Dear LAMMPS-users,
My goal is to establish a temperature gradient along the z-axis in a simulation of SPC/E water using PPPM in an orthorhombic box. I add heat to a hot region and remove it from a cold region. In order to analyse the data, I use ave/spatial to calculate the average kinetic energy per atom, average atom type and average number density. It actually works and looks like I reached a steady state (v-shaped temperature profile that does not change anymore). However, I don’t think that my way is the best way of doing it, and thus the questions:
1.) Temperature: Can avg/spatial directly calculate the temperature for rigid molecules, like water?
From ave/spatial I know the average number of atoms in the bin, <n_a>, the average kinetic energy per atom, <ek_a>, and the average atom type, , where (t(O) = 2 and t(H) = 1). Since the average type is around 1.33, I thought it would be justified to calculate T in the following way for every bin: T = 3 <ek_a> / (6/2 k_B), because a water molecule has 6 degrees of freedom.
However, as the ratio is not always exactly 1.33, I am making a slight mistake when calculating the temperature this way.
2.) Density: Can avg/spatial directly calculate the number densities of individual atom types, e.g rho_O and rho_H?
At the moment I use rho(O) = * ( - t(O)) / (t(O)-t(H)) and rho(H) = - rho(O) and apply the formula to every bin. Here is the average number density in the bin from ave/spatial. Again I have to do post processing and couldn’t find a better/direct way of calculating the densities although I am sure there must be one.
Any help would be much appreciated,
Peter
My input file:
Dear LAMMPS-users,
My goal is to establish a temperature gradient along the z-axis in a
simulation of SPC/E water using PPPM in an orthorhombic box. I add heat to a
hot region and remove it from a cold region. In order to analyse the data, I
use ave/spatial to calculate the average kinetic energy per atom, average
atom type and average number density. It actually works and looks like I
reached a steady state (v-shaped temperature profile that does not change
anymore). However, I don't think that my way is the best way of doing it,
and thus the questions:
1.) Temperature: Can avg/spatial directly calculate the temperature for
rigid molecules, like water?
From ave/spatial I know the average number of atoms in the bin, <n_a>, the
average kinetic energy per atom, <ek_a>, and the average atom type, <t>,
where (t(O) = 2 and t(H) = 1). Since the average type is around 1.33, I
thought it would be justified to calculate T in the following way for every
bin: T = 3 <ek_a> / (6/2 k_B), because a water molecule has 6 degrees of
freedom.
However, as the ratio is not always exactly 1.33, I am making a slight
mistake when calculating the temperature this way.
how about you compute spatial averages for oxygen and hydrogen atoms
independently?
since a rigid water is missing 3 DOFs vs. a flexible one, you can
subtract 1 DOF per atom.
2.) Density: Can avg/spatial directly calculate the number densities of
individual atom types, e.g rho_O and rho_H?
At the moment I use rho(O) = <rho> * (<t> - t(O)) / (t(O)-t(H)) and rho(H) =
<rho> - rho(O) and apply the formula to every bin. Here <rho> is the average
number density in the bin from ave/spatial. Again I have to do post
processing and couldn't find a better/direct way of calculating the
densities although I am sure there must be one.
same as before, why not generate a compute for each atom type separately?
if you absolutely want to do this in one sweep, you can look at the
grmask() variable function that allows you to mask out contributions
from atoms that are not part of a specific group.
of course this is not as easily applicable to more complex rigid
structures, but for water it should work.
axel.
Not 100% sure but maybe you need to subtract off a centre-of-mass term? The temperature gradient might cause the whole system to start moving in one direction, which would lead to the temperature calculated from the (uncorrected) KE not quite lining up with the expected T. I've just started doing similar work and I know we do this kind of correction, so maybe this is an issue for you too.
Chris.
Dear Axel/Chris,
thanks for your quick reply.
ad 1) Temperature
Sorry I think I misinterpreted your initial post. I thought your T calculation from KE was a bit off and you were trying to figure out why. By center-of-mass term I meant the c.o.m of the whole simulation box, not the molecular center of mass. If the system as a whole is moving, you might need to remove this part of the velocity before converting average molecular/atomic speed to temperature. But as you say the KE contribution from some motion of your whole system wouldn’t have any relation to the d.o.f. issues you are concerned with. So… carry on I guess!
@Chris>
Actually, now I understand what you mean and I am glad that you pointed it out. I will check how the centre of mass of the whole system evolves in time and see if it requires a correction. Thanks!
@Axel>
I think now I understand what you meant. There is no need to use “compute_modify”, which I wanted to use initially. I modified the example and posted it for completeness. Thank you!
Best regards
Peter
Input file:
There's also fix momentum, if you want to stop your box from accumulating net momentum.
@Niall>
Thanks for the tip, it was very useful.
Peter