From d615d8053b5e4362b383aa66d26b050281f593af Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 5 Apr 2024 14:31:53 -0600 Subject: [PATCH] support for general tri edge vectors in thermo output --- doc/src/Howto_triclinic.rst | 21 +- doc/src/thermo_style.rst | 71 ++++-- examples/README | 1 + src/domain.cpp | 5 +- src/thermo.cpp | 427 ++++++++++++++++++++++-------------- src/thermo.h | 34 ++- 6 files changed, 355 insertions(+), 204 deletions(-) diff --git a/doc/src/Howto_triclinic.rst b/doc/src/Howto_triclinic.rst index 2f2ffe85cf..9af941e4a6 100644 --- a/doc/src/Howto_triclinic.rst +++ b/doc/src/Howto_triclinic.rst @@ -94,19 +94,20 @@ restricted triclinic parallelepiped. simulation boxes in LAMMPS. Note that the :doc:`thermo_style custom ` command has -keywords for outputting the various parameters that define both -restricted and general triclinic simulation boxes. Thus you can check -the restricted triclinic box parameters LAMMPS generates to rotate a -general triclinic box to restricted triclinic form. +keywords for outputting the various parameters that define the size +and shape of orthogonal, restricted triclinic, and general triclinic +simulation boxes. + +For orthogonal boxes there these are the 6 thermo keywords +(xlo,ylo,zlo) and (xhi,yhi,zhi). For restricted triclinic boxes these are the 9 thermo keywords for -(xlo,ylo,zlo), (xhi,yhi,zhi), and the (xy,xz,yz) tilt factors. For -general triclinic boxes these are the 12 thermo keywords for +(xlo,ylo,zlo), (xhi,yhi,zhi), and the (xy,xz,yz) tilt factors. + +For general triclinic boxes these are the 12 thermo keywords for (xlo,ylo,zhi) and the components of the **A**, **B**, **C** edge -vectors. For both orthogonal and restricted triclinic boxes, the -thermo keywords lx/ly/lz refer to the box sizes, namely lx = xhi - -xlo, etc. Lx,ly,lz are the box edge vector lengths for orthogonal and -restricted/general triclinic simulation boxes. +vectors, namely (avecx,avecy,avecz), (bvecx,bvecy,bvecz), and +(cvecx,cvecy,cvecz), The remainder of this doc page explains (a) how LAMMPS operates with general triclinic simulation boxes, (b) mathematical transformations diff --git a/doc/src/thermo_style.rst b/doc/src/thermo_style.rst index f73c4baa3d..cdb39fe8a5 100644 --- a/doc/src/thermo_style.rst +++ b/doc/src/thermo_style.rst @@ -25,12 +25,18 @@ Syntax evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail, enthalpy, ecouple, econserve, - vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, - xy, xz, yz, xlat, ylat, zlat, - bonds, angles, dihedrals, impropers, - pxx, pyy, pzz, pxy, pxz, pyz, - fmax, fnorm, nbuild, ndanger, + vol, density, + xlo, xhi, ylo, yhi, zlo, zhi, + xy, xz, yz, + avecx, avecy, avecz, + bvecx, bvecy, bvecz, + cvecx, cvecy, cvecz, + lx, ly, lz, + xlat, ylat, zlat, cella, cellb, cellc, cellalpha, cellbeta, cellgamma, + pxx, pyy, pzz, pxy, pxz, pyz, + bonds, angles, dihedrals, impropers, + fmax, fnorm, nbuild, ndanger, c_ID, c_ID[I], c_ID[I][J], f_ID, f_ID[I], f_ID[I][J], v_name, v_name[I] @@ -66,18 +72,21 @@ Syntax econserve = pe + ke + ecouple = etotal + ecouple vol = volume density = mass density of system - lx,ly,lz = box lengths in x,y,z xlo,xhi,ylo,yhi,zlo,zhi = box boundaries xy,xz,yz = box tilt for restricted triclinic (non-orthogonal) simulation boxes + avecx,avecy,avecz = components of edge vector A for general triclinic simulation boxes + bvecx,bvecy,bvecz = components of edge vector B for general triclinic simulation boxes + cvecx,cvecy,cvecz = components of edge vector C for general triclinic simulation boxes + lx,ly,lz = box lengths in x,y,z xlat,ylat,zlat = lattice spacings as calculated by :doc:`lattice ` command - bonds,angles,dihedrals,impropers = # of these interactions defined + cella,cellb,cellc = periodic cell lattice constants a,b,c + cellalpha, cellbeta, cellgamma = periodic cell angles alpha,beta,gamma pxx,pyy,pzz,pxy,pxz,pyz = 6 components of pressure tensor + bonds,angles,dihedrals,impropers = # of these interactions defined fmax = max component of force on any atom in any dimension fnorm = length of force vector for all atoms nbuild = # of neighbor list builds ndanger = # of dangerous neighbor list builds - cella,cellb,cellc = periodic cell lattice constants a,b,c - cellalpha, cellbeta, cellgamma = periodic cell angles alpha,beta,gamma c_ID = global scalar value calculated by a compute with ID c_ID[I] = Ith component of global vector calculated by a compute with ID, I can include wildcard (see below) c_ID[I][J] = I,J component of global array calculated by a compute with ID @@ -248,7 +257,7 @@ and *pxx*, *pyy*, etc. ---------- Here is more information on other keywords whose meaning may not be -clear: +clear. The *step*, *elapsed*, and *elaplong* keywords refer to timestep count. *Step* is the current timestep, or iteration count when a @@ -322,6 +331,38 @@ thermostatting or barostatting to their coupling reservoirs -- that is, the NVT, NPH, or NPT ensembles, the *econserve* quantity should remain constant over time even though *etotal* may change. +In LAMMPS, the simulation box can be defined as orthogonal or +triclinic (non-orthogonal). See the :doc:`Howto_triclinic +` doc page for a detailed explanation of orthogonal, +restricted triclinic, and general triclinic simulation boxes and how +LAMMPS rotates a general triclinic box to be restricted triclinic +internally. + +The *lx*, *ly*, *lz* keywords are the extent of the simulation box in +each dimension. The *xlo*, *xhi*, *ylo*, *yhi*, *zlo*, *zhi* keywords +are the lower and upper bounds of the simulation box in each +dimension. I.e. *lx* = *xhi* - *xlo). These 9 values are the same +for all 3 kinds of boxes. I.e. for a restricted triclinic box, they +are the values as if the box were not tilted. For a general triclinic +box, they are the values after it is internally rotated to be a +restricted triclinic box. + +THe *xy*, *xz*, *yz* are the current tilt factors for a triclinic box. +They are the same for restricted and general triclinic boxes. + +The *avecx*, *avecy*, *avecz*, *bvecx*, *bvecy*, *bvecz*, *cvecx*, +*cvecy*, *cvecz* are the components of the 3 edge vectors of the +current general triclinic box. They are only non-zero if a general +triclinic box was defined when the simultion box was created. + +The *cella*, *cellb*, *cellc*, *cellalpha*, *cellbeta*, *cellgamma* +keywords correspond to the usual crystallographic quantities that +define the periodic simulation box of a crystalline system. See the +:doc:`Howto triclinic ` page for a precise definition +of these quantities in terms of the LAMMPS representation of a +restricted triclinic simulation box via *lx*, *ly*, *lz*, *yz*, *xz*, +*xy*\ . + The *pxx,pyy,pzz,pxy,pxz,pyz* keywords are the 6 components of the symmetric pressure tensor for the system. See the :doc:`compute pressure ` command doc page for details of how it is @@ -329,7 +370,7 @@ calculated. If the :doc:`thermo_modify triclinic/general ` option is set and the simulation box was created as a general triclinic box, -then the components will be output as values consistent with the +then the 6 components will be output as values consistent with the orientation of the general triclinic box relative to the standard xyz coordinate axes. If this keyword is not used, the values will be consistent with the orientation of the restricted triclinic box (which @@ -362,14 +403,6 @@ to reduce the delay factor to ensure no force interactions are missed by atoms moving beyond the neighbor skin distance before a rebuild takes place. -The keywords *cella*, *cellb*, *cellc*, *cellalpha*, *cellbeta*, -*cellgamma*, correspond to the usual crystallographic quantities that -define the periodic unit cell of a crystal. See the :doc:`Howto -triclinic ` page for a geometric description of -restricted triclinic periodic cells, including a precise definition of -these quantities in terms of the internal LAMMPS cell dimensions *lx*, -*ly*, *lz*, *yz*, *xz*, *xy*\ . - ---------- For output values from a compute or fix or variable, the bracketed diff --git a/examples/README b/examples/README index 62a09f654d..f76dced3e4 100644 --- a/examples/README +++ b/examples/README @@ -115,6 +115,7 @@ tad: temperature-accelerated dynamics of vacancy diffusion in bulk Si template: examples for using atom_style template and comparing to atom style molecular tersoff: regression test input for Tersoff variants threebody: regression test input for a variety of threebody potentials +triclinic: general and restricted triclinic examples ttm: two-temeperature model examples vashishta: models using the Vashishta potential voronoi: Voronoi tesselation via compute voronoi/atom command diff --git a/src/domain.cpp b/src/domain.cpp index 1a553061fb..adc85b68b9 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -87,7 +87,10 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp) boxlo[0] = boxlo[1] = boxlo[2] = -0.5; boxhi[0] = boxhi[1] = boxhi[2] = 0.5; xy = xz = yz = 0.0; - + avec[0] = avec[1] = avec[2] = 0.0; + bvec[0] = bvec[1] = bvec[2] = 0.0; + cvec[0] = cvec[1] = cvec[2] = 0.0; + h[3] = h[4] = h[5] = 0.0; h_inv[3] = h_inv[4] = h_inv[5] = 0.0; h_rate[0] = h_rate[1] = h_rate[2] = diff --git a/src/thermo.cpp b/src/thermo.cpp index 999932c407..4ba29bbc73 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -56,17 +56,18 @@ using namespace MathExtra; // CUSTOMIZATION: add a new keyword by adding it to this list: -// step, elapsed, elaplong, dt, time, cpu, tpcpu, spcpu, cpuremain, -// part, timeremain +// step, elapsed, elaplong, dt, time, cpu, tpcpu, spcpu, cpuremain, part, timeremain // atoms, temp, press, pe, ke, etotal // evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail // enthalpy, ecouple, econserve -// vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz +// vol, density, lx, ly, lz, +// xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz +// avecx, avecy, avecz, bvecx, bvecy, bvecz, cvecx, cvecy, cvecz, // xlat, ylat, zlat -// bonds, angles, dihedrals, impropers -// pxx, pyy, pzz, pxy, pxz, pyz -// fmax, fnorm, nbuild, ndanger // cella, cellb, cellc, cellalpha, cellbeta, cellgamma +// pxx, pyy, pzz, pxy, pxz, pyz +// bonds, angles, dihedrals, impropers +// fmax, fnorm, nbuild, ndanger // CUSTOMIZATION: add a new thermo style by adding a constant to the enumerator, // define a new string constant with the keywords and provide default formats. @@ -916,6 +917,7 @@ void Thermo::parse_fields(const std::string &str) addfield("Volume", &Thermo::compute_vol, FLOAT); } else if (word == "density") { addfield("Density", &Thermo::compute_density, FLOAT); + } else if (word == "lx") { addfield("Lx", &Thermo::compute_lx, FLOAT); } else if (word == "ly") { @@ -943,6 +945,25 @@ void Thermo::parse_fields(const std::string &str) } else if (word == "yz") { addfield("Yz", &Thermo::compute_yz, FLOAT); + } else if (word == "avecx") { + addfield("Avecx", &Thermo::compute_avecx, FLOAT); + } else if (word == "avecy") { + addfield("Avecy", &Thermo::compute_avecy, FLOAT); + } else if (word == "avecz") { + addfield("Avecz", &Thermo::compute_avecz, FLOAT); + } else if (word == "bvecx") { + addfield("Bvecx", &Thermo::compute_bvecx, FLOAT); + } else if (word == "bvecy") { + addfield("Bvecy", &Thermo::compute_bvecy, FLOAT); + } else if (word == "bvecz") { + addfield("Bvecz", &Thermo::compute_bvecz, FLOAT); + } else if (word == "cvecx") { + addfield("Cvecx", &Thermo::compute_cvecx, FLOAT); + } else if (word == "cvecy") { + addfield("Cvecy", &Thermo::compute_cvecy, FLOAT); + } else if (word == "cvecz") { + addfield("Cvecz", &Thermo::compute_cvecz, FLOAT); + } else if (word == "xlat") { addfield("Xlat", &Thermo::compute_xlat, FLOAT); } else if (word == "ylat") { @@ -950,14 +971,18 @@ void Thermo::parse_fields(const std::string &str) } else if (word == "zlat") { addfield("Zlat", &Thermo::compute_zlat, FLOAT); - } else if (word == "bonds") { - addfield("Bonds", &Thermo::compute_bonds, BIGINT); - } else if (word == "angles") { - addfield("Angles", &Thermo::compute_angles, BIGINT); - } else if (word == "dihedrals") { - addfield("Diheds", &Thermo::compute_dihedrals, BIGINT); - } else if (word == "impropers") { - addfield("Impros", &Thermo::compute_impropers, BIGINT); + } else if (word == "cella") { + addfield("Cella", &Thermo::compute_cella, FLOAT); + } else if (word == "cellb") { + addfield("Cellb", &Thermo::compute_cellb, FLOAT); + } else if (word == "cellc") { + addfield("Cellc", &Thermo::compute_cellc, FLOAT); + } else if (word == "cellalpha") { + addfield("CellAlpha", &Thermo::compute_cellalpha, FLOAT); + } else if (word == "cellbeta") { + addfield("CellBeta", &Thermo::compute_cellbeta, FLOAT); + } else if (word == "cellgamma") { + addfield("CellGamma", &Thermo::compute_cellgamma, FLOAT); } else if (word == "pxx") { if (triclinic_general) @@ -990,6 +1015,15 @@ void Thermo::parse_fields(const std::string &str) addfield("Pyz", &Thermo::compute_pyz, FLOAT); index_press_vector = add_compute(id_press, VECTOR); + } else if (word == "bonds") { + addfield("Bonds", &Thermo::compute_bonds, BIGINT); + } else if (word == "angles") { + addfield("Angles", &Thermo::compute_angles, BIGINT); + } else if (word == "dihedrals") { + addfield("Diheds", &Thermo::compute_dihedrals, BIGINT); + } else if (word == "impropers") { + addfield("Impros", &Thermo::compute_impropers, BIGINT); + } else if (word == "fmax") { addfield("Fmax", &Thermo::compute_fmax, FLOAT); } else if (word == "fnorm") { @@ -1000,19 +1034,6 @@ void Thermo::parse_fields(const std::string &str) } else if (word == "ndanger") { addfield("Ndanger", &Thermo::compute_ndanger, BIGINT); - } else if (word == "cella") { - addfield("Cella", &Thermo::compute_cella, FLOAT); - } else if (word == "cellb") { - addfield("Cellb", &Thermo::compute_cellb, FLOAT); - } else if (word == "cellc") { - addfield("Cellc", &Thermo::compute_cellc, FLOAT); - } else if (word == "cellalpha") { - addfield("CellAlpha", &Thermo::compute_cellalpha, FLOAT); - } else if (word == "cellbeta") { - addfield("CellBeta", &Thermo::compute_cellbeta, FLOAT); - } else if (word == "cellgamma") { - addfield("CellGamma", &Thermo::compute_cellgamma, FLOAT); - // compute value = c_ID, fix value = f_ID, variable value = v_ID // count trailing [] and store int arguments @@ -1316,22 +1337,6 @@ int Thermo::evaluate_keyword(const std::string &word, double *answer) compute_atoms(); dvalue = bivalue; - } else if (word == "bonds") { - compute_bonds(); - dvalue = bivalue; - - } else if (word == "angles") { - compute_angles(); - dvalue = bivalue; - - } else if (word == "dihedrals") { - compute_dihedrals(); - dvalue = bivalue; - - } else if (word == "impropers") { - compute_impropers(); - dvalue = bivalue; - } else if (word == "temp") { check_temp(word); compute_temp(); @@ -1412,6 +1417,7 @@ int Thermo::evaluate_keyword(const std::string &word, double *answer) compute_vol(); else if (word == "density") compute_density(); + else if (word == "lx") compute_lx(); else if (word == "ly") @@ -1439,6 +1445,25 @@ int Thermo::evaluate_keyword(const std::string &word, double *answer) else if (word == "yz") compute_yz(); + else if (word == "avecx") + compute_avecx(); + else if (word == "avecy") + compute_avecy(); + else if (word == "avecz") + compute_avecz(); + else if (word == "bvecx") + compute_bvecx(); + else if (word == "bvecy") + compute_bvecy(); + else if (word == "bvecz") + compute_bvecz(); + else if (word == "cvecx") + compute_cvecx(); + else if (word == "cvecy") + compute_cvecy(); + else if (word == "cvecz") + compute_cvecz(); + else if (word == "xlat") compute_xlat(); else if (word == "ylat") @@ -1446,6 +1471,19 @@ int Thermo::evaluate_keyword(const std::string &word, double *answer) else if (word == "zlat") compute_zlat(); + else if (word == "cella") + compute_cella(); + else if (word == "cellb") + compute_cellb(); + else if (word == "cellc") + compute_cellc(); + else if (word == "cellalpha") + compute_cellalpha(); + else if (word == "cellbeta") + compute_cellbeta(); + else if (word == "cellgamma") + compute_cellgamma(); + else if (word == "pxx") { check_press_vector(word); if (triclinic_general) compute_pxx_triclinic_general(); @@ -1475,9 +1513,24 @@ int Thermo::evaluate_keyword(const std::string &word, double *answer) check_press_vector(word); if (triclinic_general) compute_pyz_triclinic_general(); else compute_pyz(); - } - else if (word == "fmax") + } else if (word == "bonds") { + compute_bonds(); + dvalue = bivalue; + + } else if (word == "angles") { + compute_angles(); + dvalue = bivalue; + + } else if (word == "dihedrals") { + compute_dihedrals(); + dvalue = bivalue; + + } else if (word == "impropers") { + compute_impropers(); + dvalue = bivalue; + + } else if (word == "fmax") compute_fmax(); else if (word == "fnorm") compute_fnorm(); @@ -1490,19 +1543,6 @@ int Thermo::evaluate_keyword(const std::string &word, double *answer) dvalue = bivalue; } - else if (word == "cella") - compute_cella(); - else if (word == "cellb") - compute_cellb(); - else if (word == "cellc") - compute_cellc(); - else if (word == "cellalpha") - compute_cellalpha(); - else if (word == "cellbeta") - compute_cellbeta(); - else if (word == "cellgamma") - compute_cellgamma(); - else return 1; @@ -1774,23 +1814,6 @@ void Thermo::compute_etotal() /* ---------------------------------------------------------------------- */ -void Thermo::compute_ecouple() -{ - dvalue = modify->energy_couple(); -} - -/* ---------------------------------------------------------------------- */ - -void Thermo::compute_econserve() -{ - compute_etotal(); - double dvalue_etotal = dvalue; - compute_ecouple(); - dvalue += dvalue_etotal; -} - -/* ---------------------------------------------------------------------- */ - void Thermo::compute_evdwl() { double tmp = 0.0; @@ -1938,6 +1961,23 @@ void Thermo::compute_enthalpy() /* ---------------------------------------------------------------------- */ +void Thermo::compute_ecouple() +{ + dvalue = modify->energy_couple(); +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_econserve() +{ + compute_etotal(); + double dvalue_etotal = dvalue; + compute_ecouple(); + dvalue += dvalue_etotal; +} + +/* ---------------------------------------------------------------------- */ + void Thermo::compute_vol() { if (domain->dimension == 3) @@ -2041,6 +2081,67 @@ void Thermo::compute_yz() /* ---------------------------------------------------------------------- */ +void Thermo::compute_avecx() +{ + dvalue = domain->avec[0]; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_avecy() +{ + dvalue = domain->avec[1]; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_avecz() +{ + dvalue = domain->avec[2]; +} +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_bvecx() +{ + dvalue = domain->bvec[0]; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_bvecy() +{ + dvalue = domain->bvec[1]; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_bvecz() +{ + dvalue = domain->bvec[2]; +} +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_cvecx() +{ + dvalue = domain->cvec[0]; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_cvecy() +{ + dvalue = domain->cvec[1]; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_cvecz() +{ + dvalue = domain->cvec[2]; +} + +/* ---------------------------------------------------------------------- */ + void Thermo::compute_xlat() { dvalue = domain->lattice->xlattice; @@ -2062,30 +2163,82 @@ void Thermo::compute_zlat() /* ---------------------------------------------------------------------- */ -void Thermo::compute_bonds() +void Thermo::compute_cella() { - bivalue = atom->nbonds; + dvalue = domain->xprd; } /* ---------------------------------------------------------------------- */ -void Thermo::compute_angles() +void Thermo::compute_cellb() { - bivalue = atom->nangles; + if (!domain->triclinic) + dvalue = domain->yprd; + else { + double *h = domain->h; + dvalue = sqrt(h[1] * h[1] + h[5] * h[5]); + } } /* ---------------------------------------------------------------------- */ -void Thermo::compute_dihedrals() +void Thermo::compute_cellc() { - bivalue = atom->ndihedrals; + 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_impropers() +void Thermo::compute_cellalpha() { - bivalue = atom->nimpropers; + 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 / MY_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 / MY_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 / MY_PI; + } } /* ---------------------------------------------------------------------- */ @@ -2192,6 +2345,34 @@ void Thermo::compute_pyz_triclinic_general() /* ---------------------------------------------------------------------- */ +void Thermo::compute_bonds() +{ + bivalue = atom->nbonds; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_angles() +{ + bivalue = atom->nangles; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_dihedrals() +{ + bivalue = atom->ndihedrals; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_impropers() +{ + bivalue = atom->nimpropers; +} + +/* ---------------------------------------------------------------------- */ + void Thermo::compute_fmax() { double **f = atom->f; @@ -2235,83 +2416,3 @@ void Thermo::compute_ndanger() { bivalue = neighbor->ndanger; } - -/* ---------------------------------------------------------------------- */ - -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 / MY_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 / MY_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 / MY_PI; - } -} diff --git a/src/thermo.h b/src/thermo.h index a82d462585..31036f55d3 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -176,6 +176,7 @@ class Thermo : protected Pointers { void compute_vol(); void compute_density(); + void compute_lx(); void compute_ly(); void compute_lz(); @@ -191,14 +192,26 @@ class Thermo : protected Pointers { void compute_xz(); void compute_yz(); + void compute_avecx(); + void compute_avecy(); + void compute_avecz(); + void compute_bvecx(); + void compute_bvecy(); + void compute_bvecz(); + void compute_cvecx(); + void compute_cvecy(); + void compute_cvecz(); + void compute_xlat(); void compute_ylat(); void compute_zlat(); - void compute_bonds(); - void compute_angles(); - void compute_dihedrals(); - void compute_impropers(); + void compute_cella(); + void compute_cellb(); + void compute_cellc(); + void compute_cellalpha(); + void compute_cellbeta(); + void compute_cellgamma(); void compute_pxx(); void compute_pyy(); @@ -211,8 +224,13 @@ class Thermo : protected Pointers { void compute_pyy_triclinic_general(); void compute_pzz_triclinic_general(); void compute_pxy_triclinic_general(); - void compute_pyz_triclinic_general(); void compute_pxz_triclinic_general(); + void compute_pyz_triclinic_general(); + + void compute_bonds(); + void compute_angles(); + void compute_dihedrals(); + void compute_impropers(); void compute_fmax(); void compute_fnorm(); @@ -220,12 +238,6 @@ class Thermo : protected Pointers { void compute_nbuild(); void compute_ndanger(); - void compute_cella(); - void compute_cellb(); - void compute_cellc(); - void compute_cellalpha(); - void compute_cellbeta(); - void compute_cellgamma(); }; } // namespace LAMMPS_NS