tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
27e8d0f19cfd60ff513828af74d07d2c8f3c4451
lammps-gran-kokkos/src/BOCS/fix_bocs.cpp
Stan Gerald Moore fc7119982b whitespace
2023-10-19 12:53:13 -04:00

2524 lines
72 KiB
C++
Raw Blame History

// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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.
-------------------------------------------------------------------------
BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
from The Pennsylvania State University
------------------------------------------------------------------------- */
#include "fix_bocs.h"
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "compute.h"
#include "compute_pressure_bocs.h"
#include "domain.h"
#include "error.h"
#include "fix_deform.h"
#include "force.h"
#include "group.h"
#include "irregular.h"
#include "kspace.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "respa.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
static const char cite_user_bocs_package[] =
"BOCS package: doi:10.1021/acs.jpcb.7b09993\n\n"
"@Article{Dunn2018,\n"
" author = {N. J. H. Dunn and K. M. Lebold and M. R. {DeLyser} and\n"
" J. F. Rudzinski and W. G. Noid},\n"
" title = {{BOCS}: Bottom-Up Open-Source Coarse-Graining Software},\n"
" journal = {J.~Phys.\\ Chem.~B},\n"
" year = 2018,\n"
" volume = 122,\n"
" number = 13,\n"
" pages = {3363--3377}\n"
"}\n\n";
#define DELTAFLIP 0.1
#define TILTMAX 1.5
enum{NOBIAS,BIAS};
enum{NONE,XYZ,XY,YZ,XZ};
enum{ISO,ANISO,TRICLINIC};
const int NUM_INPUT_DATA_COLUMNS = 2; // columns in the pressure correction file
/* ----------------------------------------------------------------------
NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion
---------------------------------------------------------------------- */
FixBocs::FixBocs(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
rfix(nullptr), id_dilate(nullptr), irregular(nullptr),
id_temp(nullptr), id_press(nullptr),
eta(nullptr), eta_dot(nullptr), eta_dotdot(nullptr),
eta_mass(nullptr), etap(nullptr), etap_dot(nullptr), etap_dotdot(nullptr),
etap_mass(nullptr)
{
if (lmp->citeme) lmp->citeme->add(cite_user_bocs_package);
if (narg < 4) error->all(FLERR,"Illegal fix bocs command");
restart_global = 1;
dynamic_group_allow = 1;
time_integrate = 1;
scalar_flag = 1;
vector_flag = 1;
global_freq = 1;
extscalar = 1;
extvector = 0;
ecouple_flag = 1;
// default values
pcouple = NONE;
drag = 0.0;
allremap = 1;
id_dilate = nullptr;
mtchain = mpchain = 3;
nc_tchain = nc_pchain = 1;
mtk_flag = 1;
deviatoric_flag = 0;
nreset_h0 = 0;
eta_mass_flag = 1;
omega_mass_flag = 0;
etap_mass_flag = 0;
flipflag = 1;
dipole_flag = 0;
dlm_flag = 0;
tcomputeflag = 0;
pcomputeflag = 0;
id_temp = nullptr;
id_press = nullptr;
p_match_coeffs = nullptr;
splines = nullptr;
spline_length = 0;
// turn on tilt factor scaling, whenever applicable
dimension = domain->dimension;
scaleyz = scalexz = scalexy = 0;
if (domain->yperiodic && domain->xy != 0.0) scalexy = 1;
if (domain->zperiodic && dimension == 3) {
if (domain->yz != 0.0) scaleyz = 1;
if (domain->xz != 0.0) scalexz = 1;
}
// set fixed-point to default = center of cell
fixedpoint[0] = 0.5*(domain->boxlo[0]+domain->boxhi[0]);
fixedpoint[1] = 0.5*(domain->boxlo[1]+domain->boxhi[1]);
fixedpoint[2] = 0.5*(domain->boxlo[2]+domain->boxhi[2]);
// used by FixNVTSllod to preserve non-default value
mtchain_default_flag = 1;
tstat_flag = 0;
double t_period = 0.0;
double p_period[6];
for (int i = 0; i < 6; i++) {
p_start[i] = p_stop[i] = p_period[i] = p_target[i] = 0.0;
p_flag[i] = 0;
}
// process keywords
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg],"temp") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix bocs command");
tstat_flag = 1;
t_start = utils::numeric(FLERR,arg[iarg+1],false,lmp);
t_target = t_start;
t_stop = utils::numeric(FLERR,arg[iarg+2],false,lmp);
t_period = utils::numeric(FLERR,arg[iarg+3],false,lmp);
if (t_start <= 0.0 || t_stop <= 0.0)
error->all(FLERR,
"Target temperature for fix bocs cannot be 0.0");
iarg += 4;
} else if (strcmp(arg[iarg],"iso") == 0) {
error->all(FLERR,"Illegal fix bocs command. Pressure fix must be "
"cgiso . Use regular fix bocs for iso"); // MRD NJD
} else if (strcmp(arg[iarg],"cgiso") == 0) { // MRD NJD the whole else if
if (iarg+4 > narg)
error->all(FLERR,"Illegal fix bocs command. cgiso must be "
"followed by: P_0 P_f P_coupl");
p_match_flag = 1;
pcouple = XYZ;
p_start[0] = p_start[1] = p_start[2] =
utils::numeric(FLERR,arg[iarg+1],false,lmp);
p_stop[0] = p_stop[1] = p_stop[2] =
utils::numeric(FLERR,arg[iarg+2],false,lmp);
p_period[0] = p_period[1] = p_period[2] =
utils::numeric(FLERR,arg[iarg+3],false,lmp);
p_flag[0] = p_flag[1] = p_flag[2] = 1;
p_flag[3] = p_flag[4] = p_flag[5] = 0; // MRD
if (dimension == 2) { // Later I force 3D... MRD
p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0;
}
iarg += 4;
if (strcmp(arg[iarg], "analytic") == 0) {
if (iarg + 4 > narg) {
error->all(FLERR,"Illegal fix bocs command. Basis type analytic"
" must be followed by: avg_vol n_mol n_pmatch_coeff");
}
p_basis_type = BASIS_ANALYTIC;
vavg = utils::numeric(FLERR,arg[iarg+1],false,lmp);
N_mol = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
N_p_match = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
p_match_coeffs = (double *) (calloc(N_p_match, sizeof(double)) );
iarg += 4;
if (iarg + N_p_match > narg)
error->all(FLERR,"Illegal fix bocs command. Missing coeffs.");
for (int pmatchi = 0; pmatchi < N_p_match; pmatchi++)
p_match_coeffs[pmatchi] = utils::numeric(FLERR,arg[iarg+pmatchi],false,lmp);
iarg += (N_p_match);
} else if (strcmp(arg[iarg], "linear_spline") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bocs command. "
"Supply a file name after linear_spline.");
p_basis_type = BASIS_LINEAR_SPLINE;
spline_length = read_F_table( arg[iarg+1], p_basis_type );
iarg += 2;
} else if (strcmp(arg[iarg], "cubic_spline") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bocs command. "
"Supply a file name after cubic_spline.");
p_basis_type = BASIS_CUBIC_SPLINE;
spline_length = read_F_table( arg[iarg+1], p_basis_type );
iarg += 2;
} else {
std::string errmsg = fmt::format("CG basis type {} is not recognized\nSupported "
"basis types: analytic linear_spline cubic_spline",arg[iarg]);
error->all(FLERR,errmsg);
} // END NJD MRD
} else if (strcmp(arg[iarg],"tchain") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bocs command");
mtchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
// used by FixNVTSllod to preserve non-default value
mtchain_default_flag = 0;
if (mtchain < 1) error->all(FLERR,"Illegal fix bocs command");
iarg += 2;
} else if (strcmp(arg[iarg],"pchain") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bocs command");
mpchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
if (mpchain < 0) error->all(FLERR,"Illegal fix bocs command");
iarg += 2;
} else if (strcmp(arg[iarg],"mtk") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bocs command");
mtk_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"tloop") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bocs command");
nc_tchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
if (nc_tchain < 0) error->all(FLERR,"Illegal fix bocs command");
iarg += 2;
} else if (strcmp(arg[iarg],"ploop") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bocs command");
nc_pchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
if (nc_pchain < 0) error->all(FLERR,"Illegal fix bocs command");
iarg += 2;
} else {
std::string errmsg = fmt::format("Illegal fix bocs command: unrecognized keyword {}",
arg[iarg]);
error->all(FLERR,errmsg);
}
}
// error checks
if (dimension != 3) // MRD NJD
error->all(FLERR,"Invalid fix bocs command. Must use 3 dimensions");
// With cgiso, p_flag[0] = p_flag[1] = p_flag[2] = 1
// require periodicity in tensile dimension
if (p_flag[0] && domain->xperiodic == 0)
error->all(FLERR,"Cannot use fix bocs on a non-periodic dimension");
if (p_flag[1] && domain->yperiodic == 0)
error->all(FLERR,"Cannot use fix bocs on a non-periodic dimension");
if (p_flag[2] && domain->zperiodic == 0)
error->all(FLERR,"Cannot use fix bocs on a non-periodic dimension");
if (dipole_flag) {
if (!atom->sphere_flag)
error->all(FLERR,"Using update dipole flag requires atom style sphere");
if (!atom->mu_flag)
error->all(FLERR,"Using update dipole flag requires atom attribute mu");
}
if ((tstat_flag && t_period <= 0.0) ||
(p_flag[0] && p_period[0] <= 0.0) ||
(p_flag[1] && p_period[1] <= 0.0) ||
(p_flag[2] && p_period[2] <= 0.0) ||
(p_flag[3] && p_period[3] <= 0.0) ||
(p_flag[4] && p_period[4] <= 0.0) ||
(p_flag[5] && p_period[5] <= 0.0))
error->all(FLERR,"Fix bocs damping parameters must be > 0.0");
// set pstat_flag and box change and restart_pbc variables
pre_exchange_flag = 0;
pstat_flag = 0;
pstyle = ISO;
for (int i = 0; i < 6; i++)
if (p_flag[i]) pstat_flag = 1;
if (pstat_flag) {
if (p_flag[0]) box_change |= BOX_CHANGE_X;
if (p_flag[1]) box_change |= BOX_CHANGE_Y;
if (p_flag[2]) box_change |= BOX_CHANGE_Z;
if (p_flag[3]) box_change |= BOX_CHANGE_YZ;
if (p_flag[4]) box_change |= BOX_CHANGE_XZ;
if (p_flag[5]) box_change |= BOX_CHANGE_XY;
no_change_box = 1;
if (allremap == 0) restart_pbc = 1;
pstyle = ISO; // MRD this is the only one that can happen
// pre_exchange only required if flips can occur due to shape changes
if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5]))
pre_exchange_flag = pre_exchange_migrate = 1;
if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 ||
domain->xy != 0.0))
pre_exchange_flag = pre_exchange_migrate = 1;
}
// convert input periods to frequencies
t_freq = 0.0;
p_freq[0] = p_freq[1] = p_freq[2] = p_freq[3] = p_freq[4] = p_freq[5] = 0.0;
if (tstat_flag) t_freq = 1.0 / t_period;
if (p_flag[0]) p_freq[0] = 1.0 / p_period[0];
if (p_flag[1]) p_freq[1] = 1.0 / p_period[1];
if (p_flag[2]) p_freq[2] = 1.0 / p_period[2];
if (p_flag[3]) p_freq[3] = 1.0 / p_period[3];
if (p_flag[4]) p_freq[4] = 1.0 / p_period[4];
if (p_flag[5]) p_freq[5] = 1.0 / p_period[5];
// Nose/Hoover temp and pressure init
size_vector = 0;
if (tstat_flag) {
int ich;
eta = new double[mtchain];
// add one extra dummy thermostat, set to zero
eta_dot = new double[mtchain+1];
eta_dot[mtchain] = 0.0;
eta_dotdot = new double[mtchain];
for (ich = 0; ich < mtchain; ich++) {
eta[ich] = eta_dot[ich] = eta_dotdot[ich] = 0.0;
}
eta_mass = new double[mtchain];
size_vector += 2*2*mtchain;
}
if (pstat_flag) {
omega[0] = omega[1] = omega[2] = 0.0;
omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0;
omega_mass[0] = omega_mass[1] = omega_mass[2] = 0.0;
omega[3] = omega[4] = omega[5] = 0.0;
omega_dot[3] = omega_dot[4] = omega_dot[5] = 0.0;
omega_mass[3] = omega_mass[4] = omega_mass[5] = 0.0;
if (pstyle == ISO) size_vector += 2*2*1;
else if (pstyle == ANISO) size_vector += 2*2*3;
else if (pstyle == TRICLINIC) size_vector += 2*2*6;
if (mpchain) {
int ich;
etap = new double[mpchain];
// add one extra dummy thermostat, set to zero
etap_dot = new double[mpchain+1];
etap_dot[mpchain] = 0.0;
etap_dotdot = new double[mpchain];
for (ich = 0; ich < mpchain; ich++) {
etap[ich] = etap_dot[ich] =
etap_dotdot[ich] = 0.0;
}
etap_mass = new double[mpchain];
size_vector += 2*2*mpchain;
}
if (deviatoric_flag) size_vector += 1;
}
nrigid = 0;
rfix = nullptr;
if (pre_exchange_flag) irregular = new Irregular(lmp);
else irregular = nullptr;
// initialize vol0,t0 to zero to signal uninitialized
// values then assigned in init(), if necessary
vol0 = t0 = 0.0;
/*~ MRD I copied this from fix_npt.cpp 8/17/17 ~*/
if (!tstat_flag)
error->all(FLERR,"Temperature control must be used with fix bocs");
if (!pstat_flag)
error->all(FLERR,"Pressure control must be used with fix bocs");
// 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
id_temp = utils::strdup(std::string(id)+"_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
id_press = utils::strdup(std::string(id)+"_press");
modify->add_compute(fmt::format("{} all PRESSURE/BOCS {}",id_press,id_temp));
pcomputeflag = 1;
/*~ MRD End of stuff copied from fix_npt.cpp~*/
}
/* ---------------------------------------------------------------------- */
FixBocs::~FixBocs()
{
if (copymode) return;
delete [] id_dilate;
delete [] rfix;
delete irregular;
// delete temperature and pressure if fix created them
if (tcomputeflag) modify->delete_compute(id_temp);
delete [] id_temp;
if (tstat_flag) {
delete [] eta;
delete [] eta_dot;
delete [] eta_dotdot;
delete [] eta_mass;
}
if (pstat_flag) {
if (pcomputeflag) modify->delete_compute(id_press);
delete [] id_press;
if (mpchain) {
delete [] etap;
delete [] etap_dot;
delete [] etap_dotdot;
delete [] etap_mass;
}
}
if (p_match_coeffs) free(p_match_coeffs);
// Free splines memory structure
if (splines != nullptr) {
memory->destroy(splines);
spline_length = 0;
}
}
/* ---------------------------------------------------------------------- */
int FixBocs::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
if (pre_exchange_flag) mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixBocs::init()
{
// recheck that dilate group has not been deleted
if (allremap == 0) {
int idilate = group->find(id_dilate);
if (idilate == -1)
error->all(FLERR,"Fix bocs dilate group ID does not exist");
dilate_group_bit = group->bitmask[idilate];
}
// ensure no conflict with fix deform
if (pstat_flag)
{
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"deform") == 0) {
int *dimflag = (dynamic_cast<FixDeform *>(modify->fix[i]))->dimflag;
if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) ||
(p_flag[2] && dimflag[2]) || (p_flag[3] && dimflag[3]) ||
(p_flag[4] && dimflag[4]) || (p_flag[5] && dimflag[5]))
error->all(FLERR,"Cannot use fix bocs and fix deform on "
"same component of stress tensor");
}
}
// set temperature and pressure ptrs
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Temperature ID for fix bocs does not exist");
temperature = modify->compute[icompute];
if (temperature->tempbias) which = BIAS;
else which = NOBIAS;
if (pstat_flag) {
icompute = modify->find_compute(id_press);
if (icompute < 0)
error->all(FLERR,"Pressure ID for fix bocs does not exist");
pressure = modify->compute[icompute];
}
if (pstat_flag)
{
if (p_match_flag) // MRD NJD
{
if (pressure)
{
if (p_basis_type == BASIS_ANALYTIC)
{
(dynamic_cast<ComputePressureBocs *>(pressure))->send_cg_info(p_basis_type,
N_p_match, p_match_coeffs, N_mol, vavg);
}
else if (p_basis_type == BASIS_LINEAR_SPLINE || p_basis_type == BASIS_CUBIC_SPLINE)
{
(dynamic_cast<ComputePressureBocs *>(pressure))->send_cg_info(p_basis_type,
splines, spline_length);
}
}
else
{
error->all(FLERR,"Unable to find pressure. Are you sure you included"
" the compute bocsPress and fix_modify commands?");
}
}
}
// set timesteps and frequencies
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dthalf = 0.5 * update->dt;
dt4 = 0.25 * update->dt;
dt8 = 0.125 * update->dt;
dto = dthalf;
p_freq_max = 0.0;
if (pstat_flag) {
p_freq_max = MAX(p_freq[0],p_freq[1]);
p_freq_max = MAX(p_freq_max,p_freq[2]);
if (pstyle == TRICLINIC) {
p_freq_max = MAX(p_freq_max,p_freq[3]);
p_freq_max = MAX(p_freq_max,p_freq[4]);
p_freq_max = MAX(p_freq_max,p_freq[5]);
}
pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
}
if (tstat_flag)
tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
// tally the number of dimensions that are barostatted
// set initial volume and reference cell, if not already done
if (pstat_flag) {
pdim = p_flag[0] + p_flag[1] + p_flag[2];
if (vol0 == 0.0) {
if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd;
else vol0 = domain->xprd * domain->yprd;
h0_inv[0] = domain->h_inv[0];
h0_inv[1] = domain->h_inv[1];
h0_inv[2] = domain->h_inv[2];
h0_inv[3] = domain->h_inv[3];
h0_inv[4] = domain->h_inv[4];
h0_inv[5] = domain->h_inv[5];
}
}
boltz = force->boltz;
nktv2p = force->nktv2p;
if (force->kspace) kspace_flag = 1;
else kspace_flag = 0;
if (utils::strmatch(update->integrate_style,"^respa")) {
nlevels_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels;
step_respa = (dynamic_cast<Respa *>(update->integrate))->step;
dto = 0.5*step_respa[0];
}
// 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 = nullptr;
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;
}
}
// NJD MRD 2 functions
int FixBocs::read_F_table( char *filename, int p_basis_type )
{
std::string message;
double **data;
bool badInput = false;
int numEntries = 0;
FILE *fpi = fopen(filename,"r");
if (fpi) {
// Old code read the input file twice. Now we simply
// read all the lines from the input file into a string vector,
// then work with the data in-memory rather than do a second pass
// through the file.
// NB: LAMMPS coding guidelines prefer cstdio so we are intentionally
// foregoing reading with getline
if (comm->me == 0) {
error->message(FLERR, "INFO: About to read data file: {}", filename);
}
// Data file lines hold two floating point numbers.
// Line length we allocate should be long enough without being too long.
// 128 seems safe for a line we expect to be < 30 chars.
const int MAX_F_TABLE_LINE_LENGTH = 128;
char line[MAX_F_TABLE_LINE_LENGTH];
std::vector<std::string> inputLines;
while (fgets(line, MAX_F_TABLE_LINE_LENGTH, fpi)) {
inputLines.emplace_back(line);
}
fclose(fpi);
numEntries = inputLines.size();
if (comm->me == 0) {
error->message(FLERR, "INFO: Read {} lines from file", numEntries);
}
// Allocate memory for the two dimensional matrix
// that holds data from the input file.
memory->create(data, NUM_INPUT_DATA_COLUMNS, numEntries, "data");
double stdVolumeInterval = 0.0;
double currVolumeInterval = 0.0;
// When comparing doubles/floats, we need an Epsilon.
// The literature indicates getting this value right in the
// general case can be pretty complicated. I don't think it
// needs to be complicated here, though. At least based on the
// sample data I've seen where the volume values are fairly large.
const double volumeIntervalTolerance = 0.001;
int lineNum = 0; // this value is only for message
int numBadVolumeIntervals = 0; // count these for message
float f1, f2;
int test_sscanf;
for (int i = 0; i < (int)inputLines.size(); ++i) {
lineNum++; // count each line processed now so lineNum messages can be 1-based
test_sscanf = sscanf(inputLines.at(i).c_str()," %f , %f ",&f1, &f2);
if (test_sscanf == 2)
{
data[VOLUME][i] = (double)f1;
data[PRESSURE_CORRECTION][i] = (double)f2;
if (i == 1)
{
// second entry is used to compute the validation interval used below
stdVolumeInterval = data[VOLUME][i] - data[VOLUME][i-1];
}
else if (i > 1)
{
// after second entry, all intervals are validated
currVolumeInterval = data[VOLUME][i] - data[VOLUME][i-1];
if (fabs(currVolumeInterval - stdVolumeInterval) > volumeIntervalTolerance) {
if (comm->me == 0)
error->warning(FLERR,"Bad volume interval. Spline analysis requires uniform"
" volume distribution, found inconsistent volume"
" differential, line {} of file {}\nWARNING:\tline: {}",
lineNum, filename, inputLines.at(i));
badInput = true;
numBadVolumeIntervals++;
}
// no concluding else is intentional: i = 0, first line, no interval to validate
}
}
else
{
if (comm->me == 0)
error->warning(FLERR,"Bad input format: did not find 2 comma separated numeric"
" values in line {} of file {}\nWARNING:\tline: {}",
lineNum, filename, inputLines.at(i));
badInput = true;
}
if (badInput)
{
numBadVolumeIntervals++;
}
}
if (numBadVolumeIntervals > 0 && comm->me == 0) {
error->message(FLERR, "INFO: total number bad volume intervals = {}", numBadVolumeIntervals);
}
}
else {
error->all(FLERR,"ERROR: Unable to open file: {}", filename);
}
if (badInput && comm->me == 0) {
error->warning(FLERR,"Bad volume / pressure-correction data: {}\nSee details above", filename);
}
if (p_basis_type == BASIS_LINEAR_SPLINE)
{
spline_length = numEntries;
numEntries = build_linear_splines(data);
}
else if (p_basis_type == BASIS_CUBIC_SPLINE)
{
spline_length = numEntries;
numEntries = build_cubic_splines(data);
}
else
{
error->all(FLERR,"ERROR: invalid p_basis_type value of {} in read_F_table", p_basis_type);
}
memory->destroy(data);
return numEntries;
}
int FixBocs::build_linear_splines(double **data) {
splines = (double **) calloc(NUM_LINEAR_SPLINE_COLUMNS,sizeof(double *));
splines[VOLUME] = (double *) calloc(spline_length,sizeof(double));
splines[PRESSURE_CORRECTION] = (double *) calloc(spline_length,sizeof(double));
for (int i = 0; i < spline_length; ++i)
{
splines[VOLUME][i] = data[VOLUME][i];
splines[PRESSURE_CORRECTION][i] = data[PRESSURE_CORRECTION][i];
}
if (comm->me == 0) {
error->message(FLERR, "INFO: leaving build_linear_splines, spline_length = {}", spline_length);
}
return spline_length;
}
int FixBocs::build_cubic_splines(double **data)
{
int n = spline_length;
double *a, *b, *d, *h, *alpha, *c, *l, *mu, *z;
// 2020-07-17 ag:
// valgrind says that we read/write a[n] down in the
// for (int j=n-1; j>=0; j--) loop below
// and I agree.
// So the size of a must be n+1, not n as was found
// in the original code.
memory->create(a, n+1, "a");
memory->create(b, n+1, "b");
memory->create(c, n+1, "c");
memory->create(d, n+1, "d");
memory->create(h, n, "h");
memory->create(alpha, n, "alpha");
memory->create(l, n, "l");
memory->create(mu, n, "mu");
memory->create(z, n, "z");
for (int i=0; i<n; i++)
{
a[i] = data[1][i];
b[i] = 0.0;
d[i] = 0.0;
if (i<(n-1))
{
h[i] = (data[0][i+1] - data[0][i]);
}
double alpha_i;
if (i>1 && i<(n-1))
{
alpha_i = (3.0 / h[i]) * ( data[1][i+1] - data[1][i]) - (3.0 / h[i-1] )
* ( data[1][i] - data[1][i-1] );
alpha[i-1] = alpha_i;
}
}
l[0] = 1.0;
mu[0] = 0.0;
z[0] = 0.0;
for (int i=1; i<n-1; i++)
{
l[i] = 2*(data[0][i+1] - data[0][i-1]) - h[i-1] * mu[i-1];
mu[i] = h[i]/l[i];
z[i] = (alpha[i] - h[i-1] * z[i-1]) / l[i];
}
l[n-1] = 1.0;
mu[n-1] = 0.0;
z[n-1] = 0.0;
// 2020-07-17 ag: We've been using an uninitialized value for a[n]
// That seems like a bad idea. This may not be the right value
// but its a value.
a[n] = 0.0;
b[n] = 0.0;
c[n] = 0.0;
d[n] = 0.0;
for (int j=n-1; j>=0; j--)
{
c[j] = z[j] - mu[j]*c[j+1];
b[j] = (a[j+1]-a[j])/h[j] - h[j]*(c[j+1] + 2.0 * c[j])/3.0;
d[j] = (c[j+1]-c[j])/(3.0 * h[j]);
}
int numSplines = n - 1;
memory->create(splines, NUM_CUBIC_SPLINE_COLUMNS, numSplines, "splines");
for (int idx = 0; idx < numSplines; ++idx)
{
splines[0][idx] = data[0][idx];
splines[1][idx] = a[idx];
splines[2][idx] = b[idx];
splines[3][idx] = c[idx];
splines[4][idx] = d[idx];
}
memory->destroy(a);
memory->destroy(b);
memory->destroy(c);
memory->destroy(d);
memory->destroy(h);
memory->destroy(alpha);
memory->destroy(l);
memory->destroy(mu);
memory->destroy(z);
if (comm->me == 0) {
error->message(FLERR, "INFO: leaving build_cubic_splines, numSplines = {}", numSplines);
}
// Tell the caller how many splines we created
return numSplines;
}
// END NJD MRD 2 functions
/* ----------------------------------------------------------------------
compute T,P before integrator starts
------------------------------------------------------------------------- */
void FixBocs::setup(int /*vflag*/)
{
// tdof needed by compute_temp_target()
t_current = temperature->compute_scalar();
tdof = temperature->dof;
// t_target is needed by NVT and NPT in compute_scalar()
// If no thermostat or using fix nphug,
// t_target must be defined by other means.
if (tstat_flag && strstr(style,"nphug") == nullptr) {
compute_temp_target();
} else if (pstat_flag) {
// t0 = reference temperature for masses
// cannot be done in init() b/c temperature cannot be called there
// is b/c Modify::init() inits computes after fixes due to dof dependence
// guesstimate a unit-dependent t0 if actual T = 0.0
// if it was read in from a restart file, leave it be
if (t0 == 0.0) {
t0 = temperature->compute_scalar();
if (t0 == 0.0) {
if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0;
else t0 = 300.0;
}
}
t_target = t0;
}
if (pstat_flag) compute_press_target();
if (pstat_flag) {
if (pstyle == ISO) pressure->compute_scalar();
else pressure->compute_vector();
couple();
pressure->addstep(update->ntimestep+1);
}
// masses and initial forces on thermostat variables
if (tstat_flag) {
eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq);
for (int ich = 1; ich < mtchain; ich++)
eta_mass[ich] = boltz * t_target / (t_freq*t_freq);
for (int ich = 1; ich < mtchain; ich++) {
eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] -
boltz * t_target) / eta_mass[ich];
}
}
// masses and initial forces on barostat variables
if (pstat_flag) {
double kt = boltz * t_target;
double nkt = (atom->natoms + 1) * kt;
for (int i = 0; i < 3; i++)
if (p_flag[i])
omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
if (pstyle == TRICLINIC) {
for (int i = 3; i < 6; i++)
if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
}
// masses and initial forces on barostat thermostat variables
if (mpchain) {
etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max);
for (int ich = 1; ich < mpchain; ich++)
etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max);
for (int ich = 1; ich < mpchain; ich++)
etap_dotdot[ich] =
(etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] -
boltz * t_target) / etap_mass[ich];
}
}
}
/* ----------------------------------------------------------------------
1st half of Verlet update
------------------------------------------------------------------------- */
void FixBocs::initial_integrate(int /*vflag*/)
{
// update eta_press_dot
if (pstat_flag && mpchain) nhc_press_integrate();
// update eta_dot
if (tstat_flag) {
compute_temp_target();
nhc_temp_integrate();
}
// need to recompute pressure to account for change in KE
// t_current is up-to-date, but compute_temperature is not
// compute appropriately coupled elements of mvv_current
if (pstat_flag) {
if (pstyle == ISO) {
temperature->compute_scalar();
pressure->compute_scalar();
} else {
temperature->compute_vector();
pressure->compute_vector();
}
couple();
pressure->addstep(update->ntimestep+1);
}
if (pstat_flag) {
compute_press_target();
nh_omega_dot();
nh_v_press();
}
nve_v();
// remap simulation box by 1/2 step
if (pstat_flag) remap();
nve_x();
// remap simulation box by 1/2 step
// redo KSpace coeffs since volume has changed
if (pstat_flag) {
remap();
if (kspace_flag) force->kspace->setup();
}
}
/* ----------------------------------------------------------------------
2nd half of Verlet update
------------------------------------------------------------------------- */
void FixBocs::final_integrate()
{
nve_v();
// re-compute temp before nh_v_press()
// only needed for temperature computes with BIAS on reneighboring steps:
// b/c some biases store per-atom values (e.g. temp/profile)
// per-atom values are invalid if reneigh/comm occurred
// since temp->compute() in initial_integrate()
if (which == BIAS && neighbor->ago == 0)
t_current = temperature->compute_scalar();
if (pstat_flag) nh_v_press();
// compute new T,P after velocities rescaled by nh_v_press()
// compute appropriately coupled elements of mvv_current
t_current = temperature->compute_scalar();
tdof = temperature->dof;
if (pstat_flag) {
if (pstyle == ISO) pressure->compute_scalar();
else {
temperature->compute_vector();
pressure->compute_vector();
}
couple();
pressure->addstep(update->ntimestep+1);
}
if (pstat_flag) nh_omega_dot();
// update eta_dot
// update eta_press_dot
if (tstat_flag) nhc_temp_integrate();
if (pstat_flag && mpchain) nhc_press_integrate();
}
/* ---------------------------------------------------------------------- */
void FixBocs::initial_integrate_respa(int /*vflag*/, int ilevel, int /*iloop*/)
{
// set timesteps by level
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dthalf = 0.5 * step_respa[ilevel];
// outermost level - update eta_dot and omega_dot, apply to v
// all other levels - NVE update of v
// x,v updates only performed for atoms in group
if (ilevel == nlevels_respa-1) {
// update eta_press_dot
if (pstat_flag && mpchain) nhc_press_integrate();
// update eta_dot
if (tstat_flag) {
compute_temp_target();
nhc_temp_integrate();
}
// recompute pressure to account for change in KE
// t_current is up-to-date, but compute_temperature is not
// compute appropriately coupled elements of mvv_current
if (pstat_flag) {
if (pstyle == ISO) {
temperature->compute_scalar();
pressure->compute_scalar();
} else {
temperature->compute_vector();
pressure->compute_vector();
}
couple();
pressure->addstep(update->ntimestep+1);
}
if (pstat_flag) {
compute_press_target();
nh_omega_dot();
nh_v_press();
}
nve_v();
} else nve_v();
// innermost level - also update x only for atoms in group
// if barostat, perform 1/2 step remap before and after
if (ilevel == 0) {
if (pstat_flag) remap();
nve_x();
if (pstat_flag) remap();
}
}
/* ---------------------------------------------------------------------- */
void FixBocs::pre_force_respa(int /*vflag*/, int ilevel, int /*iloop*/)
{
// if barostat, redo KSpace coeffs at outermost level,
// since volume has changed
if (ilevel == nlevels_respa-1 && kspace_flag && pstat_flag)
force->kspace->setup();
}
/* ---------------------------------------------------------------------- */
void FixBocs::final_integrate_respa(int ilevel, int /*iloop*/)
{
// set timesteps by level
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dthalf = 0.5 * step_respa[ilevel];
// outermost level - update eta_dot and omega_dot, apply via final_integrate
// all other levels - NVE update of v
if (ilevel == nlevels_respa-1) final_integrate();
else nve_v();
}
/* ---------------------------------------------------------------------- */
void FixBocs::couple()
{
double *tensor = pressure->vector;
if (pstyle == ISO)
p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
else if (pcouple == XYZ) {
double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]);
p_current[0] = p_current[1] = p_current[2] = ave;
} else if (pcouple == XY) {
double ave = 0.5 * (tensor[0] + tensor[1]);
p_current[0] = p_current[1] = ave;
p_current[2] = tensor[2];
} else if (pcouple == YZ) {
double ave = 0.5 * (tensor[1] + tensor[2]);
p_current[1] = p_current[2] = ave;
p_current[0] = tensor[0];
} else if (pcouple == XZ) {
double ave = 0.5 * (tensor[0] + tensor[2]);
p_current[0] = p_current[2] = ave;
p_current[1] = tensor[1];
} else {
p_current[0] = tensor[0];
p_current[1] = tensor[1];
p_current[2] = tensor[2];
}
if (!std::isfinite(p_current[0]) || !std::isfinite(p_current[1]) || !std::isfinite(p_current[2]))
error->all(FLERR,"Non-numeric pressure - simulation unstable");
// switch order from xy-xz-yz to Voigt
if (pstyle == TRICLINIC) {
p_current[3] = tensor[5];
p_current[4] = tensor[4];
p_current[5] = tensor[3];
if (!std::isfinite(p_current[3]) || !std::isfinite(p_current[4]) || !std::isfinite(p_current[5]))
error->all(FLERR,"Non-numeric pressure - simulation unstable");
}
}
/* ----------------------------------------------------------------------
change box size
remap all atoms or dilate group atoms depending on allremap flag
if rigid bodies exist, scale rigid body centers-of-mass
------------------------------------------------------------------------- */
void FixBocs::remap()
{
int i;
double oldlo,oldhi;
double expfac;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *h = domain->h;
// omega is not used, except for book-keeping
for (int i = 0; i < 6; i++) omega[i] += dto*omega_dot[i];
// convert pertinent atoms and rigid bodies to lamda coords
if (allremap) domain->x2lamda(nlocal);
else {
for (i = 0; i < nlocal; i++)
if (mask[i] & dilate_group_bit)
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
// this operation corresponds to applying the
// translate and scale operations
// corresponding to the solution of the following ODE:
//
// h_dot = omega_dot * h
//
// where h_dot, omega_dot and h are all upper-triangular
// 3x3 tensors. In Voigt notation, the elements of the
// RHS product tensor are:
// h_dot = [0*0, 1*1, 2*2, 1*3+3*2, 0*4+5*3+4*2, 0*5+5*1]
//
// Ordering of operations preserves time symmetry.
double dto2 = dto/2.0;
double dto4 = dto/4.0;
double dto8 = dto/8.0;
// off-diagonal components, first half
if (pstyle == TRICLINIC) {
if (p_flag[4]) {
expfac = exp(dto8*omega_dot[0]);
h[4] *= expfac;
h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
h[4] *= expfac;
}
if (p_flag[3]) {
expfac = exp(dto4*omega_dot[1]);
h[3] *= expfac;
h[3] += dto2*(omega_dot[3]*h[2]);
h[3] *= expfac;
}
if (p_flag[5]) {
expfac = exp(dto4*omega_dot[0]);
h[5] *= expfac;
h[5] += dto2*(omega_dot[5]*h[1]);
h[5] *= expfac;
}
if (p_flag[4]) {
expfac = exp(dto8*omega_dot[0]);
h[4] *= expfac;
h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
h[4] *= expfac;
}
}
// scale diagonal components
// scale tilt factors with cell, if set
if (p_flag[0]) {
oldlo = domain->boxlo[0];
oldhi = domain->boxhi[0];
expfac = exp(dto*omega_dot[0]);
domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0];
domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0];
}
if (p_flag[1]) {
oldlo = domain->boxlo[1];
oldhi = domain->boxhi[1];
expfac = exp(dto*omega_dot[1]);
domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1];
domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1];
if (scalexy) h[5] *= expfac;
}
if (p_flag[2]) {
oldlo = domain->boxlo[2];
oldhi = domain->boxhi[2];
expfac = exp(dto*omega_dot[2]);
domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2];
domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2];
if (scalexz) h[4] *= expfac;
if (scaleyz) h[3] *= expfac;
}
// off-diagonal components, second half
if (pstyle == TRICLINIC) {
if (p_flag[4]) {
expfac = exp(dto8*omega_dot[0]);
h[4] *= expfac;
h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
h[4] *= expfac;
}
if (p_flag[3]) {
expfac = exp(dto4*omega_dot[1]);
h[3] *= expfac;
h[3] += dto2*(omega_dot[3]*h[2]);
h[3] *= expfac;
}
if (p_flag[5]) {
expfac = exp(dto4*omega_dot[0]);
h[5] *= expfac;
h[5] += dto2*(omega_dot[5]*h[1]);
h[5] *= expfac;
}
if (p_flag[4]) {
expfac = exp(dto8*omega_dot[0]);
h[4] *= expfac;
h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
h[4] *= expfac;
}
}
domain->yz = h[3];
domain->xz = h[4];
domain->xy = h[5];
// tilt factor to cell length ratio can not exceed TILTMAX in one step
if (domain->yz < -TILTMAX*domain->yprd ||
domain->yz > TILTMAX*domain->yprd ||
domain->xz < -TILTMAX*domain->xprd ||
domain->xz > TILTMAX*domain->xprd ||
domain->xy < -TILTMAX*domain->xprd ||
domain->xy > TILTMAX*domain->xprd)
error->all(FLERR,"Fix bocs has tilted box too far in one step - "
"periodic cell is too far from equilibrium state");
domain->set_global_box();
domain->set_local_box();
// convert pertinent atoms and rigid bodies back to box coords
if (allremap) domain->lamda2x(nlocal);
else {
for (i = 0; i < nlocal; i++)
if (mask[i] & dilate_group_bit)
domain->lamda2x(x[i],x[i]);
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(1);
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixBocs::write_restart(FILE *fp)
{
int nsize = size_restart_global();
double *list;
memory->create(list,nsize,"nh:list");
pack_restart_data(list);
if (comm->me == 0) {
int size = nsize * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(list,sizeof(double),nsize,fp);
}
memory->destroy(list);
}
/* ----------------------------------------------------------------------
calculate the number of data to be packed
------------------------------------------------------------------------- */
int FixBocs::size_restart_global()
{
int nsize = 2;
if (tstat_flag) nsize += 1 + 2*mtchain;
if (pstat_flag) {
nsize += 16 + 2*mpchain;
if (deviatoric_flag) nsize += 6;
}
return nsize;
}
/* ----------------------------------------------------------------------
pack restart data
------------------------------------------------------------------------- */
int FixBocs::pack_restart_data(double *list)
{
int n = 0;
list[n++] = tstat_flag;
if (tstat_flag) {
list[n++] = mtchain;
for (int ich = 0; ich < mtchain; ich++)
list[n++] = eta[ich];
for (int ich = 0; ich < mtchain; ich++)
list[n++] = eta_dot[ich];
}
list[n++] = pstat_flag;
if (pstat_flag) {
list[n++] = omega[0];
list[n++] = omega[1];
list[n++] = omega[2];
list[n++] = omega[3];
list[n++] = omega[4];
list[n++] = omega[5];
list[n++] = omega_dot[0];
list[n++] = omega_dot[1];
list[n++] = omega_dot[2];
list[n++] = omega_dot[3];
list[n++] = omega_dot[4];
list[n++] = omega_dot[5];
list[n++] = vol0;
list[n++] = t0;
list[n++] = mpchain;
if (mpchain) {
for (int ich = 0; ich < mpchain; ich++)
list[n++] = etap[ich];
for (int ich = 0; ich < mpchain; ich++)
list[n++] = etap_dot[ich];
}
list[n++] = deviatoric_flag;
if (deviatoric_flag) {
list[n++] = h0_inv[0];
list[n++] = h0_inv[1];
list[n++] = h0_inv[2];
list[n++] = h0_inv[3];
list[n++] = h0_inv[4];
list[n++] = h0_inv[5];
}
}
return n;
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixBocs::restart(char *buf)
{
int n = 0;
auto list = (double *) buf;
int flag = static_cast<int> (list[n++]);
if (flag) {
int m = static_cast<int> (list[n++]);
if (tstat_flag && m == mtchain) {
for (int ich = 0; ich < mtchain; ich++)
eta[ich] = list[n++];
for (int ich = 0; ich < mtchain; ich++)
eta_dot[ich] = list[n++];
} else n += 2*m;
}
flag = static_cast<int> (list[n++]);
if (flag) {
omega[0] = list[n++];
omega[1] = list[n++];
omega[2] = list[n++];
omega[3] = list[n++];
omega[4] = list[n++];
omega[5] = list[n++];
omega_dot[0] = list[n++];
omega_dot[1] = list[n++];
omega_dot[2] = list[n++];
omega_dot[3] = list[n++];
omega_dot[4] = list[n++];
omega_dot[5] = list[n++];
vol0 = list[n++];
t0 = list[n++];
int m = static_cast<int> (list[n++]);
if (pstat_flag && m == mpchain) {
for (int ich = 0; ich < mpchain; ich++)
etap[ich] = list[n++];
for (int ich = 0; ich < mpchain; ich++)
etap_dot[ich] = list[n++];
} else n+=2*m;
flag = static_cast<int> (list[n++]);
if (flag) {
h0_inv[0] = list[n++];
h0_inv[1] = list[n++];
h0_inv[2] = list[n++];
h0_inv[3] = list[n++];
h0_inv[4] = list[n++];
h0_inv[5] = list[n++];
}
}
}
/* ---------------------------------------------------------------------- */
int FixBocs::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
if (tcomputeflag) {
modify->delete_compute(id_temp);
tcomputeflag = 0;
}
delete [] id_temp;
id_temp = utils::strdup(arg[1]);
int icompute = modify->find_compute(arg[1]);
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 != 0 && comm->me == 0)
error->warning(FLERR,"Temperature for fix modify is not for group all");
// reset id_temp of pressure to new temperature ID
if (pstat_flag) {
icompute = modify->find_compute(id_press);
if (icompute < 0)
error->all(FLERR,"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(FLERR,"Illegal fix_modify command");
if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command");
if (pcomputeflag) {
modify->delete_compute(id_press);
pcomputeflag = 0;
}
delete [] id_press;
id_press = utils::strdup(arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");
pressure = modify->compute[icompute];
if (p_match_flag) // NJD MRD
{
if (p_basis_type == BASIS_ANALYTIC)
{
(dynamic_cast<ComputePressureBocs *>(pressure))->send_cg_info(p_basis_type, N_p_match,
p_match_coeffs, N_mol, vavg);
}
else if (p_basis_type == BASIS_LINEAR_SPLINE || p_basis_type == BASIS_CUBIC_SPLINE )
{
(dynamic_cast<ComputePressureBocs *>(pressure))->send_cg_info(p_basis_type, splines, spline_length );
}
}
if (pressure->pressflag == 0)
{
error->all(FLERR, "Fix_modify pressure ID does not compute pressure");
}
return 2;
}
return 0;
}
/* ---------------------------------------------------------------------- */
double FixBocs::compute_scalar()
{
int i;
double volume;
double energy;
double kt = boltz * t_target;
double lkt_press = 0.0;
int ich;
if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
else volume = domain->xprd * domain->yprd;
energy = 0.0;
// thermostat chain energy is equivalent to Eq. (2) in
// Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117
// Sum(0.5*p_eta_k^2/Q_k,k=1,M) + L*k*T*eta_1 + Sum(k*T*eta_k,k=2,M),
// where L = tdof
// M = mtchain
// p_eta_k = Q_k*eta_dot[k-1]
// Q_1 = L*k*T/t_freq^2
// Q_k = k*T/t_freq^2, k > 1
if (tstat_flag) {
energy += ke_target * eta[0] + 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0];
for (ich = 1; ich < mtchain; ich++)
energy += kt * eta[ich] + 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich];
}
// barostat energy is equivalent to Eq. (8) in
// Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117
// Sum(0.5*p_omega^2/W + P*V),
// where N = natoms
// p_omega = W*omega_dot
// W = N*k*T/p_freq^2
// sum is over barostatted dimensions
if (pstat_flag) {
for (i = 0; i < 3; i++) {
if (p_flag[i]) {
energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i] +
p_hydro*(volume-vol0) / (pdim*nktv2p);
lkt_press += kt;
}
}
if (pstyle == TRICLINIC) {
for (i = 3; i < 6; i++) {
if (p_flag[i]) {
energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i];
lkt_press += kt;
}
}
}
// extra contributions from thermostat chain for barostat
if (mpchain) {
energy += lkt_press * etap[0] + 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0];
for (ich = 1; ich < mpchain; ich++)
energy += kt * etap[ich] +
0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich];
}
// extra contribution from strain energy
if (deviatoric_flag) energy += compute_strain_energy();
}
return energy;
}
/* ----------------------------------------------------------------------
return a single element of the following vectors, in this order:
eta[tchain], eta_dot[tchain], omega[ndof], omega_dot[ndof]
etap[pchain], etap_dot[pchain], PE_eta[tchain], KE_eta_dot[tchain]
PE_omega[ndof], KE_omega_dot[ndof], PE_etap[pchain], KE_etap_dot[pchain]
PE_strain[1]
if no thermostat exists, related quantities are omitted from the list
if no barostat exists, related quantities are omitted from the list
ndof = 1,3,6 degrees of freedom for pstyle = ISO,ANISO,TRI
------------------------------------------------------------------------- */
double FixBocs::compute_vector(int n)
{
int ilen;
if (tstat_flag) {
ilen = mtchain;
if (n < ilen) return eta[n];
n -= ilen;
ilen = mtchain;
if (n < ilen) return eta_dot[n];
n -= ilen;
}
if (pstat_flag) {
if (pstyle == ISO) {
ilen = 1;
if (n < ilen) return omega[n];
n -= ilen;
} else if (pstyle == ANISO) {
ilen = 3;
if (n < ilen) return omega[n];
n -= ilen;
} else {
ilen = 6;
if (n < ilen) return omega[n];
n -= ilen;
}
if (pstyle == ISO) {
ilen = 1;
if (n < ilen) return omega_dot[n];
n -= ilen;
} else if (pstyle == ANISO) {
ilen = 3;
if (n < ilen) return omega_dot[n];
n -= ilen;
} else {
ilen = 6;
if (n < ilen) return omega_dot[n];
n -= ilen;
}
if (mpchain) {
ilen = mpchain;
if (n < ilen) return etap[n];
n -= ilen;
ilen = mpchain;
if (n < ilen) return etap_dot[n];
n -= ilen;
}
}
double volume;
double kt = boltz * t_target;
double lkt_press = kt;
int ich;
if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
else volume = domain->xprd * domain->yprd;
if (tstat_flag) {
ilen = mtchain;
if (n < ilen) {
ich = n;
if (ich == 0)
return ke_target * eta[0];
else
return kt * eta[ich];
}
n -= ilen;
ilen = mtchain;
if (n < ilen) {
ich = n;
if (ich == 0)
return 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0];
else
return 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich];
}
n -= ilen;
}
if (pstat_flag) {
if (pstyle == ISO) {
ilen = 1;
if (n < ilen)
return p_hydro*(volume-vol0) / nktv2p;
n -= ilen;
} else if (pstyle == ANISO) {
ilen = 3;
if (n < ilen) {
if (p_flag[n])
return p_hydro*(volume-vol0) / (pdim*nktv2p);
else
return 0.0;
}
n -= ilen;
} else {
ilen = 6;
if (n < ilen) {
if (n > 2) return 0.0;
else if (p_flag[n])
return p_hydro*(volume-vol0) / (pdim*nktv2p);
else
return 0.0;
}
n -= ilen;
}
if (pstyle == ISO) {
ilen = 1;
if (n < ilen)
return pdim*0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
n -= ilen;
} else if (pstyle == ANISO) {
ilen = 3;
if (n < ilen) {
if (p_flag[n])
return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
else return 0.0;
}
n -= ilen;
} else {
ilen = 6;
if (n < ilen) {
if (p_flag[n])
return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
else return 0.0;
}
n -= ilen;
}
if (mpchain) {
ilen = mpchain;
if (n < ilen) {
ich = n;
if (ich == 0) return lkt_press * etap[0];
else return kt * etap[ich];
}
n -= ilen;
ilen = mpchain;
if (n < ilen) {
ich = n;
if (ich == 0)
return 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0];
else
return 0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich];
}
n -= ilen;
}
if (deviatoric_flag) {
ilen = 1;
if (n < ilen)
return compute_strain_energy();
n -= ilen;
}
}
return 0.0;
}
/* ---------------------------------------------------------------------- */
void FixBocs::reset_target(double t_new)
{
t_target = t_start = t_stop = t_new;
}
/* ---------------------------------------------------------------------- */
void FixBocs::reset_dt()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dthalf = 0.5 * update->dt;
dt4 = 0.25 * update->dt;
dt8 = 0.125 * update->dt;
dto = dthalf;
// If using respa, then remap is performed in innermost level
if (utils::strmatch(update->integrate_style,"^respa"))
dto = 0.5*step_respa[0];
if (pstat_flag)
pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
if (tstat_flag)
tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
}
/* ----------------------------------------------------------------------
extract thermostat properties
------------------------------------------------------------------------- */
void *FixBocs::extract(const char *str, int &dim)
{
dim=0;
if (tstat_flag && strcmp(str,"t_target") == 0) {
return &t_target;
} else if (tstat_flag && strcmp(str,"t_start") == 0) {
return &t_start;
} else if (tstat_flag && strcmp(str,"t_stop") == 0) {
return &t_stop;
} else if (tstat_flag && strcmp(str,"mtchain") == 0) {
return &mtchain;
} else if (pstat_flag && strcmp(str,"mpchain") == 0) {
return &mtchain;
}
dim=1;
if (tstat_flag && strcmp(str,"eta") == 0) {
return &eta;
} else if (pstat_flag && strcmp(str,"etap") == 0) {
return &eta;
} else if (pstat_flag && strcmp(str,"p_flag") == 0) {
return &p_flag;
} else if (pstat_flag && strcmp(str,"p_start") == 0) {
return &p_start;
} else if (pstat_flag && strcmp(str,"p_stop") == 0) {
return &p_stop;
} else if (pstat_flag && strcmp(str,"p_target") == 0) {
return &p_target;
}
return nullptr;
}
/* ----------------------------------------------------------------------
perform half-step update of chain thermostat variables
------------------------------------------------------------------------- */
void FixBocs::nhc_temp_integrate()
{
int ich;
double expfac;
double kecurrent = tdof * boltz * t_current;
// Update masses, to preserve initial freq, if flag set
if (eta_mass_flag) {
eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq);
for (int ich = 1; ich < mtchain; ich++)
eta_mass[ich] = boltz * t_target / (t_freq*t_freq);
}
if (eta_mass[0] > 0.0)
eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0];
else eta_dotdot[0] = 0.0;
double ncfac = 1.0/nc_tchain;
for (int iloop = 0; iloop < nc_tchain; iloop++) {
for (ich = mtchain-1; ich > 0; ich--) {
expfac = exp(-ncfac*dt8*eta_dot[ich+1]);
eta_dot[ich] *= expfac;
eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4;
eta_dot[ich] *= tdrag_factor;
eta_dot[ich] *= expfac;
}
expfac = exp(-ncfac*dt8*eta_dot[1]);
eta_dot[0] *= expfac;
eta_dot[0] += eta_dotdot[0] * ncfac*dt4;
eta_dot[0] *= tdrag_factor;
eta_dot[0] *= expfac;
factor_eta = exp(-ncfac*dthalf*eta_dot[0]);
nh_v_temp();
// rescale temperature due to velocity scaling
// should not be necessary to explicitly recompute the temperature
t_current *= factor_eta*factor_eta;
kecurrent = tdof * boltz * t_current;
if (eta_mass[0] > 0.0)
eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0];
else eta_dotdot[0] = 0.0;
for (ich = 0; ich < mtchain; ich++)
eta[ich] += ncfac*dthalf*eta_dot[ich];
eta_dot[0] *= expfac;
eta_dot[0] += eta_dotdot[0] * ncfac*dt4;
eta_dot[0] *= expfac;
for (ich = 1; ich < mtchain; ich++) {
expfac = exp(-ncfac*dt8*eta_dot[ich+1]);
eta_dot[ich] *= expfac;
eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1]
- boltz * t_target)/eta_mass[ich];
eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4;
eta_dot[ich] *= expfac;
}
}
}
/* ----------------------------------------------------------------------
perform half-step update of chain thermostat variables for barostat
scale barostat velocities
------------------------------------------------------------------------- */
void FixBocs::nhc_press_integrate()
{
int ich,i,pdof;
double expfac,factor_etap,kecurrent;
double kt = boltz * t_target;
double lkt_press;
// Update masses, to preserve initial freq, if flag set
if (omega_mass_flag) {
double nkt = (atom->natoms + 1) * kt;
for (int i = 0; i < 3; i++)
if (p_flag[i])
omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
if (pstyle == TRICLINIC) {
for (int i = 3; i < 6; i++)
if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
}
}
if (etap_mass_flag) {
if (mpchain) {
etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max);
for (int ich = 1; ich < mpchain; ich++)
etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max);
for (int ich = 1; ich < mpchain; ich++)
etap_dotdot[ich] =
(etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] -
boltz * t_target) / etap_mass[ich];
}
}
kecurrent = 0.0;
pdof = 0;
for (i = 0; i < 3; i++) {
if (p_flag[i]) {
kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
pdof++;
}
}
if (pstyle == TRICLINIC) {
for (i = 3; i < 6; i++) {
if (p_flag[i]) {
kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
pdof++;
}
}
}
if (pstyle == ISO) lkt_press = kt;
else lkt_press = pdof * kt;
etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
double ncfac = 1.0/nc_pchain;
for (int iloop = 0; iloop < nc_pchain; iloop++) {
for (ich = mpchain-1; ich > 0; ich--) {
expfac = exp(-ncfac*dt8*etap_dot[ich+1]);
etap_dot[ich] *= expfac;
etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4;
etap_dot[ich] *= pdrag_factor;
etap_dot[ich] *= expfac;
}
expfac = exp(-ncfac*dt8*etap_dot[1]);
etap_dot[0] *= expfac;
etap_dot[0] += etap_dotdot[0] * ncfac*dt4;
etap_dot[0] *= pdrag_factor;
etap_dot[0] *= expfac;
for (ich = 0; ich < mpchain; ich++)
etap[ich] += ncfac*dthalf*etap_dot[ich];
factor_etap = exp(-ncfac*dthalf*etap_dot[0]);
for (i = 0; i < 3; i++)
if (p_flag[i]) omega_dot[i] *= factor_etap;
if (pstyle == TRICLINIC) {
for (i = 3; i < 6; i++)
if (p_flag[i]) omega_dot[i] *= factor_etap;
}
kecurrent = 0.0;
for (i = 0; i < 3; i++)
if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
if (pstyle == TRICLINIC) {
for (i = 3; i < 6; i++)
if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
}
etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
etap_dot[0] *= expfac;
etap_dot[0] += etap_dotdot[0] * ncfac*dt4;
etap_dot[0] *= expfac;
for (ich = 1; ich < mpchain; ich++) {
expfac = exp(-ncfac*dt8*etap_dot[ich+1]);
etap_dot[ich] *= expfac;
etap_dotdot[ich] =
(etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - boltz*t_target) /
etap_mass[ich];
etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4;
etap_dot[ich] *= expfac;
}
}
}
/* ----------------------------------------------------------------------
perform half-step barostat scaling of velocities
-----------------------------------------------------------------------*/
void FixBocs::nh_v_press()
{
double factor[3];
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
factor[0] = exp(-dt4*(omega_dot[0]+mtk_term2));
factor[1] = exp(-dt4*(omega_dot[1]+mtk_term2));
factor[2] = exp(-dt4*(omega_dot[2]+mtk_term2));
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
v[i][0] *= factor[0];
v[i][1] *= factor[1];
v[i][2] *= factor[2];
if (pstyle == TRICLINIC) {
v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]);
v[i][1] += -dthalf*v[i][2]*omega_dot[3];
}
v[i][0] *= factor[0];
v[i][1] *= factor[1];
v[i][2] *= factor[2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
v[i][0] *= factor[0];
v[i][1] *= factor[1];
v[i][2] *= factor[2];
if (pstyle == TRICLINIC) {
v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]);
v[i][1] += -dthalf*v[i][2]*omega_dot[3];
}
v[i][0] *= factor[0];
v[i][1] *= factor[1];
v[i][2] *= factor[2];
temperature->restore_bias(i,v[i]);
}
}
}
}
/* ----------------------------------------------------------------------
perform half-step update of velocities
-----------------------------------------------------------------------*/
void FixBocs::nve_v()
{
double dtfm;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
}
}
/* ----------------------------------------------------------------------
perform full-step update of positions
-----------------------------------------------------------------------*/
void FixBocs::nve_x()
{
double **x = atom->x;
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
// x update by full step only for atoms in group
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
}
/* ----------------------------------------------------------------------
perform half-step thermostat scaling of velocities
-----------------------------------------------------------------------*/
void FixBocs::nh_v_temp()
{
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
v[i][0] *= factor_eta;
v[i][1] *= factor_eta;
v[i][2] *= factor_eta;
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
v[i][0] *= factor_eta;
v[i][1] *= factor_eta;
v[i][2] *= factor_eta;
temperature->restore_bias(i,v[i]);
}
}
}
}
/* ----------------------------------------------------------------------
compute sigma tensor
needed whenever p_target or h0_inv changes
-----------------------------------------------------------------------*/
void FixBocs::compute_sigma()
{
// if nreset_h0 > 0, reset vol0 and h0_inv
// every nreset_h0 timesteps
if (nreset_h0 > 0) {
int delta = update->ntimestep - update->beginstep;
if (delta % nreset_h0 == 0) {
if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd;
else vol0 = domain->xprd * domain->yprd;
h0_inv[0] = domain->h_inv[0];
h0_inv[1] = domain->h_inv[1];
h0_inv[2] = domain->h_inv[2];
h0_inv[3] = domain->h_inv[3];
h0_inv[4] = domain->h_inv[4];
h0_inv[5] = domain->h_inv[5];
}
}
// generate upper-triangular half of
// sigma = vol0*h0inv*(p_target-p_hydro)*h0inv^t
// units of sigma are are PV/L^2 e.g. atm.A
//
// [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ]
// [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ]
// [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ]
sigma[0] =
vol0*(h0_inv[0]*((p_target[0]-p_hydro)*h0_inv[0] +
p_target[5]*h0_inv[5]+p_target[4]*h0_inv[4]) +
h0_inv[5]*(p_target[5]*h0_inv[0] +
(p_target[1]-p_hydro)*h0_inv[5]+p_target[3]*h0_inv[4]) +
h0_inv[4]*(p_target[4]*h0_inv[0]+p_target[3]*h0_inv[5] +
(p_target[2]-p_hydro)*h0_inv[4]));
sigma[1] =
vol0*(h0_inv[1]*((p_target[1]-p_hydro)*h0_inv[1] +
p_target[3]*h0_inv[3]) +
h0_inv[3]*(p_target[3]*h0_inv[1] +
(p_target[2]-p_hydro)*h0_inv[3]));
sigma[2] =
vol0*(h0_inv[2]*((p_target[2]-p_hydro)*h0_inv[2]));
sigma[3] =
vol0*(h0_inv[1]*(p_target[3]*h0_inv[2]) +
h0_inv[3]*((p_target[2]-p_hydro)*h0_inv[2]));
sigma[4] =
vol0*(h0_inv[0]*(p_target[4]*h0_inv[2]) +
h0_inv[5]*(p_target[3]*h0_inv[2]) +
h0_inv[4]*((p_target[2]-p_hydro)*h0_inv[2]));
sigma[5] =
vol0*(h0_inv[0]*(p_target[5]*h0_inv[1]+p_target[4]*h0_inv[3]) +
h0_inv[5]*((p_target[1]-p_hydro)*h0_inv[1]+p_target[3]*h0_inv[3]) +
h0_inv[4]*(p_target[3]*h0_inv[1]+(p_target[2]-p_hydro)*h0_inv[3]));
}
/* ----------------------------------------------------------------------
compute strain energy
-----------------------------------------------------------------------*/
double FixBocs::compute_strain_energy()
{
// compute strain energy = 0.5*Tr(sigma*h*h^t) in energy units
double* h = domain->h;
double d0,d1,d2;
d0 =
sigma[0]*(h[0]*h[0]+h[5]*h[5]+h[4]*h[4]) +
sigma[5]*( h[1]*h[5]+h[3]*h[4]) +
sigma[4]*( h[2]*h[4]);
d1 =
sigma[5]*( h[5]*h[1]+h[4]*h[3]) +
sigma[1]*( h[1]*h[1]+h[3]*h[3]) +
sigma[3]*( h[2]*h[3]);
d2 =
sigma[4]*( h[4]*h[2]) +
sigma[3]*( h[3]*h[2]) +
sigma[2]*( h[2]*h[2]);
double energy = 0.5*(d0+d1+d2)/nktv2p;
return energy;
}
/* ----------------------------------------------------------------------
compute deviatoric barostat force = h*sigma*h^t
-----------------------------------------------------------------------*/
void FixBocs::compute_deviatoric()
{
// generate upper-triangular part of h*sigma*h^t
// units of fdev are are PV, e.g. atm*A^3
// [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ]
// [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ]
// [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ]
double* h = domain->h;
fdev[0] =
h[0]*(sigma[0]*h[0]+sigma[5]*h[5]+sigma[4]*h[4]) +
h[5]*(sigma[5]*h[0]+sigma[1]*h[5]+sigma[3]*h[4]) +
h[4]*(sigma[4]*h[0]+sigma[3]*h[5]+sigma[2]*h[4]);
fdev[1] =
h[1]*( sigma[1]*h[1]+sigma[3]*h[3]) +
h[3]*( sigma[3]*h[1]+sigma[2]*h[3]);
fdev[2] =
h[2]*( sigma[2]*h[2]);
fdev[3] =
h[1]*( sigma[3]*h[2]) +
h[3]*( sigma[2]*h[2]);
fdev[4] =
h[0]*( sigma[4]*h[2]) +
h[5]*( sigma[3]*h[2]) +
h[4]*( sigma[2]*h[2]);
fdev[5] =
h[0]*( sigma[5]*h[1]+sigma[4]*h[3]) +
h[5]*( sigma[1]*h[1]+sigma[3]*h[3]) +
h[4]*( sigma[3]*h[1]+sigma[2]*h[3]);
}
/* ----------------------------------------------------------------------
compute target temperature and kinetic energy
-----------------------------------------------------------------------*/
void FixBocs::compute_temp_target()
{
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
t_target = t_start + delta * (t_stop-t_start);
ke_target = tdof * boltz * t_target;
}
/* ----------------------------------------------------------------------
compute hydrostatic target pressure
-----------------------------------------------------------------------*/
void FixBocs::compute_press_target()
{
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
p_hydro = 0.0;
for (int i = 0; i < 3; i++)
if (p_flag[i]) {
p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
p_hydro += p_target[i];
}
if (pdim > 0) p_hydro /= pdim;
if (pstyle == TRICLINIC)
for (int i = 3; i < 6; i++)
p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
// if deviatoric, recompute sigma each time p_target changes
if (deviatoric_flag) compute_sigma();
}
/* ----------------------------------------------------------------------
update omega_dot, omega
-----------------------------------------------------------------------*/
void FixBocs::nh_omega_dot()
{
double f_omega,volume;
if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
else volume = domain->xprd*domain->yprd;
if (deviatoric_flag) compute_deviatoric();
mtk_term1 = 0.0;
if (mtk_flag) {
if (pstyle == ISO) {
mtk_term1 = tdof * boltz * t_current;
mtk_term1 /= pdim * atom->natoms;
} else {
double *mvv_current = temperature->vector;
for (int i = 0; i < 3; i++)
if (p_flag[i])
mtk_term1 += mvv_current[i];
mtk_term1 /= pdim * atom->natoms;
}
}
for (int i = 0; i < 3; i++)
if (p_flag[i]) {
f_omega = (p_current[i]-p_hydro)*volume /
(omega_mass[i] * nktv2p) + mtk_term1 / omega_mass[i];
if (deviatoric_flag) f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= pdrag_factor;
}
mtk_term2 = 0.0;
if (mtk_flag) {
for (int i = 0; i < 3; i++)
if (p_flag[i])
mtk_term2 += omega_dot[i];
if (pdim > 0) mtk_term2 /= pdim * atom->natoms;
}
if (pstyle == TRICLINIC) {
for (int i = 3; i < 6; i++) {
if (p_flag[i]) {
f_omega = p_current[i]*volume/(omega_mass[i] * nktv2p);
if (deviatoric_flag)
f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= pdrag_factor;
}
}
}
}
/* ----------------------------------------------------------------------
if any tilt ratios exceed limits, set flip = 1 and compute new tilt values
do not flip in x or y if non-periodic (can tilt but not flip)
this is b/c the box length would be changed (dramatically) by flip
if yz tilt exceeded, adjust C vector by one B vector
if xz tilt exceeded, adjust C vector by one A vector
if xy tilt exceeded, adjust B vector by one A vector
check yz first since it may change xz, then xz check comes after
if any flip occurs, create new box in domain
image_flip() adjusts image flags due to box shape change induced by flip
remap() puts atoms outside the new box back into the new box
perform irregular on atoms in lamda coords to migrate atoms to new procs
important that image_flip comes before remap, since remap may change
image flags to new values, making eqs in doc of Domain:image_flip incorrect
------------------------------------------------------------------------- */
void FixBocs::pre_exchange()
{
double xprd = domain->xprd;
double yprd = domain->yprd;
// flip is only triggered when tilt exceeds 0.5 by DELTAFLIP
// this avoids immediate re-flipping due to tilt oscillations
double xtiltmax = (0.5+DELTAFLIP)*xprd;
double ytiltmax = (0.5+DELTAFLIP)*yprd;
int flipxy,flipxz,flipyz;
flipxy = flipxz = flipyz = 0;
if (domain->yperiodic) {
if (domain->yz < -ytiltmax) {
domain->yz += yprd;
domain->xz += domain->xy;
flipyz = 1;
} else if (domain->yz >= ytiltmax) {
domain->yz -= yprd;
domain->xz -= domain->xy;
flipyz = -1;
}
}
if (domain->xperiodic) {
if (domain->xz < -xtiltmax) {
domain->xz += xprd;
flipxz = 1;
} else if (domain->xz >= xtiltmax) {
domain->xz -= xprd;
flipxz = -1;
}
if (domain->xy < -xtiltmax) {
domain->xy += xprd;
flipxy = 1;
} else if (domain->xy >= xtiltmax) {
domain->xy -= xprd;
flipxy = -1;
}
}
int flip = 0;
if (flipxy || flipxz || flipyz) flip = 1;
if (flip) {
domain->set_global_box();
domain->set_local_box();
domain->image_flip(flipxy,flipxz,flipyz);
double **x = atom->x;
imageint *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
domain->x2lamda(atom->nlocal);
irregular->migrate_atoms();
domain->lamda2x(atom->nlocal);
}
}
/* ----------------------------------------------------------------------
memory usage of Irregular
------------------------------------------------------------------------- */
double FixBocs::memory_usage()
{
double bytes = 0.0;
if (irregular) bytes += irregular->memory_usage();
return bytes;
}
Reference in New Issue View Git Blame Copy Permalink