#!/usr/bin/env python
#
# This code was copied from the data generation program of Tencent Alchemy
# project (https://github.com/tencent-alchemy).
#
#
# Copyright 2019 Tencent America LLC. 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>
#
'''
Non-relativistic UHF analytical Hessian
'''
import numpy
import numpy as np
import scipy.linalg
from pyscf import lib
from pyscf.lib import logger
from pyscf import ao2mo
from pyscf.hessian import rhf as rhf_hess
from pyscf.hessian import uhf as uhf_hess
from pyscf.df.hessian.rhf import _load_dim0, _pinv
from pyscf.df.grad.rhf import (_int3c_wrapper, _gen_metric_solver,
LINEAR_DEP_THRESHOLD)
[docs]
def partial_hess_elec(hessobj, mo_energy=None, mo_coeff=None, mo_occ=None,
atmlst=None, max_memory=4000, verbose=None):
e1, ej, ek = _partial_hess_ejk(hessobj, mo_energy, mo_coeff, mo_occ,
atmlst, max_memory, verbose, True)
return e1 + ej - ek
def _partial_hess_ejk(hessobj, mo_energy=None, mo_coeff=None, mo_occ=None,
atmlst=None, max_memory=4000, verbose=None, with_k=True):
log = logger.new_logger(hessobj, verbose)
time0 = t1 = (logger.process_clock(), logger.perf_counter())
mol = hessobj.mol
mf = hessobj.base
if mo_energy is None: mo_energy = mf.mo_energy
if mo_occ is None: mo_occ = mf.mo_occ
if mo_coeff is None: mo_coeff = mf.mo_coeff
if atmlst is None: atmlst = range(mol.natm)
nao, nmo = mo_coeff[0].shape
mocca = mo_coeff[0][:,mo_occ[0]>0]
moccb = mo_coeff[1][:,mo_occ[1]>0]
nocca = mocca.shape[1]
noccb = moccb.shape[1]
dm0a = numpy.dot(mocca, mocca.T)
dm0b = numpy.dot(moccb, moccb.T)
dm0 = dm0a + dm0b
# Energy weighted density matrix
mo_ea = mo_energy[0][mo_occ[0]>0]
mo_eb = mo_energy[1][mo_occ[1]>0]
dme0 = numpy.einsum('pi,qi,i->pq', mocca, mocca, mo_ea)
dme0+= numpy.einsum('pi,qi,i->pq', moccb, moccb, mo_eb)
auxmol = hessobj.base.with_df.auxmol
naux = auxmol.nao
nbas = mol.nbas
auxslices = auxmol.aoslice_by_atom()
aoslices = mol.aoslice_by_atom()
aux_loc = auxmol.ao_loc
blksize = min(480, hessobj.max_memory*.3e6/8/nao**2)
aux_ranges = ao2mo.outcore.balance_partition(auxmol.ao_loc, blksize)
hcore_deriv = hessobj.hcore_generator(mol)
s1aa, s1ab, s1a = rhf_hess.get_ovlp(mol)
ftmp = lib.H5TmpFile()
get_int3c = _int3c_wrapper(mol, auxmol, 'int3c2e', 's1')
# Without RI basis response
# (20|0)(0|00)
# (11|0)(0|00)
# (10|0)(0|10)
int2c = auxmol.intor('int2c2e', aosym='s1')
solve_j2c = _gen_metric_solver(int2c)
int2c_ip1 = auxmol.intor('int2c2e_ip1', aosym='s1')
rhoj0_P = 0
if with_k:
if hessobj.max_memory*.4e6/8 < naux*nocca*(nocca+nao):
raise RuntimeError('Memory not enough. You need to increase mol.max_memory')
rhok0a_Pl_ = np.empty((naux,nao,nocca))
rhok0b_Pl_ = np.empty((naux,nao,noccb))
for i, (shl0, shl1, p0, p1) in enumerate(aoslices):
int3c = get_int3c((shl0, shl1, 0, nbas, 0, auxmol.nbas))
rhoj0_P += np.einsum('klp,kl->p', int3c, dm0[p0:p1])
if with_k:
tmp = lib.einsum('ijp,jk->pik', int3c, mocca)
tmp = solve_j2c(tmp.reshape(naux,-1))
rhok0a_Pl_[:,p0:p1] = tmp.reshape(naux,p1-p0,nocca)
tmp = lib.einsum('ijp,jk->pik', int3c, moccb)
tmp = solve_j2c(tmp.reshape(naux,-1))
rhok0b_Pl_[:,p0:p1] = tmp.reshape(naux,p1-p0,noccb)
int3c = tmp = None
rhoj0_P = solve_j2c(rhoj0_P)
get_int3c_ipip1 = _int3c_wrapper(mol, auxmol, 'int3c2e_ipip1', 's1')
vj1_diag = 0
vk1a_diag = 0
vk1b_diag = 0
for shl0, shl1, nL in aux_ranges:
shls_slice = (0, nbas, 0, nbas, shl0, shl1)
p0, p1 = aux_loc[shl0], aux_loc[shl1]
int3c_ipip1 = get_int3c_ipip1(shls_slice)
vj1_diag += np.einsum('xijp,p->xij', int3c_ipip1, rhoj0_P[p0:p1]).reshape(3,3,nao,nao)
if with_k:
tmp = lib.einsum('Plj,Jj->PlJ', rhok0a_Pl_[p0:p1], mocca)
vk1a_diag += lib.einsum('xijp,plj->xil', int3c_ipip1, tmp).reshape(3,3,nao,nao)
tmp = lib.einsum('Plj,Jj->PlJ', rhok0b_Pl_[p0:p1], moccb)
vk1b_diag += lib.einsum('xijp,plj->xil', int3c_ipip1, tmp).reshape(3,3,nao,nao)
t1 = log.timer_debug1('contracting int2e_ipip1', *t1)
int3c_ipip1 = get_int3c_ipip1 = tmp = None
get_int3c_ip1 = _int3c_wrapper(mol, auxmol, 'int3c2e_ip1', 's1')
rho_ip1 = ftmp.create_dataset('rho_ip1', (nao,nao,naux,3), 'f8')
rhoka_ip1_IkP = ftmp.create_group('rhoka_ip1_IkP')
rhokb_ip1_IkP = ftmp.create_group('rhokb_ip1_IkP')
rhoka_ip1_PkI = ftmp.create_group('rhoka_ip1_PkI')
rhokb_ip1_PkI = ftmp.create_group('rhokb_ip1_PkI')
rhoj1 = np.empty((mol.natm,naux,3))
wj1 = np.empty((mol.natm,naux,3))
for i0, ia in enumerate(atmlst):
shl0, shl1, p0, p1 = aoslices[ia]
shls_slice = (shl0, shl1, 0, nbas, 0, auxmol.nbas)
int3c_ip1 = get_int3c_ip1(shls_slice)
tmp_ip1 = solve_j2c(int3c_ip1.reshape(-1,naux).T).reshape(naux,3,p1-p0,nao)
rhoj1[i0] = np.einsum('pxij,ji->px', tmp_ip1, dm0[:,p0:p1])
wj1[i0] = np.einsum('xijp,ji->px', int3c_ip1, dm0[:,p0:p1])
rho_ip1[p0:p1] = tmp_ip1.transpose(2,3,0,1)
if with_k:
tmp = lib.einsum('pykl,li->ikpy', tmp_ip1, dm0a)
rhoka_ip1_IkP['%.4d'%ia] = tmp
rhoka_ip1_PkI['%.4d'%ia] = tmp.transpose(2,1,0,3)
tmp = None
tmp = lib.einsum('pykl,li->ikpy', tmp_ip1, dm0b)
rhokb_ip1_IkP['%.4d'%ia] = tmp
rhokb_ip1_PkI['%.4d'%ia] = tmp.transpose(2,1,0,3)
tmp = None
ej = lib.einsum('ipx,jpy->ijxy', rhoj1, wj1) * 4
ek = np.zeros_like(ej)
e1 = np.zeros_like(ej)
rhoj1 = wj1 = None
if with_k:
vk2a_buf = 0
vk2b_buf = 0
for shl0, shl1, nL in aux_ranges:
shls_slice = (0, nbas, 0, nbas, shl0, shl1)
p0, p1 = aux_loc[shl0], aux_loc[shl1]
int3c_ip1 = get_int3c_ip1(shls_slice)
vk2a_buf += lib.einsum('xijp,pkjy->xyki', int3c_ip1, _load_dim0(rhoka_ip1_PkI, p0, p1))
vk2b_buf += lib.einsum('xijp,pkjy->xyki', int3c_ip1, _load_dim0(rhokb_ip1_PkI, p0, p1))
get_int3c_ip2 = _int3c_wrapper(mol, auxmol, 'int3c2e_ip2', 's1')
wj_ip2 = np.empty((naux,3))
wka_ip2_Ipk = ftmp.create_dataset('wka_ip2', (nao,naux,3,nao), 'f8')
wkb_ip2_Ipk = ftmp.create_dataset('wkb_ip2', (nao,naux,3,nao), 'f8')
if hessobj.auxbasis_response > 1:
wka_ip2_P__ = np.empty((naux,3,nocca,nocca))
wkb_ip2_P__ = np.empty((naux,3,noccb,noccb))
for shl0, shl1, nL in aux_ranges:
shls_slice = (0, nbas, 0, nbas, shl0, shl1)
p0, p1 = aux_loc[shl0], aux_loc[shl1]
int3c_ip2 = get_int3c_ip2(shls_slice)
wj_ip2[p0:p1] = np.einsum('yklp,lk->py', int3c_ip2, dm0)
if with_k:
wka_ip2_Ipk[:,p0:p1] = lib.einsum('yklp,il->ipyk', int3c_ip2, dm0a)
wkb_ip2_Ipk[:,p0:p1] = lib.einsum('yklp,il->ipyk', int3c_ip2, dm0b)
if hessobj.auxbasis_response > 1:
wka_ip2_P__[p0:p1] = lib.einsum('xuvp,ui,vj->pxij', int3c_ip2, mocca, mocca)
wkb_ip2_P__[p0:p1] = lib.einsum('xuvp,ui,vj->pxij', int3c_ip2, moccb, moccb)
int3c_ip2 = None
if hessobj.auxbasis_response > 1:
get_int3c_ipip2 = _int3c_wrapper(mol, auxmol, 'int3c2e_ipip2', 's1')
rhok0a_P__ = lib.einsum('plj,li->pij', rhok0a_Pl_, mocca)
rhok0b_P__ = lib.einsum('plj,li->pij', rhok0b_Pl_, moccb)
rho2c_0 = lib.einsum('pij,qij->pq', rhok0a_P__, rhok0a_P__)
rho2c_0 += lib.einsum('pij,qij->pq', rhok0b_P__, rhok0b_P__)
int2c_inv = _pinv(int2c, lindep=LINEAR_DEP_THRESHOLD)
int2c_ipip1 = auxmol.intor('int2c2e_ipip1', aosym='s1')
int2c_ip_ip = lib.einsum('xpq,qr,ysr->xyps', int2c_ip1, int2c_inv, int2c_ip1)
int2c_ip_ip -= auxmol.intor('int2c2e_ip1ip2', aosym='s1').reshape(3,3,naux,naux)
int2c = solve_j2c = None
get_int3c_ipvip1 = _int3c_wrapper(mol, auxmol, 'int3c2e_ipvip1', 's1')
get_int3c_ip1ip2 = _int3c_wrapper(mol, auxmol, 'int3c2e_ip1ip2', 's1')
for i0, ia in enumerate(atmlst):
shl0, shl1, p0, p1 = aoslices[ia]
shls_slice = (shl0, shl1, 0, nbas, 0, auxmol.nbas)
# (10|0)(0|10) without response of RI basis
if with_k:
int3c_ip1 = get_int3c_ip1(shls_slice)
vk1a = lib.einsum('xijp,ikpy->xykj', int3c_ip1, _load_dim0(rhoka_ip1_IkP, p0, p1))
vk1b = lib.einsum('xijp,ikpy->xykj', int3c_ip1, _load_dim0(rhokb_ip1_IkP, p0, p1))
vk1a[:,:,:,p0:p1] += vk2a_buf[:,:,:,p0:p1]
vk1b[:,:,:,p0:p1] += vk2b_buf[:,:,:,p0:p1]
t1 = log.timer_debug1('contracting int2e_ip1ip2 for atom %d'%ia, *t1)
int3c_ip1 = None
# (11|0)(0|00) without response of RI basis
int3c_ipvip1 = get_int3c_ipvip1(shls_slice)
vj1 = np.einsum('xijp,p->xji', int3c_ipvip1, rhoj0_P).reshape(3,3,nao,p1-p0)
if with_k:
tmp = lib.einsum('pki,ji->pkj', rhok0a_Pl_, mocca[p0:p1])
vk1a += lib.einsum('xijp,pki->xjk', int3c_ipvip1, tmp).reshape(3,3,nao,nao)
tmp = lib.einsum('pki,ji->pkj', rhok0b_Pl_, moccb[p0:p1])
vk1b += lib.einsum('xijp,pki->xjk', int3c_ipvip1, tmp).reshape(3,3,nao,nao)
t1 = log.timer_debug1('contracting int2e_ipvip1 for atom %d'%ia, *t1)
int3c_ipvip1 = tmp = None
e1[i0,i0] -= numpy.einsum('xypq,pq->xy', s1aa[:,:,p0:p1], dme0[p0:p1])*2
ej[i0,i0] += numpy.einsum('xypq,pq->xy', vj1_diag[:,:,p0:p1], dm0[p0:p1])*2
if with_k:
ek[i0,i0] += numpy.einsum('xypq,pq->xy', vk1a_diag[:,:,p0:p1], dm0a[p0:p1])*2
ek[i0,i0] += numpy.einsum('xypq,pq->xy', vk1b_diag[:,:,p0:p1], dm0b[p0:p1])*2
for j0, ja in enumerate(atmlst[:i0+1]):
q0, q1 = aoslices[ja][2:]
ej[i0,j0] += numpy.einsum('xypq,pq->xy', vj1[:,:,q0:q1], dm0[q0:q1,p0:p1])*2
e1[i0,j0] -= numpy.einsum('xypq,pq->xy', s1ab[:,:,p0:p1,q0:q1], dme0[p0:p1,q0:q1])*2
if with_k:
ek[i0,j0] += numpy.einsum('xypq,pq->xy', vk1a[:,:,q0:q1], dm0a[q0:q1])*2
ek[i0,j0] += numpy.einsum('xypq,pq->xy', vk1b[:,:,q0:q1], dm0b[q0:q1])*2
h1ao = hcore_deriv(ia, ja)
e1[i0,j0] += numpy.einsum('xypq,pq->xy', h1ao, dm0)
#
# The first order RI basis response
#
# (10|1)(0|00)
# (10|0)(1|0)(0|00)
# (10|0)(0|1)(0|00)
# (10|0)(1|00)
#
if hessobj.auxbasis_response:
wk1_Pij = rho_ip1[p0:p1].transpose(2,3,0,1)
rhoj1_P = np.einsum('pxij,ji->px', wk1_Pij, dm0[:,p0:p1])
# (10|1)(0|0)(0|00)
int3c_ip1ip2 = get_int3c_ip1ip2(shls_slice)
wj11_p = np.einsum('xijp,ji->xp', int3c_ip1ip2, dm0[:,p0:p1])
# (10|0)(1|0)(0|00)
wj0_01 = np.einsum('ypq,q->yp', int2c_ip1, rhoj0_P)
if with_k:
rhok0_P_I = lib.einsum('plj,il->pji', rhok0a_Pl_, dm0a[p0:p1])
rhok0_PJI = lib.einsum('pji,Jj->pJi', rhok0_P_I, mocca)
rhok0_P_I = lib.einsum('plj,il->pji', rhok0b_Pl_, dm0b[p0:p1])
rhok0_PJI+= lib.einsum('pji,Jj->pJi', rhok0_P_I, moccb)
wk1_pJI = lib.einsum('ypq,qji->ypji', int2c_ip1, rhok0_PJI)
wk1_IpJ = lib.einsum('ipyk,kj->ipyj', wka_ip2_Ipk[p0:p1], dm0a)
wk1_IpJ+= lib.einsum('ipyk,kj->ipyj', wkb_ip2_Ipk[p0:p1], dm0b)
rho2c_PQ = lib.einsum('pxij,qji->xqp', wk1_Pij, rhok0_PJI)
for j0, (q0, q1) in enumerate(auxslices[:,2:]):
# (10|1)(0|00)
_ej = np.einsum('xp,p->x', wj11_p[:,q0:q1], rhoj0_P[q0:q1]).reshape(3,3)
# (10|0)(0|1)(0|00)
_ej -= lib.einsum('yqp,q,px->xy', int2c_ip1[:,q0:q1], rhoj0_P[q0:q1], rhoj1_P)
# (10|0)(1|0)(0|00)
_ej -= lib.einsum('px,yp->xy', rhoj1_P[q0:q1], wj0_01[:,q0:q1])
# (10|0)(1|00)
_ej += lib.einsum('px,py->xy', rhoj1_P[q0:q1], wj_ip2[q0:q1])
if hessobj.auxbasis_response > 1:
ej[i0,j0] += _ej * 2
ej[j0,i0] += _ej.T * 2
else:
ej[i0,j0] += _ej
ej[j0,i0] += _ej.T
if with_k:
_ek = lib.einsum('xijp,pji->x', int3c_ip1ip2[:,:,:,q0:q1],
rhok0_PJI[q0:q1]).reshape(3,3)
_ek -= lib.einsum('pxij,ypji->xy', wk1_Pij[q0:q1], wk1_pJI[:,q0:q1])
_ek -= lib.einsum('xqp,yqp->xy', rho2c_PQ[:,q0:q1], int2c_ip1[:,q0:q1])
_ek += lib.einsum('pxij,ipyj->xy', wk1_Pij[q0:q1], wk1_IpJ[:,q0:q1])
if hessobj.auxbasis_response > 1:
ek[i0,j0] += _ek * 2
ek[j0,i0] += _ek.T * 2
else:
ek[i0,j0] += _ek
ek[j0,i0] += _ek.T
int3c_ip1ip2 = None
#
# The second order RI basis response
#
if hessobj.auxbasis_response > 1:
# (00|2)(0|00)
# (00|0)(2|0)(0|00)
shl0, shl1, p0, p1 = auxslices[ia]
shls_slice = (0, nbas, 0, nbas, shl0, shl1)
int3c_ipip2 = get_int3c_ipip2(shls_slice)
ej[i0,i0] += np.einsum('xijp,ji,p->x', int3c_ipip2, dm0, rhoj0_P[p0:p1]).reshape(3,3)
ej[i0,i0] -= np.einsum('p,xpq,q->x', rhoj0_P[p0:p1], int2c_ipip1[:,p0:p1], rhoj0_P).reshape(3,3)
if with_k:
rhok0_PJI = lib.einsum('Pij,Jj,Ii->PJI', rhok0a_P__[p0:p1], mocca, mocca)
rhok0_PJI += lib.einsum('Pij,Jj,Ii->PJI', rhok0b_P__[p0:p1], moccb, moccb)
ek[i0,i0] += np.einsum('xijp,pij->x', int3c_ipip2, rhok0_PJI).reshape(3,3)
ek[i0,i0] -= np.einsum('pq,xpq->x', rho2c_0[p0:p1], int2c_ipip1[:,p0:p1]).reshape(3,3)
rhok0_PJI = None
# (00|0)(1|1)(0|00)
# (00|1)(1|0)(0|00)
# (00|1)(0|1)(0|00)
# (00|1)(1|00)
rhoj1 = lib.einsum('px,pq->xq', wj_ip2[p0:p1], int2c_inv[p0:p1])
# (00|0)(0|1)(1|0)(0|00)
rhoj0_01 = lib.einsum('xp,pq->xq', wj0_01[:,p0:p1], int2c_inv[p0:p1])
# (00|0)(1|0)(1|0)(0|00)
ip1_2c_2c = lib.einsum('xpq,qr->xpr', int2c_ip1[:,p0:p1], int2c_inv)
rhoj0_10 = lib.einsum('p,xpq->xq', rhoj0_P[p0:p1], ip1_2c_2c)
if with_k:
# (00|0)(0|1)(1|0)(0|00)
ip1_rho2c = .5 * lib.einsum('xpq,qr->xpr', int2c_ip1[:,p0:p1], rho2c_0)
rho2c_1 = lib.einsum('xrq,rp->xpq', ip1_rho2c, int2c_inv[p0:p1])
# (00|0)(1|0)(1|0)(0|00)
rho2c_1 += lib.einsum('xrp,rq->xpq', ip1_2c_2c, rho2c_0[p0:p1])
# (00|1)(0|1)(0|00)
# (00|1)(1|0)(0|00)
int3c_ip2 = get_int3c_ip2(shls_slice)
tmp = lib.einsum('xuvr,vj,ui,qij,rp->xpq', int3c_ip2,
mocca, mocca, rhok0a_P__, int2c_inv[p0:p1])
tmp+= lib.einsum('xuvr,vj,ui,qij,rp->xpq', int3c_ip2,
moccb, moccb, rhok0b_P__, int2c_inv[p0:p1])
rho2c_1 -= tmp
rho2c_1 -= tmp.transpose(0,2,1)
int3c_ip2 = tmp = None
for j0, (q0, q1) in enumerate(auxslices[:,2:]):
_ej = 0
# (00|0)(1|1)(0|00)
# (00|0)(1|0)(0|1)(0|00)
_ej += .5 * np.einsum('p,xypq,q->xy', rhoj0_P[p0:p1], int2c_ip_ip[:,:,p0:p1,q0:q1], rhoj0_P[q0:q1])
# (00|1)(1|0)(0|00)
_ej -= lib.einsum('xp,yp->xy', rhoj1[:,q0:q1], wj0_01[:,q0:q1])
# (00|1)(1|00)
_ej += .5 * lib.einsum('xp,py->xy', rhoj1[:,q0:q1], wj_ip2[q0:q1])
# (00|0)(0|1)(1|0)(0|00)
_ej += .5 * np.einsum('xp,yp->xy', rhoj0_01[:,q0:q1], wj0_01[:,q0:q1])
# (00|1)(0|1)(0|00)
_ej -= lib.einsum('yqp,q,xp->xy', int2c_ip1[:,q0:q1], rhoj0_P[q0:q1], rhoj1)
# (00|0)(1|0)(1|0)(0|00)
_ej += np.einsum('xp,yp->xy', rhoj0_10[:,q0:q1], wj0_01[:,q0:q1])
ej[i0,j0] += _ej
ej[j0,i0] += _ej.T
if with_k:
# (00|0)(1|1)(0|00)
# (00|0)(1|0)(0|1)(0|00)
_ek = .5 * np.einsum('pq,xypq->xy', rho2c_0[p0:p1,q0:q1], int2c_ip_ip[:,:,p0:p1,q0:q1])
# (00|1)(0|1)(0|00)
# (00|1)(1|0)(0|00)
# (00|0)(0|1)(1|0)(0|00)
# (00|0)(1|0)(1|0)(0|00)
_ek += np.einsum('xpq,ypq->xy', rho2c_1[:,q0:q1], int2c_ip1[:,q0:q1])
# (00|1)(1|00)
_ek += .5 * lib.einsum('pxij,pq,qyij->xy', wka_ip2_P__[p0:p1],
int2c_inv[p0:p1,q0:q1], wka_ip2_P__[q0:q1])
_ek += .5 * lib.einsum('pxij,pq,qyij->xy', wkb_ip2_P__[p0:p1],
int2c_inv[p0:p1,q0:q1], wkb_ip2_P__[q0:q1])
ek[i0,j0] += _ek
ek[j0,i0] += _ek.T
for i0, ia in enumerate(atmlst):
for j0 in range(i0):
e1[j0,i0] = e1[i0,j0].T
ej[j0,i0] = ej[i0,j0].T
ek[j0,i0] = ek[i0,j0].T
log.timer('UHF partial hessian', *time0)
return e1, ej, ek
[docs]
def make_h1(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None):
mol = hessobj.mol
if atmlst is None:
atmlst = range(mol.natm)
h1aoa = [None] * mol.natm
h1aob = [None] * mol.natm
for ia, h1, vj1, vk1 in _gen_jk(hessobj, mo_coeff, mo_occ, chkfile,
atmlst, verbose, True):
vk1a, vk1b = vk1
h1aoa[ia] = h1 + vj1 - vk1a
h1aob[ia] = h1 + vj1 - vk1b
return (h1aoa,h1aob)
def _gen_jk(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None,
verbose=None, with_k=True):
mol = hessobj.mol
if atmlst is None:
atmlst = range(mol.natm)
auxmol = hessobj.base.with_df.auxmol
nbas = mol.nbas
auxslices = auxmol.aoslice_by_atom()
aux_loc = auxmol.ao_loc
nao, nmo = mo_coeff[0].shape
mocca = mo_coeff[0][:,mo_occ[0]>0]
moccb = mo_coeff[1][:,mo_occ[1]>0]
nocca = mocca.shape[1]
noccb = moccb.shape[1]
dm0a = numpy.dot(mocca, mocca.T)
dm0b = numpy.dot(moccb, moccb.T)
dm0 = dm0a + dm0b
hcore_deriv = hessobj.base.nuc_grad_method().hcore_generator(mol)
get_int3c = _int3c_wrapper(mol, auxmol, 'int3c2e', 's1')
aoslices = mol.aoslice_by_atom()
naux = auxmol.nao
ftmp = lib.H5TmpFile()
rho0_Pij = ftmp.create_group('rho0_Pij')
wj_Pij = ftmp.create_group('wj_Pij')
int2c = auxmol.intor('int2c2e', aosym='s1')
solve_j2c = _gen_metric_solver(int2c)
int2c = None
int2c_ip1 = auxmol.intor('int2c2e_ip1', aosym='s1')
rhoj0_P = 0
if with_k:
rhok0a_Pl_ = np.empty((naux,nao,nocca))
rhok0b_Pl_ = np.empty((naux,nao,noccb))
for i, (shl0, shl1, p0, p1) in enumerate(aoslices):
int3c = get_int3c((shl0, shl1, 0, nbas, 0, auxmol.nbas))
coef3c = solve_j2c(int3c.reshape(-1,naux).T)
rho0_Pij['%.4d'%i] = coef3c = coef3c.reshape(naux,p1-p0,nao)
rhoj0_P += np.einsum('pkl,kl->p', coef3c, dm0[p0:p1])
if with_k:
rhok0a_Pl_[:,p0:p1] = lib.einsum('pij,jk->pik', coef3c, mocca)
rhok0b_Pl_[:,p0:p1] = lib.einsum('pij,jk->pik', coef3c, moccb)
if hessobj.auxbasis_response:
wj_Pij['%.4d'%i] = lib.einsum('xqp,pij->qixj', int2c_ip1, coef3c)
int3c = coef3c = None
get_int3c_ip1 = _int3c_wrapper(mol, auxmol, 'int3c2e_ip1', 's1')
get_int3c_ip2 = _int3c_wrapper(mol, auxmol, 'int3c2e_ip2', 's1')
aux_ranges = ao2mo.outcore.balance_partition(auxmol.ao_loc, 480)
vk1a_buf = np.zeros((3,nao,nao))
vk1b_buf = np.zeros((3,nao,nao))
vj1_buf = np.zeros((mol.natm,3,nao,nao))
for shl0, shl1, nL in aux_ranges:
shls_slice = (0, nbas, 0, nbas, shl0, shl1)
p0, p1 = aux_loc[shl0], aux_loc[shl1]
int3c_ip1 = get_int3c_ip1(shls_slice)
coef3c = _load_dim0(rho0_Pij, p0, p1)
for i, (shl0, shl1, q0, q1) in enumerate(aoslices):
wj1 = np.einsum('xijp,ji->xp', int3c_ip1[:,q0:q1], dm0[:,q0:q1])
vj1_buf[i] += np.einsum('xp,pij->xij', wj1, coef3c)
if with_k:
rhok0_PlJ = lib.einsum('plj,Jj->plJ', rhok0a_Pl_[p0:p1], mocca)
vk1a_buf += lib.einsum('xijp,plj->xil', int3c_ip1, rhok0_PlJ[p0:p1])
rhok0_PlJ = lib.einsum('plj,Jj->plJ', rhok0b_Pl_[p0:p1], moccb)
vk1b_buf += lib.einsum('xijp,plj->xil', int3c_ip1, rhok0_PlJ[p0:p1])
int3c_ip1 = None
vj1_buf = ftmp['vj1_buf'] = vj1_buf
vk1a = vk1b = np.zeros((3,nao,nao))
for i0, ia in enumerate(atmlst):
shl0, shl1, p0, p1 = aoslices[ia]
shls_slice = (shl0, shl1, 0, nbas, 0, auxmol.nbas)
int3c_ip1 = get_int3c_ip1(shls_slice)
vj1 = -np.asarray(vj1_buf[ia])
vj1[:,p0:p1] -= np.einsum('xijp,p->xij', int3c_ip1, rhoj0_P)
if with_k:
rhok0_PlJ = lib.einsum('plj,Jj->plJ', rhok0a_Pl_, mocca[p0:p1])
vk1a = -lib.einsum('xijp,pki->xkj', int3c_ip1, rhok0_PlJ)
rhok0_PlJ = lib.einsum('plj,Jj->plJ', rhok0b_Pl_, moccb[p0:p1])
vk1b = -lib.einsum('xijp,pki->xkj', int3c_ip1, rhok0_PlJ)
vk1a[:,p0:p1] -= vk1a_buf[:,p0:p1]
vk1b[:,p0:p1] -= vk1b_buf[:,p0:p1]
if hessobj.auxbasis_response:
shl0, shl1, q0, q1 = auxslices[ia]
shls_slice = (0, nbas, 0, nbas, shl0, shl1)
int3c_ip2 = get_int3c_ip2(shls_slice)
rhoj1 = np.einsum('xijp,ji->xp', int3c_ip2, dm0)
coef3c = _load_dim0(rho0_Pij, q0, q1)
Pij = _load_dim0(wj_Pij, q0, q1)
vj1 += .5 * np.einsum('pij,xp->xij', coef3c, -rhoj1)
vj1 += .5 * np.einsum('xijp,p->xij', int3c_ip2, -rhoj0_P[q0:q1])
vj1 -= .5 * lib.einsum('xpq,q,pij->xij', int2c_ip1[:,q0:q1], -rhoj0_P, coef3c)
vj1 -= .5 * lib.einsum('pixj,p->xij', Pij, -rhoj0_P[q0:q1])
if with_k:
rhok0_PlJ = lib.einsum('plj,Jj->plJ', rhok0a_Pl_[q0:q1], mocca)
vk1a -= lib.einsum('plj,xijp->xil', rhok0_PlJ, int3c_ip2)
vk1a += lib.einsum('pjxi,plj->xil', Pij, rhok0_PlJ)
rhok0_PlJ = lib.einsum('plj,Jj->plJ', rhok0b_Pl_[q0:q1], moccb)
vk1b -= lib.einsum('plj,xijp->xil', rhok0_PlJ, int3c_ip2)
vk1b += lib.einsum('pjxi,plj->xil', Pij, rhok0_PlJ)
vj1 = vj1 + vj1.transpose(0,2,1)
if with_k:
vk1a = vk1a + vk1a.transpose(0,2,1)
vk1b = vk1b + vk1b.transpose(0,2,1)
h1 = hcore_deriv(ia)
yield ia, h1, vj1, (vk1a, vk1b)
[docs]
class Hessian(uhf_hess.Hessian):
'''Non-relativistic UHF hessian'''
def __init__(self, mf):
uhf_hess.Hessian.__init__(self, mf)
auxbasis_response = 1
partial_hess_elec = partial_hess_elec
make_h1 = make_h1
#TODO: Insert into DF class
if __name__ == '__main__':
from pyscf import gto
from pyscf import scf
mol = gto.Mole()
mol.verbose = 0
mol.output = None
mol.atom = [
[1 , (1. , 0. , 0.000)],
[1 , (0. , 1. , 0.000)],
[1 , (0. , -1.517 , 1.177)],
[1 , (0. , 1.517 , 1.177)] ]
mol.basis = '631g'
mol.spin = 2
mol.unit = 'B'
mol.build()
mf = scf.UHF(mol).density_fit()
mf.conv_tol = 1e-14
mf.scf()
n3 = mol.natm * 3
hobj = Hessian(mf)
e2 = hobj.kernel()
ref = scf.UHF(mol).run().Hessian().kernel()
print(abs(e2-ref).max())
print(lib.finger(e2) - -0.23856667321975722)
e2 = hobj.set(auxbasis_response=2).kernel()
print(abs(e2-ref).max())
print(lib.finger(e2), - 0.72321237584876141)
