add dftcasci

automr/autocas.py

@@ -38,7 +38,7 @@ def do_suhf(mf):
     return mf, no, noon, nacto, (nacta, nactb), ndb, nex
 
 
-def get_uno(mf, st='st2'):
+def get_uno(mf, st='st2', uks=False):
     mol = mf.mol 
     mo = mf.mo_coeff
     S = mol.intor_symmetric('int1e_ovlp')
@@ -58,7 +58,10 @@ def get_uno(mf, st='st2'):
     nacta = (nacto + nopen)//2
     nactb = (nacto - nopen)//2
     print('nacto, nacta, nactb: %d %d %d' % (nacto, nacta, nactb))
-    mf = mf.to_rhf()
+    if uks:
+        mf = mf.to_rks()
+    else:
+        mf = mf.to_rhf()
     mf.mo_coeff = unos
     print('UNO in active space')
     dump_mat.dump_mo(mol, unos[:,ndb:ndb+nacto], ncol=10)
@@ -250,6 +253,7 @@ def cas(mf, act_user=None, crazywfn=False, max_memory=2000, natorb=True, gvb=Fal
                 mf = sort_mo(mf, sort, ndb)
         dump_mat.dump_mo(mf.mol,mf.mo_coeff[:,ndb:ndb+nacto], ncol=10)
     nopen = nacta - nactb
+
     mc = mcscf.CASSCF(mf,nacto,(nacta,nactb))
     mc.fcisolver.max_memory = max_memory // 2
     mc.max_memory = max_memory // 2

automr/dftmr.py

+from automr import autocas, mcpdft
+from pyscf import mcscf, lib
+from pyscf.dft import rks, uks
+import numpy as np
+
+def get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1):
+    '''Coulomb + XC functional for UKS.  See pyscf/dft/rks.py
+    :func:`get_veff` fore more details.
+    '''
+    if mol is None: mol = ks.mol
+    if dm is None: dm = ks.make_rdm1()
+    if not isinstance(dm, np.ndarray):
+        dm = np.asarray(dm)
+    if dm.ndim == 2:  # RHF DM
+        dm = np.asarray((dm*.5,dm*.5))
+    ground_state = (dm.ndim == 3 and dm.shape[0] == 2)
+
+    #t0 = (logger.process_clock(), logger.perf_counter())
+
+    if ks.grids.coords is None:
+        ks.grids.build(with_non0tab=True)
+        if ks.small_rho_cutoff > 1e-20 and ground_state:
+            ks.grids = rks.prune_small_rho_grids_(ks, mol, dm[0]+dm[1], ks.grids)
+        #t0 = logger.timer(ks, 'setting up grids', *t0)
+    if ks.nlc != '':
+        if ks.nlcgrids.coords is None:
+            ks.nlcgrids.build(with_non0tab=True)
+            if ks.small_rho_cutoff > 1e-20 and ground_state:
+                ks.nlcgrids = rks.prune_small_rho_grids_(ks, mol, dm[0]+dm[1], ks.nlcgrids)
+            #t0 = logger.timer(ks, 'setting up nlc grids', *t0)
+
+    ni = ks._numint
+    if hermi == 2:  # because rho = 0
+        n, exc, vxc = (0,0), 0, 0
+    else:
+        max_memory = ks.max_memory - lib.current_memory()[0]
+        n, exc, vxc = ni.nr_uks(mol, ks.grids, ks.xc, dm, max_memory=max_memory)
+        if ks.nlc:
+            assert 'VV10' in ks.nlc.upper()
+            _, enlc, vnlc = ni.nr_rks(mol, ks.nlcgrids, ks.xc+'__'+ks.nlc, dm[0]+dm[1],
+                                      max_memory=max_memory)
+            exc += enlc
+            vxc += vnlc
+        #logger.debug(ks, 'nelec by numeric integration = %s', n)
+        #t0 = logger.timer(ks, 'vxc', *t0)
+
+    #enabling range-separated hybrids
+    omega, alpha, hyb = ni.rsh_and_hybrid_coeff(ks.xc, spin=mol.spin)
+
+    if abs(hyb) < 1e-10 and abs(alpha) < 1e-10:
+        vk = None
+        if (ks._eri is None and ks.direct_scf and
+            getattr(vhf_last, 'vj', None) is not None):
+            ddm = np.asarray(dm) - np.asarray(dm_last)
+            vj = ks.get_j(mol, ddm[0]+ddm[1], hermi)
+            vj += vhf_last.vj
+        else:
+            vj = ks.get_j(mol, dm[0]+dm[1], hermi)
+        vxc += vj
+    else:
+        if (ks._eri is None and ks.direct_scf and
+            getattr(vhf_last, 'vk', None) is not None):
+            ddm = np.asarray(dm) - np.asarray(dm_last)
+            vj, vk = ks.get_jk(mol, ddm, hermi)
+            vk *= hyb
+            if abs(omega) > 1e-10:
+                vklr = ks.get_k(mol, ddm, hermi, omega)
+                vklr *= (alpha - hyb)
+                vk += vklr
+            vj = vj[0] + vj[1] + vhf_last.vj
+            vk += vhf_last.vk
+        else:
+            vj, vk = ks.get_jk(mol, dm, hermi)
+            vj = vj[0] + vj[1]
+            vk *= hyb
+            if abs(omega) > 1e-10:
+                vklr = ks.get_k(mol, dm, hermi, omega)
+                vklr *= (alpha - hyb)
+                vk += vklr
+        vxc += vj - vk
+
+        if ground_state:
+            exc -=(np.einsum('ij,ji', dm[0], vk[0]).real +
+                   np.einsum('ij,ji', dm[1], vk[1]).real) * .5
+    if ground_state:
+        ecoul = np.einsum('ij,ji', dm[0]+dm[1], vj).real * .5
+    else:
+        ecoul = None
+
+    vxc = lib.tag_array(vxc, ecoul=ecoul, exc=exc, vj=vj, vk=vk)
+    return vxc
+
+def ks_decomp(ks):
+    e_nn = ks.energy_nuc()
+    dm1 = ks.make_rdm1()
+    dm1t = dm1[0] + dm1[1]
+    e_core = np.dot(ks.get_hcore(), dm1t)
+    vhf = uks.get_veff(ks, dm=dm1)
+    e_coul = vhf.ecoul
+    #print('e_mc   : %15.8f' % e_mcscf)
+    print('e_nn   : %15.8f' % e_nn)
+    print('e_core : %15.8f' % e_core)
+    print('e_coul : %15.8f' % e_coul)
+    #    print('e_x    : %15.8f' % e_x)
+    #    print('e_otx  : %15.8f' % e_otx)
+    #    print('e_otc  : %15.8f' % e_otc)
+    #    print('e_c    : %15.8f' % e_c)
+
+
+def dftcasci(ks, act_user):
+    #mo = ks.mo_coeff
+    nacto = act_user[0]
+    nacta, nactb = act_user[1]
+    nopen = nacta - nactb
+    mf, unos, unoon, _, _, ndb, nex = autocas.get_uno(ks, uks=True)
+
+    mc = mcscf.CASCI(ks,nacto,(nacta,nactb))
+    #mc.fcisolver.max_memory = max_memory // 2
+    #mc.max_memory = max_memory // 2
+    #mc.max_cycle = 200
+    mc.fcisolver.spin = nopen
+    mc.fcisolver.max_cycle = 100
+    mc.natorb = True
+    mc.verbose = 4
+    mc.kernel()
+
+    ks_decomp(ks)
+    e_nn, e_core, e_coul, e_x, _, _, e_c = mcpdft.get_energy_decomposition(mc)
\ No newline at end of file

automr/guess.py

@@ -84,7 +84,7 @@ def from_fchk(xyz, bas, fch, cycle=2):
     mf.kernel(dm)
     return mf
 
-def mix(xyz, bas, charge=0, conv='loose', cycle=5, skipstb=False):
+def mix(xyz, bas, charge=0, conv='loose', cycle=5, skipstb=False, xc=None):
     mol = gto.Mole()
     mol.atom = xyz
     #with open(xyz, 'r') as f:
@@ -102,12 +102,13 @@ def mix(xyz, bas, charge=0, conv='loose', cycle=5, skipstb=False):
     #mf.verbose = 4
     mf.kernel() # Guess by 1e is poor,
     #dm, mo_coeff, mo_energy, mo_occ = init_guess_by_1e(mf)
-    #mf.init_guess_breaksym = True
-    #mo = (mf.mo_coeff, mf.mo_coeff)
-    #occ = (mf.mo_occ, mf.mo_occ)
     print('**** generating mix guess ****')
     dm_mix = init_guess_mixed(mf.mo_coeff, mf.mo_occ)
-    mf_mix = scf.UHF(mol)
+    if xc is None:
+        mf_mix = scf.UHF(mol)
+    else:
+        mf_mix = dft.UKS(mol)
+        mf_mix.xc = xc
     #mf_mix.verbose = 4
     if conv == 'loose':
         mf_mix.conv_tol = 1e-3

automr/mcpdft.py

@@ -12,7 +12,7 @@ print = partial(print, flush=True)
 einsum = partial(np.einsum, optimize=True)
 
 # mc is mcscf.CASSCF obj
-def get_energy_decomposition (mc, ot, mo_coeff=None, ci=None):
+def get_energy_decomposition (mc, ot=None, mo_coeff=None, ci=None):
     ''' Compute a decomposition of the MC-PDFT energy into nuclear potential, core, Coulomb, exchange,
     and correlation terms. The exchange-correlation energy at the MC-SCF level is also returned.
     This is not meant to work with MC-PDFT methods that are already hybrids. Most return arguments
@@ -25,7 +25,8 @@ def get_energy_decomposition (mc, ot, mo_coeff=None, ci=None):
     #    e_tot = mc.e_states
     #e_nuc = mc._scf.energy_nuc ()
     #h = mc.get_hcore ()
-    xfnal, cfnal = ot.split_x_c ()
+    #if ot is not None:
+    #    xfnal, cfnal = ot.split_x_c ()
     #if isinstance (mc, StateAverageMCSCFSolver):
     #    e_core = []
     #    e_coul = []
@@ -42,7 +43,7 @@ def get_energy_decomposition (mc, ot, mo_coeff=None, ci=None):
     
     #else:
     if True:
-        e_nn, e_core, e_coul, e_x, e_otx, e_otc, e_c = _get_e_decomp (mc, ot, mo_coeff, ci, e_mcscf, xfnal, cfnal)
+        e_nn, e_core, e_coul, e_x, e_otx, e_otc, e_c = _get_e_decomp (mc, ot, mo_coeff, ci, e_mcscf)
         print('e_mc   : %15.8f' % e_mcscf)
         print('e_nn   : %15.8f' % e_nn)
         print('e_core : %15.8f' % e_core)
@@ -53,7 +54,9 @@ def get_energy_decomposition (mc, ot, mo_coeff=None, ci=None):
         print('e_c    : %15.8f' % e_c)
     return e_nn, e_core, e_coul, e_x, e_otx, e_otc, e_c
 
-def _get_e_decomp (mc, ot, mo_coeff, ci, e_mcscf, xfnal, cfnal):
+def _get_e_decomp (mc, ot, mo_coeff, ci, e_mcscf):
+    if ot is not None:
+        xfnal, cfnal = ot.split_x_c ()
     ncore, ncas, nelecas = mc.ncore, mc.ncas, mc.nelecas
     _rdms = mcscf.CASCI (mc._scf, ncas, nelecas)
     _rdms.fcisolver = fci.solver (mc._scf.mol, singlet = False, symm = False)
@@ -70,11 +73,13 @@ def _get_e_decomp (mc, ot, mo_coeff, ci, e_mcscf, xfnal, cfnal):
     #print(j.shape, k.shape, dm1.shape)
     e_coul = np.tensordot (j[0] + j[1], dm1, axes=2) / 2
     e_x = -(np.tensordot (k[0], dm1s[0]) + np.tensordot (k[1], dm1s[1])) / 2
-    adm1s = np.stack (_casdms.make_rdm1s (ci, ncas, nelecas), axis=0)
-    adm2 = get_2CDM_from_2RDM (_casdms.make_rdm12 (_rdms.ci, ncas, nelecas)[1], adm1s)
-    mo_cas = mo_coeff[:,ncore:][:,:ncas]
-    e_otx = get_E_ot (xfnal, dm1s, adm2, mo_cas, max_memory=mc.max_memory)
-    e_otc = get_E_ot (cfnal, dm1s, adm2, mo_cas, max_memory=mc.max_memory)
+    e_otx = 0.0; e_otc = 0.0
+    if ot is not None:
+        adm1s = np.stack (_casdms.make_rdm1s (ci, ncas, nelecas), axis=0)
+        adm2 = get_2CDM_from_2RDM (_casdms.make_rdm12 (_rdms.ci, ncas, nelecas)[1], adm1s)
+        mo_cas = mo_coeff[:,ncore:][:,:ncas]
+        e_otx = get_E_ot (xfnal, dm1s, adm2, mo_cas, max_memory=mc.max_memory)
+        e_otc = get_E_ot (cfnal, dm1s, adm2, mo_cas, max_memory=mc.max_memory)
     e_c = e_mcscf - e_nn - e_core - e_coul - e_x
     return e_nn, e_core, e_coul, e_x, e_otx, e_otc, e_c