jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added triclinic cell keywords and tilt factor scaling options 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6432 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps 2011-06-16 16:30:04 +00:00
parent c4df1d931a
commit dfb085b321
5 changed files with 358 additions and 52 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/fix_deform.cpp
View File

@ -483,7 +483,7 @@ void FixDeform::init()
// test for WIGGLE is on min/max oscillation limit, not tilt_stop
// this is b/c the flips would induce continuous changes in xz
// in order to keep the edge vectors of the flipped shape matrix
// a linear combination of the edge vectors of the unflipped shape matrix
// an integer combination of the edge vectors of the unflipped shape matrix
if (set[3].style && set[5].style) {
int flag = 0;
@ -729,7 +729,7 @@ void FixDeform::end_of_step()
// flip will be performed on next timestep before reneighboring
// when yz flips and xy is non-zero, xz must also change
// this is to keep the edge vectors of the flipped shape matrix
// a linear combination of the edge vectors of the unflipped shape matrix
// an integer combination of the edge vectors of the unflipped shape matrix
if (triclinic) {
double xprd = set[0].hi_target - set[0].lo_target;

286
src/fix_nh.cpp
View File

@ -23,6 +23,7 @@
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "irregular.h"
#include "modify.h"
#include "fix_deform.h"
#include "compute.h"
@ -35,6 +36,9 @@
using namespace LAMMPS_NS;
#define DELTAFLIP 0.1
#define TILTMAX 1.5
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
@ -69,6 +73,17 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
deviatoric_flag = 0;
nreset_h0 = 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;
}
// Used by FixNVTSllod to preserve non-default value
mtchain_default_flag = 1;
@ -84,8 +99,6 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
// process keywords
dimension = domain->dimension;
int iarg = 3;
while (iarg < narg) {
@ -126,6 +139,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
} else if (strcmp(arg[iarg],"tri") == 0) {
if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command");
pcouple = NONE;
scalexy = scalexz = scaleyz = 0;
p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]);
p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]);
p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]);
@ -143,7 +157,6 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
p_flag[4] = 0;
}
iarg += 4;
} else if (strcmp(arg[iarg],"x") == 0) {
if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command");
p_start[0] = atof(arg[iarg+1]);
@ -173,6 +186,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
} else if (strcmp(arg[iarg],"yz") == 0) {
if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command");
scaleyz = 0;
p_start[3] = atof(arg[iarg+1]);
p_stop[3] = atof(arg[iarg+2]);
p_period[3] = atof(arg[iarg+3]);
@ -183,6 +197,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
error->all("Invalid fix nvt/npt/nph command for a 2d simulation");
} else if (strcmp(arg[iarg],"xz") == 0) {
if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command");
scalexz = 0;
p_start[4] = atof(arg[iarg+1]);
p_stop[4] = atof(arg[iarg+2]);
p_period[4] = atof(arg[iarg+3]);
@ -192,6 +207,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
if (dimension == 2)
error->all("Invalid fix nvt/npt/nph command for a 2d simulation");
} else if (strcmp(arg[iarg],"xy") == 0) {
scalexy = 0;
if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command");
p_start[5] = atof(arg[iarg+1]);
p_stop[5] = atof(arg[iarg+2]);
@ -224,7 +240,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
} else if (strcmp(arg[iarg],"tchain") == 0) {
if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command");
mtchain = atoi(arg[iarg+1]);
// Used by FixNVTSllod to preserve non-default value
// used by FixNVTSllod to preserve non-default value
mtchain_default_flag = 0;
if (mtchain < 1) error->all("Illegal fix nvt/npt/nph command");
iarg += 2;
@ -254,6 +270,24 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
nreset_h0 = atoi(arg[iarg+1]);
if (nreset_h0 < 0) error->all("Illegal fix nvt/npt/nph command");
iarg += 2;
} else if (strcmp(arg[iarg],"scalexy") == 0) {
if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command");
if (strcmp(arg[iarg+1],"yes") == 0) scalexy = 1;
else if (strcmp(arg[iarg+1],"no") == 0) scalexy = 0;
else error->all("Illegal fix nvt/npt/nph command");
iarg += 2;
} else if (strcmp(arg[iarg],"scalexz") == 0) {
if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command");
if (strcmp(arg[iarg+1],"yes") == 0) scalexz = 1;
else if (strcmp(arg[iarg+1],"no") == 0) scalexz = 0;
else error->all("Illegal fix nvt/npt/nph command");
iarg += 2;
} else if (strcmp(arg[iarg],"scaleyz") == 0) {
if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command");
if (strcmp(arg[iarg+1],"yes") == 0) scaleyz = 1;
else if (strcmp(arg[iarg+1],"no") == 0) scaleyz = 0;
else error->all("Illegal fix nvt/npt/nph command");
iarg += 2;
} else error->all("Illegal fix nvt/npt/nph command");
}
@ -263,6 +297,8 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
error->all("Invalid fix nvt/npt/nph command for a 2d simulation");
if (dimension == 2 && (pcouple == YZ || pcouple == XZ))
error->all("Invalid fix nvt/npt/nph command for a 2d simulation");
if (dimension == 2 && (scalexz == 1 || scaleyz == 1 ))
error->all("Invalid fix nvt/npt/nph command for a 2d simulation");
if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0))
error->all("Invalid fix nvt/npt/nph command pressure settings");
@ -288,6 +324,26 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
if (p_flag[5] && domain->yperiodic == 0)
error->all("Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension");
if (scaleyz == 1 && domain->zperiodic == 0)
error->all("Cannot use fix nvt/npt/nph "
"with yz dynamics when z is non-periodic dimension");
if (scalexz == 1 && domain->zperiodic == 0)
error->all("Cannot use fix nvt/npt/nph "
"with xz dynamics when z is non-periodic dimension");
if (scalexy == 1 && domain->yperiodic == 0)
error->all("Cannot use fix nvt/npt/nph "
"with xy dynamics when y is non-periodic dimension");
if (p_flag[3] && scaleyz == 1)
error->all("Cannot use fix nvt/npt/nph with"
"both yz dynamics and yz scaling");
if (p_flag[4] && scalexz == 1)
error->all("Cannot use fix nvt/npt/nph with "
"both xz dynamics and xz scaling");
if (p_flag[5] && scalexy == 1)
error->all("Cannot use fix nvt/npt/nph with "
"both xy dynamics and xy scaling");
if (!domain->triclinic && (p_flag[3] || p_flag[4] || p_flag[5]))
error->all("Can not specify Pxy/Pxz/Pyz in "
"fix nvt/npt/nph with non-triclinic box");
@ -345,6 +401,11 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
else if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO;
else pstyle = ANISO;
// reneighboring only forced if flips will occur due to shape changes
if (p_flag[3] || p_flag[4] || p_flag[5]) force_reneighbor = 1;
if (scaleyz || scalexz || scalexy) force_reneighbor = 1;
// convert input periods to frequencies
t_freq = 0.0;
@ -412,6 +473,9 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
nrigid = 0;
rfix = NULL;
if (force_reneighbor) irregular = new Irregular(lmp);
else irregular = NULL;
// initialize vol0,t0 to zero to signal uninitialized
// values then assigned in init(), if necessary
@ -424,6 +488,8 @@ FixNH::~FixNH()
{
delete [] rfix;
delete irregular;
// delete temperature and pressure if fix created them
if (tflag) modify->delete_compute(id_temp);
@ -458,6 +524,7 @@ int FixNH::setmask()
mask |= THERMO_ENERGY;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
if (force_reneighbor) mask |= PRE_EXCHANGE;
return mask;
}
@ -465,7 +532,7 @@ int FixNH::setmask()
void FixNH::init()
{
// insure no conflict with fix deform
// ensure no conflict with fix deform
if (pstat_flag)
for (int i = 0; i < modify->nfix; i++)
@ -519,7 +586,6 @@ void FixNH::init()
tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
// tally the number of dimensions that are barostatted
// also compute the initial volume and reference cell
// set initial volume and reference cell, if not already done
if (pstat_flag) {
@ -877,6 +943,9 @@ void FixNH::remap()
{
int i;
double oldlo,oldhi,ctr;
double cosalpha, cosbeta, cosgamma;
double ly, lz, singamma, yzc, lzc, clen, blen;
double expfac;
double **x = atom->x;
int *mask = atom->mask;
@ -902,7 +971,7 @@ void FixNH::remap()
// reset global and local box to new size/shape
// This operation corresponds to applying the
// this operation corresponds to applying the
// translate and scale operations
// corresponding to the solution of the following ODE:
//
@ -919,65 +988,119 @@ void FixNH::remap()
double dto4 = dto/4.0;
double dto8 = dto/8.0;
// off-diagonal components, first half
if (pstyle == TRICLINIC) {
h[4] *= exp(dto8*omega_dot[0]);
h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
h[4] *= exp(dto8*omega_dot[0]);
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;
}
h[3] *= exp(dto4*omega_dot[1]);
h[3] += dto2*(omega_dot[3]*h[2]);
h[3] *= exp(dto4*omega_dot[1]);
if (p_flag[3]) {
expfac = exp(dto4*omega_dot[1]);
h[3] *= expfac;
h[3] += dto2*(omega_dot[3]*h[2]);
h[3] *= expfac;
}
h[5] *= exp(dto4*omega_dot[0]);
h[5] += dto2*(omega_dot[5]*h[1]);
h[5] *= exp(dto4*omega_dot[0]);
if (p_flag[5]) {
expfac = exp(dto4*omega_dot[0]);
h[5] *= expfac;
h[5] += dto2*(omega_dot[5]*h[1]);
h[5] *= expfac;
}
h[4] *= exp(dto8*omega_dot[0]);
h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
h[4] *= exp(dto8*omega_dot[0]);
}
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)*exp(dto*omega_dot[i]) + ctr;
domain->boxhi[i] = (oldhi-ctr)*exp(dto*omega_dot[i]) + ctr;
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];
ctr = 0.5 * (oldlo + oldhi);
expfac = exp(dto*omega_dot[0]);
domain->boxlo[0] = (oldlo-ctr)*expfac + ctr;
domain->boxhi[0] = (oldhi-ctr)*expfac + ctr;
}
if (p_flag[1]) {
oldlo = domain->boxlo[1];
oldhi = domain->boxhi[1];
ctr = 0.5 * (oldlo + oldhi);
expfac = exp(dto*omega_dot[1]);
domain->boxlo[1] = (oldlo-ctr)*expfac + ctr;
domain->boxhi[1] = (oldhi-ctr)*expfac + ctr;
if (scalexy) h[5] *= expfac;
}
if (p_flag[2]) {
oldlo = domain->boxlo[2];
oldhi = domain->boxhi[2];
ctr = 0.5 * (oldlo + oldhi);
expfac = exp(dto*omega_dot[2]);
domain->boxlo[2] = (oldlo-ctr)*expfac + ctr;
domain->boxhi[2] = (oldhi-ctr)*expfac + ctr;
if (scalexz) h[4] *= expfac;
if (scaleyz) h[3] *= expfac;
}
// off-diagonal components, second half
if (pstyle == TRICLINIC) {
h[4] *= exp(dto8*omega_dot[0]);
h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
h[4] *= exp(dto8*omega_dot[0]);
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;
}
h[3] *= exp(dto4*omega_dot[1]);
h[3] += dto2*(omega_dot[3]*h[2]);
h[3] *= exp(dto4*omega_dot[1]);
if (p_flag[3]) {
expfac = exp(dto4*omega_dot[1]);
h[3] *= expfac;
h[3] += dto2*(omega_dot[3]*h[2]);
h[3] *= expfac;
}
h[5] *= exp(dto4*omega_dot[0]);
h[5] += dto2*(omega_dot[5]*h[1]);
h[5] *= exp(dto4*omega_dot[0]);
if (p_flag[5]) {
expfac = exp(dto4*omega_dot[0]);
h[5] *= expfac;
h[5] += dto2*(omega_dot[5]*h[1]);
h[5] *= expfac;
}
h[4] *= exp(dto8*omega_dot[0]);
h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
h[4] *= exp(dto8*omega_dot[0]);
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];
if (domain->yz < -0.5*domain->yprd || domain->yz > 0.5*domain->yprd ||
domain->xz < -0.5*domain->xprd || domain->xz > 0.5*domain->xprd ||
domain->xy < -0.5*domain->xprd || domain->xy > 0.5*domain->xprd)
error->all("Fix npt/nph has tilted box too far - "
"box flips are not yet implemented");
}
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("Fix npt/nph has tilted box too far in one step - "
"periodic cell is too far from equilibrium state");
domain->set_global_box();
domain->set_local_box();
@ -1940,3 +2063,66 @@ void FixNH::nh_omega_dot()
}
}
}
/* ----------------------------------------------------------------------
if box tilt exceeds limits,
create new box in domain
remap to put far-away atoms back into new box
perform irregular on atoms in lamda coords to get atoms to new procs
force reneighboring on next timestep
------------------------------------------------------------------------- */
void FixNH::pre_exchange()
{
double xprd = domain->xprd;
double yprd = domain->yprd;
// flip is triggered when tilt exceeds 0.5 by
// an amount DELTAFLIP that is somewhat arbitrary
double xtiltmax = (0.5+DELTAFLIP)*xprd;
double ytiltmax = (0.5+DELTAFLIP)*yprd;
int flip = 0;
if (domain->yz < -ytiltmax) {
flip = 1;
domain->yz += yprd;
domain->xz += domain->xy;
} else if (domain->yz >= ytiltmax) {
flip = 1;
domain->yz -= yprd;
domain->xz -= domain->xy;
}
if (domain->xz < -xtiltmax) {
flip = 1;
domain->xz += xprd;
} else if (domain->xz >= xtiltmax) {
flip = 1;
domain->xz -= xprd;
}
if (domain->xy < -xtiltmax) {
flip = 1;
domain->xy += xprd;
} else if (domain->xy >= xtiltmax) {
flip = 1;
domain->xy -= xprd;
}
if (flip) {
domain->set_global_box();
domain->set_local_box();
double **x = atom->x;
int *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);
}
}

7
src/fix_nh.h
View File

@ -29,6 +29,7 @@ class FixNH : public Fix {
virtual void final_integrate();
void initial_integrate_respa(int, int, int);
void final_integrate_respa(int, int);
void pre_exchange();
double compute_scalar();
double compute_vector(int);
void write_restart(FILE *);
@ -62,6 +63,7 @@ class FixNH : public Fix {
int kspace_flag; // 1 if KSpace invoked, 0 if not
int nrigid; // number of rigid fixes
int *rfix; // indices of rigid fixes
class Irregular *irregular; // for migrating atoms after box flips
int nlevels_respa;
double *step_respa;
@ -95,8 +97,13 @@ class FixNH : public Fix {
int deviatoric_flag; // 0 if target stress tensor is hydrostatic
double h0_inv[6]; // h_inv of reference (zero strain) box
int nreset_h0; // interval for resetting h0
double mtk_term1,mtk_term2; // Martyna-Tobias-Klein corrections
int scaleyz; // 1 if yz scaled with lz
int scalexz; // 1 if xz scaled with lz
int scalexy; // 1 if xy scaled with ly
void couple();
void remap();
void nhc_temp_integrate();

104
src/thermo.cpp
View File

@ -49,6 +49,7 @@ using namespace LAMMPS_NS;
// vol, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz, xlat, ylat, zlat
// pxx, pyy, pzz, pxy, pxz, pyz
// fmax, fnorm
// cella, cellb, cellc, cellalpha, cellbeta, cellgamma
// customize a new thermo style by adding a DEFINE to this list
@ -149,6 +150,8 @@ Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
format_float_user = NULL;
format_int_user = NULL;
format_bigint_user = NULL;
PI = 4.0*atan(1.0);
}
/* ---------------------------------------------------------------------- */
@ -768,6 +771,19 @@ void Thermo::parse_fields(char *str)
} else if (strcmp(word,"fnorm") == 0) {
addfield("Fnorm",&Thermo::compute_fnorm,FLOAT);
} else if (strcmp(word,"cella") == 0) {
addfield("Cella",&Thermo::compute_cella,FLOAT);
} else if (strcmp(word,"cellb") == 0) {
addfield("Cellb",&Thermo::compute_cellb,FLOAT);
} else if (strcmp(word,"cellc") == 0) {
addfield("Cellc",&Thermo::compute_cellc,FLOAT);
} else if (strcmp(word,"cellalpha") == 0) {
addfield("CellAlpha",&Thermo::compute_cellalpha,FLOAT);
} else if (strcmp(word,"cellbeta") == 0) {
addfield("CellBeta",&Thermo::compute_cellbeta,FLOAT);
} else if (strcmp(word,"cellgamma") == 0) {
addfield("CellGamma",&Thermo::compute_cellgamma,FLOAT);
// compute value = c_ID, fix value = f_ID, variable value = v_ID
// count trailing [] and store int arguments
// copy = at most 8 chars of ID to pass to addfield
@ -1292,6 +1308,13 @@ int Thermo::evaluate_keyword(char *word, double *answer)
} else if (strcmp(word,"fmax") == 0) compute_fmax();
else if (strcmp(word,"fnorm") == 0) compute_fnorm();
else if (strcmp(word,"cella") == 0) compute_cella();
else if (strcmp(word,"cellb") == 0) compute_cellb();
else if (strcmp(word,"cellc") == 0) compute_cellc();
else if (strcmp(word,"cellalpha") == 0) compute_cellalpha();
else if (strcmp(word,"cellbeta") == 0) compute_cellbeta();
else if (strcmp(word,"cellgamma") == 0) compute_cellgamma();
else return 1;
*answer = dvalue;
@ -1821,3 +1844,84 @@ void Thermo::compute_fnorm()
MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
dvalue = sqrt(dotall);
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cella()
{
dvalue = domain->xprd;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cellb()
{
if (!domain->triclinic)
dvalue = domain->yprd;
else {
double* h = domain->h;
dvalue = sqrt(h[1]*h[1]+h[5]*h[5]);
}
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cellc()
{
if (!domain->triclinic)
dvalue = domain->zprd;
else {
double* h = domain->h;
dvalue = sqrt(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]);
}
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cellalpha()
{
if (!domain->triclinic)
dvalue = 90.0;
else {
// Cos(alpha) = (xy.xz + ly.yz)/(b.c)
double* h = domain->h;
double cosalpha = (h[5]*h[4]+h[1]*h[3])/
sqrt((h[1]*h[1]+h[5]*h[5])*(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]));
dvalue = acos(cosalpha)*180.0/PI;
}
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cellbeta()
{
if (!domain->triclinic)
dvalue = 90.0;
else {
// Cos(beta) = xz/c
double* h = domain->h;
double cosbeta = h[4]/sqrt(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]);
dvalue = acos(cosbeta)*180.0/PI;
}
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cellgamma()
{
if (!domain->triclinic)
dvalue = 90.0;
else {
// Cos(gamma) = xy/b
double* h = domain->h;
double cosgamma = h[5]/sqrt(h[1]*h[1]+h[5]*h[5]);
dvalue = acos(cosgamma)*180.0/PI;
}
}

9
src/thermo.h
View File

@ -96,6 +96,8 @@ class Thermo : protected Pointers {
char **id_variable; // list of variable names
int *variables; // list of Variable indices
double PI;
// private methods
void allocate();
@ -173,6 +175,13 @@ class Thermo : protected Pointers {
void compute_fmax();
void compute_fnorm();
void compute_cella();
void compute_cellb();
void compute_cellc();
void compute_cellalpha();
void compute_cellbeta();
void compute_cellgamma();
};
}