[gmx-developers] Energy and Virial long range corrections

[Juan S. Medina Alvarez]

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


       
______________________________________________ 
¿Con Mascota por primera vez? Sé un mejor Amigo. Entra en Yahoo! Respuestas http://es.answers.yahoo.com/info/welcome