[gmx-users] truncated LJ potential

[Makoto Yoneya]

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.

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?

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.