From 4f83e46c153b823d56edcf5690c35c433d6d6268 Mon Sep 17 00:00:00 2001 From: Kristaps Ermanis Date: Sun, 6 Nov 2016 19:05:20 +0000 Subject: [PATCH] Update --- ConfPrune.pyx | 0 DP4.py | 41 +++++++++++++++++++++++++++++------------ Defaultslurm | 0 FiveConf.py | 0 Gaussian.py | 0 InchiGen.py | 0 Karplus.py | 0 LICENSE | 0 MacroModel.py | 0 NMRAnalysis.py | 5 ++++- NWChem.py | 16 +++++++++++++--- PyDP4.py | 30 +++++++++++++++++++++++------- README | 10 +++++----- TMSdata | 0 default.com | 0 default.key | 0 default1.com | 0 nmrPredictGaus.py | 0 nmrPredictNWChem.py | 0 pcs-1 | 0 pcs-2 | 0 21 files changed, 74 insertions(+), 28 deletions(-) mode change 100644 => 100755 ConfPrune.pyx mode change 100644 => 100755 DP4.py mode change 100644 => 100755 Defaultslurm mode change 100644 => 100755 FiveConf.py mode change 100644 => 100755 Gaussian.py mode change 100644 => 100755 InchiGen.py mode change 100644 => 100755 Karplus.py mode change 100644 => 100755 LICENSE mode change 100644 => 100755 MacroModel.py mode change 100644 => 100755 NMRAnalysis.py mode change 100644 => 100755 NWChem.py mode change 100644 => 100755 README mode change 100644 => 100755 TMSdata mode change 100644 => 100755 default.com mode change 100644 => 100755 default.key mode change 100644 => 100755 default1.com mode change 100644 => 100755 nmrPredictGaus.py mode change 100644 => 100755 nmrPredictNWChem.py mode change 100644 => 100755 pcs-1 mode change 100644 => 100755 pcs-2 diff --git a/ConfPrune.pyx b/ConfPrune.pyx old mode 100644 new mode 100755 diff --git a/DP4.py b/DP4.py old mode 100644 new mode 100755 index 6172ff9..9222301 --- a/DP4.py +++ b/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 diff --git a/Defaultslurm b/Defaultslurm old mode 100644 new mode 100755 diff --git a/FiveConf.py b/FiveConf.py old mode 100644 new mode 100755 diff --git a/Gaussian.py b/Gaussian.py old mode 100644 new mode 100755 diff --git a/InchiGen.py b/InchiGen.py old mode 100644 new mode 100755 diff --git a/Karplus.py b/Karplus.py old mode 100644 new mode 100755 diff --git a/LICENSE b/LICENSE old mode 100644 new mode 100755 diff --git a/MacroModel.py b/MacroModel.py old mode 100644 new mode 100755 diff --git a/NMRAnalysis.py b/NMRAnalysis.py old mode 100644 new mode 100755 index e37495c..1bb5c47 --- a/NMRAnalysis.py +++ b/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, diff --git a/NWChem.py b/NWChem.py old mode 100644 new mode 100755 index 944265d..861ce06 --- a/NWChem.py +++ b/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') diff --git a/PyDP4.py b/PyDP4.py index 944eb25..3237d60 100755 --- a/PyDP4.py +++ b/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() diff --git a/README b/README old mode 100644 new mode 100755 index 515d494..434100d --- a/README +++ b/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 diff --git a/TMSdata b/TMSdata old mode 100644 new mode 100755 diff --git a/default.com b/default.com old mode 100644 new mode 100755 diff --git a/default.key b/default.key old mode 100644 new mode 100755 diff --git a/default1.com b/default1.com old mode 100644 new mode 100755 diff --git a/nmrPredictGaus.py b/nmrPredictGaus.py old mode 100644 new mode 100755 diff --git a/nmrPredictNWChem.py b/nmrPredictNWChem.py old mode 100644 new mode 100755 diff --git a/pcs-1 b/pcs-1 old mode 100644 new mode 100755 diff --git a/pcs-2 b/pcs-2 old mode 100644 new mode 100755