#svl
//
//	pbepotential.svl    Electrostatic Potential by solving Poisson-Boltzmann Equation
//
//	21-dec-2010 (ti) keep the first 5 elements of returned value of pboltz_IonREdqC (for MOE2010.10)
//	28-jun-2010 (kt) The conditions are written in Grid Analyzer.
//      25-jun-2010 (kt) created
//
// COPYRIGHT (C) 2010 RYOKA SYSTEMS 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 title	'MOE-ElectrostaticPotential'
#set class	'MOE:simulation'
#set version	'2010.1221'
#set main 	'pbepotential'

function pboltz_AtomParameters,
pboltz_Boundary,
pboltz_ChargeDensity,
pboltz_GridShape,
pboltz_IonREdqC,
pboltz_Solve,
pboltz_ilogK,
pboltz_logD,
pboltz_vdwGrid;

function _Atoms;

local function pbepotential_panel []
    local wkey = WindowCreate [
	windowName:'Poisson-Boltzmann Electrostatic Potential'
	,title: 'Poisson-Boltzmann Electrostatic Potential', name: 'panel'
	,text: ['OK','Cancel'], onTrigger: ['return','exit']
	,Hbox:[
	    Option:[name:'atoms',title:'Atoms:',bubbleHelp:'Using atoms for calculation'],
	    Checkbox:[name:'vis_only',text:'Visible Only',bubbleHelp:'If on, hidden atoms are ignored.']
	]
	,Text:[name:'temparature', title:'Temparature:', type:'int',bubbleHelp:'Temperature of the system in Kelvin']
	,Separator:[]
	,Hbox:[
	    Text:[name:'idiele', title:'Interior Dielectric:', type:'real', bubbleHelp:'Interior dielectric of the atoms']
	    ,Text:[name:'ediele', title:'Exterior Dielectric:', type:'real', bubbleHelp:'Interior dielectric of the solvent']
	]
	,Text:[name:'offset', title:'Dielectric Offset:', type:'real', bubbleHelp:'Distance from the solute atoms of the solvent region']
	,Separator:[]
	,Hbox:[
	    Text:[name:'salt',title:'Salt:',bubbleHelp:'Chemical formula of ionic salts']
	    ,Text:[name:'salt_conc', title:'Salt Concentration:',type:'real',bubbleHelp:'Concentration of salt in \'mol/L\'']
	]
	,Separator:[]
	,Hbox:[
	    Text:[
		name:'gspace',title:'Grid Spacing:',
		bubbleHelp:
		    '   <= 5: Size of lattice\n'
		    '     > 5: Number of grid points along each axis. The value will be rounded.\n'
		    'AxBxC: Number of grid points along each axis. The each value must be integer.'
	    ]
	    ,Text:[
		name:'gextend', title:'Grid Extend:',type:'real',
		bubbleHelp:
		    'The bounding box is extended in each direction\n'
		    'by at least the value of the maximum of \'Grid Extend\',\n'
		    '3.0 and three times the maximum radius.'
	    ]
	]
	,Checkbox:[
	    name:'linearized',text:'Linearized Poisson-Boltzmann equation',
	    bubbleHelp:'If on, A linearized solution is produced'
	]
    ];
    local ctags = uniq cTag Chains[];
    ctags = ctags | app alltrue ctags;
    WindowSetAttr [wkey, [atoms:[text:cat [
	'Receptor Atoms', 'All Atoms', 'Selected Atoms', 'Unselected Atoms',
	'Selected Residues', 'Selected Chains', ctags
    ]]]];
    WindowSetAttr [wkey, [salt:[shortcut:[
	'None','NaF','NaCl','NaBr','NaI','LiF','LiCl',
	'LiBr','LiI','KF','KCl','KBr','KI','CaCl2','MgCl2'
    ]]]];
    WindowSetAttr [wkey, [gspace:[shortcut:[
	'33x33x33','65x65x65','0.25','0.5','0.75','1.0','1.5','2.0'
    ]]]];
    WindowSetData[wkey, [
	temparature:300, idiele:1, ediele:80, offset:0,
	salt:'NaCl',salt_conc:0.1,gspace:'0.75',gextend:5
    ]];
    WindowShow wkey;
    local val = first WindowWait wkey;
    WindowDestroy wkey;
    return val;
endfunction

const DEFALT_ARGS = [
    atoms:'All Atoms',
    vis_only:0,
    temparature:300,
    idiele:1,
    ediele:80,
    offset:0,
    salt:'NaCl',
    salt_conc:0.1,
    gspace:'0.75',
    gextend:5,
    linearized:0,
    return1:1
];

global function pbepotential args
    if isnull args then
	args = pbepotential_panel[];
	args.return1 = 0;
    endif

    args = cat [args, DEFALT_ARGS];
    args = args | m_uniq tags args;
    local return1 = args.return1;

    local atoms = [];
    if uniq type args.atoms === 'num' then
	if uniq oType args.atoms === 'atom' then
	    atoms = args.atoms;
	endif
    elseif args.atoms == 'Receptor Atoms' then
	atoms = _Atoms '$$receptor';
    elseif args.atoms == 'All Atoms' then
	atoms = Atoms[];
    elseif args.atoms == 'Selected Atoms' then
	atoms = SelectedAtoms[];
    elseif args.atoms == 'Unselected Atoms' then
	atoms = Atoms[];
	atoms = atoms | not aSelected atoms;
    elseif args.atoms == 'Selected Residues' then
	atoms = cat rAtoms SelectedResidues[];
    elseif args.atoms == 'Selected Chains' then
	atoms = cat cAtoms SelectedChains[];
    elseif type args.atoms == 'tok' then
	local chains = Chains[];
	chains = chains | cTag chains == args.atoms;
	atoms = cat cAtoms chains;
    endif
    if isnull atoms then
	atoms = Atoms[];
    endif
    if args.vis_only then
	atoms = atoms | not aHidden atoms;
    endif

    local pos = aPos atoms, q = aCharge atoms;
    local pos_REd = pboltz_AtomParameters atoms;

        // calculate a grid shape the encompasses the system with
	// a grid spacing of 0.75 angstroms

    local shape = pboltz_GridShape [pos, pos_REd(1), args.gextend, args.gspace];

        // calculate the log of the dielectric field

    local logD = pboltz_logD [ shape, pos, pos_REd, [
        d_in  : args.idiele,
        d_out : args.ediele,
	d_offset: args.offset
    ]];

        // calculate the ion densities

    if args.salt == 'None' then
	local salt_conc = 0;
	local ion_REdqC = [];
    else
	salt_conc = args.salt_conc;
	ion_REdqC = keep [(pboltz_IonREdqC args.salt), 5] * [1,1,1,1,salt_conc];
    endif

    local ilogK = pboltz_ilogK [shape, pos, pos_REd, ion_REdqC, logD];
    local iQ = ion_REdqC(4);

        // calculate the charge density using R/2.5 as the standard
	// deviations and then calculate the density with the folded
        // constants and 1/dielectric for the solver

    local f = pboltz_ChargeDensity [shape, pos, q, inv 2.5 * pos_REd(1)];
    local F = f * exp (- logD);

	// calculate the initial estimate including the boundary

    const AVOGADRO = 6.022045e23;
    local T = args.temparature;

    local inv_kappa = sqrt (
	  (AVOGADRO * 1e-27 * add (sqr ion_REdqC(4) * ion_REdqC(5)))
	* (4 * PI * COULOMB_SCALE) / (KBOLTZ * T)
    );

    local H = app maxE (shape - app shiftr shape);
    local H1 = cbrt mul (H / 4);
    local sigma = norm [pos_REd(1) * 2.5, H1];

    local u = pboltz_Boundary [shape, 0, pos, q, sigma, inv_kappa, logD];

    local msg_tok = tok_cat [
	'Solving ',
	select['linear','nonlinear',args.linearized],
	' Poisson-Boltzmann equation...'
    ];
    local msg = Message [0, msg_tok];

        // solve the Poisson-Boltzmann Equation
	// and then calculate the total energy and scale to kcal/mol

    u = pboltz_Solve [shape, logD, iQ, ilogK, u, F, [T:T, linear:args.linearized, verbose:1]];
    u = u * 4 * PI * COULOMB_SCALE;
    Message[msg,''];

    if return1 then
	return [u,shape];
    else
	local gridname = tok_cat[
	    'Electrostatics ',
	    token swrite [
		'Salt:{t:} {n:3.1}mol/L, Interior:{n:}, Exterior:{n:}, T:{n:}',
		args.salt, args.salt_conc, args.idiele, args.ediele, args.temparature
	    ]
	];
	local gshow_ini = [cm_level1:-5,cm_level2:0,cm_level3:5];
	run ['gridshow.svl', [[u, shape], gridname, gshow_ini]];
	return [];
    endif
endfunction

#eof