The MMSE Sandbox

Question

Similar to this sandbox, you're encouraged to write drafts of your questions here before posting them on the main site, if you think this may be beneficial. I was tempted to come up with a set of "rules" for this, or to simply say "don't over use this", but my guess is that this won't be necessary. Also, if you "over use" this, bringing more traffic into our Meta site might not be such a bad thing :)

If you have an answer to any of these questions, please leave a comment and the question can be moved from The MMSE Sandbox to the main site, if the asker wishes.

Answer 1

Can this procedure for multiply augmenting a basis set with the "basis_set_exchange" package be improved?

python,basis-sets,basis-set-exchange,quantum-chemistry,electronic-structure

1) Install the basis_set_exchange package

git clone and pip install -e are required because the default version that's installed with pip does not include the basis_set_exchange.manip.geometric_augmentation() function for multiply augmenting a basis set:

git clone https://github.com/MolSSI-BSE/basis_set_exchange.git
pip install -e basis_set_exchange

2) Create a basis set "dictionary" for the original basis set

We will multiply augment the aug-cc-pVDZ basis set for the H atom. We can also do it for other atoms by using something like elements=[6,8] or even elements=['C', 8, 'Ne', '16']:

bs_dict = basis_set_exchange.get_basis('aug-cc-pVDZ',elements=1)

3) Add the extra diffuse functions:

The 1 indicates that we're constructing a daug basis set, but we could use 2 and get a taug basis set, or 3 and get a qaug basis set:

daug=basis_set_exchange.manip.geometric_augmentation(bs_dict, 1, use_copy=True, as_component=False, steep=False)

4) Get the new basis set into the format we desire

In this case I'll use the Guassian basis set format:

daug_in_gaussian_format=basis_set_exchange.write_formatted_basis_str(daug, 'gaussian94', header=None)

4) Print the new basis set to a file!

print(daug_in_gaussian_format,  file=open('daug-cc-pVDZ.gbs', 'w')

Answer 2

Generated Fort.55 File in PySCF for MRCC, how to handle symmetry from UHF reference calculation?

Background

I am trying to build an interface between PySCF and MRCC using UHF reference. The basic idea of this interface is to generate 1e, 2e integrals using PySCF, save those integrals with ASCII format in what it is known as fort.55 file in MRCC or FCIDUMP file in PySCF, the format of this file follows the UHF standard designed by Knowles and Handy, in which after the header part, separate 4-index and 2-index integrals are written for alpha-alpha, beta-beta, and alph-beta orbitals, with a line of zeroes separating them. Unfortunately, this format and despite being used in some major quantum chemistry codes like MRCC or CFOUR, is not natively supported in PySCF.

Progress

I have implied this format to PySCF FCIDUMP file while taking care of the header part such that the user can generate and save the designated fort.55 file which can be later used as an input for MRCC in order to perform state-of-the-art post-HF calculations at the level of FCI or CC(n).

Code

The following are details about the script I am using to build the fort.55 file. Suppose that we have defined a Mol object for open-shell system, followed by running a UHF calculation. In the next step, I convert the 2e integrals using 4-fold permutation symmetry for alph-alpha, beta-beta, and alpha-beta orbitals which is saved in corresponding eri_aaaa, eri_bbbb, and eri_aabb, respectively. Moreover, the transformed 1e core Hamiltonian for alpha and beta are also calculated and stored in h_aa and h_bb, respectively.

orbs = mf.mo_coeff
nmo = orbs[0].shape[0]

eri_aaaa = pyscf.ao2mo.restore(4,pyscf.ao2mo.incore.general(mf._eri, (orbs[0],orbs[0],orbs[0],orbs[0]), compact=False),nmo)
eri_bbbb = pyscf.ao2mo.restore(4,pyscf.ao2mo.incore.general(mf._eri, (orbs[1],orbs[1],orbs[1],orbs[1]), compact=False),nmo)
eri_aabb = pyscf.ao2mo.restore(4,pyscf.ao2mo.incore.general(mf._eri, (orbs[0],orbs[0],orbs[1],orbs[1]), compact=False),nmo)

h_core = mf.get_hcore(mol)
h_aa = reduce(numpy.dot, (orbs[0].T, h_core, orbs[0]))
h_bb = reduce(numpy.dot, (orbs[1].T, h_core, orbs[1]))
nuc = mol.energy_nuc()

In the next section of the code, I try to handle the symmetry of the orbitals in a way that will match the format of fort.55, I also define the range of values for integral indices.

if mol.symmetry:
        groupname = mol.groupname
        if groupname in ('SO3', 'Dooh'):
            logger.info(mol, 'Lower symmetry from %s to D2h', groupname)
            raise RuntimeError('Lower symmetry from %s to D2h' % groupname)
        elif groupname == 'Coov':
            logger.info(mol, 'Lower symmetry from Coov to C2v')
            raise RuntimeError('''Lower symmetry from Coov to C2v''')
orbsym = pyscf.symm.label_orb_symm(mol,mol.irrep_name,mol.symm_orb,orbs[0])
orbsym = numpy.array(orbsym)
orbsym = [param.IRREP_ID_TABLE[groupname][i]+1 for i in orbsym]
a_inds = [i+1 for i in range(orbs[0].shape[0])]
b_inds = [i+1 for i in range(orbs[1].shape[1])]
nelec = mol.nelec
tol=1e-18

The last part is considered the main part of the script, which is supposed to handle the header of the generated file in a way that will match the fort.55 format, it is supposed to deal with the 4-index, and 2-index integrals in which I am assuming that we are using a 4-fold permutation symmetry. The sections after '4-fold symmetry' are in the sequence of alpha-alpha, beta-beta, alpha-beta 2e integrals, followed by alpha, beta 1e integrals, and finally the nuclear repulsion energy.

with open('fort.55', 'w') as fout:
        if not isinstance(nelec, (int, numpy.number)):
            ms    = abs(nelec[0] - nelec[1])
            nelec =     nelec[0]  + nelec[1]
        else: ms=0
        fout.write(f"{nmo:1d} {nelec:1d}\n")
        if orbsym is not None and len(orbsym) > 0:
            fout.write(f"{' '.join([str(x) for x in orbsym])}\n")
        else:
            fout.write(f"{' 1' * nmo}\n")
        fout.write(' 150000\n')
        output_format = float_format + ' %5d %5d %5d %5d\n'
        #4-fold symmetry
        kl = 0
        for l in range(nmo):
            for k in range(0, l+1):
                ij = 0
                for i in range(0, nmo):
                    for j in range(0, i+1):
                        if i >= k:
                            if abs(eri_aaaa[ij,kl]) > tol:
                                fout.write(output_format % (eri_aaaa[ij,kl], a_inds[i], a_inds[j], a_inds[k], a_inds[l]))
                        ij += 1
                kl += 1
        fout.write(' 0.00000000000000000000E+00' + '     0     0     0     0\n')
        kl = 0
        for l in range(nmo):
            for k in range(0, l+1):
                ij = 0
                for i in range(0, nmo):
                    for j in range(0, i+1):
                        if i >= k:
                            if abs(eri_bbbb[ij,kl]) > tol:
                                fout.write(output_format % (eri_bbbb[ij,kl], b_inds[i], b_inds[j], b_inds[k], b_inds[l]))
                        ij += 1
                kl += 1
        fout.write(' 0.00000000000000000000E+00' + '     0     0     0     0\n')
        ij = 0
        for j in range(nmo):
            for i in range(0, j+1):
                kl = 0
                for k in range(nmo):
                    for l in range(0, k+1):
                        if abs(eri_aabb[ij,kl]) > tol:
                            fout.write(output_format % (eri_aabb[ij,kl], a_inds[i], a_inds[j], b_inds[k], b_inds[l]))
                        kl += 1
                ij +=1

        fout.write(' 0.00000000000000000000E+00' + '     0     0     0     0\n')
        h_aa = h_aa.reshape(nmo,nmo)
        h_bb = h_bb.reshape(nmo,nmo)
        output_format = float_format + ' %5d %5d     0     0\n'
        for i in range(nmo):
            for j in range(nmo):
                fout.write(output_format % (h_aa[i,j], a_inds[i], a_inds[j]))
        fout.write(' 0.00000000000000000000E+00' + '     0     0     0     0\n')
        for i in range(nmo):
            for j in range(nmo):
                fout.write(output_format % (h_bb[i,j], b_inds[i], b_inds[j]))
        fout.write(' 0.00000000000000000000E+00' + '     0     0     0     0\n')
        output_format = float_format + '     0     0     0     0\n'
        fout.write(output_format % nuc)

Test & Debugging

In order to test the above method and further arrive at my question, I will consider two open-shell systems, namely Be(+) and Be(-). The former having 3 electrons occupying s-type orbitals, while the latter having 5 electrons with s- and p-type orbitals being occupied. by using the generated fort.55 file as described above with the following input file in MRCC:

basis=STO-XG-EMSL
iface=cfour
uncontract=off
calc=CC(5)
core=corr
itol=18
scftol=13
cctol=7
ccmaxit=999
scfmaxit=9999
scfiguess=ao
scftype=UHF
rohftype=semicanonical
rest=2
charge=-1
refdet=serialno
1,2
3

symm=0
mult=2
occ=2,0,0,0,0,0,0,1/2,0,0,0,0,0,0,0
geom
BE
tprint=0.01
verbosity=3 

You can see that I have turned off the symmetry 'symm=0', upon which, the calculation will run normally, to confirm that, I compare the energies I get with those generated using standalone MRCC, and I can confirm that the energies are identical in that case for both Be(-) and Be(+).

Those values are for Be(-) using STO-XG
mrcc-mrcc: 
CCSD        -14.274723492321
CCSD[T]     -14.274740253256
CCSD(T)     -14.274739121789
CCSDT       -14.274741793166
CCSDT[Q]    -14.274743492727
CCSDT(Q)/A  -14.274743594299
CCSDT(Q)/B  -14.274743590678
CCSDTQ      -14.274743599645
CCSDTQ[P]   -14.274743599645
CCSDTQ(P)/A -14.274743599645
CCSDTQ(P)/B -14.274743599645
CCSDTQP     -14.274743599073
pyscf-mrcc (symmetry is off): 
CCSD        -14.274723492321
CCSD[T]     -14.274740253256
CCSD(T)     -14.274739121789
CCSDT       -14.274741793166
CCSDT[Q]    -14.274743492727
CCSDT(Q)/A  -14.274743594299
CCSDT(Q)/B  -14.274743590678
CCSDTQ      -14.274743599645
CCSDTQ[P]   -14.274743599645
CCSDTQ(P)/A -14.274743599645
CCSDTQ(P)/B -14.274743599645
CCSDTQP     -14.274743599073

The Problem

Running the above calculation without using symmetry is affordable in terms of computational power and time for small systems with small basis sets, however, my attempt for such interface is to perform higher-order coupled cluster calculations such as CCSDT(Q), CCSDTQ, CCSDTQ(P), CCSDTQP, etc. with large basis sets such as aCV9Z. Thus, symmetry shall be applied as well.

I will start with Be(+) which has 3 electrons occupying s-type orbitals, using the generated fort.55 file in PySCF with STO-XG basis set, which in this case has the following ORBSYM: 'Ag' 'Ag' 'B1u' 'B2u' 'B3u' 'Ag', I should be using 'symm=1' in the input file of MRCC, I find that energies are identical to the energies generated with standalone MRCC.

Be+ (STO-XG)
pyscf-mrcc (symmetry is on): 
CCSD    -14.101017848371
CCSD[T] -14.101018583402
CCSD(T) -14.101018579146
CCSDT   -14.101018557319
mrcc-mrcc
CCSD    -14.101017848371
CCSD[T] -14.101018583402
CCSD(T) -14.101018579146
CCSDT   -14.101018557319

On the other hand, in the case of Be(-) which has 5 electrons occupying s- and p-type orbitals, the generated fort.55 file in PySCF would have the following ORBSYM: Ag' 'Ag' 'B3u' 'B1u' 'B2u' 'Ag', I should be using 'symm=8' in this case. However, the final energies are found to be different (lower) than the energies generated with standalone MRCC.

pyscf-mrcc (symmetry is on):
CCSD        -14.256521081042
CCSD[T]     -14.256521081042
CCSD(T)     -14.256521081042
CCSDT       -14.256521081190
CCSDT[Q]    -14.256521081190
CCSDT(Q)/A  -14.256521081190
CCSDT(Q)/B  -14.256521081190
CCSDTQ      -14.256521081240
CCSDTQ[P]   -14.256521081240
CCSDTQ(P)/A -14.256521081240
CCSDTQ(P)/B -14.256521081240
CCSDTQP     -14.256521081257

You can find full files for the inputs and outputs generated with PySCF and MRCC in the following Input_Output_Files.

The Question

Based on the above description, what is the reason that Be(+) works correctly with symmetry while not for Be(-)? What could possibly be wrong/missing in my code in order to correctly handle the symmetry of spatial orbitals in systems where p-type orbitals are occupied?

Answer 3

VibRot in MOLCAS doesn't find all levels

bug, molcas, spectroscopy, vibrations, software-assistance

The (7,7) isotopologue of the b-state of Li₂ supports over 90 vibrational levels, as you can see from Table IV of my 2015 paper, or lines 506 to 607 of the output of my LEVEL program (the corresponding input can be found on lines 365 to 467 of this file).

Using the VibRot program in OpenMOLCAS, I was only able to find about 80 vibrational levels (see here for an output file that contains a printout of the input file at the top). Using the Scale keyword didn't change things. The range (from 0.3 to 3000 in a.u.) of internuclear distances for numerical integration already has an extremely low lower bound and extremely high upper bound, and the Step value (Step = 0.000000001 a.u.) for energies is already extremely low. Changing the Range or the Step values also didn't improve things. I'll be more specific:

• Increasing the lower bound of the integration range from 0.3 to 0.4 (ironically, since an increase in the lower bound should make the numerical solution less accurate) gave me 2 more vibrational levels, but increasing it further to 0.5 or 0.6 didn't give me more levels.
• Decreasing the upper bound of the integration range from 3000 to 1000 gave me fewer levels, but increasing it to 4000 caused the program to crash prematurely and making it 3500 also gave me fewer levels. I could try 2000 and other values but the calculations take extremely long (for a rovibrational solver) and this type of trial-and-error work is unusual for trying to find over a dozen missing levels (perhaps I need to explicitly incorporate the C₃/r³ long-range tail of the potential like LEVEL does?).
• Strangely, increasing the Grid value usually made the numerical solution worse (leading to more convergence problems and/or fewer vibrational energies found). For example, the above-linked output file shows that I used Grid = 470 and when using 500 I got marginally more levels, but when using Grid = 510 the program found fewer levels again.
• Decreasing Step surprisingly does reduce the number of levels (marginally) but increasing it would make the cost outrageously long for the standards of a rovibrational level solver (the runtime is already extremely long with the Step value in the above calculation).

Is it possible with VibRot to get a number of vibrational levels closer to what the above-linked LEVEL output gave me (which was 98 vibrational levels including $v=0$)?

Answer 4

What are the seminal papers for ERIs?

• 1950 Boys

Based on the Rys quadrature:

• 1976,1983 Dupuis, Rys, and King x2 (Rys quadrature)
• 1978 Pople and Hehre (uses Rys quadrature, good for highly contracted s and p functions)
• 1986, 1988 Obara and Saika x 3 (reduce FLOPS by recurrence relation)
• 1988, 1989, 1990, 1991 Gill, Head-Gordon and Pople x 4 (reduce FLOPS, contraction scheme, "horizontal transfer equations")
• 1991, 1993 Lindh, Ryu, and Liu x2 (OS and GHP are less efficient than DRK when there's not much contraction)
• 1991,1996 Ishida x 2 (ACE has "attractive" FLOPS)
• 2001 Dupuis and Marquez (combination of the above to optimize FLOPS) [implementable! but notation isn't familiar to all]
• 2004 Valeev libint
• 2014 Sun (libcint, uses DRK algorithm)
• 2014 Schwenke
• 2015 Shizgal
• 2016 King

Based on the McMurchie-Davidson scheme

• 1978 McMurchie Davidson (angular momentum transfer)
• 2022 Neese (SHARK) [implementable!]

Independent/largely-unnoticed

• 2018 Milovanovic

Parallelization and hardware-related optimizations

• 1994 Blelloch et al.
• 2007 Yasuda
• 2010 Asadschev et al.
• 2011 Luehr et al.
• 2012 Asadchev and Gordon
• 2013 Titove et al.

Reviews

• 1994 M. W. Gill,Adv. Quantum Chem.1994,25, 141.
• 2002 R. Lindh, in Encylopedia Comp Chem (Ed: P. von Rague-Schleyer), JohnWiley & Sons, Chichester
• 2018 Samu (PhD thesis, notation still complicated) PDF/ertekezes.pdf

Answer 5

Who was JK Cashion?

Searching "JK Cashion" in Google seems to result in nothing but papers. I cannot even find his or her first name!

However JK Cashion remains legendary in the field of diatomic spectroscopy due to his seminal papers in the 1960s on numerically solving the Schroedinger equation to get the energies and wavefuncions of vibrational and rotational states of molecules with potentials that are more realistic than the very simple ones for which the problem is solvable analytically.

JK Cashion was mentioned in the Nobel Prize lecture of John Polyani (see this PDF). Polanyi's official website also mentions JK Cashion, saying that Cashion was Polanyi's graduate student and that their two-author 1958 paper introduced a method that "promises to provide for the first time information concerning the distribution of vibrational and possibly rotational energy among the products of a three-center reaction" which was then realized about a decade later, after the "work of a substantial group of graduate students". Interestingly, Google Scholar provides a citation to Cashion's PhD thesis from 5 years earlier (and 2 years before Polanyi joined that university, see this PDF). It would appear from his two-author paper with Donald James LeRoy (father of Robert J. Le Roy), published 1 year after the 1953 PhD thesis (see this PDF) that he did a first PhD with DJ LeRoy followed by becoming a grad student of Polanyi. Maybe Google Scholar has erroneously called his Bachelors or Masters level thesis of 1953 a PhD thesis, but in any case I was eventually able to find his 1960 PhD thesis (which also finally gives his full name: John Kenneth Cashion!). Searching "John Kenneth Cashion" basically gives 1 new hit, which is a black-and-white picture of a professor emeritus at University of Wisconsin-Parkside but not much more.

After a 1961 paper with Polanyi, Cashion seems to have joined the "Aeronical Research Laboratories" where he completed a solo-author technical report in 1962, followed by a move to what is now called the Lawrence Berkeley National Lab where in 1963 he published papers by himself and with Nobel Prize winner Dudley Herschbach's then-PhD-student Ricahrd Zare and then a papers with Herschbach himself in 1964 including an empirical potential energy surface for H3 which has been described as the Cashion-Herschbach potential in the title of a paper as recently as 2000.

When Herschbach (and Zare) moved to Harvard, it looks like Cashion came along too. Cashion and Herschbach have a 1964 paper together with Berkeley affiliation and also a 1964 paper together with Harvard affiliation, but then he stopped publishing with them and instead did an intermolecular potential for (H₂)₂ published in 1965 with Nobel Prize winner Van Vleck's Harvard PhD student Roy Gordon, followed by three solo-author papers with Harvard affiliation in 1966 (here, here, and here), seemingly nothing in 1967, then a solo-author 1968 paper with IBM Waston Lab affiliation.

I'm not able to find any 1969 papers by Cashion, but at IBM he published 4 papers with Dean E. Eastman and others from 1970-1971 including three in PRL! But Cashion's career publishing in scientific journals seems to have ended in 1971, Dean Eastman continued publishing until this seemingly final 1988 Science paper and went on to become a very successful executive at IBM and then the Arognne National lab (part of a pattern of success with Cashion's other collaborators: Nobel Prize winners Polanyi and Herschbach, extremely successful professors Zare and Gordon, and high-ranking executives LeRoy and Eastman), but Cashion seems to have fallen off the visible record after 1971.

---

After referring to some of the original papers (including Cashion's) on solving the Schroedinger equation to get rovibrational energies and wavefunctions of diatomic molecules, and not being able to find anything (not even a first name!) in a Google search, I felt like asking "who was JK Cashion?" here. In my requisite research prior to posting, I was able to trace his life from Toronto to Berkeley to Harvard to IBM, who he worked with, and when he published (almost exclusively through Google Scholar and looking at actual papers), but I wonder if someone who knows better might be able to write a biography of him in an answer here?

Answer 6

Why does our PySCF-MRCC interface work for B++ but not for Li?

Background

My colleague and I have built an interface between PySCF and MRCC which seems to work all the time for RHF references, but only sometimes for UHF references. We use PySCF to calculate the integrals and print them in FCIDUMP format, then we convert the FCIDUMP file into MRCC's fort.55 format.

The format of the FCIDUMP follows the UHF standard designed by Knowles and Handy, in which after the header part, separate 4-index and 2-index integrals are written for alpha-alpha, beta-beta, and alph-beta orbitals, with a line of zeroes separating them. Unfortunately, despite this format being used in MOLPRO, and being compatible with the fort.55 file in some major quantum chemistry codes like MRCC or CFOUR, it is not natively supported in PySCF.

The Problem

We considered two open-shell systems, Li and B(+2). Both systems have 3 electrons occupying s-type orbitals. For B(+2), the coupled cluster energies obtained from my PySCF-MRCC interface are identical to the ones obtained using MRCC-MRCC (which is how we will label "standalone MRCC" in which the integrals, SCF and AO->MO transformation are calculated using MRCC, and the post-SCF calculations are also done by MRCC).

B[+2] (6-31G) (symmetry is on)

mrcc-mrcc
CCSD    -23.372809509970
CCSD(T) -23.372810050030
CCSDT   -23.372810049178

pyscf-mrcc
CCSD    -23.372809509970
CCSD(T) -23.372810050030
CCSDT   -23.372810049178

But for Li, this is what I get:

Li (6-31G) (symmetry is on)

mrcc-mrcc
UHF.    <please add this to show that it's the same>
CCSD    -7.431553874113
CCSD(T) -7.431554248616
CCSDT   -7.431554228106

pysc-mrcc
UHF     <please add this to show that it's the same>
CCSD    -7.431274408997
CCSD(T) -7.431274447587
CCSDT   -7.431274443531

Upon some further debugging, I am finding that the generated fort.55 in the case B(+2) is similar to the one generated in MRCC in terms of having the same ORBSYM and the same indices for the integrals. However, for the case of Li, the generated fort.55 is having similar ORBSYM with that generated with MRCC but the indices and the number of integrals are not the same. You can find full files for the inputs and outputs generated with PySCF and MRCC in the following Input/Output Files.

I will point out that the energies actually do match if I turn off symmetry, but then the RAM and CPU time are significantly greater, which means that we won't be able to do very-high-order coupled cluster calculations with big basis sets.

Why is the code described below working for B++ but not for Li?

The code

Suppose that we have defined a Mol object for an open-shell system, followed by running a UHF calculation. In the next step, we convert the 2e integrals using 4-fold permutation symmetry for alpha-alpha, beta-beta, and alpha-beta orbitals which is saved in corresponding eri_aaaa, eri_bbbb, and eri_aabb, respectively. Moreover, the transformed 1e core Hamiltonian for alpha and beta are also calculated and stored in h_aa and h_bb, respectively.

orbs = mf.mo_coeff
nmo = orbs[0].shape[0]

eri_aaaa = pyscf.ao2mo.restore(4,pyscf.ao2mo.incore.general(mf._eri, (orbs[0],orbs[0],orbs[0],orbs[0]), compact=False),nmo)
eri_bbbb = pyscf.ao2mo.restore(4,pyscf.ao2mo.incore.general(mf._eri, (orbs[1],orbs[1],orbs[1],orbs[1]), compact=False),nmo)
eri_aabb = pyscf.ao2mo.restore(4,pyscf.ao2mo.incore.general(mf._eri, (orbs[0],orbs[0],orbs[1],orbs[1]), compact=False),nmo)

h_core = mf.get_hcore(mol)
h_aa = reduce(numpy.dot, (orbs[0].T, h_core, orbs[0]))
h_bb = reduce(numpy.dot, (orbs[1].T, h_core, orbs[1]))
nuc = mol.energy_nuc()

In the next section of the code, we try to handle the symmetry of the orbitals in a way that will match the format of fort.55. We also define the range of values for integral indices:

if mol.symmetry:
        groupname = mol.groupname
        if groupname in ('SO3', 'Dooh'):
            logger.info(mol, 'Lower symmetry from %s to D2h', groupname)
            raise RuntimeError('Lower symmetry from %s to D2h' % groupname)
        elif groupname == 'Coov':
            logger.info(mol, 'Lower symmetry from Coov to C2v')
            raise RuntimeError('''Lower symmetry from Coov to C2v''')
orbsym = pyscf.symm.label_orb_symm(mol,mol.irrep_name,mol.symm_orb,orbs[0])
orbsym = numpy.array(orbsym)
orbsym = [param.IRREP_ID_TABLE[groupname][i]+1 for i in orbsym]
a_inds = [i+1 for i in range(orbs[0].shape[0])]
b_inds = [i+1 for i in range(orbs[1].shape[1])]
nelec = mol.nelec
tol=1e-18

The last part is considered the main part of the script, which is supposed to handle the header of the generated file in a way that will match the fort.55 format. It is supposed to deal with the 4-index, and 2-index integrals in which we are assuming that we are using a 4-fold permutation symmetry. The sections after '4-fold symmetry' are in the sequence of alpha-alpha, beta-beta, alpha-beta 2e integrals, followed by alpha, beta 1e integrals, and finally the nuclear repulsion energy.

with open('fort.55', 'w') as fout:
        if not isinstance(nelec, (int, numpy.number)):
            ms    = abs(nelec[0] - nelec[1])
            nelec =     nelec[0]  + nelec[1]
        else: ms=0
        fout.write(f"{nmo:1d} {nelec:1d}\n")
        if orbsym is not None and len(orbsym) > 0:
            fout.write(f"{' '.join([str(x) for x in orbsym])}\n")
        else:
            fout.write(f"{' 1' * nmo}\n")
        fout.write(' 150000\n')
        output_format = float_format + ' %5d %5d %5d %5d\n'
        #4-fold symmetry
        kl = 0
        for l in range(nmo):
            for k in range(0, l+1):
                ij = 0
                for i in range(0, nmo):
                    for j in range(0, i+1):
                        if i >= k:
                            if abs(eri_aaaa[ij,kl]) > tol:
                                fout.write(output_format % (eri_aaaa[ij,kl], a_inds[i], a_inds[j], a_inds[k], a_inds[l]))
                        ij += 1
                kl += 1
        fout.write(' 0.00000000000000000000E+00' + '     0     0     0     0\n')
        kl = 0
        for l in range(nmo):
            for k in range(0, l+1):
                ij = 0
                for i in range(0, nmo):
                    for j in range(0, i+1):
                        if i >= k:
                            if abs(eri_bbbb[ij,kl]) > tol:
                                fout.write(output_format % (eri_bbbb[ij,kl], b_inds[i], b_inds[j], b_inds[k], b_inds[l]))
                        ij += 1
                kl += 1
        fout.write(' 0.00000000000000000000E+00' + '     0     0     0     0\n')
        ij = 0
        for j in range(nmo):
            for i in range(0, j+1):
                kl = 0
                for k in range(nmo):
                    for l in range(0, k+1):
                        if abs(eri_aabb[ij,kl]) > tol:
                            fout.write(output_format % (eri_aabb[ij,kl], a_inds[i], a_inds[j], b_inds[k], b_inds[l]))
                        kl += 1
                ij +=1

        fout.write(' 0.00000000000000000000E+00' + '     0     0     0     0\n')
        h_aa = h_aa.reshape(nmo,nmo)
        h_bb = h_bb.reshape(nmo,nmo)
        output_format = float_format + ' %5d %5d     0     0\n'
        for i in range(nmo):
            for j in range(nmo):
                fout.write(output_format % (h_aa[i,j], a_inds[i], a_inds[j]))
        fout.write(' 0.00000000000000000000E+00' + '     0     0     0     0\n')
        for i in range(nmo):
            for j in range(nmo):
                fout.write(output_format % (h_bb[i,j], b_inds[i], b_inds[j]))
        fout.write(' 0.00000000000000000000E+00' + '     0     0     0     0\n')
        output_format = float_format + '     0     0     0     0\n'
        fout.write(output_format % nuc)

Answer 7

Why do the PBE pseudopotentials for Zn and S from pseudo-dojo give such bad lattice constant (~10%) in siesta but it's reasonable for the PBEsol(~1%) or even LDA pseudopotentials(~3%)?

How to generate my own pseudopotential in psml format for siesta 5.2?