From d294de2a04a1a72ce3bd70b2ebbf30b321e73cce Mon Sep 17 00:00:00 2001 From: "puzhichen.996" Date: Thu, 16 Jan 2025 10:20:40 +0800 Subject: [PATCH] remove some comments and unused codes --- gpu4pyscf/tdscf/_lr_eig.py | 148 ++----------------------------- gpu4pyscf/tdscf/parameter_ris.py | 115 ------------------------ gpu4pyscf/tdscf/rhf.py | 11 +-- gpu4pyscf/tdscf/rks.py | 1 + 4 files changed, 14 insertions(+), 261 deletions(-) delete mode 100644 gpu4pyscf/tdscf/parameter_ris.py diff --git a/gpu4pyscf/tdscf/_lr_eig.py b/gpu4pyscf/tdscf/_lr_eig.py index 59f206ff..ae72899f 100644 --- a/gpu4pyscf/tdscf/_lr_eig.py +++ b/gpu4pyscf/tdscf/_lr_eig.py @@ -84,6 +84,7 @@ def eigh(aop, x0, precond, tol_residual=1e-5, lindep=1e-12, nroots=1, log = verbose else: log = logger.Logger(sys.stdout, verbose) + if isinstance(x0, np.ndarray) and x0.ndim == 1: x0 = x0[None,:] @@ -114,7 +115,7 @@ def eigh(aop, x0, precond, tol_residual=1e-5, lindep=1e-12, nroots=1, xs_ir = cp.array([], dtype=xt_ir.dtype) for icyc in range(max_cycle): - xt, xt_idx = _qr_gpu(xt, lindep) + xt, xt_idx = _qr(xt, lindep) # Generate at most space_inc trial vectors xt = xt[:space_inc] xt_idx = xt_idx[:space_inc] @@ -300,6 +301,7 @@ def eig(aop, x0, precond, tol_residual=1e-5, nroots=1, x0sym=None, pick=None, log = verbose else: log = logger.Logger(sys.stdout, verbose) + if isinstance(x0, np.ndarray) and x0.ndim == 1: x0 = x0[None,:] x0 = np.asarray(x0) @@ -636,7 +638,7 @@ def real_eig(aop, x0, precond, tol_residual=1e-5, nroots=1, x0sym=None, pick=Non pi[m0:m1,:m1] = pi_block if x0sym is None: - omega, x, y = TDDFT_subspace_eigen_solver_gpu( + omega, x, y = TDDFT_subspace_eigen_solver( a[:m1,:m1], b[:m1,:m1], sigma[:m1,:m1], pi[:m1,:m1], space_inc) else: # Diagonalize within eash symmetry sectors @@ -647,7 +649,7 @@ def real_eig(aop, x0, precond, tol_residual=1e-5, nroots=1, x0sym=None, pick=Non i1 = 0 for ir in set(xs_ir): idx = cp.nonzero(xs_ir[:m1] == ir)[0] - _w, _x, _y = TDDFT_subspace_eigen_solver_gpu( + _w, _x, _y = TDDFT_subspace_eigen_solver( a[idx[:,None],idx], b[idx[:,None],idx], sigma[idx[:,None],idx], pi[idx[:,None],idx], idx.size) i0, i1 = i1, i1 + idx.size @@ -719,7 +721,7 @@ def real_eig(aop, x0, precond, tol_residual=1e-5, nroots=1, x0sym=None, pick=Non X_new = XY_new[:,:A_size].T Y_new = XY_new[:,A_size:].T if x0sym is None: - V, W = VW_Gram_Schmidt_fill_holder_gpu( + V, W = VW_Gram_Schmidt_fill_holder( V_holder[:,:m1], W_holder[:,:m1], X_new, Y_new) else: xt_ir = xt_ir[r_index] @@ -728,7 +730,7 @@ def real_eig(aop, x0, precond, tol_residual=1e-5, nroots=1, x0sym=None, pick=Non W = [] for ir in set(xt_ir): idx = cp.nonzero(xt_ir == ir)[0] - _V, _W = VW_Gram_Schmidt_fill_holder_gpu( + _V, _W = VW_Gram_Schmidt_fill_holder( V_holder[:,:m1], W_holder[:,:m1], X_new[:,idx], Y_new[:,idx]) V.append(_V) W.append(_W) @@ -786,23 +788,6 @@ def _conj_dot(xi, xj): return xi[:,n:].dot(xj[:,:n].T) + xi[:,:n].dot(xj[:,n:].T) def _qr(xs, lindep=1e-14): - '''QR decomposition for a list of vectors (for linearly independent vectors only). - xs = (r.T).dot(qs) - ''' - nv = 0 - idx = [] - for i, xi in enumerate(xs): - for j in range(nv): - prod = xs[j].conj().dot(xi) - xi -= xs[j] * prod - norm = np.linalg.norm(xi) - if norm**2 > lindep: - xs[nv] = xi/norm - nv += 1 - idx.append(i) - return xs[:nv], idx - -def _qr_gpu(xs, lindep=1e-14): '''QR decomposition for a list of vectors (for linearly independent vectors only). xs = (r.T).dot(qs) ''' @@ -819,6 +804,7 @@ def _qr_gpu(xs, lindep=1e-14): idx.append(i) return xs[:nv], idx + def _symmetric_orth(xt, lindep=1e-6): xt = np.asarray(xt) if xt.dtype == np.float64: @@ -930,52 +916,6 @@ def TDDFT_subspace_eigen_solver(a, b, sigma, pi, nroots): ''' [ a b ] x - [ σ π] x Ω = 0 ''' ''' [ b a ] y [-π -σ] y = 0 ''' - d = abs(np.diag(sigma)) - d_mh = d**(-0.5) - - s_m_p = np.einsum('i,ij,j->ij', d_mh, sigma - pi, d_mh) - - '''LU = d^−1/2 (σ − π) d^−1/2''' - ''' A = LU ''' - L, U = scipy.linalg.lu(s_m_p, permute_l=True) - L_inv = np.linalg.inv(L) - U_inv = np.linalg.inv(U) - - '''U^-T d^−1/2 (a−b) d^-1/2 U^-1 = GG^T ''' - d_amb_d = np.einsum('i,ij,j->ij', d_mh, a-b, d_mh) - GGT = np.linalg.multi_dot([U_inv.T, d_amb_d, U_inv]) - - G = scipy.linalg.cholesky(GGT, lower=True) - G_inv = np.linalg.inv(G) - - ''' M = G^T L^−1 d^−1/2 (a+b) d^−1/2 L^−T G ''' - d_apb_d = np.einsum('i,ij,j->ij', d_mh, a+b, d_mh) - M = np.linalg.multi_dot([G.T, L_inv, d_apb_d, L_inv.T, G]) - - omega2, Z = np.linalg.eigh(M) - if np.any(omega2 <= 0): - idx = np.nonzero(omega2 > 0)[0] - omega2 = omega2[idx[:nroots]] - Z = Z[:,idx[:nroots]] - else: - omega2 = omega2[:nroots] - Z = Z[:,:nroots] - omega = omega2**0.5 - - ''' It requires Z^T Z = 1/Ω ''' - ''' x+y = d^−1/2 L^−T GZ Ω^-0.5 ''' - ''' x−y = d^−1/2 U^−1 G^−T Z Ω^0.5 ''' - x_p_y = np.einsum('i,ik,k->ik', d_mh, L_inv.T.dot(G.dot(Z)), omega**-0.5) - x_m_y = np.einsum('i,ik,k->ik', d_mh, U_inv.dot(G_inv.T.dot(Z)), omega**0.5) - - x = (x_p_y + x_m_y)/2 - y = x_p_y - x - return omega, x, y - -def TDDFT_subspace_eigen_solver_gpu(a, b, sigma, pi, nroots): - ''' [ a b ] x - [ σ π] x Ω = 0 ''' - ''' [ b a ] y [-π -σ] y = 0 ''' - d = abs(cp.diag(sigma)) d_mh = d**(-0.5) @@ -996,7 +936,6 @@ def TDDFT_subspace_eigen_solver_gpu(a, b, sigma, pi, nroots): ''' M = G^T L^−1 d^−1/2 (a+b) d^−1/2 L^−T G ''' d_apb_d = cp.einsum('i,ij,j->ij', d_mh, a+b, d_mh) - # M = cp.linalg.multi_dot([G.T, L_inv, d_apb_d, L_inv.T, G]) M = cp.dot(G.T, cp.dot(L_inv, cp.dot(d_apb_d, cp.dot(L_inv.T, G)))) omega2, Z = cp.linalg.eigh(M) @@ -1019,77 +958,8 @@ def TDDFT_subspace_eigen_solver_gpu(a, b, sigma, pi, nroots): y = x_p_y - x return omega, x, y -def VW_Gram_Schmidt_fill_holder(V_holder, W_holder, X_new, Y_new, lindep=1e-6): - ''' - QR orthogonalization for (X_new, Y_new) basis vectors, then apply symmetric - orthogonalization for {[X, Y]}, and its dual basis vectors {[Y, X]} - ''' - _x = V_holder.T.dot(X_new) - _x += W_holder.T.dot(Y_new) - _y = V_holder.T.dot(Y_new) - _y += W_holder.T.dot(X_new) - X_new -= V_holder.dot(_x) - X_new -= W_holder.dot(_y) - Y_new -= W_holder.dot(_x) - Y_new -= V_holder.dot(_y) - x0_size = X_new.shape[0] - - s11 = X_new.T.dot(X_new) - s11 += Y_new.T.dot(Y_new) - # s21 is symmetric - s21 = X_new.T.dot(Y_new) - s21 += Y_new.T.dot(X_new) - e, c = np.linalg.eigh(s11) - mask = e > lindep**2 - e = e[mask] - if e.size == 0: - return (np.zeros([0, x0_size], dtype=X_new.dtype), - np.zeros([0, x0_size], dtype=Y_new.dtype)) - c = c[:,mask] * e**-.5 - - csc = c.T.dot(s21).dot(c) - n = csc.shape[0] - for i in range(n): - w, u = np.linalg.eigh(csc[i:,i:]) - mask = 1 - abs(w) > lindep - if np.any(mask): - c = c[:,i:] - break - else: - return (np.zeros([0, x0_size], dtype=X_new.dtype), - np.zeros([0, x0_size], dtype=Y_new.dtype)) - w = w[mask] - u = u[:,mask] - c_orth = c.dot(u) - - if e[0] < 1e-6 or 1-abs(w[0]) < 1e-3: - # Rerun the orthogonalization to reduce numerical errors - e, c = np.linalg.eigh(c_orth.T.dot(s11).dot(c_orth)) - c *= e**-.5 - c_orth = c_orth.dot(c) - csc = c_orth.T.dot(s21).dot(c_orth) - w, u = np.linalg.eigh(csc) - c_orth = c_orth.dot(u) - - # Symmetric diagonalize - # [1 w] => c = [a b] - # [w 1] [b a] - # where - # a = ((1+w)**-.5 + (1-w)**-.5)/2 - # b = ((1+w)**-.5 - (1-w)**-.5)/2 - a1 = (1 + w)**-.5 - a2 = (1 - w)**-.5 - a = (a1 + a2) / 2 - b = (a1 - a2) / 2 - - x_orth = X_new.dot(c_orth * a) - x_orth += Y_new.dot(c_orth * b) - y_orth = Y_new.dot(c_orth * a) - y_orth += X_new.dot(c_orth * b) - return x_orth.T, y_orth.T - -def VW_Gram_Schmidt_fill_holder_gpu(V_holder, W_holder, X_new, Y_new, lindep=1e-6): +def VW_Gram_Schmidt_fill_holder(V_holder, W_holder, X_new, Y_new, lindep=1e-6): ''' QR orthogonalization for (X_new, Y_new) basis vectors, then apply symmetric orthogonalization for {[X, Y]}, and its dual basis vectors {[Y, X]} diff --git a/gpu4pyscf/tdscf/parameter_ris.py b/gpu4pyscf/tdscf/parameter_ris.py deleted file mode 100644 index 6b9ba05a..00000000 --- a/gpu4pyscf/tdscf/parameter_ris.py +++ /dev/null @@ -1,115 +0,0 @@ -# Description: This file contains all the parameters used in the code -# copied from https://github.com/John-zzh/pyscf_TDDFT_ris/blob/main/pyscf_TDDFT_ris/parameter.py - -''' -GB Radii -Ghosh, Dulal C and coworkers -The wave mechanical evaluation of the absolute radii of atoms. -Journal of Molecular Structure: THEOCHEM 865, no. 1-3 (2008): 60-67. -''' - -elements_106 = ['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', -'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', -'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', -'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', -'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', -'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', -'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', -'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr'] - -radii = [0.5292, 0.3113, 1.6283, 1.0855, 0.8141, 0.6513, 0.5428, 0.4652, 0.4071, 0.3618, -2.165, 1.6711, 1.3608, 1.1477, 0.9922, 0.8739, 0.7808, 0.7056, 3.293, 2.5419, -2.4149, 2.2998, 2.1953, 2.1, 2.0124, 1.9319, 1.8575, 1.7888, 1.725, 1.6654, -1.4489, 1.2823, 1.145, 1.0424, 0.9532, 0.8782, 3.8487, 2.9709, 2.8224, 2.688, -2.5658, 2.4543, 2.352, 2.2579, 2.1711, 2.0907, 2.016, 1.9465, 1.6934, 1.4986, -1.344, 1.2183, 1.1141, 1.0263, 4.2433, 3.2753, 2.6673, 2.2494, 1.9447, 1.7129, -1.5303, 1.383, 1.2615, 1.1596, 1.073, 0.9984, 0.9335, 0.8765, 0.8261, 0.7812, -0.7409, 0.7056, 0.6716, 0.6416, 0.6141, 0.589, 0.5657, 0.5443, 0.5244, 0.506, -1.867, 1.6523, 1.4818, 1.3431, 1.2283, 1.1315, 4.4479, 3.4332, 3.2615, 3.1061, -2.2756, 1.9767, 1.7473, 1.4496, 1.2915, 1.296, 1.1247, 1.0465, 0.9785, 0.9188, -0.8659, 0.8188, 0.8086] -exp = [1/(i*1.8897259885789)**2 for i in radii] - -ris_exp = dict(zip(elements_106,exp)) - -''' -range-separated hybrid functionals, (omega, alpha, beta) -''' -rsh_func = {} -rsh_func['wb97'] = (0.4, 0, 1.0) -rsh_func['wb97x'] = (0.3, 0.157706, 0.842294) # wb97 family, a+b=100% Long-range HF exchange -rsh_func['wb97x-d'] = (0.2, 0.22, 0.78) -rsh_func['wb97x-d3'] = (0.25, 0.195728, 0.804272) -rsh_func['wb97x-v'] = (0.30, 0.167, 0.833) -rsh_func['wb97x-d3bj'] = (0.30, 0.167, 0.833) -rsh_func['cam-b3lyp'] = (0.33, 0.19, 0.46) # a+b=65% Long-range HF exchange -rsh_func['lc-blyp'] = (0.33, 0, 1.0) -rsh_func['lc-PBE'] = (0.47, 0, 1.0) - -''' -pure or hybrid functionals, hybrid component a_x -''' -hbd_func = {} -hbd_func['pbe'] = 0 -hbd_func['pbe,pbe'] = 0 -hbd_func['tpss'] = 0 -hbd_func['tpssh'] = 0.1 -hbd_func['b3lyp'] = 0.2 -hbd_func['pbe0'] = 0.25 -hbd_func['bhh-lyp'] = 0.5 -hbd_func['m05-2x'] = 0.56 -hbd_func['m06'] = 0.27 -hbd_func['m06-2x'] = 0.54 -hbd_func[None] = 1 - - - - -''' for sTDA - a dictionary of chemical hardness, by mappig two lists: - list of elements 1-94 - list of hardness for elements 1-94, floats, in Hartree -''' -elements = ['H' , 'He', 'Li', 'Be', 'B' , 'C' , 'N' , 'O' , 'F' , 'Ne', -'Na', 'Mg', 'Al', 'Si', 'P' , 'S' , 'Cl', 'Ar', 'K' , 'Ca','Sc', 'Ti', -'V' , 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn','Ga', 'Ge', 'As', 'Se', -'Br', 'Kr', 'Rb', 'Sr', 'Y' , 'Zr','Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', -'Ag', 'Cd', 'In', 'Sn','Sb', 'Te', 'I' , 'Xe', 'Cs', 'Ba', 'La', 'Ce', -'Pr', 'Nd','Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', -'Lu', 'Hf', 'Ta', 'W' , 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg','Tl', 'Pb', -'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U' , 'Np', 'Pu'] -hardness = [0.47259288,0.92203391,0.17452888,0.25700733,0.33949086, -0.42195412,0.50438193,0.58691863,0.66931351,0.75191607,0.17964105, -0.22157276,0.26348578,0.30539645,0.34734014,0.38924725,0.43115670, -0.47308269,0.17105469,0.20276244,0.21007322,0.21739647,0.22471039, -0.23201501,0.23933969,0.24665638,0.25398255,0.26128863,0.26859476, -0.27592565,0.30762999,0.33931580,0.37235985,0.40273549,0.43445776, -0.46611708,0.15585079,0.18649324,0.19356210,0.20063311,0.20770522, -0.21477254,0.22184614,0.22891872,0.23598621,0.24305612,0.25013018, -0.25719937,0.28784780,0.31848673,0.34912431,0.37976593,0.41040808, -0.44105777,0.05019332,0.06762570,0.08504445,0.10247736,0.11991105, -0.13732772,0.15476297,0.17218265,0.18961288,0.20704760,0.22446752, -0.24189645,0.25932503,0.27676094,0.29418231,0.31159587,0.32902274, -0.34592298,0.36388048,0.38130586,0.39877476,0.41614298,0.43364510, -0.45104014,0.46848986,0.48584550,0.12526730,0.14268677,0.16011615, -0.17755889,0.19497557,0.21240778,0.07263525,0.09422158,0.09920295, -0.10418621,0.14235633,0.16394294,0.18551941,0.22370139] -HARDNESS = dict(zip(elements,hardness)) - -def gen_sTDA_alpha_beta_ax(a_x): - - ''' NA is for Hartree-Fork - - RSH functionals have specific a_x, beta, alpha values; - hybride fucntionals have fixed alpha12 and beta12 values, - with different a_x values, by which create beta, alpha - ''' - beta1 = 0.2 - beta2 = 1.83 - alpha1 = 1.42 - alpha2 = 0.48 - - beta = beta1 + beta2 * a_x - alpha = alpha1 + alpha2 * a_x - - return alpha, beta \ No newline at end of file diff --git a/gpu4pyscf/tdscf/rhf.py b/gpu4pyscf/tdscf/rhf.py index b05836c4..73f72040 100644 --- a/gpu4pyscf/tdscf/rhf.py +++ b/gpu4pyscf/tdscf/rhf.py @@ -15,12 +15,10 @@ import numpy as np import cupy as cp -import time -from pyscf import gto from pyscf import lib from pyscf.tdscf import rhf as tdhf_cpu from gpu4pyscf.tdscf._lr_eig import eigh as lr_eigh, real_eig -from gpu4pyscf import scf, dft +from gpu4pyscf import scf from gpu4pyscf.lib.cupy_helper import contract, tag_array from gpu4pyscf.lib import utils from gpu4pyscf.lib import logger @@ -52,7 +50,7 @@ def gen_tda_operation(mf, fock_ao=None, singlet=True, wfnsym=None): orbo2 = orbo * 2. # *2 for double occupancy e_ia = hdiag = mo_energy[viridx] - mo_energy[occidx,None] - hdiag = hdiag.ravel()#.get() + hdiag = hdiag.ravel() vresp = mf.gen_response(singlet=singlet, hermi=0) nocc, nvir = e_ia.shape @@ -65,7 +63,7 @@ def vind(zs): v1mo = contract('xpq,qo->xpo', v1ao, orbo) v1mo = contract('xpo,pv->xov', v1mo, orbv.conj()) v1mo += zs * e_ia - return v1mo.reshape(v1mo.shape[0],-1)#.get() + return v1mo.reshape(v1mo.shape[0],-1) return vind, hdiag @@ -188,7 +186,6 @@ def gen_vind(self, mf=None): '''Generate function to compute Ax''' if mf is None: mf = self._scf - return gen_tda_operation(mf, singlet=self.singlet) def init_guess(self, mf=None, nstates=None, wfnsym=None, return_symmetry=False): @@ -318,7 +315,6 @@ class TDHF(TDBase): def gen_vind(self, mf=None): if mf is None: mf = self._scf - return gen_tdhf_operation(mf, singlet=self.singlet) def init_guess(self, mf=None, nstates=None, wfnsym=None, return_symmetry=False): @@ -352,6 +348,7 @@ def kernel(self, x0=None, nstates=None): x0sym = None if x0 is None: x0 = self.init_guess() + self.converged, self.e, x1 = real_eig( vind, x0, precond, tol_residual=self.conv_tol, lindep=self.lindep, nroots=nstates, x0sym=x0sym, pick=pickeig, max_cycle=self.max_cycle, diff --git a/gpu4pyscf/tdscf/rks.py b/gpu4pyscf/tdscf/rks.py index bcaa7f37..948c7d50 100644 --- a/gpu4pyscf/tdscf/rks.py +++ b/gpu4pyscf/tdscf/rks.py @@ -100,6 +100,7 @@ def pickeig(w, v, nroots, envs): x0sym = None if x0 is None: x0 = self.init_guess() + self.converged, w2, x1 = lr_eigh( vind, x0, precond, tol_residual=self.conv_tol, lindep=self.lindep, nroots=nstates, x0sym=x0sym, pick=pickeig, max_cycle=self.max_cycle,