Skip to content

Commit bfcf788

Browse files
ss0832ss0832
authored
Add files via upload
1 parent f44482c commit bfcf788

1 file changed

Lines changed: 204 additions & 14 deletions

File tree

  • multioptpy/Calculator/ase_tools

‎multioptpy/Calculator/ase_tools/orca.py‎

Lines changed: 204 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,22 +1,212 @@
11
import os
2+
import shutil
3+
import subprocess
4+
import numpy as np
5+
from ase.calculators.orca import ORCA, OrcaProfile
6+
from ase.data import atomic_numbers
27

38
class ASE_ORCA:
9+
TARGET_ORCA_VERSION = '6.1.0'
10+
411
def __init__(self, **kwargs):
512
self.atom_obj = kwargs.get('atom_obj', None)
613
self.electric_charge_and_multiplicity = kwargs.get('electric_charge_and_multiplicity', None)
7-
self.input_file = kwargs.get('input_file', None)
14+
15+
# NOTE: self.input_file is updated in _setup_calculator to enforce 'orca.inp' in CWD.
16+
raw_input_file = kwargs.get('input_file', 'orca.inp')
17+
self.input_file = os.path.abspath(raw_input_file).replace('\\', '/')
18+
819
self.orca_path = kwargs.get('orca_path', None)
9-
self.functional = kwargs.get('functional', None)
10-
self.basis_set = kwargs.get('basis_set', None)
11-
20+
self.functional = kwargs.get('functional', 'B3LYP')
21+
self.basis_set = kwargs.get('basis_set', 'def2-SVP')
22+
23+
# Optional ORCA input blocks (raw string)
24+
self.orca_blocks = kwargs.get('orca_blocks', '')
25+
26+
# Auto-fix for unsupported Pople basis sets on heavy elements
27+
self.auto_fix_basis = kwargs.get('auto_fix_basis', True)
28+
self.heavy_atom_basis = kwargs.get('heavy_atom_basis', 'def2-SVP')
29+
30+
def _resolve_orca_exe(self, provided_path):
31+
if not provided_path:
32+
return None
33+
clean_path = os.path.normpath(provided_path)
34+
if os.path.isdir(clean_path):
35+
candidate = os.path.join(clean_path, 'orca.exe') if os.name == 'nt' else os.path.join(clean_path, 'orca')
36+
else:
37+
candidate = clean_path
38+
39+
if os.path.exists(candidate) and os.path.isfile(candidate):
40+
return os.path.abspath(candidate)
41+
42+
resolved = shutil.which(candidate)
43+
if resolved and os.path.exists(resolved):
44+
return os.path.abspath(resolved)
45+
46+
raise FileNotFoundError(f"Cannot locate ORCA executable at: {candidate}")
47+
48+
def _is_pople_basis(self, basis_name):
49+
if not basis_name: return False
50+
b = basis_name.strip().lower()
51+
return b.startswith("6-31") or b.startswith("6-311")
52+
53+
def _get_heavy_elements_for_pople(self):
54+
if self.atom_obj is None: return []
55+
symbols = self.atom_obj.get_chemical_symbols()
56+
return sorted({s for s in symbols if atomic_numbers.get(s, 0) > 30})
57+
58+
def _build_orca_blocks(self):
59+
blocks = (self.orca_blocks or "").strip()
60+
heavy_elements = []
61+
if self._is_pople_basis(self.basis_set):
62+
heavy_elements = self._get_heavy_elements_for_pople()
63+
64+
if heavy_elements:
65+
if not self.auto_fix_basis:
66+
raise ValueError("Unsupported Pople basis for heavy elements.")
67+
68+
basis_lines = ["%basis"]
69+
for elem in heavy_elements:
70+
# Per user instruction: No ECP, just GTO
71+
basis_lines.append(f' NewGTO {elem} "{self.heavy_atom_basis}"')
72+
basis_lines.append("end")
73+
74+
if blocks:
75+
blocks = blocks + "\n" + "\n".join(basis_lines)
76+
else:
77+
blocks = "\n".join(basis_lines)
78+
79+
return blocks if blocks else ""
80+
81+
def _print_orca_output_on_error(self):
82+
"""Helper to print ORCA output file content if it exists."""
83+
if not hasattr(self, 'input_file'): return
84+
85+
out_file = os.path.splitext(self.input_file)[0] + ".out"
86+
if os.path.exists(out_file):
87+
print(f"\n--- ORCA OUTPUT ({out_file}) ---")
88+
try:
89+
with open(out_file, 'r', encoding='utf-8', errors='replace') as f:
90+
content = f.read()
91+
print(content[-3000:] if len(content) > 3000 else content)
92+
except Exception as e:
93+
print(f"Failed to read output file: {e}")
94+
print("--- END ORCA OUTPUT ---\n")
95+
else:
96+
print(f"\n--- ORCA OUTPUT NOT FOUND ({out_file}) ---\n")
97+
98+
def _setup_calculator(self, task_keyword):
99+
# Force usage of Current Working Directory + "orca.inp"
100+
cwd = os.getcwd()
101+
label_path = os.path.join(cwd, 'orca').replace('\\', '/')
102+
self.input_file = label_path + '.inp'
103+
104+
print(f"DEBUG: ASE Label Path (Forced): {label_path}")
105+
106+
simple_input = f"{self.functional} {self.basis_set} {task_keyword}"
107+
108+
profile_obj = None
109+
if self.orca_path:
110+
real_exe_path = self._resolve_orca_exe(self.orca_path)
111+
orca_dir = os.path.dirname(real_exe_path)
112+
path_env = os.environ.get('PATH', '')
113+
if orca_dir not in path_env:
114+
os.environ['PATH'] = orca_dir + os.pathsep + path_env
115+
116+
ase_safe_path = real_exe_path.replace('\\', '/')
117+
profile_obj = OrcaProfile(ase_safe_path)
118+
print(f"DEBUG: ORCA Executable: {ase_safe_path}")
119+
120+
orca_blocks = self._build_orca_blocks()
121+
122+
calc = ORCA(
123+
label=label_path,
124+
profile=profile_obj,
125+
charge=int(self.electric_charge_and_multiplicity[0]),
126+
mult=int(self.electric_charge_and_multiplicity[1]),
127+
orcasimpleinput=simple_input,
128+
orcablocks=orca_blocks
129+
)
130+
131+
self.atom_obj.calc = calc
132+
return self.atom_obj
133+
12134
def run(self):
13-
14-
from ase.calculators.orca import ORCA
15-
input_dir = os.path.dirname(self.input_file)
16-
self.atom_obj.calc = ORCA(label=input_dir,
17-
profile=self.orca_path,
18-
charge=int(self.electric_charge_and_multiplicity[0]),
19-
mult=int(self.electric_charge_and_multiplicity[1]),
20-
orcasimpleinput=self.functional+' '+self.basis_set)
21-
#orcablocks='%pal nprocs 16 end')
22-
return self.atom_obj
135+
self._setup_calculator("EnGrad")
136+
print(f"--- Starting Gradient Calculation (ORCA {self.TARGET_ORCA_VERSION}) ---")
137+
try:
138+
forces = self.atom_obj.get_forces()
139+
potential_energy = self.atom_obj.get_potential_energy()
140+
print("Gradient calculation completed.")
141+
return forces, potential_energy
142+
except subprocess.CalledProcessError as e:
143+
print(f"CRITICAL: ORCA execution failed with exit code {e.returncode}")
144+
self._print_orca_output_on_error()
145+
raise e
146+
except Exception as e:
147+
print(f"CRITICAL: An unexpected error occurred: {e}")
148+
self._print_orca_output_on_error()
149+
raise e
150+
151+
def run_frequency_analysis(self):
152+
print(f"--- Starting Frequency Calculation (ORCA {self.TARGET_ORCA_VERSION}) ---")
153+
self._setup_calculator("Freq")
154+
try:
155+
self.atom_obj.get_potential_energy()
156+
print("Frequency calculation completed.")
157+
# Use self.input_file to construct hess_path instead of relying on calc.label
158+
hess_path = os.path.splitext(self.input_file)[0] + ".hess"
159+
return hess_path
160+
except subprocess.CalledProcessError as e:
161+
print(f"CRITICAL: ORCA Frequency execution failed with exit code {e.returncode}")
162+
self._print_orca_output_on_error()
163+
raise e
164+
except Exception as e:
165+
print(f"CRITICAL: An unexpected error occurred during frequency analysis: {e}")
166+
self._print_orca_output_on_error()
167+
raise e
168+
169+
def get_hessian_matrix(self, hess_file_path=None):
170+
if hess_file_path is None:
171+
# Default to orca.hess in the same dir as input_file
172+
input_dir = os.path.dirname(self.input_file)
173+
hess_file_path = os.path.join(input_dir, 'orca.hess')
174+
175+
if not os.path.exists(hess_file_path):
176+
raise FileNotFoundError(f"Hessian file not found: {hess_file_path}")
177+
178+
print(f"Reading Hessian from: {hess_file_path}")
179+
180+
with open(hess_file_path, 'r', encoding='utf-8') as f:
181+
lines = f.readlines()
182+
183+
hessian_matrix = None
184+
iterator = iter(lines)
185+
for line in iterator:
186+
if "$hessian" in line:
187+
dim_line = next(iterator).strip().split()
188+
# Parse dimensions: square matrix usually provides just one dimension
189+
if len(dim_line) == 1:
190+
n_rows = int(dim_line[0])
191+
n_cols = n_rows
192+
else:
193+
n_rows, n_cols = map(int, dim_line[:2])
194+
195+
hessian_matrix = np.zeros((n_rows, n_cols))
196+
col_pointer = 0
197+
while col_pointer < n_cols:
198+
header = next(iterator).strip()
199+
if not header: break
200+
col_indices = [int(c) for c in header.split()]
201+
for r in range(n_rows):
202+
row_data = next(iterator).strip().split()
203+
row_idx = int(row_data[0])
204+
values = [float(x) for x in row_data[1:]]
205+
for i, val in enumerate(values):
206+
hessian_matrix[row_idx, col_indices[i]] = val
207+
col_pointer += len(col_indices)
208+
break
209+
210+
if hessian_matrix is None:
211+
raise ValueError("Could not find $hessian block in the file.")
212+
return hessian_matrix

0 commit comments

Comments
 (0)