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.)

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