Creating a QChem Input File

Summary

RevQCMagic has an in-built QChem input file writer. The writer uses three things to generate an input file:

• Configuration (The Configuration object contains all information which gives a state context, i.e. it contains the nuclei involved, the geometry of the nuclei, the basis set for each atom, etc)
• MetaData (For a QChem input file this is a QChemMetadata object which contains all of the calculation parameters.)
• FilePath (The filepath to save the input file to.)

Example: Building Configuration and QChemMetadata Manually

from qcmagic.core.cspace.configurationspace import Configuration
from qcmagic.core.cspace.molecule.atom import Atom
from qcmagic.core.cspace.molecule.atomfragment import AtomFragment
from qcmagic.core.cspace.molecule.geometry import Geometry
from qcmagic.auxiliary.linearalgebra3d import Vector3D
from qcmagic.auxiliary.qcmagic_standards import ANGSTROM_TO_BOHR

from qcmagic.interfaces.input_writers.qchem_input_writer.qchem_input_writer import (
    QChemMetadata,
    QChemInputWriter,
)

import numpy as np

def generate_H2_STO3G_Qchem_input_files():
    # Each input file will have the same calculation parameters
    parameters = {
            'job_type':         'sp',
            'basis':            'STO-3G',
            'method':           'HF',
            'symmetry_ignore':  'true',
            'scf_convergence':  '11',
            'thresh':           '14',
            'unrestricted':     'false',
            'iqmol_fchk':       'true',
            'max_scf_cycles':   '500',
            'mem_total':        '5976800'
    }

    # Create a QChemMetadata object to store the calculation parameters
    # as well as the total charge (0) and the multiplicity (2*spin_proj+1)
    metadata = QChemMetadata()
    metadata.update_from_dictionary(parameters)
    metadata.state_properties = {"total_charge": 0, "spin_projection": 0}

    # Considering H-H bond lengths from 0.2 to 2,1 Angstroms
    for radius in np.linspace(0.2, 2.1, 20):
        # RevQCMagic uses Bohr whereas QChem uses Angstroms
        radius_bohr = radius * ANGSTROM_TO_BOHR

        # Every new geometry requires a new Configuration
        config = Configuration()

        # Generate a co-ordinate system for holding the molecule geometry
        coordsystem = AtomFragment(parent=None)
        H2 = AtomFragment(coordsystem)

        # Add the two hydrogen atoms, note the required use of Vector3D
        H2.add(Atom('H', H2, Vector3D([0.00, 0.00, -radius_bohr/2])))
        H2.add(Atom('H', H2, Vector3D([0.00, 0.00, +radius_bohr/2])))

        # Turn the AtomFragment object into a Geometry SubConfiguration, then
        # add to the total Configuration
        config.add_subconfiguration(Geometry([(None, None)], H2))

        # Generate the input file writer using the Configuration, and use the 
        # QChemMetadata to write the specific input file to the specified filepath.
        # If wanted to do different calculations on the same Configuration then could
        # create different QChemMetadata objects and call the writer with each of them.
        writer = QChemInputWriter(config)
        filepath = f'H2_STO-3G_radius_{radius:.1f}A.inp'
        writer(filepath, metadata)

generate_H2_STO3G_Qchem_input_files()

The above code generates input files which all look like the one given below.

$molecule
    0 1
    H    +0.00000000000000    +0.00000000000000    -0.10000000000000
    H    +0.00000000000000    +0.00000000000000    +0.10000000000000
$end
$rem
    basis STO-3G
    job_type sp
    method hf
    symmetry_ignore true
    scf_convergence 11
    thresh 14
    unrestricted false
    iqmol_fchk true
    max_scf_cycles 500
    mem_total 5976800
$end