LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Missing atoms OPLS force fied
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Missing atoms OPLS force fied

From: Axel Kohlmeyer <akohlmey@...24...>
Date: Fri, 25 May 2018 18:39:51 -0400

On Thu, May 24, 2018 at 10:56 AM, DUNDAR YILMAZ <duy42@...122...> wrote:

I am trying to simulate polymer systems using OPLS force field.
The system consists of Polyethylene chains and peroxide molecules on top of them.
I build the structure PE and peroxide parts separately and merge them using read data command twice.
Everything works fine. I minimize the structure without any problem. The system is a slab, periodic in x-y directions only. I defined a LJ wall at top and bottom of the structure using :

fix bottomwall all wall/lj93 zlo 5.0 2 2 5
variable linear equal ramp(170,105.0)
fix topwall all wall/lj93  zhi v_linear 2 2 5
run 100000

The top wall moves down to 105A in 100000 steps to compress the system in order to set the correct density.To avoid close contacts I built the structure with low density at the beginning.
Still no problem here. Bottomwall is fixed.
I start data collection part, nvt simulation at 500K for 2000000 steps.
After 243398 steps I received this error:

ERROR on proc 2: Dihedral atoms 17112 17117 17118 17032 missing on proc 2 at step 243398 (../ntopo_dihedral_al

I searched this error, it says usually bad construction of structure like close contacts of atoms etc. However I do not think thats the case here because the simulation runs more than 300 000 steps before getting this error and also minimizes without any problem. Here is my neighbor settings:

minimizing and ​running for a long time does not automatically guarantee, that close contacts won't happen. they can also happen due to bad force field parameters, e.g. missing or too weak repulsive terms or incorrect charges or missing or incorrect bond/angle/dihedral interactions.

neighbor        7.0 bin

​a large neighbor skin, is not really helpful here and slows you down a *lot*.
optimal performance is usually somewhere between a skin size of 1.0 and 2.0​
a better way to prevent "long" dihedrals to be stretched too far across the ghost area, is to increase the communication cutoff using comm_modify. however, that does not help, in the case of problematic force field settings.
neigh_modify    every 2 delay 2 check yes

​rather more helpful is probably to change this to:

neigh_modify every 1 delay 0 check yes

is the most conservative choice.

if you can reproduce the "atoms missing" error reliably, you should record a trajectory around those time steps with high resolution and visualize it to check if you see some indication of atoms "shooting" through your system. also recording thermo data and monitoring kinetic/potential energy frequently can help to identify the cause. in case of close contacts, you should see an unusual spike.


I increased the skin size to see it helps but no luck. Increased the neighborlist build frequency but still error persists.
Finally my force field parameters:

pair_style lj/cut/coul/cut 10.0 8.0
bond_style harmonic
angle_style harmonic
dihedral_style opls
pair_modify mix geometric tail yes
special_bonds lj/coul 0.0 0.0 0.5

#force field parameters

pair_coeff 1 1 0.066 3.5  # C sp3 carbon
pair_coeff 2 2 0.03 2.5   # H bonded to Carbon
pair_coeff 3 3 0.140 2.9  # Oxygen

bond_coeff 1 268.0 1.529 # C- C bond
bond_coeff 2 340.0 1.09  # C-H bond
bond_coeff 3 320.0 1.41  # C-O bond
bond_coeff 4 400.0 0.96  # O-H bond
bond_coeff 5 525.0 1.48  # O-O bond

angle_coeff 1 58.35 112.7 # C-C-C
angle_coeff 2 37.50 110.7 # C-C-H
angle_coeff 3 33.00 107.8 # H-C-H
angle_coeff 4 80.00 118   # O-C-C
angle_coeff 5 80.00 118   # DUMMY
angle_coeff 6 80.00 118   # DUMMY
angle_coeff 7 80.00 118   # O-O-C

dihedral_coeff 1 0.527 -0.186 0.900 0.0 # C- C- C - C
dihedral_coeff 2 0.0  0.00 0.150 0.0 # C- C- C - H
dihedral_coeff 3 0.0  0.00 0.150 0.0 # H- C- C - H
dihedral_coeff 4 0.0  0.00 0.820 0.0 # O- C- C - C
dihedral_coeff 5 0.0  0.00 0.468 0.0 # O- C- C - H

I suspect my non bonded parameters but the density of the system is compatible with experiment which is a nice sign for correct vdW parameters.
Is there anything that I am missing?
Thanks everyone.
Dundar Yilmaz

Check out the vibrant tech community on one of the world's most
engaging tech sites,!
lammps-users mailing list

Dr. Axel Kohlmeyer  akohlmey@...92......
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.