ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

support for general tri edge vectors in thermo output

Browse Source
This commit is contained in:
Steve Plimpton
2024-04-05 14:31:53 -06:00
parent 905ceb98f1
commit d615d8053b
6 changed files with 355 additions and 204 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
doc/src/Howto_triclinic.rst
View File

@ -94,19 +94,20 @@ restricted triclinic parallelepiped.
simulation boxes in LAMMPS.
Note that the :doc:`thermo_style custom <thermo_style>` 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

71
doc/src/thermo_style.rst
View File

@ -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 <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
<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 <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 <compute_pressure>` command doc page for details of how it is
@ -329,7 +370,7 @@ calculated.
If the :doc:`thermo_modify triclinic/general <thermo_modify>` 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 <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

1
examples/README
View File

@ -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

5
src/domain.cpp
View File

@ -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] =

427
src/thermo.cpp
View File

@ -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;
}
}

34
src/thermo.h
View File

@ -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