From: Axel Kohlmeyer (
Date: Tue Oct 09 2007 - 16:40:25 CDT

On Tue, 9 Oct 2007, Jay Shore wrote:

JS> Axel,


JS> Thank you for your comments and the copy of trimcube. They were very
JS> helpful.
JS> I guess I did not explain my mapping idea well enough. I want to take the
JS> isosurfaces as produced by VMD and then integrate the volume to determine
JS> what percent of the electron density is contained inside. I will have to
JS> decide how many surfaces to do this to. I can then write a new cube file
JS> with values based on what percent of electron density is composed in each
JS> surface.

ok. even though i still think it is not worth the effort,
here is a (crude) way of achieving what you want. i would
not recommend doing this inside of VMD for efficiency reasons,
but use trimcube as a template. please note, that this would only
work for cubefiles with only one dataset (orbcube=1), but that
kind of analysis only makes sense for density cube files anyways.

a) create an index map array that contains all indexes of the
grid points and fill it with an 1:1 mapping. i.e. something like:

int *map;
float normdata;
map = (int *)malloc(orbcube*nrdat*sizeof(int);
normdata =0.0;

for(i=0; i<nrdat*orbcube; ++i) {
   map[i] = i;
   normadata += data[i];

and then sort the dataset by swapping the indexes around.
so that when you you'd loop over data[map[i]] you would
loop from the largest value in the data array to the smallest.

c) now you can do:

float intsum;
for(i=0; i<nrdat*orbcube; ++i) {
   intsum += data[map[i]];
   data[map[i]] = intsum/normdata * 100.0;

that should give you an approximation of integrating the
data set 'from the inside out', which should be accurate enough
for visualization purposes.


JS> I have started to look at the VMD code to try to see how surfaces are
JS> defined. I will try to use trimcube for the integration part.
JS> Thanks again.
JS> Jay
JS> > On Sun, 7 Oct 2007, Jay Shore wrote:
JS> >
JS> > JS> Howdy,
JS> >
JS> > hi jay,
JS> >
JS> > JS> I am creating movies using VMD to incorporate into my General Chemistry
JS> > and
JS> > JS> physical chemistry lectures. I want to show the students electron density
JS> > JS> surfaces colored according to the electrostatic potential, but I am not
JS> > sure
JS> > JS> what is the most appropriate Isovalue to choose for the surface. It would
JS> > be
JS> >
JS> > the main purpose for mapping the electrostatic potential on the
JS> > electron density is to show the "polarity" of a molecule and how
JS> > it is perceived by other molecules. that means you want to pick an
JS> > isovalue that shows the "outside" of a molecule, i.e. a very low one.
JS> > i usually pick a value that is just large enough to still show a smooth
JS> > surface and the general shape of the molecule. you will notice, that
JS> > the colorization will not change significantly for as long as your
JS> > isovalue is small (no big surprise, there is little electron density
JS> > left).
JS> >
JS> > JS> nice to know what percent of the electron density was encompassed by a
JS> > JS> surface. To do this I think that I would have to write a program to
JS> > JS> integrate the volume inside the surfaces.
JS> >
JS> > to illustrate those porportions, you are probably better off using
JS> > a hydrogen molecule and a two dimensional graph. for most of the
JS> > time meaningful iso-surfaces will contain almost all electron density.
JS> > most of the electron density is contained in the inner shell electrons
JS> > anyways for almost all atoms/molecules.
JS> >
JS> > [...]
JS> >
JS> > JS> So, my questions:
JS> > JS>
JS> > JS> Am I missing an application that already exists that will convert my cube
JS> > JS> files into something resembling percent of electron density? I have done
JS> > a
JS> > JS> lot of searching (googling) but to no avail.
JS> >
JS> > i am not aware of any code to do what you need, but you may find the
JS> > attached code useful. it purpose is to trim away empty parts of a cube
JS> > file (and thus reduce memory consumption), but also allows to normalize
JS> > electrostatic potential maps from pseudopotential calculations.
JS> > it uses simple summing over the grid points as integration, which proved
JS> > to be accurate enough for visualization purposes. it should be
JS> > straightforward to modify it to your needs.
JS> >
JS> >
JS> > JS> If I do convert the value of the cube file to percent electron density
JS> > JS> contained (by integrating the volume of the surface) can I still use VMD
JS> > to
JS> > JS> draw the surfaces. I don't see why not, but maybe I am missing something.
JS> >
JS> i don't quite see how this mapping would work. to get the percentage of
JS> contained electron density you'll have to integrate and that will change
JS> the shape of the isosurface. particularly, if you compare multiple
JS> molecules.
JS> >
JS> > JS> What is the most meaningful way of representing an electron density
JS> > surface?
JS> > JS> Percent electron density contained or the atoms in molecules approach
JS> > (e.g.
JS> > JS> Bader). I know that it depends on exactly what I want to show, but is one
JS> > JS> method preferred over another?
JS> >
JS> > bader analyis is something different, since it computes the
JS> > _partitioning_ of the electron density. the integration of the
JS> > partitions and comparison to core charge leads numbers of how
JS> > atoms are polarized.
JS> >
JS> > JS> If I decide to write the code to do the integration and conversion, is
JS> > there
JS> > JS> a good starting point? I assume that I should do it in C or Fortran
JS> > because
JS> > JS> of the size of the cube files and a script would take too long.
JS> >
JS> > if you want to write anything, that could potentially be integrated into
JS> > VMD, i recommend writing in C/C++. that will make it easier. there are
JS> > some recent efforts to improve support for volumetric data sets akin
JS> > to what you need, but i have not followed that development closely, so
JS> > i cannot tell you how much is already available in current test release
JS> > of VMD.
JS> >
JS> > hope this helps.
JS> >
JS> > cheers,
JS> > axel.
JS> >
JS> > JS>
JS> > JS> I would appreciate any comments or suggestions.
JS> > JS> Thanks,
JS> > JS> Jay
JS> > JS>
JS> > JS>
JS> > JS>
JS> > JS>

Axel Kohlmeyer
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
If you make something idiot-proof, the universe creates a better idiot.