VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_gmail.com)
Date: Thu May 06 2010 - 19:46:26 CDT

On Thu, 2010-05-06 at 17:12 -0600, Hall, Lisa Michelle wrote:
> Hi all,

hi lisa,

> 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

right. in the current implementation the code is completely
oblivious to molecular structure. it only sees atoms.

> 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.

thanks for the compliment about efficiency. this is code that has
been tweaked and translated from fortran 77 to c to c++ (which turned
out to be the most efficient) over a period of almost 15 years now.
it is a pretty compute and memory intense problem and thus will never
be fast in script languages. wait until you see the GPU accelerated
version that my colleague ben levine has written with help from john.

> 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).

i have to think about this for a bit. outside of possibly requiring
to write a second, more complex g(r) kernel for this special case,
this should be doable to implement. if i get you right, what you would
need is some kind of exclusion that states if two atoms have the
same residue/fragment/chain/... info, they should be excluded from
the g(r) calculation. i have not looked at the code for quite some
time, but i am planning to work on some additional analysis tools
for VMD over the summer, and this could be a good "warmup". there
are also several cases, where prefixing the calculation with a grid
binning, could result in large speed ups.

> 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.

to get a speedier python code, you might try the python to c translator
(if that still exists and works).

how soon would you (or others) need this to be fully implemented
in measure gofr? realistically?

i would envision to have an additional flag (-same-exclude
residue/fragment/chain ) that indicates the desired exclusion
strategy. normalization should be automatically correct with the
mechanism that is currently in place, but i'll have to double
check.

cheers,
    axel.
>
> Lisa 

-- 
Dr. Axel Kohlmeyer  akohlmey_at_gmail.com
http://sites.google.com/site/akohlmey/
Institute for Computational Molecular Science
Temple University, Philadelphia PA, USA.