I would be very grateful if someone could recommend a way to plot RDFs for NPT ensemble molecular dynamics simulation results.

I use CP2K to run NVT simulations and usually use VMD to process results, but I don't see a way in VMD to deal with NPT results, where the volume is different at each time step. I'm hoping that someone can recommend a piece of software that can deal with my CP2K output (i.e., a position file in terms of atomic xyz coordinates and a file that contains volume output) on a Mac or Linux system.

I have spent several hours online looking for a solution, but haven't even come across discussions of this problem. This seems odd so maybe I'm using the wrong search terms. If anyone can point me to online discussions of this issue, that would be great.

I'm also not having much luck tagging this question! If anyone knows of any molecular dynamics-related tags I could use to reach people or perhaps of a more appropriate StackExchange site, I would welcome any suggestions. Thank you.

  • $\begingroup$ What happens if you attempt to use the existing processing tool to analyze the NPT data (pretending it's NVT)? $\endgroup$ – Buck Thorn Aug 26 '19 at 8:28
  • $\begingroup$ @Buck Thorn, thanks. Let me know if my explanation isn't clear and I'll rewrite it. VMD is set up to deal with NVT-type simulations where the volume is constant. There's a GUI where a sinlge volume value can be entered. So there's no way that I can see of importing my cell volume file and using a different volume value at each time step. however, VMD has a sort of coding language and can be used via a terminal. So I am wondering if there's a way to change the volume that way. I have so far averaged my volumes for each simulation and have used those, but wonder if my RDFS are meaningful. $\endgroup$ – Antst Aug 26 '19 at 9:34
  • $\begingroup$ @Buck Thorn, there don't seem to be any other tags that would be appropriate. "Molecular dynamics" and "ensemble" and related terms don't seem to be established as tags. $\endgroup$ – Antst Aug 26 '19 at 9:53
  • $\begingroup$ That's an important clarification. I'm guessing you mean this gui: ks.uiuc.edu/Research/vmd/plugins/gofrgui. The gui invokes the tcl command "measure gofr" with additional modifiers as described here: ks.uiuc.edu/Research/vmd/current/ug/node138.html While somebody else might have the script you are looking for I can recommend familiarizing yourself with the tcl scripting language in vmd as not everything is "canned" and you may want to process using your homebrew macros. $\endgroup$ – Buck Thorn Aug 26 '19 at 10:37
  • $\begingroup$ I don't know how gofr works but you can try inspecting it. It probably sets the max r cutoff based on the box size (or volume I guess). In NPT the volume will vary but you may get away with using the smallest volume sampled in the ensemble from which rmax is computed, but that is just a guess. $\endgroup$ – Buck Thorn Aug 26 '19 at 10:40

Thank you, Buck, for all your suggestions. You're right that a good way to check the influence of the volume is to try varying it. I did that by using the maximum and minimum volumes of the fluctuations in one simulation. There was some variation, so the volume does matter!

I was lucky enough to hear about "travis," an amazingly simple-to-compile code (even for Macs) available online for free. I thought I'd share that information here because other people are bound to run into this same problem. This code, travis, takes my trajectory file and my volume file and then asks a number of questions interactively. It outputs a nice CSV file that can be plotted.

Thanks again for your help, Buck. Really helped me think this through. I'd like to ask the moderators or whoever is in charge of tagging to please add some molecular dynamics-related tags if possible! Thank you.

  • $\begingroup$ It's cool that you found this software, that looks very useful for dissecting your simulations. I would be careful interpreting the effect of altering the volume setting in the gofr gui. I found the associated c routines which it presumably calls, you'll find them in the VMD source code (Measure.C). $\endgroup$ – Buck Thorn Aug 28 '19 at 8:56
  • $\begingroup$ regarding "molecular dynamics" tags, I took the freedom to add one to your question yesterday. $\endgroup$ – Buck Thorn Aug 28 '19 at 9:03
  • $\begingroup$ Reading the help page for "measure", at ks.uiuc.edu/Research/vmd/current/ug/node138.html, it says: "The size of the unitcell has to be stored in the trajectory file or has to be set manually for all frames with the molinfo command. The command uses by default only the current active frame for both selections." Not sure how this impacts your analysis. I take back a previous comment regarding pbcs etc as "usepbc=0" looks to be the default. I have to say again that you may get more mileage from posting a question at the vmd users forum. $\endgroup$ – Buck Thorn Aug 28 '19 at 9:49
  • $\begingroup$ Thanks a lot, Buck! That's good advice. I wasn't aware that there are VMD C code routines. They'll be easier for me to understand than the VMD tcl script. $\endgroup$ – Antst Aug 28 '19 at 13:12
  • $\begingroup$ Thank you also for adding the MD tags! I think your comment above means that in future, I need to set something in CP2K before running it that includes the effects of a varying unit call in my trajectory output. I think I know how to do this now. Unfortunately, I have many simulations that I can't redo, so I'm over the moon to know about travis. I was very worried I would not be able to use my existing simulations. Anyway, thank you so much for your help. Really calmed me down and put me on the right track. $\endgroup$ – Antst Aug 28 '19 at 13:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.