Source code for gaussian_n

# Gn theory.  

import PsiMod
import re
import os
import math
import warnings
from driver import *
from molutil import *
from text import *
from procutil import *
from physconst import *
# never import aliases into this file

def run_gaussian_2(name, **kwargs):

    # throw an exception for open-shells
    if (PsiMod.get_option('SCF','REFERENCE') != 'RHF' ):
        raise ValidationError("""g2 computations require "reference rhf".""")

    # stash user options:
    optstash = OptionsState(
        ['FNOCC','COMPUTE_TRIPLES'],
        ['FNOCC','COMPUTE_MP4_TRIPLES'],
        ['FREEZE_CORE'],
        ['SCF','SCF_TYPE'])

    # override default scf_type
    PsiMod.set_local_option('SCF','SCF_TYPE','OUT_OF_CORE')

    # optimize geometry at scf level
    PsiMod.clean()
    PsiMod.set_global_option('BASIS',"6-31G(D)")
    optimize('scf')
    PsiMod.clean()

    # scf frequencies for zpe
    frequency('scf')

    # thermodynamic properties
    du = PsiMod.get_variable('INTERNAL ENERGY CORRECTION')
    dh = PsiMod.get_variable('ENTHALPY CORRECTION')
    dg = PsiMod.get_variable('GIBBS FREE ENERGY CORRECTION')

    ref     = PsiMod.wavefunction()
    freqs   = ref.frequencies()
    nfreq   = freqs.dim(0)
    freqsum = 0.0
    for i in range (0,nfreq):
        freqsum += freqs.get(i)
    zpe = freqsum / psi_hartree2wavenumbers * 0.8929 * 0.5
    PsiMod.clean()

    # optimize geometry at mp2 (no frozen core) level
    # note: freeze_core isn't an option in MP2
    PsiMod.set_global_option('FREEZE_CORE',"FALSE")
    optimize('conv-mp2')
    PsiMod.clean()

    # qcisd(t)
    PsiMod.set_local_option('FNOCC','COMPUTE_MP4_TRIPLES',"TRUE")
    PsiMod.set_global_option('FREEZE_CORE',"TRUE")
    PsiMod.set_global_option('BASIS',"6-311G(D_P)")
    run_fnocc('qcisd(t)',**kwargs)

    # HLC: high-level correction based on number of valence electrons
    ref    = PsiMod.wavefunction()
    nirrep = ref.nirrep()
    frzcpi = ref.frzcpi()
    nfzc = 0
    for i in range (0,nirrep):
        nfzc += frzcpi[i]
    nalpha = ref.nalpha() - nfzc
    nbeta  = ref.nbeta() - nfzc
    # hlc of gaussian-2
    hlc = -0.00481 * nalpha -0.00019 * nbeta
    # hlc of gaussian-1
    hlc1 = -0.00614 * nalpha

    eqci_6311gdp = PsiMod.get_variable("QCISD(T) TOTAL ENERGY")
    emp4_6311gd  = PsiMod.get_variable("MP4 TOTAL ENERGY")
    emp2_6311gd  = PsiMod.get_variable("MP2 TOTAL ENERGY")
    PsiMod.clean()

    # correction for diffuse functions
    PsiMod.set_global_option('BASIS',"6-311+G(D_P)")
    energy('mp4')
    emp4_6311pg_dp = PsiMod.get_variable("MP4 TOTAL ENERGY")
    emp2_6311pg_dp = PsiMod.get_variable("MP2 TOTAL ENERGY")
    PsiMod.clean()

    # correction for polarization functions
    PsiMod.set_global_option('BASIS',"6-311G(2DF_P)")
    energy('mp4')
    emp4_6311g2dfp = PsiMod.get_variable("MP4 TOTAL ENERGY")
    emp2_6311g2dfp = PsiMod.get_variable("MP2 TOTAL ENERGY")
    PsiMod.clean()

    # big basis mp2
    PsiMod.set_global_option('BASIS',"6-311+G(3DF_2P)")
    run_fnocc('_mp2',**kwargs)
    emp2_big = PsiMod.get_variable("MP2 TOTAL ENERGY")
    PsiMod.clean()

    eqci       = eqci_6311gdp
    e_delta_g2 = emp2_big + emp2_6311gd - emp2_6311g2dfp - emp2_6311pg_dp
    e_plus     = emp4_6311pg_dp - emp4_6311gd
    e_2df      = emp4_6311g2dfp - emp4_6311gd

    eg2 = eqci + e_delta_g2 + e_plus + e_2df
    eg2_mp2_0k = eqci + (emp2_big - emp2_6311gd) + hlc + zpe

    PsiMod.print_out('\n')
    PsiMod.print_out('  ==>  G1/G2 Energy Components  <==\n')
    PsiMod.print_out('\n')
    PsiMod.print_out('        QCISD(T):            %20.12lf\n' % eqci)
    PsiMod.print_out('        E(Delta):            %20.12lf\n' % e_delta_g2)
    PsiMod.print_out('        E(2DF):              %20.12lf\n' % e_2df)
    PsiMod.print_out('        E(+):                %20.12lf\n' % e_plus)
    PsiMod.print_out('        E(G1 HLC):           %20.12lf\n' % hlc1)
    PsiMod.print_out('        E(G2 HLC):           %20.12lf\n' % hlc)
    PsiMod.print_out('        E(ZPE):              %20.12lf\n' % zpe)
    PsiMod.print_out('\n')
    PsiMod.print_out('  ==>  0 Kelvin Results  <==\n')
    PsiMod.print_out('\n')
    eg2_0k = eg2 + zpe + hlc
    PsiMod.print_out('        G1:                  %20.12lf\n' % (eqci + e_plus + e_2df + hlc1 + zpe))
    PsiMod.print_out('        G2(MP2):             %20.12lf\n' % eg2_mp2_0k)
    PsiMod.print_out('        G2:                  %20.12lf\n' % eg2_0k)

    PsiMod.set_variable("G1 TOTAL ENERGY",eqci + e_plus + e_2df + hlc1 + zpe)
    PsiMod.set_variable("G2 TOTAL ENERGY",eg2_0k)
    PsiMod.set_variable("G2(MP2) TOTAL ENERGY",eg2_mp2_0k)

    PsiMod.print_out('\n')
    T = PsiMod.get_global_option('T')
    PsiMod.print_out('  ==>  %3.0lf Kelvin Results  <==\n'% T)
    PsiMod.print_out('\n')

    internal_energy = eg2_mp2_0k + du - zpe / 0.8929
    enthalpy        = eg2_mp2_0k + dh - zpe / 0.8929
    gibbs           = eg2_mp2_0k + dg - zpe / 0.8929

    PsiMod.print_out('        G2(MP2) energy:      %20.12lf\n' % internal_energy )
    PsiMod.print_out('        G2(MP2) enthalpy:    %20.12lf\n' % enthalpy)
    PsiMod.print_out('        G2(MP2) free energy: %20.12lf\n' % gibbs)
    PsiMod.print_out('\n')

    PsiMod.set_variable("G2(MP2) INTERNAL ENERGY",internal_energy)
    PsiMod.set_variable("G2(MP2) ENTHALPY",enthalpy)
    PsiMod.set_variable("G2(MP2) FREE ENERGY",gibbs)

    internal_energy = eg2_0k + du - zpe / 0.8929
    enthalpy        = eg2_0k + dh - zpe / 0.8929
    gibbs           = eg2_0k + dg - zpe / 0.8929

    PsiMod.print_out('        G2 energy:           %20.12lf\n' % internal_energy )
    PsiMod.print_out('        G2 enthalpy:         %20.12lf\n' % enthalpy)
    PsiMod.print_out('        G2 free energy:      %20.12lf\n' % gibbs)

    PsiMod.set_variable("CURRENT ENERGY",eg2_0k)

    PsiMod.set_variable("G2 INTERNAL ENERGY",internal_energy)
    PsiMod.set_variable("G2 ENTHALPY",enthalpy)
    PsiMod.set_variable("G2 FREE ENERGY",gibbs)

    PsiMod.clean()

    optstash.restore()

    # return 0K g2 results
    return eg2_0k

# aliases for g2
procedures['energy']['gaussian-2'] = run_gaussian_2
procedures['energy']['g2']         = run_gaussian_2