#* **************************************************************************
#*
#* 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 string
import os
import pickle
import numpy as np
import CDTK.Tools.Conversion as conv
from functools import reduce
charge = 1
multiplicity = 1
symbols = []
input = ''
state = 0
QCHEM = 'qchem'
QCHEM_INP = 'qchem.inp'
QCHEM_OUT = 'qchem.out'
QCHEM_SCR = 'save'
[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)
retstates = args.get('returnStates',[])
X = X_au.copy()
X = X * conv.au2an
try:
Fileref=open(QCHEM_INP,'w')
except:
msg = 'File '+QCHEM_INP+' could not accessed for reading'
raise ValueError(msg)
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 + ' ' + QCHEM_SCR
status = os.system(QEXEC)
EHF,orbE = getHF_Energy(QCHEM_OUT)
if EHF is None:
return None
if retstates:
Epes = []
for s in retstates:
pes = EHF - float(orbE[-s-1])
Epes.append(pes)
return np.array(Epes)
elif retAll:
Epes = []
idx = np.arange(len(orbE),dtype=int)+1
for i in idx:
pes = EHF - float(orbE[-i])
Epes.append(pes)
return np.array(Epes)
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()
try:
StateF=open(QOUT,'r')
except:
msg = 'File '+QOUT+' could not be accessed for reading'
raise ValueError(msg)
Line=StateF.readline()
while Line:
if Line.find('Convergence failure') >= 0:
return None,None
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 runQchemNMode(X_au,**args):
"""
Run Qchem for given molecular coordinates.
Input:
- X_au: 3N array, molecular coordinates in bohr
Output:
- NMode: [E, ReducedMass, NCoor]
"""
retAll = args.get('returnAll',False)
retstates = args.get('returnStates',[])
X = X_au.copy()
X = X * conv.au2an
try:
Fileref=open(QCHEM_INP,'w')
except:
msg = 'File '+QCHEM_INP+' could not accessed for reading'
raise ValueError(msg)
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("\n")
Fileref.write(input)
Fileref.close()
QEXEC = QCHEM + ' ' + QCHEM_INP + ' ' + QCHEM_OUT + ' ' + QCHEM_SCR
status = os.system(QEXEC)
NMode = get_NMode(QCHEM_OUT, natoms)
return NMode
[docs]
def get_NMode(QOUT, natoms):
"""
Get normal modes from a QCHEM output file
Input
natoms: number of atoms
Output
NMode: [E, ReducedMass, NCoor]
"""
E = [] # energies of the normal modes
ReducedMass = [] # reduced masses of the normal modes
NCoord = [] # normal coordinates of the normal modes
StateF=open(QOUT,'r')
Line=StateF.readline()
while Line:
Line=StateF.readline()
if Line.find('INFRARED INTENSITIES') >= 0:
for idum in range(4):
Line = StateF.readline()
# normal modes section
while Line.find('STANDARD THERMODYNAMIC QUANTITIES') < 0:
Line = StateF.readline()
if Line.find('Frequency') >= 0:
E.append(Line.split()[1:]) # energy
Line = StateF.readline()
Line = StateF.readline() # reduced mass
ReducedMass.append(Line.split()[2:])
for idum in range(4):
Line = StateF.readline()
# normal coordinates
PreNCoord = [[], [], []]
for iatom in range(natoms):
Line = StateF.readline()
NNMode = int((len(Line.split()) - 1) / 3)
if NNMode > 0:
PreNCoord[0].append(Line.split()[1:4])
if NNMode > 1:
PreNCoord[1].append(Line.split()[4:7])
if NNMode > 2:
PreNCoord[2].append(Line.split()[7:10])
NCoord.append(PreNCoord)
# flatten
E = flatten_list(E)
ReducedMass = flatten_list(ReducedMass)
NCoord = flatten_list(NCoord)
NMode = []
for imode, e_s in enumerate(E):
# float
e_s = float(e_s) / conv.au2ic
ReducedMass[imode] = float(ReducedMass[imode]) * conv.am2au
#
NMode.append([e_s, ReducedMass[imode], flatten_float_list(NCoord[imode])])
return NMode
[docs]
def flatten_list(list):
"""
Return flattened list
"""
flattened_list = reduce(lambda x,y: x+y, list)
return flattened_list
[docs]
def flatten_float_list(list):
"""
Return flattened list of floating numbers
"""
flattened_list = (reduce(lambda x,y: x+y, list))
for i, _l_s in enumerate(flattened_list):
flattened_list[i] = float(_l_s)
return flattened_list
[docs]
def save_qchem(a_fname, ob):
"""
Save [ob] to file
"""
pickle.dump(ob, open(a_fname, 'wb'))
[docs]
def load_qchem(a_fname):
"""
Return loaded [ob]
"""
ob_load = pickle.load(open(a_fname))
return ob_load