Added Natural Transition Orbital calculation and processing.

README.md

@@ -81,7 +81,6 @@ This should open your default text editor.
 Copy-paste the following profile configuration (substitute path to ORCA and the Docker image name if needed)
 ```
 port = 8888
-default_apps = [ "aiidalab-ispg@git+https://github.com/ispg-group/aiidalab-ispg.git@main",]
 system_user = "jovyan"
 image = "ghcr.io/ispg-group/atmospec:latest"
 home_mount = "aiidalab_atmospec_home"

aiidalab_ispg/app/atmospec_steps.py

@@ -21,6 +21,7 @@
     MoleculeSettings,
     ResourceSelectionWidget,
     WignerSamplingSettings,
+    OrbitalSettings,
 )
 from .optimization_steps import OptimizationParameters
 from .steps import SubmitWorkChainStepBase, ViewWorkChainStatusStep
@@ -39,6 +40,8 @@ class AtmospecParameters(OptimizationParameters):
     tddft_functional: str
     nwigner: int
     wigner_low_freq_thr: float
+
+    calculate_orbitals: bool
 
 
 DEFAULT_ATMOSPEC_PARAMETERS = AtmospecParameters(
@@ -54,6 +57,8 @@ class AtmospecParameters(OptimizationParameters):
     tddft_functional="wB97X-D4",
     nwigner=0,
     wigner_low_freq_thr=100.0,
+
+    calculate_orbitals=True,
 )
 
 
@@ -81,11 +86,15 @@ def __init__(self, **kwargs):
         )
 
         self.wigner_settings = WignerSamplingSettings()
+
+        self.orbital_settings = OrbitalSettings()
 
         self.codes_selector = CodeSettings()
         self.resources_settings = ResourceSelectionWidget()
 
         self.codes_selector.orca.observe(self._update_state, "value")
+
+        self.codes_selector.orca_plot.observe(self._update_state, "value")
 
         # Set defaults
         self._update_ui_from_parameters(DEFAULT_ATMOSPEC_PARAMETERS)
@@ -96,6 +105,8 @@ def __init__(self, **kwargs):
             self.wigner_settings,
             self.molecule_settings,
             self.excited_state_settings,
+
+            self.orbital_settings,
         ]
         grid_layout = ipw.Layout(
             width="100%",
@@ -122,6 +133,10 @@ def _validate_input_parameters(self) -> bool:
         if self.codes_selector.orca.value is None:
             return False
         return True
+
+        if self.codes_selector.orca_plot.value is None:
+            return False
+        return True
 
     def _wigner_allowed(self):
         # Do not allow Wigner sampling for EOM-CCSD
@@ -188,6 +203,8 @@ def _update_ui_from_parameters(self, parameters: AtmospecParameters) -> None:
         self.excited_state_settings.basis.value = parameters.es_basis
         self.wigner_settings.nwigner.value = parameters.nwigner
         self.wigner_settings.wigner_low_freq_thr.value = parameters.wigner_low_freq_thr
+
+        self.orbital_settings.calculate_orbitals.value = parameters.calculate_orbitals
 
         # Infer the value of the gs_sync checkbox
         if (
@@ -213,6 +230,8 @@ def _get_parameters_from_ui(self) -> AtmospecParameters:
             nstates=self.excited_state_settings.nstates.value,
             nwigner=self.wigner_settings.nwigner.value,
             wigner_low_freq_thr=self.wigner_settings.wigner_low_freq_thr.value,
+
+            calculate_orbitals = self.orbital_settings.calculate_orbitals.value,
         )
 
     @traitlets.observe("process")
@@ -269,7 +288,7 @@ def build_base_orca_params(self, params: AtmospecParameters) -> dict:
             "input_keywords": input_keywords,
         }
 
-    def _add_mdci_orca_params(self, orca_parameters, basis, mdci_method, nroots):
+    def _add_mdci_orca_params(self, orca_parameters, basis, mdci_method, nroots, donto):
         mdci_params = deepcopy(orca_parameters)
         mdci_params["input_keywords"].append(mdci_method.value)
         mdci_params["input_keywords"].append(basis)
@@ -280,6 +299,9 @@ def _add_mdci_orca_params(self, orca_parameters, basis, mdci_method, nroots):
         mdci_params["input_blocks"]["mdci"] = {
             "nroots": nroots,
             "maxcore": MEMORY_PER_CPU,
+
+#             Calculate nto flag
+            "donto": donto,
         }
         # TODO: For efficiency reasons, in might not be necessary to calculated left-vectors
         # to obtain TDM, but we need to benchmark that first.
@@ -289,15 +311,20 @@ def _add_mdci_orca_params(self, orca_parameters, basis, mdci_method, nroots):
         return mdci_params
 
     def _add_tddft_orca_params(
-        self, base_orca_parameters, basis, es_method, functional, nroots
+        self, base_orca_parameters, basis, es_method, functional, nroots, donto
     ):
         tddft_params = deepcopy(base_orca_parameters)
         tddft_params["input_keywords"].append(functional)
         tddft_params["input_keywords"].append(basis)
         tddft_params["input_blocks"]["tddft"] = {
             "nroots": nroots,
             "maxcore": MEMORY_PER_CPU,
+
+#             Calculate nto flag
+            "donto": donto,
         }
+
+
         if es_method == ExcitedStateMethod.TDDFT:
             tddft_params["input_blocks"]["tddft"]["tda"] = "false"
         return tddft_params
@@ -323,7 +350,11 @@ def submit(self, _=None):
         builder = AtmospecWorkChain.get_builder()
 
         builder.code = load_code(self.codes_selector.orca.value)
+
+        builder.plot_code = load_code(self.codes_selector.orca_plot.value)
+
         builder.structure = self.input_structure
+
         base_orca_parameters = self.build_base_orca_params(bp)
         gs_opt_parameters = self._add_optimization_orca_params(
             base_orca_parameters, basis=bp.basis, gs_method=bp.method
@@ -338,6 +369,8 @@ def submit(self, _=None):
                 basis=bp.es_basis,
                 functional=bp.tddft_functional,
                 nroots=bp.nstates,
+
+                donto=bp.calculate_orbitals,
             )
         elif bp.excited_method in (
             ExcitedStateMethod.ADC2,
@@ -348,14 +381,16 @@ def submit(self, _=None):
                 basis=bp.es_basis,
                 mdci_method=bp.excited_method,
                 nroots=bp.nstates,
+
+                donto=bp.calculate_orbitals,
             )
         else:
             msg = f"Excited method {bp.excited_method} not implemented"
             raise NotImplementedError(msg)
 
         builder.optimize = bp.optimize
         builder.opt.orca.parameters = gs_opt_parameters
-        builder.exc.orca.parameters = es_parameters
+        builder.exc.orca.parameters = es_parameters        
 
         num_proc = self.resources_settings.num_mpi_tasks.value
         if num_proc > 1:
@@ -378,6 +413,9 @@ def submit(self, _=None):
         # Fetch GBW file from optimization step, to be used as a guess
         # for subsequent excited state calculations.
         builder.opt.orca.metadata.options.additional_retrieve_list = ["aiida.gbw"]
+
+        # Retrieve .nto and .cube files
+        builder.exc.orca.metadata.options.additional_retrieve_list = ["*.nto", "*.cube"]
 
         # Clean the remote directory by default,
         # we're copying back the main output file and gbw file anyway.

aiidalab_ispg/app/input_widgets.py

@@ -293,6 +293,33 @@ def reset(self):
         self.wigner_low_freq_thr.value = self._LOW_FREQ_THR_DEFAULT
 
 
+class OrbitalSettings(ipw.VBox):
+    title = ipw.HTML(
+        """<div style="padding-top: 0px; padding-bottom: 0px">
+        <h4>Orbital Calculation</h4>
+        </div>"""
+    )
+
+    _CALCULATE_ORBITALS_DEFAULT = True
+
+    def __init__(self, **kwargs):
+        style = {"description_width": "initial"}
+        layout = ipw.Layout(max_width="250px")
+
+        self.calculate_orbitals = ipw.Checkbox(
+            value=self._CALCULATE_ORBITALS_DEFAULT,
+            description="Calculate Natural Transition Orbitals",
+            indent=False,
+            layout=layout,
+        )
+
+        super().__init__(
+            children=[self.title, self.calculate_orbitals]
+        )
+
+    def reset(self):
+        self.charge.value = self._CALCULATE_ORBITALS_DEFAULT
+
 class CodeSettings(ipw.VBox):
     codes_title = ipw.HTML(
         """<div style="padding-top: 10px; padding-bottom: 0px">
@@ -306,17 +333,25 @@ class CodeSettings(ipw.VBox):
     # In the order of priority, we will select the default ORCA code from these
     # First, we try to use SLURM on local machine, if available
     _DEFAULT_ORCA_CODES = ("orca@slurm", "orca@localhost")
+    _DEFAULT_ORCA_PLOT_CODES = ("orca_plot@localhost")
 
     def __init__(self, **kwargs):
         self.orca = ComputationalResourcesWidget(
             default_calc_job_plugin="orca.orca",
             description="ORCA program",
         )
+
+        self.orca_plot = ComputationalResourcesWidget(
+            default_calc_job_plugin="orca.orca_plot",
+            description="ORCA plotting program",
+        )
+
         super().__init__(
             children=[
                 self.codes_title,
                 # self.codes_help,
                 self.orca,
+                self.orca_plot,
             ],
             **kwargs,
         )
@@ -341,6 +376,21 @@ def _set_default_codes(self, _=None):
 
         if not self.orca.value:
             print("WARNING: ORCA code has not been found locally")
+
+        for code_label in self._DEFAULT_ORCA_PLOT_CODES:
+            try:
+                self.orca_plot.value = load_code(code_label).uuid
+                return
+            except (NotExistent, ValueError):
+                pass
+            except tl.TraitError:
+                # This can happen if one of the code/computers is not configured/enabled or hidden
+                # In practice, this happened to me locally when importing from production DB.
+                # https://github.com/ispg-group/aiidalab-ispg/issues/240
+                pass
+
+        if not self.orca_plot.value:
+            print("WARNING: ORCA_PLOT code has not been found locally")
 
     def reset(self):
         self._set_default_codes()

aiidalab_ispg/nto/parsercalcfunction.py

@@ -0,0 +1,112 @@
+# Modified from https://github.com/radi0sus/orca_st
+
+from aiida.engine import calcfunction
+from aiida.orm import Dict
+
+import re
+
+@calcfunction
+def parse_orca_output(nto_folder, filename="aiida.out", threshold=0, states="all"):    
+
+    #Convert from Aiida nodes to python datatypes if required.
+    filename = filename.value
+    threshold = threshold.value
+    states = states.value
+
+    # global constants
+    state_string_start = 'TD-DFT/TDA EXCITED STATES'                                         #check for states from here
+    state_string_end = 'TD-DFT/TDA-EXCITATION SPECTRA'                                       #stop reading states from here
+    nto_string_start ='NATURAL TRANSITION ORBITALS'                                          #NTOs start here
+    found_nto = False                                                                        #found NTOs in orca.out
+
+    #global lists
+    orblist = list()              #list of orbital -> orbital transition
+    statedict = dict()            #dictionary of states with all transitions
+    selected_statedict = dict()   #dictionary of selected states with all transitions and or with those above the threshold
+
+    #for NTO
+    nto_orblist=list()            #list of orbital -> orbital transition for NTOs
+    nto_statedict=dict()          #dictionary of states with all transitions for NTOs
+
+    #check if threshold is between 0 and 100%, reset if not
+    if threshold:
+        if threshold > 100 or threshold < 0:
+            print("Warning! Threshold out of range. Reset to 0.")
+            threshold=0
+
+    #open a file
+    #check existence
+    try:
+        with nto_folder.open(filename) as input_file:
+            for line in input_file:
+                #start exctract text 
+                if state_string_start in line:
+                    for line in input_file:
+                        #build the state - with several transitions dict: dict[state]=list(orb1 -> orb2, xx%), list(orb3 -> orb4, xx%)
+                        #first the state
+                        if re.search(r"STATE\s{1,}(\d{1,}):",line):
+                            match_state=re.search(r"STATE\s{1,}(\d{1,}):",line)
+                        #transitions here in orblist
+                        elif re.search(r"\d{1,}[a,b]\s{1,}->\s{1,}\d{1,}[a,b]",line):
+                            match_orbs=re.search(r"(\d{1,}[a,b]\s{1,}->\s{1,}\d{1,}[a,b])\s{1,}:\s{1,}(\d{1,}.\d{1,})",line)
+                            orblist.append((match_orbs.group(1).replace("  "," "),match_orbs.group(2)))
+                        #add orblist to statedict and clear orblist for next state
+                        elif re.search(r"^\s*$",line):
+                            if orblist:
+                                statedict[match_state.group(1)] = orblist
+                                orblist=[] 
+
+                        #exit here
+                        elif nto_string_start in line:
+                            break
+                        elif state_string_end in line:
+                            break
+
+                #same for NTOs
+                if re.search(r"FOR STATE\s{1,}(\d{1,})",line):
+                    #found NTO in orca.out
+                    found_nto = True
+                    match_state_nto=re.search(r"FOR STATE\s{1,}(\d{1,})",line)
+                elif re.search(r"\d{1,}[a,b]\s{1,}->\s{1,}\d{1,}[a,b]",line):
+                    match_orbs_nto=re.search(r"(\d{1,}[a,b]\s{1,}->\s{1,}\d{1,}[a,b])\s{1,}: n=\s{1,}(\d{1,}.\d{1,})",line)
+                    nto_orblist.append((match_orbs_nto.group(1).replace(" ","").split("->"),match_orbs_nto.group(2)))
+
+                #add orblist to statedict and clear orblist for next state
+                elif re.search(r"^\s*$",line):
+                    if nto_orblist:
+                        nto_statedict[match_state_nto.group(1)] = nto_orblist
+                        nto_orblist=[] 
+
+    #file not found -> exit here
+    except IOError:
+        print(f"'{filename}'" + " not found")
+        return {}
+
+    #no NTO data in orca.out -> exit here
+    if found_nto == False:
+        print(f"'{nto_string_start}'" + " not found in " + f"'{filename}'")
+        return {}
+
+    #build selected_statedict from statedict with selected states
+    try:
+        if states == 'all':
+            selected_statedict = nto_statedict
+
+        elif re.search(r"\d",states):
+            matchstateslist=(re.findall(r"\d+",states))
+            for elements in matchstateslist:
+                selected_statedict[elements]=nto_statedict[elements]
+
+    except KeyError:
+        print("Warning! State(s) not present. Exit.")
+        return {}
+
+    #remove transitions below threshold from selected_statedict
+    for elements in selected_statedict:
+        transition_list=[]
+        for v in selected_statedict[elements]:
+            if float(v[1])*100 >= threshold:
+                transition_list.append(v)
+        selected_statedict[elements]=transition_list    
+
+    return Dict(selected_statedict)