Update

DP4.py

@@ -214,15 +214,22 @@ def DP4j(Clabels, Cvalues, Hlabels, Hvalues, Cexp, Hexp, cJvals, cJlabels,
 
 
 def CP3(Clabels, Cvalues1, Cvalues2, Hlabels, Hvalues1, Hvalues2,
-        Cexp1, Cexp2, Hexp1, Hexp2, settings):
+        Cexp1, Cexp2, Hexp1, Hexp2,
+        CalcFile1, CalcFile2, NMRFile1, NMRFile2, settings):
 
     sortedClabels1, sortedCvalues1, sortedCexp1, _ =\
         AssignExpNMR(Clabels, Cvalues1, Cexp1)
+    print "Clabels: " + str(Clabels)
+    print "Cvalues: " + str(Cvalues1)
+    print "Cexp1: " + str(Cexp1)
     Cerrors1 = [a-b for a,b in zip(sortedCvalues1,sortedCexp1)]
     sortedClabels2, sortedCvalues2, sortedCexp2, _ =\
         AssignExpNMR(Clabels, Cvalues2, Cexp2)
     Cerrors2 = [a-b for a,b in zip(sortedCvalues2,sortedCexp2)]
 
+    print "Hlabels: " + str(Hlabels)
+    print "Hvalues: " + str(Hvalues1)
+    print "Hexp1: " + str(Hexp1)
     sortedHlabels1, sortedHvalues1, sortedHexp1, _ =\
         AssignExpNMR(Hlabels, Hvalues1, Hexp1)
     Herrors1 = [a-b for a,b in zip(sortedHvalues1,sortedHexp1)]
@@ -267,21 +274,28 @@ def CP3(Clabels, Cvalues1, Cvalues2, Hlabels, Hvalues1, Hvalues2,
     CP3Cprob2 = stats.norm(CP3meanC, CP3stdevC).pdf(CP3c2)
     CP3Hprob2 = stats.norm(CP3meanH, CP3stdevH).pdf(CP3h2)
     CP3allprob2 = stats.norm(CP3meanAll, CP3stdevAll).pdf(CP3all2)
-
-    Print("Probabilities for Exp1-Calc1 & Exp2-Calc2 assignment")
-    Print("Based on C data:   " + format(CP3Cprob1, "4.1f") + "%")
-    Print("Based on H data:   " + format(CP3Hprob1, "4.1f") + "%")
-    Print("Based on all data: " + format(CP3allprob1, "4.1f") + "%")
-    Print("Probabilities for Exp1-Calc2 & Exp2-Calc1 assignment")
-    Print("Based on C data:   " + format(CP3Cprob2, "4.1f") + "%")
-    Print("Based on H data:   " + format(CP3Hprob2, "4.1f") + "%")
-    Print("Based on all data: " + format(CP3allprob2, "4.1f") + "%")
+    print NMRFile1
+    print NMRFile2
+    print CalcFile1
+    print CalcFile2
+    comb1 = NMRFile1 + '-' + CalcFile1 + ' & ' + NMRFile2 + '-' + CalcFile2
+    comb2 = NMRFile1 + '-' + CalcFile2 + ' & ' + NMRFile2 + '-' + CalcFile1
+    Print("Probabilities for " + comb1 + " assignment")
+    Print("Based on C data:   " + format(CP3Cprob1*100, "4.1f") + "%")
+    Print("Based on H data:   " + format(CP3Hprob1*100, "4.1f") + "%")
+    Print("Based on all data: " + format(CP3allprob1*100, "4.1f") + "%")
+    Print("Probabilities for " + comb2 + " assignment")
+    Print("Based on C data:   " + format(CP3Cprob2*100, "4.1f") + "%")
+    Print("Based on H data:   " + format(CP3Hprob2*100, "4.1f") + "%")
+    Print("Based on all data: " + format(CP3allprob2*100, "4.1f") + "%")
     if CP3allprob1>CP3allprob2:
-        Print("Based on all data, Exp1-Calc1 & Exp2-Calc2 assignment is\n"+
+        Print("Based on all data, " + comb1 + " assignment is\n"+
               "likely the correct one")
     else:
-        Print("Based on all data, Exp1-Calc2 & Exp2-Calc1 assignment is\n"+
+        Print("Based on all data, " + comb2 + " assignment is\n"+
               "likely the correct one")
+
+    return output
 
 
 def f3(deltaExp, deltaCalc):
@@ -290,6 +304,7 @@ def f3(deltaExp, deltaCalc):
     else:
         return deltaExp*deltaCalc
 
+
 def CalcCP3(values1, values2, exp1, exp2):
     deltaCalcs = [a - b for a,b in zip(values1, values2)]
     deltaExps = [a - b for a,b in zip(exp1, exp2)]
@@ -1025,6 +1040,8 @@ def CalculateRKDE(errors, shifts, ParamFile, t, nucleus):
     for i, e in enumerate(errors):
         region = bisect.bisect_left(regions, shifts[i])
         ep5 = ep5*float((kdes[region])(e)[0])
+        if ep5<1e-315:  #Crude underflow protection
+            ep5=1e-315
     return ep5
 
 

NMRAnalysis.py

@@ -44,6 +44,8 @@ def main(numDS, settings, *args):
         Cexp, Hexp, equivs, omits = ReadExpNMR(args[-1])
     else:
         [NMRfile1, NMRfile2] = args[-1].split(',')
+        CalcFile1 = args[1]
+        CalcFile2 = args[3]
         Cexp1, Hexp1, equivs, omits = ReadExpNMR(NMRfile1)
         Cexp2, Hexp2, equivs, omits = ReadExpNMR(NMRfile2)
 
@@ -133,7 +135,8 @@ def main(numDS, settings, *args):
         elif settings.CP3:
             CP3outp = DP4.CP3(Clabels, OptCvalues[0],OptCvalues[1],
                               Hlabels, OptHvalues[0], OptHvalues[1],
-                              Cexp1, Cexp2, Hexp1, Hexp2, settings)
+                              Cexp1, Cexp2, Hexp1, Hexp2, 
+                              CalcFile1, CalcFile2, NMRfile1, NMRfile2, settings)
             return '\n'.join(CP3outp) + '\n'
         else:
             DP4outp = DP4.DP4(Clabels, OptCvalues, Hlabels, OptHvalues, Cexp,

NWChem.py

@@ -97,6 +97,19 @@ def WriteNWChemFile(NWCheminp, conformer, atoms, charge, settings):
         basis = '6-311g*'
 
     f.write('end\n\nbasis\n  * library ' + basis + '\nend\n\n')
+    if settings.Solvent != "":
+        f.write('cosmo\n  do_cosmo_smd true\n  solvent ' + settings.Solvent + '\n')
+        f.write('end\n\n')
+    if settings.DFTOpt or settings.HFOpt:
+        f.write('driver\n  maxiter ' + str(settings.MaxDFTOptCycles)+ '\nend\n\n')
+        if settings.DFTOpt:
+            f.write('dft\n  xc b3lyp\n  mult 1\nend\n\n')
+            f.write('task dft optimize\n\n')
+        if settings.M06Opt:
+            f.write('dft\n  xc m06-2x\n  mult 1\nend\n\n')
+            f.write('task dft optimize\n\n')
+        if settings.HFOpt:
+            f.write('task scf optimize\n\n')
     if (settings.Functional).lower() == 'b3lyp':
         f.write('dft\n  xc b3lyp\n  mult 1\nend\n\n')
     elif (settings.Functional).lower() == 'm062x' or\
@@ -106,9 +119,6 @@ def WriteNWChemFile(NWCheminp, conformer, atoms, charge, settings):
         f.write('dft\n  xc mpw91 0.75 HFexch 0.25 perdew91\n  mult 1\nend\n\n')
     else:
         f.write('dft\n  xc ' + settings.Functional + '\n  mult 1\nend\n\n')
-    if settings.Solvent != "":
-        f.write('cosmo\n  do_cosmo_smd true\n  solvent ' + settings.Solvent + '\n')
-        f.write('end\n\n')
     f.write('task dft energy\n\n')
     f.write('property\n  shielding\nend\n')
     f.write('task dft property\n')

PyDP4.py

@@ -52,8 +52,8 @@
 
 #Assigning the config default values
 class Settings:
-    MMTinker = True
-    MMMacromodel = False
+    MMTinker = False
+    MMMacromodel = True
     DFT = 'z'
     Rot5Cycle = False
     Title = 'DP4molecule'
@@ -69,7 +69,7 @@ class Settings:
     PM7Opt = False
     HFOpt = False
     M06Opt = False
-    MaxOptReruns = 2
+    MaxDFTOptCycles = 50
     jKarplus = False
     jFC = False
     jJ = False
@@ -87,6 +87,7 @@ class Settings:
     RenumberFile = ''
     ScriptDir = ''
     user = 'ke291'
+    OnlyConfS = False # Stop the process after conformational search
     MMstepcount = 10000
     MMfactor = 2500  # nsteps = MMfactor*degrees of freedom
     HardConfLimit = 10000
@@ -238,6 +239,10 @@ def main(filename, ExpNMR, nfiles):
             print '\nRunning Macromodel...'
             MacroModel.RunMacromodel(len(inpfiles), settings, *mminpfiles)
 
+    if settings.OnlyConfS:
+        print "Conformational search completed, quitting as instructed."
+        quit()
+
     if (not settings.AssumeDone) and (not settings.UseExistingInputs):
         if settings.ConfPrune:
             if settings.DFT == 'z' or settings.DFT == 'g' or settings.DFT == 'd':
@@ -274,7 +279,9 @@ def main(filename, ExpNMR, nfiles):
     elif settings.DFT == 'n' or settings.DFT == 'w' or settings.DFT == 'm':
         Files2Run = NWChem.GetFiles2Run(inpfiles, settings)
     print Files2Run
+
     if len(Files2Run) == 0:
+        """
         if (settings.DFT == 'z' or settings.DFT == 'g' or settings.DFT == 'd') and\
             (settings.DFTOpt or settings.PM6Opt or settings.HFOpt or settings.M06Opt):
             print "Checking if all geometries have converged"
@@ -291,7 +298,9 @@ def main(filename, ExpNMR, nfiles):
                 QRun = True
         else:
             QRun = True
-
+        """
+        QRun = True
+
     if len(Files2Run) > settings.HardConfLimit:
         print "Hard conformation count limit exceeded, DFT calculations aborted."
         quit()
@@ -461,8 +470,8 @@ def SetTMSConstants():
     parser = argparse.ArgumentParser(description='PyDP4 script to setup\
     and run Tinker, Gaussian (on ziggy) and DP4')
     parser.add_argument('-m', '--mm', help="Select molecular mechanics program,\
-    t for tinker or m for macromodel, default is t", choices=['t', 'm'],
-    default='t')
+    t for tinker or m for macromodel, default is m", choices=['t', 'm'],
+    default='m')
     parser.add_argument('-d', '--dft', help="Select DFT program, j for Jaguar,\
     g for Gaussian, n for NWChem, z for Gaussian on ziggy, d for Gaussian on \
     Darwin, w for NWChem on ziggy, m for NWChem on Medivir cluster, default is z",
@@ -523,6 +532,8 @@ def SetTMSConstants():
     DFT inputs, avoids long conf pruning and regeneration", action="store_true")
     parser.add_argument("--NoConfPrune", help="Skip RMSD pruning, use all\
     conformers in the energy window", action="store_true")
+    parser.add_argument("--OnlyConfS", help="Quit after conformational search",
+                        action="store_true")
     parser.add_argument("--Renumber", help="Renumber the atoms in\
     diastereomers according to renumbering map in the specified file. Useful\
     when analysing manually drawn input structures")
@@ -612,6 +623,8 @@ def SetTMSConstants():
     else:
         settings.MMMacromodel = True
         settings.MMTinker = False
+    if args.OnlyConfS:
+        settings.OnlyConfS = True
     if args.DFTOpt:
         settings.DFTOpt = True
     if args.M06Opt:
@@ -665,11 +678,14 @@ def SetTMSConstants():
             print "Stats file not found, quitting."
             quit()
 
+    #settings.SCHRODINGER = os.environ['SCHRODINGER']
+    settings.SCHRODINGER = '/usr/local/shared/schrodinger/current'
+
     inpfiles = [x.split('.')[0] for x in args.StructureFiles]
 
     if len(inpfiles) == 1:
         main(inpfiles[0], args.ExpNMR, 1)
     else:
         main(inpfiles, args.ExpNMR, len(inpfiles))
-    
+
     #main()

README

@@ -3,9 +3,9 @@
 PyDP4 workflow integrating MacroModel/TINKER, Gaussian/NWChem
 and DP4 analysis
 
-version 0.4
+version 0.7
 
-Copyright (c) 2015 Kristaps Ermanis, Jonathan M. Goodman
+Copyright (c) 2016 Kristaps Ermanis, Jonathan M. Goodman
 
 distributed under MIT license
 
@@ -22,18 +22,18 @@ CONTENTS
 
 REQUIREMENTS AND SETUP
 
-All the python and Java files and one utility to convert to and from TINKER
+All the python  files and one utility to convert to and from TINKER
 nonstandard xyz file are in the attached archive. They are set up to work
 from a centralised location. It is not a requirement, but it is probably
 best to add the location to the PATH variable.
 
 The script currently is set up to use MacroModel for molecular mechanics and
-NWChem for DFT and it runs NWChem locally. Gaussian and TINKER is also supported
+Gaussian for DFT and it runs Gaussian on ziggy by default. NWChem and TINKER is also supported.
 This setup has several requirements.
 
 1) One should have MacroModel or TINKER and NWChem or Gaussian.
 The beginning PyDP4.py file contains a structure "Settings", where the location
-of the TINKER scan executable or MacroModel bmin executable should be specified.
+of the TINKER scan executable or MacroModel bmin executable should be specified. All instances of 'ke291' should be replaced with the CRSid of the user.
 
 2) The various manipulations of sdf files (renumbering, ring corner flipping)
 requires OpenBabel, including Python bindings. The following links provide