set appropriate flags for computing a vector in constructor of fix restrain (69605295) · Commits · 郑智淋 / lammps

src/fix_restrain.cpp

+25 −3

Original line number	Diff line number	Diff line
		@@ -53,6 +53,9 @@ FixRestrain::FixRestrain(LAMMPS lmp, int narg, char *arg) :
		scalar_flag = 1;
		global_freq = 1;
		extscalar = 1;
		vector_flag = 1;
		size_vector = 3;
		extvector = 1;
		respa_level_support = 1;
		ilevel_respa = 0;

		@@ -189,6 +192,10 @@ void FixRestrain::post_force(int vflag)
		{
		energy = 0.0;

		ebond = 0.0;
		eangle = 0.0;
		edihed = 0.0;

		for (int m = 0; m < nrestrain; m++)
		if (rstyle[m] == BOND) restrain_bond(m);
		else if (rstyle[m] == ANGLE) restrain_angle(m);
		@@ -271,7 +278,8 @@ void FixRestrain::restrain_bond(int m)
		if (r > 0.0) fbond = -2.0*rk/r;
		else fbond = 0.0;

		energy += rk*dr;
		ebond += rk*dr;
		energy += ebond;

		// apply force to each of 2 atoms

		@@ -378,7 +386,8 @@ void FixRestrain::restrain_angle(int m)
		dtheta = acos(c) - target[m];
		tk = k * dtheta;

		energy += tk*dtheta;
		eangle += tk*dtheta;
		energy += eangle;

		a = -2.0 * tk * s;
		a11 = a*c / rsq1;
		@@ -558,7 +567,8 @@ void FixRestrain::restrain_dihedral(int m)
		df1 *= -mult[m];
		p += 1.0;

		energy += k * p;
		edihed += k * p;
		energy += edihed;

		fg = vb1xvb2xm + vb1yvb2ym + vb1z*vb2zm;
		hg = vb3xvb2xm + vb3yvb2ym + vb3z*vb2zm;
		@@ -635,3 +645,15 @@ double FixRestrain::compute_scalar()
		MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
		return energy_all;
		}

		/* ----------------------------------------------------------------------
		return individual energy contributions
		------------------------------------------------------------------------- */

		double FixRestrain::compute_vector(int n)
		{
		if (n == 0) return ebond;
		if (n == 1) return eangle;
		if (n == 2) return edihed;
		return 0.0;
		}

src/fix_restrain.h

+2 −0

Original line number	Diff line number	Diff line
		@@ -36,6 +36,7 @@ class FixRestrain : public Fix {
		void post_force_respa(int, int, int);
		void min_post_force(int);
		double compute_scalar();
		double compute_vector(int);

		private:
		int ilevel_respa;
		@@ -46,6 +47,7 @@ class FixRestrain : public Fix {
		double kstart,kstop,*target;
		double cos_target,sin_target;
		double energy,energy_all;
		double ebond,eangle,edihed;

		void restrain_bond(int);
		void restrain_angle(int);

Admin message