search

lab-cosmo / i-pi-dev_archive

Development version of i-PI
21 stars 12 forks source link

Fix bug in temperature output for selected atoms #191

Closed mahrossi closed 6 years ago

mahrossi commented 6 years ago

Let's just say that the implementation currently in i-PI was giving different temperatures for different degrees of freedom, when one was choosing to output things like "temperature(O)" or "temperature(1)". Check what I did. I think it could be more transparent, and one has to check whether the bug is not also present for kin_cv . The way "ncount" is being counted is not transparent. What I did fixes it for PI and classical nuclei MD simulations. Before merging I would be glad if someone checked the code further as I have no time. If testing it with a molecule, please remember to set <fixcom> False </fixcom> in the input file.

mahrossi commented 6 years ago

P.S.: Sorry I started from a branch that was a bit behind master for this fix. I only changed properties.py

venkatkapil24 commented 6 years ago

The way things are implemented, the temperature of an atom when the COM is Fixed is not returned correctly. This should be discussed.

mahrossi commented 6 years ago

Thanks Venkat -- regarding the fixed COM, did you try a periodic system? Or only molecules?

venkatkapil24 commented 6 years ago

I had tried molecules only. But I just ran a simulation of 32 molecules of water with FIxcom=True and the total temperature, temperature of O atoms and temperature of H atoms is coming out to be fine. But if you ask for the temperature of lets say species j, then you will get infinity since it calculates K.E of the atom / (ncount 3 beads - nconstr), which for a classical simulation is () / (3 -3).

venkatkapil24 commented 6 years ago

Have tested it for a C H debye harmonic potential, C C debye harmonic potential and 32 molecules of qtip4pf water, for P=1,8 and computed the total temperature, temperature of j^th and k^th species, temperature of H,O,C species and all seem correct. Unless someone wants to carry out further tests or add something more, I think this branch is good for a merge.

mahrossi commented 6 years ago

Good! Merge?

pjuda commented 6 years ago

Gas phase qtip4pf water molecule at 300 K 4 kind of simulations: Classical MD (CL), Classical MD with fixed center of mass (CL-FC), PIMD with 4 beads (QN), PIMD with 4 beads and fixed center of mass (QN-FC)

We use 100 ps (200 ps in the case of CL)of simulation time to calculate the temperature of individual atoms, species and the average.

Here are the results (inputs attached):

inputs.zip --------| | SIM | H1 | H2 | O | H | AVG | | CL-FC | 312.723 +/- 5.522 | 299.390 +/- 5.401 | 300.508 +/- 0.618 | 306.057 +/- 4.859 | 304.207 +/- 3.410 | | CL | 303.842 +/- 4.565 | 294.778 +/- 4.138 | 297.008 +/- 5.064 | 299.271 +/- 3.685 | 298.498 +/- 3.223 | | QN-FC | 302.072 +/- 1.772 | 298.761 +/- 2.041 | 298.507 +/- 1.420 | 300.417 +/- 1.519 | 299.780 +/- 1.173 | | QN | 301.928 +/- 1.643 | 300.921 +/- 1.865 | 299.161 +/- 2.453 | 301.425 +/- 1.374 | 300.670 +/- 1.410 |

pjuda commented 6 years ago

Fixatoms simulation after 300 ps: input.txt (used autocorr maxlag 2000 and drop 1000)
SIM H1 H2 O1 H AVG
Fixatoms 297.547 +/- 2.947 300.512 +/- 3.315 307.670 +/- 4.488 299.969 +/- 0.084 300.059 +/- 0.077
ceriottm commented 6 years ago

OK this looks good to me. I think also it is harmless enough that it can go in even if we haven't tested it in all possible scenarios. A warning should be added in the source mentioning that it is a bit silly to change beads.p back and forth, but I was checking with Venkat and replicating the functionality without duplicating a lot of code would be a pain.