[gmx-developers] Energy and Virial long range corrections

[Berk Hess]

Juan S. Medina Alvarez wrote:

>Hello everybody
> 
> Firstly a brief summary:
>
>Juan S. Medina Alvarez wrote:
>
>  
>
>>>Dear Sir:
>>>
>>> I'm a Gromacs user, in Spain (Canary Islands). My
>>>      
>>>
>calculations are 
>  
>
>>>with water models that use intermolecular
>>>      
>>>
>interactions different from
> 
>  
>
>>>Lennard-Jones type. Specifically I am using a MCY
>>>      
>>>
>(Matsuoka Clementi 
>  
>
>>>Yoshimine) potential for the O-H, O-O, ...
>>>      
>>>
>interactions. Of course 
>  
>
>>>through the use of tables in your program. The
>>>      
>>>
>coefficients in the 
>  
>
>>>force field file (ff*nb.itp)  are, in some cases,
>>>      
>>>
>almost 10 orders of
> 
>  
>
>>>magnitude greater than the coefficients in the
>>>      
>>>
>ffgmxnb.itp force 
>  
>
>>>field, but there are not problems in the
>>>      
>>>
>calculations obtaining total
> 
>  
>
>>>energies (I compare them to some calculations
>>>      
>>>
>previously done by me 
>  
>
>>>using Moldy -another MD software-). Nevertheless
>>>      
>>>
>when I try to do the
> 
>  
>
>>>same calculations with the option "DispCorr:
>>>      
>>>
>EnerPres" in the *.mdp 
>  
>
>>>grompp configuration file I obtain very excessive
>>>      
>>>
>values of energy
> and 
>  
>
>>>pressure. For example with "DispCorr: no" the value
>>>      
>>>
>of total energy
> is 
>  
>
>>>in average -5714.39 kJ/mol while with "DispCorr:
>>>      
>>>
>EnerPres" the mean 
>  
>
>>>value is about -7.88564e+09 [sic] kJ/mol; similar
>>>      
>>>
>and unconfortable 
>  
>
>>>results occurs with pressure.
>>> I have thought that the reason for this behavior
>>>      
>>>
>is this: Gromacs 
>  
>
>>>calculates the Dispersion correction to energy and
>>>      
>>>
>pressure using the
> 
>  
>
>>>Lennard-Jones model of potential instead of
>>>      
>>>
>manipulating the tables 
>  
>
>>>for the MCY potential. (The clue: the ratio among
>>>      
>>>
>the C6 coefficients
> 
>  
>
>>>of the two models). If I am right I would like to
>>>      
>>>
>modify those lines 
>  
>
>>>of the Gromacs code which deal with the long range
>>>      
>>>
>correction in 
>  
>
>>>Energy and Virial terms. Unfortunately the code is
>>>      
>>>
>too big to do that
> 
>  
>
>>>in the dark. So here is the motivation to write
>>>      
>>>
>you. Could you point 
>  
>
>>>me out which parts of the code I might to review to
>>>      
>>>
>correct this if 
>  
>
>>>I'm really in the good way?
>>>
>>> Yours faithfully
>>>
>>> Juan S. Medina Álvarez, Ph.D.
>>> e-mail: tlazcala at yahoo.es
>>> e-mail: jmedina at iusiani.ulpgc.es
>>>      
>>>
>
>Berk Hess wrote (answer):
>
>  
>
>>I have no clue how these potentials look like.
>>You also have not told how your tables are set up.
>>    
>>
>
>  
>
>>This question is a typical question for the
>>    
>>
>gmx-developers mailing list.
>  
>
>>Could you mail a short description of the problem and
>>the details of the potentials to the gmx-developers
>>    
>>
>mailing list? I will answer there, so others might
>also profit from the discussion.
>
>  
>
>>Berk.
>>    
>>
>-------------------------------------------------------
>
> Dear Sir:
>
> A MCY potential is an extended form of the Morse
>Potential. The precise analytical expression is:
>V_MCY(r)=A exp(-br) - C exp(-dr).
>( If A=D exp(2aR), C=2D exp(aR), and b=2a, d=a we
>recover the Morse Potential, i.e.
>V_morse(r)=D{(1-exp(-a(r-R)))^2-1} )
>
> We have a water model with four sites (O,H,H,I) and
>we calculated five tables (*.xvg) for the O-O, O-H,
>O-I, H-H, H-I, interactions and an additional table
>for the I-I remaning interaction (null one). Following
>the Gromacs manual instructions the first column is
>for the value of r, the second and the third columns
>are for interpolating the Coulomb's Potential; the
>fourth and the fifth ones for the dispersion term (
>exp(-dr) and d^2*exp(-dr) ) and, finally the sixth and
>the seventh column for the repulsion term ( exp(-br)
>and the proper derivative). The coefficients are
>consigned
>in a customized force field file of type ff*nb.itp as
>C6=C and C12=A, (combinational rule 1).
> Really the first time I ran the model I didn't
>realize 
>that the dispersion correction was surely done
>following a dispersion term of r^-6 order which is
>totally unsuitable for this kind of potentials.
> I would like to correct the long range energy and
>pressure/virial; could anyone point me out where I
>must   touch the code up to do this?
>
> Yours 
>
> Juan S. Medina Álvarez
> e-mail: tlazcala at yahoo.es
> e-mail: jmedina at iusiani.ulpgc.es
>
>
>  
>
I think that the only code you have to change is in src/mdlib/sim_util.c,
lines 627-632:
      rc3  = fr->rvdw*fr->rvdw*fr->rvdw;
      rc9  = rc3*rc3*rc3;
      eners[0] += -4.0*M_PI/(3.0*rc3);
      eners[1] +=  4.0*M_PI/(9.0*rc9);
      virs[0]  +=  8.0*M_PI/rc3;
      virs[1]  += -16.0*M_PI/(3.0*rc9);

These currently give the energy and virial correction for -r^-6 and r^-12.
You should change this to the correction for your exponentials.
This will only work when b and d are the same for all potentials.

Berk.