I/O modules

docs/generate_rst.py

@@ -325,7 +325,7 @@ def render_index_rst(project: str, modules: list[ModuleInfo], out_dir: Path) ->
     out.append(".. toctree::\n   :caption: Modules\n   :hidden:\n\n")
 
 
-    module_order = ['spahm', 'fields', 'qml', 'regression'] + ['compound', 'c2mio', 'equio', 'orcaio']
+    module_order = ['spahm', 'fields', 'qml', 'regression', 'compound']
     def key(m):
         name = m.name.split('.')[1]
         i = module_order.index(name) if name in module_order else len(module_order)

examples/data/AJALIH/AJALIH.xyz

@@ -0,0 +1,54 @@
+52
+0 5
+Fe         0.14005       -0.96028        0.24860
+Br         0.04760       -2.71770       -1.35052
+Br         0.51098       -1.66925        2.48994
+P         -1.67045        0.64630        0.04154
+P          1.49137        1.00407       -0.30691
+N         -0.93908        2.01059       -0.66228
+N          0.40553        2.31997       -0.34154
+H          0.45433        2.85659        0.53275
+C         -2.55204        1.05118        1.62253
+H         -2.90953        0.06097        1.96346
+C         -1.54412        1.56608        2.63497
+H         -2.01575        1.68908        3.62514
+H         -1.15822        2.55927        2.33684
+H         -0.69342        0.87423        2.74906
+C         -3.72734        2.01525        1.48799
+H         -4.18847        2.18853        2.47612
+H         -4.51900        1.63988        0.81933
+H         -3.39664        2.99655        1.10653
+C         -3.00694        0.25147       -1.15704
+H         -3.73205        1.08518       -1.13859
+C         -2.44534        0.09384       -2.58396
+H         -3.26643       -0.15943       -3.27617
+H         -1.70660       -0.72376       -2.61226
+H         -1.95853        1.01106       -2.94640
+C         -3.71205       -1.03192       -0.70123
+H         -4.50441       -1.30055       -1.42033
+H         -4.18484       -0.92853        0.28944
+H         -2.99492       -1.86941       -0.65824
+C         -1.63840        3.16764       -1.22300
+H         -1.09847        3.52187       -2.11776
+H         -1.70279        4.00781       -0.50205
+H         -2.66091        2.88608       -1.51156
+C          2.68323        1.67170        0.94188
+H          2.01024        1.86959        1.79906
+C          3.70079        0.62535        1.40845
+H          4.24437        1.00759        2.28929
+H          4.44836        0.40440        0.63077
+H          3.20072       -0.31373        1.69785
+C          3.36116        2.98145        0.51065
+H          3.92004        3.41133        1.35970
+H          2.64130        3.73739        0.15615
+H          4.08532        2.80049       -0.30140
+C          2.31022        1.03631       -1.96046
+H          2.77156        2.03284       -2.07646
+C          3.35556       -0.07419       -2.02537
+H          3.75174       -0.15693       -3.05216
+H          2.90140       -1.04685       -1.76654
+H          4.21170        0.09757       -1.35527
+C          1.26293        0.80806       -3.03506
+H          1.73652        0.82147       -4.03221
+H          0.47748        1.58087       -3.02229
+H          0.77939       -0.17371       -2.89981

examples/data/CEHZOF/CEHZOF.xyz

@@ -0,0 +1,28 @@
+26
+0 1
+  I     -0.691564    6.065851   10.218679
+  Co    -0.896730    6.291743   12.801124
+  S     -3.111547    5.942715   12.672839
+  C     -0.900957    4.763642   14.175690
+  C     -0.461276    5.983590   14.795995
+  C      0.713137    6.405881   14.113956
+  H      1.264840    7.323863   14.311032
+  C      0.993227    5.481683   13.083596
+  H      1.802106    5.559932   12.358853
+  C      0.015965    4.473016   13.122433
+  H     -0.056529    3.635907   12.428292
+  C     -2.210760    4.059666   14.447381
+  H     -2.673227    4.493474   15.348774
+  H     -2.035736    2.991419   14.671175
+  C     -3.111610    4.202498   13.259335
+  H     -2.727972    3.649740   12.386330
+  H     -4.146425    3.867362   13.435280
+  C     -3.945153    6.793732   14.050514
+  H     -3.869637    7.852090   13.755824
+  H     -3.340371    6.686944   14.965041
+  C     -5.368941    6.350425   14.251741
+  H     -5.442055    5.313038   14.619953
+  H     -5.857396    6.994746   15.002527
+  H     -5.951474    6.421791   13.318039
+  I     -1.173397    8.852094   12.837558
+  H     -0.948251    6.502776   15.620742

examples/example_deco.py

@@ -9,10 +9,8 @@
 mol  = compound.xyz_to_mol(path+'/data/H2O.xyz', 'def2svp', charge=0, spin=0)
 dm   = fields.dm.get_converged_dm(mol,xc="pbe")
 auxmol, c = decompose(mol, dm, 'cc-pvqz jkfit')
-print("Expansion Coefficients:", c)
 
 N = moments.r2_c(auxmol, c, moments=[0])[0]
-
 print("Number of electrons after decomposition: ", N)
 
 coeffs_to_cube(auxmol, c, 'H2O.cube')
@@ -25,5 +23,6 @@
 N = moments.r2_c(auxmol, c, moments=[0])[0]
 print(N)
 
-
 c = correct_N_atomic(auxmol, np.array([8,1,1]), c)
+N = moments.r2_c(auxmol, c, moments=[0], per_atom=True)[0]
+print(N)

examples/example_hirsh.py

@@ -1,4 +1,3 @@
-import numpy as np
 import os
 from qstack import compound, fields
 from qstack.fields.decomposition import decompose

examples/example_orcaio.py

@@ -1,48 +1,48 @@
 #!/usr/bin/env python3
 
 import os
+import argparse
 import numpy as np
-import qstack
+from qstack import compound, fields
+from qstack.io import orca
+
 
 def example_orca_reader():
+    parser = argparse.ArgumentParser(description='Example of ORCA reader')
+    parser.add_argument('--name', type=str, default='AJALIH', choices=['AJALIH', 'CEHZOF'], help='job name')
+    args = parser.parse_args()
 
     path = os.path.dirname(os.path.realpath(__file__))
-    name = 'ZOFDAC_6_SPE'
-    openshell = True
-    basis_name = 'def2-TZVP'
-    version = 504
-
-    # In ORCA, def2-TZVP is not sorted wrt angular momenta
-    basis_Fe = [
-[0, [300784.84637, 0.00022806273096], [45088.970557, 0.0017681788761], [10262.516317, 0.009192708349], [2905.2897293, 0.037355495807], [946.11487137, 0.12151108426], [339.87832894, 0.28818881468], [131.94425588, 0.41126612677], [52.111494077, 0.21518583573]],
-[0, [329.48839267, -0.024745216477], [101.92332739, -0.1168308905], [16.240462745, 0.55293621136], [6.8840675801, 0.53601640182]],
-[0, [10.470693782, -0.22912708577], [1.7360039648, 0.71159319984]],
-[0, [0.72577288979, 1.0]],
-[0, [0.11595528203, 1.0]],
-[0, [0.041968227746, 1.0]],
-[1, [1585.395997, 0.0023793960179], [375.38006499, 0.019253154755], [120.31816501, 0.090021836536], [44.788749031, 0.25798172356], [17.829278584, 0.41492649744], [7.2247153786, 0.24207474784]],
-[1, [28.143219756, -0.029041755152], [3.8743241412, 0.55312260343], [1.5410752281, 0.96771136842]],
-[1, [0.5828561525, 1.0]],
-[2, [61.996675034, 0.011971972255], [17.873732552, 0.07321013541], [6.2744782934, 0.23103094314], [2.3552337175, 0.39910706494]],
-[2, [0.85432239901, 1.0]],
-[2, [0.27869254413, 1.0]],
-[1, [0.134915, 1.0]],
-[2, [0.091, 1.0]],
-[3, [1.598, 1.0]]
-]
-    mol0 = qstack.compound.xyz_to_mol(path+'/data/'+name+'.xyz', basis_name, ignore=True)
-    basis = mol0._basis
-    basis['Fe'] = basis_Fe
-    mol = qstack.compound.make_auxmol(mol0, basis)
+    basis = 'def2-TZVP'
+    version = 502
+
+    name = args.name
+    xyz = f'{path}/data/{name}/{name}.xyz'
+    mol = compound.xyz_to_mol(xyz, basis, ecp=basis, parse_comment=True)
+    name_long = f'{name}_{mol.multiplicity}_SPE'
+    gbw = f'{path}/data/{name}/{name}_{mol.multiplicity}_SPE.gbw'
+    openshell = mol.multiplicity > 1
 
     S = mol.intor('int1e_ovlp_sph')
-    D = qstack.orcaio.read_density(mol, name, version=504, openshell=openshell, directory=path+'/data')
+    C, _e, occ = orca.read_gbw(mol, gbw)
+    D1 = fields.dm.make_rdm1(np.squeeze(C), np.squeeze(occ))
+    D2 = orca.read_density(mol, name_long, reorder_dest='pyscf',
+                             version=version, openshell=openshell, directory=f'{path}/data/{name}')
 
     if openshell:
-        print(np.trace(D[0]@S))
-        print(np.trace(D[1]@S))
+        Da, Db = D1
+        print()
+        print(mol.nelec)
+        print(np.trace(Da@S), np.trace(Db@S))
+        Dp, Ds = D2
+        print()
+        print(mol.nelectron, mol.multiplicity-1)
+        print(np.trace(Dp@S), np.trace(Ds@S))
     else:
-        print(np.trace(D@S))
+        print()
+        print(mol.nelectron)
+        print(np.trace(D1@S))
+        print(np.trace(D2@S))
 
 
 if __name__ == '__main__':

examples/metal_contribution.py

@@ -0,0 +1,85 @@
+#!/usr/bin/env python3
+
+import os
+import argparse
+import numpy as np
+import scipy
+from pyscf import scf
+from pyscf.data import elements
+from qstack import compound, fields
+from qstack.io import orca
+
+
+def dipole_moment(mol, dm):
+    charges = mol.atom_charges()
+    coords  = mol.atom_coords()
+    mass = np.round(np.array(elements.MASSES)[charges], 3)
+    mass_center = np.einsum('i,ix->x', mass, coords) / sum(mass)
+    scf.hf.dip_moment(mol, dm, unit='au', origin=mass_center)
+
+
+def find_metal_ao_idx(mol, metals):
+    for iat in range(mol.natm):
+        q = mol.atom_symbol(iat)
+        if q in metals:
+            break
+    metal_ao_limits = mol.offset_ao_by_atom()[iat]
+    return np.arange(metal_ao_limits[2], metal_ao_limits[3])
+
+
+def print_metal_contribution(orb_idx, c, e, metals, do_checks=False):
+    sqrtS = scipy.linalg.sqrtm(mol.intor('int1e_ovlp_sph'))
+    center_ao_idx = find_metal_ao_idx(mol, metals)
+    if do_checks:
+        print(name)
+    else:
+        print(name, end='')
+    for set_id, orb_id in orb_idx:
+        ci = c[set_id,:,orb_id]
+        sqrtSci = sqrtS @ ci
+        fraction_S = np.linalg.norm(sqrtSci[center_ao_idx]) * 100.0
+        if do_checks:
+            print(set_id, orb_id, '\t', "%10.6f"%e[set_id,orb_id], '\t', "%6.1f%%" % fraction_S)
+        else:
+            print('\t', "%6.1f%%" % fraction_S, end='')
+    if not do_checks:
+        print()
+
+
+def get_homo_lumo_idx(occ):
+    c_homo_idx = [(i, np.where(occ[i]>0)[0][-1]) for i in range(occ.shape[0])]
+    c_lumo_idx = [(i, np.where(occ[i]==0)[0][0]) for i in range(occ.shape[0])]
+
+    e_homo = [e[i] for i in c_homo_idx]
+    e_lumo = [e[i] for i in c_lumo_idx]
+
+    homo_idx = c_homo_idx[np.argmax(e_homo)]
+    lumo_idx = c_lumo_idx[np.argmin(e_lumo)]
+    return homo_idx, lumo_idx
+
+
+if __name__=='__main__':
+
+    parser = argparse.ArgumentParser(description='Example of ORCA reader')
+    parser.add_argument('--name', type=str, default='AJALIH', choices=['AJALIH', 'CEHZOF'], help='job name')
+    parser.add_argument('--basis', type=str, default='def2-tZVP', help='basis set')
+    parser.add_argument('--no-checks', dest='do_checks', action='store_false', help='disable checks')
+    args = parser.parse_args()
+
+    metals = ['Cr', 'Mn', 'Fe', 'Co', 'Ni']
+
+    name = args.name
+    path = os.path.dirname(os.path.realpath(__file__))
+    xyz = f'{path}/data/{name}/{name}.xyz'
+    mol = compound.xyz_to_mol(xyz, args.basis, ecp=args.basis, parse_comment=True)
+    gbw = f'{path}/data/{name}/{name}_{mol.multiplicity}_SPE.gbw'
+
+    c, e, occ = orca.read_gbw(mol, gbw)
+    homo_idx, lumo_idx = get_homo_lumo_idx(occ)
+    print_metal_contribution([homo_idx, lumo_idx], c, e, metals, do_checks=args.do_checks)
+
+    if args.do_checks:
+        dm = fields.dm.make_rdm1(np.squeeze(c), np.squeeze(occ))
+        dipole_moment(mol, dm)
+        gap = e[lumo_idx]-e[homo_idx]
+        print(f'{gap=:6f} au ({gap*27.2:6f} eV)')

qstack/compound.py

@@ -6,10 +6,10 @@
 import numpy as np
 from pyscf import gto, data
 from qstack import constants
-from qstack.reorder import get_mrange
-from qstack.mathutils.array import stack_padding, loadtxtvar
+from qstack.mathutils.array import stack_padding, loadtxt_var
 from qstack.mathutils.rotation_matrix import rotate_euler
 from qstack.tools import Cursor
+from .constants import XYZ
 
 
 # detects a charge-spin line, containing only two ints (one positive or negative, the other positive and nonzero)
@@ -144,7 +144,7 @@ def xyz_to_mol(inp, basis="def2-svp", charge=None, spin=None, ignore=False, unit
     if ignore:
         if charge not in (0, None) or spin not in (0, None):
             warnings.warn("Spin and charge values are overwritten", RuntimeWarning, stacklevel=2)
-        atoms = [int(q) if q.isdigit() else data.elements.ELEMENTS_PROTON[q] for q in loadtxtvar(molxyz, dtype='str', usecols=0)]
+        atoms = [int(q) if q.isdigit() else data.elements.ELEMENTS_PROTON[q] for q in loadtxt_var(molxyz, dtype='str', usecols=0)]
         mol.spin = 0
         mol.charge = -(sum(atoms)%2)
     else:
@@ -325,6 +325,9 @@ def basis_flatten(mol, return_both=True, return_shells=False):
         If return_shells is True, also returns:
         - numpy.ndarray: starting AO indices for each shell.
     """
+    def get_mrange(l):
+        return XYZ if l==1 else range(-l, l+1)
+
     x = []
     ao_starts = []
     cursor = Cursor(action='ranger')

qstack/constants.py

@@ -24,3 +24,6 @@
 HARTREE2J = HBAR**2/(MASS_E*(BOHR2ANGS*1e-10)**2)
 HARTREE2EV = 27.21138602
 AU2DEBYE = FUND_CHARGE * BOHR2ANGS*1e-10 / DEBYE  # 2.541746
+
+# m values for p-orbitals ordered as x,y,z
+XYZ = (1, -1, 0)