Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88063545
fix_box_relax.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
Wed, Oct 16, 14:23
Size
16 KB
Mime Type
text/x-c
Expires
Fri, Oct 18, 14:23 (2 d)
Engine
blob
Format
Raw Data
Handle
21707473
Attached To
rLAMMPS lammps
fix_box_relax.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: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "fix_box_relax.h"
#include "atom.h"
#include "domain.h"
#include "update.h"
#include "comm.h"
#include "force.h"
#include "modify.h"
#include "compute.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
enum{XYZ,XY,YZ,XZ,ANISO};
/* ---------------------------------------------------------------------- */
FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 4) error->all("Illegal fix box/relax command");
box_change = 1;
if (strcmp(arg[3],"xyz") == 0) {
if (narg < 5) error->all("Illegal fix box/relax command");
press_couple = XYZ;
p_target[0] = p_target[1] = p_target[2] = atof(arg[4]);
p_flag[0] = p_flag[1] = p_flag[2] = 1;
if (domain->dimension == 2) p_flag[2] = 0;
} else {
if (strcmp(arg[3],"xy") == 0) press_couple = XY;
else if (strcmp(arg[3],"yz") == 0) press_couple = YZ;
else if (strcmp(arg[3],"xz") == 0) press_couple = XZ;
else if (strcmp(arg[3],"aniso") == 0) press_couple = ANISO;
else error->all("Illegal fix box/relax command");
if (narg < 7) error->all("Illegal fix box/relax command");
if (domain->dimension == 2 &&
(press_couple == XY || press_couple == YZ || press_couple == XZ))
error->all("Invalid fix box/relax command for a 2d simulation");
if (strcmp(arg[4],"NULL") == 0) {
p_target[0] = 0.0;
p_flag[0] = 0;
} else {
p_target[0] = atof(arg[4]);
p_flag[0] = 1;
}
if (strcmp(arg[5],"NULL") == 0) {
p_target[1] = 0.0;
p_flag[1] = 0;
} else {
p_target[1] = atof(arg[5]);
p_flag[1] = 1;
}
if (strcmp(arg[6],"NULL") == 0) {
p_target[2] = 0.0;
p_flag[2] = 0;
} else {
if (domain->dimension == 2)
error->all("Invalid fix box/relax command for a 2d simulation");
p_target[2] = atof(arg[6]);
p_flag[2] = 1;
}
}
pflagsum = p_flag[0] + p_flag[1] + p_flag[2];
// process extra keywords
allremap = 1;
vmax = 0.0001;
int iarg;
if (press_couple == XYZ) iarg = 5;
else iarg = 7;
while (iarg < narg) {
if (strcmp(arg[iarg],"dilate") == 0) {
if (iarg+2 > narg) error->all("Illegal fix box/relax command");
if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0;
else error->all("Illegal fix box/relax command");
iarg += 2;
} else if (strcmp(arg[iarg],"vmax") == 0) {
if (iarg+2 > narg) error->all("Illegal fix box/relax command");
vmax = atof(arg[iarg+1]);
iarg += 2;
} else error->all("Illegal fix box/relax command");
}
// error checks
if (press_couple == XY && (p_flag[0] == 0 || p_flag[1] == 0))
error->all("Illegal fix box/relax command");
if (press_couple == YZ && (p_flag[1] == 0 || p_flag[2] == 0))
error->all("Illegal fix box/relax command");
if (press_couple == XZ && (p_flag[0] == 0 || p_flag[2] == 0))
error->all("Illegal fix box/relax command");
if (press_couple == XY && p_target[0] != p_target[1])
error->all("Illegal fix box/relax command");
if (press_couple == YZ && p_target[1] != p_target[2])
error->all("Illegal fix box/relax command");
if (press_couple == XZ && p_target[0] != p_target[2])
error->all("Illegal fix box/relax command");
if (p_flag[0] && domain->xperiodic == 0)
error->all("Cannot use fix box/relax on a non-periodic dimension");
if (p_flag[1] && domain->yperiodic == 0)
error->all("Cannot use fix box/relax on a non-periodic dimension");
if (p_flag[2] && domain->zperiodic == 0)
error->all("Cannot use fix box/relax on a non-periodic dimension");
if (vmax <= 0.0) error->all("Illegal fix box/relax command");
// create a new compute temp style
// id = fix-ID + temp
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
tflag = 1;
// create a new compute pressure style (virial only)
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[5];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
newarg[4] = (char *) "virial";
modify->add_compute(5,newarg);
delete [] newarg;
pflag = 1;
dimension = domain->dimension;
nrigid = 0;
rfix = 0;
current_lifo = 0;
}
/* ---------------------------------------------------------------------- */
FixBoxRelax::~FixBoxRelax()
{
delete [] rfix;
// delete temperature and pressure if fix created them
if (tflag) modify->delete_compute(id_temp);
if (pflag) modify->delete_compute(id_press);
delete [] id_temp;
delete [] id_press;
}
/* ---------------------------------------------------------------------- */
int FixBoxRelax::setmask()
{
int mask = 0;
mask |= MIN_ENERGY;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixBoxRelax::init()
{
// set temperature and pressure ptrs
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all("Temperature ID for fix box/relax does not exist");
temperature = modify->compute[icompute];
icompute = modify->find_compute(id_press);
if (icompute < 0) error->all("Pressure ID for fix box/relax does not exist");
pressure = modify->compute[icompute];
// initial box dimensions
xprdinit = domain->xprd;
yprdinit = domain->yprd;
zprdinit = domain->zprd;
if (dimension == 3) volinit = domain->xprd*domain->yprd*domain->zprd;
else volinit = domain->xprd*domain->yprd;
pv2e = 1.0 / force->nktv2p;
// detect if any rigid fixes exist so rigid bodies move when box is remapped
// rfix[] = indices to each fix rigid
delete [] rfix;
nrigid = 0;
rfix = NULL;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->rigid_flag) nrigid++;
if (nrigid) {
rfix = new int[nrigid];
nrigid = 0;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i;
}
}
/* ----------------------------------------------------------------------
compute energy and force due to extra degrees of freedom
------------------------------------------------------------------------- */
double FixBoxRelax::min_energy(double *fextra)
{
double eng,scale,scalex,scaley,scalez;
double t_current = temperature->compute_scalar();
if (press_couple == XYZ) {
double tmp = pressure->compute_scalar();
} else {
temperature->compute_vector();
pressure->compute_vector();
}
couple();
// trigger virial computation on every iteration of minimizer
pressure->addstep(update->ntimestep+1);
// compute energy, forces for each extra degree of freedom
// returned eng = PV must be in units of energy
// returned fextra must likewise be in units of energy
if (press_couple == XYZ) {
scale = domain->xprd/xprdinit;
if (dimension == 3) {
eng = pv2e * p_target[0] * (scale*scale*scale-1.0)*volinit;
fextra[0] = pv2e * (p_current[0] - p_target[0])*3.0*scale*scale*volinit;
} else {
eng = pv2e * p_target[0] * (scale*scale-1.0)*volinit;
fextra[0] = pv2e * (p_current[0] - p_target[0])*2.0*scale*volinit;
}
} else {
fextra[0] = fextra[1] = fextra[2] = 0.0;
scalex = scaley = scalez = 1.0;
if (p_flag[0]) scalex = domain->xprd/xprdinit;
if (p_flag[1]) scaley = domain->yprd/yprdinit;
if (p_flag[2]) scalez = domain->zprd/zprdinit;
eng = pv2e * (p_flag[0]*p_target[0] + p_flag[1]*p_target[1] +
p_flag[2]*p_target[2])/pflagsum *
(scalex*scaley*scalez-1.0)*volinit;
if (p_flag[0])
fextra[0] = pv2e * (p_current[0] - p_target[0])*scaley*scalez*volinit;
if (p_flag[1])
fextra[1] = pv2e * (p_current[1] - p_target[1])*scalex*scalez*volinit;
if (p_flag[2])
fextra[2] = pv2e * (p_current[2] - p_target[2])*scalex*scaley*volinit;
}
return eng;
}
/* ----------------------------------------------------------------------
store extra dof values for minimization linesearch starting point
boxlo0,boxhi0 = box dimensions
s0 = ratio of current boxsize to initial boxsize
box values are pushed onto a LIFO stack so nested calls can be made
values are popped by calling min_step(0.0)
------------------------------------------------------------------------- */
void FixBoxRelax::min_store()
{
for (int i = 0; i < 3; i++) {
boxlo0[current_lifo][i] = domain->boxlo[i];
boxhi0[current_lifo][i] = domain->boxhi[i];
}
s0[0] = (boxhi0[current_lifo][0]-boxlo0[current_lifo][0])/xprdinit;
s0[1] = (boxhi0[current_lifo][1]-boxlo0[current_lifo][1])/yprdinit;
s0[2] = (boxhi0[current_lifo][2]-boxlo0[current_lifo][2])/zprdinit;
}
/* ----------------------------------------------------------------------
clear the LIFO stack for min_store
------------------------------------------------------------------------- */
void FixBoxRelax::min_clearstore()
{
current_lifo = 0;
}
/* ----------------------------------------------------------------------
push the LIFO stack for min_store
------------------------------------------------------------------------- */
void FixBoxRelax::min_pushstore()
{
if (current_lifo >= MAX_LIFO_DEPTH) {
error->all("Attempt to push beyond stack limit in fix box/relax");
return;
}
current_lifo++;
}
/* ----------------------------------------------------------------------
pop the LIFO stack for min_store
------------------------------------------------------------------------- */
void FixBoxRelax::min_popstore()
{
if (current_lifo <= 0) {
error->all("Attempt to pop empty stack in fix box/relax");
return;
}
current_lifo--;
}
/* ----------------------------------------------------------------------
change the box dimensions by fraction ds = alpha*hextra
------------------------------------------------------------------------- */
void FixBoxRelax::min_step(double alpha, double *hextra)
{
if (press_couple == XYZ) {
ds[0] = ds[1] = ds[2] = alpha*hextra[0];
} else {
if (p_flag[0]) ds[0] = alpha*hextra[0];
if (p_flag[1]) ds[1] = alpha*hextra[1];
if (p_flag[2]) ds[2] = alpha*hextra[2];
}
remap();
}
/* ----------------------------------------------------------------------
max allowed step size along hextra
------------------------------------------------------------------------- */
double FixBoxRelax::max_alpha(double *hextra)
{
double alpha = 0.0;
if (press_couple == XYZ) alpha = vmax/fabs(hextra[0]);
else {
alpha = vmax/fabs(hextra[0]);
alpha = MIN(alpha,vmax/fabs(hextra[1]));
alpha = MIN(alpha,vmax/fabs(hextra[2]));
}
return alpha;
}
/* ----------------------------------------------------------------------
return number of degrees of freedom added by this fix
------------------------------------------------------------------------- */
int FixBoxRelax::min_dof()
{
if (press_couple == XYZ) return 1;
return 3;
}
/* ----------------------------------------------------------------------
dilate the box and owned/ghost atoms around center of box
------------------------------------------------------------------------- */
void FixBoxRelax::remap()
{
int i,n;
double ctr;
// ctr = geometric center of box in a dimension
// rescale simulation box from linesearch starting point
// scale atom coords for all atoms or only for fix group atoms
double **x = atom->x;
int *mask = atom->mask;
n = atom->nlocal + atom->nghost;
// convert pertinent atoms and rigid bodies to lamda coords
if (allremap) domain->x2lamda(n);
else {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
domain->x2lamda(x[i],x[i]);
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(0);
// reset global and local box to new size/shape
for (i = 0; i < 3; i++)
if (p_flag[i]) {
double currentBoxLo0 = boxlo0[current_lifo][i];
double currentBoxHi0 = boxhi0[current_lifo][i];
ctr = 0.5 * (currentBoxLo0 + currentBoxHi0);
domain->boxlo[i] = currentBoxLo0 + (currentBoxLo0-ctr)*ds[i]/s0[i];
domain->boxhi[i] = currentBoxHi0 + (currentBoxHi0-ctr)*ds[i]/s0[i];
}
domain->set_global_box();
domain->set_local_box();
// convert pertinent atoms and rigid bodies back to box coords
if (allremap) domain->lamda2x(n);
else {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
domain->lamda2x(x[i],x[i]);
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(1);
}
/* ---------------------------------------------------------------------- */
void FixBoxRelax::couple()
{
double *tensor = pressure->vector;
if (press_couple == XYZ)
p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
else if (press_couple == XY) {
double ave = 0.5 * (tensor[0] + tensor[1]);
p_current[0] = p_current[1] = ave;
p_current[2] = tensor[2];
} else if (press_couple == YZ) {
double ave = 0.5 * (tensor[1] + tensor[2]);
p_current[1] = p_current[2] = ave;
p_current[0] = tensor[0];
} else if (press_couple == XZ) {
double ave = 0.5 * (tensor[0] + tensor[2]);
p_current[0] = p_current[2] = ave;
p_current[1] = tensor[1];
} else if (press_couple == ANISO) {
p_current[0] = tensor[0];
p_current[1] = tensor[1];
p_current[2] = tensor[2];
}
}
/* ---------------------------------------------------------------------- */
int FixBoxRelax::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all("Illegal fix_modify command");
if (tflag) {
modify->delete_compute(id_temp);
tflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all("Could not find fix_modify temperature ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all("Fix_modify temperature ID does not compute temperature");
if (temperature->igroup != 0 && comm->me == 0)
error->warning("Temperature for fix modify is not for group all");
// reset id_temp of pressure to new temperature ID
icompute = modify->find_compute(id_press);
if (icompute < 0) error->all("Pressure ID for fix modify does not exist");
modify->compute[icompute]->reset_extra_compute_fix(id_temp);
return 2;
} else if (strcmp(arg[0],"press") == 0) {
if (narg < 2) error->all("Illegal fix_modify command");
if (pflag) {
modify->delete_compute(id_press);
pflag = 0;
}
delete [] id_press;
int n = strlen(arg[1]) + 1;
id_press = new char[n];
strcpy(id_press,arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all("Could not find fix_modify pressure ID");
pressure = modify->compute[icompute];
if (pressure->pressflag == 0)
error->all("Fix_modify pressure ID does not compute pressure");
return 2;
}
return 0;
}
Event Timeline
Log In to Comment