Problem with bond or angle atom missing

Hi there,

I am simulating a polymer system with doped salts under electric field. The force field used was PCFF. I ran it with 0.5 fs timestep. 1000,0000 steps for equilibration and 4000,0000 steps for production run under electric field. The problem with bond/angle atom missing has bothered me for a long time. I know people discussed these before many times and developers give a lot of suggestions and i try:

  1. geometry. I generated several difference configurations with emc, (some with 40% smaller density than target to make sure the atoms were not too close to each other). I minimize twice with sd and cg. And it still returns error sometimes at 500,0000 steps of equilibration or 500,0000 steps during production run.
  2. cutoff. I try with a larger neighbor bin 3.0 and a larger cutoff 15.0, it delays the error but it still occurs.
  3. system size dependent. When i try with only ten chains, the error sometimes popped out. But most of the time running with small number of cores, it is successful. But when I switched to 30 chains, it just keep bugging with atoms missing.
  4. I also checked energy and density. Energy was very low and converged. Density was relaxed to the desired one.
    I know it might be caused by some bond overstretching and the proc cannot find it. It might be potential problem but i generated the parameters completely based on pcff with emc and i am not sure how to modify the params. Can people give some suggestions regarding my situation? I would really appreciate it. I attached my input script below.
# LAMMPS atomistic input script

echo		screen
units		real
atom_style	full
boundary p p p

# Variable definitions

variable	project		index	"polymer"	# project name
variable	source		index	.	# data directory
variable	params		index	.	# parameter directory
variable	temperature	index	303		# system temperature
variable	tdamp		index	50		# temperature damping
variable    pdamp           index   500             # pressure damping
variable	dielectric	index	1.0		# medium dielectric
variable	kappa		index	4		# electrostatics kappa
variable	cutoff		index	12.0		# standard cutoff
variable	charge_cutoff	index	12.0		# charge cutoff
variable	precision	index	0.00001		# kspace precision
variable	lseed		index	723853		# langevin seed
variable	vseed		index	486234		# velocity seed
variable	trun		index	40000000	# run time
variable	frestart	index	0		# 0: equil, 1: restart
variable	dtrestart	index	10000		# delta restart time
variable	dtdump		index	10000		# delta dump time
variable	dtthermo	index	10000		# delta thermo time
variable	timestep	index	0.5		# integration time step
variable	restart		index	${params}/${project}.restart

if "${frestart} != 0" then &
"variable	data		index	${restart}" &
else &
"variable	data		index	${params}/${project}.data" &

# Interaction potential definition

pair_style lj/class2/coul/cut ${cutoff} ${charge_cutoff}
bond_style	class2
angle_style	class2
dihedral_style	class2
improper_style	class2
pair_modify	mix sixthpower tail yes
special_bonds	lj/coul 0 0 1

if "${frestart} != 0" then "read_restart ${data}" else "read_data ${data}"
include		${params}/${project}.params

neighbor 2.0 bin
neigh_modify  delay 0 every 1 check yes

# Integration conditions (check)

timestep	${timestep}
dielectric	${dielectric}
fix		mom all momentum 100 linear 1 1 1 angular

if "${frestart} != 0" then "jump SELF simulate"

group liion type 1
group otfion molecule 31:60
group backbone molecule 1:30
group ion union liion otfion

# Minimization
reset_timestep 0 
thermo 10 
thermo_style custom step temp pe etotal lx ly lz vol density press pxx pyy pzz

fix min all box/relax aniso 0.0 vmax 0.001
min_style sd 1e-5
minimize 1e-25 1e-25 5000 10000 

min_style cg 1e-5
minimize 1e-25 1e-25 5000 10000
unfix min

# specify initial velocity of atoms
velocity	all create ${temperature} ${vseed} &
		dist gaussian rot yes mom yes sum yes

reset_timestep 0
thermo		${dtthermo}
dump           1 all custom 5000 equilibration.lammpstrj id type x y z

#room temp and pressure
fix npt1 all npt temp ${temperature} ${temperature} ${tdamp} aniso 1.0 1.0 ${pdamp}  
run 1000000
unfix npt1

# constant room temp
fix nvt_eq all nvt temp ${temperature} ${temperature} ${tdamp}
run 4000000
unfix nvt_eq

#room temp and pressure
fix npt2 all npt temp ${temperature} ${temperature} ${tdamp} aniso 1.0 1.0 ${pdamp}  
run 5000000
unfix npt2

undump 1
write_restart	${project}.restart

# Simulation

label		simulate

# Integrator
reset_timestep 0

pair_style lj/class2/coul/long ${cutoff} ${charge_cutoff}
bond_style	class2
angle_style	class2
dihedral_style	class2
improper_style	class2
pair_modify	mix sixthpower tail yes
special_bonds	lj/coul 0 0 1

include		${params}/${project}.params
kspace_style pppm/cg ${precision}

fix             nvt3 all nvt temp ${temperature} ${temperature} ${tdamp}
fix            3 all efield 5e-8 0.0 0.0
fix_modify 3 energy yes
fix_modify 3 virial yes

run ${trun}
1 Like

You could try visualizing the simulation, and look at what the simulation is doing when it crashes, e.g. use OVITO or VMD with a dump file, or even the in-situ visualization in LAMMPS, see dump image command — LAMMPS documentation.

It looks normal to me. I visualized the reported missing atom bonds they are all C-H bonds on the polymer backbone within the simulation box. No bond stretching is observed. I am thinking if the missing bond was caused by the too large number of processors I used. Currently I am using 2 nodes(80 cpus) for 7800 atom system. Is that one of the possible reason?

Extremely unlikely in this case. It may be inefficient to run the system with too many MPI ranks, but the communication cutoff (i.e. how far ghost atoms from a neighboring sub-domain are copied) will make certain that atoms are available. To lose atoms in a bond with a 12 angstrom cutoff with your neighbor list update settings, you would need a massive acceleration. This is in practice most commonly observed, if you have incomplete or incorrect force field parameters, or if you are doing a simulation with flexible bonds (involving hydrogens) that are expected to be constrained by SHAKE instead. Some systems are sensitive to too large timesteps, so you may want try 0.25 or 0.2 instead of 0.5, but this is less likely than just missing or unsuitable parameters.

I tried fix shake before and it returns “non-numeric pressure unstable simulation”. I think you are right, I will check the parameters. Thank you for your time!

That error hints at another possible source of a problem: your simulation cell may be too small and atoms would overlap due to periodic boundaries.

You can check for that with delete_atoms overlap 0.2 all all. If there are close contacts, atoms will be deleted. If not, the box dimensions are OK.

Thank you for the suggestion. One question I have regarding trying with the delete_atoms command is that I set special_bond lj/coul 0 0 1, which turns off the pairwise interaction for bonded pairs. According to LAMMPS documentation, “If the special_bonds command is used with a setting of 0, then a pair of bonded atoms (1-2, 1-3, or 1-4) will not appear in the neighbor list, and thus will not be considered for deletion by the overlap styles. You probably don’t want to be deleting one atom in a bonded pair anyway”.
If I just changed the lj/coul weighting coeff from 0 to 1, it will add the pairwise interaction to the bonded pair, which might make things more complex. Will instead trying generating initial configuration with less density or just simply use a larger simulation box size also indirectly check for atom overlapping?

I feel i just misunderstood the description from the documentation. It means bonded atoms in the same pair will not appear in the neighbor list of each other. I will try with delete_atoms command.

For close contacts though PBC, you don’t need to worry about the atoms that are bonded. Those have the proper distance and the bond potential will push them apart if needed anyway without creating a problem. The issues stem from the “outside” atom (say a hydrogen of a C-H bond) is in close contact with an atom from a different molecule. That can create huge forces which could result in non-numeric pressure or atoms being accelerated so much that they move more than 15 angstrom in a single time step, which will than trigger the bond atom missing error.

A different density won’t help much, but increasing the box by a small margin can do the trick. Please note that the deleting of atoms is a diagnostic, not a fix.

yeah that really makes sense! I actually thought about doing that before but i am confused about if I should increase the box size until it includes all the atoms (so that no wrapping would be carried out due to pbc at the very first), and it will make the box very large (-10 90 A lo hi compared to 0 46A lo hi). Also I am worried about if it is needed to recenter the center of mass the molecules inside the box with the new dimensions.

No, you don’t want to make a big change. That can make matters worse. Try small increments instead.

Sure i will try that. These suggestions are really helpful. Thanks so much for your time!

An update with my progress with bond atom missing: I still cannot solve the problem, even with small margin added to the simulation box. And the thing is I cannot increase the box size too much (by about 1.2 A). If I added a margin larger than 1.2/2=0.6A, the potential energy will become very high and it will return bond missing or non-numeric pressure error at the very start of the run. Increasing the box size by a tiny amount did delay the error but not eliminate it. Also tried different initial configurations with different densities, all reported the bond missing error at some point. Checked parameter file line by line compared with pcff frc file. Not sure if there is anything else I can try.

That kind of check will not tell you whether the parameters and/or force field you selected are unsuitable.

Please provide a small test input deck including input, data, and log file.

Sorry to bother you again. I tried systems with smaller size to preproof the validity of pcff under my situation( both ten chains and twenty chains works well and the property output looks reasonable to me) but when it increases to thirty chains, it keeps bugging with the atom missing error. I will upload an input deck now. (3.3 MB) (3.3 KB)
polymer.params (29.3 KB)

log.lammps (160.2 KB)
I cannot find the log file of the uploaded script since i just rerun the simulation and things were rewritten and this log file is the output of the same script as uploaded but with fix shake. I feel it is the same reason that cause the error but with fix shake the error is reported as"non-numeric pressure/atom coord"; without fix shake it is reported as “bond atom missing”.

Those two error messages have different causes.

I don’t have an answer for what is causing the problems, but the primary reason for the non-numeric pressure with fix shake is that you should not be running fix npt at that point. You have too large pressure fluctuations and that cause instabilities in the Nose-Hoove barostat which in turn will cause the error.

Since your system density is pretty close to the equilibrium density, using fix box/relax and running with fix npt early are both counterproductive. Both take your system away from the equilibrium point rather than bringing it close. You input file is probably a template and in many cases, the sequence of simulations is ok, but in yours it makes matters worse.

So here is what I found to be more stable:

  • run only cg style minimization and without fix box relax.
  • add fix shake on the type 3 bond and increase the timestep to 1fs (instead of 0.5).
  • do not start equilibration with fix npt but use fix nvt instead

But I am still getting errors of different kinds at different locations in the system.
Nothing that is easily and reliably reproducible and most of the time those seem somewhat random. There is no indication that there is a problem with the non-bonded interactions and the parameters in use look reasonable so that it is not likely anyway.

What worries me more are the many zeros with the class2 bonded interactions. I suspect that there may be some interaction where that will cause a division by zero somewhere (which would be a bug in the source code). On x86 CPUs this will usually lead to the calculation continuing, but then leading to forces, then velocities, then positions becoming NaN and that would then trigger the errors. So far, however, I have been unable to capture such an event with an instrumented binary so that I can see where it happens first.

It will take a little while extra. I would have to compile a set of differently instrumented LAMMPS executables on our local cluster (I am due to compile/install the latest stable version anyway). And then I can try a few more variants and with faster turnaround that on my tiny desktop machine at home.