mrchem icon indicating copy to clipboard operation
mrchem copied to clipboard

Rewrite read-in of SAD guess

Open robertodr opened this issue 4 years ago • 14 comments

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 file O.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 the AOBasis class to parse the basis as copy-pasted from the BSE.
  • Rework project_density such that it reads the density 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)

robertodr avatar Nov 11 '20 08:11 robertodr