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?
Initialization
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 
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