MolSSI
How to store gradients and hessians?
I'm wondering if it's possible to store gradients and hessians in a qcschema. The docs do mention this
Returned Result The “primary” return of a given computation. For energy, gradient, and Hessian quantities these are either single numbers or a array representing the derivative quantity.
I couldn't find an example of a json file like that though. Would the array be flattened or nested? What is the order? eg for gradients [gx1, gy1, gz1, ... ] or [gx1, gx2, gx3, ...]
Also, it seems the gradient or hessian would be a "return_result", not a "property" - what happens if there's more than one, for example? It seems strange that energies are both properties and return results, but gradients are only return results.
Is there a way to store this type of values as properties, eg "scf_gradient" similar to "scf_total_energy"? What is the recommended best practice?
Here's an example of a gradient input and result. QCElemental has a working implementation of QCSchema with pydantic validation, so that's where all these are using.
I couldn't find an example of a json file like that though. Would the array be flattened or nested? What is the order? eg for gradients [gx1, gy1, gz1, ... ] or [gx1, gx2, gx3, ...]
So when it's in json, arrays should be flat (not nested) and amenable to reshaping as (natoms, 3). Here's the code: https://github.com/MolSSI/QCElemental/blob/master/qcelemental/models/results.py#L290-L296 . Below, you'll see the arrays while the schemas as pydantic model instances, so they're always in nested form.
Also, it seems the gradient or hessian would be a "return_result", not a "property" - what happens if there's more than one, for example? It seems strange that energies are both properties and return results, but gradients are only return results.
There's deliberate duplication. AtomicResult.return_result is always defined and depends on the AtomicInput.driver setting. AtomicResult.properties.return_energy|gradient|hessian are defined as available. And AtomicResult.properties.scf_total_energy|gradient|hessian https://github.com/MolSSI/QCElemental/blob/master/qcelemental/models/results.py#L102-L106 are defined as available and (if running through QCEngine) depending on how thorough the harness is at collecting properties.
Let me know if there's further questions.
AtomicInput
{'driver': <DriverEnum.gradient: 'gradient'>,
'extras': {},
'id': None,
'keywords': {},
'model': {'basis': 'cc-pvdz', 'method': 'hf'},
'molecule': {'extras': {},
'fix_com': True,
'fix_orientation': True,
'geometry': array([[0., 0., 0.],
[0., 0., 6.]]),
'molecular_charge': 0.0,
'molecular_multiplicity': 1,
'name': 'Ne2',
'provenance': {'creator': 'QCElemental',
'routine': 'qcelemental.molparse.from_schema',
'version': '0.27.1'},
'schema_name': 'qcschema_molecule',
'schema_version': 2,
'symbols': array(['Ne', 'Ne'], dtype='<U2'),
'validated': True},
'protocols': {'stdout': False},
'provenance': {'creator': 'QCElemental',
'routine': 'qcelemental.models.results',
'version': '0.27.1'},
'schema_name': 'qcschema_input',
'schema_version': 1}
AtomicResult
{'driver': <DriverEnum.gradient: 'gradient'>,
'error': None,
'extras': {'qcvars': {'CURRENT DIPOLE': array([0.00000000e+00, 0.00000000e+00, 4.97379915e-14]),
'CURRENT ENERGY': -256.977621071098,
'CURRENT GRADIENT': array([[ 4.81482486e-35, 0.00000000e+00, -6.54470201e-05],
[-4.81482486e-35, 0.00000000e+00, 6.54470212e-05]]),
'CURRENT REFERENCE ENERGY': -256.977621071098,
'DD SOLVATION ENERGY': 0.0,
'HF KINETIC ENERGY': 256.97713564856093,
'HF POTENTIAL ENERGY': -513.9547567196589,
'HF TOTAL ENERGY': -256.977621071098,
'HF TOTAL GRADIENT': array([[ 4.81482486e-35, 0.00000000e+00, -6.54470201e-05],
[-4.81482486e-35, 0.00000000e+00, 6.54470212e-05]]),
'HF VIRIAL RATIO': 2.000001888971701,
'NUCLEAR REPULSION ENERGY': 16.666666666666668,
'ONE-ELECTRON ENERGY': -398.5650589792168,
'PCM POLARIZATION ENERGY': 0.0,
'PE ENERGY': 0.0,
'SCF DIPOLE': array([0.00000000e+00, 0.00000000e+00, 4.97379915e-14]),
'SCF ITERATION ENERGY': -256.977621071098,
'SCF ITERATIONS': 5.0,
'SCF TOTAL ENERGIES': array([-256.97759612, -256.97762005, -256.9776209 , -256.97762107,
-256.97762107, -256.97762107]),
'SCF TOTAL ENERGY': -256.977621071098,
'SCF TOTAL GRADIENT': array([[ 4.81482486e-35, 0.00000000e+00, -6.54470201e-05],
[-4.81482486e-35, 0.00000000e+00, 6.54470212e-05]]),
'TWO-ELECTRON ENERGY': 124.92077124145213}},
'id': None,
'keywords': {},
'model': {'basis': 'cc-pvdz', 'method': 'hf'},
'molecule': {'extras': {},
'fix_com': True,
'fix_orientation': True,
'geometry': array([[0., 0., 0.],
[0., 0., 6.]]),
'molecular_charge': 0.0,
'molecular_multiplicity': 1,
'name': 'Ne2',
'provenance': {'creator': 'QCElemental',
'routine': 'qcelemental.molparse.from_schema',
'version': '0.27.1'},
'schema_name': 'qcschema_molecule',
'schema_version': 2,
'symbols': array(['Ne', 'Ne'], dtype='<U2'),
'validated': True},
'native_files': {},
'properties': {'calcinfo_nalpha': 10,
'calcinfo_natom': 2,
'calcinfo_nbasis': 28,
'calcinfo_nbeta': 10,
'calcinfo_nmo': 28,
'nuclear_repulsion_energy': 16.666666666666668,
'return_energy': -256.977621071098,
'return_gradient': array([[ 4.81482486e-35, 0.00000000e+00, -6.54470201e-05],
[-4.81482486e-35, 0.00000000e+00, 6.54470212e-05]]),
'return_hessian': None,
'scf_dipole_moment': array([0.00000000e+00, 0.00000000e+00, 4.97379915e-14]),
'scf_iterations': 5,
'scf_one_electron_energy': -398.5650589792168,
'scf_total_energy': -256.977621071098,
'scf_total_gradient': array([[ 4.81482486e-35, 0.00000000e+00, -6.54470201e-05],
[-4.81482486e-35, 0.00000000e+00, 6.54470212e-05]]),
'scf_total_hessian': None,
'scf_two_electron_energy': 124.92077124145213},
'protocols': {'stdout': False},
'provenance': {'creator': 'Psi4',
'module': 'scf',
'routine': 'psi4.schema_runner.run_qcschema',
'version': '1.9.1'},
'return_result': array([[ 4.81482486e-35, 0.00000000e+00, -6.54470201e-05],
[-4.81482486e-35, 0.00000000e+00, 6.54470212e-05]]),
'schema_name': 'qcschema_output',
'schema_version': 1,
'stderr': None,
'stdout': None,
'success': True,
'wavefunction': None}
generating script
import pprint
from qcelemental.models import Molecule, AtomicInput, AtomicResult
import psi4
nene = Molecule(geometry= [0, 0, 0, 0, 0, 6], symbols=["Ne", "Ne"], fix_com=True, fix_orientation=True)
atin = {
"molecule": nene,
"driver": "gradient",
"model": {
"method": "hf",
"basis": "cc-pvdz",
},
"protocols": {
"stdout": False,
},
}
atin = AtomicInput(**atin)
print("\n<<< AtomicInput\n")
pprint.pprint(atin.dict(), width=100)
atres = psi4.schema_wrapper.run_qcschema(atin)
print("\n<<< AtomicResult\n")
pprint.pprint(atres.dict(), width=100)
