[CP2K-user] Dielectric constant calculation in CP2K

[Vladimir Rybkin]

Dear Nt,

I'm not sure whether the script is up-to-date as it's oldish and has been 
used for excercises 5 years ago for a course which does not exist anymore 
and by a group which does not exist anymore. It is there for the historical 
reasons. A few hints:
1) you may want to write a short scipt yourself;
2) the model used is not very accurate so you shouldn't be surprised with 
the number of epsilon far from 80. In the excercise there is a link to the 
original DFT calculation:
https://pubs.acs.org/doi/10.1021/jp4103355

As you may see even DFT values vary considerably.

Yours, 

Vladimir

пятница, 15 ноября 2019 г., 9:01:10 UTC+1 пользователь Nt написал:
>
> Dear all:
>
> I hope it's OK to repost this question, which didn't get any replies when 
> first posted.  I'm hoping that someone who has experience with using CP2K 
> to calculate the dielectric constant will see it because I'm stuck.
>
> At a temperature of 20 C and atmospheric pressure, the dielectric constant 
> of water has been measured at approximately 80.  However, when I calculate 
> it based on a FIST SPC/E simulation and a CP2K tutorial (link provided 
> below), I get a result of <1.  I am wondering if maybe there is a problem 
> with my understanding of the units from CP2K or the way I outputted the 
> dipole moment.  This is pretty new to me, so I could very well have made a 
> simple mistake somewhere.
>
> 1. I added the following in:
>
> &Force_Eval/&MM/&PRINT
>
> &DIPOLE
>   PERIODIC TRUE
>   &EACH
>     MD 1
>   &END
>   FILENAME = moments.traj
>   ADD_LAST NUMERIC
> &END DIPOLE
>
> 2. The output file "moments.traj" contains lines like these at each time 
> step with the components of the dipole moment:
>  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
>
> I extracted the lines in [Debye] units and then turned to a tutorial at 
> cp2k.org (link below), where I followed the procedure on how to calculate 
> the dielectric constant.  
>
>
> https://www.cp2k.org/exercises:2014_ethz_mmm:monte_carlo_ice#task_1calculate_the_dielectric_constant
> ("Task 1: Calculate the dielectric constant") 
>
> - First I used the tutorial to calculate a total dipole moment from the 
> three columns in the output snippet above (sum of squares of each of the 
> three components, followed by taking the square root):
>
> M = np.sqrt(np.sum(np.square(M_vec), axis=1)) # Take norm of dipole moment (from link above)
>
>
> - Then I converted the total moment in Debye units to SI units [C*m] by 
> multiplying the total moment in Debye by 3.335641e-30.
>
> 3. 
> - I executed the following line from the link, where I had converted 
> everything to SI units: 
>
> s = (e*e*4*np.pi)/(3*V*kb*T*angstrom2meter*epsilon_0)
>
>
> When I multiplied the result with the total dipole moment I calculated 
> above, the answer was far from 80.  I used another version of the equation 
> that is more common in published literature without the epsilon_0 term in 
> the denominator and without the "e" terms.  It gives me an answer much less 
> than 1.
>
> Here is how I calculated the variance of the total dipole moment:
> 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))]
>
> If anyone has any ideas on why I am not having any success, I would be 
> very grateful because I am stuck.
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20191119/0686bc43/attachment.htm>