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