From 4c0998dcf22d0079b5906191e85a28f2568b15fd Mon Sep 17 00:00:00 2001 From: Hongzhou Ye Date: Tue, 29 Aug 2023 11:03:45 -0400 Subject: [PATCH] [Ready for review] Complete implementation of pbc/tdscf (#1793) * Bug fix for RSDF due to Mole basis sorting * Bug fix for pbc df Jbuild for hermi == 0 * improve TDSCF init guess degeneracy handle * complete pbc tdscf * b3lyp -> b3lyp5 * adjust values due to v2.3 * add examples; add finalize * remove non-gamma point exception * fix flake8 error * remove numer instable tests (higher roots) * update tests * fix wrong num of states * bug fix in response * keep minimal changes for pbc _response * better handle of exxdiv, including vcut * add NotImplementedError for rsjk + non-zero kshift * fix shifted kpt for TDRKS/UKS * fix flake8 error * revert CasidaTDDFT * fix flake error --------- Co-authored-by: hongzhouye <> --- examples/pbc/22-k_points_tddft.py | 126 +++-- examples/tddft/04-tddft_missing_roots.py | 46 ++ pyscf/grad/test/test_tdrhf_grad.py | 13 +- pyscf/grad/test/test_tdrks_grad.py | 21 +- pyscf/grad/test/test_tduhf_grad.py | 9 +- pyscf/grad/test/test_tduks_grad.py | 21 +- pyscf/pbc/df/df.py | 8 +- pyscf/pbc/df/df_jk.py | 504 ++++++++++++++++++- pyscf/pbc/lib/kpts_helper.py | 15 + pyscf/pbc/scf/_response_functions.py | 68 ++- pyscf/pbc/scf/addons.py | 36 ++ pyscf/pbc/tdscf/__init__.py | 7 - pyscf/pbc/tdscf/krhf.py | 247 ++++++--- pyscf/pbc/tdscf/krks.py | 33 +- pyscf/pbc/tdscf/kuhf.py | 207 +++++--- pyscf/pbc/tdscf/kuks.py | 24 +- pyscf/pbc/tdscf/rhf.py | 27 +- pyscf/pbc/tdscf/rks.py | 23 +- pyscf/pbc/tdscf/test/test_kproxy_ks.py | 2 +- pyscf/pbc/tdscf/test/test_krhf.py | 164 ++++++ pyscf/pbc/tdscf/test/test_krhf_slow.py | 2 +- pyscf/pbc/tdscf/test/test_krhf_slow_gamma.py | 4 +- pyscf/pbc/tdscf/test/test_krks.py | 146 ++++++ pyscf/pbc/tdscf/test/test_kuhf.py | 136 +++++ pyscf/pbc/tdscf/test/test_kuks.py | 134 +++++ pyscf/pbc/tdscf/test/test_rhf.py | 229 ++++++--- pyscf/pbc/tdscf/test/test_rks.py | 309 ++++++++++++ pyscf/pbc/tdscf/test/test_uhf.py | 126 +++++ pyscf/pbc/tdscf/test/test_uks.py | 230 +++++++++ pyscf/pbc/tdscf/uks.py | 2 +- pyscf/tdscf/dhf.py | 22 +- pyscf/tdscf/ghf.py | 23 +- pyscf/tdscf/gks.py | 6 +- pyscf/tdscf/rhf.py | 64 ++- pyscf/tdscf/rks.py | 6 +- pyscf/tdscf/test/test_proxy.py | 6 +- pyscf/tdscf/test/test_rhf_slow.py | 3 +- pyscf/tdscf/test/test_tdrhf.py | 19 +- pyscf/tdscf/test/test_tdrks.py | 6 +- pyscf/tdscf/uhf.py | 24 +- pyscf/tdscf/uks.py | 7 +- 41 files changed, 2673 insertions(+), 432 deletions(-) create mode 100644 examples/tddft/04-tddft_missing_roots.py create mode 100644 pyscf/pbc/tdscf/test/test_krhf.py create mode 100644 pyscf/pbc/tdscf/test/test_krks.py create mode 100644 pyscf/pbc/tdscf/test/test_kuhf.py create mode 100644 pyscf/pbc/tdscf/test/test_kuks.py create mode 100644 pyscf/pbc/tdscf/test/test_rks.py create mode 100644 pyscf/pbc/tdscf/test/test_uhf.py create mode 100644 pyscf/pbc/tdscf/test/test_uks.py diff --git a/examples/pbc/22-k_points_tddft.py b/examples/pbc/22-k_points_tddft.py index f7d3799ab5..539ffb343c 100644 --- a/examples/pbc/22-k_points_tddft.py +++ b/examples/pbc/22-k_points_tddft.py @@ -1,65 +1,91 @@ #!/usr/bin/env python +# +# Author: Hong-Zhou Ye +# -''' -TDDFT with k-point sampling or at an individual k-point +import numpy as np +from pyscf.pbc import gto, scf, tdscf +from pyscf import scf as molscf +from pyscf import lib -(This feature is in testing. We observe numerical stability problem in TDDFT -diagonalization.) ''' +TDSCF with k-point sampling +''' +atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' +a = ''' +1.7850000000 1.7850000000 0.0000000000 +0.0000000000 1.7850000000 1.7850000000 +1.7850000000 0.0000000000 1.7850000000 +''' +basis = ''' +C S + 9.031436 -1.960629e-02 + 3.821255 -1.291762e-01 + 0.473725 5.822572e-01 +C P + 4.353457 8.730943e-02 + 1.266307 2.797034e-01 + 0.398715 5.024424e-01 +''' # a trimmed DZ basis for fast test +pseudo = 'gth-hf-rev' +cell = gto.M(atom=atom, basis=basis, a=a, pseudo=pseudo).set(verbose=3) +kmesh = [2,1,1] +kpts = cell.make_kpts(kmesh) +nkpts = len(kpts) +mf = scf.KRHF(cell, kpts).rs_density_fit().run() -from pyscf.pbc import gto -from pyscf.pbc import scf -from pyscf.pbc import df -from pyscf.pbc import tdscf +log = lib.logger.new_logger(mf) -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. +''' k-point TDSCF solutions can have non-zero momentum transfer between particle and hole. + This can be controlled by `td.kshift_lst`. By default, kshift_lst = [0] and only the + zero-momentum transfer solution (i.e., 'vertical' in k-space) will be solved, as + demonstrated in the example below. ''' -cell.basis = 'gth-szv' -cell.pseudo = 'gth-pade' -cell.build() -mf = scf.KRHF(cell, cell.make_kpts([2,2,2])) -mf.run() +td = mf.TDA().set(nstates=5).run() +log.note('RHF-TDA:') +for kshift,es in zip(td.kshift_lst,td.e): + log.note('kshift = %d Eex = %s', kshift, ' '.join([f'{e:.3f}' for e in es*27.2114])) -td = tdscf.KTDA(mf) -td.nstates = 5 -td.verbose = 5 -print(td.kernel()[0] * 27.2114) +''' If GDF/RSDF is used as the density fitting method (as in this example), solutions + with non-zero particle-hole momentum-transfer solution is also possible. The example + below demonstrates how to calculate solutions with all possible kshift. -td = tdscf.KTDDFT(mf) -td.nstates = 5 -td.verbose = 5 -print(td.kernel()[0] * 27.2114) + NOTE: if FFTDF is used, pyscf will set `kshift_lst` to [0]. +''' +td = mf.TDHF().set(nstates=5, kshift_lst=list(range(nkpts))).run() +log.note('RHF-TDHF:') +for kshift,es in zip(td.kshift_lst,td.e): + log.note('kshift = %d Eex = %s', kshift, ' '.join([f'{e:.3f}' for e in es*27.2114])) -mf = scf.RHF(cell) -mf.kernel() -td = tdscf.TDA(mf) -td.kernel() -# -# Gamma-point RKS -# -ks = scf.RKS(cell) -ks.run() +''' TDHF at a single k-point compared to molecular TDSCF +''' +atom = ''' +O 0.00000 0.00000 0.11779 +H 0.00000 0.75545 -0.47116 +H 0.00000 -0.75545 -0.47116 +''' +a = np.eye(3) * 20 # big box to match isolated molecule +basis = 'sto-3g' +auxbasis = 'weigend' +cell = gto.M(atom=atom, basis=basis, a=a).set(verbose=3) +mol = cell.to_mol() -td = tdscf.KTDDFT(ks) -td.nstates = 5 -td.verbose = 5 -print(td.kernel()[0] * 27.2114) +log = lib.logger.new_logger(cell) +xc = 'b3lyp' +# pbc +mf = scf.RKS(cell).set(xc=xc).rs_density_fit(auxbasis=auxbasis).run() +pbctda = mf.TDA().run() +pbctd = mf.TDDFT().run() +# mol +molmf = molscf.RKS(cell).set(xc=xc).density_fit(auxbasis=auxbasis).run() +moltda = molmf.TDA().run() +moltd = molmf.TDDFT().run() -# TODO: -#kpt = cell.get_abs_kpts([0.25, 0.25, 0.25]) -#mf = scf.RHF(cell, kpt=kpt) -#mf.kernel() -#td = tdscf.TDA(mf) -#td.kernel() +_format = lambda e: ' '.join([f'{x*27.2114:.3f}' for x in e]) +log.note('PBC TDA : %s', _format(pbctda.e)) +log.note('Mol TDA : %s', _format(moltda.e)) +log.note('PBC TDDFT: %s', _format(pbctd.e)) +log.note('Mol TDDFT: %s', _format(moltd.e)) diff --git a/examples/tddft/04-tddft_missing_roots.py b/examples/tddft/04-tddft_missing_roots.py new file mode 100644 index 0000000000..4c0e7eb6f7 --- /dev/null +++ b/examples/tddft/04-tddft_missing_roots.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python +# +# Author: Hong-Zhou Ye +# + + +from pyscf import gto, scf, tdscf + + +''' TDSCF selected particle-hole pairs that have low energy difference (e.g., HOMO -> LUMO, + HOMO -> LUMO+1 etc.) as init guess for the iterative solution of the TDSCF equation. + In cases where there are many near degenerate p-h pairs, their contributions to the + final TDSCF solutions can all be significnat and all such pairs should be included. + + The determination of such near degeneracies in p-h pairs is controlled by the + `deg_eia_thresh` variable: if the excitation energy of two p-h pairs differs less + than `deg_eia_thresh`, they will be considered degenerate and selected simultaneously + as init guess. + + This example demonstrates the use of this variable to reduce the chance of missing + roots in TDSCF calculations. A molecule is used here for demonstration but the same + applies to periodic TDSCF calculations. +''' + +atom = ''' +O 0.00000 0.00000 0.11779 +H 0.00000 0.75545 -0.47116 +H 0.00000 -0.75545 -0.47116 +O 4.00000 0.00000 0.11779 +H 4.00000 0.75545 -0.47116 +H 4.00000 -0.75545 -0.47116 +''' +basis = 'cc-pvdz' + +mol = gto.M(atom=atom, basis=basis) +mf = scf.RHF(mol).density_fit().run() + +''' By default, `deg_eia_thresh` = 1e-3 (Ha) as can be seen in the output if verbose >= 4. +''' +mf.TDA().set(nroots=6).run() + +''' One can check if there are missing roots by using a larger value for `deg_eia_thresh`. + Note that this will increase the number of init guess vectors and increase the cost of + solving the TDSCF equation. +''' +mf.TDA().set(nroots=6, deg_eia_thresh=0.1).run() diff --git a/pyscf/grad/test/test_tdrhf_grad.py b/pyscf/grad/test/test_tdrhf_grad.py index 6da2e61006..2f1dac782a 100644 --- a/pyscf/grad/test/test_tdrhf_grad.py +++ b/pyscf/grad/test/test_tdrhf_grad.py @@ -114,7 +114,7 @@ def fvind(x): def setUpModule(): - global mol, pmol, mf + global mol, pmol, mf, nstates mol = gto.Mole() mol.verbose = 5 mol.output = '/dev/null' @@ -126,6 +126,7 @@ def setUpModule(): mol.build() pmol = mol.copy() mf = scf.RHF(mol).set(conv_tol=1e-12).run() + nstates = 5 # to ensure the first 3 TDSCF states are converged def tearDownModule(): global mol, pmol, mf @@ -134,7 +135,7 @@ def tearDownModule(): class KnownValues(unittest.TestCase): def test_tda_singlet(self): - td = tdscf.TDA(mf).run(nstates=3) + td = tdscf.TDA(mf).run(nstates=nstates) g1ref = tda_kernel(td.nuc_grad_method(), td.xy[2]) + mf.nuc_grad_method().kernel() tdg = td.nuc_grad_method().as_scanner() @@ -155,14 +156,14 @@ def test_tda_singlet(self): mf.nuc_grad_method().kernel()).max(), 0, 8) def test_tda_triplet(self): - td = tdscf.TDA(mf).run(singlet=False, nstates=3) + td = tdscf.TDA(mf).run(singlet=False, nstates=nstates) tdg = td.nuc_grad_method() # [[ 0 0 -2.81048403e-01] # [ 0 0 2.81048403e-01]] self.assertAlmostEqual(lib.fp(tdg.kernel(state=1)), 0.19667995802487931, 6) g1 = tdg.kernel(state=3) - self.assertAlmostEqual(g1[0,2], -0.47296513687621511, 8) + self.assertAlmostEqual(g1[0,2], -0.472965206465775, 6) td_solver = td.as_scanner() e1 = td_solver(pmol.set_geom_('H 0 0 1.805; F 0 0 0', unit='B')) @@ -170,7 +171,7 @@ def test_tda_triplet(self): self.assertAlmostEqual((e1[2]-e2[2])/.002, g1[0,2], 5) def test_tdhf_singlet(self): - td = tdscf.TDDFT(mf).run(nstates=3) + td = tdscf.TDDFT(mf).run(nstates=nstates) tdg = td.nuc_grad_method() # [[ 0 0 -2.71041021e-01] # [ 0 0 2.71041021e-01]] @@ -185,7 +186,7 @@ def test_tdhf_singlet(self): self.assertAlmostEqual((e1[2]-e2[2])/.002, g1[0,2], 6) def test_tdhf_triplet(self): - td = tdscf.TDDFT(mf).run(singlet=False, nstates=3) + td = tdscf.TDDFT(mf).run(singlet=False, nstates=nstates) tdg = td.nuc_grad_method() # [[ 0 0 -2.86250870e-01] # [ 0 0 2.86250870e-01]] diff --git a/pyscf/grad/test/test_tdrks_grad.py b/pyscf/grad/test/test_tdrks_grad.py index d9d88dfa94..ef0f46c4f0 100644 --- a/pyscf/grad/test/test_tdrks_grad.py +++ b/pyscf/grad/test/test_tdrks_grad.py @@ -24,7 +24,7 @@ from pyscf.grad import tdrks as tdrks_grad def setUpModule(): - global mol, mf_lda, mf_gga + global mol, mf_lda, mf_gga, nstates mol = gto.Mole() mol.verbose = 0 mol.output = None @@ -42,6 +42,7 @@ def setUpModule(): mf_gga.grids.prune = False mf_gga.conv_tol = 1e-10 mf_gga.kernel() + nstates = 5 # to ensure the first 3 TDSCF states are converged def tearDownModule(): global mol, mf_lda, mf_gga @@ -49,13 +50,13 @@ def tearDownModule(): class KnownValues(unittest.TestCase): def test_tda_singlet_lda(self): - td = tdscf.TDA(mf_lda).run(nstates=3) + td = tdscf.TDA(mf_lda).run(nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(td.xy[2]) self.assertAlmostEqual(g1[0,2], -9.23916667e-02, 6) def test_tda_triplet_lda(self): - td = tdscf.TDA(mf_lda).run(singlet=False, nstates=3) + td = tdscf.TDA(mf_lda).run(singlet=False, nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -0.3311324654, 6) @@ -67,7 +68,7 @@ def test_tda_triplet_lda(self): self.assertAlmostEqual(abs((e1[2]-e2[2])/.002 - g1[0,2]).max(), 0, 4) def test_tda_singlet_b88(self): - td = tdscf.TDA(mf_gga).run(nstates=3) + td = tdscf.TDA(mf_gga).run(nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -9.32506535e-02, 6) @@ -81,7 +82,7 @@ def test_tda_singlet_b3lyp_xcfun(self): mf.scf() td = mf.TDA() - td.nstates = 3 + td.nstates = nstates e, z = td.kernel() tdg = td.Gradients() g1 = tdg.kernel(state=3) @@ -99,7 +100,7 @@ def test_tda_triplet_b3lyp(self): mf.xc = 'b3lyp5' mf.conv_tol = 1e-14 mf.kernel() - td = tdscf.TDA(mf).run(singlet=False, nstates=3) + td = tdscf.TDA(mf).run(singlet=False, nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -0.36333834, 6) @@ -115,7 +116,7 @@ def test_tda_singlet_mgga(self): mf.xc = 'm06l' mf.conv_tol = 1e-14 mf.kernel() - td = mf.TDA().run(nstates=3) + td = mf.TDA().run(nstates=nstates) tdg = td.Gradients() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -0.1461843283, 6) @@ -128,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=3) + td = tdscf.TDDFT(mf_lda).run(nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -1.31315477e-01, 6) @@ -145,7 +146,7 @@ def test_tddft_b3lyp_high_cost(self): mf._numint.libxc = dft.xcfun mf.grids.prune = False mf.scf() - td = tdscf.TDDFT(mf).run(nstates=3) + td = tdscf.TDDFT(mf).run(nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -1.55778110e-01, 6) @@ -154,7 +155,7 @@ def test_range_separated_high_cost(self): mol = gto.M(atom="H; H 1 1.", basis='631g', verbose=0) mf = dft.RKS(mol).set(xc='CAMB3LYP') mf._numint.libxc = dft.xcfun - td = mf.apply(tdscf.TDA) + td = mf.apply(tdscf.TDA).set(nstates=nstates) tdg_scanner = td.nuc_grad_method().as_scanner().as_scanner() g = tdg_scanner(mol, state=3)[1] self.assertAlmostEqual(lib.fp(g), 0.60109310253094916, 6) diff --git a/pyscf/grad/test/test_tduhf_grad.py b/pyscf/grad/test/test_tduhf_grad.py index f0ed046625..c76db05c64 100644 --- a/pyscf/grad/test/test_tduhf_grad.py +++ b/pyscf/grad/test/test_tduhf_grad.py @@ -24,7 +24,7 @@ def setUpModule(): - global mol, pmol, mf + global mol, pmol, mf, nstates mol = gto.Mole() mol.verbose = 5 mol.output = '/dev/null' @@ -38,6 +38,7 @@ def setUpModule(): mol.build() pmol = mol.copy() mf = scf.UHF(mol).set(conv_tol=1e-12).run() + nstates = 5 # to ensure the first 3 TDSCF states are converged def tearDownModule(): global mol, pmol, mf @@ -46,7 +47,7 @@ def tearDownModule(): class KnownValues(unittest.TestCase): def test_tda(self): - td = tdscf.TDA(mf).run(nstates=3) + td = tdscf.TDA(mf).run(nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -0.78246882668628404, 6) @@ -57,10 +58,10 @@ 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=3) + td = tdscf.TDDFT(mf).run(nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(td.xy[2]) - self.assertAlmostEqual(g1[0,2], -0.78969714300299776, 6) + self.assertAlmostEqual(g1[0,2], -0.78969714300299776, 5) td_solver = td.as_scanner() e1 = td_solver(pmol.set_geom_('H 0 0 1.805; F 0 0 0', unit='B')) diff --git a/pyscf/grad/test/test_tduks_grad.py b/pyscf/grad/test/test_tduks_grad.py index 9708ab6e76..b83e8f6a9b 100644 --- a/pyscf/grad/test/test_tduks_grad.py +++ b/pyscf/grad/test/test_tduks_grad.py @@ -24,7 +24,7 @@ from pyscf.grad import tduks as tduks_grad def setUpModule(): - global mol, pmol, mf_lda, mf_gga + global mol, pmol, mf_lda, mf_gga, nstates mol = gto.Mole() mol.verbose = 5 mol.output = '/dev/null' @@ -41,6 +41,7 @@ def setUpModule(): mf_lda.kernel() mf_gga = dft.UKS(mol).set(xc='b88,', conv_tol=1e-12) mf_gga.kernel() + nstates = 5 # to ensure the first 3 TDSCF states are converged def tearDownModule(): global mol, pmol, mf_lda, mf_gga @@ -49,10 +50,10 @@ def tearDownModule(): class KnownValues(unittest.TestCase): def test_tda_lda(self): - td = tdscf.TDA(mf_lda).run(nstates=3) + td = tdscf.TDA(mf_lda).run(nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(td.xy[2]) - self.assertAlmostEqual(g1[0,2], -0.40279473514282405, 6) + self.assertAlmostEqual(g1[0,2], -0.7944872119457362, 6) td_solver = td.as_scanner() e1 = td_solver(pmol.set_geom_('H 0 0 1.805; F 0 0 0', unit='B')) @@ -60,7 +61,7 @@ def test_tda_lda(self): self.assertAlmostEqual((e1[2]-e2[2])/.002, g1[0,2], 4) def test_tda_b88(self): - td = tdscf.TDA(mf_gga).run(nstates=3) + td = tdscf.TDA(mf_gga).run(nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -0.8120037135120326, 6) @@ -71,10 +72,10 @@ 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=3) + td = tdscf.TDDFT(mf_lda).run(nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) - self.assertAlmostEqual(g1[0,2], -0.39791714992157035, 6) + self.assertAlmostEqual(g1[0,2], -0.800487816830773, 6) td_solver = td.as_scanner() e1 = td_solver(pmol.set_geom_('H 0 0 1.805; F 0 0 0', unit='B')) @@ -87,7 +88,7 @@ def test_tda_mgga(self): mf.xc = 'm06l' mf.conv_tol = 1e-12 mf.kernel() - td = mf.TDA().run(nstates=3) + td = mf.TDA().run(nstates=nstates) tdg = td.Gradients() g1 = tdg.kernel(state=2) self.assertAlmostEqual(g1[0,2], -0.31324464083043635, 4) @@ -103,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=3) + td = tdscf.TDDFT(mf).run(nstates=nstates) tdg = td.nuc_grad_method() g1 = tdg.kernel(state=3) self.assertAlmostEqual(g1[0,2], -0.80446691153291727, 6) @@ -116,7 +117,7 @@ def test_tddft_b3lyp(self): def test_range_separated(self): mol = gto.M(atom="H; H 1 1.", basis='631g', verbose=0) mf = dft.UKS(mol).set(xc='CAMB3LYP') - td = mf.apply(tdscf.TDA) + td = mf.apply(tdscf.TDA).set(nstates=nstates) tdg_scanner = td.nuc_grad_method().as_scanner() g = tdg_scanner(mol, state=3)[1] self.assertAlmostEqual(lib.fp(g), -0.46656653988919661, 6) @@ -142,7 +143,7 @@ def test_custom_xc(self): mf.kernel() td = mf.TDA() - td.nstates = 3 + td.nstates = nstates e, z = td.kernel() tdg = td.Gradients() g1 = tdg.kernel(state=3) diff --git a/pyscf/pbc/df/df.py b/pyscf/pbc/df/df.py index 303bb42984..537aeb9f5c 100644 --- a/pyscf/pbc/df/df.py +++ b/pyscf/pbc/df/df.py @@ -308,7 +308,7 @@ def has_kpts(self, kpts): len(member(kpt, self.kpts_band))>0) for kpt in kpts) def sr_loop(self, kpti_kptj=numpy.zeros((2,3)), max_memory=2000, - compact=True, blksize=None): + compact=True, blksize=None, aux_slice=None): '''Short range part''' if self._cderi is None: self.build() @@ -350,7 +350,8 @@ def load(aux_slice): return LpqR, LpqI with _load3c(self._cderi, self._dataname, kpti_kptj) as j3c: - slices = lib.prange(0, j3c.shape[0], blksize) + if aux_slice is None: aux_slice = (0, j3c.shape[0]) + slices = lib.prange(*aux_slice, blksize) for LpqR, LpqI in lib.map_with_prefetch(load, slices): yield LpqR, LpqI, 1 LpqR = LpqI = None @@ -360,7 +361,8 @@ def load(aux_slice): # CDERI tensor of negative part. with _load3c(self._cderi, self._dataname+'-', kpti_kptj, ignore_key_error=True) as j3c: - slices = lib.prange(0, j3c.shape[0], blksize) + if aux_slice is None: aux_slice = (0, j3c.shape[0]) + slices = lib.prange(*aux_slice, blksize) for LpqR, LpqI in lib.map_with_prefetch(load, slices): yield LpqR, LpqI, -1 LpqR = LpqI = None diff --git a/pyscf/pbc/df/df_jk.py b/pyscf/pbc/df/df_jk.py index 84d95bec3a..0f60943d2e 100644 --- a/pyscf/pbc/df/df_jk.py +++ b/pyscf/pbc/df/df_jk.py @@ -30,7 +30,7 @@ from pyscf.lib import logger, zdotNN, zdotCN, zdotNC from pyscf.pbc import tools from pyscf.pbc.lib.kpts import KPoints -from pyscf.pbc.lib.kpts_helper import is_zero, gamma_point, member +from pyscf.pbc.lib.kpts_helper import is_zero, gamma_point, member, get_kconserv_ria from pyscf import __config__ DM2MO_PREC = getattr(__config__, 'pbc_gto_df_df_jk_dm2mo_prec', 1e-10) @@ -159,6 +159,114 @@ def get_j_kpts(mydf, dm_kpts, hermi=1, kpts=numpy.zeros((1,3)), kpts_band=None): return _format_jks(vj_kpts, dm_kpts, input_band, kpts) +def get_j_kpts_kshift(mydf, dm_kpts, kshift, hermi=0, kpts=numpy.zeros((1,3)), kpts_band=None): + r''' Math: + J^{k1 k1'}_{pq} + = (1/Nk) \sum_{k2} \sum_{rs} (p k1 q k1' |r k2' s k2) D_{sr}^{k2 k2'} + where k1' and k2' satisfies + (k1 - k1' - kpts[kshift]) \dot a = 2n \pi + (k2 - k2' - kpts[kshift]) \dot a = 2n \pi + For kshift = 0, :func:`get_j_kpts` is called. + ''' + if kshift == 0: + return get_j_kpts(mydf, dm_kpts, hermi=hermi, kpts=kpts, kpts_band=kpts_band) + + if kpts_band is not None: + raise NotImplementedError + + log = logger.Logger(mydf.stdout, mydf.verbose) + t0 = (logger.process_clock(), logger.perf_counter()) + if mydf._cderi is None or not mydf.has_kpts(kpts_band): + if mydf._cderi is not None: + log.warn('DF integrals for band k-points were not found %s. ' + 'DF integrals will be rebuilt to include band k-points.', + mydf._cderi) + mydf.build(kpts_band=kpts_band) + t0 = log.timer_debug1('Init get_j_kpts', *t0) + + dm_kpts = lib.asarray(dm_kpts, order='C') + dms = _format_dms(dm_kpts, kpts) + nset, nkpts, nao = dms.shape[:3] + if mydf.auxcell is None: + # If mydf._cderi is the file that generated from another calculation, + # guess naux based on the contents of the integral file. + naux = mydf.get_naoaux() + else: + naux = mydf.auxcell.nao_nr() + nao_pair = nao * (nao+1) // 2 + + 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) + + kconserv = get_kconserv_ria(mydf.cell, kpts)[kshift] + + t1 = (logger.process_clock(), logger.perf_counter()) + dmsR = dms.real.transpose(0,1,3,2).reshape(nset,nkpts,nao**2) + dmsI = dms.imag.transpose(0,1,3,2).reshape(nset,nkpts,nao**2) + rhoR = numpy.zeros((nset,naux)) + rhoI = numpy.zeros((nset,naux)) + max_memory = max(2000, (mydf.max_memory - lib.current_memory()[0])) + for k, kpt in enumerate(kpts): + kp = kconserv[k] + kptp = kpts[kp] + kptii = numpy.asarray((kptp,kpt)) + p1 = 0 + for LpqR, LpqI, sign in mydf.sr_loop(kptii, max_memory, False): + p0, p1 = p1, p1+LpqR.shape[0] + #:Lpq = (LpqR + LpqI*1j).reshape(-1,nao,nao) + #:rhoR[:,p0:p1] += numpy.einsum('Lpq,xqp->xL', Lpq, dms[:,k]).real + #:rhoI[:,p0:p1] += numpy.einsum('Lpq,xqp->xL', Lpq, dms[:,k]).imag + rhoR[:,p0:p1] += sign * numpy.einsum('Lp,xp->xL', LpqR, dmsR[:,k]) + rhoI[:,p0:p1] += sign * numpy.einsum('Lp,xp->xL', LpqR, dmsI[:,k]) + if LpqI is not None: + rhoR[:,p0:p1] -= sign * numpy.einsum('Lp,xp->xL', LpqI, dmsI[:,k]) + rhoI[:,p0:p1] += sign * numpy.einsum('Lp,xp->xL', LpqI, dmsR[:,k]) + LpqR = LpqI = None + t1 = log.timer_debug1('get_j pass 1', *t1) + + weight = 1./nkpts + rhoR *= weight + rhoI *= weight + if hermi == 0: + aos2symm = False + vjR = numpy.zeros((nset,nband,nao**2)) + vjI = numpy.zeros((nset,nband,nao**2)) + else: + aos2symm = True + vjR = numpy.zeros((nset,nband,nao_pair)) + vjI = numpy.zeros((nset,nband,nao_pair)) + for k, kpt in enumerate(kpts_band): + kp = kconserv[k] + kptp = kpts[kp] + kptii = numpy.asarray((kpt,kptp)) + p1 = 0 + for LpqR, LpqI, sign in mydf.sr_loop(kptii, max_memory, aos2symm): + p0, p1 = p1, p1+LpqR.shape[0] + #:Lpq = (LpqR + LpqI*1j)#.reshape(-1,nao,nao) + #:vjR[:,k] += numpy.dot(rho[:,p0:p1], Lpq).real + #:vjI[:,k] += numpy.dot(rho[:,p0:p1], Lpq).imag + vjR[:,k] += numpy.dot(rhoR[:,p0:p1], LpqR) + if not j_real: + vjI[:,k] += numpy.dot(rhoI[:,p0:p1], LpqR) + if LpqI is not None: + vjR[:,k] -= numpy.dot(rhoI[:,p0:p1], LpqI) + vjI[:,k] += numpy.dot(rhoR[:,p0:p1], LpqI) + LpqR = LpqI = None + t1 = log.timer_debug1('get_j pass 2', *t1) + + if j_real: + vj_kpts = vjR + else: + vj_kpts = vjR + vjI*1j + if aos2symm: + vj_kpts = lib.unpack_tril(vj_kpts.reshape(-1,nao_pair)) + vj_kpts = vj_kpts.reshape(nset,nband,nao,nao) + + log.timer('get_j', *t0) + + return _format_jks(vj_kpts, dm_kpts, input_band, kpts) + def get_k_kpts(mydf, dm_kpts, hermi=1, kpts=numpy.zeros((1,3)), kpts_band=None, exxdiv=None): cell = mydf.cell @@ -552,6 +660,400 @@ def make_kpt(ki, kj, swap_2e, inverse_idx=None): return _format_jks(vk_kpts, dm_kpts, input_band, kpts) +def get_k_kpts_kshift(mydf, dm_kpts, kshift, hermi=0, kpts=numpy.zeros((1,3)), kpts_band=None, + exxdiv=None): + r''' Math: + K^{k1 k1'}_{pq} + = (1/Nk) \sum_{k2} \sum_{rs} (p k1 s k2 | r k2' q k1') D_{sr}^{k2 k2'} + where k1' and k2' satisfies + (k1 - k1' - kpts[kshift]) \dot a = 2n \pi + (k2 - k2' - kpts[kshift]) \dot a = 2n \pi + For kshift = 0, :func:`get_k_kpts` is called. + ''' + if kshift == 0: + return get_k_kpts(mydf, dm_kpts, hermi=hermi, kpts=kpts, kpts_band=kpts_band, + exxdiv=exxdiv) + + if kpts_band is not None: + raise NotImplementedError + + cell = mydf.cell + log = logger.Logger(mydf.stdout, mydf.verbose) + + if exxdiv is not None and exxdiv != 'ewald': + log.warn('GDF does not support exxdiv %s. ' + 'exxdiv needs to be "ewald" or None', exxdiv) + raise RuntimeError('GDF does not support exxdiv %s' % exxdiv) + + t0 = (logger.process_clock(), logger.perf_counter()) + if mydf._cderi is None or not mydf.has_kpts(kpts_band): + if mydf._cderi is not None: + log.warn('DF integrals for band k-points were not found %s. ' + 'DF integrals will be rebuilt to include band k-points.', + mydf._cderi) + mydf.build(kpts_band=kpts_band) + t0 = log.timer_debug1('Init get_k_kpts', *t0) + + mo_coeff = getattr(dm_kpts, 'mo_coeff', None) + if mo_coeff is not None: + mo_occ = dm_kpts.mo_occ + + dm_kpts = lib.asarray(dm_kpts, order='C') + dms = _format_dms(dm_kpts, kpts) + nset, nkpts, nao = dms.shape[:3] + + skmoR = skmo2R = None + if not mydf.force_dm_kbuild: + if mo_coeff is not None: + if isinstance(mo_coeff[0], (list, tuple)): + mo_coeff = [mo for mo1 in mo_coeff for mo in mo1] + if len(mo_coeff) != nset*nkpts: # wrong shape + log.warn('mo_coeff from dm tag has wrong shape. ' + 'Calculating mo from dm instead.') + mo_coeff = None + elif isinstance(mo_occ[0], (list, tuple)): + mo_occ = [mo for mo1 in mo_occ for mo in mo1] + if mo_coeff is not None: + skmoR, skmoI = _format_mo(mo_coeff, mo_occ, shape=(nset,nkpts), order='F', + precision=cell.precision) + elif hermi == 1: + skmoR, skmoI = _mo_from_dm(dms.reshape(-1,nao,nao), method='eigh', + shape=(nset,nkpts), order='F', + precision=cell.precision) + if skmoR is None: + log.debug1('get_k_kpts: Eigh fails for input dm due to non-PSD. ' + 'Try SVD instead.') + # No symmetry can be explored for shifted K build; manually make it asymmetric here. + skmo2R = skmoR + skmo2I = skmoI + if skmoR is None: + skmoR, skmoI, skmo2R, skmo2I = _mo_from_dm(dms.reshape(-1,nao,nao), + method='svd', shape=(nset,nkpts), + order='F', precision=cell.precision) + if skmoR[0,0].shape[1] > nao//2: + log.debug1('get_k_kpts: rank(dm) = %d exceeds half of nao = %d. ' + 'Fall back to DM-based build.', skmoR[0,0].shape[1], nao) + skmoR = skmo2R = None + + kpts_band, input_band = _format_kpts_band(kpts_band, kpts), kpts_band + nband = len(kpts_band) + vkR = numpy.zeros((nset,nband,nao,nao)) + vkI = numpy.zeros((nset,nband,nao,nao)) + + kconserv = get_kconserv_ria(cell, kpts)[kshift] + + tspans = numpy.zeros((7,2)) + tspannames = ['buf1', 'ct11', 'ct12', 'buf2', 'ct21', 'ct22', 'load'] + + ''' math + K(p,q; k2 from k1) + = V(r k1', q k2', p k2, s k1) * D(s,r; k1 k1') + = V(L, r k1', q k2') * V(L, s k1, p k2).conj() * D(s,r; k1 k1') eqn (1) + --> in terms D's low-rank form: D^{k k'} = dot( A^k, B^k'.T.conj() ) + = ( V(L, s k1, p k2) * A(s,i; k1).conj() ).conj() + * V(L, r k1', q k2') * B(r,i; k1').conj() eqn (2) + = X(L, i k1, p k2).conj() * Y(L, i k1, q k2') eqn (3) + + if swap_2e: + K(p,q; k1 from k2) + = V(p k1, s k2, r k2', q k1') * D(s,r; k2 k2') + = V(L, p k1, s k2) * V(L, q k1', r k2').conj() * D(s,r; k2 k2') eqn (1') + --> in terms D's low-rank form: D^{k k'} = dot( A^k, B^k'.T.conj() ) + = V(L, p k1, s k2) * A(s,i; k2) + * ( V(L, q k1', r k2') * B(r,i; k2') ).conj() eqn (2') + = X(L, p k1, i k2) * Y(L, q k1, i k2).conj() eqn (3') + + Mode 1: DM-based K-build uses eqn (1) and eqn (1') + Mode 3: Asymm MO-based K-build uses eqns (2,3) and eqns (2',3') + ''' + # K_pq = ( p{k1} i{k2} | i{k2} q{k1} ) + if skmoR is None: # input dm is not Hermitian/PSD --> build K from dm + log.debug2('get_k_kpts: build K from dm') + dmsR = numpy.asarray(dms.real, order='C') + dmsI = numpy.asarray(dms.imag, order='C') + bufR = numpy.empty((mydf.blockdim*nao**2)) + bufI = numpy.empty((mydf.blockdim*nao**2)) + bufR1 = numpy.empty((mydf.blockdim*nao**2)) + bufI1 = numpy.empty((mydf.blockdim*nao**2)) + max_memory = max(2000, mydf.max_memory-lib.current_memory()[0]) + def make_kpt(ki, kj, swap_2e, inverse_idx=None): + if inverse_idx: + raise NotImplementedError + + kpti = kpts[ki] + kptj = kpts_band[kj] + kip = kconserv[ki] + kjp = kconserv[kj] + kptip = kpts[kip] + kptjp = kpts[kjp] + + ''' + K(p,q; k2 from k1) + = V(r k1', q k2', p k2, s k1) * D(s,r; k1 k1') + = (r k1' | L | q k2') (s k1 | L | p k2).conj() D(s k1 | r k1') + = [ D[k1](s | r) A(r | L | q) ] B(s | L | p).conj() + K(p,q; k1 from k2) + = V(p k1, s k2, r k2', q k1') * D(s,r; k2 k2') + = [ B(p | L | s) D[k2](s | r) ] A(q | L | r).conj() + ''' + tock = numpy.asarray((logger.process_clock(), logger.perf_counter())) + + p1 = 0 + for LpqR, LpqI, sign in mydf.sr_loop((kptip,kptjp), max_memory*0.5, compact=False): + nrow = LpqR.shape[0] + p0, p1 = p1, p1 + nrow + + tick = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[6] += tick - tock + + pLqR = numpy.ndarray((nao,nrow,nao), buffer=bufR) + pLqI = numpy.ndarray((nao,nrow,nao), buffer=bufI) + tmpR = numpy.ndarray((nao,nrow*nao), buffer=LpqR) + tmpI = numpy.ndarray((nao,nrow*nao), buffer=LpqI) + pLqR[:] = LpqR.reshape(-1,nao,nao).transpose(1,0,2) + pLqI[:] = LpqI.reshape(-1,nao,nao).transpose(1,0,2) + + tock = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[0] += tock - tick + + for LpqR1, LpqI1, sign1 in mydf.sr_loop((kpti,kptj), blksize=nrow, compact=False, + aux_slice=(p0,p1)): + nrow1 = LpqR1.shape[0] + pLqR1 = numpy.ndarray((nao,nrow1,nao), buffer=bufR1) + pLqI1 = numpy.ndarray((nao,nrow1,nao), buffer=bufI1) + pLqR1[:] = LpqR1.reshape(-1,nao,nao).transpose(1,0,2) + pLqI1[:] = LpqI1.reshape(-1,nao,nao).transpose(1,0,2) + LpqR1 = LpqI1 = None + + for i in range(nset): + zdotNN(dmsR[i,ki], dmsI[i,ki], pLqR.reshape(nao,-1), + pLqI.reshape(nao,-1), 1, tmpR, tmpI) + tick = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[1] += tick - tock + zdotCN(pLqR1.reshape(-1,nao).T, pLqI1.reshape(-1,nao).T, + tmpR.reshape(-1,nao), tmpI.reshape(-1,nao), + sign, vkR[i,kj], vkI[i,kj], 1) + tock = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[2] += tock - tick + + if swap_2e: + tmpR = tmpR.reshape(nao*nrow1,nao) + tmpI = tmpI.reshape(nao*nrow1,nao) + tick = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[3] += tick - tock + ki_tmp = ki + kj_tmp = kj + if inverse_idx: + ki_tmp = inverse_idx[0] + kj_tmp = inverse_idx[1] + for i in range(nset): + zdotNN(pLqR1.reshape(-1,nao), pLqI1.reshape(-1,nao), + dmsR[i,kj_tmp], dmsI[i,kj_tmp], 1, tmpR, tmpI) + tock = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[4] += tock - tick + zdotNC(tmpR.reshape(nao,-1), tmpI.reshape(nao,-1), + pLqR.reshape(nao,-1).T, pLqI.reshape(nao,-1).T, + sign, vkR[i,ki_tmp], vkI[i,ki_tmp], 1) + tick = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[5] += tick - tock + + tock = numpy.asarray((logger.process_clock(), logger.perf_counter())) + + LpqR = LpqI = LpqR1 = LpqI1 = pLqR = pLqI = pLqR1 = pLqI1 = tmpR = tmpI = None + else: + log.debug2('get_k_kpts: build K from mo coeff') + skmo1R = skmoR + skmo1I = skmoI + nmo = skmoR[0,0].shape[1] + log.debug2('get_k_kpts: rank(dm) = %d / %d', nmo, nao) + skmoI_mask = numpy.asarray([[max(abs(skmo1I[i,k]).max(), + abs(skmo2I[i,k]).max()) > cell.precision + for k in range(nkpts)] for i in range(nset)]) + bufR = numpy.empty((mydf.blockdim*nao**2)) + bufI = numpy.empty((mydf.blockdim*nao**2)) + bufR1 = numpy.empty((mydf.blockdim*nao**2)) + bufI1 = numpy.empty((mydf.blockdim*nao**2)) + max_memory = max(2000, mydf.max_memory-lib.current_memory()[0]) + def make_kpt(ki, kj, swap_2e, inverse_idx=None): + if inverse_idx: + raise NotImplementedError + + kip = kconserv[ki] + kjp = kconserv[kj] + kpti = kpts[ki] + kptj = kpts_band[kj] + kptip = kpts[kip] + kptjp = kpts_band[kjp] + + ''' + K(p,q; k2 from k1) + = V(r k1', q k2', p k2, s k1) * D(s,r; k1 k1') + = A(r | L | q) B(s | L | p).conj() x(s k1 | i) y(r k1' | i).conj() + = [ x(s | i) B(s | L | p).conj() ] * + [ y(r | i) A(r | L | q).conj() ].conj() + K(p,q; k1 from k2) + = V(p k1, s k2, r k2', q k1') * D(s,r; k2 k2') + = B(p | L | s) A(q | L | r).conj() x(s k2 | i) y(r k2' | i).conj() + = [ B(p | L | s) x(s | i) ] * + [ A(q | L | r) y(s | i) ].conj() + ''' + + tock = numpy.asarray((logger.process_clock(), logger.perf_counter())) + + p1 = 0 + for LpqR, LpqI, sign in mydf.sr_loop((kptip,kptjp), max_memory*0.5, compact=False): + nrow = LpqR.shape[0] + p0, p1 = p1, p1 + nrow + + tick = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[6] += tick - tock + + pLqR = numpy.ndarray((nao,nrow,nao), buffer=bufR) + pLqI = numpy.ndarray((nao,nrow,nao), buffer=bufI) + tmp1R = numpy.ndarray((nmo,nrow*nao), buffer=LpqR) + tmp1I = numpy.ndarray((nmo,nrow*nao), buffer=LpqI) + tmp2R = numpy.ndarray((nmo,nrow*nao), buffer=LpqR.reshape(-1)[tmp1R.size:]) + tmp2I = numpy.ndarray((nmo,nrow*nao), buffer=LpqI.reshape(-1)[tmp1I.size:]) + pLqR[:] = LpqR.reshape(-1,nao,nao).transpose(1,0,2) + pLqI[:] = LpqI.reshape(-1,nao,nao).transpose(1,0,2) + + tock = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[0] += tock - tick + + for LpqR1, LpqI1, sign1 in mydf.sr_loop((kpti,kptj), blksize=nrow, compact=False, + aux_slice=(p0,p1)): + nrow1 = LpqR1.shape[0] + pLqR1 = numpy.ndarray((nao,nrow1,nao), buffer=bufR1) + pLqI1 = numpy.ndarray((nao,nrow1,nao), buffer=bufI1) + pLqR1[:] = LpqR1.reshape(-1,nao,nao).transpose(1,0,2) + pLqI1[:] = LpqI1.reshape(-1,nao,nao).transpose(1,0,2) + LpqR1 = LpqI1 = None + + ''' + = [ x(s | i) B(s | L | p).conj() ] * + [ y(r | i) A(r | L | q).conj() ].conj() + ''' + for i in range(nset): + mo1R = skmo1R[i,ki] + mo2R = skmo2R[i,ki] + if skmoI_mask[i,ki]: + mo1I = skmo1I[i,ki] + mo2I = skmo2I[i,ki] + zdotNC(mo1R.T, mo1I.T, pLqR1.reshape(nao,-1), pLqI1.reshape(nao,-1), + 1, tmp1R, tmp1I) + zdotNC(mo2R.T, mo2I.T, pLqR.reshape(nao,-1), pLqI.reshape(nao,-1), + 1, tmp2R, tmp2I) + else: + lib.ddot(mo1R.T, pLqR1.reshape(nao,-1), 1, tmp1R) + lib.ddot(mo1R.T, pLqI1.reshape(nao,-1), 1, tmp1I) + lib.ddot(mo2R.T, pLqR.reshape(nao,-1), 1, tmp2R) + lib.ddot(mo2R.T, pLqI.reshape(nao,-1), 1, tmp2I) + tick = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[1] += tick - tock + zdotNC(tmp1R.reshape(-1,nao).T, tmp1I.reshape(-1,nao).T, + tmp2R.reshape(-1,nao), tmp2I.reshape(-1,nao), + sign, vkR[i,kj], vkI[i,kj], 1) + tock = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[2] += tock - tick + + if swap_2e: + tmp1R = tmp1R.reshape(nrow*nao,nmo) + tmp1I = tmp1I.reshape(nrow*nao,nmo) + tmp2R = tmp2R.reshape(nrow*nao,nmo) + tmp2I = tmp2I.reshape(nrow*nao,nmo) + tick = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[3] += tick - tock + ki_tmp = ki + kj_tmp = kj + if inverse_idx: + ki_tmp = inverse_idx[0] + kj_tmp = inverse_idx[1] + ''' + = [ B(p | L | s) x(s | i) ] * + [ A(q | L | r) y(s | i) ].conj() + ''' + for i in range(nset): + mo1R = skmo1R[i,kj_tmp] + mo2R = skmo2R[i,kj_tmp] + if skmoI_mask[i,kj_tmp]: + mo1I = skmo1I[i,kj_tmp] + mo2I = skmo2I[i,kj_tmp] + zdotNN(pLqR1.reshape(-1,nao), pLqI1.reshape(-1,nao), mo1R, mo1I, + 1, tmp1R, tmp1I) + zdotNN(pLqR.reshape(-1,nao), pLqI.reshape(-1,nao), mo2R, mo2I, + 1, tmp2R, tmp2I) + else: + lib.ddot(pLqR1.reshape(-1,nao), mo1R, 1, tmp1R) + lib.ddot(pLqI1.reshape(-1,nao), mo1R, 1, tmp1I) + lib.ddot(pLqR.reshape(-1,nao), mo2R, 1, tmp2R) + lib.ddot(pLqI.reshape(-1,nao), mo2R, 1, tmp2I) + tock = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[4] += tock - tick + zdotNC(tmp1R.reshape(nao,-1), tmp1I.reshape(nao,-1), + tmp2R.reshape(nao,-1).T, tmp2I.reshape(nao,-1).T, + sign, vkR[i,ki_tmp], vkI[i,ki_tmp], 1) + tick = numpy.asarray((logger.process_clock(), logger.perf_counter())) + tspans[5] += tick - tock + + tock = numpy.asarray((logger.process_clock(), logger.perf_counter())) + + LpqR = LpqI = pLqR = pLqI = tmp1R = tmp1I = tmp2R = tmp2I = None + + t1 = (logger.process_clock(), logger.perf_counter()) + if kpts_band is kpts: # normal k-points HF/DFT + for ki in range(nkpts): + for kj in range(ki): + make_kpt(ki, kj, True) + make_kpt(ki, ki, False) + t1 = log.timer_debug1('get_k_kpts: make_kpt ki>=kj (%d,*)'%ki, *t1) + else: + idx_in_kpts = [] + for kpt in kpts_band: + idx = member(kpt, kpts) + if len(idx) > 0: + idx_in_kpts.append(idx[0]) + else: + idx_in_kpts.append(-1) + idx_in_kpts_band = [] + for kpt in kpts: + idx = member(kpt, kpts_band) + if len(idx) > 0: + idx_in_kpts_band.append(idx[0]) + else: + idx_in_kpts_band.append(-1) + + for ki in range(nkpts): + for kj in range(nband): + if idx_in_kpts[kj] == -1 or idx_in_kpts[kj] == ki: + make_kpt(ki, kj, False) + elif idx_in_kpts[kj] < ki: + if idx_in_kpts_band[ki] == -1: + make_kpt(ki, kj, False) + else: + make_kpt(ki, kj, True, (idx_in_kpts_band[ki], idx_in_kpts[kj])) + else: + if idx_in_kpts_band[ki] == -1: + make_kpt(ki, kj, False) + t1 = log.timer_debug1('get_k_kpts: make_kpt (%d,*)'%ki, *t1) + + for tspan, tspanname in zip(tspans,tspannames): + log.debug1(' CPU time for %s %10.2f sec, wall time %10.2f sec', + tspanname, *tspan) + + if (gamma_point(kpts) and gamma_point(kpts_band) and + not numpy.iscomplexobj(dm_kpts)): + vk_kpts = vkR + else: + vk_kpts = vkR + vkI * 1j + vk_kpts *= 1./nkpts + + if exxdiv == 'ewald': + _ewald_exxdiv_for_G0(cell, kpts, dms, vk_kpts, kpts_band) + + log.timer('get_k_kpts', *t0) + + return _format_jks(vk_kpts, dm_kpts, input_band, kpts) + ################################################## # diff --git a/pyscf/pbc/lib/kpts_helper.py b/pyscf/pbc/lib/kpts_helper.py index 9188ec3008..0f393986dd 100644 --- a/pyscf/pbc/lib/kpts_helper.py +++ b/pyscf/pbc/lib/kpts_helper.py @@ -282,6 +282,21 @@ def get_kconserv(cell, kpts): kconserv[mask] = N return kconserv +def get_kconserv_ria(cell, kpts): + r''' Get the momentum conservation array for single excitation amplitudes + for a set of k-points with appropriate k-shift. + + Given k-point indices m the array kconserv[kshift,m] returns the index n that + satifies a shifted momentum conservation, + + (k(m) - k(n) - k(kshift)) \dot a = 2n\pi + + This is used in CIS/EOM-CCSD where the single excitation amplitudes + r_{i k(m)}^{a k(n)} + are zero unless n satisfies the above. + ''' + # TODO: implement it more efficiently + return get_kconserv(cell, kpts)[:,:,0].T def check_kpt_antiperm_symmetry(array, idx1, idx2, tolerance=1e-8): '''Checks antipermutational symmetry for k-point array. diff --git a/pyscf/pbc/scf/_response_functions.py b/pyscf/pbc/scf/_response_functions.py index 2926d72753..7869f0b29a 100644 --- a/pyscf/pbc/scf/_response_functions.py +++ b/pyscf/pbc/scf/_response_functions.py @@ -59,7 +59,7 @@ def _gen_rhf_response(mf, mo_coeff=None, mo_occ=None, max_memory = max(2000, mf.max_memory*.8-mem_now) if singlet is None: # Without specify singlet, general case - def vind(dm1): + def vind(dm1, kshift=0): # The singlet hessian if hermi == 2: v1 = numpy.zeros_like(dm1) @@ -68,16 +68,16 @@ def vind(dm1): rho0, vxc, fxc, kpts, max_memory=max_memory) if hybrid: if hermi != 2: - vj, vk = mf.get_jk(cell, dm1, hermi=hermi, kpts=kpts) + vj, vk = _get_jk(mf, cell, dm1, hermi, kpts, kshift) v1 += vj - .5 * hyb * vk else: - v1 -= .5 * hyb * mf.get_k(cell, dm1, hermi=hermi, kpts=kpts) + v1 -= .5 * hyb * _get_k(mf, cell, dm1, hermi, kpts, kshift) elif hermi != 2: - v1 += mf.get_j(cell, dm1, hermi=hermi, kpts=kpts) + v1 += _get_j(mf, cell, dm1, hermi, kpts, kshift) return v1 elif singlet: - def vind(dm1): + def vind(dm1, kshift=0): if hermi == 2: v1 = numpy.zeros_like(dm1) else: @@ -88,15 +88,15 @@ def vind(dm1): v1 *= .5 if hybrid: if hermi != 2: - vj, vk = mf.get_jk(cell, dm1, hermi=hermi, kpts=kpts) + vj, vk = _get_jk(mf, cell, dm1, hermi, kpts, kshift) v1 += vj - .5 * hyb * vk else: - v1 -= .5 * hyb * mf.get_k(cell, dm1, hermi=hermi, kpts=kpts) + v1 -= .5 * hyb * _get_k(mf, cell, dm1, hermi, kpts, kshift) elif hermi != 2: - v1 += mf.get_j(cell, dm1, hermi=hermi, kpts=kpts) + v1 += _get_j(mf, cell, dm1, hermi, kpts, kshift) return v1 else: # triplet - def vind(dm1): + def vind(dm1, kshift=0): if hermi == 2: v1 = numpy.zeros_like(dm1) else: @@ -106,17 +106,17 @@ def vind(dm1): max_memory=max_memory) v1 *= .5 if hybrid: - v1 += -.5 * hyb * mf.get_k(cell, dm1, hermi=hermi, kpts=kpts) + v1 += -.5 * hyb * _get_k(mf, cell, dm1, hermi, kpts, kshift) return v1 else: # HF if (singlet is None or singlet) and hermi != 2: - def vind(dm1): - vj, vk = mf.get_jk(cell, dm1, hermi=hermi, kpts=kpts) + def vind(dm1, kshift=0): + vj, vk = _get_jk(mf, cell, dm1, hermi, kpts, kshift) return vj - .5 * vk else: - def vind(dm1): - return -.5 * mf.get_k(cell, dm1, hermi=hermi, kpts=kpts) + def vind(dm1, kshift=0): + return -.5 * _get_k(mf, cell, dm1, hermi, kpts, kshift) return vind @@ -152,7 +152,7 @@ def _gen_uhf_response(mf, mo_coeff=None, mo_occ=None, mem_now = lib.current_memory()[0] max_memory = max(2000, mf.max_memory*.8-mem_now) - def vind(dm1): + def vind(dm1, kshift=0): if hermi == 2: v1 = numpy.zeros_like(dm1) else: @@ -160,25 +160,25 @@ def vind(dm1): rho0, vxc, fxc, kpts, max_memory=max_memory) if not hybrid: if with_j: - vj = mf.get_j(cell, dm1, hermi=hermi, kpts=kpts) + vj = _get_j(mf, cell, dm1, hermi, kpts, kshift) v1 += vj[0] + vj[1] else: if with_j: - vj, vk = mf.get_jk(cell, dm1, hermi=hermi, kpts=kpts) + vj, vk = _get_jk(mf, cell, dm1, hermi, kpts, kshift) v1 += vj[0] + vj[1] - vk * hyb else: - v1 -= hyb * mf.get_k(cell, dm1, hermi=hermi, kpts=kpts) + v1 -= hyb * _get_k(mf, cell, dm1, hermi, kpts, kshift) return v1 elif with_j: - def vind(dm1): - vj, vk = mf.get_jk(cell, dm1, hermi=hermi, kpts=kpts) + def vind(dm1, kshift=0): + vj, vk = _get_jk(mf, cell, dm1, hermi, kpts, kshift) v1 = vj[0] + vj[1] - vk return v1 else: - def vind(dm1): - return -mf.get_k(cell, dm1, hermi=hermi, kpts=kpts) + def vind(dm1, kshift=0): + return -_get_k(mf, cell, dm1, hermi, kpts, kshift) return vind @@ -189,6 +189,29 @@ def _gen_ghf_response(mf, mo_coeff=None, mo_occ=None, ''' raise NotImplementedError +def _get_jk_kshift(mf, dm_kpts, hermi, kpts, kshift, with_j=True, with_k=True): + from pyscf.pbc.df.df_jk import get_j_kpts_kshift, get_k_kpts_kshift + vj = vk = None + if with_j: + vj = get_j_kpts_kshift(mf.with_df, dm_kpts, kshift, hermi=hermi, kpts=kpts) + if with_k: + vk = get_k_kpts_kshift(mf.with_df, dm_kpts, kshift, hermi=hermi, kpts=kpts, + exxdiv=mf.exxdiv) + return vj, vk +def _get_jk(mf, cell, dm1, hermi, kpts, kshift, with_j=True, with_k=True): + from pyscf.pbc import df + if kshift == 0: + return mf.get_jk(cell, dm1, hermi=hermi, kpts=kpts, with_j=with_j, with_k=with_k) + elif mf.rsjk is not None or not isinstance(mf.with_df, df.df.DF): + lib.logger.error(mf, 'Non-zero kshift is only supported by GDF/RSDF.') + raise NotImplementedError + else: + return _get_jk_kshift(mf, dm1, hermi, kpts, kshift, with_j=with_j, with_k=with_k) +def _get_j(mf, cell, dm1, hermi, kpts, kshift): + return _get_jk(mf, cell, dm1, hermi, kpts, kshift, True, False)[0] +def _get_k(mf, cell, dm1, hermi, kpts, kshift): + return _get_jk(mf, cell, dm1, hermi, kpts, kshift, False, True)[1] + khf.KRHF.gen_response = _gen_rhf_response kuhf.KUHF.gen_response = _gen_uhf_response @@ -372,4 +395,3 @@ def _gen_ghf_response_gam(mf, mo_coeff=None, mo_occ=None, uhf.UHF.gen_response = _gen_uhf_response_gam ghf.GHF.gen_response = _gen_ghf_response_gam rohf.ROHF.gen_response = _gen_uhf_response_gam - diff --git a/pyscf/pbc/scf/addons.py b/pyscf/pbc/scf/addons.py index d486a66fde..28ab6d72e5 100644 --- a/pyscf/pbc/scf/addons.py +++ b/pyscf/pbc/scf/addons.py @@ -579,6 +579,42 @@ def convert_to_khf(mf, out=None): out.__dict__.update(mf.__dict__) return out +def mo_energy_with_exxdiv_none(mf, mo_coeff=None): + ''' compute mo energy from the diagonal of fock matrix with exxdiv=None + ''' + from pyscf.pbc import scf, dft + from pyscf.lib import logger + log = logger.new_logger(mf) + + if mo_coeff is None: mo_coeff = mf.mo_coeff + + if mf.exxdiv is None and mf.mo_coeff is mo_coeff: + return mf.mo_energy + + with lib.temporary_env(mf, exxdiv=None): + dm = mf.make_rdm1(mo_coeff) + vhf = mf.get_veff(mf.mol, dm) + fockao = mf.get_fock(vhf=vhf, dm=dm) + + def _get_moe1(C, fao): + return lib.einsum('pi,pq,qi->i', C.conj(), fao, C) + def _get_moek(kC, kfao): + return [_get_moe1(C, fao) for C,fao in zip(kC, kfao)] + + # avoid using isinstance as some are other's derived class + if mf.__class__ in [scf.rhf.RHF, scf.ghf.GHF, dft.rks.RKS, dft.gks.GKS]: + return _get_moe1(mo_coeff, fockao) + elif mf.__class__ in [scf.uhf.UHF, dft.uks.UKS]: + return _get_moek(mo_coeff, fockao) + elif mf.__class__ in [scf.krhf.KRHF, scf.kghf.KGHF, dft.krks.KRKS, dft.kgks.KGKS]: + return _get_moek(mo_coeff, fockao) + elif mf.__class__ in [scf.kuhf.KUHF, dft.kuks.KUKS]: + return [_get_moek(kC, kfao) for kC,kfao in zip(mo_coeff,fockao)] + else: + log.error(f'Unknown SCF type {mf.__class__.__name__}') + raise NotImplementedError + + del (SMEARING_METHOD) diff --git a/pyscf/pbc/tdscf/__init__.py b/pyscf/pbc/tdscf/__init__.py index cd0e354d67..a06fe8834b 100644 --- a/pyscf/pbc/tdscf/__init__.py +++ b/pyscf/pbc/tdscf/__init__.py @@ -34,8 +34,6 @@ def TDHF(mf): import numpy if isinstance(mf, scf.khf.KSCF): return KTDHF(mf) - if numpy.abs(getattr(mf, 'kpt', 0)).max() > 1e-9: - raise NotImplementedError('Only supports gamma-point TDHF') if isinstance(mf, scf.hf.KohnShamDFT): raise RuntimeError('TDHF does not support DFT object %s' % mf) #TODO: mf = mf.remove_soscf() @@ -48,8 +46,6 @@ def TDA(mf): import numpy if isinstance(mf, scf.khf.KSCF): return KTDA(mf) - if numpy.abs(getattr(mf, 'kpt', 0)).max() > 1e-9: - raise NotImplementedError('Only supports gamma-point TDA') #TODO: mf = mf.remove_soscf() if isinstance(mf, scf.rohf.ROHF): if isinstance(mf, scf.hf.KohnShamDFT): @@ -62,8 +58,6 @@ def TDDFT(mf): import numpy if isinstance(mf, scf.khf.KSCF): return KTDDFT(mf) - if numpy.abs(getattr(mf, 'kpt', 0)).max() > 1e-9: - raise NotImplementedError('Only supports gamma-point TDDFT') if isinstance(mf, scf.hf.KohnShamDFT): #TODO: mf = mf.remove_soscf() if isinstance(mf, scf.rohf.ROHF): @@ -99,4 +93,3 @@ def KTDDFT(mf): return KTDHF(mf) KTD = KTDDFT - diff --git a/pyscf/pbc/tdscf/krhf.py b/pyscf/pbc/tdscf/krhf.py index e7740a6652..7ff6f791bc 100644 --- a/pyscf/pbc/tdscf/krhf.py +++ b/pyscf/pbc/tdscf/krhf.py @@ -29,39 +29,93 @@ from pyscf.pbc import scf from pyscf.pbc.tdscf.rhf import TDMixin from pyscf.pbc.scf import _response_functions # noqa -from pyscf.pbc.lib.kpts_helper import gamma_point +from pyscf.pbc.lib.kpts_helper import gamma_point, get_kconserv_ria from pyscf.pbc.df.df_ao2mo import warn_pbc2d_eri +from pyscf.pbc import df as pbcdf +from pyscf.data import nist from pyscf import __config__ REAL_EIG_THRESHOLD = getattr(__config__, 'pbc_tdscf_rhf_TDDFT_pick_eig_threshold', 1e-3) class KTDMixin(TDMixin): - def __init__(self, mf): + def __init__(self, mf, kshift_lst=None): assert isinstance(mf, scf.khf.KSCF) TDMixin.__init__(self, mf) 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 + self._keys.update(['kconserv', 'kshift_lst']) + + def dump_flags(self, verbose=None): + log = logger.new_logger(self, verbose) + log.info('\n') + log.info('******** %s for %s ********', + self.__class__, self._scf.__class__) + if self.singlet is None: + log.info('nstates = %d', self.nstates) + elif self.singlet: + log.info('nstates = %d singlet', self.nstates) + else: + log.info('nstates = %d triplet', self.nstates) + log.info('deg_eia_thresh = %.3e', self.deg_eia_thresh) + log.info('kshift_lst = %s', self.kshift_lst) + log.info('wfnsym = %s', self.wfnsym) + 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)', + self.max_memory, lib.current_memory()[0]) + if not self._scf.converged: + log.warn('Ground state SCF is not converged') + log.info('\n') + + def check_sanity(self): + TDMixin.check_sanity(self) + mf = self._scf + if any([k != 0 for k in self.kshift_lst]): + if mf.rsjk is not None or not isinstance(mf.with_df, pbcdf.df.DF): + logger.error(self, 'Solutions with non-zero kshift for %s are ' + 'only supported by GDF/RSDF') + raise NotImplementedError + + def _finalize(self): + '''Hook for dumping results and clearing up the object.''' + for k,kshift in enumerate(self.kshift_lst): + if not all(self.converged[k]): + logger.note(self, 'kshift = %d TD-SCF states %s not converged.', + kshift, [i for i, x in enumerate(self.converged[k]) if not x]) + logger.note(self, 'kshift = %d Excited State energies (eV)\n%s', + kshift, self.e[k] * nist.HARTREE2EV) + return self + get_nto = lib.invalid_method('get_nto') + class TDA(KTDMixin): conv_tol = getattr(__config__, 'pbc_tdscf_rhf_TDA_conv_tol', 1e-6) - def gen_vind(self, mf): + 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 singlet = self.singlet + kconserv = self.kconserv[kshift] mo_coeff = mf.mo_coeff - mo_energy = mf.mo_energy mo_occ = mf.mo_occ nkpts = len(mo_occ) nao, nmo = mo_coeff[0].shape 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[k][:,viridx[k]] for k in range(nkpts)] - e_ia = _get_e_ia(mo_energy, mo_occ) + orbv = [mo_coeff[kconserv[k]][:,viridx[kconserv[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]) mem_now = lib.current_memory()[0] @@ -70,7 +124,7 @@ def gen_vind(self, mf): def vind(zs): nz = len(zs) - z1s = [_unpack(z, mo_occ) for z in 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): @@ -79,30 +133,36 @@ def vind(zs): dmov[i,k] = reduce(numpy.dot, (orbo[k], dm1, orbv[k].conj().T)) with lib.temporary_env(mf, exxdiv=None): - v1ao = vresp(dmov) + v1ao = vresp(dmov, kshift) v1s = [] 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])) v1vo += e_ia[k] * dm1[k] - v1s.append(v1vo.ravel()) + v1.append(v1vo.ravel()) + v1s.append( numpy.concatenate(v1) ) return lib.asarray(v1s).reshape(nz,-1) return vind, hdiag - def init_guess(self, mf, nstates=None): + def init_guess(self, mf, kshift, nstates=None): if nstates is None: nstates = self.nstates mo_energy = mf.mo_energy mo_occ = mf.mo_occ - e_ia = numpy.hstack([x.ravel() for x in _get_e_ia(mo_energy, mo_occ)]) + kconserv = self.kconserv[kshift] + e_ia = numpy.concatenate( [x.reshape(-1) for x in + _get_e_ia(mo_energy, mo_occ, kconserv)] ) + e_ia = numpy.ceil(e_ia / self.deg_eia_thresh) * self.deg_eia_thresh + e_ia_uniq = numpy.unique(e_ia) e_ia_max = e_ia.max() nov = e_ia.size nstates = min(nstates, nov) - e_threshold = min(e_ia_max, e_ia[numpy.argsort(e_ia)[nstates-1]]) - # Handle degeneracy, include all degenerated states in initial guess - e_threshold += 1e-6 + nstates_thresh = min(nstates, e_ia_uniq.size) + e_threshold = min(e_ia_max, e_ia_uniq[numpy.argsort(e_ia_uniq)[nstates_thresh-1]]) + e_threshold += self.deg_eia_thresh idx = numpy.where(e_ia <= e_threshold)[0] x0 = numpy.zeros((idx.size, nov)) @@ -112,12 +172,20 @@ def init_guess(self, mf, nstates=None): def kernel(self, x0=None): '''TDA diagonalization solver + + Args: + x0: list of init guess arrays for each k-shift specified in :attr:`self.kshift_lst` + [x0_1, x0_2, ..., x0_nshift] + x0_i ~ (nstates, nkpts*nocc*nvir) ''' + cpu0 = (logger.process_clock(), logger.perf_counter()) self.check_sanity() self.dump_flags() - vind, hdiag = self.gen_vind(self._scf) - precond = self.get_precond(hdiag) + log = logger.new_logger(self) + + mf = self._scf + mo_occ = mf.mo_occ def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] @@ -126,41 +194,58 @@ def pickeig(w, v, nroots, envs): log = logger.Logger(self.stdout, self.verbose) precision = self.cell.precision * 1e-2 - 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=self.nstates, lindep=self.lindep, - max_space=self.max_space, pick=pickeig, - fill_heff=purify_krlyov_heff(precision, 0, log), - verbose=self.verbose) - - mo_occ = self._scf.mo_occ -# 1/sqrt(2) because self.x is for alpha excitation amplitude and 2(X^+*X) = 1 - self.xy = [(_unpack(xi*numpy.sqrt(.5), mo_occ), 0) for xi in x1] + self.converged = [] + self.e = [] + self.xy = [] + for i,kshift in enumerate(self.kshift_lst): + kconserv = self.kconserv[kshift] + + vind, hdiag = self.gen_vind(self._scf, kshift) + precond = self.get_precond(hdiag) + + if x0 is None: + x0k = self.init_guess(self._scf, kshift, self.nstates) + x0k = self.trunc_workspace(vind, x0k, nstates=self.nstates, pick=pickeig)[1] + 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) + self.converged.append( converged ) + self.e.append( e ) + # 1/sqrt(2) because self.x is for alpha excitation amplitude and 2(X^+*X) = 1 + self.xy.append( [(_unpack(xi*numpy.sqrt(.5), mo_occ, kconserv), 0) for xi in x1] ) + + log.timer(self.__class__.__name__, *cpu0) + self._finalize() return self.e, self.xy CIS = KTDA = TDA class TDHF(TDA): - def gen_vind(self, mf): + def gen_vind(self, mf, kshift): ''' [ A B ][X] [-B* -A*][Y] ''' singlet = self.singlet + kconserv = self.kconserv[kshift] mo_coeff = mf.mo_coeff - mo_energy = mf.mo_energy mo_occ = mf.mo_occ nkpts = len(mo_occ) nao, nmo = mo_coeff[0].shape 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[k][:,viridx[k]] for k in range(nkpts)] - e_ia = _get_e_ia(mo_energy, mo_occ) + orbv = [mo_coeff[kconserv[k]][:,viridx[kconserv[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]) tot_x = hdiag.size hdiag = numpy.hstack((hdiag, -hdiag)) @@ -171,8 +256,8 @@ def gen_vind(self, mf): def vind(xys): nz = len(xys) - z1xs = [_unpack(xy[:tot_x], mo_occ) for xy in xys] - z1ys = [_unpack(xy[tot_x:], mo_occ) for xy in 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) for i in range(nz): for k in range(nkpts): @@ -183,7 +268,7 @@ def vind(xys): dmov[i,k]+= reduce(numpy.dot, (orbv[k], dmy.T, orbo[k].T.conj())) with lib.temporary_env(mf, exxdiv=None): - v1ao = vresp(dmov) + v1ao = vresp(dmov, kshift) v1s = [] for i in range(nz): dmx = z1xs[i] @@ -197,29 +282,27 @@ def vind(xys): v1y+= e_ia[k] * dmy[k] v1xs.append(v1x.ravel()) v1ys.append(-v1y.ravel()) - v1s += v1xs + v1ys + # v1s += v1xs + v1ys + v1s.append( numpy.concatenate(v1xs + v1ys) ) return lib.asarray(v1s).reshape(nz,-1) return vind, hdiag - def init_guess(self, mf, nstates=None): - x0 = TDA.init_guess(self, mf, nstates) + 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()]])) def kernel(self, x0=None): '''TDHF diagonalization with non-Hermitian eigenvalue solver ''' - logger.warn(self, 'PBC-TDDFT is an experimental feature. ' - 'It is numerically sensitive to the accuracy of integrals ' - '(relating to cell.precision).') - + cpu0 = (logger.process_clock(), logger.perf_counter()) self.check_sanity() self.dump_flags() - vind, hdiag = self.gen_vind(self._scf) - precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) + log = logger.new_logger(self) + + mf = self._scf + mo_occ = mf.mo_occ real_system = (gamma_point(self._scf.kpts) and self._scf.mo_coeff[0].dtype == numpy.double) @@ -234,42 +317,65 @@ def pickeig(w, v, nroots, envs): precision = self.cell.precision * 1e-2 hermi = 0 - self.converged, w, x1 = \ - lib.davidson_nosym1(vind, x0, precond, - tol=self.conv_tol, - nroots=self.nstates, lindep=self.lindep, - max_space=self.max_space, pick=pickeig, - fill_heff=purify_krlyov_heff(precision, hermi, log), - verbose=self.verbose) - mo_occ = self._scf.mo_occ - self.e = w - def norm_xy(z): + def norm_xy(z, kconserv): x, y = z.reshape(2,-1) norm = 2*(lib.norm(x)**2 - lib.norm(y)**2) norm = 1/numpy.sqrt(norm) x *= norm y *= norm - return _unpack(x, mo_occ), _unpack(y, mo_occ) - self.xy = [norm_xy(z) for z in x1] - + return _unpack(x, mo_occ, kconserv), _unpack(y, mo_occ, kconserv) + + self.converged = [] + self.e = [] + self.xy = [] + for i,kshift in enumerate(self.kshift_lst): + kconserv = self.kconserv[kshift] + + vind, hdiag = self.gen_vind(self._scf, kshift) + precond = self.get_precond(hdiag) + + if x0 is None: + x0k = self.init_guess(self._scf, kshift, self.nstates) + x0k = self.trunc_workspace(vind, x0k, nstates=self.nstates*2, pick=pickeig)[1] + 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) + self.converged.append( converged ) + self.e.append( e ) + self.xy.append( [norm_xy(z, kconserv) for z in x1] ) + + log.timer(self.__class__.__name__, *cpu0) + self._finalize() return self.e, self.xy RPA = KTDHF = TDHF - -def _get_e_ia(mo_energy, mo_occ): +def _get_e_ia(mo_energy, mo_occ, kconserv=None): e_ia = [] - for k, occ in enumerate(mo_occ): - occidx = occ > 0 - viridx = occ == 0 - e_ia.append(mo_energy[k][viridx] - mo_energy[k][occidx,None]) + nkpts = len(mo_occ) + if kconserv is None: kconserv = numpy.arange(nkpts) + for k in range(nkpts): + kp = kconserv[k] + moeocc = mo_energy[k][mo_occ[k] > 1e-6] + moevir = mo_energy[kp][mo_occ[kp] < 1e-6] + e_ia.append( -moeocc[:,None] + moevir ) return e_ia -def _unpack(vo, mo_occ): +def _unpack(vo, mo_occ, kconserv): z = [] p1 = 0 for k, occ in enumerate(mo_occ): no = numpy.count_nonzero(occ > 0) - nv = occ.size - no + no1 = numpy.count_nonzero(mo_occ[kconserv[k]] > 0) + nv = occ.size - no1 p0, p1 = p1, p1 + no * nv z.append(vo[p0:p1].reshape(no,nv)) return z @@ -330,7 +436,7 @@ def fill_heff(heff, xs, ax, xt, axt, dot): td = TDA(mf) td.verbose = 5 - print(td.kernel()[0] * 27.2114) + 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] @@ -345,9 +451,8 @@ def fill_heff(heff, xs, ax, xt, axt, dot): td = TDHF(mf) td.verbose = 5 - print(td.kernel()[0] * 27.2114) + 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 33ed716df8..b38c8d309a 100644 --- a/pyscf/pbc/tdscf/krks.py +++ b/pyscf/pbc/tdscf/krks.py @@ -22,11 +22,34 @@ from pyscf import lib from pyscf.pbc import dft +from pyscf.pbc import df from pyscf.pbc.tdscf import krhf -KTDA = TDA = krhf.TDA -RPA = KTDDFT = TDDFT = krhf.TDHF +class TDA(krhf.TDA): + def kernel(self, x0=None): + _rebuild_df(self) + return krhf.TDA.kernel(self, x0=x0) + +KTDA = TDA + +class TDDFT(krhf.TDHF): + def kernel(self, x0=None): + _rebuild_df(self) + return krhf.TDHF.kernel(self, x0=x0) + +RPA = KTDDFT = TDDFT + +def _rebuild_df(td): + log = lib.logger.new_logger(td) + mf = td._scf + if any([k != 0 for k in td.kshift_lst]): + if isinstance(mf.with_df, df.df.DF): + if mf.with_df._j_only: + log.warn(f'Non-zero kshift is requested for {td.__class__.__name__}, ' + f'recomputing DF integrals with _j_only = False') + mf.with_df._j_only = False + mf.with_df.build() dft.krks.KRKS.TDA = lib.class_as_method(KTDA) dft.krks.KRKS.TDHF = None @@ -38,7 +61,7 @@ if __name__ == '__main__': from pyscf.pbc import gto - from pyscf.pbc import dft, df + from pyscf.pbc import dft cell = gto.Cell() cell.unit = 'B' cell.atom = ''' @@ -68,13 +91,13 @@ td = TDDFT(mf) td.nstates = 5 td.verbose = 5 - print(td.kernel()[0] * 27.2114) + 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] * 27.2114) + 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 e6f88ba03e..5452722788 100644 --- a/pyscf/pbc/tdscf/kuhf.py +++ b/pyscf/pbc/tdscf/kuhf.py @@ -32,10 +32,10 @@ class TDA(KTDMixin): conv_tol = getattr(__config__, 'pbc_tdscf_rhf_TDA_conv_tol', 1e-6) - def gen_vind(self, mf): + def gen_vind(self, mf, kshift): '''Compute Ax''' + kconserv = self.kconserv[kshift] mo_coeff = mf.mo_coeff - mo_energy = mf.mo_energy mo_occ = mf.mo_occ nkpts = len(mo_occ[0]) nao, nmo = mo_coeff[0][0].shape @@ -45,11 +45,12 @@ def gen_vind(self, mf): 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][k][:,viridxa[k]] for k in range(nkpts)] - orbvb = [mo_coeff[1][k][:,viridxb[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)] - e_ia_a = _get_e_ia(mo_energy[0], mo_occ[0]) - e_ia_b = _get_e_ia(mo_energy[1], mo_occ[1]) + 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) hdiag = numpy.hstack([x.ravel() for x in (e_ia_a + e_ia_b)]) mem_now = lib.current_memory()[0] @@ -58,7 +59,7 @@ def gen_vind(self, mf): def vind(zs): nz = len(zs) - zs = [_unpack(z, mo_occ) for z in zs] + zs = [_unpack(z, mo_occ, kconserv) for z in zs] dmov = numpy.empty((2,nz,nkpts,nao,nao), dtype=numpy.complex128) for i in range(nz): dm1a, dm1b = zs[i] @@ -67,8 +68,8 @@ def vind(zs): dmov[1,i,k] = reduce(numpy.dot, (orbob[k], dm1b[k], orbvb[k].conj().T)) with lib.temporary_env(mf, exxdiv=None): - dmov = dmov.reshape(2*nz,nkpts,nao,nao) - v1ao = vresp(dmov) + dmov = dmov.reshape(2,nz,nkpts,nao,nao) + v1ao = vresp(dmov, kshift) v1ao = v1ao.reshape(2,nz,nkpts,nao,nao) v1s = [] @@ -83,26 +84,30 @@ def vind(zs): v1b += e_ia_b[k] * dm1b[k] v1as.append(v1a.ravel()) v1bs.append(v1b.ravel()) - v1s += v1as + v1bs - return numpy.hstack(v1s).reshape(nz,-1) + v1s.append( numpy.concatenate(v1as + v1bs) ) + return lib.asarray(v1s).reshape(nz,-1) return vind, hdiag - def init_guess(self, mf, nstates=None): + def init_guess(self, mf, kshift, nstates=None): if nstates is None: nstates = self.nstates mo_energy = mf.mo_energy mo_occ = mf.mo_occ - e_ia_a = _get_e_ia(mo_energy[0], mo_occ[0]) - e_ia_b = _get_e_ia(mo_energy[1], mo_occ[1]) + kconserv = self.kconserv[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)]) + e_ia = numpy.ceil(e_ia / self.deg_eia_thresh) * self.deg_eia_thresh + e_ia_uniq = numpy.unique(e_ia) e_ia_max = e_ia.max() nov = e_ia.size nstates = min(nstates, nov) - e_threshold = min(e_ia_max, e_ia[numpy.argsort(e_ia)[nstates-1]]) - # Handle degeneracy - e_threshold += 1e-6 + nstates_thresh = min(nstates, e_ia_uniq.size) + e_threshold = min(e_ia_max, e_ia_uniq[numpy.argsort(e_ia_uniq)[nstates_thresh-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): @@ -112,13 +117,14 @@ def init_guess(self, mf, nstates=None): def kernel(self, x0=None): '''TDA diagonalization solver ''' + cpu0 = (logger.process_clock(), logger.perf_counter()) self.check_sanity() self.dump_flags() - vind, hdiag = self.gen_vind(self._scf) - precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) + log = logger.new_logger(self) + + mf = self._scf + mo_occ = mf.mo_occ def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] @@ -128,27 +134,46 @@ def pickeig(w, v, nroots, envs): precision = self.cell.precision * 1e-2 hermi = 1 - self.converged, self.e, x1 = \ - lib.davidson1(vind, x0, precond, - tol=self.conv_tol, - nroots=self.nstates, lindep=self.lindep, - max_space=self.max_space, pick=pickeig, - fill_heff=purify_krlyov_heff(precision, hermi, log), - verbose=self.verbose) - - mo_occ = self._scf.mo_occ - self.xy = [(_unpack(xi, mo_occ), # (X_alpha, X_beta) - (0, 0)) # (Y_alpha, Y_beta) - for xi in x1] + self.converged = [] + self.e = [] + self.xy = [] + for i,kshift in enumerate(self.kshift_lst): + kconserv = self.kconserv[kshift] + + vind, hdiag = self.gen_vind(self._scf, kshift) + precond = self.get_precond(hdiag) + + if x0 is None: + x0k = self.init_guess(self._scf, kshift, self.nstates) + x0k = self.trunc_workspace(vind, x0k, nstates=self.nstates, pick=pickeig)[1] + 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) + self.converged.append( converged ) + self.e.append( e ) + self.xy.append( [(_unpack(xi, mo_occ, kconserv), # (X_alpha, X_beta) + (0, 0)) # (Y_alpha, Y_beta) + for xi in x1] ) #TODO: analyze CIS wfn point group symmetry + log.timer(self.__class__.__name__, *cpu0) + self._finalize() return self.e, self.xy CIS = KTDA = TDA class TDHF(TDA): - def gen_vind(self, mf): + def gen_vind(self, mf, kshift): + kconserv = self.kconserv[kshift] mo_coeff = mf.mo_coeff - mo_energy = mf.mo_energy mo_occ = mf.mo_occ nkpts = len(mo_occ[0]) nao, nmo = mo_coeff[0][0].shape @@ -158,11 +183,12 @@ def gen_vind(self, mf): 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][k][:,viridxa[k]] for k in range(nkpts)] - orbvb = [mo_coeff[1][k][:,viridxb[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)] - e_ia_a = _get_e_ia(mo_energy[0], mo_occ[0]) - e_ia_b = _get_e_ia(mo_energy[1], mo_occ[1]) + 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) hdiag = numpy.hstack([x.ravel() for x in (e_ia_a + e_ia_b)]) hdiag = numpy.hstack((hdiag, -hdiag)) tot_x_a = sum(x.size for x in e_ia_a) @@ -175,8 +201,8 @@ def gen_vind(self, mf): def vind(xys): nz = len(xys) - x1s = [_unpack(x[:tot_x], mo_occ) for x in xys] - y1s = [_unpack(x[tot_x:], mo_occ) for x in xys] + x1s = [_unpack(x[:tot_x], mo_occ, kconserv) for x in xys] + y1s = [_unpack(x[tot_x:], mo_occ, kconserv) for x in xys] dmov = numpy.empty((2,nz,nkpts,nao,nao), dtype=numpy.complex128) for i in range(nz): xa, xb = x1s[i] @@ -190,8 +216,8 @@ def vind(xys): dmov[1,i,k] = dmx + dmy # AX + BY with lib.temporary_env(mf, exxdiv=None): - dmov = dmov.reshape(2*nz,nkpts,nao,nao) - v1ao = vresp(dmov) + dmov = dmov.reshape(2,nz,nkpts,nao,nao) + v1ao = vresp(dmov, kshift) v1ao = v1ao.reshape(2,nz,nkpts,nao,nao) v1s = [] @@ -215,26 +241,27 @@ def vind(xys): v1xsb.append(v1xb.ravel()) v1ysa.append(-v1ya.ravel()) v1ysb.append(-v1yb.ravel()) - v1s += v1xsa + v1xsb + v1ysa + v1ysb + v1s.append( numpy.concatenate(v1xsa + v1xsb + v1ysa + v1ysb) ) return numpy.hstack(v1s).reshape(nz,-1) return vind, hdiag - def init_guess(self, mf, nstates=None, wfnsym=None): - x0 = TDA.init_guess(self, mf, nstates) + 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()]])) def kernel(self, x0=None): '''TDHF diagonalization with non-Hermitian eigenvalue solver ''' + cpu0 = (logger.process_clock(), logger.perf_counter()) self.check_sanity() self.dump_flags() - vind, hdiag = self.gen_vind(self._scf) - precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) + log = logger.new_logger(self) + + mf = self._scf + mo_occ = mf.mo_occ real_system = (gamma_point(self._scf.kpts) and self._scf.mo_coeff[0][0].dtype == numpy.double) @@ -248,44 +275,66 @@ def pickeig(w, v, nroots, envs): log = logger.Logger(self.stdout, self.verbose) precision = self.cell.precision * 1e-2 - self.converged, w, x1 = \ - lib.davidson_nosym1(vind, x0, precond, - tol=self.conv_tol, - nroots=self.nstates, lindep=self.lindep, - max_space=self.max_space, pick=pickeig, - fill_heff=purify_krlyov_heff(precision, 0, log), - verbose=self.verbose) - - mo_occ = self._scf.mo_occ - e = [] - xy = [] - for i, z in enumerate(x1): - xs, ys = z.reshape(2,-1) - norm = lib.norm(xs)**2 - lib.norm(ys)**2 - if norm > 0: - norm = 1/numpy.sqrt(norm) - xs *= norm - ys *= norm - e.append(w[i]) - xy.append((_unpack(xs, mo_occ), _unpack(ys, mo_occ))) - self.e = numpy.array(e) - self.xy = xy + self.converged = [] + self.e = [] + self.xy = [] + for i,kshift in enumerate(self.kshift_lst): + kconserv = self.kconserv[kshift] + + vind, hdiag = self.gen_vind(self._scf, kshift) + precond = self.get_precond(hdiag) + + if x0 is None: + x0k = self.init_guess(self._scf, kshift, self.nstates) + x0k = self.trunc_workspace(vind, x0k, nstates=self.nstates, pick=pickeig)[1] + 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) + self.converged.append( converged ) + + e = [] + xy = [] + for i, z in enumerate(x1): + xs, ys = z.reshape(2,-1) + norm = lib.norm(xs)**2 - lib.norm(ys)**2 + if norm > 0: + norm = 1/numpy.sqrt(norm) + xs *= norm + ys *= norm + e.append(w[i]) + xy.append((_unpack(xs, mo_occ, kconserv), _unpack(ys, mo_occ, kconserv))) + self.e.append( numpy.array(e) ) + self.xy.append( xy ) + + log.timer(self.__class__.__name__, *cpu0) + self._finalize() return self.e, self.xy RPA = KTDHF = TDHF -def _unpack(vo, mo_occ): +def _unpack(vo, mo_occ, kconserv): za = [] zb = [] p1 = 0 for k, occ in enumerate(mo_occ[0]): no = numpy.count_nonzero(occ > 0) - nv = occ.size - no + no1 = numpy.count_nonzero(mo_occ[0][kconserv[k]] > 0) + nv = occ.size - no1 p0, p1 = p1, p1 + no * nv za.append(vo[p0:p1].reshape(no,nv)) for k, occ in enumerate(mo_occ[1]): no = numpy.count_nonzero(occ > 0) - nv = occ.size - no + no1 = numpy.count_nonzero(mo_occ[1][kconserv[k]] > 0) + nv = occ.size - no1 p0, p1 = p1, p1 + no * nv zb.append(vo[p0:p1].reshape(no,nv)) return za, zb @@ -325,12 +374,12 @@ def _unpack(vo, mo_occ): td = TDA(mf) td.verbose = 5 td.nstates = 5 - print(td.kernel()[0] * 27.2114) + print(td.kernel()[0][0] * 27.2114) td = TDHF(mf) td.verbose = 5 td.nstates = 5 - print(td.kernel()[0] * 27.2114) + print(td.kernel()[0][0] * 27.2114) cell.spin = 2 mf = scf.KUHF(cell, cell.make_kpts([2,1,1])).set(exxdiv=None) @@ -339,9 +388,9 @@ def _unpack(vo, mo_occ): td = TDA(mf) td.verbose = 5 td.nstates = 5 - print(td.kernel()[0] * 27.2114) + print(td.kernel()[0][0] * 27.2114) td = TDHF(mf) td.verbose = 5 td.nstates = 5 - print(td.kernel()[0] * 27.2114) + print(td.kernel()[0][0] * 27.2114) diff --git a/pyscf/pbc/tdscf/kuks.py b/pyscf/pbc/tdscf/kuks.py index ebfd237d8e..1f76d12c8f 100644 --- a/pyscf/pbc/tdscf/kuks.py +++ b/pyscf/pbc/tdscf/kuks.py @@ -19,10 +19,22 @@ from pyscf import lib from pyscf.pbc import dft from pyscf.pbc.tdscf import kuhf +from pyscf.pbc.tdscf.krks import _rebuild_df -KTDA = TDA = kuhf.TDA -RPA = KTDDFT = TDDFT = kuhf.TDHF +class TDA(kuhf.TDA): + def kernel(self, x0=None): + _rebuild_df(self) + kuhf.TDA.kernel(self, x0=x0) + +KTDA = TDA + +class TDDFT(kuhf.TDHF): + def kernel(self, x0=None): + _rebuild_df(self) + kuhf.TDHF.kernel(self, x0=x0) + +RPA = KTDDFT = TDDFT dft.kuks.KUKS.TDA = lib.class_as_method(KTDA) dft.kuks.KUKS.TDHF = None @@ -59,14 +71,14 @@ td = TDDFT(mf) td.verbose = 5 td.nstates = 5 - print(td.kernel()[0] * 27.2114) + 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] * 27.2114) + 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') @@ -75,11 +87,11 @@ td = TDDFT(mf) td.verbose = 5 td.nstates = 5 - print(td.kernel()[0] * 27.2114) + 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] * 27.2114) + print(td.kernel()[0][0] * 27.2114) diff --git a/pyscf/pbc/tdscf/rhf.py b/pyscf/pbc/tdscf/rhf.py index a20846e2b0..09f132fb6b 100644 --- a/pyscf/pbc/tdscf/rhf.py +++ b/pyscf/pbc/tdscf/rhf.py @@ -22,6 +22,7 @@ from pyscf import lib from pyscf.tdscf import rhf +from pyscf.pbc import scf from pyscf import __config__ class TDMixin(rhf.TDMixin): @@ -56,27 +57,9 @@ class TDA(TDMixin): _gen_vind = rhf.TDA.gen_vind def gen_vind(self, mf): - # gen_vind calls get_jk functions to compute the contraction between - # two-electron integrals and X,Y amplitudes. There are two choices for - # the treatment of exxdiv. - # 1. Removing exxdiv corrections in both orbital energies (in hdiag) and - # get_jk functions (in vind function). This treatment can make TDA the - # same to CIS method in which exxdiv was completely excluded when - # constructing Hamiltonians. - # 2. Excluding exxdiv corrections from get_jk only. Keep its correction - # to orbital energy. This treatment can make the TDDFT excitation - # energies closed to the relevant DFT orbital energy gaps. - # DFT orbital energy gaps can be used as a good estimation for - # excitation energies. Current implementation takes the second choice so - # as to make the TDDFT excitation energies agree to DFT orbital energy gaps. - # - # There might be a third treatment: Taking treatment 1 first then adding - # certain types of corrections to the excitation energy at last. - # I'm not sure how to do this properly. - # - # See also issue https://github.com/pyscf/pyscf/issues/1187 - - vind, hdiag = self._gen_vind(mf) + moe = scf.addons.mo_energy_with_exxdiv_none(mf) + with lib.temporary_env(mf, mo_energy=moe): + vind, hdiag = self._gen_vind(mf) def vindp(x): with lib.temporary_env(mf, exxdiv=None): return vind(x) @@ -95,9 +78,7 @@ class TDHF(TDA): RPA = TDRHF = TDHF -from pyscf.pbc import scf scf.hf.RHF.TDA = lib.class_as_method(TDA) scf.hf.RHF.TDHF = lib.class_as_method(TDHF) scf.rohf.ROHF.TDA = None scf.rohf.ROHF.TDHF = None - diff --git a/pyscf/pbc/tdscf/rks.py b/pyscf/pbc/tdscf/rks.py index 9995b45bd1..55da0c0911 100644 --- a/pyscf/pbc/tdscf/rks.py +++ b/pyscf/pbc/tdscf/rks.py @@ -19,29 +19,30 @@ from pyscf import lib from pyscf.pbc import dft from pyscf.tdscf import rks -from pyscf.pbc.tdscf.rhf import TDA -from pyscf.pbc.tdscf.rhf import TDHF as TDDFT +from pyscf.pbc.tdscf import rhf +from pyscf.pbc.lib.kpts_helper import gamma_point -RPA = TDRKS = TDDFT +RPA = TDRKS = TDDFT = rhf.TDHF -class CasidaTDDFT(TDA): +class CasidaTDDFT(rhf.TDA): _gen_vind = rks.TDDFTNoHybrid.gen_vind - gen_vind = TDA.gen_vind + gen_vind = rhf.TDA.gen_vind + kernel = rks.TDDFTNoHybrid.kernel TDDFTNoHybrid = CasidaTDDFT def tddft(mf): '''Driver to create TDDFT or CasidaTDDFT object''' - if mf._numint.libxc.is_hybrid_xc(mf.xc): - return TDDFT(mf) - else: + kpt = getattr(mf, 'kpt', None) + if not mf._numint.libxc.is_hybrid_xc(mf.xc) and gamma_point(kpt): return CasidaTDDFT(mf) + else: + return TDDFT(mf) -dft.rks.RKS.TDA = lib.class_as_method(TDA) +dft.rks.RKS.TDA = lib.class_as_method(rhf.TDA) dft.rks.RKS.TDHF = None -dft.rks.RKS.TDDFTNoHybrid = lib.class_as_method(TDDFTNoHybrid) +dft.rks.RKS.TDDFTNoHybrid = tddft dft.rks.RKS.CasidaTDDFT = lib.class_as_method(CasidaTDDFT) dft.rks.RKS.TDDFT = tddft #dft.rks.RKS.dTDA = lib.class_as_method(dTDA) #dft.rks.RKS.dRPA = lib.class_as_method(dRPA) - diff --git a/pyscf/pbc/tdscf/test/test_kproxy_ks.py b/pyscf/pbc/tdscf/test/test_kproxy_ks.py index 8182451522..4679408ff2 100644 --- a/pyscf/pbc/tdscf/test/test_kproxy_ks.py +++ b/pyscf/pbc/tdscf/test/test_kproxy_ks.py @@ -98,7 +98,7 @@ def test_class(self): for k in range(self.k): # Prepare indexes r1, r2, c1, c2 = krhf_slow.get_block_k_ix(model.eri, k) - r = k2k(r1, r2) + # r = k2k(r1, r2) c = k2k(c1, c2) # Select roots diff --git a/pyscf/pbc/tdscf/test/test_krhf.py b/pyscf/pbc/tdscf/test/test_krhf.py new file mode 100644 index 0000000000..eef979ed31 --- /dev/null +++ b/pyscf/pbc/tdscf/test/test_krhf.py @@ -0,0 +1,164 @@ +#!/usr/bin/env python +# Copyright 2021 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 unittest +import numpy as np +from pyscf import __config__ +from pyscf.pbc import gto, scf, tdscf, cc +from pyscf import gto as molgto, scf as molscf, tdscf as moltdscf +from pyscf.pbc.cc.eom_kccsd_rhf import EOMEESinglet +from pyscf.data.nist import HARTREE2EV as unitev + + +class Diamond(unittest.TestCase): + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-hf-rev' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + kpts = cell.make_kpts((2,1,1)) + mf = scf.KRHF(cell, kpts=kpts).rs_density_fit(auxbasis='weigend').run() + cls.cell = cell + cls.mf = mf + + cls.nstates = 4 # make sure first `nstates_test` states are converged + cls.nstates_test = 1 + @classmethod + 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)), + 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) + + def test_tda_singlet_eomccs(self): + ''' Brute-force solution to the KTDA equation. Compared to the brute-force + implementation of KTDA from EOM-EE-CCSD. + ''' + td = tdscf.KTDA(self.mf).set(kshift_lst=np.arange(len(self.mf.kpts))) + ecis_k = [] + for kshift in td.kshift_lst: + vind, hdiag = td.gen_vind(td._scf, kshift) + heff = vind(np.eye(hdiag.size)) + ecis_k.append(np.linalg.eigh(heff)[0]) + + mcc = cc.KCCSD(self.mf) + mcc.kernel() + meom = EOMEESinglet(mcc) + ecc_k = meom.cis(nroots=ecis_k[0].size, kptlist=td. kshift_lst)[0] + + for e,ecc in zip(ecis_k, ecc_k): + self.assertAlmostEqual(abs(e * unitev - ecc.real * unitev).max(), 0, 3) + + def test_tda_singlet(self): + ref = [[10.9573977036], + [11.0418311085]] + self.kernel(tdscf.KTDA, ref) + + def test_tda_triplet(self): + ref = [[6.4440137833], + [7.4264899075]] + self.kernel(tdscf.KTDA, ref, singlet=False) + + def test_tdhf_singlet(self): + ref = [[10.7665673889], + [10.8485048947]] + self.kernel(tdscf.KTDHF, ref) + + def test_tdhf_triplet(self): + ref = [[5.9794378466], + [6.1703909932]] + self.kernel(tdscf.KTDHF, ref, singlet=False) + + +class WaterBigBox(unittest.TestCase): + ''' Match molecular CIS + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = ''' + O 0.00000 0.00000 0.11779 + H 0.00000 0.75545 -0.47116 + H 0.00000 -0.75545 -0.47116 + ''' + cell.a = np.eye(3) * 15 + cell.basis = 'sto-3g' + cell.build() + kpts = cell.make_kpts((2,1,1)) + mf = scf.KRHF(cell, kpts=kpts).rs_density_fit(auxbasis='weigend') + mf.with_df.omega = 0.1 + mf.kernel() + cls.cell = cell + cls.mf = mf + + mol = molgto.Mole() + for key in ['verbose','output','atom','basis']: + setattr(mol, key, getattr(cell, key)) + mol.build() + molmf = molscf.RHF(mol).density_fit(auxbasis=mf.with_df.auxbasis).run() + cls.mol = mol + cls.molmf = molmf + + cls.nstates = 5 # make sure first `nstates_test` states are converged + cls.nstates_test = 3 + + @classmethod + def tearDownClass(cls): + cls.cell.stdout.close() + cls.mol.stdout.close() + del cls.cell, cls.mf + del cls.mol, cls.molmf + + def kernel(self, TD, MOLTD, **kwargs): + td = TD(self.mf).set(nstates=self.nstates, **kwargs).run() + moltd = MOLTD(self.molmf).set(nstates=self.nstates, **kwargs).run() + ref = moltd.e[:self.nstates_test] + for kshift,e in enumerate(td.e): + self.assertAlmostEqual(abs(e[:self.nstates_test] * unitev - ref * unitev).max(), 0, 2) + + def test_tda_singlet(self): + self.kernel(tdscf.KTDA, moltdscf.TDA) + + def test_tda_triplet(self): + self.kernel(tdscf.KTDA, moltdscf.TDA, singlet=False) + + def test_tdhf_singlet(self): + self.kernel(tdscf.KTDHF, moltdscf.TDHF) + + def test_tdhf_triplet(self): + self.kernel(tdscf.KTDHF, moltdscf.TDHF, singlet=False) + + +if __name__ == "__main__": + print("Full Tests for krhf-TDA and krhf-TDHF") + unittest.main() diff --git a/pyscf/pbc/tdscf/test/test_krhf_slow.py b/pyscf/pbc/tdscf/test/test_krhf_slow.py index b6b026d3de..7040dfa204 100644 --- a/pyscf/pbc/tdscf/test/test_krhf_slow.py +++ b/pyscf/pbc/tdscf/test/test_krhf_slow.py @@ -157,7 +157,7 @@ def test_eig_kernel(self): for k in range(self.k): # Prepare indexes r1, r2, c1, c2 = ktd.get_block_k_ix(model.eri, k) - r = k2k(r1, r2) + # r = k2k(r1, r2) c = k2k(c1, c2) # Select roots diff --git a/pyscf/pbc/tdscf/test/test_krhf_slow_gamma.py b/pyscf/pbc/tdscf/test/test_krhf_slow_gamma.py index 1e7b6f2204..445bd38c02 100644 --- a/pyscf/pbc/tdscf/test/test_krhf_slow_gamma.py +++ b/pyscf/pbc/tdscf/test/test_krhf_slow_gamma.py @@ -136,9 +136,9 @@ def test_class(self): model.kernel() mask_o, mask_v = tdhf_frozen_mask(model.eri, kind="o,v") testing.assert_allclose(model.e, self.td_model_krhf.e, atol=1e-3) - assert_vectors_close(model.xy, numpy.array(self.td_model_krhf.xy)[..., mask_o, :][..., mask_v], atol=1e-2) + assert_vectors_close(model.xy, + numpy.array(self.td_model_krhf.xy)[..., mask_o, :][..., mask_v], atol=1e-2) except Exception: print("When testing class with frozen={} the following exception occurred:".format(repr(frozen))) raise - diff --git a/pyscf/pbc/tdscf/test/test_krks.py b/pyscf/pbc/tdscf/test/test_krks.py new file mode 100644 index 0000000000..c99315a043 --- /dev/null +++ b/pyscf/pbc/tdscf/test/test_krks.py @@ -0,0 +1,146 @@ +#!/usr/bin/env python +# Copyright 2021 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 unittest +import numpy as np +from pyscf import __config__ +from pyscf.pbc import gto, scf, tdscf, cc +from pyscf import gto as molgto, scf as molscf, tdscf as moltdscf +from pyscf.pbc.cc.eom_kccsd_rhf import EOMEESinglet +from pyscf.data.nist import HARTREE2EV as unitev + + +class DiamondPBE(unittest.TestCase): + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-hf-rev' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + kpts = cell.make_kpts((2,1,1)) + + xc = 'pbe' + mf = scf.KRKS(cell, kpts=kpts).set(xc=xc).rs_density_fit(auxbasis='weigend').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)), + 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) + + def test_tda_singlet(self): + ref = [[7.7172854747, 7.7173219160], + [8.3749594280, 8.3749980463]] + self.kernel('TDA', ref) + + def test_tda_triplet(self): + ref = [[5.7465112548, 5.7465526327], + [6.9888184993, 6.9888609925]] + self.kernel('TDA', ref, singlet=False) + + def test_tdhf_singlet(self): + ref = [[7.5824302393, 7.5824675688], + [8.3438648420, 8.3439036271]] + self.kernel('TDDFT', ref) + + def test_tdhf_triplet(self): + ref = [[5.5659966435, 5.5660393021], + [6.7992845776, 6.7993290046]] + self.kernel('TDDFT', ref, singlet=False) + + +class DiamondPBE0(unittest.TestCase): + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-hf-rev' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + kpts = cell.make_kpts((2,1,1)) + + xc = 'pbe0' + mf = scf.KRKS(cell, kpts=kpts).set(xc=xc).rs_density_fit(auxbasis='weigend').run() + + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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)), + 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) + + def test_tda_singlet(self): + ref = [[9.3936718451, 9.4874866060], + [10.0697605303, 10.0697862958]] + self.kernel('TDA', ref) + + def test_tda_triplet(self): + ref = [[6.6703797643, 6.6704110631], + [7.4081863259, 7.4082204017]] + self.kernel('TDA', ref, singlet=False) + + def test_tdhf_singlet(self): + ref = [[9.2519208050, 9.3208025447], + [9.9609751875, 9.9610015227]] + self.kernel('TDDFT', ref) + + def test_tdhf_triplet(self): + ref = [[6.3282716764, 6.3283051217], + [7.0656766298, 7.0657111705]] + self.kernel('TDDFT', ref, singlet=False) + + +if __name__ == "__main__": + print("Full Tests for krks-TDA and krks-TDDFT") + unittest.main() diff --git a/pyscf/pbc/tdscf/test/test_kuhf.py b/pyscf/pbc/tdscf/test/test_kuhf.py new file mode 100644 index 0000000000..903090c8d2 --- /dev/null +++ b/pyscf/pbc/tdscf/test/test_kuhf.py @@ -0,0 +1,136 @@ +#!/usr/bin/env python +# Copyright 2021 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 unittest +import numpy as np +from pyscf import __config__ +from pyscf.pbc import gto, scf, tdscf +from pyscf import gto as molgto, scf as molscf, tdscf as moltdscf +from pyscf.data.nist import HARTREE2EV as unitev + + +class Diamond(unittest.TestCase): + ''' Reproduce KRHF-TDSCF + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-hf-rev' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + kpts = cell.make_kpts((2,1,1)) + mf = scf.KUHF(cell, kpts=kpts).rs_density_fit(auxbasis='weigend').run() + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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 = TD(self.mf).set(kshift_lst=np.arange(len(self.mf.kpts)), + 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) + + 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) + + 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) + + +class WaterBigBox(unittest.TestCase): + ''' Match molecular CIS + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = ''' + O 0.00000 0.00000 0.11779 + H 0.00000 0.75545 -0.47116 + H 0.00000 -0.75545 -0.47116 + ''' + cell.spin = 4 # TOTAL spin in the corresponding supercell + cell.a = np.eye(3) * 15 + cell.basis = 'sto-3g' + cell.build() + kpts = cell.make_kpts((2,1,1)) + mf = scf.KUHF(cell, kpts=kpts).rs_density_fit(auxbasis='weigend') + mf.with_df.omega = 0.1 + mf.kernel() + cls.cell = cell + cls.mf = mf + + mol = molgto.Mole() + for key in ['verbose','output','atom','basis']: + setattr(mol, key, getattr(cell, key)) + mol.spin = cell.spin // len(kpts) + mol.build() + molmf = molscf.UHF(mol).density_fit(auxbasis=mf.with_df.auxbasis).run() + cls.mol = mol + cls.molmf = molmf + + cls.nstates = 5 # make sure first `nstates_test` states are converged + cls.nstates_test = 2 + + @classmethod + def tearDownClass(cls): + cls.cell.stdout.close() + cls.mol.stdout.close() + del cls.cell, cls.mf + 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() + moltd = MOLTD(self.molmf).set(nstates=self.nstates, **kwargs).run() + ref = moltd.e + for kshift,e in enumerate(td.e): + self.assertAlmostEqual(abs(e[:self.nstates_test] * unitev - + ref[:self.nstates_test] * unitev).max(), 0, 2) + + def test_tda(self): + self.kernel(tdscf.KTDA, moltdscf.TDA) + + def test_tdhf(self): + self.kernel(tdscf.KTDHF, moltdscf.TDHF) + + +if __name__ == "__main__": + print("Full Tests for kuhf-TDA and kuhf-TDHF") + unittest.main() diff --git a/pyscf/pbc/tdscf/test/test_kuks.py b/pyscf/pbc/tdscf/test/test_kuks.py new file mode 100644 index 0000000000..aa05bab624 --- /dev/null +++ b/pyscf/pbc/tdscf/test/test_kuks.py @@ -0,0 +1,134 @@ +#!/usr/bin/env python +# Copyright 2021 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 unittest +import numpy as np +from pyscf import __config__ +from pyscf.pbc import gto, scf, tdscf, cc +from pyscf import gto as molgto, scf as molscf, tdscf as moltdscf +from pyscf.pbc.cc.eom_kccsd_rhf import EOMEESinglet +from pyscf.data.nist import HARTREE2EV as unitev + + +class DiamondPBE(unittest.TestCase): + ''' Reproduce KRKS-TDSCF + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-hf-rev' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + kpts = cell.make_kpts((2,1,1)) + + xc = 'pbe' + mf = scf.KUKS(cell, kpts=kpts).set(xc=xc).rs_density_fit(auxbasis='weigend').run() + + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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)), + 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) + + def test_tda(self): + # same as lowest roots in DiamondPBE->test_tda_singlet/triplet in test_krks.py + ref = [[5.7465112548, 5.7465526327], + [6.9888184993, 6.9888609925]] + self.kernel('TDA', ref, singlet=False) + + 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]] + self.kernel('TDDFT', ref, singlet=False) + + +class DiamondPBE0(unittest.TestCase): + ''' Reproduce KRKS-TDSCF + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-hf-rev' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + kpts = cell.make_kpts((2,1,1)) + + xc = 'pbe0' + mf = scf.KUKS(cell, kpts=kpts).set(xc=xc).rs_density_fit(auxbasis='weigend').run() + + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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)), + 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) + + def test_tda(self): + # same as lowest roots in DiamondPBE0->test_tda_singlet/triplet in test_krks.py + ref = [[6.6703797643, 6.6704110631], + [7.4081863259, 7.4082204017]] + self.kernel('TDA', ref, singlet=False) + + 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]] + self.kernel('TDDFT', ref, singlet=False) + + +if __name__ == "__main__": + print("Full Tests for kuks-TDA and kuks-TDDFT") + unittest.main() diff --git a/pyscf/pbc/tdscf/test/test_rhf.py b/pyscf/pbc/tdscf/test/test_rhf.py index 957a534630..566e2687fb 100644 --- a/pyscf/pbc/tdscf/test/test_rhf.py +++ b/pyscf/pbc/tdscf/test/test_rhf.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -# Copyright 2014-2018 The PySCF Developers. All Rights Reserved. +# Copyright 2021 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. @@ -13,65 +13,176 @@ # See the License for the specific language governing permissions and # limitations under the License. # -# Author: Qiming Sun -# import unittest -from pyscf import lib -from pyscf.pbc.gto import Cell -from pyscf.pbc.scf import RHF, KRHF -from pyscf.pbc import tdscf - -def setUpModule(): - global cell - cell = Cell() - cell.atom = ''' - C 0.000000000000 0.000000000000 0.000000000000 - C 1.685068664391 1.685068664391 1.685068664391 +import numpy as np +from pyscf import __config__ +from pyscf import gto as molgto, scf as molscf, tdscf as moltdscf +from pyscf.pbc import gto, scf, tdscf +from pyscf.data.nist import HARTREE2EV as unitev + + +class Diamond(unittest.TestCase): + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-hf-rev' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + + mf = scf.RHF(cell).rs_density_fit(auxbasis='weigend').run() + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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(nstates=self.nstates, **kwargs).run() + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + + def test_tda_singlet(self): + ref = [9.6425852906, 9.6425852906] + self.kernel('TDA', ref) + + def test_tda_triplet(self): + ref = [4.7209460258, 5.6500725912] + self.kernel('TDA', ref, singlet=False) + + def test_tdhf_singlet(self): + ref = [9.2573219105, 9.2573219106] + self.kernel('TDHF', ref) + + def test_tdhf_triplet(self): + ref = [3.0396052214, 3.0396052214] + self.kernel('TDHF', ref, singlet=False) + + +class DiamondShifted(unittest.TestCase): + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-hf-rev' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + + kpt = np.asarray([0.3721, 0.2077, 0.1415]) + + mf = scf.RHF(cell, kpt).rs_density_fit(auxbasis='weigend').run() + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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(nstates=self.nstates, **kwargs).run() + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + + def test_tda_singlet(self): + ref = [12.7166510188, 13.5460934688] + self.kernel('TDA', ref) + + def test_tda_triplet(self): + ref = [8.7359078688, 9.2565904604] + self.kernel('TDA' , ref, singlet=False) + + def test_tdhf_singlet(self): + ref = [12.6104811733, 13.4160717812] + self.kernel('TDHF', ref) + + def test_tdhf_triplet(self): + ref = [3.8940277713, 7.9448161493] + self.kernel('TDHF', ref, singlet=False) + + +class WaterBigBox(unittest.TestCase): + ''' Match molecular CIS ''' - cell.basis = {'C': [[0, (0.8, 1.0)], - [1, (1.0, 1.0)]]} - # cell.basis = 'gth-dzvp' - cell.pseudo = 'gth-pade' - cell.a = ''' - 0.000000000, 3.370137329, 3.370137329 - 3.370137329, 0.000000000, 3.370137329 - 3.370137329, 3.370137329, 0.000000000''' - cell.unit = 'B' - cell.verbose = 0 - cell.build() - -def tearDownModule(): - global cell - del cell - -class KnownValues(unittest.TestCase): - def test_tda_gamma_point(self): - mf = RHF(cell).run(conv_tol=1e-10) - td_model = tdscf.TDA(mf) - td_model.kernel() - e1 = td_model.e - - kmf = KRHF(cell, cell.make_kpts([1, 1, 1])).run(conv_tol=1e-10) - td_model = tdscf.KTDA(kmf) - td_model.kernel() - e2 = td_model.e - self.assertAlmostEqual(abs(e1-e2).max(), 0, 6) - self.assertAlmostEqual(abs(e1 - 1.0329858545904074).max(), 0, 5) - - def test_tdhf_gamma_point(self): - mf = RHF(cell).run(conv_tol=1e-10) - td_model = tdscf.TDHF(mf) - td_model.kernel() - e1 = td_model.e - - kmf = KRHF(cell, cell.make_kpts([1, 1, 1])).run(conv_tol=1e-10) - td_model = tdscf.KTDHF(kmf) - td_model.kernel() - e2 = td_model.e - self.assertAlmostEqual(abs(e1-e2).max(), 0, 5) - self.assertAlmostEqual(abs(e1 - 1.0301736485136344).max(), 0, 5) - -if __name__ == '__main__': - print("Tests for pbc.tdscf.rhf") + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = ''' + O 0.00000 0.00000 0.11779 + H 0.00000 0.75545 -0.47116 + H 0.00000 -0.75545 -0.47116 + ''' + cell.a = np.eye(3) * 15 + cell.basis = 'sto-3g' + cell.build() + mf = scf.RHF(cell).rs_density_fit(auxbasis='weigend') + mf.with_df.omega = 0.1 + mf.kernel() + cls.cell = cell + cls.mf = mf + + mol = molgto.Mole() + for key in ['verbose','output','atom','basis']: + setattr(mol, key, getattr(cell, key)) + mol.build() + molmf = molscf.RHF(mol).density_fit(auxbasis=mf.with_df.auxbasis).run() + cls.mol = mol + cls.molmf = molmf + + cls.nstates = 5 # make sure first `nstates_test` states are converged + cls.nstates_test = 2 + + @classmethod + def tearDownClass(cls): + cls.cell.stdout.close() + cls.mol.stdout.close() + del cls.cell, cls.mf + del cls.mol, cls.molmf + + def kernel(self, TD, **kwargs): + td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + moltd = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - + moltd.e[:self.nstates_test] * unitev).max(), 0, 2) + + def test_tda_singlet(self): + self.kernel('TDA') + + def test_tda_triplet(self): + self.kernel('TDA', singlet=False) + + def test_tdhf_singlet(self): + self.kernel('TDHF') + + def test_tdhf_triplet(self): + self.kernel('TDHF', singlet=False) + + +if __name__ == "__main__": + print("Full Tests for rhf-TDA and rhf-TDHF") unittest.main() diff --git a/pyscf/pbc/tdscf/test/test_rks.py b/pyscf/pbc/tdscf/test/test_rks.py new file mode 100644 index 0000000000..674649e208 --- /dev/null +++ b/pyscf/pbc/tdscf/test/test_rks.py @@ -0,0 +1,309 @@ +#!/usr/bin/env python +# Copyright 2021 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 unittest +import numpy as np +from pyscf import __config__ +from pyscf import gto as molgto, scf as molscf, tdscf as moltdscf +from pyscf.pbc import gto, scf, tdscf +from pyscf.data.nist import HARTREE2EV as unitev + + +class DiamondPBE(unittest.TestCase): + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-pbe' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + + xc = 'pbe' + mf = scf.RKS(cell).set(xc=xc).rs_density_fit(auxbasis='weigend').run() + + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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(nstates=self.nstates, **kwargs).run() + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + + def test_tda_singlet(self): + ref = [9.2625251659, 9.2625251659] + self.kernel('TDA', ref) + + def test_tda_triplet(self): + ref = [4.8364248190, 4.8364248190] + self.kernel('TDA', ref, singlet=False) + + def test_tddft_singlet(self): + ref = [8.8773512896, 8.8773512896] + self.kernel('TDDFT', ref) + + def test_tddft_triplet(self): + ref = [4.7694490556, 4.7694490556] + self.kernel('TDDFT', ref, singlet=False) + + +class DiamondPBEShifted(unittest.TestCase): + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-pbe' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + + kpt = np.asarray([0.3721, 0.2077, 0.1415]) + + xc = 'pbe' + mf = scf.RKS(cell, kpt).set(xc=xc).rs_density_fit(auxbasis='weigend').run() + + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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(nstates=self.nstates, **kwargs).run() + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + + def test_tda_singlet(self): + ref = [11.9664870288, 12.7605699008] + self.kernel('TDA', ref) + + def test_tda_triplet(self): + ref = [8.5705015296, 9.3030273411] + self.kernel('TDA', ref, singlet=False) + + def test_tddft_singlet(self): + ref = [11.8322851619, 12.6207316217] + self.kernel('TDDFT', ref) + + def test_tddft_triplet(self): + ref = [8.4227532516, 9.1695913993] + self.kernel('TDDFT', ref, singlet=False) + + +class WaterBigBoxPBE(unittest.TestCase): + ''' Match molecular CIS + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = ''' + O 0.00000 0.00000 0.11779 + H 0.00000 0.75545 -0.47116 + H 0.00000 -0.75545 -0.47116 + ''' + cell.a = np.eye(3) * 15 + cell.basis = 'sto-3g' + cell.build() + + xc = 'pbe' + mf = scf.RKS(cell).set(xc=xc).rs_density_fit(auxbasis='weigend') + mf.with_df.omega = 0.1 + mf.kernel() + cls.cell = cell + cls.mf = mf + + mol = molgto.Mole() + for key in ['verbose','output','atom','basis']: + setattr(mol, key, getattr(cell, key)) + mol.build() + molmf = molscf.RKS(mol).set(xc=xc).density_fit(auxbasis=mf.with_df.auxbasis).run() + cls.mol = mol + cls.molmf = molmf + + cls.nstates = 5 # make sure first `nstates_test` states are converged + cls.nstates_test = 2 + + @classmethod + def tearDownClass(cls): + cls.cell.stdout.close() + cls.mol.stdout.close() + del cls.cell, cls.mf + del cls.mol, cls.molmf + + def kernel(self, TD, **kwargs): + td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + moltd = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + # larger tol (0.1 eV) potentially due to DFT integration error + self.assertTrue((abs(td.e[:self.nstates_test] * unitev - + moltd.e[:self.nstates_test] * unitev).max() < 0.1)) + + def test_tda_singlet(self): + self.kernel('TDA') + + def test_tda_triplet(self): + self.kernel('TDA', singlet=False) + + def test_tddft_singlet(self): + self.kernel('TDDFT') + + def test_tddft_triplet(self): + self.kernel('TDDFT', singlet=False) + + +class DiamondPBE0(unittest.TestCase): + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-pbe' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + + xc = 'pbe0' + mf = scf.RKS(cell).set(xc=xc).rs_density_fit(auxbasis='weigend').run() + + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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(nstates=self.nstates, **kwargs).run() + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + + def test_tda_singlet(self): + ref = [9.6202884134, 9.6202884134] + self.kernel('TDA', ref) + + def test_tda_triplet(self): + ref = [5.1659937745, 5.1659937745] + self.kernel('TDA', ref, singlet=False) + + def test_tddft_singlet(self): + ref = [9.2585491075, 9.2585491075] + self.kernel('TDDFT', ref) + + def test_tddft_triplet(self): + ref = [4.5570089807, 4.5570089807] + self.kernel('TDDFT', ref, singlet=False) + + +class WaterBigBoxPBE0(unittest.TestCase): + ''' Match molecular CIS + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = ''' + O 0.00000 0.00000 0.11779 + H 0.00000 0.75545 -0.47116 + H 0.00000 -0.75545 -0.47116 + ''' + cell.a = np.eye(3) * 15 + cell.basis = 'sto-3g' + cell.build() + + xc = 'pbe0' + mf = scf.RKS(cell).set(xc=xc).rs_density_fit(auxbasis='weigend') + mf.with_df.omega = 0.1 + mf.kernel() + cls.cell = cell + cls.mf = mf + + mol = molgto.Mole() + for key in ['verbose','output','atom','basis']: + setattr(mol, key, getattr(cell, key)) + mol.build() + molmf = molscf.RKS(mol).set(xc=xc).density_fit(auxbasis=mf.with_df.auxbasis).run() + cls.mol = mol + cls.molmf = molmf + + cls.nstates = 5 # make sure first `nstates_test` states are converged + cls.nstates_test = 2 + + @classmethod + def tearDownClass(cls): + cls.cell.stdout.close() + cls.mol.stdout.close() + del cls.cell, cls.mf + del cls.mol, cls.molmf + + def kernel(self, TD, **kwargs): + td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + moltd = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + # larger tol (0.1 eV) potentially due to DFT integration error + self.assertTrue((abs(td.e[:self.nstates_test] * unitev - + moltd.e[:self.nstates_test] * unitev).max() < 0.1)) + + def test_tda_singlet(self): + self.kernel('TDA') + + def test_tda_triplet(self): + self.kernel('TDA', singlet=False) + + def test_tddft_singlet(self): + self.kernel('TDDFT') + + def test_tddft_triplet(self): + self.kernel('TDDFT', singlet=False) + + +if __name__ == "__main__": + print("Full Tests for rks-TDA and rks-TDDFT") + unittest.main() diff --git a/pyscf/pbc/tdscf/test/test_uhf.py b/pyscf/pbc/tdscf/test/test_uhf.py new file mode 100644 index 0000000000..39d4d8042a --- /dev/null +++ b/pyscf/pbc/tdscf/test/test_uhf.py @@ -0,0 +1,126 @@ +#!/usr/bin/env python +# Copyright 2021 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 unittest +import numpy as np +from pyscf import __config__ +from pyscf import gto as molgto, scf as molscf, tdscf as moltdscf +from pyscf.pbc import gto, scf, tdscf +from pyscf.data.nist import HARTREE2EV as unitev + + +class Diamond(unittest.TestCase): + ''' Reproduce RHF-TDSCF results + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-hf-rev' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + mf = scf.UHF(cell).rs_density_fit(auxbasis='weigend').run() + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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(nstates=self.nstates, **kwargs).run() + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + + def test_tda(self): + # same as lowest roots in test_rhf/Diamond -> test_tda_singlet/triplet + ref = [4.7209460258, 5.6500725912] + self.kernel('TDA', ref) + + def test_tdhf(self): + # same as lowest roots in test_rhf/Diamond -> test_tdhf_singlet/triplet + ref = [3.0396052214, 3.0396052214] + self.kernel('TDHF', ref) + + +class WaterBigBox(unittest.TestCase): + ''' Match molecular CIS + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = ''' + O 0.00000 0.00000 0.11779 + H 0.00000 0.75545 -0.47116 + H 0.00000 -0.75545 -0.47116 + ''' + cell.spin = 2 + cell.a = np.eye(3) * 15 + cell.basis = 'sto-3g' + cell.build() + mf = scf.UHF(cell).rs_density_fit(auxbasis='weigend') + mf.with_df.omega = 0.1 + mf.kernel() + cls.cell = cell + cls.mf = mf + + mol = molgto.Mole() + for key in ['verbose','output','atom','spin','basis']: + setattr(mol, key, getattr(cell, key)) + mol.build() + molmf = molscf.UHF(mol).density_fit(auxbasis=mf.with_df.auxbasis).run() + cls.mol = mol + cls.molmf = molmf + + cls.nstates = 5 # make sure first `nstates_test` states are converged + cls.nstates_test = 2 + + @classmethod + def tearDownClass(cls): + cls.cell.stdout.close() + cls.mol.stdout.close() + del cls.cell, cls.mf + del cls.mol, cls.molmf + + def kernel(self, TD, **kwargs): + td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + moltd = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - + moltd.e[:self.nstates_test] * unitev).max(), 0, 2) + + def test_tda(self): + self.kernel('TDA') + + def test_tdhf(self): + self.kernel('TDHF') + + +if __name__ == "__main__": + print("Full Tests for uhf-TDA and uhf-TDHF") + unittest.main() diff --git a/pyscf/pbc/tdscf/test/test_uks.py b/pyscf/pbc/tdscf/test/test_uks.py new file mode 100644 index 0000000000..19caa353d4 --- /dev/null +++ b/pyscf/pbc/tdscf/test/test_uks.py @@ -0,0 +1,230 @@ +#!/usr/bin/env python +# Copyright 2021 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 unittest +import numpy as np +from pyscf import __config__ +from pyscf import gto as molgto, scf as molscf, tdscf as moltdscf +from pyscf.pbc import gto, scf, tdscf +from pyscf.data.nist import HARTREE2EV as unitev + + +class DiamondPBE(unittest.TestCase): + ''' Reproduce RKS-TDSCF results + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-hf-rev' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + + xc = 'pbe' + mf = scf.UKS(cell).set(xc=xc).rs_density_fit(auxbasis='weigend').run() + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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(nstates=self.nstates, **kwargs).run() + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + + def test_tda(self): + ref = [4.8137763755, 4.8137763755] + self.kernel('TDA', ref) + + def test_tdhf(self): + ref = [4.7455726874, 4.7455726874] + self.kernel('TDDFT', ref) + + +class WaterBigBoxPBE(unittest.TestCase): + ''' Match molecular CIS + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = ''' + O 0.00000 0.00000 0.11779 + H 0.00000 0.75545 -0.47116 + H 0.00000 -0.75545 -0.47116 + ''' + cell.spin = 2 + cell.a = np.eye(3) * 15 + cell.basis = 'sto-3g' + cell.build() + + xc = 'pbe' + mf = scf.UKS(cell).set(xc=xc).rs_density_fit(auxbasis='weigend') + mf.with_df.omega = 0.1 + mf.kernel() + cls.cell = cell + cls.mf = mf + + mol = molgto.Mole() + for key in ['verbose','output','atom','spin','basis']: + setattr(mol, key, getattr(cell, key)) + mol.build() + molmf = molscf.UKS(mol).set(xc=xc).density_fit(auxbasis=mf.with_df.auxbasis).run() + cls.mol = mol + cls.molmf = molmf + + cls.nstates = 5 # make sure first `nstates_test` states are converged + cls.nstates_test = 2 + + @classmethod + def tearDownClass(cls): + cls.cell.stdout.close() + cls.mol.stdout.close() + del cls.cell, cls.mf + del cls.mol, cls.molmf + + def kernel(self, TD, **kwargs): + td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + moltd = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + self.assertTrue(abs(td.e[:self.nstates_test] * unitev - + moltd.e[:self.nstates_test] * unitev).max() < 0.1) + + def test_tda(self): + self.kernel('TDA') + + def test_tdhf(self): + self.kernel('TDDFT') + + +class DiamondPBE0(unittest.TestCase): + ''' Reproduce RKS-TDSCF results + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = 'C 0 0 0; C 0.8925000000 0.8925000000 0.8925000000' + cell.a = ''' + 1.7850000000 1.7850000000 0.0000000000 + 0.0000000000 1.7850000000 1.7850000000 + 1.7850000000 0.0000000000 1.7850000000 + ''' + cell.pseudo = 'gth-hf-rev' + cell.basis = {'C': [[0, (0.8, 1.0)], + [1, (1.0, 1.0)]]} + cell.precision = 1e-10 + cell.build() + + xc = 'pbe0' + mf = scf.UKS(cell).set(xc=xc).rs_density_fit(auxbasis='weigend').run() + mf.grids.level = 4 + cls.cell = cell + cls.mf = mf + + cls.nstates = 5 # 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(nstates=self.nstates, **kwargs).run() + self.assertAlmostEqual(abs(td.e[:self.nstates_test] * unitev - ref).max(), 0, 4) + + def test_tda(self): + ref = [5.1435121896, 5.1435121896] + self.kernel('TDA', ref) + + def test_tdhf(self): + ref = [4.5331029147, 4.5331029147] + self.kernel('TDDFT', ref) + + +class WaterBigBoxPBE0(unittest.TestCase): + ''' Match molecular CIS + ''' + @classmethod + def setUpClass(cls): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = ''' + O 0.00000 0.00000 0.11779 + H 0.00000 0.75545 -0.47116 + H 0.00000 -0.75545 -0.47116 + ''' + cell.spin = 2 + cell.a = np.eye(3) * 15 + cell.basis = 'sto-3g' + cell.build() + + xc = 'pbe0' + mf = scf.UKS(cell).set(xc=xc).rs_density_fit(auxbasis='weigend') + mf.with_df.omega = 0.1 + mf.kernel() + cls.cell = cell + cls.mf = mf + + mol = molgto.Mole() + for key in ['verbose','output','atom','spin','basis']: + setattr(mol, key, getattr(cell, key)) + mol.build() + molmf = molscf.UKS(mol).set(xc=xc).density_fit(auxbasis=mf.with_df.auxbasis).run() + cls.mol = mol + cls.molmf = molmf + + cls.nstates = 5 # make sure first `nstates_test` states are converged + cls.nstates_test = 2 + + @classmethod + def tearDownClass(cls): + cls.cell.stdout.close() + cls.mol.stdout.close() + del cls.cell, cls.mf + del cls.mol, cls.molmf + + def kernel(self, TD, **kwargs): + td = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + moltd = getattr(self.mf, TD)().set(nstates=self.nstates, **kwargs).run() + self.assertTrue(abs(td.e[:self.nstates_test] * unitev - + moltd.e[:self.nstates_test] * unitev).max() < 0.1) + + def test_tda(self): + self.kernel('TDA') + + def test_tdhf(self): + self.kernel('TDDFT') + + +if __name__ == "__main__": + print("Full Tests for uks-TDA and uks-TDDFT") + unittest.main() diff --git a/pyscf/pbc/tdscf/uks.py b/pyscf/pbc/tdscf/uks.py index 05a518ea06..c8b1fc7e34 100644 --- a/pyscf/pbc/tdscf/uks.py +++ b/pyscf/pbc/tdscf/uks.py @@ -28,6 +28,7 @@ class CasidaTDDFT(TDA): _gen_vind = uks.TDDFTNoHybrid.gen_vind gen_vind = TDA.gen_vind + kernel = uks.TDDFTNoHybrid.kernel TDDFTNoHybrid = CasidaTDDFT @@ -45,4 +46,3 @@ def tddft(mf): dft.uks.UKS.TDDFT = tddft #dft.rks.RKS.dTDA = lib.class_as_method(dTDA) #dft.rks.RKS.dRPA = lib.class_as_method(dRPA) - diff --git a/pyscf/tdscf/dhf.py b/pyscf/tdscf/dhf.py index 6d9946e39f..fe09ef5e37 100644 --- a/pyscf/tdscf/dhf.py +++ b/pyscf/tdscf/dhf.py @@ -411,14 +411,17 @@ def init_guess(self, mf, nstates=None, wfnsym=None): occidx = n2c + numpy.where(mo_occ[n2c:] == 1)[0] viridx = n2c + numpy.where(mo_occ[n2c:] == 0)[0] e_ia = mo_energy[viridx] - mo_energy[occidx,None] + # make degenerate excitations equal for later selection by energy + e_ia = numpy.ceil(e_ia / self.deg_eia_thresh) * self.deg_eia_thresh + e_ia_uniq = numpy.unique(e_ia) e_ia_max = e_ia.max() nov = e_ia.size nstates = min(nstates, nov) + nstates_thresh = min(nstates, e_ia_uniq.size) e_ia = e_ia.ravel() - e_threshold = min(e_ia_max, e_ia[numpy.argsort(e_ia)[nstates-1]]) - # Handle degeneracy, include all degenerated states in initial guess - e_threshold += 1e-6 + e_threshold = min(e_ia_max, e_ia_uniq[numpy.argsort(e_ia_uniq)[nstates_thresh-1]]) + e_threshold += self.deg_eia_thresh idx = numpy.where(e_ia <= e_threshold)[0] x0 = numpy.zeros((idx.size, nov)) @@ -442,13 +445,14 @@ def kernel(self, x0=None, nstates=None): vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) - def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + if x0 is None: + x0 = self.init_guess(self._scf, self.nstates) + x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] + # FIXME: Is it correct to call davidson1 for complex integrals? self.converged, self.e, x1 = \ lib.davidson1(vind, x0, precond, @@ -549,8 +553,6 @@ def kernel(self, x0=None, nstates=None): vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) def pickeig(w, v, nroots, envs): realidx = numpy.where((abs(w.imag) < REAL_EIG_THRESHOLD) & @@ -559,6 +561,10 @@ def pickeig(w, v, nroots, envs): return lib.linalg_helper._eigs_cmplx2real(w, v, realidx, real_eigenvectors=False) + if x0 is None: + x0 = self.init_guess(self._scf, self.nstates) + x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] + self.converged, w, x1 = \ lib.davidson_nosym1(vind, x0, precond, tol=self.conv_tol, diff --git a/pyscf/tdscf/ghf.py b/pyscf/tdscf/ghf.py index 00cc0a672b..1cf0424e02 100644 --- a/pyscf/tdscf/ghf.py +++ b/pyscf/tdscf/ghf.py @@ -387,6 +387,8 @@ def init_guess(self, mf, nstates=None, wfnsym=None): occidx = numpy.where(mo_occ==1)[0] viridx = numpy.where(mo_occ==0)[0] e_ia = mo_energy[viridx] - mo_energy[occidx,None] + # make degenerate excitations equal for later selection by energy + e_ia = numpy.ceil(e_ia / self.deg_eia_thresh) * self.deg_eia_thresh e_ia_max = e_ia.max() if wfnsym is not None and mf.mol.symmetry: @@ -397,12 +399,14 @@ def init_guess(self, mf, nstates=None, wfnsym=None): 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_uniq = numpy.unique(e_ia) + nov = e_ia.size nstates = min(nstates, nov) + nstates_thresh = min(nstates, e_ia_uniq.size) e_ia = e_ia.ravel() - e_threshold = min(e_ia_max, e_ia[numpy.argsort(e_ia)[nstates-1]]) - # Handle degeneracy, include all degenerated states in initial guess - e_threshold += 1e-6 + e_threshold = min(e_ia_max, e_ia_uniq[numpy.argsort(e_ia_uniq)[nstates_thresh-1]]) + e_threshold += self.deg_eia_thresh idx = numpy.where(e_ia <= e_threshold)[0] x0 = numpy.zeros((idx.size, nov)) @@ -426,13 +430,14 @@ def kernel(self, x0=None, nstates=None): vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) - def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + if x0 is None: + x0 = self.init_guess(self._scf, self.nstates) + x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] + # FIXME: Is it correct to call davidson1 for complex integrals self.converged, self.e, x1 = \ lib.davidson1(vind, x0, precond, @@ -562,8 +567,6 @@ def kernel(self, x0=None, nstates=None): vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) ensure_real = self._scf.mo_coeff.dtype == numpy.double def pickeig(w, v, nroots, envs): @@ -572,6 +575,10 @@ 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) + if x0 is None: + x0 = self.init_guess(self._scf, self.nstates) + x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] + self.converged, w, x1 = \ lib.davidson_nosym1(vind, x0, precond, tol=self.conv_tol, diff --git a/pyscf/tdscf/gks.py b/pyscf/tdscf/gks.py index 96bc0e77a4..7681abba27 100644 --- a/pyscf/tdscf/gks.py +++ b/pyscf/tdscf/gks.py @@ -116,13 +116,15 @@ def kernel(self, x0=None, nstates=None): vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + if x0 is None: + x0 = self.init_guess(self._scf, self.nstates) + x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] + self.converged, w2, x1 = \ lib.davidson1(vind, x0, precond, tol=self.conv_tol, diff --git a/pyscf/tdscf/rhf.py b/pyscf/tdscf/rhf.py index cf84135118..dadfc73b7f 100644 --- a/pyscf/tdscf/rhf.py +++ b/pyscf/tdscf/rhf.py @@ -24,6 +24,7 @@ from functools import reduce import numpy +import scipy.linalg from pyscf import lib from pyscf import gto from pyscf import scf @@ -49,7 +50,7 @@ def gen_tda_operation(mf, fock_ao=None, singlet=True, wfnsym=None): ''' mol = mf.mol mo_coeff = mf.mo_coeff - assert (mo_coeff.dtype == numpy.double) + # assert (mo_coeff.dtype == numpy.double) mo_energy = mf.mo_energy mo_occ = mf.mo_occ nao, nmo = mo_coeff.shape @@ -114,7 +115,7 @@ def get_ab(mf, mo_energy=None, mo_coeff=None, mo_occ=None): if mo_energy is None: mo_energy = mf.mo_energy if mo_coeff is None: mo_coeff = mf.mo_coeff if mo_occ is None: mo_occ = mf.mo_occ - assert (mo_coeff.dtype == numpy.double) + # assert (mo_coeff.dtype == numpy.double) mol = mf.mol nao, nmo = mo_coeff.shape @@ -646,6 +647,8 @@ class TDMixin(lib.StreamObject): 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) + # Threshold to handle degeneracy in init guess + deg_eia_thresh = getattr(__config__, 'tdscf_rhf_TDDFT_deg_eia_thresh', 1e-3) def __init__(self, mf): self.verbose = mf.verbose @@ -690,6 +693,7 @@ def dump_flags(self, verbose=None): log.info('nstates = %d singlet', self.nstates) else: log.info('nstates = %d triplet', self.nstates) + log.info('deg_eia_thresh = %.3e', self.deg_eia_thresh) log.info('wfnsym = %s', self.wfnsym) log.info('conv_tol = %g', self.conv_tol) log.info('eigh lindep = %g', self.lindep) @@ -729,6 +733,26 @@ def precond(x, e, x0): return x/diagd return precond + def trunc_workspace(self, vind, x0, nstates=None, pick=None): + log = logger.new_logger(self) + if nstates is None: nstates = self.nstates + if x0.shape[0] <= nstates: + return None, x0 + else: + heff = lib.einsum('xa,ya->xy', x0.conj(), vind(x0)) + e, u = scipy.linalg.eig(heff) # heff not necessarily Hermitian + e = e.real + order = numpy.argsort(e) + e = e[order] + u = u[:,order] + if callable(pick): + e, u = pick(e, u, None, None)[:2] + e = e[:nstates] + log.debug('truncating %d states to %d states with excitation energy (eV): %s', + x0.shape[0], e.size, e*27.211399) + x1 = numpy.dot(x0.T, u[:,:nstates]).T + return e, x1 + analyze = analyze get_nto = get_nto oscillator_strength = oscillator_strength @@ -793,6 +817,8 @@ def init_guess(self, mf, nstates=None, wfnsym=None): occidx = numpy.where(mo_occ==2)[0] viridx = numpy.where(mo_occ==0)[0] e_ia = mo_energy[viridx] - mo_energy[occidx,None] + # make degenerate excitations equal for later selection by energy + e_ia = numpy.ceil(e_ia / self.deg_eia_thresh) * self.deg_eia_thresh e_ia_max = e_ia.max() if wfnsym is not None and mf.mol.symmetry: @@ -803,12 +829,14 @@ def init_guess(self, mf, nstates=None, wfnsym=None): 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_uniq = numpy.unique(e_ia) + nov = e_ia.size nstates = min(nstates, nov) + nstates_thresh = min(nstates, e_ia_uniq.size) e_ia = e_ia.ravel() - e_threshold = min(e_ia_max, e_ia[numpy.argsort(e_ia)[nstates-1]]) - # Handle degeneracy, include all degenerated states in initial guess - e_threshold += 1e-6 + e_threshold = min(e_ia_max, e_ia_uniq[numpy.argsort(e_ia_uniq)[nstates_thresh-1]]) + e_threshold += self.deg_eia_thresh idx = numpy.where(e_ia <= e_threshold)[0] x0 = numpy.zeros((idx.size, nov)) @@ -832,13 +860,14 @@ def kernel(self, x0=None, nstates=None): vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) - def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + if x0 is None: + x0 = self.init_guess(self._scf, self.nstates) + x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] + self.converged, self.e, x1 = \ lib.davidson1(vind, x0, precond, tol=self.conv_tol, @@ -872,7 +901,7 @@ def gen_tdhf_operation(mf, fock_ao=None, singlet=True, wfnsym=None): ''' mol = mf.mol mo_coeff = mf.mo_coeff - assert (mo_coeff.dtype == numpy.double) + # assert (mo_coeff.dtype == numpy.double) mo_energy = mf.mo_energy mo_occ = mf.mo_occ nao, nmo = mo_coeff.shape @@ -987,8 +1016,14 @@ def kernel(self, x0=None, nstates=None): vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) + + # 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): @@ -997,8 +1032,11 @@ def pickeig(w, v, nroots, envs): # If the complex eigenvalue has small imaginary part, both the # real part and the imaginary part of the eigenvector can # approximately be used as the "real" eigen solutions. - return lib.linalg_helper._eigs_cmplx2real(w, v, realidx, - real_eigenvectors=True) + return lib.linalg_helper._eigs_cmplx2real(w, v, realidx, real_system) + + if x0 is None: + x0 = self.init_guess(self._scf, self.nstates) + x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] self.converged, w, x1 = \ lib.davidson_nosym1(vind, x0, precond, diff --git a/pyscf/tdscf/rks.py b/pyscf/tdscf/rks.py index 82486b6215..fc9860b8fd 100644 --- a/pyscf/tdscf/rks.py +++ b/pyscf/tdscf/rks.py @@ -120,13 +120,15 @@ def kernel(self, x0=None, nstates=None): vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + if x0 is None: + x0 = self.init_guess(self._scf, self.nstates) + x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] + self.converged, w2, x1 = \ lib.davidson1(vind, x0, precond, tol=self.conv_tol, diff --git a/pyscf/tdscf/test/test_proxy.py b/pyscf/tdscf/test/test_proxy.py index d182c5fe66..37cbb3a517 100644 --- a/pyscf/tdscf/test/test_proxy.py +++ b/pyscf/tdscf/test/test_proxy.py @@ -79,7 +79,8 @@ def test_class_frozen(self): model.kernel() mask_o, mask_v = tdhf_frozen_mask(model.eri, kind="o,v") testing.assert_allclose(model.e, self.td_model_rks.e, atol=1e-3) - assert_vectors_close(model.xy, numpy.array(self.td_model_rks.xy)[..., mask_o, :][..., mask_v], atol=1e-3) + assert_vectors_close(model.xy, + numpy.array(self.td_model_rks.xy)[..., mask_o, :][..., mask_v], atol=1e-3) except Exception: print("When testing class with frozen={} the following exception occurred:".format(repr(frozen))) @@ -193,7 +194,8 @@ def test_class_frozen(self): model.kernel() mask_o, mask_v = tdhf_frozen_mask(model.eri, kind="o,v") testing.assert_allclose(model.e, self.td_model_rhf.e, atol=1e-3) - assert_vectors_close(model.xy, numpy.array(self.td_model_rhf.xy)[..., mask_o, :][..., mask_v], atol=1e-3) + assert_vectors_close(model.xy, + numpy.array(self.td_model_rhf.xy)[..., mask_o, :][..., mask_v], atol=1e-3) except Exception: print("When testing class with frozen={} the following exception occurred:".format(repr(frozen))) diff --git a/pyscf/tdscf/test/test_rhf_slow.py b/pyscf/tdscf/test/test_rhf_slow.py index c104335215..44dcb0daa8 100644 --- a/pyscf/tdscf/test/test_rhf_slow.py +++ b/pyscf/tdscf/test/test_rhf_slow.py @@ -126,7 +126,8 @@ def test_class_frozen(self): model.kernel() mask_o, mask_v = tdhf_frozen_mask(model.eri, kind="o,v") testing.assert_allclose(model.e, self.td_model_rhf.e, atol=1e-3) - assert_vectors_close(model.xy, numpy.array(self.td_model_rhf.xy)[..., mask_o, :][..., mask_v], atol=1e-3) + assert_vectors_close(model.xy, + numpy.array(self.td_model_rhf.xy)[..., mask_o, :][..., mask_v], atol=1e-3) except Exception: print("When testing class with frozen={} the following exception occurred:".format(repr(frozen))) diff --git a/pyscf/tdscf/test/test_tdrhf.py b/pyscf/tdscf/test/test_tdrhf.py index 3d4d3621b6..243d897bb2 100644 --- a/pyscf/tdscf/test/test_tdrhf.py +++ b/pyscf/tdscf/test/test_tdrhf.py @@ -18,7 +18,7 @@ from pyscf import lib, gto, scf, tdscf def setUpModule(): - global mol, mf + global mol, mf, nstates mol = gto.Mole() mol.verbose = 7 mol.output = '/dev/null' @@ -28,6 +28,7 @@ def setUpModule(): mol.basis = '631g' mol.build() mf = scf.RHF(mol).run() + nstates = 5 # make sure first 3 states are converged def tearDownModule(): global mol, mf @@ -36,30 +37,30 @@ def tearDownModule(): class KnownValues(unittest.TestCase): def test_tda_singlet(self): - td = mf.TDA() + td = mf.TDA().set(nstates=nstates) e = td.kernel()[0] ref = [11.90276464, 11.90276464, 16.86036434] - self.assertAlmostEqual(abs(e * 27.2114 - ref).max(), 0, 5) + self.assertAlmostEqual(abs(e[:len(ref)] * 27.2114 - ref).max(), 0, 5) def test_tda_triplet(self): - td = mf.TDA() + td = mf.TDA().set(nstates=nstates) td.singlet = False e = td.kernel()[0] ref = [11.01747918, 11.01747918, 13.16955056] - self.assertAlmostEqual(abs(e * 27.2114 - ref).max(), 0, 5) + self.assertAlmostEqual(abs(e[:len(ref)] * 27.2114 - ref).max(), 0, 5) def test_tdhf_singlet(self): - td = mf.TDHF() + td = mf.TDHF().set(nstates=nstates) e = td.kernel()[0] ref = [11.83487199, 11.83487199, 16.66309285] - self.assertAlmostEqual(abs(e * 27.2114 - ref).max(), 0, 5) + self.assertAlmostEqual(abs(e[:len(ref)] * 27.2114 - ref).max(), 0, 5) def test_tdhf_triplet(self): - td = mf.TDHF() + td = mf.TDHF().set(nstates=nstates) td.singlet = False e = td.kernel()[0] ref = [10.8919234, 10.8919234, 12.63440705] - self.assertAlmostEqual(abs(e * 27.2114 - ref).max(), 0, 5) + self.assertAlmostEqual(abs(e[:len(ref)] * 27.2114 - ref).max(), 0, 5) if __name__ == "__main__": diff --git a/pyscf/tdscf/test/test_tdrks.py b/pyscf/tdscf/test/test_tdrks.py index a59421d3ad..fef61007cc 100644 --- a/pyscf/tdscf/test/test_tdrks.py +++ b/pyscf/tdscf/test/test_tdrks.py @@ -161,7 +161,7 @@ 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) + self.assertAlmostEqual(lib.fp(es), -39.74044291202006, 5) ref = [9.0139312, 9.0139312, 12.42444659] self.assertAlmostEqual(abs(es[:3] - ref).max(), 0, 4) @@ -197,7 +197,7 @@ def test_tda_rsh(self): e = mf.kernel() self.assertAlmostEqual(e, -1.14670613191817, 8) - e_td = mf.TDA().kernel()[0] + e_td = mf.TDA().set(nstates=5).kernel()[0] ref = [16.25021865, 27.93720198, 49.4665691] self.assertAlmostEqual(abs(e_td*nist.HARTREE2EV - ref).max(), 0, 4) @@ -286,7 +286,7 @@ def test_ab_mgga(self): def test_nto(self): mf = scf.RHF(mol).run() - td = rks.TDA(mf).run() + td = rks.TDA(mf).run(nstates=5) w, nto = td.get_nto(state=3) self.assertAlmostEqual(w[0], 0.98655300613468389, 7) self.assertAlmostEqual(lib.fp(w), 0.98625701534112464, 7) diff --git a/pyscf/tdscf/uhf.py b/pyscf/tdscf/uhf.py index d9e3d2c127..92aabf0258 100644 --- a/pyscf/tdscf/uhf.py +++ b/pyscf/tdscf/uhf.py @@ -619,6 +619,10 @@ def init_guess(self, mf, nstates=None, wfnsym=None): 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 + # make degenerate excitations equal for later selection by energy + e_ia_a = numpy.ceil(e_ia_a / self.deg_eia_thresh) * self.deg_eia_thresh + e_ia_b = numpy.ceil(e_ia_b / self.deg_eia_thresh) * self.deg_eia_thresh + e_ia_max = max(numpy.max(e_ia_a), numpy.max(e_ia_b)) if wfnsym is not None and mol.symmetry: if isinstance(wfnsym, str): @@ -631,12 +635,12 @@ def init_guess(self, mf, nstates=None, wfnsym=None): 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())) - e_ia_max = e_ia.max() + e_ia_uniq = numpy.unique(e_ia) nov = e_ia.size nstates = min(nstates, nov) - e_threshold = min(e_ia_max, e_ia[numpy.argsort(e_ia)[nstates-1]]) - # Handle degeneracy - e_threshold += 1e-6 + nstates_thresh = min(nstates, e_ia_uniq.size) + e_threshold = min(e_ia_max, e_ia_uniq[numpy.argsort(e_ia_uniq)[nstates_thresh-1]]) + e_threshold += self.deg_eia_thresh idx = numpy.where(e_ia <= e_threshold)[0] x0 = numpy.zeros((idx.size, nov)) @@ -658,13 +662,15 @@ def kernel(self, x0=None, nstates=None): vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + if x0 is None: + x0 = self.init_guess(self._scf, self.nstates) + x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] + self.converged, self.e, x1 = \ lib.davidson1(vind, x0, precond, tol=self.conv_tol, @@ -812,8 +818,6 @@ def kernel(self, x0=None, nstates=None): vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) # We only need positive eigenvalues def pickeig(w, v, nroots, envs): @@ -822,6 +826,10 @@ def pickeig(w, v, nroots, envs): return lib.linalg_helper._eigs_cmplx2real(w, v, realidx, real_eigenvectors=True) + if x0 is None: + x0 = self.init_guess(self._scf, self.nstates) + x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] + self.converged, w, x1 = \ lib.davidson_nosym1(vind, x0, precond, tol=self.conv_tol, diff --git a/pyscf/tdscf/uks.py b/pyscf/tdscf/uks.py index e2d92a62a4..d4301584f9 100644 --- a/pyscf/tdscf/uks.py +++ b/pyscf/tdscf/uks.py @@ -138,13 +138,14 @@ def kernel(self, x0=None, nstates=None): vind, hdiag = self.gen_vind(self._scf) precond = self.get_precond(hdiag) - if x0 is None: - x0 = self.init_guess(self._scf, self.nstates) - def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + if x0 is None: + x0 = self.init_guess(self._scf, self.nstates) + x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] + self.converged, w2, x1 = \ lib.davidson1(vind, x0, precond, tol=self.conv_tol,