sharc-md/sharc

Issues with overlap calculation in Bagel Interface

MartinPeschel opened this issue · 15 comments

Hi!

I have encountered an issue when running initial conditions with the Bagel interface.
The Bagel calculation runs smoothly. The calculation of the initial overlaps crashes.

QM.out:

>>>>>>>>>>>>> Starting the BAGEL job execution
copying from file  /home.io/map/sharc/sharc_3.0.1/sharc/tests_rdv/traj_bagel/ICOND_00001/SAVE/archive.1.old.archive  to file  /scratch/map/81820/ICOND_00001/master_1/archive.archive

Running with OpenMP (1 cores)
START:  /scratch/map/81820/ICOND_00001/master_1         2023-10-31 15:56:25.713577      "/opt/adm/bagel/rev-210120/bin/BAGEL     "
FINISH: /scratch/map/81820/ICOND_00001/master_1         2023-10-31 15:58:45.997774      Runtime: 0:02:20.284197 Error Code: 0
>>>>>>>>>>>>> Saving files
/home.io/map...01/SAVE/orbitals.molden.1
  .. transforming to cartesian AO basis
/home.io/map...el/ICOND_00001/SAVE/mos.1
/home.io/map...l/ICOND_00001/SAVE/dets.1
Saving Runtime: 0:00:00.317126
>>>>>>>>>>>>> Starting the WFOVERLAP job execution
START:  /scratch/map/81820/ICOND_00001/WFOVL_1_1        2023-10-31 15:58:48.329854      "/home.io/map...ap.x -m 1000 -f wfovl.inp"
FINISH: /scratch/map/81820/ICOND_00001/WFOVL_1_1        2023-10-31 15:58:48.428781      Runtime: 0:00:00.098927 Error Code: 1
Error Codes:
        WFOVL_1_1       1
Some subprocesses did not finish successfully!

Turning debug mode on, it seems that the number of orbitals in the mo-files and the aoovl file
that are passed to wfoverlap are incompatible.

wfovl.err:
read_col_moheader: MOcoef file states different number of atomic orbitals than AOoverlap file. naotmp = 120, nao = 114 1

The difference seems to be the number of frozen core orbitals.
No frozen core orbitals in the MO-files, 6 frozen core orbitals in the aoovl file.
Any idea what might be going wrong?

Hi Martin,
frozen core orbitals should not affect the number of AOs, only the number of MOs. My suspicion would be that there might be a mix up between Cartesian versus spherical harmonics. What basis set are you using? Can you replace it by another one that is about the same quality (e.g., if you used 6-31G*, try with def2-SVP or cc-pvdz instead) and try whether it works?

Best,
Sebastian

Hi Sebastian,
I replaced "cc-pvdz" by "6-31g". With "6-31g" it works correctly.
Best, Martin

But be aware that 6-31g does not have polarization functions (or does it in Bagel?). If your calculations are just for testing, that is not a problem. But for a serious simulation with the goal of a publication, I strongly recommend against using 6-31g.

Also, one of the examples in $SHARC/../examples uses cc-pvdz, so I would expect that it works with cc-pvdz. Could you send me an input to reproduce your error (either here or via email)?

I appended the ICOND_00000 and ICOND_00001 folders.
ICOND_00000.zip
ICOND_00001.zip
You will only have to change the basis and df_basis parts in BAGEL.template. I checked that my version of the cc-pvdz
basis-file is identical with the one provided in the examples folder.

Hi Martin,
I checked your inputs and can reproduce your error. The reason is some code in make_mos_from_Molden that was not checked and adapted when converting from PyQuante to PySCF. This code converts the MO coefficients read from the Molden file from mixed/spherical representation to Cartesian (because PyQuante could only work in Cartesian AOs). Thus, in your example with cc-pvdz, Bagel uses 114 spherical AOs, but the interface was producing MOs based on 120 basis functions.

The error disappears if you comment out lines 2718 to 2771 in SHARC_BAGEL.py. However, after this change the obtained wave function overlaps do not look correct for me and I will investigate this further. Here, it would be good to have the archive file with your initial orbitals, as the orbitals that I see in the Molden files from Bagel do not look reasonable.
Best,
Sebastian

Hi Sebastian,
Thanks you! That means, PyQuante is no longer needed for BAGEL calculations with the new interface?

I just checked the molden files in the SAVE folder, and they seem reasonable... the guess of bagel would have to be pretty poor to not get the pi space of benzene correct. I attached my archive and molden files again, so you can check...

In my experience, moving archive files from one cluster to another (and one bagel version to another) does not work.
Therefore I have attached a Bagel input to generate the orbitals.

orbitals.zip
archive.zip
BAGEL.json

Edit: Updated the orbitals file, there was an issue with the .zip...

Best, Martin

Dear Sebastian and Martin,

I'm encountering the same issue in my calculations. To note, I'm using the 6-31++G(d) basis set.

Was this issue ever fully resolved?

Thanks,

Andrew

To follow on from my previous comment.

Hardcoding the following into the writeBAGELinput(QMin) function, just after it prints the dkh flag, seems to have done the trick (i.e. forcing the BAGEL calculation to use cartesian basis functions).

string += '    "cartesian": "true",\n'

Is it true to assume that wfoverlap cannot use spherical basis functions or is it something in the conversion that SHARC performs?

Thanks,

Andrew

Hi all,
thanks for bringing up the topic again. Indeed, we were never able to resolve this issue. Enforcing Cartesian basis functions does seem to help to let the job run through, but the resulting overlap matrix does not look correct. I will have another look, maybe I can figure it out.

Best,
Sebastian

Coming closer... Apparently, the MO coefficients in the Molden file are garbage if Bagel writes the Molden after doing any gradient/NAC calculation. For a single-point without grad/nac, the orbitals look good. Now I can try to make the MO coefficients and the AO overlaps consistent...

Hmm....that's rather strange, thanks for letting me know of this.

So just to be clear this is an issue for all sharc-bagel wfoverlap calculations and isn't related to the discrepancy between basis functions?

Yes, this seems to be a general issue with all wfoverlap calculations in SHARC_BAGEL since SHARC 3.0, when we ditched PyQuante and ported to PySCF. This was apparently not fully tested. Sorry for the inconvenience.

Update: I now reorder the MO coefficients so that the shell order matches between the Molden file and the PySCF-derived AO overlaps (i.e., xx, yy, zz, xy, xz, yz VS xx, xy, xz, yy, yz, zz). But the overlaps still do not look good. There seems to be something else, possibly sign conventions of the basis functions or basis function normalization.

Ok, thanks for the update. In the meantime I'll revert back to sharc-2.1.

Thanks,

Andrew

It might be possible to use SHARC 3 together with the SHARC_BAGEL.py from SHARC 2.1, if you want to use other features of SHARC 3.0 (like options in sharc.x or some analysis tools). You just need to be careful with mixing Python 2 and Python 3.

In the sharc3preview branch, I now pushed a fix. It has the following changes:

  • Molden files are saved before any force calculations. This makes it slightly slower, because one has to calculate the energy first, then partially redo the energy for the forces. This might jeopardize phase-consistency for NACs, but I did not test this. As long as you do not use NACs, this is fine.
  • Cartesian basis functions are always enforced
  • AO overlap integrals are properly rescaled (PySCF does not give AO overlaps where all basis functions are normalized to 1).
  • MO coefficients read from Molden are resorted to match the PySCF order. This works up to f functions.

For your test example (benzene, CAS(6,6)/cc-pvdz) it works. You can check this for yourself by doing the same SHARC singlepoint calculation twice, (once with "init", once with "overlap" in QM.in), and you should obtain a unit overlap matrix.

Can you please test it and confirm that it works. Then I can merge into the main branch