Dear LAMMPS users,
I want to calculate the heat flux through a group of water molecules. A rigid water model, TIP4P, is used in my simulations. Angular velocities and torque are supposed to be considered for rigid molecules. The LAMMPS command “compute heat/flux” computes the heat flux based on contributions from atoms rather than molecules. I wonder whether it is appropriate to use this command to calculate the heat flux through a group of water molecules.
Regards,
Han Hu
Dear LAMMPS users,
I want to calculate the heat flux through a group of water molecules. A
rigid water model, TIP4P, is used in my simulations. Angular velocities and
torque are supposed to be considered for rigid molecules. The LAMMPS command
"compute heat/flux" computes the heat flux based on contributions from
atoms rather than molecules. I wonder whether it is appropriate to use this
command to calculate the heat flux through a group of water molecules.
please note that a rigid molecule integrator separates the total motion
into a motion of the center of mass and the rotation around that center
of mass. you can project that motion back into the motion of individual
atoms that would comprise such a rigid molecule. thus i don't see a
principal problem with that.
*however*, there have been repeated discussions on the mailing list
about using fix heat/flux on water and it has been pointed out that
you cannot use it for as long as the correct representation of your
system requires the use of long-range electrostatics, since those
cannot be easily decomposed into per atom energy contributions.
i suggest you dig into the mailing list archives and read through
some of the discussions in order to see, if this scenario applies
to your system as well and what alternatives are available, if yes.
cheers,
axel.
Green-Kubo and compute heat/flux are equilibrium
methods. I think they can be applied to any system,
regardless of whether there are molecules are not.
But as Axel says, if there is long-range Coulombics
then it isn't included in the per-atom virial used
by compute heat/flux.
Steve
Dear Axel and Steve,
Thanks. In my simulation, the long-range Coulombics interaction between water molecules is taken into consideration by using the command “kspace pppm/tip4p”. In the previous discussion of calculating heat flux, the LAMMPS user German provided a modified compute heat/flux code which worked with pair potentials and Ewald summation. I will try this code first. I am grateful for German’s contribution.
I have one more question. Should the bonds be accounted for in the calculation of potential per atom? In order to compute the heat flux using the Green-Kubo formula, the potential per atom should be calculated first using the command “compute pe/atom”. In the command “compute pe/atom”, we could choose which contribution is considered. In my simulations of water molecules, TIP4P pair potential is applied. Harmonic bond and harmonic angle is applied for intramolecular interaction. Should the bond energy be accounted for in the calculation of potential energy per atom? In my opinion, the bonds are the intramolecular interaction so it should be be taken into consideration in the computation of pe/atom. What do you think? Thanks a lot.
Sincerely,
Han Hu
Dear Axel and Steve,
Thanks. In my simulation, the long-range Coulombics interaction between
water molecules is taken into consideration by using the command "kspace
pppm/tip4p". In the previous discussion of calculating heat flux, the LAMMPS
user German provided a modified compute heat/flux code which worked with
pair potentials and Ewald summation. I will try this code first. I am
grateful for German's contribution.
that won't work on TIP4P water, though. all /tip4p styles compute the
position of the charge point M of the TIP4P water on the fly. hence using
any plain long-range style will produce incorrect results.
I have one more question. Should the bonds be accounted for in the
calculation of potential per atom? In order to compute the heat flux using
the Green-Kubo formula, the potential per atom should be calculated first
using the command "compute pe/atom". In the command "compute pe/atom", we
could choose which contribution is considered. In my simulations of water
molecules, TIP4P pair potential is applied. Harmonic bond and harmonic
angle is applied for intramolecular interaction. Should the bond energy be
for rigid TIP4P water, the bond/angle potentials are only needed
to get the equilibrium distrance/angle.
accounted for in the calculation of potential energy per atom? In my
opinion, the bonds are the intramolecular interaction so it should be be
taken into consideration in the computation of pe/atom. What do you
think? Thanks a lot.
if you use fix shake to keep the water rigid, bond/angle
forces of the water molecules are not computed.
axel.
Dear Axel,
I’ve read German’s codes and try to modify the codes following his example. But I encounter a problem that the potential energy and virial are calculated at grids rather than atoms. In German’s code, the long-range Coulombic interaction is calculated with a standard Ewald summation(using ewald.cpp). The potential energy and virial are calculated at atoms.
I am using pppm/tip4p ( and some functions defined in pppm.cpp are also called). In the code pppm,cpp, the energy and virial are calculated at the grids. The total potential energy is calculated by summing the potential energy of each grid through all processors. The potential energy and virial of each atom are not calculated at all. Do you think it possible to calculate potential energy per atom and virial per atom using current potential energy per grid and virial per grid?
Thanks,
Han Hu
Dear Axel,
I've read German's codes and try to modify the codes following his example.
But I encounter a problem that the potential energy and virial are
calculated at grids rather than atoms. In German's code, the long-range
Coulombic interaction is calculated with a standard Ewald summation(using
ewald.cpp). The potential energy and virial are calculated at atoms.I am using pppm/tip4p ( and some functions defined in pppm.cpp are also
called). In the code pppm,cpp, the energy and virial are calculated at the
grids. The total potential energy is calculated by summing the potential
energy of each grid through all processors. The potential energy and virial
of each atom are not calculated at all. Do you think it possible to
calculate potential energy per atom and virial per atom using current
potential energy per grid and virial per grid?
this topic has come up quite a few times before.
so far nobody on this list know how to handle this.
your only (potentially) viable option would be to
write a similarly modified ewald/tip4p kspace style.
axel.
Do you think it possible to
calculate potential energy per atom and virial per atom using current
potential energy per grid and virial per grid?
How? Since a grid point contains contributions to many
nearby atoms.
Steve
Han…
I suggest you use Ewald methodology for the long range calculations. This methodology has been reported multiple times in the literature of polymers, ionic salts and other systems. Please see one of our papers (Varshney et al. Polymer 50, 2009, 3378-3385) and the referencs there in (refs 27,28). I used the same methodology (however I wrote my own post processing code) before German provided his (with works with lammps). I know ewald is slower but better to go with something that works.
Best Regards,
Vikas
Dear Steve,
Thanks. To tell the truth, I don’t know how to make it yet. The code ‘pppm.cpp’ is complex and I have not understood every function.
1). A charge density function has been established to map to the charges of the atoms to the grid. Could this function be used to interpolate the potential and virial stress from the grid back to the atoms?
2). The electric filed and forces are obtained by interpolating from grid to atoms. Could the similar method be used to interpolate the potential and virial stress from the grid back to the atoms? Or could we just calculate the potential or virial stress using the forces that have been calculated.
Thanks a lot.
Best wishes,
Han Hu
Dear Steve,
Thanks. To tell the truth, I don't know how to make it yet. The code
'pppm.cpp' is complex and I have not understood every function.1). A charge density function has been established to map to the charges of
the atoms to the grid. Could this function be used to interpolate the
potential and virial stress from the grid back to the atoms?
no.
2). The electric filed and forces are obtained by interpolating from grid to
atoms. Could the similar method be used to interpolate the potential and
virial stress from the grid back to the atoms? Or could we just calculate
the potential or virial stress using the forces that have been calculated.
no.
axel.
Dear Alex,
Thanks. Then my last hope will be to write a code of standard ewald summation for TIP4P water. The approximate procedures are shown below.
(1) compute the position of M site of each molecule as the position of the negative charge;
(2) use ewald summation to calculate the potential and virial stress. redistribute the potential and virial to the physical atoms.
(3 ) a 1d array *eatom (in the class kspace) is defined to store the potential per atom and a 2d array **vatom is defined to store the virial per atom.
(4) the kspace->eatom and kspace->vatom will be used in the calculation of pe/atom and stress/atom.
then the heat flux through a group of TIP4P water with long-range Coulombic interaction could be calculated.
Do you think it possible to redistribute the virial stress tensor back to the physical atoms?
Sincerely,
Han Hu
Dear Vikas,
Thanks a lot. I have browsed your paper and ref.27 and ref.28 of your paper. Now I am determined to use standard Ewald summation rather than the pppm method.
In the email, you mentioned that you used your post processing code to compute heat flux. Could you briefly tell me your idea about your post processing code because I really have an idea about calculating heat flux with a post progressing code. In the formula of calculating heat flux (e.g., Eq. (6) in your paper), the velocities and the stress of each atom are needed. The velocities could be output in the dump files. But how could we output the stress tensor of each atom? And I think a large amount of storage will be needed. In fact, I prefer to write a post processing code rather than modifying the functions of LAMMPS because I cannot modify the source files LAMMPS in the supercomputers and most of my cases are performed in supercomputers
Thank you very much!
Sincerely,
Han Hu
Dear Han,
As you see from equation 6, the stress tensor is not for a single atom but for a pair S_ij. As you will see from equation 7, this is calculated in reciprocal space. What you need is just the coordinates of the atoms and the charges. Rest is just coding to calculate S^{alpha,beta}{i,j}. I will email the code that takes care of Sij from my other email account (I am having problems downloading the attachment on my machine). Please go through it and see if that helps. It will see a lot of familiarity with ewald.cpp in lammps as I used their framework. It essentially calculates equation 7. I hope this helps.
Best Regards,
Vikas
Dear Vikas,
Many thanks. I will try to run the code later.
Bests,
Han Hu