Changes to the NH fix enabling Cauchy stress control (Cauhchystat) due to Miller, Tadmor, Gibson, Bernstein and Pavia, J Chem Phys,

144, 184107 (2016).

doc/src/fix_nh.txt

@@ -54,6 +54,9 @@ keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {y
   {scaleyz} value = {yes} or {no} = scale yz with lz
   {scalexz} value = {yes} or {no} = scale xz with lz
   {flip} value = {yes} or {no} = allow or disallow box flips when it becomes highly skewed
+  {cauchystat} cauchystat values = alpha continue
+    alpha = strength of Cauchystat control parameter
+    continue = {yes} or {no} = whether of not to continue from a previous run
   {fixedpoint} values = x y z
     x,y,z = perform barostat dilation/contraction around this point (distance units)
   {update} value = {dipole} or {dipole/dlm}
@@ -606,6 +609,39 @@ can only be used if the 2nd dimension in the keyword is periodic,
 and if the tilt factor is not coupled to the barostat via keywords
 {tri}, {yz}, {xz}, and {xy}.
 
+Without the {cauchystat} keyword, the barostat algorithm
+controls the Second-Piola Kirchhoff stress, which is a stress measure
+referred to the undeformed (initial) simulation box.  If the box
+deforms substantially during the equilibration, the difference between
+the set values and the final true (Cauchy) stresses can be
+considerable.
+
+The {cauchystat} keyword modifies the barostat as per Miller et
+al. (Miller)_"#nh-Miller" so that the Cauchy stress is controlled.
+{alpha} is the non-dimensional parameter, typically set to 0.001 or
+0.01 that determines how aggresively the algorithm drives the system
+towards the set Cauchy stresses.  Larger values of {alpha} will modify
+the system more quickly, but can lead to instabilities.  Smaller
+values will lead to longer convergence time.  Since {alpha} also
+influences how much the stress fluctuations deviate from the
+equilibrium fluctuations, it should be set as small as possible.
+
+A {continue} value of {yes} indicates that the fix is subsequent to a
+previous run with the Cauchystat fix, and the intention is to continue
+from the converged stress state at the end of the previous run.  This
+may be required, for example, when implementing a multi-step loading/unloading
+sequence over several fixes.
+
+Setting {alpha} to zero is not permitted.  To "turn off" the
+Cauchystat control and thus restore the equilibrium stress
+fluctuations, two subsequent fixes should be used.  In the first, the
+Cauchystat is used and the simulation box equilibrates to the correct
+shape for the desired stresses.  In the second, the {fix} statement is
+identical except that the {cauchystat} keyword is removed (along with
+related {alpha} and {continue} values). This restores the original
+Parrinello-Rahman algorithm, but now with the correct simulation box
+shape from the first fix.
+
 These fixes can be used with dynamic groups as defined by the
 "group"_group.html command.  Likewise they can be used with groups to
 which atoms are added or deleted over time, e.g. a deposition
@@ -623,6 +659,7 @@ over time or the atom count becomes very small.
 
 The keyword defaults are tchain = 3, pchain = 3, mtk = yes, tloop =
 ploop = 1, nreset = 0, drag = 0.0, dilate = all, couple = none,
+cauchystat = no,
 scaleyz = scalexz = scalexy = yes if periodic in 2nd dimension and
 not coupled to barostat, otherwise no.
 
@@ -644,3 +681,7 @@ Martyna, J Phys A: Math Gen, 39, 5629 (2006).
 :link(nh-Dullweber)
 [(Dullweber)] Dullweber, Leimkuhler and McLachlan, J Chem Phys, 107,
 5840 (1997).
+
+:link(nh-Miller)
+[(Miller)] Miller, Tadmor, Gibson, Bernstein and Pavia, J Chem Phys,
+144, 184107 (2016).

examples/cauchystat/Al_shear_CS/lammps.in

@@ -0,0 +1,69 @@
+units		metal
+atom_style	atomic
+atom_modify	map array
+
+processors      1 1 1
+
+# Box and atom positions:
+boundary p p p
+
+# Defining lattice and creating simulation
+# box with atoms inside
+lattice          fcc 4.05
+region           simbox prism 0 5 0 5 0 5 0 0 0 units lattice
+create_box       2 simbox
+create_atoms     2 box
+
+# Atomic mass:
+mass 1 58.69
+mass 2 26.98154
+
+# Potential, Al fcc crystal
+pair_style eam/alloy
+pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
+
+
+thermo 100
+thermo_style custom step temp pxx pyy pzz pxy pxz pyz
+compute cna all cna/atom 2.8
+
+fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 couple none cauchystat 0.001 no
+
+dump 1 all cfg 1000 test*.cfg mass type xs ys zs type c_cna
+
+timestep 0.002
+
+variable px equal pxx
+variable py equal pyy
+variable pz equal pzz
+variable sxy equal pxy
+variable sxz equal pxz
+variable syz equal pyz
+variable t equal temp
+
+fix avg all ave/time 1 100 100 v_t v_px v_py v_pz v_sxy v_sxz v_syz file avg.txt
+
+variable lx equal lx
+variable ly equal ly
+variable lz equal ly
+variable xy equal xy
+variable xz equal xz
+variable yz equal yz
+
+fix box all ave/time 1 100 100 v_lx v_ly v_lz v_xy v_xz v_yz file box.txt
+
+velocity all create 1200 4928459 rot yes dist gaussian
+
+run 10000
+
+unfix 1
+
+fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 xy -10000.0 -10000.0 0.1 couple none cauchystat 0.001 yes
+
+run 10000
+
+unfix 1
+
+fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1  xy -10000.0 -10000.0 0.1 couple none
+
+run 10000

examples/cauchystat/Al_shear_PR/lammps.in

@@ -0,0 +1,69 @@
+units		metal
+atom_style	atomic
+atom_modify	map array
+
+processors      1 1 1
+
+# Box and atom positions:
+boundary p p p
+
+# Defining lattice and creating simulation
+# box with atoms inside
+lattice          fcc 4.05
+region           simbox prism 0 5 0 5 0 5 0 0 0 units lattice
+create_box       2 simbox
+create_atoms     2 box
+
+# Atomic mass:
+mass 1 58.69
+mass 2 26.98154
+
+# Potential, Al fcc crystal
+pair_style eam/alloy
+pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
+
+
+thermo 100
+thermo_style custom step temp pxx pyy pzz pxy pxz pyz
+compute cna all cna/atom 2.8
+
+fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 couple none
+
+dump 1 all cfg 1000 test*.cfg mass type xs ys zs type c_cna
+
+timestep 0.002
+
+variable px equal pxx
+variable py equal pyy
+variable pz equal pzz
+variable sxy equal pxy
+variable sxz equal pxz
+variable syz equal pyz
+variable t equal temp
+
+fix avg all ave/time 1 100 100 v_t v_px v_py v_pz v_sxy v_sxz v_syz file avg.txt
+
+variable lx equal lx
+variable ly equal ly
+variable lz equal ly
+variable xy equal xy
+variable xz equal xz
+variable yz equal yz
+
+fix box all ave/time 1 100 100 v_lx v_ly v_lz v_xy v_xz v_yz file box.txt
+
+velocity all create 1200 4928459 rot yes dist gaussian
+
+run 10000
+
+unfix 1
+
+fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 xy -10000.0 -10000.0 0.1 couple none
+
+run 10000
+
+unfix 1
+
+fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1  xy -10000.0 -10000.0 0.1 couple none
+
+run 10000

examples/cauchystat/NiAl_cool/lammps.in

@@ -0,0 +1,45 @@
+units		metal
+atom_style	atomic
+atom_modify	map array
+
+processors      1 1 1
+
+# Box and atom positions:
+boundary p p p
+read_data input.dat
+
+# Atomic mass:
+mass 1 58.69
+mass 2 26.98154
+
+# Potential, Al fcc crystal
+pair_style eam/alloy
+pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
+
+thermo 100
+thermo_style custom step temp pxx pyy pzz lx ly lz
+compute cna all cna/atom 2.8
+
+velocity all create 2400.0 4928459 rot yes dist gaussian
+
+fix 1 all npt temp 1200.0 1200.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z -2800.0 -2800.0 0.1 couple none cauchystat 0.001 no
+
+dump 1 all cfg 1000 test*.cfg mass type xs ys zs type c_cna
+
+timestep 0.002
+
+variable l equal lz
+variable p equal pzz
+variable t equal temp
+
+fix avg all ave/time 1 1000 1000 v_t v_l v_p file avg.txt
+
+velocity all create 2400 4928459 rot yes dist gaussian
+
+run 10000
+
+unfix 1
+
+fix 1 all npt temp 1200.0 300.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z -2800.0 -2800.0 0.1 couple none cauchystat 0.001 yes
+
+run 100000

examples/cauchystat/NiAl_slow_stretch_CS/lammps.in

@@ -0,0 +1,45 @@
+units		metal
+atom_style	atomic
+atom_modify	map array
+
+processors      1 1 1
+
+# Box and atom positions:
+boundary p p p
+read_data input.dat
+
+# Atomic mass:
+mass 1 58.69
+mass 2 26.98154
+
+# Potential, Al fcc crystal
+pair_style eam/alloy
+pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
+
+thermo 100
+thermo_style custom step temp pxx pyy pzz lx ly lz
+compute cna all cna/atom 2.8
+
+velocity all create 1000.0 4928459 rot yes dist gaussian
+
+fix 1 all npt temp 500.0 500.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 couple none cauchystat 0.001 no
+
+dump 1 all cfg 5000 test*.cfg mass type xs ys zs type c_cna
+
+timestep 0.002
+
+variable l equal lz
+variable p equal pzz
+variable t equal temp
+
+fix avg all ave/time 1 1000 1000 v_t v_l v_p file avg.txt
+
+velocity all create 2400 4928459 rot yes dist gaussian
+
+run 10000
+
+unfix 1
+
+fix 1 all npt temp 500.0 500.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 -14000.0 0.1 couple none cauchystat 0.001 yes
+
+run 350000

examples/cauchystat/NiAl_slow_stretch_PR/lammps.in

@@ -0,0 +1,45 @@
+units		metal
+atom_style	atomic
+atom_modify	map array
+
+processors      1 1 1
+
+# Box and atom positions:
+boundary p p p
+read_data input.dat
+
+# Atomic mass:
+mass 1 58.69
+mass 2 26.98154
+
+# Potential, Al fcc crystal
+pair_style eam/alloy
+pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
+
+thermo 100
+thermo_style custom step temp pxx pyy pzz lx ly lz
+compute cna all cna/atom 2.8
+
+velocity all create 1000.0 4928459 rot yes dist gaussian
+
+fix 1 all npt temp 500.0 500.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 couple none
+
+dump 1 all cfg 5000 test*.cfg mass type xs ys zs type c_cna
+
+timestep 0.002
+
+variable l equal lz
+variable p equal pzz
+variable t equal temp
+
+fix avg all ave/time 1 1000 1000 v_t v_l v_p file avg.txt
+
+velocity all create 2400 4928459 rot yes dist gaussian
+
+run 10000
+
+unfix 1
+
+fix 1 all npt temp 500.0 500.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 -14000.0 0.1 couple none
+
+run 350000

examples/cauchystat/NiAl_stretch/lammps.in

@@ -0,0 +1,48 @@
+units		metal
+atom_style	atomic
+atom_modify	map array
+
+processors      1 1 1
+
+# Box and atom positions:
+boundary p p p
+read_data input.dat
+
+# Atomic mass:
+mass 1 58.69
+mass 2 26.98154
+
+# Potential, Al fcc crystal
+pair_style eam/alloy
+pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
+
+
+thermo 100
+thermo_style custom step temp pxx pyy pzz lx ly lz
+compute cna all cna/atom 2.8
+
+velocity all create 2400.0 4928459 rot yes dist gaussian
+
+fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 couple none cauchystat 0.001 no
+
+dump 1 all cfg 1000 test*.cfg mass type xs ys zs type c_cna
+
+timestep 0.002
+
+variable l equal lz
+variable p equal pzz
+variable t equal temp
+
+fix avg all ave/time 1 1000 1000 v_t v_l v_p file avg.txt
+
+#restart 1000 restart.*.test
+
+velocity all create 2400 4928459 rot yes dist gaussian
+
+run 10000
+
+unfix 1
+
+fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 -6000.0 0.1 couple none cauchystat 0.001 yes
+
+run 100000

examples/cauchystat/README

@@ -0,0 +1,42 @@
+There are 6 examples in this directory.  All are run by executing:
+
+% lmp_serial < lammps.in
+
+from the respective home directory.
+
+Note that these examples use an EAM potential, and therefore must be
+run with a LAMMPS executable built with the MANYBODY package.
+
+Al_shear_PR: Application of a simple shear stress, illustrating that
+there is a non-zero direct stress in the equilibrated system due to
+the fact that Parinello-Rahman controls the second Piola-Kirchhoff
+stress instead of the Cauchy stress.  First fix equilibrates the
+temperature at zero stress, second fix applies the shear.  The third
+fix corrects the problem because each new run resets the reference
+simulation cell (H0) to the current cell shape.
+
+Al_shear_CS: Cauchystat example. Same as above, using the Cauchystat
+instead of PR control.  In this case, the equilibrated stresses are as
+prescribed (ie., the set stresses control the Cauchy stress, not Piola
+Kirchhoff)
+
+NiAl_cool: Cauchystat example. First fix equilibrates the system to
+uniaxial tension of 2800 bars, 1200K.  Second fix cools to 300 K,
+causing a phase transformation.  Cauchy stress remains at the set
+level after the transformation
+
+NiAl_stretch: Cauchystat example.  After a brief equilibration
+at 600 K and zero stress, quickly ramp up the applied tensile stress
+from 0 to 6000 bars over 100000 time steps.  After the phase
+transformation, Cauchy stress continues to follow the set value.
+
+NiAl_slow_stretch_CS: Cauchystat example.  After a brief equilibration
+at 500 K and zero stress, gradually ramp up the applied tensile stress
+from 0 to 14000 bars over 350000 time steps.  After the phase
+transformation, Cauchy stress continues to follow the set value.
+
+NiAl_slow_stretch_PR: After a brief equilibration at 500 K and zero
+stress, gradually ramp up the applied tensile stress from 0 to 14000
+bars over 350000 time steps.  After the phase transformation, Cauchy
+stress departs substantially from the set value, because the PR
+control is used instead of the Cauchystat.

src/fix_nh.cpp

@@ -46,6 +46,10 @@ enum{NOBIAS,BIAS};
 enum{NONE,XYZ,XY,YZ,XZ};
 enum{ISO,ANISO,TRICLINIC};
 
+int FixNH::restartPK;
+int FixNH::restart_stored=0;
+double FixNH::setPKinit[6];
+
 /* ----------------------------------------------------------------------
    NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion
  ---------------------------------------------------------------------- */
@@ -59,6 +63,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg < 4) error->all(FLERR,"Illegal fix nvt/npt/nph command");
 
+  usePK = 1;
   restart_global = 1;
   dynamic_group_allow = 1;
   time_integrate = 1;
@@ -341,6 +346,78 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
         dlm_flag = 1;
       } else error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"cauchystat") == 0) {
+      if (iarg+3 > narg) error->all(FLERR,"Illegal fix npt cauchystat command: wrong number of arguments");
+      if (strcmp(arg[iarg+2],"yes") != 0 && strcmp(arg[iarg+2],"no") != 0) error->all(FLERR,"Illegal cauchystat continue value.  Must be 'yes' or 'no'");
+      alpha = force->numeric(FLERR,arg[iarg+1]);
+      restartPK = !strcmp(arg[iarg+2],"yes");
+
+
+      if (comm->me == 0) {
+	if (screen) {
+	  fprintf(screen,"Using the Cauchystat fix with alpha=%f\n",alpha);
+	  if(restartPK==1) {
+	    fprintf(screen,"   (this is a continuation fix)\n");
+	  }
+	  else {
+	    fprintf(screen,"   (this is NOT a continuation fix)\n");
+	  }
+	}
+	if (logfile) {
+	  fprintf(logfile,"Using the Cauchystat with alpha=%f\n",alpha);
+	  if(restartPK==1) {
+	    fprintf(logfile,"   this is a continuation run\n");
+	  }
+	  else {
+	    fprintf(logfile,"   this is NOT a continuation run\n");
+	  }
+	}
+      }
+
+
+      if(restartPK==1 && restart_stored==0) {
+	error->all(FLERR,"Illegal cauchystat command.  Continuation run must follow a previously equilibrated Cauchystat run");
+      }
+
+      if(alpha<=0.0) {
+	error->all(FLERR,"Illegal cauchystat command.  Alpha cannot be zero or negative");
+      }
+
+      initRUN = 0;
+      initPK = 1;
+      usePK = 0;
+
+      #define H0(row,col) (H0[(row-1)][(col-1)])
+      #define invH0(row,col) (invH0[(row-1)][(col-1)])
+
+
+      // initialize H0 to current box shape
+
+      int triclinic = domain->triclinic;
+      double *h = domain->h;
+      double *h_inv = domain->h_inv;
+      double *boxhi = domain->boxhi;
+      double *boxlo = domain->boxlo;
+      double yz = domain->yz;
+      double xz = domain->xz;
+      double xy = domain->xy;
+
+
+      H0(1,1)=h[0];  H0(1,2)=h[5];  H0(1,3)=h[4];
+      H0(2,1)=0.0;   H0(2,2)=h[1];  H0(2,3)=h[3];
+      H0(3,1)=0.0;   H0(3,2)=0.0;   H0(3,3)=h[2];
+
+      invH0(1,1)=h_inv[0];  invH0(1,2)=h_inv[5];  invH0(1,3)=h_inv[4];
+      invH0(2,1)=0.0;       invH0(2,2)=h_inv[1];  invH0(2,3)=h_inv[3];
+      invH0(3,1)=0.0;       invH0(3,2)=0.0;       invH0(3,3)=h_inv[2];
+
+      myvol0 = abs(MathExtra::det3(H0));   //find reference volume
+
+      #undef H0
+      #undef invH0
+
+      iarg += 3;
+
     } else if (strcmp(arg[iarg],"fixedpoint") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       fixedpoint[0] = force->numeric(FLERR,arg[iarg+1]);
@@ -2230,8 +2307,170 @@ void FixNH::compute_press_target()
     for (int i = 3; i < 6; i++)
       p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
 
-  // if deviatoric, recompute sigma each time p_target changes
+  // CauchyStat: call CauchyStat to modify p_target[i] and p_hydro, they are used in compute_sigma()
+  //             if CauchyStat enabled and pressure->vector computation has been initiated
+  if( (usePK==0) && (initRUN==1) ) {
 
+    double* h = domain->h;           // shape matrix in Voigt notation
+    double* h_rate = domain->h_rate; // rate of box size/shape change in Voigt notation
+    double H[3][3];
+    double Hdot[3][3];
+
+    #define H(row,col) (H[(row-1)][(col-1)])
+    #define Hdot(row,col) (Hdot[(row-1)][(col-1)])
+    #define F(row,col) (F[(row-1)][(col-1)])
+    #define Fi(row,col) (Fi[(row-1)][(col-1)])
+    #define Fdot(row,col) (Fdot[(row-1)][(col-1)])
+    #define invH0(row,col) (invH0[(row-1)][(col-1)])
+    #define cauchy(row,col) (cauchy[(row-1)][(col-1)])
+    #define setcauchy(row,col) (setcauchy[(row-1)][(col-1)])
+    #define setPK(row,col) (setPK[(row-1)][(col-1)])
+    #define sigmahat(row,col) (sigmahat[(row-1)][(col-1)])
+
+    H(1,1)=h[0]; H(1,2)=0.0;  H(1,3)=0.0;
+    H(2,1)=0.0;  H(2,2)=h[1];  H(2,3)=0.0;
+    H(3,1)=0.0;  H(3,2)=0.0;   H(3,3)=h[2];
+
+    for (int i=0;i<6;i++) {
+      h_rate[i]=(h[i]-h_old[i])/update->dt;
+      h_old[i]=h[i];
+    }
+
+    Hdot(1,1)=h_rate[0]; Hdot(1,2)=0.0;        Hdot(1,3)=0.0;
+    Hdot(2,1)=0.0;       Hdot(2,2)=h_rate[1];  Hdot(2,3)=0.0;
+    Hdot(3,1)=0.0;       Hdot(3,2)=0.0;        Hdot(3,3)=h_rate[2];
+
+    if (domain->triclinic) {
+      H(1,2)=h[5];  H(1,3)=h[4]; H(2,3)=h[3];
+      Hdot(1,2)=h_rate[5];  Hdot(1,3)=h_rate[4]; Hdot(2,3)=h_rate[3];
+    }
+
+    double F[3][3]={0.0};
+    double FT[3][3]={0.0};
+    double Fi[3][3]={0.0};
+    double Fdot[3][3]={0.0};
+    double J,vol;
+
+    MathExtra::times3(H,invH0,F);       //find F
+    MathExtra::times3(Hdot,invH0,Fdot); //find Fdot
+
+    MathExtra::invert3(F,Fi);
+    J = MathExtra::det3(F);   //find J
+    vol=myvol0*J;             // actual volume
+
+    double deltat;
+    deltat = update->dt; //increment of time step
+
+    // Current pressure on the cell boundaries:
+    //tensor:  0    1    2    3    4    5
+    //         x    y    z    xy   xz   yz
+
+    double *tensor = pressure->vector;
+
+    double cauchy[3][3]; //stress
+    cauchy(1,1)=-tensor[0]; cauchy(1,2)=0.0;        cauchy(1,3)=0.0;
+    cauchy(2,1)=0.0;        cauchy(2,2)=-tensor[1]; cauchy(2,3)=0.0;
+    cauchy(3,1)=0.0;        cauchy(3,2)=0.0;        cauchy(3,3)=-tensor[2];
+
+    if (domain->triclinic) {
+      cauchy(1,2)=-tensor[3]; cauchy(1,3)=-tensor[4];
+      cauchy(2,1)=-tensor[3]; cauchy(2,3)=-tensor[5];
+      cauchy(3,1)=-tensor[4]; cauchy(3,2)=-tensor[5];
+    }
+
+    // target pressure on the cell boundaries:
+    //p_target:          0          1          2          3          4          5
+    //                   x          y          z         yz         xz         xy
+
+    double setcauchy[3][3]; //stress
+    setcauchy(1,1)=-p_target[0]; setcauchy(1,2)=0.0;          setcauchy(1,3)=0.0;
+    setcauchy(2,1)=0.0;          setcauchy(2,2)=-p_target[1]; setcauchy(2,3)=0.0;
+    setcauchy(3,1)=0.0;          setcauchy(3,2)=0.0;          setcauchy(3,3)=-p_target[2];
+
+    if (domain->triclinic) {
+      setcauchy(1,2)=-p_target[5]; setcauchy(1,3)=-p_target[4];
+      setcauchy(2,1)=-p_target[5]; setcauchy(2,3)=-p_target[3];
+      setcauchy(3,1)=-p_target[4]; setcauchy(3,2)=-p_target[3];
+    }
+
+    //initialize:
+    if(initPK==1) {
+      if(restartPK==1) {
+        setPK(1,1)=setPKinit[0]; setPK(1,2)=setPKinit[1]; setPK(1,3)=setPKinit[2];
+        setPK(2,1)=setPKinit[1]; setPK(2,2)=setPKinit[3]; setPK(2,3)=setPKinit[4];
+        setPK(3,1)=setPKinit[2]; setPK(3,2)=setPKinit[4]; setPK(3,3)=setPKinit[5];
+      }else {
+        setPK(1,1)=cauchy(1,1); setPK(1,2)=cauchy(1,2); setPK(1,3)=cauchy(1,3);
+        setPK(2,1)=cauchy(2,1); setPK(2,2)=cauchy(2,2); setPK(2,3)=cauchy(2,3);
+        setPK(3,1)=cauchy(3,1); setPK(3,2)=cauchy(3,2); setPK(3,3)=cauchy(3,3);
+      }
+
+      initPK=0;
+    }
+
+    //cauchystat:
+    bigint step = update->ntimestep;
+
+
+
+
+
+    CauchyStat(step,F,Fi,Fdot,cauchy,setcauchy,setPK,vol,myvol0,deltat,alpha);
+
+    // use currentPK as new target:
+    //p_target:          0          1          2          3          4          5
+    //                   x          y          z         yz         xz         xy
+
+    p_target[0]=-setPK(1,1);
+    p_target[1]=-setPK(2,2);
+    p_target[2]=-setPK(3,3);
+    if (pstyle == TRICLINIC) {
+      p_target[3]=-setPK(2,3);
+      p_target[4]=-setPK(1,3);
+      p_target[5]=-setPK(1,2);
+    }
+
+    p_hydro = 0.0;
+    for (int i = 0; i < 3; i++)
+      if (p_flag[i]) {
+        p_hydro += p_target[i];
+      }
+    p_hydro /= pdim;
+
+    // save information for Cauchystat restart
+    setPKinit[0] = setcauchy(1,1);
+    setPKinit[1] = setcauchy(1,2);
+    setPKinit[2] = setcauchy(1,3);
+    setPKinit[3] = setcauchy(2,2);
+    setPKinit[4] = setcauchy(2,3);
+    setPKinit[5] = setcauchy(3,3);
+    restart_stored=1;
+
+    #undef H
+    #undef Hdot
+    #undef F
+    #undef Fi
+    #undef Fdot
+    #undef invH0
+    #undef cauchy
+    #undef setcauchy
+    #undef setPK
+    #undef sigmahat
+
+  }
+
+  if(initRUN==0){
+    double* h = domain->h;           // shape matrix in Voigt notation
+    for (int i=0;i<6;i++) {
+      h_old[i]=h[i];
+    }
+  }
+
+
+  initRUN=1; // when run is initialized tensor[] not available (pressure on cell wall)
+
+
+  // if deviatoric, recompute sigma each time p_target changes
   if (deviatoric_flag) compute_sigma();
 }
 
@@ -2381,3 +2620,142 @@ double FixNH::memory_usage()
   if (irregular) bytes += irregular->memory_usage();
   return bytes;
 }
+
+/* ----------------------------------------------------------------------
+   CauchyStat
+
+   Inputs:
+   step             :  current step number
+   F(3,3)           :  current deformation gradient
+   Fi(3,3)          :  inverse of the deformation gradient
+   Fdot(3,3)        :  time derivative of the deformation gradient (velocity)
+   cauchy(3,3)      :  current cauchy stress tensor
+   setcauchy(3,3)   :  target cauchy stress tensor
+   alpha            :  parameter =0.01
+   nstat            :  =1, flag for mod(step,nstat).ne.0
+   setPK(3,3)       :  current PK stress tensor, at initialization (time=0) setPK=cauchy
+   volume           :  current volume of the periodic box
+   volume0          :  initial volume of the periodic box
+   deltat           :  simulation step increment
+   alpha            :  parameter
+
+   Outputs:
+   setPK(3,3)       :  PK stress tensor at the next time step
+------------------------------------------------------------------------- */
+
+void FixNH::CauchyStat(bigint step, double (&F)[3][3], double (&Fi)[3][3], double (&Fdot)[3][3], double (&cauchy)[3][3], double (&setcauchy)[3][3], double (&setPK)[3][3], double volume, double volume0, double deltat, double alpha)
+{
+
+  int nstat=1;
+
+  //macros to go from c to fortran style for arrays:
+  #define F(row,col) (F[(row-1)][(col-1)])
+  #define Fi(row,col) (Fi[(row-1)][(col-1)])
+  #define Fdot(row,col) (Fdot[(row-1)][(col-1)])
+  #define cauchy(row,col) (cauchy[(row-1)][(col-1)])
+  #define setcauchy(row,col) (setcauchy[(row-1)][(col-1)])
+  #define setPK(row,col) (setPK[(row-1)][(col-1)])
+
+  #define uv(row,col) (uv[(row-1)][(col-1)])
+  #define deltastress(row) (deltastress[(row-1)])
+  #define fdotvec(row) (fdotvec[(row-1)])
+  #define dsdf(row,col) (dsdf[(row-1)][(col-1)])
+  #define dsds(row,col) (dsds[(row-1)][(col-1)])
+  #define deltaF(row) (deltaF[(row-1)])
+  #define deltaPK(row) (deltaPK[(row-1)])
+
+  //initialize local variables:
+  int i,j,m,n;
+
+  double deltastress[6],fdotvec[6];
+  double dsdf[6][6]={0.0}; //zeroed, used in incremental form
+  double dsds[6][6]={0.0}; //zeroed, used in incremental form
+  double jac;
+  double deltaF[6]={0.0}; //zeroed, used in incremental form
+  double deltaPK[6]={0.0}; //zeroed, used in incremental form
+
+  int uv[6][2];
+  uv(1,1)=1; uv(1,2)=1;
+  uv(2,1)=2; uv(2,2)=2;
+  uv(3,1)=3; uv(3,2)=3;
+  uv(4,1)=2; uv(4,2)=3;
+  uv(5,1)=1; uv(5,2)=3;
+  uv(6,1)=1; uv(6,2)=2;
+
+  if( (step % nstat) != 0) {
+    //do nothing!
+  }else {
+    for(int ii = 1;ii <= 6;ii++) {
+      i=uv(ii,1);
+      j=uv(ii,2);
+      deltastress(ii)=setcauchy(i,j)-cauchy(i,j);
+      if(ii>3) deltastress(ii)=deltastress(ii)*2.0;
+      fdotvec(ii)=Fdot(i,j)*deltat;
+    }
+
+    for(int ii = 1;ii <= 6;ii++) {
+      i=uv(ii,1);
+      j=uv(ii,2);
+      for(int jj = 1;jj <= 6;jj++) {
+        m=uv(jj,1);
+        n=uv(jj,2);
+        dsds(ii,jj) = Fi(i,m)*Fi(j,n) + Fi(i,n)*Fi(j,m) + Fi(j,m)*Fi(i,n) + Fi(j,n)*Fi(i,m);
+        for(int l = 1;l <= 3;l++) {
+          for(int k = 1;k <= 3;k++) {
+            dsdf(ii,jj) = dsdf(ii,jj) + cauchy(k,l)*( Fi(i,k)*Fi(j,l)*Fi(n,m) - Fi(i,m)*Fi(j,l)*Fi(n,k) - Fi(i,k)*Fi(j,m)*Fi(n,l) );
+          }//k
+        }//l
+      }//jj
+    }//ii
+
+    jac=volume/volume0;
+    for(int ii = 1;ii <= 6;ii++) {
+      for(int jj = 1;jj <= 6;jj++) {
+        dsds(ii,jj)=dsds(ii,jj)*jac/4.0;
+      dsdf(ii,jj)=dsdf(ii,jj)*jac;
+      }//jj
+    }//ii
+
+    for(int ii = 1;ii <= 6;ii++) {
+      for(int jj = 1;jj <= 6;jj++) {
+        deltaF(ii)=deltaF(ii)+dsdf(ii,jj)*fdotvec(jj); // deltaF=matmul(dsdf,fdotvec) in the fortran implementation
+      }//jj
+    }//ii
+
+    for(int ii = 1;ii <= 6;ii++) {
+      for(int jj = 1;jj <= 6;jj++) {
+        deltaPK(ii)=deltaPK(ii)+alpha*dsds(ii,jj)*deltastress(jj); // deltaPK=alpha*matmul(dsds,deltastress) + deltaF in the fortran implementation
+      }//jj
+      // this line applies alpha to both terms instead of just the delta sigma term.
+      deltaPK(ii)=deltaPK(ii)+alpha*deltaF(ii);
+      // this is the old version
+      //deltaPK(ii)=deltaPK(ii)+deltaF(ii);
+    }//ii
+
+    setPK(1,1)=setPK(1,1)+deltaPK(1); //equation (4) in SD-notes.pdf
+    setPK(2,2)=setPK(2,2)+deltaPK(2);
+    setPK(3,3)=setPK(3,3)+deltaPK(3);
+    setPK(2,3)=setPK(2,3)+deltaPK(4);
+    setPK(3,2)=setPK(3,2)+deltaPK(4);
+    setPK(1,3)=setPK(1,3)+deltaPK(5);
+    setPK(3,1)=setPK(3,1)+deltaPK(5);
+    setPK(1,2)=setPK(1,2)+deltaPK(6);
+    setPK(2,1)=setPK(2,1)+deltaPK(6);
+  }
+
+  //undefine macros:
+  #undef F
+  #undef Fi
+  #undef Fdot
+  #undef cauchy
+  #undef setcauchy
+  #undef setPK
+  #undef uv
+  #undef deltastress
+  #undef fdotvec
+  #undef dsdf
+  #undef dsds
+  #undef deltaF
+  #undef deltaPK
+
+}

src/fix_nh.h

@@ -139,6 +139,23 @@ class FixNH : public Fix {
   double compute_strain_energy();
   void compute_press_target();
   void nh_omega_dot();
+
+  // Implementation of CauchyStat
+  double H0[3][3];               //shape matrix for the undeformed cell
+  double h_old[6];               //previous time step shape matrix for the undeformed cell
+  double invH0[3][3];            //inverse of H0;
+  double myvol0;
+  double setPK[3][3];
+  static double setPKinit[6];
+  double alpha;                  //integration parameter cauchystat
+  int initPK;                    // 1 if setPK needs to be initialized either from cauchy or restart, else 0
+  int usePK;                     // 0 if use CauchyStat else 1
+  static int restartPK;          // Read PK stress from the previous step
+  static int restart_stored;          // Read PK stress from the previous step
+  int initRUN;                   // 0 if run not initialized (pressure->vector not computed yet), else 1 (pressure->vector available)
+
+  virtual void CauchyStat(bigint step, double (&F)[3][3], double (&Fi)[3][3], double (&Fdot)[3][3], double (&cauchy)[3][3], double (&setcauchy)[3][3], double (&setPK)[3][3], double volume, double volume0, double deltat, double alpha);
+
 };
 
 }