<div dir="ltr">Dear all,<div><br></div><div><div>I am using this link to learn how to calculate the static dielectric constant/relative permittivity from CP2K "moments" output files:</div><div><span style="color: rgb(0, 0, 255);"><br></span></div><div><span style="color: rgb(0, 0, 255);">https://www.cp2k.org/exercises:2014_ethz_mmm:monte_carlo_ice#task_1calculate_the_dielectric_constant</span><br></div><div>("Task 1: Calculate the dielectric constant")</div><div><br></div><div>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).</div><div><br></div><div>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:</div></div><div><br></div><div>- CP2K SPC/E simulations of water at 20 C give me: "moments.traj" files. Here is a snippet:</div><div><br></div><div><div class="prettyprint" style="background-color: rgb(250, 250, 250); border-color: rgb(187, 187, 187); border-style: solid; border-width: 1px; overflow-wrap: break-word;"><code class="prettyprint"><div class="subprettyprint"><span style="color: #000;" class="styled-by-prettify"> MM DIPOLE </span><span style="color: #660;" class="styled-by-prettify">[</span><span style="color: #000;" class="styled-by-prettify">BERRY PHASE</span><span style="color: #660;" class="styled-by-prettify">](</span><span style="color: #000;" class="styled-by-prettify">A</span><span style="color: #660;" class="styled-by-prettify">.</span><span style="color: #000;" class="styled-by-prettify">U</span><span style="color: #660;" class="styled-by-prettify">.)|</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #660;" class="styled-by-prettify">-</span><span style="color: #066;" class="styled-by-prettify">3.333675</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">0.611131</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">1.550260</span><span style="color: #000;" class="styled-by-prettify"><br> MM DIPOLE </span><span style="color: #660;" class="styled-by-prettify">[</span><span style="color: #000;" class="styled-by-prettify">BERRY PHASE</span><span style="color: #660;" class="styled-by-prettify">](</span><span style="color: #606;" class="styled-by-prettify">Debye</span><span style="color: #660;" class="styled-by-prettify">)|</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #660;" class="styled-by-prettify">-</span><span style="color: #066;" class="styled-by-prettify">8.473356</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">1.553341</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">3.940367</span><span style="color: #000;" class="styled-by-prettify"><br> MM DIPOLE </span><span style="color: #660;" class="styled-by-prettify">[</span><span style="color: #000;" class="styled-by-prettify">BERRY PHASE</span><span style="color: #660;" class="styled-by-prettify">]</span><span style="color: #000;" class="styled-by-prettify"> DERIVATIVE</span><span style="color: #660;" class="styled-by-prettify">(</span><span style="color: #000;" class="styled-by-prettify">A</span><span style="color: #660;" class="styled-by-prettify">.</span><span style="color: #000;" class="styled-by-prettify">U</span><span style="color: #660;" class="styled-by-prettify">.)|</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">0.000767</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #660;" class="styled-by-prettify">-</span><span style="color: #066;" class="styled-by-prettify">0.000211</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">0.001268</span><span style="color: #000;" class="styled-by-prettify"><br> MM DIPOLE </span><span style="color: #660;" class="styled-by-prettify">[</span><span style="color: #000;" class="styled-by-prettify">BERRY PHASE</span><span style="color: #660;" class="styled-by-prettify">](</span><span style="color: #000;" class="styled-by-prettify">A</span><span style="color: #660;" class="styled-by-prettify">.</span><span style="color: #000;" class="styled-by-prettify">U</span><span style="color: #660;" class="styled-by-prettify">.)|</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #660;" class="styled-by-prettify">-</span><span style="color: #066;" class="styled-by-prettify">3.600074</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">0.921496</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">1.793125</span><span style="color: #000;" class="styled-by-prettify"><br> MM DIPOLE </span><span style="color: #660;" class="styled-by-prettify">[</span><span style="color: #000;" class="styled-by-prettify">BERRY PHASE</span><span style="color: #660;" class="styled-by-prettify">](</span><span style="color: #606;" class="styled-by-prettify">Debye</span><span style="color: #660;" class="styled-by-prettify">)|</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #660;" class="styled-by-prettify">-</span><span style="color: #066;" class="styled-by-prettify">9.150475</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">2.342208</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">4.557669</span><span style="color: #000;" class="styled-by-prettify"><br> MM DIPOLE </span><span style="color: #660;" class="styled-by-prettify">[</span><span style="color: #000;" class="styled-by-prettify">BERRY PHASE</span><span style="color: #660;" class="styled-by-prettify">]</span><span style="color: #000;" class="styled-by-prettify"> DERIVATIVE</span><span style="color: #660;" class="styled-by-prettify">(</span><span style="color: #000;" class="styled-by-prettify">A</span><span style="color: #660;" class="styled-by-prettify">.</span><span style="color: #000;" class="styled-by-prettify">U</span><span style="color: #660;" class="styled-by-prettify">.)|</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">0.000754</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #660;" class="styled-by-prettify">-</span><span style="color: #066;" class="styled-by-prettify">0.000212</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #066;" class="styled-by-prettify">0.001278</span></div></code></div><div><br></div></div><div>- I keep the lines in [Debye] units and dump everything else. </div>
<style type="text/css">
p.p1 {margin: 0.0px 0.0px 0.0px 0.0px; font: 12.0px Menlo; color: #272ad8; background-color: #ffffff}
</style>
<div>- I then calculate a "total moment" for each line by summing the squares of each of the three components and taking the square root:</div><div><br></div><div><pre class="code python" style="padding: 0.7em 1em; font-family: Consolas, "Andale Mono WT", "Andale Mono", "Bitstream Vera Sans Mono", "Nimbus Mono L", Monaco, "Courier New", monospace; font-size: 14px; direction: ltr; background-color: rgb(251, 250, 249); color: rgb(51, 51, 51); box-shadow: rgb(204, 204, 204) -4px -4px 0.5em -0.3em inset; border-radius: 2px; overflow: auto; overflow-wrap: normal; border-width: 1px; border-style: solid; border-color: rgb(204, 204, 204);">M <span class="sy0" style="color: rgb(102, 204, 102);">=</span> np.<span class="me1" style="color: rgb(0, 102, 0);">sqrt</span><span class="br0" style="color: rgb(102, 204, 102);">(</span>np.<span class="kw2" style="color: rgb(0, 0, 0); font-weight: bold;">sum</span><span class="br0" style="color: rgb(102, 204, 102);">(</span>np.<span class="me1" style="color: rgb(0, 102, 0);">square</span><span class="br0" style="color: rgb(102, 204, 102);">(</span>M_vec<span class="br0" style="color: rgb(102, 204, 102);">)</span><span class="sy0" style="color: rgb(102, 204, 102);">,</span> axis<span class="sy0" style="color: rgb(102, 204, 102);">=</span><span class="nu0" style="color: rgb(204, 102, 204);">1</span><span class="br0" style="color: rgb(102, 204, 102);">)</span><span class="br0" style="color: rgb(102, 204, 102);">)</span> <span class="co1" style="color: rgb(128, 128, 128); font-style: italic;"># Take norm of dipole moment (from link above)</span>
</pre></div><div><span class="co1" style="color: rgb(128, 128, 128); font-style: italic;"><br></span></div><div><font color="#000000"><i>- </i></font>I then convert the total moment in Debye units to SI units [C*m] by multiplying the total moment in Debye by <span style="color: rgb(39, 42, 216); font-family: Menlo; font-size: 12px;">3.335641e-30</span>.</div><div>- I then execute the following line from the link: </div><div><br></div><div><pre class="code python" style="padding: 0.7em 1em; font-family: Consolas, "Andale Mono WT", "Andale Mono", "Bitstream Vera Sans Mono", "Nimbus Mono L", Monaco, "Courier New", monospace; font-size: 14px; direction: ltr; background-color: rgb(251, 250, 249); color: rgb(51, 51, 51); box-shadow: rgb(204, 204, 204) -4px -4px 0.5em -0.3em inset; border-radius: 2px; overflow: auto; overflow-wrap: normal; border-width: 1px; border-style: solid; border-color: rgb(204, 204, 204);">s <span class="sy0" style="color: rgb(102, 204, 102);">=</span> <span class="br0" style="color: rgb(102, 204, 102);">(</span>e*e*<span class="nu0" style="color: rgb(204, 102, 204);">4</span>*np.<span class="me1" style="color: rgb(0, 102, 0);">pi</span><span class="br0" style="color: rgb(102, 204, 102);">)</span>/<span class="br0" style="color: rgb(102, 204, 102);">(</span><span class="nu0" style="color: rgb(204, 102, 204);">3</span>*V*kb*T*angstrom2meter<span class="br0" style="color: rgb(102, 204, 102);">)</span>
</pre></div><div><br></div><div>Note: </div><div>* 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." </div><div>* 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.<br></div><div>* 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:</div><div><br></div><div><div class="prettyprint" style="background-color: rgb(250, 250, 250); border-color: rgb(187, 187, 187); border-style: solid; border-width: 1px; overflow-wrap: break-word;"><code class="prettyprint"><div class="subprettyprint"><span style="color: #008;" class="styled-by-prettify">var</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #660;" class="styled-by-prettify">=</span><span style="color: #000;" class="styled-by-prettify"> np</span><span style="color: #660;" class="styled-by-prettify">.</span><span style="color: #000;" class="styled-by-prettify">average</span><span style="color: #660;" class="styled-by-prettify">(</span><span style="color: #000;" class="styled-by-prettify">np</span><span style="color: #660;" class="styled-by-prettify">.</span><span style="color: #000;" class="styled-by-prettify">square</span><span style="color: #660;" class="styled-by-prettify">(</span><span style="color: #000;" class="styled-by-prettify">totalMoment</span><span style="color: #660;" class="styled-by-prettify">))</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #660;" class="styled-by-prettify">-</span><span style="color: #000;" class="styled-by-prettify"> np</span><span style="color: #660;" class="styled-by-prettify">.</span><span style="color: #000;" class="styled-by-prettify">square</span><span style="color: #660;" class="styled-by-prettify">(</span><span style="color: #000;" class="styled-by-prettify">np</span><span style="color: #660;" class="styled-by-prettify">.</span><span style="color: #000;" class="styled-by-prettify">average</span><span style="color: #660;" class="styled-by-prettify">(</span><span style="color: #000;" class="styled-by-prettify">totalMoment</span><span style="color: #660;" class="styled-by-prettify">))</span><span style="color: #000;" class="styled-by-prettify"><br></span></div></code></div><br>I also tried Python's variance and it seemed to generate a similar result:</div><div><br></div><div><div class="prettyprint" style="background-color: rgb(250, 250, 250); border-color: rgb(187, 187, 187); border-style: solid; border-width: 1px; overflow-wrap: break-word;"><code class="prettyprint"><div class="subprettyprint"><span style="color: #008;" class="styled-by-prettify">var</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #660;" class="styled-by-prettify">=</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #660;" class="styled-by-prettify">[</span><span style="color: #000;" class="styled-by-prettify">np</span><span style="color: #660;" class="styled-by-prettify">.</span><span style="color: #008;" class="styled-by-prettify">var</span><span style="color: #660;" class="styled-by-prettify">(</span><span style="color: #000;" class="styled-by-prettify">totalMoment</span><span style="color: #660;" class="styled-by-prettify">[:</span><span style="color: #000;" class="styled-by-prettify">i</span><span style="color: #660;" class="styled-by-prettify">+</span><span style="color: #066;" class="styled-by-prettify">1</span><span style="color: #660;" class="styled-by-prettify">])</span><span style="color: #000;" class="styled-by-prettify"> </span><span style="color: #008;" class="styled-by-prettify">for</span><span style="color: #000;" class="styled-by-prettify"> i </span><span style="color: #008;" class="styled-by-prettify">in</span><span style="color: #000;" class="styled-by-prettify"> range</span><span style="color: #660;" class="styled-by-prettify">(</span><span style="color: #000;" class="styled-by-prettify">len</span><span style="color: #660;" class="styled-by-prettify">(</span><span style="color: #000;" class="styled-by-prettify">totalMoment</span><span style="color: #660;" class="styled-by-prettify">))]</span><span style="color: #000;" class="styled-by-prettify"> </span></div></code></div><br>Thank you.</div></div>