Bug? The conserved quantity is NOT conserved after flexible cell MD is restarted

[Rustam]

Hi guys,

As Amandine Chouanard mentioned in a recent post the conserved quantity shows discontinuities when a flexible cell MD simulation is restarted. I used her input file and found that restart works fine with NVT and NPT_I simulations, but NOT with NPT_F (see attached figure). The discontinuities in NVT_F the conserved quantity are larger if restarts are less frequent (i.e. trajectory pieces are longer).

It seems that the conserved quantities in NPT_I and NPT_F simulations are computed by the same subroutine in motion/md_conserved_quantities.F:390. So it seems that the restart file contains ALL data on the current state of the barostat and thermostat to ensure a smooth restart.

Is it possible that this has something to do with the CELL_REF setting? Something along the lines that the grids before and after the restart are different? I tried printing the GRID info but it seems to keep the number of points constant...

I will appreciate any suggestions/comments,

Rustam

[Matthias Krack]

Hi Rustam

in the first place I would suggest to use a compatible pair of basis set and potential for Mg. The Mg TZV2P-GTH basis set was generated for a Mg GTH-*-q10, but a GTH-BLYP-q2 is employed. Unfortunately, the GTH_BASIS_SET file provides only generic labels. I will fix that. Please, run the case with a valid combinations of potential and basis set and report if the error persists. Note, that the Mg q10 potential requires a very large cutoff and for the Mg q2 potential there is only a DZVP basis set available in the file BASIS_SET.

Best regards

Matthias

PS: The basis sets in GTH_BASIS_SETS are usually not well suited for (dense) solid state systems as they have been generated mainly for molecular systems.

[Rustam]

Hi Matthias,

Thanks for the great advice! I tried changing the basis set to the one compatible with Mg q2 and low cutoff. Unfortunately, the problem remains as shown in the new figure.

On a different note, what basis sets would you recommend for solid state systems? All-electron basis sets? Or do you think I should generate a pair of potential/basis for the problem at hand?

Thanks again,

Rustam

PS. The input file is also attached.

[Matthias Krack]

Dear Rustam

I would recommend to try a Mg-q10 MOLOPT basis set.

Best regards

Matthias

[Matt W]

Hi,

(I could well be wrong here, this is sort of anecdotal rather than rigorously checking the details) I think that the conserved quantity for the CSVR thermostat is quite a complicated thing compared to other thermostats. To my memory, the proper conserved quantity requires sums over information from all previous states in a non-trivial way that is tough to implement efficiently (i.e. adding new information at every step without complete recomputation).

After all that waffle, I suspect the 'problem' is the CSVR thermostat, and it is not actually an indication of bad behaviour, but a question of efficiency in calculating the conserved quantity (i.e. dynamics are correct, but it is hard to be sure without a conserved quantity). Maybe try Nose-Hoover after equilibriating with CSVR to be sure?

Matt

[Rustam]

Hi Matt,

I hope you are right and the dynamics is correct even though there are discontinuities in the conserved quantity after restarts. Unfortunately, Nose thermostat has the same problem (see the new figure) so there is no way so far to tell if the dynamics is correct.

Rustam