#!/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 gradients for 1-electron spin-free x2c method
Ref.
JCP 135, 084114 (2011); DOI:10.1063/1.3624397
'''
from functools import reduce
import numpy
import scipy.linalg
from pyscf import lib
from pyscf import gto
from pyscf.x2c import x2c
[docs]
def hcore_grad_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(atm_id):
h1 = get_h1_xmol(atm_id)
if contr_coeff is not None:
h1 = lib.einsum('pi,xpq,qj->xij', 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 = _get_h0_s0(mol)
e0, c0 = scipy.linalg.eigh(h0, s0)
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
s_nesc0 = s0[:nao,:nao] + reduce(numpy.dot, (x0.T, s0[nao:,nao:], x0))
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)
get_h1_etc = _gen_first_order_quantities(mol, e0, c0, x0, approx)
def hcore_deriv(ia):
h1_ao, s1_ao, e1, c1, x1, s_nesc1, R1, c_fw1 = get_h1_etc(ia)
hfw1 = lib.einsum('xpi,pj->xij', c_fw1, h0_fw_half)
hfw1 = hfw1 + hfw1.transpose(0,2,1)
hfw1+= lib.einsum('pi,xpq,qj->xij', c_fw0, h1_ao, c_fw0)
return hfw1
return hcore_deriv
def _get_h0_s0(mol):
c = lib.param.LIGHT_SPEED
s = mol.intor_symmetric('int1e_ovlp')
t = mol.intor_symmetric('int1e_kin')
v = mol.intor_symmetric('int1e_nuc')
w = mol.intor_symmetric('int1e_pnucp')
nao = s.shape[0]
n2 = nao * 2
h = numpy.zeros((n2,n2), dtype=v.dtype)
m = numpy.zeros((n2,n2), dtype=v.dtype)
h[:nao,:nao] = v
h[:nao,nao:] = t
h[nao:,:nao] = t
h[nao:,nao:] = w * (.25/c**2) - t
m[:nao,:nao] = s
m[nao:,nao:] = t * (.5/c**2)
return h, m
def _gen_h1_s1(mol):
c = lib.param.LIGHT_SPEED
s1 = mol.intor('int1e_ipovlp', comp=3)
t1 = mol.intor('int1e_ipkin', comp=3)
v1 = mol.intor('int1e_ipnuc', comp=3)
w1 = mol.intor('int1e_ippnucp', comp=3)
aoslices = mol.aoslice_by_atom()
nao = s1.shape[1]
n2 = nao * 2
def get_h1_s1(ia):
h1 = numpy.zeros((3,n2,n2), dtype=v1.dtype)
m1 = numpy.zeros((3,n2,n2), dtype=v1.dtype)
ish0, ish1, i0, i1 = aoslices[ia]
with mol.with_rinv_origin(mol.atom_coord(ia)):
z = mol.atom_charge(ia)
rinv1 = -z*mol.intor('int1e_iprinv', comp=3)
prinvp1 = -z*mol.intor('int1e_ipprinvp', comp=3)
rinv1 [:,i0:i1,:] -= v1[:,i0:i1]
prinvp1[:,i0:i1,:] -= w1[:,i0:i1]
for i in range(3):
s1cc = numpy.zeros((nao,nao))
t1cc = numpy.zeros((nao,nao))
s1cc[i0:i1,:] =-s1[i,i0:i1]
s1cc[:,i0:i1]-= s1[i,i0:i1].T
t1cc[i0:i1,:] =-t1[i,i0:i1]
t1cc[:,i0:i1]-= t1[i,i0:i1].T
v1cc = rinv1[i] + rinv1[i].T
w1cc = prinvp1[i] + prinvp1[i].T
h1[i,:nao,:nao] = v1cc
h1[i,:nao,nao:] = t1cc
h1[i,nao:,:nao] = t1cc
h1[i,nao:,nao:] = w1cc * (.25/c**2) - t1cc
m1[i,:nao,:nao] = s1cc
m1[i,nao:,nao:] = t1cc * (.5/c**2)
return h1, m1
return get_h1_s1
def _gen_first_order_quantities(mol, e0, c0, x0, approx='1E'):
c = lib.param.LIGHT_SPEED
nao = e0.size // 2
n2 = nao * 2
epq = e0[:,None] - e0
degen_mask = abs(epq) < 1e-7
epq[degen_mask] = 1e200
cl0 = c0[:nao,nao:]
# cs0 = c0[nao:,nao:]
s0 = mol.intor('int1e_ovlp')
t0 = mol.intor('int1e_kin')
t0x0 = numpy.dot(t0, x0) * (.5/c**2)
s_nesc0 = s0[:nao,:nao] + numpy.dot(x0.T, t0x0)
w_s, v_s = scipy.linalg.eigh(s0)
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))
get_h1_s1 = _gen_h1_s1(mol)
def get_first_order(ia):
h1ao, s1ao = get_h1_s1(ia)
h1mo = lib.einsum('pi,xpq,qj->xij', c0.conj(), h1ao, c0)
s1mo = lib.einsum('pi,xpq,qj->xij', c0.conj(), s1ao, c0)
if 'ATOM' in approx:
e1 = c1_ao = x1 = None
s_nesc1 = lib.einsum('pi,xpq,qj->xij', x0, s1ao[:,nao:,nao:], x0)
s_nesc1+= s1ao[:,:nao,:nao]
else:
f1 = h1mo[:,:,nao:] - s1mo[:,:,nao:] * e0[nao:]
c1 = f1 / -epq[:,nao:]
e1 = f1[:,nao:]
e1[:,~degen_mask[nao:,nao:]] = 0
c1_ao = lib.einsum('pq,xqi->xpi', c0, c1)
cl1 = c1_ao[:,:nao]
cs1 = c1_ao[:,nao:]
tmp = cs1 - lib.einsum('pq,xqi->xpi', x0, cl1)
x1 = scipy.linalg.solve(cl0.T, tmp.reshape(-1,nao).T)
x1 = x1.T.reshape(3,nao,nao)
s_nesc1 = lib.einsum('xpi,pj->xij', x1, t0x0)
s_nesc1 = s_nesc1 + s_nesc1.transpose(0,2,1)
s_nesc1+= lib.einsum('pi,xpq,qj->xij', x0, s1ao[:,nao:,nao:], x0)
s_nesc1+= s1ao[:,:nao,:nao]
R1 = numpy.empty((3,nao,nao))
c_fw1 = numpy.empty((3,n2,nao))
for i in range(3):
R1[i] = _get_r1((w_sqrt,v_s), s_nesc0_vbas,
s1ao[i,:nao,:nao], s_nesc1[i], (wr0_sqrt,vr0))
c_fw1[i,:nao] = R1[i]
c_fw1[i,nao:] = numpy.dot(x0, R1[i])
if 'ATOM' not in approx:
c_fw1[i,nao:] += numpy.dot(x1[i], R0)
return h1ao, s1ao, e1, c1_ao, x1, s_nesc1, R1, c_fw1
return get_first_order
def _get_r1(s0_roots, s_nesc0, s1, s_nesc1, r0_roots):
# See JCP 135, 084114 (2011); DOI:10.1063/1.3624397, Eq (34)
w_sqrt, v_s = s0_roots
w_invsqrt = 1. / w_sqrt
wr0_sqrt, vr0 = r0_roots
wr0_invsqrt = 1. / wr0_sqrt
s1 = reduce(numpy.dot, (v_s.T, s1, v_s))
s_nesc1 = reduce(numpy.dot, (v_s.T, s_nesc1, v_s))
s1_sqrt = s1 / (w_sqrt[:,None] + w_sqrt)
s1_invsqrt = (numpy.einsum('i,ij,j->ij', w_invsqrt**2, s1, w_invsqrt**2)
/ -(w_invsqrt[:,None] + w_invsqrt))
R1_mid = numpy.dot(s1_invsqrt, s_nesc0) * w_invsqrt
R1_mid = R1_mid + R1_mid.T
R1_mid += numpy.einsum('i,ij,j->ij', w_invsqrt, s_nesc1, w_invsqrt)
R1_mid = reduce(numpy.dot, (vr0.T, R1_mid, vr0))
R1_mid /= -(wr0_invsqrt[:,None] + wr0_invsqrt)
R1_mid = numpy.einsum('i,ij,j->ij', wr0_invsqrt**2, R1_mid, wr0_invsqrt**2)
vr0_wr0_sqrt = vr0 * wr0_invsqrt
vr0_s0_sqrt = vr0.T * w_sqrt
vr0_s0_invsqrt = vr0.T * w_invsqrt
R1 = reduce(numpy.dot, (vr0_s0_invsqrt.T, R1_mid, vr0_s0_sqrt))
R1 += reduce(numpy.dot, (s1_invsqrt, vr0_wr0_sqrt, vr0_s0_sqrt))
R1 += reduce(numpy.dot, (vr0_s0_invsqrt.T, vr0_wr0_sqrt.T, s1_sqrt))
R1 = reduce(numpy.dot, (v_s, R1, v_s.T))
return R1
if __name__ == '__main__':
bak = lib.param.LIGHT_SPEED
lib.param.LIGHT_SPEED = 10
def get_h(mol):
c = lib.param.LIGHT_SPEED
t = mol.intor_symmetric('int1e_kin')
v = mol.intor_symmetric('int1e_nuc')
s = mol.intor_symmetric('int1e_ovlp')
w = mol.intor_symmetric('int1e_pnucp')
return x2c._x2c1e_get_hcore(t, v, w, s, c)
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',
)
h_1 = get_h(mol)
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',
)
h_2 = get_h(mol)
h_ref = (h_1 - h_2) / 0.0002 * lib.param.BOHR
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',
)
hcore_deriv = gen_sf_hfw(mol)
h1 = hcore_deriv(0)
print(abs(h1[2]-h_ref).max())
lib.param.LIGHT_SPEED = bak
print(lib.finger(h1) - -1.4618392662849411)
hcore_deriv = gen_sf_hfw(mol, approx='atom1e')
h1 = hcore_deriv(0)
print(lib.finger(h1) - -1.3596826558976405)
