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

Two computes that computes the local mechanical pressure tensor in Cartesian and spherical coordinates

Browse Source
This commit is contained in:
Olav Galteland
2022-03-04 16:04:18 +01:00
parent 798975b809
commit 760d85b9c4
6 changed files with 1214 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

63
doc/src/compute_pressure_cartesian.rst Normal file
View File

@ -0,0 +1,63 @@
.. index:: compute pressure/cartesian
compute pressure/cartesian command
==================================
Syntax
""""""
.. parsed-literal::
compute ID group-ID pressure/cartesian dim bin_width
* ID, group-ID are documented in :doc:`compute <compute>` command
* pressure/cartesian = style name of this compute command
* one or two dim/bin_width pairs may be appended
* dim = x, y, or z
* bin_width = width of the bin
Examples
""""""""
.. code-block:: LAMMPS
compute 1 all pressure/cartesian x 0.1
compute 1 all pressure/cartesian y 0.25 z 0.1
Description
"""""""""""
This computes the Cartesian pressure tensor in one or two dimensions, as described in :ref:`(Ikeshoji)<Ikeshoji2>`. This can for example be used to calculate the local pressure tensor of flat liquid-vapor interfaces. This compute obeys momentum balance, such that the pressure tensor component normal to a surface is constant. This compute uses the Irving-Kirkwood contour, which is the straight line between particle pairs. The pressure tensor is split into a kinetic contribution :math:`P^k` and a virial contribution :math:`P^v`. The sum gives the total pressure tensor :math:`P = P^k+P^v`.
Output info
"""""""""""
The output is a global array with 8 columns when only one dimension is specified and 9 columns when two dimensions are specified. There are (L1*L2)/(bin_width1*bin_width2) rows, where L1 and L2 are the size of the simulation box in the relevant dimensions. The output columns are position of the center of the local volume in the first and second dimension (when two dimensions are specified), number density, :math:`P^k_{xx}`, :math:`P^k_{yy}`, :math:`P^k_{zz}`, :math:`P^v_{xx}`, :math:`P^v_{yy}`, and :math:`P^v_{zz}`.
This array can be output with the :doc:`fix ave/time <fix_ave_time>`,
.. code-block:: LAMMPS
compute p all pressure/cartesian x 0.1
fix 2 all ave/time 100 1 100 c_p[*] file dump_p.out mode vector
Restrictions
""""""""""""
This compute is part of the EXTRA-COMPUTE package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
Related commands
""""""""""""""""
:doc:`compute stress/atom <compute_stress_atom>`, :doc:`compute pressure/cylinder <compute_pressure_cylinder>`, :doc:`compute pressure/spherical <compute_pressure_spherical>`
Default
"""""""
none
----------
.. _Ikeshoji2:
**(Ikeshoji)** Ikeshoji, Hafskjold, Furuholt, Mol Sim, 29, 101-109, (2003).

61
doc/src/compute_pressure_spherical.rst Normal file
View File

@ -0,0 +1,61 @@
.. index:: compute pressure/spherical
compute pressure/spherical command
==================================
Syntax
""""""
.. parsed-literal::
compute ID group-ID pressure/spherical x0 y0 z0 bin_width Rmax
* ID, group-ID are documented in :doc:`compute <compute>` command
* pressure/cartesian = style name of this compute command
* x0, y0, z0 = origin of the spherical coordinate system
* bin_width = width of spherical shells
* Rmax = maximum radius of spherical shells
Examples
""""""""
.. code-block:: LAMMPS
compute 1 all pressure/spherical 0 0 0 0.1 10
Description
"""""""""""
This computes the diagonal components of the spherical pressure tensor in spherical shells, as described in :ref:`(Ikeshoji)<Ikeshoji2>`. This can be used to calculate the local pressure tensor components of for example droplets and bubbles. This compute obeys momentum balance. This compute uses the Irving-Kirkwood contour, whici is the straight line between particle pairs. The pressure tensor is split into a kinetic contribution :math:`P^k` and a virial contribution :math:`P^v`. The sum gives the pressure :math:`P = P^k+P^v`.
Output info
"""""""""""
This compute outputs a global array with 8 columns and Rmax/bin_width. The output columns are position of the center of the local spherical shell, number density, :math:`P^k_{rr}`, :math:`P^k_{\theta\theta}`, :math:`P^k_{\phi\phi}`, :math:`P^v_{rr}`, :math:`P^v_{\theta\theta}`, and :math:`P^v_{\phi\phi}`.
This array can be output with :doc:`fix ave/time <fix_ave_time>`,
.. code-block:: LAMMPS
compute p all pressure/spherical 0 0 0 0.1 10
fix 2 all ave/time 100 1 100 c_p[*] file dump_p.out mode vector
Restrictions
""""""""""""
This compute is part of the EXTRA-COMPUTE package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
Related commands
""""""""""""""""
:doc:`compute stress/atom <compute_stress_atom>`, :doc:`compute pressure/cylinder <compute_pressure_cylinder>`, :doc:`compute pressure/cartesian <compute_pressure_cartesian>`
Default
"""""""
none
----------
.. _Ikeshoji2:
**(Ikeshoji)** Ikeshoji, Hafskjold, Furuholt, Mol Sim, 29, 101-109, (2003).

537
src/EXTRA-COMPUTE/compute_pressure_cartesian.cpp Normal file
View File

@ -0,0 +1,537 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS dir1ectory.
------------------------------------------------------------------------- */
#include "compute_pressure_cartesian.h"
#include <mpi.h>
#include <cmath>
#include <cstring>
#include "atom.h"
#include "force.h"
#include "pair.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "citeme.h"
#include "domain.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
static const char cite_compute_pressure_cartesian[] =
"compute pressure/cartesian:\n\n"
"@article{ikeshoji2003molecular,\n"
"title={Molecular-level calculation scheme for pressure in inhomogeneous systems of flat and spherical layers},\n"
"author={Ikeshoji, Tamio and Hafskjold, Bjørn and Furuholt, Hilde},\n"
"journal={Molecular Simulation},\n"
"volume={29},\n"
"number={2},\n"
"pages={101--109},\n"
"year={2003},\n"
"publisher={Taylor & Francis}\n"
"}\n\n";
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
ComputePressureCart::ComputePressureCart(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
dens(NULL),
pkxx(NULL), pkyy(NULL), pkzz(NULL),
pcxx(NULL), pcyy(NULL), pczz(NULL),
tdens(NULL),
tpkxx(NULL), tpkyy(NULL), tpkzz(NULL),
tpcxx(NULL), tpcyy(NULL), tpczz(NULL)
{
if (lmp->citeme)
lmp->citeme->add(cite_compute_pressure_cartesian);
// narg == 5 for one-dimensional and narg == 7 for two-dimensional
if (narg == 5) dims = 1;
else if (narg == 7) dims = 2;
else error->all(FLERR,"Illegal compute pressure/cartesian command. Illegal number of arguments.");
if (strcmp(arg[3], "x") == 0) dir1 = 0;
else if (strcmp(arg[3], "y") == 0) dir1 = 1;
else if (strcmp(arg[3], "z") == 0) dir1 = 2;
else error->all(FLERR,"Illegal compute pressure/cartesian command.");
dir2 = 0;
bin_width1 = utils::numeric(FLERR, arg[4], false, lmp);
bin_width2 = 0.0;
nbins1 = (int)((domain->boxhi[dir1]-domain->boxlo[dir1])/bin_width1);
nbins2 = 1;
invV = bin_width1;
if (bin_width1 <= 0.0)
error->all(FLERR,"Illegal compute pressure/cartesian command. Bin width must be > 0");
else if (bin_width1 > domain->boxhi[dir1]-domain->boxlo[dir1])
error->all(FLERR,"Illegal compute pressure/cartesian command. Bin width must be smaller than box size.");
if (dims == 2){
if (strcmp(arg[5], "x") == 0) dir2 = 0;
else if (strcmp(arg[5], "y") == 0) dir2 = 1;
else if (strcmp(arg[5], "z") == 0) dir2 = 2;
else error->all(FLERR,"Illegal compute pressure/cartesian command.");
bin_width2 = utils::numeric(FLERR, arg[6], false, lmp);
nbins2 = (int)((domain->boxhi[dir2]-domain->boxlo[dir2])/bin_width2);
invV *= bin_width2;
if (bin_width2 <= 0.0)
error->all(FLERR,"Illegal compute pressure/cartesian command. Bin width must be > 0");
else if (bin_width2 > domain->boxhi[dir2]-domain->boxlo[dir2])
error->all(FLERR,"Illegal compute pressure/cartesian command. Bin width must be smaller than box size.");
}
for (int i = 0; i < 3; i++)
if ((dims == 1 && i != dir1) || (dims == 2 && (i != dir1 && i != dir2)))
invV *= domain->boxhi[i] - domain->boxlo[i];
invV = 1.0/invV;
array_flag = 1;
vector_flag = 0;
extarray = 0;
size_array_cols = 7+dims; // dir1, dir2, number density, pkxx, pkyy, pkzz, pcxx, pcyy, pczz
size_array_rows = nbins1*nbins2;
memory->create(dens, nbins1*nbins2, "dens" );
memory->create(pkxx, nbins1*nbins2, "pkxx" );
memory->create(pkyy, nbins1*nbins2, "pkyy" );
memory->create(pkzz, nbins1*nbins2, "pkzz" );
memory->create(pcxx, nbins1*nbins2, "pcxx" );
memory->create(pcyy, nbins1*nbins2, "pcyy" );
memory->create(pczz, nbins1*nbins2, "pczz" );
memory->create(tdens, nbins1*nbins2, "tdens");
memory->create(tpkxx, nbins1*nbins2, "tpkxx");
memory->create(tpkyy, nbins1*nbins2, "tpkyy");
memory->create(tpkzz, nbins1*nbins2, "tpkzz");
memory->create(tpcxx, nbins1*nbins2, "tpcxx");
memory->create(tpcyy, nbins1*nbins2, "tpcyy");
memory->create(tpczz, nbins1*nbins2, "tpczz");
memory->create(array, size_array_rows, size_array_cols, "pressure:cartesian:output");
}
/* ---------------------------------------------------------------------- */
ComputePressureCart::~ComputePressureCart()
{
memory->destroy(dens);
memory->destroy(pkxx);
memory->destroy(pkyy);
memory->destroy(pkzz);
memory->destroy(pcxx);
memory->destroy(pcyy);
memory->destroy(pczz);
memory->destroy(tdens);
memory->destroy(tpkxx);
memory->destroy(tpkyy);
memory->destroy(tpkzz);
memory->destroy(tpcxx);
memory->destroy(tpcyy);
memory->destroy(tpczz);
memory->destroy(array);
}
/* ---------------------------------------------------------------------- */
void ComputePressureCart::init()
{
if (force->pair == NULL)
error->all(FLERR,"No pair style is defined for compute pressure/cartesian");
if (force->pair->single_enable == 0)
error->all(FLERR,"Pair style does not support compute pressure/cartesian");
// need an occasional half neighbor list.
int irequest = neighbor->request(this, instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
neighbor->requests[irequest]->occasional = 1;
}
/* ---------------------------------------------------------------------- */
void ComputePressureCart::init_list(int /* id */, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
count pairs and compute pair info on this proc
only count pair once if newton_pair is off
both atom I,J must be in group
if flag is set, compute requested info about pair
------------------------------------------------------------------------- */
void ComputePressureCart::compute_array()
{
//invoked_array = update->ntimestep;
int me;
MPI_Comm_rank(world, &me);
int i, j, ii, jj, inum, jnum, itype, jtype;
int bin, bin1, bin2, bin3;
tagint itag, jtag;
double xtmp, ytmp, ztmp, delx, dely, delz;
double rsq, fpair, factor_coul, factor_lj;
int *ilist, *jlist, *numneigh, **firstneigh;
double **x = atom->x;
double **v = atom->v;
double *mass = atom->mass;
tagint *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
// invoke half neighbor list (will copy or build if necessary)
neighbor->build_one(list);
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// Zero arrays
for (bin = 0; bin < nbins1*nbins2; bin++) {
tdens[bin] = 0;
tpkxx[bin] = 0;
tpkyy[bin] = 0;
tpkzz[bin] = 0;
tpcxx[bin] = 0;
tpcyy[bin] = 0;
tpczz[bin] = 0;
}
// calculate number density and kinetic contribution to pressure
for (i = 0; i < nlocal; i++){
bin1 = (int)(x[i][dir1]/bin_width1)%nbins1;
bin2 = 0;
if (dims == 2)
bin2 = (int)(x[i][dir2]/bin_width2)%nbins2;
j = bin1+bin2*nbins1;
tdens[j] += 1;
// TODO: Subtract streaming velocity
tpkxx[j] += mass[type[i]]*v[i][0]*v[i][0];
tpkyy[j] += mass[type[i]]*v[i][1]*v[i][1];
tpkzz[j] += mass[type[i]]*v[i][2]*v[i][2];
}
// loop over neighbors of my atoms
// skip if I or J are not in group
// for newton = 0 and J = ghost atom,
// need to insure I,J pair is only output by one proc
// use same itag,jtag logic as in Neighbor::neigh_half_nsq()
// for flag = 0, just count pair interactions within force cutoff
// for flag = 1, calculate requested output fields
Pair *pair = force->pair;
double **cutsq = force->pair->cutsq;
double risq, rjsq;
double xi1, xi2, xj1, xj2;
for (ii = 0; ii < inum; ii++){
i = ilist[ii];
if (!(mask[i] & groupbit)) continue;
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
xi1 = x[i][dir1];
xi2 = x[i][dir2];
itag = tag[i];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
if (!(mask[j] & groupbit)) continue;
// itag = jtag is possible for long cutoffs that include images of self
// do calculation only on appropriate processor
if (newton_pair == 0 && j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
}
xj1 = x[j][dir1];
xj2 = x[j][dir2];
delx = x[j][0]-xtmp;
dely = x[j][1]-ytmp;
delz = x[j][2]-ztmp;
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
// Check if inside cut-off
if (rsq >= cutsq[itype][jtype]) continue;
pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
if (dims == 1)
compute_pressure_1d(fpair, xi1, xj1, delx, dely, delz);
if (dims == 2)
compute_pressure_2d(fpair, xi1, xi2, xj1, xj2, delx, dely, delz);
}
}
// normalize pressure
for (bin = 0; bin < nbins1*nbins2; bin++) {
tdens[bin] *= invV;
tpkxx[bin] *= invV;
tpkyy[bin] *= invV;
tpkzz[bin] *= invV;
tpcxx[bin] *= invV;
tpcyy[bin] *= invV;
tpczz[bin] *= invV;
}
// communicate across processors
MPI_Allreduce(tdens, dens, nbins1*nbins2, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpkxx, pkxx, nbins1*nbins2, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpkyy, pkyy, nbins1*nbins2, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpkzz, pkzz, nbins1*nbins2, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpcxx, pcxx, nbins1*nbins2, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpcyy, pcyy, nbins1*nbins2, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpczz, pczz, nbins1*nbins2, MPI_DOUBLE, MPI_SUM, world);
// populate array to output.
for (bin = 0; bin < nbins1*nbins2; bin++) {
array[bin][0] = (bin%nbins1+0.5)*bin_width1;
if (dims == 2)
array[bin][1] = ((int)(bin/nbins1)+0.5)*bin_width2;
array[bin][0+dims] = dens[bin];
array[bin][1+dims] = pkxx[bin];
array[bin][2+dims] = pkyy[bin];
array[bin][3+dims] = pkzz[bin];
array[bin][4+dims] = pcxx[bin];
array[bin][5+dims] = pcyy[bin];
array[bin][6+dims] = pczz[bin];
}
}
void ComputePressureCart::compute_pressure_1d(double fpair, double xi, double xj, double delx, double dely, double delz)
{
int bin_s, bin_e, bin_step, bin, bin_limit;
double xa, xb, l_sum;
if (xi < domain->boxlo[dir1])
xi += (domain->boxhi[dir1]-domain->boxlo[dir1]);
else if (xi > domain->boxhi[dir1])
xi -= (domain->boxhi[dir1]-domain->boxlo[dir1]);
if (xj < domain->boxlo[dir1])
xj += (domain->boxhi[dir1]-domain->boxlo[dir1]);
else if (xj > domain->boxhi[dir1])
xj -= (domain->boxhi[dir1]-domain->boxlo[dir1]);
// Integrating contour from bin_s to bin_e
bin_s = lround((int)((xi - domain->boxlo[dir1])/bin_width1)%nbins1);
bin_e = lround((int)((xj - domain->boxlo[dir1])/bin_width1)%nbins1);
// If not periodic in dir1
if (domain->periodicity[dir1] == 0){
bin_s = lround((int)((xi - domain->boxlo[dir1])/bin_width1));
bin_e = lround((int)((xj - domain->boxlo[dir1])/bin_width1));
if (bin_e == nbins1)
bin_e--;
if (bin_s == nbins1)
bin_s--;
}
bin_step = 1;
if (domain->periodicity[dir1] == 1){
if (bin_e - bin_s > 0.5*nbins1)
bin_step = -1;
else if (bin_s - bin_e > 0.5*nbins1)
bin_step = 1;
else if (bin_s > bin_e)
bin_step = -1;
}
else {
if (bin_s > bin_e)
bin_step = -1;
}
if (domain->periodicity[dir1] == 1)
bin_limit = (bin_e+bin_step)%nbins1 < 0 ? (bin_e+bin_step)%nbins1+nbins1 : (bin_e+bin_step)%nbins1;
else
bin_limit = bin_e+bin_step;
bin = bin_s;
// Integrate from bin_s to bin_e with step bin_step.
while (bin != bin_limit){
// Calculating exit and entry point (xa, xb). Checking if inside current bin.
if (bin == bin_s){
if (domain->periodicity[dir1] == 1)
xa = fmod(xi, domain->boxhi[dir1])+domain->boxlo[dir1];
else
xa = xi;
}
else
xa = (bin_step == 1) ? bin*bin_width1 : (bin+1)*bin_width1;
if (bin == bin_e){
if (domain->periodicity[dir1] == 1)
xb = fmod(xj, domain->boxhi[dir1])+domain->boxlo[dir1];
else
xb = xj;
}
else
xb = (bin_step == 1) ? (bin+1)*bin_width1 : bin*bin_width1;
if (bin < 0 || bin >= nbins1)
error->all(FLERR,"ERROR: Bin outside simulation.");
if (bin_s != bin_e){
if (dir1 == 0){
tpcxx[bin] += (fpair*delx*delx) * (xb - xa)/delx;
tpcyy[bin] += (fpair*dely*dely) * (xb - xa)/delx;
tpczz[bin] += (fpair*delz*delz) * (xb - xa)/delx;
}
else if (dir1 == 1){
tpcxx[bin] += (fpair*delx*delx) * (xb - xa)/dely;
tpcyy[bin] += (fpair*dely*dely) * (xb - xa)/dely;
tpczz[bin] += (fpair*delz*delz) * (xb - xa)/dely;
}
else if (dir1 == 2){
tpcxx[bin] += (fpair*delx*delx) * (xb - xa)/delz;
tpcyy[bin] += (fpair*dely*dely) * (xb - xa)/delz;
tpczz[bin] += (fpair*delz*delz) * (xb - xa)/delz;
}
}
// Particle i and j in same bin. Avoiding zero divided by zero.
else {
tpcxx[bin] += fpair*delx*delx;
tpcyy[bin] += fpair*dely*dely;
tpczz[bin] += fpair*delz*delz;
}
// Stepping bin to next bin
if (domain->periodicity[dir1] == 1)
bin = (bin+bin_step)%nbins1 < 0 ? (bin+bin_step)%nbins1+nbins1 : (bin+bin_step)%nbins1;
else
bin = bin+bin_step;
}
}
void ComputePressureCart::compute_pressure_2d(double fpair, double xi, double yi, double xj, double yj, double delx, double dely, double delz)
{
int bin1, bin2, next_bin1, next_bin2;
double la = 0.0, lb = 0.0, l_sum = 0.0;
double rij[3] = {delx, dely, delz};
double l1 = 0.0, l2, rij1, rij2;
double SMALL = 10e-5;
rij1 = rij[dir1];
rij2 = rij[dir2];
next_bin1 = (int)floor(xi/bin_width1);
next_bin2 = (int)floor(yi/bin_width2);
// Integrating along line
while (lb != 1.0){
bin1 = next_bin1;
bin2 = next_bin2;
if (rij1 > 0)
l1 = ((bin1+1)*bin_width1 - xi)/rij1;
else
l1 = (bin1*bin_width1 - xi)/rij1;
if (rij2 > 0)
l2 = ((bin2+1)*bin_width2 - yi)/rij2;
else
l2 = (bin2*bin_width2 - yi)/rij2;
if ((l1 < l2 || l2 < lb+SMALL) && l1 <= 1.0 && l1 > lb){
lb = l1;
next_bin1 = bin1+(int)(rij1/fabs(rij1));
}
else if (l2 <= 1.0 && l2 > lb){
lb = l2;
next_bin2 = bin2+(int)(rij2/fabs(rij2));
}
else
lb = 1.0;
// Periodic boundary conditions
if (domain->periodicity[dir1] == 1){
if (bin1 < 0)
bin1 += nbins1;
else if (bin1 >= nbins1)
bin1 -= nbins1;
}
else if (bin1 < 0)
bin1 = 0;
else if (bin1 >= nbins1)
bin1 = nbins1-1;
if (domain->periodicity[dir2] == 1){
if (bin2 < 0)
bin2 += nbins2;
else if (bin2 >= nbins2)
bin2 -= nbins2;
}
else if (bin2 < 0)
bin2 = 0;
else if (bin2 >= nbins2)
bin2 = nbins2-1;
if (bin1+bin2*nbins1 >= nbins1*nbins2)
error->all(FLERR, "Bin outside");
tpcxx[bin1+bin2*nbins1] += (fpair*delx*delx*(lb-la));
tpcyy[bin1+bin2*nbins1] += (fpair*dely*dely*(lb-la));
tpczz[bin1+bin2*nbins1] += (fpair*delz*delz*(lb-la));
l_sum += lb-la;
la = lb;
}
if (l_sum > 1.0+SMALL || l_sum < 1.0-SMALL)
error->all(FLERR,"Sum of fractional line segments does not equal 1.");
}
/* ----------------------------------------------------------------------
memory usage of data
------------------------------------------------------------------------- */
double ComputePressureCart::memory_usage()
{
return (14.0+dims+7)*(double)(nbins1*nbins2)*sizeof(double);
}

71
src/EXTRA-COMPUTE/compute_pressure_cartesian.h Normal file
View File

@ -0,0 +1,71 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(pressure/cartesian,ComputePressureCart)
#else
#ifndef LMP_COMPUTE_PRESSURE_CARTESIAN
#define LMP_COMPUTE_PRESSURE_CARTESIAN
#include "compute.h"
namespace LAMMPS_NS {
class ComputePressureCart : public Compute {
public:
ComputePressureCart(class LAMMPS *, int, char **);
~ComputePressureCart();
void init();
void init_list(int, class NeighList *);
void compute_array();
double memory_usage();
private:
int nbins1, nbins2, dir1, dir2, dims;
double bin_width1, bin_width2, invV;
// Number density, kinetic and configurational contribution to the pressure.
double *dens, *pkxx, *pkyy, *pkzz, *pcxx, *pcyy, *pczz;
double *tdens, *tpkxx, *tpkyy, *tpkzz, *tpcxx, *tpcyy, *tpczz;
class NeighList *list;
void compute_pressure_1d(double, double, double, double, double, double);
void compute_pressure_2d(double, double, double, double, double, double, double, double);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: No pair style is defined for compute pressure/cartesian
Self-explanatory.
E: Pair style does not support compute pressure/cartesian
The pair style does not have a single() function, so it can
not be invoked by compute pressure/cartesian.
*/

413
src/EXTRA-COMPUTE/compute_pressure_spherical.cpp Normal file
View File

@ -0,0 +1,413 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS dir1ectory.
------------------------------------------------------------------------- */
#include "compute_pressure_spherical.h"
#include <mpi.h>
#include <cmath>
#include <cstring>
#include "atom.h"
#include "update.h"
#include "force.h"
#include "pair.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "citeme.h"
#include "domain.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
static const char cite_compute_pressure_sphere[] =
"compute pressure/spherical:\n\n"
"@article{ikeshoji2003molecular,\n"
"title={Molecular-level calculation scheme for pressure in inhomogeneous systems of flat and spherical layers},\n"
"author={Ikeshoji, Tamio and Hafskjold, Bjørn and Furuholt, Hilde},\n"
"journal={Molecular Simulation},\n"
"volume={29},\n"
"number={2},\n"
"pages={101--109},\n"
"year={2003},\n"
"publisher={Taylor & Francis}\n"
"}\n\n";
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
ComputePressureSpherical::ComputePressureSpherical(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
dens(nullptr),
pkrr(nullptr), pktt(nullptr), pkpp(nullptr),
pcrr(nullptr), pctt(nullptr), pcpp(nullptr),
tdens(nullptr),
tpkrr(nullptr), tpktt(nullptr), tpkpp(nullptr),
tpcrr(nullptr), tpctt(nullptr), tpcpp(nullptr)
{
if (lmp->citeme)
lmp->citeme->add(cite_compute_pressure_sphere);
if (narg != 8)
error->all(FLERR,"Illegal compute pressure/spherical command. Illegal number of arguments.");
x0 = utils::numeric(FLERR, arg[3], false, lmp);
y0 = utils::numeric(FLERR, arg[4], false, lmp);
z0 = utils::numeric(FLERR, arg[5], false, lmp);
bin_width = utils::numeric(FLERR, arg[6], false, lmp);
Rmax = utils::numeric(FLERR, arg[7], false, lmp);
nbins = (int)(Rmax/bin_width)+1;
if (bin_width <= 0.0)
error->all(FLERR,"Illegal compute pressure/spherical command. Bin width must be > 0");
array_flag = 1;
vector_flag = 0;
extarray = 0;
size_array_cols = 8; // r, dens, pkrr, pktt, pkpp, pcrr, pctt, pcpp
size_array_rows = nbins;
memory->create(invV, nbins, "compute/pressure/spherical:invV");
memory->create(dens, nbins, "compute/pressure/spherical:dens");
memory->create(pkrr, nbins, "compute/pressure/spherical:pkrr");
memory->create(pktt, nbins, "compute/pressure/spherical:pktt");
memory->create(pkpp, nbins, "compute/pressure/spherical:pkpp");
memory->create(pcrr, nbins, "compute/pressure/spherical:pcrr");
memory->create(pctt, nbins, "compute/pressure/spherical:pctt");
memory->create(pcpp, nbins, "compute/pressure/spherical:pcpp");
memory->create(tdens, nbins, "compute/pressure/spherical:tdens");
memory->create(tpkrr, nbins, "compute/pressure/spherical:tpkrr");
memory->create(tpktt, nbins, "compute/pressure/spherical:tpktt");
memory->create(tpkpp, nbins, "compute/pressure/spherical:tpkpp");
memory->create(tpcrr, nbins, "compute/pressure/spherical:tpcrr");
memory->create(tpctt, nbins, "compute/pressure/spherical:tpctt");
memory->create(tpcpp, nbins, "compute/pressure/spherical:tpcpp");
memory->create(array, size_array_rows, size_array_cols, "compute/pressure/spherical:array");
}
/* ---------------------------------------------------------------------- */
ComputePressureSpherical::~ComputePressureSpherical()
{
memory->destroy(invV);
memory->destroy(dens);
memory->destroy(pkrr);
memory->destroy(pktt);
memory->destroy(pkpp);
memory->destroy(pcrr);
memory->destroy(pctt);
memory->destroy(pcpp);
memory->destroy(tdens);
memory->destroy(tpkrr);
memory->destroy(tpktt);
memory->destroy(tpkpp);
memory->destroy(tpcrr);
memory->destroy(tpctt);
memory->destroy(tpcpp);
memory->destroy(array);
}
/* ---------------------------------------------------------------------- */
void ComputePressureSpherical::init()
{
if (force->pair == nullptr)
error->all(FLERR,"No pair style is defined for compute pressure/spherical");
if (force->pair->single_enable == 0)
error->all(FLERR,"Pair style does not support compute pressure/spherical");
for (int bin = 0; bin < nbins; bin++)
invV[bin] = 3.0 / (4.0*MY_PI*(pow((bin+1)*bin_width, 3.0) - pow(bin*bin_width, 3.0)));
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
neighbor->requests[irequest]->occasional = 1;
}
/* ---------------------------------------------------------------------- */
void ComputePressureSpherical::init_list(int /* id */, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
count pairs and compute pair info on this proc
only count pair once if newton_pair is off
both atom I,J must be in group
if flag is set, compute requested info about pair
------------------------------------------------------------------------- */
void ComputePressureSpherical::compute_array()
{
invoked_array = update->ntimestep;
int me;
MPI_Comm_rank(world, &me);
int bin;
// Zero arrays
for (int bin = 0; bin < nbins; bin++) {
tdens[bin] = 0.0;
tpkrr[bin] = 0.0;
tpktt[bin] = 0.0;
tpkpp[bin] = 0.0;
tpcrr[bin] = 0.0;
tpctt[bin] = 0.0;
tpcpp[bin] = 0.0;
}
int i, j, ii, jj, inum, jnum, itype, jtype;
tagint itag, jtag;
double ri[3], xtmp, ytmp, ztmp, tmp;
double rsq, fpair, factor_coul, factor_lj;
int *ilist, *jlist, *numneigh, **firstneigh;
double r, vr, vt, vp, theta;
double **x = atom->x;
double **v = atom->v;
double *mass = atom->mass;
tagint *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
// invoke half neighbor list (will copy or build if necessary)
neighbor->build_one(list);
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// calculate number density and kinetic contribution to pressure
for (i = 0; i < nlocal; i++){
ri[0] = x[i][0]-x0;
ri[1] = x[i][1]-y0;
ri[2] = x[i][2]-z0;
for (j = 0; j < 3; j++){
tmp = domain->boxhi[j]-domain->boxlo[j];
if (ri[j] > 0.5*tmp)
ri[j] -= tmp;
else if (ri[j] < -0.5*tmp)
ri[j] += tmp;
}
r = sqrt(ri[0]*ri[0]+ri[1]*ri[1]+ri[2]*ri[2]);
if (r >= Rmax) continue;
bin = (int)(r/bin_width);
// Avoiding division by zero
if (r != 0){
theta = acos(ri[2]/r);
tdens[bin] += 1.0;
vr = (ri[0]*v[i][0]+ri[1]*v[i][1]+ri[2]*v[i][2]) / r;
vt = r*sin(theta)/(pow(ri[1]/ri[0], 2.0) + 1.0)*((ri[0]*v[i][1] - ri[1]*v[i][0])/(ri[0]*ri[0]));
vp = (ri[2]*vr - r*v[i][2])/(r*sqrt(1.0 - pow(ri[2]/r, 2.0)));
tpkrr[bin] += vr*vr;
tpktt[bin] += vt*vt;
tpkpp[bin] += vp*vp;
}
}
// loop over neighbors of my atoms
// skip if I or J are not in group
// for newton = 0 and J = ghost atom,
// need to insure I,J pair is only output by one proc
// use same itag,jtag logic as in Neighbor::neigh_half_nsq()
// for flag = 0, just count pair interactions within force cutoff
// for flag = 1, calculate requested output fields
Pair *pair = force->pair;
double **cutsq = force->pair->cutsq;
double risq, rjsq;
double dir1i, dir2i, dir3i, dir1j, dir2j, dir3j;
double xi, yi, zi, xj, yj, zj;
double qi[3], l1, l2, l3, l4, R1, R2, Fa, Fb, l_sum;
double rij, f, ririj, sqr, la, lb, sql0, lambda0;
double SMALL = 10e-12;
double rsqxy, ririjxy, sqrixy, sqlxy0, A, B, C;
int end_bin;
for (ii = 0; ii < inum; ii++){
i = ilist[ii];
if (!(mask[i] & groupbit)) continue;
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itag = tag[i];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
if (!(mask[j] & groupbit)) continue;
// itag = jtag is possible for long cutoffs that include images of self
// do calculation only on appropriate processor
if (newton_pair == 0 && j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
}
ri[0] = x[j][0]-xtmp;
ri[1] = x[j][1]-ytmp;
ri[2] = x[j][2]-ztmp;
rsq = ri[0]*ri[0]+ri[1]*ri[1]+ri[2]*ri[2];
jtype = type[j];
// Check if inside cut-off
if (rsq >= cutsq[itype][jtype]) continue;
pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
qi[0] = xtmp-x0;
qi[1] = ytmp-y0;
qi[2] = ztmp-z0;
for (int k = 0; k < 3; k++){
tmp = domain->boxhi[k]-domain->boxlo[k];
if (qi[k] > 0.5*tmp)
qi[k] -= tmp;
else if (qi[k] < -0.5*tmp)
qi[k] += tmp;
}
l_sum = 0.0;
rij = sqrt(rsq);
f = -rij*fpair;
ririj = qi[0]*ri[0] + qi[1]*ri[1] + qi[2]*ri[2];
sqr = qi[0]*qi[0] + qi[1]*qi[1] + qi[2]*qi[2];
la = 0.0;
lb = 0.0;
sql0 = sqr - ririj*ririj/rsq;
lambda0 = -ririj/rsq;
rsqxy = ri[0]*ri[0]+ri[1]*ri[1];
ririjxy = qi[0]*ri[0]+qi[1]*ri[1];
sqrixy = qi[0]*qi[0]+qi[1]*qi[1];
sqlxy0 = sqrixy - ririjxy*ririjxy / rsqxy;
A = pow(qi[0]*ri[1] - qi[1]*ri[0], 2.0);
C = sqrt(rsqxy*sqrixy - ririjxy*ririjxy);
B = sqrt(rsq*sqr - ririj*ririj);
end_bin = (int)(sqrt(rsq+2.0*ririj+sqr)/bin_width);
while (lb < 1.0) {
l1 = lb+SMALL;
bin = (int)floor(sqrt(rsq*l1*l1+2.0*ririj*l1+sqr)/bin_width);
R1 = (bin)*bin_width;
R2 = (bin+1)*bin_width;
l1 = lambda0 + sqrt((R1*R1 - sql0)/rsq);
l2 = lambda0 - sqrt((R1*R1 - sql0)/rsq);
l3 = lambda0 + sqrt((R2*R2 - sql0)/rsq);
l4 = lambda0 - sqrt((R2*R2 - sql0)/rsq);
if (l4 >= 0.0 && l4 <= 1.0 && l4 > la)
lb = l4;
else if (l3 >= 0.0 && l3 <= 1.0 && l3 > la)
lb = l3;
else if (l2 >= 0.0 && l2 <= 1.0 && l2 > la)
lb = l2;
else if (l1 >= 0.0 && l1 <= 1.0 && l1 > la)
lb = l1;
else
lb = 1.0;
if (bin == end_bin)
lb = 1.0;
if (la > lb)
error->all(FLERR,"Error: la > lb\n");
if (bin >= 0 && bin < nbins){
if (sql0 > SMALL){
Fa = -B*atan2(rsq*la+ririj, B);
Fb = -B*atan2(rsq*lb+ririj, B);
tpcrr[bin] -= (f/rij)*( rsq*(lb - la) + Fb - Fa );
tpcpp[bin] -= (f/rij)*( Fa - Fb )/2.0;
}
else
tpcrr[bin] -= f*rij*(lb-la);
if (sqlxy0 > SMALL){
tpctt[bin] -= (f/rij)*A/C*(atan2(rsqxy*lb + ririjxy, C) - atan2(rsqxy*la + ririjxy, C));
}
}
l_sum += lb-la;
la = lb;
}
// Error check
if (fabs(l_sum - 1.0) > SMALL)
error->all(FLERR,"ERROR: The sum of the fractional line segments, calculated in spherical pressure tensor, does not sum up to 1.");
}
}
// normalize pressure
for (bin = 0; bin < nbins; bin++) {
tdens[bin] *= invV[bin];
tpkrr[bin] *= invV[bin];
tpktt[bin] *= invV[bin];
tpkpp[bin] *= invV[bin];
tpcrr[bin] *= invV[bin];
tpctt[bin] *= invV[bin];
tpcpp[bin] *= invV[bin];
}
// communicate across processors
MPI_Allreduce(tdens, dens, nbins, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpkrr, pkrr, nbins, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpktt, pktt, nbins, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpkpp, pkpp, nbins, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpcrr, pcrr, nbins, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpctt, pctt, nbins, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(tpcpp, pcpp, nbins, MPI_DOUBLE, MPI_SUM, world);
// populate array to output.
for (bin = 0; bin < nbins; bin++) {
array[bin][0] = (bin+0.5)*bin_width;
array[bin][1] = dens[bin];
array[bin][2] = pkrr[bin];
array[bin][3] = pktt[bin];
array[bin][4] = pkpp[bin];
array[bin][5] = pcrr[bin];
array[bin][6] = pctt[bin];
array[bin][7] = pcpp[bin];
}
}
double ComputePressureSpherical::memory_usage()
{
return 15.0*(double)(nbins+size_array_rows*size_array_cols)*sizeof(double);
}

69
src/EXTRA-COMPUTE/compute_pressure_spherical.h Normal file
View File

@ -0,0 +1,69 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(pressure/spherical,ComputePressureSpherical)
#else
#ifndef LMP_COMPUTE_PRESSURE_SPHERICAL
#define LMP_COMPUTE_PRESSURE_SPHERICAL
#include "compute.h"
namespace LAMMPS_NS {
class ComputePressureSpherical : public Compute {
public:
ComputePressureSpherical(class LAMMPS *, int, char **);
~ComputePressureSpherical();
void init();
void init_list(int, class NeighList *);
void compute_array();
double memory_usage();
private:
int nbins;
double bin_width, x0, y0, z0, Rmax;
// Number density, kinetic and configurational contribution to the pressure.
double *invV, *dens, *pkrr, *pktt, *pkpp, *pcrr, *pctt, *pcpp;
double *tdens, *tpkrr, *tpktt, *tpkpp, *tpcrr, *tpctt, *tpcpp;
class NeighList *list;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: No pair style is defined for compute pressure/cartesian
Self-explanatory.
E: Pair style does not support compute pressure/cartesian
The pair style does not have a single() function, so it can
not be invoked by compute pressure/cartesian.
*/