[CP2K-user] Dielectric constant calculation in CP2K

[Nt]

Dear all:

I hope it's OK to repost this question, which didn't get any replies when I 
asked it recently, in case someone who is familiar with calculating the 
dielectric constant in CP2K sees it.  If the question is unclear, please 
feel free to let me know and I will fix it.

I am using the link below to learn how to calculate the static dielectric 
constant/relative permittivity from CP2K "moments" output files.  I am 
doing something wrong because I can't calculate the dielectric constant of 
water at 20 C, which has been measured at ~80.  I have wondered if there is 
a problem with the example at the link or if there's something I don't 
understand correctly about the units required for the equation or outputted 
by CP2K.

https://www.cp2k.org/exercises:2014_ethz_mmm:monte_carlo_ice#task_1calculate_the_dielectric_constant
("Task 1: Calculate the dielectric constant")

I'm using total dipole moments for the "M" parameter from CP2K "moments" 
output files outputted by classical (SPC/E) simulations.  I have written my 
own Python script to deal with the equation at the link that is based on 
the Python script at the link (I can't figure out how to make the script 
provided there to work with my CP2K dipole moment output files).

I have tried using two versions of the equation at the link--the one shown 
and the and one without the epsilon_0 term in the denominator (which is 
identical to a publication I'm interested in).  Neither gives me anything 
like the right answer.  I'm pretty new to this and wonder if maybe there's 
something I don't understand about the CP2K moments output or units.  Here 
is a brief summary of what I have done.  I would be extremely appreciative 
of any advice on possible mistakes:

- CP2K SPC/E simulations of water at 20 C give me: "moments.traj" files. 
Here is a snippet:

     MM DIPOLE [BERRY PHASE](A.U.)|                    -3.333675    0.611131
   1.550260
     MM DIPOLE [BERRY PHASE](Debye)|                   -8.473356    1.553341
   3.940367
     MM DIPOLE [BERRY PHASE] DERIVATIVE(A.U.)|          0.000767   -0.000211
   0.001268
     MM DIPOLE [BERRY PHASE](A.U.)|                    -3.600074    0.921496
   1.793125
     MM DIPOLE [BERRY PHASE](Debye)|                   -9.150475    2.342208
   4.557669
     MM DIPOLE [BERRY PHASE] DERIVATIVE(A.U.)|          0.000754   -0.000212
   0.001278

- I keep the lines in [Debye] units and dump everything else.  
- I then calculate a "total moment" for each line by summing the squares of 
each of the three components and taking the square root:

M = np.sqrt(np.sum(np.square(M_vec), axis=1)) # Take norm of dipole moment (from link above)


*- *I then convert the total moment in Debye units to SI units [C*m] by 
multiplying the total moment in Debye by 3.335641e-30.
- I then execute the following line from the link: 

s = (e*e*4*np.pi)/(3*V*kb*T*angstrom2meter)


Note: 
* I don't know where the "e" terms come from--they are not part of the Task 
1 equation and I have also not seen them anywhere else in publications.  I 
tried both excluding them and including them and excluding and including 
the "epsilon_0" term in the denominator; nothing works to give me a result 
like "80."  
* Also, I think the "angstrom2meter" term should be "1e-30" instead of 
"1e-10" as shown at the link because it is meant to operate on the volume 
term.  I did try using 1e-10 in case I misunderstood something and the 
result is still wrong.
* Finally, here is how I calculated variance.  I don't understand what 
"weights" in the code at the link is, so left that term out:

var = np.average(np.square(totalMoment)) - np.square(np.average(totalMoment
))

I also tried Python's variance and it seemed to generate a similar result:

var  = [np.var(totalMoment[:i+1]) for i in range(len(totalMoment))]

Thank you.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20191115/0600be67/attachment.htm>