Source code for pyscf.symm.sph

#!/usr/bin/env python
# Copyright 2014-2019 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.

'''
Spherical harmonics
'''

import numpy
import scipy.linalg
from pyscf.symm.cg import cg_spin

[docs] def real_sph_vec(r, lmax, reorder_p=False): '''Computes (all) real spherical harmonics up to the angular momentum lmax''' #:import scipy.special #:ngrid = r.shape[0] #:cosphi = r[:,2] #:sinphi = (1-cosphi**2)**.5 #:costheta = numpy.ones(ngrid) #:sintheta = numpy.zeros(ngrid) #:costheta[sinphi!=0] = r[sinphi!=0,0] / sinphi[sinphi!=0] #:sintheta[sinphi!=0] = r[sinphi!=0,1] / sinphi[sinphi!=0] #:costheta[costheta> 1] = 1 #:costheta[costheta<-1] =-1 #:sintheta[sintheta> 1] = 1 #:sintheta[sintheta<-1] =-1 #:varphi = numpy.arccos(cosphi) #:theta = numpy.arccos(costheta) #:theta[sintheta<0] = 2*numpy.pi - theta[sintheta<0] #:ylms = [] #:for l in range(lmax+1): #: ylm = numpy.empty((l*2+1,ngrid)) #: ylm[l] = scipy.special.sph_harm(0, l, theta, varphi).real #: for m in range(1, l+1): #: f1 = scipy.special.sph_harm(-m, l, theta, varphi) #: f2 = scipy.special.sph_harm( m, l, theta, varphi) #: # complex to real spherical functions #: if m % 2 == 1: #: ylm[l-m] = (-f1.imag - f2.imag) / numpy.sqrt(2) #: ylm[l+m] = ( f1.real - f2.real) / numpy.sqrt(2) #: else: #: ylm[l-m] = (-f1.imag + f2.imag) / numpy.sqrt(2) #: ylm[l+m] = ( f1.real + f2.real) / numpy.sqrt(2) #: ylms.append(ylm) #:return ylms # When r is a normalized vector: norm = 1./numpy.linalg.norm(r, axis=1) return multipoles(r*norm.reshape(-1,1), lmax, reorder_p)
[docs] def multipoles(r, lmax, reorder_dipole=True): ''' Compute all multipoles upto lmax rad = numpy.linalg.norm(r, axis=1) ylms = real_ylm(r/rad.reshape(-1,1), lmax) pol = [rad**l*y for l, y in enumerate(ylms)] Kwargs: reorder_p : bool sort dipole to the order (x,y,z) ''' from pyscf import gto # libcint cart2sph transformation provide the capability to compute # multipole directly. cart2sph function is fast for low angular moment. ngrid = r.shape[0] xs = numpy.ones((lmax+1,ngrid)) ys = numpy.ones((lmax+1,ngrid)) zs = numpy.ones((lmax+1,ngrid)) for i in range(1,lmax+1): xs[i] = xs[i-1] * r[:,0] ys[i] = ys[i-1] * r[:,1] zs[i] = zs[i-1] * r[:,2] ylms = [] for l in range(lmax+1): nd = (l+1)*(l+2)//2 c = numpy.empty((nd,ngrid)) k = 0 for lx in reversed(range(0, l+1)): for ly in reversed(range(0, l-lx+1)): lz = l - lx - ly c[k] = xs[lx] * ys[ly] * zs[lz] k += 1 ylm = gto.cart2sph(l, c.T).T ylms.append(ylm) # when call libcint, p functions are ordered as px,py,pz # reorder px,py,pz to p(-1),p(0),p(1) if (not reorder_dipole) and lmax >= 1: ylms[1] = ylms[1][[1,2,0]] return ylms
[docs] def sph_pure2real(l, reorder_p=True): r''' Transformation matrix: from the pure spherical harmonic functions Y_m to the real spherical harmonic functions O_m:: O_m = \sum Y_m' * U(m',m) Y(-1) = 1/\sqrt(2){-iO(-1) + O(1)}; Y(1) = 1/\sqrt(2){-iO(-1) - O(1)} Y(-2) = 1/\sqrt(2){-iO(-2) + O(2)}; Y(2) = 1/\sqrt(2){iO(-2) + O(2)} O(-1) = i/\sqrt(2){Y(-1) + Y(1)}; O(1) = 1/\sqrt(2){Y(-1) - Y(1)} O(-2) = i/\sqrt(2){Y(-2) - Y(2)}; O(2) = 1/\sqrt(2){Y(-2) + Y(2)} Kwargs: reorder_p (bool): Whether the p functions are in the (x,y,z) order. Returns: 2D array U_{complex,real} ''' n = 2 * l + 1 u = numpy.zeros((n,n), dtype=complex) sqrthfr = numpy.sqrt(.5) sqrthfi = numpy.sqrt(.5)*1j if reorder_p and l == 1: u[1,2] = 1 u[0,1] = sqrthfi u[2,1] = sqrthfi u[0,0] = sqrthfr u[2,0] = -sqrthfr else: u[l,l] = 1 for m in range(1, l+1, 2): u[l-m,l-m] = sqrthfi u[l+m,l-m] = sqrthfi u[l-m,l+m] = sqrthfr u[l+m,l+m] = -sqrthfr for m in range(2, l+1, 2): u[l-m,l-m] = sqrthfi u[l+m,l-m] = -sqrthfi u[l-m,l+m] = sqrthfr u[l+m,l+m] = sqrthfr return u
[docs] def sph_real2pure(l, reorder_p=True): ''' Transformation matrix: from real spherical harmonic functions to the pure spherical harmonic functions. Kwargs: reorder_p (bool): Whether the real p functions are in the (x,y,z) order. ''' # numpy.linalg.inv(sph_pure2real(l)) return sph_pure2real(l, reorder_p).conj().T
# |spinor> = (|real_sph>, |real_sph>) * / u_alpha \ # \ u_beta / # Return 2D array U_{sph,spinor}
[docs] def sph2spinor(l, reorder_p=True): if l == 0: return numpy.array((0., 1.)).reshape(1,-1), \ numpy.array((1., 0.)).reshape(1,-1) else: u1 = sph_real2pure(l, reorder_p) ua = numpy.zeros((2*l+1,4*l+2),dtype=complex) ub = numpy.zeros((2*l+1,4*l+2),dtype=complex) j = l * 2 - 1 mla = l + (-j-1)//2 mlb = l + (-j+1)//2 for k,mj in enumerate(range(-j, j+1, 2)): ua[:,k] = u1[:,mla] * cg_spin(l, j, mj, 1) ub[:,k] = u1[:,mlb] * cg_spin(l, j, mj,-1) mla += 1 mlb += 1 j = l * 2 + 1 mla = l + (-j-1)//2 mlb = l + (-j+1)//2 for k,mj in enumerate(range(-j, j+1, 2)): if mla < 0: ua[:,l*2+k] = 0 else: ua[:,l*2+k] = u1[:,mla] * cg_spin(l, j, mj, 1) if mlb >= 2*l+1: ub[:,l*2+k] = 0 else: ub[:,l*2+k] = u1[:,mlb] * cg_spin(l, j, mj,-1) mla += 1 mlb += 1 return ua, ub
real2spinor = sph2spinor # Returns 2D array U_{sph,spinor}
[docs] def sph2spinor_coeff(mol): '''Transformation matrix that transforms real-spherical GTOs to spinor GTOs for all basis functions Examples:: >>> from pyscf import gto >>> from pyscf.symm import sph >>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvtz') >>> ca, cb = sph.sph2spinor_coeff(mol) >>> s0 = mol.intor('int1e_ovlp_spinor') >>> s1 = ca.conj().T.dot(mol.intor('int1e_ovlp_sph')).dot(ca) >>> s1+= cb.conj().T.dot(mol.intor('int1e_ovlp_sph')).dot(cb) >>> print(abs(s1-s0).max()) >>> 6.66133814775e-16 ''' lmax = max([mol.bas_angular(i) for i in range(mol.nbas)]) ualst = [] ublst = [] for l in range(lmax+1): u1, u2 = sph2spinor(l, reorder_p=True) ualst.append(u1) ublst.append(u2) ca = [] cb = [] for ib in range(mol.nbas): l = mol.bas_angular(ib) kappa = mol.bas_kappa(ib) if kappa == 0: ua = ualst[l] ub = ublst[l] elif kappa < 0: ua = ualst[l][:,l*2:] ub = ublst[l][:,l*2:] else: ua = ualst[l][:,:l*2] ub = ublst[l][:,:l*2] nctr = mol.bas_nctr(ib) ca.extend([ua]*nctr) cb.extend([ub]*nctr) return numpy.stack([scipy.linalg.block_diag(*ca), scipy.linalg.block_diag(*cb)])
real2spinor_whole = sph2spinor_coeff
[docs] def cart2spinor(l): '''Cartesian to spinor for angular moment l''' from pyscf import gto return gto.cart2spinor_l(l)
if __name__ == '__main__': for l in range(3): print(sph_pure2real(l)) print(sph_real2pure(l)) for l in range(3): print(sph2spinor(l)[0]) print(sph2spinor(l)[1])