Calculation of the static dielectric constant of a 3D periodic water system
Nathalie Smith
Dear CP2K users,
I am having trouble with a problem that I have seen sporadically in this forum, but didn't really find an answer to.
I am trying to calculate the static dielectric constant of water using the variance of the total dipole moment of the system, like shown in several papers and CP2K tutorials (e.g. https://www.cp2k.org/exercises:2014_ethz_mmm:monte_carlo_ice or https://www.cp2k.org/exercises:2018_uzh_acpc2:mol_sol). Currently, I am tyring to reproduce the MD simulation shown in the last link, just for longer simulation times (in the nanosecond scale) and in the NVT instead of the NVE ensemble. Therefore, I use the TIP3/Fw water model, PBCs in all directions, the spme method to deal with electrostatic interactions, and the berry phase formula for the calculation of the dipole moments. You can see all the parameters in the attached input file.
When I got the file with the dipole moments, I used the linux terminal command „cat filename.dat | grep Debye > newfilename.dat“ to extract the dipole moments in Debye, removed the „MM DIPOLE [BERRY PHASE](Debye)|“ section and used a python skript to calculate the static dielectric constant from it (also attached to this message). The problem: The value I get is way too low, around 15 and not around 80 as I would expect. In other simulations that I set up (using SPC/E water and a larger cell), the problem was even worse and I got a value around 8.
I checked my units and formulas several times and I honestly don’t think they are the reason I get wrong values. From what I have seen in the CP2K forum, this is a rather common problem; e.g. here https://lists.cp2k.org/archives/cp2k-user/2020-March/013043.html someone only got the right result after calculating the total dipole moment by themselves (which is something I could do for most MM methods, but not quite as easily for QM methods I guess, which I want to use later on. I know I could calculate molecular dipole moments using maximally localized Wannier centers, but this seems to be computationally more expensive).
So my question is: Is there something wrong with my calculations that I don’t see? Or is it a problem with the CP2K dipole moments (e.g. wrong unit, or wrong values altogether)? I am convinced that CP2K can be used for the calculation of static dielectric constants using the variance of the total dipole moment because it was already used for this purpose in the literature (e.g. https://pubs.acs.org/doi/10.1021/jp4103355; please note that I am aware of the differences between their and my calculations (MC instead of MD, ice instead of water, DFT instead of a MM water model etc.)).
Thanks in advance for anyone who reads my question, and to anyone who also gives an answer.
Best wishes,
Nathalie Smith