Unverified Commit 2775b937 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer Committed by GitHub
Browse files

Merge pull request #2282 from akohlmey/more-unit-tests 

Add more unit tests for pair and kspace styles
parents aa393f35 54b93316

src/USER-CGSDK/pair_lj_sdk_coul_long.cpp

+25 −10
Original line number Diff line number Diff line
@@ -593,17 +593,32 @@ double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype, 


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

    const int ljt = lj_type[itype][jtype];

    const double ljpow1 = lj_pow1[ljt];

    const double ljpow2 = lj_pow2[ljt];

    const double ljpref = lj_prefact[ljt];



    const double ratio = sigma[itype][jtype]/sqrt(rsq);

    const double eps = epsilon[itype][jtype];

    if (ljt == LJ12_4) {

      const double r4inv=r2inv*r2inv;

      forcelj = r4inv*(lj1[itype][jtype]*r4inv*r4inv

                       - lj2[itype][jtype]);



      philj = r4inv*(lj3[itype][jtype]*r4inv*r4inv

                     - lj4[itype][jtype]) - offset[itype][jtype];



    } else if (ljt == LJ9_6) {

      const double r3inv = r2inv*sqrt(r2inv);

      const double r6inv = r3inv*r3inv;

      forcelj = r6inv*(lj1[itype][jtype]*r3inv

                       - lj2[itype][jtype]);

      philj = r6inv*(lj3[itype][jtype]*r3inv

                     - lj4[itype][jtype]) - offset[itype][jtype];



    forcelj = factor_lj * ljpref*eps * (ljpow1*pow(ratio,ljpow1)

                          - ljpow2*pow(ratio,ljpow2))/rsq;

    philj = factor_lj * (ljpref*eps * (pow(ratio,ljpow1) - pow(ratio,ljpow2))

                         - offset[itype][jtype]);

    } else if (ljt == LJ12_6) {

      const double r6inv = r2inv*r2inv*r2inv;

      forcelj = r6inv*(lj1[itype][jtype]*r6inv

                       - lj2[itype][jtype]);

      philj = r6inv*(lj3[itype][jtype]*r6inv

                     - lj4[itype][jtype]) - offset[itype][jtype];

    }

    forcelj *= factor_lj;

    philj *= factor_lj;

  }



  fforce = (forcecoul + forcelj) * r2inv;

src/USER-CGSDK/pair_lj_sdk_coul_msm.cpp

+25 −10
Original line number Diff line number Diff line
@@ -257,17 +257,32 @@ double PairLJSDKCoulMSM::single(int i, int j, int itype, int jtype, 


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

    const int ljt = lj_type[itype][jtype];

    const double ljpow1 = lj_pow1[ljt];

    const double ljpow2 = lj_pow2[ljt];

    const double ljpref = lj_prefact[ljt];



    const double ratio = sigma[itype][jtype]/sqrt(rsq);

    const double eps = epsilon[itype][jtype];

    if (ljt == LJ12_4) {

      const double r4inv=r2inv*r2inv;

      forcelj = r4inv*(lj1[itype][jtype]*r4inv*r4inv

                       - lj2[itype][jtype]);



      philj = r4inv*(lj3[itype][jtype]*r4inv*r4inv

                     - lj4[itype][jtype]) - offset[itype][jtype];



    } else if (ljt == LJ9_6) {

      const double r3inv = r2inv*sqrt(r2inv);

      const double r6inv = r3inv*r3inv;

      forcelj = r6inv*(lj1[itype][jtype]*r3inv

                       - lj2[itype][jtype]);

      philj = r6inv*(lj3[itype][jtype]*r3inv

                     - lj4[itype][jtype]) - offset[itype][jtype];



    forcelj = factor_lj * ljpref*eps * (ljpow1*pow(ratio,ljpow1)

                          - ljpow2*pow(ratio,ljpow2))/rsq;

    philj = factor_lj * (ljpref*eps * (pow(ratio,ljpow1) - pow(ratio,ljpow2))

                         - offset[itype][jtype]);

    } else if (ljt == LJ12_6) {

      const double r6inv = r2inv*r2inv*r2inv;

      forcelj = r6inv*(lj1[itype][jtype]*r6inv

                       - lj2[itype][jtype]);

      philj = r6inv*(lj3[itype][jtype]*r6inv

                     - lj4[itype][jtype]) - offset[itype][jtype];

    }

    forcelj *= factor_lj;

    philj *= factor_lj;

  }



  fforce = (forcecoul + forcelj) * r2inv;

src/USER-MISC/pair_coul_diel.cpp

+2 −0
Original line number Diff line number Diff line
@@ -281,6 +281,8 @@ void PairCoulDiel::read_restart(FILE *fp) 
  int me = comm->me;

  for (i = 1; i <= atom->ntypes; i++)

    for (j = i; j <= atom->ntypes; j++) {

      if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,NULL,error);

      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);

      if (setflag[i][j]) {

        if (me == 0) {

          utils::sfread(FLERR,&rme[i][j],sizeof(double),1,fp,NULL,error);

src/USER-MISC/pair_coul_shield.cpp

+33 −24
Original line number Diff line number Diff line
@@ -115,8 +115,8 @@ void PairCoulShield::compute(int eflag, int vflag) 


        // turn on/off taper function

        if (tap_flag) {

          Tap = calc_Tap(r,sqrt(cutsq[itype][jtype]));

          dTap = calc_dTap(r,sqrt(cutsq[itype][jtype]));

          Tap = calc_Tap(r,cut[itype][jtype]);

          dTap = calc_dTap(r,cut[itype][jtype]);

        } else {Tap = 1.0; dTap = 0.0;}



        forcecoul = qqrd2e*qtmp*q[j]*r*depsdr;
@@ -297,6 +297,8 @@ void PairCoulShield::read_restart(FILE *fp) 
  int me = comm->me;

  for (i = 1; i <= atom->ntypes; i++)

    for (j = i; j <= atom->ntypes; j++) {

      if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,NULL,error);

      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);

      if (setflag[i][j]) {

        if (me == 0) {

          utils::sfread(FLERR,&sigmae[i][j],sizeof(double),1,fp,NULL,error);
@@ -346,6 +348,13 @@ double PairCoulShield::single(int i, int j, int itype, int jtype, 
  double *q = atom->q;

  double qqrd2e = force->qqrd2e;



  // only computed between different layers as indicated by different molecule ids.



  if (atom->molecule[i] == atom->molecule[j]) {

    fforce = 0.0;

    return 0.0;

  }



  r = sqrt(rsq);

  r3 = rsq*r;

  rarg = 1.0/sigmae[itype][jtype];
@@ -357,15 +366,15 @@ double PairCoulShield::single(int i, int j, int itype, int jtype, 


  // turn on/off taper function

  if (tap_flag) {

     Tap = calc_Tap(r,sqrt(cutsq[itype][jtype]));

     dTap = calc_dTap(r,sqrt(cutsq[itype][jtype]));

    Tap = calc_Tap(r,cut[itype][jtype]);

    dTap = calc_dTap(r,cut[itype][jtype]);

  } else {Tap = 1.0; dTap = 0.0;}



  forcecoul = qqrd2e*q[i]*q[j]*r*depsdr;

  fvc = forcecoul*Tap - Vc*dTap/r;

  fforce = factor_coul*fvc;



  if (tap_flag) phishieldec = factor_coul*Vc*Tap;

  if (tap_flag) phishieldec = Vc*Tap;

  else phishieldec = Vc - offset[itype][jtype];

  return factor_coul*phishieldec;

}

src/USER-MISC/pair_meam_spline.cpp

+16 −3
Original line number Diff line number Diff line
@@ -348,6 +348,7 @@ void PairMEAMSpline::allocate() 
  zero_atom_energies = new double[n];



  map = new int[n+1];

  for (int i=0; i <= n; ++i) map[i] = -1;

}



/* ----------------------------------------------------------------------
@@ -400,7 +401,7 @@ void PairMEAMSpline::coeff(int narg, char **arg) 
        if (strcmp(arg[i],elements[j]) == 0)

          break;

      if (j < nelements) map[i-2] = j;

      else error->all(FLERR,"No matching element in EAM potential file");

      else error->all(FLERR,"No matching element in meam/spline potential file");

    }

  }

  // clear setflag since coeff() called once with I,J = * *
@@ -419,8 +420,17 @@ void PairMEAMSpline::coeff(int narg, char **arg) 
        setflag[i][j] = 1;

        count++;

      }



  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");



  // check that each element is mapped to exactly one atom type



  for (int i = 0; i < nelements; i++) {

    count = 0;

    for (int j = 1; j <= n; j++)

      if (map[j] == i) count++;

    if (count != 1)

      error->all(FLERR,"Pair style meam/spline requires one atom type per element");

  } 

}



#define MAXLINE 1024
@@ -480,7 +490,10 @@ void PairMEAMSpline::read_file(const char* filename) 
    }



    nmultichoose2 = ((nelements+1)*nelements)/2;

    // allocate!!



    if (nelements != atom->ntypes)

      error->all(FLERR,"Pair style meam/spline requires one atom type per element");



    allocate();



    // Parse spline functions.

CRA Git | Maintained and supported by SUSTech CRA and CCSE