Binning heat capacity

LAMMPS LAMMPS Beginners
1

Hi everybody,

I need to calculate the heat capacity of water in a series of cylindrical bins.
This calculation is supposed to give information on the local heat capacity, thus the thermodynamic variables should be calculated for each bin where I use equation 2.82 from Allen-Tildesley as follows:

C = (3/2) (N k_B) / ( 1 - ( <KE^2> - < KE >^2 )*2/(3 N k^2_B T^2) )

I read some of the previous threads and came up with an approach but I am not sure what I calculate is correct or not. Here is the way that I do binning:

compute dvd all chunk/atom bin/cylinder x lower 50 40 40 0.0 40.0 20 discard yes

variable np equal atoms # number of atoms
compute ke_peratom all ke/atom # KE per atom: ke_i, i=1, 2, …, N where N is the number of atoms

variable meanke atom c_ke_peratom/v_np # KE/N per atom: ke_i/N i=1, 2, …, N where N is the number of atoms
variable meanke2 atom v_meanke^2 # (KE/N)^2 per atom: (ke_i/N)^2 i=1, 2, …, N where N is the number of atoms

variable meank2e1 atom c_ke_peratom^2 # KE^2 per atom: ke^2_i, i=1, 2, …, N where N is the number of atoms
variable meank2e2 atom v_meank2e1/v_np # KE^2/N per atom: ke^2_i/N, i=1, 2, …, N where N is the number of atoms
variable kestd_v1 atom v_meank2e1^2/v_np-v_meanke2 # KE^2/N - (KE/N)^2 per atom: ke^2_i/N - (sum_i ke_i/N)^2

variable kb equal 0.00198722886 # kcal/mol/K
variable np1 atom density/mass
variable npc atom v_np1/mass

variable hc atom (1.5 * v_kb * v_npc)/(1-v_kestd_v1 * 2/(3 * v_kb^2 * temp^2 * v_npc))
fix heat_capacity_profile all ave/chunk 10 5 100 dvd v_hc file heat_capacity.profile # in kcal/mol/K

I think I do not calculate the number of atoms per bin correctly.
I appreciate it if anyone could take a look at this and give me some advice.

2

Hey,

I tried to read a bit about the chunk-related commands you are using (I never used them myself) and checked the equation you are referring to in the book. I never calculated Cv using this equation but I have calculated it using a similar expression, derived for the dynamics in the NVT ensemble (pages 65 and 66 of Chandler’s book “Introduction to Modern Statistical Mechanics”). I think there are two problems: the first is that the variables you are defining have too many divisions by N on it. For example, if I understood correctly, your variables np, ke_peratom, meank2e1, meank2e2 and kestd_v1 would have the values shown below at each time instant (and the fact that they are being used within an ave/chunk command would make sure only atoms within a given cylinder bin are considered in the calculation for each corresponding bin, so this should be okay). So you can see that, firstly, you have too many divisions by N compared to the expression from which Cv is calculated from. Secondly, I think the quantity <delta_H^2> involved in Allen-Tildesley formula would require a different calculation (but I could be wrong since I didnt read it very carefully and am thinking comparetively to the expression I have used). As shown below for <delta_X^2> (where X is just a general variable I am inventing to make the point), you would need to do the square of <X>, which requires determining the average first.

image
image1920×1725 109 KB

I think that the most proper way to determine the quantity you want to determine would be to first compute, separately, the value of <delta_H^2> (which is an extensive quantity) and the average value of N at each bin and use them in the formula. From quickly looking into the book, I didnt understand what H is supposed to mean in that formula, but I suppose it is the kinetic energy, given the calculation you are doing. You could define a variable that does the sum of ke_peratomin a bin for example and calculate the due average values you need (which should be extensive quantities), shown for the variable X I invented, at each bin. At the same time, you could calculate the average value of N at each bin. And then, as I said, you can plug these values in the expression for Cv derived in the NVE ensemble, which is the one you seem to be using.

3

In addition to all those specific concerns, there is also the more general question whether there is a meaningful property to be computed at such a small scale.

Take for example temperature. While you can compute “a temperature” even for very small systems, its meaning can be very limited and is based on the assumption of equipartitioning, i.e. that every degree of freedom contains the same amount of energy. This really is only true in the limit of large systems and only if your system is in equilibrium, you can average over time in addition to averaging over volume. Yet still for some properties there are finite size effects and thus the interpretation of data computed at such level of detail difficult if not impossible. With heat capacities being a far more complicated to determine property compared to temperature, I would be much more concerned about what you can learn from this, even if the math works out.

4

Many thanks for checking the formula and suggestions. To be honest, I am totally confused if ave/chunk lumps the per-atoms variables and algebraically sum or it multiplies a factor of 1/N to the result. I tried various ways to test, and came to the conclusion that it just sums the numbers without dividing the result by N. So if this is the case, for example, considering ki_peratom to be a series of k1/N, k2/N, … then at the end, ave/chunk gives (k1/N + k2/N + k3/N + …) = sum_i ki/N that is the average . By this logic, I came up with that script. I, however, may be wrong. Then, I agree, I should remove 1/N. Any comment on this? And yes, KE is the kinetic energy.

5

I am glad to see we share the same concerns as given in your general comments and the validity of such model calculation at small scales. What I tend to observe, is how strong the fluctuations deviate the liquids away from simplistic and ideal “text-book” models. The systems under study are strong thermal spikes induced by the passage of high-energy charged particles such as cosmic rays. The strong expansion of thermal spikes makes the formation of nano-cavities, similar to what I showed in these videos (link is given below). By calculation of “local” heat capacity, I tend to examine how the behavior of shock waves, in a more realistic system, deviates from the adiabatic expansion of ideal gas - a text-book model known as Sedov solutions, a model developed initially for supernova instabilities in astrophysics. In this model, the inside of the cavity behaves like an ideal gas, and the shock wave propagates adiabatically. I want to know if in the atomistic models, this is indeed the case, or if we must deal with a more complex system and fully calculate hydrodynamic equations such as Navier-stokes equations. Finally, here is the link to the system of interest: (289) 2023 08 09 07 52 30 167 2SWs v1 FB - YouTube

6

So, in the manual I found this:

Each input value can also be averaged over the atoms in each chunk. The way the averaging is done across the Nrepeat time steps to produce output on the Nfreq time steps, and across multiple Nfreq outputs, is determined by the norm and ave keyword settings, as discussed below.

The norm keyword affects how averaging is done for the per-chunk values that are output every Nfreq time steps.

It the norm setting is all, which is the default, a chunk value is summed over all atoms in all Nrepeat samples, as is the count of atoms in the chunk. The averaged output value for the chunk on the Nfreq time steps is Total-sum / Total-count. (…)

If the norm setting is sample, the chunk value is summed over atoms for each sample, as is the count, and an “average sample value” is computed for each sample (i.e., Sample-sum / Sample-count). The output value for the chunk on the Nfreq time steps is the average of the Nrepeat “average sample values” (i.e., the sum of Nrepeat “average sample values” divided by Nrepeat).

If the norm setting is none, a similar computation as for the sample setting is done, except the individual “average sample values” are “summed sample values”. A summed sample value is simply the chunk value summed over atoms in the sample, without dividing by the number of atoms in the sample.

So if my understanding is correct, I think the keyword all will sum the value of the quantity (let’s say KE of each atom) over all atoms in a chunk (Ni) and over all the L microstates concerned in the fix ave/chunk and divide the final value by Y, where Y would be the summation of Ni from i = 1 to i = L. The number of microstates concerned by the fix ave/chunk (i.e., L) would be determined by Nevery, Nrepeat and Nfreq together.

The keyword sample will sum the value of the quantity in a given microstate i concerned by the ave/chunk and divide it by Ni (ie., number of atoms in teh chunk in that given microstate) and then average this (which is an average on its own) over all the L microstates concerned in the fix ave/chunk. Finally, the keyword none would be the average kinetic energy (extensive) over the L microstates concerned in the ave/chunk. So this would not be normalized by the number atoms (which is why I am saying it is extensive).

7

By the way, do you have any suggestions on how I can calculate the number of particles in these bins where their sum gives the total number of particles? N1+N2+N3+ … = N? Where Ni is the number of particles in the bin ith?

8

I think that what Axel meant (not sure) is that these formulas for calculating Cv derived within statistical mechanics’ formalism should hold for very big systems. If you have a very small system (such as the cylinder bins), maybe this formula no longer holds and you are calculating something without physical meaning (so it would be complicated to compare with the reference value of the Cv for the system at a given thermodynamic state). But maybe I misunderstood what you meant ^^

9

Yes, I agree with Axel’s comment if the expression for Cv on small scales is valid or not. However, people have tried to model the same small-scale systems using the standard models in thermodynamics in which we may raise the same concern on the validity of the analytical equations and proposing propagation under adiabatic expansion.

10

I suspect that what you are trying to do is beyond the capabilities of what can be easily done with variables and scripting in LAMMPS. The information of atoms/chunk is available in compute chunk/atom, but the more practical approach would be to implement your computation as a compute cv/chunk style or something along those lines. The chunk mechanism is quite powerful and flexible, but the more complex a property is to compute, the more likely is the need for a piece of custom C++ code.

11

Is this something you can contribute to our projects and be part of the authors if you have time? Our focus is mainly on modeling (FLASH) radiotherapy for treating cancer and we have facilities to check experimentally what we predict.

12

Sorry, but no. I am significantly underfunded for working on LAMMPS. Technically, only 10% of my time is currently covered. I have been fortunate that COVID has generated a lot of additional available time due to reducing other demands on my time, but that has now mostly gone away. I also have a long list of projects on my mind that I am more comfortable doing apart from all the normal effort of responding to questions here and managing pull requests on GitHub presence and so on that is consuming most of my time.

I have been waiting for some time now for a miracle to happen where someone says: I like you to keep working on LAMMPS full time and will pay you for it so you can quit your current job.
But the nature of miracles is that they very rarely happen.

1 Like
13

I think I had an idea on how you can achieve this calculation with the currently implemented chunk/ave!!!

I am not sure about the syntax of the “variable” command because I dont have it fresh in my mind right now, but if you have only water, maybe you can do the following:

variable my_mass atom mass[i]

and use this as input in the ave/chunk command together with the keyword none. Since the ave/chunk command will only consider the atoms that are in each respective bin for the calculation, it should output you the average mass found at each cylinder bin (not normalized by N, since you are using the keyword ‘none’). And then, if your system has only water, you can figure out the corresponding value of N in each bin by answering the questions: if I have this many grams (or whatever unit) of water in this bin, how many water molecules do I have? I know that it is in the realm of possibilities that part of the molecule is in a bin and another part in another, but maybe this gives you enough of an esitmate.

So, afterwards, given that the formula you need to use involve <Ek> (average kinetic energy, extensive), <Ek^2> (average kinetic energy (also extensive) square) for each bin, you can define also

compute ke_peratom all ke/atom
variable meanke atom c_ke_peratom
variable meanke2 atom v_meanke^2

And then, similarly to what you have written, you output all the things you need in something like:

fix heat_capacity_profile all ave/chunk 10 5 100 dvd v_mymass v_meanke v_meanke2 file heat_capacity.profile norm none

You then do the square of the output column corresponding to v_meanke manually for each of the bins. This should give you all the info you need to compute the Cv for each bin (together with boltzmann constant and etc).

14

Thanks a lot for the script and the time you have spent on this. I ran your script but it seems the numbers are growing with bins. Especially the mass. Maybe I need to share the entire script and we work on it in more detail. We may go offline and if we turn the results into a scientifically satisfactory work, we will turn it into an extra publication. Does it work for you? If so, please send me an email so we carry on.

15

On top of everything you’ve discussed: are you using rigid water models? If you are, your “heat capacity” will be even further off. Firstly, a rigid water molecule only has 6 spatial degrees of freedom, where a fully-flexible 3-point water molecule would have 9. Secondly, those degrees of freedom are not equally partitioned between atoms because hydrogen and oxygen atoms have different masses.

I have a colleague whom you might be interested to talk to. You can PM me for my email.

16

No, it is ReaxFF

17

My email: [email protected]

18

The preprint of this work is now available: https://arxiv.org/pdf/2403.05880.pdf
Please give me your critical feedback. Thank you all for your help and responses.

1 Like