From: "Eric B. Isaacs" <[email protected]...>
To: [email protected]
Date: Tue, 25 May 2010 00:53:48 -0700
Subject: [lammps-users] Thermal conductivity
Hi,
Alternatively, is it possible to use the Muller-Plathe method for my system given that there are multiple atom types? (I am concerned by the "the masses of all exchanged atom pairs should be the same to use this fix" in the manual.)
For convenience, robustness, and better statistics, I prefer to do M-P
thermal conductivity calculations swapping between all atoms. Below is
a patch I wrote several months ago that modifies
fix_thermal_conductivity.cpp to allow this while still conserving
energy and momentum. See "Non-equilibrium momentum exchange algorithm
for molecular dynamics simulation of heat flow in multicomponent
systems" (Nieto-Draghi, Molecular Physics 20, 2003) for more info.
--- fix_thermal_conductivity.cpp.any_mass 2010-02-27
06:10:23.000000000 -0500
+++ ../lammps-8May10/src/fix_thermal_conductivity.cpp 2010-01-04
15:25:40.000000000 -0500
@@ -225,8 +225,7 @@
// exchange kinetic energy between the 2 particles
// if I own both particles just swap, else point2point comm of velocities
- double sbuf[4],rbuf[4],vcm[3];
- double eswap = 0.0;
+ double sbuf[3],rbuf[3];
mine[0].proc = mine[1].proc = me;
int ilo = 0;
@@ -240,80 +239,45 @@
MPI_Allreduce(mine,all,2,MPI_DOUBLE_INT,MPI_MINLOC,world);
if (all[0].value == BIG || all[1].value == BIG) continue;
+ all[0].value = -all[0].value;
+ e_exchange += force->mvv2e * (all[0].value - all[1].value);
if (me == all[0].proc && me == all[1].proc) {
i = index_lo[ilo++];
j = index_hi[ihi++];
- sbuf[0] = v[j][0];
- sbuf[1] = v[j][1];
- sbuf[2] = v[j][2];
- if (rmass) sbuf[3] = rmass[j];
- else sbuf[3] = mass[type[j]];
rbuf[0] = v[i][0];
rbuf[1] = v[i][1];
rbuf[2] = v[i][2];
- if (rmass) rbuf[3] = rmass[i];
- else rbuf[3] = mass[type[i]];
- vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]);
- vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]);
- vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]);
- v[j][0] = 2.0 * vcm[0] - sbuf[0];
- v[j][1] = 2.0 * vcm[1] - sbuf[1];
- v[j][2] = 2.0 * vcm[2] - sbuf[2];
- eswap += sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) +
- vcm[1] * (vcm[1] - sbuf[1]) +
- vcm[2] * (vcm[2] - sbuf[2]));
- v[i][0] = 2.0 * vcm[0] - rbuf[0];
- v[i][1] = 2.0 * vcm[1] - rbuf[1];
- v[i][2] = 2.0 * vcm[2] - rbuf[2];
- eswap -= rbuf[3] * (vcm[0] * (vcm[0] - rbuf[0]) +
- vcm[1] * (vcm[1] - rbuf[1]) +
- vcm[2] * (vcm[2] - rbuf[2]));
+ v[i][0] = v[j][0];
+ v[i][1] = v[j][1];
+ v[i][2] = v[j][2];
+ v[j][0] = rbuf[0];
+ v[j][1] = rbuf[1];
+ v[j][2] = rbuf[2];
} else if (me == all[0].proc) {
- j = index_lo[ilo++];
- sbuf[0] = v[j][0];
- sbuf[1] = v[j][1];
- sbuf[2] = v[j][2];
- if (rmass) sbuf[3] = rmass[j];
- else sbuf[3] = mass[type[j]];
- MPI_Sendrecv(sbuf,4,MPI_DOUBLE,all[1].proc,0,
- rbuf,4,MPI_DOUBLE,all[1].proc,0,world,&status);
- vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]);
- vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]);
- vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]);
- v[j][0] = 2.0 * vcm[0] - sbuf[0];
- v[j][1] = 2.0 * vcm[1] - sbuf[1];
- v[j][2] = 2.0 * vcm[2] - sbuf[2];
- eswap -= sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) +
- vcm[1] * (vcm[1] - sbuf[1]) +
- vcm[2] * (vcm[2] - sbuf[2]));
+ i = index_lo[ilo++];
+ sbuf[0] = v[i][0];
+ sbuf[1] = v[i][1];
+ sbuf[2] = v[i][2];
+ MPI_Sendrecv(sbuf,3,MPI_DOUBLE,all[1].proc,0,
+ rbuf,3,MPI_DOUBLE,all[1].proc,0,world,&status);
+ v[i][0] = rbuf[0];
+ v[i][1] = rbuf[1];
+ v[i][2] = rbuf[2];
} else if (me == all[1].proc) {
j = index_hi[ihi++];
sbuf[0] = v[j][0];
sbuf[1] = v[j][1];
sbuf[2] = v[j][2];
- if (rmass) sbuf[3] = rmass[j];
- else sbuf[3] = mass[type[j]];
- MPI_Sendrecv(sbuf,4,MPI_DOUBLE,all[0].proc,0,
- rbuf,4,MPI_DOUBLE,all[0].proc,0,world,&status);
- vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]);
- vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]);
- vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]);
- v[j][0] = 2.0 * vcm[0] - sbuf[0];
- v[j][1] = 2.0 * vcm[1] - sbuf[1];
- v[j][2] = 2.0 * vcm[2] - sbuf[2];
- eswap += sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) +
- vcm[1] * (vcm[1] - sbuf[1]) +
- vcm[2] * (vcm[2] - sbuf[2]));
+ MPI_Sendrecv(sbuf,3,MPI_DOUBLE,all[0].proc,0,
+ rbuf,3,MPI_DOUBLE,all[0].proc,0,world,&status);
+ v[j][0] = rbuf[0];
+ v[j][1] = rbuf[1];
+ v[j][2] = rbuf[2];
}
}
- // tally energy exchange from all swaps