From 27809e8a244164e27a3ea625c65541f645eac696 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 2 Oct 2012 20:08:19 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8865 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- .../pair_lj_charmm_coul_long_proxy_omp.cpp | 259 +++++++++ .../pair_lj_charmm_coul_long_proxy_omp.h | 55 ++ .../pair_lj_class2_coul_long_proxy_omp.cpp | 226 ++++++++ .../pair_lj_class2_coul_long_proxy_omp.h | 55 ++ .../pair_lj_cut_coul_long_proxy_omp.cpp | 244 ++++++++ .../pair_lj_cut_coul_long_proxy_omp.h | 54 ++ .../pair_lj_cut_coul_long_proxy_tip4p_omp.cpp | 540 ++++++++++++++++++ .../pair_lj_cut_coul_long_proxy_tip4p_omp.h | 57 ++ 8 files changed, 1490 insertions(+) create mode 100644 src/USER-OMP/pair_lj_charmm_coul_long_proxy_omp.cpp create mode 100644 src/USER-OMP/pair_lj_charmm_coul_long_proxy_omp.h create mode 100644 src/USER-OMP/pair_lj_class2_coul_long_proxy_omp.cpp create mode 100644 src/USER-OMP/pair_lj_class2_coul_long_proxy_omp.h create mode 100644 src/USER-OMP/pair_lj_cut_coul_long_proxy_omp.cpp create mode 100644 src/USER-OMP/pair_lj_cut_coul_long_proxy_omp.h create mode 100644 src/USER-OMP/pair_lj_cut_coul_long_proxy_tip4p_omp.cpp create mode 100644 src/USER-OMP/pair_lj_cut_coul_long_proxy_tip4p_omp.h diff --git a/src/USER-OMP/pair_lj_charmm_coul_long_proxy_omp.cpp b/src/USER-OMP/pair_lj_charmm_coul_long_proxy_omp.cpp new file mode 100644 index 0000000000..afdfe05571 --- /dev/null +++ b/src/USER-OMP/pair_lj_charmm_coul_long_proxy_omp.cpp @@ -0,0 +1,259 @@ +/* ---------------------------------------------------------------------- + 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_long_proxy_omp.h" +#include "pppm_proxy.h" +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "update.h" + +#include + +#include "suffix.h" +using namespace LAMMPS_NS; + +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairLJCharmmCoulLongProxyOMP::PairLJCharmmCoulLongProxyOMP(LAMMPS *lmp) : + PairLJCharmmCoulLong(lmp), ThrOMP(lmp, THR_PAIR|THR_PROXY) +{ + proxyflag = 1; + suffix_flag |= Suffix::OMP; + respa_enable = 0; + nproxy=1; + + kspace = NULL; +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCharmmCoulLongProxyOMP::init_style() +{ + if (comm->nthreads < 2) + error->all(FLERR,"need at least two threads per MPI task for this pair style"); + + kspace = static_cast(force->kspace); + + PairLJCharmmCoulLong::init_style(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCharmmCoulLongProxyOMP::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, nproxy); + ThrData *thr = fix->get_thr(tid); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr); + + // thread id 0 runs pppm, the rest the pair style + if (tid < nproxy) { + kspace->compute_proxy(eflag,vflag); + } else { + 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); + } + } + + sync_threads(); + reduce_thr(this, eflag, vflag, thr, nproxy); + } // end of omp parallel region +} + +/* ---------------------------------------------------------------------- */ + +template +void PairLJCharmmCoulLongProxyOMP::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 grij,expm2,prefactor,t,erfc; + double philj,switch1,switch2; + 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); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = qqrd2e * qtmp*q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + 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) { + r6inv = r2inv*r2inv*r2inv; + jtype = type[j]; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq > cut_lj_innersq) { + switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + switch2 = 12.0*rsq * (cut_ljsq-rsq) * + (rsq-cut_lj_innersq) / denom_lj; + philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); + forcelj = forcelj*switch1 + philj*switch2; + } + 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*erfc; + 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) { + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + if (rsq > cut_lj_innersq) { + switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + evdwl *= switch1; + } + 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 PairLJCharmmCoulLongProxyOMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += PairLJCharmmCoulLong::memory_usage(); + + return bytes; +} diff --git a/src/USER-OMP/pair_lj_charmm_coul_long_proxy_omp.h b/src/USER-OMP/pair_lj_charmm_coul_long_proxy_omp.h new file mode 100644 index 0000000000..b7fba4f6cb --- /dev/null +++ b/src/USER-OMP/pair_lj_charmm_coul_long_proxy_omp.h @@ -0,0 +1,55 @@ +/* -*- 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/long/proxy/omp,PairLJCharmmCoulLongProxyOMP) + +#else + +#ifndef LMP_PAIR_LJ_CHARMM_COUL_LONG_PROXY_OMP_H +#define LMP_PAIR_LJ_CHARMM_COUL_LONG_PROXY_OMP_H + +#include "pair_lj_charmm_coul_long.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairLJCharmmCoulLongProxyOMP : + public PairLJCharmmCoulLong, public ThrOMP { + + public: + PairLJCharmmCoulLongProxyOMP(class LAMMPS *); + virtual ~PairLJCharmmCoulLongProxyOMP() {}; + + virtual void init_style(); + + virtual void compute(int, int); + virtual double memory_usage(); + + private: + template + void eval(int ifrom, int ito, ThrData * const thr); + + class PPPMProxy *kspace; + int nproxy; +}; + +} + +#endif +#endif diff --git a/src/USER-OMP/pair_lj_class2_coul_long_proxy_omp.cpp b/src/USER-OMP/pair_lj_class2_coul_long_proxy_omp.cpp new file mode 100644 index 0000000000..c49e4abe1f --- /dev/null +++ b/src/USER-OMP/pair_lj_class2_coul_long_proxy_omp.cpp @@ -0,0 +1,226 @@ +/* ---------------------------------------------------------------------- + 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_class2_coul_long_proxy_omp.h" +#include "pppm_proxy.h" +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "update.h" + +#include + +#include "suffix.h" +using namespace LAMMPS_NS; + +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairLJClass2CoulLongProxyOMP::PairLJClass2CoulLongProxyOMP(LAMMPS *lmp) : + PairLJClass2CoulLong(lmp), ThrOMP(lmp, THR_PAIR|THR_PROXY) +{ + proxyflag = 1; + suffix_flag |= Suffix::OMP; + respa_enable = 0; + nproxy=1; + + kspace = NULL; +} + +/* ---------------------------------------------------------------------- */ + +void PairLJClass2CoulLongProxyOMP::init_style() +{ + if (comm->nthreads < 2) + error->all(FLERR,"need at least two threads per MPI task for this pair style"); + + kspace = static_cast(force->kspace); + + PairLJClass2CoulLong::init_style(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJClass2CoulLongProxyOMP::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, nproxy); + ThrData *thr = fix->get_thr(tid); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr); + + // thread id 0 runs pppm, the rest the pair style + if (tid < nproxy) { + kspace->compute_proxy(eflag,vflag); + } else { + 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); + } + } + + sync_threads(); + reduce_thr(this, eflag, vflag, thr, nproxy); + } // end of omp parallel region +} + +/* ---------------------------------------------------------------------- */ + +template +void PairLJClass2CoulLongProxyOMP::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 r,rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double grij,expm2,prefactor,t,erfc; + 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) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = qqrd2e * qtmp*q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + rinv = sqrt(r2inv); + r3inv = r2inv*rinv; + r6inv = r3inv*r3inv; + forcelj = r6inv * (lj1[itype][jtype]*r3inv - 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) { + ecoul = prefactor*erfc; + 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]*r3inv-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 PairLJClass2CoulLongProxyOMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += PairLJClass2CoulLong::memory_usage(); + + return bytes; +} diff --git a/src/USER-OMP/pair_lj_class2_coul_long_proxy_omp.h b/src/USER-OMP/pair_lj_class2_coul_long_proxy_omp.h new file mode 100644 index 0000000000..df2aa72043 --- /dev/null +++ b/src/USER-OMP/pair_lj_class2_coul_long_proxy_omp.h @@ -0,0 +1,55 @@ +/* -*- 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/class2/coul/long/proxy/omp,PairLJClass2CoulLongProxyOMP) + +#else + +#ifndef LMP_PAIR_LJ_CLASS2_COUL_LONG_PROXY_OMP_H +#define LMP_PAIR_LJ_CLASS2_COUL_LONG_PROXY_OMP_H + +#include "pair_lj_class2_coul_long.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairLJClass2CoulLongProxyOMP : + public PairLJClass2CoulLong, public ThrOMP { + + public: + PairLJClass2CoulLongProxyOMP(class LAMMPS *); + virtual ~PairLJClass2CoulLongProxyOMP() {}; + + virtual void init_style(); + + virtual void compute(int, int); + virtual double memory_usage(); + + private: + template + void eval(int ifrom, int ito, ThrData * const thr); + + class PPPMProxy *kspace; + int nproxy; +}; + +} + +#endif +#endif diff --git a/src/USER-OMP/pair_lj_cut_coul_long_proxy_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_long_proxy_omp.cpp new file mode 100644 index 0000000000..b1e7ba1839 --- /dev/null +++ b/src/USER-OMP/pair_lj_cut_coul_long_proxy_omp.cpp @@ -0,0 +1,244 @@ +/* ---------------------------------------------------------------------- + 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_long_proxy_omp.h" +#include "pppm_proxy.h" +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "update.h" + +#include + +#include "suffix.h" +using namespace LAMMPS_NS; + +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairLJCutCoulLongProxyOMP::PairLJCutCoulLongProxyOMP(LAMMPS *lmp) : + PairLJCutCoulLong(lmp), ThrOMP(lmp, THR_PAIR|THR_PROXY) +{ + proxyflag = 1; + suffix_flag |= Suffix::OMP; + respa_enable = 0; + nproxy=1; + + kspace = NULL; +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLongProxyOMP::init_style() +{ + if (comm->nthreads < 2) + error->all(FLERR,"need at least two threads per MPI task for this pair style"); + + kspace = static_cast(force->kspace); + + PairLJCutCoulLong::init_style(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLongProxyOMP::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, nproxy); + ThrData *thr = fix->get_thr(tid); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr); + + // thread id 0 runs pppm, the rest the pair style + if (tid < nproxy) { + kspace->compute_proxy(eflag,vflag); + } else { + 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); + } + } + + sync_threads(); + reduce_thr(this, eflag, vflag, thr, nproxy); + } // end of omp parallel region +} +/* ---------------------------------------------------------------------- */ + +template +void PairLJCutCoulLongProxyOMP::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 grij,expm2,prefactor,t,erfc; + 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); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = qqrd2e * qtmp*q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + 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*erfc; + 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 PairLJCutCoulLongProxyOMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += PairLJCutCoulLong::memory_usage(); + + return bytes; +} diff --git a/src/USER-OMP/pair_lj_cut_coul_long_proxy_omp.h b/src/USER-OMP/pair_lj_cut_coul_long_proxy_omp.h new file mode 100644 index 0000000000..7d653b3d2e --- /dev/null +++ b/src/USER-OMP/pair_lj_cut_coul_long_proxy_omp.h @@ -0,0 +1,54 @@ +/* -*- 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/long/proxy/omp,PairLJCutCoulLongProxyOMP) + +#else + +#ifndef LMP_PAIR_LJ_CUT_COUL_LONG_PROXY_OMP_H +#define LMP_PAIR_LJ_CUT_COUL_LONG_PROXY_OMP_H + +#include "pair_lj_cut_coul_long.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairLJCutCoulLongProxyOMP : public PairLJCutCoulLong, public ThrOMP { + + public: + PairLJCutCoulLongProxyOMP(class LAMMPS *); + virtual ~PairLJCutCoulLongProxyOMP() {}; + + virtual void init_style(); + + virtual void compute(int, int); + virtual double memory_usage(); + + private: + template + void eval(int ifrom, int ito, ThrData * const thr); + + class PPPMProxy *kspace; + int nproxy; +}; + +} + +#endif +#endif diff --git a/src/USER-OMP/pair_lj_cut_coul_long_proxy_tip4p_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_long_proxy_tip4p_omp.cpp new file mode 100644 index 0000000000..64173d72d6 --- /dev/null +++ b/src/USER-OMP/pair_lj_cut_coul_long_proxy_tip4p_omp.cpp @@ -0,0 +1,540 @@ +/* ---------------------------------------------------------------------- + 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_long_proxy_tip4p_omp.h" +#include "pppm_tip4p_proxy.h" +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "update.h" + +#include + +#include "suffix.h" +using namespace LAMMPS_NS; + +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairLJCutCoulLongProxyTIP4POMP::PairLJCutCoulLongProxyTIP4POMP(LAMMPS *lmp) : + PairLJCutCoulLongTIP4P(lmp), ThrOMP(lmp, THR_PAIR|THR_PROXY) +{ + proxyflag = 1; + suffix_flag |= Suffix::OMP; + respa_enable = 0; + nproxy=1; + + kspace = NULL; + + // TIP4P cannot compute virial as F dot r + // due to finding bonded H atoms which are not near O atom + + no_virial_fdotr_compute = 1; +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLongProxyTIP4POMP::init_style() +{ + if (comm->nthreads < 2) + error->all(FLERR,"need at least two threads per MPI task for this pair style"); + + kspace = static_cast(force->kspace); + + PairLJCutCoulLongTIP4P::init_style(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutCoulLongProxyTIP4POMP::compute(int eflag, int vflag) +{ + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + const int nlocal = atom->nlocal; + const int nall = nlocal + atom->nghost; + + // reallocate hneigh & newsite if necessary + // initialize hneigh[0] to -1 on steps when reneighboring occurred + // initialize hneigh[2] to 0 every step + + if (atom->nmax > nmax) { + nmax = atom->nmax; + memory->destroy(hneigh); + memory->create(hneigh,nmax,3,"pair:hneigh"); + memory->destroy(newsite); + memory->create(newsite,nmax,3,"pair:newsite"); + } + + // XXX: this could be threaded, too. + int i; + if (neighbor->ago == 0) + for (i = 0; i < nall; i++) hneigh[i][0] = -1; + for (i = 0; i < nall; i++) hneigh[i][2] = 0; + + 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, nproxy); + ThrData *thr = fix->get_thr(tid); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr); + + // thread id 0 runs pppm, the rest the pair style + if (tid < nproxy) { + kspace->compute_proxy(eflag,vflag); + } else { + if (!ncoultablebits) { + if (evflag) { + if (eflag) { + if (vflag) eval<1,1,1,1>(ifrom, ito, thr); + else eval<1,1,1,0>(ifrom, ito, thr); + } else { + if (vflag) eval<1,1,0,1>(ifrom, ito, thr); + else eval<1,1,0,0>(ifrom, ito, thr); + } + } else eval<1,0,0,0>(ifrom, ito, thr); + } else { + if (evflag) { + if (eflag) { + if (vflag) eval<0,1,1,1>(ifrom, ito, thr); + else eval<0,1,1,0>(ifrom, ito, thr); + } else { + if (vflag) eval<0,1,0,1>(ifrom, ito, thr); + else eval<0,1,0,0>(ifrom, ito, thr); + } + } else eval<0,0,0,0>(ifrom, ito, thr); + } + } + + sync_threads(); + reduce_thr(this, eflag, vflag, thr, nproxy); + } // end of omp parallel region +} + +/* ---------------------------------------------------------------------- */ + +template +void PairLJCutCoulLongProxyTIP4POMP::eval(int iifrom, int iito, ThrData * const thr) +{ + int i,j,ii,jj,jnum,itype,jtype,itable; + int n,vlist[6]; + int iH1,iH2,jH1,jH2; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul; + double fraction,table; + double delxOM,delyOM,delzOM; + double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce; + double factor_coul,factor_lj; + double grij,expm2,prefactor,t,erfc,ddotf; + double v[6],xH1[3],xH2[3]; + double fdx,fdy,fdz,f1x,f1y,f1z,fOx,fOy,fOz,fHx,fHy,fHz; + const double *x1,*x2; + 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 int tid = thr->get_tid(); + const double * const special_coul = force->special_coul; + const double * const special_lj = force->special_lj; + const double qqrd2e = force->qqrd2e; + const double cut_coulsqplus = (cut_coul+2.0*qdist) * (cut_coul+2.0*qdist); + + 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]; + + // if atom I = water O, set x1 = offset charge site + // else x1 = x of atom I + // NOTE: to make this part thread safe, we need to + // make sure that the hneigh[][] entries only get + // updated, when all data is in place. worst case, + // some calculation is repeated, but since the results + // will be the same, there is no race condition. + if (itype == typeO) { + if (hneigh[i][0] < 0) { + iH1 = atom->map(atom->tag[i] + 1); + iH2 = atom->map(atom->tag[i] + 2); + if (iH1 == -1 || iH2 == -1) + error->one(FLERR,"TIP4P hydrogen is missing"); + if (atom->type[iH1] != typeH || atom->type[iH2] != typeH) + error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); + compute_newsite_thr(x[i],x[iH1],x[iH2],newsite[i]); + hneigh[i][2] = 1; + hneigh[i][1] = iH2; + hneigh[i][0] = iH1; + } else { + iH1 = hneigh[i][0]; + iH2 = hneigh[i][1]; + if (hneigh[i][2] == 0) { + compute_newsite_thr(x[i],x[iH1],x[iH2],newsite[i]); + hneigh[i][2] = 1; + } + } + x1 = newsite[i]; + } else x1 = x[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]; + + // LJ interaction based on true rsq + + if (rsq < cut_ljsq[itype][jtype]) { + r2inv = 1.0/rsq; + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + forcelj *= factor_lj * r2inv; + + fxtmp += delx*forcelj; + fytmp += dely*forcelj; + fztmp += delz*forcelj; + f[j][0] -= delx*forcelj; + f[j][1] -= dely*forcelj; + f[j][2] -= delz*forcelj; + + if (EFLAG) { + 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 = */ 1, + evdwl,0.0,forcelj,delx,dely,delz,thr); + } + + // adjust rsq and delxyz for off-site O charge(s) if necessary + // but only if they are within reach + // NOTE: to make this part thread safe, we need to + // make sure that the hneigh[][] entries only get + // updated, when all data is in place. worst case, + // some calculation is repeated, but since the results + // will be the same, there is no race condition. + if (rsq < cut_coulsqplus) { + if (itype == typeO || jtype == typeO) { + + // if atom J = water O, set x2 = offset charge site + // else x2 = x of atom J + + if (jtype == typeO) { + if (hneigh[j][0] < 0) { + jH1 = atom->map(atom->tag[j] + 1); + jH2 = atom->map(atom->tag[j] + 2); + if (jH1 == -1 || jH2 == -1) + error->one(FLERR,"TIP4P hydrogen is missing"); + if (atom->type[jH1] != typeH || atom->type[jH2] != typeH) + error->one(FLERR,"TIP4P hydrogen has incorrect atom type"); + compute_newsite_thr(x[j],x[jH1],x[jH2],newsite[j]); + hneigh[j][2] = 1; + hneigh[j][1] = jH2; + hneigh[j][0] = jH1; + } else { + jH1 = hneigh[j][0]; + jH2 = hneigh[j][1]; + if (hneigh[j][2] == 0) { + compute_newsite_thr(x[j],x[jH1],x[jH2],newsite[j]); + hneigh[j][2] = 1; + } + } + x2 = newsite[j]; + } else x2 = x[j]; + + delx = x1[0] - x2[0]; + dely = x1[1] - x2[1]; + delz = x1[2] - x2[2]; + rsq = delx*delx + dely*dely + delz*delz; + } + + // Coulombic interaction based on modified rsq + + if (rsq < cut_coulsq) { + r2inv = 1 / rsq; + if (CTABLE || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = qqrd2e * qtmp*q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + 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; + } + } + + cforce = forcecoul * r2inv; + + // if i,j are not O atoms, force is applied directly + // if i or j are O atoms, force is on fictitious atom & partitioned + // force partitioning due to Feenstra, J Comp Chem, 20, 786 (1999) + // f_f = fictitious force, fO = f_f (1 - 2 alpha), fH = alpha f_f + // preserves total force and torque on water molecule + // virial = sum(r x F) where each water's atoms are near xi and xj + // vlist stores 2,4,6 atoms whose forces contribute to virial + + n = 0; + + if (itype != typeO) { + fxtmp += delx * cforce; + fytmp += dely * cforce; + fztmp += delz * cforce; + + if (VFLAG) { + v[0] = x[i][0] * delx * cforce; + v[1] = x[i][1] * dely * cforce; + v[2] = x[i][2] * delz * cforce; + v[3] = x[i][0] * dely * cforce; + v[4] = x[i][0] * delz * cforce; + v[5] = x[i][1] * delz * cforce; + vlist[n++] = i; + } + + } else { + + fdx = delx*cforce; + fdy = dely*cforce; + fdz = delz*cforce; + + delxOM = x[i][0] - x1[0]; + delyOM = x[i][1] - x1[1]; + delzOM = x[i][2] - x1[2]; + + ddotf = (delxOM * fdx + delyOM * fdy + delzOM * fdz) / + (qdist*qdist); + + f1x = alpha * (fdx - ddotf * delxOM); + f1y = alpha * (fdy - ddotf * delyOM); + f1z = alpha * (fdz - ddotf * delzOM); + + fOx = fdx - f1x; + fOy = fdy - f1y; + fOz = fdz - f1z; + + fHx = 0.5 * f1x; + fHy = 0.5 * f1y; + fHz = 0.5 * f1z; + + fxtmp += fOx; + fytmp += fOy; + fztmp += fOz; + + f[iH1][0] += fHx; + f[iH1][1] += fHy; + f[iH1][2] += fHz; + + f[iH2][0] += fHx; + f[iH2][1] += fHy; + f[iH2][2] += fHz; + + if (VFLAG) { + domain->closest_image(x[i],x[iH1],xH1); + domain->closest_image(x[i],x[iH2],xH2); + + v[0] = x[i][0]*fOx + xH1[0]*fHx + xH2[0]*fHx; + v[1] = x[i][1]*fOy + xH1[1]*fHy + xH2[1]*fHy; + v[2] = x[i][2]*fOz + xH1[2]*fHz + xH2[2]*fHz; + v[3] = x[i][0]*fOy + xH1[0]*fHy + xH2[0]*fHy; + v[4] = x[i][0]*fOz + xH1[0]*fHz + xH2[0]*fHz; + v[5] = x[i][1]*fOz + xH1[1]*fHz + xH2[1]*fHz; + + vlist[n++] = i; + vlist[n++] = iH1; + vlist[n++] = iH2; + } + } + + if (jtype != typeO) { + f[j][0] -= delx * cforce; + f[j][1] -= dely * cforce; + f[j][2] -= delz * cforce; + + if (VFLAG) { + v[0] -= x[j][0] * delx * cforce; + v[1] -= x[j][1] * dely * cforce; + v[2] -= x[j][2] * delz * cforce; + v[3] -= x[j][0] * dely * cforce; + v[4] -= x[j][0] * delz * cforce; + v[5] -= x[j][1] * delz * cforce; + vlist[n++] = j; + } + + } else { + + fdx = -delx*cforce; + fdy = -dely*cforce; + fdz = -delz*cforce; + + delxOM = x[j][0] - x2[0]; + delyOM = x[j][1] - x2[1]; + delzOM = x[j][2] - x2[2]; + + ddotf = (delxOM * fdx + delyOM * fdy + delzOM * fdz) / + (qdist*qdist); + + f1x = alpha * (fdx - ddotf * delxOM); + f1y = alpha * (fdy - ddotf * delyOM); + f1z = alpha * (fdz - ddotf * delzOM); + + fOx = fdx - f1x; + fOy = fdy - f1y; + fOz = fdz - f1z; + + fHx = 0.5 * f1x; + fHy = 0.5 * f1y; + fHz = 0.5 * f1z; + + f[j][0] += fOx; + f[j][1] += fOy; + f[j][2] += fOz; + + f[jH1][0] += fHx; + f[jH1][1] += fHy; + f[jH1][2] += fHz; + + f[jH2][0] += fHx; + f[jH2][1] += fHy; + f[jH2][2] += fHz; + + if (VFLAG) { + domain->closest_image(x[j],x[jH1],xH1); + domain->closest_image(x[j],x[jH2],xH2); + + v[0] += x[j][0]*fOx + xH1[0]*fHx + xH2[0]*fHx; + v[1] += x[j][1]*fOy + xH1[1]*fHy + xH2[1]*fHy; + v[2] += x[j][2]*fOz + xH1[2]*fHz + xH2[2]*fHz; + v[3] += x[j][0]*fOy + xH1[0]*fHy + xH2[0]*fHy; + v[4] += x[j][0]*fOz + xH1[0]*fHz + xH2[0]*fHz; + v[5] += x[j][1]*fOz + xH1[1]*fHz + xH2[1]*fHz; + + vlist[n++] = j; + vlist[n++] = jH1; + vlist[n++] = jH2; + } + } + + if (EFLAG) { + if (CTABLE || rsq <= tabinnersq) + ecoul = prefactor*erfc; + 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 (EVFLAG) ev_tally_list_thr(this,n,vlist,ecoul,v,thr); + } + } + } + f[i][0] += fxtmp; + f[i][1] += fytmp; + f[i][2] += fztmp; + } +} + +/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + compute position xM of fictitious charge site for O atom and 2 H atoms + return it as xM +------------------------------------------------------------------------- */ + +void PairLJCutCoulLongProxyTIP4POMP::compute_newsite_thr(const double * xO, + const double * xH1, + const double * xH2, + double * xM) const +{ + double delx1 = xH1[0] - xO[0]; + double dely1 = xH1[1] - xO[1]; + double delz1 = xH1[2] - xO[2]; + domain->minimum_image(delx1,dely1,delz1); + + double delx2 = xH2[0] - xO[0]; + double dely2 = xH2[1] - xO[1]; + double delz2 = xH2[2] - xO[2]; + domain->minimum_image(delx2,dely2,delz2); + + const double prefac = alpha * 0.5; + xM[0] = xO[0] + prefac * (delx1 + delx2); + xM[1] = xO[1] + prefac * (dely1 + dely2); + xM[2] = xO[2] + prefac * (delz1 + delz2); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJCutCoulLongProxyTIP4POMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += PairLJCutCoulLongTIP4P::memory_usage(); + return bytes; +} diff --git a/src/USER-OMP/pair_lj_cut_coul_long_proxy_tip4p_omp.h b/src/USER-OMP/pair_lj_cut_coul_long_proxy_tip4p_omp.h new file mode 100644 index 0000000000..c66bbe86e8 --- /dev/null +++ b/src/USER-OMP/pair_lj_cut_coul_long_proxy_tip4p_omp.h @@ -0,0 +1,57 @@ +/* -*- 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/long/proxy/tip4p/omp,PairLJCutCoulLongProxyTIP4POMP) + +#else + +#ifndef LMP_PAIR_LJ_CUT_COUL_LONG_PROXY_TIP4P_OMP_H +#define LMP_PAIR_LJ_CUT_COUL_LONG_PROXY_TIP4P_OMP_H + +#include "pair_lj_cut_coul_long_tip4p.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairLJCutCoulLongProxyTIP4POMP : + public PairLJCutCoulLongTIP4P, public ThrOMP { + + public: + PairLJCutCoulLongProxyTIP4POMP(class LAMMPS *); + virtual ~PairLJCutCoulLongProxyTIP4POMP() {}; + + virtual void init_style(); + + virtual void compute(int, int); + virtual double memory_usage(); + + private: + template + void eval(int ifrom, int ito, ThrData * const thr); + void compute_newsite_thr(const double *, const double *, + const double *, double *) const; + + class PPPMTIP4PProxy *kspace; + int nproxy; +}; + +} + +#endif +#endif