mrchem
mrchem copied to clipboard
Rewrite read-in of SAD guess
Currently the SAD guess is read from 2 files:
- The basis set file in DALTON format (I believe) in
O.bas
- The atomic density in
O.dens
The two can be unified into a single JSON fileO.json
that contains both the basis and the corresponding atomic density. The basis can be taken directly from the Basis Set Exchange (BSE), see here for 3-21G of oxygen. The density can be computed with any code. The JSON would have to keys:
{
"basis": { ... }, # <--- copy-paste from the BSE repository,
"density": [ ... ] # <--- 1D array with the atomic density.
}
In the C++ code one needs to:
- Add a method
from_json
to theAOBasis
class to parse the basis as copy-pasted from the BSE. - Rework
project_density
such that it reads thedensity
field of the JSON file.
The advantages of this approach are that:
- Easier for developers and it's in a standard, human-readable format
- We are sure that the basis set is correct, since it comes from the BSE
- One can build its own SAD guesses on top of the ones shipped with MRChem, for example, to use a larger atomic basis. I am not sure how useful this is in general though.
On top of the C++ changes there should be machinery to generate the SAD guesses in JSON format. This can be nicely achieved with Psi4 and the BSE Python library. Note that I haven't tested the following script:
import json
import numpy as np
import basis_set_exchange as bse
import psi4
mol = psi4.geometry("""
O 0.0 0.0 0.0
symmetry c1
noreorient
no_com
""")
psi4.core.set_active_molecule(mol)
# set options
options = {
'SCF_TYPE':'PK',
'E_CONVERGENCE':1e-12,
'D_CONVERGENCE':1e-12
}
# build custom basis from BSE
def define_custom_basis(mol, role):
basstrings = {}
mol.set_basis_all_atoms("3-21G", role=role)
basstrings['3-21g'] = bse.get_basis("3-21G", elements=[8], fmt="psi4")
return basstrings
psi4.qcdb.libmintsbasisset.basishorde['3-21G_BSE'] = define_custom_basis
psi4.core.set_global_option("BASIS", "3-21G_BSE")
psi4.set_options(options)
e, wfn = psi4.energy('SCF', return_wfn=True)
# density as NumPy array
density = wfn.Da().to_array()
output = Path("O.json")
with output.open("w") as f:
json.dump({"basis": bse.get_basis("3-21G", elements=[8], fmt="json"), "density": density.tolist()} , f)