From 9ed023d697c787f712398a5b8a384655b5fa239c Mon Sep 17 00:00:00 2001 From: ifilot Date: Mon, 9 Dec 2024 22:03:09 +0100 Subject: [PATCH 1/5] Integration of automated FD schemes --- examples/h2.py | 2 +- pydft/atomicgrid.py | 26 +++++++++++++++++--------- pydft/dft.py | 10 ++++------ 3 files changed, 22 insertions(+), 16 deletions(-) diff --git a/examples/h2.py b/examples/h2.py index b9a6560..9aabaa8 100644 --- a/examples/h2.py +++ b/examples/h2.py @@ -14,5 +14,5 @@ CO = MoleculeBuilder().from_name("H2") dft = DFT(CO, basis='sto3g') -en = dft.scf(1e-4) +en = dft.scf(1e-4, verbose=True) print("Total electronic energy: %f Ht" % en) \ No newline at end of file diff --git a/pydft/atomicgrid.py b/pydft/atomicgrid.py index 144d4c0..fabe5d0 100644 --- a/pydft/atomicgrid.py +++ b/pydft/atomicgrid.py @@ -5,6 +5,7 @@ from . import bragg_slater from .angulargrid import AngularGrid from .spherical_harmonics import spherical_harmonic +import math class AtomicGrid: def __init__(self, at, nshells:int=32, nangpts:int=110, lmax:int=8): @@ -434,18 +435,25 @@ def __build_finite_difference_matrix(self, r, rm): A[i,i] = 1.0 continue - c1 /= 180.0 * h * h - c2 /= 60.0 * h - A[i,i-3] = 2.0 * c1 - 1.0 * c2 - A[i,i-2] = -27.0 * c1 + 9.0 * c2 - A[i,i-1] = 270.0 * c1 - 45.0 * c2 - A[i,i] = -490.0 * c1 - A[i,i+1] = 270.0 * c1 + 45.0 * c2 - A[i,i+2] = -27.0 * c1 - 9.0 * c2 - A[i,i+3] = 2.0 * c1 + 1.0 * c2 + 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 return A + def __calculate_fd_coeff(self, x, k): + """ + Calculate finite difference coefficients + + x: sampling points + k: derivative order + """ + T = np.vander(x, increasing=True).transpose() + b = np.zeros(len(x)) + b[k] = math.factorial(k) + + return np.linalg.solve(T,b) + def __set_bragg_slater_radius(self): """ Set the Bragg-Slater radius for the atom diff --git a/pydft/dft.py b/pydft/dft.py index c01b5b1..e4d703f 100644 --- a/pydft/dft.py +++ b/pydft/dft.py @@ -19,8 +19,7 @@ def __init__(self, nshells:int = 32, nangpts:int = 110, lmax:int = 8, - normalize:bool = True, - verbose:bool = False): + normalize:bool = True): """ Constructs the DFT class @@ -47,7 +46,6 @@ def __init__(self, self.__mol = mol self.__integrator = PyQInt() self.__basis = basis - self.__verbose = verbose self.__time_stats = {} self.__itermax = 100 self.__nshells = nshells @@ -155,7 +153,7 @@ def get_gradient_at_points(self, spoints:np.ndarray) -> np.ndarray: return self.__molgrid.get_gradient_at_points(spoints, self.__P) - def scf(self, tol:float=1e-5) -> float: + def scf(self, tol:float=1e-5, verbose:bool=False) -> float: """ Perform the self-consistent field procedure @@ -187,7 +185,7 @@ def scf(self, tol:float=1e-5) -> float: itertime = stop - start self.__time_stats['iterations'].append(itertime) - if self.__verbose: + if verbose: print('%03i | Energy: %12.6f | %0.4f ms' % (niter+1, energy, itertime)) if niter > 2: @@ -200,7 +198,7 @@ def scf(self, tol:float=1e-5) -> float: continue # terminate self-convergence cycle - if self.__verbose: + if verbose: print("Stopping SCF cycle, convergence reached.") # update density matrix from last found coefficient matrix From 1e7035134dc6957ecb1ffae4c51ab0e8110564fb Mon Sep 17 00:00:00 2001 From: ifilot Date: Sat, 14 Dec 2024 20:05:08 +0100 Subject: [PATCH 2/5] 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() From a68efc4b2eb2cae3f5260b913542073042efa667 Mon Sep 17 00:00:00 2001 From: ifilot Date: Sat, 14 Dec 2024 20:53:40 +0100 Subject: [PATCH 3/5] Adding time analysis --- examples/benzene.py | 19 +++++++++++++++++++ examples/co.py | 5 +++-- meta.yaml | 2 +- pydft/dft.py | 43 ++++++++++++++++++++++++++++++++++++++---- pydft/moleculargrid.py | 4 ++-- pyproject.toml | 2 +- 6 files changed, 65 insertions(+), 10 deletions(-) create mode 100644 examples/benzene.py diff --git a/examples/benzene.py b/examples/benzene.py new file mode 100644 index 0000000..c40f368 --- /dev/null +++ b/examples/benzene.py @@ -0,0 +1,19 @@ +# -*- coding: utf-8 -*- +import os,sys + +# add a reference to load the module +ROOT = os.path.dirname(__file__) +sys.path.insert(1, os.path.join(ROOT, '..')) + +from pydft import MoleculeBuilder, DFT + +# +# Example: Calculate total electronic energy for CO using standard +# settings. +# + +CO = MoleculeBuilder().from_name("benzene") +dft = DFT(CO, basis='sto3g') +en = dft.scf(1e-4, verbose=True) +print("Total electronic energy: %f Ht" % en) +dft.print_time_statistics() \ No newline at end of file diff --git a/examples/co.py b/examples/co.py index f00b1c5..e11b93b 100644 --- a/examples/co.py +++ b/examples/co.py @@ -14,5 +14,6 @@ CO = MoleculeBuilder().from_name("CO") dft = DFT(CO, basis='sto3g') -en = dft.scf(1e-4) -print("Total electronic energy: %f Ht" % en) \ No newline at end of file +en = dft.scf(1e-4, verbose=True) +print("Total electronic energy: %f Ht" % en) +dft.print_time_statistics() \ No newline at end of file diff --git a/meta.yaml b/meta.yaml index 4dd9d02..4786f34 100644 --- a/meta.yaml +++ b/meta.yaml @@ -1,6 +1,6 @@ package: name: "pydft" - version: "0.6.4" + version: "0.7.0" source: path: . diff --git a/pydft/dft.py b/pydft/dft.py index 4ccc59a..c688c23 100644 --- a/pydft/dft.py +++ b/pydft/dft.py @@ -57,6 +57,13 @@ def __init__(self, self.__lmax = lmax self.__functional = functional self.__normalize = normalize + + # keep track of time + self.calctimes = { + 'density_hartree': [], + 'calculate_J': [], + 'calculate_XC': [], + } def get_data(self) -> dict: """ @@ -87,6 +94,7 @@ def get_data(self) -> dict: * :code:`orbc`: coeffient matrix (duplicate) * :code:`orbe`: molecular orbital eigenvalues * :code:`enucrep`: electrostatic repulsion of the nuclei + * :code:`timedata`: computation time for various parts of the calculation """ data = { @@ -106,6 +114,10 @@ def get_data(self) -> dict: 'orbc': self.__C, # coeffient matrix (duplicate) 'orbe': self.__e, # molecular orbital eigenvalues 'enucrep': self.__enuc, # electrostatic repulsion of the nuclei + + # also output computation times + 'timedata': {'construct_times': self.__molgrid.construct_times, + 'calc_times': self.calctimes}, } return data @@ -179,6 +191,7 @@ def scf(self, tol:float=1e-5, verbose:bool=False) -> float: # start SCF iterative procedure nitfin = 0 + ediff = 0 for niter in range(0, self.__itermax): start = time.time() energy = self.__iterate(niter, @@ -190,15 +203,15 @@ def scf(self, tol:float=1e-5, verbose:bool=False) -> float: self.__time_stats['iterations'].append(itertime) if verbose: - print('%03i | Energy: %12.6f | %0.4f ms' % (niter+1, energy, itertime)) + print('%03i | E = %12.6f | dE = %5.4e | %0.4f ms' % (niter+1, energy, ediff, itertime)) - if niter > 2: + if niter > 0: ediff = np.abs(energy - self.__energies[-2]) + if niter > 2: if ediff < tol: # terminate giis self-convergence and continue with mixing nitfin += 1 - - if nitfin < 3: + if nitfin < 2: continue # terminate self-convergence cycle @@ -211,6 +224,19 @@ def scf(self, tol:float=1e-5, verbose:bool=False) -> float: return energy + def print_time_statistics(self): + print('-- Construction times --') + print('Atomic grids: %.4f s' % self.__molgrid.construct_times['atomic_grids']) + print('Fuzzy cell decomposition: %.4f s' % self.__molgrid.construct_times['fuzzy_cell_decomposition']) + print('Spherical harmonics: %.4f s' % self.__molgrid.construct_times['spherical_harmonics']) + print('Nuclear distance and potential: %.4f s' % self.__molgrid.construct_times['nuclear_distance_and_potential']) + print('Basis set amplitudes: %.4f s' % self.__molgrid.construct_times['basis_function_amplitudes']) + print() + print('-- Calculation times --') + print('Classical e-e repulsion matrix (J): %.4f s' % np.average(self.calctimes['calculate_J'])) + print('Electron density and Hartree potential (U): %.4f s' % np.average(self.calctimes['density_hartree'])) + print('Exchange-correlation matrices (XC): %.4f s' % np.average(self.calctimes['calculate_XC'])) + def get_construction_times(self) -> dict: """ Return construct times dictionary from MolecularGrid to get insights @@ -243,9 +269,18 @@ 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): + st = time.time() self.__molgrid.build_density(self.__P, normalize=self.__normalize) + self.calctimes['density_hartree'].append(time.time() - st) + + st = time.time() self.__J = self.__calculate_J() + self.calctimes['calculate_J'].append(time.time() - st) + + st = time.time() self.__XC, self.__Exc = self.__calculate_XC() + self.calctimes['calculate_XC'].append(time.time() - st) + # calculate Fock matrix self.__F = self.__H + self.__J + self.__XC diff --git a/pydft/moleculargrid.py b/pydft/moleculargrid.py index aad6f6a..7ad8e19 100644 --- a/pydft/moleculargrid.py +++ b/pydft/moleculargrid.py @@ -845,9 +845,9 @@ def __build_molecular_grid(self): self.__ylmgpts = np.ndarray((len(self.__atoms), (self.__lmax+1)**2, np.prod(self.__mweights.shape))) - for i,at in enumerate(self.__atoms): + for i,at in enumerate(self.__atoms): # loop over atoms lmctr = 0 - for l in range(0, self.__lmax+1): + for l in range(0, self.__lmax+1): # loop over spherical harmonics for m in range(-l, l+1): self.__ylmgpts[i,lmctr,:] = spherical_harmonic(l, m, \ self.__theta_gridpoints[i,:], diff --git a/pyproject.toml b/pyproject.toml index 33a75cb..1c47938 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "pydft" -version = "0.6.4" +version = "0.7.0" authors = [ { name="Ivo Filot", email="ivo@ivofilot.nl" } ] From 6753e5ca06be862416da5246f11912278e551081 Mon Sep 17 00:00:00 2001 From: ifilot Date: Sat, 14 Dec 2024 21:06:32 +0100 Subject: [PATCH 4/5] Improving computation speed --- pydft/atomicgrid.py | 35 ++++++++++++++++------ pydft/dft.py | 2 +- pydft/moleculargrid.py | 66 ++++++++++++++++++++---------------------- 3 files changed, 58 insertions(+), 45 deletions(-) diff --git a/pydft/atomicgrid.py b/pydft/atomicgrid.py index 68e9ae7..dc6eca3 100644 --- a/pydft/atomicgrid.py +++ b/pydft/atomicgrid.py @@ -187,14 +187,17 @@ def get_gradient_squared(self): """ grad = self.get_gradient() #return np.linalg.norm(grad, axis=1) - return np.einsum('ij,ij->i', grad, grad) + return np.einsum('ij,ij->i', grad, grad, optimize=True) def count_electrons(self): """ Sum over the electron density to obtain the number of electrons for this atomic cell """ - return np.einsum('ij,ij,ij', self.__edens, self.__wgrid, self.__mweights) + return np.einsum('ij,ij,ij', self.__edens, + self.__wgrid, + self.__mweights, + optimize=True) def get_dfa_exchange(self): """ @@ -202,14 +205,20 @@ def get_dfa_exchange(self): """ alpha = 2.0 / 3.0 fac = -2.25 * alpha * np.power(0.75 / np.pi, 1.0 / 3.0) - return fac * np.einsum('ij,ij,ij', np.power(self.__edens, 4/3), self.__wgrid, self.__mweights) + return fac * np.einsum('ij,ij,ij', np.power(self.__edens, 4/3), + self.__wgrid, + self.__mweights, + optimize=True) def get_dfa_kinetic(self): """ Get density functional approximation of the kinetic energy for this atomic cell """ Ckin = 3.0 / 40.0 * (3 / np.pi)**(2/3) * (2 * np.pi)**2 - return Ckin * np.einsum('ij,ij,ij', np.power(self.__edens, 5/3), self.__wgrid, self.__mweights) + return Ckin * np.einsum('ij,ij,ij', np.power(self.__edens, 5/3), + self.__wgrid, + self.__mweights, + optimize=True) def get_dfa_nuclear_local(self): """ @@ -231,7 +240,8 @@ def perform_spherical_harmonic_expansion(self, f): self.__ylm, f, self.__wangpts, - self.__mweights) * 4.0 * np.pi + self.__mweights, + optimize=True) * 4.0 * np.pi def get_bragg_slater_radius(self) -> float: """ @@ -307,9 +317,12 @@ def calculate_coulomb_energy(self): """ # build Hartree potential pot = 1.0 / self.__rr - self.__htpot = np.einsum('ij,jk,i,ik->ik', self.__ulm, self.__ylm, pot, self.__edens) + self.__htpot = np.einsum('ij,jk,i,ik->ik', self.__ulm, + self.__ylm, pot, + self.__edens, + optimize=True) - return np.einsum('ij,ij', self.__htpot, self.__wgrid) + return np.einsum('ij,ij', self.__htpot, self.__wgrid, optimize=True) def calculate_coulomb_energy_interpolation(self): """ @@ -323,9 +336,13 @@ def calculate_coulomb_energy_interpolation(self): # build Hartree potential pot = 1.0 / self.__rr ulmintp = self.calculate_interpolated_ulm(self.__rr).transpose() - self.__htpot = np.einsum('ij,jk,i,ik->ik', ulmintp, self.__ylm, pot, self.__edens) + self.__htpot = np.einsum('ij,jk,i,ik->ik', ulmintp, + self.__ylm, + pot, + self.__edens, + optimize=True) - return np.einsum('ij,ij', self.__htpot, self.__wgrid) + return np.einsum('ij,ij', self.__htpot, self.__wgrid, optimize=True) def __build_weight_grid(self): """ diff --git a/pydft/dft.py b/pydft/dft.py index c688c23..eaf6933 100644 --- a/pydft/dft.py +++ b/pydft/dft.py @@ -203,7 +203,7 @@ def scf(self, tol:float=1e-5, verbose:bool=False) -> float: self.__time_stats['iterations'].append(itertime) if verbose: - print('%03i | E = %12.6f | dE = %5.4e | %0.4f ms' % (niter+1, energy, ediff, itertime)) + print('%03i | E = %12.6f | dE = %5.4e | %0.4f s' % (niter+1, energy, ediff, itertime)) if niter > 0: ediff = np.abs(energy - self.__energies[-2]) diff --git a/pydft/moleculargrid.py b/pydft/moleculargrid.py index 7ad8e19..2eacdbc 100644 --- a/pydft/moleculargrid.py +++ b/pydft/moleculargrid.py @@ -87,13 +87,15 @@ def build_density(self, P:np.ndarray, normalize:bool=True): self.__densities = np.einsum('ijk,jl,ilk->ik', self.__amplitudes, P, - self.__amplitudes) + self.__amplitudes, + optimize=True) # also build the gradient of the density self.__gradients = 2.0 * np.einsum('ijk,jl,ilkm->ikm', self.__amplitudes, P, - self.__ampgrads) + self.__ampgrads, + optimize=True) # perform optional normalization if normalize: @@ -251,7 +253,8 @@ def count_electrons_from_rho_lm(self) -> float: return np.einsum('ijk,ijk,ik', rho_lm_gpts, ylm, - np.array([an.get_weights() for an in self.__atomgrids])) + np.array([an.get_weights() for an in self.__atomgrids]), + optimize=True) def calculate_dfa_nuclear_attraction_local(self) -> float: """ @@ -302,7 +305,10 @@ def calculate_dfa_coulomb(self) -> float: # for each grid point determine the Hartree potential from the # spherical harmonic coefficients with respect to each atomic center - self.__ugpts = np.einsum('ijk,ik,ijk->k', ulmgpts, self.__rigridpoints, self.__ylmgpts) + self.__ugpts = np.einsum('ijk,ik,ijk->k', ulmgpts, + self.__rigridpoints, + self.__ylmgpts, + optimize=True) return 0.5 * np.einsum('i,i,i', self.__ugpts, np.array([an.get_density() for an in self.__atomgrids]).flatten(), @@ -337,20 +343,18 @@ def calculate_coulombic_matrix(self) -> np.array: # for each grid point determine the Hartree potential from the # spherical harmonic coefficients with respect to each atomic center - self.__ugpts = np.einsum('ijk,ik,ijk->k', ulmgpts, self.__rigridpoints, self.__ylmgpts) + self.__ugpts = np.einsum('ijk,ik,ijk->k', ulmgpts, + self.__rigridpoints, + self.__ylmgpts, + optimize=True) # construct coulombic repulsion matrix by integrating the interaction # of the hartree potential with the basis function amplitudes - N = len(self.__basis) - J = np.zeros((N,N)) - for i in range(0, N): - for j in range(i, N): - J[i,j] = np.einsum('i,i,i,i', self.__ugpts, - self.__fullgrid_amplitudes[i,:], - self.__fullgrid_amplitudes[j,:], - self.__mgw) - if i != j: - J[j,i] = J[i,j] + J = np.einsum('k,ik,jk,k->ij', self.__ugpts, + self.__fullgrid_amplitudes, + self.__fullgrid_amplitudes, + self.__mgw, + optimize=True) return J @@ -400,15 +404,11 @@ def calculate_exchange(self) -> (np.ndarray,float): X = np.zeros((len(self.__basis), len(self.__basis))) # exchange parameters - for i in range(0, len(self.__basis)): - for j in range(i, len(self.__basis)): - X[i,j] = np.einsum('i,i,i,i', vfx, - self.__fullgrid_amplitudes[i,:], - self.__fullgrid_amplitudes[j,:], - self.__mgw) - - if i != j: - X[j,i] = X[i,j] + X = np.einsum('k,ik,jk,k->ij', vfx, + self.__fullgrid_amplitudes, + self.__fullgrid_amplitudes, + self.__mgw, + optimize=True) return X, ex @@ -435,15 +435,11 @@ def calculate_correlation(self) -> (np.ndarray, float): C = np.zeros((len(self.__basis), len(self.__basis))) # exchange parameters - for i in range(0, len(self.__basis)): - for j in range(i, len(self.__basis)): - C[i,j] = np.einsum('i,i,i,i', vfc, - self.__fullgrid_amplitudes[i,:], - self.__fullgrid_amplitudes[j,:], - self.__mgw) - - if i != j: - C[j,i] = C[i,j] + C = np.einsum('k,ik,jk,k->ij', vfc, + self.__fullgrid_amplitudes, + self.__fullgrid_amplitudes, + self.__mgw, + optimize=True) return C, ec @@ -474,7 +470,7 @@ def get_density_at_points(self, for i,cgf in enumerate(self.__basis): amps[i,:] = np.array([cgf.get_amp(p) for p in spoints]) - dens = np.einsum('ik,ij,jk->k', amps, P, amps) + dens = np.einsum('ik,ij,jk->k', amps, P, amps, optimize=True) return dens @@ -504,7 +500,7 @@ def get_amplitude_at_points(self, for i,cgf in enumerate(self.__basis): amps[i,:] = np.array([cgf.get_amp(p) for p in spoints]) - wfamp = np.einsum('ik,i->k', amps, c) + wfamp = np.einsum('ik,i->k', amps, c, optimize=True) return wfamp From 5ddb4efff96faa3467d47479b6d60707427873ed Mon Sep 17 00:00:00 2001 From: ifilot Date: Sat, 14 Dec 2024 22:22:31 +0100 Subject: [PATCH 5/5] Adding multiprocessing features --- pydft/dft.py | 9 ++- pydft/moleculargrid.py | 85 +++++++++++++++++++++------- tests/test_finite_difference_grid.py | 6 +- tests/test_parallel.py | 29 ++++++++++ 4 files changed, 106 insertions(+), 23 deletions(-) create mode 100644 tests/test_parallel.py diff --git a/pydft/dft.py b/pydft/dft.py index eaf6933..cc1b7f5 100644 --- a/pydft/dft.py +++ b/pydft/dft.py @@ -20,7 +20,8 @@ def __init__(self, nangpts:int = 110, lmax:int = 8, fdpts:int=7, - normalize:bool = True): + normalize:bool = True, + parallel:bool = True): """ Constructs the DFT class @@ -45,6 +46,8 @@ def __init__(self, density verbose : bool, optional whether to provide verbose output, by default False + parallel : bool, optional + whether to use multiprocessing features, by default False """ self.__mol = mol self.__integrator = PyQInt() @@ -57,6 +60,7 @@ def __init__(self, self.__lmax = lmax self.__functional = functional self.__normalize = normalize + self.__parallel = parallel # keep track of time self.calctimes = { @@ -352,7 +356,8 @@ def __setup(self): nangpts=self.__nangpts, lmax=self.__lmax, fdpts=self.__fdpts, - functional=self.__functional) + functional=self.__functional, + parallel=self.__parallel) self.__molgrid.initialize() # molecular grid uses late initialization # build one-electron matrices; because these matrices are Hermetian, diff --git a/pydft/moleculargrid.py b/pydft/moleculargrid.py index 2eacdbc..fe61fb9 100644 --- a/pydft/moleculargrid.py +++ b/pydft/moleculargrid.py @@ -7,7 +7,9 @@ from pyqint import PyQInt import pyqint as pq import time +import multiprocessing from .xcfunctionals import Functionals +from functools import partial class MolecularGrid: def __init__(self, @@ -17,7 +19,8 @@ def __init__(self, nangpts:int=110, lmax:int=8, fdpts:int=7, - functional:str='svwn5'): + functional:str='svwn5', + parallel:bool=False): """ Construct MolecularGrid @@ -37,6 +40,8 @@ def __init__(self, number of grid point in finite difference scheme, by default 7 functional : str, optional exchange-correlation functional, by default 'svwn5' + parallel : bool, optional + whether to use multiprocessing features, by default False """ # keep track of build times self.construct_times = {} @@ -51,6 +56,7 @@ def __init__(self, self.__basis = cgfs self.__functionals = Functionals(functional) self.__is_initialized = False + self.__enable_parallel = True def initialize(self): """ @@ -62,7 +68,7 @@ def initialize(self): if self.__is_initialized: return - self.__build_molecular_grid() + self.__build_molecular_grid(self.__enable_parallel) self.__build_amplitudes() self.__is_initialized = True @@ -344,17 +350,17 @@ def calculate_coulombic_matrix(self) -> np.array: # for each grid point determine the Hartree potential from the # spherical harmonic coefficients with respect to each atomic center self.__ugpts = np.einsum('ijk,ik,ijk->k', ulmgpts, - self.__rigridpoints, - self.__ylmgpts, - optimize=True) + self.__rigridpoints, + self.__ylmgpts, + optimize=True) # construct coulombic repulsion matrix by integrating the interaction # of the hartree potential with the basis function amplitudes J = np.einsum('k,ik,jk,k->ij', self.__ugpts, - self.__fullgrid_amplitudes, - self.__fullgrid_amplitudes, - self.__mgw, - optimize=True) + self.__fullgrid_amplitudes, + self.__fullgrid_amplitudes, + self.__mgw, + optimize=True) return J @@ -541,7 +547,7 @@ def get_gradient_at_points(self, #t2 = np.einsum('ij,jkl,ik->kl', P, grads, amps) #return t1 + t2 - return 2.0 * np.einsum('ij,ikl,jk->kl', P, grads, amps) + return 2.0 * np.einsum('ij,ikl,jk->kl', P, grads, amps, optimize=True) def calculate_coulomb_potential_at_points(self, pts:np.ndarray) -> np.ndarray: @@ -729,7 +735,7 @@ def calculate_weights_at_points(self, return mweights - def __build_molecular_grid(self): + def __build_molecular_grid(self, build_parallel=False): """ Build the molecular grid from the atomic grids """ @@ -841,20 +847,61 @@ def __build_molecular_grid(self): self.__ylmgpts = np.ndarray((len(self.__atoms), (self.__lmax+1)**2, np.prod(self.__mweights.shape))) - for i,at in enumerate(self.__atoms): # loop over atoms - lmctr = 0 - for l in range(0, self.__lmax+1): # loop over spherical harmonics - for m in range(-l, l+1): - self.__ylmgpts[i,lmctr,:] = spherical_harmonic(l, m, \ - self.__theta_gridpoints[i,:], - self.__phi_gridpoints[i,:]) - lmctr += 1 + + # perform parallellized calculation of spherical harmonics + if build_parallel: + atoms = list(range(len(self.__atoms))) + inputs = zip([self.__theta_gridpoints[i,:] for i in atoms], + [self.__phi_gridpoints[i,:] for i in atoms]) + calculate_partial = partial(self.build_ylmgpts_atom, + lmax=self.__lmax, + npts=np.prod(self.__mweights.shape)) + with multiprocessing.Pool(processes=multiprocessing.cpu_count()) as pool: + result = pool.map(calculate_partial, inputs) + for i,res in enumerate(result): + self.__ylmgpts[i,:,:] = res + else: # non-parallelized evaluation of spherical harmonics + for i,at in enumerate(self.__atoms): # loop over atoms + self.__ylmgpts[i,:,:] = self.__build_ylmgpts_atom(i) self.construct_times['spherical_harmonics'] = time.time() - st # collect complete weights for the full molecular grid weights = np.array([an.get_weights() for an in self.__atomgrids]).flatten() self.__mgw = np.multiply(weights, self.__mweights.flatten()) + @staticmethod + def build_ylmgpts_atom(Omega, lmax, npts): + """ + Get the spherical harmonic parameters for a single atom + + This function is 'pickable' and can be used in a multiprocessing pool + """ + lmctr = 0 + theta = Omega[0] + phi = Omega[1] + res = np.ndarray(((lmax+1)**2, npts)) + for l in range(0, lmax+1): # loop over spherical harmonics + for m in range(-l, l+1): + res[lmctr,:] = spherical_harmonic(l, m, \ + theta, + phi) + lmctr += 1 + return res + + def __build_ylmgpts_atom(self, atidx): + """ + Get the spherical harmonic parameters for a single atom + """ + lmctr = 0 + res = np.ndarray(((self.__lmax+1)**2, np.prod(self.__mweights.shape))) + for l in range(0, self.__lmax+1): # loop over spherical harmonics + for m in range(-l, l+1): + res[lmctr,:] = spherical_harmonic(l, m, \ + self.__theta_gridpoints[atidx,:], + self.__phi_gridpoints[atidx,:]) + lmctr += 1 + return res + def __build_amplitudes(self): """ Precalculate the basis function amplitudes at the selected grid diff --git a/tests/test_finite_difference_grid.py b/tests/test_finite_difference_grid.py index 29c07e2..159995d 100644 --- a/tests/test_finite_difference_grid.py +++ b/tests/test_finite_difference_grid.py @@ -8,11 +8,13 @@ from pydft import MoleculeBuilder, DFT -class TestDFT(unittest.TestCase): +class TestFiniteDifferenceSchemes(unittest.TestCase): def test_co(self): """ - Test DFT calculation of CO molecule + Test DFT calculation of CO molecule using various finite difference + schemes (always diagonal matrices, but with different number of + non-zero diagonals) """ mol_builder = MoleculeBuilder() mol = mol_builder.from_name('CO') diff --git a/tests/test_parallel.py b/tests/test_parallel.py new file mode 100644 index 0000000..cfe8e2f --- /dev/null +++ b/tests/test_parallel.py @@ -0,0 +1,29 @@ +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 TestParallel(unittest.TestCase): + + + def test_benzene(self): + """ + Test DFT calculation of benzene molecule + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('benzene') + + # construct dft object + dft = DFT(mol, basis='sto3g', parallel=True) + energy = dft.scf() + + answer = -228.06536332917378 + np.testing.assert_almost_equal(energy, answer, 4) + +if __name__ == '__main__': + unittest.main()