Commit 6ae15d08 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

first batch of pair styles converted to double precision erfc()

parent 2826449b

src/CLASS2/pair_lj_class2_coul_long.cpp

+10 −18
Original line number Diff line number Diff line
@@ -22,21 +22,15 @@ 
#include "kspace.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "math_special.h"

#include "math_const.h"

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



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



PairLJClass2CoulLong::PairLJClass2CoulLong(LAMMPS *lmp) : Pair(lmp)
@@ -76,7 +70,7 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag) 
  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;

  double fraction,table;

  double rsq,r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj;

  double grij,expm2,prefactor,t,erfc;

  double grij,expm2,prefactor,erfc;

  double factor_coul,factor_lj;

  int *ilist,*jlist,*numneigh,**firstneigh;
@@ -130,11 +124,10 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag) 
          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;
@@ -484,7 +477,7 @@ double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, 
                                    double factor_coul, double factor_lj,

                                    double &fforce)

{

  double r2inv,r,rinv,r3inv,r6inv,grij,expm2,t,erfc,prefactor;

  double r2inv,r,rinv,r3inv,r6inv,grij,expm2,erfc,prefactor;

  double fraction,table,forcecoul,forcelj,phicoul,philj;

  int itable;
@@ -493,11 +486,10 @@ double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, 
    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 = force->qqrd2e * atom->q[i]*atom->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/DIPOLE/pair_lj_cut_dipole_long.cpp

+162 −169
Original line number Diff line number Diff line
@@ -22,6 +22,7 @@ 
#include "neigh_list.h"

#include "force.h"

#include "kspace.h"

#include "math_special.h"

#include "math_const.h"

#include "memory.h"

#include "error.h"
@@ -30,16 +31,9 @@ 




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



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



PairLJCutDipoleLong::PairLJCutDipoleLong(LAMMPS *lmp) : Pair(lmp)
@@ -82,7 +76,7 @@ void PairLJCutDipoleLong::compute(int eflag, int vflag) 
  double tixcoul,tiycoul,tizcoul,tjxcoul,tjycoul,tjzcoul;

  double fx,fy,fz,fdx,fdy,fdz,fax,fay,faz;

  double pdotp,pidotr,pjdotr,pre1,pre2,pre3;

  double grij,expm2,t,erfc;

  double grij,expm2,erfc;

  double g0,g1,g2,b0,b1,b2,b3,d0,d1,d2,d3;

  double zdix,zdiy,zdiz,zdjx,zdjy,zdjz,zaix,zaiy,zaiz,zajx,zajy,zajz;

  double g0b1_g1b2_g2b3,g0d1_g1d2_g2d3;
@@ -112,14 +106,14 @@ void PairLJCutDipoleLong::compute(int eflag, int vflag) 
  firstneigh = list->firstneigh;



  pre1 = 2.0 * g_ewald / MY_PIS;

  pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS;

  pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS;

  pre2 = 4.0 * g_ewald*g_ewald*g_ewald / MY_PIS;

  pre3 = 8.0 * g_ewald*g_ewald*g_ewald*g_ewald*g_ewald / MY_PIS;



  // loop over neighbors of my atoms



  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    qtmp = atom->q[i];

    qtmp = q[i];

    xtmp = x[i][0];

    ytmp = x[i][1];

    ztmp = x[i][2];
@@ -146,9 +140,8 @@ void PairLJCutDipoleLong::compute(int eflag, int vflag) 
        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;



          pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2];

          pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;

src/KSPACE/pair_born_coul_long.cpp

+10 −18
Original line number Diff line number Diff line
@@ -26,21 +26,15 @@ 
#include "kspace.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "math_special.h"

#include "math_const.h"

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



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



PairBornCoulLong::PairBornCoulLong(LAMMPS *lmp) : Pair(lmp)
@@ -82,7 +76,7 @@ void PairBornCoulLong::compute(int eflag, int vflag) 
  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;

  double fraction,table;

  double rsq,r2inv,r6inv,forcecoul,forceborn,factor_coul,factor_lj;

  double grij,expm2,prefactor,t,erfc;

  double grij,expm2,prefactor,erfc;

  double r,rexp;

  int *ilist,*jlist,*numneigh,**firstneigh;
@@ -136,11 +130,10 @@ void PairBornCoulLong::compute(int eflag, int vflag) 
          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;
@@ -514,7 +507,7 @@ double PairBornCoulLong::single(int i, int j, int itype, int jtype, 
                                double factor_coul, double factor_lj,

                                double &fforce)

{

  double r2inv,r6inv,r,rexp,grij,expm2,t,erfc,prefactor;

  double r2inv,r6inv,r,rexp,grij,expm2,erfc,prefactor;

  double fraction,table,forcecoul,forceborn,phicoul,phiborn;

  int itable;
@@ -523,11 +516,10 @@ double PairBornCoulLong::single(int i, int j, int itype, int jtype, 
    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 = force->qqrd2e * atom->q[i]*atom->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/KSPACE/pair_buck_coul_long.cpp

+11 −18
Original line number Diff line number Diff line
@@ -22,21 +22,15 @@ 
#include "kspace.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "math_special.h"

#include "math_const.h"

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



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



PairBuckCoulLong::PairBuckCoulLong(LAMMPS *lmp) : Pair(lmp)
@@ -77,7 +71,7 @@ void PairBuckCoulLong::compute(int eflag, int vflag) 
  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;

  double fraction,table;

  double rsq,r2inv,r6inv,forcecoul,forcebuck,factor_coul,factor_lj;

  double grij,expm2,prefactor,t,erfc;

  double grij,expm2,prefactor,erfc;

  double r,rexp;

  int *ilist,*jlist,*numneigh,**firstneigh;
@@ -126,15 +120,15 @@ void PairBuckCoulLong::compute(int eflag, int vflag) 


      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;

            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;
@@ -486,7 +480,7 @@ double PairBuckCoulLong::single(int i, int j, int itype, int jtype, 
                                double factor_coul, double factor_lj,

                                double &fforce)

{

  double r2inv,r6inv,r,rexp,grij,expm2,t,erfc,prefactor;

  double r2inv,r6inv,r,rexp,grij,expm2,erfc,prefactor;

  double fraction,table,forcecoul,forcebuck,phicoul,phibuck;

  int itable;
@@ -495,11 +489,10 @@ double PairBuckCoulLong::single(int i, int j, int itype, int jtype, 
    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 = force->qqrd2e * atom->q[i]*atom->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/KSPACE/pair_coul_long.cpp

+14 −20
Original line number Diff line number Diff line
@@ -24,23 +24,19 @@ 
#include "comm.h"

#include "force.h"

#include "kspace.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "update.h"

#include "integrate.h"

#include "respa.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "math_special.h"

#include "math_const.h"

#include "memory.h"

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



/* ---------------------------------------------------------------------- */
@@ -74,7 +70,7 @@ void PairCoulLong::compute(int eflag, int vflag) 
  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;

  double fraction,table;

  double r,r2inv,forcecoul,factor_coul;

  double grij,expm2,prefactor,t,erfc;

  double grij,expm2,prefactor,erfc;

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

  double rsq;
@@ -124,11 +120,10 @@ void PairCoulLong::compute(int eflag, int vflag) 
        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;
@@ -343,7 +338,7 @@ double PairCoulLong::single(int i, int j, int itype, int jtype, 
                            double factor_coul, double factor_lj,

                            double &fforce)

{

  double r2inv,r,grij,expm2,t,erfc,prefactor;

  double r2inv,r,grij,expm2,erfc,prefactor;

  double fraction,table,forcecoul,phicoul;

  int itable;
@@ -351,11 +346,10 @@ double PairCoulLong::single(int i, int j, int itype, int jtype, 
  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 = force->qqrd2e * atom->q[i]*atom->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;

CRA Git | Maintained and supported by SUSTech CRA and CCSE