[AMBER] pbsa EDISPER with radiopt=0 vs radiopt=1

From: Paulo Costa <pjfmcosta.gmail.com>
Date: Thu, 23 Feb 2017 14:57:40 +0000

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

Sample single point PB calculation
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

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
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
Received on Thu Feb 23 2017 - 07:00:03 PST
Custom Search