1#!/usr/bin/env python 2# Copyright 2014-2018,2021 The PySCF Developers. All Rights Reserved. 3# 4# Licensed under the Apache License, Version 2.0 (the "License"); 5# you may not use this file except in compliance with the License. 6# You may obtain a copy of the License at 7# 8# http://www.apache.org/licenses/LICENSE-2.0 9# 10# Unless required by applicable law or agreed to in writing, software 11# distributed under the License is distributed on an "AS IS" BASIS, 12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13# See the License for the specific language governing permissions and 14# limitations under the License. 15# 16# Author: Qiming Sun <osirpt.sun@gmail.com> 17# 18 19'''Analytic PP integrals. See also pyscf/pbc/gto/pesudo/pp.py 20 21For GTH/HGH PPs, see: 22 Goedecker, Teter, Hutter, PRB 54, 1703 (1996) 23 Hartwigsen, Goedecker, and Hutter, PRB 58, 3641 (1998) 24''' 25 26import ctypes 27import copy 28import numpy 29import scipy.special 30from pyscf import lib 31from pyscf import gto 32 33libpbc = lib.load_library('libpbc') 34 35def get_pp_loc_part1(cell, kpts=None): 36 '''PRB, 58, 3641 Eq (1), integrals associated to erf 37 ''' 38 raise NotImplementedError 39 40def get_gth_vlocG_part1(cell, Gv): 41 '''PRB, 58, 3641 Eq (5) first term 42 ''' 43 from pyscf.pbc import tools 44 coulG = tools.get_coulG(cell, Gv=Gv) 45 G2 = numpy.einsum('ix,ix->i', Gv, Gv) 46 G0idx = numpy.where(G2==0)[0] 47 48 if cell.dimension != 2 or cell.low_dim_ft_type == 'inf_vacuum': 49 vlocG = numpy.zeros((cell.natm, len(G2))) 50 for ia in range(cell.natm): 51 Zia = cell.atom_charge(ia) 52 symb = cell.atom_symbol(ia) 53 # Note the signs -- potential here is positive 54 vlocG[ia] = Zia * coulG 55 if symb in cell._pseudo: 56 pp = cell._pseudo[symb] 57 rloc, nexp, cexp = pp[1:3+1] 58 vlocG[ia] *= numpy.exp(-0.5*rloc**2 * G2) 59 # alpha parameters from the non-divergent Hartree+Vloc G=0 term. 60 vlocG[ia,G0idx] = -2*numpy.pi*Zia*rloc**2 61 62 elif cell.dimension == 2: 63 # The following 2D ewald summation is taken from: 64 # Minary, Tuckerman, Pihakari, Martyna J. Chem. Phys. 116, 5351 (2002) 65 vlocG = numpy.zeros((cell.natm,len(G2))) 66 b = cell.reciprocal_vectors() 67 inv_area = numpy.linalg.norm(numpy.cross(b[0], b[1]))/(2*numpy.pi)**2 68 lzd2 = cell.vol * inv_area / 2 69 lz = lzd2*2. 70 71 G2[G0idx] = 1e200 72 Gxy = numpy.linalg.norm(Gv[:,:2],axis=1) 73 Gz = abs(Gv[:,2]) 74 75 for ia in range(cell.natm): 76 Zia = cell.atom_charge(ia) 77 symb = cell.atom_symbol(ia) 78 if symb not in cell._pseudo: 79 vlocG[ia] = Zia * coulG 80 continue 81 82 pp = cell._pseudo[symb] 83 rloc, nexp, cexp = pp[1:3+1] 84 85 ew_eta = 1./numpy.sqrt(2)/rloc 86 JexpG2 = 4*numpy.pi / G2 * numpy.exp(-G2/(4*ew_eta**2)) 87 fac = 4*numpy.pi / G2 * numpy.cos(Gz*lzd2) 88 JexpG2 -= fac * numpy.exp(-Gxy*lzd2) 89 eta_z1 = (ew_eta**2 * lz + Gxy) / (2.*ew_eta) 90 eta_z2 = (ew_eta**2 * lz - Gxy) / (2.*ew_eta) 91 JexpG2 += fac * 0.5*(numpy.exp(-eta_z1**2)*scipy.special.erfcx(eta_z2) + 92 numpy.exp(-eta_z2**2)*scipy.special.erfcx(eta_z1)) 93 vlocG[ia,:] = Zia * JexpG2 94 95 JexpG0 = ( - numpy.pi * lz**2 / 2. * scipy.special.erf( ew_eta * lzd2 ) 96 + numpy.pi/ew_eta**2 * scipy.special.erfc(ew_eta*lzd2) 97 - numpy.sqrt(numpy.pi)*lz/ew_eta * numpy.exp( - (ew_eta*lzd2)**2 ) ) 98 vlocG[ia,G0idx] = -2*numpy.pi*Zia*rloc**2 + Zia*JexpG0 99 else: 100 raise NotImplementedError('Low dimension ft_type %s' 101 ' not implemented for dimension %d' % 102 (cell.low_dim_ft_type, cell.dimension)) 103 return vlocG 104 105# part2 Vnuc - Vloc 106def get_pp_loc_part2(cell, kpts=None): 107 '''PRB, 58, 3641 Eq (1), integrals associated to C1, C2, C3, C4 108 ''' 109 from pyscf.pbc.df import incore 110 if kpts is None: 111 kpts_lst = numpy.zeros((1,3)) 112 else: 113 kpts_lst = numpy.reshape(kpts, (-1,3)) 114 nkpts = len(kpts_lst) 115 116 intors = ('int3c2e', 'int3c1e', 'int3c1e_r2_origk', 117 'int3c1e_r4_origk', 'int3c1e_r6_origk') 118 kptij_lst = numpy.hstack((kpts_lst,kpts_lst)).reshape(-1,2,3) 119 buf = 0 120 for cn in range(1, 5): 121 fakecell = fake_cell_vloc(cell, cn) 122 if fakecell.nbas > 0: 123 v = incore.aux_e2(cell, fakecell, intors[cn], aosym='s2', comp=1, 124 kptij_lst=kptij_lst) 125 buf += numpy.einsum('...i->...', v) 126 127 if isinstance(buf, int): 128 if any(cell.atom_symbol(ia) in cell._pseudo for ia in range(cell.natm)): 129 pass 130 else: 131 lib.logger.warn(cell, 'cell.pseudo was specified but its elements %s ' 132 'were not found in the system.', cell._pseudo.keys()) 133 vpploc = [0] * nkpts 134 else: 135 buf = buf.reshape(nkpts,-1) 136 vpploc = [] 137 for k, kpt in enumerate(kpts_lst): 138 v = lib.unpack_tril(buf[k]) 139 if abs(kpt).sum() < 1e-9: # gamma_point: 140 v = v.real 141 vpploc.append(v) 142 if kpts is None or numpy.shape(kpts) == (3,): 143 vpploc = vpploc[0] 144 return vpploc 145 146 147def get_pp_nl(cell, kpts=None): 148 if kpts is None: 149 kpts_lst = numpy.zeros((1,3)) 150 else: 151 kpts_lst = numpy.reshape(kpts, (-1,3)) 152 nkpts = len(kpts_lst) 153 154 fakecell, hl_blocks = fake_cell_vnl(cell) 155 ppnl_half = _int_vnl(cell, fakecell, hl_blocks, kpts_lst) 156 nao = cell.nao_nr() 157 buf = numpy.empty((3*9*nao), dtype=numpy.complex128) 158 159 # We set this equal to zeros in case hl_blocks loop is skipped 160 # and ppnl is returned 161 ppnl = numpy.zeros((nkpts,nao,nao), dtype=numpy.complex128) 162 for k, kpt in enumerate(kpts_lst): 163 offset = [0] * 3 164 for ib, hl in enumerate(hl_blocks): 165 l = fakecell.bas_angular(ib) 166 nd = 2 * l + 1 167 hl_dim = hl.shape[0] 168 ilp = numpy.ndarray((hl_dim,nd,nao), dtype=numpy.complex128, buffer=buf) 169 for i in range(hl_dim): 170 p0 = offset[i] 171 ilp[i] = ppnl_half[i][k][p0:p0+nd] 172 offset[i] = p0 + nd 173 ppnl[k] += numpy.einsum('ilp,ij,jlq->pq', ilp.conj(), hl, ilp) 174 175 if abs(kpts_lst).sum() < 1e-9: # gamma_point: 176 ppnl = ppnl.real 177 178 if kpts is None or numpy.shape(kpts) == (3,): 179 ppnl = ppnl[0] 180 return ppnl 181 182 183def fake_cell_vloc(cell, cn=0): 184 '''Generate fake cell for V_{loc}. 185 186 Each term of V_{loc} (erf, C_1, C_2, C_3, C_4) is a gaussian type 187 function. The integral over V_{loc} can be transfered to the 3-center 188 integrals, in which the auxiliary basis is given by the fake cell. 189 190 The kwarg cn indiciates which term to generate for the fake cell. 191 If cn = 0, the erf term is generated. C_1,..,C_4 are generated with cn = 1..4 192 ''' 193 fake_env = [cell.atom_coords().ravel()] 194 fake_atm = cell._atm.copy() 195 fake_atm[:,gto.PTR_COORD] = numpy.arange(0, cell.natm*3, 3) 196 ptr = cell.natm * 3 197 fake_bas = [] 198 half_sph_norm = .5/numpy.sqrt(numpy.pi) 199 for ia in range(cell.natm): 200 if cell.atom_charge(ia) == 0: # pass ghost atoms 201 continue 202 203 symb = cell.atom_symbol(ia) 204 if cn == 0: 205 if symb in cell._pseudo: 206 pp = cell._pseudo[symb] 207 rloc, nexp, cexp = pp[1:3+1] 208 alpha = .5 / rloc**2 209 else: 210 alpha = 1e16 211 norm = half_sph_norm / gto.gaussian_int(2, alpha) 212 fake_env.append([alpha, norm]) 213 fake_bas.append([ia, 0, 1, 1, 0, ptr, ptr+1, 0]) 214 ptr += 2 215 elif symb in cell._pseudo: 216 pp = cell._pseudo[symb] 217 rloc, nexp, cexp = pp[1:3+1] 218 if cn <= nexp: 219 alpha = .5 / rloc**2 220 norm = cexp[cn-1]/rloc**(cn*2-2) / half_sph_norm 221 fake_env.append([alpha, norm]) 222 fake_bas.append([ia, 0, 1, 1, 0, ptr, ptr+1, 0]) 223 ptr += 2 224 225 fakecell = copy.copy(cell) 226 fakecell._atm = numpy.asarray(fake_atm, dtype=numpy.int32) 227 fakecell._bas = numpy.asarray(fake_bas, dtype=numpy.int32) 228 fakecell._env = numpy.asarray(numpy.hstack(fake_env), dtype=numpy.double) 229 return fakecell 230 231# sqrt(Gamma(l+1.5)/Gamma(l+2i+1.5)) 232_PLI_FAC = 1/numpy.sqrt(numpy.array(( 233 (1, 3.75 , 59.0625 ), # l = 0, 234 (1, 8.75 , 216.5625 ), # l = 1, 235 (1, 15.75, 563.0625 ), # l = 2, 236 (1, 24.75, 1206.5625), # l = 3, 237 (1, 35.75, 2279.0625), # l = 4, 238 (1, 48.75, 3936.5625), # l = 5, 239 (1, 63.75, 6359.0625), # l = 6, 240 (1, 80.75, 9750.5625))))# l = 7, 241 242def fake_cell_vnl(cell): 243 '''Generate fake cell for V_{nl}. 244 245 gaussian function p_i^l Y_{lm} 246 ''' 247 fake_env = [cell.atom_coords().ravel()] 248 fake_atm = cell._atm.copy() 249 fake_atm[:,gto.PTR_COORD] = numpy.arange(0, cell.natm*3, 3) 250 ptr = cell.natm * 3 251 fake_bas = [] 252 hl_blocks = [] 253 for ia in range(cell.natm): 254 if cell.atom_charge(ia) == 0: # pass ghost atoms 255 continue 256 257 symb = cell.atom_symbol(ia) 258 if symb in cell._pseudo: 259 pp = cell._pseudo[symb] 260 # nproj_types = pp[4] 261 for l, (rl, nl, hl) in enumerate(pp[5:]): 262 if nl > 0: 263 alpha = .5 / rl**2 264 norm = gto.gto_norm(l, alpha) 265 fake_env.append([alpha, norm]) 266 fake_bas.append([ia, l, 1, 1, 0, ptr, ptr+1, 0]) 267 268# 269# Function p_i^l (PRB, 58, 3641 Eq 3) is (r^{2(i-1)})^2 square normalized to 1. 270# But here the fake basis is square normalized to 1. A factor ~ p_i^l / p_1^l 271# is attached to h^l_ij (for i>1,j>1) so that (factor * fake-basis * r^{2(i-1)}) 272# is normalized to 1. The factor is 273# r_l^{l+(4-1)/2} sqrt(Gamma(l+(4-1)/2)) sqrt(Gamma(l+3/2)) 274# ------------------------------------------ = ---------------------------------- 275# r_l^{l+(4i-1)/2} sqrt(Gamma(l+(4i-1)/2)) sqrt(Gamma(l+2i-1/2)) r_l^{2i-2} 276# 277 fac = numpy.array([_PLI_FAC[l,i]/rl**(i*2) for i in range(nl)]) 278 hl = numpy.einsum('i,ij,j->ij', fac, numpy.asarray(hl), fac) 279 hl_blocks.append(hl) 280 ptr += 2 281 282 fakecell = copy.copy(cell) 283 fakecell._atm = numpy.asarray(fake_atm, dtype=numpy.int32) 284 fakecell._bas = numpy.asarray(fake_bas, dtype=numpy.int32) 285 fakecell._env = numpy.asarray(numpy.hstack(fake_env), dtype=numpy.double) 286 return fakecell, hl_blocks 287 288def _int_vnl(cell, fakecell, hl_blocks, kpts): 289 '''Vnuc - Vloc''' 290 rcut = max(cell.rcut, fakecell.rcut) 291 Ls = cell.get_lattice_Ls(rcut=rcut) 292 nimgs = len(Ls) 293 expkL = numpy.asarray(numpy.exp(1j*numpy.dot(kpts, Ls.T)), order='C') 294 nkpts = len(kpts) 295 296 fill = getattr(libpbc, 'PBCnr2c_fill_ks1') 297 intopt = lib.c_null_ptr() 298 299 def int_ket(_bas, intor): 300 if len(_bas) == 0: 301 return [] 302 intor = cell._add_suffix(intor) 303 atm, bas, env = gto.conc_env(cell._atm, cell._bas, cell._env, 304 fakecell._atm, _bas, fakecell._env) 305 atm = numpy.asarray(atm, dtype=numpy.int32) 306 bas = numpy.asarray(bas, dtype=numpy.int32) 307 env = numpy.asarray(env, dtype=numpy.double) 308 natm = len(atm) 309 nbas = len(bas) 310 shls_slice = (cell.nbas, nbas, 0, cell.nbas) 311 ao_loc = gto.moleintor.make_loc(bas, intor) 312 ni = ao_loc[shls_slice[1]] - ao_loc[shls_slice[0]] 313 nj = ao_loc[shls_slice[3]] - ao_loc[shls_slice[2]] 314 out = numpy.empty((nkpts,ni,nj), dtype=numpy.complex128) 315 comp = 1 316 317 fintor = getattr(gto.moleintor.libcgto, intor) 318 319 drv = libpbc.PBCnr2c_drv 320 drv(fintor, fill, out.ctypes.data_as(ctypes.c_void_p), 321 ctypes.c_int(nkpts), ctypes.c_int(comp), ctypes.c_int(nimgs), 322 Ls.ctypes.data_as(ctypes.c_void_p), 323 expkL.ctypes.data_as(ctypes.c_void_p), 324 (ctypes.c_int*4)(*(shls_slice[:4])), 325 ao_loc.ctypes.data_as(ctypes.c_void_p), intopt, lib.c_null_ptr(), 326 atm.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(natm), 327 bas.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(nbas), 328 env.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(env.size)) 329 return out 330 331 hl_dims = numpy.asarray([len(hl) for hl in hl_blocks]) 332 out = (int_ket(fakecell._bas[hl_dims>0], 'int1e_ovlp'), 333 int_ket(fakecell._bas[hl_dims>1], 'int1e_r2_origi'), 334 int_ket(fakecell._bas[hl_dims>2], 'int1e_r4_origi')) 335 return out 336 337