[gmx-developers] Question about ewald_LRcorrection() function
Berk Hess
hess at kth.se
Wed Jun 8 11:21:57 CEST 2016
Hi,
No, that code is correct.
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.
Cheers,
Berk
On 06/07/2016 05:03 PM, Sheng Bi wrote:
>
> Hi Berk,
>
> Thanks for your response!
>
> 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.
>
> Actually, I have printed dipcorrA for both gmx331 and gmx467 and find
> that there is no difference. So it wouldn’t be the PBC trouble.
>
> I think *you are right about ewald_LRcorrection().* I check this
> function further and think it works fine.
>
> However, I think there is *something wrong in forcerec_set_excl_load()
> function (in forcerec.c file).* We know in ewald_LRcorrection(), for
> each thread, we loop over (i = fr->excl_load[t]; i <
> fr->excl_load[t+1]) to do dipole correction. The excl_load array is
> calculated in forcerec_set_excl_load(), as shown in following:
>
> / ****** line 3018 ********* /
>
> for (t = 1; t <= fr->nthreads; t++)
>
> {
>
> ntarget = (ntot*t)/fr->nthreads;
>
> *while (i < top->excls.nr && n < ntarget)*
>
> {
>
> for (j = ind[i]; j < ind[i+1]; j++)
>
> {
>
> if (a[j] > i)
>
> {
>
> n++;
>
> }
>
> }
>
> i++;
>
> }
>
> fr->excl_load[t] = i;
>
> }
>
> Suppose we have a system consisted by *2 types of molecules*. The
> whole system contains 2 *atoms*, 1 for molecule *A(positively
> charged)* and 1 for molecule *B(negatively charged)*. We perform
> simulation in one core with one thread. So, for code above, we can
> deduce parameters like following:
>
> ***top->excls.nr = 2*
>
> *fr->nthreads = 1 *
>
> *ntarget = ntot =1*
>
> Since it *use while() statement here*, we can deduce that:
>
> *fr->excl_load[0] =0*
>
> * fr->excl_load[1] =1*
>
> However, *fr->excl_load[1] should be 2*. This wrong value lead to
> *missing a dipole correction on the 2^th atom* in
> ewald_LRcorrection(). Because in ewald_LRcorrection(), we loop from
> fr->excl_load[0]=0 to fr->excl_load[1] =1 using for statement.
>
> 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 system dipole.
>
> Please correct me if I am wrong.
>
> Best regards!
>
> Sheng
>
> --
> PhD student
> State Key Laboratory of Coal Combustion,
> School of Energy and Power Engineering,
> Huazhong University of Science and Technology (HUST),
> Room 302 Power Building,
> 1037 Luoyu Road, Wuhan, Hubei 430074 China
> Email: chrishengbee at hust.edu.cn
> Phone: (86)27-87548122
> http://itp.energy.hust.edu.cn
>
> *From: *Berk Hess <mailto:hess at kth.se>
> *Sent: *Tuesday, June 7, 2016 5:49 PM
> *To: *gmx-developers at gromacs.org <mailto:gmx-developers at gromacs.org>
> *Subject: *Re: [gmx-developers] Question about ewald_LRcorrection()
> function
>
> Hi,
>
> I think you misinterpreted the code in ewald_LRcorrection. It does the
> same things in 3.3.1 and 4.6.7.
> 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.
>
> Cheers,
>
> Berk
>
> On 2016-06-07 10:12, Sheng Bi wrote:
>
> Sorry, some type-offs in previous question.
>
> In gmx-4.6.7, ewald_LRcorrection()*.* We loop over excluded
> neighbours only in condition *calc_excl_corr=1*. But when using
> Verlet cut-off scheme, *calc_excl_corr will be 0*.
>
> --
> PhD student
> State Key Laboratory of Coal Combustion,
> School of Energy and Power Engineering,
> Huazhong University of Science and Technology (HUST),
> Room 302 Power Building,
> 1037 Luoyu Road, Wuhan, Hubei 430074 China
> Email: chrishengbee at hust.edu.cn <mailto:chrishengbee at hust.edu.cn>
> Phone: (86)27-87548122
> http://itp.energy.hust.edu.cn
>
> *From: *Sheng Bi <mailto:chrishengbee at hust.edu.cn>
> *Sent: *Tuesday, June 7, 2016 1:28 AM
> *To: *gmx-developers at gromacs.org <mailto:gmx-developers at gromacs.org>
> *Subject: *[gmx-developers] Question about ewald_LRcorrection()
> function
>
> Dear GMX developers:
>
> I find a abnormal phenomenon when using
> *gromacs-4.6.7*. 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. *Whole system seems
> like below*:
>
> ****************************************
>
> * - H + *
>
> * - \ +
> vacuum *
>
> * - O + *
>
> * - / + *
>
> * - H
> + *
>
> ****************************************
>
> I use *Verlet cut-off scheme, PME method, and
> ewald-geometry=3DC*. Temperature of whole system is at *0K*,
> controlled by *v-rescale method*. I run whole simulation *only in
> OPENMP parallel*.
>
> One can have following expected behavior of water
> molecule without MD simulation: Since it can be seen as an
> infinite panel capacitor, *there is constant electric field
> between two walls. So the water will change it orientation, but
> won’t change it’s position(Since 0K)*.
>
> This phenomenon can be *confirmed by gromcas-3.3.1*
> (although there is position shift, but quite small). Also, if one
> use higher version of gromacs, such as 4.6.7 or higher *with
> ewald-geometry=3D, either no problem*. However, if one use *3DC*
> (and it should be), you can see *an obvious force on water molecule*.
>
> I hacked into source code, and find this might due
> to *ewald_LRcorrection()* function:
>
> Firstly, in gmx-4.6.7, ewald_LRcorrection()*.* We
> *loop over excluded neighbours only in condition calc_excl_corr=0.
> But when using verlet cut-off scheme, calc_excl_corr will be 1*.
> In gmx-3.3.1, we do loop over excluded neighbours without
> additional if() statement.
>
> Secondly, in gmx-4.6.7, , ewald_LRcorrection().
> When do dipole correction on force, we loop over *excl_load*
> (Exclusion load distribution over the threads). In fact, most
> water model will have nrexel=2. So, as a result, we *only do
> dipole correction on O atom, but do nothing on H atoms* and we
> will see an obvious force on the whole water molecule. I think
> this is not purpose of 3DC。
>
> Best regards！
>
> Sheng Bi
>
> **
>
> ---
> PhD student
> State Key Laboratory of Coal Combustion,
> School of Energy and Power Engineering,
> Huazhong University of Science and Technology (HUST),
> Room 302 Power Building,
> 1037 Luoyu Road, Wuhan, Hubei 430074 China
> Email: chrishengbee at hust.edu.cn <mailto:chrishengbee at hust.edu.cn>
> Phone: (86)27-87548122
> http://itp.energy.hust.edu.cn
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20160608/999d3380/attachment-0001.html>
More information about the gromacs.org_gmx-developers
mailing list