Skip to content

Solid Harmonics module

contracted_norm(a, w, l)

Compute normalization factor of contracted basis function

Source code in braketlab/solid_harmonics.py 
169
170
171
172
173
174
def contracted_norm(a, w, l):
    """
    Compute normalization factor 
    of contracted basis function
    """ 
    return np.sum(w[:,None]*w[None,:]*(np.sqrt(4.0*a[:,None]*a[None,:])/(a[:,None]+a[None,:]))**(1.5 + l))

f(m)

factorial m!

Source code in braketlab/solid_harmonics.py 
144
145
146
147
148
def f(m):
    """
    factorial m!
    """
    return np.float(np.max([np.prod(np.arange(m)+1), 1]))

get_Nao(a, l, m)

return normalized AO in sympy-format a = exponent l = angular quantum number m = magnetic quantum number

Source code in braketlab/solid_harmonics.py 
122
123
124
125
126
127
128
129
130
131
def get_Nao(a,l,m):
    """
    return normalized AO in sympy-format
    a = exponent
    l = angular quantum number
    m = magnetic quantum number
    """

    a = np.float(a)
    return get_ao(a,l,m)*get_Npi(a,l)*norm_extra(l)

get_Nao_at(pos, a, l, m)

return normalized AO in sympy-format a = exponent l = angular quantum number m = magnetic quantum number

Source code in braketlab/solid_harmonics.py 
133
134
135
136
137
138
139
140
141
def get_Nao_at(pos, a,l,m):
    """
    return normalized AO in sympy-format
    a = exponent
    l = angular quantum number
    m = magnetic quantum number
    """
    a = np.float(a)
    return get_ao_at(pos, a,l,m)*get_Npi(a,l)*norm_extra(l)

get_Nao_lambda(a, l, m)

return a normalized solid harmonic gaussian in numpy lambda format, for convenient evaluation.

Note that every function is centered in (0,0,0) translations should be performed retrospectively

Source code in braketlab/solid_harmonics.py 
157
158
159
160
161
162
163
164
165
166
def get_Nao_lambda(a,l,m):
    """
    return a normalized solid harmonic gaussian
    in numpy lambda format, for convenient evaluation.

    Note that every function is centered in (0,0,0)
    translations should be performed retrospectively
    """
    x,y,z = sp.symbols("x y z")
    return sp.lambdify([x,y,z], get_Nao(a,l,m), "numpy")

get_Npi(a_i, l)

Returns the normalization prefactor for S_lm(a_i, r) a_i = exponent l = angular quantum number

Source code in braketlab/solid_harmonics.py 
114
115
116
117
118
119
120
def get_Npi(a_i, l):
    """
    Returns the normalization prefactor for S_lm(a_i, r)
    a_i = exponent
    l = angular quantum number
    """
    return (2.0*np.pi)**(-.75) * (4.0*a_i)**(0.75 + l/2.0)

get_Slm(l, m)

return the sympy real solid harmonic gaussian S_{lm}(r) as presented in table 6.3 of Helgaker, Jørgensen and Olsen (page 211)

Source code in braketlab/solid_harmonics.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
def get_Slm(l, m):
    """
    return the sympy real solid harmonic gaussian S_{lm}(r) 
    as presented in table 6.3 of Helgaker, Jørgensen and Olsen
    (page 211)
    """
    x,y,z = sp.symbols("x y z")
    r = sp.sqrt(x**2.0 + y**2.0 + z**2.0)

    assert(l<5), "Only l<=4 permitted"
    assert(l>=0), "Invalid l value"
    assert(np.abs(m)<=l), "Invalid m value"

    if l==0:
        if m== 0:
            return x**0

    if l==1:
        if m==1:
            return x
        if m==0:
            return z
        if m==-1:
            return y
    if l==2:
        if m==2:
            return (x**2.0 - y**2.0)*sp.sqrt(3.0)/2.0
        if m==1:
            return x*z*sp.sqrt(3.0)
        if m==0:
            return (3.0*z**2.0 - r**2.0)/2.0
        if m==-1:
            return y*z*sp.sqrt(3.0)
        if m==-2:
            return x*y*sp.sqrt(3.0)

    if l==3:
        if m==3:
            return x*(x**2.0 - 3.0*y**2.0)*sp.sqrt(5.0/2.0)/2.0
        if m==2:
            return z*(x**2.0 - y**2.0)*sp.sqrt(15.0)/2
        if m==1:
            return x*(5.0*z**2.0 - r**2.0)*sp.sqrt(3.0/2.0)/2.0
        if m==0:
            return z*(2.0*z**2.0 - 3.0*x**2.0 - 3.0*y**2.0)/2.0
        if m==-1:
            return y*(5.0*z**2.0 - r**2.0)*sp.sqrt(3.0/2.0)/2.0
        if m==-2:
            return x*y*z*sp.sqrt(15.0) 
        if m==-3:
            return .5*sp.sqrt(5.0/2.0)*(3*x**2 - y**2)*y


    if l==4:
        if m==4:
            return (x**4.0 - 6.0*x**2.0*y**2.0 + y**4.0)*sp.sqrt(35.0)/8.0
        if m==3:
            return (x**2.0 - 3.0*y**2.0)*x*z*sp.sqrt(35.0/2.0)/2.0
        if m==2:
            return (7.0*z**2.0 - r**2.0)*(x**2.0 - y**2.0)*sp.sqrt(5.0)/4.0
        if m==1:
            return (7.0*z**2.0 - 3.0*r**2.0)*x*z*sp.sqrt(5.0/2.0)/2.0
        if m==0:
            return (35.0*z**4.0 - 30.0*z**2.0*r**2.0 + 3.0*r**4.0)/8.0
        if m==-1:
            return (7.0*z**2.0 - 3.0*r**2.0)*y*z*sp.sqrt(5.0/2.0)/2.0
        if m==-2:
            return (7.0*z**2.0 - r**2.0)*x*y*sp.sqrt(5.0)/2.0
        if m==-3:
            return (3.0*x**2.0- y**2.0)*y*z*sp.sqrt(35.0/2.0)/2.0
        if m==-4:
            return (x**2.0 - y**2.0)*x*y*sp.sqrt(35.0)/2.0

get_ao(a, l, m)

return unnormalized solid harmonic gaussian for quantum numbers l, m a = exponent

Source code in braketlab/solid_harmonics.py 
84
85
86
87
88
89
90
91
92
93
94
def get_ao(a, l, m):
    """
    return unnormalized 
    solid harmonic gaussian
    for quantum numbers l, m
    a = exponent
    """
    a = np.float(a)
    x,y,z = sp.symbols("x y z")
    slm = get_Slm(l,m)
    return slm*sp.exp(-sp.UnevaluatedExpr(a)*(x**2.0 + y**2.0 + z**2.0))

get_ao_at(pos, a, l, m)

return unnormalized solid harmonic gaussian for quantum numbers l, m a = exponent

Source code in braketlab/solid_harmonics.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
def get_ao_at(pos, a, l, m):
    """
    return unnormalized 
    solid harmonic gaussian
    for quantum numbers l, m
    a = exponent
    """
    a = np.float(a)
    x,y,z = sp.symbols("x y z")
    slm = get_Slm(l,m)

    chi = slm*sp.exp(-sp.UnevaluatedExpr(a)*(x**2.0 + y**2.0 + z**2.0))
    chi = chi.subs(x, x-pos[0])
    chi = chi.subs(y, y-pos[1])
    chi = chi.subs(z, z-pos[2])

    return chi

get_contracted(a, w, l, m, representation='numeric')

Generates Solid Harmonic Gaussian lambda functions a = exponent

Source code in braketlab/solid_harmonics.py 
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
def get_contracted(a,w,l,m, representation = "numeric"):
    """
    Generates Solid Harmonic Gaussian lambda functions
    a = exponent
    """
    a = np.float(a)
    S = contracted_norm(a,w,l)
    CGO = 0
    for i in np.arange(a.shape[0]):
        CGO += w[i]*get_Nao(a[i],l,m)/np.sqrt(S)

    if representation == "numeric":
        x,y,z = sp.symbols("x y z")

        return sp.lambdify([x,y,z], CGO, "numpy")
    if representation == "sympy":
        return CGO

get_contracted_at(pos, a, w, l, m)

Generates Solid Harmonic Gaussian lambda functions a = exponents

Source code in braketlab/solid_harmonics.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
def get_contracted_at(pos, a,w,l,m):
    """
    Generates Solid Harmonic Gaussian lambda functions
    a = exponents

    """
    S = contracted_norm(a,w,l)
    CGO = 0
    for i in np.arange(a.shape[0]):
        CGO += w[i]*get_Nao_at(pos, a[i],l,m)/np.sqrt(S)
    #print(CGO)
    x,y,z = sp.symbols("x y z")


    return CGO #sp.lambdify([x,y,z], c*CGO, "numpy")

get_contracted_sympy(a, w, l, m)

Generates Solid Harmonic Gaussian lambda functions a = exponents

Source code in braketlab/solid_harmonics.py 
213
214
215
216
217
218
219
220
221
222
223
224
225
def get_contracted_sympy(a,w,l,m):
    """
    Generates Solid Harmonic Gaussian lambda functions
    a = exponents

    """
    a = np.float(a)
    S = contracted_norm(a,w,l)
    CGO = 0
    for i in np.arange(a.shape[0]):
        CGO += w[i]*get_Nao(a[i],l,m)/np.sqrt(S)

    return CGO

norm_extra(l)

Factor required that is not accounted for in eq. 3.3 in LSDalton manual

Source code in braketlab/solid_harmonics.py 
150
151
152
153
154
155
def norm_extra(l):
    """
    Factor required that is _not_ accounted for
    in eq. 3.3 in LSDalton manual
    """
    return (np.array([1.0,1.0,3.0,15.0,105.0])**-.5)[l]