(539fb97f) · Commits · 郑智淋 / lammps

src/fix_langevin.cpp

+89 −37

Original line number	Diff line number	Diff line
		@@ -31,6 +31,7 @@
		using namespace LAMMPS_NS;

		enum{NOBIAS,BIAS};
		enum{MASS,RMASS};

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

		@@ -105,22 +106,25 @@ int FixLangevin::setmask()

		void FixLangevin::init()
		{
		if (atom->mass == NULL)
		error->all("Cannot use fix langevin without per-type mass defined");

		// set force prefactors

		if (atom->mass) {
		for (int i = 1; i <= atom->ntypes; i++) {
		gfactor1[i] = -atom->mass[i] / t_period / force->ftm2v;
		gfactor2[i] = sqrt(atom->mass[i]) *
		sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) / force->ftm2v;
		sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
		force->ftm2v;
		gfactor1[i] *= 1.0/ratio[i];
		gfactor2[i] *= 1.0/sqrt(ratio[i]);
		}
		}

		if (temperature && temperature->tempbias) which = BIAS;
		else which = NOBIAS;

		if (atom->mass) massflag = MASS;
		else massflag = RMASS;

		if (strcmp(update->integrate_style,"respa") == 0)
		nlevels_respa = ((Respa *) update->integrate)->nlevels;
		}
		@@ -158,6 +162,7 @@ void FixLangevin::post_force(int vflag)

		// apply damping and thermostat to appropriate atoms

		if (massflag == MASS) {
		if (which == NOBIAS) {
		for (int i = 0; i < nlocal; i++) {
		if (mask[i] & groupbit) {
		@@ -170,7 +175,7 @@ void FixLangevin::post_force(int vflag)
		}

		// invoke temperature since some computes require it to remove bias
		// test on v = 0 since some computes mask non-participating atoms via v = 0
		// test v = 0 since some computes mask non-participating atoms via v = 0

		} else if (which == BIAS) {
		double tmp = temperature->compute_scalar();
		@@ -189,7 +194,51 @@ void FixLangevin::post_force(int vflag)
		temperature->restore_bias(v[i]);
		}
		}
		}

		} else {
		double *rmass = atom->rmass;
		double boltz = force->boltz;
		double dt = update->dt;
		double mvv2e = force->mvv2e;
		double ftm2v = force->ftm2v;

		if (which == NOBIAS) {
		for (int i = 0; i < nlocal; i++) {
		if (mask[i] & groupbit) {
		gamma1 = -rmass[i] / t_period / ftm2v;
		gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
		gamma1 *= 1.0/ratio[type[i]];
		gamma2 = 1.0/sqrt(ratio[type[i]]) tsqrt;
		f[i][0] += gamma1v[i][0] + gamma2(random->uniform()-0.5);
		f[i][1] += gamma1v[i][1] + gamma2(random->uniform()-0.5);
		f[i][2] += gamma1v[i][2] + gamma2(random->uniform()-0.5);
		}
		}

		// invoke temperature since some computes require it to remove bias
		// test v = 0 since some computes mask non-participating atoms via v = 0

		} else if (which == BIAS) {
		double tmp = temperature->compute_scalar();

		for (int i = 0; i < nlocal; i++) {
		if (mask[i] & groupbit) {
		gamma1 = -rmass[i] / t_period / ftm2v;
		gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
		gamma1 *= 1.0/ratio[type[i]];
		gamma2 = 1.0/sqrt(ratio[type[i]]) tsqrt;
		temperature->remove_bias(i,v[i]);
		if (v[i][0] != 0.0)
		f[i][0] += gamma1v[i][0] + gamma2(random->uniform()-0.5);
		if (v[i][1] != 0.0)
		f[i][1] += gamma1v[i][1] + gamma2(random->uniform()-0.5);
		if (v[i][2] != 0.0)
		f[i][2] += gamma1v[i][2] + gamma2(random->uniform()-0.5);
		temperature->restore_bias(v[i]);
		}
		}
		}
		}
		}

		@@ -211,12 +260,15 @@ void FixLangevin::reset_target(double t_new)

		void FixLangevin::reset_dt()
		{
		if (atom->mass) {
		for (int i = 1; i <= atom->ntypes; i++) {
		gfactor2[i] = sqrt(atom->mass[i]) *
		sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) / force->ftm2v;
		sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
		force->ftm2v;
		gfactor2[i] *= 1.0/sqrt(ratio[i]);
		}
		}
		}

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

src/fix_langevin.h

+1 −1

Original line number	Diff line number	Diff line
		@@ -32,7 +32,7 @@ class FixLangevin : public Fix {
		int modify_param(int, char **);

		private:
		int which;
		int which,massflag;
		double t_start,t_stop,t_period;
		double gfactor1,gfactor2,*ratio;

Admin message