Merge branch 'merge-pull-153' into lammps-icms

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

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

:link(ld,Manual.html)

:link(lc,Section_commands.html#comm)

:line

fix flow/gauss command :h3

[Syntax:]

fix ID group-ID flow/gauss xflag yflag zflag keyword :pre

ID, group-ID are documented in "fix"_fix.html command :ulb,l

flow/gauss = style name of this fix command :l

xflag,yflag,zflag = 0 or 1 :l

    0 = do not conserve current in this dimension

    1 = conserve current in this dimension :pre

zero or more keyword/value pairs may be appended :l

keyword = {energy} :l

  {energy} value = no or yes

    no = do not compute work done by this fix

    yes = compute work done by this fix :pre

:ule

[Examples:]

fix GD fluid flow/gauss 1 0 0 

fix GD fluid flow/gauss 1 1 1 energy yes :pre

[Description:]

This fix implements the Gaussian dynamics (GD) method to simulate a system at 

constant mass flux "(Strong)"_#Strong. GD is a nonequilibrium molecular 

dynamics simulation method that can be used to study fluid flows through  

pores, pipes, and channels. In its original implementation GD was used to 

compute the pressure required to achieve a fixed mass flux through an opening. 

The flux can be conserved in any combination of the directions, x, y, or z, 

using xflag,yflag,zflag. This fix does not initialize a net flux through 

a system, it only conserves the center-of-mass momentum that is present 

when the fix is declared in the input script. Use the "velocity"_velocity.html 

command to generate an initial center-of-mass momentum.

GD applies an external fluctuating gravitational field that acts as a driving 

force to keep the system away from equilibrium. To maintain steady state, a 

profile-unbiased thermostat must be implemented to dissipate the heat that is 

added by the driving force. "Compute temp/profile"_compute_temp_profile.html 

can be used to implement a profile-unbiased thermostat. 

A common use of this fix is to compute a pressure drop across a pipe, pore, or 

membrane. The pressure profile can be computed in LAMMPS with "compute 

stress/atom"_compute_stress_atom.html and "fix ave/chunk"_fix_ave_chunk.html, 

or with the hardy method in "fix atc"_fix_atc.html. Note that the simple 

"compute stress/atom"_compute_stress_atom.html method is only accurate away 

from inhomogeneities in the fluid, such as fixed wall atoms. Further, the 

computed pressure profile must be corrected for the acceleration applied by 

GD before computing a pressure drop or comparing it to other methods, such as 

the pump method "(Zhu)"_#Zhu. The pressure correction is discussed and 

described in "(Strong)"_#Strong. 

NOTE: For a complete example including the considerations discussed above, see 

the examples/USER/flow_gauss directory. 

NOTE: Only the flux of the atoms in group-ID will be conserved. If the 

velocities of the group-ID atoms are coupled to the velocities of other atoms 

in the simulation, the flux will not be conserved. For example, in a 

simulation with fluid atoms and harmonically constrained wall atoms, if a 

single thermostat is applied to group {all}, the fluid atom velocities will be 

coupled to the wall atom velocities, and the flux will not be conserved. This 

issue can be avoided by thermostatting the fluid and wall groups separately.

Adding an acceleration to atoms does work on the system. This added energy 

 can be optionally subtracted from the potential energy for the thermodynamic 

output (see below) to check that the timestep is small enough to conserve 

energy. Since the applied acceleration is fluctuating in time, the work cannot 

be computed from a potential. As a result, computing the work is slightly more 

computationally expensive than usual, so it is not performed by default. To 

invoke the work calculation, use the {energy} keyword. The 

"fix_modify"_fix_modify.html {energy} option also invokes the work 

calculation, and overrides an {energy no} setting here. If neither {energy yes}

or {fix_modify energy yes} are set, the global scalar computed by the fix 

will return zero.

NOTE: In order to check energy conservation, any other fixes that do work on 

the system must have {fix_modify energy yes} set as well. This includes 

thermostat fixes and any constraints that hold the positions of wall atoms 

fixed, such as "fix spring/self"_fix_spring_self.html.

:line

[Restart, fix_modify, output, run start/stop, minimize info:]

This fix is part of the USER-MISC package.  It is only enabled if

LAMMPS was built with that package.  See the "Making

LAMMPS"_Section_start.html#start_3 section for more info.

No information about this fix is written to "binary restart

files"_restart.html.

The "fix_modify"_fix_modify.html {energy} option is supported by this

fix to subtract the work done from the

system's potential energy as part of "thermodynamic

output"_thermo_style.html.

This fix computes a global scalar and a global 3-vector of forces,

which can be accessed by various "output

commands"_Section_howto.html#howto_15.  The scalar is the negative of the 

work done on the system, see above discussion.  The vector is the total force 

that this fix applied to the group of atoms on the current timestep.

The scalar and vector values calculated by this fix are "extensive".

No parameter of this fix can be used with the {start/stop} keywords of

the "run"_run.html command.

[Restrictions:] none

[Related commands:]

"fix addforce"_fix_addforce.html, "compute temp/profile"_compute_temp_profile.html, "velocity"_velocity.html

[Default:]

The option default for the {energy} keyword is energy = no.

:line

:link(Strong)

[(Strong)] Strong and Eaves, J. Phys. Chem. Lett. 7, 1907 (2016).

:link(Evans)

[(Evans)] Evans and Morriss, Phys. Rev. Lett. 56, 2172 (1986).

:link(Zhu)

[(Zhu)] Zhu, Tajkhorshid, and Schulten, Biophys. J. 83, 154 (2002).

:line

The input script in.GD is an example simulation using Gaussian dynamics (GD).

The simulation is of a simple 2d Lennard-Jones fluid flowing through a pipe.

For details see online LAMMPS documentation and 

Strong and Eaves, J. Phys. Chem. Lett. 7(10) 2016, p. 1907.

Note that the run times and box size are chosen to allow a fast example run. 

They are not adequate for a real simulation. 

The script has the following parts:

1) initialize variables

	These can be modified to customize the simulation. Note that if the 

	pipe dimensions L or d are changed, the geometry should be checked 

	by visualizing the coordinates in all.init.lammpstrj. 

2) create box

3) set up potential

4) create atoms

5) set up profile-unbiased thermostat (PUT)

	see Evans and Morriss, Phys. Rev. Lett. 56(20) 1986, p. 2172

	By default, this uses boxes which contain on average 8 molecules.

6) equilibrate without GD

7) initialize the center-of-mass velocity and run to achieve steady-state

	The system is initialized with a uniform velocity profile, which 

	relaxes over the course of the simulation.

8) collect data

	The data is output in several files:

	GD.out contains the force that GD applies, and the flux in the x- and 

		y- directions. The output Jx should be equal to the value of 

		J set in section 1, which is 0.1 by default.

	x_profiles contains the velocity, density, and pressure profiles in 

		the x-direction. The pressure profile is given by 

		(-1/2V)*(c_spa[1] + c_spa[2]), where V is the volume of a 

		slice. The pressure profile is computed with IK1, see 

		Todd, Evans, and Davis, Phys. Rev. E 52(2) 1995, p. 1627.

		Note that to compare with the pump method, or to 

		compute a pressure drop, you must correct this pressure 

		profile as described in Strong 2016 above.

	Vy_profile is the velocity profile inside the pipe along the 

		y-direction, u_x(y).

#LAMMPS input script

#in.GD

#see README for details

###############################################################################

#initialize variables

clear

#frequency for outputting info (timesteps)

variable	dump_rate 	equal 50

variable	thermo_rate	equal 10

#equilibration time (timesteps)

variable	equil		equal 1000

#stabilization time (timesteps to reach steady-state)

variable	stabil		equal 1000

#data collection time (timesteps)

variable	run		equal 2000

#length of pipe

variable	L		equal 30

#width of pipe

variable	d		equal 20

#flux (mass/sigma*tau)

variable	J		equal 0.1

#simulation box dimensions

variable	Lx		equal 100

variable	Ly		equal 40

#bulk fluid density

variable	dens		equal 0.8

#lattice spacing for wall atoms

variable	aWall		equal 1.0 #1.7472

#timestep

variable	ts		equal 0.001

#temperature

variable	T		equal 2.0

#thermostat damping constant

variable	tdamp		equal ${ts}*100

units 		lj

dimension	2

atom_style	atomic

###############################################################################

#create box

#create lattice with the spacing aWall

variable	rhoWall		equal ${aWall}^(-2)

lattice		sq ${rhoWall}

#modify input dimensions to be multiples of aWall

variable	L1 equal round($L/${aWall})*${aWall}

variable	d1 equal round($d/${aWall})*${aWall}

variable	Ly1 equal round(${Ly}/${aWall})*${aWall}

variable	Lx1 equal round(${Lx}/${aWall})*${aWall}

#create simulation box

variable	lx2 equal ${Lx1}/2

variable	ly2 equal ${Ly1}/2

region		simbox block -${lx2} ${lx2} -${ly2} ${ly2} 0 0.1 units box

create_box	2 simbox

#####################################################################

#set up potential

mass 		1 1.0  		#fluid atoms

mass 		2 1.0		#wall atoms

pair_style	lj/cut 2.5

pair_modify	shift yes 

pair_coeff	1 1 1.0 1.0 2.5

pair_coeff	1 2 1.0 1.0 1.12246 

pair_coeff	2 2 0.0 0.0 0.0

timestep 	${ts}

#####################################################################

#create atoms

#create wall atoms everywhere

create_atoms    2 box

#define region which is "walled off"

variable 	dhalf equal ${d1}/2

variable 	Lhalf equal ${L1}/2

region		walltop block -${Lhalf} ${Lhalf} ${dhalf} EDGE -0.1 0.1 &

		units box

region		wallbot block -${Lhalf} ${Lhalf} EDGE -${dhalf} -0.1 0.1 &

		units box

region 		outsidewall union 2 walltop wallbot side out

#remove wall atoms outside wall region

group 		outside region outsidewall

delete_atoms 	group outside

#remove wall atoms that aren't on edge of wall region

variable	x1 equal ${Lhalf}-${aWall}

variable	y1 equal ${dhalf}+${aWall}

region 		insideTop block -${x1} ${x1} ${y1} EDGE -0.1 0.1 units box

region 		insideBot block -${x1} ${x1} EDGE -${y1} -0.1 0.1 units box

region		insideWall union 2 insideTop insideBot

group 		insideWall region insideWall

delete_atoms	group insideWall

#define new lattice, to give correct fluid density

#y lattice const must be a multiple of aWall

variable	atrue equal ${dens}^(-1/2)

variable	ay equal round(${atrue}/${aWall})*${aWall}

#choose x lattice const to give correct density

variable	ax equal (${ay}*${dens})^(-1)

#change Lx to be multiple of ax

variable	Lx1 equal round(${Lx}/${ax})*${ax}

variable	lx2 equal ${Lx1}/2

change_box 	all x final -${lx2} ${lx2} units box

#define new lattice

lattice 	custom ${dens} &

		a1 ${ax} 0.0 0.0 a2 0.0 ${ay} 0.0 a3 0.0 0.0 1.0 &

		basis 0.0 0.0 0.0

#fill in rest of box with bulk particles

variable	delta equal 0.001

variable	Ldelt equal ${Lhalf}+${delta}

variable	dDelt equal ${dhalf}-${delta}

region		left block EDGE -${Ldelt} EDGE EDGE -0.1 0.1 units box

region 		right block ${Ldelt} EDGE EDGE EDGE -0.1 0.1 units box

region		pipe block -${Ldelt} ${Ldelt} -${dDelt} ${dDelt} -0.1 0.1 &

		units box

region 		bulk union 3 left pipe right 

create_atoms	1 region bulk

group 		bulk type 1

group		wall type 2

#remove atoms that are too close to wall

delete_atoms    overlap 0.9 bulk wall

neighbor	0.3 bin

neigh_modify	delay 0 every 1 check yes

neigh_modify    exclude group wall wall

velocity 	bulk create $T 78915 dist gaussian rot yes mom yes loop geom

#####################################################################

#set up PUT

#see Evans and Morriss, Phys. Rev. Lett. 56(20) 1986, p. 2172

#average number of particles per box, Evans and Morriss used 2.0

variable	NperBox equal 8.0

#calculate box sizes

variable	boxSide equal sqrt(${NperBox}/${dens})

variable	nX equal round(lx/${boxSide})

variable	nY equal round(ly/${boxSide})

variable	dX equal lx/${nX}

variable	dY equal ly/${nY}

#temperature of fluid (excluding wall)

compute 	myT bulk temp

#profile-unbiased temperature of fluid

compute 	myTp bulk temp/profile 1 1 0 xy ${nX} ${nY} 

#thermo setup

thermo		${thermo_rate}	

thermo_style 	custom step c_myT c_myTp etotal press 

#dump initial configuration

dump 		55 all custom 1 all.init.lammpstrj id type x y z vx vy vz

dump 		56 wall custom 1 wall.init.lammpstrj id type x y z

dump_modify   	55 sort id

dump_modify   	56 sort id

run 		0

undump 		55

undump 		56

#####################################################################

#equilibrate without GD

fix 		nvt bulk nvt temp $T $T ${tdamp}

fix_modify	nvt temp myTp

fix		2 bulk enforce2d

run		${equil}

#####################################################################

#initialize the COM velocity and run to achieve steady-state

#calculate velocity to add: V=J/rho_total

variable	Vadd		equal $J*lx*ly/count(bulk)

#first remove any COM velocity, then add back the streaming velocity

velocity        bulk zero linear

velocity 	bulk set ${Vadd} 0.0 0.0 units box sum yes mom no

fix		GD bulk flow/gauss 1 0 0 #energy yes

#fix_modify	GD energy yes

run		${stabil}

#####################################################################

#collect data

#print the applied force and total flux to ensure conservation of Jx

variable	Fapp equal f_GD[1]

compute 	vxBulk bulk reduce sum vx

compute 	vyBulk bulk reduce sum vy

variable        invVol equal 1.0/(lx*ly)

variable	jx equal c_vxBulk*${invVol}

variable	jy equal c_vyBulk*${invVol}

variable	curr_step equal step

fix		print_vCOM all print ${dump_rate} &

		"${curr_step} ${Fapp} ${jx} ${jy}" file GD.out screen no &

		title "timestep Fapp Jx Jy"

#compute IK1 pressure profile 

#see Todd, Evans, and Davis, Phys. Rev. E 52(2) 1995, p. 1627

#use profile-unbiased temperature to remove the streaming velocity

#from the kinetic part of the pressure

compute		spa bulk stress/atom myTp

#for the pressure profile, use the same grid as the PUT

compute		chunkX bulk chunk/atom bin/1d x lower ${dX} units box

#output pressure profile and other profiles

#the pressure profile is (-1/2V)*(c_spa[1] + c_spa[2]), where

#V is the volume of a slice

fix		profiles bulk ave/chunk 1 1 ${dump_rate} chunkX &

		vx density/mass  c_spa[1] c_spa[2] &

		file x_profiles ave running overwrite

#compute velocity profile across the pipe with a finer grid

variable	dYnew equal ${dY}/10

compute		chunkY bulk chunk/atom bin/1d y center ${dYnew} units box & 

		region pipe

fix		velYprof bulk ave/chunk 1 1 ${dump_rate} chunkY &

		vx file Vy_profile ave running overwrite

#full trajectory

dump 		7 bulk custom ${dump_rate} bulk.lammpstrj &

		id type x y z 

dump_modify	7 sort id

run 		${run}

dihedral_style table, Andrew Jewett, jewett.aij@gmail.com, 10 Jan 12

fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11

fix ave/correlate/long, Jorge Ramirez (UPM Madrid), jorge.ramirez at upm.es, 21 Oct 2015

fix flow/gauss, Joel D. Eaves (CU Boulder), Joel.Eaves@Colorado.edu, 23 Aug 2016

fix gle, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014

fix imd, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 2009

fix ipi, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014

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

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

#include "atom.h"

#include "force.h"

#include "group.h"

#include "comm.h"

#include "update.h"

#include "domain.h"

#include "error.h"

#include "citeme.h"

using namespace LAMMPS_NS;

using namespace FixConst;

static const char cite_flow_gauss[] =

  "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";

FixFlowGauss::FixFlowGauss(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)

{

  if (lmp->citeme) lmp->citeme->add(cite_flow_gauss);

  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 flow/gauss 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 FixFlowGauss::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 FixFlowGauss::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 FixFlowGauss::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 FixFlowGauss::compute_scalar()

{

  return -pe_tot*dt;

}

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

   return components of applied force

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

double FixFlowGauss::compute_vector(int n)

{

  return -f_tot[n];

}