Commit 5c5b4ffa authored by sjplimp's avatar sjplimp Committed by GitHub
Browse files

Merge pull request #522 from akohlmey/tip4p-cleanup-refactor 

Refactor and bugfix for some TIP4P pair styles
parents 30177c4e 2e728972

src/KSPACE/pair_lj_long_tip4p_long.cpp

+390 −363
Original line number Diff line number Diff line
@@ -158,7 +158,6 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag) 
        hneigh[i][0] = iH1;

        hneigh[i][1] = iH2;

        hneigh[i][2] = 1;



      } else {

        iH1 = hneigh[i][0];

        iH2 = hneigh[i][1];
@@ -272,7 +271,6 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag) 
              hneigh[j][0] = jH1;

              hneigh[j][1] = jH2;

              hneigh[j][2] = 1;



            } else {

              jH1 = hneigh[j][0];

              jH2 = hneigh[j][1];
@@ -473,7 +471,7 @@ void PairLJLongTIP4PLong::compute_inner() 
  int iH1,iH2,jH1,jH2;

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

  double r2inv,forcecoul,forcelj,cforce;

  double fO[3],fH[3],fd[3];// f1[3];

  double fO[3],fH[3],fd[3];

  double *x1,*x2;

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

  double rsq, qri;
@@ -534,14 +532,19 @@ void PairLJLongTIP4PLong::compute_inner() 
    itype = type[i];

    if (itype == typeO && order1) {

      if (hneigh[i][0] < 0) {

        hneigh[i][0] = iH1 = atom->map(tag[i] + 1);

        hneigh[i][1] = iH2 = atom->map(tag[i] + 2);

        hneigh[i][2] = 1;

        iH1 = atom->map(tag[i] + 1);

        iH2 = atom->map(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");

        // set iH1,iH2 to closest image to O

        iH1 = domain->closest_image(i,iH1);

        iH2 = domain->closest_image(i,iH2);

        compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);

        hneigh[i][0] = iH1;

        hneigh[i][1] = iH2;

        hneigh[i][2] = 1;

      } else {

        iH1 = hneigh[i][0];

        iH2 = hneigh[i][1];
@@ -599,14 +602,19 @@ void PairLJLongTIP4PLong::compute_inner() 
        if (itype == typeO || jtype == typeO) {

          if (jtype == typeO) {

            if (hneigh[j][0] < 0) {

              hneigh[j][0] = jH1 = atom->map(tag[j] + 1);

              hneigh[j][1] = jH2 = atom->map(tag[j] + 2);

              hneigh[j][2] = 1;

              jH1 = atom->map(tag[j] + 1);

              jH2 = atom->map(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");

              // set jH1,jH2 to closest image to O

              jH1 = domain->closest_image(j,jH1);

              jH2 = domain->closest_image(j,jH2);

              compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);

              hneigh[j][0] = jH1;

              hneigh[j][1] = jH2;

              hneigh[j][2] = 1;

            } else {

              jH1 = hneigh[j][0];

              jH2 = hneigh[j][1];
@@ -777,14 +785,19 @@ void PairLJLongTIP4PLong::compute_middle() 
    itype = type[i];

    if (itype == typeO && order1) {

      if (hneigh[i][0] < 0) {

        hneigh[i][0] = iH1 = atom->map(tag[i] + 1);

        hneigh[i][1] = iH2 = atom->map(tag[i] + 2);

        hneigh[i][2] = 1;

        iH1 = atom->map(tag[i] + 1);

        iH2 = atom->map(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");

        // set iH1,iH2 to closest image to O

        iH1 = domain->closest_image(i,iH1);

        iH2 = domain->closest_image(i,iH2);

        compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);

        hneigh[i][0] = iH1;

        hneigh[i][1] = iH2;

        hneigh[i][2] = 1;

      } else {

        iH1 = hneigh[i][0];

        iH2 = hneigh[i][1];
@@ -846,14 +859,19 @@ void PairLJLongTIP4PLong::compute_middle() 
        if (itype == typeO || jtype == typeO) {

          if (jtype == typeO) {

            if (hneigh[j][0] < 0) {

              hneigh[j][0] = jH1 = atom->map(tag[j] + 1);

              hneigh[j][1] = jH2 = atom->map(tag[j] + 2);

              hneigh[j][2] = 1;

              jH1 = atom->map(tag[j] + 1);

              jH2 = atom->map(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");

              // set jH1,jH2 to closest image to O

              jH1 = domain->closest_image(j,jH1);

              jH2 = domain->closest_image(j,jH2);

              compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);

              hneigh[j][0] = jH1;

              hneigh[j][1] = jH2;

              hneigh[j][2] = 1;

            } else {

              jH1 = hneigh[j][0];

              jH2 = hneigh[j][1];
@@ -979,8 +997,8 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag) 
  int iH1,iH2,jH1,jH2;

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

  double r2inv,forcecoul,forcelj,cforce, respa_coul, respa_lj, frespa,fvirial;

  double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];// f1[3];

  double *x1,*x2;

  double fO[3],fH[3],fd[3],v[6];

  double *x1,*x2,*xH1,*xH2;

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

  double rsq,qri;

  int respa_flag;
@@ -1048,14 +1066,19 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag) 
    itype = type[i];

    if (itype == typeO) {

      if (hneigh[i][0] < 0) {

        hneigh[i][0] = iH1 = atom->map(tag[i] + 1);

        hneigh[i][1] = iH2 = atom->map(tag[i] + 2);

        hneigh[i][2] = 1;

        iH1 = atom->map(tag[i] + 1);

        iH2 = atom->map(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");

        // set iH1,iH2 to closest image to O

        iH1 = domain->closest_image(i,iH1);

        iH2 = domain->closest_image(i,iH2);

        compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);

        hneigh[i][0] = iH1;

        hneigh[i][1] = iH2;

        hneigh[i][2] = 1;

      } else {

        iH1 = hneigh[i][0];

        iH2 = hneigh[i][1];
@@ -1167,14 +1190,19 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag) 
        if (itype == typeO || jtype == typeO) {

          if (jtype == typeO) {

            if (hneigh[j][0] < 0) {

              hneigh[j][0] = jH1 = atom->map(tag[j] + 1);

              hneigh[j][1] = jH2 = atom->map(tag[j] + 2);

              hneigh[j][2] = 1;

              jH1 = atom->map(tag[j] + 1);

              jH2 = atom->map(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");

              // set jH1,jH2 to closest image to O

              jH1 = domain->closest_image(j,jH1);

              jH2 = domain->closest_image(j,jH2);

              compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);

              hneigh[j][0] = jH1;

              hneigh[j][1] = jH2;

              hneigh[j][2] = 1;

            } else {

              jH1 = hneigh[j][0];

              jH2 = hneigh[j][1];
@@ -1309,9 +1337,8 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag) 
              fH[1] = 0.5 * alpha * fd[1];

              fH[2] = 0.5 * alpha * fd[2];



	      domain->closest_image(x[i],x[iH1],xH1);

	      domain->closest_image(x[i],x[iH2],xH2);



              xH1 = x[jH1];

              xH2 = x[jH2];

              v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];

              v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];

              v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -1380,8 +1407,8 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag) 
              fH[1] = 0.5 * alpha * fd[1];

              fH[2] = 0.5 * alpha * fd[2];



	      domain->closest_image(x[j],x[jH1],xH1);

	      domain->closest_image(x[j],x[jH2],xH2);

              xH1 = x[jH1];

              xH2 = x[jH2];



              v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];

              v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];

src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp

+13 −16
Original line number Diff line number Diff line
@@ -104,7 +104,7 @@ void PairLJCutTIP4PLongOMP::compute(int eflag, int vflag) 
    thr->timer(Timer::START);

    ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);



  if (!ncoultablebits) {

  if (ncoultablebits) {

    if (evflag) {

      if (eflag) {

        if (vflag) eval<1,1,1,1>(ifrom, ito, thr);
@@ -156,6 +156,7 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
  dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];

  const double * _noalias const q = atom->q;

  const int * _noalias const type = atom->type;

  const tagint * _noalias const tag = atom->tag;

  const int nlocal = atom->nlocal;

  const double * _noalias const special_coul = force->special_coul;

  const double * _noalias const special_lj = force->special_lj;
@@ -187,8 +188,8 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
    // will be the same, there is no race condition.

    if (itype == typeO) {

      if (hneigh_thr[i].a < 0) {

        iH1 = atom->map(atom->tag[i] + 1);

        iH2 = atom->map(atom->tag[i] + 2);

        iH1 = atom->map(tag[i] + 1);

        iH2 = atom->map(tag[i] + 2);

        if (iH1 == -1 || iH2 == -1)

          error->one(FLERR,"TIP4P hydrogen is missing");

        if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
@@ -267,8 +268,8 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr) 


          if (jtype == typeO) {

            if (hneigh_thr[j].a < 0) {

              jH1 = atom->map(atom->tag[j] + 1);

              jH2 = atom->map(atom->tag[j] + 2);

              jH1 = atom->map(tag[j] + 1);

              jH2 = atom->map(tag[j] + 2);

              if (jH1 == -1 || jH2 == -1)

                error->one(FLERR,"TIP4P hydrogen is missing");

              if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
@@ -301,7 +302,7 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr) 


        if (rsq < cut_coulsq) {

          r2inv = 1 / rsq;

          if (CTABLE || rsq <= tabinnersq) {

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

            r = sqrt(rsq);

            grij = g_ewald * r;

            expm2 = exp(-grij*grij);
@@ -337,7 +338,7 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
          // 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



          if (EVFLAG) {

          if (VFLAG) {

            n = 0;

            key = 0;

          }
@@ -354,11 +355,11 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
              v[3] = x[i].x * dely * cforce;

              v[4] = x[i].x * delz * cforce;

              v[5] = x[i].y * delz * cforce;

              vlist[n++] = i;

            }

            if (EVFLAG) vlist[n++] = i;



          } else {

            if (EVFLAG) key++;

            if (VFLAG) key++;



            fdx = delx*cforce;

            fdy = dely*cforce;
@@ -393,8 +394,6 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
              v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;

              v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;

              v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;

            }

            if (EVFLAG) {

              vlist[n++] = i;

              vlist[n++] = iH1;

              vlist[n++] = iH2;
@@ -413,11 +412,11 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
              v[3] -= x[j].x * dely * cforce;

              v[4] -= x[j].x * delz * cforce;

              v[5] -= x[j].y * delz * cforce;

              vlist[n++] = j;

            }

            if (EVFLAG) vlist[n++] = j;



          } else {

            if (EVFLAG) key += 2;

            if (VFLAG) key += 2;



            fdx = -delx*cforce;

            fdy = -dely*cforce;
@@ -452,8 +451,6 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
              v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;

              v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;

              v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;

            }

            if (EVFLAG) {

              vlist[n++] = j;

              vlist[n++] = jH1;

              vlist[n++] = jH2;
@@ -461,7 +458,7 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
          }



          if (EFLAG) {

            if (CTABLE || rsq <= tabinnersq)

            if (!CTABLE || rsq <= tabinnersq)

              ecoul = prefactor*erfc;

            else {

              table = etable[itable] + fraction*detable[itable];

src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp

+630 −566

File changed.

Preview size limit exceeded, changes collapsed.

CRA Git | Maintained and supported by SUSTech CRA and CCSE