From 486f77d580bcab00ce6a3830875702b74fc79dc8 Mon Sep 17 00:00:00 2001 From: Fabian Bohle <49951809+fabothch@users.noreply.github.com> Date: Thu, 25 Mar 2021 15:36:55 +0100 Subject: [PATCH] symmetry in G_mRRHO (#10) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit • better time printout, which also includes restarts • symmetry in the G_mRRHO contribution is always detected, but its contribution to S_rot can be switched on and off by CENSO -consider_sym without the need to recalculate • if new censorc remote configuration files are created and an existing .censorc is found the program paths can be copied to the new censorc • printout of OR data to file • printout of NMR data to file • new option to creat a censo.inp file from all settings set by command line --- README.rst | 2 +- censo_qm/censo.py | 74 ++++++++++--- censo_qm/cfg.py | 39 ++++++- censo_qm/cheapscreening.py | 13 ++- censo_qm/datastructure.py | 180 +++++++++++++++++++++++-------- censo_qm/ensembledata.py | 13 +++ censo_qm/inputhandling.py | 199 ++++++++++++++++++++++++---------- censo_qm/nmrproperties.py | 104 ++++++++++++------ censo_qm/opticalrotation.py | 66 ++++++++++-- censo_qm/optimization.py | 112 +++++++++++++++---- censo_qm/prescreening.py | 207 +++++++++++++++++++++++++++--------- censo_qm/qm_job.py | 155 +++++++++++++++++---------- censo_qm/refinement.py | 109 +++++++++++++++---- censo_qm/setupcenso.py | 54 +++++++++- censo_qm/tm_job.py | 1 - censo_qm/utilities.py | 20 +++- 16 files changed, 1033 insertions(+), 315 deletions(-) diff --git a/README.rst b/README.rst index 08c16ae..248c21a 100644 --- a/README.rst +++ b/README.rst @@ -203,7 +203,7 @@ a replacement solvent. E.g. if CCl4 is not available choose CHCl3. .. raw:: html

- censo_solvents.json + censo_solvents.json

The solvent file is directly used in `CENSO` and typos will cause calculations to crash! diff --git a/censo_qm/censo.py b/censo_qm/censo.py index ba200e2..5b7b3fa 100644 --- a/censo_qm/censo.py +++ b/censo_qm/censo.py @@ -62,6 +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"] print(f"Ran part0 in {ensembledata.part_info['part0']:0.4f} seconds") # RUNNING PART1 @@ -85,6 +86,8 @@ def main(argv=None): sys.exit(1) 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"] print(f"Ran part1 in {ensembledata.part_info['part1']:0.4f} seconds") # RUNNING PART2 @@ -108,6 +111,8 @@ def main(argv=None): sys.exit(1) 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"] print(f"Ran part2 in {ensembledata.part_info['part2']:0.4f} seconds") # RUNNING PART3 @@ -131,6 +136,7 @@ def main(argv=None): sys.exit(1) toc = perf_counter() ensembledata.part_info["part3"] = toc - tic + ensembledata.previous_part_info["part3"] += ensembledata.part_info["part3"] print(f"Ran part3 in {ensembledata.part_info['part3']:0.4f} seconds") # RUNNING PART4 @@ -154,6 +160,7 @@ def main(argv=None): sys.exit(1) toc = perf_counter() ensembledata.part_info["part4"] = toc - tic + ensembledata.previous_part_info["part4"] += ensembledata.part_info["part4"] print(f"Ran part4 in {ensembledata.part_info['part4']:0.4f} seconds") # RUNNING PART5 @@ -177,6 +184,7 @@ def main(argv=None): sys.exit(1) toc = perf_counter() ensembledata.part_info["part5"] = toc - tic + ensembledata.previous_part_info["part5"] += ensembledata.part_info["part5"] print(f"Ran part5 in {ensembledata.part_info['part5']:0.4f} seconds") # save current data to jsonfile @@ -190,51 +198,91 @@ def main(argv=None): # END of CENSO timings = 0.0 + prev_timings = 0.0 if len(str(config.nconf)) > 5: conflength = len(str(config.nconf)) else: conflength = 5 - print(f"\n\n{'Part':20}: {'#conf':>{conflength}} time") - print("".ljust(int(PLENGTH / 2), "-")) - print(f"{'Input':20}: {ensembledata.nconfs_per_part['starting']:{conflength}} -") + try: + tmp = [] + tmp_prev = [] + if config.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']) + 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']) + if config.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']) + 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 + 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}}") if config.part0: print( - f"{'Part0_all':20}: {ensembledata.nconfs_per_part['part0']:{conflength}} {ensembledata.part_info['part0']:.2f} s" + 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'] if config.part1: print( - f"{'Part1_initial_sort':20}: {ensembledata.nconfs_per_part['part1_firstsort']:{conflength}} -" + 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" ) print( - f"{'Part1_all':20}: {ensembledata.nconfs_per_part['part1_firstsort']:{conflength}} {ensembledata.part_info['part1']:.2f} s" + 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'] if config.part2: print( - f"{'Part2_opt':20}: {ensembledata.nconfs_per_part['part2_opt']:{conflength}} -" + 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" ) print( - f"{'Part2_all':20}: {ensembledata.nconfs_per_part['part2']:{conflength}} {ensembledata.part_info['part2']:.2f} s" + 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'] if config.part3: print( - f"{'Part3_all':20}: {ensembledata.nconfs_per_part['part3']:{conflength}} {ensembledata.part_info['part3']:.2f} s" + 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'] if config.part4: print( - f"{'Part4':20}: {'':{conflength}} {ensembledata.part_info['part4']:.2f} s" + 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'] if config.optical_rotation: print( - f"{'Part5':20}: {'':{conflength}} {ensembledata.part_info['part5']:.2f} s" + 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"] - print("".ljust(int(PLENGTH / 2), "-")) - print(f"{'All parts':20}: {'':{conflength}} {timings:.2f} s") + 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("\nCENSO all done!") return 0 diff --git a/censo_qm/cfg.py b/censo_qm/cfg.py index 5414460..69023f9 100644 --- a/censo_qm/cfg.py +++ b/censo_qm/cfg.py @@ -5,7 +5,7 @@ """ import os -__version__ = "1.0.7" +__version__ = "1.0.8" DESCR = f""" ______________________________________________________________ @@ -37,6 +37,7 @@ PLENGTH = 100 AU2J = 4.3597482e-18 # a.u.(hartree/mol) to J KB = 1.3806485279e-23 # J/K +R = 1.987203585e-03 # kcal/(mol*K) AU2KCAL = 627.50947428 BOHR2ANG = 0.52917721067 WARNLEN = max([len(i) for i in ['WARNING:', 'ERROR:', 'INFORMATION:']])+1 @@ -2654,3 +2655,39 @@ def dict_to_NMRRef(self, dictionary): 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 +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 +} diff --git a/censo_qm/cheapscreening.py b/censo_qm/cheapscreening.py index d809bb4..09ed79c 100644 --- a/censo_qm/cheapscreening.py +++ b/censo_qm/cheapscreening.py @@ -259,6 +259,7 @@ 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["info"] = "calculated" conf.cheap_prescreening_gsolv_info["method"] = conf.job["method2"] conf.cheap_prescreening_gsolv_info["gas-energy"] = conf.job[ @@ -345,7 +346,13 @@ def part0(config, conformers, ensembledata): rrho = None energy = "cheap_prescreening_sp_info" for conf in calculate: - conf.calc_free_energy(e=energy, solv=solvation, rrho=rrho) + conf.calc_free_energy( + e=energy, + solv=solvation, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) try: minfree = min([i.free_energy for i in calculate if i is not None]) except ValueError: @@ -428,8 +435,8 @@ def part0(config, conformers, ensembledata): (5, 2), (5, 2), ] - columndescription[1] = f"{config.part0_gfnv.upper()}-xTB[{config.sm_rrho.upper()}]" - columndescription[2] = f"{config.part0_gfnv.upper()}-xTB[{config.sm_rrho.upper()}]" + columndescription[1] = f"{config.part0_gfnv.upper()}-xTB[{config.sm_rrho}]" + columndescription[2] = f"{config.part0_gfnv.upper()}-xTB[{config.sm_rrho}]" columndescription[3] = instruction["method"] columndescription[4] = instruction["method2"] if config.solvent == "gas": diff --git a/censo_qm/datastructure.py b/censo_qm/datastructure.py index 3c7ecc4..3478579 100755 --- a/censo_qm/datastructure.py +++ b/censo_qm/datastructure.py @@ -4,6 +4,8 @@ """ from collections import OrderedDict from .utilities import print +from .cfg import R, AU2KCAL +from math import log class MoleculeData: @@ -23,6 +25,8 @@ def __init__( rel_xtb_energy=None, rel_xtb_free_energy=None, sym="c1", + linear=False, + symnum=1, gi=1.0, removed=False, free_energy=0.0, @@ -37,6 +41,7 @@ def __init__( "energy": None, "gas-energy": None, "solv-energy": None, + "range": None, "info": "not_calculated", "method": None, "prev_methods": None, @@ -61,6 +66,7 @@ def __init__( }, prescreening_grrho_info={ "energy": None, + "range": None, "info": "not_calculated", "method": None, "fuzzythr": 0.0, @@ -101,6 +107,7 @@ def __init__( }, prescreening_gsolv_info={ "energy": None, + "range": None, "gas-energy": None, "info": "not_calculated", "method": None, @@ -231,8 +238,8 @@ def __init__( temperature_info["range"] = [] if xtb_energy is None: xtb_energy = 100.0 - if xtb_energy_unbiased is None: - xtb_energy_unbiased = 100.0 + #if xtb_energy_unbiased is None: + # xtb_energy_unbiased = 100.0 if xtb_free_energy is None: xtb_free_energy = 100.0 if rel_xtb_energy is None: @@ -243,6 +250,7 @@ def __init__( cheap_prescreening_sp_info["energy"] = 0.0 if cheap_prescreening_gsolv_info.get("energy") is None: cheap_prescreening_gsolv_info["energy"] = 0.0 + self._initialize(cheap_prescreening_gsolv_info) if prescreening_sp_info.get("energy") is None: prescreening_sp_info["energy"] = 0.0 if lowlevel_sp_info.get("energy") is None: @@ -251,6 +259,7 @@ def __init__( highlevel_sp_info["energy"] = 0.0 if prescreening_grrho_info.get("energy") is None: prescreening_grrho_info["energy"] = 0.0 + self._initialize(prescreening_grrho_info) if lowlevel_grrho_info.get("energy") is None: lowlevel_grrho_info["energy"] = 0.0 self._initialize(lowlevel_grrho_info) @@ -259,6 +268,7 @@ def __init__( self._initialize(lowlevel_hrrho_info) if prescreening_gsolv_info.get("energy") is None: prescreening_gsolv_info["energy"] = 0.0 + self._initialize(prescreening_gsolv_info) if lowlevel_gsolv_info.get("energy") is None: lowlevel_gsolv_info["energy"] = 0.0 self._initialize(lowlevel_gsolv_info) @@ -313,6 +323,8 @@ def __init__( ) if not isinstance(sym, str): raise TypeError("Please provide symmetry as string.") + if not isinstance(linear, bool): + raise TypeError("Please provide information on linear molecules as bool.") if not isinstance(gi, float): try: gi = float(gi) @@ -367,6 +379,8 @@ def __init__( self.rel_xtb_energy = rel_xtb_energy self.rel_xtb_free_energy = rel_xtb_free_energy self.sym = sym + self.linear = linear + self.symnum = symnum self.gi = gi self.cheap_prescreening_gsolv_info = cheap_prescreening_gsolv_info self.cheap_prescreening_sp_info = cheap_prescreening_sp_info @@ -423,6 +437,8 @@ def reset_range_info(self, trange=None): attributes = [ "lowlevel_grrho_info", "lowlevel_hrrho_info", + "highlevel_grrho_info", + "highlevel_hrrho_info", "lowlevel_gsolv_info", "highlevel_gsolv_info", ] @@ -491,7 +507,59 @@ def provide_runinfo(self): runinfo.append((key, getattr(self, key))) return OrderedDict(runinfo) - def calc_free_energy(self, e=None, solv=None, rrho=None, t=None, out=False): + # def calc_free_energy(self, e=None, solv=None, rrho=None, t=None, out=False): + # """ + # 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 t is None: + # try: + # f = 0.0 + # if e is not None: + # if e in ("xtb_energy", "xtb_energy_unbiased"): + # f += getattr(self, e, 0.0) + # else: + # f += getattr(self, e, {"energy": 0.0})["energy"] + # if solv is not None: + # f += getattr(self, solv, {"energy": 0.0})["energy"] + # if rrho is not None: + # f += getattr(self, rrho, {"energy": 0.0})["energy"] + # if not out: + # self.free_energy = f + # else: + # return f + # except Exception as error: + # print("ERROR in _calc_free_energy: ", error) + # if not out: + # self.free_energy = None + # else: + # return f + # else: + # try: + # f = 0.0 + # if e is not None: + # if e in ("xtb_energy", "xtb_energy_unbiased"): + # f += getattr(self, e, 0.0) + # else: + # f += getattr(self, e)["energy"] + # if solv is not None: + # f += getattr(self, solv)["range"].get(t, 0.0) + # if rrho is not None: + # f += getattr(self, rrho)["range"].get(t, 0.0) + # if not out: + # self.free_energy = f + # else: + # return f + # except (Exception, KeyError) as error: + # print("ERROR in _calc_free_energy: ", error) + # if not out: + # self.free_energy = None + # else: + # return f + + 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. @@ -499,46 +567,72 @@ def calc_free_energy(self, e=None, solv=None, rrho=None, t=None, out=False): if out=True free energy is simply returned """ if t is None: - try: - f = 0.0 - if e is not None: - if e in ("xtb_energy", "xtb_energy_unbiased"): - f += getattr(self, e, 0.0) - else: - f += getattr(self, e, {"energy": 0.0})["energy"] - if solv is not None: - f += getattr(self, solv, {"energy": 0.0})["energy"] - if rrho is not None: - f += getattr(self, rrho, {"energy": 0.0})["energy"] - if not out: - self.free_energy = f + print("xxxxxxxxxxxxx No temperature provided in calc_free_energy xxxxxxxxxxxxxxx") + try: + f = 0.0 + if e is not None: + if e in ("xtb_energy", "xtb_energy_unbiased"): + f += getattr(self, e, 0.0) else: - return f - except Exception as error: - print("ERROR in _calc_free_energy: ", error) - if not out: - self.free_energy = None + 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"]) else: - return f - else: - try: - f = 0.0 - if e is not None: - if e in ("xtb_energy", "xtb_energy_unbiased"): - f += getattr(self, e, 0.0) - else: - f += getattr(self, e)["energy"] - if solv is not None: f += getattr(self, solv)["range"].get(t, 0.0) - if rrho is not None: - f += getattr(self, rrho)["range"].get(t, 0.0) - if not out: - self.free_energy = f - else: - return f - except (Exception, KeyError) as error: - print("ERROR in _calc_free_energy: ", error) - if not out: - self.free_energy = None - else: - return f + if rrho is not None: + f += self.get_mrrho( + t, + rrho=rrho, + consider_sym=consider_sym, + symnum=self.symnum, + ) + if not out: + self.free_energy = f + else: + return f + except (Exception, KeyError) as error: + print("ERROR in _calc_free_energy: ", error) + if not out: + self.free_energy = None + else: + return f + + def calc_entropy_sym(self, temperature, symnum=None): + """ RTln(sigma) rotational entropy""" + if symnum is None: + 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 + ): + """ return mRRHO with or without RTln(sigma) (rot entropy)""" + f = 0.0 + if rrho is not None: + if rrho == "direct_input": + f = direct_input + elif rrho in ("prescreening_grrho_info",): + f += ( + getattr(self, rrho) + .get("range", {}) + .get(temperature, getattr(self, rrho, {"energy": 0.0})["energy"]) + ) + elif rrho in ("rrho_optimization",): + 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: + consider_sym = True + if not consider_sym: + # sym is considered but consider_sym off + if symnum is None: + symnum = self.symnum + f += -self.calc_entropy_sym(temperature, symnum=symnum) + return f + else: + return f diff --git a/censo_qm/ensembledata.py b/censo_qm/ensembledata.py index a53df7a..9c8850f 100644 --- a/censo_qm/ensembledata.py +++ b/censo_qm/ensembledata.py @@ -10,6 +10,18 @@ def __init__( "part2_opt": None, "part2": None, "part3": None, + "part4": None, + "part5": None, + }, + previous_part_info={ + "part0": 0.0, + "part1_firstsort": 0.0, + "part1": 0.0, + "part2_opt": 0.0, + "part2": 0.0, + "part3": 0.0, + "part4": 0.0, + "part5": 0.0, }, avGcorrection=None, comment=None, @@ -42,6 +54,7 @@ def __init__( self.id = id self.filename = filename self.part_info = part_info + self.previous_part_info = previous_part_info self.avGcorrection = avGcorrection self.comment = comment self.bestconf = bestconf diff --git a/censo_qm/inputhandling.py b/censo_qm/inputhandling.py index 490ec2f..eeee45a 100755 --- a/censo_qm/inputhandling.py +++ b/censo_qm/inputhandling.py @@ -118,18 +118,18 @@ def cml(startup_description, options, argv=None): action="store", required=False, metavar="", - help="Applies structure constraint to input/DFT geometry for mRRHO calcuation." + help="Uses SPH and applies structure constraint to input/DFT geometry for mRRHO calcuation. " "Options are: ['on' or 'off'].", ) group1.add_argument( "-consider_sym", - "---consider_sym", - dest="-consider_sym", + "--consider_sym", + dest="consider_sym", choices=["on", "off"], action="store", required=False, metavar="", - help="Consider symmetry in mRRHO calcuation (based on desy xtb threshold)." + help="Consider symmetry in mRRHO calcuation (based on desy xtb threshold). " "Options are: ['on' or 'off'].", ) group1.add_argument( @@ -140,7 +140,7 @@ def cml(startup_description, options, argv=None): action="store", required=False, metavar="", - help="Applies constraint to rmsdpot.xyz to be consistent to CREST." + help="Applies constraint to rmsdpot.xyz to be consistent to CREST. " "Options are: ['on' or 'off'].", ) group1.add_argument( @@ -927,6 +927,16 @@ def cml(startup_description, options, argv=None): help="Write new configuration file, which is placed into the current " "directory.", ) + group8.add_argument( + "-copyinput", + "--copyinput", + dest="copyinput", + default=False, + action="store_true", + required=False, + help="Write all current settings to a censo.inp configuration file, which is placed into the current " + "directory.", + ) group8.add_argument( "-progress", "--progress", @@ -1356,6 +1366,7 @@ class internal_settings: "dmf", "dmso", "ether", + "ethanol", "ethylacetate", "furane", "hexadecane", @@ -1502,7 +1513,7 @@ def __init__(self): ("trange", {"default": [273.15, 378.15, 5], "type": list}), ("multitemp", {"default": True, "type": bool}), ("evaluate_rrho", {"default": True, "type": bool}), - ("consider_sym", {"default": False, "type": bool}), + ("consider_sym", {"default": True, "type": bool}), ("bhess", {"default": True, "type": bool}), ("imagthr", {"default": "automatic", "type": str}), ("sthr", {"default": "automatic", "type": str}), @@ -1739,7 +1750,6 @@ def __init__(self): "sm2", "nat", "radsize", - "consider_sym", "cosmorsparam", ] # might be changed, but data may be lost/overwritten @@ -1760,12 +1770,13 @@ def __init__(self): "freq_or": False, "func3":False, "basis3":False, - # "consider_sym": False, # --> reset all rrho values! - # funcJ - # basisJ - # funcS - # basisS - # sm4_ j s + "func_j":False, + "basis_j":False, + "sm4_j":False, + "func_s":False, + "basis_s":False, + "sm4_s":False, + # "consider_sym": calculated on the fly } @@ -2791,7 +2802,8 @@ def print_parameters(self): ) info.append(["justprint", "".ljust(int(PLENGTH / 2), "-")]) info.append(["part1", "part1"]) - info.append(["nconf", "starting number of considered conformers"]) + if not self.part0: + info.append(["nconf", "starting number of considered conformers"]) info.append(["prog", "program for part1"]) info.append(["func", "functional for initial evaluation"]) info.append(["basis", "basis set for initial evaluation"]) @@ -3450,7 +3462,7 @@ def processQMpaths(self, requirements, error_logical): external_paths.update(self.external_paths) return error_logical - def write_rcfile(self, pathtofile): + def write_rcfile(self, pathtofile, usepaths=False): """ write new global configruation file into the current directroy. """ @@ -3459,20 +3471,33 @@ def write_rcfile(self, pathtofile): outdata.write("$CENSO global configuration file: .censorc\n") outdata.write(f"$VERSION:{__version__} \n") outdata.write("\n") - outdata.write("ORCA: /path/excluding/binary/\n") - outdata.write("ORCA version: 4.2.1\n") - outdata.write("GFN-xTB: /path/including/binary/xtb-binary\n") - outdata.write("CREST: /path/including/binary/crest-binary\n") - outdata.write("mpshift: /path/including/binary/mpshift-binary\n") - outdata.write("escf: /path/including/binary/escf-binary\n") - outdata.write("\n") - outdata.write("#COSMO-RS\n") - outdata.write( - "ctd = BP_TZVP_C30_1601.ctd cdir = " - '"/software/cluster/COSMOthermX16/COSMOtherm/CTDATA-FILES" ldir = ' - '"/software/cluster/COSMOthermX16/COSMOtherm/CTDATA-FILES"\n' - ) - # outdata.write("cosmothermversion: 16\n") + if usepaths: + # write stored program paths to file + outdata.write(f"ORCA: {self.external_paths['orcapath']}\n") + outdata.write(f"ORCA version: {self.external_paths['orcaversion']}\n") + outdata.write(f"GFN-xTB: {self.external_paths['xtbpath']}\n") + outdata.write(f"CREST: {self.external_paths['crestpath']}\n") + outdata.write(f"mpshift: {self.external_paths['mpshiftpath']}\n") + outdata.write(f"escf: {self.external_paths['escfpath']}\n") + outdata.write("\n") + outdata.write("#COSMO-RS\n") + outdata.write(f"{self.external_paths['cosmorssetup']}\n") + # outdata.write("cosmothermversion: 16\n") + else: + outdata.write("ORCA: /path/excluding/binary/\n") + outdata.write("ORCA version: 4.2.1\n") + outdata.write("GFN-xTB: /path/including/binary/xtb-binary\n") + outdata.write("CREST: /path/including/binary/crest-binary\n") + outdata.write("mpshift: /path/including/binary/mpshift-binary\n") + outdata.write("escf: /path/including/binary/escf-binary\n") + outdata.write("\n") + outdata.write("#COSMO-RS\n") + outdata.write( + "ctd = BP_TZVP_C30_1601.ctd cdir = " + '"/software/cluster/COSMOthermX16/COSMOtherm/CTDATA-FILES" ldir = ' + '"/software/cluster/COSMOthermX16/COSMOtherm/CTDATA-FILES"\n' + ) + # outdata.write("cosmothermversion: 16\n") outdata.write("$ENDPROGRAMS\n\n") outdata.write("$CRE SORTING SETTINGS:\n") outdata.write("$GENERAL SETTINGS:\n") @@ -3546,38 +3571,104 @@ def write_rcfile(self, pathtofile): outdata.write(format_line(key, value, options)) outdata.write("$END CENSORC\n") - def write_enso_inp(self, path=None): + def write_censo_inp(self, path=None): """ - Write file "enso.inp" which is the control file for the calculation + Write file "censo.inp" which stores current settings in the .censorc format """ if path is None: path = self.cwd - with open(os.path.join(path, "enso.inp"), "w", newline=None) as inp: - inp.write("$ File: enso.inp settings of current calculation\n") - data = self.provide_runinfo(extend=False) - for key in data.keys(): - value = data[key] - options = self.value_options.get(key, "possibilities") + args_key = {v: k for k, v in self.key_args_dict.items()} + data = self.provide_runinfo(extend=False) + with open(os.path.join(path, "censo.inp"), "w", newline=None) as outdata: + outdata.write("$File: censo.inp settings of current calculation\n") + outdata.write(f"$VERSION:{__version__} \n") + outdata.write("\n") + # write stored program paths to file + outdata.write(f"ORCA: {self.external_paths['orcapath']}\n") + outdata.write(f"ORCA version: {self.external_paths['orcaversion']}\n") + outdata.write(f"GFN-xTB: {self.external_paths['xtbpath']}\n") + outdata.write(f"CREST: {self.external_paths['crestpath']}\n") + outdata.write(f"mpshift: {self.external_paths['mpshiftpath']}\n") + outdata.write(f"escf: {self.external_paths['escfpath']}\n") + outdata.write("\n") + outdata.write("#COSMO-RS\n") + outdata.write(f"{self.external_paths['cosmorssetup']}\n") + # outdata.write("cosmothermversion: 16\n") + outdata.write("$ENDPROGRAMS\n\n") + outdata.write("$CRE SORTING SETTINGS:\n") + outdata.write("$GENERAL SETTINGS:\n") + for key in OrderedDict(self.defaults_refine_ensemble_general): + value = self._exchange_onoff( + data.get( + key, + OrderedDict(self.defaults_refine_ensemble_general)[key]["default"]), + reverse=True, + ) if key == "nconf" and value is None: value = "all" - value = self._exchange_onoff(value, reverse=True) - # limit printout of possibilities - if len(str(options)) > 80: - length = 0 - reduced = [] - for item in options: - length += len(item) + 2 - if length < 80: - reduced.append(item) - reduced.append("...") - options = reduced - length = 0 - inp.write( - "{}: {:{digits}} # {}\n".format( - str(key), str(value), options, digits=30 - len(str(key)) - ) + key = args_key.get(key, key) + outdata.write(format_line(key, value, "")) + outdata.write("\n$PART0 - CHEAP-PRESCREENING - SETTINGS:\n") + for key in OrderedDict(self.defaults_refine_ensemble_part0): + value = self._exchange_onoff( + data.get( + key, + OrderedDict(self.defaults_refine_ensemble_part0)[key]["default"]), + reverse=True, + ) + key = args_key.get(key, key) + outdata.write(format_line(key, value, "")) + 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"]), + reverse=True, + ) + key = args_key.get(key, key) + outdata.write(format_line(key, value, "")) + 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"]), + 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"]), + reverse=True, + ) + key = args_key.get(key, key) + outdata.write(format_line(key, value, "")) + 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"]), + reverse=True, + ) + key = args_key.get(key, key) + outdata.write(format_line(key, value, "")) + 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"]), + reverse=True, ) - inp.write("$end\n") + 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 3e54fed..61558fd 100644 --- a/censo_qm/nmrproperties.py +++ b/censo_qm/nmrproperties.py @@ -113,7 +113,13 @@ 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, rrho=rrho) + conf.calc_free_energy( + e=energy, + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) conf.free_energy += normalvariate( 0.0, get_std_dev(conf) ) @@ -134,7 +140,13 @@ def average_shieldings(config, calculate, element_ref_shield, energy, solv, rrho if wanted: for _ in range(1000): for conf in calculate: - conf.calc_free_energy(e=energy, solv=solv, rrho=rrho) + conf.calc_free_energy( + e=energy, + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) conf.free_energy += normalvariate( 0.0, (0.4/AU2KCAL) ) @@ -151,34 +163,42 @@ 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]) - - print("# in coord element σ(sigma) SD(σ based on SD Gsolv) SD(σ by 0.4 kcal/mol) shift σ_ref") - print("".ljust(int(100), "-")) - 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) - else "None" - ) - for atom, sigma in averaged.items(): - try: - std_dev = calc_std_dev(sigma_std_dev[atom]) - except Exception as e: - #print(e) - std_dev = 0.0 - try: - std_dev_const = calc_std_dev(sigma_std_dev_const[atom]) - except Exception as e: - #print(e) - std_dev_const = 0.0 - try: - print( - f"{atom:< {10}} {element[atom]:^{7}} {sigma:> {maxsigma}.2f} " - f"{std_dev:^ 24.6f} {std_dev_const:^ 24.6f} {make_shift(atom):>5} {element_ref_shield.get(element[atom], 0.0)}" - ) - except Exception as e: - print(f"{atom:< {10}} {element[atom]:^{7}} {sigma:> {maxsigma}.2f}") - print("".ljust(int(80), "-")) + 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), "-")) + print(line) + 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) + else "None" + ) + for atom, sigma in averaged.items(): + try: + std_dev = calc_std_dev(sigma_std_dev[atom]) + except Exception: + std_dev = 0.0 + try: + std_dev_const = calc_std_dev(sigma_std_dev_const[atom]) + except Exception: + std_dev_const = 0.0 + try: + line = ( + f"{atom:< {10}} {element[atom]:^{7}} {sigma:> {maxsigma}.2f} " + 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") + except Exception: + line=(f"{atom:< {10}} {element[atom]:^{7}} {sigma:> {maxsigma}.2f}") + print(line) + out.write(line+"\n") + line = ("".ljust(int(105), "-")) + print(line) + out.write(line+"\n") def write_anmrrc(config): @@ -648,7 +668,13 @@ def part4(config, conformers, store_confs, ensembledata): for conf in calculate: - conf.calc_free_energy(e=energy, solv=gsolv, rrho=rrho) + conf.calc_free_energy( + e=energy, + solv=gsolv, + rrho=rrho, + t=config.temperature, + 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]) @@ -666,7 +692,10 @@ 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: 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, @@ -1050,7 +1079,13 @@ def part4(config, conformers, store_confs, ensembledata): # write anmr_enso output! print("\nGenerating file anmr_enso for processing with the ANMR program.") for conf in calculate: - conf.calc_free_energy(e=energy, solv=gsolv, rrho=rrho) + conf.calc_free_energy( + e=energy, + solv=gsolv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) try: length = max([str(i.id) for i in calculate]) @@ -1081,7 +1116,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"{getattr(conf, str(rrho), {'energy': 0.0})['energy']: {fmtrrho-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.gi:.3f}\n" ) diff --git a/censo_qm/opticalrotation.py b/censo_qm/opticalrotation.py index 1e51fd8..95a1b45 100644 --- a/censo_qm/opticalrotation.py +++ b/censo_qm/opticalrotation.py @@ -292,7 +292,13 @@ def part5(config, conformers, store_confs, ensembledata): ) for conf in calculate: - conf.calc_free_energy(e=energy, solv=gsolv, rrho=rrho) + conf.calc_free_energy( + e=energy, + solv=gsolv, + rrho=rrho, + t=config.temperature, + 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]) @@ -310,7 +316,12 @@ 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: 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, @@ -540,15 +551,36 @@ def part5(config, conformers, store_confs, ensembledata): f" populated to {conf.bm_weight*100:.2f} % " ) calculate.append(prev_calculated.pop(prev_calculated.index(conf))) + + outlen = config.lenconfx+1 + if len('ENSEMBLE') >= config.lenconfx: + outlen = len('ENSEMBLE')+1 + for freq in config.freq_or: averaged_or = 0.0 - for conf in calculate: - averaged_or += conf.bm_weight * conf.optical_rotation_info["range"].get( - freq, 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") + 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" + ) + # ENSEMBLE + outdata.write( + 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: .3f} in deg*[dm(g/cc)]^(-1)" + f"{averaged_or:> {max_fmt}.3f} in deg*[dm(g/cc)]^(-1)" ) @@ -560,7 +592,13 @@ def part5(config, conformers, store_confs, ensembledata): for _ in range(1000): averaged_or = 0.0 for conf in calculate: - conf.calc_free_energy(e=energy, solv=gsolv, rrho=rrho) + conf.calc_free_energy( + e=energy, + solv=gsolv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) conf.free_energy += normalvariate( 0.0, conf.lowlevel_gsolv_compare_info["std_dev"] ) @@ -590,7 +628,7 @@ def part5(config, conformers, store_confs, ensembledata): print(e) max_fmt = 16 print( - f" SD based on SD of Gsolv (part2) " + f" SD based on SD of Gsolv (part2) " f": {calc_std_dev(all_or):> {max_fmt}.3f} in deg*[dm(g/cc)]^(-1)" ) @@ -600,7 +638,13 @@ def part5(config, conformers, store_confs, ensembledata): for _ in range(1000): averaged_or = 0.0 for conf in calculate: - conf.calc_free_energy(e=energy, solv=gsolv, rrho=rrho) + conf.calc_free_energy( + e=energy, + solv=gsolv, + rrho=rrho, + t=config.consider_sym, + consider_sym=config.consider_sym + ) conf.free_energy += normalvariate( 0.0, (0.4/AU2KCAL) ) @@ -630,7 +674,7 @@ def part5(config, conformers, store_confs, ensembledata): print(e) max_fmt = 16 print( - f" SD based on SD in G of 0.4 kcal/mol " + 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)" ) diff --git a/censo_qm/optimization.py b/censo_qm/optimization.py index 51074bd..4ec5574 100755 --- a/censo_qm/optimization.py +++ b/censo_qm/optimization.py @@ -369,11 +369,11 @@ def part2(config, conformers, store_confs, ensembledata): conf.reset_job_info() # *************************************************************************** # NEW ENSEMBLE OPTIMIZER: + start_opt_time = time.perf_counter() if calculate: 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 @@ -525,12 +525,22 @@ def part2(config, conformers, store_confs, ensembledata): # check if too many calculations failed ### for conf in list(calculate): + # 'rrho_optimization' used in get_mrrho print( f"The G_mRRHO calculation on crudely optimized DFT " f"geometry @ {conf.job['symmetry']} " f"{check[conf.job['success']]} for " f"{last_folders(conf.job['workdir'], 3):>{pl}}: " - f"{conf.job['energy']:>.8f}" + f"{conf.job['energy']:>.8f} " + f"""S_rot(sym)= {conf.calc_entropy_sym( + config.temperature, + symnum=conf.job['symnum']):>.7f} """ + f"""using= {conf.get_mrrho( + config.temperature, + rrho='direct_input', + consider_sym=config.consider_sym, + symnum=conf.job['symnum'], + direct_input=conf.job["energy"]):>.7f}""" ) if not conf.job["success"]: store_confs.append(calculate.pop(calculate.index(conf))) @@ -540,6 +550,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'] for conf in list(calculate): if conf.id in converged_run1: @@ -580,7 +592,7 @@ def part2(config, conformers, store_confs, ensembledata): print("Spearman rank evaluation is performed in the next cycle.") cycle_spearman.append("") run_spearman = False - # if run == 1 and not gesc: + # 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("") @@ -603,18 +615,24 @@ 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 = "axqzv" # to get 0.0 contribution + rrho_energy = None # to get 0.0 contribution for i in range(maxecyc): try: minecyc.append( min( [ conf.job["ecyc"][i] - + getattr(conf, "optimization_info").get( - rrho_energy, 0.0 + + conf.get_mrrho( + config.temperature, + rrho=rrho_energy, + consider_sym=config.consider_sym, ) + #+ getattr(conf, "optimization_info").get( + # rrho_energy, 0.0 + #) for conf in calculate + prev_calculated if conf.job["ecyc"][i] is not None ] @@ -632,9 +650,14 @@ def part2(config, conformers, store_confs, ensembledata): conf.job["decyc"].append( ( conf.job["ecyc"][i] - + getattr(conf, "optimization_info").get( - rrho_energy, 0.0 + + conf.get_mrrho( + config.temperature, + rrho=rrho_energy, + consider_sym=config.consider_sym, ) + #+ getattr(conf, "optimization_info").get( + # rrho_energy, 0.0 + #) - minecyc[i] ) * AU2KCAL @@ -868,6 +891,8 @@ def part2(config, conformers, store_confs, ensembledata): print("sum: {:>.2f}".format(sum(timings))) print("\nCONVERGED optimizations for the following remaining conformers:") prev_calculated.sort(key=lambda x: int(x.id)) + end_opt_time = time.perf_counter() + ensembledata.part_info["part2_opt"] = end_opt_time - start_opt_time # end if calculate-- for conf in list(prev_calculated): print( @@ -1340,13 +1365,24 @@ def part2(config, conformers, store_confs, ensembledata): f"The lowlevel G_mRRHO calculation @ {conf.job['symmetry']} " f"{check[conf.job['success']]} for " f"{last_folders(conf.job['workdir'], 2):>{pl}}: " - f"{conf.job['energy']:>.8f}" + f"{conf.job['energy']:>.8f} " + f"""S_rot(sym)= {conf.calc_entropy_sym( + config.temperature, + symnum=conf.job['symnum']):>.7f} """ + f"""using= {conf.get_mrrho( + config.temperature, + rrho='direct_input', + consider_sym=config.consider_sym, + symnum=conf.job['symnum'], + direct_input=conf.job["energy"]):>.7f}""" ) if not conf.job["success"]: conf.lowlevel_grrho_info["info"] = "failed" store_confs.append(calculate.pop(calculate.index(conf))) else: conf.sym = conf.job["symmetry"] + conf.linear = conf.job["linear"] + conf.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" @@ -1377,7 +1413,12 @@ def part2(config, conformers, store_confs, ensembledata): f"The lowlevel G_mRRHO calculation @ {conf.sym} " f"{check[conf.job['success']]} for " f"{last_folders(conf.job['workdir'], 2):>{pl}}: " - f"{conf.lowlevel_grrho_info['energy']:>.8f}" + f"{conf.lowlevel_grrho_info['energy']:>.8f} " + f"S_rot(sym)= {conf.calc_entropy_sym(config.temperature):>.7f}" + f""" using= {conf.get_mrrho( + config.temperature, + rrho='lowlevel_grrho_info', + consider_sym=config.consider_sym):>.7f}""" ) calculate.append(prev_calculated.pop(prev_calculated.index(conf))) if not calculate: @@ -1395,7 +1436,10 @@ 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 + ), lambda conf: getattr(conf, "free_energy"), lambda conf: getattr(conf, "rel_free_energy"), lambda conf: getattr(conf, "bm_weight") * 100, @@ -1469,7 +1513,13 @@ def part2(config, conformers, store_confs, ensembledata): else: solv = "lowlevel_gsolv_info" e = "lowlevel_sp_info" - conf.calc_free_energy(e=e, solv=solv, rrho=rrho) + conf.calc_free_energy( + e=e, + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) try: @@ -1716,7 +1766,13 @@ def part2(config, conformers, store_confs, ensembledata): else: solv = "lowlevel_gsolv_info" e = "lowlevel_sp_info" - conf.calc_free_energy(e=e, solv=solv, rrho=rrho, t=temperature) + conf.calc_free_energy( + 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] @@ -1733,9 +1789,11 @@ def part2(config, conformers, store_confs, ensembledata): for conf in calculate: avG += conf.bm_weight * conf.free_energy avE += conf.bm_weight * conf.lowlevel_sp_info["energy"] - avGRRHO += conf.bm_weight * conf.lowlevel_grrho_info["range"].get( - temperature, 0.0 - ) + avGRRHO += conf.bm_weight * conf.get_mrrho( + temperature, + "lowlevel_grrho_info", + config.consider_sym + ) avGsolv += conf.bm_weight * conf.lowlevel_gsolv_info["range"].get( temperature, 0.0 ) @@ -1791,7 +1849,13 @@ def part2(config, conformers, store_confs, ensembledata): else: solv = "lowlevel_gsolv_info" e = "lowlevel_sp_info" - conf.calc_free_energy(e=e, solv=solv, rrho=rrho) + conf.calc_free_energy( + e=e, + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) # ensembledata is used to store avGcorrection ensembledata.comment = [ lowestconf, @@ -1843,6 +1907,8 @@ def part2(config, conformers, store_confs, ensembledata): config.func, config.nat, "free_energy", + config.temperature, + config.consider_sym, overwrite=True, **kwargs, ) @@ -1900,6 +1966,8 @@ def part2(config, conformers, store_confs, ensembledata): config.func, config.nat, "free_energy", + config.temperature, + config.consider_sym, **kwargs, ) @@ -1919,7 +1987,13 @@ def part2(config, conformers, store_confs, ensembledata): ) as best: best.write( "$coord # {} {} !CONF{} \n".format( - conf.free_energy, conf.lowlevel_grrho_info["energy"], conf.id + conf.free_energy, + conf.get_mrrho( + config.temperature, + rrho="lowlevel_grrho_info", + consider_sym=config.consider_sym + ), + conf.id ) ) for line in coord[1:]: diff --git a/censo_qm/prescreening.py b/censo_qm/prescreening.py index 0ef114a..14649dc 100755 --- a/censo_qm/prescreening.py +++ b/censo_qm/prescreening.py @@ -5,6 +5,7 @@ import os import sys import math +import time from multiprocessing import JoinableQueue as Queue from .cfg import PLENGTH, DIGILEN, AU2KCAL, CODING, WARNLEN from .parallel import run_in_parallel @@ -77,20 +78,19 @@ def part1(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": @@ -280,6 +280,7 @@ def part1(config, conformers, store_confs, ensembledata): folder = instruction["func"] check = {True: "was successful", False: "FAILED"} + start_firstsort = time.perf_counter() if calculate: print(f"The {name} is calculated for:") print_block(["CONF" + str(i.id) for i in calculate]) @@ -313,7 +314,7 @@ def part1(config, conformers, store_confs, ensembledata): calculate, instruction, config.balance, - folder + folder, ) for conf in list(calculate): @@ -335,6 +336,7 @@ 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["info"] = "calculated" conf.prescreening_gsolv_info["method"] = conf.job["method2"] elif instruction["jobtype"] in ( @@ -362,6 +364,7 @@ def part1(config, conformers, store_confs, ensembledata): conf.prescreening_sp_info["method"] = conf.job["method"] conf.prescreening_gsolv_info["gas-energy"] = conf.job["energy"] conf.prescreening_gsolv_info["energy"] = conf.job["energy2"] + conf.prescreening_gsolv_info["range"] = conf.job["erange1"] conf.prescreening_gsolv_info["info"] = "calculated" conf.prescreening_gsolv_info["method"] = conf.job["method2"] else: @@ -404,6 +407,8 @@ def part1(config, conformers, store_confs, ensembledata): f"{conf.prescreening_gsolv_info['energy']:>.8f}" ) calculate.append(prev_calculated.pop(prev_calculated.index(conf))) + end_firstsort = time.perf_counter() + ensembledata.part_info["part1_firstsort"] = end_firstsort -start_firstsort for conf in calculate: conf.reset_job_info() if not calculate: @@ -425,7 +430,13 @@ def part1(config, conformers, store_confs, ensembledata): else: solv = "prescreening_gsolv_info" e = "prescreening_sp_info" - conf.calc_free_energy(e=e, solv=solv, rrho=rrho) + conf.calc_free_energy( + e=e, + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) try: minfree = min([i.free_energy for i in calculate if i is not None]) except ValueError: @@ -435,7 +446,9 @@ def part1(config, conformers, store_confs, 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")), @@ -494,14 +507,16 @@ def part1(config, conformers, store_confs, ensembledata): if conf.rel_free_energy > (config.part1_threshold): store_confs.append(calculate.pop(calculate.index(conf))) if calculate: - print(f"Below the g_thr(1) threshold of {config.part1_threshold} kcal/mol.\n") + print( + f"Below the g_thr(1) threshold of {config.part1_threshold} kcal/mol.\n" + ) print_block(["CONF" + str(i.id) for i in calculate]) else: print(f"{'ERROR:':{WARNLEN}}There are no more conformers left!") else: print( "\nAll relative (free) energies are below the g_thr(1) threshold " - f"of ({config.part1_threshold} kcal/mol.\nAll conformers are " + f"of {config.part1_threshold} kcal/mol.\nAll conformers are " "considered further." ) ensembledata.nconfs_per_part["part1_firstsort"] = len(calculate) @@ -564,7 +579,7 @@ def part1(config, conformers, store_confs, ensembledata): "success": False, "imagthr": config.imagthr, "sthr": config.sthr, - "scale":config.scale, + "scale": config.scale, } instruction_prerrho["method"], _ = config.get_method_name( @@ -601,20 +616,33 @@ def part1(config, conformers, store_confs, ensembledata): # check if too many calculations failed ### + for conf in list(calculate): print( - f"The prescreening G_mRRHO calculation @ {conf.job['symmetry']} " + f"The prescreening G_mRRHO run @ {conf.job['symmetry']} " f"{check[conf.job['success']]} for " f"{last_folders(conf.job['workdir'], 2):>{pl}}: " - f"{conf.job['energy']:>.8f}" + f"{conf.job['energy']:>.8f} " + f"""S_rot(sym)= {conf.calc_entropy_sym( + config.temperature, + symnum=conf.job['symnum']):>.7f} """ + f"""using= {conf.get_mrrho( + config.temperature, + rrho='direct_input', + consider_sym=config.consider_sym, + symnum=conf.job['symnum'], + direct_input=conf.job['energy']):>.7f}""" ) if not conf.job["success"]: conf.prescreening_grrho_info["info"] = "failed" store_confs.append(calculate.pop(calculate.index(conf))) else: conf.sym = conf.job["symmetry"] + conf.symnum = conf.job["symnum"] + conf.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["info"] = "calculated" conf.prescreening_grrho_info["method"] = instruction_prerrho[ "method" @@ -639,10 +667,15 @@ def part1(config, conformers, store_confs, ensembledata): os.path.join(config.cwd, "CONF" + str(conf.id), folderrho) ) print( - f"The prescreening G_mRRHO calculation @ {conf.sym} " + f"The prescreening G_mRRHO run @ {conf.sym} " f"{check[conf.job['success']]} for " f"{last_folders(conf.job['workdir'], 2):>{pl}}: " - f"{conf.prescreening_grrho_info['energy']:>.8f}" + f"{conf.prescreening_grrho_info['energy']:>.8f} " + f"S_rot(sym)= {conf.calc_entropy_sym(config.temperature):>.7f}" + f""" using= {conf.get_mrrho( + config.temperature, + rrho='prescreening_grrho_info', + consider_sym=config.consider_sym):>.7f}""" ) calculate.append(prev_calculated.pop(prev_calculated.index(conf))) if not calculate: @@ -670,7 +703,10 @@ def part1(config, conformers, store_confs, ensembledata): lambda conf: getattr(conf, "rel_xtb_free_energy"), lambda conf: getattr(conf, "prescreening_sp_info")["energy"], lambda conf: getattr(conf, "prescreening_gsolv_info")["energy"], - lambda conf: getattr(conf, "prescreening_grrho_info")["energy"], + # lambda conf: getattr(conf, "prescreening_grrho_info")["energy"], + lambda conf: conf.get_mrrho( + config.temperature, "prescreening_grrho_info", config.consider_sym + ), lambda conf: getattr(conf, "free_energy"), lambda conf: getattr(conf, "rel_free_energy"), ] @@ -729,9 +765,20 @@ def part1(config, conformers, store_confs, ensembledata): solv = None else: solv = "prescreening_gsolv_info" - conf.calc_free_energy(e="prescreening_sp_info", solv=solv, rrho=rrho) + conf.calc_free_energy( + e="prescreening_sp_info", + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) conf.xtb_free_energy = conf.calc_free_energy( - e="xtb_energy", solv=None, rrho=rrho, out=True + e="xtb_energy", + solv=None, + rrho=rrho, + out=True, + t=config.temperature, + consider_sym=config.consider_sym ) try: @@ -745,7 +792,9 @@ def part1(config, conformers, store_confs, 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 calculate.sort(key=lambda x: int(x.id)) printout( @@ -782,9 +831,17 @@ def part1(config, conformers, store_confs, ensembledata): else: std_dev = calc_std_dev( [ - conf.prescreening_grrho_info["energy"] * AU2KCAL + conf.get_mrrho( + config.temperature, + rrho="prescreening_grrho_info", + consider_sym=config.consider_sym + ) * AU2KCAL for conf in calculate - if conf.prescreening_grrho_info["energy"] is not None + if conf.get_mrrho( + config.temperature, + rrho="prescreening_grrho_info", + consider_sym=config.consider_sym + ) is not None ] ) max_fuzzy = 1 @@ -810,7 +867,13 @@ def part1(config, conformers, store_confs, ensembledata): else: solv = "prescreening_gsolv_info" e = "prescreening_sp_info" - conf.calc_free_energy(e=e, solv=solv, rrho=rrho) + conf.calc_free_energy( + e=e, + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) try: minfree = min([i.free_energy for i in calculate if i is not None]) except ValueError: @@ -828,7 +891,13 @@ def part1(config, conformers, store_confs, ensembledata): else: solv = "prescreening_gsolv_info" e = "prescreening_sp_info" - conf.calc_free_energy(e=e, solv=solv, rrho=rrho) + conf.calc_free_energy( + e=e, + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) try: minfree = min([i.free_energy for i in calculate if i is not None]) except ValueError: @@ -874,7 +943,9 @@ def part1(config, conformers, store_confs, 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["part1"] = "passed" @@ -904,7 +975,13 @@ def part1(config, conformers, store_confs, ensembledata): solv = None else: solv = "prescreening_gsolv_info" - conf.calc_free_energy(e="prescreening_sp_info", solv=solv, rrho=rrho) + conf.calc_free_energy( + e="prescreening_sp_info", + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) # write coord.enso_best for conf in calculate: @@ -923,7 +1000,11 @@ def part1(config, conformers, store_confs, ensembledata): best.write( "$coord # {} {} !CONF{} \n".format( conf.free_energy, - conf.prescreening_grrho_info["energy"], + conf.get_mrrho( + config.temperature, + rrho="prescreening_grrho_info", + consider_sym=config.consider_sym + ), conf.id, ) ) @@ -975,7 +1056,13 @@ def part1(config, conformers, store_confs, ensembledata): else: solv = "prescreening_gsolv_info" e = "prescreening_sp_info" - conf.calc_free_energy(e=e, solv=solv, rrho=rrho) + conf.calc_free_energy( + e=e, + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) avG = 0.0 @@ -985,7 +1072,12 @@ 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 + ) avGsolv += conf.bm_weight * conf.prescreening_gsolv_info["energy"] # printout: @@ -1011,25 +1103,33 @@ def part1(config, conformers, store_confs, ensembledata): print("".ljust(int(PLENGTH), "-")) print("") - ################################################################################ + ############################################################################ + # Calculate unbiased GFNn-xTB energy: + for conf in list(calculate): + if conf.xtb_energy_unbiased is None: + pass + else: + conf = calculate.pop(calculate.index(conf)) + conf.job["success"] = True + prev_calculated.append(conf) - print("\nCalculating unbiased GFNn-xTB energy") - instruction_gfn = { - "jobtype": "xtb_sp", - "func": getattr(config, "part1_gfnv"), - "charge": config.charge, - "unpaired": config.unpaired, - "solvent": config.solvent, - "progpath": config.external_paths["xtbpath"], - "sm": config.sm_rrho, - "rmsdbias": config.rmsdbias, - "temperature": config.temperature, - "gfn_version": config.part1_gfnv, - "energy": 0.0, - "energy2": 0.0, - "success": False, - } if calculate: + print("\nCalculating unbiased GFNn-xTB energy") + instruction_gfn = { + "jobtype": "xtb_sp", + "func": getattr(config, "part1_gfnv"), + "charge": config.charge, + "unpaired": config.unpaired, + "solvent": config.solvent, + "progpath": config.external_paths["xtbpath"], + "sm": config.sm_rrho, + "rmsdbias": config.rmsdbias, + "temperature": config.temperature, + "gfn_version": config.part1_gfnv, + "energy": 0.0, + "energy2": 0.0, + "success": False, + } folder_gfn = "GFN_unbiased" save_errors, store_confs, calculate = new_folders( config.cwd, calculate, folder_gfn, save_errors, store_confs @@ -1055,6 +1155,11 @@ def part1(config, conformers, store_confs, ensembledata): conf.xtb_energy_unbiased = conf.xtb_energy else: conf.xtb_energy_unbiased = conf.job["energy"] + if prev_calculated: + for conf in list(prev_calculated): + calculate.append(prev_calculated.pop(prev_calculated.index(conf))) + ############################################################################ + # write ensemble move_recursively(config.cwd, "enso_ensemble_part1.xyz") if config.evaluate_rrho: @@ -1068,6 +1173,8 @@ def part1(config, conformers, store_confs, ensembledata): config.func, config.nat, "free_energy", + config.temperature, + config.consider_sym, **kwargs, ) @@ -1078,17 +1185,11 @@ def part1(config, conformers, store_confs, ensembledata): conf.bm_weight = 0.0 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 Part1")) / 2) print("\n" + "".ljust(tmp, ">") + "END of Part1" + "".rjust(tmp, "<")) diff --git a/censo_qm/qm_job.py b/censo_qm/qm_job.py index db97aa1..994e6c8 100644 --- a/censo_qm/qm_job.py +++ b/censo_qm/qm_job.py @@ -11,7 +11,7 @@ import time import subprocess import json -from .cfg import ENVIRON, CODING, WARNLEN, censo_solvent_db +from .cfg import ENVIRON, CODING, WARNLEN, censo_solvent_db, rot_sym_num from .utilities import last_folders, print from .datastructure import MoleculeData @@ -87,6 +87,7 @@ def reset_job_info(self): "internal_error": [], # "cosmorsparam": "", # normal/fine + "symnum":1, } def _sp(self, silent=False): @@ -109,6 +110,25 @@ def _genericoutput(self): """ pass + 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 _xtb_sp(self): """ Get single-point energy from xTB @@ -141,9 +161,16 @@ def _xtb_sp(self): [ "--" + str(self.job["sm"]), censo_solvent_db[self.job["solvent"]]["xtb"][1], + "-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: @@ -293,6 +320,12 @@ def _xtb_gsolv(self): # run single-point in solution: # ``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: @@ -308,6 +341,8 @@ def _xtb_gsolv(self): "--chrg", str(self.job["charge"]), "--norestart", + "-I", + "xcontrol-inp", ], shell=False, stdin=None, @@ -424,8 +459,12 @@ def _xtbrrho(self): if self.job["consider_sym"]: # xcout.write(" desy=0.1\n") # taken from xtb defaults xcout.write(" maxat=1000\n") - else: - xcout.write(" desy=0.0\n") + # always consider symmetry! + #else: + # xcout.write(" desy=0.0\n") + if self.job["solvent"] != "gas": + xcout.write("$gbsa\n") + xcout.write(" gbsagrid=tight\n") xcout.write("$end\n") if self.job["bhess"]: # set ohess or bhess @@ -497,47 +536,49 @@ def _xtbrrho(self): print(self.job["errormessage"][-1]) return # start reading output! - if self.job["trange"]: - if not os.path.isfile(os.path.join(self.job["workdir"], "ohess.out")): - self.job["energy"] = 0.0 - self.job["success"] = False - self.job["errormessage"].append( - f"{'ERROR:':{WARNLEN}}file {self.job['workdir']}/ohess.out could not be found!" - ) - print(self.job["errormessage"][-1]) - return - gt = {} - ht = {} - rotS = {} - with open( + if not os.path.isfile(os.path.join(self.job["workdir"], "ohess.out")): + self.job["energy"] = 0.0 + self.job["success"] = False + self.job["errormessage"].append( + f"{'ERROR:':{WARNLEN}}file {self.job['workdir']}/ohess.out could not be found!" + ) + 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() - # rotational entropy: - for line in store: - if "VIB" in line: - try: - realline = store.index(line) + 1 - T = float(line.split()[0]) - rotS[T] = float(store[realline].split()[4]) - except (KeyError, ValueError): - pass - for line in store: - if "T/K" in line: - start = store.index(line) - for line in store[start + 2 :]: - if "----------------------------------" in line: - break - else: - try: - T = float(line.split()[0]) - gt[T] = float(line.split()[4]) - ht[T] = float(line.split()[2]) - except (ValueError, KeyError): - print(f"{'ERROR:':{WARNLEN}}can not convert G(T)") + + if self.job["trange"]: + gt = {} + ht = {} + rotS = {} + # rotational entropy: + for line in store: + if "VIB" in line: + try: + realline = store.index(line) + 1 + T = float(line.split()[0]) + rotS[T] = float(store[realline].split()[4]) + except (KeyError, ValueError): + pass + for line in store: + if "T/K" in line: + start = store.index(line) + for line in store[start + 2 :]: + if "----------------------------------" in line: + break + else: + try: + T = float(line.split()[0]) + gt[T] = float(line.split()[4]) + ht[T] = float(line.split()[2]) + except (ValueError, KeyError): + print(f"{'ERROR:':{WARNLEN}}can not convert G(T)") if len(self.job["trange"]) == len(gt): self.job["success"] = True self.job["erange1"] = gt @@ -547,21 +588,22 @@ def _xtbrrho(self): self.job["success"] = False return # end self.trange - if self.job["bhess"]: - # read rmsd_info - with open( - os.path.join(self.job["workdir"], "ohess.out"), - "r", - encoding=CODING, - newline=None, - ) as inp: - store = inp.readlines() - for line in store: - if "final rmsd / " in line: - try: - self.job["rmsd"] = float(line.split()[3]) - except (ValueError): - self.job["rmsd"] = 0.0 + for line in store: + if "final rmsd / " in line and self.job['bhess']: + try: + self.job["rmsd"] = float(line.split()[3]) + except (ValueError): + self.job["rmsd"] = 0.0 + if ": linear? " in line: + # linear needed for symmetry and S_rot (only if 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 + print(e) if os.path.isfile(os.path.join(self.job["workdir"], "xtb_enso.json")): with open( os.path.join(self.job["workdir"], "xtb_enso.json"), @@ -590,12 +632,16 @@ def _xtbrrho(self): self.job["success"] = True if "point group" in data: self.job["symmetry"] = data["point group"] + if "linear" in data: + self.job["linear"] = data.get("linear", self.job["linear"]) else: print( f"{'ERROR:':{WARNLEN}}while converting mRRHO in: {last_folders(self.job['workdir'], 2)}" ) 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"]) else: print( f"{'WARNING:':{WARNLEN}}File " @@ -604,6 +650,7 @@ 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 b8513fa..d1de651 100755 --- a/censo_qm/refinement.py +++ b/censo_qm/refinement.py @@ -662,7 +662,16 @@ def part3(config, conformers, store_confs, ensembledata): f"The G_mRRHO calculation @ {conf.job['symmetry']} " f"{check[conf.job['success']]} for " f"{last_folders(conf.job['workdir'], 2):>{pl}}: " - f"{conf.job['energy']:>.8f}" + f"{conf.job['energy']:>.8f} " + f"""S_rot(sym)= {conf.calc_entropy_sym( + config.temperature, + symnum=conf.job['symnum']):>.7f} """ + f"""using= {conf.get_mrrho( + config.temperature, + rrho='direct_input', + consider_sym=config.consider_sym, + symnum=conf.job['symnum'], + direct_input=conf.job["energy"]):>.7f}""" ) if not conf.job["success"]: conf.part_info["part3"] = "refused" @@ -671,6 +680,7 @@ def part3(config, conformers, store_confs, ensembledata): store_confs.append(calculate.pop(calculate.index(conf))) else: conf.sym = conf.job["symmetry"] + conf.linear = conf.job["linear"] conf.highlevel_grrho_info["rmsd"] = conf.job["rmsd"] conf.highlevel_grrho_info["energy"] = conf.job["energy"] conf.highlevel_grrho_info["info"] = "calculated" @@ -701,7 +711,12 @@ def part3(config, conformers, store_confs, ensembledata): f"The G_mRRHO calculation @ {conf.sym} " f"{check[conf.job['success']]} for " f"{last_folders(conf.job['workdir'], 2):>{pl}}: " - f"{conf.highlevel_grrho_info['energy']:>.8f}" + f"{conf.highlevel_grrho_info['energy']:>.8f} " + f"S_rot(sym)= {conf.calc_entropy_sym(config.temperature):>.7f}" + f""" using= {conf.get_mrrho( + config.temperature, + rrho='highlevel_grrho_info', + consider_sym=config.consider_sym):>.7f}""" ) calculate.append(prev_calculated.pop(prev_calculated.index(conf))) if not calculate: @@ -727,9 +742,13 @@ 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 + ), lambda conf: getattr(conf, "free_energy"), lambda conf: getattr(conf, "rel_free_energy"), + lambda conf: getattr(conf, "bm_weight") * 100, ] columnheader = [ "CONF#", @@ -740,17 +759,38 @@ def part3(config, conformers, store_confs, ensembledata): "GmRRHO [Eh]", "Gtot", "ΔGtot", + "Boltzmannweight", ] - columndescription = ["", "[a.u.]", "[kcal/mol]", "", "", "", "[Eh]", "[kcal/mol]"] - columndescription2 = ["", "", "", "", "", "", "", ""] - columnformat = ["", (12, 7), (5, 2), (12, 7), (12, 7), (12, 7), (12, 7), (5, 3)] + columndescription = [ + "", + "[a.u.]", + "[kcal/mol]", + "", + "", + "", + "[Eh]", + "[kcal/mol]", + f" % at {config.temperature:.2f} K", + ] + columndescription2 = ["", "", "", "", "", "", "", "", ""] + columnformat = [ + "", + (12, 7), + (5, 2), + (12, 7), + (12, 7), + (12, 7), + (12, 7), + (5, 3), + (5, 2) + ] if config.solvent == "gas": columndescription[3] = instruction["method"] elif config.solvent != "gas": columndescription[3] = instruction["method"] columndescription[4] = instruction["method2"] if config.evaluate_rrho: - columndescription[5] = str(instruction_rrho["method"]).upper() # Grrho + columndescription[5] = instruction_rrho["method"] if not config.evaluate_rrho or config.solvent == "gas": if not config.evaluate_rrho: # ignore rrho in printout @@ -776,7 +816,13 @@ def part3(config, conformers, store_confs, ensembledata): else: solv = "highlevel_gsolv_info" e = "highlevel_sp_info" - conf.calc_free_energy(e=e, solv=solv, rrho=rrho) + conf.calc_free_energy( + e=e, + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) try: @@ -858,13 +904,19 @@ def part3(config, conformers, store_confs, ensembledata): else: solv = "highlevel_gsolv_info" e = "highlevel_sp_info" - conf.calc_free_energy(e=e, solv=solv, rrho=rrho, t=temperature) - try: - minfreeT = min( - [conf.free_energy for conf in calculate if conf.free_energy is not None] - ) - except ValueError: - raise ValueError + conf.calc_free_energy( + 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 calculate = calc_boltzmannweights(calculate, "free_energy", temperature) avG = 0.0 @@ -875,8 +927,13 @@ 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( - temperature, 0.0 + #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 ) avGsolv += conf.bm_weight * conf.highlevel_gsolv_info["range"].get( temperature, 0.0 @@ -936,7 +993,13 @@ def part3(config, conformers, store_confs, ensembledata): else: solv = "highlevel_gsolv_info" e = "highlevel_sp_info" - conf.calc_free_energy(e=e, solv=solv, rrho=rrho) + conf.calc_free_energy( + e=e, + solv=solv, + rrho=rrho, + t=config.temperature, + consider_sym=config.consider_sym + ) calculate = calc_boltzmannweights(calculate, "free_energy", config.temperature) # SORTING for the next part: @@ -991,6 +1054,8 @@ def part3(config, conformers, store_confs, ensembledata): folder, config.nat, "free_energy", + config.temperature, + config.consider_sym, **kwargs, ) @@ -1010,7 +1075,13 @@ def part3(config, conformers, store_confs, ensembledata): ) as best: best.write( "$coord # {} {} !CONF{} \n".format( - conf.free_energy, conf.highlevel_grrho_info["energy"], conf.id + conf.free_energy, + conf.get_mrrho( + config.temperature, + rrho="highlevel_grrho_info", + consider_sym=config.consider_sym + ), + conf.id ) ) for line in coord[1:]: diff --git a/censo_qm/setupcenso.py b/censo_qm/setupcenso.py index 24699cf..10b518c 100755 --- a/censo_qm/setupcenso.py +++ b/censo_qm/setupcenso.py @@ -55,8 +55,29 @@ def enso_startup(cwd, args): print("Removed files and going to exit!") sys.exit(0) + configfname = ".censorc" if args.writeconfig: + #check if .censorc in local or home dir and ask user it the program + # paths should be copied + tmp = None newconfigfname = "censorc_new" + if os.path.isfile(os.path.join(config.cwd, configfname)): + tmp = os.path.join(config.cwd, configfname) + elif os.path.isfile(os.path.join(os.path.expanduser("~"), configfname)): + 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("Please type 'yes' or 'no':") + user_input = input() + if user_input.strip() in ('y', 'yes'): + config.read_program_paths(tmp) + config.write_rcfile(os.path.join(config.cwd, newconfigfname), usepaths=True) + print("") + 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)) print( "A new ensorc was written into the current directory file: " f"{newconfigfname}!\nYou have to adjust the settings to your needs" @@ -65,10 +86,8 @@ def enso_startup(cwd, args): "and place it either in your /home/$USER/ or current directory.\n" "\nAll done!" ) - config.write_rcfile(os.path.join(config.cwd, newconfigfname)) sys.exit(0) - configfname = ".censorc" if args.inprcpath: try: tmp_path = os.path.abspath(os.path.expanduser(args.inprcpath)) @@ -179,6 +198,14 @@ def enso_startup(cwd, args): 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("\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) @@ -376,6 +403,7 @@ 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), avGcorrection=save_data[conf].get("avGcorrection"), comment=save_data[conf].get("comment"), bestconf=save_data[conf].get("bestconf"), @@ -398,8 +426,10 @@ def enso_startup(cwd, args): xtb_free_energy=save_data[conf].get("xtb_free_energy"), 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"), - gi=save_data[conf].get("gi"), + 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")), + gi=save_data[conf].get("gi", getattr(MoleculeData(0), "sym")), removed=save_data[conf].get( "removed", getattr(MoleculeData(0), "removed") ), @@ -638,6 +668,22 @@ def enso_startup(cwd, args): func2=config.func_or_scf, ) molecule.load_prev("optical_rotation_info", method) + ##### + 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" + ): + # reset to default, load is not available + molecule.nmr_shielding_info=getattr(MoleculeData(0), "nmr_shielding_info") + ##### elif value and key in ("func3", "basis3"): # save calculated to molecule.save_prev( diff --git a/censo_qm/tm_job.py b/censo_qm/tm_job.py index 7eaee43..c54b90c 100644 --- a/censo_qm/tm_job.py +++ b/censo_qm/tm_job.py @@ -746,7 +746,6 @@ def _cosmors(self): ) as inp: stor = inp.readlines() for line in stor: - #vwork = 0 if "T=" in line: T = float(line.split()[5]) vwork = R * T * math.log(videal * T) diff --git a/censo_qm/utilities.py b/censo_qm/utilities.py index 2d303bb..d7b20f9 100755 --- a/censo_qm/utilities.py +++ b/censo_qm/utilities.py @@ -172,7 +172,7 @@ def x2t(path, infile="inp.xyz"): def write_trj( - results, cwd, outpath, optfolder, nat, attribute, 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,8 +199,13 @@ def write_trj( ### coordinates in xyz out.write(" {}\n".format(nat)) xtbfree = conf.calc_free_energy( - e=energy, solv=None, rrho=rrho, out=True - ) + e=energy, + solv=None, + rrho=rrho, + out=True, + t=temperature, + consider_sym=consider_sym + ) if xtbfree is not None: xtbfree = f"{xtbfree:20.8f}" out.write( @@ -779,8 +784,13 @@ def crest_routine(config, conformers, func, store_confs, prev_calculated=None): def format_line(key, value, options, optionlength=70, dist_to_options=30): """ used in print_parameters + 'key: value ---space--- #(separator) option1 option2 ...' """ # limit printout of possibilities + if len(options) > 0: + separator = '#' + else: + separator = '' if len(str(options)) > optionlength: length = 0 reduced = [] @@ -790,8 +800,8 @@ def format_line(key, value, options, optionlength=70, dist_to_options=30): reduced.append(item) reduced.append("...") options = reduced - line = "{}: {:{digits}} # {} \n".format( - key, str(value), options, digits=dist_to_options - len(key) + line = "{}: {:{digits}} {} {} \n".format( + key, str(value), separator, options, digits=dist_to_options - len(key) ) return line