How to store gradients and hessians?

[aizvorski]

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?

[loriab]

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)