Bug in NBFix input routines in X-Plor
Mitch Miller
MITCH at BRAGG1.BCHS.UH.EDU
Fri Aug 21 18:06:37 EST 1998
Hi,
I found an error in the X-PLOR v3.851 code (& 3.843, although I have not
had a chance to check the MSI version yet). The problem is with the
non-bonded parameter input file parmio.s When non-bonded parameters are
input via the NBFix statement, the value of (Rmin/2) is stored in the
van der Waals radius array, CNBVR, rather than Rmin. This error has
very serious consequences for the repel vdW potential and more minor
effects on the Lennard-Jones potential.
In summary, NBFix incorrectly stores Rmin/2 in the vdW array CNBVR.
Line 972 of parmio.s needs to be changed from
CNBVR(I1,J1)=SIGI*TWO56
to
CNBVR(I1,J1)=TWO*SIGI*TWO56
so that Rmin is stored in CNBVR.
This only has an effect if NBFix is used to input parameters. This
causes an error in 2 of the 4 test cases outlined below.
Without this fix in X-PLOR 3.851:
1. vdW energy and forces using a Lennard-Jones potential WITHOUT nbfix
works correctly
2. vdW energy and forces using the repel WITHOUT nbfix works correctly
3. vdW energy and forces using a Lennard-Jones potential WITH nbfix
works many cases. The NBFix statement properly stores the A and B
parameters for the L-J potential which is used for calculation of the
vdW energy. The value of Rmin from the array CNBVR is only used to
calculate the truncation radius. Since NBFix parameter input is storing
Rmin/2 instead of Rmin, the truncation radius (INHIbit * Rmin) is cut in
half. So the energy term is correct EXCEPT when there is a very close
contact. Specifically, truncation of the energy functions does not
occur until R < (Rmin/2)*INHIbit, so the energy is wrong whenever
R < Rmin* INHIbit
Since the default value for INHIbit is 0.25, this only changes the vdW
energy when two atoms are closer than 1/4 of the sum of their van der
Waals radii (Rmin) or around 0.75 Angstroms in most cases.
4. The most egregious errors occur when the Repel van der Waals energy
function [fVDW(R) = 16*(Rmin - R)^4] is used WITH nbfix. Since Rmin/2
is stored in the CNBVR array for all NBFix input interactions, the
effective van der Waals radii are cut in half! The vdW energy is zero
until the atoms are Rmin/2 apart. This allows atoms to get much closer
than they should.
Again this only affects interactions modified with the NBFix statement.
For example, NBFIX is used in param19.sol to modify water-water interactions.
-Mitch
========================================================================
| Mitchell D. Miller | mitch at bragg1.bchs.uh.edu |
| Department of Biochemistry | Office: (713) 743-8371 |
| University of Houston | X-ray Lab: (713) 743-8387 |
| Houston, TX 77204-5934 | FAX: (713) 743-8373 |
========================================================================
For those who are interested, here are the details about what I found:
I found this error when trying to modify the non-bonded potential
between a Cys-SG and Zn ion with an NBFIX statement. In testing to see
if I input the parameters correctly, I found that the repel vdW potential
was zero until R = Rmin / 2.
However, if I used the L-J vdW potential, the function had a minimum of
-eps at R = Rmin.
To understand these differences, I did a series of tests of vdW energy
as a function of distance for several different 2 atom models. And then
looked into the code to try to find an explanation for the behavior I
found.
Careful inspection of the L-J potential plotted as a function of
distance for a 2 atom model also revealed that the function was not
quite the same when parameters were input via NBFix. That is, the
maximum of the potential function was at a shorter distance (Rmin/8
instead of Rmin/4, which in turn makes the maximum much higher).
(In X-PLOR 3.1, CNBVR was a linear array with one entry / atom to store
its vdW radii. In X-PLOR 3.857, this array was made 2D with ability to
patch Rmin for specific interactions via NBFix.)
The following relations are defined in the X-PLOR 3.1 manual (p.76),
Source code (parmio.s), and / or comments in the protein_rep.param
parameter file.
Rmin = sigma * 2^(1/6)
Emin = -eps
A = 4 * sigma^12 * eps
B = 4 * sigma^6 * eps
TWO56 = (1/2)^(5/6)
When parameters are input via the NONBonded statement, Rmin is placed in
the CNBVR array. The NONBonded statement creates the diagonal elemnets
in the array, and the off-diagonal elements are filled in via the
combination rules. However, the when the NBFIx statement is used to
modify the off-diagonal terms, ( Rmin / 2 ) is placed in this same CNBVR
array.
Line 746 (NBON parameter input: values from the NBON statement are read
in as sigma and eps. The diagonal elements of the CNBVR array are
correctly assigns Rmin to the arrary)
CNBVR(NCN,NCN)=TWO*SIG*TWO56
The when the cross terms are calculated the values are correct according
to the combination rules in the manual. (see p. 76 of X-PLOR 3.1 manual)
CNBVR(I,NCN)=(SIGI+SIG)*TWO56
(Recall that sigma is Rmin*(1/2)^1/6 which when multiplied by
(1/2)^(5/6) gives Rmin/2, so that (SigI + SigJ)*TWO56 = RminI/2 + RminJ/2
In the NBfix statement, (see parmio.s lines 891 ff)A and B are read in
(see line 960 and manual p. 35)and placed in the arrays CNBA and CNBB
CNBA(I1,J1)=A
CNBB(I1,J1)=B
CNBA(J1,I1)=CNBA(I1,J1)
CNBB(J1,I1)=CNBB(I1,J1)
The subroutine LJESAB is called to extract EPS & SIGMA from A & B and
then add these terms to the CNBVR array. (line 971) (Note the error is
line 972)
CALL LJESAB(EPSI,SIGI,A,B)
CNBVR(I1,J1)=SIGI*TWO56
CNBVR(J1,I1)=CNBVR(I1,J1)
The subroutine LJESAB returns EPS and Sigma given A and B
(see line 1207)
EPSI = 1/4 * B * B / A = 1/4 ( 2 * Rmin^6 * eps)^2 / (Rmin^12 * eps)
= 4 * Rmin^12 * eps^2 / (4 * Rmin^12 * eps) = eps
SIGI = (A/B)^(1/6) = ((Rmin^12 * eps)/(2 * Rmin^6 * eps))^(1/6)
= (Rmin^6 / 2)^(1/6) = Rmin * (1/2)^(1/6)
The value which is stored in CNBVR array is
[Recall: TWO56=ONE/TWO**(FIVE/SIX) ]
CNBVR = SIGI * (1/2)^(5/6) = Rmin * (1/2)^(1/6) * (1/2)^(5/6) = Rmin / 2
However, the value stored in the CNBVR array should be Rmin, thus line
972 should read:
CNBVR(I1,J1)=TWO*SIGI*TWO56
or CNBVR = 2 * [Rmin * (1/2)^(1/6)] * (1/2)^(5/6) = Rmin
This explains why I did not see any vdW energy with the repel function
until I moved the two atoms closer than (Rmin / 2). One can compensate
for this in the repel function by using an NBFix statement where A and
B are calculated using 2*Rmin. But the code should be fixed since using
the same 2*Rmin with the L-J vdW potential will give very wrong answers.
========================================================================
I would be happy to send the X-plor input script file I used to
verify the error to anyone who would like it.
========================================================================
-Mitch
More information about the X-plor
mailing list