#svl
//
//	unshared_vol.svl
//
//	24-nov-2020 (kk) fixed the year of COPYRIGHT
//	17-jan-2020 (kk) added db_AtomRegion function and related functions
//	30-sep-2011 (na) added "None" to Atoms2 menu
//	27-oct-2010 (ti) output volume information into CLI
//	30-sep-2010 (ti) bug fixed (Unselected Residues -> Unselected Chains)
//	30-sep-2010 (ti) support dummy atoms
//	08-apr-2010 (ti) added VDW allowance option
//	17-apr-2009 (ti) added TR2 menu
//	16-apr-2009 (ti) added "Near Atoms1" to Atoms2 menu
//	04-feb-2009 (ti) modified default menu for amoe
//	22-jan-2009 (kt) GUI
//	16-jan-2009 (kt) bug fix (translation of positions)
//	15-jan-2009 (kt) added 'unshared_region'
//      14-jan-2009 (kt) create
//
// COPYRIGHT (C) 2009-2020 MOLSIS INC.  ALL RIGHTS RESERVED.
//
// PERMISSION TO USE, COPY, MODIFY AND DISTRIBUTE THIS SOFTWARE IS HEREBY
// GRANTED PROVIDED THAT: (1) UNMODIFIED OR FUNCTIONALLY EQUIVALENT CODE
// DERIVED FROM THIS SOFTWARE MUST CONTAIN THIS NOTICE; (2) ALL CODE DERIVED
// FROM THIS SOFTWARE MUST ACKNOWLEDGE THE AUTHOR(S) AND INSTITUTION(S); (3)
// THE NAMES OF THE AUTHOR(S) AND INSTITUTION(S) NOT BE USED IN ADVERTISING
// OR PUBLICITY PERTAINING TO THE DISTRIBUTION OF THE SOFTWARE WITHOUT
// SPECIFIC, WRITTEN PRIOR PERMISSION; (4) ALL CODE DERIVED FROM THIS SOFTWARE
// BE EXECUTED WITH THE MOLECULAR OPERATING ENVIRONMENT (MOE) LICENSED FROM
// RYOKA SYSTEMS INC.
//
// RYOKA SYSTEMS INC. DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS
// SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS,
// AND IN NO EVENT SHALL RYOKA SYSTEMS INC. BE LIABLE FOR ANY
// SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER
// RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF
// CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
// CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

#set main 	'AtomRegion'
#set title 	'AtomRegion'
#set class 	'MOE:interactive'
#set version    'MOLSIS 2020.11.24'

function _Atoms;

const SOLID = 1;
const DOT = 1;
const G_SPACE = 0.5; // Grid space default: 0.5

static gkey;
static gkeydot;

local function grid_SphereCalc [[p1, rad1], [p2, rad2], operation, gspace, trade_flag, db]

    local pos = app cat tr [p1,p2];
    local radius = cat [rad1, rad2];
    local center = app add pos * invz l_length pos;
    pos = pos - center;					// translate
    local [e,Q] = rot3d_CovarianceEigenSystem pos;
    pos = app add (Q * [pos]);				// rotate

    [p1, p2] = tr apt split [pos, [app length [p1(1), p2(1)]]];

    if length radius === 0 then return 0; endif

	// calculate the grid shape from the origin, half dimensions
	// of the box and the grid spacing

    local max_radius = maxE radius;
    local maxpos = app maxE pos, minpos = app minE pos;
    local org = 0.5 * (maxpos + minpos);

    local h   = gspace;					// grid spacing
    local dim = ceil ( (0.5 * (maxpos - minpos) + max_radius + 0.001) / h );
    local shape = org + (app igen inc (2*dim) - inc dim) * h;

	// get the grid coordinates and use as a probe for the molecule
	// intersection test

    local gpos = grid_coord [shape, igen mul app length shape];

    local pkey = prox_open [max_radius+1, p1, rad1];
    local px1 = first prox_find [pkey, gpos, 0];
    prox_close pkey;

    pkey = prox_open [max_radius+1, p2, rad2];
    local px2 = first prox_find [pkey, gpos, 0];
    prox_close pkey;

    if istrue db then return [px1, px2]; endif

    if operation == 'AND' then
	local px = px1 and px2;
    elseif operation == 'DIFF' then
	px = px1 and not px2;
    else
	px = px1 or px2;
    endif

    gpos = app add (matinv Q * [gpos]);
    gpos = gpos + center;

    if istrue trade_flag then
	[px2, px1] = [px1, px2]; 
    endif
    write ['----------------------\n'];
    write ['Atoms1 = {n:} A^3\n', add (app istrue px1) * cube gspace];
    write ['Atoms2 = {n:} A^3\n', add (app istrue px2) * cube gspace];
    write ['Total = {n:} A^3\n', add px * cube gspace];
    write ['----------------------\n'];

    return [px, gpos];
endfunction


local function  grid_SphereReg [[p1, rad1], [p2, rad2], operation, gspace, opt]
    local pos = app cat tr [p1,p2];
    local radius = cat [rad1, rad2];
    local center = app add pos * invz l_length pos;
    pos = pos - center;					// translate
    local [e,Q] = rot3d_CovarianceEigenSystem pos;
    pos = app add (Q * [pos]);				// rotate
    [p1, p2] = tr apt split [pos, [app length [p1(1), p2(1)]]];

    if length radius === 0 then return 0; endif

	// calculate the grid shape from the origin, half dimensions
	// of the box and the grid spacing

    local max_radius = maxE radius;
    local maxpos = app maxE pos, minpos = app minE pos;
    local org = 0.5 * (maxpos + minpos);

    local h   = gspace;					// grid spacing
    local dim = ceil ( (0.5 * (maxpos - minpos) + max_radius + 0.001) / h );
    local shape = org + (app igen inc (2*dim) - inc dim) * h;

	// get the grid coordinates and use as a probe for the molecule
	// intersection test

    local gpos = grid_coord [shape, igen mul app length shape];

    local pkey = prox_open [max_radius+1, p1, rad1];
    local px1 = first prox_find [pkey, gpos, 0];
    prox_close pkey;

    pkey = prox_open [max_radius+1, p2, rad2];
    local px2 = first prox_find [pkey, gpos, 0];
    prox_close pkey;

    if operation == 'AND' then
	local px = px1 and px2;
    elseif operation == 'DIFF' then
	px = px1 and not px2;
    else
	px = px1 or px2;
    endif
    local color = opt.color_sol;
    local dims = [1,2,3];
    local gtype = 3;
    local level = 0.5;

    local [isovtx, isoidx] = grid_iso [px, shape, level, dims];
    local points = grid_isovtx_coord [px, shape, isovtx, level, dims];
    local normals = - grid_isovtx_grad [px, shape, isovtx, level, dims];
    local idx = grid_isoidx_emit [isoidx, gtype, -1];

    points = app add (matinv Q * [points]);
    points = points + center;
    normals = app add (matinv Q * [normals]);

    if opt.sol then
	gkey = GCreateUnique opt.gname;
    	GSetTransparency [ gkey, [opt.tr1]  ];
	GVertex  [gkey, gtype, idx, color, points, normals];
    endif

    if opt.dot then
	gkeydot = GCreateUnique tok_cat [opt.gname, '_dot'];
	GSetTransparency [ gkeydot, [opt.tr2]  ];
	gpos = app add (matinv Q * [gpos]);
	gpos = gpos + center;
	gpos = gpos || [px];
	GVertex  [gkeydot , 1, igen l_length gpos, opt.color_dot, gpos];
    endif
endfunction

local function diff_SphereVolume [[p1, rad1],[p2, rad2]]

    local pos = app cat tr [p1,p2];
    local radius = cat [rad1, rad2];

    if length radius === 0 then return 0; endif

	// calculate the grid shape from the origin, half dimensions
	// of the box and the grid spacing

    local max_radius = maxE radius;
    local maxpos = app maxE pos, minpos = app minE pos;
    local org = 0.5 * (maxpos + minpos);

    local h   = G_SPACE;					// grid spacing
    local dim = ceil ( (0.5 * (maxpos - minpos) + max_radius + 0.001) / h );
    local shape = org + (app igen inc (2*dim) - inc dim) * h;

	// get the grid coordinates and use as a probe for the molecule
	// intersection test

    local [x,y,z] = grid_coord [shape, igen mul app length shape];

    local pkey = prox_open [max_radius+1, p1, rad1];
    local px1 = first prox_find [pkey, [x,y,z], 0];
    prox_close pkey;

    pkey = prox_open [max_radius+1, p2, rad2];
    local px2 = first prox_find [pkey, [x,y,z], 0];
    prox_close pkey;

    local npoints = add andE[px1, xorE[px1,px2]];
    return npoints * cube h;
endfunction


global function unshared_vol [atoms1, atoms2]

    if andE app isnull [atoms1, atoms2] then
	atoms1 = SelectedAtoms[];
	atoms2 = diff [Atoms[], atoms1];
    endif

    local el1 = aElement atoms1;
    local el2 = aElement atoms2;
    local p1 = aPos atoms1;
    local p2 = aPos atoms2;

    local pos = app cat tr [p1, p2];
    pos = pos - app add pos * invz l_length pos;	// translate
    local [e,Q] = rot3d_CovarianceEigenSystem pos;
    pos = app add (Q * [pos]);				// rotate

    [p1, p2] = tr apt split [pos, [app length [p1(1),p2(1)]]];

    local volume = diff_SphereVolume [
	[p1, el_VDW_Radius el1],[p2, el_VDW_Radius el2]
    ];

    return volume;
endfunction

const OPT_DEFALT = [color: 0xcccccc, gname: 'Unshared region' ];

global function unshared_region [atoms1, atoms2, opt]
    opt = cat [opt, OPT_DEFALT];
    opt = opt | m_uniq tags opt;

    if andE app isnull [atoms1, atoms2] then
	atoms1 = SelectedAtoms[];
	atoms2 = diff [Atoms[], atoms1];
    endif

    local el1 = aElement atoms1;
    local el2 = aElement atoms2;
    local rad1 = el_VDW_Radius el1;
    local rad2 = el_VDW_Radius el2;
    local p1 = aPos atoms1;
    local p2 = aPos atoms2;

    local pos = app cat tr [p1, p2];
    local trpos = app add pos * invz l_length pos;
    pos = pos - trpos;					// translate
    local [e,Q] = rot3d_CovarianceEigenSystem pos;
    pos = app add (Q * [pos]);				// rotate

    [p1, p2] = tr apt split [pos, [app length [p1(1),p2(1)]]];

    local radius = cat app el_VDW_Radius [el1, el2];

	// calculate the grid shape from the origin, half dimensions
	// of the box and the grid spacing

    local max_radius = maxE radius;
    local maxpos = app maxE pos, minpos = app minE pos;
    local org = 0.5 * (maxpos + minpos);

    local h   = G_SPACE;					// grid spacing
    local dim = ceil ( (0.5 * (maxpos - minpos) + max_radius + 0.001) / h );
    local shape = org + (app igen inc (2*dim) - inc dim) * h;

	// get the grid coordinates and use as a probe for the molecule
	// intersection test

    local gpos = grid_coord [shape, igen mul app length shape];

    local pkey = prox_open [max_radius+1, p1, rad1];
    local px1 = first prox_find [pkey, gpos, 0];
    prox_close pkey;

    pkey = prox_open [max_radius+1, p2, rad2];
    local px2 = first prox_find [pkey, gpos, 0];
    prox_close pkey;

    local px = andE[px1, xorE[px1,px2]];
    local color = opt.color;
    local dims = [1,2,3];
    local gtype = 3;
    local level = 0.5;

    local [isovtx, isoidx] = grid_iso [px, shape, level, dims];
    local points = grid_isovtx_coord [px, shape, isovtx, level, dims];
    local normals = - grid_isovtx_grad [px, shape, isovtx, level, dims];
    local idx = grid_isoidx_emit [isoidx, gtype, -1];

    points = app add (matinv Q * [points]);
    points = points + trpos;
    normals = app add (matinv Q * [normals]);

#if SOLID
    local gname = tok_cat [opt.gname, ' (sol)'];
    GVertex  [GCreateUnique gname, gtype, idx, color, points, normals];
#endif

#if DOT
    gname = tok_cat [opt.gname, ' (dot)'];
    gpos = app add (matinv Q * [gpos]);
    gpos = gpos + trpos;
    gpos = gpos || [px];
    GVertex  [GCreateUnique 'Unshared region (dot)' , 1, igen l_length gpos, 0x00FF00, gpos];
#endif
endfunction

const ATOMS_OPTION_DEF = [
	'Selected atoms','Unselected atoms', 'Selected residues', 'Unselected residues', 'Selected chains', 'Unselected chains',
	'Dummies(B factor)', 'Dummies(HYD:C, LPA:O)'
];


global function AtomRegion[]

    if WindowShow ['AtomRegion', 1] then return; endif

    local atoms_option_def = cat ['Ligand Atoms', 'Receptor and Solvent', ATOMS_OPTION_DEF];

    local wkey = WindowCreate [
	windowName: 'AtomRegion',
	name: 'shell',
	title: 'AtomRegion',
	text:['Apply', 'Cancel'],
	onTrigger: ['return', 'exit'],

	Text: [name:'gname', title: 'Graphic Name:', len:30],
	Hbox:[
	    title: 'Graphics:',
	    Hbox:[
		shadow:'etched-out',
		Mbox: [
		    columns: 4,
		    Checkbox:[name:'sol', text: 'Solid', font:'mediumBold', onTrigger:'return'],
		    Separator:[vertical:1, shadow:'none'],
		    Color:[name:'color_sol', minWidth:3, title:'Color:', onTrigger:'return'],
		    Slider :[name:'tr1', title:'TR:', range: [0,255,1],width:5, len:5, onTrigger:'return'],

		    Checkbox:[name:'dot', text: 'Dot', font:'mediumBold', onTrigger:'return'],
		    Separator:[vertical:1, shadow:'none'],
		    Color:[name:'color_dot', minWidth:3, title:'Color:', onTrigger:'return'],
		    Slider :[name:'tr2', title:'TR:', range: [0,255,1],width:5, len:5, onTrigger:'return']
		]
	        //Slider :[name:'tr1', title:'TR:', range: [0,255,1],width:5, len:5]
	    ]
	],
	Hbox:[
	    Text: [name:'gspace', title: 'Grid space', onTrigger:'return', type:'real', len:6, shortcut:['0.5','0.4','0.3','0.2','0.1']],
	    Slider :[name:'vdw_allow', title:' VDW radius:', range: [0, 200, 1], width:5, onTrigger:'return'],
	    Text:[name:'vdw_val', len:4, type:'real', onTrigger:'return'],
	    Label:[text:'%']	    
	],
	Separator: [],
	Hbox: [
	    Option:[name:'atoms1', title:'Atoms1:', text:atoms_option_def] //onTrigger:'return']
	],
	Hbox: [    
	    Option:[name:'atoms2', title:'Atoms2:', text:cat ['Near Atoms1', atoms_option_def, 'None'], onTrigger:'return'],
	    Text: [name:'area', title: ' Within:', onTrigger:'return', type:'real', len:6, shortcut:['3.0','3.5','4.0','4.5','5.0','5.5','6.0']]
	],
	Hbox: [
	    Option:[
		name: 'oper',
		title: 'Operation:', 
//		text: ['Atoms1 - Atoms2', 'Atoms2 - Atoms1', 'Atoms1 AND Atoms2', 'Atoms1 OR Atoms2'], //original
		text: ['Atoms1 AND Atoms2', 'Atoms1 OR Atoms2', 'Atoms1 - Atoms2', 'Atoms2 - Atoms1'], //modify default
		onTrigger: 'return'
	    ]
	],
	Hbox:[
	    Text: [name:'volume', title: 'Volume:', len:10, type:'real'],
	    Label:[text:'Å ^ 3', font:'largeBold']
	]
    ];

    WindowSetAttr [wkey, [volume:[sensitive:0]]];
    local atoms_opt = app tok_cat tr ['Chain #', totok x_id Chains[], '- ', cName Chains []];
    WindowSetAttr [wkey, tag[['atoms1', 'atoms2'],[[text:cat[atoms_option_def, atoms_opt]]]]]; 
    WindowSetAttr [wkey, [atoms1:[text:cat[atoms_option_def, atoms_opt]]]];
    WindowSetAttr [wkey, [atoms2:[text:cat['Near Atoms1', atoms_option_def, atoms_opt, 'None']]]];

    WindowSetData [wkey, [vdw_allow:100, vdw_val:100]];
    WindowSetData [wkey, [gname:'Atom Region', sol:1, color_sol:0xCCCCCC, dot:1, color_dot:0x00FF00, 
			gspace:0.5, atoms2:'Unselected atoms']];
    WindowShow [wkey, 1];

    local function get_atoms opt
	local atoms = Atoms[];
	local res = Residues[];
	local chains = Chains[];
	if opt == 'Selected atoms' then
	    return atoms | aSelected atoms;
	elseif opt == 'Unselected atoms' then
	    return atoms | not aSelected atoms;
	elseif opt == 'Selected residues' then
	    return cat rAtoms (res | rSelected res);
	elseif opt == 'Unselected residues' then
	    return cat rAtoms (res | not rSelected res);
	elseif opt == 'Selected chains' then
	    return cat cAtoms (chains | cSelected chains);
	elseif opt == 'Unselected chains' then
	    return cat cAtoms (chains | not cSelected chains);
	elseif opt == 'Ligand Atoms' then
	    return _Atoms '$$ligand';
	elseif opt == 'Receptor and Solvent' then
	    return cat [_Atoms '$$receptor', _Atoms '$$solvent'];
	elseif opt == 'Dummies(B factor)' or opt == 'Dummies(HYD:C, LPA:O)' then
	    return atoms | m_join [aElement atoms, 'LP'];
	else
	    local cnum = atoi tok_drop [opt, 7];
	    return cat cAtoms cat chains[indexof [cnum, cNumber chains]];
	endif
    endfunction

    local function _grid_SphereCalc val
	local atoms1 = get_atoms val.atoms1;
	local atoms2;
	if val.atoms2 === 'Near Atoms1' then
	    local r = val.area;
	    local pkey = prox_open [r, aPos Atoms [], r/2];
	    local idx = second prox_find [pkey, aPos atoms1, r/2];
	    prox_close pkey;
	    atoms2 = get [Atoms [], idx];
	    atoms2 = diff [atoms2, atoms1];
	elseif val.atoms2 === 'None' then	// 30-sep-2011 (na) added Atoms2:'None'
		atoms2 = atoms1;
	    val.oper = 'Atoms1 AND Atoms2';
	else
	    atoms2 = get_atoms val.atoms2;
	endif

	local [target1, target2] = [val.atoms1, val.atoms2];
	if val.atoms2 === 'None' then target2 = target1; endif	// 30-sep-2011 (na) added Atoms2:'None'

	local trade = 0; 

	if orE app isnull [atoms1, atoms2] then
	    return[];
	endif
	if val.oper == 'Atoms1 - Atoms2' then
	    local oper = 'DIFF';
	elseif val.oper == 'Atoms2 - Atoms1' then
	    [atoms1, atoms2] = [atoms2, atoms1];
	    [target1, target2] = [target2, target1];
	    oper = 'DIFF';
	    trade = 1;	   // make a trade of atoms1 for atoms2
	elseif val.oper == 'Atoms1 AND Atoms2' then
	    oper = 'AND';
	else
	    oper = 'OR';
	endif

	local p1 = aPos atoms1;
	local rad1 = (el_VDW_Radius aElement atoms1);
	local aname1 = aName atoms1;

	if target1 == 'Dummies(HYD:C, LPA:O)' then
	    rad1 | aname1 == 'LPA' = el_VDW_Radius 'O';
	    rad1 | aname1 == 'HYD' = el_VDW_Radius 'C';
	elseif target1 == 'Dummies(B factor)' then
	    rad1 = aTempFactor atoms1;
	    rad1 | rad1 < 0 = 0;
	endif
	rad1 = rad1 * (val.vdw_val/100);

	local p2 = aPos atoms2;
	local rad2 = (el_VDW_Radius aElement atoms2);
	local aname2 = aName atoms2;

	if target2 == 'Dummies(HYD:C, LPA:O)' then
	    rad2 | aname2 == 'LPA' = el_VDW_Radius 'O';
	    rad2 | aname2 == 'HYD' = el_VDW_Radius 'C';
	elseif target2 == 'Dummies(B factor)' then

	    rad2 = aTempFactor atoms2;
	    rad2 | rad2 < 0 = 0;
	endif
	rad2 = rad2 * (val.vdw_val/100);

	local gspace = val.gspace;
	return grid_SphereCalc [[p1, rad1], [p2, rad2], oper, gspace, trade];
    endfunction

    local function _grid_SphereReg val
	local atoms1 = get_atoms val.atoms1;
//	local atoms2 = get_atoms val.atoms2;
	local atoms2;
	if val.atoms2 == 'Near Atoms1' then
	    local r = val.area;
	    local pkey = prox_open [r, aPos Atoms [], r/2];
	    local idx = second prox_find [pkey, aPos atoms1, r/2];
	    prox_close pkey;
	    atoms2 = get [Atoms [], idx];
	    atoms2 = diff [atoms2, atoms1];
	elseif val.atoms2 === 'None' then	// 30-sep-2011 (na) added Atoms2:'None'
		atoms2 = atoms1;
	    val.oper = 'Atoms1 AND Atoms2';
	else
	    atoms2 = get_atoms val.atoms2;
	endif

	if orE app isnull [atoms1, atoms2] then
	    return[];
	endif

	local [target1, target2] = [val.atoms1, val.atoms2];
	if val.atoms2 === 'None' then target2 = target1; endif	// 30-sep-2011 (na) added Atoms2:'None'

	if val.oper === 'Atoms1 - Atoms2' then
	    local oper = 'DIFF';
	elseif val.oper === 'Atoms2 - Atoms1' then
	    [atoms1, atoms2] = [atoms2, atoms1];
	    [target1, target2] = [target2, target1];
	    oper = 'DIFF';
	elseif val.oper === 'Atoms1 AND Atoms2' then
	    oper = 'AND';
	else
	    oper = 'OR';
	endif

	local p1 = aPos atoms1;
	local rad1 = (el_VDW_Radius aElement atoms1);
	local aname1 = aName atoms1;
	if target1 === 'Dummies(HYD:C, LPA:O)' then
	    rad1 | aname1 == 'LPA' = el_VDW_Radius 'O';
	    rad1 | aname1 == 'HYD' = el_VDW_Radius 'C';
	elseif target1 === 'Dummies(B factor)' then
	    rad1 = aTempFactor atoms1;
	    rad1 | rad1 < 0 = 0;
	endif
	rad1 = rad1 * (val.vdw_val/100);

	local p2 = aPos atoms2;
	local rad2 = (el_VDW_Radius aElement atoms2);
	local aname2 = aName atoms2;
	if target2 === 'Dummies(HYD:C, LPA:O)' then
	    rad2 | aname2 == 'LPA' = el_VDW_Radius 'O';
	    rad2 | aname2 == 'HYD' = el_VDW_Radius 'C';
	elseif target2 === 'Dummies(B factor)' then
	    rad2 = aTempFactor atoms2;
	    rad2 | rad2 < 0 = 0;
	endif
	rad2 = rad2 * (val.vdw_val/100);

	local gspace = val.gspace;
	return grid_SphereReg [[p1, rad1], [p2, rad2], oper, gspace, val];
    endfunction

    WindowSetData [wkey, [
	atoms2:['Near Atoms1'],
	area:[4.5]
    ]];

    loop
	local [val,trg] = WindowWait wkey;
	
	if trg == 'shell' and val.shell == 'Apply' then
	    UndoSnapshot[];
	    local [px, pos] = _grid_SphereCalc val;
	    if isnull px then
		WindowSetData [wkey, [volume:[]]];
		continue;
	    endif
	    WindowSetData [wkey, [volume:add px * cube val.gspace]];
	    if andE not [val.sol, val.dot] then
		continue;
	    else
	    	_grid_SphereReg val;
	    endif
	elseif trg == 'atoms2' then
	    WindowSetAttr [wkey, [area:[sensitive:[val.atoms2 == 'Near Atoms1']]]];
	    WindowSetAttr [wkey, [oper:[sensitive:[val.atoms2 <> 'None']]]];
	elseif trg == 'sol' or trg == 'dot' then
	    WindowSetAttr [wkey, [gname:[sensitive:orE [val.sol, val.dot]]]];
	elseif trg == 'tr1' then
	    if istrue gkey and istrue m_join [gkey, GKeyList []] then
		GSetTransparency [ gkey, [val.tr1]  ];
	    else
		continue;
	    endif
	elseif trg == 'tr2' then
	    if istrue gkeydot and istrue m_join [gkeydot, GKeyList []] then
	    	GSetTransparency [ gkeydot, [val.tr2]  ];		
	    else
		continue;
	    endif
	elseif trg == 'vdw_allow' then
	    WindowSetData [wkey, [vdw_val:val.vdw_allow / 100]];
	    WindowSetData [wkey, [vdw_val:val.vdw_allow]];
	elseif trg == 'vdw_val' then
	    //WindowSetData [wkey, [vdw_allow:val.vdw_val * 100]];
	    WindowSetData [wkey, [vdw_allow:val.vdw_val]];
 	elseif trg == 'color_sol' then
	    if istrue gkey and istrue m_join [gkey, GKeyList []] then
	    	local gdata = GVertexData gkey;
	    	gdata(1)(3) = val.color_sol;
	    	GClearVertexData gkey;
	    	GVertex prepend [cat gdata, gkey];  
	    else
	    	continue;
	    endif
	elseif trg == 'color_dot' then
	    if istrue gkeydot and istrue m_join [gkeydot, GKeyList []] then
	    	gdata = GVertexData gkeydot;
	    	gdata(1)(3) = val.color_dot;
	    	GClearVertexData gkeydot;
	    	GVertex prepend [cat gdata, gkeydot];  
	    else
	    	continue;
	    endif
#if 0
	elseif m_join [trg, ['atoms1','atoms2','oper']] then
	    [px, pos] = _grid_SphereCalc val;
	    if isnull px then
		WindowSetData [wkey, [volume:[]]];
		continue;
	    endif
	    WindowSetData [wkey, [volume:add px * cube val.gspace]];
#endif
	endif
    endloop

endfunction


local function mdb_fileset [wkey, mdb_file]
    local fn = FilePrompt [
	name: mdb_file,
	title: 'Select MDB File',
	filter: '*.mdb',
	mode: 'open'
    ];
    if istrue wkey and istrue fn then
	WindowSetData [wkey, [mdb_file: fn]];
    endif

    return fn;
endfunction

local function fieldset [wkey, val]
    if istrue ('file' === ftype val.mdb_file) then
	local dbkey = db_Open [val.mdb_file, 'read'];
	local [flds, typs] = db_Fields dbkey;
	local mfld_lst = flds | m_join [typs, 'molecule'];

	if isfalse mfld_lst then
	    Warning 'There is no molecular type field.';
	    mfld_lst = '';
	    val.mfield = '';
	else
	    val.mfield = first mfld_lst;
	endif
	
	db_Close dbkey;
    else
	mfld_lst = '';
	val.mfield = '';
    endif

    if istrue wkey and istrue mfld_lst then
	WindowSetAttr [wkey, [mfield: [text: mfld_lst]]];
	WindowSetData [wkey, val];
    endif
    return val;
endfunction

local function fileview [wkey, val]
    if isfalse val.mdb_file then
	val.mdb_file = mdb_fileset [wkey, val.mdb_file];
	fieldset [wkey, val];
    endif
    if istrue ('file' === ftype val.mdb_file) then
	local dbkey = dbv_Open val.mdb_file;
	db_Close dbkey;
    else
	Warning 'This is not a valid file.';
    endif
endfunction

local function get_atoms opt
    local atoms = Atoms[];
    local res = Residues[];
    local chains = Chains[];
    if opt == 'Selected atoms' then
	return atoms | aSelected atoms;
    elseif opt == 'Unselected atoms' then
	return atoms | not aSelected atoms;
    elseif opt == 'Selected residues' then
	return cat rAtoms (res | rSelected res);
    elseif opt == 'Unselected residues' then
	return cat rAtoms (res | not rSelected res);
    elseif opt == 'Selected chains' then
	return cat cAtoms (chains | cSelected chains);
    elseif opt == 'Unselected chains' then
	return cat cAtoms (chains | not cSelected chains);
    elseif opt == 'Ligand Atoms' then
	return _Atoms '$$ligand';
    elseif opt == 'Receptor and Solvent' then
	return cat [_Atoms '$$receptor', _Atoms '$$solvent'];
    elseif opt == 'Dummies(B factor)' or opt == 'Dummies(HYD:C, LPA:O)' then
	return atoms | m_join [aElement atoms, 'LP'];
    else
	local cnum = atoi tok_drop [opt, 7];
	return cat cAtoms cat chains[indexof [cnum, cNumber chains]];
    endif
endfunction

local function get_SphereRadius [atoms, vdw_val, expr]
    local rad = (el_VDW_Radius aElement atoms);
    local aname = aName atoms;

    if istrue expr then
	if expr == 'Dummies(HYD:C, LPA:O)' then
	    rad | aname == 'LPA' = el_VDW_Radius 'O';
	    rad | aname == 'HYD' = el_VDW_Radius 'C';
	elseif expr == 'Dummies(B factor)' then
	    rad = aTempFactor atoms;
	    rad | rad < 0 = 0;
	endif
    endif

    return rad * (vdw_val/100);
endfunction

const OUTPUT_DEFAULTS = [
    tanimoto: 1,
    intersect: 0,
    union: 0,
    diff12: 0,
    diff21: 0
];

const DB_ATOMREGION_DEFAULTS = [
    mdb_file: '',
    esel: 0,
    mfield: '',
    atoms2: '',
    gspace: G_SPACE,
    vdw_val: 100,
    output: OUTPUT_DEFAULTS
];

global function _db_AtomRegion opt
    opt = tagcat [opt, DB_ATOMREGION_DEFAULTS];
    local output = tagcat [opt.output, OUTPUT_DEFAULTS];

    if isfalse opt.mdb_file then
	opt.mdb_file = mdb_fileset [];
	if isfalse opt.mdb_file then exit []; endif
    endif
    local mdb = db_Open [opt.mdb_file, 'read-write'];

    local [flds, typs] = db_Fields mdb;
    local mflds = flds | (typs == 'molecule');
    if isfalse opt.mfield then
	opt.mfield = first mflds;
    elseif isfalse (opt.mfield == flds) then
	exit twrite ['{} does not exist in {}', opt.mfield, ftail opt.mdb_file];
    elseif isfalse (opt.mfield == mflds) then
	exit twrite ['{} is not a molecule field', opt.mfield];
    endif

    if opt.esel then
	if isfalse dbv_nSelectedEntries mdb then
	    exit 'There are no selected entries';
	endif
	local ents = dbv_SelectedEntries mdb;
    else
	ents = db_Entries mdb;
    endif

    if isfalse opt.atoms2 then
	local atoms2 = SelectedAtoms [];
	if isfalse atoms2 then atoms2 = Atoms []; endif
	if isfalse atoms2 then
	    exit 'There are no atoms in MOE window';
	endif
    else
	atoms2 = get_atoms opt.atoms2;
	if isfalse atoms2 then
	    exit twrite ['{} is invalid expression', opt.atoms2];
	endif
    endif
    local apos2 = aPos atoms2;
    local rad2 = get_SphereRadius [atoms2, opt.vdw_val, opt.atoms2];

    output = tr output;
    local outflds = apt tok_cat ['v_', first output | second output];
    apt db_EnsureField [mdb, outflds, 'float'];

    local pdata = SystemPush atoms2;
    local db_flag = 1;
    local i;
    for i = 1, length ents loop
	local mol = cat db_ReadFields [mdb, ents(i), opt.mfield];
	local chains1 = mol_Create mol;
	local atoms1 = cat cAtoms chains1;
	local apos1 = aPos atoms1;
	local rad1 = get_SphereRadius [atoms1, opt.vdw_val];

	local [px1, px2] = grid_SphereCalc [
	    [apos1, rad1], [apos2, rad2], [], opt.gspace, 0, db_flag
	];

	local data = [];
	data.v_intersect = add (px1 and px2) * cube opt.gspace;
	data.v_union = add (px1 or px2) * cube opt.gspace;
	data.v_diff12 = add (px1 and not px2) * cube opt.gspace;
	data.v_diff21 = add (px2 and not px1) * cube opt.gspace;
	data.v_tanimoto = data.v_intersect / data.v_union;

	data = tr data;
	data = data || [m_join [first data, outflds]];
	data = tr data;
	db_Write [mdb, ents(i), data];

	oDestroy chains1;
    endloop

    SystemPop pdata;
    db_Close mdb;
    write 'done.\n';
endfunction

global function db_AtomRegion[]

    if WindowShow ['db_AtomRegion', 1] then return; endif

    local atoms_option_def = cat ['Ligand Atoms', 'Receptor and Solvent', ATOMS_OPTION_DEF];

    local wkey = WindowCreate [
	windowName: 'db_AtomRegion',
	name: 'shell',
	title: 'db_AtomRegion',
	text: ['Calculate', 'Cancel'],
	onTrigger: ['return', 'exit'],

	Hbox: [
	    FSBText: [
		name: 'mdb_file', title: 'Database:', len: 30, 
		mode:'open', onTrigger: 'return'
	    ],
	    Button: [name: 'mdb_brws'],
	    Button: [name: 'mdb_vw']
	],
	Checkbox: [name: 'esel', text: 'Selected Entries Only'],
	Hbox: [
	    Option: [name: 'mfield', title:'Molecule Field:'],
	    Label: [text: '   Used as Atoms1']
	],
	Hbox: [
	    Option:[name: 'atoms2', title: 'Atoms2:', text: atoms_option_def]
	],
	Hbox: [
	    Text: [
		name: 'gspace', title: 'Grid space:', onTrigger: 'return', 
		type: 'real', len: 6, shortcut: ['0.5','0.4','0.3','0.2','0.1']
	    ],
	    Slider: [
		name: 'vdw_allow', title: '  VDW radius:', 
		range: [0, 200, 1], width: 5, onTrigger: 'return'
	    ],
	    Text: [name: 'vdw_val', len: 4, type: 'real', onTrigger: 'return'],
	    Label: [text: '%']	    
	],
	Checkbox: [title: 'Output:', name: 'tanimoto', text: 'tanimoto (v_tanimoto)'],
	Checkbox: [name: 'intersect', text: 'Atoms1 AND Atoms2 (v_intersect)'],
	Checkbox: [name: 'union', text: 'Atoms1 OR Atoms2 (v_union)'],
	Checkbox: [name: 'diff12', text: 'Atoms1 - Atoms2 (v_diff12)'],
	Checkbox: [name: 'diff21', text: 'Atoms2 - Atoms1 (v_diff21)']
    ];

    local atoms_opt = app tok_cat tr ['Chain #', totok x_id Chains[], '- ', cName Chains []];
    local browse_icon = gr_icon [ 'open.ico', 'medium', 'foreground'];
    local view_icon = gr_icon ['fileview.ico', 'medium', 'forground'];
    WindowSetAttr [wkey, [
	atoms2: [text: cat [atoms_option_def, atoms_opt]],
	mdb_brws: [graphics: browse_icon],
	mdb_vw: [graphics: view_icon]
    ]];
    WindowSetData [wkey, [vdw_allow: 100, vdw_val: 100, gspace: G_SPACE]];
    WindowSetData [wkey, OUTPUT_DEFAULTS];

    WindowShow [wkey, 1];

    loop
	local [val,trg] = WindowWait wkey;
	
	if trg == 'shell' and val.shell == 'Calculate' then
	    local opt = val;
	    opt.output.tanimoto = val.tanimoto;
	    opt.output.intersect = val.intersect;
	    opt.output.union = val.union;
	    opt.output.diff12 = val.diff12;
	    opt.output.diff21 = val.diff21;
	    if isfalse opt.mdb_file then
		Warning 'Specify Database!';
		continue;
	    endif
	    _db_AtomRegion opt;
//	    continue;
	    exit [];
	elseif trg == 'mdb_brws' then
	    val.mdb_file = mdb_fileset [wkey, val.mdb_file];
	    val = fieldset [wkey, val];
	elseif trg == 'mdb_vw' then
	    fileview [wkey, val];
	elseif trg == 'vdw_allow' then
	    WindowSetData [wkey, [vdw_val:val.vdw_allow / 100]];
	    WindowSetData [wkey, [vdw_val:val.vdw_allow]];
	elseif trg == 'vdw_val' then
	    WindowSetData [wkey, [vdw_allow:val.vdw_val]];
	endif
    endloop

endfunction
#eof