Commit d0235dd0 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

remove cauchystat changes from fix_nh.cpp/.h/.txt (again). delete misplaced... 

remove cauchystat changes from fix_nh.cpp/.h/.txt (again). delete misplaced example that came back from the dead, too.
parent a4835fa7

doc/src/fix_nh.txt

+0 −41
Original line number Diff line number Diff line
@@ -54,9 +54,6 @@ 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}
@@ -609,39 +606,6 @@ 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
@@ -659,7 +623,6 @@ 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.
@@ -681,7 +644,3 @@ 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/Mishin-Ni-Al-2009.eam.alloy

deleted100644 → 0
+0 −70007

File deleted.

Preview size limit exceeded, changes collapsed.

examples/cauchystat/README

deleted100644 → 0
+0 −10
Original line number Diff line number Diff line 
Run this example by executing:



% lmp_serial < lammps.in



Note that this example use an EAM potential, and therefore must be

run with a LAMMPS executable built with the MANYBODY package.



The first cauchystat fix equilibrates the temperature at zero stress,

the second fix applies a shear stress. Output in avg.txt shows

convergence to correct values.

examples/cauchystat/lammps.in

deleted100644 → 0
+0 −63
Original line number Diff line number Diff line 
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

src/fix_nh.cpp

+1 −365
Original line number Diff line number Diff line
@@ -59,6 +59,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : 
{

  if (narg < 4) error->all(FLERR,"Illegal fix nvt/npt/nph command");



  restart_global = 1;

  dynamic_group_allow = 1;

  time_integrate = 1;

  scalar_flag = 1;
@@ -67,12 +68,6 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : 
  extscalar = 1;

  extvector = 0;



  // for CauchyStat



  usePK = 1;

  restart_global = 1;

  restart_stored = 0;



  // default values



  pcouple = NONE;
@@ -346,16 +341,6 @@ 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");

      CauchyStat_init();

      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]);
@@ -2245,20 +2230,6 @@ 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]);



  // CauchyStat: call CauchyStat to modify p_target[i] and p_hydro,

  // if CauchyStat enabled and pressure->vector computation has been initiated



  if ((usePK == 0) && (initRUN == 1)) CauchyStat();

  if (initRUN == 0) {

    for (int i=0; i < 6; i++) {

      h_old[i]=domain->h[i];

    }

  }



  // when run is initialized tensor[] not available (pressure on cell wall)



  initRUN=1;



  // if deviatoric, recompute sigma each time p_target changes



  if (deviatoric_flag) compute_sigma();
@@ -2410,338 +2381,3 @@ double FixNH::memory_usage() 
  if (irregular) bytes += irregular->memory_usage();

  return bytes;

}
 


/* ----------------------------------------------------------------------

   initialize Cauchystat

------------------------------------------------------------------------- */



void FixNH::CauchyStat_init()

{

  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



  double *h = domain->h;

  double *h_inv = domain->h_inv;



  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];



  CSvol0 = abs(MathExtra::det3(H0));   //find reference volume



#undef H0

#undef invH0

}







/* ----------------------------------------------------------------------

   Cauchystat

------------------------------------------------------------------------- */



void FixNH::CauchyStat()

{

  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,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};

  double Fi[3][3]={{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};

  double Fdot[3][3]={{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};

  double J,vol;



  MathExtra::times3(H,invH0,F);

  MathExtra::times3(Hdot,invH0,Fdot);

  MathExtra::invert3(F,Fi);

  J = MathExtra::det3(F);

  vol=CSvol0*J;



  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];

  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];

  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_Step(Fi,Fdot,cauchy,setcauchy,setPK,vol,CSvol0,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



}



/* ----------------------------------------------------------------------

   CauchyStat_Step



   Inputs:

   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

   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_Step(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)

{



  //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;



  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) );

	}

      }

    }

  }



  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;

    }

  }



  for(int ii = 1;ii <= 6;ii++) {

    for(int jj = 1;jj <= 6;jj++) {

      deltaF(ii)=deltaF(ii)+dsdf(ii,jj)*fdotvec(jj);

    }

  }



  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(ii)=deltaPK(ii)+alpha*deltaF(ii);

  }



  setPK(1,1)=setPK(1,1)+deltaPK(1);

  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);



#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



}

CRA Git | Maintained and supported by SUSTech CRA and CCSE