From 1e7035134dc6957ecb1ffae4c51ab0e8110564fb Mon Sep 17 00:00:00 2001 From: ifilot Date: Sat, 14 Dec 2024 20:05:08 +0100 Subject: [PATCH] Allowing user to set finite difference scheme --- pydft/_version.py | 2 +- pydft/atomicgrid.py | 86 ++++++++++++---------------- pydft/dft.py | 7 ++- pydft/moleculargrid.py | 7 ++- tests/test_finite_difference_grid.py | 36 ++++++++++++ 5 files changed, 87 insertions(+), 51 deletions(-) create mode 100644 tests/test_finite_difference_grid.py diff --git a/pydft/_version.py b/pydft/_version.py index 02f8497..a71c5c7 100644 --- a/pydft/_version.py +++ b/pydft/_version.py @@ -1 +1 @@ -__version__ = '0.6.4' +__version__ = '0.7.0' diff --git a/pydft/atomicgrid.py b/pydft/atomicgrid.py index fabe5d0..68e9ae7 100644 --- a/pydft/atomicgrid.py +++ b/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 diff --git a/pydft/dft.py b/pydft/dft.py index e4d703f..4ccc59a 100644 --- a/pydft/dft.py +++ b/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 diff --git a/pydft/moleculargrid.py b/pydft/moleculargrid.py index a863e9d..aad6f6a 100644 --- a/pydft/moleculargrid.py +++ b/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 diff --git a/tests/test_finite_difference_grid.py b/tests/test_finite_difference_grid.py new file mode 100644 index 0000000..29c07e2 --- /dev/null +++ b/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()