diff --git a/doc/compute_voronoi_atom.txt b/doc/compute_voronoi_atom.txt
index 73dc80e92..b97106aa0 100644
--- a/doc/compute_voronoi_atom.txt
+++ b/doc/compute_voronoi_atom.txt
@@ -1,183 +1,213 @@
"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
compute voronoi/atom command :h3
[Syntax:]
compute ID group-ID voronoi/atom keyword arg ... :pre
ID, group-ID are documented in "compute"_compute.html command :ulb,l
voronoi/atom = style name of this compute command :l
zero or more keyword/value pairs may be appended :l
keyword = {only_group} or {surface} or {radius} or {edge_histo} or {edge_threshold} or {face_threshold} :l
{only_group} = no arg
{occupation} = no arg
{surface} arg = sgroup-ID
sgroup-ID = compute the dividing surface between group-ID and sgroup-ID
this keyword adds a third column to the compute output
{radius} arg = v_r
v_r = radius atom style variable for a poly-disperse Voronoi tessellation
{edge_histo} arg = maxedge
maxedge = maximum number of Voronoi cell edges to be accounted in the histogram
{edge_threshold} arg = minlength
minlength = minimum length for an edge to be counted
{face_threshold} arg = minarea
- minarea = minimum area for a face to be counted :pre
+ minarea = minimum area for a face to be counted
+ {neighbors} value = {yes} or {no} = store list of all neighbors or no :pre
:ule
[Examples:]
compute 1 all voronoi/atom
compute 2 precipitate voronoi/atom surface matrix
compute 3b precipitate voronoi/atom radius v_r
compute 4 solute voronoi/atom only_group :pre
compute 5 defects voronoi/atom occupation :pre
+compute 6 all voronoi/atom neighbors yes
[Description:]
Define a computation that calculates the Voronoi tessellation of the
atoms in the simulation box. The tessellation is calculated using all
atoms in the simulation, but non-zero values are only stored for atoms
in the group.
By default two quantities per atom are calculated by this compute.
The first is the volume of the Voronoi cell around each atom. Any
point in an atom's Voronoi cell is closer to that atom than any other.
The second is the number of faces of the Voronoi cell. This is
equal to the number of nearest neighbors of the central atom,
plus any exterior faces (see note below).
:line
If the {only_group} keyword is specified the tessellation is performed
only with respect to the atoms contained in the compute group. This is
equivalent to deleting all atoms not contained in the group prior to
evaluating the tessellation.
If the {surface} keyword is specified a third quantity per atom is
computed: the Voronoi cell surface of the given atom. {surface} takes
a group ID as an argument. If a group other than {all} is specified,
only the Voronoi cell facets facing a neighbor atom from the specified
group are counted towards the surface area.
In the example above, a precipitate embedded in a matrix, only atoms
at the surface of the precipitate will have non-zero surface area, and
only the outward facing facets of the Voronoi cells are counted (the
hull of the precipitate). The total surface area of the precipitate
can be obtained by running a "reduce sum" compute on c_2\[3\]
If the {radius} keyword is specified with an atom style variable as
the argument, a poly-disperse Voronoi tessellation is
performed. Examples for radius variables are
variable r1 atom (type==1)*0.1+(type==2)*0.4
compute radius all property/atom radius
variable r2 atom c_radius :pre
Here v_r1 specifies a per-type radius of 0.1 units for type 1 atoms
and 0.4 units for type 2 atoms, and v_r2 accesses the radius property
present in atom_style sphere for granular models.
The {edge_histo} keyword activates the compilation of a histogram of
number of edges on the faces of the Voronoi cells in the compute
-group. The argument maxedge of the this keyword is the largest number
+group. The argument {maxedge} of the this keyword is the largest number
of edges on a single Voronoi cell face expected to occur in the
sample. This keyword adds the generation of a global vector with
-maxedge+1 entries. The last entry in the vector contains the number of
-faces with with more than maxedge edges. Since the polygon with the
+{maxedge}+1 entries. The last entry in the vector contains the number of
+faces with with more than {maxedge} edges. Since the polygon with the
smallest amount of edges is a triangle, entries 1 and 2 of the vector
will always be zero.
The {edge_threshold} and {face_threshold} keywords allow the
suppression of edges below a given minimum length and faces below a
given minimum area. Ultra short edges and ultra small faces can occur
as artifacts of the Voronoi tessellation. These keywords will affect
the neighbor count and edge histogram outputs.
If the {occupation} keyword is specified the tessellation is only
performed for the first invocation of the compute and then stored.
For all following invocations of the compute the number of atoms in
each Voronoi cell in the stored tessellation is counted. In this mode
the compute returns a per-atom array with 2 columns. The first column
is the number of atoms currently in the Voronoi volume defined by this
atom at the time of the first invocation of the compute (note that the
atom may have moved significantly). The second column contains the
total number of atoms sharing the Voronoi cell of the stored
tessellation at the location of the current atom. Numbers in column
one can be any positive integer including zero, while column two
values will always be greater than zero. Column one data can be used
to locate vacancies (the coordinates are given by the atom coordinates
at the time step when the compute was first invoked), while column two
data can be used to identify interstitial atoms.
+If the {neighbors} value is set to yes, then
+this compute creates a local array with 3 columns. There
+is one row for each face of each Voronoi cell. The
+3 columns are the atom ID of the atom that owns the cell,
+the atom ID of the atom in the neighboring cell
+(or zero if the face is external), and the area of the face.
+The array can be accessed by any command that
+uses local values from a compute as input. See "this
+section"_Section_howto.html#howto_15 for an overview of LAMMPS output
+options. More specifically, the array can be accessed by a
+"dump local"_dump.html command to write a file containing
+all the Voronoi neighbors in a system:
+
+compute 6 all voronoi/atom neighbors yes
+dump d2 all local 1 dump.neighbors index c_6\[1\] c_6\[2\] c_6\[3\] :pre
+
+If the {face_threshold} keyword is used, then only faces
+with areas greater than the threshold are stored.
+
:line
The Voronoi calculation is performed by the freely available "Voro++
package"_voronoi, written by Chris Rycroft at UC Berkeley and LBL,
which must be installed on your system when building LAMMPS for use
with this compute. See instructions on obtaining and installing the
Voro++ software in the src/VORONOI/README file.
:link(voronoi,http://math.lbl.gov/voro++)
NOTE: The calculation of Voronoi volumes is performed by each
processor for the atoms it owns, and includes the effect of ghost
atoms stored by the processor. This assumes that the Voronoi cells of
owned atoms are not affected by atoms beyond the ghost atom cut-off
distance. This is usually a good assumption for liquid and solid
systems, but may lead to underestimation of Voronoi volumes in low
density systems. By default, the set of ghost atoms stored by each
processor is determined by the cutoff used for
"pair_style"_pair_style.html interactions. The cutoff can be set
explicitly via the "comm_modify cutoff"_comm_modify.html command.
The Voronoi cells for atoms adjacent to empty regions will extend
into those regions up to the communication cutoff in x, y, or z.
In that situation, an exterior
face is created at the cutoff distance normal to the x, y, or z
direction. For triclinic systems, the exterior face is
parallel to the corresponding reciprocal lattice vector.
NOTE: The Voro++ package performs its calculation in 3d. This will
still work for a 2d LAMMPS simulation, provided all the atoms have
the same z coordinate. The Voronoi cell of each atom will be
a columnar polyhedron with constant cross-sectional area
along the z direction and two exterior faces at the top and bottom of the
simulation box. If the atoms do not all
have the same z coordinate, then the columnar cells will be
accordingly distorted. The cross-sectional area of each Voronoi
cell can be obtained by dividing its volume by the z extent
of the simulation box. Note
that you define the z extent of the simulation box for 2d simulations
when using the "create_box"_create_box.html or
"read_data"_read_data.html commands.
[Output info:]
This compute calculates a per-atom array with 2 columns. In regular
dynamic tessellation mode the first column is the Voronoi volume, the
second is the neighbor count, as described above (read above for the
output data in case the {occupation} keyword is specified).
These values can be accessed by any command that
uses per-atom values from a compute as input. See "Section_howto
15"_Section_howto.html#howto_15 for an overview of LAMMPS output
options.
+If the {edge_histo} keyword is used, then this compute
+generates a global vector of length {maxedge}+1, containing
+a histogram of the number of edges per face.
+
+If the {neighbors} value is set to yes, then
+this compute calculates a local array with 3 columns. There
+is one row for each face of each Voronoi cell.
+
The Voronoi cell volume will be in distance "units"_units.html cubed.
+The Voronoi face area will be in distance "units"_units.html squared.
[Restrictions:]
This compute is part of the VORONOI 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.
[Related commands:]
-"dump custom"_dump.html
+"dump custom"_dump.html, "dump local"_dump.html
[Default:] none
diff --git a/examples/voronoi/in.voronoi.2d b/examples/voronoi/in.voronoi.2d
index 72f973b6f..64059ca5e 100644
--- a/examples/voronoi/in.voronoi.2d
+++ b/examples/voronoi/in.voronoi.2d
@@ -1,54 +1,56 @@
# Test volume definitions for 2d and finite systems
variable rcut equal 10.0
variable rskin equal 2.0
variable rcomm equal 20.0
variable len equal 4.0
variable lenz equal 10.0
dimension 2
units metal
boundary p p p
#lattice sq 1.0 origin 0.5 0.5 0.0
lattice hex 1.0 origin 0.5 0.5 0.0
atom_style atomic
region box block 0 ${len} 0 ${len} 0.0 ${lenz}
region atoms block 0 ${len} 0 ${len} 0.0 0.0
create_box 1 box
create_atoms 1 region atoms
mass 1 1.0
pair_style lj/cut ${rcut}
pair_coeff 1 1 0.0 1.0
neighbor ${rskin} nsq
neigh_modify delay 10
# set the minimum communication cut-off
comm_modify cutoff ${rcomm}
thermo 1
#
# TEST 1: Volume check for 2d bulk system
#
-compute v1 all voronoi/atom
+compute v1 all voronoi/atom neighbors yes
dump d1 all custom 1 dump.voro id type x y z c_v1[1] c_v1[2]
+dump d2 all local 1 dump.neighbors index c_v1[1] c_v1[2] c_v1[3]
+
compute volvor all reduce sum c_v1[1]
variable volsys equal lz*lx*ly
variable err equal c_volvor-v_volsys
thermo_style custom c_volvor v_volsys vol v_err
run 0
#
# TEST 2: Volume check for 2d finite system
# add margins in x and y directions
#
change_box all boundary f f p
run 0
diff --git a/src/VORONOI/compute_voronoi_atom.cpp b/src/VORONOI/compute_voronoi_atom.cpp
index 052260895..305a98716 100644
--- a/src/VORONOI/compute_voronoi_atom.cpp
+++ b/src/VORONOI/compute_voronoi_atom.cpp
@@ -1,571 +1,624 @@
/* ----------------------------------------------------------------------
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: Daniel Schwen
------------------------------------------------------------------------- */
#include <mpi.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "compute_voronoi_atom.h"
#include "atom.h"
#include "group.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
#include "comm.h"
#include "variable.h"
#include "input.h"
#include "force.h"
#include <vector>
using namespace LAMMPS_NS;
using namespace voro;
+#define FACESDELTA 10000
+
/* ---------------------------------------------------------------------- */
ComputeVoronoi::ComputeVoronoi(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
int sgroup;
size_peratom_cols = 2;
peratom_flag = 1;
comm_forward = 1;
+ faces_flag = 0;
surface = VOROSURF_NONE;
maxedge = 0;
fthresh = ethresh = 0.0;
radstr = NULL;
onlyGroup = false;
occupation = false;
con_mono = NULL;
con_poly = NULL;
tags = NULL;
occvec = sendocc = lroot = lnext = NULL;
+ faces = NULL;
int iarg = 3;
while ( iarg<narg ) {
if (strcmp(arg[iarg], "occupation") == 0) {
occupation = true;
iarg++;
}
else if (strcmp(arg[iarg], "only_group") == 0) {
onlyGroup = true;
iarg++;
}
else if (strcmp(arg[iarg], "radius") == 0) {
if (iarg + 2 > narg || strstr(arg[iarg+1],"v_") != arg[iarg+1] )
error->all(FLERR,"Illegal compute voronoi/atom command");
int n = strlen(&arg[iarg+1][2]) + 1;
radstr = new char[n];
strcpy(radstr,&arg[iarg+1][2]);
iarg += 2;
}
else if (strcmp(arg[iarg], "surface") == 0) {
if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
// group all is a special case where we just skip group testing
if(strcmp(arg[iarg+1], "all") == 0) {
surface = VOROSURF_ALL;
} else {
sgroup = group->find(arg[iarg+1]);
if (sgroup == -1) error->all(FLERR,"Could not find compute/voronoi surface group ID");
sgroupbit = group->bitmask[sgroup];
surface = VOROSURF_GROUP;
}
size_peratom_cols = 3;
iarg += 2;
} else if (strcmp(arg[iarg], "edge_histo") == 0) {
if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
maxedge = force->inumeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg], "face_threshold") == 0) {
if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
fthresh = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg], "edge_threshold") == 0) {
if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
ethresh = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
+ } else if (strcmp(arg[iarg], "neighbors") == 0) {
+ if (iarg + 2 > narg) error->all(FLERR,"Illegal compute voronoi/atom command");
+ if (strcmp(arg[iarg+1],"yes") == 0) faces_flag = 1;
+ else if (strcmp(arg[iarg+1],"no") == 0) faces_flag = 0;
+ else error->all(FLERR,"Illegal compute voronoi/atom command");
+ iarg += 2;
}
else error->all(FLERR,"Illegal compute voronoi/atom command");
}
if (occupation && ( surface!=VOROSURF_NONE || maxedge>0 ) )
error->all(FLERR,"Illegal compute voronoi/atom command (occupation and (surface or edges))");
nmax = rmax = 0;
edge = rfield = sendvector = NULL;
voro = NULL;
if ( maxedge > 0 ) {
vector_flag = 1;
size_vector = maxedge+1;
memory->create(edge,maxedge+1,"voronoi/atom:edge");
memory->create(sendvector,maxedge+1,"voronoi/atom:sendvector");
vector = edge;
}
+
+ // store local face data: i, j, area
+
+ if (faces_flag) {
+ local_flag = 1;
+ size_local_cols = 3;
+ nfacesmax = 0;
+ }
}
/* ---------------------------------------------------------------------- */
ComputeVoronoi::~ComputeVoronoi()
{
memory->destroy(edge);
memory->destroy(rfield);
memory->destroy(sendvector);
memory->destroy(voro);
delete[] radstr;
// voro++ container classes
delete con_mono;
delete con_poly;
// occupation analysis stuff
memory->destroy(lroot);
memory->destroy(lnext);
memory->destroy(occvec);
#ifdef NOTINPLACE
memory->destroy(sendocc);
#endif
memory->destroy(tags);
+ memory->destroy(faces);
}
/* ---------------------------------------------------------------------- */
void ComputeVoronoi::init()
{
}
/* ----------------------------------------------------------------------
gather compute vector data from other nodes
------------------------------------------------------------------------- */
void ComputeVoronoi::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow per atom array if necessary
int nlocal = atom->nlocal;
if (nlocal > nmax) {
memory->destroy(voro);
nmax = atom->nmax;
memory->create(voro,nmax,size_peratom_cols,"voronoi/atom:voro");
array_atom = voro;
}
// decide between occupation or per-frame tesselation modes
if (occupation) {
// build cells only once
int i, nall = nlocal + atom->nghost;
if (con_mono==NULL && con_poly==NULL) {
// generate the voronoi cell network for the initial structure
buildCells();
// save tags of atoms (i.e. of each voronoi cell)
memory->create(tags,nall,"voronoi/atom:tags");
for (i=0; i<nall; i++) tags[i] = atom->tag[i];
// linked list structure for cell occupation count on the atoms
oldnall= nall;
memory->create(lroot,nall,"voronoi/atom:lroot"); // point to first atom index in cell (or -1 for empty cell)
lnext = NULL;
lmax = 0;
// build the occupation buffer
oldnatoms = atom->natoms;
memory->create(occvec,oldnatoms,"voronoi/atom:occvec");
#ifdef NOTINPLACE
memory->create(sendocc,oldnatoms,"voronoi/atom:sendocc");
#endif
}
// get the occupation of each original voronoi cell
checkOccupation();
} else {
// build cells for each output
buildCells();
loopCells();
}
}
void ComputeVoronoi::buildCells()
{
int i;
const double e = 0.01;
int nlocal = atom->nlocal;
int dim = domain->dimension;
// in the onlyGroup mode we are not setting values for all atoms later in the voro loop
// initialize everything to zero here
if (onlyGroup) {
if (surface == VOROSURF_NONE)
for (i = 0; i < nlocal; i++) voro[i][0] = voro[i][1] = 0.0;
else
for (i = 0; i < nlocal; i++) voro[i][0] = voro[i][1] = voro[i][2] = 0.0;
}
double *sublo = domain->sublo, *sublo_lamda = domain->sublo_lamda, *boxlo = domain->boxlo;
double *subhi = domain->subhi, *subhi_lamda = domain->subhi_lamda;
double *cut = comm->cutghost;
double sublo_bound[3], subhi_bound[3], cut_bound[3];
double **x = atom->x;
// setup bounds for voro++ domain for orthogonal and triclinic simulation boxes
if( domain->triclinic ) {
// triclinic box: embed parallelepiped into orthogonal voro++ domain
// cutghost is in lamda coordinates for triclinic boxes, use subxx_lamda
double *h = domain->h;
for( i=0; i<3; ++i ) {
sublo_bound[i] = sublo[i]-cut[i]-e;
subhi_bound[i] = subhi[i]+cut[i]+e;
if (domain->periodicity[i]==0) {
sublo_bound[i] = MAX(sublo_bound[i],0.0);
subhi_bound[i] = MIN(subhi_bound[i],1.0);
}
}
if (dim == 2) {
sublo_bound[2] = 0.0;
subhi_bound[2] = 1.0;
}
sublo_bound[0] = h[0]*sublo_bound[0] + h[5]*sublo_bound[1] + h[4]*sublo_bound[2] + boxlo[0];
sublo_bound[1] = h[1]*sublo_bound[1] + h[3]*sublo_bound[2] + boxlo[1];
sublo_bound[2] = h[2]*sublo_bound[2] + boxlo[2];
subhi_bound[0] = h[0]*subhi_bound[0] + h[5]*subhi_bound[1] + h[4]*subhi_bound[2] + boxlo[0];
subhi_bound[1] = h[1]*subhi_bound[1] + h[3]*subhi_bound[2] + boxlo[1];
subhi_bound[2] = h[2]*subhi_bound[2] + boxlo[2];
} else {
// orthogonal box
for( i=0; i<3; ++i ) {
sublo_bound[i] = sublo[i]-cut[i]-e;
subhi_bound[i] = subhi[i]+cut[i]+e;
if (domain->periodicity[i]==0) {
sublo_bound[i] = MAX(sublo_bound[i],domain->boxlo[i]);
subhi_bound[i] = MIN(subhi_bound[i],domain->boxhi[i]);
}
}
if (dim == 2) {
sublo_bound[2] = sublo[2];
subhi_bound[2] = subhi[2];
}
}
// n = # of voro++ spatial hash cells (with approximately cubic cells)
int nall = nlocal + atom->nghost;
double n[3], V;
for( i=0; i<3; ++i ) n[i] = subhi_bound[i] - sublo_bound[i];
V = n[0]*n[1]*n[2];
for( i=0; i<3; ++i ) {
n[i] = round( n[i]*pow( double(nall)/(V*8.0), 0.333333 ) );
n[i] = n[i]==0 ? 1 : n[i];
}
// clear edge statistics
for (i = 0; i < maxedge; ++i) edge[i]=0;
// initialize voro++ container
// preallocates 8 atoms per cell
// voro++ allocates more memory if needed
int *mask = atom->mask;
if (radstr) {
// check and fetch atom style variable data
int radvar = input->variable->find(radstr);
if (radvar < 0)
error->all(FLERR,"Variable name for voronoi radius does not exist");
if (!input->variable->atomstyle(radvar))
error->all(FLERR,"Variable for voronoi radius is not atom style");
// prepare destination buffer for variable evaluation
if (nlocal > rmax) {
memory->destroy(rfield);
rmax = atom->nmax;
memory->create(rfield,rmax,"voronoi/atom:rfield");
}
// compute atom style radius variable
input->variable->compute_atom(radvar,0,rfield,1,0);
// communicate values to ghost atoms of neighboring nodes
comm->forward_comm_compute(this);
// polydisperse voro++ container
delete con_poly;
con_poly = new container_poly(sublo_bound[0],
subhi_bound[0],
sublo_bound[1],
subhi_bound[1],
sublo_bound[2],
subhi_bound[2],
int(n[0]),int(n[1]),int(n[2]),
false,false,false,8);
// pass coordinates for local and ghost atoms to voro++
for (i = 0; i < nall; i++) {
if( !onlyGroup || (mask[i] & groupbit) )
con_poly->put(i,x[i][0],x[i][1],x[i][2],rfield[i]);
}
} else {
// monodisperse voro++ container
delete con_mono;
con_mono = new container(sublo_bound[0],
subhi_bound[0],
sublo_bound[1],
subhi_bound[1],
sublo_bound[2],
subhi_bound[2],
int(n[0]),int(n[1]),int(n[2]),
false,false,false,8);
// pass coordinates for local and ghost atoms to voro++
for (i = 0; i < nall; i++)
if( !onlyGroup || (mask[i] & groupbit) )
con_mono->put(i,x[i][0],x[i][1],x[i][2]);
}
}
void ComputeVoronoi::checkOccupation()
{
// clear occupation vector
memset(occvec, 0, oldnatoms*sizeof(*occvec));
int i, j, k,
nlocal = atom->nlocal,
nall = atom->nghost + nlocal;
double rx, ry, rz,
**x = atom->x;
// prepare destination buffer for variable evaluation
if (nall > lmax) {
memory->destroy(lnext);
lmax = atom->nmax;
memory->create(lnext,lmax,"voronoi/atom:lnext");
}
// clear lroot
for (i=0; i<oldnall; ++i) lroot[i] = -1;
// clear lnext
for (i=0; i<nall; ++i) lnext[i] = -1;
// loop over all local atoms and find out in which of the local first frame voronoi cells the are in
// (need to loop over ghosts, too, to get correct occupation numbers for the second column)
for (i=0; i<nall; ++i) {
// again: find_voronoi_cell() should be in the common base class. Why it is not, I don't know. Ask the voro++ author.
if (( radstr && con_poly->find_voronoi_cell(x[i][0], x[i][1], x[i][2], rx, ry, rz, k)) ||
( !radstr && con_mono->find_voronoi_cell(x[i][0], x[i][1], x[i][2], rx, ry, rz, k) )) {
// increase occupation count of this particular cell
// only for local atoms, as we do an MPI reduce sum later
if (i<nlocal) occvec[tags[k]-1]++;
// add this atom to the linked list of cell j
if (lroot[k]<0)
lroot[k]=i;
else {
j = lroot[k];
while (lnext[j]>=0) j=lnext[j];
lnext[j] = i;
}
}
}
// MPI sum occupation
#ifdef NOTINPLACE
memcpy(sendocc, occvec, oldnatoms*sizeof(*occvec));
MPI_Allreduce(sendocc, occvec, oldnatoms, MPI_INT, MPI_SUM, world);
#else
MPI_Allreduce(MPI_IN_PLACE, occvec, oldnatoms, MPI_INT, MPI_SUM, world);
#endif
// determine the total number of atoms in this atom's currently occupied cell
int c;
for (i=0; i<oldnall; i++) { // loop over lroot (old voronoi cells)
// count
c = 0;
j = lroot[i];
while (j>=0) {
c++;
j = lnext[j];
}
// set
j = lroot[i];
while (j>=0) {
voro[j][1] = c;
j = lnext[j];
}
}
// cherry pick currently owned atoms
for (i=0; i<nlocal; i++) {
// set the new atom count in the atom's first frame voronoi cell
voro[i][0] = occvec[atom->tag[i]-1];
}
}
void ComputeVoronoi::loopCells()
{
// invoke voro++ and fetch results for owned atoms in group
voronoicell_neighbor c;
int i;
+ if (faces_flag) nfaces = 0;
if (radstr) {
c_loop_all cl(*con_poly);
if (cl.start()) do if (con_poly->compute_cell(c,cl)) {
i = cl.pid();
processCell(c,i);
} while (cl.inc());
} else {
c_loop_all cl(*con_mono);
if (cl.start()) do if (con_mono->compute_cell(c,cl)) {
i = cl.pid();
processCell(c,i);
} while (cl.inc());
}
+ if (faces_flag) size_local_rows = nfaces;
+
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
void ComputeVoronoi::processCell(voronoicell_neighbor &c, int i)
{
int j,k, *mask = atom->mask;
std::vector<int> neigh, norder, vlist;
std::vector<double> narea, vcell;
bool have_narea = false;
// zero out surface area if surface computation was requested
if (surface != VOROSURF_NONE && !onlyGroup) voro[i][2] = 0.0;
if (i < atom->nlocal && (mask[i] & groupbit)) {
// cell volume
voro[i][0] = c.volume();
// number of cell faces
c.neighbors(neigh);
int neighs = neigh.size();
-// // DEBUG CODE!!
-// // loop over all faces (neighbors)
-// for (j=0; j<neighs; ++j)
-// printf("neighbors %d %d %d\n",i,j,neigh[j]);
-// // END DEBUG CODE!!
-
if (fthresh > 0) {
// count only faces above area threshold
c.face_areas(narea);
have_narea = true;
voro[i][1] = 0.0;
for (j=0; j<narea.size(); ++j)
if (narea[j] > fthresh) voro[i][1] += 1.0;
} else {
// unthresholded face count
voro[i][1] = neighs;
}
// cell surface area
if (surface == VOROSURF_ALL) {
voro[i][2] = c.surface_area();
} else if (surface == VOROSURF_GROUP) {
if (!have_narea) c.face_areas(narea);
voro[i][2] = 0.0;
// each entry in neigh should correspond to an entry in narea
if (neighs != narea.size())
- error->all(FLERR,"Voro++ error: narea and neigh have a different size");
+ error->one(FLERR,"Voro++ error: narea and neigh have a different size");
// loop over all faces (neighbors) and check if they are in the surface group
for (j=0; j<neighs; ++j)
if (neigh[j] >= 0 && mask[neigh[j]] & sgroupbit)
voro[i][2] += narea[j];
}
// histogram of number of face edges
+
if (maxedge>0) {
if (ethresh > 0) {
// count only edges above length threshold
c.vertices(vcell);
c.face_vertices(vlist); // for each face: vertex count followed list of vertex indices (n_1,v1_1,v2_1,v3_1,..,vn_1,n_2,v2_1,...)
double dx, dy, dz, r2, t2 = ethresh*ethresh;
for( j=0; j<vlist.size(); j+=vlist[j]+1 ) {
int a, b, nedge = 0;
// vlist[j] contains number of vertex indices for the current face
for( k=0; k<vlist[j]; ++k ) {
a = vlist[j+1+k]; // first vertex in edge
b = vlist[j+1+(k+1)%vlist[j]]; // second vertex in edge (possible wrap around to first vertex in list)
dx = vcell[a*3] - vcell[b*3];
dy = vcell[a*3+1] - vcell[b*3+1];
dz = vcell[a*3+2] - vcell[b*3+2];
r2 = dx*dx+dy*dy+dz*dz;
if (r2 > t2) nedge++;
}
// counted edges above threshold, now put into the correct bin
if (nedge>0) {
if (nedge<=maxedge)
edge[nedge-1]++;
else
edge[maxedge]++;
}
}
} else {
// unthresholded edge counts
c.face_orders(norder);
for (j=0; j<voro[i][1]; ++j)
if (norder[j]>0) {
if (norder[j]<=maxedge)
edge[norder[j]-1]++;
else
edge[maxedge]++;
}
}
}
+
+ // store info for local faces
+
+ if (faces_flag) {
+ if (nfaces+voro[i][1] > nfacesmax) {
+ while (nfacesmax < nfaces+voro[i][1]) nfacesmax += FACESDELTA;
+ memory->grow(faces,nfacesmax,size_local_cols,"compute/voronoi/atom:faces");
+ array_local = faces;
+ }
+
+ if (!have_narea) c.face_areas(narea);
+
+ if (neighs != narea.size())
+ error->one(FLERR,"Voro++ error: narea and neigh have a different size");
+ tagint itag, jtag;
+ tagint *tag = atom->tag;
+ itag = tag[i];
+ for (j=0; j<neighs; ++j)
+ if (narea[j] > fthresh) {
+
+ // external faces assigned the tag 0
+
+ int jj = neigh[j];
+ if (jj >= 0) jtag = tag[jj];
+ else jtag = 0;
+
+ faces[nfaces][0] = itag;
+ faces[nfaces][1] = jtag;
+ faces[nfaces][2] = narea[j];
+ nfaces++;
+ }
+ }
+
+
} else if (i < atom->nlocal) voro[i][0] = voro[i][1] = 0.0;
}
double ComputeVoronoi::memory_usage()
{
double bytes = size_peratom_cols * nmax * sizeof(double);
+ // estimate based on average coordination of 12
+ if (faces_flag) bytes += 12 * size_local_cols * nmax * sizeof(double);
return bytes;
}
void ComputeVoronoi::compute_vector()
{
invoked_vector = update->ntimestep;
if( invoked_peratom < invoked_vector ) compute_peratom();
for( int i=0; i<size_vector; ++i ) sendvector[i] = edge[i];
MPI_Allreduce(sendvector,edge,size_vector,MPI_DOUBLE,MPI_SUM,world);
}
/* ---------------------------------------------------------------------- */
int ComputeVoronoi::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,m=0;
for (i = 0; i < n; ++i) buf[m++] = rfield[list[i]];
return m;
}
/* ---------------------------------------------------------------------- */
void ComputeVoronoi::unpack_forward_comm(int n, int first, double *buf)
{
int i,last,m=0;
last = first + n;
for (i = first; i < last; ++i) rfield[i] = buf[m++];
}
diff --git a/src/VORONOI/compute_voronoi_atom.h b/src/VORONOI/compute_voronoi_atom.h
index ca58eca57..086b0e293 100644
--- a/src/VORONOI/compute_voronoi_atom.h
+++ b/src/VORONOI/compute_voronoi_atom.h
@@ -1,94 +1,96 @@
/* -*- 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 COMPUTE_CLASS
ComputeStyle(voronoi/atom,ComputeVoronoi)
#else
#ifndef LMP_COMPUTE_VORONOI_H
#define LMP_COMPUTE_VORONOI_H
#include "compute.h"
#include "voro++.hh"
namespace LAMMPS_NS {
class ComputeVoronoi : public Compute {
public:
ComputeVoronoi(class LAMMPS *, int, char **);
~ComputeVoronoi();
void init();
void compute_peratom();
void compute_vector();
double memory_usage();
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
private:
voro::container *con_mono;
voro::container_poly *con_poly;
void buildCells();
void checkOccupation();
void loopCells();
void processCell(voro::voronoicell_neighbor&, int);
int nmax, rmax, maxedge, sgroupbit;
char *radstr;
double fthresh, ethresh;
double **voro;
double *edge, *sendvector, *rfield;
enum { VOROSURF_NONE, VOROSURF_ALL, VOROSURF_GROUP } surface;
bool onlyGroup, occupation;
tagint *tags;
int *occvec, *sendocc, *lroot, *lnext, lmax, oldnatoms, oldnall;
+ int faces_flag, nfaces, nfacesmax;
+ double **faces;
};
}
#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: Could not find compute/voronoi surface group ID
Self-explanatory.
E: Illegal compute voronoi/atom command (occupation and (surface or edges))
Self-explanatory.
E: Variable name for voronoi radius does not exist
Self-explanatory.
E: Variable for voronoi radius is not atom style
Self-explanatory.
E: Voro++ error: narea and neigh have a different size
This error is returned by the Voro++ library.
*/
diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp
index c71ffbc3f..cf3161e8a 100644
--- a/src/fix_langevin.cpp
+++ b/src/fix_langevin.cpp
@@ -1,919 +1,916 @@
/* ----------------------------------------------------------------------
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: Carolyn Phillips (U Mich), reservoir energy tally
Aidan Thompson (SNL) GJF formulation
------------------------------------------------------------------------- */
#include <mpi.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "fix_langevin.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "comm.h"
#include "input.h"
#include "variable.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
#include "group.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{NOBIAS,BIAS};
enum{CONSTANT,EQUAL,ATOM};
#define SINERTIA 0.4 // moment of inertia prefactor for sphere
#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
/* ---------------------------------------------------------------------- */
FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) error->all(FLERR,"Illegal fix langevin command");
dynamic_group_allow = 1;
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
nevery = 1;
tstr = NULL;
if (strstr(arg[3],"v_") == arg[3]) {
int n = strlen(&arg[3][2]) + 1;
tstr = new char[n];
strcpy(tstr,&arg[3][2]);
} else {
t_start = force->numeric(FLERR,arg[3]);
t_target = t_start;
tstyle = CONSTANT;
}
t_stop = force->numeric(FLERR,arg[4]);
t_period = force->numeric(FLERR,arg[5]);
seed = force->inumeric(FLERR,arg[6]);
if (t_period <= 0.0) error->all(FLERR,"Fix langevin period must be > 0.0");
if (seed <= 0) error->all(FLERR,"Illegal fix langevin command");
// 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];
ratio = new double[atom->ntypes+1];
// optional args
for (int i = 1; i <= atom->ntypes; i++) ratio[i] = 1.0;
ascale = 0.0;
gjfflag = 0;
oflag = 0;
tallyflag = 0;
zeroflag = 0;
int iarg = 7;
while (iarg < narg) {
if (strcmp(arg[iarg],"angmom") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
if (strcmp(arg[iarg+1],"no") == 0) ascale = 0.0;
else ascale = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"gjf") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
if (strcmp(arg[iarg+1],"no") == 0) gjfflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) gjfflag = 1;
else error->all(FLERR,"Illegal fix langevin command");
iarg += 2;
} else if (strcmp(arg[iarg],"omega") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
if (strcmp(arg[iarg+1],"no") == 0) oflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) oflag = 1;
else error->all(FLERR,"Illegal fix langevin command");
iarg += 2;
} else if (strcmp(arg[iarg],"scale") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix langevin command");
int itype = force->inumeric(FLERR,arg[iarg+1]);
double scale = force->numeric(FLERR,arg[iarg+2]);
if (itype <= 0 || itype > atom->ntypes)
error->all(FLERR,"Illegal fix langevin command");
ratio[itype] = scale;
iarg += 3;
} else if (strcmp(arg[iarg],"tally") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
if (strcmp(arg[iarg+1],"no") == 0) tallyflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) tallyflag = 1;
else error->all(FLERR,"Illegal fix langevin command");
iarg += 2;
} else if (strcmp(arg[iarg],"zero") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
if (strcmp(arg[iarg+1],"no") == 0) zeroflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) zeroflag = 1;
else error->all(FLERR,"Illegal fix langevin command");
iarg += 2;
} else error->all(FLERR,"Illegal fix langevin command");
}
// set temperature = NULL, user can override via fix_modify if wants bias
id_temp = NULL;
temperature = NULL;
- // flangevin is unallocated until first call to setup()
- // compute_scalar checks for this and returns 0.0 if flangevin is NULL
-
energy = 0.0;
flangevin = NULL;
franprev = NULL;
tforce = NULL;
maxatom1 = maxatom2 = 0;
// Setup atom-based array for franprev
// register with Atom class
// No need to set peratom_flag
// as this data is for internal use only
if (gjfflag) {
nvalues = 3;
grow_arrays(atom->nmax);
atom->add_callback(0);
// initialize franprev to zero
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
franprev[i][0] = 0.0;
franprev[i][1] = 0.0;
franprev[i][2] = 0.0;
}
}
if (tallyflag && zeroflag && comm->me == 0)
error->warning(FLERR,"Energy tally does not account for 'zero yes'");
}
/* ---------------------------------------------------------------------- */
FixLangevin::~FixLangevin()
{
delete random;
delete [] tstr;
delete [] gfactor1;
delete [] gfactor2;
delete [] ratio;
delete [] id_temp;
memory->destroy(flangevin);
memory->destroy(tforce);
if (gjfflag) {
memory->destroy(franprev);
atom->delete_callback(id,0);
}
}
/* ---------------------------------------------------------------------- */
int FixLangevin::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= END_OF_STEP;
mask |= THERMO_ENERGY;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixLangevin::init()
{
if (oflag && !atom->sphere_flag)
error->all(FLERR,"Fix langevin omega requires atom style sphere");
if (ascale && !atom->ellipsoid_flag)
error->all(FLERR,"Fix langevin angmom requires atom style ellipsoid");
// check variable
if (tstr) {
tvar = input->variable->find(tstr);
if (tvar < 0)
error->all(FLERR,"Variable name for fix langevin does not exist");
if (input->variable->equalstyle(tvar)) tstyle = EQUAL;
else if (input->variable->atomstyle(tvar)) tstyle = ATOM;
else error->all(FLERR,"Variable for fix langevin is invalid style");
}
// if oflag or ascale set, check that all group particles are finite-size
if (oflag) {
double *radius = atom->radius;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (radius[i] == 0.0)
error->one(FLERR,"Fix langevin omega requires extended particles");
}
if (ascale) {
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec)
error->all(FLERR,"Fix langevin angmom requires atom style ellipsoid");
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (ellipsoid[i] < 0)
error->one(FLERR,"Fix langevin angmom requires extended particles");
}
// set force prefactors
if (!atom->rmass) {
for (int i = 1; i <= atom->ntypes; i++) {
gfactor1[i] = -atom->mass[i] / t_period / force->ftm2v;
gfactor2[i] = sqrt(atom->mass[i]) *
sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v;
gfactor1[i] *= 1.0/ratio[i];
gfactor2[i] *= 1.0/sqrt(ratio[i]);
}
}
if (temperature && temperature->tempbias) tbiasflag = BIAS;
else tbiasflag = NOBIAS;
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (gjfflag) gjffac = 1.0/(1.0+update->dt/2.0/t_period);
}
/* ---------------------------------------------------------------------- */
void FixLangevin::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 FixLangevin::post_force(int vflag)
{
double *rmass = atom->rmass;
// enumerate all 2^6 possibilities for template parameters
// this avoids testing them inside inner loop:
// TSTYLEATOM, GJF, TALLY, BIAS, RMASS, ZERO
#ifdef TEMPLATED_FIX_LANGEVIN
if (tstyle == ATOM)
if (gjfflag)
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,1,1,1,1,1>();
else post_force_templated<1,1,1,1,1,0>();
else
if (zeroflag) post_force_templated<1,1,1,1,0,1>();
else post_force_templated<1,1,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,1,1,0,1,1>();
else post_force_templated<1,1,1,0,1,0>();
else
if (zeroflag) post_force_templated<1,1,1,0,0,1>();
else post_force_templated<1,1,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,1,0,1,1,1>();
else post_force_templated<1,1,0,1,1,0>();
else
if (zeroflag) post_force_templated<1,1,0,1,0,1>();
else post_force_templated<1,1,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,1,0,0,1,1>();
else post_force_templated<1,1,0,0,1,0>();
else
if (zeroflag) post_force_templated<1,1,0,0,0,1>();
else post_force_templated<1,1,0,0,0,0>();
else
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,0,1,1,1,1>();
else post_force_templated<1,0,1,1,1,0>();
else
if (zeroflag) post_force_templated<1,0,1,1,0,1>();
else post_force_templated<1,0,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,0,1,0,1,1>();
else post_force_templated<1,0,1,0,1,0>();
else
if (zeroflag) post_force_templated<1,0,1,0,0,1>();
else post_force_templated<1,0,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,0,0,1,1,1>();
else post_force_templated<1,0,0,1,1,0>();
else
if (zeroflag) post_force_templated<1,0,0,1,0,1>();
else post_force_templated<1,0,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,0,0,0,1,1>();
else post_force_templated<1,0,0,0,1,0>();
else
if (zeroflag) post_force_templated<1,0,0,0,0,1>();
else post_force_templated<1,0,0,0,0,0>();
else
if (gjfflag)
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,1,1,1,1,1>();
else post_force_templated<0,1,1,1,1,0>();
else
if (zeroflag) post_force_templated<0,1,1,1,0,1>();
else post_force_templated<0,1,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,1,1,0,1,1>();
else post_force_templated<0,1,1,0,1,0>();
else
if (zeroflag) post_force_templated<0,1,1,0,0,1>();
else post_force_templated<0,1,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,1,0,1,1,1>();
else post_force_templated<0,1,0,1,1,0>();
else
if (zeroflag) post_force_templated<0,1,0,1,0,1>();
else post_force_templated<0,1,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,1,0,0,1,1>();
else post_force_templated<0,1,0,0,1,0>();
else
if (zeroflag) post_force_templated<0,1,0,0,0,1>();
else post_force_templated<0,1,0,0,0,0>();
else
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,0,1,1,1,1>();
else post_force_templated<0,0,1,1,1,0>();
else
if (zeroflag) post_force_templated<0,0,1,1,0,1>();
else post_force_templated<0,0,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,0,1,0,1,1>();
else post_force_templated<0,0,1,0,1,0>();
else
if (zeroflag) post_force_templated<0,0,1,0,0,1>();
else post_force_templated<0,0,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,0,0,1,1,1>();
else post_force_templated<0,0,0,1,1,0>();
else
if (zeroflag) post_force_templated<0,0,0,1,0,1>();
else post_force_templated<0,0,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,0,0,0,1,1>();
else post_force_templated<0,0,0,0,1,0>();
else
if (zeroflag) post_force_templated<0,0,0,0,0,1>();
else post_force_templated<0,0,0,0,0,0>();
#else
post_force_untemplated(int(tstyle==ATOM), gjfflag, tallyflag,
int(tbiasflag==BIAS), int(rmass!=NULL), zeroflag);
#endif
}
/* ---------------------------------------------------------------------- */
void FixLangevin::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ----------------------------------------------------------------------
modify forces using one of the many Langevin styles
------------------------------------------------------------------------- */
#ifdef TEMPLATED_FIX_LANGEVIN
template < int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO >
void FixLangevin::post_force_templated()
#else
void FixLangevin::post_force_untemplated
(int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO)
#endif
{
double gamma1,gamma2;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// apply damping and thermostat to atoms in group
// for Tp_TSTYLEATOM:
// use per-atom per-coord target temperature
// for Tp_GJF:
// use Gronbech-Jensen/Farago algorithm
// else use regular algorithm
// for Tp_TALLY:
// store drag plus random forces in flangevin[nlocal][3]
// for Tp_BIAS:
// calculate temperature since some computes require temp
// computed on current nlocal atoms to remove bias
// test v = 0 since some computes mask non-participating atoms via v = 0
// and added force has extra term not multiplied by v = 0
// for Tp_RMASS:
// use per-atom masses
// else use per-type masses
// for Tp_ZERO:
// sum random force over all atoms in group
// subtract sum/count from each atom in group
double fdrag[3],fran[3],fsum[3],fsumall[3];
bigint count;
double fswap;
double boltz = force->boltz;
double dt = update->dt;
double mvv2e = force->mvv2e;
double ftm2v = force->ftm2v;
compute_target();
if (Tp_ZERO) {
fsum[0] = fsum[1] = fsum[2] = 0.0;
count = group->count(igroup);
if (count == 0)
error->all(FLERR,"Cannot zero Langevin force of 0 atoms");
}
// reallocate flangevin if necessary
if (Tp_TALLY) {
if (atom->nlocal > maxatom1) {
memory->destroy(flangevin);
maxatom1 = atom->nmax;
memory->create(flangevin,maxatom1,3,"langevin:flangevin");
}
}
if (Tp_BIAS) temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (Tp_TSTYLEATOM) tsqrt = sqrt(tforce[i]);
if (Tp_RMASS) {
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
} else {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
}
fran[0] = gamma2*(random->uniform()-0.5);
fran[1] = gamma2*(random->uniform()-0.5);
fran[2] = gamma2*(random->uniform()-0.5);
if (Tp_BIAS) {
temperature->remove_bias(i,v[i]);
fdrag[0] = gamma1*v[i][0];
fdrag[1] = gamma1*v[i][1];
fdrag[2] = gamma1*v[i][2];
if (v[i][0] == 0.0) fran[0] = 0.0;
if (v[i][1] == 0.0) fran[1] = 0.0;
if (v[i][2] == 0.0) fran[2] = 0.0;
temperature->restore_bias(i,v[i]);
} else {
fdrag[0] = gamma1*v[i][0];
fdrag[1] = gamma1*v[i][1];
fdrag[2] = gamma1*v[i][2];
}
if (Tp_GJF) {
fswap = 0.5*(fran[0]+franprev[i][0]);
franprev[i][0] = fran[0];
fran[0] = fswap;
fswap = 0.5*(fran[1]+franprev[i][1]);
franprev[i][1] = fran[1];
fran[1] = fswap;
fswap = 0.5*(fran[2]+franprev[i][2]);
franprev[i][2] = fran[2];
fran[2] = fswap;
fdrag[0] *= gjffac;
fdrag[1] *= gjffac;
fdrag[2] *= gjffac;
fran[0] *= gjffac;
fran[1] *= gjffac;
fran[2] *= gjffac;
f[i][0] *= gjffac;
f[i][1] *= gjffac;
f[i][2] *= gjffac;
}
f[i][0] += fdrag[0] + fran[0];
f[i][1] += fdrag[1] + fran[1];
f[i][2] += fdrag[2] + fran[2];
if (Tp_TALLY) {
flangevin[i][0] = fdrag[0] + fran[0];
flangevin[i][1] = fdrag[1] + fran[1];
flangevin[i][2] = fdrag[2] + fran[2];
}
if (Tp_ZERO) {
fsum[0] += fran[0];
fsum[1] += fran[1];
fsum[2] += fran[2];
}
}
}
// set total force to zero
if (Tp_ZERO) {
MPI_Allreduce(fsum,fsumall,3,MPI_DOUBLE,MPI_SUM,world);
fsumall[0] /= count;
fsumall[1] /= count;
fsumall[2] /= count;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
f[i][0] -= fsumall[0];
f[i][1] -= fsumall[1];
f[i][2] -= fsumall[2];
}
}
}
// thermostat omega and angmom
if (oflag) omega_thermostat();
if (ascale) angmom_thermostat();
}
/* ----------------------------------------------------------------------
set current t_target and t_sqrt
------------------------------------------------------------------------- */
void FixLangevin::compute_target()
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
// if variable temp, evaluate variable, wrap with clear/add
// reallocate tforce array if necessary
if (tstyle == CONSTANT) {
t_target = t_start + delta * (t_stop-t_start);
tsqrt = sqrt(t_target);
} else {
modify->clearstep_compute();
if (tstyle == EQUAL) {
t_target = input->variable->compute_equal(tvar);
if (t_target < 0.0)
error->one(FLERR,"Fix langevin variable returned negative temperature");
tsqrt = sqrt(t_target);
} else {
if (nlocal > maxatom2) {
maxatom2 = atom->nmax;
memory->destroy(tforce);
memory->create(tforce,maxatom2,"langevin:tforce");
}
input->variable->compute_atom(tvar,igroup,tforce,1,0);
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (tforce[i] < 0.0)
error->one(FLERR,
"Fix langevin variable returned negative temperature");
}
modify->addstep_compute(update->ntimestep + 1);
}
}
/* ----------------------------------------------------------------------
thermostat rotational dof via omega
------------------------------------------------------------------------- */
void FixLangevin::omega_thermostat()
{
double gamma1,gamma2;
double boltz = force->boltz;
double dt = update->dt;
double mvv2e = force->mvv2e;
double ftm2v = force->ftm2v;
double **torque = atom->torque;
double **omega = atom->omega;
double *radius = atom->radius;
double *rmass = atom->rmass;
int *mask = atom->mask;
int *type = atom->type;
int nlocal = atom->nlocal;
// rescale gamma1/gamma2 by 10/3 & sqrt(10/3) for spherical particles
// does not affect rotational thermosatting
// gives correct rotational diffusivity behavior
double tendivthree = 10.0/3.0;
double tran[3];
double inertiaone;
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && (radius[i] > 0.0)) {
inertiaone = SINERTIA*radius[i]*radius[i]*rmass[i];
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
gamma1 = -tendivthree*inertiaone / t_period / ftm2v;
gamma2 = sqrt(inertiaone) * sqrt(80.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
tran[0] = gamma2*(random->uniform()-0.5);
tran[1] = gamma2*(random->uniform()-0.5);
tran[2] = gamma2*(random->uniform()-0.5);
torque[i][0] += gamma1*omega[i][0] + tran[0];
torque[i][1] += gamma1*omega[i][1] + tran[1];
torque[i][2] += gamma1*omega[i][2] + tran[2];
}
}
}
/* ----------------------------------------------------------------------
thermostat rotational dof via angmom
------------------------------------------------------------------------- */
void FixLangevin::angmom_thermostat()
{
double gamma1,gamma2;
double boltz = force->boltz;
double dt = update->dt;
double mvv2e = force->mvv2e;
double ftm2v = force->ftm2v;
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
double **torque = atom->torque;
double **angmom = atom->angmom;
double *rmass = atom->rmass;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int *type = atom->type;
int nlocal = atom->nlocal;
// rescale gamma1/gamma2 by ascale for aspherical particles
// does not affect rotational thermosatting
// gives correct rotational diffusivity behavior if (nearly) spherical
// any value will be incorrect for rotational diffusivity if aspherical
double inertia[3],omega[3],tran[3];
double *shape,*quat;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
shape = bonus[ellipsoid[i]].shape;
inertia[0] = EINERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
inertia[1] = EINERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
inertia[2] = EINERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
quat = bonus[ellipsoid[i]].quat;
MathExtra::mq_to_omega(angmom[i],quat,inertia,omega);
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
gamma1 = -ascale / t_period / ftm2v;
gamma2 = sqrt(ascale*24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
tran[0] = sqrt(inertia[0])*gamma2*(random->uniform()-0.5);
tran[1] = sqrt(inertia[1])*gamma2*(random->uniform()-0.5);
tran[2] = sqrt(inertia[2])*gamma2*(random->uniform()-0.5);
torque[i][0] += inertia[0]*gamma1*omega[0] + tran[0];
torque[i][1] += inertia[1]*gamma1*omega[1] + tran[1];
torque[i][2] += inertia[2]*gamma1*omega[2] + tran[2];
}
}
}
/* ----------------------------------------------------------------------
tally energy transfer to thermal reservoir
------------------------------------------------------------------------- */
void FixLangevin::end_of_step()
{
if (!tallyflag) return;
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
energy_onestep = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2];
energy += energy_onestep*update->dt;
}
/* ---------------------------------------------------------------------- */
void FixLangevin::reset_target(double t_new)
{
t_target = t_start = t_stop = t_new;
}
/* ---------------------------------------------------------------------- */
void FixLangevin::reset_dt()
{
if (atom->mass) {
for (int i = 1; i <= atom->ntypes; i++) {
gfactor2[i] = sqrt(atom->mass[i]) *
sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v;
gfactor2[i] *= 1.0/sqrt(ratio[i]);
}
}
}
/* ---------------------------------------------------------------------- */
int FixLangevin::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Could not find fix_modify temperature ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all(FLERR,
"Fix_modify temperature ID does not compute temperature");
if (temperature->igroup != igroup && comm->me == 0)
error->warning(FLERR,"Group for fix_modify temp != fix group");
return 2;
}
return 0;
}
/* ---------------------------------------------------------------------- */
double FixLangevin::compute_scalar()
{
- if (!tallyflag || flangevin == NULL) return 0.0;
+ if (!tallyflag) return 0.0;
// capture the very first energy transfer to thermal reservoir
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (update->ntimestep == update->beginstep) {
energy_onestep = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2];
energy = 0.5*energy_onestep*update->dt;
}
// convert midstep energy back to previous fullstep energy
double energy_me = energy - 0.5*energy_onestep*update->dt;
double energy_all;
MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
return -energy_all;
}
/* ----------------------------------------------------------------------
extract thermostat properties
------------------------------------------------------------------------- */
void *FixLangevin::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"t_target") == 0) {
return &t_target;
}
return NULL;
}
/* ----------------------------------------------------------------------
memory usage of tally array
------------------------------------------------------------------------- */
double FixLangevin::memory_usage()
{
double bytes = 0.0;
if (gjfflag) bytes += atom->nmax*3 * sizeof(double);
if (tallyflag) bytes += atom->nmax*3 * sizeof(double);
if (tforce) bytes += atom->nmax * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
allocate atom-based array for franprev
------------------------------------------------------------------------- */
void FixLangevin::grow_arrays(int nmax)
{
memory->grow(franprev,nmax,3,"fix_langevin:franprev");
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
------------------------------------------------------------------------- */
void FixLangevin::copy_arrays(int i, int j, int delflag)
{
for (int m = 0; m < nvalues; m++)
franprev[j][m] = franprev[i][m];
}
/* ----------------------------------------------------------------------
pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */
int FixLangevin::pack_exchange(int i, double *buf)
{
for (int m = 0; m < nvalues; m++) buf[m] = franprev[i][m];
return nvalues;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */
int FixLangevin::unpack_exchange(int nlocal, double *buf)
{
for (int m = 0; m < nvalues; m++) franprev[nlocal][m] = buf[m];
return nvalues;
}