Using WCA potential in gromacs

User discussions
1

GROMACS version: 2019
GROMACS modification: Yes/No
Here post your question

Hi, I am trying to use Weeks-Chandler-Anderson potential in gromacs. I have a system containing fullerenes in water. I want to treat the fullerenes as hard-spheres. So, that ful-ful,ful-water interactions are computed using WCA potential. I have gone through the gromacs manual about how to use tabulated potential. It says to define function as: u®=f®+Cj®+Ah®, where f®-for electrostatics,j®-for dispersion and h®-for repulsive interactions. A,C parameters are going to be defined in topology file as A=4epsilonsigma^12,C=4epsilonsigma^6, depending on atom types.

WCA=4*epsilon((sigma/r)^12 - (sigma/r)^6) + epsilon

Now for a typical LJ type function, we can easily use j®=-1/r^6,h®=1/r^12; but for WCA there is an additional epsilon(potential well depth) parameter, which is typically defined in topology file. How to include that additional epsilon in the tabular form, what would be h® and j® if we use WCA model? I am struggling with functional form that must be used in order to compute the tabular form of potential energy function for WCA potential. Any help will be much appreciated thank you

2

Hi,
The potential function is defined as
V(rij) = (qiqj/4πϵ0) * f(rij) + C6 * g(rij) + C12* h(rij)

where C6, C12 and the charges are read from the topology.
To generate the tabulated potential for WCA potential
I suggest is to assign a value of 1 to C6 and C12 (in the topology), then you have
g(rij) + h(rij) = 4*epsilon((sigma/r)^12 - (sigma/r)^6) + epsilon
now you can decide which function assign to g(rij) and to h(rij).

It is always good to run a small test to check that the potential is implemented correct.

Best regards
Alessandra

3

Thank you for your reply Alessandra.

The functional form for g(rij)+h(rij) will still contain two additional variables: sigma and epsilon. My system consists many fullerenes in water. Lets say I want to model fullerene molecules as hard spheres. Now for FUL-FUL interaction I can compute the table for potential energy function by putting sigma and epsilon values of Carbon atom prior to performing the simulation. However, in case of FUL-SOL interaction, can we proceed with the same table? Here we will have non-bonding interactions between C-atoms of FUL and O-atoms of SOL. Now in order to compute table for this interaction how do we compute corresponding sigma and epsilon? As far as I know for using tabulated potential combination rule is set to 1. If you could this, that would be really helpful, thank you

4

Hi,
I suggest you generate a table for each interaction type.
For a system with two particles types (C and O), you can generate 3 tables, each one for each interactions (C-C, C-O, O-O).
Does this answer your question?
Alessandra

5

Thank you for your response. For similar atom types like C-C or O-O interactions, we can compute the corresponding table using sigma and epsilon values of C or O respectively. But, for C-O interaction, how do we calculate the corresponding tabular energy file? I mean what value of sigma and epsilon must be used in that case?

6

Hi
That is independent from the tabulated potential.
This depends on which rules applies to your potential. If you use a combination rule to combine sigma and epsilon, then you use the same rule to generate the C-O potential.
I hope it helps
Best regards
Alessandra

7

I was thinking we can calculate corresponding sigma and epsilon values for C-O interactions in the following way:

according the combination rule=1, we get C6 and C12 for C-O interactions i.e. C as: C6(CO)=sqrt(C6©C6(O)) and C12(CO)=sqrt(C12©C12(O)); and C6=4epsilonsigma^6 and C12=4epsilonsigma^12; thus from these two equations we can get the sigma and epsilon values for C-O interactions.

Is the above mentioned a correct wayout?