Chem Snap (d3d67ceb) · Commits · 郑智淋 / lammps

src/SNAP/compute_snap.cpp

+6 −4

Changes for src/SNAP/compute_snap.cpp: 6 added lines, 4 removed lines.
Original line number	Diff line number	Diff line
		@@ -40,7 +40,7 @@ ComputeSnap::ComputeSnap(LAMMPS lmp, int narg, char *arg) :
		extarray = 0;

		double rfac0, rmin0;
		int twojmax, switchflag, bzeroflag, bnormflag;
		int twojmax, switchflag, bzeroflag;
		radelem = NULL;
		wjelem = NULL;

		@@ -54,7 +54,6 @@ ComputeSnap::ComputeSnap(LAMMPS lmp, int narg, char *arg) :
		rmin0 = 0.0;
		switchflag = 1;
		bzeroflag = 1;
		bnormflag = 0;
		quadraticflag = 0;
		alloyflag = 0;
		wselfallflag = 0;
		@@ -113,10 +112,9 @@ ComputeSnap::ComputeSnap(LAMMPS lmp, int narg, char *arg) :
		quadraticflag = atoi(arg[iarg+1]);
		iarg += 2;
		} else if (strcmp(arg[iarg],"alloyflag") == 0) {
		if (iarg+2+ntypes > narg)
		if (iarg+2 > narg)
		error->all(FLERR,"Illegal compute snap command");
		alloyflag = 1;
		bnormflag = 1;
		memory->create(map,ntypes+1,"compute_snap:map");
		nelements = force->inumeric(FLERR,arg[iarg+1]);
		for(int i = 0; i < ntypes; i++) {
		@@ -168,6 +166,8 @@ ComputeSnap::~ComputeSnap()
		memory->destroy(wjelem);
		memory->destroy(cutsq);
		delete snaptr;

		if (alloyflag) memory->destroy(map);
		}

		/* ---------------------------------------------------------------------- */
		@@ -541,6 +541,8 @@ double ComputeSnap::memory_usage()
		sizeof(double); // snapall
		bytes += nmaxsize_peratom sizeof(double); // snap_peratom
		bytes += snaptr->memory_usage(); // SNA object
		int n = atom->ntypes+1;
		bytes += n*sizeof(int); // map

		return bytes;
		}

src/SNAP/sna.cpp

+7 −64

Changes for src/SNAP/sna.cpp: 7 added lines, 64 removed lines.
Original line number	Diff line number	Diff line
		@@ -277,8 +277,6 @@ void SNA::build_indexlist()
		idxz[idxz_count].mb1min = MAX(0, (2 * mb - j - j2 + j1) / 2);
		idxz[idxz_count].mb2max = (2 * mb - j - (2 * idxz[idxz_count].mb1min - j1) + j2) / 2;
		idxz[idxz_count].nb = MIN(j1, (2 * mb - j + j2 + j1) / 2) - idxz[idxz_count].mb1min + 1;
		idxz[idxz_count].ma = ma;
		idxz[idxz_count].mb = mb;
		// apply to z(j1,j2,j,ma,mb) to unique element of y(j)

		const int jju = idxu_block[j] + (j+1)*mb + ma;
		@@ -336,7 +334,6 @@ void SNA::compute_ui(int jnum, int ielem)
		// utot(j,ma,mb) += u(r0;j,ma,mb) for all j,ma,mb

		zero_uarraytot(ielem);
		addself_uarraytot(wself, ielem);

		for(int j = 0; j < jnum; j++) {
		x = rij[j][0];
		@@ -438,8 +435,6 @@ void SNA::compute_yi(const double* beta)
		int jju;
		double betaj;
		int itriple;
		int jelem;
		double temp;

		for(int ielem1 = 0; ielem1 < nelements; ielem1++)
		for(int j = 0; j <= twojmax; j++) {
		@@ -464,8 +459,6 @@ void SNA::compute_yi(const double* beta)
		const int mb1min = idxz[jjz].mb1min;
		const int mb2max = idxz[jjz].mb2max;
		const int nb = idxz[jjz].nb;
		const int ma = idxz[jjz].ma;
		const int mb = idxz[jjz].mb;

		const double *cgblock = cglist + idxcg_block[j1][j2][j];

		@@ -518,65 +511,29 @@ void SNA::compute_yi(const double* beta)
		jju = idxz[jjz].jju;
		for(int elem3 = 0; elem3 < nelements; elem3++) {
		// pick out right beta value
		if (alloy_flag) {
		if (j >= j1) {
		const int jjb = idxb_block[j1][j2][j];
		itriple = ((elem3 * nelements + elem2) * nelements + elem1) * idxb_max + jjb;
		if (j1 ==j && j2 == j && elem1 == elem2 && elem1 == elem3) betaj = 3 * beta[itriple];
		else if (j1 == j && j2 == j && elem2 == elem3 && elem1!=elem3)
		betaj =
		2 * (beta[itriple] + beta[((elem1 * nelements + elem2) * nelements + elem3) * idxb_max + jjb]);
		else if (j1 == j && j2 == j && elem1 == elem2 && elem1 != elem3)
		betaj = beta[itriple] + beta[((elem2 * nelements + elem1) * nelements + elem3) * idxb_max + jjb] +
		beta[((elem1 * nelements + elem3) * nelements + elem2) * idxb_max + jjb];
		else if (j1 ==j && (elem1 == elem3 \|\| elem2 == elem3)) // this line covers quite a few cases
		betaj = 2 * beta[itriple];
		else if (j1 == j && j2 != j && elem2 == elem1 && elem1 != elem3)
		betaj = beta[itriple] + beta[((elem1 * nelements + elem2) * nelements + elem3) * idxb_max + jjb];
		else
		betaj = beta[((elem1 * nelements + elem2) * nelements + elem3) * idxb_max + jjb];
		itriple = ((elem1 * nelements + elem2) * nelements + elem3) * idxb_max + jjb;
		if (j1 == j) {
		if (j2 == j) betaj = 3*beta[itriple];
		else betaj = 2*beta[itriple];
		} else betaj = beta[itriple];
		} else if (j >= j2) {
		const int jjb = idxb_block[j][j2][j1];
		itriple = ((elem3 * nelements + elem2) * nelements + elem1) * idxb_max + jjb;
		if (j2 == j) {
		if (elem3 == elem2)
		betaj = 2 * beta[itriple];
		else if (elem1 == elem2)
		betaj = beta[itriple] + beta[((elem1 * nelements + elem3) * nelements + elem2) * idxb_max + jjb];
		else if (elem1 == elem3)
		betaj = beta[itriple] + beta[((elem2 * nelements + elem1) * nelements + elem3) * idxb_max + jjb];
		} else
		betaj = beta[itriple];
		if (j2 == j) betaj = 2*beta[itriple];
		else betaj = beta[itriple];
		} else {
		const int jjb = idxb_block[j2][j][j1];
		itriple = ((elem2 * nelements + elem3) * nelements + elem1) * idxb_max + jjb;
		betaj = beta[itriple];
		}
		} else {
		if (j >= j1) {
		const int jjb = idxb_block[j1][j2][j];
		if (j1 == j) {
		if (j2 == j) betaj = 3*beta[jjb];
		else betaj = 2*beta[jjb];
		} else betaj = beta[jjb];
		} else if (j >= j2) {
		const int jjb = idxb_block[j][j2][j1];
		if (j2 == j) betaj = 2*beta[jjb];
		else betaj = beta[jjb];
		} else {
		const int jjb = idxb_block[j2][j][j1];
		betaj = beta[jjb];
		}
		}

		if (!bnorm_flag && j1 > j)
		betaj *= (j1 + 1) / (j + 1.0);

		ylist_r[elem3 * idxu_max + jju] += betaj * ztmp_r;
		ylist_i[elem3 * idxu_max + jju] += betaj * ztmp_i;
		//if (elem3==0 && j ==4 && ma==2 && mb==0)
		// fprintf(screen, "%i %i %i %i %i %i %f %f %f\n", elem1, elem2, elem3, j, j1, j2, betaj, ztmp_r, ylist_r[elem3 * idxu_max + jju]);

		}
		} // end loop over jjz
		}
		@@ -1023,20 +980,6 @@ void SNA::zero_uarraytot(int ielem)
		}
		}

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

		void SNA::addself_uarraytot(double wself_in, int ielem)
		{
		// for (int j = 0; j <= twojmax; j++) {
		// int jju = idxu_block[j];
		// for (int ma = 0; ma <= j; ma++) {
		// ulisttot_r[ielem*idxu_max+jju] = wself_in;
		// ulisttot_i[ielem*idxu_max+jju] = 0.0;
		// jju += j+2;
		// }
		// }
		}

		/* ----------------------------------------------------------------------
		add Wigner U-functions for one neighbor to the total
		------------------------------------------------------------------------- */

src/SNAP/sna.h

+1 −2

Changes for src/SNAP/sna.h: 1 added line, 2 removed lines.
Original line number	Diff line number	Diff line
		@@ -23,7 +23,7 @@
		namespace LAMMPS_NS {

		struct SNA_ZINDICES {
		int j1, j2, j, ma1min, ma2max, mb1min, mb2max, na, nb, jju, ma, mb;
		int j1, j2, j, ma1min, ma2max, mb1min, mb2max, na, nb, jju;
		};

		struct SNA_BINDICES {
		@@ -110,7 +110,6 @@ private:
		void print_clebsch_gordan();
		void init_rootpqarray();
		void zero_uarraytot(int);
		void addself_uarraytot(double, int);
		void add_uarraytot(double, double, double, int, int);
		void compute_uarray(double, double, double,
		double, double, int);

Admin message