diff --git a/example/kbe_polyacetylene.py b/example/kbe_polyacetylene.py index 1ce7dbb2..22f0f613 100644 --- a/example/kbe_polyacetylene.py +++ b/example/kbe_polyacetylene.py @@ -2,7 +2,7 @@ # A supercell with 4 carbon & 4 hydrogen atoms is defined as unit cell in # pyscf's periodic HF calculation -import numpy +import numpy as np from pyscf.pbc import df, gto, scf from quemb.kbe import BE, fragpart @@ -14,7 +14,7 @@ b = 8.0 c = 2.455 * 2.0 -lat = numpy.eye(3) +lat = np.eye(3) lat[0, 0] = a lat[1, 1] = b lat[2, 2] = c diff --git a/example/molbe_dmrg_block2.py b/example/molbe_dmrg_block2.py index 2f06c52e..d1e99e32 100644 --- a/example/molbe_dmrg_block2.py +++ b/example/molbe_dmrg_block2.py @@ -3,7 +3,7 @@ # Garnet-Chan group at Caltech: https://block2.readthedocs.io/en/latest/index.html import matplotlib.pyplot as plt -import numpy +import numpy as np from pyscf import cc, fci, gto, scf from quemb.molbe import BE, fragpart @@ -11,7 +11,7 @@ # We'll consider the dissociation curve for a 1D chain of 8 H-atoms: num_points = 3 -seps = numpy.linspace(0.60, 1.6, num=num_points) +seps = np.linspace(0.60, 1.6, num=num_points) fci_ecorr, ccsd_ecorr, ccsdt_ecorr, bedmrg_ecorr = [], [], [], [] for a in seps: diff --git a/src/quemb/kbe/autofrag.py b/src/quemb/kbe/autofrag.py index 384d5628..fc62dcb3 100644 --- a/src/quemb/kbe/autofrag.py +++ b/src/quemb/kbe/autofrag.py @@ -1,7 +1,7 @@ # Author(s): Oinam Romesh Meitei -import numpy +from numpy import arange, asarray, where from numpy.linalg import norm from pyscf import lib @@ -83,10 +83,10 @@ def sidefunc( rlist=[], ): if ext_list == []: - main_list.extend(unit2[numpy.where(unit1 == Idx)[0]]) - sub_list.extend(unit2[numpy.where(unit1 == Idx)[0]]) + main_list.extend(unit2[where(unit1 == Idx)[0]]) + sub_list.extend(unit2[where(unit1 == Idx)[0]]) else: - for sub_i in unit2[numpy.where(unit1 == Idx)[0]]: + for sub_i in unit2[where(unit1 == Idx)[0]]: if sub_i in rlist: continue if sub_i in ext_list: @@ -104,7 +104,7 @@ def sidefunc( close_be3 = [] if be_type == "be3" or be_type == "be4": - for lmin1 in unit2[numpy.where(unit1 == Idx)[0]]: + for lmin1 in unit2[where(unit1 == Idx)[0]]: for jdx, j in enumerate(coord): if ( jdx not in unit1 @@ -396,11 +396,9 @@ def autogen( lnk2 = 1 lattice_vector = cell.lattice_vectors() - Ts = lib.cartesian_prod( - (numpy.arange(lkpt[0]), numpy.arange(lkpt[1]), numpy.arange(lkpt[2])) - ) + Ts = lib.cartesian_prod((arange(lkpt[0]), arange(lkpt[1]), arange(lkpt[2]))) - Ls = numpy.dot(Ts, lattice_vector) + Ls = Ts @ lattice_vector # 1-2-(1-2)-1-2 # * * @@ -469,12 +467,12 @@ def autogen( kmsites = [] ktsites = [] - lunit = numpy.asarray(lunit) - runit = numpy.asarray(runit) - uunit = numpy.asarray(uunit) - dunit = numpy.asarray(dunit) - munit = numpy.asarray(munit) - tunit = numpy.asarray(tunit) + lunit = asarray(lunit) + runit = asarray(runit) + uunit = asarray(uunit) + dunit = asarray(dunit) + munit = asarray(munit) + tunit = asarray(tunit) inter_dist = 1000.0 if twoD and interlayer: @@ -1307,7 +1305,7 @@ def autogen( continue if be_type == "be3" or be_type == "be4": if jdx in lunit: - lmin1 = runit[numpy.where(lunit == jdx)[0]] + lmin1 = runit[where(lunit == jdx)[0]] if not twoD: flist.extend(lmin1) lsts.extend(lmin1) @@ -1387,7 +1385,7 @@ def autogen( bond=bond, ) if jdx in runit: - rmin1 = lunit[numpy.where(runit == jdx)[0]] + rmin1 = lunit[where(runit == jdx)[0]] if not twoD: flist.extend(rmin1) rsts.extend(rmin1) @@ -1467,7 +1465,7 @@ def autogen( ) if jdx in uunit: - umin1 = dunit[numpy.where(uunit == jdx)[0]] + umin1 = dunit[where(uunit == jdx)[0]] add_check_k(umin1, flist, usts, kusts, 2) if be_type == "be4": for kdx, k in enumerate(coord): @@ -1534,7 +1532,7 @@ def autogen( bond=bond, ) if jdx in dunit: - dmin1 = uunit[numpy.where(dunit == jdx)[0]] + dmin1 = uunit[where(dunit == jdx)[0]] add_check_k(dmin1, flist, dsts, kdsts, nk1) if be_type == "be4": @@ -1602,7 +1600,7 @@ def autogen( bond=bond, ) if jdx in munit: # - mmin1 = tunit[numpy.where(munit == jdx)[0]] + mmin1 = tunit[where(munit == jdx)[0]] add_check_k(mmin1, flist, msts, kmsts, nk1 * nk2) if be_type == "be4": for kdx, k in enumerate(coord): @@ -1656,7 +1654,7 @@ def autogen( ) if jdx in tunit: - tmin1 = munit[numpy.where(tunit == jdx)[0]] + tmin1 = munit[where(tunit == jdx)[0]] add_check_k(tmin1, flist, tsts, ktsts, nk1 + 2) if be_type == "be4": @@ -1724,7 +1722,7 @@ def autogen( pedg.append(kdx) if be_type == "be4": if kdx in lunit: - lmin1 = runit[numpy.where(lunit == kdx)[0]] + lmin1 = runit[where(lunit == kdx)[0]] for zdx in lmin1: if ( zdx in lsts @@ -1740,7 +1738,7 @@ def autogen( else: klsts.append(nk1) if kdx in runit: - rmin1 = lunit[numpy.where(runit == kdx)[0]] + rmin1 = lunit[where(runit == kdx)[0]] for zdx in rmin1: if ( zdx in rsts @@ -1755,7 +1753,7 @@ def autogen( else: krsts.append(2) if kdx in uunit: - umin1 = dunit[numpy.where(uunit == kdx)[0]] + umin1 = dunit[where(uunit == kdx)[0]] for zdx in umin1: if ( zdx in usts @@ -1767,7 +1765,7 @@ def autogen( usts.append(zdx) kusts.append(2) if kdx in dunit: - dmin1 = uunit[numpy.where(dunit == kdx)[0]] + dmin1 = uunit[where(dunit == kdx)[0]] for zdx in dmin1: if ( zdx in dsts @@ -1779,7 +1777,7 @@ def autogen( dsts.append(zdx) kdsts.append(nk1) if kdx in munit: - mmin1 = tunit[numpy.where(munit == kdx)[0]] + mmin1 = tunit[where(munit == kdx)[0]] for zdx in mmin1: if ( zdx in msts @@ -1791,7 +1789,7 @@ def autogen( msts.append(zdx) kmsts.append(nk1 * nk2) if kdx in tunit: - tmin1 = munit[numpy.where(tunit == kdx)[0]] + tmin1 = munit[where(tunit == kdx)[0]] for zdx in tmin1: if ( zdx in tsts diff --git a/src/quemb/kbe/helper.py b/src/quemb/kbe/helper.py index 1a8860e5..dd3c52b6 100644 --- a/src/quemb/kbe/helper.py +++ b/src/quemb/kbe/helper.py @@ -1,6 +1,6 @@ # Author(s): Oinam Romesh Meitei -import numpy +from numpy import asarray, complex128, float64, zeros from numpy.linalg import multi_dot from pyscf import scf @@ -33,21 +33,21 @@ def get_veff(eri_, dm, S, TA, hf_veff, return_veff0=False): # construct rdm nk, nao, neo = TA.shape unused(nao) - P_ = numpy.zeros((neo, neo), dtype=numpy.complex128) + P_ = zeros((neo, neo), dtype=complex128) for k in range(nk): - Cinv = numpy.dot(TA[k].conj().T, S[k]) + Cinv = TA[k].conj().T @ S[k] P_ += multi_dot((Cinv, dm[k], Cinv.conj().T)) P_ /= float(nk) - P_ = numpy.asarray(P_.real, dtype=numpy.double) + P_ = asarray(P_.real, dtype=float64) - eri_ = numpy.asarray(eri_, dtype=numpy.double) + eri_ = asarray(eri_, dtype=float64) vj, vk = scf.hf.dot_eri_dm(eri_, P_, hermi=1, with_j=True, with_k=True) Veff_ = vj - 0.5 * vk # remove core contribution from hf_veff - Veff0 = numpy.zeros((neo, neo), dtype=numpy.complex128) + Veff0 = zeros((neo, neo), dtype=complex128) for k in range(nk): Veff0 += multi_dot((TA[k].conj().T, hf_veff[k], TA[k])) Veff0 /= float(nk) diff --git a/src/quemb/kbe/lo.py b/src/quemb/kbe/lo.py index 89b85445..4fca07b4 100644 --- a/src/quemb/kbe/lo.py +++ b/src/quemb/kbe/lo.py @@ -3,8 +3,21 @@ # import os -import numpy +import numpy as np from libdmet.lo import pywannier90 +from numpy import ( + allclose, + array, + asarray, + complex128, + diag, + eye, + hstack, + sqrt, + where, + zeros, + zeros_like, +) from numpy.linalg import eigh, multi_dot, svd from quemb.kbe.lo_k import ( @@ -59,52 +72,52 @@ def localize( if lo_method == "lowdin": # Lowdin orthogonalization with k-points - W = numpy.zeros_like(self.S) + W = zeros_like(self.S) nk, nao, nmo = self.C.shape if self.frozen_core: - W_nocore = numpy.zeros_like(self.S[:, :, self.ncore :]) - lmo_coeff = numpy.zeros_like(self.C[:, self.ncore :, self.ncore :]) - cinv_ = numpy.zeros((nk, nmo - self.ncore, nao), dtype=numpy.complex128) + W_nocore = zeros_like(self.S[:, :, self.ncore :]) + lmo_coeff = zeros_like(self.C[:, self.ncore :, self.ncore :]) + cinv_ = zeros((nk, nmo - self.ncore, nao), dtype=complex128) else: - lmo_coeff = numpy.zeros_like(self.C) - cinv_ = numpy.zeros((nk, nmo, nao), dtype=numpy.complex128) + lmo_coeff = zeros_like(self.C) + cinv_ = zeros((nk, nmo, nao), dtype=complex128) for k in range(self.nkpt): es_, vs_ = eigh(self.S[k]) edx = es_ > 1.0e-14 - W[k] = (vs_[:, edx] / numpy.sqrt(es_[edx])) @ vs_[:, edx].conj().T + W[k] = (vs_[:, edx] / sqrt(es_[edx])) @ vs_[:, edx].conj().T for i in range(W[k].shape[1]): if W[k][i, i] < 0: W[:, i] *= -1 if self.frozen_core: - pcore = numpy.eye(W[k].shape[0]) - (self.P_core[k] @ self.S[k]) - C_ = numpy.dot(pcore, W[k]) + pcore = eye(W[k].shape[0]) - (self.P_core[k] @ self.S[k]) + C_ = pcore @ W[k] # PYSCF has basis in 1s2s3s2p2p2p3p3p3p format # fix no_core_idx - use population for now # C_ = C_[:,self.no_core_idx] Cpop = multi_dot((C_.conj().T, self.S[k], C_)) - Cpop = numpy.diag(Cpop.real) + Cpop = diag(Cpop.real) - no_core_idx = numpy.where(Cpop > 0.7)[0] + no_core_idx = where(Cpop > 0.7)[0] C_ = C_[:, no_core_idx] S_ = multi_dot((C_.conj().T, self.S[k], C_)) es_, vs_ = eigh(S_) edx = es_ > 1.0e-14 - W_ = (vs_[:, edx] / numpy.sqrt(es_[edx])) @ vs_[:, edx].conj().T - W_nocore[k] = numpy.dot(C_, W_) + W_ = (vs_[:, edx] / sqrt(es_[edx])) @ vs_[:, edx].conj().T + W_nocore[k] = C_ @ W_ lmo_coeff[k] = multi_dot( (W_nocore[k].conj().T, self.S[k], self.C[k][:, self.ncore :]), ) - cinv_[k] = numpy.dot(W_nocore[k].conj().T, self.S[k]) + cinv_[k] = W_nocore[k].conj().T @ self.S[k] else: lmo_coeff[k] = multi_dot((W[k].conj().T, self.S[k], self.C[k])) - cinv_[k] = numpy.dot(W[k].conj().T, self.S[k]) + cinv_[k] = W[k].conj().T @ self.S[k] if self.frozen_core: self.W = W_nocore else: @@ -120,7 +133,7 @@ def localize( # tmp - aos are not rearrange and so below is not necessary nk, nao, nlo = ciao_.shape - Ciao_ = numpy.zeros((nk, nao, nlo), dtype=numpy.complex128) + Ciao_ = zeros((nk, nao, nlo), dtype=complex128) for k in range(self.nkpt): aoind_by_atom = get_aoind_by_atom(self.cell) ctmp, iaoind_by_atom = reorder_by_atom_( @@ -136,7 +149,7 @@ def localize( ) nk, nao, nlo = cpao_.shape - Cpao_ = numpy.zeros((nk, nao, nlo), dtype=numpy.complex128) + Cpao_ = zeros((nk, nao, nlo), dtype=complex128) for k in range(self.nkpt): aoind_by_atom = get_aoind_by_atom(self.cell) ctmp, paoind_by_atom = reorder_by_atom_( @@ -147,9 +160,7 @@ def localize( nk, nao, nlo = Ciao_.shape if self.frozen_core: nk, nao, nlo = Ciao_.shape - Ciao_nocore = numpy.zeros( - (nk, nao, nlo - self.ncore), dtype=numpy.complex128 - ) + Ciao_nocore = zeros((nk, nao, nlo - self.ncore), dtype=complex128) for k in range(nk): Ccore = self.C[k][:, : self.ncore] Ciao_nocore[k] = remove_core_mo_k(Ciao_[k], Ccore, self.S[k]) @@ -162,17 +173,17 @@ def localize( s12_core_, s2_core = get_xovlp_k(self.cell, self.kpts, basis=core_basis) C_core_ = self.C[:, :, : self.ncore].copy() nk_, nao_, nmo_ = C_core_.shape - s1_core = numpy.zeros((nk_, nmo_, nmo_), dtype=self.S.dtype) - s12_core = numpy.zeros( + s1_core = zeros((nk_, nmo_, nmo_), dtype=self.S.dtype) + s12_core = zeros( (nk_, nmo_, s12_core_.shape[-1]), dtype=s12_core_.dtype ) - C_core = numpy.zeros((nk_, self.ncore, self.ncore), dtype=C_core_.dtype) + C_core = zeros((nk_, self.ncore, self.ncore), dtype=C_core_.dtype) for k in range(nk_): C_core[k] = C_core_[k].conj().T @ self.S[k] @ C_core_[k] s1_core[k] = C_core_[k].conj().T @ self.S[k] @ C_core_[k] s12_core[k] = C_core_[k].conj().T @ s12_core_[k] ciao_core_ = get_iao_k(C_core, s12_core, s1_core, s2_core, ortho=False) - ciao_core = numpy.zeros( + ciao_core = zeros( (nk_, nao_, ciao_core_.shape[-1]), dtype=ciao_core_.dtype ) for k in range(nk_): @@ -186,11 +197,9 @@ def localize( C_nocore = self.C[:, :, self.ncore :].copy() C_nocore_occ_ = C_nocore[:, :, : self.Nocc].copy() nk_, nao_, nmo_ = C_nocore.shape - s1_val = numpy.zeros((nk_, nmo_, nmo_), dtype=self.S.dtype) - s12_val = numpy.zeros( - (nk_, nmo_, s12_val_.shape[-1]), dtype=s12_val_.dtype - ) - C_nocore_occ = numpy.zeros( + s1_val = zeros((nk_, nmo_, nmo_), dtype=self.S.dtype) + s12_val = zeros((nk_, nmo_, s12_val_.shape[-1]), dtype=s12_val_.dtype) + C_nocore_occ = zeros( (nk_, nao_ - self.ncore, C_nocore_occ_.shape[-1]), dtype=C_nocore_occ_.dtype, ) @@ -203,20 +212,18 @@ def localize( ciao_val_ = get_iao_k( C_nocore_occ, s12_val, s1_val, s2_val, ortho=False ) - Ciao_ = numpy.zeros( - (nk_, nao_, ciao_val_.shape[-1]), dtype=ciao_val_.dtype - ) + Ciao_ = zeros((nk_, nao_, ciao_val_.shape[-1]), dtype=ciao_val_.dtype) for k in range(nk_): Ciao_[k] = C_nocore[k] @ ciao_val_[k] Ciao_[k] = symm_orth_k(Ciao_[k], ovlp=self.S[k]) # stack core|val nao = self.S.shape[-1] - c_core_val = numpy.zeros( + c_core_val = zeros( (nk_, nao, Ciao_.shape[-1] + self.ncore), dtype=Ciao_.dtype ) for k in range(nk_): - c_core_val[k] = numpy.hstack((ciao_core[k], Ciao_[k])) + c_core_val[k] = hstack((ciao_core[k], Ciao_[k])) # tmp - aos are not rearrange and so below is not necessary # (iaoind_by_atom is used to stack iao|pao later) @@ -231,7 +238,7 @@ def localize( c_core_val, self.S, self.cell, valence_basis, self.kpts, ortho=True ) nk, nao, nlo = cpao_.shape - Cpao_ = numpy.zeros((nk, nao, nlo), dtype=numpy.complex128) + Cpao_ = zeros((nk, nao, nlo), dtype=complex128) for k in range(self.nkpt): aoind_by_atom = get_aoind_by_atom(self.cell) ctmp, paoind_by_atom = reorder_by_atom_( @@ -254,7 +261,7 @@ def localize( self.mol, kpts=self.kpts, mo_coeff=Ciao_, mo_energy=mo_energy_ ) - num_wann = numpy.asarray(iaomf.mo_coeff).shape[2] + num_wann = asarray(iaomf.mo_coeff).shape[2] keywords = """ num_iter = 5000 dis_num_iter = 0 @@ -270,33 +277,31 @@ def localize( iaomf, self.kmesh, num_wann, other_keywords=keywords ) - A_matrix = numpy.zeros( - (self.nkpt, num_wann, num_wann), dtype=numpy.complex128 - ) + A_matrix = zeros((self.nkpt, num_wann, num_wann), dtype=complex128) for k in range(self.nkpt): - A_matrix[k] = numpy.eye(num_wann, dtype=numpy.complex128) + A_matrix[k] = eye(num_wann, dtype=complex128) A_matrix = A_matrix.transpose(1, 2, 0) w90.kernel(A_matrix=A_matrix) - u_mat = numpy.array( - w90.U_matrix.transpose(2, 0, 1), order="C", dtype=numpy.complex128 + u_mat = array( + w90.U_matrix.transpose(2, 0, 1), order="C", dtype=complex128 ) os.system("cp wannier90.wout wannier90_iao.wout") os.system("rm wannier90.*") nk, nao, nlo = Ciao_.shape - Ciao = numpy.zeros((nk, nao, nlo), dtype=numpy.complex128) + Ciao = zeros((nk, nao, nlo), dtype=complex128) for k in range(self.nkpt): - Ciao[k] = numpy.dot(Ciao_[k], u_mat[k]) + Ciao[k] = Ciao_[k] @ u_mat[k] # Stack Ciao - Wstack = numpy.zeros( + Wstack = zeros( (self.nkpt, Ciao.shape[1], Ciao.shape[2] + Cpao.shape[2]), - dtype=numpy.complex128, + dtype=complex128, ) if self.frozen_core: for k in range(self.nkpt): @@ -335,18 +340,16 @@ def localize( nlo = self.W.shape[2] nao = self.S.shape[2] - lmo_coeff = numpy.zeros((self.nkpt, nlo, nmo), dtype=numpy.complex128) - cinv_ = numpy.zeros((self.nkpt, nlo, nao), dtype=numpy.complex128) + lmo_coeff = zeros((self.nkpt, nlo, nmo), dtype=complex128) + cinv_ = zeros((self.nkpt, nlo, nao), dtype=complex128) if nmo > nlo: Co_nocore = self.C[:, :, self.ncore : self.Nocc] Cv = self.C[:, :, self.Nocc :] # Ensure that the LOs span the occupied space for k in range(self.nkpt): - assert numpy.allclose( - numpy.sum( - (self.W[k].conj().T @ self.S[k] @ Co_nocore[k]) ** 2.0 - ), + assert allclose( + np.sum((self.W[k].conj().T @ self.S[k] @ Co_nocore[k]) ** 2.0), self.Nocc - self.ncore, ) # Find virtual orbitals that lie in the span of LOs @@ -355,23 +358,21 @@ def localize( ) unused(u) nvlo = nlo - self.Nocc - self.ncore - assert numpy.allclose(numpy.sum(l[:nvlo]), nvlo) - C_ = numpy.hstack([Co_nocore[k], Cv[k] @ vt[:nvlo].conj().T]) + assert allclose(np.sum(l[:nvlo]), nvlo) + C_ = hstack([Co_nocore[k], Cv[k] @ vt[:nvlo].conj().T]) lmo_ = self.W[k].conj().T @ self.S[k] @ C_ - assert numpy.allclose( - lmo_.conj().T @ lmo_, numpy.eye(lmo_.shape[1]) - ) + assert allclose(lmo_.conj().T @ lmo_, eye(lmo_.shape[1])) lmo_coeff.append(lmo_) else: for k in range(self.nkpt): lmo_coeff[k] = multi_dot( (self.W[k].conj().T, self.S[k], self.C[k][:, self.ncore :]), ) - cinv_[k] = numpy.dot(self.W[k].conj().T, self.S[k]) + cinv_[k] = self.W[k].conj().T @ self.S[k] - assert numpy.allclose( + assert allclose( lmo_coeff[k].conj().T @ lmo_coeff[k], - numpy.eye(lmo_coeff[k].shape[1]), + eye(lmo_coeff[k].shape[1]), ) self.lmo_coeff = lmo_coeff @@ -379,16 +380,12 @@ def localize( elif lo_method == "wannier": nk, nao, nmo = self.C.shape - lorb = numpy.zeros((nk, nao, nmo), dtype=numpy.complex128) - lorb_nocore = numpy.zeros( - (nk, nao, nmo - self.ncore), dtype=numpy.complex128 - ) + lorb = zeros((nk, nao, nmo), dtype=complex128) + lorb_nocore = zeros((nk, nao, nmo - self.ncore), dtype=complex128) for k in range(nk): es_, vs_ = eigh(self.S[k]) edx = es_ > 1.0e-14 - lorb[k] = numpy.dot( - vs_[:, edx] / numpy.sqrt(es_[edx]), vs_[:, edx].conj().T - ) + lorb[k] = vs_[:, edx] / sqrt(es_[edx]) @ vs_[:, edx].conj().T if self.frozen_core: Ccore = self.C[k][:, : self.ncore] @@ -428,39 +425,33 @@ def localize( """ w90 = pywannier90.W90(lmf, self.kmesh, num_wann, other_keywords=keywords) - A_matrix = numpy.zeros( - (self.nkpt, num_wann, num_wann), dtype=numpy.complex128 - ) + A_matrix = zeros((self.nkpt, num_wann, num_wann), dtype=complex128) # Using A=I + lowdin orbital and A= + |psi> is the same for k in range(self.nkpt): - A_matrix[k] = numpy.eye(num_wann, dtype=numpy.complex128) + A_matrix[k] = eye(num_wann, dtype=complex128) A_matrix = A_matrix.transpose(1, 2, 0) w90.kernel(A_matrix=A_matrix) - u_mat = numpy.array( - w90.U_matrix.transpose(2, 0, 1), order="C", dtype=numpy.complex128 - ) + u_mat = array(w90.U_matrix.transpose(2, 0, 1), order="C", dtype=complex128) nk, nao, nlo = lmf.mo_coeff.shape - W = numpy.zeros((nk, nao, nlo), dtype=numpy.complex128) + W = zeros((nk, nao, nlo), dtype=complex128) for k in range(nk): - W[k] = numpy.dot(lmf.mo_coeff[k], u_mat[k]) + W[k] = lmf.mo_coeff[k] @ u_mat[k] self.W = W - lmo_coeff = numpy.zeros( - (self.nkpt, nlo, nmo - self.ncore), dtype=numpy.complex128 - ) - cinv_ = numpy.zeros((self.nkpt, nlo, nao), dtype=numpy.complex128) + lmo_coeff = zeros((self.nkpt, nlo, nmo - self.ncore), dtype=complex128) + cinv_ = zeros((self.nkpt, nlo, nao), dtype=complex128) for k in range(nk): lmo_coeff[k] = multi_dot( (self.W[k].conj().T, self.S[k], self.C[k][:, self.ncore :]), ) - cinv_[k] = numpy.dot(self.W[k].conj().T, self.S[k]) - assert numpy.allclose( + cinv_[k] = self.W[k].conj().T @ self.S[k] + assert allclose( lmo_coeff[k].conj().T @ lmo_coeff[k], - numpy.eye(lmo_coeff[k].shape[1]), + eye(lmo_coeff[k].shape[1]), ) self.lmo_coeff = lmo_coeff self.cinv = cinv_ diff --git a/src/quemb/kbe/lo_k.py b/src/quemb/kbe/lo_k.py index 62992f00..b174f9e7 100644 --- a/src/quemb/kbe/lo_k.py +++ b/src/quemb/kbe/lo_k.py @@ -2,8 +2,19 @@ # Oinam Meitei # -import numpy +import numpy as np import scipy +from numpy import ( + allclose, + array, + asarray, + complex128, + diag, + eye, + where, + zeros, + zeros_like, +) from numpy.linalg import eigh, inv, multi_dot, norm from pyscf.pbc import gto as pgto @@ -35,13 +46,13 @@ def cano_orth(A, thr=1.0e-7, ovlp=None): def get_symm_orth_mat_k(A, thr=1.0e-7, ovlp=None): S = dot_gen(A, A, ovlp) e, u = scipy.linalg.eigh(S) - if int(numpy.sum(e < thr)) > 0: + if (e < thr).any(): raise ValueError( "Linear dependence is detected in the column space of A: " "smallest eigenvalue (%.3E) is less than thr (%.3E). " - "Please use 'cano_orth' instead." % (numpy.min(e), thr) + "Please use 'cano_orth' instead." % (np.min(e), thr) ) - return u @ numpy.diag(e**-0.5) @ u.conj().T + return u @ diag(e**-0.5) @ u.conj().T def symm_orth_k(A, thr=1.0e-7, ovlp=None): @@ -69,27 +80,25 @@ def get_xovlp_k(cell, kpts, basis="sto-3g"): cell_alt.basis = basis cell_alt.build() - S22 = numpy.array( - cell_alt.pbc_intor("int1e_ovlp", hermi=1, kpts=kpts), dtype=numpy.complex128 - ) - S12 = numpy.array( + S22 = array(cell_alt.pbc_intor("int1e_ovlp", hermi=1, kpts=kpts), dtype=complex128) + S12 = array( pgto.cell.intor_cross("int1e_ovlp", cell, cell_alt, kpts=kpts), - dtype=numpy.complex128, + dtype=complex128, ) return (S12, S22) def remove_core_mo_k(Clo, Ccore, S, thr=0.5): - assert numpy.allclose(Clo.conj().T @ S @ Clo, numpy.eye(Clo.shape[1])) - assert numpy.allclose(Ccore.conj().T @ S @ Ccore, numpy.eye(Ccore.shape[1])) + assert allclose(Clo.conj().T @ S @ Clo, eye(Clo.shape[1])) + assert allclose(Ccore.conj().T @ S @ Ccore, eye(Ccore.shape[1])) n, nlo = Clo.shape ncore = Ccore.shape[1] Pcore = Ccore @ Ccore.conj().T @ S - Clo1 = (numpy.eye(n) - Pcore) @ Clo - pop = numpy.diag(Clo1.conj().T @ S @ Clo1) - idx_keep = numpy.where(pop > thr)[0] + Clo1 = (eye(n) - Pcore) @ Clo + pop = diag(Clo1.conj().T @ S @ Clo1) + idx_keep = where(pop > thr)[0] assert len(idx_keep) == nlo - ncore Clo2 = symm_orth_k(Clo1[:, idx_keep], ovlp=S) @@ -118,29 +127,29 @@ def get_iao_k(Co, S12, S1, S2=None, ortho=True): nk, nao, nmo = S12.shape unused(nmo) - P1 = numpy.zeros_like(S1, dtype=numpy.complex128) - P2 = numpy.zeros_like(S2, dtype=numpy.complex128) + P1 = zeros_like(S1, dtype=complex128) + P2 = zeros_like(S2, dtype=complex128) for k in range(nk): P1[k] = scipy.linalg.inv(S1[k]) P2[k] = scipy.linalg.inv(S2[k]) - Ciao = numpy.zeros((nk, nao, S12.shape[-1]), dtype=numpy.complex128) + Ciao = zeros((nk, nao, S12.shape[-1]), dtype=complex128) for k in range(nk): # Cotil = P1[k] @ S12[k] @ P2[k] @ S12[k].conj().T @ Co[k] Cotil = multi_dot((P1[k], S12[k], P2[k], S12[k].conj().T, Co[k])) - ptil = numpy.dot(P1[k], S12[k]) + ptil = P1[k] @ S12[k] Stil = multi_dot((Cotil.conj().T, S1[k], Cotil)) - Po = numpy.dot(Co[k], Co[k].conj().T) + Po = Co[k] @ Co[k].conj().T Stil_inv = inv(Stil) Potil = multi_dot((Cotil, Stil_inv, Cotil.conj().T)) Ciao[k] = ( - numpy.eye(nao, dtype=numpy.complex128) - - numpy.dot((Po + Potil - 2.0 * multi_dot((Po, S1[k], Potil))), S1[k]) + eye(nao, dtype=complex128) + - (Po + Potil - 2.0 * multi_dot([Po, S1[k], Potil])) @ S1[k] ) @ ptil if ortho: Ciao[k] = symm_orth_k(Ciao[k], ovlp=S1[k]) @@ -174,13 +183,13 @@ def get_pao_k(Ciao, S, S12): Cpao = [] for k in range(nk): s12 = scipy.linalg.inv(S[k]) @ S12[k] - nonval = numpy.eye(nao) - s12 @ s12.conj().T + nonval = eye(nao) - s12 @ s12.conj().T Piao = Ciao[k] @ Ciao[k].conj().T @ S[k] - cpao_ = (numpy.eye(nao) - Piao) @ nonval + cpao_ = (eye(nao) - Piao) @ nonval Cpao.append(cano_orth(cpao_, ovlp=S[k])) - return numpy.asarray(Cpao) + return asarray(Cpao) def get_pao_native_k(Ciao, S, mol, valence_basis, ortho=True): @@ -219,10 +228,10 @@ def get_pao_native_k(Ciao, S, mol, valence_basis, ortho=True): ] niao = len(vir_idx) - Cpao = numpy.zeros((nk, nao, niao), dtype=numpy.complex128) + Cpao = zeros((nk, nao, niao), dtype=complex128) for k in range(nk): Piao = multi_dot((Ciao[k], Ciao[k].conj().T, S[k])) - cpao_ = (numpy.eye(nao) - Piao)[:, vir_idx] + cpao_ = (eye(nao) - Piao)[:, vir_idx] if ortho: try: Cpao[k] = symm_orth_k(cpao_, ovlp=S[k]) diff --git a/src/quemb/kbe/pbe.py b/src/quemb/kbe/pbe.py index d97f19ef..583cbc64 100644 --- a/src/quemb/kbe/pbe.py +++ b/src/quemb/kbe/pbe.py @@ -5,9 +5,9 @@ from multiprocessing import Pool import h5py -import numpy +import numpy as np from libdmet.basis_transform.eri_transform import get_emb_eri_fast_gdf -from numpy import array, floating +from numpy import array, einsum, floating, result_type, zeros, zeros_like from pyscf import ao2mo, pbc from pyscf.pbc import df, gto from pyscf.pbc.df.df_jk import _ewald_exxdiv_for_G0 @@ -145,7 +145,7 @@ def __init__( self.enuc = mf.energy_nuc() self.hcore = mf.get_hcore() self.S = mf.get_ovlp() - self.C = numpy.array(mf.mo_coeff) + self.C = array(mf.mo_coeff) self.hf_dm = mf.make_rdm1() self.hf_veff = mf.get_veff( self.fobj.mol, @@ -209,21 +209,17 @@ def __init__( nk, nao = self.hf_dm.shape[:2] - dm_nocore = numpy.zeros( - (nk, nao, nao), dtype=numpy.result_type(self.C, self.C) - ) - C_core = numpy.zeros((nk, nao, self.ncore), dtype=self.C.dtype) - P_core = numpy.zeros( - (nk, nao, nao), dtype=numpy.result_type(self.C, self.C) - ) + dm_nocore = zeros((nk, nao, nao), dtype=result_type(self.C, self.C)) + C_core = zeros((nk, nao, self.ncore), dtype=self.C.dtype) + P_core = zeros((nk, nao, nao), dtype=result_type(self.C, self.C)) for k in range(nk): - dm_nocore[k] += 2.0 * numpy.dot( - self.C[k][:, self.ncore : self.ncore + self.Nocc], - self.C[k][:, self.ncore : self.ncore + self.Nocc].conj().T, + dm_nocore[k] += 2.0 * ( + self.C[k][:, self.ncore : self.ncore + self.Nocc] + @ self.C[k][:, self.ncore : self.ncore + self.Nocc].conj().T ) C_core[k] += self.C[k][:, : self.ncore] - P_core[k] += numpy.dot(C_core[k], C_core[k].conj().T) + P_core[k] += C_core[k] @ C_core[k].conj().T self.C_core = C_core self.P_core = P_core @@ -239,23 +235,20 @@ def __init__( ecore_h1 = 0.0 ecore_veff = 0.0 for k in range(nk): - ecore_h1 += numpy.einsum( - "ij,ji", self.hcore[k], 2.0 * self.P_core[k] - ) + ecore_h1 += einsum("ij,ji", self.hcore[k], 2.0 * self.P_core[k]) ecore_veff += ( - numpy.einsum("ij,ji", 2.0 * self.P_core[k], self.core_veff[k]) - * 0.5 + einsum("ij,ji", 2.0 * self.P_core[k], self.core_veff[k]) * 0.5 ) ecore_h1 /= float(nk) ecore_veff /= float(nk) E_core = ecore_h1 + ecore_veff - if numpy.abs(E_core.imag).max() < 1.0e-10: + if np.abs(E_core.imag).max() < 1.0e-10: self.E_core = E_core.real else: raise ValueError( - f"Imaginary density in E_core {numpy.abs(E_core.imag).max()}" + f"Imaginary density in E_core {np.abs(E_core.imag).max()}" ) for k in range(nk): @@ -443,7 +436,7 @@ def ewald_sum(self): dm_ = self.mf.make_rdm1() nk, nao = dm_.shape[:2] - vk_kpts = numpy.zeros(dm_.shape) * 1j + vk_kpts = zeros(dm_.shape) * 1j _ewald_exxdiv_for_G0( self.mf.cell, self.kpts, @@ -451,7 +444,7 @@ def ewald_sum(self): vk_kpts.reshape(-1, nk, nao, nao), kpts_band=self.kpts, ) - e_ = numpy.einsum("kij,kji->", vk_kpts, dm_) * 0.25 + e_ = einsum("kij,kji->", vk_kpts, dm_) * 0.25 e_ /= float(nk) return e_.real @@ -517,7 +510,7 @@ def initialize(self, compute_hf: bool, restart: bool = False) -> None: ) fobjs_.cons_h1(self.hcore) - fobjs_.heff = numpy.zeros_like(fobjs_.h1) + fobjs_.heff = zeros_like(fobjs_.h1) fobjs_.dm_init = fobjs_.get_nsocc( self.S, self.C, self.Nocc, ncore=self.ncore ) @@ -589,11 +582,11 @@ def initialize(self, compute_hf: bool, restart: bool = False) -> None: for frg in range(self.fobj.Nfrag): veff0, veff_ = veffs[frg] - if numpy.abs(veff_.imag).max() < 1.0e-6: + if np.abs(veff_.imag).max() < 1.0e-6: self.Fobjs[frg].veff = veff_.real self.Fobjs[frg].veff0 = veff0.real else: - raise ValueError(f"Imaginary Veff {numpy.abs(veff_.imag).max()}") + raise ValueError(f"Imaginary Veff {np.abs(veff_.imag).max()}") self.Fobjs[frg].fock = self.Fobjs[frg].h1 + veff_.real del veffs @@ -624,12 +617,9 @@ def initialize(self, compute_hf: bool, restart: bool = False) -> None: self.Fobjs[frg]._mo_coeffs = mo_coeffs[frg] for frg in range(self.fobj.Nfrag): - self.Fobjs[frg].dm0 = ( - numpy.dot( - self.Fobjs[frg]._mo_coeffs[:, : self.Fobjs[frg].nsocc], - self.Fobjs[frg]._mo_coeffs[:, : self.Fobjs[frg].nsocc].conj().T, - ) - * 2.0 + self.Fobjs[frg].dm0 = 2.0 * ( + self.Fobjs[frg]._mo_coeffs[:, : self.Fobjs[frg].nsocc] + @ self.Fobjs[frg]._mo_coeffs[:, : self.Fobjs[frg].nsocc].conj().T ) if compute_hf: diff --git a/src/quemb/kbe/pfrag.py b/src/quemb/kbe/pfrag.py index 67035d6b..0b8fc387 100644 --- a/src/quemb/kbe/pfrag.py +++ b/src/quemb/kbe/pfrag.py @@ -2,7 +2,19 @@ import h5py -import numpy +import numpy as np +from numpy import ( + abs, + complex128, + diag_indices, + einsum, + outer, + result_type, + trace, + tril_indices, + zeros, + zeros_like, +) from numpy.linalg import multi_dot from quemb.kbe.helper import get_veff @@ -142,20 +154,18 @@ def sd( Number of k-points in each lattice vector direction """ nk, nao, nlo = lao.shape - rdm1_lo_k = numpy.zeros((nk, nlo, nlo), dtype=numpy.result_type(lmo, lmo)) + rdm1_lo_k = zeros((nk, nlo, nlo), dtype=result_type(lmo, lmo)) for k in range(nk): - rdm1_lo_k[k] += numpy.dot(lmo[k][:, :nocc], lmo[k][:, :nocc].conj().T) + rdm1_lo_k[k] += lmo[k][:, :nocc] @ lmo[k][:, :nocc].conj().T self.rdm1_lo_k = rdm1_lo_k phase = get_phase(cell, kpts, kmesh) - supcell_rdm = numpy.einsum("Rk,kuv,Sk->RuSv", phase, rdm1_lo_k, phase.conj()) + supcell_rdm = einsum("Rk,kuv,Sk->RuSv", phase, rdm1_lo_k, phase.conj()) supcell_rdm = supcell_rdm.reshape(nk * nlo, nk * nlo) - if numpy.abs(supcell_rdm.imag).max() < 1.0e-6: + if (max_val := np.abs(supcell_rdm.imag).max()) < 1.0e-6: supcell_rdm = supcell_rdm.real else: - raise ValueError( - f"Imaginary density in Full SD {numpy.abs(supcell_rdm.imag).max()}" - ) + raise ValueError(f"Imaginary density in Full SD {max_val}") Sites = [i + (nlo * 0) for i in self.fsites] if not frag_type == "autogen": @@ -166,31 +176,27 @@ def sd( TA_R = TA_R.reshape(nk, nlo, teo) phase1 = get_phase1(cell, kpts, kmesh) - TA_k = numpy.einsum("Rim, Rk -> kim", TA_R, phase1) + TA_k = einsum("Rim, Rk -> kim", TA_R, phase1) self.TA_lo_eo = TA_k - TA_ao_eo_k = numpy.zeros( - (nk, nao, teo), dtype=numpy.result_type(lao.dtype, TA_k.dtype) - ) + TA_ao_eo_k = zeros((nk, nao, teo), dtype=result_type(lao.dtype, TA_k.dtype)) for k in range(nk): - TA_ao_eo_k[k] = numpy.dot(lao[k], TA_k[k]) + TA_ao_eo_k[k] = lao[k] @ TA_k[k] self.TA = TA_ao_eo_k self.nao = TA_ao_eo_k.shape[-1] # useful for debugging -- - rdm1_eo = numpy.zeros((teo, teo), dtype=numpy.complex128) + rdm1_eo = zeros((teo, teo), dtype=complex128) for k in range(nk): rdm1_eo += multi_dot((TA_k[k].conj().T, rdm1_lo_k[k], TA_k[k])) rdm1_eo /= float(nk) - h1_eo = numpy.zeros((teo, teo), dtype=numpy.complex128) + h1_eo = zeros((teo, teo), dtype=complex128) for k in range(nk): h1_eo += multi_dot((self.TA[k].conj().T, h1[k], self.TA[k])) h1_eo /= float(nk) - e1 = 2.0 * numpy.einsum( - "ij,ij->i", h1_eo[: self.nfsites], rdm1_eo[: self.nfsites] - ) + e1 = 2.0 * einsum("ij,ij->i", h1_eo[: self.nfsites], rdm1_eo[: self.nfsites]) e_h1 = 0.0 for i in self.efac[1]: e_h1 += self.efac[0] * e1[i] @@ -207,15 +213,15 @@ def cons_h1(self, h1): nk, nao, teo = self.TA.shape unused(nao) - h1_eo = numpy.zeros((teo, teo), dtype=numpy.complex128) + h1_eo = zeros((teo, teo), dtype=complex128) for k in range(nk): h1_eo += multi_dot((self.TA[k].conj().T, h1[k], self.TA[k])) h1_eo /= float(nk) - if numpy.abs(h1_eo.imag).max() < 1.0e-7: + if np.abs(h1_eo.imag).max() < 1.0e-7: self.h1 = h1_eo.real else: - raise ValueError(f"Imaginary Hcore {numpy.abs(h1_eo.imag).max()}") + raise ValueError(f"Imaginary Hcore {np.abs(h1_eo.imag).max()}") def cons_fock(self, hf_veff, S, dm, eri_=None): """ @@ -239,11 +245,11 @@ def cons_fock(self, hf_veff, S, dm, eri_=None): ) veff0, veff_ = get_veff(eri_, dm, S, self.TA, hf_veff, return_veff0=True) - if numpy.abs(veff_.imag).max() < 1.0e-6: + if np.abs(veff_.imag).max() < 1.0e-6: self.veff = veff_.real self.veff0 = veff0.real else: - raise ValueError(f"Imaginary Veff {numpy.abs(veff_.imag).max()}") + raise ValueError(f"Imaginary Veff {abs(veff_.imag).max()}") self.fock = self.h1 + veff_.real @@ -269,26 +275,24 @@ def get_nsocc(self, S, C, nocc, ncore=0): """ nk, nao, neo = self.TA.shape - dm_ = numpy.zeros((nk, nao, nao), dtype=numpy.result_type(C, C)) + dm_ = zeros((nk, nao, nao), dtype=result_type(C, C)) for k in range(nk): - dm_[k] = 2.0 * numpy.dot( - C[k][:, ncore : ncore + nocc], C[k][:, ncore : ncore + nocc].conj().T + dm_[k] = 2.0 * ( + C[k][:, ncore : ncore + nocc] @ C[k][:, ncore : ncore + nocc].conj().T ) - P_ = numpy.zeros((neo, neo), dtype=numpy.complex128) + P_ = zeros((neo, neo), dtype=complex128) for k in range(nk): - Cinv = numpy.dot(self.TA[k].conj().T, S[k]) + Cinv = self.TA[k].conj().T @ S[k] P_ += multi_dot((Cinv, dm_[k], Cinv.conj().T)) P_ /= float(nk) - if numpy.abs(P_.imag).max() < 1.0e-6: + if np.abs(P_.imag).max() < 1.0e-6: P_ = P_.real else: - raise ValueError( - f"Imaginary density in get_nsocc {numpy.abs(P_.imag).max()}" - ) + raise ValueError(f"Imaginary density in get_nsocc {abs(P_.imag).max()}") - nsocc_ = numpy.trace(P_) - nsocc = int(numpy.round(nsocc_.real) / 2) + nsocc_ = trace(P_) + nsocc = int(round(nsocc_.real) / 2) self.nsocc = nsocc return P_ @@ -326,12 +330,9 @@ def scf( eri = get_eri(self.dname, self.nao, eri_file=self.eri_file) if dm0 is None: - dm0 = ( - numpy.dot( - self._mo_coeffs[:, : self.nsocc], - self._mo_coeffs[:, : self.nsocc].conj().T, - ) - * 2.0 + dm0 = 2.0 * ( + self._mo_coeffs[:, : self.nsocc] + @ self._mo_coeffs[:, : self.nsocc].conj().T ) mf_ = get_scfObj( @@ -358,7 +359,7 @@ def update_heff( only_chem=False, ): """Update the effective Hamiltonian for the fragment.""" - heff_ = numpy.zeros_like(self.h1) + heff_ = zeros_like(self.h1) if cout is None: cout = self.udim @@ -398,22 +399,18 @@ def update_ebe_hf( mo_coeffs = self._mo_coeffs if rdm_hf is None: - rdm_hf = numpy.dot( - mo_coeffs[:, : self.nsocc], mo_coeffs[:, : self.nsocc].conj().T - ) + rdm_hf = mo_coeffs[:, : self.nsocc] @ mo_coeffs[:, : self.nsocc].conj().T unrestricted = 1.0 if unrestricted else 2.0 - e1 = unrestricted * numpy.einsum( + e1 = unrestricted * einsum( "ij,ij->i", self.h1[: self.nfsites], rdm_hf[: self.nfsites] ) ec = ( 0.5 * unrestricted - * numpy.einsum( - "ij,ij->i", self.veff[: self.nfsites], rdm_hf[: self.nfsites] - ) + * einsum("ij,ij->i", self.veff[: self.nfsites], rdm_hf[: self.nfsites]) ) if self.TA.ndim == 3: @@ -424,16 +421,16 @@ def update_ebe_hf( with h5py.File(self.eri_file, "r") as r: eri = r[self.dname][()] - e2 = numpy.zeros_like(e1) + e2 = zeros_like(e1) for i in range(self.nfsites): for j in range(jmax): ij = i * (i + 1) // 2 + j if i > j else j * (j + 1) // 2 + i - Gij = (2.0 * rdm_hf[i, j] * rdm_hf - numpy.outer(rdm_hf[i], rdm_hf[j]))[ + Gij = (2.0 * rdm_hf[i, j] * rdm_hf - outer(rdm_hf[i], rdm_hf[j]))[ :jmax, :jmax ] - Gij[numpy.diag_indices(jmax)] *= 0.5 + Gij[diag_indices(jmax)] *= 0.5 Gij += Gij.T - e2[i] += 0.5 * unrestricted * Gij[numpy.tril_indices(jmax)] @ eri[ij] + e2[i] += 0.5 * unrestricted * Gij[tril_indices(jmax)] @ eri[ij] e_ = e1 + e2 + ec etmp = 0.0 diff --git a/src/quemb/kbe/solver.py b/src/quemb/kbe/solver.py index 26f9bd7e..4a384db1 100644 --- a/src/quemb/kbe/solver.py +++ b/src/quemb/kbe/solver.py @@ -1,7 +1,7 @@ # Author(s): Oinam Romesh Meitei -import numpy import scipy.linalg +from numpy import array, complex128, eye, zeros from quemb.shared.helper import unused @@ -31,15 +31,15 @@ def schmidt_decomp_svd(rdm, Frag_sites): Fragsites = [i if i >= 0 else Tot_sites + i for i in Frag_sites] - Env_sites1 = numpy.array([i for i in range(Tot_sites) if i not in Fragsites]) + Env_sites1 = array([i for i in range(Tot_sites) if i not in Fragsites]) nfs = len(Frag_sites) Denv = rdm[Env_sites1][:, Fragsites] U, sigma, V = scipy.linalg.svd(Denv, full_matrices=False, lapack_driver="gesvd") unused(V) nbath = (sigma >= thres).sum() - TA = numpy.zeros((Tot_sites, nfs + nbath), dtype=numpy.complex128) - TA[Fragsites, :nfs] = numpy.eye(nfs) + TA = zeros((Tot_sites, nfs + nbath), dtype=complex128) + TA[Fragsites, :nfs] = eye(nfs) TA[Env_sites1, nfs:] = U[:, :nbath] return TA diff --git a/src/quemb/molbe/be_parallel.py b/src/quemb/molbe/be_parallel.py index e221823a..7dcc9679 100644 --- a/src/quemb/molbe/be_parallel.py +++ b/src/quemb/molbe/be_parallel.py @@ -3,8 +3,7 @@ import os from multiprocessing import Pool -import numpy -from numpy import float64 +from numpy import asarray, diag_indices, einsum, float64, zeros_like from numpy.linalg import multi_dot from pyscf import ao2mo, fci, mcscf @@ -164,7 +163,7 @@ def run_solver( h1_ = multi_dot((mf_.mo_coeff.T, h1, mf_.mo_coeff)) eci, civec = ci_.kernel(h1_, eri, nmo, nelec) unused(eci) - civec = numpy.asarray(civec) + civec = asarray(civec) (rdm1a_, rdm1b_), (rdm2aa, rdm2ab, rdm2bb) = ci_.make_rdm12s(civec, nmo, nelec) rdm1_tmp = rdm1a_ + rdm1b_ @@ -236,19 +235,19 @@ def run_solver( elif solver == "FCI": rdm2s = mc_.make_rdm2(civec, mc_.norb, mc_.nelec) if use_cumulant: - hf_dm = numpy.zeros_like(rdm1_tmp) - hf_dm[numpy.diag_indices(nocc)] += 2.0 + hf_dm = zeros_like(rdm1_tmp) + hf_dm[diag_indices(nocc)] += 2.0 del_rdm1 = rdm1_tmp.copy() - del_rdm1[numpy.diag_indices(nocc)] -= 2.0 + del_rdm1[diag_indices(nocc)] -= 2.0 nc = ( - numpy.einsum("ij,kl->ijkl", hf_dm, hf_dm) - + numpy.einsum("ij,kl->ijkl", hf_dm, del_rdm1) - + numpy.einsum("ij,kl->ijkl", del_rdm1, hf_dm) + einsum("ij,kl->ijkl", hf_dm, hf_dm) + + einsum("ij,kl->ijkl", hf_dm, del_rdm1) + + einsum("ij,kl->ijkl", del_rdm1, hf_dm) ) nc -= ( - numpy.einsum("ij,kl->iklj", hf_dm, hf_dm) - + numpy.einsum("ij,kl->iklj", hf_dm, del_rdm1) - + numpy.einsum("ij,kl->iklj", del_rdm1, hf_dm) + einsum("ij,kl->iklj", hf_dm, hf_dm) + + einsum("ij,kl->iklj", hf_dm, del_rdm1) + + einsum("ij,kl->iklj", del_rdm1, hf_dm) ) * 0.5 rdm2s -= nc e_f = get_frag_energy( diff --git a/src/quemb/molbe/eri_onthefly.py b/src/quemb/molbe/eri_onthefly.py index d169f9d7..364fba76 100644 --- a/src/quemb/molbe/eri_onthefly.py +++ b/src/quemb/molbe/eri_onthefly.py @@ -1,6 +1,6 @@ # Author(s): Minsik Cho, Hong-Zhou Ye -import numpy +from numpy import moveaxis, transpose, zeros from pyscf import lib from pyscf.ao2mo.addons import restore from pyscf.df.addons import make_auxmol @@ -106,7 +106,7 @@ def block_step_size(nfrag, naux, nao): j2c = auxmol.intor(mf.mol._add_suffix("int2c2e"), hermi=1) # (L|M) low = cholesky(j2c, lower=True) pqL_frag = [ - numpy.zeros((auxmol.nao, fragobj.nao, fragobj.nao)) for fragobj in Fobjs + zeros((auxmol.nao, fragobj.nao, fragobj.nao)) for fragobj in Fobjs ] # place to store fragment (pq|L) end = 0 atm, bas, env = mole.conc_env( @@ -132,9 +132,9 @@ def block_step_size(nfrag, naux, nao): for fragidx in range(len(Fobjs)): if settings.PRINT_LEVEL > 10: print("(μν|P) -> (ij|P) for frag #", fragidx, flush=True) - Lqp = numpy.transpose(ints, axes=(2, 1, 0)) + Lqp = transpose(ints, axes=(2, 1, 0)) Lqi = Lqp @ Fobjs[fragidx].TA - Liq = numpy.moveaxis(Lqi, 2, 1) + Liq = moveaxis(Lqi, 2, 1) pqL_frag[fragidx][start:end, :, :] = Liq @ Fobjs[fragidx].TA # Fit to get B_{ij}^{L} for fragidx in range(len(Fobjs)): diff --git a/src/quemb/molbe/helper.py b/src/quemb/molbe/helper.py index 1f8a6b9d..752058cc 100644 --- a/src/quemb/molbe/helper.py +++ b/src/quemb/molbe/helper.py @@ -3,7 +3,17 @@ import h5py -import numpy +from numpy import ( + array, + asarray, + diag_indices, + einsum, + eye, + float64, + tril_indices, + zeros, + zeros_like, +) from numpy.linalg import multi_dot from pyscf import ao2mo, gto, lib, scf @@ -38,12 +48,12 @@ def get_veff(eri_, dm, S, TA, hf_veff): """ # Transform the density matrix - ST = numpy.dot(S, TA) + ST = S @ TA P_ = multi_dot((ST.T, dm, ST)) # Ensure the transformed density matrix and ERI are real and double-precision - P_ = numpy.asarray(P_.real, dtype=numpy.double) - eri_ = numpy.asarray(eri_, dtype=numpy.double) + P_ = asarray(P_.real, dtype=float64) + eri_ = asarray(eri_, dtype=float64) # Compute the Coulomb (J) and exchange (K) integrals vj, vk = scf.hf.dot_eri_dm(eri_, P_, hermi=1, with_j=True, with_k=True) @@ -89,7 +99,7 @@ def get_scfObj( nao = h1.shape[0] # Initialize a dummy molecule with the required number of electrons - S = numpy.eye(nao) + S = eye(nao) mol = gto.M() mol.nelectron = nocc * 2 mol.incore_anyway = True @@ -165,7 +175,7 @@ def get_eri(i_frag, Nao, symm=8, ignore_symm=False, eri_file="eri_file.h5"): """ # Open the HDF5 file and read the ERI for the specified fragment with h5py.File(eri_file, "r") as r: - eri__ = numpy.array(r.get(i_frag)) + eri__ = array(r.get(i_frag)) # Optionally restore the symmetry of the ERI if not ignore_symm: @@ -264,20 +274,20 @@ def get_frag_energy( rdm1s_rot = mo_coeffs @ rdm1 @ mo_coeffs.T * 0.5 # Construct the Hartree-Fock 1-RDM - hf_1rdm = numpy.dot(mo_coeffs[:, :nsocc], mo_coeffs[:, :nsocc].conj().T) + hf_1rdm = mo_coeffs[:, :nsocc] @ mo_coeffs[:, :nsocc].conj().T if use_cumulant: # Compute the difference between the rotated RDM1 and the Hartree-Fock 1-RDM delta_rdm1 = 2 * (rdm1s_rot - hf_1rdm) # Calculate the one-electron contributions - e1 = numpy.einsum("ij,ij->i", h1[:nfsites], delta_rdm1[:nfsites]) - ec = numpy.einsum("ij,ij->i", veff0[:nfsites], delta_rdm1[:nfsites]) + e1 = einsum("ij,ij->i", h1[:nfsites], delta_rdm1[:nfsites]) + ec = einsum("ij,ij->i", veff0[:nfsites], delta_rdm1[:nfsites]) else: # Calculate the one-electron and effective potential energy contributions - e1 = 2 * numpy.einsum("ij,ij->i", h1[:nfsites], rdm1s_rot[:nfsites]) - ec = numpy.einsum("ij,ij->i", veff[:nfsites], rdm1s_rot[:nfsites]) + e1 = 2 * einsum("ij,ij->i", h1[:nfsites], rdm1s_rot[:nfsites]) + ec = einsum("ij,ij->i", veff[:nfsites], rdm1s_rot[:nfsites]) if TA.ndim == 3: jmax = TA[0].shape[1] @@ -289,21 +299,21 @@ def get_frag_energy( eri = r[dname][()] # Rotate the RDM2 into the MO basis - rdm2s = numpy.einsum( + rdm2s = einsum( "ijkl,pi,qj,rk,sl->pqrs", 0.5 * rdm2s, *([mo_coeffs] * 4), optimize=True ) # Initialize the two-electron energy contribution - e2 = numpy.zeros_like(e1) + e2 = zeros_like(e1) # Calculate the two-electron energy contribution for i in range(nfsites): for j in range(jmax): ij = i * (i + 1) // 2 + j if i > j else j * (j + 1) // 2 + i Gij = rdm2s[i, j, :jmax, :jmax].copy() - Gij[numpy.diag_indices(jmax)] *= 0.5 + Gij[diag_indices(jmax)] *= 0.5 Gij += Gij.T - e2[i] += Gij[numpy.tril_indices(jmax)] @ eri[ij] + e2[i] += Gij[tril_indices(jmax)] @ eri[ij] # Sum the energy contributions e_ = e1 + e2 + ec @@ -388,7 +398,7 @@ def get_frag_energy_u( # Construct the Hartree-Fock RDM1 for both spin the the Schmidt space hf_1rdm = [ - numpy.dot(mo_coeffs[s][:, : nsocc[s]], mo_coeffs[s][:, : nsocc[s]].conj().T) + mo_coeffs[s][:, : nsocc[s]] @ mo_coeffs[s][:, : nsocc[s]].conj().T for s in [0, 1] ] @@ -407,11 +417,11 @@ def get_frag_energy_u( # Calculate the one-electron and effective potential energy contributions e1 = [ - numpy.einsum("ij,ij->i", h1[s][: nfsites[s]], delta_rdm1[s][: nfsites[s]]) + einsum("ij,ij->i", h1[s][: nfsites[s]], delta_rdm1[s][: nfsites[s]]) for s in [0, 1] ] ec = [ - numpy.einsum("ij,ij->i", veff0[s][: nfsites[s]], delta_rdm1[s][: nfsites[s]]) + einsum("ij,ij->i", veff0[s][: nfsites[s]], delta_rdm1[s][: nfsites[s]]) for s in [0, 1] ] @@ -423,7 +433,7 @@ def get_frag_energy_u( # Rotate the RDM2 into the MO basis rdm2s_k = [ - numpy.einsum( + einsum( "ijkl,pi,qj,rk,sl->pqrs", rdm2s[s], *([mo_coeffs[s12[0]]] * 2 + [mo_coeffs[s12[1]]] * 2), @@ -433,11 +443,11 @@ def get_frag_energy_u( ] # Initialize the two-electron energy contribution - e2 = [numpy.zeros(h1[0].shape[0]), numpy.zeros(h1[1].shape[0])] + e2 = [zeros(h1[0].shape[0]), zeros(h1[1].shape[0])] # Calculate the two-electron energy contribution for alpha and beta def contract_2e(jmaxs, rdm2_, V_, s, sym): - e2_ = numpy.zeros(nfsites[s]) + e2_ = zeros(nfsites[s]) jmax1, jmax2 = [jmaxs] * 2 if isinstance(jmaxs, int) else jmaxs for i in range(nfsites[s]): for j in range(jmax1): @@ -448,9 +458,9 @@ def contract_2e(jmaxs, rdm2_, V_, s, sym): else: Gij = rdm2_[:jmax2, :jmax2, i, j].copy() Vij = V_[:, ij] - Gij[numpy.diag_indices(jmax2)] *= 0.5 + Gij[diag_indices(jmax2)] *= 0.5 Gij += Gij.T - e2_[i] += Gij[numpy.tril_indices(jmax2)] @ Vij + e2_[i] += Gij[tril_indices(jmax2)] @ Vij e2_ *= 0.5 return e2_ diff --git a/src/quemb/molbe/lo.py b/src/quemb/molbe/lo.py index 44ec2213..6ab763e0 100644 --- a/src/quemb/molbe/lo.py +++ b/src/quemb/molbe/lo.py @@ -1,8 +1,8 @@ # Author(s): Henry Tran, Oinam Meitei, Shaun Weatherly # -import numpy -from numpy import allclose, eye +import numpy as np +from numpy import allclose, diag, eye, sqrt, where, zeros from numpy.linalg import eigh, inv, multi_dot, norm, svd from pyscf.gto import intor_cross from pyscf.gto.mole import Mole @@ -45,9 +45,9 @@ def get_symm_orth_mat( raise ValueError( "Linear dependence is detected in the column space of A: " "smallest eigenvalue (%.3E) is less than thr (%.3E). " - "Please use 'cano_orth' instead." % (numpy.min(e), thr) + "Please use 'cano_orth' instead." % (np.min(e), thr) ) - return u @ numpy.diag(e**-0.5) @ u.T + return u @ diag(e**-0.5) @ u.T def symm_orth(A: Matrix, thr: float = 1.0e-6, ovlp: Matrix | None = None) -> Matrix: @@ -63,8 +63,8 @@ def remove_core_mo(Clo: Matrix, Ccore: Matrix, S: Matrix, thr: float = 0.5) -> M ncore = Ccore.shape[1] Pcore = Ccore @ Ccore.T @ S Clo1 = (eye(n) - Pcore) @ Clo - pop = numpy.diag(Clo1.T @ S @ Clo1) - idx_keep = numpy.where(pop > thr)[0] + pop = diag(Clo1.T @ S @ Clo1) + idx_keep = where(pop > thr)[0] assert len(idx_keep) == nlo - ncore return symm_orth(Clo1[:, idx_keep], ovlp=S) @@ -283,41 +283,40 @@ def localize( if lo_method == "lowdin": es_, vs_ = eigh(self.S) edx = es_ > 1.0e-15 - self.W = numpy.dot(vs_[:, edx] / numpy.sqrt(es_[edx]), vs_[:, edx].T) + self.W = vs_[:, edx] / sqrt(es_[edx]) @ vs_[:, edx].T if self.frozen_core: if self.unrestricted: P_core = [ - eye(self.W.shape[0]) - numpy.dot(self.P_core[s], self.S) - for s in [0, 1] + eye(self.W.shape[0]) - (self.P_core[s] @ self.S) for s in [0, 1] ] - C_ = numpy.dot(P_core, self.W) + C_ = P_core @ self.W Cpop = [multi_dot((C_[s].T, self.S, C_[s])) for s in [0, 1]] - Cpop = [numpy.diag(Cpop[s]) for s in [0, 1]] - no_core_idx = [numpy.where(Cpop[s] > 0.7)[0] for s in [0, 1]] + Cpop = [diag(Cpop[s]) for s in [0, 1]] + no_core_idx = [where(Cpop[s] > 0.7)[0] for s in [0, 1]] C_ = [C_[s][:, no_core_idx[s]] for s in [0, 1]] S_ = [multi_dot((C_[s].T, self.S, C_[s])) for s in [0, 1]] W_ = [] for s in [0, 1]: es_, vs_ = eigh(S_[s]) - s_ = numpy.sqrt(es_) - s_ = numpy.diag(1.0 / s_) + s_ = sqrt(es_) + s_ = diag(1.0 / s_) W_.append(multi_dot((vs_, s_, vs_.T))) - self.W = [numpy.dot(C_[s], W_[s]) for s in [0, 1]] + self.W = [C_[s] @ W_[s] for s in [0, 1]] else: - P_core = eye(self.W.shape[0]) - numpy.dot(self.P_core, self.S) - C_ = numpy.dot(P_core, self.W) + P_core = eye(self.W.shape[0]) - self.P_core @ self.S + C_ = P_core @ self.W # NOTE: PYSCF has basis in 1s2s3s2p2p2p3p3p3p format # fix no_core_idx - use population for now Cpop = multi_dot((C_.T, self.S, C_)) - Cpop = numpy.diag(Cpop) - no_core_idx = numpy.where(Cpop > 0.7)[0] + Cpop = diag(Cpop) + no_core_idx = where(Cpop > 0.7)[0] C_ = C_[:, no_core_idx] S_ = multi_dot((C_.T, self.S, C_)) es_, vs_ = eigh(S_) - s_ = numpy.sqrt(es_) - s_ = numpy.diag(1.0 / s_) + s_ = sqrt(es_) + s_ = diag(1.0 / s_) W_ = multi_dot((vs_, s_, vs_.T)) - self.W = numpy.dot(C_, W_) + self.W = C_ @ W_ if self.unrestricted: if self.frozen_core: @@ -341,24 +340,24 @@ def localize( elif lo_method in ["pipek-mezey", "pipek", "PM"]: es_, vs_ = eigh(self.S) edx = es_ > 1.0e-15 - self.W = numpy.dot(vs_[:, edx] / numpy.sqrt(es_[edx]), vs_[:, edx].T) + self.W = vs_[:, edx] / sqrt(es_[edx]) @ vs_[:, edx].T es_, vs_ = eigh(self.S) edx = es_ > 1.0e-15 - W_ = numpy.dot(vs_[:, edx] / numpy.sqrt(es_[edx]), vs_[:, edx].T) + W_ = vs_[:, edx] / sqrt(es_[edx]) @ vs_[:, edx].T if self.frozen_core: - P_core = eye(W_.shape[0]) - numpy.dot(self.P_core, self.S) - C_ = numpy.dot(P_core, W_) + P_core = eye(W_.shape[0]) - self.P_core @ self.S + C_ = P_core @ W_ Cpop = multi_dot((C_.T, self.S, C_)) - Cpop = numpy.diag(Cpop) - no_core_idx = numpy.where(Cpop > 0.55)[0] + Cpop = diag(Cpop) + no_core_idx = where(Cpop > 0.55)[0] C_ = C_[:, no_core_idx] S_ = multi_dot((C_.T, self.S, C_)) es_, vs_ = eigh(S_) - s_ = numpy.sqrt(es_) - s_ = numpy.diag(1.0 / s_) + s_ = sqrt(es_) + s_ = diag(1.0 / s_) W_ = multi_dot((vs_, s_, vs_.T)) - W_ = numpy.dot(C_, W_) + W_ = C_ @ W_ self.W = get_loc( self.mol, W_, "PM", pop_method=pop_method, init_guess=init_guess @@ -409,11 +408,11 @@ def localize( shift = 0 ncore = 0 if not valence_only: - Wstack = numpy.zeros( + Wstack = zeros( (Ciao.shape[0], Ciao.shape[1] + Cpao.shape[1]) ) # -self.ncore)) else: - Wstack = numpy.zeros((Ciao.shape[0], Ciao.shape[1])) + Wstack = zeros((Ciao.shape[0], Ciao.shape[1])) if self.frozen_core: for ix in range(self.mol.natm): @@ -440,7 +439,7 @@ def localize( ] shift += npao else: - Wstack = numpy.hstack((Ciao, Cpao)) + Wstack = hstack((Ciao, Cpao)) if not nosave: self.W = Wstack assert allclose(self.W.T @ self.S @ self.W, eye(self.W.shape[1])) @@ -456,15 +455,15 @@ def localize( Cv = self.C[:, self.Nocc :] # Ensure that the LOs span the occupied space assert allclose( - numpy.sum((self.W.T @ self.S @ Co_nocore) ** 2.0), + np.sum((self.W.T @ self.S @ Co_nocore) ** 2.0), self.Nocc - self.ncore, ) # Find virtual orbitals that lie in the span of LOs u, l, vt = svd(self.W.T @ self.S @ Cv, full_matrices=False) unused(u) nvlo = nlo - self.Nocc - self.ncore - assert allclose(numpy.sum(l[:nvlo]), nvlo) - C_ = numpy.hstack([Co_nocore, Cv @ vt[:nvlo].T]) + assert allclose(np.sum(l[:nvlo]), nvlo) + C_ = hstack([Co_nocore, Cv @ vt[:nvlo].T]) self.lmo_coeff = self.W.T @ self.S @ C_ else: self.lmo_coeff = self.W.T @ self.S @ self.C[:, self.ncore :] @@ -474,20 +473,20 @@ def localize( elif lo_method == "boys": es_, vs_ = eigh(self.S) edx = es_ > 1.0e-15 - W_ = numpy.dot(vs_[:, edx] / numpy.sqrt(es_[edx]), vs_[:, edx].T) + W_ = vs_[:, edx] / sqrt(es_[edx]) @ vs_[:, edx].T if self.frozen_core: - P_core = eye(W_.shape[0]) - numpy.dot(self.P_core, self.S) - C_ = numpy.dot(P_core, W_) + P_core = eye(W_.shape[0]) - self.P_core @ self.S + C_ = P_core @ W_ Cpop = multi_dot((C_.T, self.S, C_)) - Cpop = numpy.diag(Cpop) - no_core_idx = numpy.where(Cpop > 0.55)[0] + Cpop = diag(Cpop) + no_core_idx = where(Cpop > 0.55)[0] C_ = C_[:, no_core_idx] S_ = multi_dot((C_.T, self.S, C_)) es_, vs_ = eigh(S_) - s_ = numpy.sqrt(es_) - s_ = numpy.diag(1.0 / s_) + s_ = sqrt(es_) + s_ = diag(1.0 / s_) W_ = multi_dot((vs_, s_, vs_.T)) - W_ = numpy.dot(C_, W_) + W_ = C_ @ W_ self.W = get_loc(self.mol, W_, "BOYS") diff --git a/src/quemb/molbe/mbe.py b/src/quemb/molbe/mbe.py index ffd05f12..73c77b30 100644 --- a/src/quemb/molbe/mbe.py +++ b/src/quemb/molbe/mbe.py @@ -5,7 +5,7 @@ import h5py import numpy from attrs import define -from numpy import floating +from numpy import array, diag_indices, einsum, float64, floating, zeros, zeros_like from pyscf import ao2mo, scf from quemb.molbe.be_parallel import be_func_parallel @@ -159,7 +159,7 @@ def __init__( self.hcore = mf.get_hcore() self.S = mf.get_ovlp() - self.C = numpy.array(mf.mo_coeff) + self.C = array(mf.mo_coeff) self.hf_dm = mf.make_rdm1() self.hf_veff = mf.get_veff() self.hf_etot = mf.e_tot @@ -193,14 +193,14 @@ def __init__( if not restart: self.Nocc -= self.ncore - self.hf_dm = 2.0 * numpy.dot( - self.C[:, self.ncore : self.ncore + self.Nocc], - self.C[:, self.ncore : self.ncore + self.Nocc].T, + self.hf_dm = 2.0 * ( + self.C[:, self.ncore : self.ncore + self.Nocc] + @ self.C[:, self.ncore : self.ncore + self.Nocc].T ) self.C_core = self.C[:, : self.ncore] - self.P_core = numpy.dot(self.C_core, self.C_core.T) + self.P_core = self.C_core @ self.C_core.T self.core_veff = mf.get_veff(dm=self.P_core * 2.0) - self.E_core = numpy.einsum( + self.E_core = einsum( "ji,ji->", 2.0 * self.hcore + self.core_veff, self.P_core ) self.hf_veff -= self.core_veff @@ -311,20 +311,20 @@ def rdm1_fullbasis( nao = C_mo.shape[0] # Initialize density matrices for atomic orbitals (AO) - rdm1AO = numpy.zeros((nao, nao)) - rdm2AO = numpy.zeros((nao, nao, nao, nao)) + rdm1AO = zeros((nao, nao)) + rdm2AO = zeros((nao, nao, nao, nao)) for fobjs in self.Fobjs: if return_RDM2: # Adjust the one-particle reduced density matrix (RDM1) drdm1 = fobjs.rdm1__.copy() - drdm1[numpy.diag_indices(fobjs.nsocc)] -= 2.0 + drdm1[diag_indices(fobjs.nsocc)] -= 2.0 # Compute the two-particle reduced density matrix (RDM2) and subtract # non-connected component - dm_nc = numpy.einsum( + dm_nc = einsum( "ij,kl->ijkl", drdm1, drdm1, dtype=numpy.float64, optimize=True - ) - 0.5 * numpy.einsum( + ) - 0.5 * einsum( "ij,kl->iklj", drdm1, drdm1, dtype=numpy.float64, optimize=True ) fobjs.rdm2__ -= dm_nc @@ -350,13 +350,13 @@ def rdm1_fullbasis( if not only_rdm1: # Transform RDM2 to AO basis - rdm2s = numpy.einsum( + rdm2s = einsum( "ijkl,pi,qj,rk,sl->pqrs", fobjs.rdm2__, *([fobjs.mo_coeffs] * 4), optimize=True, ) - rdm2_ao = numpy.einsum( + rdm2_ao = einsum( "xi,ijkl,px,qj,rk,sl->pqrs", Pc_, rdm2s, @@ -373,14 +373,14 @@ def rdm1_fullbasis( rdm2AO = (rdm2AO + rdm2AO.T) / 2.0 if return_RDM2: nc_AO = ( - numpy.einsum( + einsum( "ij,kl->ijkl", rdm1AO, rdm1AO, dtype=numpy.float64, optimize=True, ) - - numpy.einsum( + - einsum( "ij,kl->iklj", rdm1AO, rdm1AO, @@ -394,7 +394,7 @@ def rdm1_fullbasis( # Transform RDM2 to the molecular orbital (MO) basis if needed if not return_ao: CmoT_S = self.C.T @ self.S - rdm2MO = numpy.einsum( + rdm2MO = einsum( "ijkl,pi,qj,rk,sl->pqrs", rdm2AO, CmoT_S, @@ -407,7 +407,7 @@ def rdm1_fullbasis( # Transform RDM2 to the localized orbital (LO) basis if needed if return_lo: CloT_S = self.W.T @ self.S - rdm2LO = numpy.einsum( + rdm2LO = einsum( "ijkl,pi,qj,rk,sl->pqrs", rdm2AO, CloT_S, @@ -431,9 +431,9 @@ def rdm1_fullbasis( if return_RDM2 and print_energy: # Compute and print energy contributions - Eh1 = numpy.einsum("ij,ij", self.hcore, rdm1AO, optimize=True) + Eh1 = einsum("ij,ij", self.hcore, rdm1AO, optimize=True) eri = ao2mo.restore(1, self.mf._eri, self.mf.mo_coeff.shape[1]) - E2 = 0.5 * numpy.einsum("pqrs,pqrs", eri, rdm2AO, optimize=True) + E2 = 0.5 * einsum("pqrs,pqrs", eri, rdm2AO, optimize=True) print(flush=True) print("-----------------------------------------------------", flush=True) print(" BE ENERGIES with cumulant-based expression", flush=True) @@ -517,12 +517,8 @@ def compute_energy_full( if return_rdm: # Construct the full RDM2 from RDM1 RDM2_full = ( - numpy.einsum( - "ij,kl->ijkl", rdm1f, rdm1f, dtype=numpy.float64, optimize=True - ) - - numpy.einsum( - "ij,kl->iklj", rdm1f, rdm1f, dtype=numpy.float64, optimize=True - ) + einsum("ij,kl->ijkl", rdm1f, rdm1f, dtype=float64, optimize=True) + - einsum("ij,kl->iklj", rdm1f, rdm1f, dtype=float64, optimize=True) * 0.5 ) @@ -539,30 +535,30 @@ def compute_energy_full( veff = scf.hf.get_veff(self.fobj.mol, rdm1f, hermi=0) # Compute the one-electron energy - Eh1 = numpy.einsum("ij,ij", self.hcore, rdm1f, optimize=True) + Eh1 = einsum("ij,ij", self.hcore, rdm1f, optimize=True) # Compute the energy due to the effective potential - EVeff = numpy.einsum("ij,ij", veff, rdm1f, optimize=True) + EVeff = einsum("ij,ij", veff, rdm1f, optimize=True) # Compute the change in the one-electron energy - Eh1_dg = numpy.einsum("ij,ij", self.hcore, del_gamma, optimize=True) + Eh1_dg = einsum("ij,ij", self.hcore, del_gamma, optimize=True) # Compute the change in the effective potential energy - Eveff_dg = numpy.einsum("ij,ij", self.hf_veff, del_gamma, optimize=True) + Eveff_dg = einsum("ij,ij", self.hf_veff, del_gamma, optimize=True) # Restore the electron repulsion integrals (ERI) eri = ao2mo.restore(1, self.mf._eri, self.mf.mo_coeff.shape[1]) # Compute the cumulant part of the two-electron energy - EKumul = numpy.einsum("pqrs,pqrs", eri, Kumul, optimize=True) + EKumul = einsum("pqrs,pqrs", eri, Kumul, optimize=True) if not approx_cumulant: # Compute the true two-electron energy if not using approximate cumulant - EKumul_T = numpy.einsum("pqrs,pqrs", eri, Kumul_T, optimize=True) + EKumul_T = einsum("pqrs,pqrs", eri, Kumul_T, optimize=True) if use_full_rdm and return_rdm: # Compute the full two-electron energy using the full RDM2 - E2 = numpy.einsum("pqrs,pqrs", eri, RDM2_full, optimize=True) + E2 = einsum("pqrs,pqrs", eri, RDM2_full, optimize=True) # Compute the approximate BE total energy EKapprox = self.ebe_hf + Eh1_dg + Eveff_dg + EKumul / 2.0 @@ -699,7 +695,7 @@ def optimize( if method == "QN": # Prepare the initial Jacobian matrix if only_chem: - J0 = numpy.array([[0.0]]) + J0 = array([[0.0]]) J0 = self.get_be_error_jacobian(jac_solver="HF") J0 = J0[-1:, -1:] else: @@ -847,7 +843,7 @@ def initialize(self, eri_, compute_hf, restart=False) -> None: for fobjs_ in self.Fobjs: # Process each fragment - eri = numpy.array(file_eri.get(fobjs_.dname)) + eri = array(file_eri.get(fobjs_.dname)) _ = fobjs_.get_nsocc(self.S, self.C, self.Nocc, ncore=self.ncore) fobjs_.cons_h1(self.hcore) @@ -857,7 +853,7 @@ def initialize(self, eri_, compute_hf, restart=False) -> None: fobjs_.cons_fock(self.hf_veff, self.S, self.hf_dm, eri_=eri) - fobjs_.heff = numpy.zeros_like(fobjs_.h1) + fobjs_.heff = zeros_like(fobjs_.h1) fobjs_.scf(fs=True, eri=eri) fobjs_.dm0 = 2.0 * ( diff --git a/src/quemb/molbe/misc.py b/src/quemb/molbe/misc.py index 3167d96c..aad0eefc 100644 --- a/src/quemb/molbe/misc.py +++ b/src/quemb/molbe/misc.py @@ -4,7 +4,7 @@ import time import h5py -import numpy +from numpy import einsum, ix_, loadtxt from pyscf import ao2mo, df, gto, qmmm, scf from pyscf.lib import chkfile from pyscf.tools import fcidump @@ -74,7 +74,7 @@ def libint2pyscf( raise ValueError("Input core Hamiltonian file does not exist") mol = gto.M(atom=xyzfile, basis=basis, spin=spin, charge=charge) - hcore_libint = numpy.loadtxt(hcore, skiprows=hcore_skiprows) + hcore_libint = loadtxt(hcore, skiprows=hcore_skiprows) libint2pyscf = [] for labelidx, label in enumerate(mol.ao_labels()): @@ -90,7 +90,7 @@ def libint2pyscf( elif "z" in label.split()[2]: libint2pyscf.append(labelidx - 1) - hcore_pyscf = hcore_libint[numpy.ix_(libint2pyscf, libint2pyscf)] + hcore_pyscf = hcore_libint[ix_(libint2pyscf, libint2pyscf)] mol.incore_anyway = True if use_df: @@ -130,10 +130,10 @@ def be2fcidump(be_obj, fcidump_prefix, basis): h2e = eri elif basis == "fragment_mo": frag.scf() # make sure that we have mo coefficients - h1e = numpy.einsum( + h1e = einsum( "ij,ia,jb->ab", frag.fock, frag.mo_coeffs, frag.mo_coeffs, optimize=True ) - h2e = numpy.einsum( + h2e = einsum( "ijkl,ia,jb,kc,ld->abcd", eri, frag.mo_coeffs, @@ -180,10 +180,10 @@ def ube2fcidump(be_obj, fcidump_prefix, basis): h2e = eri elif basis == "fragment_mo": frag.scf() # make sure that we have mo coefficients - h1e = numpy.einsum( + h1e = einsum( "ij,ia,jb->ab", frag.fock, frag.mo_coeffs, frag.mo_coeffs, optimize=True ) - h2e = numpy.einsum( + h2e = einsum( "ijkl,ia,jb,kc,ld->abcd", eri, frag.mo_coeffs, @@ -214,10 +214,10 @@ def ube2fcidump(be_obj, fcidump_prefix, basis): h2e = eri elif basis == "fragment_mo": frag.scf() # make sure that we have mo coefficients - h1e = numpy.einsum( + h1e = einsum( "ij,ia,jb->ab", frag.fock, frag.mo_coeffs, frag.mo_coeffs, optimize=True ) - h2e = numpy.einsum( + h2e = einsum( "ijkl,ia,jb,kc,ld->abcd", eri, frag.mo_coeffs, @@ -343,18 +343,14 @@ def be2puffin( elif "z" in label.split()[2]: libint2pyscf.append(labelidx - 1) - hcore_pyscf = hcore[numpy.ix_(libint2pyscf, libint2pyscf)] + hcore_pyscf = hcore[ix_(libint2pyscf, libint2pyscf)] else: # Input hcore is in PySCF format hcore_pyscf = hcore if jk is not None: jk_pyscf = ( - jk[0][ - numpy.ix_(libint2pyscf, libint2pyscf, libint2pyscf, libint2pyscf) - ], - jk[1][ - numpy.ix_(libint2pyscf, libint2pyscf, libint2pyscf, libint2pyscf) - ], + jk[0][ix_(libint2pyscf, libint2pyscf, libint2pyscf, libint2pyscf)], + jk[1][ix_(libint2pyscf, libint2pyscf, libint2pyscf, libint2pyscf)], ) mol.incore_anyway = True diff --git a/src/quemb/molbe/opt.py b/src/quemb/molbe/opt.py index dea701f0..cfda6f99 100644 --- a/src/quemb/molbe/opt.py +++ b/src/quemb/molbe/opt.py @@ -1,7 +1,6 @@ # Author(s): Oinam Romesh Meitei -import numpy from attrs import Factory, define from numpy import array, float64 @@ -172,7 +171,7 @@ def optimize(self, method, J0=None, trust_region=False): # Initialize the Quasi-Newton optimizer optQN = FrankQN( - self.objfunc, numpy.array(self.pot), f0, J0, max_space=self.max_space + self.objfunc, array(self.pot), f0, J0, max_space=self.max_space ) if self.err < self.conv_tol: diff --git a/src/quemb/molbe/pfrag.py b/src/quemb/molbe/pfrag.py index 7561a926..3d6767bb 100644 --- a/src/quemb/molbe/pfrag.py +++ b/src/quemb/molbe/pfrag.py @@ -1,9 +1,21 @@ # Author(s): Oinam Romesh Meitei import h5py -import numpy +import numpy as np import scipy.linalg -from numpy.linalg import multi_dot +from numpy import ( + argsort, + array, + diag_indices, + einsum, + eye, + outer, + trace, + tril_indices, + zeros, + zeros_like, +) +from numpy.linalg import eigh, multi_dot from quemb.molbe.helper import get_eri, get_scfObj, get_veff @@ -136,7 +148,7 @@ def sd(self, lao, lmo, nocc, norb=None, return_orb_count=False): else: TA = schmidt_decomposition(lmo, nocc, self.fsites) self.C_lo_eo = TA - TA = numpy.dot(lao, TA) + TA = lao @ TA self.nao = TA.shape[1] self.TA = TA if return_orb_count: @@ -202,9 +214,9 @@ def get_nsocc(self, S, C, nocc, ncore=0): Projected density matrix. """ C_ = multi_dot((self.TA.T, S, C[:, ncore : ncore + nocc])) - P_ = numpy.dot(C_, C_.T) - nsocc_ = numpy.trace(P_) - nsocc = int(numpy.round(nsocc_)) + P_ = C_ @ C_.T + nsocc_ = trace(P_) + nsocc = int(round(nsocc_)) try: mo_coeffs = scipy.linalg.svd(C_)[0] except scipy.linalg.LinAlgError: @@ -251,12 +263,9 @@ def scf( eri = get_eri(dname, self.nao, eri_file=self.eri_file) if dm0 is None: - dm0 = ( - numpy.dot( - self._mo_coeffs[:, : self.nsocc], - self._mo_coeffs[:, : self.nsocc].conj().T, - ) - * 2.0 + dm0 = 2.0 * ( + self._mo_coeffs[:, : self.nsocc] + @ self._mo_coeffs[:, : self.nsocc].conj().T ) mf_ = get_scfObj(self.fock + heff, eri, self.nsocc, dm0=dm0) @@ -269,7 +278,7 @@ def scf( def update_heff(self, u, cout=None, only_chem=False): """Update the effective Hamiltonian for the fragment.""" - heff_ = numpy.zeros_like(self.h1) + heff_ = zeros_like(self.h1) if cout is None: cout = self.udim @@ -317,22 +326,18 @@ def update_ebe_hf( mo_coeffs = self._mo_coeffs if rdm_hf is None: - rdm_hf = numpy.dot( - mo_coeffs[:, : self.nsocc], mo_coeffs[:, : self.nsocc].conj().T - ) + rdm_hf = mo_coeffs[:, : self.nsocc] @ mo_coeffs[:, : self.nsocc].conj().T unrestricted_fac = 1.0 if unrestricted else 2.0 - e1 = unrestricted_fac * numpy.einsum( + e1 = unrestricted_fac * einsum( "ij,ij->i", self.h1[: self.nfsites], rdm_hf[: self.nfsites] ) ec = ( 0.5 * unrestricted_fac - * numpy.einsum( - "ij,ij->i", self.veff[: self.nfsites], rdm_hf[: self.nfsites] - ) + * einsum("ij,ij->i", self.veff[: self.nfsites], rdm_hf[: self.nfsites]) ) if self.TA.ndim == 3: @@ -346,14 +351,14 @@ def update_ebe_hf( else: eri = r[self.dname][()] - e2 = numpy.zeros_like(e1) + e2 = zeros_like(e1) for i in range(self.nfsites): for j in range(jmax): ij = i * (i + 1) // 2 + j if i > j else j * (j + 1) // 2 + i - Gij = (2.0 * rdm_hf[i, j] * rdm_hf - numpy.outer(rdm_hf[i], rdm_hf[j]))[ + Gij = (2.0 * rdm_hf[i, j] * rdm_hf - outer(rdm_hf[i], rdm_hf[j]))[ :jmax, :jmax ] - Gij[numpy.diag_indices(jmax)] *= 0.5 + Gij[diag_indices(jmax)] *= 0.5 Gij += Gij.T if ( unrestricted @@ -361,13 +366,11 @@ def update_ebe_hf( e2[i] += ( 0.5 * unrestricted_fac - * Gij[numpy.tril_indices(jmax)] + * Gij[tril_indices(jmax)] @ eri[spin_ind][ij] ) else: - e2[i] += ( - 0.5 * unrestricted_fac * Gij[numpy.tril_indices(jmax)] @ eri[ij] - ) + e2[i] += 0.5 * unrestricted_fac * Gij[tril_indices(jmax)] @ eri[ij] e_ = e1 + e2 + ec etmp = 0.0 @@ -433,7 +436,7 @@ def schmidt_decomposition( if mo_coeff is not None: C = mo_coeff[:, :nocc] if rdm is None: - Dhf = numpy.dot(C, C.T) + Dhf = C @ C.T if cinv is not None: Dhf = multi_dot((cinv, Dhf, cinv.conj().T)) else: @@ -443,15 +446,15 @@ def schmidt_decomposition( Tot_sites = Dhf.shape[0] # Identify environment sites (indices not in Frag_sites) - Env_sites1 = numpy.array([i for i in range(Tot_sites) if i not in Frag_sites]) - Env_sites = numpy.array([[i] for i in range(Tot_sites) if i not in Frag_sites]) - Frag_sites1 = numpy.array([[i] for i in Frag_sites]) + Env_sites1 = array([i for i in range(Tot_sites) if i not in Frag_sites]) + Env_sites = array([[i] for i in range(Tot_sites) if i not in Frag_sites]) + Frag_sites1 = array([[i] for i in Frag_sites]) # Compute the environment part of the density matrix Denv = Dhf[Env_sites, Env_sites.T] # Perform eigenvalue decomposition on the environment density matrix - Eval, Evec = numpy.linalg.eigh(Denv) + Eval, Evec = eigh(Denv) # Identify significant environment orbitals based on eigenvalue threshold Bidx = [] @@ -461,19 +464,19 @@ def schmidt_decomposition( if norb is not None: n_frag_ind = len(Frag_sites1) n_bath_ind = norb - n_frag_ind - ind_sort = numpy.argsort(numpy.abs(Eval)) + ind_sort = argsort(np.abs(Eval)) first_el = [x for x in ind_sort if x < 1.0 - thres][-1 * n_bath_ind] for i in range(len(Eval)): - if numpy.abs(Eval[i]) >= first_el: + if np.abs(Eval[i]) >= first_el: Bidx.append(i) else: for i in range(len(Eval)): - if thres < numpy.abs(Eval[i]) < 1.0 - thres: + if thres < np.abs(Eval[i]) < 1.0 - thres: Bidx.append(i) # Initialize the transformation matrix (TA) - TA = numpy.zeros([Tot_sites, len(Frag_sites) + len(Bidx)]) - TA[Frag_sites, : len(Frag_sites)] = numpy.eye(len(Frag_sites)) # Fragment part + TA = zeros([Tot_sites, len(Frag_sites) + len(Bidx)]) + TA[Frag_sites, : len(Frag_sites)] = eye(len(Frag_sites)) # Fragment part TA[Env_sites1, len(Frag_sites) :] = Evec[:, Bidx] # Environment part if return_orb_count: diff --git a/src/quemb/molbe/solver.py b/src/quemb/molbe/solver.py index 77be5855..636675c2 100644 --- a/src/quemb/molbe/solver.py +++ b/src/quemb/molbe/solver.py @@ -4,8 +4,18 @@ from abc import ABC from typing import Final -import numpy from attrs import Factory, define, field +from numpy import ( + allclose, + array, + asarray, + diag, + diag_indices, + einsum, + mean, + ndarray, + zeros_like, +) from numpy.linalg import multi_dot from pyscf import ao2mo, cc, fci, mcscf, mp from pyscf.cc.ccsd_rdm import make_rdm2 @@ -330,7 +340,7 @@ def be_func( h1_ = multi_dot((fobj._mf.mo_coeff.T, h1_, fobj._mf.mo_coeff)) eci, civec = ci_.kernel(h1_, eri, nmo, nelec) unused(eci) - civec = numpy.asarray(civec) + civec = asarray(civec) (rdm1a_, rdm1b_), (rdm2aa, rdm2ab, rdm2bb) = ci_.make_rdm12s( civec, nmo, nelec @@ -440,19 +450,19 @@ def be_func( elif solver == "FCI": rdm2s = mc.make_rdm2(civec, mc.norb, mc.nelec) if use_cumulant: - hf_dm = numpy.zeros_like(rdm1_tmp) - hf_dm[numpy.diag_indices(fobj.nsocc)] += 2.0 + hf_dm = zeros_like(rdm1_tmp) + hf_dm[diag_indices(fobj.nsocc)] += 2.0 del_rdm1 = rdm1_tmp.copy() - del_rdm1[numpy.diag_indices(fobj.nsocc)] -= 2.0 + del_rdm1[diag_indices(fobj.nsocc)] -= 2.0 nc = ( - numpy.einsum("ij,kl->ijkl", hf_dm, hf_dm) - + numpy.einsum("ij,kl->ijkl", hf_dm, del_rdm1) - + numpy.einsum("ij,kl->ijkl", del_rdm1, hf_dm) + einsum("ij,kl->ijkl", hf_dm, hf_dm) + + einsum("ij,kl->ijkl", hf_dm, del_rdm1) + + einsum("ij,kl->ijkl", del_rdm1, hf_dm) ) nc -= ( - numpy.einsum("ij,kl->iklj", hf_dm, hf_dm) - + numpy.einsum("ij,kl->iklj", hf_dm, del_rdm1) - + numpy.einsum("ij,kl->iklj", del_rdm1, hf_dm) + einsum("ij,kl->iklj", hf_dm, hf_dm) + + einsum("ij,kl->iklj", hf_dm, del_rdm1) + + einsum("ij,kl->iklj", del_rdm1, hf_dm) ) * 0.5 rdm2s -= nc fobj.rdm2__ = rdm2s.copy() @@ -647,7 +657,7 @@ def solve_error(Fobjs, Nocc, only_chem=False): err_chempot /= Fobjs[0].unitcell_nkpt err = err_chempot - Nocc - return abs(err), numpy.asarray([err]) + return abs(err), asarray([err]) # Compute edge and chemical potential errors for fobj in Fobjs: @@ -678,14 +688,14 @@ def solve_error(Fobjs, Nocc, only_chem=False): err_cen.append(Fobjs[fobj.center[cindx]]._rdm1[cens[j_], cens[k_]]) err_cen.append(Nocc) - err_edge = numpy.array(err_edge) - err_cen = numpy.array(err_cen) + err_edge = array(err_edge) + err_cen = array(err_cen) # Compute the error vector err_vec = err_edge - err_cen # Compute the norm of the error vector - norm_ = numpy.mean(err_vec * err_vec) ** 0.5 + norm_ = mean(err_vec * err_vec) ** 0.5 return norm_, err_vec @@ -809,7 +819,7 @@ def solve_ccsd( # Prepare the integrals and Fock matrix eris = mycc.ao2mo() eris.mo_energy = mo_energy - eris.fock = numpy.diag(mo_energy) + eris.fock = diag(mo_energy) # Solve the CCSD equations try: @@ -828,8 +838,8 @@ def solve_ccsd( # Compute and return the density matrices if requested if rdm_return: if not relax: - l1 = numpy.zeros_like(t1) - l2 = numpy.zeros_like(t2) + l1 = zeros_like(t1) + l2 = zeros_like(t2) rdm1a = cc.ccsd_rdm.make_rdm1(mycc, t1, t2, l1, l2) else: rdm1a = mycc.make_rdm1(with_frozen=False) @@ -910,20 +920,20 @@ def solve_block2( # Subtract off non-cumulant contribution to correlated 2RDM. if use_cumulant: - hf_dm = numpy.zeros_like(rdm1) - hf_dm[numpy.diag_indices(nocc)] += 2.0 + hf_dm = zeros_like(rdm1) + hf_dm[diag_indices(nocc)] += 2.0 del_rdm1 = rdm1.copy() - del_rdm1[numpy.diag_indices(nocc)] -= 2.0 + del_rdm1[diag_indices(nocc)] -= 2.0 nc = ( - numpy.einsum("ij,kl->ijkl", hf_dm, hf_dm) - + numpy.einsum("ij,kl->ijkl", hf_dm, del_rdm1) - + numpy.einsum("ij,kl->ijkl", del_rdm1, hf_dm) + einsum("ij,kl->ijkl", hf_dm, hf_dm) + + einsum("ij,kl->ijkl", hf_dm, del_rdm1) + + einsum("ij,kl->ijkl", del_rdm1, hf_dm) ) nc -= ( - numpy.einsum("ij,kl->iklj", hf_dm, hf_dm) - + numpy.einsum("ij,kl->iklj", hf_dm, del_rdm1) - + numpy.einsum("ij,kl->iklj", del_rdm1, hf_dm) + einsum("ij,kl->iklj", hf_dm, hf_dm) + + einsum("ij,kl->iklj", hf_dm, del_rdm1) + + einsum("ij,kl->iklj", del_rdm1, hf_dm) ) * 0.5 rdm2 -= nc @@ -983,7 +993,7 @@ def solve_uccsd( Vos = eris_inp[-1] def ao2mofn(moish): - if isinstance(moish, numpy.ndarray): + if isinstance(moish, ndarray): # Since inside '_make_eris_incore' it does not differentiate spin # for the two same-spin components, we here brute-forcely determine # what spin component we are dealing with by comparing the first @@ -992,7 +1002,7 @@ def ao2mofn(moish): moish_feature = moish[:2, :2] s = -1 for ss in [0, 1]: - if numpy.allclose(moish_feature, C[ss][:2, :2]): + if allclose(moish_feature, C[ss][:2, :2]): s = ss break if s < 0: @@ -1008,8 +1018,8 @@ def ao2mofn(moish): for s in [0, 1]: Cs_feature = C[s][:2, :2] if not ( - numpy.allclose(moish_feature[2 * s], Cs_feature) - and numpy.allclose(moish_feature[2 * s + 1], Cs_feature) + allclose(moish_feature[2 * s], Cs_feature) + and allclose(moish_feature[2 * s + 1], Cs_feature) ): raise ValueError( "Expect a list/tuple of 4 numpy arrays in the order " @@ -1019,7 +1029,7 @@ def ao2mofn(moish): return ao2mo.incore.general(Vos, moish, compact=False) except NotImplementedError: # ao2mo.incore.general is not implemented for complex numbers - return numpy.einsum( + return einsum( "ijkl,ip,jq,kr,ls->pqrs", Vos, moish[0], diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 648cea6c..83d24f10 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -15,7 +15,7 @@ from pathlib import Path import h5py -import numpy +from numpy import array, einsum, zeros_like from pyscf import ao2mo from pyscf.scf.uhf import UHF @@ -74,7 +74,7 @@ def __init__( self.hcore = mf.get_hcore() self.S = mf.get_ovlp() - self.C = [numpy.array(mf.mo_coeff[0]), numpy.array(mf.mo_coeff[1])] + self.C = [array(mf.mo_coeff[0]), array(mf.mo_coeff[1])] self.hf_dm = [mf.make_rdm1()[0], mf.make_rdm1()[1]] self.hf_veff = [mf.get_veff()[0], mf.get_veff()[1]] @@ -112,20 +112,18 @@ def __init__( self.Nocc[1] -= self.ncore self.hf_dm = [ - numpy.dot( - self.C[s][:, self.ncore : self.ncore + self.Nocc[s]], - self.C[s][:, self.ncore : self.ncore + self.Nocc[s]].T, - ) + self.C[s][:, self.ncore : self.ncore + self.Nocc[s]] + @ self.C[s][:, self.ncore : self.ncore + self.Nocc[s]].T for s in [0, 1] ] self.C_core = [self.C[s][:, : self.ncore] for s in [0, 1]] - self.P_core = [numpy.dot(self.C_core[s], self.C_core[s].T) for s in [0, 1]] + self.P_core = [self.C_core[s] @ self.C_core[s].T for s in [0, 1]] self.core_veff = 1.0 * mf.get_veff(dm=self.P_core) self.E_core = ( sum( [ - numpy.einsum( + einsum( "ji,ji->", 2 * self.hcore + self.core_veff[s], self.P_core[s], @@ -137,8 +135,8 @@ def __init__( ) # iao ignored for now - self.C_a = numpy.array(mf.mo_coeff[0]) - self.C_b = numpy.array(mf.mo_coeff[1]) + self.C_a = array(mf.mo_coeff[0]) + self.C_b = array(mf.mo_coeff[1]) del self.C self.localize( @@ -293,11 +291,11 @@ def initialize(self, eri_, compute_hf): fobj_a.cons_fock(self.hf_veff[0], self.S, self.hf_dm[0] * 2.0, eri_=eri_a) fobj_a.hf_veff = self.hf_veff[0] - fobj_a.heff = numpy.zeros_like(fobj_a.h1) + fobj_a.heff = zeros_like(fobj_a.h1) fobj_a.scf(fs=True, eri=eri_a) - fobj_a.dm0 = numpy.dot( - fobj_a._mo_coeffs[:, : fobj_a.nsocc], - fobj_a._mo_coeffs[:, : fobj_a.nsocc].conj().T, + fobj_a.dm0 = ( + fobj_a._mo_coeffs[:, : fobj_a.nsocc] + @ fobj_a._mo_coeffs[:, : fobj_a.nsocc].conj().T ) if compute_hf: @@ -315,12 +313,12 @@ def initialize(self, eri_, compute_hf): eri_b = ao2mo.restore(8, eri_b, fobj_b.nao) fobj_b.cons_fock(self.hf_veff[1], self.S, self.hf_dm[1] * 2.0, eri_=eri_b) fobj_b.hf_veff = self.hf_veff[1] - fobj_b.heff = numpy.zeros_like(fobj_b.h1) + fobj_b.heff = zeros_like(fobj_b.h1) fobj_b.scf(fs=True, eri=eri_b) - fobj_b.dm0 = numpy.dot( - fobj_b._mo_coeffs[:, : fobj_b.nsocc], - fobj_b._mo_coeffs[:, : fobj_b.nsocc].conj().T, + fobj_b.dm0 = ( + fobj_b._mo_coeffs[:, : fobj_b.nsocc] + @ fobj_b._mo_coeffs[:, : fobj_b.nsocc].conj().T ) if compute_hf: diff --git a/src/quemb/shared/external/ccsd_rdm.py b/src/quemb/shared/external/ccsd_rdm.py index bb8b5318..a1ea7286 100644 --- a/src/quemb/shared/external/ccsd_rdm.py +++ b/src/quemb/shared/external/ccsd_rdm.py @@ -3,17 +3,17 @@ # The code has been slightly modified. # -import numpy +from numpy import asarray, diag_indices, einsum, zeros, zeros_like from pyscf.cc.uccsd_rdm import make_rdm1, make_rdm2 def make_rdm1_ccsd_t1(t1): nocc, nvir = t1.shape nmo = nocc + nvir - dm = numpy.zeros((nmo, nmo), dtype=t1.dtype) + dm = zeros((nmo, nmo), dtype=t1.dtype) dm[:nocc, nocc:] = t1 dm[nocc:, :nocc] = t1.T - dm[numpy.diag_indices(nocc)] += 2.0 + dm[diag_indices(nocc)] += 2.0 return dm @@ -22,12 +22,12 @@ def make_rdm2_urlx(t1, t2, with_dm1=True): nocc, nvir = t1.shape nmo = nocc + nvir - goovv = (numpy.einsum("ia,jb->ijab", t1, t1) + t2) * 0.5 + goovv = (einsum("ia,jb->ijab", t1, t1) + t2) * 0.5 dovov = goovv.transpose(0, 2, 1, 3) * 2 - goovv.transpose(1, 2, 0, 3) - dm2 = numpy.zeros([nmo, nmo, nmo, nmo], dtype=t1.dtype) + dm2 = zeros([nmo, nmo, nmo, nmo], dtype=t1.dtype) - dovov = numpy.asarray(dovov) + dovov = asarray(dovov) dm2[:nocc, nocc:, :nocc, nocc:] = dovov dm2[:nocc, nocc:, :nocc, nocc:] += dovov.transpose(2, 3, 0, 1) dm2[nocc:, :nocc, nocc:, :nocc] = ( @@ -37,7 +37,7 @@ def make_rdm2_urlx(t1, t2, with_dm1=True): if with_dm1: dm1 = make_rdm1_ccsd_t1(t1) - dm1[numpy.diag_indices(nocc)] -= 2 + dm1[diag_indices(nocc)] -= 2 for i in range(nocc): dm2[i, i, :, :] += dm1 * 2 @@ -57,8 +57,8 @@ def make_rdm1_uccsd(ucc, relax=False): if relax: rdm1 = make_rdm1(ucc, ucc.t1, ucc.t2, ucc.l1, ucc.l2) else: - l1 = [numpy.zeros_like(ucc.t1[s]) for s in [0, 1]] - l2 = [numpy.zeros_like(ucc.t2[s]) for s in [0, 1, 2]] + l1 = [zeros_like(ucc.t1[s]) for s in [0, 1]] + l2 = [zeros_like(ucc.t2[s]) for s in [0, 1, 2]] rdm1 = make_rdm1(ucc, ucc.t1, ucc.t2, l1, l2) return rdm1 @@ -67,7 +67,7 @@ def make_rdm2_uccsd(ucc, relax=False, with_dm1=True): if relax: rdm2 = make_rdm2(ucc, ucc.t1, ucc.t2, ucc.l1, ucc.l2, with_dm1=with_dm1) else: - l1 = [numpy.zeros_like(ucc.t1[s]) for s in [0, 1]] - l2 = [numpy.zeros_like(ucc.t2[s]) for s in [0, 1, 2]] + l1 = [zeros_like(ucc.t1[s]) for s in [0, 1]] + l2 = [zeros_like(ucc.t2[s]) for s in [0, 1, 2]] rdm2 = make_rdm2(ucc, ucc.t1, ucc.t2, l1, l2, with_dm1=with_dm1) return rdm2 diff --git a/src/quemb/shared/external/lo_helper.py b/src/quemb/shared/external/lo_helper.py index 55af2927..0e8c1276 100644 --- a/src/quemb/shared/external/lo_helper.py +++ b/src/quemb/shared/external/lo_helper.py @@ -4,7 +4,8 @@ # The code has been slightly modified. # -import numpy +import numpy as np +from numpy import diag, where from numpy.linalg import eigh, matrix_power, norm @@ -21,7 +22,7 @@ def get_symm_mat_pow(A, p, check_symm=True, thresh=1.0e-8): assert norm(A - A.conj().T) < thresh e, u = eigh(A) - Ap = u @ numpy.diag(e**p) @ u.conj().T + Ap = u @ diag(e**p) @ u.conj().T return Ap @@ -60,8 +61,8 @@ def reorder_by_atom_(Clo, aoind_by_atom, S, thr=0.5): loshift = 0 for ia in range(natom): ra = aoind_by_atom[ia] - poplo_by_atom = numpy.sum(Clo_soao[ra] ** 2.0, axis=0) - loind_a = numpy.where(poplo_by_atom > thr)[0].tolist() + poplo_by_atom = np.sum(Clo_soao[ra] ** 2.0, axis=0) + loind_a = where(poplo_by_atom > thr)[0].tolist() loind_reorder += loind_a nlo_a = len(loind_a) loind_by_atom[ia] = list(range(loshift, loshift + nlo_a)) diff --git a/src/quemb/shared/external/optqn.py b/src/quemb/shared/external/optqn.py index c03b6e09..f68c1d5f 100644 --- a/src/quemb/shared/external/optqn.py +++ b/src/quemb/shared/external/optqn.py @@ -6,7 +6,7 @@ # The code has been slightly modified. # -import numpy +from numpy import array, empty, outer, zeros from numpy.linalg import inv, norm, pinv from quemb.molbe.helper import get_eri, get_scfObj @@ -121,7 +121,7 @@ def trustRegion(func, xold, fold, Binv, c=0.5): ) tdx_sd = t * dx_sd diff = dx_gn - tdx_sd - # s = (-dx_sd.T@diff + numpy.sqrt((dx_sd.T@diff)**2 - + # s = (-dx_sd.T@diff + sqrt((dx_sd.T@diff)**2 - # norm(diff)**2*(norm(dx_sd)**2-(c ** microiter)**2))) # / (norm(dx_sd))**2 # s is largest value in [0, 1] s.t. ||dx|| \le trust radius @@ -130,7 +130,7 @@ def trustRegion(func, xold, fold, Binv, c=0.5): while norm(dx) > c**microiter and s > 0: s -= 0.001 dx = tdx_sd + s * diff - if prevdx is None or not numpy.all(dx == prevdx): + if prevdx is None or not all(dx == prevdx): # Actual Reduction := f(x_k) - f(x_k + dx) fnew = func(xold + dx) ared = 0.5 * (norm(fold) ** 2 - norm(fnew) ** 2) @@ -140,7 +140,7 @@ def trustRegion(func, xold, fold, Binv, c=0.5): # r = ared/pred \le rho ratio = ared / pred microiter += 1 - if prevdx is None or not numpy.all(dx == prevdx) and settings.PRINT_LEVEL > 2: + if prevdx is None or not all(dx == prevdx) and settings.PRINT_LEVEL > 2: print(" ||δx||: ", norm(dx), flush=True) print( " Reduction Ratio (Actual / Predicted): ", @@ -180,10 +180,10 @@ def __init__(self, func, x0, f0, J0, trust=0.5, max_space=500): self.fnew = None # new jacobian? self.fold = None # old jacobian? self.max_subspace = max_space - self.dxs = numpy.empty([self.max_subspace, self.n]) - self.fs = numpy.empty([self.max_subspace, self.n]) - self.us = numpy.empty([self.max_subspace, self.n]) # u_m = B_m @ f_m - self.vs = numpy.empty([self.max_subspace, self.n]) # v_m = B_0 @ f_{m+1} + self.dxs = empty([self.max_subspace, self.n]) + self.fs = empty([self.max_subspace, self.n]) + self.us = empty([self.max_subspace, self.n]) # u_m = B_m @ f_m + self.vs = empty([self.max_subspace, self.n]) # v_m = B_0 @ f_{m+1} self.B = None self.trust = trust @@ -192,7 +192,7 @@ def next_step(self, trust_region=False): self.xnew = self.x0 self.fnew = self.func(self.xnew) if self.f0 is None else self.f0 self.fs[0] = self.fnew.copy() - self.us[0] = numpy.dot(self.B0, self.fnew) + self.us[0] = self.B0 @ self.fnew self.Binv = self.B0.copy() # Book keeping @@ -204,7 +204,7 @@ def next_step(self, trust_region=False): self.fold = self.fnew.copy() if not self.iter_ == 0: - tmp__ = numpy.outer(dx_i - self.Binv @ df_i, dx_i @ self.Binv) / ( + tmp__ = outer(dx_i - self.Binv @ df_i, dx_i @ self.Binv) / ( dx_i @ self.Binv @ df_i ) self.Binv += tmp__ @@ -221,7 +221,7 @@ def next_step(self, trust_region=False): ) # udpate vs, dxs, and fs - self.vs[self.iter_] = numpy.dot(self.B0, self.fnew) + self.vs[self.iter_] = self.B0 @ self.fnew self.dxs[self.iter_] = self.xnew - self.xold self.fs[self.iter_ + 1] = self.fnew.copy() @@ -285,12 +285,12 @@ def get_be_error_jacobian(Nfrag, Fobjs, jac_solver="HF"): M4 C2-2 E3 """ N_ = sum(Ncout) - J = numpy.zeros((N_ + 1, N_ + 1)) + J = zeros((N_ + 1, N_ + 1)) cout = 0 for findx, fobj in enumerate(Fobjs): J[cout : Ncout[findx] + cout, cout : Ncout[findx] + cout] = Jes[findx] - J[cout : Ncout[findx] + cout, N_:] = numpy.array(xes[findx]).reshape(-1, 1) + J[cout : Ncout[findx] + cout, N_:] = array(xes[findx]).reshape(-1, 1) J[N_:, cout : Ncout[findx] + cout] = ys[findx] coutc = 0 @@ -300,7 +300,7 @@ def get_be_error_jacobian(Nfrag, Fobjs, jac_solver="HF"): start_ = sum(Ncout[: fobj.center[cindx]]) end_ = start_ + Ncout[fobj.center[cindx]] J[cout + coutc_ : cout + coutc, start_:end_] += Jcs[fobj.center[cindx]] - J[cout + coutc_ : cout + coutc, N_:] += numpy.array( + J[cout + coutc_ : cout + coutc, N_:] += array( xcs[fobj.center[cindx]] ).reshape(-1, 1) coutc_ = coutc @@ -313,10 +313,7 @@ def get_be_error_jacobian(Nfrag, Fobjs, jac_solver="HF"): def get_atbe_Jblock_frag(fobj, res_func): vpots = get_vpots_frag(fobj.nao, fobj.edge_idx, fobj.fsites) eri_ = get_eri(fobj.dname, fobj.nao, eri_file=fobj.eri_file) - dm0 = ( - numpy.dot(fobj._mo_coeffs[:, : fobj.nsocc], fobj._mo_coeffs[:, : fobj.nsocc].T) - * 2.0 - ) + dm0 = 2.0 * (fobj._mo_coeffs[:, : fobj.nsocc] @ fobj._mo_coeffs[:, : fobj.nsocc].T) mf_ = get_scfObj(fobj.fock + fobj.heff, eri_, fobj.nsocc, dm0=dm0) dPs, dP_mu = res_func(mf_, vpots, eri_, fobj.nsocc) @@ -370,8 +367,8 @@ def get_atbe_Jblock_frag(fobj, res_func): # edge xe.append(dP_mu[edge[j_], edge[k_]]) cout += 1 - Je = numpy.array(Je).T - Jc = numpy.array(Jc).T + Je = array(Je).T + Jc = array(Jc).T alpha = 0.0 for fidx, _ in enumerate(fobj.fsites): @@ -403,12 +400,12 @@ def get_be_error_jacobian_selffrag(self, jac_solver="HF"): Jes, _, xes, xcs, ys, alphas, Ncout = get_atbe_Jblock_frag(self.Fobjs[0], res_func) N_ = Ncout - J = numpy.zeros((N_ + 1, N_ + 1)) + J = zeros((N_ + 1, N_ + 1)) J[:Ncout, :Ncout] = Jes - J[:Ncout, N_:] = numpy.array(xes).reshape(-1, 1) + J[:Ncout, N_:] = array(xes).reshape(-1, 1) J[N_:, :Ncout] = ys - J[:Ncout, N_:] += numpy.array([*xcs, *xcs]).reshape(-1, 1) + J[:Ncout, N_:] += array([*xcs, *xcs]).reshape(-1, 1) J[N_:, N_:] = alphas return J @@ -434,7 +431,7 @@ def mp2res_func(mf, vpots, eri, nsocc): no = nsocc dPs_an = get_dPmp2_batch_r(C, moe, eri, no, vpots, aorep=True) - dPs_an = numpy.array([dp_ * 0.5 for dp_ in dPs_an]) + dPs_an = array([dp_ * 0.5 for dp_ in dPs_an]) dP_mu = dPs_an[-1] return dPs_an[:-1], dP_mu @@ -463,13 +460,13 @@ def get_vpots_frag(nao, edge_idx, fsites): if j__ > k__: continue - tmppot = numpy.zeros((nao, nao)) + tmppot = zeros((nao, nao)) tmppot[edge_[j__], edge_[k__]] = tmppot[edge_[k__], edge_[j__]] = 1 vpots.append(tmppot) # only the centers # outer edges not included - tmppot = numpy.zeros((nao, nao)) + tmppot = zeros((nao, nao)) for fidx, fval in enumerate(fsites): if not any(fidx in sublist for sublist in edge_idx): tmppot[fidx, fidx] = -1 diff --git a/src/quemb/shared/external/uccsd_eri.py b/src/quemb/shared/external/uccsd_eri.py index 8c513cfd..41646e21 100644 --- a/src/quemb/shared/external/uccsd_eri.py +++ b/src/quemb/shared/external/uccsd_eri.py @@ -1,4 +1,4 @@ -import numpy +from numpy import einsum from pyscf.cc.uccsd import _make_eris_incore @@ -12,13 +12,12 @@ def make_eris_incore(mycc, Vss, Vos, mo_coeff=None, ao2mofn=None, frozen=False): def frank_get_veff(dm, Vss, Vos): veffss = [ - numpy.einsum("pqrs,sr->pq", Vss[s], dm[s]) - - numpy.einsum("psrq,sr->pq", Vss[s], dm[s]) + einsum("pqrs,sr->pq", Vss[s], dm[s]) - einsum("psrq,sr->pq", Vss[s], dm[s]) for s in [0, 1] ] veffos = [ - numpy.einsum("pqrs,sr->pq", Vos, dm[1]), - numpy.einsum("pqrs,qp->rs", Vos, dm[0]), + einsum("pqrs,sr->pq", Vos, dm[1]), + einsum("pqrs,qp->rs", Vos, dm[0]), ] veff = [veffss[s] + veffos[s] for s in [0, 1]] diff --git a/src/quemb/shared/external/unrestricted_utils.py b/src/quemb/shared/external/unrestricted_utils.py index a94588d1..ac6a3736 100644 --- a/src/quemb/shared/external/unrestricted_utils.py +++ b/src/quemb/shared/external/unrestricted_utils.py @@ -3,7 +3,7 @@ import ctypes import h5py -import numpy +from numpy import empty from numpy.linalg import multi_dot from pyscf import ao2mo, lib, scf @@ -139,11 +139,11 @@ def _convert_eri_gen(origsym, targetsym, eri, norb1, norb2): fn = getattr(libao2mo, "AO2MOrestore_nr%sto%s_gen" % (origsym, targetsym)) if targetsym == 1: - eri_out = numpy.empty([norb1, norb1, norb2, norb2]) + eri_out = empty([norb1, norb1, norb2, norb2]) elif targetsym == 4: npair1 = norb1 * (norb1 + 1) // 2 npair2 = norb2 * (norb2 + 1) // 2 - eri_out = numpy.empty([npair1, npair2]) + eri_out = empty([npair1, npair2]) fn( eri.ctypes.data_as(ctypes.c_void_p), diff --git a/tests/chem_dm_kBE_test.py b/tests/chem_dm_kBE_test.py index dff2ceea..9f720de1 100644 --- a/tests/chem_dm_kBE_test.py +++ b/tests/chem_dm_kBE_test.py @@ -8,7 +8,7 @@ import unittest -import numpy +from numpy import eye from pyscf.pbc import df, gto, scf from quemb.kbe import BE, fragpart @@ -29,7 +29,7 @@ def test_kc2_sto3g_be1_chempot(self) -> None: b = 1.0 c = 12.0 - lat = numpy.eye(3) + lat = eye(3) lat[0, 0] = a lat[1, 1] = b lat[2, 2] = c @@ -54,7 +54,7 @@ def test_kc4_sto3g_be2_density(self) -> None: b = 1.0 c = 12.0 - lat = numpy.eye(3) + lat = eye(3) lat[0, 0] = a lat[1, 1] = b lat[2, 2] = c