changed computes names to stress/mop and stress/mop/profile

"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c

:link(lws,http://lammps.sandia.gov)

:link(ld,Manual.html)

:link(lc,Section_commands.html#comm)

:line

compute stress/mop command :h3

compute stress/mop/profile command :h3

[Syntax:]

compute ID group-ID style dir args keywords ... :pre

ID, group-ID are documented in "compute"_compute.html command

style = {stress/mop} or {stress/mop/profile}

dir = {x} or {y} or {z} is the direction normal to the plane

args = argument specific to the compute style

keywords = {kin} or {conf} or {total} (one of more can be specified) :ul

  {stress/mop} args = pos

    pos = {lower} or {center} or {upper} or coordinate value (distance units) is the position of the plane

  {stress/mop/profile} args = origin delta

    origin = {lower} or {center} or {upper} or coordinate value (distance units) is the position of the first plane

    delta = value (distance units) is the distance between planes :pre

compute 1 all stress/mop x lower total

compute 1 liquid stress/mop z 0.0 kin conf

fix 1 all ave/time 10 1000 10000 c_1\[*\] file mop.time

fix 1 all ave/time 10 1000 10000 c_1\[2\] file mop.time :pre

compute 1 all stress/mop/profile x lower 0.1 total

compute 1 liquid stress/mop/profile z 0.0 0.25 kin conf

fix 1 all ave/time 500 20 10000 c_1\[*\] ave running overwrite file mopp.time mode vector :pre

[Description:]

Compute {stress/mop} and compute {stress/mop/profile} define computations that

calculate components of the local stress tensor using the method of

planes "(Todd)"_#mop-todd.  Specifically in compute {stress/mop} calculates 3

components are computed in directions {dir},{x}; {dir},{y}; and

{dir},{z}; where {dir} is the direction normal to the plane, while

in compute {stress/mop/profile} the profile of the stress is computed.

Contrary to methods based on histograms of atomic stress (i.e. using

"compute stress/atom"_compute_stress_atom.html), the method of planes is

compatible with mechanical balance in heterogeneous systems and at

interfaces "(Todd)"_#mop-todd.

The stress tensor is the sum of a kinetic term and a configurational

term, which are given respectively by Eq. (21) and Eq. (16) in

"(Todd)"_#mop-todd. For the kinetic part, the algorithm considers that

atoms have crossed the plane if their positions at times t-dt and t are

one on either side of the plane, and uses the velocity at time t-dt/2

given by the velocity-Verlet algorithm.

Between one and three keywords can be used to indicate which

contributions to the stress must be computed: kinetic stress (kin),

configurational stress (conf), and/or total stress (total).

NOTE 1: The configurational stress is computed considering all pairs of atoms where at least one atom belongs to group group-ID. 

NOTE 2: The local stress does not include any Lennard-Jones tail

corrections to the pressure added by the "pair_modify tail

yes"_pair_modify.html command, since those are contributions to the global system pressure.

[Output info:]

Compute {stress/mop} calculates a global vector (indices starting at 1), with 3

values for each declared keyword (in the order the keywords have been

declared). For each keyword, the stress tensor components are ordered as

follows: stress_dir,x, stress_dir,y, and stress_dir,z.

Compute {stress/mop/profile} instead calculates a global array, with 1 column

giving the position of the planes where the stress tensor was computed,

and with 3 columns of values for each declared keyword (in the order the

keywords have been declared). For each keyword, the profiles of stress

tensor components are ordered as follows: stress_dir,x; stress_dir,y;

and stress_dir,z.

The values are in pressure "units"_units.html. 

The values produced by this compute can be accessed by various "output commands"_Howto_output.html. For instance, the results can be written to a file using the "fix ave/time"_fix_ave_time.html command. Please see the example in the examples/USER/mop folder.

[Restrictions:] 

This style is part of the USER-MOP package. It is only enabled if LAMMPS

is built with that package. See the "Build package"_Build_package.html

doc page on for more info.

The method is only implemented for 3d orthogonal simulation boxes whose

size does not change in time, and axis-aligned planes.

The method only works with two-body pair interactions, because it

requires the class method pair->single() to be implemented. In

particular, it does not work with more than two-body pair interactions,

intra-molecular interactions, and long range (kspace) interactions.

[Related commands:]

"compute stress/atom"_compute_stress_atom.html

[Default:] none

:line

:link(mop-todd)

[(Todd)] B. D. Todd, Denis J. Evans, and Peter J. Daivis: "Pressure tensor for inhomogeneous fluids", 

Phys. Rev. E 52, 1627 (1995).

variable T equal 0.8

variable p_solid equal 0.05

lattice fcc 1.0

region box block 0.0 6.0 0.0 6.0 -2.0 12.0

create_box 2 box

mass * 1.0

pair_style lj/cut 2.5

pair_coeff * * 1.0 1.0

pair_coeff 1 2 0.5 1.0 

pair_coeff 2 2 0.0 0.0

neigh_modify delay 0

region solid_bottom block INF INF INF INF -1.1 0.1

region liquid block INF INF INF INF 1.1 8.9

region solid_up block INF INF INF INF 9.9 11.1

create_atoms 1 region liquid

delete_atoms porosity liquid 0.26 88765

group liquid region liquid

create_atoms 2 region solid_bottom

group solid_bottom region solid_bottom

create_atoms 2 region solid_up

group solid_up region solid_up

group solid union solid_bottom solid_up

variable faSolid equal ${p_solid}*lx*ly/count(solid_up)

fix piston_up solid_up aveforce NULL NULL -${faSolid} 

fix freeze_up solid_up setforce 0.0 0.0 NULL 

fix freeze_bottom solid_bottom setforce 0.0 0.0 0.0

fix nvesol solid nve

compute Tliq liquid temp

fix nvtliq liquid nvt temp $T $T 0.5

fix_modify nvtliq temp Tliq

thermo 10000

thermo_modify flush yes temp Tliq

dump 1 all atom 10000 dump.lammpstrj

fix fxbal all balance 1000 1.05 shift z 10 1.05

velocity liquid create $T 47298 dist gaussian rot yes

run 50000

undump 1

reset_timestep 0

compute bin_z liquid chunk/atom bin/1d z 0.0 0.1 units box

compute liquidStress_ke liquid stress/atom NULL ke

compute liquidStress_vir liquid stress/atom NULL virial

fix profile_z liquid ave/chunk 10 1000 10000 bin_z density/number temp c_liquidStress_ke[1] c_liquidStress_ke[2] c_liquidStress_ke[3] c_liquidStress_ke[4] c_liquidStress_ke[5] c_liquidStress_ke[6] c_liquidStress_vir[1] c_liquidStress_vir[2] c_liquidStress_vir[3] c_liquidStress_vir[4] c_liquidStress_vir[5] c_liquidStress_vir[6] ave running overwrite file profile.z

compute mopz0 all stress/mop z center kin conf

fix mopz0t all ave/time 10 1000 10000 c_mopz0[*] file mopz0.time

compute moppz liquid stress/mop/profile z 0.0 0.1 kin conf

fix moppzt all ave/time 100 100 10000 c_moppz[*] ave running overwrite file moppz.time mode vector

run 40000

USER-MOP

============

This package provides the mop and mop/profile compute styles.

See the doc page for compute mop or compute mop/profile command for how to use

This package provides the stress/mop and stress/mop/profile compute styles.

See the doc page for compute stress/mop or compute stress/mop/profile command for how to use

them. 

PACKAGE DESCRIPTION

-------------------

This is a LAMMPS (http://lammps.sandia.gov/) compute style that calculates 

These are LAMMPS (http://lammps.sandia.gov/) compute styles that calculate 

components of the local stress tensor using the method of planes as described in

the paper by Todd et al. (B. D. Todd, Denis J. Evans, and Peter J. Daivis: "Pressure 

tensor for inhomogeneous fluids", Phys. Rev. E 52, 1627 (1995)). 

This package contains the source files of two compute styles:

* compute mop calculates components of the local stress tensor using the method of planes, applied on one plane.

* compute mop/profile calculates profiles of components of the local stress tensor using the method of planes, applied on an array of regularly spaced planes.

* compute stress/mop calculates components of the local stress tensor using the method of planes, applied on one plane.

* compute stress/mop/profile calculates profiles of components of the local stress tensor using the method of planes, applied on an array of regularly spaced planes.

The persons who created these files are Laurent Joly at University of Lyon 1 and Romain Vermorel at University of Pau and Pays de l'Adour.

Contact them directly if you have questions.

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

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.

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

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

  Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)

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

#include <mpi.h>

#include <cmath>

#include <cstring>

#include <cstdlib>

#include "compute_stress_mop.h"

#include "atom.h"

#include "update.h"

#include "domain.h"

#include "group.h"

#include "modify.h"

#include "fix.h"

#include "neighbor.h"

#include "force.h"

#include "pair.h"

#include "neigh_request.h"

#include "neigh_list.h"

#include "error.h"

#include "memory.h"

using namespace LAMMPS_NS;

enum{X,Y,Z};

enum{TOTAL,CONF,KIN};

#define BIG 1000000000

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

ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) :

  Compute(lmp, narg, arg)

{

  if (narg < 6) error->all(FLERR,"Illegal compute stress/mop command");

  MPI_Comm_rank(world,&me);

  // set compute mode and direction of plane(s) for pressure calculation

  if (strcmp(arg[3],"x")==0) {

    dir = X;

  } else if (strcmp(arg[3],"y")==0) {

    dir = Y;

  } else if (strcmp(arg[3],"z")==0) {

    dir = Z;

  } else error->all(FLERR,"Illegal compute stress/mop command");

  // Position of the plane

  if (strcmp(arg[4],"lower")==0) {

    pos = domain->boxlo[dir];

  } else if (strcmp(arg[4],"upper")==0) {

    pos = domain->boxhi[dir];

  } else if (strcmp(arg[4],"center")==0) {

    pos = 0.5*(domain->boxlo[dir]+domain->boxhi[dir]);

  } else pos = force->numeric(FLERR,arg[4]);

  if ( pos < (domain->boxlo[dir]+domain->prd_half[dir]) ) {

    pos1 = pos + domain->prd[dir];

  } else {

    pos1 = pos - domain->prd[dir];

  }

  // parse values until one isn't recognized

  which = new int[3*(narg-5)];

  nvalues = 0;

  int i;

  int iarg=5;

  while (iarg < narg) {

    if (strcmp(arg[iarg],"conf") == 0) {

      for (i=0; i<3; i++) {

        which[nvalues] = CONF;

        nvalues++;

      }

    } else if (strcmp(arg[iarg],"kin") == 0) {

      for (i=0; i<3; i++) {

        which[nvalues] = KIN;

        nvalues++;

      }

    } else if (strcmp(arg[iarg],"total") == 0) {

      for (i=0; i<3; i++) {

        which[nvalues] = TOTAL;

        nvalues++;

      }

    } else error->all(FLERR, "Illegal compute stress/mop command"); //break;

    iarg++;

  }

  // Error check

  // 3D only

  if (domain->dimension < 3)

    error->all(FLERR, "Compute stress/mop incompatible with simulation dimension");

  // orthogonal simulation box

  if (domain->triclinic != 0)

    error->all(FLERR, "Compute stress/mop incompatible with triclinic simulation box");

  // plane inside the box

  if (pos >domain->boxhi[dir] || pos <domain->boxlo[dir])

    error->all(FLERR, "Plane for compute stress/mop is out of bounds");

  // Initialize some variables

  values_local = values_global = vector = NULL;

  // this fix produces a global vector

  memory->create(vector,nvalues,"stress/mop:vector");

  memory->create(values_local,nvalues,"stress/mop/spatial:values_local");

  memory->create(values_global,nvalues,"stress/mop/spatial:values_global");

  size_vector = nvalues;

  vector_flag = 1;

  extvector = 0;

}

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

ComputeStressMop::~ComputeStressMop()

{

  delete [] which;

  memory->destroy(values_local);

  memory->destroy(values_global);

  memory->destroy(vector);

}

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

void ComputeStressMop::init()

{

  // Conversion constants

  nktv2p = force->nktv2p;

  ftm2v = force->ftm2v;

  // Plane area

  area = 1;

  int i;

  for (i=0; i<3; i++){

    if (i!=dir) area = area*domain->prd[i];

  }

  // Timestep Value

  dt = update->dt;

  // Error check

  // Compute stress/mop requires fixed simulation box

  if (domain->box_change_size || domain->box_change_shape || domain->deform_flag)

    error->all(FLERR, "Compute stress/mop requires a fixed simulation box");

  // This compute requires a pair style with pair_single method implemented

  if (force->pair == NULL)

    error->all(FLERR,"No pair style is defined for compute stress/mop");

  if (force->pair->single_enable == 0)

    error->all(FLERR,"Pair style does not support compute stress/mop");

  // Warnings

  if (me==0){

    //Compute stress/mop only accounts for pair interactions.

    // issue a warning if any intramolecular potential or Kspace is defined.

    if (force->bond!=NULL)

      error->warning(FLERR,"compute stress/mop does not account for bond potentials");

    if (force->angle!=NULL)

      error->warning(FLERR,"compute stress/mop does not account for angle potentials");

    if (force->dihedral!=NULL)

      error->warning(FLERR,"compute stress/mop does not account for dihedral potentials");

    if (force->improper!=NULL)

      error->warning(FLERR,"compute stress/mop does not account for improper potentials");

    if (force->kspace!=NULL)

      error->warning(FLERR,"compute stress/mop does not account for kspace contributions");

  }

  // need an occasional half neighbor list

  int irequest = neighbor->request((void *) this);

  neighbor->requests[irequest]->pair = 0;

  neighbor->requests[irequest]->compute = 1;

  neighbor->requests[irequest]->occasional = 1;

}

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

void ComputeStressMop::init_list(int id, NeighList *ptr)

{

  list = ptr;

}

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

   compute output vector

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

void ComputeStressMop::compute_vector()

{

  invoked_array = update->ntimestep;

  //Compute pressures on separate procs

  compute_pairs();

  // sum pressure contributions over all procs

  MPI_Allreduce(values_local,values_global,nvalues,

                MPI_DOUBLE,MPI_SUM,world);

  int m;

  for (m=0; m<nvalues; m++) {

    vector[m] = values_global[m];

  }

}

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

  compute pressure contribution of local proc

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

void ComputeStressMop::compute_pairs()

{

  int i,j,m,n,ii,jj,inum,jnum,itype,jtype;

  double delx,dely,delz;

  double rsq,eng,fpair,factor_coul,factor_lj;

  int *ilist,*jlist,*numneigh,**firstneigh;

  double *mass = atom->mass;

  int *type = atom->type;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  double *special_coul = force->special_coul;

  double *special_lj = force->special_lj;

  int newton_pair = force->newton_pair;

  // zero out arrays for one sample

  for (i = 0; i < nvalues; i++) values_local[i] = 0.0;

  // invoke half neighbor list (will copy or build if necessary)

  neighbor->build_one(list);

  inum = list->inum;

  ilist = list->ilist;

  numneigh = list->numneigh;

  firstneigh = list->firstneigh;

  // loop over neighbors of my atoms

  Pair *pair = force->pair;

  double **cutsq = force->pair->cutsq;

  // Parse values

  double xi[3];

  double vi[3];

  double fi[3];

  double xj[3];

  m = 0;

  while (m<nvalues) {

    if (which[m] == CONF || which[m] == TOTAL) {

      //Compute configurational contribution to pressure

      for (ii = 0; ii < inum; ii++) {

        i = ilist[ii];

        xi[0] = atom->x[i][0];

        xi[1] = atom->x[i][1];

        xi[2] = atom->x[i][2];

        itype = type[i];

        jlist = firstneigh[i];

        jnum = numneigh[i];

        for (jj = 0; jj < jnum; jj++) {

          j = jlist[jj];

          factor_lj = special_lj[sbmask(j)];

          factor_coul = special_coul[sbmask(j)];

          j &= NEIGHMASK;

          // skip if neither I nor J are in group

          if (!(mask[i] & groupbit || mask[j] & groupbit)) continue;

          xj[0] = atom->x[j][0];

          xj[1] = atom->x[j][1];

          xj[2] = atom->x[j][2];

          delx = xi[0] - xj[0];

          dely = xi[1] - xj[1];

          delz = xi[2] - xj[2];

          rsq = delx*delx + dely*dely + delz*delz;

          jtype = type[j];

          if (rsq >= cutsq[itype][jtype]) continue;

          if (newton_pair || j < nlocal) {

            //check if ij pair is accross plane, add contribution to pressure

            if ( ((xi[dir]>pos) && (xj[dir]<pos)) || ((xi[dir]>pos1) && (xj[dir]<pos1)) ) {

              pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);

              values_local[m] += fpair*(xi[0]-xj[0])/area*nktv2p;

              values_local[m+1] += fpair*(xi[1]-xj[1])/area*nktv2p;

              values_local[m+2] += fpair*(xi[2]-xj[2])/area*nktv2p;

            }

            else if ( ((xi[dir]<pos) && (xj[dir]>pos)) || ((xi[dir]<pos1) && (xj[dir]>pos1)) ){

              pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);

              values_local[m] -= fpair*(xi[0]-xj[0])/area*nktv2p;

              values_local[m+1] -= fpair*(xi[1]-xj[1])/area*nktv2p;

              values_local[m+2] -= fpair*(xi[2]-xj[2])/area*nktv2p;

            }

          } else {

            if ( ((xi[dir]>pos) && (xj[dir]<pos)) || ((xi[dir]>pos1) && (xj[dir]<pos1)) ) {

              pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);

              values_local[m] += fpair*(xi[0]-xj[0])/area*nktv2p;

              values_local[m+1] += fpair*(xi[1]-xj[1])/area*nktv2p;

              values_local[m+2] += fpair*(xi[2]-xj[2])/area*nktv2p;

            }

          }

        }

      }

    }

    // Compute kinetic contribution to pressure

    // counts local particles transfers across the plane

    if (which[m] == KIN || which[m] == TOTAL){

      double vcm[3];

      double masstotal,sgn;

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

        // skip if I is not in group

        if (mask[i] & groupbit){

          itype = type[i];

          //coordinates at t

          xi[0] = atom->x[i][0];

          xi[1] = atom->x[i][1];

          xi[2] = atom->x[i][2];

          //velocities at t

          vi[0] = atom->v[i][0];

          vi[1] = atom->v[i][1];

          vi[2] = atom->v[i][2];

          //forces at t

          fi[0] = atom->f[i][0];

          fi[1] = atom->f[i][1];

          fi[2] = atom->f[i][2];

          //coordinates at t-dt (based on Velocity-Verlet alg.)

          xj[0] = xi[0]-vi[0]*dt+fi[0]/2/mass[itype]*dt*dt*ftm2v;

          xj[1] = xi[1]-vi[1]*dt+fi[1]/2/mass[itype]*dt*dt*ftm2v;

          xj[2] = xi[2]-vi[2]*dt+fi[2]/2/mass[itype]*dt*dt*ftm2v;

          // because LAMMPS does not put atoms back in the box

          // at each timestep, must check atoms going through the

          // image of the plane that is closest to the box

          double pos_temp = pos+copysign(1,domain->prd_half[dir]-pos)*domain->prd[dir];

          if (fabs(xi[dir]-pos)<fabs(xi[dir]-pos_temp)) pos_temp = pos;

          if (((xi[dir]-pos_temp)*(xj[dir]-pos_temp)<0)){

            //sgn = copysign(1,vi[dir]-vcm[dir]);

            sgn = copysign(1,vi[dir]);

            //approximate crossing velocity by v(t-dt/2) (based on Velocity-Verlet alg.)

            double vcross[3];

            vcross[0] = vi[0]-fi[0]/mass[itype]/2*ftm2v*dt;

            vcross[1] = vi[1]-fi[1]/mass[itype]/2*ftm2v*dt;

            vcross[2] = vi[2]-fi[2]/mass[itype]/2*ftm2v*dt;

            values_local[m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;

            values_local[m+1] += mass[itype]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;

            values_local[m+2] += mass[itype]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;

          }

        }

      }

    }

    m+=3;

  }

}

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

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.

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

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

  Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon)

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

#ifdef COMPUTE_CLASS

ComputeStyle(stress/mop,ComputeStressMop)

#else

#ifndef LMP_COMPUTE_STRESS_MOP_H

#define LMP_COMPUTE_STRESS_MOP_H

#include "compute.h"

namespace LAMMPS_NS {

  class ComputeStressMop : public Compute {

  public:

    ComputeStressMop(class LAMMPS *, int, char **);

    virtual ~ComputeStressMop();

    void init();

    void init_list(int, class NeighList *);

    void compute_vector();

  private:

    void compute_pairs();

    int me,nvalues,dir;

    int *which;

    double *values_local,*values_global;

    double pos,pos1,dt,nktv2p,ftm2v;

    double area;

    class NeighList *list;

  };

}

#endif

#endif

/* ERROR/WARNING messages:

 E: Illegal ... command

 Self-explanatory.  Check the input script syntax and compare to the

 documentation for the command.  You can use -echo screen as a

 command-line option when running LAMMPS to see the offending line.

 E: Compute stress/mop incompatible with simulation dimension

 Compute stress/mop only works with 3D simulations.

 E: Compute stress/mop incompatible with triclinic simulation box

 Self-explanatory.

 E: Compute stress/mop requires a fixed simulation box

 Compute stress/mop is not compatible with any change of volume or shape

 or boundary conditions of the simulation box.

 E: No pair style is defined for compute stress/mop

 Self-explanatory. Compute stress/mop requires the definition of a pair style.

 E: Pair style does not support compute stress/mop

 The pair style does not have a single() function, so it can

 not be invoked by compute stress/mop.

 W: compute stress/mop does not account for bond potentials

 W: compute stress/mop does not account for angle potentials

 W: compute stress/mop does not account for dihedral potentials

 W: compute stress/mop does not account for improper potentials

 W: compute stress/mop does not account for kspace contributions

 Compute stress/mop only accounts for pairwise additive interactions for

 the computation of local stress tensor components.

 */