#!/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.
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#
import numpy
from pyscf import gto
from pyscf.df import incore
# (ij|L)
[docs]
def aux_e2(mol, auxmol, intor='int3c2e_spinor', aosym='s1', comp=None, hermi=0):
intor, comp = gto.moleintor._get_intor_and_comp(mol._add_suffix(intor), comp)
atm, bas, env = gto.mole.conc_env(mol._atm, mol._bas, mol._env,
auxmol._atm, auxmol._bas, auxmol._env)
shls_slice = (0, mol.nbas, 0, mol.nbas, mol.nbas, mol.nbas+auxmol.nbas)
ao_loc1 = mol.ao_loc_2c()
ao_loc2 = auxmol.ao_loc_nr('ssc' in intor)
nao = ao_loc1[-1]
ao_loc = numpy.append(ao_loc1, ao_loc2[1:]+nao)
out = gto.moleintor.getints3c(intor, atm, bas, env, shls_slice,
comp, aosym, ao_loc=ao_loc)
return out
# (L|ij)
[docs]
def aux_e1(mol, auxmol, intor='int3c2e_spinor', aosym='s1', comp=1, hermi=0):
raise NotImplementedError
[docs]
def cholesky_eri(mol, auxbasis='weigend+etb', auxmol=None,
int3c='int3c2e_spinor', aosym='s1', int2c='int2c2e_sph', comp=1,
verbose=0):
return incore.cholesky_eri_debug(mol, auxbasis, auxmol, int3c, aosym, int2c,
comp, verbose, aux_e2)
if __name__ == '__main__':
from pyscf import lib
from pyscf import scf
mol = gto.Mole()
mol.build(
verbose = 0,
atom = [["O" , (0. , 0. , 0.)],
[1 , (0. , -0.757 , 0.587)],
[1 , (0. , 0.757 , 0.587)] ],
basis = 'ccpvdz',
)
cderi = (cholesky_eri(mol, int3c='int3c2e_spinor', verbose=5),
cholesky_eri(mol, int3c='int3c2e_spsp1_spinor', verbose=5))
n2c = mol.nao_2c()
c2 = .5 / lib.param.LIGHT_SPEED
def fjk(mol, dm, *args, **kwargs):
# dm is 4C density matrix
cderi_ll = cderi[0].reshape(-1,n2c,n2c)
cderi_ss = cderi[1].reshape(-1,n2c,n2c)
vj = numpy.zeros((n2c*2,n2c*2), dtype=dm.dtype)
vk = numpy.zeros((n2c*2,n2c*2), dtype=dm.dtype)
rho =(numpy.dot(cderi[0], dm[:n2c,:n2c].T.reshape(-1)) +
numpy.dot(cderi[1], dm[n2c:,n2c:].T.reshape(-1)*c2**2))
vj[:n2c,:n2c] = numpy.dot(rho, cderi[0]).reshape(n2c,n2c)
vj[n2c:,n2c:] = numpy.dot(rho, cderi[1]).reshape(n2c,n2c) * c2**2
v1 = numpy.einsum('pij,jk->pik', cderi_ll, dm[:n2c,:n2c])
vk[:n2c,:n2c] = numpy.einsum('pik,pkj->ij', v1, cderi_ll)
v1 = numpy.einsum('pij,jk->pik', cderi_ss, dm[n2c:,n2c:])
vk[n2c:,n2c:] = numpy.einsum('pik,pkj->ij', v1, cderi_ss) * c2**4
v1 = numpy.einsum('pij,jk->pik', cderi_ll, dm[:n2c,n2c:])
vk[:n2c,n2c:] = numpy.einsum('pik,pkj->ij', v1, cderi_ss) * c2**2
vk[n2c:,:n2c] = vk[:n2c,n2c:].T.conj()
return vj, vk
mf = scf.DHF(mol)
mf.get_jk = fjk
mf.direct_scf = False
ehf1 = mf.scf()
print(ehf1, -76.08073868516945)
cderi = cderi[0].reshape(-1,n2c,n2c)
print(numpy.allclose(cderi, cderi.transpose(0,2,1).conj()))
