[lammps-users] fix npt xyz

LAMMPS LAMMPS Mailing List Mirror
1

Hi. I have a system in which the pressure pzz is different from pxx and pyy because the system is anisotropic. However, their sum is zero, and I am trying to do a run in which the ratios of the cell lengths along x, y, and z are held constant and the total pressure P = (pxx + pyy + pzz)/3 is held at zero. I thought that the following command

fix 1 all npt 0.4 0.4 1 xyz 0 0 10

would do the trick, but for some reason the ratios of Lx to Ly to Lz change (greatly) during the run.

Here are snippets of the log file

units lj
atom_style angle
boundary p p p
processors 4 4 1
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes
dump 1 all custom 100000 configs.holdiso tag x y z ix iy iz
read_restart heatup.res.10000
   orthogonal box = (-41.1127 -42.0356 -5.19833) to (43.7383 42.8154 5.00473)
   4 by 4 by 1 processor grid
   70000 atoms
   69800 bonds
   69600 angles
   2 = max # of 1-2 neighbors
   2 = max # of special neighbors
reset_timestep 0
group ends type 2
400 atoms in group ends
group notends type 1
69600 atoms in group notends
bond_style fene
bond_coeff 1 30 1.5 1.0 1
special_bonds 0 1 1
angle_style cosine
angle_coeff 1 0.75
pair_style lj/cut 1.5
pair_modify shift yes
pair_coeff * * 1.0 1.0 1.5
fix 1 all npt 0.4 0.4 1 xyz 0 0 10
Resetting global state of Fix 1 Style npt from restart file info
timestep .01
thermo_style custom step temp epair ebond etotal press pxx pyy pzz lx ly lz vol
thermo 100
restart 500000 heatup.res
run 10000000
Memory usage per processor = 3.9225 Mbytes
Step Temp E_pair E_bond TotEng Press Pxx Pyy Pzz Lx Ly Lz Volume
        0 0.3987103 -2.4078843 20.357496 18.838813 0.018431475 -0.13542447 -0.15594064 0.34665954 84.851011 84.851011 10.203057 73458.892
      100 0.40108556 -2.4081637 20.355631 18.8401 -0.030417656 -0.25440208 -0.12787674 0.29102585 84.885421 84.885421 10.207537 73550.766
      200 0.3995556 -2.4056779 20.355978 18.840843 -0.027694001 -0.20666648 -0.18260281 0.30618728 84.885885 84.885885 10.207936 73554.44
      300 0.40283211 -2.4083831 20.354151 18.841343 0.0071930841 -0.18056975 -0.13478428 0.33693328 84.874065 84.874065 10.206857 73526.183

(snip)

9999800 0.40074939 -2.3910793 20.3509 18.890757 0.062947775 0.096553237 0.088241034 0.0040490548 27.80659 27.80659 95.77338 74052.595
9999900 0.40154465 -2.3907971 20.348669 18.891117 -0.0041623688 0.068006223 0.0082873043 -0.088780634 27.808561 27.808561 95.783381 74070.825
10000000 0.40053655 -2.3903391 20.349659 18.891075 -0.031163708 0.015608897 0.0042295935 -0.11332961 27.809911 27.809911 95.791244 74084.099
Loop time of 87932.2 on 16 procs for 10000000 steps with 70000 atoms

So, as you see, the initial Lx Ly Lz are: 84.851011 84.851011 10.203057
and the final Lx Ly Lz are: 27.809911 27.809911 95.791244

So the shape of the box is changing, which I do not want, and I did not think would happen with the "xyz" option.

Any clues as to what I did wrong and how to correct it?

Thanks,
Rob

2

Hi,
From what i understand, since you are running NPT, volume is is free to change. Now to keep the two ratio’s fixed leaves you with one cell length to be varied and rest to be adjusted as per the given ratios. I think you can use combination of deform and variable command to deal with this along with the npt fix to force the ratios as per requirement.

Thanks
Arnab

Arnab Chakrabarty
Graduate Student
Texas A & M University

3

If you have no other fix that changes the volume, then you're
right, the fix npt xyz style should keep the aspect ratio
of the box constant. Can you print out the 3 values
of the dilation[] variable in fix_npt.cpp as you run and
verify they are all the same each step? If so, each box length
is scaled by the same factor and their ratios should remain
constant.

Steve

4

(SEE THE "P. S." AT THE BOTTOM OF THIS MESSAGE FIRST)

Somehow, in fix_npt.cpp, a constant is being added to omega_dot[2]. Adding a line

for (i = 0; i < 3; i++) { p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
omega_dot[i] += f_omega*dthalf; omega_dot[i] *= drag_factor;
omega[i] += dtv*omega_dot[i];
factor[i] = exp(-dthalf*(eta_dot+omega_dot[i]));
dilation[i] = exp(dthalf*omega_dot[i]);
/* ADDED */ if (comm->me == 0 && update->ntimestep % 25 == 0) fprintf(screen, "Step %d, direction %d, f_omega %f, omega_dot %f, dilation-1 %e\n", update->ntimestep, i, f_omega, omega_dot[i], dilation[i]-1.0); }

in initial_integrate() produces in the logfile

Step Temp E_pair E_bond TotEng Press Pxx Pyy Pzz Lx Ly Lz Volume
        0 0.3987103 -2.4078843 20.357496 18.838813 0.018431475 -0.13542447 -0.15594064 0.34665954 84.851011 84.851011 10.203057 73458.892
Step 25, direction 0, f_omega -0.000901, omega_dot 0.000309, dilation-1 1.546047e-06
Step 25, direction 1, f_omega -0.000901, omega_dot 0.000309, dilation-1 1.546047e-06
Step 25, direction 2, f_omega -0.000901, omega_dot 0.000343, dilation-1 1.713797e-06
Step 50, direction 0, f_omega 0.000443, omega_dot 0.000436, dilation-1 2.180021e-06
Step 50, direction 1, f_omega 0.000443, omega_dot 0.000436, dilation-1 2.180021e-06
Step 50, direction 2, f_omega 0.000443, omega_dot 0.000470, dilation-1 2.347771e-06
Step 75, direction 0, f_omega -0.000509, omega_dot 0.000464, dilation-1 2.321444e-06
Step 75, direction 1, f_omega -0.000509, omega_dot 0.000464, dilation-1 2.321444e-06
Step 75, direction 2, f_omega -0.000509, omega_dot 0.000498, dilation-1 2.489194e-06
Step 100, direction 0, f_omega -0.000655, omega_dot 0.000324, dilation-1 1.620370e-06
Step 100, direction 1, f_omega -0.000655, omega_dot 0.000324, dilation-1 1.620370e-06
Step 100, direction 2, f_omega -0.000655, omega_dot 0.000358, dilation-1 1.788120e-06
      100 0.40108556 -2.4081637 20.355631 18.8401 -0.030417656 -0.25440208 -0.12787674 0.29102585 84.885421 84.885421 10.207537 73550.766

So, for some reason omega_dot[2] is always larger by about .000034. Puzzling, since in the constructor (line 192) there is

   omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0;

and I verified that this line of code executes. Somewhere between the constructor and the beginning of init(), the constant is added to omega_dot[2].

Constructor Omega_dot: x 0.000000 y 0.000000 z 0.000000
Begin Init Omega_dot: x 0.000302 y 0.000302 z 0.000335

I don't know how this could happen, since I use only one fix. Putting in the "sledgehammer" fix

omega_dot[2] = omega_dot[1] = omega_dot[0];

later fixes the problem, but I don't like this, since I don't know whether .000302 or .000335 is the correct value. I attach the fix_npt code with my extra output commands.

Thanks,
Rob

fix_npt.cpp (24.6 KB)

5

P. S. I found a workaround which eliminates the need for the
"sledgehammer" fix. If I run for 0 steps, unfix npt, refix npt, and
run, omega_dot is reset to zero and all three components remain the
same (as they change with time). So I guess the offset came from
reading in the restart file, which was generated by some older,
nonstandard code. But this still (might) indicate an undesirable
"leak".

The 3 omega and 3 omega_dot factors are saved in a restart file, so a fix npt
specified the same way can pick up where it left off. As described
on the read_restart doc page, the fix-ID and group-ID of the fix npt command
in your new script would have to match the old one and LAMMPS will
print a message at start-up indicating the fix is being re-enabled.

If this is happening, then it is likely the source of the anisotropy. If
not, and there is still an issue, let me know.

Steve