Allowing user to set finite difference scheme

pydft/_version.py

@@ -1 +1 @@
-__version__ = '0.6.4'
+__version__ = '0.7.0'

pydft/atomicgrid.py

@@ -8,7 +8,8 @@
 import math
 
 class AtomicGrid:
-    def __init__(self, at, nshells:int=32, nangpts:int=110, lmax:int=8):
+    def __init__(self, at, nshells:int=32, nangpts:int=110, lmax:int=8,
+                 fdpts:int=7):
         """
         Initialize the class
         
@@ -17,9 +18,16 @@ def __init__(self, at, nshells:int=32, nangpts:int=110, lmax:int=8):
         nangpts: number of angular points
         """
         self.__atom = at
-        self.__cp = at[0]   # center of the atomic grid (= position of the atom)
-        self.__aidx = at[1] # element id
-        self.__nshells = nshells
+        self.__cp = at[0]           # center of the atomic grid (= position of the atom)
+        self.__aidx = at[1]         # element id
+        self.__nshells = nshells    # number of radial shells
+        self.__fdpts = fdpts        # number of grid points used in finite difference scheme
+
+        # perform some parameter checking
+        if self.__fdpts < 3:
+            raise Exception('Number of grid points in stencil must at least be 3')
+        if self.__fdpts % 2 != 1:
+            raise Exception('Only odd number of grid points are allowed')
 
         # set coefficients and parameters
         self.__set_bragg_slater_radius()
@@ -370,12 +378,20 @@ def __build_finite_difference_matrix(self, r, rm):
         """
         Construct the finite difference matrix to solve for Ulm
         
+        The discretization points are automatically generated using the
+        '__calculate_fd_coeff' function. The user can specify the number of
+        grid points used via npts. The script always constructs a
+        npts-diagonal matrix, e.g. a heptadiagonal matrix when n=7 or a
+        pentadiagonal matrix when n-5.
+        
         r : vector of distances
         rm: half of the Bragg-Slater radius
+        nrpts: number of stencil points
         """
         N = len(r)
         A = np.zeros((N+2, N+2))
-        h = 1.0 / float(N+1)
+        h = 1.0 / float(N+1)       
+        Np = self.__fdpts//2
 
         for i in range(0, N+2):
             c1 = 0.0
@@ -385,59 +401,33 @@ def __build_finite_difference_matrix(self, r, rm):
                 c1 = self.__dzdrsq(r[i-1], rm)
                 c2 = self.__d2zdr2(r[i-1], rm)
 
+            # first point
             if i == 0:
                 A[0,0] = 1.0
                 continue
 
-            if i == 1:
-                c1 /= 12.0 * h * h
-                c2 /= 12.0 * h
-                A[i,0] = 11.0 * c1 -3.0 * c2
-                A[i,1] = -20.0 * c1 - 10.0 * c2
-                A[i,2] = 6.0 * c1 + 18.0 * c2
-                A[i,3] = 4.0 * c1 - 6.0 * c2
-                A[i,4] = -1.0 * c1 + 1.0 * c2
-                continue
-
-            if i == 2:
-                c1 /= 12.0 * h * h
-                c2 /= 60.0 * h
-                A[i,0] = -1.0 * c1 + 3.0 * c2
-                A[i,1] = 16.0 * c1 - 30.0 * c2
-                A[i,2] = -30.0 * c1 - 20.0 * c2
-                A[i,3] = 16.0 * c1 + 60.0 * c2
-                A[i,4] = -1.0 * c1 - 15.0 * c2
-                A[i,5] = 0.0 * c1 + 2.0 * c2
-                continue
-
-            if i == N-1:
-                c1 /= 12.0 * h * h
-                c2 /= 60.0 * h
-                A[i,N-4] = 0.0 * c1 - 2.0 * c2
-                A[i,N-3] = -1.0 * c1 + 15.0 * c2
-                A[i,N-2] = 16.0 * c1 - 60.0 * c2
-                A[i,N-1] = -30.0 * c1 + 20.0 * c2
-                A[i,N] = 16.0 * c1 + 30.0 * c2
-                A[i,N+1] = -1.0 * c1 - 3.0 * c2
+            # last point
+            if i == N+1:
+                A[i,i] = 1.0
                 continue
 
-            if i == N:
-                c1 /= 12.0 * h * h
-                c2 /= 12.0 * h
-                A[i,N-3] = -1.0 * c1 - 1.0 * c2
-                A[i,N-2] = 4.0 * c1 + 6.0 * c2
-                A[i,N-1] = 6.0 * c1 - 18.0 * c2
-                A[i,N] = -20.0 * c1 + 10.0 * c2
-                A[i,N+1] = 11.0 * c1 + 3.0 * c2
+            # left edge, excluding first point
+            if i < Np:
+                c1 *= self.__calculate_fd_coeff(np.arange(0,i+Np+1)-i, 2)
+                c2 *= self.__calculate_fd_coeff(np.arange(0,i+Np+1)-i, 1)
+                A[i,0:Np+i+1] = c1 / h**2 + c2 / h
                 continue
 
-            if i == N+1:
-                A[i,i] = 1.0
+            # right edge, excluding last point
+            if i > (N - Np + 1):
+                c1 *= self.__calculate_fd_coeff(np.arange(i-Np,N+2)-i, 2)
+                c2 *= self.__calculate_fd_coeff(np.arange(i-Np,N+2)-i, 1)
+                A[i,i-Np:N+2] = c1 / h**2 + c2 / h
                 continue
 
-            c1 *= self.__calculate_fd_coeff(np.arange(-3,4), 2)
-            c2 *= self.__calculate_fd_coeff(np.arange(-3,4), 1)
-            A[i,(i-3):(i+4)] = c1 / h**2 + c2 / h
+            c1 *= self.__calculate_fd_coeff(np.arange(-Np,Np+1), 2)
+            c2 *= self.__calculate_fd_coeff(np.arange(-Np,Np+1), 1)
+            A[i,(i-Np):(i+Np+1)] = c1 / h**2 + c2 / h
 
         return A
 

pydft/dft.py

@@ -19,6 +19,7 @@ def __init__(self,
                  nshells:int = 32,
                  nangpts:int = 110,
                  lmax:int = 8,
+                 fdpts:int=7,
                  normalize:bool = True):
         """
         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
+        fdpts: int, optional
+            number of grid point in finite difference scheme, by default 7
         normalize: whether to perform intermediary normalization of the electron 
             density
         verbose : bool, optional
@@ -50,6 +53,7 @@ def __init__(self,
         self.__itermax = 100
         self.__nshells = nshells
         self.__nangpts = nangpts
+        self.__fdpts = fdpts
         self.__lmax = lmax
         self.__functional = functional
         self.__normalize = normalize
@@ -310,8 +314,9 @@ def __setup(self):
         self.__molgrid = MolecularGrid(self.__nuclei, 
                                        self.__cgfs, 
                                        nshells=self.__nshells, 
-                                       nangpts=self.__nangpts, 
+                                       nangpts=self.__nangpts,
                                        lmax=self.__lmax,
+                                       fdpts=self.__fdpts,
                                        functional=self.__functional)
         self.__molgrid.initialize() # molecular grid uses late initialization
 

pydft/moleculargrid.py

@@ -16,6 +16,7 @@ def __init__(self,
                  nshells:int=32, 
                  nangpts:int=110, 
                  lmax:int=8,
+                 fdpts:int=7,
                  functional:str='svwn5'):
         """
         Construct MolecularGrid
@@ -32,6 +33,8 @@ def __init__(self,
             number of angular sampling points per shell, by default 110
         lmax : int, optional
             maximum value for l for projection of spherical harmonics, by default 8
+        fdpts: int, optional
+            number of grid point in finite difference scheme, by default 7
         functional : str, optional
             exchange-correlation functional, by default 'svwn5'
         """
@@ -42,6 +45,7 @@ def __init__(self,
         self.__atoms = atoms
         self.__nelec = np.sum([nuc[1] for nuc in self.__atoms])
         self.__lmax = lmax
+        self.__fdpts = fdpts
         self.__nshells = nshells
         self.__nangpts = nangpts
         self.__basis = cgfs
@@ -741,7 +745,8 @@ def __build_molecular_grid(self):
             self.__atomgrids.append(AtomicGrid(atom, 
                                                self.__nshells, 
                                                self.__nangpts,
-                                               self.__lmax)
+                                               self.__lmax,
+                                               self.__fdpts)
                                     )
         self.construct_times['atomic_grids'] = time.time() - st
 

tests/test_finite_difference_grid.py

@@ -0,0 +1,36 @@
+import unittest
+import sys
+import os
+import numpy as np
+
+# add a reference to load the pyDFT module
+sys.path.append(os.path.join(os.path.dirname(__file__), '..'))
+
+from pydft import MoleculeBuilder, DFT
+
+class TestDFT(unittest.TestCase):
+
+    def test_co(self):
+        """
+        Test DFT calculation of CO molecule
+        """
+        mol_builder = MoleculeBuilder()
+        mol = mol_builder.from_name('CO')
+
+        answers = [
+            -110.96476654905626,
+            -111.14751240002687,
+            -111.14683799873168,
+            -111.15048475412650,
+            -111.14985895386218,
+            -111.14831640666813,
+        ]
+
+        for fdpts,answer in zip([3,5,7,9,11,13], answers):
+            dft = DFT(mol, basis='sto3g', fdpts=fdpts)
+            energy = dft.scf()
+
+            np.testing.assert_almost_equal(energy, answer, 4)
+
+if __name__ == '__main__':
+    unittest.main()