Source code for pyscf.fci.selected_ci

#!/usr/bin/env python
# Copyright 2014-2021 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

'''
Selected CI.

This is an inefficient dialect of Selected CI using the same structure as
determinant based FCI algorithm. For the efficient Selected CI programs,
Dice program (https://github.com/sanshar/Dice.git) is a good candidate.

Simple usage::

    >>> from pyscf import gto, scf, ao2mo, fci
    >>> mol = gto.M(atom='C 0 0 0; C 0 0 1')
    >>> mf = scf.RHF(mol).run()
    >>> h1 = mf.mo_coeff.T.dot(mf.get_hcore()).dot(mf.mo_coeff)
    >>> h2 = ao2mo.kernel(mol, mf.mo_coeff)
    >>> e = fci.selected_ci.kernel(h1, h2, mf.mo_coeff.shape[1], mol.nelectron)[0]
'''

import ctypes
import numpy
from pyscf import lib
from pyscf.lib import logger
from pyscf import ao2mo
from pyscf.fci import cistring
from pyscf.fci import direct_spin1
from pyscf.fci import rdm
from pyscf import __config__

libfci = direct_spin1.libfci

[docs] @lib.with_doc(direct_spin1.contract_2e.__doc__) def contract_2e(eri, civec_strs, norb, nelec, link_index=None): ci_coeff, nelec, ci_strs = _unpack(civec_strs, nelec) if link_index is None: link_index = _all_linkstr_index(ci_strs, norb, nelec) cd_indexa, dd_indexa, cd_indexb, dd_indexb = link_index na, nlinka = cd_indexa.shape[:2] nb, nlinkb = cd_indexb.shape[:2] eri = ao2mo.restore(1, eri, norb) eri1 = eri.transpose(0,2,1,3) - eri.transpose(0,2,3,1) idx,idy = numpy.tril_indices(norb, -1) idx = idx * norb + idy eri1 = lib.take_2d(eri1.reshape(norb**2,-1), idx, idx) * 2 fcivec = ci_coeff.reshape(na,nb) # (bb|bb) if nelec[1] > 1: mb, mlinkb = dd_indexb.shape[:2] fcivecT = lib.transpose(fcivec) ci1T = numpy.zeros((nb,na)) libfci.SCIcontract_2e_aaaa(eri1.ctypes.data_as(ctypes.c_void_p), fcivecT.ctypes.data_as(ctypes.c_void_p), ci1T.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(nb), ctypes.c_int(na), ctypes.c_int(mb), ctypes.c_int(mlinkb), dd_indexb.ctypes.data_as(ctypes.c_void_p)) ci1 = lib.transpose(ci1T, out=fcivecT) else: ci1 = numpy.zeros_like(fcivec) # (aa|aa) if nelec[0] > 1: ma, mlinka = dd_indexa.shape[:2] libfci.SCIcontract_2e_aaaa(eri1.ctypes.data_as(ctypes.c_void_p), fcivec.ctypes.data_as(ctypes.c_void_p), ci1.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(na), ctypes.c_int(nb), ctypes.c_int(ma), ctypes.c_int(mlinka), dd_indexa.ctypes.data_as(ctypes.c_void_p)) # Adding h_ps below to because contract_2e function computes the # contraction "E_{pq}E_{rs} V_{pqrs} |CI>" (~ p^+ q r^+ s |CI>) while # the actual contraction for (aa|aa) and (bb|bb) part is # "p^+ r^+ s q V_{pqrs} |CI>". To make (aa|aa) and (bb|bb) code reproduce # "p^+ q r^+ s |CI>", we employ the identity # p^+ q r^+ s = p^+ r^+ s q + delta(qr) p^+ s # the second term is the source of h_ps h_ps = numpy.einsum('pqqs->ps', eri) eri1 = eri * 2 for k in range(norb): eri1[:,:,k,k] += h_ps/nelec[0] eri1[k,k,:,:] += h_ps/nelec[1] eri1 = ao2mo.restore(4, eri1, norb) # (bb|aa) libfci.SCIcontract_2e_bbaa(eri1.ctypes.data_as(ctypes.c_void_p), fcivec.ctypes.data_as(ctypes.c_void_p), ci1.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(na), ctypes.c_int(nb), ctypes.c_int(nlinka), ctypes.c_int(nlinkb), cd_indexa.ctypes.data_as(ctypes.c_void_p), cd_indexb.ctypes.data_as(ctypes.c_void_p)) return _as_SCIvector(ci1.reshape(ci_coeff.shape), ci_strs)
[docs] def select_strs(myci, eri, eri_pq_max, civec_max, strs, norb, nelec): strs = numpy.asarray(strs, dtype=numpy.int64) nstrs = len(strs) nvir = norb - nelec strs_add = numpy.empty((nstrs*((nelec+1)*(nvir+1))**2//4), dtype=numpy.int64) libfci.SCIselect_strs.restype = ctypes.c_int nadd = libfci.SCIselect_strs(strs_add.ctypes.data_as(ctypes.c_void_p), strs.ctypes.data_as(ctypes.c_void_p), eri.ctypes.data_as(ctypes.c_void_p), eri_pq_max.ctypes.data_as(ctypes.c_void_p), civec_max.ctypes.data_as(ctypes.c_void_p), ctypes.c_double(myci.select_cutoff), ctypes.c_int(norb), ctypes.c_int(nelec), ctypes.c_int(nstrs)) strs_add = sorted(set(strs_add[:nadd]) - set(strs)) return numpy.asarray(strs_add, dtype=numpy.int64)
[docs] def enlarge_space(myci, civec_strs, eri, norb, nelec): if isinstance(civec_strs, (tuple, list)): nelec, (strsa, strsb) = _unpack(civec_strs[0], nelec)[1:] ci_coeff = lib.asarray(civec_strs) else: ci_coeff, nelec, (strsa, strsb) = _unpack(civec_strs, nelec) na = len(strsa) nb = len(strsb) ci0 = ci_coeff.reshape(-1,na,nb) civec_a_max = lib.norm(ci0, axis=2).max(axis=0) civec_b_max = lib.norm(ci0, axis=1).max(axis=0) ci_aidx = numpy.where(civec_a_max > myci.ci_coeff_cutoff)[0] ci_bidx = numpy.where(civec_b_max > myci.ci_coeff_cutoff)[0] civec_a_max = civec_a_max[ci_aidx] civec_b_max = civec_b_max[ci_bidx] strsa = strsa[ci_aidx] strsb = strsb[ci_bidx] eri = ao2mo.restore(1, eri, norb) eri_pq_max = abs(eri.reshape(norb**2,-1)).max(axis=1).reshape(norb,norb) strsa_add = select_strs(myci, eri, eri_pq_max, civec_a_max, strsa, norb, nelec[0]) strsb_add = select_strs(myci, eri, eri_pq_max, civec_b_max, strsb, norb, nelec[1]) strsa = numpy.append(strsa, strsa_add) strsb = numpy.append(strsb, strsb_add) aidx = numpy.argsort(strsa) bidx = numpy.argsort(strsb) ci_strs = (strsa[aidx], strsb[bidx]) aidx = numpy.where(aidx < len(ci_aidx))[0] bidx = numpy.where(bidx < len(ci_bidx))[0] ma = len(strsa) mb = len(strsb) cs = [] for i in range(ci0.shape[0]): ci1 = numpy.zeros((ma,mb)) tmp = lib.take_2d(ci0[i], ci_aidx, ci_bidx) lib.takebak_2d(ci1, tmp, aidx, bidx) cs.append(_as_SCIvector(ci1, ci_strs)) if not isinstance(civec_strs, (tuple, list)) and civec_strs.ndim < 3: cs = cs[0] return cs
[docs] def cre_des_linkstr(strs, norb, nelec, tril=False): '''Given intermediates, the link table to generate input strs ''' strs = numpy.asarray(strs, dtype=numpy.int64) nvir = norb - nelec nstrs = len(strs) link_index = numpy.zeros((nstrs,nelec+nelec*nvir,4), dtype=numpy.int32) libfci.SCIcre_des_linkstr(link_index.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(nstrs), ctypes.c_int(nelec), strs.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(tril)) return link_index
[docs] def cre_des_linkstr_tril(strs, norb, nelec): '''Given intermediates, the link table to generate input strs ''' return cre_des_linkstr(strs, norb, nelec, True)
[docs] def des_des_linkstr(strs, norb, nelec, tril=False): '''Given intermediates, the link table to generate input strs ''' if nelec < 2: return None strs = numpy.asarray(strs, dtype=numpy.int64) nvir = norb - nelec nstrs = len(strs) inter1 = numpy.empty((nstrs*nelec), dtype=numpy.int64) libfci.SCIdes_uniq_strs.restype = ctypes.c_int ninter = libfci.SCIdes_uniq_strs(inter1.ctypes.data_as(ctypes.c_void_p), strs.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(nelec), ctypes.c_int(nstrs)) inter1 = numpy.asarray(sorted(set(inter1[:ninter])), dtype=numpy.int64) ninter = len(inter1) inter = numpy.empty((ninter*nelec), dtype=numpy.int64) ninter = libfci.SCIdes_uniq_strs(inter.ctypes.data_as(ctypes.c_void_p), inter1.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(nelec-1), ctypes.c_int(ninter)) inter = numpy.asarray(sorted(set(inter[:ninter])), dtype=numpy.int64) ninter = len(inter) nvir += 2 link_index = numpy.zeros((ninter,nvir*nvir,4), dtype=numpy.int32) libfci.SCIdes_des_linkstr(link_index.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(nelec), ctypes.c_int(nstrs), ctypes.c_int(ninter), strs.ctypes.data_as(ctypes.c_void_p), inter.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(tril)) return link_index
[docs] def des_des_linkstr_tril(strs, norb, nelec): '''Given intermediates, the link table to generate input strs ''' return des_des_linkstr(strs, norb, nelec, True)
[docs] def gen_des_linkstr(strs, norb, nelec): '''Given intermediates, the link table to generate input strs ''' if nelec < 1: return None strs = numpy.asarray(strs, dtype=numpy.int64) nvir = norb - nelec nstrs = len(strs) inter = numpy.empty((nstrs*nelec), dtype=numpy.int64) libfci.SCIdes_uniq_strs.restype = ctypes.c_int ninter = libfci.SCIdes_uniq_strs(inter.ctypes.data_as(ctypes.c_void_p), strs.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(nelec), ctypes.c_int(nstrs)) inter = numpy.asarray(sorted(set(inter[:ninter])), dtype=numpy.int64) ninter = len(inter) nvir += 1 link_index = numpy.zeros((ninter,nvir,4), dtype=numpy.int32) libfci.SCIdes_linkstr(link_index.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(nelec), ctypes.c_int(nstrs), ctypes.c_int(ninter), strs.ctypes.data_as(ctypes.c_void_p), inter.ctypes.data_as(ctypes.c_void_p)) return link_index
[docs] def gen_cre_linkstr(strs, norb, nelec): '''Given intermediates, the link table to generate input strs ''' if nelec == norb: return None strs = numpy.asarray(strs, dtype=numpy.int64) nvir = norb - nelec nstrs = len(strs) inter = numpy.empty((nstrs*nvir), dtype=numpy.int64) libfci.SCIcre_uniq_strs.restype = ctypes.c_int ninter = libfci.SCIcre_uniq_strs(inter.ctypes.data_as(ctypes.c_void_p), strs.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(nelec), ctypes.c_int(nstrs)) inter = numpy.asarray(sorted(set(inter[:ninter])), dtype=numpy.int64) ninter = len(inter) link_index = numpy.zeros((ninter,nelec+1,4), dtype=numpy.int32) libfci.SCIcre_linkstr(link_index.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(nelec), ctypes.c_int(nstrs), ctypes.c_int(ninter), strs.ctypes.data_as(ctypes.c_void_p), inter.ctypes.data_as(ctypes.c_void_p)) return link_index
[docs] def make_hdiag(h1e, eri, ci_strs, norb, nelec, compress=False): ci_coeff, nelec, ci_strs = _unpack(None, nelec, ci_strs) na = len(ci_strs[0]) nb = len(ci_strs[1]) hdiag = numpy.empty(na*nb) h1e = numpy.asarray(h1e, order='C') eri = ao2mo.restore(1, eri, norb) jdiag = numpy.asarray(numpy.einsum('iijj->ij',eri), order='C') kdiag = numpy.asarray(numpy.einsum('ijji->ij',eri), order='C') c_h1e = h1e.ctypes.data_as(ctypes.c_void_p) c_jdiag = jdiag.ctypes.data_as(ctypes.c_void_p) c_kdiag = kdiag.ctypes.data_as(ctypes.c_void_p) occslsta = cistring._strs2occslst(ci_strs[0], norb) occslstb = cistring._strs2occslst(ci_strs[1], norb) libfci.FCImake_hdiag_uhf(hdiag.ctypes.data_as(ctypes.c_void_p), c_h1e, c_h1e, c_jdiag, c_jdiag, c_jdiag, c_kdiag, c_kdiag, ctypes.c_int(norb), ctypes.c_int(na), ctypes.c_int(nb), ctypes.c_int(nelec[0]), ctypes.c_int(nelec[1]), occslsta.ctypes.data_as(ctypes.c_void_p), occslstb.ctypes.data_as(ctypes.c_void_p)) return hdiag
[docs] def kernel_fixed_space(myci, h1e, eri, norb, nelec, ci_strs, ci0=None, tol=None, lindep=None, max_cycle=None, max_space=None, nroots=None, davidson_only=None, max_memory=None, verbose=None, ecore=0, **kwargs): log = logger.new_logger(myci, verbose) if tol is None: tol = myci.conv_tol if lindep is None: lindep = myci.lindep if max_cycle is None: max_cycle = myci.max_cycle if max_space is None: max_space = myci.max_space if max_memory is None: max_memory = myci.max_memory if nroots is None: nroots = myci.nroots if myci.verbose >= logger.WARN: myci.check_sanity() nelec = direct_spin1._unpack_nelec(nelec, myci.spin) ci0, nelec, ci_strs = _unpack(ci0, nelec, ci_strs) na = len(ci_strs[0]) nb = len(ci_strs[1]) h2e = direct_spin1.absorb_h1e(h1e, eri, norb, nelec, .5) h2e = ao2mo.restore(1, h2e, norb) link_index = _all_linkstr_index(ci_strs, norb, nelec) hdiag = myci.make_hdiag(h1e, eri, ci_strs, norb, nelec, compress=True) if isinstance(ci0, SCIvector): if ci0.size == na*nb: ci0 = [ci0.ravel()] else: ci0 = [x.ravel() for x in ci0] else: ci0 = myci.get_init_guess(ci_strs, norb, nelec, nroots, hdiag) cpu0 = [logger.process_clock(), logger.perf_counter()] def hop(c): hc = myci.contract_2e(h2e, _as_SCIvector(c, ci_strs), norb, nelec, link_index) cpu0[:] = log.timer_debug1('contract_2e', *cpu0) return hc.reshape(-1) precond = lambda x, e, *args: x/(hdiag-e+1e-4) #e, c = lib.davidson(hop, ci0, precond, tol=myci.conv_tol) e, c = myci.eig(hop, ci0, precond, tol=tol, lindep=lindep, max_cycle=max_cycle, max_space=max_space, nroots=nroots, max_memory=max_memory, verbose=log, **kwargs) if nroots > 1: return e+ecore, [_as_SCIvector(ci.reshape(na,nb),ci_strs) for ci in c] else: return e+ecore, _as_SCIvector(c.reshape(na,nb), ci_strs)
[docs] def kernel_float_space(myci, h1e, eri, norb, nelec, ci0=None, tol=None, lindep=None, max_cycle=None, max_space=None, nroots=None, davidson_only=None, max_memory=None, verbose=None, ecore=0, **kwargs): log = logger.new_logger(myci, verbose) if tol is None: tol = myci.conv_tol if lindep is None: lindep = myci.lindep if max_cycle is None: max_cycle = myci.max_cycle if max_space is None: max_space = myci.max_space if max_memory is None: max_memory = myci.max_memory if nroots is None: nroots = myci.nroots if myci.verbose >= logger.WARN: myci.check_sanity() nelec = direct_spin1._unpack_nelec(nelec, myci.spin) h2e = direct_spin1.absorb_h1e(h1e, eri, norb, nelec, .5) h2e = ao2mo.restore(1, h2e, norb) # TODO: initial guess from CISD if isinstance(ci0, SCIvector): if ci0.size == len(ci0._strs[0])*len(ci0._strs[1]): ci0 = [ci0.ravel()] else: ci0 = [x.ravel() for x in ci0] else: ci_strs = (numpy.asarray([int('1'*nelec[0], 2)]), numpy.asarray([int('1'*nelec[1], 2)])) ci0 = _as_SCIvector(numpy.ones((1,1)), ci_strs) ci0 = myci.enlarge_space(ci0, h2e, norb, nelec) if ci0.size < nroots: log.warn(''' Selected-CI space generated from HF ground state (by double exciting) is not enough for excited states. HOMO->LUMO excitations are included in the initial guess. NOTE: This may introduce excited states of different symmetry.\n''') corea = '1' * (nelec[0]-1) coreb = '1' * (nelec[1]-1) ci_strs = (numpy.asarray([int('1'+corea, 2), int('10'+corea, 2)]), numpy.asarray([int('1'+coreb, 2), int('10'+coreb, 2)])) ci0 = _as_SCIvector(numpy.ones((2,2)), ci_strs) ci0 = myci.enlarge_space(ci0, h2e, norb, nelec) if ci0.size < nroots: raise RuntimeError('Not enough selected-CI space for %d states' % nroots) ci_strs = ci0._strs hdiag = myci.make_hdiag(h1e, eri, ci_strs, norb, nelec, compress=True) ci0 = myci.get_init_guess(ci_strs, norb, nelec, nroots, hdiag) def hop(c): hc = myci.contract_2e(h2e, _as_SCIvector(c, ci_strs), norb, nelec, link_index) return hc.ravel() precond = lambda x, e, *args: x/(hdiag-e+myci.level_shift) namax = cistring.num_strings(norb, nelec[0]) nbmax = cistring.num_strings(norb, nelec[1]) e_last = 0 float_tol = myci.start_tol tol_decay_rate = myci.tol_decay_rate # conv = False for icycle in range(norb): ci_strs = ci0[0]._strs float_tol = max(float_tol*tol_decay_rate, tol*1e2) log.debug('cycle %d ci.shape %s float_tol %g', icycle, (len(ci_strs[0]), len(ci_strs[1])), float_tol) ci0 = [c.ravel() for c in ci0] link_index = _all_linkstr_index(ci_strs, norb, nelec) hdiag = myci.make_hdiag(h1e, eri, ci_strs, norb, nelec, compress=True) #e, ci0 = lib.davidson(hop, ci0.reshape(-1), precond, tol=float_tol) e, ci0 = myci.eig(hop, ci0, precond, tol=float_tol, lindep=lindep, max_cycle=max_cycle, max_space=max_space, nroots=nroots, max_memory=max_memory, verbose=log, **kwargs) if nroots > 1: ci0 = [_as_SCIvector(c, ci_strs) for c in ci0] de, e_last = min(e)-e_last, min(e) log.info('cycle %d E = %s dE = %.8g', icycle, e+ecore, de) else: ci0 = [_as_SCIvector(ci0, ci_strs)] de, e_last = e-e_last, e log.info('cycle %d E = %.15g dE = %.8g', icycle, e+ecore, de) if ci0[0].shape == (namax,nbmax) or abs(de) < tol*1e3: # conv = True break last_ci0_size = float(len(ci_strs[0])), float(len(ci_strs[1])) ci0 = myci.enlarge_space(ci0, h2e, norb, nelec) na = len(ci0[0]._strs[0]) nb = len(ci0[0]._strs[1]) if ((.99 < na/last_ci0_size[0] < 1.01) and (.99 < nb/last_ci0_size[1] < 1.01)): # conv = True break ci_strs = ci0[0]._strs log.debug('Extra CI in selected space %s', (len(ci_strs[0]), len(ci_strs[1]))) ci0 = [c.ravel() for c in ci0] link_index = _all_linkstr_index(ci_strs, norb, nelec) hdiag = myci.make_hdiag(h1e, eri, ci_strs, norb, nelec, compress=True) e, c = myci.eig(hop, ci0, precond, tol=tol, lindep=lindep, max_cycle=max_cycle, max_space=max_space, nroots=nroots, max_memory=max_memory, verbose=log, **kwargs) na = len(ci_strs[0]) nb = len(ci_strs[1]) if nroots > 1: for i, ei in enumerate(e+ecore): log.info('Selected CI state %d E = %.15g', i, ei) return e+ecore, [_as_SCIvector(ci.reshape(na,nb),ci_strs) for ci in c] else: log.info('Selected CI E = %.15g', e+ecore) return e+ecore, _as_SCIvector(c.reshape(na,nb), ci_strs)
[docs] def kernel(h1e, eri, norb, nelec, ci0=None, level_shift=1e-3, tol=1e-10, lindep=1e-14, max_cycle=50, max_space=12, nroots=1, davidson_only=False, pspace_size=400, orbsym=None, wfnsym=None, select_cutoff=1e-3, ci_coeff_cutoff=1e-3, ecore=0, **kwargs): return direct_spin1._kfactory(SelectedCI, h1e, eri, norb, nelec, ci0, level_shift, tol, lindep, max_cycle, max_space, nroots, davidson_only, pspace_size, select_cutoff=select_cutoff, ci_coeff_cutoff=ci_coeff_cutoff, ecore=ecore, **kwargs)
[docs] def make_rdm1s(civec_strs, norb, nelec, link_index=None): r'''Spin separated 1-particle density matrices. The return values include two density matrices: (alpha,alpha), (beta,beta) dm1[p,q] = <q^\dagger p> The convention is based on McWeeney's book, Eq (5.4.20). The contraction between 1-particle Hamiltonian and rdm1 is E = einsum('pq,qp', h1, rdm1) ''' ci_coeff, nelec, ci_strs = _unpack(civec_strs, nelec) if link_index is None: cd_indexa = cre_des_linkstr(ci_strs[0], norb, nelec[0]) cd_indexb = cre_des_linkstr(ci_strs[1], norb, nelec[1]) else: cd_indexa, dd_indexa, cd_indexb, dd_indexb = link_index rdm1a = rdm.make_rdm1_spin1('FCImake_rdm1a', ci_coeff, ci_coeff, norb, nelec, (cd_indexa,cd_indexb)) rdm1b = rdm.make_rdm1_spin1('FCImake_rdm1b', ci_coeff, ci_coeff, norb, nelec, (cd_indexa,cd_indexb)) return rdm1a, rdm1b
[docs] def make_rdm1(civec_strs, norb, nelec, link_index=None): r'''Spin-traced 1-particle density matrix. dm1[p,q] = <q_alpha^\dagger p_alpha> + <q_beta^\dagger p_beta> The convention is based on McWeeney's book, Eq (5.4.20) The contraction between 1-particle Hamiltonian and rdm1 is E = einsum('pq,qp', h1, rdm1) ''' rdm1a, rdm1b = make_rdm1s(civec_strs, norb, nelec, link_index) return rdm1a + rdm1b
# dm[p,q,r,s] = <|p^+ q r^+ s|>
[docs] def make_rdm2s(civec_strs, norb, nelec, link_index=None, **kwargs): r'''Spin separated 2-particle density matrices. The return values include three density matrices: (alpha,alpha,alpha,alpha), (alpha,alpha,beta,beta), (beta,beta,beta,beta) 2pdm[p,q,r,s] = :math:`\langle p^\dagger r^\dagger s q\rangle` ''' ci_coeff, nelec, ci_strs = _unpack(civec_strs, nelec) if link_index is None: cd_indexa = cre_des_linkstr(ci_strs[0], norb, nelec[0]) dd_indexa = des_des_linkstr(ci_strs[0], norb, nelec[0]) cd_indexb = cre_des_linkstr(ci_strs[1], norb, nelec[1]) dd_indexb = des_des_linkstr(ci_strs[1], norb, nelec[1]) else: cd_indexa, dd_indexa, cd_indexb, dd_indexb = link_index na, nlinka = cd_indexa.shape[:2] nb, nlinkb = cd_indexb.shape[:2] fcivec = ci_coeff.reshape(na,nb) # (bb|aa) and (aa|bb) dm2ab = rdm.make_rdm12_spin1('FCItdm12kern_ab', fcivec, fcivec, norb, nelec, (cd_indexa,cd_indexb), 0)[1] # (aa|aa) dm2aa = numpy.zeros([norb]*4) if nelec[0] > 1: ma, mlinka = dd_indexa.shape[:2] libfci.SCIrdm2_aaaa(libfci.SCIrdm2kern_aaaa, dm2aa.ctypes.data_as(ctypes.c_void_p), fcivec.ctypes.data_as(ctypes.c_void_p), fcivec.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(na), ctypes.c_int(nb), ctypes.c_int(ma), ctypes.c_int(mlinka), dd_indexa.ctypes.data_as(ctypes.c_void_p)) # (bb|bb) dm2bb = numpy.zeros([norb]*4) if nelec[1] > 1: mb, mlinkb = dd_indexb.shape[:2] fcivecT = lib.transpose(fcivec) libfci.SCIrdm2_aaaa(libfci.SCIrdm2kern_aaaa, dm2bb.ctypes.data_as(ctypes.c_void_p), fcivecT.ctypes.data_as(ctypes.c_void_p), fcivecT.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(norb), ctypes.c_int(nb), ctypes.c_int(na), ctypes.c_int(mb), ctypes.c_int(mlinkb), dd_indexb.ctypes.data_as(ctypes.c_void_p)) return dm2aa, dm2ab, dm2bb
[docs] def make_rdm2(civec_strs, norb, nelec, link_index=None, **kwargs): r'''Spin-traced two-particle density matrix. 2pdm[p,q,r,s] = :math:`\langle p_\alpha^\dagger r_\alpha^\dagger s_\alpha q_\alpha\rangle + \langle p_\beta^\dagger r_\alpha^\dagger s_\alpha q_\beta\rangle + \langle p_\alpha^\dagger r_\beta^\dagger s_\beta q_\alpha\rangle + \langle p_\beta^\dagger r_\beta^\dagger s_\beta q_\beta\rangle`. ''' dm2aa, dm2ab, dm2bb = make_rdm2s(civec_strs, norb, nelec, link_index) dm2aa += dm2bb dm2aa += dm2ab dm2aa += dm2ab.transpose(2,3,0,1) return dm2aa
[docs] def trans_rdm1s(cibra_strs, ciket_strs, norb, nelec, link_index=None): r'''Spin separated transition 1-particle density matrices. See also function :func:`make_rdm1s` 1pdm[p,q] = :math:`\langle q^\dagger p \rangle` ''' cibra, nelec, ci_strs = _unpack(cibra_strs, nelec) ciket, nelec1, ci_strs1 = _unpack(ciket_strs, nelec) assert (all(ci_strs[0] == ci_strs1[0]) and all(ci_strs[1] == ci_strs1[1])) if link_index is None: cd_indexa = cre_des_linkstr(ci_strs[0], norb, nelec[0]) cd_indexb = cre_des_linkstr(ci_strs[1], norb, nelec[1]) else: cd_indexa, dd_indexa, cd_indexb, dd_indexb = link_index rdm1a = rdm.make_rdm1_spin1('FCItrans_rdm1a', cibra, ciket, norb, nelec, (cd_indexa,cd_indexb)) rdm1b = rdm.make_rdm1_spin1('FCItrans_rdm1b', cibra, ciket, norb, nelec, (cd_indexa,cd_indexb)) return rdm1a, rdm1b
[docs] def trans_rdm1(cibra_strs, ciket_strs, norb, nelec, link_index=None): r'''Spin traced transition 1-particle density matrices. See also function :func:`make_rdm1` 1pdm[p,q] = :math:`\langle q_\alpha^\dagger p_\alpha \rangle + \langle q_\beta^\dagger p_\beta \rangle` ''' rdm1a, rdm1b = trans_rdm1s(cibra_strs, ciket_strs, norb, nelec, link_index) return rdm1a + rdm1b
[docs] def spin_square(civec_strs, norb, nelec): '''Spin square for RHF-FCI CI wfn only (obtained from spin-degenerated Hamiltonian)''' ci1 = contract_ss(civec_strs, norb, nelec) ss = numpy.einsum('ij,ij->', civec_strs.reshape(ci1.shape), ci1) s = numpy.sqrt(ss+.25) - .5 multip = s*2+1 return ss, multip
[docs] def contract_ss(civec_strs, norb, nelec): r''' S^2 |\Psi\rangle ''' ci_coeff, nelec, ci_strs = _unpack(civec_strs, nelec) strsa, strsb = ci_strs neleca, nelecb = nelec ci_coeff = ci_coeff.reshape(len(strsa),len(strsb)) def gen_map(fstr_index, strs, nelec, des=True): a_index = fstr_index(strs, norb, nelec) amap = numpy.zeros((a_index.shape[0],norb,2), dtype=numpy.int32) if des: for k, tab in enumerate(a_index): sign = tab[:,3] tab = tab[sign!=0] amap[k,tab[:,1]] = tab[:,2:] else: for k, tab in enumerate(a_index): sign = tab[:,3] tab = tab[sign!=0] amap[k,tab[:,0]] = tab[:,2:] return amap if neleca > 0: ades = gen_map(gen_des_linkstr, strsa, neleca) else: ades = None if nelecb > 0: bdes = gen_map(gen_des_linkstr, strsb, nelecb) else: bdes = None if neleca < norb: acre = gen_map(gen_cre_linkstr, strsa, neleca, False) else: acre = None if nelecb < norb: bcre = gen_map(gen_cre_linkstr, strsb, nelecb, False) else: bcre = None def trans(ci1, aindex, bindex): if aindex is None or bindex is None: return None ma = len(aindex) mb = len(bindex) t1 = numpy.zeros((ma,mb)) for i in range(norb): signa = aindex[:,i,1] signb = bindex[:,i,1] maska = numpy.where(signa!=0)[0] maskb = numpy.where(signb!=0)[0] addra = aindex[maska,i,0] addrb = bindex[maskb,i,0] citmp = lib.take_2d(ci_coeff, addra, addrb) citmp *= signa[maska].reshape(-1,1) citmp *= signb[maskb] #: t1[addra.reshape(-1,1),addrb] += citmp lib.takebak_2d(t1, citmp, maska, maskb) for i in range(norb): signa = aindex[:,i,1] signb = bindex[:,i,1] maska = numpy.where(signa!=0)[0] maskb = numpy.where(signb!=0)[0] addra = aindex[maska,i,0] addrb = bindex[maskb,i,0] citmp = lib.take_2d(t1, maska, maskb) citmp *= signa[maska].reshape(-1,1) citmp *= signb[maskb] #: ci1[maska.reshape(-1,1), maskb] += citmp lib.takebak_2d(ci1, citmp, addra, addrb) ci1 = numpy.zeros_like(ci_coeff) trans(ci1, ades, bcre) # S+*S- trans(ci1, acre, bdes) # S-*S+ ci1 *= .5 ci1 += (neleca-nelecb)**2*.25*ci_coeff return _as_SCIvector(ci1, ci_strs)
[docs] def to_fci(civec_strs, norb, nelec): ci_coeff, nelec, ci_strs = _unpack(civec_strs, nelec) addrsa = [cistring.str2addr(norb, nelec[0], x) for x in ci_strs[0]] addrsb = [cistring.str2addr(norb, nelec[1], x) for x in ci_strs[1]] na = cistring.num_strings(norb, nelec[0]) nb = cistring.num_strings(norb, nelec[1]) ci0 = numpy.zeros((na,nb)) lib.takebak_2d(ci0, ci_coeff, addrsa, addrsb) return ci0
[docs] def from_fci(fcivec, ci_strs, norb, nelec): fcivec, nelec, ci_strs = _unpack(fcivec, nelec, ci_strs) addrsa = [cistring.str2addr(norb, nelec[0], x) for x in ci_strs[0]] addrsb = [cistring.str2addr(norb, nelec[1], x) for x in ci_strs[1]] na = cistring.num_strings(norb, nelec[0]) nb = cistring.num_strings(norb, nelec[1]) fcivec = fcivec.reshape(na,nb) civec = lib.take_2d(fcivec, addrsa, addrsb) return _as_SCIvector(civec, ci_strs)
[docs] class SelectedCI(direct_spin1.FCISolver): ci_coeff_cutoff = getattr(__config__, 'fci_selected_ci_SCI_ci_coeff_cutoff', .5e-3) select_cutoff = getattr(__config__, 'fci_selected_ci_SCI_select_cutoff', .5e-3) conv_tol = getattr(__config__, 'fci_selected_ci_SCI_conv_tol', 1e-9) start_tol = getattr(__config__, 'fci_selected_ci_SCI_start_tol', 3e-4) tol_decay_rate = getattr(__config__, 'fci_selected_ci_SCI_tol_decay_rate', 0.3) _keys = { 'ci_coeff_cutoff', 'select_cutoff', 'conv_tol', 'start_tol', 'tol_decay_rate', } def __init__(self, mol=None): direct_spin1.FCISolver.__init__(self, mol) ################################################## # don't modify the following attributes, they are not input options #self.converged = False #self.ci = None self._strs = None
[docs] def dump_flags(self, verbose=None): direct_spin1.FCISolver.dump_flags(self, verbose) logger.info(self, 'ci_coeff_cutoff %g', self.ci_coeff_cutoff) logger.info(self, 'select_cutoff %g', self.select_cutoff)
[docs] def contract_2e(self, eri, civec_strs, norb, nelec, link_index=None, **kwargs): # The argument civec_strs is a CI vector in function FCISolver.contract_2e. # Save and patch self._strs to make this contract_2e function compatible to # FCISolver.contract_2e. if getattr(civec_strs, '_strs', None) is not None: self._strs = civec_strs._strs else: assert (civec_strs.size == len(self._strs[0])*len(self._strs[1])) civec_strs = _as_SCIvector(civec_strs, self._strs) return contract_2e(eri, civec_strs, norb, nelec, link_index)
[docs] def get_init_guess(self, ci_strs, norb, nelec, nroots, hdiag): '''Initial guess is the single Slater determinant ''' na = len(ci_strs[0]) nb = len(ci_strs[1]) ci0 = direct_spin1._get_init_guess(na, nb, nroots, hdiag, nelec) return [_as_SCIvector(x, ci_strs) for x in ci0]
make_hdiag = staticmethod(make_hdiag) enlarge_space = enlarge_space kernel = kernel_float_space kernel_fixed_space = kernel_fixed_space # def approx_kernel(self, h1e, eri, norb, nelec, ci0=None, link_index=None, # tol=None, lindep=None, max_cycle=None, # max_memory=None, verbose=None, **kwargs): # ci_strs = getattr(ci0, '_strs', self._strs) # return self.kernel_fixed_space(h1e, eri, norb, nelec, ci_strs, # ci0, link_index, tol, lindep, 6, # max_memory, verbose, **kwargs)
[docs] @lib.with_doc(spin_square.__doc__) def spin_square(self, civec_strs, norb, nelec): nelec = direct_spin1._unpack_nelec(nelec, self.spin) return spin_square(_as_SCIvector_if_not(civec_strs, self._strs), norb, nelec)
[docs] def large_ci(self, civec_strs, norb, nelec, tol=.1, return_strs=True): nelec = direct_spin1._unpack_nelec(nelec, self.spin) ci, _, (strsa, strsb) = _unpack(civec_strs, nelec, self._strs) addra, addrb = numpy.where(abs(ci) > tol) if return_strs: strsa = [bin(x) for x in strsa[addra]] strsb = [bin(x) for x in strsb[addrb]] return list(zip(ci[addra,addrb], strsa, strsb)) else: occslsta = cistring._strs2occslst(strsa[addra], norb) occslstb = cistring._strs2occslst(strsb[addrb], norb) return list(zip(ci[addra,addrb], occslsta, occslstb))
[docs] def contract_ss(self, fcivec, norb, nelec): return contract_ss(fcivec, norb, nelec)
[docs] @lib.with_doc(make_rdm1s.__doc__) def make_rdm1s(self, civec_strs, norb, nelec, link_index=None): nelec = direct_spin1._unpack_nelec(nelec, self.spin) civec_strs = _as_SCIvector_if_not(civec_strs, self._strs) return make_rdm1s(civec_strs, norb, nelec, link_index)
[docs] @lib.with_doc(make_rdm1.__doc__) def make_rdm1(self, civec_strs, norb, nelec, link_index=None): nelec = direct_spin1._unpack_nelec(nelec, self.spin) rdm1a, rdm1b = self.make_rdm1s(civec_strs, norb, nelec, link_index) return rdm1a + rdm1b
[docs] @lib.with_doc(make_rdm2s.__doc__) def make_rdm2s(self, civec_strs, norb, nelec, link_index=None, **kwargs): nelec = direct_spin1._unpack_nelec(nelec, self.spin) civec_strs = _as_SCIvector_if_not(civec_strs, self._strs) return make_rdm2s(civec_strs, norb, nelec, link_index)
[docs] @lib.with_doc(make_rdm2.__doc__) def make_rdm2(self, civec_strs, norb, nelec, link_index=None, **kwargs): nelec = direct_spin1._unpack_nelec(nelec, self.spin) civec_strs = _as_SCIvector_if_not(civec_strs, self._strs) return make_rdm2(civec_strs, norb, nelec, link_index)
[docs] def make_rdm12s(self, civec_strs, norb, nelec, link_index=None, **kwargs): neleca, nelecb = nelec = direct_spin1._unpack_nelec(nelec, self.spin) civec_strs = _as_SCIvector_if_not(civec_strs, self._strs) dm2aa, dm2ab, dm2bb = make_rdm2s(civec_strs, norb, nelec, link_index) if neleca > 1 and nelecb > 1: dm1a = numpy.einsum('iikl->kl', dm2aa) / (neleca-1) dm1b = numpy.einsum('iikl->kl', dm2bb) / (nelecb-1) else: dm1a, dm1b = make_rdm1s(civec_strs, norb, nelec, link_index) return (dm1a, dm1b), (dm2aa, dm2ab, dm2bb)
[docs] def make_rdm12(self, civec_strs, norb, nelec, link_index=None, **kwargs): nelec = direct_spin1._unpack_nelec(nelec, self.spin) nelec_tot = sum(nelec) civec_strs = _as_SCIvector_if_not(civec_strs, self._strs) dm2 = make_rdm2(civec_strs, norb, nelec, link_index) if nelec_tot > 1: dm1 = numpy.einsum('iikl->kl', dm2) / (nelec_tot-1) else: dm1 = make_rdm1(civec_strs, norb, nelec, link_index) return dm1, dm2
[docs] @lib.with_doc(trans_rdm1s.__doc__) def trans_rdm1s(self, cibra, ciket, norb, nelec, link_index=None): nelec = direct_spin1._unpack_nelec(nelec, self.spin) cibra = _as_SCIvector_if_not(cibra, self._strs) ciket = _as_SCIvector_if_not(ciket, self._strs) return trans_rdm1s(cibra, ciket, norb, nelec, link_index)
[docs] @lib.with_doc(trans_rdm1.__doc__) def trans_rdm1(self, cibra, ciket, norb, nelec, link_index=None): nelec = direct_spin1._unpack_nelec(nelec, self.spin) cibra = _as_SCIvector_if_not(cibra, self._strs) ciket = _as_SCIvector_if_not(ciket, self._strs) return trans_rdm1(cibra, ciket, norb, nelec, link_index)
[docs] def gen_linkstr(self, norb, nelec, tril=True, spin=None, ci_strs=None): if spin is None: spin = self.spin if ci_strs is None: ci_strs = self._strs neleca, nelecb = direct_spin1._unpack_nelec(nelec, spin) if tril: cd_indexa = cre_des_linkstr_tril(ci_strs[0], norb, neleca) dd_indexa = des_des_linkstr_tril(ci_strs[0], norb, neleca) cd_indexb = cre_des_linkstr_tril(ci_strs[1], norb, nelecb) dd_indexb = des_des_linkstr_tril(ci_strs[1], norb, nelecb) else: cd_indexa = cre_des_linkstr(ci_strs[0], norb, neleca) dd_indexa = des_des_linkstr(ci_strs[0], norb, neleca) cd_indexb = cre_des_linkstr(ci_strs[1], norb, nelecb) dd_indexb = des_des_linkstr(ci_strs[1], norb, nelecb) return cd_indexa, dd_indexa, cd_indexb, dd_indexb
SCI = SelectedCI def _unpack(civec_strs, nelec, ci_strs=None, spin=None): neleca, nelecb = direct_spin1._unpack_nelec(nelec, spin) ci_strs = getattr(civec_strs, '_strs', ci_strs) if ci_strs is not None: strsa, strsb = ci_strs strsa = numpy.asarray(strsa) strsb = numpy.asarray(strsb) ci_strs = (strsa, strsb) return civec_strs, (neleca, nelecb), ci_strs def _all_linkstr_index(ci_strs, norb, nelec): cd_indexa = cre_des_linkstr_tril(ci_strs[0], norb, nelec[0]) dd_indexa = des_des_linkstr_tril(ci_strs[0], norb, nelec[0]) cd_indexb = cre_des_linkstr_tril(ci_strs[1], norb, nelec[1]) dd_indexb = des_des_linkstr_tril(ci_strs[1], norb, nelec[1]) return cd_indexa, dd_indexa, cd_indexb, dd_indexb # numpy.ndarray does not allow to attach attributes. Overwrite the # numpy.ndarray class to tag the ._strs attribute
[docs] class SCIvector(numpy.ndarray): '''An 2D np array for selected CI coefficients''' def __array_finalize__(self, obj): self._strs = getattr(obj, '_strs', None) # Special cases for ndarray when the array was modified (through ufunc) def __array_wrap__(self, out): if out.shape == self.shape: return out elif out.shape == (): # if ufunc returns a scalar return out[()] else: return out.view(numpy.ndarray)
def _as_SCIvector(civec, ci_strs): civec = civec.view(SCIvector) civec._strs = ci_strs return civec def _as_SCIvector_if_not(civec, ci_strs): if getattr(civec, '_strs', None) is None: civec = _as_SCIvector(civec, ci_strs) return civec if __name__ == '__main__': from functools import reduce from pyscf import gto from pyscf import scf from pyscf.fci import spin_op from pyscf.fci import addons mol = gto.Mole() mol.verbose = 0 mol.output = None mol.atom = [ ['H', ( 1.,-1. , 0. )], ['H', ( 0.,-1. ,-1. )], ['H', ( 1.,-0.5 ,-1. )], ['H', ( 0.,-0. ,-1. )], ['H', ( 1.,-0.5 , 0. )], ['H', ( 0., 1. , 1. )], ['H', ( 1., 2. , 3. )], ['H', ( 1., 2. , 4. )], ] mol.basis = 'sto-3g' mol.build() m = scf.RHF(mol) m.kernel() norb = m.mo_coeff.shape[1] nelec = mol.nelectron h1e = reduce(numpy.dot, (m.mo_coeff.T, m.get_hcore(), m.mo_coeff)) eri = ao2mo.kernel(m._eri, m.mo_coeff, compact=False) eri = eri.reshape(norb,norb,norb,norb) e1, c1 = kernel(h1e, eri, norb, nelec) e2, c2 = direct_spin1.kernel(h1e, eri, norb, nelec) print(e1, e1 - -11.894559902235565, 'diff to FCI', e1-e2) print(c1.shape, c2.shape) dm1_1 = make_rdm1(c1, norb, nelec) dm1_2 = direct_spin1.make_rdm1(c2, norb, nelec) print(abs(dm1_1 - dm1_2).sum()) dm2_1 = make_rdm2(c1, norb, nelec) dm2_2 = direct_spin1.make_rdm12(c2, norb, nelec)[1] print(abs(dm2_1 - dm2_2).sum()) myci = SelectedCI() e, c = kernel_fixed_space(myci, h1e, eri, norb, nelec, c1._strs) print(e - -11.894559902235565) print(myci.large_ci(c1, norb, nelec)) print(myci.spin_square(c1, norb, nelec)[0] - spin_op.spin_square0(to_fci(c1, norb, nelec), norb, nelec)[0]) myci = SelectedCI() myci = addons.fix_spin_(myci) e1, c1 = myci.kernel(h1e, eri, norb, nelec) print(e1, e1 - -11.89467612053687) print(myci.spin_square(c1, norb, nelec))