Repository navigation
Open-shell calculation with BLOCK. Reading input issues. #55
Description
Activity
- The DMRG code is (by default) spin adapted. To do the triplet state, you can take the restricted integrals and set the spin=1 in the Block input.…On Fri, Dec 4, 2020 at 2:31 AM sivkov ***@***.***> wrote: Hello, I tried to simulate open-shell system, but BLOCK crashes on reading the FCIDUMP file. For example, I want H2 molecule as open-shell. FCIDUMP file is exported as follows &FCI NORB= 2,NELEC= 2,MS2= 1, ORBSYM=1,1, IUHF=.TRUE. , ISYM=0, / 4.6190068036607312E-01 1 1 1 1 ... 2.0693987712918233E-01 4 4 4 4 ... -2.5214903899133018E-01 4 4 0 0 4.3708714060572051E-01 0 0 0 0 it catches an exeption about wrong indexation in https://github.com/sanshar/StackBlock/blob/f95317b08043b7c531289576d59ad74a6d920741/input.C#L2714 because it creates matrix 2x2, since I have NORB=2 but not 4x4. But, according to agreements NORB should define number of orbitals for one spin direction. In case of open shell the total number of spin orbitals will be 2*NORB. IF I change it to 4 orbitals, I get the error Initial HF occupancy guess: 1 0 0 0 1 0 0 0 Checking input for errors Orbital cannot have an irreducible representation of 0 with c1 symmetry and we see, that the total number of orbitals became 8, so NORB=4 is multiplied somewhere by 2. What is wrong here? Also, I found in the code that o searches for "IUHF", not for "UHF" in the FCIDUMP file. Could you help to resolve these problems? Thanks. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#55>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAN5NRCOKLKXPIRE2ZSPZTTSTC3ADANCNFSM4UNH7QWQ> .
I would take the restricted integrals but the core matrix in my case remains spin-dependent due to some external and periodic effects. So I need UHF kind of FCIDUMP. I see in the code that BLOCK is able to read it but there are some strange bugs or misunderstanding from my side.
The DMRG code is (by default) spin adapted. To do the triplet state, you can take the restricted integrals and set the spin=1 in the Block input.
…
On Fri, Dec 4, 2020 at 2:31 AM sivkov @.**> wrote: Hello, I tried to simulate open-shell system, but BLOCK crashes on reading the FCIDUMP file. For example, I want H2 molecule as open-shell. FCIDUMP file is exported as follows &FCI NORB= 2,NELEC= 2,MS2= 1, ORBSYM=1,1, IUHF=.TRUE. , ISYM=0, / 4.6190068036607312E-01 1 1 1 1 ... 2.0693987712918233E-01 4 4 4 4 ... -2.5214903899133018E-01 4 4 0 0 4.3708714060572051E-01 0 0 0 0 it catches an exeption about wrong indexation in https://github.com/sanshar/StackBlock/blob/f95317b08043b7c531289576d59ad74a6d920741/input.C#L2714 because it creates matrix 2x2, since I have NORB=2 but not 4x4. But, according to agreements NORB should define number of orbitals for one spin direction. In case of open shell the total number of spin orbitals will be 2NORB. IF I change it to 4 orbitals, I get the error Initial HF occupancy guess: 1 0 0 0 1 0 0 0 Checking input for errors Orbital cannot have an irreducible representation of 0 with c1 symmetry and we see, that the total number of orbitals became 8, so NORB=4 is multiplied somewhere by 2. What is wrong here? Also, I found in the code that o searches for "IUHF", not for "UHF" in the FCIDUMP file. Could you help to resolve these problems? Thanks. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#55>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN5NRCOKLKXPIRE2ZSPZTTSTC3ADANCNFSM4UNH7QWQ .(1) It looks like to run Block you need at least 4 spatial orbitals, 2 is too small.
(2) You can use the following python script to generate an UHF FCIDUMP:import numpy as np def write_fcidump(filename, h1e, h2e, nmo, nelec, nuc, ms, isym=1, orbsym=None, tol=1E-15): with open(filename, 'w') as fout: fout.write(' &FCI NORB=%4d,NELEC=%2d,MS2=%d,\n' % (nmo, nelec, ms)) if orbsym is not None and len(orbsym) > 0: fout.write(' ORBSYM=%s,\n' % ','.join([str(x) for x in orbsym])) else: fout.write(' ORBSYM=%s\n' % ('1,' * nmo)) fout.write(' ISYM=%d,\n' % isym) if isinstance(h1e, tuple) and len(h1e) == 2: fout.write(' IUHF=1,\n') assert isinstance(h2e, tuple) if h2e[0].ndim == 4: assert False elif h2e.ndim == 4: assert False fout.write(' &END\n') output_format = '%20.16f%4d%4d%4d%4d\n' npair = nmo * (nmo + 1) // 2 def write_eri(fout, eri): assert eri.ndim in [1, 2, 4] if eri.ndim == 2: # 4-fold symmetry assert(eri.size == npair ** 2) ij = 0 for i in range(nmo): for j in range(0, i + 1): kl = 0 for k in range(0, nmo): for l in range(0, k + 1): if abs(eri[ij, kl]) > tol: fout.write(output_format % (eri[ij, kl], i + 1, j + 1, k + 1, l + 1)) kl += 1 ij += 1 else: # 8-fold symmetry assert(eri.size == npair * (npair + 1) // 2) ij = 0 ijkl = 0 for i in range(nmo): for j in range(0, i + 1): kl = 0 for k in range(0, i + 1): for l in range(0, k + 1): if ij >= kl: if abs(eri[ijkl]) > tol: fout.write(output_format % (eri[ijkl], i + 1, j + 1, k + 1, l + 1)) ijkl += 1 kl += 1 ij += 1 def write_h1e(fout, hx): h = hx.reshape(nmo, nmo) for i in range(nmo): for j in range(0, i + 1): if abs(h[i, j]) > tol: fout.write(output_format % (h[i, j], i + 1, j + 1, 0, 0)) if isinstance(h2e, tuple): assert len(h2e) == 3 vaa, vab, vbb = h2e assert vaa.ndim == vbb.ndim == 1 and vab.ndim == 2 for eri in [vaa, vbb, vab]: write_eri(fout, eri) fout.write(output_format % (0, 0, 0, 0, 0)) assert len(h1e) == 2 for hx in h1e: write_h1e(fout, hx) fout.write(output_format % (0, 0, 0, 0, 0)) fout.write(output_format % (nuc, 0, 0, 0, 0)) else: write_eri(fout, h2e) write_h1e(fout, h1e) fout.write(output_format % (nuc, 0, 0, 0, 0)) if __name__ == "__main__": from pyscf import gto, scf, symm, ao2mo, cc # H chain N = 4 hf_type = 'UHF' BOHR = 0.52917721092 # Angstroms R = 1.8 * BOHR mol = gto.M(atom=[['H', (0, 0, i * R)] for i in range(N)], basis='sto6g', verbose=0, symmetry='d2h') pg = mol.symmetry.lower() if hf_type == "RHF": mf = scf.RHF(mol) elif hf_type == "UHF": mf = scf.UHF(mol) ener = mf.kernel() print("SCF Energy = %20.15f" % ener) if pg == 'd2h': fcidump_sym = ["Ag", "B3u", "B2u", "B1g", "B1u", "B2g", "B3g", "Au"] elif pg == 'c1': fcidump_sym = ["A"] else: assert False mo_coeff = mf.mo_coeff if hf_type == "RHF": n_mo = mo_coeff.shape[1] orb_sym_str = symm.label_orb_symm( mol, mol.irrep_name, mol.symm_orb, mo_coeff) orb_sym = np.array([fcidump_sym.index(i) + 1 for i in orb_sym_str]) h1e = mo_coeff.T @ mf.get_hcore() @ mo_coeff g2e = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff), n_mo) ecore = mol.energy_nuc() na = nb = mol.nelectron // 2 else: mo_coeff_a, mo_coeff_b = mo_coeff[0], mo_coeff[1] n_mo = mo_coeff_b.shape[1] orb_sym_str_a = symm.label_orb_symm( mol, mol.irrep_name, mol.symm_orb, mo_coeff_a) orb_sym_a = np.array([fcidump_sym.index(i) + 1 for i in orb_sym_str_a]) orb_sym = orb_sym_a h1ea = mo_coeff_a.T @ mf.get_hcore() @ mo_coeff_a h1eb = mo_coeff_b.T @ mf.get_hcore() @ mo_coeff_b g2eaa = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff_a), n_mo) g2ebb = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff_b), n_mo) g2eab = ao2mo.kernel( mol, [mo_coeff_a, mo_coeff_a, mo_coeff_b, mo_coeff_b]) h1e = (h1ea, h1eb) g2e = (g2eaa, g2eab, g2ebb) ecore = mol.energy_nuc() na, nb = mol.nelec write_fcidump('FCIDUMP', h1e, g2e, n_mo, na + nb, ecore, na - nb, isym=1, orbsym=orb_sym, tol=1E-15)
- @Sivkov are you running Block via pyscf? [as per the last reply]…On Thu, Dec 10, 2020 at 7:01 PM hczhai ***@***.***> wrote: (1) It looks like to run Block you need at least 4 spatial orbitals, 2 is too small. (2) You can use the following python script to generate an UHF FCIDUMP: import numpy as np def write_fcidump(filename, h1e, h2e, nmo, nelec, nuc, ms, isym=1, orbsym=None, tol=1E-15): with open(filename, 'w') as fout: fout.write(' &FCI NORB=%4d,NELEC=%2d,MS2=%d,\n' % (nmo, nelec, ms)) if orbsym is not None and len(orbsym) > 0: fout.write(' ORBSYM=%s,\n' % ','.join([str(x) for x in orbsym])) else: fout.write(' ORBSYM=%s\n' % ('1,' * nmo)) fout.write(' ISYM=%d,\n' % isym) if isinstance(h1e, tuple) and len(h1e) == 2: fout.write(' IUHF=1,\n') assert isinstance(h2e, tuple) if h2e[0].ndim == 4: assert False elif h2e.ndim == 4: assert False fout.write(' &END\n') output_format = '%20.16f%4d%4d%4d%4d\n' npair = nmo * (nmo + 1) // 2 def write_eri(fout, eri): assert eri.ndim in [1, 2, 4] if eri.ndim == 2: # 4-fold symmetry assert(eri.size == npair ** 2) ij = 0 for i in range(nmo): for j in range(0, i + 1): kl = 0 for k in range(0, nmo): for l in range(0, k + 1): if abs(eri[ij, kl]) > tol: fout.write(output_format % (eri[ij, kl], i + 1, j + 1, k + 1, l + 1)) kl += 1 ij += 1 else: # 8-fold symmetry assert(eri.size == npair * (npair + 1) // 2) ij = 0 ijkl = 0 for i in range(nmo): for j in range(0, i + 1): kl = 0 for k in range(0, i + 1): for l in range(0, k + 1): if ij >= kl: if abs(eri[ijkl]) > tol: fout.write(output_format % (eri[ijkl], i + 1, j + 1, k + 1, l + 1)) ijkl += 1 kl += 1 ij += 1 def write_h1e(fout, hx): h = hx.reshape(nmo, nmo) for i in range(nmo): for j in range(0, i + 1): if abs(h[i, j]) > tol: fout.write(output_format % (h[i, j], i + 1, j + 1, 0, 0)) if isinstance(h2e, tuple): assert len(h2e) == 3 vaa, vab, vbb = h2e assert vaa.ndim == vbb.ndim == 1 and vab.ndim == 2 for eri in [vaa, vbb, vab]: write_eri(fout, eri) fout.write(output_format % (0, 0, 0, 0, 0)) assert len(h1e) == 2 for hx in h1e: write_h1e(fout, hx) fout.write(output_format % (0, 0, 0, 0, 0)) fout.write(output_format % (nuc, 0, 0, 0, 0)) else: write_eri(fout, h2e) write_h1e(fout, h1e) fout.write(output_format % (nuc, 0, 0, 0, 0)) if __name__ == "__main__": from pyscf import gto, scf, symm, ao2mo, cc # H chain N = 4 hf_type = 'UHF' BOHR = 0.52917721092 # Angstroms R = 1.8 * BOHR mol = gto.M(atom=[['H', (0, 0, i * R)] for i in range(N)], basis='sto6g', verbose=0, symmetry='d2h') pg = mol.symmetry.lower() if hf_type == "RHF": mf = scf.RHF(mol) elif hf_type == "UHF": mf = scf.UHF(mol) ener = mf.kernel() print("SCF Energy = %20.15f" % ener) if pg == 'd2h': fcidump_sym = ["Ag", "B3u", "B2u", "B1g", "B1u", "B2g", "B3g", "Au"] elif pg == 'c1': fcidump_sym = ["A"] else: assert False mo_coeff = mf.mo_coeff if hf_type == "RHF": n_mo = mo_coeff.shape[1] orb_sym_str = symm.label_orb_symm( mol, mol.irrep_name, mol.symm_orb, mo_coeff) orb_sym = np.array([fcidump_sym.index(i) + 1 for i in orb_sym_str]) h1e = mo_coeff.T @ mf.get_hcore() @ mo_coeff g2e = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff), n_mo) ecore = mol.energy_nuc() na = nb = mol.nelectron // 2 else: mo_coeff_a, mo_coeff_b = mo_coeff[0], mo_coeff[1] n_mo = mo_coeff_b.shape[1] orb_sym_str_a = symm.label_orb_symm( mol, mol.irrep_name, mol.symm_orb, mo_coeff_a) orb_sym_a = np.array([fcidump_sym.index(i) + 1 for i in orb_sym_str_a]) orb_sym = orb_sym_a h1ea = mo_coeff_a.T @ mf.get_hcore() @ mo_coeff_a h1eb = mo_coeff_b.T @ mf.get_hcore() @ mo_coeff_b g2eaa = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff_a), n_mo) g2ebb = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff_b), n_mo) g2eab = ao2mo.kernel( mol, [mo_coeff_a, mo_coeff_a, mo_coeff_b, mo_coeff_b]) h1e = (h1ea, h1eb) g2e = (g2eaa, g2eab, g2ebb) ecore = mol.energy_nuc() na, nb = mol.nelec write_fcidump('FCIDUMP', h1e, g2e, n_mo, na + nb, ecore, na - nb, isym=1, orbsym=orb_sym, tol=1E-15) — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#55 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAN5NRGPZV4TVVIWZG724ZTSUGDRXANCNFSM4UNH7QWQ> .
I generate FCIDUMP by myself in the HF/DFT code. Then I try to read the FCIDUMP either via PySCF or as standalone. Ideally I dont need PySCF.
@Sivkov are you running Block via pyscf? [as per the last reply]
…
On Thu, Dec 10, 2020 at 7:01 PM hczhai @.***> wrote: (1) It looks like to run Block you need at least 4 spatial orbitals, 2 is too small. (2) You can use the following python script to generate an UHF FCIDUMP: import numpy as np def write_fcidump(filename, h1e, h2e, nmo, nelec, nuc, ms, isym=1, orbsym=None, tol=1E-15): with open(filename, 'w') as fout: fout.write(' &FCI NORB=%4d,NELEC=%2d,MS2=%d,\n' % (nmo, nelec, ms)) if orbsym is not None and len(orbsym) > 0: fout.write(' ORBSYM=%s,\n' % ','.join([str(x) for x in orbsym])) else: fout.write(' ORBSYM=%s\n' % ('1,' * nmo)) fout.write(' ISYM=%d,\n' % isym) if isinstance(h1e, tuple) and len(h1e) == 2: fout.write(' IUHF=1,\n') assert isinstance(h2e, tuple) if h2e[0].ndim == 4: assert False elif h2e.ndim == 4: assert False fout.write(' &END\n') output_format = '%20.16f%4d%4d%4d%4d\n' npair = nmo * (nmo + 1) // 2 def write_eri(fout, eri): assert eri.ndim in [1, 2, 4] if eri.ndim == 2: # 4-fold symmetry assert(eri.size == npair ** 2) ij = 0 for i in range(nmo): for j in range(0, i + 1): kl = 0 for k in range(0, nmo): for l in range(0, k + 1): if abs(eri[ij, kl]) > tol: fout.write(output_format % (eri[ij, kl], i + 1, j + 1, k + 1, l + 1)) kl += 1 ij += 1 else: # 8-fold symmetry assert(eri.size == npair * (npair + 1) // 2) ij = 0 ijkl = 0 for i in range(nmo): for j in range(0, i + 1): kl = 0 for k in range(0, i + 1): for l in range(0, k + 1): if ij >= kl: if abs(eri[ijkl]) > tol: fout.write(output_format % (eri[ijkl], i + 1, j + 1, k + 1, l + 1)) ijkl += 1 kl += 1 ij += 1 def write_h1e(fout, hx): h = hx.reshape(nmo, nmo) for i in range(nmo): for j in range(0, i + 1): if abs(h[i, j]) > tol: fout.write(output_format % (h[i, j], i + 1, j + 1, 0, 0)) if isinstance(h2e, tuple): assert len(h2e) == 3 vaa, vab, vbb = h2e assert vaa.ndim == vbb.ndim == 1 and vab.ndim == 2 for eri in [vaa, vbb, vab]: write_eri(fout, eri) fout.write(output_format % (0, 0, 0, 0, 0)) assert len(h1e) == 2 for hx in h1e: write_h1e(fout, hx) fout.write(output_format % (0, 0, 0, 0, 0)) fout.write(output_format % (nuc, 0, 0, 0, 0)) else: write_eri(fout, h2e) write_h1e(fout, h1e) fout.write(output_format % (nuc, 0, 0, 0, 0)) if name == "main": from pyscf import gto, scf, symm, ao2mo, cc # H chain N = 4 hf_type = 'UHF' BOHR = 0.52917721092 # Angstroms R = 1.8 * BOHR mol = gto.M(atom=[['H', (0, 0, i * R)] for i in range(N)], basis='sto6g', verbose=0, symmetry='d2h') pg = mol.symmetry.lower() if hf_type == "RHF": mf = scf.RHF(mol) elif hf_type == "UHF": mf = scf.UHF(mol) ener = mf.kernel() print("SCF Energy = %20.15f" % ener) if pg == 'd2h': fcidump_sym = ["Ag", "B3u", "B2u", "B1g", "B1u", "B2g", "B3g", "Au"] elif pg == 'c1': fcidump_sym = ["A"] else: assert False mo_coeff = mf.mo_coeff if hf_type == "RHF": n_mo = mo_coeff.shape[1] orb_sym_str = symm.label_orb_symm( mol, mol.irrep_name, mol.symm_orb, mo_coeff) orb_sym = np.array([fcidump_sym.index(i) + 1 for i in orb_sym_str]) h1e = mo_coeff.T @ mf.get_hcore() @ mo_coeff g2e = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff), n_mo) ecore = mol.energy_nuc() na = nb = mol.nelectron // 2 else: mo_coeff_a, mo_coeff_b = mo_coeff[0], mo_coeff[1] n_mo = mo_coeff_b.shape[1] orb_sym_str_a = symm.label_orb_symm( mol, mol.irrep_name, mol.symm_orb, mo_coeff_a) orb_sym_a = np.array([fcidump_sym.index(i) + 1 for i in orb_sym_str_a]) orb_sym = orb_sym_a h1ea = mo_coeff_a.T @ mf.get_hcore() @ mo_coeff_a h1eb = mo_coeff_b.T @ mf.get_hcore() @ mo_coeff_b g2eaa = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff_a), n_mo) g2ebb = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff_b), n_mo) g2eab = ao2mo.kernel( mol, [mo_coeff_a, mo_coeff_a, mo_coeff_b, mo_coeff_b]) h1e = (h1ea, h1eb) g2e = (g2eaa, g2eab, g2ebb) ecore = mol.energy_nuc() na, nb = mol.nelec write_fcidump('FCIDUMP', h1e, g2e, n_mo, na + nb, ecore, na - nb, isym=1, orbsym=orb_sym, tol=1E-15) — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#55 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN5NRGPZV4TVVIWZG724ZTSUGDRXANCNFSM4UNH7QWQ .I generate FCIDUMP by myself. I think I need some example of UHF fcidump. But I will try to take more orbitals.
(1) It looks like to run Block you need at least 4 spatial orbitals, 2 is too small.
(2) You can use the following python script to generate an UHF FCIDUMP:import numpy as np def write_fcidump(filename, h1e, h2e, nmo, nelec, nuc, ms, isym=1, orbsym=None, tol=1E-15): with open(filename, 'w') as fout: fout.write(' &FCI NORB=%4d,NELEC=%2d,MS2=%d,\n' % (nmo, nelec, ms)) if orbsym is not None and len(orbsym) > 0: fout.write(' ORBSYM=%s,\n' % ','.join([str(x) for x in orbsym])) else: fout.write(' ORBSYM=%s\n' % ('1,' * nmo)) fout.write(' ISYM=%d,\n' % isym) if isinstance(h1e, tuple) and len(h1e) == 2: fout.write(' IUHF=1,\n') assert isinstance(h2e, tuple) if h2e[0].ndim == 4: assert False elif h2e.ndim == 4: assert False fout.write(' &END\n') output_format = '%20.16f%4d%4d%4d%4d\n' npair = nmo * (nmo + 1) // 2 def write_eri(fout, eri): assert eri.ndim in [1, 2, 4] if eri.ndim == 2: # 4-fold symmetry assert(eri.size == npair ** 2) ij = 0 for i in range(nmo): for j in range(0, i + 1): kl = 0 for k in range(0, nmo): for l in range(0, k + 1): if abs(eri[ij, kl]) > tol: fout.write(output_format % (eri[ij, kl], i + 1, j + 1, k + 1, l + 1)) kl += 1 ij += 1 else: # 8-fold symmetry assert(eri.size == npair * (npair + 1) // 2) ij = 0 ijkl = 0 for i in range(nmo): for j in range(0, i + 1): kl = 0 for k in range(0, i + 1): for l in range(0, k + 1): if ij >= kl: if abs(eri[ijkl]) > tol: fout.write(output_format % (eri[ijkl], i + 1, j + 1, k + 1, l + 1)) ijkl += 1 kl += 1 ij += 1 def write_h1e(fout, hx): h = hx.reshape(nmo, nmo) for i in range(nmo): for j in range(0, i + 1): if abs(h[i, j]) > tol: fout.write(output_format % (h[i, j], i + 1, j + 1, 0, 0)) if isinstance(h2e, tuple): assert len(h2e) == 3 vaa, vab, vbb = h2e assert vaa.ndim == vbb.ndim == 1 and vab.ndim == 2 for eri in [vaa, vbb, vab]: write_eri(fout, eri) fout.write(output_format % (0, 0, 0, 0, 0)) assert len(h1e) == 2 for hx in h1e: write_h1e(fout, hx) fout.write(output_format % (0, 0, 0, 0, 0)) fout.write(output_format % (nuc, 0, 0, 0, 0)) else: write_eri(fout, h2e) write_h1e(fout, h1e) fout.write(output_format % (nuc, 0, 0, 0, 0)) if __name__ == "__main__": from pyscf import gto, scf, symm, ao2mo, cc # H chain N = 4 hf_type = 'UHF' BOHR = 0.52917721092 # Angstroms R = 1.8 * BOHR mol = gto.M(atom=[['H', (0, 0, i * R)] for i in range(N)], basis='sto6g', verbose=0, symmetry='d2h') pg = mol.symmetry.lower() if hf_type == "RHF": mf = scf.RHF(mol) elif hf_type == "UHF": mf = scf.UHF(mol) ener = mf.kernel() print("SCF Energy = %20.15f" % ener) if pg == 'd2h': fcidump_sym = ["Ag", "B3u", "B2u", "B1g", "B1u", "B2g", "B3g", "Au"] elif pg == 'c1': fcidump_sym = ["A"] else: assert False mo_coeff = mf.mo_coeff if hf_type == "RHF": n_mo = mo_coeff.shape[1] orb_sym_str = symm.label_orb_symm( mol, mol.irrep_name, mol.symm_orb, mo_coeff) orb_sym = np.array([fcidump_sym.index(i) + 1 for i in orb_sym_str]) h1e = mo_coeff.T @ mf.get_hcore() @ mo_coeff g2e = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff), n_mo) ecore = mol.energy_nuc() na = nb = mol.nelectron // 2 else: mo_coeff_a, mo_coeff_b = mo_coeff[0], mo_coeff[1] n_mo = mo_coeff_b.shape[1] orb_sym_str_a = symm.label_orb_symm( mol, mol.irrep_name, mol.symm_orb, mo_coeff_a) orb_sym_a = np.array([fcidump_sym.index(i) + 1 for i in orb_sym_str_a]) orb_sym = orb_sym_a h1ea = mo_coeff_a.T @ mf.get_hcore() @ mo_coeff_a h1eb = mo_coeff_b.T @ mf.get_hcore() @ mo_coeff_b g2eaa = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff_a), n_mo) g2ebb = ao2mo.restore(8, ao2mo.kernel(mol, mo_coeff_b), n_mo) g2eab = ao2mo.kernel( mol, [mo_coeff_a, mo_coeff_a, mo_coeff_b, mo_coeff_b]) h1e = (h1ea, h1eb) g2e = (g2eaa, g2eab, g2ebb) ecore = mol.energy_nuc() na, nb = mol.nelec write_fcidump('FCIDUMP', h1e, g2e, n_mo, na + nb, ecore, na - nb, isym=1, orbsym=orb_sym, tol=1E-15)
Here is an example of (H4-d2h symmetry) UHF FCIDUMP.
It has 6 sections for the integrals, V-alpha-alpha, V-beta-beta, V-alpha-beta, T-alpha, T-beta, and Const, separated by "0.0 0 0 0 0". See the previous script for details.&FCI NORB= 4,NELEC= 4,MS2=0, ORBSYM=1,5,1,5, ISYM=1, IUHF=1, &END 0.5083007829564030 1 1 1 1 0.1574060863159039 2 1 2 1 0.4461364581632330 2 2 1 1 0.4643541373866802 2 2 2 2 -0.0835460048314562 3 1 1 1 0.0084520980677123 3 1 2 2 0.1076814510738056 3 1 3 1 0.0993405121123775 3 2 2 1 0.1374351924433553 3 2 3 2 0.4572693492091709 3 3 1 1 0.4580462798315313 3 3 2 2 -0.0100607334199195 3 3 3 1 0.4784591094144346 3 3 3 3 0.0438882947078274 4 1 2 1 -0.0506041964088772 4 1 3 2 0.0961125852495129 4 1 4 1 0.0863400560690246 4 2 1 1 0.0061930898726832 4 2 2 2 -0.0976978440469843 4 2 3 1 0.0057796879060008 4 2 3 3 0.1041641131899388 4 2 4 2 -0.1499542552157829 4 3 2 1 -0.1001455350087904 4 3 3 2 -0.0417868235594825 4 3 4 1 0.1620733641164259 4 3 4 3 0.5356532301946455 4 4 1 1 0.4759994469592082 4 4 2 2 -0.0883961794249322 4 4 3 1 0.4937433605905979 4 4 3 3 0.0964548805904054 4 4 4 2 0.5984395584964084 4 4 4 4 0.0000000000000000 0 0 0 0 0.5083004232237522 1 1 1 1 0.1574062438620759 2 1 2 1 0.4461363655544462 2 2 1 1 0.4643541214908927 2 2 2 2 -0.0835458279360813 3 1 1 1 0.0084522362699454 3 1 2 2 0.1076816092801142 3 1 3 1 0.0993406193092624 3 2 2 1 0.1374351070956478 3 2 3 2 0.4572695074154798 3 3 1 1 0.4580462542172687 3 3 2 2 -0.0100609424382043 3 3 3 1 0.4784591527344686 3 3 3 3 0.0438881181472035 4 1 2 1 -0.0506042608952912 4 1 3 2 0.0961125516109522 4 1 4 1 0.0863397882936577 4 2 1 1 0.0061929487198046 4 2 2 2 -0.0976978686202398 4 2 3 1 0.0057798753338584 4 2 3 3 0.1041639973511403 4 2 4 2 -0.1499543197021969 4 3 2 1 -0.1001453349270675 4 3 3 2 -0.0417868812506004 4 3 4 1 0.1620733255565223 4 3 4 3 0.5356531506914772 4 4 1 1 0.4759993311204100 4 4 2 2 -0.0883963499205894 4 4 3 1 0.4937435583168170 4 4 3 3 0.0964549357030432 4 4 4 2 0.5984398060697933 4 4 4 4 0.0000000000000000 0 0 0 0 0.5083006030898282 1 1 1 1 0.4461363473576081 1 1 2 2 -0.0835460597641479 1 1 3 1 0.4572695290757468 1 1 3 3 0.0863399986276825 1 1 4 2 0.5356533410002735 1 1 4 4 0.1574061650888555 2 1 2 1 0.0993404388951413 2 1 3 2 0.0438882342939144 2 1 4 1 -0.1499542387144289 2 1 4 3 0.4461364763598014 2 2 1 1 0.4643541294387007 2 2 2 2 0.0084521108880433 2 2 3 1 0.4580462616349640 2 2 3 3 0.0061930824000966 2 2 4 2 0.4759994549071901 2 2 4 4 -0.0835457730035601 3 1 1 1 0.0084522234495031 3 1 2 2 0.1076815301769181 3 1 3 1 -0.0100609652478155 3 1 3 3 -0.0976977819011402 3 1 4 2 -0.0883963048067246 3 1 4 4 0.0993406925265315 3 2 2 1 0.1374351497694806 3 2 3 2 -0.0506042404656113 3 2 4 1 -0.1001453923461528 3 2 4 3 0.4572693275494037 3 3 1 1 0.4580462724141067 3 3 2 2 -0.0100607106101379 3 3 3 1 0.4784591310742020 3 3 3 3 0.0057796649998342 3 3 4 2 0.4937433680080252 3 3 4 4 0.0438881785611532 4 1 2 1 -0.0506042168385719 4 1 3 2 0.0961125684302098 4 1 4 1 -0.0417869594918433 4 1 4 3 0.0863398457348663 4 2 1 1 0.0061929561923169 4 2 2 2 -0.0976979307661215 4 2 3 1 0.0057798982401591 4 2 3 3 0.1041640552705168 4 2 4 2 0.0964550142707728 4 2 4 4 -0.1499543362037927 4 3 2 1 -0.1001454775897001 4 3 3 2 -0.0417867453182301 4 3 4 1 0.1620733448363552 4 3 4 3 0.5356530398861216 4 4 1 1 0.4759993231726017 4 4 2 2 -0.0883962245386873 4 4 3 1 0.4937435508991220 4 4 3 3 0.0964548020227513 4 4 4 2 0.5984396822830154 4 4 4 4 0.0000000000000000 0 0 0 0 -1.8949792637704244 1 1 0 0 -1.5938622714227224 2 2 0 0 0.1659819525071990 3 1 0 0 -1.2671440860563838 3 3 0 0 -0.1349843015574928 4 2 0 0 -0.8813389499224603 4 4 0 0 0.0000000000000000 0 0 0 0 -1.8949789064265117 1 1 0 0 -1.5938620981885263 2 2 0 0 0.1659826283412012 3 1 0 0 -1.2671444434002974 3 3 0 0 -0.1349847587705340 4 2 0 0 -0.8813391231566573 4 4 0 0 0.0000000000000000 0 0 0 0 2.4074074074074074 0 0 0 0Thank you for the answer. The format is that the index of spin-down orbital is NOT 2x of the index of a spin-up orbital. The blocks are just followed one after the other. Ok, this is clear, there are different agreements how to export UHF FCIDUMP and I used indices of spin-down orbitals as 2 x indices of spin-up orbitals. It was wrong.
Here is an example of (H4-d2h symmetry) UHF FCIDUMP.
It has 6 sections for the integrals, V-alpha-alpha, V-beta-beta, V-alpha-beta, T-alpha, T-beta, and Const, separated by "0.0 0 0 0 0". See the previous script for details.&FCI NORB= 4,NELEC= 4,MS2=0, ORBSYM=1,5,1,5, ISYM=1, IUHF=1, &END 0.5083007829564030 1 1 1 1 0.1574060863159039 2 1 2 1 0.4461364581632330 2 2 1 1 0.4643541373866802 2 2 2 2 -0.0835460048314562 3 1 1 1 0.0084520980677123 3 1 2 2 0.1076814510738056 3 1 3 1 0.0993405121123775 3 2 2 1 0.1374351924433553 3 2 3 2 0.4572693492091709 3 3 1 1 0.4580462798315313 3 3 2 2 -0.0100607334199195 3 3 3 1 0.4784591094144346 3 3 3 3 0.0438882947078274 4 1 2 1 -0.0506041964088772 4 1 3 2 0.0961125852495129 4 1 4 1 0.0863400560690246 4 2 1 1 0.0061930898726832 4 2 2 2 -0.0976978440469843 4 2 3 1 0.0057796879060008 4 2 3 3 0.1041641131899388 4 2 4 2 -0.1499542552157829 4 3 2 1 -0.1001455350087904 4 3 3 2 -0.0417868235594825 4 3 4 1 0.1620733641164259 4 3 4 3 0.5356532301946455 4 4 1 1 0.4759994469592082 4 4 2 2 -0.0883961794249322 4 4 3 1 0.4937433605905979 4 4 3 3 0.0964548805904054 4 4 4 2 0.5984395584964084 4 4 4 4 0.0000000000000000 0 0 0 0 0.5083004232237522 1 1 1 1 0.1574062438620759 2 1 2 1 0.4461363655544462 2 2 1 1 0.4643541214908927 2 2 2 2 -0.0835458279360813 3 1 1 1 0.0084522362699454 3 1 2 2 0.1076816092801142 3 1 3 1 0.0993406193092624 3 2 2 1 0.1374351070956478 3 2 3 2 0.4572695074154798 3 3 1 1 0.4580462542172687 3 3 2 2 -0.0100609424382043 3 3 3 1 0.4784591527344686 3 3 3 3 0.0438881181472035 4 1 2 1 -0.0506042608952912 4 1 3 2 0.0961125516109522 4 1 4 1 0.0863397882936577 4 2 1 1 0.0061929487198046 4 2 2 2 -0.0976978686202398 4 2 3 1 0.0057798753338584 4 2 3 3 0.1041639973511403 4 2 4 2 -0.1499543197021969 4 3 2 1 -0.1001453349270675 4 3 3 2 -0.0417868812506004 4 3 4 1 0.1620733255565223 4 3 4 3 0.5356531506914772 4 4 1 1 0.4759993311204100 4 4 2 2 -0.0883963499205894 4 4 3 1 0.4937435583168170 4 4 3 3 0.0964549357030432 4 4 4 2 0.5984398060697933 4 4 4 4 0.0000000000000000 0 0 0 0 0.5083006030898282 1 1 1 1 0.4461363473576081 1 1 2 2 -0.0835460597641479 1 1 3 1 0.4572695290757468 1 1 3 3 0.0863399986276825 1 1 4 2 0.5356533410002735 1 1 4 4 0.1574061650888555 2 1 2 1 0.0993404388951413 2 1 3 2 0.0438882342939144 2 1 4 1 -0.1499542387144289 2 1 4 3 0.4461364763598014 2 2 1 1 0.4643541294387007 2 2 2 2 0.0084521108880433 2 2 3 1 0.4580462616349640 2 2 3 3 0.0061930824000966 2 2 4 2 0.4759994549071901 2 2 4 4 -0.0835457730035601 3 1 1 1 0.0084522234495031 3 1 2 2 0.1076815301769181 3 1 3 1 -0.0100609652478155 3 1 3 3 -0.0976977819011402 3 1 4 2 -0.0883963048067246 3 1 4 4 0.0993406925265315 3 2 2 1 0.1374351497694806 3 2 3 2 -0.0506042404656113 3 2 4 1 -0.1001453923461528 3 2 4 3 0.4572693275494037 3 3 1 1 0.4580462724141067 3 3 2 2 -0.0100607106101379 3 3 3 1 0.4784591310742020 3 3 3 3 0.0057796649998342 3 3 4 2 0.4937433680080252 3 3 4 4 0.0438881785611532 4 1 2 1 -0.0506042168385719 4 1 3 2 0.0961125684302098 4 1 4 1 -0.0417869594918433 4 1 4 3 0.0863398457348663 4 2 1 1 0.0061929561923169 4 2 2 2 -0.0976979307661215 4 2 3 1 0.0057798982401591 4 2 3 3 0.1041640552705168 4 2 4 2 0.0964550142707728 4 2 4 4 -0.1499543362037927 4 3 2 1 -0.1001454775897001 4 3 3 2 -0.0417867453182301 4 3 4 1 0.1620733448363552 4 3 4 3 0.5356530398861216 4 4 1 1 0.4759993231726017 4 4 2 2 -0.0883962245386873 4 4 3 1 0.4937435508991220 4 4 3 3 0.0964548020227513 4 4 4 2 0.5984396822830154 4 4 4 4 0.0000000000000000 0 0 0 0 -1.8949792637704244 1 1 0 0 -1.5938622714227224 2 2 0 0 0.1659819525071990 3 1 0 0 -1.2671440860563838 3 3 0 0 -0.1349843015574928 4 2 0 0 -0.8813389499224603 4 4 0 0 0.0000000000000000 0 0 0 0 -1.8949789064265117 1 1 0 0 -1.5938620981885263 2 2 0 0 0.1659826283412012 3 1 0 0 -1.2671444434002974 3 3 0 0 -0.1349847587705340 4 2 0 0 -0.8813391231566573 4 4 0 0 0.0000000000000000 0 0 0 0 2.4074074074074074 0 0 0 0
Hello,
I tried to simulate open-shell system, but BLOCK crashes on reading the FCIDUMP file.
For example, I want H2 molecule as open-shell. FCIDUMP file is exported as follows
it catches an exeption about wrong indexation in
https://github.com/sanshar/StackBlock/blob/f95317b08043b7c531289576d59ad74a6d920741/input.C#L2714
because it creates matrix 2x2, since I have NORB=2 but not 4x4.
But, according to agreements NORB should define number of orbitals for one spin direction. In case of open shell the total number of spin orbitals will be 2*NORB.
IF I change it to 4 orbitals, I get the error
and we see, that the total number of orbitals became 8, so NORB=4 is multiplied somewhere by 2. What is wrong here?
Also, I found in the code that o searches for "IUHF", not for "UHF" in the FCIDUMP file.
Could you help to resolve these problems?
Thanks.