From d2b0abf03d01e95471024dd643072ebdb332aa79 Mon Sep 17 00:00:00 2001 From: Fabian Bohle <49951809+fabothch@users.noreply.github.com> Date: Mon, 12 Apr 2021 18:35:43 +0200 Subject: [PATCH] devel (#12) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit • work in progress: automatic supporting information (SI) generation (with bib entries) • new printout showing the number of conformers within an energy window • symmetry number correction (for part1) because of (DFT) unoptimized geometries • the configuration file location information is saved and automatically reloaded on restart • printing command line argument call at program startup • experimental “read only” feature to fill the enso.json information by only reading the output of previous calculations (the program input arguments have to be exactly the same!) • fixed NMR shielding constant calculations and MORead • NMR new printout of min and max calculated shielding constants per investigated nuclei • OR new printout of min and max OR value and corresponding conformers • fixed balance option and wrong xtb call • improved restarting of optical rotation calculations • tutorial option has details on solvents and naming within CENSO, and information on functionals for OR is updated • improved optimization routine --- censo_qm/adf_job.py | 2 +- censo_qm/censo.py | 74 +++-- censo_qm/cfg.py | 530 +++++++++++++++++++++++++++------ censo_qm/cheapscreening.py | 133 +++++---- censo_qm/datastructure.py | 63 ++-- censo_qm/ensembledata.py | 57 ++++ censo_qm/inputhandling.py | 568 ++++++++++++++++++++++++++---------- censo_qm/nmrproperties.py | 304 +++++++++++++------ censo_qm/opticalrotation.py | 171 +++++++---- censo_qm/optimization.py | 277 +++++++++++------- censo_qm/orca_job.py | 303 +++++++++---------- censo_qm/parallel.py | 33 ++- censo_qm/prescreening.py | 229 ++++++++++----- censo_qm/qm_job.py | 412 ++++++++++---------------- censo_qm/refinement.py | 256 ++++++++++------ censo_qm/setupcenso.py | 118 +++++--- censo_qm/tm_job.py | 501 ++++++++++++++++--------------- censo_qm/tutorial.py | 44 ++- censo_qm/utilities.py | 92 +++++- setup.cfg | 2 +- 20 files changed, 2646 insertions(+), 1523 deletions(-) diff --git a/censo_qm/adf_job.py b/censo_qm/adf_job.py index 714b9d0..d798563 100644 --- a/censo_qm/adf_job.py +++ b/censo_qm/adf_job.py @@ -1,3 +1,3 @@ # adf_job.py -# for J and S calculation only \ No newline at end of file +# for J and S calculation only diff --git a/censo_qm/censo.py b/censo_qm/censo.py index 5b7b3fa..dcf8cab 100644 --- a/censo_qm/censo.py +++ b/censo_qm/censo.py @@ -62,7 +62,7 @@ def main(argv=None): sys.exit(1) toc = perf_counter() ensembledata.part_info["part0"] = toc - tic - ensembledata.previous_part_info["part0"] += ensembledata.part_info["part0"] + ensembledata.previous_part_info["part0"] += ensembledata.part_info["part0"] print(f"Ran part0 in {ensembledata.part_info['part0']:0.4f} seconds") # RUNNING PART1 @@ -87,7 +87,9 @@ def main(argv=None): toc = perf_counter() ensembledata.part_info["part1"] = toc - tic ensembledata.previous_part_info["part1"] += ensembledata.part_info["part1"] - ensembledata.previous_part_info["part1_firstsort"] += ensembledata.part_info["part1_firstsort"] + ensembledata.previous_part_info["part1_firstsort"] += ensembledata.part_info[ + "part1_firstsort" + ] print(f"Ran part1 in {ensembledata.part_info['part1']:0.4f} seconds") # RUNNING PART2 @@ -112,7 +114,9 @@ def main(argv=None): toc = perf_counter() ensembledata.part_info["part2"] = toc - tic ensembledata.previous_part_info["part2"] += ensembledata.part_info["part2"] - ensembledata.previous_part_info["part2_opt"] += ensembledata.part_info["part2_opt"] + ensembledata.previous_part_info["part2_opt"] += ensembledata.part_info[ + "part2_opt" + ] print(f"Ran part2 in {ensembledata.part_info['part2']:0.4f} seconds") # RUNNING PART3 @@ -196,6 +200,9 @@ def main(argv=None): config.provide_runinfo(), ) + if args.create_si: + config.create_SI(ensembledata) + # END of CENSO timings = 0.0 prev_timings = 0.0 @@ -208,44 +215,47 @@ def main(argv=None): tmp = [] tmp_prev = [] if config.part0: - tmp.append(ensembledata.part_info['part0']) - tmp_prev.append(ensembledata.previous_part_info['part0']) + tmp.append(ensembledata.part_info["part0"]) + tmp_prev.append(ensembledata.previous_part_info["part0"]) if config.part1: - tmp.append(ensembledata.part_info['part1']) - tmp_prev.append(ensembledata.previous_part_info['part1']) - tmp_prev.append(ensembledata.previous_part_info['part1_firstsort']) + tmp.append(ensembledata.part_info["part1"]) + tmp_prev.append(ensembledata.previous_part_info["part1"]) + tmp_prev.append(ensembledata.previous_part_info["part1_firstsort"]) if config.part2: - tmp.append(ensembledata.part_info['part2']) - tmp_prev.append(ensembledata.previous_part_info['part2']) - tmp.append(ensembledata.part_info['part2_opt']) - tmp_prev.append(ensembledata.previous_part_info['part2_opt']) + tmp.append(ensembledata.part_info["part2"]) + tmp_prev.append(ensembledata.previous_part_info["part2"]) + tmp.append(ensembledata.part_info["part2_opt"]) + tmp_prev.append(ensembledata.previous_part_info["part2_opt"]) if config.part3: - tmp.append(ensembledata.part_info['part3']) - tmp_prev.append(ensembledata.previous_part_info['part3']) + tmp.append(ensembledata.part_info["part3"]) + tmp_prev.append(ensembledata.previous_part_info["part3"]) if config.part4: - tmp.append(ensembledata.part_info['part4']) - tmp_prev.append(ensembledata.previous_part_info['part4']) + tmp.append(ensembledata.part_info["part4"]) + tmp_prev.append(ensembledata.previous_part_info["part4"]) if config.optical_rotation: - tmp.append(ensembledata.part_info['part5']) - tmp_prev.append(ensembledata.previous_part_info['part5']) - timelen = max([len(f"{float(value):.2f}") for value in tmp]) +2 - prev_timelen = max([len(f"{float(value):.2f}") for value in tmp_prev]) +2 + tmp.append(ensembledata.part_info["part5"]) + tmp_prev.append(ensembledata.previous_part_info["part5"]) + timelen = max([len(f"{float(value):.2f}") for value in tmp]) + 2 + prev_timelen = max([len(f"{float(value):.2f}") for value in tmp_prev]) + 2 if timelen < 7: timelen = 7 except Exception: timelen = 20 prev_timelen = 20 - - print(f"\n\n{'Part':20}: {'#conf':>{conflength}} {'time': >{timelen}} time (including restarts)") - print("".ljust(int(PLENGTH/1.4), "-")) - print(f"{'Input':20}: {ensembledata.nconfs_per_part['starting']:{conflength}} {'-':^{timelen+2}} {'-':^{timelen+2}}") + print( + f"\n\n{'Part':20}: {'#conf':>{conflength}} {'time': >{timelen}} time (including restarts)" + ) + print("".ljust(int(PLENGTH / 1.4), "-")) + print( + f"{'Input':20}: {ensembledata.nconfs_per_part['starting']:{conflength}} {'-':^{timelen+2}} {'-':^{timelen+2}}" + ) if config.part0: print( f"{'Part0_all':20}: {ensembledata.nconfs_per_part['part0']:{conflength}} {ensembledata.part_info['part0']:{timelen}.2f} s {ensembledata.previous_part_info['part0']:>{prev_timelen}.2f} s" ) timings += ensembledata.part_info["part0"] - prev_timings += ensembledata.previous_part_info['part0'] + prev_timings += ensembledata.previous_part_info["part0"] if config.part1: print( f"{'Part1_initial_sort':20}: {ensembledata.nconfs_per_part['part1_firstsort']:{conflength}} {ensembledata.part_info['part1_firstsort']:{timelen}.2f} s {ensembledata.previous_part_info['part1_firstsort']:>{prev_timelen}.2f} s" @@ -254,7 +264,7 @@ def main(argv=None): f"{'Part1_all':20}: {ensembledata.nconfs_per_part['part1_firstsort']:{conflength}} {ensembledata.part_info['part1']:{timelen}.2f} s {ensembledata.previous_part_info['part1']:>{prev_timelen}.2f} s" ) timings += ensembledata.part_info["part1"] - prev_timings += ensembledata.previous_part_info['part1'] + prev_timings += ensembledata.previous_part_info["part1"] if config.part2: print( f"{'Part2_opt':20}: {ensembledata.nconfs_per_part['part2_opt']:{conflength}} {ensembledata.part_info['part2_opt']:{timelen}.2f} s {ensembledata.previous_part_info['part2_opt']:>{prev_timelen}.2f} s" @@ -263,26 +273,28 @@ def main(argv=None): f"{'Part2_all':20}: {ensembledata.nconfs_per_part['part2']:{conflength}} {ensembledata.part_info['part2']:{timelen}.2f} s {ensembledata.previous_part_info['part2']:>{prev_timelen}.2f} s" ) timings += ensembledata.part_info["part2"] - prev_timings += ensembledata.previous_part_info['part2'] + prev_timings += ensembledata.previous_part_info["part2"] if config.part3: print( f"{'Part3_all':20}: {ensembledata.nconfs_per_part['part3']:{conflength}} {ensembledata.part_info['part3']:{timelen}.2f} s {ensembledata.previous_part_info['part3']:>{prev_timelen}.2f} s" ) timings += ensembledata.part_info["part3"] - prev_timings += ensembledata.previous_part_info['part3'] + prev_timings += ensembledata.previous_part_info["part3"] if config.part4: print( f"{'Part4':20}: {'':{conflength}} {ensembledata.part_info['part4']:{timelen}.2f} s {ensembledata.previous_part_info['part4']:>{prev_timelen}.2f} s" ) timings += ensembledata.part_info["part4"] - prev_timings += ensembledata.previous_part_info['part4'] + prev_timings += ensembledata.previous_part_info["part4"] if config.optical_rotation: print( f"{'Part5':20}: {'':{conflength}} {ensembledata.part_info['part5']:{timelen}.2f} s {ensembledata.previous_part_info['part5']:>{prev_timelen}.2f} s" ) timings += ensembledata.part_info["part5"] - prev_timings += ensembledata.previous_part_info['part5'] + prev_timings += ensembledata.previous_part_info["part5"] print("".ljust(int(PLENGTH / 1.4), "-")) - print(f"{'All parts':20}: {'-':>{conflength}} {timings:{timelen}.2f} s {prev_timings:{prev_timelen}.2f} s") + print( + f"{'All parts':20}: {'-':>{conflength}} {timings:{timelen}.2f} s {prev_timings:{prev_timelen}.2f} s" + ) print("\nCENSO all done!") return 0 diff --git a/censo_qm/cfg.py b/censo_qm/cfg.py index 69023f9..258a6d8 100644 --- a/censo_qm/cfg.py +++ b/censo_qm/cfg.py @@ -5,7 +5,7 @@ """ import os -__version__ = "1.0.8" +__version__ = "1.0.9" DESCR = f""" ______________________________________________________________ @@ -40,7 +40,7 @@ R = 1.987203585e-03 # kcal/(mol*K) AU2KCAL = 627.50947428 BOHR2ANG = 0.52917721067 -WARNLEN = max([len(i) for i in ['WARNING:', 'ERROR:', 'INFORMATION:']])+1 +WARNLEN = max([len(i) for i in ["WARNING:", "ERROR:", "INFORMATION:"]]) + 1 # definitions: composite_method_basis = { @@ -53,15 +53,16 @@ "r2scan-3c": "def2-mTZVPP", } composite_dfa = tuple(composite_method_basis.keys()) -gga_dfa = ("tpss", "pbe", "kt2", "b97-d3") +gga_dfa = ("tpss", "pbe", "kt1", "kt2", "b97-d3", "b97-d") hybrid_dfa = ( "pbe0", "pw6b95", "wb97x-d3", - 'wb97x-d3bj', + "wb97x-d3bj", "wb97x-v", "cam-b3lyp", - "b3-lyp", + "b3-lyp", # tm + "b3lyp", # orca "pbeh-3c", "m06x", "bh-lyp", @@ -69,11 +70,12 @@ ) dh_dfa = ("dsd-blyp",) disp_already_included_in_func = composite_dfa + ( - 'b97-d3', - 'wb97x-d3', - 'wb97x-d3bj', + "b97-d3", # orca + "b97-d", # tm + "wb97x-d3", + "wb97x-d3bj", "wb97x-v", - ) +) # program paths: external_paths = {} @@ -330,7 +332,8 @@ }, } -class NmrRef(): + +class NmrRef: """nmrreference data in the format: [reference-molecule][func-geometry][funcS][solvent] h_tm_shieldings c_tm_shieldings @@ -343,6 +346,7 @@ class NmrRef(): si_orca_shieldings p_orca_shieldings """ + h_tm_shieldings = { "TMS": { "r2scan-3c": { @@ -356,7 +360,7 @@ class NmrRef(): "h2o": 31.51254989166667, "methanol": 31.520660391666667, "thf": 31.509407100000004, - "toluene": 31.489292383333336 + "toluene": 31.489292383333336, }, "pbeh-3c": { "gas": 32.174366075, @@ -368,7 +372,7 @@ class NmrRef(): "h2o": 32.142147875000006, "methanol": 32.14963604166667, "thf": 32.13958754166667, - "toluene": 32.121203466666664 + "toluene": 32.121203466666664, }, "b97-3c": { "gas": 32.174913100000005, @@ -380,7 +384,7 @@ class NmrRef(): "h2o": 32.14978633333334, "methanol": 32.15600726666666, "thf": 32.14655236666667, - "toluene": 32.128767625 + "toluene": 32.128767625, }, "tpss": { "gas": 31.898731900000012, @@ -392,7 +396,7 @@ class NmrRef(): "h2o": 31.87253611666667, "methanol": 31.87929998333333, "thf": 31.870970700000004, - "toluene": 31.851654483333338 + "toluene": 31.851654483333338, }, "pbe0": { "gas": 31.66924141666666, @@ -404,7 +408,7 @@ class NmrRef(): "h2o": 31.637786341666665, "methanol": 31.644428750000003, "thf": 31.633562058333336, - "toluene": 31.6153103 + "toluene": 31.6153103, }, "kt2": { "gas": 31.667160058333337, @@ -416,8 +420,8 @@ class NmrRef(): "h2o": 31.635154825, "methanol": 31.642592158333333, "thf": 31.631994574999997, - "toluene": 31.612610516666663 - } + "toluene": 31.612610516666663, + }, }, "pbeh-3c": { "tpss": { @@ -774,7 +778,7 @@ class NmrRef(): "h2o": 194.73630830000002, "methanol": 194.47241222500003, "thf": 194.5933908, - "toluene": 194.76792010000003 + "toluene": 194.76792010000003, }, "pbeh-3c": { "gas": 197.51118245, @@ -786,7 +790,7 @@ class NmrRef(): "h2o": 198.22405855, "methanol": 197.911550875, "thf": 198.02789195, - "toluene": 198.28027294999998 + "toluene": 198.28027294999998, }, "b97-3c": { "gas": 184.57099812500002, @@ -798,7 +802,7 @@ class NmrRef(): "h2o": 185.317873325, "methanol": 185.01346424999997, "thf": 185.100088675, - "toluene": 185.3094498 + "toluene": 185.3094498, }, "tpss": { "gas": 185.706348925, @@ -810,7 +814,7 @@ class NmrRef(): "h2o": 186.333736775, "methanol": 186.07168465, "thf": 186.13180945, - "toluene": 186.32699157500002 + "toluene": 186.32699157500002, }, "pbe0": { "gas": 187.894285475, @@ -822,20 +826,20 @@ class NmrRef(): "h2o": 188.5966455, "methanol": 188.355007975, "thf": 188.4689729, - "toluene": 188.594034275 + "toluene": 188.594034275, }, "kt2": { - "gas": 189.78494644999998, - "acetone": 190.329502875, - "chcl3": 190.204013175, - "acetonitrile": 190.397052075, - "ch2cl2": 190.06665505, - "dmso": 190.4107424, - "h2o": 190.40970589999998, - "methanol": 190.188391875, - "thf": 190.2872299, - "toluene": 190.41299607500002 - } + "gas": 189.78494644999998, + "acetone": 190.329502875, + "chcl3": 190.204013175, + "acetonitrile": 190.397052075, + "ch2cl2": 190.06665505, + "dmso": 190.4107424, + "h2o": 190.40970589999998, + "methanol": 190.188391875, + "thf": 190.2872299, + "toluene": 190.41299607500002, + }, }, "pbeh-3c": { "tpss": { @@ -2220,7 +2224,7 @@ class NmrRef(): "h2o": 357.9281381, "methanol": 357.9240386, "thf": 357.8718386, - "toluene": 357.8388333 + "toluene": 357.8388333, }, "pbeh-3c": { "gas": 425.427401, @@ -2232,7 +2236,7 @@ class NmrRef(): "h2o": 425.1994123, "methanol": 425.1168623, "thf": 425.1081863, - "toluene": 425.1184486 + "toluene": 425.1184486, }, "b97-3c": { "gas": 352.2983555, @@ -2244,7 +2248,7 @@ class NmrRef(): "h2o": 351.7346671, "methanol": 351.6729095, "thf": 351.6321795, - "toluene": 351.648296 + "toluene": 351.648296, }, "tpss": { "gas": 334.0278062, @@ -2256,7 +2260,7 @@ class NmrRef(): "h2o": 333.555991, "methanol": 333.5402484, "thf": 333.5146132, - "toluene": 333.4768581 + "toluene": 333.4768581, }, "pbe0": { "gas": 331.8236884, @@ -2268,7 +2272,7 @@ class NmrRef(): "h2o": 331.3437693, "methanol": 331.3449225, "thf": 331.3305993, - "toluene": 331.3160006 + "toluene": 331.3160006, }, "kt2": { "gas": 340.923509, @@ -2280,8 +2284,8 @@ class NmrRef(): "h2o": 340.4180652, "methanol": 340.3808603, "thf": 340.348664, - "toluene": 340.3406705 - } + "toluene": 340.3406705, + }, }, "pbeh-3c": { "tpss": { @@ -2625,6 +2629,7 @@ class NmrRef(): }, } } + def NMRRef_to_dict(self): """Convert NMRRef data to a dict object""" dict_ret = dict( @@ -2638,56 +2643,415 @@ def NMRRef_to_dict(self): f_orca_shieldings=self.f_orca_shieldings, si_orca_shieldings=self.si_orca_shieldings, p_orca_shieldings=self.p_orca_shieldings, - ) + ) return dict_ret def dict_to_NMRRef(self, dictionary): """Convert dict object to NMRRef data """ NmrRef_object = NmrRef() - NmrRef_object.h_tm_shieldings = dictionary.get('h_tm_shieldings', NmrRef_object.h_tm_shieldings) - NmrRef_object.c_tm_shieldings = dictionary.get('c_tm_shieldings', NmrRef_object.c_tm_shieldings) - NmrRef_object.f_tm_shieldings = dictionary.get('f_tm_shieldings', NmrRef_object.f_tm_shieldings) - NmrRef_object.si_tm_shieldings = dictionary.get('si_tm_shieldings', NmrRef_object.si_tm_shieldings) - NmrRef_object.p_tm_shieldings = dictionary.get('p_tm_shieldings', NmrRef_object.p_tm_shieldings) - NmrRef_object.h_orca_shieldings = dictionary.get('h_orca_shieldings', NmrRef_object.h_orca_shieldings) - NmrRef_object.c_orca_shieldings = dictionary.get('c_orca_shieldings', NmrRef_object.c_orca_shieldings) - NmrRef_object.f_orca_shieldings = dictionary.get('f_orca_shieldings', NmrRef_object.f_orca_shieldings) - NmrRef_object.si_orca_shieldings = dictionary.get('si_orca_shieldings', NmrRef_object.si_orca_shieldings) - NmrRef_object.p_orca_shieldings = dictionary.get('p_orca_shieldings', NmrRef_object.p_orca_shieldings) + NmrRef_object.h_tm_shieldings = dictionary.get( + "h_tm_shieldings", NmrRef_object.h_tm_shieldings + ) + NmrRef_object.c_tm_shieldings = dictionary.get( + "c_tm_shieldings", NmrRef_object.c_tm_shieldings + ) + NmrRef_object.f_tm_shieldings = dictionary.get( + "f_tm_shieldings", NmrRef_object.f_tm_shieldings + ) + NmrRef_object.si_tm_shieldings = dictionary.get( + "si_tm_shieldings", NmrRef_object.si_tm_shieldings + ) + NmrRef_object.p_tm_shieldings = dictionary.get( + "p_tm_shieldings", NmrRef_object.p_tm_shieldings + ) + NmrRef_object.h_orca_shieldings = dictionary.get( + "h_orca_shieldings", NmrRef_object.h_orca_shieldings + ) + NmrRef_object.c_orca_shieldings = dictionary.get( + "c_orca_shieldings", NmrRef_object.c_orca_shieldings + ) + NmrRef_object.f_orca_shieldings = dictionary.get( + "f_orca_shieldings", NmrRef_object.f_orca_shieldings + ) + NmrRef_object.si_orca_shieldings = dictionary.get( + "si_orca_shieldings", NmrRef_object.si_orca_shieldings + ) + NmrRef_object.p_orca_shieldings = dictionary.get( + "p_orca_shieldings", NmrRef_object.p_orca_shieldings + ) return NmrRef_object + # rotational entropy from symmetry -#https://cccbdb.nist.gov/thermo.asp +# https://cccbdb.nist.gov/thermo.asp rot_sym_num = { - 'c1':1, - 'ci':1, - 'cs':1, - 'c2':2, - 'c3':3, - 'c4':4, - 'c5':5, - 'c6':6, - 'c7':7, - 'c8':8, - 'c9':9, - 'c10':10, - 'c11':11, - 's4':2, - 's6':3, - 's8':4, - 'd2':4, - 'd3':6, - 'd4':8, - 'd5':10, - 'd6':12, - 'd7':14, - 'd8':16, - 'd9':18, - 'd10':20, - 't':12, - 'th': 12, - 'td': 12, - 'o': 24, - 'oh':24, - 'ih':60 + "c1": 1, + "ci": 1, + "cs": 1, + "c2": 2, + "c3": 3, + "c4": 4, + "c5": 5, + "c6": 6, + "c7": 7, + "c8": 8, + "c9": 9, + "c10": 10, + "c11": 11, + "s4": 2, + "s6": 3, + "s8": 4, + "d2": 4, + "d3": 6, + "d4": 8, + "d5": 10, + "d6": 12, + "d7": 14, + "d8": 16, + "d9": 18, + "d10": 20, + "t": 12, + "th": 12, + "td": 12, + "o": 24, + "oh": 24, + "ih": 60, +} + + +si_bib = { + "tm": [ + r"@misc{TURBOMOLE,", + r" title = {{TURBOMOLE V7.5 2020}, a development of {University of Karlsruhe} and", + r" {Forschungszentrum Karlsruhe GmbH}, 1989-2007,", + r" {TURBOMOLE GmbH}, since 2007; available from \\", + r" {\tt https://www.turbomole.org}.}", + r"}", + r"@Article{TURBOMOLE.2020", + r" author = {Balasubramani, Sree Ganesh and Chen, Guo P. ", + r" and Coriani, Sonia and Diedenhofen, Michael and ", + r" Frank, Marius S. and Franzke, Yannick J. and ", + r" Furche, Filipp and Grotjahn, Robin and Harding, Michael E. ", + r" and H\"attig, Christof and Hellweg, Arnim and ", + r" Helmich-Paris, Benjamin and Holzer, Christof and Huniar, Uwe", + r" and Kaupp, Martin and Marefat Khah, Alireza ", + r" and Karbalaei Khani, Sarah and M\"uller, Thomas and Mack, Fabian", + r" and Nguyen, Brian D. and Parker, Shane M. and Perlt, Eva ", + r" and Rappoport, Dmitrij and Reiter, Kevin and Roy, Saswata and", + r" R\"uckert, Matthias and Schmitz, Gunnar and Sierka, Marek", + r" and Tapavicza, Enrico and Tew, David P. and van W\"ullen, Christoph", + r" and Voora, Vamsee K. and Weigend, Florian and", + r" Wody{\’n}ski, Artur and Yu, Jason M.},", + r" title = {TURBOMOLE: Modular program suite for \textit{ab initio}", + r" quantum-chemical and condensed-matter simulations},", + r" journal = {J. Chem. Phys.},", + r" volume = {152},", + r" issue = {18},", + r" pages = {184107},", + r" year = {2020},", + r" url = { https://doi.org/10.1063/5.0004635},", + r" DOI = {10.1063/5.0004635}", + r"}", + ], + "orca": [ + r"@article{ORCA_generic,", + r" author = {Neese, Frank},", + r" title = {The ORCA program system},", + r" journal = {WIREs Computational Molecular Science},", + r" volume = {2},", + r" number = {1},", + r" pages = {73-78},", + r" doi = {https://doi.org/10.1002/wcms.81},", + r" url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/wcms.81},", + r" eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.81},", + r" year = {2012}", + r"}", + r"@article{ORCA4.0,", + r" author = {Neese, Frank},", + r" title = {Software update: the ORCA program system, version 4.0},", + r" journal = {WIREs Computational Molecular Science},", + r" volume = {8},", + r" number = {1},", + r" pages = {e1327},", + r" doi = {https://doi.org/10.1002/wcms.1327},", + r" url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/wcms.1327},", + r" eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1327},", + r" year = {2018}", + r"}", + ], + "cosmors": [ + r"@article{Klamt1995,", + r" author = {Klamt, Andreas},", + r" doi = {10.1021/j100007a062},", + r" journal = {The Journal of Physical Chemistry},", + r" number = {7},", + r" pages = {2224--2235},", + r" title = {{Conductor-like Screening Model for Real Solvents: A New Approach to the Quantitative Calculation of Solvation Phenomena}},", + r" url = {https://pubs.acs.org/doi/abs/10.1021/j100007a062},", + r" volume = {99},", + r" year = {1995}", + r"}", + r"@article{Klamt1998,", + r" author = {Klamt, Andreas and Jonas, Volker and B{\"{u}}rger, Thorsten and Lohrenz, John C. W.},", + r" doi = {10.1021/jp980017s},", + r" journal = {The Journal of Physical Chemistry A},", + r" number = {26},", + r" pages = {5074--5085},", + r" title = {{Refinement and Parametrization of COSMO-RS}},", + r" url = {https://pubs.acs.org/doi/10.1021/jp980017s},", + r" volume = {102},", + r" year = {1998}", + r"}", + r"@article{Eckert2002,", + r" author = {Eckert, Frank and Klamt, Andreas},", + r" doi = {10.1002/aic.690480220},", + r" journal = {AIChE Journal},", + r" number = {2},", + r" pages = {369--385},", + r" title = {{Fast solvent screening via quantum chemistry: COSMO-RS approach}},", + r" url = {http://doi.wiley.com/10.1002/aic.690480220},", + r" volume = {48},", + r" year = {2002}", + r"}", + r"@misc{COSMOtherm,", + r" title = {COSMOtherm, Release 19; COSMOlogic GmbH & Co. KG, {\tt http://www.cosmologic.de}.}", + r"}", + ], + "xtb": [ + r"@article{xtb_generic,", + r" title={Extended tight‐binding quantum chemistry methods},", + r" author={Bannwarth, Christoph and Caldeweyher, Eike and Ehlert, Sebastian and Hansen, Andreas and Pracht, Philipp and Seibert, Jakob and Spicher, Spicher and Grimme, Stefan},", + r" journal={{WIREs} Comput{.} Mol{.} Sci{.}},", + r" volume = {11},", + r" year={2020},", + r" pages={e01493},", + r" doi={10.1002/wcms.1493},", + r" url={https://dx.doi.org/10.1002/wcms.1493}", + r"}", + r"@article{GFN2,", + r" title={GFN2-xTB—An accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions},", + r" author={Bannwarth, Christoph and Ehlert, Sebastian and Grimme, Stefan},", + r" journal={J{.} Chem{.} Theory Comput{.}},", + r" volume={15},", + r" number={3},", + r" pages={1652--1671},", + r" year={2019},", + r" doi={10.1021/acs.jctc.8b01176},", + r" url={https://dx.doi.org/10.1021/acs.jctc.8b01176},", + r"}", + r"@article{GFN1,", + r" title={A robust and accurate tight-binding quantum chemical method for structures, vibrational frequencies, and noncovalent interactions of large molecular systems parametrized for all spd-block elements (Z=1--86)},", + r" author={Grimme, Stefan and Bannwarth, Christoph and Shushkov, Philip},", + r" journal={J{.} Chem{.} Theory Comput{.}},", + r" volume={13},", + r" number={5},", + r" pages={1989--2009},", + r" year={2017},", + r" doi={10.1021/acs.jctc.7b00118},", + r" url={https://dx.doi.org/10.1021/acs.jctc.7b00118},", + r"}", + ], + "censo": [ + r"@article{CENSO", + r" author = {Grimme, Stefan and Bohle, Fabian and Hansen, Andreas and Pracht, Philipp and Spicher, Sebastian and Stahn, Marcel},", + r" title = {Efficient Quantum Chemical Calculation of Structure Ensembles and Free Energies for Nonrigid Molecules},", + r" journal = {The Journal of Physical Chemistry A},", + r" volume = {0},", + r" number = {0},", + r" pages = {null},", + r" year = {0},", + r" doi = {10.1021/acs.jpca.1c00971},", + r" note ={PMID: 33688730},", + r" URL = {https://doi.org/10.1021/acs.jpca.1c00971},", + r" eprint = {https://doi.org/10.1021/acs.jpca.1c00971}", + r"}", + ], + # selected functionals: + "r2scan-3c": [ + r"@article{r2scan-3c,", + r" author = {Grimme, Stefan and Hansen, Andreas and Ehlert, Sebastian and Mewes, Jan-Michael},", + r" doi = {10.1063/5.0040021},", + r" journal = {The Journal of Chemical Physics},", + r" number = {6},", + r" pages = {064103},", + r" title = {{r 2 SCAN-3c: A “Swiss army knife” composite electronic-structure method}},", + r" url = {https://doi.org/10.1063/5.0040021},", + r" volume = {154},", + r" year = {2021}", + r"}", + ], + "pbeh-3c": [ + r"@article{Grimme2015,", + r" author = {Grimme, Stefan and Brandenburg, Jan Gerit and Bannwarth, Christoph and Hansen, Andreas},", + r" doi = {10.1063/1.4927476},", + r" journal = {The Journal of Chemical Physics},", + r" number = {5},", + r" pages = {054107},", + r" pmid = {26254642},", + r" title = {{Consistent structures and interactions by density functional theory with small atomic orbital basis sets}},", + r" url = {http://scitation.aip.org/content/aip/journal/jcp/143/5/10.1063/1.4927476},", + r" volume = {143},", + r" year = {2015}", + r"}", + ], + "b97-3c": [ + r"@article{b97-3c,", + r" author = {Brandenburg, Jan Gerit and Bannwarth, Christoph and Hansen, Andreas and Grimme, Stefan},", + r" doi = {10.1063/1.5012601},", + r" journal = {The Journal of Chemical Physics},", + r" number = {6},", + r" pages = {064104},", + r" title = {{B97-3c: A revised low-cost variant of the B97-D density functional method}},", + r" url = {http://aip.scitation.org/doi/10.1063/1.5012601},", + r" volume = {148},", + r" year = {2018}", + r"}", + ], + "pbe0": [ + r"@article{pbe0_one,", + r" author = {Perdew, John P. and Burke, Kieron and Ernzerhof, Matthias},", + r" doi = {10.1103/PhysRevLett.77.3865},", + r" journal = {Physical Review Letters},", + r" number = {18},", + r" pages = {3865--3868},", + r" title = {{Generalized Gradient Approximation Made Simple}},", + r" url = {https://link.aps.org/doi/10.1103/PhysRevLett.77.3865},", + r" volume = {77},", + r" year = {1996}", + r"}", + r"@article{pbe0_one_erratum,", + r" title = {Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)]},", + r" author = {Perdew, John P. and Burke, Kieron and Ernzerhof, Matthias},", + r" journal = {Phys. Rev. Lett.},", + r" volume = {78},", + r" issue = {7},", + r" pages = {1396--1396},", + r" numpages = {0},", + r" year = {1997},", + r" month = {Feb},", + r" publisher = {American Physical Society},", + r" doi = {10.1103/PhysRevLett.78.1396},", + r" url = {https://link.aps.org/doi/10.1103/PhysRevLett.78.1396}", + r"}", + r"@article{PBE0_two,", + r" author = {Adamo,Carlo and Barone,Vincenzo},", + r" title = {Toward reliable density functional methods without adjustable parameters: The PBE0 model},", + r" journal = {The Journal of Chemical Physics},", + r" volume = {110},", + r" number = {13},", + r" pages = {6158-6170},", + r" year = {1999},", + r" doi = {10.1063/1.478522},", + r" url = {https://doi.org/10.1063/1.478522},", + r" eprint = {https://doi.org/10.1063/1.478522}", + r"}", + ], + "b3-lyp": [ + r"@article{b3lyp_one,", + r" author = {Becke,Axel D. },", + r" title = {A new mixing of Hartree–Fock and local density‐functional theories},", + r" journal = {The Journal of Chemical Physics},", + r" volume = {98},", + r" number = {2},", + r" pages = {1372-1377},", + r" year = {1993},", + r" doi = {10.1063/1.464304},", + r" url = {https://doi.org/10.1063/1.464304},", + r" eprint = {https://doi.org/10.1063/1.464304}", + r"}", + r"@article{b3lyp_two,", + r" author = {Stephens, P. J. and Devlin, F. J. and Chabalowski, C. F. and Frisch, M. J.},", + r" title = {Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields},", + r" journal = {The Journal of Physical Chemistry},", + r" volume = {98},", + r" number = {45},", + r" pages = {11623-11627},", + r" year = {1994},", + r" doi = {10.1021/j100096a001},", + r" url = {https://doi.org/10.1021/j100096a001},", + r" eprint = {https://doi.org/10.1021/j100096a001}", + r"}", + ], + "def2_basis": [ + r"@article{Weigend2005,", + r" author = {Weigend, Florian and Ahlrichs, Reinhart},", + r" doi = {10.1039/b508541a},", + r" journal = {Physical Chemistry Chemical Physics},", + r" number = {18},", + r" pages = {3297},", + r" title = {{Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy}},", + r" url = {http://xlink.rsc.org/?DOI=b508541a},", + r" volume = {7},", + r" year = {2005}", + r"}", + ], + "def2_auxbasis": [ + r"@article{Eichkorn1997,", + r"author = {Eichkorn, Karin and Weigend, Florian and Treutler, Oliver and Ahlrichs, Reinhart},", + r"doi = {10.1007/s002140050244},", + r"journal = {Theoretical Chemistry Accounts: Theory, Computation, and Modeling (Theoretica Chimica Acta)},", + r"number = {1-4},", + r"pages = {119--124},", + r"publisher = {Springer-Verlag},", + r"title = {{Auxiliary basis sets for main row atoms and transition metals and their use to approximate Coulomb potentials}},", + r"url = {http://link.springer.com/10.1007/s002140050244},", + r"volume = {97},", + r"year = {1997}", + r"}", + r"@article{Weigend2006,", + r" author = {Weigend, Florian},", + r" doi = {10.1039/b515623h},", + r" journal = {Physical Chemistry Chemical Physics},", + r" number = {9},", + r" pages = {1057},", + r" title = {{Accurate Coulomb-fitting basis sets for H to Rn}},", + r" url = {http://xlink.rsc.org/?DOI=b515623h},", + r" volume = {8},", + r" year = {2006}", + r"}", + ], + "ri": [ + r"@article{RI,", + r"title = {Integral approximations for LCAO-SCF calculations},", + r"journal = {Chemical Physics Letters},", + r"volume = {213},", + r"number = {5},", + r"pages = {514-518},", + r"year = {1993},", + r"issn = {0009-2614},", + r"doi = {https://doi.org/10.1016/0009-2614(93)89151-7},", + r"url = {https://www.sciencedirect.com/science/article/pii/0009261493891517},", + r"author = {O. Vahtras and J. Almlöf and M.W. Feyereisen},", + r"}", + ], + "sph": [ + r"@article{SPH,", + r" author = {Spicher, Sebastian and Grimme, Stefan},", + r" doi = {10.1021/acs.jctc.0c01306},", + r" journal = {Journal of Chemical Theory and Computation},", + r" number = {3},", + r" pages = {1701--1714},", + r" pmid = {33554604},", + r" title = {{Single-Point Hessian Calculations for Improved Vibrational Frequencies and Rigid-Rotor-Harmonic-Oscillator Thermodynamics}},", + r" url = {https://pubs.acs.org/doi/10.1021/acs.jctc.0c01306},", + r" volume = {17},", + r" year = {2021}", + r"}", + ], +} + +# qm_prepinfo: grid and scfconv settings for ORCA and TM +qm_prepinfo = { + "orca": { + "low": ["grid4 nofinalgrid", "loosescf"], + "low+": ["grid4 nofinalgrid", "scfconv6"], + "high": ["grid4 nofinalgrid", "scfconv7"], + "high+": ["grid5 nofinalgrid", "scfconv7"], + }, + "tm": { + "low": ["-grid", "m3", "-scfconv", "6"], + "low+": ["-grid", "m4", "-scfconv", "6"], + "high": ["-grid", "m4", "-scfconv", "7"], + "high+": ["-grid", "m5", "-scfconv", "7"], + }, } diff --git a/censo_qm/cheapscreening.py b/censo_qm/cheapscreening.py index 09ed79c..87ffc85 100644 --- a/censo_qm/cheapscreening.py +++ b/censo_qm/cheapscreening.py @@ -5,7 +5,7 @@ import os import sys from multiprocessing import JoinableQueue as Queue -from .cfg import PLENGTH, DIGILEN, AU2KCAL, WARNLEN +from .cfg import PLENGTH, DIGILEN, AU2KCAL, WARNLEN, qm_prepinfo from .parallel import run_in_parallel from .orca_job import OrcaJob from .tm_job import TmJob @@ -20,6 +20,7 @@ print, print_errors, calc_boltzmannweights, + conf_in_interval, ) @@ -52,16 +53,16 @@ def part0(config, conformers, ensembledata): max_len_digilen = 0 for item in info: - if item[0] == 'justprint': + if item[0] == "justprint": if "short-notation" in item[1]: - tmp = len(item[1]) -len('short-notation:') + tmp = len(item[1]) - len("short-notation:") else: tmp = len(item[1]) else: tmp = len(item[1]) if tmp > max_len_digilen: max_len_digilen = tmp - max_len_digilen +=1 + max_len_digilen += 1 if max_len_digilen < DIGILEN: max_len_digilen = DIGILEN @@ -147,30 +148,10 @@ def part0(config, conformers, ensembledata): "energy": 0.0, "energy2": 0.0, "success": False, - "xtb_driver_path": config.external_paths["xtbpath"], + "temperature": config.temperature, + "onlyread": config.onlyread, } - tmp_disp = "" - if config.prog == "tm": - instruction["prepinfo"] = ["clear", "-grid", "1", "-scfconv", "5", "DOGCP"] - if config.func0 == "b97-d": - instruction["prepinfo"].append("-zero") - tmp_disp = "D3(0)" - - elif config.prog == "orca": - instruction["progpath"] = config.external_paths["orcapath"] - instruction["prepinfo"] = ["low", "DOGCP"] - - instruction["method"], instruction["method2"], = config.get_method_name( - instruction["jobtype"], - func=instruction["func"], - basis=instruction["basis"], - sm=instruction["sm"], - solvent=instruction["solvent"], - prog=config.prog, - disp=tmp_disp, - gfn_version=instruction["gfn_version"], - ) elif config.solvent == "gas": instruction = { "jobtype": "sp", @@ -187,29 +168,37 @@ def part0(config, conformers, ensembledata): "energy": 0.0, "energy2": 0.0, "success": False, + "onlyread": config.onlyread, } - tmp_disp = "" - if config.prog == "tm": - instruction["prepinfo"] = ["clear", "-grid", "1", "-scfconv", "5", "DOGCP"] - if config.func0 == "b97-d": - instruction["prepinfo"].append("-zero") - tmp_disp = "D3(0)" - - elif config.prog == "orca": - instruction["progpath"] = config.external_paths["orcapath"] - instruction["prepinfo"] = ["low", "DOGCP"] + tmp_disp = "" + if config.prog == "tm": + instruction["prepinfo"] = ["clear", "-grid", "1", "-scfconv", "5", "DOGCP"] + ensembledata.si["part0"]["Energy_settings"] = "scfconv 5, grid 1" + if config.func0 == "b97-d": + instruction["prepinfo"].append("-zero") + tmp_disp = "D3(0)" + ensembledata.si["part0"][ + "Energy_settings" + ] += " using D3(0) instead of D3(BJ)" - instruction["method"], instruction["method2"], = config.get_method_name( - instruction["jobtype"], - func=instruction["func"], - basis=instruction["basis"], - sm=instruction["sm"], - solvent=instruction["solvent"], - prog=config.prog, - disp=tmp_disp, + elif config.prog == "orca": + instruction["prepinfo"] = ["low", "DOGCP"] + ensembledata.si["part0"]["Energy_settings"] = " ".join( + qm_prepinfo["orca"][instruction["prepinfo"][0]] ) + instruction["method"], instruction["method2"], = config.get_method_name( + instruction["jobtype"], + func=instruction["func"], + basis=instruction["basis"], + sm=instruction["sm"], + solvent=instruction["solvent"], + prog=config.prog, + disp=tmp_disp, + gfn_version=instruction.get("gfn_version", ""), + ) + name = "efficient gas-phase single-point" # folder = "part0_sp" check = {True: "was successful", False: "FAILED"} @@ -235,7 +224,7 @@ def part0(config, conformers, ensembledata): calculate, instruction, config.balance, - folder + folder, ) for conf in list(calculate): @@ -259,7 +248,9 @@ def part0(config, conformers, ensembledata): conf.cheap_prescreening_sp_info["info"] = "calculated" conf.cheap_prescreening_sp_info["method"] = conf.job["method"] conf.cheap_prescreening_gsolv_info["energy"] = conf.job["energy2"] - conf.cheap_prescreening_gsolv_info["range"] = {config.temperature: conf.job['energy2'],} + conf.cheap_prescreening_gsolv_info["range"] = { + config.temperature: conf.job["energy2"] + } conf.cheap_prescreening_gsolv_info["info"] = "calculated" conf.cheap_prescreening_gsolv_info["method"] = conf.job["method2"] conf.cheap_prescreening_gsolv_info["gas-energy"] = conf.job[ @@ -289,6 +280,7 @@ def part0(config, conformers, ensembledata): f"{'ERROR:':{WARNLEN}}UNEXPECTED BEHAVIOUR: {conf.job['success']} {conf.job['jobtype']}", save_errors, ) + # save current data to jsonfile config.write_json( config.cwd, @@ -301,6 +293,28 @@ def part0(config, conformers, ensembledata): check_tasks(calculate, config.check) else: print("No conformers are considered additionally.") + + # create data for possible SI generation:----------------------------------- + # Energy: + ensembledata.si["part0"]["Energy"] = instruction["method"] + if "DOGCP" in instruction["prepinfo"]: + ensembledata.si["part0"]["Energy"] += " + GCP" + # Energy_settings are set above! + # GmRRHO: + ensembledata.si["part0"]["G_mRRHO"] = "not included" + # Solvation: + if config.solvent == "gas": + ensembledata.si["part0"]["G_solv"] = "gas-phase" + else: + ensembledata.si["part0"]["G_solv"] = instruction["method2"] + # Geometry: + ensembledata.si["part0"]["Geometry"] = "GFNn-xTB (input geometry)" + # QM-CODE: + ensembledata.si["part0"]["main QM code"] = str(config.prog).upper() + # Threshold: + ensembledata.si["part0"]["Threshold"] = f"{config.part0_threshold} kcal/mol" + # END SI generation -------------------------------------------------------- + if prev_calculated: # adding conformers calculated before: for conf in list(prev_calculated): @@ -351,8 +365,8 @@ def part0(config, conformers, ensembledata): solv=solvation, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) try: minfree = min([i.free_energy for i in calculate if i is not None]) except ValueError: @@ -389,7 +403,9 @@ def part0(config, conformers, ensembledata): try: maxreldft = max([i.rel_free_energy for i in calculate if i is not None]) except ValueError: - print_errors(f"{'ERROR:':{WARNLEN}}No conformer left or error in maxreldft!", save_errors) + print_errors( + f"{'ERROR:':{WARNLEN}}No conformer left or error in maxreldft!", save_errors + ) # print sorting columncall = [ lambda conf: "CONF" + str(getattr(conf, "id")), @@ -460,6 +476,9 @@ def part0(config, conformers, ensembledata): ) print("".ljust(int(PLENGTH), "-")) # -------------------------------------------------------------------------- + calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) + conf_in_interval(calculate, full_free_energy=False) + # -------------------------------------------------------------------------- # write to enso.json config.write_json( @@ -489,7 +508,9 @@ def part0(config, conformers, ensembledata): ) print_block(["CONF" + str(i.id) for i in calculate]) else: - print_errors(f"{'ERROR:':{WARNLEN}}There are no more conformers left!", save_errors) + print_errors( + f"{'ERROR:':{WARNLEN}}There are no more conformers left!", save_errors + ) else: for conf in list(calculate): conf.part_info["part0"] = "passed" @@ -543,17 +564,11 @@ def part0(config, conformers, ensembledata): ) if save_errors: - print( - "\n***---------------------------------------------------------***" - ) - print( - "Printing most relevant errors again, just for user convenience:" - ) + print("\n***---------------------------------------------------------***") + print("Printing most relevant errors again, just for user convenience:") for _ in list(save_errors): print(save_errors.pop()) - print( - "***---------------------------------------------------------***" - ) + print("***---------------------------------------------------------***") tmp = int((PLENGTH - len("END of Part0")) / 2) print("\n" + "".ljust(tmp, ">") + "END of Part0" + "".rjust(tmp, "<")) diff --git a/censo_qm/datastructure.py b/censo_qm/datastructure.py index 3478579..9baa570 100755 --- a/censo_qm/datastructure.py +++ b/censo_qm/datastructure.py @@ -4,7 +4,7 @@ """ from collections import OrderedDict from .utilities import print -from .cfg import R, AU2KCAL +from .cfg import R, AU2KCAL, rot_sym_num from math import log @@ -72,6 +72,9 @@ def __init__( "fuzzythr": 0.0, "rmsd": None, "prev_methods": None, + "sym": "c1", + "linear": False, + "symnum": 1, }, lowlevel_grrho_info={ "energy": None, @@ -80,6 +83,9 @@ def __init__( "method": None, "rmsd": None, "prev_methods": None, + "sym": "c1", + "linear": False, + "symnum": 1, }, lowlevel_hrrho_info={ "energy": None, @@ -96,6 +102,9 @@ def __init__( "method": None, "rmsd": None, "prev_methods": None, + "sym": "c1", + "linear": False, + "symnum": 1, }, highlevel_hrrho_info={ "energy": None, @@ -238,7 +247,7 @@ def __init__( temperature_info["range"] = [] if xtb_energy is None: xtb_energy = 100.0 - #if xtb_energy_unbiased is None: + # if xtb_energy_unbiased is None: # xtb_energy_unbiased = 100.0 if xtb_free_energy is None: xtb_free_energy = 100.0 @@ -512,7 +521,7 @@ def provide_runinfo(self): # Calculate free energy for molecule either at normal temperature, # or if the temperature is not None from the range of temperatures. # if out=False free energy is written to self.free_energy - # if out=True free energy is simply returned + # if out=True free energy is simply returned # """ # if t is None: # try: @@ -559,7 +568,9 @@ def provide_runinfo(self): # else: # return f - def calc_free_energy(self, e=None, solv=None, rrho=None, t=None, out=False, consider_sym=None): + def calc_free_energy( + self, e=None, solv=None, rrho=None, t=None, out=False, consider_sym=None + ): """ Calculate free energy for molecule either at normal temperature, or if the temperature is not None from the range of temperatures. @@ -567,7 +578,9 @@ def calc_free_energy(self, e=None, solv=None, rrho=None, t=None, out=False, cons if out=True free energy is simply returned """ if t is None: - print("xxxxxxxxxxxxx No temperature provided in calc_free_energy xxxxxxxxxxxxxxx") + print( + "xxxxxxxxxxxxx No temperature provided in calc_free_energy xxxxxxxxxxxxxxx" + ) try: f = 0.0 if e is not None: @@ -577,16 +590,17 @@ def calc_free_energy(self, e=None, solv=None, rrho=None, t=None, out=False, cons f += getattr(self, e)["energy"] if solv is not None: if solv in ("cheap_prescreening_gsolv_info", "prescreening_gsolv_info"): - f += getattr(self, solv).get("range", {}).get(t,getattr(self, solv, {"energy": 0.0})["energy"]) + f += ( + getattr(self, solv) + .get("range", {}) + .get(t, getattr(self, solv, {"energy": 0.0})["energy"]) + ) else: f += getattr(self, solv)["range"].get(t, 0.0) if rrho is not None: f += self.get_mrrho( - t, - rrho=rrho, - consider_sym=consider_sym, - symnum=self.symnum, - ) + t, rrho=rrho, consider_sym=consider_sym, symnum=self.symnum + ) if not out: self.free_energy = f else: @@ -598,13 +612,31 @@ def calc_free_energy(self, e=None, solv=None, rrho=None, t=None, out=False, cons else: return f + def _get_sym_num(self, sym=None, linear=None): + """Get rotational symmetry number from Schoenfließ symbol""" + if sym is None: + sym = "c1" + if linear is None: + linear = False + symnum = 1 + if linear and "c" in sym.lower()[0]: + symnum = 1 + return symnum + elif linear and "d" in sym.lower()[0]: + symnum = 2 + return symnum + for key in rot_sym_num.keys(): + if key in sym.lower(): + symnum = rot_sym_num.get(key, 1) + break + return symnum + def calc_entropy_sym(self, temperature, symnum=None): """ RTln(sigma) rotational entropy""" if symnum is None: - symnum=self.symnum + symnum = self.symnum return R / AU2KCAL * temperature * log(symnum) - def get_mrrho( self, temperature, rrho=None, consider_sym=None, symnum=None, direct_input=0.0 ): @@ -620,10 +652,7 @@ def get_mrrho( .get(temperature, getattr(self, rrho, {"energy": 0.0})["energy"]) ) elif rrho in ("rrho_optimization",): - f += ( - getattr(self, "optimization_info") - .get("energy_rrho", 0.0) - ) + f += getattr(self, "optimization_info").get("energy_rrho", 0.0) else: f += getattr(self, rrho)["range"].get(temperature, 0.0) if consider_sym is None: diff --git a/censo_qm/ensembledata.py b/censo_qm/ensembledata.py index 9c8850f..575520c 100644 --- a/censo_qm/ensembledata.py +++ b/censo_qm/ensembledata.py @@ -35,6 +35,62 @@ def __init__( "part2": None, "part3": None, }, + supporting_info={ + "part0": { + "Energy": None, + "Energy_settings": None, + "G_mRRHO": None, + "G_solv": None, + "Geometry": None, + "Threshold": None, + "main QM code": None, + }, + "part1": { + "Energy": None, + "Energy_settings": None, + "G_mRRHO": None, + "G_solv": None, + "Geometry": None, + "Threshold": None, + "main QM code": None, + }, + "part2": { + "Energy": None, + "Energy_settings": None, + "G_mRRHO": None, + "G_solv": None, + "Geometry": None, + "Threshold": None, + "main QM code": None, + }, + "part3": { + "Energy": None, + "Energy_settings": None, + "G_mRRHO": None, + "G_solv": None, + "Geometry": None, + "Threshold": None, + "main QM code": None, + }, + "part4": { + "Energy": None, + "Energy_settings": None, + "G_mRRHO": None, + "G_solv": None, + "Geometry": None, + "Threshold": None, + "main QM code": None, + }, + "part5": { + "Energy": None, + "Energy_settings": None, + "G_mRRHO": None, + "G_solv": None, + "Geometry": None, + "Threshold": None, + "main QM code": None, + }, + }, ): """ ensemble_data: Creates an object where data @@ -59,3 +115,4 @@ def __init__( self.comment = comment self.bestconf = bestconf self.nconfs_per_part = nconfs_per_part + self.si = supporting_info diff --git a/censo_qm/inputhandling.py b/censo_qm/inputhandling.py index eeee45a..6b9a0fc 100755 --- a/censo_qm/inputhandling.py +++ b/censo_qm/inputhandling.py @@ -23,6 +23,8 @@ external_paths, __version__, cosmors_param, + composite_method_basis, + disp_already_included_in_func, ) from .utilities import frange, format_line, print @@ -444,7 +446,7 @@ def cml(startup_description, options, argv=None): ) group4.add_argument( "-ancopt", - choices=["on"], # there is no other option right now! + choices=["on"], # there is no other option right now! dest="ancopt", required=False, metavar="", @@ -727,7 +729,7 @@ def cml(startup_description, options, argv=None): required=False, metavar="", help="Investigates hydrogen nuclei in coupling and shielding calculations." - "choices=['on', 'off']", + "choices=['on', 'off']", ) group6.add_argument( "-cactive", @@ -737,7 +739,7 @@ def cml(startup_description, options, argv=None): required=False, metavar="", help="Investigates carbon nuclei in coupling and shielding calculations." - "choices=['on', 'off']", + "choices=['on', 'off']", ) group6.add_argument( "-factive", @@ -747,7 +749,7 @@ def cml(startup_description, options, argv=None): required=False, metavar="", help="Investigates fluorine nuclei in coupling and shielding calculations." - "choices=['on', 'off']", + "choices=['on', 'off']", ) group6.add_argument( "-siactive", @@ -757,7 +759,7 @@ def cml(startup_description, options, argv=None): required=False, metavar="", help="Investigates silicon nuclei in coupling and shielding calculations." - "choices=['on', 'off']", + "choices=['on', 'off']", ) group6.add_argument( "-pactive", @@ -767,7 +769,7 @@ def cml(startup_description, options, argv=None): required=False, metavar="", help="Investigates phosophorus nuclei in coupling and shielding calculations." - "choices=['on', 'off']", + "choices=['on', 'off']", ) group9 = parser.add_argument_group("OPTICAL ROTATION MODE") group9.add_argument( @@ -854,7 +856,7 @@ def cml(startup_description, options, argv=None): action="store", metavar="", help="Automatically balance the number of threads and cores when a low number" - "of conformers is left. (never exceed O*P cores).", + "of conformers is left. (never exceed O*P cores).", ) group11 = parser.add_argument_group("Concerning overall mRRHO calculations") group11.add_argument( @@ -916,6 +918,14 @@ def cml(startup_description, options, argv=None): help="Delete all unneeded files from current working directory. " "Stronger than -cleanup !", ) + group8.add_argument( + "--readonly", + "-readonly", + dest="onlyread", + choices=["on", "off"], + action="store", + help="Create new enso.json from the output of previous calculations.", + ) group8.add_argument( "-newconfig", "-write_censorc", @@ -945,7 +955,7 @@ def cml(startup_description, options, argv=None): choices=["on", "off"], default="off", help="Print progress to stderr when starting and finishing a sorting Part." - "Choices are 'on' or 'off'.", + "Choices are 'on' or 'off'.", ) group8.add_argument( "-inprc", @@ -962,6 +972,14 @@ def cml(startup_description, options, argv=None): action="store_true", help="Start interactive CENSO documentation.", ) + group8.add_argument( + "-create_SI", + "--create_SI", + dest="create_si", + required=False, + action="store_true", + help="Start automatic SI generation after CENSO run.", + ) args = parser.parse_args(argv) # apply logK settings but don't override user input! @@ -1088,6 +1106,7 @@ class internal_settings: "cosmorsparam": "cosmorsparam", "progress": "progress", "balance": "balance", + "onlyread": "onlyread", } knownbasissets3 = [ "SVP", @@ -1122,28 +1141,28 @@ class internal_settings: ] # information on functionals: - composite_method_basis = { - "pbeh-3c": "def2-mSVP", - "pbeh3c": "def2-mSVP", - "b97-3c": "def2-mTZVP", - "b973c": "def2-mTZVP", - "hf3c": "minix", - "hf-3c": "minix", - "r2scan-3c": "def2-mTZVPP", - } - gga_dfa = ("tpss", "pbe", "kt2") - hybrid_dfa = ( - "pbe0", - "pw6b95", - "wb97x-d3", - "cam-b3lyp", - "b3-lyp", - "pbeh-3c", - "m06x", - "bh-lyp", - "tpssh", - ) - dh_dfa = ("dsd-blyp",) + # composite_method_basis = { + # "pbeh-3c": "def2-mSVP", + # "pbeh3c": "def2-mSVP", + # "b97-3c": "def2-mTZVP", + # "b973c": "def2-mTZVP", + # "hf3c": "minix", + # "hf-3c": "minix", + # "r2scan-3c": "def2-mTZVPP", + # } + # gga_dfa = ("tpss", "pbe", "kt2") + # hybrid_dfa = ( + # "pbe0", + # "pw6b95", + # "wb97x-d3", + # "cam-b3lyp", + # "b3-lyp", + # "pbeh-3c", + # "m06x", + # "bh-lyp", + # "tpssh", + # ) + # dh_dfa = ("dsd-blyp",) knownbasissetsJ = knownbasissets3 + ["pcJ-0", "pcJ-1", "pcJ-2"] knownbasissetsS = knownbasissets3 + [ @@ -1156,12 +1175,38 @@ class internal_settings: ] func_orca = ["pbeh-3c", "b97-3c", "tpss", "b97-d3", "pbe"] func_tm = ["pbeh-3c", "b97-3c", "tpss", "r2scan-3c", "b97-d", "pbe"] - func3_orca = ["pw6b95", "pbe0","b97-d3", "wb97x-d3", "wb97x-d3bj", "wb97x-v", "dsd-blyp", "b97-3c", "pbeh-3c", "tpss", "pbe"] - func3_tm = ["pw6b95", "pbe0", "b97-d", "r2scan-3c", "b97-3c", "wb97x-v", "pbeh-3c", "tpss", "pbe",] - func_j_tm = ["tpss", "pbe0", "pbeh-3c", "r2scan-3c"] - func_j_orca = ["tpss", "pbe0", "pbeh-3c"] - func_s_tm = ["tpss", "pbe0", "pbeh-3c", "b97-3c", "kt1", "kt2", "r2scan-3c"] - func_s_orca = ["tpss", "pbe0", "dsd-blyp", "pbeh-3c", "kt2"] + func3_orca = [ + "pw6b95", + "pbe0", + "b97-d3", + "wb97x-d3", + "wb97x-d3bj", + "wb97x-v", + "dsd-blyp", + "b97-3c", + "pbeh-3c", + "tpss", + "pbe", + "b3lyp", + ] + func3_tm = [ + "pw6b95", + "pbe0", + "b97-d", + "r2scan-3c", + "b97-3c", + "wb97x-v", + "pbeh-3c", + "tpss", + "pbe", + "b3-lyp", + ] + func_j_tm = ["tpss", "pbe0", "pbeh-3c", "r2scan-3c", "pbe"] + func_j_orca = ["tpss", "pbe0", "pbeh-3c", "pbe"] + func_s_tm = ["tpss", "pbe0", "pbeh-3c", "b97-3c", "kt1", "kt2", "r2scan-3c", "pbe"] + func_s_orca = ["tpss", "pbe0", "dsd-blyp", "pbeh-3c", "kt2", "pbe"] + func_or_scf_tm = list(set(func_tm + func3_tm)) + func_or_tm = list(set(func_tm + func3_tm)) impgfnv = ["gfn1", "gfn2", "gfnff"] tmp_smd_solvents = [ "1,1,1-TRICHLOROETHANE", @@ -1529,7 +1574,7 @@ def __init__(self): ("omp", {"default": 1, "type": int}), ("balance", {"default": False, "type": bool}), ("cosmorsparam", {"default": "automatic", "type": str}), - + ("onlyread", {"default": False, "type": bool}), ] self.defaults_refine_ensemble_part0 = [ # part0 @@ -1650,7 +1695,7 @@ def __init__(self): "part4": ["on", "off"], "optical_rotation": ["on", "off"], "prog3": ["tm", "orca", "prog"], - "ancopt": ["on",], #, "off"], + "ancopt": ["on"], # , "off"], "opt_spearman": ["on", "off"], "evaluate_rrho": ["on", "off"], "consider_sym": ["on", "off"], @@ -1734,6 +1779,7 @@ def __init__(self): "sthr": ["automatic or e.g., 50 # in cm-1"], "scale": ["automatic or e.g., 1.0 "], "cosmorsparam": ["automatic"] + sorted(list(cosmors_param.keys())), + "onlyread": ["on", "off"], } # must not be changed if restart(concerning optimization) self.restart_unchangeable = [ @@ -1768,15 +1814,15 @@ def __init__(self): "basis_or": False, "func_or_scf": False, "freq_or": False, - "func3":False, - "basis3":False, - "func_j":False, - "basis_j":False, - "sm4_j":False, - "func_s":False, - "basis_s":False, - "sm4_s":False, - # "consider_sym": calculated on the fly + "func3": False, + "basis3": False, + "func_j": False, + "basis_j": False, + "sm4_j": False, + "func_s": False, + "basis_s": False, + "sm4_s": False, + # "consider_sym": calculated on the fly } @@ -1811,7 +1857,7 @@ def __init__(self, path=os.getcwd(), *args, **kwargs): self.basis = "automatic" self.maxthreads = 1 self.omp = 1 - self.balance=False + self.balance = False self.cosmorsparam = "automatic" # part0 self.part0 = False @@ -1885,7 +1931,7 @@ def __init__(self, path=os.getcwd(), *args, **kwargs): self.save_errors = [] self.save_infos = [] - self.startupinfokeys = ["nat", "md5", "maxconf", "run"] + self.startupinfokeys = ["nat", "md5", "maxconf", "run", "configpath"] self.nat = 0 self.md5 = "" self.maxconf = 0 @@ -1921,6 +1967,16 @@ def cleanup_run(self, complete=False): if int(file.split(".")[2]) > 1: print(f"Removing: {file}") os.remove(os.path.join(self.cwd, file)) + # get files like amat.tmp from COSMO calculation (can be several Mb) + files_in_sub_cwd = glob.glob(self.cwd + "/**/**/**/*mat.tmp", recursive=True) + size = 0 + for tmpfile in files_in_sub_cwd: + if any(x in tmpfile for x in ["a3mat.tmp", "a2mat.tmp", "amat.tmp"]): + if os.path.isfile(tmpfile): + size += os.path.getsize(tmpfile) + # print(f"Removing {tmpfile} {os.path.getsize(tmpfile)} byte") + os.remove(tmpfile) + print(f"Removed {size/(1024*1024): .2f} Mb") if complete: if "enso.json" in files_in_cwd: print(f"Removing: {'enso.json'}") @@ -1931,16 +1987,6 @@ def cleanup_run(self, complete=False): if os.path.isdir(os.path.join(self.cwd, "conformer_rotamer_check")): print("Removing conformer_rotamer_check") shutil.rmtree(os.path.join(self.cwd, "conformer_rotamer_check")) - # get files like amat.tmp from COSMO calculation (can be several Mb) - files_in_sub_cwd = glob.glob(self.cwd + '/**/**/**/*mat.tmp', recursive=True) - size = 0 - for tmpfile in files_in_sub_cwd: - if any(x in tmpfile for x in ['a3mat.tmp', 'a2mat.tmp', 'amat.tmp']): - if os.path.isfile(tmpfile): - size += os.path.getsize(tmpfile) - #print(f"Removing {tmpfile} {os.path.getsize(tmpfile)} byte") - os.remove(tmpfile) - print(f"Removed {size/1024} Mb") # ask if CONF folders should be removed def get_method_name( @@ -1962,17 +2008,14 @@ def get_method_name( --> method2 gsolv """ if func is not None and basis is not None: - if func in self.composite_method_basis.keys(): + if func in composite_method_basis.keys(): if basis == self.func_basis_default.get(func, None): # composite method (e.g. r2scan-3c) tmp_func_basis = func - # elif disp is not None: - # # FUNC-DISP/BASIS - # tmp_func_basis = f"{func}-{disp}/{basis}" else: # FUNC/BASIS tmp_func_basis = f"{func}/{basis}" - elif disp is not None: + elif disp is not None: #and func not in disp_already_included_in_func: # FUNC-DISP/BASIS tmp_func_basis = f"{func}-{disp}/{basis}" else: @@ -2259,9 +2302,13 @@ def check_logic(self, error_logical=False, silent=False): # Handle basis0 for func0: if self.basis0 == "None" or self.basis0 is None or self.basis0 == "automatic": if self.prog == "tm": - default = self.internal_defaults_tm.get("basis0", {'default':"def2-SV(P)"})['default'] + default = self.internal_defaults_tm.get( + "basis0", {"default": "def2-SV(P)"} + )["default"] elif self.prog == "orca": - default = self.internal_defaults_orca.get("basis0",{'default':"def2-SV(P)"})['default'] + default = self.internal_defaults_orca.get( + "basis0", {"default": "def2-SV(P)"} + )["default"] else: default = "def2-SV(P)" self.basis0 = self.func_basis_default.get(self.func0, default) @@ -2284,9 +2331,13 @@ def check_logic(self, error_logical=False, silent=False): # Handle basis for func: if self.basis == "None" or self.basis is None or self.basis == "automatic": if self.prog == "tm": - default = self.internal_defaults_tm.get("basis", {'default':"def2-TZVP"})['default'] + default = self.internal_defaults_tm.get( + "basis", {"default": "def2-TZVP"} + )["default"] elif self.prog == "orca": - default = self.internal_defaults_orca.get("basis", {'default':"def2-TZVP"})['default'] + default = self.internal_defaults_orca.get( + "basis", {"default": "def2-TZVP"} + )["default"] else: default = "def2-TZVP" self.basis = self.func_basis_default.get(self.func, default) @@ -2294,29 +2345,33 @@ def check_logic(self, error_logical=False, silent=False): # Handle basis3: if self.basis3 == "None" or self.basis3 is None or self.basis3 == "automatic": if self.prog3 == "tm": - default = self.internal_defaults_tm.get("basis3", {'default':"def2-TZVP"})['default'] + default = self.internal_defaults_tm.get( + "basis3", {"default": "def2-TZVP"} + )["default"] elif self.prog3 == "orca": - default = self.internal_defaults_orca.get("basis3", {'default':"def2-TZVP"})['default'] + default = self.internal_defaults_orca.get( + "basis3", {"default": "def2-TZVP"} + )["default"] else: default = "def2-TZVP" self.basis3 = self.func_basis_default.get(self.func3, default) # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Handle func3 - if self.part3 and self.func3 in self.composite_method_basis.keys(): - if self.basis3 != self.composite_method_basis[self.func3]: + if self.part3 and self.func3 in composite_method_basis.keys(): + if self.basis3 != composite_method_basis[self.func3]: self.save_errors.append( f"{'INFORMATION:':{WARNLEN}}You are using a basis " f"set (basis3) {self.basis3} different to the original composite method" - f" basis set ({self.composite_method_basis[self.func3]})!" + f" basis set ({composite_method_basis[self.func3]})!" ) - if self.part3 and self.prog3 == 'orca' and self.func3 not in self.func3_orca: + if self.part3 and self.prog3 == "orca" and self.func3 not in self.func3_orca: self.save_errors.append( f"{'ERROR:':{WARNLEN}}The functional " f"(func3) {self.func3} is not implemented with the {self.prog3} program package.\n" f"{'':{WARNLEN}}Options are: {self.func3_orca}" ) error_logical = True - elif self.part3 and self.prog3 == 'tm' and self.func3 not in self.func3_tm: + elif self.part3 and self.prog3 == "tm" and self.func3 not in self.func3_tm: self.save_errors.append( f"{'ERROR:':{WARNLEN}}The functional " f"(func3) {self.func3} is not implemented with the {self.prog3} program package.\n" @@ -2349,7 +2404,7 @@ def check_logic(self, error_logical=False, silent=False): getattr(self, "basis_j", None) is None or getattr(self, "basis_j", None) == "automatic" ): - default_basis_j = self.composite_method_basis.get(self.func_j, "def2-TZVP") + default_basis_j = composite_method_basis.get(self.func_j, "def2-TZVP") setattr(self, "basis_j", default_basis_j) # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Handle func_s @@ -2374,7 +2429,7 @@ def check_logic(self, error_logical=False, silent=False): getattr(self, "basis_s", None) is None or getattr(self, "basis_s", None) == "automatic" ): - default_basis_s = self.composite_method_basis.get(self.func_s, "def2-TZVP") + default_basis_s = composite_method_basis.get(self.func_s, "def2-TZVP") setattr(self, "basis_s", default_basis_s) # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # no unpaired electrons in coupling or shielding calculations! @@ -2649,7 +2704,9 @@ def check_logic(self, error_logical=False, silent=False): f" for solventmodel/program {key} can not be checked but is used anyway." ) else: - if censo_solvent_db[self.solvent].get(key, "nothing_found")[1] not in getattr(self, lookup[key]): + if censo_solvent_db[self.solvent].get(key, "nothing_found")[ + 1 + ] not in getattr(self, lookup[key]): self.save_errors.append( f"{'WARNING:':{WARNLEN}}The solvent " f"{censo_solvent_db[self.solvent].get(key, 'nothing_found')[1]} " @@ -2657,12 +2714,16 @@ def check_logic(self, error_logical=False, silent=False): ) # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # prog3 if smgsolv3 cosmors - if self.part3 and self.solvent != 'gas' and self.smgsolv3 in ('cosmors', 'cosmors-fine'): - if self.prog3 == 'orca': + if ( + self.part3 + and self.solvent != "gas" + and self.smgsolv3 in ("cosmors", "cosmors-fine") + ): + if self.prog3 == "orca": self.save_errors.append( - f"{'ERROR:':{WARNLEN}}Part3 when used with COSMO-RS can only be prog3 = TM!\n" - f"{'':{WARNLEN}}Set prog3 to tm !" - ) + f"{'ERROR:':{WARNLEN}}Part3 when used with COSMO-RS can only be prog3 = TM!\n" + f"{'':{WARNLEN}}Set prog3 to tm !" + ) error_logical = True # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Handle optlevel2: @@ -2703,16 +2764,16 @@ def check_logic(self, error_logical=False, silent=False): ) # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # optical rotation - if self.optical_rotation and self.prog == 'orca': + if self.optical_rotation and self.prog == "orca": self.save_errors.append( - f"{'ERROR:':{WARNLEN}}Optical rotation calculations are only possible with TM." - ) + f"{'ERROR:':{WARNLEN}}Optical rotation calculations are only possible with TM." + ) error_logical = True if silent: self.save_errors = store_errors return error_logical - def print_parameters(self): + def print_parameters(self, cmlcall=None): """ print settings at startup """ @@ -2722,6 +2783,8 @@ def print_parameters(self): print("PARAMETERS".center(PLENGTH, " ")) print("".ljust(PLENGTH, "-") + "\n") + if cmlcall is not None: + print(f"program call: {' '.join(cmlcall)}") print( f"The configuration file {os.path.basename(self.configpath)} is read " f"from {self.configpath}." @@ -2730,6 +2793,14 @@ def print_parameters(self): if self.save_infos: for _ in list(self.save_infos): print(self.save_infos.pop(0)) + if self.onlyread: + self.save_errors.append( + f"{'WARNING:':{WARNLEN}}Using the option ``readonly`` to re-read" + + f" data from old outputs. This option is experimental therefore " + + f"check results carefully.\n" + + f"{' ':{WARNLEN}}It can only succeed if exactly the same " + + f"input commands (.censorc + command line) are supplied!" + ) if self.save_errors: print("") for _ in list(self.save_errors): @@ -2826,7 +2897,12 @@ def print_parameters(self): if self.solvent != "gas": info.append(["sm_rrho", "solvent model applied with xTB"]) info.append(["printoption", "evaluate at different temperatures", "off"]) - info.append(["part1_threshold", "threshold g_thr(1) and G_thr(1) for sorting in part1"]) + info.append( + [ + "part1_threshold", + "threshold g_thr(1) and G_thr(1) for sorting in part1", + ] + ) if self.solvent != "gas": info.append( ["smgsolv1", "solvent model for Gsolv contribution of part1"] @@ -2881,7 +2957,7 @@ def print_parameters(self): info.append( [ "opt_limit", - "optimize all conformers below this G_thr(opt,2) threshold" + "optimize all conformers below this G_thr(opt,2) threshold", ] ) info.append(["printoption", "spearmanthr", f"{self.spearmanthr:.3f}"]) @@ -2893,7 +2969,10 @@ def print_parameters(self): info.append(["smgsolv2", "solvent model for Gsolv contribution"]) info.append(["multitemp", "evaluate at different temperatures"]) info.append( - ["part2_threshold", "Boltzmann sum threshold G_thr(2) for sorting in part2"] + [ + "part2_threshold", + "Boltzmann sum threshold G_thr(2) for sorting in part2", + ] ) info.append(["evaluate_rrho", "calculate mRRHO contribution"]) if self.evaluate_rrho: @@ -2906,8 +2985,8 @@ def print_parameters(self): info.append( [ "bhess", - "Apply constraint to input geometry "+ - "during mRRHO calculation", + "Apply constraint to input geometry " + + "during mRRHO calculation", ] ) if self.solvent != "gas": @@ -3054,8 +3133,25 @@ def print_parameters(self): if getattr(self, "p_active"): info.append(["p_active", "Calculating phosphorus spectrum"]) info.append(["p_ref", "reference for 31P"]) - if all([ not getattr(self, active) for active in ("h_active", "c_active", "f_active", "si_active", "p_active")]): - info.append(["printoption", "Calculating spectrum for all nuclei", "H, C, F, Si, P"]) + if all( + [ + not getattr(self, active) + for active in ( + "h_active", + "c_active", + "f_active", + "si_active", + "p_active", + ) + ] + ): + info.append( + [ + "printoption", + "Calculating spectrum for all nuclei", + "H, C, F, Si, P", + ] + ) info.append(["resonance_frequency", "resonance frequency"]) # short notation: @@ -3080,17 +3176,17 @@ def print_parameters(self): max_len_digilen = 0 for item in info: - if item[0] == 'justprint': + if item[0] == "justprint": if "short-notation" in item[1]: continue - #tmp = len(item[1]) -len('short-notation:') + # tmp = len(item[1]) -len('short-notation:') else: tmp = len(item[1]) else: tmp = len(item[1]) if tmp > max_len_digilen: max_len_digilen = tmp - max_len_digilen +=1 + max_len_digilen += 1 if max_len_digilen < DIGILEN: max_len_digilen = DIGILEN @@ -3126,6 +3222,134 @@ def print_parameters(self): ) print("END of parameters\n") + def create_SI(self, ensembledata): + """Automatically create a supporting information block""" + si_out = [] + l_length = 25 + si_out.append(f"\nAutomatic CENSO-SI creation:\n") + si_out.append( + "The preparation of this supporting information (SI) is by no means complete and " + + "does not relieve the\nuser of the responsibility to verify/complete the SI." + ) + si_out.append(f"In this CENSO run the following settings were employed:\n") + # General information: + tmp = {"Temperature": self.temperature, "Solvent": self.solvent} + si_out.append(f"{'':^{l_length}} | {'General information:'}") + si_out.append("".ljust(int(PLENGTH), "-")) + for key, value in tmp.items(): + si_out.append(f"{key:^{l_length}} | {value}") + si_out.append("".ljust(int(PLENGTH), "-")) + # ----------------------------PART0-------------------------------------- + if self.part0: + # part0: E = fast single-point with normally b97-d/def2-SV(P) + # GRRHO is omitted for speed! + # Gsolv = ALPB_GSOLV + si_out.append(f"{'':^{l_length}} | {'PART0 - cheap prescreening'}") + si_out.append("".ljust(int(PLENGTH), "-")) + for key, value in getattr(ensembledata, "si")["part0"].items(): + si_out.append(f"{key:^{l_length}} | {value}") + si_out.append("".ljust(int(PLENGTH), "-")) + # ==========================END PART0==================================== + # ----------------------------PART1-------------------------------------- + if self.part1: + # part1: E normally with r2scan-3c + # Gsolv can be anything + # GmRRHO xtb with GFNn-xTB + ALPB GBSA gas-phase + ### + si_out.append(f"{'':^{l_length}} | {'PART1 - prescreening'}") + si_out.append("".ljust(int(PLENGTH), "-")) + for key, value in getattr(ensembledata, "si")["part1"].items(): + si_out.append(f"{key:^{l_length}} | {value}") + si_out.append("".ljust(int(PLENGTH), "-")) + # ==========================END PART1==================================== + # ----------------------------PART2-------------------------------------- + if self.part2: + # part2: E normally with r2scan-3c + # Gsolv can be anything + # GmRRHO xtb with GFNn-xTB + ALPB GBSA gas-phase + ### + si_out.append(f"{'':^{l_length}} | {'PART2 - optimization'}") + si_out.append("".ljust(int(PLENGTH), "-")) + for key, value in getattr(ensembledata, "si")["part2"].items(): + si_out.append(f"{key:^{l_length}} | {value}") + si_out.append("".ljust(int(PLENGTH), "-")) + # ==========================END PART2==================================== + # ----------------------------PART3-------------------------------------- + if self.part3: + # # part3: E normally hybrid dft + # # Gsolv can be anything + # # GmRRHO xtb with GFNn-xTB + ALPB GBSA gas-phase + si_out.append(f"{'':^{l_length}} | {'PART3 - refinement'}") + si_out.append("".ljust(int(PLENGTH), "-")) + for key, value in getattr(ensembledata, "si")["part3"].items(): + si_out.append(f"{key:^{l_length}} | {value}") + si_out.append("".ljust(int(PLENGTH), "-")) + # ==========================END PART3==================================== + # ----------------------------PART4-------------------------------------- + if self.part4: + ### + si_out.append(f"{'':^{l_length}} | {'PART4 - NMR mode'}") + si_out.append("".ljust(int(PLENGTH), "-")) + for key, value in getattr(ensembledata, "si")["part4"].items(): + si_out.append(f"{key:^{l_length}} | {value}") + si_out.append("".ljust(int(PLENGTH), "-")) + # ==========================END PART4==================================== + # ----------------------------PART5-------------------------------------- + if self.optical_rotation: + ### + si_out.append(f"{'':^{l_length}} | {'PART5 - OR mode'}") + si_out.append("".ljust(int(PLENGTH), "-")) + for key, value in getattr(ensembledata, "si")["part5"].items(): + si_out.append(f"{key:^{l_length}} | {value}") + si_out.append("".ljust(int(PLENGTH), "-")) + # ==========================END PART5==================================== + # PRINT SI PART information + for line in si_out: + print(line) + si_text = """ + Within CENSO free energies (G) are calculated from: + + G(T) = E + G_mRRHO(T) + G_solv(T) (1) + + The input ensemble (originating from CREST or crest_combi) are sorted + in part0 and part1 at DFT level. This provides an improved description of + the electronic energy and solvation contributions compared to the input + SQM level. Thermostatistical contributions, including ZPVE are evaluated + at the GFNn-xTB level in part1. This first filtering removes structures + high lying in free energy and reduces the computational cost by passing + fewer conformers to the computational costly DFT optimizations. + + COSMO-RS calculations are performed with the energies and densities from + the functional and basis set combination of the 'Energy' evaluation in the + respective part. + + """ + print(si_text) + print("\nSome citations in bib style are provided below:") + # PRINT bib of employed programs + from .cfg import si_bib + + out_bib = [] + requirements = self.needed_external_programs() + bib_need = { + "needtm": "tm", + "needxtb": "xtb", + "needcosmors": "cosmors", + "needorca": "orca", + } + out_bib.extend(si_bib.get("censo")) + for item in ("needtm", "needxtb", "needcosmors", "needorca"): + if requirements.get(item, False): + out_bib.extend(si_bib.get(bib_need[item], [])) + for item in set([self.func0, self.func, self.func3, self.func_j, self.func_s]): + if si_bib.get(item, None) is not None: + out_bib.extend(si_bib.get(item, [])) + if self.evaluate_rrho and self.bhess: + out_bib.extend(si_bib.get("sph", [])) + print("\nBib entries:") + for line in out_bib: + print(line) + def read_program_paths(self, configpath): """ Get absolute paths of external programs employed in enso @@ -3164,7 +3388,9 @@ def read_program_paths(self, configpath): try: self.external_paths["orcapath"] = str(line.split()[1]) except Exception: - print(f"{'WARNING:':{WARNLEN}}Could not read path for ORCA from .censorc!.") + print( + f"{'WARNING:':{WARNLEN}}Could not read path for ORCA from .censorc!." + ) if "ORCA version:" in line: try: tmp = line.split()[2] @@ -3173,12 +3399,16 @@ def read_program_paths(self, configpath): tmp = "".join(tmp) self.external_paths["orcaversion"] = tmp except Exception: - print(f"{'WARNING:':{WARNLEN}}Could not read ORCA version from .censorc!") + print( + f"{'WARNING:':{WARNLEN}}Could not read ORCA version from .censorc!" + ) if "GFN-xTB:" in line: try: self.external_paths["xtbpath"] = str(line.split()[1]) except Exception: - print(f"{'WARNING:':{WARNLEN}}Could not read path for GFNn-xTB from .censorc!") + print( + f"{'WARNING:':{WARNLEN}}Could not read path for GFNn-xTB from .censorc!" + ) if shutil.which("xtb") is not None: self.external_paths["xtbpath"] = shutil.which("xtb") print( @@ -3188,7 +3418,9 @@ def read_program_paths(self, configpath): try: self.external_paths["crestpath"] = str(line.split()[1]) except Exception: - print(f"{'WARNING:':{WARNLEN}}Could not read path for CREST from .censorc!") + print( + f"{'WARNING:':{WARNLEN}}Could not read path for CREST from .censorc!" + ) if shutil.which("crest") is not None: self.external_paths["crestpath"] = shutil.which("crest") print( @@ -3198,16 +3430,20 @@ def read_program_paths(self, configpath): try: self.external_paths["mpshiftpath"] = str(line.split()[1]) except Exception: - print(f"{'WARNING:':{WARNLEN}}Could not read path for mpshift from .censorc!") + print( + f"{'WARNING:':{WARNLEN}}Could not read path for mpshift from .censorc!" + ) if "escf:" in line: try: self.external_paths["escfpath"] = str(line.split()[1]) except Exception: - print(f"{'WARNING:':{WARNLEN}}Could not read path for escf from .censorc!") + print( + f"{'WARNING:':{WARNLEN}}Could not read path for escf from .censorc!" + ) if "$ENDPROGRAMS" in line: break - def needed_external_programs(self, config): + def needed_external_programs(self): """ Automatically checks which external programs are required for the current run. @@ -3215,56 +3451,56 @@ def needed_external_programs(self, config): requirements = {} # xTB if ( - config.prog_rrho == "xtb" - or config.part0 - or config.ancopt - or config.smgsolv2 in ("gbsa_gsolv", "alpb_gsolv") + self.prog_rrho == "xtb" + or self.part0 + or self.ancopt + or self.smgsolv2 in ("gbsa_gsolv", "alpb_gsolv") ): requirements["needxtb"] = True # TM if ( - config.prog == "tm" - or config.prog3 == "tm" - or config.prog4_j == "tm" - or config.prog4_s == "tm" - or config.smgsolv1 in ("cosmors", "cosmors-fine") - or config.smgsolv2 in ("cosmors", "cosmors-fine") - or config.smgsolv3 in ("cosmors", "cosmors-fine") + self.prog == "tm" + or self.prog3 == "tm" + or self.prog4_j == "tm" + or self.prog4_s == "tm" + or self.smgsolv1 in ("cosmors", "cosmors-fine") + or self.smgsolv2 in ("cosmors", "cosmors-fine") + or self.smgsolv3 in ("cosmors", "cosmors-fine") ): requirements["needtm"] = True requirements["needcefine"] = True - if config.part4 and (config.prog4_j == "tm" or config.prog4_s == "tm"): - if config.couplings: + if self.part4 and (self.prog4_j == "tm" or self.prog4_s == "tm"): + if self.couplings: requirements["needescf"] = True - if config.shieldings: + if self.shieldings: requirements["needmpshift"] = True # COSMORS if ( - (config.part1 and config.smgsolv1 == "cosmors") - or (config.part2 and config.smgsolv2 == "cosmors") - or (config.part3 and config.smgsolv3 == "cosmors") + (self.part1 and self.smgsolv1 == "cosmors") + or (self.part2 and self.smgsolv2 == "cosmors") + or (self.part3 and self.smgsolv3 == "cosmors") ): requirements["needcosmors"] = True requirements["needcosmors-normal"] = True if ( - (config.part1 and config.smgsolv1 == "cosmors-fine") - or (config.part2 and config.smgsolv2 == "cosmors-fine") - or (config.part3 and config.smgsolv3 == "cosmors-fine") + (self.part1 and self.smgsolv1 == "cosmors-fine") + or (self.part2 and self.smgsolv2 == "cosmors-fine") + or (self.part3 and self.smgsolv3 == "cosmors-fine") ): requirements["needcosmors"] = True requirements["needcosmors-fine"] = True # ORCA if ( - config.prog == "orca" - or config.prog3 == "orca" - or config.prog4_j == "orca" - or config.prog4_s == "orca" - or config.smgsolv1 == "smd_gsolv" - or config.smgsolv2 == "smd_gsolv" - or config.smgsolv3 == "smd_gsolv" + self.prog == "orca" + or self.prog3 == "orca" + or self.prog4_j == "orca" + or self.prog4_s == "orca" + or self.smgsolv1 == "smd_gsolv" + or self.smgsolv2 == "smd_gsolv" + or self.smgsolv3 == "smd_gsolv" ): requirements["needorca"] = True - if config.run: + if self.run: requirements["startenso"] = True return requirements @@ -3448,10 +3684,14 @@ def processQMpaths(self, requirements, error_logical): # COSMORS if requirements.get("needcosmors", False): if self.external_paths["cosmorssetup"] is None: - print(f"{'ERROR:':{WARNLEN}}Set up for COSMO-RS has to be written to .censorc!") + print( + f"{'ERROR:':{WARNLEN}}Set up for COSMO-RS has to be written to .censorc!" + ) error_logical = True if self.external_paths["cosmothermversion"] is None: - print(f"{'ERROR:':{WARNLEN}}Version of COSMO-RS has to be written to .censorc!") + print( + f"{'ERROR:':{WARNLEN}}Version of COSMO-RS has to be written to .censorc!" + ) error_logical = True if shutil.which("cosmotherm") is not None: print(" Using COSMOtherm from {}".format(shutil.which("cosmotherm"))) @@ -3601,7 +3841,10 @@ def write_censo_inp(self, path=None): value = self._exchange_onoff( data.get( key, - OrderedDict(self.defaults_refine_ensemble_general)[key]["default"]), + OrderedDict(self.defaults_refine_ensemble_general)[key][ + "default" + ], + ), reverse=True, ) if key == "nconf" and value is None: @@ -3613,7 +3856,10 @@ def write_censo_inp(self, path=None): value = self._exchange_onoff( data.get( key, - OrderedDict(self.defaults_refine_ensemble_part0)[key]["default"]), + OrderedDict(self.defaults_refine_ensemble_part0)[key][ + "default" + ], + ), reverse=True, ) key = args_key.get(key, key) @@ -3621,9 +3867,13 @@ def write_censo_inp(self, path=None): outdata.write("\n$PART1 - PRESCREENING - SETTINGS:\n") outdata.write("# func and basis is set under GENERAL SETTINGS\n") for key in OrderedDict(self.defaults_refine_ensemble_part1): - value = self._exchange_onoff(data.get( - key, - OrderedDict(self.defaults_refine_ensemble_part1)[key]["default"]), + value = self._exchange_onoff( + data.get( + key, + OrderedDict(self.defaults_refine_ensemble_part1)[key][ + "default" + ], + ), reverse=True, ) key = args_key.get(key, key) @@ -3631,18 +3881,26 @@ def write_censo_inp(self, path=None): outdata.write("\n$PART2 - OPTIMIZATION - SETTINGS:\n") outdata.write("# func and basis is set under GENERAL SETTINGS\n") for key in OrderedDict(self.defaults_refine_ensemble_part2): - value = self._exchange_onoff(data.get( - key, - OrderedDict(self.defaults_refine_ensemble_part2)[key]["default"]), + value = self._exchange_onoff( + data.get( + key, + OrderedDict(self.defaults_refine_ensemble_part2)[key][ + "default" + ], + ), reverse=True, ) key = args_key.get(key, key) outdata.write(format_line(key, value, "")) outdata.write("\n$PART3 - REFINEMENT - SETTINGS:\n") for key in OrderedDict(self.defaults_refine_ensemble_part3): - value = self._exchange_onoff(data.get( - key, - OrderedDict(self.defaults_refine_ensemble_part3)[key]["default"]), + value = self._exchange_onoff( + data.get( + key, + OrderedDict(self.defaults_refine_ensemble_part3)[key][ + "default" + ], + ), reverse=True, ) key = args_key.get(key, key) @@ -3650,9 +3908,10 @@ def write_censo_inp(self, path=None): outdata.write("\n$NMR PROPERTY SETTINGS:\n") outdata.write("$PART4 SETTINGS:\n") for key in OrderedDict(self.defaults_nmrprop_part4): - value = self._exchange_onoff(data.get( - key, - OrderedDict(self.defaults_nmrprop_part4)[key]["default"]), + value = self._exchange_onoff( + data.get( + key, OrderedDict(self.defaults_nmrprop_part4)[key]["default"] + ), reverse=True, ) key = args_key.get(key, key) @@ -3660,15 +3919,20 @@ def write_censo_inp(self, path=None): outdata.write("\n$OPTICAL ROTATION PROPERTY SETTINGS:\n") outdata.write("$PART5 SETTINGS:\n") for key in OrderedDict(self.defaults_optical_rotation_part5): - value = self._exchange_onoff(data.get( - key, - OrderedDict(self.defaults_optical_rotation_part5)[key]["default"]), + value = self._exchange_onoff( + data.get( + key, + OrderedDict(self.defaults_optical_rotation_part5)[key][ + "default" + ], + ), reverse=True, ) key = args_key.get(key, key) outdata.write(format_line(key, value, "")) outdata.write("\n$END censo.inp\n") -########################## + + ########################## def read_json(self, path, silent=False): """ diff --git a/censo_qm/nmrproperties.py b/censo_qm/nmrproperties.py index 61558fd..37a9e8c 100644 --- a/censo_qm/nmrproperties.py +++ b/censo_qm/nmrproperties.py @@ -8,14 +8,7 @@ import json from random import normalvariate from multiprocessing import JoinableQueue as Queue -from .cfg import ( - PLENGTH, - DIGILEN, - AU2KCAL, - WARNLEN, - CODING, - NmrRef -) +from .cfg import PLENGTH, DIGILEN, AU2KCAL, WARNLEN, CODING, NmrRef, qm_prepinfo from .parallel import run_in_parallel from .orca_job import OrcaJob from .tm_job import TmJob @@ -114,16 +107,16 @@ def average_shieldings(config, calculate, element_ref_shield, energy, solv, rrho for _ in range(1000): for conf in calculate: conf.calc_free_energy( - e=energy, - solv=solv, + e=energy, + solv=solv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) - conf.free_energy += normalvariate( - 0.0, get_std_dev(conf) + consider_sym=config.consider_sym, ) - calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) + conf.free_energy += normalvariate(0.0, get_std_dev(conf)) + calculate = calc_boltzmannweights( + calculate, "free_energy", config.temperature + ) tmp_sigma = {} for i in range(1, config.nat + 1): tmp_sigma[i] = 0 @@ -145,12 +138,12 @@ def average_shieldings(config, calculate, element_ref_shield, energy, solv, rrho solv=solv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) - conf.free_energy += normalvariate( - 0.0, (0.4/AU2KCAL) + consider_sym=config.consider_sym, ) - calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) + conf.free_energy += normalvariate(0.0, (0.4 / AU2KCAL)) + calculate = calc_boltzmannweights( + calculate, "free_energy", config.temperature + ) tmp_sigma = {} for i in range(1, config.nat + 1): tmp_sigma[i] = 0 @@ -163,14 +156,16 @@ def average_shieldings(config, calculate, element_ref_shield, energy, solv, rrho for atom in conf.shieldings.keys(): sigma_std_dev_const[atom].append(tmp_sigma[atom]) - with open(os.path.join(config.cwd, 'averaged_shift.dat'), "w", newline=None) as out: - line=("# in coord element σ(sigma) SD(σ based on SD Gsolv) SD(σ by 0.4 kcal/mol) shift σ_ref") + with open(os.path.join(config.cwd, "averaged_shift.dat"), "w", newline=None) as out: + line = "# in coord element σ(sigma) SD(σ based on SD Gsolv) SD(σ by 0.4 kcal/mol) shift σ_ref" print(line) - out.write(line+"\n") - line= ("".ljust(int(105), "-")) + out.write(line + "\n") + line = "".ljust(int(105), "-") print(line) - out.write(line+"\n") - maxsigma = max([len(str(sigma).split(".")[0]) for sigma in averaged.values()]) + 5 + out.write(line + "\n") + maxsigma = ( + max([len(str(sigma).split(".")[0]) for sigma in averaged.values()]) + 5 + ) make_shift = ( lambda atom: f"{-sigma+element_ref_shield.get(element[atom], 0.0):> {maxsigma}.2f}" if (element_ref_shield.get(element[atom], None) is not None) @@ -191,14 +186,53 @@ def average_shieldings(config, calculate, element_ref_shield, energy, solv, rrho f"{std_dev:^ 24.6f} {std_dev_const:^ 24.6f} {make_shift(atom):>5} {float(element_ref_shield.get(element[atom], 0.0)):> 10.3f}" ) print(line) - out.write(line+"\n") + out.write(line + "\n") except Exception: - line=(f"{atom:< {10}} {element[atom]:^{7}} {sigma:> {maxsigma}.2f}") + line = f"{atom:< {10}} {element[atom]:^{7}} {sigma:> {maxsigma}.2f}" print(line) - out.write(line+"\n") - line = ("".ljust(int(105), "-")) + out.write(line + "\n") + line = "".ljust(int(105), "-") print(line) - out.write(line+"\n") + out.write(line + "\n") + + ############################################################################ + # print min and max shielding within all conformers with averaging, but no boltzmann weighting + atom_sigma = {} + for atom in conf.shieldings.keys(): + atom_sigma[atom] = {} + for conf in calculate: + # get shielding constants + atom_sigma[atom][conf.id] = sum( + [conf.shieldings.get(eq_atom, 0.0) for eq_atom in chemeq[atom]] + ) / len(chemeq[atom]) + tmppl = maxsigma + 1 + if tmppl <= 8: + tmppl = 8 + line = f"\n{'# in coord':<{10}} {'element':^{7}} {'σ(sigma)':>{tmppl}} {'min(σ)*':>{tmppl}} CONFX {'max(σ)*':>{tmppl}} CONFX Δ(max-min)" + print(line) + line = "".ljust(int(105), "-") + print(line) + for atom in conf.shieldings.keys(): + try: + minx = min(list(atom_sigma[atom].values())) + minid = next( + key for key, value in atom_sigma[atom].items() if value == minx + ) + maxx = max(list(atom_sigma[atom].values())) + maxid = next( + key for key, value in atom_sigma[atom].items() if value == maxx + ) + line = (f"{atom:< {10}} {element[atom]:^{7}} {averaged[atom]:> {maxsigma}.2f} " + + f"{minx :> {maxsigma}.2f} {'CONF'+str(minid)} {maxx :> {maxsigma}.2f} " + + f"{'CONF'+str(maxid)} {maxx-minx:> {maxsigma}.2f}") + print(line) + except Exception: + pass + line = "".ljust(int(105), "-") + print(line) + print( + f"* min(σ) and max(σ) are averaged over the chemical equivalent atoms, but not Boltzmann weighted." + ) def write_anmrrc(config): @@ -379,7 +413,6 @@ def write_anmrrc(config): return refs - def part4(config, conformers, store_confs, ensembledata): """ Calculate nmr properties: shielding and coupling constants on the populated @@ -425,23 +458,30 @@ def part4(config, conformers, store_confs, ensembledata): if getattr(config, "p_active"): info.append(["p_active", "Calculating phosphorus spectrum"]) info.append(["p_ref", "reference for 31P"]) - if all([ not getattr(config, active) for active in ("h_active", "c_active", "f_active", "si_active", "p_active")]): - info.append(["printoption", "Calculating spectrum for all nuclei", "H, C, F, Si, P"]) + if all( + [ + not getattr(config, active) + for active in ("h_active", "c_active", "f_active", "si_active", "p_active") + ] + ): + info.append( + ["printoption", "Calculating spectrum for all nuclei", "H, C, F, Si, P"] + ) info.append(["resonance_frequency", "spectrometer frequency"]) # active nuclei max_len_digilen = 0 for item in info: - if item[0] == 'justprint': + if item[0] == "justprint": if "short-notation" in item[1]: - tmp = len(item[1]) -len('short-notation:') + tmp = len(item[1]) - len("short-notation:") else: tmp = len(item[1]) else: tmp = len(item[1]) if tmp > max_len_digilen: max_len_digilen = tmp - max_len_digilen +=1 + max_len_digilen += 1 if max_len_digilen < DIGILEN: max_len_digilen = DIGILEN @@ -514,7 +554,7 @@ def part4(config, conformers, store_confs, ensembledata): boltzmannthr = config.part2_threshold elif config.part1: # part2 is not calculated use Boltzmann weights directly from part1 - #--> misappropriate config.part2_threshold + # --> misappropriate config.part2_threshold # This means starting from not DFT optimized geometries! energy = "prescreening_sp_info" rrho = "prescreening_grrho_info" @@ -530,19 +570,31 @@ def part4(config, conformers, store_confs, ensembledata): elif getattr(conf, rrho)["info"] != "calculated" and config.evaluate_rrho: store_confs.append(mol) continue - elif ( - getattr(conf, gsolv)["info"] != "calculated" and config.solvent != "gas" - ): + elif getattr(conf, gsolv)["info"] != "calculated" and config.solvent != "gas": store_confs.append(mol) continue else: calculate.append(mol) if unoptimized_warning: - print_errors(f"{'INFORMATION:':{WARNLEN}} Conformers have not been optimized at DFT level!!!\n" - f"{'':{WARNLEN}}Use results with care!\n", save_errors + print_errors( + f"{'INFORMATION:':{WARNLEN}} Conformers have not been optimized at DFT level!!!\n" + f"{'':{WARNLEN}}Use results with care!\n", + save_errors, ) - + # SI geometry: + if unoptimized_warning: + ensembledata.si["part4"]["Geometry"] = "GFNn-xTB (input geometry)" + else: + ensembledata.si["part4"]["Geometry"], _ = config.get_method_name( + "xtbopt", + func=config.func, + basis=config.basis, + solvent=config.solvent, + sm=config.sm2, + ) + ensembledata.si["part4"]["Geometry"] += f" @optlevel: {config.optlevel2}" + # if not calculate and not prev_calculated: print(f"{'ERROR:':{WARNLEN}}No conformers left!") print("Going to exit!") @@ -554,6 +606,7 @@ def part4(config, conformers, store_confs, ensembledata): # Calculate Boltzmann weight for confs: if config.part3: + using_part = "part3 - refinement" if not config.evaluate_rrho: rrho = None rrho_method = None @@ -591,6 +644,7 @@ def part4(config, conformers, store_confs, ensembledata): solvent=config.solvent, ) elif config.part2: + using_part = "part2 - optimization" if not config.evaluate_rrho: rrho = None rrho_method = None @@ -628,6 +682,7 @@ def part4(config, conformers, store_confs, ensembledata): solvent=config.solvent, ) elif config.part1: + using_part = "part1 - prescreening" # on DFT unoptimized geometries! if not config.evaluate_rrho: rrho = None @@ -665,7 +720,19 @@ def part4(config, conformers, store_confs, ensembledata): gfn_version=config.part1_gfnv, solvent=config.solvent, ) - + # SI information:------------------------------------------------------- + ensembledata.si["part4"]["Energy"] = f"using Energy from {using_part}" + ensembledata.si["part4"]["Energy_settings"] = f"see {using_part}" + ensembledata.si["part4"]["G_mRRHO"] = f"using G_mRRHO from {using_part}" + ensembledata.si["part4"]["G_solv"] = f"using G_solv from {using_part}" + ensembledata.si["part4"]["Threshold"] = f"Boltzmann sum threshold: {boltzmannthr} %" + if config.prog4_j == config.prog4_s: + ensembledata.si["part4"]["main QM code"] = str(config.prog4_j).upper() + else: + ensembledata.si["part4"]["main QM code"] = ( + str(config.prog4_j).upper() + " " + str(config.prog4_s).upper() + ) + # END SI information-------------------------------------------------------- for conf in calculate: conf.calc_free_energy( @@ -673,8 +740,8 @@ def part4(config, conformers, store_confs, ensembledata): solv=gsolv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) try: minfree = min([i.free_energy for i in calculate if i is not None]) @@ -692,10 +759,8 @@ def part4(config, conformers, store_confs, ensembledata): lambda conf: "CONF" + str(getattr(conf, "id")), lambda conf: getattr(conf, energy)["energy"], lambda conf: getattr(conf, gsolv)["energy"], - #lambda conf: getattr(conf, rrho)["energy"], - lambda conf: conf.get_mrrho( - config.temperature, rrho, config.consider_sym - ), + # lambda conf: getattr(conf, rrho)["energy"], + lambda conf: conf.get_mrrho(config.temperature, rrho, config.consider_sym), lambda conf: getattr(conf, "free_energy"), lambda conf: getattr(conf, "rel_free_energy"), lambda conf: getattr(conf, "bm_weight") * 100, @@ -781,7 +846,7 @@ def part4(config, conformers, store_confs, ensembledata): print(f"{'ERROR:':{WARNLEN}}can't copy geometry!") store_confs.append(calculate.pop(calculate.index(conf))) elif config.part1: - # do not use coord from folder config.func it could be optimized if + # do not use coord from folder config.func it could be optimized if # part2 has ever been run, take coord from ensemble file # write coord to folder calculate, store_confs, save_errors = ensemble2coord( @@ -794,7 +859,9 @@ def part4(config, conformers, store_confs, ensembledata): if getattr(conf, "nmr_coupling_info")["info"] == "calculated": prev_calculated.append(calculate.pop(calculate.index(conf))) elif getattr(conf, "nmr_coupling_info")["info"] == "failed": - print(f"{'INFORMATION:':{WARNLEN}}The calculation failed for CONF{conf.id} in the previous run.") + print( + f"{'INFORMATION:':{WARNLEN}}The calculation failed for CONF{conf.id} in the previous run." + ) store_confs.append(calculate.pop(calculate.index(conf))) else: # still in calculate @@ -821,11 +888,11 @@ def part4(config, conformers, store_confs, ensembledata): "f_active": config.f_active, "p_active": config.p_active, "si_active": config.si_active, + "onlyread": config.onlyread, } if config.prog4_j == "orca": job = OrcaJob - instruction_j["progpath"] = config.external_paths["orcapath"] - instruction_j["prepinfo"].extend(['nmrJ',]) + instruction_j["prepinfo"].extend(["nmrJ"]) instruction_j["method"], _ = config.get_method_name( instruction_j["jobtype"], func=instruction_j["func"], @@ -836,7 +903,6 @@ def part4(config, conformers, store_confs, ensembledata): ) elif config.prog4_j == "tm": job = TmJob - instruction_j["progpath"] = config.external_paths["escfpath"] instruction_j["method"], _ = config.get_method_name( instruction_j["jobtype"], func=instruction_j["func"], @@ -855,6 +921,7 @@ def part4(config, conformers, store_confs, ensembledata): solvent=instruction_j["solvent"], prog=config.prog4_j, ) + check = {True: "was successful", False: "FAILED"} pl = config.lenconfx + 4 + len(str("/" + folder)) if calculate: @@ -956,6 +1023,8 @@ def part4(config, conformers, store_confs, ensembledata): "f_active": config.f_active, "p_active": config.p_active, "si_active": config.si_active, + "onlyread": config.onlyread, + "moread": None, } if config.basis_j != config.basis_s: @@ -968,9 +1037,12 @@ def part4(config, conformers, store_confs, ensembledata): instruction_s["prepinfo"] = ["high+"] if (config.basis_j == config.basis_s) and (config.prog4_j == config.prog4_s): if config.func_j == config.func_s and config.couplings: - instruction_s["prepinfo"] = [] # don't do single-point + if config.prog4_s == "tm": + instruction_s["prepinfo"] = [] instruction_s["jobtype"] = "shieldings" + if config.prog4_s == "orca": + instruction_s["moread"] = ["! MORead", '%moinp "inpJ.gbw"'] elif config.func_j != config.func_s and config.couplings: instruction_s["prepinfo"] = ["high+"] # use already converged mos as start mos @@ -980,8 +1052,7 @@ def part4(config, conformers, store_confs, ensembledata): if config.prog4_s == "orca": job = OrcaJob - instruction_s["prepinfo"].extend(['nmrS',]) - instruction_s["progpath"] = config.external_paths["orcapath"] + instruction_s["prepinfo"].extend(["nmrS"]) instruction_s["method"], _ = config.get_method_name( instruction_s["jobtype"], func=instruction_s["func"], @@ -992,7 +1063,6 @@ def part4(config, conformers, store_confs, ensembledata): ) elif config.prog4_s == "tm": job = TmJob - instruction_s["progpath"] = config.external_paths["mpshiftpath"] instruction_s["method"], _ = config.get_method_name( instruction_s["jobtype"], func=instruction_s["func"], @@ -1071,11 +1141,48 @@ def part4(config, conformers, store_confs, ensembledata): print(line) calculate.append(prev_calculated.pop(prev_calculated.index(conf))) + # SI update----------------------------------------------------------------- + if config.couplings: + ensembledata.si["part4"]["Coupling constants"] = instruction_j["method"] + if config.prog4_j == "tm": + ensembledata.si["part4"]["Coupling constants"] += " " + " ".join( + qm_prepinfo["tm"][instruction_j["prepinfo"][0]] + ).replace("-", "") + elif config.prog4_j == "orca": + ensembledata.si["part4"]["Coupling constants"] += " " + " ".join( + qm_prepinfo["orca"][instruction_j["prepinfo"][0]] + ) + else: + ensembledata.si["part4"]["Coupling constants"] = "not calculated" + if config.shieldings: + ensembledata.si["part4"]["Shielding constants"] = instruction_s["method"] + if config.prog4_s == "tm": + if ( + config.prog4_j == "tm" + and config.func_s == config.func_j + and config.basis_s == config.basis_j + and config.couplings + ): + ensembledata.si["part4"]["Shielding constants"] += " " + " ".join( + qm_prepinfo["tm"][instruction_j["prepinfo"][0]] + ).replace("-", "") + else: + ensembledata.si["part4"]["Shielding constants"] += " " + " ".join( + qm_prepinfo["tm"][instruction_s["prepinfo"][0]] + ).replace("-", "") + elif config.prog4_s == "orca": + ensembledata.si["part4"]["Shielding constants"] += " " + " ".join( + qm_prepinfo["orca"][instruction_s["prepinfo"][0]] + ) + else: + ensembledata.si["part4"]["Shielding constants"] = "not calculated" + # END SI update-------------------------------------------------------------- + if not calculate: print(f"{'ERROR:':{WARNLEN}}No conformers left!") print("Going to exit!") sys.exit(1) - + # write anmr_enso output! print("\nGenerating file anmr_enso for processing with the ANMR program.") for conf in calculate: @@ -1084,20 +1191,30 @@ def part4(config, conformers, store_confs, ensembledata): solv=gsolv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) try: length = max([str(i.id) for i in calculate]) if int(length) < 4: length = 4 fmtenergy = max([len("{: .7f}".format(i.free_energy)) for i in calculate]) - if config.solvent != 'gas': - fmtsolv = max([len("{: .7f}".format(getattr(i, gsolv,{'energy': 0.0})["energy"])) for i in calculate]) + if config.solvent != "gas": + fmtsolv = max( + [ + len("{: .7f}".format(getattr(i, gsolv, {"energy": 0.0})["energy"])) + for i in calculate + ] + ) else: fmtsolv = 10 if config.evaluate_rrho: - fmtrrho = max([len("{: .7f}".format(getattr(i, rrho,{'energy': 0.0})["energy"])) for i in calculate]) + fmtrrho = max( + [ + len("{: .7f}".format(getattr(i, rrho, {"energy": 0.0})["energy"])) + for i in calculate + ] + ) else: fmtrrho = 10 except Exception as e: @@ -1116,8 +1233,8 @@ def part4(config, conformers, store_confs, ensembledata): f"{1:<5} {conf.id:^{length}} {conf.id:^{length}} " f"{conf.bm_weight: {2}.4f} {getattr(conf, energy)['energy']: {fmtenergy}.{7}f} " f"{getattr(conf, str(gsolv), {'energy': 0.0})['energy']: {fmtsolv-7}.{7}f} " - f"{conf.get_mrrho(config.temperature, rrho, config.consider_sym)} " - #f"{getattr(conf, str(rrho), {'energy': 0.0})['energy']: {fmtrrho-7}.{7}f} " + f"{conf.get_mrrho(config.temperature, rrho, config.consider_sym): {fmtrrho-7}.{7}f} " + # f"{getattr(conf, str(rrho), {'energy': 0.0})['energy']: {fmtrrho-7}.{7}f} " f"{conf.gi:.3f}\n" ) @@ -1125,23 +1242,28 @@ def part4(config, conformers, store_confs, ensembledata): print("\nWriting .anmrrc!") element_ref_shield = write_anmrrc(config) - print("\nGenerating plain nmrprop.dat files for each populated conformer.") - print("These files contain all calculated shielding and coupling constants.") - print("The files can be read by ANMR using the keyword '-plain'.\n") - # write generic: - instructgeneric = {"jobtype": "genericout", "nat": int(config.nat)} - calculate = run_in_parallel( - config, - q, - resultq, - job, - config.maxthreads, - config.omp, - calculate, - instructgeneric, - False, - folder - ) + if config.prog4_s == config.prog4_j: + print("\nGenerating plain nmrprop.dat files for each populated conformer.") + print("These files contain all calculated shielding and coupling constants.") + print("The files can be read by ANMR using the keyword '-plain'.\n") + # write generic: + instructgeneric = { + "jobtype": "genericout", + "nat": int(config.nat), + "onlyread": config.onlyread, + } + calculate = run_in_parallel( + config, + q, + resultq, + job, + config.maxthreads, + config.omp, + calculate, + instructgeneric, + False, + folder, + ) # printout the averaged shielding constants if config.shieldings: @@ -1152,17 +1274,11 @@ def part4(config, conformers, store_confs, ensembledata): conf.reset_job_info() if save_errors: - print( - "\n***---------------------------------------------------------***" - ) - print( - "Printing most relevant errors again, just for user convenience:" - ) + print("\n***---------------------------------------------------------***") + print("Printing most relevant errors again, just for user convenience:") for _ in list(save_errors): print(save_errors.pop()) - print( - "***---------------------------------------------------------***" - ) + print("***---------------------------------------------------------***") # end printout for part4 tmp = int((PLENGTH - len("END of Part4")) / 2) print("\n" + "".ljust(tmp, ">") + "END of Part4" + "".rjust(tmp, "<")) diff --git a/censo_qm/opticalrotation.py b/censo_qm/opticalrotation.py index 95a1b45..c9069e3 100644 --- a/censo_qm/opticalrotation.py +++ b/censo_qm/opticalrotation.py @@ -56,20 +56,19 @@ def part5(config, conformers, store_confs, ensembledata): max_len_digilen = 0 for item in info: - if item[0] == 'justprint': + if item[0] == "justprint": if "short-notation" in item[1]: - tmp = len(item[1]) -len('short-notation:') + tmp = len(item[1]) - len("short-notation:") else: tmp = len(item[1]) else: tmp = len(item[1]) if tmp > max_len_digilen: max_len_digilen = tmp - max_len_digilen +=1 + max_len_digilen += 1 if max_len_digilen < DIGILEN: max_len_digilen = DIGILEN - optionsexchange = {True: "on", False: "off"} for item in info: if item[0] == "justprint": @@ -140,7 +139,7 @@ def part5(config, conformers, store_confs, ensembledata): boltzmannthr = config.part2_threshold elif config.part1: # part2 is not calculated use boltzmann weights directly from part1 - #--> misappropriate config.part2_threshold + # --> misappropriate config.part2_threshold # This means starting from not DFT optimized geometries! energy = "prescreening_sp_info" rrho = "prescreening_grrho_info" @@ -155,19 +154,31 @@ def part5(config, conformers, store_confs, ensembledata): elif getattr(conf, rrho)["info"] != "calculated" and config.evaluate_rrho: store_confs.append(mol) continue - elif ( - getattr(conf, gsolv)["info"] != "calculated" and config.solvent != "gas" - ): + elif getattr(conf, gsolv)["info"] != "calculated" and config.solvent != "gas": store_confs.append(mol) continue else: calculate.append(mol) if unoptimized_warning: - print_errors(f"{'INFORMATION:':{WARNLEN}}Conformers have not been optimized at DFT level!!!\n" - f"{'':{WARNLEN}}Use results with care!\n", save_errors + print_errors( + f"{'INFORMATION:':{WARNLEN}}Conformers have not been optimized at DFT level!!!\n" + f"{'':{WARNLEN}}Use results with care!\n", + save_errors, ) - + # SI geometry: + if unoptimized_warning: + ensembledata.si["part5"]["Geometry"] = "GFNn-xTB (input geometry)" + else: + ensembledata.si["part5"]["Geometry"], _ = config.get_method_name( + "xtbopt", + func=config.func, + basis=config.basis, + solvent=config.solvent, + sm=config.sm2, + ) + ensembledata.si["part5"]["Geometry"] += f" @optlevel: {config.optlevel2}" + # if not calculate and not prev_calculated: print(f"{'ERROR:':{WARNLEN}}No conformers left!") print("Going to exit!") @@ -179,6 +190,7 @@ def part5(config, conformers, store_confs, ensembledata): # Calculate Boltzmann weight for confs: if config.part3: + using_part = "part3 - refinement" if not config.evaluate_rrho: rrho = None rrho_method = None @@ -216,6 +228,7 @@ def part5(config, conformers, store_confs, ensembledata): solvent=config.solvent, ) elif config.part2: + using_part = "part2 - optimization" if not config.evaluate_rrho: rrho = None rrho_method = None @@ -253,6 +266,7 @@ def part5(config, conformers, store_confs, ensembledata): solvent=config.solvent, ) elif config.part1: + using_part = "part1 - prescreening" # on DFT unoptimized geometries! if not config.evaluate_rrho: rrho = None @@ -291,14 +305,23 @@ def part5(config, conformers, store_confs, ensembledata): solvent=config.solvent, ) + # SI information:------------------------------------------------------- + ensembledata.si["part5"]["Energy"] = f"using Energy from {using_part}" + ensembledata.si["part5"]["Energy_settings"] = f"see {using_part}" + ensembledata.si["part5"]["G_mRRHO"] = f"using G_mRRHO from {using_part}" + ensembledata.si["part5"]["G_solv"] = f"using G_solv from {using_part}" + ensembledata.si["part5"]["Threshold"] = f"Boltzmann sum threshold: {boltzmannthr} %" + ensembledata.si["part5"]["main QM code"] = str("tm").upper() + # END SI information-------------------------------------------------------- + for conf in calculate: conf.calc_free_energy( e=energy, solv=gsolv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) try: minfree = min([i.free_energy for i in calculate if i is not None]) @@ -316,12 +339,8 @@ def part5(config, conformers, store_confs, ensembledata): lambda conf: "CONF" + str(getattr(conf, "id")), lambda conf: getattr(conf, energy)["energy"], lambda conf: getattr(conf, gsolv)["energy"], - #lambda conf: getattr(conf, rrho)["energy"], - lambda conf: conf.get_mrrho( - config.temperature, - rrho, - config.consider_sym - ), + # lambda conf: getattr(conf, rrho)["energy"], + lambda conf: conf.get_mrrho(config.temperature, rrho, config.consider_sym), lambda conf: getattr(conf, "free_energy"), lambda conf: getattr(conf, "rel_free_energy"), lambda conf: getattr(conf, "bm_weight") * 100, @@ -408,7 +427,7 @@ def part5(config, conformers, store_confs, ensembledata): print(f"{'ERROR:':{WARNLEN}}can't copy geometry!") store_confs.append(calculate.pop(calculate.index(conf))) elif config.part1: - # do not use coord from folder config.func it could be optimized if + # do not use coord from folder config.func it could be optimized if # part2 has ever been run, take coord from ensemble file # write coord to folder calculate, store_confs, save_errors = ensemble2coord( @@ -425,7 +444,9 @@ def part5(config, conformers, store_confs, ensembledata): if getattr(conf, "optical_rotation_info")["info"] == "calculated": prev_calculated.append(calculate.pop(calculate.index(conf))) elif getattr(conf, "optical_rotation_info")["info"] == "failed": - print(f"{'INFORMATION:':{WARNLEN}}The calculation failed for CONF{conf.id} in the previous run.") + print( + f"{'INFORMATION:':{WARNLEN}}The calculation failed for CONF{conf.id} in the previous run." + ) store_confs.append(calculate.pop(calculate.index(conf))) instruction_or = { @@ -445,19 +466,19 @@ def part5(config, conformers, store_confs, ensembledata): "freq_or": config.freq_or, } if config.prog == "orca": - print(f"{'ERROR:':{WARNLEN}}Can't calculate OR with ORCA! You can use TM instead.") + print( + f"{'ERROR:':{WARNLEN}}Can't calculate OR with ORCA! You can use TM instead." + ) print("Going to exit!") sys.exit(1) # ORCA can't calculate optical rotation!!! --> job = OrcaJob - instruction_or["progpath"] = config.external_paths["orcapath"] if config.solvent != "gas": instruction_or["solvent"] = config.solvent instruction_or["sm"] = "cpcm" if config.prog == "tm": job = TmJob instruction_or["prepinfo"] = ["clear", "-grid", "2", "-scfconv", "6"] - instruction_or["progpath"] = config.external_paths["escfpath"] if config.basis == config.basis_or: instruction_or["copymos"] = config.func instruction_or["jobtype"] = "opt-rot" @@ -474,6 +495,12 @@ def part5(config, conformers, store_confs, ensembledata): func2=instruction_or["func"], sm=instruction_or["sm"], ) + # SI update----------------------------------------------------------------- + ensembledata.si["part5"]["Optical rotation"] = ( + instruction_or["method"] + + f" using {' '.join(instruction_or['prepinfo']).replace('-', '').replace('clear', '')}" + ) + # END SI update-------------------------------------------------------------- print(f"\nOptical-rotation is calculated at {instruction_or['method']} level.\n") check = {True: "was successful", False: "FAILED"} @@ -552,41 +579,50 @@ def part5(config, conformers, store_confs, ensembledata): ) calculate.append(prev_calculated.pop(prev_calculated.index(conf))) - outlen = config.lenconfx+1 - if len('ENSEMBLE') >= config.lenconfx: - outlen = len('ENSEMBLE')+1 + outlen = config.lenconfx + 1 + if len("ENSEMBLE") >= config.lenconfx: + outlen = len("ENSEMBLE") + 1 for freq in config.freq_or: averaged_or = 0.0 - with open(os.path.join(config.cwd, f"OR_{int(freq)}.dat"), "w", newline=None) as outdata: - outdata.write(f"{'#label':{outlen}} {'#unmod_alpha':>{max_fmt}} {'#%pop':^7} {'#pop_alpha':{max_fmt}} \n") + with open( + os.path.join(config.cwd, f"OR_{int(freq)}.dat"), "w", newline=None + ) as outdata: + outdata.write( + f"{'#label':{outlen}} {'#unmod_alpha':>{max_fmt}} {'#%pop':^7} {'#pop_alpha':{max_fmt}} \n" + ) for conf in calculate: averaged_or += conf.bm_weight * conf.optical_rotation_info["range"].get( freq, 0.0 ) # CONFX outdata.write( - f"{'CONF'+str(conf.id):{outlen}} "+ - f"{conf.optical_rotation_info['range'].get(freq, 0.0):> {max_fmt}.7f} "+ - f"{conf.bm_weight*100:> 6.2f} "+ - f"{conf.optical_rotation_info['range'].get(freq, 0.0)*conf.bm_weight:> {max_fmt}.7f}\n" - ) + f"{'CONF'+str(conf.id):{outlen}} " + + f"{conf.optical_rotation_info['range'].get(freq, 0.0):> {max_fmt}.7f} " + + f"{conf.bm_weight*100:> 6.2f} " + + f"{conf.optical_rotation_info['range'].get(freq, 0.0)*conf.bm_weight:> {max_fmt}.7f}\n" + ) # ENSEMBLE outdata.write( - f"{'ENSEMBLE':{outlen}} "+ - f"{'-':^{max_fmt-1}} "+ - f"{100.00:> 6.2f} "+ - f"{averaged_or:> {max_fmt}.7f}\n" - ) + f"{'ENSEMBLE':{outlen}} " + + f"{'-':^{max_fmt-1}} " + + f"{100.00:> 6.2f} " + + f"{averaged_or:> {max_fmt}.7f}\n" + ) print( f"\nAveraged specific rotation at {freq} nm : " f"{averaged_or:> {max_fmt}.3f} in deg*[dm(g/cc)]^(-1)" ) - - if all( - [conf.lowlevel_gsolv_compare_info["std_dev"] is not None for conf in calculate] - ) and calculate: + if ( + all( + [ + conf.lowlevel_gsolv_compare_info["std_dev"] is not None + for conf in calculate + ] + ) + and calculate + ): for freq in config.freq_or: all_or = [] for _ in range(1000): @@ -597,8 +633,8 @@ def part5(config, conformers, store_confs, ensembledata): solv=gsolv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) conf.free_energy += normalvariate( 0.0, conf.lowlevel_gsolv_compare_info["std_dev"] ) @@ -639,15 +675,13 @@ def part5(config, conformers, store_confs, ensembledata): averaged_or = 0.0 for conf in calculate: conf.calc_free_energy( - e=energy, - solv=gsolv, + e=energy, + solv=gsolv, rrho=rrho, t=config.consider_sym, - consider_sym=config.consider_sym - ) - conf.free_energy += normalvariate( - 0.0, (0.4/AU2KCAL) + consider_sym=config.consider_sym, ) + conf.free_energy += normalvariate(0.0, (0.4 / AU2KCAL)) calculate = calc_boltzmannweights( calculate, "free_energy", config.temperature ) @@ -677,7 +711,28 @@ def part5(config, conformers, store_confs, ensembledata): f" SD based on SD in G of 0.4 kcal/mol" f": {calc_std_dev(all_or):> {max_fmt}.3f} in deg*[dm(g/cc)]^(-1)" ) - + if calculate: + # print min and max OR values and corresponding conformers, but no Boltzmann weighting + print("") + output = [] + output.append(f"{'frequency'} {'min(OR)'} {'CONF'} {'max(OR)'} {'CONF'}") + output.append("".ljust(int(50), "-")) + for freq in config.freq_or: + tmp_or = {} + for conf in calculate: + tmp_or[conf.id] = conf.optical_rotation_info["range"].get(freq) + minx = min(list(tmp_or.values())) + minid = next(key for key, value in tmp_or.items() if value == minx) + maxx = max(list(tmp_or.values())) + maxid = next(key for key, value in tmp_or.items() if value == maxx) + + output.append( + f"{float(freq): ^12} {minx: .2f} {'CONF'+str(minid)} {maxx: .2f} {'CONF'+str(maxid)}" + ) + output.append("".ljust(int(50), "-")) + output.append("* min and max values are not Boltzmann weighted.") + for line in output: + print(line) for conf in calculate: conf.reset_job_info() @@ -693,17 +748,11 @@ def part5(config, conformers, store_confs, ensembledata): ) if save_errors: - print( - "\n***---------------------------------------------------------***" - ) - print( - "Printing most relevant errors again, just for user convenience:" - ) + print("\n***---------------------------------------------------------***") + print("Printing most relevant errors again, just for user convenience:") for _ in list(save_errors): print(save_errors.pop()) - print( - "***---------------------------------------------------------***" - ) + print("***---------------------------------------------------------***") if not calculate: print("ERROR: No conformers left!") print("Going to exit!") diff --git a/censo_qm/optimization.py b/censo_qm/optimization.py index 4ec5574..30fa945 100755 --- a/censo_qm/optimization.py +++ b/censo_qm/optimization.py @@ -8,7 +8,7 @@ import os import sys from copy import deepcopy -from .cfg import PLENGTH, CODING, AU2KCAL, DIGILEN, WARNLEN +from .cfg import PLENGTH, CODING, AU2KCAL, DIGILEN, WARNLEN, qm_prepinfo from .utilities import ( check_for_folder, print_block, @@ -26,6 +26,7 @@ print, print_errors, calc_weighted_std_dev, + conf_in_interval, ) from .orca_job import OrcaJob from .tm_job import TmJob @@ -88,7 +89,9 @@ def part2(config, conformers, store_confs, ensembledata): ], ] ) - info.append(["part2_threshold", "Boltzmann sum threshold G_thr(2) for sorting in part2"]) + info.append( + ["part2_threshold", "Boltzmann sum threshold G_thr(2) for sorting in part2"] + ) info.append(["evaluate_rrho", "calculate mRRHO contribution"]) if config.evaluate_rrho: info.append(["prog_rrho", "program for mRRHO contribution"]) @@ -103,16 +106,16 @@ def part2(config, conformers, store_confs, ensembledata): ) max_len_digilen = 0 for item in info: - if item[0] == 'justprint': + if item[0] == "justprint": if "short-notation" in item[1]: - tmp = len(item[1]) -len('short-notation:') + tmp = len(item[1]) - len("short-notation:") else: tmp = len(item[1]) else: tmp = len(item[1]) if tmp > max_len_digilen: max_len_digilen = tmp - max_len_digilen +=1 + max_len_digilen += 1 if max_len_digilen < DIGILEN: max_len_digilen = DIGILEN @@ -228,6 +231,7 @@ def part2(config, conformers, store_confs, ensembledata): "energy": 0.0, "energy2": 0.0, "success": False, + "onlyread": config.onlyread, } # INSTRUCTION OPT !!!! @@ -247,6 +251,7 @@ def part2(config, conformers, store_confs, ensembledata): "energy": 0.0, "energy2": 0.0, "success": False, + "onlyread": config.onlyread, } instruction_rrho_crude = { @@ -257,18 +262,17 @@ def part2(config, conformers, store_confs, ensembledata): "charge": config.charge, "unpaired": config.unpaired, "solvent": config.solvent, - "progpath": config.external_paths["xtbpath"], "bhess": config.bhess, "sm_rrho": config.sm_rrho, "rmsdbias": config.rmsdbias, "cwd": config.cwd, - "consider_sym": config.consider_sym, "energy": 0.0, "energy2": 0.0, "success": False, "imagthr": config.imagthr, "sthr": config.sthr, - "scale":config.scale, + "scale": config.scale, + "onlyread": config.onlyread, } instruction_rrho_crude["method"], _ = config.get_method_name( "rrhoxtb", @@ -312,7 +316,17 @@ def part2(config, conformers, store_confs, ensembledata): solvent=instruction_opt["solvent"], sm=instruction_opt["sm"], ) - + # SI + if config.prog == "tm": + ensembledata.si["part2"]["Geometry"] = ( + instruction_opt["method"] + + f" @optlevel: {config.optlevel2} using {' '.join(qm_prepinfo['tm'][instruction_opt['prepinfo'][0]]).replace('-', '')}" + ) + elif config.prog == "orca": + ensembledata.si["part2"]["Geometry"] = ( + instruction_opt["method"] + + f" @optlevel: {config.optlevel2} using {' '.join(qm_prepinfo['orca'][instruction_opt['prepinfo'][0]])}" + ) check = {True: "was successful", False: "FAILED"} if calculate: print("The optimization is calculated for:") @@ -374,38 +388,39 @@ def part2(config, conformers, store_confs, ensembledata): print("Starting optimizations".center(70, "*")) ### settings in instruction_opt are overwriting conf.job everytime,(while loop) ### therefore dont write information which has to be reaccessed to it! - run = 1 timings = [] # time per cycle cycle_spearman = [] # spearmanthr used in evaluation per cycle nconf_cycle = [] # number of conformers at end of each cycle - do_increase = 0.6 - if config.opt_limit * do_increase >= 1.5: - ewin_increase = config.opt_limit * do_increase - print( - f"\nStarting threshold is set to {config.opt_limit} + " - f"{do_increase*100} % = {config.opt_limit + ewin_increase} kcal/mol\n" - ) - else: - ewin_increase = 1.5 - print( - f"\nStarting threshold is set to {config.opt_limit} + " - f"{ewin_increase} kcal/mol = {config.opt_limit + ewin_increase} kcal/mol\n" - ) - ewin_initial = config.opt_limit + ewin_increase - ewin = config.opt_limit + ewin_increase - - print(f"Lower limit is set to G_thr(opt,2) = {config.opt_limit} kcal/mol\n") - lower_limit = config.opt_limit - maxecyc_prev = 1 - maxecyc = 0 - converged_run1 = [] - if config.nat > 200: - # stopcycle = don't optimize more than stopcycle cycles - stopcycle = config.nat * 2 - else: - stopcycle = 200 if config.opt_spearman: + run = 1 + do_increase = 0.6 + if config.opt_limit * do_increase >= 1.5: + ewin_increase = config.opt_limit * do_increase + print( + f"\nStarting threshold is set to {config.opt_limit} + " + f"{do_increase*100} % = {config.opt_limit + ewin_increase} kcal/mol\n" + ) + else: + ewin_increase = 1.5 + print( + f"\nStarting threshold is set to {config.opt_limit} + " + f"{ewin_increase} kcal/mol = {config.opt_limit + ewin_increase} kcal/mol\n" + ) + ewin_initial = config.opt_limit + ewin_increase + ewin = config.opt_limit + ewin_increase + + print(f"Lower limit is set to G_thr(opt,2) = {config.opt_limit} kcal/mol\n") + lower_limit = config.opt_limit + maxecyc_prev = 1 + maxecyc = 0 + converged_run1 = [] + if config.nat > 200: + # stopcycle = don't optimize more than stopcycle cycles + stopcycle = config.nat * 2 + else: + stopcycle = 200 + # start batch calculations: while calculate: tic = time.perf_counter() print(f"CYCLE {str(run)}".center(70, "*")) @@ -499,13 +514,18 @@ def part2(config, conformers, store_confs, ensembledata): if not os.path.isfile(os.path.join(tmp_from, "coord")): print_errors( f"{'ERROR:':{WARNLEN}}while copying the coord file from {tmp_from}! " - "The corresponding file does not exist.", save_errors + "The corresponding file does not exist.", + save_errors, ) elif not os.path.isdir(tmp_to): print_errors( - f"{'ERROR:':{WARNLEN}}Could not create folder {tmp_to}!", save_errors + f"{'ERROR:':{WARNLEN}}Could not create folder {tmp_to}!", + save_errors, ) - print_errors(f"{'ERROR:':{WARNLEN}}Removing conformer {conf.name}!", save_errors) + print_errors( + f"{'ERROR:':{WARNLEN}}Removing conformer {conf.name}!", + save_errors, + ) conf.lowlevel_grrho_info["info"] = "prep-failed" store_confs.append(calculate.pop(calculate.index(conf))) # parallel execution: @@ -550,8 +570,8 @@ def part2(config, conformers, store_confs, ensembledata): "method_rrho" ] = instruction_rrho_crude["method"] conf.optimization_info["info_rrho"] = "calculated" - conf.symnum = conf.job['symnum'] - conf.sym = conf.job['symmetry'] + conf.symnum = conf.job["symnum"] + conf.sym = conf.job["symmetry"] for conf in list(calculate): if conf.id in converged_run1: @@ -592,17 +612,6 @@ def part2(config, conformers, store_confs, ensembledata): print("Spearman rank evaluation is performed in the next cycle.") cycle_spearman.append("") run_spearman = False - # elif run == 1 and not gesc: - # # only evaluate spearman starting from second cycle - # print("Spearman rank evaluation is performed in the next cycle.") - # cycle_spearman.append("") - # run_spearman = False - # run += 1 - # toc = time.perf_counter() - # timings.append(toc - tic) - # nconf_cycle.append(len(calculate) + len(prev_calculated)) - # print(f"CYCLE {run} performed in {toc -tic:0.4f} seconds") - # continue # lists of equal lenght: for conf in sorted( @@ -615,7 +624,7 @@ def part2(config, conformers, store_confs, ensembledata): # calculate min of each cycle: minecyc = [] if config.evaluate_rrho: - #rrho_energy = "energy_rrho" + # rrho_energy = "energy_rrho" rrho_energy = "rrho_optimization" else: rrho_energy = None # to get 0.0 contribution @@ -630,9 +639,9 @@ def part2(config, conformers, store_confs, ensembledata): rrho=rrho_energy, consider_sym=config.consider_sym, ) - #+ getattr(conf, "optimization_info").get( + # + getattr(conf, "optimization_info").get( # rrho_energy, 0.0 - #) + # ) for conf in calculate + prev_calculated if conf.job["ecyc"][i] is not None ] @@ -655,9 +664,9 @@ def part2(config, conformers, store_confs, ensembledata): rrho=rrho_energy, consider_sym=config.consider_sym, ) - #+ getattr(conf, "optimization_info").get( + # + getattr(conf, "optimization_info").get( # rrho_energy, 0.0 - #) + # ) - minecyc[i] ) * AU2KCAL @@ -740,7 +749,7 @@ def part2(config, conformers, store_confs, ensembledata): f"Final averaged Spearman correlation coefficient: " f"{(sum(evalspearman)/2):>.4f}" ) - + cycle_spearman.append(f"{sum(evalspearman)/2:.3f}") if ( len(evalspearman) >= 2 and sum(evalspearman) / 2 >= config.spearmanthr @@ -759,7 +768,6 @@ def part2(config, conformers, store_confs, ensembledata): print( f"Current optimization threshold: {ewin:.2f} kcal/mol" ) - cycle_spearman.append(f"{sum(evalspearman)/2:.3f}") for conf in list(calculate): if conf.job["decyc"][-1] > ewin and conf.job["grad_norm"] < 0.01: @@ -814,7 +822,9 @@ def part2(config, conformers, store_confs, ensembledata): conf.optimization_info["decyc"] = conf.job["decyc"] conf.optimization_info[ "convergence" - ] = "converged" #### THIS IS NOT CORRECT /or can lead to errors! + ] = ( + "converged" + ) #### THIS IS NOT CORRECT /or can lead to errors! conf.optimization_info["method"] = instruction_opt["method"] prev_calculated.append(calculate.pop(calculate.index(conf))) # @@ -865,9 +875,7 @@ def part2(config, conformers, store_confs, ensembledata): timings.append(toc - tic) # ********************end standard optimization ********************* print("Finished optimizations!".center(70, "*")) - if config.opt_spearman and ( - len(timings) == len(cycle_spearman) == len(nconf_cycle) - ): + if config.opt_spearman and (len(timings) == len(nconf_cycle)): try: tl = max([len(f"{i: .2f}") for i in timings]) if tl > 7: @@ -877,9 +885,13 @@ def part2(config, conformers, store_confs, ensembledata): print("Timings:") print(f"Cycle: [s] {'#nconfs':^{tmp1}} Spearman coeff.") for i in range(len(timings)): - print( - f"{i+1:>4} {timings[i]:> {tl}.2f} {nconf_cycle[i]:^{tmp1}} {cycle_spearman[i]}" - ) + try: + print( + f"{i+1:>4} {timings[i]:> {tl}.2f} {nconf_cycle[i]:^{tmp1}} {cycle_spearman[i]}" + ) + except Exception as e: + print(i, len(nconf_cycle), len(cycle_spearman)) + print(f"{i+1:>4} {timings[i]:> {tl}.2f} {'-':^{tmp1}} {'-'}") print("sum: {:> .2f}".format(sum(timings))) except Exception as e: print(e) @@ -939,6 +951,7 @@ def part2(config, conformers, store_confs, ensembledata): "energy2": 0.0, "success": False, "gfn_version": config.part2_gfnv, + "onlyread": config.onlyread, } if config.multitemp: instruction_gsolv["trange"] = [ @@ -965,7 +978,6 @@ def part2(config, conformers, store_confs, ensembledata): # COSMO-RS if "cosmors" in config.smgsolv2 and config.smgsolv2 != "dcosmors": job = TmJob - instruction_gsolv["prepinfo"] = ["low+"] exc_fine = {"cosmors": "normal", "cosmors-fine": "fine"} tmp = { "jobtype": "cosmors", @@ -989,9 +1001,6 @@ def part2(config, conformers, store_confs, ensembledata): # GBSA-Gsolv / ALPB-Gsolv elif config.smgsolv2 in ("gbsa_gsolv", "alpb_gsolv"): instruction_gsolv["jobtype"] = instruction_gsolv["sm"] - if config.prog == "orca": - instruction_gsolv["progpath"] = config.external_paths["orcapath"] - instruction_gsolv["xtb_driver_path"] = config.external_paths["xtbpath"] instruction_gsolv["method"], instruction_gsolv[ "method2" ] = config.get_method_name( @@ -1015,7 +1024,6 @@ def part2(config, conformers, store_confs, ensembledata): elif config.smgsolv2 == "smd_gsolv": job = OrcaJob instruction_gsolv["jobtype"] = "smd_gsolv" - instruction_gsolv["progpath"] = config.external_paths["orcapath"] instruction_gsolv["method"], instruction_gsolv[ "method2" ] = config.get_method_name( @@ -1029,8 +1037,6 @@ def part2(config, conformers, store_confs, ensembledata): else: # with implicit solvation instruction_gsolv["jobtype"] = "sp_implicit" - if config.prog == "orca": - instruction_gsolv["progpath"] = config.external_paths["orcapath"] instruction_gsolv["method"], instruction_gsolv[ "method2" ] = config.get_method_name( @@ -1144,7 +1150,7 @@ def part2(config, conformers, store_confs, ensembledata): ) for conf in list(calculate): - if instruction_gsolv["jobtype"] in("sp", "sp_implicit"): + if instruction_gsolv["jobtype"] in ("sp", "sp_implicit"): line = ( f"{name} calculation {check[conf.job['success']]}" f" for {last_folders(conf.job['workdir'], 2):>{pl}}: " @@ -1232,6 +1238,36 @@ def part2(config, conformers, store_confs, ensembledata): f"{conf.lowlevel_gsolv_info['energy']:>.7f}" ) calculate.append(prev_calculated.pop(prev_calculated.index(conf))) + + # create data for possible SI generation:----------------------------------- + # Energy: + ensembledata.si["part2"]["Energy"] = instruction_gsolv["method"] + if "DOGCP" in instruction_gsolv["prepinfo"]: + ensembledata.si["part2"]["Energy"] += " + GCP" + # Energy_settings: + if job == TmJob: + ensembledata.si["part2"]["Energy_settings"] = " ".join( + qm_prepinfo["tm"][instruction_gsolv["prepinfo"][0]] + ).replace("-", "") + elif job == OrcaJob: + ensembledata.si["part2"]["Energy_settings"] = " ".join( + qm_prepinfo["orca"][instruction_gsolv["prepinfo"][0]] + ) + # GmRRHO is done below! + # Solvation: + if config.solvent == "gas": + ensembledata.si["part2"]["G_solv"] = "gas-phase" + else: + ensembledata.si["part2"]["G_solv"] = instruction_gsolv["method2"] + # Geometry is done above: + # QM-CODE: + ensembledata.si["part2"]["main QM code"] = str(config.prog).upper() + # Threshold: + ensembledata.si["part2"][ + "Threshold" + ] = f"Opt_limit: {config.opt_limit} kcal/mol, Boltzmann sum threshold: {config.part2_threshold} %" + # END SI generation -------------------------------------------------------- + # reset for conf in calculate: conf.reset_job_info() @@ -1265,6 +1301,9 @@ def part2(config, conformers, store_confs, ensembledata): conf = calculate.pop(calculate.index(conf)) conf.__class__ = job conf.job["success"] = True + conf.sym = conf.lowlevel_grrho_info.get("sym", conf.sym) + conf.linear = conf.lowlevel_grrho_info.get("linear", conf.linear) + conf.symnum = conf.lowlevel_grrho_info.get("symnum", conf.symnum) prev_calculated.append(conf) if not calculate and not prev_calculated: @@ -1285,7 +1324,6 @@ def part2(config, conformers, store_confs, ensembledata): "charge": config.charge, "unpaired": config.unpaired, "solvent": config.solvent, - "progpath": config.external_paths["xtbpath"], "bhess": config.bhess, "sm_rrho": config.sm_rrho, "rmsdbias": config.rmsdbias, @@ -1296,7 +1334,8 @@ def part2(config, conformers, store_confs, ensembledata): "success": False, "imagthr": config.imagthr, "sthr": config.sthr, - "scale":config.scale, + "scale": config.scale, + "onlyread": config.onlyread, } instruction_rrho["method"], _ = config.get_method_name( "rrhoxtb", @@ -1305,6 +1344,16 @@ def part2(config, conformers, store_confs, ensembledata): sm=instruction_rrho["sm_rrho"], solvent=instruction_rrho["solvent"], ) + + # GmRRHO for SI information: + if config.evaluate_rrho: + ensembledata.si["part2"]["G_mRRHO"] = instruction_rrho["method"] + if "bhess" in ensembledata.si["part2"]["G_mRRHO"]: + ensembledata.si["part2"]["G_mRRHO"] += " SPH" + else: + ensembledata.si["part2"]["G_mRRHO"] = "not included" + # END SI + if config.multitemp: instruction_rrho["trange"] = [ i for i in frange(config.trange[0], config.trange[1], config.trange[2]) @@ -1336,11 +1385,17 @@ def part2(config, conformers, store_confs, ensembledata): print_errors( f"{'ERROR:':{WARNLEN}}while copying the coord file from {tmp_from}! " "The corresponding file does not exist.", - save_errors + save_errors, ) elif not os.path.isdir(tmp_to): - print_errors(f"{'ERROR:':{WARNLEN}}Could not create folder {tmp_to}!", save_errors) - print_errors(f"{'ERROR:':{WARNLEN}}Removing conformer {conf.name}!", save_errors) + print_errors( + f"{'ERROR:':{WARNLEN}}Could not create folder {tmp_to}!", + save_errors, + ) + print_errors( + f"{'ERROR:':{WARNLEN}}Removing conformer {conf.name}!", + save_errors, + ) conf.lowlevel_grrho_info["info"] = "prep-failed" store_confs.append(calculate.pop(calculate.index(conf))) # parallel execution: @@ -1383,6 +1438,9 @@ def part2(config, conformers, store_confs, ensembledata): conf.sym = conf.job["symmetry"] conf.linear = conf.job["linear"] conf.symnum = conf.job["symnum"] + conf.lowlevel_grrho_info["sym"] = conf.job["symmetry"] + conf.lowlevel_grrho_info["linear"] = conf.job["linear"] + conf.lowlevel_grrho_info["symnum"] = conf.job["symnum"] conf.lowlevel_grrho_info["rmsd"] = conf.job["rmsd"] conf.lowlevel_grrho_info["energy"] = conf.job["energy"] conf.lowlevel_grrho_info["info"] = "calculated" @@ -1436,7 +1494,7 @@ def part2(config, conformers, store_confs, ensembledata): lambda conf: getattr(conf, "rel_xtb_energy"), lambda conf: getattr(conf, "lowlevel_sp_info")["energy"], lambda conf: getattr(conf, "lowlevel_gsolv_info")["energy"], - #lambda conf: getattr(conf, "lowlevel_grrho_info")["energy"], + # lambda conf: getattr(conf, "lowlevel_grrho_info")["energy"], lambda conf: conf.get_mrrho( config.temperature, "lowlevel_grrho_info", config.consider_sym ), @@ -1514,12 +1572,12 @@ def part2(config, conformers, store_confs, ensembledata): solv = "lowlevel_gsolv_info" e = "lowlevel_sp_info" conf.calc_free_energy( - e=e, - solv=solv, + e=e, + solv=solv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) try: @@ -1544,7 +1602,8 @@ def part2(config, conformers, store_confs, ensembledata): columndescription2=columndescription2, ) # printout for part2 ------------------------------------------------------- - + conf_in_interval(calculate) + # -------------------------------------------------------------------------- # *************************************************************************** # SD on solvation: # DCOSMO-RS_GSOLV @@ -1560,10 +1619,10 @@ def part2(config, conformers, store_confs, ensembledata): "energy2": 0.0, "success": False, "gfn_version": config.part2_gfnv, + "onlyread": config.onlyread, } instruction_gsolv_compare["trange"] = [] instruction_gsolv_compare["prepinfo"] = [] - instruction_gsolv_compare["xtb_driver_path"] = config.external_paths["xtbpath"] instruction_gsolv_compare["jobtype"] = instruction_gsolv_compare["sm"] _, instruction_gsolv_compare["method"] = config.get_method_name( instruction_gsolv_compare["jobtype"], @@ -1596,7 +1655,8 @@ def part2(config, conformers, store_confs, ensembledata): ) except (TypeError, KeyError): print_errors( - f"{'ERROR:':{WARNLEN}}Can't calculate DCOSMO-RS_gsolv. Skipping SD of Gsolv!", save_errors + f"{'ERROR:':{WARNLEN}}Can't calculate DCOSMO-RS_gsolv. Skipping SD of Gsolv!", + save_errors, ) dorun = False break @@ -1634,7 +1694,10 @@ def part2(config, conformers, store_confs, ensembledata): try: shutil.copy(tmp1, tmp2) except FileNotFoundError: - print_errors(f"{'ERROR:':{WARNLEN}}can't copy optimized geometry!", save_errors) + print_errors( + f"{'ERROR:':{WARNLEN}}can't copy optimized geometry!", + save_errors, + ) if calculate: calculate = run_in_parallel( @@ -1771,8 +1834,8 @@ def part2(config, conformers, store_confs, ensembledata): solv=solv, rrho=rrho, t=temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) try: minfreeT = min( [conf.free_energy for conf in calculate if conf.free_energy is not None] @@ -1790,10 +1853,8 @@ def part2(config, conformers, store_confs, ensembledata): avG += conf.bm_weight * conf.free_energy avE += conf.bm_weight * conf.lowlevel_sp_info["energy"] avGRRHO += conf.bm_weight * conf.get_mrrho( - temperature, - "lowlevel_grrho_info", - config.consider_sym - ) + temperature, "lowlevel_grrho_info", config.consider_sym + ) avGsolv += conf.bm_weight * conf.lowlevel_gsolv_info["range"].get( temperature, 0.0 ) @@ -1850,12 +1911,12 @@ def part2(config, conformers, store_confs, ensembledata): solv = "lowlevel_gsolv_info" e = "lowlevel_sp_info" conf.calc_free_energy( - e=e, - solv=solv, + e=e, + solv=solv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) # ensembledata is used to store avGcorrection ensembledata.comment = [ lowestconf, @@ -1989,11 +2050,11 @@ def part2(config, conformers, store_confs, ensembledata): "$coord # {} {} !CONF{} \n".format( conf.free_energy, conf.get_mrrho( - config.temperature, - rrho="lowlevel_grrho_info", - consider_sym=config.consider_sym + config.temperature, + rrho="lowlevel_grrho_info", + consider_sym=config.consider_sym, ), - conf.id + conf.id, ) ) for line in coord[1:]: @@ -2010,17 +2071,11 @@ def part2(config, conformers, store_confs, ensembledata): conf.reset_job_info() if save_errors: - print( - "\n***---------------------------------------------------------***" - ) - print( - "Printing most relevant errors again, just for user convenience:" - ) + print("\n***---------------------------------------------------------***") + print("Printing most relevant errors again, just for user convenience:") for _ in list(save_errors): print(save_errors.pop()) - print( - "***---------------------------------------------------------***" - ) + print("***---------------------------------------------------------***") tmp = int((PLENGTH - len("END of Part2")) / 2) print("\n" + "".ljust(tmp, ">") + "END of Part2" + "".rjust(tmp, "<")) diff --git a/censo_qm/orca_job.py b/censo_qm/orca_job.py index 98a8289..bf471d8 100644 --- a/censo_qm/orca_job.py +++ b/censo_qm/orca_job.py @@ -5,7 +5,6 @@ import os import time import subprocess -import shutil from .cfg import ( CODING, ENVIRON, @@ -14,7 +13,7 @@ external_paths, composite_method_basis, composite_dfa, - #gga_dfa, + # gga_dfa, hybrid_dfa, dh_dfa, disp_already_included_in_func, @@ -50,10 +49,11 @@ def _prep_input(self, xyzfile=False, returndict=False): # understood commands in prepinfo: # DOGCP = uses gcp if basis known for gcp - # + # orcainput_start = OrderedDict( [ + ("moread", None), ("functional", None), ("disp", None), ("basis", None), @@ -63,16 +63,18 @@ def _prep_input(self, xyzfile=False, returndict=False): ("scfconv", None), ("frozencore", None), ("mp2", None), - ("default", [ - "! smallprint printgap noloewdin", - "! NOSOSCF", - "%MaxCore 8000", - "%output", - " print[P_BondOrder_M] 1", - " print[P_Mayer] 1", - " print[P_basis] 2", - "end", - ] + ( + "default", + [ + "! smallprint printgap noloewdin", + "! NOSOSCF", + "%MaxCore 8000", + "%output", + " print[P_BondOrder_M] 1", + " print[P_Mayer] 1", + " print[P_basis] 2", + "end", + ], ), ("job", None), ("optthreshold", None), @@ -88,19 +90,12 @@ def _prep_input(self, xyzfile=False, returndict=False): else: nmrprop = False - # definitions taken from .cfg: - # composite_dfa - # gga_dfa - # hybrid_dfa - # dh_dfa - # disp_already_included_in_func - - # build up call: orcainput = orcainput_start.copy() # set functional - if (self.job["func"] in composite_dfa and - self.job['basis'] == composite_method_basis.get(self.job["func"],'NONE')): + if self.job["func"] in composite_dfa and self.job[ + "basis" + ] == composite_method_basis.get(self.job["func"], "NONE"): orcainput["functional"] = [f"! {self.job['func']}"] else: if self.job["func"] == "kt2": @@ -117,20 +112,24 @@ def _prep_input(self, xyzfile=False, returndict=False): # set basis set orcainput["basis"] = [f"! {self.job['basis']}"] # set gcp: - if "DOGCP" in self.job['prepinfo']: - gcp_keywords = { - 'minis': "MINIS", - "sv": "SV", - "6-31g(d)": "631GD", - 'def2-sv(p)': "SV(P)", - 'def2-svp': "SVP", - 'def2-tzvp': "TZ", + if "DOGCP" in self.job["prepinfo"]: + gcp_keywords = { + "minis": "MINIS", + "sv": "SV", + "6-31g(d)": "631GD", + "def2-sv(p)": "SV(P)", + "def2-svp": "SVP", + "def2-tzvp": "TZ", } - if self.job['basis'].lower() in gcp_keywords.keys(): - if self.job['func'] in composite_dfa and self.job['basis'] == composite_method_basis.get(self.job["func"],'NONE'): + if self.job["basis"].lower() in gcp_keywords.keys(): + if self.job["func"] in composite_dfa and self.job[ + "basis" + ] == composite_method_basis.get(self.job["func"], "NONE"): pass else: - orcainput["gcp"] = [f"! GCP(DFT/{gcp_keywords[self.job['basis'].lower()]})"] + orcainput["gcp"] = [ + f"! GCP(DFT/{gcp_keywords[self.job['basis'].lower()]})" + ] # set RI def2/J, RIJCOSX def2/J gridx6 NOFINALGRIDX, RIJK def2/JK if self.job["func"] in dh_dfa: if nmrprop: @@ -169,7 +168,10 @@ def _prep_input(self, xyzfile=False, returndict=False): orcainput["grid"] = ["! grid5 nofinalgrid"] else: orcainput["grid"] = ["! grid4 nofinalgrid"] - + if self.job["moread"] is not None: + #! MORead + #%moinp "jobname2.gbw" + orcainput["moread"] = self.job["moread"] orcainput["scfconv"] = ["! scfconv6"] # set scfconv or convergence threshold e.g. loosescf or scfconv6 @@ -184,8 +186,6 @@ def _prep_input(self, xyzfile=False, returndict=False): if self.job["prepinfo"][0] in extension.keys(): orcainput["grid"] = extension[self.job["prepinfo"][0]]["grid"] orcainput["scfconv"] = extension[self.job["prepinfo"][0]]["scfconv"] - else: - pass # add dispersion if self.job["func"] not in disp_already_included_in_func: @@ -269,7 +269,7 @@ def _prep_input(self, xyzfile=False, returndict=False): orcainput["geom"] = [f"*xyz {self.job['charge']} {self.job['unpaired']+1}"] orcainput["geom"].extend(geom) orcainput["geom"].append("*") - #nmr kt2 disp + # nmr kt2 disp if self.job["func"] == "kt2" and ("nmrJ" or "nmrS" in self.job["prepinfo"]): orcainput["disp"] = [] # couplings @@ -328,10 +328,11 @@ def _prep_input(self, xyzfile=False, returndict=False): else: return tmp - def _sp(self, silent=False): + def _sp(self, silent=False, filename="sp.out"): """ ORCA input generation and single-point calculation """ + outputpath = os.path.join(self.job["workdir"], filename) if not self.job["onlyread"]: with open( os.path.join(self.job["workdir"], "inp"), "w", newline=None @@ -344,9 +345,7 @@ def _sp(self, silent=False): if not silent: print(f"Running single-point in {last_folders(self.job['workdir'], 2)}") # start SP calculation - with open( - os.path.join(self.job["workdir"], "sp.out"), "w", newline=None - ) as outputfile: + with open(outputpath, "w", newline=None) as outputfile: call = [os.path.join(external_paths["orcapath"], "orca"), "inp"] subprocess.call( call, @@ -358,14 +357,10 @@ def _sp(self, silent=False): stdout=outputfile, ) time.sleep(0.05) + # read output # check if scf is converged: - if os.path.isfile(os.path.join(self.job["workdir"], "sp.out")): - with open( - os.path.join(self.job["workdir"], "sp.out"), - "r", - encoding=CODING, - newline=None, - ) as inp: + if os.path.isfile(outputpath): + with open(outputpath, "r", encoding=CODING, newline=None) as inp: stor = inp.readlines() for line in stor: if "FINAL SINGLE POINT ENERGY" in line: @@ -382,10 +377,7 @@ def _sp(self, silent=False): else: self.job["energy"] = 0.0 self.job["success"] = False - print( - f"{'WARNING:':{WARNLEN}}{os.path.join(self.job['workdir'], 'sp.out')} " - "doesn't exist!" - ) + print(f"{'WARNING:':{WARNLEN}}{outputpath} doesn't exist!") return def _smd_gsolv(self): @@ -409,10 +401,9 @@ def _smd_gsolv(self): keepsm = self.job["sm"] self.job["solvent"] = "gas" self.job["sm"] = "gas-phase" - self._sp(silent=True) + self._sp(silent=True, filename="sp_gas.out") if self.job["success"] == False: - self.job["success"] = False self.job["energy"] = 0.0 self.job["energy2"] = 0.0 print( @@ -423,46 +414,21 @@ def _smd_gsolv(self): else: energy_gas = self.job["energy"] self.job["energy"] = 0.0 - # mv inp inp_solv sp.out sp_solv.out - try: - shutil.move( - os.path.join(self.job["workdir"], "inp"), - os.path.join(self.job["workdir"], "inp_gas"), - ) - shutil.move( - os.path.join(self.job["workdir"], "sp.out"), - os.path.join(self.job["workdir"], "sp_gas.out"), - ) - except FileNotFoundError: - pass # calculate in solution self.job["solvent"] = keepsolv self.job["sm"] = keepsm - self._sp(silent=True) + self._sp(silent=True, filename="sp_solv.out") if self.job["success"] == False: - self.job["success"] = False self.job["energy"] = 0.0 self.job["energy2"] = 0.0 print( - f"{'ERROR:':{WARNLEN}}in gas solution phase single-point " + f"{'ERROR:':{WARNLEN}}in gas solution phase single-point " f"of {last_folders(self.job['workdir'], 2):18}" ) return else: energy_solv = self.job["energy"] self.job["energy"] = 0.0 - # mv inp inp_solv sp.out sp_solv.out - try: - shutil.move( - os.path.join(self.job["workdir"], "inp"), - os.path.join(self.job["workdir"], "inp_solv"), - ) - shutil.move( - os.path.join(self.job["workdir"], "sp.out"), - os.path.join(self.job["workdir"], "sp_solv.out"), - ) - except FileNotFoundError: - pass if self.job["success"]: if energy_solv is None or energy_gas is None: self.job["energy"] = 0.0 @@ -475,7 +441,9 @@ def _smd_gsolv(self): else: self.job["energy"] = energy_gas self.job["energy2"] = energy_solv - energy_gas - self.job["erange1"] = {self.job["temperature"] :energy_solv - energy_gas} + self.job["erange1"] = { + self.job["temperature"]: energy_solv - energy_gas + } self.job["success"] = True return @@ -502,6 +470,7 @@ def _xtbopt(self): output = "opt-part2.out" else: output = "opt-part1.out" + outputpath = os.path.join(self.job["workdir"], output) if not self.job["onlyread"]: print(f"Running optimization in {last_folders(self.job['workdir'], 2):18}") files = [ @@ -515,7 +484,6 @@ def _xtbopt(self): for file in files: if os.path.isfile(os.path.join(self.job["workdir"], file)): os.remove(os.path.join(self.job["workdir"], file)) - if not self.job["onlyread"]: # convert coord to xyz, write inp.xyz t2x(self.job["workdir"], writexyz=True, outfile="inp.xyz") # add inputfile information to coord (xtb as a driver) @@ -530,7 +498,9 @@ def _xtbopt(self): newcoord.write(line) newcoord.write("$external\n") newcoord.write(" orca input file= inp\n") - newcoord.write(f" orca bin= {os.path.join(self.job['progpath'], 'orca')}") + newcoord.write( + f" orca bin= {os.path.join(external_paths['orcapath'], 'orca')}" + ) newcoord.write("$end") with open( @@ -540,13 +510,13 @@ def _xtbopt(self): inp.write(line + "\n") # Done writing input! callargs = [ - self.job["xtb_driver_path"], + external_paths["xtbpath"], "coord", "--opt", self.job["optlevel"], "--orca", "-I", - "opt.inp" + "opt.inp", ] with open( os.path.join(self.job["workdir"], "opt.inp"), "w", newline=None @@ -565,12 +535,12 @@ def _xtbopt(self): out.write("engine=lbfgs\n") out.write("$external\n") out.write(" orca input file= inp\n") - out.write(f" orca bin= {os.path.join(self.job['progpath'], 'orca')}") + out.write( + f" orca bin= {os.path.join(external_paths['orcapath'], 'orca')}" + ) out.write("$end \n") time.sleep(0.02) - with open( - os.path.join(self.job["workdir"], output), "w", newline=None - ) as outputfile: + with open(outputpath, "w", newline=None) as outputfile: returncode = subprocess.call( callargs, shell=False, @@ -588,14 +558,10 @@ def _xtbopt(self): f"in {last_folders(self.job['workdir'], 2):18} not converged" ) time.sleep(0.02) - # check if optimization finished correctly: - if os.path.isfile(os.path.join(self.job["workdir"], output)): - with open( - os.path.join(self.job["workdir"], output), - "r", - encoding=CODING, - newline=None, - ) as inp: + # read output + # check if optimization finished correctly: + if os.path.isfile(outputpath): + with open(outputpath, "r", encoding=CODING, newline=None) as inp: stor = inp.readlines() for line in stor: if ( @@ -615,23 +581,24 @@ def _xtbopt(self): elif "*** GEOMETRY OPTIMIZATION CONVERGED AFTER " in line: self.job["cycles"] += int(line.split()[5]) self.job["converged"] = True - with open( - os.path.join(self.job["workdir"], output), - "r", - encoding=CODING, - newline=None, - ) as inp: - for line in inp: - if "av. E: " in line: - # self.job["ecyc"].append(float(line.split("Eh")[0].split()[-1])) - self.job["ecyc"].append(float(line.split("->")[-1])) + ################ + for line in stor: + if "av. E: " in line and "->" in line: + try: + self.job["ecyc"].append(float(line.split("->")[-1])) + except ValueError as e: + error_logical = True + print(f"{'ERROR:':{WARNLEN}}in CONF{self.id} calculation:\n{e}") + break if " :: gradient norm " in line: - self.job["grad_norm"] = float(line.split()[3]) + try: + self.job["grad_norm"] = float(line.split()[3]) + except ValueError as e: + error_logical = True + print(f"{'ERROR:':{WARNLEN}}in CONF{self.id} calculation:\n{e}") + break else: - print( - f"{'WARNING:':{WARNLEN}}" - f"{os.path.join(self.job['workdir'], output)} doesn't exist!" - ) + print(f"{'WARNING:':{WARNLEN}}{outputpath} doesn't exist!") error_logical = True if not error_logical: try: @@ -646,38 +613,42 @@ def _xtbopt(self): self.job["ecyc"] = [] self.job["grad_norm"] = 10.0 - # convert optimized xyz to coord file - x2t(self.job["workdir"], infile="inp.xyz") + if not self.job["onlyread"]: + # convert optimized xyz to coord file + x2t(self.job["workdir"], infile="inp.xyz") return def _nmrS(self): """ ORCA NMR shielding constant calculation """ - with open(os.path.join(self.job["workdir"], "inpS"), "w", newline=None) as inp: - for line in self._prep_input(): - inp.write(line + "\n") - # Done input! - # shielding calculation - print( - "Running shielding calculation in {:18}".format( - last_folders(self.job["workdir"], 2) - ) - ) - with open( - os.path.join(self.job["workdir"], "orcaS.out"), "w", newline=None - ) as outputfile: - call = [os.path.join(self.job["progpath"], "orca"), "inpS"] - subprocess.call( - call, - shell=False, - stdin=None, - stderr=subprocess.STDOUT, - universal_newlines=False, - cwd=self.job["workdir"], - stdout=outputfile, + if not self.job["onlyread"]: + with open( + os.path.join(self.job["workdir"], "inpS"), "w", newline=None + ) as inp: + for line in self._prep_input(): + inp.write(line + "\n") + # Done input! + # shielding calculation + print( + "Running shielding calculation in {:18}".format( + last_folders(self.job["workdir"], 2) + ) ) - time.sleep(0.1) + with open( + os.path.join(self.job["workdir"], "orcaS.out"), "w", newline=None + ) as outputfile: + call = [os.path.join(external_paths["orcapath"], "orca"), "inpS"] + subprocess.call( + call, + shell=False, + stdin=None, + stderr=subprocess.STDOUT, + universal_newlines=False, + cwd=self.job["workdir"], + stdout=outputfile, + ) + time.sleep(0.1) # check if calculation was successfull: with open( os.path.join(self.job["workdir"], "orcaS.out"), @@ -707,31 +678,34 @@ def _nmrJ(self): progpath success """ - # generate input # double hybrids not implemented - with open(os.path.join(self.job["workdir"], "inpJ"), "w", newline=None) as inp: - for line in self._prep_input(): - inp.write(line + "\n") - # Done input! - # start coupling calculation - print( - "Running coupling calculation in {}".format( - last_folders(self.job["workdir"], 2) - ) - ) - with open( - os.path.join(self.job["workdir"], "orcaJ.out"), "w", newline=None - ) as outputfile: - call = [os.path.join(self.job["progpath"], "orca"), "inpJ"] - subprocess.call( - call, - shell=False, - stdin=None, - stderr=subprocess.STDOUT, - universal_newlines=False, - cwd=self.job["workdir"], - stdout=outputfile, + if not self.job["onlyread"]: + # generate input # double hybrids not implemented + with open( + os.path.join(self.job["workdir"], "inpJ"), "w", newline=None + ) as inp: + for line in self._prep_input(): + inp.write(line + "\n") + # Done input! + # start coupling calculation + print( + "Running coupling calculation in {}".format( + last_folders(self.job["workdir"], 2) + ) ) - time.sleep(0.1) + with open( + os.path.join(self.job["workdir"], "orcaJ.out"), "w", newline=None + ) as outputfile: + call = [os.path.join(external_paths["orcapath"], "orca"), "inpJ"] + subprocess.call( + call, + shell=False, + stdin=None, + stderr=subprocess.STDOUT, + universal_newlines=False, + cwd=self.job["workdir"], + stdout=outputfile, + ) + time.sleep(0.1) # check if calculation was successfull: with open( os.path.join(self.job["workdir"], "orcaJ.out"), @@ -823,6 +797,7 @@ def _genericoutput(self): ) as out: s = sorted(zip(atom, sigma)) atom, sigma = map(list, zip(*s)) + self.shieldings = dict(zip(atom, sigma)) for i in range(len(atom)): out.write("{:{digits}} {}\n".format(atom[i], sigma[i], digits=4)) for i in range(self.job["nat"] - len(atom)): @@ -867,7 +842,7 @@ def execute(self): self._smd_gsolv() elif self.job["jobtype"] == "couplings_sp": self._nmrJ() - elif self.job["jobtype"] == "shieldings_sp": + elif self.job["jobtype"] in ("shieldings_sp", "shieldings"): self._nmrS() elif self.job["jobtype"] == "genericout": self._genericoutput() diff --git a/censo_qm/parallel.py b/censo_qm/parallel.py index 6e56aa0..f66b367 100644 --- a/censo_qm/parallel.py +++ b/censo_qm/parallel.py @@ -11,21 +11,22 @@ from .utilities import print from .cfg import ENVIRON -def balance_load(P,O,nconf, do_change): + +def balance_load(P, O, nconf, do_change): """Balance calculation load between threads (P) and number of cores per thread (O) """ changed = False - max_cores = P*O + max_cores = P * O if do_change: if nconf < P: try: P = nconf - O=1 + O = 1 while True: - if P*O <= max_cores: - if P*(O+1) <=max_cores: - O+=1 + if P * O <= max_cores: + if P * (O + 1) <= max_cores: + O += 1 else: break else: @@ -34,9 +35,12 @@ def balance_load(P,O,nconf, do_change): except Exception: pass if changed: - print(f"Adjusting the number of threads (P) = {P} and number of cores per thread (O) = {O}") + print( + f"Adjusting the number of threads (P) = {P} and number of cores per thread (O) = {O}" + ) return P, O, changed + def execute_data(q, resultq): """ code that the worker has to execute @@ -66,7 +70,16 @@ def execute_data(q, resultq): def run_in_parallel( - config, q, resultq, job, maxthreads, omp, loopover, instructdict, balance=False,foldername="" + config, + q, + resultq, + job, + maxthreads, + omp, + loopover, + instructdict, + balance=False, + foldername="", ): """Run jobs in parallel q = queue to put assemble tasks @@ -105,7 +118,7 @@ def run_in_parallel( njobs = q.qsize() time.sleep(0.02) if changed: - ENVIRON['PARNODES'] = str(omp) + ENVIRON["PARNODES"] = str(omp) if instructdict.get("onlyread", False): print(f"\nReading data from {njobs} conformers calculated in previous run.") @@ -166,5 +179,5 @@ def run_in_parallel( print(f"ERROR some conformers were lost!") if changed: - ENVIRON['PARNODES'] = str(omp_initial) + ENVIRON["PARNODES"] = str(omp_initial) return results diff --git a/censo_qm/prescreening.py b/censo_qm/prescreening.py index 14649dc..860f0fd 100755 --- a/censo_qm/prescreening.py +++ b/censo_qm/prescreening.py @@ -7,7 +7,7 @@ import math import time from multiprocessing import JoinableQueue as Queue -from .cfg import PLENGTH, DIGILEN, AU2KCAL, CODING, WARNLEN +from .cfg import PLENGTH, DIGILEN, AU2KCAL, CODING, WARNLEN, qm_prepinfo from .parallel import run_in_parallel from .orca_job import OrcaJob from .tm_job import TmJob @@ -26,6 +26,7 @@ print, print_errors, calc_boltzmannweights, + conf_in_interval, ) @@ -194,6 +195,7 @@ def part1(config, conformers, store_confs, ensembledata): "energy": 0.0, "energy2": 0.0, "success": False, + "onlyread": config.onlyread, } if config.solvent == "gas": @@ -203,8 +205,6 @@ def part1(config, conformers, store_confs, ensembledata): ) name = "prescreening_single-point" folder = instruction["func"] - if config.prog == "orca": - instruction["progpath"] = config.external_paths["orcapath"] else: if config.smgsolv1 in config.smgsolv_1: # additive Gsolv @@ -232,9 +232,6 @@ def part1(config, conformers, store_confs, ensembledata): elif instruction["sm"] in ("gbsa_gsolv", "alpb_gsolv"): # do DFT gas phase sp and additive Gsolv instruction["jobtype"] = instruction["sm"] - if config.prog == "orca": - instruction["progpath"] = config.external_paths["orcapath"] - instruction["xtb_driver_path"] = config.external_paths["xtbpath"] instruction["method"], instruction["method2"] = config.get_method_name( instruction["jobtype"], func=instruction["func"], @@ -267,9 +264,6 @@ def part1(config, conformers, store_confs, ensembledata): else: # with implicit solvation instruction["jobtype"] = "sp_implicit" - # instruction["prepinfo"] = ["low+"] - if config.prog == "orca": - instruction["progpath"] = config.external_paths["orcapath"] instruction["method"], instruction["method2"] = config.get_method_name( "sp_implicit", func=instruction["func"], @@ -279,6 +273,18 @@ def part1(config, conformers, store_confs, ensembledata): name = "prescreening_single-point" folder = instruction["func"] + # Supporting info update: + try: + if job == TmJob: + ensembledata.si["part1"]["Energy_settings"] = " ".join( + qm_prepinfo["tm"][instruction["prepinfo"][0]] + ).replace("-", "") + elif job == OrcaJob: + ensembledata.si["part1"]["Energy_settings"] = " ".join( + qm_prepinfo["orca"][instruction["prepinfo"][0]] + ) + except TypeError as e: + print(e) check = {True: "was successful", False: "FAILED"} start_firstsort = time.perf_counter() if calculate: @@ -302,20 +308,69 @@ def part1(config, conformers, store_confs, ensembledata): calculate, store_confs, save_errors = ensemble2coord( config, folder, calculate, store_confs, save_errors ) - - # parallel calculation: - calculate = run_in_parallel( - config, - q, - resultq, - job, - config.maxthreads, - config.omp, - calculate, - instruction, - config.balance, - folder, - ) + tmppath = os.path.join(config.cwd, "part1preG.dat") + if config.onlyread and not os.path.isfile(tmppath): + print( + f"Reading data from {tmppath} is not possible since it can not be found." + ) + print( + f"{'ERROR:':{WARNLEN}}The re-read data in part1 can be the data from part2" + " since files have been overwritten!" + ) + if config.onlyread and os.path.isfile(tmppath): + #################################################################### + # this is only a temporary fix, and needed because in part1 + # COSMO-RS (and single-point) data is overwritten in part2. + # This has to be changed when backward compability is broken. + with open(tmppath, "r", newline=None, encoding=CODING) as inp: + data = inp.readlines() + if instruction["jobtype"] in ("sp", "sp_implicit"): + for line in data[3:]: + confid = int(line.strip().split()[0][4:]) + try: + energy = float(line.strip().split()[3]) + except (ValueError, IndexError) as e: + print(e) + for conf in calculate: + if conf.id == confid: + conf.job["energy"] = energy + conf.job["success"] = True + conf.job["workdir"] = f"CONF{confid}/{folder}" + elif instruction["jobtype"] in ( + "cosmors", + "smd_gsolv", + "gbsa_gsolv", + "alpb_gsolv", + ): + for line in data[3:]: + confid = int(line.strip().split()[0][4:]) + try: + energy = float(line.strip().split()[3]) + gsolv = float(line.strip().split()[4]) + except (ValueError, IndexError) as e: + print(e) + for conf in calculate: + if conf.id == confid: + conf.job["energy"] = energy + conf.job["energy2"] = gsolv + conf.job["erange1"] = {instruction["temperature"]: gsolv} + conf.job["success"] = True + conf.job["workdir"] = f"CONF{confid}/{folder}" + #################################################################### + else: + # parallel calculation: + calculate = run_in_parallel( + config, + q, + resultq, + job, + config.maxthreads, + config.omp, + calculate, + instruction, + config.balance, + folder, + ) for conf in list(calculate): if instruction["jobtype"] in ("sp", "sp_implicit"): @@ -336,7 +391,9 @@ def part1(config, conformers, store_confs, ensembledata): conf.prescreening_sp_info["method"] = conf.job["method"] if instruction["jobtype"] == "sp_implicit": conf.prescreening_gsolv_info["energy"] = 0.0 - conf.prescreening_gsolv_info["range"] = {conf.job['temperature']: 0.0,} + conf.prescreening_gsolv_info["range"] = { + conf.job["temperature"]: 0.0 + } conf.prescreening_gsolv_info["info"] = "calculated" conf.prescreening_gsolv_info["method"] = conf.job["method2"] elif instruction["jobtype"] in ( @@ -408,7 +465,28 @@ def part1(config, conformers, store_confs, ensembledata): ) calculate.append(prev_calculated.pop(prev_calculated.index(conf))) end_firstsort = time.perf_counter() - ensembledata.part_info["part1_firstsort"] = end_firstsort -start_firstsort + ensembledata.part_info["part1_firstsort"] = end_firstsort - start_firstsort + + # create data for possible SI generation:----------------------------------- + # Energy: + ensembledata.si["part1"]["Energy"] = instruction["method"] + if "DOGCP" in instruction["prepinfo"]: + ensembledata.si["part1"]["Energy"] += " + GCP" + # Energy_settings are set above! + # GmRRHO is done below! + # Solvation: + if config.solvent == "gas": + ensembledata.si["part1"]["G_solv"] = "gas-phase" + else: + ensembledata.si["part1"]["G_solv"] = instruction["method2"] + # Geometry: + ensembledata.si["part1"]["Geometry"] = "GFNn-xTB (input geometry)" + # QM-CODE: + ensembledata.si["part1"]["main QM code"] = str(config.prog).upper() + # Threshold: + ensembledata.si["part1"]["Threshold"] = f"{config.part1_threshold} kcal/mol" + # END SI generation -------------------------------------------------------- + for conf in calculate: conf.reset_job_info() if not calculate: @@ -431,12 +509,12 @@ def part1(config, conformers, store_confs, ensembledata): solv = "prescreening_gsolv_info" e = "prescreening_sp_info" conf.calc_free_energy( - e=e, + e=e, solv=solv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) try: minfree = min([i.free_energy for i in calculate if i is not None]) except ValueError: @@ -545,6 +623,9 @@ def part1(config, conformers, store_confs, ensembledata): conf = calculate.pop(calculate.index(conf)) conf.__class__ = job conf.job["success"] = True + conf.sym = conf.prescreening_grrho_info.get("sym", conf.sym) + conf.linear = conf.prescreening_grrho_info.get("linear", conf.linear) + conf.symnum = conf.prescreening_grrho_info.get("symnum", conf.symnum) prev_calculated.append(conf) if not calculate and not prev_calculated: @@ -565,7 +646,6 @@ def part1(config, conformers, store_confs, ensembledata): "charge": config.charge, "unpaired": config.unpaired, "solvent": config.solvent, - "progpath": config.external_paths["xtbpath"], "bhess": config.bhess, "consider_sym": config.consider_sym, "sm_rrho": config.sm_rrho, @@ -580,6 +660,7 @@ def part1(config, conformers, store_confs, ensembledata): "imagthr": config.imagthr, "sthr": config.sthr, "scale": config.scale, + "onlyread": config.onlyread, } instruction_prerrho["method"], _ = config.get_method_name( @@ -589,6 +670,14 @@ def part1(config, conformers, store_confs, ensembledata): sm=instruction_prerrho["sm_rrho"], solvent=instruction_prerrho["solvent"], ) + # GmRRHO for SI information: + if config.evaluate_rrho: + ensembledata.si["part1"]["G_mRRHO"] = instruction_prerrho["method"] + if "bhess" in ensembledata.si["part1"]["G_mRRHO"]: + ensembledata.si["part1"]["G_mRRHO"] += " SPH" + else: + ensembledata.si["part1"]["G_mRRHO"] = "not included" + if calculate: print("The prescreening G_mRRHO calculation is now performed for:") print_block(["CONF" + str(i.id) for i in calculate]) @@ -640,9 +729,14 @@ def part1(config, conformers, store_confs, ensembledata): conf.sym = conf.job["symmetry"] conf.symnum = conf.job["symnum"] conf.linear = conf.job["linear"] + conf.prescreening_grrho_info["sym"] = conf.job["symmetry"] + conf.prescreening_grrho_info["symnum"] = conf.job["symnum"] + conf.prescreening_grrho_info["linear"] = conf.job["linear"] conf.prescreening_grrho_info["rmsd"] = conf.job["rmsd"] conf.prescreening_grrho_info["energy"] = conf.job["energy"] - conf.prescreening_grrho_info["range"] = {conf.job['temperature']: conf.job["energy"],} + conf.prescreening_grrho_info["range"] = { + conf.job["temperature"]: conf.job["energy"] + } conf.prescreening_grrho_info["info"] = "calculated" conf.prescreening_grrho_info["method"] = instruction_prerrho[ "method" @@ -766,19 +860,19 @@ def part1(config, conformers, store_confs, ensembledata): else: solv = "prescreening_gsolv_info" conf.calc_free_energy( - e="prescreening_sp_info", - solv=solv, + e="prescreening_sp_info", + solv=solv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) conf.xtb_free_energy = conf.calc_free_energy( e="xtb_energy", solv=None, - rrho=rrho, + rrho=rrho, out=True, t=config.temperature, - consider_sym=config.consider_sym + consider_sym=config.consider_sym, ) try: @@ -808,6 +902,9 @@ def part1(config, conformers, store_confs, ensembledata): columndescription2=columndescription2, ) # -------------------------------------------------------------------------- + calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) + conf_in_interval(calculate) + # -------------------------------------------------------------------------- for conf in calculate: if conf.free_energy == minfree: ensembledata.bestconf["part1"] = conf.id @@ -834,14 +931,16 @@ def part1(config, conformers, store_confs, ensembledata): conf.get_mrrho( config.temperature, rrho="prescreening_grrho_info", - consider_sym=config.consider_sym - ) * AU2KCAL + consider_sym=config.consider_sym, + ) + * AU2KCAL for conf in calculate if conf.get_mrrho( config.temperature, rrho="prescreening_grrho_info", - consider_sym=config.consider_sym - ) is not None + consider_sym=config.consider_sym, + ) + is not None ] ) max_fuzzy = 1 @@ -868,12 +967,12 @@ def part1(config, conformers, store_confs, ensembledata): solv = "prescreening_gsolv_info" e = "prescreening_sp_info" conf.calc_free_energy( - e=e, - solv=solv, - rrho=rrho, + e=e, + solv=solv, + rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) try: minfree = min([i.free_energy for i in calculate if i is not None]) except ValueError: @@ -892,12 +991,12 @@ def part1(config, conformers, store_confs, ensembledata): solv = "prescreening_gsolv_info" e = "prescreening_sp_info" conf.calc_free_energy( - e=e, - solv=solv, + e=e, + solv=solv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) try: minfree = min([i.free_energy for i in calculate if i is not None]) except ValueError: @@ -976,12 +1075,12 @@ def part1(config, conformers, store_confs, ensembledata): else: solv = "prescreening_gsolv_info" conf.calc_free_energy( - e="prescreening_sp_info", - solv=solv, + e="prescreening_sp_info", + solv=solv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) # write coord.enso_best for conf in calculate: @@ -1001,9 +1100,9 @@ def part1(config, conformers, store_confs, ensembledata): "$coord # {} {} !CONF{} \n".format( conf.free_energy, conf.get_mrrho( - config.temperature, - rrho="prescreening_grrho_info", - consider_sym=config.consider_sym + config.temperature, + rrho="prescreening_grrho_info", + consider_sym=config.consider_sym, ), conf.id, ) @@ -1057,12 +1156,12 @@ def part1(config, conformers, store_confs, ensembledata): solv = "prescreening_gsolv_info" e = "prescreening_sp_info" conf.calc_free_energy( - e=e, - solv=solv, + e=e, + solv=solv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) avG = 0.0 @@ -1072,12 +1171,10 @@ def part1(config, conformers, store_confs, ensembledata): for conf in calculate: avG += conf.bm_weight * conf.free_energy avE += conf.bm_weight * conf.prescreening_sp_info["energy"] - #avGRRHO += conf.bm_weight * conf.prescreening_grrho_info["energy"] + # avGRRHO += conf.bm_weight * conf.prescreening_grrho_info["energy"] avGRRHO += conf.bm_weight * conf.get_mrrho( - config.temperature, - "prescreening_grrho_info", - config.consider_sym - ) + config.temperature, "prescreening_grrho_info", config.consider_sym + ) avGsolv += conf.bm_weight * conf.prescreening_gsolv_info["energy"] # printout: @@ -1121,7 +1218,6 @@ def part1(config, conformers, store_confs, ensembledata): "charge": config.charge, "unpaired": config.unpaired, "solvent": config.solvent, - "progpath": config.external_paths["xtbpath"], "sm": config.sm_rrho, "rmsdbias": config.rmsdbias, "temperature": config.temperature, @@ -1129,6 +1225,7 @@ def part1(config, conformers, store_confs, ensembledata): "energy": 0.0, "energy2": 0.0, "success": False, + "onlyread": config.onlyread, } folder_gfn = "GFN_unbiased" save_errors, store_confs, calculate = new_folders( diff --git a/censo_qm/qm_job.py b/censo_qm/qm_job.py index 994e6c8..bcb0379 100644 --- a/censo_qm/qm_job.py +++ b/censo_qm/qm_job.py @@ -4,6 +4,7 @@ code. (xTB always available) """ import os + try: from math import isclose except ImportError: @@ -11,7 +12,7 @@ import time import subprocess import json -from .cfg import ENVIRON, CODING, WARNLEN, censo_solvent_db, rot_sym_num +from .cfg import ENVIRON, CODING, WARNLEN, censo_solvent_db, rot_sym_num, external_paths from .utilities import last_folders, print from .datastructure import MoleculeData @@ -36,6 +37,7 @@ def reset_job_info(self): "method2": "", # description of the method "workdir": "", "copymos": "", + "moread": None, "omp": 1, "charge": 0, "unpaired": 0, @@ -87,7 +89,7 @@ def reset_job_info(self): "internal_error": [], # "cosmorsparam": "", # normal/fine - "symnum":1, + "symnum": 1, } def _sp(self, silent=False): @@ -113,7 +115,7 @@ def _genericoutput(self): def _get_sym_num(self, sym=None, linear=None): """Get rotational symmetry number from Schoenfließ symbol""" if sym is None: - sym = 'c1' + sym = "c1" if linear is None: linear = False symnum = 1 @@ -129,77 +131,77 @@ def _get_sym_num(self, sym=None, linear=None): break return symnum - def _xtb_sp(self): + def _xtb_sp(self, filename="sp.out", silent=False): """ - Get single-point energy from xTB + Get single-point energy from GFNn-xTB + --> return self.job["energy"] + --> return self.job["success"] """ - files = [ - "xtbrestart", - "xtbtopo.mol", - "xcontrol-inp", - "wbo", - "charges", - "gfnff_topo", - "sp.out", - ] - for file in files: - if os.path.isfile(os.path.join(self.job["workdir"], file)): - os.remove(os.path.join(self.job["workdir"], file)) - # run single-point: - - call = [ - self.job["progpath"], - "coord", - "--" + self.job["gfn_version"], - "--sp", - "--chrg", - str(self.job["charge"]), - "--norestart", - ] - if self.job["solvent"] != "gas": - call.extend( - [ - "--" + str(self.job["sm"]), - censo_solvent_db[self.job["solvent"]]["xtb"][1], - "-I", - "xcontrol-inp", - ] - ) - with open( + outputpath = os.path.join(self.job["workdir"], filename) + if not self.job["onlyread"]: + files = [ + "xtbrestart", + "xtbtopo.mol", + "xcontrol-inp", + "wbo", + "charges", + "gfnff_topo", + filename, + ] + for file in files: + if os.path.isfile(os.path.join(self.job["workdir"], file)): + os.remove(os.path.join(self.job["workdir"], file)) + # run single-point: + call = [ + external_paths["xtbpath"], + "coord", + "--" + self.job["gfn_version"], + "--sp", + "--chrg", + str(self.job["charge"]), + "--norestart", + "--parallel", + str(self.job["omp"]), + ] + if self.job["solvent"] != "gas": + call.extend( + [ + "--" + str(self.job["sm"]), + censo_solvent_db[self.job["solvent"]]["xtb"][1], + "reference", + "-I", + "xcontrol-inp", + ] + ) + with open( os.path.join(self.job["workdir"], "xcontrol-inp"), "w", newline=None ) as xcout: - xcout.write("$gbsa\n") - xcout.write(" gbsagrid=tight\n") - xcout.write("$end\n") - with open( - os.path.join(self.job["workdir"], "sp.out"), "w", newline=None - ) as outputfile: - returncode = subprocess.call( - call, - shell=False, - stdin=None, - stderr=subprocess.STDOUT, - universal_newlines=False, - cwd=self.job["workdir"], - stdout=outputfile, - env=ENVIRON, - ) - if returncode != 0: - self.job["energy"] = 0.0 - self.job["success"] = False - print( - f"{'ERROR:':{WARNLEN}}{self.job['gfn_version']}-xTB error in " - f"{last_folders(self.job['workdir'], 2)}" - ) - return - # read gas phase energy: - if os.path.isfile(os.path.join(self.job["workdir"], "sp.out")): - with open( - os.path.join(self.job["workdir"], "sp.out"), - "r", - encoding=CODING, - newline=None, - ) as inp: + xcout.write("$gbsa\n") + xcout.write(" gbsagrid=tight\n") + xcout.write("$end\n") + with open(outputpath, "w", newline=None) as outputfile: + returncode = subprocess.call( + call, + shell=False, + stdin=None, + stderr=subprocess.STDOUT, + universal_newlines=False, + cwd=self.job["workdir"], + stdout=outputfile, + env=ENVIRON, + ) + if returncode != 0: + self.job["energy"] = 0.0 + self.job["success"] = False + print( + f"{'ERROR:':{WARNLEN}}{self.job['gfn_version']}-xTB error in " + f"{last_folders(self.job['workdir'], 2)}" + ) + return + time.sleep(0.02) + # read energy: + if os.path.isfile(outputpath): + with open(outputpath, "r", encoding=CODING, newline=None) as inp: store = inp.readlines() for line in store: if "| TOTAL ENERGY" in line: @@ -207,25 +209,28 @@ def _xtb_sp(self): self.job["energy"] = float(line.split()[3]) self.job["success"] = True except Exception: - print( - f"{'ERROR:':{WARNLEN}}while converting " - f"single-point in: {last_folders(self.job['workdir'], 2)}" - ) + if not silent: + print( + f"{'ERROR:':{WARNLEN}}while converting " + f"single-point in: {last_folders(self.job['workdir'], 2)}" + ) self.job["energy"] = 0.0 self.job["success"] = False return else: self.job["energy"] = 0.0 self.job["success"] = False + if not silent: + print( + f"{'ERROR:':{WARNLEN}}{self.job['gfn_version']}-xTB error in " + f"{last_folders(self.job['workdir'], 2)}" + ) + return + if not silent: print( - f"{'ERROR:':{WARNLEN}}{self.job['gfn_version']}-xTB error in " - f"{last_folders(self.job['workdir'], 2)}" + f"{self.job['gfn_version']}-xTB energy for {last_folders(self.job['workdir'], 2)}" + f" = {self.job['energy']: >.7f}" ) - return - print( - f"{self.job['gfn_version']}-xTB energy for {last_folders(self.job['workdir'], 2)}" - f" = {self.job['energy']: >.7f}" - ) def _xtb_gsolv(self): """ @@ -234,82 +239,28 @@ def _xtb_gsolv(self): using xTB and the GFNn or GFN-FF hamiltonian. --> return gsolv at energy2 """ - tmp_gas = 0 - tmp_solv = 0 - if self.job["jobtype"] == "gbsa_gsolv": - xtbsm = "gbsa" - elif self.job["jobtype"] == "alpb_gsolv": - xtbsm = "alpb" print( f"Running {self.job['jobtype'].upper()} calculation in " f"{last_folders(self.job['workdir'], 3)}" ) - files = [ - "xtbrestart", - "xtbtopo.mol", - "xcontrol-inp", - "wbo", - "charges", - "gfnff_topo", - "gas.out", - "solv.out", - ] - for file in files: - if os.path.isfile(os.path.join(self.job["workdir"], file)): - os.remove(os.path.join(self.job["workdir"], file)) - # run gas phase single-point: - with open( - os.path.join(self.job["workdir"], "gas.out"), "w", newline=None - ) as outputfile: - returncode = subprocess.call( - [ - self.job["xtb_driver_path"], - "coord", - "--" + self.job["gfn_version"], - "--sp", - "--chrg", - str(self.job["charge"]), - "--norestart", - ], - shell=False, - stdin=None, - stderr=subprocess.STDOUT, - universal_newlines=False, - cwd=self.job["workdir"], - stdout=outputfile, - env=ENVIRON, - ) - if returncode != 0: - self.job["energy2"] = 0.0 + tmp_gas = None + tmp_solv = None + # keep possible DFT energy: + keep_dft_energy = self.job["energy"] + self.job["energy"] = 0.0 + keep_solvent = self.job["solvent"] + self.job["solvent"] = "gas" + self.job["sm"] = "gas-phase" + + # run gas-phase GFNn-xTB single-point + self._xtb_sp(filename="gas.out", silent=True) + if self.job["success"]: + tmp_gas = self.job["energy"] + # reset for next calculation + self.job["energy"] = 0.0 self.job["success"] = False - print( - f"{'ERROR:':{WARNLEN}}Gas phase {self.job['gfn_version']}-xTB error in " - f"{last_folders(self.job['workdir'], 3)}" - ) - return - # read gas phase energy: - if os.path.isfile(os.path.join(self.job["workdir"], "gas.out")): - with open( - os.path.join(self.job["workdir"], "gas.out"), - "r", - encoding=CODING, - newline=None, - ) as inp: - store = inp.readlines() - for line in store: - if "| TOTAL ENERGY" in line: - try: - tmp_gas = float(line.split()[3]) - self.job["success"] = True - except Exception: - print( - f"{'ERROR:':{WARNLEN}}while converting gas phase " - f"single-point in: {last_folders(self.job['workdir'], 3)}" - ) - self.job["energy2"] = 0.0 - self.job["success"] = False - return else: + self.job["energy"] = keep_dft_energy self.job["energy2"] = 0.0 self.job["success"] = False print( @@ -317,73 +268,22 @@ def _xtb_gsolv(self): f"{last_folders(self.job['workdir'], 3)}" ) return + # run gas-phase GFNn-xTB single-point # run single-point in solution: - # ``reference'' corresponds to 1\;bar of ideal gas and 1\;mol/L of liquid + # ''reference'' corresponds to 1\;bar of ideal gas and 1\;mol/L of liquid # solution at infinite dilution, - with open( - os.path.join(self.job["workdir"], "xcontrol-inp"), "w", newline=None - ) as xcout: - xcout.write("$gbsa\n") - xcout.write(" gbsagrid=tight\n") - xcout.write("$end\n") - with open( - os.path.join(self.job["workdir"], "solv.out"), "w", newline=None - ) as outputfile: - returncode = subprocess.call( - [ - self.job["xtb_driver_path"], - "coord", - "--" + self.job["gfn_version"], - "--sp", - "--" + xtbsm, - censo_solvent_db[self.job["solvent"]]["xtb"][1], - "reference", - "--chrg", - str(self.job["charge"]), - "--norestart", - "-I", - "xcontrol-inp", - ], - shell=False, - stdin=None, - stderr=subprocess.STDOUT, - universal_newlines=False, - cwd=self.job["workdir"], - stdout=outputfile, - env=ENVIRON, - ) - if returncode != 0: - self.job["energy2"] = 0.0 - self.job["success"] = False - print( - f"{'ERROR:':{WARNLEN}}Solution phase {self.job['gfn_version']}-xTB error in " - f"{last_folders(self.job['workdir'], 3)}" - ) - return - time.sleep(0.05) - # #read solv.out - if os.path.isfile(os.path.join(self.job["workdir"], "solv.out")): - with open( - os.path.join(self.job["workdir"], "solv.out"), - "r", - encoding=CODING, - newline=None, - ) as inp: - store = inp.readlines() - for line in store: - if "| TOTAL ENERGY" in line: - try: - tmp_solv = float(line.split()[3]) - self.job["success"] = True - except Exception: - print( - f"{'ERROR:':{WARNLEN}}while converting solution phase " - f"single-point in: {last_folders(self.job['workdir'], 3)}" - ) - self.job["energy2"] = 0.0 - self.job["success"] = False - return + if self.job["jobtype"] == "gbsa_gsolv": + self.job["sm"] = "gbsa" + elif self.job["jobtype"] == "alpb_gsolv": + self.job["sm"] = "alpb" + self.job["solvent"] = keep_solvent + self._xtb_sp(filename="solv.out", silent=True) + if self.job["success"]: + tmp_solv = self.job["energy"] + # reset + self.job["energy"] = keep_dft_energy else: + self.job["energy"] = keep_dft_energy self.job["energy2"] = 0.0 self.job["success"] = False print( @@ -397,15 +297,16 @@ def _xtb_gsolv(self): self.job["success"] = False else: self.job["energy2"] = tmp_solv - tmp_gas - self.job["erange1"] = {self.job["temperature"] :tmp_solv - tmp_gas} - self.job["success"] = True + self.job["erange1"] = {self.job["temperature"]: tmp_solv - tmp_gas} self.job["energy_xtb_gas"] = tmp_gas self.job["energy_xtb_solv"] = tmp_solv - def _xtbrrho(self): + def _xtbrrho(self, filename="ohess.out"): """ - mRRHO contribution with GFNn/GFN-FF-XTB + mRRHO contribution with GFNn/GFN-FF-xTB """ + outputpath = os.path.join(self.job["workdir"], filename) + xcontrolpath = os.path.join(self.job["workdir"], "xcontrol-inp") if not self.job["onlyread"]: print( f"Running {str(self.job['gfn_version']).upper()}-xTB mRRHO in " @@ -427,9 +328,7 @@ def _xtbrrho(self): if isclose(self.job["temperature"], t, abs_tol=0.6): self.job["trange"].pop(self.job["trange"].index(t)) self.job["trange"].append(self.job["temperature"]) - with open( - os.path.join(self.job["workdir"], "xcontrol-inp"), "w", newline=None - ) as xcout: + with open(xcontrolpath, "w", newline=None) as xcout: xcout.write("$thermo\n") if self.job["trange"]: xcout.write( @@ -438,33 +337,29 @@ def _xtbrrho(self): ) else: xcout.write(" temp={}\n".format(self.job["temperature"])) - if self.job.get("sthr", "automatic") == 'automatic': + if self.job.get("sthr", "automatic") == "automatic": xcout.write(" sthr=50.0\n") else: - xcout.write(" sthr={}\n".format(self.job['sthr'])) - if self.job.get("imagthr", "automatic") == 'automatic': + xcout.write(" sthr={}\n".format(self.job["sthr"])) + if self.job.get("imagthr", "automatic") == "automatic": if self.job["bhess"]: xcout.write(" imagthr={}\n".format("-100")) else: xcout.write(" imagthr={}\n".format("-50")) else: xcout.write(" imagthr={}\n".format(self.job["imagthr"])) - if self.job.get("scale", "automatic") == 'automatic': - #xcout.write(" scale={}\n".format("1.0")) - pass # is method dependant leave it to xTB e.g. GFNFF has a + if self.job.get("scale", "automatic") != "automatic": + # for automatic --> is method dependant leave it to xTB e.g. GFNFF has a # different scaling factor than GFN2 - else: - xcout.write(" scale={}\n".format(self.job['scale'])) + xcout.write(" scale={}\n".format(self.job["scale"])) xcout.write("$symmetry\n") - if self.job["consider_sym"]: - # xcout.write(" desy=0.1\n") # taken from xtb defaults - xcout.write(" maxat=1000\n") - # always consider symmetry! - #else: - # xcout.write(" desy=0.0\n") + xcout.write(" maxat=1000\n") + # always consider symmetry + # xcout.write(" desy=0.1\n") # taken from xtb defaults + # xcout.write(" desy=0.0\n") if self.job["solvent"] != "gas": xcout.write("$gbsa\n") - xcout.write(" gbsagrid=tight\n") + xcout.write(" gbsagrid=tight\n") xcout.write("$end\n") if self.job["bhess"]: # set ohess or bhess @@ -473,13 +368,11 @@ def _xtbrrho(self): else: dohess = "--ohess" olevel = "vtight" - time.sleep(0.05) - with open( - os.path.join(self.job["workdir"], "ohess.out"), "w", newline=None - ) as outputfile: + time.sleep(0.02) + with open(outputpath, "w", newline=None) as outputfile: if self.job["solvent"] != "gas": callargs = [ - self.job["progpath"], + external_paths["xtbpath"], "coord", "--" + str(self.job["gfn_version"]), dohess, @@ -492,10 +385,12 @@ def _xtbrrho(self): "--norestart", "-I", "xcontrol-inp", + "--parallel", + str(self.job["omp"]), ] else: callargs = [ - self.job["progpath"], + external_paths["xtbpath"], "coord", "--" + str(self.job["gfn_version"]), dohess, @@ -506,6 +401,8 @@ def _xtbrrho(self): "--norestart", "-I", "xcontrol-inp", + "--parallel", + str(self.job["omp"]), ] if self.job["rmsdbias"]: callargs.extend( @@ -524,7 +421,6 @@ def _xtbrrho(self): stdout=outputfile, env=ENVIRON, ) - time.sleep(0.05) # check if converged: if returncode != 0: self.job["energy"] = 0.0 @@ -536,7 +432,7 @@ def _xtbrrho(self): print(self.job["errormessage"][-1]) return # start reading output! - if not os.path.isfile(os.path.join(self.job["workdir"], "ohess.out")): + if not os.path.isfile(outputpath): self.job["energy"] = 0.0 self.job["success"] = False self.job["errormessage"].append( @@ -545,14 +441,8 @@ def _xtbrrho(self): print(self.job["errormessage"][-1]) return # start reading file: - with open( - os.path.join(self.job["workdir"], "ohess.out"), - "r", - encoding=CODING, - newline=None, - ) as inp: - store = inp.readlines() - + with open(outputpath, "r", encoding=CODING, newline=None) as inp: + store = inp.readlines() if self.job["trange"]: gt = {} ht = {} @@ -587,22 +477,21 @@ def _xtbrrho(self): else: self.job["success"] = False return - # end self.trange + # end self.job["trange"] for line in store: - if "final rmsd / " in line and self.job['bhess']: + if "final rmsd / " in line and self.job["bhess"]: try: self.job["rmsd"] = float(line.split()[3]) - except (ValueError): + except (ValueError, IndexError): self.job["rmsd"] = 0.0 if ": linear? " in line: - # linear needed for symmetry and S_rot (only if turned off) + # linear needed for symmetry and S_rot (only if considersym is turned off) try: if line.split()[2] == "false": self.job["linear"] = False elif line.split()[2] == "true": self.job["linear"] = True - except (IndexError,Exception) as e: - #pass + except (IndexError, Exception) as e: print(e) if os.path.isfile(os.path.join(self.job["workdir"], "xtb_enso.json")): with open( @@ -641,7 +530,9 @@ def _xtbrrho(self): self.job["energy"] = 0.0 self.job["success"] = False # calculate symnum - self.job["symnum"] = self._get_sym_num(sym=self.job['symmetry'], linear=self.job["linear"]) + self.job["symnum"] = self._get_sym_num( + sym=self.job["symmetry"], linear=self.job["linear"] + ) else: print( f"{'WARNING:':{WARNLEN}}File " @@ -650,7 +541,6 @@ def _xtbrrho(self): self.job["energy"] = 0.0 self.job["success"] = False - def execute(self): """ Chooses which function to execute based on jobtype. diff --git a/censo_qm/refinement.py b/censo_qm/refinement.py index d1de651..2388f81 100755 --- a/censo_qm/refinement.py +++ b/censo_qm/refinement.py @@ -7,7 +7,7 @@ import shutil import os import sys -from .cfg import PLENGTH, CODING, AU2KCAL, DIGILEN, WARNLEN +from .cfg import PLENGTH, CODING, AU2KCAL, DIGILEN, WARNLEN, qm_prepinfo from .utilities import ( check_for_folder, print_block, @@ -22,6 +22,7 @@ print, print_errors, ensemble2coord, + conf_in_interval, ) from .parallel import run_in_parallel from .orca_job import OrcaJob @@ -77,19 +78,19 @@ def part3(config, conformers, store_confs, ensembledata): info.append( ["bhess", "Apply constraint to input geometry during mRRHO calculation"] ) - + max_len_digilen = 0 for item in info: - if item[0] == 'justprint': + if item[0] == "justprint": if "short-notation" in item[1]: - tmp = len(item[1]) -len('short-notation:') + tmp = len(item[1]) - len("short-notation:") else: tmp = len(item[1]) else: tmp = len(item[1]) if tmp > max_len_digilen: max_len_digilen = tmp - max_len_digilen +=1 + max_len_digilen += 1 if max_len_digilen < DIGILEN: max_len_digilen = DIGILEN @@ -191,21 +192,43 @@ def part3(config, conformers, store_confs, ensembledata): geometries_from_input = False # check if calculated on unoptimized geometries: - if any([conf.optimization_info["info"] == "not_calculated" for conf in calculate+prev_calculated]): + if any( + [ + conf.optimization_info["info"] == "not_calculated" + for conf in calculate + prev_calculated + ] + ): if config.part2: - print_errors(f"{'ERROR:':{WARNLEN}}Calculating (free) energies on " - f"DFT unoptimized geometries!\n" - f"{'':{WARNLEN}}Even though part2 is calculated!\n" - f"{'':{WARNLEN}}Calculation on mixture of optimized " - f"and unoptimized geometries is not advised!", - save_errors) + print_errors( + f"{'ERROR:':{WARNLEN}}Calculating (free) energies on " + f"DFT unoptimized geometries!\n" + f"{'':{WARNLEN}}Even though part2 is calculated!\n" + f"{'':{WARNLEN}}Calculation on mixture of optimized " + f"and unoptimized geometries is not advised!", + save_errors, + ) print("Going to exit!") sys.exit(1) else: - print_errors(f"{'WARNING:':{WARNLEN}} Calculating (free) energies " - f"on DFT unoptimized geometries!", save_errors) + print_errors( + f"{'WARNING:':{WARNLEN}} Calculating (free) energies " + f"on DFT unoptimized geometries!", + save_errors, + ) geometries_from_input = True - + # SI geometry: + if geometries_from_input: + ensembledata.si["part3"]["Geometry"] = "GFNn-xTB (input geometry)" + else: + ensembledata.si["part3"]["Geometry"], _ = config.get_method_name( + "xtbopt", + func=config.func, + basis=config.basis, + solvent=config.solvent, + sm=config.sm2, + ) + ensembledata.si["part3"]["Geometry"] += f" @optlevel: {config.optlevel2}" + # if config.solvent == "gas": print("\nCalculating single-point energies!") @@ -227,6 +250,7 @@ def part3(config, conformers, store_confs, ensembledata): "energy": 0.0, "energy2": 0.0, "success": False, + "onlyread": config.onlyread, } if config.multitemp: instruction["trange"] = [ @@ -240,8 +264,6 @@ def part3(config, conformers, store_confs, ensembledata): instruction["method"], _ = config.get_method_name( "sp", func=instruction["func"], basis=instruction["basis"] ) - if config.prog3 == "orca": - instruction["progpath"] = config.external_paths["orcapath"] folder = "part3" name = "highlevel single-point" else: @@ -271,9 +293,6 @@ def part3(config, conformers, store_confs, ensembledata): elif instruction["sm"] in ("gbsa_gsolv", "alpb_gsolv"): # do DFT gas phase sp and additive Gsolv instruction["jobtype"] = instruction["sm"] - if config.prog3 == "orca": - instruction["progpath"] = config.external_paths["orcapath"] - instruction["xtb_driver_path"] = config.external_paths["xtbpath"] instruction["method"], instruction["method2"] = config.get_method_name( instruction["jobtype"], func=instruction["func"], @@ -297,7 +316,6 @@ def part3(config, conformers, store_confs, ensembledata): job = OrcaJob instruction["prepinfo"] = ["high"] instruction["jobtype"] = instruction["sm"] - instruction["progpath"] = config.external_paths["orcapath"] instruction["method"], instruction["method2"] = config.get_method_name( "smd_gsolv", func=instruction["func"], @@ -310,8 +328,6 @@ def part3(config, conformers, store_confs, ensembledata): # with implicit solvation instruction["jobtype"] = "sp_implicit" instruction["prepinfo"] = ["high"] - if config.prog3 == "orca": - instruction["progpath"] = config.external_paths["orcapath"] instruction["method"], instruction["method2"] = config.get_method_name( "sp_implicit", func=instruction["func"], @@ -320,7 +336,7 @@ def part3(config, conformers, store_confs, ensembledata): ) name = "high level single-point" folder = "part3" - + if prev_calculated: check_for_folder(config.cwd, [i.id for i in prev_calculated], folder) print("The calculation was performed before for:") @@ -337,10 +353,10 @@ def part3(config, conformers, store_confs, ensembledata): config.cwd, calculate, folder, save_errors, store_confs ) # write coord to folder - if config.smgsolv3 in ('cosmors', 'cosmors-fine'): + if config.smgsolv3 in ("cosmors", "cosmors-fine"): cp_to = ["part3", folder] else: - cp_to = [folder,] + cp_to = [folder] for folder in cp_to: if geometries_from_input: # working on DFT unoptimized geometries @@ -350,13 +366,18 @@ def part3(config, conformers, store_confs, ensembledata): else: # need to copy optimized coord to folder for conf in list(calculate): - tmp1 = os.path.join(config.cwd, "CONF" + str(conf.id), config.func, "coord") + tmp1 = os.path.join( + config.cwd, "CONF" + str(conf.id), config.func, "coord" + ) tmp2 = os.path.join("CONF" + str(conf.id), folder, "coord") try: shutil.copy(tmp1, tmp2) except FileNotFoundError: - print_errors(f"{'ERROR:':{WARNLEN}}can't copy optimized " - f"geometry of CONF{conf.id}!", save_errors) + print_errors( + f"{'ERROR:':{WARNLEN}}can't copy optimized " + f"geometry of CONF{conf.id}!", + save_errors, + ) store_confs.append(calculate.pop(calculate.index(conf))) if config.solvent == "gas": @@ -375,7 +396,7 @@ def part3(config, conformers, store_confs, ensembledata): calculate, instruction, config.balance, - folder + folder, ) # check if too many calculations failed @@ -471,6 +492,36 @@ def part3(config, conformers, store_confs, ensembledata): f"{conf.highlevel_gsolv_info['energy']:>.8f}" ) calculate.append(prev_calculated.pop(prev_calculated.index(conf))) + + # create data for possible SI generation:----------------------------------- + # Energy: + ensembledata.si["part3"]["Energy"] = instruction["method"] + if "DOGCP" in instruction["prepinfo"]: + ensembledata.si["part3"]["Energy"] += " + GCP" + # Energy_settings: + if job == TmJob: + ensembledata.si["part3"]["Energy_settings"] = " ".join( + qm_prepinfo["tm"][instruction["prepinfo"][0]] + ).replace("-", "") + elif job == OrcaJob: + ensembledata.si["part3"]["Energy_settings"] = " ".join( + qm_prepinfo["orca"][instruction["prepinfo"][0]] + ) + # GmRRHO is done below! + # Solvation: + if config.solvent == "gas": + ensembledata.si["part3"]["G_solv"] = "gas-phase" + else: + ensembledata.si["part3"]["G_solv"] = instruction["method2"] + # Geometry is done above: + # QM-CODE: + ensembledata.si["part3"]["main QM code"] = str(config.prog3).upper() + # Threshold: + ensembledata.si["part3"][ + "Threshold" + ] = f"Boltzmann sum threshold: {config.part3_threshold} %" + # END SI generation -------------------------------------------------------- + for conf in calculate: conf.reset_job_info() @@ -499,7 +550,8 @@ def part3(config, conformers, store_confs, ensembledata): "progpath": config.external_paths["xtbpath"], "imagthr": config.imagthr, "sthr": config.sthr, - "scale":config.scale, + "scale": config.scale, + "onlyread": config.onlyread, } folder_rrho = "rrho_part3" instruction_rrho["method"], _ = config.get_method_name( @@ -509,6 +561,16 @@ def part3(config, conformers, store_confs, ensembledata): sm=instruction_rrho["sm_rrho"], solvent=instruction_rrho["solvent"], ) + + # GmRRHO for SI information: + if config.evaluate_rrho: + ensembledata.si["part3"]["G_mRRHO"] = instruction_rrho["method"] + if "bhess" in ensembledata.si["part3"]["G_mRRHO"]: + ensembledata.si["part3"]["G_mRRHO"] += " SPH" + else: + ensembledata.si["part3"]["G_mRRHO"] = "not included" + # END SI + if config.multitemp: instruction_rrho["trange"] = [ i for i in frange(config.trange[0], config.trange[1], config.trange[2]) @@ -516,7 +578,6 @@ def part3(config, conformers, store_confs, ensembledata): else: instruction_rrho["trange"] = [] - # check if calculated for conf in list(calculate): if conf.removed: @@ -526,16 +587,16 @@ def part3(config, conformers, store_confs, ensembledata): if conf.highlevel_grrho_info["info"] == "calculated": conf = calculate.pop(calculate.index(conf)) conf.job["success"] = True + conf.sym = conf.highlevel_grrho_info.get("sym", conf.sym) + conf.linear = conf.highlevel_grrho_info.get("linear", conf.linear) + conf.symnum = conf.highlevel_grrho_info.get("symnum", conf.symnum) prev_calculated.append(conf) elif conf.highlevel_grrho_info["info"] == "failed": conf = calculate.pop(calculate.index(conf)) conf.part_info["part3"] = "refused" store_confs.append(conf) print(f"Calculation of CONF{conf.id} failed in the previous run!") - elif conf.highlevel_grrho_info["info"] in ( - "not_calculated", - "prep-failed", - ): + elif conf.highlevel_grrho_info["info"] in ("not_calculated", "prep-failed"): # stay in calculate (e.g not_calculated or prep-failed) # check if method has been calculated in part2 if instruction_rrho["method"] == conf.lowlevel_grrho_info["method"]: @@ -549,8 +610,13 @@ def part3(config, conformers, store_confs, ensembledata): getattr(conf, "highlevel_grrho_info").update(tmp) # if rrho is taken from part2 still make folder for rrho_part3 save_errors, store_confs, calculate = new_folders( - config.cwd, [conf,], folder_rrho, save_errors, store_confs, - silent=True) + config.cwd, + [conf], + folder_rrho, + save_errors, + store_confs, + silent=True, + ) prev_calculated.append(calculate.pop(calculate.index(conf))) elif ( instruction_rrho["method"] @@ -565,8 +631,13 @@ def part3(config, conformers, store_confs, ensembledata): ) # if rrho is taken from part2 still make folder for rrho_part3 save_errors, store_confs, calculate = new_folders( - config.cwd, [conf,], folder_rrho, save_errors, store_confs, - silent=True) + config.cwd, + [conf], + folder_rrho, + save_errors, + store_confs, + silent=True, + ) prev_calculated.append(calculate.pop(calculate.index(conf))) else: # stays in calculate and has to be calculated now @@ -624,7 +695,8 @@ def part3(config, conformers, store_confs, ensembledata): config.cwd, "CONF" + str(conf.id), folder_rrho ) shutil.copy( - os.path.join(tmp_from, "coord"), os.path.join(tmp_to, "coord") + os.path.join(tmp_from, "coord"), + os.path.join(tmp_to, "coord"), ) except shutil.SameFileError: pass @@ -635,11 +707,15 @@ def part3(config, conformers, store_confs, ensembledata): "The corresponding file does not exist." ) elif not os.path.isdir(tmp_to): - print(f"{'ERROR:':{WARNLEN}}Could not create folder {tmp_to}!") + print( + f"{'ERROR:':{WARNLEN}}Could not create folder {tmp_to}!" + ) print(f"{'ERROR:':{WARNLEN}}Removing conformer {conf.name}!") conf.highlevel_grrho_info["info"] = "prep-failed" store_confs.append(calculate.pop(calculate.index(conf))) - save_errors.append(f"CONF{conf.id} was removed, because IO failed!") + save_errors.append( + f"CONF{conf.id} was removed, because IO failed!" + ) # parallel execution: calculate = run_in_parallel( config, @@ -681,6 +757,9 @@ def part3(config, conformers, store_confs, ensembledata): else: conf.sym = conf.job["symmetry"] conf.linear = conf.job["linear"] + conf.highlevel_grrho_info["sym"] = conf.job["symmetry"] + conf.highlevel_grrho_info["linear"] = conf.job["linear"] + conf.highlevel_grrho_info["symnum"] = conf.job["symnum"] conf.highlevel_grrho_info["rmsd"] = conf.job["rmsd"] conf.highlevel_grrho_info["energy"] = conf.job["energy"] conf.highlevel_grrho_info["info"] = "calculated" @@ -742,7 +821,7 @@ def part3(config, conformers, store_confs, ensembledata): lambda conf: getattr(conf, "rel_xtb_energy"), lambda conf: getattr(conf, "highlevel_sp_info")["energy"], lambda conf: getattr(conf, "highlevel_gsolv_info")["energy"], - #lambda conf: getattr(conf, "highlevel_grrho_info")["energy"], + # lambda conf: getattr(conf, "highlevel_grrho_info")["energy"], lambda conf: conf.get_mrrho( config.temperature, "highlevel_grrho_info", config.consider_sym ), @@ -764,14 +843,14 @@ def part3(config, conformers, store_confs, ensembledata): columndescription = [ "", "[a.u.]", - "[kcal/mol]", + "[kcal/mol]", "", "", "", - "[Eh]", + "[Eh]", "[kcal/mol]", f" % at {config.temperature:.2f} K", - ] + ] columndescription2 = ["", "", "", "", "", "", "", "", ""] columnformat = [ "", @@ -782,8 +861,8 @@ def part3(config, conformers, store_confs, ensembledata): (12, 7), (12, 7), (5, 3), - (5, 2) - ] + (5, 2), + ] if config.solvent == "gas": columndescription[3] = instruction["method"] elif config.solvent != "gas": @@ -817,12 +896,12 @@ def part3(config, conformers, store_confs, ensembledata): solv = "highlevel_gsolv_info" e = "highlevel_sp_info" conf.calc_free_energy( - e=e, - solv=solv, + e=e, + solv=solv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) try: @@ -843,20 +922,16 @@ def part3(config, conformers, store_confs, ensembledata): columndescription2=columndescription2, ) # end printout for part3 - + conf_in_interval(calculate) + # -------------------------------------------------------------------------- for conf in calculate: if conf.free_energy == minfree: ensembledata.bestconf["part3"] = conf.id -################################################################################ + ################################################################################ # calculate average G correction print("\nCalculating Boltzmann averaged free energy of ensemble!\n") - avGcorrection = { - "avG": {}, - "avE": {}, - "avGsolv": {}, - "avGRRHO": {}, - } + avGcorrection = {"avG": {}, "avE": {}, "avGsolv": {}, "avGRRHO": {}} if config.multitemp: trange = [ i for i in frange(config.trange[0], config.trange[1], config.trange[2]) @@ -905,18 +980,12 @@ def part3(config, conformers, store_confs, ensembledata): solv = "highlevel_gsolv_info" e = "highlevel_sp_info" conf.calc_free_energy( - e=e, - solv=solv, - rrho=rrho, + e=e, + solv=solv, + rrho=rrho, t=temperature, - consider_sym=config.consider_sym - ) - # try: - # minfreeT = min( - # [conf.free_energy for conf in calculate if conf.free_energy is not None] - # ) - # except ValueError: - # raise ValueError + consider_sym=config.consider_sym, + ) calculate = calc_boltzmannweights(calculate, "free_energy", temperature) avG = 0.0 @@ -927,13 +996,11 @@ def part3(config, conformers, store_confs, ensembledata): for conf in calculate: avG += conf.bm_weight * conf.free_energy avE += conf.bm_weight * conf.highlevel_sp_info["energy"] - #avGRRHO += conf.bm_weight * conf.highlevel_grrho_info["range"].get( + # avGRRHO += conf.bm_weight * conf.highlevel_grrho_info["range"].get( # temperature, 0.0 - #) + # ) avGRRHO += conf.bm_weight * conf.get_mrrho( - temperature, - "highlevel_grrho_info", - config.consider_sym + temperature, "highlevel_grrho_info", config.consider_sym ) avGsolv += conf.bm_weight * conf.highlevel_gsolv_info["range"].get( temperature, 0.0 @@ -973,13 +1040,15 @@ def part3(config, conformers, store_confs, ensembledata): else: print(line) print("".ljust(int(PLENGTH), "-")) - print("avGcorrection(T)* = Correction to free energy, which can be added (manually) to avG(T).") + print( + "avGcorrection(T)* = Correction to free energy, which can be added (manually) to avG(T)." + ) print(" If only a small ensemble is evaluated in part3 ") print(" this can incorporate information of higher lying ") print(" conformers of a more complete ensemble from part2.") print("") -################################################################################ + ################################################################################ # -----------------------------Trange Ouput END------------------------------ # reset boltzmannweights to correct temperature # get free energy at (T) @@ -994,12 +1063,12 @@ def part3(config, conformers, store_confs, ensembledata): solv = "highlevel_gsolv_info" e = "highlevel_sp_info" conf.calc_free_energy( - e=e, + e=e, solv=solv, rrho=rrho, t=config.temperature, - consider_sym=config.consider_sym - ) + consider_sym=config.consider_sym, + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) # SORTING for the next part: @@ -1075,13 +1144,13 @@ def part3(config, conformers, store_confs, ensembledata): ) as best: best.write( "$coord # {} {} !CONF{} \n".format( - conf.free_energy, + conf.free_energy, conf.get_mrrho( - config.temperature, - rrho="highlevel_grrho_info", - consider_sym=config.consider_sym + config.temperature, + rrho="highlevel_grrho_info", + consider_sym=config.consider_sym, ), - conf.id + conf.id, ) ) for line in coord[1:]: @@ -1098,18 +1167,11 @@ def part3(config, conformers, store_confs, ensembledata): conf.reset_job_info() if save_errors: - print( - "\n***---------------------------------------------------------***" - ) - print( - "Printing most relevant errors again, just for user convenience:" - ) + print("\n***---------------------------------------------------------***") + print("Printing most relevant errors again, just for user convenience:") for _ in list(save_errors): print(save_errors.pop()) - print( - "***---------------------------------------------------------***" - - ) + print("***---------------------------------------------------------***") tmp = int((PLENGTH - len("END of Part3")) / 2) print("\n" + "".ljust(tmp, ">") + "END of Part3" + "".rjust(tmp, "<")) if config.progress: diff --git a/censo_qm/setupcenso.py b/censo_qm/setupcenso.py index 10b518c..a7a37af 100755 --- a/censo_qm/setupcenso.py +++ b/censo_qm/setupcenso.py @@ -6,15 +6,7 @@ import sys import json from collections import OrderedDict -from .cfg import ( - CODING, - PLENGTH, - DESCR, - WARNLEN, - censo_solvent_db, - __version__, - NmrRef, -) +from .cfg import CODING, PLENGTH, DESCR, WARNLEN, censo_solvent_db, __version__, NmrRef from .inputhandling import config_setup, internal_settings from .datastructure import MoleculeData from .qm_job import QmJob @@ -57,7 +49,7 @@ def enso_startup(cwd, args): configfname = ".censorc" if args.writeconfig: - #check if .censorc in local or home dir and ask user it the program + # check if .censorc in local or home dir and ask user it the program # paths should be copied tmp = None newconfigfname = "censorc_new" @@ -67,14 +59,18 @@ def enso_startup(cwd, args): tmp = os.path.join(os.path.expanduser("~"), configfname) if tmp is not None: print(f"An existing .censorc has been found in {tmp}") - print(f"Do you want to copy existing program path information to the new remote configuration file?") + print( + f"Do you want to copy existing program path information to the new remote configuration file?" + ) print("Please type 'yes' or 'no':") user_input = input() - if user_input.strip() in ('y', 'yes'): + if user_input.strip() in ("y", "yes"): config.read_program_paths(tmp) - config.write_rcfile(os.path.join(config.cwd, newconfigfname), usepaths=True) + config.write_rcfile( + os.path.join(config.cwd, newconfigfname), usepaths=True + ) print("") - elif user_input.strip() in ('n', 'no'): + elif user_input.strip() in ("n", "no"): config.write_rcfile(os.path.join(config.cwd, newconfigfname)) else: config.write_rcfile(os.path.join(config.cwd, newconfigfname)) @@ -96,9 +92,26 @@ def enso_startup(cwd, args): else: raise FileNotFoundError except FileNotFoundError: - print(f"{'ERROR:':{WARNLEN}}Could not find the configuration file: {configfname}.") + print( + f"{'ERROR:':{WARNLEN}}Could not find the configuration file: {configfname}." + ) print("\nGoing to exit!") sys.exit(1) + elif ( + args.restart + and os.path.isfile(os.path.join(config.cwd, "enso.json")) + and os.path.isfile( + config.read_json(os.path.join(config.cwd, "enso.json"), silent=True) + .get("settings", {}) + .get("configpath", "") + ) + ): + # read configpath from previous enso.json if restart is requested + config.configpath = ( + config.read_json(os.path.join(config.cwd, "enso.json"), silent=True) + .get("settings", {}) + .get("configpath", "") + ) elif os.path.isfile(os.path.join(config.cwd, configfname)): # local configuration file before remote configuration file config.configpath = os.path.join(config.cwd, configfname) @@ -191,21 +204,26 @@ def enso_startup(cwd, args): myfile.seek(0) # reset reader except (ValueError, KeyError, AttributeError) as e: print(e) - print(f"{'ERROR:':{WARNLEN}}Please create a new .censorc --> 'censo -newconfig'") + print( + f"{'ERROR:':{WARNLEN}}Please create a new .censorc --> 'censo -newconfig'" + ) sys.exit(1) if not startread in myfile.read(): - print(f"{'ERROR:':{WARNLEN}}You are using a corrupted .censorc. Create a new one!") + print( + f"{'ERROR:':{WARNLEN}}You are using a corrupted .censorc. Create a new one!" + ) sys.exit(1) config.read_config(config.configpath, startread, args) if args.copyinput: config.read_program_paths(config.configpath) config.write_censo_inp(config.cwd) - print("The file censo.inp with the current settings has been written to the current working directory.") + print( + "The file censo.inp with the current settings has been written to the current working directory." + ) print("\nGoing to exit!") sys.exit(1) - # read inputfile: if os.path.isfile(os.path.join(config.cwd, "enso.json")): tmp = config.read_json(os.path.join(config.cwd, "enso.json"), silent=True) @@ -278,9 +296,9 @@ def enso_startup(cwd, args): error_logical = config.check_logic() # printing parameters - config.print_parameters() + config.print_parameters(cmlcall=sys.argv) config.read_program_paths(config.configpath) - requirements = config.needed_external_programs(config) + requirements = config.needed_external_programs() error_logical = config.processQMpaths(requirements, error_logical) if error_logical and not args.debug and args.checkinput: @@ -291,8 +309,7 @@ def enso_startup(cwd, args): print("\nGoing to exit!") sys.exit(1) elif error_logical and not args.debug: - print(f"\n{'ERROR:':{WARNLEN}}CENSO can not continue due to input errors!" - ) + print(f"\n{'ERROR:':{WARNLEN}}CENSO can not continue due to input errors!") print("\nGoing to exit!") sys.exit(1) @@ -403,7 +420,9 @@ def enso_startup(cwd, args): id=save_data[conf].get("id"), filename=save_data[conf].get("filename"), part_info=save_data[conf].get("part_info"), - previous_part_info=save_data[conf].get("previous_part_info", EnsembleData().previous_part_info), + previous_part_info=save_data[conf].get( + "previous_part_info", EnsembleData().previous_part_info + ), avGcorrection=save_data[conf].get("avGcorrection"), comment=save_data[conf].get("comment"), bestconf=save_data[conf].get("bestconf"), @@ -427,8 +446,20 @@ def enso_startup(cwd, args): rel_xtb_energy=save_data[conf].get("rel_xtb_energy"), rel_xtb_free_energy=save_data[conf].get("rel_xtb_free_energy"), sym=save_data[conf].get("sym", getattr(MoleculeData(0), "sym")), - linear=save_data[conf].get("linear", getattr(MoleculeData(0), "linear")), - symnum=save_data[conf].get("symnum", getattr(MoleculeData(0), "symnum")), + linear=save_data[conf].get( + "linear", getattr(MoleculeData(0), "linear") + ), + symnum=save_data[conf].get( + "symnum", + MoleculeData(0)._get_sym_num( + sym=save_data[conf].get( + "sym", getattr(MoleculeData(0), "sym") + ), + linear=save_data[conf].get( + "linear", getattr(MoleculeData(0), "linear") + ), + ), + ), gi=save_data[conf].get("gi", getattr(MoleculeData(0), "sym")), removed=save_data[conf].get( "removed", getattr(MoleculeData(0), "removed") @@ -665,35 +696,31 @@ def enso_startup(cwd, args): prog=config.prog, basis=config.basis_or, func=config.func_or, + solvent=config.solvent, + sm="cosmo", func2=config.func_or_scf, ) molecule.load_prev("optical_rotation_info", method) ##### - elif value and key in ( - "func_j", - "basis_j", - "sm4_j" - ): + elif value and key in ("func_j", "basis_j", "sm4_j"): # reset to default, load is not available - molecule.nmr_coupling_info=getattr(MoleculeData(0), "nmr_coupling_info") - elif value and key in ( - "func_s", - "basis_s", - "sm4_s" - ): + molecule.nmr_coupling_info = getattr( + MoleculeData(0), "nmr_coupling_info" + ) + elif value and key in ("func_s", "basis_s", "sm4_s"): # reset to default, load is not available - molecule.nmr_shielding_info=getattr(MoleculeData(0), "nmr_shielding_info") + molecule.nmr_shielding_info = getattr( + MoleculeData(0), "nmr_shielding_info" + ) ##### elif value and key in ("func3", "basis3"): # save calculated to molecule.save_prev( "highlevel_sp_info", - getattr(molecule, "highlevel_sp_info").get( - "method" - ), + getattr(molecule, "highlevel_sp_info").get("method"), ) # load new if available - if config.smgsolv3 in ('cpcm', 'smd', 'cosmo', 'dcosmors'): + if config.smgsolv3 in ("cpcm", "smd", "cosmo", "dcosmors"): expected = "sp_implicit" else: expected = "sp" @@ -702,7 +729,7 @@ def enso_startup(cwd, args): prog=config.prog3, basis=config.basis3, func=config.func3, - sm=config.smgsolv3 + sm=config.smgsolv3, ) molecule.load_prev("highlevel_sp_info", method) # finally add molecule to list @@ -724,7 +751,9 @@ def enso_startup(cwd, args): elif not args.checkinput: # don't create enso.json on checkinput # enso.json does not exist, create new conformers - print(f"{'INFORMATION:':{WARNLEN}}No restart information exists and is created during this run!\n") + print( + f"{'INFORMATION:':{WARNLEN}}No restart information exists and is created during this run!\n" + ) conformers = [] for i in range(1, config.nconf + 1): conformers.append(QmJob(i)) @@ -746,6 +775,7 @@ def enso_startup(cwd, args): print(f"{'ERROR:':{WARNLEN}}No conformers are considered!") print("\nGoing to exit!") sys.exit(1) + # formatting information: config.lenconfx = max([len(str(i.id)) for i in conformers]) conformers.sort(key=lambda x: int(x.id)) diff --git a/censo_qm/tm_job.py b/censo_qm/tm_job.py index c54b90c..f4b726d 100644 --- a/censo_qm/tm_job.py +++ b/censo_qm/tm_job.py @@ -20,7 +20,8 @@ external_paths, cosmors_param, composite_dfa, - composite_method_basis + composite_method_basis, + hybrid_dfa, ) from .utilities import last_folders, print from .qm_job import QmJob @@ -242,7 +243,9 @@ def _prep_cefine(self): self.job["internal_error"].append(line) if self.job["solvent"] != "gas" and self.job["sm"] in ("cosmo", "dcosmors"): if solvent_dcosmors.get(self.job["solvent"], "not found!") == "not found!": - print(f"{'ERROR:':{WARNLEN}} Solvent {self.job['solvent']} is not known for cefine!") + print( + f"{'ERROR:':{WARNLEN}} Solvent {self.job['solvent']} is not known for cefine!" + ) self.job["success"] = False self.job["internal_error"].append("prep-failed") return @@ -269,7 +272,9 @@ def _prep_cefine(self): controlappend.append("$rpaconv 4") if dogcp: - if self.job['func'] in composite_dfa and self.job['basis'] == composite_method_basis.get(self.job["func"],'NONE'): + if self.job["func"] in composite_dfa and self.job[ + "basis" + ] == composite_method_basis.get(self.job["func"], "NONE"): pass else: if self.job["basis"] == "def2-SV(P)": @@ -412,18 +417,18 @@ def _prep_cefine(self): # ****************************end cefine************************************ - def _sp(self, silent=False): + def _sp(self, silent=False, outfile="ridft.out"): """ Turbomole single-point calculation, needs previous cefine run """ + fullpath = os.path.join(self.job["workdir"], outfile) if not self.job["onlyread"]: + # perform calculation: if not silent: print( f"Running single-point in {last_folders(self.job['workdir'], 2):18}" ) - with open( - os.path.join(self.job["workdir"], "ridft.out"), "w", newline=None - ) as outputfile: + with open(fullpath, "w", newline=None) as outputfile: subprocess.call( ["ridft"], shell=False, @@ -434,15 +439,11 @@ def _sp(self, silent=False): stdout=outputfile, env=ENVIRON, ) + # process output: time.sleep(0.02) # check if scf is converged: - if os.path.isfile(os.path.join(self.job["workdir"], "ridft.out")): - with open( - os.path.join(self.job["workdir"], "ridft.out"), - "r", - encoding=CODING, - newline=None, - ) as inp: + if os.path.isfile(fullpath): + with open(fullpath, "r", encoding=CODING, newline=None) as inp: stor = inp.readlines() if " ENERGY CONVERGED !\n" not in stor: print( @@ -453,36 +454,33 @@ def _sp(self, silent=False): self.job["energy"] = 0.0 self.job["energy2"] = 0.0 return + # get energy: + for line in stor: + if "| total energy = " in line: + try: + self.job["energy"] = float(line.strip().split()[4]) + self.job["success"] = True + if self.job["jobtype"] == "sp_implicit": + self.job["energy2"] = 0.0 + except (ValueError, IndexError): + print( + f"{'ERROR:':{WARNLEN}}while converting energy " + f"in: {last_folders(self.job['workdir'], 2):18}" + ) + self.job["success"] = False + self.job["energy"] = 0.0 + self.job["energy2"] = 0.0 + return + if not self.job["success"]: + self.job["energy"] = 0.0 + self.job["energy2"] = 0.0 + return else: - print( - f"{'WARNING:':{WARNLEN}}" - f"{os.path.join(self.job['workdir'], 'ridft.out')} doesn't exist!" - ) + print(f"{'WARNING:':{WARNLEN}}{fullpath} doesn't exist!") self.job["success"] = False self.job["energy"] = 0.0 self.job["energy2"] = 0.0 return - if os.path.isfile(os.path.join(self.job["workdir"], "energy")): - with open( - os.path.join(self.job["workdir"], "energy"), - "r", - encoding=CODING, - newline=None, - ) as energy: - storage = energy.readlines() - try: - self.job["energy"] = float(storage[-2].split()[1]) - self.job["success"] = True - except ValueError: - print( - f"{'ERROR:':{WARNLEN}}while converting energy " - f"in: {last_folders(self.job['workdir'], 2):18}" - ) - if self.job["jobtype"] == "sp_implicit": - self.job["energy2"] = 0.0 - else: - self.job["energy"] = 0.0 - self.job["success"] = False # ****************************end _sp*************************************** @@ -503,7 +501,7 @@ def _cosmors(self): f"{last_folders(self.job['workdir'], 3):18}" ) # parametrisation - # ctd and database have to match fine -> fine or non-fine -> nofine + # ctd and database have to match fine -> fine or non-fine -> non-fine if self.job["ctd-param"] == "automatic": if self.job["cosmorsparam"] == "fine": if ( @@ -511,11 +509,23 @@ def _cosmors(self): and "1201" not in self.job["cosmorssetup"] and "1301" not in self.job["cosmorssetup"] ): - self.job["cosmorssetup"] = self.job["cosmorssetup"].replace("BP_TZVP", "BP_TZVPD_FINE") - elif "FINE" not in self.job["cosmorssetup"] and "1201" in self.job["cosmorssetup"]: - self.job["cosmorssetup"] = self.job["cosmorssetup"].replace("BP_TZVP", "BP_TZVPD_FINE_HB2012") - elif "FINE" not in self.job["cosmorssetup"] and "1301" in self.job["cosmorssetup"]: - self.job["cosmorssetup"] = self.job["cosmorssetup"].replace("BP_TZVP", "BP_TZVPD_FINE_HB2012") + self.job["cosmorssetup"] = self.job["cosmorssetup"].replace( + "BP_TZVP", "BP_TZVPD_FINE" + ) + elif ( + "FINE" not in self.job["cosmorssetup"] + and "1201" in self.job["cosmorssetup"] + ): + self.job["cosmorssetup"] = self.job["cosmorssetup"].replace( + "BP_TZVP", "BP_TZVPD_FINE_HB2012" + ) + elif ( + "FINE" not in self.job["cosmorssetup"] + and "1301" in self.job["cosmorssetup"] + ): + self.job["cosmorssetup"] = self.job["cosmorssetup"].replace( + "BP_TZVP", "BP_TZVPD_FINE_HB2012" + ) elif self.job["cosmorsparam"] == "normal": if "FINE" in self.job["cosmorssetup"]: ## normal cosmors @@ -537,6 +547,7 @@ def _cosmors(self): tmp_new = tmp_new + " ".join(tmp_dat[tmp_dat.index(item) + 1 :]) self.job["cosmorssetup"] = tmp_new # run two single-points: + # copy mos first if self.job["copymos"]: if self.job["unpaired"] > 0: molist = ["alpha", "beta"] @@ -575,6 +586,12 @@ def _cosmors(self): f"{last_folders(self.job['workdir'], 3):18}" ) return + # try: + # tmp_from = os.path.join(self.job["workdir"], "ridft.out") + # tmp_to = os.path.join(self.job["workdir"], "gas_"+ self.job["partx"] +".out") + # shutil.copy(tmp_from, tmp_to) + # except (FileNotFoundError, KeyError): + # pass with open( os.path.join(self.job["workdir"], "out.energy"), "w", newline=None ) as out: @@ -730,72 +747,72 @@ def _cosmors(self): env=ENVIRON, ) time.sleep(0.1) - # get T and Gsolv for version > cosmothermX16 - ## volumework: - R = 1.987203585e-03 # kcal/(mol*K) - videal = ( - 24.789561955 / 298.15 - ) # molar volume for ideal gas at 298.15 K 100.0 kPa - gsolvt = {} - try: - with open( - os.path.join(self.job["workdir"], "cosmotherm.tab"), - "r", - encoding=CODING, - newline=None, - ) as inp: - stor = inp.readlines() - for line in stor: - if "T=" in line: - T = float(line.split()[5]) - vwork = R * T * math.log(videal * T) - elif " out " in line: - gsolvt[T] = float(line.split()[-1]) / AU2KCAL + vwork / AU2KCAL - self.job["erange1"] = gsolvt - except (FileNotFoundError, ValueError): - print( - f"{'ERROR:':{WARNLEN}}cosmotherm.tab was not written, this error can be " - "due to a missing licensefile information, or wrong path " - "to the COSMO-RS Database." - ) - self.job["energy"] = 0.0 - self.job["energy2"] = 0.0 - self.job["erange1"] = {} - self.job["success"] = False - return - except IndexError: - print(f"{'ERROR:':{WARNLEN}}IndexERROR in cosmotherm.tab!") - self.job["energy"] = 0.0 - self.job["energy2"] = 0.0 - self.job["erange1"] = {} + # calculations done, get output + if self.job["onlyread"]: + # if self.job.get("partx", None) is None: + # get gas phase energy: + gaspath = os.path.join(self.job["workdir"], "out.energy") + if os.path.exists(gaspath) and os.path.getsize(gaspath) > 0: + with open(gaspath, "r", newline=None) as inp: + try: + gas_phase_energy = float(inp.readline().strip()) + except Exception: + self.job["success"] = False + return + else: self.job["success"] = False return - # cosmothermrd - if ( - os.stat(os.path.join(self.job["workdir"], "cosmotherm.tab")).st_size - == 0 - ): - print( - f"{'ERROR:':{WARNLEN}}cosmotherm.tab was not written, this error can be " - "due to a missing licensefile information, or wrong path " - "to the COSMO-RS Database." - ) - self.job["energy"] = 0.0 - self.job["energy2"] = 0.0 - self.job["erange1"] = {} - self.job["success"] = False - gsolv_out = 0.0 - for temp in gsolvt.keys(): - if isclose(self.job["temperature"], temp, abs_tol=0.6): - gsolv_out = gsolvt[temp] - temp = float(self.job["temperature"]) - ## volumework: - R = 1.987203585e-03 # kcal/(mol*K) - videal = ( - 24.789561955 / 298.15 - ) # molar volume for ideal gas at 298.15 K 100.0 kPa - volwork = R * temp * math.log(videal * temp) + # get T and Gsolv for version > cosmothermX16 + ## volumework: + R = 1.987203585e-03 # kcal/(mol*K) + videal = ( + 24.789561955 / 298.15 + ) # molar volume for ideal gas at 298.15 K 100.0 kPa + gsolvt = {} + + cosmothermtab = os.path.join(self.job["workdir"], "cosmotherm.tab") + try: + with open(cosmothermtab, "r", encoding=CODING, newline=None) as inp: + stor = inp.readlines() + for line in stor: + if "T=" in line: + T = float(line.split()[5]) + vwork = R * T * math.log(videal * T) + elif " out " in line: + gsolvt[T] = float(line.split()[-1]) / AU2KCAL + vwork / AU2KCAL + self.job["erange1"] = gsolvt + except (FileNotFoundError, ValueError): + print( + f"{'ERROR:':{WARNLEN}}cosmotherm.tab was not written, this error can be " + "due to a missing licensefile information, or wrong path " + "to the COSMO-RS Database." + ) + self.job["energy"] = 0.0 + self.job["energy2"] = 0.0 + self.job["erange1"] = {} + self.job["success"] = False + return + except IndexError: + print(f"{'ERROR:':{WARNLEN}}IndexERROR in cosmotherm.tab!") + self.job["energy"] = 0.0 + self.job["energy2"] = 0.0 + self.job["erange1"] = {} + self.job["success"] = False + return + # cosmothermrd + gsolv_out = 0.0 + for temp in gsolvt.keys(): + if isclose(self.job["temperature"], temp, abs_tol=0.6): + gsolv_out = gsolvt[temp] + temp = float(self.job["temperature"]) + ## volumework: + R = 1.987203585e-03 # kcal/(mol*K) + videal = ( + 24.789561955 / 298.15 + ) # molar volume for ideal gas at 298.15 K 100.0 kPa + volwork = R * temp * math.log(videal * temp) + if not self.job["onlyread"]: with open( os.path.join(os.path.dirname(self.job["workdir"]), "cosmors.out"), "w", @@ -820,10 +837,10 @@ def _cosmors(self): ) ) time.sleep(0.01) - self.job["energy"] = gas_phase_energy - self.job["energy2"] = gsolv_out # VOLWORK INCLUDED - self.job["erange1"][self.job["temperature"]] = gsolv_out # VOLWORK INCLUDED - self.job["success"] = True + self.job["energy"] = gas_phase_energy + self.job["energy2"] = gsolv_out # VOLWORK INCLUDED + self.job["erange1"][self.job["temperature"]] = gsolv_out # VOLWORK INCLUDED + self.job["success"] = True # ********************************end _cosmors*********************************** def _xtbopt(self): @@ -835,6 +852,7 @@ def _xtbopt(self): output = "opt-part2.out" else: output = "opt-part1.out" + outputpath = os.path.join(self.job["workdir"], output) if not self.job["onlyread"]: print(f"Running optimization in {last_folders(self.job['workdir'], 2):18}") files = [ @@ -850,7 +868,7 @@ def _xtbopt(self): os.remove(os.path.join(self.job["workdir"], file)) callargs = [ - self.job["xtb_driver_path"], + external_paths["xtbpath"], "coord", "--opt", self.job["optlevel"], @@ -875,9 +893,7 @@ def _xtbopt(self): callargs.append("-I") callargs.append("opt.inp") time.sleep(0.02) - with open( - os.path.join(self.job["workdir"], output), "w", newline=None - ) as outputfile: + with open(outputpath, "w", newline=None) as outputfile: returncode = subprocess.call( callargs, shell=False, @@ -895,14 +911,10 @@ def _xtbopt(self): f"{last_folders(self.job['workdir'], 2):18} not converged" ) time.sleep(0.02) + # read output: # check if converged: - if os.path.isfile(os.path.join(self.job["workdir"], output)): - with open( - os.path.join(self.job["workdir"], output), - "r", - encoding=CODING, - newline=None, - ) as inp: + if os.path.isfile(outputpath): + with open(outputpath, "r", encoding=CODING, newline=None) as inp: stor = inp.readlines() for line in stor: if ( @@ -923,23 +935,24 @@ def _xtbopt(self): elif "*** GEOMETRY OPTIMIZATION CONVERGED AFTER " in line: self.job["cycles"] += int(line.split()[5]) self.job["converged"] = True - with open( - os.path.join(self.job["workdir"], output), - "r", - encoding=CODING, - newline=None, - ) as inp: - for line in inp: - if "av. E: " in line: - # self.job["ecyc"].append(float(line.split("Eh")[0].split()[-1])) - self.job["ecyc"].append(float(line.split("->")[-1])) + ##### + for line in stor: + if "av. E: " in line and "->" in line: + try: + self.job["ecyc"].append(float(line.split("->")[-1])) + except ValueError as e: + error_logical = True + print(f"{'ERROR:':{WARNLEN}}in CONF{self.id} calculation:\n{e}") + break if " :: gradient norm " in line: - self.job["grad_norm"] = float(line.split()[3]) + try: + self.job["grad_norm"] = float(line.split()[3]) + except ValueError as e: + error_logical = True + print(f"{'ERROR:':{WARNLEN}}in CONF{self.id} calculation:\n{e}") + break else: - print( - f"{'WARNING:':{WARNLEN}}" - f"{os.path.join(self.job['workdir'], output)} doesn't exist!" - ) + print(f"{'WARNING:':{WARNLEN}}{outputpath} doesn't exist!") error_logical = True if not error_logical: try: @@ -1108,50 +1121,49 @@ def _nmr_coupling(self): print( f"Running couplings calculation in {last_folders(self.job['workdir'], 2)}" ) - # escf doesnot allow for mgrids! - with open( - os.path.join(self.job["workdir"], "control"), "r", newline=None - ) as inp: - tmp = inp.readlines() - with open( - os.path.join(self.job["workdir"], "control"), "w", newline=None - ) as out: - for line in tmp: - if "gridsize" in line: - out.write(f" gridsize {5} \n") - else: - out.write(line + "\n") + outputpath = os.path.join(self.job["workdir"], "escf.out") + if not self.job["onlyread"]: + # escf doesnot allow for mgrids! + with open( + os.path.join(self.job["workdir"], "control"), "r", newline=None + ) as inp: + tmp = inp.readlines() + with open( + os.path.join(self.job["workdir"], "control"), "w", newline=None + ) as out: + for line in tmp: + if "gridsize" in line: + out.write(f" gridsize {5} \n") + else: + out.write(line + "\n") - with open( - os.path.join(self.job["workdir"], "escf.out"), "w", newline=None - ) as outputfile: - subprocess.call( - [self.job["progpath"], "-smpcpus", str(self.job["omp"])], - shell=False, - stdin=None, - stderr=subprocess.STDOUT, - universal_newlines=False, - cwd=self.job["workdir"], - stdout=outputfile, - env=ENVIRON, - ) - time.sleep(0.02) - # check for convergence - with open( - os.path.join(self.job["workdir"], "escf.out"), - "r", - encoding=CODING, - newline=None, - ) as inp: - stor = inp.readlines() - if " **** escf : all done ****\n" in stor: - self.job["success"] = True - else: - print( - f"{'ERROR:':{WARNLEN}}coupling calculation failed in " - f"{last_folders(self.job['workdir'], 1):18}" + with open(outputpath, "w", newline=None) as outputfile: + subprocess.call( + [external_paths["escfpath"], "-smpcpus", str(self.job["omp"])], + shell=False, + stdin=None, + stderr=subprocess.STDOUT, + universal_newlines=False, + cwd=self.job["workdir"], + stdout=outputfile, + env=ENVIRON, ) - self.job["success"] = False + time.sleep(0.02) + # check for convergence + if os.path.isfile(outputpath): + with open(outputpath, "r", encoding=CODING, newline=None) as inp: + stor = inp.readlines() + if " **** escf : all done ****\n" in stor: + self.job["success"] = True + else: + print( + f"{'ERROR:':{WARNLEN}}coupling calculation failed in " + f"{last_folders(self.job['workdir'], 1):18}" + ) + self.job["success"] = False + else: + print(f"The file escf.out could not be found!") + self.job["success"] = False def _nmr_shielding(self): """ @@ -1162,44 +1174,40 @@ def _nmr_shielding(self): last_folders(self.job["workdir"], 2) ) ) - # update grid to m5! - with open( - os.path.join(self.job["workdir"], "control"), "r", newline=None - ) as inp: - tmp = inp.readlines() - with open( - os.path.join(self.job["workdir"], "control"), "w", newline=None - ) as out: - for line in tmp: - if "gridsize" in line: - out.write(f" gridsize {'m5'} \n") - if "$disp" in line and self.job["func"] in ("kt2", "kt1"): - pass - else: - out.write(line + "\n") - time.sleep(0.02) - - with open( - os.path.join(self.job["workdir"], "mpshift.out"), "w", newline=None - ) as outputfile: - subprocess.call( - [self.job["progpath"], "-smpcpus", str(self.job["omp"])], - shell=False, - stdin=None, - stderr=subprocess.STDOUT, - universal_newlines=False, - cwd=self.job["workdir"], - stdout=outputfile, - env=ENVIRON, - ) - time.sleep(0.02) - # check if shift calculation is converged: + outputpath = os.path.join(self.job["workdir"], "mpshift.out") + if not self.job["onlyread"]: + # update grid to m5! with open( - os.path.join(self.job["workdir"], "mpshift.out"), - "r", - encoding=CODING, - newline=None, + os.path.join(self.job["workdir"], "control"), "r", newline=None ) as inp: + tmp = inp.readlines() + with open( + os.path.join(self.job["workdir"], "control"), "w", newline=None + ) as out: + for line in tmp: + if "gridsize" in line: + out.write(f" gridsize {'m5'} \n") + if "$disp" in line and self.job["func"] in ("kt2", "kt1"): + pass + else: + out.write(line + "\n") + time.sleep(0.02) + + with open(outputpath, "w", newline=None) as outputfile: + subprocess.call( + [external_paths["mpshiftpath"], "-smpcpus", str(self.job["omp"])], + shell=False, + stdin=None, + stderr=subprocess.STDOUT, + universal_newlines=False, + cwd=self.job["workdir"], + stdout=outputfile, + env=ENVIRON, + ) + time.sleep(0.02) + # check if shift calculation is converged: + if os.path.isfile(outputpath): + with open(outputpath, "r", encoding=CODING, newline=None) as inp: stor = inp.readlines() found = False for line in stor: @@ -1212,6 +1220,9 @@ def _nmr_shielding(self): f"{last_folders(self.job['workdir'], 2):18}" ) self.job["success"] = False + else: + print(f"The file mpshift.out could not be found!") + self.job["success"] = False def _genericoutput(self): """ @@ -1244,7 +1255,9 @@ def _genericoutput(self): ) self.job["success"] = False except ValueError: - print(f"{'ERROR:':{WARNLEN}}ValueError in generic_output, nmrprop.dat can be flawed !") + print( + f"{'ERROR:':{WARNLEN}}ValueError in generic_output, nmrprop.dat can be flawed !" + ) self.job["success"] = False self.job["success"] = True fnamecoupl = "escf.out" @@ -1278,7 +1291,9 @@ def _genericoutput(self): ) self.job["success"] = False except ValueError: - print(f"{'ERROR:':{WARNLEN}}ValueError in generic_output, nmrprop.dat can be flawed") + print( + f"{'ERROR:':{WARNLEN}}ValueError in generic_output, nmrprop.dat can be flawed" + ) self.job["success"] = False self.job["success"] = True with open( @@ -1304,6 +1319,7 @@ def _optrot(self, silent=False): calculate optical rotation """ if not self.job["onlyread"]: + # update control file with open( os.path.join(self.job["workdir"], "control"), "r", newline=None ) as inp: @@ -1331,7 +1347,7 @@ def _optrot(self, silent=False): os.path.join(self.job["workdir"], "escf.out"), "w", newline=None ) as outputfile: subprocess.call( - [self.job["progpath"]], + [external_paths["escfpath"]], shell=False, stdin=None, stderr=subprocess.STDOUT, @@ -1340,7 +1356,8 @@ def _optrot(self, silent=False): stdout=outputfile, env=ENVIRON, ) - time.sleep(0.02) + time.sleep(0.02) + # Read output # check if scf is converged: if os.path.isfile(os.path.join(self.job["workdir"], "escf.out")): with open( @@ -1356,18 +1373,7 @@ def _optrot(self, silent=False): velocity_rep = False # (velocity representation) do_read_length = False do_read_velocity = True - hybriddfa = ( - "pbe0", - "pw6b95", - "wb97x-d3", - "cam-b3lyp", - "b3-lyp", - "pbeh-3c", - "m06x", - "bh-lyp", - "tpssh", - ) - if self.job["func2"] in hybriddfa: + if self.job["func2"] in hybrid_dfa: do_read_length = True dum = 0 frequencies = self.job["freq_or"] @@ -1431,11 +1437,17 @@ def execute(self): elif self.job["jobtype"] == "xtb_sp": self._xtb_sp() elif self.job["jobtype"] == "prep": - self._prep_cefine() + if not self.job["onlyread"]: + self._prep_cefine() + else: + self.job["success"] = True elif self.job["jobtype"] in ("sp", "sp_implicit"): if self.job["prepinfo"]: # do cefine first - self._prep_cefine() + if not self.job["onlyread"]: + self._prep_cefine() + else: + self.job["success"] = True if not self.job["success"]: return self._sp() @@ -1447,7 +1459,10 @@ def execute(self): self._genericoutput() elif self.job["jobtype"] in ("couplings", "couplings_sp"): if self.job["prepinfo"]: - self._prep_cefine() + if not self.job["onlyread"]: + self._prep_cefine() + else: + self.job["success"] = True if not self.job["success"]: return if self.job["jobtype"] == "couplings_sp": @@ -1464,8 +1479,10 @@ def execute(self): self._nmr_coupling() elif self.job["jobtype"] in ("shieldings", "shieldings_sp"): if self.job["prepinfo"]: - self._prep_cefine() - # print('performed cefine') + if not self.job["onlyread"]: + self._prep_cefine() + else: + self.job["success"] = True if not self.job["success"]: return if self.job["copymos"]: @@ -1489,7 +1506,10 @@ def execute(self): if self.job["prepinfo"]: tmp_solvent = self.job["solvent"] self.job["solvent"] = "gas" - self._prep_cefine() + if not self.job["onlyread"]: + self._prep_cefine() + else: + self.job["success"] = True if not self.job["success"]: return self._sp() @@ -1499,7 +1519,10 @@ def execute(self): self._xtb_gsolv() elif self.job["jobtype"] in ("opt-rot", "opt-rot_sp"): if self.job["prepinfo"]: - self._prep_cefine() + if not self.job["onlyread"]: + self._prep_cefine() + else: + self.job["success"] = True if not self.job["success"]: return if self.job["jobtype"] == "opt-rot_sp": diff --git a/censo_qm/tutorial.py b/censo_qm/tutorial.py index 1fa259a..6ec1f8e 100644 --- a/censo_qm/tutorial.py +++ b/censo_qm/tutorial.py @@ -1,14 +1,15 @@ """ Guide for setting up CENSO and performing a CENSO calculation """ -from .cfg import DESCR +from .cfg import DESCR, censo_solvent_db from .inputhandling import internal_settings + def interactiv_doc(): """Guide for interactive explaination of CENSO""" default_settings = internal_settings() - + general = """ Commandline Energetic Sorting (CENSO) is a sorting algorithm for efficient evaluation of Structure Ensembles (SE). The input ensemble (or single structure) @@ -613,7 +614,16 @@ def interactiv_doc(): * ALPB_Gsolv [xtb] * GBSA_Gsolv [xtb] - """ + Available solvents and naming convention employed in CENSO: + +""" + tmp = " " + for count, solvent in enumerate(list(censo_solvent_db.keys())): + tmp += f"{solvent}, " + if (count + 1) % 5 == 0: + solvation += tmp + "\n" + tmp = " " + solvation += tmp example_applications = """ CENSO can be used for several applications / target quantities. Some are @@ -661,6 +671,12 @@ def interactiv_doc(): Functionals that can be employed by TM for funcS {', '.join(internal_settings.func_s_tm)} + Functionals that can be employed by TM for funcOR + {', '.join(internal_settings.func_or_tm)} + + Functionals that can be employed by TM for funcOR_SCF + {', '.join(internal_settings.func_or_scf_tm)} + Functionals that can be employed by ORCA for func: {', '.join(internal_settings.func_orca)} @@ -674,20 +690,16 @@ def interactiv_doc(): {', '.join(internal_settings.func_s_orca)} """ -# CONFORMER numbering kept from crest input - -# folders which are created - #part0_sp - #GFN_unbiased/ - #rrho_part1/ - #b97-3c/ folder of func name - + # CONFORMER numbering kept from crest input + # folders which are created + # part0_sp + # GFN_unbiased/ + # rrho_part1/ + # b97-3c/ folder of func name # which functionals are available - - tutorial_data = { "general": general, "censorc": censorc, @@ -699,7 +711,7 @@ def interactiv_doc(): "functionals": functionals, "jobscript": jobscript, } - tutorial_data['everything'] = '\n\n'.join(tutorial_data.values()) + tutorial_data["everything"] = "\n\n".join(tutorial_data.values()) options = list(tutorial_data.keys()) print(DESCR) @@ -711,7 +723,7 @@ def interactiv_doc(): print("\nPlease input your information request:") while True: user_input = input() - if user_input.strip() in ('q', 'exit'): + if user_input.strip() in ("q", "exit"): break if user_input.strip() not in options: print(f"Options are: {options}") @@ -722,7 +734,7 @@ def interactiv_doc(): # print(everything) # if user_input == 'jobscript': # print(jobscript) - #break + # break print("\n\nDo you want information on any of the other topics?") print("\nIf you want to exit please type one of the following: exit or q") print("\n****CENSO TUTORIAL END****") diff --git a/censo_qm/utilities.py b/censo_qm/utilities.py index d7b20f9..fbfa5c8 100755 --- a/censo_qm/utilities.py +++ b/censo_qm/utilities.py @@ -106,7 +106,7 @@ def t2x(path, writexyz=False, outfile="original.xyz"): if not os.path.basename(path) == "coord": if os.path.isfile(path): with open(path, "r", encoding=CODING, newline=None) as f: - if not '$coord' in f.readline(): + if not "$coord" in f.readline(): path = os.path.join(path, "coord") else: path = os.path.join(path, "coord") @@ -172,7 +172,17 @@ def x2t(path, infile="inp.xyz"): def write_trj( - results, cwd, outpath, optfolder, nat, attribute, temperature, consider_sym, overwrite=False, *args, **kwargs + results, + cwd, + outpath, + optfolder, + nat, + attribute, + temperature, + consider_sym, + overwrite=False, + *args, + **kwargs, ): """ Write trajectory (multiple xyz geometries) to file. @@ -199,13 +209,13 @@ def write_trj( ### coordinates in xyz out.write(" {}\n".format(nat)) xtbfree = conf.calc_free_energy( - e=energy, + e=energy, solv=None, rrho=rrho, out=True, t=temperature, - consider_sym=consider_sym - ) + consider_sym=consider_sym, + ) if xtbfree is not None: xtbfree = f"{xtbfree:20.8f}" out.write( @@ -216,9 +226,10 @@ def write_trj( for line in conf_xyz: out.write(line + "\n") except (FileExistsError, ValueError): - print(f"{'WARNING:':{WARNLEN}}Could not write trajectory: " - f"{last_folders(outpath, 1)}." - ) + print( + f"{'WARNING:':{WARNLEN}}Could not write trajectory: " + f"{last_folders(outpath, 1)}." + ) def check_for_float(line): @@ -396,7 +407,9 @@ def calc_boltzmannweights(confs, property, T): T = float(T) except ValueError: T = 298.15 # K - print(f"{'WARNING:':{WARNLEN}}Temperature can not be converted and is therfore set to T = {T} K.") + print( + f"{'WARNING:':{WARNLEN}}Temperature can not be converted and is therfore set to T = {T} K." + ) if T == 0: T += 0.00001 # avoid division by zero try: @@ -436,8 +449,12 @@ def new_folders(cwd, conflist, foldername, save_errors, store_confs, silent=Fals print(e) if not os.path.isdir(tmp_dir): print(f"{'ERROR:':{WARNLEN}}Could not create folder for CONF{conf.id}!") - print(f"{'ERROR:':{WARNLEN}}CONF{conf.id} is removed, because IO failed!") - save_errors.append(f"{'ERROR:':{WARNLEN}}CONF{conf.id} was removed, because IO failed!") + print( + f"{'ERROR:':{WARNLEN}}CONF{conf.id} is removed, because IO failed!" + ) + save_errors.append( + f"{'ERROR:':{WARNLEN}}CONF{conf.id} was removed, because IO failed!" + ) store_confs.append(conflist.pop(conflist.index(conf))) if not silent: print("Constructed folders!") @@ -752,7 +769,9 @@ def crest_routine(config, conformers, func, store_confs, prev_calculated=None): store = inp.readlines() except (Exception) as e: print(f"{'ERROR:':{WARNLEN}}{e}") - print(f"{'ERROR:':{WARNLEN}}output file (enso.tags) of CREST routine does not exist!") + print( + f"{'ERROR:':{WARNLEN}}output file (enso.tags) of CREST routine does not exist!" + ) keep = [] if config.crestcheck: try: @@ -788,9 +807,9 @@ def format_line(key, value, options, optionlength=70, dist_to_options=30): """ # limit printout of possibilities if len(options) > 0: - separator = '#' + separator = "#" else: - separator = '' + separator = "" if len(str(options)) > optionlength: length = 0 reduced = [] @@ -822,7 +841,8 @@ def check_tasks(results, check=False, thresh=0.25): sys.exit(1) if fail_rate >= thresh and check: print( - f"{'ERROR:':{WARNLEN}}{fail_rate*100} % of the calculations failed!" "\nGoing to exit!" + f"{'ERROR:':{WARNLEN}}{fail_rate*100} % of the calculations failed!" + "\nGoing to exit!" ) sys.exit(1) elif fail_rate >= thresh: @@ -874,10 +894,50 @@ def calc_weighted_std_dev(data, weights=None): std_dev = math.sqrt(variance) return std_dev + def print_errors(line, save_errors): """print line and append to list save_errors""" print(line) try: save_errors.append(line) except Exception: - pass \ No newline at end of file + pass + + +def conf_in_interval(conformers, full_free_energy=True, bm=True): + """ number of conformers (and Boltzmann weights) within free energy window intervals """ + basins = {} + max_rel_free = max([getattr(conf, "rel_free_energy") for conf in conformers]) + for i in frange(0.5, 6.0, 0.5): + if i <= max_rel_free + 0.5: + basins[i] = {"nconf": 0, "bmweight": 0.0} + + for conf in conformers: + for ewin in basins.keys(): + if getattr(conf, "rel_free_energy") <= float(ewin): + basins[ewin]["nconf"] += 1 + if bm: + basins[ewin]["bmweight"] += getattr(conf, "bm_weight") + out = [] + if full_free_energy: + outg = "G" + else: + outg = "g" + out.append(f"\nNumber of conformers observed within the following Δ{outg} windows:") + if bm: + out.append(f"Δ{outg} [kcal/mol] #CONF sum(Boltzmann_weights)") + else: + out.append(f"Δ{outg} [kcal/mol] #CONF") + out.append("".ljust(int(45), "-")) + for ewin in basins.keys(): + if bm: + tmp = f"0 - {ewin}" + out.append( + f"{tmp:^{14}} {basins[ewin]['nconf']} {basins[ewin]['bmweight']: .2f}" + ) + else: + tmp = f"0 - {ewin}" + out.append(f"{tmp:^{14}} {basins[ewin]['nconf']}") + out.append("".ljust(int(45), "-")) + for line in out: + print(line) diff --git a/setup.cfg b/setup.cfg index 70979fb..40e7742 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = censo-QM -version = 1.0.6 +version = 1.0.9 description = CENSO - Comandline ENergetic SOrting for conformer rotamer ensembles long_description = file: README.rst long_description_content_type = text/x-rst