============================
Brief description
===========================

"bornRadius" program computes the total solvation free energy and its
components, using either the canonical GB of Still et al. 
or CHA-GB (Mukhopadhyay et al.
"Introducing Charge Hydration Asymmetry into the Generalized Born Model").  The
input molecule is specified by coordinates and charges of each atom, either by
its AMBER format topology/coordinate files or the corresponding PQR. 


The nonpolar component of the solvation energy is calculated using the "surface
area + VDW" approach  described by Gallicchio and Levy.  (Journal of
computational Chemistry (2004), 25 479). In our implementation, the non-polar
terms are based directly on the R6 integrals already computed for the effective
radii needed for the polar part.

With either canonical GB or CHAGB,  the accuracy of the polar component of the
solvation energy can be improved via the ALPB model, which re-introduces the
correct dependence on dielectrics into the GB model.  The correction is
recommended in all cases except small molecules with decidedly non-spherical
topology (e.g. rings ) or structures that a topologically not singly-connected (e.g. 
two molecules not in contact with each other).  
For the ALPB reference see G. Sigalov, A. Fenley, and
A. Onufriev, "Analytical Linearized Poisson--Boltzman n Approach: Beyond the
Generalized Born Approximation", J. Chem. Phys., 124, 124 902 (2006)

The effective Born radii are estimated via the "R6" method as integrals over an
approximation to molecular sutface.  For a detailed description of the methods
and algorithms to compute the R6 effective radii please refer to: B. Aguilar,
R. Shadrach, and A. V. Onufriev, ``Reducing the Secondary Structure Bias in the
Generalized Born Model via R6 Effective Radii", J. Chem. Theory and Comput. 6,
3613 (2010), and B. Aguilar and A.V. Onufriev. ``Efficient Computation of the
Total Solvation Energy of Small Molecules via the R6 Generalized Born Model".
J. Chem. Theory and Comput., 8, 2404-2411 (2012).

Effective Born radii are computed by either of the two methods: NSR6 
(numerical) and  AR6 (analytical).

RECOMMENDED: NSR6 is the default for performing R6 calculation.  
The only required
input for NSR6 is the coordinate file in PQR format.
 
(NO LONGER RECOMMENDED: The version of AR6 in this suite 
is intended either for development or for fast estimation
of solvation energies of snapshots from very long MD trajectories. Requires both
pqr and topoly files.)

For accuracy benchmarks see B. Aguilar, R. Shadrach, and A. V. Onufriev,
``Reducing the Secondary Structure Bias in the Generalized Born Model via R6
Effective Radii", J. Chem. Theory and Comput. 6, 3613(2010), and B. Aguilar and
A.V. Onufriev. ``Efficient Computation of the Total Solvation Energy of Small
Molecules via the R6 Generalized Born Model". J. Chem. Theory and Comput., 8,
2404-2411 (2012). 

==================================================================
Recommended citations: Mukhopadhyay et. al., Sigalov et al., Aguilar et al., 
see above.  

==================================================================
Authors: Boris Aguilar, Abhishek Mukhapodhyay, Igor Tolokh, Alexey V. Onufriev. 

==================================================================
License: Free for any academic/educational/any other non-profit use. 
In which case you can do whatever you want with this code. No
restrictions of any kind, explicit or implied. 

For commercial purposes it is still 
free for up to a month (for evaluation).  
Afterwards, please contact the authors. 
==================================================================

Contact:  baguilar@cs.vt.edu 
          alexey@cs.vt.edu
          
==================================================================
System Requirements
==================================================================

   -- technically any Unix system
   -- gcc and g++ compilers 
   -- make
   -- MSMS (install first)

============
Installation
============

- Install MSMS. Make sure it is in your PATH. The
binary name should be set as "msms" 
MSMS is needed to compute surface triangulation, this program
can be downloaded from : http://mgltools.scripps.edu/packages/MSMS/

- do ./make

- the following program will be generated :
  bonrRadius        : program that computes the R6 Born radii, and solvation energies

- The following files are also provided with this ditribution:
  LennarJones.data  : parameters used in the computation of the nonpolar solvation energy 
  temp.neck         : values of Aij and Bij for Bondi radii set  
  

=================
Details of usage 
=================

If you do not have your structure in pqr format, which is
PDB + charges + atomic radii, you can generate it via  
H++ webserver sever (biophysics.cs.vt.edu/H++).  
It will generate you pqr files and AMBER topology files from pdb inputs. 
The protonation 
(charge) state of your structure will be properly set at the input pH. 
This method is probably best for macromolecules. You may also use 
ambpdb utility of AMBER to get PQR directly from inprd/prmtop files of 
your structure. For small molecules, ZAP9 radii set is recommended. 
See example. 

SPECIFY MODEL (GREEN FUNCTION)
-------------
                (The DEFAULT model for solvation free energy calculations is the Analytical
                Linear Poisson Boltzmann (ALPB). The other two models incorporated are the
                canonical GB without the ALPB correction
                and the new Charge Hydration Asymmetric Generalized Born (CHAGB).)

    -canonical : The polar solvation energy is calculated via the canonical GB model.

    -chagb : The polar solvation energy is calculated via CHAGB model, see "Introducing
             Charge Hydration asymmetry to the Generalized Born Model".
             The optimum parameters specific to CHAGB, including atomic radii,
             are provided in Mukhopadhyay et. al. The user defined intrinsic
             atomic radii in the pqr files are, by default, replaced  by the 
	     optimum radii set.  However the "-inradii inpqr" option (below) 
	     may be used for overriding the default, in which case the 	
             atomic radii are read from the input pqr file.
             JUST BE CAREFUL: Atomic radii suitable for CHAGB are quite different 
             from "standard".
	         
	     By default, the ALPB correction is used with the -chagb (recommended), 
	     but -chagb can also be run with the -canonical option (untested).



INPUT PARAMETERS
----------------
	-pqr filname : File containing your molecule structure in pqr format. Required. 
 

	-prmtop filename : File containing information on the topology of the molecule,
			   AMBER format. Used for the ZAP9 radii set or AR6.
			   Not needed otherwise.  

	-surf_density value : density parameter (vertices per A^2)) used for the
			      computation if the surface integral corresponding to
			      the R6 integral via MSMS. The default is 6, which is
			      a good trade-off between speed and accuracy. Use higher
			      density for slightly more accurate result, however the
			      calculation will become more expensive. Densities higher
			      than 20 are not recommended, especially for larger 
			      structures.

	-A value : Electrostatic size of the molecule. If not specified, the program
		   will automatically compute this value, in which case it is kept 
		   constant for every snaphot. The option is useful if the structure 
		   does not undergo large conformational changes and so you can save
		   a little bit of computational time by keeping the value fixed. The
		   savings are pretty small though.

	-B value : Empirical offset to be added to the inverse radii computation
		   at the very end. The default = 0.028 is optimal for agreement 
		   with the PB. However, B = 0 is recommended for small molecules
		   to improve agreement with explicit solvent free energies. 


        -gamma value : Surface tension parameter. Default=0.01 optimized 
			for small molecules.  

	-probe : probe radius for the calculation of the dielectric boundary (molecular 
		 surface). Default=1.4 for both the canonical and CHAGB. But see 
		 "Introducing Charge Hydration asymmetry to the Generalized Born Model" 
                 for details of how the dielectric boundary is constructed in the case
		 of CHAGB (all of the modifications are done internally by the code). 
 

        -epsin value : Internal (solute) dielectric constant, DEFAULT = 1
        
	-epsout value : External (solvent) dielectric constant, DEFAULT = 80
   
        -kappa value : Inverse of Debye-Huckel length, DEFAULT = 0.0; 
		       Use 0.1 to mimic physiological salt conditions. 

	-rinv_in filename : file containing the inverse of the effective 
                            Born radii per atom (second column of filname). 
                            The number of columns must be equal to the number 
                            of atoms in the pqr file. The first column 
                            especifies the atom index. Please refer to 
                            "PQR/beta.bondi.test.rinv" file for an example.

                        For CHAGB the inverse effective Born radii are "charge-symmetric", 
			just like in the canonical GB. With -chagb flag, the input
			effective radii will be 
                        further scaled internally by the CHA scaling factor (see
                        "Introducing Charge Hydration asymmetry to the Generalized
                         Born Model" for details).

        -inradii option : intrinsic atomic radii option can be zap9, bestgb, inpqr. 
			 Internally overrides the radii found in the input pqr.

            zap9      The atomic radii of the pqr file is replaced by the 
                      ZAP9 radii set and used in all the calculations, see 
                      Nicholls et. al., J. Med. Chem., 2008, 51 (4), pp 769–779. 
                      The option requires the AMBER topology file loaded with 
                      the -prmtop input command.
                      THIS RADII SET IS NOT RECOMMENDED FOR CHAGB
 
            bestgb    The atomic radii of the pqr file is replaced by the
                      a new radii set obtained by optimization agains TIP3P
                      solvation energies of small molecules, see
                      our paper "Introducing Charge Hydration asymmetry 
                      to the Generalized Born Model". Amber topology file 
                      is not required in this case. 
                      THIS RADII SET DOES NOT WORK WITH CHAGB 
 
            inpqr     Forces to use intrinsic radii from the pqr file for CHAGB.
		      Only useful if you want to use your own input atomic radii in CHAGB, 
		      and not have them automatically replaced by the optimal values. 
		      Expert only option. 


         
         CHAGB specific flags:
	     ---------------------
    
         -ROH value : The ROH value used with -chagb. 
			In CHAGB, the propensity
			of a water model being emulated is controlled by
			the water model-dependent parameter ROH, 
			which characterizes asymmetry of the molecule's 
			charge distribution with respect to its center
			(oxygen). 
                        The default is 0.58588 (Angstroms), ROH of TIP3P. 


         -tau value : The tau value used for the CHAGB model, 
                      tau is a positive constant
                      that controls the effective range of the 
                      neighboring charges (j) affecting the CHA of
                      atom (i). Optimized default = 1.47,  see Mukhopadhyay et. al.

EFFECTIVE BORN RADII
--------------------
(The RECOMMENDED DEFAULT for computing the effective Born radii is NSR6. 
	The other option is AR6 described below.)

	-AR6 : Flag to compute the inverse effective Born radii by analytic
           integration over approximate molecular volume. pqr and topology 
           file are mandatory inputs. This method is generally NOT RECOMMENDED 
	   unless speed of calculations is paramount for a large number of
           snapshots (of the same molecule). 
           Note that this is an alternative to the default NSR6, 
           and it has not been tested with CHAGB.
               
            AR6 specific flags:
	        -------------------

	    -neck filename : File to read in neck parameters.  Uses temp.neck by default.     
	
	    -prmtop filename : File containing information on the
                           topology of the molecule, AMBER format.
                           Used for the ZAP9 radii set and AR6 only.


             -beta1 value 
             -beta2 value : Parameters intended to improve analytical effective radii 
			    of buried atoms in larger structures. The default 
		            values are beta1 = 18.4377, beta2  = 313.1717. 

	     -Svdw value : Pairwise descreening weight. 
		           Optimized value 0.6211.

	     -Sneck value : Neck weight.  
		            Optimized value 0.4058.  Default.
             	

	     -neckRw : Probe radius for Necks. Default 1.4 A. You 
	               should delete temp.neck first if you are not
                       using the default.

	     -psf filename : File containing information on the topology 
		             of the molecule in PSF format.  
			     Will be overwritten by prmtop if 
			     this is available as well.
                             This option is experimental, tested on a
                             very limited number of cases. 

	     -prmtop filename : File containing information on the 
		      	        topology of the molecule, AMBER format. 
        				A must for for chunk computations in the AR6.  
			            This takes priority over the psf flag. 

	     -chunk filename depth : File to save chunk data 
			   	     to (or read it from if previously computed).  
				     depth is either a nonnegative integer 
				     or the string "infinity".
				     depth 0 is used to switch chunks off. 
				     Default depth = 3. If filename is not specified
				     "temp.chunk" will be used as chunk file name.
OUTPUT VALUES
-----------------

	-radii filename : print the inverse of the effective Born radii in filename

	-energy total_solv  : print the total solvation energy and its components	

	-energy polar :  print the polar component of the solvation energy

	-decomp filename : print the solvation energie for each pair of atoms 

	-energy nonpolar : print nonpolar component of the solvation energy
			   which is the sum of the solute-solvent 
			   VDW and "hydrophobic" parts.   

        -energy coulomb : print the Coulomb energy. In this version, all
			charge-charge interactions are computed, including 
			neareast neighbors. This is diffeernt from AMBER 
			convention where 1-2 and 1-3 interactions are 
			excluded. 

	-energy all : print the total electrostatic energy 
                      (nonpolar part excluded), coulomb energy +
                      and polar solvaion energy.

        -nolabel : this flag allows no to write out labels. The order in which 
                   the different energies appear is not changed.
  
	 (If the flag -energy is not especified, the program will 
	 compute and output all energies.)

	-decom_tot filename : print the total electrostatic energy 
         (polar solvation energy + coulomb energy)  between all pairs 
         of atoms. The first two columns especify the indices of the
         two atoms. Units in kcal/mol

	-decom_solv filename : print the polar componenet of the solvation energy
         between pairs of atoms. The output format is similar to the previous option. 


======================================================
Examples:

Compute all energies with labels using NSR6 (all parameters = defaults), using ALPB
./bornRadius -pqr PQR/beta.bondi.pqr -energy all

Compute the exact R6 (NSR6) inverse radii with surface density 6, using ALPB
./bornRadius -pqr PQR/beta.bondi.pqr -surf_density 6 -radii beta.bondi.rinv

Compute all energies without labels using AR6 with default parameters, using ALPB 
./bornRadius -pqr PQR/beta.bondi.pqr -prmtop PQR/beta.bondi.top -AR6 -energy all

Compute the analytic R6 (AR6) inverse effective radii with chunk depth 3, using specified 
parameters Sneck Svdw beta1 beta2:
./bornRadius -pqr PQR/beta.bondi.pqr -prmtop PQR/beta.bondi.top -chunk PQR/beta.bondi.d3.chunk 3 -AR6  -Sneck 0.4058 -Svdw 0.6211 -beta1 18.4377 -beta2 313.1717 -B 0.028 -radii beta.bondi.rinv

Compute solvation energies using ALPB for methylcyclopentane with specified probe radius 
./bornRadius -pqr PQR/methylcyclopentane.pqr -probe 1.4 -energy total_solv -B 0.0

Compute all energies using previously calculated effective Born radii:
./bornRadius -pqr PQR/beta.bondi.pqr -rinv_in PQR/beta.bondi.test.rinv

Compute all energies with labels using NSR6 (all parameters = defaults), and print the energies between pairs of atoms in the file "decom.dat"
./bornRadius -pqr PQR/beta.bondi.pqr -energy all -decom_tot decom.dat

Compute solvation energies for methylcyclopentane using the ZAP9 radii set
./bornRadius -pqr PQR/methylcyclopentane.pqr -prmtop PQR/methylcyclopentane.prmtop -inradii zap9 -energy total_solv -B 0.0

Compute solvation energies with canonical GB used for the polar part
./bornRadius -pqr PQR/methylcyclopentane.pqr -energy total_solv -canonical -B 0.0 

Compute solvation energies using the CHAGB solvation model (NEW) for the polar part
./bornRadius -pqr PQR/methylcyclopentane.pqr -energy total_solv -chagb -B 0.0 

Compute solvation energies using the CHAGB solvation model with specified inverse Born radii
./bornRadius -pqr PQR/methylcyclopentane.pqr -energy total_solv -chagb -rinv_in PQR/methycyclopentane.rinv

Compute polar solv  with CHAGB, and print the energies between pairs of atoms into the file "decom.dat"
./bornRadius -pqr PQR/methylcyclopentane.pqr -chagb  -energy polar  -decom_solv  decom.dat -B 0.0

(Check : awk '{sum = sum + $3} END{print sum }'  decom.dat  
should give you the computed polar solvation energy)

