Hi everyone,
I ran into some trouble with my latest simulations and since I’m quite new to the LAMMPS code, it would be great to get some help and maybe even a suggestions of improvement for my simulations.
Just to get a basic idea what the simulation aims for: My goal is to use the TraPPE UA force field of ethane (also the one of propane and the CO2 model) for adsorption processes on different metal surfaces. The simulations with CO2 worked in most cases very well, but ethane and propane caused almost always „nan“ errors.
you have very high initial potential energy. that means, your initial
configuration is not optimal and there may be overlaps (possibly
through periodic boundaries), this will result in unstable time
integration and hence the NaNs.
this can be remedied by either running a minimization first (crank up
your bond/angle/etc force constants to mimic rigid bonds, if needed)
or by addressing the geometry issue directly when creating it and
converting it to a data file.
That’s why I tried to just simulate the ethane and propane model without any other compounds at different densities (from densities in the gas phase to densities within the 2-phase region). The results of these tests showed that only simulations with low densities or low box dimensions worked. So, for example a simulation with a density of 60 kg/m^3 (gas phase, close to phase boundary) worked with a cubic box with an edge length of 50 Å but not with 120 Å. As a different example, a simulation with a density of 150 kg/m^3 (within 2-phase region) ran in a 20 Å box but not in a 50 Å box. I attached the input and a shortened data file of the latter one to this message.
please note, that these are not LAMMPS specific issues, but rather
practical MD simulation knowledge. this is best learned from somebody
local, that can tutor you. there are multiple research groups at your
institution that have this knowledge.
A few more things, I’m a little wondering about, are:
The Ewald summation is supposed to be used for the model, but I can’t make it work the way I tried to implement it so far
using ewald summation is straightforward (you probably want to use
pppm for better performance, but that is besides the point). there are
plenty of example inputs and information in the manual.
however, your input data file has no (partial) charges assigned to any
atom. if there are no charges, there are no coulomb interactions and
there is not point in using a kspace solver.
I never got the minimize command to run, is it maybe connected to the rigid model I use? I thought it might speed up the equilibration of my simulations
you cannot do rigid body position updates during minimization. that is
a "hard" problem. but you can do a minimization *before* you define
fix rigid and then just increase the bond/angle/dihedral force
constants (temporarily) to have the same effect.
I wasn’t able to set the temperature damping parameter in all my simulations to lower values than 20 fs. Is it because of the neigh_modify command? I was actually used to work with temperature damping parameters of around 1-10 fs.
the tdamp factor should be *much* larger than your time step. the
purpose of a nose-hoover thermostat is to have a *weak* coupling
between the thermostat DOFs and the simulation DOFs. so you want
something of the order of 100-1000fs for a simulation with a timestep
of 1fs. this is explained in the manual.
no, this has nothing to do with neighbor list updates. your settings
for that are quite aggressive, but seem to be ok for this specific
system. please carefully check the "dangerous builds" output property.
if this is different from zero, you have to change those settings to
be less aggressive.
axel.