diff --git a/.gitignore b/.gitignore index 8cfc943f..f020e4db 100644 --- a/.gitignore +++ b/.gitignore @@ -16,9 +16,20 @@ Makefile *.molden *.npy *.log +*.log.* *.dat *.out *.err *.swp *.chk +# assorted garbage +test_all.sh +*.status +*.rasscf.h5 +*.inp +*.timing +*.ods +*.tab +*.txt + diff --git a/__init__.py b/__init__.py index 940ec539..63a3bc11 100644 --- a/__init__.py +++ b/__init__.py @@ -1,2 +1,3 @@ # Let's see how complicated this has to get +from mrh.lib import patch_pyscf_sys_info diff --git a/debug/fci/1414.py b/debug/fci/1414.py new file mode 100644 index 00000000..35c22e50 --- /dev/null +++ b/debug/fci/1414.py @@ -0,0 +1,16 @@ +import numpy as np +from pyscf import lib, gto, ao2mo +from mrh.my_pyscf.fci import csf_solver +mol = gto.M (atom='Si 0 0 0', spin=14, charge=0, output='1414.log', verbose=lib.logger.DEBUG1) +fci = csf_solver (mol, smult=1) +norb = 14 +nelec = (7,7) +h1 = np.random.rand (14,14) +h1 += h1.T +h2 = np.random.rand (14,14,14,14) +h2 = ao2mo.restore (8, h2, norb) +h2 = ao2mo.restore (1, h2, norb) +e, ci = fci.kernel (h1, h2, norb, nelec) +print (e) + + diff --git a/debug/fci/lih_casscf25_sto3g.py b/debug/fci/lih_casscf25_sto3g.py new file mode 100644 index 00000000..16326157 --- /dev/null +++ b/debug/fci/lih_casscf25_sto3g.py @@ -0,0 +1,17 @@ +from pyscf import gto, scf, mcscf, lib + +mol = gto.M (atom='Li 0 0 0\nH 1.5 0 0', basis='sto-3g', symmetry=True, + output=__file__+'.log', verbose=lib.logger.INFO) +mf = scf.RHF (mol).run () +print (getattr (mf.mo_coeff, 'degen_mapping', None)) +mc0 = mcscf.CASSCF (mf, 5, 2) +mc0.kernel () + +print (getattr (mc0.mo_coeff, 'degen_mapping', None)) +print (getattr (mc0.fcisolver.orbsym, 'degen_mapping', None)) +mc1 = mcscf.CASCI (mf, 5, 2) +mc1.kernel (mc0.mo_coeff, ci0=mc0.ci) + + + + diff --git a/debug/lasscf/async_input.py b/debug/lasscf/async_input.py new file mode 100644 index 00000000..5ac5d0c2 --- /dev/null +++ b/debug/lasscf/async_input.py @@ -0,0 +1,25 @@ +import pyscf +#from geometry_generator import generator +from pyscf import gto, scf, tools, mcscf,lib +from mrh.my_pyscf.mcscf.lasscf_async import LASSCF +#from mrh.my_pyscf.mcscf.lasscf_sync_o0 import LASSCF +from pyscf.mcscf import avas +lib.logger.TIMER_LEVEL = lib.logger.INFO + +nfrags=2 +basis='sto-3g' +outputfile='async_1_6-31g_out.log' +#xyz=generator(nfrags) +xyz='butadiene.xyz' +mol=gto.M(atom=xyz,basis=basis,verbose=4,output=outputfile) +mf=scf.RHF(mol) +mf=mf.density_fit() +mf.run() +ncas,nelecas,guess_mo_coeff = avas.kernel(mf, ['C 2p']) +las=LASSCF(mf, list((2,)*nfrags),list((2,)*nfrags)) +frag_atom_list=[list(range(1+4*nfrag,3+4*nfrag)) for nfrag in range(nfrags)] +mo_coeff=las.set_fragments_(frag_atom_list, guess_mo_coeff) +las.kernel(mo_coeff) +#mf.mo_coeff=las.mo_coeff +#mf.analyze() + diff --git a/debug/lasscf/butadiene.xyz b/debug/lasscf/butadiene.xyz new file mode 100644 index 00000000..9f5a1f88 --- /dev/null +++ b/debug/lasscf/butadiene.xyz @@ -0,0 +1,12 @@ +10 + +H -0.943967690125 0.545000000000 0.000000000000 +C 0.000000000000 0.000000000000 0.000000000000 +C 1.169134295109 0.675000000000 0.000000000000 +H 0.000000000000 -1.090000000000 0.000000000000 +H 1.169134295109 1.765000000000 0.000000000000 +C 2.459512146748 -0.070000000000 0.000000000000 +C 3.628646441857 0.605000000000 0.000000000000 +H 2.459512146748 -1.160000000000 0.000000000000 +H 3.628646441857 1.695000000000 0.000000000000 +H 4.572614131982 0.060000000000 0.000000000000 diff --git a/debug/lasscf/c2h6n4.xyz b/debug/lasscf/c2h6n4.xyz new file mode 100644 index 00000000..45b39bb4 --- /dev/null +++ b/debug/lasscf/c2h6n4.xyz @@ -0,0 +1,12 @@ +H -0.145525 2.534624 1.196098 +N 0.586282 2.685058 0.454251 +N 0.586282 1.768962 -0.376701 +C -0.376894 0.666619 -0.222182 +H -0.986448 0.788749 0.689250 +H -1.039496 0.688384 -1.095107 +C 0.376894 -0.666619 -0.222182 +H 0.986448 -0.788749 0.689250 +H 1.039496 -0.688384 -1.095107 +N -0.586282 -1.768962 -0.376701 +N -0.586282 -2.685058 0.454251 +H 0.145525 -2.534624 1.196098 diff --git a/debug/lasscf/c2h6n4_struct.py b/debug/lasscf/c2h6n4_struct.py new file mode 100644 index 00000000..1d996f68 --- /dev/null +++ b/debug/lasscf/c2h6n4_struct.py @@ -0,0 +1,28 @@ +import numpy as np +from pyscf import gto + +def structure (dnn1=0, dnn2=0, basis='6-31g', symmetry=False): + with open ('c2h6n4.xyz', 'r') as f: + equilgeom = f.read () + mol = gto.M (atom = equilgeom, basis = basis, symmetry=True, unit='au') + atoms = tuple(mol.atom_symbol (i) for i in range (mol.natm)) + coords = mol.atom_coords () + + idx_nn = np.asarray ([[0, 1], [10, 11]]) + nn_vec = coords[[1,10],:] - coords[[2,9],:] + nn_equil = np.mean (np.linalg.norm (nn_vec, axis=1)) + nn_vec /= np.linalg.norm (nn_vec, axis=1)[:,np.newaxis] + delta_coords = np.zeros_like (coords) + delta_coords[idx_nn[0],:] = np.broadcast_to (nn_vec[0,:], (2,3)) + delta_coords[idx_nn[1],:] = np.broadcast_to (nn_vec[1,:], (2,3)) + scale = np.zeros (12, dtype=coords.dtype) + scale[idx_nn[0]] = dnn1 + scale[idx_nn[1]] = dnn2 + newcoords = coords + scale[:,np.newaxis] * delta_coords + carts = [[atoms[i]] + list(newcoords[i,:]) for i in range(12)] + dummymol = gto.M (atom = carts, basis=basis, symmetry=True) + newcoords = dummymol.atom_coords () + carts = [[atoms[i]] + list(newcoords[i,:]) for i in range(12)] + return gto.M (atom = carts, basis=basis, symmetry=symmetry, unit='au') + + diff --git a/debug/lassi/debug_c2h4n4.py b/debug/lassi/debug_c2h4n4.py new file mode 100644 index 00000000..24ade981 --- /dev/null +++ b/debug/lassi/debug_c2h4n4.py @@ -0,0 +1,186 @@ +#!/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. +import os +import copy +import unittest +import numpy as np +from scipy import linalg +from pyscf import lib, gto, scf, dft, fci, mcscf, df +from pyscf.tools import molden +from mrh.tests.lasscf.c2h4n4_struct import structure as struct +from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF +from mrh.my_pyscf.lassi.lassi import roots_make_rdm12s, root_make_rdm12s, make_stdm12s, ham_2q +from mrh.my_pyscf.lassi import LASSI +topdir = os.path.abspath (os.path.join (__file__, '../../../tests/lassi/')) + +def setUpModule (): + global las, lsi, rdm1s, rdm2s + dr_nn = 2.0 + mol = struct (dr_nn, dr_nn, '6-31g', symmetry=False) + mol.verbose = 0 #lib.logger.DEBUG + mol.output = '/dev/null' #'test_c2h4n4.log' + mol.spin = 0 + mol.build () + mf = scf.RHF (mol).run () + las = LASSCF (mf, (4,4), (4,4), spin_sub=(1,1)) + las.state_average_(weights=[1.0/5.0,]*5, + spins=[[0,0],[0,0],[2,-2],[-2,2],[2,2]], + smults=[[1,1],[3,3],[3,3],[3,3],[3,3]]) + las.frozen = list (range (las.mo_coeff.shape[-1])) + ugg = las.get_ugg () + las.mo_coeff = np.loadtxt (os.path.join (topdir, 'test_c2h4n4_mo.dat')) + las.ci = ugg.unpack (np.loadtxt (os.path.join (topdir, 'test_c2h4n4_ci.dat')))[1] + lroots = 2 * np.ones ((2,7), dtype=int) + lroots[:,1] = 1 # <- that rootspace is responsible for spin contamination + las.conv_tol_grad = 1e-8 + las.lasci (lroots=lroots) + # TODO: potentially save and load CI vectors + # requires extending features of ugg + #las.set (conv_tol_grad=1e-8).run () + #np.savetxt ('test_c2h4n4_mo.dat', las.mo_coeff) + #np.savetxt ('test_c2h4n4_ci.dat', ugg.pack (las.mo_coeff, las.ci)) + #las.e_states = las.energy_nuc () + las.states_energy_elec () + lsi = LASSI (las).run () + rdm1s, rdm2s = roots_make_rdm12s (las, las.ci, lsi.si) + las = LASSCF (mf, (4,2,4), ((4,0),(1,1),(0,4)), spin_sub=(5,3,5)) + ci_quin = np.array ([[1.0,],]) + ci_sing = np.diag ([1,-1])[::-1] / np.sqrt (2) + las.ci = [[ci_quin], [ci_sing], [ci_quin.copy ()]] + las.e_states = las.energy_nuc () + las.states_energy_elec () + +def tearDownModule(): + global las, lsi, rdm1s, rdm2s + mol = lsi._las.mol + mol.stdout.close () + del las, lsi, rdm1s, rdm2s + +class KnownValues(unittest.TestCase): + #def test_evals (self): + # self.assertAlmostEqual (lib.fp (lsi.e_roots), 71.48678303162376, 6) + + #def test_si (self): + # # Arbitrary signage in both the SI and CI vector, sadly + # # Actually this test seems really inconsistent overall... + # dens = lsi.si.conj () * lsi.si + # self.assertAlmostEqual (lib.fp (np.diag (dens)), 1.2321935065030327, 4) + + #def test_nelec (self): + # for ix, ne in enumerate (lsi.nelec): + # with self.subTest(ix): + # if ix in (1,5,8,11): + # self.assertEqual (ne, (6,2)) + # else: + # self.assertEqual (ne, (4,4)) + + #def test_s2 (self): + # s2_array = np.zeros (16) + # quintets = [1,2,5,8,11] + # for ix in quintets: s2_array[ix] = 6 + # triplets = [3,6,7,9,10,12,13] + # for ix in triplets: s2_array[ix] = 2 + # self.assertAlmostEqual (lib.fp (lsi.s2), lib.fp (s2_array), 3) + + #def test_tdms (self): + # las, si = lsi._las, lsi.si + # stdm1s, stdm2s = make_stdm12s (las) + # nelec = float (sum (las.nelecas)) + # for ix in range (stdm1s.shape[0]): + # d1 = stdm1s[ix,...,ix].sum (0) + # d2 = stdm2s[ix,...,ix].sum ((0,3)) + # with self.subTest (root=ix): + # self.assertAlmostEqual (np.trace (d1), nelec, 9) + # self.assertAlmostEqual (np.einsum ('ppqq->',d2), nelec*(nelec-1), 9) + # rdm1s_test = np.einsum ('ar,asijb,br->rsij', si.conj (), stdm1s, si) + # rdm2s_test = np.einsum ('ar,asijtklb,br->rsijtkl', si.conj (), stdm2s, si) + # self.assertAlmostEqual (lib.fp (rdm1s_test), lib.fp (rdm1s), 9) + # self.assertAlmostEqual (lib.fp (rdm2s_test), lib.fp (rdm2s), 9) + + #def test_rdms (self): + # las, e_roots = lsi._las, lsi.e_roots + # h0, h1, h2 = ham_2q (las, las.mo_coeff) + # d1_r = rdm1s.sum (1) + # d2_r = rdm2s.sum ((1, 4)) + # nelec = float (sum (las.nelecas)) + # for ix, (d1, d2) in enumerate (zip (d1_r, d2_r)): + # with self.subTest (root=ix): + # self.assertAlmostEqual (np.trace (d1), nelec, 9) + # self.assertAlmostEqual (np.einsum ('ppqq->',d2), nelec*(nelec-1), 9) + # e_roots_test = h0 + np.tensordot (d1_r, h1, axes=2) + np.tensordot (d2_r, h2, axes=4) / 2 + # for e1, e0 in zip (e_roots_test, e_roots): + # self.assertAlmostEqual (e1, e0, 8) + + #def test_singles_constructor (self): + # from mrh.my_pyscf.lassi.states import all_single_excitations + # las2 = all_single_excitations (lsi._las) + # las2.check_sanity () + # # Meaning of tuple: (na+nb,smult) + # # from state 0, (1+2,2)&(3+2,2) & a<->b & l<->r : 4 states + # # from state 1, smults=((2,4),(4,2),(4,4)) permutations of above: 12 additional states + # # from states 2 and 3, (4+1,4)&(0+3,4) & a<->b & l<->r : 4 additional states + # # from state 4, (2+1,2)&(4+1,4) & (2+1,4)&(4+1,4) + # # & (3+0,4)&(3+2,2) & (3+0,4)&(3+2,4) & l<->r : 8 additional states + # # 5 + 4 + 12 + 4 + 8 = 33 + # self.assertEqual (las2.nroots, 33) + + #def test_spin_shuffle (self): + # from mrh.my_pyscf.lassi.states import spin_shuffle, spin_shuffle_ci + # mf = lsi._las._scf + # las3 = spin_shuffle (las) + # las3.check_sanity () + # # The number of states is the number of graphs connecting one number + # # in each row which sum to zero: + # # -2 -1 0 +1 +2 + # # -1 0 +1 + # # -2 -1 0 +1 +2 + # # For the first two rows, paths which sum to -3 and +3 are immediately + # # excluded. Two paths connecting the first two rows each sum to -2 and +2 + # # and three paths each sum to -1, 0, +1. Each partial sum then has one + # # remaining option to complete the path, so + # # 2 + 3 + 3 + 3 + 2 = 13 + # with self.subTest ("state construction"): + # self.assertEqual (las3.nroots, 13) + # las3.ci = spin_shuffle_ci (las3, las3.ci) + # lsi2 = LASSI (las3).run () + # errvec = lsi2.s2 - np.around (lsi2.s2) + # with self.subTest ("CI vector rotation"): + # self.assertLess (np.amax (np.abs (errvec)), 1e-8) + + #def test_lassis (self): + # from mrh.my_pyscf.lassi.lassis import LASSIS + # las1 = LASSCF (las._scf, (4,4), (4,4), spin_sub=(1,1)) + # las1.mo_coeff = las.mo_coeff + # las1.lasci () + # lsis = LASSIS (las1).run (max_cycle_macro=1) + + def test_lassis_slow (self): + from mrh.my_pyscf.lassi.lassis import LASSIS + # TODO: optimize implementation and eventually merge with test_lassis + mol = struct (2.0, 2.0, '6-31g', symmetry=False) + mol.output = 'debug_c2h4n4.log' + mol.verbose = lib.logger.DEBUG + mol.spin = 8 + mol.build () + mf = scf.RHF (mol).run () + las1 = LASSCF (mf, (5,5), ((3,2),(2,3)), spin_sub=(2,2)) + mo_coeff = las1.localize_init_guess ((list (range (5)), list (range (5,10)))) + las1.kernel (mo_coeff) + lsis = LASSIS (las1, opt=1).run () + self.assertAlmostEqual (lsis.e_roots[0], -295.52103109, 7) + self.assertTrue (lsis.converged) + +if __name__ == "__main__": + print("Full Tests for SA-LASSI of c2h4n4 molecule") + unittest.main() + diff --git a/debug/lassi/debug_lassis_targets_slow.py b/debug/lassi/debug_lassis_targets_slow.py new file mode 100644 index 00000000..75181a52 --- /dev/null +++ b/debug/lassi/debug_lassis_targets_slow.py @@ -0,0 +1,264 @@ +#!/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. +import os +import copy +import unittest +import numpy as np +from pyscf import gto, scf, mcscf +from pyscf.lib import chkfile +from pyscf.mcscf import avas +from pyscf.data import nist +from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF +from mrh.my_pyscf.mcscf import lasscf_async +from mrh.tests.lasscf.c2h4n4_struct import structure as struct +from mrh.my_pyscf import lassi +topdir = os.path.abspath (os.path.join (__file__, '..')) + +au2cm = nist.HARTREE2J / nist.PLANCK / nist.LIGHT_SPEED_SI * 1e-2 +def yamaguchi (e_roots, s2, nelec=6): + '''The Yamaguchi formula for nelec unpaired electrons''' + hs = (nelec/2.0) * ((nelec/2.0)+1.0) + ls = hs - np.floor (hs) + idx = np.argsort (e_roots) + e_roots = e_roots[idx] + s2 = s2[idx] + idx_hs = (np.around (s2) == np.around (hs)) + assert (np.count_nonzero (idx_hs)), 'high-spin ground state not found' + idx_hs = np.where (idx_hs)[0][0] + e_hs = e_roots[idx_hs] + idx_ls = (np.around (s2) == np.around (ls)) + assert (np.count_nonzero (idx_ls)), 'low-spin ground state not found' + idx_ls = np.where (idx_ls)[0][0] + e_ls = e_roots[idx_ls] + j = (e_ls - e_hs) / (hs - ls) + return j*au2cm + +def setUpModule (): + pass + +def tearDownModule(): + pass + +class KnownValues(unittest.TestCase): + + def test_c2h4n4_3frag (self): + # str + mol = struct (2.0, 2.0, '6-31g', symmetry=False) + mol.output = 'c2h4n4_3frag_str.log'#'/dev/null' + mol.verbose = 5#0 + mol.spin = 8 + mol.build () + mf = scf.RHF (mol).run () + las = LASSCF (mf, (4,2,4), ((2,2),(1,1),(2,2)), spin_sub=(1,1,1)) + mo_coeff = las.localize_init_guess ([[0,1,2],[3,4,5,6],[7,8,9]]) + las.kernel (mo_coeff) + lsi = lassi.LASSIS (las).run () + e_str = lsi.e_roots[0] + with self.subTest ('str converged'): + self.assertTrue (lsi.converged) + # equil + mol = struct (0.0, 0.0, '6-31g', symmetry=False) + mol.spin = 0 + mol.output = 'c2h4n4_3frag_equil.log'#'/dev/null' + mol.verbose = 5#0 + mol.build () + mf = scf.RHF (mol).run () + las = LASSCF (mf, (4,2,4), ((2,2),(1,1),(2,2)), spin_sub=(1,1,1)) + mo_coeff = las.sort_mo ([7,8,16,18,22,23,24,26,33,34]) + mo_coeff = las.localize_init_guess ([[0,1,2],[3,4,5,6],[7,8,9]], mo_coeff=mo_coeff) + las.kernel (mo_coeff) + lsi = lassi.LASSIS (las).run () + e_equil = lsi.e_roots[0] + with self.subTest ('equil converged'): + self.assertTrue (lsi.converged) + # test + de = 1000 * (e_str - e_equil) + self.assertAlmostEqual (de, 208.21775437138967, 1) + + def test_c2h4n4_2frag (self): + # str + mol = struct (2.0, 2.0, '6-31g', symmetry=False) + mol.output = 'c2h4n4_2frag_str.log'#'/dev/null' + mol.verbose = 5#0 + mol.spin = 8 + mol.build () + mf = scf.RHF (mol).run () + las = lasscf_async.LASSCF (mf, (5,5), ((3,2),(2,3))) + las = las.state_average ([0.5,0.5], + spins=[[1,-1],[-1,1]], + smults=[[2,2],[2,2]], + charges=[[0,0],[0,0]], + wfnsyms=[[0,0],[0,0]]) + mo = las.set_fragments_((list (range (5)), list (range (5,10)))) + las.kernel (mo) + mo_coeff = las.mo_coeff + las = lasscf_async.LASSCF (mf, (5,5), ((3,2),(2,3))) + las.lasci_(mo_coeff) + lsi = lassi.LASSIS (las).run () + e_str = lsi.e_roots[0] + with self.subTest ('str converged'): + self.assertTrue (lsi.converged) + # equil + mol = struct (0.0, 0.0, '6-31g', symmetry=False) + mol.spin = 0 + mol.output = 'c2h4n4_2frag_equil.log'#'/dev/null' + mol.verbose = 5#0 + mol.build () + mf = scf.RHF (mol).run () + las = lasscf_async.LASSCF (mf, (5,5), ((3,2),(2,3))) + las = las.state_average ([0.5,0.5], + spins=[[1,-1],[-1,1]], + smults=[[2,2],[2,2]], + charges=[[0,0],[0,0]], + wfnsyms=[[0,0],[0,0]]) + mo = las.sort_mo ([7,8,16,18,22,23,24,26,33,34]) + mo = las.set_fragments_((list (range (5)), list (range (5,10))), mo) + las.kernel (mo) + mo_coeff = las.mo_coeff + las = lasscf_async.LASSCF (mf, (5,5), ((3,2),(2,3))) + las.lasci_(mo_coeff) + lsi = lassi.LASSIS (las).run () + e_equil = lsi.e_roots[0] + with self.subTest ('equil converged'): + self.assertTrue (lsi.converged) + # test + de = 1000 * (e_str - e_equil) + self.assertAlmostEqual (de, 191.185141573740, 1) + + def test_kremer_cr2_model (self): + xyz='''Cr -1.320780000000 0.000050000000 -0.000070000000 + Cr 1.320770000000 0.000050000000 -0.000070000000 + O 0.000000000000 -0.165830000000 1.454680000000 + O 0.000000000000 1.342770000000 -0.583720000000 + O 0.000000000000 -1.176830000000 -0.871010000000 + H 0.000020000000 0.501280000000 2.159930000000 + H 0.000560000000 1.618690000000 -1.514480000000 + H -0.000440000000 -2.120790000000 -0.644130000000 + N -2.649800000000 -1.445690000000 0.711420000000 + H -2.186960000000 -2.181980000000 1.244400000000 + H -3.053960000000 -1.844200000000 -0.136070000000 + H -3.367270000000 -1.005120000000 1.287210000000 + N -2.649800000000 1.339020000000 0.896300000000 + N -2.649800000000 0.106770000000 -1.607770000000 + H -3.367270000000 -0.612160000000 -1.514110000000 + H -3.053960000000 0.804320000000 1.665160000000 + N 2.649800000000 -1.445680000000 0.711420000000 + N 2.649790000000 1.339030000000 0.896300000000 + N 2.649800000000 0.106780000000 -1.607770000000 + H -2.186970000000 2.168730000000 1.267450000000 + H -3.367270000000 1.617370000000 0.226860000000 + H -2.186960000000 0.013340000000 -2.511900000000 + H -3.053970000000 1.039980000000 -1.529140000000 + H 2.186960000000 -2.181970000000 1.244400000000 + H 3.053960000000 -1.844190000000 -0.136080000000 + H 3.367270000000 -1.005100000000 1.287200000000 + H 2.186950000000 2.168740000000 1.267450000000 + H 3.053960000000 0.804330000000 1.665160000000 + H 3.367260000000 1.617380000000 0.226850000000 + H 2.186960000000 0.013350000000 -2.511900000000 + H 3.053960000000 1.039990000000 -1.529140000000 + H 3.367270000000 -0.612150000000 -1.514110000000''' + basis = {'C': 'sto-3g','H': 'sto-3g','O': 'sto-3g','N': 'sto-3g','Cr': 'cc-pvdz'} + mol = gto.M (atom=xyz, spin=6, charge=3, basis=basis, + verbose=5, output='kremer.log') + #verbose=0, output='/dev/null') + mf = scf.ROHF(mol) + mf.chkfile = 'test_lassis_targets_slow.kremer_cr2_model.chk' + mf.kernel () + las = LASSCF(mf,(6,6),((3,0),(0,3)),spin_sub=(4,4)) + try: + mo_coeff = chkfile.load (las.chkfile, 'las')['mo_coeff'] + except (OSError, TypeError, KeyError) as e: + ncas_avas, nelecas_avas, mo_coeff = avas.kernel (mf, ['Cr 3d', 'Cr 4d'], minao=mol.basis) + mc_avas = mcscf.CASCI (mf, ncas_avas, nelecas_avas) + mo_list = mc_avas.ncore + np.array ([5,6,7,8,9,10,15,16,17,18,19,20]) + mo_coeff = las.sort_mo (mo_list, mo_coeff) + mo_coeff = las.localize_init_guess (([0],[1]), mo_coeff) + las = lassi.states.spin_shuffle (las) # generate direct-exchange states + las.weights = [1.0/las.nroots,]*las.nroots # set equal weights + nroots_ref = las.nroots + las.kernel (mo_coeff) # optimize orbitals + mo_coeff = las.mo_coeff + las = LASSCF(mf,(6,6),((3,0),(0,3)),spin_sub=(4,4)) + las.lasci_(mo_coeff) + lsi = lassi.LASSIS (las).run () + with self.subTest('convergence'): + self.assertTrue (lsi.converged) + self.assertAlmostEqual (yamaguchi (lsi.e_roots, lsi.s2, 6), -12.45, 2) + + def test_alfefe (self): + xyz='''O -2.2201982441 0.3991903003 1.6944716989 + O -1.6855532446 -1.7823063217 1.4313995072 + C -2.2685178651 -0.8319550379 1.983951274 + H -2.9133420167 -1.0767285885 2.8437868053 + O 1.3882880809 2.0795561562 -1.3470856753 + O -0.7599088595 2.5809236347 -0.849227704 + C 0.3465674686 2.753895325 -1.438835699 + H 0.3723796145 3.6176996598 -2.1230911503 + O 1.0469609081 -2.6425342 0.9613246722 + O -0.3991903003 2.2201982441 1.6944716989 + O 1.7823063217 1.6855532446 1.4313995072 + C 0.8319550379 2.2685178651 1.983951274 + H 1.0767285885 2.9133420167 2.8437868053 + O -2.0795561562 -1.3882880809 -1.3470856753 + O -2.5809236347 0.7599088595 -0.849227704 + C -2.753895325 -0.3465674686 -1.438835699 + H -3.6176996598 -0.3723796145 -2.1230911503 + Fe -0.3798741893 -1.7926033629 -0.2003377303 + O 0.6422231803 -2.237796553 -1.8929145176 + Fe 1.7926033629 0.3798741893 -0.2003377303 + O 2.237796553 -0.6422231803 -1.8929145176 + O 2.6425342 -1.0469609081 0.9613246722 + Al -1.2362716421 1.2362716421 0.350650148 + O -0.0449501954 0.0449501954 0.0127621413 + C 1.645528685 -1.645528685 -2.3571124654 + H 2.0569062984 -2.0569062984 -3.2945712916 + C 2.1609242811 -2.1609242811 1.2775841868 + H 2.7960805433 -2.7960805433 1.9182443024''' + basis = {'C': 'sto-3g','H': 'sto-3g','O': 'sto-3g','Al': 'cc-pvdz','Fe': 'cc-pvdz'} + mol = gto.M (atom=xyz, spin=9, charge=0, basis=basis, max_memory=10000, + verbose=5, output='alfefe.log') + #verbose=0, output='/dev/null') + mf = scf.ROHF(mol) + mf.init_guess='chk' + mf.chkfile='test_lassis_targets_slow.alfefe.chk' + mf = mf.density_fit() + mf.max_cycle=100 + mf.kernel() + las = LASSCF (mf, (5,5), ((5,1),(5,0)), spin_sub=(5,6)) + try: + las.load_chk_() + las.kernel () + except (OSError, TypeError, KeyError) as e: + ncas, nelecas, mo_coeff = avas.kernel (mf, ['Fe 3d'], openshell_option=3) + mo_coeff = las.localize_init_guess (([17],[19]), mo_coeff) + las.kernel (mo_coeff) + with self.subTest ('LASSCF convergence'): + self.assertTrue (las.converged) + with self.subTest ('same LASSCF result'): + self.assertAlmostEqual (las.e_tot, -3955.98148841553, 6) + mf.chkfile = None # prevent the spins flips down there from messing things up + las2 = LASSCF (mf, (5,5), ((1,5),(5,0)), spin_sub=(5,6)) + las2.lasci_(las.mo_coeff) + lsi = lassi.LASSIS (las2).run () + with self.subTest('LASSI convergence'): + self.assertTrue (lsi.converged) + self.assertAlmostEqual (yamaguchi (lsi.e_roots, lsi.s2, 9), -4.40, 2) + + +if __name__ == "__main__": + print("Full Tests for SA-LASSI of c2h4n4 molecule") + unittest.main() + diff --git a/debug/lassi/debug_opt57_slow.py b/debug/lassi/debug_opt57_slow.py new file mode 100644 index 00000000..ad484d8f --- /dev/null +++ b/debug/lassi/debug_opt57_slow.py @@ -0,0 +1,163 @@ +#!/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. + +import copy +import unittest +import numpy as np +from scipy import linalg +from copy import deepcopy +from itertools import product +from pyscf import lib, gto, scf, dft, fci, mcscf, df +from pyscf.tools import molden +from pyscf.fci import cistring +from pyscf.fci.direct_spin1 import _unpack_nelec +from mrh.tests.lasscf.c2h4n4_struct import structure as struct +from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF +from mrh.my_pyscf.mcscf.lasci import get_space_info +from mrh.my_pyscf.lassi.lassi import roots_make_rdm12s, make_stdm12s, ham_2q +from mrh.my_pyscf.lassi.citools import get_lroots, get_rootaddr_fragaddr +from mrh.my_pyscf.lassi import op_o0 +from mrh.my_pyscf.lassi import op_o1 + +def setUpModule (): + global mol, mf, las, nstates, nelec_frs, si, orbsym, wfnsym + # Build crazy state list + states = {'charges': [[0,0,0],], + 'spins': [[0,0,0],], + 'smults': [[1,1,1],], + 'wfnsyms': [[0,0,0],]} + states1 = {'charges': [[-1,1,0],[-1,1,0],[1,-1,0],[1,-1,0],[0,1,-1],[0,1,-1],[0,-1,1],[0,-1,1]], + 'spins': [[-1,1,0],[1,-1,0],[-1,1,0],[1,-1,0],[0,1,-1],[0,-1,1],[0,1,-1],[0,-1,1]], + 'smults': [[2,2,1], [2,2,1], [2,2,1], [2,2,1], [1,2,2], [1,2,2], [1,2,2], [1,2,2]], + 'wfnsyms': [[1,1,0], [1,1,0], [1,1,0], [1,1,0], [0,1,1], [0,1,1], [0,1,1], [0,1,1]]} + states2 = {'charges': [[0,0,0],]*6, + 'spins': [[2,-2,0],[0,0,0],[-2,2,0],[0,2,-2],[0,0,0],[0,-2,2]], + 'smults': [[3,3,1], [3,3,1],[3,3,1], [1,3,3], [1,3,3],[1,3,3]], + 'wfnsyms': [[0,0,0],]*6} + states3 = {'charges': [[-1,2,-1],[-1,2,-1],[1,-2,1],[1,-2,1],[-1,0,1],[-1,0,1],[1,0,-1],[1,0,-1]], + 'spins': [[1,0,-1], [-1,0,1], [1,0,-1],[-1,0,1],[1,0,-1],[-1,0,1],[1,0,-1],[-1,0,1]], + 'smults': [[2,1,2], [2,1,2], [2,1,2], [2,1,2], [2,1,2], [2,1,2], [2,1,2], [2,1,2]], + 'wfnsyms': [[1,0,1],]*8} + states4 = {'charges': [[0,0,0],]*10, + 'spins': [[-2,0,2],[0,0,0],[2,0,-2],[-2,0,2],[0,0,0],[2,0,-2],[2,-2,0],[-2,2,0],[0,2,-2],[0,-2,2]], + 'smults': [[3,1,3], [3,1,3],[3,1,3], [3,3,3], [3,3,3],[3,3,3], [3,3,3], [3,3,3], [3,3,3], [3,3,3]], + 'wfnsyms': [[0,0,0],]*10} + states5 = {'charges': [[-1,1,0],[-1,1,0], [-1,1,0],[-1,1,0],[1,-1,0],[1,-1,0], [1,-1,0],[1,-1,0]], + 'spins': [[1,1,-2],[-1,-1,2],[1,-1,0],[-1,1,0],[1,1,-2],[-1,-1,2],[1,-1,0],[-1,1,0]], + 'smults': [[2,2,3], [2,2,3], [2,2,3], [2,2,3], [2,2,3], [2,2,3], [2,2,3], [2,2,3]], + 'wfnsyms': [[1,1,0],]*8} + states6 = deepcopy (states5) + states7 = deepcopy (states5) + lroots = np.ones ((3, 57), dtype=int) + offs = [0,] + for field in ('charges', 'spins', 'smults', 'wfnsyms'): + states6[field] = [[row[1], row[2], row[0]] for row in states5[field]] + states7[field] = [[row[2], row[0], row[1]] for row in states5[field]] + for d in [states1, states2, states3, states4, states5, states6, states7]: + offs.append (len (states['charges'])) + for field in ('charges', 'spins', 'smults', 'wfnsyms'): + states[field] = states[field] + d[field] + lroots[:,offs[0]] = [2, 2, 2] + lroots[:,offs[1]] = [2, 2, 2] + lroots[:,offs[2]] = [2, 1, 2] + lroots[:,offs[3]] = [2, 1, 2] + lroots[:,offs[4]] = [2, 2, 2] + lroots[:,offs[5]] = [2, 2, 2] + weights = [1.0,] + [0.0,]*56 + nroots = 57 + nstates = 91 + # End building crazy state list + + dr_nn = 2.0 + mol = struct (dr_nn, dr_nn, '6-31g', symmetry='Cs') + mol.verbose = 0 #lib.logger.INFO + mol.output = '/dev/null' #'test_lassi_op.log' + mol.spin = 0 + mol.build () + mf = scf.RHF (mol).run () + las = LASSCF (mf, (4,2,4), (4,2,4)) + las.state_average_(weights=weights, **states) + las.mo_coeff = las.localize_init_guess ((list (range (3)), + list (range (3,7)), list (range (7,10))), mf.mo_coeff) + las.ci = las.get_init_guess_ci (las.mo_coeff, las.get_h2eff (las.mo_coeff)) + nelec_frs = np.array ( + [[_unpack_nelec (fcibox._get_nelec (solver, nelecas)) for solver in fcibox.fcisolvers] + for fcibox, nelecas in zip (las.fciboxes, las.nelecas_sub)] + ) + ndet_frs = np.array ( + [[[cistring.num_strings (las.ncas_sub[ifrag], nelec_frs[ifrag,iroot,0]), + cistring.num_strings (las.ncas_sub[ifrag], nelec_frs[ifrag,iroot,1])] + for iroot in range (las.nroots)] for ifrag in range (las.nfrags)] + ) + np.random.seed (1) + for iroot in range (las.nroots): + for ifrag in range (las.nfrags): + lroots_r = lroots[ifrag,iroot] + ndet_s = ndet_frs[ifrag,iroot] + ci = np.random.rand (lroots_r, ndet_s[0], ndet_s[1]) + ci /= linalg.norm (ci.reshape (lroots_r,-1), axis=1)[:,None,None] + if lroots_r==1: ci=ci[0] + las.ci[ifrag][iroot] = ci + orbsym = getattr (las.mo_coeff, 'orbsym', None) + if orbsym is None and callable (getattr (las, 'label_symmetry_', None)): + orbsym = las.label_symmetry_(las.mo_coeff).orbsym + if orbsym is not None: + orbsym = orbsym[las.ncore:las.ncore+las.ncas] + wfnsym = 0 + #las.lasci (lroots=lroots) + rand_mat = np.random.rand (nstates,nstates) + rand_mat += rand_mat.T + e, si = linalg.eigh (rand_mat) + +def tearDownModule(): + global mol, mf, las, nstates, nelec_frs, si, orbsym, wfnsym + mol.stdout.close () + del mol, mf, las, nstates, nelec_frs, si, orbsym, wfnsym + +class KnownValues(unittest.TestCase): + def test_stdm12s (self): + d12_o0 = make_stdm12s (las, opt=0) + d12_o1 = make_stdm12s (las, opt=1) + rootaddr, fragaddr = get_rootaddr_fragaddr (get_lroots (las.ci)) + for r in range (2): + for i, j in product (range (nstates), repeat=2): + with self.subTest (rank=r+1, idx=(i,j), spaces=(rootaddr[i], rootaddr[j]), + envs=(list(fragaddr[:,i]),list(fragaddr[:,j]))): + self.assertAlmostEqual (lib.fp (d12_o0[r][i,...,j]), + lib.fp (d12_o1[r][i,...,j]), 9) + + def test_ham_s2_ovlp (self): + h1, h2 = ham_2q (las, las.mo_coeff, veff_c=None, h2eff_sub=None)[1:] + lbls = ('ham','s2','ovlp') + mats_o0 = op_o0.ham (las, h1, h2, las.ci, nelec_frs, orbsym=orbsym, wfnsym=wfnsym) + fps_o0 = [lib.fp (mat) for mat in mats_o0] + mats_o1 = op_o1.ham (las, h1, h2, las.ci, nelec_frs, orbsym=orbsym, wfnsym=wfnsym) + for lbl, mat, fp in zip (lbls, mats_o1, fps_o0): + with self.subTest(matrix=lbl): + self.assertAlmostEqual (lib.fp (mat), fp, 9) + + def test_rdm12s (self): + d12_o0 = op_o0.roots_make_rdm12s (las, las.ci, nelec_frs, si, orbsym=orbsym, wfnsym=wfnsym) + d12_o1 = op_o1.roots_make_rdm12s (las, las.ci, nelec_frs, si, orbsym=orbsym, wfnsym=wfnsym) + for r in range (2): + for i in range (nstates): + with self.subTest (rank=r+1, root=i): + self.assertAlmostEqual (lib.fp (d12_o0[r][i]), + lib.fp (d12_o1[r][i]), 9) + +if __name__ == "__main__": + print("Full Tests for LASSI matrix elements of 57-space (91-state) manifold") + unittest.main() + diff --git a/debug/mcpdft/lih_ftlda22_sto3g.py b/debug/mcpdft/lih_ftlda22_sto3g.py new file mode 100644 index 00000000..f63d7dd0 --- /dev/null +++ b/debug/mcpdft/lih_ftlda22_sto3g.py @@ -0,0 +1,91 @@ +import numpy as np +from pyscf import gto, scf +from pyscf.lib import logger +from pyscf import mcpdft +from pyscf.mcpdft.otfnal import transfnal +from pyscf.mcpdft.tfnal_derivs import _ftGGA_jT_op, _unpack_sigma_vector + +mol = gto.M (atom='Li 0 0 0; H 1.2 0 0', basis='sto-3g', verbose=logger.DEBUG2, + output='lih_ftlda22_sto3g.log') +mf = scf.RHF (mol).run () +mc = mcpdft.CASSCF (mf, 'ftLDA,VWN3', 2, 2, grids_level=1).run () +mc_grad = mc.nuc_grad_method ().run () + +def num_vrho (rho, Pi, d): + rhop = rho.copy () + rhop[:,0,:] *= 1+d + rhom = rho.copy () + rhom[:,0,:] *= 1-d + denom = 2*d*rho[:,0,0].sum () + eotp = mc.otfnal.eval_ot (rhop, Pi)[0][0] + eotm = mc.otfnal.eval_ot (rhom, Pi)[0][0] + return (eotp-eotm)/denom + +def num_vPi (rho, Pi, d): + Pip = Pi.copy () + Pip[0,:] *= 1+d + Pim = Pi.copy () + Pim[0,:] *= 1-d + denom = 2*d*Pi[0,0] + eotp = mc.otfnal.eval_ot (rho, Pip)[0][0] + eotm = mc.otfnal.eval_ot (rho, Pim)[0][0] + return (eotp-eotm)/denom + +def debug_jT_op (ot, rho, Pi, vxc): + R = ot.get_ratio (Pi, rho/2) + zeta = ot.get_zeta (R[0], fn_deriv=2) + vot = _ftGGA_jT_op (vxc, rho, Pi, R, zeta, _incl_tGGA=True) + vot = _unpack_sigma_vector (vot, deriv1=rho[1:4], deriv2=Pi[1:4]) + vrho, vPi = vot[0][0,0], vot[1][0,0] + return vrho, vPi + +for R in (0.1, 0.95, 1, 1.05, 1.16): + #if R != 1: continue + print ("R =",R) + rho = np.ones ((2,4,1)) / 2 + rho[:,0,:] *= np.random.rand (1)[0] + rho_tot = rho.sum (0) + Pi = np.ones ((4,1)) * R * 0.25 * rho_tot[0,0] * rho_tot[0,0] + rho[:,1:4,:] *= 1-(2*np.random.rand (3))[None,:,None] + Pi[1:4,:] *= 1-(2*np.random.rand (3))[:,None] + eot, vot, fot = mc.otfnal.eval_ot (rho, Pi) + vrho = vot[0][0,0] + vPi = vot[1][0,0] + # subtract tGGA contributions to debug _ftGGA_jT_op specifically + rho_t = mc.otfnal.get_rho_translated (Pi, rho) + xc_grid = mc.otfnal._numint.eval_xc ( + mc.otfnal.otxc, (rho_t[0], rho_t[1]), spin=1, relativity=0, + deriv=1)[:2] + exc = xc_grid[0] * rho_t[:,0,:].sum (0) + vxcr, vxcs = xc_grid[1][:2] + vxc = list (vxcr.T) + list (vxcs.T) + vot_tGGA = transfnal.jT_op (mc.otfnal, vxc, rho, Pi) + vrho_tGGA, vPi_tGGA = vot_tGGA[:2,0] + for p in range (20): + d = 1.0/(2**p) + vrho1 = num_vrho (rho, Pi, d) + denom = vrho1-vrho_tGGA + vrho_err = vrho-vrho1 + vrho_rerr = 0 + if abs (denom)>0 and abs (vrho_err)>1e-8: + vrho_rerr = vrho_err/denom + vPi1 = num_vPi (rho, Pi, d) + denom = vPi1-vPi_tGGA + vPi_err = vPi-vPi1 + vPi_rerr = 0 + if abs (denom)>0 and abs (vPi_err)>1e-8: + vPi_rerr = vPi_err/denom + print (p, vrho_err, vrho_rerr, vPi_err, vPi_rerr, rho_tot[0,0]*vrho_err/(-2*R), + rho_tot[0,0]*rho_tot[0,0]*vPi_err/4) + rho = rho.sum (0).ravel () + Pi = Pi.ravel () + srr = rho[1:4].dot (rho[1:4]) + srP = rho[1:4].dot (Pi[1:4]) + sPP = Pi[1:4].dot (Pi[1:4]) + x0 = np.ravel (vxc[2:]) + x1 = np.zeros_like (x0) + x1[0] = (x0[0] + x0[2] + x0[1]) / 4.0 + x1[1] = (x0[0] - x0[2]) / 2.0 + x1[2] = (x0[0] + x0[2] - x0[1]) / 4.0 + + diff --git a/debug/mcpdft/lih_ftpbe22_sto3g.py b/debug/mcpdft/lih_ftpbe22_sto3g.py new file mode 100644 index 00000000..e30df8ce --- /dev/null +++ b/debug/mcpdft/lih_ftpbe22_sto3g.py @@ -0,0 +1,91 @@ +import numpy as np +from pyscf import gto, scf +from pyscf.lib import logger +from pyscf import mcpdft +from pyscf.mcpdft.otfnal import transfnal +from pyscf.mcpdft.tfnal_derivs import _ftGGA_jT_op, _unpack_sigma_vector + +mol = gto.M (atom='Li 0 0 0; H 1.2 0 0', basis='sto-3g', verbose=logger.DEBUG2, + output='lih_ftpbe22_sto3g.log') +mf = scf.RHF (mol).run () +mc = mcpdft.CASSCF (mf, 'ftPBE', 2, 2, grids_level=1).run () +mc_grad = mc.nuc_grad_method ().run () + +def num_vrho (rho, Pi, d): + rhop = rho.copy () + rhop[:,0,:] *= 1+d + rhom = rho.copy () + rhom[:,0,:] *= 1-d + denom = 2*d*rho[:,0,0].sum () + eotp = mc.otfnal.eval_ot (rhop, Pi)[0][0] + eotm = mc.otfnal.eval_ot (rhom, Pi)[0][0] + return (eotp-eotm)/denom + +def num_vPi (rho, Pi, d): + Pip = Pi.copy () + Pip[0,:] *= 1+d + Pim = Pi.copy () + Pim[0,:] *= 1-d + denom = 2*d*Pi[0,0] + eotp = mc.otfnal.eval_ot (rho, Pip)[0][0] + eotm = mc.otfnal.eval_ot (rho, Pim)[0][0] + return (eotp-eotm)/denom + +def debug_jT_op (ot, rho, Pi, vxc): + R = ot.get_ratio (Pi, rho/2) + zeta = ot.get_zeta (R[0], fn_deriv=2) + vot = _ftGGA_jT_op (vxc, rho, Pi, R, zeta, _incl_tGGA=True) + vot = _unpack_sigma_vector (vot, deriv1=rho[1:4], deriv2=Pi[1:4]) + vrho, vPi = vot[0][0,0], vot[1][0,0] + return vrho, vPi + +for R in (0.1, 0.95, 1, 1.05, 1.16): + #if R != 1: continue + print ("R =",R) + rho = np.ones ((2,4,1)) / 2 + rho[:,0,:] *= np.random.rand (1)[0] + rho_tot = rho.sum (0) + Pi = np.ones ((4,1)) * R * 0.25 * rho_tot[0,0] * rho_tot[0,0] + rho[:,1:4,:] *= 1-(2*np.random.rand (3))[None,:,None] + Pi[1:4,:] *= 1-(2*np.random.rand (3))[:,None] + eot, vot, fot = mc.otfnal.eval_ot (rho, Pi) + vrho = vot[0][0,0] + vPi = vot[1][0,0] + # subtract tGGA contributions to debug _ftGGA_jT_op specifically + rho_t = mc.otfnal.get_rho_translated (Pi, rho) + xc_grid = mc.otfnal._numint.eval_xc ( + mc.otfnal.otxc, (rho_t[0], rho_t[1]), spin=1, relativity=0, + deriv=1)[:2] + exc = xc_grid[0] * rho_t[:,0,:].sum (0) + vxcr, vxcs = xc_grid[1][:2] + vxc = list (vxcr.T) + list (vxcs.T) + vot_tGGA = transfnal.jT_op (mc.otfnal, vxc, rho, Pi) + vrho_tGGA, vPi_tGGA = vot_tGGA[:2,0] + for p in range (20): + d = 1.0/(2**p) + vrho1 = num_vrho (rho, Pi, d) + denom = vrho1-vrho_tGGA + vrho_err = vrho-vrho1 + vrho_rerr = 0 + if abs (denom)>0 and abs (vrho_err)>1e-8: + vrho_rerr = vrho_err/denom + vPi1 = num_vPi (rho, Pi, d) + denom = vPi1-vPi_tGGA + vPi_err = vPi-vPi1 + vPi_rerr = 0 + if abs (denom)>0 and abs (vPi_err)>1e-8: + vPi_rerr = vPi_err/denom + print (p, vrho_err, vrho_rerr, vPi_err, vPi_rerr, rho_tot[0,0]*vrho_err/(-2*R), + rho_tot[0,0]*rho_tot[0,0]*vPi_err/4) + rho = rho.sum (0).ravel () + Pi = Pi.ravel () + srr = rho[1:4].dot (rho[1:4]) + srP = rho[1:4].dot (Pi[1:4]) + sPP = Pi[1:4].dot (Pi[1:4]) + x0 = np.ravel (vxc[2:]) + x1 = np.zeros_like (x0) + x1[0] = (x0[0] + x0[2] + x0[1]) / 4.0 + x1[1] = (x0[0] - x0[2]) / 2.0 + x1[2] = (x0[0] + x0[2] - x0[1]) / 4.0 + + diff --git a/debug/mcpdft/lih_tpbe22_sto3g.py b/debug/mcpdft/lih_tpbe22_sto3g.py new file mode 100644 index 00000000..cac3a1d0 --- /dev/null +++ b/debug/mcpdft/lih_tpbe22_sto3g.py @@ -0,0 +1,13 @@ +import numpy as np +from pyscf import gto, scf +from pyscf.lib import logger +from mrh.my_pyscf import mcpdft + +mol = gto.M (atom='Li 0 0 0; H 1.2 0 0', basis='sto-3g', verbose=logger.DEBUG2, + output='lih_tpbe22_sto3g.log') +mf = scf.RHF (mol).run () +mc = mcpdft.CASSCF (mf, 'tPBE', 2, 2, grids_level=1).run () +mc_grad = mc.nuc_grad_method ().run () + + + diff --git a/debug/mcpdft/proc_v_err.py b/debug/mcpdft/proc_v_err.py new file mode 100644 index 00000000..d001db1e --- /dev/null +++ b/debug/mcpdft/proc_v_err.py @@ -0,0 +1,69 @@ +import numpy as np +import argparse +from subprocess import Popen, PIPE +from shlex import split +from io import StringIO +from itertools import product + +parser = argparse.ArgumentParser (description="tabulate the output of v_err_report for a logfile") +parser.add_argument ("-l", "--logfile", default='lih_ftpbe22_sto3g.log', help='log file (default: %(default)s)') +args = parser.parse_args () + +logfile=str(args.logfile) +fnals=['xc','ot'] +cases_fnals=[["rhoa","rhob","rhoa'","rhob'"], + ["rho","Pi","rho'","Pi'"]] + +for fnal, cases in zip (fnals, cases_fnals): + for case in cases: + p1 = Popen (split ('grep "gradient debug {}:" {}'.format (case, logfile)), + stdout=PIPE) + p2 = Popen (split ('grep "^eval_{}"'.format (fnal, logfile)), + stdin=p1.stdout, stdout=PIPE) + p3 = Popen (split ("awk '{ print $2,$6,$12,$13 }'"), + stdin=p2.stdout, stdout=PIPE) + stdout, stderr = p3.communicate () + rawtab = stdout.decode ("utf-8").replace (")", "").replace ("(","")[:-1] + rawtab = np.loadtxt (StringIO(rawtab)) + if not rawtab.size: continue + rawtab = rawtab.reshape (-1,20,4).transpose (1,0,2) + idx = rawtab[:,0,0].astype (int) + val = rawtab[:,:,1] + err1 = rawtab[:,:,2] + err2 = rawtab[:,:,3] + fmt_str = ' '.join (['{}',]*(val.shape[1]+1)) + #relerr = np.abs (err1/val) + #print (fnal, case, 'grad') + #for ix, row in zip (idx, relerr): + # print (fmt_str.format (2.0**(-ix), *row)) + relerr = np.abs (err2/val) + print (fnal, case, 'grad+hess') + for ix, row in zip (idx, relerr): + print (fmt_str.format (2.0**(-ix), *row)) + relerr = relerr[1:] / relerr[:-1] + print (fnal, case, 'grad+hess conv') + for ix, row in zip (idx, relerr): + print (fmt_str.format (2.0**(-ix), *row)) + relerr = relerr[1:] / relerr[:-1] + for c1, c2 in product (cases, repeat=2): + p1 = Popen (split ('grep "Hessian debug (H.x_{})_{}:" {}'.format (c1, c2, logfile)), + stdout=PIPE) + p2 = Popen (split ('grep "^eval_{}"'.format (fnal, logfile)), + stdin=p1.stdout, stdout=PIPE) + p3 = Popen (split ("awk '{ print $2,$6,$10 }'"), + stdin=p2.stdout, stdout=PIPE) + stdout, stderr = p3.communicate () + rawtab = np.loadtxt (StringIO(stdout.decode ("utf-8")[:-1])) + if not rawtab.size: continue + rawtab = rawtab.reshape (-1,20,3).transpose (1,0,2) + idx = rawtab[:,0,0].astype (int) + val = rawtab[:,:,1] + err = rawtab[:,:,2] + val[(val==0)&(err==0)] = 1 + relerr = np.abs (err/val) + fmt_str = ' '.join (['{}',]*(relerr.shape[1]+1)) + print (fnal, c1, c2) + for ix, row in zip (idx, relerr): + print (fmt_str.format (2.0**(-ix), *row)) + + diff --git a/examples/laspdft/c2h4n4_si_laspdft.py b/examples/laspdft/c2h4n4_si_laspdft.py index ec1dcbd6..944d1695 100755 --- a/examples/laspdft/c2h4n4_si_laspdft.py +++ b/examples/laspdft/c2h4n4_si_laspdft.py @@ -1,7 +1,7 @@ from pyscf import gto, scf, lib, mcscf from mrh.my_pyscf.fci import csf_solver from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF -from mrh.my_pyscf.mcscf import lassi +from mrh.my_pyscf import lassi from mrh.my_pyscf import mcpdft from mrh.my_pyscf.tools.molden import from_lasscf from c2h4n4_struct import structure as struct @@ -22,16 +22,27 @@ guess_mo = las.sort_mo([16,18,22,23,24,26]) mo0 = las.localize_init_guess((list(range (5)), list(range (5,10))), guess_mo) las.kernel(mo0) - + +las = lassi.states.all_single_excitations (las) +las.lasci () + +lsi = lassi.LASSI(las) +lsi.kernel() + # LASSI-PDFT -mc = mcpdft.LASSI(las, 'tPBE', (3, 3), ((2,1),(1,2))) -mc = all_single_excitations(mc) # Level of charge transfer -mc.kernel(las.mo_coeff) # SA-LAS orbitals +mc = mcpdft.LASSI(lsi, 'tPBE', (3, 3), ((2,1),(1,2))) +mc.kernel() + +# CASCI-PDFT in las orbitals +from pyscf import mcpdft +mc_ci = mcpdft.CASCI(mf, 'tPBE', 6, 6) +mc_ci.kernel(las.mo_coeff) # Results print("\n----Results-------\n") -#print("State",' \t', "LASSCF Energy",'\t\t',"LASSI Energy",'\t\t', "LASSI-PDFT Energy") -#[print(sn,'\t',x,'\t', y,'\t', z) for sn, x, y, z in zip(list(range(mc.nroots)), mc.e_mcscf, mc.e_lassi, mc.e_tot)] -print("LASSI state-0 =", mc.e_mcscf[0]) +print("CASCI state-0 =", mc_ci.e_mcscf) +print("LASSI state-0 =", lsi.e_roots[0]) + +print("CASCI-PDFT state-0 =", mc_ci.e_tot) print("LASSI-PDFT state-0 =", mc.e_tot[0]) diff --git a/examples/mcscf_mcpdft_dferi_geomopts/h2co_ftpbe66_631g_grad.py b/examples/mcscf_mcpdft_dferi_geomopts/h2co_ftpbe66_631g_grad.py new file mode 100644 index 00000000..558adf6a --- /dev/null +++ b/examples/mcscf_mcpdft_dferi_geomopts/h2co_ftpbe66_631g_grad.py @@ -0,0 +1,89 @@ +#from rhf import monkeypatch_setup +#monkeypatch_teardown = monkeypatch_setup () +import math +import numpy as np +from scipy import linalg +from pyscf import scf, gto, lib, mcscf, df +from pyscf import mcpdft +from pyscf.df.grad import mcpdft as mcpdft_grad +from mrh.my_pyscf.grad import numeric as numeric_grad + +# Convenience functions to get the internal coordinates for human inspection +def bond_length (carts, i, j): + return linalg.norm (carts[i] - carts[j]) +def bond_angle (carts, i, j, k): + rij = carts[i] - carts[j] + rkj = carts[k] - carts[j] + res = max (min (1.0, np.dot (rij, rkj) / linalg.norm (rij) / linalg.norm (rkj)), -1.0) + return math.acos (res) * 180 / math.pi +def out_of_plane_angle (carts, i, j, k, l): + eji = carts[j] - carts[i] + eki = carts[k] - carts[i] + eli = carts[l] - carts[i] + eji /= linalg.norm (eji) + eki /= linalg.norm (eki) + eli /= linalg.norm (eli) + return -math.asin (np.dot (eji, (np.cross (eki, eli) / math.sin (bond_angle (carts, j, i, k) * math.pi / 180)))) * 180 / math.pi +def h2co_geom_analysis (carts): + print ("rCO = {:.4f} Angstrom".format (bond_length (carts, 1, 0))) + print ("rCH1 = {:.4f} Angstrom".format (bond_length (carts, 2, 0))) + print ("rCH2 = {:.4f} Angstrom".format (bond_length (carts, 3, 0))) + print ("tOCH1 = {:.2f} degrees".format (bond_angle (carts, 1, 0, 2))) + print ("tOCH2 = {:.2f} degrees".format (bond_angle (carts, 1, 0, 3))) + print ("tHCH = {:.2f} degrees".format (bond_angle (carts, 3, 0, 2))) + print ("eta = {:.2f} degrees".format (out_of_plane_angle (carts, 0, 2, 3, 1))) + + +h2co_casscf66_631g_xyz = '''C 0.534004 0.000000 0.000000 +O -0.676110 0.000000 0.000000 +H 1.102430 0.000000 0.920125 +H 1.102430 0.000000 -0.920125''' +mol = gto.M (atom = h2co_casscf66_631g_xyz, basis = '6-31g', symmetry = False, verbose = lib.logger.INFO, output = 'h2co_ftpbe66_631g_grad.log') +mf_conv = scf.RHF (mol).run () +mc_conv = mcpdft.CASSCF (mf_conv, 'ftPBE', 6, 6, grids_level=6) +mc_conv.conv_tol = 1e-10 +mc_conv.kernel () + + +mf = scf.RHF (mol).density_fit (auxbasis = df.aug_etb (mol)).run () +mc_df = mcpdft.CASSCF (mf, 'tPBE', 6, 6, grids_level=6) +mc_df.conv_tol = 1e-10 +mc_df.kernel () + +try: + de_num = np.load ('h2co_ftpbe66_631g_grad_num.npy') + de_conv_num, de_df_num = list (de_num) +except OSError as e: + de_conv_num = numeric_grad.Gradients (mc_conv).kernel () + de_df_num = numeric_grad.Gradients (mc_df).kernel () + np.save ('h2co_ftpbe66_631g_grad_num.npy', np.stack ((de_conv_num, de_df_num), axis=0)) +def my_call (env): + carts = env['mol'].atom_coords () * lib.param.BOHR + h2co_geom_analysis (carts) + +conv_params = { + 'convergence_energy': 1e-6, # Eh + 'convergence_grms': 5.0e-5, # Eh/Bohr + 'convergence_gmax': 7.5e-5, # Eh/Bohr + 'convergence_drms': 1.0e-4, # Angstrom + 'convergence_dmax': 1.5e-4, # Angstrom +} + +de_conv = mc_conv.nuc_grad_method ().kernel () +de_df = mcpdft_grad.Gradients (mc_df).kernel () + +def printable_grad (arr): + arr[np.abs (arr)<1e-10] = 0.0 + line_fmt = ' '.join (('{:13.9f}',)*3) + return '\n'.join ((line_fmt.format (*row) for row in arr)) + +print ("Analytic gradient with conventional ERIs:\n", printable_grad (de_conv)[1:]) +print ("Numeric gradient with conventional ERIs:\n", printable_grad (de_conv_num)[1:]) +print ("Gradient error with conventional ERIs = {:.6e}".format (linalg.norm (de_conv - de_conv_num))) +print ("Analytic gradient with DF ERIs:\n", printable_grad (de_df)[1:]) +print ("Numeric gradient with DF ERIs:\n", printable_grad (de_df_num)[1:]) +print ("Gradient error with DF ERIs = {:.6e}".format (linalg.norm (de_df - de_df_num))) +print ("Gradient disagreements of state 0: Analytic = {:.6e} ; Numeric = {:.6e}".format ( + linalg.norm (de_df - de_conv), linalg.norm (de_df_num - de_conv_num))) + + diff --git a/exploratory/citools/dms.py b/exploratory/citools/dms.py new file mode 100644 index 00000000..1b19a7f5 --- /dev/null +++ b/exploratory/citools/dms.py @@ -0,0 +1,110 @@ +import numpy as np +from scipy import linalg +from pyscf.fci import cistring +from pyscf.fci.addons import _unpack_nelec +from pyscf.fci.direct_spin1 import make_rdm12s +from mrh.my_pyscf.mcpdft._dms import dm2_cumulant +from mrh.my_pyscf.fci.csfstring import CSFTransformer + +def random_ci_vector (norb, nelec, csf=False, csf_off=0): + neleca, nelecb = _unpack_nelec (nelec) + na = cistring.num_strings (norb, neleca) + nb = cistring.num_strings (norb, nelecb) + ci0 = np.random.rand (na,nb) + ci0 = ci0 / linalg.norm (ci0) + if csf: + try: + smult = (neleca - nelecb) + 1 + 2*csf_off + t = CSFTransformer (norb, neleca, nelecb, smult) + ci0 = t.vec_det2csf (ci0) + ci0 = t.vec_csf2det (ci0) + except AssertionError as e: + return np.zeros ((na,nb)) + return ci0 + +def random_1det_rdms (norb, nelec): + neleca, nelecb = _unpack_nelec (nelec) + dm1s = [] + for n in (neleca, nelecb): + kappa = np.random.rand (norb, norb) + kappa -= kappa.T + umat = linalg.expm (kappa) [:,:n] + dm1s.append (umat @ umat.T) + assert (abs (np.trace (dm1s[0])-neleca) < 1e-8), dm1s[0] + assert (abs (np.trace (dm1s[1])-nelecb) < 1e-8), dm1s[1] + dm1 = dm1s[0] + dm1s[1] + dm2 = np.multiply.outer (dm1, dm1) + for d in dm1s: + dm2 -= np.multiply.outer (d, d).transpose (0,3,2,1) + return dm1s, dm2 + +def random_rohf_rdms (norb, nelec): + neleca, nelecb = _unpack_nelec (nelec) + kappa = np.random.rand (norb, norb) + kappa -= kappa.T + umat = linalg.expm (kappa) + umata, umatb = umat[:,:neleca], umat[:,:nelecb] + dm1s = [umata @ umata.T, umatb @ umatb.T] + assert (abs (np.trace (dm1s[0])-neleca) < 1e-8), dm1s[0] + assert (abs (np.trace (dm1s[1])-nelecb) < 1e-8), dm1s[1] + dm1 = dm1s[0] + dm1s[1] + dm2 = np.multiply.outer (dm1, dm1) + for d in dm1s: + dm2 -= np.multiply.outer (d, d).transpose (0,3,2,1) + return dm1s, dm2 + +def daniels_fn (dm1s, dm2, nelec): + neleca, nelecb = _unpack_nelec (nelec) + dm1 = dm1s[0] + dm1s[1] + dm1n = (2-(neleca+nelecb)/2.) * dm1 - np.einsum('pkkq->pq', dm2) + dm1n *= 2./(neleca-nelecb+2) + return dm1n + +if __name__=='__main__': + for norb in range (1,7): + for neleca in range (1,norb+1): + for nelecb in range (neleca+1): + #if norb > 2: break + if neleca==norb and nelecb==norb: continue + ci0 = random_ci_vector (norb, (neleca, nelecb), csf=True) + if linalg.norm (ci0) == 0: + continue + dm1s, dm2s = make_rdm12s (ci0, norb, (neleca, nelecb)) + dm2 = dm2s[0] + dm2s[1] + dm2s[2] + dm2s[1].transpose (2,3,0,1) + sm1_ref = dm1s[0] - dm1s[1] + sm1_test = daniels_fn (dm1s, dm2, (neleca, nelecb)) + try: + sm1_err = linalg.norm (sm1_test-sm1_ref) + if np.abs (sm1_err) < 1e-8: sm1_err = "yes" + except ValueError as e: + sm1_err = str (e) + is_1det = (nelecb == 0) or (neleca == norb) + print ("general", norb, (neleca,nelecb), is_1det, sm1_err) + #dm1 = dm1s[0] + dm1s[1] + #cm2 = dm2_cumulant (dm2, dm1s) + #print ("2D_pq =", 2*dm1) + #print ("l_prrq =", np.einsum ('prrq->pq', cm2)) + #print ("D_pr D_qr =", dm1 @ dm1.T) + #print ("2D_pq - l_prrq - D_pr D_qr =", 2*dm1 - dm1 @ dm1.T - np.einsum ('prrq->pq', cm2)) + #print ("M_pq test =", sm1_test) + #print ("M_pq ref =", sm1_ref) + dm1s, dm2 = random_1det_rdms (norb, (neleca, nelecb)) + sm1_ref = dm1s[0] - dm1s[1] + sm1_test = daniels_fn (dm1s, dm2, (neleca, nelecb)) + try: + sm1_err = linalg.norm (sm1_test-sm1_ref) + if np.abs (sm1_err) < 1e-8: sm1_err = "yes" + except ValueError as e: + sm1_err = str (e) + #print ("1det", norb, (neleca,nelecb), sm1_err) + dm1s, dm2 = random_rohf_rdms (norb, (neleca, nelecb)) + sm1_ref = dm1s[0] - dm1s[1] + sm1_test = daniels_fn (dm1s, dm2, (neleca, nelecb)) + try: + sm1_err = linalg.norm (sm1_test-sm1_ref) + if np.abs (sm1_err) < 1e-8: sm1_err = "yes" + except ValueError as e: + sm1_err = str (e) + #print ("ROHF", norb, (neleca,nelecb), sm1_err) + + diff --git a/lib/patch_pyscf_sys_info.py b/lib/patch_pyscf_sys_info.py new file mode 100644 index 00000000..37bb6d31 --- /dev/null +++ b/lib/patch_pyscf_sys_info.py @@ -0,0 +1,20 @@ +import os +import pyscf.lib.misc + +if not hasattr (pyscf.lib.misc, 'mrh_patched') and hasattr (pyscf.lib.misc, 'format_sys_info'): + from pyscf.lib.misc import format_sys_info as pyscf_format_sys_info + from pyscf.lib.misc import repo_info + + def format_sys_info (): + result = pyscf_format_sys_info () + mrh_info = repo_info (os.path.join (__file__, '..')) + result.append (f'mrh path {mrh_info["path"]}') + if 'git' in mrh_info: + result.append (mrh_info['git']) + return result + + + pyscf.lib.misc.format_sys_info = format_sys_info + pyscf.lib.misc.mrh_patched = True + + diff --git a/my_dmet/pyscf_rhf.py b/my_dmet/pyscf_rhf.py index 477ad1ca..0dd56b2b 100644 --- a/my_dmet/pyscf_rhf.py +++ b/my_dmet/pyscf_rhf.py @@ -115,14 +115,14 @@ def __init__(self, my_mf): self.__dict__.update (my_mf.__dict__) def get_fock (self, h1e=None, s1e=None, vhf=None, dm=None, cycle=1, diis=None, - diis_start_cycle=None, level_shift_factor=None, damp_factor=None): + diis_start_cycle=None, level_shift_factor=None, damp_factor=None, fock_last=None): if vhf is None: vhf = self.get_veff (self.mol, dm) vhf[0] += h1e_s vhf[1] -= h1e_s return super().get_fock (h1e=h1e, s1e=s1e, vhf=vhf, dm=dm, cycle=cycle, diis=diis, diis_start_cycle=diis_start_cycle, level_shift_factor=level_shift_factor, - damp_factor=damp_factor) + damp_factor=damp_factor, fock_last=fock_last) def energy_elec (self, dm=None, h1e=None, vhf=None): if dm is None: dm = self.make_rdm1 () diff --git a/my_pyscf/fci/csf.py b/my_pyscf/fci/csf.py index a6dba6d3..138fb87e 100644 --- a/my_pyscf/fci/csf.py +++ b/my_pyscf/fci/csf.py @@ -92,6 +92,10 @@ def make_hdiag_csf (h1e, eri, norb, nelec, transformer, hdiag_det=None, max_memo csd_offset = npair_csd_offset[ipair] det_addr = transformer.csd_mask[csd_offset:][:nconf*ndet] # mem safety + # Issue #54: PySCF wants "max_memory" on entrance to FCI to be "remaining memory". However, + # the first few lines of this function consume some memory, so that's difficult to + # implement consistently. "max_memory" here is currently the config parameter for the whole + # calculation. mem_remaining = max_memory - lib.current_memory ()[0] safety_factor = 1.2 nfloats = nconf*ndet*ndet + det_addr.size @@ -284,6 +288,9 @@ def pspace (fci, h1e, eri, norb, nelec, transformer, hdiag_det=None, hdiag_csf=N safety_factor = 1.2 nfloats_h0 = (npsp_det+npsp)**2 mem_h0 = safety_factor * nfloats_h0 * np.dtype (float).itemsize / 1e6 + # Issue #54: PySCF wants "max_memory" on entrance to FCI to be "remaining memory". However, + # the earlier lines of this function consume some memory, so that's difficult to implement + # consistently. "max_memory" here is currently the config parameter for the whole calculation. mem_remaining = max_memory - lib.current_memory ()[0] memstr = ("pspace_size of {} CSFs -> {} determinants requires {} MB, cf {} MB " "remaining memory").format (npsp, npsp_det, mem_h0, mem_remaining) diff --git a/my_pyscf/grad/mspdft_nacs.py b/my_pyscf/grad/mspdft_nacs.py index acbde816..1721972a 100644 --- a/my_pyscf/grad/mspdft_nacs.py +++ b/my_pyscf/grad/mspdft_nacs.py @@ -1,183 +1,3 @@ -import numpy as np -from pyscf import lib -from pyscf.lib import logger -from pyscf.fci import direct_spin1 -from pyscf.mcscf import newton_casscf -from pyscf.grad import casscf as casscf_grad -from pyscf.grad import sacasscf as sacasscf_grad -from pyscf.grad import mspdft as mspdft_grad -from mrh.my_pyscf.grad import sacasscf_nacs -from functools import reduce - -_unpack_state = mspdft_grad._unpack_state -_nac_csf = sacasscf_nacs._nac_csf - -def nac_model (mc_grad, mo_coeff=None, ci=None, si_bra=None, si_ket=None, - mf_grad=None, atmlst=None): - '''Compute the "model-state contribution" to the MS-PDFT NAC''' - mc = mc_grad.base - mol = mc.mol - ci_bra = np.tensordot (si_bra, ci, axes=1) - ci_ket = np.tensordot (si_ket, ci, axes=1) - ncore, ncas, nelecas = mc.ncore, mc.ncas, mc.nelecas - castm1 = direct_spin1.trans_rdm1 (ci_bra, ci_ket, ncas, nelecas) - # if PySCF commentary is to be trusted, trans_rdm1[p,q] is - # . I want . - castm1 = castm1.conj ().T - castm1 - mo_cas = mo_coeff[:,ncore:][:,:ncas] - tm1 = reduce (np.dot, (mo_cas, castm1, mo_cas.conj ().T)) - return _nac_csf (mol, mf_grad, tm1, atmlst) - -class NonAdiabaticCouplings (mspdft_grad.Gradients): - '''MS-PDFT non-adiabatic couplings (NACs) between states - - kwargs/attributes: - - state : tuple of length 2 - The NACs returned are . - In other words, state = (ket, bra). - mult_ediff : logical - If True, returns NACs multiplied by the energy difference. - Useful near conical intersections to avoid numerical problems. - use_etfs : logical - If True, use the ``electron translation factors'' of Fatehi and - Subotnik [JPCL 3, 2039 (2012)], which guarantee conservation of - total electron + nuclear momentum when the nuclei are moving - (i.e., in non-adiabatic molecular dynamics). This corresponds - to omitting the ``model state contribution''. - ''' - - def __init__(self, mc, state=None, mult_ediff=False, use_etfs=False): - self.mult_ediff = mult_ediff - self.use_etfs = use_etfs - self.state = state - mspdft_grad.Gradients.__init__(self, mc) - - def get_wfn_response (self, si_bra=None, si_ket=None, state=None, si=None, - verbose=None, **kwargs): - g_all = mspdft_grad.Gradients.get_wfn_response ( - self, si_bra=si_bra, si_ket=si_ket, state=state, si=si, - verbose=verbose, **kwargs - ) - g_orb, g_ci, g_is = self.unpack_uniq_var (g_all) - if state is None: state = self.state - ket, bra = _unpack_state (state) - if si is None: si = self.base.si - if si_bra is None: si_bra = si[:,bra] - if si_ket is None: si_ket = si[:,ket] - nroots = self.nroots - log = lib.logger.new_logger (self, verbose) - g_model = np.multiply.outer (si_bra.conj (), si_ket) - g_model -= g_model.T - g_model *= self.base.e_states[bra]-self.base.e_states[ket] - g_model = g_model[np.tril_indices (nroots, k=-1)] - log.debug ("NACs g_is additional component:\n{}".format (g_model)) - return self.pack_uniq_var (g_orb, g_ci, g_is+g_model) - - def get_ham_response (self, **kwargs): - nac = mspdft_grad.Gradients.get_ham_response (self, **kwargs) - use_etfs = kwargs.get ('use_etfs', self.use_etfs) - if not use_etfs: - verbose = kwargs.get ('verbose', self.verbose) - log = lib.logger.new_logger (self, verbose) - nac_model = self.nac_model (**kwargs) - log.info ('NACs model-state contribution:\n{}'.format (nac_model)) - nac += nac_model - return nac - - def nac_model (self, mo_coeff=None, ci=None, si=None, si_bra=None, - si_ket=None, state=None, mf_grad=None, atmlst=None, - **kwargs): - if state is None: state = self.state - ket, bra = _unpack_state (state) - if mo_coeff is None: mo_coeff = self.base.mo_coeff - if ci is None: ci = self.base.ci - if si is None: si = self.base.si - if si_bra is None: si_bra = si[:,bra] - if si_ket is None: si_ket = si[:,ket] - if mf_grad is None: mf_grad = self.base._scf.nuc_grad_method () - if atmlst is None: atmlst = self.atmlst - nac = nac_model (self, mo_coeff=mo_coeff, ci=ci, si_bra=si_bra, - si_ket=si_ket, mf_grad=mf_grad, atmlst=atmlst) - e_bra = self.base.e_states[bra] - e_ket = self.base.e_states[ket] - nac *= e_bra - e_ket - return nac - - def kernel (self, *args, **kwargs): - mult_ediff = kwargs.get ('mult_ediff', self.mult_ediff) - state = kwargs.get ('state', self.state) - nac = mspdft_grad.Gradients.kernel (self, *args, **kwargs) - if not mult_ediff: - ket, bra = _unpack_state (state) - e_bra = self.base.e_states[bra] - e_ket = self.base.e_states[ket] - nac /= e_bra - e_ket - return nac - -if __name__=='__main__': - from pyscf import gto, scf, mcscf, mcpdft - from mrh.my_pyscf.dft.openmolcas_grids import quasi_ultrafine - from scipy import linalg - mol = gto.M (atom = 'Li 0 0 0; H 0 0 1.5', basis='sto-3g', - output='mspdft_nacs.log', verbose=lib.logger.INFO) - mf = scf.RHF (mol).run () - mc = mcpdft.CASSCF (mf, 'ftLDA,VWN3', 2, 2, grids_attr=quasi_ultrafine) - mc = mc.fix_spin_(ss=0).multi_state ([0.5,0.5], 'cms').run (conv_tol=1e-10) - #openmolcas_energies = np.array ([-7.85629118, -7.72175252]) - print ("energies:",mc.e_states) - #print ("disagreement w openmolcas:", np.around (mc.e_states-openmolcas_energies, 8)) - mc_nacs = NonAdiabaticCouplings (mc) - print ("no csf contr") - print ("antisym") - nac_01 = mc_nacs.kernel (state=(0,1), use_etfs=True) - nac_10 = mc_nacs.kernel (state=(1,0), use_etfs=True) - print (nac_01) - print (nac_10) - print ("checking antisym:",linalg.norm(nac_01+nac_10)) - print ("sym") - nac_01_mult = mc_nacs.kernel (state=(0,1), use_etfs=True, mult_ediff=True) - nac_10_mult = mc_nacs.kernel (state=(1,0), use_etfs=True, mult_ediff=True) - print (nac_01_mult) - print (nac_10_mult) - print ("checking sym:",linalg.norm(nac_01_mult-nac_10_mult)) - - - print ("incl csf contr") - print ("antisym") - nac_01 = mc_nacs.kernel (state=(0,1), use_etfs=False) - nac_10 = mc_nacs.kernel (state=(1,0), use_etfs=False) - print (nac_01) - print ("checking antisym:",linalg.norm(nac_01+nac_10)) - print ("sym") - nac_01_mult = mc_nacs.kernel (state=(0,1), use_etfs=False, mult_ediff=True) - nac_10_mult = mc_nacs.kernel (state=(1,0), use_etfs=False, mult_ediff=True) - print (nac_01_mult) - print ("checking sym:",linalg.norm(nac_01_mult-nac_10_mult)) - - print ("Check gradients") - mc_grad = mc.nuc_grad_method () - de_0 = mc_grad.kernel (state=0) - print (de_0) - de_1 = mc_grad.kernel (state=1) - print (de_1) - - #from mrh.my_pyscf.tools.molcas2pyscf import * - #mol = get_mol_from_h5 ('LiH_sa2casscf22_sto3g.rasscf.h5', - # output='sacasscf_nacs_fromh5.log', - # verbose=lib.logger.INFO) - #mo = get_mo_from_h5 (mol, 'LiH_sa2casscf22_sto3g.rasscf.h5') - #nac_etfs_ref = np.array ([9.14840490109073E-02, -9.14840490109074E-02]) - #nac_ref = np.array ([1.83701578323929E-01, -6.91459741744125E-02]) - #mf = scf.RHF (mol).run () - #mc = mcscf.CASSCF (mol, 2, 2).fix_spin_(ss=0).state_average ([0.5,0.5]) - #mc.run (mo, natorb=True, conv_tol=1e-10) - #mc_nacs = NonAdiabaticCouplings (mc) - #nac = mc_nacs.kernel (state=(0,1)) - #print (nac) - #print (nac_ref) - #nac_etfs = mc_nacs.kernel (state=(0,1), use_etfs=True) - #print (nac_etfs) - #print (nac_etfs_ref) - - +from pyscf.nac.mspdft import * +from mrh.util.io import mcpdft_removal_warn +mcpdft_removal_warn () diff --git a/my_pyscf/grad/sacasscf_nacs.py b/my_pyscf/grad/sacasscf_nacs.py.old similarity index 100% rename from my_pyscf/grad/sacasscf_nacs.py rename to my_pyscf/grad/sacasscf_nacs.py.old diff --git a/my_pyscf/lassi/__init__.py b/my_pyscf/lassi/__init__.py index 55550acf..c032291d 100644 --- a/my_pyscf/lassi/__init__.py +++ b/my_pyscf/lassi/__init__.py @@ -1,3 +1,4 @@ from mrh.my_pyscf.lassi.lassi import LASSI from mrh.my_pyscf.lassi.lassis import LASSIS +from mrh.my_pyscf.lassi.lassirq import LASSIrq, LASSIrqCT diff --git a/my_pyscf/lassi/excitations.py b/my_pyscf/lassi/excitations.py index 6275d942..4e662776 100644 --- a/my_pyscf/lassi/excitations.py +++ b/my_pyscf/lassi/excitations.py @@ -15,18 +15,8 @@ LOWEST_REFOVLP_EIGVAL_THRESH = getattr (__config__, 'lassi_excitations_refovlp_eigval_thresh', 1e-9) IMAG_SHIFT = getattr (__config__, 'lassi_excitations_imag_shift', 1e-6) - -def only_ground_states (ci0): - '''For a list of sequences of CI vectors in the same Hilbert space, - generate a list in which all but the first element of each sequence - is discarded.''' - # TODO: examine whether this should be considered a limitation - ci1 = [] - for c in ci0: - c = np.asarray (c) - if c.ndim==3: c = c[0] - ci1.append (c) - return ci1 +MAX_CYCLE_E0 = getattr (__config__, 'lassi_excitations_max_cycle_e0', 1) +CONV_TOL_E0 = getattr (__config__, 'lassi_excitations_conv_tol_e0', 1e-8) def lowest_refovlp_eigval (ham_pq, ovlp_thresh=LOWEST_REFOVLP_EIGVAL_THRESH): ''' Return the lowest eigenvalue of the matrix ham_pq, whose corresponding @@ -46,6 +36,9 @@ def sort_ci0 (obj, ham_pq, ci0): (h_pp + h_pq (e0 - e_q)^-1 h_qp) |p> = e0|p> + NOTE: in the current LASSIS method, this isn't how the converged vectors should be + sorted: this arrangement only for the initial guess. + Args: ham_pq: ndarray of shape (p+q,p+q) Hamiltonian matrix in model-space basis. In the p-space, ENVs ascend in @@ -94,6 +87,7 @@ def project_1p (ip): e_p_arr = e_p.reshape (*lroots[::-1]).T h_pp = h_pp.reshape (*(list(lroots[::-1])*2)) h_pq = ham_pq[:p,p:].reshape (*(list(lroots[::-1])+[q,])) + ci1 = [c.copy () for c in ci0] for ifrag in range (nfrag): if lroots[ifrag]<2: continue e_p_slice = e_p_arr @@ -103,7 +97,7 @@ def project_1p (ip): e_p_slice = e_p_slice[:,addr[jfrag]] sort_idx = np.argsort (e_p_slice) assert (sort_idx[0] == addr[ifrag]) - ci0[ifrag] = np.stack ([ci0[ifrag][i] for i in sort_idx], axis=0) + ci1[ifrag] = np.stack ([ci1[ifrag][i] for i in sort_idx], axis=0) dimorder = list(range(h_pp.ndim)) dimorder.insert (0, dimorder.pop (nfrag-(1+ifrag))) dimorder.insert (1, dimorder.pop (2*nfrag-(1+ifrag))) @@ -120,7 +114,7 @@ def project_1p (ip): h_pp = h_pp.reshape (p, p) h_pq = h_pq.reshape (p, q) ham_pq = np.block ([[h_pp, h_pq],[h_pq.conj ().T, h_qq]]) - return ci0, e0_p, ham_pq + return ci1, e0_p, ham_pq class _vrvloop_env (object): def __init__(self, fciobj, vrvsolvers, e_q, si_q): @@ -219,7 +213,7 @@ def get_excited_h (self, h0, h1, h2): norb_ref = [norb_excited,] + [n for ifrag, n in enumerate (self.norb_ref) if not (ifrag in self.excited_frags)] h1eff, h0eff = self.project_hfrag (h1, h2, self.ci_ref, norb_ref, self.nelec_ref, - ecore=h0, dm1s=dm1s, dm2=dm2) + ecore=h0, dm1s=dm1s, dm2=dm2)[:2] h0, h1 = h0eff[0], h1eff[0] h2 = h2[:norb_excited][:,:norb_excited][:,:,:norb_excited][:,:,:,:norb_excited] return h0, h1, h2 @@ -256,9 +250,9 @@ def set_excited_fragment_(self, ifrag, nelec, smult, weights=None): self.excited_frags = [self.excited_frags[i] for i in idx] self.fcisolvers = [self.fcisolvers[i] for i in idx] - def kernel (self, h1, h2, ecore=0, - conv_tol_grad=1e-4, conv_tol_self=1e-10, max_cycle_macro=50, - serialfrag=False, **kwargs): + def kernel (self, h1, h2, ecore=0, ci0=None, + conv_tol_grad=1e-4, conv_tol_self=1e-6, max_cycle_macro=50, + serialfrag=False, _add_vrv_energy=False, **kwargs): h0, h1, h2 = self.get_excited_h (ecore, h1, h2) norb_f = np.asarray ([self.norb_ref[ifrag] for ifrag in self.excited_frags]) nelec_f = np.asarray ([self.nelec_ref[ifrag] for ifrag in self.excited_frags]) @@ -267,74 +261,76 @@ def kernel (self, h1, h2, ecore=0, idx = self.get_excited_orb_idx () orbsym = [orbsym[iorb] for iorb in range (norb_tot) if idx[iorb]] # TODO: point group symmetry; I probably also have to do something to wfnsym - ci0, vrvsolvers, e_q, si_q = self.prepare_vrvsolvers_(h0, h1, h2) + ci0, vrvsolvers, e_q, si_q = self.prepare_vrvsolvers_(h0, h1, h2, ci0=ci0) with _vrvloop_env (self, vrvsolvers, e_q, si_q): converged, energy_elec, ci1_active = ProductStateFCISolver.kernel ( self, h1, h2, norb_f, nelec_f, ecore=h0, ci0=ci0, orbsym=orbsym, conv_tol_grad=conv_tol_grad, conv_tol_self=conv_tol_self, max_cycle_macro=max_cycle_macro, serialfrag=serialfrag, **kwargs ) + if _add_vrv_energy: # for a sanity check in unittests only + energy_elec += self._energy_vrv (h1, h2, ci1_active) ci1 = [c for c in self.ci_ref] for ifrag, c in zip (self.excited_frags, ci1_active): ci1[ifrag] = np.asarray (c) return converged, energy_elec, ci1 def project_hfrag (self, h1, h2, ci, norb_f, nelec_f, ecore=0, dm1s=None, dm2=None, **kwargs): - h1eff, h0eff = ProductStateFCISolver.project_hfrag ( + h1eff, h0eff, ci = ProductStateFCISolver.project_hfrag ( self, h1, h2, ci, norb_f, nelec_f, ecore=ecore, dm1s=dm1s, dm2=dm2, **kwargs ) + nj = np.cumsum (norb_f) + ni = nj - norb_f # Project the part coupling the p and q rootspaces + ci1 = [c for c in ci] if len (self._e_q) and not self._deactivate_vrv: ci0 = [np.asarray (c) for c in ci] - ham_pq = self.get_ham_pq (ecore, h1, h2, ci0) - # TODO: optimize this so that we only need one op_ham_pq_ref call each time we land here, - # and not get_ham_pq - ci0, e0 = self.sort_ci0 (ham_pq, ci0) hci_f_pabq = self.op_ham_pq_ref (h1, h2, ci0) - for ifrag, (c, hci_pabq, solver) in enumerate (zip (ci0, hci_f_pabq, self.fcisolvers)): + zipper = zip (hci_f_pabq, self.fcisolvers, ci0, norb_f, nelec_f, h1eff, h0eff, ni, nj) + for ifrag, (hci_pabq, solver, c, norb, nelec, h1e, h0e, i, j) in enumerate (zipper): + h2e = h2[i:j,i:j,i:j,i:j] + ne = self._get_nelec (solver, nelec) solver.v_qpab = np.tensordot (self._si_q, hci_pabq, axes=((0),(-1))) - # The slight disagreement here is causing massive convergence problems! - ci[ifrag] = c - return h1eff, h0eff - - def energy_elec (self, h1, h2, ci, norb_f, nelec_f, ecore=0, **kwargs): - energy_tot = ProductStateFCISolver.energy_elec ( - self, h1, h2, ci, norb_f, nelec_f, ecore=ecore, **kwargs - ) - # Also compute the vrv perturbation energy - if len (self._e_q) and not self._deactivate_vrv: - ci0 = [] - # extract this from the putatively converged solver cycles - # if you attempt to recalculate it, then it has to be reconverged i guess - denom_q = 0 - for c, solver in zip (ci, self.fcisolvers): - t = solver.transformer - c = np.asarray (c).reshape (-1, t.ndeta, t.ndetb) - ci0.append (c) - denom_q += solver.denom_q - denom_q /= len (ci) - lroots = get_lroots (ci0) - p = np.prod (lroots) - ham_pq = self.get_ham_pq (ecore, h1, h2, ci0) - idx = np.ones (len (ham_pq), dtype=bool) - idx[1:p] = False - e0, si = linalg.eigh (ham_pq[idx,:][:,idx]) - h_pp = ham_pq[:p,:p] - h_pq = ham_pq[:p,p:] - h_qq = ham_pq[p:,p:] - h_qq = self._si_q.conj ().T @ h_qq @ self._si_q - h_pq = np.dot (h_pq, self._si_q) - ham_pq = ham_pq[idx,:][:,idx] - idx = np.abs (denom_q) > 1e-16 - e_p = np.diag (np.dot (h_pq[:,idx].conj () / denom_q[None,idx], h_pq[:,idx].T)) - e_p = e_p.reshape (*lroots[::-1]).T - for solver in self.fcisolvers: - if hasattr (getattr (solver, 'weights', None), '__len__'): - e_p = np.dot (solver.weights, e_p) - else: - e_p = e_p[0] - energy_tot += e_p - return energy_tot + e0, ci1[ifrag] = solver.sort_ci (h0e, h1e, h2e, norb, ne, c) + solver.denom_q = e0 - solver.e_q + # The reason I do the above two lines here and not in the fragment-solver kernel is + # that between this function and the start of the fragment-solver kernel, the + # wrapper needs to compute the current gradient for the whole system. The fragment + # solvers need to already know the correct e0 for that gradient calculation. + return h1eff, h0eff, ci1 + + def _energy_vrv (self, h1, h2, ci): + # The only purpose this serves is for a sanity test in test_excitations.py + # It can be reattached to energy_elec if I can figure out how to define a whole-system + # energy for this step in LASSIS that makes sense. + # If so, I should rewrite it so that it can be called outside of this class's kernel + # function + if (len (self._e_q) == 0) or self._deactivate_vrv: return 0 + ci0 = [] + denom_q = 0 + for c, solver in zip (ci, self.fcisolvers): + t = solver.transformer + c = np.asarray (c).reshape (-1, t.ndeta, t.ndetb) + ci0.append (c) + denom_q += solver.denom_q + denom_q /= len (ci) + lroots = get_lroots (ci0) + p = np.prod (lroots) + ham_pq = self.get_ham_pq (0, h1, h2, ci0) + h_pp = ham_pq[:p,:p] + h_pq = ham_pq[:p,p:] + h_qq = ham_pq[p:,p:] + h_qq = self._si_q.conj ().T @ h_qq @ self._si_q + h_pq = np.dot (h_pq, self._si_q) + idx = np.abs (denom_q) > 1e-16 + e_p = np.diag (np.dot (h_pq[:,idx].conj () / denom_q[None,idx], h_pq[:,idx].T)) + e_p = e_p.reshape (*lroots[::-1]).T + for solver in self.fcisolvers: + if hasattr (getattr (solver, 'weights', None), '__len__'): + e_p = np.dot (solver.weights, e_p) + else: + e_p = e_p[0] + return e_p def get_ham_pq (self, h0, h1, h2, ci_p): '''Build the model-space Hamiltonian matrix for the current state of the P-space. @@ -427,18 +423,20 @@ def op_ham_pq_ref (self, h1, h2, ci): def sort_ci0 (self, ham_pq, ci0): return sort_ci0 (self, ham_pq, ci0)[:2] - def prepare_vrvsolvers_(self, h0, h1, h2): + def prepare_vrvsolvers_(self, h0, h1, h2, ci0=None): + do_sort_ci0 = (ci0 is None) norb_f = np.asarray ([self.norb_ref[ifrag] for ifrag in self.excited_frags]) nelec_f = np.asarray ([self.nelec_ref[ifrag] for ifrag in self.excited_frags]) - ci0 = self.get_init_guess (None, norb_f, nelec_f, h1, h2) + ci0 = self.get_init_guess (ci0, norb_f, nelec_f, h1, h2) ham_pq = self.get_ham_pq (h0, h1, h2, ci0) p = np.prod (get_lroots (ci0)) h_qq = ham_pq[p:,p:] e_q, si_q = linalg.eigh (h_qq) - ci0, e0 = self.sort_ci0 (ham_pq, ci0) + ci0_sorted, e0 = self.sort_ci0 (ham_pq, ci0) + if do_sort_ci0: ci0 = ci0_sorted vrvsolvers = [] for ix, solver in enumerate (self.fcisolvers): - vrvsolvers.append (vrv_fcisolver (solver, e0, e_q, None, max_cycle_e0=1, + vrvsolvers.append (vrv_fcisolver (solver, e0, e_q, None, max_cycle_e0=MAX_CYCLE_E0, crash_locmin=self.crash_locmin)) return ci0, vrvsolvers, e_q, si_q @@ -493,8 +491,8 @@ class VRVDressedFCISolver (object): ''' _keys = {'contract_vrv', 'base', 'v_qpab', 'denom_q', 'e_q', 'max_cycle_e0', 'conv_tol_e0', 'charge', 'crash_locmin', 'imag_shift'} - def __init__(self, fcibase, my_vrv, my_eq, my_e0, max_cycle_e0=100, conv_tol_e0=1e-8, - crash_locmin=False): + def __init__(self, fcibase, my_vrv, my_eq, my_e0, max_cycle_e0=MAX_CYCLE_E0, + conv_tol_e0=CONV_TOL_E0, crash_locmin=False): self.base = copy.copy (fcibase) if isinstance (fcibase, StateAverageFCISolver): self._undressed_class = fcibase._base_class @@ -510,12 +508,14 @@ def __init__(self, fcibase, my_vrv, my_eq, my_e0, max_cycle_e0=100, conv_tol_e0= self.crash_locmin = crash_locmin self.davidson_only = self.base.davidson_only = True # TODO: Relaxing this ^ requires accounting for pspace, precond, and/or hdiag - def contract_2e(self, eri, fcivec, norb, nelec, link_index=None, **kwargs): + def contract_2e(self, eri, fcivec, norb, nelec, link_index=None, v_qpab=None, denom_q=None, + **kwargs): ci0 = self.undressed_contract_2e (eri, fcivec, norb, nelec, link_index, **kwargs) - ci0 += self.contract_vrv (fcivec) + ci0 += self.contract_vrv (fcivec, v_qpab=v_qpab, denom_q=denom_q) return ci0 - def contract_vrv (self, ket): - v_qpab, denom_q = self.v_qpab, self.denom_q + def contract_vrv (self, ket, v_qpab=None, denom_q=None): + if v_qpab is None: v_qpab = self.v_qpab + if denom_q is None: denom_q = self.denom_q if v_qpab is None: return np.zeros_like (ket) ket_shape = ket.shape idx = np.abs (denom_q) > 1e-16 @@ -547,8 +547,9 @@ def test_locmin (self, e0, ci, norb, nelec, h0e, h1e, h2e, warntag='Apparent loc e_pq = np.append ([e_p,], e_q) h_diagmin = np.amin (e_pq) if e0-h_diagmin > 1e-8: - log.warn ("%s in VRVSolver: min (hdiag) = %.6f < e0 = %.6f", - warntag, np.amin (e_pq), e0) + if self.verbose >= lib.logger.DEBUG: + log.warn ("%s in VRVSolver: min (hdiag) = %.6f < e0 = %.6f", + warntag, np.amin (e_pq), e0) log.debug ('e_p = %.6f ; vrv = %.6f', e_p, vrv) log.debug ('e_q = {}'.format (e_q)) log.debug ('valid de_pq = {}'.format (de_pq[idx])) @@ -588,6 +589,22 @@ def solve_e0 (self, h0e, h1e, h2e, norb, nelec, ket): log.debug2 ('e_pq = {}'.format (e_pq)) e0 = lowest_refovlp_eigval (ham_pq) return e0 + def sort_ci (self, h0e, h1e, h2e, norb, nelec, ci): + if self.nroots==1: ci = [ci] + e0 = [self.solve_e0 (h0e, h1e, h2e, norb, nelec, ket) for ket in ci] + idx = np.argsort (e0) + e0 = [e0[ix] for ix in idx] + ci = [ci[ix] for ix in idx] + den = e0[0] - self.e_q + h2eff = self.absorb_h1e (h1e, h2e, norb, nelec, 0.5) + e = [np.dot (ket.ravel(), self.contract_2e (h2eff, ket, norb, nelec, denom_q=den).ravel()) + for ket in ci] + if self.nroots > 1: + idx = np.argsort (e[1:]) + ci = [ci[0]] + [ci[1:][ix] for ix in idx] + else: + ci = ci[0] + return e0[0], ci def kernel (self, h1e, h2e, norb, nelec, ecore=0, ci0=None, orbsym=None, **kwargs): log = lib.logger.new_logger (self, self.verbose) max_cycle_e0 = self.max_cycle_e0 @@ -599,45 +616,40 @@ def kernel (self, h1e, h2e, norb, nelec, ecore=0, ci0=None, orbsym=None, **kwarg ci1 = ci0 self.denom_q = e0 - self.e_q log.debug ("Self-energy singularities in VRVSolver: {}".format (self.e_q)) + log.debug ("e0 = %.8g", e0) log.debug ("Denominators in VRVSolver: {}".format (self.denom_q)) self.test_locmin (e0, ci1, norb, nelec, ecore, h1e, h2e, warntag='Saddle-point initial guess') - warn_swap = False # annoying loud warning not necessary - #print (lib.fp (ci0), self.denom_q) + h2eff = self.absorb_h1e (h1e, h2e, norb, nelec, 0.5) for it in range (max_cycle_e0): e, ci1 = self.undressed_kernel ( h1e, h2e, norb, nelec, ecore=ecore, ci0=ci1, orbsym=orbsym, **kwargs ) - e0_last = e0 + # Subtract the vrv energy so that agreement between different fragments can + # be checked in the impure-state case if isinstance (e, (list,tuple,np.ndarray)): - delta_e0 = np.array (e) - e0 - if warn_swap and np.argmin (np.abs (delta_e0)) != 0: - log.warn ("Possible rootswapping in H(E)|Psi> = E|Psi> fixed-point iteration") - warn_swap = False - ket, e0 = ci1[0], e[0] + for i in range (len (e)): + hci = self.undressed_contract_2e (h2eff, ci1[i], norb, nelec) + e[i] = ecore + np.dot (ci1[i].ravel (), hci.ravel ()) else: - ket, e0 = ci1, e + hci = self.undressed_contract_2e (h2eff, ci1, norb, nelec) + e = ecore + np.dot (ci1.ravel (), hci.ravel ()) + e0_last = e0 e0 = self.solve_e0 (ecore, h1e, h2e, norb, nelec, ket) self.denom_q = e0 - self.e_q + log.debug ("e0 = %.8g", e0) log.debug ("Denominators in VRVSolver: {}".format (self.denom_q)) - #hket = self.contract_2e (self.absorb_h1e (h1e, h2e, norb, nelec, 0.5), ket, norb, nelec) - #brahket = np.dot (ket.ravel (), hket.ravel ()) - #e0_test = ecore + brahket - #g_test = hket - ket*brahket - #print (e, e0, e0_test, linalg.norm (g_test)) if abs(e0-e0_last)1: @@ -220,7 +244,10 @@ def cisolve (sm, nroots): ifrag, nelec, norb, smult-2) smults1_i.extend ([smult-2,]*(smult-2)) spins1_i.extend (list (range (smult-3, -(smult-3)-1, -2))) - ci1_i.extend (cisolve (smult-2, ndn0[ifrag])) + ci0 = lsi.ci_spin_flips.get ((ifrag,'d'), None) + ci1_i_down = cisolve (smult-2, ndn0[ifrag], ci0) + lsi.ci_spin_flips[(ifrag,'d')] = ci1_i_down[0] + ci1_i.extend (ci1_i_down) min_npair = max (0, nelec-norb) max_smult = (nelec - 2*min_npair) + 1 if smult < max_smult: # spin-raised @@ -228,20 +255,26 @@ def cisolve (sm, nroots): ifrag, nelec, norb, smult+2) smults1_i.extend ([smult+2,]*(smult+2)) spins1_i.extend (list (range (smult+1, -(smult+1)-1, -2))) - ci1_i.extend (cisolve (smult+2, nup0[ifrag])) + ci0 = lsi.ci_spin_flips.get ((ifrag,'u'), None) + ci1_i_up = cisolve (smult+2, nup0[ifrag], ci0) + lsi.ci_spin_flips[(ifrag,'u')] = ci1_i_up[0] + ci1_i.extend (ci1_i_up) smults1.append (smults1_i) spins1.append (spins1_i) ci1.append (ci1_i) i = j - spin_halfexcs = [SpinHalfexcitations (c,m,s) for c, m, s in zip (ci1, spins1, smults1)] - return spin_halfexcs + spin_flips = [SpinFlips (c,m,s) for c, m, s in zip (ci1, spins1, smults1)] + return spin_flips -def _spin_halfexcitation_products (spaces, spin_halfexcs, nroots_ref=1, frozen_frags=None): - if spin_halfexcs is None or len (spin_halfexcs)==0: return spaces +def _spin_flip_products (spaces, spin_flips, nroots_ref=1, frozen_frags=None): + # NOTE: this actually only uses the -first- rootspace in las, so it can be done before + # the initial spin shuffle + '''Combine spin-flip excitations in all symmetrically permissible ways''' + if spin_flips is None or len (spin_flips)==0: return spaces spaces_ref = spaces[:nroots_ref] - spins3 = [she.spins for she in spin_halfexcs] - smults3 = [she.smults for she in spin_halfexcs] - ci3 = [she.ci for she in spin_halfexcs] + spins3 = [she.spins for she in spin_flips] + smults3 = [she.smults for she in spin_flips] + ci3 = [she.ci for she in spin_flips] nelec0 = spaces[0].nelec smults0 = spaces[0].smults nfrags = spaces[0].nfrag @@ -252,6 +285,9 @@ def _spin_halfexcitation_products (spaces, spin_halfexcs, nroots_ref=1, frozen_f new_spaces = [] m3, s3, c3 = spins3[ifrag], smults3[ifrag], ci3[ifrag] for space in spaces: + # I want to inject the spin-flip into all distinct references, + # but if two references differ only in ifrag then this would + # generate duplicates. The two lines below filter this case. if space.nelec[ifrag] != nelec0[ifrag]: continue if space.smults[ifrag] != smults0[ifrag]: continue for m3i, s3i, c3i in zip (m3, s3, c3): @@ -268,11 +304,12 @@ def _spin_halfexcitation_products (spaces, spin_halfexcs, nroots_ref=1, frozen_f spaces = [space for space in spaces if not ((space in seen) or seen.add (space))] return spaces -def spin_halfexcitation_products (las2, spin_halfexcs, nroots_ref=1): +def spin_flip_products (las2, spin_flips, nroots_ref=1): + '''Inject spin-flips into las2 in all possible ways''' log = logger.new_logger (las2, las2.verbose) spaces = [SingleLASRootspace (las2, m, s, c, las2.weights[ix], ci=[c[ix] for c in las2.ci]) for ix, (c, m, s, w) in enumerate (zip (*get_space_info (las2)))] - spaces = _spin_halfexcitation_products (spaces, spin_halfexcs, nroots_ref=nroots_ref) + spaces = _spin_flip_products (spaces, spin_flips, nroots_ref=nroots_ref) nfrags = spaces[0].nfrag weights = [space.weight for space in spaces] charges = [space.charges for space in spaces] @@ -299,9 +336,68 @@ def charge_excitation_products (las2, las1): # TODO: direct product of single-electron hops raise NotImplementedError (">3-frag LASSIS") +def as_scanner(lsi): + '''Generating a scanner for LASSIS PES. + + The returned solver is a function. This function requires one argument + "mol" as input and returns total LASSIS energy. + + The solver will automatically use the results of last calculation as the + initial guess of the new calculation. All parameters of LASSIS object + are automatically applied in the solver. + + Note scanner has side effects. It may change many underlying objects + (_scf, with_df, with_x2c, ...) during calculation. + ''' + if isinstance(lsi, lib.SinglePointScanner): + return lsi + + logger.info(lsi, 'Create scanner for %s', lsi.__class__) + name = lsi.__class__.__name__ + LASSIS_Scanner.__name_mixin__ + return lib.set_class(LASSIS_Scanner(lsi), (LASSIS_Scanner, lsi.__class__), name) + +class LASSIS_Scanner(lib.SinglePointScanner): + def __init__(self, lsi, state=0): + self.__dict__.update(lsi.__dict__) + self._las = lsi._las.as_scanner() + self._scan_state = state + + def __call__(self, mol_or_geom, **kwargs): + if isinstance(mol_or_geom, gto.MoleBase): + mol = mol_or_geom + else: + mol = self.mol.set_geom_(mol_or_geom, inplace=False) + + self.reset (mol) + for key in ('with_df', 'with_x2c', 'with_solvent', 'with_dftd3'): + sub_mod = getattr(self, key, None) + if sub_mod: + sub_mod.reset(mol) + + las_scanner = self._las + las_scanner(mol) + self.mol = mol + self.mo_coeff = las_scanner.mo_coeff + e_tot = self.kernel()[0][self._scan_state] + if hasattr (e_tot, '__len__'): + e_tot = np.average (e_tot) + return e_tot + class LASSIS (LASSI): def __init__(self, las, ncharge='s', nspin='s', sa_heff=True, deactivate_vrv=False, crash_locmin=False, opt=1, **kwargs): + ''' + Key attributes: + _las : instance of class `LASCINoSymm` + The encapsulated LASSCF wave function. The CI vectors of the reference state are, + i.e., _las.get_single_state_las (state=0).ci. + ci_spin_flips : dict + Keys are (i,s) with integer i and string s = 'u' or 'd'. Values are the spin-up + (s='u') or spin-down (s='d') CI vectors of the ith fragment. + ci_charge_hops: dict + Keys are hash strings for particular charge-hop rootspaces, and values are the CI + vectors of the involved fragments. + ''' self.ncharge = ncharge self.nspin = nspin self.sa_heff = sa_heff @@ -309,9 +405,13 @@ def __init__(self, las, ncharge='s', nspin='s', sa_heff=True, deactivate_vrv=Fal self.crash_locmin = crash_locmin self.e_states_meaningless = True # a tag to silence an invalid warning LASSI.__init__(self, las, opt=opt, **kwargs) + self.max_cycle_macro = 50 + self.conv_tol_self = 1e-6 + self.ci_spin_flips = {} + self.ci_charge_hops = {} if las.nroots>1: - logger.warn (self, ("LASSIS builds the model space for you! I don't know what will " - "happen if you build a model space by hand!")) + logger.warn (self, ("Only the first LASSCF state is used by LASSIS! " + "Other states are discarded!")) def kernel (self, ncharge=None, nspin=None, sa_heff=None, deactivate_vrv=None, crash_locmin=None, **kwargs): @@ -320,9 +420,21 @@ def kernel (self, ncharge=None, nspin=None, sa_heff=None, deactivate_vrv=None, if sa_heff is None: sa_heff = self.sa_heff if deactivate_vrv is None: deactivate_vrv = self.deactivate_vrv if crash_locmin is None: crash_locmin = self.crash_locmin + t0 = (logger.process_clock (), logger.perf_counter ()) + log = logger.new_logger (self, self.verbose) self.converged, las = prepare_states (self, ncharge=ncharge, nspin=nspin, sa_heff=sa_heff, deactivate_vrv=deactivate_vrv, crash_locmin=crash_locmin) - self.__dict__.update(las.__dict__) - return LASSI.kernel (self, **kwargs) + #self.__dict__.update(las.__dict__) # Unsafe + self.fciboxes = las.fciboxes + self.ci = las.ci + self.nroots = las.nroots + self.weights = las.weights + self.e_lexc = las.e_lexc + self.e_states = las.e_states + self.e_roots, self.si = LASSI.kernel (self, **kwargs) + log.timer ("LASSIS", *t0) + return self.e_roots, self.si + + as_scanner = as_scanner diff --git a/my_pyscf/lassi/s2.py b/my_pyscf/lassi/s2.py new file mode 100644 index 00000000..1e905fbf --- /dev/null +++ b/my_pyscf/lassi/s2.py @@ -0,0 +1,57 @@ +import numpy as np + +def gencoup_table (local_smults, global_smult): + ''' Genealogical coupling table for states of 2S+1 = global_smult given len (local_smults) + subunits with local 2S+1 = local_smults ''' + + nfrags = len (local_smults) + local_2s = np.array (local_smults) - 1 + assert (np.sum (local_2s) % 2 != global_smult % 2), 'illegal parity' + assert (np.sum (local_2s)+1 >= global_smult), 'total spin too high' + left_ceiling = np.cumsum (local_2s) + right_ceiling = np.empty_like (left_ceiling) + right_ceiling[:] = global_smult-1 + right_ceiling[1:] += np.cumsum (local_2s[::-1])[:-1] + right_ceiling = right_ceiling[::-1] + ceiling = np.minimum (left_ceiling, right_ceiling) + + def find_lowerable_nodes (path): + path = np.array (path) + idx = np.zeros (nfrags, dtype=np.bool_) + idx[1:-1] = True + idx = idx & (path>1) + left_min = np.array (path) + left_min[1:] = path[:-1] - local_2s[1:] + right_min = np.array (path) + right_min[:-1] = (path - local_2s)[1:] + idx = idx & (path > left_min) + idx = idx & (path > right_min) + return idx + + current_paths = [tuple(ceiling),] + all_paths = set () + for i in range (np.prod (local_smults)): + if not len (current_paths): break + all_paths.update (current_paths) + next_paths = [] + for path in current_paths: + idx_lowerable = find_lowerable_nodes (path) + n_lowerable = np.count_nonzero (idx_lowerable) + if not n_lowerable: continue + idx_lowerable = np.diag (idx_lowerable)[idx_lowerable] + new_paths = np.tile (path, (n_lowerable, 1)) + new_paths[idx_lowerable] -= 2 + next_paths.extend (list (new_paths)) + next_paths = [tuple (next_path) for next_path in next_paths] + next_paths = list (set (next_paths)) + current_paths = next_paths + + all_paths = np.array (list (all_paths)) + for ix in range (2,nfrags): + all_paths = all_paths[all_paths[:,-ix].argsort (kind='mergesort')] + return all_paths + + +if __name__=='__main__': + print (gencoup_table ([5,3,5],10)) + diff --git a/my_pyscf/lassi/sitools.py b/my_pyscf/lassi/sitools.py index 3a180461..d1834a05 100644 --- a/my_pyscf/lassi/sitools.py +++ b/my_pyscf/lassi/sitools.py @@ -2,8 +2,8 @@ from pyscf import lib, symm from scipy import linalg from mrh.my_pyscf.mcscf.lasci import get_space_info -from mrh.my_pyscf.lassi.lassi import ham_2q -from mrh.my_pyscf.lassi.citools import get_lroots +from mrh.my_pyscf.lassi.lassi import ham_2q, root_make_rdm12s, LASSI +from mrh.my_pyscf.lassi.citools import get_lroots, get_rootaddr_fragaddr def decompose_sivec_by_rootspace (las, si, ci=None): '''Decompose a set of LASSI vectors as @@ -189,15 +189,26 @@ def analyze (las, si, ci=None, state=0, print_all_but=1e-8, lbasis='primitive', ) if ci is None: ci = las.ci ci0 = ci + states = np.atleast_1d (state) + nstates = len (states) log = lib.logger.new_logger (las, las.verbose) + log.info ("Natural-orbital analysis for state(s) %s", str (state)) + casdm1s = root_make_rdm12s (las, ci, si, state=state)[0] + if nstates > 1: + casdm1s = casdm1s.sum (0) / nstates + casdm1 = casdm1s.sum (0) + if isinstance (las, LASSI): + no_coeff, no_ene, no_occ = las._las.canonicalize (natorb_casdm1=casdm1)[:3] + else: + no_coeff, no_ene, no_occ = las.canonicalize (natorb_casdm1=casdm1)[:3] + log.info ("Analyzing LASSI vectors for states = %s",str(state)) log.info ("Average quantum numbers:") space_weights, state_coeffs, idx_space = decompose_sivec_by_rootspace ( las, si ) - states = np.atleast_1d (state) nelelas = np.array ([sum (n) for n in las.nelecas_sub]) c, m, smults = get_rootspace_central_moment (las, space_weights[:,states]) neleca = .5*(nelelas[None,:]-c+m) @@ -211,7 +222,6 @@ def analyze (las, si, ci=None, state=0, print_all_but=1e-8, lbasis='primitive', log.info (("Analyzing rootspace fragment density matrices for LASSI " "states %s averaged together"), str (states)) log.info ("Continue until 1-%e of wave function(s) accounted for", print_all_but) - nstates = len (states) avg_weights = space_weights[:,states].sum (1) / nstates lroots = get_lroots (ci0).T running_weight = 1 @@ -445,5 +455,91 @@ def _print_hdiag (log, iroot, state_coeffs, lroots, hdiag0, hdiag_i, log.info ("All %d ENVs in rootspace %d accounted for", nprods, iroot) return +def sivec_fermion_spin_shuffle (si0, nelec_frs, lroots): + '''Permute the signs of rows of the SI vector to account for the difference in convention + between + + ... a2' a1' a0' ... b2' b1' b0' |vac> (spin-major order) + + and + + ... a2' b2' a1' b1' a0' b0' |vac> (fragment-major order) + + Operators and SI vectors in LASSI use spin-major order by default. However, fragment-major + order is sometimes more useful when computing various properties. + + Args: + si0: ndarray of shape (ndim,*) + Contains SI vecotr + nelec_frs: ndarray of shape (nfrags, nroots, 2) + Number of electrons of each spin in each fragment in each rootspace + lroots: ndarray of shape (nfrags, nroots) + Number of states in each fragment in each rootspace + + + Returns: + si1: ndarray of shape (ndim,*) + si0 with permuted row signs + ''' + from mrh.my_pyscf.lassi.op_o1 import fermion_spin_shuffle + nelec_rsf = np.asarray (nelec_frs).transpose (1,2,0) + rootaddr = get_rootaddr_fragaddr (lroots)[0] + si1 = si0.copy () + for iroot, nelec_sf in enumerate (nelec_rsf): + idx = (rootaddr==iroot) + si1[idx,:] *= fermion_spin_shuffle (nelec_sf[0], nelec_sf[1]) + return si1 + +# TODO: resort to spin-major order within the new vacuum! Will probably be +# necessary to make sense of excitons! +def sivec_vacuum_shuffle (si0, nelec_frs, lroots, nelec_vac=None, state=None): + '''Define a particular number of electrons in each fragment as the vacuum, + set the signs of the LAS basis functions accordingly and in fragment-major + order, and return the correspondingly-modified SI vector. + + Args: + si0: ndarray of shape (ndim,*) + Contains SI vecotr + nelec_frs: ndarray of shape (nfrags, nroots, 2) + Number of electrons of each spin in each fragment in each rootspace + lroots: ndarray of shape (nfrags, nroots) + Number of states in each fragment in each rootspace + + Kwargs: + nelec_vac: ndarray of shape (nfrags) + Number of electrons (spinless) in each fragment in the new vacuum. + Defaults to nelec_frs[:,state,:].sum (1) + state: integer + Index of the rootspace identified as the new vacuum. Required if + nelec_vac is unset. + + Returns: + si1: ndarray of shape (ndim,*) + si0 with permuted row signs corresponding to nelec_vac electrons in + each fragment in the vacuum and the fermion creation operators in + fragment-major order + ''' + assert ((nelec_vac is not None) or (state is not None)) + nfrags, nroots = nelec_frs.shape[:2] + si1 = sivec_fermion_spin_shuffle (si0, nelec_frs, lroots) + nelec_rf = nelec_frs.sum (2).T + if nelec_vac is None: nelec_vac = nelec_rf[state] + nelec_vac = np.asarray (nelec_vac) + dparity_rf = np.remainder (nelec_rf - nelec_vac[None,:], 2) + parity_vac = np.remainder (nelec_vac, 2) + rootaddr = get_rootaddr_fragaddr (lroots)[0] + for iroot, dparity_f in enumerate (dparity_rf): + odd_frags = np.where (dparity_f)[0] + nperms = 0 + # Remember: c2' c1' c0' |vac> + for ifrag in odd_frags[::-1]: + if ifrag == nfrags-1: continue + nperms += sum (parity_vac[ifrag+1:]) + idx = (rootaddr==iroot) + si1[idx,:] *= (1,-1)[nperms%2] + return si1 + + + diff --git a/my_pyscf/lassi/states.py b/my_pyscf/lassi/states.py index 799662d7..8862afe0 100644 --- a/my_pyscf/lassi/states.py +++ b/my_pyscf/lassi/states.py @@ -229,24 +229,20 @@ def describe_single_excitation (self, other): e_spin = 'a' if np.any (self.neleca!=other.neleca) else 'b' src_ds = 'u' if self.smults[src_frag]>other.smults[src_frag] else 'd' dest_ds = 'u' if self.smults[dest_frag]>other.smults[dest_frag] else 'd' - return src_frag, dest_frag, e_spin, src_ds, dest_ds + lroots_s = min (other.nelecu[src_frag], other.nholed[dest_frag]) + return src_frag, dest_frag, e_spin, src_ds, dest_ds, lroots_s def single_excitation_description_string (self, other): - src, dest, e_spin, src_ds, dest_ds = self.describe_single_excitation (other) - fmt_str = '{:d}({:s}) --{:s}--> {:d}({:s})' - return fmt_str.format (src, src_ds, e_spin, dest, dest_ds) + src, dest, e_spin, src_ds, dest_ds, lroots_s = self.describe_single_excitation (other) + fmt_str = '{:d}({:s}) --{:s}--> {:d}({:s}) ({:d} lroots)' + return fmt_str.format (src, src_ds, e_spin, dest, dest_ds, lroots_s) def compute_single_excitation_lroots (self, ref): if isinstance (ref, (list, tuple)): lroots = np.array ([self.compute_single_excitation_lroots (r) for r in ref]) return np.amax (lroots) assert (self.is_single_excitation_of (ref)) - src, dest, e_spin = self.describe_single_excitation (ref)[:3] - if e_spin == 'a': - nelec, nhole = ref.neleca, ref.nholea - else: - nelec, nhole = ref.nelecb, ref.nholeb - return min (nelec[src], nhole[dest]) + return self.describe_single_excitation (ref)[5] def is_spin_shuffle_of (self, other): if np.any (self.nelec != other.nelec): return False diff --git a/my_pyscf/mcpdft/__init__.py b/my_pyscf/mcpdft/__init__.py index 05d2a98a..beb29124 100644 --- a/my_pyscf/mcpdft/__init__.py +++ b/my_pyscf/mcpdft/__init__.py @@ -1,5 +1,6 @@ # Lahh dee dah import copy +import numpy as np from pyscf.mcpdft.mcpdft import get_mcpdft_child_class from pyscf import mcscf, gto from pyscf.lib import logger @@ -48,7 +49,6 @@ def CASCIPDFT (mc_or_mf_or_mol, ot, ncas, nelecas, ncore=None, **kwargs): CASSCF=CASSCFPDFT CASCI=CASCIPDFT - # LAS-PDFT def _laspdftEnergy(mc_class, mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, DoLASSI=False, ncore=None, spin_sub=None, frozen=None, **kwargs): @@ -64,16 +64,43 @@ def _laspdftEnergy(mc_class, mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, DoLASSI if frozen is not None: mc1 = mc_class(mf_or_mol, ncas_sub, nelecas_sub, ncore=ncore, spin_sub = spin_sub, frozen=frozen) else: mc1 = mc_class(mf_or_mol, ncas_sub, nelecas_sub, ncore=ncore, spin_sub = spin_sub) - + from mrh.my_pyscf.mcpdft.laspdft import get_mcpdft_child_class mc2 = get_mcpdft_child_class(mc1, ot, DoLASSI=DoLASSI, **kwargs) if mc0 is not None: - mc2.mo_coeff = mc_or_mf_or_mol.mo_coeff.copy () + mc2.mo_coeff = mc_or_mf_or_mol.mo_coeff.copy () mc2.ci = copy.deepcopy (mc_or_mf_or_mol.ci) mc2.converged = mc0.converged return mc2 - + + +def _lassipdftEnergy(mc_class, mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, DoLASSI=False, ncore=None, spin_sub=None, + frozen=None, **kwargs): + + from mrh.my_pyscf.lassi import lassi + + if isinstance(mc_or_mf_or_mol, lassi.LASSI): + mc0 = mc_or_mf_or_mol._las + mf_or_mol = mc_or_mf_or_mol._las._scf + else: + raise "Requires lassi instance" + + mc1 = mc_class(mf_or_mol, ncas_sub, nelecas_sub, ncore=ncore, spin_sub=spin_sub) + + from mrh.my_pyscf.mcpdft.laspdft import get_mcpdft_child_class + mc2 = get_mcpdft_child_class(mc1, ot, DoLASSI=DoLASSI, **kwargs) + + if mc0 is not None: + mc2.mo_coeff = mc_or_mf_or_mol.mo_coeff.copy() + mc2.ci = copy.deepcopy(mc_or_mf_or_mol.ci) + mc2.converged = mc0.converged + _keys = mc_or_mf_or_mol._keys.copy() + mc2.__dict__.update(mc_or_mf_or_mol.__dict__) + mc2._keys = mc2._keys.union(_keys) + mc2.e_mcscf = np.average(mc_or_mf_or_mol.e_roots).copy() + return mc2 + def LASSCFPDFT(mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, ncore=None, spin_sub=None, frozen=None, **kwargs): from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF @@ -83,11 +110,13 @@ def LASSCFPDFT(mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, ncore=None, spin_sub def LASSIPDFT(mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, ncore=None, spin_sub=None, frozen=None, **kwargs): from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF - return _laspdftEnergy(LASSCF, mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, DoLASSI=True, ncore=ncore, + return _lassipdftEnergy(LASSCF, mc_or_mf_or_mol, ot, ncas_sub, nelecas_sub, DoLASSI=True, ncore=ncore, spin_sub=spin_sub, frozen=frozen, **kwargs) + LASSCF = LASSCFPDFT LASSI = LASSIPDFT +LASSIS = LASSIPDFT # Not Sure acc. to issue #34 of mrh def CIMCPDFT (*args, **kwargs): from mrh.my_pyscf.mcpdft.var_mcpdft import CIMCPDFT as fn diff --git a/my_pyscf/mcpdft/laspdft.py b/my_pyscf/mcpdft/laspdft.py index 77a10b98..e9a4e42b 100644 --- a/my_pyscf/mcpdft/laspdft.py +++ b/my_pyscf/mcpdft/laspdft.py @@ -68,13 +68,13 @@ def optimize_mcscf_(self, mo_coeff=None, ci0=None, **kwargs): '''Optimize the MC-SCF wave function underlying an MC-PDFT calculation. Has the same calling signature as the parent kernel method. ''' with _mcscf_env(self): - self.e_mcscf, self.e_cas, self.ci, self.mo_coeff, self.mo_energy = \ - self._mc_class.kernel(self, mo_coeff, ci0=ci0, **kwargs)[:-2] - self.fcisolver.nroots = self.nroots if self.DoLASSI: - self.e_states, self.si = self.lassi() - return self.e_mcscf, self.e_cas, self.ci, self.mo_coeff, self.mo_energy - + self.fcisolver.nroots = len(self.e_states) + self.e_states = self.e_roots + else: + self.e_mcscf, self.e_cas, self.ci, self.mo_coeff, self.mo_energy = \ + self._mc_class.kernel(self, mo_coeff, ci0=ci0, **kwargs)[:-2] + self.fcisolver.nroots = self.nroots pdft = PDFT(mc._scf, mc.ncas_sub, mc.nelecas_sub, my_ot=ot, **kwargs) _keys = pdft._keys.copy() diff --git a/my_pyscf/mcscf/addons.py b/my_pyscf/mcscf/addons.py index 3331c520..18abb18b 100644 --- a/my_pyscf/mcscf/addons.py +++ b/my_pyscf/mcscf/addons.py @@ -60,6 +60,29 @@ def get_h1e_zipped_fcisolver (fcisolver): class FCISolver (fcisolver.__class__, H1EZipFCISolver): + # This is a hack designed to, i.e., set + # self.fcisolvers[i].nroots = 1 for all i + # in the context of StateAverageMixFCISolver._loop_civecs, so that local state-averaging + # in the ith Hilbert space can be totally abstracted away from the HIEZipFCISolver, + # StateAverageNMixFCISolver layer + def _loop_civecs (self, *args, **kwargs): + _solver_args = self._solver_args + _state_args = self._state_args + for i, solver in enumerate (self.fcisolvers): + my_args = [] + for arg in args: + if isinstance (arg, (_state_args, _solver_args)): + my_args.append (arg[i]) + else: + my_args.append (arg) + my_kwargs = {} + for key, item in kwargs.items (): + if isinstance (arg, (_state_args, _solver_args)): + my_kwargs[key] = item[i] + else: + my_kwargs[key] = item + yield solver, my_args, my_kwargs + def kernel(self, h1, h2, norb, nelec, ci0=None, verbose=0, ecore=0, orbsym=None, **kwargs): # Note self.orbsym is initialized lazily in mc1step_symm.kernel function log = logger.new_logger(self, verbose) diff --git a/my_pyscf/mcscf/lasci.py b/my_pyscf/mcscf/lasci.py index 1e3b2806..4aaa897d 100644 --- a/my_pyscf/mcscf/lasci.py +++ b/my_pyscf/mcscf/lasci.py @@ -1,6 +1,7 @@ from pyscf.scf.rohf import get_roothaan_fock from pyscf.fci import cistring from pyscf.mcscf import casci, casci_symm, df +from pyscf.tools import dump_mat from pyscf import symm, gto, scf, ao2mo, lib from pyscf.fci.direct_spin1 import _unpack_nelec from mrh.my_pyscf.mcscf.addons import state_average_n_mix, get_h1e_zipped_fcisolver, las2cas_civec @@ -8,7 +9,7 @@ from mrh.my_pyscf.fci import csf_solver from mrh.my_pyscf.df.sparse_df import sparsedf_array from mrh.my_pyscf.mcscf import chkfile -from mrh.my_pyscf.mcscf.productstate import ImpureProductStateFCISolver +from mrh.my_pyscf.mcscf.productstate import ImpureProductStateFCISolver, state_average_fcisolver from mrh.util.la import matrix_svd_control_options from itertools import combinations from scipy.sparse import linalg as sparse_linalg @@ -180,7 +181,6 @@ def get_grad_ci (las, mo_coeff=None, ci=None, h1eff_sub=None, h2eff_sub=None, ve for isub, (fcibox, h1e, ci0, ncas, nelecas) in enumerate (zip ( las.fciboxes, h1eff_sub, ci, las.ncas_sub, las.nelecas_sub)): eri_cas = las.get_h2eff_slice (h2eff_sub, isub, compact=8) - max_memory = max(400, las.max_memory-lib.current_memory()[0]) linkstrl = fcibox.states_gen_linkstr (ncas, nelecas, True) linkstr = fcibox.states_gen_linkstr (ncas, nelecas, False) h2eff = fcibox.states_absorb_h1e(h1e, eri_cas, ncas, nelecas, .5) @@ -409,7 +409,8 @@ def canonicalize (las, mo_coeff=None, ci=None, casdm1fs=None, natorb_casdm1=None i = sum (ncas_sub[:ix]) j = i + ncas check_diag[i:j,i:j] = 0.0 - if np.amax (np.abs (check_diag)) < 1e-8: + is_block_diag = np.amax (np.abs (check_diag)) < 1e-8 + if is_block_diag: # No off-diagonal RDM elements -> extra effort to prevent diagonalizer from breaking frags for isub, (ncas, nelecas) in enumerate (zip (ncas_sub, nelecas_sub)): i = sum (ncas_sub[:isub]) @@ -460,12 +461,31 @@ def canonicalize (las, mo_coeff=None, ci=None, casdm1fs=None, natorb_casdm1=None h2eff_sub = np.tensordot (ucas, h2eff_sub, axes=((0),(3))).transpose (1,2,3,0) h2eff_sub = h2eff_sub.reshape (nmo*las.ncas, las.ncas, las.ncas) h2eff_sub = lib.numpy_helper.pack_tril (h2eff_sub).reshape (nmo, -1) + + # I/O + log = lib.logger.new_logger (las, las.verbose) + if las.verbose >= lib.logger.INFO: + if is_block_diag: + for isub, nlas in enumerate (ncas_sub): + log.info ("Fragment %d natural orbitals", isub) + i = ncore + sum (ncas_sub[:isub]) + j = i + nlas + log.info ('Natural occ %s', str (mo_occ[i:j])) + log.info ('Natural orbital (expansion on AOs) in CAS space') + label = las.mol.ao_labels() + mo_las = mo_coeff[:,i:j] + dump_mat.dump_rec(log.stdout, mo_las, label, start=1) + else: + log.info ("Delocalized natural orbitals do not reflect LAS fragmentation") + log.info ('Natural occ %s', str (mo_occ[ncore:nocc])) + log.info ('Natural orbital (expansion on AOs) in CAS space') + label = las.mol.ao_labels() + mo_las = mo_coeff[:,ncore:nocc] + dump_mat.dump_rec(log.stdout, mo_las, label, start=1) + return mo_coeff, mo_ene, mo_occ, ci, h2eff_sub def get_init_guess_ci (las, mo_coeff=None, h2eff_sub=None, ci0=None): - # TODO: come up with a better algorithm? This might be working better than what I had before - # but it omits inter-active Coulomb and exchange interactions altogether. Is there a - # non-outer-product algorithm for finding the lowest-energy single product of CSFs? if mo_coeff is None: mo_coeff = las.mo_coeff if ci0 is None: ci0 = [[None for i in range (las.nroots)] for j in range (las.nfrags)] if h2eff_sub is None: h2eff_sub = las.get_h2eff (mo_coeff) @@ -487,14 +507,37 @@ def get_init_guess_ci (las, mo_coeff=None, h2eff_sub=None, ci0=None): eri = eri_cas[i:j,i:j,i:j,i:j] for iy, solver in enumerate (fcibox.fcisolvers): nelec = fcibox._get_nelec (solver, nelecas) - ndet = tuple ([cistring.num_strings (norb, n) for n in nelec]) - if isinstance (ci0[ix][iy], np.ndarray) and ci0[ix][iy].size==ndet[0]*ndet[1]: continue if hasattr (mo_coeff, 'orbsym'): solver.orbsym = mo_coeff.orbsym[ncore+i:ncore+j] hdiag_csf = solver.make_hdiag_csf (h1e, eri, norb, nelec, max_memory=las.max_memory) - ci0[ix][iy] = solver.get_init_guess (norb, nelec, solver.nroots, hdiag_csf)[0] + ci0g = solver.get_init_guess (norb, nelec, solver.nroots, hdiag_csf) + ci0[ix][iy] = las._combine_init_guess_ci (ci0[ix][iy], ci0g, norb, nelec, + solver.nroots) return ci0 +def _combine_init_guess_ci (las, ci0i, ci0g, norb, nelec, nroots): + ''' Function to handle the combination of a generated set of guess CI vectors (ci0g) with + existing CI vectors (ci0i) already stored. ''' + nroots0 = 0 + ndet = tuple ([cistring.num_strings (norb, n) for n in nelec]) + ci0g = np.asarray (ci0g) # TODO: this leads to deprecated behavior in lasscf_rdm + if isinstance (ci0i, np.ndarray) and ci0i.size % ndet[0]*ndet[1] == 0: + ci0i = ci0i.reshape (-1, ndet[0], ndet[1]) + nroots0 = ci0i.shape[0] + if nroots0 == 0: + ci0i = ci0g + elif nroots0 < nroots: + cg = ci0g.reshape (nroots,-1) + cs = ci0i.reshape (nroots0,-1) + ovlp = cg.conj () @ cs.T + cg -= ovlp @ cs + ovlp = cg.conj () @ cg.T + Q, R = linalg.qr (ovlp) + ci0g = Q.T @ ci0g + ci0i = np.append (ci0i, ci0g, axis=0)[:nroots] + if nroots == 1: ci0i = ci0i[0] + return ci0i + def get_space_info (las): ''' Retrieve the quantum numbers defining the states of a LASSCF calculation ''' nfrags, nroots = las.nfrags, las.nroots @@ -505,9 +548,9 @@ def get_space_info (las): nelec = fcibox._get_nelec (solver, las.nelecas_sub[ifrag]) charges[iroot,ifrag] = np.sum (las.nelecas_sub[ifrag]) - np.sum (nelec) spins[iroot,ifrag] = nelec[0]-nelec[1] - smults[iroot,ifrag] = solver.smult + smults[iroot,ifrag] = getattr (solver, 'smult', abs(spins[iroot,ifrag])+1) try: - wfnsyms[iroot,ifrag] = solver.wfnsym or 0 + wfnsyms[iroot,ifrag] = getattr (solver, 'wfnsym', None) or 0 except ValueError as e: wfnsyms[iroot,ifrag] = symm.irrep_name2id (las.mol.groupname, solver.wfnsym) return charges, spins, smults, wfnsyms @@ -540,7 +583,8 @@ def assert_no_duplicates (las, tab=None): raise e from None def state_average_(las, weights=[0.5,0.5], charges=None, spins=None, - smults=None, wfnsyms=None, assert_no_dupes=True): + smults=None, wfnsyms=None, lroots=None, lweights=None, + assert_no_dupes=True): ''' Transform LASCI/LASSCF object into state-average LASCI/LASSCF Args: @@ -569,6 +613,15 @@ def state_average_(las, weights=[0.5,0.5], charges=None, spins=None, wfnsyms[i][j] identifies the point-group irreducible representation Defaults to all zeros (i.e., the totally-symmetric irrep) + lroots: ndarray of shape (nroots,nfrags) + Number of local roots in each fragment for each global state. + The corresponding local weights are set to [1/n,1/n,1/n,...]. + Note that this is transposed compared to the kwarg of run_lasci, + in order to be consistent with the other kwargs of this function! + lweights: list of length nfrags of list of length nroots of sequence + Weights of local roots in each fragment for each global state. + Passing lweights is incompatible with passing lroots. Defaults + to, i.e., np.ones (las.nfrags, las.nroots, 1).tolist () Returns: las: LASCI/LASSCF instance @@ -576,6 +629,7 @@ def state_average_(las, weights=[0.5,0.5], charges=None, spins=None, state-averaged LASCI/LASSCF instance. ''' + # TODO: incorporate lroots counting into get_space_info and simplify this! old_states = np.stack (get_space_info (las), axis=-1) nroots = len (weights) nfrags = las.nfrags @@ -583,6 +637,23 @@ def state_average_(las, weights=[0.5,0.5], charges=None, spins=None, if wfnsyms is None: wfnsyms = np.zeros ((nroots, nfrags), dtype=np.int32) if spins is None: spins = np.asarray ([[n[0]-n[1] for n in las.nelecas_sub] for i in weights]) if smults is None: smults = np.abs (spins)+1 + if lroots is not None and lweights is not None: + raise RuntimeError ("lroots sets lweights: pass either or none but not both") + elif lweights is None: + if lroots is None: + lroots = np.ones ((nfrags, nroots), dtype=int) + else: + lroots = np.asarray (lroots).T + lweights = [] + for i in range (nfrags): + lwi = [] + for j in range (nroots): + lwij = np.ones (lroots[i,j], dtype=float) / lroots[i,j] + lwi.append (lwij) + lweights.append (lwi) + else: + lroots = np.array ([[len (lwij) for lwij in lwi] for lwi in lweights]) + if np.any (lroots>1): raise NotImplementedError ("Local state-averaged LASSCF") charges = np.asarray (charges) wfnsyms = np.asarray (wfnsyms) @@ -611,6 +682,12 @@ def state_average_(las, weights=[0.5,0.5], charges=None, spins=None, las.e_states = np.zeros (nroots) las.nroots = nroots las.weights = weights + for ifrag, fcibox in enumerate (las.fciboxes): + for iroot in range (las.nroots): + if lroots[ifrag][iroot] == 1: continue + fcibox.fcisolvers[iroot] = state_average_fcisolver ( + fcibox.fcisolvers[iroot], weights=lweights[ifrag][iroot] + ) if las.ci is not None: log = lib.logger.new_logger(las, las.verbose) @@ -630,6 +707,7 @@ def state_average_(las, weights=[0.5,0.5], charges=None, spins=None, elif np.count_nonzero (idx) > 1: raise RuntimeError ("Duplicate states specified?\n{}".format (idx)) las.ci = ci0 + las.converged = False return las @@ -638,7 +716,9 @@ def state_average_(las, weights=[0.5,0.5], charges=None, spins=None, See lasci.state_average_ docstring below:\n\n''' + state_average_.__doc__) def state_average (las, weights=[0.5,0.5], charges=None, spins=None, - smults=None, wfnsyms=None, assert_no_dupes=True): + smults=None, wfnsyms=None, lroots=None, lweights=None, assert_no_dupes=True): + is_scanner = isinstance (las, lib.SinglePointScanner) + if is_scanner: las = las.undo_scanner () new_las = las.__class__(las._scf, las.ncas_sub, las.nelecas_sub) new_las.__dict__.update (las.__dict__) new_las.mo_coeff = las.mo_coeff.copy () @@ -649,8 +729,22 @@ def state_average (las, weights=[0.5,0.5], charges=None, spins=None, if las.ci is not None: new_las.ci = [[c2.copy () if isinstance (c2, np.ndarray) else None for c2 in c1] for c1 in las.ci] - return state_average_(new_las, weights=weights, charges=charges, spins=spins, - smults=smults, wfnsyms=wfnsyms, assert_no_dupes=assert_no_dupes) + las = state_average_(new_las, weights=weights, charges=charges, spins=spins, + smults=smults, wfnsyms=wfnsyms, lroots=lroots, lweights=lweights, + assert_no_dupes=assert_no_dupes) + if is_scanner: las = las.as_scanner () + return las + +def get_single_state_las (las, state=0): + ''' Quickly extract a state-specific las calculation from a state-average one ''' + charges, spins, smults, wfnsyms = get_space_info (las) + charges = charges[state:state+1] + spins = spins[state:state+1] + smults = smults[state:state+1] + wfnsyms = wfnsyms[state:state+1] + weights = [1,] + return state_average (las, weights=weights, charges=charges, spins=spins, smults=smults, + wfnsyms=wfnsyms) def run_lasci (las, mo_coeff=None, ci0=None, lroots=None, lweights=None, verbose=0, assert_no_dupes=False, _dry_run=False): @@ -737,8 +831,11 @@ def run_lasci (las, mo_coeff=None, ci0=None, lroots=None, lweights=None, verbose ci0_i = [c[state] for c in ci0] solver = ImpureProductStateFCISolver (fcisolvers, stdout=las.stdout, lweights=[l[state] for l in lweights], verbose=verbose) - # TODO: better handling of CSF symmetry quantum numbers in general for ix, s in enumerate (solver.fcisolvers): + # Set the calling las objects local bottom-layer fcisolvers to the + # locally-state-averaged ones I just made so that I can more easily get + # locally-state-averaged density matrices after this function exits + las.fciboxes[ix].fcisolvers[state] = s i = sum (ncas_sub[:ix]) j = i + ncas_sub[ix] if orbsym is not None: s.orbsym = orbsym[i:j] @@ -954,7 +1051,22 @@ def kernel(self, mo_coeff=None, ci0=None, casdm0_fr=None, conv_tol_grad=None, return self.e_tot, self.e_cas, self.ci, self.mo_coeff, self.mo_energy, h2eff_sub, veff def states_make_casdm1s_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, **kwargs): - ''' Spin-separated 1-RDMs in the MO basis for each subspace in sequence ''' + ''' Get spin-separated, rootspace-separated, locally state-averaged 1-RDMs of the active + orbitals for each fragment in sequence. + + Kwargs: + ci: list of list of ndarrays + Contains CI vectors + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + + Returns: + casdm1frs: list of ndarray of shape (nroots, 2, ncas_sub[i], ncas_sub[i]) + Spin-separated 1-body reduced density matrices for each fragment for each + rootspace, locally state-averaged if applicable. + ''' if ci is None: ci = self.ci if ncas_sub is None: ncas_sub = self.ncas_sub if nelecas_sub is None: nelecas_sub = self.nelecas_sub @@ -971,6 +1083,22 @@ def states_make_casdm1s_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, **k def make_casdm1s_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, casdm1frs=None, w=None, **kwargs): + ''' Get spin-separated, rootspace- and state-averaged 1-RDMs of the active orbitals for + each fragment in sequence. + + Kwargs: + ci: list of list of ndarrays + Contains CI vectors + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + + Returns: + casdm1fs: list of ndarray of shape (2, ncas_sub[i], ncas_sub[i]) + Spin-separated 1-body reduced density matrices for each fragment, + rootspace-averaged and locally state-averaged if applicable. + ''' if casdm1frs is None: casdm1frs = self.states_make_casdm1s_sub (ci=ci, ncas_sub=ncas_sub, nelecas_sub=nelecas_sub, **kwargs) if w is None: w = self.weights @@ -978,6 +1106,24 @@ def make_casdm1s_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, def states_make_casdm1s (self, ci=None, ncas_sub=None, nelecas_sub=None, casdm1frs=None, **kwargs): + ''' Get spin-separated, rootspace-separated, locally state-averaged 1-RDMs of all active + orbitals collectively. + + Kwargs: + ci: list of list of ndarrays + Contains CI vectors + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + casdm1frs: sequence of ndarrays of shape (nroots, 2, ncas_sub[i], ncas_sub[i]) + output of states_make_casdm1s_sub + + Returns: + casdm1rs: ndarray of shape (nroots, 2, ncas, ncas) + Spin-separated 1-body reduced density matrix for the active orbitals for each + rootspace, locally state-averaged if applicable. + ''' if casdm1frs is None: casdm1frs = self.states_make_casdm1s_sub (ci=ci, ncas_sub=ncas_sub, nelecas_sub=nelecas_sub, **kwargs) return np.stack ([np.stack ([linalg.block_diag (*[dm1rs[iroot][ispin] @@ -986,7 +1132,22 @@ def states_make_casdm1s (self, ci=None, ncas_sub=None, nelecas_sub=None, for iroot in range (self.nroots)], axis=0) def states_make_casdm2_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, **kwargs): - ''' Spin-separated 1-RDMs in the MO basis for each subspace in sequence ''' + ''' Get spin-summed, rootspace-separated, locally state-averaged 2-RDMs of the active + orbitals for each fragment in sequence. + + Kwargs: + ci: list of list of ndarrays + Contains CI vectors + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + + Returns: + casdm2fr: list of ndarray of shape [nroots,]+[ncas_sub[i]]*4 + Spin-summed 2-body reduced density matrices for each fragment for each + rootspace, locally state-averaged if applicable. + ''' if ci is None: ci = self.ci if ncas_sub is None: ncas_sub = self.ncas_sub if nelecas_sub is None: nelecas_sub = self.nelecas_sub @@ -996,6 +1157,24 @@ def states_make_casdm2_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, **kw return casdm2 def make_casdm2_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, casdm2fr=None, **kwargs): + ''' Get spin-summed, rootspace- and state-averaged 2-RDMs of the active orbitals for each + fragment in sequence. + + Kwargs: + ci: list of list of ndarrays + Contains CI vectors + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + casdm2fr: sequence of ndarrays of shape [nroots,] + [ncas_sub[i],]*4 + output of states_make_casdm2_sub + + Returns: + casdm2f: list of ndarray of shape [nroots,]+[ncas_sub[i]]*4 + Spin-summed 2-body reduced density matrices for each fragment, rootspace- and + state-averaged if applicable. + ''' if casdm2fr is None: casdm2fr = self.states_make_casdm2_sub (ci=ci, ncas_sub=ncas_sub, nelecas_sub=nelecas_sub, **kwargs) return [np.einsum ('rijkl,r->ijkl', dm2, box.weights) @@ -1003,6 +1182,26 @@ def make_casdm2_sub (self, ci=None, ncas_sub=None, nelecas_sub=None, casdm2fr=No def states_make_rdm1s (self, mo_coeff=None, ci=None, ncas_sub=None, nelecas_sub=None, casdm1rs=None, casdm1frs=None, **kwargs): + ''' Get spin-separated, rootspace-separated, locally state-averaged 1-RDMs in the AO basis. + + Kwargs: + mo_coeff: ndarray of shape (nao,nmo) + Contains MO coefficients + ci: list of list of ndarrays + Contains CI vectors + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + casdm1rs: ndarray of shape (nroots,2,ncas,ncas) + output of states_make_casdm1s + casdm1frs: list of ndarray of shape (nroots,2,ncas_sub[i],ncas_sub[i]) + output of states_make_casdm1s_sub + + Returns: + dm1rs: ndarray of shape (nroots,2,nao,nao) + Rootspace-separated, locally state-averaged, spin-separated 1-RDMs + ''' if mo_coeff is None: mo_coeff = self.mo_coeff if ci is None: ci = self.ci if ncas_sub is None: ncas_sub = self.ncas_sub @@ -1019,13 +1218,35 @@ def states_make_rdm1s (self, mo_coeff=None, ci=None, ncas_sub=None, def make_rdm1s_sub (self, mo_coeff=None, ci=None, ncas_sub=None, nelecas_sub=None, include_core=False, casdm1s_sub=None, **kwargs): + ''' Get spin-separated, rootspace- and state-averaged 1-RDMs of the active orbitals iun the + AO basis for each fragment in sequence. This function is suboptimal but retained for + compatibility with the old DMET-based LASSCF implementation. + + Kwargs: + mo_coeff: ndarray of shape (nao,nmo) + Contains MO coefficients + ci: list of list of ndarrays + Contains CI vectors + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + include_core: logical + If true, density of inactive orbitals is included on return + casdm1s_sub: list of ndarray of shape (2, ncas_sub[i], ncas_sub[i]) + Output of make_casdm1s_sub + + Returns: + dm1fs: list of ndarray of shape (2, nao, nao) + Spin-separated 1-body reduced density matrices for each fragment in the AO basis, + rootspace-averaged and locally state-averaged if applicable. + ''' if mo_coeff is None: mo_coeff = self.mo_coeff if ci is None: ci = self.ci if ncas_sub is None: ncas_sub = self.ncas_sub if nelecas_sub is None: nelecas_sub = self.nelecas_sub if casdm1s_sub is None: casdm1s_sub = self.make_casdm1s_sub (ci=ci, ncas_sub=ncas_sub, nelecas_sub=nelecas_sub, **kwargs) - ''' Same as make_casdm1s_sub, but in the ao basis ''' rdm1s = [] for idx, casdm1s in enumerate (casdm1s_sub): mo = self.get_mo_slice (idx, mo_coeff=mo_coeff) @@ -1040,9 +1261,44 @@ def make_rdm1s_sub (self, mo_coeff=None, ci=None, ncas_sub=None, return rdm1s def make_rdm1_sub (self, **kwargs): + ''' Get spin-summed, rootspace- and state-averaged 1-RDMs of the active orbitals iun the + AO basis for each fragment in sequence. This function is suboptimal but retained for + compatibility with the old DMET-based LASSCF implementation. + + Kwargs: + mo_coeff: ndarray of shape (nao,nmo) + Contains MO coefficients + ci: list of list of ndarrays + Contains CI vectors + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + include_core: logical + If true, density of inactive orbitals is included on return + casdm1fs: list of ndarray of shape (2, ncas_sub[i], ncas_sub[i]) + Output of make_casdm1s_sub + + Returns: + dm1f: list of ndarray of shape (nao, nao) + Spin-summed, 1-body reduced density matrices for each fragment in the AO basis, + rootspace-averaged and locally state-averaged if applicable. + ''' return self.make_rdm1s_sub (**kwargs).sum (1) def make_rdm1s (self, mo_coeff=None, ncore=None, **kwargs): + ''' Get spin-separated, rootspace- and locally state-averaged 1-RDMs in the AO basis. + + Kwargs: + mo_coeff: ndarray of shape (nao,nmo) + Contains MO coefficients + ncore: integer + Number of core orbitals + + Returns: + dm1s: ndarray of shape (2,nao,nao) + Rootspace- and locally state-averaged, spin-separated 1-RDM + ''' if mo_coeff is None: mo_coeff = self.mo_coeff if ncore is None: ncore = self.ncore mo = mo_coeff[:,:ncore] @@ -1052,26 +1308,70 @@ def make_rdm1s (self, mo_coeff=None, ncore=None, **kwargs): return dm_core[None,:,:] + dm_cas def make_rdm1 (self, mo_coeff=None, ci=None, **kwargs): + ''' Get spin-summed, rootspace- and locally state-averaged 1-RDM in the AO basis. + + Kwargs: + mo_coeff: ndarray of shape (nao,nmo) + Contains MO coefficients + ci: list of list of ndarrays + Contains CI vectors + + Returns: + dm1: ndarray of shape (nao,nao) + Rootspace- and locally state-averaged, spin-summed 1-RDM + ''' return self.make_rdm1s (mo_coeff=mo_coeff, ci=ci, **kwargs).sum (0) def make_casdm1s (self, ci=None, **kwargs): - ''' Make the full-dimensional casdm1s spanning the collective active space ''' + ''' Get spin-separated, rootspace- and state-averaged 1-RDMs of all active + orbitals collectively. + + Kwargs: + ci: list of list of ndarrays + Contains CI vectors + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + + Returns: + casdm1s: ndarray of shape (2, ncas, ncas) + Spin-separated, rootspace- and state-averaged 1-body reduced density matrix for the + active orbitals. + ''' casdm1s_sub = self.make_casdm1s_sub (ci=ci, **kwargs) casdm1a = linalg.block_diag (*[dm[0] for dm in casdm1s_sub]) casdm1b = linalg.block_diag (*[dm[1] for dm in casdm1s_sub]) return np.stack ([casdm1a, casdm1b], axis=0) def make_casdm1 (self, ci=None, **kwargs): - ''' Spin-sum make_casdm1s ''' + ''' Get spin-summed, rootspace- and state-averaged 1-RDMs of all active + orbitals collectively. + + Kwargs: + ci: list of list of ndarrays + Contains CI vectors + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + + Returns: + casdm1s: ndarray of shape (ncas, ncas) + Spin-summed, rootspace- and state-averaged 1-body reduced density matrix for the + active orbitals. + ''' return self.make_casdm1s (ci=ci, **kwargs).sum (0) def states_make_casdm2 (self, ci=None, ncas_sub=None, nelecas_sub=None, casdm1frs=None, casdm2fr=None, **kwargs): - ''' Make the full-dimensional casdm2 spanning the collective active space ''' - log = lib.logger.new_logger (self, self.verbose) - log.warn (("You have found yourself in states_make_casdm2, which is " - "a very bad piece of code that Matt should be avoiding. " - "Please yell at him about this at earliest convenience.")) + ''' Make the full-dimensional casdm2 spanning the collective active space. ILLEGAL FUNCTION + DO NOT USE. ''' + raise DeprecationWarning ( + ("states_make_casdm2 is BANNED. There is no reason to EVER make an array this huge.\n" + "Use states_make_casdm*_sub instead, and substitute the factorization into your " + "expressions.") + ) if ci is None: ci = self.ci if ncas_sub is None: ncas_sub = self.ncas_sub if nelecas_sub is None: nelecas_sub = self.nelecas_sub @@ -1106,6 +1406,26 @@ def states_make_casdm2 (self, ci=None, ncas_sub=None, nelecas_sub=None, def state_make_casdm1s(self, ci=None, state=0, ncas_sub=None, nelecas_sub=None, casdm1frs=None, **kwargs): + ''' Get spin-separated, locally state-averaged 1-RDM of all active orbitals collectively + for a single rootspace. + + Kwargs: + ci: list of list of ndarrays + Contains CI vectors + state: integer + Index of rootspace to return + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + casdm1frs: sequence of ndarrays of shape (nroots, 2, ncas_sub[i], ncas_sub[i]) + output of states_make_casdm1s_sub + + Returns: + casdm1s: ndarray of shape (2, ncas, ncas) + Spin-separated 1-body reduced density matrix for the active orbitals for the chosen + rootspace, locally state-averaged if applicable. + ''' if casdm1frs is None: casdm1frs = self.states_make_casdm1s_sub (ci=ci, ncas_sub=ncas_sub, nelecas_sub=nelecas_sub, **kwargs) casdm1s = np.stack([np.stack ([linalg.block_diag (*[dm1rs[iroot][ispin] @@ -1116,16 +1436,97 @@ def state_make_casdm1s(self, ci=None, state=0, ncas_sub=None, nelecas_sub=None, def state_make_casdm2(self, ci=None, state=0, ncas_sub=None, nelecas_sub=None, casdm1frs=None, casdm2fr=None, **kwargs): - ''' State wise casdm2 spanning the collective active space. ''' - # This is producing the casdm2 for all states, but need to generate only for one state - casdm2r = self.states_make_casdm2(ci=ci, ncas_sub=ncas_sub, nelecas_sub=nelecas_sub, - casdm1frs=casdm1frs, casdm2fr=casdm2fr, **kwargs) - return casdm2r[state] + ''' Get spin-summed, locally state-averaged 2-RDM of all active orbitals for a single + rootspace. The formulas for nonzero elements are (state=r): + + casdm2r[p1,p2,p3,p4] = casdm2fr[p][r,p1,p2,p3,p4] + casdm2r[p1,p2,q1,q2] = casdm1fr[p][r,p1,p2] * casdm1fr[q][r,q1,q2] + casdm2r[p1,q2,q1,p2] = -(casdm1frs[p][r,0,p1,p2] * casdm1frs[q][r,0,q1,q2] + + casdm1frs[p][r,1,p1,p2] * casdm1frs[q][r,1,q1,q2]) + + Kwargs: + ci: list of list of ndarrays + Contains CI vectors + state: integer + Index of the chosen rootspace + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + casdm1frs: sequence of ndarrays of shape [nroots,] + [ncas_sub[i],]*2 + output of states_make_casdm1s_sub + casdm2fr: sequence of ndarrays of shape [nroots,] + [ncas_sub[i],]*4 + output of states_make_casdm2_sub + + Returns: + casdm2: ndarray of shape [ncas,]*4 + Spin-summed 2-body reduced density matrices for all active orbitals for the + rootspace identified by "state", locally state-averaged if applicable. + ''' + if ci is None: ci = self.ci + if ncas_sub is None: ncas_sub = self.ncas_sub + if nelecas_sub is None: nelecas_sub = self.nelecas_sub + if casdm1frs is None: casdm1frs = self.states_make_casdm1s_sub (ci=ci) + if casdm2fr is None: casdm2fr = self.states_make_casdm2_sub (ci=ci, + ncas_sub=ncas_sub, nelecas_sub=nelecas_sub, **kwargs) + ncas = sum (ncas_sub) + ncas_cum = np.cumsum ([0] + ncas_sub.tolist ()) + casdm2 = np.zeros ((ncas,ncas,ncas,ncas)) + # Diagonal + for isub, dm2_r in enumerate (casdm2fr): + i = ncas_cum[isub] + j = ncas_cum[isub+1] + casdm2[i:j, i:j, i:j, i:j] = dm2_r[state] + # Off-diagonal + for (isub1, dm1s1_r), (isub2, dm1s2_r) in combinations (enumerate (casdm1frs), 2): + i = ncas_cum[isub1] + j = ncas_cum[isub1+1] + k = ncas_cum[isub2] + l = ncas_cum[isub2+1] + dma1, dmb1 = dm1s1_r[state][0], dm1s1_r[state][1] + dma2, dmb2 = dm1s2_r[state][0], dm1s2_r[state][1] + # Coulomb slice + casdm2[i:j, i:j, k:l, k:l] = np.multiply.outer (dma1+dmb1, dma2+dmb2) + casdm2[k:l, k:l, i:j, i:j] = casdm2[i:j, i:j, k:l, k:l].transpose (2,3,0,1) + # Exchange slice + casdm2[i:j, k:l, k:l, i:j] = -(np.multiply.outer (dma1, dma2) + +np.multiply.outer (dmb1, dmb2)).transpose (0,3,2,1) + casdm2[k:l, i:j, i:j, k:l] = casdm2[i:j, k:l, k:l, i:j].transpose (1,0,3,2) + return casdm2 def make_casdm2 (self, ci=None, ncas_sub=None, nelecas_sub=None, casdm2r=None, casdm2f=None, casdm1frs=None, casdm2fr=None, **kwargs): - ''' Make the full-dimensional casdm2 spanning the collective active space ''' + ''' Get spin-summed, rootspace- and state-averaged 2-RDM of all active orbitals. The + formulas for nonzero elements are + + casdm2r[r,p1,p2,p3,p4] = casdm2fr[p][r,p1,p2,p3,p4] + casdm2r[r,p1,p2,q1,q2] = casdm1fr[p][r,p1,p2] * casdm1fr[q][r,q1,q2] + casdm2r[r,p1,q2,q1,p2] = -(casdm1frs[p][r,0,p1,p2] * casdm1frs[q][r,0,q1,q2] + + casdm1frs[p][r,1,p1,p2] * casdm1frs[q][r,1,q1,q2]) + casdm2 = sum_r casdm2r[r] + + Kwargs: + ci: list of list of ndarrays + Contains CI vectors + ncas_sub: sequence of integers + number of orbitals in each fragment + nelecas_sub: sequence of integers + number of orbitals in each fragment in the reference rootspace + casdm2r: ndarray of shape [nroots,] + [ncas,]*4 + output of states_make_casdm2 + casdm2f: sequence of ndarrays of shape [ncas_sub[i],]*4 + output of make_casdm2_sub + casdm1frs: sequence of ndarrays of shape [nroots,] + [ncas_sub[i],]*2 + output of states_make_casdm1s_sub + casdm2fr: sequence of ndarrays of shape [nroots,] + [ncas_sub[i],]*4 + output of states_make_casdm2_sub + + Returns: + casdm2: ndarray of shape [ncas,]*4 + Spin-summed 2-body reduced density matrices for all active orbitals, rootspace- + and state-averaged if applicable. + ''' if casdm2r is not None: return np.einsum ('rijkl,r->ijkl', casdm2r, self.weights) if ci is None: ci = self.ci @@ -1407,28 +1808,33 @@ def lasci (self, mo_coeff=None, ci0=None, lroots=None, lweights=None, verbose=No return self.converged, self.e_tot, self.e_states, self.e_cas, e_lexc, self.ci @lib.with_doc(run_lasci.__doc__) - def lasci_(self, mo_coeff=None, ci0=None, verbose=None, - assert_no_dupes=False, _dry_run=False): + def lasci_(self, mo_coeff=None, ci0=None, lroots=None, lweights=None, verbose=None, + assert_no_dupes=False, _dry_run=False): if mo_coeff is not None: self.mo_coeff = mo_coeff - return self.lasci (mo_coeff=mo_coeff, ci0=ci0, verbose=verbose, - assert_no_dupes=assert_no_dupes, _dry_run=_dry_run) + return self.lasci (mo_coeff=mo_coeff, ci0=ci0, lroots=lroots, lweights=lweights, + verbose=verbose, assert_no_dupes=assert_no_dupes, _dry_run=_dry_run) state_average = state_average state_average_ = state_average_ + get_single_state_las = get_single_state_las - def lassi(self, **kwargs): + def lassi(self, mo_coeff=None, ci=None, veff_c=None, h2eff_sub=None, orbsym=None, + soc=False, break_symmetry=False, opt=1, **kwargs): #import warnings #lassi_kernel_warn = "Now LASSI have kernel, which takes las instance as input. This [las.lassi()] function " \ # "will be removed soon." #warnings.warn(lassi_kernel_warn, stacklevel=3) from mrh.my_pyscf.lassi import lassi - mylassi = lassi.LASSI(self, **kwargs) - return mylassi.kernel(**kwargs) + mylassi = lassi.LASSI(self, mo_coeff=mo_coeff, ci=ci, soc=soc, opt=opt, + break_symmetry=break_symmetry, **kwargs) + return mylassi.kernel(mo_coeff=mo_coeff, ci=ci, veff_c=veff_c, h2eff_sub=h2eff_sub, + orbsym=orbsym) las2cas_civec = las2cas_civec assert_no_duplicates = assert_no_duplicates get_init_guess_ci = get_init_guess_ci + _combine_init_guess_ci = _combine_init_guess_ci localize_init_guess=lasscf_guess.localize_init_guess def _svd (self, mo_lspace, mo_rspace, s=None, **kwargs): if s is None: s = self._scf.get_ovlp () diff --git a/my_pyscf/mcscf/lasci_sync.py b/my_pyscf/mcscf/lasci_sync.py index e2181768..601c0780 100644 --- a/my_pyscf/mcscf/lasci_sync.py +++ b/my_pyscf/mcscf/lasci_sync.py @@ -263,6 +263,8 @@ def ci_cycle (las, mo, ci0, veff, h2eff_sub, casdm1frs, log): las.nelecas_sub, h1eff_sub, ci0)): eri_cas = las.get_h2eff_slice (h2eff_sub, isub, compact=8) + #max_memory = max(400, las.max_memory-lib.current_memory()[0]) + # Issue #54: compute max_memory here, or in fcisolver? orbsym = getattr (mo, 'orbsym', None) if orbsym is not None: i = ncas_cum[isub] @@ -288,6 +290,7 @@ def ci_cycle (las, mo, ci0, veff, h2eff_sub, casdm1frs, log): e_sub, fcivec = fcibox.kernel(h1e, eri_cas, ncas, nelecas, ci0=fcivec, verbose=log, + #max_memory = max_memory issue #54 ecore=e0, orbsym=orbsym) e_cas.append (e_sub) ci1.append (fcivec) @@ -508,6 +511,7 @@ def _init_df_(h_op): h_op.las.cderi_ao2mo (h_op.mo_coeff, h_op.mo_coeff[:,:h_op.nocc], compact=False)) +# TODO: local state-average generalization class LASCI_HessianOperator (sparse_linalg.LinearOperator): ''' The Hessian-vector product for a `LASCI' energy minimization, implemented as a linear operator from the scipy.sparse.linalg module. `LASCI' here means that the CAS is frozen diff --git a/my_pyscf/mcscf/lasscf_async/crunch.py b/my_pyscf/mcscf/lasscf_async/crunch.py index 52978cbf..651ba3f2 100644 --- a/my_pyscf/mcscf/lasscf_async/crunch.py +++ b/my_pyscf/mcscf/lasscf_async/crunch.py @@ -86,6 +86,13 @@ def skip_value(dic): class ImpuritySCF (scf.hf.SCF): + def _is_mem_enough (self, df_naux=None): + nao = self.mol.nao () + if df_naux is not None: + return 4*df_naux*nao*(nao+1)/1e6+lib.current_memory()[0] < self.max_memory*.95 + else: + return 2*(nao**4)/1e6+lib.current_memory()[0] < self.max_memory*.95 + def _update_space_(self, imporb_coeff, nelec_imp): '''Syntactic sugar for updating the impurity orbital subspace in the encapsulated ImpurityMole object.''' @@ -115,10 +122,15 @@ def _update_impham_1_(self, veff, dm1s, e_tot=None): # Two-electron integrals log = logger.new_logger (self, self.verbose) t0 = (logger.process_clock(), logger.perf_counter()) - if getattr (mf, '_eri', None) is not None: - self._eri = ao2mo.full (mf._eri, imporb_coeff, 4) + conv_eris_mem_error = MemoryError (("Conventional two-electron integrals in asynchronous " + "LASSCF (integral-direct algorithm is not yet " + "supported)")) + df_eris_mem_error = MemoryError (("Density-fitted two-electron integrals in asynchronous " + "LASSCF (outcore algorithm is not yet supported")) if getattr (mf, 'with_df', None) is not None: # TODO: impurity outcore cderi + if not self._is_mem_enough (df_naux = mf.with_df.get_naoaux ()): + raise df_eris_mem_error self.with_df._cderi = np.empty ((mf.with_df.get_naoaux (), nimp*(nimp+1)//2), dtype=imporb_coeff.dtype) ijmosym, mij_pair, moij, ijslice = ao2mo.incore._conc_mos (imporb_coeff, imporb_coeff, @@ -130,6 +142,14 @@ def _update_impham_1_(self, veff, dm1s, e_tot=None): eri2 = ao2mo._ao2mo.nr_e2 (eri1, moij, ijslice, aosym='s2', mosym=ijmosym, out=eri2) b0 = b1 + else: + if getattr (mf, '_eri', None) is None: + if not mf._is_mem_enough (): + raise conv_eris_mem_error + mf._eri = mf.mol.intor('int2e', aosym='s8') + if not self._is_mem_enough (): + raise conv_eris_mem_error + self._eri = ao2mo.full (mf._eri, imporb_coeff, 4) t0 = log.timer ("Two-electron integrals in embedding subspace", *t0) # External mean-field; potentially spin-broken h1s = mf.get_hcore ()[None,:,:] + veff @@ -292,10 +312,11 @@ def casci_kernel(casci, mo_coeff=None, ci0=None, verbose=logger.NOTE, envs=None) (mo_coeff.shape[1], casci.ncore, ncas)) # FCI - max_memory = max(400, casci.max_memory-lib.current_memory()[0]) + #max_memory = max(400, casci.max_memory-lib.current_memory()[0]) + # Issue #54: count memory here, or in FCISolver? e_tot, fcivec = casci.fcisolver.kernel(h1eff, eri_cas, ncas, nelecas, ci0=ci0, verbose=log, - max_memory=max_memory, + #max_memory=max_memory, ecore=energy_core) t1 = log.timer('FCI solver', *t1) diff --git a/my_pyscf/mcscf/lasscf_rdm.py b/my_pyscf/mcscf/lasscf_rdm.py index 7e179416..65fcf736 100644 --- a/my_pyscf/mcscf/lasscf_rdm.py +++ b/my_pyscf/mcscf/lasscf_rdm.py @@ -232,6 +232,7 @@ def kernel (las, mo_coeff=None, casdm1frs=None, casdm2fr=None, conv_tol_grad=1e- if ((norm_gorb < conv_tol_grad) or (norm_gorb < norm_gx/10)) and rdmjk_conv: converged = True break + las.dump_chk (mo_coeff=mo_coeff, ci=[casdm1frs, casdm2fr]) H_op._init_eri_() # Take this part out of the true initialization b/c # if I'm already converged I don't want to waste the cycles t1 = log.timer ('LASSCF Hessian constructor', *t1) @@ -390,7 +391,6 @@ def make_hdiag_csf (self, h1s, h2, norb, nelec, max_memory): return (h1s, h2) def get_init_guess (self, norb, nelec, nroots, ham): - ''' Important: zeroth item selected in get_init_guess_ci ''' h1s, h2 = ham if callable (self._get_init_guess): dm1s, dm2 = self._get_init_guess (norb, nelec, nroots, h1s, h2) @@ -399,7 +399,7 @@ def get_init_guess (self, norb, nelec, nroots, ham): hdiag = fci.make_hdiag_csf (h1s, h2, norb, nelec) ci = fci.get_init_guess (norb, nelec, nroots, hdiag) dm1s, dm2 = self._ci2rdm (fci, ci, norb, nelec) - return [dm1s, dm2], None + return dm1s, dm2 def kernel (self, norb, nelec, h0, h1s, h2): h2 = ao2mo.restore (1, h2, norb) @@ -468,6 +468,15 @@ def make_fcibox (mol, kernel=None, get_init_guess=None, spin=None): s.spin = spin return FCIBox ([s]) +def _combine_init_guess_ci (las, ci0i, ci0g, norb, nelec, nroots): + if nroots>1: raise NotImplementedError ("Multiple local roots for lasscf_rdm") + if getattr (ci0i, '__len__', None) is None: return ci0g + if len (ci0i) != 2: return ci0g + for n, ci0in in enumerate (ci0i): + if not isinstance (ci0in, np.ndarray): return ci0g + if ci0in.size != norb ** (2*n): return ci0g + return ci0i + class LASSCFNoSymm (lasscf_sync_o0.LASSCFNoSymm): def __init__(self, *args, **kwargs): @@ -484,6 +493,7 @@ def __init__(self, *args, **kwargs): _hop = LASSCF_HessianOperator canonicalize = canonicalize rdm_cycle = rdm_cycle + _combine_init_guess_ci = _combine_init_guess_ci def _init_fcibox (self, smult, nel): return make_fcibox (self.mol, spin=nel[0]-nel[1]) diff --git a/my_pyscf/mcscf/mc1step_csf.py b/my_pyscf/mcscf/mc1step_csf.py index c11dc311..5a5b81e8 100644 --- a/my_pyscf/mcscf/mc1step_csf.py +++ b/my_pyscf/mcscf/mc1step_csf.py @@ -58,9 +58,10 @@ def contract_2e(c): if mc.ci_response_space > 7: logger.debug(mc, 'CI step by full response') # full response - max_memory = max(400, mc.max_memory-lib.current_memory()[0]) + #max_memory = max(400, mc.max_memory-lib.current_memory()[0]) + # Issue #54: count memory here, or in fcisolver? e, ci1 = mc.fcisolver.kernel(h1, h2, ncas, nelecas, ecore=ecore, - ci0=ci0, tol=tol, max_memory=max_memory) + ci0=ci0, tol=tol)#, max_memory=max_memory) else: # MRH 03/24/2019: this is where I intervene to enforce CSFs fci = mc.fcisolver diff --git a/my_pyscf/mcscf/productstate.py b/my_pyscf/mcscf/productstate.py index 158f55f6..fd437109 100644 --- a/my_pyscf/mcscf/productstate.py +++ b/my_pyscf/mcscf/productstate.py @@ -34,14 +34,17 @@ def kernel (self, h1, h2, norb_f, nelec_f, ecore=0, ci0=None, orbsym=None, log.info ('Entering product-state fixed-point CI iteration') for it in range (max_cycle_macro): ci0 = ci1 - h1eff, h0eff = self.project_hfrag (h1, h2, ci0, norb_f, nelec_f, + h1eff, h0eff, ci0 = self.project_hfrag (h1, h2, ci0, norb_f, nelec_f, ecore=ecore, **kwargs) grad = self._get_grad (h1eff, h2, ci0, norb_f, nelec_f, **kwargs) grad_max = np.amax (np.abs (grad)) - log.info ('Cycle %d: max grad = %e ; sigma = %e', it, grad_max, - e_sigma) - log.debug ('e vector = {}'.format (e)) solvers_converged = [np.all (np.asarray (s.converged)) for s in self.fcisolvers] + nconv = sum ([int (c) for c in solvers_converged]) + log.info ('Cycle %d: max grad = %e ; sigma = %e ; %d/%d fragment CI solvers converged', + it, grad_max, e_sigma, nconv, len (self.fcisolvers)) + log.debug ('e vector = {}'.format (e)) + if nconv0): converged = True @@ -63,12 +66,25 @@ def get_init_guess (self, ci0, norb_f, nelec_f, h1, h2): if h1.ndim < 3: h1 = np.stack ([h1, h1], axis=0) for ix, (no, ne, solver) in enumerate (zip (norb_f, nelec_f, self.fcisolvers)): nelec = self._get_nelec (solver, ne) + i = sum (norb_f[:ix]) + j = i + norb_f[ix] + hdiag_csf = solver.make_hdiag_csf (h1[:,i:j,i:j], h2[i:j,i:j,i:j,i:j], + no, nelec) + ci1_guess = solver.get_init_guess (no, nelec, solver.nroots, hdiag_csf) + na = cistring.num_strings (no, nelec[0]) + nb = cistring.num_strings (no, nelec[1]) if ci1[ix] is None: - i = sum (norb_f[:ix]) - j = i + norb_f[ix] - hdiag_csf = solver.make_hdiag_csf (h1[:,i:j,i:j], h2[i:j,i:j,i:j,i:j], - no, nelec) - ci1[ix] = solver.get_init_guess (no, nelec, solver.nroots, hdiag_csf) + ci1[ix] = ci1_guess + elif np.asarray (ci1[ix]).reshape (-1,na*nb).shape[0] < solver.nroots: + ci1_inp = np.asarray (ci1[ix]).reshape (-1,na*nb) + ci1_guess = np.asarray (ci1_guess).reshape (-1,na*nb) + ovlp = ci1_inp.conj () @ ci1_guess.T + ci1_guess -= ovlp.T @ ci1_inp + ovlp = ci1_guess.conj () @ ci1_guess.T + Q, R = linalg.qr (ovlp) + ci1_guess = Q.T @ ci1_guess + ci1[ix] = np.append (ci1_inp, ci1_guess, axis=0)[:solver.nroots].reshape ( + solver.nroots, na, nb) return self._check_init_guess (ci1, norb_f, nelec_f) def _check_init_guess (self, ci0, norb_f, nelec_f): @@ -223,7 +239,8 @@ def project_hfrag (self, h1, h2, ci, norb_f, nelec_f, ecore=0, dm1s=None, dm2=No e_i -= (np.tensordot (h1_i, dm1s_i, axes=3) + 0.5*np.tensordot (h2_i, dm2_i, axes=4)) h0eff.append (e_i) - return h1eff, h0eff + return h1eff, h0eff, ci + # A child class will modify ci in this function def make_rdm1s (self, ci, norb_f, nelec_f, **kwargs): norb = sum (norb_f) diff --git a/my_pyscf/tools/molden.py b/my_pyscf/tools/molden.py index 0a8d9a92..30414e6a 100644 --- a/my_pyscf/tools/molden.py +++ b/my_pyscf/tools/molden.py @@ -46,12 +46,16 @@ def from_lasscf (las, fname, state=None, natorb_casdm1=None, **kwargs): mo_coeff, mo_ene, mo_occ = las.canonicalize (natorb_casdm1=natorb_casdm1)[:3] return from_mo (las.mol, fname, mo_coeff, occ=mo_occ, ene=mo_ene, **kwargs) -def from_lassi (las, fname, state=0, si=None, opt=1, **kwargs): - if si is None: si = getattr (las, 'si', None) +def from_lassi (lsi, fname, state=0, si=None, opt=1, **kwargs): + if si is None: si = getattr (lsi, 'si', None) from mrh.my_pyscf.lassi.lassi import root_make_rdm12s - natorb_casdm1 = root_make_rdm12s (las, las.ci, si, state=state, opt=opt)[0].sum (0) - mo_coeff, mo_ene, mo_occ = las.canonicalize (natorb_casdm1=natorb_casdm1)[:3] - return from_mo (las.mol, fname, mo_coeff, occ=mo_occ, ene=mo_ene, **kwargs) + from mrh.my_pyscf.lassi import LASSI + natorb_casdm1 = root_make_rdm12s (lsi, lsi.ci, si, state=state, opt=opt)[0].sum (0) + if isinstance (lsi, LASSI): + mo_coeff, mo_ene, mo_occ = lsi._las.canonicalize (natorb_casdm1=natorb_casdm1)[:3] + else: + mo_coeff, mo_ene, mo_occ = lsi.canonicalize (natorb_casdm1=natorb_casdm1)[:3] + return from_mo (lsi.mol, fname, mo_coeff, occ=mo_occ, ene=mo_ene, **kwargs) def from_mcscf (mc, filename, ignore_h=IGNORE_H, cas_natorb=False): ncore, ncas = mc.ncore, mc.ncas diff --git a/pyscf-forge_version.txt b/pyscf-forge_version.txt index 1e589b06..4790f30f 100644 --- a/pyscf-forge_version.txt +++ b/pyscf-forge_version.txt @@ -1 +1 @@ -git+https://github.com/pyscf/pyscf-forge.git@1077e48d5efab407eb6116701ba699ccaa107192 +git+https://github.com/pyscf/pyscf-forge.git@f817911cb0c984207d9309a1d7347f6dbb46a9c6 diff --git a/pyscf_version.txt b/pyscf_version.txt index e2efaa5d..22f2bd9e 100644 --- a/pyscf_version.txt +++ b/pyscf_version.txt @@ -1 +1 @@ -git+https://github.com/pyscf/pyscf.git@0bc590ffddec98eec5d2513a3d85274b7a791bc6 +git+https://github.com/pyscf/pyscf.git@10f89c376371ed55b99075cd77ba209181629ca2 diff --git a/tests/lasscf/test_lasci_mod.py b/tests/lasscf/test_lasci_mod.py index 9bf1f118..a02c81a9 100644 --- a/tests/lasscf/test_lasci_mod.py +++ b/tests/lasscf/test_lasci_mod.py @@ -70,30 +70,32 @@ def build (mf, m1=0, m2=0, ir1=0, ir2=0, CASlist=None, active_first=False, calcn # -------------------------------------------------------------------------------------------------------------------- return c2h4n4_dmet -dr_nn = 2.0 -mol = struct (dr_nn, dr_nn, '6-31g', symmetry=False) -mol.verbose = lib.logger.DEBUG -mol.output = '/dev/null' -mol.spin = 0 -mol.build () -mf = scf.RHF (mol).run () -dmet = build (mf, 1, -1, active_first=False) -dmet.conv_tol_grad = 1e-5 -try: - e_tot = dmet.doselfconsistent () -except RuntimeError as e: - e_tot = 0.0 - -#np.savetxt ('test_lasci_mo.dat', dmet.las.mo_coeff) -#np.savetxt ('test_lasci_ci0.dat', dmet.las.ci0) -#np.savetxt ('test_lasci_ci1.dat', dmet.las.ci[1]) -dmet.las.mo_coeff = np.loadtxt (os.path.join (topdir, 'test_lasci_mo.dat')) -dmet.las.ci[0] = [np.loadtxt (os.path.join (topdir, 'test_lasci_ci0.dat'))] -dmet.las.ci[1] = [-np.loadtxt (os.path.join (topdir, 'test_lasci_ci1.dat')).T] -ugg = LASCI_UnitaryGroupGenerators (dmet.las, dmet.las.mo_coeff, dmet.las.ci) -h_op = LASCI_HessianOperator (dmet.las, ugg) -np.random.seed (0) -x = np.random.rand (ugg.nvar_tot) +def setUpModule(): + global mol, mf, dmet, ugg, h_op, x + dr_nn = 2.0 + mol = struct (dr_nn, dr_nn, '6-31g', symmetry=False) + mol.verbose = lib.logger.DEBUG + mol.output = '/dev/null' + mol.spin = 0 + mol.build () + mf = scf.RHF (mol).run () + dmet = build (mf, 1, -1, active_first=False) + dmet.conv_tol_grad = 1e-5 + try: + e_tot = dmet.doselfconsistent () + except RuntimeError as e: + e_tot = 0.0 + + #np.savetxt ('test_lasci_mo.dat', dmet.las.mo_coeff) + #np.savetxt ('test_lasci_ci0.dat', dmet.las.ci0) + #np.savetxt ('test_lasci_ci1.dat', dmet.las.ci[1]) + dmet.las.mo_coeff = np.loadtxt (os.path.join (topdir, 'test_lasci_mo.dat')) + dmet.las.ci[0] = [np.loadtxt (os.path.join (topdir, 'test_lasci_ci0.dat'))] + dmet.las.ci[1] = [-np.loadtxt (os.path.join (topdir, 'test_lasci_ci1.dat')).T] + ugg = LASCI_UnitaryGroupGenerators (dmet.las, dmet.las.mo_coeff, dmet.las.ci) + h_op = LASCI_HessianOperator (dmet.las, ugg) + np.random.seed (0) + x = np.random.rand (ugg.nvar_tot) def tearDownModule(): global mol, mf, dmet, ugg, h_op, x diff --git a/tests/lasscf/test_lasscf_rdm.py b/tests/lasscf/test_lasscf_rdm.py index 9f361757..a1eb217e 100644 --- a/tests/lasscf/test_lasscf_rdm.py +++ b/tests/lasscf/test_lasscf_rdm.py @@ -20,16 +20,20 @@ from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF as LASSCFRef from mrh.my_pyscf.mcscf.lasscf_rdm import LASSCF as LASSCFTest -xyz = '''H 0.0 0.0 0.0 - H 1.0 0.0 0.0 - H 0.2 3.9 0.1 - H 1.159166 4.1 -0.1''' -mol = gto.M (atom = xyz, basis = '6-31g', output='lasscf_rdm.log', - verbose=lib.logger.INFO) -mf = scf.RHF (mol).run () -las_ref = LASSCFRef (mf, (2,2), (2,2), spin_sub=(1,1)) -las_test = LASSCFTest (mf, (2,2), (2,2), spin_sub=(1,1)) -mo_loc = las_ref.localize_init_guess (((0,1),(2,3)), mf.mo_coeff) +def setUpModule(): + global mol, mf, las_test, las_ref, mo_loc + xyz = '''H 0.0 0.0 0.0 + H 1.0 0.0 0.0 + H 0.2 3.9 0.1 + H 1.159166 4.1 -0.1''' + mol = gto.M (atom = xyz, basis = '6-31g', output='lasscf_rdm.log', + verbose=lib.logger.INFO) + mf = scf.RHF (mol).run () + las_ref = LASSCFRef (mf, (2,2), (2,2), spin_sub=(1,1)) + las_ref.chkfile = 'lasscf_rdm_ref.chk' + las_test = LASSCFTest (mf, (2,2), (2,2), spin_sub=(1,1)) + las_test.chkfile = 'lasscf_rdm_test.chk' + mo_loc = las_ref.localize_init_guess (((0,1),(2,3)), mf.mo_coeff) def tearDownModule(): global mol, mf, las_test, las_ref, mo_loc diff --git a/tests/lasscf/test_old_lasscf.py b/tests/lasscf/test_old_lasscf.py index eeeb990d..9fa0c730 100644 --- a/tests/lasscf/test_old_lasscf.py +++ b/tests/lasscf/test_old_lasscf.py @@ -66,22 +66,24 @@ def build (mf, m1=0, m2=0, ir1=0, ir2=0, CASlist=None, active_first=False, calcn # -------------------------------------------------------------------------------------------------------------------- return c2h4n4_dmet -dr_nn = 3.0 -mol = struct (dr_nn, dr_nn, '6-31g', symmetry=False) -mol.verbose = lib.logger.DEBUG -mol.output = '/dev/null' -mol.spin = 8 -mol.build () -mf = scf.RHF (mol).run () -dmet = build (mf, 1, -1, active_first=True) -dmet.conv_tol_grad = 1e-6 -e_tot = dmet.doselfconsistent () +def setUpModule(): + global mol, mf, dmet, e_tot + dr_nn = 3.0 + mol = struct (dr_nn, dr_nn, '6-31g', symmetry=False) + mol.verbose = lib.logger.DEBUG + mol.output = '/dev/null' + mol.spin = 8 + mol.build () + mf = scf.RHF (mol).run () + dmet = build (mf, 1, -1, active_first=True) + dmet.conv_tol_grad = 1e-6 + e_tot = dmet.doselfconsistent () def tearDownModule(): - global mol, mf, dmet + global mol, mf, dmet, e_tot mol.stdout.close () dmet.lasci_log.close () - del mol, mf, dmet + del mol, mf, dmet, e_tot class KnownValues(unittest.TestCase): diff --git a/tests/lasscf/test_utilities.py b/tests/lasscf/test_utilities.py new file mode 100644 index 00000000..5837ef40 --- /dev/null +++ b/tests/lasscf/test_utilities.py @@ -0,0 +1,35 @@ +import unittest +from pyscf import gto, scf + +def setUpModule(): + pass + +def tearDownModule(): + pass + +class KnownValues (unittest.TestCase): + + def test_get_single_state_las (self): + from mrh.my_pyscf.mcscf.lasscf_async import LASSCF + from mrh.my_pyscf.lassi import states as lassi_states + xyz='''N 0 0 0, + N 3 0 0''' + mol = gto.M (atom=xyz, basis='cc-pvdz', symmetry=False, verbose=0, output='/dev/null') + mf = scf.RHF (mol).run () + las = LASSCF (mf, (6,6), ((3,0),(0,3))) + + mo = las.set_fragments_(([0],[1])) + las = lassi_states.spin_shuffle (las) + las.weights = [1.0/las.nroots,]*las.nroots + las.kernel (mo) + for i in range (len (las.e_states)): + las1 = las.get_single_state_las (state=i) + las1.lasci () + self.assertAlmostEqual (las1.e_tot, las.e_states[i], 6) + e_lasci = las1.e_tot + las1.kernel () + self.assertLessEqual (las1.e_tot, las.e_states[i]) + +if __name__ == "__main__": + print("Full Tests for LASSCF/LASCI miscellaneous") + unittest.main() diff --git a/tests/lassi/test_22.py b/tests/lassi/test_22.py index 408f4222..64481215 100644 --- a/tests/lassi/test_22.py +++ b/tests/lassi/test_22.py @@ -18,14 +18,14 @@ from scipy import linalg from pyscf import lib, gto, scf, mcscf, ao2mo from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF -from mrh.my_pyscf.lassi import LASSI +from mrh.my_pyscf.lassi import LASSI, LASSIrq, LASSIrqCT from mrh.my_pyscf.lassi.lassi import root_make_rdm12s, make_stdm12s from mrh.my_pyscf.lassi.states import all_single_excitations, SingleLASRootspace from mrh.my_pyscf.mcscf.lasci import get_space_info from mrh.my_pyscf.lassi import op_o0, op_o1, lassis def setUpModule (): - global mol, mf, lsi, las, op + global mol, mf, lsi, las, mc, op xyz='''H 0 0 0 H 1 0 0 H 3 0 0 @@ -55,12 +55,15 @@ def setUpModule (): lsi = LASSI (las1) lsi.kernel (opt=0) + # CASCI limit + mc = mcscf.CASCI (mf, 4, 4).run () + op = (op_o0, op_o1) def tearDownModule(): - global mol, mf, lsi, las, op + global mol, mf, lsi, las, mc, op mol.stdout.close () - del mol, mf, lsi, las, op + del mol, mf, lsi, las, mc, op class KnownValues(unittest.TestCase): @@ -73,7 +76,6 @@ def test_op_o1 (self): def test_casci_limit (self): # CASCI limit - mc = mcscf.CASCI (mf, 4, 4).run () casdm1, casdm2 = mc.fcisolver.make_rdm12 (mc.ci, mc.ncas, mc.nelecas) # LASSI in the CASCI limit @@ -98,6 +100,14 @@ def test_casci_limit (self): self.assertAlmostEqual (lib.fp (stdm1s[0,:,2:,2:,0]), lib.fp (stdm1s[1,:,2:,2:,1])) + def test_lassirq (self): + lsi1 = LASSIrq (las, 2, 3).run () + self.assertAlmostEqual (lsi1.e_roots[0], mc.e_tot, 8) + + def test_lassirqct (self): + lsi1 = LASSIrqCT (las, 2, 3).run () + self.assertAlmostEqual (lsi1.e_roots[0], -4.2879945248402445, 8) + def test_contract_hlas_ci (self): e_roots, si, las = lsi.e_roots, lsi.si, lsi._las h0, h1, h2 = lsi.ham_2q () diff --git a/tests/lassi/test_c2h4n4.py b/tests/lassi/test_c2h4n4.py index cb539487..3ce330de 100644 --- a/tests/lassi/test_c2h4n4.py +++ b/tests/lassi/test_c2h4n4.py @@ -181,7 +181,7 @@ def test_lassis_slow (self): for opt in (0,1): with self.subTest (opt=opt): lsis = LASSIS (las1).run (opt=opt) - self.assertAlmostEqual (lsis.e_roots[0], -295.52103109, 7) + self.assertAlmostEqual (lsis.e_roots[0], -295.5210783894406, 7) self.assertTrue (lsis.converged) if __name__ == "__main__": diff --git a/tests/lassi/test_excitations.py b/tests/lassi/test_excitations.py index f5d61f90..6730e475 100644 --- a/tests/lassi/test_excitations.py +++ b/tests/lassi/test_excitations.py @@ -12,6 +12,7 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. +import sys import copy import unittest import numpy as np @@ -23,10 +24,21 @@ from mrh.my_pyscf.lassi import LASSI, op_o0, op_o1 from mrh.my_pyscf.lassi.lassi import root_make_rdm12s, make_stdm12s from mrh.my_pyscf.lassi.states import all_single_excitations -from mrh.my_pyscf.lassi.excitations import ExcitationPSFCISolver, only_ground_states +from mrh.my_pyscf.lassi.excitations import ExcitationPSFCISolver from mrh.my_pyscf.mcscf.lasci import get_space_info from mrh.my_pyscf.mcscf.productstate import ImpureProductStateFCISolver +def only_ground_states (ci0): + '''For a list of sequences of CI vectors in the same Hilbert space, + generate a list in which all but the first element of each sequence + is discarded.''' + ci1 = [] + for c in ci0: + c = np.asarray (c) + if c.ndim==3: c = c[0] + ci1.append (c) + return ci1 + def setUpModule (): global mol, mf, lsi, op xyz='''He 0 0 0 @@ -34,7 +46,7 @@ def setUpModule (): H 11 0 0 H 10 2 0 H 11 2 0''' - mol = gto.M (atom=xyz, basis='sto3g', symmetry=False, verbose=0, output='/dev/null') + mol = gto.M (atom=xyz, basis='sto3g', symmetry=False, verbose=5, output='test_excitations.log')#0, output='/dev/null') mf = scf.RHF (mol).run () # LASSCF with CASCI-limit model space @@ -120,7 +132,7 @@ def lassi_ref (ci1, iroot): weights[0] = 1 psexc.set_excited_fragment_(1+i, (neleca[iroot,i], nelecb[iroot,i]), smults[iroot,i], weights=weights) - conv, energy_tot, ci1 = psexc.kernel (h1, h2, ecore=h0) + conv, energy_tot, ci1 = psexc.kernel (h1, h2, ecore=h0, _add_vrv_energy=True) self.assertTrue (conv) e_roots1, si1 = lassi_ref (ci1, iroot) idx_match = np.argmin (np.abs (e_roots1-energy_tot)) @@ -136,7 +148,7 @@ def lassi_ref (ci1, iroot): weights[0] = 1 psexc.set_excited_fragment_(1+i, (neleca[iroot,i], nelecb[iroot,i]), smults[iroot,i], weights=weights) - conv, energy_tot, ci1 = psexc.kernel (h1, h2, ecore=h0) + conv, energy_tot, ci1 = psexc.kernel (h1, h2, ecore=h0, _add_vrv_energy=True) self.assertTrue (conv) self.assertAlmostEqual (energy_tot, lsi._las.e_states[iroot], 8) @@ -167,7 +179,8 @@ def test_multiref (self): for i in range (3)] ci_ref = las.ci nelec_ref = [[1,1] for i in range (3)] - psexc = ExcitationPSFCISolver (psref, ci_ref, las.ncas_sub, nelec_ref) + psexc = ExcitationPSFCISolver (psref, ci_ref, las.ncas_sub, nelec_ref, + stdout=mol.stdout, verbose=mol.verbose) charges, spins, smults, wfnsyms = get_space_info (lsi._las) dneleca = (spins - charges) // 2 dnelecb = -(charges + spins) // 2 @@ -178,6 +191,10 @@ def test_multiref (self): neleca = dneleca + neleca[None,1:] nelecb = dnelecb + nelecb[None,1:] lroots = lsi.get_lroots () + lroots[:,:] = 1 + # Convergence failures if lroots>1 because I'm not state-averaging, + # but if I state-average then the equality between energy_tot and the + # ref is broken, so I can only have 1 lroot here smults_rf = dsmults + 1 ci0_ref = [only_ground_states (c) for c in ci_ref] @@ -207,11 +224,11 @@ def lassi_ref (ci1, iroot): for iroot in range (1, 5): with self.subTest (rootspace=iroot): for i in range (2): - weights = np.zeros (lroots[i,iroot]) - weights[0] = 1 + weights = np.ones (lroots[i,iroot]) / lroots[i,iroot] psexc.set_excited_fragment_(1+i, (neleca[iroot,i], nelecb[iroot,i]), smults[iroot,i], weights=weights) - conv, energy_tot, ci1 = psexc.kernel (h1, h2, ecore=h0) + conv, energy_tot, ci1 = psexc.kernel (h1, h2, ecore=h0, + _add_vrv_energy=True) self.assertTrue (conv) e_roots1, si1 = lassi_ref (ci1, iroot) idx_match = np.argmin (np.abs (e_roots1-energy_tot)) diff --git a/tests/lassi/test_lassis_targets_slow.py b/tests/lassi/test_lassis_targets_slow.py new file mode 100644 index 00000000..7e68c511 --- /dev/null +++ b/tests/lassi/test_lassis_targets_slow.py @@ -0,0 +1,262 @@ +#!/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. +import os +import copy +import unittest +import numpy as np +from pyscf import gto, scf, mcscf +from pyscf.lib import chkfile +from pyscf.mcscf import avas +from pyscf.data import nist +from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF +from mrh.my_pyscf.mcscf import lasscf_async +from mrh.tests.lasscf.c2h4n4_struct import structure as struct +from mrh.my_pyscf import lassi +topdir = os.path.abspath (os.path.join (__file__, '..')) + +au2cm = nist.HARTREE2J / nist.PLANCK / nist.LIGHT_SPEED_SI * 1e-2 +def yamaguchi (e_roots, s2, nelec=6): + '''The Yamaguchi formula for nelec unpaired electrons''' + hs = (nelec/2.0) * ((nelec/2.0)+1.0) + ls = hs - np.floor (hs) + idx = np.argsort (e_roots) + e_roots = e_roots[idx] + s2 = s2[idx] + idx_hs = (np.around (s2) == np.around (hs)) + assert (np.count_nonzero (idx_hs)), 'high-spin ground state not found' + idx_hs = np.where (idx_hs)[0][0] + e_hs = e_roots[idx_hs] + idx_ls = (np.around (s2) == np.around (ls)) + assert (np.count_nonzero (idx_ls)), 'low-spin ground state not found' + idx_ls = np.where (idx_ls)[0][0] + e_ls = e_roots[idx_ls] + j = (e_ls - e_hs) / (hs - ls) + return j*au2cm + +def setUpModule (): + pass + +def tearDownModule(): + pass + +class KnownValues(unittest.TestCase): + + def test_c2h4n4_3frag (self): + # str + mol = struct (2.0, 2.0, '6-31g', symmetry=False) + mol.output = '/dev/null' + mol.verbose = 0 + mol.spin = 8 + mol.build () + mf = scf.RHF (mol).run () + las = LASSCF (mf, (4,2,4), ((2,2),(1,1),(2,2)), spin_sub=(1,1,1)) + mo_coeff = las.localize_init_guess ([[0,1,2],[3,4,5,6],[7,8,9]]) + las.kernel (mo_coeff) + lsi = lassi.LASSIS (las).run () + e_str = lsi.e_roots[0] + with self.subTest ('str converged'): + self.assertTrue (lsi.converged) + # equil + mol = struct (0.0, 0.0, '6-31g', symmetry=False) + mol.spin = 0 + mol.verbose = 0 + mol.output = '/dev/null' + mol.build () + mf = scf.RHF (mol).run () + las = LASSCF (mf, (4,2,4), ((2,2),(1,1),(2,2)), spin_sub=(1,1,1)) + mo_coeff = las.sort_mo ([7,8,16,18,22,23,24,26,33,34]) + mo_coeff = las.localize_init_guess ([[0,1,2],[3,4,5,6],[7,8,9]], mo_coeff=mo_coeff) + las.kernel (mo_coeff) + lsi = lassi.LASSIS (las).run () + e_equil = lsi.e_roots[0] + with self.subTest ('equil converged'): + self.assertTrue (lsi.converged) + # test + de = 1000 * (e_str - e_equil) + self.assertAlmostEqual (de, 208.21775437138967, 1) + + def test_c2h4n4_2frag (self): + # str + mol = struct (2.0, 2.0, '6-31g', symmetry=False) + mol.output = '/dev/null' + mol.verbose = 0 + mol.spin = 8 + mol.build () + mf = scf.RHF (mol).run () + las = lasscf_async.LASSCF (mf, (5,5), ((3,2),(2,3))) + las = las.state_average ([0.5,0.5], + spins=[[1,-1],[-1,1]], + smults=[[2,2],[2,2]], + charges=[[0,0],[0,0]], + wfnsyms=[[0,0],[0,0]]) + mo = las.set_fragments_((list (range (5)), list (range (5,10)))) + las.kernel (mo) + mo_coeff = las.mo_coeff + las = lasscf_async.LASSCF (mf, (5,5), ((3,2),(2,3))) + las.lasci_(mo_coeff) + lsi = lassi.LASSIS (las).run () + e_str = lsi.e_roots[0] + with self.subTest ('str converged'): + self.assertTrue (lsi.converged) + # equil + mol = struct (0.0, 0.0, '6-31g', symmetry=False) + mol.spin = 0 + mol.verbose = 0 + mol.output = '/dev/null' + mol.build () + mf = scf.RHF (mol).run () + las = lasscf_async.LASSCF (mf, (5,5), ((3,2),(2,3))) + las = las.state_average ([0.5,0.5], + spins=[[1,-1],[-1,1]], + smults=[[2,2],[2,2]], + charges=[[0,0],[0,0]], + wfnsyms=[[0,0],[0,0]]) + mo = las.sort_mo ([7,8,16,18,22,23,24,26,33,34]) + mo = las.set_fragments_((list (range (5)), list (range (5,10))), mo) + las.kernel (mo) + mo_coeff = las.mo_coeff + las = lasscf_async.LASSCF (mf, (5,5), ((3,2),(2,3))) + las.lasci_(mo_coeff) + lsi = lassi.LASSIS (las).run () + e_equil = lsi.e_roots[0] + with self.subTest ('equil converged'): + self.assertTrue (lsi.converged) + # test + de = 1000 * (e_str - e_equil) + self.assertAlmostEqual (de, 191.185141573740, 1) + + def test_kremer_cr2_model (self): + xyz='''Cr -1.320780000000 0.000050000000 -0.000070000000 + Cr 1.320770000000 0.000050000000 -0.000070000000 + O 0.000000000000 -0.165830000000 1.454680000000 + O 0.000000000000 1.342770000000 -0.583720000000 + O 0.000000000000 -1.176830000000 -0.871010000000 + H 0.000020000000 0.501280000000 2.159930000000 + H 0.000560000000 1.618690000000 -1.514480000000 + H -0.000440000000 -2.120790000000 -0.644130000000 + N -2.649800000000 -1.445690000000 0.711420000000 + H -2.186960000000 -2.181980000000 1.244400000000 + H -3.053960000000 -1.844200000000 -0.136070000000 + H -3.367270000000 -1.005120000000 1.287210000000 + N -2.649800000000 1.339020000000 0.896300000000 + N -2.649800000000 0.106770000000 -1.607770000000 + H -3.367270000000 -0.612160000000 -1.514110000000 + H -3.053960000000 0.804320000000 1.665160000000 + N 2.649800000000 -1.445680000000 0.711420000000 + N 2.649790000000 1.339030000000 0.896300000000 + N 2.649800000000 0.106780000000 -1.607770000000 + H -2.186970000000 2.168730000000 1.267450000000 + H -3.367270000000 1.617370000000 0.226860000000 + H -2.186960000000 0.013340000000 -2.511900000000 + H -3.053970000000 1.039980000000 -1.529140000000 + H 2.186960000000 -2.181970000000 1.244400000000 + H 3.053960000000 -1.844190000000 -0.136080000000 + H 3.367270000000 -1.005100000000 1.287200000000 + H 2.186950000000 2.168740000000 1.267450000000 + H 3.053960000000 0.804330000000 1.665160000000 + H 3.367260000000 1.617380000000 0.226850000000 + H 2.186960000000 0.013350000000 -2.511900000000 + H 3.053960000000 1.039990000000 -1.529140000000 + H 3.367270000000 -0.612150000000 -1.514110000000''' + basis = {'C': 'sto-3g','H': 'sto-3g','O': 'sto-3g','N': 'sto-3g','Cr': 'cc-pvdz'} + mol = gto.M (atom=xyz, spin=6, charge=3, basis=basis, + verbose=0, output='/dev/null') + mf = scf.ROHF(mol) + mf.chkfile = 'test_lassis_targets_slow.kremer_cr2_model.chk' + mf.kernel () + las = LASSCF(mf,(6,6),((3,0),(0,3)),spin_sub=(4,4)) + try: + mo_coeff = chkfile.load (las.chkfile, 'las')['mo_coeff'] + except (OSError, TypeError, KeyError) as e: + ncas_avas, nelecas_avas, mo_coeff = avas.kernel (mf, ['Cr 3d', 'Cr 4d'], minao=mol.basis) + mc_avas = mcscf.CASCI (mf, ncas_avas, nelecas_avas) + mo_list = mc_avas.ncore + np.array ([5,6,7,8,9,10,15,16,17,18,19,20]) + mo_coeff = las.sort_mo (mo_list, mo_coeff) + mo_coeff = las.localize_init_guess (([0],[1]), mo_coeff) + las = lassi.states.spin_shuffle (las) # generate direct-exchange states + las.weights = [1.0/las.nroots,]*las.nroots # set equal weights + nroots_ref = las.nroots + las.kernel (mo_coeff) # optimize orbitals + mo_coeff = las.mo_coeff + las = LASSCF(mf,(6,6),((3,0),(0,3)),spin_sub=(4,4)) + las.lasci_(mo_coeff) + lsi = lassi.LASSIS (las).run () + with self.subTest('convergence'): + self.assertTrue (lsi.converged) + self.assertAlmostEqual (yamaguchi (lsi.e_roots, lsi.s2, 6), -12.45, 2) + + def test_alfefe (self): + xyz='''O -2.2201982441 0.3991903003 1.6944716989 + O -1.6855532446 -1.7823063217 1.4313995072 + C -2.2685178651 -0.8319550379 1.983951274 + H -2.9133420167 -1.0767285885 2.8437868053 + O 1.3882880809 2.0795561562 -1.3470856753 + O -0.7599088595 2.5809236347 -0.849227704 + C 0.3465674686 2.753895325 -1.438835699 + H 0.3723796145 3.6176996598 -2.1230911503 + O 1.0469609081 -2.6425342 0.9613246722 + O -0.3991903003 2.2201982441 1.6944716989 + O 1.7823063217 1.6855532446 1.4313995072 + C 0.8319550379 2.2685178651 1.983951274 + H 1.0767285885 2.9133420167 2.8437868053 + O -2.0795561562 -1.3882880809 -1.3470856753 + O -2.5809236347 0.7599088595 -0.849227704 + C -2.753895325 -0.3465674686 -1.438835699 + H -3.6176996598 -0.3723796145 -2.1230911503 + Fe -0.3798741893 -1.7926033629 -0.2003377303 + O 0.6422231803 -2.237796553 -1.8929145176 + Fe 1.7926033629 0.3798741893 -0.2003377303 + O 2.237796553 -0.6422231803 -1.8929145176 + O 2.6425342 -1.0469609081 0.9613246722 + Al -1.2362716421 1.2362716421 0.350650148 + O -0.0449501954 0.0449501954 0.0127621413 + C 1.645528685 -1.645528685 -2.3571124654 + H 2.0569062984 -2.0569062984 -3.2945712916 + C 2.1609242811 -2.1609242811 1.2775841868 + H 2.7960805433 -2.7960805433 1.9182443024''' + basis = {'C': 'sto-3g','H': 'sto-3g','O': 'sto-3g','Al': 'cc-pvdz','Fe': 'cc-pvdz'} + mol = gto.M (atom=xyz, spin=9, charge=0, basis=basis, max_memory=10000, verbose=0, + output='/dev/null') + mf = scf.ROHF(mol) + mf.init_guess='chk' + mf.chkfile='test_lassis_targets_slow.alfefe.chk' + mf = mf.density_fit() + mf.max_cycle=100 + mf.kernel() + las = LASSCF (mf, (5,5), ((5,1),(5,0)), spin_sub=(5,6)) + try: + las.load_chk_() + las.kernel () + except (OSError, TypeError, KeyError) as e: + ncas, nelecas, mo_coeff = avas.kernel (mf, ['Fe 3d'], openshell_option=3) + mo_coeff = las.localize_init_guess (([17],[19]), mo_coeff) + las.kernel (mo_coeff) + with self.subTest ('LASSCF convergence'): + self.assertTrue (las.converged) + with self.subTest ('same LASSCF result'): + self.assertAlmostEqual (las.e_tot, -3955.98148841553, 6) + mf.chkfile = None # prevent the spins flips down there from messing things up + las2 = LASSCF (mf, (5,5), ((1,5),(5,0)), spin_sub=(5,6)) + las2.lasci_(las.mo_coeff) + lsi = lassi.LASSIS (las2).run () + with self.subTest('LASSI convergence'): + self.assertTrue (lsi.converged) + self.assertAlmostEqual (yamaguchi (lsi.e_roots, lsi.s2, 9), -4.40, 2) + + +if __name__ == "__main__": + print("Full Tests for SA-LASSI of c2h4n4 molecule") + unittest.main() + diff --git a/tests/mcpdft/test_lassipdft.py b/tests/mcpdft/test_lassipdft.py index f6ee8f94..4ecfa31c 100644 --- a/tests/mcpdft/test_lassipdft.py +++ b/tests/mcpdft/test_lassipdft.py @@ -1,36 +1,48 @@ -import sys -from pyscf import gto, scf, tools, dft, lib -from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF -from mrh.my_dmet import localintegrals, dmet, fragments -from mrh.my_pyscf import mcpdft import unittest +from pyscf import lib, gto, scf +from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF +from mrh.my_pyscf.lassi import LASSI +from mrh.my_pyscf.lassi.states import all_single_excitations +from mrh.my_pyscf.mcscf.lasci import get_space_info class KnownValues(unittest.TestCase): def test_ethene (self): - mol = gto.M() - mol.atom = '''C -0.662958 0.000000 0.000000; - C 0.662958 0.000000 0.000000; - H -1.256559 -0.924026 0.000000; - H 1.256559 -0.924026 0.000000; - H -1.256559 0.924026 0.000000; - H 1.256559 0.924026 0.000000''' - mol.basis='sto3g' - mol.verbose = 0 - mol.output = '/dev/null' - mol.build() - mf = scf.ROHF(mol).newton().run() - mc = mcpdft.LASSI(mf, 'tPBE', (2, ), (2, ), grid_level=1).state_average( - [1, 0], spins=[[0,], [2, ]], smults=[[1, ], [3, ]], charges=[[0, ],[0, ]]) - mo0 = mc.localize_init_guess(([0, 1],), mc.sort_mo([8, 9])) - mc.kernel(mo0) - elassi = mc.e_mcscf[0] - epdft = mc.e_tot[0] - self.assertAlmostEqual (elassi , -77.1154672717181, 7) # Reference values of CASSCF and CAS-PDFT - self.assertAlmostEqual (epdft , -77.49805221093968, 7) + xyz='''H 0 0 0 + H 1 0 0 + H 3 0 0 + H 4 0 0''' + + mol = gto.M (atom=xyz, basis='sto3g', symmetry=False, verbose=0, output='/dev/null') + mf = scf.RHF (mol).run () + + # LASSCF and LASSI + las = LASSCF (mf, (2,2), (2,2), spin_sub=(1,1)) + las.lasci () + las1 = las + for i in range (2): las1 = all_single_excitations (las1) + charges, spins, smults, wfnsyms = get_space_info (las1) + lroots = 4 - smults + idx = (charges!=0) & (lroots==3) + lroots[idx] = 1 + las1.conv_tol_grad = las.conv_tol_self = 9e99 + las1.lasci (lroots=lroots.T) + las1.dump_spaces () + lsi = LASSI (las1) + lsi.kernel (opt=0) + + from mrh.my_pyscf import mcpdft + lsipdft = mcpdft.LASSI(lsi, 'tPBE', (2,2), (2,2)) + lsipdft.kernel() + + # CASCI limit + from pyscf import mcpdft + mc = mcpdft.CASCI (mf, 'tPBE', 4, 4).run () + + self.assertAlmostEqual (lsi.e_roots[0], mc.e_mcscf, 7) + self.assertAlmostEqual (lsipdft.e_tot[0], mc.e_tot, 7) if __name__ == "__main__": print("Full Tests for LASSI-PDFT") unittest.main() - diff --git a/tests/nacs/test_cmspdft_nacs.py b/tests/nacs/test_cmspdft_nacs.py index 7390de86..24149051 100644 --- a/tests/nacs/test_cmspdft_nacs.py +++ b/tests/nacs/test_cmspdft_nacs.py @@ -1,9 +1,6 @@ import numpy as np -from pyscf import gto, scf, df, mcscf, lib, mcpdft -from mrh.my_pyscf.fci import csf_solver -from mrh.my_pyscf.tools.molcas2pyscf import * +from pyscf import gto, scf, mcpdft from mrh.my_pyscf.dft.openmolcas_grids import quasi_ultrafine -from mrh.my_pyscf.grad.mspdft_nacs import NonAdiabaticCouplings import unittest @@ -46,9 +43,7 @@ def diatomic(atom1, atom2, r, basis, ncas, nelecas, nstates, mo = mc.sort_mo_by_irrep(cas_irrep) mc.kernel(mo) - - - return NonAdiabaticCouplings(mc) + return mc.nac_method() def tearDownModule(): diff --git a/tests/nacs/test_sacasscf_nacs.py b/tests/nacs/test_sacasscf_nacs.py index 32d00ff5..cdd44c76 100644 --- a/tests/nacs/test_sacasscf_nacs.py +++ b/tests/nacs/test_sacasscf_nacs.py @@ -1,6 +1,5 @@ import numpy as np -from pyscf import gto, scf, df, mcscf, lib -from mrh.my_pyscf.grad.sacasscf_nacs import NonAdiabaticCouplings +from pyscf import gto, scf, mcscf import unittest @@ -31,7 +30,7 @@ def diatomic(atom1, atom2, r, basis, ncas, nelecas, nstates, mc.kernel(mo) - return NonAdiabaticCouplings(mc) + return mc.nac_method() def tearDownModule(): diff --git a/util/la.py b/util/la.py index f15eedad..0e58bb90 100644 --- a/util/la.py +++ b/util/la.py @@ -73,38 +73,47 @@ def safe_svd (a, full_matrices=True): except scipy.linalg.LinAlgError as err: warn ("Encountered scipy.linalg.LinAlgError during SVD: {}".format (err)) warn ("attempting fake SVD using eigh") - aTa = a.conj ().T @ a - evals_aTa, v = scipy.linalg.eigh (-aTa) - N = len (evals_aTa) - aaT = a @ a.conj ().T - evals_aaT, u = scipy.linalg.eigh (-aaT) - M = len (evals_aaT) - K = min (M,N) - svals = np.sqrt (np.maximum ((evals_aaT[:K] + evals_aTa[:K])*-.5, 0)) - # Try to sort out degenerate-sval and sign issues by a second pass to - # linalg.svd with the null space removed - try: - svals_rel = svals/np.amax(svals) - K1 = np.count_nonzero (svals_rel>1e-7) - a1 = u[:,:K1].conj ().T @ a @ v[:,:K1] - u1, svals1, v1h = scipy.linalg.svd (a1, full_matrices=False) - svals[:K1] = svals1 - u[:,:K1] = u[:,:K1] @ u1 - v[:,:K1] = v[:,:K1] @ v1h.conj ().T - except scipy.linalg.LinAlgError as err1: - warn ("Second pass to scipy.linalg also failed: {}".format (err1)) - warn ("Degenerate manifolds are mixed randomly!") - # Assign sign flips to u arbitrarily! - idx_signflip = ((a @ v[:,:K]) * u[:,:K].conj ()).sum (0) < 0 - u[:,:K][:,idx_signflip] = -u[:,:K][:,idx_signflip] - vH = v.conj ().T - a_test = (u[:,:K] * svals[None,:]) @ vH[:K] - a_err = np.amax (np.abs (a_test-a)) - warn ("Checking fake svd accuracy") - warn ("max(abs(u.s.vh - a)) = %e", a_err) - if not full_matrices: - u, v = u[:,:K], v[:,:K] - return u, svals, vH + return fake_svd (a, full_matrices=full_matrices) + except ValueError as err: + warn ("Encountered ValueError during SVD: {}".format (err)) + warn ("Matrix contains {} infs and {} nans".format ( + np.count_nonzero (np.isinf (a)), np.count_nonzero (np.isnan (a)))) + warn ("attempting fake SVD using eigh") + return fake_svd (a, full_matrices=full_matrices) + + def fake_svd (a, full_matrices=True): + aTa = a.conj ().T @ a + evals_aTa, v = scipy.linalg.eigh (-aTa) + N = len (evals_aTa) + aaT = a @ a.conj ().T + evals_aaT, u = scipy.linalg.eigh (-aaT) + M = len (evals_aaT) + K = min (M,N) + svals = np.sqrt (np.maximum ((evals_aaT[:K] + evals_aTa[:K])*-.5, 0)) + # Try to sort out degenerate-sval and sign issues by a second pass to + # linalg.svd with the null space removed + try: + svals_rel = svals/np.amax(svals) + K1 = np.count_nonzero (svals_rel>1e-7) + a1 = u[:,:K1].conj ().T @ a @ v[:,:K1] + u1, svals1, v1h = scipy.linalg.svd (a1, full_matrices=False) + svals[:K1] = svals1 + u[:,:K1] = u[:,:K1] @ u1 + v[:,:K1] = v[:,:K1] @ v1h.conj ().T + except scipy.linalg.LinAlgError as err1: + warn ("Second pass to scipy.linalg also failed: {}".format (err1)) + warn ("Degenerate manifolds are mixed randomly!") + # Assign sign flips to u arbitrarily! + idx_signflip = ((a @ v[:,:K]) * u[:,:K].conj ()).sum (0) < 0 + u[:,:K][:,idx_signflip] = -u[:,:K][:,idx_signflip] + vH = v.conj ().T + a_test = (u[:,:K] * svals[None,:]) @ vH[:K] + a_err = np.amax (np.abs (a_test-a)) + warn ("Checking fake svd accuracy") + warn ("max(abs(u.s.vh - a)) = %e", a_err) + if not full_matrices: + u, v = u[:,:K], v[:,:K] + return u, svals, vH return safe_svd safe_svd = safe_svd_warner ()