I am running a simple system of one methanol molecule in 500 water molecules (SPC/E)
My simulation keeps crashing.
Per MPI rank memory allocation (min/avg/max) = 14.09 | 14.09 | 14.09 Mbytes
Step Temp PotEng KinEng Press Density
1 300 19902.939 898.71423 260270.24 0.96067272
ERROR on proc 0: Out of range atoms - cannot compute PPPM (src/src/KSPACE/pppm.cpp:1839)
I have tried the exact same input structure with GROMACS and the simulation works perfectly fine. I can see that energy minimization does nothing to the structure in LAMMPS but the post EM structure I get in GROMACS is starkingly different than the initial structure.
I don’t, obviously, expect that post-EM structure match in LAMMPS and GROAMCS, but I expected LAMMPS to change the structure to some extent.
I have attatched the files I am using to simulate. Please let me know if you need any further information.
Kindly help.
Thanks and regards,
Ved
GROAMCS .itp file:
[ defaults ]
; nbfunc→ comb-rule gen-pairs fudgeLJ fudgeQQ
1 3 yes 0.5 0.5
[ atomtypes ]
opls_116 OW 8 15.999 -0.8476 A 3.16557e-01 6.50194e-01
opls_117 HW 1 1.008 0.4238 A 0.00000e+00 0.00000e+00
MTH0 C 12.0110 0.000 A 0.35000 0.27614
MTH1 O 15.9990 0.000 A 0.31200 0.71128
MTH2 HA 1.0080 0.000 A 0.25000 0.12552
MTH3 HB 1.0080 0.000 A 0.25000 0.12552
MTH4 HC 1.0080 0.000 A 0.25000 0.12552
MTH5 H 1.0080 0.000 A 0.00000 0.00000
That is highly unlikely, since your initial geometry is highly symmetric. For the minimization to have some effect you need to break the symmetry, e.g. by adding a random displacement to atom positions.
For instance with:
displace_atoms all random 0.5 0.5 0.5 21344
There are multiple other problematic choices in your input:
fix myshake spc shake 0.0001 10 100 b 4 a 4
This applies SHAKE only to the hydrogen bond/angle in water, but you also have hydrogen atoms in your methanol molecule, so you need to shake their bonds as well or reduce the timestep. I you apply fix shake to group all and also shake bond types 1 and 3
neighbor 0.0 bin
neigh_modify every 2 delay 10 check yes
A neighbor list skin of zero is a very bad choice. This means that the neighborlists would become invalid after every MD step and would need to be rebuilt. This is in stark contrast with your neighbor list build setting. A delay of 10 is completely bogus for a timestep of 2 fs. I would suggest to stay with the default settings of 2.0 skin, delay 0, every 1, check yes.
A timestep of 2.0 can work for your system when it is well equilibrated, but not so well when you have a high potential energy geometry. So reducing it to 1.0fs during the initial equilibration seems like a good idea.
Since you have such a high initial potential energy (and pressure), it is a bad idea to start with fix npt. This will needlessly expand the box and shrinking it back to the “normal” size can take a very, very long time. Better to run with fix nvt until most of the excess potential energy has been removed by the thermostat and the temperature and pressure gone down.
Thank you so much for such a detailed response. It’s very helpful and I highly appreciate it. I have certain questions about your suggestions:
Regarding EM, is the a certain feature of LAMMPS EM? GROMACS is able to minimize the energy when given the same structure.
Applying SHAKE to OPLS-AA methanol model is going to change the FF for methanol, right? The I would be running a different model than OPLS-AA. Am I right?
I had tried exactly this but the simulation still kept crashing.
If I now, incorporate your suggestions into my input script:
# -- Initialise system --
units real
atom_style full
# -- read data --
read_data system.data
displace_atoms all random 0.5 0.5 0.5 21344
# -- set charges --
# water
set type 1 charge -0.8476
set type 2 charge 0.4238
# methanol
set type 3 charge -0.683
set type 4 charge 0.418
set type 5 charge 0.04
set type 6 charge 0.145
# -- non bonded interactions --
include "lj.settings"
#include "ffnonbonded.settings"
# -- kspace --
kspace_style pppm 0.00001
# -- bonded interaction
include "ffbonded.settings"
# -- make groups --
group spc type 1:2
group ow type 1
group met type 3:6
# -- write ff data --
write_data ff.data pair ij
# -- water shake --
fix myshake spc shake 0.0001 10 100 b 4 a 4
# -- neighbour list --
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
# -- simulation protocol --
# EM
timestep 2
minimize 0.0 0.0 1000 10000
write_data postem.data
# equilibriate pressure
velocity all create 300.0 21423623
fix relax1 all npt temp 300 300 100.0 iso 0.986923 0.986923 2000.0
thermo_style custom step temp pe ke press density
thermo 10
run 5000000
dump myDump1 all custom 100 npt_eq.lammpstrj id mol type x y z
unfix relax1
# prod
fix relax2 all npt temp 300 300 200.0 iso 0.986923 0.986923 2000.0
thermo_style custom step temp pe ke press density
thermo 10
run 5000000
dump myDump2 all custom 50 prod.lammpstrj id mol type x y z
dump myDump3 ow h5md 100 prod_ow.h5md position box yes
unfix relax2
write_data prod.data nocoeff
The temperature, pressure and density are all looking good!
What is EM? electron microscopy?? It is bad to use ambiguous acronyms in a public forum. Always keep in mind that it is you that wants good advice and that requires that you communicate well and with as little ambiguity as possible.
Minimizer algorithms can be very unpredictable. There are a whole bunch of parameters that can be tuned. Given the highly symmetric nature of your input geometry, I would consider the behavior of LAMMPS more along the expected behavior. I don’t know what is different in Gromacs.
That said, starting from a highly symmetric geometry is generally a bad idea since it can take a lot of simulation time to remove the “memory” of that highly ordered structure. Since the implementation of the “overlap” keyword in create_atoms random, I prefer that in combination with a molecule file to create water boxes (. An even better choice would be to start from a pre-equilibrated water box and then just delete the overlapping waters.
Not entirely. Running OPLS-AA with a timestep of 2fs and without fix shake on bonds involving hydrogens is definitely wrong. Using shake on hydrogens doesn’t really change the force field itself, it just removes some degrees of freedom from your simulation, so you can see it as some simple form of coarse graining (as an intermediate step to OPLS-UA).
Please note that using SPC/E water is not consistent with OPLS-AA, either, and probably a much bigger deviation than using fix shake. This should be TIP3P or TIP4P water.
This is really not a LAMMPS issue but a science question, so it would be better to discuss with your adviser/tutor/supervisor/colleagues.
It is a mistake to assume that only one change can “fix” a problem. Unfortunately, humans seem to be hooked into making this kind of assumption (not only in science).
Thank you so very much for the detailed answer. I read up on some things and I now understand what you mean. I fully agree with your suggestions. I incorporate all of them and I get an RDF that matches with my target RDF.