I prepare the Tersoff potential (TP) for GaP/GaN/GaAs crystal based on zinc-blende structure. The initial TP values are given by a fitting procedure based on results of DFT calculations . MD simulations based on it give rather far from expected values of lattice parameter, cohesive energy and a few others.
So I developed a protocol for generating thousands of samples from which one can select the most promising to further evaluations. Samples are differ only by TP parameters. There are randomly changed 6 parameters, i.e.: A, l1, B, l2, R, D. Every sample uses on input well relaxed initial model.
Simulation of a single sample is based on NPT ensemble which controls 1000 steps of relaxation for T= 300K and next the raising of temperature from 300K to 1300K and constant 1Ba pressure. In the most cases one run (11 000 steps) takes about 5 s. However and quite often the simulation time is greater than 60s and what worse, I’ve also noticed that during relaxation temperature exceeds very high values (above 2500K). The problem is that NPT doesn’t adjust the temperature, doesn’t reduce to expected 300K value even for very long relaxation time (I’ve tested one case where relaxation took 100 000 steps without reducing the temperature). Temperature is “frozen”. The same phenomena is present in NVT ensemble.
Also I created in the case of NVT the thermostats, one for each type of atom. And similarly to above temperature is frozen but for low values ( about 170-180 K); NVTs are not able to increase value of temperature to 300K.
I 've noticed that problems arise only when R, D parameters are the subjects of adjusting. The R/D params mainly regulate the range of potential. The LAMMPS docs don’t give details on implemented NPT/NVT algorithms and how they are altered by TP (or other) potentials. It bothers me for a last few days and stopped potential developing/testing.
So question is why ensembles don’t work properly for some R/D params ?
First off, please don’t use the term “ensemble” in this context. When using fix npt or fix nvt you are not automatically simulating an NPT or NVT ensemble. What fix nvt does is a) time integration and b) coupling of the system to a Nose-Hoover thermostat. For practicality reasons both steps are combined to a single fix. To get an NVT ensemble, multiple other conditions are needed. Details are in any decent text book on statistical mechanics. Similarly, the fix npt command adds a Nose-Hoover barostat to the time integration and thermostat.
This is not so much a question of the “ensembles” but rather a question of reasonable behavior of the potential and the properties of numerical time integration. If you have a system in equilibrium, there should be not much need for a thermostat, but rather the system should have a stable time integration for suitable energy conservation and fluctuations of the temperature around the equilibrium value with a magnitude correlated with the system size (the smaller the system, the larger the fluctuations). However, if your potential parameters lead to very large forces (e.g. due to lack of repulsion), then you can easily run into problems with the numerical time integration. The larger the forces (and the lighter the atoms), the shorter a timestep is needed or else the simulation will diverge, i.e. the temperature will increase in an uncontrolled fashion. Same as with the topic of ensembles, this is fundamental text book knowledge of MD simulations.
As for why the thermostat does not guarantee the desired temperature. Well, the Nose-Hoover algorithm is meant to be used for a system in equilibrium and designed to mimic the fluctuations that are required to replicate a larger system so that sampling the NVT statistical mechanical ensemble is possible. There are other thermostats that are more effective and removing excess kinetic energy, but that would be just hiding the problem. MD simulations follow the “garbage in, garbage out” (GI-GO) principle and thus when you have bad simulation settings (timestep, cutoff, potential parameters) then you cannot “correct” this with a thermostat.
As for variations in the simulation time. This can happen when your system collapses due to lack of proper repulsion. In that scenario the number of neighbors can grow very much (in cubic order with the average distance of atoms) and that can slow down your calculations substantially. Again, this is the GI-GO principle at work.
The LAMMPS documentation to fix nvt/npt does quote all relevant publications to the Nose-Hoover algorithms as they are used. However, the choice of potential is a completely different thing. Time integration and thermostat only “see” atom positions, velocities, and forces and the thermostat specific settings as well as the time step. They pay no attention to the details of how you compute forces. If the forces are bogus, the trajectory will be bogus and time integration will diverge or worse.
Thanks a lot for your answer, especially for explanation of differences between statistical ensembles and Lammps NVT/NPT algorithms. You inspired me, so yesterday I’ve done some other tests. I found that reduction timestep 2x from 0.001 to 0.0005 restored expected behaviour of barostat/thermostat and temperature T=300K. I have to think about it more, how it works ,etc… There are still some details which I need to understand.