[gmx-users] truncated LJ potential

[Mark Abraham]

On 27/01/2011 12:11 PM, Makoto Yoneya wrote:
> Dear Berk and all:
>
> I'd tried to rewrite the routine
>
>>>> src/gmxlib/nonbonded/nb_generic.c
> to modify the LJ potential to the shifted and truncated one.
>
> First, I'd add the new switch(ivdw) as case 4, but when I'd tried;
> [ defaults ]
> ; nbfunc comb-rule   gen-pairs   fudgeLJ  fudgeQQ
>    4      1     no    1.0   1.0,
> I'd got,
> ERROR 1 [file topol.top, line 8]:
>    Invalid nonbond function selector '4' using LJ.

So you have to update the machinery that parses the .top to recognise 
that the value of 4 is now legal.

> Then, I'd rewrite the standard switch(ivdw) case 1 as follows
> (with real tc12c5 and int factor).
>
> #ifdef TRUNCATED_LJ
> 	tc12c6 = 2.0*c12/c6;
> 	factor = ceil(floor(tc12c6*rinvsix)/(tc12c6*rinvsix));
> 	fscal += factor*(12.0*Vvdw_rep-6.0*Vvdw_disp)*rinvsq;
> 	Vvdwtot = Vvdwtot+factor*(Vvdw_rep-Vvdw_disp+0.5*c6/tc12c6);
> #else
> 	fscal += (12.0*Vvdw_rep-6.0*Vvdw_disp)*rinvsq;
> 	Vvdwtot = Vvdwtot+Vvdw_rep-Vvdw_disp;
> #endif
>
> However with this modification (recompile with #define TRUNCATED_LJ),
> the simulation result was exactly the same.
>
> I'd appreciate to tell me what kind of mistake I'd done?

Did you set the environment variable to actually call the generic 
nonbonded lists?

Mark

> Makoto Yoneya, Dr.
> AIST
> http://staff.aist.go.jp/makoto-yoneya/
>
>> Yes.
>>
>> The pow function is expensive though. The code will run much faster
>> if you can use rinvsix, such as check for  2*rinvsix>  c6/c12.
>> (I forgot the factor 2 in my previous mail).
>>
>> Berk
>>
>> From: makoto-yoneya at aist.go.jp
>> To: gmx-users at gromacs.org
>> Date: Tue, 11 Jan 2011 10:10:56 +0900
>> Subject: [gmx-users] truncated LJ potential
>>
>> Dear Berk:
>>
>> Thanks again for the further reply.
>>
>>>> The LJ potential and force code in the above looks like in the c6-c12
>> form
>>>> not in epsilon-sigma one.
>>>> The LJ potential modification I'd like to try is based on the epsilon-
>>>> sigma form and the mixing rule is the Lorents-Bertelot's one.
>>>> Could you kindly tell me the LJ potential and force routine in the
>> epsilon-
>>>> Sigma form.
>>> There is no such code.
>>>
>>> You can simply check for rinvsix>  c6/c12
>> Is it means that the LJ potential in epsion-sigma form (with the
>> Lorents-Bertelot
>> mixing rule) is evaluated after the coversion into c6-c12 form in GROMACS?
>> Then, may I evaluate the modified LJ potential:
>>
>>>> V(r)
>>>>     = 4*epsilon*{ (sigma/r)^(12) - (sigma/r)^6 + (1/4) } for r<=
>>>> 2^(1/6)*sigma
>>>>     = 0 for r>   2^(1/6)*sigma
>> with translated into the equivalent c6-c12 form:
>>
>> V(r)
>>     = (c12/r)^(12) - (c6/r)^6 + (c6/2)*(c6/2*c12) for r<= (2*c12/c6)^(1/6)
>>     = 0 for r>  (2*c12/c6)^(1/6)
>>
>> in the following routines.
>>>> You can also set the environment variable nb_generic.c and modify
>>>> src/gmxlib/nonbonded/nb_generic.c, but might lead to somewhat
>>>> slower simulations.
>> If my understanding in the above would correct, I'll try that.