[gmx-users] Implementing WCA with OPLS-aa

[jtomlins at purdue.edu]

Hi,
I've been trying to implement WCA potential by the user defined potential
function as Berk suggested in the following thread:
http://www.gromacs.org/component/option,com_wrapper/Itemid,84/

I think I have it working for a simple case of argon in tip4p water, but am
having trouble going to the more complex solute of methane.

I'm using OPLS-aa force field, and Gromacs 3.3.1.  I've made my user defined
tables by specifying the entire potential as either the g(r) or h(r), where g(r)
is the entire attractive part of wca and h(r) is the repulsive part of wca. 
Thus my C6/C12 in my topology are just 1 or some other scaled value between 0
and 1.  I also changed my combination rule to 1.

1) Only the oxygen in tip4p water has sigma/eplison values, so I just used the
sigma/eplison of oxygen for the sol interaction.  In the topology I figured that
hydrogen and dummy point both had C6/C12 of zero so they would only
contribute energies of 0.  Is this reasonable? 
I did test this and all of the calculated energies are almost exactly
equivalent to using the cutoff option and sigma/eplison for calculation of the
vdw.  If I made some error here and the numbers just happened to still work, I
figure it might help answer my next question.  

2) Next I wanted to do a more complex solute, methane. Since both the carbon and
hydrogen atoms have different sigma/eplison then my thought was that I need a
table for each interaction (c-sol, h-sol, etc.).  Then the total interaction
potential energy of methane with water would be LJ(c-sol)+LJ(h-sol). As
suggested in the discussion:
http://www.gromacs.org/component/option,com_wrapper/Itemid,165/

I made an index file with the system, sol, and
[ CA ]
   1 
[ HA ]
   2    3    4    5

And I've defined the energy tables and groups in my .mdp as follows:
energygrps          =  CA HA SOL
energygrp_tables    =  CA SOL HA SOL SOL SOL

However when I try to do this simulation I get the grompp error:
Fatal error:
atoms 1 and 2 in charge group 1 are in different energy groups

Atoms 1 and 2 are my C and H in the methane molecule.  I can see how using a
united atom potential where I just have one sigma/eplison for each charge group
might solve this problem, but I'd really like to do this with the all atom
potential.  Any suggestions as to where I'm going wrong?

Thanks a ton,
Jill

Ben-Amotz group
Department of Chemistry
Purdue University
West Lafayette, IN