OPLS-AA naphthalene

I’m using the OPLS-AA force field for simulating a system of 100 naphthalene molecules and have observed unexpectedly high potential energy values. I have set up my system with rigid naphthalene molecules and am aware that minimization is not applicable for rigid bodies. Could the high potential energy be a result of the way I’ve configured the force field parameters, or might it relate to the treatment of naphthalene as a rigid body? What steps can I take to troubleshoot and resolve this issue in my simulation?


units real
atom_style full
boundary p p p
special_bonds lj/coul 0.0 0.0 0.5
pair_style lj/cut/coul/long 12.0
kspace_style pppm 1.0e-4

Read in initial configuration

read_data 10naphthalene.data

Define Lennard-Jones coefficients

pair_coeff 1 1 0.070 3.550
pair_coeff 1 2 0.055 3.000
pair_coeff 2 2 0.030 2.420

Define group for

group naphthalene type 1 2

Output settings

compute mype all pe
compute mymsd all msd
compute myTemp all temp
compute mypre all pressure myTemp

thermo 1000
thermo_style custom step time temp density press vol pe evdwl ecoul

dump 2 all atom 10 naphthalene_1.dump

Apply rigid fix to naphthalene molecules

fix myrigid naphthalene rigid/nvt molecule temp 300 300 100

Set the timestep

timestep 1.0

Run the simulation

run 100000

PotEng         E_vdwl         E_coul 

44719076 44720091 -724.14662

First off, please have a look at this post to learn how to correct quote text in the forum so that special characters are not interpreted as markup.

Sometimes this can be worked around by temporarily using extremely stiff force constants for bonds/angles/dihedrals/impropers to keep the molecules “near rigid”.

Could be by it is not likely.

very unlikely

Most likely candidate for a problem is your data file, e.g. the box dimensions are incorrect not large enough or molecules are overlapping. This can be tested by first creating a system with few molecules (like 1 or 2) and checking what the energy per molecule would be like. Then you can gradually increase to number of molecules until you see a significant difference. That is where you need to look at the file more carefully.

Another option would be to use a different strategy to create the system, e.g. through using a molecule file and then create_atoms from within LAMMPS with the “overlap” parameter set to guarantee a minimum distance between molecules.

Here are some tips for your simulation.

Firstly, you can check the density through thermo_style to determine if the simulation box is too crowded.

Secondly, try to turn off the kspace by kspace_style none (also change the pair_style) at the beginning of the simulation. After a fast relaxation, add the kspace through pair_modify and kspace_style pppm 1.0e-4. This could save a lot of time for a large amount of atoms.

Hope you will find these useful. My first reply :slight_smile:

1 Like

To me a high E_vdwl points to overlapping atoms in a crowded box. You are also using a NVT integrator, so the potential energy should not decrease significantly. As @leo-lyy pointed out, printing the density will easily verify this hypothesis.
Also, please quote the code. You are not sending an SMS.

1 Like