mirror of https://github.com/lammps/lammps.git
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9280 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
b078dc5e0b
commit
56a3ae0eca
|
|
@ -0,0 +1,244 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing authors: Axel Kohlmeyer (Temple U), Stan Moore (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "msm_omp.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "math_const.h"
|
||||
|
||||
#include <string.h>
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#include <omp.h>
|
||||
#endif
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
MSMOMP::MSMOMP(LAMMPS *lmp, int narg, char **arg) :
|
||||
MSM(lmp, narg, arg), ThrOMP(lmp, THR_KSPACE)
|
||||
{
|
||||
suffix_flag |= Suffix::OMP;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
run the regular toplevel compute method from plain PPPPM
|
||||
which will have individual methods replaced by our threaded
|
||||
versions and then call the obligatory force reduction.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void MSMOMP::compute(int eflag, int vflag)
|
||||
{
|
||||
|
||||
MSM::compute(eflag,vflag);
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp parallel default(none) shared(eflag,vflag)
|
||||
#endif
|
||||
{
|
||||
#if defined(_OPENMP)
|
||||
const int tid = omp_get_thread_num();
|
||||
#else
|
||||
const int tid = 0;
|
||||
#endif
|
||||
ThrData *thr = fix->get_thr(tid);
|
||||
reduce_thr(this, eflag, vflag, thr);
|
||||
} // end of omp parallel region
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
MSM direct part procedure for intermediate grid levels
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void MSMOMP::direct(int n)
|
||||
{
|
||||
// zero out electric potential
|
||||
|
||||
memset(&(egrid[n][nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
|
||||
|
||||
// zero out virial
|
||||
|
||||
if (vflag_atom) {
|
||||
memset(&(v0grid[n][nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
|
||||
memset(&(v1grid[n][nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
|
||||
memset(&(v2grid[n][nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
|
||||
memset(&(v3grid[n][nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
|
||||
memset(&(v4grid[n][nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
|
||||
memset(&(v5grid[n][nzlo_out[n]][nylo_out[n]][nxlo_out[n]]),0,ngrid[n]*sizeof(double));
|
||||
}
|
||||
|
||||
if (eflag_global) {
|
||||
if (vflag_global) {
|
||||
if (vflag_atom)
|
||||
direct_eval<1,1,1>(n);
|
||||
else
|
||||
direct_eval<1,1,0>(n);
|
||||
} else {
|
||||
if (vflag_atom)
|
||||
direct_eval<1,0,1>(n);
|
||||
else
|
||||
direct_eval<1,0,0>(n);
|
||||
}
|
||||
} else { // !eflag_global
|
||||
if (vflag_global) {
|
||||
if (vflag_atom)
|
||||
direct_eval<0,1,1>(n);
|
||||
else
|
||||
direct_eval<0,1,0>(n);
|
||||
} else {
|
||||
if (vflag_atom)
|
||||
direct_eval<0,0,1>(n);
|
||||
else
|
||||
direct_eval<0,0,0>(n);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <int EFLAG_GLOBAL, int VFLAG_GLOBAL, int VFLAG_ATOM>
|
||||
void MSMOMP::direct_eval(const int n)
|
||||
{
|
||||
double * const * const * const egridn = egrid[n];
|
||||
const double * const * const * const qgridn = qgrid[n];
|
||||
const double * const g_directn = g_direct[n];
|
||||
const double * const v0_directn = v0_direct[n];
|
||||
const double * const v1_directn = v1_direct[n];
|
||||
const double * const v2_directn = v2_direct[n];
|
||||
const double * const v3_directn = v3_direct[n];
|
||||
const double * const v4_directn = v4_direct[n];
|
||||
const double * const v5_directn = v5_direct[n];
|
||||
|
||||
double v0,v1,v2,v3,v4,v5,emsm;
|
||||
v0 = v1 = v2 = v3 = v4 = v5 = emsm = 0.0;
|
||||
const int alphan = alpha[n];
|
||||
const int betaxn = betax[n];
|
||||
const int betayn = betay[n];
|
||||
const int betazn = betaz[n];
|
||||
|
||||
const int nx = nxhi_direct - nxlo_direct + 1;
|
||||
const int ny = nyhi_direct - nylo_direct + 1;
|
||||
|
||||
// merge three outer loops into one for better threading
|
||||
|
||||
const int nzlo_inn = nzlo_in[n];
|
||||
const int nylo_inn = nylo_in[n];
|
||||
const int nxlo_inn = nxlo_in[n];
|
||||
const int numz = nzhi_in[n] - nzlo_inn + 1;
|
||||
const int numy = nyhi_in[n] - nylo_inn + 1;
|
||||
const int numx = nxhi_in[n] - nxlo_inn + 1;
|
||||
const int inum = numz*numy*numx;
|
||||
|
||||
const int zper = domain->zperiodic;
|
||||
const int yper = domain->yperiodic;
|
||||
const int xper = domain->xperiodic;
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp parallel default(none) reduction(+:v0,v1,v2,v3,v4,v5,emsm)
|
||||
#endif
|
||||
{
|
||||
double qtmp,esum,v0sum,v1sum,v2sum,v3sum,v4sum,v5sum;
|
||||
int i,ifrom,ito,tid,icx,icy,icz,ix,iy,iz,k;
|
||||
|
||||
loop_setup_thr(ifrom, ito, tid, inum, comm->nthreads);
|
||||
|
||||
for (i = ifrom; i < ito; ++i) {
|
||||
|
||||
// infer outer loop indices icx, icy, icz from master loop index
|
||||
|
||||
icz = i/(numy*numx);
|
||||
icy = (i - icz*numy*numx) / numx;
|
||||
icx = i - icz*numy*numx - icy*numx;
|
||||
icz += nzlo_inn;
|
||||
icy += nylo_inn;
|
||||
icx += nxlo_inn;
|
||||
|
||||
const int kmin = zper ? nzlo_direct : MAX(nzlo_direct,alphan - icz);
|
||||
const int kmax = zper ? nzhi_direct : MIN(nzhi_direct,betazn - icz);
|
||||
const int jmin = yper ? nylo_direct : MAX(nylo_direct,alphan - icy);
|
||||
const int jmax = yper ? nyhi_direct : MIN(nyhi_direct,betayn - icy);
|
||||
const int imin = xper ? nxlo_direct : MAX(nxlo_direct,alphan - icx);
|
||||
const int imax = xper ? nxhi_direct : MIN(nxhi_direct,betaxn - icx);
|
||||
|
||||
esum = 0.0;
|
||||
if (VFLAG_GLOBAL || VFLAG_ATOM)
|
||||
v0sum = v1sum = v2sum = v3sum = v4sum = v5sum = 0.0;
|
||||
|
||||
for (iz = kmin; iz <= kmax; iz++) {
|
||||
const int kk = icz+iz;
|
||||
const int zk = (iz + nzhi_direct)*ny;
|
||||
for (iy = jmin; iy <= jmax; iy++) {
|
||||
const int jj = icy+iy;
|
||||
const int zyk = (zk + iy + nyhi_direct)*nx;
|
||||
for (ix = imin; ix <= imax; ix++) {
|
||||
qtmp = qgridn[kk][jj][icx+ix];
|
||||
k = zyk + ix + nxhi_direct;
|
||||
esum += g_directn[k] * qtmp;
|
||||
|
||||
if (VFLAG_GLOBAL || VFLAG_ATOM) {
|
||||
v0sum += v0_directn[k] * qtmp;
|
||||
v1sum += v1_directn[k] * qtmp;
|
||||
v2sum += v2_directn[k] * qtmp;
|
||||
v3sum += v3_directn[k] * qtmp;
|
||||
v4sum += v4_directn[k] * qtmp;
|
||||
v5sum += v5_directn[k] * qtmp;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
egridn[icz][icy][icx] = esum;
|
||||
|
||||
if (VFLAG_ATOM) {
|
||||
v0grid[n][icz][icy][icx] = v0sum;
|
||||
v1grid[n][icz][icy][icx] = v1sum;
|
||||
v2grid[n][icz][icy][icx] = v2sum;
|
||||
v3grid[n][icz][icy][icx] = v3sum;
|
||||
v4grid[n][icz][icy][icx] = v4sum;
|
||||
v5grid[n][icz][icy][icx] = v5sum;
|
||||
}
|
||||
|
||||
if (EFLAG_GLOBAL || VFLAG_GLOBAL) {
|
||||
qtmp = qgridn[icz][icy][icx];
|
||||
if (EFLAG_GLOBAL) emsm += esum * qtmp;
|
||||
if (VFLAG_GLOBAL) {
|
||||
v0 += v0sum * qtmp;
|
||||
v1 += v1sum * qtmp;
|
||||
v2 += v2sum * qtmp;
|
||||
v3 += v3sum * qtmp;
|
||||
v4 += v4sum * qtmp;
|
||||
v5 += v5sum * qtmp;
|
||||
}
|
||||
}
|
||||
}
|
||||
} // end of omp parallel region
|
||||
|
||||
if (EFLAG_GLOBAL || VFLAG_GLOBAL) {
|
||||
if (EFLAG_GLOBAL) energy += emsm;
|
||||
if (VFLAG_GLOBAL) {
|
||||
virial[0] += v0;
|
||||
virial[1] += v1;
|
||||
virial[2] += v2;
|
||||
virial[3] += v3;
|
||||
virial[4] += v4;
|
||||
virial[5] += v5;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,45 @@
|
|||
/* -*- 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 KSPACE_CLASS
|
||||
|
||||
KSpaceStyle(msm/omp,MSMOMP)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_MSM_OMP_H
|
||||
#define LMP_MSM_OMP_H
|
||||
|
||||
#include "msm.h"
|
||||
#include "thr_omp.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class MSMOMP : public MSM, public ThrOMP {
|
||||
public:
|
||||
MSMOMP(class LAMMPS *, int, char **);
|
||||
virtual ~MSMOMP () {};
|
||||
|
||||
protected:
|
||||
virtual void direct(int);
|
||||
virtual void compute(int,int);
|
||||
|
||||
private:
|
||||
template <int, int, int> void direct_eval(int);
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,188 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 "pair_born_coul_msm_omp.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairBornCoulMSMOMP::PairBornCoulMSMOMP(LAMMPS *lmp) :
|
||||
PairBornCoulMSM(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
{
|
||||
suffix_flag |= Suffix::OMP;
|
||||
respa_enable = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairBornCoulMSMOMP::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;
|
||||
|
||||
#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 PairBornCoulMSMOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
{
|
||||
int i,j,ii,jj,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double rsq,r2inv,r6inv,r,rexp,forcecoul,forceborn,factor_coul,factor_lj;
|
||||
double egamma,fgamma,prefactor;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
|
||||
const double * const * const x = atom->x;
|
||||
double * const * const f = thr->get_f();
|
||||
const double * const q = atom->q;
|
||||
const int * const type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
double *special_lj = force->special_lj;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
double fxtmp,fytmp,fztmp;
|
||||
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = iifrom; ii < iito; ++ii) {
|
||||
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[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];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
r = sqrt(rsq);
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
|
||||
fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);
|
||||
forcecoul = prefactor * fgamma;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]);
|
||||
forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv
|
||||
+ born3[itype][jtype]*r2inv*r6inv;
|
||||
} else forceborn = 0.0;
|
||||
|
||||
fpair = (forcecoul + factor_lj*forceborn)*r2inv;
|
||||
|
||||
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) {
|
||||
if (rsq < cut_coulsq) {
|
||||
ecoul = prefactor*egamma;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv
|
||||
+ d[itype][jtype]*r6inv*r2inv - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
|
||||
evdwl,ecoul,fpair,delx,dely,delz,thr);
|
||||
}
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairBornCoulMSMOMP::memory_usage()
|
||||
{
|
||||
double bytes = memory_usage_thr();
|
||||
bytes += PairBornCoulMSM::memory_usage();
|
||||
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
/* -*- 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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(born/coul/msm/omp,PairBornCoulMSMOMP)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_BORN_COUL_MSM_OMP_H
|
||||
#define LMP_PAIR_BORN_COUL_MSM_OMP_H
|
||||
|
||||
#include "pair_born_coul_msm.h"
|
||||
#include "thr_omp.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairBornCoulMSMOMP : public PairBornCoulMSM, public ThrOMP {
|
||||
|
||||
public:
|
||||
PairBornCoulMSMOMP(class LAMMPS *);
|
||||
|
||||
virtual void compute(int, int);
|
||||
virtual double memory_usage();
|
||||
|
||||
private:
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||
void eval(int ifrom, int ito, ThrData * const thr);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,188 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 "pair_buck_coul_msm_omp.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairBuckCoulMSMOMP::PairBuckCoulMSMOMP(LAMMPS *lmp) :
|
||||
PairBuckCoulMSM(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
{
|
||||
suffix_flag |= Suffix::OMP;
|
||||
respa_enable = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairBuckCoulMSMOMP::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;
|
||||
|
||||
#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 per thread forces into global force array.
|
||||
reduce_thr(this, eflag, vflag, thr);
|
||||
} // end of omp parallel region
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||
void PairBuckCoulMSMOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
{
|
||||
int i,j,ii,jj,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double rsq,r2inv,r6inv,r,rexp,forcecoul,forcebuck,factor_coul,factor_lj;
|
||||
double egamma,fgamma,prefactor;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
|
||||
const double * const * const x = atom->x;
|
||||
double * const * const f = thr->get_f();
|
||||
const double * const q = atom->q;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
double *special_lj = force->special_lj;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
double fxtmp,fytmp,fztmp;
|
||||
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = iifrom; ii < iito; ++ii) {
|
||||
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[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];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
r = sqrt(rsq);
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
|
||||
fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);
|
||||
forcecoul = prefactor * fgamma;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
rexp = exp(-r*rhoinv[itype][jtype]);
|
||||
forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv;
|
||||
} else forcebuck = 0.0;
|
||||
|
||||
fpair = (forcecoul + factor_lj*forcebuck)*r2inv;
|
||||
|
||||
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) {
|
||||
if (rsq < cut_coulsq) {
|
||||
ecoul = prefactor*egamma;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv -
|
||||
offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
|
||||
evdwl,ecoul,fpair,delx,dely,delz,thr);
|
||||
}
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairBuckCoulMSMOMP::memory_usage()
|
||||
{
|
||||
double bytes = memory_usage_thr();
|
||||
bytes += PairBuckCoulMSM::memory_usage();
|
||||
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
/* -*- 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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(buck/coul/msm/omp,PairBuckCoulMSMOMP)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_BUCK_COUL_MSM_OMP_H
|
||||
#define LMP_PAIR_BUCK_COUL_MSM_OMP_H
|
||||
|
||||
#include "pair_buck_coul_msm.h"
|
||||
#include "thr_omp.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairBuckCoulMSMOMP : public PairBuckCoulMSM, public ThrOMP {
|
||||
|
||||
public:
|
||||
PairBuckCoulMSMOMP(class LAMMPS *);
|
||||
|
||||
virtual void compute(int, int);
|
||||
virtual double memory_usage();
|
||||
|
||||
private:
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||
void eval(int ifrom, int ito, ThrData * const thr);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,191 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 "pair_coul_msm_omp.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulMSMOMP::PairCoulMSMOMP(LAMMPS *lmp) :
|
||||
PairCoulMSM(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
{
|
||||
suffix_flag |= Suffix::OMP;
|
||||
respa_enable = 0;
|
||||
cut_respa = NULL;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulMSMOMP::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;
|
||||
|
||||
#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 PairCoulMSMOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
{
|
||||
int i,j,ii,jj,jnum,itable,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,r2inv,rsq,forcecoul,factor_coul;
|
||||
double egamma,fgamma,prefactor;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
ecoul = 0.0;
|
||||
|
||||
const double * const * const x = atom->x;
|
||||
double * const * const f = thr->get_f();
|
||||
const double * const q = atom->q;
|
||||
const int * const type = atom->type;
|
||||
const int nlocal = atom->nlocal;
|
||||
const double * const special_coul = force->special_coul;
|
||||
const double qqrd2e = force->qqrd2e;
|
||||
double fxtmp,fytmp,fztmp;
|
||||
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = iifrom; ii < iito; ++ii) {
|
||||
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[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];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
r2inv = 1.0/rsq;
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
|
||||
egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
|
||||
fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);
|
||||
forcecoul = prefactor * fgamma;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
}
|
||||
|
||||
fpair = forcecoul * r2inv;
|
||||
|
||||
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) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq)
|
||||
ecoul = prefactor*egamma;
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
ecoul = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
}
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
|
||||
if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
|
||||
0.0,ecoul,fpair,delx,dely,delz,thr);
|
||||
}
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairCoulMSMOMP::memory_usage()
|
||||
{
|
||||
double bytes = memory_usage_thr();
|
||||
bytes += PairCoulMSM::memory_usage();
|
||||
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
/* -*- 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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(coul/msm/omp,PairCoulMSMOMP)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_COUL_MSM_OMP_H
|
||||
#define LMP_PAIR_COUL_MSM_OMP_H
|
||||
|
||||
#include "pair_coul_msm.h"
|
||||
#include "thr_omp.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairCoulMSMOMP : public PairCoulMSM, public ThrOMP {
|
||||
|
||||
public:
|
||||
PairCoulMSMOMP(class LAMMPS *);
|
||||
|
||||
virtual void compute(int, int);
|
||||
virtual double memory_usage();
|
||||
|
||||
private:
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||
void eval(int ifrom, int ito, ThrData * const thr);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,216 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 "pair_lj_charmm_coul_msm_omp.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCharmmCoulMSMOMP::PairLJCharmmCoulMSMOMP(LAMMPS *lmp) :
|
||||
PairLJCharmmCoulMSM(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
{
|
||||
suffix_flag |= Suffix::OMP;
|
||||
respa_enable = 0;
|
||||
cut_respa = NULL;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCharmmCoulMSMOMP::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;
|
||||
|
||||
#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 PairLJCharmmCoulMSMOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
{
|
||||
|
||||
const double * const * const x = atom->x;
|
||||
double * const * const f = thr->get_f();
|
||||
const double * const q = atom->q;
|
||||
const int * const type = atom->type;
|
||||
const double * const special_coul = force->special_coul;
|
||||
const double * const special_lj = force->special_lj;
|
||||
const double qqrd2e = force->qqrd2e;
|
||||
const double inv_denom_lj = 1.0/denom_lj;
|
||||
|
||||
const int * const ilist = list->ilist;
|
||||
const int * const numneigh = list->numneigh;
|
||||
const int * const * const firstneigh = list->firstneigh;
|
||||
const int nlocal = atom->nlocal;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (int ii = iifrom; ii < iito; ++ii) {
|
||||
|
||||
const int i = ilist[ii];
|
||||
const int itype = type[i];
|
||||
const double qtmp = q[i];
|
||||
const double xtmp = x[i][0];
|
||||
const double ytmp = x[i][1];
|
||||
const double ztmp = x[i][2];
|
||||
double fxtmp,fytmp,fztmp;
|
||||
fxtmp=fytmp=fztmp=0.0;
|
||||
|
||||
const int * const jlist = firstneigh[i];
|
||||
const int jnum = numneigh[i];
|
||||
|
||||
for (int jj = 0; jj < jnum; jj++) {
|
||||
double forcecoul, forcelj, evdwl, ecoul;
|
||||
forcecoul = forcelj = evdwl = ecoul = 0.0;
|
||||
|
||||
const int sbindex = sbmask(jlist[jj]);
|
||||
const int j = jlist[jj] & NEIGHMASK;
|
||||
|
||||
const double delx = xtmp - x[j][0];
|
||||
const double dely = ytmp - x[j][1];
|
||||
const double delz = ztmp - x[j][2];
|
||||
const double rsq = delx*delx + dely*dely + delz*delz;
|
||||
const int jtype = type[j];
|
||||
|
||||
if (rsq < cut_bothsq) {
|
||||
const double r2inv = 1.0/rsq;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
|
||||
const double r = sqrt(rsq);
|
||||
const double prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
const double egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
|
||||
const double fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);
|
||||
forcecoul = prefactor * (fgamma - 1.0);
|
||||
|
||||
if (EFLAG) ecoul = prefactor*egamma;
|
||||
if (sbindex) {
|
||||
const double adjust = (1.0-special_coul[sbindex])*prefactor;
|
||||
forcecoul -= adjust;
|
||||
if (EFLAG) ecoul -= adjust;
|
||||
}
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
|
||||
const double fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
const double table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (EFLAG) ecoul = qtmp*q[j] * (etable[itable] + fraction*detable[itable]);
|
||||
if (sbindex) {
|
||||
const double table2 = ctable[itable] + fraction*dctable[itable];
|
||||
const double prefactor = qtmp*q[j] * table2;
|
||||
const double adjust = (1.0-special_coul[sbindex])*prefactor;
|
||||
forcecoul -= adjust;
|
||||
if (EFLAG) ecoul -= adjust;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (rsq < cut_ljsq) {
|
||||
const double r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
if (EFLAG) evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
|
||||
|
||||
if (rsq > cut_lj_innersq) {
|
||||
const double drsq = cut_ljsq - rsq;
|
||||
const double cut2 = (rsq - cut_lj_innersq) * drsq;
|
||||
const double switch1 = drsq * (drsq*drsq + 3.0*cut2) * inv_denom_lj;
|
||||
const double switch2 = 12.0*rsq * cut2 * inv_denom_lj;
|
||||
if (EFLAG) {
|
||||
forcelj = forcelj*switch1 + evdwl*switch2;
|
||||
evdwl *= switch1;
|
||||
} else {
|
||||
const double philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
|
||||
forcelj = forcelj*switch1 + philj*switch2;
|
||||
}
|
||||
}
|
||||
|
||||
if (sbindex) {
|
||||
const double factor_lj = special_lj[sbindex];
|
||||
forcelj *= factor_lj;
|
||||
if (EFLAG) evdwl *= factor_lj;
|
||||
}
|
||||
|
||||
}
|
||||
const double fpair = (forcecoul + forcelj) * r2inv;
|
||||
|
||||
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 (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR,
|
||||
evdwl,ecoul,fpair,delx,dely,delz,thr);
|
||||
}
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCharmmCoulMSMOMP::memory_usage()
|
||||
{
|
||||
double bytes = memory_usage_thr();
|
||||
bytes += PairLJCharmmCoulMSM::memory_usage();
|
||||
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
/* -*- 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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(lj/charmm/coul/msm/omp,PairLJCharmmCoulMSMOMP)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_LJ_CHARMM_COUL_MSM_OMP_H
|
||||
#define LMP_PAIR_LJ_CHARMM_COUL_MSM_OMP_H
|
||||
|
||||
#include "pair_lj_charmm_coul_msm.h"
|
||||
#include "thr_omp.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairLJCharmmCoulMSMOMP : public PairLJCharmmCoulMSM, public ThrOMP {
|
||||
|
||||
public:
|
||||
PairLJCharmmCoulMSMOMP(class LAMMPS *);
|
||||
|
||||
virtual void compute(int, int);
|
||||
virtual double memory_usage();
|
||||
|
||||
private:
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||
void eval(int ifrom, int ito, ThrData * const thr);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,210 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 "pair_lj_cut_coul_msm_omp.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include "suffix.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulMSMOMP::PairLJCutCoulMSMOMP(LAMMPS *lmp) :
|
||||
PairLJCutCoulMSM(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
{
|
||||
suffix_flag |= Suffix::OMP;
|
||||
respa_enable = 0;
|
||||
cut_respa = NULL;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulMSMOMP::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;
|
||||
|
||||
#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 PairLJCutCoulMSMOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
{
|
||||
int i,j,ii,jj,jnum,itype,jtype,itable;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double egamma,fgamma,prefactor;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
|
||||
const double * const * const x = atom->x;
|
||||
double * const * const f = thr->get_f();
|
||||
const double * const q = atom->q;
|
||||
const int * const type = atom->type;
|
||||
const int nlocal = atom->nlocal;
|
||||
const double * const special_coul = force->special_coul;
|
||||
const double * const special_lj = force->special_lj;
|
||||
const double qqrd2e = force->qqrd2e;
|
||||
double fxtmp,fytmp,fztmp;
|
||||
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = iifrom; ii < iito; ++ii) {
|
||||
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[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];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
|
||||
fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);
|
||||
forcecoul = prefactor * fgamma;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = qtmp*q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
}
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
forcelj *= factor_lj;
|
||||
} else forcelj = 0.0;
|
||||
|
||||
fpair = (forcecoul + forcelj) * r2inv;
|
||||
|
||||
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) {
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq)
|
||||
ecoul = prefactor*egamma;
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
ecoul = qtmp*q[j] * table;
|
||||
}
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
|
||||
evdwl,ecoul,fpair,delx,dely,delz,thr);
|
||||
}
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutCoulMSMOMP::memory_usage()
|
||||
{
|
||||
double bytes = memory_usage_thr();
|
||||
bytes += PairLJCutCoulMSM::memory_usage();
|
||||
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
/* -*- 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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(lj/cut/coul/msm/omp,PairLJCutCoulMSMOMP)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_LJ_CUT_COUL_MSM_OMP_H
|
||||
#define LMP_PAIR_LJ_CUT_COUL_MSM_OMP_H
|
||||
|
||||
#include "pair_lj_cut_coul_msm.h"
|
||||
#include "thr_omp.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairLJCutCoulMSMOMP : public PairLJCutCoulMSM, public ThrOMP {
|
||||
|
||||
public:
|
||||
PairLJCutCoulMSMOMP(class LAMMPS *);
|
||||
|
||||
virtual void compute(int, int);
|
||||
virtual double memory_usage();
|
||||
|
||||
private:
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||
void eval(int ifrom, int ito, ThrData * const thr);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
Loading…
Reference in New Issue