VMD-L Mailing List

From: Hall, Lisa Michelle (lhall_at_sandia.gov)
Date: Thu May 06 2010 - 18:12:49 CDT

Hi all,

I want to share a script in case it is useful to others, and also in case I might get some comments on how to improve it.

In the liquid state for non-biomolecules such as polymers, you often don't want to include intramolecular correlations in g(r). You will get a big spike, for instance, at the average distance of a backbone bond, and that doesn't tell you anything that you didn't already know but could obscure other interesting things. There's no easy way to exclude intramolecular terms from g(r) in measure gofr, which has been brought up on this list before. I want to use VMD to calculate g(r) because, at least when I do "measure gofr" for the full g(r) including intramolecular correlations, it is much, much faster than my own python script. This is partly because VMD doesn't have to think about the intramolecular connectivity, partly because VMD automatically does threading and is otherwise pretty efficient about things, and partly because my scripts are not efficient.

So, I've written a tcl script which selects a certain type of atom on a particular 'residue', and another type of atom on all other residues, (my residues are my polymer molecules) and uses measure gofr between those selections. [Actually, my types aren't read in correctly for now so it is selecting the types by their different charges of 0 and -1.] Then it loops over all of the first 100 residues summing up the individual molecule - other molecules g(r)s, and renormalizes at the end. Note that if you feed measure gofr selections of residue 1 and residue !=1, and there are N residues, each of n atoms, it will normalize g(r) by 1/(n*1)/(n*(N-1)). However, just as you would include the atom at the origin when normalizing an atomic g(r), you do want to include the molecule at the origin. So after summing over the residues, I divide by N and multiply by (N-1)/N so that overall it is effectively normalized by 1/(n*n*N*N).

To use the script, I put it in my home folder and type "play molec_gofr.tcl" into VMD. It isn't as slow as you would think; even though it has to call measure gofr 300 times, it is calling it on small selections.

Since I don't actually have any idea how to write a tcl script, and am not very proficient in VMD, it is not general and anyone else who wants to use it would need to modify it for their own systems. But it would give you somewhere to start if you want to do the same thing. If anyone else knows a better way (short of learning to write my python code to be more efficient and include MPI), I'd love to hear it.

Lisa