I’m exploring options for estimating molecular volumes. One method is, essentially, to calculate density on a 3D grid and then count the number of points that are above a cutoff. I know I could do this by outputting a cube file of the total density and then work on that but is there a simple way to do so without going through an intermediate cube file (and the associated IO operations)?
This would be through the python interface. The mechanics of how to do this are clearly there but are they accessible through the python API? I can see that there is compute_esp_over_grid_in_memory() but I’m after just the density on a 3D “cubic” grid.
IDK about psi4’s capabilities for in-memory features for this. However, to avoid IO it is always possible to make a RAM disk (tmpfs).
Do you have experience if a standard cubic grid used for visualization is accurate enough for this purpose?
Ooo… tmpfs might be the simplest way to go. Thanks for the suggestion.
In terms of how well it works, good question! There is literature on using this method: https://onlinelibrary.wiley.com/doi/10.1002/jcc.540160312 (I believe this is how the “volume” keyword in Gaussian works.)
I did have a play about doing this a while back and needed to use quite fine grids for it, hence not wanting to write to disk and back. I’m exploring whether there’s anything to be learned from comparing molecular volumes of a range of different conformers of a set of molecules. I.e., I don’t really care how well the calculated values match experimental values as long as I can see a difference between different conformers. No idea whether it will work or not. Yet!
I used the Monte Carlo method for volumes with simple vdW radii and it gets slow for larger molecules. At least my implementation
In principle it would be straightforward to use the optimized DFT grids instead of the simple regular grid of the cube files. (https://github.com/psi4/psi4numpy/blob/master/Tutorials/04_Density_Functional_Theory/4a_DFT_Grid.ipynb)
The freely available Multiwfn program (http://sobereva.com/multiwfn) fully supports the fchk file generated by PSI4, using Multiwfn to calculate molecular volume is fairly easy and quick. Boot up Multiwfn, then input commands like below
h2o.fchk // Example file
12 // Main function 12: Quantitative analysis of molecular surface
6 // Start analysis without considering mapped function
Then you will see below output:
Volume: 177.09406 Bohr^3 ( 26.24264 Angstrom^3)
Estimated density according to mass and volume (M/V): 1.1399 g/cm^3
Overall surface area: 154.10915 Bohr^2 ( 43.15496 Angstrom^2)
Note that the volume estimated in this way (see J. Mol. Graph. Model., 38, 314 (2012) for detail) has evidently higher numerical accuracy than the “volume” keyword in Gaussian, which employs Monte Carlo method.
By default, isosurface of electron density of 0.001 a.u. is used for defining molecular surface, you can alter this setting by using option “1 Select the way to define surface” in main function 12 of Multiwfn.
Multiwfn also supports calculating volume via Monte Carlo method based on molecular structure and atomic van der Waals radii (via subfunction 3 of main function 100).
If you simply need to calculate grid data, you can use main function 5 of Multiwfn. Multiwfn supports more than 100 real space functions, including electron density. After calculation, the grid data can be exported to cube file via corresponding option in the menu.
Ooo… Good suggestion. I’ve use Multiwfn in the past for various things but hadn’t thought to try it for volumes. What I’m really after is an algorithm that I can incorporate into my own code because I’m planning on running it against several hundred (if not thousands!) of different structures to build up a vague sort of database for later analysis of some other data. I’ll have a play…
I think this is a good item for our devs todo list.
(I assume you mean the in-RAM access to cubic grid data rather than volumes!)
Currently, I’m generating cubefiles of density, writing to (RAM) disk, then reading them back in to a numpy array for analysis. It would vastly simplify things if the gridded data were available as, e.g., a numpy array. I might have a look into that myself but I have not yet looked into how the C++ parts interface with the python interface so I probably won’t get that far. I’m sure it will be a good learning experience, though.
A much simpler speed up (for me, at least) would be being able to request just the total density as a cubefile, rather than alpha, beta, total, and spin density all being created. That could readily be achieved by adding another option to CUBEPROP_TASKS (TOTAL_DENSITY?). Most of the DENSITY parts can be used but only one out of four cubefiles would actually be written to disk (or available through the python interface!). I might have a go at implementing that in the next few days…
The volume part could also be added I think.
You would want to look at the libcubeprop directory for the C++ part.
The computations are done in-core and then written to a file (e.g csg.cc line 593 and line l231).
The forum does not get full attention from all the core devs. You can join the PSI4 slack (see github page) for a more immediate response on questions or open up an early work-in-progress pull request to get code pointers and opinions.