Commit 8778452e authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

convert a few more styles to use fast full DP erfc()

parent 6ae15d08

src/OPT/pair_lj_charmm_coul_long_opt.cpp

+13 −25
Original line number Diff line number Diff line
@@ -23,17 +23,13 @@ 
#include "pair_lj_charmm_coul_long_opt.h"

#include "atom.h"

#include "force.h"

#include "math_special.h"

#include "math_const.h"

#include "neigh_list.h"



using namespace LAMMPS_NS;



#define EWALD_F   1.12837917

#define EWALD_P   0.3275911

#define EWALD_A1  0.254829592

#define EWALD_A2 -0.284496736

#define EWALD_A3  1.421413741

#define EWALD_A4 -1.453152027

#define EWALD_A5  1.061405429

using namespace MathSpecial;

using namespace MathConst;



/* ---------------------------------------------------------------------- */
@@ -76,7 +72,7 @@ void PairLJCharmmCoulLongOpt::eval() 
  int i,j,ii,jj,inum,jnum,itype,jtype,itable,sbindex;

  double fraction,table;

  double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;

  double grij,expm2,prefactor,t,erfc;

  double grij,expm2,prefactor,erfc;

  double philj,switch1,switch2;



  double rsq;
@@ -156,13 +152,10 @@ void PairLJCharmmCoulLongOpt::eval() 
            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 *

                (EWALD_A1+t*(EWALD_A2+t*(EWALD_A3+t*(EWALD_A4+t*EWALD_A5)))) *

                expm2;

              expm2 = expmsq(grij);

              erfc = my_erfcx(grij) * expm2;

              prefactor = qqrd2e * tmp_coef3/r;

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

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

            } else {

              union_int_float_t rsq_lookup;

              rsq_lookup.f = rsq;
@@ -245,16 +238,11 @@ void PairLJCharmmCoulLongOpt::eval() 
            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 *

                (EWALD_A1+t*(EWALD_A2+t*(EWALD_A3+t*(EWALD_A4+t*EWALD_A5)))) *

                expm2;

              prefactor = qqrd2e * tmp_coef3/r;

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

              if (factor_coul < 1.0) {

                forcecoul -= (1.0-factor_coul)*prefactor;

              }

              expm2 = expmsq(grij);

              erfc = my_erfcx(grij) * expm2;

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

              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;

src/OPT/pair_lj_cut_coul_long_opt.cpp

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

#include "atom.h"

#include "force.h"

#include "math_special.h"

#include "math_const.h"

#include "neigh_list.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;



/* ---------------------------------------------------------------------- */
@@ -79,7 +75,7 @@ void PairLJCutCoulLongOpt::eval() 
  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;

  double fraction,table;

  double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;

  double grij,expm2,prefactor,t,erfc;

  double grij,expm2,prefactor,erfc;

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

  double rsq;
@@ -132,11 +128,10 @@ void PairLJCutCoulLongOpt::eval() 
          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;

            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;

src/OPT/pair_lj_cut_tip4p_long_opt.cpp

+8 −13
Original line number Diff line number Diff line
@@ -22,18 +22,14 @@ 
#include "force.h"

#include "error.h"

#include "memory.h"

#include "math_special.h"

#include "math_const.h"

#include "neighbor.h"

#include "neigh_list.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;



/* ---------------------------------------------------------------------- */
@@ -110,7 +106,7 @@ void PairLJCutTIP4PLongOpt::eval() 
  double fraction,table;

  double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce;

  double factor_coul,factor_lj;

  double grij,expm2,prefactor,t,erfc;

  double grij,expm2,prefactor,erfc;

  double v[6];

  double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;

  const double *x1,*x2,*xH1,*xH2;
@@ -272,11 +268,10 @@ void PairLJCutTIP4PLongOpt::eval() 
          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;

            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;

            }

src/pair_coul_dsf.cpp

+13 −18
Original line number Diff line number Diff line
@@ -27,20 +27,14 @@ 
#include "neighbor.h"

#include "neigh_list.h"

#include "memory.h"

#include "math_special.h"

#include "math_const.h"

#include "error.h"



using namespace LAMMPS_NS;

using namespace MathSpecial;

using namespace MathConst;



#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



/* ---------------------------------------------------------------------- */



PairCoulDSF::PairCoulDSF(LAMMPS *lmp) : Pair(lmp) {}
@@ -64,7 +58,7 @@ void PairCoulDSF::compute(int eflag, int vflag) 
  int i,j,ii,jj,inum,jnum;

  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;

  double r,rsq,r2inv,forcecoul,factor_coul;

  double prefactor,erfcc,erfcd,t;

  double prefactor,erfcc,erfcd,grij;

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



  ecoul = 0.0;
@@ -115,9 +109,9 @@ void PairCoulDSF::compute(int eflag, int vflag) 


        r = sqrt(rsq);

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

        erfcd = exp(-alpha*alpha*rsq);

        t = 1.0 / (1.0 + EWALD_P*alpha*r);

        erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;

        grij = alpha*r;

        erfcd = expmsq(grij);

        erfcc = my_erfcx(grij) * erfcd;

        forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +

                                 r*f_shift) * r;

        if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
@@ -295,7 +289,7 @@ double PairCoulDSF::single(int i, int j, int itype, int jtype, double rsq, 
                           double factor_coul, double factor_lj,

                           double &fforce)

{

  double r2inv,r,erfcc,erfcd,prefactor,t;

  double r2inv,r,erfcc,erfcd,prefactor,grij;

  double forcecoul,phicoul;



  r2inv = 1.0/rsq;
@@ -303,15 +297,16 @@ double PairCoulDSF::single(int i, int j, int itype, int jtype, double rsq, 
  double eng = 0.0;

  if (rsq < cut_coulsq) {

    r = sqrt(rsq);

    prefactor = factor_coul * force->qqrd2e * atom->q[i]*atom->q[j]/r;

    erfcd = exp(-alpha*alpha*rsq);

    t = 1.0 / (1.0 + EWALD_P*alpha*r);

    erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;



    prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;

    grij = alpha * r;

    erfcd = expmsq(grij);

    erfcc = my_erfcx(grij) * erfcd;

    forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS*erfcd +

      r*f_shift) * r;

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



    phicoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);

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

    eng += phicoul;

  } else forcecoul = 0.0;

src/pair_lj_cut_coul_dsf.cpp

+6 −12
Original line number Diff line number Diff line
@@ -27,20 +27,14 @@ 
#include "neighbor.h"

#include "neigh_list.h"

#include "memory.h"

#include "math_special.h"

#include "math_const.h"

#include "error.h"



using namespace LAMMPS_NS;

using namespace MathSpecial;

using namespace MathConst;



#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



/* ---------------------------------------------------------------------- */



PairLJCutCoulDSF::PairLJCutCoulDSF(LAMMPS *lmp) : Pair(lmp)
@@ -77,7 +71,7 @@ void PairLJCutCoulDSF::compute(int eflag, int vflag) 
  int i,j,ii,jj,inum,jnum,itype,jtype;

  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;

  double r,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;

  double prefactor,erfcc,erfcd,t;

  double prefactor,erfcc,erfcd,grij;

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



  evdwl = ecoul = 0.0;
@@ -139,9 +133,9 @@ void PairLJCutCoulDSF::compute(int eflag, int vflag) 
        if (rsq < cut_coulsq) {

          r = sqrt(rsq);

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

          erfcd = exp(-alpha*alpha*r*r);

          t = 1.0 / (1.0 + EWALD_P*alpha*r);

          erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;

          grij = alpha*r;

          erfcd = expmsq(grij);

          erfcc = my_erfcx(grij) * erfcd;

          forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +

                                   r*f_shift) * r;

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

CRA Git | Maintained and supported by SUSTech CRA and CCSE