update compute heat/flux to work with compute centroid/stress/atom (9301fc6b) · Commits · 郑智淋 / lammps

src/compute_heat_flux.cpp

+45 −14

Original line number	Diff line number	Diff line
		@@ -65,9 +65,10 @@ ComputeHeatFlux::ComputeHeatFlux(LAMMPS lmp, int narg, char *arg) :
		error->all(FLERR,"Compute heat/flux compute ID does not compute ke/atom");
		if (modify->compute[ipe]->peatomflag == 0)
		error->all(FLERR,"Compute heat/flux compute ID does not compute pe/atom");
		if (modify->compute[istress]->pressatomflag == 0)
		if (modify->compute[istress]->pressatomflag != 1
		&& modify->compute[istress]->pressatomflag != 2)
		error->all(FLERR,
		"Compute heat/flux compute ID does not compute stress/atom");
		"Compute heat/flux compute ID does not compute stress/atom or centroid/stress/atom");

		vector = new double[size_vector];
		}
		@@ -137,6 +138,35 @@ void ComputeHeatFlux::compute_vector()
		double jv[3] = {0.0,0.0,0.0};
		double eng;

		// heat flux via centroid atomic stress
		if (c_stress->pressatomflag == 2) {
		for (int i = 0; i < nlocal; i++) {
		if (mask[i] & groupbit) {
		eng = pe[i] + ke[i];
		jc[0] += eng*v[i][0];
		jc[1] += eng*v[i][1];
		jc[2] += eng*v[i][2];
		// stress[0]: rijx*fijx
		// stress[1]: rijy*fijy
		// stress[2]: rijz*fijz
		// stress[3]: rijx*fijy
		// stress[4]: rijx*fijz
		// stress[5]: rijy*fijz
		// stress[6]: rijy*fijx
		// stress[7]: rijz*fijx
		// stress[8]: rijz*fijy
		// jv[0] = rijx fijx vjx + rijx fijy vjy + rijx fijz vjz
		jv[0] -= stress[i][0]v[i][0] + stress[i][3]v[i][1] +
		stress[i][4]*v[i][2];
		// jv[1] = rijy fijx vjx + rijy fijy vjy + rijy fijz vjz
		jv[1] -= stress[i][6]v[i][0] + stress[i][1]v[i][1] +
		stress[i][5]*v[i][2];
		// jv[2] = rijz fijx vjx + rijz fijy vjy + rijz fijz vjz
		jv[2] -= stress[i][7]v[i][0] + stress[i][8]v[i][1] +
		stress[i][2]*v[i][2];
		}
		}
		} else {
		for (int i = 0; i < nlocal; i++) {
		if (mask[i] & groupbit) {
		eng = pe[i] + ke[i];
		@@ -151,6 +181,7 @@ void ComputeHeatFlux::compute_vector()
		stress[i][2]*v[i][2];
		}
		}
		}

		// convert jv from stress*volume to energy units via nktv2p factor

Admin message