Source code for pyscf.x2c.sfx2c1e_hess

#!/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.

'''
Analytical nuclear hessian for 1-electron spin-free x2c method

Ref.
JCP 135, 244104 (2011); DOI:10.1063/1.3667202
JCTC 8, 2617 (2012); DOI:10.1021/ct300127e
'''

from functools import reduce
import numpy
import scipy.linalg
from pyscf import lib
from pyscf import gto
from pyscf.x2c import x2c
from pyscf.x2c import sfx2c1e_grad

[docs]def hcore_hess_generator(x2cobj, mol=None): '''nuclear gradients of 1-component X2c hcore Hamiltonian (spin-free part only) ''' if mol is None: mol = x2cobj.mol xmol, contr_coeff = x2cobj.get_xmol(mol) if x2cobj.basis is not None: s22 = xmol.intor_symmetric('int1e_ovlp') s21 = gto.intor_cross('int1e_ovlp', xmol, mol) contr_coeff = lib.cho_solve(s22, s21) get_h1_xmol = gen_sf_hfw(xmol, x2cobj.approx) def hcore_deriv(ia, ja): h1 = get_h1_xmol(ia, ja) if contr_coeff is not None: h1 = lib.einsum('pi,xypq,qj->xyij', contr_coeff, h1, contr_coeff) return numpy.asarray(h1) return hcore_deriv
[docs]def gen_sf_hfw(mol, approx='1E'): approx = approx.upper() c = lib.param.LIGHT_SPEED h0, s0 = sfx2c1e_grad._get_h0_s0(mol) e0, c0 = scipy.linalg.eigh(h0, s0) c0[:,c0[1]<0] *= -1 aoslices = mol.aoslice_by_atom() nao = mol.nao_nr() if 'ATOM' in approx: x0 = numpy.zeros((nao,nao)) for ia in range(mol.natm): ish0, ish1, p0, p1 = aoslices[ia] shls_slice = (ish0, ish1, ish0, ish1) t1 = mol.intor('int1e_kin', shls_slice=shls_slice) s1 = mol.intor('int1e_ovlp', shls_slice=shls_slice) with mol.with_rinv_at_nucleus(ia): z = -mol.atom_charge(ia) v1 = z * mol.intor('int1e_rinv', shls_slice=shls_slice) w1 = z * mol.intor('int1e_prinvp', shls_slice=shls_slice) x0[p0:p1,p0:p1] = x2c._x2c1e_xmatrix(t1, v1, w1, s1, c) else: cl0 = c0[:nao,nao:] cs0 = c0[nao:,nao:] x0 = scipy.linalg.solve(cl0.T, cs0.T).T t0x0 = numpy.dot(s0[nao:,nao:], x0) s_nesc0 = s0[:nao,:nao] + numpy.dot(x0.T, t0x0) w_s, v_s = scipy.linalg.eigh(s0[:nao,:nao]) w_sqrt = numpy.sqrt(w_s) s_nesc0_vbas = reduce(numpy.dot, (v_s.T, s_nesc0, v_s)) R0_mid = numpy.einsum('i,ij,j->ij', 1./w_sqrt, s_nesc0_vbas, 1./w_sqrt) wr0, vr0 = scipy.linalg.eigh(R0_mid) wr0_sqrt = numpy.sqrt(wr0) # R0 in v_s basis R0 = numpy.dot(vr0/wr0_sqrt, vr0.T) R0 *= w_sqrt R0 /= w_sqrt[:,None] # Transform R0 back R0 = reduce(numpy.dot, (v_s, R0, v_s.T)) R0 = x2c._get_r(s0[:nao,:nao], s_nesc0) c_fw0 = numpy.vstack((R0, numpy.dot(x0, R0))) h0_fw_half = numpy.dot(h0, c_fw0) epq = e0[:,None] - e0 degen_mask = abs(epq) < 1e-7 epq[degen_mask] = 1e200 s2aa = mol.intor('int1e_ipipovlp', comp=9).reshape(3,3,nao,nao) t2aa = mol.intor('int1e_ipipkin', comp=9).reshape(3,3,nao,nao) v2aa = mol.intor('int1e_ipipnuc', comp=9).reshape(3,3,nao,nao) w2aa = mol.intor('int1e_ipippnucp', comp=9).reshape(3,3,nao,nao) s2ab = mol.intor('int1e_ipovlpip', comp=9).reshape(3,3,nao,nao) t2ab = mol.intor('int1e_ipkinip', comp=9).reshape(3,3,nao,nao) v2ab = mol.intor('int1e_ipnucip', comp=9).reshape(3,3,nao,nao) w2ab = mol.intor('int1e_ippnucpip', comp=9).reshape(3,3,nao,nao) n2 = nao * 2 h2ao = numpy.zeros((3,3,n2,n2), dtype=v2aa.dtype) s2ao = numpy.zeros((3,3,n2,n2), dtype=v2aa.dtype) get_h1_etc = sfx2c1e_grad._gen_first_order_quantities(mol, e0, c0, x0, approx) def hcore_deriv(ia, ja): ish0, ish1, i0, i1 = aoslices[ia] jsh0, jsh1, j0, j1 = aoslices[ja] s2cc = numpy.zeros_like(s2aa) t2cc = numpy.zeros_like(s2aa) v2cc = numpy.zeros_like(s2aa) w2cc = numpy.zeros_like(s2aa) if ia == ja: with mol.with_rinv_origin(mol.atom_coord(ia)): z = mol.atom_charge(ia) rinv2aa = z*mol.intor('int1e_ipiprinv', comp=9).reshape(3,3,nao,nao) rinv2ab = z*mol.intor('int1e_iprinvip', comp=9).reshape(3,3,nao,nao) prinvp2aa = z*mol.intor('int1e_ipipprinvp', comp=9).reshape(3,3,nao,nao) prinvp2ab = z*mol.intor('int1e_ipprinvpip', comp=9).reshape(3,3,nao,nao) s2cc[:,:,i0:i1 ] = s2aa[:,:,i0:i1 ] s2cc[:,:,i0:i1,j0:j1]+= s2ab[:,:,i0:i1,j0:j1] t2cc[:,:,i0:i1 ] = t2aa[:,:,i0:i1 ] t2cc[:,:,i0:i1,j0:j1]+= t2ab[:,:,i0:i1,j0:j1] v2cc -= rinv2aa + rinv2ab v2cc[:,:,i0:i1 ]+= v2aa[:,:,i0:i1 ] v2cc[:,:,i0:i1,j0:j1]+= v2ab[:,:,i0:i1,j0:j1] v2cc[:,:,i0:i1 ]+= rinv2aa[:,:,i0:i1] v2cc[:,:,i0:i1 ]+= rinv2ab[:,:,i0:i1] v2cc[:,:,: ,i0:i1]+= rinv2aa[:,:,i0:i1].transpose(0,1,3,2) v2cc[:,:,: ,i0:i1]+= rinv2ab[:,:,:,i0:i1] w2cc -= prinvp2aa + prinvp2ab w2cc[:,:,i0:i1 ]+= w2aa[:,:,i0:i1 ] w2cc[:,:,i0:i1,j0:j1]+= w2ab[:,:,i0:i1,j0:j1] w2cc[:,:,i0:i1 ]+= prinvp2aa[:,:,i0:i1] w2cc[:,:,i0:i1 ]+= prinvp2ab[:,:,i0:i1] w2cc[:,:,: ,i0:i1]+= prinvp2aa[:,:,i0:i1].transpose(0,1,3,2) w2cc[:,:,: ,i0:i1]+= prinvp2ab[:,:,:,i0:i1] else: s2cc[:,:,i0:i1,j0:j1] = s2ab[:,:,i0:i1,j0:j1] t2cc[:,:,i0:i1,j0:j1] = t2ab[:,:,i0:i1,j0:j1] v2cc[:,:,i0:i1,j0:j1] = v2ab[:,:,i0:i1,j0:j1] w2cc[:,:,i0:i1,j0:j1] = w2ab[:,:,i0:i1,j0:j1] zi = mol.atom_charge(ia) zj = mol.atom_charge(ja) with mol.with_rinv_at_nucleus(ia): shls_slice = (jsh0, jsh1, 0, mol.nbas) rinv2aa = mol.intor('int1e_ipiprinv', comp=9, shls_slice=shls_slice) rinv2ab = mol.intor('int1e_iprinvip', comp=9, shls_slice=shls_slice) prinvp2aa = mol.intor('int1e_ipipprinvp', comp=9, shls_slice=shls_slice) prinvp2ab = mol.intor('int1e_ipprinvpip', comp=9, shls_slice=shls_slice) rinv2aa = zi * rinv2aa.reshape(3,3,j1-j0,nao) rinv2ab = zi * rinv2ab.reshape(3,3,j1-j0,nao) prinvp2aa = zi * prinvp2aa.reshape(3,3,j1-j0,nao) prinvp2ab = zi * prinvp2ab.reshape(3,3,j1-j0,nao) v2cc[:,:,j0:j1] += rinv2aa v2cc[:,:,j0:j1] += rinv2ab.transpose(1,0,2,3) w2cc[:,:,j0:j1] += prinvp2aa w2cc[:,:,j0:j1] += prinvp2ab.transpose(1,0,2,3) with mol.with_rinv_at_nucleus(ja): shls_slice = (ish0, ish1, 0, mol.nbas) rinv2aa = mol.intor('int1e_ipiprinv', comp=9, shls_slice=shls_slice) rinv2ab = mol.intor('int1e_iprinvip', comp=9, shls_slice=shls_slice) prinvp2aa = mol.intor('int1e_ipipprinvp', comp=9, shls_slice=shls_slice) prinvp2ab = mol.intor('int1e_ipprinvpip', comp=9, shls_slice=shls_slice) rinv2aa = zj * rinv2aa.reshape(3,3,i1-i0,nao) rinv2ab = zj * rinv2ab.reshape(3,3,i1-i0,nao) prinvp2aa = zj * prinvp2aa.reshape(3,3,i1-i0,nao) prinvp2ab = zj * prinvp2ab.reshape(3,3,i1-i0,nao) v2cc[:,:,i0:i1] += rinv2aa v2cc[:,:,i0:i1] += rinv2ab w2cc[:,:,i0:i1] += prinvp2aa w2cc[:,:,i0:i1] += prinvp2ab s2cc = s2cc + s2cc.transpose(0,1,3,2) t2cc = t2cc + t2cc.transpose(0,1,3,2) v2cc = v2cc + v2cc.transpose(0,1,3,2) w2cc = w2cc + w2cc.transpose(0,1,3,2) h2ao[:,:,:nao,:nao] = v2cc h2ao[:,:,:nao,nao:] = t2cc h2ao[:,:,nao:,:nao] = t2cc h2ao[:,:,nao:,nao:] = w2cc * (.25/c**2) - t2cc s2ao[:,:,:nao,:nao] = s2cc s2ao[:,:,nao:,nao:] = t2cc * (.5/c**2) h1i, s1i, e1i, c1i, x1i, s_nesc1i, R1i, c_fw1i = get_h1_etc(ia) h1j, s1j, e1j, c1j, x1j, s_nesc1j, R1j, c_fw1j = get_h1_etc(ja) if 'ATOM' not in approx: f2 = lib.einsum('xypq,qj->xypj', h2ao, c0[:,nao:]) f2+= lib.einsum('xpq,yqj->xypj', h1i, c1j) f2+= lib.einsum('ypq,xqj->xypj', h1j, c1i) sc2 = lib.einsum('xypq,qj->xypj', s2ao, c0[:,nao:]) sc2+= lib.einsum('xpq,yqj->xypj', s1i, c1j) sc2+= lib.einsum('ypq,xqj->xypj', s1j, c1i) f2-= sc2 * e0[nao:] sc1i = lib.einsum('xpq,qj->xpj', s1i, c0[:,nao:]) sc1j = lib.einsum('xpq,qj->xpj', s1j, c0[:,nao:]) sc1i+= lib.einsum('pq,xqj->xpj', s0, c1i) sc1j+= lib.einsum('pq,xqj->xpj', s0, c1j) f2-= lib.einsum('xpq,yqj->xypj', sc1i, e1j) f2-= lib.einsum('ypq,xqj->xypj', sc1j, e1i) c2 = lib.einsum('pi,xypj->xyij', c0.conj(), f2) / -epq[:,nao:] c2_ao = lib.einsum('pq,xyqi->xypi', c0, c2) cl2 = c2_ao[:,:,:nao] cs2 = c2_ao[:,:,nao:] tmp = cs2 - lib.einsum('pq,xyqi->xypi', x0, cl2) tmp-= lib.einsum('xpq,yqi->xypi', x1i, c1j[:,:nao]) tmp-= lib.einsum('ypq,xqi->xypi', x1j, c1i[:,:nao]) x2 = scipy.linalg.solve(cl0.T, tmp.reshape(-1,nao).T).T.reshape(3,3,nao,nao) hfw2 = numpy.empty((3,3,nao,nao)) for i in range(3): for j in range(3): if 'ATOM' in approx: s_nesc2 = reduce(numpy.dot, (x0.T, s2ao[i,j,nao:,nao:], x0)) s_nesc2 += s2ao[i,j,:nao,:nao] R2 = _get_r2((w_sqrt,v_s), s_nesc0, s1i[i,:nao,:nao], s_nesc1i[i], s1j[j,:nao,:nao], s_nesc1j[j], s2ao[i,j,:nao,:nao], s_nesc2, (wr0_sqrt,vr0)) c_fw2 = numpy.vstack((R2, numpy.dot(x0, R2))) else: s_nesc2 = numpy.dot(x2[i,j].T, t0x0) s_nesc2 += reduce(numpy.dot, (x1i[i].T, s1j[j,nao:,nao:], x0)) s_nesc2 += reduce(numpy.dot, (x0.T, s1i[i,nao:,nao:], x1j[j])) s_nesc2 += reduce(numpy.dot, (x1i[i].T, s0[nao:,nao:], x1j[j])) s_nesc2 = s_nesc2 + s_nesc2.T s_nesc2 += reduce(numpy.dot, (x0.T, s2ao[i,j,nao:,nao:], x0)) s_nesc2 += s2ao[i,j,:nao,:nao] R2 = _get_r2((w_sqrt,v_s), s_nesc0, s1i[i,:nao,:nao], s_nesc1i[i], s1j[j,:nao,:nao], s_nesc1j[j], s2ao[i,j,:nao,:nao], s_nesc2, (wr0_sqrt,vr0)) c_fw_s = (numpy.dot(x0, R2) + numpy.dot(x1i[i], R1j[j]) + numpy.dot(x1j[j], R1i[i]) + numpy.dot(x2[i,j], R0)) c_fw2 = numpy.vstack((R2, c_fw_s)) tmp = numpy.dot(c_fw2.T, h0_fw_half) tmp += reduce(numpy.dot, (c_fw1i[i].T, h1j[j], c_fw0)) tmp += reduce(numpy.dot, (c_fw0.T, h1i[i], c_fw1j[j])) tmp += reduce(numpy.dot, (c_fw1i[i].T, h0, c_fw1j[j])) hfw2[i,j] = tmp + tmp.T hfw2[i,j]+= reduce(numpy.dot, (c_fw0.T, h2ao[i,j], c_fw0)) return hfw2 return hcore_deriv
def _get_r2(s0_roots, sa0, s1i, sa1i, s1j, sa1j, s2, sa2, r0_roots): w_sqrt, v_s = s0_roots w_invsqrt = 1. / w_sqrt wr0_sqrt, vr0 = r0_roots wr0_invsqrt = 1. / wr0_sqrt sa0 = lib.einsum('pi,pq,qj->ij', v_s, sa0 , v_s) s1i = lib.einsum('pi,pq,qj->ij', v_s, s1i , v_s) s1j = lib.einsum('pi,pq,qj->ij', v_s, s1j , v_s) s2 = lib.einsum('pi,pq,qj->ij', v_s, s2 , v_s) sa1i = lib.einsum('pi,pq,qj->ij', v_s, sa1i, v_s) sa1j = lib.einsum('pi,pq,qj->ij', v_s, sa1j, v_s) sa2 = lib.einsum('pi,pq,qj->ij', v_s, sa2 , v_s) s1i_sqrt = s1i / (w_sqrt[:,None] + w_sqrt) s1i_invsqrt = (numpy.einsum('i,ij,j->ij', w_invsqrt**2, s1i, w_invsqrt**2) / -(w_invsqrt[:,None] + w_invsqrt)) s1j_sqrt = s1j / (w_sqrt[:,None] + w_sqrt) s1j_invsqrt = (numpy.einsum('i,ij,j->ij', w_invsqrt**2, s1j, w_invsqrt**2) / -(w_invsqrt[:,None] + w_invsqrt)) tmp = numpy.dot(s1i_sqrt, s1j_sqrt) s2_sqrt = (s2 - tmp - tmp.T) / (w_sqrt[:,None] + w_sqrt) tmp = numpy.dot(s1i*w_invsqrt**2, s1j) tmp = s2 - tmp - tmp.T tmp = -numpy.einsum('i,ij,j->ij', w_invsqrt**2, tmp, w_invsqrt**2) tmp1 = numpy.dot(s1i_invsqrt, s1j_invsqrt) s2_invsqrt = (tmp - tmp1 - tmp1.T) / (w_invsqrt[:,None] + w_invsqrt) R1i_mid = lib.einsum('ip,pj,j->ij', s1i_invsqrt, sa0, w_invsqrt) R1i_mid = R1i_mid + R1i_mid.T R1i_mid+= numpy.einsum('i,ij,j->ij', w_invsqrt, sa1i, w_invsqrt) R1i_mid = tmpi = lib.einsum('pi,pq,qj->ij', vr0, R1i_mid, vr0) R1i_mid = (numpy.einsum('i,ij,j->ij', wr0_invsqrt**2, R1i_mid, wr0_invsqrt**2) / -(wr0_invsqrt[:,None] + wr0_invsqrt)) R1j_mid = lib.einsum('ip,pj,j->ij', s1j_invsqrt, sa0, w_invsqrt) R1j_mid = R1j_mid + R1j_mid.T R1j_mid+= numpy.einsum('i,ij,j->ij', w_invsqrt, sa1j, w_invsqrt) R1j_mid = tmpj = lib.einsum('pi,pq,qj->ij', vr0, R1j_mid, vr0) R1j_mid = (numpy.einsum('i,ij,j->ij', wr0_invsqrt**2, R1j_mid, wr0_invsqrt**2) / -(wr0_invsqrt[:,None] + wr0_invsqrt)) # second derivative of (s_invsqrt * sa * s_invsqrt), 9 terms R2_mid = lib.einsum('ip,pj,j->ij', s2_invsqrt , sa0 , w_invsqrt) R2_mid+= lib.einsum('ip,pj,j->ij', s1i_invsqrt, sa1j, w_invsqrt) R2_mid+= lib.einsum('i,ip,pj->ij', w_invsqrt , sa1i, s1j_invsqrt) R2_mid+= lib.einsum('ip,pq,qj->ij', s1i_invsqrt, sa0 , s1j_invsqrt) R2_mid = R2_mid + R2_mid.T R2_mid+= numpy.einsum('i,ij,j->ij', w_invsqrt, sa2, w_invsqrt) R2_mid = lib.einsum('pi,pq,qj->ij', vr0, R2_mid, vr0) tmp = numpy.dot(tmpi*wr0_invsqrt**2, tmpj) tmp = R2_mid - tmp - tmp.T tmp = -numpy.einsum('i,ij,j->ij', wr0_invsqrt**2, tmp, wr0_invsqrt**2) tmp1 = numpy.dot(R1i_mid, R1j_mid) R2_mid = (tmp - tmp1 - tmp1.T) / (wr0_invsqrt[:,None] + wr0_invsqrt) R0_mid = numpy.dot(vr0*wr0_invsqrt, vr0.T) R1i_mid = reduce(numpy.dot, (vr0, R1i_mid, vr0.T)) R1j_mid = reduce(numpy.dot, (vr0, R1j_mid, vr0.T)) R2_mid = reduce(numpy.dot, (vr0, R2_mid, vr0.T)) R2 = lib.einsum('ip,pj,j->ij' , s2_invsqrt , R0_mid , w_sqrt) R2 += lib.einsum('ip,pj,j->ij' , s1i_invsqrt, R1j_mid, w_sqrt) R2 += lib.einsum('ip,pq,qj->ij', s1i_invsqrt, R0_mid , s1j_sqrt) R2 += lib.einsum('ip,pj,j->ij' , s1j_invsqrt, R1i_mid, w_sqrt) R2 += numpy.einsum('i,ij,j->ij', w_invsqrt , R2_mid , w_sqrt) R2 += lib.einsum('i,iq,qj->ij' , w_invsqrt , R1i_mid, s1j_sqrt) R2 += lib.einsum('ip,pq,qj->ij', s1j_invsqrt, R0_mid , s1i_sqrt) R2 += lib.einsum('i,iq,qj->ij' , w_invsqrt , R1j_mid, s1i_sqrt) R2 += lib.einsum('i,iq,qj->ij' , w_invsqrt , R0_mid , s2_sqrt) R2 = reduce(numpy.dot, (v_s, R2, v_s.T)) return R2 if __name__ == '__main__': bak = lib.param.LIGHT_SPEED lib.param.LIGHT_SPEED = 10 mol = gto.M( verbose = 0, atom = [["O" , (0. , 0. , 0.0001)], [1 , (0. , -0.757 , 0.587)], [1 , (0. , 0.757 , 0.587)]], basis = '3-21g', ) h1_deriv_1 = sfx2c1e_grad.gen_sf_hfw(mol, approx='1E') mol = gto.M( verbose = 0, atom = [["O" , (0. , 0. ,-0.0001)], [1 , (0. , -0.757 , 0.587)], [1 , (0. , 0.757 , 0.587)]], basis = '3-21g', ) h1_deriv_2 = sfx2c1e_grad.gen_sf_hfw(mol, approx='1E') mol = gto.M( verbose = 0, atom = [["O" , (0. , 0. , 0. )], [1 , (0. , -0.757 , 0.587)], [1 , (0. , 0.757 , 0.587)]], basis = '3-21g', ) h2_deriv = gen_sf_hfw(mol) h2 = h2_deriv(0,0) h2_ref = (h1_deriv_1(0)[2] - h1_deriv_2(0)[2]) / 0.0002 * lib.param.BOHR print(abs(h2[2,2]-h2_ref).max()) print(lib.finger(h2) - 33.71188112440316) h2 = h2_deriv(1,0) h2_ref = (h1_deriv_1(1)[2] - h1_deriv_2(1)[2]) / 0.0002 * lib.param.BOHR print(abs(h2[2,2]-h2_ref).max()) print(lib.finger(h2) - -23.609411428378138) lib.param.LIGHT_SPEED = bak