Homec4science
Diffusion lammps (efficient_neuronet) Browse

lammps/src/USER-OMP/pair_peri_lps_omp.cpp
1eb5be159143efficient_neuronet

Last Change
Raw File

pair_peri_lps_omp.cpp
View Options

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "math.h"
#include "float.h"
#include "pair_peri_lps_omp.h"
#include "fix.h"
#include "fix_peri_neigh.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "force.h"
#include "memory.h"
#include "lattice.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "math_const.h"
#include "suffix.h"
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairPeriLPSOMP::PairPeriLPSOMP(LAMMPS *lmp) :
PairPeriLPS(lmp), ThrOMP(lmp, THR_PAIR)
{
suffix_flag |= Suffix::OMP;
respa_enable = 0;
}
/* ---------------------------------------------------------------------- */
void PairPeriLPSOMP::compute(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
const int nall = atom->nlocal + atom->nghost;
const int nthreads = comm->nthreads;
const int inum = list->inum;
// grow bond forces array if necessary
if (atom->nmax > nmax) {
memory->destroy(s0_new);
memory->destroy(theta);
nmax = atom->nmax;
memory->create(s0_new,nmax,"pair:s0_new");
memory->create(theta,nmax,"pair:theta");
}
#if defined(_OPENMP)
#pragma omp parallel default(none) shared(eflag,vflag)
#endif
{
int ifrom, ito, tid;
loop_setup_thr(ifrom, ito, tid, inum, nthreads);
ThrData *thr = fix->get_thr(tid);
ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
if (evflag) {
if (eflag) {
if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
else eval<1,1,0>(ifrom, ito, thr);
} else {
if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
else eval<1,0,0>(ifrom, ito, thr);
}
} else {
if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
else eval<0,0,0>(ifrom, ito, thr);
}
reduce_thr(this, eflag, vflag, thr);
} // end of omp parallel region
}
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairPeriLPSOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,ii,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz;
double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0,rsq0;
double rsq,r,dr,rk,evdwl,fpair,fbond;
int *ilist,*jlist,*numneigh,**firstneigh;
double d_ij,delta,stretch;
evdwl = 0.0;
const double * const * const x = atom->x;
double * const * const f = thr->get_f();
const int * const type = atom->type;
const int nlocal = atom->nlocal;
double fxtmp,fytmp,fztmp;
double *vfrac = atom->vfrac;
double *s0 = atom->s0;
double **x0 = atom->x0;
double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0;
tagint **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume;
// lc = lattice constant
// init_style guarantees it's the same in x, y, and z
double lc = domain->lattice->xlattice;
double half_lc = 0.5*lc;
double vfrac_scale = 1.0;
// short-range forces
int periodic = (domain->xperiodic || domain->yperiodic || domain->zperiodic);
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
// need minimg() for x0 difference since not ghosted
for (ii = iifrom; ii < iito; ++ii) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
xtmp0 = x0[i][0];
ytmp0 = x0[i][1];
ztmp0 = x0[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
fxtmp=fytmp=fztmp=0.0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
delx0 = xtmp0 - x0[j][0];
dely0 = ytmp0 - x0[j][1];
delz0 = ztmp0 - x0[j][2];
if (periodic) domain->minimum_image(delx0,dely0,delz0);
rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
jtype = type[j];
r = sqrt(rsq);
// short-range interaction distance based on initial particle position
// 0.9 and 1.35 are constants
d_ij = MIN(0.9*sqrt(rsq0),1.35*lc);
// short-range contact forces
// 15 is constant taken from the EMU Theory Manual
// Silling, 12 May 2005, p 18
if (r < d_ij) {
dr = r - d_ij;
// kshort based upon short-range force constant
// of the bond-based theory used in PMB model
double kshort = (15.0 * 18.0 * bulkmodulus[itype][itype]) /
(MY_PI * cutsq[itype][jtype] * cutsq[itype][jtype]);
rk = (kshort * vfrac[j]) * (dr / cut[itype][jtype]);
if (r > 0.0) fpair = -(rk/r);
else fpair = 0.0;
fxtmp += delx*fpair;
fytmp += dely*fpair;
fztmp += delz*fpair;
if (NEWTON_PAIR || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (EFLAG) evdwl = 0.5*rk*dr;
if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR,evdwl,0.0,
fpair*vfrac[i],delx,dely,delz,thr);
}
}
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
}
// wait until all threads are done since we
// need to distribute the work differently.
sync_threads();
#if defined(_OPENMP)
// each thread works on a fixed chunk of atoms.
const int idelta = 1 + nlocal/comm->nthreads;
iifrom = thr->get_tid()*idelta;
iito = ((iifrom + idelta) > nlocal) ? nlocal : (iifrom + idelta);
#else
iifrom = 0;
iito = nlocal;
#endif
// Compute the dilatation on each particle
if (iifrom < nlocal)
compute_dilatation_thr(iifrom, iito);
// wait until all threads are done before communication
sync_threads();
#if defined(_OPENMP)
#pragma omp master
#endif
{ // communicate dilatation (theta) of each particle
comm->forward_comm_pair(this);
// communicate weighted volume (wvolume) upon every reneighbor
if (neighbor->ago == 0)
comm->forward_comm_fix(modify->fix[ifix_peri]);
}
sync_threads();
// Volume-dependent part of the energy
if (EFLAG) {
for (i = iifrom; i < iito; i++) {
itype = type[i];
e_tally_thr(this, i, i, nlocal, NEWTON_PAIR,
0.5 * bulkmodulus[itype][itype] * (theta[i] * theta[i]), 0.0, thr);
}
}
// loop over my particles and their partners
// partner list contains all bond partners, so I-J appears twice
// if bond already broken, skip this partner
// first = true if this is first neighbor of particle i
bool first;
double omega_minus, omega_plus;
for (i = iifrom; i < iito; ++i) {
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
xtmp0 = x0[i][0];
ytmp0 = x0[i][1];
ztmp0 = x0[i][2];
itype = type[i];
jnum = npartner[i];
first = true;
for (jj = 0; jj < jnum; jj++) {
if (partner[i][jj] == 0) continue;
j = atom->map(partner[i][jj]);
// check if lost a partner without first breaking bond
if (j < 0) {
partner[i][jj] = 0;
continue;
}
// compute force density, add to PD equation of motion
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
if (periodic) domain->minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz;
delx0 = xtmp0 - x0[j][0];
dely0 = ytmp0 - x0[j][1];
delz0 = ztmp0 - x0[j][2];
if (periodic) domain->minimum_image(delx0,dely0,delz0);
jtype = type[j];
delta = cut[itype][jtype];
r = sqrt(rsq);
dr = r - r0[i][jj];
// avoid roundoff errors
if (fabs(dr) < 2.2204e-016) dr = 0.0;
// scale vfrac[j] if particle j near the horizon
if ((fabs(r0[i][jj] - delta)) <= half_lc)
vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) +
(1.0 + ((delta - half_lc)/(2*half_lc) ) );
else vfrac_scale = 1.0;
omega_plus = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0);
omega_minus = influence_function(delx0,dely0,delz0);
rk = ( (3.0 * bulkmodulus[itype][itype]) -
(5.0 * shearmodulus[itype][itype]) ) * vfrac[j] * vfrac_scale *
( (omega_plus * theta[i] / wvolume[i]) +
( omega_minus * theta[j] / wvolume[j] ) ) * r0[i][jj];
rk += 15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) *
( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr;
if (r > 0.0) fbond = -(rk/r);
else fbond = 0.0;
f[i][0] += delx*fbond;
f[i][1] += dely*fbond;
f[i][2] += delz*fbond;
// since I-J is double counted, set newton off & use 1/2 factor and I,I
double deviatoric_extension = dr - (theta[i]* r0[i][jj] / 3.0);
if (EFLAG) evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
omega_plus*(deviatoric_extension * deviatoric_extension) *
vfrac[j] * vfrac_scale;
if (EVFLAG) ev_tally_thr(this,i,i,nlocal,0,0.5*evdwl,0.0,
0.5*fbond*vfrac[i],delx,dely,delz,thr);
// find stretch in bond I-J and break if necessary
// use s0 from previous timestep
stretch = dr / r0[i][jj];
if (stretch > MIN(s0[i],s0[j])) partner[i][jj] = 0;
// update s0 for next timestep
if (first)
s0_new[i] = s00[itype][jtype] - (alpha[itype][jtype] * stretch);
else
s0_new[i] = MAX(s0_new[i],s00[itype][jtype] - (alpha[itype][jtype] * stretch));
first = false;
}
}
sync_threads();
// store new s0 (in parallel)
if (iifrom < nlocal)
for (i = iifrom; i < iito; i++) s0[i] = s0_new[i];
}
/* ---------------------------------------------------------------------- */
void PairPeriLPSOMP::compute_dilatation_thr(int ifrom, int ito)
{
int i,j,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz;
double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0;
double rsq,r,dr;
double delta;
double **x = atom->x;
int *type = atom->type;
double **x0 = atom->x0;
double *vfrac = atom->vfrac;
double vfrac_scale = 1.0;
double lc = domain->lattice->xlattice;
double half_lc = 0.5*lc;
double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0;
tagint **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume;
int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic;
// compute the dilatation theta
for (i = ifrom; i < ito; i++) {
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
xtmp0 = x0[i][0];
ytmp0 = x0[i][1];
ztmp0 = x0[i][2];
jnum = npartner[i];
theta[i] = 0.0;
itype = type[i];
for (jj = 0; jj < jnum; jj++) {
// if bond already broken, skip this partner
if (partner[i][jj] == 0) continue;
// Look up local index of this partner particle
j = atom->map(partner[i][jj]);
// Skip if particle is "lost"
if (j < 0) continue;
// Compute force density and add to PD equation of motion
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
if (periodic) domain->minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz;
delx0 = xtmp0 - x0[j][0];
dely0 = ytmp0 - x0[j][1];
delz0 = ztmp0 - x0[j][2];
if (periodic) domain->minimum_image(delx0,dely0,delz0);
r = sqrt(rsq);
dr = r - r0[i][jj];
if (fabs(dr) < 2.2204e-016) dr = 0.0;
jtype = type[j];
delta = cut[itype][jtype];
// scale vfrac[j] if particle j near the horizon
if ((fabs(r0[i][jj] - delta)) <= half_lc)
vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) +
(1.0 + ((delta - half_lc)/(2*half_lc) ) );
else vfrac_scale = 1.0;
theta[i] += influence_function(delx0, dely0, delz0) * r0[i][jj] * dr *
vfrac[j] * vfrac_scale;
}
// if wvolume[i] is zero, then particle i has no bonds
// therefore, the dilatation is set to
if (wvolume[i] != 0.0) theta[i] = (3.0/wvolume[i]) * theta[i];
else theta[i] = 0;
}
}
/* ---------------------------------------------------------------------- */
double PairPeriLPSOMP::memory_usage()
{
double bytes = memory_usage_thr();
bytes += PairPeriLPS::memory_usage();
return bytes;
}
c4science · Help