psi4-data 1:0.3-5 / usr / share / psi / python / wrappers_cfour.py

#
#@BEGIN LICENSE
#
# PSI4: an ab initio quantum chemistry software package
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
#@END LICENSE
#

"""Module with functions for Psi4/Cfour interface. Portions that require
calls to Boost Python psi4 module are here, otherwise in qcdb module.
Also calls to qcdb module are here and not elsewhere in driver.
Organizationally, this module isolates qcdb code from psi4 code.

"""
from __future__ import print_function
import shutil
import os
import subprocess
import re
import inspect
import glob
import shelve
import datetime
import difflib
#CUimport psi4
#CUimport p4const
#CUimport p4util
#CUimport qcdb
#CUfrom p4regex import *
#from extend_Molecule import *
#CUfrom molutil import *
#CUfrom functional import *
from driver import *
from p4xcpt import *
# never import driver, wrappers, or aliases into this file


def run_cfour_module(xmod):
    # Find environment by merging PSIPATH and PATH environment variables
    lenv = os.environ
    lenv['PATH'] = ':'.join([os.path.abspath(x) for x in os.environ.get('PSIPATH', '').split(':')]) + ':' + lenv.get('PATH') + ':' + psi4.Process.environment["PSIDATADIR"] + '/basis' + ':' + psi4.psi_top_srcdir() + '/lib/basis'

    # Call executable xcfour, directing cfour output to the psi4 output file
    try:
        retcode = subprocess.Popen([xmod], bufsize=0, stdout=subprocess.PIPE, env=lenv)
    except OSError as e:
        sys.stderr.write('Program %s not found in path or execution failed: %s\n' % (cfour_executable, e.strerror))
        #p4out.write('Program %s not found in path or execution failed: %s\n' % (cfour_executable, e.strerror))
        message = ('Program %s not found in path or execution failed: %s\n' % (cfour_executable, e.strerror))
        raise ValidationError(message)        

    c4out = ''
    while True:
        data = retcode.stdout.readline()
        if not data:
            break
        #if psi4.outfile_name() == 'stdout':
        #    sys.stdout.write(data)
        #else:
        #    p4out.write(data)
        #    p4out.flush()
        c4out += data
    #internal_p4c4_info['output'] = c4out
    return c4out


def vpt2(name, **kwargs):
    """Perform vibrational second-order perturbation computation through
    Cfour to get anharmonic frequencies. This version uses c4 for the disp
    and pt2 but gets gradients from p4.

    :type c4full: :ref:`boolean <op_py_boolean>`
    :param c4full: ``'on'`` || |dl| ``'off'`` |dr|

        Indicates whether when *name* indicates a Cfour method and *mode*
        indicates a sow/reap approach, sown files are direct ZMAT files
        and FJOBARC files are expected to reap, so that Cfour only, not
        Cfour-through-Psi4, is needed for distributed jobs.

    .. caution:: Some features are not yet implemented. Buy a developer a coffee.

       - Presently uses all gradients. Could mix in analytic 2nd-derivs.

       - Collect resutls.

       - Manage scratch / subdir better.

       - Untangle CCSD(T) vs CCSD[T] and FJOBARC issue

       - Allow CFOUR_BASIS

       - Consider forcing some tighter convcrit, c4 and p4

       - sow/reap

       - mixed ang/bohr signals

       - error by converting to ang in psi?

       - Expand CURRENT DIPOLE XYZ beyond SCF

       - Remember additional FJOBARC record TOTENER2 if EXCITE .ne. NONE

       - S/R P4grad

       - S/R C4grad

       - C P4grad

       - C C4grad

       - switch C --> S/R with recovery using shelf

       - pure C mode where only need P4 for wrapper

    """
    lowername = name.lower()
    kwargs = p4util.kwargs_lower(kwargs)

    optstash = p4util.OptionsState(
        ['BASIS'])

    # Option mode of operation- whether vpt2 run in one job or files farmed out
    if not('vpt2_mode' in kwargs):
        if ('mode' in kwargs):
            kwargs['vpt2_mode'] = kwargs['mode']
            del kwargs['mode']
        else:
            kwargs['vpt2_mode'] = 'continuous'

    # Switches for route through code- S/R or continuous & Psi4 or Cfour gradients
    isSowReap = True if kwargs['vpt2_mode'].lower() == 'sowreap' else False
    isC4notP4 = bool(re.match('cfour', lowername)) or bool(re.match('c4-', lowername))
    isC4fully = True if ('c4full' in kwargs and yes.match(str(kwargs['c4full'])) and isC4notP4 and isSowReap) else False

    # Save submission directory and basis set
    current_directory = os.getcwd()
    user_basis = psi4.get_global_option('BASIS')

    # Open data persistence shelf- vital for sowreap, checkpoint for continuouw
    shelf = shelve.open(current_directory + '/' + os.path.splitext(psi4.outfile_name())[0] + '.shelf', writeback=True)

    # Cfour keywords to request vpt2 analysis through findif gradients
    psi4.set_local_option('CFOUR', 'CFOUR_VIBRATION', 'FINDIF')
    psi4.set_local_option('CFOUR', 'CFOUR_FREQ_ALGORITHM', 'PARALLEL')
    psi4.set_local_option('CFOUR', 'CFOUR_ANH_ALGORITHM', 'PARALLEL')
    psi4.set_local_option('CFOUR', 'CFOUR_ANHARMONIC', 'VPT2')
    psi4.set_local_option('CFOUR', 'CFOUR_FD_PROJECT', 'OFF')

    # When a Psi4 method is requested for vpt2, a skeleton of
    #   computations in Cfour is still required to hang the gradients
    #   upon. The skeleton is as cheap as possible (integrals only
    #   & sto-3g) and set up here.
    if isC4notP4:
        skelname = lowername
    else:
        skelname = 'c4-scf'
        psi4.set_global_option('BASIS', 'STO-3G')
    #    P4  'c4-scf'/'cfour'CALC_LEVEL      lowername  # temporary
    #    C4  lowername                       cfour{}  # temporary

    if 'status' not in shelf:
        shelf['status'] = 'initialized'
        shelf['linkage'] = os.getpid()
        shelf['zmat'] = {}  # Cfour-generated ZMAT files with finite difference geometries
        shelf['fjobarc'] = {}  # Cfour- or Psi4-generated ascii files with packaged gradient results
        shelf.sync()
    else:
        pass
        # how decide whether to use. keep precedent of intco.dat in mind

    # Construct and move into directory job scratch / cfour scratch / harm
    psioh = psi4.IOManager.shared_object()
    psio = psi4.IO.shared_object()
    os.chdir(psioh.get_default_path())  # psi_scratch
    cfour_tmpdir = kwargs['path'] if 'path' in kwargs else \
        'psi.' + str(os.getpid()) + '.' + psio.get_default_namespace() + \
        '.cfour.' + str(random.randint(0, 99999))
    if not os.path.exists(cfour_tmpdir):
        os.mkdir(cfour_tmpdir)
    os.chdir(cfour_tmpdir)  # psi_scratch/cfour
    if not os.path.exists('harm'):
        os.mkdir('harm')
    os.chdir('harm')  # psi_scratch/cfour/harm

    psioh.set_specific_retention(32, True)  # temporary, to track p4 scratch
    #shelf['status'] = 'anharm_jobs_sown'  # temporary to force backtrack
    print('STAT', shelf['status'])  # temporary

    # Generate the ZMAT input file in scratch
    with open('ZMAT', 'w') as handle:
        cfour_infile = write_zmat(skelname, 1)
        handle.write(cfour_infile)
    print('\n====== Begin ZMAT input for CFOUR ======')
    print(open('ZMAT', 'r').read())
    print('======= End ZMAT input for CFOUR =======\n')
    shelf['genbas'] = open('GENBAS', 'r').read()

    # Check existing shelf consistent with generated ZMAT, store
    if ('000-000' in shelf['zmat']) and (shelf['zmat']['000-000'] != cfour_infile):
        diff = difflib.Differ().compare(shelf['zmat']['000-000'].splitlines(), cfour_infile.splitlines())
        raise ValidationError("""Input file translated to Cfour ZMAT does not match ZMAT stored in shelf.\n\n""" +
            '\n'.join(list(diff)))
    shelf['zmat']['000-000'] = cfour_infile
    shelf.sync()

    # Reset basis after Cfour skeleton seeded
    psi4.set_global_option('BASIS', user_basis)

    if shelf['status'] == 'initialized':
        p4util.banner(' VPT2 Setup: Harmonic ')

        # Generate the displacements that will form the harmonic freq
        os.chdir(psioh.get_default_path() + cfour_tmpdir + '/harm')  # psi_scratch/cfour/harm
        with open('partial.out', 'w') as handle:
            handle.write(run_cfour_module('xjoda'))
            handle.write(run_cfour_module('xsymcor'))

        # Read the displacements that will form the harmonic freq
        zmats0N = ['000-' + item[-3:] for item in sorted(glob.glob('zmat*'))]
        for zm12 in zmats0N:
            zm1, zm2 = zm12.split('-')
            with open('zmat' + zm2, 'r') as handle:
                shelf['zmat'][zm12] = handle.read()
                shelf.sync()
            psi4.print_out('  CFOUR scratch file %s for %s-%s has been read\n' % ('zmat' + zm2, zm1, zm2))
            psi4.print_out('%s\n' % shelf['zmat'][zm12])

        # S/R: Write distributed input files for harmonic freq
        if isSowReap:
            os.chdir(current_directory)
            inputSansMol = p4util.format_currentstate_for_input(gradient, lowername, allButMol=True, **kwargs)
            for zm12 in zmats0N:
                zm1, zm2 = zm12.split('-')

                ifile = vpt2_sow_files(zm12, shelf['linkage'], isC4notP4, isC4fully,
                    shelf['zmat'][zm12], inputSansMol, shelf['genbas'])

                with open('VPT2-' + zm12 + '.in', 'w') as handle:
                    handle.write(ifile)

            msg = vpt2_instructions('harmonic', current_directory, zmats0N)
            psi4.print_out(msg)
            print(msg)

        shelf['status'] = 'harm_jobs_sown'

        # S/R: Pause for distributed calculations
        if isSowReap:
            shelf.close()
            return 0.0

    if shelf['status'] == 'harm_jobs_sown':
        zmats0N = [item for item in sorted(shelf['zmat'].keys()) if (item[:3] == '000' and item[-3:] != '000')]

        # S/R: Check that distributed calcs all completed correctly
        if isSowReap:
            msg = vpt2_instructions('harmonic', current_directory, zmats0N)
            psi4.print_out(msg)
            isOk, msg = sown_jobs_status(current_directory, 'VPT2', zmats0N, reap_job_validate,
                shelf['linkage'], ['CURRENT ENERGY', 'CURRENT DIPOLE', 'CURRENT GRADIENT'])
            psi4.print_out(msg)
            print(msg)
            if not isOk:
                shelf.close()
                return 0.0

        # Collect all results from gradients forming the harmonic freq
        for zm12 in zmats0N:
            zm1, zm2 = zm12.split('-')
            if zm12 not in shelf['fjobarc']:
                p4util.banner(' VPT2 Computation: %s ' % (zm12))
                print(' VPT2 Computation: %s ' % (zm12))

                fjobarc = vpt2_reaprun_files(zm12, shelf['linkage'], isSowReap, isC4notP4, isC4fully,
                    shelf['zmat'][zm12], current_directory, psioh.get_default_path(), cfour_tmpdir,
                    lowername, kwargs)
                shelf['fjobarc'][zm12] = fjobarc
                shelf.sync()
        shelf['status'] = 'harm_jobs_reaped'

    if shelf['status'] == 'harm_jobs_reaped':
        zmats0N = [item for item in sorted(shelf['zmat'].keys()) if (item[:3] == '000' and item[-3:] != '000')]

        p4util.banner(' VPT2 Results: Harmonic ')

        # Process the gradients into harmonic freq
        os.chdir(psioh.get_default_path() + cfour_tmpdir + '/harm')  # psi_scratch/cfour/harm
        harmout = run_cfour_module('xjoda')
        harmout += run_cfour_module('xsymcor')
        for zm12 in zmats0N:
            zm1, zm2 = zm12.split('-')
            with open('FJOBARC', 'w') as handle:
                handle.write(shelf['fjobarc'][zm12])
            harmout += run_cfour_module('xja2fja')
            harmout += run_cfour_module('xsymcor')
            shutil.move('FJOBARC', 'fja.' + zm12)
            try:
                os.remove('zmat' + zm2)
            except OSError:
                pass
        harmout += run_cfour_module('xjoda')
        harmout += run_cfour_module('xcubic')

        psi4.print_out(harmout)
        with open('harm.out', 'w') as handle:
            handle.write(harmout)

        # Generate displacements along harmonic normal modes
        zmatsN0 = [item[-3:] for item in sorted(glob.glob('zmat*'))]
        os.chdir('..')  # psi_scratch/cfour
        for zm1 in zmatsN0:
            zm12 = zm1 + '-000'
            with open(psioh.get_default_path() + cfour_tmpdir + '/harm/zmat' + zm1, 'r') as handle:
                shelf['zmat'][zm12] = handle.read()
                shelf.sync()
                psi4.print_out('  CFOUR scratch file %s for %s has been read\n' % ('zmat' + zm1, zm12))
                psi4.print_out('%s\n' % shelf['zmat'][zm12])

            # Collect displacements along the normal coordinates generated by the harmonic freq.
            #   Further harmonic freqs are to be run at each of these to produce quartic force field.
            #   To carry these out, generate displacements for findif by gradient at each displacement.
            if os.path.exists(zm1):
                shutil.rmtree(zm1)
            os.mkdir(zm1)
            os.chdir(zm1)  # psi_scratch/cfour/004
            with open('ZMAT', 'w') as handle:
                handle.write(shelf['zmat'][zm12])
            shutil.copy2('../harm/GENBAS', 'GENBAS')  # ln -s $ecpdir/ECPDATA $j/ECPDATA
            with open('partial.out', 'w') as handle:
                handle.write(run_cfour_module('xjoda'))
                handle.write(run_cfour_module('xsymcor'))

            # Read the displacements that will form the anharmonic freq
            zmatsNN = [item[-3:] for item in sorted(glob.glob('zmat*'))]
            for zm2 in zmatsNN:
                zm12 = zm1 + '-' + zm2
                with open(psioh.get_default_path() + cfour_tmpdir + '/' + zm1 + '/zmat' + zm2, 'r') as handle:
                    shelf['zmat'][zm12] = handle.read()
                    shelf.sync()
                    psi4.print_out('  CFOUR scratch file %s for %s has been read\n' % ('zmat' + zm2, zm12))
                    psi4.print_out('%s\n' % shelf['zmat'][zm12])
            os.chdir('..')  # psi_scratch/cfour

        zmatsNN = [item for item in sorted(shelf['zmat'].keys()) if (item[:3] != '000' and item[-3:] != '000')]

        # S/R: Write distributed input files for anharmonic freq
        if isSowReap:
            os.chdir(current_directory)
            inputSansMol = p4util.format_currentstate_for_input(gradient, lowername, allButMol=True, **kwargs)
            for zm12 in zmatsNN:
                zm1, zm2 = zm12.split('-')

                ifile = vpt2_sow_files(zm12, shelf['linkage'], isC4notP4, isC4fully,
                    shelf['zmat'][zm12], inputSansMol, shelf['genbas'])
                # GENBAS needed here

                with open('VPT2-' + zm12 + '.in', 'w') as handle:
                    handle.write(ifile)

            msg = vpt2_instructions('anharmonic', current_directory, zmatsNN)
            psi4.print_out(msg)
            print(msg)

        shelf['status'] = 'anharm_jobs_sown'

        # S/R: Pause for distributed calculations
        if isSowReap:
            shelf.close()
            return 0.0

    if shelf['status'] == 'anharm_jobs_sown':
        zmatsNN = [item for item in sorted(shelf['zmat'].keys()) if (item[:3] != '000' and item[-3:] != '000')]

        # S/R: Check that distributed calcs all completed correctly
        if isSowReap:
            msg = vpt2_instructions('anharmonic', current_directory, zmatsNN)
            psi4.print_out(msg)
            isOk, msg = sown_jobs_status(current_directory, 'VPT2', zmatsNN, 
                reap_job_validate, shelf['linkage'], 
                ['CURRENT ENERGY', 'CURRENT DIPOLE', 'CURRENT GRADIENT'])
            psi4.print_out(msg)
            print(msg)
            if not isOk:
                shelf.close()
                return 0.0

        # Collect all results from gradients forming the anharmonic freq
        for zm12 in zmatsNN:
            zm1, zm2 = zm12.split('-')
            if zm12 not in shelf['fjobarc']:
                p4util.banner(' VPT2 Computation: %s ' % (zm12))
                print(' VPT2 Computation: %s ' % (zm12))

                fjobarc = vpt2_reaprun_files(zm12, shelf['linkage'], isSowReap, isC4notP4, isC4fully,
                    shelf['zmat'][zm12], current_directory, psioh.get_default_path(), cfour_tmpdir,
                    lowername, kwargs)
                shelf['fjobarc'][zm12] = fjobarc
                shelf.sync()
        shelf['status'] = 'anharm_jobs_reaped'

    if shelf['status'] == 'anharm_jobs_reaped':
        zmats0N = [item for item in sorted(shelf['zmat'].keys()) if (item[:3] == '000' and item[-3:] != '000')]
        zmatsN0 = [item for item in sorted(shelf['zmat'].keys()) if (item[:3] != '000' and item[-3:] == '000')]
        zmatsNN = [item for item in sorted(shelf['zmat'].keys()) if (item[:3] != '000' and item[-3:] != '000')]

        p4util.banner(' VPT2 Results: Harmonic ')

        # Process the gradients into harmonic freq
        os.chdir(psioh.get_default_path() + cfour_tmpdir)  # psi_scratch/cfour
        if os.path.exists('anharm'):
            shutil.rmtree('anharm')
        os.mkdir('anharm')
        os.chdir('harm')  # psi_scratch/cfour/harm
        run_cfour_module('xclean')
        anharmout = run_cfour_module('xjoda')
        anharmout += run_cfour_module('xsymcor')

        for zm12 in zmats0N:
            zm1, zm2 = zm12.split('-')
            with open('FJOBARC', 'w') as handle:
                handle.write(shelf['fjobarc'][zm12])
            anharmout += run_cfour_module('xja2fja')
            anharmout += run_cfour_module('xsymcor')
            shutil.move('FJOBARC', 'fja.' + zm12)
        anharmout += run_cfour_module('xjoda')
        anharmout += run_cfour_module('xcubic')

        psi4.print_out(anharmout)
        with open('harm.out', 'w') as handle:
            handle.write(anharmout)

        # Process the gradients into harmonic freq at each normco displaced point
        os.chdir('..')  # psi_scratch/cfour
        for zm11 in zmatsN0:
            zm1 = zm11[:3]
            if os.path.exists(zm1):
                shutil.rmtree(zm1)
            os.mkdir(zm1)
            os.chdir(zm1)  # psi_scratch/cfour/004

            run_cfour_module('xclean')
            with open('ZMAT', 'w') as handle:
                handle.write(shelf['zmat'][zm11])
            shutil.copy2('../harm/GENBAS', 'GENBAS')
            anharmout = run_cfour_module('xjoda')
            anharmout += run_cfour_module('xsymcor')

            for zm22 in [item for item in zmatsNN if (item[:3] == zm1 and item[-3:] != '000')]:
                zm2 = zm22[-3:]
                zm12 = zm1 + '-' + zm2
                print(zm12)
                with open('FJOBARC', 'w') as handle:
                    handle.write(shelf['fjobarc'][zm12])
                anharmout += run_cfour_module('xja2fja')
                anharmout += run_cfour_module('xsymcor')
                shutil.move('FJOBARC', 'fja.' + zm12)
            anharmout += run_cfour_module('xjoda')
            anharmout += run_cfour_module('xja2fja')
            with open('FJOBARC', 'r') as handle:
                shelf['fjobarc'][zm11] = handle.read()
                shelf.sync()

            psi4.print_out(anharmout)
            with open('partial.out', 'w') as handle:
                handle.write(anharmout)

            os.chdir('..')  # psi_scratch/cfour

        # Process the harmonic freqs at normco displacements into anharmonic freq
        p4util.banner(' VPT2 Results: Anharmonic ')

        os.chdir('anharm')  # psi_scratch/cfour/anharm
        shutil.copy2('../harm/JOBARC', 'JOBARC')
        shutil.copy2('../harm/JAINDX', 'JAINDX')

        for zm12 in zmatsN0:
            with open('FJOBARC', 'w') as handle:
                handle.write(shelf['fjobarc'][zm12])
            anharmout = run_cfour_module('xja2fja')
            anharmout += run_cfour_module('xcubic')
            shutil.move('FJOBARC', 'fja.' + zm12)

        psi4.print_out(anharmout)
        with open('anharm.out', 'w') as handle:
            handle.write(anharmout)

        shelf['status'] = 'vpt2_completed'

    # Finish up
    os.chdir(current_directory)
    shelf.close()
    optstash.restore()


def vpt2_sow_files(item, linkage, isC4notP4, isC4fully, zmat, inputSansMol, inputGenbas):
    """Provided with the particular displacement number *item* and the
    associated *zmat* file contents and *linkage*, and common contents
    *inputSansMol*, returns contents of input file to be sown.

    """
    inputReapOrders = r"""
print_variables()

print_out('VPT2 RESULT: linkage {0} for item {1} yields CURRENT ENERGY being %r\n' % (get_variable('CURRENT ENERGY')))
print_out('VPT2 RESULT: linkage {0} for item {1} yields CURRENT GRADIENT being %r\n' % (p4util.mat2arr(psi4.get_gradient())))
print_out('VPT2 RESULT: linkage {0} for item {1} yields CURRENT DIPOLE being [%r, %r, %r]\n' % (get_variable('CURRENT DIPOLE X'), get_variable('CURRENT DIPOLE Y'), get_variable('CURRENT DIPOLE Z')))
""".format(linkage, item)

    # Direct Cfour for gradients
    if isC4fully:
        inputString = zmat
        with open('VPT2-GENBAS', 'w') as handle:
            handle.write(inputGenbas)

    # Cfour for gradients
    elif isC4notP4:
        # GENBAS needed here
        inputString = 'extracted_genbas = """\n' + inputGenbas.replace('\n\n', '\nblankline\n') + '\n"""\n\n'
        inputString += """cfour {\n%s\n}\n\nenergy('cfour', genbas=extracted_genbas)\n\n""" % (zmat)
        inputString += inputReapOrders
        inputString += r"""
print_out('VPT2 RESULT: linkage {0} for item {1} yields CURRENT MOLECULE being %r\n' % (get_active_molecule().create_psi4_string_from_molecule()))
""".format(linkage, item)

    # Psi4 for gradients
    else:
        inputString = p4util.format_molecule_for_input(
            qcdb.cfour.harvest_zmat(zmat).create_psi4_string_from_molecule(),
            name='disp' + item[:3] + item[-3:])
        inputString += inputSansMol
        inputString += inputReapOrders

    return inputString


def vpt2_reaprun_files(item, linkage, isSowReap, isC4notP4, isC4fully, zmat, outdir, scrdir, c4scrdir, lowername, kwargs):
    """Provided with the particular displacement number *item* and the
    associated *zmat* file with geometry and *linkage*, returns the
    FJOBARC contents. Depending on the mode settings of *isC4notP4*,
    *isSowReap*, and *isC4fully*, either runs (using *lowername* and
    *kwargs*) or reaps contents. *outdir* is where psi4 was invoked,
    *scrdir* is the psi4 scratch directory, and *c4scrdir* is Cfour
    scratch directory within.

    """
    os.chdir(outdir)  # current_directory
    # Extract qcdb.Molecule at findif orientation
    zmmol = qcdb.cfour.harvest_zmat(zmat)

    # Cfour S/R Direct for gradients
    if isC4fully:
        with open('VPT2-' + item + '.fja', 'r') as handle:
            fjobarc = handle.read()

    # Cfour for gradients
    elif isC4notP4:

        # S/R: Reap results from output file
        if isSowReap:
            isOk, msg, results = reap_job_validate(outdir, 'VPT2', item, linkage,
                ['CURRENT ENERGY', 'CURRENT DIPOLE', 'CURRENT GRADIENT', 'CURRENT MOLECULE'])
            if not isOk:
                raise ValidationError(msg)

            fje = results['CURRENT ENERGY']
            fjgrd = results['CURRENT GRADIENT']
            fjdip = [item / p4const.psi_dipmom_au2debye for item in results['CURRENT DIPOLE']]
            c4mol = qcdb.Molecule(results['CURRENT MOLECULE'])
            c4mol.update_geometry()

        # C: Run the job and collect results
        else:
            # Prepare Cfour skeleton calc directory
            os.chdir(scrdir + c4scrdir)  # psi_scratch/cfour
            if os.path.exists('scr.' + item):
                shutil.rmtree('scr.' + item)
            os.mkdir('scr.' + item)
            os.chdir('scr.' + item)  # psi_scratch/cfour/scr.000-004
            with open('ZMAT', 'w') as handle:
                handle.write(zmat)
            shutil.copy2('../harm/GENBAS', 'GENBAS')

            #os.chdir(scrdir + '/scr.' + item)
            #run_cfour_module('xja2fja')
            #with open('FJOBARC', 'r') as handle:
            #    fjobarc = handle.read()

            # Run Cfour calc using ZMAT & GENBAS in scratch, outdir redirects to outfile
            os.chdir(outdir)  # current_directory
            psi4.get_active_molecule().set_name('blank_molecule_psi4_yo')
            energy('cfour', path=c4scrdir + '/scr.' + item)
#            os.chdir(scrdir + '/scr.' + item)

            fje = psi4.get_variable('CURRENT ENERGY')
            fjgrd = p4util.mat2arr(psi4.get_gradient())
            fjdip = [psi4.get_variable('CURRENT DIPOLE X') / p4const.psi_dipmom_au2debye,
                     psi4.get_variable('CURRENT DIPOLE Y') / p4const.psi_dipmom_au2debye,
                     psi4.get_variable('CURRENT DIPOLE Z') / p4const.psi_dipmom_au2debye]
            c4mol = qcdb.Molecule(psi4.get_active_molecule().create_psi4_string_from_molecule())
            c4mol.update_geometry()

        # Get map btwn ZMAT and C4 orientation, then use it, grad and dipole to forge FJOBARC file
        fjobarc = qcdb.cfour.format_fjobarc(fje,
            *qcdb.cfour.backtransform(chgeMol=zmmol, permMol=c4mol), gradient=fjgrd, dipole=fjdip)

    # Psi4 for gradients
    else:
        # Prepare Cfour skeleton calc directory
        os.chdir(scrdir + c4scrdir)  # psi_scratch/cfour
        if os.path.exists('scr.' + item):
            shutil.rmtree('scr.' + item)
        os.mkdir('scr.' + item)
        os.chdir('scr.' + item)  # psi_scratch/cfour/scr.000-004
        with open('ZMAT', 'w') as handle:
            handle.write(zmat)
        shutil.copy2('../harm/GENBAS', 'GENBAS')

        # Run Cfour skeleton calc and extract qcdb.Molecule at needed C4 orientation
        with open('partial.out', 'w') as handle:
            handle.write(run_cfour_module('xjoda'))
            handle.write(run_cfour_module('xvmol'))
            handle.write(run_cfour_module('xvmol2ja'))
        psi4.print_out('  CFOUR scratch file %s for %s has been read\n' % ('JOBARC (binary)', item))
        c4mol = qcdb.cfour.jajo2mol(qcdb.jajo.getrec(['COORD   ', 'ATOMCHRG', 'MAP2ZMAT']))

        # S/R: Reap results from output file
        if isSowReap:
            isOk, msg, results = reap_job_validate(outdir, 'VPT2', item, linkage,
                ['CURRENT ENERGY', 'CURRENT DIPOLE', 'CURRENT GRADIENT'])
            if not isOk:
                raise ValidationError(msg)

            fje = results['CURRENT ENERGY']
            fjgrd = results['CURRENT GRADIENT']
            fjdip = [item / p4const.psi_dipmom_au2debye for item in results['CURRENT DIPOLE']]

        # C: Run the job and collect results
        else:
            psi4.IO.set_default_namespace(item)
            molecule = geometry(zmmol.create_psi4_string_from_molecule(), 'disp-' + item)
            molecule.update_geometry()
            gradient(lowername, **kwargs)

            fje = psi4.get_variable('CURRENT ENERGY')
            fjgrd = p4util.mat2arr(psi4.get_gradient())
            fjdip = [psi4.get_variable('CURRENT DIPOLE X') / p4const.psi_dipmom_au2debye,
                     psi4.get_variable('CURRENT DIPOLE Y') / p4const.psi_dipmom_au2debye,
                     psi4.get_variable('CURRENT DIPOLE Z') / p4const.psi_dipmom_au2debye]

        # Transform results into C4 orientation (defined by c4mol) & forge FJOBARC file
        fjobarc = qcdb.cfour.format_fjobarc(fje,
            *qcdb.cfour.backtransform(chgeMol=zmmol, permMol=c4mol, chgeGrad=fjgrd, chgeDip=fjdip))

    return fjobarc


def vpt2_instructions(stage, dir, zmats):
    """Stores all the instructions to the user for running
    :py:func:`~wrappers_cfour.vpt2` in sowreap mode. Depending on the
    *stage*, Pieces together instruction strings for the appropriate
    *stage* individualized by working directory *dir* and sown inputs
    *zmats* information.

    """
    stepFiles = ''
    for zm12 in sorted(zmats):
        stepFiles += """             psi4 %-27s %-27s\n""" % ('VPT2-' + zm12 + '.in', 'VPT2-' + zm12 + '.out')

    step0 = """
    The vpt2 sow/reap procedure has been selected through mode='sowreap'. This
    output file, the corresponding input file, and the data persistence file
    must not be edited by the user over the course of the sow/reap procedure.
    Throughout, psi4 can be invoked to move to the next stage of the procedure
    or to tally up the 'sown' jobs. This output file is overwritten each time
    psi4 is invoked, but all results and instructions accumulate.

    This procedure involves two stages of distributed calculations, harmonic and
    anharmonic, and a mimimum of three invokations of psi4 on the original input
    file (including the one that initially generated this text). From the input
    geometry (0), displacements are generated for which gradients are required.
    Input files for these are 'sown' in the current directory (1). Upon
    completion, their output files are 'reaped' into a harmonic force field (2).
    At displacements along the normal coordinates, further displacements are
    generated for which gradients are required. Input files for these are again
    'sown' in the current directory (3). Upon completion, their output files are
    'reaped' into an anharmonic force field (4), terminating the vpt2 procedure.
    Follow the instructions below to continue.

    (0)  Read Only
    --------------
       %s
       %s
       %s

""" % (dir + '/' + os.path.splitext(psi4.outfile_name())[0] + '.in',
       dir + '/' + psi4.outfile_name(),
       dir + '/' + os.path.splitext(psi4.outfile_name())[0] + '.shelf')
    step1 = """
    (1)  Sow
    --------
       Run all of the VPT2-000-*.in input files on any variety of computer
       architecture. The output file names must be as given below (default).

"""
    step2 = """
    (2)  Reap
    ---------
       Gather all the resulting output files in this directory along with the
       three read-only files from (0). Invoke psi4 again. The job will be
       trivial in length (unless sto-3g integrals on the molecule are costly)
       and give results for the harmonic frequency stage in this output file. It
       will also supply the next set of instructions.

             psi4 %-27s %-27s

""" % (os.path.splitext(psi4.outfile_name())[0] + '.in', psi4.outfile_name())
    step3 = """
    (3)  Sow
    --------
       Run all of the VPT2-*-*.in input files on any variety of computer
       architecture. The output file names must be as given below (default).

"""
    step4 = """
    (4)  Reap
    ---------
       Gather all the resulting output files in this directory along with the
       three read-only files from (0). Invoke psi4 again. The job will be
       trivial in length (unless sto-3g integrals on the molecule are costly)
       and give results for the harmonic and anharmonic frequency stages in this
       output file.

             psi4 %-27s %-27s

""" % (os.path.splitext(psi4.outfile_name())[0] + '.in', psi4.outfile_name())

    if stage == 'harmonic':
        instructions = step0 + step1 + stepFiles + step2
    elif stage == 'anharmonic':
        instructions = step0 + step3 + stepFiles + step4

    return instructions


def sown_jobs_status(dir, prefix, zmats, validate_func=None, linkage=None, keys=None):
    """Evaluate the output file status of jobs in *zmats* which should
    exist at *dir* + '/' + prefix + '-' + job + '.out'. Returns string with
    formatted summary of job status and boolean of whether all complete.
    Return boolean *isOk* signals whether all *zmats* have completed and,
    if *validate_func* present, are validated.

    """
    isOk = True
    msgError = ''
    instructions = '\n'
    instructions += p4util.banner(prefix + ' Status: ' + datetime.datetime.now().strftime("%Y-%m-%d %H:%M"), strNotOutfile=True)
    instructions += '\n'
    for job in sorted(zmats):
        outfile = dir + '/' + prefix + '-' + job + '.out'
        fjafile = dir + '/' + prefix + '-' + job + '.fja'
        formatArgs = [prefix + '-' + job, '', '', '', '']
        if os.path.isfile(outfile):
            with open(outfile, 'r') as handle:
                for line in handle:
                    if line.find('Buy a developer a beer!') > -1:
                        formatArgs[3] = 'Completed'
                        if reap_job_validate is not None:
                            isOkJob, msg, temp = reap_job_validate(dir, prefix, job, linkage, keys)
                            if isOkJob:
                                formatArgs[4] = '& Validated'
                            else:
                                isOk = False
                                msgError += msg
                                formatArgs[4] = 'INVALID'
                        break
                else:
                    isOk = False
                    formatArgs[2] = 'Running'
        elif os.path.isfile(fjafile):
            formatArgs[3] = 'Completed'
        else:
            isOk = False
            formatArgs[1] = 'Waiting'
        instructions += """             {0:<27} {1:^10} {2:^10} {3:^10} {4:^10}\n""".format(*formatArgs)
    instructions += '\n' + msgError + '\n\n'

    return isOk, instructions


def reap_job_validate(dir, prefix, item, linkage, keys):
    """For a given output file whose path is constructed with
    *dir* + '/' + *prefix* + '-' + *item* + '.out', tests that the file
    exists and has *prefix* RESULTS lines for each piece of information
    requested in list *keys* and that those lines correspond to the
    appropriate *linkage* and *item*. Returns *keys* along with their
    scanned values in dict *reapings*, along with error and success
    messages in *instructions* and a boolean *isOk* indicating whether
    all *keys* reaped sucessfully.

    """
    isOk = True
    instructions = ''
    reapings = {}
    outfile = dir + '/' + prefix + '-' + item + '.out'

    try:
        with open(outfile, 'r') as handle:
            for line in handle:
                if line.find(prefix + ' RESULT:') == 0:
                    sline = line.split()
                    if sline[2:7] == ['linkage', str(linkage), 'for', 'item', item]:
                        yieldsAt = line.find('yields')
                        beingAt = line.find('being')
                        if beingAt > yieldsAt > -1:
                            key = line[yieldsAt + 6:beingAt].strip()
                            val = line[beingAt + 5:].strip()
                            if key in keys:
                                reapings[key] = eval(val)
                                #psi4.print_out('  CFOUR scratch file %s for %s has been read\n' % ('JOBARC', zm12))
                        else:
                            isOk = False
                            instructions += """Outfile file %s
    has corrupted sowreap result line:\n%s\n\n""" % (outfile, line)
                    else:
                        isOk = False
                        instructions += """Outfile file %s
    has sowreap result of either incompatible linkage (observed: %s, expected: %s)
    or incompatible job affiliation (observed: %s, expected: %s).\n\n""" % \
                            (outfile, sline[3], linkage, sline[6], item)
            else:
                if len(reapings) != len(keys):
                    isOk = False
                    instructions += """Output file %s
    has missing results (observed: %s, expected: %s).\n\n""" % \
                        (outfile, reapings.keys(), keys)
    except IOError:
        isOk = False
        instructions += """Output file %s
    that was judged present and complete at the beginning of this
    job is now missing. Replace it and invoke psi4 again.\n\n""" % (outfile)

    # return file contents in instructions
    return isOk, instructions, reapings