Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F71957513
pair_peri_pmb_omp.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Jul 13, 23:07
Size
8 KB
Mime Type
text/x-c++
Expires
Mon, Jul 15, 23:07 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
19018450
Attached To
rLAMMPS lammps
pair_peri_pmb_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_pmb_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 "suffix.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairPeriPMBOMP::PairPeriPMBOMP(LAMMPS *lmp) :
PairPeriPMB(lmp), ThrOMP(lmp, THR_PAIR)
{
suffix_flag |= Suffix::OMP;
respa_enable = 0;
}
/* ---------------------------------------------------------------------- */
void PairPeriPMBOMP::compute(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else evflag = vflag_fdotr = 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);
nmax = atom->nmax;
memory->create(s0_new,nmax,"pair:s0_new");
}
#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);
thr->timer(Timer::START);
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);
}
thr->timer(Timer::PAIR);
reduce_thr(this, eflag, vflag, thr);
} // end of omp parallel region
}
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairPeriPMBOMP::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;
// 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;
rk = (15.0 * kspring[itype][jtype] * 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
// 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;
for (i = iifrom; i < iito; ++i) {
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jnum = npartner[i];
s0_new[i] = DBL_MAX;
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;
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;
stretch = dr / r0[i][jj];
rk = (kspring[itype][jtype] * vfrac[j]) * vfrac_scale * stretch;
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
if (EFLAG) evdwl = 0.5*rk*dr;
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
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];
}
/* ---------------------------------------------------------------------- */
double PairPeriPMBOMP::memory_usage()
{
double bytes = memory_usage_thr();
bytes += PairPeriPMB::memory_usage();
return bytes;
}
Event Timeline
Log In to Comment