<html>
  <head>
    <meta content="text/html; charset=utf-8" http-equiv="Content-Type">
  </head>
  <body text="#000000" bgcolor="#FFFFFF">
    <div class="moz-cite-prefix">Hi,<br>
      <br>
      No, that code is correct.<br>
      <br>
      I ran 4.6.7 on a water slab in the middle of a box with vacuum on
      top and bottom and I don't observe any issue with OpenMP only
      parallelization.<br>
      <br>
      Cheers,<br>
      <br>
      Berk<br>
      <br>
      On 06/07/2016 05:03 PM, Sheng Bi wrote:<br>
    </div>
    <blockquote cite="mid:5756E235.01A905.28023@mail.hust.edu.cn"
      type="cite">
      <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
      <meta name="Generator" content="Microsoft Word 15 (filtered
        medium)">
      <style><!--
/* Font Definitions */
@font-face
        {font-family:"Cambria Math";
        panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
        {font-family:DengXian;
        panose-1:2 1 6 0 3 1 1 1 1 1;}
@font-face
        {font-family:Calibri;
        panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
        {font-family:"\@DengXian";
        panose-1:2 1 6 0 3 1 1 1 1 1;}
@font-face
        {font-family:"Times New Roman \,serif";}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
        {margin:0in;
        margin-bottom:.0001pt;
        font-size:11.0pt;
        font-family:"Calibri",sans-serif;}
h1
        {mso-style-priority:9;
        mso-style-link:"Heading 1 Char";
        margin-top:12.0pt;
        margin-right:0in;
        margin-bottom:0in;
        margin-left:0in;
        margin-bottom:.0001pt;
        page-break-after:avoid;
        font-size:16.0pt;
        font-family:"Calibri Light",sans-serif;
        color:#2E74B5;
        font-weight:normal;}
a:link, span.MsoHyperlink
        {mso-style-priority:99;
        color:blue;
        text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
        {mso-style-priority:99;
        color:#954F72;
        text-decoration:underline;}
span.Heading1Char
        {mso-style-name:"Heading 1 Char";
        mso-style-priority:9;
        mso-style-link:"Heading 1";
        font-family:"Calibri Light",sans-serif;
        color:#2E74B5;}
.MsoChpDefault
        {mso-style-type:export-only;
        font-family:"Calibri",sans-serif;}
@page WordSection1
        {size:8.5in 11.0in;
        margin:1.0in 1.25in 1.0in 1.25in;}
div.WordSection1
        {page:WordSection1;}
--></style>
      <div class="WordSection1">
        <p class="MsoNormal"><span style="font-size:12.0pt">Hi Berk,<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">           
            Thanks for your response!<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">           
            However, my question is still there. In fact, I run
            simulation in one core (with multi-threadings), so domain
            decomposition could not be the reason for this phenomenon.
            Meanwhile, we use a quite large box (due to 3DC, Z length
            usually need a huge vacuum), so I think the water is within
            the box.<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">           
            Actually, I have printed dipcorrA for both gmx331 and gmx467
            and find that there is no difference. So it wouldn’t be the
          </span><span style="font-size:12.0pt">PBC trouble.<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">           
            I think <b>you are right about ewald_LRcorrection().</b> I
            check this function further and think it works fine. <o:p></o:p></span></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">           
            However, I think there is <b>something wrong in
              forcerec_set_excl_load() function (in forcerec.c file).</b>
            We know in ewald_LRcorrection(), for each thread, we loop
            over (i = fr-&gt;excl_load[t]; i &lt; fr-&gt;excl_load[t+1])
             to do dipole correction. The excl_load array is calculated
            in forcerec_set_excl_load(), as shown in following:<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">/ ******
            line 3018 ********* /<o:p></o:p></span></p>
        <div>
          <div>
            <p class="MsoNormal"><span style="font-size:12.0pt">    for
                (t = 1; t &lt;= fr-&gt;nthreads; t++)<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">    {<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">       
                ntarget = (ntot*t)/fr-&gt;nthreads;<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">       
                <b>while (i &lt; top-&gt;excls.nr &amp;&amp; n &lt;
                  ntarget)<o:p></o:p></b></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">       
                {<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">           
                for (j = ind[i]; j &lt; ind[i+1]; j++)<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">           
                {<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">               
                if (a[j] &gt; i)<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">               
                {<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">                   
                n++;<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">               
                }<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">           
                }<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">           
                i++;<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">       
                }<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">       
                fr-&gt;excl_load[t] = i;<o:p></o:p></span></p>
            <p class="MsoNormal"><span style="font-size:12.0pt">    }<o:p></o:p></span></p>
          </div>
        </div>
        <p class="MsoNormal"><span style="font-size:12.0pt">           
                        </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">           
            Suppose we have a system consisted by <b>2 types of
              molecules</b>. The whole system contains 2 <b>atoms</b>,
            1 for molecule <b>A(positively charged)</b> and 1 for
            molecule <b>B(negatively charged)</b>. We perform
            simulation in one core with one thread. So, for code above,
            we can deduce parameters like following:<o:p></o:p></span></p>
        <p class="MsoNormal"><b><span style="font-size:12.0pt">           
            </span></b><b><span style="font-size:12.0pt">top-&gt;excls.nr
              = 2<o:p></o:p></span></b></p>
        <p class="MsoNormal"><b><span style="font-size:12.0pt">           
              fr-&gt;nthreads = 1 <o:p></o:p></span></b></p>
        <p class="MsoNormal"><b><span style="font-size:12.0pt">           
              ntarget = ntot =1<o:p></o:p></span></b></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">           
            Since it <b>use while() statement here</b>, we can deduce
            that:<o:p></o:p></span></p>
        <p class="MsoNormal"><b><span style="font-size:12.0pt">           
              fr-&gt;excl_load[0] =0<o:p></o:p></span></b></p>
        <p class="MsoNormal"><b><span style="font-size:12.0pt">
                         fr-&gt;excl_load[1] =1<o:p></o:p></span></b></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">           
            However, <b>fr-&gt;excl_load[1] should be 2</b>. </span><span
            style="font-size:12.0pt">This wrong value lead to <b>missing
              a dipole correction on the 2<sup>th</sup> atom</b> in </span><span
            style="font-size:12.0pt">ewald_LRcorrection()</span><span
            style="font-size:12.0pt">. Because in </span><span
            style="font-size:12.0pt">ewald_LRcorrection()</span><span
            style="font-size:12.0pt">, we loop from </span><span
            style="font-size:12.0pt">fr-&gt;excl_load[0]</span><span
            style="font-size:12.0pt">=0 to </span><span
            style="font-size:12.0pt">fr-&gt;excl_load[1] =1 </span><span
            style="font-size:12.0pt">using for statement.<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">           
            I think this is the direct reason for wrong additional force
            on water molecule in my system. This error might result big
            trouble when someone using umbrella sampling in a system
            with </span><span style="font-size:12.0pt">system dipole.<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">           
          </span><span style="font-size:12.0pt">Please correct me if I
            am wrong.<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">Best
            regards!<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="font-size:12.0pt">Sheng <o:p></o:p></span></p>
        <p class="MsoNormal">--<br>
          PhD student<br>
          State Key Laboratory of Coal Combustion,<br>
          School of Energy and Power Engineering,<br>
          Huazhong University of Science and Technology (HUST),<br>
          Room 302 Power Building, <br>
          1037 Luoyu Road, Wuhan, Hubei 430074 China<br>
          Email: <a class="moz-txt-link-abbreviated" href="mailto:chrishengbee@hust.edu.cn">chrishengbee@hust.edu.cn</a><br>
          Phone: (86)27-87548122 <br>
          <a class="moz-txt-link-freetext" href="http://itp.energy.hust.edu.cn">http://itp.energy.hust.edu.cn</a> </p>
        <p class="MsoNormal"><span
            style="font-size:12.0pt;font-family:&quot;Times New
            Roman&quot;,serif"><o:p> </o:p></span></p>
        <div
          style="mso-element:para-border-div;border:none;border-top:solid
          #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in">
          <p class="MsoNormal" style="border:none;padding:0in"><b>From:
            </b><a moz-do-not-send="true" href="mailto:hess@kth.se">Berk
              Hess</a><br>
            <b>Sent: </b>Tuesday, June 7, 2016 5:49 PM<br>
            <b>To: </b><a moz-do-not-send="true"
              href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a><br>
            <b>Subject: </b>Re: [gmx-developers] Question about
            ewald_LRcorrection() function</p>
        </div>
        <p class="MsoNormal"><span
            style="font-size:12.0pt;font-family:&quot;Times New
            Roman&quot;,serif"><o:p> </o:p></span></p>
        <p class="MsoNormal"><span
            style="font-size:12.0pt;font-family:&quot;Times New
            Roman&quot;,serif">Hi,<br>
            <br>
            I think you misinterpreted the code in ewald_LRcorrection.
            It does the same things in 3.3.1 and 4.6.7.<br>
            My guess is that the issue you are observing is due to
            periodic boundary conditions. With group cut-off scheme
            charge groups are whole (not split by PBC), so water
            molecules are usually whole. Then the 3DC correction always
            works correctly when you have charge groups with zero net
            charge. In the Verlet cut-off scheme with domain
            decomposition, or with charged charge groups, (partial)
            charges moving over PBC change the system dipole. To get
            your system to run correctly you need to ensure that no
            water molecule moves over the periodic boundary in along the
            z-dimension. I think grompp will print a warning for this
            issue.<br>
            <br>
            Cheers,<br>
            <br>
            Berk<br>
            <br>
            On 2016-06-07 10:12, Sheng Bi wrote:</span><span
            style="font-size:12.0pt;font-family:&quot;Times New
            Roman&quot;,serif"><o:p></o:p></span></p>
        <blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
          <p class="MsoNormal">Sorry, some type-offs in previous
            question.<o:p></o:p></p>
          <p class="MsoNormal">In gmx-4.6.7, ewald_LRcorrection()<b>.</b>
            We loop over excluded neighbours only in condition <b>calc_excl_corr=1</b>.
            But when using Verlet cut-off scheme, <b>calc_excl_corr
              will be 0</b>.<o:p></o:p></p>
          <p class="MsoNormal"> <o:p></o:p></p>
          <p class="MsoNormal"> </p>
          <p class="MsoNormal">--<br>
            PhD student<br>
            State Key Laboratory of Coal Combustion,<br>
            School of Energy and Power Engineering,<br>
            Huazhong University of Science and Technology (HUST),<br>
            Room 302 Power Building, <br>
            1037 Luoyu Road, Wuhan, Hubei 430074 China<br>
            Email: <a moz-do-not-send="true"
              href="mailto:chrishengbee@hust.edu.cn">chrishengbee@hust.edu.cn</a><br>
            Phone: (86)27-87548122 <br>
            <a moz-do-not-send="true"
              href="http://itp.energy.hust.edu.cn">http://itp.energy.hust.edu.cn</a>
            <o:p></o:p></p>
          <p class="MsoNormal"><span
              style="font-size:12.0pt;font-family:&quot;Times New Roman
              ,serif&quot;,serif"> </span></p>
          <div style="border:none;border-top:solid #E1E1E1
            1.0pt;padding:3.0pt 0in 0in 0in">
            <p class="MsoNormal"><b>From: </b><a moz-do-not-send="true"
                href="mailto:chrishengbee@hust.edu.cn">Sheng Bi</a><br>
              <b>Sent: </b>Tuesday, June 7, 2016 1:28 AM<br>
              <b>To: </b><a moz-do-not-send="true"
                href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a><br>
              <b>Subject: </b>[gmx-developers] Question about
              ewald_LRcorrection() function<o:p></o:p></p>
          </div>
          <p class="MsoNormal"><span
              style="font-size:12.0pt;font-family:&quot;Times New Roman
              ,serif&quot;,serif"> </span></p>
          <p class="MsoNormal">Dear GMX developers:</p>
          <p class="MsoNormal">               I find a abnormal
            phenomenon when using <b>gromacs-4.6.7</b>. I set up a
            system containing two walls(consist of freeze atoms), one
            wall has positively charges on atoms, another has negatively
            charges.  I put an water molecule(rigid spc/spce/tip3p
            model) between these two walls. <b>Whole system seems like
              below</b>:</p>
          <p class="MsoNormal">              
            ****************************************</p>
          <p class="MsoNormal">               * -         
            H                          
            +                                          *             </p>
          <p class="MsoNormal">               * -             \         
                           +            
            vacuum                              *</p>
          <p class="MsoNormal">              * -            
              O                    
            +                                          *</p>
          <p class="MsoNormal">               * -           
             /                       
            +                                          *</p>
          <p class="MsoNormal">               * -          H           
                           +                                          *</p>
          <p class="MsoNormal">              
            ****************************************</p>
          <p class="MsoNormal">               I use <b>Verlet cut-off
              scheme, PME method, and ewald-geometry=3DC</b>.
            Temperature of whole system is at <b>0K</b>, controlled by
            <b>v-rescale method</b>. I run whole simulation <b>only in
              OPENMP parallel</b>.</p>
          <p class="MsoNormal">               One can have following
            expected behavior of water molecule without MD simulation:
            Since it can be seen as an infinite panel capacitor, <b>there
              is constant electric field between two walls. So the water
              will change it orientation, but won’t change it’s
              position(Since 0K)</b>.</p>
          <p class="MsoNormal">               This phenomenon can be <b>confirmed
              by gromcas-3.3.1</b> (although there is position shift,
            but quite small). Also, if one use higher version of
            gromacs, such as 4.6.7 or higher <b>with ewald-geometry=3D,
              either no problem</b>. However, if one use <b>3DC</b>
            (and it should be), you can see <b>an obvious force on
              water molecule</b>.</p>
          <p class="MsoNormal">               I hacked into source code,
            and find this might due to <b>ewald_LRcorrection()</b>
            function:</p>
          <p class="MsoNormal">               Firstly, in gmx-4.6.7,
            ewald_LRcorrection()<b>.</b> We <b>loop over excluded
              neighbours only in condition calc_excl_corr=0. But when
              using verlet cut-off scheme, calc_excl_corr will be 1</b>.
            In gmx-3.3.1, we do loop over excluded neighbours without
            additional if() statement.</p>
          <p class="MsoNormal">               Secondly, in gmx-4.6.7, ,
            ewald_LRcorrection(). When do dipole correction on force, we
            loop over <b>excl_load</b> (Exclusion load distribution
            over the threads). In fact, most water model will have
            nrexel=2. So, as a result, we <b>only do dipole correction
              on O atom, but do nothing on H atoms</b> and we will see
            an obvious force on the whole water molecule. I think this
            is not purpose of 3DC<span style="font-family:DengXian"
              lang="ZH-CN">。</span></p>
          <p class="MsoNormal">               Best regards<span
              style="font-family:DengXian" lang="ZH-CN">!</span></p>
          <p class="MsoNormal">Sheng Bi</p>
          <p class="MsoNormal"><b>               </b></p>
          <p class="MsoNormal">---<br>
            PhD student<br>
            State Key Laboratory of Coal Combustion,<br>
            School of Energy and Power Engineering,<br>
            Huazhong University of Science and Technology (HUST),<br>
            Room 302 Power Building, <br>
            1037 Luoyu Road, Wuhan, Hubei 430074 China<br>
            Email: <a moz-do-not-send="true"
              href="mailto:chrishengbee@hust.edu.cn">chrishengbee@hust.edu.cn</a><br>
            Phone: (86)27-87548122 <br>
            <a moz-do-not-send="true"
              href="http://itp.energy.hust.edu.cn">http://itp.energy.hust.edu.cn</a>
            <o:p></o:p></p>
          <p class="MsoNormal"><span
              style="font-size:12.0pt;font-family:&quot;Times New Roman
              ,serif&quot;,serif"> </span><o:p></o:p></p>
          <p class="MsoNormal"> </p>
          <p class="MsoNormal"><span style="font-family:&quot;Times New
              Roman&quot;,serif"><br>
              <br>
              <o:p></o:p></span></p>
        </blockquote>
        <p class="MsoNormal"><span style="font-family:&quot;Times New
            Roman&quot;,serif"><o:p> </o:p></span></p>
        <p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
      </div>
    </blockquote>
    <br>
  </body>
</html>