Dielectric constant calculation in CP2K

[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.

[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 написал:

[Ant Stone]

Dear Vladimir:

Thank you very much indeed for your reply.  I haven't managed to find anyone who has performed dielectric constant calculations in CP2K so far and my results have been so strange (always ~1 for water no matter which T I use in my simulations) that I have been wondering if maybe CP2K is not appropriate for this purpose for some reason or whether there's a problem with the way I outputted my dipole moments files.  I just can't figure out what is going wrong. 

Thank you for your suggestion; I did write a script myself, but keep getting values of ~1 (instead of ~80) at 20 C.  I have tried several versions of the fluctuation expression I provided in my initial post, I have checked my units repeatedly, and nothing helps.

For the past few days I have been trying to understand the tutorial script you referred to in case there's a hint there as to why my own script doesn't generate reasonable results.  However, there are a few points in the tutorial script that I don't understand, such as what the "weights" parameter is and where it comes from (it looks like the dipole moment files generated by CP2K and used in the tutorial script must have been formatted differently to current-day moment files so that makes it tough to figure out what is going on--i.e., exactly which data the tutorial script is using).  It's also not possible to tell exactly what M_vec is made up of.  Thank you very much for the link to the paper, but unfortunately, it doesn't provide the kind of nitty gritty I need to follow the example in the tutorial and my own calculations.

I have been hoping that someone here who has run these calculations successfully might see my post and give me some hints or maybe even let me know if CP2K is not the method to use for this.  Thanks a lot anyway.

[Vladimir Rybkin]

Dear Nt, 

you write:

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.

So, if I understand correctly, you used to formulas: one from the script and another one, that you find more proper. The latter gives less than 1. What is the number from the first one? 

GENERALLY:

1. Run the trajectory;

2. Extract dipoles step by step (not ONE dipole moment as you say, but all of them, one at a step);

3. Compute dipole moment autocorrelation function;

4. Apply the Kubo formula.

It is clear from your original message, that you do not do it. Rather, you tell us about ONE dipole moment you have computed.

Yours,

Vladimir

вторник, 19 ноября 2019 г., 16:51:52 UTC+1 пользователь Nt написал:

[Ant Stone]

Dear Vladimir:

Thank you for asking and for your interest.  Both formulae, which are similar to each other (based on the "fluctuation" method), give me answers of <1.  Both formulae are of the form, static dielectric constant = 1 + x and my problem is that the x term is far too small.

I have searched through the literature quite extensively and don't see any alternatives to the formula I provided here or the one that the CP2K tutorial uses.  I just can't figure out why my result is so far off what it should be.  My water density is consistent with experiments and I have checked and rechecked my units.  I'm wondering if there is something fundamental I don't understand about the CP2K dipole moments output.

On Wed, 20 Nov 2019 at 09:46, Vladimir Rybkin <rybk...@gmail.com> wrote:

Dear Nt, 

you write:

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.

So, if I understand correctly, you used to formulas: one from the script and another one, that you find more proper. The latter gives less than 1. What is the number from the first one? 

Yours,

Vladimir

вторник, 19 ноября 2019 г., 16:51:52 UTC+1 пользователь Nt написал:

Dear Vladimir:

Thank you very much indeed for your reply.  I haven't managed to find anyone who has performed dielectric constant calculations in CP2K so far and my results have been so strange (always ~1 for water no matter which T I use in my simulations) that I have been wondering if maybe CP2K is not appropriate for this purpose for some reason or whether there's a problem with the way I outputted my dipole moments files.  I just can't figure out what is going wrong. 

Thank you for your suggestion; I did write a script myself, but keep getting values of ~1 (instead of ~80) at 20 C.  I have tried several versions of the fluctuation expression I provided in my initial post, I have checked my units repeatedly, and nothing helps.

For the past few days I have been trying to understand the tutorial script you referred to in case there's a hint there as to why my own script doesn't generate reasonable results.  However, there are a few points in the tutorial script that I don't understand, such as what the "weights" parameter is and where it comes from (it looks like the dipole moment files generated by CP2K and used in the tutorial script must have been formatted differently to current-day moment files so that makes it tough to figure out what is going on--i.e., exactly which data the tutorial script is using).  It's also not possible to tell exactly what M_vec is made up of.  Thank you very much for the link to the paper, but unfortunately, it doesn't provide the kind of nitty gritty I need to follow the example in the tutorial and my own calculations.

I have been hoping that someone here who has run these calculations successfully might see my post and give me some hints or maybe even let me know if CP2K is not the method to use for this.  Thanks a lot anyway.

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/abab9c1d-6d8a-45b0-9165-6b2c30be811b%40googlegroups.com.

[Vladimir Rybkin]

Dear Nt,

1) it would be helpful to see your script;

2) please double check the units;

3) have you tried the very script from the excercise (not only some formulae from it)?

Yours,

Vladimir

среда, 20 ноября 2019 г., 19:07:19 UTC+1 пользователь Nt написал:

Dear Vladimir:

Thank you for asking and for your interest.  Both formulae, which are similar to each other (based on the "fluctuation" method), give me answers of <1.  Both formulae are of the form, static dielectric constant = 1 + x and my problem is that the x term is far too small.

I have searched through the literature quite extensively and don't see any alternatives to the formula I provided here or the one that the CP2K tutorial uses.  I just can't figure out why my result is so far off what it should be.  My water density is consistent with experiments and I have checked and rechecked my units.  I'm wondering if there is something fundamental I don't understand about the CP2K dipole moments output.

On Wed, 20 Nov 2019 at 09:46, Vladimir Rybkin <rybk...@gmail.com> wrote:

Dear Nt, 

you write:

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.

So, if I understand correctly, you used to formulas: one from the script and another one, that you find more proper. The latter gives less than 1. What is the number from the first one? 

Yours,

Vladimir

вторник, 19 ноября 2019 г., 16:51:52 UTC+1 пользователь Nt написал:

Dear Vladimir:

Thank you very much indeed for your reply.  I haven't managed to find anyone who has performed dielectric constant calculations in CP2K so far and my results have been so strange (always ~1 for water no matter which T I use in my simulations) that I have been wondering if maybe CP2K is not appropriate for this purpose for some reason or whether there's a problem with the way I outputted my dipole moments files.  I just can't figure out what is going wrong. 

Thank you for your suggestion; I did write a script myself, but keep getting values of ~1 (instead of ~80) at 20 C.  I have tried several versions of the fluctuation expression I provided in my initial post, I have checked my units repeatedly, and nothing helps.

For the past few days I have been trying to understand the tutorial script you referred to in case there's a hint there as to why my own script doesn't generate reasonable results.  However, there are a few points in the tutorial script that I don't understand, such as what the "weights" parameter is and where it comes from (it looks like the dipole moment files generated by CP2K and used in the tutorial script must have been formatted differently to current-day moment files so that makes it tough to figure out what is going on--i.e., exactly which data the tutorial script is using).  It's also not possible to tell exactly what M_vec is made up of.  Thank you very much for the link to the paper, but unfortunately, it doesn't provide the kind of nitty gritty I need to follow the example in the tutorial and my own calculations.

I have been hoping that someone here who has run these calculations successfully might see my post and give me some hints or maybe even let me know if CP2K is not the method to use for this.  Thanks a lot anyway.

--
You received this message because you are subscribed to the Google Groups "cp2k" group.

To unsubscribe from this group and stop receiving emails from it, send an email to cp...@googlegroups.com.