new USER-MOP package submitted

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

them. 

PACKAGE DESCRIPTION

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

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

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.

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.

Laurent Joly

Institut Lumière Matière

Université Lyon 1 et CNRS

43 boulevard du 11 novembre 1918

69622 Villeurbanne Cedex

e-mail: laurent.joly@univ-lyon1.fr

Romain Vermorel

Laboratory of Complex Fluids and their Reservoirs (LFCR)

Université de Pau et des Pays de l'Adour

Avenue de l'Université 

64012 Pau

e-mail: romain.vermorel@univ-pau.fr

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

   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 "math.h"

#include "mpi.h"

#include "string.h"

#include "stdlib.h"

#include "compute_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"

#include <iostream>

using namespace LAMMPS_NS;

enum{X,Y,Z};

enum{TOTAL,CONF,KIN};

#define BIG 1000000000

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

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

  Compute(lmp, narg, arg)

{

  if (narg < 6) error->all(FLERR,"Illegal compute 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 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 mop command"); //break;

    iarg++;

  }

  //Error check

    // 3D only

    if (domain->dimension<3)

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

    // orthogonal simulation box

    if (domain->triclinic != 0)

        error->all(FLERR, "Compute 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 mop is out of bounds");

  // Initialize some variables

  values_local = values_global = vector = NULL;

  // this fix produces a global vector

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

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

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

  size_vector = nvalues;

  vector_flag = 1;

  extvector = 0;

}

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

ComputeMop::~ComputeMop()

{

  delete [] which;

  memory->destroy(values_local);

  memory->destroy(values_global);

  memory->destroy(vector);

}

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

void ComputeMop::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 mop requires fixed simulation box

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

    error->all(FLERR, "Compute 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 mop");

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

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

// Warnings

if (me==0){

  //Compute 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 mop does not account for bond potentials");

  if (force->angle!=NULL)

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

  if (force->dihedral!=NULL)

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

  if (force->improper!=NULL)

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

  if (force->kspace!=NULL)

        error->warning(FLERR,"compute 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 ComputeMop::init_list(int id, NeighList *ptr)

{

  list = ptr;

}

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

   compute output vector

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

void ComputeMop::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 ComputeMop::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;

            }

          }

        }

      }

    }

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

      //Compute kinetic contribution to pressure

      // counts local particles transfers across the plane

      double vcm[3];

      double masstotal,sgn;

      //Velocity of the center of mass

      //int ivcm,vcmbit;

      //ivcm = igroup;

      //vcmbit = groupbit;

      //masstotal = group->mass(ivcm);

      //group->vcm(ivcm,masstotal,vcm);

      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;

  }

}

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 mop z center kin conf

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

compute moppz liquid 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

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

   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(mop,ComputeMop)

#else

#ifndef LMP_COMPUTE_MOP_H

#define LMP_COMPUTE_MOP_H

#include "compute.h"

namespace LAMMPS_NS {

  class ComputeMop : public Compute {

  public:

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

    virtual ~ComputeMop();

    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