More On The Virtues Of VirtualBox – ACID (or AICD) Under Ubuntu 14.04 (By Way Of OpenSuse 11.2)

“Stop that!” – George Carlin

If you’ve obtained source code from an academic lab that was last developed some time ago and you spent a whole day installing libraries and symbolic links and redefining variables in your .bashrc and downgrading libraries and redefining paths and have 20 tabs open in your browser that all go to 20 different obscure error discussions on Stack Overflow and it’s late and you’re tired and you think you might not need the program after all if you do a bunch of other workaround things instead – what’s below is for you.

Academics have been developing small code for (nearly) millions of years to make their lives easier – and we all benefit when that code is made available to others (esp. when it helps in data analysis). When that code is a series of perl or python scripts, there’s generally little reason why you should have any run issues. When they call on external libraries or specific tools, generally that information is available in the README somewhere. Generally speaking, there’s no reason why a code shouldn’t work in a straightforward manner when the developer doesn’t make it known that something else needs be installed to make it work.

So, why doesn’t code A work on your linux box? A few possibilities.

0. You didn’t read the README right.

1. The developer skipped some documentation and you’re missing a library or program that they installed months earlier but forgot about (and they forgot to acknowledge that that was a needed item to run their program).

2. Your linux distro is too new, they developed it on something much older, have never updated their own machine, and something fundamental changed between, for instance, Ubuntu 8.04 and 14.04 (I’ve run into a few in that category – depending on the code, this can sometimes just be chalked up to a compiler difference issue that went unnoticed by the developers).

3. The one you learn the harder way – Ubuntu is not Fedora is not OpenSuse is not [insert distro here]. Half of the Ubuntu posts on this website stem from my having needed to do something under Ubuntu to make it look a little tiny bit more like whatever distro a developer developed their code with.

Let’s hammer on point 3 for a minute – I know precious few academics who use computers as part of research who are both (1) experts in coding and (2) experts in interpreting the output of some code. A properly working code in the hands of a productive researcher can mean the difference between an OK interpretation of data and a very descriptive interpretation of data that anyone else could follow (which, let’s not lie to ourselves, may also mean the difference between a journal with an impact factor of 0.5 and one of 5.0). Disregarding any “it was hard to code, it should be hard to use!” mentality, finding out that a program runs painlessly on distro A with NO effort after trying in vain on distro B with CONSIDERABLE effort is a rather important piece of information.

For the record – this isn’t about a user’s ability to program, install, admin their machine, etc. It’s simply about ease-of-use for ALL users regardless of their computational aptitude (which is why many of the Ubuntu posts on this website are exhaustive in their reporting of errors as well as work-arounds).

Which brings me to one of the most versatile (and storage-consuming) programs on my Ubuntu 14.04 machine – VirtualBox. VirtualBox (and related) are absolutely awesome for their abilities to let you try – very quickly – the quick fix that is the building or running of a strange code “somewhere else.” For the record, I ALWAYS try something new using a clean Ubuntu 14.04 image under VirtualBox before I ever let it touch my main machine – and that’s because I’ve had more than one instance in the past of trying to build a piece of software that then destroyed (well, sort of) my machine. When you destroy an image with a bad installation attempt – all you have to do is trash the image and copy a fresh one over from some “clean image” folder. Could not be easier. Furthermore, I’ve still ONE app that works in Windows XP that doesn’t work anywhere else. And I’ve two or three under Windows that I use generally – all ready for use under VirtualBox without having to reboot anything.

Additionally importantly – Virtual Machines allow you to “sample in time” – working around new & improved features and functionality in new compilers and libraries by instead building code on what was available at the time of the code’s release.

Anisotropy of the Current (Induced) Density – ACID

2016mar29_benzol.nmr_40000_0.050_1_0_0_Aniso_4

Sample output (from a successful Gaussian09 benzene tutorial run).

From the website of Prof. Dr. Rainer Herges

The “density of delocalized electrons” is a concept that is intuitively used to explain the electronic structure of conjugated systems. Unfortunately, however, there is no rigorous way to separate the total electron density in a density of localized and delocalized electrons. Like aromaticity, bond order, point charge, and other important concepts of chemistry, a definition for the density of delocalized electrons has to be derived from more fundamental quantum theoretical parameters.

As observables, magnetic properties are a suitable starting point for a general description of of delocalization and conjugation. In analogy to the anisotropy of the magnetic susceptibility, which is a powerful measure of aromaticity, we investigate the anisotropy of the current (induced) density (ACID). Similar to the square of the wavefunction which defines the total electron density, the ACID scalar field defines the density of delocalized electrons.

The ACID method proved to be an extremely versatile and descriptive tool to analyze delocalization in ground-state molecules, excited states, transition states, organometallics, hyperconjugation and other through-bond and through-space interactions.

Having spent quite some time determining the best way to account for electron delocalization in an upcoming project, the ACID method (“Delocalization of Electrons in Molecules,” J. Phys. Chem. A 2001, 105, 3214-3220. dx.doi.org/10.1021/jp0034426) struck me as near-ideal, leading to my request for the code from the developers.

With the code graciously provided, my first run attempt under Ubuntu 14.04 (and Cygwin, and Ubuntu 10.04 because it was waaaaay back when) went less-than successfully after the successful Gaussian run.

Ubuntu 14.04 make Results (Note: Using gcc 4.8)

user@hostname:~/AICD-2.0.0/$ make

sh modify_shellpath_in_scripts.sh AICD AICD-extract.pl ModifyPov.pl AICD-convert.pl
wd=`pwd` ; \
	sed -i -e "2s+^AICD_BaseDir=.*+AICD_BaseDir=$wd+" AICD
g++ -O2 -ffast-math -march=native -o AICD-smooth_isosurface.o -c AICD-smooth_isosurface.cpp
g++ AICD-smooth_isosurface.o -Wl,-s -o AICD-smooth_isosurface
g++ -O2 -ffast-math -march=native -o AICD-isosurface.o -c AICD-isosurface.cpp
g++ AICD-isosurface.o -Wl,-s -o AICD-isosurface
g++ -O2 -ffast-math -march=native -o AICD-extract.o -c AICD-extract.cpp
g++ AICD-extract.o -Wl,-s -o AICD-extract
g++ -O2 -ffast-math -march=native -o AICD-cube.o -c AICD-cube.cpp
g++ AICD-cube.o -Wl,-s -o AICD-cube
g++ -O2 -ffast-math -march=native -o AICD-remap.o -c AICD-remap.cpp
g++ AICD-remap.o -Wl,-s -o AICD-remap
g++ -O2 -ffast-math -march=native -o AICD-opt_remap.o -c AICD-opt_remap.cpp
g++ AICD-opt_remap.o -Wl,-s -o AICD-opt_remap
g++ -O2 -ffast-math -march=native -o AICD-rotate_mol.o -c AICD-rotate_mol.cpp
g++ AICD-rotate_mol.o -Wl,-s -o AICD-rotate_mol
g++ -O2 -ffast-math -march=native -o AICD-isocut.o -c AICD-isocut.cpp
g++ AICD-isocut.o -Wl,-s -o AICD-isocut
cc -O2 -ffast-math -march=native -o povchem/povchem.o -c povchem/povchem.c
povchem/povchem.c: In function ‘Get_User_Input’:
povchem/povchem.c:1173:10: warning: ignoring return value of ‘scanf’, declared with attribute warn_unused_result [-Wunused-result]
     scanf("%c",&string[i]);
          ^
povchem/povchem.c: In function ‘main’:
povchem/povchem.c:2958:13: warning: ignoring return value of ‘system’, declared with attribute warn_unused_result [-Wunused-result]
       system(tmp);
             ^
povchem/povchem.c:2965:17: warning: ignoring return value of ‘system’, declared with attribute warn_unused_result [-Wunused-result]
           system(tmp);
                 ^
povchem/povchem.c:2985:15: warning: ignoring return value of ‘system’, declared with attribute warn_unused_result [-Wunused-result]
         system(tmp);
               ^
povchem/povchem.c:2989:15: warning: ignoring return value of ‘system’, declared with attribute warn_unused_result [-Wunused-result]
         system(tmp);
               ^
povchem/povchem.c:2994:15: warning: ignoring return value of ‘system’, declared with attribute warn_unused_result [-Wunused-result]
         system(tmp);
               ^
povchem/povchem.c:3003:15: warning: ignoring return value of ‘system’, declared with attribute warn_unused_result [-Wunused-result]
         system(tmp);
               ^
g++ povchem/povchem.o -Wl,-s -o povchem/povchem

Warnings are not errors, so meh. Google searching indicated a code issue (version issue) that someone on Stack Overflow proposed fixing by turning off optimization (removing -O2 from the Makefile). The error goes away.

sh modify_shellpath_in_scripts.sh AICD AICD-extract.pl ModifyPov.pl AICD-convert.pl
wd=`pwd` ; \
	sed -i -e "2s+^AICD_BaseDir=.*+AICD_BaseDir=$wd+" AICD
g++ -O2 -ffast-math -march=native -o AICD-smooth_isosurface.o -c AICD-smooth_isosurface.cpp
g++ AICD-smooth_isosurface.o -Wl,-s -o AICD-smooth_isosurface
g++ -O2 -ffast-math -march=native -o AICD-isosurface.o -c AICD-isosurface.cpp
g++ AICD-isosurface.o -Wl,-s -o AICD-isosurface
g++ -O2 -ffast-math -march=native -o AICD-extract.o -c AICD-extract.cpp
g++ AICD-extract.o -Wl,-s -o AICD-extract
g++ -O2 -ffast-math -march=native -o AICD-cube.o -c AICD-cube.cpp
g++ AICD-cube.o -Wl,-s -o AICD-cube
g++ -O2 -ffast-math -march=native -o AICD-remap.o -c AICD-remap.cpp
g++ AICD-remap.o -Wl,-s -o AICD-remap
g++ -O2 -ffast-math -march=native -o AICD-opt_remap.o -c AICD-opt_remap.cpp
g++ AICD-opt_remap.o -Wl,-s -o AICD-opt_remap
g++ -O2 -ffast-math -march=native -o AICD-rotate_mol.o -c AICD-rotate_mol.cpp
g++ AICD-rotate_mol.o -Wl,-s -o AICD-rotate_mol
g++ -O2 -ffast-math -march=native -o AICD-isocut.o -c AICD-isocut.cpp
g++ AICD-isocut.o -Wl,-s -o AICD-isocut
cc -O2 -ffast-math -march=native -o povchem/povchem.o -c povchem/povchem.c
g++ povchem/povchem.o -Wl,-s -o povchem/povchem

And the tutorial run output itself is below – complete with error.

$ AICD -m 2 -rot 90 0 0 -b 1 0 0 -runpov benzol.nmr.log

Datei benzol.nmr.log wird bearbeitet.
rufe jetzt AICDcube auf

1. Dateienpaar: benzol.nmr.icd40000 und benzol.nmr.pdb

Molek?linformationen:
6 0 0 0 2.6192
6 0 0 2.2677 1.3096
6 0 0 2.2677 -1.3096
6 0 0 0 -2.6192
6 0 0 -2.2677 -1.3096
6 0 0 -2.2677 1.3096
1 0 0 0 4.6525
1 0 0 4.0289 2.3263
1 0 0 4.0289 -2.3263
1 0 0 0 -4.6525
1 0 0 -4.0289 -2.3263
1 0 0 -4.0289 2.3263

40128 = 19 * 44 * 48
Grenzen:
-0.066811 < = Isotropie <= 0.066811
0 <= Anisotropie <= 1.1042
-2.0667 <= X <= 2.0667
-4.7203 <= Y <= 4.7203
-5.1407 <= Z <= 5.1407

Cube-file wird geschrieben. Bitte warten.

1. Input-file: benzol.nmr.icd40000
40128 = 19 * 44 * 48
Grenzen:
-0.0668113 <= Isotropie <= 0.0668113
0 <= Anisotropie <= 1.10423
-2.06672 <= X <= 2.06672
-4.72031 <= Y <= 4.72032
-5.14072 <= Z <= 5.14072

Punkte auf der Isooberfl?che
Z = 0
...................
...................
[…]
...................
...................

Es wurden 2604 Pfeile generiert.
Davon zeigen 1302 in die Isooberfl?che hinein (50%).

Mittelwert der Pfeill?nge: 3.65747
Pfeilstatistik:
0  0  0  0  0  0  0  0  0  0  0  2604  
Dreiecke werden generiert. Bitte warten.
Povray-input wird geschrieben. Bitte warten.
Bildgrenzen:
0 <= Isotropie <= 0
0.05 <= Anisotropie <= 0.05
-1.4586 <= X <= 1.4586
-4.29131 <= Y <= 4.29131
-4.64666 <= Z <= 4.64666

producing benzol.nmr.pov

Welcome to PovChem! This is version 1.00. For general
instructions see http://ludwig.scs.uiuc.edu/~paul/Manual.html

Read configuration file "/home/damianallis/AICD-2.0.0/povray-AICD-templates/povchem.cfg".
Read periodic table "/home/damianallis/AICD-2.0.0/povray-AICD-templates/periodic.tab".

Found 12 atoms...
...loaded into memory.
Wrote 12 spheres...
...and read 24 bonds into memory.

Wrote 0 single, 12 double, 0 triple, 0 higher order, and 0 hydrogen bonds.

benzol.nmr.inc created.
Writing color and atom definitions...
Writing bond definitions...
Pov header benzol.nmr.pov created.


modifying benzol.nmr.pov to benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Molekuel.pov
*** Error in `AICD-rotate_mol': free(): invalid next size (fast): 0x0000000001f22160 ***
/opt/AICD-2.0.0/AICD: line 656: 17519 Aborted                 (core dumped) $Animation $inputbasename.pdb 
$Atom1 $Atom2 $Atom3 < "$RenderMich" > $Pov_RenderMich_pov
Replacing Camera Postion to \<0, 0, -250\>
Inserting #version directive in RenderMich.pov
calling Replace_Pov_Filenames benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Isoober.inc
calling Replace_Pov_Filenames benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Molekuel.inc
calling Replace_Pov_Filenames benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Molekuel.pov
calling Replace_Pov_Filenames benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.RenderMich.pov
calling Replace_Pov_Filenames benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Rotate.inc
calling AICD-isocut -m -1 -M -1 -r 10000000 < benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Isoober.inc.noncut > benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Isoober.inc
This is AICD-isocut....
 arrow_maxlength = -1
 arrow_minlength = -1
 Point of origin:  0 0 0
 Maximum distance: 1e+07
Number of deleted arrows: 0
povray: cannot open the user configuration file /home/damianallis/.povray/3.7/povray.conf: No such file or directory
Persistence of Vision(tm) Ray Tracer Version 3.7.0.unofficial (g++ 4.8 @
x86_64-pc-linux-gnu)
This is an unofficial version compiled by:
Felix Geyer  for Debian 
The POV-Ray Team is not responsible for supporting this version.

POV-Ray is based on DKBTrace 2.12 by David K. Buck & Aaron A. Collins
Copyright 1991-2013 Persistence of Vision Raytracer Pty. Ltd.

Primary POV-Ray 3.7 Architects/Developers: (Alphabetically)
 Chris Cason         Thorsten Froehlich  Christoph Lipka   

...

Support libraries used by POV-Ray:
 ZLib 1.2.8, Copyright 1995-2012 Jean-loup Gailly and Mark Adler
 LibPNG 1.2.50, Copyright 1998-2012 Glenn Randers-Pehrson
 LibJPEG 80, Copyright 1991-2013 Thomas G. Lane, Guido Vollbeding
 LibTIFF 4.0.3, Copyright 1988-1997 Sam Leffler, 1991-1997 SGI
 Boost 1.54, http://www.boost.org/
 OpenEXR, Copyright (c) 2004-2007, Industrial Light & Magic.

Parser Options
 Input file: benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.RenderMich.pov
 Remove bounds........On 
 Split unions.........Off
 Library paths:
   /usr/share/povray-3.7
   /usr/share/povray-3.7/ini
   /usr/share/povray-3.7/include
   /usr/lib/povray3
   /usr/lib/povray3/include
 Clock value:    0.000  (Animation off)
Image Output Options
 Image resolution.....1024 by 768 (rows 1 to 768, columns 1 to 1024).
 Output file..........benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.png, 24 bpp PNG
 Dithering............Off
 Graphic display......Off
 Mosaic preview.......Off
 Continued trace......Off
Information Output Options
 All Streams to console..........On 
 Debug Stream to console.........On 
 Fatal Stream to console.........On 
 Render Stream to console........On 
 Statistics Stream to console....On 
 Warning Stream to console.......On 
==== [Parsing...] ==========================================================
Parse Warning: No objects in scene.
Parse Warning: assumed_gamma not specified, so gamma_correction is turned off
for compatibility with this pre POV-Ray 3.7 scene. See the documentation for
more details.
Parse Warning: This scene did not contain a #version directive. Please be aware
that as of POV-Ray 3.7, unless already specified via an INI option, a #version
is expected as the first declaration in a scene file. POV-Ray may apply
settings to some features that are intended to maintain compatibility with
pre-3.7 scenes. You are strongly encouraged to add a #version statement to the
scene to make your intent clear. Future versions of POV-Ray may make the
presence of a #version statement mandatory.
----------------------------------------------------------------------------
Parser Statistics
----------------------------------------------------------------------------
Finite Objects:            0
Infinite Objects:          0
Light Sources:             0
Total:                     0
----------------------------------------------------------------------------
Parser Time
 Parse Time:       0 hours  0 minutes  0 seconds (0.000 seconds)
             using 1 thread(s) with 0.000 CPU-seconds total
 Bounding Time:    0 hours  0 minutes  0 seconds (0.000 seconds)
             using 1 thread(s) with 0.000 CPU-seconds total
----------------------------------------------------------------------------
Render Options
 Quality:  9
 Bounding boxes.......On   Bounding threshold: 3
 Antialiasing.........Off
==== [Rendering...] ========================================================
Rendered 786432 of 786432 pixels (100%)
----------------------------------------------------------------------------
Render Statistics
Image Resolution 1024 x 768
----------------------------------------------------------------------------
Pixels:           786432   Samples:               0   Smpls/Pxl: 0.00
Rays:             786432   Saved:                 0   Max Level: 1/5
----------------------------------------------------------------------------
Ray->Shape Intersection          Tests       Succeeded  Percentage
----------------------------------------------------------------------------
----------------------------------------------------------------------------
----------------------------------------------------------------------------
----------------------------------------------------------------------------
Render Time:
 Photon Time:      No photons
 Radiosity Time:   No radiosity
 Trace Time:       0 hours  0 minutes  0 seconds (0.031 seconds)
             using 12 thread(s) with 0.143 CPU-seconds total
POV-Ray finished

*******************************************************************************************************
*******************************************************************************************************

Not a particularly Ubuntu-looking problem, and not one for which anything online was available that helped at all. The .pov file is written, but the file is empty. Googling the error above…

*** Error in `AICD-rotate_mol’: free(): invalid next size (fast): 0x0000000001f22160 ***
/opt/AICD-2.0.0/AICD: line 656: 17519 Aborted … (core dumped) $Animation $inputbasename.pdb

… sends you to, among other pages, one Stack Overflow discussion that starts with the error and ends with recommendations of finding a good debugger – not something I’d any interest in pursuing as part of getting this code to work.

I did not find any version-dependence discussion with the errors above, but opted to try using an older gcc – that I wasn’t about to install on my new Ubuntu box. Ergo, find an older distro that installs an older gcc. I threw it on an OpenSUSE 11.2 image because the developer group is in Germany – that’s basically why – and because I couldn’t find an old enough Ubuntu distro that would update/upgrade/install apps from the current repositories.

OpenSUSE 11.2 make Results (Note: Using gcc 4.5)

From a fresh OpenSUSE 11.2 install, I added the appropriate developer tools to get the ACID code made (note the specific request for gcc 4.5):

sudo zypper install autoconf automake binutils-gold gcc gcc45 glibc-devel kernel-desktop-devel kernel-devel kernel-source linux-glibc-devel m4 make gcc-c++

Then added povray (and nano, but that’s because I hate vi)

sudo zypper install povray nano

The build process then goes as follows – with no issues.

sh modify_shellpath_in_scripts.sh AICD AICD-extract.pl ModifyPov.pl AICD-convert.pl
wd=`pwd` ; \
	sed -i -e "2s+^AICD_BaseDir=.*+AICD_BaseDir=$wd+" AICD
g++ -O2 -ffast-math -march=native -o AICD-smooth_isosurface.o -c AICD-smooth_isosurface.cpp
g++ AICD-smooth_isosurface.o -Wl,-s -o AICD-smooth_isosurface
g++ -O2 -ffast-math -march=native -o AICD-isosurface.o -c AICD-isosurface.cpp
g++ AICD-isosurface.o -Wl,-s -o AICD-isosurface
g++ -O2 -ffast-math -march=native -o AICD-extract.o -c AICD-extract.cpp
g++ AICD-extract.o -Wl,-s -o AICD-extract
g++ -O2 -ffast-math -march=native -o AICD-cube.o -c AICD-cube.cpp
g++ AICD-cube.o -Wl,-s -o AICD-cube
g++ -O2 -ffast-math -march=native -o AICD-remap.o -c AICD-remap.cpp
g++ AICD-remap.o -Wl,-s -o AICD-remap
g++ -O2 -ffast-math -march=native -o AICD-opt_remap.o -c AICD-opt_remap.cpp
g++ AICD-opt_remap.o -Wl,-s -o AICD-opt_remap
g++ -O2 -ffast-math -march=native -o AICD-rotate_mol.o -c AICD-rotate_mol.cpp
g++ AICD-rotate_mol.o -Wl,-s -o AICD-rotate_mol
g++ -O2 -ffast-math -march=native -o AICD-isocut.o -c AICD-isocut.cpp
g++ AICD-isocut.o -Wl,-s -o AICD-isocut
cc -O2 -ffast-math -march=native -o povchem/povchem.o -c povchem/povchem.c
g++ povchem/povchem.o -Wl,-s -o povchem/povchem

And then there’s the run through the tutorial again…

Datei benzol.nmr.log wird bearbeitet.
Weder benzol.nmr.icd40000 noch benzol.nmr.icd40000.gz existiert.
Weder benzol.nmr.icd noch benzol.nmr.icd.gz existieren.
Starting AICDextractsh
Separate AICD-output in file: test.txt
Extrahiere_Gitterpunkte

1. Input-file: benzol.nmr.icd-raw
Es wurden 84726 Stromdichtetensoren eingelesen.

Grenzen:
-0.0668113 < = Isotropie <= 0.0668113
0 <= Anisotropie <= 1.10423
-20.3137 <= X <= 20.3137
-23.6411 <= Y <= 23.6411
-24.2643 <= Z <= 24.2643

Atomposition:
0.0214824  
0  
2.61948  

Radius: 26.8838
Datei benzol.nmr wird bearbeitet.
Weder benzol.nmr.icd40000 noch benzol.nmr.icd40000.gz existiert.

Input-file: benzol.nmr.icd

Output-file: benzol.nmr.icd40000
========================
 1. Iteration
========================
AICD-remap 40000 1 - benzol.nmr.icd40000 benzol.nmr.icd 

Input-file: benzol.nmr.icd
Anzahl der Stromdichtetensoren: 84726
Grenzen:
-0.0668113 <= Isotropie <= 0.0668113
0 <= Anisotropie <= 1.10423
-10.3137 <= X <= 10.3137
-13.6411 <= Y <= 13.6411
-14.2643 <= Z <= 14.2643


Grenzen des kartesischen Netzes:
-0.0668113 <= Isotropie <= 0.0668113
0 <= Anisotropie <= 1.10423
-10.3137 <= X <= 10.3137
-13.6411 <= Y <= 13.6411
-14.2643 <= Z <= 14.2643

Globales Netz: 40404 = 28 * 37 * 39

Input-file: benzol.nmr.icd

1 %
2 %
3 %

...

Z = 47
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................
...................

Es wurden 2604 Pfeile generiert.
Davon zeigen 1302 in die Isooberfl‰che hinein (50%).

Mittelwert der Pfeill‰nge: 0.0642683
Pfeilstatistik:
0  32  80  52  124  268  476  716  412  260  144  40  
Dreiecke werden generiert. Bitte warten.
Povray-input wird geschrieben. Bitte warten.
Bildgrenzen:
0 <= Isotropie <= 0
0.05 <= Anisotropie <= 0.05
-1.4586 <= X <= 1.4586
-4.29131 <= Y <= 4.29131
-4.64666 <= Z <= 4.64666

producing benzol.nmr.pov

Welcome to PovChem! This is version 1.00. For general
instructions see http://ludwig.scs.uiuc.edu/~paul/Manual.html

Read configuration file "/home/damianallis/ACID/AICD-2.0.0/povray-AICD-templates/povchem.cfg".
Read periodic table "/home/damianallis/ACID/AICD-2.0.0/povray-AICD-templates/periodic.tab".

Found 12 atoms...
 ...loaded into memory.
Wrote 12 spheres...
 ...and read 24 bonds into memory.

Wrote 0 single, 12 double, 0 triple, 0 higher order, and 0 hydrogen bonds.

benzol.nmr.inc created.
Writing color and atom definitions...
Writing bond definitions...
Pov header benzol.nmr.pov created.


modifying benzol.nmr.pov to benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Molekuel.pov
Replacing Camera Postion to \<0, 0, -250\>
Inserting #version directive in RenderMich.pov
calling Replace_Pov_Filenames benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Isoober.inc
calling Replace_Pov_Filenames benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Molekuel.inc
calling Replace_Pov_Filenames benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Molekuel.pov
calling Replace_Pov_Filenames benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.RenderMich.pov
calling Replace_Pov_Filenames benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Rotate.inc
calling AICD-isocut -m -1 -M -1 -r 10000000 < benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Isoober.inc.noncut > benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Isoober.inc
This is AICD-isocut....
  arrow_maxlength = -1
  arrow_minlength = -1
  Point of origin:  0 0 0
  Maximum distance: 1e+07
Number of deleted arrows: 0
povray: cannot open the user configuration file /home/damianallis/.povray/3.6/povray.conf: No such file or directory
Persistence of Vision(tm) Ray Tracer Version 3.6.1 (g++ @ x86_64-unknown-linux-g
nu)
This is an unofficial version compiled by:
 SUSE LINUX Products GmbH, Nuernberg, Germany
 The POV-Ray Team(tm) is not responsible for supporting this version.
POV-Ray is based on DKBTrace 2.12 by David K. Buck & Aaron A. Collins
Copyright 1991-2003 Persistence of Vision Team
Copyright 2003-2004 Persistence of Vision Raytracer Pty. Ltd.

...

Support libraries used by POV-Ray:
  ZLib 1.2.5, Copyright 1995-1998 Jean-loup Gailly and Mark Adler
  LibPNG 1.4.4, Copyright 1998-2002 Glenn Randers-Pehrson
  LibJPEG 6b, Copyright 1998 Thomas G. Lane
  LibTIFF 3.9.4, Copyright 1988-1997 Sam Leffler, 1991-1997 SGI
Redirecting Options
  All Streams to console..........On 
  Debug Stream to console.........On 
  Fatal Stream to console.........On 
  Render Stream to console........On 
  Statistics Stream to console....On 
  Warning Stream to console.......On 
Parsing Options
  Input file: benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.RenderMich.pov (compatible
 to version 3.61)
  Remove bounds........On 
  Split unions.........Off
  Library paths:
    /usr/share/povray-3.6
    /usr/share/povray-3.6/ini
    /usr/share/povray-3.6/include
    /usr/lib/povray3
    /usr/lib/povray3/include
Output Options
  Image resolution 1024 by 768 (rows 1 to 768, columns 1 to 1024).
  Output file: /home/damianallis/ACID/AICD-2.0.0/tutorial-data/benzol.nmr_40000_
0.050_1_0_0_Aniso_4.2.png, 24 bpp PNG
  Graphic display......Off
  Mosaic preview.......Off
  CPU usage histogram..Off
  Continued trace......Off
Tracing Options
  Quality:  9
  Bounding boxes.......On   Bounding threshold: 3
  Light Buffer.........On 
  Vista Buffer.........On   Draw Vista Buffer....Off
  Antialiasing.........Off
  Clock value:    0.000  (Animation off)

  0:00:00 Parsing
File: /usr/share/povray-3.6/include/glass_old.inc  Line: 18
Parse Warning:  Due to changes in version 3.1, you must add interior {I_Glass}
 to all objects calling glass_old.inc textures and finishes... 
File: benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Molekuel.pov  Line: 1
Possible Parse Error: All #version and #declares of float, vector, and color
 require semi-colon ';' at end.
File: benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Molekuel.pov  Line: 2
Possible Parse Error: All #version and #declares of float, vector, and color
 require semi-colon ';' at end.
File: benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Molekuel.pov  Line: 21
Possible Parse Error: All #version and #declares of float, vector, and color
 require semi-colon ';' at end.
File: benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Isoober.inc  Line: 9
Possible Parse Error: Text may not be displayed as expected. Please refer to the
 user manual regarding changes in POV-Ray 3.5 and later.
File: benzol.nmr_40000_0.050_1_0_0_Aniso_4.2.Isoober.inc  Line: 17
Possible Parse Error: Text may not be displayed as expected. Please refer to the
 user manual regarding changes in POV-Ray 3.5 and later.

  0:00:01 Creating bounding slabs
  0:00:01 Creating vista buffer
  0:00:01 Creating light buffers
  0:00:01 Creating light buffers
  0:00:01 Creating light buffers 508K tokens
Scene Statistics
  Finite objects:         5233
  Infinite objects:          0
  Light sources:             4
  Total:                  5237

  0:00:00 Rendering line 1 of 768
  0:00:01 Rendering line 237 of 768
  0:00:02 Rendering line 292 of 768
  0:00:03 Rendering line 383 of 768
  0:00:04 Rendering line 479 of 768
  0:00:05 Rendering line 538 of 768
  0:00:05 Done Tracing
Render Statistics
Image Resolution 1024 x 768
----------------------------------------------------------------------------
Pixels:           786432   Samples:          786432   Smpls/Pxl: 1.00
Rays:             786432   Saved:                 0   Max Level: 1/5
----------------------------------------------------------------------------
Ray->Shape Intersection          Tests       Succeeded  Percentage
----------------------------------------------------------------------------
Cone/Cylinder                  1927296          164391      8.53
Mesh                            862314          164798     19.11
Sphere                           21764            5562     25.56
Bounding Box                  39215172        13715953     34.98
Light Buffer                  60713588        16726656     27.55
Vista Buffer                  11930298         6504493     54.52
----------------------------------------------------------------------------
Calls to Noise:                   0   Calls to DNoise:              10
----------------------------------------------------------------------------
Shadow Ray Tests:           1215156   Succeeded:                109903
----------------------------------------------------------------------------
Smallest Alloc:                   9 bytes
Largest  Alloc:             2097160 bytes
Peak memory used:          24172005 bytes
Total Scene Processing Times
  Parse Time:    0 hours  0 minutes  1 seconds (1 seconds)
  Photon Time:   0 hours  0 minutes  0 seconds (0 seconds)
  Render Time:   0 hours  0 minutes  5 seconds (5 seconds)
  Total Time:    0 hours  0 minutes  6 seconds (6 seconds)

Producing the beautiful and informative image shown above (despite the Parse Warning and Possible Parse Errors above).

For the record, compiler version differences can be subtle beasts – and my whole discussion above is not directly related to getting ACID to work. The above is a mix of brute force and resource consumption that provided a much faster turnaround in positive results than any amount of work spent in Ubuntu itself to make ACID run successfully – and does not require modifying anyone’s code (an added bonus for those who might have to ask questions of the developers).

So, I would posit that the issue turns out to be a 4.5 vs. 4.8 thing, but I’ve not explored the issue further (frankly, no time – and it’s working – so I don’t know at what point in gcc development that the make process began to produce errors in the ACID code). This isn’t technically the way to run it natively in Ubuntu, but is a way to get it to run in Ubuntu, which still counts in my book. In theory, a VirtualBox install and some book work on finding the dates of certain distros provides yet another way of getting stuff to just work.

That said, getting OpenSUSE to play with the Networking in VirtualBox appears to be a real pain. My solution for that for getting files back-and-forth to process?

DropBox.

Generating Molecular Orbitals (And Visualizing Assorted Properties) With The Gaussian09 cubegen Utility

To begin, this post owes its existence to the efforts of Dr. Douglas Fox at Gaussian, Inc., who provided me with an alternative explanation of how the cubegen utility works. After much wailing and gnashing of teeth, I intend on taking Dr. Fox’s advice and asking Gaussian Support for assistance earlier in my endeavors. What follows below, I hope, will save you some significant frustration (and, given how little there is online that really describes the extra workings of cubegen in a clear and example’ed way, it is my expectation that this page appeared early in your search list).

What I wanted out of cubegen that I couldn’t figure out how to get:

The situation was simple. I wanted my molecule centered and bound within an arbitrarily-sized box (X,Z,Y) for making images and doing additional post-processing. Specifically, I wanted to be able to take many different molecules (from hydrogen gas to big biomolecules) defined within the same-sized box for layering and presentation (different boxes for each, but all the same size).

I am assuming for this that you’re using cubegen from a terminal (not within GaussView or the like) to produce .cub/.cube files for use in some kind of rendering-capable program (like VESTA or VMD) and that cubegen and formchk are in your PATH (either properly placed or by running the Gaussian install script). I’ll be demonstrating usage with benzene (C6H6) and the benzene cation (C6H6+).

1. The Checkpoint File

To extract any kind of data for making .cub/.cube files, you need a checkpoint file (.chk) from your run. This is performed by adding a %chk=FILENAME.chk line to the top of the input file (which, if you’re a Gaussian user, you likely already know). If you want additional properties cube’d, check the Gaussian Tech Document, specifically looking at the Pop keyword for most of the properties you’d want visualized (this data gets placed into the .chk file for .cub/.cube generation after the run). For the standard molecular orbitals, they’re already saved in the .chk file (or their coefficients, anyway).

For benzene.gjf:

%chk=benzene.chk
# b3lyp/6-31G(d,p)

Benzene

0 1
 C                  1.20809735    0.69749533   -0.00000000
 C                  0.00000000    1.39499067   -0.00000000
 C                 -1.20809735    0.69749533   -0.00000000
 C                 -1.20809735   -0.69749533   -0.00000000
 C                  0.00000000   -1.39499067   -0.00000000
 C                  1.20809735   -0.69749533   -0.00000000
 H                  2.16038781    1.24730049   -0.00000000
 H                  0.00000000    2.49460097   -0.00000000
 H                 -2.16038781    1.24730049   -0.00000000
 H                 -2.16038781   -1.24730049   -0.00000000
 H                  0.00000000   -2.49460097   -0.00000000
 H                  2.16038781   -1.24730049   -0.00000000

For benzenecation.gjf:

%chk=benzenecation.chk
# b3lyp/6-31G(d,p)

Benzene cation

1 2
 C                  1.20809735    0.69749533   -0.00000000
 C                  0.00000000    1.39499067   -0.00000000
 C                 -1.20809735    0.69749533   -0.00000000
 C                 -1.20809735   -0.69749533   -0.00000000
 C                  0.00000000   -1.39499067   -0.00000000
 C                  1.20809735   -0.69749533   -0.00000000
 H                  2.16038781    1.24730049   -0.00000000
 H                  0.00000000    2.49460097   -0.00000000
 H                 -2.16038781    1.24730049   -0.00000000
 H                 -2.16038781   -1.24730049   -0.00000000
 H                  0.00000000   -2.49460097   -0.00000000
 H                  2.16038781   -1.24730049   -0.00000000

2. Convert The .chk To .fchk With formchk

As per the Gaussian Tech Doc:

formchk converts the data in a Gaussian checkpoint file into a formatted form which is suitable for input into a variety of visualization software.

Basically, making the .chk file something that cubegen can manipulate to generate .cub/.cube files of orbitals, densities, electrostatic potentials, etc. This run is simple for most users (for the rest, see formchk).

formchk benzene.chk benzene.fchk
formchk benzenecation.chk benzenecation.fchk

3. Using cubegen

And now the fun begins. A typical cubegen run looks like the following:

cubegen 0 MO=HOMO benzene.fchk benzene_HOMO.cub 0 h

cubegen – run cubegen
0 – an old memory flag (must be there, but not important)
MO=HOMO – generate the highest occupied molecular orbital
benzene.fchk – the .fchk file
benzene_HOMO.cub – the generated .cub file
0 – use the default grid point specification (80*80*80 points total in the whole cube file)
h – write out the .cub file with headers

The output you find summarized in VESTA is below for this case.

DEFAULT:
OpenGL version: 2.1 INTEL-8.26.34
Video configuration: Intel HD Graphics 4000 OpenGL Engine
Maximum supported width and height of the viewport: 16384 x 16384
OpenGL depth buffer bit: 16

/Users/damianallis/benzene_HOMO_default_0.cub
====================================================================================
Title Benzene MO=HOMO
Dimensions 87 91 65

Lattice parameters

a b c alpha beta gamma
9.39704 9.82909 7.02078 90.0000 90.0000 90.0000

Unit-cell volume = 648.469273 Ã…^3

Total number of polygons and unique vertices on slices;
(1 0 0): 0 ( 0), 0 ( 0)
(0 1 0): 0 ( 0), 0 ( 0)
(0 0 1): 0 ( 0), 0 ( 0)
====================================================================================

====================================================================================
Title Benzene

Lattice type P
Space group name P 1
Space group number 1
Setting number 1

Lattice parameters

a b c alpha beta gamma
1.00000 1.00000 1.00000 90.0000 90.0000 90.0000

Unit-cell volume = 1.000000 Ã…^3

Structure parameters

x y z Occ. B Site Sym.
1 C C1 4.65450 6.23638 3.44640 1.000 1.000 1 –
2 C C2 5.86259 5.53889 3.44640 1.000 1.000 1 –
3 C C3 5.86259 4.14389 3.44640 1.000 1.000 1 –
4 C C4 4.65450 3.44640 3.44640 1.000 1.000 1 –
5 C C5 3.44640 4.14389 3.44640 1.000 1.000 1 –
6 C C6 3.44640 5.53889 3.44640 1.000 1.000 1 –
7 H H1 4.65450 7.33599 3.44640 1.000 1.000 1 –
8 H H2 6.81488 6.08869 3.44640 1.000 1.000 1 –
9 H H3 6.81488 3.59409 3.44640 1.000 1.000 1 –
10 H H4 4.65450 2.34679 3.44640 1.000 1.000 1 –
11 H H5 2.49411 3.59409 3.44640 1.000 1.000 1 –
12 H H6 2.49411 6.08869 3.44640 1.000 1.000 1 –
====================================================================================

Number of polygons and unique vertices on isosurface = 16904 (8460)
12 atoms, 12 bonds, 0 polyhedra; CPU time = 39 ms

For the coarse grid (-2) case:

cubegen 0 MO=HOMO benzene.fchk benzene_HOMO_default_m2.cub -2 h

The output you find summarized in VESTA is below for this case.

/Users/damianallis/benzene_HOMO_default_m2.cub
====================================================================================
Title Benzene MO=HOMO
Dimensions 54 56 40

Lattice parameters

a b c alpha beta gamma
9.52518 9.87796 7.05569 90.0000 90.0000 90.0000

Unit-cell volume = 663.865482 Ã…^3

Total number of polygons and unique vertices on slices;
(1 0 0): 0 ( 0), 0 ( 0)
(0 1 0): 0 ( 0), 0 ( 0)
(0 0 1): 0 ( 0), 0 ( 0)
====================================================================================

====================================================================================
Title Benzene

Lattice type P
Space group name P 1
Space group number 1
Setting number 1

Lattice parameters

a b c alpha beta gamma
1.00000 1.00000 1.00000 90.0000 90.0000 90.0000

Unit-cell volume = 1.000000 Ã…^3

Structure parameters

x y z Occ. B Site Sym.
1 C C1 4.65450 6.23638 3.44640 1.000 1.000 1 –
2 C C2 5.86259 5.53889 3.44640 1.000 1.000 1 –
3 C C3 5.86259 4.14389 3.44640 1.000 1.000 1 –
4 C C4 4.65450 3.44640 3.44640 1.000 1.000 1 –
5 C C5 3.44640 4.14389 3.44640 1.000 1.000 1 –
6 C C6 3.44640 5.53889 3.44640 1.000 1.000 1 –
7 H H1 4.65450 7.33599 3.44640 1.000 1.000 1 –
8 H H2 6.81488 6.08869 3.44640 1.000 1.000 1 –
9 H H3 6.81488 3.59409 3.44640 1.000 1.000 1 –
10 H H4 4.65450 2.34679 3.44640 1.000 1.000 1 –
11 H H5 2.49411 3.59409 3.44640 1.000 1.000 1 –
12 H H6 2.49411 6.08869 3.44640 1.000 1.000 1 –
====================================================================================

Number of polygons and unique vertices on isosurface = 6516 (3266)
12 atoms, 12 bonds, 0 polyhedra; CPU time = 10 ms

For the medium grid (-3) case:

cubegen 0 MO=HOMO benzene.fchk benzene_HOMO_default_m3.cub -3 h

The output you find summarized in VESTA is below for this case.

/Users/damianallis/benzene_HOMO_default_m3.cub
====================================================================================
Title Benzene MO=HOMO
Dimensions 107 111 79

Lattice parameters

a b c alpha beta gamma
9.43701 9.78980 6.96751 90.0000 90.0000 90.0000

Unit-cell volume = 643.703858 Ã…^3

Total number of polygons and unique vertices on slices;
(1 0 0): 0 ( 0), 0 ( 0)
(0 1 0): 0 ( 0), 0 ( 0)
(0 0 1): 0 ( 0), 0 ( 0)
====================================================================================

====================================================================================
Title Benzene

Lattice type P
Space group name P 1
Space group number 1
Setting number 1

Lattice parameters

a b c alpha beta gamma
1.00000 1.00000 1.00000 90.0000 90.0000 90.0000

Unit-cell volume = 1.000000 Ã…^3

Structure parameters

x y z Occ. B Site Sym.
1 C C1 4.65450 6.23638 3.44640 1.000 1.000 1 –
2 C C2 5.86259 5.53889 3.44640 1.000 1.000 1 –
3 C C3 5.86259 4.14389 3.44640 1.000 1.000 1 –
4 C C4 4.65450 3.44640 3.44640 1.000 1.000 1 –
5 C C5 3.44640 4.14389 3.44640 1.000 1.000 1 –
6 C C6 3.44640 5.53889 3.44640 1.000 1.000 1 –
7 H H1 4.65450 7.33599 3.44640 1.000 1.000 1 –
8 H H2 6.81488 6.08869 3.44640 1.000 1.000 1 –
9 H H3 6.81488 3.59409 3.44640 1.000 1.000 1 –
10 H H4 4.65450 2.34679 3.44640 1.000 1.000 1 –
11 H H5 2.49411 3.59409 3.44640 1.000 1.000 1 –
12 H H6 2.49411 6.08869 3.44640 1.000 1.000 1 –
====================================================================================

Number of polygons and unique vertices on isosurface = 25532 (12774)
12 atoms, 12 bonds, 0 polyhedra; CPU time = 51 ms

For the fine grid (-4) case:

cubegen 0 MO=HOMO benzene.fchk benzene_HOMO_default_m4.cub -4 h

The output you find summarized in VESTA is below for this case.

/Users/damianallis/benzene_HOMO_default_m4.cub
====================================================================================
Title Benzene MO=HOMO
Dimensions 212 221 157

Lattice parameters

a b c alpha beta gamma
9.34876 9.74564 6.92337 90.0000 90.0000 90.0000

Unit-cell volume = 630.786281 Ã…^3

Total number of polygons and unique vertices on slices;
(1 0 0): 0 ( 0), 0 ( 0)
(0 1 0): 0 ( 0), 0 ( 0)
(0 0 1): 0 ( 0), 0 ( 0)
====================================================================================

====================================================================================
Title Benzene

Lattice type P
Space group name P 1
Space group number 1
Setting number 1

Lattice parameters

a b c alpha beta gamma
1.00000 1.00000 1.00000 90.0000 90.0000 90.0000

Unit-cell volume = 1.000000 Ã…^3

Structure parameters

x y z Occ. B Site Sym.
1 C C1 4.65450 6.23638 3.44640 1.000 1.000 1 –
2 C C2 5.86259 5.53889 3.44640 1.000 1.000 1 –
3 C C3 5.86259 4.14389 3.44640 1.000 1.000 1 –
4 C C4 4.65450 3.44640 3.44640 1.000 1.000 1 –
5 C C5 3.44640 4.14389 3.44640 1.000 1.000 1 –
6 C C6 3.44640 5.53889 3.44640 1.000 1.000 1 –
7 H H1 4.65450 7.33599 3.44640 1.000 1.000 1 –
8 H H2 6.81488 6.08869 3.44640 1.000 1.000 1 –
9 H H3 6.81488 3.59409 3.44640 1.000 1.000 1 –
10 H H4 4.65450 2.34679 3.44640 1.000 1.000 1 –
11 H H5 2.49411 3.59409 3.44640 1.000 1.000 1 –
12 H H6 2.49411 6.08869 3.44640 1.000 1.000 1 –
====================================================================================

Number of polygons and unique vertices on isosurface = 100680 (50348)
12 atoms, 12 bonds, 0 polyhedra; CPU time = 155 ms

These all generate a file containing the highest occupied molecular orbital (or one of the degenerate HOMO’s in this case. Do I have to qualify that this doesn’t mean what 99.5% of the people coming to this page thinks this means?). The box is generated by something in cubegen to be 9.3ish x 9.7ish x 6.9ish Angstroms on a side and containing X points per Angstrom (and you can change the fineness of the grid points). The image below shows the four cases for the benzene HOMO. Click to see larger versions if you want to see the influence of grid fineness on the final image.

benzene_homo_gaussian_defaults_small

Click for a larger view.

Now, then, while the boxes are almost all identical, the same molecule and input gives four slightly different results. Fine for individual images, but not ideal for the obsessive-compulsive image maker. Also, you see how a box simply bounds the molecule, meaning no standardization of size if you needed that standardization for some reason.

   a        b        c       alpha    beta     gamma
 9.39704  9.82909  7.02078  90.0000  90.0000  90.0000 < - default (0)
 9.52518  9.87796  7.05569  90.0000  90.0000  90.0000 <- coarse (-2)
 9.43701  9.78980  6.96751  90.0000  90.0000  90.0000 <- medium (-3)
 9.34876  9.74564  6.92337  90.0000  90.0000  90.0000 <- fine (-4)

So, for a specific case - suppose I wanted this orbital in a box exactly 15 x 20 x 25 Angstroms on a side with the molecule offset from the center by -1.0 Angstrom in each direction.

I was pleased to finally discover that cubegen allows for that, although you have to ask Gaussian Support to find out how (until now, that is) and you need to do a little bit of math to get the placement right (or use the excel file I've linked in a .zip file found at 2014june7_cubegen_excel_file.xlsx).

You begin with the following:

cubegen 0 MO=HOMO benzene.fchk benzene_HOMO.cub -1 h

But for -1, Where do the numbers go?

From the Gaussian Tech Doc:

A value of -1 says to read the cube specification from the input stream, according to the following format:

IFlag, X0, Y0, Z0 Output unit number and initial point.
N1, X1, Y1, Z1 Number of points and step-size in the X-direction.
N2, X2, Y2, Z2 Number of points and step-size in the Y-direction.
N3, X3, Y3, Z3 Number of points and step-size in the Z-direction.

IFlag is the output unit number. If IFlag is less than 0, then a formatted file will be produced; otherwise, an unformatted file will be written.

Admittedly, “input stream” made no sense to me upon first and second read. I just knew that the program didn’t do anything when I ran it. Now obvious, this means you input the cube specifications by typing in (or, better, pasting in) the 16 numbers it asks for.

Continuing…

The -1 tells cubegen to “expect more input.” In this case, without explanation first, my input would look as below:

-12  -6.50000  -9.00000 -11.50000
 60   0.25000   0.00000   0.00000
 80   0.00000   0.25000   0.00000
100   0.00000   0.00000   0.25000

Which you just paste into your terminal at the new line (having pressed ENTER after typing out the cubegen line above).

How this works (and note the use of minus signs!):

-[# Atoms] -[Start Point For Box In X] -[Start Point For Box In Y] -[Start Point For Box In Z]
[Number of Points In X]   [Grid Fineness In X]   [Grid Fineness In Y]   [Grid Fineness In Z]
[Number of Points In Y]   [Grid Fineness In X]   [Grid Fineness In Y]   [Grid Fineness In Z]
[Number of Points In Y]   [Grid Fineness In X]   [Grid Fineness In Y]   [Grid Fineness In Z]

Assuming orthogonality in your box, the off-diagonals for the grid fineness matrix are zero.

-[# Atoms] -[Start Point For Box In X] -[Start Point For Box In Y] -[Start Point For Box In Z]
[Number of Points In X]  [Grid Fineness In X]   0.000   0.000
[Number of Points In Y]  0.000   [Grid Fineness In Y]   0.000
[Number of Points In Y]  0.000   0.000   [Grid Fineness In Z]

4. -6.5, -9, -11.5?

You build the box around your molecule in cubegen, which means you combine (1) where you want the molecule positions with (2) the number of grid points along each direction and (3) the fineness of the grid to generate the box. Here, I’m starting my hypothetical box at -6.5 in X, -9 in Y, and -11.5 in Z, then building out my molecule 121*.25 points in X, 161*.25 in Y, 201*.25 in Z. This will produce the intended box size with the molecule technically centered at the origin in the box (0,0,0), but the generation of all 121, 161, and 201 points in X, Y, and Z will result in the box going from -6.5 to 8.5, -9 to 11, and -11.5 to 13.5 (and there’s your asymmetry in the box). Alternatively, you could think of it as generating a box 15 x 20 x 25, then placing the center of the molecule at 6.5, 9, 11.5 (but you don’t specify the box size directly, instead relying on the relative position of the molecule and the fineness of the grid to determine the position (from which you could work back to get the number of points you needed in each direction if you knew the size of the box you wanted. Yes, you might have to re-read that a few times).

I demonstrate this below for a benzene orbital “walk” along X using direct output from VESTA. The rest of the numbers in my matrix above are the same except for the “-[Start Point For The Box In X]” value.

benzene_homo_walk

The benzene walk (numbers show the spacing based on the cubegen input above).

5. Formula For Boxes And Grid Points

You can, in fact, work from the box size you want and relative position of the molecule in that box with some simple math. That looks like the table below:

-(# Atoms)           -(X Position)  -(Y Position)  -(Z Position)
(Box Size / X Mesh)    X Mesh         0.00000        0.00000
(Box Size / Y Mesh)    0.00000        Y Mesh         0.00000
(Box Size / Z Mesh)    0.00000        0.00000        Z Mesh

You specify # Atoms, X Position, Y Position, Z Position, X Mesh, Y Mesh and Z Mesh, then decide on how big your box is going to be. Also, note that X Position, Y Position, Z Position all need to be 1/2 the size of your box if you want the molecule centered. A way to help force this is to force the molecule to have its center of mass shifted to the origin using Symm=COM in your input file.

As mentioned above, a simple excel file for performing this task is provided for download at 2014june7_cubegen_excel_file.xlsx.

6. Lastly, A Procedure For Scripting The Generation Of Many Orbitals

That first stone passed, everything about making custom .cub/.cube files finally made sense. But it lead to another problem. Suppose I want to generate many molecular orbitals. Does one have to paste in the IFlag…Z3 block each time?

Thankfully, this process can be scripted to automation as well, although it’s not just a matter of pasting IFlag…Z3 below each run of cubegen. Doing that produces the following…

Example:

This isn’t a cubegen problem, it’s a Linux issue with the interpretation of stdin. The cubegen script needs to be fed in the matrix in a file (say cubegen.dat if you always want the same .cub/.cube file generated) or via the use of an EOF call.

Cubegen.dat:

cubegen 0 MO=1 benzene.fchk benzene_MO1.cub -1 h < cubegen.dat
cubegen 0 MO=2 benzene.fchk benzene_MO2.cub -1 h < cubegen.dat
cubegen 0 MO=3 benzene.fchk benzene_MO3.cub -1 h < cubegen.dat
...

EOF

cubegen 0 MO=1 benzene.fchk benzene_MO1.cub -1 h << EOF
-12  -6.50000  -9.00000 -11.50000
 60   0.25000   0.00000   0.00000
 80   0.00000   0.25000   0.00000
100   0.00000   0.00000   0.25000
EOF
cubegen 0 MO=2 benzene.fchk benzene_MO2.cub -1 h << EOF
-12  -6.50000  -9.00000 -11.50000
 60   0.25000   0.00000   0.00000
 80   0.00000   0.25000   0.00000
100   0.00000   0.00000   0.25000
EOF
cubegen 0 MO=3 benzene.fchk benzene_MO3.cub -1 h << EOF
-12  -6.50000  -9.00000 -11.50000
 60   0.25000   0.00000   0.00000
 80   0.00000   0.25000   0.00000
100   0.00000   0.00000   0.25000
EOF
...

7. What’s The Deal With The Benzene Cation?

Nothing, except I saw a question in my perusing of cubegen problems and found one related to UHF wavefunctions. How do you render alpha spin orbitals and beta spin orbitals? The answer is you dig into the .log file for the orbital energies and count (to the best of my knowledge).

Benzene (21 alpha/beta-occupied)

 Alpha  occ. eigenvalues -- -10.18955 -10.18928 -10.18928 -10.18872 -10.18872
 Alpha  occ. eigenvalues -- -10.18845  -0.84761  -0.73971  -0.73971  -0.59595
 Alpha  occ. eigenvalues --  -0.59595  -0.51588  -0.45423  -0.43943  -0.41518
 Alpha  occ. eigenvalues --  -0.41518  -0.36090  -0.33862  -0.33862  -0.24750
 Alpha  occ. eigenvalues --  -0.24750
 Alpha virt. eigenvalues --   0.00266   0.00266   0.08636   0.14126   0.14126
 Alpha virt. eigenvalues --   0.16238   0.17957   0.17957   0.18681   0.29989
 Alpha virt. eigenvalues --   0.29989   0.31908   0.31908   0.46637   0.52628
 Alpha virt. eigenvalues --   0.54782   0.55099   0.56222   0.59294   0.60077
 Alpha virt. eigenvalues --   0.60077   0.60084   0.60084   0.62384   0.62384
 Alpha virt. eigenvalues --   0.66653   0.66653   0.74180   0.81178   0.81178
 Alpha virt. eigenvalues --   0.82134   0.83694   0.83694   0.91676   0.93745
 Alpha virt. eigenvalues --   0.93745   0.95812   1.08054   1.08054   1.12992
 Alpha virt. eigenvalues --   1.12992   1.20098   1.26111   1.30051   1.40786
 Alpha virt. eigenvalues --   1.40786   1.42585   1.42585   1.42914   1.42914
 Alpha virt. eigenvalues --   1.74102   1.76078   1.80542   1.87583   1.90680
 Alpha virt. eigenvalues --   1.90680   1.97195   1.97195   1.97924   1.97924
 Alpha virt. eigenvalues --   2.02762   2.07664   2.07664   2.29609   2.29609
 Alpha virt. eigenvalues --   2.34429   2.34429   2.35491   2.39944   2.40328
 Alpha virt. eigenvalues --   2.40328   2.44636   2.44636   2.48731   2.48731
 Alpha virt. eigenvalues --   2.50802   2.58538   2.58538   2.60300   2.65987
 Alpha virt. eigenvalues --   2.75521   2.80103   2.80103   3.03123   3.03123
 Alpha virt. eigenvalues --   3.18490   3.20485   3.21867   3.21867   3.37166
 Alpha virt. eigenvalues --   3.48298   3.48298   3.93339   4.13215   4.16289
 Alpha virt. eigenvalues --   4.16289   4.43754   4.43754   4.82384

RHF wave functions are easy as the alpha and beta spin orbitals are identical (so you just call one).

Benzene Cation (21 alpha occ, 20 beta occ)

 Alpha  occ. eigenvalues --  -10.44746 -10.44745 -10.44690 -10.44689 -10.41307
 Alpha  occ. eigenvalues --  -10.41306  -1.09893  -0.99649  -0.97270  -0.83278
 Alpha  occ. eigenvalues --   -0.83268  -0.74423  -0.68358  -0.67574  -0.65278
 Alpha  occ. eigenvalues --   -0.63494  -0.61047  -0.56837  -0.56618  -0.51141
 Alpha  occ. eigenvalues --   -0.47878
 Alpha virt. eigenvalues --   -0.25225  -0.22671  -0.10624  -0.07758  -0.05310
 Alpha virt. eigenvalues --   -0.04280  -0.01821  -0.00871   0.00401   0.08260
 Alpha virt. eigenvalues --    0.08579   0.09642   0.10056   0.25206   0.29439
 Alpha virt. eigenvalues --    0.31399   0.31852   0.34121   0.36475   0.36906
 Alpha virt. eigenvalues --    0.37451   0.38343   0.38500   0.39459   0.40284
 Alpha virt. eigenvalues --    0.43576   0.45334   0.52549   0.60260   0.60770
 Alpha virt. eigenvalues --    0.61287   0.62929   0.64337   0.70989   0.71650
 Alpha virt. eigenvalues --    0.71731   0.74333   0.85713   0.86949   0.90112
 Alpha virt. eigenvalues --    0.90952   0.98816   1.00856   1.05831   1.15646
 Alpha virt. eigenvalues --    1.17792   1.17972   1.18789   1.20601   1.20854
 Alpha virt. eigenvalues --    1.49713   1.52475   1.57000   1.65756   1.66784
 Alpha virt. eigenvalues --    1.68337   1.73545   1.74011   1.74167   1.74723
 Alpha virt. eigenvalues --    1.80258   1.82880   1.84586   2.04024   2.06015
 Alpha virt. eigenvalues --    2.12117   2.12667   2.14025   2.17682   2.18940
 Alpha virt. eigenvalues --    2.19096   2.22084   2.22451   2.24748   2.25480
 Alpha virt. eigenvalues --    2.28544   2.35165   2.36888   2.39005   2.41062
 Alpha virt. eigenvalues --    2.52629   2.57091   2.57724   2.79730   2.80863
 Alpha virt. eigenvalues --    2.95189   2.99029   2.99731   3.01110   3.14403
 Alpha virt. eigenvalues --    3.25310   3.26537   3.70063   3.88553   3.90763
 Alpha virt. eigenvalues --    3.92953   4.18629   4.20462   4.58339
  Beta  occ. eigenvalues --  -10.44304 -10.44303 -10.44252 -10.44250 -10.41463
  Beta  occ. eigenvalues --  -10.41462  -1.08758  -0.97673  -0.97028  -0.82708
  Beta  occ. eigenvalues --   -0.82377  -0.74165  -0.67883  -0.67164  -0.64793
  Beta  occ. eigenvalues --   -0.63478  -0.57727  -0.56637  -0.56323  -0.47270
  Beta virt. eigenvalues --   -0.41639  -0.21435  -0.21139  -0.10438  -0.05496
  Beta virt. eigenvalues --   -0.05056  -0.04232  -0.01054  -0.00739   0.00754
  Beta virt. eigenvalues --    0.08748   0.08784   0.10027   0.10356   0.25410
  Beta virt. eigenvalues --    0.30875   0.31655   0.33033   0.34430   0.37599
  Beta virt. eigenvalues --    0.38243   0.38423   0.38827   0.38857   0.40471
  Beta virt. eigenvalues --    0.40510   0.45633   0.45687   0.53548   0.60543
  Beta virt. eigenvalues --    0.61003   0.61366   0.63303   0.64325   0.71163
  Beta virt. eigenvalues --    0.71910   0.72371   0.74501   0.86611   0.87153
  Beta virt. eigenvalues --    0.90721   0.90982   0.99163   1.02443   1.07028
  Beta virt. eigenvalues --    1.17547   1.18130   1.19642   1.19672   1.20955
  Beta virt. eigenvalues --    1.21374   1.51458   1.52709   1.57335   1.66396
  Beta virt. eigenvalues --    1.67580   1.68460   1.73895   1.74747   1.75260
  Beta virt. eigenvalues --    1.75568   1.80924   1.84865   1.84936   2.06229
  Beta virt. eigenvalues --    2.06582   2.12479   2.12665   2.14334   2.18350
  Beta virt. eigenvalues --    2.18883   2.19283   2.22289   2.22978   2.25783
  Beta virt. eigenvalues --    2.25938   2.29233   2.36212   2.37068   2.39062
  Beta virt. eigenvalues --    2.42549   2.53376   2.57824   2.57840   2.79980
  Beta virt. eigenvalues --    2.80952   2.95964   2.99101   2.99875   3.01115
  Beta virt. eigenvalues --    3.14561   3.25632   3.26592   3.70353   3.89317
  Beta virt. eigenvalues --    3.92008   3.93146   4.19813   4.20623   4.58989

In the case of UHF wave functions, you specify alpha or beta using AMO= or BMO= when you run cubegen.