Hi,
I’m encountering an error that I haven’t seen for the past decade. I have a model which contains 4 water reservoirs (and two silica slabs and some molecules irrelevant to this error). Each water data file (with bonds and angles) is read and shifted so that they are far from each other like 4 separate drops of water.
The minimization energies look like this:
Energy initial, next-to-last, final =
3072175.56291282 -6.27252401814407e+16 -6.27252401814407e+16
I noticed that after outputting the energy contribution only the pairwise energies are orders of magnitude off:
After dumping the atoms with their pairwise energy I found two hydrogen atoms that have the nonphysical energy values.
These are two separate H atoms from two different water molecules (however water molecules IDs are consecutive)
They have the exact same energy values.
They are 35 Angstroms apart and are relatively in the middle of the simulation box (so that the images are also too far to generate such a large energy values)
Weird thing is that they exist in the third time the water data file in read, the other ones have regular energy values.
IMPORTANT: I plotted the distance vector VS time during minimization and noticed that the vector flips
2.118 -2.4147 34.7487
-2.117 2.4168 -34.7483
This happens without the atoms moving too much (they don’t go through the PBC).
Forcefield:
pair_style lj/charmm/coul/long 10 12
bond_style harmonic
angle_style charmm
dihedral_style charmm
special_bonds charmm
improper_style harmonic (there is no improper)
What is causing this error and how to fix it?
Please let me know if you need more details, files, etc.
@Caro This is all too vague to make any specific assessment and give suggestions.
I feel like I sound like a broken record, but here we go…
For starters, we need a minimal (i.e. small, fast to run, and as few commands as possible) input deck that can reproduce the issue. We need to know which LAMMPS version you are using and what platform you are running one with which command line.
Bottom line, the simpler it is to reproduce what you are seeing, the faster we can figure out what is going wrong or where there is a mistake.
All your other observations are helpful, but useless unless we can reproduce what you see.
Thank you for your response. I wanted to see if this was something trivial that you might have seen before, before I bothered you with files.
Here are the files attached.
Another observation: When I remove the molecules (Linoleic acid data file) it works, but I have no idea why would it make hydrogen atoms in the water molecules to have high energies?
lammps-29Aug2024
Ubuntu 22.04.5 LTS test.zip (837.9 KB)
There are multiple warnings about image flags when I run your input. That is not a good sign.
Also when I add thermo 1 before minimization and change the output to show more details about energies, I can observe that something very bad is starting to happen with the Coulomb energy around minimization step 42. This is usually an indication of a bad force field, one or more pairs of atoms with opposite sign charges can get too close due to an insufficient LJ repulsion. The most common case of that is when mixing and matching parameters from different studies and for different materials. This is something you will have to work out for yourself.
Thank you. Yes, that is, I think, exactly the problem. when I set the charges on the molecule atom types to zero, it started to work.
I guess my best option is either to find a different set of parameters, or do a quick DFT on the molecules and fit an LJ on it.
Note that you are running your water molecules with SHAKE fixing the bond lengths, but not angles. Please double check that this is intended. Also, if you have linoleic acid – O-H bonds and angles are generally tricky and often worth making rigid in organic molecules, not just in water.
This is the first time I’m working with Lineloic Acid, Are you suggesting to use SHAKE on the Acid molecules? Or fix rigid? because I could make the molecules RIGID, but it would make it less “realistic”.
You can use SHAKE to freeze OH bond lengths and X-O-H angles by specifying the relevant types.
Note that you will have to verify that your simulation runs correctly. I first started freezing hydroxyls when running a nonequilibrium graphene oxide simulation – since we were trying to simulate fairly high shear rates, it wasn’t going to be a particularly energy-conserving simulation anyway, and I visualised my simulation frame-by-frame and verified that other than those few O-H groups everything else in the simulation was behaving as intended.
It is up to you to know what your simulation “should” do – in particular, what molecular structures or dynamics correspond to the experimental observables you are trying to model. Only then will you know whether a particular shortcut is acceptable or not.
In particular, if your molecules are showing unphysical dynamics because they are interacting with other particles using a force field that has not been parameterised for that interaction then you are in uncharted waters.
I would be particularly worried about the parameters for silica. There are many parameter sets available that are only suitable for bulk silica and not for slabs or interaction with other materials.
Please also not that popular force fields like OPLS-AA, Amber, or CHARMM all have very different parameterization strategies and thus different balances between partial charges and LJ terms. Generic force fields have low accuracy and Dreiding is only parameterized on representing molecular structures, not non-bonded interactions. So there are lots of pitfalls.