added an interface for the inverse c.o.b. matrix to UEF_utils (c3bf7d09) · Commits · 郑智淋 / lammps

src/USER-UEF/uef_utils.cpp

+24 −5

Original line number	Diff line number	Diff line
		@@ -55,6 +55,7 @@ UEFBox::UEFBox()
		{
		l[k][j] = l0[k][j];
		r[j][k]=(j==k);
		ri[j][k]=(j==k);
		}

		// get the initial rotation and upper triangular matrix
		@@ -119,6 +120,14 @@ void UEFBox::get_rot(double x[3][3])
		x[k][j]=rot[k][j];
		}

		// get inverse change of basis matrix
		void UEFBox::get_inverse_cob(double x[3][3])
		{
		for (int k=0;k<3;k++)
		for (int j=0;j<3;j++)
		x[k][j]=ri[k][j];
		}

		// diagonal, incompressible deformation
		void UEFBox::step_deform(const double ex, const double ey)
		{
		@@ -167,19 +176,30 @@ bool UEFBox::reduce()
		if (f2 < 0) for (int k=0;k<-f2;k++) mul_m2 (a2i,r0);

		// robust reduction to the box defined by Dobson
		double lc[3][3];
		for (int k=0;k<3;k++)
		{
		double eps = exp(theta[0]w1[k]+theta[1]w2[k]);
		l[k][0] = eps*l0[k][0];
		l[k][1] = eps*l0[k][1];
		l[k][2] = eps*l0[k][2];
		lc[k][0] = eps*l0[k][0];
		lc[k][1] = eps*l0[k][1];
		lc[k][2] = eps*l0[k][2];
		}

		// further reduce the box using greedy reduction and check
		// if it changed from the last step using the change of basis
		// matrices r and r0
		int ri[3][3];

		greedy(l,r,ri);
		mul_m2(r,ri);

		// multiplying the inverse by the old change of basis matrix gives
		// the inverse of the transformation itself (should be identity if
		// no reduction takes place). This is used for image flags only.

		mul_m1(ri,r0);

		rotation_matrix(rot,lrot, l);
		return !mat_same(r,r0);
		}
		@@ -197,7 +217,6 @@ void UEFBox::set_strain(const double ex, const double ey)
		l[k][1] = eps*l0[k][1];
		l[k][2] = eps*l0[k][2];
		}
		int ri[3][3];
		greedy(l,r,ri);
		rotation_matrix(rot,lrot, l);
		}

+9 −8

Original line number	Diff line number	Diff line
		@@ -27,12 +27,13 @@ class UEFBox
		bool reduce();
		void get_box(double[3][3], double);
		void get_rot(double[3][3]);
		void get_inverse_cob(double[3][3]);
		private:
		double l0[3][3]; // initial basis
		double w1[3],w2[3],winv[3][3];//omega1 and omega2 (spectra of automorphisms)
		double theta[2];
		double l[3][3], rot[3][3], lrot[3][3];
		int r[3][3],a1[3][3],a2[3][3],a1i[3][3],a2i[3][3];
		int r[3][3],ri[3][3],a1[3][3],a2[3][3],a1i[3][3],a2i[3][3];
		};


		@@ -112,10 +113,10 @@ void transpose(T m[3][3])
		}
		}

		template<typename T>
		void mul_m1(T m1[3][3], const T m2[3][3])
		template<typename T1,typename T2>
		void mul_m1(T1 m1[3][3], const T2 m2[3][3])
		{
		T t[3][3];
		T1 t[3][3];
		for (int k=0;k<3;k++)
		for (int j=0;j<3;j++)
		t[k][j]=m1[k][j];
		@@ -125,10 +126,10 @@ void mul_m1(T m1[3][3], const T m2[3][3])
		m1[k][j] = t[k][0]m2[0][j] + t[k][1]m2[1][j] + t[k][2]*m2[2][j];
		}

		template<typename T>
		void mul_m2(const T m1[3][3], T m2[3][3])
		template<typename T1, typename T2>
		void mul_m2(const T1 m1[3][3], T2 m2[3][3])
		{
		T t[3][3];
		T2 t[3][3];
		for (int k=0;k<3;k++)
		for (int j=0;j<3;j++)
		t[k][j]=m2[k][j];