From dfb085b3215f4e146533c3a7d85d1a6c17803031 Mon Sep 17 00:00:00 2001 From: athomps Date: Thu, 16 Jun 2011 16:30:04 +0000 Subject: [PATCH] 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 --- src/fix_deform.cpp | 4 +- src/fix_nh.cpp | 286 +++++++++++++++++++++++++++++++++++++-------- src/fix_nh.h | 7 ++ src/thermo.cpp | 104 +++++++++++++++++ src/thermo.h | 9 ++ 5 files changed, 358 insertions(+), 52 deletions(-) diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index ac8a64ba01..75db867299 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -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; diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index cf4d5535fb..77650a9784 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -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); + } +} + diff --git a/src/fix_nh.h b/src/fix_nh.h index 58d4f6fb06..bdbef59824 100644 --- a/src/fix_nh.h +++ b/src/fix_nh.h @@ -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(); diff --git a/src/thermo.cpp b/src/thermo.cpp index afb32ac9a2..d128d5bb1f 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -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; + } +} + diff --git a/src/thermo.h b/src/thermo.h index 49e2f79e57..352209cb35 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -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(); }; }