#* **************************************************************************
#*
#* CDTK, Chemical Dynamics Toolkit
#* A modular system for chemical dynamics applications and more
#*
#* Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016
#* Oriol Vendrell, DESY, <oriol.vendrell@desy.de>
#*
#* Copyright (C) 2017, 2018, 2019
#* Ralph Welsch, DESY, <ralph.welsch@desy.de>
#*
#* Copyright (C) 2020, 2021, 2022, 2023
#* Ludger Inhester, DESY, ludger.inhester@cfel.de
#*
#* This file is part of CDTK.
#*
#* CDTK 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 3 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, see <http://www.gnu.org/licenses/>.
#*
#* **************************************************************************
#! /usr/bin/env python
import sys
import string
import os
import numpy as np
import CDTK.Tools.Conversion as conv
charge = 1
multiplicity = 1
symbols = []
input = ''
QCHEM = 'qchem'
QCHEM_INP = 'qchem.inp'
QCHEM_OUT = 'qchem.out'
[docs]
def runQchem(X_au,**args):
"""
Run Qchen for given molecular coordinates.
Input:
- X_au: 3N array, molecular coordinates in bohr
Output: if state == 0, EHF (Hartree Fock)
if state == 1, EHF - e1
...
"""
retAll = args.get('returnAll',False)
X = X_au.copy()
X = X * conv.au2an
Fileref=open(QCHEM_INP,'w')
natoms = X.size/3
X.shape = (natoms,3)
Fileref.write("$molecule\n")
Fileref.write('%i %i\n' % (charge,multiplicity) )
for i,atom in enumerate(X):
Fileref.write('%s %12.7f %12.7f %12.7f\n' % (symbols[i],atom[0],atom[1],atom[2]) )
Fileref.write("$end\n")
Fileref.write(input)
Fileref.close()
QEXEC = QCHEM + ' ' + QCHEM_INP + ' > ' + QCHEM_OUT
os.system(QEXEC)
EHF,orbE = getHF_Energy(QCHEM_OUT)
if retAll:
return EHF,orbE
if state == 0:
return EHF
if state > 0:
return EHF - float(orbE[-state])
[docs]
def getHF_Energy(QOUT):
"""
Get HF and orbital energies from a QCHEM output file
"""
SCF_Energy = ''
NOCC = ''
NOCCN = ''
OrbEi=''
StateF=open(QOUT,'r')
Line=StateF.readline()
while Line:
if Line.find('SCF energy in the final basis set = ')>=0:
StateF.readline()
SCF_Energy=SCF_Energy+Line.split()[8]
E0=float(SCF_Energy)
Line=StateF.readline()
StateF=open(QOUT,'r')
Line=StateF.readline()
while Line:
if Line.find('alpha and ')>=0:
StateF.readline()
NOCC=NOCC+string.split(Line)[2]
NOCCN=int(NOCC)
if Line.find('Final Alpha MO Eigenvalues')>=0:
OrbEi=[]
StateF.readline()
while Line.find('SCF energy in the final basis set =') < 0:
Line=StateF.readline()
if Line.find('.')>=0:
OrbEi=OrbEi+ Line.split()[1:]
Line=StateF.readline()
return (E0,OrbEi[0:NOCCN])
Line=StateF.readline()
[docs]
def getHartreeFockEnergy():
"""
Return the HF energy of a QChem calculation
Raise exception if HF energy not found
"""
energy = None
outfile=open(QOUT,'r')
line=outfile.readline()
while line:
if line.find('SCF energy in the final basis set = ')>=0:
outfile.readline()
energy=float(line.split()[8])
break
line=outfile.readline()
if energy is None:
raise ValueError
return energy