<html><body><div style="color:#000; background-color:#fff; font-family:tahoma, new york, times, serif;font-size:14pt"><div><span>Dear </span>Julius,</div>
<div>That's really cool! Thank you very much for sharing this solution! <br></div><div>I thing mdlib/domdec.h and mdlib/domdec.c should be officially patched to allow such things.</div><div><br></div><div>Regards,</div><div>Semen<br></div><div><br><blockquote style="border-left: 2px solid rgb(16, 16, 255); margin-left: 5px; padding-left: 5px;"><div style="font-family: tahoma, new york, times, serif; font-size: 14pt;"><div style="font-family: times new roman, new york, times, serif; font-size: 12pt;"><font size="2" face="Arial"><hr size="1"><b><span style="font-weight:bold;">From:</span></b> Julius Su &lt;jsu@caltech.edu&gt;<br><b><span style="font-weight: bold;">To:</span></b> gmx-developers@gromacs.org<br><b><span style="font-weight: bold;">Sent:</span></b> Saturday, October 15, 2011 12:50 AM<br><b><span style="font-weight: bold;">Subject:</span></b> [gmx-developers] adding external forces to all atoms<br></font><br>
Hi everyone,<br><br>I would like to present a solution I've arrived at to implement a feature -- the ability to add a set of external forces F(x1 ... xN) to all the atoms of a parallel GROMACS simulation -- which has been previously discussed:<br><br>http://lists.gromacs.org/pipermail/gmx-developers/2011-May/005193.html<br>http://lists.gromacs.org/pipermail/gmx-developers/2011-April/005141.html<br><br>This could useful for applying biasing potentials that span decomposition domains, such as ones acting on collective coordinates, or for applying novel experimental restraints which include large groups of atoms in the system. There is existing functionality in GROMACS to do this sort of thing, but the code below makes it straightforward to try out different potentials, without needing to perform a more extensive integration with the existing code base. So it may be useful for prototyping purposes ..<br><br>Note that in the approach below, all the work for
 evaluating F(x1 ... xN) is localized to one processor (the Master), so it should not a time consuming potential to evaluate -- else the parallel benefits of GROMACS are negated. Also the synchronization required to evaluate and distribute the forces may cause the parallel performance to suffer slightly.<br><br>In any case, since there appears to be some (scattered) interest, I'll present my code modifications. They may be of use to others in the community. Also if I've done anything egregiously (or even slightly) wrong, please let me know!<br><br>Sincerely,<br><br>Julius<br><br>-----<br>(1) In kernel/md.c, modify such that f_global is always initialized:<br><br>&nbsp; &nbsp; &nbsp; &nbsp; //if (DDMASTER(cr-&gt;dd) &amp;&amp; ir-&gt;nstfout) {<br>&nbsp; &nbsp; &nbsp; &nbsp; if (1) {&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;  // change so that f_global is always initialized<br>&nbsp;
 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; snew(f_global,state_global-&gt;natoms);<br><br>(2) In kernel/md.c, add lines in<br><br>&nbsp; &nbsp; &nbsp;  /*&nbsp; ################## END TRAJECTORY OUTPUT ################ */<br><br>&nbsp; &nbsp; &nbsp;  // Collect all<br>&nbsp; &nbsp; &nbsp;  if (DOMAINDECOMP(cr))<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; dd_collect_vec(cr-&gt;dd, state, f, f_global);<br><br>&nbsp; &nbsp; &nbsp; &nbsp; if (MASTER(cr))<br>&nbsp; &nbsp; &nbsp; &nbsp; {<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;  // Fetch positions of all atoms<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; for (i = 0; i &lt; state_global-&gt;natoms; i++)<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; {<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; current_x[i] = state_global-&gt;x[i][0];<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; current_y[i] = state_global-&gt;x[i][1];<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; current_z[i] = state_global-&gt;x[i][2];<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;
 }<br><br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; // Compute biasing potential (ADD NEW POTENTIAL HERE)<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; ... Given current_x, current_y, current_z --&gt; returns delta_fx, delta_fy, delta_fz.<br><br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; // Apply biasing potential to all atoms<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; for (i = 0; i &lt; state_global-&gt;natoms; i++)<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; {<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;  f_global[i][0] += delta_fx[i];<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;  f_global[i][1] += delta_fy[i];<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;  f_global[i][2] += delta_fz[i];<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; }<br>&nbsp; &nbsp; &nbsp; &nbsp; }<br><br>&nbsp; &nbsp; &nbsp; &nbsp; // Distribute modified f_global back to individual processors<br>&nbsp; &nbsp; &nbsp; &nbsp; if (DOMAINDECOMP(cr))<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; dd_distribute_vec(cr-&gt;dd,
 dd_charge_groups_global(cr-&gt;dd), f_global, f);<br><br>&nbsp; &nbsp; &nbsp; &nbsp; // Make sure all processes are synchronized before continuing<br>&nbsp; &nbsp; &nbsp; &nbsp; MPI_Barrier(cr-&gt;mpi_comm_mygroup);<br><br>(3) In mdlib/domdec.c, modify dd_distribute_vec to have external scope:<br><br>//static void dd_distribute_vec(gmx_domdec_t *dd,t_block *cgs,rvec *v,rvec *lv)<br>void dd_distribute_vec(gmx_domdec_t *dd,t_block *cgs,rvec *v,rvec *lv) // allow external functions access<br><br>(4) In include/domdec.h, expose dd_distribute_vec for use by other functions:<br><br>void dd_distribute_vec(gmx_domdec_t *dd,t_block *cgs,rvec *v,rvec *lv);<br><br>----<br><br>-- gmx-developers mailing list<br><a ymailto="mailto:gmx-developers@gromacs.org" href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a><br>http://lists.gromacs.org/mailman/listinfo/gmx-developers<br>Please don't post (un)subscribe requests to the list. Use the www interface
 or send it to <a ymailto="mailto:gmx-developers-request@gromacs.org" href="mailto:gmx-developers-request@gromacs.org">gmx-developers-request@gromacs.org</a>.<br><br><br></div></div></blockquote></div></div></body></html>