[CP2K-user] Dielectric constant calculation in CP2K

[kris Gu]

Dear Vladimir,
I just came across the same problem as you.
First, when you calculate the variance of the totoal moment, you need to 
use the norm of dipole vector, rather than directly calculate the average 
number like this:  np.average(totalMoment).

Second, I recalculate the dipole by myself using sum(coordination X point 
charge), it gives me a  different dipole number compared to the CP2K output 
dipole moment, which finally gives me a correct dielectric constant, I 
don't know why there's such a difference as I can not find the way of CP2K 
calculating dipole moment.

I hope it also works in your case.

Best,

Chris

在 2019年11月15日星期五 UTC+8下午4:01:10，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/20200314/43b7d7eb/attachment.htm>