search

sanshar / Block

Block implements the density matrix renormalization group (DMRG) algorithm for quantum chemistry.
GNU General Public License v3.0
31 stars 33 forks source link

Open-shell calculation with BLOCK. Reading input issues. #55

Open dithillobothrium opened 3 years ago

dithillobothrium commented 3 years ago

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.

gkc1000 commented 3 years ago

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 notifications@github.com 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 https://github.com/sanshar/Block/issues/55, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN5NRCOKLKXPIRE2ZSPZTTSTC3ADANCNFSM4UNH7QWQ .

dithillobothrium commented 3 years ago

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 .

hczhai commented 3 years ago

(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)
gkc1000 commented 3 years ago

@sivkov are you running Block via pyscf? [as per the last reply]

On Thu, Dec 10, 2020 at 7:01 PM hczhai notifications@github.com 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 https://github.com/sanshar/Block/issues/55#issuecomment-742934364, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN5NRGPZV4TVVIWZG724ZTSUGDRXANCNFSM4UNH7QWQ .

dithillobothrium commented 3 years ago

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 .

dithillobothrium commented 3 years ago

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)
hczhai commented 3 years ago

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
dithillobothrium commented 3 years ago

Thank 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