Fix shake

Dear all, I am using lammps to simulate a carbon chain with more than 100 atoms in vacuum, using CVFF force field. In NVT ensemble, the distance between C-C bond is often too large, resulting in simulation failure. However, the fix shake command cannot be applied to C-C bond because there are more than 4 atoms. I don’t know what to do. Is it necessary to increase the force constant of the C-C bond? How many times does it generally need to be increased?

You don’t solve problems with a simulation by suppressing the symptoms that are an indication of the problem. You have to first understand where the issues are originating and only then can you think about addressing the situation.

What exactly are the errors you are facing and when do they happen?
Why do you conclude that the bonds are “too long”?

I’m sorry I didn’t make it clear. In the NVT ensemble, after about 188 ps, bond atoms missing result in simulation failure, I want to do something to keep the simulation going, details attached
nohup.out (283.9 KB)

  • It makes no sense to run a simulation with less than 4000 atoms on 96 CPUs. You are wasting CPUs and you have serious load balancing issues due to that.
  • Second, have a look at the timing summary. You are spending far too much time on Kspace.
  • If a simulation can run for so many steps and then crashes, something unusual happens. That needs to be captured in a dump file and visualized and carefully checked.
  • Does your system contain hydrogen atoms? There is crucial information missing to say for certain, but it seems likely. In that case, a timestep of 1fs is too large. That could cause a situation where atoms accelerate rapidly and cause a bond atom missing error.
  • Please note that you have shake issues before the bond atom missing error. It also can be that the shake issue is causing the bond issue. However, the shake issue can also be caused by a too large timestep. Unfortunately, there also is no information about the fix shake setting and how well it can maintain the shake constraints
  • You say you are using nvt integator. Does the system conserve energy well, when run with nve after the first chunk of simulation with fix nvt? Again, this can be an indication of a too large timestep or other inaccuracies of your simulation settings.

Additional question: is this using the data file and input you posted here?

Data file and input file are attached (680.5 KB) (2.0 KB)

How did you create the data file and more specifically how did you assign the force field parameters?

There are several atoms where the Lennard-Jones parameters are zero. That includes mobile atoms with a non-zero charge. That is very bad and could be causing the problems you are seeing. Once these atoms come closer to an atom with an opposite charge, they can collapse and cause extreme acceleration or a crash due to division by zero if they - by chance - can come exactly on top of each other. Usually, the LJ parameters prohibit that because its repulsion is stronger than the Coulomb attraction.

The atoms with L-J parameter 0 is the Hydroxyl hydrogen atom, this is set by the force field. The force field is assigned by Materials Studio software with CVFF

Does Materials Studio directly output LAMMPS data files? That would surprise me.
So you likely used msi2lmp, right? What command line exactly did you use for that?
Did you get any warnings when running it?

Yes, I did as you said. I used Materials Studio to build two models, and output data file by mis2lmp with CLAYFF and CVFF, respectively. Then use the “read_data” command in LAMMPS to combine the two models. Details attached. (1.4 KB) (334.1 KB) (54.1 KB)

I need to know the exact msi2lmp command that you used in both cases.

that is “./msi2lmp.exe Otay-MMT_12x7 -frc ./clayff.frc -i -s 0.163558254 0.344586616 -0.0006982185” and “./msi2lmp.exe PAAS_Iso_0.4 -frc ./cvff.frc -i -s -0.59066686 7.049548922 0.111241514”

That is what I suspected. What happens, if you leave out the “-i” flag?

If leave out the “-i” flag, the “CLAYFF” unable to output normally, while “CVFF” output normally. Details attached.
CLAYFF_Out.dat (530 Bytes)
CVFF_Out.dat (548 Bytes)

Well, first of all, mixing different force fields is rarely a good idea.
But the fact that you are missing parameters for one part of your system is a problem that must not be ignored. You have an incomplete description.

But I suspect the issue of not having any LJ repulsion on the hydroxyl hydrogens may cause the problems, because it looks like the shake and missing bond atom issue happens when the molecule approaches the surface. Basically, you have an imbalance/incompatibility there, that cannot easily be fixed with any LAMMPS commands. You will likely need to get a better, more complete and consistent set of potential parameters.