Implemented the ECP part and added several missing factors in the TREXIO interface.

[kousuke-nakano]

@matthew-hennefarth

Sorry for my late response. Today, I have implemented lines to write ECP information to a TREXIO file. Please have a look at my pull request.

The ECP part works. However, a problem I found is that the cc-pVXZ basis sets accompanied with ccECP, which are often used for QMC calculations, did not pass the following assertion line in trexio.py.

  assert all(mol._bas[:,gto.NCTR_OF] == 1)

Would you give me a clue for solving this problem?

#60

[matthew-hennefarth]

The ECP part works. However, a problem I found is that the cc-pVXZ basis sets accompanied with ccECP, which are often used for QMC calculations, did not pass the following assertion line in trexio.py.

  assert all(mol._bas[:,gto.NCTR_OF] == 1)

It seems that mol._bas is the internal representation for the libcint library, and NCTR_OF is a constant being 3 (which I believe represents number of cartesians; ie x, y, z). Unfortunately, I do not know what the internal representation of _bas is supposed to be.

Perhaps, since @sunqm wrote the assert statement and knows the details of libcint, he can provide you with more guidance.

[kousuke-nakano]

@matthew-hennefarth, thank you very much. I will wait for @sunqm's reply.

[kousuke-nakano]

@matthew-hennefarth

I found that mol._bas[:,gto.NCTR_OF] refers to the numbers of generalized contractions (i.e., one set of exponents and multiple sets of coefficients, like molpro and molcas). See this issue. Do you know how to switch it off with some option when we build a mol?

[matthew-hennefarth]

I found that mol._bas[:,gto.NCTR_OF] refers to the numbers of generalized contractions (i.e., one set of exponents and multiple sets of coefficients, like molpro and molcas). See this issue. Do you know how to switch it off with some option when we build a mol?

I am not sure how to turn it off, sorry. From the discussion you linked it seems like the best bet would be to just manually expand the contraction set.

[kousuke-nakano]

Dear @matthew-hennefarth, Sure. Thank you! Anyway, this is nothing with the ECP part.

Dear @matthew-hennefarth and @sunqm, I have completed implementing a function to write the ECP part to a TREXIO file. The read function is my ToDO work. In addition, I have added several missing factors, such as AO_normalization (#80). This pull request is ready for merge. I would appreciate it if you could review and merge it.

[kousuke-nakano]

Sorry @sunqm, I committed new changes (6949296) by mistake. I have reverted it (db9cde9).

[kousuke-nakano]

@sunqm, I have resolved a trivial conflict.

[kousuke-nakano]

Dear @sunqm, Thank you for running the workflows. One of the workflows failed at the part unrelated to this pull request.

[q-posev]

@kousuke-nakano are you done with this PR?

I left one comment about the UHF case, let me know if you prefer to fix it here or if you prefer me to do it in another PR. It should not be too complicated, just a question of extracting the number of alpha and beta electrons from the UHF object.

[kousuke-nakano]

Dear @q-posev, @sunqm has already implemented both RHF and UHF in #60.

[q-posev]

@kousuke-nakano yes, but on line 63 in trexio.py in this PR you added the following line:

    electron_up_num, electron_dn_num = mol.nelec

which is not true for UHF.

The original code of @sunqm was not producing the number of electrons.

[q-posev]

The electron number writing can be probably moved to _scf_to_trexio function where we have an if/else block depending on the output of isinstance(mf, scf.uhf.UHF) (I mean this line).

[kousuke-nakano]

Dear @q-posev,

@kousuke-nakano yes, but on line 63 in trexio.py in this PR you added the following line:

    electron_up_num, electron_dn_num = mol.nelec

which is not true for UHF.

The following code gives the correct numbers of up and down electrons.

import pyscf

mol_S_0 = pyscf.M(atom='H 0 0 0; F 0 0 1', basis='6-31g*', cart=True)
electron_up_num, electron_dn_num = mol_S_0.nelec
print(electron_up_num, electron_dn_num) # -> 5 5

mol_S_2 = pyscf.M(atom='H 0 0 0; F 0 0 1', basis='6-31g*', cart=True)
mol_S_2.spin = 2
electron_up_num, electron_dn_num = mol_S_2.nelec
print(electron_up_num, electron_dn_num) # -> 6 4

[q-posev]

@kousuke-nakano you are totally right, my bad! Please ignore my comments above.