[CP2K-user] Problem in calculating liquid static dielectric constant/relative permittivity

[Nt]

Dear all,

I am using this link to learn how to calculate the static dielectric 
constant/relative permittivity from CP2K "moments" output files:

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

I have started by trying to calculate the dielectric constant of water at 
20 C so that I know I am doing things correctly.  I'm using CP2K "moments" 
output files from 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/20191113/921079c9/attachment.htm>