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

@@ -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,14 +1819,17 @@ 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):
             setattr(datum, 'atm_{}'.format(i+1), int(atm_num))
     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:

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

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)
+        ]
+    )
+
+

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)