[gmx-users] Re: WCA and free energy calculations

jtomlins at purdue.edu jtomlins at purdue.edu
Mon Jun 25 16:31:09 CEST 2007


> Hi,
> 
> I'm CCing the users list on this since it might be of general interest.
> 
> >  I can't begin to tell you already how much help you've given me with your
> >  gromacs user-list posts and free energy tutorial.
> 
> Glad I can help!
> 
> > I was searching the Gromacs user list and came across a thread by you
> > titled:tabulated nonbonded interactions and free energy calculations.
> > 
> > I'm trying to do exactly the same thing that you described, using the WCA
> > LJ potential and free energy calculations.  No one ever responded to this
> > thread, but I was wondering if you figured out the answer?  Also, just to 
> > confirm that I'm not spinning my wheels did you get the WCA LJ potential to 
> > work using the user specified potential function? 
> 
> So, no one ever responded to that thread, as you noticed. But my end
> conclusion was that the tabulated nonbonded interactions know nothing
> about lambda. So while one could implement the WCA separation using
> tabulated nonbonded interactions, one couldn't use it to do free
> energy calculations in a straightforward way, at least not
> simultaneously with the soft core potentials I wanted to use (where
> the effective radius is modified as a function of lambda). So I've put
> it on the back burner for now. I guess you could in principle get WCA
> to work with linear lambda scaing, but there is a lot of work in the
> literature indicating that linear scaling is bad for transformations
> involving insertion or deletion of LJ sites, which is presumably the
> point of using the WCA separation in the first place.
> 
> Anyway -- I think that all means that one would actually need to
> implement WCA into GROMACS to be able to use it for free energy
> calculations.
> 
> I hope this makes sense. Please let me know if any of it needs further
> explanation. Developers are welcome to weigh in if I'm going wrong
> here anywhere...
> 
> Thanks,
> David
(Sorry so long, wanted to get all the info in.)
I was thinking of another way to possibly implement WCA of a sorts into GROMACS 
without adding it into the code and since I'm still an baby in simulation work I 
wanted some comments on where I might run into problems...
I'm using the OPLS potential with TIP4P water and other small molecules, GROMACS 
v. 3.3.1.

The thought is that on page 97 of the v. 3.3 manual it says that the LJ 
parameters can be interpolated during FE calculations.  For the dummy atoms I 
leave the charge at full and change the elpison value to zero, do a cutoff at 
rmin or 2^(1/6)*sigma, and then do a shift of eplison.  Is this equilivant to 
removing the attractive term of the LJ potential?  So any configurations found 
would be due to the replusive term only.  Thus to find the attractive part of 
the potential I can do a mdrun -rerun and change the cutoff back to 9A and 
remove the shift then get my full LJ, so then attractive part is just the full 
LJ-replusive LJ.  

One problem I know that I will have with this is that each atom has a different 
eplison value, so then what eplison value would be approiate to use?

I've tried to implement this idea with some sucess, using the average eplison.  
For example a "normal" simulation of methane in water (no FEP, Coul cutoff of 
9A, vdw switch 8.5A and cutoff at 9A) has
Coul=-0.078  LJ=-13.083
When I use this method for lambda=0 (which by my definitions in the topology 
should have full LJ) the mdrun -rerun I find
Coul=-0.180  LJ=-7.601
This seems kinda logical because the cutoff was shorter so there is more 
emphasis put on the Coulomb term and decreased LJ.  When I do lambda=0.2 I get
Coul=-0.190  LJ=-7.371
and for lambda=0.4 I get
Coul=-0.213  LJ=-6.924
Again, the numbers are mostly changing the way I expect them to, I'm just 
worried that there is a major problem I'm missing since the energies are so 
different from that of a "normal" run.

Ideally the goal is to separate the different terms of the LJ potential so I can 
identify what interaction energies are due to the attractive and repulsive terms 
of the LJ. 

Thanks,
Jill

Jill Tomlinson-Phillips
Graduate Student
Department of Chemistry
Purdue University
560 Oval Drive Box 791
West Lafayette, IN 47906
jtomlins at purdue.edu



More information about the gromacs.org_gmx-users mailing list