diff --git a/examples/pprpa/01-pprpa_total_energy.py b/examples/pprpa/01-pprpa_total_energy.py deleted file mode 100644 index 000bb4a05..000000000 --- a/examples/pprpa/01-pprpa_total_energy.py +++ /dev/null @@ -1,42 +0,0 @@ -from pyscf import gto, dft - -# For ppRPA total energy of N-electron system -# mean field is also N-electron - -mol = gto.Mole() -mol.verbose = 5 -mol.atom = [ - ["O", (0.0, 0.0, 0.0)], - ["H", (0.0, -0.7571, 0.5861)], - ["H", (0.0, 0.7571, 0.5861)]] -mol.basis = 'def2-svp' -mol.build() - -# =====> Part I. Restricted ppRPA <===== -mf = dft.RKS(mol) -mf.xc = "b3lyp" -mf.kernel() - -from pyscf.pprpa.rpprpa_direct import RppRPADirect -pp = RppRPADirect(mf) -ec = pp.get_correlation() -etot, ehf, ec = pp.energy_tot() -print("H2O Hartree-Fock energy = %.8f" % ehf) -print("H2O ppRPA correlation energy = %.8f" % ec) -print("H2O ppRPA total energy = %.8f" % etot) - -# =====> Part II. Unrestricted ppRPA <===== -# unrestricted KS-DFT calculation as starting point of UppRPA -umf = dft.UKS(mol) -umf.xc = "b3lyp" -umf.kernel() - -# direct diagonalization, N6 scaling -from pyscf.pprpa.upprpa_direct import UppRPADirect -pp = UppRPADirect(umf) -ec = pp.get_correlation() -etot, ehf, ec = pp.energy_tot() -print("H2O Hartree-Fock energy = %.8f" % ehf) -print("H2O ppRPA correlation energy = %.8f" % ec) -print("H2O ppRPA total energy = %.8f" % etot) - diff --git a/examples/pprpa/02-pprpa_excitation_energy.py b/examples/pprpa/02-pprpa_excitation_energy.py deleted file mode 100644 index a94660e88..000000000 --- a/examples/pprpa/02-pprpa_excitation_energy.py +++ /dev/null @@ -1,64 +0,0 @@ -from pyscf import gto, dft - -# For ppRPA excitation energy of N-electron system in particle-particle channel -# mean field is (N-2)-electron - -mol = gto.Mole() -mol.verbose = 5 -mol.atom = [ - ["O", (0.0, 0.0, 0.0)], - ["H", (0.0, -0.7571, 0.5861)], - ["H", (0.0, 0.7571, 0.5861)]] -mol.basis = 'def2-svp' -# create a (N-2)-electron system for charged-neutral H2O -mol.charge = 2 -mol.build() - -# =====> Part I. Restricted ppRPA <===== -# restricted KS-DFT calculation as starting point of RppRPA -rmf = dft.RKS(mol) -rmf.xc = "b3lyp" -rmf.kernel() - -# direct diagonalization, N6 scaling -from pyscf.pprpa.rpprpa_direct import RppRPADirect -# ppRPA can be solved in an active space -pp = RppRPADirect(rmf, nocc_act=None, nvir_act=10) -# number of two-electron addition states to print -pp.pp_state = 10 -# solve for singlet states -pp.kernel("s") -# solve for triplet states -pp.kernel("t") -pp.analyze() - -# Davidson algorithm, N4 scaling -from pyscf.pprpa.rpprpa_davidson import RppRPADavidson -# ppRPA can be solved in an active space -pp = RppRPADavidson(rmf, nocc_act=3, nvir_act=None, nroot=10) -# solve for singlet states -pp.kernel("s") -# solve for triplet states -pp.kernel("t") -pp.analyze() - -# =====> Part II. Unrestricted ppRPA <===== -# unrestricted KS-DFT calculation as starting point of UppRPA -umf = dft.UKS(mol) -umf.xc = "b3lyp" -umf.kernel() - -# direct diagonalization, N6 scaling -from pyscf.pprpa.upprpa_direct import UppRPADirect -# ppRPA can be solved in an active space -pp = UppRPADirect(umf, nocc_act=None, nvir_act=10) -# number of two-electron addition states to print -pp.pp_state = 10 -# solve ppRPA in the (alpha alpha, alpha alpha) subspace -pp.kernel(subspace=['aa']) -# solve ppRPA in the (alpha beta, alpha beta) subspace -pp.kernel(subspace=['ab']) -# solve ppRPA in the (beta beta, beta beta) subspace -pp.kernel(subspace=['bb']) -pp.analyze() - diff --git a/examples/pprpa/03-hhrpa_excitation_energy.py b/examples/pprpa/03-hhrpa_excitation_energy.py deleted file mode 100644 index 409893c4b..000000000 --- a/examples/pprpa/03-hhrpa_excitation_energy.py +++ /dev/null @@ -1,63 +0,0 @@ -from pyscf import gto, dft - -# For ppRPA excitation energy of N-electron system in hole-hole channel -# mean field is (N+2)-electron - -mol = gto.Mole() -mol.verbose = 5 -mol.atom = [ - ["O", (0.0, 0.0, 0.0)], - ["H", (0.0, -0.7571, 0.5861)], - ["H", (0.0, 0.7571, 0.5861)]] -mol.basis = 'def2-svp' -# create a (N+2)-electron system for charged-neutral H2O -mol.charge = -2 -mol.build() - -# =====> Part I. Restricted ppRPA <===== -mf = dft.RKS(mol) -mf.xc = "b3lyp" -mf.kernel() - -# direct diagonalization, N6 scaling -# ppRPA can be solved in an active space -from pyscf.pprpa.rpprpa_direct import RppRPADirect -pp = RppRPADirect(mf, nvir_act=10, nelec="n+2") -# number of two-electron removal states to print -pp.hh_state = 10 -# solve for singlet states -pp.kernel("s") -# solve for triplet states -pp.kernel("t") -pp.analyze() - -# Davidson algorithm, N4 scaling -# ppRPA can be solved in an active space -from pyscf.pprpa.rpprpa_davidson import RppRPADavidson -pp = RppRPADavidson(mf, nvir_act=10, channel="hh") -# solve for singlet states -pp.kernel("s") -# solve for triplet states -pp.kernel("t") -pp.analyze() - -# =====> Part II. Unrestricted ppRPA <===== -# unrestricted KS-DFT calculation as starting point of UppRPA -umf = dft.UKS(mol) -umf.xc = "b3lyp" -umf.kernel() - -# direct diagonalization, N6 scaling -from pyscf.pprpa.upprpa_direct import UppRPADirect -# ppRPA can be solved in an active space -pp = UppRPADirect(umf, nvir_act=10, nelec='n+2') -# number of two-electron addition states to print -pp.pp_state = 10 -# solve ppRPA in the (alpha alpha, alpha alpha) subspace -pp.kernel(subspace=['aa']) -# solve ppRPA in the (alpha beta, alpha beta) subspace -pp.kernel(subspace=['ab']) -# solve ppRPA in the (beta beta, beta beta) subspace -pp.kernel(subspace=['bb']) -pp.analyze() - diff --git a/examples/pprpa/04-gamma_pprpa_excitation_energy.py b/examples/pprpa/04-gamma_pprpa_excitation_energy.py deleted file mode 100644 index bf37ae854..000000000 --- a/examples/pprpa/04-gamma_pprpa_excitation_energy.py +++ /dev/null @@ -1,166 +0,0 @@ -import os - -from pyscf.pbc import df, dft, gto, lib - -# For ppRPA excitation energy of N-electron system in particle-particle channel -# mean field is (N-2)-electron - -# carbon vacancy in diamond -# see Table.1 in https://doi.org/10.1021/acs.jctc.4c00829 -cell = gto.Cell() -cell.build(unit='angstrom', - a=[[7.136000, 0.000000, 0.000000], - [0.000000, 7.136000, 0.000000], - [0.000000, 0.000000, 7.136000]], - atom=[ - ["C", (0.001233209267, 0.001233209267, 1.777383281725)], - ["C", (0.012225089066, 1.794731051862, 3.561871956432)], - ["C", (1.782392880032, 1.782392880032, 5.362763822515)], - ["C", (1.780552680464, -0.013167298675, -0.008776569968)], - ["C", (3.557268948138, 0.012225089066, 1.790128043568)], - ["C", (5.339774910934, 1.794731051862, 1.790128043568)], - ["C", (-0.013167298675, 1.780552680464, -0.008776569968)], - ["C", (1.782392880032, 3.569607119968, -0.010763822515)], - ["C", (1.794731051862, 5.339774910934, 1.790128043568)], - ["C", (2.671194343578, 0.900046784093, 0.881569117555)], - ["C", (0.887735886768, 0.887735886768, 6.244403380954)], - ["C", (0.900046784093, 2.680805656422, 4.470430882445)], - ["C", (2.680805656422, 4.451953215907, 0.881569117555)], - ["C", (0.895786995277, 6.238213500378, 0.896853305635)], - ["C", (0.930821705350, 0.930821705350, 2.673641624787)], - ["C", (4.421178294650, 0.930821705350, 2.678358375213)], - ["C", (0.900046784093, 2.671194343578, 0.881569117555)], - ["C", (6.238213500378, 0.895786995277, 0.896853305635)], - ["C", (2.680805656422, 0.900046784093, 4.470430882445)], - ["C", (4.451953215907, 2.680805656422, 0.881569117555)], - ["C", (0.930821705350, 4.421178294650, 2.678358375213)], - ["C", (1.794731051862, 0.012225089066, 3.561871956432)], - ["C", (0.012225089066, 3.557268948138, 1.790128043568)], - ["C", (3.569607119968, 1.782392880032, -0.010763822515)], - ["C", (1.736746319267, 1.736746319267, 1.671367479693)], - ["C", (5.351404126874, 0.000595873126, 0.004129648157)], - ["C", (0.000595873126, 5.351404126874, 0.004129648157)], - ["C", (2.676000000000, 2.676000000000, 6.244000000000)], - ["C", (6.244000000000, 2.676000000000, 2.676000000000)], - ["C", (2.676000000000, 6.244000000000, 2.676000000000)], - ["C", (0.000595873126, 0.000595873126, 5.347870351843)], - ["C", (0.001233209267, 5.350766790733, 3.574616718275)], - ["C", (1.780552680464, 5.365167298675, 5.360776569968)], - ["C", (3.571447319536, -0.013167298675, 5.360776569968)], - ["C", (5.365167298675, 1.780552680464, 5.360776569968)], - ["C", (5.365167298675, 3.571447319536, -0.008776569968)], - ["C", (5.350766790733, 5.350766790733, 1.777383281725)], - ["C", (4.464264113232, 0.887735886768, 6.243596619046)], - ["C", (4.451953215907, 2.671194343578, 4.470430882445)], - ["C", (2.671194343578, 4.451953215907, 4.470430882445)], - ["C", (0.895786995277, 6.249786499622, 4.455146694365)], - ["C", (4.421178294650, 4.421178294650, 2.673641624787)], - ["C", (6.249786499622, 4.456213004723, 0.896853305635)], - ["C", (0.887735886768, 4.464264113232, 6.243596619046)], - ["C", (6.249786499622, 0.895786995277, 4.455146694365)], - ["C", (4.456213004723, 6.249786499622, 0.896853305635)], - ["C", (5.350766790733, 0.001233209267, 3.574616718275)], - ["C", (-0.013167298675, 3.571447319536, 5.360776569968)], - ["C", (3.571447319536, 5.365167298675, -0.008776569968)], - ["C", (3.615253680733, 1.736746319267, 3.680632520307)], - ["C", (3.615253680733, 3.615253680733, 1.671367479693)], - ["C", (6.244000000000, 2.676000000000, 6.244000000000)], - ["C", (6.244000000000, 6.244000000000, 2.676000000000)], - ["C", (1.736746319267, 3.615253680733, 3.680632520307)], - ["C", (2.676000000000, 6.244000000000, 6.244000000000)], - ["C", (3.557268948138, 5.339774910934, 3.561871956432)], - ["C", (5.351404126874, 5.351404126874, 5.347870351843)], - ["C", (3.569607119968, 3.569607119968, 5.362763822515)], - ["C", (5.339774910934, 3.557268948138, 3.561871956432)], - ["C", (4.464264113232, 4.464264113232, 6.244403380954)], - ["C", (4.456213004723, 6.238213500378, 4.455146694365)], - ["C", (6.238213500378, 4.456213004723, 4.455146694365)], - ["C", (6.244000000000, 6.244000000000, 6.244000000000)], - ], - dimension=3, - max_memory=90000, - verbose=5, - basis='cc-pvdz', - # create a (N-2)-electron system - charge=2, - precision=1e-12) - -gdf = df.RSDF(cell) -gdf.auxbasis = "cc-pvdz-ri" -gdf_fname = 'gdf_ints.h5' -gdf._cderi_to_save = gdf_fname -if not os.path.isfile(gdf_fname): - gdf.build() - -# =====> Part I. Restricted ppRPA <===== -# After SCF, PySCF might fail in Makov-Payne correction -# save chkfile to restart -chkfname = 'scf.chk' -if os.path.isfile(chkfname): - kmf = dft.RKS(cell).rs_density_fit() - kmf.xc = "b3lyp" - kmf.exxdiv = None - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - data = lib.chkfile.load(chkfname, 'scf') - kmf.__dict__.update(data) -else: - kmf = dft.RKS(cell).rs_density_fit() - kmf.xc = "b3lyp" - kmf.exxdiv = None - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - kmf.chkfile = chkfname - kmf.kernel() - -# direct diagonalization, N6 scaling -# ppRPA can be solved in an active space -from pyscf.pprpa.rpprpa_direct import RppRPADirect -pp = RppRPADirect(kmf, nocc_act=50, nvir_act=50) -# number of two-electron addition states to print -pp.pp_state = 50 -# solve for singlet states -pp.kernel("s") -# solve for triplet states -pp.kernel("t") -pp.analyze() - -# Davidson algorithm, N4 scaling -# ppRPA can be solved in an active space -from pyscf.pprpa.rpprpa_davidson import RppRPADavidson -pp = RppRPADavidson(kmf, nocc_act=50, nvir_act=50, nroot=50) -# solve for singlet states -pp.kernel("s") -# solve for triplet states -pp.kernel("t") -pp.analyze() - -# =====> Part II. Unrestricted ppRPA <===== -chkfname = 'uscf.chk' -if os.path.isfile(chkfname): - kmf = dft.UKS(cell).rs_density_fit() - kmf.xc = "b3lyp" - kmf.exxdiv = None - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - data = lib.chkfile.load(chkfname, 'scf') - kmf.__dict__.update(data) -else: - kmf = dft.UKS(cell).rs_density_fit() - kmf.xc = "b3lyp" - kmf.exxdiv = None - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - kmf.chkfile = chkfname - kmf.kernel() - -# direct diagonalization, N6 scaling -# ppRPA can be solved in an active space -from pyscf.pprpa.upprpa_direct import UppRPADirect -pp = UppRPADirect(kmf, nocc_act=50, nvir_act=50) -# number of two-electron addition states to print -pp.pp_state = 50 -# solve ppRPA -pp.kernel() -pp.analyze() - diff --git a/examples/pprpa/05-gamma_hhrpa_excitation_energy.py b/examples/pprpa/05-gamma_hhrpa_excitation_energy.py deleted file mode 100644 index 20d80fac1..000000000 --- a/examples/pprpa/05-gamma_hhrpa_excitation_energy.py +++ /dev/null @@ -1,156 +0,0 @@ -import os - -from pyscf.pbc import df, dft, gto, lib - -# For ppRPA excitation energy of N-electron system in hole-hole channel -# mean field is (N+2)-electron - -# nitrogen vacancy in diamond (NV-) -# see in Table.1 https://doi.org/10.1021/acs.jpclett.4c00184 -cell = gto.Cell() -cell.build(unit='angstrom', - a=[[7.136000, 0.000000, 0.000000], - [0.000000, 7.136000, 0.000000], - [0.000000, 0.000000, 7.136000]], - atom=[ - ["C", (0.906073512160, 2.664323563008, 2.694339584224)], - ["C", (0.894090626784, 2.684866651264, 6.241909316128)], - ["C", (0.899084378176, 6.235829736704, 2.679298965664)], - ["C", (0.892434989152, 6.244294788160, 6.243564953760)], - ["C", (4.459048614208, 2.719457426496, 2.676951150304)], - ["C", (4.441660180288, 2.664323563008, 6.229926430752)], - ["C", (4.460040375488, 6.254734827520, 2.675959389024)], - ["C", (4.456700798848, 6.235829736704, 6.236915564736)], - ["C", (0.000121369088, 0.000121369088, -0.000121369088)], - ["C", (-0.010102670688, -0.010102670688, 3.575517319296)], - ["C", (-0.010102670688, 3.560482680704, 0.010102670688)], - ["C", (0.008425625056, 3.562468772224, 3.573531227776)], - ["C", (3.560482680704, -0.010102670688, 0.010102670688)], - ["C", (3.562468772224, 0.008425625056, 3.573531227776)], - ["C", (3.562468772224, 3.562468772224, -0.008425625056)], - ["C", (0.892434989152, 0.892434989152, 0.891705154752)], - ["C", (0.894090626784, 0.894090626784, 4.451133113248)], - ["C", (0.899084378176, 4.456700798848, 0.900170206208)], - ["C", (0.906073512160, 4.441660180288, 4.471676201504)], - ["C", (4.456700798848, 0.899084378176, 0.900170206208)], - ["C", (4.441660180288, 0.906073512160, 4.471676201504)], - ["C", (4.460040375488, 4.460040375488, 0.881265115392)], - ["C", (4.459048614208, 4.459048614208, 4.416542338016)], - ["C", (-0.006901660896, 1.786905622208, 1.780569760480)], - ["C", (0.002512956672, 1.784927308928, 5.351072569760)], - ["C", (0.006719393184, 5.349206041920, 1.786793836768)], - ["C", (-0.006901660896, 5.355430118208, 5.349094256480)], - ["C", (3.557699819104, 1.801090513056, 1.799927287968)], - ["C", (3.623443694336, 1.707240281568, 5.428759597120)], - ["C", (3.570268719936, 5.363769868640, 1.772230010048)], - ["C", (3.557699819104, 5.336072590720, 5.334909365632)], - ["C", (2.675027363200, 2.675027363200, 0.892428602432)], - ["C", (2.675779925760, 6.243436648480, 0.892563294432)], - ["C", (2.675027363200, 6.243571340480, 4.460972401312)], - ["C", (6.243436648480, 2.675779925760, 0.892563294432)], - ["C", (6.243571340480, 2.675027363200, 4.460972401312)], - ["C", (6.244306769504, 6.244306769504, 0.891693173408)], - ["C", (6.243436648480, 6.243436648480, 4.460219838752)], - ["C", (1.786905622208, -0.006901660896, 1.780569760480)], - ["C", (1.784927308928, 0.002512956672, 5.351072569760)], - ["C", (1.801090513056, 3.557699819104, 1.799927287968)], - ["C", (1.707240281568, 3.623443694336, 5.428759597120)], - ["C", (5.349206041920, 0.006719393184, 1.786793836768)], - ["C", (5.355430118208, -0.006901660896, 5.349094256480)], - ["C", (5.363769868640, 3.570268719936, 1.772230010048)], - ["C", (5.336072590720, 3.557699819104, 5.334909365632)], - ["C", (2.664323563008, 0.906073512160, 2.694339584224)], - ["C", (2.684866651264, 0.894090626784, 6.241909316128)], - ["C", (2.719457426496, 4.459048614208, 2.676951150304)], - ["C", (2.664323563008, 4.441660180288, 6.229926430752)], - ["C", (6.235829736704, 0.899084378176, 2.679298965664)], - ["C", (6.244294788160, 0.892434989152, 6.243564953760)], - ["C", (6.254734827520, 4.460040375488, 2.675959389024)], - ["C", (6.235829736704, 4.456700798848, 6.236915564736)], - ["C", (1.784927308928, 1.784927308928, -0.002512956672)], - ["C", (1.707240281568, 1.707240281568, 3.512556305664)], - ["C", (1.786905622208, 5.355430118208, 0.006901660896)], - ["C", (1.801090513056, 5.336072590720, 3.578300180896)], - ["C", (5.355430118208, 1.786905622208, 0.006901660896)], - ["C", (5.336072590720, 1.801090513056, 3.578300180896)], - ["C", (5.349206041920, 5.349206041920, -0.006719393184)], - ["C", (5.363769868640, 5.363769868640, 3.565731280064)], - ["N", (3.661895773440, 3.661895773440, 3.474104226560)], - ], - dimension=3, - max_memory=90000, - verbose=5, - basis='cc-pvdz', - # create a (N+2)-electron system - charge=-3, - precision=1e-12) - -gdf = df.RSDF(cell) -gdf.auxbasis = "cc-pvdz-ri" -gdf_fname = 'gdf_ints.h5' -gdf._cderi_to_save = gdf_fname -if not os.path.isfile(gdf_fname): - gdf.build() - -# =====> Part I. Restricted ppRPA <===== -# After SCF, PySCF might fail in Makov-Payne correction -# save chkfile to restart -chkfname = 'scf.chk' -if os.path.isfile(chkfname): - kmf = dft.RKS(cell).rs_density_fit() - kmf.xc = "b3lyp" - kmf.exxdiv = None - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - data = lib.chkfile.load(chkfname, 'scf') - kmf.__dict__.update(data) -else: - kmf = dft.RKS(cell).rs_density_fit() - kmf.xc = "b3lyp" - kmf.exxdiv = None - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - kmf.chkfile = chkfname - kmf.kernel() - -# direct diagonalization, N6 scaling -# ppRPA can be solved in an active space -from pyscf.pprpa.rpprpa_direct import RppRPADirect -pp = RppRPADirect(kmf, nocc_act=50, nvir_act=50, nelec="n+2") -# number of two-electron removal states to print -pp.hh_state = 50 -# solve for singlet states -pp.kernel("s") -# solve for triplet states -pp.kernel("t") -pp.analyze() - -# =====> Part II. Unrestricted ppRPA <===== -chkfname = 'uscf.chk' -if os.path.isfile(chkfname): - kmf = dft.UKS(cell).rs_density_fit() - kmf.xc = "b3lyp" - kmf.exxdiv = None - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - data = lib.chkfile.load(chkfname, 'scf') - kmf.__dict__.update(data) -else: - kmf = dft.UKS(cell).rs_density_fit() - kmf.xc = "b3lyp" - kmf.exxdiv = None - kmf.with_df = gdf - kmf.with_df._cderi = gdf_fname - kmf.chkfile = chkfname - kmf.kernel() - -# direct diagonalization, N6 scaling -# ppRPA can be solved in an active space -from pyscf.pprpa.upprpa_direct import UppRPADirect -pp = UppRPADirect(kmf, nocc_act=50, nvir_act=50, nelec="n+2") -# number of two-electron addition states to print -pp.pp_state = 50 -# solve ppRPA -pp.kernel() -pp.analyze() - diff --git a/pyscf/pprpa/__init__.py b/pyscf/pprpa/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/pyscf/pprpa/rpprpa_davidson.py b/pyscf/pprpa/rpprpa_davidson.py deleted file mode 100644 index d0a438195..000000000 --- a/pyscf/pprpa/rpprpa_davidson.py +++ /dev/null @@ -1,693 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2020 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. -# -# Author: Jiachen Li -# Author: Jincheng Yu -# - -import numpy -import scipy - -from pyscf.lib import einsum, logger, StreamObject -from pyscf.mp.mp2 import get_nocc, get_nmo - -from pyscf.pprpa.rpprpa_direct import ao2mo, inner_product, ij2index, \ - get_chemical_potential, pprpa_orthonormalize_eigenvector, pprpa_print_a_pair - - -def kernel(pprpa): - # initialize trial vector and product matrix - tri_vec = numpy.zeros(shape=[pprpa.max_vec, pprpa.full_dim], dtype=numpy.double) - tri_vec_sig = numpy.zeros(shape=[pprpa.max_vec], dtype=numpy.double) - if pprpa.channel == "pp": - ntri = min(pprpa.nroot * 4, pprpa.vv_dim) - else: - ntri = min(pprpa.nroot * 4, pprpa.oo_dim) - tri_vec[:ntri], tri_vec_sig[:ntri] = pprpa_get_trial_vector(pprpa, ntri) - mv_prod = numpy.zeros_like(tri_vec) - - iter = 0 - nprod = 0 # number of contracted vectors - while iter < pprpa.max_iter: - logger.info( - pprpa, "\nppRPA Davidson %d-th iteration, ntri= %d , nprod= %d .", - iter + 1, ntri, nprod) - mv_prod[nprod:ntri] = pprpa_contraction(pprpa, tri_vec[nprod:ntri]) - nprod = ntri - - # get ppRPA matrix and metric matrix in subspace - m_tilde = numpy.dot(tri_vec[:ntri], mv_prod[:ntri].T) - w_tilde = numpy.zeros_like(m_tilde) - for i in range(ntri): - if inner_product(tri_vec[i], tri_vec[i], pprpa.oo_dim) > 0: - w_tilde[i, i] = 1 - else: - w_tilde[i, i] = -1 - - # diagonalize subspace matrix - alphar, _, beta, _, v_tri, _, _ = scipy.linalg.lapack.dggev( - m_tilde, w_tilde, compute_vl=0) - e_tri = alphar / beta - v_tri = v_tri.T # Fortran matrix to Python order - - if pprpa.channel == "pp": - # sort eigenvalues and eigenvectors by ascending order - idx = e_tri.argsort() - e_tri = e_tri[idx] - v_tri = v_tri[idx, :] - - # re-order all states by signs, first hh then pp - sig = numpy.zeros(shape=[ntri], dtype=int) - for i in range(ntri): - if numpy.sum((v_tri[i] ** 2) * tri_vec_sig[: ntri]) > 0: - sig[i] = 1 - else: - sig[i] = -1 - - hh_index = numpy.where(sig < 0)[0] - pp_index = numpy.where(sig > 0)[0] - e_tri_hh = e_tri[hh_index] - e_tri_pp = e_tri[pp_index] - e_tri[:len(hh_index)] = e_tri_hh - e_tri[len(hh_index):] = e_tri_pp - v_tri_hh = v_tri[hh_index] - v_tri_pp = v_tri[pp_index] - v_tri[:len(hh_index)] = v_tri_hh - v_tri[len(hh_index):] = v_tri_pp - - # get only two-electron addition energy - first_state=len(hh_index) - pprpa.exci = e_tri[first_state:first_state+pprpa.nroot] - else: - # sort eigenvalues and eigenvectors by descending order - idx = e_tri.argsort()[::-1] - e_tri = e_tri[idx] - v_tri = v_tri[idx, :] - - # re-order all states by signs, first pp then hh - sig = numpy.zeros(shape=[ntri], dtype=int) - for i in range(ntri): - if numpy.sum((v_tri[i] ** 2) * tri_vec_sig[:ntri]) > 0: - sig[i] = 1 - else: - sig[i] = -1 - - hh_index = numpy.where(sig < 0)[0] - pp_index = numpy.where(sig > 0)[0] - e_tri_hh = e_tri[hh_index] - e_tri_pp = e_tri[pp_index] - e_tri[:len(pp_index)] = e_tri_pp - e_tri[len(pp_index):] = e_tri_hh - v_tri_hh = v_tri[hh_index] - v_tri_pp = v_tri[pp_index] - v_tri[:len(pp_index)] = v_tri_pp - v_tri[len(pp_index):] = v_tri_hh - - # get only two-electron removal energy - first_state=len(pp_index) - pprpa.exci = e_tri[first_state:first_state+pprpa.nroot] - - ntri_old = ntri - conv, ntri = pprpa_expand_space( - pprpa=pprpa, first_state=first_state, tri_vec=tri_vec, - tri_vec_sig=tri_vec_sig, mv_prod=mv_prod, v_tri=v_tri) - logger.info(pprpa, "add %d new trial vectors.", ntri - ntri_old) - - iter += 1 - if conv is True: - break - - assert conv is True, "ppRPA Davidson is not converged!" - logger.info( - pprpa, "\nppRPA Davidson converged in %d iterations, final subspace size = %d", - iter, nprod) - - pprpa_orthonormalize_eigenvector(pprpa.multi, pprpa.nocc_act, pprpa.exci, pprpa.xy) - - return - - -# Davidson algorithm functions -def pprpa_get_trial_vector(pprpa, ntri): - """Generate initial trial vectors in particle-particle or hole-hole channel. - The order is determined by the pair orbital energy summation. - The initial trial vectors are diagonal, and signatures are all 1 or -1. - - Args: - pprpa (ppRPA_Davidson): ppRPA_Davidson object. - ntri (int): the number of initial trial vectors. - - Returns: - tri_vec (double ndarray): initial trial vectors. - tri_vec_sig (double ndarray): signature of initial trial vectors. - """ - # for convenience, use XX as XX_act in this function - nocc, nvir, nmo = pprpa.nocc_act, pprpa.nvir_act, pprpa.nmo_act - mo_energy = pprpa.mo_energy_act - - is_singlet = 1 if pprpa.multi == "s" else 0 - - max_orb_sum = 1.0e15 - - class pair(): - def __init__(self): - self.p = -1 - self.q = -1 - self.eig_sum = max_orb_sum if pprpa.channel == "pp" else -max_orb_sum - - pairs = [] - for r in range(ntri): - t = pair() - pairs.append(t) - - if pprpa.channel == "pp": - # find particle-particle pairs with lowest orbital energy summation - for r in range(ntri): - for p in range(nocc, nmo): - for q in range(nocc, p + is_singlet): - valid = True - for rr in range(r): - if pairs[rr].p == p and pairs[rr].q == q: - valid = False - break - if (valid is True - and (mo_energy[p] + mo_energy[q]) < pairs[r].eig_sum): - pairs[r].p, pairs[r].q = p, q - pairs[r].eig_sum = mo_energy[p] + mo_energy[q] - - # sort pairs by ascending energy order - for i in range(ntri-1): - for j in range(i+1, ntri): - if pairs[i].eig_sum > pairs[j].eig_sum: - p_tmp, q_tmp, eig_sum_tmp = \ - pairs[i].p, pairs[i].q, pairs[i].eig_sum - pairs[i].p, pairs[i].q, pairs[i].eig_sum = \ - pairs[j].p, pairs[j].q, pairs[j].eig_sum - pairs[j].p, pairs[j].q, pairs[j].eig_sum = \ - p_tmp, q_tmp, eig_sum_tmp - - assert pairs[ntri-1].eig_sum < max_orb_sum, \ - "cannot find enough pairs for trial vectors" - - tri_vec = numpy.zeros(shape=[ntri, pprpa.full_dim], dtype=numpy.double) - tri_vec_sig = numpy.zeros(shape=[ntri], dtype=numpy.double) - tri_row_v, tri_col_v = numpy.tril_indices(nvir, is_singlet-1) - for r in range(ntri): - p, q = pairs[r].p, pairs[r].q - pq = ij2index(p - nocc, q - nocc, tri_row_v, tri_col_v) - tri_vec[r, pprpa.oo_dim + pq] = 1.0 - tri_vec_sig[r] = 1.0 - else: - # find hole-hole pairs with highest orbital energy summation - for r in range(ntri): - for p in range(nocc-1, -1, -1): - for q in range(nocc-1, p - is_singlet, -1): - valid = True - for rr in range(r): - if pairs[rr].p == q and pairs[rr].q == p: - valid = False - break - if (valid is True - and (mo_energy[p] + mo_energy[q]) > pairs[r].eig_sum): - pairs[r].p, pairs[r].q = q, p - pairs[r].eig_sum = mo_energy[p] + mo_energy[q] - - # sort pairs by descending energy order - for i in range(ntri-1): - for j in range(i+1, ntri): - if pairs[i].eig_sum < pairs[j].eig_sum: - p_tmp, q_tmp, eig_sum_tmp = \ - pairs[i].p, pairs[i].q, pairs[i].eig_sum - pairs[i].p, pairs[i].q, pairs[i].eig_sum = \ - pairs[j].p, pairs[j].q, pairs[j].eig_sum - pairs[j].p, pairs[j].q, pairs[j].eig_sum = \ - p_tmp, q_tmp, eig_sum_tmp - - assert pairs[ntri-1].eig_sum < max_orb_sum, \ - "cannot find enough pairs for trial vectors" - - tri_vec = numpy.zeros(shape=[ntri, pprpa.full_dim], dtype=numpy.double) - tri_vec_sig = numpy.zeros(shape=[ntri], dtype=numpy.double) - tri_row_o, tri_col_o = numpy.tril_indices(nocc, is_singlet-1) - for r in range(ntri): - p, q = pairs[r].p, pairs[r].q - pq = ij2index(p, q, tri_row_o, tri_col_o) - tri_vec[r, pq] = 1.0 - tri_vec_sig[r] = -1.0 - - return tri_vec, tri_vec_sig - - -def pprpa_contraction(pprpa, tri_vec): - """ppRPA contraction. - - Args: - pprpa (ppRPA_Davidson): ppRPA_Davidson object. - tri_vec (double ndarray): trial vector. - - Returns: - mv_prod (double ndarray): product between ppRPA matrix and trial vectors. - """ - # for convenience, use XX as XX_act in this function - nocc, nvir, nmo = pprpa.nocc_act, pprpa.nvir_act, pprpa.nmo_act - mo_energy = pprpa.mo_energy_act - Lpi = pprpa.Lpi - Lpa = pprpa.Lpa - naux = Lpi.shape[0] - - ntri = tri_vec.shape[0] - mv_prod = numpy.zeros(shape=[ntri, pprpa.full_dim], dtype=numpy.double) - - is_singlet = 1 if pprpa.multi == "s" else 0 - tri_row_o, tri_col_o = numpy.tril_indices(nocc, is_singlet-1) - tri_row_v, tri_col_v = numpy.tril_indices(nvir, is_singlet-1) - - z_oo = numpy.zeros(shape=[nocc, nocc], dtype=numpy.double) - z_vv = numpy.zeros(shape=[nvir, nvir], dtype=numpy.double) - for ivec in range(ntri): - z_oo[tri_row_o, tri_col_o] = tri_vec[ivec][:pprpa.oo_dim] - z_oo[numpy.diag_indices(nocc)] *= 1.0 / numpy.sqrt(2) - z_vv[tri_row_v, tri_col_v] = tri_vec[ivec][pprpa.oo_dim:] - z_vv[numpy.diag_indices(nvir)] *= 1.0 / numpy.sqrt(2) - - # Lpqz_{L,pr} = \sum_s Lpq_{L,ps} z_{rs} - Lpq_z = numpy.zeros(shape=[naux * nmo, nmo], dtype=numpy.double) - Lpq_z[:, :nocc] = Lpi.reshape(naux * nmo, nocc) @ z_oo.T - Lpq_z[:, nocc:] = Lpa.reshape(naux * nmo, nvir) @ z_vv.T - - # transpose and reshape for faster multiplication - Lpq_z = Lpq_z.reshape(naux, nmo, nmo).transpose(1, 0, 2) - Lpq_z = Lpq_z.reshape(nmo, naux * nmo) - - # MV_{pq} = \sum_{Lr} Lpq_{L,pr} Lpqz_{L,qr} \pm Lpq_{L,qr} Lpqz_{L,pr} - # NOTE: here assuming Lpq[L,p,q] = Lpq[L,q,p] for real orbitals - mv_prod_full = numpy.zeros(shape=[nmo, nmo], dtype=numpy.double) - mv_prod_full[:nocc, :nocc] = Lpq_z[:nocc] @ Lpi.reshape(naux * nmo, nocc) - mv_prod_full[nocc:, nocc:] = Lpq_z[nocc:] @ Lpa.reshape(naux * nmo, nvir) - if pprpa.multi == "s": - mv_prod_full += mv_prod_full.T - else: - mv_prod_full -= mv_prod_full.T - mv_prod_full = mv_prod_full.T - mv_prod_full[numpy.diag_indices(nmo)] *= 1.0 / numpy.sqrt(2) - mv_prod[ivec][: pprpa.oo_dim] =\ - mv_prod_full[: nocc, : nocc][tri_row_o, tri_col_o] - mv_prod[ivec][pprpa.oo_dim:] = \ - mv_prod_full[nocc:, nocc:][tri_row_v, tri_col_v] - - orb_sum_oo = (mo_energy[None, : nocc] + mo_energy[: nocc, None]) - 2.0 * pprpa.mu - orb_sum_oo = orb_sum_oo[tri_row_o, tri_col_o] - orb_sum_vv = (mo_energy[None, nocc:] + mo_energy[nocc:, None]) - 2.0 * pprpa.mu - orb_sum_vv = orb_sum_vv[tri_row_v, tri_col_v] - for ivec in range(ntri): - oz_oo = -orb_sum_oo * tri_vec[ivec][:pprpa.oo_dim] - mv_prod[ivec][:pprpa.oo_dim] += oz_oo - oz_vv = orb_sum_vv * tri_vec[ivec][pprpa.oo_dim:] - mv_prod[ivec][pprpa.oo_dim:] += oz_vv - - return mv_prod - - -def pprpa_expand_space( - pprpa, first_state, tri_vec, tri_vec_sig, mv_prod, v_tri): - """Expand trial vector space in Davidson algorithm. - - Args: - pprpa (ppRPA_Davidson): ppRPA_Davidson object. - first_state (int): index of first particle-particle or hole-hole state. - tri_vec (double ndarray): trial vector. - tri_vec_sig (int array): signature of trial vector. - mv_prod (double ndarray): product matrix of ppRPA matrix and trial vector. - v_tri (double ndarray): eigenvector of subspace matrix. - - Returns: - conv (bool): if Davidson algorithm is converged. - ntri (int): updated number of trial vectors. - """ - nocc_act, nvir_act = pprpa.nocc_act, pprpa.nvir_act - mo_energy_act = pprpa.mo_energy_act - nroot = pprpa.nroot - exci = pprpa.exci - max_vec = pprpa.max_vec - residue_thresh = pprpa.residue_thresh - - is_singlet = 1 if pprpa.multi == "s" else 0 - - tri_row_o, tri_col_o = numpy.tril_indices(nocc_act, is_singlet-1) - tri_row_v, tri_col_v = numpy.tril_indices(nvir_act, is_singlet-1) - - # take only nRoot vectors, starting from first pp channel - tmp = v_tri[first_state:(first_state+nroot)] - - # get the eigenvectors in the original space - ntri = v_tri.shape[0] - pprpa.xy = numpy.dot(tmp, tri_vec[:ntri]) - - # compute residue vectors - residue = numpy.dot(tmp, mv_prod[:ntri]) - for i in range(nroot): - residue[i][:pprpa.oo_dim] -= -exci[i] * pprpa.xy[i][:pprpa.oo_dim] - residue[i][pprpa.oo_dim:] += -exci[i] * pprpa.xy[i][pprpa.oo_dim:] - - # check convergence - conv_record = numpy.zeros(shape=[nroot], dtype=bool) - max_residue = 0 - for i in range(nroot): - max_residue = max(max_residue, abs(numpy.max(residue[i]))) - if len(residue[i][abs(residue[i]) > residue_thresh]) == 0: - conv_record[i] = True - else: - conv_record[i] = False - nconv = len(conv_record[conv_record is True]) - logger.info(pprpa, "max residue = %.6e", max_residue) - if nconv == nroot: - return True, ntri - - orb_sum_oo = mo_energy_act[None, :nocc_act] + mo_energy_act[:nocc_act, None] - orb_sum_oo -= 2.0 * pprpa.mu - orb_sum_oo = orb_sum_oo[tri_row_o, tri_col_o] - orb_sum_vv = mo_energy_act[None, nocc_act:] + mo_energy_act[nocc_act:, None] - orb_sum_vv -= 2.0 * pprpa.mu - orb_sum_vv = orb_sum_vv[tri_row_v, tri_col_v] - - # Schmidt orthogonalization - ntri_old = ntri - for iroot in range(nroot): - if conv_record[iroot] is True: - continue - - # convert residuals - residue[iroot][:pprpa.oo_dim] /= -(exci[iroot] - orb_sum_oo) - residue[iroot][pprpa.oo_dim:] /= (exci[iroot] - orb_sum_vv) - - for ivec in range(ntri): - # compute product between new vector and old vector - inp = -inner_product(residue[iroot], tri_vec[ivec], pprpa.oo_dim) - # eliminate parallel part - if tri_vec_sig[ivec] < 0: - inp = -inp - residue[iroot] += inp * tri_vec[ivec] - - # add a new trial vector - if len(residue[iroot][abs(residue[iroot]) > residue_thresh]) > 0: - assert ntri < max_vec, ( - "ppRPA Davidson expansion failed! ntri %d exceeds max_vec %d!" % - (ntri, max_vec)) - inp = inner_product(residue[iroot], residue[iroot], pprpa.oo_dim) - tri_vec_sig[ntri] = 1 if inp > 0 else -1 - tri_vec[ntri] = residue[iroot] / numpy.sqrt(abs(inp)) - ntri = ntri + 1 - - conv = True if ntri_old == ntri else False - return conv, ntri - - -# analysis functions -def pprpa_davidson_print_eigenvector(pprpa, multi, exci0, exci, xy): - """Print dominant components of an eigenvector. - - Args: - pprpa (RppRPADavidson): ppRPA object. - multi (string): multiplicity. - exci0 (double): lowest eigenvalue. - exci (double array): ppRPA eigenvalue. - xy (double ndarray): ppRPA eigenvector. - """ - nocc = pprpa.nocc - nocc_act, nvir_act = pprpa.nocc_act, pprpa.nvir_act - if multi == "s": - oo_dim = int((nocc_act + 1) * nocc_act / 2) - is_singlet = 1 - logger.info(pprpa, "\n print ppRPA excitations: singlet\n") - elif multi == "t": - oo_dim = int((nocc_act - 1) * nocc_act / 2) - is_singlet = 0 - logger.info(pprpa, "\n print ppRPA excitations: triplet\n") - - tri_row_o, tri_col_o = numpy.tril_indices(nocc_act, is_singlet-1) - tri_row_v, tri_col_v = numpy.tril_indices(nvir_act, is_singlet-1) - - nroot = len(exci) - au2ev = 27.211386 - if pprpa.channel == "pp": - for iroot in range(nroot): - logger.info(pprpa, "#%-d %s excitation: exci= %-12.4f eV 2e= %-12.4f eV", - iroot + 1, multi, (exci[iroot] - exci0) * au2ev, exci[iroot] * au2ev) - if nocc_act > 0: - full = numpy.zeros(shape=[nocc_act, nocc_act], dtype=numpy.double) - full[tri_row_o, tri_col_o] = xy[iroot][:oo_dim] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > pprpa.print_thresh) - for i, j in pairs: - pprpa_print_a_pair( - pprpa, is_pp=False, p=i, q=j, percentage=full[i, j]) - - full = numpy.zeros(shape=[nvir_act, nvir_act], dtype=numpy.double) - full[tri_row_v, tri_col_v] = xy[iroot][oo_dim:] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > pprpa.print_thresh) - for a, b in pairs: - pprpa_print_a_pair( - pprpa, is_pp=True, p=a+nocc, q=b+nocc, percentage=full[a, b]) - logger.info(pprpa, "") - else: - for iroot in range(nroot): - logger.info(pprpa, "#%-d %s de-excitation: exci= %-12.4f eV 2e= %-12.4f eV", - iroot + 1, multi, (exci[iroot] - exci0) * au2ev, exci[iroot] * au2ev) - full = numpy.zeros(shape=[nocc_act, nocc_act], dtype=numpy.double) - full[tri_row_o, tri_col_o] = xy[iroot][:oo_dim] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > pprpa.print_thresh) - for i, j in pairs: - pprpa_print_a_pair( - pprpa, is_pp=False, p=i, q=j, percentage=full[i, j]) - - if nvir_act > 0: - full = numpy.zeros(shape=[nvir_act, nvir_act], dtype=numpy.double) - full[tri_row_v, tri_col_v] = xy[iroot][oo_dim:] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > pprpa.print_thresh) - for a, b in pairs: - pprpa_print_a_pair( - pprpa, is_pp=True, p=a+nocc, q=b+nocc, percentage=full[a, b]) - logger.info(pprpa, "") - - return - - -def analyze_pprpa_davidson(pprpa): - """Analyze ppRPA (Davidson algorithm) excited states. - - Args: - pprpa (RppRPADavidson): ppRPA object. - """ - logger.info(pprpa, "\nanalyze ppRPA results.") - - if pprpa.exci_s is not None and pprpa.exci_t is not None: - logger.info(pprpa, "both singlet and triplet results found.") - if pprpa.channel == "pp": - exci0 = min(pprpa.exci_s[0], pprpa.exci_t[0]) - else: - exci0 = max(pprpa.exci_s[0], pprpa.exci_t[0]) - pprpa_davidson_print_eigenvector( - pprpa, multi="s", exci0=exci0, exci=pprpa.exci_s, xy=pprpa.xy_s) - pprpa_davidson_print_eigenvector( - pprpa, multi="t", exci0=exci0, exci=pprpa.exci_t, xy=pprpa.xy_t) - else: - if pprpa.exci_s is not None: - logger.info(pprpa, "only singlet results found.") - pprpa_davidson_print_eigenvector( - pprpa, multi="s", exci0=pprpa.exci_s[0], exci=pprpa.exci_s, xy=pprpa.xy_s) - else: - logger.info(pprpa, "only triplet results found.") - pprpa_davidson_print_eigenvector( - pprpa, multi="t", exci0=pprpa.exci_t[0], exci=pprpa.exci_t, xy=pprpa.xy_t) - return - - -class RppRPADavidson(StreamObject): - def __init__( - self, mf, nocc_act=None, nvir_act=None, channel="pp", nroot=5, - max_vec=200, max_iter=100, residue_thresh=1.0e-7, print_thresh=0.1, - auxbasis=None): - # necessary input - self.mol = mf.mol - self._scf = mf - self.verbose = self.mol.verbose - self.stdout = self.mol.stdout - self.max_memory = mf.max_memory - - # options - self.nocc_act = nocc_act # number of active occupied orbitals - self.nvir_act = nvir_act # number of active virtual orbitals - self.channel = channel # channel of desired states, particle-particle or hole-hole - self.nroot = nroot # number of desired roots - self.max_vec = max_vec # max size of trial vectors - self.max_iter = max_iter # max iteration - self.residue_thresh = residue_thresh # residue threshold - self.print_thresh = print_thresh # threshold to print component - self.auxbasis = auxbasis # auxiliary basis set - - # internal flags - self.multi = None # multiplicity - self.is_singlet = None # multiplicity is singlet - self.mu = None # chemical potential - self.nmo_act = None # number of active orbitals - self.mo_energy_act = None # orbital energy in active space - self.oo_dim = None # particle-particle block dimension - self.vv_dim = None # hole-hole block dimension - self.full_dim = None # full matrix dimension - - # results - self.exci = None # two-electron addition energy - self.xy = None # ppRPA eigenvector - self.exci_s = None # singlet two-electron addition energy - self.xy_s = None # singlet two-electron addition eigenvector - self.exci_t = None # triplet two-electron addition energy - self.xy_t = None # triplet two-electron addition eigenvector - - ################################################## - # don't modify the following attributes, they are not input options - self._nocc = None # number of occupied orbitals - self._nmo = None # number of molecular orbitals - self.nvir = None # number of virtual orbitals - self.mo_energy = numpy.array(self._scf.mo_energy) # orbital energy - self.Lpq = None # three-center density-fitting matrix in MO - self.Lpi = None # three-center density-fitting matrix in MO, [naux, nmo, nocc] - self.Lpa = None # three-center density-fitting matrix in MO, [naux, nmo, nvir] - self.mo_occ = self._scf.mo_occ # used in get_nocc() and get_nmo(), not in ppRPA - self.frozen = 0 # used in get_nocc() and get_nmo(), not in ppRPA - - return - - @property - def nocc(self): - return self.get_nocc() - @nocc.setter - def nocc(self, n): - self._nocc = n - - @property - def nmo(self): - return self.get_nmo() - @nmo.setter - def nmo(self, n): - self._nmo = n - - get_nocc = get_nocc - get_nmo = get_nmo - - analyze = analyze_pprpa_davidson - ao2mo = ao2mo - - def check_parameter(self): - """Initialize and check options. - """ - assert self.channel in ["pp", "hh"] - assert self.multi in ["s", "t"] - - self.nvir = self.nmo - self.nocc - - # adjust active space - if self.nocc_act is None: - self.nocc_act = self.nocc - else: - self.nocc_act = min(self.nocc_act, self.nocc) - if self.nvir_act is None: - self.nvir_act = self.nvir - else: - self.nvir_act = min(self.nvir_act, self.nvir) - - self.nmo_act = self.nocc_act + self.nvir_act - self.mo_energy_act = self.mo_energy[self.nocc-self.nocc_act:self.nocc+self.nvir_act] - - if self.multi == "s": - self.oo_dim = int((self.nocc_act + 1) * self.nocc_act / 2) - self.vv_dim = int((self.nvir_act + 1) * self.nvir_act / 2) - elif self.multi == "t": - self.oo_dim = int((self.nocc_act - 1) * self.nocc_act / 2) - self.vv_dim = int((self.nvir_act - 1) * self.nvir_act / 2) - self.full_dim = self.oo_dim + self.vv_dim - - self.max_vec = min(self.max_vec, self.full_dim) - - assert self.residue_thresh > 0 - assert 0.0 < self.print_thresh < 1.0 - - if self.mu is None: - self.mu = get_chemical_potential(self.nocc, self.mo_energy) - - return - - def dump_flags(self): - log = logger.Logger(self.stdout, self.verbose) - log.info('') - log.info('\n******** %s ********' % self.__class__) - multiplicity = "singlet" if self.multi == "s" else "triplet" - log.info('multiplicity = %s', multiplicity) - log.info('state channel = %s' % self.channel) - log.info('nmo = %d' % self.nmo) - log.info('nocc = %d nvir = %d', self.nocc, self.nvir) - log.info('nocc_act = %d nvir_act = %d', self.nocc_act, self.nvir_act) - log.info('occ-occ dimension = %d vir-vir dimension = %d', self.oo_dim, self.vv_dim) - log.info('full dimension = %d', self.full_dim) - log.info('number of roots = %d', self.nroot) - log.info('max subspace size = %d', self.max_vec) - log.info('max iteration = %d', self.max_iter) - log.info('residue threshold = %.3e', self.residue_thresh) - log.info('print threshold = %.2f%%', self.print_thresh*100) - log.info('') - return - - def check_memory(self): - """Check required memory. - """ - # intermediate in contraction; mv_prod, tri_vec, xy - mem = (3 * self.max_vec * self.full_dim) * 8 / 1.0e6 - if mem < 1000: - logger.info(self, "ppRPA needs at least %d MB memory.", mem) - else: - logger.info(self, "ppRPA needs at least %.1f GB memory.", mem/1.0e3) - return - - def kernel(self, multi): - """Run ppRPA direct diagonalization. - - Args: - multi (char): multiplicity. - """ - self.multi = multi - self.check_parameter() - self.check_memory() - - cput0 = (logger.process_clock(), logger.perf_counter()) - if self.Lpi is None or self.Lpa is None: - Lpq = self.ao2mo() - self.Lpi = numpy.ascontiguousarray(Lpq[:, :, :self.nocc_act]) - self.Lpa = numpy.ascontiguousarray(Lpq[:, :, self.nocc_act:]) - logger.timer(self, "ppRPA integral transformation: %s" % multi, *cput0) - self.dump_flags() - kernel(pprpa=self) - logger.timer(self, "ppRPA Davidson: %s" % multi, *cput0) - - if self.multi == "s": - self.exci_s = self.exci.copy() - self.xy_s = self.xy.copy() - else: - self.exci_t = self.exci.copy() - self.xy_t = self.xy.copy() - self.exci = self.xy = None - return - diff --git a/pyscf/pprpa/rpprpa_direct.py b/pyscf/pprpa/rpprpa_direct.py deleted file mode 100644 index e222245db..000000000 --- a/pyscf/pprpa/rpprpa_direct.py +++ /dev/null @@ -1,760 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2020 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. -# -# Author: Jiachen Li -# Author: Jincheng Yu -# - -import numpy -import scipy - -from pyscf import df, dft, pbc, scf -from pyscf.ao2mo._ao2mo import nr_e2 -from pyscf.lib import einsum, logger, StreamObject, current_memory -from pyscf.mp.mp2 import get_nocc, get_nmo -from pyscf.pbc.df.fft_ao2mo import _format_kpts -from pyscf.scf.hf import KohnShamDFT - - -def diagonalize_pprpa_singlet(nocc, mo_energy, Lpq, mu=None): - """Diagonalize singlet ppRPA matrix. - - Reference: - [1] https://doi.org/10.1063/1.4828728 - - Args: - nocc (int): number of occupied orbitals. - mo_energy (double array): orbital energy. - Lpq (double ndarray): three-center density-fitting matrix in MO. - mu (double, optional): chemical potential. Defaults to None. - - Returns: - exci (double array): ppRPA eigenvalue. - xy (double ndarray): ppRPA eigenvector. - ec (double): singlet correlation energy. - """ - nmo = len(mo_energy) - nvir = nmo - nocc - if mu is None: - mu = get_chemical_potential(nocc=nocc, mo_energy=mo_energy) - - oo_dim = int((nocc + 1) * nocc / 2) # number of hole-hole pairs - - # low triangular index (including diagonal) - tri_row_o, tri_col_o = numpy.tril_indices(nocc) - tri_row_v, tri_col_v = numpy.tril_indices(nvir) - - # A matrix: particle-particle block - # two-electron integral part, + - A = einsum("Pac,Pbd->abcd", Lpq[:, nocc:, nocc:], Lpq[:, nocc:, nocc:]) - A += einsum("Pad,Pbc->abcd", Lpq[:, nocc:, nocc:], Lpq[:, nocc:, nocc:]) - # scale the diagonal elements - A[numpy.diag_indices(nvir)] *= 1.0 / numpy.sqrt(2) # a=b - A = A.transpose(2, 3, 0, 1) # A_{ab,cd} to A_{cd,ab} - A[numpy.diag_indices(nvir)] *= 1.0 / numpy.sqrt(2) # c=d - A = A.transpose(2, 3, 0, 1) # A_{cd,ab} to A_{ab,cd} - # orbital energy part - A = A.reshape(nvir*nvir, nvir*nvir) - orb_sum = numpy.asarray(mo_energy[nocc:, None] + mo_energy[None, nocc:]) - orb_sum = orb_sum.reshape(-1) - 2.0 * mu - numpy.fill_diagonal(A, A.diagonal() + orb_sum) - A = A.reshape(nvir, nvir, nvir, nvir) - # take only low-triangular part - A = A[tri_row_v, tri_col_v, ...] - A = A[..., tri_row_v, tri_col_v] - trace_A = numpy.trace(A) - - # B matrix: particle-hole block - # two-electron integral part, + - B = einsum("Pai,Pbj->abij", Lpq[:, nocc:, :nocc], Lpq[:, nocc:, :nocc]) - B += einsum("Paj,Pbi->abij", Lpq[:, nocc:, :nocc], Lpq[:, nocc:, :nocc]) - # scale the diagonal elements - B[numpy.diag_indices(nvir)] *= 1.0 / numpy.sqrt(2) # a=b - B = B.transpose(2, 3, 0, 1) # B_{ab,ij} to B_{ij,ab} - B[numpy.diag_indices(nocc)] *= 1.0 / numpy.sqrt(2) # i=j - B = B.transpose(2, 3, 0, 1) # B_{ij,ab} to B_{ab,ij} - # take only low-triangular part - B = B[tri_row_v, tri_col_v, ...] - B = B[..., tri_row_o, tri_col_o] - - # C matrix: hole-hole block - # two-electron integral part, + - C = einsum("Pik,Pjl->ijkl", Lpq[:, :nocc, :nocc], Lpq[:, :nocc, :nocc]) - C += einsum("Pil,Pjk->ijkl", Lpq[:, :nocc, :nocc], Lpq[:, :nocc, :nocc]) - # scale the diagonal elements - C[numpy.diag_indices(nocc)] *= 1.0 / numpy.sqrt(2) # i=j - C = C.transpose(2, 3, 0, 1) # C_{ij,kl} to C_{kl,ij} - C[numpy.diag_indices(nocc)] *= 1.0 / numpy.sqrt(2) # k=l - C = C.transpose(2, 3, 0, 1) # C_{kl,ij} to C_{ij,kl} - # orbital energy part - C = C.reshape(nocc*nocc, nocc*nocc) - orb_sum = numpy.asarray(mo_energy[:nocc, None] + mo_energy[None, :nocc]) - orb_sum = orb_sum.reshape(-1) - 2.0 * mu - numpy.fill_diagonal(C, C.diagonal() - orb_sum) - C = C.reshape(nocc, nocc, nocc, nocc) - # take only low-triangular part - C = C[tri_row_o, tri_col_o, ...] - C = C[..., tri_row_o, tri_col_o] - - # combine A, B and C matrix as - # | C B^T | - # | B A | - M_upper = numpy.concatenate((C, B.T), axis=1) - M_lower = numpy.concatenate((B, A), axis=1) - M = numpy.concatenate((M_upper, M_lower), axis=0) - del A, B, C - # M to WM, W is the metric matrix [[-I,0],[0,I]] - M[:oo_dim][:] *= -1.0 - - # diagonalize ppRPA matrix - exci, xy = scipy.linalg.eig(M) - exci = exci.real - xy = xy.T # Fortran to Python order - - # sort eigenvalue and eigenvectors by ascending order - idx = exci.argsort() - exci = exci[idx] - xy = xy[idx, :] - - pprpa_orthonormalize_eigenvector(multi="s", nocc=nocc, exci=exci, xy=xy) - - sum_exci = numpy.sum(exci[oo_dim:]) - ec = sum_exci - trace_A - - return exci, xy, ec - - -def diagonalize_pprpa_triplet(nocc, mo_energy, Lpq, mu=None): - """Diagonalize triplet ppRPA matrix. - - Reference: - [1] https://doi.org/10.1063/1.4828728 - - Args: - nocc (int): number of occupied orbitals. - mo_energy (double array): orbital energy. - Lpq (double ndarray): three-center density-fitting matrix in MO. - mu (double, optional): chemical potential. Defaults to None. - - Returns: - exci (double array): ppRPA eigenvalue. - xy (double ndarray): ppRPA eigenvector. - ec (double): triplet correlation energy, with a factor of 3. - """ - nmo = len(mo_energy) - nvir = nmo - nocc - if mu is None: - mu = get_chemical_potential(nocc=nocc, mo_energy=mo_energy) - - oo_dim = int((nocc - 1) * nocc / 2) # number of hole-hole pairs - - # low triangular index (not including diagonal) - tri_row_o, tri_col_o = numpy.tril_indices(nocc, -1) - tri_row_v, tri_col_v = numpy.tril_indices(nvir, -1) - - # A matrix: particle-particle block - # two-electron integral part, - - A = einsum("Pac,Pbd->abcd", Lpq[:, nocc:, nocc:], Lpq[:, nocc:, nocc:]) - A -= einsum("Pad,Pbc->abcd", Lpq[:, nocc:, nocc:], Lpq[:, nocc:, nocc:]) - # orbital energy part - A = A.reshape(nvir*nvir, nvir*nvir) - orb_sum = numpy.asarray(mo_energy[nocc:, None] + mo_energy[None, nocc:]) - orb_sum = orb_sum.reshape(-1) - 2.0 * mu - numpy.fill_diagonal(A, A.diagonal() + orb_sum) - A = A.reshape(nvir, nvir, nvir, nvir) - # take only low-triangular part - A = A[tri_row_v, tri_col_v, ...] - A = A[..., tri_row_v, tri_col_v] - trace_A = numpy.trace(A) - - # B matrix: particle-hole block - # two-electron integral part, - - B = einsum("Pai,Pbj->abij", Lpq[:, nocc:, :nocc], Lpq[:, nocc:, :nocc]) - B -= einsum("Paj,Pbi->abij", Lpq[:, nocc:, :nocc], Lpq[:, nocc:, :nocc]) - # take only low-triangular part - B = B[tri_row_v, tri_col_v, ...] - B = B[..., tri_row_o, tri_col_o] - - # C matrix: hole-hole block - # two-electron integral part, - - C = einsum("Pik,Pjl->ijkl", Lpq[:, :nocc, :nocc], Lpq[:, :nocc, :nocc]) - C -= einsum("Pil,Pjk->ijkl", Lpq[:, :nocc, :nocc], Lpq[:, :nocc, :nocc]) - # orbital energy part - C = C.reshape(nocc*nocc, nocc*nocc) - orb_sum = numpy.asarray(mo_energy[:nocc, None] + mo_energy[None, :nocc]) - orb_sum = orb_sum.reshape(-1) - 2.0 * mu - numpy.fill_diagonal(C, C.diagonal() - orb_sum) - C = C.reshape(nocc, nocc, nocc, nocc) - # take only low-triangular part - C = C[tri_row_o, tri_col_o, ...] - C = C[..., tri_row_o, tri_col_o] - - # combine A, B and C matrix as - # | C B^T | - # |B A | - M_upper = numpy.concatenate((C, B.T), axis=1) - M_lower = numpy.concatenate((B, A), axis=1) - M = numpy.concatenate((M_upper, M_lower), axis=0) - del A, B, C - # M to WM, W is the metric matrix [[-I,0],[0,I]] - M[:oo_dim][:] *= -1.0 - - # diagonalize ppRPA matrix - exci, xy = scipy.linalg.eig(M) - exci = exci.real - xy = xy.T # Fortran to Python order - - # sort eigenvalue and eigenvectors by ascending order - idx = exci.argsort() - exci = exci[idx] - xy = xy[idx, :] - - pprpa_orthonormalize_eigenvector(multi="t", nocc=nocc, exci=exci, xy=xy) - - sum_exci = numpy.sum(exci[oo_dim:]) - ec = (sum_exci - trace_A) * 3.0 - - return exci, xy, ec - - -# utility function -def ij2index(r, c, row, col): - """Get index of a row and column in a square matrix in a lower triangular matrix. - - Args: - r (int): row index in s square matrix. - c (int): column index in s square matrix. - row (int array): row index array of a lower triangular matrix. - col (int array): column index array of a lower triangular matrix. - - Returns: - i (int): index in the lower triangular matrix. - """ - for i in range(len(row)): - if r == row[i] and c == col[i]: - return i - - raise ValueError("cannot find the index!") - - -def inner_product(u, v, oo_dim): - """Calculate inner product between two ppRPA eigenvectors. - product = - , where X is occ-occ block, Y is vir-vir block. - - Args: - u (double array): first vector. - v (double array): second vector - oo_dim (int): occ-occ block dimension - - Returns: - inp (double): inner product. - """ - product = numpy.sum(u[oo_dim:] * v[oo_dim:]) - product -= numpy.sum(u[:oo_dim] * v[:oo_dim]) - return product - - -def get_chemical_potential(nocc, mo_energy): - """Get chemical potential as the average between HOMO and LUMO. - In the case there is no occupied or virtual orbital, return 0. - - Args: - nocc (int): number of occupied orbitals. - mo_energy (double array/list): orbital energy. - - Returns: - mu (double): chemical potential. - """ - if nocc == 0: - mu = 0.0 - else: - mu = (mo_energy[nocc-1] + mo_energy[nocc]) * 0.5 - return mu - - -def ao2mo(pprpa): - """Get three-center density-fitting matrix in MO active space. - - Args: - pprpa: ppRPA object. - - Returns: - Lpq (double ndarray): three-center DF matrix in MO active space. - """ - mf = pprpa._scf - mo_coeff = mf.mo_coeff - nocc = pprpa.nocc - nocc_act, nvir_act, nmo_act = pprpa.nocc_act, pprpa.nvir_act, pprpa.nmo_act - - nao = mo_coeff.shape[0] - mo = numpy.asarray(mo_coeff, order='F') - ijslice = (nocc-nocc_act, nocc+nvir_act, nocc-nocc_act, nocc+nvir_act) - - if isinstance(mf, (scf.rhf.RHF, dft.rks.RKS)): - # molecule - if getattr(mf, 'with_df', None): - pprpa.with_df = mf.with_df - else: - pprpa.with_df = df.DF(mf.mol) - if pprpa.auxbasis is not None: - pprpa.with_df.auxbasis = pprpa.auxbasis - else: - pprpa.with_df.auxbasis = df.make_auxbasis( - mf.mol, mp2fit=True) - pprpa._keys.update(['with_df']) - - naux = pprpa.with_df.get_naoaux() - mem_incore = (2*nao**2*naux) * 8/1e6 - mem_now = current_memory()[0] - - Lpq = None - if (mem_incore + mem_now < pprpa.max_memory) or pprpa.mol.incore_anyway: - Lpq = nr_e2(pprpa.with_df._cderi, mo, ijslice, aosym='s2', out=Lpq) - return Lpq.reshape(naux, nmo_act, nmo_act) - else: - logger.warn(pprpa, 'Memory may not be enough!') - raise NotImplementedError - elif isinstance(mf, (pbc.scf.rhf.RHF, pbc.dft.rks.RKS)): - # supercell - if getattr(mf, 'with_df', None): - pprpa.with_df = mf.with_df - else: - pprpa.with_df = df.DF(mf.mol) - if pprpa.auxbasis is not None: - pprpa.with_df.auxbasis = pprpa.auxbasis - else: - pprpa.with_df.auxbasis = pbc.df.make_auxbasis( - mf.mol, mp2fit=True) - pprpa._keys.update(['with_df']) - - naux = pprpa.with_df.get_naoaux() - mem_incore = (nao**2*naux) * 8/1e6 - mem_now = current_memory()[0] - max_memory = max(4000, pprpa.max_memory - mem_now - mem_incore) - - kptijkl = _format_kpts(mf.with_df.kpts) - Lpq = [] - for LpqR, _, _ in mf.with_df.sr_loop(kptijkl[:2], - max_memory=0.8*max_memory, - compact=False): - LpqR = LpqR.reshape(-1, nao, nao) - tmp = None - tmp = nr_e2(LpqR, mo, ijslice, aosym='s1', mosym='s1', out=tmp) - Lpq.append(tmp) - Lpq = numpy.vstack(Lpq).reshape(naux, nmo_act, nmo_act) - return Lpq - - -def pprpa_orthonormalize_eigenvector(multi, nocc, exci, xy): - """Orthonormalize ppRPA eigenvector. - The eigenvector is normalized as Y^2 - X^2 = 1. - This function will rewrite input exci and xy, after calling this function, - exci and xy will be re-ordered as [hole-hole, particle-particle]. - - Args: - multi (string): multiplicity. - nocc (int): number of occupied orbitals. - exci (double array): ppRPA eigenvalue. - xy (double ndarray): ppRPA eigenvector. - """ - nroot = xy.shape[0] - - if multi == "s": - oo_dim = int((nocc + 1) * nocc / 2) - elif multi == "t": - oo_dim = int((nocc - 1) * nocc / 2) - - # determine the vector is pp or hh - sig = numpy.zeros(shape=[nroot], dtype=numpy.double) - for i in range(nroot): - sig[i] = 1 if inner_product(xy[i], xy[i], oo_dim) > 0 else -1 - - # eliminate parallel component - for i in range(nroot): - for j in range(i): - if abs(exci[i] - exci[j]) < 1.0e-7: - inp = inner_product(xy[i], xy[j], oo_dim) - xy[i] -= sig[j] * xy[j] * inp - - # normalize - for i in range(nroot): - inp = inner_product(xy[i], xy[i], oo_dim) - inp = numpy.sqrt(abs(inp)) - xy[i] /= inp - - # re-order all states by signs, first hh then pp - hh_index = numpy.where(sig < 0)[0] - pp_index = numpy.where(sig > 0)[0] - exci_hh = exci[hh_index] - exci_pp = exci[pp_index] - exci[:len(hh_index)] = exci_hh - exci[len(hh_index):] = exci_pp - xy_hh = xy[hh_index] - xy_pp = xy[pp_index] - xy[:len(hh_index)] = xy_hh - xy[len(hh_index):] = xy_pp - - return - - -# analysis functions -def pprpa_print_direct_eigenvector(pprpa, multi, exci0, exci, xy): - """Print dominant components of an eigenvector. - - Args: - pprpa (RppRPADirect): ppRPA object. - multi (string): multiplicity. - exci0 (double): lowest eigenvalue. - exci (double array): ppRPA eigenvalue. - xy (double ndarray): ppRPA eigenvector. - """ - nocc = pprpa.nocc - nocc_act, nvir_act = pprpa.nocc_act, pprpa.nvir_act - nocc_fro = nocc - nocc_act - if multi == "s": - oo_dim = int((nocc_act + 1) * nocc_act / 2) - vv_dim = int((nvir_act + 1) * nvir_act / 2) - is_singlet = 1 - logger.info(pprpa, "\n print ppRPA excitations: singlet\n") - elif multi == "t": - oo_dim = int((nocc_act - 1) * nocc_act / 2) - vv_dim = int((nvir_act - 1) * nvir_act / 2) - is_singlet = 0 - logger.info(pprpa, "\n print ppRPA excitations: triplet\n") - - tri_row_o, tri_col_o = numpy.tril_indices(nocc_act, is_singlet-1) - tri_row_v, tri_col_v = numpy.tril_indices(nvir_act, is_singlet-1) - - au2ev = 27.211386 - - for istate in range(min(pprpa.hh_state, oo_dim)): - logger.info(pprpa, "#%-d %s de-excitation: exci= %-12.4f eV 2e= %-12.4f eV", - istate + 1, multi, (exci[oo_dim-istate-1] - exci0) * au2ev, - exci[oo_dim-istate-1] * au2ev) - full = numpy.zeros(shape=[nocc_act, nocc_act], dtype=numpy.double) - full[tri_row_o, tri_col_o] = xy[oo_dim-istate-1][:oo_dim] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > pprpa.print_thresh) - for i, j in pairs: - pprpa_print_a_pair( - pprpa, is_pp=False, p=i+nocc_fro, q=j+nocc_fro, - percentage=full[i, j]) - - full = numpy.zeros(shape=[nvir_act, nvir_act], dtype=numpy.double) - full[tri_row_v, tri_col_v] = xy[oo_dim-istate-1][oo_dim:] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > pprpa.print_thresh) - for a, b in pairs: - pprpa_print_a_pair( - pprpa, s_pp=True, p=a+nocc_fro+nocc_act, q=b+nocc_fro+nocc_act, - percentage=full[a, b]) - - logger.info(pprpa, "") - - for istate in range(min(pprpa.pp_state, vv_dim)): - logger.info(pprpa, "#%-d %s excitation: exci= %-12.4f eV 2e= %-12.4f eV", - istate + 1, multi, (exci[oo_dim+istate] - exci0) * au2ev, - exci[oo_dim+istate] * au2ev) - full = numpy.zeros(shape=[nocc_act, nocc_act], dtype=numpy.double) - full[tri_row_o, tri_col_o] = xy[oo_dim+istate][:oo_dim] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > pprpa.print_thresh) - for i, j in pairs: - pprpa_print_a_pair( - pprpa, is_pp=False, p=i+nocc_fro, q=j+nocc_fro, - percentage=full[i, j]) - - full = numpy.zeros(shape=[nvir_act, nvir_act], dtype=numpy.double) - full[tri_row_v, tri_col_v] = xy[oo_dim+istate][oo_dim:] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > pprpa.print_thresh) - for a, b in pairs: - pprpa_print_a_pair( - pprpa, is_pp=True, p=a+nocc_fro+nocc_act, q=b+nocc_fro+nocc_act, - percentage=full[a, b]) - - logger.info(pprpa, "") - - return - - -def analyze_pprpa_direct(pprpa): - """Analyze ppRPA (direct diagonalization) excited states. - - Args: - pprpa (RppRPADirect): ppRPA object. - """ - oo_dim_s = int((pprpa.nocc_act + 1) * pprpa.nocc_act / 2) - oo_dim_t = int((pprpa.nocc_act - 1) * pprpa.nocc_act / 2) - - logger.info(pprpa, "\nanalyze ppRPA results.") - if pprpa.exci_s is not None and pprpa.exci_t is not None: - logger.info(pprpa, "both singlet and triplet results found.") - if pprpa.nelec == "n-2": - exci0 = min(pprpa.exci_s[oo_dim_s], pprpa.exci_t[oo_dim_t]) - else: - exci0 = max(pprpa.exci_s[oo_dim_s-1], pprpa.exci_t[oo_dim_t-1]) - pprpa_print_direct_eigenvector( - pprpa, multi="s", exci0=exci0, exci=pprpa.exci_s, xy=pprpa.xy_s) - pprpa_print_direct_eigenvector( - pprpa, multi="t", exci0=exci0, exci=pprpa.exci_t, xy=pprpa.xy_t) - else: - if pprpa.exci_s is not None: - logger.info(pprpa, "only singlet results found.") - exci0 = pprpa.exci_s[oo_dim_s if pprpa.nelec == "n-2" else oo_dim_s-1] - pprpa_print_direct_eigenvector( - pprpa, multi="s", exci0=exci0, exci=pprpa.exci_s, xy=pprpa.xy_s) - else: - logger.info(pprpa, "only triplet results found.") - exci0 = pprpa.exci_s[oo_dim_t if pprpa.nelec == "n-2" else oo_dim_t-1] - pprpa_print_direct_eigenvector( - pprpa, multi="t", exci0=exci0, exci=pprpa.exci_t, xy=pprpa.xy_t) - return - - -def pprpa_print_a_pair(pprpa, is_pp, p, q, percentage): - """Print the percentage of a pair in the eigenvector. - - Args: - pprpa: ppRPA object. - is_pp (bool): the eigenvector is in particle-particle channel. - p (int): MO index of the first orbital. - q (int): MO index of the second orbital. - percentage (double): the percentage of this pair. - """ - if is_pp: - logger.info(pprpa, " particle-particle pair: %5d %5d %5.2f%%", - p + 1, q + 1, percentage * 100) - else: - logger.info(pprpa, " hole-hole pair: %5d %5d %5.2f%%", - p + 1, q + 1, percentage * 100) - return - - -class RppRPADirect(StreamObject): - def __init__( - self, mf, nocc_act=None, nvir_act=None, hh_state=5, pp_state=5, - nelec="n-2", print_thresh=0.1, auxbasis=None): - self.mol = mf.mol - self._scf = mf - self.verbose = self.mol.verbose - self.stdout = self.mol.stdout - self.max_memory = mf.max_memory - - # options - self.nocc_act = nocc_act # number of active occupied orbitals - self.nvir_act = nvir_act # number of active virtual orbitals - self.hh_state = hh_state # number of hole-hole states to print - self.pp_state = pp_state # number of particle-particle states to print - self.nelec = nelec # "n-2" or "n+2" for system is an N-2 or N+2 system - self.print_thresh = print_thresh # threshold to print component - self.auxbasis = auxbasis # auxiliary basis set - - # internal flags - self.multi = None # multiplicity - self.mu = None # chemical potential - self.nmo_act = None # number of active orbitals - self.mo_energy_act = None # orbital energy in active space - - # results - self.ec = None # correlation energy - self.ec_s = None # singlet correlation energy - self.ec_t = None # triplet correlation energy - self.exci = None # two-electron addition energy - self.xy = None # ppRPA eigenvector - self.exci_s = None # singlet two-electron addition energy - self.xy_s = None # singlet two-electron addition eigenvector - self.exci_t = None # triplet two-electron addition energy - self.xy_t = None # triplet two-electron addition eigenvector - - ################################################## - # don't modify the following attributes, they are not input options - self._nocc = None # number of occupied orbitals - self._nmo = None # number of molecular orbitals - self.nvir = None # number of virtual orbitals - self.mo_energy = numpy.array(self._scf.mo_energy) # orbital energy - self.Lpq = None # three-center density-fitting matrix in MO - self.mo_occ = self._scf.mo_occ # used in get_nocc() and get_nmo(), not in ppRPA - self.frozen = 0 # used in get_nocc() and get_nmo(), not in ppRPA - - return - - @property - def nocc(self): - return self.get_nocc() - @nocc.setter - def nocc(self, n): - self._nocc = n - - @property - def nmo(self): - return self.get_nmo() - @nmo.setter - def nmo(self, n): - self._nmo = n - - get_nocc = get_nocc - get_nmo = get_nmo - - ao2mo = ao2mo - analyze = analyze_pprpa_direct - - def check_parameter(self): - """Initialize and check options. - """ - assert 0.0 < self.print_thresh < 1.0 - assert self.nelec in ["n-2", "n+2"] - - self.nvir = self.nmo - self.nocc - - # adjust active space - if self.nocc_act is None: - self.nocc_act = self.nocc - else: - self.nocc_act = min(self.nocc_act, self.nocc) - if self.nvir_act is None: - self.nvir_act = self.nvir - else: - self.nvir_act = min(self.nvir_act, self.nvir) - - self.nmo_act = self.nocc_act + self.nvir_act - self.mo_energy_act = self.mo_energy[self.nocc-self.nocc_act:self.nocc+self.nvir_act] - - if self.mu is None: - self.mu = get_chemical_potential(self.nocc, self.mo_energy) - - return - - def dump_flags(self): - log = logger.Logger(self.stdout, self.verbose) - log.info('') - log.info('\n******** %s ********', self.__class__) - if self.multi == "s": - oo_dim = int((self.nocc_act + 1) * self.nocc_act / 2) - vv_dim = int((self.nvir_act + 1) * self.nvir_act / 2) - elif self.multi == "t": - oo_dim = int((self.nocc_act - 1) * self.nocc_act / 2) - vv_dim = int((self.nvir_act - 1) * self.nvir_act / 2) - full_dim = oo_dim + vv_dim - multiplicity = "singlet" if self.multi == "s" else "triplet" - log.info('multiplicity = %s', multiplicity) - log.info('nmo = %d', self.nmo) - log.info('nocc = %d nvir = %d', self.nocc, self.nvir) - log.info('nocc_act = %d nvir_act = %d', self.nocc_act, self.nvir_act) - log.info('occ-occ dimension = %d vir-vir dimension = %d', oo_dim, vv_dim) - log.info('full dimension = %d', full_dim) - log.info('interested hh state = %d', self.hh_state) - log.info('interested pp state = %d', self.pp_state) - log.info('ground state = %s', self.nelec) - log.info('print threshold = %.2f%%', self.print_thresh*100) - log.info('') - return - - def check_memory(self): - """Check required memory. - In direct diagonalization, dominant memory cost is saving A and full - ppRPA matrix. - """ - log = logger.Logger(self.stdout, self.verbose) - if self.multi == "s": - oo_dim = int((self.nocc + 1) * self.nocc / 2) - vv_dim = int((self.nvir + 1) * self.nvir / 2) - elif self.multi == "t": - oo_dim = int((self.nocc - 1) * self.nocc / 2) - vv_dim = int((self.nvir - 1) * self.nvir / 2) - full_dim = oo_dim + vv_dim - - # ppRPA matrix: A block and full matrix, eigenvector - mem = (3 * full_dim * full_dim) * 8 / 1.0e6 - if mem < 1000: - log.info("ppRPA needs at least %d MB memory.", mem) - else: - log.info("ppRPA needs at least %.1f GB memory.", mem / 1.0e3) - return - - def kernel(self, multi): - """Run ppRPA direct diagonalization. - - Args: - multi (char): multiplicity. - """ - self.multi = multi - self.check_parameter() - self.check_memory() - - cput0 = (logger.process_clock(), logger.perf_counter()) - self.dump_flags() - if self.Lpq is None: - self.Lpq = self.ao2mo() - if self.multi == "s": - self.exci_s, self.xy_s, self.ec_s = diagonalize_pprpa_singlet( - nocc=self.nocc_act, mo_energy=self.mo_energy_act, Lpq=self.Lpq, - mu=self.mu) - elif multi == "t": - self.exci_t, self.xy_t, self.ec_t = diagonalize_pprpa_triplet( - nocc=self.nocc_act, mo_energy=self.mo_energy_act, Lpq=self.Lpq, - mu=self.mu) - logger.timer(self, "ppRPA direct: %s" % multi, *cput0) - return - - def get_correlation(self): - """Get ppRPA correlation energy. - Triplet contribution is multiplied by a factor of 3. - - Reference: - [1] https://doi.org/10.1063/1.4828728 - - Returns: - ec (double): ppRPA correlation energy. - """ - self.check_parameter() - - if self.Lpq is None: - self.Lpq = self.ao2mo() - if self.ec_s is None: - cput0 = (logger.process_clock(), logger.perf_counter()) - self.exci_s, self.xy_s, self.ec_s = diagonalize_pprpa_singlet( - nocc=self.nocc_act, mo_energy=self.mo_energy_act, Lpq=self.Lpq, - mu=self.mu) - logger.timer(self, "ppRPA correlation energy: singlet", *cput0) - - if self.ec_t is None: - cput0 = (logger.process_clock(), logger.perf_counter()) - self.exci_t, self.xy_t, self.ec_t = diagonalize_pprpa_triplet( - nocc=self.nocc_act, mo_energy=self.mo_energy_act, Lpq=self.Lpq, - mu=self.mu) - logger.timer(self, "ppRPA correlation energy: triplet", *cput0) - - self.ec = self.ec_s + self.ec_t - - return self.ec - - def energy_tot(self): - """Get ppRPA total energy. - Total energy = Hartree-Fock energy + ppRPA correlation energy. - - Returns: - e_tot (double): ppRPA total energy. - """ - mf = self._scf - assert mf.converged - hf_obj = mf if not isinstance(mf, KohnShamDFT) else mf.to_hf() - - dm = hf_obj.make_rdm1() - e_hf = hf_obj.energy_nuc() + hf_obj.energy_elec(dm=dm)[0] - ec = self.get_correlation() - e_tot = e_hf + ec - - return e_tot, e_hf, ec diff --git a/pyscf/pprpa/tests/test_rpprpa.py b/pyscf/pprpa/tests/test_rpprpa.py deleted file mode 100644 index 899c805bc..000000000 --- a/pyscf/pprpa/tests/test_rpprpa.py +++ /dev/null @@ -1,121 +0,0 @@ -#!/usr/bin/env python -# Copyright 2024 The PySCF Developers. All Rights Reserved. -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -import unittest -from pyscf import gto, dft, lib -from pyscf.pprpa.rpprpa_direct import RppRPADirect -from pyscf.pprpa.rpprpa_davidson import RppRPADavidson -from pyscf.pprpa.upprpa_direct import UppRPADirect - -def setUpModule(): - global mol, rmf, umf - mol = gto.Mole() - mol.verbose = 5 - mol.output = '/dev/null' - mol.atom = [ - ["O", (0.0, 0.0, 0.0)], - ["H", (0.0, -0.7571, 0.5861)], - ["H", (0.0, 0.7571, 0.5861)]] - mol.basis = 'def2-svp' - # create a (N-2)-electron system for charged-neutral H2O - mol.charge = 2 - mol.build() - rmf = mol.RKS(xc = "b3lyp").run() - umf = mol.UKS(xc = "b3lyp").run() - -def tearDownModule(): - global mol, rmf, umf - mol.stdout.close() - del mol, rmf, umf - -class KnownValues(unittest.TestCase): - def test_rpprpa_correlation_energy(self): - pp = RppRPADirect(rmf, nocc_act=None, nvir_act=10) - ec = pp.get_correlation() - etot, ehf, ec = pp.energy_tot() - self.assertAlmostEqual(ec, -0.0450238550202-0.0159221841934, 8) - - pp = RppRPADirect(rmf) - ec = pp.get_correlation() - etot, ehf, ec = pp.energy_tot() - self.assertAlmostEqual(ec, -0.11242089840288827) - - def test_upprpa_correlation_energy(self): - pp = UppRPADirect(umf) - ec = pp.get_correlation() - etot, ehf, ec = pp.energy_tot() - self.assertAlmostEqual(ec, -0.11242089840288827) - - def test_rpprpa_direct(self): - pp = RppRPADirect(rmf, nocc_act=None, nvir_act=10) - pp.pp_state = 10 - pp.kernel("s") - pp.kernel("t") - self.assertAlmostEqual(pp.ec_s, -0.0450238550202, 8) - self.assertAlmostEqual(pp.ec_t, -0.0159221841934, 8) - - ref = [0.92727944, 1.18456136, 1.25715794, 1.66117646, 1.72616615] - self.assertAlmostEqual(abs(pp.exci_s[10:15] - ref).max(), 0, 7) - self.assertAlmostEqual(lib.fp(pp.exci_s), -21.19709249893, 8) - - ref = [1.16289368, 1.24603019, 1.64840497, 1.69022118, 1.74312208] - self.assertAlmostEqual(abs(pp.exci_t[6:11] - ref).max(), 0, 7) - self.assertAlmostEqual(lib.fp(pp.exci_t), -19.45490052780, 8) - pp.analyze() - - def test_upprpa_direct(self): - pp = UppRPADirect(umf, nocc_act=None, nvir_act=10) - pp.pp_state = 10 - pp.kernel() - self.assertAlmostEqual(pp.ec[0], -0.0159221841934/3, 8) - self.assertAlmostEqual(pp.ec[1], -0.0159221841934/3, 8) - self.assertAlmostEqual(pp.ec[2], -0.0450238550202-0.0159221841934/3, 7) - self.assertAlmostEqual(lib.fp(pp.exci[0]), -19.45490020712, 8) - self.assertAlmostEqual(lib.fp(pp.exci[1]), -19.45490025537, 8) - self.assertAlmostEqual(lib.fp(pp.exci[2]), -35.07654958367, 8) - pp.analyze() - - def test_rpprpa_davidson(self): - pp = RppRPADavidson(rmf, nocc_act=None, nvir_act=10, nroot=5) - ref = [0.92727944, 1.18456136, 1.25715794, 1.66117646, 1.72616615] - pp.kernel("s") - self.assertAlmostEqual(abs(pp.exci_s - ref).max(), 0, 7) - ref = [1.16289368, 1.24603019, 1.64840497, 1.69022118, 1.74312208] - pp.kernel("t") - self.assertAlmostEqual(abs(pp.exci_t - ref).max(), 0, 7) - pp.analyze() - - def test_hhrpa(self): - mol = gto.Mole() - mol.verbose = 5 - mol.atom = [ - ["O", (0.0, 0.0, 0.0)], - ["H", (0.0, -0.7571, 0.5861)], - ["H", (0.0, 0.7571, 0.5861)]] - mol.basis = 'def2-svp' - mol.charge = -2 - mol.build() - mf = mol.RKS(xc = "b3lyp").run() - pp = RppRPADirect(mf, nvir_act=10, nelec="n+2") - pp.hh_state = 10 - pp.kernel("s") - pp.kernel("t") - self.assertAlmostEqual(pp.ec_s, -0.0837534201778, 8) - self.assertAlmostEqual(pp.ec_t, -0.0633019675263, 8) - pp.analyze() - -if __name__ == "__main__": - print('Full Tests for ppRPA') - unittest.main() diff --git a/pyscf/pprpa/upprpa_direct.py b/pyscf/pprpa/upprpa_direct.py deleted file mode 100644 index fdcb7d8d6..000000000 --- a/pyscf/pprpa/upprpa_direct.py +++ /dev/null @@ -1,751 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2020 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. -# -# Author: Jiachen Li -# Author: Jincheng Yu -# - -import numpy -import scipy -from pyscf import df, dft, pbc, scf -from pyscf.lib import einsum, logger, StreamObject, current_memory -from pyscf.mp.ump2 import get_nocc, get_nmo -from pyscf.pbc.df.fft_ao2mo import _format_kpts -from pyscf.ao2mo._ao2mo import nr_e2 -from pyscf.scf.hf import KohnShamDFT -from pyscf.pprpa.rpprpa_direct import inner_product, ij2index, \ - pprpa_orthonormalize_eigenvector, \ - pprpa_print_a_pair, diagonalize_pprpa_triplet - - -def upprpa_orthonormalize_eigenvector(subspace, nocc, exci, xy): - """Orthonormalize U-ppRPA eigenvectors. - The eigenvector is normalized as Y^2 - X^2 = 1. - This function will rewrite input exci and xy, after calling this function, - exci and xy will be re-ordered as [hole-hole, particle-particle]. - - Args: - subspace (str): subspace, 'aaaa', 'bbbb', or 'abab'. - nocc (int/tuple of int): number of occupied orbitals. - exci (double array): ppRPA eigenvalue. - xy (double ndarray): ppRPA eigenvector. - """ - nroot = xy.shape[0] - - if subspace == 'abab': - oo_dim = int(nocc[0] * nocc[1]) - - # determine the vector is pp or hh - sig = numpy.zeros(shape=[nroot], dtype=numpy.double) - for i in range(nroot): - sig[i] = 1 if inner_product(xy[i], xy[i], oo_dim) > 0 else -1 - - # eliminate parallel component - for i in range(nroot): - for j in range(i): - if abs(exci[i] - exci[j]) < 1.0e-7: - inp = inner_product(xy[i], xy[j], oo_dim) - xy[i] -= sig[j] * xy[j] * inp - - # normalize - for i in range(nroot): - inp = inner_product(xy[i], xy[i], oo_dim) - inp = numpy.sqrt(abs(inp)) - xy[i] /= inp - - # re-order all states by signs, first hh then pp - hh_index = numpy.where(sig < 0)[0] - pp_index = numpy.where(sig > 0)[0] - exci_hh = exci[hh_index] - exci_pp = exci[pp_index] - exci[:len(hh_index)] = exci_hh - exci[len(hh_index):] = exci_pp - xy_hh = xy[hh_index] - xy_pp = xy[pp_index] - xy[:len(hh_index)] = xy_hh - xy[len(hh_index):] = xy_pp - - # change |X -Y> to |X Y> - xy[:][:oo_dim] *= -1 - - else: - pprpa_orthonormalize_eigenvector('t', nocc, exci, xy) - - return - - -def get_chemical_potential(nocc, mo_energy): - """Get chemical potential as the average between HOMO and LUMO. - In the case there is no occupied or virtual orbital, return0. - - Args: - nocc (list/tuple of int): number of occupied orbitals, [alpha, beta]. - mo_energy (double ndarrya/list): orbital energy. - - Returns: - mu (double): chemical potential. - """ - nmo = (len(mo_energy[0]), len(mo_energy[1])) - if (nocc[0] == nocc[1] == 0) or (nocc[0] == nmo[0] and nocc[1] == nmo[1]): - mu = 0.0 - else: - assert nocc[0] >= nocc[1] - if nocc[1] == 0: - homo = mo_energy[0][nocc[0]-1] - else: - homo = max(mo_energy[0][nocc[0]-1], mo_energy[1][nocc[1]-1]) - lumo = min(mo_energy[0][nocc[0]], mo_energy[1][nocc[1]]) - mu = (homo + lumo) * 0.5 - - return mu - - -def ao2mo(pprpa): - """Get three-center density-fitting matrix in MO active space. - - Args: - pprpa: ppRPA object. - - Returns: - Lpq (double ndarray): three-center DF matrices in MO active space. - """ - mf = pprpa._scf - mo_coeff = numpy.asarray(mf.mo_coeff) - nocc = pprpa.nocc - nocc_act, nvir_act, nmo_act = pprpa.nocc_act, pprpa.nvir_act, pprpa.nmo_act - - nao = mo_coeff[0].shape[0] - ijslice = [ - ( - nocc[0]-nocc_act[0], nocc[0]+nvir_act[0], - nocc[0]-nocc_act[0], nocc[0]+nvir_act[0]), - ( - nocc[1]-nocc_act[1], nocc[1]+nvir_act[1], - nocc[1]-nocc_act[1], nocc[1]+nvir_act[1]), - ] - - if isinstance(mf, (scf.uhf.UHF, dft.uks.UKS)): - if getattr(mf, 'with_df', None): - pprpa.with_df = mf.with_df - else: - pprpa.with_df = df.DF(mf.mol) - if pprpa.auxbasis is not None: - pprpa.with_df.auxbasis = pprpa.auxbasis - else: - pprpa.with_df.auxbasis = df.make_auxbasis( - mf.mol, mp2fit=True) - pprpa._keys.update(['with_df']) - - naux = pprpa.with_df.get_naoaux() - mem_incore = (2*nao**2*naux) * 8/1e6 - mem_now = current_memory()[0] - - Lpq_a = None - Lpq_b = None - if (mem_incore + mem_now < pprpa.max_memory) or pprpa.mol.incore_anyway: - Lpq_a = nr_e2( - pprpa.with_df._cderi, mo_coeff[0], ijslice[0], - aosym='s2', out=Lpq_a) - Lpq_b = nr_e2( - pprpa.with_df._cderi, mo_coeff[1], ijslice[1], - aosym='s2', out=Lpq_b) - Lpq_a = Lpq_a.reshape(naux, nmo_act[0], nmo_act[0]) - Lpq_b = Lpq_b.reshape(naux, nmo_act[1], nmo_act[1]) - return numpy.asarray([Lpq_a, Lpq_b]) - else: - logger.warn(pprpa, 'Memory may not be enough!') - raise NotImplementedError - elif isinstance(mf, (pbc.scf.uhf.UHF, pbc.dft.uks.UKS)): - if getattr(mf, 'with_df', None): - pprpa.with_df = mf.with_df - else: - pprpa.with_df = df.DF(mf.mol) - if pprpa.auxbasis is not None: - pprpa.with_df.auxbasis = pprpa.auxbasis - else: - pprpa.with_df.auxbasis = pbc.df.make_auxbasis( - mf.mol, mp2fit=True) - pprpa._keys.update(['with_df']) - - naux = pprpa.with_df.get_naoaux() - mem_incore = (nao**2*naux) * 8/1e6 - mem_now = current_memory()[0] - max_memory = max(4000, pprpa.max_memory - mem_now - mem_incore) - - kptijkl = _format_kpts(mf.with_df.kpts) - mo_a = numpy.asarray(mo_coeff[0], order='F') - mo_b = numpy.asarray(mo_coeff[1], order='F') - eri_3d = [] - Lpq_a, Lpq_b = [], [] - - for LpqR, _, _ in mf.with_df.sr_loop(kptijkl[:2], - max_memory=0.3*max_memory, - compact=False): - tmp_a, tmp_b = None, None - tmp_a = nr_e2( - LpqR.reshape(-1,nao,nao), mo_a, ijslice[0], aosym='s1', - mosym='s1', out=tmp_a) - tmp_b = nr_e2( - LpqR.reshape(-1,nao,nao), mo_b, ijslice[1], aosym='s1', - mosym='s1', out=tmp_b) - Lpq_a.append(tmp_a) - Lpq_b.append(tmp_b) - Lpq_a = numpy.vstack(Lpq_a).reshape(-1, nmo_act[0], nmo_act[0]) - Lpq_b = numpy.vstack(Lpq_b).reshape(-1, nmo_act[1], nmo_act[1]) - eri_3d = numpy.asarray([Lpq_a, Lpq_b]) - - return eri_3d - - -def diagonalize_pprpa_subspace_same_spin(nocc, mo_energy, Lpq, mu=None): - """Diagonalize UppRPA matrix in subspace (alpha alpha, alpha alpha) - or (beta beta, beta beta). - - See function `pprpa_direct.diagonalize_pprpa_triplet`. - - """ - exci, xy, ec = diagonalize_pprpa_triplet(nocc, mo_energy, Lpq, mu=mu) - - return exci, xy, ec/3.0 - - -def diagonalize_pprpa_subspace_diff_spin(nocc, mo_energy, Lpq, mu=None): - """Diagonalize UppRPA matrix in subspace (alpha beta, alpha beta). - - Reference: - [1] https://doi.org/10.1063/1.4828728 (equation 14) - - Args: - nocc(tuple of int): number of occupied orbitals, (nalpha, nbeta). - mo_energy (list of double array): orbital energies. - Lpq (list of double ndarray): three-center RI matrices in MO space. - - Kwarg: - mu (double): chemical potential. - - Returns: - exci (double array): ppRPA eigenvalue. - xy (double ndarray): ppRPA eigenvector. - ec (double): correlation energy from one subspace. - """ - nmo = (len(mo_energy[0]), len(mo_energy[1])) - nvir = (nmo[0]-nocc[0], nmo[1]-nocc[1]) - if mu is None: - mu = get_chemical_potential(nocc, mo_energy) - - # ===========================> A matrix <============================ - # - A = einsum( - 'Pac,Pbd->abcd', Lpq[0][:, nocc[0]:, nocc[0]:], - Lpq[1][:, nocc[1]:, nocc[1]:], optimize=True) - # delta_ac delta_bd (e_a + e_b - 2 * mu) - A = A.reshape(nvir[0]*nvir[1], nvir[0]*nvir[1]) - orb_sum = numpy.asarray( - mo_energy[0][nocc[0]:, None] + mo_energy[1][None, nocc[1]:] - ).reshape(-1) - orb_sum -= 2.0 * mu - numpy.fill_diagonal(A, A.diagonal() + orb_sum) - trace_A = numpy.trace(A) - - # ===========================> B matrix <============================ - # - B = einsum( - 'Pai,Pbj->abij', Lpq[0][:, nocc[0]:, :nocc[0]], - Lpq[1][:, nocc[1]:, :nocc[1]], optimize=True) - B = B.reshape(nvir[0]*nvir[1], nocc[0]*nocc[1]) - - # ===========================> C matrix <============================ - # - C = einsum( - 'Pik,Pjl->ijkl', Lpq[0][:, :nocc[0], :nocc[0]], - Lpq[1][:, :nocc[1], :nocc[1]], optimize=True) - # - delta_ik delta_jl (e_i + e_j - 2 * mu) - C = C.reshape(nocc[0]*nocc[1], nocc[0]*nocc[1]) - orb_sum = numpy.asarray( - mo_energy[0][:nocc[0], None] + mo_energy[1][None, :nocc[1]] - ).reshape(-1) - orb_sum -= 2.0 * mu - numpy.fill_diagonal(C, C.diagonal() - orb_sum) - - # ==================> whole matrix in the subspace<================== - # C B^T - # B A - M_upper = numpy.concatenate((C, B.T), axis=1) - M_lower = numpy.concatenate((B, A), axis=1) - M = numpy.concatenate((M_upper, M_lower), axis=0) - del A, B, C - # M to WM, where W is the metric matrix [[-I, 0], [0, I]] - M[:nocc[0]*nocc[1]][:] *= -1.0 - - # =====================> solve for eigenpairs <====================== - exci, xy = scipy.linalg.eig(M) - exci = exci.real - xy = xy.T # Fortran to Python order - - # sort eigenpairs - idx = exci.argsort() - exci = exci[idx] - xy = xy[idx, :] - upprpa_orthonormalize_eigenvector('abab', nocc, exci, xy) - - sum_exci = numpy.sum(exci[nocc[0]*nocc[1]:]) - ec = sum_exci - trace_A - - return exci, xy, ec - - -def pprpa_print_direct_eigenvector( - pprpa, subspace, nocc, nvir, nocc_fro, thresh, hh_state, - pp_state, exci0, exci, xy): - """Print components of an eigenvector. - - Args: - pprpa (UppRPADirect object): unrestricted pprpa object. - subspace (str): subspace, 'aaaa', 'bbbb', or 'abab'. - nocc (int/tuple of int): number of occupied orbitals. - nvir (int/tuple of int): number of virtual orbitals. - nocc_fro (int/tuple of int): number of frozen occupied orbitals. - thresh (double): threshold to print a pair. - hh_state (int): number of interested hole-hole states. - pp_state (int): number of interested particle-particle states. - exci0 (double): lowest eigenvalue. - exci (double array): ppRPA eigenvalue. - xy (double ndarray): ppRPA eigenvector. - """ - if subspace == 'aaaa': - oo_dim = int((nocc - 1) * nocc / 2) - vv_dim = int((nvir - 1) * nvir / 2) - print("\n print U-ppRPA excitations: (alpha alpha, alpha alpha)\n") - elif subspace == 'bbbb': - oo_dim = int((nocc - 1) * nocc / 2) - vv_dim = int((nvir - 1) * nvir / 2) - print("\n print U-ppRPA excitations: (beta beta, beta beta)\n") - elif subspace == 'abab': - oo_dim = int(nocc[0] * nocc[1]) - vv_dim = int(nvir[0] * nvir[1]) - print("\n print U-ppRPA excitations: (alpha beta, alpha beta)\n") - else: - raise ValueError("Not recognized subspace: %s." % subspace) - - if subspace == 'aaaa' or subspace == 'bbbb': - tri_row_o, tri_col_o = numpy.tril_indices(nocc, -1) - tri_row_v, tri_col_v = numpy.tril_indices(nvir, -1) - - au2ev = 27.211386 - - # =====================> two-electron removal <====================== - for istate in range(min(hh_state, oo_dim)): - print("#%-d %s de-excitation: exci= %-12.4f eV 2e= %-12.4f eV" % - (istate + 1, subspace[:2], (exci[oo_dim-istate-1] - exci0) * au2ev, - exci[oo_dim-istate-1] * au2ev)) - if subspace == 'aaaa' or subspace == 'bbbb': - full = numpy.zeros(shape=[nocc, nocc], dtype=numpy.double) - full[tri_row_o, tri_col_o] = xy[oo_dim-istate-1][:oo_dim] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > thresh) - for i, j in pairs: - pprpa_print_a_pair(pprpa, is_pp=False, p=i+nocc_fro, q=j+nocc_fro, - percentage=full[i, j]) - - full = numpy.zeros(shape=[nvir, nvir], dtype=numpy.double) - full[tri_row_v, tri_col_v] = xy[oo_dim-istate-1][oo_dim:] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > thresh) - for a, b in pairs: - pprpa_print_a_pair(pprpa, is_pp=True, p=a+nocc_fro+nocc, - q=b+nocc_fro+nocc, percentage=full[a, b]) - - else: - full = xy[oo_dim-istate-1][:oo_dim].reshape(nocc[0], nocc[1]) - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > thresh) - for i, j in pairs: - pprpa_print_a_pair(pprpa, is_pp=False, p=i+nocc_fro[0], q=j+nocc_fro[1], - percentage=full[i, j]) - - full = xy[oo_dim-istate-1][oo_dim:].reshape(nvir[0], nvir[1]) - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > thresh) - for a, b in pairs: - pprpa_print_a_pair(pprpa, is_pp=True, p=a+nocc_fro[0]+nocc[0], - q=b+nocc_fro[1]+nocc[1], percentage=full[a, b]) - print("") - - # =====================> two-electron addition <===================== - for istate in range(min(pp_state, vv_dim)): - print("#%-d %s excitation: exci= %-12.4f eV 2e= %-12.4f eV" % - (istate + 1, subspace[:2], (exci[oo_dim+istate] - exci0) * au2ev, - exci[oo_dim+istate] * au2ev)) - if subspace == 'aaaa' or subspace == 'bbbb': - full = numpy.zeros(shape=[nocc, nocc], dtype=numpy.double) - full[tri_row_o, tri_col_o] = xy[oo_dim+istate][:oo_dim] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > thresh) - for i, j in pairs: - pprpa_print_a_pair(pprpa, is_pp=False, p=i+nocc_fro, q=j+nocc_fro, - percentage=full[i, j]) - - full = numpy.zeros(shape=[nvir, nvir], dtype=numpy.double) - full[tri_row_v, tri_col_v] = xy[oo_dim+istate][oo_dim:] - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > thresh) - for a, b in pairs: - pprpa_print_a_pair(pprpa, is_pp=True, p=a+nocc_fro+nocc, - q=b+nocc_fro+nocc, percentage=full[a, b]) - - else: - full = xy[oo_dim+istate][:oo_dim].reshape(nocc[0], nocc[1]) - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > thresh) - for i, j in pairs: - pprpa_print_a_pair( - pprpa, is_pp=False, p=i+nocc_fro[0], q=j+nocc_fro[1], - percentage=full[i, j]) - - full = xy[oo_dim+istate][oo_dim:].reshape(nvir[0], nvir[1]) - full = numpy.power(full, 2) - pairs = numpy.argwhere(full > thresh) - for a, b in pairs: - pprpa_print_a_pair( - pprpa, is_pp=True, p=a+nocc_fro[0]+nocc[0], - q=b+nocc_fro[1]+nocc[1], percentage=full[a, b]) - print("") - - return - - -def analyze_pprpa_direct(pprpa): - """Analyze ppRPA (direct diagonalization) excited states. - - Args: - pprpa (UppRPADirect): ppRPA object. - """ - logger.info(pprpa, '\nanalyze U-ppRPA results.') - nocc_fro = ( - pprpa.nocc[0] - pprpa.nocc_act[0], - pprpa.nocc[1] - pprpa.nocc_act[1]) - oo_dim_aa = int((pprpa.nocc_act[0] - 1) * pprpa.nocc_act[0] / 2) - oo_dim_bb = int((pprpa.nocc_act[1] - 1) * pprpa.nocc_act[1] / 2) - oo_dim_ab = int(pprpa.nocc_act[0] * pprpa.nocc_act[1]) - - exci_aa = pprpa.exci[0] - exci_bb = pprpa.exci[1] - exci_ab = pprpa.exci[2] - - exci0_list = [] - if exci_aa is not None: - logger.info(pprpa, '(alpha alpha, alpha alpha) results found.') - if pprpa.nelec == 'n-2': - exci0_list.append(exci_aa[oo_dim_aa]) - else: - exci0_list.append(exci_aa[oo_dim_aa - 1]) - if exci_bb is not None: - logger.info(pprpa, '(beta beta, beta beta) results found.') - if pprpa.nelec == 'n-2': - exci0_list.append(exci_bb[oo_dim_bb]) - else: - exci0_list.append(exci_bb[oo_dim_bb - 1]) - if exci_ab is not None: - logger.info(pprpa, '(alpha beta, alpha beta) results found.') - if pprpa.nelec == 'n-2': - exci0_list.append(exci_ab[oo_dim_ab]) - else: - exci0_list.append(exci_ab[oo_dim_ab - 1]) - - if pprpa.nelec == 'n-2': - exci0 = min(exci0_list) - else: - exci0 = max(exci0_list) - - if exci_aa is not None: - pprpa_print_direct_eigenvector( - pprpa, 'aaaa', pprpa.nocc_act[0], pprpa.nvir_act[0], nocc_fro[0], - pprpa.print_thresh, pprpa.hh_state, - pprpa.pp_state, exci0, exci_aa, pprpa.xy[0]) - if exci_bb is not None: - pprpa_print_direct_eigenvector( - pprpa, 'bbbb', pprpa.nocc_act[1], pprpa.nvir_act[1], nocc_fro[1], - pprpa.print_thresh, pprpa.hh_state, - pprpa.pp_state, exci0, exci_bb, pprpa.xy[1]) - if exci_ab is not None: - pprpa_print_direct_eigenvector( - pprpa, 'abab', pprpa.nocc_act, pprpa.nvir_act, nocc_fro, - pprpa.print_thresh, pprpa.hh_state, pprpa.pp_state, - exci0, exci_ab, pprpa.xy[2]) - - -class UppRPADirect(StreamObject): - def __init__( - self, mf, nocc_act=None, nvir_act=None, hh_state=5, pp_state=5, - nelec='n-2', print_thresh=0.1, auxbasis=None): - self.mol = mf.mol - self._scf = mf - self.verbose = self.mol.verbose - self.stdout = self.mol.stdout - self.max_memory = mf.max_memory - - # options - self.nocc_act = nocc_act # number of active occupied orbitals (alpha) - self.nvir_act = nvir_act # number of active virtual orbitals (alpha) - self.hh_state = hh_state # number of hole-hole states to print - self.pp_state = pp_state # number of particle-particle states to print - self.nelec = nelec # 'n-2' for N-2 system, 'n+2' for N+2 system - self.print_thresh = print_thresh # threshold to print component - self.auxbasis = auxbasis # auxiliary basis set to construct Lpq - - # ======================> internal flags <======================= - self.mu = None # chemical potential - self.nmo_act = None # number of active orbitals - self.mo_energy_act = None # orbital energy in active space - self.subspace = None # subspace(s) to perform ppRPA calculations - - # =========================> results <=========================== - self.ec = [None, None, None] # correlation energy [aaaa, bbbb, abab] - self.exci = [None, None, None] # two-electron addition energy [aaaa, bbbb, abab] - self.xy = [None, None, None] # ppRPA eigenvector [aaaa, bbbb, abab] - - # =============================================================== - # don't modify the following attributes, they are not input options - self._nocc = None # number of occupied orbitals - self._nmo = None # number of molecular orbitals - self.nvir = None # number of virtual orbitals - self.mo_energy = numpy.asarray(self._scf.mo_energy) # MO energies - self.Lpq = None # three-center density-fitting matrix in MO - self.mo_occ = self._scf.mo_occ - self.frozen = 0 - - return - - @property - def nocc(self): - return self.get_nocc() - @nocc.setter - def nocc(self, n): - self._nocc = n - - @property - def nmo(self): - return self.get_nmo() - @nmo.setter - def nmo(self, n): - self._nmo = n - - get_nocc = get_nocc - get_nmo = get_nmo - ao2mo = ao2mo - - analyze = analyze_pprpa_direct - - def check_parameter(self): - "Initialize and check options." - assert 0.0 < self.print_thresh < 1.0 - assert self.nelec in ["n-2", "n+2"] - - self.nvir = (self.nmo[0] - self.nocc[0], self.nmo[1] - self.nocc[1]) - - # adjust active space - if self.nocc_act is None: - self.nocc_act = self.nocc - else: - if isinstance(self.nocc_act, (int, numpy.integer)): - self.nocc_act = (self.nocc_act, self.nocc_act) - self.nocc_act = ( - min(self.nocc_act[0], self.nocc[0]), - min(self.nocc_act[1], self.nocc[1])) - if self.nvir_act is None: - self.nvir_act = self.nvir - else: - if isinstance(self.nvir_act, (int, numpy.integer)): - self.nvir_act = (self.nvir_act, self.nvir_act) - self.nvir_act = ( - min(self.nvir_act[0], self.nvir[0]), - min(self.nvir_act[1], self.nvir[1])) - - self.nmo_act = (self.nocc_act[0] + self.nvir_act[0], - self.nocc_act[1] + self.nvir_act[1]) - nocc = self.nocc - nocc_act = self.nocc_act - nvir_act = self.nvir_act - self.mo_energy_act = ( - self.mo_energy[0][nocc[0]-nocc_act[0]:nocc[0]+nvir_act[0]], - self.mo_energy[0][nocc[1]-nocc_act[1]:nocc[1]+nvir_act[1]]) - if self.mu is None: - self.mu = get_chemical_potential(self.nocc, self.mo_energy) - return - - def dump_flags(self): - log = logger.Logger(self.stdout, self.verbose) - # ====================> calculate dimensions <=================== - # (alpha, alpha) subspace - aavv_dim = int(self.nvir_act[0] * (self.nvir_act[0] + 1) / 2) - aaoo_dim = int(self.nocc_act[0] * (self.nocc_act[0] + 1) / 2) - # (alpha, beta) subspace - abvv_dim = int(self.nvir_act[0] * self.nvir_act[1]) - aboo_dim = int(self.nocc_act[0] * self.nocc_act[1]) - # (beta, beta) subspace - bbvv_dim = int(self.nvir_act[1] * (self.nvir_act[1] + 1) / 2) - bboo_dim = int(self.nocc_act[1] * (self.nocc_act[1] + 1) / 2) - - log.info('\n******** %s ********' % self.__class__) - log.info('nmo = %d (%d alpha, %d beta)', - self.nmo[0]+self.nmo[1], self.nmo[0], self.nmo[1]) - log.info('nocc = %d (%d alpha, %d beta), nvir = %d (%d alpha, %d beta)', - self.nocc[0] + self.nocc[1], self.nocc[0], self.nocc[1], - self.nvir[0] + self.nvir[1], self.nvir[0], self.nvir[1]) - log.info('nocc_act = %d (%d alpha, %d beta)', - self.nocc_act[0] + self.nocc_act[1], - self.nocc_act[0], self.nocc_act[1]) - log.info('nvir_act = %d (%d alpha, %d beta)', - self.nvir_act[0] + self.nvir_act[1], - self.nvir_act[0], self.nvir_act[1]) - log.info('for (alpha alpha, alpha alpha) subspace:') - log.info(' occ-occ dimension = %d vir-vir dimension = %d', - aaoo_dim, aavv_dim) - log.info('for (beta beta, beta beta) subspace:') - log.info(' occ-occ dimension = %d vir-vir dimension = %d', - bboo_dim, bbvv_dim) - log.info('for (alpha beta, alpha beta) subspace:') - log.info(' occ-occ dimension = %d vir-vir dimension = %d', - aboo_dim, abvv_dim) - log.info('interested hh state = %d', self.hh_state) - log.info('interested pp state = %d', self.pp_state) - log.info('ground state = %s', self.nelec) - log.info('print threshold = %.2f%%', self.print_thresh*100) - log.info('') - return - - def check_memory(self): - log = logger.Logger(self.stdout, self.verbose) - # ====================> calculate dimensions <=================== - # (alpha, alpha) subspace - aavv_dim = int(self.nvir_act[0] * (self.nvir_act[0] + 1) / 2) - aaoo_dim = int(self.nocc_act[0] * (self.nocc_act[0] + 1) / 2) - aafull_dim = aavv_dim + aaoo_dim - # (alpha, beta) subspace - abvv_dim = int(self.nvir_act[0] * self.nvir_act[1]) - aboo_dim = int(self.nocc_act[0] * self.nocc_act[1]) - abfull_dim = abvv_dim + aboo_dim - # (beta, beta) subspace - bbvv_dim = int(self.nvir_act[1] * (self.nvir_act[1] + 1) / 2) - bboo_dim = int(self.nocc_act[1] * (self.nocc_act[1] + 1) / 2) - bbfull_dim = bbvv_dim + bboo_dim - - full_dim = max(aafull_dim, abfull_dim, bbfull_dim) - - mem = (3 * full_dim * full_dim) * 8 / 1.0e6 - if mem < 1000: - log.info("U-ppRPA needs at least %d MB memory." % mem) - else: - log.info("U-ppRPA needs at least %.1f GB memory." % (mem / 1.0e3)) - return - - def kernel(self, subspace=['aa', 'bb', 'ab']): - """Run ppRPA direct diagonalization. - - Kwargs: - subspace (list of string): subspace(s) to run diagonalization, - 'aa' for (alpha, alpha, alpha, alpha), - 'bb' for (beta, beta, beta, beta), - and 'ab' for (alpha, beta, alpha, beta). - """ - self.subspace = subspace - self.check_parameter() - self.check_memory() - - cput0 = (logger.process_clock(), logger.perf_counter()) - self.dump_flags() - if self.Lpq is None: - self.Lpq = self.ao2mo() - if 'aa' in subspace: - aa_exci, aa_xy, aa_ec = diagonalize_pprpa_subspace_same_spin( - self.nocc_act[0], self.mo_energy_act[0], - self.Lpq[0], mu=self.mu) - else: - aa_exci = aa_xy = aa_ec = None - - if 'bb' in subspace: - bb_exci, bb_xy, bb_ec = diagonalize_pprpa_subspace_same_spin( - self.nocc_act[1], self.mo_energy_act[1], - self.Lpq[1], mu=self.mu) - else: - bb_exci = bb_xy = bb_ec = None - - if 'ab' in subspace: - ab_exci, ab_xy, ab_ec = diagonalize_pprpa_subspace_diff_spin( - self.nocc_act, self.mo_energy_act, self.Lpq, mu=self.mu) - else: - ab_exci = ab_xy = ab_ec = None - logger.timer(self, "ppRPA direct: %s" % subspace, *cput0) - - self.ec = [aa_ec, bb_ec, ab_ec] - self.exci = [aa_exci, bb_exci, ab_exci] - self.xy = [aa_xy, bb_xy, ab_xy] - - return - - def get_correlation(self): - """Get ppRPA correlation energy. - - Returns: - ec (double): ppRPA correlation energy. - """ - self.check_parameter() - - if self.Lpq is None: - self.Lpq = self.ao2mo() - if self.ec[0] is None: - cput0 = (logger.process_clock(), logger.perf_counter()) - self.exci[0], self.xy[0], self.ec[0] = \ - diagonalize_pprpa_subspace_same_spin( - self.nocc[0], self.mo_energy[0], self.Lpq[0], mu=self.mu) - logger.timer( - self, "ppRPA correlation energy: (alpha alpha, alpha alpha)", - *cput0) - if self.ec[1] is None: - cput0 = (logger.process_clock(), logger.perf_counter()) - self.exci[1], self.xy[1], self.ec[1] = \ - diagonalize_pprpa_subspace_same_spin( - self.nocc[1], self.mo_energy[1], self.Lpq[1], mu=self.mu) - logger.timer( - self, "ppRPA correlation energy: (beta beta, beta beta)", - *cput0) - if self.ec[2] is None: - cput0 = (logger.process_clock(), logger.perf_counter()) - self.exci[2], self.xy[2], self.ec[2] = \ - diagonalize_pprpa_subspace_diff_spin( - self.nocc, self.mo_energy, self.Lpq, mu=self.mu) - logger.timer( - self, "ppRPA correlation energy: (alpha beta, alpha beta)", - *cput0) - - return self.ec[0] + self.ec[1] + self.ec[2] - - def energy_tot(self): - """Get ppRPA total energy. - Totoal energy = Hartree-Fock energy + correlation energy from ppRPA. - - Returns: - e_tot (double); ppRPA total energy. - """ - mf = self._scf - assert mf.converged - hf_obj = mf if not isinstance(mf, KohnShamDFT) else mf.to_hf() - - dm = hf_obj.make_rdm1() - e_hf = hf_obj.energy_nuc() + hf_obj.energy_elec(dm=dm)[0] - ec = self.get_correlation() - e_tot = e_hf + ec - - return e_tot, e_hf, ec