Commit b39a1d95 authored by hebrewsnabla's avatar hebrewsnabla
Browse files

update dftcasci

parent a9d6a9a2
...@@ -2,6 +2,10 @@ from automr import autocas, mcpdft ...@@ -2,6 +2,10 @@ from automr import autocas, mcpdft
from pyscf import mcscf, lib from pyscf import mcscf, lib
from pyscf.dft import rks, uks from pyscf.dft import rks, uks
import numpy as np import numpy as np
from functools import partial
print = partial(print, flush=True)
einsum = partial(np.einsum, optimize=True)
def get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1): 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 '''Coulomb + XC functional for UKS. See pyscf/dft/rks.py
...@@ -78,37 +82,39 @@ def get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1): ...@@ -78,37 +82,39 @@ def get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1):
vklr *= (alpha - hyb) vklr *= (alpha - hyb)
vk += vklr vk += vklr
vxc += vj - vk vxc += vj - vk
if ground_state: if ground_state:
exc -=(np.einsum('ij,ji', dm[0], vk[0]).real + exc0 = exc
ek = -(np.einsum('ij,ji', dm[0], vk[0]).real +
np.einsum('ij,ji', dm[1], vk[1]).real) * .5 np.einsum('ij,ji', dm[1], vk[1]).real) * .5
if ground_state: if ground_state:
ecoul = np.einsum('ij,ji', dm[0]+dm[1], vj).real * .5 ecoul = np.einsum('ij,ji', dm[0]+dm[1], vj).real * .5
else: else:
ecoul = None ecoul = None
vxc = lib.tag_array(vxc, ecoul=ecoul, exc=exc, vj=vj, vk=vk) #vxc = lib.tag_array(vxc, ecoul=ecoul, exc=exc, vj=vj, vk=vk)
return vxc return vxc, exc0, ecoul, ek
def ks_decomp(ks): def ks_decomp(ks):
e_nn = ks.energy_nuc() e_nn = ks.energy_nuc()
dm1 = ks.make_rdm1() dm1 = ks.make_rdm1()
dm1t = dm1[0] + dm1[1] dm1t = dm1[0] + dm1[1]
e_core = np.dot(ks.get_hcore(), dm1t) e_core = einsum('ij,ji->', ks.get_hcore(), dm1t)
vhf = uks.get_veff(ks, dm=dm1) vhf, e_xcdft, e_coul, e_xhf = get_veff(ks)
e_coul = vhf.ecoul #e_coul = vhf.ecoul
#print('e_mc : %15.8f' % e_mcscf) #print('e_mc : %15.8f' % e_mcscf)
print('e_nn : %15.8f' % e_nn) print('e_nn : %15.8f' % e_nn)
print('e_core : %15.8f' % e_core) print('e_core : %15.8f' % e_core)
print('e_coul : %15.8f' % e_coul) print('e_coul : %15.8f' % e_coul)
# print('e_x : %15.8f' % e_x) print('e_xhf : %15.8f' % e_xhf)
# print('e_otx : %15.8f' % e_otx) print('e_xcdft : %15.8f' % e_xcdft)
# print('e_otc : %15.8f' % e_otc) # print('e_otc : %15.8f' % e_otc)
# print('e_c : %15.8f' % e_c) # print('e_c : %15.8f' % e_c)
def dftcasci(ks, act_user): def dftcasci(ks, act_user):
#mo = ks.mo_coeff #mo = ks.mo_coeff
ks_decomp(ks)
nacto = act_user[0] nacto = act_user[0]
nacta, nactb = act_user[1] nacta, nactb = act_user[1]
nopen = nacta - nactb nopen = nacta - nactb
...@@ -124,5 +130,4 @@ def dftcasci(ks, act_user): ...@@ -124,5 +130,4 @@ def dftcasci(ks, act_user):
mc.verbose = 4 mc.verbose = 4
mc.kernel() mc.kernel()
ks_decomp(ks)
e_nn, e_core, e_coul, e_x, _, _, e_c = mcpdft.get_energy_decomposition(mc) e_nn, e_core, e_coul, e_x, _, _, e_c = mcpdft.get_energy_decomposition(mc)
\ No newline at end of file
...@@ -9,7 +9,7 @@ from automr import stability, dump_mat ...@@ -9,7 +9,7 @@ from automr import stability, dump_mat
import time import time
import copy import copy
def gen(xyz, bas, charge, spin, conv='tight', level_shift=0): def gen(xyz, bas, charge, spin, conv='tight', level_shift=0, xc=None):
'''for states other than singlets''' '''for states other than singlets'''
mol = gto.Mole() mol = gto.Mole()
mol.atom = xyz mol.atom = xyz
...@@ -19,7 +19,11 @@ def gen(xyz, bas, charge, spin, conv='tight', level_shift=0): ...@@ -19,7 +19,11 @@ def gen(xyz, bas, charge, spin, conv='tight', level_shift=0):
mol.verbose = 4 mol.verbose = 4
mol.build() mol.build()
mf = scf.UHF(mol) if xc is None:
mf = scf.UHF(mol)
else:
mf = dft.UKS(mol)
mf.xc = xc
if conv == 'loose': if conv == 'loose':
mf.conv_tol = 1e-6 mf.conv_tol = 1e-6
mf.max_cycle = 10 mf.max_cycle = 10
......
...@@ -43,6 +43,7 @@ def get_energy_decomposition (mc, ot=None, mo_coeff=None, ci=None): ...@@ -43,6 +43,7 @@ def get_energy_decomposition (mc, ot=None, mo_coeff=None, ci=None):
#else: #else:
if True: if True:
print('energy decomposition of MCPDFT')
e_nn, e_core, e_coul, e_x, e_otx, e_otc, e_c = _get_e_decomp (mc, ot, mo_coeff, ci, e_mcscf) 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_mc : %15.8f' % e_mcscf)
print('e_nn : %15.8f' % e_nn) print('e_nn : %15.8f' % e_nn)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment