Homec4science
Diffusion lammps (lammm-master-deprecated) Browse

lammps/src/RIGID/fix_rigid_nh_small.cpp
e2790b1c630blammm-master-deprecated

Last Change
Raw File

fix_rigid_nh_small.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: Trung Dac Nguyen (ORNL)
references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005)
Miller et al., J Chem Phys. 116, 8649-8659 (2002)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "fix_rigid_nh_small.h"
#include "math_extra.h"
#include "atom.h"
#include "compute.h"
#include "domain.h"
#include "update.h"
#include "modify.h"
#include "fix_deform.h"
#include "group.h"
#include "comm.h"
#include "force.h"
#include "kspace.h"
#include "output.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathExtra;
enum{NONE,XYZ,XY,YZ,XZ}; // same as in FixRigid
enum{ISO,ANISO,TRICLINIC}; // same as in FixRigid
#define EPSILON 1.0e-7
enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF};
/* ---------------------------------------------------------------------- */
FixRigidNHSmall::FixRigidNHSmall(LAMMPS *lmp, int narg, char **arg) :
FixRigidSmall(lmp, narg, arg)
{
// error checks
if ((p_flag[0] == 1 && p_period[0] <= 0.0) ||
(p_flag[1] == 1 && p_period[1] <= 0.0) ||
(p_flag[2] == 1 && p_period[2] <= 0.0))
error->all(FLERR,"Fix rigid/small npt/nph period must be > 0.0");
dimension = domain->dimension;
if (dimension == 2 && p_flag[2])
error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"for a 2d simulation");
if (dimension == 2 && (pcouple == YZ || pcouple == XZ))
error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"for a 2d simulation");
if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0))
error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"pressure settings");
if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0)
error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"pressure settings");
if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0))
error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"pressure settings");
if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0))
error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"pressure settings");
if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0))
error->all(FLERR,"Invalid fix rigid/small npt/nph command "
"pressure settings");
// require periodicity in tensile dimension
if (p_flag[0] && domain->xperiodic == 0)
error->all(FLERR,
"Cannot use fix rigid/small npt/nph on a "
"non-periodic dimension");
if (p_flag[1] && domain->yperiodic == 0)
error->all(FLERR,
"Cannot use fix rigid/small npt/nph on a "
"non-periodic dimension");
if (p_flag[2] && domain->zperiodic == 0)
error->all(FLERR,
"Cannot use fix rigid/small npt/nph on a "
"non-periodic dimension");
if (pcouple == XYZ && dimension == 3 &&
(p_start[0] != p_start[1] || p_start[0] != p_start[2] ||
p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] ||
p_period[0] != p_period[1] || p_period[0] != p_period[2]))
error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings");
if (pcouple == XYZ && dimension == 2 &&
(p_start[0] != p_start[1] || p_stop[0] != p_stop[1] ||
p_period[0] != p_period[1]))
error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings");
if (pcouple == XY &&
(p_start[0] != p_start[1] || p_stop[0] != p_stop[1] ||
p_period[0] != p_period[1]))
error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings");
if (pcouple == YZ &&
(p_start[1] != p_start[2] || p_stop[1] != p_stop[2] ||
p_period[1] != p_period[2]))
error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings");
if (pcouple == XZ &&
(p_start[0] != p_start[2] || p_stop[0] != p_stop[2] ||
p_period[0] != p_period[2]))
error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings");
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))
error->all(FLERR,"Fix rigid/small nvt/npt/nph damping parameters "
"must be > 0.0");
// memory allocation and initialization
if (tstat_flag || pstat_flag) {
allocate_chain();
allocate_order();
}
if (tstat_flag) {
eta_t[0] = eta_r[0] = 0.0;
eta_dot_t[0] = eta_dot_r[0] = 0.0;
f_eta_t[0] = f_eta_r[0] = 0.0;
for (int i = 1; i < t_chain; i++) {
eta_t[i] = eta_r[i] = 0.0;
eta_dot_t[i] = eta_dot_r[i] = 0.0;
}
}
if (pstat_flag) {
epsilon_dot[0] = epsilon_dot[1] = epsilon_dot[2] = 0.0;
eta_b[0] = eta_dot_b[0] = f_eta_b[0] = 0.0;
for (int i = 1; i < p_chain; i++)
eta_b[i] = eta_dot_b[i] = 0.0;
}
// rigid body pointers
nrigidfix = 0;
rfix = NULL;
vol0 = 0.0;
t0 = 1.0;
tcomputeflag = 0;
pcomputeflag = 0;
}
/* ---------------------------------------------------------------------- */
FixRigidNHSmall::~FixRigidNHSmall()
{
if (tstat_flag || pstat_flag) {
deallocate_chain();
deallocate_order();
}
if (rfix) delete [] rfix;
if (tcomputeflag) {
modify->delete_compute(id_temp);
delete [] id_temp;
}
// delete pressure if fix created it
if (pstat_flag) {
if (pcomputeflag) modify->delete_compute(id_press);
delete [] id_press;
}
}
/* ---------------------------------------------------------------------- */
int FixRigidNHSmall::setmask()
{
int mask = 0;
mask = FixRigidSmall::setmask();
if (tstat_flag || pstat_flag) mask |= THERMO_ENERGY;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::init()
{
FixRigidSmall::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 rigid npt/nph dilate group ID does not exist");
dilate_group_bit = group->bitmask[idilate];
}
// initialize thermostats
// set timesteps, constants
// store Yoshida-Suzuki integrator parameters
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dtq = 0.5 * update->dt;
boltz = force->boltz;
nktv2p = force->nktv2p;
mvv2e = force->mvv2e;
if (force->kspace) kspace_flag = 1;
else kspace_flag = 0;
// see Table 1 in Kamberaj et al
if (tstat_flag || pstat_flag) {
if (t_order == 3) {
w[0] = 1.0 / (2.0 - pow(2.0, 1.0/3.0));
w[1] = 1.0 - 2.0*w[0];
w[2] = w[0];
} else if (t_order == 5) {
w[0] = 1.0 / (4.0 - pow(4.0, 1.0/3.0));
w[1] = w[0];
w[2] = 1.0 - 4.0 * w[0];
w[3] = w[0];
w[4] = w[0];
}
}
int icompute;
if (tcomputeflag) {
icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Temperature ID for fix rigid nvt/npt/nph does not exist");
temperature = modify->compute[icompute];
}
if (pstat_flag) {
if (domain->triclinic)
error->all(FLERR,"Fix rigid npt/nph does not yet allow triclinic box");
// ensure no conflict with fix deform
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"deform") == 0) {
int *dimflag = ((FixDeform *) modify->fix[i])->dimflag;
if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) ||
(p_flag[2] && dimflag[2]))
error->all(FLERR,"Cannot use fix rigid npt/nph and fix deform on "
"same component of stress tensor");
}
// set frequency
p_freq_max = 0.0;
p_freq_max = MAX(p_freq[0],p_freq[1]);
p_freq_max = MAX(p_freq_max,p_freq[2]);
// tally the number of dimensions that are barostatted
// set initial volume and reference cell, if not already done
pdim = p_flag[0] + p_flag[1] + p_flag[2];
if (vol0 == 0.0) {
if (dimension == 2) vol0 = domain->xprd * domain->yprd;
else vol0 = domain->xprd * domain->yprd * domain->zprd;
}
// set pressure compute ptr
icompute = modify->find_compute(id_press);
if (icompute < 0)
error->all(FLERR,"Pressure ID for fix rigid npt/nph does not exist");
pressure = modify->compute[icompute];
// detect if any rigid fixes exist so rigid bodies move on remap
// rfix[] = indices to each fix rigid
// this will include self
if (rfix) delete [] rfix;
nrigidfix = 0;
rfix = NULL;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->rigid_flag) nrigidfix++;
if (nrigidfix) {
rfix = new int[nrigidfix];
nrigidfix = 0;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->rigid_flag) rfix[nrigidfix++] = i;
}
}
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::setup(int vflag)
{
FixRigidSmall::setup(vflag);
// total translational and rotational degrees of freedom
nf_t = dimension * nlocal_body;
if (dimension == 3) {
nf_r = dimension * nlocal_body;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
for (int k = 0; k < dimension; k++)
if (fabs(b->inertia[k]) < EPSILON) nf_r--;
}
} else if (dimension == 2) {
nf_r = nlocal_body;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
if (fabs(b->inertia[2]) < EPSILON) nf_r--;
}
}
double nf[2], nfall[2];
nf[0] = nf_t;
nf[1] = nf_r;
MPI_Allreduce(nf,nfall,2,MPI_DOUBLE,MPI_SUM,world);
nf_t = nfall[0];
nf_r = nfall[1];
g_f = nf_t + nf_r;
onednft = 1.0 + (double)(dimension) / (double)g_f;
onednfr = (double) (dimension) / (double)g_f;
double mbody[3];
akin_t = akin_r = 0.0;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space,
b->angmom,mbody);
MathExtra::quatvec(b->quat,mbody,b->conjqm);
b->conjqm[0] *= 2.0;
b->conjqm[1] *= 2.0;
b->conjqm[2] *= 2.0;
b->conjqm[3] *= 2.0;
if (tstat_flag || pstat_flag) {
akin_t += b->mass*(b->vcm[0]*b->vcm[0] + b->vcm[1]*b->vcm[1] +
b->vcm[2]*b->vcm[2]);
akin_r += b->angmom[0]*b->omega[0] + b->angmom[1]*b->omega[1] +
b->angmom[2]*b->omega[2];
}
}
// accumulate translational and rotational kinetic energies
if (tstat_flag || pstat_flag) {
double ke[2],keall[2];
ke[0] = akin_t;
ke[1] = akin_r;
MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world);
akin_t = keall[0];
akin_r = keall[1];
}
// compute target temperature
if (tstat_flag) compute_temp_target();
else if (pstat_flag) {
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;
}
// compute target pressure
// compute current pressure
// trigger virial computation on next timestep
if (pstat_flag) {
compute_press_target();
if (pstyle == ISO) {
temperature->compute_scalar();
pressure->compute_scalar();
} else {
temperature->compute_vector();
pressure->compute_vector();
}
couple();
pressure->addstep(update->ntimestep+1);
}
// initialize thermostat/barostat settings
double kt, t_mass, tb_mass;
kt = boltz * t_target;
if (tstat_flag) {
t_mass = kt / (t_freq*t_freq);
q_t[0] = nf_t * t_mass;
q_r[0] = nf_r * t_mass;
for (int i = 1; i < t_chain; i++)
q_t[i] = q_r[i] = t_mass;
for (int i = 1; i < t_chain; i++) {
f_eta_t[i] = (q_t[i-1] * eta_dot_t[i-1] * eta_dot_t[i-1] - kt)/q_t[i];
f_eta_r[i] = (q_r[i-1] * eta_dot_r[i-1] * eta_dot_r[i-1] - kt)/q_r[i];
}
}
// initial forces on barostat thermostat variables
if (pstat_flag) {
for (int i = 0; i < 3; i++)
if (p_flag[i]) {
epsilon_mass[i] = (g_f + dimension) * kt / (p_freq[i]*p_freq[i]);
epsilon[i] = log(vol0)/dimension;
}
tb_mass = kt / (p_freq_max * p_freq_max);
q_b[0] = dimension * dimension * tb_mass;
for (int i = 1; i < p_chain; i++) {
q_b[i] = tb_mass;
f_eta_b[i] = (q_b[i] * eta_dot_b[i-1] * eta_dot_b[i-1] - kt)/q_b[i];
}
}
// update order/timestep dependent coefficients
if (tstat_flag || pstat_flag) {
for (int i = 0; i < t_order; i++) {
wdti1[i] = w[i] * dtv / t_iter;
wdti2[i] = wdti1[i] / 2.0;
wdti4[i] = wdti1[i] / 4.0;
}
}
if (pstat_flag) {
compute_press_target();
nh_epsilon_dot();
}
}
/* ----------------------------------------------------------------------
perform preforce velocity Verlet integration
see Kamberaj paper for step references
------------------------------------------------------------------------- */
void FixRigidNHSmall::initial_integrate(int vflag)
{
double tmp,scale_r,scale_t[3],scale_v[3];
double dtfm,mbody[3],tbody[3],fquat[4];
double dtf2 = dtf * 2.0;
// compute scale variables
scale_t[0] = scale_t[1] = scale_t[2] = 1.0;
scale_v[0] = scale_v[1] = scale_v[2] = 1.0;
scale_r = 1.0;
if (tstat_flag) {
tmp = exp(-dtq * eta_dot_t[0]);
scale_t[0] = scale_t[1] = scale_t[2] = tmp;
tmp = exp(-dtq * eta_dot_r[0]);
scale_r = tmp;
}
if (pstat_flag) {
scale_t[0] *= exp(-dtq * (epsilon_dot[0] + mtk_term2));
scale_t[1] *= exp(-dtq * (epsilon_dot[1] + mtk_term2));
scale_t[2] *= exp(-dtq * (epsilon_dot[2] + mtk_term2));
scale_r *= exp(-dtq * (pdim * mtk_term2));
tmp = dtq * epsilon_dot[0];
scale_v[0] = dtv * exp(tmp) * maclaurin_series(tmp);
tmp = dtq * epsilon_dot[1];
scale_v[1] = dtv * exp(tmp) * maclaurin_series(tmp);
tmp = dtq * epsilon_dot[2];
scale_v[2] = dtv * exp(tmp) * maclaurin_series(tmp);
}
// update xcm, vcm, quat, conjqm and angmom
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
// step 1.1 - update vcm by 1/2 step
dtfm = dtf / b->mass;
b->vcm[0] += dtfm * b->fcm[0];
b->vcm[1] += dtfm * b->fcm[1];
b->vcm[2] += dtfm * b->fcm[2];
if (tstat_flag || pstat_flag) {
b->vcm[0] *= scale_t[0];
b->vcm[1] *= scale_t[1];
b->vcm[2] *= scale_t[2];
}
// step 1.2 - update xcm by full step
if (!pstat_flag) {
b->xcm[0] += dtv * b->vcm[0];
b->xcm[1] += dtv * b->vcm[1];
b->xcm[2] += dtv * b->vcm[2];
} else {
b->xcm[0] += scale_v[0] * b->vcm[0];
b->xcm[1] += scale_v[1] * b->vcm[1];
b->xcm[2] += scale_v[2] * b->vcm[2];
}
// step 1.3 - apply torque (body coords) to quaternion momentum
MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space,
b->torque,tbody);
MathExtra::quatvec(b->quat,tbody,fquat);
b->conjqm[0] += dtf2 * fquat[0];
b->conjqm[1] += dtf2 * fquat[1];
b->conjqm[2] += dtf2 * fquat[2];
b->conjqm[3] += dtf2 * fquat[3];
if (tstat_flag || pstat_flag) {
b->conjqm[0] *= scale_r;
b->conjqm[1] *= scale_r;
b->conjqm[2] *= scale_r;
b->conjqm[3] *= scale_r;
}
// step 1.4 to 1.13 - use no_squish rotate to update p and q
no_squish_rotate(3,b->conjqm,b->quat,b->inertia,dtq);
no_squish_rotate(2,b->conjqm,b->quat,b->inertia,dtq);
no_squish_rotate(1,b->conjqm,b->quat,b->inertia,dtv);
no_squish_rotate(2,b->conjqm,b->quat,b->inertia,dtq);
no_squish_rotate(3,b->conjqm,b->quat,b->inertia,dtq);
// update exyz_space
// transform p back to angmom
// update angular velocity
MathExtra::q_to_exyz(b->quat,b->ex_space,b->ey_space,
b->ez_space);
MathExtra::invquatvec(b->quat,b->conjqm,mbody);
MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,
mbody,b->angmom);
b->angmom[0] *= 0.5;
b->angmom[1] *= 0.5;
b->angmom[2] *= 0.5;
MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
b->ez_space,b->inertia,b->omega);
}
// forward communicate updated info of all bodies
commflag = INITIAL;
comm->forward_comm_fix(this,26);
// accumulate translational and rotational kinetic energies
if (tstat_flag || pstat_flag) {
akin_t = akin_r = 0.0;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
akin_t += b->mass*(b->vcm[0]*b->vcm[0] + b->vcm[1]*b->vcm[1] +
b->vcm[2]*b->vcm[2]);
akin_r += b->angmom[0]*b->omega[0] + b->angmom[1]*b->omega[1] +
b->angmom[2]*b->omega[2];
}
double ke[2],keall[2];
ke[0] = akin_t;
ke[1] = akin_r;
MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world);
akin_t = keall[0];
akin_r = keall[1];
}
// compute target temperature
// update thermostat chains using akin_t and akin_r
// refer to update_nhcp() in Kamberaj et al.
if (tstat_flag) {
compute_temp_target();
nhc_temp_integrate();
}
// update thermostat chains coupled with barostat
// refer to update_nhcb() in Kamberaj et al.
if (pstat_flag) {
nhc_press_integrate();
}
// virial setup before call to set_xv
if (vflag) v_setup(vflag);
else evflag = 0;
// remap simulation box by 1/2 step
if (pstat_flag) remap();
// set coords/orient and velocity/rotation of atoms in rigid bodies
// from quarternion and omega
set_xv();
// remap simulation box by full step
// redo KSpace coeffs since volume has changed
if (pstat_flag) {
remap();
if (kspace_flag) force->kspace->setup();
}
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::final_integrate()
{
int i,ibody;
double tmp,scale_t[3],scale_r;
double dtfm;
double mbody[3],tbody[3],fquat[4];
double dtf2 = dtf * 2.0;
// compute scale variables
scale_t[0] = scale_t[1] = scale_t[2] = 1.0;
scale_r = 1.0;
if (tstat_flag) {
tmp = exp(-1.0 * dtq * eta_dot_t[0]);
scale_t[0] = scale_t[1] = scale_t[2] = tmp;
scale_r = exp(-1.0 * dtq * eta_dot_r[0]);
}
if (pstat_flag) {
scale_t[0] *= exp(-dtq * (epsilon_dot[0] + mtk_term2));
scale_t[1] *= exp(-dtq * (epsilon_dot[1] + mtk_term2));
scale_t[2] *= exp(-dtq * (epsilon_dot[2] + mtk_term2));
scale_r *= exp(-dtq * (pdim * mtk_term2));
}
// sum over atoms to get force and torque on rigid body
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;
double dx,dy,dz;
double unwrap[3];
double *xcm,*fcm,*tcm;
for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
fcm = body[ibody].fcm;
fcm[0] = fcm[1] = fcm[2] = 0.0;
tcm = body[ibody].torque;
tcm[0] = tcm[1] = tcm[2] = 0.0;
}
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
fcm = b->fcm;
fcm[0] += f[i][0];
fcm[1] += f[i][1];
fcm[2] += f[i][2];
domain->unmap(x[i],xcmimage[i],unwrap);
xcm = b->xcm;
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
tcm = b->torque;
tcm[0] += dy*f[i][2] - dz*f[i][1];
tcm[1] += dz*f[i][0] - dx*f[i][2];
tcm[2] += dx*f[i][1] - dy*f[i][0];
}
// extended particles add their torque to torque of body
if (extended) {
double **torque = atom->torque;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
if (eflags[i] & TORQUE) {
tcm = body[atom2body[i]].torque;
tcm[0] += torque[i][0];
tcm[1] += torque[i][1];
tcm[2] += torque[i][2];
}
}
}
// reverse communicate fcm, torque of all bodies
commflag = FORCE_TORQUE;
comm->reverse_comm_fix(this,6);
// include Langevin thermostat forces and torques
if (langflag) {
for (int ibody = 0; ibody < nlocal_body; ibody++) {
fcm = body[ibody].fcm;
fcm[0] += langextra[ibody][0];
fcm[1] += langextra[ibody][1];
fcm[2] += langextra[ibody][2];
tcm = body[ibody].torque;
tcm[0] += langextra[ibody][3];
tcm[1] += langextra[ibody][4];
tcm[2] += langextra[ibody][5];
}
}
// update vcm and angmom
// include Langevin thermostat forces
// fflag,tflag = 0 for some dimensions in 2d
for (ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
// update vcm by 1/2 step
dtfm = dtf / b->mass;
if (tstat_flag || pstat_flag) {
b->vcm[0] *= scale_t[0];
b->vcm[1] *= scale_t[1];
b->vcm[2] *= scale_t[2];
}
b->vcm[0] += dtfm * b->fcm[0];
b->vcm[1] += dtfm * b->fcm[1];
b->vcm[2] += dtfm * b->fcm[2];
// update conjqm, then transform to angmom, set velocity again
// virial is already setup from initial_integrate
MathExtra::transpose_matvec(b->ex_space,b->ey_space,
b->ez_space,b->torque,tbody);
MathExtra::quatvec(b->quat,tbody,fquat);
if (tstat_flag || pstat_flag) {
b->conjqm[0] = scale_r * b->conjqm[0] + dtf2 * fquat[0];
b->conjqm[1] = scale_r * b->conjqm[1] + dtf2 * fquat[1];
b->conjqm[2] = scale_r * b->conjqm[2] + dtf2 * fquat[2];
b->conjqm[3] = scale_r * b->conjqm[3] + dtf2 * fquat[3];
} else {
b->conjqm[0] += dtf2 * fquat[0];
b->conjqm[1] += dtf2 * fquat[1];
b->conjqm[2] += dtf2 * fquat[2];
b->conjqm[3] += dtf2 * fquat[3];
}
MathExtra::invquatvec(b->quat,b->conjqm,mbody);
MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,mbody,b->angmom);
b->angmom[0] *= 0.5;
b->angmom[1] *= 0.5;
b->angmom[2] *= 0.5;
MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
b->ez_space,b->inertia,b->omega);
}
// forward communicate updated info of all bodies
commflag = FINAL;
comm->forward_comm_fix(this,10);
// accumulate translational and rotational kinetic energies
if (pstat_flag) {
akin_t = akin_r = 0.0;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
akin_t += b->mass*(b->vcm[0]*b->vcm[0] + b->vcm[1]*b->vcm[1] +
b->vcm[2]*b->vcm[2]);
akin_r += b->angmom[0]*b->omega[0] + b->angmom[1]*b->omega[1] +
b->angmom[2]*b->omega[2];
}
double ke[2],keall[2];
ke[0] = akin_t;
ke[1] = akin_r;
MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world);
akin_t = keall[0];
akin_r = keall[1];
}
// set velocity/rotation of atoms in rigid bodies
// virial is already setup from initial_integrate
set_v();
// compute current temperature
if (tcomputeflag) t_current = temperature->compute_scalar();
// compute current and target pressures
// update epsilon dot using akin_t and akin_r
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);
compute_press_target();
nh_epsilon_dot();
}
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::nhc_temp_integrate()
{
int i,j,k;
double kt,gfkt_t,gfkt_r,tmp,ms,s,s2;
kt = boltz * t_target;
gfkt_t = nf_t * kt;
gfkt_r = nf_r * kt;
// update thermostat masses
double t_mass = boltz * t_target / (t_freq * t_freq);
q_t[0] = nf_t * t_mass;
q_r[0] = nf_r * t_mass;
for (i = 1; i < t_chain; i++)
q_t[i] = q_r[i] = t_mass;
// update force of thermostats coupled to particles
f_eta_t[0] = (akin_t * mvv2e - gfkt_t) / q_t[0];
f_eta_r[0] = (akin_r * mvv2e - gfkt_r) / q_r[0];
// multiple timestep iteration
for (i = 0; i < t_iter; i++) {
for (j = 0; j < t_order; j++) {
// update thermostat velocities half step
eta_dot_t[t_chain-1] += wdti2[j] * f_eta_t[t_chain-1];
eta_dot_r[t_chain-1] += wdti2[j] * f_eta_r[t_chain-1];
for (k = 1; k < t_chain; k++) {
tmp = wdti4[j] * eta_dot_t[t_chain-k];
ms = maclaurin_series(tmp);
s = exp(-1.0 * tmp);
s2 = s * s;
eta_dot_t[t_chain-k-1] = eta_dot_t[t_chain-k-1] * s2 +
wdti2[j] * f_eta_t[t_chain-k-1] * s * ms;
tmp = wdti4[j] * eta_dot_r[t_chain-k];
ms = maclaurin_series(tmp);
s = exp(-1.0 * tmp);
s2 = s * s;
eta_dot_r[t_chain-k-1] = eta_dot_r[t_chain-k-1] * s2 +
wdti2[j] * f_eta_r[t_chain-k-1] * s * ms;
}
// update thermostat positions a full step
for (k = 0; k < t_chain; k++) {
eta_t[k] += wdti1[j] * eta_dot_t[k];
eta_r[k] += wdti1[j] * eta_dot_r[k];
}
// update thermostat forces
for (k = 1; k < t_chain; k++) {
f_eta_t[k] = q_t[k-1] * eta_dot_t[k-1] * eta_dot_t[k-1] - kt;
f_eta_t[k] /= q_t[k];
f_eta_r[k] = q_r[k-1] * eta_dot_r[k-1] * eta_dot_r[k-1] - kt;
f_eta_r[k] /= q_r[k];
}
// update thermostat velocities a full step
for (k = 0; k < t_chain-1; k++) {
tmp = wdti4[j] * eta_dot_t[k+1];
ms = maclaurin_series(tmp);
s = exp(-1.0 * tmp);
s2 = s * s;
eta_dot_t[k] = eta_dot_t[k] * s2 + wdti2[j] * f_eta_t[k] * s * ms;
tmp = q_t[k] * eta_dot_t[k] * eta_dot_t[k] - kt;
f_eta_t[k+1] = tmp / q_t[k+1];
tmp = wdti4[j] * eta_dot_r[k+1];
ms = maclaurin_series(tmp);
s = exp(-1.0 * tmp);
s2 = s * s;
eta_dot_r[k] = eta_dot_r[k] * s2 + wdti2[j] * f_eta_r[k] * s * ms;
tmp = q_r[k] * eta_dot_r[k] * eta_dot_r[k] - kt;
f_eta_r[k+1] = tmp / q_r[k+1];
}
eta_dot_t[t_chain-1] += wdti2[j] * f_eta_t[t_chain-1];
eta_dot_r[t_chain-1] += wdti2[j] * f_eta_r[t_chain-1];
}
}
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::nhc_press_integrate()
{
int i,j,k;
double tmp,s,s2,ms,kecurrent;
double kt = boltz * t_target;
double lkt_press = kt;
// update thermostat masses
double tb_mass = kt / (p_freq_max * p_freq_max);
q_b[0] = dimension * dimension * tb_mass;
for (int i = 1; i < p_chain; i++) {
q_b[i] = tb_mass;
f_eta_b[i] = q_b[i-1] * eta_dot_b[i-1] * eta_dot_b[i-1] - kt;
f_eta_b[i] /= q_b[i];
}
// update forces acting on thermostat
kecurrent = 0.0;
for (i = 0; i < 3; i++)
if (p_flag[i]) {
epsilon_mass[i] = (g_f + dimension) * kt / (p_freq[i] * p_freq[i]);
kecurrent += epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i];
}
kecurrent /= pdim;
f_eta_b[0] = (kecurrent - lkt_press) / q_b[0];
// multiple timestep iteration
for (i = 0; i < t_iter; i++) {
for (j = 0; j < t_order; j++) {
// update thermostat velocities a half step
eta_dot_b[p_chain-1] += wdti2[j] * f_eta_b[p_chain-1];
for (k = 1; k < p_chain; k++) {
tmp = wdti4[j] * eta_dot_b[p_chain-k];
ms = maclaurin_series(tmp);
s = exp(-0.5 * tmp);
s2 = s * s;
eta_dot_b[p_chain-k-1] = eta_dot_b[p_chain-k-1] * s2 +
wdti2[j] * f_eta_b[p_chain-k-1] * s * ms;
}
// update thermostat positions
for (k = 0; k < p_chain; k++)
eta_b[k] += wdti1[j] * eta_dot_b[k];
// update thermostat forces
for (k = 1; k < p_chain; k++) {
f_eta_b[k] = q_b[k-1] * eta_dot_b[k-1] * eta_dot_b[k-1] - kt;
f_eta_b[k] /= q_b[k];
}
// update thermostat velocites a full step
for (k = 0; k < p_chain-1; k++) {
tmp = wdti4[j] * eta_dot_b[k+1];
ms = maclaurin_series(tmp);
s = exp(-0.5 * tmp);
s2 = s * s;
eta_dot_b[k] = eta_dot_b[k] * s2 + wdti2[j] * f_eta_b[k] * s * ms;
tmp = q_b[k] * eta_dot_b[k] * eta_dot_b[k] - kt;
f_eta_b[k+1] = tmp / q_b[k+1];
}
eta_dot_b[p_chain-1] += wdti2[j] * f_eta_b[p_chain-1];
}
}
}
/* ----------------------------------------------------------------------
compute kinetic energy in the extended Hamiltonian
conserved quantity = sum of returned energy and potential energy
-----------------------------------------------------------------------*/
double FixRigidNHSmall::compute_scalar()
{
int i,k;
double kt = boltz * t_target;
double energy,ke_t,ke_q,tmp,Pkq[4];
double *vcm,*quat;
// compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114)
// translational and rotational kinetic energies
ke_t = 0.0;
ke_q = 0.0;
for (int i = 0; i < nlocal_body; i++) {
vcm = body[i].vcm;
quat = body[i].quat;
ke_t += body[i].mass * (vcm[0]*vcm[0] + vcm[1]*vcm[1] +
vcm[2]*vcm[2]);
for (k = 1; k < 4; k++) {
if (k == 1) {
Pkq[0] = -quat[1];
Pkq[1] = quat[0];
Pkq[2] = quat[3];
Pkq[3] = -quat[2];
} else if (k == 2) {
Pkq[0] = -quat[2];
Pkq[1] = -quat[3];
Pkq[2] = quat[0];
Pkq[3] = quat[1];
} else if (k == 3) {
Pkq[0] = -quat[3];
Pkq[1] = quat[2];
Pkq[2] = -quat[1];
Pkq[3] = quat[0];
}
tmp = body[i].conjqm[0]*Pkq[0] + body[i].conjqm[1]*Pkq[1] +
body[i].conjqm[2]*Pkq[2] + body[i].conjqm[3]*Pkq[3];
tmp *= tmp;
if (fabs(body[i].inertia[k-1]) < 1e-6) tmp = 0.0;
else tmp /= (8.0 * body[i].inertia[k-1]);
ke_q += tmp;
}
}
double ke[2],keall[2];
ke[0] = ke_t;
ke[1] = ke_q;
MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world);
ke_t = keall[0];
ke_q = keall[1];
energy = (ke_t + ke_q) * mvv2e;
if (tstat_flag) {
// thermostat chain energy: from equation 12 in Kameraj et al (JCP 2005)
energy += kt * (nf_t * eta_t[0] + nf_r * eta_r[0]);
for (i = 1; i < t_chain; i++)
energy += kt * (eta_t[i] + eta_r[i]);
for (i = 0; i < t_chain; i++) {
energy += 0.5 * q_t[i] * (eta_dot_t[i] * eta_dot_t[i]);
energy += 0.5 * q_r[i] * (eta_dot_r[i] * eta_dot_r[i]);
}
}
if (pstat_flag) {
// using equation 22 in Kameraj et al for H_NPT
double e = 0.0;
for (i = 0; i < 3; i++)
if (p_flag[i])
e += epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i];
energy += e*(0.5/pdim);
double vol;
if (dimension == 2) vol = domain->xprd * domain->yprd;
else vol = domain->xprd * domain->yprd * domain->zprd;
double p0 = (p_target[0] + p_target[1] + p_target[2]) / 3.0;
energy += p0 * vol / nktv2p;
for (i = 0; i < p_chain; i++) {
energy += kt * eta_b[i];
energy += 0.5 * q_b[i] * (eta_dot_b[i] * eta_dot_b[i]);
}
}
return energy;
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::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];
}
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::remap()
{
int i;
double oldlo,oldhi,ctr,expfac;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// epsilon is not used, except for book-keeping
for (i = 0; i < 3; i++) epsilon[i] += dtq * epsilon_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 (nrigidfix)
for (i = 0; i < nrigidfix; 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]) {
oldlo = domain->boxlo[i];
oldhi = domain->boxhi[i];
ctr = 0.5 * (oldlo + oldhi);
expfac = exp(dtq * epsilon_dot[i]);
domain->boxlo[i] = (oldlo-ctr)*expfac + ctr;
domain->boxhi[i] = (oldhi-ctr)*expfac + ctr;
}
}
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 (nrigidfix)
for (i = 0; i< nrigidfix; i++)
modify->fix[rfix[i]]->deform(1);
}
/* ----------------------------------------------------------------------
compute target temperature and kinetic energy
-----------------------------------------------------------------------*/
void FixRigidNHSmall::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);
}
/* ----------------------------------------------------------------------
compute hydrostatic target pressure
-----------------------------------------------------------------------*/
void FixRigidNHSmall::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];
}
p_hydro /= pdim;
}
/* ----------------------------------------------------------------------
update epsilon_dot
-----------------------------------------------------------------------*/
void FixRigidNHSmall::nh_epsilon_dot()
{
int i;
double volume,scale,f_epsilon;
if (dimension == 2) volume = domain->xprd*domain->yprd;
else volume = domain->xprd*domain->yprd*domain->zprd;
// MTK terms
mtk_term1 = (akin_t + akin_r) * mvv2e / g_f;
scale = exp(-1.0 * dtq * eta_dot_b[0]);
for (i = 0; i < 3; i++)
if (p_flag[i]) {
f_epsilon = (p_current[i]-p_hydro)*volume / nktv2p + mtk_term1;
f_epsilon /= epsilon_mass[i];
epsilon_dot[i] += dtq * f_epsilon;
epsilon_dot[i] *= scale;
}
mtk_term2 = 0.0;
for (i = 0; i < 3; i++)
if (p_flag[i]) mtk_term2 += epsilon_dot[i];
mtk_term2 /= g_f;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixRigidNHSmall::write_restart(FILE *fp)
{
if (tstat_flag == 0 && pstat_flag == 0) return;
int nsize = 2; // tstat_flag and pstat_flag
if (tstat_flag) {
nsize += 1; // t_chain
nsize += 4*t_chain; // eta_t, eta_r, eta_dot_t, eta_dot_r
}
if (pstat_flag) {
nsize += 7; // p_chain, epsilon(3) and epsilon_dot(3)
nsize += 2*p_chain;
}
double *list;
memory->create(list,nsize,"rigid_nh:list");
int n = 0;
list[n++] = tstat_flag;
if (tstat_flag) {
list[n++] = t_chain;
for (int i = 0; i < t_chain; i++) {
list[n++] = eta_t[i];
list[n++] = eta_r[i];
list[n++] = eta_dot_t[i];
list[n++] = eta_dot_r[i];
}
}
list[n++] = pstat_flag;
if (pstat_flag) {
list[n++] = epsilon[0];
list[n++] = epsilon[1];
list[n++] = epsilon[2];
list[n++] = epsilon_dot[0];
list[n++] = epsilon_dot[1];
list[n++] = epsilon_dot[2];
list[n++] = p_chain;
for (int i = 0; i < p_chain; i++) {
list[n++] = eta_b[i];
list[n++] = eta_dot_b[i];
}
}
if (comm->me == 0) {
int size = (nsize)*sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(list,sizeof(double),nsize,fp);
}
memory->destroy(list);
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixRigidNHSmall::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
int flag = static_cast<int> (list[n++]);
if (flag) {
int m = static_cast<int> (list[n++]);
if (tstat_flag && m == t_chain) {
for (int i = 0; i < t_chain; i++) {
eta_t[i] = list[n++];
eta_r[i] = list[n++];
eta_dot_t[i] = list[n++];
eta_dot_r[i] = list[n++];
}
} else n += 4*m;
}
flag = static_cast<int> (list[n++]);
if (flag) {
epsilon[0] = list[n++];
epsilon[1] = list[n++];
epsilon[2] = list[n++];
epsilon_dot[0] = list[n++];
epsilon_dot[1] = list[n++];
epsilon_dot[2] = list[n++];
int m = static_cast<int> (list[n++]);
if (pstat_flag && m == p_chain) {
for (int i = 0; i < p_chain; i++) {
eta_b[i] = list[n++];
eta_dot_b[i] = list[n++];
}
} else n += 2*m;
}
}
/* ---------------------------------------------------------------------- */
int FixRigidNHSmall::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command");
if (tcomputeflag) {
modify->delete_compute(id_temp);
tcomputeflag = 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(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;
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(FLERR,"Could not find fix_modify pressure ID");
pressure = modify->compute[icompute];
if (pressure->pressflag == 0)
error->all(FLERR,"Fix_modify pressure ID does not compute pressure");
return 2;
}
return 0;
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::allocate_chain()
{
if (tstat_flag) {
q_t = new double[t_chain];
q_r = new double[t_chain];
eta_t = new double[t_chain];
eta_r = new double[t_chain];
eta_dot_t = new double[t_chain];
eta_dot_r = new double[t_chain];
f_eta_t = new double[t_chain];
f_eta_r = new double[t_chain];
}
if (pstat_flag) {
q_b = new double[p_chain];
eta_b = new double[p_chain];
eta_dot_b = new double[p_chain];
f_eta_b = new double[p_chain];
}
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::reset_target(double t_new)
{
t_start = t_stop = t_new;
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::allocate_order()
{
w = new double[t_order];
wdti1 = new double[t_order];
wdti2 = new double[t_order];
wdti4 = new double[t_order];
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::deallocate_chain()
{
if (tstat_flag) {
delete [] q_t;
delete [] q_r;
delete [] eta_t;
delete [] eta_r;
delete [] eta_dot_t;
delete [] eta_dot_r;
delete [] f_eta_t;
delete [] f_eta_r;
}
if (pstat_flag) {
delete [] q_b;
delete [] eta_b;
delete [] eta_dot_b;
delete [] f_eta_b;
}
}
/* ---------------------------------------------------------------------- */
void FixRigidNHSmall::deallocate_order()
{
delete [] w;
delete [] wdti1;
delete [] wdti2;
delete [] wdti4;
}
c4science · Help