Source code for qcdb.interface_gcp

# Psi4: an open-source quantum chemistry software package
# Copyright (c) 2007-2019 The Psi4 Developers.
# The copyrights for code used from other parties are included in
# the corresponding files.
# This file is part of Psi4.
# Psi4 is free software; you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, version 3.
# Psi4 is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# GNU Lesser General Public License for more details.
# You should have received a copy of the GNU Lesser General Public License along
# with Psi4; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

"""Module with functions that interface with Grimme's GCP code."""
import os
import re
import uuid
import shutil
import socket
import subprocess

    from psi4.driver.p4util.exceptions import *
    from psi4 import core
    isP4regime = True
except ImportError:
    from .exceptions import *
    isP4regime = False
from .util import parse_dertype
from .molecule import Molecule

[docs]def run_gcp(self, func=None, dertype=None, verbose=False): # dashlvl=None, dashparam=None """Function to call Grimme's GCP program to compute an a posteriori geometrical BSSE correction to *self* for several HF, generic DFT, and specific HF-3c and PBEh-3c method/basis combinations, *func*. Returns energy if *dertype* is 0, gradient if *dertype* is 1, else tuple of energy and gradient if *dertype* unspecified. The gcp executable must be independently compiled and found in :envvar:`PATH` or :envvar:`PSIPATH`. *self* may be either a qcdb.Molecule (sensibly) or a psi4.Molecule (works b/c psi4.Molecule has been extended by this method py-side and only public interface fns used) or a string that can be instantiated into a qcdb.Molecule. """ # Create (if necessary) and update qcdb.Molecule if isinstance(self, Molecule): # called on a qcdb.Molecule pass elif isinstance(self, core.Molecule): # called on a python export of a psi4.core.Molecule (py-side through Psi4's driver) self.create_psi4_string_from_molecule() elif isinstance(self, str): # called on a string representation of a psi4.Molecule (c-side through psi4.Dispersion) self = Molecule(self) else: raise ValidationError("""Argument mol must be psi4string or qcdb.Molecule""") self.update_geometry() # # Validate arguments # dashlvl = dashlvl.lower() # dashlvl = dash_alias['-' + dashlvl][1:] if ('-' + dashlvl) in dash_alias.keys() else dashlvl # if dashlvl not in dashcoeff.keys(): # raise ValidationError("""-D correction level %s is not available. Choose among %s.""" % (dashlvl, dashcoeff.keys())) if dertype is None: derint, derdriver = -1, 'gradient' else: derint, derdriver = parse_dertype(dertype, max_derivative=1) # if func is None: # if dashparam is None: # # defunct case # raise ValidationError("""Parameters for -D correction missing. Provide a func or a dashparam kwarg.""") # else: # # case where all param read from dashparam dict (which must have all correct keys) # func = 'custom' # dashcoeff[dashlvl][func] = {} # dashparam = dict((k.lower(), v) for k, v in dashparam.iteritems()) # for key in dashcoeff[dashlvl]['b3lyp'].keys(): # if key in dashparam.keys(): # dashcoeff[dashlvl][func][key] = dashparam[key] # else: # raise ValidationError("""Parameter %s is missing from dashparam dict %s.""" % (key, dashparam)) # else: # func = func.lower() # if func not in dashcoeff[dashlvl].keys(): # raise ValidationError("""Functional %s is not available for -D level %s.""" % (func, dashlvl)) # if dashparam is None: # # (normal) case where all param taken from dashcoeff above # pass # else: # # case where items in dashparam dict can override param taken from dashcoeff above # dashparam = dict((k.lower(), v) for k, v in dashparam.iteritems()) # for key in dashcoeff[dashlvl]['b3lyp'].keys(): # if key in dashparam.keys(): # dashcoeff[dashlvl][func][key] = dashparam[key] # TODO temp until figure out paramfile allowed_funcs = ['HF/MINIS', 'DFT/MINIS', 'HF/MINIX', 'DFT/MINIX', 'HF/SV', 'DFT/SV', 'HF/def2-SV(P)', 'DFT/def2-SV(P)', 'HF/def2-SVP', 'DFT/def2-SVP', 'HF/DZP', 'DFT/DZP', 'HF/def-TZVP', 'DFT/def-TZVP', 'HF/def2-TZVP', 'DFT/def2-TZVP', 'HF/631Gd', 'DFT/631Gd', 'HF/def2-TZVP', 'DFT/def2-TZVP', 'HF/cc-pVDZ', 'DFT/cc-pVDZ', 'HF/aug-cc-pVDZ', 'DFT/aug-cc-pVDZ', 'DFT/SV(P/h,c)', 'DFT/LANL', 'DFT/pobTZVP', 'TPSS/def2-SVP', 'PW6B95/def2-SVP', # specials 'hf3c', 'pbeh3c'] allowed_funcs = [f.lower() for f in allowed_funcs] if func.lower() not in allowed_funcs: raise Dftd3Error("""bad gCP func: %s. need one of: %r""" % (func, allowed_funcs)) # Move ~/.dftd3par.<hostname> out of the way so it won't interfere defaultfile = os.path.expanduser('~') + '/.dftd3par.' + socket.gethostname() defmoved = False if os.path.isfile(defaultfile): os.rename(defaultfile, defaultfile + '_hide') defmoved = True # Find environment by merging PSIPATH and PATH environment variables lenv = { 'PATH': ':'.join([os.path.abspath(x) for x in os.environ.get('PSIPATH', '').split(':') if x != '']) + \ ':' + os.environ.get('PATH'), 'LD_LIBRARY_PATH': os.environ.get('LD_LIBRARY_PATH') } # Filter out None values as subprocess will fault on them lenv = {k: v for k, v in lenv.items() if v is not None} # Find out if running from Psi4 for scratch details and such try: import psi4 except ImportError as err: isP4regime = False else: isP4regime = True # Setup unique scratch directory and move in current_directory = os.getcwd() if isP4regime: psioh = core.IOManager.shared_object() psio = core.IO.shared_object() os.chdir(psioh.get_default_path()) gcp_tmpdir = 'psi.' + str(os.getpid()) + '.' + psio.get_default_namespace() + \ '.gcp.' + str(uuid.uuid4())[:8] else: gcp_tmpdir = os.environ['HOME'] + os.sep + 'gcp_' + str(uuid.uuid4())[:8] if os.path.exists(gcp_tmpdir) is False: os.mkdir(gcp_tmpdir) os.chdir(gcp_tmpdir) # Write gcp_parameters file that governs cp correction # paramcontents = gcp_server(func, dashlvl, 'dftd3') # paramfile1 = 'dftd3_parameters' # older patched name # with open(paramfile1, 'w') as handle: # handle.write(paramcontents) # paramfile2 = '.gcppar' # with open(paramfile2, 'w') as handle: # handle.write(paramcontents) ###Two kinds of parameter files can be read in: A short and an extended version. Both are read from ###$HOME/.gcppar.$HOSTNAME by default. If the option -local is specified the file is read in from ###the current working directory: .gcppar ###The short version reads in: basis-keywo # Write dftd3_geometry file that supplies geometry to dispersion calc numAtoms = self.natom() geom = self.save_string_xyz() reals = [] for line in geom.splitlines(): lline = line.split() if len(lline) != 4: continue if lline[0] == 'Gh': numAtoms -= 1 else: reals.append(line) geomtext = str(numAtoms) + '\n\n' for line in reals: geomtext += line.strip() + '\n' geomfile = './' with open(geomfile, 'w') as handle: handle.write(geomtext) # TODO somehow the variations on save_string_xyz and # whether natom and chgmult does or doesn't get written # have gotten all tangled. I fear this doesn't work # the same btwn libmints and qcdb or for ghosts # Call gcp program command = ['gcp', geomfile] command.extend(['-level', func]) if derint != 0: command.append('-grad') try: #print('command', command) dashout = subprocess.Popen(command, stdout=subprocess.PIPE, env=lenv) except OSError as e: raise ValidationError('Program gcp not found in path. %s' % e) out, err = dashout.communicate() # Parse output success = False for line in out.splitlines(): line = line.decode('utf-8') if re.match(' Egcp:', line): sline = line.split() dashd = float(sline[1]) if re.match(' normal termination of gCP', line): success = True if not success: os.chdir(current_directory) raise Dftd3Error("""Unsuccessful gCP run.""") # Parse grad output if derint != 0: derivfile = './gcp_gradient' dfile = open(derivfile, 'r') dashdderiv = [] for line in geom.splitlines(): lline = line.split() if len(lline) != 4: continue if lline[0] == 'Gh': dashdderiv.append([0.0, 0.0, 0.0]) else: dashdderiv.append([float(x.replace('D', 'E')) for x in dfile.readline().split()]) dfile.close() if len(dashdderiv) != self.natom(): raise ValidationError('Program gcp gradient file has %d atoms- %d expected.' % \ (len(dashdderiv), self.natom())) # Prepare results for Psi4 if isP4regime and derint != 0: core.set_variable('GCP CORRECTION ENERGY', dashd) psi_dashdderiv = core.Matrix.from_list(dashdderiv) # Print program output to file if verbose if not verbose and isP4regime: verbose = True if core.get_option('SCF', 'PRINT') >= 3 else False if verbose: text = '\n ==> GCP Output <==\n' text += out.decode('utf-8') if derint != 0: with open(derivfile, 'r') as handle: text +='D', 'E') text += '\n' if isP4regime: core.print_out(text) else: print(text) # # Clean up files and remove scratch directory # os.unlink(paramfile1) # os.unlink(paramfile2) # os.unlink(geomfile) # if derint != 0: # os.unlink(derivfile) # if defmoved is True: # os.rename(defaultfile + '_hide', defaultfile) # clean up files and remove scratch directory os.chdir('..') try: shutil.rmtree(gcp_tmpdir) except OSError as err: raise OSError('Unable to remove gcp temporary directory: {}'.format(gcp_tmpdir)) from err os.chdir(current_directory) # return -D & d(-D)/dx if derint == -1: return dashd, dashdderiv elif derint == 0: return dashd elif derint == 1: return psi_dashdderiv
try: # Attach method to libmints psi4.Molecule class core.Molecule.run_gcp = run_gcp except (NameError, AttributeError): # But don't worry if that doesn't work b/c # it'll get attached to qcdb.Molecule class pass