Bug in NBFix input routines in X-Plor

[Mitch Miller]

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