Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F76659044
fix_npt.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Aug 9, 16:21
Size
23 KB
Mime Type
text/x-c
Expires
Sun, Aug 11, 16:21 (2 d)
Engine
blob
Format
Raw Data
Handle
19748561
Attached To
rLAMMPS lammps
fix_npt.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: Mark Stevens (SNL)
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "math.h"
#include "fix_npt.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "modify.h"
#include "fix_deform.h"
#include "compute.h"
#include "kspace.h"
#include "update.h"
#include "respa.h"
#include "domain.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
/* ---------------------------------------------------------------------- */
FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) error->all("Illegal fix npt command");
restart_global = 1;
pressure_every = 1;
box_change = 1;
scalar_flag = 1;
scalar_vector_freq = 1;
extensive = 1;
t_start = atof(arg[3]);
t_stop = atof(arg[4]);
double t_period = atof(arg[5]);
if (t_start < 0.0 || t_stop <= 0.0)
error->all("Target T for fix npt cannot be 0.0");
double p_period[3];
if (strcmp(arg[6],"xyz") == 0) {
if (narg < 10) error->all("Illegal fix npt command");
press_couple = 0;
p_start[0] = p_start[1] = p_start[2] = atof(arg[7]);
p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[8]);
p_period[0] = p_period[1] = p_period[2] = atof(arg[9]);
p_flag[0] = p_flag[1] = p_flag[2] = 1;
if (domain->dimension == 2) {
p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0;
}
} else {
if (strcmp(arg[6],"xy") == 0) press_couple = 1;
else if (strcmp(arg[6],"yz") == 0) press_couple = 2;
else if (strcmp(arg[6],"xz") == 0) press_couple = 3;
else if (strcmp(arg[6],"aniso") == 0) press_couple = 4;
else error->all("Illegal fix npt command");
if (narg < 14) error->all("Illegal fix npt command");
if (domain->dimension == 2 &&
(press_couple == 1 || press_couple == 2 || press_couple == 3))
error->all("Invalid fix npt command for a 2d simulation");
if (strcmp(arg[7],"NULL") == 0) {
p_start[0] = p_stop[0] = p_period[0] = 0.0;
p_flag[0] = 0;
} else {
p_start[0] = atof(arg[7]);
p_stop[0] = atof(arg[8]);
p_flag[0] = 1;
}
if (strcmp(arg[9],"NULL") == 0) {
p_start[1] = p_stop[1] = p_period[1] = 0.0;
p_flag[1] = 0;
} else {
p_start[1] = atof(arg[9]);
p_stop[1] = atof(arg[10]);
p_flag[1] = 1;
}
if (strcmp(arg[11],"NULL") == 0) {
p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0;
} else {
if (domain->dimension == 2)
error->all("Invalid fix npt command for a 2d simulation");
p_start[2] = atof(arg[11]);
p_stop[2] = atof(arg[12]);
p_flag[2] = 1;
}
double period = atof(arg[13]);
if (p_flag[0]) p_period[0] = period;
if (p_flag[1]) p_period[1] = period;
if (p_flag[2]) p_period[2] = period;
}
// process extra keywords
drag = 0.0;
allremap = 1;
int iarg;
if (press_couple == 0) iarg = 10;
else iarg = 14;
while (iarg < narg) {
if (strcmp(arg[iarg],"drag") == 0) {
if (iarg+2 > narg) error->all("Illegal fix npt command");
drag = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"dilate") == 0) {
if (iarg+2 > narg) error->all("Illegal fix npt command");
if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0;
else error->all("Illegal fix npt command");
iarg += 2;
} else error->all("Illegal fix npt command");
}
// error checks
if (p_flag[0] && domain->xperiodic == 0)
error->all("Cannot use fix npt on a non-periodic dimension");
if (p_flag[1] && domain->yperiodic == 0)
error->all("Cannot use fix npt on a non-periodic dimension");
if (p_flag[2] && domain->zperiodic == 0)
error->all("Cannot use fix npt on a non-periodic dimension");
// convert input periods to frequencies
if (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("Fix npt periods must be > 0.0");
t_freq = 1.0 / t_period;
p_freq[0] = p_freq[1] = p_freq[2] = 0.0;
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];
// create a new compute temp style
// id = fix-ID + temp
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
if (strcmp(style,"npt") == 0) newarg[2] = (char *) "temp";
else if (strcmp(style,"npt/asphere") == 0)
newarg[2] = (char *) "temp/asphere";
modify->add_compute(3,newarg);
delete [] newarg;
tflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
pflag = 1;
// Nose/Hoover temp and pressure init
eta = eta_dot = 0.0;
omega[0] = omega[1] = omega[2] = 0.0;
omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0;
nrigid = 0;
rfix = NULL;
}
/* ---------------------------------------------------------------------- */
FixNPT::~FixNPT()
{
delete [] rfix;
// delete temperature and pressure if fix created them
if (tflag) modify->delete_compute(id_temp);
if (pflag) modify->delete_compute(id_press);
delete [] id_temp;
delete [] id_press;
}
/* ---------------------------------------------------------------------- */
int FixNPT::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
mask |= THERMO_ENERGY;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixNPT::init()
{
if (domain->triclinic) error->all("Cannot use fix npt with triclinic box");
if (atom->mass == NULL)
error->all("Cannot use fix npt without per-type mass defined");
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("Cannot use fix npt and fix deform on same dimension");
}
// set temperature and pressure ptrs
// set ptemperature only if pressure's id_pre[0] is not id_temp
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all("Temp ID for fix npt does not exist");
temperature = modify->compute[icompute];
icompute = modify->find_compute(id_press);
if (icompute < 0) error->all("Press ID for fix npt does not exist");
pressure = modify->compute[icompute];
if (strcmp(id_temp,pressure->id_pre[0]) == 0) ptemperature = NULL;
else {
icompute = modify->find_compute(pressure->id_pre[0]);
if (icompute < 0)
error->all("Temp ID of press ID for fix npt does not exist");
ptemperature = modify->compute[icompute];
}
// set timesteps and frequencies
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dtq = 0.5 * update->dt;
dthalf = 0.5 * update->dt;
double freq = MAX(p_freq[0],p_freq[1]);
freq = MAX(freq,p_freq[2]);
drag_factor = 1.0 - (update->dt * freq * drag);
boltz = force->boltz;
nktv2p = force->nktv2p;
dimension = domain->dimension;
if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd;
else vol0 = domain->xprd * domain->yprd;
if (force->kspace) kspace_flag = 1;
else kspace_flag = 0;
if (strcmp(update->integrate_style,"respa") == 0) {
nlevels_respa = ((Respa *) update->integrate)->nlevels;
step_respa = ((Respa *) update->integrate)->step;
}
// detect if any rigid fixes exist so rigid bodies move when box is remapped
// rfix[] = indices to each fix rigid
delete [] rfix;
nrigid = 0;
rfix = NULL;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->rigid_flag) nrigid++;
if (nrigid) {
rfix = new int[nrigid];
nrigid = 0;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i;
}
}
/* ----------------------------------------------------------------------
compute T,P before integrator starts
------------------------------------------------------------------------- */
void FixNPT::setup()
{
t_target = t_start; // used by compute_scalar()
p_target[0] = p_start[0];
p_target[1] = p_start[1];
p_target[2] = p_start[2];
t_current = temperature->compute_scalar();
if (press_couple == 0) {
if (ptemperature) double ptmp = ptemperature->compute_scalar();
double tmp = pressure->compute_scalar();
} else {
temperature->compute_vector();
if (ptemperature) ptemperature->compute_vector();
pressure->compute_vector();
}
couple();
}
/* ----------------------------------------------------------------------
1st half of Verlet update
------------------------------------------------------------------------- */
void FixNPT::initial_integrate()
{
int i;
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
// update eta_dot
t_target = t_start + delta * (t_stop-t_start);
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dthalf;
eta_dot *= drag_factor;
eta += dtv*eta_dot;
// update omega_dot
// for non-varying dims, p_freq is 0.0, so omega_dot doesn't change
double f_omega,volume;
if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
else volume = domain->xprd*domain->yprd;
double denskt = atom->natoms*boltz*t_target / volume * nktv2p;
for (i = 0; i < 3; i++) {
p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= drag_factor;
omega[i] += dtv*omega_dot[i];
factor[i] = exp(-dthalf*(eta_dot+omega_dot[i]));
dilation[i] = exp(dthalf*omega_dot[i]);
}
// v update only for atoms in NPT group
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double dtfm;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
}
}
// remap simulation box and all owned atoms by 1/2 step
remap(0);
// x update by full step only for atoms in NPT group
for (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];
}
}
// remap simulation box and all owned atoms by 1/2 step
// redo KSpace coeffs since volume has changed
remap(0);
if (kspace_flag) force->kspace->setup();
}
/* ----------------------------------------------------------------------
2nd half of Verlet update
------------------------------------------------------------------------- */
void FixNPT::final_integrate()
{
int i;
// v update only for atoms in NPT group
double **v = atom->v;
double **f = atom->f;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double dtfm;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
}
}
// compute new T,P
t_current = temperature->compute_scalar();
if (press_couple == 0) {
if (ptemperature) double ptmp = ptemperature->compute_scalar();
double tmp = pressure->compute_scalar();
} else {
temperature->compute_vector();
if (ptemperature) ptemperature->compute_vector();
pressure->compute_vector();
}
couple();
// update eta_dot
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dthalf;
eta_dot *= drag_factor;
// update omega_dot
// for non-varying dims, p_freq is 0.0, so omega_dot doesn't change
double f_omega,volume;
if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
else volume = domain->xprd*domain->yprd;
double denskt = atom->natoms*boltz*t_target / volume * nktv2p;
for (i = 0; i < 3; i++) {
f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= drag_factor;
}
}
/* ---------------------------------------------------------------------- */
void FixNPT::initial_integrate_respa(int ilevel, int flag)
{
// if flag = 1, then is 2nd call at outermost level from rRESPA
// perform 2nd half of box remap on own + ghost atoms and return
// redo KSpace coeffs since volume has changed
if (flag == 1) {
remap(1);
if (kspace_flag) force->kspace->setup();
return;
}
// set timesteps by level
double dtfm;
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dthalf = 0.5 * step_respa[ilevel];
// atom quantities
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// outermost level - update eta_dot and omega_dot, apply to v, remap box
// all other levels - NVE update of v
// x,v updates only performed for atoms in NPT group
if (ilevel == nlevels_respa-1) {
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
// update eta_dot
t_target = t_start + delta * (t_stop-t_start);
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dthalf;
eta_dot *= drag_factor;
eta += dtv*eta_dot;
// update omega_dot
// for non-varying dims, p_freq is 0.0, so omega_dot doesn't change
double f_omega,volume;
if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
else volume = domain->xprd*domain->yprd;
double denskt = atom->natoms*boltz*t_target / volume * nktv2p;
for (int i = 0; i < 3; i++) {
p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= drag_factor;
omega[i] += dtv*omega_dot[i];
factor[i] = exp(-dthalf*(eta_dot+omega_dot[i]));
dilation[i] = exp(dthalf*omega_dot[i]);
}
// v update only for atoms in NPT group
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
}
}
// remap simulation box and all owned atoms by 1/2 step
remap(0);
} else {
// v update only for atoms in NPT group
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];
}
}
}
// innermost level - also update x only for atoms in NPT group
if (ilevel == 0) {
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];
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixNPT::final_integrate_respa(int ilevel)
{
double dtfm;
// 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 to v via final_integrate()
// all other levels - NVE update of v
// v update only performed for atoms in NPT group
if (ilevel == nlevels_respa-1) final_integrate();
else {
// v update only for atoms in NPT group
double **v = atom->v;
double **f = atom->f;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
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];
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixNPT::couple()
{
double *tensor = pressure->vector;
if (press_couple == 0)
p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
else if (press_couple == 1) {
double ave = 0.5 * (tensor[0] + tensor[1]);
p_current[0] = p_current[1] = ave;
p_current[2] = tensor[2];
} else if (press_couple == 2) {
double ave = 0.5 * (tensor[1] + tensor[2]);
p_current[1] = p_current[2] = ave;
p_current[0] = tensor[0];
} else if (press_couple == 3) {
double ave = 0.5 * (tensor[0] + tensor[2]);
p_current[0] = p_current[2] = ave;
p_current[1] = tensor[1];
} if (press_couple == 4) {
p_current[0] = tensor[0];
p_current[1] = tensor[1];
p_current[2] = tensor[2];
}
}
/* ----------------------------------------------------------------------
change box size
remap owned or owned+ghost atoms depending on flag
remap all atoms or fix group atoms depending on allremap flag
if rigid bodies exist, scale rigid body centers-of-mass
------------------------------------------------------------------------- */
void FixNPT::remap(int flag)
{
int i,n;
double oldlo,oldhi,ctr;
double **x = atom->x;
int *mask = atom->mask;
if (flag) n = atom->nlocal + atom->nghost;
else n = atom->nlocal;
// convert pertinent atoms and rigid bodies to lamda coords
if (allremap) domain->x2lamda(n);
else {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
domain->x2lamda(x[i],x[i]);
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(0);
// reset global and local box to new size/shape
for (i = 0; i < 3; i++) {
if (p_flag[i]) {
oldlo = domain->boxlo[i];
oldhi = domain->boxhi[i];
ctr = 0.5 * (oldlo + oldhi);
domain->boxlo[i] = (oldlo-ctr)*dilation[i] + ctr;
domain->boxhi[i] = (oldhi-ctr)*dilation[i] + ctr;
}
}
domain->set_global_box();
domain->set_local_box();
// convert pertinent atoms and rigid bodies back to box coords
if (allremap) domain->lamda2x(n);
else {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
domain->lamda2x(x[i],x[i]);
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(1);
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixNPT::write_restart(FILE *fp)
{
int n = 0;
double list[8];
list[n++] = eta;
list[n++] = eta_dot;
list[n++] = omega[0];
list[n++] = omega[1];
list[n++] = omega[2];
list[n++] = omega_dot[0];
list[n++] = omega_dot[1];
list[n++] = omega_dot[2];
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(&list,sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixNPT::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
eta = list[n++];
eta_dot = list[n++];
omega[0] = list[n++];
omega[1] = list[n++];
omega[2] = list[n++];
omega_dot[0] = list[n++];
omega_dot[1] = list[n++];
omega_dot[2] = list[n++];
}
/* ---------------------------------------------------------------------- */
int FixNPT::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all("Illegal fix_modify command");
if (tflag) {
modify->delete_compute(id_temp);
tflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all("Could not find fix_modify temp ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all("Fix_modify temp ID does not compute temperature");
if (temperature->igroup != 0 && comm->me == 0)
error->warning("Temperature for NPT is not for group all");
// reset id_pre[0] of pressure to new temp ID
icompute = modify->find_compute(id_press);
if (icompute < 0) error->all("Press ID for fix npt does not exist");
delete [] modify->compute[icompute]->id_pre[0];
modify->compute[icompute]->id_pre[0] = new char[n];
strcpy(modify->compute[icompute]->id_pre[0],id_temp);
return 2;
} else if (strcmp(arg[0],"press") == 0) {
if (narg < 2) error->all("Illegal fix_modify command");
if (pflag) {
modify->delete_compute(id_press);
pflag = 0;
}
delete [] id_press;
int n = strlen(arg[1]) + 1;
id_press = new char[n];
strcpy(id_press,arg[1]);
int icompute = modify->find_compute(id_press);
if (icompute < 0) error->all("Could not find fix_modify press ID");
pressure = modify->compute[icompute];
if (pressure->pressflag == 0)
error->all("Fix_modify press ID does not compute pressure");
return 2;
}
return 0;
}
/* ---------------------------------------------------------------------- */
double FixNPT::compute_scalar()
{
double ke = temperature->dof * boltz * t_target;
double keplus = atom->natoms * boltz * t_target;
double volume;
if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
else volume = domain->xprd * domain->yprd;
int pdim = p_flag[0] + p_flag[1] + p_flag[2];
double energy = ke * (eta + 0.5*eta_dot*eta_dot/(t_freq*t_freq));
for (int i = 0; i < 3; i++)
if (p_freq[i] > 0.0)
energy += 0.5*keplus*omega_dot[i]*omega_dot[i] /
(p_freq[i]*p_freq[i]) + p_target[i]*(volume-vol0) / (pdim*nktv2p);
return energy;
}
/* ---------------------------------------------------------------------- */
void FixNPT::reset_dt()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dtq = 0.5 * update->dt;
dthalf = 0.5 * update->dt;
double freq = MAX(p_freq[0],p_freq[1]);
freq = MAX(freq,p_freq[2]);
drag_factor = 1.0 - (update->dt * freq * drag);
}
Event Timeline
Log In to Comment