Impossible to import basis set for second-row elements

[NastaMauger]

Hello,

This minimal example failed to be imported into QP2 with the following error:

Nuclei		...	 OK
Electrons	...	 5 5
Basis		...	 Gaussian
AOS		...	 Traceback (most recent call last):
  File "/opt/qp2_dev_stable/scripts/qp_import_trexio.py", line 576, in <module>
    write_ezfio(trexio_filename, EZFIO_FILE)
  File "/opt/qp2_dev_stable/scripts/qp_import_trexio.py", line 317, in write_ezfio
    assert (len(coefficient) == ao_num)

I realized that the issue is related to the basis set of the carbon atoms, not hydrogen. If I use H or He it works
I used the cc-pVDZ basis set in an uncontracted form from BSE. I am not sure if this is related to issue #42, so I opened a new issue. Feel free to close it if you think it might be redundant

from pyscf import gto, scf, symm, mcscf, fci, ao2mo, lib                                                                                                                                                                                                                                                                   
import pyscf.mcscf.chkfile as chkfile
import numpy as np
import os

# Define the molecule
# Trexio needs to have the uncontracted general basis from BSE
# Basis set: cc-pvdz
mol = gto.M(
   atom = ''' 
     C 0.0000   0.0000  0.0000
     H 0.6266   0.6266  0.6266
     H 0.6266  -0.6266 -0.6266
     H -0.6266  0.6266 -0.6266
     H -0.6266 -0.6266  0.6266
   ''',
   charge = 0,
   spin = 0,
   symmetry = True,
#    symmetry_subgroup = 'D2h',
   verbose = 4,
   basis ={
     'H': gto.parse('''
           H    S
                 1.301000E+01           1.968500E-02
                 1.962000E+00           1.379770E-01
                 4.446000E-01           4.781480E-01
                 1.220000E-01           5.012400E-01
           H    S
                 1.220000E-01           1.000000E+00
           H    P
                 7.270000E-01           1.0000000
     '''),
     'C': gto.parse('''
C    S
     6.665000E+03           6.920000E-04
     1.000000E+03           5.329000E-03
     2.280000E+02           2.707700E-02
     6.471000E+01           1.017180E-01
     2.106000E+01           2.747400E-01
     7.495000E+00           4.485640E-01
     2.797000E+00           2.850740E-01
     5.215000E-01           1.520400E-02
     1.596000E-01          -3.191000E-03
C    S
     6.665000E+03          -1.460000E-04
     1.000000E+03          -1.154000E-03
     2.280000E+02          -5.725000E-03
     6.471000E+01          -2.331200E-02
     2.106000E+01          -6.395500E-02
     7.495000E+00          -1.499810E-01
     2.797000E+00          -1.272620E-01
     5.215000E-01           5.445290E-01
     1.596000E-01           5.804960E-01
C    S
     1.596000E-01           1.000000E+00
C    P
     9.439000E+00           3.810900E-02
     2.002000E+00           2.094800E-01
     5.456000E-01           5.085570E-01
     1.517000E-01           4.688420E-01
C    P
     1.517000E-01           1.000000E+00
C    D
     5.500000E-01           1.0000000
      ''')},
)

# Perform SCF calculation
print("------------------------------------------------")
print("----------------------RHF-----------------------")
print("------------------------------------------------")
rhf = scf.RHF(mol)
rhf.verbose = 4
rhf.chkfile = 'rhf.chk'
rhf.kernel()

filename = 'rhf.hdf5'
# Convert PySCF checkpoint to TREXIO
os.system('trexio convert-from -t pyscf -i rhf.chk -b hdf5 rhf.hdf5')

[NastaMauger]

I encountered this issue as soon as I used a basis set for a second-row element. Not only C

[q-posev]

@NastaMauger can you try to use pyscf-forge for producing TREXIO instead of the converter from this repo? Does the error persists there?

[NastaMauger]

@q-posev The error is slighly different now with the lastest pyscf-forge

Nuclei		...	 OK
Electrons	...	 2 2
Basis		...	 Traceback (most recent call last):
  File "/opt/qp2_dev_stable/scripts/qp_import_trexio.py", line 576, in <module>
    write_ezfio(trexio_filename, EZFIO_FILE)
  File "/opt/qp2_dev_stable/scripts/qp_import_trexio.py", line 154, in write_ezfio
    ao_shell    = trexio.read_ao_shell(trexio_file)
  File "/opt/trexio/python/trexio.py", line 9474, in read_ao_shell
    raise Error(rc)
trexio.Error: Dataset does not exist in the file

which is related to our previous discussions. This happens no matter the molecule (first or second row).

Now if I use chk file and convert it using trexio-tools, I have the AO error only for second row elements.

Best

[q-posev]

Are you sure that you are using the latest commit of pyscf-forge? Our friend and colleague @kousuke-nakano implemented the AO shell writing in MR #68, which appears on the master branch as commit 9564956cb7443394f0aacb6a51db03b5da2f69f2.

Also, I just noticed that Qiming Sun has modified the trexio interface to work for generally contracted basis sets too. It's in MR #84, you may want to try it out.

[kousuke-nakano]

Dear @NastaMauger, thank you for the report. It seems that you are using the pyscf-to-trexio converter of this repo.

# Convert PySCF checkpoint to TREXIO
os.system('trexio convert-from -t pyscf -i rhf.chk -b hdf5 rhf.hdf5')

It is no longer maintained because we are now implementing the converter directly inside pyscf. Would you try the converter implemented in pyscf-forge?

Dear @q-posev, Shall we remove the pyscf-to-trexio converter of this repo, though PBC implementations in pyscf-forge are still in progress?

[q-posev]

Dear @kousuke-nakano ,
I have nothing against that if that's OK with you. Tagging @scemama and @NastaMauger to confirm with them too as they are the main users of the pyscf-to-trexio converter in this repo.

[NastaMauger]

Hello

I am confirming that I am using lastest commit from pyscf-forge.
Pyscf-forge leads to this error doing qp_import_trexio.py data.hdf5 :

Nuclei		...	 OK
Electrons	...	 2 2
Basis		...	 Traceback (most recent call last):
  File "/opt/qp2_dev_stable/scripts/qp_import_trexio.py", line 576, in <module>
    write_ezfio(trexio_filename, EZFIO_FILE)
  File "/opt/qp2_dev_stable/scripts/qp_import_trexio.py", line 154, in write_ezfio
    ao_shell    = trexio.read_ao_shell(trexio_file)
  File "/opt/trexio/python/trexio.py", line 9474, in read_ao_shell
    raise Error(rc)
trexio.Error: Dataset does not exist in the file

And this error comes no matter the elemet (first/second row)

Here a minimal example which reproduces the problem:

# Basis set: cc-pvdz
mol = gto.M(
    atom = ''' 
       He 0.0000   0.0000  0.0000
       He 0.0000   0.0000  3.5000
    ''',
    charge = 0,
    spin = 0,    
    symmetry = True, 
    symmetry_subgroup = 'D2h',
    verbose = 4,
    max_memory = 700000,
    basis ={'He': gto.parse('''
            #BASIS SET: (10s,5p,2d) -> [4s,3p,2d]
        He    S
              3.836000E+01           2.380900E-02
              5.770000E+00           1.548910E-01
              1.240000E+00           4.699870E-01
              2.976000E-01           5.130270E-01
        He    S
              2.976000E-01           1.000000E+00
        He    P
              1.275000E+00           1.0000000
              ''')},
)
rhf = scf.RHF(mol)
rhf.verbose = 4 
rhf.kernel()

trexio_file = 'data.hdf5'
to_trexio(rhf, trexio_file)
eri = mol.intor('int2e')
write_eri(eri, trexio_file)

Best

[kousuke-nakano]

Dear @NastaMauger

I cannot reproduce your result with the latest pyscf-forge repo.

import trexio
from pyscf import gto, scf
from pyscf.tools import trexio as pyscf_trexio

mol = gto.M(
atom = '''
He 0.0000 0.0000 0.0000
He 0.0000 0.0000 3.5000
''',
charge = 0,
spin = 0,
symmetry = True,
symmetry_subgroup = 'D2h',
verbose = 4,
max_memory = 700000,
basis ={'He': gto.parse('''
#BASIS SET: (10s,5p,2d) -> [4s,3p,2d]
He S
3.836000E+01 2.380900E-02
5.770000E+00 1.548910E-01
1.240000E+00 4.699870E-01
2.976000E-01 5.130270E-01
He S
2.976000E-01 1.000000E+00
He P
1.275000E+00 1.0000000
''')},
)
rhf = scf.RHF(mol)
rhf.verbose = 4
rhf.kernel()
eri = mol.intor('int2e')

trexio_file='data_rhf.hdf5'
pyscf_trexio.to_trexio(rhf, trexio_file)
pyscf_trexio.write_eri(eri, trexio_file)

file_r = trexio.File(trexio_file, mode="r", back_end=trexio.TREXIO_HDF5)
ao_shell = trexio.read_ao_shell(file_r)
print(ao_shell)
[0 1 2 2 2 3 4 5 5 5]

This works.

I am sure I have implemented ao_shell in this pull request.

Would you let us know the detail of your error?

[NastaMauger]

Hello @kousuke-nakano

In your last commit, you also added the PySCF version to the metadata of TREXIO file.
Here is my data.hdf5.zip. The PySCF version in the metadata confirms that I am using the correct PR. However, this HDF5 file leads to the error I previously mentioned with QP2 (qp_import_trexio.py data.hdf5).

Best

[NastaMauger]

The data are in TREXIO. The issue might then arise from the QP2 import
@scemama might know what the issue is here.

[kousuke-nakano]

Dear @NastaMauger, thank you! Let's wait for @q-posev and @scemama replies. I have no idea about this converter, sorry.

[q-posev]

@NastaMauger the PySCF version has nothing to do with the version of pyscf-forge. Please install the latest commit from the master branch in your software environment (ideally in an isolated Python env like venv or conda where you can control the version of both pyscf and pyscf-forge). You can install the latest pyscf-forge with: pip install git+https://github.com/pyscf/pyscf-forge.

[q-posev]

@kousuke-nakano thanks for testing!

[q-posev]

@NastaMauger sorry, I missed your comment. Indeed, you are right, the ao_shell dataset is present in the HDF5 file. I don't understand how it's possible that the Python interface reports it as missing. Are you sure that you provide the same file to the QP2 import function?