Non-zero pressure component in contrast to fix NPT

Dear Lammps users and developers,

I have a difficulty with fix NPT to keep the xy shear component constant while all other
components are zero. After reducing the system complexity I ended up with a simple
periodic box of Cu atoms. Relaxation with all pressure components set to zero is working
correct. However, when I set one of the shear component (say Pxy) to a non-zero value,
one of the normal components (Pxx in this case) is also oscillates around a non-zero value
in contrast to the fact that it should oscillate around zero by fix NPT.

I already checked this with Cu anisotropy and there should not be a coupling between xy
and xx components. Besides the system is fully periodic with triclinic box, so any stress
state should be possible. Changing the potential has also no effect on the problem.

I tried to solve this by changing Tdamp, Pdamp, Tchain and Pchain values of fix NPT,
but no success. At the end, I found that the problem can be solved by using the same
fix NPT twice separated by a short run. First run can be even 100 steps, then the second
run will work fine which is strange. I do not understand why this is happening. Issuing
the same fix again, shouldn’t change anything. Is it a bug or some kind of stress locking
is happening in NPT for solids? In the first run, Pxx is not oscillating much and it is stuck
at around -255. Using unfix and different name for two fixed has the same result.

Here are the input file and parts of the log file:

##------ INPUT File--------------

units metal
boundary p p p
neigh_modify every 1 delay 2

lattice fcc 3.615

region whole block 0 6 0 6 0 6
create_box 1 whole
change_box all triclinic
create_atoms 1 box

pair_style eam
pair_coeff * * Cu_u6.eam

thermo 10
thermo_style custom step temp pxx pyy pzz pxy pxz pyz press vol

variable Tdmp equal 0.1
variable Pdmp equal 1

velocity all create 4 7654321

fix nsmb all npt temp 2 2 {Tdmp} x 0 0 {Pdmp} y 0 0 {Pdmp} z 0 0 {Pdmp} xy 10000 10000 {Pdmp} yz 0 0 {Pdmp} xz 0 0 ${Pdmp}
run 5000

fix nsmb all npt temp 2 2 {Tdmp} x 0 0 {Pdmp} y 0 0 {Pdmp} z 0 0 {Pdmp} xy 10000 10000 {Pdmp} yz 0 0 {Pdmp} xz 0 0 ${Pdmp}
run 5000

##------ Log File --------------

LAMMPS (12 Jan 2013)
Lattice spacing in x,y,z = 3.615 3.615 3.615
Created orthogonal box = (0 0 0) to (21.69 21.69 21.69)
2 by 2 by 2 MPI processor grid
Changing box …
triclinic box = (0 0 0) to (21.69 21.69 21.69) with tilt (0 0 0)
Created 864 atoms
Setting up run …
Memory usage per processor = 2.25883 Mbytes
Step Temp Pxx Pyy Pzz Pxy Pxz Pyz Press Volume

500 2.2111239 -252.98756 13.397256 -11.230803 10090.326 -0.33264225 -0.23636272 -83.607035 10207.053

5000 2.0185036 -255.04601 0.39089872 -0.65204267 10006.385 0.14726901 -0.22797605 -85.102385 10207.135
Loop time of 6.99789 on 8 procs for 5000 steps with 864 atoms

Setting up run …
Memory usage per processor = 2.25883 Mbytes
Step Temp Pxx Pyy Pzz Pxy Pxz Pyz Press Volume
5000 2.0185036 -255.04601 0.39089872 -0.65204267 10006.385 0.14726901 -0.22797605 -85.102385 10207.135
10000 2.0333702 -70.565601 1.1895423 -5.4735159 9999.8132 -0.068677885 -0.0026304791 -24.949858 10206.642

Thanks and regards,
Jaber Rezaei Mianroodi

Jaber,

Since you are assigning a relatively large pressure to the xy tilt dimension instantaneously via the Pstart argument, the simulation box dimensions change significantly (and almost instantaneously) so that you should reset your reference state.

From the fix_nh doc page, “The keyword nreset controls how often the reference dimensions used to define the strain energy are reset. If this keyword is not used, or is given a value of zero, then the reference dimensions are set to those of the initial simulation domain and are never changed. If the simulation domain changes significantly during the simulation, then the final average pressure tensor will differ significantly from the specified values of the external stress tensor. A value of nstep means that every nstep timesteps, the reference dimensions are set to those of the current simulation domain.”

When the second npt runs, the reference state is reset to the tilted box, so you obtain a better match between the final stress tensor and the specified external stress tensor. A nreset value of 2000 for the first npt run (where significant dimension change occurs) should work fine. Also, 5000 steps (5 ps) is far from enough to obtain good statistics, especially the Pdamp is in the order of 1 ps. You should run this much longer.

Cheers,
Ray