Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85687643
fix_ehex.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
Tue, Oct 1, 05:33
Size
15 KB
Mime Type
text/x-c
Expires
Thu, Oct 3, 05:33 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21242743
Attached To
rLAMMPS lammps
fix_ehex.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 author: Peter Wirnsberger (University of Cambridge)
This source file implements the asymmetric version of the enhanced heat
exchange (eHEX/a) algorithm. The paper is available for download on
arXiv: http://arxiv.org/pdf/1507.07081.pdf.
This file is based on fix_heat.cpp written by Paul Crozier (SNL)
which implements the heat exchange (HEX) algorithm.
------------------------------------------------------------------------- */
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "fix_ehex.h"
#include "atom.h"
#include "domain.h"
#include "region.h"
#include "group.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
#include "fix_shake.h"
#include "neighbor.h"
#include "comm.h"
#include "timer.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{CONSTANT,EQUAL,ATOM};
/* ---------------------------------------------------------------------- */
FixEHEX::FixEHEX(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
MPI_Comm_rank(world, &me);
// check
if (narg < 4) error->all(FLERR,"Illegal fix ehex command: wrong number of parameters ");
scalar_flag = 1;
global_freq = 1;
extscalar = 0;
// apply fix every nevery timesteps
nevery = force->inumeric(FLERR,arg[3]);
if (nevery <= 0) error->all(FLERR,"Illegal fix ehex command");
// heat flux into the reservoir
heat_input = force->numeric(FLERR,arg[4]);
// optional args
iregion = -1;
idregion = NULL;
// NOTE: constraints are deactivated by default
constraints = 0;
// NOTE: cluster rescaling is deactivated by default
cluster = 0;
// NOTE: hex = 1 means that no coordinate correction is applied in which case eHEX reduces to HEX
hex = 0;
int iarg = 5;
while (iarg < narg) {
if (strcmp(arg[iarg],"region") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ehex command: wrong number of parameters ");
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1)
error->all(FLERR,"Region ID for fix ehex does not exist");
int n = strlen(arg[iarg+1]) + 1;
idregion = new char[n];
strcpy(idregion,arg[iarg+1]);
iarg += 2;
}
// apply constraints (shake/rattle) at the end of the timestep
else if (strcmp(arg[iarg], "constrain") == 0) {
constraints = 1;
iarg += 1;
}
// rescale only if the entire molecule is contained within the region
else if (strcmp(arg[iarg], "com") == 0) {
cluster = 1;
iarg += 1;
}
// don't apply a coordinate correction if this keyword is specified
else if (strcmp(arg[iarg], "hex") == 0) {
hex = 1;
iarg+= 1;
}
else
error->all(FLERR, "Illegal fix ehex keyword ");
}
// check options
if (cluster && !constraints)
error->all(FLERR, "You can only use the keyword 'com' together with the keyword 'constrain' ");
scale = 1.0;
scalingmask = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
}
/* ---------------------------------------------------------------------- */
void FixEHEX::grow_arrays(int nmax) {
memory->grow(scalingmask, nmax,"ehex:scalingmask");
}
/* ---------------------------------------------------------------------- */
FixEHEX::~FixEHEX()
{
atom->delete_callback(id,0);
delete [] idregion;
memory->destroy(scalingmask);
}
/* ---------------------------------------------------------------------- */
int FixEHEX::setmask()
{
int mask = 0;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixEHEX::init()
{
// set index and check validity of region
if (iregion >= 0) {
iregion = domain->find_region(idregion);
if (iregion == -1)
error->all(FLERR,"Region ID for fix ehex does not exist");
}
// cannot have 0 atoms in group
if (group->count(igroup) == 0)
error->all(FLERR,"Fix ehex group has no atoms");
fshake = NULL;
if (constraints) {
// check if constraining algorithm is used (FixRattle inherits from FixShake)
int cnt_shake = 0;
int id_shake;
for (int i = 0; i < modify->nfix; i++) {
if (strcmp("rattle", modify->fix[i]->style) == 0 ||
strcmp("shake", modify->fix[i]->style) == 0) {
cnt_shake++;
id_shake = i;
}
}
if (cnt_shake > 1)
error->all(FLERR,"Multiple instances of fix shake/rattle detected (not supported yet)");
else if (cnt_shake == 1) {
fshake = ((FixShake*) modify->fix[id_shake]);
}
else if (cnt_shake == 0)
error->all(FLERR, "Fix ehex was configured with keyword constrain, but shake/rattle was not defined");
}
}
/* ---------------------------------------------------------------------- */
void FixEHEX::end_of_step() {
// store local pointers
x = atom->x;
f = atom->f;
v = atom->v;
mass = atom->mass;
rmass = atom->rmass;
type = atom->type;
nlocal = atom->nlocal;
// determine which sites are to be rescaled
update_scalingmask();
// rescale velocities
rescale();
// if required use shake/rattle to correct coordinates and velocities
if (constraints && fshake)
fshake->shake_end_of_step(0);
}
/* ----------------------------------------------------------------------
Iterate over all atoms, rescale the velocities and apply coordinate
corrections.
------------------------------------------------------------------------- */
void FixEHEX::rescale() {
double heat,ke,massone;
double Kr, Ke, escale;
double vsub[3],vcm[3], sfr[3];
double dvcmdt[3];
double mi;
double dt;
double F, mr, Fo2Kr, epsr_ik, sfvr, eta_ik;
dt = update->dt;
// calculate centre of mass properties
com_properties(vcm, sfr, &sfvr, &Ke, &Kr, &masstotal);
// heat flux into the reservoir
F = heat_input * force->ftm2v * nevery;
// total mass
mr = masstotal;
Fo2Kr = F / (2.*Kr);
// energy scaling factor
escale = 1. + (F*dt)/Kr;
// safety check for kinetic energy
if (escale < 0.0) error->all(FLERR,"Fix ehex kinetic energy went negative");
scale = sqrt(escale);
vsub[0] = (scale-1.0) * vcm[0];
vsub[1] = (scale-1.0) * vcm[1];
vsub[2] = (scale-1.0) * vcm[2];
for (int i = 0; i < nlocal; i++){
if (scalingmask[i]) {
mi = (rmass) ? rmass[i] : mass[type[i]];
for (int k=0; k<3; k++) {
// apply coordinate correction unless running in hex mode
if (!hex) {
// epsr_ik implements Eq. (20) in the paper
eta_ik = mi * F/(2.*Kr) * (v[i][k] - vcm[k]);
epsr_ik = eta_ik / (mi*Kr) * (F/48. + sfvr/6.*force->ftm2v) - F/(12.*Kr) * (f[i][k]/mi - sfr[k]/mr)*force->ftm2v;
x[i][k] -= dt*dt*dt * epsr_ik;
}
// rescale the velocity
v[i][k] = scale*v[i][k] - vsub[k];
}
}
}
}
/* ---------------------------------------------------------------------- */
double FixEHEX::compute_scalar()
{
return scale;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixEHEX::memory_usage()
{
double bytes = 0.0;
bytes += atom->nmax * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
Update the array scalingmask depending on which individual atoms
will be rescaled or not.
------------------------------------------------------------------------- */
void FixEHEX::update_scalingmask() {
int m;
int lid;
bool stat;
int nsites;
// prematch region
Region *region = NULL;
if (iregion >= 0) {
region = domain->regions[iregion];
region->prematch();
}
// only rescale molecules whose center of mass if fully contained in the region
if (cluster) {
// loop over all clusters
for (int i=0; i < fshake->nlist; i++) {
// cluster id
m = fshake->list[i];
// check if the centre of mass of the cluster is inside the region
// if region == NULL, just check the group information of all sites
if (fshake->shake_flag[m] == 1) nsites = 3;
else if (fshake->shake_flag[m] == 2) nsites = 2;
else if (fshake->shake_flag[m] == 3) nsites = 3;
else if (fshake->shake_flag[m] == 4) nsites = 4;
else nsites = 0;
if (nsites == 0) {
error->all(FLERR,"Internal error: shake_flag[m] has to be between 1 and 4 for m in nlist");
}
stat = check_cluster(fshake->shake_atom[m], nsites, region);
for (int l=0; l < nsites; l++) {
lid = atom->map(fshake->shake_atom[m][l]);
scalingmask[lid] = stat;
}
}
// check atoms that do not belong to any cluster
for (int i=0; i<atom->nlocal; i++) {
if (fshake->shake_flag[i] == 0)
scalingmask[i] = rescale_atom(i,region);
}
}
// no clusters, just individual sites (e.g. monatomic system or flexible molecules)
else {
for (int i=0; i<atom->nlocal; i++)
scalingmask[i] = rescale_atom(i,region);
}
}
/* ----------------------------------------------------------------------
Check if the centre of mass of the cluster to be constrained is
inside the region.
------------------------------------------------------------------------- */
bool FixEHEX::check_cluster(int *shake_atom, int n, Region * region) {
// IMPORTANT NOTE: If any site of the cluster belongs to a group
// which should not be rescaled than all of the sites
// will be ignored!
double **x = atom->x;
double * rmass = atom->rmass;
double * mass = atom->mass;
int * type = atom->type;
int * mask = atom->mask;
double xcom[3], xtemp[3];
double mcluster, mi;
bool stat;
int lid[4];
// accumulate mass and centre of mass position
stat = true;
xcom[0] = 0.;
xcom[1] = 0.;
xcom[2] = 0.;
mcluster = 0;
for (int i = 0; i < n; i++) {
// get local id
lid[i] = atom->map(shake_atom[i]);
// check if all sites of the cluster belong to the correct group
stat = stat && (mask[lid[i]] & groupbit);
if (region && stat) {
// check if reduced mass is used
mi = (rmass) ? rmass[lid[i]] : mass[type[lid[i]]];
mcluster += mi;
// accumulate centre of mass
// NOTE: you can either use unwrapped coordinates or take site x[lid[0]] as reference,
// i.e. reconstruct the molecule around this site and calculate the com.
for (int k=0; k<3; k++)
xtemp[k] = x[lid[i]][k] - x[lid[0]][k];
// take into account pbc
domain->minimum_image(xtemp);
for (int k=0; k<3; k++)
xcom[k] += mi * (x[lid[0]][k] + xtemp[k]) ;
}
}
// check if centre of mass is inside the region (if specified)
if (region && stat) {
// check mass
if (mcluster < 1.e-14) {
error->all(FLERR, "Fix ehex shake cluster has almost zero mass.");
}
// divide by total mass
for (int k=0; k<3; k++)
xcom[k] = xcom[k]/mcluster;
// apply periodic boundary conditions (centre of mass could be outside the box)
// and check if molecule is inside the region
domain->remap(xcom);
stat = stat && region->match(xcom[0], xcom[1], xcom[2]);
}
return stat;
}
/* ----------------------------------------------------------------------
Check if atom i has the correct group and is inside the region.
------------------------------------------------------------------------- */
bool FixEHEX::rescale_atom(int i, Region*region) {
bool stat;
double x_r[3];
// check mask and group
stat = (atom->mask[i] & groupbit);
if (region) {
x_r[0] = atom->x[i][0];
x_r[1] = atom->x[i][1];
x_r[2] = atom->x[i][2];
// apply periodic boundary conditions
domain->remap(x_r);
// check if the atom is in the group/region
stat = stat && region->match(x_r[0],x_r[1],x_r[2]);
}
return stat;
}
/* ----------------------------------------------------------------------
Calculate global properties of the atoms inside the reservoir.
(e.g. com velocity, kinetic energy, total mass,...)
------------------------------------------------------------------------- */
void FixEHEX::com_properties(double * vr, double * sfr, double *sfvr, double *K, double *Kr, double *mr) {
double ** f = atom->f;
double ** v = atom->v;
int nlocal = atom->nlocal;
double *rmass= atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
double l_vr[3];
double l_mr;
double l_sfr[3];
double l_sfvr;
double l_K;
double mi;
double l_buf[9];
double buf[9];
// calculate partial sums
l_vr[0] = l_vr[1] = l_vr[2] = 0;
l_sfr[0] = l_sfr[1] = l_sfr[2] = 0;
l_sfvr = 0;
l_mr = 0;
l_K = 0;
for (int i = 0; i < nlocal; i++) {
if (scalingmask[i]) {
// check if reduced mass is used
mi = (rmass) ? rmass[i] : mass[type[i]];
// accumulate total mass
l_mr += mi;
// accumulate kinetic energy
l_K += mi/2. * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
// sum_j f_j * v_j
l_sfvr += f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2];
// accumulate com velocity and sum of forces
for (int k=0; k<3; k++) {
l_vr[k] += mi * v[i][k];
l_sfr[k]+= f[i][k];
}
}
}
// reduce sums
l_buf[0] = l_vr[0];
l_buf[1] = l_vr[1];
l_buf[2] = l_vr[2];
l_buf[3] = l_K;
l_buf[4] = l_mr;
l_buf[5] = l_sfr[0];
l_buf[6] = l_sfr[1];
l_buf[7] = l_sfr[2];
l_buf[8] = l_sfvr;
MPI_Allreduce(l_buf, buf, 9, MPI_DOUBLE, MPI_SUM, world);
// total mass of region
*mr = buf[4];
if (*mr < 1.e-14) {
error->all(FLERR, "Fix ehex error mass of region is close to zero");
}
// total kinetic energy of region
*K = buf[3];
// centre of mass velocity of region
vr[0] = buf[0]/(*mr);
vr[1] = buf[1]/(*mr);
vr[2] = buf[2]/(*mr);
// sum of forces
sfr[0] = buf[5];
sfr[1] = buf[6];
sfr[2] = buf[7];
// calculate non-translational kinetic energy
*Kr = *K - 0.5* (*mr) * (vr[0]*vr[0]+vr[1]*vr[1]+vr[2]*vr[2]);
// calculate sum_j f_j * (v_j - v_r) = sum_j f_j * v_j - v_r * sum_j f_j
*sfvr = buf[8] - (vr[0]*sfr[0] + vr[1]*sfr[1] + vr[2]*sfr[2]);
}
Event Timeline
Log In to Comment