[lammps-users] Problems with fix npt - could it be a bug or am I doing it wrong?

Dear LAMMPS users,

I have for quite some time been doing simulations of cracked iron crystals that are deformed (tensile tests). In most of these simulations I have been using a "fix npt" to control pressure and temperature, and a "fix deform" to deform the system in one direction. In addition, I have used a "fix temp/rescale" to keep the temperature stable. I have run many simulations at 5 K (or other very cold temperature values), and things have looked okay.

Lately, I have started doing simulations at higher temperatures (300-400 K), using the same kind of fixes, and I have observed very strange behavior (e.g. that the material is more brittle at high temperature, which is certainly not the case). This made me look more closely at how the pressure was varying in the undeformed directions. My fixes (after some initial relaxation) look like this:

fix 1 all npt temp 300.0 300.0 100.0 x 0.0 0.0 1000.0 z 0.0 0.0 1000.0 couple none
fix 2 all temp/rescale 100 300.0 300.0 0.05 1.0
fix 3 all deform 10 y erate 0.0005

So I should have 0.0 pressure in the x and z directions (at least after a while, or on average). However, the pressures in these directions are far from zero, and grow as the sample is deformed. When the crack finally advances it gets really crazy. This increasing pressure happens also at low temperatures, but the effect is smaller here.

After noticing this I started to play around with different thing in order to debug the problem. I made a sample without a crack, and tried also the "nreset" keyword in the "fix npt" command in order to see wether it helped. It did not. I also tried all of this without the "fix temp/rescale", and it did not help.

As this seemed to be linked to temperature (the effect was larger at higher temperatures), I ran simulations of small perfect crystals, using the "fix npt" command and no temperature rescaling. What I notice is that the temperature is not at all stable, no matter what combination of nreset, drag og other keywords I try. An example of the fix I have used (on a perfect crystal with all periodic boundaries) is this:

fix 1 all npt temp 400 400 1.0 drag 2.0 iso 0.0 0.0 10.0 nreset 1

What I see in all cases is that the system temperature immediately drops from 400 to about 200, something whichs happens in all simulations I have tried without using the "fix temp/rescale" command. (More generally, it drops from T to about T/2 in all cases.) After a looong time (e.g. about 300000 timesteps in some cases) the system approaches the desired temperature with some combinations of parameters. Other combinations make the waiting time shorter, but I am not able to avoid the initial sudden drop from T to T/2. Also, the temperature fluctuates quite a bit (+- 50 K) around the target temperature once this has been reached.

My question is: Could this be some kind of bug, or is it just a whole study in itself to find the right parameters to use in a "fix npt" command?

Sorry for the (too) long mail, but I really want to find out of this. It is either something fundamentally wrong with LAMMPS or with how I have been using it all this time. Any comments and suggestions are welcome!

Sincerely,
Christer H. Ersland

PS: I use metal units, bcc lattice and an EAM iron potential, in case anyone was wondering.

dear christer,

Dear LAMMPS users,

I have for quite some time been doing simulations of cracked iron crystals that are deformed (tensile tests). In most of these simulations I have been using a "fix npt" to control pressure and temperature, and a "fix deform" to deform the system in one direction. In addition, I have used a "fix temp/rescale" to keep the temperature stable. I have run many simulations at 5 K (or other very cold temperature values), and things have looked okay.

that low a temperature is easy. nothing really moves. :wink:

Lately, I have started doing simulations at higher temperatures (300-400 K), using the same kind of fixes, and I have observed very strange behavior (e.g. that the material is more brittle at high temperature, which is certainly not the case). This made me look more closely at how the pressure was varying in the undeformed directions. My fixes (after some initial relaxation) look like this:

fix 1 all npt temp 300.0 300.0 100.0 x 0.0 0.0 1000.0 z 0.0 0.0 1000.0 couple none
fix 2 all temp/rescale 100 300.0 300.0 0.05 1.0
fix 3 all deform 10 y erate 0.0005

So I should have 0.0 pressure in the x and z directions (at least after a while, or on average). However, the pressures in these directions are far from zero, and grow as the sample is deformed. When the crack finally advances it gets really crazy. This increasing pressure happens also at low temperatures, but the effect is smaller here.

After noticing this I started to play around with different thing in order to debug the problem. I made a sample without a crack, and tried also the "nreset" keyword in the "fix npt" command in order to see wether it helped. It did not. I also tried all of this without the "fix temp/rescale", and it did not help.

what i am missing here is that you describe tests using a shorter time step.
did you try that?

As this seemed to be linked to temperature (the effect was larger at higher temperatures), I ran simulations of small perfect crystals, using the "fix npt" command and no temperature rescaling. What I notice is that the temperature is not at all stable, no matter what combination of nreset, drag og other keywords I try. An example of the fix I have used (on a perfect crystal with all periodic boundaries) is this:

fix 1 all npt temp 400 400 1.0 drag 2.0 iso 0.0 0.0 10.0 nreset 1

What I see in all cases is that the system temperature immediately drops from 400 to about 200, something whichs happens in all simulations I have tried without using the "fix temp/rescale" command. (More generally, it drops from T to about T/2 in all cases.) After a looong time (e.g. about 300000

this would be perfectly normal, if you start from a perfect crystal
and are not running an equilibration first. at equilibrium you have
to have about half of the initial kinetic energy converted into potential
energy.

timesteps in some cases) the system approaches the desired temperature with some combinations of parameters. Other combinations make the waiting time shorter, but I am not able to avoid the initial

the npt thermostat is not so good at transferring kinetic energy
and particularly guaranteeing equipartitioning. i would recommend
to first run for a while with fix nve plus fix langevin with a reasonably
short time constant and only then switch to nvt and then ntp.

sudden drop from T to T/2. Also, the temperature fluctuates quite a bit (+- 50 K) around the target temperature once this has been reached.

you may have some resonances with the lattice due to the initial
"push" you give your system. to some degree those fluctuations
are to be expected and depend on the system size. fix langevin
at the beginning should help to dissipate very large initial fluctuations.

My question is: Could this be some kind of bug, or is it just a whole study in itself to find the right parameters to use in a "fix npt" command?

Sorry for the (too) long mail, but I really want to find out of this. It is either something fundamentally wrong with LAMMPS or with how I have been using it all this time. Any comments and suggestions are welcome!

one recommendation on top of all.
please don't run fix npt or fix nvt at the same time as fix temp/rescale
this is a bad combination and if you feel you need it, then something
else is wrong. you are just trying to use a sledgehammer fix to hide
the real problem.

cheers,
    axel.