From a6a12ad4868ce4ab6a2d5b56a1c836114719e101 Mon Sep 17 00:00:00 2001 From: Wesley Sanchez Date: Mon, 29 Jan 2024 09:41:05 -0500 Subject: [PATCH 1/9] create new functions for evaluating density matrices on a grid --- gbasis/evals/density.py | 99 +++++++++++++++++++++++++++++++++++++++++ tests/test_density.py | 29 ++++++++++++ 2 files changed, 128 insertions(+) diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py index c0712924..8c6dfffd 100644 --- a/gbasis/evals/density.py +++ b/gbasis/evals/density.py @@ -102,6 +102,105 @@ def evaluate_density(one_density_matrix, basis, points, transform=None, threshol raise ValueError(f"Found negative density <= {-threshold}, got {min_output}.") return output.clip(min=0.0) +def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): + """Return the evaluation of the density matrix given the evaluated orbitals. + + Parameters + ---------- + one_density_matrix : np.ndarray(K_orb, K_orb) + One-electron density matrix in terms of the given orbitals. + orb_eval_list : list of orbitals evaluated on grid points [np.ndarray(K_orb, N) ...] + Orbitals evaluated at :math:`N` grid points. + The set of orbitals must be the same basis set used to build the one-electron density + matrix. + + Returns + ------- + density : np.ndarray(N1,N2,K_orb,K_orb) + Density Matrix evaluated at `N1,N2` grid points. + + Raises + ------ + TypeError + If `orb_eval` is not a 2-dimensional `numpy` array with `dtype` float. + If `one_density_matrix` is not a 2-dimensional `numpy` array with `dtype` float. + + +""" + if not ( + isinstance(one_density_matrix, np.ndarray) + and one_density_matrix.ndim == 2 + and one_density_matrix.dtype == float + ): + raise TypeError( + "One-electron density matrix must be a two-dimensional `numpy` array with `dtype`" + " float." + ) + for orb in orb_eval_list: + if not (isinstance(orb, np.ndarray) and orb.ndim == 2 and orb.dtype == float): + raise TypeError( + "Evaluation of orbitals must be a two-dimensional `numpy` array with `dtype` float." + ) + if one_density_matrix.shape[0] != one_density_matrix.shape[1]: + raise ValueError("One-electron density matrix must be a square matrix.") + + if not np.allclose(one_density_matrix, one_density_matrix.T): + raise ValueError("One-electron density matrix must be symmetric.") + + + #Tensor product for \gamma(\mathbf{r}_1,\mathbf{r}_2) = \sum_{pq} \gamma_{pq} \chi_p(\mathbf{r}_1) \chi_q(\mathbf{r}_2) + tensor_product = np.einsum('ij,ik,jl->klij',one_density_matrix, orb_eval_list[0],orb_eval_list[1]) + + #returns dm evaluated on each grid point + return tensor_product + + +def evaluate_dm_density(one_density_matrix, basis, points_list, transform=None): + r"""Return the density of the given basis set at the given points. + + Parameters + ---------- + one_density_matrix : np.ndarray(K_orbs, K_orbs) + One-electron density matrix in terms of the given basis set. + If the basis is transformed using `transform` keyword, then the density matrix is assumed to + be expressed with respect to the transformed basis set. + basis : list/tuple of GeneralizedContractionShell + Shells of generalized contractions. + points_list : list of points [np.ndarray(N, 3)...] + Cartesian coordinates of the points in space (in atomic units) where the basis functions + are evaluated. + Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` + components. + This function can take a list of points at which basis functions are evaluated. If only one + set of points is given, it will be duplicated. + transform : np.ndarray(K_orbs, K_cont) + Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear + combinations of contractions (e.g. MO). + Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` + and index 0 of the array for contractions. + Default is no transformation. + + + Returns + ------- + dm_on_grid : np.ndarray(N1,N2,K_orb,K_orb) + Density Matrix evaluated at `N1,N2` grid points. + + """ + + orb_evals = [] + #evaluate basi(e)s on the point set(s) + for grid in points_list: + orb_eval = evaluate_basis(basis, grid, transform=transform) + orb_evals.append(orb_eval) + #if only one set of points is supplied, it is duplicated + if len(points_list)==1: + orb_evals.append(orb_eval) + + #Calulated performed using the evaluated orbitals + dm_on_grid = evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_evals) + + return dm_on_grid def evaluate_deriv_reduced_density_matrix( orders_one, diff --git a/tests/test_density.py b/tests/test_density.py index 0b12eeda..03147775 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -1,6 +1,8 @@ """Test gbasis.evals.density.""" from gbasis.evals.density import ( evaluate_density, + evaluate_dm_using_evaluated_orbs, + evaluate_dm_density, evaluate_density_gradient, evaluate_density_hessian, evaluate_density_laplacian, @@ -73,6 +75,32 @@ def test_evaluate_density(): assert np.all(dens >= 0.0) assert np.allclose(dens, np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) +def test_evaluate_dm_using_evaluated_orbs(): + """Test gbasis.evals.density.evaluate_density_using_evaluated_orbs.""" + density_mat = np.array([[1.0, 2.0], [2.0, 3.0]]) + orb_eval = np.array([[1.0], [2.0]]) + dens = evaluate_dm_using_evaluated_orbs(density_mat, [orb_eval,orb_eval]) + assert np.all(dens >= 0.0) + assert np.allclose(np.einsum('iikl->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) + + density_mat = np.array([[1.0, 2.0], [2.0, 3.0]]) + orb_eval = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) + dens = evaluate_dm_using_evaluated_orbs(density_mat, [orb_eval,orb_eval]) + assert np.all(dens >= 0) + assert np.allclose(np.einsum('iikl->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) + +def test_evaluate_dm_density(): + """Test gbasis.evals.density.evaluate_density.""" + basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) + basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") + transform = np.random.rand(14, 18) + density = np.random.rand(14, 14) + density += density.T + points = np.random.rand(100, 3) + + evaluate_orbs = evaluate_basis(basis, points, transform) + dens = evaluate_dm_density(density, basis, [points], transform) + assert np.allclose(np.einsum('iikl->i',dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) def test_evaluate_deriv_density(): """Test gbasis.evals.density.evaluate_deriv_density.""" @@ -448,3 +476,4 @@ def test_evaluate_general_kinetic_energy_density(): evaluate_posdef_kinetic_energy_density(np.identity(40), basis, points, np.identity(40)) + evaluate_density_laplacian(np.identity(40), basis, points, np.identity(40)), ) + From e0b738ae2cb18dd6b039723b7b42237626602363 Mon Sep 17 00:00:00 2001 From: Wesley Sanchez Date: Mon, 29 Jan 2024 09:58:00 -0500 Subject: [PATCH 2/9] create function for evaluating two-body hole correlation function on a grid --- gbasis/evals/density.py | 60 ++++++++++++++++++++++++++++++++++- tests/data_rdm_h2o_sto3g.npy | Bin 0 -> 520 bytes tests/data_rdm_kr_sto6g.npy | Bin 0 -> 2720 bytes tests/test_density.py | 17 ++++++++++ 4 files changed, 76 insertions(+), 1 deletion(-) create mode 100644 tests/data_rdm_h2o_sto3g.npy create mode 100644 tests/data_rdm_kr_sto6g.npy diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py index 8c6dfffd..d24d2ccb 100644 --- a/gbasis/evals/density.py +++ b/gbasis/evals/density.py @@ -147,7 +147,7 @@ def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): if not np.allclose(one_density_matrix, one_density_matrix.T): raise ValueError("One-electron density matrix must be symmetric.") - + #Tensor product for \gamma(\mathbf{r}_1,\mathbf{r}_2) = \sum_{pq} \gamma_{pq} \chi_p(\mathbf{r}_1) \chi_q(\mathbf{r}_2) tensor_product = np.einsum('ij,ik,jl->klij',one_density_matrix, orb_eval_list[0],orb_eval_list[1]) @@ -202,6 +202,64 @@ def evaluate_dm_density(one_density_matrix, basis, points_list, transform=None): return dm_on_grid +def evaluate_hole_x2(one_density_matrix, basis, points_list, transform=None): + r"""Return the two-body hole correlation function. + + .. math :: + + \left. + h(\mathbf{r}_1,\mathbf{r}_2) = + \right. + -\frac{|\gamma(\mathbf{r}_1,\mathbf{r}_2)|^2} + {\rho({\mathbf{r}_1})\rho({\mathbf{r}_2})} + + Parameters + ---------- + one_density_matrix : np.ndarray(K_orbs, K_orbs) + One-electron density matrix in terms of the given basis set. + If the basis is transformed using `transform` keyword, then the density matrix is assumed to + be expressed with respect to the transformed basis set. + basis : list/tuple of GeneralizedContractionShell + Shells of generalized contractions. + points_list : list of points [np.ndarray(N, 3)] + Cartesian coordinates of the points in space (in atomic units) where the basis functions + are evaluated. + Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` + components. + This function can take a list of points at which basis functions are evaluated. If only one + set of points is given, it will be duplicated. + transform : np.ndarray(K_orbs, K_cont) + Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear + combinations of contractions (e.g. MO). + Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` + and index 0 of the array for contractions. + Default is no transformation. + + + Returns + ------- + hole_x2 : np.ndarray(N1,N2,K_orb,K_orb) + Two-body Exchange Hole evaluated at `N1,N2` grid points. + + """ + dens_evals = [] + + for grid in points_list: + dens_eval = evaluate_density(one_density_matrix,basis, grid, transform=transform) + dens_evals.append(dens_eval) + + if len(points_list)==1: + dens_evals.append(dens_eval) + + #build density matrix on grid + dm_eval = evaluate_dm_density(one_density_matrix, basis, points_list, transform=transform) + + #evaluate hole function + numerator = np.einsum('ijkl,jikl->ijkl',dm_eval,dm_eval) + hole_x2 = -1*np.einsum('ijkl,i,j->ijkl',numerator,1/dens_evals[0],1/dens_evals[1]) + + return hole_x2 + def evaluate_deriv_reduced_density_matrix( orders_one, orders_two, diff --git a/tests/data_rdm_h2o_sto3g.npy b/tests/data_rdm_h2o_sto3g.npy new file mode 100644 index 0000000000000000000000000000000000000000..0e2b60e593483d87f43ab0824aeca3dfef7d35e3 GIT binary patch literal 520 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(4=E~51qv5u zBo?Fsxf4R5Z4NzYJj1?U zaoNhatW_I#Z*ovCG=49n>s^WeYzgd2y8zLn3}lf-Io8t&X? zf7YBYRaD1zkI2+bM+Df*?Vs7^PmoOCxd&!GO#O)~e>P~A|F>VXUfi!=W%-^xrVIHW zu;uPaNo|R!Z@6h=5ctBkE2(-v%snvkVd_()G?k_2|F@6TJ21&zA=8G-)vror(Nz1i zZMMrdUAnXf=6;xaVCKWrZ=cg(?)&q<{j~Z2Ow{f6?74UAM$tRvvbEtPL!DVE%=e2Ky8eXJ@0BNd6(-(>11po#xAMyYI literal 0 HcmV?d00001 diff --git a/tests/test_density.py b/tests/test_density.py index 03147775..4fe47ca3 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -3,6 +3,7 @@ evaluate_density, evaluate_dm_using_evaluated_orbs, evaluate_dm_density, + evaluate_hole_x2, evaluate_density_gradient, evaluate_density_hessian, evaluate_density_laplacian, @@ -102,6 +103,21 @@ def test_evaluate_dm_density(): dens = evaluate_dm_density(density, basis, [points], transform) assert np.allclose(np.einsum('iikl->i',dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) + +def test_evaluate_hole_x2(): + basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) + basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") + transform = np.ones((18, 18)) + density = np.load(find_datafile("data_rdm_kr_sto6g.npy")) + + points = np.random.rand(10, 3) + + eh = evaluate_hole_x2(density, basis, points_list=[points,points],transform=transform) + ed = evaluate_density(density,basis,points,transform=transform) + + assert np.allclose(np.einsum('ijkl,j->',eh,ed)/np.sum(ed)/len(points)*len(transform), -1) + + def test_evaluate_deriv_density(): """Test gbasis.evals.density.evaluate_deriv_density.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) @@ -477,3 +493,4 @@ def test_evaluate_general_kinetic_energy_density(): + evaluate_density_laplacian(np.identity(40), basis, points, np.identity(40)), ) +test_evaluate_hole_x2() \ No newline at end of file From 21cb294c468552cf10771db82d675b8c5109aecf Mon Sep 17 00:00:00 2001 From: Wesley Sanchez Date: Mon, 29 Jan 2024 10:26:10 -0500 Subject: [PATCH 3/9] fixed various instances of docustring var_names not matching var_names in return statements in density.py --- gbasis/evals/density.py | 58 ++++++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py index d24d2ccb..88d3bd33 100644 --- a/gbasis/evals/density.py +++ b/gbasis/evals/density.py @@ -92,15 +92,15 @@ def evaluate_density(one_density_matrix, basis, points, transform=None, threshol ------- density : np.ndarray(N,) Density evaluated at `N` grid points. - + """ orb_eval = evaluate_basis(basis, points, transform=transform) - output = evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval) + density = evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval) # Fix #117: check magnitude of small negative density values, then use clip to remove them - min_output = np.min(output) - if min_output < 0.0 and abs(min_output) > threshold: - raise ValueError(f"Found negative density <= {-threshold}, got {min_output}.") - return output.clip(min=0.0) + min_density = np.min(density) + if min_density < 0.0 and abs(min_density) > threshold: + raise ValueError(f"Found negative density <= {-threshold}, got {min_density}.") + return density.clip(min=0.0) def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): """Return the evaluation of the density matrix given the evaluated orbitals. @@ -122,8 +122,8 @@ def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): Raises ------ TypeError - If `orb_eval` is not a 2-dimensional `numpy` array with `dtype` float. - If `one_density_matrix` is not a 2-dimensional `numpy` array with `dtype` float. + If `orb_eval_list` does not consist of 2-dimensional `numpy` arrays with `dtype` float. + If `one_density_matrix` is not a symmetric, 2-dimensional `numpy` array with `dtype` float. """ @@ -149,10 +149,10 @@ def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): #Tensor product for \gamma(\mathbf{r}_1,\mathbf{r}_2) = \sum_{pq} \gamma_{pq} \chi_p(\mathbf{r}_1) \chi_q(\mathbf{r}_2) - tensor_product = np.einsum('ij,ik,jl->klij',one_density_matrix, orb_eval_list[0],orb_eval_list[1]) + density = np.einsum('ij,ik,jl->klij',one_density_matrix, orb_eval_list[0],orb_eval_list[1]) #returns dm evaluated on each grid point - return tensor_product + return density def evaluate_dm_density(one_density_matrix, basis, points_list, transform=None): @@ -345,8 +345,8 @@ def evaluate_deriv_reduced_density_matrix( ) density = one_density_matrix.dot(deriv_orb_eval_two) density *= deriv_orb_eval_one - density = np.sum(density, axis=0) - return density + deriv_reduced_density_matrix = np.sum(density, axis=0) + return deriv_reduced_density_matrix def evaluate_deriv_density( @@ -395,7 +395,7 @@ def evaluate_deriv_density( # pylint: disable=R0914 total_l_x, total_l_y, total_l_z = orders - output = np.zeros(points.shape[0]) + density_deriv = np.zeros(points.shape[0]) for l_x in range(total_l_x // 2 + 1): # prevent double counting for the middle of the even total_l_x # e.g. If total_l_x == 4, then l_x is in [0, 1, 2, 3, 4]. Exploiting symmetry we only need @@ -430,8 +430,8 @@ def evaluate_deriv_density( transform=transform, deriv_type=deriv_type, ) - output += factor * num_occurence * density - return output + density_deriv += factor * num_occurence * density + return density_deriv def evaluate_density_gradient( @@ -475,7 +475,7 @@ def evaluate_density_gradient( """ orders_one = np.array(([1, 0, 0], [0, 1, 0], [0, 0, 1])) - output = np.zeros((3, len(points))) + density_gradient = np.zeros((3, len(points))) # Evaluation of generalized contraction shell for zeroth order = 0,0,0 zeroth_deriv = evaluate_deriv_basis( basis, @@ -497,8 +497,8 @@ def evaluate_density_gradient( # output[ind] = 2*(np.einsum('ij,ik,jk -> k',one_density_matrix, zeroth_deriv, deriv_comp)) density = one_density_matrix.dot(zeroth_deriv) density *= deriv_comp - output[ind] = 2 * 1 * np.sum(density, axis=0) - return output.T + density_gradient[ind] = 2 * 1 * np.sum(density, axis=0) + return density_gradient.T def evaluate_density_laplacian( @@ -544,7 +544,7 @@ def evaluate_density_laplacian( orders_one_second = np.array(([2, 0, 0], [0, 2, 0], [0, 0, 2])) orders_one_first = np.array(([1, 0, 0], [0, 1, 0], [0, 0, 1])) orders_two = np.array(([1, 0, 0], [0, 1, 0], [0, 0, 1])) - output = np.zeros(points.shape[0]) + density_laplacian = np.zeros(points.shape[0]) # Evaluation of generalized contraction shell for zeroth order = 0,0,0 zeroth_deriv = evaluate_deriv_basis( basis, @@ -566,7 +566,7 @@ def evaluate_density_laplacian( density = one_density_matrix.dot(zeroth_deriv) density *= deriv_one - output += 2 * 1 * np.sum(density, axis=0) + density_laplacian += 2 * 1 * np.sum(density, axis=0) for orders in zip(orders_one_first, orders_two): deriv_one = evaluate_deriv_basis( @@ -587,9 +587,9 @@ def evaluate_density_laplacian( # output[ind] = 2*(np.einsum('ij,ik,jk -> k',one_density_matrix, zeroth_deriv, deriv_comp)) density = one_density_matrix.dot(deriv_two) density *= deriv_one - output += 2 * 1 * np.sum(density, axis=0) + density_laplacian += 2 * 1 * np.sum(density, axis=0) - return output + return density_laplacian def evaluate_density_hessian( @@ -702,13 +702,13 @@ def evaluate_density_hessian( density_2 = np.einsum("ijkm,ijmk -> ijkm", one_two_arr_1, raw_density_2) # factors and sum over basis functions - output = 2 * 1 * np.sum(density_1, axis=2) - output += 2 * 1 * np.sum(density_2, axis=2) + density_hessian = 2 * 1 * np.sum(density_1, axis=2) + density_hessian += 2 * 1 * np.sum(density_2, axis=2) # copying lower matrix to upper matrix - upp = np.swapaxes(output, 0, 1) + upp = np.swapaxes(density_hessian, 0, 1) upp = np.triu(upp.T, 1) - return output.T + upp + return density_hessian.T + upp def evaluate_posdef_kinetic_energy_density( @@ -767,9 +767,9 @@ def evaluate_posdef_kinetic_energy_density( `N` grid points. """ - output = np.zeros(points.shape[0]) + posdef_kindetic_energy_density = np.zeros(points.shape[0]) for orders in np.identity(3, dtype=int): - output += evaluate_deriv_reduced_density_matrix( + posdef_kindetic_energy_density += evaluate_deriv_reduced_density_matrix( orders, orders, one_density_matrix, @@ -782,7 +782,7 @@ def evaluate_posdef_kinetic_energy_density( min_output = np.min(output) if min_output < 0.0 and abs(min_output) > threshold: raise ValueError(f"Found negative density <= {-threshold}, got {min_output}.") - return (0.5 * output).clip(min=0.0) + return (0.5 * posdef_kindetic_energy_density).clip(min=0.0) # TODO: test against a reference From 06bfe3ad82a4bbcca15b746d93b75d1df5dd34ec Mon Sep 17 00:00:00 2001 From: Wesley Sanchez Date: Mon, 29 Jan 2024 12:54:00 -0500 Subject: [PATCH 4/9] fix typo --- gbasis/evals/density.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py index 88d3bd33..50c3b7ee 100644 --- a/gbasis/evals/density.py +++ b/gbasis/evals/density.py @@ -255,7 +255,7 @@ def evaluate_hole_x2(one_density_matrix, basis, points_list, transform=None): dm_eval = evaluate_dm_density(one_density_matrix, basis, points_list, transform=transform) #evaluate hole function - numerator = np.einsum('ijkl,jikl->ijkl',dm_eval,dm_eval) + numerator = np.einsum('ijkl,ijkl->ijkl',dm_eval,dm_eval) hole_x2 = -1*np.einsum('ijkl,i,j->ijkl',numerator,1/dens_evals[0],1/dens_evals[1]) return hole_x2 From 96cd22699d34b1375dd8a60080612d70d0aaa082 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 29 Jan 2024 21:36:37 +0000 Subject: [PATCH 5/9] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- gbasis/evals/density.py | 43 +++++++++++++++++++++-------------------- tests/test_density.py | 29 ++++++++++++++++++--------- 2 files changed, 42 insertions(+), 30 deletions(-) diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py index 50c3b7ee..a086f08f 100644 --- a/gbasis/evals/density.py +++ b/gbasis/evals/density.py @@ -92,7 +92,7 @@ def evaluate_density(one_density_matrix, basis, points, transform=None, threshol ------- density : np.ndarray(N,) Density evaluated at `N` grid points. - + """ orb_eval = evaluate_basis(basis, points, transform=transform) density = evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval) @@ -102,6 +102,7 @@ def evaluate_density(one_density_matrix, basis, points, transform=None, threshol raise ValueError(f"Found negative density <= {-threshold}, got {min_density}.") return density.clip(min=0.0) + def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): """Return the evaluation of the density matrix given the evaluated orbitals. @@ -125,8 +126,7 @@ def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): If `orb_eval_list` does not consist of 2-dimensional `numpy` arrays with `dtype` float. If `one_density_matrix` is not a symmetric, 2-dimensional `numpy` array with `dtype` float. - -""" + """ if not ( isinstance(one_density_matrix, np.ndarray) and one_density_matrix.ndim == 2 @@ -143,15 +143,14 @@ def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): ) if one_density_matrix.shape[0] != one_density_matrix.shape[1]: raise ValueError("One-electron density matrix must be a square matrix.") - + if not np.allclose(one_density_matrix, one_density_matrix.T): raise ValueError("One-electron density matrix must be symmetric.") - - #Tensor product for \gamma(\mathbf{r}_1,\mathbf{r}_2) = \sum_{pq} \gamma_{pq} \chi_p(\mathbf{r}_1) \chi_q(\mathbf{r}_2) - density = np.einsum('ij,ik,jl->klij',one_density_matrix, orb_eval_list[0],orb_eval_list[1]) + # Tensor product for \gamma(\mathbf{r}_1,\mathbf{r}_2) = \sum_{pq} \gamma_{pq} \chi_p(\mathbf{r}_1) \chi_q(\mathbf{r}_2) + density = np.einsum("ij,ik,jl->klij", one_density_matrix, orb_eval_list[0], orb_eval_list[1]) - #returns dm evaluated on each grid point + # returns dm evaluated on each grid point return density @@ -171,7 +170,7 @@ def evaluate_dm_density(one_density_matrix, basis, points_list, transform=None): are evaluated. Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` components. - This function can take a list of points at which basis functions are evaluated. If only one + This function can take a list of points at which basis functions are evaluated. If only one set of points is given, it will be duplicated. transform : np.ndarray(K_orbs, K_cont) Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear @@ -189,26 +188,27 @@ def evaluate_dm_density(one_density_matrix, basis, points_list, transform=None): """ orb_evals = [] - #evaluate basi(e)s on the point set(s) + # evaluate basi(e)s on the point set(s) for grid in points_list: orb_eval = evaluate_basis(basis, grid, transform=transform) orb_evals.append(orb_eval) - #if only one set of points is supplied, it is duplicated - if len(points_list)==1: + # if only one set of points is supplied, it is duplicated + if len(points_list) == 1: orb_evals.append(orb_eval) - #Calulated performed using the evaluated orbitals + # Calulated performed using the evaluated orbitals dm_on_grid = evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_evals) return dm_on_grid + def evaluate_hole_x2(one_density_matrix, basis, points_list, transform=None): r"""Return the two-body hole correlation function. .. math :: \left. - h(\mathbf{r}_1,\mathbf{r}_2) = + h(\mathbf{r}_1,\mathbf{r}_2) = \right. -\frac{|\gamma(\mathbf{r}_1,\mathbf{r}_2)|^2} {\rho({\mathbf{r}_1})\rho({\mathbf{r}_2})} @@ -226,7 +226,7 @@ def evaluate_hole_x2(one_density_matrix, basis, points_list, transform=None): are evaluated. Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` components. - This function can take a list of points at which basis functions are evaluated. If only one + This function can take a list of points at which basis functions are evaluated. If only one set of points is given, it will be duplicated. transform : np.ndarray(K_orbs, K_cont) Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear @@ -245,21 +245,22 @@ def evaluate_hole_x2(one_density_matrix, basis, points_list, transform=None): dens_evals = [] for grid in points_list: - dens_eval = evaluate_density(one_density_matrix,basis, grid, transform=transform) + dens_eval = evaluate_density(one_density_matrix, basis, grid, transform=transform) dens_evals.append(dens_eval) - if len(points_list)==1: + if len(points_list) == 1: dens_evals.append(dens_eval) - #build density matrix on grid + # build density matrix on grid dm_eval = evaluate_dm_density(one_density_matrix, basis, points_list, transform=transform) - #evaluate hole function - numerator = np.einsum('ijkl,ijkl->ijkl',dm_eval,dm_eval) - hole_x2 = -1*np.einsum('ijkl,i,j->ijkl',numerator,1/dens_evals[0],1/dens_evals[1]) + # evaluate hole function + numerator = np.einsum("ijkl,ijkl->ijkl", dm_eval, dm_eval) + hole_x2 = -1 * np.einsum("ijkl,i,j->ijkl", numerator, 1 / dens_evals[0], 1 / dens_evals[1]) return hole_x2 + def evaluate_deriv_reduced_density_matrix( orders_one, orders_two, diff --git a/tests/test_density.py b/tests/test_density.py index 4fe47ca3..03e6f5d6 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -76,19 +76,25 @@ def test_evaluate_density(): assert np.all(dens >= 0.0) assert np.allclose(dens, np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) + def test_evaluate_dm_using_evaluated_orbs(): """Test gbasis.evals.density.evaluate_density_using_evaluated_orbs.""" density_mat = np.array([[1.0, 2.0], [2.0, 3.0]]) orb_eval = np.array([[1.0], [2.0]]) - dens = evaluate_dm_using_evaluated_orbs(density_mat, [orb_eval,orb_eval]) + dens = evaluate_dm_using_evaluated_orbs(density_mat, [orb_eval, orb_eval]) assert np.all(dens >= 0.0) - assert np.allclose(np.einsum('iikl->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) + assert np.allclose( + np.einsum("iikl->i", dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval) + ) density_mat = np.array([[1.0, 2.0], [2.0, 3.0]]) orb_eval = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) - dens = evaluate_dm_using_evaluated_orbs(density_mat, [orb_eval,orb_eval]) + dens = evaluate_dm_using_evaluated_orbs(density_mat, [orb_eval, orb_eval]) assert np.all(dens >= 0) - assert np.allclose(np.einsum('iikl->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) + assert np.allclose( + np.einsum("iikl->i", dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval) + ) + def test_evaluate_dm_density(): """Test gbasis.evals.density.evaluate_density.""" @@ -101,7 +107,9 @@ def test_evaluate_dm_density(): evaluate_orbs = evaluate_basis(basis, points, transform) dens = evaluate_dm_density(density, basis, [points], transform) - assert np.allclose(np.einsum('iikl->i',dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) + assert np.allclose( + np.einsum("iikl->i", dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs) + ) def test_evaluate_hole_x2(): @@ -112,10 +120,12 @@ def test_evaluate_hole_x2(): points = np.random.rand(10, 3) - eh = evaluate_hole_x2(density, basis, points_list=[points,points],transform=transform) - ed = evaluate_density(density,basis,points,transform=transform) + eh = evaluate_hole_x2(density, basis, points_list=[points, points], transform=transform) + ed = evaluate_density(density, basis, points, transform=transform) - assert np.allclose(np.einsum('ijkl,j->',eh,ed)/np.sum(ed)/len(points)*len(transform), -1) + assert np.allclose( + np.einsum("ijkl,j->", eh, ed) / np.sum(ed) / len(points) * len(transform), -1 + ) def test_evaluate_deriv_density(): @@ -493,4 +503,5 @@ def test_evaluate_general_kinetic_energy_density(): + evaluate_density_laplacian(np.identity(40), basis, points, np.identity(40)), ) -test_evaluate_hole_x2() \ No newline at end of file + +test_evaluate_hole_x2() From 59c34a5dde57b49ab779bba0029b92e42d4a8340 Mon Sep 17 00:00:00 2001 From: Wesley Sanchez Date: Thu, 1 Feb 2024 10:37:06 -0500 Subject: [PATCH 6/9] code changes+update functions --- gbasis/evals/density.py | 116 +++++++++++++++++++++++++--------------- tests/test_density.py | 60 ++++++++++++++++----- 2 files changed, 121 insertions(+), 55 deletions(-) diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py index 50c3b7ee..9eb85dc6 100644 --- a/gbasis/evals/density.py +++ b/gbasis/evals/density.py @@ -102,21 +102,25 @@ def evaluate_density(one_density_matrix, basis, points, transform=None, threshol raise ValueError(f"Found negative density <= {-threshold}, got {min_density}.") return density.clip(min=0.0) -def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): +def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_1,orb_eval_2): """Return the evaluation of the density matrix given the evaluated orbitals. Parameters ---------- one_density_matrix : np.ndarray(K_orb, K_orb) One-electron density matrix in terms of the given orbitals. - orb_eval_list : list of orbitals evaluated on grid points [np.ndarray(K_orb, N) ...] + orb_eval_1 : orbital evaluated on grid points np.ndarray(K_orb, N) + Orbitals evaluated at :math:`N` grid points. + The set of orbitals must be the same basis set used to build the one-electron density + matrix. + orb_eval_2 : orbital evaluated on grid points np.ndarray(K_orb, N) Orbitals evaluated at :math:`N` grid points. The set of orbitals must be the same basis set used to build the one-electron density matrix. Returns ------- - density : np.ndarray(N1,N2,K_orb,K_orb) + dm_on_grid : np.ndarray(N1,N2) Density Matrix evaluated at `N1,N2` grid points. Raises @@ -136,11 +140,15 @@ def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): "One-electron density matrix must be a two-dimensional `numpy` array with `dtype`" " float." ) - for orb in orb_eval_list: - if not (isinstance(orb, np.ndarray) and orb.ndim == 2 and orb.dtype == float): - raise TypeError( - "Evaluation of orbitals must be a two-dimensional `numpy` array with `dtype` float." - ) + + if not (isinstance(orb_eval_1, np.ndarray) and orb_eval_1.ndim == 2 and orb_eval_1.dtype == float): + raise TypeError( + "Evaluation of orbitals must be a two-dimensional `numpy` array with `dtype` float." + ) + if not (isinstance(orb_eval_2, np.ndarray) and orb_eval_2.ndim == 2 and orb_eval_2.dtype == float): + raise TypeError( + "Evaluation of orbitals must be a two-dimensional `numpy` array with `dtype` float." + ) if one_density_matrix.shape[0] != one_density_matrix.shape[1]: raise ValueError("One-electron density matrix must be a square matrix.") @@ -148,15 +156,18 @@ def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): raise ValueError("One-electron density matrix must be symmetric.") - #Tensor product for \gamma(\mathbf{r}_1,\mathbf{r}_2) = \sum_{pq} \gamma_{pq} \chi_p(\mathbf{r}_1) \chi_q(\mathbf{r}_2) - density = np.einsum('ij,ik,jl->klij',one_density_matrix, orb_eval_list[0],orb_eval_list[1]) + #Evaluation of \gamma(\mathbf{r}_1,\mathbf{r}_2) = \sum_{pq} \gamma_{pq} \chi_p(\mathbf{r}_1) \chi_q(\mathbf{r}_2) + dm_on_grid=0 + for ii, orb1 in enumerate(one_density_matrix): + for jj, orb2 in enumerate(one_density_matrix): + dm_on_grid += one_density_matrix[ii][jj]* np.outer(orb_eval_1[ii],orb_eval_2[jj]) - #returns dm evaluated on each grid point - return density + #returns dm evaluated on each grid point summed over the orbitals + return dm_on_grid -def evaluate_dm_density(one_density_matrix, basis, points_list, transform=None): - r"""Return the density of the given basis set at the given points. +def evaluate_dm_density(one_density_matrix, basis, point1, point2 = np.array(None), transform=None): + r"""Return the density matrix of the given basis set evaluated on the given points. Parameters ---------- @@ -166,7 +177,14 @@ def evaluate_dm_density(one_density_matrix, basis, points_list, transform=None): be expressed with respect to the transformed basis set. basis : list/tuple of GeneralizedContractionShell Shells of generalized contractions. - points_list : list of points [np.ndarray(N, 3)...] + point1 : set of points np.ndarray(N, 3) + Cartesian coordinates of the points in space (in atomic units) where the basis functions + are evaluated. + Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` + components. + This function can take a list of points at which basis functions are evaluated. If only one + set of points is given, it will be duplicated. + point2 : optional set of points np.ndarray(N, 3) --> default None Cartesian coordinates of the points in space (in atomic units) where the basis functions are evaluated. Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` @@ -183,26 +201,31 @@ def evaluate_dm_density(one_density_matrix, basis, points_list, transform=None): Returns ------- - dm_on_grid : np.ndarray(N1,N2,K_orb,K_orb) - Density Matrix evaluated at `N1,N2` grid points. + dm_on_grid : np.ndarray(N1,N2) + Density Matrix evaluated at `N1,N2` grid points summed over orbitals. """ - orb_evals = [] #evaluate basi(e)s on the point set(s) - for grid in points_list: - orb_eval = evaluate_basis(basis, grid, transform=transform) - orb_evals.append(orb_eval) - #if only one set of points is supplied, it is duplicated - if len(points_list)==1: - orb_evals.append(orb_eval) + orb_eval_1 = evaluate_basis(basis, point1, transform=transform) + #if only one set of points is supplied, it is duplicated + if point2.all() == None: + orb_eval_2 = orb_eval_1 + else: + orb_eval_2 = evaluate_basis(basis, point2, transform=transform) + + #Calulated performed using the evaluated orbitals - dm_on_grid = evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_evals) + dm_on_grid = evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_1,orb_eval_2) return dm_on_grid -def evaluate_hole_x2(one_density_matrix, basis, points_list, transform=None): + + + + +def evaluate_hole_x2(one_density_matrix, basis, point1, point2=np.array(None), transform=None): r"""Return the two-body hole correlation function. .. math :: @@ -221,7 +244,14 @@ def evaluate_hole_x2(one_density_matrix, basis, points_list, transform=None): be expressed with respect to the transformed basis set. basis : list/tuple of GeneralizedContractionShell Shells of generalized contractions. - points_list : list of points [np.ndarray(N, 3)] + point1 : set of points np.ndarray(N, 3) + Cartesian coordinates of the points in space (in atomic units) where the basis functions + are evaluated. + Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` + components. + This function can take a list of points at which basis functions are evaluated. If only one + set of points is given, it will be duplicated. + point2 : optional set of points np.ndarray(N, 3) --> default None Cartesian coordinates of the points in space (in atomic units) where the basis functions are evaluated. Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` @@ -238,25 +268,27 @@ def evaluate_hole_x2(one_density_matrix, basis, points_list, transform=None): Returns ------- - hole_x2 : np.ndarray(N1,N2,K_orb,K_orb) + hole_x2 : np.ndarray(N1,N2) Two-body Exchange Hole evaluated at `N1,N2` grid points. """ - dens_evals = [] - for grid in points_list: - dens_eval = evaluate_density(one_density_matrix,basis, grid, transform=transform) - dens_evals.append(dens_eval) + dens_eval_1 = evaluate_density(one_density_matrix,basis, point1, transform=transform) - if len(points_list)==1: - dens_evals.append(dens_eval) + #if only one set of points is supplied, density evaluation is duplicated + if point2.all() == None: + dens_eval_2 = dens_eval_1 + else: + dens_eval_2 = evaluate_density(one_density_matrix,basis, point2, transform=transform) #build density matrix on grid - dm_eval = evaluate_dm_density(one_density_matrix, basis, points_list, transform=transform) + dm_eval = evaluate_dm_density(one_density_matrix, basis, point1, point2, transform=transform) #evaluate hole function - numerator = np.einsum('ijkl,ijkl->ijkl',dm_eval,dm_eval) - hole_x2 = -1*np.einsum('ijkl,i,j->ijkl',numerator,1/dens_evals[0],1/dens_evals[1]) + numerator = dm_eval*dm_eval + + hole_x2 = -1*np.einsum('ij,i,j->ij',numerator,1/dens_eval_1,1/dens_eval_2) + return hole_x2 @@ -762,14 +794,14 @@ def evaluate_posdef_kinetic_energy_density( Returns ------- - posdef_kindetic_energy_density : np.ndarray(N,) + posdef_kinetic_energy_density : np.ndarray(N,) Positive definite kinetic energy density of the given transformed basis set evaluated at `N` grid points. """ - posdef_kindetic_energy_density = np.zeros(points.shape[0]) + posdef_kinetic_energy_density = np.zeros(points.shape[0]) for orders in np.identity(3, dtype=int): - posdef_kindetic_energy_density += evaluate_deriv_reduced_density_matrix( + posdef_kinetic_energy_density += evaluate_deriv_reduced_density_matrix( orders, orders, one_density_matrix, @@ -779,10 +811,10 @@ def evaluate_posdef_kinetic_energy_density( deriv_type=deriv_type, ) # Fix #117: check magnitude of small negative density values, then use clip to remove them - min_output = np.min(output) + min_output = np.min(posdef_kinetic_energy_density) if min_output < 0.0 and abs(min_output) > threshold: raise ValueError(f"Found negative density <= {-threshold}, got {min_output}.") - return (0.5 * posdef_kindetic_energy_density).clip(min=0.0) + return (0.5 * posdef_kinetic_energy_density).clip(min=0.0) # TODO: test against a reference diff --git a/tests/test_density.py b/tests/test_density.py index 4fe47ca3..ef205a31 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -78,30 +78,51 @@ def test_evaluate_density(): def test_evaluate_dm_using_evaluated_orbs(): """Test gbasis.evals.density.evaluate_density_using_evaluated_orbs.""" + density_mat = np.array([[1.0, 2.0], [2.0, 3.0]]) orb_eval = np.array([[1.0], [2.0]]) - dens = evaluate_dm_using_evaluated_orbs(density_mat, [orb_eval,orb_eval]) - assert np.all(dens >= 0.0) - assert np.allclose(np.einsum('iikl->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) + dens = evaluate_dm_using_evaluated_orbs(density_mat, orb_eval,orb_eval) + assert np.all(dens >= 0.0) + assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) + density_mat = np.array([[1.0, 2.0], [2.0, 3.0]]) orb_eval = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) - dens = evaluate_dm_using_evaluated_orbs(density_mat, [orb_eval,orb_eval]) + dens = evaluate_dm_using_evaluated_orbs(density_mat, orb_eval,orb_eval) + + assert np.all(dens >= 0) - assert np.allclose(np.einsum('iikl->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) + assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) def test_evaluate_dm_density(): """Test gbasis.evals.density.evaluate_density.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") transform = np.random.rand(14, 18) + density = np.random.rand(14, 14) density += density.T points = np.random.rand(100, 3) + evaluate_orbs = evaluate_basis(basis, points, transform) - dens = evaluate_dm_density(density, basis, [points], transform) - assert np.allclose(np.einsum('iikl->i',dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) + + dens = evaluate_dm_density(density, basis, points, transform=transform) + assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) + + basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) + basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") + transform = np.random.rand(14, 18) + + density = np.random.rand(14, 14) + density += density.T + points = np.random.rand(100, 3) + + + evaluate_orbs = evaluate_basis(basis, points, transform) + + dens = evaluate_dm_density(density, basis, points, point2=points, transform=transform) + assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) def test_evaluate_hole_x2(): @@ -110,12 +131,27 @@ def test_evaluate_hole_x2(): transform = np.ones((18, 18)) density = np.load(find_datafile("data_rdm_kr_sto6g.npy")) - points = np.random.rand(10, 3) + points1 = np.random.rand(10, 3) + points2 = np.random.rand(12, 3) - eh = evaluate_hole_x2(density, basis, points_list=[points,points],transform=transform) - ed = evaluate_density(density,basis,points,transform=transform) + eh = evaluate_hole_x2(density, basis, points1,point2=points2,transform=transform) + ed = evaluate_density(density,basis,points1,transform=transform) - assert np.allclose(np.einsum('ijkl,j->',eh,ed)/np.sum(ed)/len(points)*len(transform), -1) + assert np.allclose(eh, np.ones((len(points1),len(points2)))*-1) + + basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) + basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") + transform = np.ones((18, 18)) + density = np.load(find_datafile("data_rdm_kr_sto6g.npy")) + + points1 = np.random.rand(10, 3) + + + eh = evaluate_hole_x2(density, basis, points1,transform=transform) + ed = evaluate_density(density,basis,points1,transform=transform) + + + assert np.allclose(eh, np.ones((len(points1),len(points1)))*-1) def test_evaluate_deriv_density(): @@ -492,5 +528,3 @@ def test_evaluate_general_kinetic_energy_density(): evaluate_posdef_kinetic_energy_density(np.identity(40), basis, points, np.identity(40)) + evaluate_density_laplacian(np.identity(40), basis, points, np.identity(40)), ) - -test_evaluate_hole_x2() \ No newline at end of file From a0434ef01a116d56e206208f7b433a58281b0478 Mon Sep 17 00:00:00 2001 From: Wesley Sanchez Date: Thu, 1 Feb 2024 11:43:29 -0500 Subject: [PATCH 7/9] fix conflict --- gbasis/evals/density.py | 120 ++++++++++++++++++++++++---------------- tests/test_density.py | 70 ++++++++++++++++------- 2 files changed, 123 insertions(+), 67 deletions(-) diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py index a086f08f..bd837a0d 100644 --- a/gbasis/evals/density.py +++ b/gbasis/evals/density.py @@ -102,22 +102,25 @@ def evaluate_density(one_density_matrix, basis, points, transform=None, threshol raise ValueError(f"Found negative density <= {-threshold}, got {min_density}.") return density.clip(min=0.0) - -def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): +def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_1,orb_eval_2): """Return the evaluation of the density matrix given the evaluated orbitals. Parameters ---------- one_density_matrix : np.ndarray(K_orb, K_orb) One-electron density matrix in terms of the given orbitals. - orb_eval_list : list of orbitals evaluated on grid points [np.ndarray(K_orb, N) ...] + orb_eval_1 : orbital evaluated on grid points np.ndarray(K_orb, N) + Orbitals evaluated at :math:`N` grid points. + The set of orbitals must be the same basis set used to build the one-electron density + matrix. + orb_eval_2 : orbital evaluated on grid points np.ndarray(K_orb, N) Orbitals evaluated at :math:`N` grid points. The set of orbitals must be the same basis set used to build the one-electron density matrix. Returns ------- - density : np.ndarray(N1,N2,K_orb,K_orb) + dm_on_grid : np.ndarray(N1,N2) Density Matrix evaluated at `N1,N2` grid points. Raises @@ -136,26 +139,33 @@ def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_list): "One-electron density matrix must be a two-dimensional `numpy` array with `dtype`" " float." ) - for orb in orb_eval_list: - if not (isinstance(orb, np.ndarray) and orb.ndim == 2 and orb.dtype == float): - raise TypeError( - "Evaluation of orbitals must be a two-dimensional `numpy` array with `dtype` float." - ) + + if not (isinstance(orb_eval_1, np.ndarray) and orb_eval_1.ndim == 2 and orb_eval_1.dtype == float): + raise TypeError( + "Evaluation of orbitals must be a two-dimensional `numpy` array with `dtype` float." + ) + if not (isinstance(orb_eval_2, np.ndarray) and orb_eval_2.ndim == 2 and orb_eval_2.dtype == float): + raise TypeError( + "Evaluation of orbitals must be a two-dimensional `numpy` array with `dtype` float." + ) if one_density_matrix.shape[0] != one_density_matrix.shape[1]: raise ValueError("One-electron density matrix must be a square matrix.") if not np.allclose(one_density_matrix, one_density_matrix.T): raise ValueError("One-electron density matrix must be symmetric.") - # Tensor product for \gamma(\mathbf{r}_1,\mathbf{r}_2) = \sum_{pq} \gamma_{pq} \chi_p(\mathbf{r}_1) \chi_q(\mathbf{r}_2) - density = np.einsum("ij,ik,jl->klij", one_density_matrix, orb_eval_list[0], orb_eval_list[1]) + #Evaluation of \gamma(\mathbf{r}_1,\mathbf{r}_2) = \sum_{pq} \gamma_{pq} \chi_p(\mathbf{r}_1) \chi_q(\mathbf{r}_2) + dm_on_grid=0 + for ii, orb1 in enumerate(one_density_matrix): + for jj, orb2 in enumerate(one_density_matrix): + dm_on_grid += one_density_matrix[ii][jj]* np.outer(orb_eval_1[ii],orb_eval_2[jj]) - # returns dm evaluated on each grid point - return density + #returns dm evaluated on each grid point summed over the orbitals + return dm_on_grid -def evaluate_dm_density(one_density_matrix, basis, points_list, transform=None): - r"""Return the density of the given basis set at the given points. +def evaluate_dm_density(one_density_matrix, basis, point1, point2 = np.array(None), transform=None): + r"""Return the density matrix of the given basis set evaluated on the given points. Parameters ---------- @@ -165,7 +175,14 @@ def evaluate_dm_density(one_density_matrix, basis, points_list, transform=None): be expressed with respect to the transformed basis set. basis : list/tuple of GeneralizedContractionShell Shells of generalized contractions. - points_list : list of points [np.ndarray(N, 3)...] + point1 : set of points np.ndarray(N, 3) + Cartesian coordinates of the points in space (in atomic units) where the basis functions + are evaluated. + Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` + components. + This function can take a list of points at which basis functions are evaluated. If only one + set of points is given, it will be duplicated. + point2 : optional set of points np.ndarray(N, 3) --> default None Cartesian coordinates of the points in space (in atomic units) where the basis functions are evaluated. Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` @@ -182,27 +199,28 @@ def evaluate_dm_density(one_density_matrix, basis, points_list, transform=None): Returns ------- - dm_on_grid : np.ndarray(N1,N2,K_orb,K_orb) - Density Matrix evaluated at `N1,N2` grid points. + dm_on_grid : np.ndarray(N1,N2) + Density Matrix evaluated at `N1,N2` grid points summed over orbitals. """ - orb_evals = [] - # evaluate basi(e)s on the point set(s) - for grid in points_list: - orb_eval = evaluate_basis(basis, grid, transform=transform) - orb_evals.append(orb_eval) - # if only one set of points is supplied, it is duplicated - if len(points_list) == 1: - orb_evals.append(orb_eval) + #evaluate basi(e)s on the point set(s) + orb_eval_1 = evaluate_basis(basis, point1, transform=transform) - # Calulated performed using the evaluated orbitals - dm_on_grid = evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_evals) + #if only one set of points is supplied, it is duplicated + if point2.all() == None: + orb_eval_2 = orb_eval_1 + else: + orb_eval_2 = evaluate_basis(basis, point2, transform=transform) + + + #Calulated performed using the evaluated orbitals + dm_on_grid = evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_1,orb_eval_2) return dm_on_grid -def evaluate_hole_x2(one_density_matrix, basis, points_list, transform=None): +def evaluate_hole_x2(one_density_matrix, basis, point1, point2=np.array(None), transform=None): r"""Return the two-body hole correlation function. .. math :: @@ -221,7 +239,14 @@ def evaluate_hole_x2(one_density_matrix, basis, points_list, transform=None): be expressed with respect to the transformed basis set. basis : list/tuple of GeneralizedContractionShell Shells of generalized contractions. - points_list : list of points [np.ndarray(N, 3)] + point1 : set of points np.ndarray(N, 3) + Cartesian coordinates of the points in space (in atomic units) where the basis functions + are evaluated. + Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` + components. + This function can take a list of points at which basis functions are evaluated. If only one + set of points is given, it will be duplicated. + point2 : optional set of points np.ndarray(N, 3) --> default None Cartesian coordinates of the points in space (in atomic units) where the basis functions are evaluated. Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` @@ -238,25 +263,26 @@ def evaluate_hole_x2(one_density_matrix, basis, points_list, transform=None): Returns ------- - hole_x2 : np.ndarray(N1,N2,K_orb,K_orb) + hole_x2 : np.ndarray(N1,N2) Two-body Exchange Hole evaluated at `N1,N2` grid points. """ - dens_evals = [] - for grid in points_list: - dens_eval = evaluate_density(one_density_matrix, basis, grid, transform=transform) - dens_evals.append(dens_eval) + dens_eval_1 = evaluate_density(one_density_matrix,basis, point1, transform=transform) + + #if only one set of points is supplied, density evaluation is duplicated + if point2.all() == None: + dens_eval_2 = dens_eval_1 + else: + dens_eval_2 = evaluate_density(one_density_matrix,basis, point2, transform=transform) - if len(points_list) == 1: - dens_evals.append(dens_eval) + #build density matrix on grid + dm_eval = evaluate_dm_density(one_density_matrix, basis, point1, point2, transform=transform) - # build density matrix on grid - dm_eval = evaluate_dm_density(one_density_matrix, basis, points_list, transform=transform) + #evaluate hole function + numerator = dm_eval*dm_eval - # evaluate hole function - numerator = np.einsum("ijkl,ijkl->ijkl", dm_eval, dm_eval) - hole_x2 = -1 * np.einsum("ijkl,i,j->ijkl", numerator, 1 / dens_evals[0], 1 / dens_evals[1]) + hole_x2 = -1*np.einsum('ij,i,j->ij',numerator,1/dens_eval_1,1/dens_eval_2) return hole_x2 @@ -763,14 +789,14 @@ def evaluate_posdef_kinetic_energy_density( Returns ------- - posdef_kindetic_energy_density : np.ndarray(N,) + posdef_kinetic_energy_density : np.ndarray(N,) Positive definite kinetic energy density of the given transformed basis set evaluated at `N` grid points. """ - posdef_kindetic_energy_density = np.zeros(points.shape[0]) + posdef_kinetic_energy_density = np.zeros(points.shape[0]) for orders in np.identity(3, dtype=int): - posdef_kindetic_energy_density += evaluate_deriv_reduced_density_matrix( + posdef_kinetic_energy_density += evaluate_deriv_reduced_density_matrix( orders, orders, one_density_matrix, @@ -780,10 +806,10 @@ def evaluate_posdef_kinetic_energy_density( deriv_type=deriv_type, ) # Fix #117: check magnitude of small negative density values, then use clip to remove them - min_output = np.min(output) + min_output = np.min(posdef_kinetic_energy_density) if min_output < 0.0 and abs(min_output) > threshold: raise ValueError(f"Found negative density <= {-threshold}, got {min_output}.") - return (0.5 * posdef_kindetic_energy_density).clip(min=0.0) + return (0.5 * posdef_kinetic_energy_density).clip(min=0.0) # TODO: test against a reference diff --git a/tests/test_density.py b/tests/test_density.py index 03e6f5d6..a7c1d326 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -79,37 +79,51 @@ def test_evaluate_density(): def test_evaluate_dm_using_evaluated_orbs(): """Test gbasis.evals.density.evaluate_density_using_evaluated_orbs.""" + density_mat = np.array([[1.0, 2.0], [2.0, 3.0]]) orb_eval = np.array([[1.0], [2.0]]) - dens = evaluate_dm_using_evaluated_orbs(density_mat, [orb_eval, orb_eval]) - assert np.all(dens >= 0.0) - assert np.allclose( - np.einsum("iikl->i", dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval) - ) + dens = evaluate_dm_using_evaluated_orbs(density_mat, orb_eval,orb_eval) + assert np.all(dens >= 0.0) + assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) + density_mat = np.array([[1.0, 2.0], [2.0, 3.0]]) orb_eval = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) - dens = evaluate_dm_using_evaluated_orbs(density_mat, [orb_eval, orb_eval]) - assert np.all(dens >= 0) - assert np.allclose( - np.einsum("iikl->i", dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval) - ) + dens = evaluate_dm_using_evaluated_orbs(density_mat, orb_eval,orb_eval) + assert np.all(dens >= 0) + assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) + def test_evaluate_dm_density(): """Test gbasis.evals.density.evaluate_density.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") transform = np.random.rand(14, 18) + density = np.random.rand(14, 14) density += density.T points = np.random.rand(100, 3) + evaluate_orbs = evaluate_basis(basis, points, transform) - dens = evaluate_dm_density(density, basis, [points], transform) - assert np.allclose( - np.einsum("iikl->i", dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs) - ) + + dens = evaluate_dm_density(density, basis, points, transform=transform) + assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) + + basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) + basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") + transform = np.random.rand(14, 18) + + density = np.random.rand(14, 14) + density += density.T + points = np.random.rand(100, 3) + + + evaluate_orbs = evaluate_basis(basis, points, transform) + + dens = evaluate_dm_density(density, basis, points, point2=points, transform=transform) + assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) def test_evaluate_hole_x2(): @@ -118,14 +132,27 @@ def test_evaluate_hole_x2(): transform = np.ones((18, 18)) density = np.load(find_datafile("data_rdm_kr_sto6g.npy")) - points = np.random.rand(10, 3) + points1 = np.random.rand(10, 3) + points2 = np.random.rand(12, 3) - eh = evaluate_hole_x2(density, basis, points_list=[points, points], transform=transform) - ed = evaluate_density(density, basis, points, transform=transform) + eh = evaluate_hole_x2(density, basis, points1,point2=points2,transform=transform) + ed = evaluate_density(density,basis,points1,transform=transform) - assert np.allclose( - np.einsum("ijkl,j->", eh, ed) / np.sum(ed) / len(points) * len(transform), -1 - ) + assert np.allclose(eh, np.ones((len(points1),len(points2)))*-1) + + basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) + basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") + transform = np.ones((18, 18)) + density = np.load(find_datafile("data_rdm_kr_sto6g.npy")) + + points1 = np.random.rand(10, 3) + + + eh = evaluate_hole_x2(density, basis, points1,transform=transform) + ed = evaluate_density(density,basis,points1,transform=transform) + + + assert np.allclose(eh, np.ones((len(points1),len(points1)))*-1) def test_evaluate_deriv_density(): @@ -502,6 +529,9 @@ def test_evaluate_general_kinetic_energy_density(): evaluate_posdef_kinetic_energy_density(np.identity(40), basis, points, np.identity(40)) + evaluate_density_laplacian(np.identity(40), basis, points, np.identity(40)), ) +<<<<<<< HEAD test_evaluate_hole_x2() +======= +>>>>>>> 59c34a5 (code changes+update functions) From 14c973cc654a6371c9f3ad3ea9aaa952d037b418 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 1 Feb 2024 16:48:37 +0000 Subject: [PATCH 8/9] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- gbasis/evals/density.py | 60 ++++++++++++++++++++--------------------- tests/test_density.py | 44 +++++++++++++++--------------- 2 files changed, 53 insertions(+), 51 deletions(-) diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py index 1710b22f..1f658e25 100644 --- a/gbasis/evals/density.py +++ b/gbasis/evals/density.py @@ -102,7 +102,8 @@ def evaluate_density(one_density_matrix, basis, points, transform=None, threshol raise ValueError(f"Found negative density <= {-threshold}, got {min_density}.") return density.clip(min=0.0) -def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_1,orb_eval_2): + +def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_1, orb_eval_2): """Return the evaluation of the density matrix given the evaluated orbitals. Parameters @@ -140,11 +141,15 @@ def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_1,orb_eval_2): " float." ) - if not (isinstance(orb_eval_1, np.ndarray) and orb_eval_1.ndim == 2 and orb_eval_1.dtype == float): + if not ( + isinstance(orb_eval_1, np.ndarray) and orb_eval_1.ndim == 2 and orb_eval_1.dtype == float + ): raise TypeError( "Evaluation of orbitals must be a two-dimensional `numpy` array with `dtype` float." ) - if not (isinstance(orb_eval_2, np.ndarray) and orb_eval_2.ndim == 2 and orb_eval_2.dtype == float): + if not ( + isinstance(orb_eval_2, np.ndarray) and orb_eval_2.ndim == 2 and orb_eval_2.dtype == float + ): raise TypeError( "Evaluation of orbitals must be a two-dimensional `numpy` array with `dtype` float." ) @@ -154,17 +159,17 @@ def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_1,orb_eval_2): if not np.allclose(one_density_matrix, one_density_matrix.T): raise ValueError("One-electron density matrix must be symmetric.") - #Evaluation of \gamma(\mathbf{r}_1,\mathbf{r}_2) = \sum_{pq} \gamma_{pq} \chi_p(\mathbf{r}_1) \chi_q(\mathbf{r}_2) - dm_on_grid=0 + # Evaluation of \gamma(\mathbf{r}_1,\mathbf{r}_2) = \sum_{pq} \gamma_{pq} \chi_p(\mathbf{r}_1) \chi_q(\mathbf{r}_2) + dm_on_grid = 0 for ii, orb1 in enumerate(one_density_matrix): for jj, orb2 in enumerate(one_density_matrix): - dm_on_grid += one_density_matrix[ii][jj]* np.outer(orb_eval_1[ii],orb_eval_2[jj]) + dm_on_grid += one_density_matrix[ii][jj] * np.outer(orb_eval_1[ii], orb_eval_2[jj]) - #returns dm evaluated on each grid point summed over the orbitals + # returns dm evaluated on each grid point summed over the orbitals return dm_on_grid -def evaluate_dm_density(one_density_matrix, basis, point1, point2 = np.array(None), transform=None): +def evaluate_dm_density(one_density_matrix, basis, point1, point2=np.array(None), transform=None): r"""Return the density matrix of the given basis set evaluated on the given points. Parameters @@ -180,14 +185,14 @@ def evaluate_dm_density(one_density_matrix, basis, point1, point2 = np.array(Non are evaluated. Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` components. - This function can take a list of points at which basis functions are evaluated. If only one + This function can take a list of points at which basis functions are evaluated. If only one set of points is given, it will be duplicated. point2 : optional set of points np.ndarray(N, 3) --> default None Cartesian coordinates of the points in space (in atomic units) where the basis functions are evaluated. Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` components. - This function can take a list of points at which basis functions are evaluated. If only one + This function can take a list of points at which basis functions are evaluated. If only one set of points is given, it will be duplicated. point2 : optional set of points np.ndarray(N, 3) --> default None Cartesian coordinates of the points in space (in atomic units) where the basis functions @@ -211,23 +216,19 @@ def evaluate_dm_density(one_density_matrix, basis, point1, point2 = np.array(Non """ - #evaluate basi(e)s on the point set(s) + # evaluate basi(e)s on the point set(s) orb_eval_1 = evaluate_basis(basis, point1, transform=transform) - #if only one set of points is supplied, it is duplicated + # if only one set of points is supplied, it is duplicated if point2.all() == None: orb_eval_2 = orb_eval_1 else: - orb_eval_2 = evaluate_basis(basis, point2, transform=transform) - - - #Calulated performed using the evaluated orbitals - dm_on_grid = evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_1,orb_eval_2) - - return dm_on_grid - + orb_eval_2 = evaluate_basis(basis, point2, transform=transform) + # Calulated performed using the evaluated orbitals + dm_on_grid = evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_1, orb_eval_2) + return dm_on_grid def evaluate_hole_x2(one_density_matrix, basis, point1, point2=np.array(None), transform=None): @@ -254,14 +255,14 @@ def evaluate_hole_x2(one_density_matrix, basis, point1, point2=np.array(None), t are evaluated. Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` components. - This function can take a list of points at which basis functions are evaluated. If only one + This function can take a list of points at which basis functions are evaluated. If only one set of points is given, it will be duplicated. point2 : optional set of points np.ndarray(N, 3) --> default None Cartesian coordinates of the points in space (in atomic units) where the basis functions are evaluated. Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` components. - This function can take a list of points at which basis functions are evaluated. If only one + This function can take a list of points at which basis functions are evaluated. If only one set of points is given, it will be duplicated. point2 : optional set of points np.ndarray(N, 3) --> default None Cartesian coordinates of the points in space (in atomic units) where the basis functions @@ -285,22 +286,21 @@ def evaluate_hole_x2(one_density_matrix, basis, point1, point2=np.array(None), t """ - dens_eval_1 = evaluate_density(one_density_matrix,basis, point1, transform=transform) + dens_eval_1 = evaluate_density(one_density_matrix, basis, point1, transform=transform) - #if only one set of points is supplied, density evaluation is duplicated + # if only one set of points is supplied, density evaluation is duplicated if point2.all() == None: dens_eval_2 = dens_eval_1 else: - dens_eval_2 = evaluate_density(one_density_matrix,basis, point2, transform=transform) + dens_eval_2 = evaluate_density(one_density_matrix, basis, point2, transform=transform) - #build density matrix on grid + # build density matrix on grid dm_eval = evaluate_dm_density(one_density_matrix, basis, point1, point2, transform=transform) - #evaluate hole function - numerator = dm_eval*dm_eval - - hole_x2 = -1*np.einsum('ij,i,j->ij',numerator,1/dens_eval_1,1/dens_eval_2) + # evaluate hole function + numerator = dm_eval * dm_eval + hole_x2 = -1 * np.einsum("ij,i,j->ij", numerator, 1 / dens_eval_1, 1 / dens_eval_2) return hole_x2 diff --git a/tests/test_density.py b/tests/test_density.py index 6dd84c1b..6c0b5f0c 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -79,21 +79,25 @@ def test_evaluate_density(): def test_evaluate_dm_using_evaluated_orbs(): """Test gbasis.evals.density.evaluate_density_using_evaluated_orbs.""" - + density_mat = np.array([[1.0, 2.0], [2.0, 3.0]]) orb_eval = np.array([[1.0], [2.0]]) - dens = evaluate_dm_using_evaluated_orbs(density_mat, orb_eval,orb_eval) + dens = evaluate_dm_using_evaluated_orbs(density_mat, orb_eval, orb_eval) assert np.all(dens >= 0.0) - assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) - + assert np.allclose( + np.einsum("ii->i", dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval) + ) + density_mat = np.array([[1.0, 2.0], [2.0, 3.0]]) orb_eval = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) - dens = evaluate_dm_using_evaluated_orbs(density_mat, orb_eval,orb_eval) - + dens = evaluate_dm_using_evaluated_orbs(density_mat, orb_eval, orb_eval) assert np.all(dens >= 0) - assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval)) + assert np.allclose( + np.einsum("ii->i", dens), np.einsum("ij,ik,jk->k", density_mat, orb_eval, orb_eval) + ) + def test_evaluate_dm_density(): """Test gbasis.evals.density.evaluate_density.""" @@ -105,11 +109,12 @@ def test_evaluate_dm_density(): density += density.T points = np.random.rand(100, 3) - evaluate_orbs = evaluate_basis(basis, points, transform) dens = evaluate_dm_density(density, basis, points, transform=transform) - assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) + assert np.allclose( + np.einsum("ii->i", dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs) + ) basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") @@ -119,11 +124,12 @@ def test_evaluate_dm_density(): density += density.T points = np.random.rand(100, 3) - evaluate_orbs = evaluate_basis(basis, points, transform) dens = evaluate_dm_density(density, basis, points, point2=points, transform=transform) - assert np.allclose(np.einsum('ii->i',dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs)) + assert np.allclose( + np.einsum("ii->i", dens), np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs) + ) def test_evaluate_hole_x2(): @@ -135,10 +141,10 @@ def test_evaluate_hole_x2(): points1 = np.random.rand(10, 3) points2 = np.random.rand(12, 3) - eh = evaluate_hole_x2(density, basis, points1,point2=points2,transform=transform) - ed = evaluate_density(density,basis,points1,transform=transform) + eh = evaluate_hole_x2(density, basis, points1, point2=points2, transform=transform) + ed = evaluate_density(density, basis, points1, transform=transform) - assert np.allclose(eh, np.ones((len(points1),len(points2)))*-1) + assert np.allclose(eh, np.ones((len(points1), len(points2))) * -1) basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") @@ -147,12 +153,10 @@ def test_evaluate_hole_x2(): points1 = np.random.rand(10, 3) + eh = evaluate_hole_x2(density, basis, points1, transform=transform) + ed = evaluate_density(density, basis, points1, transform=transform) - eh = evaluate_hole_x2(density, basis, points1,transform=transform) - ed = evaluate_density(density,basis,points1,transform=transform) - - - assert np.allclose(eh, np.ones((len(points1),len(points1)))*-1) + assert np.allclose(eh, np.ones((len(points1), len(points1))) * -1) def test_evaluate_deriv_density(): @@ -529,5 +533,3 @@ def test_evaluate_general_kinetic_energy_density(): evaluate_posdef_kinetic_energy_density(np.identity(40), basis, points, np.identity(40)) + evaluate_density_laplacian(np.identity(40), basis, points, np.identity(40)), ) - - From fe9fef95cdc9ad358704c0a5ca2a6d04f6bcd098 Mon Sep 17 00:00:00 2001 From: leila-pujal Date: Wed, 3 Jul 2024 16:47:32 -0400 Subject: [PATCH 9/9] Change np.all(None) to None for numpy 2.0 --- gbasis/evals/density.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py index 1f658e25..8c6e2069 100644 --- a/gbasis/evals/density.py +++ b/gbasis/evals/density.py @@ -169,7 +169,7 @@ def evaluate_dm_using_evaluated_orbs(one_density_matrix, orb_eval_1, orb_eval_2) return dm_on_grid -def evaluate_dm_density(one_density_matrix, basis, point1, point2=np.array(None), transform=None): +def evaluate_dm_density(one_density_matrix, basis, point1, point2=None, transform=None): r"""Return the density matrix of the given basis set evaluated on the given points. Parameters @@ -216,11 +216,11 @@ def evaluate_dm_density(one_density_matrix, basis, point1, point2=np.array(None) """ - # evaluate basi(e)s on the point set(s) + # Evaluate basis on the point set orb_eval_1 = evaluate_basis(basis, point1, transform=transform) - # if only one set of points is supplied, it is duplicated - if point2.all() == None: + # If only one set of points is supplied, it is duplicated + if point2 is None: orb_eval_2 = orb_eval_1 else: orb_eval_2 = evaluate_basis(basis, point2, transform=transform) @@ -231,8 +231,8 @@ def evaluate_dm_density(one_density_matrix, basis, point1, point2=np.array(None) return dm_on_grid -def evaluate_hole_x2(one_density_matrix, basis, point1, point2=np.array(None), transform=None): - r"""Return the two-body hole correlation function. +def evaluate_hole_x2(one_density_matrix, basis, point1, point2=None, transform=None): + r"""Return the two-body exchange hole correlation function. .. math :: @@ -289,7 +289,7 @@ def evaluate_hole_x2(one_density_matrix, basis, point1, point2=np.array(None), t dens_eval_1 = evaluate_density(one_density_matrix, basis, point1, transform=transform) # if only one set of points is supplied, density evaluation is duplicated - if point2.all() == None: + if point2 is None: dens_eval_2 = dens_eval_1 else: dens_eval_2 = evaluate_density(one_density_matrix, basis, point2, transform=transform)