#!/usr/bin/env python3
'''
.. code-block::
!---------------------------------------------------------------------------!
! gaussian: interface to the Gaussian program !
! Implementations by: Pavlo O. Dral, Peikun Zheng, Yi-Fan Hou !
!---------------------------------------------------------------------------!
'''
import os, sys, subprocess, math
import numpy as np
pythonpackage = True
from mlatom import constants
from mlatom import data
from mlatom import models
from mlatom import stopper
from mlatom import simulations
from mlatom.decorators import doc_inherit
[文档]
class gaussian_methods(models.model):
'''
Gaussian interface
Arguments:
method (str): Method to use
nthreads (int): equivalent to %proc in Gaussian input file
save_files_in_current_directory (bool): whether to keep input and output files, default ``'False'``
working_directory (str): path to the directory where the program output files and other tempory files are saved, default ``'None'``
.. note::
The format of method should be the same as that in Gaussian, e.g., ``'B3LYP/6-31G*'``
'''
def __init__(self,method='B3LYP/6-31G*',nthreads=None,save_files_in_current_directory=False, working_directory=None, **kwargs):
if not "GAUSS_EXEDIR" in os.environ:
if pythonpackage: raise ValueError('enviromental variable GAUSS_EXEDIR is not set')
else: stopper.stopMLatom('enviromental variable GAUSS_EXEDIR is not set')
self.method = method
self.find_energy_to_read_in_Gaussian()
if nthreads is None:
from multiprocessing import cpu_count
self.nthreads = cpu_count()
else:
self.nthreads = nthreads
self.save_files_in_current_directory = save_files_in_current_directory
self.working_directory = working_directory
if 'writechk' in kwargs:
self.writechk = kwargs['writechk']
else:
self.writechk = False
[文档]
@doc_inherit
def predict(self,molecular_database=None,molecule=None,
calculate_energy=True,
calculate_energy_gradients=False,
calculate_hessian=False,
gaussian_keywords=None,):
'''
gaussian_keywords (some type): ``# needs to be documented``.
'''
molDB = super().predict(molecular_database=molecular_database, molecule=molecule)
method = self.method
#self.energy_to_read = energy_to_read
self.calculate_energy_gradients = calculate_energy_gradients
self.calculate_energy = calculate_energy
self.calculate_hessian = calculate_hessian
# self.writechk = False
if gaussian_keywords == None:
self.gaussian_keywords = ''
else:
self.gaussian_keywords = gaussian_keywords
if calculate_energy_gradients:
if not calculate_hessian:
if 'force'.casefold() in method.casefold():
pass
else:
method += ' force(nostep)'
if calculate_hessian:
self.writechk = True
if 'freq'.casefold() in method.casefold():
pass
else:
method += ' freq'
import tempfile, subprocess
with tempfile.TemporaryDirectory() as tmpdirname:
if self.save_files_in_current_directory: tmpdirname = '.'
if self.working_directory is not None:
tmpdirname = self.working_directory
if not os.path.exists(tmpdirname):
os.makedirs(tmpdirname)
tmpdirname = os.path.abspath(tmpdirname)
for imol in range(len(molDB.molecules)):
imolecule = molDB.molecules[imol]
# Run Gaussian job
if self.gaussian_keywords != '':
run_gaussian_job(filename='molecule'+str(imol)+'.com',molecule=imolecule,gaussian_keywords=self.gaussian_keywords,nthreads=self.nthreads,method=method,cwd=tmpdirname,writechk=self.writechk)
else:
run_gaussian_job(filename='molecule'+str(imol)+'.com',molecule=imolecule,nthreads=self.nthreads,method=method,cwd=tmpdirname,writechk=self.writechk)
# Read Gaussian output file
self.parse_gaussian_output(os.path.join(tmpdirname,'molecule'+str(imol)+'.log'),imolecule)
def parse_gaussian_output(self, filename, molecule):
assistant = Gaussian_output_reading_assistant(self.calculate_energy,self.calculate_energy_gradients,self.calculate_hessian,self.energy_to_read)
assistant.read_output(filename,molecule)
def find_energy_to_read_in_Gaussian(self):
HF = data.properties_tree_node(name='HF')
MP2 = data.properties_tree_node(name='MP2')
MP3 = data.properties_tree_node(name='MP3')
MP4D = data.properties_tree_node(name='MP4D')
MP4DQ = data.properties_tree_node(name='MP4DQ')
MP4SDQ = data.properties_tree_node(name='MP4SDQ')
MP4SDTQ = data.properties_tree_node(name='MP4SDTQ')
MP5 = data.properties_tree_node(name='MP5')
CISD = data.properties_tree_node(name='CISD')
QCISD = data.properties_tree_node(name='QCISD')
QCISD_T = data.properties_tree_node(name='QCISD(T)')
CCSD = data.properties_tree_node(name='CCSD')
CCSD_T = data.properties_tree_node(name='CCSD(T)')
children_properties = []
method = self.method
parse_method = method.split()
for each in method.split():
if '/' in each:
parse_method = each.split('/')
functional = parse_method[0]
if 'HF'.casefold() in functional.casefold():
children_properties = [HF]
elif 'B3LYP'.casefold() in functional.casefold():
children_properties = [HF]
# MP2
elif 'MP2'.casefold() in functional.casefold():
children_properties = [HF,MP2]
# MP3
elif 'MP3'.casefold() in functional.casefold():
children_properties = [HF,MP2,MP3]
# MP4(SDQ) & MP4(SDTQ) (by default latter)
elif 'MP4'.casefold() in functional.casefold():
if 'SDQ'.casefold() in functional.casefold(): # MP4(SDQ)
children_properties = [HF,MP2,MP3,MP4D,MP4DQ,MP4SDQ]
elif 'SDTQ'.casefold() in functional.casefold(): # MP4(SDTQ)
children_properties = [HF,MP2,MP3,MP4D,MP4DQ,MP4SDQ,MP4SDTQ]
else: # MP4
children_properties = [HF,MP2,MP3,MP4D,MP4DQ,MP4SDQ,MP4SDTQ]
# CISD
elif 'CISD'.casefold() in functional.casefold():
children_properties = [HF,MP2,MP3,CISD]
# QCISD & QCISD(T)
elif 'QCISD'.casefold() in functional.casefold():
if 'QCISD(T'.casefold() in functional.casefold(): # QCISD(T)
children_properties = [HF,MP2,MP3,MP4D,MP4DQ,MP4SDQ,QCISD,QCISD_T]
else: # QCISD
children_properties = [HF,MP2,MP3,MP4D,MP4DQ,MP4SDQ,QCISD]
# CCSD & CCSD(T)
elif 'CCSD'.casefold() in functional.casefold():
if 'CCSD(T'.casefold() in functional.casefold(): # CCSD(T)
children_properties = [HF,MP2,MP3,MP4D,MP4DQ,MP4SDQ,CCSD,CCSD_T]
else: # CCSD
children_properties = [HF,MP2,MP3,MP4D,MP4DQ,MP4SDQ,CCSD]
else:
children_properties = [HF]
self.energy_to_read = data.properties_tree_node(name='energy',children=children_properties)
def run_gaussian_job(**kwargs):
if 'filename' in kwargs:
filename = kwargs['filename']
if 'molecule' in kwargs: molecule = kwargs['molecule']
gaussian_keywords = ''
if 'gaussian_keywords' in kwargs:
gaussian_keywords = kwargs['gaussian_keywords']
if 'nthreads' in kwargs:
nthreads = kwargs['nthreads']
else:
nthreads = 1
memory = ''
if 'memory' in kwargs:
memory = f"%mem={kwargs['memory']}\n"
gaussian_keywords = f'{memory}%nproc={nthreads}\n' + gaussian_keywords
if 'model_predict_kwargs' in kwargs:
model_predict_kwargs_str = str(kwargs['model_predict_kwargs'])
else:
model_predict_kwargs_str = "{}"
model_predict_kwargs_str_file = 'model_predict_kwargs'
with open(model_predict_kwargs_str_file, 'w') as f:
f.write(model_predict_kwargs_str)
if 'external_task' in kwargs:
pythonbin = sys.executable
path_to_this_file=os.path.abspath(__file__)
external_task = kwargs['external_task']
if 'gaussian_keywords' in kwargs:
gaussian_keywords += "\nexternal='%s %s'\n" % (pythonbin, path_to_this_file)
elif external_task.lower() == 'opt':
gaussian_keywords += "#p opt(nomicro) external='%s %s %s'\n" % (pythonbin, path_to_this_file, model_predict_kwargs_str_file)
elif external_task.lower() == 'freq':
gaussian_keywords += "#p freq external='%s %s %s'\n" % (pythonbin, path_to_this_file, model_predict_kwargs_str_file)
elif external_task.lower() == 'freq(anharmonic)':
if 'frequency_keywords' in kwargs:
if len(kwargs['frequency_keywords']) !=0:
extra_keywords = ','.join(kwargs['frequency_keywords'])
extra_keywords = ','+extra_keywords
else:
extra_keywords = ''
else:
extra_keywords = ''
gaussian_keywords += f"#p freq(anharmonic{extra_keywords}) external='%s %s %s'\n" % (pythonbin, path_to_this_file, model_predict_kwargs_str_file)
elif external_task.lower() == 'ts':
gaussian_keywords += "#p opt(ts,calcfc,noeigen,nomicro) external='%s %s %s'\n" % (pythonbin, path_to_this_file, model_predict_kwargs_str_file)
elif external_task.lower() == 'irc':
gaussian_keywords += "#p irc(calcfc) external='%s %s %s'\n" % (pythonbin, path_to_this_file, model_predict_kwargs_str_file)
else:
if 'gaussian_keywords' in kwargs:
gaussian_keywords += '\n'
if 'method' in kwargs:
if 'freq' in kwargs['method']:
kwargs['method'] = 'p ' + kwargs['method']
gaussian_keywords += '# '+kwargs['method']+'\n'
if 'cwd' in kwargs:
cwd = kwargs['cwd']
else:
cwd='.'
if 'writechk' in kwargs:
writechk = kwargs['writechk']
else:
writechk = False
if writechk:
gaussian_keywords = f'%chk={filename[:-4]}.chk\n'+gaussian_keywords
if 'additional_input' in kwargs:
additional_input = kwargs['additional_input']
else:
additional_input = ''
write_gaussian_input_file(filename=os.path.join(cwd,filename), molecule=molecule, gaussian_keywords=gaussian_keywords,additional_input=additional_input)
Gaussianbin, _ = check_gaussian()
proc = subprocess.Popen([Gaussianbin, filename], stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=cwd, universal_newlines=True)
proc.communicate()
def check_gaussian():
status = os.popen('echo $GAUSS_EXEDIR').read().strip()
if len(status) != 0:
Gaussianroot = status.split('bsd')[0]
if 'g16' in Gaussianroot:
version = 'g16'
elif 'g09' in Gaussianroot:
version = 'g09'
Gaussianbin = os.path.join(Gaussianroot,version)
else:
stopper.stopMLatom('Cannot find Gaussian software in the environment, set $GAUSS_EXEDIR environmental variable')
version = version.replace('g', '')
return Gaussianbin, version
def write_gaussian_input_file(**kwargs):
if 'filename' in kwargs: filename = kwargs['filename']
if 'molecule' in kwargs: molecule = kwargs['molecule']
if 'gaussian_keywords' in kwargs: gaussian_keywords = kwargs['gaussian_keywords']
if 'additional_input' in kwargs:
additional_input = kwargs['additional_input']
else:
additional_input = ''
if 'comment' in molecule.__dict__:
if molecule.comment != '':
title = molecule.comment
elif molecule.id != '':
title = molecule.id
else:
title = 'Gaussian calculations from MLatom interface'
with open(filename, 'w') as f:
f.writelines(gaussian_keywords)
f.writelines(f'\n{title}\n')
f.writelines(f'\n{molecule.charge} {molecule.multiplicity}\n')
for atom in molecule.atoms:
f.writelines('%-3s %25.13f %25.13f %25.13f\n' % (atom.element_symbol,
atom.xyz_coordinates[0], atom.xyz_coordinates[1], atom.xyz_coordinates[2]))
f.writelines('\n')
f.writelines(additional_input)
def gaussian_external(EIn_file, EOu_file, model_predict_kwargs):
# write new coordinate into 'xyz_temp.dat'
derivs, molecule = read_gaussian_EIn(EIn_file)
# calculate energy, gradients, hessian for new coordinates
# import json
# with open('model.json', 'r') as fjson:
# model_dict = json.load(fjson)
# if 'method' in model_dict:
# kwargs = {}
# if 'kwargs' in model_dict:
# kwargs = model_dict['kwargs']
# del model_dict['kwargs']
# model = models.methods(**model_dict, **kwargs)
model = models.load('model.json')
calc_hessian = False
if derivs == 2: calc_hessian = True
model.predict(molecule=molecule, calculate_energy=True, calculate_energy_gradients=True, calculate_hessian=calc_hessian, **model_predict_kwargs)
if not 'energy' in molecule.__dict__:
if pythonpackage: raise ValueError('model did not return any energy')
else: stopper.stopMLatom('model did not return any energy')
write_gaussian_EOu(EOu_file, derivs, molecule)
if os.path.exists('gaussian_opttraj.json'):
opttraj = data.molecular_trajectory()
opttraj.load(filename='gaussian_opttraj.json', format='json')
nsteps = len(opttraj.steps)
opttraj.steps.append(data.molecular_trajectory_step(step=nsteps, molecule=molecule))
opttraj.dump(filename='gaussian_opttraj.json', format='json')
if os.path.exists('gaussian_freq_mol.json'):
molecule.dump(filename='gaussian_freq_mol.json', format='json')
def read_gaussian_EIn(EIn_file):
molecule = data.molecule()
with open(EIn_file, 'r') as fEIn:
lines = fEIn.readlines()
line0 = lines[0].strip().split()
natoms = int(line0[0]); derivs = int(line0[1])
molecule.charge = int(line0[2]); molecule.multiplicity = int(line0[3])
for i in range(1, natoms+1):
xx = lines[i].strip().split()
atom = data.atom(atomic_number=int(xx[0]),
xyz_coordinates=np.array(xx[1:-1]).astype('float')*constants.Bohr2Angstrom)
molecule.atoms.append(atom)
return derivs, molecule
def write_gaussian_EOu(EOu_file, derivs, molecule):
import fortranformat
with open(EOu_file, 'w') as fEOu:
# energy, dipole-moment (xyz) E, Dip(I), I=1,3
if 'dipole_moment' in molecule.__dict__.keys():
dp = molecule.dipole_moment
else:
dp = [0.0,0.0,0.0]
writer = fortranformat.FortranRecordWriter('(4D20.12)')
output = writer.write([molecule.energy, dp[0], dp[1], dp[2]])
fEOu.write(output)
fEOu.write('\n')
writer = fortranformat.FortranRecordWriter('(3D20.12)')
# gradient on atom (xyz) FX(J,I), J=1,3; I=1,NAtoms
output = writer.write(molecule.get_energy_gradients().flatten()*constants.Bohr2Angstrom)
fEOu.write(output)
fEOu.write('\n')
if derivs == 2:
natoms = len(molecule.atoms)
# polarizability Polar(I), I=1,6
polor = np.zeros(6)
output = writer.write(polor)
fEOu.write(output)
fEOu.write('\n')
# dipole derivatives DDip(I), I=1,9*NAtoms
ddip = np.zeros(9*natoms)
if 'dipole_derivatives' in molecule.__dict__.keys():
# Rescale the dipole derivatives by a factor of 0.1 so that anharmonic
# infrared intensities can be printed properly when using AIQM1
# Later the intensities will be rescaled by a factor of 100
ddip = molecule.dipole_derivatives * constants.Bohr2Angstrom * constants.Debye2au / 10.0
output = writer.write(ddip)
fEOu.write(output)
fEOu.write('\n')
# force constants FFX(I), I=1,(3*NAtoms*(3*NAtoms+1))/2
output = writer.write(molecule.hessian[np.tril_indices(natoms*3)]*constants.Bohr2Angstrom**2)
fEOu.write(output)
def read_freq_thermochemistry_from_Gaussian_output(outputfile, molecule):
successful = False
if os.path.exists('gaussian_freq_mol.json'):
molecule.load(filename='gaussian_freq_mol.json', format='json')
lines = open(outputfile, 'r').readlines()
natoms = len(molecule.atoms)
molecule.frequencies = []
molecule.force_constants = []
molecule.reduced_masses = []
molecule.infrared_intensities = []
for atom in molecule.atoms:
atom.normal_modes = []
# Ugly fix: read from last job step
termination_list = [0]
for iline in range(len(lines)):
if 'termination' in lines[iline]:
termination_list.append(iline+1)
istart = termination_list[-2]
for iline in range(istart,len(lines)):
if 'Frequencies --' in lines[iline]: molecule.frequencies += [float(xx) for xx in lines[iline].split()[2:]]
if 'Red. masses --' in lines[iline]: molecule.reduced_masses += [float(xx) for xx in lines[iline].split()[3:]]
if 'Frc consts --' in lines[iline]: molecule.force_constants += [float(xx) for xx in lines[iline].split()[3:]]
if 'IR Inten --' in lines[iline]:
if 'dipole_derivatives' in molecule.__dict__.keys(): molecule.infrared_intensities += [float(xx)*100 for xx in lines[iline].split()[3:]]
if 'Atom AN X Y Z' in lines[iline]:
for iatom in range(natoms):
xyzs = np.array([float(xx) for xx in lines[iline+iatom+1].split()[2:]])
nxyzs = len(xyzs) / 3
for xyz in np.array_split(xyzs, nxyzs):
molecule.atoms[iatom].normal_modes.append(list(xyz))
if 'Zero-point correction=' in lines[iline]: molecule.ZPE = float(lines[iline].split()[2])
if 'Thermal correction to Energy=' in lines[iline]: molecule.DeltaE2U = float(lines[iline].split()[-1])
if 'Thermal correction to Enthalpy=' in lines[iline]: molecule.DeltaE2H = float(lines[iline].split()[-1])
if 'Thermal correction to Gibbs Free Energy=' in lines[iline]: molecule.DeltaE2G = float(lines[iline].split()[-1])
if 'Sum of electronic and zero-point Energies=' in lines[iline]:
molecule.U0 = float(lines[iline].split()[-1])
molecule.H0 = molecule.U0
if 'Sum of electronic and thermal Energies=' in lines[iline]: molecule.U = float(lines[iline].split()[-1])
if 'Sum of electronic and thermal Enthalpies=' in lines[iline]: molecule.H = float(lines[iline].split()[-1])
if 'Sum of electronic and thermal Free Energies=' in lines[iline]: molecule.G = float(lines[iline].split()[-1])
if 'E (Thermal) CV S' in lines[iline]:
molecule.S = float(lines[iline+2].split()[-1]) * constants.kcalpermol2Hartree / 1000.0
if 'Normal termination of Gaussian' in lines[iline]:
successful = True
for atom in molecule.atoms:
atom.normal_modes = np.array(atom.normal_modes)
return successful
def read_anharmonic_frequencies(outputfile,molecule):
freq_len = len(molecule.frequencies)
anharmonic_frequencies = []
harmonic_frequencies = []
anharmonic_overtones = []
harmonic_overtones = []
anharmonic_combination_bands = []
harmonic_combination_bands = []
fundamental_bands_flag = True
overtones_flag = True
combination_bands_flag = True
anharmonic_infrared_intensities = []
anharmonic_overtones_infrared_intensities = []
anharmonic_combination_bands_infrared_intensities = []
IRflags = []
with open(outputfile,'r') as f:
lines = f.readlines()
for iline in range(len(lines)):
if 'Fundamental Bands' in lines[iline] and fundamental_bands_flag:
flag = iline+3
for ii in range(len(molecule.frequencies)):
templine = lines[flag+ii]
anharmonic_frequencies.append(eval(templine[38:48]))
harmonic_frequencies.append(eval(templine[24:38]))
# temp = lines[flag+ii].strip().split()
# anharmonic_frequencies.append(eval(temp[-4]))
# harmonic_frequencies.append(eval(temp[-5]))
fundamental_bands_flag = False
if 'Overtones' in lines[iline] and overtones_flag:
flag = iline+3
read_overtone = True
icount = 0
while read_overtone:
templine = lines[flag+icount]
temp = lines[flag+icount].strip().split()
if temp == []:
read_overtone = False
else:
anharmonic_overtones.append(eval(templine[38:48]))
harmonic_overtones.append(eval(templine[24:38]))
icount += 1
overtones_flag = False
if 'Combination Bands' in lines[iline] and combination_bands_flag:
flag = iline+3
read_combinaion_band = True
icount = 0
while read_combinaion_band:
templine = lines[flag+icount]
temp = lines[flag+icount].strip().split()
if temp == []:
read_combinaion_band = False
else:
anharmonic_combination_bands.append(eval(templine[38:48]))
harmonic_combination_bands.append(eval(templine[24:38]))
icount += 1
combination_bands_flag = False
if 'ZPE(anh)' in lines[iline]:
# molecule.anharmonic_ZPE = eval(lines[iline].strip().split()[-2].replace('D','E')) * constants.kJpermol2Hartree
molecule.anharmonic_ZPE = eval(lines[iline][48:61].replace('D','E')) * constants.kJpermol2Hartree
if 'Input values of T(K) and P(atm):' in lines[iline]:
temperature = eval(lines[iline].strip().split()[-2])
E = molecule.U - molecule.DeltaE2U
molecule.anharmonic_U0 = E + molecule.anharmonic_ZPE
molecule.anharmonic_H0 = E + molecule.anharmonic_ZPE
molecule.anharmonic_DeltaE2U = eval(lines[iline+4].strip().split()[-2].replace('D','E')) * constants.kJpermol2Hartree
molecule.anharmonic_U = E + molecule.anharmonic_DeltaE2U
molecule.anharmonic_DeltaE2H = eval(lines[iline+5].strip().split()[-2].replace('D','E')) * constants.kJpermol2Hartree
molecule.anharmonic_H = E + molecule.anharmonic_DeltaE2H
molecule.anharmonic_S = eval(lines[iline+6].strip().split()[-3].replace('D','E')) * constants.kJpermol2Hartree / 1000.0
molecule.anharmonic_G = molecule.anharmonic_H - temperature * molecule.anharmonic_S
molecule.anharmonic_DeltaE2G = molecule.anharmonic_G - E
# Read infrared intensities if found
if 'I(anharm)' in lines[iline]:
if 'dipole_derivatives' in molecule.__dict__.keys():
IRflags.append(iline+1)
# Read infrared intensities if found
if len(IRflags) == 3:
# Read fundamental bands intensities
flag = IRflags[0]
for ii in range(len(molecule.frequencies)):
templine = lines[flag+ii]
try:
intensity = eval(templine.strip().split()[-1]) * 100 # the intensities are rescaled by a factor of 100
anharmonic_infrared_intensities.append(intensity)
except:
anharmonic_infrared_intensities.append(math.nan)
# Read overtones intensities
flag = IRflags[1]
read_overtone = True
icount = 0
while read_overtone:
templine = lines[flag+icount]
temp = lines[flag+icount].strip().split()
if temp == []:
read_overtone = False
else:
try:
intensity = eval(temp[-1]) * 100 # the intensities are rescaled by a factor of 100
anharmonic_overtones_infrared_intensities.append(intensity)
except:
anharmonic_overtones_infrared_intensities.append(math.nan)
icount += 1
# Read combination bands intensities
flag = IRflags[2]
read_combinaion_band = True
icount = 0
while read_combinaion_band:
templine = lines[flag+icount]
temp = lines[flag+icount].strip().split()
if temp == []:
read_combinaion_band = False
else:
try:
intensity = eval(temp[-1]) * 100 # the intensities are rescaled by a factor of 100
anharmonic_combination_bands_infrared_intensities.append(intensity)
except:
anharmonic_combination_bands_infrared_intensities.append(math.nan)
icount += 1
molecule.anharmonic_frequencies = []
order = np.argsort(harmonic_frequencies)
for each in order:
molecule.anharmonic_frequencies.append(anharmonic_frequencies[each])
if anharmonic_infrared_intensities != []:
molecule.anharmonic_infrared_intensities = []
for each in order:
molecule.anharmonic_infrared_intensities.append(anharmonic_infrared_intensities[each])
molecule.anharmonic_overtones_infrared_intensities = anharmonic_overtones_infrared_intensities
molecule.anharmonic_combination_bands_infrared_intensities = anharmonic_combination_bands_infrared_intensities
molecule.anharmonic_overtones = anharmonic_overtones
molecule.harmonic_overtones = harmonic_overtones
molecule.anharmonic_combination_bands = anharmonic_combination_bands
molecule.harmonic_conmination_bands = harmonic_combination_bands
def check_Gaussian_job_status(lines,flag,mol):
if flag == -1:
pass
mol.Gaussian_job_status = False
# for line in lines:
# if 'Normal termination of Gaussian' in line:
# successful = True
if 'Normal termination of Gaussian' in lines[flag]:
mol.Gaussian_job_status = True
#return successful
def read_energy_gradients_from_Gaussian_output(lines,flag,mol):
if flag == -1:
return
Natoms = len(mol.atoms)
for iatom in range(Natoms):
force = np.array([eval(each) for each in lines[flag+iatom].strip().split()[2:]])
mol.atoms[iatom].energy_gradients = - force / constants.Bohr2Angstrom
def read_mulliken_charges_from_Gaussian_output(lines,flag,mol):
if flag == -1:
return
Natoms = len(mol.atoms)
for iatom in range(Natoms):
charge = eval(lines[flag+iatom].strip().split()[2])
mol.atoms[iatom].mulliken_charge = charge
def read_dipole_moment_from_Gaussian_output(lines,flag,mol):
if flag == -1:
return
try:
raw = lines[flag].strip().split()
dipole = np.array([eval(raw[1]),eval(raw[3]),eval(raw[5]),eval(raw[7])])
mol.dipole_moment = dipole
except:
return
def read_electronic_energy_from_Gaussian_output(lines,flag,mol,energy_to_read):
if flag == -1:
return
tmpline = lines[flag]
archive = ''
while(tmpline!='\n'):
archive += tmpline[1:-1]
flag += 1
tmpline = lines[flag]
archive_split = archive.split('\\') # split with \
for ii in range(len(archive_split)):
for each_children in energy_to_read.children:
if each_children.name.casefold()+'=' in archive_split[ii].casefold():
energy = eval(archive_split[ii].split('=')[-1])
mol.__dict__[each_children.name] = energy
mol.energy = mol.__dict__[energy_to_read.children[-1].name]
def read_Hessian_matrix_from_Gaussian_output(filename,molecule):
natoms = len(molecule.atoms)
hessian = np.zeros((3*natoms,3*natoms))
lines = open(filename,'r').readlines()
for iline in range(len(lines)):
if 'The second derivative matrix:' in lines[iline]:
flag = iline + 1
nblocks = int((3*natoms-0.5)//5 + 1)
icount = 0
for iblock in range(nblocks):
for ii in range(3*natoms-5*iblock):
temp = [eval(each) for each in lines[flag+icount+ii+1].strip().split()[1:]]
for jj in range(min(5,3*natoms-5*iblock)):
if ii >= jj:
hessian[ii+5*iblock,jj+5*iblock] = temp[jj]
icount += 3*natoms-5*iblock+1
for ii in range(3*natoms):
for jj in range(ii+1,3*natoms):
hessian[ii,jj] = hessian[jj,ii]
hessian /= constants.Bohr2Angstrom**2
molecule.hessian = hessian
def read_xyz_coordinates_from_Gaussian_output(filename:str):
lines = open(filename,'r').readlines()
# Choose the last standard orientation
for iline in range(len(lines)):
if 'Standard orientation:' in lines[iline]:
flag = iline+5
xyz_string = ''
icount = 0
while True:
xyz = lines[flag+icount].strip().split()
if len(xyz) > 1:
xyz_string += f'{xyz[1]} {xyz[3]} {xyz[4]} {xyz[5]}\n'
else:
break
icount += 1
natoms = icount + 1
xyz_string = f'{natoms}\n\n' + xyz_string[:-1]
return data.molecule.from_xyz_string(string=xyz_string)
def read_Hessian_matrix_from_Gaussian_chkfile(filename,molecule):
natoms = len(molecule.atoms)
hessian = np.zeros((3*natoms,3*natoms))
cmdargs = ['chkchk','-p',filename]
cmd = f'chkchk -p {filename}'
proc = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE, universal_newlines=True, shell=True)
outs,errs = proc.communicate()
stdout = outs.split('\n')
stderr = errs.split('\n')
outputs = []
for readable in stdout:
outputs.append(readable)
flag = None
for iline in range(len(outputs)):
if 'Input for Opt=FCCards' in outputs[iline]:
flag = iline+1
# If there is only 1 atom, chk does not have hessian matrix
if flag is None:
molecule.hessian = hessian
else:
flag += 1 + int((3*natoms-0.5)//6 + 1)# Skip energy and forces
hessian_len = int(((3*natoms)*(3*natoms+1)//2-0.5)//6 + 1)
hessian_array = []
for ii in range(hessian_len):
hessian_array += [eval(each) for each in outputs[flag+ii].strip().split()]
#print(hessian_array)
for ii in range(3*natoms):
hessian[:ii+1,ii] = hessian_array[(ii+1)*ii//2:(ii+2)*(ii+1)//2]
for ii in range(3*natoms):
for jj in range(ii+1,3*natoms):
hessian[jj,ii] = hessian[ii,jj]
#print(hessian[:,-2])
hessian /= constants.Bohr2Angstrom**2
molecule.hessian = hessian
class Gaussian_output_reading_assistant():
def __init__(self,calculate_energy=True,calculate_energy_gradients=False,calculate_hessian=False,
energy_to_read=None):
# Some arguments of Gaussian job
self.calculate_energy = calculate_energy
self.calculate_energy_gradients = calculate_energy_gradients
self.calculate_hessian = calculate_hessian
self.energy_to_read = energy_to_read
# What to read in Gaussian output file
self.read_energy = False
self.read_energy_gradients = False
self.read_hessian = False
self.read_Mulliken_charges = True
self.read_dipole_moment = True
self.read_job_status = True
if self.calculate_energy:
self.read_energy = True
if self.calculate_energy_gradients:
self.read_energy_gradients = True
if self.calculate_hessian:
self.read_hessian = True
def read_output(self,filename,molecule):
job_status_flag = -1
mulliken_charge_flag = -1
dipole_moment_flag = -1
forces_flag = -1
energy_flag = -1
with open(filename,'r') as f:
lines = f.readlines()
for iline in range(len(lines)):
if self.read_job_status:
if 'termination' in lines[iline]:
job_status_flag = iline
check_Gaussian_job_status(lines,job_status_flag,molecule)
self.read_job_status = False
if self.read_energy:
if '1\\1\\' in lines[iline]:
energy_flag = iline
read_electronic_energy_from_Gaussian_output(lines,energy_flag,molecule,energy_to_read=self.energy_to_read)
self.read_energy = False
if self.read_Mulliken_charges:
if 'Mulliken charges' in lines[iline]:
mulliken_charge_flag = iline+2
read_mulliken_charges_from_Gaussian_output(lines,mulliken_charge_flag,molecule)
self.read_Mulliken_charges = False
if self.read_dipole_moment:
if 'Dipole moment' in lines[iline]:
dipole_moment_flag = iline+1
read_dipole_moment_from_Gaussian_output(lines,dipole_moment_flag,molecule)
self.read_dipole_moment = False
if self.calculate_energy_gradients:
if self.read_energy_gradients:
if 'Forces (Hartrees/Bohr)' in lines[iline]:
forces_flag = iline+3
read_energy_gradients_from_Gaussian_output(lines,forces_flag,molecule)
if self.calculate_hessian:
newmolecule = molecule.copy()
freq_flag = read_freq_thermochemistry_from_Gaussian_output(filename,newmolecule)
#read_Hessian_matrix_from_Gaussian_output(filename,molecule)
if os.path.exists(filename[:-4]+'.chk'):
read_Hessian_matrix_from_Gaussian_chkfile(filename[:-4]+'.chk',molecule)
# !!!
# need to be fixed later in method **read_freq_thermochemistry_from_Gaussian_output** above
# now just move properties of new mol to old mol
# !!!
if freq_flag:
molecule.frequencies = newmolecule.frequencies
molecule.force_constants = newmolecule.force_constants
molecule.reduced_masses = newmolecule.reduced_masses
for ii in range(len(newmolecule.atoms)):
molecule.atoms[ii].normal_modes = newmolecule.atoms[ii].normal_modes
molecule.ZPE = newmolecule.ZPE
molecule.DeltaE2U = newmolecule.DeltaE2U
molecule.DeltaE2H = newmolecule.DeltaE2H
molecule.DeltaE2G = newmolecule.DeltaE2G
molecule.U0 = newmolecule.U0
molecule.H0 = newmolecule.H0
molecule.U = newmolecule.U
molecule.H = newmolecule.H
molecule.G = newmolecule.G
molecule.S = newmolecule.S
if __name__ == '__main__':
model_predict_kwargs_str_file, _, EIn_file, EOu_file, _, _, _ = sys.argv[1:]
with open(model_predict_kwargs_str_file) as f:
model_predict_kwargs = eval(f.read())
gaussian_external(EIn_file, EOu_file, model_predict_kwargs)