Dear AMBER users
I'm sorry for the long post ...
I'm performing some test calculations with the pbsa solver (AmberTools16)
to get the solvation free energies of small molecules. These small
molecules are built with antechamber using GAFF atom types.
I started with a simple single point calculation using the example from the
manual
Sample single point PB calculation
&cntrl
/
inp=2
&pb
npbverb=1, istrng=150, fillratio=1.5, saopt=1, radiopt=1
/
I got the non-polar contribution modeled by 2 terms, cavity + dispersion
(inp=2) and I since GAFF atom types are employed, radii from the prmtop
file will be used according to the manual (radipot=1). This calculation
yielded
EPB = -12.3409 | ECAVITY= 16.4597 | EDISPER = -15.1553
So, let's imagine I want to optimize the radii of the elements of my
molecule. The manual states that one should change the RADII section of the
topology and read those radii via RADIOPT=0. But if I use the exact same
topology as before, the radii should be the same. Therefore, I repeated the
calculation with
Sample single point PB calculation
&cntrl
/
inp=2
&pb
npbverb=1, istrng=150, fillratio=1.5, saopt=1, radiopt=0
/
and I got
EPB = -12.3409, ECAVITY= 16.4597, EDISPER = -17.0143
the value for EDISPER is completely different. Using the verbose option in
both runs, I confirmed in the output that all radii (PB and NP) are exactly
the same for both runs.
Besides EDISPER, the outputs differ in the following
< Use Tan, Yang, and Luo optimized cavity radii definition
---
> Use cavity radii in the prmtop file
Indeed, if one changes the RADII section in the topology for one or more
atoms, this only affects the EPB values as these are used for for the
electrostatic contribution calculation. Apparently the NP radii are always
the Rmin values of the FF. If so, it makes sense that the ECAVITY term is
equal regardless if we use radiopt = 0 or 1 since, in both cases, Rmin
values are used. On the other hand, why are EDISPER values different?
If we optimize a radii for an element (via the topology), we should run a
radiopt = 0 but how to compare it with the non-optimized raddii? The
default calculations are performed with radiopt = 1 but that will change te
contribution from EDISPER as well ...
If one wishes to optimize the NP radii, used for the calculation of the
non-polar contributions, we have to "hack" the source code? Or is there any
other way to do it (e.g. via topology)?
Again, sorry for the long post and thanks in advance for your help!
Paulo Costa
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Feb 24 2017 - 02:00:03 PST
