Homec4science
Diffusion lammps (lammm-devel) Browse

lammps/src/fix_gaussflow.cpp
2e7404c96c42lammm-devel

Last Change
Raw File

fix_gaussflow.cpp
View Options

/* ----------------------------------------------------------------------
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: Steven E. Strong and Joel D. Eaves
Joel.Eaves@Colorado.edu
------------------------------------------------------------------------- */
#include "stdlib.h"
#include "string.h"
#include "fix_gaussflow.h"
#include "atom.h"
#include "force.h"
#include "group.h"
#include "comm.h"
#include "modify.h"
#include "update.h"
#include "domain.h"
#include "error.h"
#include "citeme.h"
#include "velocity.h"
using namespace LAMMPS_NS;
using namespace FixConst;
static const char cite_gaussflow[] =
"Gaussian dynamics package:\n\n"
"@Article{strong_atomistic_2016,\n"
"title = {Atomistic Hydrodynamics and the Dynamical Hydrophobic Effect in Porous Graphene},\n"
"volume = {7},\n"
"number = {10},\n"
"issn = {1948-7185},\n"
"url = {http://dx.doi.org/10.1021/acs.jpclett.6b00748},\n"
"doi = {10.1021/acs.jpclett.6b00748},\n"
"urldate = {2016-05-10},\n"
"journal = {J. Phys. Chem. Lett.},\n"
"author = {Strong, Steven E. and Eaves, Joel D.},\n"
"year = {2016},\n"
"pages = {1907--1912}\n"
"}\n\n";
// ID is arg[0]
//fix ID group gaussFlow Xf Yf Zf (opt: energy)
//Conserve a pre-applied flux in the X Y and Z directions depending on flags.
//if Pcorr keyword is given, then correct pressure tensor as discussed in ref
FixGaussFlow::FixGaussFlow(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
if (lmp->citeme) lmp->citeme->add(cite_gaussflow);
if (narg < 6) error->all(FLERR,"Not enough input arguments");
dynamic_group_allow = 0; //group which conserves momentum must also conserve particle number
scalar_flag = 1;
vector_flag = 1;
extscalar = 1;
extvector = 1;
size_vector = 3;
global_freq = 1; //data available every timestep
dimension = domain->dimension;
//get inputs
int tmpFlag;
for (int ii=0; ii<3; ii++)
{
tmpFlag=force->inumeric(FLERR,arg[3+ii]);
if (tmpFlag==1 || tmpFlag==0)
flow[ii]=tmpFlag;
else
error->all(FLERR,"Constraint flags must be 1 or 0");
}
//by default, do not compute work done
workflag=0;
// process optional keyword
int iarg = 6;
while (iarg < narg) {
if ( strcmp(arg[iarg],"energy") == 0 ) {
if ( iarg+2 > narg ) error->all(FLERR,"Illegal energy keyword");
if ( strcmp(arg[iarg+1],"yes") == 0 ) workflag = 1;
else if ( strcmp(arg[iarg+1],"no") == 1 ) error->all(FLERR,"Illegal energy keyword");
iarg += 2;
} else error->all(FLERR,"Illegal fix gaussFlow command");
}
//error checking
if (dimension == 2) {
if (flow[2])
error->all(FLERR,"Can't constrain z flow in 2d simulation");
}
dt=update->dt;
pe_tot=0.0;
}
/* ---------------------------------------------------------------------- */
int FixGaussFlow::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= THERMO_ENERGY;
return mask;
}
/* ----------------------------------------------------------------------
setup is called after the initial evaluation of forces before a run, so we
must remove the total force here too
------------------------------------------------------------------------- */
void FixGaussFlow::setup(int vflag)
{
//need to compute work done if set fix_modify energy yes
if (thermo_energy)
workflag=1;
//get total mass of group
mTot=group->mass(igroup);
post_force(vflag);
}
/* ----------------------------------------------------------------------
this is where Gaussian dynamics constraint is applied
------------------------------------------------------------------------- */
void FixGaussFlow::post_force(int vflag)
{
double **f = atom->f;
double **x = atom->x;
double **v = atom->v;
int *mask = atom->mask;
int *type = atom->type;
double *mass = atom->mass;
int nlocal = atom->nlocal;
int ii,jj;
//find the total force on all atoms
//initialize to zero
double f_thisProc[3];
for (ii=0; ii<3; ii++)
f_thisProc[ii]=0.0;
//add all forces on each processor
for(ii=0; ii<nlocal; ii++)
if (mask[ii] & groupbit)
for (jj=0; jj<3; jj++)
if (flow[jj])
f_thisProc[jj] += f[ii][jj];
//add the processor sums together
MPI_Allreduce(f_thisProc, f_tot, 3, MPI_DOUBLE, MPI_SUM, world);
//compute applied acceleration
for (ii=0; ii<3; ii++)
a_app[ii] = -f_tot[ii] / mTot;
//apply added accelleration to each atom
double f_app[3];
peAdded=0.0;
for( ii = 0; ii<nlocal; ii++)
if (mask[ii] & groupbit) {
f_app[0] = a_app[0]*mass[type[ii]];
f_app[1] = a_app[1]*mass[type[ii]];
f_app[2] = a_app[2]*mass[type[ii]];
f[ii][0] += f_app[0]; //f_app[jj] is 0 if flow[jj] is false
f[ii][1] += f_app[1];
f[ii][2] += f_app[2];
//calculate added energy, since more costly, only do this if requested
if (workflag)
peAdded += f_app[0]*v[ii][0] + f_app[1]*v[ii][1] + f_app[2]*v[ii][2];
}
//finish calculation of work done, sum over all procs
if (workflag) {
double pe_tmp=0.0;
MPI_Allreduce(&peAdded,&pe_tmp,1,MPI_DOUBLE,MPI_SUM,world);
pe_tot += pe_tmp;
}
}
/* ----------------------------------------------------------------------
negative of work done by this fix
This is only computed if requested, either with fix_modify energy yes, or with the energy keyword. Otherwise returns 0.
------------------------------------------------------------------------- */
double FixGaussFlow::compute_scalar()
{
return -pe_tot*dt;
}
/* ----------------------------------------------------------------------
return components of applied force
------------------------------------------------------------------------- */
double FixGaussFlow::compute_vector(int n)
{
return -f_tot[n];
}
c4science · Help