README update (#6)

* more functionals for func3 available * SD on OR and NMR available for a fixed value of 0.4 kcal/mol
fabothch · Mar 4, 2021 · c5d9354 · c5d9354
1 parent 80bb98f
commit c5d9354
Show file tree

Hide file tree

Showing 8 changed files with 154 additions and 66 deletions.
diff --git a/README.rst b/README.rst
@@ -1,22 +1,35 @@
 [![Latest Version](https://img.shields.io/github/v/release/grimme-lab/CENSO)](https://github.com/grimme-lab/CENSO/releases/latest)
 
+====================================================================
 CENSO - Commandline ENergetic SOrting of Conformer Rotamer Ensembles
 ====================================================================
 
-This repository hosts the `CENSO` code for the refinement of Conformer Rotamer 
-Ensembles (CRE) as obtained from `CREST`.
+.. figure:: docs/src/censo_logo_300dpi.png
+	:scale: 100%
+	:align: center
+	:alt: CENSO logo
 
+This repository hosts the ``CENSO`` code for the refinement of Conformer Rotamer 
+Ensembles (CRE) as obtained from ``CREST``.
 
 
 Installation
-------------
+============
 
 There are several options possible. The easiest is to use the packaged censo programs
 (by use of Pyinstaller) which can be found at the release section. The packaged 
 censo is linked against GLIBC version 2.19 and will work for GLIBC version 2.19 and above.
 
-Other options to use censo are shown below:
+For the packaged "binary" of CENSO, download it from the 
+`release <https://github.com/grimme-lab/CENSO/releases/>`_ site, 
+copy to your bin and make executable:
 
+.. code::
+
+    $ cp censo ~/bin/censo
+    $ chmod u+x ~/bin/censo
+
+Other options to use censo are shown below:
 
 Download the git repository and run:
 
@@ -45,7 +58,7 @@ Flexible Invocation
 
 
 Getting started:
-----------------
+================
 
 Create the remote configuration file .censorc where the user can adjust default
 settings and provide paths to the external programs e.g. `xtb`, `crest`, `orca` ...
@@ -58,6 +71,11 @@ settings and provide paths to the external programs e.g. `xtb`, `crest`, `orca`
     vi /home/$USER/.ensorc
 
 
+Documentation:
+==============
+
+Can be found following: https://fbohle.gitbook.io/censo/
+
 **Interactive Documentation can be accessed:**
 
 .. code::
@@ -70,11 +88,39 @@ Explainations on the commandline arguments can be printed by:
 
     $ censo --help
 
-Online Documentation:
----------------------
 
-Can be found following: https://fbohle.gitbook.io/censo/
+Requirements:
+-------------
 
+* newest xtb (currently: https://github.com/grimme-lab/xtb/releases/tag/v6.4.0 )
+* newest cefine https://github.com/grimme-lab/cefine/releases
+* ORCA > version 4.1
+
+
+Further information (will be ordered later on):
+
+* the file .censorc can be used in the current working directory and will be preferred to 
+  the global configuration file in ~/.censorc
+* a folder ~/.censo_assets/ will be created upon usage of censo
+* ORCA has not been tested extensively so please be careful, test calculations
+  and report possible "bad" settings
+* To be efficient COSMO-RS calculations are not performed with BP86 but whith the functionals
+  for energy evaluation.
+
+
+Usage:
+------
+
+.. code::
+
+    # check if setting-combinations match:
+    $ censo -inp structure_ensemble.xyz -part2 on -solvent h2o --checkinput
+    # start the calculation:
+    $ censo -inp structure_ensemble.xyz -part2 on -solvent h2o > censo.out 2> error.censo &
+
+
+Short overview:
+---------------
 
 The molecule numbering from the input structure ensemble is kept throughout the 
 entire program. There are several program parts which can be used to filter a structure 
@@ -100,42 +146,17 @@ ensemble:
     for the populated ensemble.
 
 
-Usage:
-------
-
-.. code::
-
-    # check if setting-combinations match:
-    $ censo -inp structure_ensemble.xyz -part2 on -solvent h2o --checkinput
-    # start the calculation:
-    $ censo -inp structure_ensemble.xyz -part2 on -solvent h2o > censo.out 2> error.censo &
-
-Requirements:
--------------
-
-* newest xtb (currently: https://github.com/grimme-lab/xtb/releases/tag/v6.4.0 )
-* newest cefine https://github.com/grimme-lab/cefine/releases
-* ORCA > version 4.1
-
-
-Further information (will be ordered later on):
-
-* the file .censorc can be used in the current working directory and will be preferred to 
-  the global configuration file in ~/.censorc
-* a folder ~/.censo_assets/ will be created upon usage of censo
-* ORCA has not been used extensively so please be careful, test calculations
-  and report possible "bad" settings
-* To be efficient COSMO-RS calculations are not performed with BP86 but whith the functionals
-  for energy evaluation.
-
-
-
+For Turbomole user:
+-------------------
 
-License
--------
+The amount of *ricore* for each calculation can be set in your `.cefinerc`. The same
+holds for *maxcor* and/or *rpacor*.
 
-LGPL3
+.. code::
 
+    $ echo "ricore  4000" > .cefinerc
+    $ echo "maxcor  4000" >> .cefinerc
+    $ echo "rpacor  4000" >> .cefinerc
 
 Available solvation models:
 ---------------------------
@@ -157,20 +178,6 @@ to free energy (Gibbs energy) of the choosen geometry (SMGSOLV).
     "","ALPB_Gsolv","(SMGSOLV)"
 
 
-
-For Turbomole user:
--------------------
-
-The amount of *ricore* for each calculation can be set in your `.cefinerc`. The same
-holds for *maxcor* and/or *rpacor*.
-
-.. code::
-
-    $ echo "ricore  4000" > .cefinerc
-    $ echo "maxcor  4000" >> .cefinerc
-    $ echo "rpacor  4000" >> .cefinerc
-
-
 Solvents:
 ---------
 
@@ -192,3 +199,16 @@ a replacement solvent. E.g. if CCl4 is not available choose CHCl3.
 
 The solvent file is directly used in `CENSO` and typos will cause calculations to crash!
 Adding a new solvent is as easy as adding a new dictionary to the file.
+
+License
+-------
+
+``CENSO`` is free software: you can redistribute it and/or modify it under
+the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+``CENSO`` is distributed in the hope that it will be useful,
+but without any warranty; without even the implied warranty of
+merchantability or fitness for a particular purpose. See the
+GNU Lesser General Public License for more details.
diff --git a/censo_qm/cfg.py b/censo_qm/cfg.py
@@ -5,7 +5,7 @@
 """
 import os
 
-__version__ = "1.0.5"
+__version__ = "1.0.6"
 
 DESCR = f"""
          ______________________________________________________________

diff --git a/censo_qm/inputhandling.py b/censo_qm/inputhandling.py
@@ -1145,8 +1145,8 @@ 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"]
-    func3_tm = ["pw6b95", "pbe0", "b97-d", "r2scan-3c", "b97-3c", "wb97x-v"]
+    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"]

diff --git a/censo_qm/nmrproperties.py b/censo_qm/nmrproperties.py
@@ -83,13 +83,15 @@ def average_shieldings(config, calculate, element_ref_shield, energy, solv, rrho
     if not os.path.isfile(path_anmr_nucinfo):
         print(f"File anmr_nucinfo (generated by CREST) is not present!")
         return
+    print("\nAveraged shielding constants:")
     chemeq = read_chemeq(path_anmr_nucinfo)
     averaged = {}
     element = {}
-
     sigma_std_dev = {}
+    sigma_std_dev_const = {}
     for i in range(1, config.nat + 1):
         sigma_std_dev[i] = []
+        sigma_std_dev_const[i] = []
 
     for conf in calculate:
         # get averaged shielding constants
@@ -127,11 +129,31 @@ def average_shieldings(config, calculate, element_ref_shield, energy, solv, rrho
                     tmp_sigma[atom] += conf.bm_weight * sigma
             for atom in conf.shieldings.keys():
                 sigma_std_dev[atom].append(tmp_sigma[atom])
+    # SD with 0.25 kcal/mol
+    wanted = True
+    if wanted:
+        for _ in range(1000):
+            for conf in calculate:
+                conf.calc_free_energy(e=energy, solv=solv, rrho=rrho)
+                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
+            for conf in calculate:
+                for atom in conf.shieldings.keys():
+                    sigma = sum(
+                        [conf.shieldings.get(eq_atom, 0.0) for eq_atom in chemeq[atom]]
+                    ) / len(chemeq[atom])
+                    tmp_sigma[atom] += conf.bm_weight * sigma
+            for atom in conf.shieldings.keys():
+                sigma_std_dev_const[atom].append(tmp_sigma[atom])
 
-
-    print("\nAveraged shielding constants:")
-    print("# in coord  element  σ(sigma)  SD(σ based on SD Gsolv)   shift    σ_ref")
-    print("".ljust(int(80), "-"))
+
+    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}"
@@ -144,10 +166,15 @@ def average_shieldings(config, calculate, element_ref_shield, energy, solv, rrho
         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} {make_shift(atom):>5}    {element_ref_shield.get(element[atom], 0.0)}"
+                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}")

diff --git a/censo_qm/opticalrotation.py b/censo_qm/opticalrotation.py
@@ -594,6 +594,47 @@ def part5(config, conformers, store_confs, ensembledata):
                 f": {calc_std_dev(all_or):> {max_fmt}.3f}   in deg*[dm(g/cc)]^(-1)"
             )
 
+    if calculate:
+        for freq in config.freq_or:
+            all_or = []
+            for _ in range(1000):
+                averaged_or = 0.0
+                for conf in calculate:
+                    conf.calc_free_energy(e=energy, solv=gsolv, rrho=rrho)
+                    conf.free_energy += normalvariate(
+                        0.0, (0.4/AU2KCAL)
+                    )
+                calculate = calc_boltzmannweights(
+                    calculate, "free_energy", config.temperature
+                )
+                for conf in calculate:
+                    averaged_or += conf.bm_weight * conf.optical_rotation_info[
+                        "range"
+                    ].get(freq)
+                all_or.append(averaged_or)
+            try:
+                max_fmt = max(
+                    [
+                        len(
+                            str(
+                                item.optical_rotation_info["range"].get(
+                                    config.freq_or[0]
+                                )
+                            ).split(".")[0]
+                        )
+                        for item in calculate
+                    ]
+                )
+                max_fmt += 9
+            except Exception as e:
+                print(e)
+                max_fmt = 16
+            print(
+                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)"
+            )
+
+
     for conf in calculate:
         conf.reset_job_info()
 

diff --git a/censo_qm/setupcenso.py b/censo_qm/setupcenso.py
@@ -166,7 +166,7 @@ def enso_startup(cwd, args):
                 if int(censorc_version.split(".")[0]) < int(__version__.split(".")[0]):
                     print(
                         f"{'ERROR:':{WARNLEN}}There has been an API break and you have to "
-                        "create a new .censorc.\n{'':{WARNLEN}}E.g. 'censo -newconfig'"
+                        f"create a new .censorc.\n{'':{WARNLEN}}E.g. 'censo -newconfig'"
                     )
                     sys.exit(1)
                 myfile.seek(0)  # reset reader

diff --git a/docs/src/censo_logo_300dpi.png b/docs/src/censo_logo_300dpi.png
diff --git a/setup.cfg b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = censo-QM
-version = 1.0.5
+version = 1.0.6
 description = CENSO - Comandline ENergetic SOrting for conformer rotamer ensembles
 long_description = file: README.rst
 long_description_content_type = text/x-rst