From 6a6d08e18ef9f47e679dd79b841bdc1faf05d849 Mon Sep 17 00:00:00 2001 From: Julien Devemy Date: Tue, 25 Jun 2019 12:01:29 +0200 Subject: [PATCH] Better compute_pressure hybrid and doc --- doc/src/compute_pressure.txt | 8 ++++-- src/compute_pressure.cpp | 50 ++++++++++++++++++++++++++---------- 2 files changed, 42 insertions(+), 16 deletions(-) diff --git a/doc/src/compute_pressure.txt b/doc/src/compute_pressure.txt index bd6e38e392..fd9bfb7aff 100644 --- a/doc/src/compute_pressure.txt +++ b/doc/src/compute_pressure.txt @@ -16,12 +16,13 @@ ID, group-ID are documented in "compute"_compute.html command pressure = style name of this compute command temp-ID = ID of compute that calculates temperature, can be NULL if not needed zero or more keywords may be appended -keyword = {ke} or {pair} or {bond} or {angle} or {dihedral} or {improper} or {kspace} or {fix} or {virial} :ul +keyword = {ke} or {pair} or {bond} or {angle} or {dihedral} or {improper} or {kspace} or {fix} or {virial} or {hybridpair} :ul [Examples:] compute 1 all pressure thermo_temp -compute 1 all pressure NULL pair bond :pre +compute 1 all pressure NULL pair bond +compute 1 all pressure NULL hybridpair lj/cut :pre [Description:] @@ -67,6 +68,9 @@ extra keywords are listed, then only those components are summed to compute temperature or ke and/or the virial. The {virial} keyword means include all terms except the kinetic energy {ke}. +The {hybridpair} keyword means to only include contribution +from a subpair in a {hybrid} or {hybrid/overlay} pair style. + Details of how LAMMPS computes the virial efficiently for the entire system, including for many-body potentials and accounting for the effects of periodic boundary conditions are discussed in diff --git a/src/compute_pressure.cpp b/src/compute_pressure.cpp index 8381e04fa8..96f577219d 100644 --- a/src/compute_pressure.cpp +++ b/src/compute_pressure.cpp @@ -81,8 +81,37 @@ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) : int iarg = 4; while (iarg < narg) { if (strcmp(arg[iarg],"ke") == 0) keflag = 1; - else if (strcmp(arg[iarg],"hybridpair") == 0) - hybridpairflag = force->inumeric(FLERR, arg[++iarg]); + else if (strcmp(arg[iarg],"hybridpair") == 0) { + int n = strlen(arg[++iarg]) + 1; + if (lmp->suffix) n += strlen(lmp->suffix) + 1; + pstyle = new char[n]; + strcpy(pstyle,arg[iarg++]); + + nsub = 0; + + if (narg > iarg) { + if (isdigit(arg[iarg][0])) { + nsub = force->inumeric(FLERR,arg[iarg]); + ++iarg; + if (nsub <= 0) + error->all(FLERR,"Illegal compute pressure command"); + } + } + + // check if pair style with and without suffix exists + + hybridpair = (Pair *) force->pair_match(pstyle,1,nsub); + if (!hybridpair && lmp->suffix) { + strcat(pstyle,"/"); + strcat(pstyle,lmp->suffix); + hybridpair = (Pair *) force->pair_match(pstyle,1,nsub); + } + + if (!hybridpair) + error->all(FLERR,"Unrecognized pair style in compute pressure command"); + + hybridpairflag = 1; + } else if (strcmp(arg[iarg],"pair") == 0) pairflag = 1; else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1; else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1; @@ -144,12 +173,7 @@ void ComputePressure::init() nvirial = 0; vptr = NULL; - if (hybridpairflag > 0 && force->pair) { - if (strstr(force->pair_style, "hybrid")) { - PairHybrid *ph = (PairHybrid *) force->pair; - if (hybridpairflag <= ph->nstyles) nvirial++; - } - } + if (hybridpairflag && force->pair) nvirial++; if (pairflag && force->pair) nvirial++; if (bondflag && atom->molecular && force->bond) nvirial++; if (angleflag && atom->molecular && force->angle) nvirial++; @@ -162,12 +186,10 @@ void ComputePressure::init() if (nvirial) { vptr = new double*[nvirial]; nvirial = 0; - if (hybridpairflag > 0 && force->pair) { - if (strstr(force->pair_style, "hybrid")) { - PairHybrid *ph = (PairHybrid *) force->pair; - if (hybridpairflag <= ph->nstyles) - vptr[nvirial++] = ph->styles[hybridpairflag-1]->virial; - } + if (hybridpairflag && force->pair) { + PairHybrid *ph = (PairHybrid *) force->pair; + ph->no_virial_fdotr_compute = 1; + vptr[nvirial++] = hybridpair->virial; } if (pairflag && force->pair) vptr[nvirial++] = force->pair->virial; if (bondflag && force->bond) vptr[nvirial++] = force->bond->virial;