port corrected triclinic handling from pppm/tip4p to pppm/tip4p/omp (f8364342) · Commits · 郑智淋 / lammps

src/USER-OMP/pppm_tip4p_omp.cpp

+94 −14

Original line number	Diff line number	Diff line
		@@ -48,7 +48,7 @@ using namespace MathSpecial;
		PPPMTIP4POMP::PPPMTIP4POMP(LAMMPS lmp, int narg, char *arg) :
		PPPMTIP4P(lmp, narg, arg), ThrOMP(lmp, THR_KSPACE)
		{
		triclinic_support = 0;
		triclinic_support = 1;
		suffix_flag \|= Suffix::OMP;
		}

		@@ -735,6 +735,8 @@ void PPPMTIP4POMP::fieldforce_ad()

		void PPPMTIP4POMP::find_M_thr(int i, int &iH1, int &iH2, dbl3_t &xM)
		{
		double **x = atom->x;

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

		@@ -742,6 +744,83 @@ void PPPMTIP4POMP::find_M_thr(int i, int &iH1, int &iH2, dbl3_t &xM)
		if (atom->type[iH1] != typeH \|\| atom->type[iH2] != typeH)
		error->one(FLERR,"TIP4P hydrogen has incorrect atom type");

		if (triclinic) {

		// need to use custom code to find the closest image for triclinic,
		// since local atoms are in lambda coordinates, but ghosts are not.

		int *sametag = atom->sametag;
		double xo[3],xh1[3],xh2[3];

		domain->lamda2x(x[i],xo);
		domain->lamda2x(x[iH1],xh1);
		domain->lamda2x(x[iH2],xh2);

		double delx = xo[0] - xh1[0];
		double dely = xo[1] - xh1[1];
		double delz = xo[2] - xh1[2];
		double rsqmin = delxdelx + delydely + delz*delz;
		double rsq;
		int closest = iH1;

		while (sametag[iH1] >= 0) {
		iH1 = sametag[iH1];
		delx = xo[0] - x[iH1][0];
		dely = xo[1] - x[iH1][1];
		delz = xo[2] - x[iH1][2];
		rsq = delxdelx + delydely + delz*delz;
		if (rsq < rsqmin) {
		rsqmin = rsq;
		closest = iH1;
		xh1[0] = x[iH1][0];
		xh1[1] = x[iH1][1];
		xh1[2] = x[iH1][2];
		}
		}
		iH1 = closest;

		closest = iH2;
		delx = xo[0] - xh2[0];
		dely = xo[1] - xh2[1];
		delz = xo[2] - xh2[2];
		rsqmin = delxdelx + delydely + delz*delz;

		while (sametag[iH2] >= 0) {
		iH2 = sametag[iH2];
		delx = xo[0] - x[iH2][0];
		dely = xo[1] - x[iH2][1];
		delz = xo[2] - x[iH2][2];
		rsq = delxdelx + delydely + delz*delz;
		if (rsq < rsqmin) {
		rsqmin = rsq;
		closest = iH2;
		xh2[0] = x[iH2][0];
		xh2[1] = x[iH2][1];
		xh2[2] = x[iH2][2];
		}
		}
		iH2 = closest;

		// finally compute M in real coordinates ...

		double delx1 = xh1[0] - xo[0];
		double dely1 = xh1[1] - xo[1];
		double delz1 = xh1[2] - xo[2];

		double delx2 = xh2[0] - xo[0];
		double dely2 = xh2[1] - xo[1];
		double delz2 = xh2[2] - xo[2];

		xM.x = xo[0] + alpha * 0.5 * (delx1 + delx2);
		xM.y = xo[1] + alpha * 0.5 * (dely1 + dely2);
		xM.z = xo[2] + alpha * 0.5 * (delz1 + delz2);

		// ... and convert M to lamda space for PPPM

		domain->x2lamda((double )&xM,(double )&xM);

		} else {

		// set iH1,iH2 to index of closest image to O

		iH1 = domain->closest_image(i,iH1);
		@@ -761,6 +840,7 @@ void PPPMTIP4POMP::find_M_thr(int i, int &iH1, int &iH2, dbl3_t &xM)
		xM.y = x[i].y + alpha * 0.5 * (dely1 + dely2);
		xM.z = x[i].z + alpha * 0.5 * (delz1 + delz2);
		}
		}


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

Admin message