I have a linear coarse-grained polymer system, with nonbonded, bond, angle, and dihedral interactions as a tabulated potential. My simulations consist of 1000 polymer chains, each chain has 200 monomers. The simulations sometimes crash because of dihedral atom missing errors in some processors. I am running the latest version of lammps on 96 CPUs. This error is more likely to happen if I have larger systems(like 10000 polymer chains). The timestep I am using is 5 femtoseconds. The simulation ran with no problem with a lower time step (1 femtoseconds), but this is the timestep for the all-atom(AA) model and the CG model I have should be able to run with a higher timestep than the original AA model. If I remove the dihedral, the simulations run with a timestep of 12 with no problem. Even so the dihedral energy or forces is not high at all. I also get the same error if I keep the dihedrals in my data file and have zero forces for dihedral in the tabulated potential.
I will add some simulation settings in case could help to find the problem.
Units are real.
neigh_modify every 1 delay 0
neighbor 2.0 bin cutoff = 15 A
comm_modify cutoff 20.0
T = 400 K
boundary p p p
special_bonds lj 0.0 0.0 1.0
How coarse is your mapping? When I was working with CG, I found that I had to diminish the timestep when I was working at low coarsening degrees. The three mappings I investigated with iterative Boltzmann inversion and force matching models cherished timesteps of 5fs, 2fs and 1fs and it diminished together with the degree of coarsening. I dont see why it would be a problem to work with small timesteps in CG. The meaning of time would still be absent and it would rather be the parameter for the numerical approximation of the derivatives involved in the equations of motion regardless. From the top of my head, I remember this paper: Systematic study of temperature and density variations in effective potentials for coarse-grained models of molecular liquids | The Journal of Chemical Physics | AIP Publishing in wihch people used 1 fs as timestep in their CG simulation. I suppose they also could not (like me ([2312.05192] Force Matching and Iterative Boltzmann Inversion Coarse Grained Force Fields for ZIF-8 ;D) and you) get it to “stable-ly” (?) run with large ones. I think I saw others, but I dont remember so readily.
The bead mass is 70 u, it’s not a problem to use 1fs, but I believe that we can run it using higher timesteps, for longer times also. I forgot to mention that even with higher timestep if I include the dihedrals the simulations run for a few 100 000 steps and then without temperature changes (temperature fluctuation is quite low and stays on 400K+/-1K) or energy the simulations crash. I am trying to understand what could be wrong with higher timesteps when I include the dihedrals.
I forgot to mention that I am doing systematic coarse-grading ibi and the average bond length is 3.68 A.
We need a fully functional example input deck (it doesn’t have to be your production input. In fact something smaller that still reproduces the issue would be beneficial) sIt would also help if you report what exact version of LAMMPS you are using (“latest” useless as it keeps changing and depends on what people consider “latest”) and what platform you are running on with which command line.
You should have a careful look at the log file for any warnings (or best attach one of those as well).
Do you get the same error, if you use a simple harmonic potential instead of the tables?
module purge
module load intel/env/2018
module load intel/mpi/2018
srun /home/snorouzi/programs/lammps/lammps-2Aug2023/build/lmp -in input.in
on our CentOS cluster.
I tested two simulations one with 2000 atoms (10 chain with 200 monomers) and the other one 800 chains with athe same number of atoms(160000). the output file for system with 160000 atoms is attached here out.160000, and the output file for 2000 atoms is called out.2000
I have also the input file and data files here.
The problem happnes for large systems than small ones. the system with 2000 atoms runs with time step 10 fs and the other one crashes. if I make the much bigger system the even crashes with time step 2fs. all systems run with timesteps of 1 fs. out.2000 (521.2 KB)
I think the actual error is because of the dihedrals, because if I replace the bond with a harmonic one, I get the error for a dihedral atom missing, and If I remove the dihedrals everything are file.
I have not tried the harmonic with the dihedral.
I cannot attach the data file for 160000 atoms. it seems to be large file to upload here.
even if I choose a bigger cutoff range I have the same error. for small systems no problem with timesteps but for the larger ones immediately crashes with timesteps higher than 1fs.
That is why I urged you to provide a small input example for testing that still (eventually) reproduces the crash. I suspect that the problem may be due to the tabulation of the potential, but that is just a guess.
here is my tabulated potentials, I checked the forces for dihedrals evey this seems file even if I replace the forces with 0 means no force on dihedrals but the dihedrals exist in the data file I have the same error. angles.dat (8.9 KB) bonds.dat (61.1 KB) dihedral.dat (10.4 KB) pairs.dat (58.1 KB)
Ok. This is it. This is still the huge system, you have failed to provide the additional information I have asked about and are making me ask for stuff again and again. I have lost my interest.
If you want help from people that volunteer their time to help you, you have to respect their effort by making it easy to help you. What you are doing is the opposite.
Perhaps you get lucky and someone else is willing to spend the effort and help you anyway.
I don’t have the time at the moment (neither possibly the competence) to try to help you find your answers, but if I was going to come with a hypothesis, I would tend to think that the more bonded potentials you put, the more forces restraining the movement of the superatoms you have. It is not necessarily true always, but if i am right and the dynamics is more restrained, it makes sense that the simulation crashes due to numerically instability as you add more and more potentials restraining the range of motion. So it would explain why it crashes when you add the potentials up to dihedrals but not necessarily only with up to angles. And the crashing doesn’t need to happen right away: in principle, there could be specific “phase space trajectories” that lead to collapsing, and they may or may not happen in a more readily fashion.
You can find out if the force field becomes “less allowing” for larger range of motions if adding dihedral potentials by making histograms for bonds, angles, dihedrals and also rdfs and seeing how broad the peaks become. However it is true that it doesn’t prove that this is the cause foe the crashing. Yet, I guess it’s a valid hypothesis.
You asked me to send you the smallest possible system that shows the problem and I did. Systems smaller than that do not show the problem and run fine. Thanks for your effort and sorry for taking your time.
I fixed the problem. it was the wrong spline interpolations for the inner cutoff of the pair potential. The created tables by lammps for the potential energy for the inner cutoff was very high for two atoms closer than 1 A, so they were experiencing a very strong forces. I fixed it by using larger inner cutoff and also using the linear style for tables.
