I am running simulations in the dreaded LJ units. At this point, I have a system of a single molecule of type ABBB with no charges, but LJ, bond, angular and dihedral interactions are present.
The simulation is being run in a NVT ensemble after minimization, using the standard Nose-Hoover thermostat. However, for some reason, my temperature is not equilibrating. What could be the reason for this?
I have attached my sys.lipid.data, sys.lipid.settings, and in.lipid files to this message.
However, for some reason, my temperature is not equilibrating. What could be the reason for this?
This line:
fix nvt_dynamics all nvt temp ${T0} ${Tf} 5
does not follow the usual rule of thumb for initial Tdamp values. As stated in the doc:
A Nose-Hoover thermostat will not work well for arbitrary values of Tdamp. If Tdamp is too small, the temperature can fluctuate wildly; if it is too large, the temperature will take a very long time to equilibrate. A good choice for many models is a Tdamp of around 100 timesteps. Note that this is NOT the same as 100 time units for most units settings. A simple way to ensure this, is via using an immediate variable expression accessing the thermo property âdtâ, which is the length of the time step.
Yet your timestep is defined as 7.5e-4 so taking a value of 5 is not likely to work.
A couple of remarks concerning your input files:
LJ units are dimensionless. That is expressed using other physical quantities as references. For example taking your lennard-jones energy eps, you can express a bond constant k in eps reduced units by taking k/eps. If eps is 4 kcal/mol (arbitrary value) and your bond constant 40 kcal/mol, then your bond constant becomes 10 in reduced units. As a ratio this is a unitless value. I still see units written in your sys.lipid.settings constants. If you are using lj units, this makes no sense.
In LAMMPS, lj units are expressed with regards to a mass, a distance, an energy and a kb unit. As stated in the doc, time is expressed in sqrt(eps/m*sig^2) units. A comment like timestep 0.00075 # real-time units does not make much sense.
You could play around with pen and paper on dimensional analysis trying to get the how and why of units conversion (for example why is the unit of the PV product analogous to an energy and what is its conversion factor?) to get a better feel on the lj units you are working with. You can take a look in the examples/UNITS directory, lj unit simulations and try to compute expected values from other units system using the conversion factors before running the other simulations.
Using a too large thermostat time constant can delay reaching equilibrium, but is not a general problem. The thermostat will eventually reach the desired target temperature unless there are other causes that may lead to rising energy. ⌠and that is the case here, too.
There are multiple settings in your input that are rather dubious:
you are using a LJ cutoff of 12 which is extremely large for reduced units. Since a typical sigma for molecular systems would be between 3 and 4 angstrom, that would correspond to a cutoff of 40 angstrom. reducing it to a more typical value of 3 or 4 exposes other problems with your input deck, but also significantly speeds up your calculation
you are updating your neighbor lists only every 100 steps. that will likely lead to problems. lowering it exposes your problems even more
your bond and angle coefficients look more like they are in real units instead of reduced units
Now after adjusting your settings to more sane values, you get all kinds of bond/angle atom missing crashes, especially when running in parallel. This is alleviated to some degree by using a shorter thermostat time constant. But it is important to observe what is happening to your system. If you look at the temperature or total energy, you see that there are âspikesâ and those are an indication that your force field parameters are bad: atoms get too close, then forces get too large, then atoms are accelerated too much and then you âlooseâ atoms or have âdangerous buildsâ for the neighbor list.
even though your timestep is already reduced from the default, it is still far too large to get a stable numeric integration of the equations of motion. you would need to reduce it by another order of magnitude to resolve this, and you will still see spikes (albeit softer) in temperature and energy unless you reduce the thermostat timeconstant to âdrainâ the energy faster than it can be created.
bottom line. bad parameters â bad simulation (or in short GI-GO = garbage in, garbage out).
Thank you for your commens, @akohlmey, @Germain! I have redone my simulation. I admit, I am not particularly versed with LJ units, and I do not have a good understanding of what parameters seem âreasonableâ.
About the bonded parameters (bond, angle, dihedral), I was following this paper.
After changing up the LJ/Coul cutoff lengths, and updating neighborlists sooner, and reducing the damping time for the thermostat, I have gotten the simulation to run stably.
I revamped the simulation per the paper, and ran it using the NVE+Langevin thermostatting mechanism, but all in all, things have worked per your comments.
Sure. I have added a charge to atom type 2, I was playing with the temperature and is now at 1.8. In my VMD, they look simply like worms on a lattice, without any critical blowupsâŚ
Thanks. The simulation runs without crash, but a running simulation doesnât mean it is a correct one.
I still think it is bogus for three reasons:
the structure after minimization is very different from your initial structure. Bonds that had around 1.5 sigma length became extremely short (because of your bond_coeff setting).
you need far too small a timestep for a system of this kind. that means you have some very steep potential somewhere or something is too soft. with charges there is often a problem that the coulomb interaction is balanced with a barrier, but when that barrier is too soft you need to suppress fluctuations to avoid the charges overcoming those barriers and then interactions collapse.
you need too short a thermostat time constant. after equilibration you should be able to run without a thermostat, but that is not the case. so the thermostat is required to suppress normal fluctuations.
If, for example, I change your equilibrium bond length to 1.5, it is possible to run a much larger time step, even the default timestep. Also it is possible to use a much larger time thermostat time constant, but not without a thermostat and for constants like 0.1 tau the energy conservation is already pretty bad.
In summary, your force field is still bogus, just not quite as extremely so as before.
Thank you for your insight.
I have some questions about your arguments:
Isnât that the point of minimization, though? If I make a bad structure using some molecule-making software like Avogadro, minimization helps me get to an optimal point where energies are not too crazy. Going from 1.5 to 0.2 is what I expect from the minimization process. Would you say that is not a sound way of using minimization?
The way 1 fs is a âgood first guessâ for a real units MD simulation; what would you classify as a âreasonable time stepâ for an LJ simulation of this density? From the de Pablo paper, he mentions that 0.00075 is sometimes what you gotta do. I see that my simulation breaks at timestep = 0.005, but runs at timestep = 0.001.
I have actually never heard of being able to remove the thermostat at equilibrium. If I donât have a thermostat, doesnât that simply mean I am integrating Hamiltonian equation of motion => running an NVE simulation? I do not understand this point⌠Letâs take changing bond-length as an example. Why is changing bond-length affecting stability/timestep? I know bond length affects how close atoms are, but if I have turned off nonbonded interactions using special_bonds, why should increase in bond length allow for a larger time step?
This ended up being longer than expected. I am genuinely curious, and I would really appreciate your word on this!
Not automatically, but it then questions your way of generating the initial structure.
If 0.2 sigma is the intended distance, then why create them at 1.5 sigma? That is a huuuge difference and all kinds of bad things can happen when starting so far from the intended geometry.
My statement was based on the assumption that your initial geometry is representative of what you want to simulate.
This is not so easy to tell. If you have only LJ with epsilon/sigma/mass 1.0, the 0.005 tau is a reasonable start. But the addition of bonds and charges complicates things. And particularly, if those bonds are supposed to be that short relative to the diameter of the LJ interactions.
This is a common technique to check on energy conservation for (purely) molecular systems. Yes you get an NVE ensemble instead of NVT, but that distinction becomes smaller the larger the system gets.
Things get a bit different for systems where you use a langevin thermostat to represent a solvent in addition to the molecules. In those cases the time constant of the thermostat determines the viscosity of the implicit solvent and thus is an input parameter.
But for cases where you use a nose-hoover thermostat, you want to âdisruptâ the simulation as little as possible, so you usually try to have the thermostat couple to the system very weakly and ideally to the border of some very wide bands in the IR spectrum. If you couple too strongly, you can get artifacts or suppress desired fluctuations thus not get the desired statistical sampling. Since you were initially using a nose-hoover thermostat, I ruled out the implicit solvent case and thus concluded that it should be possible to run without a thermostat.
If you exclude bonds at a length of 0.2 sigma, you maximally exclude bonded interactions up to a distance of 0.6 sigma (3 bonds in a line). The next atom in the line (at 0.8 sigma) is still interacting strongly via the LJ potential. With a longer bond, those atoms are much further apart and interacting much weaker. And this is assuming the equilibrium distances. Depending on the strength of the bond potential (which doesnât seem very strong compared to molecular all-atom systems) you can also have significant fluctuations on top of that, which can bring the 1-5 pair unexpectedly close, same also when the chain is not stretched, but bending. If you are that high up on the repulsive branch of an LJ potential, you need to use a very short timestep for a stable time integration.