1#!/usr/bin/env python 2# Copyright 2014-2018 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'''radii grids''' 20 21import numpy 22from pyscf.data import radii 23from pyscf.data.elements import charge as elements_proton 24 25BRAGG_RADII = radii.BRAGG 26COVALENT_RADII = radii.COVALENT 27 28# P.M.W. Gill, B.G. Johnson, J.A. Pople, Chem. Phys. Letters 209 (1993) 506-512 29SG1RADII = numpy.array(( 30 0, 31 1.0000, 0.5882, 32 3.0769, 2.0513, 1.5385, 1.2308, 1.0256, 0.8791, 0.7692, 0.6838, 33 4.0909, 3.1579, 2.5714, 2.1687, 1.8750, 1.6514, 1.4754, 1.3333)) 34 35 36# Murray, N.C. Handy, G.J. Laming, Mol. Phys. 78, 997(1993) 37def murray(n, *args, **kwargs): 38 raise RuntimeError('Not implemented') 39 40# Gauss-Chebyshev of the first kind, and the transformed interval [0,\infty) 41def becke(n, charge, *args, **kwargs): 42 '''Becke, JCP 88, 2547 (1988); DOI:10.1063/1.454033''' 43 if charge == 1: 44 rm = BRAGG_RADII[charge] 45 else: 46 rm = BRAGG_RADII[charge] * .5 47 t, w = numpy.polynomial.chebyshev.chebgauss(n) 48 r = (1+t)/(1-t) * rm 49 w *= 2/(1-t)**2 * rm 50 return r, w 51 52# scale rad and rad_weight if necessary 53# gauss-legendre 54def delley(n, *args, **kwargs): 55 '''B. Delley radial grids. Ref. JCP 104, 9848 (1996); DOI:10.1063/1.471749. log2 algorithm''' 56 r = numpy.empty(n) 57 dr = numpy.empty(n) 58 r_outer = 12. 59 step = 1. / (n+1) 60 rfac = r_outer / numpy.log(1 - (n*step)**2) 61 for i in range(1, n+1): 62 xi = rfac * numpy.log(1-(i*step)**2) 63 r[i-1] = xi 64 dri = rfac * (-2.0*i*(step)**2) / ((1-(i*step)**2)) # d xi / dr 65 dr[i-1] = dri 66 return r, dr 67gauss_legendre = delley 68 69def mura_knowles(n, charge=None, *args, **kwargs): 70 '''Mura-Knowles [JCP 104, 9848 (1996); DOI:10.1063/1.471749] log3 quadrature radial grids''' 71 r = numpy.empty(n) 72 dr = numpy.empty(n) 73# 7 for Li, Be, Na, Mg, K, Ca, otherwise 5 74 if charge in (3, 4, 11, 12, 19, 20): 75 far = 7 76 else: 77 far = 5.2 78 for i in range(n): 79 x = (i+.5) / n 80 r[i] = -far * numpy.log(1-x**3) 81 dr[i] = far * 3*x*x/((1-x**3)*n) 82 return r, dr 83 84# Gauss-Chebyshev of the second kind, and the transformed interval [0,\infty) 85# Ref Matthias Krack and Andreas M. Koster, J. Chem. Phys. 108 (1998), 3226 86def gauss_chebyshev(n, *args, **kwargs): 87 '''Gauss-Chebyshev [JCP 108, 3226 (1998); DOI:10.1063/1.475719) radial grids''' 88 ln2 = 1 / numpy.log(2) 89 fac = 16./3 / (n+1) 90 x1 = numpy.arange(1,n+1) * numpy.pi / (n+1) 91 xi = ((n-1-numpy.arange(n)*2) / (n+1.) + 92 (1+2./3*numpy.sin(x1)**2) * numpy.sin(2*x1) / numpy.pi) 93 xi = (xi - xi[::-1])/2 94 r = 1 - numpy.log(1+xi) * ln2 95 dr = fac * numpy.sin(x1)**4 * ln2/(1+xi) 96 return r, dr 97 98 99def treutler_ahlrichs(n, *args, **kwargs): 100 ''' 101 Treutler-Ahlrichs [JCP 102, 346 (1995); DOI:10.1063/1.469408] (M4) radial grids 102 ''' 103 r = numpy.empty(n) 104 dr = numpy.empty(n) 105 step = numpy.pi / (n+1) 106 ln2 = 1 / numpy.log(2) 107 for i in range(n): 108 x = numpy.cos((i+1)*step) 109 r [i] = -ln2*(1+x)**.6 * numpy.log((1-x)/2) 110 dr[i] = step * numpy.sin((i+1)*step) \ 111 * ln2*(1+x)**.6 *(-.6/(1+x)*numpy.log((1-x)/2)+1/(1-x)) 112 return r[::-1], dr[::-1] 113treutler = treutler_ahlrichs 114 115 116 117 118def becke_atomic_radii_adjust(mol, atomic_radii): 119 '''Becke atomic radii adjust function''' 120# Becke atomic size adjustment. J. Chem. Phys. 88, 2547 121# i > j 122# fac(i,j) = \frac{1}{4} ( \frac{ra(j)}{ra(i)} - \frac{ra(i)}{ra(j)} 123# fac(j,i) = -fac(i,j) 124 charges = [elements_proton(x) for x in mol.elements] 125 rad = atomic_radii[charges] + 1e-200 126 rr = rad.reshape(-1,1) * (1./rad) 127 a = .25 * (rr.T - rr) 128 a[a<-.5] = -.5 129 a[a>0.5] = 0.5 130 #:return lambda i,j,g: g + a[i,j]*(1-g**2) 131 def fadjust(i, j, g): 132 g1 = g**2 133 g1 -= 1. 134 g1 *= -a[i,j] 135 g1 += g 136 return g1 137 return fadjust 138 139def treutler_atomic_radii_adjust(mol, atomic_radii): 140 '''Treutler atomic radii adjust function: [JCP 102, 346 (1995); DOI:10.1063/1.469408]''' 141# JCP 102, 346 (1995) 142# i > j 143# fac(i,j) = \frac{1}{4} ( \frac{ra(j)}{ra(i)} - \frac{ra(i)}{ra(j)} 144# fac(j,i) = -fac(i,j) 145 charges = [elements_proton(x) for x in mol.elements] 146 rad = numpy.sqrt(atomic_radii[charges]) + 1e-200 147 rr = rad.reshape(-1,1) * (1./rad) 148 a = .25 * (rr.T - rr) 149 a[a<-.5] = -.5 150 a[a>0.5] = 0.5 151 #:return lambda i,j,g: g + a[i,j]*(1-g**2) 152 def fadjust(i, j, g): 153 g1 = g**2 154 g1 -= 1. 155 g1 *= -a[i,j] 156 g1 += g 157 return g1 158 return fadjust 159 160