Add files via upload

### 3. Integration of Multiple PES Information and Model Function Optimization (BITSS, etc.)

#### Overview
A feature to perform geometry optimization on an "effective potential" constructed by combining energies or gradients from multiple electronic states (PES).
Specifically, version 1.20.3 implements the Binary-Image Transition State Search (BITSS) method.

#### Implementation Details
*   **Independent Calculator Instances:** `ModelFunctionHandler` creates independent directories and calculators for State1 and State2 to prevent interference between states.
*   **BITSS Implementation:**
    *   **6N-Dimensional Expansion:** Concatenates two structures (Image 1 & 2) and treats them as a $6N \times 6N$ Hessian and a $2N$ atom system.
    *   **Second Derivative (Hessian):** Mathematically derives the second derivatives for distance and energy equality constraint terms and implements them as the `calc_hess` method.

#### Command Line Arguments (`-mf`)
Use the `-mf` argument to specify the type of model function and accompanying parameters (the file path of the reference structure in the case of BITSS).

#### Usage Example:
```bash
# Search for a pathway connecting two points using the BITSS method, with target.xyz as the target structure.
multioptpy start.xyz -mf bitss target.xyz
```

```bash
# Minimum energy seam of crossing (MECI) using Seam Model Function (between charge 0, multiplicity 1 and 3 states)
multioptpy input.xyz -mf opt_meci 0 3 -elec 0 -spin 3
```

Author: ss0832

multioptpy/ModelFunction/opt_meci.py

@@ -2,49 +2,116 @@
 
 class OptMECI:
     def __init__(self):
-        # ref.:https://doi.org/10.1021/ct1000268
+        # ref.: J. Am. Chem. Soc. 2015, 137, 3433-3445
+        # MECI optimization using GP method with Branching Plane Updating (BPU) 
 
-        self.switch_threshold = 5e-4
-        self.alpha = 1e-3
-        self.approx_cdv_vec = None
-        self.prev_dgv_vec = None
-        self.dgv_vec = None
-        
+        self.approx_cdv_vec = None # Represents 'y' vector (orthogonal to dgv inside BP)
+        self.prev_dgv_vec = None   # Represents 'x_{k-1}'
+        self.prev_y_vec = None     # Represents 'y_{k-1}'
         return
 
     def calc_energy(self, energy_1, energy_2):
         tot_energy = (energy_1 + energy_2) / 2.0 
         print("energy_1:", energy_1, "hartree")
         print("energy_2:", energy_2, "hartree")
-        print("energy_1 - energy_2:", abs(energy_1 - energy_2), "hartree")
+        print("|energy_1 - energy_2|:", abs(energy_1 - energy_2), "hartree")
         return tot_energy
 
     def calc_grad(self, energy_1, energy_2, grad_1, grad_2):
-        if self.approx_cdv_vec is None:
-            self.approx_cdv_vec = np.ones((len(grad_1)*3, 1))
-
-        delta_grad = grad_1 - grad_2
-        dgv_vec = delta_grad / np.linalg.norm(delta_grad)
-        dgv_vec = dgv_vec.reshape(-1, 1)
+        # Reshape inputs
+        grad_1_flat = grad_1.reshape(-1, 1)
+        grad_2_flat = grad_2.reshape(-1, 1)
 
-        if self.prev_dgv_vec is None:
-            self.prev_dgv_vec = dgv_vec
+        # 1. Calculate Difference Gradient Vector (x_k) 
+        delta_grad = grad_1_flat - grad_2_flat
+        norm_delta_grad = np.linalg.norm(delta_grad)
+        if norm_delta_grad < 1e-8:
+            dgv_vec = np.zeros_like(delta_grad) # Avoid division by zero
+        else:
+            dgv_vec = delta_grad / norm_delta_grad # x_k
 
-        self.approx_cdv_vec = (np.dot(self.approx_cdv_vec.T, dgv_vec) * self.prev_dgv_vec -1 * np.dot(self.prev_dgv_vec.T, dgv_vec) * self.approx_cdv_vec) / np.sqrt(np.dot(self.approx_cdv_vec.T, dgv_vec) ** 2 + np.dot(self.prev_dgv_vec.T, dgv_vec) ** 2) 
+        # 2. Determine Approximate Coupling Vector (y_k) using BPU 
+        if self.prev_dgv_vec is None:
+            # Initialization Step 
+            # "A plane made of x0 and the mean energy gradient vector was used as an initial BP."
+            mean_grad = 0.5 * (grad_1_flat + grad_2_flat)
+            
+            # Project mean_grad to be orthogonal to dgv_vec (Gram-Schmidt)
+            overlap = np.dot(mean_grad.T, dgv_vec)
+            ortho_vec = mean_grad - overlap * dgv_vec
+            
+            norm_ortho = np.linalg.norm(ortho_vec)
+            if norm_ortho < 1e-8:
+                 # Fallback if mean grad is parallel to diff grad (unlikely)
+                 ortho_vec = np.random.rand(*dgv_vec.shape)
+                 ortho_vec = ortho_vec - np.dot(ortho_vec.T, dgv_vec) * dgv_vec
+                 norm_ortho = np.linalg.norm(ortho_vec)
 
-        P_matrix = np.eye((len(dgv_vec))) -1 * np.dot(dgv_vec, dgv_vec.T) -1 * np.dot(self.approx_cdv_vec, self.approx_cdv_vec.T)
-        P_matrix = 0.5 * (P_matrix + P_matrix.T)
-        gp_grad =  2 * (energy_1 - energy_2) * dgv_vec + np.dot(P_matrix, 0.5 * (grad_1.reshape(-1, 1) + grad_2.reshape(-1, 1)))
-        
-        self.prev_dgv_vec = dgv_vec
-        gp_grad = gp_grad.reshape(len(grad_1), 3)
-        return gp_grad
+            self.approx_cdv_vec = ortho_vec / norm_ortho # Initial y_0
+            
+        else:
+            # Update Step using Eq 4 
+            # y_k = [ (y_{k-1}.x_k) * x_{k-1} - (x_{k-1}.x_k) * y_{k-1} ] / normalization
+            
+            x_k = dgv_vec
+            x_prev = self.prev_dgv_vec
+            y_prev = self.prev_y_vec
+            
+            dot_yx = np.dot(y_prev.T, x_k)
+            dot_xx = np.dot(x_prev.T, x_k)
+            
+            numerator = dot_yx * x_prev - dot_xx * y_prev
+            norm_num = np.linalg.norm(numerator)
+            
+            if norm_num < 1e-8:
+                 # If x_k didn't change much, keep y_prev orthogonalized to x_k
+                 numerator = y_prev - np.dot(y_prev.T, x_k) * x_k
+                 norm_num = np.linalg.norm(numerator)
 
+            self.approx_cdv_vec = numerator / norm_num # y_k
+
+        # Store vectors for next step
+        self.prev_dgv_vec = dgv_vec.copy()
+        self.prev_y_vec = self.approx_cdv_vec.copy()
+        
+        # 3. Construct Projection Matrix P for MECI 
+        # Projects out BOTH dgv (x) and approx_cdv (y) directions
+        P_matrix = np.eye(len(dgv_vec)) \
+                   - np.dot(dgv_vec, dgv_vec.T) \
+                   - np.dot(self.approx_cdv_vec, self.approx_cdv_vec.T)
+        
+        # 4. Compose Gradient Projection (GP) Gradient
+        # Force to reduce energy gap: 2 * (E1 - E2) * dgv
+        # Force to minimize mean energy on intersection space (N-2 dim): P * mean_grad
+        mean_grad = 0.5 * (grad_1_flat + grad_2_flat)
+        
+        gap_force = 2.0 * (energy_1 - energy_2) * dgv_vec
+        seam_force = np.dot(P_matrix, mean_grad)
+        
+        gp_grad = gap_force + seam_force
+        
+        return gp_grad.reshape(len(grad_1), 3)
 
     def calc_hess(self, hess_1, hess_2):
+        # Approximate Hessian for GP method
+        # Projects the mean Hessian onto the intersection space
         mean_hess = 0.5 * (hess_1 + hess_2)
-        P_matrix = np.eye((len(self.prev_dgv_vec))) -1 * np.dot(self.prev_dgv_vec, self.prev_dgv_vec.T) -1 * np.dot(self.approx_cdv_vec, self.approx_cdv_vec.T)
-        P_matrix = 0.5 * (P_matrix + P_matrix.T)
 
-        gp_hess = np.dot(P_matrix, np.dot(mean_hess, P_matrix))
+        # Need current P_matrix. Reconstruct it from stored vectors.
+        if self.approx_cdv_vec is None or self.prev_dgv_vec is None:
+             # Should not happen if calc_grad is called first
+             return mean_hess
+             
+        dgv_vec = self.prev_dgv_vec
+        cdv_vec = self.approx_cdv_vec
+        
+        P_matrix = np.eye(len(dgv_vec)) \
+                   - np.dot(dgv_vec, dgv_vec.T) \
+                   - np.dot(cdv_vec, cdv_vec.T)
+                   
+        # Projected Mean Hessian + Gap Penalty Curvature
+        proj_hess = np.dot(P_matrix, np.dot(mean_hess, P_matrix))
+        gap_curvature = 2.0 * np.dot(dgv_vec, dgv_vec.T)
+        
+        gp_hess = proj_hess + gap_curvature
         return gp_hess

multioptpy/ModelFunction/opt_mesx.py

@@ -2,46 +2,78 @@
 
 class OptMESX:
     def __init__(self):
-        #ref.: Chemical Physics Letters 674 (2017) 141-145
-       
-        self.switch_threshold = 5e-4
-        self.alpha = 1e-3
+        # ref.: J. Am. Chem. Soc. 2015, 137, 3433-3445 
+        # MESX optimization using Gradient Projection (GP) method.
+        # Only the difference gradient vector (DG or f) is projected out.
         return
 
     def calc_energy(self, energy_1, energy_2):
+        # The objective is to minimize the mean energy on the seam.
         tot_energy = (energy_1 + energy_2) / 2.0 
         print("energy_1:", energy_1, "hartree")
         print("energy_2:", energy_2, "hartree")
         print("|energy_1 - energy_2|:", abs(energy_1 - energy_2), "hartree")
         return tot_energy
 
     def calc_grad(self, energy_1, energy_2, grad_1, grad_2):
+        # 1. Calculate Difference Gradient Vector (DGV) / f vector
         delta_grad = grad_1 - grad_2
         norm_delta_grad = np.linalg.norm(delta_grad)
+        
         if norm_delta_grad < 1e-8:
             dgv_vec = np.zeros_like(delta_grad)
         else:
             dgv_vec = delta_grad / norm_delta_grad
+        
+        # Ensure correct shape for matrix operations
         dgv_vec = dgv_vec.reshape(-1, 1)
+        grad_1 = grad_1.reshape(-1, 1)
+        grad_2 = grad_2.reshape(-1, 1)
+        
+        # 2. Define Projection Matrix P for MESX
+        # Projects out the component along the difference vector (degenerate lifting direction)
+        # P = I - v * v.T
+        P_matrix = np.eye(len(dgv_vec)) - np.dot(dgv_vec, dgv_vec.T)
 
-        P_matrix = np.eye((len(dgv_vec))) -1 * np.dot(dgv_vec, dgv_vec.T)
-        P_matrix = 0.5 * (P_matrix + P_matrix.T)
-        gp_grad =  2 * (energy_1 - energy_2) * dgv_vec + np.dot(P_matrix, 0.5 * (grad_1.reshape(-1, 1) + grad_2.reshape(-1, 1)))
+        # 3. Calculate Mean Gradient
+        mean_grad = 0.5 * (grad_1 + grad_2)
 
-        gp_grad = gp_grad.reshape(len(grad_1), 3)
-        return gp_grad
+        # 4. Compose Gradient Projection (GP) Gradient 
+        # Force to reduce energy gap: 2 * (E1 - E2) * dgv
+        # Force to minimize mean energy on seam: P * mean_grad
+        gap_force = 2.0 * (energy_1 - energy_2) * dgv_vec
+        seam_force = np.dot(P_matrix, mean_grad)
+        
+        gp_grad = gap_force + seam_force
+        
+        return gp_grad.reshape(-1, 3)
 
     def calc_hess(self, grad_1, grad_2, hess_1, hess_2):
+        # Approximate Hessian for GP method
         delta_grad = grad_1 - grad_2
         norm_delta_grad = np.linalg.norm(delta_grad)
+        
         if norm_delta_grad < 1e-8:
             dgv_vec = np.zeros_like(delta_grad)
         else:
             dgv_vec = delta_grad / norm_delta_grad
-        delta_grad = delta_grad.reshape(-1, 1)
+            
         dgv_vec = dgv_vec.reshape(-1, 1)
-        P_matrix = np.eye((len(dgv_vec))) -1 * np.dot(dgv_vec, dgv_vec.T)
-        P_matrix = 0.5 * (P_matrix + P_matrix.T)
-        gp_hess = 2.0 * np.dot(delta_grad, dgv_vec.T) + np.dot(P_matrix, 0.5 * (hess_1 + hess_2))
-        return gp_hess
-        
+        
+        # Projection Matrix
+        P_matrix = np.eye(len(dgv_vec)) - np.dot(dgv_vec, dgv_vec.T)
+        
+        # Mean Hessian
+        mean_hess = 0.5 * (hess_1 + hess_2)
+        
+        # Projected Mean Hessian
+        # This describes curvature along the seam.
+        proj_hess = np.dot(P_matrix, np.dot(mean_hess, P_matrix))
+        
+        # Gap Curvature (Penalty term approximation)
+        # Adds large curvature along the difference vector to enforce the gap constraint strongly.
+        gap_curvature = 2.0 * np.dot(dgv_vec, dgv_vec.T)
+        
+        gp_hess = proj_hess + gap_curvature
+        
+        return gp_hess

multioptpy/ModelFunction/opt_mesx_2.py

@@ -2,48 +2,65 @@
 
 class OptMESX2:
     def __init__(self):
-        #ref.: Theor Chem Acc 99, 95–99 (1998)
-       
-        self.switch_threshold = 5e-4
-        self.alpha = 1e-3
+        # ref.: Theor Chem Acc 99, 95–99 (1998)
+        # This reference describes the Gradient Projection method.
+        # The implementation has been corrected to follow the standard GP formulation
+        # as described in J. Am. Chem. Soc. 2015, 137, 3433-3445 .
         return
 
     def calc_energy(self, energy_1, energy_2):
         tot_energy = (energy_1 + energy_2) / 2.0 
         print("energy_1:", energy_1, "hartree")
         print("energy_2:", energy_2, "hartree")
-        print("energy_1 - energy_2:", abs(energy_1 - energy_2), "hartree")
+        print("|energy_1 - energy_2|:", abs(energy_1 - energy_2), "hartree")
         return tot_energy
 
     def calc_grad(self, energy_1, energy_2, grad_1, grad_2):
-        grad_1 = grad_1.reshape(-1, 1)
-        grad_2 = grad_2.reshape(-1, 1)
-        
+        # 1. Difference Vector (normalized)
         delta_grad = grad_1 - grad_2
         norm_delta_grad = np.linalg.norm(delta_grad)
+        
         if norm_delta_grad < 1e-8:
-            projection = np.zeros_like(delta_grad)
+            dgv_vec = np.zeros_like(delta_grad)
         else:
-            projection = np.sum(grad_1 * delta_grad) / norm_delta_grad
+            dgv_vec = delta_grad / norm_delta_grad
 
-        parallel = grad_1 - delta_grad * projection / norm_delta_grad
+        dgv_vec = dgv_vec.reshape(-1, 1)
+        grad_1_flat = grad_1.reshape(-1, 1)
+        grad_2_flat = grad_2.reshape(-1, 1)
+
+        # 2. Projection Matrix (P = I - v v^T)
+        P_matrix = np.eye(len(dgv_vec)) - np.dot(dgv_vec, dgv_vec.T)
 
-        gp_grad = (energy_1 - energy_2) * 140 * delta_grad + 1.0 * parallel
+        # 3. Mean Gradient
+        mean_grad = 0.5 * (grad_1_flat + grad_2_flat)
 
-        gp_grad = gp_grad.reshape(-1, 3)
-        return gp_grad
+        # 4. Recomposed Gradient 
+        # Replaces the arbitrary '140' factor with the analytical gap force 2(E1-E2)
+        gap_force = 2.0 * (energy_1 - energy_2) * dgv_vec
+        seam_force = np.dot(P_matrix, mean_grad)
+        
+        gp_grad = gap_force + seam_force
+        
+        return gp_grad.reshape(-1, 3)
 
     def calc_hess(self, grad_1, grad_2, hess_1, hess_2):
+        # Robust Hessian construction for GP
         delta_grad = grad_1 - grad_2
         norm_delta_grad = np.linalg.norm(delta_grad)
         if norm_delta_grad < 1e-8:
             dgv_vec = np.zeros_like(delta_grad)
         else:
             dgv_vec = delta_grad / norm_delta_grad
-        delta_grad = delta_grad.reshape(-1, 1)
+            
         dgv_vec = dgv_vec.reshape(-1, 1)
-        P_matrix = np.eye((len(dgv_vec))) -1 * np.dot(dgv_vec, dgv_vec.T)
-        P_matrix = 0.5 * (P_matrix + P_matrix.T)
 
-        gp_hess = 2.0 * np.dot(delta_grad, dgv_vec.T) + np.dot(P_matrix, 0.5 * (hess_1 + hess_2))
+        P_matrix = np.eye(len(dgv_vec)) - np.dot(dgv_vec, dgv_vec.T)
+        mean_hess = 0.5 * (hess_1 + hess_2)
+        
+        # Projected Mean Hessian + Gap Curvature
+        proj_hess = np.dot(P_matrix, np.dot(mean_hess, P_matrix))
+        gap_curvature = 2.0 * np.dot(dgv_vec, dgv_vec.T)
+        
+        gp_hess = proj_hess + gap_curvature
         return gp_hess

multioptpy/fileio.py

@@ -399,9 +399,9 @@ def print_geometry_list(self, new_geometry, element_list, electric_charge_and_mu
             for num, geometry in enumerate(new_geometry):    
                 element = element_list[num]
                 print(f"{element:2}   {float(geometry[0]):>17.12f}   {float(geometry[1]):>17.12f}   {float(geometry[2]):>17.12f}")
-        
+            print("\n")
+            
         geometry_list = [[electric_charge_and_multiplicity, *formatted_geometries]]
-        print("\n")
 
         return geometry_list
 

multioptpy/optimization.py

@@ -562,11 +562,11 @@ def _run_calc(self, calc_inst, geom, elems, chg_mult, label, iter_idx):
 
     def finalize_bitss_trajectory(self):
         if not self.is_bitss or not self.bitss_geom1_history: return
-        filename = os.path.join(self.base_dir, f"{self.file_io.NOEXT_START_FILE}_bitss_path.xyz")
+        filename = os.path.join(self.base_dir, f"{self.file_io.NOEXT_START_FILE}_traj.xyz")
         full_seq = self.bitss_geom1_history + self.bitss_geom2_history[::-1]
         with open(filename, 'w') as f:
             for s, g in enumerate(full_seq):
-                f.write(f"{len(g)}\nStep {s}\n")
+                f.write(f"{len(g)}\nBITSS_Step {s}\n")
                 for i, atom in enumerate(g):
                     f.write(f"{self.single_element_list[i]:2s} {atom[0]:12.8f} {atom[1]:12.8f} {atom[2]:12.8f}\n")