diff --git a/pyscf/adc/test/test_radc/test_N2_radc_ea.py b/pyscf/adc/test/test_radc/test_N2_radc_ea.py index 518c339f57..ef2441d38e 100644 --- a/pyscf/adc/test/test_radc/test_N2_radc_ea.py +++ b/pyscf/adc/test/test_radc/test_N2_radc_ea.py @@ -110,7 +110,7 @@ def test_ea_adc2x(self): self.assertAlmostEqual(p[3],1.9689940350592559, 6) dm1_exc = myadcea.make_rdm1() - self.assertAlmostEqual(rdms_test(dm1_exc[0]), 83.53454941982915, 6) + self.assertAlmostEqual(rdms_test(dm1_exc[0]), 83.53454941982915, 5) self.assertAlmostEqual(rdms_test(dm1_exc[1]), 66.8646378314392, 6) self.assertAlmostEqual(rdms_test(dm1_exc[2]), 66.29435797091572, 6) self.assertAlmostEqual(rdms_test(dm1_exc[3]), 66.29435797091572, 6) diff --git a/pyscf/grad/test/test_tdrks_grad.py b/pyscf/grad/test/test_tdrks_grad.py index b7e11cbb5b..f38c44ada3 100644 --- a/pyscf/grad/test/test_tdrks_grad.py +++ b/pyscf/grad/test/test_tdrks_grad.py @@ -129,7 +129,7 @@ def test_tda_singlet_mgga(self): self.assertAlmostEqual(abs((e1[2]-e2[2])/.002 - g1[0,2]).max(), 0, 3) def test_tddft_lda(self): - td = tdscf.TDDFT(mf_lda).run(nstates=nstates) + td = tdscf.TDDFT(mf_lda).run(nstates=nstates, conv_tol=1e-8) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -1.31315477e-01, 6) diff --git a/pyscf/grad/test/test_tduhf_grad.py b/pyscf/grad/test/test_tduhf_grad.py index cce0ed8c3a..7866b20b56 100644 --- a/pyscf/grad/test/test_tduhf_grad.py +++ b/pyscf/grad/test/test_tduhf_grad.py @@ -58,7 +58,7 @@ def test_tda(self): self.assertAlmostEqual((e1[2]-e2[2])/.002, g1[0,2], 4) def test_tdhf(self): - td = tdscf.TDDFT(mf).run(nstates=nstates) + td = tdscf.TDDFT(mf).run(nstates=nstates, conv_tol=1e-6) tdg = td.nuc_grad_method() g1 = tdg.kernel(td.xy[2]) self.assertAlmostEqual(g1[0,2], -0.78969714300299776, 5) diff --git a/pyscf/grad/test/test_tduks_grad.py b/pyscf/grad/test/test_tduks_grad.py index b83e8f6a9b..90b8037e1b 100644 --- a/pyscf/grad/test/test_tduks_grad.py +++ b/pyscf/grad/test/test_tduks_grad.py @@ -72,7 +72,7 @@ def test_tda_b88(self): self.assertAlmostEqual((e1[2]-e2[2])/.002, g1[0,2], 4) def test_tddft_lda(self): - td = tdscf.TDDFT(mf_lda).run(nstates=nstates) + td = tdscf.TDDFT(mf_lda).run(nstates=nstates, conv_tol=1e-7) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -0.800487816830773, 6) @@ -104,7 +104,7 @@ def test_tddft_b3lyp(self): mf = dft.UKS(mol).set(conv_tol=1e-12) mf.xc = '.2*HF + .8*b88, vwn' mf.scf() - td = tdscf.TDDFT(mf).run(nstates=nstates) + td = tdscf.TDDFT(mf).run(nstates=nstates, conv_tol=1e-6) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -0.80446691153291727, 6) diff --git a/pyscf/lib/linalg_helper.py b/pyscf/lib/linalg_helper.py index 59a78b5c68..161a09463b 100644 --- a/pyscf/lib/linalg_helper.py +++ b/pyscf/lib/linalg_helper.py @@ -402,7 +402,6 @@ def davidson1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=12, v = None conv = numpy.zeros(nroots, dtype=bool) emin = None - level_shift = 0 for icyc in range(max_cycle): if fresh_start: @@ -475,9 +474,9 @@ def davidson1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=12, else: e = w[:nroots] v = v[:,:nroots] - conv = numpy.zeros(nroots, dtype=bool) - elast, conv_last = _sort_elast(elast, conv_last, vlast, v, - fresh_start, log) + conv = numpy.zeros(e.size, dtype=bool) + if not fresh_start: + elast, conv_last = _sort_elast(elast, conv_last, vlast, v, log) if elast is None: de = e @@ -495,16 +494,15 @@ def davidson1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=12, else: ax0 = _gen_x0(v, ax) - dx_norm = numpy.zeros(nroots) + dx_norm = numpy.zeros(e.size) xt = [None] * nroots for k, ek in enumerate(e): - if not conv[k]: - xt[k] = ax0[k] - ek * x0[k] - dx_norm[k] = numpy.sqrt(dot(xt[k].conj(), xt[k]).real) - conv[k] = abs(de[k]) < tol and dx_norm[k] < toloose - if conv[k] and not conv_last[k]: - log.debug('root %d converged |r|= %4.3g e= %s max|de|= %4.3g', - k, dx_norm[k], ek, de[k]) + xt[k] = ax0[k] - ek * x0[k] + dx_norm[k] = numpy.sqrt(dot(xt[k].conj(), xt[k]).real) + conv[k] = abs(de[k]) < tol and dx_norm[k] < toloose + if conv[k] and not conv_last[k]: + log.debug('root %d converged |r|= %4.3g e= %s max|de|= %4.3g', + k, dx_norm[k], ek, de[k]) ax0 = None max_dx_norm = max(dx_norm) ide = numpy.argmax(abs(de)) @@ -528,7 +526,7 @@ def davidson1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=12, # remove subspace linear dependency for k, ek in enumerate(e): if (not conv[k]) and dx_norm[k]**2 > lindep: - xt[k] = precond(xt[k], e[0]-level_shift, x0[k]) + xt[k] = precond(xt[k], e[0], x0[k]) xt[k] *= dot(xt[k].conj(), xt[k]).real ** -.5 elif not conv[k]: # Remove linearly dependent vector @@ -537,21 +535,13 @@ def davidson1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=12, else: xt[k] = None - xt, ill_precond = _project_xt_(xt, xs, e, lindep, dot, precond) - if ill_precond: - # Manually adjust the precond because precond function may not be - # able to generate linearly dependent basis vectors. e.g. issue 1362 - log.warn('Matrix may be already a diagonal matrix. ' - 'level_shift is applied to precond') - level_shift = 0.1 - - xt, norm_min = _normalize_xt_(xt, lindep, dot) + xt, norm_min = _normalize_xt_(xt, xs, lindep, dot) log.debug('davidson %d %d |r|= %4.3g e= %s max|de|= %4.3g lindep= %4.3g', icyc, space, max_dx_norm, e, de[ide], norm_min) if len(xt) == 0: log.debug('Linear dependency in trial subspace. |r| for each state %s', dx_norm) - conv[dx_norm < toloose] = True + conv = dx_norm < toloose break max_dx_last = max_dx_norm @@ -570,17 +560,16 @@ def davidson1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=12, # can be generated. # 2. The initial guess sits in the subspace which is smaller than the # required number of roots. - msg = 'Not enough eigenvectors (len(x0)=%d, nroots=%d)' % (len(x0), nroots) - warnings.warn(msg) + log.warn(f'Not enough eigenvectors (len(x0)={len(x0)}, nroots={nroots})') return numpy.asarray(conv), e, x0 -def make_diag_precond(diag, level_shift=0): +def make_diag_precond(diag, level_shift=1e-3): '''Generate the preconditioner function with the diagonal function.''' # For diagonal matrix A, precond (Ax-x*e)/(diag(A)-e) is not able to - # generate linearly independent basis. Use level_shift to break the - # correlation between Ax-x*e and diag(A)-e. + # generate linearly independent basis (see issue 1362). Use level_shift to + # break the correlation between Ax-x*e and diag(A)-e. def precond(dx, e, *args): diagd = diag - (e - level_shift) diagd[abs(diagd)<1e-8] = 1e-8 @@ -787,7 +776,6 @@ def davidson_nosym1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=20, v = None conv = numpy.zeros(nroots, dtype=bool) emin = None - level_shift = 0 for icyc in range(max_cycle): if fresh_start: @@ -822,10 +810,7 @@ def davidson_nosym1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=20, head, space = space, space+rnow if dtype is None: - try: - dtype = numpy.result_type(axt[0], xt[0]) - except IndexError: - dtype = numpy.result_type(ax[0].dtype, xs[0].dtype) + dtype = numpy.result_type(*axt, *xt) if heff is None: # Lazy initialize heff to determine the dtype heff = numpy.empty((max_space+nroots,max_space+nroots), dtype=dtype) else: @@ -849,9 +834,9 @@ def davidson_nosym1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=20, else: e = w[:nroots] v = v[:,:nroots] - conv = numpy.zeros(nroots, dtype=bool) - elast, conv_last = _sort_elast(elast, conv_last, vlast, v, - fresh_start, log) + conv = numpy.zeros(e.size, dtype=bool) + if not fresh_start: + elast, conv_last = _sort_elast(elast, conv_last, vlast, v, log) if elast is None: de = e @@ -869,7 +854,7 @@ def davidson_nosym1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=20, else: ax0 = _gen_x0(v, ax) - dx_norm = numpy.zeros(nroots) + dx_norm = numpy.zeros(e.size) xt = [None] * nroots for k, ek in enumerate(e): if not conv[k]: @@ -879,7 +864,6 @@ def davidson_nosym1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=20, if conv[k] and not conv_last[k]: log.debug('root %d converged |r|= %4.3g e= %s max|de|= %4.3g', k, dx_norm[k], ek, de[k]) - conv[k] = True ax0 = None max_dx_norm = max(dx_norm) ide = numpy.argmax(abs(de)) @@ -903,7 +887,7 @@ def davidson_nosym1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=20, # remove subspace linear dependency for k, ek in enumerate(e): if (not conv[k]) and dx_norm[k]**2 > lindep: - xt[k] = precond(xt[k], e[0]-level_shift, x0[k]) + xt[k] = precond(xt[k], e[0], x0[k]) xt[k] *= dot(xt[k].conj(), xt[k]).real ** -.5 elif not conv[k]: # Remove linearly dependent vector @@ -912,21 +896,13 @@ def davidson_nosym1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=20, else: xt[k] = None - xt, ill_precond = _project_xt_(xt, xs, e, lindep, dot, precond) - if ill_precond: - # Manually adjust the precond because precond function may not be - # able to generate linearly dependent basis vectors. e.g. issue 1362 - log.warn('Matrix may be already a diagonal matrix. ' - 'level_shift is applied to precond') - level_shift = 0.1 - - xt, norm_min = _normalize_xt_(xt, lindep, dot) + xt, norm_min = _normalize_xt_(xt, xs, lindep, dot) log.debug('davidson %d %d |r|= %4.3g e= %s max|de|= %4.3g lindep= %4.3g', icyc, space, max_dx_norm, e, de[ide], norm_min) if len(xt) == 0: log.debug('Linear dependency in trial subspace. |r| for each state %s', dx_norm) - conv = [conv[k] or (norm < toloose) for k,norm in enumerate(dx_norm)] + conv = dx_norm < toloose break max_dx_last = max_dx_norm @@ -1103,7 +1079,6 @@ def dgeev1(abop, x0, precond, type=1, tol=1e-12, max_cycle=50, max_space=12, seff = numpy.empty((max_space,max_space), dtype=x0[0].dtype) fresh_start = True conv = False - level_shift = 0 for icyc in range(max_cycle): if fresh_start: @@ -1193,7 +1168,7 @@ def dgeev1(abop, x0, precond, type=1, tol=1e-12, max_cycle=50, max_space=12, conv = True break - dx_norm = numpy.zeros(nroots) + dx_norm = numpy.zeros(e.size) xt = [None] * nroots for k, ek in enumerate(e): if type == 1: @@ -1212,27 +1187,19 @@ def dgeev1(abop, x0, precond, type=1, tol=1e-12, max_cycle=50, max_space=12, # remove subspace linear dependency for k, ek in enumerate(e): if dx_norm[k]**2 > lindep: - xt[k] = precond(xt[k], e[0]-level_shift, x0[k]) + xt[k] = precond(xt[k], e[0], x0[k]) xt[k] *= dot(xt[k].conj(), xt[k]).real ** -.5 else: log.debug1('Drop eigenvector %d, norm=%4.3g', k, dx_norm[k]) xt[k] = None - xt, ill_precond = _project_xt_(xt, xs, e, lindep, dot, precond) - if ill_precond: - # Manually adjust the precond because precond function may not be - # able to generate linearly dependent basis vectors. e.g. issue 1362 - log.warn('Matrix may be already a diagonal matrix. ' - 'level_shift is applied to precond') - level_shift = 0.1 - - xt, norm_min = _normalize_xt_(xt, lindep, dot) + xt, norm_min = _normalize_xt_(xt, xs, lindep, dot) log.debug('davidson %d %d |r|= %4.3g e= %s max|de|= %4.3g lindep= %4.3g', icyc, space, max(dx_norm), e, de[ide], norm_min) if len(xt) == 0: log.debug('Linear dependency in trial subspace. |r| for each state %s', dx_norm) - conv = all(norm < toloose for norm in dx_norm) + conv = all(dx_norm < toloose) break fresh_start = fresh_start or (space+len(xt) > max_space) @@ -1495,58 +1462,47 @@ def _sort_by_similarity(w, v, nroots, conv, vlast, emin=None, heff=None): c = v[:,sorted_idx] return e, c -def _sort_elast(elast, conv_last, vlast, v, fresh_start, log): +def _sort_elast(elast, conv_last, vlast, v, log): ''' Eigenstates may be flipped during the Davidson iterations. Reorder the eigenvalues of last iteration to make them comparable to the eigenvalues of the current iterations. ''' - if fresh_start: - return elast, conv_last head, nroots = vlast.shape ovlp = abs(numpy.dot(v[:head].conj().T, vlast)) - idx = numpy.argmax(ovlp, axis=1) + mapping = numpy.argmax(ovlp, axis=1) + found = numpy.any(ovlp > .5, axis=1) if log.verbose >= logger.DEBUG: - ordering_diff = (idx != numpy.arange(len(idx))) - if numpy.any(ordering_diff): + ordering_diff = (mapping != numpy.arange(len(mapping))) + if any(ordering_diff & found): log.debug('Old state -> New state') for i in numpy.where(ordering_diff)[0]: - log.debug(' %3d -> %3d ', idx[i], i) + log.debug(' %3d -> %3d ', mapping[i], i) - return elast[idx], conv_last[idx] + conv = conv_last[mapping] + e = elast[mapping] + conv[~found] = False + e[~found] = 0. + return e, conv -def _project_xt_(xt, xs, e, threshold, dot, precond): +def _normalize_xt_(xt, xs, threshold, dot): '''Projects out existing basis vectors xs. Also checks whether the precond function is ill-conditioned''' - ill_precond = False - for i, xsi in enumerate(xs): + xt = [x for x in xt if x is not None] + for xsi in xs: xsi = numpy.asarray(xsi) - for k, xi in enumerate(xt): - if xi is None: - continue - ovlp = dot(xsi.conj(), xi) - # xs[i] == xt[k] - if abs(1 - ovlp)**2 < threshold: - ill_precond = True - # rebuild xt[k] to remove correlation between xt[k] and xs[i] - xi[:] = precond(xi, e[k], xi) - ovlp = dot(xsi.conj(), xi) - xi -= xsi * ovlp - xsi = None - return xt, ill_precond - -def _normalize_xt_(xt, threshold, dot): + for xi in xt: + xi -= xsi * dot(xsi.conj(), xi) + norm_min = 1 out = [] - for i, xi in enumerate(xt): - if xi is None: - continue - norm = dot(xi.conj(), xi).real ** .5 + for xi in xt: + norm = dot(xi.conj(), xi).real**.5 if norm**2 > threshold: - xt[i] *= 1/norm + xi *= 1/norm + out.append(xi) norm_min = min(norm_min, norm) - out.append(xt[i]) return out, norm_min @@ -1592,5 +1548,3 @@ def __len__(self): def pop(self, index): key = self.index.pop(index) del (self.scr_h5[str(key)]) - -del (SAFE_EIGH_LINDEP, DAVIDSON_LINDEP, DSOLVE_LINDEP, MAX_MEMORY) diff --git a/pyscf/pbc/df/df_jk.py b/pyscf/pbc/df/df_jk.py index 68bf20a010..57aad2949b 100644 --- a/pyscf/pbc/df/df_jk.py +++ b/pyscf/pbc/df/df_jk.py @@ -193,7 +193,8 @@ def get_j_kpts_kshift(mydf, dm_kpts, kshift, hermi=0, kpts=numpy.zeros((1,3)), k kpts_band, input_band = _format_kpts_band(kpts_band, kpts), kpts_band nband = len(kpts_band) - j_real = gamma_point(kpts_band) and not numpy.iscomplexobj(dms) + j_real = (gamma_point(kpts_band) and gamma_point(kpts[kshift]) and + not numpy.iscomplexobj(dms)) kconserv = get_kconserv_ria(mydf.cell, kpts)[kshift] diff --git a/pyscf/pbc/df/fft_jk.py b/pyscf/pbc/df/fft_jk.py index 7785d30f4a..18e52293a6 100644 --- a/pyscf/pbc/df/fft_jk.py +++ b/pyscf/pbc/df/fft_jk.py @@ -54,10 +54,10 @@ def get_j_kpts(mydf, dm_kpts, hermi=1, kpts=np.zeros((1,3)), kpts_band=None): assert cell.dimension > 1 ni = mydf._numint - make_rho, nset, nao = ni._gen_rho_evaluator(cell, dm_kpts, hermi) dm_kpts = lib.asarray(dm_kpts, order='C') dms = _format_dms(dm_kpts, kpts) nset, nkpts, nao = dms.shape[:3] + make_rho, nset, nao = ni._gen_rho_evaluator(cell, dms, hermi) coulG = tools.get_coulG(cell, mesh=mesh) ngrids = len(coulG) @@ -120,10 +120,10 @@ def get_j_e1_kpts(mydf, dm_kpts, kpts=np.zeros((1,3)), kpts_band=None): assert cell.dimension > 1 ni = mydf._numint - make_rho, nset, nao = ni._gen_rho_evaluator(cell, dm_kpts, hermi=1) dm_kpts = lib.asarray(dm_kpts, order='C') dms = _format_dms(dm_kpts, kpts) nset, nkpts, nao = dms.shape[:3] + make_rho, nset, nao = ni._gen_rho_evaluator(cell, dms, hermi=1) coulG = tools.get_coulG(cell, mesh=mesh) ngrids = len(coulG) diff --git a/pyscf/pbc/dft/numint.py b/pyscf/pbc/dft/numint.py index 7332b45240..e77b57e088 100644 --- a/pyscf/pbc/dft/numint.py +++ b/pyscf/pbc/dft/numint.py @@ -276,7 +276,8 @@ def eval_rho2(cell, ao, mo_coeff, mo_occ, non0tab=None, xctype='LDA', rho[4] -= rho5 * 2 rho[tau_idx] -= rho5 * .5 else: - rho = numint.eval_rho2(cell, ao, mo_coeff, mo_occ, non0tab, xctype, verbose) + rho = numint.eval_rho2(cell, ao, mo_coeff, mo_occ, non0tab, xctype, + with_lapl, verbose) return rho @@ -477,7 +478,7 @@ def nr_uks(ni, cell, grids, xc_code, dms, spin=1, relativity=0, hermi=1, vmata[i] += ni._vxc_mat(cell, ao_k1, wv[0], mask, xctype, shls_slice, ao_loc, v_hermi) vmatb[i] += ni._vxc_mat(cell, ao_k1, wv[1], mask, xctype, - shls_slice, ao_loc, hermi) + shls_slice, ao_loc, v_hermi) vmat = numpy.stack([vmata, vmatb]) # call swapaxes method to swap last two indices because vmat may be a 3D @@ -636,11 +637,12 @@ def nr_rks_fxc(ni, cell, grids, xc_code, dm0, dms, relativity=0, hermi=0, elif isinstance(kpts, KPoints): kpts = kpts.kpts_ibz + v_hermi = 0 if is_zero(kpts): if isinstance(dms, numpy.ndarray) and dms.dtype == numpy.double: # for real orbitals and real matrix, K_{ia,bj} = K_{ia,jb} # The output matrix v = K*x_{ia} is symmetric - hermi = 1 + v_hermi = 1 if xctype in ('LDA', 'GGA', 'MGGA'): make_rho, nset, nao = ni._gen_rho_evaluator(cell, dms, hermi, False) @@ -660,10 +662,10 @@ def nr_rks_fxc(ni, cell, grids, xc_code, dm0, dms, relativity=0, hermi=0, vxc1 = numpy.einsum('xg,xyg->yg', rho1, _fxc) wv = weight * vxc1 vmat[i] += ni._vxc_mat(cell, ao_k1, wv, mask, xctype, - shls_slice, ao_loc, hermi) + shls_slice, ao_loc, v_hermi) vmat = numpy.stack(vmat) - if hermi == 1: + if v_hermi == 1: # call swapaxes method to swap last two indices because vmat may be a 3D # array (nset,nao,nao) in single k-point mode or a 4D array # (nset,nkpts,nao,nao) in k-points mode @@ -759,10 +761,11 @@ def nr_uks_fxc(ni, cell, grids, xc_code, dm0, dms, relativity=0, hermi=0, kpts = kpts.kpts_ibz dma, dmb = _format_uks_dm(dms) + v_hermi = 0 if is_zero(kpts) and dma.dtype == numpy.double: # for real orbitals and real matrix, K_{ia,bj} = K_{ia,jb} # The output matrix v = K*x_{ia} is symmetric - hermi = 1 + v_hermi = 1 nao = dma.shape[-1] make_rhoa, nset = ni._gen_rho_evaluator(cell, dma, hermi, False)[:2] @@ -782,19 +785,19 @@ def nr_uks_fxc(ni, cell, grids, xc_code, dm0, dms, relativity=0, hermi=0, for i in range(nset): rho1a = make_rhoa(i, ao_k1, mask, xctype) rho1b = make_rhob(i, ao_k1, mask, xctype) - rho1 = numpy.stack([rho1a, rho1b]) if xctype == 'LDA': - vxc1 = numpy.einsum('ag,abyg->byg', rho1, _fxc[:,0,:]) + wv = rho1a * _fxc[0,0] + rho1b * _fxc[1,0] else: - vxc1 = numpy.einsum('axg,axbyg->byg', rho1, _fxc) - wv = weight * vxc1 + wv = numpy.einsum('xg,xbyg->byg', rho1a, _fxc[0]) + wv += numpy.einsum('xg,xbyg->byg', rho1b, _fxc[1]) + wv *= weight vmata[i] += ni._vxc_mat(cell, ao_k1, wv[0], mask, xctype, - shls_slice, ao_loc, hermi) + shls_slice, ao_loc, v_hermi) vmatb[i] += ni._vxc_mat(cell, ao_k1, wv[1], mask, xctype, - shls_slice, ao_loc, hermi) + shls_slice, ao_loc, v_hermi) vmat = numpy.stack([vmata, vmatb]) - if hermi == 1: + if v_hermi == 1: vmat = vmat + vmat.conj().swapaxes(-2,-1) if nset == 1: vmat = vmat.reshape((2,) + dma.shape) @@ -852,12 +855,13 @@ def cache_xc_kernel(ni, cell, grids, xc_code, mo_coeff, mo_occ, spin=0, else: is_rhf = mo_coeff[0].ndim == 1 + with_lapl = False nao = cell.nao_nr() if is_rhf: rho = [] for ao_k1, ao_k2, mask, weight, coords \ in ni.block_loop(cell, grids, nao, ao_deriv, kpts, None, max_memory): - rho.append(ni.eval_rho2(cell, ao_k1, mo_coeff, mo_occ, mask, xctype)) + rho.append(ni.eval_rho2(cell, ao_k1, mo_coeff, mo_occ, mask, xctype, with_lapl)) rho = numpy.hstack(rho) if spin == 1: rho *= .5 @@ -868,8 +872,8 @@ def cache_xc_kernel(ni, cell, grids, xc_code, mo_coeff, mo_occ, spin=0, rhob = [] for ao_k1, ao_k2, mask, weight, coords \ in ni.block_loop(cell, grids, nao, ao_deriv, kpts, None, max_memory): - rhoa.append(ni.eval_rho2(cell, ao_k1, mo_coeff[0], mo_occ[0], mask, xctype)) - rhob.append(ni.eval_rho2(cell, ao_k1, mo_coeff[1], mo_occ[1], mask, xctype)) + rhoa.append(ni.eval_rho2(cell, ao_k1, mo_coeff[0], mo_occ[0], mask, xctype, with_lapl)) + rhob.append(ni.eval_rho2(cell, ao_k1, mo_coeff[1], mo_occ[1], mask, xctype, with_lapl)) rho = numpy.stack([numpy.hstack(rhoa), numpy.hstack(rhob)]) vxc, fxc = ni.eval_xc_eff(xc_code, rho, deriv=2, xctype=xctype)[1:3] return rho, vxc, fxc @@ -898,12 +902,13 @@ def cache_xc_kernel1(ni, cell, grids, xc_code, dm, spin=0, # 0th order density matrix must be hermitian hermi = 1 + with_lapl = False nao = cell.nao_nr() if is_rhf: rho = [] for ao_k1, ao_k2, mask, weight, coords \ in ni.block_loop(cell, grids, nao, ao_deriv, kpts, None, max_memory): - rho.append(ni.eval_rho1(cell, ao_k1, dm, mask, xctype, hermi)) + rho.append(ni.eval_rho1(cell, ao_k1, dm, mask, xctype, hermi, with_lapl)) rho = numpy.hstack(rho) if spin == 1: rho *= .5 @@ -914,8 +919,8 @@ def cache_xc_kernel1(ni, cell, grids, xc_code, dm, spin=0, rhob = [] for ao_k1, ao_k2, mask, weight, coords \ in ni.block_loop(cell, grids, nao, ao_deriv, kpts, None, max_memory): - rhoa.append(ni.eval_rho1(cell, ao_k1, dm[0], mask, xctype, hermi)) - rhob.append(ni.eval_rho1(cell, ao_k1, dm[1], mask, xctype, hermi)) + rhoa.append(ni.eval_rho1(cell, ao_k1, dm[0], mask, xctype, hermi, with_lapl)) + rhob.append(ni.eval_rho1(cell, ao_k1, dm[1], mask, xctype, hermi, with_lapl)) rho = numpy.stack([numpy.hstack(rhoa), numpy.hstack(rhob)]) vxc, fxc = ni.eval_xc_eff(xc_code, rho, deriv=2, xctype=xctype)[1:3] return rho, vxc, fxc @@ -985,11 +990,13 @@ def nr_uks(self, cell, grids, xc_code, dms, relativity=0, hermi=1, hermi, kpt, kpts_band, max_memory, verbose) def _vxc_mat(self, cell, ao, wv, mask, xctype, shls_slice, ao_loc, hermi): + if xctype == 'MGGA': + tau_idx = 4 + wv[tau_idx] *= .5 # *.5 for 1/2 in tau if hermi == 1: # *.5 because mat + mat.T in the caller when hermi=1 wv[0] *= .5 if xctype == 'MGGA': - tau_idx = 4 wv[tau_idx] *= .5 return _vxc_mat(cell, ao, wv, mask, xctype, shls_slice, ao_loc, hermi) @@ -1190,11 +1197,13 @@ def _vxc_mat(self, cell, ao_kpts, wv, mask, xctype, shls_slice, ao_loc, hermi): nao = ao_kpts[0].shape[-1] dtype = numpy.result_type(wv, *ao_kpts) mat = numpy.empty((nkpts,nao,nao), dtype=dtype) + if xctype == 'MGGA': + tau_idx = 4 + wv[tau_idx] *= .5 # *.5 for 1/2 in tau if hermi == 1: # *.5 because mat + mat.T in the caller when hermi=1 wv[0] *= .5 if xctype == 'MGGA': - tau_idx = 4 wv[tau_idx] *= .5 for k in range(nkpts): diff --git a/pyscf/pbc/lib/kpts_helper.py b/pyscf/pbc/lib/kpts_helper.py index 2fbfb76b87..51ed89a098 100644 --- a/pyscf/pbc/lib/kpts_helper.py +++ b/pyscf/pbc/lib/kpts_helper.py @@ -30,7 +30,7 @@ def is_zero(kpt): return abs(np.asarray(kpt)).sum() < KPT_DIFF_TOL -gamma_point = is_zero +is_gamma_point = gamma_point = is_zero def round_to_fbz(kpts, wrap_around=False, tol=KPT_DIFF_TOL): ''' diff --git a/pyscf/pbc/scf/addons.py b/pyscf/pbc/scf/addons.py index 64add227fc..3c09fbabef 100644 --- a/pyscf/pbc/scf/addons.py +++ b/pyscf/pbc/scf/addons.py @@ -526,7 +526,7 @@ def mo_energy_with_exxdiv_none(mf, mo_coeff=None): fockao = mf.get_fock(vhf=vhf, dm=dm) def _get_moe1(C, fao): - return lib.einsum('pi,pq,qi->i', C.conj(), fao, C) + return lib.einsum('pi,pq,qi->i', C.conj(), fao, C).real def _get_moek(kC, kfao): return [_get_moe1(C, fao) for C,fao in zip(kC, kfao)] diff --git a/pyscf/pbc/scf/hf.py b/pyscf/pbc/scf/hf.py index 6b8818c876..510d1b35f4 100644 --- a/pyscf/pbc/scf/hf.py +++ b/pyscf/pbc/scf/hf.py @@ -552,28 +552,12 @@ def kpts(self, x): def build(self, cell=None): # To handle the attribute kpt or kpts loaded from chkfile - if 'kpts' in self.__dict__: - self.kpts = self.__dict__.pop('kpts') - elif 'kpt' in self.__dict__: + if 'kpt' in self.__dict__: self.kpt = self.__dict__.pop('kpt') - # "vcut_ws" precomputing is triggered by pbc.tools.pbc.get_coulG - #if self.exxdiv == 'vcut_ws': - # if self.exx_built is False: - # self.precompute_exx() - # logger.info(self, 'WS alpha = %s', self.exx_alpha) - - kpts = self.kpts if self.rsjk: if not np.all(self.rsjk.kpts == self.kpt): - self.rsjk = self.rsjk.__class__(cell, kpts) - - # for GDF and MDF - with_df = self.with_df - if len(kpts) > 1 and getattr(with_df, '_j_only', False): - logger.warn(self, 'df.j_only cannot be used with k-point HF') - with_df._j_only = False - with_df.reset() + self.rsjk = self.rsjk.__class__(cell, self.kpt) if self.verbose >= logger.WARN: self.check_sanity() diff --git a/pyscf/pbc/scf/khf.py b/pyscf/pbc/scf/khf.py index 8024a45953..631e9d28d9 100644 --- a/pyscf/pbc/scf/khf.py +++ b/pyscf/pbc/scf/khf.py @@ -474,6 +474,33 @@ def mo_coeff_kpts(self): def mo_occ_kpts(self): return self.mo_occ + def build(self, cell=None): + # To handle the attribute kpt or kpts loaded from chkfile + if 'kpts' in self.__dict__: + self.kpts = self.__dict__.pop('kpts') + + # "vcut_ws" precomputing is triggered by pbc.tools.pbc.get_coulG + #if self.exxdiv == 'vcut_ws': + # if self.exx_built is False: + # self.precompute_exx() + # logger.info(self, 'WS alpha = %s', self.exx_alpha) + + kpts = self.kpts + if self.rsjk: + if not np.all(self.rsjk.kpts == kpts): + self.rsjk = self.rsjk.__class__(cell, kpts) + + # for GDF and MDF + with_df = self.with_df + if len(kpts) > 1 and getattr(with_df, '_j_only', False): + logger.warn(self, 'df.j_only cannot be used with k-point HF') + with_df._j_only = False + with_df.reset() + + if self.verbose >= logger.WARN: + self.check_sanity() + return self + def dump_flags(self, verbose=None): mol_hf.SCF.dump_flags(self, verbose) logger.info(self, '\n') diff --git a/pyscf/pbc/tdscf/krhf.py b/pyscf/pbc/tdscf/krhf.py index 7c49538ca8..9323e877c1 100644 --- a/pyscf/pbc/tdscf/krhf.py +++ b/pyscf/pbc/tdscf/krhf.py @@ -26,10 +26,11 @@ from pyscf.lib import linalg_helper from pyscf.lib import logger from pyscf.tdscf import rhf +from pyscf.tdscf._lr_eig import eigh as lr_eigh, eig as lr_eig from pyscf.pbc import scf from pyscf.pbc.tdscf.rhf import TDBase from pyscf.pbc.scf import _response_functions # noqa -from pyscf.pbc.lib.kpts_helper import gamma_point, get_kconserv_ria +from pyscf.pbc.lib.kpts_helper import is_gamma_point, get_kconserv_ria, conj_mapping from pyscf.pbc.df.df_ao2mo import warn_pbc2d_eri from pyscf.pbc import df as pbcdf from pyscf.data import nist @@ -37,8 +38,139 @@ REAL_EIG_THRESHOLD = getattr(__config__, 'pbc_tdscf_rhf_TDDFT_pick_eig_threshold', 1e-3) +def get_ab(mf, kshift=0): + r'''A and B matrices for TDDFT response function. + + A[i,a,j,b] = \delta_{ab}\delta_{ij}(E_a - E_i) + (ia||bj) + B[i,a,j,b] = (ia||jb) + + Ref: Chem Phys Lett, 256, 454 + + Kwargs: + kshift : integer + The index of the k-point that represents the transition between + k-points in the excitation coefficients. + ''' + cell = mf.cell + mo_energy = scf.addons.mo_energy_with_exxdiv_none(mf) + mo = numpy.asarray(mf.mo_coeff) + mo_occ = numpy.asarray(mf.mo_occ) + kpts = mf.kpts + nkpts, nao, nmo = mo.shape + noccs = numpy.count_nonzero(mo_occ==2, axis=1) + nocc = noccs[0] + nvir = nmo - nocc + assert all(noccs == nocc) + orbo = mo[:,:,:nocc] + orbv = mo[:,:,nocc:] + + kconserv = get_kconserv_ria(cell, kpts)[kshift] + e_ia = numpy.asarray(_get_e_ia(mo_energy, mo_occ, kconserv)).astype(mo.dtype) + a = numpy.diag(e_ia.ravel()).reshape(nkpts,nocc,nvir,nkpts,nocc,nvir) + b = numpy.zeros_like(a) + weight = 1./nkpts + + def add_hf_(a, b, hyb=1): + eri = mf.with_df.ao2mo_7d([orbo,mo,mo,mo], kpts) + eri *= weight + eri = eri.reshape(nkpts,nkpts,nkpts,nocc,nmo,nmo,nmo) + for ki, ka in enumerate(kconserv): + for kj, kb in enumerate(kconserv): + a[ki,:,:,kj] += numpy.einsum('iabj->iajb', eri[ki,ka,kb,:nocc,nocc:,nocc:,:nocc]) * 2 + a[ki,:,:,kj] -= numpy.einsum('ijba->iajb', eri[ki,kj,kb,:nocc,:nocc,nocc:,nocc:]) * hyb + + for kb, kj in enumerate(kconserv): + b[ki,:,:,kj] += numpy.einsum('iajb->iajb', eri[ki,ka,kj,:nocc,nocc:,:nocc,nocc:]) * 2 + b[ki,:,:,kj] -= numpy.einsum('ibja->iajb', eri[ki,kb,kj,:nocc,nocc:,:nocc,nocc:]) * hyb + + if isinstance(mf, scf.hf.KohnShamDFT): + assert kshift == 0 + ni = mf._numint + omega, alpha, hyb = ni.rsh_and_hybrid_coeff(mf.xc, cell.spin) + + add_hf_(a, b, hyb) + if omega != 0: # For RSH + raise NotImplementedError + + xctype = ni._xc_type(mf.xc) + dm0 = mf.make_rdm1(mo, mo_occ) + make_rho = ni._gen_rho_evaluator(cell, dm0, hermi=1, with_lapl=False)[0] + mem_now = lib.current_memory()[0] + max_memory = max(2000, mf.max_memory*.8-mem_now) + cmap = conj_mapping(cell, kpts) + + if xctype == 'LDA': + ao_deriv = 0 + for ao, _, mask, weight, coords \ + in ni.block_loop(cell, mf.grids, nao, ao_deriv, kpts, None, max_memory): + rho = make_rho(0, ao, mask, xctype) + fxc = ni.eval_xc_eff(mf.xc, rho, deriv=2, xctype=xctype)[2] + wfxc = fxc[0,0] * weight + + rho_o = lib.einsum('krp,kpi->kri', ao, orbo) + rho_v = lib.einsum('krp,kpi->kri', ao, orbv) + rho_ov = numpy.einsum('kri,kra->kria', rho_o, rho_v) + w_ov = numpy.einsum('kria,r->kria', rho_ov, wfxc) * (2/nkpts) + rho_vo = rho_ov.conj()[cmap] + a += lib.einsum('kria,lrjb->kialjb', w_ov, rho_vo) + b += lib.einsum('kria,lrjb->kialjb', w_ov, rho_ov) + + elif xctype == 'GGA': + ao_deriv = 1 + for ao, _, mask, weight, coords \ + in ni.block_loop(cell, mf.grids, nao, ao_deriv, kpts, None, max_memory): + rho = make_rho(0, ao, mask, xctype) + fxc = ni.eval_xc_eff(mf.xc, rho, deriv=2, xctype=xctype)[2] + wfxc = fxc * weight + rho_o = lib.einsum('kxrp,kpi->kxri', ao, orbo) + rho_v = lib.einsum('kxrp,kpi->kxri', ao, orbv) + rho_ov = numpy.einsum('kxri,kra->kxria', rho_o, rho_v[:,0]) + rho_ov[:,1:4] += numpy.einsum('kri,kxra->kxria', rho_o[:,0], rho_v[:,1:4]) + w_ov = numpy.einsum('xyr,kxria->kyria', wfxc, rho_ov) * (2/nkpts) + rho_vo = rho_ov.conj()[cmap] + a += lib.einsum('kxria,lxrjb->kialjb', w_ov, rho_vo) + b += lib.einsum('kxria,lxrjb->kialjb', w_ov, rho_ov) + + elif xctype == 'HF': + pass + + elif xctype == 'NLC': + raise NotImplementedError('NLC') + + elif xctype == 'MGGA': + ao_deriv = 1 + for ao, _, mask, weight, coords \ + in ni.block_loop(cell, mf.grids, nao, ao_deriv, kpts, None, max_memory): + rho = make_rho(0, ao, mask, xctype) + fxc = ni.eval_xc_eff(mf.xc, rho, deriv=2, xctype=xctype)[2] + wfxc = fxc * weight + rho_o = lib.einsum('kxrp,kpi->kxri', ao, orbo) + rho_v = lib.einsum('kxrp,kpi->kxri', ao, orbv) + rho_ov = numpy.einsum('kxri,kra->kxria', rho_o, rho_v[:,0]) + rho_ov[:,1:4] += numpy.einsum('kri,kxra->kxria', rho_o[:,0], rho_v[:,1:4]) + tau_ov = numpy.einsum('kxri,kxra->kria', rho_o[:,1:4], rho_v[:,1:4]) * .5 + rho_ov = numpy.concatenate([rho_ov, tau_ov[:,numpy.newaxis]], axis=1) + w_ov = numpy.einsum('xyr,kxria->kyria', wfxc, rho_ov) * (2/nkpts) + rho_vo = rho_ov.conj()[cmap] + a += lib.einsum('kxria,lxrjb->kialjb', w_ov, rho_vo) + b += lib.einsum('kxria,lxrjb->kialjb', w_ov, rho_ov) + else: + add_hf_(a, b) + + return a, b + class KTDBase(TDBase): - _keys = {'kconserv', 'kshift_lst'} + ''' + Attributes: + kshift_lst : list of integers + Each element in the list is the index of the k-point that + represents the transition between k-points in the excitation + coefficients. + ''' + + conv_tol = getattr(__config__, 'pbc_tdscf_rhf_TDA_conv_tol', 1e-4) + + _keys = {'kshift_lst'} def __init__(self, mf, kshift_lst=None): assert isinstance(mf, scf.khf.KSCF) @@ -46,8 +178,6 @@ def __init__(self, mf, kshift_lst=None): warn_pbc2d_eri(mf) if kshift_lst is None: kshift_lst = [0] - - self.kconserv = get_kconserv_ria(mf.cell, mf.kpts) self.kshift_lst = kshift_lst def dump_flags(self, verbose=None): @@ -67,7 +197,6 @@ def dump_flags(self, verbose=None): log.info('conv_tol = %g', self.conv_tol) log.info('eigh lindep = %g', self.lindep) log.info('eigh level_shift = %g', self.level_shift) - log.info('eigh max_space = %d', self.max_space) log.info('eigh max_cycle = %d', self.max_cycle) log.info('chkfile = %s', self.chkfile) log.info('max_memory %d MB (current use %d MB)', @@ -98,14 +227,22 @@ def _finalize(self): get_nto = lib.invalid_method('get_nto') class TDA(KTDBase): - conv_tol = getattr(__config__, 'pbc_tdscf_rhf_TDA_conv_tol', 1e-6) - def gen_vind(self, mf, kshift): - # exxdiv corrections are kept in hdiag while excluding them when calling - # the contractions between two-electron integrals and X/Y amplitudes. - # See also the relevant comments in function pbc.tdscf.rhf.TDA.gen_vind + @lib.with_doc(get_ab.__doc__) + def get_ab(self, mf=None, kshift=0): + if mf is None: mf = self._scf + return get_ab(mf, kshift) + + def gen_vind(self, mf, kshift=0): + '''Compute Ax + + Kwargs: + kshift : integer + The index of the k-point that represents the transition between + k-points in the excitation coefficients. + ''' singlet = self.singlet - kconserv = self.kconserv[kshift] + kconserv = get_kconserv_ria(mf.cell, mf.kpts)[kshift] mo_coeff = mf.mo_coeff mo_occ = mf.mo_occ @@ -114,7 +251,7 @@ def gen_vind(self, mf, kshift): occidx = [numpy.where(mo_occ[k]==2)[0] for k in range(nkpts)] viridx = [numpy.where(mo_occ[k]==0)[0] for k in range(nkpts)] orbo = [mo_coeff[k][:,occidx[k]] for k in range(nkpts)] - orbv = [mo_coeff[kconserv[k]][:,viridx[kconserv[k]]] for k in range(nkpts)] + orbv = [mo_coeff[k][:,viridx[k]] for k in range(nkpts)] e_ia = _get_e_ia(scf.addons.mo_energy_with_exxdiv_none(mf), mo_occ, kconserv) hdiag = numpy.hstack([x.ravel() for x in e_ia]) @@ -127,10 +264,10 @@ def vind(zs): z1s = [_unpack(z, mo_occ, kconserv) for z in zs] dmov = numpy.empty((nz,nkpts,nao,nao), dtype=numpy.complex128) for i in range(nz): - for k in range(nkpts): + for k, kp in enumerate(kconserv): # *2 for double occupancy dm1 = z1s[i][k] * 2 - dmov[i,k] = reduce(numpy.dot, (orbo[k], dm1, orbv[k].conj().T)) + dmov[i,k] = reduce(numpy.dot, (orbo[k], dm1, orbv[kp].conj().T)) with lib.temporary_env(mf, exxdiv=None): v1ao = vresp(dmov, kshift) @@ -138,8 +275,8 @@ def vind(zs): for i in range(nz): dm1 = z1s[i] v1 = [] - for k in range(nkpts): - v1vo = reduce(numpy.dot, (orbo[k].conj().T, v1ao[i,k], orbv[k])) + for k, kp in enumerate(kconserv): + v1vo = reduce(numpy.dot, (orbo[k].conj().T, v1ao[i,k], orbv[kp])) v1vo += e_ia[k] * dm1[k] v1.append(v1vo.ravel()) v1s.append( numpy.concatenate(v1) ) @@ -151,13 +288,13 @@ def init_guess(self, mf, kshift, nstates=None): mo_energy = mf.mo_energy mo_occ = mf.mo_occ - kconserv = self.kconserv[kshift] + kconserv = get_kconserv_ria(mf.cell, mf.kpts)[kshift] e_ia = numpy.concatenate( [x.reshape(-1) for x in _get_e_ia(mo_energy, mo_occ, kconserv)] ) nov = e_ia.size nstates = min(nstates, nov) - e_threshold = numpy.sort(e_ia)[nstates-1] + e_threshold = numpy.partition(e_ia, nstates-1)[nstates-1] e_threshold += self.deg_eia_thresh idx = numpy.where(e_ia <= e_threshold)[0] @@ -188,13 +325,12 @@ def pickeig(w, v, nroots, envs): return w[idx], v[:,idx], idx log = logger.Logger(self.stdout, self.verbose) - precision = self.cell.precision * 1e-2 self.converged = [] self.e = [] self.xy = [] for i,kshift in enumerate(self.kshift_lst): - kconserv = self.kconserv[kshift] + kconserv = get_kconserv_ria(mf.cell, mf.kpts)[kshift] vind, hdiag = self.gen_vind(self._scf, kshift) precond = self.get_precond(hdiag) @@ -204,14 +340,10 @@ def pickeig(w, v, nroots, envs): else: x0k = x0[i] - converged, e, x1 = \ - lib.davidson1(vind, x0k, precond, - tol=self.conv_tol, - max_cycle=self.max_cycle, - nroots=self.nstates, lindep=self.lindep, - max_space=self.max_space, pick=pickeig, - fill_heff=purify_krlyov_heff(precision, 0, log), - verbose=self.verbose) + converged, e, x1 = lr_eigh( + vind, x0k, precond, tol_residual=self.conv_tol, lindep=self.lindep, + nroots=self.nstates, pick=pickeig, max_cycle=self.max_cycle, + max_memory=self.max_memory, verbose=log) self.converged.append( converged ) self.e.append( e ) # 1/sqrt(2) because self.x is for alpha excitation amplitude and 2(X^+*X) = 1 @@ -223,14 +355,21 @@ def pickeig(w, v, nroots, envs): CIS = KTDA = TDA -class TDHF(TDA): - def gen_vind(self, mf, kshift): +class TDHF(KTDBase): + + @lib.with_doc(get_ab.__doc__) + def get_ab(self, mf=None, kshift=0): + if mf is None: mf = self._scf + return get_ab(mf, kshift) + + def gen_vind(self, mf, kshift=0): ''' [ A B ][X] [-B* -A*][Y] ''' + assert kshift == 0 + singlet = self.singlet - kconserv = self.kconserv[kshift] mo_coeff = mf.mo_coeff mo_occ = mf.mo_occ @@ -239,7 +378,9 @@ def gen_vind(self, mf, kshift): occidx = [numpy.where(mo_occ[k]==2)[0] for k in range(nkpts)] viridx = [numpy.where(mo_occ[k]==0)[0] for k in range(nkpts)] orbo = [mo_coeff[k][:,occidx[k]] for k in range(nkpts)] - orbv = [mo_coeff[kconserv[k]][:,viridx[kconserv[k]]] for k in range(nkpts)] + orbv = [mo_coeff[k][:,viridx[k]] for k in range(nkpts)] + + kconserv = numpy.arange(nkpts) e_ia = _get_e_ia(scf.addons.mo_energy_with_exxdiv_none(mf), mo_occ, kconserv) hdiag = numpy.hstack([x.ravel() for x in e_ia]) tot_x = hdiag.size @@ -253,14 +394,14 @@ def vind(xys): nz = len(xys) z1xs = [_unpack(xy[:tot_x], mo_occ, kconserv) for xy in xys] z1ys = [_unpack(xy[tot_x:], mo_occ, kconserv) for xy in xys] - dmov = numpy.empty((nz,nkpts,nao,nao), dtype=numpy.complex128) + dmov = numpy.zeros((nz,nkpts,nao,nao), dtype=numpy.complex128) for i in range(nz): for k in range(nkpts): # *2 for double occupancy dmx = z1xs[i][k] * 2 dmy = z1ys[i][k] * 2 - dmov[i,k] = reduce(numpy.dot, (orbo[k], dmx, orbv[k].T.conj())) - dmov[i,k]+= reduce(numpy.dot, (orbv[k], dmy.T, orbo[k].T.conj())) + dmov[i,k] += reduce(numpy.dot, (orbo[k], dmx , orbv[k].T.conj())) + dmov[i,k] += reduce(numpy.dot, (orbv[k], dmy.T, orbo[k].T.conj())) with lib.temporary_env(mf, exxdiv=None): v1ao = vresp(dmov, kshift) # = Xjb + Yjb @@ -268,8 +409,8 @@ def vind(xys): for i in range(nz): dmx = z1xs[i] dmy = z1ys[i] - v1xs = [] - v1ys = [] + v1xs = [0] * nkpts + v1ys = [0] * nkpts for k in range(nkpts): # AX + BY # = Xjb + Yjb @@ -279,11 +420,10 @@ def vind(xys): # = Xjb + Yjb # = ( Xjb + Yjb) Cma* Cni v1y = reduce(numpy.dot, (orbv[k].T.conj(), v1ao[i,k], orbo[k])).T - v1x+= e_ia[k] * dmx[k] - v1y+= e_ia[k].conj() * dmy[k] - v1xs.append(v1x.ravel()) - v1ys.append(-v1y.ravel()) - # v1s += v1xs + v1ys + v1x += e_ia[k] * dmx[k] + v1y += e_ia[k].conj() * dmy[k] + v1xs[k] += v1x.ravel() + v1ys[k] -= v1y.ravel() v1s.append( numpy.concatenate(v1xs + v1ys) ) return lib.asarray(v1s).reshape(nz,-1) return vind, hdiag @@ -291,7 +431,9 @@ def vind(xys): def init_guess(self, mf, kshift, nstates=None): x0 = TDA.init_guess(self, mf, kshift, nstates) y0 = numpy.zeros_like(x0) - return numpy.asarray(numpy.block([[x0, y0], [y0, x0.conj()]])) + return numpy.hstack([x0, y0]) + + get_precond = rhf.TDHF.get_precond def kernel(self, x0=None): '''TDHF diagonalization with non-Hermitian eigenvalue solver @@ -305,9 +447,21 @@ def kernel(self, x0=None): mf = self._scf mo_occ = mf.mo_occ - real_system = (gamma_point(self._scf.kpts) and + real_system = (is_gamma_point(self._scf.kpts) and self._scf.mo_coeff[0].dtype == numpy.double) + if any(k != 0 for k in self.kshift_lst): + # It's not clear how to define the Y matrix for kshift!=0 . + # When the A tensor is constructed against the X(kshift) matrix, + # the diagonal terms e_ia are calculated as e_i[k] - e_k[k+kshift]. + # Given the k-conserve relation in the A tensor, the j-b indices in + # the A tensor should follow j[k'], b[k'+kshift]. This leads to the + # j-b indices in the B tensor being defined as (j[k'+shift], b[k']). + # To form the square A-B-B-A matrix, the diagonal terms for the + # -A* part need to be constructed as e_i[k+kshift] - e_a[k], which + # conflict to the diagonal terms of the A tensor. + raise RuntimeError('kshift != 0 for TDHF') + # We only need positive eigenvalues def pickeig(w, v, nroots, envs): realidx = numpy.where((abs(w.imag) < REAL_EIG_THRESHOLD) & @@ -315,8 +469,6 @@ def pickeig(w, v, nroots, envs): return lib.linalg_helper._eigs_cmplx2real(w, v, realidx, real_system) log = logger.Logger(self.stdout, self.verbose) - precision = self.cell.precision * 1e-2 - hermi = 0 def norm_xy(z, kconserv): x, y = z.reshape(2,-1) @@ -330,7 +482,7 @@ def norm_xy(z, kconserv): self.e = [] self.xy = [] for i,kshift in enumerate(self.kshift_lst): - kconserv = self.kconserv[kshift] + kconserv = get_kconserv_ria(mf.cell, mf.kpts)[kshift] vind, hdiag = self.gen_vind(self._scf, kshift) precond = self.get_precond(hdiag) @@ -340,15 +492,10 @@ def norm_xy(z, kconserv): else: x0k = x0[i] - converged, e, x1 = \ - lib.davidson_nosym1(vind, x0k, precond, - tol=self.conv_tol, - max_cycle=self.max_cycle, - nroots=self.nstates, - lindep=self.lindep, - max_space=self.max_space, pick=pickeig, - fill_heff=purify_krlyov_heff(precision, hermi, log), - verbose=self.verbose) + converged, e, x1 = lr_eig( + vind, x0k, precond, tol_residual=self.conv_tol, lindep=self.lindep, + nroots=self.nstates, pick=pickeig, max_cycle=self.max_cycle, + max_memory=self.max_memory, verbose=log) self.converged.append( converged ) self.e.append( e ) self.xy.append( [norm_xy(z, kconserv) for z in x1] ) @@ -372,87 +519,16 @@ def _get_e_ia(mo_energy, mo_occ, kconserv=None): def _unpack(vo, mo_occ, kconserv): z = [] p1 = 0 - for k, occ in enumerate(mo_occ): - no = numpy.count_nonzero(occ > 0) - no1 = numpy.count_nonzero(mo_occ[kconserv[k]] > 0) - nv = occ.size - no1 + no_kpts = [numpy.count_nonzero(occ) for occ in mo_occ] + for k, no in enumerate(no_kpts): + kp = kconserv[k] + nv = mo_occ[kp].size - no_kpts[kp] p0, p1 = p1, p1 + no * nv z.append(vo[p0:p1].reshape(no,nv)) return z -def purify_krlyov_heff(precision, hermi, log): - def fill_heff(heff, xs, ax, xt, axt, dot): - if hermi == 1: - heff = linalg_helper._fill_heff_hermitian(heff, xs, ax, xt, axt, dot) - else: - heff = linalg_helper._fill_heff(heff, xs, ax, xt, axt, dot) - space = len(axt) - # TODO: PBC integrals has larger errors than molecule systems. - # purify the effective Hamiltonian with symmetry and other - # possible conditions. - if abs(heff[:space,:space].imag).max() < precision: - log.debug('Remove imaginary part of the Krylov space effective Hamiltonian') - heff[:space,:space].imag = 0 - return heff - return fill_heff - scf.khf.KRHF.TDA = lib.class_as_method(KTDA) scf.khf.KRHF.TDHF = lib.class_as_method(KTDHF) scf.krohf.KROHF.TDA = None scf.krohf.KROHF.TDHF = None - - -if __name__ == '__main__': - from pyscf.pbc import gto - from pyscf.pbc import scf - from pyscf.pbc import df - cell = gto.Cell() - cell.unit = 'B' - cell.atom = ''' - C 0. 0. 0. - C 1.68506879 1.68506879 1.68506879 - ''' - cell.a = ''' - 0. 3.37013758 3.37013758 - 3.37013758 0. 3.37013758 - 3.37013758 3.37013758 0. - ''' - - cell.basis = 'gth-szv' - cell.pseudo = 'gth-pade' - cell.mesh = [25]*3 - cell.build() - mf = scf.KRHF(cell, cell.make_kpts([2,1,1])).set(exxdiv=None) - #mf.with_df = df.MDF(cell, cell.make_kpts([2,1,1])) - #mf.with_df.auxbasis = 'weigend' - #mf.with_df._cderi = 'eri3d-mdf.h5' - #mf.with_df.build(with_j3c=False) - mf.run() - #mesh=9 -8.65192427146353 - #mesh=12 -8.65192352289817 - #mesh=15 -8.6519235231529 - #MDF mesh=5 -8.6519301815144 - - td = TDA(mf) - td.verbose = 5 - print(td.kernel()[0][0] * 27.2114) - #mesh=9 [ 6.0073749 6.09315355 6.3479901 ] - #mesh=12 [ 6.00253282 6.09317929 6.34799109] - #mesh=15 [ 6.00253396 6.09317949 6.34799109] - #MDF mesh=5 [ 6.09317489 6.09318265 6.34798637] - - #from pyscf.pbc import tools - #scell = tools.super_cell(cell, [2,1,1]) - #mf = scf.RHF(scell).run() - #td = rhf.TDA(mf) - #td.verbose = 5 - #print(td.kernel()[0] * 27.2114) - - td = TDHF(mf) - td.verbose = 5 - print(td.kernel()[0][0] * 27.2114) - #mesh=9 [ 6.03860914 6.21664545 8.20305225] - #mesh=12 [ 6.03868259 6.03860343 6.2167623 ] - #mesh=15 [ 6.03861321 6.03861324 6.21675868] - #MDF mesh=5 [ 6.03861693 6.03861775 6.21675694] diff --git a/pyscf/pbc/tdscf/krks.py b/pyscf/pbc/tdscf/krks.py index 2c734b3589..336b47d9a4 100644 --- a/pyscf/pbc/tdscf/krks.py +++ b/pyscf/pbc/tdscf/krks.py @@ -57,47 +57,3 @@ def _rebuild_df(td): dft.kroks.KROKS.TDA = None dft.kroks.KROKS.TDHF = None dft.kroks.KROKS.TDDFT = None - - -if __name__ == '__main__': - from pyscf.pbc import gto - from pyscf.pbc import dft - cell = gto.Cell() - cell.unit = 'B' - cell.atom = ''' - C 0. 0. 0. - C 1.68506879 1.68506879 1.68506879 - ''' - cell.a = ''' - 0. 3.37013758 3.37013758 - 3.37013758 0. 3.37013758 - 3.37013758 3.37013758 0. - ''' - cell.basis = 'gth-szv' - cell.pseudo = 'gth-pade' - cell.mesh = [25]*3 - cell.build() - - mf = dft.KRKS(cell, cell.make_kpts([2,1,1])) - #mf.with_df = df.MDF(cell, cell.make_kpts([2,1,1])) - #mf.with_df.auxbasis = 'weigend' - #mf.with_df._cderi = 'eri3d-mdf.h5' - #mf.with_df.build(with_j3c=False) - mf.xc = 'lda,' - mf.kernel() -#mesh=12 -10.3077341607895 -#mesh=5 -10.3086623157515 - - td = TDDFT(mf) - td.nstates = 5 - td.verbose = 5 - print(td.kernel()[0][0] * 27.2114) -#mesh=12 [ 6.08108297 6.10231481 6.10231478 6.38355803 6.38355804] -#MDF mesh=5 [ 6.07919157 6.10251718 6.10253961 6.37202499 6.37565246] - - td = TDA(mf) - td.singlet = False - td.verbose = 5 - print(td.kernel()[0][0] * 27.2114) -#mesh=12 [ 4.01539192 5.1750807 5.17508071] -#MDF mesh=5 [ 4.01148649 5.18043397 5.18043459] diff --git a/pyscf/pbc/tdscf/kuhf.py b/pyscf/pbc/tdscf/kuhf.py index 7e98bdbd33..30e476635b 100644 --- a/pyscf/pbc/tdscf/kuhf.py +++ b/pyscf/pbc/tdscf/kuhf.py @@ -21,20 +21,29 @@ from pyscf import lib from pyscf.lib import logger from pyscf.tdscf import uhf +from pyscf.tdscf._lr_eig import eigh as lr_eigh, eig as lr_eig from pyscf.pbc import scf -from pyscf.pbc.tdscf.krhf import KTDBase, _get_e_ia, purify_krlyov_heff -from pyscf.pbc.lib.kpts_helper import gamma_point +from pyscf.pbc.tdscf.krhf import KTDBase, _get_e_ia +from pyscf.pbc.lib.kpts_helper import is_gamma_point, get_kconserv_ria from pyscf.pbc.scf import _response_functions # noqa from pyscf import __config__ REAL_EIG_THRESHOLD = getattr(__config__, 'pbc_tdscf_uhf_TDDFT_pick_eig_threshold', 1e-3) class TDA(KTDBase): - conv_tol = getattr(__config__, 'pbc_tdscf_rhf_TDA_conv_tol', 1e-6) - def gen_vind(self, mf, kshift): - '''Compute Ax''' - kconserv = self.kconserv[kshift] + def get_ab(self, mf=None, kshift=0): + raise NotImplementedError + + def gen_vind(self, mf, kshift=0): + '''Compute Ax + + Kwargs: + kshift : integer + The index of the k-point that represents the transition between + k-points in the excitation coefficients. + ''' + kconserv = get_kconserv_ria(mf.cell, mf.kpts)[kshift] mo_coeff = mf.mo_coeff mo_occ = mf.mo_occ nkpts = len(mo_occ[0]) @@ -94,7 +103,7 @@ def init_guess(self, mf, kshift, nstates=None): mo_energy = mf.mo_energy mo_occ = mf.mo_occ - kconserv = self.kconserv[kshift] + kconserv = get_kconserv_ria(mf.cell, mf.kpts)[kshift] e_ia_a = _get_e_ia(mo_energy[0], mo_occ[0], kconserv) e_ia_b = _get_e_ia(mo_energy[1], mo_occ[1], kconserv) e_ia = numpy.hstack([x.ravel() for x in (e_ia_a + e_ia_b)]) @@ -127,14 +136,12 @@ def pickeig(w, v, nroots, envs): return w[idx], v[:,idx], idx log = logger.Logger(self.stdout, self.verbose) - precision = self.cell.precision * 1e-2 - hermi = 1 self.converged = [] self.e = [] self.xy = [] for i,kshift in enumerate(self.kshift_lst): - kconserv = self.kconserv[kshift] + kconserv = get_kconserv_ria(mf.cell, mf.kpts)[kshift] vind, hdiag = self.gen_vind(self._scf, kshift) precond = self.get_precond(hdiag) @@ -144,15 +151,10 @@ def pickeig(w, v, nroots, envs): else: x0k = x0[i] - converged, e, x1 = \ - lib.davidson1(vind, x0k, precond, - tol=self.conv_tol, - max_cycle=self.max_cycle, - nroots=self.nstates, - lindep=self.lindep, - max_space=self.max_space, pick=pickeig, - fill_heff=purify_krlyov_heff(precision, hermi, log), - verbose=self.verbose) + converged, e, x1 = lr_eigh( + vind, x0k, precond, tol_residual=self.conv_tol, lindep=self.lindep, + nroots=self.nstates, pick=pickeig, max_cycle=self.max_cycle, + max_memory=self.max_memory, verbose=log) self.converged.append( converged ) self.e.append( e ) self.xy.append( [(_unpack(xi, mo_occ, kconserv), # (X_alpha, X_beta) @@ -165,9 +167,14 @@ def pickeig(w, v, nroots, envs): CIS = KTDA = TDA -class TDHF(TDA): - def gen_vind(self, mf, kshift): - kconserv = self.kconserv[kshift] +class TDHF(KTDBase): + + def get_ab(self, mf=None, kshift=0): + raise NotImplementedError + + def gen_vind(self, mf, kshift=0): + assert kshift == 0 + mo_coeff = mf.mo_coeff mo_occ = mf.mo_occ nkpts = len(mo_occ[0]) @@ -178,9 +185,10 @@ def gen_vind(self, mf, kshift): viridxb = [numpy.where(mo_occ[1][k]==0)[0] for k in range(nkpts)] orboa = [mo_coeff[0][k][:,occidxa[k]] for k in range(nkpts)] orbob = [mo_coeff[1][k][:,occidxb[k]] for k in range(nkpts)] - orbva = [mo_coeff[0][kconserv[k]][:,viridxa[kconserv[k]]] for k in range(nkpts)] - orbvb = [mo_coeff[1][kconserv[k]][:,viridxb[kconserv[k]]] for k in range(nkpts)] + orbva = [mo_coeff[0][k][:,viridxa[k]] for k in range(nkpts)] + orbvb = [mo_coeff[1][k][:,viridxb[k]] for k in range(nkpts)] + kconserv = numpy.arange(nkpts) moe = scf.addons.mo_energy_with_exxdiv_none(mf) e_ia_a = _get_e_ia(moe[0], mo_occ[0], kconserv) e_ia_b = _get_e_ia(moe[1], mo_occ[1], kconserv) @@ -219,23 +227,23 @@ def vind(xys): for i in range(nz): xa, xb = x1s[i] ya, yb = y1s[i] - v1xsa = [] - v1xsb = [] - v1ysa = [] - v1ysb = [] + v1xsa = [0] * nkpts + v1xsb = [0] * nkpts + v1ysa = [0] * nkpts + v1ysb = [0] * nkpts for k in range(nkpts): v1xa = reduce(numpy.dot, (orboa[k].conj().T, v1ao[0,i,k], orbva[k])) v1xb = reduce(numpy.dot, (orbob[k].conj().T, v1ao[1,i,k], orbvb[k])) v1ya = reduce(numpy.dot, (orbva[k].conj().T, v1ao[0,i,k], orboa[k])).T v1yb = reduce(numpy.dot, (orbvb[k].conj().T, v1ao[1,i,k], orbob[k])).T - v1xa+= e_ia_a[k] * xa[k] - v1xb+= e_ia_b[k] * xb[k] - v1ya+= e_ia_a[k].conj() * ya[k] - v1yb+= e_ia_b[k].conj() * yb[k] - v1xsa.append(v1xa.ravel()) - v1xsb.append(v1xb.ravel()) - v1ysa.append(-v1ya.ravel()) - v1ysb.append(-v1yb.ravel()) + v1xa += e_ia_a[k] * xa[k] + v1xb += e_ia_b[k] * xb[k] + v1ya += e_ia_a[k].conj() * ya[k] + v1yb += e_ia_b[k].conj() * yb[k] + v1xsa[k] += v1xa.ravel() + v1xsb[k] += v1xb.ravel() + v1ysa[k] += -v1ya.ravel() + v1ysb[k] += -v1yb.ravel() v1s.append( numpy.concatenate(v1xsa + v1xsb + v1ysa + v1ysb) ) return numpy.hstack(v1s).reshape(nz,-1) @@ -244,7 +252,9 @@ def vind(xys): def init_guess(self, mf, kshift, nstates=None, wfnsym=None): x0 = TDA.init_guess(self, mf, kshift, nstates) y0 = numpy.zeros_like(x0) - return numpy.asarray(numpy.block([[x0, y0], [y0, x0.conj()]])) + return numpy.hstack([x0, y0]) + + get_precond = uhf.TDHF.get_precond def kernel(self, x0=None): '''TDHF diagonalization with non-Hermitian eigenvalue solver @@ -258,23 +268,23 @@ def kernel(self, x0=None): mf = self._scf mo_occ = mf.mo_occ - real_system = (gamma_point(self._scf.kpts) and + real_system = (is_gamma_point(self._scf.kpts) and self._scf.mo_coeff[0][0].dtype == numpy.double) + if any(k != 0 for k in self.kshift_lst): + raise RuntimeError('kshift != 0 for TDHF') + # We only need positive eigenvalues def pickeig(w, v, nroots, envs): realidx = numpy.where((abs(w.imag) < REAL_EIG_THRESHOLD) & (w.real > self.positive_eig_threshold))[0] return lib.linalg_helper._eigs_cmplx2real(w, v, realidx, real_system) - log = logger.Logger(self.stdout, self.verbose) - precision = self.cell.precision * 1e-2 - self.converged = [] self.e = [] self.xy = [] for i,kshift in enumerate(self.kshift_lst): - kconserv = self.kconserv[kshift] + kconserv = get_kconserv_ria(mf.cell, mf.kpts)[kshift] vind, hdiag = self.gen_vind(self._scf, kshift) precond = self.get_precond(hdiag) @@ -284,15 +294,10 @@ def pickeig(w, v, nroots, envs): else: x0k = x0[i] - converged, w, x1 = \ - lib.davidson_nosym1(vind, x0k, precond, - tol=self.conv_tol, - max_cycle=self.max_cycle, - nroots=self.nstates, - lindep=self.lindep, - max_space=self.max_space, pick=pickeig, - fill_heff=purify_krlyov_heff(precision, 0, log), - verbose=self.verbose) + converged, w, x1 = lr_eig( + vind, x0k, precond, tol_residual=self.conv_tol, lindep=self.lindep, + nroots=self.nstates, pick=pickeig, max_cycle=self.max_cycle, + max_memory=self.max_memory, verbose=log) self.converged.append( converged ) e = [] @@ -336,55 +341,3 @@ def _unpack(vo, mo_occ, kconserv): scf.kuhf.KUHF.TDA = lib.class_as_method(KTDA) scf.kuhf.KUHF.TDHF = lib.class_as_method(KTDHF) - - -if __name__ == '__main__': - from pyscf.pbc import gto - from pyscf.pbc import scf - from pyscf.pbc import df - cell = gto.Cell() - cell.unit = 'B' - cell.atom = ''' - C 0. 0. 0. - C 1.68506879 1.68506879 1.68506879 - ''' - cell.a = ''' - 0. 3.37013758 3.37013758 - 3.37013758 0. 3.37013758 - 3.37013758 3.37013758 0. - ''' - - cell.basis = 'gth-szv' - cell.pseudo = 'gth-pade' - cell.mesh = [37]*3 - cell.build() - mf = scf.KUHF(cell, cell.make_kpts([2,1,1])).set(exxdiv=None) -# mf.with_df = df.DF(cell, cell.make_kpts([2,1,1])) -# mf.with_df.auxbasis = 'weigend' -# mf.with_df._cderi = 'eri3d-df.h5' -# mf.with_df.build(with_j3c=False) - mf.run() - - td = TDA(mf) - td.verbose = 5 - td.nstates = 5 - print(td.kernel()[0][0] * 27.2114) - - td = TDHF(mf) - td.verbose = 5 - td.nstates = 5 - print(td.kernel()[0][0] * 27.2114) - - cell.spin = 2 - mf = scf.KUHF(cell, cell.make_kpts([2,1,1])).set(exxdiv=None) - mf.run() - - td = TDA(mf) - td.verbose = 5 - td.nstates = 5 - print(td.kernel()[0][0] * 27.2114) - - td = TDHF(mf) - td.verbose = 5 - td.nstates = 5 - print(td.kernel()[0][0] * 27.2114) diff --git a/pyscf/pbc/tdscf/kuks.py b/pyscf/pbc/tdscf/kuks.py index 1f76d12c8f..aee05f0f61 100644 --- a/pyscf/pbc/tdscf/kuks.py +++ b/pyscf/pbc/tdscf/kuks.py @@ -39,59 +39,3 @@ def kernel(self, x0=None): dft.kuks.KUKS.TDA = lib.class_as_method(KTDA) dft.kuks.KUKS.TDHF = None dft.kuks.KUKS.TDDFT = lib.class_as_method(TDDFT) - - -if __name__ == '__main__': - from pyscf.pbc import gto - from pyscf.pbc import dft - from pyscf.pbc import df - cell = gto.Cell() - cell.unit = 'B' - cell.atom = ''' - C 0. 0. 0. - C 1.68506879 1.68506879 1.68506879 - ''' - cell.a = ''' - 0. 3.37013758 3.37013758 - 3.37013758 0. 3.37013758 - 3.37013758 3.37013758 0. - ''' - - cell.basis = 'gth-szv' - cell.pseudo = 'gth-pade' - cell.mesh = [37]*3 - cell.build() - mf = dft.KUKS(cell, cell.make_kpts([2,1,1])).set(exxdiv=None, xc='b88,p86') - #mf.with_df = df.MDF(cell, cell.make_kpts([2,1,1])) - #mf.with_df.auxbasis = 'weigend' - #mf.with_df._cderi = 'eri3d-mdf.h5' - #mf.with_df.build(with_j3c=False) - mf.run() - - td = TDDFT(mf) - td.verbose = 5 - td.nstates = 5 - print(td.kernel()[0][0] * 27.2114) - - mf.xc = 'lda,vwn' - mf.run() - td = TDA(mf) - td.verbose = 5 - td.nstates = 5 - print(td.kernel()[0][0] * 27.2114) - - cell.spin = 2 - mf = dft.KUKS(cell, cell.make_kpts([2,1,1])).set(exxdiv=None, xc='b88,p86') - mf.run() - - td = TDDFT(mf) - td.verbose = 5 - td.nstates = 5 - print(td.kernel()[0][0] * 27.2114) - - mf.xc = 'lda,vwn' - mf.run() - td = TDA(mf) - td.verbose = 5 - td.nstates = 5 - print(td.kernel()[0][0] * 27.2114) diff --git a/pyscf/pbc/tdscf/rhf.py b/pyscf/pbc/tdscf/rhf.py index 9a311849b7..1a787c6640 100644 --- a/pyscf/pbc/tdscf/rhf.py +++ b/pyscf/pbc/tdscf/rhf.py @@ -20,11 +20,117 @@ # J. Mol. Struct. THEOCHEM, 914, 3 # +import numpy as np from pyscf import lib from pyscf.tdscf import rhf from pyscf.pbc import scf from pyscf import __config__ +def get_ab(mf): + r'''A and B matrices for TDDFT response function. + + A[i,a,j,b] = \delta_{ab}\delta_{ij}(E_a - E_i) + (ia||bj) + B[i,a,j,b] = (ia||jb) + + Ref: Chem Phys Lett, 256, 454 + ''' + cell = mf.cell + mo_energy = scf.addons.mo_energy_with_exxdiv_none(mf) + mo = np.asarray(mf.mo_coeff) + mo_occ = np.asarray(mf.mo_occ) + kpt = mf.kpt + nao, nmo = mo.shape + nocc = np.count_nonzero(mo_occ==2) + nvir = nmo - nocc + orbo = mo[:,:nocc] + orbv = mo[:,nocc:] + + e_ia = mo_energy[nocc:] - mo_energy[:nocc,None] + a = np.diag(e_ia.ravel()).reshape(nocc,nvir,nocc,nvir).astype(mo.dtype) + b = np.zeros_like(a) + + def add_hf_(a, b, hyb=1): + eri = mf.with_df.ao2mo([orbo,mo,mo,mo], kpt, compact=False) + eri = eri.reshape(nocc,nmo,nmo,nmo) + a += np.einsum('iabj->iajb', eri[:nocc,nocc:,nocc:,:nocc]) * 2 + a -= np.einsum('ijba->iajb', eri[:nocc,:nocc,nocc:,nocc:]) * hyb + b += np.einsum('iajb->iajb', eri[:nocc,nocc:,:nocc,nocc:]) * 2 + b -= np.einsum('ibja->iajb', eri[:nocc,nocc:,:nocc,nocc:]) * hyb + + if isinstance(mf, scf.hf.KohnShamDFT): + ni = mf._numint + omega, alpha, hyb = ni.rsh_and_hybrid_coeff(mf.xc, cell.spin) + + add_hf_(a, b, hyb) + if omega != 0: # For RSH + raise NotImplementedError + + xctype = ni._xc_type(mf.xc) + dm0 = mf.make_rdm1(mo, mo_occ) + make_rho = ni._gen_rho_evaluator(cell, dm0, hermi=1, with_lapl=False)[0] + mem_now = lib.current_memory()[0] + max_memory = max(2000, mf.max_memory*.8-mem_now) + + if xctype == 'LDA': + ao_deriv = 0 + for ao, _, mask, weight, coords \ + in ni.block_loop(cell, mf.grids, nao, ao_deriv, kpt, None, max_memory): + rho = make_rho(0, ao, mask, xctype) + fxc = ni.eval_xc_eff(mf.xc, rho, deriv=2, xctype=xctype)[2] + wfxc = fxc[0,0] * weight + + rho_o = lib.einsum('rp,pi->ri', ao, orbo) + rho_v = lib.einsum('rp,pi->ri', ao, orbv) + rho_ov = np.einsum('ri,ra->ria', rho_o, rho_v) + w_ov = np.einsum('ria,r->ria', rho_ov, wfxc) * 2 + rho_vo = rho_ov.conj() + a += lib.einsum('ria,rjb->iajb', w_ov, rho_vo) + b += lib.einsum('ria,rjb->iajb', w_ov, rho_ov) + + elif xctype == 'GGA': + ao_deriv = 1 + for ao, _, mask, weight, coords \ + in ni.block_loop(cell, mf.grids, nao, ao_deriv, kpt, None, max_memory): + rho = make_rho(0, ao, mask, xctype) + fxc = ni.eval_xc_eff(mf.xc, rho, deriv=2, xctype=xctype)[2] + wfxc = fxc * weight + rho_o = lib.einsum('xrp,pi->xri', ao, orbo) + rho_v = lib.einsum('xrp,pi->xri', ao, orbv) + rho_ov = np.einsum('xri,ra->xria', rho_o, rho_v[0]) + rho_ov[1:4] += np.einsum('ri,xra->xria', rho_o[0], rho_v[1:4]) + w_ov = np.einsum('xyr,xria->yria', wfxc, rho_ov) * 2 + rho_vo = rho_ov.conj() + a += lib.einsum('xria,xrjb->iajb', w_ov, rho_vo) + b += lib.einsum('xria,xrjb->iajb', w_ov, rho_ov) + + elif xctype == 'HF': + pass + + elif xctype == 'NLC': + raise NotImplementedError('NLC') + + elif xctype == 'MGGA': + ao_deriv = 1 + for ao, _, mask, weight, coords \ + in ni.block_loop(cell, mf.grids, nao, ao_deriv, kpt, None, max_memory): + rho = make_rho(0, ao, mask, xctype) + fxc = ni.eval_xc_eff(mf.xc, rho, deriv=2, xctype=xctype)[2] + wfxc = fxc * weight + rho_o = lib.einsum('xrp,pi->xri', ao, orbo) + rho_v = lib.einsum('xrp,pi->xri', ao, orbv) + rho_ov = np.einsum('xri,ra->xria', rho_o, rho_v[0]) + rho_ov[1:4] += np.einsum('ri,xra->xria', rho_o[0], rho_v[1:4]) + tau_ov = np.einsum('xri,xra->ria', rho_o[1:4], rho_v[1:4]) * .5 + rho_ov = np.vstack([rho_ov, tau_ov[np.newaxis]]) + w_ov = np.einsum('xyr,xria->yria', wfxc, rho_ov) * 2 + rho_vo = rho_ov.conj() + a += lib.einsum('xria,xrjb->iajb', w_ov, rho_vo) + b += lib.einsum('xria,xrjb->iajb', w_ov, rho_ov) + else: + add_hf_(a, b) + + return a, b + class TDBase(rhf.TDBase): _keys = {'cell'} @@ -33,7 +139,8 @@ def __init__(self, mf): self.cell = mf.cell def get_ab(self, mf=None): - raise NotImplementedError + if mf is None: mf = self._scf + return get_ab(mf) def nuc_grad_method(self): raise NotImplementedError @@ -57,7 +164,8 @@ class TDA(TDBase): kernel = rhf.TDA.kernel _gen_vind = rhf.TDA.gen_vind - def gen_vind(self, mf): + def gen_vind(self, mf=None): + if mf is None: mf = self._scf moe = scf.addons.mo_energy_with_exxdiv_none(mf) with lib.temporary_env(mf, mo_energy=moe): vind, hdiag = self._gen_vind(mf) @@ -69,7 +177,7 @@ def vindp(x): CIS = TDA -class TDHF(TDA): +class TDHF(TDBase): init_guess = rhf.TDHF.init_guess kernel = rhf.TDHF.kernel diff --git a/pyscf/pbc/tdscf/rks.py b/pyscf/pbc/tdscf/rks.py index 55da0c0911..348bc2e780 100644 --- a/pyscf/pbc/tdscf/rks.py +++ b/pyscf/pbc/tdscf/rks.py @@ -24,7 +24,8 @@ RPA = TDRKS = TDDFT = rhf.TDHF -class CasidaTDDFT(rhf.TDA): +class CasidaTDDFT(TDDFT): + init_guess = rhf.TDA.init_guess _gen_vind = rks.TDDFTNoHybrid.gen_vind gen_vind = rhf.TDA.gen_vind kernel = rks.TDDFTNoHybrid.kernel diff --git a/pyscf/pbc/tdscf/test/test_krhf.py b/pyscf/pbc/tdscf/test/test_krhf.py index eef979ed31..dbaf765e56 100644 --- a/pyscf/pbc/tdscf/test/test_krhf.py +++ b/pyscf/pbc/tdscf/test/test_krhf.py @@ -22,12 +22,22 @@ from pyscf.pbc.cc.eom_kccsd_rhf import EOMEESinglet from pyscf.data.nist import HARTREE2EV as unitev +def diagonalize(a, b, nroots=4): + nkpts, nocc, nvir = a.shape[:3] + a = a.reshape(nkpts*nocc*nvir, -1) + b = b.reshape(nkpts*nocc*nvir, -1) + h = np.block([[a , b ], + [-b.conj(),-a.conj()]]) + e = np.linalg.eigvals(np.asarray(h)) + lowest_e = np.sort(e[e.real > 0].real) + lowest_e = lowest_e[lowest_e > 1e-3][:nroots] + return lowest_e class Diamond(unittest.TestCase): @classmethod def setUpClass(cls): cell = gto.Cell() - cell.verbose = 4 + cell.verbose = 5 cell.output = '/dev/null' cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' cell.a = ''' @@ -52,11 +62,13 @@ def tearDownClass(cls): cls.cell.stdout.close() del cls.cell, cls.mf - def kernel(self, TD, ref, **kwargs): - td = TD(self.mf).set(kshift_lst=np.arange(len(self.mf.kpts)), + def kernel(self, TD, ref, kshift_lst, **kwargs): + td = TD(self.mf).set(kshift_lst=kshift_lst, nstates=self.nstates, **kwargs).run() for kshift,e in enumerate(td.e): - self.assertAlmostEqual(abs(e[:self.nstates_test] * unitev - ref[kshift]).max(), 0, 4) + n = len(ref[kshift]) + self.assertAlmostEqual(abs(e[:n] * unitev - ref[kshift]).max(), 0, 4) + return td def test_tda_singlet_eomccs(self): ''' Brute-force solution to the KTDA equation. Compared to the brute-force @@ -80,22 +92,52 @@ def test_tda_singlet_eomccs(self): def test_tda_singlet(self): ref = [[10.9573977036], [11.0418311085]] - self.kernel(tdscf.KTDA, ref) + td = self.kernel(tdscf.KTDA, ref, np.arange(2)) + a0, _ = td.get_ab(kshift=0) + nk, no, nv = a0.shape[:3] + a0 = a0.reshape(nk*no*nv,-1) + eref0 = np.linalg.eigvalsh(a0)[:4] + a1, _ = td.get_ab(kshift=1) + nk, no, nv = a1.shape[:3] + a1 = a1.reshape(nk*no*nv,-1) + eref1 = np.linalg.eigvalsh(a1)[:4] + self.assertAlmostEqual(abs(td.e[0][:2] - eref0[:2]).max(), 0, 5) + self.assertAlmostEqual(abs(td.e[1][:2] - eref1[:2]).max(), 0, 5) + + vind, hdiag = td.gen_vind(td._scf, kshift=0) + z = a0[:1] + self.assertAlmostEqual(abs(vind(z) - a0.dot(z[0])).max(), 0, 10) + vind, hdiag = td.gen_vind(td._scf, kshift=1) + self.assertAlmostEqual(abs(vind(z) - a1.dot(z[0])).max(), 0, 10) def test_tda_triplet(self): ref = [[6.4440137833], [7.4264899075]] - self.kernel(tdscf.KTDA, ref, singlet=False) + self.kernel(tdscf.KTDA, ref, np.arange(2), singlet=False) def test_tdhf_singlet(self): - ref = [[10.7665673889], - [10.8485048947]] - self.kernel(tdscf.KTDHF, ref) + # The first state can be obtained by solving nstates=9 + #ref = [[10.60761946, 10.76654619, 10.76654619]] + ref = [[10.76654619, 10.76654619]] + td = self.kernel(tdscf.KTDHF, ref, [0]) + a0, b0 = td.get_ab(kshift=0) + eref0 = diagonalize(a0, b0) + self.assertAlmostEqual(abs(td.e[0][:2] - eref0[1:3]).max(), 0, 5) + + nk, no, nv = a0.shape[:3] + nkov = nk * no * nv + a0 = a0.reshape(nkov,-1) + b0 = b0.reshape(nkov,-1) + z = np.hstack([a0[:1], a0[1:2]]) + + h = np.block([[ a0 , b0 ], + [-b0.conj(),-a0.conj()]]) + vind, hdiag = td.gen_vind(td._scf, kshift=0) + self.assertAlmostEqual(abs(vind(z) - h.dot(z[0])).max(), 0, 10) def test_tdhf_triplet(self): - ref = [[5.9794378466], - [6.1703909932]] - self.kernel(tdscf.KTDHF, ref, singlet=False) + ref = [[5.9794378466]] + self.kernel(tdscf.KTDHF, ref, [0], singlet=False) class WaterBigBox(unittest.TestCase): diff --git a/pyscf/pbc/tdscf/test/test_krks.py b/pyscf/pbc/tdscf/test/test_krks.py index c99315a043..e16737f9c8 100644 --- a/pyscf/pbc/tdscf/test/test_krks.py +++ b/pyscf/pbc/tdscf/test/test_krks.py @@ -23,6 +23,17 @@ from pyscf.data.nist import HARTREE2EV as unitev +def diagonalize(a, b, nroots=4): + nkpts, nocc, nvir = a.shape[:3] + a = a.reshape(nkpts*nocc*nvir, -1) + b = b.reshape(nkpts*nocc*nvir, -1) + h = np.block([[a , b ], + [-b.conj(),-a.conj()]]) + e = np.linalg.eigvals(np.asarray(h)) + lowest_e = np.sort(e[e.real > 0].real)[:nroots] + lowest_e = lowest_e[lowest_e > 1e-3] + return lowest_e + class DiamondPBE(unittest.TestCase): @classmethod def setUpClass(cls): @@ -43,28 +54,37 @@ def setUpClass(cls): kpts = cell.make_kpts((2,1,1)) xc = 'pbe' - mf = scf.KRKS(cell, kpts=kpts).set(xc=xc).rs_density_fit(auxbasis='weigend').run() + mf = scf.KRKS(cell, kpts=kpts).set(xc=xc).rs_density_fit(auxbasis='weigend') + mf.with_df._j_only = False + mf.with_df.build() + mf.run() cls.cell = cell cls.mf = mf cls.nstates = 4 # make sure first `nstates_test` states are converged cls.nstates_test = 2 + @classmethod def tearDownClass(cls): cls.cell.stdout.close() del cls.cell, cls.mf def kernel(self, TD, ref, **kwargs): - td = getattr(self.mf, TD)().set(kshift_lst=np.arange(len(self.mf.kpts)), + td = getattr(self.mf, TD)().set(kshift_lst=np.arange(len(ref)), nstates=self.nstates, **kwargs).run() for kshift,e in enumerate(td.e): - self.assertAlmostEqual(abs(e[:self.nstates_test] * unitev - ref[kshift]).max(), 0, 4) + self.assertAlmostEqual(abs(e[:self.nstates_test] * unitev - ref[kshift]).max(), 0, 5) + return td def test_tda_singlet(self): - ref = [[7.7172854747, 7.7173219160], + ref = [[7.7172857896, 7.7173222336], [8.3749594280, 8.3749980463]] - self.kernel('TDA', ref) + td = self.kernel('TDA', ref) + a0, _ = td.get_ab(kshift=0) + nk, no, nv = a0.shape[:3] + eref = np.linalg.eigvalsh(a0.reshape(nk*no*nv,-1))[:4] + self.assertAlmostEqual(abs(td.e[0][:2] - eref[:2]).max(), 0, 7) def test_tda_triplet(self): ref = [[5.7465112548, 5.7465526327], @@ -72,13 +92,14 @@ def test_tda_triplet(self): self.kernel('TDA', ref, singlet=False) def test_tdhf_singlet(self): - ref = [[7.5824302393, 7.5824675688], - [8.3438648420, 8.3439036271]] - self.kernel('TDDFT', ref) + ref = [[7.58243026, 7.58246786]] + td = self.kernel('TDDFT', ref) + a, b = td.get_ab(kshift=0) + eref = diagonalize(a, b) + self.assertAlmostEqual(abs(td.e[0][:2] - eref[:2]).max(), 0, 7) def test_tdhf_triplet(self): - ref = [[5.5659966435, 5.5660393021], - [6.7992845776, 6.7993290046]] + ref = [[5.56599665, 5.56603980]] self.kernel('TDDFT', ref, singlet=False) @@ -115,15 +136,20 @@ def tearDownClass(cls): del cls.cell, cls.mf def kernel(self, TD, ref, **kwargs): - td = getattr(self.mf, TD)().set(kshift_lst=np.arange(len(self.mf.kpts)), + td = getattr(self.mf, TD)().set(kshift_lst=np.arange(len(ref)), nstates=self.nstates, **kwargs).run() for kshift,e in enumerate(td.e): - self.assertAlmostEqual(abs(e[:self.nstates_test] * unitev - ref[kshift]).max(), 0, 4) + self.assertAlmostEqual(abs(e[:self.nstates_test] * unitev - ref[kshift]).max(), 0, 5) + return td def test_tda_singlet(self): ref = [[9.3936718451, 9.4874866060], [10.0697605303, 10.0697862958]] - self.kernel('TDA', ref) + td = self.kernel('TDA', ref) + a0, _ = td.get_ab(kshift=0) + nk, no, nv = a0.shape[:3] + eref = np.linalg.eigvalsh(a0.reshape(nk*no*nv,-1))[:4] + self.assertAlmostEqual(abs(td.e[0][:2] - eref[:2]).max(), 0, 8) def test_tda_triplet(self): ref = [[6.6703797643, 6.6704110631], @@ -131,13 +157,14 @@ def test_tda_triplet(self): self.kernel('TDA', ref, singlet=False) def test_tdhf_singlet(self): - ref = [[9.2519208050, 9.3208025447], - [9.9609751875, 9.9610015227]] - self.kernel('TDDFT', ref) + ref = [[9.25192096, 9.32080304]] + td = self.kernel('TDDFT', ref, conv_tol=1e-6) + a, b = td.get_ab(kshift=0) + eref = diagonalize(a, b) + self.assertAlmostEqual(abs(td.e[0][:2] - eref[:2]).max(), 0, 8) def test_tdhf_triplet(self): - ref = [[6.3282716764, 6.3283051217], - [7.0656766298, 7.0657111705]] + ref = [[6.3282716764, 6.3283051217]] self.kernel('TDDFT', ref, singlet=False) diff --git a/pyscf/pbc/tdscf/test/test_kuhf.py b/pyscf/pbc/tdscf/test/test_kuhf.py index 903090c8d2..4f306ee6cc 100644 --- a/pyscf/pbc/tdscf/test/test_kuhf.py +++ b/pyscf/pbc/tdscf/test/test_kuhf.py @@ -53,8 +53,8 @@ def tearDownClass(cls): cls.cell.stdout.close() del cls.cell, cls.mf - def kernel(self, TD, ref, **kwargs): - td = TD(self.mf).set(kshift_lst=np.arange(len(self.mf.kpts)), + def kernel(self, TD, ref, kshift_lst, **kwargs): + td = TD(self.mf).set(kshift_lst=kshift_lst, nstates=self.nstates, **kwargs).run() for kshift,e in enumerate(td.e): self.assertAlmostEqual(abs(e[:self.nstates_test] * unitev - ref[kshift]).max(), 0, 4) @@ -63,13 +63,12 @@ def test_tda(self): # same as lowest roots in Diamond->test_tda_singlet/triplet in test_krhf.py ref = [[6.4440137833, 7.5317890777], [7.4264899075, 7.6381352853]] - self.kernel(tdscf.KTDA, ref) + self.kernel(tdscf.KTDA, ref, np.arange(len(ref))) def test_tdhf(self): # same as lowest roots in Diamond->test_tdhf_singlet/triplet in test_krhf.py - ref = [[5.9794378466, 5.9794378466], - [6.1703909932, 6.1703909932]] - self.kernel(tdscf.KTDHF, ref) + ref = [[5.9794378466, 5.9794378466]] + self.kernel(tdscf.KTDHF, ref, [0]) class WaterBigBox(unittest.TestCase): @@ -116,8 +115,7 @@ def tearDownClass(cls): del cls.mol, cls.molmf def kernel(self, TD, MOLTD, **kwargs): - td = TD(self.mf).set(kshift_lst=np.arange(len(self.mf.kpts)), - nstates=self.nstates, **kwargs).run() + td = TD(self.mf).set(nstates=self.nstates, **kwargs).run() moltd = MOLTD(self.molmf).set(nstates=self.nstates, **kwargs).run() ref = moltd.e for kshift,e in enumerate(td.e): diff --git a/pyscf/pbc/tdscf/test/test_kuks.py b/pyscf/pbc/tdscf/test/test_kuks.py index aa05bab624..15a93792c3 100644 --- a/pyscf/pbc/tdscf/test/test_kuks.py +++ b/pyscf/pbc/tdscf/test/test_kuks.py @@ -58,7 +58,7 @@ def tearDownClass(cls): del cls.cell, cls.mf def kernel(self, TD, ref, **kwargs): - td = getattr(self.mf, TD)().set(kshift_lst=np.arange(len(self.mf.kpts)), + td = getattr(self.mf, TD)().set(kshift_lst=np.arange(len(ref)), nstates=self.nstates, **kwargs).run() for kshift,e in enumerate(td.e): self.assertAlmostEqual(abs(e[:self.nstates_test] * unitev - ref[kshift]).max(), 0, 4) @@ -71,8 +71,7 @@ def test_tda(self): def test_tdhf(self): # same as lowest roots in DiamondPBE->test_tdhf_singlet/triplet in test_krks.py - ref = [[5.5659966435, 5.5660393021], - [6.7992845776, 6.7993290046]] + ref = [[5.56602, 5.56602]] self.kernel('TDDFT', ref, singlet=False) @@ -111,7 +110,7 @@ def tearDownClass(cls): del cls.cell, cls.mf def kernel(self, TD, ref, **kwargs): - td = getattr(self.mf, TD)().set(kshift_lst=np.arange(len(self.mf.kpts)), + td = getattr(self.mf, TD)().set(kshift_lst=np.arange(len(ref)), nstates=self.nstates, **kwargs).run() for kshift,e in enumerate(td.e): self.assertAlmostEqual(abs(e[:self.nstates_test] * unitev - ref[kshift]).max(), 0, 4) @@ -124,8 +123,7 @@ def test_tda(self): def test_tdhf(self): # same as lowest roots in DiamondPBE0->test_tdhf_singlet/triplet in test_krks.py - ref = [[6.3282716764, 6.3283051217], - [7.0656766298, 7.0657111705]] + ref = [[6.3282716764, 6.3283051217]] self.kernel('TDDFT', ref, singlet=False) diff --git a/pyscf/pbc/tdscf/test/test_rhf.py b/pyscf/pbc/tdscf/test/test_rhf.py index 8bac9829a2..32d385c388 100644 --- a/pyscf/pbc/tdscf/test/test_rhf.py +++ b/pyscf/pbc/tdscf/test/test_rhf.py @@ -55,6 +55,7 @@ def tearDownClass(cls): def kernel(self, TD, ref, **kwargs): td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + return td def test_tda_singlet(self): ref = [9.6425852906, 9.6425852906] @@ -66,7 +67,15 @@ def test_tda_triplet(self): def test_tdhf_singlet(self): ref = [9.2573219105, 9.2573219106] - self.kernel('TDHF', ref) + td = self.kernel('TDHF', ref) + a, b = td.get_ab() + no, nv = a.shape[:2] + a = a.reshape(no*nv,-1) + b = b.reshape(no*nv,-1) + h = np.block([[a, b], [-b.conj(), -a.conj()]]) + eref = np.sort(np.linalg.eigvals(h).real) + eref = eref[eref > 1e-3] + self.assertAlmostEqual(abs(td.e[:2] - eref[:2]).max(), 0, 5) def test_tdhf_triplet(self): ref = [3.0396052214, 3.0396052214] @@ -108,6 +117,7 @@ def tearDownClass(cls): def kernel(self, TD, ref, **kwargs): td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + return td def test_tda_singlet(self): ref = [12.7166510188, 13.5460934688] @@ -119,7 +129,15 @@ def test_tda_triplet(self): def test_tdhf_singlet(self): ref = [12.6104811733, 13.4160717812] - self.kernel('TDHF', ref) + td = self.kernel('TDHF', ref) + a, b = td.get_ab() + no, nv = a.shape[:2] + a = a.reshape(no*nv,-1) + b = b.reshape(no*nv,-1) + h = np.block([[a, b], [-b.conj(), -a.conj()]]) + eref = np.sort(np.linalg.eigvals(h).real) + eref = eref[eref > 1e-3] + self.assertAlmostEqual(abs(td.e[:2] - eref[:2]).max(), 0, 5) def test_tdhf_triplet(self): ref = [3.8940277713, 7.9448161493] diff --git a/pyscf/pbc/tdscf/test/test_rks.py b/pyscf/pbc/tdscf/test/test_rks.py index 6ebbbc5343..4be1045424 100644 --- a/pyscf/pbc/tdscf/test/test_rks.py +++ b/pyscf/pbc/tdscf/test/test_rks.py @@ -22,6 +22,17 @@ from pyscf.data.nist import HARTREE2EV as unitev +def diagonalize(a, b, nroots=4): + nocc, nvir = a.shape[:2] + a = a.reshape(nocc*nvir, -1) + b = b.reshape(nocc*nvir, -1) + h = np.block([[a , b ], + [-b.conj(),-a.conj()]]) + e = np.linalg.eigvals(np.asarray(h)) + lowest_e = np.sort(e[e.real > 0].real) + lowest_e = lowest_e[lowest_e > 1e-3][:nroots] + return lowest_e + class DiamondPBE(unittest.TestCase): @classmethod def setUpClass(cls): @@ -40,7 +51,7 @@ def setUpClass(cls): cell.precision = 1e-10 cell.build() - xc = 'pbe' + xc = 'm06' mf = scf.RKS(cell).set(xc=xc).rs_density_fit(auxbasis='weigend').run() cls.cell = cell @@ -55,22 +66,30 @@ def tearDownClass(cls): def kernel(self, TD, ref, **kwargs): td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() - self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 5) + return td def test_tda_singlet(self): - ref = [9.26752544, 9.26752544] - self.kernel('TDA', ref) + ref = [14.68587442, 14.68589929] + td = self.kernel('TDA', ref) + a, b = td.get_ab() + no, nv = a.shape[:2] + eref = np.linalg.eigvalsh(a.reshape(no*nv,-1)) + self.assertAlmostEqual(abs(td.e[:2] - eref[:2]).max(), 0, 8) def test_tda_triplet(self): - ref = [4.79532598, 4.79532598] + ref = [11.10049832, 11.59365532] self.kernel('TDA', ref, singlet=False) def test_tddft_singlet(self): - ref = [8.8773512896, 8.8773512896] - self.kernel('TDDFT', ref) + ref = [14.42819773, 14.42822009] + td = self.kernel('TDDFT', ref) + a, b = td.get_ab() + eref = diagonalize(a, b) + self.assertAlmostEqual(abs(td.e[:2] - eref[:2]).max(), 0, 8) def test_tddft_triplet(self): - ref = [4.7694490556, 4.7694490556] + ref = [ 9.09496456, 11.53650896] self.kernel('TDDFT', ref, singlet=False) @@ -109,7 +128,7 @@ def tearDownClass(cls): def kernel(self, TD, ref, **kwargs): td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() - self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 5) def test_tda_singlet(self): ref = [11.9664870288, 12.7605699008] @@ -223,22 +242,30 @@ def tearDownClass(cls): def kernel(self, TD, ref, **kwargs): td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() - self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 5) + return td def test_tda_singlet(self): - ref = [9.61243337, 9.61243337] - self.kernel('TDA', ref) + ref = [9.62238067, 9.62238067] + td = self.kernel('TDA', ref) + a, b = td.get_ab() + no, nv = a.shape[:2] + eref = np.linalg.eigvalsh(a.reshape(no*nv,-1)) + self.assertAlmostEqual(abs(td.e[:2] - eref[:2]).max(), 0, 8) def test_tda_triplet(self): - ref = [5.13034084, 5.13034084] + ref = [5.39995144, 5.39995144] self.kernel('TDA', ref, singlet=False) def test_tddft_singlet(self): - ref = [9.2585491075, 9.2585491075] - self.kernel('TDDFT', ref) + ref = [9.26011401, 9.26011401] + td = self.kernel('TDDFT', ref) + a, b = td.get_ab() + eref = diagonalize(a, b) + self.assertAlmostEqual(abs(td.e[:2] - eref[:2]).max(), 0, 8) def test_tddft_triplet(self): - ref = [4.5570089807, 4.5570089807] + ref = [4.68905023, 4.81439580] self.kernel('TDDFT', ref, singlet=False) diff --git a/pyscf/pbc/tdscf/test/test_uks.py b/pyscf/pbc/tdscf/test/test_uks.py index c483813cfc..6c232bd975 100644 --- a/pyscf/pbc/tdscf/test/test_uks.py +++ b/pyscf/pbc/tdscf/test/test_uks.py @@ -22,7 +22,28 @@ from pyscf.data.nist import HARTREE2EV as unitev -class DiamondPBE(unittest.TestCase): +def diagonalize(a, b, nroots=4): + a_aa, a_ab, a_bb = a + b_aa, b_ab, b_bb = b + nocc_a, nvir_a, nocc_b, nvir_b = a_ab.shape + a_aa = a_aa.reshape((nocc_a*nvir_a,nocc_a*nvir_a)) + a_ab = a_ab.reshape((nocc_a*nvir_a,nocc_b*nvir_b)) + a_bb = a_bb.reshape((nocc_b*nvir_b,nocc_b*nvir_b)) + b_aa = b_aa.reshape((nocc_a*nvir_a,nocc_a*nvir_a)) + b_ab = b_ab.reshape((nocc_a*nvir_a,nocc_b*nvir_b)) + b_bb = b_bb.reshape((nocc_b*nvir_b,nocc_b*nvir_b)) + a = np.block([[ a_aa , a_ab], + [ a_ab.T, a_bb]]) + b = np.block([[ b_aa , b_ab], + [ b_ab.T, b_bb]]) + abba = np.block([[a , b ], + [-b.conj(),-a.conj()]]) + e = np.linalg.eig(abba)[0] + lowest_e = np.sort(e[e.real > 0].real) + lowest_e = lowest_e[lowest_e > 1e-3][:nroots] + return lowest_e + +class DiamondM06(unittest.TestCase): ''' Reproduce RKS-TDSCF results ''' @classmethod @@ -42,7 +63,7 @@ def setUpClass(cls): cell.precision = 1e-10 cell.build() - xc = 'pbe' + xc = 'm06' mf = scf.UKS(cell).set(xc=xc).rs_density_fit(auxbasis='weigend').run() cls.cell = cell cls.mf = mf @@ -56,15 +77,29 @@ def tearDownClass(cls): def kernel(self, TD, ref, **kwargs): td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() - self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 5) + return td def test_tda(self): - ref = [4.77090949, 4.77090949] - self.kernel('TDA', ref) + ref = [11.09731427, 11.57079413] + td = self.kernel('TDA', ref) + a, b = td.get_ab() + a_aa, a_ab, a_bb = a + nocc_a, nvir_a, nocc_b, nvir_b = a_ab.shape + a_aa = a_aa.reshape((nocc_a*nvir_a,nocc_a*nvir_a)) + a_ab = a_ab.reshape((nocc_a*nvir_a,nocc_b*nvir_b)) + a_bb = a_bb.reshape((nocc_b*nvir_b,nocc_b*nvir_b)) + a = np.block([[ a_aa , a_ab], + [ a_ab.T, a_bb]]) + eref = np.linalg.eigvalsh(a) + self.assertAlmostEqual(abs(td.e[:4] - eref[:4]).max(), 0, 8) def test_tdhf(self): - ref = [4.7455726874, 4.7455726874] - self.kernel('TDDFT', ref) + ref = [9.09165361, 11.51362009] + td = self.kernel('TDDFT', ref) + a, b = td.get_ab() + eref = diagonalize(a, b) + self.assertAlmostEqual(abs(td.e[:4] - eref[:4]).max(), 0, 8) class WaterBigBoxPBE(unittest.TestCase): @@ -151,6 +186,7 @@ def setUpClass(cls): cls.nstates = 5 # make sure first `nstates_test` states are converged cls.nstates_test = 2 + @classmethod def tearDownClass(cls): cls.cell.stdout.close() @@ -158,15 +194,30 @@ def tearDownClass(cls): def kernel(self, TD, ref, **kwargs): td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() - self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 5) + return td def test_tda(self): - ref = [5.1069789, 5.10697989] - self.kernel('TDA', ref) + ref = [5.37745381, 5.37745449] + td = self.kernel('TDA', ref) + a, b = td.get_ab() + a_aa, a_ab, a_bb = a + nocc_a, nvir_a, nocc_b, nvir_b = a_ab.shape + a_aa = a_aa.reshape((nocc_a*nvir_a,nocc_a*nvir_a)) + a_ab = a_ab.reshape((nocc_a*nvir_a,nocc_b*nvir_b)) + a_bb = a_bb.reshape((nocc_b*nvir_b,nocc_b*nvir_b)) + a = np.block([[ a_aa , a_ab], + [ a_ab.T, a_bb]]) + eref = np.linalg.eigvalsh(a) + self.assertAlmostEqual(abs(td.e[:2] - eref[:2]).max(), 0, 8) def test_tdhf(self): - ref = [4.5331029147, 4.5331029147] - self.kernel('TDDFT', ref) + # nstates=6 is required to derive the lowest state + ref = [4.6851639, 4.79043398, 4.79043398] + td = self.kernel('TDDFT', ref[1:3]) + a, b = td.get_ab() + eref = diagonalize(a, b) + self.assertAlmostEqual(abs(td.e[:2] - eref[1:3]).max(), 0, 8) class WaterBigBoxPBE0(unittest.TestCase): diff --git a/pyscf/pbc/tdscf/uhf.py b/pyscf/pbc/tdscf/uhf.py index 42e7599acc..ec8b86977e 100644 --- a/pyscf/pbc/tdscf/uhf.py +++ b/pyscf/pbc/tdscf/uhf.py @@ -16,14 +16,212 @@ # Author: Qiming Sun # +import numpy as np from pyscf import lib from pyscf.tdscf import uhf from pyscf.pbc.tdscf import rhf as td_rhf from pyscf.pbc.tdscf.rhf import TDBase +def get_ab(mf): + r'''A and B matrices for TDDFT response function. + + A[i,a,j,b] = \delta_{ab}\delta_{ij}(E_a - E_i) + (ia||bj) + B[i,a,j,b] = (ia||jb) + + Spin symmetry is considered in the returned A, B lists. List A has three + items: (A_aaaa, A_aabb, A_bbbb). A_bbaa = A_aabb.transpose(2,3,0,1). + B has three items: (B_aaaa, B_aabb, B_bbbb). + B_bbaa = B_aabb.transpose(2,3,0,1). + ''' + cell = mf.cell + nao = cell.nao_nr() + mo_energy = scf.addons.mo_energy_with_exxdiv_none(mf) + mo = np.asarray(mf.mo_coeff) + mo_occ = np.asarray(mf.mo_occ) + kpt = mf.kpt + + occidx_a = np.where(mo_occ[0]==1)[0] + viridx_a = np.where(mo_occ[0]==0)[0] + occidx_b = np.where(mo_occ[1]==1)[0] + viridx_b = np.where(mo_occ[1]==0)[0] + orbo_a = mo[0][:,occidx_a] + orbv_a = mo[0][:,viridx_a] + orbo_b = mo[1][:,occidx_b] + orbv_b = mo[1][:,viridx_b] + nocc_a = orbo_a.shape[1] + nvir_a = orbv_a.shape[1] + nocc_b = orbo_b.shape[1] + nvir_b = orbv_b.shape[1] + mo_a = np.hstack((orbo_a,orbv_a)) + mo_b = np.hstack((orbo_b,orbv_b)) + nmo_a = nocc_a + nvir_a + nmo_b = nocc_b + nvir_b + + e_ia_a = (mo_energy[0][viridx_a,None] - mo_energy[0][occidx_a]).T + e_ia_b = (mo_energy[1][viridx_b,None] - mo_energy[1][occidx_b]).T + a_aa = np.diag(e_ia_a.ravel()).reshape(nocc_a,nvir_a,nocc_a,nvir_a) + a_bb = np.diag(e_ia_b.ravel()).reshape(nocc_b,nvir_b,nocc_b,nvir_b) + a_ab = np.zeros((nocc_a,nvir_a,nocc_b,nvir_b)) + b_aa = np.zeros_like(a_aa) + b_ab = np.zeros_like(a_ab) + b_bb = np.zeros_like(a_bb) + a = (a_aa, a_ab, a_bb) + b = (b_aa, b_ab, b_bb) + + def add_hf_(a, b, hyb=1): + eri_aa = mf.with_df.ao2mo([orbo_a,mo_a,mo_a,mo_a], kpt, compact=False) + eri_ab = mf.with_df.ao2mo([orbo_a,mo_a,mo_b,mo_b], kpt, compact=False) + eri_bb = mf.with_df.ao2mo([orbo_b,mo_b,mo_b,mo_b], kpt, compact=False) + eri_aa = eri_aa.reshape(nocc_a,nmo_a,nmo_a,nmo_a) + eri_ab = eri_ab.reshape(nocc_a,nmo_a,nmo_b,nmo_b) + eri_bb = eri_bb.reshape(nocc_b,nmo_b,nmo_b,nmo_b) + a_aa, a_ab, a_bb = a + b_aa, b_ab, b_bb = b + + a_aa += np.einsum('iabj->iajb', eri_aa[:nocc_a,nocc_a:,nocc_a:,:nocc_a]) + a_aa -= np.einsum('ijba->iajb', eri_aa[:nocc_a,:nocc_a,nocc_a:,nocc_a:]) * hyb + b_aa += np.einsum('iajb->iajb', eri_aa[:nocc_a,nocc_a:,:nocc_a,nocc_a:]) + b_aa -= np.einsum('jaib->iajb', eri_aa[:nocc_a,nocc_a:,:nocc_a,nocc_a:]) * hyb + + a_bb += np.einsum('iabj->iajb', eri_bb[:nocc_b,nocc_b:,nocc_b:,:nocc_b]) + a_bb -= np.einsum('ijba->iajb', eri_bb[:nocc_b,:nocc_b,nocc_b:,nocc_b:]) * hyb + b_bb += np.einsum('iajb->iajb', eri_bb[:nocc_b,nocc_b:,:nocc_b,nocc_b:]) + b_bb -= np.einsum('jaib->iajb', eri_bb[:nocc_b,nocc_b:,:nocc_b,nocc_b:]) * hyb + + a_ab += np.einsum('iabj->iajb', eri_ab[:nocc_a,nocc_a:,nocc_b:,:nocc_b]) + b_ab += np.einsum('iajb->iajb', eri_ab[:nocc_a,nocc_a:,:nocc_b,nocc_b:]) + + if isinstance(mf, scf.hf.KohnShamDFT): + ni = mf._numint + omega, alpha, hyb = ni.rsh_and_hybrid_coeff(mf.xc, cell.spin) + + add_hf_(a, b, hyb) + if omega != 0: # For RSH + raise NotImplementedError + + xctype = ni._xc_type(mf.xc) + dm0 = mf.make_rdm1(mo, mo_occ) + make_rho = ni._gen_rho_evaluator(cell, dm0, hermi=1, with_lapl=False)[0] + mem_now = lib.current_memory()[0] + max_memory = max(2000, mf.max_memory*.8-mem_now) + + if xctype == 'LDA': + ao_deriv = 0 + for ao, _, mask, weight, coords \ + in ni.block_loop(cell, mf.grids, nao, ao_deriv, kpt, None, max_memory): + rho0a = make_rho(0, ao, mask, xctype) + rho0b = make_rho(1, ao, mask, xctype) + rho = (rho0a, rho0b) + fxc = ni.eval_xc_eff(mf.xc, rho, deriv=2, xctype=xctype)[2] + wfxc = fxc[:,0,:,0] * weight + + rho_o_a = lib.einsum('rp,pi->ri', ao, orbo_a) + rho_v_a = lib.einsum('rp,pi->ri', ao, orbv_a) + rho_o_b = lib.einsum('rp,pi->ri', ao, orbo_b) + rho_v_b = lib.einsum('rp,pi->ri', ao, orbv_b) + rho_ov_a = np.einsum('ri,ra->ria', rho_o_a, rho_v_a) + rho_ov_b = np.einsum('ri,ra->ria', rho_o_b, rho_v_b) + rho_vo_a = rho_ov_a.conj() + rho_vo_b = rho_ov_b.conj() + w_ov_aa = np.einsum('ria,r->ria', rho_ov_a, wfxc[0,0]) + w_ov_ab = np.einsum('ria,r->ria', rho_ov_a, wfxc[0,1]) + w_ov_bb = np.einsum('ria,r->ria', rho_ov_b, wfxc[1,1]) + + a_aa += lib.einsum('ria,rjb->iajb', w_ov_aa, rho_vo_a) + b_aa += lib.einsum('ria,rjb->iajb', w_ov_aa, rho_ov_a) + + a_ab += lib.einsum('ria,rjb->iajb', w_ov_ab, rho_vo_b) + b_ab += lib.einsum('ria,rjb->iajb', w_ov_ab, rho_ov_b) + + a_bb += lib.einsum('ria,rjb->iajb', w_ov_bb, rho_vo_b) + b_bb += lib.einsum('ria,rjb->iajb', w_ov_bb, rho_ov_b) + + elif xctype == 'GGA': + ao_deriv = 1 + for ao, _, mask, weight, coords \ + in ni.block_loop(cell, mf.grids, nao, ao_deriv, kpt, None, max_memory): + rho0a = make_rho(0, ao, mask, xctype) + rho0b = make_rho(1, ao, mask, xctype) + rho = (rho0a, rho0b) + fxc = ni.eval_xc_eff(mf.xc, rho, deriv=2, xctype=xctype)[2] + wfxc = fxc * weight + rho_o_a = lib.einsum('xrp,pi->xri', ao, orbo_a) + rho_v_a = lib.einsum('xrp,pi->xri', ao, orbv_a) + rho_o_b = lib.einsum('xrp,pi->xri', ao, orbo_b) + rho_v_b = lib.einsum('xrp,pi->xri', ao, orbv_b) + rho_ov_a = np.einsum('xri,ra->xria', rho_o_a, rho_v_a[0]) + rho_ov_b = np.einsum('xri,ra->xria', rho_o_b, rho_v_b[0]) + rho_ov_a[1:4] += np.einsum('ri,xra->xria', rho_o_a[0], rho_v_a[1:4]) + rho_ov_b[1:4] += np.einsum('ri,xra->xria', rho_o_b[0], rho_v_b[1:4]) + rho_vo_a = rho_ov_a.conj() + rho_vo_b = rho_ov_b.conj() + w_ov_aa = np.einsum('xyr,xria->yria', wfxc[0,:,0], rho_ov_a) + w_ov_ab = np.einsum('xyr,xria->yria', wfxc[0,:,1], rho_ov_a) + w_ov_bb = np.einsum('xyr,xria->yria', wfxc[1,:,1], rho_ov_b) + + a_aa += lib.einsum('xria,xrjb->iajb', w_ov_aa, rho_vo_a) + b_aa += lib.einsum('xria,xrjb->iajb', w_ov_aa, rho_ov_a) + + a_ab += lib.einsum('xria,xrjb->iajb', w_ov_ab, rho_vo_b) + b_ab += lib.einsum('xria,xrjb->iajb', w_ov_ab, rho_ov_b) + + a_bb += lib.einsum('xria,xrjb->iajb', w_ov_bb, rho_vo_b) + b_bb += lib.einsum('xria,xrjb->iajb', w_ov_bb, rho_ov_b) + + elif xctype == 'HF': + pass + + elif xctype == 'NLC': + raise NotImplementedError('NLC') + + elif xctype == 'MGGA': + ao_deriv = 1 + for ao, _, mask, weight, coords \ + in ni.block_loop(cell, mf.grids, nao, ao_deriv, kpt, None, max_memory): + rho0a = make_rho(0, ao, mask, xctype) + rho0b = make_rho(1, ao, mask, xctype) + rho = (rho0a, rho0b) + fxc = ni.eval_xc_eff(mf.xc, rho, deriv=2, xctype=xctype)[2] + wfxc = fxc * weight + rho_oa = lib.einsum('xrp,pi->xri', ao, orbo_a) + rho_ob = lib.einsum('xrp,pi->xri', ao, orbo_b) + rho_va = lib.einsum('xrp,pi->xri', ao, orbv_a) + rho_vb = lib.einsum('xrp,pi->xri', ao, orbv_b) + rho_ov_a = np.einsum('xri,ra->xria', rho_oa, rho_va[0]) + rho_ov_b = np.einsum('xri,ra->xria', rho_ob, rho_vb[0]) + rho_ov_a[1:4] += np.einsum('ri,xra->xria', rho_oa[0], rho_va[1:4]) + rho_ov_b[1:4] += np.einsum('ri,xra->xria', rho_ob[0], rho_vb[1:4]) + tau_ov_a = np.einsum('xri,xra->ria', rho_oa[1:4], rho_va[1:4]) * .5 + tau_ov_b = np.einsum('xri,xra->ria', rho_ob[1:4], rho_vb[1:4]) * .5 + rho_ov_a = np.vstack([rho_ov_a, tau_ov_a[np.newaxis]]) + rho_ov_b = np.vstack([rho_ov_b, tau_ov_b[np.newaxis]]) + rho_vo_a = rho_ov_a.conj() + rho_vo_b = rho_ov_b.conj() + w_ov_aa = np.einsum('xyr,xria->yria', wfxc[0,:,0], rho_ov_a) + w_ov_ab = np.einsum('xyr,xria->yria', wfxc[0,:,1], rho_ov_a) + w_ov_bb = np.einsum('xyr,xria->yria', wfxc[1,:,1], rho_ov_b) + + a_aa += lib.einsum('xria,xrjb->iajb', w_ov_aa, rho_vo_a) + b_aa += lib.einsum('xria,xrjb->iajb', w_ov_aa, rho_ov_a) + + a_ab += lib.einsum('xria,xrjb->iajb', w_ov_ab, rho_vo_b) + b_ab += lib.einsum('xria,xrjb->iajb', w_ov_ab, rho_ov_b) + + a_bb += lib.einsum('xria,xrjb->iajb', w_ov_bb, rho_vo_b) + b_bb += lib.einsum('xria,xrjb->iajb', w_ov_bb, rho_ov_b) + + else: + add_hf_(a, b) + + return a, b + class TDA(TDBase): + def get_ab(self, mf=None): + if mf is None: mf = self._scf + return get_ab(mf) + singlet = None init_guess = uhf.TDA.init_guess @@ -34,7 +232,11 @@ class TDA(TDBase): CIS = TDA -class TDHF(TDA): +class TDHF(TDBase): + + def get_ab(self, mf=None): + if mf is None: mf = self._scf + return get_ab(mf) singlet = None diff --git a/pyscf/pbc/tdscf/uks.py b/pyscf/pbc/tdscf/uks.py index c8b1fc7e34..5f1099922e 100644 --- a/pyscf/pbc/tdscf/uks.py +++ b/pyscf/pbc/tdscf/uks.py @@ -25,7 +25,8 @@ RPA = TDUKS = TDDFT -class CasidaTDDFT(TDA): +class CasidaTDDFT(TDDFT): + init_guess = TDA.init_guess _gen_vind = uks.TDDFTNoHybrid.gen_vind gen_vind = TDA.gen_vind kernel = uks.TDDFTNoHybrid.kernel diff --git a/pyscf/tdscf/_lr_eig.py b/pyscf/tdscf/_lr_eig.py new file mode 100644 index 0000000000..98595a2526 --- /dev/null +++ b/pyscf/tdscf/_lr_eig.py @@ -0,0 +1,505 @@ +#!/usr/bin/env python +# Copyright 2024 The PySCF Developers. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import sys +import numpy as np +import scipy.linalg +from pyscf.lib.parameters import MAX_MEMORY +from pyscf.lib import logger +from pyscf.lib.linalg_helper import _sort_elast, _outprod_to_subspace + +# Add at most 20 states in each iteration +MAX_SPACE_INC = 20 + +def eigh(aop, x0, precond, tol_residual=1e-5, lindep=1e-12, nroots=1, + x0sym=None, pick=None, max_cycle=50, max_memory=MAX_MEMORY, + verbose=logger.WARN): + ''' + Solve symmetric eigenvalues. + + This solver is similar to the `linalg_helper.davidson` solver, with + optimizations for performance and wavefunction symmetry, specifically + tailored for linear response methods. + + Args: + aop : function(x) => array_like_x + The matrix-vector product operator. + x0 : 1D array or a list of 1D array + Initial guess. + precond : function(dx, e) => array_like_dx + Preconditioner to generate new trial vector. The argument dx is a + residual vector ``A*x0-e*x0``; e is the eigenvalue. + + Kwargs: + tol_residual : float + Convergence tolerance for the norm of residual vector ``A*x0-e*x0``. + lindep : float + Linear dependency threshold. The function is terminated when the + smallest eigenvalue of the metric of the trial vectors is lower + than this threshold. + nroots : int + Number of eigenvalues to be computed. + x0sym: + The symmetry label for each initial guess vectors. + pick : function(w,v,nroots) => (e[idx], w[:,idx], idx) + Function to filter eigenvalues and eigenvectors. + max_cycle : int + max number of iterations. + max_memory : int or float + Allowed memory in MB. + + Returns: + e : list of floats + Eigenvalues. + c : list of 1D arrays + Eigenvectors. + ''' + + assert callable(pick) + assert callable(precond) + + if isinstance(verbose, logger.Logger): + 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) + + space_inc = nroots if MAX_SPACE_INC is None else MAX_SPACE_INC + x0_size = x0.shape[1] + max_space = int(max_memory*1e6/8/x0_size / 2 - nroots - space_inc) + if max_space < nroots * 2 < x0_size: + log.warn('Not enough memory to store trial space in _lr_eig.eigh') + max_space = max(max_space, nroots * 2) + max_space = min(max_space, x0_size) + log.debug(f'Set max_space {max_space}, space_inc {space_inc}') + + xs = np.zeros((0, x0_size)) + ax = np.zeros((0, x0_size)) + e = w = v = None + conv_last = conv = np.zeros(nroots, dtype=bool) + xt = x0 + + if x0sym is not None: + xt_ir = np.asarray(x0sym) + xs_ir = np.array([], dtype=xt_ir.dtype) + + for icyc in range(max_cycle): + xt, xt_idx = _qr(xt, lindep) + # Generate at most space_inc trial vectors + xt = xt[:space_inc] + xt_idx = xt_idx[:space_inc] + + row0 = len(xs) + axt = aop(xt) + xs = np.vstack([xs, xt]) + ax = np.vstack([ax, axt]) + if x0sym is not None: + xs_ir = np.hstack([xs_ir, xt_ir[xt_idx]]) + + # Compute heff = xs.conj().dot(ax.T) + if w is None: + heff = xs.conj().dot(ax.T) + else: + hsub = xt.conj().dot(ax.T) + heff = np.block([[np.diag(w), hsub[:,:row0].conj().T], + [hsub[:,:row0], hsub[:,row0:]]]) + + if x0sym is None: + w, v = scipy.linalg.eigh(heff) + else: + # Diagonalize within eash symmetry sectors + row1 = len(xs) + w = np.empty(row1) + v = np.zeros((row1, row1)) + v_ir = [] + i1 = 0 + for ir in set(xs_ir): + idx = np.where(xs_ir == ir)[0] + i0, i1 = i1, i1 + idx.size + w_sub, v_sub = scipy.linalg.eigh(heff[idx[:,None],idx]) + w[i0:i1] = w_sub + v[idx,i0:i1] = v_sub + v_ir.append([ir] * idx.size) + w_idx = np.argsort(w) + w = w[w_idx] + v = v[:,w_idx] + xs_ir = np.hstack(v_ir)[w_idx] + + w, v, idx = pick(w, v, nroots, locals()) + if x0sym is not None: + xs_ir = xs_ir[idx] + if len(w) == 0: + raise RuntimeError('Not enough eigenvalues') + + e, elast = w[:nroots], e + if elast is None: + de = e + elif elast.size != e.size: + log.debug('Number of roots different from the previous step (%d,%d)', + e.size, elast.size) + de = e + else: + # mapping to previous eigenvectors + vlast = np.eye(nroots) + elast, conv_last = _sort_elast(elast, conv, vlast, + v[:nroots,:nroots], log) + de = e - elast + + xs = v.T.dot(xs) + ax = v.T.dot(ax) + if len(xs) * 2 > max_space: + row0 = max(nroots, max_space-space_inc) + xs = xs[:row0] + ax = ax[:row0] + w = w[:row0] + if x0sym is not None: + xs_ir = xs_ir[:row0] + + t_size = max(nroots, max_space-len(xs)) + xt = -w[:t_size,None] * xs[:t_size] + xt += ax[:t_size] + if x0sym is not None: + xt_ir = xs_ir[:t_size] + + dx_norm = np.linalg.norm(xt, axis=1) + max_dx_norm = max(dx_norm[:nroots]) + conv = dx_norm[:nroots] < tol_residual + for k, ek in enumerate(e[:nroots]): + if conv[k] and not conv_last[k]: + log.debug('root %d converged |r|= %4.3g e= %s max|de|= %4.3g', + k, dx_norm[k], ek, de[k]) + ide = np.argmax(abs(de)) + if all(conv): + log.debug('converged %d %d |r|= %4.3g e= %s max|de|= %4.3g', + icyc, len(xs), max_dx_norm, e, de[ide]) + break + + # remove subspace linear dependency + for k, xk in enumerate(xt): + if dx_norm[k] > tol_residual: + xt[k] = precond(xk, e[0]) + xt -= xs.conj().dot(xt.T).T.dot(xs) + xt_norm = np.linalg.norm(xt, axis=1) + + remaining = [] + for k, xk in enumerate(xt): + if dx_norm[k] > tol_residual and xt_norm[k]**2 > lindep: + xt[k] /= xt_norm[k] + remaining.append(k) + if len(remaining) == 0: + log.debug(f'Linear dependency in trial subspace. |r| for each state {dx_norm}') + break + + xt = xt[remaining] + if x0sym is not None: + xt_ir = xt_ir[remaining] + norm_min = xt_norm[remaining].min() + log.debug('davidson %d %d |r|= %4.3g e= %s max|de|= %4.3g lindep= %4.3g', + icyc, len(xs), max_dx_norm, e, de[ide], norm_min) + + x0 = xs[:nroots] + # Check whether the solver finds enough eigenvectors. + if len(x0) < min(x0_size, nroots): + log.warn(f'Not enough eigenvectors (len(x0)={len(x0)}, nroots={nroots})') + + return conv, e, x0 + +def eig(aop, x0, precond, tol_residual=1e-5, nroots=1, x0sym=None, pick=None, + max_cycle=50, max_memory=MAX_MEMORY, lindep=1e-12, verbose=logger.WARN): + ''' + Solver for linear response eigenvalues + [ A B] [X] = w [X] + [-B* -A*] [Y] [Y] + + subject to normalization X^2 - Y^2 = 1 + + Reference: + Olsen, Jensen, and Jorgenson, J Comput Phys, 74, 265, + DOI: 10.1016/0021-9991(88)90081-2 + + Args: + aop : function(x) => array_like_x + The matrix-vector product operator. + x0 : 1D array or a list of 1D array + Initial guess. + precond : function(dx, e) => array_like_dx + Preconditioner to generate new trial vector. + The argument dx is a residual vector ``A*x0-e*x0``; e is the eigenvalue. + + Kwargs: + tol_residual : float + Convergence tolerance for the norm of residual vector ``A*x0-e*x0``. + lindep : float + Linear dependency threshold. + nroots : int + Number of eigenvalues to be computed. + x0sym: + The symmetry label for each initial guess vectors. + pick : function(w,v,nroots) => (e[idx], w[:,idx], idx) + Function to filter eigenvalues and eigenvectors. + max_cycle : int + max number of iterations. + max_memory : int or float + Allowed memory in MB. + + Returns: + e : list of floats + Eigenvalues. + c : list of 1D arrays + Eigenvectors. + ''' + + assert callable(pick) + assert callable(precond) + + if isinstance(verbose, logger.Logger): + 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) + space_inc = nroots + 2 + + x0_size = x0.shape[1] + max_space = int(max_memory*1e6/8/(2*x0_size) / 2 - space_inc) + if max_space < nroots * 2 < x0_size: + log.warn('Not enough memory to store trial space in _lr_eig.eig') + max_space = space_inc * 2 + max_space = max(max_space, nroots * 2) + max_space = min(max_space, x0_size) + log.debug(f'Set max_space {max_space}') + + heff = None + e = None + v = None + conv_last = conv = np.zeros(nroots, dtype=bool) + + if x0sym is not None: + x0_ir = np.asarray(x0sym) + + half_size = x0[0].size // 2 + fresh_start = True + for icyc in range(max_cycle): + if fresh_start: + xs = np.zeros((0, x0_size)) + ax = np.zeros((0, x0_size)) + row1 = 0 + xt = x0 + if x0sym is not None: + xs_ir = [] + xt_ir = x0_ir + + if x0sym is None: + xt = _symmetric_orth(xt) + else: + xt_orth = [] + xt_orth_ir = [] + for ir in set(xt_ir): + idx = np.where(xt_ir == ir)[0] + xt_sub = _symmetric_orth(xt[idx]) + xt_orth.append(xt_sub) + xt_orth_ir.append([ir] * len(xt_sub)) + xt = np.vstack(xt_orth) + xs_ir = np.hstack([xs_ir, *xt_orth_ir]) + xt_orth = xt_orth_ir = None + + axt = aop(xt) + xs = np.vstack([xs, xt]) + ax = np.vstack([ax, axt]) + row0, row1 = row1, row1+len(xt) + + if heff is None: + dtype = np.result_type(axt, xt) + heff = np.empty((max_space*2,max_space*2), dtype=dtype) + + h11 = xs[:row0].conj().dot(axt.T) + h21 = _conj_dot(xs[:row0], axt) + heff[0:row0*2:2, row0*2+0:row1*2:2] = h11 + heff[1:row0*2:2, row0*2+0:row1*2:2] = h21 + heff[0:row0*2:2, row0*2+1:row1*2:2] = -h21.conj() + heff[1:row0*2:2, row0*2+1:row1*2:2] = -h11.conj() + + h11 = xt.conj().dot(ax.T) + h21 = _conj_dot(xt, ax) + heff[row0*2+0:row1*2:2, 0:row1*2:2] = h11 + heff[row0*2+1:row1*2:2, 0:row1*2:2] = h21 + heff[row0*2+0:row1*2:2, 1:row1*2:2] = -h21.conj() + heff[row0*2+1:row1*2:2, 1:row1*2:2] = -h11.conj() + + if x0sym is None: + w, v = scipy.linalg.eig(heff[:row1*2,:row1*2]) + else: + # Diagonalize within eash symmetry sectors + xs_ir2 = np.repeat(xs_ir, 2) + w = np.empty(row1*2, dtype=np.complex128) + v = np.zeros((row1*2, row1*2), dtype=np.complex128) + v_ir = [] + i1 = 0 + for ir in set(xs_ir): + idx = np.where(xs_ir2 == ir)[0] + i0, i1 = i1, i1 + idx.size + w_sub, v_sub = scipy.linalg.eig(heff[idx[:,None],idx]) + w[i0:i1] = w_sub + v[idx,i0:i1] = v_sub + v_ir.append([ir] * idx.size) + v_ir = np.hstack(v_ir) + + w, v, idx = pick(w, v, nroots, locals()) + if x0sym is not None: + v_ir = v_ir[idx] + if len(w) == 0: + raise RuntimeError('Not enough eigenvalues') + + w, e, elast = w[:space_inc], w[:nroots], e + v, vlast = v[:,:space_inc], v[:,:nroots] + if not fresh_start: + elast, conv_last = _sort_elast(elast, conv, vlast, v[:,:nroots], log) + + if elast is None: + de = e + elif elast.size != e.size: + log.debug('Number of roots different from the previous step (%d,%d)', + e.size, elast.size) + de = e + else: + de = e - elast + + x0 = _gen_x0(v, xs) + ax0 = xt = _gen_ax0(v, ax) + xt -= w[:,None] * x0 + ax0 = None + if x0sym is not None: + xt_ir = x0_ir = v_ir[:space_inc] + + dx_norm = np.linalg.norm(xt, axis=1) + max_dx_norm = max(dx_norm[:nroots]) + conv = dx_norm[:nroots] < tol_residual + for k, ek in enumerate(e[:nroots]): + if conv[k] and not conv_last[k]: + log.debug('root %d converged |r|= %4.3g e= %s max|de|= %4.3g', + k, dx_norm[k], ek, de[k]) + ide = np.argmax(abs(de)) + if all(conv): + log.debug('converged %d %d |r|= %4.3g e= %s max|de|= %4.3g', + icyc, len(xs), max_dx_norm, e, de[ide]) + break + + # remove subspace linear dependency + for k, xk in enumerate(xt): + if dx_norm[k] > tol_residual: + xt[k] = precond(xk, e[0]) + + xt -= xs.conj().dot(xt.T).T.dot(xs) + c = _conj_dot(xs, xt) + xt[:,:half_size] -= c.T.dot(xs[:,half_size:].conj()) + xt[:,half_size:] -= c.T.dot(xs[:,:half_size].conj()) + xt_norm = np.linalg.norm(xt, axis=1) + + remaining = [] + for k, xk in enumerate(xt): + if (dx_norm[k] > tol_residual and + xt_norm[k] > tol_residual and xt_norm[k]**2 > lindep): + xt[k] /= xt_norm[k] + remaining.append(k) + if len(remaining) == 0: + log.debug(f'Linear dependency in trial subspace. |r| for each state {dx_norm}') + break + + xt = xt[remaining] + if x0sym is not None: + xt_ir = xt_ir[remaining] + norm_min = xt_norm[remaining].min() + log.debug('davidson %d %d |r|= %4.3g e= %s max|de|= %4.3g lindep= %4.3g', + icyc, len(xs), max_dx_norm, e, de[ide], norm_min) + + fresh_start = len(xs)+space_inc > max_space + + # Check whether the solver finds enough eigenvectors. + h_dim = x0[0].size + if len(x0) < min(h_dim, nroots): + log.warn(f'Not enough eigenvectors (len(x0)={len(x0)}, nroots={nroots})') + + return conv[:nroots], e[:nroots], x0[:nroots] + +def _gen_x0(v, xs): + out = _outprod_to_subspace(v[::2], xs) + out_conj = _outprod_to_subspace(v[1::2].conj(), xs) + n = out.shape[1] // 2 + # v[1::2] * xs.conj() = (v[1::2].conj() * xs).conj() + out[:,:n] += out_conj[:,n:].conj() + out[:,n:] += out_conj[:,:n].conj() + return out + +def _gen_ax0(v, xs): + out = _outprod_to_subspace(v[::2], xs) + out_conj = _outprod_to_subspace(v[1::2].conj(), xs) + n = out.shape[1] // 2 + out[:,:n] -= out_conj[:,n:].conj() + out[:,n:] -= out_conj[:,:n].conj() + return out + +def _conj_dot(xi, xj): + '''Dot product between the conjugated basis of xi and xj. + The conjugated basis of xi is np.hstack([xi[half:], xi[:half]]).conj() + ''' + n = xi.shape[-1] // 2 + 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 _symmetric_orth(xt, lindep=1e-6): + xt = np.asarray(xt) + n, m = xt.shape + if n == 0: + raise RuntimeError('Linear dependency in trial bases') + m = m // 2 + # The conjugated basis np.hstack([xt[:,m:], xt[:,:m]]).conj() + s11 = xt.conj().dot(xt.T) + s21 = _conj_dot(xt, xt) + s = np.block([[s11, s21.conj().T], + [s21, s11.conj() ]]) + e, c = scipy.linalg.eigh(s) + if e[0] < lindep: + if n == 1: + return xt + return _symmetric_orth(xt[:-1], lindep) + + c_orth = (c * e**-.5).dot(c[:n].conj().T) + x_orth = c_orth[:n].T.dot(xt) + # Contribution from the conjugated basis + x_orth[:,:m] += c_orth[n:].T.dot(xt[:,m:].conj()) + x_orth[:,m:] += c_orth[n:].T.dot(xt[:,:m].conj()) + return x_orth diff --git a/pyscf/tdscf/dhf.py b/pyscf/tdscf/dhf.py index 6e54e0a664..00cf71a06a 100644 --- a/pyscf/tdscf/dhf.py +++ b/pyscf/tdscf/dhf.py @@ -23,11 +23,11 @@ from functools import reduce import numpy from pyscf import lib -from pyscf import dft +from pyscf import scf from pyscf import ao2mo from pyscf.lib import logger from pyscf.tdscf import rhf -from pyscf.data import nist +from pyscf.tdscf._lr_eig import eigh as lr_eigh, eig as lr_eig from pyscf import __config__ OUTPUT_THRESHOLD = getattr(__config__, 'tdscf_uhf_get_nto_threshold', 0.3) @@ -120,7 +120,7 @@ def add_hf_(a, b, hyb=1): b = b - numpy.einsum('jaib->iajb', eri_mo[:nocc,nocc:,:nocc,nocc:]) * hyb return a, b - if isinstance(mf, dft.KohnShamDFT): + if isinstance(mf, scf.hf.KohnShamDFT): from pyscf.dft import xc_deriv ni = mf._numint ni.libxc.test_deriv_order(mf.xc, 2, raise_error=True) @@ -398,6 +398,9 @@ def nuc_grad_method(self): @lib.with_doc(rhf.TDA.__doc__) class TDA(TDBase): + + singlet = None + def gen_vind(self, mf=None): '''Generate function to compute Ax''' if mf is None: @@ -449,14 +452,10 @@ def pickeig(w, v, nroots, envs): if x0 is None: x0 = self.init_guess(self._scf, self.nstates) - # FIXME: Is it correct to call davidson1 for complex integrals? - self.converged, self.e, x1 = \ - lib.davidson1(vind, x0, precond, - tol=self.conv_tol, - nroots=nstates, lindep=self.lindep, - max_cycle=self.max_cycle, - max_space=self.max_space, pick=pickeig, - verbose=log) + self.converged, self.e, x1 = lr_eigh( + vind, x0, precond, tol_residual=self.conv_tol, lindep=self.lindep, + nroots=nstates, pick=pickeig, max_cycle=self.max_cycle, + max_memory=self.max_memory, verbose=log) nocc = (self._scf.mo_occ>0).sum() nmo = self._scf.mo_occ.size @@ -525,6 +524,9 @@ def vind(xys): class TDHF(TDBase): + + singlet = None + @lib.with_doc(gen_tdhf_operation.__doc__) def gen_vind(self, mf=None): if mf is None: @@ -534,7 +536,9 @@ def gen_vind(self, mf=None): def init_guess(self, mf, nstates=None, wfnsym=None): x0 = TDA.init_guess(self, mf, nstates, wfnsym) y0 = numpy.zeros_like(x0) - return numpy.asarray(numpy.block([[x0, y0], [y0, x0.conj()]])) + return numpy.hstack([x0, y0]) + + get_precond = rhf.TDHF.get_precond def kernel(self, x0=None, nstates=None): '''TDHF diagonalization with non-Hermitian eigenvalue solver @@ -562,13 +566,10 @@ def pickeig(w, v, nroots, envs): if x0 is None: x0 = self.init_guess(self._scf, self.nstates) - self.converged, w, x1 = \ - lib.davidson_nosym1(vind, x0, precond, - tol=self.conv_tol, - nroots=nstates, lindep=self.lindep, - max_cycle=self.max_cycle, - max_space=self.max_space, pick=pickeig, - verbose=log) + self.converged, w, x1 = lr_eig( + vind, x0, precond, tol_residual=self.conv_tol, lindep=self.lindep, + nroots=nstates, pick=pickeig, max_cycle=self.max_cycle, + max_memory=self.max_memory, verbose=log) nocc = (self._scf.mo_occ>0).sum() nmo = self._scf.mo_occ.size @@ -590,7 +591,6 @@ def norm_xy(z): self._finalize() return self.e, self.xy -from pyscf import scf scf.dhf.DHF.TDA = scf.dhf.RDHF.TDA = lib.class_as_method(TDA) scf.dhf.DHF.TDHF = scf.dhf.RDHF.TDHF = lib.class_as_method(TDHF) diff --git a/pyscf/tdscf/dks.py b/pyscf/tdscf/dks.py index 1f7bdf0563..fec0ad1e85 100644 --- a/pyscf/tdscf/dks.py +++ b/pyscf/tdscf/dks.py @@ -17,12 +17,8 @@ # -import numpy from pyscf import lib from pyscf.tdscf import dhf -from pyscf.data import nist -from pyscf.dft.rks import KohnShamDFT -from pyscf import __config__ class TDA(dhf.TDA): diff --git a/pyscf/tdscf/ghf.py b/pyscf/tdscf/ghf.py index 69f70ea13b..9bd02b74ec 100644 --- a/pyscf/tdscf/ghf.py +++ b/pyscf/tdscf/ghf.py @@ -25,14 +25,13 @@ from functools import reduce import numpy from pyscf import lib -from pyscf import dft -from pyscf.dft import numint +from pyscf import scf from pyscf import ao2mo from pyscf import symm from pyscf.lib import logger from pyscf.tdscf import rhf +from pyscf.tdscf._lr_eig import eigh as lr_eigh, eig as lr_eig from pyscf.scf import ghf_symm -from pyscf.data import nist from pyscf import __config__ OUTPUT_THRESHOLD = getattr(__config__, 'tdscf_rhf_get_nto_threshold', 0.3) @@ -62,9 +61,7 @@ def gen_tda_operation(mf, fock_ao=None, wfnsym=None): if isinstance(wfnsym, str): wfnsym = symm.irrep_name2id(mol.groupname, wfnsym) wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsym = ghf_symm.get_orbsym(mol, mo_coeff) - orbsym_in_d2h = numpy.asarray(orbsym) % 10 # convert to D2h irreps - sym_forbid = (orbsym_in_d2h[occidx,None] ^ orbsym_in_d2h[viridx]) != wfnsym + sym_forbid = _get_x_sym_table(mf) != wfnsym if fock_ao is None: e_ia = hdiag = mo_energy[viridx] - mo_energy[occidx,None] @@ -102,6 +99,14 @@ def vind(zs): return vind, hdiag gen_tda_hop = gen_tda_operation +def _get_x_sym_table(mf): + '''Irrep (up to D2h symmetry) of each coefficient in X[nocc,nvir]''' + mol = mf.mol + mo_occ = mf.mo_occ + orbsym = ghf_symm.get_orbsym(mol, mf.mo_coeff) + orbsym = numpy.asarray(orbsym) % 10 # convert to D2h irreps + return orbsym[mo_occ==1,None] ^ orbsym[mo_occ==0] + def get_ab(mf, mo_energy=None, mo_coeff=None, mo_occ=None): r'''A and B matrices for TDDFT response function. @@ -151,7 +156,7 @@ def add_hf_(a, b, hyb=1): b -= numpy.einsum('jaib->iajb', eri_mo[:nocc,nocc:,:nocc,nocc:]) * hyb return a, b - if isinstance(mf, dft.KohnShamDFT): + if isinstance(mf, scf.hf.KohnShamDFT): from pyscf.dft import xc_deriv ni = mf._numint ni.libxc.test_deriv_order(mf.xc, 2, raise_error=True) @@ -379,7 +384,7 @@ def gen_vind(self, mf=None): mf = self._scf return gen_tda_hop(mf, wfnsym=self.wfnsym) - def init_guess(self, mf, nstates=None, wfnsym=None): + def init_guess(self, mf, nstates=None, wfnsym=None, return_symmetry=False): if nstates is None: nstates = self.nstates if wfnsym is None: wfnsym = self.wfnsym @@ -387,27 +392,36 @@ def init_guess(self, mf, nstates=None, wfnsym=None): mo_occ = mf.mo_occ occidx = numpy.where(mo_occ==1)[0] viridx = numpy.where(mo_occ==0)[0] - e_ia = mo_energy[viridx] - mo_energy[occidx,None] - - if wfnsym is not None and mf.mol.symmetry: - if isinstance(wfnsym, str): - wfnsym = symm.irrep_name2id(mf.mol.groupname, wfnsym) - wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsym = ghf_symm.get_orbsym(mf.mol, mf.mo_coeff) - orbsym_in_d2h = numpy.asarray(orbsym) % 10 # convert to D2h irreps - e_ia[(orbsym_in_d2h[occidx,None] ^ orbsym_in_d2h[viridx]) != wfnsym] = 1e99 - + e_ia = (mo_energy[viridx] - mo_energy[occidx,None]).ravel() nov = e_ia.size nstates = min(nstates, nov) - e_ia = e_ia.ravel() - e_threshold = numpy.sort(e_ia)[nstates-1] + + if (wfnsym is not None or return_symmetry) and mf.mol.symmetry: + x_sym = _get_x_sym_table(mf).ravel() + if wfnsym is not None: + if isinstance(wfnsym, str): + wfnsym = symm.irrep_name2id(mf.mol.groupname, wfnsym) + wfnsym = wfnsym % 10 # convert to D2h subgroup + e_ia[x_sym != wfnsym] = 1e99 + nov_allowed = numpy.count_nonzero(x_sym == wfnsym) + nstates = min(nstates, nov_allowed) + + e_threshold = numpy.partition(e_ia, nstates-1)[nstates-1] e_threshold += self.deg_eia_thresh idx = numpy.where(e_ia <= e_threshold)[0] x0 = numpy.zeros((idx.size, nov)) for i, j in enumerate(idx): x0[i, j] = 1 # Koopmans' excitations - return x0 + + if return_symmetry: + if mf.mol.symmetry: + x0sym = x_sym[idx] + else: + x0sym = None + return x0, x0sym + else: + return x0 def kernel(self, x0=None, nstates=None): '''TDA diagonalization solver @@ -419,6 +433,7 @@ def kernel(self, x0=None, nstates=None): nstates = self.nstates else: self.nstates = nstates + mol = self.mol log = logger.Logger(self.stdout, self.verbose) @@ -429,17 +444,18 @@ def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + x0sym = None if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) + x0, x0sym = self.init_guess( + self._scf, self.nstates, return_symmetry=True) + elif mol.symmetry: + x_sym = _get_x_sym_table(self._scf).ravel() + x0sym = [rhf._guess_wfnsym_id(self, x_sym, x) for x in x0] - # FIXME: Is it correct to call davidson1 for complex integrals - self.converged, self.e, x1 = \ - lib.davidson1(vind, x0, precond, - tol=self.conv_tol, - nroots=nstates, lindep=self.lindep, - max_cycle=self.max_cycle, - max_space=self.max_space, pick=pickeig, - verbose=log) + self.converged, self.e, 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, + max_memory=self.max_memory, verbose=log) nocc = (self._scf.mo_occ>0).sum() nmo = self._scf.mo_occ.size @@ -479,9 +495,7 @@ def gen_tdhf_operation(mf, fock_ao=None, wfnsym=None): if isinstance(wfnsym, str): wfnsym = symm.irrep_name2id(mol.groupname, wfnsym) wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsym = ghf_symm.get_orbsym(mol, mo_coeff) - orbsym_in_d2h = numpy.asarray(orbsym) % 10 # convert to D2h irreps - sym_forbid = (orbsym_in_d2h[occidx,None] ^ orbsym_in_d2h[viridx]) != wfnsym + sym_forbid = _get_x_sym_table(mf) != wfnsym assert fock_ao is None @@ -538,10 +552,15 @@ def gen_vind(self, mf=None): mf = self._scf return gen_tdhf_operation(mf, wfnsym=self.wfnsym) - def init_guess(self, mf, nstates=None, wfnsym=None): - x0 = TDA.init_guess(self, mf, nstates, wfnsym) - y0 = numpy.zeros_like(x0) - return numpy.asarray(numpy.block([[x0, y0], [y0, x0.conj()]])) + def init_guess(self, mf, nstates=None, wfnsym=None, return_symmetry=False): + if return_symmetry: + x0, x0sym = TDA.init_guess(self, mf, nstates, wfnsym, return_symmetry) + y0 = numpy.zeros_like(x0) + return numpy.hstack([x0, y0]), x0sym + else: + x0 = TDA.init_guess(self, mf, nstates, wfnsym, return_symmetry) + y0 = numpy.zeros_like(x0) + return numpy.hstack([x0, y0]) def kernel(self, x0=None, nstates=None): '''TDHF diagonalization with non-Hermitian eigenvalue solver @@ -553,6 +572,7 @@ def kernel(self, x0=None, nstates=None): nstates = self.nstates else: self.nstates = nstates + mol = self.mol log = logger.Logger(self.stdout, self.verbose) @@ -566,16 +586,19 @@ def pickeig(w, v, nroots, envs): # FIXME: Should the amplitudes be real? It also affects x2c-tdscf return lib.linalg_helper._eigs_cmplx2real(w, v, realidx, ensure_real) + x0sym = None if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) - - self.converged, w, x1 = \ - lib.davidson_nosym1(vind, x0, precond, - tol=self.conv_tol, - nroots=nstates, lindep=self.lindep, - max_cycle=self.max_cycle, - max_space=self.max_space, pick=pickeig, - verbose=log) + x0, x0sym = self.init_guess( + self._scf, self.nstates, return_symmetry=True) + elif mol.symmetry: + x_sym = y_sym = _get_x_sym_table(self._scf).ravel() + x_sym = numpy.append(x_sym, y_sym) + x0sym = [rhf._guess_wfnsym_id(self, x_sym, x) for x in x0] + + self.converged, w, x1 = lr_eig( + vind, x0, precond, tol_residual=self.conv_tol, lindep=self.lindep, + nroots=nstates, x0sym=x0sym, pick=pickeig, max_cycle=self.max_cycle, + max_memory=self.max_memory, verbose=log) nocc = (self._scf.mo_occ>0).sum() nmo = self._scf.mo_occ.size @@ -598,6 +621,5 @@ def norm_xy(z): RPA = TDGHF = TDHF -from pyscf import scf scf.ghf.GHF.TDA = lib.class_as_method(TDA) scf.ghf.GHF.TDHF = lib.class_as_method(TDHF) diff --git a/pyscf/tdscf/gks.py b/pyscf/tdscf/gks.py index 461fc60635..3d372a6a08 100644 --- a/pyscf/tdscf/gks.py +++ b/pyscf/tdscf/gks.py @@ -20,10 +20,8 @@ import numpy from pyscf import lib from pyscf import symm -from pyscf.tdscf import ghf -from pyscf.scf import ghf_symm -from pyscf.scf import _response_functions # noqa -from pyscf.data import nist +from pyscf.tdscf import ghf, rhf +from pyscf.tdscf._lr_eig import eigh as lr_eigh from pyscf.dft.rks import KohnShamDFT from pyscf import __config__ @@ -39,6 +37,7 @@ class TDDFT(ghf.TDHF): class CasidaTDDFT(TDDFT, TDA): '''Solve the Casida TDDFT formula (A-B)(A+B)(X+Y) = (X+Y)w^2 ''' + init_guess = TDA.init_guess def gen_vind(self, mf=None): @@ -62,9 +61,7 @@ def gen_vind(self, mf=None): if isinstance(wfnsym, str): wfnsym = symm.irrep_name2id(mol.groupname, wfnsym) wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsym = ghf_symm.get_orbsym(mol, mo_coeff) - orbsym_in_d2h = numpy.asarray(orbsym) % 10 # convert to D2h irreps - sym_forbid = (orbsym_in_d2h[occidx,None] ^ orbsym_in_d2h[viridx]) != wfnsym + sym_forbid = ghf._get_x_sym_table(mf) != wfnsym e_ia = (mo_energy[viridx].reshape(-1,1) - mo_energy[occidx]).T if wfnsym is not None and mol.symmetry: @@ -101,6 +98,7 @@ def kernel(self, x0=None, nstates=None): '''TDDFT diagonalization solver ''' cpu0 = (lib.logger.process_clock(), lib.logger.perf_counter()) + mol = self.mol mf = self._scf if mf._numint.libxc.is_hybrid_xc(mf.xc): raise RuntimeError('%s cannot be used with hybrid functional' @@ -121,16 +119,18 @@ def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + x0sym = None if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) - - self.converged, w2, x1 = \ - lib.davidson1(vind, x0, precond, - tol=self.conv_tol, - nroots=nstates, lindep=self.lindep, - max_cycle=self.max_cycle, - max_space=self.max_space, pick=pickeig, - verbose=log) + x0, x0sym = self.init_guess( + self._scf, self.nstates, return_symmetry=True) + elif mol.symmetry: + x_sym = ghf._get_x_sym_table(self._scf).ravel() + x0sym = [rhf._guess_wfnsym_id(self, x_sym, x) for x in x0] + + 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, + max_memory=self.max_memory, verbose=log) mo_energy = self._scf.mo_energy mo_occ = self._scf.mo_occ diff --git a/pyscf/tdscf/rhf.py b/pyscf/tdscf/rhf.py index d2cb63e086..b083184de9 100644 --- a/pyscf/tdscf/rhf.py +++ b/pyscf/tdscf/rhf.py @@ -32,8 +32,9 @@ from pyscf import symm from pyscf.lib import logger from pyscf.scf import hf_symm -from pyscf.scf import _response_functions +from pyscf.scf import _response_functions # noqa from pyscf.data import nist +from pyscf.tdscf._lr_eig import eigh as lr_eigh, eig as lr_eig from pyscf import __config__ OUTPUT_THRESHOLD = getattr(__config__, 'tdscf_rhf_get_nto_threshold', 0.3) @@ -65,9 +66,8 @@ def gen_tda_operation(mf, fock_ao=None, singlet=True, wfnsym=None): if isinstance(wfnsym, str): wfnsym = symm.irrep_name2id(mol.groupname, wfnsym) wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsym = hf_symm.get_orbsym(mol, mo_coeff) - orbsym_in_d2h = numpy.asarray(orbsym) % 10 # convert to D2h irreps - sym_forbid = (orbsym_in_d2h[occidx,None] ^ orbsym_in_d2h[viridx]) != wfnsym + x_sym = _get_x_sym_table(mf) + sym_forbid = x_sym != wfnsym if fock_ao is None: e_ia = hdiag = mo_energy[viridx] - mo_energy[occidx,None] @@ -106,6 +106,14 @@ def vind(zs): return vind, hdiag gen_tda_hop = gen_tda_operation +def _get_x_sym_table(mf): + '''Irrep (up to D2h symmetry) of each coefficient in X[nocc,nvir]''' + mol = mf.mol + mo_occ = mf.mo_occ + orbsym = hf_symm.get_orbsym(mol, mf.mo_coeff) + orbsym = orbsym % 10 # convert to D2h irreps + return orbsym[mo_occ==2,None] ^ orbsym[mo_occ==0] + def get_ab(mf, mo_energy=None, mo_coeff=None, mo_occ=None): r'''A and B matrices for TDDFT response function. @@ -128,7 +136,6 @@ def get_ab(mf, mo_energy=None, mo_coeff=None, mo_occ=None): nvir = orbv.shape[1] nocc = orbo.shape[1] mo = numpy.hstack((orbo,orbv)) - nmo = nocc + nvir e_ia = lib.direct_sum('a-i->ia', mo_energy[viridx], mo_energy[occidx]) a = numpy.diag(e_ia.ravel()).reshape(nocc,nvir,nocc,nvir) @@ -144,7 +151,6 @@ def add_hf_(a, b, hyb=1): b -= numpy.einsum('jaib->iajb', eri_mo[:nocc,nocc:,:nocc,nocc:]) * hyb if isinstance(mf, scf.hf.KohnShamDFT): - from pyscf.dft import xc_deriv ni = mf._numint ni.libxc.test_deriv_order(mf.xc, 2, raise_error=True) if mf.do_nlc(): @@ -390,7 +396,7 @@ def analyze(tdobj, verbose=None): log.note('Excited State %3d: %12.5f eV %9.2f nm f=%.4f', i+1, e_ev[i], wave_length[i], f_oscillator[i]) else: - wfnsym = analyze_wfnsym(tdobj, x_sym, x) + wfnsym = _analyze_wfnsym(tdobj, x_sym, x) log.note('Excited State %3d: %4s %12.5f eV %9.2f nm f=%.4f', i+1, wfnsym, e_ev[i], wave_length[i], f_oscillator[i]) @@ -428,23 +434,23 @@ def analyze(tdobj, verbose=None): i+1, m[0], m[1], m[2]) return tdobj -def analyze_wfnsym(tdobj, x_sym, x): - '''Guess the wfn symmetry of TDDFT X amplitude''' - possible_sym = x_sym[(x > 1e-7) | (x < -1e-7)] - if numpy.all(possible_sym == symm.MULTI_IRREPS): - if tdobj.wfnsym is None: - wfnsym = '???' - else: - wfnsym = tdobj.wfnsym +def _analyze_wfnsym(tdobj, x_sym, x): + '''Guess the wfn symmetry of TDDFT X amplitude. Return a label''' + wfnsym = _guess_wfnsym_id(tdobj, x_sym, x) + if wfnsym == symm.MULTI_IRREPS: + wfnsym = '???' else: - ids = possible_sym[possible_sym != symm.MULTI_IRREPS] - ids = numpy.unique(ids) - if ids.size == 1: - wfnsym = symm.irrep_id2name(tdobj.mol.groupname, ids[0]) - else: - wfnsym = '???' + wfnsym = symm.irrep_id2name(tdobj.mol.groupname, wfnsym) return wfnsym +def _guess_wfnsym_id(tdobj, x_sym, x): + '''Guess the wfn symmetry of TDDFT X amplitude. Return an ID''' + possible_sym = x_sym[(x > 1e-7) | (x < -1e-7)] + wfnsym = symm.MULTI_IRREPS + ids = possible_sym[possible_sym != symm.MULTI_IRREPS] + if len(ids) > 0 and all(ids == ids[0]): + wfnsym = ids[0] + return wfnsym def transition_dipole(tdobj, xy=None): '''Transition dipole moments in the length gauge''' @@ -653,12 +659,11 @@ def __call__(self, mol_or_geom, **kwargs): class TDBase(lib.StreamObject): - conv_tol = getattr(__config__, 'tdscf_rhf_TDA_conv_tol', 1e-9) + conv_tol = getattr(__config__, 'tdscf_rhf_TDA_conv_tol', 1e-5) nstates = getattr(__config__, 'tdscf_rhf_TDA_nstates', 3) singlet = getattr(__config__, 'tdscf_rhf_TDA_singlet', True) lindep = getattr(__config__, 'tdscf_rhf_TDA_lindep', 1e-12) level_shift = getattr(__config__, 'tdscf_rhf_TDA_level_shift', 0) - max_space = getattr(__config__, 'tdscf_rhf_TDA_max_space', 50) max_cycle = getattr(__config__, 'tdscf_rhf_TDA_max_cycle', 100) # Low excitation filter to avoid numerical instability positive_eig_threshold = getattr(__config__, 'tdscf_rhf_TDDFT_positive_eig_threshold', 1e-3) @@ -666,7 +671,7 @@ class TDBase(lib.StreamObject): deg_eia_thresh = getattr(__config__, 'tdscf_rhf_TDDFT_deg_eia_thresh', 1e-3) _keys = { - 'conv_tol', 'nstates', 'singlet', 'lindep', 'level_shift', 'max_space', + 'conv_tol', 'nstates', 'singlet', 'lindep', 'level_shift', 'max_cycle', 'mol', 'chkfile', 'wfnsym', 'converged', 'e', 'xy', } @@ -714,7 +719,6 @@ def dump_flags(self, verbose=None): log.info('conv_tol = %g', self.conv_tol) log.info('eigh lindep = %g', self.lindep) log.info('eigh level_shift = %g', self.level_shift) - log.info('eigh max_space = %d', self.max_space) log.info('eigh max_cycle = %d', self.max_cycle) log.info('chkfile = %s', self.chkfile) log.info('max_memory %d MB (current use %d MB)', @@ -743,7 +747,7 @@ def get_ab(self, mf=None): return get_ab(mf) def get_precond(self, hdiag): - def precond(x, e, x0): + def precond(x, e, *args): diagd = hdiag - (e-self.level_shift) diagd[abs(diagd)<1e-8] = 1e-8 return x/diagd @@ -801,13 +805,25 @@ class TDA(TDBase): excited state. (X,Y) are normalized to 1/2 in RHF/RKS methods and normalized to 1 for UHF/UKS methods. In the TDA calculation, Y = 0. ''' + def gen_vind(self, mf=None): '''Generate function to compute Ax''' if mf is None: mf = self._scf return gen_tda_hop(mf, singlet=self.singlet, wfnsym=self.wfnsym) - def init_guess(self, mf, nstates=None, wfnsym=None): + def init_guess(self, mf, nstates=None, wfnsym=None, return_symmetry=False): + ''' + Generate initial guess for TDA + + Kwargs: + nstates : int + The number of initial guess vectors. + wfnsym : int or str + The irrep label or ID of the wavefunction. + return_symmetry : bool + Whether to return symmetry labels for initial guess vectors. + ''' if nstates is None: nstates = self.nstates if wfnsym is None: wfnsym = self.wfnsym @@ -815,27 +831,37 @@ def init_guess(self, mf, nstates=None, wfnsym=None): mo_occ = mf.mo_occ occidx = numpy.where(mo_occ==2)[0] viridx = numpy.where(mo_occ==0)[0] - e_ia = mo_energy[viridx] - mo_energy[occidx,None] - - if wfnsym is not None and mf.mol.symmetry: - if isinstance(wfnsym, str): - wfnsym = symm.irrep_name2id(mf.mol.groupname, wfnsym) - wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsym = hf_symm.get_orbsym(mf.mol, mf.mo_coeff) - orbsym_in_d2h = numpy.asarray(orbsym) % 10 # convert to D2h irreps - e_ia[(orbsym_in_d2h[occidx,None] ^ orbsym_in_d2h[viridx]) != wfnsym] = 1e99 - + e_ia = (mo_energy[viridx] - mo_energy[occidx,None]).ravel() nov = e_ia.size nstates = min(nstates, nov) - e_ia = e_ia.ravel() - e_threshold = numpy.sort(e_ia)[nstates-1] + + if (wfnsym is not None or return_symmetry) and mf.mol.symmetry: + x_sym = _get_x_sym_table(mf).ravel() + if wfnsym is not None: + if isinstance(wfnsym, str): + wfnsym = symm.irrep_name2id(mf.mol.groupname, wfnsym) + wfnsym = wfnsym % 10 # convert to D2h subgroup + e_ia[x_sym != wfnsym] = 1e99 + nov_allowed = numpy.count_nonzero(x_sym == wfnsym) + nstates = min(nstates, nov_allowed) + + # Find the nstates-th lowest energy gap + e_threshold = numpy.partition(e_ia, nstates-1)[nstates-1] e_threshold += self.deg_eia_thresh idx = numpy.where(e_ia <= e_threshold)[0] x0 = numpy.zeros((idx.size, nov)) for i, j in enumerate(idx): x0[i, j] = 1 # Koopmans' excitations - return x0 + + if return_symmetry: + if mf.mol.symmetry: + x0sym = x_sym[idx] + else: + x0sym = None + return x0, x0sym + else: + return x0 def kernel(self, x0=None, nstates=None): '''TDA diagonalization solver @@ -847,6 +873,7 @@ def kernel(self, x0=None, nstates=None): nstates = self.nstates else: self.nstates = nstates + mol = self.mol log = logger.Logger(self.stdout, self.verbose) @@ -857,21 +884,23 @@ def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + x0sym = None if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) + x0, x0sym = self.init_guess( + self._scf, self.nstates, return_symmetry=True) + elif mol.symmetry: + x_sym = _get_x_sym_table(self._scf).ravel() + x0sym = [_guess_wfnsym_id(self, x_sym, x) for x in x0] - self.converged, self.e, x1 = \ - lib.davidson1(vind, x0, precond, - tol=self.conv_tol, - nroots=nstates, lindep=self.lindep, - max_cycle=self.max_cycle, - max_space=self.max_space, pick=pickeig, - verbose=log) + self.converged, self.e, 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, + max_memory=self.max_memory, verbose=log) nocc = (self._scf.mo_occ>0).sum() nmo = self._scf.mo_occ.size nvir = nmo - nocc -# 1/sqrt(2) because self.x is for alpha excitation amplitude and 2(X^+*X) = 1 + # 1/sqrt(2) because self.x is for alpha excitation and 2(X^+*X) = 1 self.xy = [(xi.reshape(nocc,nvir)*numpy.sqrt(.5),0) for xi in x1] if self.chkfile: @@ -910,9 +939,7 @@ def gen_tdhf_operation(mf, fock_ao=None, singlet=True, wfnsym=None): if isinstance(wfnsym, str): wfnsym = symm.irrep_name2id(mol.groupname, wfnsym) wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsym = hf_symm.get_orbsym(mol, mo_coeff) - orbsym_in_d2h = numpy.asarray(orbsym) % 10 # convert to D2h irreps - sym_forbid = (orbsym_in_d2h[occidx,None] ^ orbsym_in_d2h[viridx]) != wfnsym + sym_forbid = _get_x_sym_table(mf) != wfnsym assert fock_ao is None @@ -960,12 +987,12 @@ def vind(xys): return vind, hdiag -class TDHF(TDA): +class TDHF(TDBase): '''Time-dependent Hartree-Fock Attributes: conv_tol : float - Diagonalization convergence tolerance. Default is 1e-9. + Diagonalization convergence tolerance. Default is 1e-4. nstates : int Number of TD states to be computed. Default is 3. @@ -981,16 +1008,22 @@ class TDHF(TDA): excited state. (X,Y) are normalized to 1/2 in RHF/RKS methods and normalized to 1 for UHF/UKS methods. In the TDA calculation, Y = 0. ''' + @lib.with_doc(gen_tdhf_operation.__doc__) def gen_vind(self, mf=None): if mf is None: mf = self._scf return gen_tdhf_operation(mf, singlet=self.singlet, wfnsym=self.wfnsym) - def init_guess(self, mf, nstates=None, wfnsym=None): - x0 = TDA.init_guess(self, mf, nstates, wfnsym) - y0 = numpy.zeros_like(x0) - return numpy.asarray(numpy.block([[x0, y0], [y0, x0.conj()]])) + def init_guess(self, mf, nstates=None, wfnsym=None, return_symmetry=False): + if return_symmetry: + x0, x0sym = TDA.init_guess(self, mf, nstates, wfnsym, return_symmetry) + y0 = numpy.zeros_like(x0) + return numpy.hstack([x0, y0]), x0sym + else: + x0 = TDA.init_guess(self, mf, nstates, wfnsym, return_symmetry) + y0 = numpy.zeros_like(x0) + return numpy.hstack([x0, y0]) def kernel(self, x0=None, nstates=None): '''TDHF diagonalization with non-Hermitian eigenvalue solver @@ -1002,6 +1035,7 @@ def kernel(self, x0=None, nstates=None): nstates = self.nstates else: self.nstates = nstates + mol = self.mol log = logger.Logger(self.stdout, self.verbose) @@ -1025,16 +1059,19 @@ def pickeig(w, v, nroots, envs): # approximately be used as the "real" eigen solutions. return lib.linalg_helper._eigs_cmplx2real(w, v, realidx, real_system) + x0sym = None if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) - - self.converged, w, x1 = \ - lib.davidson_nosym1(vind, x0, precond, - tol=self.conv_tol, - nroots=nstates, lindep=self.lindep, - max_cycle=self.max_cycle, - max_space=self.max_space, pick=pickeig, - verbose=log) + x0, x0sym = self.init_guess( + self._scf, self.nstates, return_symmetry=True) + elif mol.symmetry: + x_sym = y_sym = _get_x_sym_table(self._scf).ravel() + x_sym = numpy.append(x_sym, y_sym) + x0sym = [_guess_wfnsym_id(self, x_sym, x) for x in x0] + + self.converged, w, x1 = lr_eig( + vind, x0, precond, tol_residual=self.conv_tol, lindep=self.lindep, + nroots=nstates, x0sym=x0sym, pick=pickeig, max_cycle=self.max_cycle, + max_memory=self.max_memory, verbose=log) nocc = (self._scf.mo_occ>0).sum() nmo = self._scf.mo_occ.size diff --git a/pyscf/tdscf/rks.py b/pyscf/tdscf/rks.py index 70073b1085..5b22e06449 100644 --- a/pyscf/tdscf/rks.py +++ b/pyscf/tdscf/rks.py @@ -25,8 +25,7 @@ from pyscf import lib from pyscf import symm from pyscf.tdscf import rhf -from pyscf.scf import hf_symm -from pyscf.data import nist +from pyscf.tdscf._lr_eig import eigh as lr_eigh from pyscf.dft.rks import KohnShamDFT from pyscf import __config__ @@ -46,6 +45,7 @@ def nuc_grad_method(self): class CasidaTDDFT(TDDFT, TDA): '''Solve the Casida TDDFT formula (A-B)(A+B)(X+Y) = (X+Y)w^2 ''' + init_guess = TDA.init_guess def gen_vind(self, mf=None): @@ -71,9 +71,7 @@ def gen_vind(self, mf=None): if isinstance(wfnsym, str): wfnsym = symm.irrep_name2id(mol.groupname, wfnsym) wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsym = hf_symm.get_orbsym(mol, mo_coeff) - orbsym_in_d2h = numpy.asarray(orbsym) % 10 # convert to D2h irreps - sym_forbid = (orbsym_in_d2h[occidx,None] ^ orbsym_in_d2h[viridx]) != wfnsym + sym_forbid = rhf._get_x_sym_table(mf) != wfnsym e_ia = (mo_energy[viridx].reshape(-1,1) - mo_energy[occidx]).T if wfnsym is not None and mol.symmetry: @@ -105,6 +103,7 @@ def kernel(self, x0=None, nstates=None): '''TDDFT diagonalization solver ''' cpu0 = (lib.logger.process_clock(), lib.logger.perf_counter()) + mol = self.mol mf = self._scf if mf._numint.libxc.is_hybrid_xc(mf.xc): raise RuntimeError('%s cannot be used with hybrid functional' @@ -125,16 +124,18 @@ def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + x0sym = None if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) - - self.converged, w2, x1 = \ - lib.davidson1(vind, x0, precond, - tol=self.conv_tol, - nroots=nstates, lindep=self.lindep, - max_cycle=self.max_cycle, - max_space=self.max_space, pick=pickeig, - verbose=log) + x0, x0sym = self.init_guess( + self._scf, self.nstates, return_symmetry=True) + elif mol.symmetry: + x_sym = rhf._get_x_sym_table(mf).ravel() + x0sym = [rhf._guess_wfnsym_id(self, x_sym, x) for x in x0] + + 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, + max_memory=self.max_memory, verbose=log) mo_energy = self._scf.mo_energy mo_occ = self._scf.mo_occ diff --git a/pyscf/tdscf/test/test_tddks.py b/pyscf/tdscf/test/test_tddks.py index 7b8b093330..e548377def 100644 --- a/pyscf/tdscf/test/test_tddks.py +++ b/pyscf/tdscf/test/test_tddks.py @@ -114,13 +114,15 @@ def test_mcol_mgga_ab_ks(self): mcol_m06l = dft.UDKS(mol).set(xc='m06l', collinear='mcol') mcol_m06l._numint.spin_samples = 6 mcol_m06l.__dict__.update(scf.chkfile.load(mf_lda.chkfile, 'scf')) - self._check_against_ab_ks(mcol_m06l.TDDFT(), 14.934345395514491, 9.539340104227188) + self._check_against_ab_ks(mcol_m06l.TDDFT()) - def _check_against_ab_ks(self, td, refa, refb): + def _check_against_ab_ks(self, td, refa=None, refb=None): mf = td._scf a, b = td.get_ab() - self.assertAlmostEqual(lib.fp(abs(a)), refa, 4) - self.assertAlmostEqual(lib.fp(abs(b)), refb, 4) + if refa is not None: + self.assertAlmostEqual(lib.fp(abs(a)), refa, 4) + if refb is not None: + self.assertAlmostEqual(lib.fp(abs(b)), refb, 4) ftda = mf.TDA().gen_vind()[0] ftdhf = td.gen_vind()[0] n2c = mf.mo_occ.size // 2 diff --git a/pyscf/tdscf/test/test_tdrhf.py b/pyscf/tdscf/test/test_tdrhf.py index b2b5ed6218..80b30d2600 100644 --- a/pyscf/tdscf/test/test_tdrhf.py +++ b/pyscf/tdscf/test/test_tdrhf.py @@ -27,6 +27,7 @@ def setUpModule(): ['H' , (0. , 0. , .917)], ['F' , (0. , 0. , 0.)], ] mol.basis = '631g' + mol.symmetry = True mol.build() mf = scf.RHF(mol).run() nstates = 5 # make sure first 3 states are converged @@ -39,6 +40,7 @@ def tearDownModule(): class KnownValues(unittest.TestCase): def test_tda_singlet(self): td = mf.TDA().set(nstates=nstates) + td.max_memory = 1e-3 e = td.kernel()[0] ref = [11.90276464, 11.90276464, 16.86036434] self.assertAlmostEqual(abs(e[:len(ref)] * 27.2114 - ref).max(), 0, 5) @@ -71,6 +73,29 @@ def test_tdhf_triplet(self): dip = td.transition_dipole() self.assertAlmostEqual(abs(dip).max(), 0, 8) + def test_symmetry_init_guess(self): + mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='631g', symmetry='D2h') + mf = mol.RHF.run() + td = mf.TDA().run(nstates=1) + self.assertAlmostEqual(td.e[0], 0.22349707455528, 7) + mo_coeff = mf.mo_coeff + mo_occ = mf.mo_occ + orbsym = tdscf.rhf.hf_symm.get_orbsym(mol, mo_coeff) + x_sym = tdscf.rhf.symm.direct_prod(orbsym[mo_occ==2], orbsym[mo_occ==0], mol.groupname) + wfnsym = tdscf.rhf._analyze_wfnsym(td, x_sym, td.xy[0][0]) + self.assertEqual(wfnsym, 'Au') + + mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='631g', symmetry=True) + mf = mol.RHF.run() + td = mf.TDA().run(nstates=1, singlet=False) + self.assertAlmostEqual(td.e[0], 0.14147328219131602, 7) + mo_coeff = mf.mo_coeff + mo_occ = mf.mo_occ + orbsym = tdscf.rhf.hf_symm.get_orbsym(mol, mo_coeff) + x_sym = tdscf.rhf.symm.direct_prod(orbsym[mo_occ==2], orbsym[mo_occ==0], mol.groupname) + wfnsym = tdscf.rhf._analyze_wfnsym(td, x_sym, td.xy[0][0]) + self.assertEqual(wfnsym, 'A1u') + if __name__ == "__main__": print("Full Tests for rhf-TDA and rhf-TDHF") diff --git a/pyscf/tdscf/test/test_tdrks.py b/pyscf/tdscf/test/test_tdrks.py index 2ca7cbb540..7bc7a228d4 100644 --- a/pyscf/tdscf/test/test_tdrks.py +++ b/pyscf/tdscf/test/test_tdrks.py @@ -37,7 +37,7 @@ def setUpModule(): mol.build() mf = scf.RHF(mol).run() - td_hf = tdscf.TDHF(mf).run() + td_hf = tdscf.TDHF(mf).run(conv_tol=1e-6) mf_lda = dft.RKS(mol) mf_lda.xc = 'lda, vwn' @@ -84,20 +84,21 @@ def test_nohbrid_lda(self): def test_nohbrid_b88p86(self): td = rks.CasidaTDDFT(mf_bp86) es = td.kernel(nstates=5)[0] * 27.2114 - self.assertAlmostEqual(lib.fp(es), -40.462005239920558, 4) + self.assertAlmostEqual(lib.fp(es), -40.4619799852133, 6) a, b = td.get_ab() ref = diagonalize(a, b, nroots=5) * 27.2114 - self.assertAlmostEqual(abs(es - ref).max(), 0, 7) + self.assertAlmostEqual(abs(es - ref).max(), 0, 6) def test_tddft_lda(self): td = rks.TDDFT(mf_lda) es = td.kernel(nstates=5)[0] * 27.2114 - self.assertAlmostEqual(lib.fp(es), -41.100806721759945, 4) + self.assertAlmostEqual(lib.fp(es), -41.100806721759945, 5) def test_tddft_b88p86(self): td = rks.TDDFT(mf_bp86) + td.conv_tol = 1e-5 es = td.kernel(nstates=5)[0] * 27.2114 - self.assertAlmostEqual(lib.fp(es), -40.462005239920558, 4) + self.assertAlmostEqual(lib.fp(es), -40.4619799852133, 6) def test_tddft_b3lyp(self): td = rks.TDDFT(mf_b3lyp) @@ -105,15 +106,16 @@ def test_tddft_b3lyp(self): self.assertAlmostEqual(lib.fp(es), -41.29609453661341, 4) a, b = td.get_ab() ref = diagonalize(a, b, nroots=5) * 27.2114 - self.assertAlmostEqual(abs(es - ref).max(), 0, 7) + self.assertAlmostEqual(abs(es - ref).max(), 0, 6) def test_tddft_camb3lyp(self): mf = mol.RKS(xc='camb3lyp').run() td = mf.TDDFT() + td.conv_tol = 1e-5 es = td.kernel(nstates=4)[0] a,b = td.get_ab() e_ref = diagonalize(a, b, 5) - self.assertAlmostEqual(abs(es[:3]-e_ref[:3]).max(), 0, 8) + self.assertAlmostEqual(abs(es[:3]-e_ref[:3]).max(), 0, 7) self.assertAlmostEqual(lib.fp(es[:3]*27.2114), 9.0054057603534, 4) def test_tda_b3lypg(self): @@ -168,10 +170,10 @@ def test_tda_b3lyp_triplet(self): def test_tda_lda_triplet(self): td = rks.TDA(mf_lda) td.singlet = False - es = td.kernel(nstates=5)[0] * 27.2114 - self.assertAlmostEqual(lib.fp(es), -39.988118769202416, 5) - ref = [9.0139312, 9.0139312, 12.42444659] - self.assertAlmostEqual(abs(es[:3] - ref).max(), 0, 4) + es = td.kernel(nstates=6)[0] * 27.2114 + self.assertAlmostEqual(lib.fp(es[[0,1,2,4,5]]), -39.988118769202416, 5) + ref = [9.0139312, 9.0139312, 12.42444659, 29.38040677, 29.63058493, 29.63058493] + self.assertAlmostEqual(abs(es - ref).max(), 0, 4) def test_tddft_b88p86_triplet(self): td = rks.TDDFT(mf_bp86) @@ -294,7 +296,7 @@ def test_ab_mgga(self): def test_nto(self): mf = scf.RHF(mol).run() - td = rks.TDA(mf).run(nstates=5) + td = rks.TDA(mf).run(conv_tol=1e-6, nstates=5) w, nto = td.get_nto(state=3) self.assertAlmostEqual(w[0], 0.98655300613468389, 7) self.assertAlmostEqual(lib.fp(w), 0.98625701534112464, 7) @@ -307,7 +309,7 @@ def test_nto(self): pmol.symmetry = True pmol.build(0, 0) mf = scf.RHF(pmol).run() - td = rks.TDA(mf).run(nstates=3) + td = rks.TDA(mf).run(conv_tol=1e-6, nstates=3) w, nto = td.get_nto(state=-1) self.assertAlmostEqual(w[0], 0.98655300613468389, 7) self.assertAlmostEqual(lib.fp(w), 0.98625701534112464, 7) @@ -364,7 +366,7 @@ def test_init(self): def test_tda_with_wfnsym(self): pmol = mol.copy() - pmol.symmetry = True + pmol.symmetry = 'C2v' pmol.build(0, 0) mf = dft.RKS(pmol).run() @@ -471,12 +473,6 @@ def test_custom_rsh(self): ref = [16.14837289, 28.01968627, 49.00854076] self.assertAlmostEqual(abs(e_td*nist.HARTREE2EV - ref).max(), 0, 4) - def test_symmetry_init_guess(self): - mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='631g', output='/dev/null', symmetry=True) - td = mol.RHF.run().TDA().run(nstates=1) - self.assertAlmostEqual(td.e[0], 0.22349707455528, 7) - # TODO: verify symmetry of td.x == A1u - if __name__ == "__main__": print("Full Tests for TD-RKS") unittest.main() diff --git a/pyscf/tdscf/test/test_tduhf.py b/pyscf/tdscf/test/test_tduhf.py index 34eaa0587b..08f8734bab 100644 --- a/pyscf/tdscf/test/test_tduhf.py +++ b/pyscf/tdscf/test/test_tduhf.py @@ -15,7 +15,7 @@ # import unittest -from pyscf import lib, gto, scf, tdscf +from pyscf import lib, gto, scf, tdscf, symm def setUpModule(): global mol, mol1, mf, mf1 @@ -25,6 +25,7 @@ def setUpModule(): ['H' , (0. , 0. , .917)], ['F' , (0. , 0. , 0.)], ] mol.basis = '631g' + mol.symmetry = True mol.build() mf = scf.UHF(mol).run(conv_tol=1e-10) @@ -56,6 +57,7 @@ def test_tdhf(self): td = mf.TDHF() td.nstates = 5 td.singlet = False + td.conv_tol = 1e-5 e = td.kernel()[0] ref = [10.89192986, 10.89192986, 11.83487865, 11.83487865, 12.6344099] self.assertAlmostEqual(abs(e * 27.2114 - ref).max(), 0, 4) @@ -74,6 +76,20 @@ def test_tdhf_triplet(self): ref = [3.31267103, 18.4954748, 20.84935404, 21.54808392] self.assertAlmostEqual(abs(e * 27.2114 - ref).max(), 0, 4) + def test_symmetry_init_guess(self): + mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='631g', symmetry=True, verbose=0) + mf = mol.UHF.run() + td = mf.TDA().run(nstates=1) + self.assertAlmostEqual(td.e[0], 0.14147328219131602, 7) + mo_coeff = mf.mo_coeff + mo_occa, mo_occb = mf.mo_occ + orbsyma, orbsymb = scf.uhf_symm.get_orbsym(mol, mo_coeff) + x_syma = symm.direct_prod(orbsyma[mo_occa==1], orbsyma[mo_occa==0], mol.groupname) + x_symb = symm.direct_prod(orbsymb[mo_occb==1], orbsymb[mo_occb==0], mol.groupname) + wfnsyma = tdscf.rhf._analyze_wfnsym(td, x_syma, td.xy[0][0][0]) + wfnsymb = tdscf.rhf._analyze_wfnsym(td, x_symb, td.xy[0][0][1]) + self.assertAlmostEqual(wfnsyma, 'A1u') + self.assertAlmostEqual(wfnsymb, 'A1u') if __name__ == "__main__": print("Full Tests for uhf-TDA and uhf-TDHF") diff --git a/pyscf/tdscf/test/test_tduks.py b/pyscf/tdscf/test/test_tduks.py index 61a70172ff..3bd42317af 100644 --- a/pyscf/tdscf/test/test_tduks.py +++ b/pyscf/tdscf/test/test_tduks.py @@ -18,7 +18,7 @@ import unittest import numpy -from pyscf import lib, gto, scf, dft +from pyscf import lib, gto, scf, dft, symm from pyscf import tdscf def diagonalize(a, b, nroots=4): @@ -65,7 +65,7 @@ def setUpModule(): mol1.build() mf_uhf = scf.UHF(mol).run() - td_hf = tdscf.TDHF(mf_uhf).run(conv_tol=1e-12) + td_hf = tdscf.TDHF(mf_uhf).run(conv_tol=1e-6) mf_lda = dft.UKS(mol).set(xc='lda', conv_tol=1e-12) mf_lda.grids.prune = None @@ -86,35 +86,35 @@ def tearDownModule(): class KnownValues(unittest.TestCase): def test_nohybrid_lda(self): - td = tdscf.uks.CasidaTDDFT(mf_lda).set(conv_tol=1e-12) + td = tdscf.uks.CasidaTDDFT(mf_lda) es = td.kernel(nstates=4)[0] a,b = td.get_ab() e_ref = diagonalize(a, b, 5) - self.assertAlmostEqual(abs(es[:3]-e_ref[:3]).max(), 0, 8) + self.assertAlmostEqual(abs(es[:3]-e_ref[:3]).max(), 0, 6) self.assertAlmostEqual(lib.fp(es[:3]*27.2114), 1.294630966929489, 4) mf = dft.UKS(mol1).run(xc='lda, vwn_rpa').run() td = mf.CasidaTDDFT() td.nstates = 5 es = td.kernel()[0] * 27.2114 - ref = [6.94083826, 7.61492553, 8.55550045, 9.36308859, 9.84896499] + ref = [6.94083826, 7.61492553, 8.55550045, 9.36308859, 9.49948318] self.assertAlmostEqual(abs(es - ref).max(), 0, 4) def test_nohybrid_b88p86(self): - td = tdscf.uks.CasidaTDDFT(mf_bp86).set(conv_tol=1e-12) + td = tdscf.uks.CasidaTDDFT(mf_bp86) es = td.kernel(nstates=4)[0] a,b = td.get_ab() e_ref = diagonalize(a, b, 5) - self.assertAlmostEqual(abs(es[:3]-e_ref[:3]).max(), 0, 8) + self.assertAlmostEqual(abs(es[:3]-e_ref[:3]).max(), 0, 6) self.assertAlmostEqual(lib.fp(es[:3]*27.2114), 1.4624730971221087, 4) def test_tddft_lda(self): - td = tdscf.uks.TDDFT(mf_lda).set(conv_tol=1e-12) + td = tdscf.uks.TDDFT(mf_lda) es = td.kernel(nstates=4)[0] * 27.2114 self.assertAlmostEqual(lib.fp(es[:3]), 1.2946309669294163, 4) def test_tddft_b88p86(self): - td = tdscf.uks.TDDFT(mf_bp86).set(conv_tol=1e-12) + td = tdscf.uks.TDDFT(mf_bp86) es = td.kernel(nstates=5)[0] * 27.2114 self.assertAlmostEqual(lib.fp(es[:3]), 1.4624730971221087, 4) ref = [2.45700922, 2.93224712, 6.19693767, 12.22264487, 13.40445012] @@ -122,11 +122,11 @@ def test_tddft_b88p86(self): mf = dft.UKS(mol1).run(xc='b88,p86').run() es = mf.TDDFT().kernel(nstates=5)[0] * 27.2114 - ref = [6.96397206, 7.70955605, 8.59882964, 9.35357180, 9.92828610] + ref = [6.96396398, 7.70954799, 8.59882244, 9.35356454, 9.69774071] self.assertAlmostEqual(abs(es - ref).max(), 0, 4) def test_tddft_b3lyp(self): - td = tdscf.uks.TDDFT(mf_b3lyp).set(conv_tol=1e-12) + td = tdscf.uks.TDDFT(mf_b3lyp) es = td.kernel(nstates=4)[0] * 27.2114 self.assertAlmostEqual(lib.fp(es[:3]), 1.2984822994759448, 4) @@ -136,16 +136,16 @@ def test_tddft_camb3lyp(self): es = td.kernel(nstates=4)[0] a,b = td.get_ab() e_ref = diagonalize(a, b, 5) - self.assertAlmostEqual(abs(es[:3]-e_ref[:3]).max(), 0, 8) + self.assertAlmostEqual(abs(es[:3]-e_ref[:3]).max(), 0, 6) self.assertAlmostEqual(lib.fp(es[:3]*27.2114), 7.69383202636, 4) def test_tda_b3lyp(self): - td = tdscf.TDA(mf_b3lyp).set(conv_tol=1e-12) + td = tdscf.TDA(mf_b3lyp) es = td.kernel(nstates=4)[0] * 27.2114 self.assertAlmostEqual(lib.fp(es[:3]), 1.4303636271767162, 4) def test_tda_lda(self): - td = tdscf.TDA(mf_lda).set(conv_tol=1e-12) + td = tdscf.TDA(mf_lda) es = td.kernel(nstates=5)[0] * 27.2114 self.assertAlmostEqual(lib.fp(es[:3]), 1.4581538269747121, 4) ref = [2.14644585, 3.27738191, 5.90913787, 12.14980714, 13.15535042] @@ -155,16 +155,15 @@ def test_tda_lda(self): td = mf.TDA() td.nstates = 5 es = td.kernel()[0] * 27.2114 - ref = [6.88046608, 7.58244885, 8.49961771, 9.30209259, 9.79775972] + ref = [6.88046608, 7.58244885, 8.49961771, 9.30209259, 9.53368005] self.assertAlmostEqual(abs(es - ref).max(), 0, 4) def test_tda_m06l(self): td = mf_m06l.TDA() es = td.kernel(nstates=5)[0] * 27.2114 - self.assertAlmostEqual(lib.fp(es), -20.70191947889884, 4) + self.assertAlmostEqual(lib.fp(es), -20.49388623318, 4) ref = [2.74346804, 3.10082138, 6.87321246, 12.8332282, 14.30085068, 14.61913328] - self.assertAlmostEqual(abs(es[:4] - ref[:4]).max(), 0, 4) - self.assertAlmostEqual(abs(es[4] - ref[5]), 0, 4) + self.assertAlmostEqual(abs(es - ref[:5]).max(), 0, 4) def test_ab_hf(self): mf = mf_uhf @@ -354,10 +353,10 @@ def test_nto(self): mf = mf_uhf td = tdscf.TDA(mf).run() w, nto = td.get_nto(state=1) - self.assertAlmostEqual(w[0][0], 0.00018520143461015, 9) - self.assertAlmostEqual(w[1][0], 0.99963372674044326, 9) - self.assertAlmostEqual(lib.fp(w[0]), 0.00027305600430816, 9) - self.assertAlmostEqual(lib.fp(w[1]), 0.99964370569529093, 9) + self.assertAlmostEqual(w[0][0], 0.00018520143461015, 7) + self.assertAlmostEqual(w[1][0], 0.99963372674044326, 7) + self.assertAlmostEqual(lib.fp(w[0]), 0.00027305600430816, 7) + self.assertAlmostEqual(lib.fp(w[1]), 0.99964370569529093, 7) pmol = mol.copy(deep=False) pmol.symmetry = True @@ -470,13 +469,6 @@ def test_tddft_with_wfnsym(self): self.assertAlmostEqual(lib.fp(es), 0.15403661700414412, 6) td.analyze() - def test_symmetry_init_guess(self): - mol = gto.M(atom='N 0 0 0; N 0 0 1.2', basis='631g', symmetry=True, verbose=0) - td = mol.UHF.run().TDA().run(nstates=1) - self.assertAlmostEqual(td.e[0], 0.14147328219131602, 7) - # TODO: verify symmetry of td.x == A1u, close to triplet state - - if __name__ == "__main__": print("Full Tests for TD-UKS") unittest.main() diff --git a/pyscf/tdscf/uhf.py b/pyscf/tdscf/uhf.py index bed59eae6b..39ddb137f6 100644 --- a/pyscf/tdscf/uhf.py +++ b/pyscf/tdscf/uhf.py @@ -23,9 +23,9 @@ from pyscf import symm from pyscf import ao2mo from pyscf.lib import logger -from pyscf.tdscf import rhf from pyscf.scf import uhf_symm -from pyscf.scf import _response_functions +from pyscf.tdscf import rhf +from pyscf.tdscf._lr_eig import eigh as lr_eigh, eig as lr_eig from pyscf.data import nist from pyscf import __config__ @@ -63,13 +63,9 @@ def gen_tda_operation(mf, fock_ao=None, wfnsym=None): if wfnsym is not None and mol.symmetry: if isinstance(wfnsym, str): wfnsym = symm.irrep_name2id(mol.groupname, wfnsym) - orbsyma, orbsymb = uhf_symm.get_orbsym(mol, mo_coeff) wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsyma_in_d2h = numpy.asarray(orbsyma) % 10 - orbsymb_in_d2h = numpy.asarray(orbsymb) % 10 - sym_forbida = (orbsyma_in_d2h[occidxa,None] ^ orbsyma_in_d2h[viridxa]) != wfnsym - sym_forbidb = (orbsymb_in_d2h[occidxb,None] ^ orbsymb_in_d2h[viridxb]) != wfnsym - sym_forbid = numpy.hstack((sym_forbida.ravel(), sym_forbidb.ravel())) + x_sym_a, x_sym_b = _get_x_sym_table(mf) + sym_forbid = numpy.append(x_sym_a.ravel(), x_sym_b.ravel()) != wfnsym e_ia_a = mo_energy[0][viridxa] - mo_energy[0][occidxa,None] e_ia_b = mo_energy[1][viridxb] - mo_energy[1][occidxb,None] @@ -83,13 +79,14 @@ def gen_tda_operation(mf, fock_ao=None, wfnsym=None): vresp = mf.gen_response(hermi=0, max_memory=max_memory) def vind(zs): + nz = len(zs) zs = numpy.asarray(zs) if wfnsym is not None and mol.symmetry: zs = numpy.copy(zs) zs[:,sym_forbid] = 0 - za = zs[:,:nocca*nvira].reshape(-1,nocca,nvira) - zb = zs[:,nocca*nvira:].reshape(-1,noccb,nvirb) + za = zs[:,:nocca*nvira].reshape(nz,nocca,nvira) + zb = zs[:,nocca*nvira:].reshape(nz,noccb,nvirb) dmova = lib.einsum('xov,qv,po->xpq', za, orbva.conj(), orboa) dmovb = lib.einsum('xov,qv,po->xpq', zb, orbvb.conj(), orbob) @@ -100,7 +97,6 @@ def vind(zs): v1a += numpy.einsum('xia,ia->xia', za, e_ia_a) v1b += numpy.einsum('xia,ia->xia', zb, e_ia_b) - nz = zs.shape[0] hx = numpy.hstack((v1a.reshape(nz,-1), v1b.reshape(nz,-1))) if wfnsym is not None and mol.symmetry: hx[:,sym_forbid] = 0 @@ -109,6 +105,17 @@ def vind(zs): return vind, hdiag gen_tda_hop = gen_tda_operation +def _get_x_sym_table(mf): + '''Irrep (up to D2h symmetry) of each coefficient in X amplitude''' + mol = mf.mol + mo_occa, mo_occb = mf.mo_occ + orbsyma, orbsymb = uhf_symm.get_orbsym(mol, mf.mo_coeff) + orbsyma = orbsyma % 10 + orbsymb = orbsymb % 10 + x_sym_a = orbsyma[mo_occa>0,None] ^ orbsyma[mo_occa==0] + x_sym_b = orbsymb[mo_occb>0,None] ^ orbsymb[mo_occb==0] + return x_sym_a, x_sym_b + def get_ab(mf, mo_energy=None, mo_coeff=None, mo_occ=None): r'''A and B matrices for TDDFT response function. @@ -521,8 +528,8 @@ def analyze(tdobj, verbose=None): log.note('Excited State %3d: %12.5f eV %9.2f nm f=%.4f', i+1, e_ev[i], wave_length[i], f_oscillator[i]) else: - wfnsyma = rhf.analyze_wfnsym(tdobj, x_syma, x[0]) - wfnsymb = rhf.analyze_wfnsym(tdobj, x_symb, x[1]) + wfnsyma = rhf._analyze_wfnsym(tdobj, x_syma, x[0]) + wfnsymb = rhf._analyze_wfnsym(tdobj, x_symb, x[1]) if wfnsyma == wfnsymb: wfnsym = wfnsyma else: @@ -619,7 +626,7 @@ def gen_vind(self, mf=None): mf = self._scf return gen_tda_hop(mf, wfnsym=self.wfnsym) - def init_guess(self, mf, nstates=None, wfnsym=None): + def init_guess(self, mf, nstates=None, wfnsym=None, return_symmetry=False): if nstates is None: nstates = self.nstates if wfnsym is None: wfnsym = self.wfnsym @@ -630,30 +637,42 @@ def init_guess(self, mf, nstates=None, wfnsym=None): occidxb = numpy.where(mo_occ[1]>0)[0] viridxa = numpy.where(mo_occ[0]==0)[0] viridxb = numpy.where(mo_occ[1]==0)[0] - e_ia_a = (mo_energy[0][viridxa,None] - mo_energy[0][occidxa]).T - e_ia_b = (mo_energy[1][viridxb,None] - mo_energy[1][occidxb]).T - - if wfnsym is not None and mol.symmetry: - if isinstance(wfnsym, str): - wfnsym = symm.irrep_name2id(mol.groupname, wfnsym) - orbsyma, orbsymb = uhf_symm.get_orbsym(mol, mf.mo_coeff) - wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsyma_in_d2h = numpy.asarray(orbsyma) % 10 - orbsymb_in_d2h = numpy.asarray(orbsymb) % 10 - e_ia_a[(orbsyma_in_d2h[occidxa,None] ^ orbsyma_in_d2h[viridxa]) != wfnsym] = 1e99 - e_ia_b[(orbsymb_in_d2h[occidxb,None] ^ orbsymb_in_d2h[viridxb]) != wfnsym] = 1e99 - - e_ia = numpy.hstack((e_ia_a.ravel(), e_ia_b.ravel())) - nov = e_ia.size + e_ia_a = mo_energy[0][viridxa] - mo_energy[0][occidxa,None] + e_ia_b = mo_energy[1][viridxb] - mo_energy[1][occidxb,None] + nov = e_ia_a.size + e_ia_b.size nstates = min(nstates, nov) - e_threshold = numpy.sort(e_ia)[nstates-1] + + if (wfnsym is not None or return_symmetry) and mf.mol.symmetry: + x_sym_a, x_sym_b = _get_x_sym_table(mf) + if wfnsym is not None: + if isinstance(wfnsym, str): + wfnsym = symm.irrep_name2id(mol.groupname, wfnsym) + wfnsym = wfnsym % 10 # convert to D2h subgroup + e_ia_a[x_sym_a != wfnsym] = 1e99 + e_ia_b[x_sym_b != wfnsym] = 1e99 + nov_allowed = (numpy.count_nonzero(x_sym_a == wfnsym) + + numpy.count_nonzero(x_sym_b == wfnsym)) + nstates = min(nstates, nov_allowed) + + e_ia = numpy.append(e_ia_a.ravel(), e_ia_b.ravel()) + # Find the nstates-th lowest energy gap + e_threshold = numpy.partition(e_ia, nstates-1)[nstates-1] e_threshold += self.deg_eia_thresh idx = numpy.where(e_ia <= e_threshold)[0] x0 = numpy.zeros((idx.size, nov)) for i, j in enumerate(idx): x0[i, j] = 1 # Koopmans' excitations - return x0 + + if return_symmetry: + if mf.mol.symmetry: + x_sym = numpy.append(x_sym_a.ravel(), x_sym_b.ravel()) + x0sym = x_sym[idx] + else: + x0sym = None + return x0, x0sym + else: + return x0 def kernel(self, x0=None, nstates=None): '''TDA diagonalization solver @@ -666,6 +685,7 @@ def kernel(self, x0=None, nstates=None): else: self.nstates = nstates log = logger.Logger(self.stdout, self.verbose) + mol = self.mol vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) @@ -674,16 +694,19 @@ def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + x0sym = None if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) - - self.converged, self.e, x1 = \ - lib.davidson1(vind, x0, precond, - tol=self.conv_tol, - nroots=nstates, lindep=self.lindep, - max_cycle=self.max_cycle, - max_space=self.max_space, pick=pickeig, - verbose=log) + x0, x0sym = self.init_guess( + self._scf, self.nstates, return_symmetry=True) + elif mol.symmetry: + x_sym_a, x_sym_b = _get_x_sym_table(self._scf) + x_sym = numpy.append(x_sym_a.ravel(), x_sym_b.ravel()) + x0sym = [rhf._guess_wfnsym_id(self, x_sym, x) for x in x0] + + self.converged, self.e, 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, + max_memory=self.max_memory, verbose=log) nmo = self._scf.mo_occ[0].size nocca = (self._scf.mo_occ[0]>0).sum() @@ -716,7 +739,6 @@ def gen_tdhf_operation(mf, fock_ao=None, singlet=True, wfnsym=None): ''' mol = mf.mol mo_coeff = mf.mo_coeff - assert (mo_coeff[0].dtype == numpy.double) mo_energy = mf.mo_energy mo_occ = mf.mo_occ nao, nmo = mo_coeff[0].shape @@ -736,13 +758,9 @@ def gen_tdhf_operation(mf, fock_ao=None, singlet=True, wfnsym=None): if wfnsym is not None and mol.symmetry: if isinstance(wfnsym, str): wfnsym = symm.irrep_name2id(mol.groupname, wfnsym) - orbsyma, orbsymb = uhf_symm.get_orbsym(mol, mo_coeff) wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsyma_in_d2h = numpy.asarray(orbsyma) % 10 - orbsymb_in_d2h = numpy.asarray(orbsymb) % 10 - sym_forbida = (orbsyma_in_d2h[occidxa,None] ^ orbsyma_in_d2h[viridxa]) != wfnsym - sym_forbidb = (orbsymb_in_d2h[occidxb,None] ^ orbsymb_in_d2h[viridxb]) != wfnsym - sym_forbid = numpy.hstack((sym_forbida.ravel(), sym_forbidb.ravel())) + x_sym_a, x_sym_b = _get_x_sym_table(mf) + sym_forbid = numpy.append(x_sym_a.ravel(), x_sym_b.ravel()) != wfnsym e_ia_a = mo_energy[0][viridxa] - mo_energy[0][occidxa,None] e_ia_b = mo_energy[1][viridxb] - mo_energy[1][occidxb,None] @@ -793,7 +811,7 @@ def vind(xys): return vind, hdiag -class TDHF(TDA): +class TDHF(TDBase): singlet = None @@ -803,10 +821,17 @@ def gen_vind(self, mf=None): mf = self._scf return gen_tdhf_operation(mf, singlet=self.singlet, wfnsym=self.wfnsym) - def init_guess(self, mf, nstates=None, wfnsym=None): - x0 = TDA.init_guess(self, mf, nstates, wfnsym) - y0 = numpy.zeros_like(x0) - return numpy.asarray(numpy.block([[x0, y0], [y0, x0.conj()]])) + def init_guess(self, mf, nstates=None, wfnsym=None, return_symmetry=False): + if return_symmetry: + x0, x0sym = TDA.init_guess(self, mf, nstates, wfnsym, return_symmetry) + y0 = numpy.zeros_like(x0) + return numpy.hstack([x0, y0]), x0sym + else: + x0 = TDA.init_guess(self, mf, nstates, wfnsym, return_symmetry) + y0 = numpy.zeros_like(x0) + return numpy.hstack([x0, y0]) + + get_precond = rhf.TDHF.get_precond def kernel(self, x0=None, nstates=None): '''TDHF diagonalization with non-Hermitian eigenvalue solver @@ -818,29 +843,41 @@ def kernel(self, x0=None, nstates=None): nstates = self.nstates else: self.nstates = nstates + mol = self.mol log = logger.Logger(self.stdout, self.verbose) vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) + # handle single kpt PBC SCF + if getattr(self._scf, 'kpt', None) is not None: + from pyscf.pbc.lib.kpts_helper import gamma_point + real_system = (gamma_point(self._scf.kpt) and + self._scf.mo_coeff[0].dtype == numpy.double) + else: + real_system = True + # We only need positive eigenvalues def pickeig(w, v, nroots, envs): realidx = numpy.where((abs(w.imag) < REAL_EIG_THRESHOLD) & (w.real > self.positive_eig_threshold))[0] - return lib.linalg_helper._eigs_cmplx2real(w, v, realidx, - real_eigenvectors=True) + return lib.linalg_helper._eigs_cmplx2real(w, v, realidx, real_system) + x0sym = None if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) - - self.converged, w, x1 = \ - lib.davidson_nosym1(vind, x0, precond, - tol=self.conv_tol, - nroots=nstates, lindep=self.lindep, - max_cycle=self.max_cycle, - max_space=self.max_space, pick=pickeig, - verbose=log) + x0, x0sym = self.init_guess( + self._scf, self.nstates, return_symmetry=True) + elif mol.symmetry: + x_sym_a, x_sym_b = _get_x_sym_table(self._scf) + x_sym = y_sym = numpy.append(x_sym_a.ravel(), x_sym_b.ravel()) + x_sym = numpy.append(x_sym, y_sym) + x0sym = [rhf._guess_wfnsym_id(self, x_sym, x) for x in x0] + + self.converged, w, x1 = lr_eig( + vind, x0, precond, tol_residual=self.conv_tol, lindep=self.lindep, + nroots=nstates, x0sym=x0sym, pick=pickeig, max_cycle=self.max_cycle, + max_memory=self.max_memory, verbose=log) nmo = self._scf.mo_occ[0].size nocca = (self._scf.mo_occ[0]>0).sum() diff --git a/pyscf/tdscf/uks.py b/pyscf/tdscf/uks.py index 7b16a3d2d1..25c49e2ffc 100644 --- a/pyscf/tdscf/uks.py +++ b/pyscf/tdscf/uks.py @@ -21,9 +21,8 @@ from pyscf import symm from pyscf import lib from pyscf.lib import logger -from pyscf.tdscf import uhf -from pyscf.scf import uhf_symm -from pyscf.data import nist +from pyscf.tdscf import uhf, rhf +from pyscf.tdscf._lr_eig import eigh as lr_eigh from pyscf.dft.rks import KohnShamDFT from pyscf import __config__ @@ -74,13 +73,9 @@ def gen_vind(self, mf=None): if wfnsym is not None and mol.symmetry: if isinstance(wfnsym, str): wfnsym = symm.irrep_name2id(mol.groupname, wfnsym) - orbsyma, orbsymb = uhf_symm.get_orbsym(mol, mo_coeff) wfnsym = wfnsym % 10 # convert to D2h subgroup - orbsyma_in_d2h = numpy.asarray(orbsyma) % 10 - orbsymb_in_d2h = numpy.asarray(orbsymb) % 10 - sym_forbida = (orbsyma_in_d2h[occidxa,None] ^ orbsyma_in_d2h[viridxa]) != wfnsym - sym_forbidb = (orbsymb_in_d2h[occidxb,None] ^ orbsymb_in_d2h[viridxb]) != wfnsym - sym_forbid = numpy.hstack((sym_forbida.ravel(), sym_forbidb.ravel())) + x_sym_a, x_sym_b = uhf._get_x_sym_table(mf) + sym_forbid = numpy.append(x_sym_a.ravel(), x_sym_b.ravel()) != wfnsym e_ia_a = (mo_energy[0][viridxa,None] - mo_energy[0][occidxa]).T e_ia_b = (mo_energy[1][viridxb,None] - mo_energy[1][occidxb]).T @@ -123,6 +118,7 @@ def kernel(self, x0=None, nstates=None): '''TDDFT diagonalization solver ''' cpu0 = (logger.process_clock(), logger.perf_counter()) + mol = self.mol mf = self._scf if mf._numint.libxc.is_hybrid_xc(mf.xc): raise RuntimeError('%s cannot be used with hybrid functional' @@ -142,16 +138,19 @@ def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + x0sym = None if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) - - self.converged, w2, x1 = \ - lib.davidson1(vind, x0, precond, - tol=self.conv_tol, - nroots=nstates, lindep=self.lindep, - max_cycle=self.max_cycle, - max_space=self.max_space, pick=pickeig, - verbose=log) + x0, x0sym = self.init_guess( + self._scf, self.nstates, return_symmetry=True) + elif mol.symmetry: + x_sym_a, x_sym_b = uhf._get_x_sym_table(self._scf) + x_sym = numpy.append(x_sym_a.ravel(), x_sym_b.ravel()) + x0sym = [rhf._guess_wfnsym_id(self, x_sym, x) for x in x0] + + 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, + max_memory=self.max_memory, verbose=log) mo_energy = self._scf.mo_energy mo_occ = self._scf.mo_occ diff --git a/pyscf/x2c/tdscf.py b/pyscf/x2c/tdscf.py index 905b795d0f..77f6a0be46 100644 --- a/pyscf/x2c/tdscf.py +++ b/pyscf/x2c/tdscf.py @@ -290,9 +290,9 @@ def gen_vind(self, mf=None): mf = self._scf return gen_tda_hop(mf) - def init_guess(self, mf, nstates=None, wfnsym=None): + def init_guess(self, mf, nstates=None, wfnsym=None, return_symmetry=False): assert self.wfnsym is None - return ghf.TDA.init_guess(self, mf, nstates, None) + return ghf.TDA.init_guess(self, mf, nstates, None, return_symmetry) kernel = ghf.TDA.kernel @@ -313,9 +313,9 @@ def gen_vind(self, mf=None): mf = self._scf return gen_tdhf_operation(mf) - def init_guess(self, mf, nstates=None, wfnsym=None): + def init_guess(self, mf, nstates=None, wfnsym=None, return_symmetry=False): assert self.wfnsym is None - return ghf.TDHF.init_guess(self, mf, nstates, None) + return ghf.TDHF.init_guess(self, mf, nstates, None, return_symmetry) kernel = ghf.TDHF.kernel diff --git a/pyscf/x2c/test/test_tdscf.py b/pyscf/x2c/test/test_tdscf.py index ad14d4c6f4..4048219568 100644 --- a/pyscf/x2c/test/test_tdscf.py +++ b/pyscf/x2c/test/test_tdscf.py @@ -115,7 +115,7 @@ def test_mcol_mgga_ab_ks(self): mcol_m06l = dft.UKS(mol).set(xc='m06l', collinear='mcol') mcol_m06l._numint.spin_samples = 6 mcol_m06l.__dict__.update(scf.chkfile.load(mf_lda.chkfile, 'scf')) - self._check_against_ab_ks(mcol_m06l.TDDFT(), -0.6984240332038076, 2.0192987108288794) + self._check_against_ab_ks(mcol_m06l.TDDFT(), -0.6984240332038076, 2.0192987108288794, 5) def _check_against_ab_ks(self, td, refa, refb, places=6): mf = td._scf @@ -144,7 +144,7 @@ def _check_against_ab_ks(self, td, refa, refb, places=6): def test_mcol_vs_gks(self): with lib.temporary_env(lib.param, LIGHT_SPEED=20): mol = gto.M(atom='C', basis='6-31g') - ref = dft.RKS(mol) + ref = dft.UKS(mol) ref.xc = 'pbe' ref.collinear = 'mcol' ref._numint.spin_samples = 6 @@ -165,7 +165,7 @@ def test_mcol_vs_gks(self): td = mf.TDA() td.positive_eig_threshold = -10 es = td.kernel(nstates=5)[0] - self.assertAlmostEqual(abs(es - eref).max(), 0, 7) + self.assertAlmostEqual(abs(es - eref).max(), 0, 6) if __name__ == "__main__":