Subject: SUMMARY: Volume of a Molecule within an Isodensity Surface

From: "C.F. Matta" <[email-masked]>
Date: Sat, 3 Mar 2001 12:48:04 -0500 (EST)

Hi everyone,

A couple of weeks ago I posted the following question:

"   I wonder is someone knows of a program that calculates the volume
of a molecule contained within a particular density envelope (e.g. the
0.001 a.u. envelope) from a GAUSSIAN94 wavefunction or cube ? (The volume
I want to calculate is the INTEGRATED volume within the envelope of choice
AND NOT the volume computed from intersecting spheres).  "

Cherif F. Matta

	SUMMARY OF REPLIES:
	===================

Several very useful answers were received.  

Dr. Curt Breneman (Reply4) suggested a nice and simple algorithm,
Dr. Hoefinger (Reply7) sent a useful awk script that I tried and obtained
volumes that agree with those calculated from GAUSSIAN94 by invoking 
SCRF=IPCM as suggested by Dr. Doug Fox (Reply6).  I would like to thank
others who sent helpful comments and web pointers: Dr. N.O.J. Malcolm
(Reply1: MORPHY program);  Dr. Paulo Couto (Reply2: Moelekel 
program);  Dr. Laurence Cuffe (Reply3: "volume" keyword);  and  Dr. TJ
O'Donnell (Reply5: Codes to calculate molecular volumes and surfaces).

	INDIVIDUAL REPLIES:
	===================

Reply 1
======
From: noj malcolm <[email-masked]>

Our in-house code MORPHY (a general AIM program) does what you
require, for details see our web-pages
http://www.ch.umist.ac.uk/morphy/default.htm.

Dr. N.O.J. Malcolm                       [email-masked]
Chemistry Department
UMIST, Sackville Street
Manchester, England. M60 1QD
-------------------------------------------------------------------------------------------
Reply2
======
From: Paulo Couto <[email-masked]>

Moelekel visualization package allows you to
plot isodensity surfaces from a Gaussian or
Gamess output file. Furthermore, the volume
enclosed by a given isodensity surface is computed.
Please check http://www.cscs.ch/molekel/
-------------------------------------------------------------------------------------------
Reply3
======
From: Laurence Cuffe <[email-masked]>

I am probably mising something here but have you used the
"volume" keyword? It seems to do what you want via a monte carlo
integration procedure.
-------------------------------------------------------------------------------------------
Reply4
======
From: Curt M. Breneman <[email-masked]>

One thing we used to do in the Wiberg group (back in my late 80's
post-doc days) was to run the cube file through a little algorithm that
counted the number of points with greater than or equal to 0.001 (or
0.002 e/au3) density values, then multiply that number by the cubic
volume increment to get the volume.  This works pretty well,
particularly for small molecules and cubes that are around 81x81x81
points.  More coarse grids also work, but with slightly lower accuracy.

Curt Breneman
RPI Chemistry
-------------------------------------------------------------------------------------------
Reply5
======
From: tj <[email-masked]>
Reply-To: [email-masked]
Cc: Jim Kress <[email-masked]>

Most any 3D contouring software, for example those using the
marching cubes algorithm, should output the volume (and area) of
the "envelope".  I have some old code that can do this.  I
could dig it up if you can't find anything else.  It would most likely
have to be modified to read G94 files.  I have a vague recollection
that I had added the feature to read G9x cube files.
I had sent my code to Jim Kress who was using it with G94.
I haven't heard from him in a while.  I've copied him on this message
and you might try getting in touch with him. and you might try getting in
touch with him.
My agreement with him was that he could use the code, but not
distribute it.
I have no problem with you using my code, as long as proper
acknowledgement is given.
TJ O'Donnell
[email-masked]
-------------------------------------------------------------------------------------------
Reply6
======
From: Cust. Service Doug <[email-masked]>
Reply-To: [email-masked]

If you turn on #p running SCRF=IPCM you will get the integrated
volume and surface area from the isodensity contour specified.  Remember
that this is supplied after the structure, dielectric and then the
cutoff.  The volume is printed after the first SCF when it gets into
l117.

Douglas J. Fox
Technical Support
Gaussian, Inc.
[email-masked]
-------------------------------------------------------------------------------------------
Reply7
======
From: Siegfried Hoefinger <[email-masked]>

If you've got a *.wfn file, then you could easily get a *.qub file with
AIMPAC and the below attatched awk script should give you the volume.
Note that the file that specifies the *.qub is called tmp.inf below
and threshold becomes an isodensity value.  ( example shown below )

Siegfried Hoefinger, M.S.,Ph.D.      Universita' degli Studi di Bologna
phone: [0039]/051 2099498            Dipartimento di Chimica
"G. Ciamician"
fax:   [0039]/051 2099456            Via F. Selmi 2, I-40126, Bologna,
Italy
http://poseidon.gup.uni-linz.ac.at/sh/sh_frames.html

The awk script should be cited as:
**********************************

	Hoefinger, S. (2000):  Aimpacvol, a program to calculate molecular
        volumes, [unpublished, private communication].

Here is the script in full and notes on how to use it:
======================================================

1.) Create a *.qub with AIMPAC from a *.wfn file: the qub specifying file
    is called 'tmp.inf' and the *.wfn file 'out.wfn' and the resulting
    *.qub fiel will be termed 'tmp.qub'.
/SOME_PATH/AIMPAC/cubev tmp out
***********
* tmp.inf *
***********
  TITLE:  ch3co-/631+g
  INCRM:  10.0 0.3 0
  CENTR:  0.0 0.0 0.0
  PLANE:  0 1 2
  NFUNC:  1
#
#
# Specfies a box of 10 Bohr x 10 Bohr x 10 Bohr, centered at 0/0/0
# with increments of 0.3 Bohr, computing the electron density ( NFUNC )
# considering no symmetry and leaving atoms 1 and 2 in the xy plane
# of interest. The whole thing will result in a tmp.qub file.
# NOTE: use increments not smaller than 1/100 of the boxlength to
#       avoid core dumps.
#

2.) Estimate the volume with the help of awk script 'aimpacvol.awk'

*********************
*  aimpacvol.awk    *
*********************
#
#  PURPOSE: Counts the number of entries greater than a user defined
#           threshold value within a *.qub file produced from AIMPAC.
#           Based on this infinitesimal partial volume of one of these
#           grid cells, an estimation of the molecular volume inside
#           the specified threshold may be performed and the result
#           is presented in the end. In addition, all the coordinates
#           of these gridpoints having an electron density value
#           between exact the threshold and 15% above, are written
#           to another file called 'tmp.xyz' and may be used for
#           further processing, probably from awk-script 'fndclspnt.awk'.
#  USAGE:   awk -f aimpacvol.awk tmp.qub tmp.inf threshold
#  AUTHOR:  Siegfried Hoefinger, [email-masked]
#  DATE:    24.05.2000
#
BEGIN{
#       for( i=1; i<ARGC; i++ )
#         {
#            printf "%s %s \n", "Data file used is:", ARGV[i]
#         }
        THRESH=ARGV[3]
        CUBEINFFILE=ARGV[2]
        CUBEFILE=ARGV[1]
        getline < CUBEFILE
        NX=$1
        NY=$2
        NZ=$3
        getline < CUBEINFFILE
        getline < CUBEINFFILE
        LENGTH=$2
        INCR=$3
        getline < CUBEINFFILE
        CX=$2
        CY=$3
        CZ=$4
        CNT=0
        SUCCESSCNT=0
        STX=CX-(LENGTH/2.0)
        STY=CY-(LENGTH/2.0)
        STZ=CZ-(LENGTH/2.0)
        NXNY=NX*NY
        ARGC-=2
     }
{
  while( (getline < CUBEFILE) > 0 )
      {
         for( i=1; i<=NF; i++ )
         {
            CNT++
            IZ=int( (CNT-1)/NXNY )
            IY=int( (CNT-(IZ*NXNY)-1)/NX )
            IX=CNT-(IZ*NXNY)-(IY*NX)-1
            if( $i >= THRESH )
             {
                SUCCESSCNT++
                if( $i <= (THRESH*1.15) )
                 {
                    X=STX+(IX*INCR)
                    Y=STY+(IY*INCR)
                    Z=STZ+(IZ*INCR)
                    printf "%12.6lf %12.6lf %12.6lf\n", X, Y, Z > "tmp.xyz"
                 }
             }
         }
      }
}
END{
        printf "%10s %dx%dx%d %12s\n", "There were", NX, NY, NZ,"gridpoints"
        printf "%d %18s\n", CNT, "actually counted"
        printf "%d %18s\n", SUCCESSCNT, "inside threshold"
        printf "%12s %12.6lf %9s %12.6lf %7s\n", \
               "gridvolume::", INCR*INCR*INCR, " Bohr^3 =", \
               INCR*INCR*INCR*0.1481847435, " Angs^3"
        printf "%25s %12.6lf %9s %12.6lf %7s\n", "volume inside
threshold::",
               SUCCESSCNT*INCR*INCR*INCR, " Bohr^3 =", \
               SUCCESSCNT*INCR*INCR*INCR*0.1481847435, " Angs^3"

   }
-------------------------------------------------------------------------------------------

			END OF REPLIES
			==============

Cheers,

Cherif
___________________________________________________________________________

  Cherif F. Matta		    	  tel. (905) 525-9140 ext. 22502
  Chemistry Department                    fax  (905) 522-2509
  McMaster University                     
  Hamilton, Ontario, CANADA L8S 4M1.        
___________________________________________________________________________