jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10366 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2013-07-25 22:26:18 +00:00
parent dea3df948e
commit 87b941f250
12 changed files with 0 additions and 3214 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

508
src/fix_deposit.cpp
View File

@ -1,508 +0,0 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "fix_deposit.h"
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "comm.h"
#include "domain.h"
#include "lattice.h"
#include "region.h"
#include "random_park.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) error->all(FLERR,"Illegal fix deposit command");
restart_global = 1;
time_depend = 1;
// required args
ninsert = force->inumeric(FLERR,arg[3]);
ntype = force->inumeric(FLERR,arg[4]);
nfreq = force->inumeric(FLERR,arg[5]);
seed = force->inumeric(FLERR,arg[6]);
if (seed <= 0) error->all(FLERR,"Illegal fix deposit command");
// set defaults
iregion = -1;
idregion = NULL;
idnext = 0;
globalflag = localflag = 0;
lo = hi = deltasq = 0.0;
nearsq = 0.0;
maxattempt = 10;
rateflag = 0;
vxlo = vxhi = vylo = vyhi = vzlo = vzhi = 0.0;
scaleflag = 1;
targetflag = 0;
// read options from end of input line
options(narg-7,&arg[7]);
// error checks on region and its extent being inside simulation box
if (iregion == -1) error->all(FLERR,"Must specify a region in fix deposit");
if (domain->regions[iregion]->bboxflag == 0)
error->all(FLERR,"Fix deposit region does not support a bounding box");
if (domain->regions[iregion]->dynamic_check())
error->all(FLERR,"Fix deposit region cannot be dynamic");
xlo = domain->regions[iregion]->extent_xlo;
xhi = domain->regions[iregion]->extent_xhi;
ylo = domain->regions[iregion]->extent_ylo;
yhi = domain->regions[iregion]->extent_yhi;
zlo = domain->regions[iregion]->extent_zlo;
zhi = domain->regions[iregion]->extent_zhi;
if (domain->triclinic == 0) {
if (xlo < domain->boxlo[0] || xhi > domain->boxhi[0] ||
ylo < domain->boxlo[1] || yhi > domain->boxhi[1] ||
zlo < domain->boxlo[2] || zhi > domain->boxhi[2])
error->all(FLERR,"Deposition region extends outside simulation box");
} else {
if (xlo < domain->boxlo_bound[0] || xhi > domain->boxhi_bound[0] ||
ylo < domain->boxlo_bound[1] || yhi > domain->boxhi_bound[1] ||
zlo < domain->boxlo_bound[2] || zhi > domain->boxhi_bound[2])
error->all(FLERR,"Deposition region extends outside simulation box");
}
// setup scaling
double xscale,yscale,zscale;
if (scaleflag) {
xscale = domain->lattice->xlattice;
yscale = domain->lattice->ylattice;
zscale = domain->lattice->zlattice;
}
else xscale = yscale = zscale = 1.0;
// apply scaling to all input parameters with dist/vel units
if (domain->dimension == 2) {
lo *= yscale;
hi *= yscale;
rate *= yscale;
} else {
lo *= zscale;
hi *= zscale;
rate *= zscale;
}
deltasq *= xscale*xscale;
nearsq *= xscale*xscale;
vxlo *= xscale;
vxhi *= xscale;
vylo *= yscale;
vyhi *= yscale;
vzlo *= zscale;
vzhi *= zscale;
tx *= xscale;
ty *= yscale;
tz *= zscale;
// maxtag_all = current max tag for all atoms
if (idnext) {
int *tag = atom->tag;
int nlocal = atom->nlocal;
int maxtag = 0;
for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]);
MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_INT,MPI_MAX,world);
}
// random number generator, same for all procs
random = new RanPark(lmp,seed);
// set up reneighboring
force_reneighbor = 1;
next_reneighbor = update->ntimestep + 1;
nfirst = next_reneighbor;
ninserted = 0;
}
/* ---------------------------------------------------------------------- */
FixDeposit::~FixDeposit()
{
delete random;
delete [] idregion;
}
/* ---------------------------------------------------------------------- */
int FixDeposit::setmask()
{
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixDeposit::init()
{
// set index and check validity of region
iregion = domain->find_region(idregion);
if (iregion == -1)
error->all(FLERR,"Region ID for fix deposit does not exist");
}
/* ----------------------------------------------------------------------
perform particle insertion
------------------------------------------------------------------------- */
void FixDeposit::pre_exchange()
{
int i,j;
int flag,flagall;
double coord[3],lamda[3],delx,dely,delz,rsq;
double *newcoord;
// just return if should not be called on this timestep
if (next_reneighbor != update->ntimestep) return;
// compute current offset = bottom of insertion volume
double offset = 0.0;
if (rateflag) offset = (update->ntimestep - nfirst) * update->dt * rate;
double *sublo,*subhi;
if (domain->triclinic == 0) {
sublo = domain->sublo;
subhi = domain->subhi;
} else {
sublo = domain->sublo_lamda;
subhi = domain->subhi_lamda;
}
// attempt an insertion until successful
int nfix = modify->nfix;
Fix **fix = modify->fix;
int success = 0;
int attempt = 0;
while (attempt < maxattempt) {
attempt++;
// choose random position for new atom within region
coord[0] = xlo + random->uniform() * (xhi-xlo);
coord[1] = ylo + random->uniform() * (yhi-ylo);
coord[2] = zlo + random->uniform() * (zhi-zlo);
while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) {
coord[0] = xlo + random->uniform() * (xhi-xlo);
coord[1] = ylo + random->uniform() * (yhi-ylo);
coord[2] = zlo + random->uniform() * (zhi-zlo);
}
// adjust vertical coord by offset
if (domain->dimension == 2) coord[1] += offset;
else coord[2] += offset;
// if global, reset vertical coord to be lo-hi above highest atom
// if local, reset vertical coord to be lo-hi above highest "nearby" atom
// local computation computes lateral distance between 2 particles w/ PBC
if (globalflag || localflag) {
int dim;
double max,maxall,delx,dely,delz,rsq;
if (domain->dimension == 2) {
dim = 1;
max = domain->boxlo[1];
} else {
dim = 2;
max = domain->boxlo[2];
}
double **x = atom->x;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (localflag) {
delx = coord[0] - x[i][0];
dely = coord[1] - x[i][1];
delz = 0.0;
domain->minimum_image(delx,dely,delz);
if (domain->dimension == 2) rsq = delx*delx;
else rsq = delx*delx + dely*dely;
if (rsq > deltasq) continue;
}
if (x[i][dim] > max) max = x[i][dim];
}
MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world);
if (domain->dimension == 2)
coord[1] = maxall + lo + random->uniform()*(hi-lo);
else
coord[2] = maxall + lo + random->uniform()*(hi-lo);
}
// now have final coord
// if distance to any atom is less than near, try again
double **x = atom->x;
int nlocal = atom->nlocal;
flag = 0;
for (i = 0; i < nlocal; i++) {
delx = coord[0] - x[i][0];
dely = coord[1] - x[i][1];
delz = coord[2] - x[i][2];
domain->minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < nearsq) flag = 1;
}
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
if (flagall) continue;
// insertion will proceed
// choose random velocity for new atom
double vxtmp = vxlo + random->uniform() * (vxhi-vxlo);
double vytmp = vylo + random->uniform() * (vyhi-vylo);
double vztmp = vzlo + random->uniform() * (vzhi-vzlo);
// if we have a sputter target change velocity vector accordingly
if (targetflag) {
double vel = sqrt(vxtmp*vxtmp + vytmp*vytmp + vztmp*vztmp);
delx = tx - coord[0];
dely = ty - coord[1];
delz = tz - coord[2];
double rsq = delx*delx + dely*dely + delz*delz;
if (rsq > 0.0) {
double rinv = sqrt(1.0/rsq);
vxtmp = delx*rinv*vel;
vytmp = dely*rinv*vel;
vztmp = delz*rinv*vel;
}
}
// check if new atom is in my sub-box or above it if I'm highest proc
// if so, add to my list via create_atom()
// initialize info about the atoms
// set group mask to "all" plus fix group
if (domain->triclinic) {
domain->x2lamda(coord,lamda);
newcoord = lamda;
} else newcoord = coord;
flag = 0;
if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] &&
newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1;
else if (domain->dimension == 3 && newcoord[2] >= domain->boxhi[2] &&
comm->myloc[2] == comm->procgrid[2]-1 &&
newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1;
else if (domain->dimension == 2 && newcoord[1] >= domain->boxhi[1] &&
comm->myloc[1] == comm->procgrid[1]-1 &&
newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1;
if (flag) {
atom->avec->create_atom(ntype,coord);
int m = atom->nlocal - 1;
atom->type[m] = ntype;
atom->mask[m] = 1 | groupbit;
atom->v[m][0] = vxtmp;
atom->v[m][1] = vytmp;
atom->v[m][2] = vztmp;
for (j = 0; j < nfix; j++)
if (fix[j]->create_attribute) fix[j]->set_arrays(m);
}
MPI_Allreduce(&flag,&success,1,MPI_INT,MPI_MAX,world);
break;
}
// warn if not successful b/c too many attempts or no proc owned particle
if (!success && comm->me == 0)
error->warning(FLERR,"Particle deposition was unsuccessful",0);
// reset global natoms
// if idnext, set new atom ID to incremented maxtag_all
// else set new atom ID to value beyond all current atoms
// if global map exists, reset it now instead of waiting for comm
// since adding an atom messes up ghosts
if (success) {
atom->natoms += 1;
if (atom->tag_enable) {
if (idnext) {
maxtag_all++;
if (atom->nlocal && atom->tag[atom->nlocal-1] == 0)
atom->tag[atom->nlocal-1] = maxtag_all;
} else atom->tag_extend();
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
}
}
// next timestep to insert
// next_reneighbor = 0 if done
if (success) ninserted++;
if (ninserted < ninsert) next_reneighbor += nfreq;
else next_reneighbor = 0;
}
/* ----------------------------------------------------------------------
parse optional parameters at end of input line
------------------------------------------------------------------------- */
void FixDeposit::options(int narg, char **arg)
{
if (narg < 0) error->all(FLERR,"Illegal fix indent command");
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"region") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1)
error->all(FLERR,"Region ID for fix deposit does not exist");
int n = strlen(arg[iarg+1]) + 1;
idregion = new char[n];
strcpy(idregion,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"id") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
if (strcmp(arg[iarg+1],"max") == 0) idnext = 0;
else if (strcmp(arg[iarg+1],"next") == 0) idnext = 1;
else error->all(FLERR,"Illegal fix deposit command");
iarg += 2;
} else if (strcmp(arg[iarg],"global") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
globalflag = 1;
localflag = 0;
lo = force->numeric(FLERR,arg[iarg+1]);
hi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"local") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command");
localflag = 1;
globalflag = 0;
lo = force->numeric(FLERR,arg[iarg+1]);
hi = force->numeric(FLERR,arg[iarg+2]);
deltasq = force->numeric(FLERR,arg[iarg+3])*force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"near") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
nearsq = force->numeric(FLERR,arg[iarg+1])*force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"attempt") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
maxattempt = force->inumeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"rate") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
rateflag = 1;
rate = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"vx") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
vxlo = force->numeric(FLERR,arg[iarg+1]);
vxhi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"vy") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
vylo = force->numeric(FLERR,arg[iarg+1]);
vyhi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"vz") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
vzlo = force->numeric(FLERR,arg[iarg+1]);
vzhi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"units") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
else error->all(FLERR,"Illegal fix deposit command");
iarg += 2;
} else if (strcmp(arg[iarg],"target") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command");
tx = force->numeric(FLERR,arg[iarg+1]);
ty = force->numeric(FLERR,arg[iarg+2]);
tz = force->numeric(FLERR,arg[iarg+3]);
targetflag = 1;
iarg += 4;
} else error->all(FLERR,"Illegal fix deposit command");
}
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixDeposit::write_restart(FILE *fp)
{
int n = 0;
double list[4];
list[n++] = random->state();
list[n++] = ninserted;
list[n++] = nfirst;
list[n++] = next_reneighbor;
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(list,sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixDeposit::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
seed = static_cast<int> (list[n++]);
ninserted = static_cast<int> (list[n++]);
nfirst = static_cast<int> (list[n++]);
next_reneighbor = static_cast<int> (list[n++]);
random->reset(seed);
}

98
src/fix_deposit.h
View File

@ -1,98 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(deposit,FixDeposit)
#else
#ifndef LMP_FIX_DEPOSIT_H
#define LMP_FIX_DEPOSIT_H
#include "stdio.h"
#include "fix.h"
namespace LAMMPS_NS {
class FixDeposit : public Fix {
public:
FixDeposit(class LAMMPS *, int, char **);
~FixDeposit();
int setmask();
void init();
void pre_exchange();
void write_restart(FILE *);
void restart(char *);
private:
int ninsert,ntype,nfreq,seed;
int iregion,globalflag,localflag,maxattempt,rateflag,scaleflag,targetflag;
char *idregion;
double lo,hi,deltasq,nearsq,rate;
double vxlo,vxhi,vylo,vyhi,vzlo,vzhi;
double xlo,xhi,ylo,yhi,zlo,zhi;
double tx,ty,tz;
int nfirst,ninserted;
int idnext,maxtag_all;
class RanPark *random;
void options(int, char **);
};
}
#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: Must specify a region in fix deposit
The region keyword must be specified with this fix.
E: Fix deposit region does not support a bounding box
Not all regions represent bounded volumes. You cannot use
such a region with the fix deposit command.
E: Fix deposit region cannot be dynamic
Only static regions can be used with fix deposit.
E: Deposition region extends outside simulation box
Self-explanatory.
E: Region ID for fix deposit does not exist
Self-explanatory.
W: Particle deposition was unsuccessful
The fix deposit command was not able to insert as many atoms as
needed. The requested volume fraction may be too high, or other atoms
may be in the insertion region.
U: Use of fix deposit with undefined lattice
Must use lattice command with compute fix deposit command if units
option is set to lattice.
*/

406
src/fix_efield.cpp
View File

@ -1,406 +0,0 @@
/* ----------------------------------------------------------------------
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 author: Christina Payne (Vanderbilt U)
Stan Moore (Sandia) for dipole terms
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "fix_efield.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "comm.h"
#include "modify.h"
#include "force.h"
#include "respa.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{NONE,CONSTANT,EQUAL,ATOM};
/* ---------------------------------------------------------------------- */
FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg != 6) error->all(FLERR,"Illegal fix efield command");
vector_flag = 1;
scalar_flag = 1;
size_vector = 3;
global_freq = 1;
extvector = 1;
extscalar = 1;
qe2f = force->qe2f;
xstr = ystr = zstr = NULL;
if (strstr(arg[3],"v_") == arg[3]) {
int n = strlen(&arg[3][2]) + 1;
xstr = new char[n];
strcpy(xstr,&arg[3][2]);
} else {
ex = qe2f * force->numeric(FLERR,arg[3]);
xstyle = CONSTANT;
}
if (strstr(arg[4],"v_") == arg[4]) {
int n = strlen(&arg[4][2]) + 1;
ystr = new char[n];
strcpy(ystr,&arg[4][2]);
} else {
ey = qe2f * force->numeric(FLERR,arg[4]);
ystyle = CONSTANT;
}
if (strstr(arg[5],"v_") == arg[5]) {
int n = strlen(&arg[5][2]) + 1;
zstr = new char[n];
strcpy(zstr,&arg[5][2]);
} else {
ez = qe2f * force->numeric(FLERR,arg[5]);
zstyle = CONSTANT;
}
// optional args
estr = NULL;
int iarg = 6;
while (iarg < narg) {
if (strcmp(arg[iarg],"energy") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix efield command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
int n = strlen(&arg[iarg+1][2]) + 1;
estr = new char[n];
strcpy(estr,&arg[iarg+1][2]);
} else error->all(FLERR,"Illegal fix efield command");
iarg += 2;
} else error->all(FLERR,"Illegal fix efield command");
}
force_flag = 0;
fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0;
maxatom = 0;
efield = NULL;
}
/* ---------------------------------------------------------------------- */
FixEfield::~FixEfield()
{
delete [] xstr;
delete [] ystr;
delete [] zstr;
delete [] estr;
memory->destroy(efield);
}
/* ---------------------------------------------------------------------- */
int FixEfield::setmask()
{
int mask = 0;
mask |= THERMO_ENERGY;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixEfield::init()
{
qflag = muflag = 0;
if (atom->q_flag) qflag = 1;
if (atom->mu_flag && atom->torque_flag) muflag = 1;
if (!qflag && !muflag)
error->all(FLERR,"Fix efield requires atom attribute q or mu");
// check variables
if (xstr) {
xvar = input->variable->find(xstr);
if (xvar < 0)
error->all(FLERR,"Variable name for fix efield does not exist");
if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
else if (input->variable->atomstyle(xvar)) xstyle = ATOM;
else error->all(FLERR,"Variable for fix efield is invalid style");
}
if (ystr) {
yvar = input->variable->find(ystr);
if (yvar < 0)
error->all(FLERR,"Variable name for fix efield does not exist");
if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
else if (input->variable->atomstyle(yvar)) ystyle = ATOM;
else error->all(FLERR,"Variable for fix efield is invalid style");
}
if (zstr) {
zvar = input->variable->find(zstr);
if (zvar < 0)
error->all(FLERR,"Variable name for fix efield does not exist");
if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
else if (input->variable->atomstyle(zvar)) zstyle = ATOM;
else error->all(FLERR,"Variable for fix efield is invalid style");
}
if (estr) {
evar = input->variable->find(estr);
if (evar < 0)
error->all(FLERR,"Variable name for fix efield does not exist");
if (input->variable->atomstyle(evar)) estyle = ATOM;
else error->all(FLERR,"Variable for fix efield is invalid style");
} else estyle = NONE;
if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)
varflag = ATOM;
else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL)
varflag = EQUAL;
else varflag = CONSTANT;
if (muflag && varflag == ATOM)
error->all(FLERR,"Fix efield with dipoles cannot use atom-style variables");
if (muflag && update->whichflag == 2 && comm->me == 0)
error->warning(FLERR,
"The minimizer does not re-orient dipoles "
"when using fix efield");
if (varflag == CONSTANT && estyle != NONE)
error->all(FLERR,"Cannot use variable energy with "
"constant efield in fix efield");
if ((varflag == EQUAL || varflag == ATOM) &&
update->whichflag == 2 && estyle == NONE)
error->all(FLERR,"Must use variable energy with fix efield");
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixEfield::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}
/* ---------------------------------------------------------------------- */
void FixEfield::min_setup(int vflag)
{
post_force(vflag);
}
/* ----------------------------------------------------------------------
apply F = qE
------------------------------------------------------------------------- */
void FixEfield::post_force(int vflag)
{
double **f = atom->f;
double *q = atom->q;
int *mask = atom->mask;
tagint *image = atom->image;
int nlocal = atom->nlocal;
// reallocate efield array if necessary
if (varflag == ATOM && nlocal > maxatom) {
maxatom = atom->nmax;
memory->destroy(efield);
memory->create(efield,maxatom,4,"efield:efield");
}
// fsum[0] = "potential energy" for added force
// fsum[123] = extra force added to atoms
fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0;
force_flag = 0;
double **x = atom->x;
double fx,fy,fz;
// constant efield
if (varflag == CONSTANT) {
double unwrap[3];
// charge interactions
// force = qE, potential energy = F dot x in unwrapped coords
if (qflag) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
fx = q[i]*ex;
fy = q[i]*ey;
fz = q[i]*ez;
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
domain->unmap(x[i],image[i],unwrap);
fsum[0] -= fx*unwrap[0]+fy*unwrap[1]+fz*unwrap[2];
fsum[1] += fx;
fsum[2] += fy;
fsum[3] += fz;
}
}
// dipole interactions
// no force, torque = mu cross E, potential energy = -mu dot E
if (muflag) {
double **mu = atom->mu;
double **t = atom->torque;
double tx,ty,tz;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
tx = ez*mu[i][1] - ey*mu[i][2];
ty = ex*mu[i][2] - ez*mu[i][0];
tz = ey*mu[i][0] - ex*mu[i][1];
t[i][0] += tx;
t[i][1] += ty;
t[i][2] += tz;
fsum[0] -= mu[i][0]*ex + mu[i][1]*ey + mu[i][2]*ez;
}
}
// variable efield, wrap with clear/add
// potential energy = evar if defined, else 0.0
} else {
modify->clearstep_compute();
if (xstyle == EQUAL) ex = qe2f * input->variable->compute_equal(xvar);
else if (xstyle == ATOM && efield)
input->variable->compute_atom(xvar,igroup,&efield[0][0],3,0);
if (ystyle == EQUAL) ey = qe2f * input->variable->compute_equal(yvar);
else if (ystyle == ATOM && efield)
input->variable->compute_atom(yvar,igroup,&efield[0][1],3,0);
if (zstyle == EQUAL) ez = qe2f * input->variable->compute_equal(zvar);
else if (zstyle == ATOM && efield)
input->variable->compute_atom(zvar,igroup,&efield[0][2],3,0);
if (estyle == ATOM && efield)
input->variable->compute_atom(evar,igroup,&efield[0][3],4,0);
modify->addstep_compute(update->ntimestep + 1);
// charge interactions
// force = qE
if (qflag) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (xstyle == ATOM) fx = qe2f * q[i]*efield[i][0];
else fx = q[i]*ex;
f[i][0] += fx;
fsum[1] += fx;
if (ystyle == ATOM) fy = qe2f * q[i]*efield[i][1];
else fy = q[i]*ey;
f[i][1] += fy;
fsum[2] += fy;
if (zstyle == ATOM) fz = qe2f * q[i]*efield[i][2];
else fz = q[i]*ez;
f[i][2] += fz;
fsum[3] += fz;
if (estyle == ATOM) fsum[0] += efield[0][3];
}
}
// dipole interactions
// no force, torque = mu cross E
if (muflag) {
double **mu = atom->mu;
double **t = atom->torque;
double tx,ty,tz;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
tx = ez*mu[i][1] - ey*mu[i][2];
ty = ex*mu[i][2] - ez*mu[i][0];
tz = ey*mu[i][0] - ex*mu[i][1];
t[i][0] += tx;
t[i][1] += ty;
t[i][2] += tz;
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixEfield::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixEfield::min_post_force(int vflag)
{
post_force(vflag);
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double FixEfield::memory_usage()
{
double bytes = 0.0;
if (varflag == ATOM) bytes = atom->nmax*4 * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
return energy added by fix
------------------------------------------------------------------------- */
double FixEfield::compute_scalar(void)
{
if (force_flag == 0) {
MPI_Allreduce(fsum,fsum_all,4,MPI_DOUBLE,MPI_SUM,world);
force_flag = 1;
}
return fsum_all[0];
}
/* ----------------------------------------------------------------------
return total extra force due to fix
------------------------------------------------------------------------- */
double FixEfield::compute_vector(int n)
{
if (force_flag == 0) {
MPI_Allreduce(fsum,fsum_all,4,MPI_DOUBLE,MPI_SUM,world);
force_flag = 1;
}
return fsum_all[n+1];
}

101
src/fix_efield.h
View File

@ -1,101 +0,0 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(efield,FixEfield)
#else
#ifndef LMP_FIX_EFIELD_H
#define LMP_FIX_EFIELD_H
#include "fix.h"
namespace LAMMPS_NS {
class FixEfield : public Fix {
public:
FixEfield(class LAMMPS *, int, char **);
~FixEfield();
int setmask();
void init();
void setup(int);
void min_setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
double memory_usage();
double compute_scalar();
double compute_vector(int);
private:
double ex,ey,ez;
int varflag;
char *xstr,*ystr,*zstr,*estr;
int xvar,yvar,zvar,evar,xstyle,ystyle,zstyle,estyle;
int nlevels_respa;
double qe2f;
double fdotx;
int qflag,muflag;
int maxatom;
double **efield;
int force_flag;
double fsum[4],fsum_all[4];
};
}
#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: Fix efield requires atom attribute q or mu
Self-explanatory.
E: Variable name for fix efield does not exist
Self-explanatory.
E: Variable for fix efield is invalid style
Only equal-style or atom-style variables can be used.
E: Fix efield with dipoles cannot use atom-style variables
This feature is not yet supported.
W: The minimizer does not re-orient dipoles when using fix efield
Self-explanatory.
E: Cannot use variable energy with constant efield in fix efield
Self-explanatory.
E: Must use variable energy with fix efield
One or more variables are defined for fix efield, which require
variable energy when using the minimizer.
*/

386
src/fix_evaporate.cpp
View File

@ -1,386 +0,0 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "fix_evaporate.h"
#include "atom.h"
#include "atom_vec.h"
#include "update.h"
#include "domain.h"
#include "region.h"
#include "comm.h"
#include "force.h"
#include "group.h"
#include "random_park.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixEvaporate::FixEvaporate(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) error->all(FLERR,"Illegal fix evaporate command");
scalar_flag = 1;
global_freq = 1;
extscalar = 0;
nevery = force->inumeric(FLERR,arg[3]);
nflux = force->inumeric(FLERR,arg[4]);
iregion = domain->find_region(arg[5]);
int n = strlen(arg[5]) + 1;
idregion = new char[n];
strcpy(idregion,arg[5]);
int seed = force->inumeric(FLERR,arg[6]);
if (nevery <= 0 || nflux <= 0)
error->all(FLERR,"Illegal fix evaporate command");
if (iregion == -1)
error->all(FLERR,"Region ID for fix evaporate does not exist");
if (seed <= 0) error->all(FLERR,"Illegal fix evaporate command");
// random number generator, same for all procs
random = new RanPark(lmp,seed);
// optional args
molflag = 0;
int iarg = 7;
while (iarg < narg) {
if (strcmp(arg[iarg],"molecule") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix evaporate command");
if (strcmp(arg[iarg+1],"no") == 0) molflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) molflag = 1;
else error->all(FLERR,"Illegal fix evaporate command");
iarg += 2;
} else error->all(FLERR,"Illegal fix evaporate command");
}
// set up reneighboring
force_reneighbor = 1;
next_reneighbor = (update->ntimestep/nevery)*nevery + nevery;
ndeleted = 0;
nmax = 0;
list = NULL;
mark = NULL;
}
/* ---------------------------------------------------------------------- */
FixEvaporate::~FixEvaporate()
{
delete [] idregion;
delete random;
memory->destroy(list);
memory->destroy(mark);
}
/* ---------------------------------------------------------------------- */
int FixEvaporate::setmask()
{
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixEvaporate::init()
{
// set index and check validity of region
iregion = domain->find_region(idregion);
if (iregion == -1)
error->all(FLERR,"Region ID for fix evaporate does not exist");
// check that no deletable atoms are in atom->firstgroup
// deleting such an atom would not leave firstgroup atoms first
if (atom->firstgroup >= 0) {
int *mask = atom->mask;
int nlocal = atom->nlocal;
int firstgroupbit = group->bitmask[atom->firstgroup];
int flag = 0;
for (int i = 0; i < nlocal; i++)
if ((mask[i] & groupbit) && (mask[i] && firstgroupbit)) flag = 1;
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall)
error->all(FLERR,"Cannot evaporate atoms in atom_modify first group");
}
// if molflag not set, warn if any deletable atom has a mol ID
if (molflag == 0 && atom->molecule_flag) {
int *molecule = atom->molecule;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int flag = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (molecule[i]) flag = 1;
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall && comm->me == 0)
error->warning(FLERR,
"Fix evaporate may delete atom with non-zero molecule ID");
}
if (molflag && atom->molecule_flag == 0)
error->all(FLERR,
"Fix evaporate molecule requires atom attribute molecule");
}
/* ----------------------------------------------------------------------
perform particle deletion
done before exchange, borders, reneighbor
so that ghost atoms and neighbor lists will be correct
------------------------------------------------------------------------- */
void FixEvaporate::pre_exchange()
{
int i,j,m,iwhichglobal,iwhichlocal;
int ndel,ndeltopo[4];
if (update->ntimestep != next_reneighbor) return;
// grow list and mark arrays if necessary
if (atom->nlocal > nmax) {
memory->destroy(list);
memory->destroy(mark);
nmax = atom->nmax;
memory->create(list,nmax,"evaporate:list");
memory->create(mark,nmax,"evaporate:mark");
}
// ncount = # of deletable atoms in region that I own
// nall = # on all procs
// nbefore = # on procs before me
// list[ncount] = list of local indices of atoms I can delete
double **x = atom->x;
int *mask = atom->mask;
int *tag = atom->tag;
int nlocal = atom->nlocal;
int ncount = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
list[ncount++] = i;
int nall,nbefore;
MPI_Allreduce(&ncount,&nall,1,MPI_INT,MPI_SUM,world);
MPI_Scan(&ncount,&nbefore,1,MPI_INT,MPI_SUM,world);
nbefore -= ncount;
// ndel = total # of atom deletions, in or out of region
// ndeltopo[1,2,3,4] = ditto for bonds, angles, dihedrals, impropers
// mark[] = 1 if deleted
ndel = 0;
for (i = 0; i < nlocal; i++) mark[i] = 0;
// atomic deletions
// choose atoms randomly across all procs and mark them for deletion
// shrink eligible list as my atoms get marked
// keep ndel,ncount,nall,nbefore current after each atom deletion
if (molflag == 0) {
while (nall && ndel < nflux) {
iwhichglobal = static_cast<int> (nall*random->uniform());
if (iwhichglobal < nbefore) nbefore--;
else if (iwhichglobal < nbefore + ncount) {
iwhichlocal = iwhichglobal - nbefore;
mark[list[iwhichlocal]] = 1;
list[iwhichlocal] = list[ncount-1];
ncount--;
}
ndel++;
nall--;
}
// molecule deletions
// choose one atom in one molecule randomly across all procs
// bcast mol ID and delete all atoms in that molecule on any proc
// update deletion count by total # of atoms in molecule
// shrink list of eligible candidates as any of my atoms get marked
// keep ndel,ndeltopo,ncount,nall,nbefore current after each mol deletion
} else {
int me,proc,iatom,imolecule,ndelone,ndelall;
int *molecule = atom->molecule;
ndeltopo[0] = ndeltopo[1] = ndeltopo[2] = ndeltopo[3] = 0;
while (nall && ndel < nflux) {
// pick an iatom,imolecule on proc me to delete
iwhichglobal = static_cast<int> (nall*random->uniform());
if (iwhichglobal >= nbefore && iwhichglobal < nbefore + ncount) {
iwhichlocal = iwhichglobal - nbefore;
iatom = list[iwhichlocal];
imolecule = molecule[iatom];
me = comm->me;
} else me = -1;
// bcast mol ID to delete all atoms from
// if mol ID > 0, delete any atom in molecule and decrement counters
// if mol ID == 0, delete single iatom
// be careful to delete correct # of bond,angle,etc for newton on or off
MPI_Allreduce(&me,&proc,1,MPI_INT,MPI_MAX,world);
MPI_Bcast(&imolecule,1,MPI_INT,proc,world);
ndelone = 0;
for (i = 0; i < nlocal; i++) {
if (imolecule && molecule[i] == imolecule) {
mark[i] = 1;
ndelone++;
if (atom->avec->bonds_allow) {
if (force->newton_bond) ndeltopo[0] += atom->num_bond[i];
else {
for (j = 0; j < atom->num_bond[i]; j++) {
if (tag[i] < atom->bond_atom[i][j]) ndeltopo[0]++;
}
}
}
if (atom->avec->angles_allow) {
if (force->newton_bond) ndeltopo[1] += atom->num_angle[i];
else {
for (j = 0; j < atom->num_angle[i]; j++) {
m = atom->map(atom->angle_atom2[i][j]);
if (m >= 0 && m < nlocal) ndeltopo[1]++;
}
}
}
if (atom->avec->dihedrals_allow) {
if (force->newton_bond) ndeltopo[2] += atom->num_dihedral[i];
else {
for (j = 0; j < atom->num_dihedral[i]; j++) {
m = atom->map(atom->dihedral_atom2[i][j]);
if (m >= 0 && m < nlocal) ndeltopo[2]++;
}
}
}
if (atom->avec->impropers_allow) {
if (force->newton_bond) ndeltopo[3] += atom->num_improper[i];
else {
for (j = 0; j < atom->num_improper[i]; j++) {
m = atom->map(atom->improper_atom2[i][j]);
if (m >= 0 && m < nlocal) ndeltopo[3]++;
}
}
}
} else if (me == proc && i == iatom) {
mark[i] = 1;
ndelone++;
}
}
// remove any atoms marked for deletion from my eligible list
i = 0;
while (i < ncount) {
if (mark[list[i]]) {
list[i] = list[ncount-1];
ncount--;
} else i++;
}
// update ndel,ncount,nall,nbefore
// ndelall is total atoms deleted on this iteration
// ncount is already correct, so resum to get nall and nbefore
MPI_Allreduce(&ndelone,&ndelall,1,MPI_INT,MPI_SUM,world);
ndel += ndelall;
MPI_Allreduce(&ncount,&nall,1,MPI_INT,MPI_SUM,world);
MPI_Scan(&ncount,&nbefore,1,MPI_INT,MPI_SUM,world);
nbefore -= ncount;
}
}
// delete my marked atoms
// loop in reverse order to avoid copying marked atoms
AtomVec *avec = atom->avec;
for (i = nlocal-1; i >= 0; i--) {
if (mark[i]) {
avec->copy(atom->nlocal-1,i,1);
atom->nlocal--;
}
}
// reset global natoms and bonds, angles, etc
// if global map exists, reset it now instead of waiting for comm
// since deleting atoms messes up ghosts
atom->natoms -= ndel;
if (molflag) {
int all[4];
MPI_Allreduce(ndeltopo,all,4,MPI_INT,MPI_SUM,world);
atom->nbonds -= all[0];
atom->nangles -= all[1];
atom->ndihedrals -= all[2];
atom->nimpropers -= all[3];
}
if (ndel && atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
// statistics
ndeleted += ndel;
next_reneighbor = update->ntimestep + nevery;
}
/* ----------------------------------------------------------------------
return number of deleted particles
------------------------------------------------------------------------- */
double FixEvaporate::compute_scalar()
{
return 1.0*ndeleted;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixEvaporate::memory_usage()
{
double bytes = 2*nmax * sizeof(int);
return bytes;
}

80
src/fix_evaporate.h
View File

@ -1,80 +0,0 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(evaporate,FixEvaporate)
#else
#ifndef LMP_FIX_EVAPORATE_H
#define LMP_FIX_EVAPORATE_H
#include "fix.h"
namespace LAMMPS_NS {
class FixEvaporate : public Fix {
public:
FixEvaporate(class LAMMPS *, int, char **);
~FixEvaporate();
int setmask();
void init();
void pre_exchange();
double compute_scalar();
double memory_usage();
private:
int nevery,nflux,iregion;
int molflag;
int ndeleted;
char *idregion;
int nmax;
int *list,*mark;
class RanPark *random;
};
}
#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: Region ID for fix evaporate does not exist
Self-explanatory.
E: Cannot evaporate atoms in atom_modify first group
This is a restriction due to the way atoms are organized in
a list to enable the atom_modify first command.
W: Fix evaporate may delete atom with non-zero molecule ID
This is probably an error, since you should not delete only one atom
of a molecule.
E: Fix evaporate molecule requires atom attribute molecule
The atom style being used does not define a molecule ID.
*/

330
src/fix_thermal_conductivity.cpp
View File

@ -1,330 +0,0 @@
/* ----------------------------------------------------------------------
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 author: Craig Tenney (UND) added support
for swapping atoms of different masses
------------------------------------------------------------------------- */
#include "math.h"
#include "mpi.h"
#include "string.h"
#include "stdlib.h"
#include "fix_thermal_conductivity.h"
#include "atom.h"
#include "force.h"
#include "domain.h"
#include "modify.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
#define BIG 1.0e10
/* ---------------------------------------------------------------------- */
FixThermalConductivity::FixThermalConductivity(LAMMPS *lmp,
int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 6) error->all(FLERR,"Illegal fix thermal/conductivity command");
MPI_Comm_rank(world,&me);
nevery = force->inumeric(FLERR,arg[3]);
if (nevery <= 0) error->all(FLERR,"Illegal fix thermal/conductivity command");
scalar_flag = 1;
global_freq = nevery;
extscalar = 0;
if (strcmp(arg[4],"x") == 0) edim = 0;
else if (strcmp(arg[4],"y") == 0) edim = 1;
else if (strcmp(arg[4],"z") == 0) edim = 2;
else error->all(FLERR,"Illegal fix thermal/conductivity command");
nbin = force->inumeric(FLERR,arg[5]);
if (nbin % 2 || nbin <= 2)
error->all(FLERR,"Illegal fix thermal/conductivity command");
// optional keywords
nswap = 1;
int iarg = 6;
while (iarg < narg) {
if (strcmp(arg[iarg],"swap") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix thermal/conductivity command");
nswap = force->inumeric(FLERR,arg[iarg+1]);
if (nswap <= 0)
error->all(FLERR,
"Fix thermal/conductivity swap value must be positive");
iarg += 2;
} else error->all(FLERR,"Illegal fix thermal/conductivity command");
}
// initialize array sizes to nswap+1 so have space to shift values down
index_lo = new int[nswap+1];
index_hi = new int[nswap+1];
ke_lo = new double[nswap+1];
ke_hi = new double[nswap+1];
e_exchange = 0.0;
}
/* ---------------------------------------------------------------------- */
FixThermalConductivity::~FixThermalConductivity()
{
delete [] index_lo;
delete [] index_hi;
delete [] ke_lo;
delete [] ke_hi;
}
/* ---------------------------------------------------------------------- */
int FixThermalConductivity::setmask()
{
int mask = 0;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixThermalConductivity::init()
{
// warn if any fix ave/spatial comes after this fix
// can cause glitch in averaging since ave will happen after swap
int foundme = 0;
for (int i = 0; i < modify->nfix; i++) {
if (modify->fix[i] == this) foundme = 1;
if (foundme && strcmp(modify->fix[i]->style,"ave/spatial") == 0 && me == 0)
error->warning(FLERR,
"Fix thermal/conductivity comes before fix ave/spatial");
}
// set bounds of 2 slabs in edim
// only necessary for static box, else re-computed in end_of_step()
// lo bin is always bottom bin
// hi bin is just above half height
if (domain->box_change == 0) {
prd = domain->prd[edim];
boxlo = domain->boxlo[edim];
boxhi = domain->boxhi[edim];
double binsize = (boxhi-boxlo) / nbin;
slablo_lo = boxlo;
slablo_hi = boxlo + binsize;
slabhi_lo = boxlo + (nbin/2)*binsize;
slabhi_hi = boxlo + (nbin/2+1)*binsize;
}
periodicity = domain->periodicity[edim];
}
/* ---------------------------------------------------------------------- */
void FixThermalConductivity::end_of_step()
{
int i,j,m,insert;
double coord,ke;
MPI_Status status;
struct {
double value;
int proc;
} mine[2],all[2];
// if box changes, recompute bounds of 2 slabs in edim
if (domain->box_change) {
prd = domain->prd[edim];
boxlo = domain->boxlo[edim];
boxhi = domain->boxhi[edim];
double binsize = (boxhi-boxlo) / nbin;
slablo_lo = boxlo;
slablo_hi = boxlo + binsize;
slabhi_lo = boxlo + (nbin/2)*binsize;
slabhi_hi = boxlo + (nbin/2+1)*binsize;
}
// make 2 lists of up to nswap atoms
// hottest atoms in lo slab, coldest atoms in hi slab (really mid slab)
// lo slab list is sorted by hottest, hi slab is sorted by coldest
// map atoms back into periodic box if necessary
// insert = location in list to insert new atom
double **x = atom->x;
double **v = atom->v;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
nhi = nlo = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
coord = x[i][edim];
if (coord < boxlo && periodicity) coord += prd;
else if (coord >= boxhi && periodicity) coord -= prd;
if (coord >= slablo_lo && coord < slablo_hi) {
ke = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
if (rmass) ke *= 0.5*rmass[i];
else ke *= 0.5*mass[type[i]];
if (nlo < nswap || ke > ke_lo[nswap-1]) {
for (insert = nlo-1; insert >= 0; insert--)
if (ke < ke_lo[insert]) break;
insert++;
for (m = nlo-1; m >= insert; m--) {
ke_lo[m+1] = ke_lo[m];
index_lo[m+1] = index_lo[m];
}
ke_lo[insert] = ke;
index_lo[insert] = i;
if (nlo < nswap) nlo++;
}
}
if (coord >= slabhi_lo && coord < slabhi_hi) {
ke = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
if (rmass) ke *= 0.5*rmass[i];
else ke *= 0.5*mass[type[i]];
if (nhi < nswap || ke < ke_hi[nswap-1]) {
for (insert = nhi-1; insert >= 0; insert--)
if (ke > ke_hi[insert]) break;
insert++;
for (m = nhi-1; m >= insert; m--) {
ke_hi[m+1] = ke_hi[m];
index_hi[m+1] = index_hi[m];
}
ke_hi[insert] = ke;
index_hi[insert] = i;
if (nhi < nswap) nhi++;
}
}
}
// loop over nswap pairs
// pair 2 global atoms at beginning of sorted lo/hi slab lists via Allreduce
// BIG values are for procs with no atom to contribute
// use negative of hottest KE since is doing a MINLOC
// MINLOC also communicates which procs own them
// exchange kinetic energy between the 2 particles
// if I own both particles just swap, else point2point comm of velocities
double sbuf[4],rbuf[4],vcm[3];
double eswap = 0.0;
mine[0].proc = mine[1].proc = me;
int ilo = 0;
int ihi = 0;
for (m = 0; m < nswap; m++) {
if (ilo < nlo) mine[0].value = -ke_lo[ilo];
else mine[0].value = BIG;
if (ihi < nhi) mine[1].value = ke_hi[ihi];
else mine[1].value = BIG;
MPI_Allreduce(mine,all,2,MPI_DOUBLE_INT,MPI_MINLOC,world);
if (all[0].value == BIG || all[1].value == BIG) continue;
if (me == all[0].proc && me == all[1].proc) {
i = index_lo[ilo++];
j = index_hi[ihi++];
sbuf[0] = v[j][0];
sbuf[1] = v[j][1];
sbuf[2] = v[j][2];
if (rmass) sbuf[3] = rmass[j];
else sbuf[3] = mass[type[j]];
rbuf[0] = v[i][0];
rbuf[1] = v[i][1];
rbuf[2] = v[i][2];
if (rmass) rbuf[3] = rmass[i];
else rbuf[3] = mass[type[i]];
vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]);
vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]);
vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]);
v[j][0] = 2.0 * vcm[0] - sbuf[0];
v[j][1] = 2.0 * vcm[1] - sbuf[1];
v[j][2] = 2.0 * vcm[2] - sbuf[2];
eswap += sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) +
vcm[1] * (vcm[1] - sbuf[1]) +
vcm[2] * (vcm[2] - sbuf[2]));
v[i][0] = 2.0 * vcm[0] - rbuf[0];
v[i][1] = 2.0 * vcm[1] - rbuf[1];
v[i][2] = 2.0 * vcm[2] - rbuf[2];
eswap -= rbuf[3] * (vcm[0] * (vcm[0] - rbuf[0]) +
vcm[1] * (vcm[1] - rbuf[1]) +
vcm[2] * (vcm[2] - rbuf[2]));
} else if (me == all[0].proc) {
j = index_lo[ilo++];
sbuf[0] = v[j][0];
sbuf[1] = v[j][1];
sbuf[2] = v[j][2];
if (rmass) sbuf[3] = rmass[j];
else sbuf[3] = mass[type[j]];
MPI_Sendrecv(sbuf,4,MPI_DOUBLE,all[1].proc,0,
rbuf,4,MPI_DOUBLE,all[1].proc,0,world,&status);
vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]);
vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]);
vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]);
v[j][0] = 2.0 * vcm[0] - sbuf[0];
v[j][1] = 2.0 * vcm[1] - sbuf[1];
v[j][2] = 2.0 * vcm[2] - sbuf[2];
eswap -= sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) +
vcm[1] * (vcm[1] - sbuf[1]) +
vcm[2] * (vcm[2] - sbuf[2]));
} else if (me == all[1].proc) {
j = index_hi[ihi++];
sbuf[0] = v[j][0];
sbuf[1] = v[j][1];
sbuf[2] = v[j][2];
if (rmass) sbuf[3] = rmass[j];
else sbuf[3] = mass[type[j]];
MPI_Sendrecv(sbuf,4,MPI_DOUBLE,all[0].proc,0,
rbuf,4,MPI_DOUBLE,all[0].proc,0,world,&status);
vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]);
vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]);
vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]);
v[j][0] = 2.0 * vcm[0] - sbuf[0];
v[j][1] = 2.0 * vcm[1] - sbuf[1];
v[j][2] = 2.0 * vcm[2] - sbuf[2];
eswap += sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) +
vcm[1] * (vcm[1] - sbuf[1]) +
vcm[2] * (vcm[2] - sbuf[2]));
}
}
// tally energy exchange from all swaps
double eswap_all;
MPI_Allreduce(&eswap,&eswap_all,1,MPI_DOUBLE,MPI_SUM,world);
e_exchange += force->mvv2e * eswap_all;
}
/* ---------------------------------------------------------------------- */
double FixThermalConductivity::compute_scalar()
{
return e_exchange;
}

75
src/fix_thermal_conductivity.h
View File

@ -1,75 +0,0 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(thermal/conductivity,FixThermalConductivity)
#else
#ifndef LMP_FIX_THERMAL_CONDUCTIVITY_H
#define LMP_FIX_THERMAL_CONDUCTIVITY_H
#include "fix.h"
namespace LAMMPS_NS {
class FixThermalConductivity : public Fix {
public:
FixThermalConductivity(class LAMMPS *, int, char **);
~FixThermalConductivity();
int setmask();
void init();
void end_of_step();
double compute_scalar();
private:
int me;
int edim,nbin,periodicity;
int nswap;
double prd,boxlo,boxhi;
double slablo_lo,slablo_hi,slabhi_lo,slabhi_hi;
double e_exchange;
int nlo,nhi;
int *index_lo,*index_hi;
double *ke_lo,*ke_hi;
};
}
#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: Fix thermal/conductivity swap value must be positive
Self-explanatory.
W: Fix thermal/conductivity comes before fix ave/spatial
The order of these 2 fixes in your input script is such that fix
thermal/conductivity comes first. If you are using fix ave/spatial to
measure the temperature profile induced by fix viscosity, then this
may cause a glitch in the profile since you are averaging immediately
after swaps have occurred. Flipping the order of the 2 fixes
typically helps.
*/

685
src/fix_ttm.cpp
View File

@ -1,685 +0,0 @@
/* ----------------------------------------------------------------------
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: Paul Crozier (SNL)
Carolyn Phillips (University of Michigan)
------------------------------------------------------------------------- */
#include "lmptype.h"
#include "mpi.h"
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "fix_ttm.h"
#include "atom.h"
#include "force.h"
#include "update.h"
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "comm.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
#define MAXLINE 1024
/* ---------------------------------------------------------------------- */
FixTTM::FixTTM(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 15) error->all(FLERR,"Illegal fix ttm command");
vector_flag = 1;
size_vector = 2;
global_freq = 1;
extvector = 1;
nevery = 1;
restart_peratom = 1;
restart_global = 1;
seed = force->inumeric(FLERR,arg[3]);
electronic_specific_heat = force->numeric(FLERR,arg[4]);
electronic_density = force->numeric(FLERR,arg[5]);
electronic_thermal_conductivity = force->numeric(FLERR,arg[6]);
gamma_p = force->numeric(FLERR,arg[7]);
gamma_s = force->numeric(FLERR,arg[8]);
v_0 = force->numeric(FLERR,arg[9]);
nxnodes = force->inumeric(FLERR,arg[10]);
nynodes = force->inumeric(FLERR,arg[11]);
nznodes = force->inumeric(FLERR,arg[12]);
fpr = fopen(arg[13],"r");
if (fpr == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",arg[13]);
error->one(FLERR,str);
}
nfileevery = force->inumeric(FLERR,arg[14]);
if (nfileevery) {
if (narg != 16) error->all(FLERR,"Illegal fix ttm command");
MPI_Comm_rank(world,&me);
if (me == 0) {
fp = fopen(arg[15],"w");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open fix ttm file %s",arg[15]);
error->one(FLERR,str);
}
}
}
// error check
if (seed <= 0) error->all(FLERR,"Invalid random number seed in fix ttm command");
if (electronic_specific_heat <= 0.0)
error->all(FLERR,"Fix ttm electronic_specific_heat must be > 0.0");
if (electronic_density <= 0.0)
error->all(FLERR,"Fix ttm electronic_density must be > 0.0");
if (electronic_thermal_conductivity < 0.0)
error->all(FLERR,"Fix ttm electronic_thermal_conductivity must be >= 0.0");
if (gamma_p <= 0.0) error->all(FLERR,"Fix ttm gamma_p must be > 0.0");
if (gamma_s < 0.0) error->all(FLERR,"Fix ttm gamma_s must be >= 0.0");
if (v_0 < 0.0) error->all(FLERR,"Fix ttm v_0 must be >= 0.0");
if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0)
error->all(FLERR,"Fix ttm number of nodes must be > 0");
v_0_sq = v_0*v_0;
// initialize Marsaglia RNG with processor-unique seed
random = new RanMars(lmp,seed + comm->me);
// allocate per-type arrays for force prefactors
gfactor1 = new double[atom->ntypes+1];
gfactor2 = new double[atom->ntypes+1];
// allocate 3d grid variables
total_nnodes = nxnodes*nynodes*nznodes;
memory->create(nsum,nxnodes,nynodes,nznodes,"ttm:nsum");
memory->create(nsum_all,nxnodes,nynodes,nznodes,"ttm:nsum_all");
memory->create(T_initial_set,nxnodes,nynodes,nznodes,"ttm:T_initial_set");
memory->create(sum_vsq,nxnodes,nynodes,nznodes,"ttm:sum_vsq");
memory->create(sum_mass_vsq,nxnodes,nynodes,nznodes,"ttm:sum_mass_vsq");
memory->create(sum_vsq_all,nxnodes,nynodes,nznodes,"ttm:sum_vsq_all");
memory->create(sum_mass_vsq_all,nxnodes,nynodes,nznodes,
"ttm:sum_mass_vsq_all");
memory->create(T_electron_old,nxnodes,nynodes,nznodes,"ttm:T_electron_old");
memory->create(T_electron,nxnodes,nynodes,nznodes,"ttm:T_electron");
memory->create(net_energy_transfer,nxnodes,nynodes,nznodes,
"TTM:net_energy_transfer");
memory->create(net_energy_transfer_all,nxnodes,nynodes,nznodes,
"TTM:net_energy_transfer_all");
flangevin = NULL;
grow_arrays(atom->nmax);
// zero out the flangevin array
for (int i = 0; i < atom->nmax; i++) {
flangevin[i][0] = 0;
flangevin[i][1] = 0;
flangevin[i][2] = 0;
}
atom->add_callback(0);
atom->add_callback(1);
// set initial electron temperatures from user input file
if (me == 0) read_initial_electron_temperatures();
MPI_Bcast(&T_electron[0][0][0],total_nnodes,MPI_DOUBLE,0,world);
}
/* ---------------------------------------------------------------------- */
FixTTM::~FixTTM()
{
if (nfileevery && me == 0) fclose(fp);
delete random;
delete [] gfactor1;
delete [] gfactor2;
memory->destroy(nsum);
memory->destroy(nsum_all);
memory->destroy(T_initial_set);
memory->destroy(sum_vsq);
memory->destroy(sum_mass_vsq);
memory->destroy(sum_vsq_all);
memory->destroy(sum_mass_vsq_all);
memory->destroy(T_electron_old);
memory->destroy(T_electron);
memory->destroy(flangevin);
memory->destroy(net_energy_transfer);
memory->destroy(net_energy_transfer_all);
}
/* ---------------------------------------------------------------------- */
int FixTTM::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixTTM::init()
{
if (domain->dimension == 2)
error->all(FLERR,"Cannot use fix ttm with 2d simulation");
if (domain->nonperiodic != 0)
error->all(FLERR,"Cannot use nonperiodic boundares with fix ttm");
if (domain->triclinic)
error->all(FLERR,"Cannot use fix ttm with triclinic box");
// set force prefactors
for (int i = 1; i <= atom->ntypes; i++) {
gfactor1[i] = - gamma_p / force->ftm2v;
gfactor2[i] =
sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
}
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
net_energy_transfer_all[ixnode][iynode][iznode] = 0;
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixTTM::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force_setup(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa_setup(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}
/* ---------------------------------------------------------------------- */
void FixTTM::post_force(int vflag)
{
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double gamma1,gamma2;
// apply damping and thermostat to all atoms in fix group
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
int ixnode = static_cast<int>(xscale*nxnodes);
int iynode = static_cast<int>(yscale*nynodes);
int iznode = static_cast<int>(zscale*nznodes);
while (ixnode > nxnodes-1) ixnode -= nxnodes;
while (iynode > nynodes-1) iynode -= nynodes;
while (iznode > nznodes-1) iznode -= nznodes;
while (ixnode < 0) ixnode += nxnodes;
while (iynode < 0) iynode += nynodes;
while (iznode < 0) iznode += nznodes;
if (T_electron[ixnode][iynode][iznode] < 0)
error->all(FLERR,"Electronic temperature dropped below zero");
double tsqrt = sqrt(T_electron[ixnode][iynode][iznode]);
gamma1 = gfactor1[type[i]];
double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
if (vsq > v_0_sq) gamma1 *= (gamma_p + gamma_s)/gamma_p;
gamma2 = gfactor2[type[i]] * tsqrt;
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixTTM::post_force_setup(int vflag)
{
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// apply langevin forces that have been stored from previous run
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixTTM::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixTTM::post_force_respa_setup(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force_setup(vflag);
}
/* ---------------------------------------------------------------------- */
void FixTTM::reset_dt()
{
for (int i = 1; i <= atom->ntypes; i++)
gfactor2[i] =
sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
}
/* ----------------------------------------------------------------------
read in initial electron temperatures from a user-specified file
only called by proc 0
------------------------------------------------------------------------- */
void FixTTM::read_initial_electron_temperatures()
{
char line[MAXLINE];
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
T_initial_set[ixnode][iynode][iznode] = 0;
// read initial electron temperature values from file
int ixnode,iynode,iznode;
double T_tmp;
while (1) {
if (fgets(line,MAXLINE,fpr) == NULL) break;
sscanf(line,"%d %d %d %lg",&ixnode,&iynode,&iznode,&T_tmp);
if (T_tmp < 0.0) error->one(FLERR,"Fix ttm electron temperatures must be > 0.0");
T_electron[ixnode][iynode][iznode] = T_tmp;
T_initial_set[ixnode][iynode][iznode] = 1;
}
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
if (T_initial_set[ixnode][iynode][iznode] == 0)
error->one(FLERR,"Initial temperatures not all set in fix ttm");
// close file
fclose(fpr);
}
/* ---------------------------------------------------------------------- */
void FixTTM::end_of_step()
{
double **x = atom->x;
double **v = atom->v;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
net_energy_transfer[ixnode][iynode][iznode] = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
int ixnode = static_cast<int>(xscale*nxnodes);
int iynode = static_cast<int>(yscale*nynodes);
int iznode = static_cast<int>(zscale*nznodes);
while (ixnode > nxnodes-1) ixnode -= nxnodes;
while (iynode > nynodes-1) iynode -= nynodes;
while (iznode > nznodes-1) iznode -= nznodes;
while (ixnode < 0) ixnode += nxnodes;
while (iynode < 0) iynode += nynodes;
while (iznode < 0) iznode += nznodes;
net_energy_transfer[ixnode][iynode][iznode] +=
(flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2]);
}
MPI_Allreduce(&net_energy_transfer[0][0][0],
&net_energy_transfer_all[0][0][0],
total_nnodes,MPI_DOUBLE,MPI_SUM,world);
double dx = domain->xprd/nxnodes;
double dy = domain->yprd/nynodes;
double dz = domain->zprd/nznodes;
double del_vol = dx*dy*dz;
// num_inner_timesteps = # of inner steps (thermal solves)
// required this MD step to maintain a stable explicit solve
int num_inner_timesteps = 1;
double inner_dt = update->dt;
double stability_criterion = 1.0 -
2.0*inner_dt/(electronic_specific_heat*electronic_density) *
(electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
if (stability_criterion < 0.0) {
inner_dt = 0.5*(electronic_specific_heat*electronic_density) /
(electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
num_inner_timesteps = static_cast<int>(update->dt/inner_dt) + 1;
inner_dt = update->dt/double(num_inner_timesteps);
if (num_inner_timesteps > 1000000)
error->warning(FLERR,"Too many inner timesteps in fix ttm",0);
}
for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps;
ith_inner_timestep++) {
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
T_electron_old[ixnode][iynode][iznode] =
T_electron[ixnode][iynode][iznode];
// compute new electron T profile
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
int right_xnode = ixnode + 1;
int right_ynode = iynode + 1;
int right_znode = iznode + 1;
if (right_xnode == nxnodes) right_xnode = 0;
if (right_ynode == nynodes) right_ynode = 0;
if (right_znode == nznodes) right_znode = 0;
int left_xnode = ixnode - 1;
int left_ynode = iynode - 1;
int left_znode = iznode - 1;
if (left_xnode == -1) left_xnode = nxnodes - 1;
if (left_ynode == -1) left_ynode = nynodes - 1;
if (left_znode == -1) left_znode = nznodes - 1;
T_electron[ixnode][iynode][iznode] =
T_electron_old[ixnode][iynode][iznode] +
inner_dt/(electronic_specific_heat*electronic_density) *
(electronic_thermal_conductivity *
((T_electron_old[right_xnode][iynode][iznode] +
T_electron_old[left_xnode][iynode][iznode] -
2*T_electron_old[ixnode][iynode][iznode])/dx/dx +
(T_electron_old[ixnode][right_ynode][iznode] +
T_electron_old[ixnode][left_ynode][iznode] -
2*T_electron_old[ixnode][iynode][iznode])/dy/dy +
(T_electron_old[ixnode][iynode][right_znode] +
T_electron_old[ixnode][iynode][left_znode] -
2*T_electron_old[ixnode][iynode][iznode])/dz/dz) -
(net_energy_transfer_all[ixnode][iynode][iznode])/del_vol);
}
}
// output nodal temperatures for current timestep
if ((nfileevery) && !(update->ntimestep % nfileevery)) {
// compute atomic Ta for each grid point
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
nsum[ixnode][iynode][iznode] = 0;
nsum_all[ixnode][iynode][iznode] = 0;
sum_vsq[ixnode][iynode][iznode] = 0.0;
sum_mass_vsq[ixnode][iynode][iznode] = 0.0;
sum_vsq_all[ixnode][iynode][iznode] = 0.0;
sum_mass_vsq_all[ixnode][iynode][iznode] = 0.0;
}
double massone;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
int ixnode = static_cast<int>(xscale*nxnodes);
int iynode = static_cast<int>(yscale*nynodes);
int iznode = static_cast<int>(zscale*nznodes);
while (ixnode > nxnodes-1) ixnode -= nxnodes;
while (iynode > nynodes-1) iynode -= nynodes;
while (iznode > nznodes-1) iznode -= nznodes;
while (ixnode < 0) ixnode += nxnodes;
while (iynode < 0) iynode += nynodes;
while (iznode < 0) iznode += nznodes;
double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
nsum[ixnode][iynode][iznode] += 1;
sum_vsq[ixnode][iynode][iznode] += vsq;
sum_mass_vsq[ixnode][iynode][iznode] += massone*vsq;
}
MPI_Allreduce(&nsum[0][0][0],&nsum_all[0][0][0],total_nnodes,
MPI_INT,MPI_SUM,world);
MPI_Allreduce(&sum_vsq[0][0][0],&sum_vsq_all[0][0][0],total_nnodes,
MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&sum_mass_vsq[0][0][0],&sum_mass_vsq_all[0][0][0],
total_nnodes,MPI_DOUBLE,MPI_SUM,world);
if (me == 0) {
fprintf(fp,BIGINT_FORMAT,update->ntimestep);
double T_a;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
T_a = 0;
if (nsum_all[ixnode][iynode][iznode] > 0)
T_a = sum_mass_vsq_all[ixnode][iynode][iznode]/
(3.0*force->boltz*nsum_all[ixnode][iynode][iznode]/force->mvv2e);
fprintf(fp," %f",T_a);
}
fprintf(fp,"\t");
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
fprintf(fp,"%f ",T_electron[ixnode][iynode][iznode]);
fprintf(fp,"\n");
}
}
}
/* ----------------------------------------------------------------------
memory usage of 3d grid
------------------------------------------------------------------------- */
double FixTTM::memory_usage()
{
double bytes = 0.0;
bytes += 5*total_nnodes * sizeof(int);
bytes += 14*total_nnodes * sizeof(double);
return bytes;
}
/* ---------------------------------------------------------------------- */
void FixTTM::grow_arrays(int ngrow)
{
memory->grow(flangevin,ngrow,3,"TTM:flangevin");
}
/* ----------------------------------------------------------------------
return the energy of the electronic subsystem or the net_energy transfer
between the subsystems
------------------------------------------------------------------------- */
double FixTTM::compute_vector(int n)
{
double e_energy = 0.0;
double transfer_energy = 0.0;
double dx = domain->xprd/nxnodes;
double dy = domain->yprd/nynodes;
double dz = domain->zprd/nznodes;
double del_vol = dx*dy*dz;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
e_energy +=
T_electron[ixnode][iynode][iznode]*electronic_specific_heat*
electronic_density*del_vol;
transfer_energy +=
net_energy_transfer_all[ixnode][iynode][iznode]*update->dt;
}
if (n == 0) return e_energy;
if (n == 1) return transfer_energy;
return 0.0;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixTTM::write_restart(FILE *fp)
{
double *rlist;
memory->create(rlist,nxnodes*nynodes*nznodes+1,"TTM:rlist");
int n = 0;
rlist[n++] = seed;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
rlist[n++] = T_electron[ixnode][iynode][iznode];
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(rlist,sizeof(double),n,fp);
}
memory->destroy(rlist);
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixTTM::restart(char *buf)
{
int n = 0;
double *rlist = (double *) buf;
// the seed must be changed from the initial seed
seed = static_cast<int> (0.5*rlist[n++]);
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
T_electron[ixnode][iynode][iznode] = rlist[n++];
delete random;
random = new RanMars(lmp,seed+comm->me);
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for restart file
------------------------------------------------------------------------- */
int FixTTM::pack_restart(int i, double *buf)
{
buf[0] = 4;
buf[1] = flangevin[i][0];
buf[2] = flangevin[i][1];
buf[3] = flangevin[i][2];
return 4;
}
/* ----------------------------------------------------------------------
unpack values from atom->extra array to restart the fix
------------------------------------------------------------------------- */
void FixTTM::unpack_restart(int nlocal, int nth)
{
double **extra = atom->extra;
// skip to Nth set of extra values
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;
flangevin[nlocal][0] = extra[nlocal][m++];
flangevin[nlocal][1] = extra[nlocal][m++];
flangevin[nlocal][2] = extra[nlocal][m++];
}
/* ----------------------------------------------------------------------
maxsize of any atom's restart data
------------------------------------------------------------------------- */
int FixTTM::maxsize_restart()
{
return 4;
}
/* ----------------------------------------------------------------------
size of atom nlocal's restart data
------------------------------------------------------------------------- */
int FixTTM::size_restart(int nlocal)
{
return 4;
}

156
src/fix_ttm.h
View File

@ -1,156 +0,0 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(ttm,FixTTM)
#else
#ifndef LMP_FIX_TTM_H
#define LMP_FIX_TTM_H
#include "fix.h"
namespace LAMMPS_NS {
class FixTTM : public Fix {
public:
FixTTM(class LAMMPS *, int, char **);
~FixTTM();
int setmask();
void init();
void setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void post_force_setup(int);
void post_force_respa_setup(int, int, int);
void end_of_step();
void reset_dt();
void write_restart(FILE *);
void restart(char *);
int pack_restart(int, double *);
void unpack_restart(int, int);
int size_restart(int);
int maxsize_restart();
double memory_usage();
void grow_arrays(int);
double compute_vector(int);
private:
int me;
int nfileevery;
int nlevels_respa;
int seed;
class RanMars *random;
FILE *fp,*fpr;
int nxnodes,nynodes,nznodes,total_nnodes;
int ***nsum;
int ***nsum_all,***T_initial_set;
double *gfactor1,*gfactor2,*ratio;
double **flangevin;
double ***T_electron,***T_electron_old;
double ***sum_vsq,***sum_mass_vsq;
double ***sum_vsq_all,***sum_mass_vsq_all;
double ***net_energy_transfer,***net_energy_transfer_all;
double electronic_specific_heat,electronic_density;
double electronic_thermal_conductivity;
double gamma_p,gamma_s,v_0,v_0_sq;
void read_initial_electron_temperatures();
};
}
#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: Cannot open file %s
The specified file cannot be opened. Check that the path and name are
correct.
E: Cannot open fix ttm file %s
The output file for the fix ttm command cannot be opened. Check that
the path and name are correct.
E: Invalid random number seed in fix ttm command
Random number seed must be > 0.
E: Fix ttm electronic_specific_heat must be > 0.0
Self-explanatory.
E: Fix ttm electronic_density must be > 0.0
Self-explanatory.
E: Fix ttm electronic_thermal_conductivity must be >= 0.0
Self-explanatory.
E: Fix ttm gamma_p must be > 0.0
Self-explanatory.
E: Fix ttm gamma_s must be >= 0.0
Self-explanatory.
E: Fix ttm v_0 must be >= 0.0
Self-explanatory.
E: Fix ttm number of nodes must be > 0
Self-explanatory.
E: Cannot use fix ttm with 2d simulation
This is a current restriction of this fix due to the grid it creates.
E: Cannot use nonperiodic boundares with fix ttm
This fix requires a fully periodic simulation box.
E: Cannot use fix ttm with triclinic box
This is a current restriction of this fix due to the grid it creates.
E: Electronic temperature dropped below zero
Something has gone wrong with the fix ttm electron temperature model.
E: Fix ttm electron temperatures must be > 0.0
Self-explanatory.
E: Initial temperatures not all set in fix ttm
Self-explantory.
W: Too many inner timesteps in fix ttm
Self-explanatory.
*/

309
src/fix_viscosity.cpp
View File

@ -1,309 +0,0 @@
/* ----------------------------------------------------------------------
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 author: Craig Tenney (UND) added support
for swapping atoms of different masses
------------------------------------------------------------------------- */
#include "math.h"
#include "mpi.h"
#include "string.h"
#include "stdlib.h"
#include "fix_viscosity.h"
#include "atom.h"
#include "domain.h"
#include "modify.h"
#include "error.h"
#include "force.h"
using namespace LAMMPS_NS;
using namespace FixConst;
// needs to be big, but not so big that lose precision when subtract velocity
#define BIG 1.0e10
/* ---------------------------------------------------------------------- */
FixViscosity::FixViscosity(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) error->all(FLERR,"Illegal fix viscosity command");
MPI_Comm_rank(world,&me);
nevery = force->inumeric(FLERR,arg[3]);
if (nevery <= 0) error->all(FLERR,"Illegal fix viscosity command");
scalar_flag = 1;
global_freq = nevery;
extscalar = 0;
if (strcmp(arg[4],"x") == 0) vdim = 0;
else if (strcmp(arg[4],"y") == 0) vdim = 1;
else if (strcmp(arg[4],"z") == 0) vdim = 2;
else error->all(FLERR,"Illegal fix viscosity command");
if (strcmp(arg[5],"x") == 0) pdim = 0;
else if (strcmp(arg[5],"y") == 0) pdim = 1;
else if (strcmp(arg[5],"z") == 0) pdim = 2;
else error->all(FLERR,"Illegal fix viscosity command");
nbin = force->inumeric(FLERR,arg[6]);
if (nbin % 2 || nbin <= 2) error->all(FLERR,"Illegal fix viscosity command");
// optional keywords
nswap = 1;
vtarget = BIG;
int iarg = 7;
while (iarg < narg) {
if (strcmp(arg[iarg],"swap") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix viscosity command");
nswap = force->inumeric(FLERR,arg[iarg+1]);
if (nswap <= 0)
error->all(FLERR,"Fix viscosity swap value must be positive");
iarg += 2;
} else if (strcmp(arg[iarg],"vtarget") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix viscosity command");
if (strcmp(arg[iarg+1],"INF") == 0) vtarget = BIG;
else vtarget = force->numeric(FLERR,arg[iarg+1]);
if (vtarget <= 0.0)
error->all(FLERR,"Fix viscosity vtarget value must be positive");
iarg += 2;
} else error->all(FLERR,"Illegal fix viscosity command");
}
// initialize array sizes to nswap+1 so have space to shift values down
pos_index = new int[nswap+1];
neg_index = new int[nswap+1];
pos_delta = new double[nswap+1];
neg_delta = new double[nswap+1];
p_exchange = 0.0;
}
/* ---------------------------------------------------------------------- */
FixViscosity::~FixViscosity()
{
delete [] pos_index;
delete [] neg_index;
delete [] pos_delta;
delete [] neg_delta;
}
/* ---------------------------------------------------------------------- */
int FixViscosity::setmask()
{
int mask = 0;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixViscosity::init()
{
// warn if any fix ave/spatial comes after this fix
// can cause glitch in averaging since ave will happen after swap
int foundme = 0;
for (int i = 0; i < modify->nfix; i++) {
if (modify->fix[i] == this) foundme = 1;
if (foundme && strcmp(modify->fix[i]->style,"ave/spatial") == 0 && me == 0)
error->warning(FLERR,"Fix viscosity comes before fix ave/spatial");
}
// set bounds of 2 slabs in pdim
// only necessary for static box, else re-computed in end_of_step()
// lo bin is always bottom bin
// hi bin is just above half height
if (domain->box_change == 0) {
prd = domain->prd[pdim];
boxlo = domain->boxlo[pdim];
boxhi = domain->boxhi[pdim];
double binsize = (boxhi-boxlo) / nbin;
slablo_lo = boxlo;
slablo_hi = boxlo + binsize;
slabhi_lo = boxlo + (nbin/2)*binsize;
slabhi_hi = boxlo + (nbin/2+1)*binsize;
}
periodicity = domain->periodicity[pdim];
}
/* ---------------------------------------------------------------------- */
void FixViscosity::end_of_step()
{
int i,m,insert;
double coord,delta;
MPI_Status status;
struct {
double value;
int proc;
} mine[2],all[2];
// if box changes, recompute bounds of 2 slabs in pdim
if (domain->box_change) {
prd = domain->prd[pdim];
boxlo = domain->boxlo[pdim];
boxhi = domain->boxhi[pdim];
double binsize = (boxhi-boxlo) / nbin;
slablo_lo = boxlo;
slablo_hi = boxlo + binsize;
slabhi_lo = boxlo + (nbin/2)*binsize;
slabhi_hi = boxlo + (nbin/2+1)*binsize;
}
// make 2 lists of up to nswap atoms with velocity closest to +/- vtarget
// lists are sorted by closeness to vtarget
// only consider atoms in the bottom/middle slabs
// map atoms back into periodic box if necessary
// insert = location in list to insert new atom
double **x = atom->x;
double **v = atom->v;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
npositive = nnegative = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
coord = x[i][pdim];
if (coord < boxlo && periodicity) coord += prd;
else if (coord >= boxhi && periodicity) coord -= prd;
if (coord >= slablo_lo && coord < slablo_hi) {
if (v[i][vdim] < 0.0) continue;
delta = fabs(v[i][vdim] - vtarget);
if (npositive < nswap || delta < pos_delta[nswap-1]) {
for (insert = npositive-1; insert >= 0; insert--)
if (delta > pos_delta[insert]) break;
insert++;
for (m = npositive-1; m >= insert; m--) {
pos_delta[m+1] = pos_delta[m];
pos_index[m+1] = pos_index[m];
}
pos_delta[insert] = delta;
pos_index[insert] = i;
if (npositive < nswap) npositive++;
}
}
if (coord >= slabhi_lo && coord < slabhi_hi) {
if (v[i][vdim] > 0.0) continue;
delta = fabs(v[i][vdim] + vtarget);
if (nnegative < nswap || delta < neg_delta[nswap-1]) {
for (insert = nnegative-1; insert >= 0; insert--)
if (delta > neg_delta[insert]) break;
insert++;
for (m = nnegative-1; m >= insert; m--) {
neg_delta[m+1] = neg_delta[m];
neg_index[m+1] = neg_index[m];
}
neg_delta[insert] = delta;
neg_index[insert] = i;
if (nnegative < nswap) nnegative++;
}
}
}
// loop over nswap pairs
// find 2 global atoms with smallest delta in bottom/top slabs
// BIG values are for procs with no atom to contribute
// MINLOC also communicates which procs own them
// exchange momenta between the 2 particles
// if I own both particles just swap, else point2point comm of vel,mass
double *mass = atom->mass;
double *rmass = atom->rmass;
int ipos,ineg;
double sbuf[2],rbuf[2],vcm;
double pswap = 0.0;
mine[0].proc = mine[1].proc = me;
int ipositive = 0;
int inegative = 0;
for (m = 0; m < nswap; m++) {
if (ipositive < npositive) mine[0].value = pos_delta[ipositive];
else mine[0].value = BIG;
if (inegative < nnegative) mine[1].value = neg_delta[inegative];
else mine[1].value = BIG;
MPI_Allreduce(mine,all,2,MPI_DOUBLE_INT,MPI_MINLOC,world);
if (all[0].value == BIG || all[1].value == BIG) continue;
if (me == all[0].proc && me == all[1].proc) {
ipos = pos_index[ipositive++];
ineg = neg_index[inegative++];
rbuf[0] = v[ipos][vdim];
if (rmass) rbuf[1] = rmass[ipos];
else rbuf[1] = mass[type[ipos]];
sbuf[0] = v[ineg][vdim];
if (rmass) sbuf[1] = rmass[ineg];
else sbuf[1] = mass[type[ineg]];
vcm = (sbuf[1]*sbuf[0] + rbuf[1]*rbuf[0]) / (sbuf[1] + rbuf[1]);
v[ineg][vdim] = 2.0 * vcm - sbuf[0];
v[ipos][vdim] = 2.0 * vcm - rbuf[0];
pswap += rbuf[1] * (vcm - rbuf[0]) - sbuf[1] * (vcm - sbuf[0]);
} else if (me == all[0].proc) {
ipos = pos_index[ipositive++];
sbuf[0] = v[ipos][vdim];
if (rmass) sbuf[1] = rmass[ipos];
else sbuf[1] = mass[type[ipos]];
MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[1].proc,0,
rbuf,2,MPI_DOUBLE,all[1].proc,0,world,&status);
vcm = (sbuf[1]*sbuf[0] + rbuf[1]*rbuf[0]) / (sbuf[1] + rbuf[1]);
v[ipos][vdim] = 2.0 * vcm - sbuf[0];
pswap += sbuf[1] * (vcm - sbuf[0]);
} else if (me == all[1].proc) {
ineg = neg_index[inegative++];
sbuf[0] = v[ineg][vdim];
if (rmass) sbuf[1] = rmass[ineg];
else sbuf[1] = mass[type[ineg]];
MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[0].proc,0,
rbuf,2,MPI_DOUBLE,all[0].proc,0,world,&status);
vcm = (sbuf[1]*sbuf[0] + rbuf[1]*rbuf[0]) / (sbuf[1] + rbuf[1]);
v[ineg][vdim] = 2.0 * vcm - sbuf[0];
pswap -= sbuf[1] * (vcm - sbuf[0]);
}
}
// tally momentum exchange from all swaps
double pswap_all;
MPI_Allreduce(&pswap,&pswap_all,1,MPI_DOUBLE,MPI_SUM,world);
p_exchange += pswap_all;
}
/* ---------------------------------------------------------------------- */
double FixViscosity::compute_scalar()
{
return p_exchange;
}

80
src/fix_viscosity.h
View File

@ -1,80 +0,0 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(viscosity,FixViscosity)
#else
#ifndef LMP_FIX_VISCOSITY_H
#define LMP_FIX_VISCOSITY_H
#include "fix.h"
namespace LAMMPS_NS {
class FixViscosity : public Fix {
public:
FixViscosity(class LAMMPS *, int, char **);
~FixViscosity();
int setmask();
void init();
void end_of_step();
double compute_scalar();
private:
int me;
int vdim,pdim,nbin,periodicity;
int nswap;
double vtarget;
double prd,boxlo,boxhi;
double slablo_lo,slablo_hi,slabhi_lo,slabhi_hi;
double p_exchange;
int npositive,nnegative;
int *pos_index,*neg_index;
double *pos_delta,*neg_delta;
};
}
#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: Fix viscosity swap value must be positive
Self-explanatory.
E: Fix viscosity vtarget value must be positive
Self-explanatory.
W: Fix viscosity comes before fix ave/spatial
The order of these 2 fixes in your input script is such that
fix viscosity comes first. If you are using fix ave/spatial
to measure the velocity profile induced by fix viscosity, then
this may cause a glitch in the profile since you are averaging
immediately after swaps have occurred. Flipping the order
of the 2 fixes typically helps.
*/