Polishing class and example script

examples/blender_render.py

@@ -1,26 +1,120 @@
 from pyqint import Molecule, PyQInt, FosterBoys, GeometryOptimization, HF
-from pyqint import BlenderRender
+from pyqint import MoleculeBuilder, BlenderRender
 import pyqint
 import numpy as np
 import matplotlib.pyplot as plt
 import os
+import subprocess
 
 #
-# Plot the isosurfaces for the CO molecule
+# Plot the isosurfaces for a number of molecules, prior and after
+# orbital localization. Note that this script has been designed to be executed
+# in a Linux Ubuntu 22.04 LTS (WSL) container.
 #
 
-print(pyqint.__version__)
 outpath = os.path.dirname(__file__)
 
-def main():   
+def main():
+    # build_orbitals_co()
+    build_orbitals_ch4()
+
+def build_orbitals_co():
+    """
+    Build a montage image of the canonical and localized molecular orbitals
+    of CO
+    """
+    molname = 'CO'
     mol = Molecule('CO')
     mol.add_atom('C', 0, 0, -1.08232106)
     mol.add_atom('O', 0, 0, 1.08232106)
     res = HF().rhf(mol, 'sto3g')
     resfb = FosterBoys(res).run()
 
+    build(molname, res, resfb, nrows=2)
+
+def build_orbitals_ch4():
+    """
+    Build a montage image of the canonical and localized molecular orbitals
+    of CH4
+    """
+    molname = 'CH4'
+    mol = MoleculeBuilder().from_name('ch4')
+    res = GeometryOptimization().run(mol, 'sto3g')['data']
+    resfb = FosterBoys(res).run()
+
+    build(molname, res, resfb, nrows=3)
+
+def build(molname, res, resfb, nrows=2):
+    """
+    Build isosurfaces, montage and print energies to a file
+
+    :param      molname:  Name of the molecule
+    :type       molname:  string
+    :param      res:      Results of a Hartree-Fock calculation
+    :type       res:      dictionary
+    :param      resfb:    Results of a Foster-Boys localization
+    :type       resfb:    dictionary
+    :param      nrows:    Number of rows in the montage
+    :type       nrows:    int
+    """
+    #build_isosurfaces(molname, res, resfb)
+    montage(molname, nrows)
+    store_energies(os.path.join(os.path.dirname(__file__), 'MO_%s_energies.txt' % molname), res['orbe'], resfb['orbe'])
+
+def build_isosurfaces(molname, res, resfb):
+    """
+    Builds isosurfaces.
+
+    :param      molname:  Name of the molecule
+    :type       molname:  string
+    :param      res:      Result of a Hartree-Fock calculation
+    :type       res:      dictionary
+    :param      resfb:    Result of a Foster-Boys localization
+    :type       resfb:    dictionary
+    """
     br = BlenderRender()
-    br.render_molecular_orbitals(res['mol'], res['cgfs'], resfb['orbc'], outpath)
+    br.render_molecular_orbitals(res['mol'], res['cgfs'], res['orbc'], outpath,
+        prefix='MO_CAN_%s' % molname)
+
+    br.render_molecular_orbitals(resfb['mol'], res['cgfs'], resfb['orbc'], outpath,
+        prefix='MO_FB_%s' % molname)
+
+def montage(molname, nrows=2):
+    """
+    Produce a montage of the images
+
+    :param      molname:  Name of the molecule
+    :type       molname:  string
+    :param      nrows:    Number of rows in the montage
+    :type       nrows:    int
+    """
+    out = subprocess.check_output(
+        ['montage', 'MO_CAN_%s_????.png' % molname, '-tile', 'x%i' % nrows, '-geometry', '128x128+2+2', 'MO_%s_CAN.png' % molname],
+        cwd=os.path.dirname(__file__)
+    )
+
+    out = subprocess.check_output(
+        ['montage', 'MO_FB_%s_????.png' % molname, '-tile', 'x%i' % nrows, '-geometry', '128x128+2+2', 'MO_%s_FB.png' % molname],
+        cwd=os.path.dirname(__file__)
+    )
+
+def store_energies(filename, orbe, orbe_fb):
+    """
+    Stores energies.
+
+    :param      filename:  Name of the molecule
+    :type       filename:  string
+    :param      orbe:      Array of the canonical molecular orbital energies
+    :type       orbe:      list or numpy array
+    :param      orbe_fb:   Array of the localized molecular orbital energies
+    :type       orbe_fb:   list or numpy array
+    """
+    f = open(filename, 'w')
+    f.write('id    localized    canonical\n')
+    f.write('----------------------------\n')
+    for i in range(len(orbe)):
+        f.write("%02i %12.4f %12.4f\n" % (i+1, orbe[i], orbe_fb[i]))
+    f.close()
 
 if __name__ == '__main__':
     main()

pyqint/blenderrender.py

@@ -40,7 +40,8 @@ def __init__(self):
 
     def render_molecular_orbitals(self, molecule, cgfs, orbc, outpath,
                                   mo_indices=None, sz=5, isovalue=0.03,
-                                  npts=100, negcol='E72F65', poscol='3F9EE7'):
+                                  prefix='MO', npts=100,
+                                  negcol='E72F65', poscol='3F9EE7'):
         if mo_indices is None: # render all orbitals
             mo_indices = np.arange(0, len(orbc))
 
@@ -54,12 +55,12 @@ def render_molecular_orbitals(self, molecule, cgfs, orbc, outpath,
         for idx in pbar:
             # build isosurfaces
             pbar.set_description('Producing isosurfaces (#%i)' % (idx+1))
-            plyfile = os.path.join(tempdir, 'MO_%04i' % (idx+1))
+            plyfile = os.path.join(tempdir, '%s_%04i' % (prefix,idx+1))
             plypos, plyneg = self.__build_isosurface(plyfile, cgfs, orbc[:,idx], isovalue, sz, npts)
 
             # execute blender
             pbar.set_description('Producing molecular orbital (#%i)' % (idx+1))
-            outfile = os.path.join(outpath, 'MO_%04i.png' % (idx+1))
+            outfile = os.path.join(outpath, '%s_%04i.png' % (prefix,idx+1))
             logoutput = self.__run_blender(plypos, plyneg, xyzfile, outfile, tempdir, negcol, poscol)
 
             self.log.append("### START LOG: MOLECULAR ORBITAL %i ###" % (idx+1))