<jupyter_start><jupyter_text>AO to SO transformation<jupyter_code>import psi4
import numpy as np
import sys
sys.path.append('../lib')
from P4toC4_aux import *
BASIS='def2-SV'
#BASIS='cc-pvdz'
psi4.set_memory('500 MB')
psi4.core.set_global_option("BASIS", BASIS)
psi4.core.set_global_option("SCF_TYPE", "pk")
psi4.core.set_global_option("REFERENCE", "RHF")
psi4.core.set_global_option("D_CONVERGENCE", 1e-8)
psi4.core.set_global_option("PUREAM", "True")
psi4.core.set_output_file('output.dat', False)

# O
# H 1 0.96
# H 1 0.96 2 104.5

h2o = psi4.geometry("""
0 1
O           -0.005580556816     0.124770196240     0.000000000000
H           -1.409609476816    -1.024054213760     0.000000000000
H            1.498176963184    -0.956139343760     0.000000000000
units Bohr
no_reorient
symmetry cs
""")

E, wf = psi4.energy('scf', return_wfn=True)
E
mol=wf.molecule()
ptgr=mol.point_group()
print(ptgr.symbol(), ptgr.order())
#psi4.core.Vector.array_interface(wf.epsilon_a())
n_irrep=wf.nirrep()
g=wf.nmopi()
n_mo_pi=g.to_tuple()
n_mo_pi
# Psi4 MOs in SO basis
Ca=wf.Ca()
Ca.shape
C_AO=wf.Ca_subset('AO','ALL')
print(C_AO.shape)
C_AO_by_E=psi4.core.Matrix.array_interface(C_AO)[0]
#C_AO[:,0:6].round(3)<jupyter_output>(13, 13)<jupyter_text>Reorder C_AO by symmetry* establish j_by_e by argsort of symmetry ordered epsilons* reverse j_by_e mapping to j_by_sym* use j_by_sym to reorder C_AO_by_E* make a Psi4.Matrix from the symmetry blocks of C_AO_by_sym  * C_AO_by_sym.nph[irrep] has two indices  * first index is AO  * second index is MO of irrep<jupyter_code>eps_pi=wf.epsilon_a()
#print(n_irrep, n_mo_pi)
eps_by_sym=np.concatenate(eps_pi.nph)
#print(eps_by_sym)
j_by_e=np.argsort(eps_by_sym)
print('j_by_energy:', j_by_e)
#print('Reordered orbital energies')
#print(all_eps[j_by_e])
n = len(j_by_e)
j_by_sym=np.zeros(n, int)
for i in range(n):
    j_by_sym[j_by_e[i]]=i
print('j_by_symmetry:', j_by_sym)
#
#  Test by sorting epsilsons by energy and applying j_by_sym vs all
#
eps_by_e=np.sort(eps_by_sym)
np.max(eps_by_e[j_by_sym]-eps_by_sym)
C_AO_by_sym=C_AO_by_E[:,j_by_sym]
n_irrep=wf.nirrep()
n_mo_pi=wf.nmopi()
N=n_mo_pi.to_tuple()
irrep_lst = []
offset = 0
for i in range(n_irrep):
    dim = N[i]
    irrep_lst.append(C_AO_by_sym[:,offset:offset+dim])
    offset += dim
C_AO_pi=psi4.core.Matrix.from_array(irrep_lst)
print(C_AO_pi.shape)
#C_AO_pi.nph[1]<jupyter_output>((13, 11), (13, 2))<jupyter_text>AO to SO transformation* I want C_SO[irrep] = L[irrep] * C_AO[irrep]* The left matrix L has an all-AO and an SO[irrep] index.* The C_AO matrix are has an all-AO and an MO[irrep] index.* Options:  * `wf.sobasisset()`: two matrices, which are transposes of each other.  * `wf.aotoso()`: one matrix, looks like a left matrix.<jupyter_code># no clue what a petite_list contains
# seems clear for STO-3G, but eludes me for def2-SV
# I=wf.sobasisset()
# pl=I.petite_list()
# Rs=pl.aotoso()
# Ls=pl.sotoao()
# print(Ls.shape)
# Ls.nph[1]
#
#  this looks right for def2-SV
#
Ls=wf.aotoso()
print(Ls.shape)
#Ls.nph[1]
for sym in range(n_irrep):
    L=Ls.nph[sym]
    C=C_AO_pi.nph[sym]
    C_SO=np.matmul(np.transpose(L),C)
    #np.round(C_SO,4)
    print(f'{sym}   delta={np.max(C_SO-Ca.nph[sym]):10.3e}')<jupyter_output>0   delta= 0.000e+00
1   delta= 0.000e+00
{'fraction_non_alphanumeric': 0.1485906862745098, 'fraction_numerical': 0.06004901960784314, 'mean_word_length': 5.183712121212121, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '122186', 'license': 'MIT', 'n_tokens_mistral': 1654, 'n_tokens_neox': 1492, 'n_words': 280, 'path': 'dev/.ipynb_checkpoints/AO_TO_SO-checkpoint.ipynb', 'repo_name': 'tsommerfeld/Psi4_to_OLDMOS', 'size': '3264'}