diff --git a/pygsm/level_of_theories/qchem.py b/pygsm/level_of_theories/qchem.py index bb3d050..cc0e788 100644 --- a/pygsm/level_of_theories/qchem.py +++ b/pygsm/level_of_theories/qchem.py @@ -55,6 +55,25 @@ def write_preamble(self,geom,multiplicity,tempfilename,jobtype='FORCE'): tempfile.write(line) tempfile.write('{} {}\n'.format(self.charge,multiplicity)) + # QMMM + if os.path.isfile("MM_region.txt"): + with open("MM_region.txt") as MM_region: + MM_geom = MM_region.readlines() + if not os.path.isfile("link.txt"): + raise FileNotFoundError('QMMM calculation needs link.txt file') + else: + with open("link.txt") as link: + link_lines = link.readlines() + # calling copy to avoid modifing the geomtry + from copy import deepcopy + geom = deepcopy(geom) + nHcap = len(geom) + len(MM_geom) - len(link_lines) + geom = deepcopy(geom[:-nHcap]) + for line in MM_geom: + line = line.split() + geom.append(line) + # the general format for the $molecule section in QMMM is + # if os.path.isfile("link.txt"): with open("link.txt") as link: link_lines = link.readlines() @@ -111,17 +130,18 @@ def run(self,geom,multiplicity): temp += 1 with open(efilepath) as efile: - gradlines = efile.readlines() - temp = 0 - tmp=[] - for lines in gradlines: - if '$' in lines: - temp+=1 - elif temp == 2: - tmpline = lines.split() - tmp.append([float(i) for i in tmpline]) - elif temp == 3: - break + tmp=[] + line = efile.readline() + while line != '': + if line.startswith('$gradient'): + line = efile.readline() + natoms = len(self.geom) + for i in range(natoms): + tmpline = line.split() + tmp.append([float(i) for i in tmpline]) + line = efile.readline() + break + line = efile.readline() self.grada.append((multiplicity,tmp)) else: raise NotImplementedError @@ -179,3 +199,4 @@ def copy(cls,lot,options,copy_wavefunction=True): os.system(cmd) return cls(lot.options.copy().set_values(options)) + diff --git a/scripts/f2s.py b/scripts/f2s.py new file mode 100644 index 0000000..b857e67 --- /dev/null +++ b/scripts/f2s.py @@ -0,0 +1,191 @@ +# -*- coding: utf-8 -*- +""" +Create QM_region.xyz, MM_region.txt, link.txt, and frozen.txt for GSM job in QMMM system from a given freq output file. + +python f2s.py qchem_freq.out +""" + +import argparse + +class QChemLog(): + """ + Represent an output file from QChem. The attribute `path` refers to the + location on disk of the QChem output file of interest. + """ + def __init__(self, path): + self.path = path + self.check() + + def check(self): + is_freq = False + if not self.is_QM_MM_INTERFACE: + raise ValueError('This file is not for QMMM system') + with open(self.path) as f: + line = f.readline() + while line != '': + if 'VIBRATIONAL ANALYSIS' in line: + is_freq = True + line = f.readline() + if not is_freq: + raise ValueError('This file is not freq job output file') + + def get_QM_ATOMS(self): + """ + Return the index of QM_atoms. + """ + QM_atoms = [] + + with open(self.path, 'r') as f: + line = f.readline() + while line != '': + if '$QM_ATOMS' in line.upper(): + line = f.readline() + while '$end' not in line and line != '\n': + QM_atoms.append(line.strip()) + line = f.readline() + break + line = f.readline() + + return QM_atoms + + def is_QM_MM_INTERFACE(self): + """ + Return the bool value. + """ + is_QM_MM_INTERFACE = False + + with open(self.path, 'r') as f: + line = f.readline() + while line != '': + if 'QM_MM_INTERFACE' in line.upper(): + is_QM_MM_INTERFACE = True + line = f.readline() + + return is_QM_MM_INTERFACE + + def get_USER_CONNECT(self): + """ + Return a list of string with " " . + """ + USER_CONNECTS = [] + + with open(self.path, 'r') as f: + line = f.readline() + while line != '': + if '$MOLECULE' in line.upper(): + + for i in range(2): + line = f.readline() + while '$end' not in line and line != '\n': + _, _, _, _, MM_atom_type, Bond_1, Bond_2, Bond_3, Bond_4 = line.split() + USER_CONNECT = '{} {} {} {} {}\n'.format(MM_atom_type, Bond_1, Bond_2, Bond_3, Bond_4) + USER_CONNECTS.append(USER_CONNECT) + line = f.readline() + break + line = f.readline() + + return USER_CONNECTS + + def get_fixed_molecule(self): + """ + Return the string of fixed part in $moledule + """ + # + # O 7.256000 1.298000 9.826000 + # O 6.404000 1.114000 12.310000 + # O 4.077000 1.069000 0.082000 + # H 1.825000 1.405000 12.197000 + # H 2.151000 1.129000 9.563000 + # ----------------------------------- + + fixed_molecule_string = '' + + n_atoms = len(self.get_QM_ATOMS()) + with open(self.path, 'r') as f: + line = f.readline() + while line != '': + if '$MOLECULE' in line.upper(): + for i in range(2): + line = f.readline() + for i in range(n_atoms): + line = f.readline() + while '$end' not in line and line != '\n': + a, x, y, z, _, _, _, _, _ = line.split() + line = '{:<2s} {:.10f} {:.10f} {:.10f}\n'.format(a, float(x), float(y), float(z)) + fixed_molecule_string += line + line = f.readline() + break + line = f.readline() + + return fixed_molecule_string[:-1] + + def load_geometry(self): + with open(self.path) as f: + log = f.readlines() + for line in reversed(log): + if 'VIBRATIONAL ANALYSIS' in line: + end_ind = log.index(line) + start_ind = end_ind + for line in reversed(log[:end_ind]): + start_ind -= 1 + if 'Standard Nuclear Orientation' in line: + break + atom, coord = [], [] + geometry_flag = False + # Now look for the geometry. + # Will return the final geometry in the file under Standard Nuclear Orientation. + for line in log[start_ind+3:end_ind]: + if '------------' not in line: + line.strip() + data = line.split() + atom.append(data[1]) + coord.append([float(c) for c in data[2:]]) + else: + break + self.natom = len(atom) + self.nHcap = len(atom) - len(self.get_QM_ATOMS()) + return atom, coord + + def getXYZ(self): + """ + Return a string of the molecule in the XYZ file format. + """ + atoms, cart_coords = self.load_geometry() + natom = len(atoms) + xyz = str(natom) + '\n\n' + for i in range(natom): + xyz += '{:<2s} {:.10f} {:.10f} {:.10f}\n'.format(atoms[i],cart_coords[i][0],cart_coords[i][1],cart_coords[i][2]) + return xyz + + def write_files(self): + """ + Create files for GSM simulation in QMMM system. + """ + geometry = self.getXYZ() + g = open('QM_region.xyz', 'w') + g.write(geometry) + g.close() + + MM_xyz = self.get_fixed_molecule() + m = open('MM_region.txt', 'w') + m.write(MM_xyz) + m.close() + + link = self.get_USER_CONNECT() + f = open('link.txt', 'w') + f.writelines(link) + f.close() + + h = open('frozen.txt', 'w') + for i in range(self.natom-self.nHcap, self.natom): + h.write(str(i)) + h.write('\n') + h.close() + +parser = argparse.ArgumentParser(description='Qchem frequency files to SSM starting files') +parser.add_argument('file', metavar='FILE', type=str, nargs=1, + help='a Q-Chem freq job output file') +args = parser.parse_args() +input_file = args.file[0] +log = QChemLog(path=input_file) +log.write_files()