Merge pull request #91 from davidkastner/correct-molden

Fix ECP in molden

Author: David W. Kastner

qp/analyze/checkup.py

@@ -53,22 +53,27 @@ def check_submit_record(submit_record_path, delete_queued):
     with open(submit_record_path, 'r') as f:
         content = f.read()
 
-        author = extract_author(content)
-        queue_time = "Queue Time:" in content
-        run_start_time = "Run Start Time:" in content
-        run_end_time = "Run End Time:" in content
-
-        # if queue_time and not run_start_time:
-        if run_start_time and not run_end_time:
-            if delete_queued:
-                print(f"Deleting queued job record: {submit_record_path}")
-                os.remove(submit_record_path)
-            return "queue", author
-        elif run_start_time and not run_end_time:
-            return "running", author
-        elif run_end_time:
-            return "done", author
-        
+    author = extract_author(content)
+    queue_time     = "Queue Time:" in content
+    run_start_time = "Run Start Time:" in content
+    run_end_time   = "Run End Time:" in content
+
+    # Queued but never started
+    if queue_time and not run_start_time:
+        if delete_queued:
+            print(f"Deleting queued job record: {submit_record_path}")
+            os.remove(submit_record_path)
+        return "queue", author
+
+    # Started but not finished
+    if run_start_time and not run_end_time:
+        return "running", author
+
+    # Finished
+    if run_end_time:
+        return "done", author
+
+    # Fallback (record present but missing expected markers)
     return "backlog", author
 
 

qp/analyze/molden.py

@@ -1,4 +1,56 @@
 import numpy as np
+import os
+
+VALENCE_ELECTRONS = {
+    'CA': 2, 'CD': 12, 'CO': 9, 'CU': 11, 'FE': 8, 'K': 1,
+    'MN': 7, 'NI': 10, 'PB': 14, 'V': 5, 'ZN': 12
+}
+
+def correct_ecp_charges_inplace(molden_file):
+    """
+    Corrects the nuclear charge in a Molden file for elements with ECPs,
+    overwriting the original file.
+    """
+    try:
+        with open(molden_file, 'r') as f:
+            lines = f.readlines()
+    except FileNotFoundError:
+        print(f"> ERROR: Molden file not found at {molden_file}")
+        return
+
+    modified_lines = []
+    in_atoms_section = False
+    was_modified = False
+    for line in lines:
+        if line.strip() == "[Atoms] Angs":
+            in_atoms_section = True
+            modified_lines.append(line)
+            continue
+        elif in_atoms_section and line.strip().startswith("["):
+            in_atoms_section = False
+        
+        if in_atoms_section:
+            parts = line.strip().split()
+            if len(parts) >= 6:
+                element_symbol = parts[0].upper()
+                # Check if the element needs correction and if it hasn't been corrected already
+                if element_symbol in VALENCE_ELECTRONS and parts[2] != str(VALENCE_ELECTRONS[element_symbol]):
+                    parts[2] = str(VALENCE_ELECTRONS[element_symbol])
+                    line_to_write = '     '.join(parts) + '\n'
+                    modified_lines.append(line_to_write)
+                    was_modified = True
+                else:
+                    modified_lines.append(line)
+            else:
+                modified_lines.append(line)
+        else:
+            modified_lines.append(line)
+    
+    if was_modified:
+        with open(molden_file, 'w') as f:
+            f.writelines(modified_lines)
+        print(f"> Overwrote {os.path.basename(molden_file)} with ECP-corrected nuclear charges.")
+
 
 def translate_to_origin(coordinates, center_of_mass):
     """Translates the molecule so the center of mass is at the origin."""
@@ -92,4 +144,4 @@ def center_molden(molden_file, com):
     process_molden(molden_file, output_molden, com)
     print(f"> Centered coordinates written to: {output_molden}")
 
-    return output_molden
+    return output_molden

qp/analyze/multiwfn.py

@@ -54,38 +54,34 @@ def iterate_qm_output(pdb_all, method, base_output_dir, multiwfn_path, settings_
                 # Copy the atmrad dir and settings.ini into the current working directory
                 atmrad_dest_path = os.path.join(scr_dir_path, 'atmrad/')
                 if not os.path.exists(atmrad_dest_path):
-                    shutil.copytree(atmrad_path, atmrad_dest_path, dirs_exist_ok=True)
+                    shutil.copytree(atmrad_path, atmrad_dest_path)
 
                 settings_ini_dest_path = os.path.join(scr_dir_path, 'settings.ini')
-                shutil.copy(settings_ini_path, settings_ini_dest_path)
+                # MODIFICATION: Add a check to avoid re-copying settings.ini if it already exists
+                if not os.path.exists(settings_ini_dest_path):
+                    shutil.copy(settings_ini_path, settings_ini_dest_path)
+                
                 cpu_count = get_cpu_count(settings_ini_dest_path)
-                print(f"> Using {cpu_count} threads for Multiwfn")
+                print(f"> Using {cpu_count} threads for Multiwfn", flush=True)
 
                 if os.path.exists(scr_dir_path):
+                    # Move into the QM scr directory
                     original_dir = os.getcwd()
                     os.chdir(scr_dir_path)
-                    try:
-                        scr_dir_path = os.getcwd()
-
-                        molden_files = glob.glob(f'{chain_dir}.molden')
-                        if molden_files:
-                            molden_file = os.path.basename(molden_files[0])
-                            try:
-                                task_function(scr_dir_path, molden_file, multiwfn_path, charge_scheme)
-                            except Exception as e:
-                                print(f"> ERROR: {current_pdb_dir} {chain_dir} failed: {e}. Skipping.")
-                        else:
-                            print(f"> WARNING: No molden file in {scr_dir_path}")
-                    finally:
-                        shutil.rmtree('atmrad/', ignore_errors=True)
-                        try:
-                            os.remove('settings.ini')
-                        except FileNotFoundError:
-                            pass
-                        os.chdir(original_dir)
+                    scr_dir_path = os.getcwd()
+
+                    molden_files = glob.glob(f'{chain_dir}.molden')
+                    if molden_files:
+                        molden_file = os.path.basename(molden_files[0])
+                        task_function(scr_dir_path, molden_file, multiwfn_path, charge_scheme)
+                        # MODIFICATION: Comment out the deletion lines to prevent race conditions
+                        # shutil.rmtree('atmrad/')  # Delete the atmrad directory when done
+                        # os.remove('settings.ini')  # Delete the settings.ini file when done
+                    else:
+                        print(f"> WARNING: No molden file in {scr_dir_path}", flush=True)
+                    os.chdir(original_dir)
                 else:
-                    print(f"> WARNING: No scr directory in {method_dir_path}.")
-
+                    print(f"> WARNING: No scr directory in {method_dir_path}.", flush=True)
         else:
             continue
 
@@ -187,7 +183,7 @@ def get_coc(scr_dir_path, molden_file, fragment_atom_indices, selected_charge_sc
 
     # Check if the .chg file exists
     if not os.path.exists(charge_file):
-        print(f"> Charge file {charge_file} not found. Generating it using charge_scheme function.")
+        print(f"> Charge file {charge_file} not found. Generating it using charge_scheme function.", flush=True)
         # Generate the .chg file using the charge_scheme function
         charge_scheme(scr_dir_path, molden_file, multiwfn_path, selected_charge_scheme)
 
@@ -293,9 +289,9 @@ def get_cpu_count(settings_ini_dest_path):
                     # Extract the first matching group (the number after 'nthreads=')
                     return int(match.group(1))
     except FileNotFoundError:
-        print(f"Error: The file {settings_ini_dest_path} was not found.")
+        print(f"Error: The file {settings_ini_dest_path} was not found.", flush=True)
     except ValueError:
-        print(f"Error: Unable to parse 'nthreads' in {settings_ini_dest_path}.")
+        print(f"Error: Unable to parse 'nthreads' in {settings_ini_dest_path}.", flush=True)
 
     # Default to 1 if nthreads isn't found or there is an error
     return 1
@@ -364,7 +360,7 @@ def charge_scheme(scr_dir_path, molden_file, multiwfn_path, charge_scheme):
 
     # If the charge file already exists, skip
     if os.path.exists(new_charge_file):
-        print(f"> Charge file {new_charge_file} already exists. Skipping {charge_scheme}...")
+        print(f"> Charge file {new_charge_file} already exists. Skipping {charge_scheme}...", flush=True)
         return
 
     # Run Multiwfn with the user-selected charge scheme
@@ -376,9 +372,9 @@ def charge_scheme(scr_dir_path, molden_file, multiwfn_path, charge_scheme):
     charge_output_file = Path(f"{molden_base_name}.chg")
     if charge_output_file.exists():
         charge_output_file.rename(new_charge_file)
-        print(f"> {pdb_name.upper()} {chain_name} {charge_scheme} scheme computed and saved to {new_charge_file}")
+        print(f"> {pdb_name.upper()} {chain_name} {charge_scheme} scheme computed and saved to {new_charge_file}", flush=True)
     else:
-        print(f"> ERROR: Multiwfn failed for {pdb_name} {chain_name} using {charge_scheme} scheme.")
+        print(f"> ERROR: Multiwfn failed for {pdb_name} {chain_name} using {charge_scheme} scheme.", flush=True)
 
 
 
@@ -434,4 +430,3 @@ def calc_dipole(scr_dir_path, molden_file, multiwfn_path, charge_scheme):
                 csv_writer.writerow(['Molecular dipole moment'] + dipole_a_u + [magnitude_a_u])
 
     print(f"> {pdb_name.upper()} {chain_name} dipole computed using center of mass\n", flush=True)
-