From 7ec4a5818d50a91453a15168ac6a7a1ae35d6304 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Sun, 18 Feb 2018 16:38:11 +0100
Subject: [PATCH] add USER-OMP versions of lj/cut/coul/wolf
---
doc/src/pair_lj.txt | 1 +
src/USER-OMP/pair_lj_cut_coul_wolf_omp.cpp | 204 +++++++++++++++++++++
src/USER-OMP/pair_lj_cut_coul_wolf_omp.h | 48 +++++
3 files changed, 253 insertions(+)
create mode 100644 src/USER-OMP/pair_lj_cut_coul_wolf_omp.cpp
create mode 100644 src/USER-OMP/pair_lj_cut_coul_wolf_omp.h
diff --git a/doc/src/pair_lj.txt b/doc/src/pair_lj.txt
index 1f3805643a..e297d479bc 100644
--- a/doc/src/pair_lj.txt
+++ b/doc/src/pair_lj.txt
@@ -33,6 +33,7 @@ pair_style lj/cut/coul/msm command :h3
pair_style lj/cut/coul/msm/gpu command :h3
pair_style lj/cut/coul/msm/omp command :h3
pair_style lj/cut/coul/wolf command :h3
+pair_style lj/cut/coul/wolf/omp command :h3
pair_style lj/cut/tip4p/cut command :h3
pair_style lj/cut/tip4p/cut/omp command :h3
pair_style lj/cut/tip4p/long command :h3
diff --git a/src/USER-OMP/pair_lj_cut_coul_wolf_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_wolf_omp.cpp
new file mode 100644
index 0000000000..f5dc03e368
--- /dev/null
+++ b/src/USER-OMP/pair_lj_cut_coul_wolf_omp.cpp
@@ -0,0 +1,204 @@
+/* ----------------------------------------------------------------------
+ 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_wolf_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+
+#include "suffix.h"
+#include "math_const.h"
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+PairLJCutCoulWolfOMP::PairLJCutCoulWolfOMP(LAMMPS *lmp) :
+ PairLJCutCoulWolf(lmp), ThrOMP(lmp, THR_PAIR)
+{
+ suffix_flag |= Suffix::OMP;
+ respa_enable = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJCutCoulWolfOMP::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);
+ 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 PairLJCutCoulWolfOMP::eval(int iifrom, int iito, ThrData * const thr)
+{
+ int i,j,ii,jj,jnum,itype,jtype;
+ double qitmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
+ double r,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
+ double prefactor,erfcc,erfcd,v_sh,dvdrr,e_self,qisq;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = ecoul = 0.0;
+
+ const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
+ dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
+ const double * _noalias const q = atom->q;
+ const int * _noalias const type = atom->type;
+ const int nlocal = atom->nlocal;
+ const double * _noalias const special_coul = force->special_coul;
+ const double * _noalias const special_lj = force->special_lj;
+ const double qqrd2e = force->qqrd2e;
+ double fxtmp,fytmp,fztmp;
+
+ // self and shifted Coulombic energy
+
+ e_self = v_sh = 0.0;
+ const double e_shift = erfc(alf*cut_coul)/cut_coul;
+ const double f_shift = -(e_shift+ 2.0*alf/MY_PIS
+ * exp(-alf*alf*cut_coul*cut_coul)) / cut_coul;
+
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = iifrom; ii < iito; ++ii) {
+
+ i = ilist[ii];
+ qitmp = q[i];
+ xtmp = x[i].x;
+ ytmp = x[i].y;
+ ztmp = x[i].z;
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+ fxtmp=fytmp=fztmp=0.0;
+
+ if (EFLAG) {
+ e_self = -(e_shift/2.0 + alf/MY_PIS) * qitmp*qitmp*qqrd2e;
+ ev_tally_thr(this,i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0,thr);
+ }
+
+ 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].x;
+ dely = ytmp - x[j].y;
+ delz = ztmp - x[j].z;
+ rsq = delx*delx + dely*dely + delz*delz;
+ jtype = type[j];
+
+ if (rsq < cutsq[itype][jtype]) {
+ r2inv = 1.0/rsq;
+
+ if (rsq < cut_coulsq) {
+ r = sqrt(rsq);
+ prefactor = qqrd2e*qitmp*q[j]/r;
+ erfcc = erfc(alf*r);
+ erfcd = exp(-alf*alf*r*r);
+ v_sh = (erfcc - e_shift*r) * prefactor;
+ dvdrr = (erfcc/rsq + 2.0*alf/MY_PIS * erfcd/r) + f_shift;
+ forcecoul = dvdrr*rsq*prefactor;
+ 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;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ } else forcelj = 0.0;
+
+ fpair = (forcecoul + factor_lj*forcelj) * r2inv;
+
+ fxtmp += delx*fpair;
+ fytmp += dely*fpair;
+ fztmp += delz*fpair;
+ if (NEWTON_PAIR || j < nlocal) {
+ f[j].x -= delx*fpair;
+ f[j].y -= dely*fpair;
+ f[j].z -= delz*fpair;
+ }
+
+ if (EFLAG) {
+ 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 (rsq < cut_coulsq) {
+ ecoul = v_sh;
+ if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+ } else ecoul = 0.0;
+ }
+
+ if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
+ evdwl,ecoul,fpair,delx,dely,delz,thr);
+ }
+ }
+ f[i].x += fxtmp;
+ f[i].y += fytmp;
+ f[i].z += fztmp;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLJCutCoulWolfOMP::memory_usage()
+{
+ double bytes = memory_usage_thr();
+ bytes += PairLJCutCoulWolf::memory_usage();
+
+ return bytes;
+}
diff --git a/src/USER-OMP/pair_lj_cut_coul_wolf_omp.h b/src/USER-OMP/pair_lj_cut_coul_wolf_omp.h
new file mode 100644
index 0000000000..4e56808c7b
--- /dev/null
+++ b/src/USER-OMP/pair_lj_cut_coul_wolf_omp.h
@@ -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/wolf/omp,PairLJCutCoulWolfOMP)
+
+#else
+
+#ifndef LMP_PAIR_LJ_CUT_COUL_WOLF_OMP_H
+#define LMP_PAIR_LJ_CUT_COUL_WOLF_OMP_H
+
+#include "pair_lj_cut_coul_wolf.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairLJCutCoulWolfOMP : public PairLJCutCoulWolf, public ThrOMP {
+
+ public:
+ PairLJCutCoulWolfOMP(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
--
GitLab