Develop (#23)

* Initial commit docs

* Adding deployment docs

* Adding deployment docs

* Expanding documentation

* Building docker image for docs

* Adding Dockerfile

* Adding example script

* Expanding documentation

* Expanding documentation

* Redesigning README.md

* Expanding documentation

* Initial commit Foster-Boys procedure

* Fixing requirements

* Adding scipy in build and host

* Adding scipy to environment.yml

* Adding Foster-Boys script

* Adding testing

* Initial commit COHP and adding RNG seed for FB localization

* Adding seed

* Add multiple runners

* Decreasing accuracy to two decimals for FB

* Putting RNG as class member

* Fixing nuclear derivatives

* Fixing meta.yaml

* Adding derivatives tests

* Adding error catching in diis

* Adding LiH example for COHP

* Adding molecule builder

* Adding molecules

* Fixing manifest file

* Fixing type

* Making test more lenient

* Update cohp_lih.py

* Adding bounds for Foster-Boys

* Refactoring HF derivatives

* Pre-seeding calculation

* Initial commit steepest descent

* Initial commit geometric derivatives parallel parser

* Update cohp_lih.py

* Initial commit OpenMP calculation of derivatives

* Adding steepest descent

* More stringent conversion criterium

* Reducing number of digits for MacOS

* Adding conjugate gradient algorithm

* Making test more lenient

* Fixing issue #17

* Tracking history in geometry optimization class

* Fixing bugs in HF algo

* Fixing bugs in GeometryOptimization and in MoleculeBuilder

* Fixing molecule builder test

* Forcing CH4

* Adding .xyz file to unit test

* Adding example scripts

* Adding nelec() method

* Adding charges

* Adding molecule rendering

* Incrementing version number

* Adding safeguard for PyTessel

* Removing tqdm and adding mendeleev in gaurds

* Changing versions of dependencies

* Adding tqdm and mendeleev as dependencies

* Removing mendeleev as package

* Fixing test environment

* Removing potential troublesome mendeleev package

* Testing compilation Python 3.11

* Specifying versions in conda_build_config.yml

COMPILATION.md

@@ -99,11 +99,11 @@ python3 setup.py build
 
 and install it locally
 ```
-pip install -e .
+pip3 install -e .
 ```
 
 and finally test it
 
 ```
-pytest tests/*
+pytest-3 tests/*
 ```

conda_build_config.yaml

@@ -2,3 +2,12 @@ python:
   - 3.8
   - 3.9
   - 3.10
+  - 3.11
+numpy:
+  - 1.22
+  - 1.22
+  - 1.22
+  - 1.24
+zip_keys:
+  - python
+  - numpy

environment.yml

@@ -5,3 +5,4 @@ dependencies:
   - conda-build
   - conda-verify
   - anaconda-client
+  - tqdm

examples/foster_boys.py

@@ -1,50 +1,69 @@
-from pyqint import Molecule, HF, PyQInt, FosterBoys
-import pyqint
+from pyqint import Molecule, PyQInt, FosterBoys, GeometryOptimization
 import numpy as np
-from pytessel import PyTessel
+import matplotlib.pyplot as plt
 
 #
 # Plot the isosurfaces for the CO molecule
 #
 
 def main():   
-    res = calculate_co(1.145414)
+    res = optimize_co()
     resfb = FosterBoys(res).run()
+    energies = resfb['orbe']
+    coeff = resfb['orbc']
+    cgfs = res['cgfs']
 
-    for i in range(len(res['cgfs'])):
-        build_isosurface('MO_%03i' % (i+1), 
-                         res['cgfs'],
-                         resfb['orbc'][:,i],
-                         0.1)
+    fig, ax = plt.subplots(2, 5, dpi=144, figsize=(12,5))           
+    for i,(chi,e) in enumerate(zip(coeff.transpose(), energies)):
+        res, x, y = build_contourplot(cgfs, chi, sz=3, plane='xz')
+        vmax = np.max(np.abs(res))
+        vmin = -vmax
+        ax[i//5,i%5].contourf(x, y, res, levels=15, cmap='PiYG',
+                              vmin=vmin, vmax=vmax)
+        ax[i//5,i%5].contour(x, y, res, levels=15, colors='black',
+                             vmin=vmin, vmax=vmax)
+        ax[i//5,i%5].set_xlabel('x [Bohr]')
+        ax[i//5,i%5].set_ylabel('z [Bohr]')
+        ax[i//5,i%5].set_title('Energy = %.4f Ht' % e)
+        ax[i//5,i%5].grid(linestyle='--', color='black', alpha=0.5)
+
+    plt.tight_layout()
+    plt.savefig('co_fb.jpg')
 
-def calculate_co(d):
+def optimize_co():
     """
-    Full function for evaluation
+    Optimization function for scipy.optimize.minimize
     """
     mol = Molecule()
-    mol.add_atom('C', 0.0, 0.0, -d/2, unit='angstrom')
-    mol.add_atom('O', 0.0, 0.0,  d/2, unit='angstrom')
-
-    result = HF().rhf(mol, 'sto3g')
-
-    return result
+    mol.add_atom('C', 0.0, 0.0, -0.6, unit='angstrom')
+    mol.add_atom('O', 0.0, 0.0,  0.6, unit='angstrom')
+    
+    res = GeometryOptimization().run(mol, 'sto3g')
+    
+    return res['data']
 
-def build_isosurface(filename, cgfs, coeff, isovalue, sz=5, npts=100):
-    # generate some data
-    isovalue = np.abs(isovalue)
+def build_contourplot(cgfs, coeff, sz=2, npts=50, plane='xy'):
     integrator = PyQInt()
-    grid = integrator.build_rectgrid3d(-sz, sz, npts)
-    scalarfield = np.reshape(integrator.plot_wavefunction(grid, coeff, cgfs), (npts, npts, npts))
-    unitcell = np.diag(np.ones(3) * 2 * sz)
-
-    pytessel = PyTessel()
-    vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), isovalue)
-    fname = filename + '_pos.ply'
-    pytessel.write_ply(fname, vertices, normals, indices)
 
-    vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), -isovalue)
-    fname = filename + '_neg.ply'
-    pytessel.write_ply(fname, vertices, normals, indices)
+    # build grid
+    x = np.linspace(-sz, sz, 50)
+    y = np.linspace(-sz, sz, 50)
+    xx, yy = np.meshgrid(x,y)
+    zz = np.zeros(len(x) * len(y))
+
+    if plane == 'xy':
+        points = [xx.flatten(), yy.flatten(), zz]
+    elif plane == 'xz':
+        points = [xx.flatten(), zz, yy.flatten()]
+    elif plane == 'yz':
+        points = [zz, xx.flatten(), yy.flatten()]
+    else:
+        raise Exception('Unknown plane: %s' % plane)
+
+    grid = np.vstack(points).reshape(3,-1).T
+    res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(y), len(x)))
+
+    return res, x, y
 
 if __name__ == '__main__':
     main()

examples/molecular_orbitals_co.py

@@ -1,80 +1,68 @@
-from pyqint import Molecule, HF, PyQInt
-import numpy as np
-import matplotlib.pyplot as plt
-
-#
-# Plot the isosurfaces for the CO molecule
-#
-
-def main():
-    # calculate sto-3g coefficients for co
-    result = calculate_co(1.145414)
-    energies = result['orbe']
-    coeff = result['orbc']
-
-    fig, ax = plt.subplots(2, 5, dpi=144, figsize=(12,5))           
-    for i,(chi,e) in enumerate(zip(coeff.transpose(), energies)):
-        res, x, y = build_contourplot(result['cgfs'], chi, sz=3, plane='xz')
-        vmax = np.max(np.abs(res))
-        vmin = -vmax
-        ax[i//5,i%5].contourf(x, y, res, levels=15, cmap='PiYG',
-                              vmin=vmin, vmax=vmax)
-        ax[i//5,i%5].contour(x, y, res, levels=15, colors='black',
-                             vmin=vmin, vmax=vmax)
-        ax[i//5,i%5].set_xlabel('x [Bohr]')
-        ax[i//5,i%5].set_ylabel('z [Bohr]')
-        ax[i//5,i%5].set_title('Energy = %.4f Ht' % e)
-        ax[i//5,i%5].grid(linestyle='--', color='black', alpha=0.5)
-
-    plt.tight_layout()
-    plt.savefig('co.jpg')
-
-def optimize_co(d):
-    """
-    Optimization function for scipy.optimize.minimize
-    """
-    mol = Molecule()
-    mol.add_atom('C', 0.0, 0.0, -d[0]/2, unit='angstrom')
-    mol.add_atom('O', 0.0, 0.0,  d[0]/2, unit='angstrom')
-
-    result = HF().rhf(mol, 'sto3g')
-
-    return result['energy']
-
-def calculate_co(d):
-    """
-    Full function for evaluation
-    """
-    mol = Molecule()
-    mol.add_atom('C', 0.0, 0.0, -d/2, unit='angstrom')
-    mol.add_atom('O', 0.0, 0.0,  d/2, unit='angstrom')
-
-    result = HF().rhf(mol, 'sto3g')
-
-    return result
-
-def build_contourplot(cgfs, coeff, sz=2, npts=50, plane='xy'):
-    integrator = PyQInt()
-
-    # build grid
-    x = np.linspace(-sz, sz, 50)
-    y = np.linspace(-sz, sz, 50)
-    xx, yy = np.meshgrid(x,y)
-    zz = np.zeros(len(x) * len(y))
-
-    if plane == 'xy':
-        points = [xx.flatten(), yy.flatten(), zz]
-    elif plane == 'xz':
-        points = [xx.flatten(), zz, yy.flatten()]
-    elif plane == 'yz':
-        points = [zz, xx.flatten(), yy.flatten()]
-    else:
-        raise Exception('Unknown plane: %s' % plane)
-
-    grid = np.vstack(points).reshape(3,-1).T
-    res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(y), len(x)))
-
-    return res, x, y
-
-if __name__ == '__main__':
-    main()
+from pyqint import Molecule, HF, PyQInt, GeometryOptimization
+import numpy as np
+import matplotlib.pyplot as plt
+
+#
+# Plot the isosurfaces for the CO molecule
+#
+
+def main():
+    # calculate sto-3g coefficients for co
+    result = optimize_co()
+    energies = result['orbe']
+    coeff = result['orbc']
+
+    fig, ax = plt.subplots(2, 5, dpi=144, figsize=(12,5))           
+    for i,(chi,e) in enumerate(zip(coeff.transpose(), energies)):
+        res, x, y = build_contourplot(result['cgfs'], chi, sz=3, plane='xz')
+        vmax = np.max(np.abs(res))
+        vmin = -vmax
+        ax[i//5,i%5].contourf(x, y, res, levels=15, cmap='PiYG',
+                              vmin=vmin, vmax=vmax)
+        ax[i//5,i%5].contour(x, y, res, levels=15, colors='black',
+                             vmin=vmin, vmax=vmax)
+        ax[i//5,i%5].set_xlabel('x [Bohr]')
+        ax[i//5,i%5].set_ylabel('z [Bohr]')
+        ax[i//5,i%5].set_title('Energy = %.4f Ht' % e)
+        ax[i//5,i%5].grid(linestyle='--', color='black', alpha=0.5)
+
+    plt.tight_layout()
+    plt.savefig('co.jpg')
+
+def optimize_co():
+    """
+    Optimization function for scipy.optimize.minimize
+    """
+    mol = Molecule()
+    mol.add_atom('C', 0.0, 0.0, -0.6, unit='angstrom')
+    mol.add_atom('O', 0.0, 0.0,  0.6, unit='angstrom')
+
+    res = GeometryOptimization().run(mol, 'sto3g')
+
+    return res['data']
+
+def build_contourplot(cgfs, coeff, sz=2, npts=50, plane='xy'):
+    integrator = PyQInt()
+
+    # build grid
+    x = np.linspace(-sz, sz, 50)
+    y = np.linspace(-sz, sz, 50)
+    xx, yy = np.meshgrid(x,y)
+    zz = np.zeros(len(x) * len(y))
+
+    if plane == 'xy':
+        points = [xx.flatten(), yy.flatten(), zz]
+    elif plane == 'xz':
+        points = [xx.flatten(), zz, yy.flatten()]
+    elif plane == 'yz':
+        points = [zz, xx.flatten(), yy.flatten()]
+    else:
+        raise Exception('Unknown plane: %s' % plane)
+
+    grid = np.vstack(points).reshape(3,-1).T
+    res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(y), len(x)))
+
+    return res, x, y
+
+if __name__ == '__main__':
+    main()

examples/molecular_orbitals_hco.py

@@ -0,0 +1,74 @@
+from pyqint import Molecule, HF, PyQInt, GeometryOptimization
+import numpy as np
+import matplotlib.pyplot as plt
+
+#
+# Plot the isosurfaces for the CO molecule
+#
+
+def main():
+    # calculate sto-3g coefficients for co
+    result = optimize_hco()
+    energies = result['orbe']
+    coeff = result['orbc']
+
+    fig, ax = plt.subplots(2, 5, dpi=144, figsize=(12,5))           
+    for i,(chi,e) in enumerate(zip(coeff.transpose(), energies)):
+
+        if i >= 10:
+            break
+
+        res, x, y = build_contourplot(result['cgfs'], chi, sz=3, plane='xz')
+        vmax = np.max(np.abs(res))
+        vmin = -vmax
+        ax[i//5,i%5].contourf(x, y, res, levels=15, cmap='PiYG',
+                              vmin=vmin, vmax=vmax)
+        ax[i//5,i%5].contour(x, y, res, levels=15, colors='black',
+                             vmin=vmin, vmax=vmax)
+        ax[i//5,i%5].set_xlabel('x [Bohr]')
+        ax[i//5,i%5].set_ylabel('z [Bohr]')
+        ax[i//5,i%5].set_title('Energy = %.4f Ht' % e)
+        ax[i//5,i%5].grid(linestyle='--', color='black', alpha=0.5)
+
+    plt.tight_layout()
+    plt.savefig('hco.jpg')
+
+def optimize_hco():
+    """
+    Optimization function for scipy.optimize.minimize
+    """
+    mol = Molecule()
+    mol.add_atom('O', -0.12770367,  -0.00000000,   1.23419284)
+    mol.add_atom('C', -0.50625665,   0.00000000,  -1.09149431)
+    mol.add_atom('H',  1.57882331,  -0.00000000,  -2.06681794)
+    mol.set_charge(-1)
+
+    res = GeometryOptimization().run(mol, 'p321')
+
+    return res['data']
+
+def build_contourplot(cgfs, coeff, sz=2, npts=50, plane='xy'):
+    integrator = PyQInt()
+
+    # build grid
+    x = np.linspace(-sz, sz, 50)
+    y = np.linspace(-sz, sz, 50)
+    xx, yy = np.meshgrid(x,y)
+    zz = np.zeros(len(x) * len(y))
+
+    if plane == 'xy':
+        points = [xx.flatten(), yy.flatten(), zz]
+    elif plane == 'xz':
+        points = [xx.flatten(), zz, yy.flatten()]
+    elif plane == 'yz':
+        points = [zz, xx.flatten(), yy.flatten()]
+    else:
+        raise Exception('Unknown plane: %s' % plane)
+
+    grid = np.vstack(points).reshape(3,-1).T
+    res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(y), len(x)))
+
+    return res, x, y
+
+if __name__ == '__main__':
+    main()

meta.yaml

@@ -1,6 +1,6 @@
 package:
   name: "pyqint"
-  version: "0.12.2.0"
+  version: "0.13.0.0"
 
 source:
   path: .
@@ -23,6 +23,7 @@ requirements:
     - python
     - numpy
     - scipy
+    - tqdm
 
 test:
   requires:

pyqint/__init__.py

@@ -7,5 +7,6 @@
 from .cohp import COHP
 from .molecule_builder import MoleculeBuilder
 from .geometry_optimization import GeometryOptimization
+from .blenderrender import BlenderRender
 
 from ._version import __version__

pyqint/_version.py

@@ -1,2 +1,2 @@
-__version__ = "0.12.2.0"
+__version__ = "0.14.0.0"