Homec4science
Diffusion lammps (master) Browse

lammps/src/USER-OMP/thr_omp.cpp
9a9291036d88master

Last Change
Raw File

thr_omp.cpp
View Options

/* -------------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
OpenMP based threading support for LAMMPS
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "thr_omp.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "math_const.h"
#include <string.h>
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
ThrOMP::ThrOMP(LAMMPS *ptr, int style)
: lmp(ptr), fix(NULL), thr_style(style), thr_error(0)
{
// register fix omp with this class
int ifix = lmp->modify->find_fix("package_omp");
if (ifix < 0)
lmp->error->all(FLERR,"The 'package omp' command is required for /omp styles");
fix = static_cast<FixOMP *>(lmp->modify->fix[ifix]);
}
/* ---------------------------------------------------------------------- */
ThrOMP::~ThrOMP()
{
// nothing to do?
}
/* ----------------------------------------------------------------------
Hook up per thread per atom arrays into the tally infrastructure
---------------------------------------------------------------------- */
void ThrOMP::ev_setup_thr(int eflag, int vflag, int nall, double *eatom,
double **vatom, ThrData *thr)
{
const int tid = thr->get_tid();
if (tid == 0) thr_error = 0;
if (thr_style & THR_PAIR) {
if (eflag & 2) {
thr->eatom_pair = eatom + tid*nall;
if (nall > 0)
memset(&(thr->eatom_pair[0]),0,nall*sizeof(double));
}
if (vflag & 4) {
thr->vatom_pair = vatom + tid*nall;
if (nall > 0)
memset(&(thr->vatom_pair[0][0]),0,nall*6*sizeof(double));
}
}
if (thr_style & THR_BOND) {
if (eflag & 2) {
thr->eatom_bond = eatom + tid*nall;
if (nall > 0)
memset(&(thr->eatom_bond[0]),0,nall*sizeof(double));
}
if (vflag & 4) {
thr->vatom_bond = vatom + tid*nall;
if (nall > 0)
memset(&(thr->vatom_bond[0][0]),0,nall*6*sizeof(double));
}
}
if (thr_style & THR_ANGLE) {
if (eflag & 2) {
thr->eatom_angle = eatom + tid*nall;
if (nall > 0)
memset(&(thr->eatom_angle[0]),0,nall*sizeof(double));
}
if (vflag & 4) {
thr->vatom_angle = vatom + tid*nall;
if (nall > 0)
memset(&(thr->vatom_angle[0][0]),0,nall*6*sizeof(double));
}
}
if (thr_style & THR_DIHEDRAL) {
if (eflag & 2) {
thr->eatom_dihed = eatom + tid*nall;
if (nall > 0)
memset(&(thr->eatom_dihed[0]),0,nall*sizeof(double));
}
if (vflag & 4) {
thr->vatom_dihed = vatom + tid*nall;
if (nall > 0)
memset(&(thr->vatom_dihed[0][0]),0,nall*6*sizeof(double));
}
}
if (thr_style & THR_IMPROPER) {
if (eflag & 2) {
thr->eatom_imprp = eatom + tid*nall;
if (nall > 0)
memset(&(thr->eatom_imprp[0]),0,nall*sizeof(double));
}
if (vflag & 4) {
thr->vatom_imprp = vatom + tid*nall;
if (nall > 0)
memset(&(thr->vatom_imprp[0][0]),0,nall*6*sizeof(double));
}
}
// nothing to do for THR_KSPACE
}
/* ----------------------------------------------------------------------
Reduce per thread data into the regular structures
Reduction of global properties is serialized with a "critical"
directive, so that only one thread at a time will access the
global variables. Since we are not synchronized, this should
come with little overhead. The reduction of per-atom properties
in contrast is parallelized over threads in the same way as forces.
---------------------------------------------------------------------- */
void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag,
ThrData *const thr, const int nproxy)
{
const int nlocal = lmp->atom->nlocal;
const int nghost = lmp->atom->nghost;
const int nall = nlocal + nghost;
const int nfirst = lmp->atom->nfirst;
const int nthreads = lmp->comm->nthreads;
const int evflag = eflag | vflag;
const int tid = thr->get_tid();
double **f = lmp->atom->f;
double **x = lmp->atom->x;
int need_force_reduce = 1;
if (evflag)
sync_threads();
switch (thr_style) {
case THR_PAIR: {
Pair * const pair = lmp->force->pair;
if (pair->vflag_fdotr) {
if (style == fix->last_pair_hybrid) {
// pair_style hybrid will compute fdotr for us
// but we first need to reduce the forces
data_reduce_thr(&(f[0][0]), nall, nthreads, 3, tid);
need_force_reduce = 0;
}
// this is a non-hybrid pair style. compute per thread fdotr
if (fix->last_pair_hybrid == NULL) {
if (lmp->neighbor->includegroup == 0)
thr->virial_fdotr_compute(x, nlocal, nghost, -1);
else
thr->virial_fdotr_compute(x, nlocal, nghost, nfirst);
}
}
if (evflag) {
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if (eflag & 1) {
pair->eng_vdwl += thr->eng_vdwl;
pair->eng_coul += thr->eng_coul;
thr->eng_vdwl = 0.0;
thr->eng_coul = 0.0;
}
if (vflag & 3)
for (int i=0; i < 6; ++i) {
pair->virial[i] += thr->virial_pair[i];
thr->virial_pair[i] = 0.0;
}
}
if (eflag & 2) {
data_reduce_thr(&(pair->eatom[0]), nall, nthreads, 1, tid);
}
if (vflag & 4) {
data_reduce_thr(&(pair->vatom[0][0]), nall, nthreads, 6, tid);
}
}
}
break;
case THR_PAIR|THR_PROXY: {
Pair * const pair = lmp->force->pair;
if (tid >= nproxy && pair->vflag_fdotr) {
if (fix->last_pair_hybrid) {
if (tid == nproxy)
lmp->error->all(FLERR,
"Cannot use hybrid pair style with kspace proxy");
else return;
}
// this is a non-hybrid pair style. compute per thread fdotr
if (fix->last_pair_hybrid == NULL) {
if (lmp->neighbor->includegroup == 0)
thr->virial_fdotr_compute(x, nlocal, nghost, -1);
else
thr->virial_fdotr_compute(x, nlocal, nghost, nfirst);
}
}
if (evflag) {
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if (tid < nproxy) {
// nothing to collect for kspace proxy threads
// just reset pair accumulators to 0.0.
if (eflag & 1) {
thr->eng_vdwl = 0.0;
thr->eng_coul = 0.0;
}
if (vflag & 3)
for (int i=0; i < 6; ++i) {
thr->virial_pair[i] = 0.0;
}
} else {
if (eflag & 1) {
pair->eng_vdwl += thr->eng_vdwl;
pair->eng_coul += thr->eng_coul;
thr->eng_vdwl = 0.0;
thr->eng_coul = 0.0;
}
if (vflag & 3)
for (int i=0; i < 6; ++i) {
pair->virial[i] += thr->virial_pair[i];
thr->virial_pair[i] = 0.0;
}
}
}
if (eflag & 2) {
data_reduce_thr(&(pair->eatom[0]), nall, nthreads, 1, tid);
}
if (vflag & 4) {
data_reduce_thr(&(pair->vatom[0][0]), nall, nthreads, 6, tid);
}
}
}
break;
case THR_BOND:
if (evflag) {
Bond * const bond = lmp->force->bond;
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if (eflag & 1) {
bond->energy += thr->eng_bond;
thr->eng_bond = 0.0;
}
if (vflag & 3) {
for (int i=0; i < 6; ++i) {
bond->virial[i] += thr->virial_bond[i];
thr->virial_bond[i] = 0.0;
}
}
}
if (eflag & 2) {
data_reduce_thr(&(bond->eatom[0]), nall, nthreads, 1, tid);
}
if (vflag & 4) {
data_reduce_thr(&(bond->vatom[0][0]), nall, nthreads, 6, tid);
}
}
break;
case THR_ANGLE:
if (evflag) {
Angle * const angle = lmp->force->angle;
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if (eflag & 1) {
angle->energy += thr->eng_angle;
thr->eng_angle = 0.0;
}
if (vflag & 3) {
for (int i=0; i < 6; ++i) {
angle->virial[i] += thr->virial_angle[i];
thr->virial_angle[i] = 0.0;
}
}
}
if (eflag & 2) {
data_reduce_thr(&(angle->eatom[0]), nall, nthreads, 1, tid);
}
if (vflag & 4) {
data_reduce_thr(&(angle->vatom[0][0]), nall, nthreads, 6, tid);
}
}
break;
case THR_DIHEDRAL:
if (evflag) {
Dihedral * const dihedral = lmp->force->dihedral;
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if (eflag & 1) {
dihedral->energy += thr->eng_dihed;
thr->eng_dihed = 0.0;
}
if (vflag & 3) {
for (int i=0; i < 6; ++i) {
dihedral->virial[i] += thr->virial_dihed[i];
thr->virial_dihed[i] = 0.0;
}
}
}
if (eflag & 2) {
data_reduce_thr(&(dihedral->eatom[0]), nall, nthreads, 1, tid);
}
if (vflag & 4) {
data_reduce_thr(&(dihedral->vatom[0][0]), nall, nthreads, 6, tid);
}
}
break;
case THR_DIHEDRAL|THR_CHARMM: // special case for CHARMM dihedrals
if (evflag) {
Dihedral * const dihedral = lmp->force->dihedral;
Pair * const pair = lmp->force->pair;
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if (eflag & 1) {
dihedral->energy += thr->eng_dihed;
pair->eng_vdwl += thr->eng_vdwl;
pair->eng_coul += thr->eng_coul;
thr->eng_dihed = 0.0;
thr->eng_vdwl = 0.0;
thr->eng_coul = 0.0;
}
if (vflag & 3) {
for (int i=0; i < 6; ++i) {
dihedral->virial[i] += thr->virial_dihed[i];
pair->virial[i] += thr->virial_pair[i];
thr->virial_dihed[i] = 0.0;
thr->virial_pair[i] = 0.0;
}
}
}
if (eflag & 2) {
data_reduce_thr(&(dihedral->eatom[0]), nall, nthreads, 1, tid);
data_reduce_thr(&(pair->eatom[0]), nall, nthreads, 1, tid);
}
if (vflag & 4) {
data_reduce_thr(&(dihedral->vatom[0][0]), nall, nthreads, 6, tid);
data_reduce_thr(&(pair->vatom[0][0]), nall, nthreads, 6, tid);
}
}
break;
case THR_IMPROPER:
if (evflag) {
Improper *improper = lmp->force->improper;
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if (eflag & 1) {
improper->energy += thr->eng_imprp;
thr->eng_imprp = 0.0;
}
if (vflag & 3) {
for (int i=0; i < 6; ++i) {
improper->virial[i] += thr->virial_imprp[i];
thr->virial_imprp[i] = 0.0;
}
}
}
if (eflag & 2) {
data_reduce_thr(&(improper->eatom[0]), nall, nthreads, 1, tid);
}
if (vflag & 4) {
data_reduce_thr(&(improper->vatom[0][0]), nall, nthreads, 6, tid);
}
}
break;
case THR_KSPACE|THR_PROXY: // fallthrough
case THR_KSPACE:
// nothing to do
break;
default:
printf("tid:%d unhandled thr_style case %d\n", tid, thr_style);
break;
}
if (style == fix->last_omp_style) {
if (need_force_reduce)
data_reduce_thr(&(f[0][0]), nall, nthreads, 3, tid);
if (lmp->atom->torque)
data_reduce_thr(&(lmp->atom->torque[0][0]), nall, nthreads, 3, tid);
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and eng_coul into per thread global and per-atom accumulators
------------------------------------------------------------------------- */
void ThrOMP::e_tally_thr(Pair * const pair, const int i, const int j,
const int nlocal, const int newton_pair,
const double evdwl, const double ecoul, ThrData * const thr)
{
if (pair->eflag_global) {
if (newton_pair) {
thr->eng_vdwl += evdwl;
thr->eng_coul += ecoul;
} else {
const double evdwlhalf = 0.5*evdwl;
const double ecoulhalf = 0.5*ecoul;
if (i < nlocal) {
thr->eng_vdwl += evdwlhalf;
thr->eng_coul += ecoulhalf;
}
if (j < nlocal) {
thr->eng_vdwl += evdwlhalf;
thr->eng_coul += ecoulhalf;
}
}
}
if (pair->eflag_atom) {
const double epairhalf = 0.5 * (evdwl + ecoul);
if (newton_pair || i < nlocal) thr->eatom_pair[i] += epairhalf;
if (newton_pair || j < nlocal) thr->eatom_pair[j] += epairhalf;
}
}
/* helper functions */
static void v_tally(double * const vout, const double * const vin)
{
vout[0] += vin[0];
vout[1] += vin[1];
vout[2] += vin[2];
vout[3] += vin[3];
vout[4] += vin[4];
vout[5] += vin[5];
}
static void v_tally(double * const vout, const double scale, const double * const vin)
{
vout[0] += scale*vin[0];
vout[1] += scale*vin[1];
vout[2] += scale*vin[2];
vout[3] += scale*vin[3];
vout[4] += scale*vin[4];
vout[5] += scale*vin[5];
}
/* ----------------------------------------------------------------------
tally virial into per thread global and per-atom accumulators
------------------------------------------------------------------------- */
void ThrOMP::v_tally_thr(Pair * const pair, const int i, const int j,
const int nlocal, const int newton_pair,
const double * const v, ThrData * const thr)
{
if (pair->vflag_global) {
double * const va = thr->virial_pair;
if (newton_pair) {
v_tally(va,v);
} else {
if (i < nlocal) v_tally(va,0.5,v);
if (j < nlocal) v_tally(va,0.5,v);
}
}
if (pair->vflag_atom) {
if (newton_pair || i < nlocal) {
double * const va = thr->vatom_pair[i];
v_tally(va,0.5,v);
}
if (newton_pair || j < nlocal) {
double * const va = thr->vatom_pair[j];
v_tally(va,0.5,v);
}
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into per thread global and per-atom accumulators
need i < nlocal test since called by bond_quartic and dihedral_charmm
------------------------------------------------------------------------- */
void ThrOMP::ev_tally_thr(Pair * const pair, const int i, const int j, const int nlocal,
const int newton_pair, const double evdwl, const double ecoul,
const double fpair, const double delx, const double dely,
const double delz, ThrData * const thr)
{
if (pair->eflag_either)
e_tally_thr(pair, i, j, nlocal, newton_pair, evdwl, ecoul, thr);
if (pair->vflag_either) {
double v[6];
v[0] = delx*delx*fpair;
v[1] = dely*dely*fpair;
v[2] = delz*delz*fpair;
v[3] = delx*dely*fpair;
v[4] = delx*delz*fpair;
v[5] = dely*delz*fpair;
v_tally_thr(pair, i, j, nlocal, newton_pair, v, thr);
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into global and per-atom accumulators
for virial, have delx,dely,delz and fx,fy,fz
------------------------------------------------------------------------- */
void ThrOMP::ev_tally_xyz_thr(Pair * const pair, const int i, const int j,
const int nlocal, const int newton_pair,
const double evdwl, const double ecoul,
const double fx, const double fy, const double fz,
const double delx, const double dely, const double delz,
ThrData * const thr)
{
if (pair->eflag_either)
e_tally_thr(pair, i, j, nlocal, newton_pair, evdwl, ecoul, thr);
if (pair->vflag_either) {
double v[6];
v[0] = delx*fx;
v[1] = dely*fy;
v[2] = delz*fz;
v[3] = delx*fy;
v[4] = delx*fz;
v[5] = dely*fz;
v_tally_thr(pair, i, j, nlocal, newton_pair, v, thr);
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into global and per-atom accumulators
called by SW and hbond potentials, newton_pair is always on
virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk
------------------------------------------------------------------------- */
void ThrOMP::ev_tally3_thr(Pair * const pair, const int i, const int j, const int k,
const double evdwl, const double ecoul,
const double * const fj, const double * const fk,
const double * const drji, const double * const drki,
ThrData * const thr)
{
if (pair->eflag_either) {
if (pair->eflag_global) {
thr->eng_vdwl += evdwl;
thr->eng_coul += ecoul;
}
if (pair->eflag_atom) {
const double epairthird = THIRD * (evdwl + ecoul);
thr->eatom_pair[i] += epairthird;
thr->eatom_pair[j] += epairthird;
thr->eatom_pair[k] += epairthird;
}
}
if (pair->vflag_either) {
double v[6];
v[0] = drji[0]*fj[0] + drki[0]*fk[0];
v[1] = drji[1]*fj[1] + drki[1]*fk[1];
v[2] = drji[2]*fj[2] + drki[2]*fk[2];
v[3] = drji[0]*fj[1] + drki[0]*fk[1];
v[4] = drji[0]*fj[2] + drki[0]*fk[2];
v[5] = drji[1]*fj[2] + drki[1]*fk[2];
if (pair->vflag_global) v_tally(thr->virial_pair,v);
if (pair->vflag_atom) {
v_tally(thr->vatom_pair[i],THIRD,v);
v_tally(thr->vatom_pair[j],THIRD,v);
v_tally(thr->vatom_pair[k],THIRD,v);
}
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into global and per-atom accumulators
called by AIREBO potential, newton_pair is always on
------------------------------------------------------------------------- */
void ThrOMP::ev_tally4_thr(Pair * const pair, const int i, const int j,
const int k, const int m, const double evdwl,
const double * const fi, const double * const fj,
const double * const fk, const double * const drim,
const double * const drjm, const double * const drkm,
ThrData * const thr)
{
double v[6];
if (pair->eflag_either) {
if (pair->eflag_global) thr->eng_vdwl += evdwl;
if (pair->eflag_atom) {
const double epairfourth = 0.25 * evdwl;
thr->eatom_pair[i] += epairfourth;
thr->eatom_pair[j] += epairfourth;
thr->eatom_pair[k] += epairfourth;
thr->eatom_pair[m] += epairfourth;
}
}
if (pair->vflag_atom) {
v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]);
v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]);
v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]);
v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]);
v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]);
v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]);
v_tally(thr->vatom_pair[i],v);
v_tally(thr->vatom_pair[j],v);
v_tally(thr->vatom_pair[k],v);
v_tally(thr->vatom_pair[m],v);
}
}
/* ----------------------------------------------------------------------
tally ecoul and virial into each of n atoms in list
called by TIP4P potential, newton_pair is always on
changes v values by dividing by n
------------------------------------------------------------------------- */
void ThrOMP::ev_tally_list_thr(Pair * const pair, const int n,
const int * const list, const double ecoul,
const double * const v, ThrData * const thr)
{
if (pair->eflag_either) {
if (pair->eflag_global) thr->eng_coul += ecoul;
if (pair->eflag_atom) {
double epairatom = ecoul/static_cast<double>(n);
for (int i = 0; i < n; i++) thr->eatom_pair[list[i]] += epairatom;
}
}
if (pair->vflag_either) {
if (pair->vflag_global)
v_tally(thr->virial_pair,v);
if (pair->vflag_atom) {
const double s = 1.0/static_cast<double>(n);
double vtmp[6];
vtmp[0] = s * v[0];
vtmp[1] = s * v[1];
vtmp[2] = s * v[2];
vtmp[3] = s * v[3];
vtmp[4] = s * v[4];
vtmp[5] = s * v[5];
for (int i = 0; i < n; i++) {
const int j = list[i];
v_tally(thr->vatom_pair[j],vtmp);
}
}
}
}
/* ----------------------------------------------------------------------
tally energy and virial into global and per-atom accumulators
------------------------------------------------------------------------- */
void ThrOMP::ev_tally_thr(Bond * const bond, const int i, const int j, const int nlocal,
const int newton_bond, const double ebond, const double fbond,
const double delx, const double dely, const double delz,
ThrData * const thr)
{
if (bond->eflag_either) {
const double ebondhalf = 0.5*ebond;
if (newton_bond) {
if (bond->eflag_global)
thr->eng_bond += ebond;
if (bond->eflag_atom) {
thr->eatom_bond[i] += ebondhalf;
thr->eatom_bond[j] += ebondhalf;
}
} else {
if (bond->eflag_global) {
if (i < nlocal) thr->eng_bond += ebondhalf;
if (j < nlocal) thr->eng_bond += ebondhalf;
}
if (bond->eflag_atom) {
if (i < nlocal) thr->eatom_bond[i] += ebondhalf;
if (j < nlocal) thr->eatom_bond[j] += ebondhalf;
}
}
}
if (bond->vflag_either) {
double v[6];
v[0] = delx*delx*fbond;
v[1] = dely*dely*fbond;
v[2] = delz*delz*fbond;
v[3] = delx*dely*fbond;
v[4] = delx*delz*fbond;
v[5] = dely*delz*fbond;
if (bond->vflag_global) {
if (newton_bond)
v_tally(thr->virial_bond,v);
else {
if (i < nlocal)
v_tally(thr->virial_bond,0.5,v);
if (j < nlocal)
v_tally(thr->virial_bond,0.5,v);
}
}
if (bond->vflag_atom) {
v[0] *= 0.5;
v[1] *= 0.5;
v[2] *= 0.5;
v[3] *= 0.5;
v[4] *= 0.5;
v[5] *= 0.5;
if (newton_bond) {
v_tally(thr->vatom_bond[i],v);
v_tally(thr->vatom_bond[j],v);
} else {
if (j < nlocal)
v_tally(thr->vatom_bond[i],v);
if (j < nlocal)
v_tally(thr->vatom_bond[j],v);
}
}
}
}
/* ----------------------------------------------------------------------
tally energy and virial into global and per-atom accumulators
virial = r1F1 + r2F2 + r3F3 = (r1-r2) F1 + (r3-r2) F3 = del1*f1 + del2*f3
------------------------------------------------------------------------- */
void ThrOMP::ev_tally_thr(Angle * const angle, const int i, const int j, const int k,
const int nlocal, const int newton_bond, const double eangle,
const double * const f1, const double * const f3,
const double delx1, const double dely1, const double delz1,
const double delx2, const double dely2, const double delz2,
ThrData * const thr)
{
if (angle->eflag_either) {
const double eanglethird = THIRD*eangle;
if (newton_bond) {
if (angle->eflag_global)
thr->eng_angle += eangle;
if (angle->eflag_atom) {
thr->eatom_angle[i] += eanglethird;
thr->eatom_angle[j] += eanglethird;
thr->eatom_angle[k] += eanglethird;
}
} else {
if (angle->eflag_global) {
if (i < nlocal) thr->eng_angle += eanglethird;
if (j < nlocal) thr->eng_angle += eanglethird;
if (k < nlocal) thr->eng_angle += eanglethird;
}
if (angle->eflag_atom) {
if (i < nlocal) thr->eatom_angle[i] += eanglethird;
if (j < nlocal) thr->eatom_angle[j] += eanglethird;
if (k < nlocal) thr->eatom_angle[k] += eanglethird;
}
}
}
if (angle->vflag_either) {
double v[6];
v[0] = delx1*f1[0] + delx2*f3[0];
v[1] = dely1*f1[1] + dely2*f3[1];
v[2] = delz1*f1[2] + delz2*f3[2];
v[3] = delx1*f1[1] + delx2*f3[1];
v[4] = delx1*f1[2] + delx2*f3[2];
v[5] = dely1*f1[2] + dely2*f3[2];
if (angle->vflag_global) {
if (newton_bond) {
v_tally(thr->virial_angle,v);
} else {
int cnt = 0;
if (i < nlocal) ++cnt;
if (j < nlocal) ++cnt;
if (k < nlocal) ++cnt;
v_tally(thr->virial_angle,cnt*THIRD,v);
}
}
if (angle->vflag_atom) {
v[0] *= THIRD;
v[1] *= THIRD;
v[2] *= THIRD;
v[3] *= THIRD;
v[4] *= THIRD;
v[5] *= THIRD;
if (newton_bond) {
v_tally(thr->vatom_angle[i],v);
v_tally(thr->vatom_angle[j],v);
v_tally(thr->vatom_angle[k],v);
} else {
if (j < nlocal) v_tally(thr->vatom_angle[i],v);
if (j < nlocal) v_tally(thr->vatom_angle[j],v);
if (k < nlocal) v_tally(thr->vatom_angle[k],v);
}
}
}
}
/* ----------------------------------------------------------------------
tally energy and virial from 1-3 repulsion of SDK angle into accumulators
------------------------------------------------------------------------- */
void ThrOMP::ev_tally13_thr(Angle * const angle, const int i1, const int i3,
const int nlocal, const int newton_bond,
const double epair, const double fpair,
const double delx, const double dely,
const double delz, ThrData * const thr)
{
if (angle->eflag_either) {
const double epairhalf = 0.5 * epair;
if (angle->eflag_global) {
if (newton_bond || i1 < nlocal)
thr->eng_angle += epairhalf;
if (newton_bond || i3 < nlocal)
thr->eng_angle += epairhalf;
}
if (angle->eflag_atom) {
if (newton_bond || i1 < nlocal) thr->eatom_angle[i1] += epairhalf;
if (newton_bond || i3 < nlocal) thr->eatom_angle[i3] += epairhalf;
}
}
if (angle->vflag_either) {
double v[6];
v[0] = delx*delx*fpair;
v[1] = dely*dely*fpair;
v[2] = delz*delz*fpair;
v[3] = delx*dely*fpair;
v[4] = delx*delz*fpair;
v[5] = dely*delz*fpair;
if (angle->vflag_global) {
double * const va = thr->virial_angle;
if (newton_bond || i1 < nlocal) v_tally(va,0.5,v);
if (newton_bond || i3 < nlocal) v_tally(va,0.5,v);
}
if (angle->vflag_atom) {
if (newton_bond || i1 < nlocal) {
double * const va = thr->vatom_angle[i1];
v_tally(va,0.5,v);
}
if (newton_bond || i3 < nlocal) {
double * const va = thr->vatom_angle[i3];
v_tally(va,0.5,v);
}
}
}
}
/* ----------------------------------------------------------------------
tally energy and virial into global and per-atom accumulators
virial = r1F1 + r2F2 + r3F3 + r4F4 = (r1-r2) F1 + (r3-r2) F3 + (r4-r2) F4
= (r1-r2) F1 + (r3-r2) F3 + (r4-r3 + r3-r2) F4
= vb1*f1 + vb2*f3 + (vb3+vb2)*f4
------------------------------------------------------------------------- */
void ThrOMP::ev_tally_thr(Dihedral * const dihed, const int i1, const int i2,
const int i3, const int i4, const int nlocal,
const int newton_bond, const double edihedral,
const double * const f1, const double * const f3,
const double * const f4, const double vb1x,
const double vb1y, const double vb1z, const double vb2x,
const double vb2y, const double vb2z, const double vb3x,
const double vb3y, const double vb3z, ThrData * const thr)
{
if (dihed->eflag_either) {
if (dihed->eflag_global) {
if (newton_bond) {
thr->eng_dihed += edihedral;
} else {
const double edihedralquarter = 0.25*edihedral;
int cnt = 0;
if (i1 < nlocal) ++cnt;
if (i2 < nlocal) ++cnt;
if (i3 < nlocal) ++cnt;
if (i4 < nlocal) ++cnt;
thr->eng_dihed += static_cast<double>(cnt)*edihedralquarter;
}
}
if (dihed->eflag_atom) {
const double edihedralquarter = 0.25*edihedral;
if (newton_bond) {
thr->eatom_dihed[i1] += edihedralquarter;
thr->eatom_dihed[i2] += edihedralquarter;
thr->eatom_dihed[i3] += edihedralquarter;
thr->eatom_dihed[i4] += edihedralquarter;
} else {
if (i1 < nlocal) thr->eatom_dihed[i1] += edihedralquarter;
if (i2 < nlocal) thr->eatom_dihed[i2] += edihedralquarter;
if (i3 < nlocal) thr->eatom_dihed[i3] += edihedralquarter;
if (i4 < nlocal) thr->eatom_dihed[i4] += edihedralquarter;
}
}
}
if (dihed->vflag_either) {
double v[6];
v[0] = vb1x*f1[0] + vb2x*f3[0] + (vb3x+vb2x)*f4[0];
v[1] = vb1y*f1[1] + vb2y*f3[1] + (vb3y+vb2y)*f4[1];
v[2] = vb1z*f1[2] + vb2z*f3[2] + (vb3z+vb2z)*f4[2];
v[3] = vb1x*f1[1] + vb2x*f3[1] + (vb3x+vb2x)*f4[1];
v[4] = vb1x*f1[2] + vb2x*f3[2] + (vb3x+vb2x)*f4[2];
v[5] = vb1y*f1[2] + vb2y*f3[2] + (vb3y+vb2y)*f4[2];
if (dihed->vflag_global) {
if (newton_bond) {
v_tally(thr->virial_dihed,v);
} else {
int cnt = 0;
if (i1 < nlocal) ++cnt;
if (i2 < nlocal) ++cnt;
if (i3 < nlocal) ++cnt;
if (i4 < nlocal) ++cnt;
v_tally(thr->virial_dihed,0.25*static_cast<double>(cnt),v);
}
}
v[0] *= 0.25;
v[1] *= 0.25;
v[2] *= 0.25;
v[3] *= 0.25;
v[4] *= 0.25;
v[5] *= 0.25;
if (dihed->vflag_atom) {
if (newton_bond) {
v_tally(thr->vatom_dihed[i1],v);
v_tally(thr->vatom_dihed[i2],v);
v_tally(thr->vatom_dihed[i3],v);
v_tally(thr->vatom_dihed[i4],v);
} else {
if (i1 < nlocal) v_tally(thr->vatom_dihed[i1],v);
if (i2 < nlocal) v_tally(thr->vatom_dihed[i2],v);
if (i3 < nlocal) v_tally(thr->vatom_dihed[i3],v);
if (i4 < nlocal) v_tally(thr->vatom_dihed[i4],v);
}
}
}
}
/* ----------------------------------------------------------------------
tally energy and virial into global and per-atom accumulators
virial = r1F1 + r2F2 + r3F3 + r4F4 = (r1-r2) F1 + (r3-r2) F3 + (r4-r2) F4
= (r1-r2) F1 + (r3-r2) F3 + (r4-r3 + r3-r2) F4
= vb1*f1 + vb2*f3 + (vb3+vb2)*f4
------------------------------------------------------------------------- */
void ThrOMP::ev_tally_thr(Improper * const imprp, const int i1, const int i2,
const int i3, const int i4, const int nlocal,
const int newton_bond, const double eimproper,
const double * const f1, const double * const f3,
const double * const f4, const double vb1x,
const double vb1y, const double vb1z, const double vb2x,
const double vb2y, const double vb2z, const double vb3x,
const double vb3y, const double vb3z, ThrData * const thr)
{
if (imprp->eflag_either) {
if (imprp->eflag_global) {
if (newton_bond) {
thr->eng_imprp += eimproper;
} else {
const double eimproperquarter = 0.25*eimproper;
int cnt = 0;
if (i1 < nlocal) ++cnt;
if (i2 < nlocal) ++cnt;
if (i3 < nlocal) ++cnt;
if (i4 < nlocal) ++cnt;
thr->eng_imprp += static_cast<double>(cnt)*eimproperquarter;
}
}
if (imprp->eflag_atom) {
const double eimproperquarter = 0.25*eimproper;
if (newton_bond) {
thr->eatom_imprp[i1] += eimproperquarter;
thr->eatom_imprp[i2] += eimproperquarter;
thr->eatom_imprp[i3] += eimproperquarter;
thr->eatom_imprp[i4] += eimproperquarter;
} else {
if (i1 < nlocal) thr->eatom_imprp[i1] += eimproperquarter;
if (i2 < nlocal) thr->eatom_imprp[i2] += eimproperquarter;
if (i3 < nlocal) thr->eatom_imprp[i3] += eimproperquarter;
if (i4 < nlocal) thr->eatom_imprp[i4] += eimproperquarter;
}
}
}
if (imprp->vflag_either) {
double v[6];
v[0] = vb1x*f1[0] + vb2x*f3[0] + (vb3x+vb2x)*f4[0];
v[1] = vb1y*f1[1] + vb2y*f3[1] + (vb3y+vb2y)*f4[1];
v[2] = vb1z*f1[2] + vb2z*f3[2] + (vb3z+vb2z)*f4[2];
v[3] = vb1x*f1[1] + vb2x*f3[1] + (vb3x+vb2x)*f4[1];
v[4] = vb1x*f1[2] + vb2x*f3[2] + (vb3x+vb2x)*f4[2];
v[5] = vb1y*f1[2] + vb2y*f3[2] + (vb3y+vb2y)*f4[2];
if (imprp->vflag_global) {
if (newton_bond) {
v_tally(thr->virial_imprp,v);
} else {
int cnt = 0;
if (i1 < nlocal) ++cnt;
if (i2 < nlocal) ++cnt;
if (i3 < nlocal) ++cnt;
if (i4 < nlocal) ++cnt;
v_tally(thr->virial_imprp,0.25*static_cast<double>(cnt),v);
}
}
v[0] *= 0.25;
v[1] *= 0.25;
v[2] *= 0.25;
v[3] *= 0.25;
v[4] *= 0.25;
v[5] *= 0.25;
if (imprp->vflag_atom) {
if (newton_bond) {
v_tally(thr->vatom_imprp[i1],v);
v_tally(thr->vatom_imprp[i2],v);
v_tally(thr->vatom_imprp[i3],v);
v_tally(thr->vatom_imprp[i4],v);
} else {
if (i1 < nlocal) v_tally(thr->vatom_imprp[i1],v);
if (i2 < nlocal) v_tally(thr->vatom_imprp[i2],v);
if (i3 < nlocal) v_tally(thr->vatom_imprp[i3],v);
if (i4 < nlocal) v_tally(thr->vatom_imprp[i4],v);
}
}
}
}
/* ----------------------------------------------------------------------
tally virial into per-atom accumulators
called by AIREBO potential, newton_pair is always on
fpair is magnitude of force on atom I
------------------------------------------------------------------------- */
void ThrOMP::v_tally2_thr(const int i, const int j, const double fpair,
const double * const drij, ThrData * const thr)
{
double v[6];
v[0] = 0.5 * drij[0]*drij[0]*fpair;
v[1] = 0.5 * drij[1]*drij[1]*fpair;
v[2] = 0.5 * drij[2]*drij[2]*fpair;
v[3] = 0.5 * drij[0]*drij[1]*fpair;
v[4] = 0.5 * drij[0]*drij[2]*fpair;
v[5] = 0.5 * drij[1]*drij[2]*fpair;
v_tally(thr->vatom_pair[i],v);
v_tally(thr->vatom_pair[j],v);
}
/* ----------------------------------------------------------------------
tally virial into per-atom accumulators
called by AIREBO and Tersoff potential, newton_pair is always on
------------------------------------------------------------------------- */
void ThrOMP::v_tally3_thr(const int i, const int j, const int k,
const double * const fi, const double * const fj,
const double * const drik, const double * const drjk,
ThrData * const thr)
{
double v[6];
v[0] = THIRD * (drik[0]*fi[0] + drjk[0]*fj[0]);
v[1] = THIRD * (drik[1]*fi[1] + drjk[1]*fj[1]);
v[2] = THIRD * (drik[2]*fi[2] + drjk[2]*fj[2]);
v[3] = THIRD * (drik[0]*fi[1] + drjk[0]*fj[1]);
v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]);
v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]);
v_tally(thr->vatom_pair[i],v);
v_tally(thr->vatom_pair[j],v);
v_tally(thr->vatom_pair[k],v);
}
/* ----------------------------------------------------------------------
tally virial into per-atom accumulators
called by AIREBO potential, newton_pair is always on
------------------------------------------------------------------------- */
void ThrOMP::v_tally4_thr(const int i, const int j, const int k, const int m,
const double * const fi, const double * const fj,
const double * const fk, const double * const drim,
const double * const drjm, const double * const drkm,
ThrData * const thr)
{
double v[6];
v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]);
v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]);
v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]);
v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]);
v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]);
v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]);
v_tally(thr->vatom_pair[i],v);
v_tally(thr->vatom_pair[j],v);
v_tally(thr->vatom_pair[k],v);
v_tally(thr->vatom_pair[m],v);
}
/* ---------------------------------------------------------------------- */
double ThrOMP::memory_usage_thr()
{
double bytes=0.0;
return bytes;
}
c4science · Help