<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=Windows-1252">
</head>
<body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;">
So i did as you advised previously and what i found is that everything diverges when comparing the serial case vs the parallel one. I checked ke and pe. It takes several thousand integration steps (1.5fs time step) but the differences in these values become
quite large. The values i am comparing are <span style="font-family: Menlo; font-size: 11px;">enerd->term[F_EPOT] and </span><font face="Menlo"><span style="font-size: 11px;">enerd->term[F_EKIN]. These are only compared on an energy calculation step. During
a certain part of the simulation i change ir->nstcalcenergy equal to 1 such that it will calculate energy every step. When I’m done i set </span></font><span style="font-family: Menlo; font-size: 11px;">ir->nstcalcenergy back to its original value. </span><font face="Menlo"><span style="font-size: 11px;">Furthermore,
i am scaling v(t-1/2dt). My equation for computing scaling_factor considers the current forces f(t), v(t-1/2dt),pe(t) and what v(t+1/2dt) will be after scaling. This goes as follows</span></font>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">Target = pe(t) + ke^(t)</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">here ke^ implies that velocity scaling has already been done and Target is just some target energy we hope to have after scaling v(t-1/2dt) and computing v(t+1/2dt).</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">with leapfrog we have ke(t) = .5ke(t-1/2dt) + .5ke(t+1/2dt). Thus</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">Target - pe(t) = </span></font><span style="font-family: Menlo; font-size: 11px;">.5ke^(t-1/2dt) + .5ke^(t+1/2dt)</span></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">Target - pe(t) = 1/2*SUM[ 1/2*m*v^(t-1/2dt)*</span></font><span style="font-family: Menlo; font-size: 11px;">v^(t-1/2dt)</span><font face="Menlo"><span style="font-size: 11px;"> + 1/2*m*v^(t+1/2dt)**</span></font><span style="font-family: Menlo; font-size: 11px;">v^(t+1/2dt)
]</span></div>
<div><span style="font-family: Menlo; font-size: 11px;"><br>
</span></div>
<div><span style="font-family: Menlo; font-size: 11px;">Here SUM just means we are summing over all the atoms in the system. </span></div>
<div><span style="font-family: Menlo; font-size: 11px;"><br>
</span></div>
<div><font face="Menlo"><span style="font-size: 11px;">with leapfrog v(t+1/2dt) = v(t-1/2dt) + f(t)*dt/m </span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">Thus, we have the following</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><span style="font-family: Menlo; font-size: 11px;">Target - pe(t) = 1/4*SUM[ </span><font face="Menlo"><span style="font-size: 11px;">m*v^(t-1/2dt)*</span></font><span style="font-family: Menlo; font-size: 11px;">v^(t-1/2dt) ] + 1/4*SUM[ m*(v^(t-1/2dt)
+ f(t)*dt/mi)*</span><span style="font-family: Menlo; font-size: 11px;">(v^(t-1/2dt) + f(t)*dt/mi) ]</span></div>
<div><span style="font-family: Menlo; font-size: 11px;"><br>
</span></div>
<div><span style="font-family: Menlo; font-size: 11px;">Now we substitute v^(t-1/2dt) for K*v(t-1/2dt). Here v(t-1/2dt) is the original unscaled velocity and k is the scaling factor we are trying to solve for. This gives</span></div>
<div><span style="font-family: Menlo; font-size: 11px;"><br>
</span></div>
<div>
<div><span style="font-family: Menlo; font-size: 11px;">Target - pe(t) = 1/4*K*K*SUM[ </span><font face="Menlo"><span style="font-size: 11px;">m*v(t-1/2dt)*</span></font><span style="font-family: Menlo; font-size: 11px;">v(t-1/2dt) ] + 1/4*SUM[ m*(K*v(t-1/2dt)
+ f(t)*dt/mi)*</span><span style="font-family: Menlo; font-size: 11px;">(K*v(t-1/2dt) + f(t)*dt/mi) ]</span></div>
</div>
<div><span style="font-family: Menlo; font-size: 11px;"><br>
</span></div>
<div><font face="Menlo"><span style="font-size: 11px;">The first term on the right hand side is just K*K*ke(t-1/2dt). The second term is more troublesome. To get things in terms of K on the right hand side we now must perform a vector sum between f(t) and v(t-1/2dt).
We then perform a dot product of the resulting vector. And, we can’t forget that we have a summation sign in front of all of this. If you do the math you end up with terms that have K*K*(a number) + K*(another number) + (yet another number) = 0. </span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">So to solve for K we merely solve quadratic formula. This is easy to do with a computer. We simply put in a few loops etc. </span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">So… i do this. I let gromacs get past do_forec() such that we have pe(t) and f(t). Then, I run a loop over all the atoms and solve for ke(t-1/2dt). I then solve for k. Next, I multiply every component of
my velocities v(t-1/2dt) by this factor. I also multiply my previously computed ke(t-1/2dt) by this factor twice to get the new ke(t-1/2dt). Then Gromacs does what it does. It computes v(t+1/2dt). With v(t+1/2dt), i run a loop over all the atoms again and
compute ke(t+1/2dt). And, finally Etot(t) = pe(t) + .5ke(t-1/2dt) + .5ke(t+1/2dt) and that value should equal Target. </span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">I realize that more efficient means of doing this are possible. However, this only needs to be done rarely. Say once every 20,000 steps or so. Thus, i am not concerned with that right now. Futhermore, when
i do this using single processor it works perfectly. But, my problem is that when using multiple processors it fails. </span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">Any further advices is greatly appreciated.</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">Thanks</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">Nathan</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font>
<div>
<div>On Mar 19, 2015, at 3:26 AM, Berk Hess <<a href="mailto:hess@kth.se">hess@kth.se</a>> wrote:</div>
<br class="Apple-interchange-newline">
<blockquote type="cite">
<div bgcolor="#FFFFFF" text="#000000">
<div class="moz-cite-prefix">Hi,<br>
<br>
You can simply print the kinetic energy and scaling factor you calculate and compare a serial against a parallel run and you will see where things go off.<br>
<br>
But calculating kinetic energies and scaling is not trivial in MD. With the leap-frog integrator the reported kinetic energy is the average of that at t-0.5*delta_t and t+0.5*delta_t. And then the result will also depend when during the step you scale the velocities
(before/after which half-step Ekin calculation).<br>
<br>
Collecting the state is a very expensive operation that you should avoid. Better is to calculate the rank contribution to Ekin and do a MPI_Allreduce on that. But since mdrun already does this for energy reporting and temperature scaling, why not simply use
the Ekin mdrun calculates already? I would suggest to look at one of the temperature coupling algorithms, e.g. v-rescale, and reuse most of that machinery.<br>
<br>
Cheers,<br>
<br>
Berk<br>
<br>
On 03/19/2015 09:01 AM, Bernhardt, Nathan A. wrote:<br>
</div>
<blockquote cite="mid:FB2E3680-B205-44DE-BD2E-58EE16A1F292@ou.edu" type="cite">Hello everybody,
<div><br>
</div>
<div>This is my first time to post a message so feel free to comment accordingly. With that said however, let me get right to it.</div>
<div><br>
</div>
<div>i am performing md simulations (GROMACS 4.6.5) using the leap frog integrator. At this point i am not using LINCS but hope to be able to later. Lets just assume that at some point in the simulation i need to scale velocities such that the total energy
after scaling is equal to some target value. I can do this easily so long as a single processor is used in the simulation. Unfortunately, if i attempt to use parallel processing (domain decomp) my equation for calculating the scaling factor seems to no longer
work. This should not be the case and leads me to believe i may have misunderstood how gromacs is working.</div>
<div><br>
</div>
<div>I have tried two ways of doing this when domain decomp is used. the first uses dd_partition_system() and is as follows.</div>
<div><br>
</div>
<div> <span style="font-family: Menlo; font-size: 11px; color:
rgb(206, 121, 36);">if</span><span style="font-family: Menlo;
font-size: 11px;">(DOMAINDECOMP(cr))</span></div>
<div><span style="font-family: Menlo; font-size: 11px;">{</span></div>
<div><span style="font-family: Menlo; font-size: 11px;"><span class="Apple-tab-span" style="white-space:pre"></span></span><span style="font-family: Menlo; font-size: 11px;">dd_collect_state(cr->dd, state, state_global);</span></div>
<div><font face="Menlo"><span style="font-size: 11px;">}</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">if(group_rank == 0)</span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">{</span></font></div>
<div><font face="Menlo"><span class="Apple-tab-span" style="font-size: 11px; white-space: pre;"></span><span style="font-size: 11px;">calc_scaling_factor(state_global, energy_data); //this function calculates the scaling factor. i believe it is correct. </span></font></div>
<div><font face="Menlo"><span style="font-size: 11px;">}</span></font></div>
<div><br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><span style="color: rgb(206, 121, 36);">for</span>(i=<span style="color: rgb(195, 55, 32);">0</span>; i<nr_atoms; i++)</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;">{</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><span class="Apple-tab-span" style="white-space:pre"></span>state_global->v[i][<span style="color: rgb(195, 55, 32);">0</span>] = scaling_factor*state_global->v[i][<span style="color:
rgb(195, 55, 32);">0</span>];</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><span class="Apple-tab-span" style="white-space:pre"></span>state_global->v[i][<font color="#c33720">1</font>] = scaling_factor*state_global->v[i][<font color="#c33720">1</font>];</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><span class="Apple-tab-span" style="white-space:pre"></span>state_global->v[i][<font color="#c33720">2</font>] = scaling_factor*state_global->v[i][<font color="#c33720">2</font>];</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;">}</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;">if(DOMAINDECOMP(cr))</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;">{</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><span class="Apple-tab-span" style="white-space:pre"></span>dd_partition_system(fplog, step, cr, TRUE,
<span style="color: rgb(195, 55, 32);">1</span>,</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"> <span class="Apple-tab-span" style="white-space:pre">
</span>state_global, top_global, ir,</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"> state, &f, mdatoms, top, fr,</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"> vsite, shellfc, constr,</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"> nrnb, wcycle, FALSE);<br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;">}</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;">I will admit i do not fully understand what all dd_partition_system() does and the parameters it takes in. Thus, i tried to avoid using it by simply having each node scale the velocities for its
home atoms. Below is sample code from that attempt.</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;">
<div style="font-family: Helvetica; font-size: 12px;"> <span style="font-family: Menlo; font-size: 11px; color: rgb(206,
121, 36);">if</span><span style="font-family: Menlo;
font-size: 11px;">(DOMAINDECOMP(cr))</span></div>
<div style="font-family: Helvetica; font-size: 12px;"><span style="font-family: Menlo; font-size: 11px;">{</span></div>
<div style="font-family: Helvetica; font-size: 12px;"><span style="font-family: Menlo; font-size: 11px;"><span class="Apple-tab-span" style="white-space: pre;"></span></span><span style="font-family: Menlo; font-size: 11px;">dd_collect_state(cr->dd, state,
state_global);</span></div>
<div style="font-family: Helvetica; font-size: 12px;"><font face="Menlo"><span style="font-size: 11px;">}</span></font></div>
<div style="font-family: Helvetica; font-size: 12px;"><font face="Menlo"><span style="font-size: 11px;"><br>
</span></font></div>
<div style="font-family: Helvetica; font-size: 12px;"><font face="Menlo"><span style="font-size: 11px;">if(group_rank == 0)</span></font></div>
<div style="font-family: Helvetica; font-size: 12px;"><font face="Menlo"><span style="font-size: 11px;">{</span></font></div>
<div style="font-family: Helvetica; font-size: 12px;"><font face="Menlo"><span class="Apple-tab-span" style="font-size:
11px; white-space: pre;"></span><span style="font-size:
11px;">calc_scaling_factor(state_global, energy_data);
//this function calculates the scaling factor. i believe it is correct. </span></font></div>
<div style="font-family: Helvetica; font-size: 12px;"><font face="Menlo"><span style="font-size: 11px;">}</span></font></div>
<div><font face="Menlo"><br>
</font></div>
<div><font face="Menlo">broadcast_scaling_factor(scaling_factor); //This sends the scaling factor to all nodes so they can scale their home atoms velocities accordingly. </font></div>
<div><font face="Menlo"><br>
</font></div>
<div>
<div style="margin: 0px;">if(DOMAINDECOMP(cr))</div>
<div style="margin: 0px;">{</div>
</div>
<div style="margin: 0px;"><span class="Apple-tab-span" style="white-space:pre"></span><span style="color: rgb(206,
121, 36);">for</span>(i=<span style="color: rgb(195, 55,
32);">0</span>; i<mdatoms[<span style="color: rgb(195,
55, 32);">0</span>].homenr;
i++)</div>
<div style="margin: 0px;"><span class="Apple-tab-span" style="white-space:pre"></span>{</div>
<div style="margin: 0px;"><span class="Apple-tab-span" style="white-space:pre"></span>state->v[i][<span style="color: rgb(195, 55, 32);">0</span>] = scaling_factor*state->v[i][<span style="color: rgb(195, 55,
32);">0</span>];</div>
<div style="margin: 0px;"><span class="Apple-tab-span" style="white-space:pre"></span>state->v[i][<font color="#c33720">1</font>] = scaling_factor*state->v[i][<font color="#c33720">1</font>];</div>
<div style="margin: 0px;"><span class="Apple-tab-span" style="white-space:pre"></span>state->v[i][<font color="#c33720">2</font>] = scaling_factor*state->v[i][<font color="#c33720">2</font>];</div>
<div style="margin: 0px;"><span class="Apple-tab-span" style="white-space:pre"></span>}</div>
<div style="margin: 0px;">}</div>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;">Neither of these methods have proved successful. In both cases the total energy after scaling is the same and off slightly from the target value by about 10kj/mol. In the second method I have tested
that broadcasting of the scaling factor is successful. Since both methods require a global state to be collected and used for calculating the scaling factor i see no reason why the function that calculates this factor would work for the single processor case
and not the domain decomp one. This leads me to believe that either the dd_collect_state() function is not doing exactly what i think it is, the dd_partition_system() may be working in ways i do not understand or that state->v[][] may not be the array I need
to modify. Any help is appreciated. </div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><br>
</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;">Thanks</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;">Nathan Bernhardt</div>
<div style="margin: 0px; font-size: 11px; font-family: Menlo;"><br>
</div>
<br>
<fieldset class="mimeAttachmentHeader"></fieldset> <br>
</blockquote>
<br>
</div>
-- <br>
Gromacs Developers mailing list<br>
<br>
* Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List">
http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List</a> before posting!<br>
<br>
* Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists">http://www.gromacs.org/Support/Mailing_Lists</a><br>
<br>
* For (un)subscribe requests visit<br>
<a href="https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers">https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers</a> or send a mail to
<a href="mailto:gmx-developers-request@gromacs.org">gmx-developers-request@gromacs.org</a>.</blockquote>
</div>
<br>
</div>
</body>
</html>