#* **************************************************************************
#*
#* 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/>.
#*
#* **************************************************************************
"""Provides functions to access gaussian log files and retrieve information:
-- getNormalMode
-- getAtomNumber
-- getMasses
-- getCoordinates
-- tester
"""
import os
import sys
[docs]
def getNormalMode(f,i):
"""Get the coordinates of a mode from a gaussian log file
f -- full path to the gaussian frequency file
i -- index of the normal mode to retrieve
"""
mode=[]
log = file(f,'r')
n = getAtomNumber(f)
if hasHighPrecisionModes(f):
col = (i-1)%5 + 1
ind = col + 2
row = (i+4)/5
j=0
while 1:
line=log.readline()
if not line: break
if line.find(" Frequencies --") != -1:
j=j+1
if j == row:
while 1:
line=log.readline()
if line.find("Coord") != -1:
break
for nn in range(3*n):
lvec = log.readline().split()
mode.append(lvec[ind])
else:
col = (i-1)%3 + 1
ind = 3*(col-1) + 2
row = (i+2)/3
j=0
while 1:
line=log.readline()
if not line: break
if line.find(" Frequencies --") != -1:
j=j+1
if j == row:
log.readline()
log.readline()
log.readline()
log.readline()
for nn in range(n):
lvec=log.readline().split()
mode= mode + lvec[ind:ind+3]
log.close()
mode = map(float,mode)
return mode
[docs]
def getNormalModeFreq(f,i):
"""
Get normal mode frequency of mode i from gaussian log file
f -- full path to the gaussian frequency file
i -- index of the normal mode to retrieve
"""
log = file(f,'r')
freqs = []
while 1:
line = log.readline()
if not line: break
if line.find(' Frequencies -- ') == 0:
frline = line.split()[2:]
for fr in frline:
freqs.append(fr)
return float(freqs[i-1])
[docs]
def getAtomNumber(f):
"""Get the number of atoms of a calculation stored in a log file
f -- full path to the gaussian file
"""
num = 0
log = file(f,'r')
while 1:
line=log.readline()
if not line: break
if line.find(" Number Number") != -1:
line=log.readline()
while 1:
line=log.readline()
if not line: break
if line.find(" ---") != -1: break
else:
lvec = line.split()
num = lvec[0]
log.close()
return int(num)
log.close()
sys.exit("Error in getAtomNumber, no number of atoms found")
[docs]
def getMasses(f):
"""Get the masses from a gaussian freq file.
In this implementation the "Thermochemistry" section is used.
f -- full path to the gaussian file
"""
masses=[]
log=file(f,'r')
for line in log:
if line.find("and mass") != -1:
lvec=line.split()
masses.append(lvec[8])
masses.append(lvec[8])
masses.append(lvec[8])
if line.find("Molecular mass:") != -1:
break
masses = map(float,masses)
return masses
[docs]
def getCoordinates(f):
"""Get the standard orientation coordinates from a
gaussian freq file
"""
coords=[]
log=file(f,'r')
while 1:
line=log.readline()
if not line: break
if line.find("Standard orientation:") != -1:
line=log.readline()
line=log.readline()
line=log.readline()
line=log.readline()
while 1:
line=log.readline()
if line.find(" ---") != -1: break
lvec=line.split()
coords.append(lvec[3])
coords.append(lvec[4])
coords.append(lvec[5])
break
coords = map(float,coords)
return coords
[docs]
def getSymbols(f):
"""Get the list of atomic symbols from a gaussian log file
Their are taken from the mulliken population analysis,
which is almost always there
"""
s=[]
log=file(f,'r')
while 1:
line=log.readline()
if not line: break
if line.find("Mulliken atomic charges") != -1:
line=log.readline()
while 1:
line=log.readline()
if line.find("Sum") != -1: break
lvec = line.split()
s.append(lvec[1])
break
return s
[docs]
def hasHighPrecisionModes(f):
"""Return True if the Gaussian file contains high precision modes"""
log=file(f,'r')
for l in log:
low = l.lower()
if (low.find('freq') != -1) and (low.find('#') == 1):
if low.find('hpmodes') != -1:
log.close()
return True
else:
log.close()
return False
log.close()
return None
[docs]
def getSCFEnergy(f):
"""
Get last SCF electronic energy in log file f
"""
log=file(f,'r')
ener = None
for line in log:
if (line.find('SCF Done') == 1):
ener = float(line.split()[4])
return ener
[docs]
def getBOMD(f):
"""
Get trajectory from a gaussian BOMD file
Return a list of list, one for each geometry
"""
au2an = 0.529177249
natoms = len(getSymbols(f))
log=file(f,'r')
traj = []
while 1:
line=log.readline()
if not line: break
if line.find(' Summary information for step') == 0:
geom = []
for i in range(8):
line=log.readline()
for i in range(natoms):
atomxyz = []
sline=log.readline().split()
atomxyz.append(float('e'.join(sline[3].split('D')))*au2an)
atomxyz.append(float('e'.join(sline[5].split('D')))*au2an)
atomxyz.append(float('e'.join(sline[7].split('D')))*au2an)
geom.append(atomxyz)
traj.append(geom)
return traj
[docs]
def g03toVTF(f,vtf=None):
"""
Generate a trajectory in .vtf format from a gaussian log file
"""
if not vtf:
out = sys.stdout
else:
out = file(vtf,'w')
symb = getSymbols(f)
traj = getBOMD(f)
for i,s in enumerate(symb):
out.write( 'atom %i radius 1.0 name %s\n' % (i,s) )
out.write('\n')
for geom in traj:
out.write('timestep\n')
for atom in geom:
out.write(' %12.7f %12.7f %12.7f\n' % (atom[0],atom[1],atom[2]) )
out.write('\n')
if vtf:
out.close()
[docs]
def tester():
"""Tests the module
"""
f1="/home/oriol/Code/python/packages/testfiles/gaussian.2.test"
print("atom numbers")
print(getAtomNumber(f1))
print("mode 8")
print(getNormalMode(f1,8))
print("masses")
print(getMasses(f1))
print("coordinates")
print(getCoordinates(f1))
print("symbols")
print(getSymbols(f1))
if __name__=='__main__':
tester()