#!/usr/bin/env python
# Copyright 2014-2018 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>
#
'''
Effective core potential (ECP)
This module exposes some ecp integration functions from the C implementation.
Reference for ecp integral computation
* Analytical integration
J. Chem. Phys. 65, 3826 (1976); https://doi.org/10.1063/1.432900
J. Chem. Phys. 111, 8778 (1999); https://doi.org/10.1063/1.480225
J. Comput. Phys. 44, 289 (1981); https://doi.org/10.1016/0021-9991(81)90053-X
* Numerical integration
J. Comput. Chem. 27, 1009 (2006); https://doi.org/10.1002/jcc.20410
Chem. Phys. Lett. 296, 445 (1998); https://doi.org/10.1016/S0009-2614(98)01077-X
'''
import ctypes
import numpy
from pyscf import lib
from pyscf.gto import moleintor
from pyscf.data.elements import ELEMENTS
libecp = moleintor.libcgto
libecp.ECPscalar_cache_size.restype = ctypes.c_int
[docs]
def type1_by_shell(mol, shls, cart=False):
li = mol.bas_angular(shls[0])
lj = mol.bas_angular(shls[1])
if cart:
fn = libecp.ECPtype1_cart
di = (li+1)*(li+2)//2 * mol.bas_nctr(shls[0])
dj = (lj+1)*(lj+2)//2 * mol.bas_nctr(shls[1])
else:
fn = libecp.ECPtype1_sph
di = (li*2+1) * mol.bas_nctr(shls[0])
dj = (lj*2+1) * mol.bas_nctr(shls[1])
cache_size = libecp.ECPscalar_cache_size(
ctypes.c_int(1), (ctypes.c_int*2)(*shls),
mol._atm.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.natm),
mol._bas.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.nbas),
mol._env.ctypes.data_as(ctypes.c_void_p))
cache = numpy.empty(cache_size)
buf = numpy.zeros((di,dj), order='F')
fn(buf.ctypes.data_as(ctypes.c_void_p),
(ctypes.c_int*2)(*shls),
mol._ecpbas.ctypes.data_as(ctypes.c_void_p),
ctypes.c_int(len(mol._ecpbas)),
mol._atm.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.natm),
mol._bas.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.nbas),
mol._env.ctypes.data_as(ctypes.c_void_p), lib.c_null_ptr(),
cache.ctypes.data_as(ctypes.c_void_p))
return buf
[docs]
def type2_by_shell(mol, shls, cart=False):
li = mol.bas_angular(shls[0])
lj = mol.bas_angular(shls[1])
if cart:
fn = libecp.ECPtype2_cart
di = (li+1)*(li+2)//2 * mol.bas_nctr(shls[0])
dj = (lj+1)*(lj+2)//2 * mol.bas_nctr(shls[1])
else:
fn = libecp.ECPtype2_sph
di = (li*2+1) * mol.bas_nctr(shls[0])
dj = (lj*2+1) * mol.bas_nctr(shls[1])
cache_size = libecp.ECPscalar_cache_size(
ctypes.c_int(1), (ctypes.c_int*2)(*shls),
mol._atm.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.natm),
mol._bas.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.nbas),
mol._env.ctypes.data_as(ctypes.c_void_p))
cache = numpy.empty(cache_size)
buf = numpy.zeros((di,dj), order='F')
fn(buf.ctypes.data_as(ctypes.c_void_p),
(ctypes.c_int*2)(*shls),
mol._ecpbas.ctypes.data_as(ctypes.c_void_p),
ctypes.c_int(len(mol._ecpbas)),
mol._atm.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.natm),
mol._bas.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.nbas),
mol._env.ctypes.data_as(ctypes.c_void_p), lib.c_null_ptr(),
cache.ctypes.data_as(ctypes.c_void_p))
return buf
AS_ECPBAS_OFFSET= 18
AS_NECPBAS = 19
[docs]
def so_by_shell(mol, shls):
'''Spin-orbit coupling ECP in spinor basis
i/2 <Pauli_matrix dot l U(r)>
'''
li = mol.bas_angular(shls[0])
lj = mol.bas_angular(shls[1])
di = (li*4+2) * mol.bas_nctr(shls[0])
dj = (lj*4+2) * mol.bas_nctr(shls[1])
bas = numpy.vstack((mol._bas, mol._ecpbas))
mol._env[AS_ECPBAS_OFFSET] = len(mol._bas)
mol._env[AS_NECPBAS] = len(mol._ecpbas)
buf = numpy.zeros((di,dj), order='F', dtype=numpy.complex128)
cache = numpy.empty(buf.size*48+100000)
fn = libecp.ECPso_spinor
fn(buf.ctypes.data_as(ctypes.c_void_p),
(ctypes.c_int*2)(di, dj),
(ctypes.c_int*2)(*shls),
mol._atm.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.natm),
bas.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.nbas),
mol._env.ctypes.data_as(ctypes.c_void_p), lib.c_null_ptr(),
cache.ctypes.data_as(ctypes.c_void_p))
return buf
[docs]
def core_configuration(nelec_core, atom_symbol=None):
conf_dic = {
0 : '0s0p0d0f',
2 : '1s0p0d0f',
10: '2s1p0d0f',
18: '3s2p0d0f',
28: '3s2p1d0f',
36: '4s3p1d0f',
46: '4s3p2d0f',
54: '5s4p2d0f',
60: '4s3p2d1f',
68: '5s4p2d1f',
78: '5s4p3d1f',
92: '5s4p3d2f',
}
# Core configurations for f-in-core ECPs defined in the following references
# 10.1007/BF00528565 , 10.1007/s00214-005-0629-0 , 10.1007/s00214-009-0584-2
elements_4f = ELEMENTS[57:71]
elements_5f = ELEMENTS[89:103]
if atom_symbol in elements_4f:
# TODO: fix La3+ and Ce4+ ECP46MWB f-in-core ECPs
for i in range(47, 59):
conf_dic[i] = '4s3p2d1f'
if atom_symbol in elements_5f:
# TODO: fix Ac3+, Th4+, Pa5+ and U6+ ECP78MWB f-in-core ECPs
for i in range(79, 91):
conf_dic[i] = '5s4p3d2f'
if nelec_core not in conf_dic:
raise RuntimeError('Core configuration for %d core electrons is not available.' % nelec_core)
coreshell = [int(x) for x in conf_dic[nelec_core][::2]]
return coreshell
if __name__ == '__main__':
from pyscf import gto, scf
mol = gto.M(atom='''
Cu 0. 0. 0.
H 0. 0. -1.56
H 0. 0. 1.56
''',
basis={'Cu':'lanl2dz', 'H':'sto3g'},
ecp = {'cu':'lanl2dz'},
#basis={'Cu':'crenbs', 'H':'sto3g'},
#ecp = {'cu':'crenbs'},
charge=-1,
verbose=4)
mf = scf.RHF(mol)
print(mf.kernel(), -196.09477546034623)
mol = gto.M(atom='''
Na 0. 0. 0.
H 0. 0. 1.
''',
basis={'Na':'lanl2dz', 'H':'sto3g'},
ecp = {'Na':'lanl2dz'},
verbose=0)
mf = scf.RHF(mol)
print(mf.kernel(), -0.45002315562861461)
