port nh fixes to USER-BOCS package (e37ee02e) · Commits · 郑智淋 / lammps

src/USER-BOCS/fix_bocs.cpp

+28 −13

Original line number	Diff line number	Diff line
		@@ -846,7 +846,7 @@ void FixBocs::setup(int vflag)

		if (pstat_flag) {
		double kt = boltz * t_target;
		double nkt = atom->natoms * kt;
		double nkt = (atom->natoms + 1) * kt;

		for (int i = 0; i < 3; i++)
		if (p_flag[i])
		@@ -1508,7 +1508,7 @@ double FixBocs::compute_scalar()
		double volume;
		double energy;
		double kt = boltz * t_target;
		double lkt_press = kt;
		double lkt_press = 0.0;
		int ich;
		if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
		else volume = domain->xprd * domain->yprd;
		@@ -1539,15 +1539,21 @@ double FixBocs::compute_scalar()
		// sum is over barostatted dimensions

		if (pstat_flag) {
		for (i = 0; i < 3; i++)
		if (p_flag[i])
		for (i = 0; i < 3; i++) {
		if (p_flag[i]) {
		energy += 0.5omega_dot[i]omega_dot[i]*omega_mass[i] +
		p_hydro(volume-vol0) / (pdimnktv2p);
		lkt_press += kt;
		}
		}

		if (pstyle == TRICLINIC) {
		for (i = 3; i < 6; i++)
		if (p_flag[i])
		for (i = 3; i < 6; i++) {
		if (p_flag[i]) {
		energy += 0.5omega_dot[i]omega_dot[i]*omega_mass[i];
		lkt_press += kt;
		}
		}
		}

		// extra contributions from thermostat chain for barostat
		@@ -1880,15 +1886,14 @@ void FixBocs::nhc_temp_integrate()

		void FixBocs::nhc_press_integrate()
		{
		int ich,i;
		int ich,i,pdof;
		double expfac,factor_etap,kecurrent;
		double kt = boltz * t_target;
		double lkt_press = kt;

		// Update masses, to preserve initial freq, if flag set

		if (omega_mass_flag) {
		double nkt = atom->natoms * kt;
		double nkt = (atom->natoms + 1) * kt;
		for (int i = 0; i < 3; i++)
		if (p_flag[i])
		omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
		@@ -1912,14 +1917,24 @@ void FixBocs::nhc_press_integrate()
		}

		kecurrent = 0.0;
		for (i = 0; i < 3; i++)
		if (p_flag[i]) kecurrent += omega_mass[i]omega_dot[i]omega_dot[i];
		pdof = 0;
		for (i = 0; i < 3; i++) {
		if (p_flag[i]) {
		kecurrent += omega_mass[i]omega_dot[i]omega_dot[i];
		pdof++;
		}
		}

		if (pstyle == TRICLINIC) {
		for (i = 3; i < 6; i++)
		if (p_flag[i]) kecurrent += omega_mass[i]omega_dot[i]omega_dot[i];
		for (i = 3; i < 6; i++) {
		if (p_flag[i]) {
		kecurrent += omega_mass[i]omega_dot[i]omega_dot[i];
		pdof++;
		}
		}
		}

		double lkt_press = pdof * kt;
		etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];

		double ncfac = 1.0/nc_pchain;

Admin message