Hi,
I have found and I believe solved a bug causing the dispersion term to
be incorrect when the radiopt=0 option is set. This applies to pbsa but
also to sander and MMPBSA.py as the latter use the pbsa code. A patch to
be applied on AmberTools19 is provided.
I previously reported (
http://archive.ambermd.org/201905/0152.html and
201906/0090.html) that the dispersion term is different between radiopt=
0 and 1, which is a problem, as radiopt only concerns the PB radii.
I further investigated the problem on a single atom (Na+), for which the
dispersion term can be calculated analytically.
E_disp = 4*pi*rhow*(A/(9*sigma^9) - B/(3*sigma^3)) when r <= sigma
with r = rmin_Na+ + sprob
rmin_Na+ = 1.8680
eps_Na+ = 0.00277
rmin_w = 1.7683
eps_w = 0.1520
sprob = 0.557
rhow = 1.129*0.03333
I found that the correct value (-0.2932) is obtained with radiopt=1,
while a value of -0.4005 is obtained with radiopt=0.
(The radiopt=0 value corresponds to an inconsistent situation where
eps_Na+ is set to B_Na+/(256*rpb_Na+^6), where rpb_Na+ is the radius
read from the prmtop (e.g. 1.5 for mbondi radii, giving an incorrect
eps_Na+ of 0.00517), while rmin_Na+ is correct.)
The bug is caused by the following line in the "if ( radiopt == 0 )
then" section of AmberTools/src/pbsa/pb_init.F90:
mdsig = rin ! for np
I think there is no reason to redefine mdsig here, as it is not used for
PB and the radiopt option should not affect NP calculations.
The problem could have been avoided though as sigma is recalculated
for each atom from the prmtop A and B coefficients in the dispersion
calculation performed in AmberTools/src/pbsa/np_force.F90:
mdsig_iatm = (cn1(ic)/cn2(ic))**SIXTH/2
Unfortunately, on the following line, epsilon is calculated from the
mdsig array (incorrectly set to rin as seen above) instead of from the
just calculated mdsig_iatm:
epsln_iatm = cn2(ic)/(256.0d0*mdsig(iatm)**6)
Please find attached a patch solving this problem.
The problem can be reproduced as follows:
cat > pbsa.radiopt$radiopt.in <<EOF
PBSADI radiopt=$radiopt
&cntrl
ipb=2,
/
&pb
radiopt=$radiopt,
/
EOF
pbsa -O -i pbsa.radiopt$radiopt.in -o pbsa.radiopt$radiopt.out -p
Na+.top -c Na+.crd
Na+ .top and .crd files are attached.
--
Thomas Gaillard
Professeur assistant
Laboratoire de Biochimie
Ecole Polytechnique
91128 Palaiseau cedex
tel: +33 1 69 33 48 62
fax: +33 1 69 33 49 09
thomas.gaillard.polytechnique.edu
http://thomasgaillard.fr
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Jun 17 2019 - 06:30:04 PDT
