From 5322ca527cf9ec3b14bc5a037b1352fe9271ca69 Mon Sep 17 00:00:00 2001 From: arosale4 Date: Fri, 30 Nov 2018 12:03:56 -0500 Subject: [PATCH] Esp fit (#26) * ESP RMS for a charge scheme This is code that was borrowed from Anna Tomberg and modified to fit the Q2MM organization. Simply put I added two new commands "-mgESP" and "-mjESP" which take two files filename.mae,filename.chk and filename.mae,filename.in, respectively. I have not finished coding the "-mjESP" section, as I do not know of anyone who really wants the Jaguar end of this. For a loop, use the following syntax: RDAT -r some_ref_text_file.txt CDAT -mgESP filename.mae,filename.chk This will take filename.mae and determine the partial charges from the force field and then write them into a gaussian com file (filename.ESP.q2mm.com) which will calculate the RMS when compared to the ESP found in filename.chk. I added a few things into constants.py that should be expanded. 1) Gaussian methods with just a few common functions that we have used in the lab 2) CHELPG radii that are used. This is important for metals that do not have built in radii. A GaussCom class was added to filetypes, and the GaussLog class was modified to read in the RMS from the resulting filename.ESP.q2mm.log file. I have not yet checked to see how compatible this code is with more data types. I imagine there is probably some sorting errors that happen. * Bug: Early Return in GaussLog class After adding in the ESP_fit commit previously I included an early return argument in the read_out function of the GaussLog class. This had been remedied so it shouldn't be a problem anymore. * Will add comments later. * Minor bug fixes There were bugs where atom objects would not have consistent properties, such as an atom type corresponding to Ru but still contained old information from the merging process such as the atomic number. --- q2mm/calculate.py | 98 ++++++++++++++++++++++++++-- q2mm/compare.py | 10 ++- q2mm/constants.py | 35 +++++++++- q2mm/datatypes.py | 30 +++++++++ q2mm/filetypes.py | 158 +++++++++++++++++++++++++++++++++++++++------ q2mm/gradient.py | 4 +- q2mm/loop.py | 34 +++++----- screen/merge.py | 23 ++++--- tools/clean_mae.py | 1 + tools/setup_qm.py | 23 +++---- 10 files changed, 351 insertions(+), 65 deletions(-) diff --git a/q2mm/calculate.py b/q2mm/calculate.py index e1fc47e..f448d10 100755 --- a/q2mm/calculate.py +++ b/q2mm/calculate.py @@ -55,7 +55,7 @@ 'ma', 'mb', 'mt', 'me', 'meo', 'mea', 'meao', 'mh', 'mjeig', 'mgeig', - 'mp'] + 'mp', 'mgESP', 'mjESP'] # Commands related to Tinker. COM_TINKER = ['ta','tao', 'tb', 'tbo', 'tt','tto', 'te', 'teo', @@ -64,7 +64,7 @@ # Commands related to Amber. COM_AMBER = ['ae'] # All other commands. -COM_OTHER = ['r'] +COM_OTHER = ['r'] # All possible commands. COM_ALL = COM_GAUSSIAN + COM_JAGUAR + COM_MACROMODEL + COM_TINKER + \ COM_AMBER + COM_OTHER @@ -147,9 +147,23 @@ def main(args): os.path.join(opts.directory, filename)) inps[filename].commands = commands_for_filename inps[filename].write_com(sometext=opts.append) + #Has to be here even though this is a Gaussian Job. + if os.path.splitext(filename)[1] == '.chk': + # The generated com file will be used as the input filename. It + # also seems best to do the gaussian calculation in the + # collect_data function since we need to collect the force + # fields partial charges. + com_filename = os.path.splitext(filename)[0] + '.ESP.q2mm.com' + inps[com_filename] = filetypes.GaussCom( + os.path.join(opts.directory, com_filename)) + inps[com_filename].commands = commands_for_filename + inps[com_filename].read_newzmat(filename) + + + elif any(x in COM_TINKER for x in commands_for_filename): if os.path.splitext(filename)[1] == '.xyz': - inps[filename] = filetypes.Tinker_xyz( + inps[filename] = filetypes.TinkerXYZ( os.path.join(opts.directory, filename)) inps[filename].commands = commands_for_filename elif any(x in COM_AMBER for x in commands_for_filename): @@ -441,6 +455,17 @@ def return_calculate_parser(add_help=True, parents=None): help='Uses a MM3* FF file (somename.fld) and a parameter file ' '(somename.txt) to use the current FF parameter values as data. This ' 'is used for harmonic parameter tethering.') + mm_args.add_argument( + '-mgESP', type=str, nargs='+', action='append', + default=[], metavar='somename.mae,somename.chk', + help='Uses the partial charges obtained from the FF and *mae file to ' + 'determine the RMS of electrostatic fitting from a gaussain *chk file.') + mm_args.add_argument( + '-mjESP', type=str, nargs='+', action='append', + default=[], metavar='somename.mae,somename.in', + help='Uses the partial charges obtained from the FF and *mae file to ' + 'determine the RMS of electrostatic fitting from a schrodinger *in ' + 'file.') # TINKER OPTIONS tin_args = parser.add_argument_group("tinker data types") tin_args.add_argument( @@ -1228,6 +1253,68 @@ def collect_data(coms, inps, direc='.', sub_names=['OPT'], invert=None): src_1=filename, idx_1=idx_1 + 1, atm_1=atom.index)) + # MACROMODEL+GUASSIAN ESP + filenames = chain.from_iterable(coms['mgESP']) + for comma_filenames in filenames: + charges_list = [] + filename_mae, name_gau_chk = comma_filenames.split(',') + #Filename of the output *mae file (i.e. filename.q2mm.mae) + name_mae = inps[filename_mae].name_mae + mae = check_outs(name_mae, outs, filetypes.Mae, direc) + structures = filetypes.select_structures( + mae.structures, inps[filename_mae]._index_output_mae, 'pre') + for idx_1, structure in structures: + for atom in structure.atoms: + ### I think we want all the charges, right? + #if not 'b_q_use_charge' in atom.props or \ + # atom.props['b_q_use_charge']: + if atom.atomic_num > 0: + charges_list.append(atom.partial_charge) + com_filename = os.path.splitext(name_gau_chk)[0] + '.ESP.q2mm.com' + inps[com_filename].charge_list = charges_list + inps[com_filename].write_com() + inps[com_filename].run_gaussian() + name_gauss_log = inps[com_filename].name_log + gauss = check_outs(name_gauss_log, outs, filetypes.GaussLog, direc) + esp_rms = gauss.esp_rms + if esp_rms < 0.0: + raise Exception('A negative RMS was obtained for the ESP fitting ' + 'which indicates an error occured. Look at the ' + 'following file: {}'.format(name_gauss_log)) + data.append(datatypes.Datum( + val=esp_rms, + com='mgESP', + typ='esp', + src_1= name_mae, + src_2='gaussian', + idx_1 = 1)) + # MACROMODEL+JAGUAR ESP + ## This does not work, I still need to write code to support Jaguaer. -TR + filenames = chain.from_iterable(coms['mjESP']) + for comma_filenames in filenames: + charges_list = [] + name_mae, name_jag_chk = comma_filenames.split(',') + mae = check_outs(name_mae, outs, filetypes.Mae, direc) + structures = filetypes.select_structures( + mae.structures, inps[name_mae]._index_output_mae, 'pre') + for idx_1, structure in structures: + for atom in structure.atoms: + if not 'b_q_use_charge' in atom.props or \ + atom.props['b_q_use_charge']: + charges_list.append(atom.partial_charge) + ###Filler for ESP calculations#### + ### This is what is used in anna's code + current_RMS = run_ChelpG_inp.run_JCHelpG(charges_list,name_jag_chk) + + ### End of filler + if current_RMS < 0: + sys.exit("Error while computing RMS. Exiting") + data.append(datatypes.Datum( + val=current_RMS, + com='mjESP', + typ='esp', + src_1=name_mae, + idx_1=1)) # JAGUAR CHARGES EXCLUDING ALIPHATIC HYDROGENS filenames = chain.from_iterable(coms['jqh']) for filename in filenames: @@ -1732,7 +1819,7 @@ def lbl_to_data_attrs(datum, lbl): idxs = parts[-1] # if len(parts) == 4: if datum.typ in ['b','t','a']: - idxs = parts [-2] + idxs = parts[-2] atm_nums = parts[-1] atm_nums = atm_nums.split('-') for i, atm_num in enumerate(atm_nums): @@ -1740,6 +1827,9 @@ def lbl_to_data_attrs(datum, lbl): if datum.typ in ['p']: datum.src_1 = parts[1] idxs = parts[-1] + if datum.typ in ['esp']: + datum.src_1 = parts[1] + idxs = parts[-1] idxs = idxs.split('-') datum.idx_1 = int(idxs[0]) if len(idxs) == 2: diff --git a/q2mm/compare.py b/q2mm/compare.py index f3ccdd4..c919d24 100755 --- a/q2mm/compare.py +++ b/q2mm/compare.py @@ -130,6 +130,10 @@ def compare_data(r_dict, c_dict, output=None, doprint=False): for typ in r_dict: data_types.append(typ) data_types.sort() + total_num_energy = 0 + for typ in data_types: + if typ in ['e','eo','ea','eao']: + total_num_energy += len(r_dict[typ]) for typ in data_types: total_num += int(len(r_dict[typ])) if typ in ['e','eo','ea','eao']: @@ -142,7 +146,11 @@ def compare_data(r_dict, c_dict, output=None, doprint=False): diff = 360. - diff else: diff = r.val - c.val - score = (r.wht**2 * diff**2)/len(r_dict[typ]) + #score = (r.wht**2 * diff**2) + if typ in ['e', 'eo', 'ea', 'eao']: + score = (r.wht**2 * diff**2)/total_num_energy + else: + score = (r.wht**2 * diff**2)/len(r_dict[typ]) score_tot += score score_typ[c.typ] += score num_typ[c.typ] += 1 diff --git a/q2mm/constants.py b/q2mm/constants.py index 9dc6fca..522928e 100755 --- a/q2mm/constants.py +++ b/q2mm/constants.py @@ -83,7 +83,9 @@ 'imp2': 0.2, 'sb': 0.2, 'q': 0.1, - 'q_p': 0.05 + 'q_p': 0.05, + 'vdwr': 0.1, + 'vdwfc': 0.02 } # WEIGHTS @@ -102,6 +104,7 @@ 'q': 10.00, 'qh': 10.00, 'qa': 10.00, + 'esp': 10.00, 'p': 10.00 } @@ -281,3 +284,33 @@ ('Lr', 262.109692000) ] ) + + +# ELECTRONIC STRUCTURE METHODS +gaussian_methods = [ 'b3lyp', + 'm06', + 'm062x', + 'm06L'] + + + + +# CHELPG NEEDED RADII +# These are a combination of the default values in Gaussian09 (first and second +# row?) and values we have used in the past for metals. +CHELPG_RADII = OrderedDict( + [ + ('H', 1.45), + ('C', 1.50), + ('N', 1.70), + ('O', 1.70), + ('F', 1.70), + ('Pd', 2.40), + ('Ir', 2.40), + ('Ru', 2.40), + ('Rh', 2.40), + ('S', 2.00) + ] + ) + + diff --git a/q2mm/datatypes.py b/q2mm/datatypes.py index dd24bd5..1ab01b0 100755 --- a/q2mm/datatypes.py +++ b/q2mm/datatypes.py @@ -848,6 +848,33 @@ def import_ff(self, path=None, sub_search='OPT'): mm3_label = line[:2], value = parm_cols[1]))) continue + # Bonds. + elif match_mm3_vdw(line): + logger.log( + 5, '[L{}] Found vdw:\n{}'.format( + i + 1, line.strip('\n'))) + if section_sub: + atm_lbls = [line[4:6], line[8:10]] + atm_typs = self.convert_to_types( + atm_lbls, self.atom_types[-1]) + parm_cols = line[P_1_START:P_3_END] + parm_cols = map(float, parm_cols.split()) + self.params.extend(( + ParamMM3(atom_labels = atm_lbls, + atom_types = atm_typs, + ptype = 'vdwr', + mm3_col = 1, + mm3_row = i + 1, + mm3_label = line[:2], + value = parm_cols[0]), + ParamMM3(atom_labels = atm_lbls, + atom_types = atm_typs, + ptype = 'vdwfc', + mm3_col = 2, + mm3_row = i + 1, + mm3_label = line[:2], + value = parm_cols[1]))) + continue # The Van der Waals are stored in annoying way. if line.startswith('-6'): section_vdw = True @@ -1215,6 +1242,9 @@ def match_mm3_label(mm3_label): in a Schrodinger mm3.fld file. """ return re.match('[\s5a-z][1-5]', mm3_label) +def match_mm3_vdw(mm3_label): + """Matches MM3* label for bonds.""" + return re.match('[\sa-z]6', mm3_label) def match_mm3_bond(mm3_label): """Matches MM3* label for bonds.""" return re.match('[\sa-z]1', mm3_label) diff --git a/q2mm/filetypes.py b/q2mm/filetypes.py index 0cf5d43..c196c70 100755 --- a/q2mm/filetypes.py +++ b/q2mm/filetypes.py @@ -466,6 +466,107 @@ def run(self,check_tokens=False): \n======================\n") os.chdir(current_directory) +class GaussCom(File): + """ + Used to write Gaussian command files to run quick calculations. The ony use + thus far is for checking the RMS of the electrostatic potential using + partial charges obtained from a force field. + """ + def __init__(self, path): + super(GaussCom, self).__init__(path) + self.name = os.path.splitext(self.filename)[0] + self.commands = None + self.old_chk = None + self.name_chk = self.name + '.chk' + self.name_log = self.name + '.log' + self.atom_and_coords = [] + self.charge = None + self.charge_list = None + self.multiplicity = 1 + self.memory = 24 + self.procs = 24 + self.method = 'M06' + + def read_newzmat(self, old_check): + self.old_chk = old_check + if os.path.isfile('TEMP.com'): + os.remove('TEMP.com') + sp.call('newzmat -ichk -ocart {} {}'.format(self.old_chk, 'TEMP'), + shell=True) + with open('TEMP.com', 'r') as f: + for line in f: + if re.search('[A-Z][a-z]?\s+{0}\s+{0}\s+{0}'.format( + co.RE_FLOAT),line): + ele, x, y, z = line.split() + self.atom_and_coords.append([ele,float(x),float(y),float(z)]) + if re.search('[+-]?\d,\d',line): + charge, multiplicity = line.split(',') + self.charge = charge + self.multiplicity = multiplicity + if re.search('^[#]',line): + split = line.split() + for keyword in split: + for method in co.gaussian_methods: + if method.upper() in keyword.upper(): + self.method = method.upper() + #os.remove('TEMP.com') + + def write_com(self): + """ + Gaussian has a utility, newzmat, which could write the com file, but + in my experience it is rather slow especially for larger chkpoint files. + It is easier just writting the file ourselves. -Tony + """ + elements_in_structure = [] + with open(self.filename, 'w') as f: + if self.old_chk: + f.write('%oldchk={}\n'.format(self.old_chk)) + f.write('%chk={}\n'.format(self.name_chk)) + f.write('%mem={}GB\n'.format(self.memory)) + f.write('%nprocshared={}\n'.format(self.procs)) + f.write('#N Guess=TCheck SCRF=Check GenChk pop=(chelpg,readradii) ' + 'IOp(6/20=30133) chkbasis {}\n\n'.format(self.method)) + f.write('SP with fixed charges\n\n') + f.write('{} {}'.format(self.charge,self.multiplicity)) + if self.charge_list: + try: + if len(self.charge_list) != len(self.atom_and_coords): + raise ValueError("Length of the charge list and atom " + "list are not the same: {}".format( + self.filename)) + except: + print("Just print something") + for atom, charge in itertools.izip(self.atom_and_coords,self.charge_list): + if atom[0] not in elements_in_structure: + elements_in_structure.append(atom[0]) + f.write(' '.join((' {}--{:8.6f}'.format(atom[0],charge), + '{:14.10f}'.format(atom[1]), + '{:14.10f}'.format(atom[2]), + '{:14.10f}'.format(atom[3]), + '\n' + ))) + f.write('\n') + for atom in elements_in_structure: + f.write('{0} {1}\n'.format(atom,co.CHELPG_RADII[atom])) + f.write('\n') + + def run_gaussian(self): + """ + This runs the gaussian job. I have this labeled differently than our + typical "run" functions because I don't want to use this function until + after we have calculated and collected partial charge data. + """ + logger.log(5, 'RUNNING: {}'.format(self.filename)) + self._index_output_log = [] + current_directory = os.getcwd() + os.chdir(self.directory) + if os.path.isfile(self.name_log): + os.remove(self.name_log) + if os.path.isfile(self.name_chk): + os.remove(self.name_chk) + sp.call('g09 {}'.format(self.filename), shell=True) + os.chdir(current_directory) + class GaussFormChk(File): """ Used to retrieve data from Gaussian formatted checkpoint files. @@ -531,6 +632,7 @@ def __init__(self, path): self._evals = None self._evecs = None self._structures = None + self._esp_rms = None @property def evecs(self): if self._evecs is None: @@ -547,6 +649,12 @@ def structures(self): # self.read_out() self.read_archive() return self._structures + @property + def esp_rms(self): + if self._esp_rms is None: + self._esp_rms = -1 + self.read_out() + return self._esp_rms def read_out(self): """ Read force constant and eigenvector data from a frequency @@ -575,8 +683,12 @@ def read_out(self): except: # End of file. break + if 'Charges from ESP fit' in line: + pattern = re.compile('RMS=\s+({0})'.format(co.RE_FLOAT)) + match = pattern.search(line) + self._esp_rms = float(match.group(1)) # Gathering some geometric information. - if 'orientation:' in line: + elif 'orientation:' in line: self._structures.append(Structure()) file_iterator.next() file_iterator.next() @@ -736,24 +848,26 @@ def read_out(self): # We know we're done if this is in the line. if 'Harmonic' in line: break - for evec in self._evecs: - # Each evec is a single eigenvector. - # Add up the sum of squares over an eigenvector. - sum_of_squares = 0. - # Appropriately named, element is an element of that single - # eigenvector. - for element in evec: - sum_of_squares += element * element - # Now x is the inverse of the square root of the sum of squares - # for an individual eigenvector. - element = 1 / np.sqrt(sum_of_squares) - for i in range(len(evec)): - evec[i] *= element - self._evals = np.array(self._evals) - self._evecs = np.array(self._evecs) - logger.log(1, '>>> self._evals: {}'.format(self._evals)) - logger.log(1, '>>> self._evecs: {}'.format(self._evecs)) - logger.log(5, ' -- {} structures found.'.format(len(self.structures))) + if self._evals and self._evecs: + for evec in self._evecs: + # Each evec is a single eigenvector. + # Add up the sum of squares over an eigenvector. + sum_of_squares = 0. + # Appropriately named, element is an element of that single + # eigenvector. + for element in evec: + sum_of_squares += element * element + # Now x is the inverse of the square root of the sum of squares + # for an individual eigenvector. + element = 1 / np.sqrt(sum_of_squares) + for i in range(len(evec)): + evec[i] *= element + self._evals = np.array(self._evals) + self._evecs = np.array(self._evecs) + logger.log(1, '>>> self._evals: {}'.format(self._evals)) + logger.log(1, '>>> self._evecs: {}'.format(self._evecs)) + logger.log(5, ' -- {} structures found.'.format( + len(self.structures))) # May want to move some attributes assigned to the structure class onto # this filetype class. def read_archive(self): @@ -1634,7 +1748,8 @@ def get_com_opts(self): com_opts['strs'] = True if any(x in ['jb', 'ja', 'jt'] for x in self.commands): com_opts['sp_mmo'] = True - if any(x in ['me', 'mea', 'mq', 'mqh', 'mqa'] for x in self.commands): + if any(x in ['me', 'mea', 'mq', 'mqh', 'mqa', 'mgESP', 'mjESP'] + for x in self.commands): com_opts['sp'] = True # Command meig is depreciated. if any(x in ['mh', 'meig', 'mjeig', 'mgeig'] for x in self.commands): @@ -1704,6 +1819,7 @@ def write_com(self, sometext=None): else: com += co.COM_FORM.format('MMOD', 0, 1, 0, 0, 0, 0, 0, 0) # May want to turn on/off arg2 (continuum solvent). + #com += co.COM_FORM.format('FFLD', 2, 0, 0, 0, 36.7, 0, 0, 0) com += co.COM_FORM.format('FFLD', 2, 0, 0, 0, 0, 0, 0, 0) # Also may want to turn on/off cutoffs using BDCO. ## We have noticed there are some oddities for electrostatic @@ -1713,7 +1829,7 @@ def write_com(self, sometext=None): ## EXNB is used to set all vdW and electrostatic cutoffs to 99.0 ## ensuring all interactions are gathered. The seventh column is the ## cutoff for hydrogen bonds, and 4 is the default value. - #com += co.COM_FORM.format('EXNB', 0, 0, 0, 0, 99., 99., 4., 0) + com += co.COM_FORM.format('EXNB', 0, 0, 0, 0, 99., 99., 4., 0) if com_opts['strs']: com += co.COM_FORM.format('BGIN', 0, 0, 0, 0, 0, 0, 0, 0) # Look into differences. diff --git a/q2mm/gradient.py b/q2mm/gradient.py index 3075b49..34f32b8 100755 --- a/q2mm/gradient.py +++ b/q2mm/gradient.py @@ -82,9 +82,9 @@ def __init__(self, # Whether or not to generate parameters with these methods. self.do_lstsq = False - self.do_lagrange = True + self.do_lagrange = False self.do_levenberg = False - self.do_newton = True + self.do_newton = False self.do_svd = False # Particular settings for each method. diff --git a/q2mm/loop.py b/q2mm/loop.py index 2677649..066f572 100755 --- a/q2mm/loop.py +++ b/q2mm/loop.py @@ -186,7 +186,7 @@ def run_loop_input(self, lines, score=None): if "lstsq" in col: g_args = col.split('=')[1].split(',') for arg in g_args: - if arg == True: + if arg == "True": grad.do_lstsq=True elif arg == False: grad.do_lstsq=False @@ -215,7 +215,7 @@ def run_loop_input(self, lines, score=None): elif "newton" in col: g_args = col.split('=')[1].split(',') for arg in g_args: - if arg == True: + if arg == "True": grad.do_newton=True elif arg == False: grad.do_newton=False @@ -244,7 +244,7 @@ def run_loop_input(self, lines, score=None): elif "levenberg" in col: g_args = col.split('=')[1].split(',') for arg in g_args: - if arg == True: + if arg == "True": grad.do_levenberg=True elif arg == False: grad.do_levenberg=False @@ -279,48 +279,48 @@ def run_loop_input(self, lines, score=None): else: for val in factor_vals: grad.levenberg_factor.append(float(val)) - elif "legrange" in col: + elif "lagrange" in col: g_args = col.split('=')[1].split(',') for arg in g_args: - if arg == True: - grad.do_legrange=True + if arg == "True": + grad.do_lagrange=True elif arg == False: - grad.do_legrange=False + grad.do_lagrange=False if 'radii' in arg: - grad.legrange_radii = [] + grad.lagrange_radii = [] radii_vals = re.search( r"\[(.+)\]",arg).group(1).split('/') if radii_vals=='None': - grad.legrange_radii = None + grad.lagrange_radii = None else: for val in radii_vals: - grad.legrange_radii.append(float(val)) + grad.lagrange_radii.append(float(val)) if 'cutoff' in arg: - grad.legrange_cutoff = [] + grad.lagrange_cutoff = [] cutoff_vals = re.search( r"\[(.+)\]",arg).group(1).split('/') if cutoff_vals=='None': - grad.legrange_cutoff = None + grad.lagrange_cutoff = None else: if len(cutoff_vals) > 2 or \ len(cutoff_vals) < 2: raise Exception("Cutoff values must " \ "be between two numbers.") for val in cutoff_vals: - grad.legrange_cutoff.append(float(val)) + grad.lagrange_cutoff.append(float(val)) if 'factor' in arg: - grad.legrange_factors = [] + grad.lagrange_factors = [] factor_vals = re.search( r"\[(.+)\]",arg).group(1).split('/') if factor_vals=='None': - grad.legrange_factors = None + grad.lagrange_factors = None else: for val in factor_vals: - grad.legrange_factors.append(float(val)) + grad.lagrange_factors.append(float(val)) elif "svd" in col: g_args = col.split('=')[1].split(',') for arg in g_args: - if arg == True: + if arg == "True": grad.do_svd=True elif arg == False: grad.do_svd=False diff --git a/screen/merge.py b/screen/merge.py index 381fd52..1925db2 100755 --- a/screen/merge.py +++ b/screen/merge.py @@ -554,8 +554,15 @@ def merge_structures_from_matching_atoms(struct_1, match_1, struct_2, match_2): if common_atom_1.atom_type == 64: common_atom_1.atom_type = common_atom_2.atom_type common_atom_1.color = common_atom_2.color - common_atom_1.atom_type = common_atom_2.atom_type + #common_atom_1.atom_type = common_atom_2.atom_type common_atom_1.color = common_atom_2.color + #For somereason the rest of the properties of the atom class aren't + #updated when the atom_type is changed. Here we set the atom type to + #what the atom type is, which is redundant, but it resets the atomic + #number and weight. The AtomName isn't important, but it makes it look + #pretty in the final mae file. + common_atom_1._setAtomType(common_atom_1.atom_type) + common_atom_1._setAtomName(str(common_atom_1.element) + str(common_atom_1.index)) # Below are alternatives options for which atoms to keep. Currently, the # coordinates of the atoms from struct_1 are kept. @@ -666,7 +673,7 @@ def merge_structures_from_matching_atoms(struct_1, match_1, struct_2, match_2): reordered_atoms.extend(interest) merge = build.reorder_atoms(merge,reordered_atoms) - print(reordered_atoms) +# print(reordered_atoms) for bond in merge.bond: for prop in bond.property: initial_index = bond.property[prop] @@ -691,11 +698,11 @@ def merge_structures_from_matching_atoms(struct_1, match_1, struct_2, match_2): frozen_atoms=frozen_atoms, fix_torsions=fix_torsions)[0] # Short conformational sampling. - merge = mcmm( - [merge], - #frozen_atoms=range(1, num_atoms + 1))[0] - frozen_atoms=frozen_atoms)[0] - # Do another minimization, this time without frozen atoms. +# merge = mcmm( +# [merge], +# #frozen_atoms=range(1, num_atoms + 1))[0] +# frozen_atoms=frozen_atoms)[0] +# # Do another minimization, this time without frozen atoms. #if fix_torsions: # merge = mini( # [merge], @@ -703,7 +710,7 @@ def merge_structures_from_matching_atoms(struct_1, match_1, struct_2, match_2): merge = mini( [merge], fix_torsions=fix_torsions)[0] - merge = mcmm([merge])[0] +# merge = mcmm([merge])[0] return merge def add_chirality(structure): diff --git a/tools/clean_mae.py b/tools/clean_mae.py index 5e6f557..28cdd40 100755 --- a/tools/clean_mae.py +++ b/tools/clean_mae.py @@ -48,6 +48,7 @@ 'i_m_formal_charge', 's_m_color_rgb', 's_m_atom_name', + 's_m_atomic_number', 's_m_label_format', 'i_m_label_color', 's_m_label_user_text', diff --git a/tools/setup_qm.py b/tools/setup_qm.py index ffdd91b..712a738 100644 --- a/tools/setup_qm.py +++ b/tools/setup_qm.py @@ -142,7 +142,7 @@ def determine_method(self): if self.frozen_atoms: self.opt += ' geom=modredundant' if self.calculation_type == 'FZTS': - self.opt += ' opt=(calcfc,maxcycle=50)' + self.opt += ' opt=(calcfc,maxcycle=500)' if self.calculation_type == 'TS': # Do I need this part? Will self.frequency ever be a value other # than None or freq=noraman? @@ -203,14 +203,14 @@ def get_frozen_coords(self, dict_of_patterns, from_command_line=False): dict_of_patterns = self.get_dictionary_of_frozen_coords( dict_of_patterns) for pattern in dict_of_patterns: - atom_indicies = analyze.evaluate_substructure(self.mae_struct, + matches = analyze.evaluate_substructure(self.mae_struct, pattern, - first_match_only=False)[0] - str_indicies = [] - for index in atom_indicies: - str_indicies.append(str(index)) - - frozen_coord_lines.append('{} {} F'.format( + first_match_only=False) + for atom_indicies in matches: + str_indicies = [] + for index in atom_indicies: + str_indicies.append(str(index)) + frozen_coord_lines.append('{} {} F'.format( dict_of_patterns[pattern], ' '.join(str_indicies))) return frozen_coord_lines @@ -263,13 +263,14 @@ def get_atoms_from_schrodinger_struct(structures_list): structures = {} for struct in structures_list: atom_list = [] - dummy = None + dummy = [] for atom in struct.atom: + atom._setAtomType(atom.atom_type) if atom.atomic_number < 0: - dummy = atom + dummy.append(atom.index) atom_list.append((atom.element, atom.x, atom.y, atom.z)) if dummy: - struct.deleteAtoms((dummy.index),renumber_map=False) + struct.deleteAtoms(dummy,renumber_map=False) if 'r_mmod_Potential_Energy-MM3*' in struct.property: energy = struct.property['r_mmod_Potential_Energy-MM3*'] structures[struct] = [struct.formal_charge, atom_list, energy]