Merge pull request #2 from ifilot/develop

Patch 0.6.3

README.md

@@ -22,7 +22,8 @@ packages offer. It **does** offer a unique insight into a working code and a
 considerable effort was made in documenting everything.
 
 > [!TIP]  
-> Interested in other **education** quantum chemical codes? Have a look at the packages below.
+> Interested in other **education** quantum chemical codes? Have a look at the
+> packages below.
 > * [PyQInt](https://github.com/ifilot/pyqint) is a hybrid C++/Python (Cython)
 >   code for performing Hartree-Fock calculations. This code contains many
 >   relevant features, such a geometry optimization, orbital localization and
@@ -91,10 +92,10 @@ print("Total electronic energy: %f Ht" % en)
 * Contributions to PyDFT are always welcome and appreciated. Before doing so,
   please first read the [CONTRIBUTING](CONTRIBUTING.md) guide.
 * For reporting issues or problems with the software, you are kindly invited to
-  to open a [new issue with the bug label](https://github.com/ifilot/pydft/issues/new?labels=bug).
-* If you seek support in using PyDFT, please 
-  [open an issue with the question](https://github.com/ifilot/pydft/issues/new?labels=question)
-  label.
+  to open a [new issue with the bug
+  label](https://github.com/ifilot/pydft/issues/new?labels=bug).
+* If you seek support in using PyDFT, please [open an issue with the
+  question](https://github.com/ifilot/pydft/issues/new?labels=question) label.
 * If you wish to contact the developers, please send an e-mail to ivo@ivofilot.nl.
 
 ## License

meta.yaml

@@ -1,6 +1,6 @@
 package:
   name: "pydft"
-  version: "0.6.2"
+  version: "0.6.3"
 
 source:
   path: .

pydft/_version.py

@@ -1 +1 @@
-__version__ = '0.6.2'
+__version__ = '0.6.3'

pydft/dft.py

@@ -19,6 +19,7 @@ def __init__(self,
                  nshells:int = 32,
                  nangpts:int = 110,
                  lmax:int = 8,
+                 normalize:bool = True,
                  verbose:bool = False):
         """
         Constructs the DFT class
@@ -38,6 +39,8 @@ def __init__(self,
             number of angular sampling points, by default 110
         lmax : int, optional
             maximum value of l in the spherical harmonic expansion, by default 8
+        normalize: whether to perform intermediary normalization of the electron 
+            density
         verbose : bool, optional
             whether to provide verbose output, by default False
         """
@@ -51,6 +54,7 @@ def __init__(self,
         self.__nangpts = nangpts
         self.__lmax = lmax
         self.__functional = functional
+        self.__normalize = normalize
 
     def get_data(self) -> dict:
         """
@@ -237,7 +241,7 @@ def __iterate(self, niter, giis=True, mix=0.9):
         # calculate J and XC matrices based on the current electron
         # density estimate as captured in the density matrix P
         if np.any(self.__P):
-            self.__molgrid.build_density(self.__P, normalize=True)
+            self.__molgrid.build_density(self.__P, normalize=self.__normalize)
             self.__J = self.__calculate_J()
             self.__XC, self.__Exc = self.__calculate_XC()
 

pydft/moleculargrid.py

@@ -62,7 +62,7 @@ def initialize(self):
         self.__build_amplitudes()
         self.__is_initialized = True
 
-    def build_density(self, P:np.ndarray, normalize:bool=False):
+    def build_density(self, P:np.ndarray, normalize:bool=True):
         """
         Build the electron density from the density matrix
 
@@ -72,7 +72,7 @@ def build_density(self, P:np.ndarray, normalize:bool=False):
             density matrix
         normalize : bool, optional
             whether to normalize the density matrix based on the total number
-            of electrons, by default False
+            of electrons, by default True
         """
         # this function requires the amplitudes and molecular grid
         # to be built
@@ -94,7 +94,7 @@ def build_density(self, P:np.ndarray, normalize:bool=False):
         # perform optional normalization
         if normalize:
             nelec = np.sum(self.__mgw * self.__densities.flatten())
-            cn = nelec / self.__nelec # normalization constant
+            cn = self.__nelec / nelec  # normalization constant
             self.__densities *= cn
             self.__gradients *= cn
 

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pydft"
-version = "0.6.2"
+version = "0.6.3"
 authors = [
   { name="Ivo Filot", email="i.a.w.filot@tue.nl" }
 ]

tests/test_dft.py

@@ -29,7 +29,7 @@ def test_h2o(self):
         dft = DFT(mol, basis='sto3g')
         energy = dft.scf()
 
-        answer = -74.9393412503195
+        answer = -74.93750649033662
         np.testing.assert_almost_equal(energy, answer, 4)
 
     def test_co(self):
@@ -43,7 +43,7 @@ def test_co(self):
         dft = DFT(mol, basis='sto3g')
         energy = dft.scf()
 
-        answer = -111.1458218142885
+        answer = -111.14683799873168
         np.testing.assert_almost_equal(energy, answer, 4)
 
     def test_bf3(self):
@@ -57,7 +57,7 @@ def test_bf3(self):
         dft = DFT(mol, basis='sto3g')
         energy = dft.scf()
 
-        answer = -318.2583234959651
+        answer = -318.2753601057624
         np.testing.assert_almost_equal(energy, answer, 4)
 
     def test_ch4(self):
@@ -71,7 +71,7 @@ def test_ch4(self):
         dft = DFT(mol, basis='sto3g')
         energy = dft.scf()
 
-        answer = -39.801509613192096
+        answer = -39.80546943950668
         np.testing.assert_almost_equal(energy, answer, 4)
 
     def test_co2(self):
@@ -85,7 +85,7 @@ def test_co2(self):
         dft = DFT(mol, basis='sto3g')
         energy = dft.scf()
 
-        answer = -185.0106604454506
+        answer = -185.01476596292474
         np.testing.assert_almost_equal(energy, answer, 4)
 
     def test_ethylene(self):
@@ -99,7 +99,7 @@ def test_ethylene(self):
         dft = DFT(mol, basis='sto3g')
         energy = dft.scf()
 
-        answer = -77.16054913818746
+        answer = -77.16253628136911
         np.testing.assert_almost_equal(energy, answer, 4)
 
     def test_h2(self):
@@ -127,7 +127,7 @@ def test_lih(self):
         dft = DFT(mol, basis='sto3g')
         energy = dft.scf()
 
-        answer = -7.86536764070593
+        answer = -7.865828015750928
         np.testing.assert_almost_equal(energy, answer, 4)
 
     # def test_benzene(self):

tests/test_normalization.py

@@ -0,0 +1,57 @@
+import unittest
+import numpy as np
+from pydft import MoleculeBuilder, DFT
+
+class TestNormalization(unittest.TestCase):
+    """
+    Test whether normalization yields **exactly** the number of electrons of
+    the system
+    """
+
+    def test_helium(self):
+        """
+        Test DFT calculation of Helium atom
+        """
+        mol_builder = MoleculeBuilder()
+        mol = mol_builder.from_name('He')
+        Nelec = 2.0
+
+        # construct dft object
+        dft = DFT(mol, basis='sto3g', normalize=False)
+        dft.scf()
+        nelec = dft.get_molgrid_copy().count_electrons()
+
+        err_no_norm = np.abs(Nelec - nelec)
+
+        # construct dft object
+        dft = DFT(mol, basis='sto3g', normalize=True)
+        dft.scf()
+        nelec = dft.get_molgrid_copy().count_electrons()
+
+        err_norm = np.abs(Nelec - nelec)
+        self.assertTrue(err_norm < err_no_norm)
+        np.testing.assert_almost_equal(nelec, Nelec, 8)
+
+    def test_methane(self):
+        """
+        Test DFT calculation of methane molecule
+        """
+        mol_builder = MoleculeBuilder()
+        mol = mol_builder.from_name('CH4')
+        Nelec = 10.0
+
+        # construct dft object
+        dft = DFT(mol, basis='sto3g', normalize=False)
+        dft.scf()
+        nelec = dft.get_molgrid_copy().count_electrons()
+
+        err_no_norm = np.abs(Nelec - nelec)
+
+        # construct dft object
+        dft = DFT(mol, basis='sto3g', normalize=True)
+        dft.scf()
+        nelec = dft.get_molgrid_copy().count_electrons()
+
+        err_norm = np.abs(Nelec - nelec)
+        self.assertTrue(err_norm < err_no_norm)
+        np.testing.assert_almost_equal(nelec, Nelec, 8)