<div dir="ltr"><div><div><div><div><div>Hi Justin,<br><br></div>Thanks for the recommendation, but the filter implementation is not really my question. I am trying to update mdrun to include this filtering so that it can be done on a signal with the proper sampling rate.<br><br>I'm thinking that I should update md.c with some sort of do_filter_update() and do_filter_writing() function called right before do_md_trajectory_writing(). I can maintain the value of the cumulative sum for the filter with arrays managed by md.c, and then I need to:<br>(1) make sure that each DD node computes only their own filtered positions independently every time step<br>(2) they communicate their results when appropriate<br></div>(3) cumulative sums for atoms move with the atoms when atoms move across domains (not sure how to do this)<br></div>(4) I add nstfilterxtcout or something like that to the options<br><br></div>Does this sound like a sufficiently canonical approach? <br><br>Just to be clear on my intention, low pass filtering is only useful/rigorously correct if the original signal being filtered is not polluted with aliasing, i.e. it is sampled at the nyquist frequency of the system, e.g. the sampled signal from every time step. I think it is a reasonable approximation to say that all of the high frequency noise is "white" noise relative to conformational dynamics and just do the filtering on the 1ps-10ps sampled frames, but that is not really rigorous. If it doesn't cost additional memory and only adds O(N) computation at each step, it seems like a useful thing to have in the program, and in particular a thing that I am quite interested in.<br><br></div><div>Thanks for all the timely responses.<br><br></div>John<br></div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Apr 12, 2017 at 3:25 PM, Justin Lemkul <span dir="ltr"><<a href="mailto:jalemkul@vt.edu" target="_blank">jalemkul@vt.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
It sounds like you may want to start by looking at gmx filter -ol and the underlying code.<br>
<br>
-Justin<span class=""><br>
<br>
On 4/12/17 1:50 PM, John Haberstroh wrote:<br>
</span><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class="">
Hello Mark,<br>
<br>
Thank you for the suggestion, but as you suggested it kind of misses the point<br>
of what I was asking. Maybe *I* misunderstand the -aver option in g_energy if<br>
that's the case. I'll give an example of what I want.<br>
<br>
I want an additional trajectory output that averages over all steps (every 2fs),<br>
but I want it to be output every 10ps. This should result in a fairly smooth<br>
trajectory without sacrificing temporal resolution below 10ps. I am referring to<br>
this averaging procedure as "lowpass" filtering, though the square window form<br>
would need to be more complicated to be strictly a lowpass filter.<br>
<br>
Alongside that, I want the instantaneous samples at every 10ps (i.e. the current<br>
.xtc output), which includes fluctuations with random phases from frequency<br>
components with periods much higher than 10ps, but in exchange obeys boltzmann<br>
statistics.<br>
<br>
If I could do it "in post", I would be fine, but I can't do it in post without<br>
writing out every single time step, which is far too resource intensive for my<br>
needs.<br>
<br>
Thanks,<br>
John<br>
<br>
ps: I was subscribed to digest so I apologize if this does not remain within the<br>
original thread, but I have fixed it now, so there should be no further issues.<br>
<br>
Date: Wed, 12 Apr 2017 16:55:01 +0100<br></span>
From: David van der Spoel <<a href="mailto:spoel@xray.bmc.uu.se" target="_blank">spoel@xray.bmc.uu.se</a> <mailto:<a href="mailto:spoel@xray.bmc.uu.se" target="_blank">spoel@xray.bmc.uu.se</a>>><br>
To: <a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs.org</a> <mailto:<a href="mailto:gmx-developers@gromacs.org" target="_blank">gmx-developers@gromacs<wbr>.org</a>><span class=""><br>
Subject: Re: [gmx-developers] Lowpass trajectory implementation<br>
Message-ID: <<a href="mailto:1d2ce3f6-2a0f-6137-e351-6a84b5eab007@xray.bmc.uu.se" target="_blank">1d2ce3f6-2a0f-6137-e351-6a84b<wbr>5eab007@xray.bmc.uu.se</a><br></span>
<mailto:<a href="mailto:1d2ce3f6-2a0f-6137-e351-6a84b5eab007@xray.bmc.uu.se" target="_blank">1d2ce3f6-2a0f-6137-e35<wbr>1-6a84b5eab007@xray.bmc.uu.se</a>><wbr>><div><div class="h5"><br>
Content-Type: text/plain; charset=windows-1252; format=flowed<br>
<br>
On 12/04/17 16:30, John Haberstroh wrote:<br>
> Hello developers,<br>
><br>
> I am interested in implementing a lowpass filter system for position and<br>
> dihedral trajectories akin to the -aver option in g_energy ("Also print<br>
> the exact average and rmsd stored in the energy frames (only when 1 term<br>
> is requested)") . I have found that a lot of computational effort goes<br>
> into eliminating aliased noise from very high frequencies in these<br>
> position coordinates, and it would be nice to turn this noise down a few<br>
> decibels. As it stands, it is not even possible to output dihedrals<br>
> during simulation though they are ostensibly be computed, let alone to<br>
> print moving averages alongside samples. That said, my first step would<br>
> be doing these averages for atomic positions as it would be (maybe)<br>
> pretty benign to extract dihedrals from the filtered atomic positions.<br>
><br>
> I understand that averaged positions will not necessarily be physical<br>
> configurations and do not obey boltzmann statistics. I am interested in<br>
> this procedure for the sake of getting a degree of spectral<br>
> separability. I also realize that the fluctuations will have to be<br>
> captured in a covariance matrix to be appropriately described. These<br>
> concerns are secondary to getting a signal with some semblance of low<br>
> pass filtering.<br>
><br>
> Any recommendations of where I should start for storing a rolling mean<br>
> of atomic positions? I can probably piggy back on the trajectory dump<br>
> procedures without issue if I can update this rolling mean every<br>
> timestep, so if you know where I can create a new array that is<br>
> persistent throughout the simulation, that would be the right kind of tip.<br>
><br>
I must say it is not clear to me what you would get out of this, but<br>
making a variant of gmx trjconv that reads the first N frames and writes<br>
out the time averaged coordinates to a new trajectory could do the<br>
trick. A couple of lines of code should suffice.<br>
<br>
<br>
<br>
> Thanks!<br>
> John<br>
><br>
> PS: I have edited gromacs source before, but only "brute forcing" a<br>
> modification to the pull.c function for alternative biased sampling,<br>
> which did not require generating any new persistent data variables. I do<br>
> have some experience with the gromacs internals, though!<br>
><br>
><br>
<br>
<br>
--<br>
David van der Spoel, Ph.D., Professor of Biology<br>
Head of Department, Cell & Molecular Biology, Uppsala University.<br></div></div>
Box 596, SE-75124 Uppsala, Sweden. Phone: <a href="tel:%2B46184714205" value="+46184714205" target="_blank">+46184714205</a> <tel:%2B46184714205>.<br>
<a href="http://www.icm.uu.se" rel="noreferrer" target="_blank">http://www.icm.uu.se</a><br>
<br>
<br>
<br>
<br>
</blockquote>
<br>
-- <br>
==============================<wbr>====================<br>
<br>
Justin A. Lemkul, Ph.D.<br>
Ruth L. Kirschstein NRSA Postdoctoral Fellow<br>
<br>
Department of Pharmaceutical Sciences<br>
School of Pharmacy<br>
Health Sciences Facility II, Room 629<br>
University of Maryland, Baltimore<br>
20 Penn St.<br>
Baltimore, MD 21201<br>
<br>
<a href="mailto:jalemkul@outerbanks.umaryland.edu" target="_blank">jalemkul@outerbanks.umaryland.<wbr>edu</a> | <a href="tel:%28410%29%20706-7441" value="+14107067441" target="_blank">(410) 706-7441</a><br>
<a href="http://mackerell.umaryland.edu/~jalemkul" rel="noreferrer" target="_blank">http://mackerell.umaryland.edu<wbr>/~jalemkul</a><br>
<br>
==============================<wbr>====================<span class="HOEnZb"><font color="#888888"><br>
-- <br>
Gromacs Developers mailing list<br>
<br>
* Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List" rel="noreferrer" target="_blank">http://www.gromacs.org/Support<wbr>/Mailing_Lists/GMX-developers_<wbr>List</a> before posting!<br>
<br>
* Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists" rel="noreferrer" target="_blank">http://www.gromacs.org/Support<wbr>/Mailing_Lists</a><br>
<br>
* For (un)subscribe requests visit<br>
<a href="https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers" rel="noreferrer" target="_blank">https://maillist.sys.kth.se/ma<wbr>ilman/listinfo/gromacs.org_gmx<wbr>-developers</a> or send a mail to <a href="mailto:gmx-developers-request@gromacs.org" target="_blank">gmx-developers-request@gromacs<wbr>.org</a>.<br>
</font></span></blockquote></div><br></div>