Commit 381d87b7 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

convert pair styles in USER-OMP to use fast DP analytical coulomb

parent 8778452e

src/USER-OMP/pair_born_coul_long_omp.cpp

+35 −19
Original line number Diff line number Diff line
@@ -17,19 +17,15 @@ 
#include "atom.h"

#include "comm.h"

#include "force.h"

#include "math_special.h"

#include "math_const.h"

#include "neighbor.h"

#include "neigh_list.h"



#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

using namespace MathSpecial;

using namespace MathConst;



/* ---------------------------------------------------------------------- */
@@ -89,8 +85,8 @@ void PairBornCoulLongOMP::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 grij,expm2,prefactor,t,erfc;

  int *ilist,*jlist,*numneigh,**firstneigh;

  double grij,expm2,prefactor,erfc,table,fraction;

  int *ilist,*jlist,*numneigh,**firstneigh,itable;



  evdwl = ecoul = 0.0;
@@ -136,19 +132,34 @@ void PairBornCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr) 


      if (rsq < cutsq[itype][jtype]) {

        r2inv = 1.0/rsq;

        r = sqrt(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;

            expm2 = expmsq(grij);

            erfc = my_erfcx(grij) * expm2;

            prefactor = qqrd2e * qtmp*q[j]/r;

          forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);

            forcecoul = prefactor * (erfc + MY_ISPI4*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]) {

          r = sqrt(rsq);

          r6inv = r2inv*r2inv*r2inv;

          rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]);

          forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv
@@ -168,7 +179,12 @@ void PairBornCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr) 


        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]) {

src/USER-OMP/pair_buck_coul_long_omp.cpp

+35 −19
Original line number Diff line number Diff line
@@ -17,19 +17,15 @@ 
#include "atom.h"

#include "comm.h"

#include "force.h"

#include "math_special.h"

#include "math_const.h"

#include "neighbor.h"

#include "neigh_list.h"



#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

using namespace MathSpecial;

using namespace MathConst;



/* ---------------------------------------------------------------------- */
@@ -90,8 +86,8 @@ void PairBuckCoulLongOMP::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 grij,expm2,prefactor,t,erfc;

  int *ilist,*jlist,*numneigh,**firstneigh;

  double grij,expm2,prefactor,erfc,fraction,table;

  int *ilist,*jlist,*numneigh,**firstneigh,itable;



  evdwl = ecoul = 0.0;
@@ -137,19 +133,34 @@ void PairBuckCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr) 


      if (rsq < cutsq[itype][jtype]) {

        r2inv = 1.0/rsq;

        r = sqrt(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;

            expm2 = expmsq(grij);

            erfc = my_erfcx(grij) * expm2;

            prefactor = qqrd2e * qtmp*q[j]/r;

          forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);

            forcecoul = prefactor * (erfc + MY_ISPI4*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]) {

          r = sqrt(rsq);

          r6inv = r2inv*r2inv*r2inv;

          rexp = exp(-r*rhoinv[itype][jtype]);

          forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv;
@@ -168,7 +179,12 @@ void PairBuckCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr) 


        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]) {

src/USER-OMP/pair_coul_long_omp.cpp

+8 −13
Original line number Diff line number Diff line
@@ -17,19 +17,15 @@ 
#include "atom.h"

#include "comm.h"

#include "force.h"

#include "math_special.h"

#include "math_const.h"

#include "neighbor.h"

#include "neigh_list.h"



#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

using namespace MathSpecial;

using namespace MathConst;



/* ---------------------------------------------------------------------- */
@@ -91,7 +87,7 @@ void PairCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;

  double fraction,table;

  double r,r2inv,rsq,forcecoul,factor_coul;

  double grij,expm2,prefactor,t,erfc;

  double grij,expm2,prefactor,erfc;

  int *ilist,*jlist,*numneigh,**firstneigh;



  ecoul = 0.0;
@@ -139,11 +135,10 @@ void PairCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
        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;

          expm2 = expmsq(grij);

          erfc = my_erfcx(grij) * expm2;

          prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;

          forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);

          forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);

          if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;

        } else {

          union_int_float_t rsq_lookup;

src/USER-OMP/pair_lj_charmm_coul_long_omp.cpp

+7 −12
Original line number Diff line number Diff line
@@ -17,11 +17,15 @@ 
#include "atom.h"

#include "comm.h"

#include "force.h"

#include "math_special.h"

#include "math_const.h"

#include "neighbor.h"

#include "neigh_list.h"



#include "suffix.h"

using namespace LAMMPS_NS;

using namespace MathSpecial;

using namespace MathConst;



/* ---------------------------------------------------------------------- */
@@ -132,21 +136,12 @@ void PairLJCharmmCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr) 


        if (rsq < cut_coulsq) {

          if (!ncoultablebits || rsq <= tabinnersq) {

            const double A1 =  0.254829592;

            const double A2 = -0.284496736;

            const double A3 =  1.421413741;

            const double A4 = -1.453152027;

            const double A5 =  1.061405429;

            const double EWALD_F = 1.12837917;

            const double INV_EWALD_P = 1.0/0.3275911;



            const double r = sqrt(rsq);

            const double grij = g_ewald * r;

            const double expm2 = exp(-grij*grij);

            const double t = INV_EWALD_P / (INV_EWALD_P + grij);

            const double erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;

            const double expm2 = expmsq(grij);

            const double erfc = my_erfcx(grij) * expm2;

            const double prefactor = qqrd2e * qtmp*q[j]/r;

            forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);

            forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);

            if (EFLAG) ecoul = prefactor*erfc;

            if (sbindex) {

              const double adjust = (1.0-special_coul[sbindex])*prefactor;

src/USER-OMP/pair_lj_class2_coul_long_omp.cpp

+8 −13
Original line number Diff line number Diff line
@@ -17,19 +17,15 @@ 
#include "atom.h"

#include "comm.h"

#include "force.h"

#include "math_special.h"

#include "math_const.h"

#include "neighbor.h"

#include "neigh_list.h"



#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

using namespace MathSpecial;

using namespace MathConst;



/* ---------------------------------------------------------------------- */
@@ -89,7 +85,7 @@ void PairLJClass2CoulLongOMP::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;

  double grij,expm2,prefactor,erfc;

  int *ilist,*jlist,*numneigh,**firstneigh;



  evdwl = ecoul = 0.0;
@@ -140,11 +136,10 @@ void PairLJClass2CoulLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
        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;

          expm2 = expmsq(grij);

          erfc = my_erfcx(grij) * expm2;

          prefactor = qqrd2e * qtmp*q[j]/r;

          forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);

          forcecoul = prefactor * (erfc + MY_ISPI4*grij*expm2);

          if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;

        } else forcecoul = 0.0;

CRA Git | Maintained and supported by SUSTech CRA and CCSE