Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F87224129
fix_gaussflow.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Oct 11, 09:17
Size
6 KB
Mime Type
text/x-c
Expires
Sun, Oct 13, 09:17 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
21551481
Attached To
rLAMMPS lammps
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];
}
Event Timeline
Log In to Comment