Add files via upload

- new option to align image distances (-aid N, nebmain.py)

Author: ss0832

multioptpy/interface.py

@@ -197,7 +197,9 @@ def nebparser(parser):
     parser.add_argument("-ewbneb", "--EWBNEB", action='store_true', help='Energy-weighted Nudged elastic band method')
 
     parser.add_argument("-idpp", "--use_image_dependent_pair_potential", action='store_true', help='use image dependent pair potential (IDPP) method (ref. arXiv:1406.1512v1)')
-    parser.add_argument("-ad", "--align_distances", type=int, default=0, help='distribute images at equal intervals linearly')
+    parser.add_argument("-ad", "--align_distances", type=int, default=0, help='distribute images at equal intervals on the reaction coordinate')
+    parser.add_argument("-aid", "--align_image_distances", type=int, default=0, help='distribute images at equal intervals between images')
+
     parser.add_argument("-nd", "--node_distance", type=float, default=None, help='distribute images at equal intervals linearly based ont specific distance (ex.) [distance (ang.)] (default: None)')
     parser.add_argument("-p", "--partition",  type=int, default='0', help='number of nodes')
     parser.add_argument("-core", "--N_THREAD",  type=int, default='8', help='threads')

multioptpy/neb.py

@@ -104,6 +104,7 @@ def __init__(self, args):
         # Flags
         self.IDPP_flag = args.use_image_dependent_pair_potential
         self.align_distances = args.align_distances
+        self.align_image_distances = args.align_image_distances
         self.excited_state = args.excited_state
         self.unrestrict = args.unrestrict
         self.save_pict = args.save_pict
@@ -465,7 +466,12 @@ def execute(self):
                 tmp_new_geometry = distribute_geometry(np.array(new_geometry))
                 for k in range(len(new_geometry)):
                     new_geometry[k] = copy.copy(tmp_new_geometry[k])
-            
+            if self.config.align_image_distances >= 1 and optimize_num % self.config.align_image_distances == 0 and optimize_num > 0:
+                print("Aligning images...")
+                tmp_new_geometry = align_geometry_by_image_distance(np.array(new_geometry))
+                for k in range(len(new_geometry)):
+                    new_geometry[k] = copy.copy(tmp_new_geometry[k])
+
             # Save analysis files
             tmp_instance_fileio = FileIO(file_directory + "/", "dummy.txt")
             tmp_instance_fileio.argrelextrema_txt_save(biased_energy_list, "approx_TS_node", "max")
@@ -896,4 +902,145 @@ def distribute_geometry_by_length(geometry_list, angstrom_spacing):
         new_geometry_list.append(new_geometry)
 
     new_geometry_list.append(geometry_list[-1])
-    return new_geometry_list
+    return new_geometry_list
+
+
+def align_geometry_by_image_distance(geometry_list, spacing_dist=None, max_iterations=1000, tolerance=1e-6):
+    """Distribute geometries by specified distance spacing"""
+
+    for i in range(len(geometry_list)-1):
+        geom_i = geometry_list[i]
+        geom_j = geometry_list[i+1]
+        geom_i_after_kabsch, geom_j_after_kabsch = Calculationtools().kabsch_algorithm(geom_i, geom_j)
+        geometry_list[i] = geom_i_after_kabsch
+    geometry_list[-1] = geom_j_after_kabsch
+
+    # 1. Create path length list
+    path_length_list = calc_path_length_list(geometry_list)
+    total_length = path_length_list[-1]
+    nnode = len(geometry_list)
+    #print("Number of nodes : ", nnode)
+    #print("Path length list (before the process) : ", path_length_list)
+    #print("Total path length : ", total_length)
+
+    # 2. Calculate equal spacing distance
+    if spacing_dist is not None:
+        node_distance = spacing_dist
+        n_new_nodes = int(total_length / node_distance) + 1
+    else:
+        n_new_nodes = nnode
+        node_distance = total_length / (n_new_nodes - 1)
+    #print("Node distance : ", node_distance)
+
+    # 3. Create target distance list with equal spacing
+    target_distances = [i * node_distance for i in range(n_new_nodes)]
+    
+    # 4. Generate nodes with accurate equal spacing using iterative method
+    new_geometry_list = []
+    
+    # Fix the first node
+    new_geometry_list.append(geometry_list[0])
+    
+    # Generate remaining nodes (excluding first and last nodes)
+    for i in range(1, n_new_nodes-1):
+        # Target distance
+        target_distance = node_distance
+        
+        # Initial estimated position
+        # Find which segment it belongs to
+        for j in range(len(path_length_list) - 1):
+            if path_length_list[j] <= target_distances[i] <= path_length_list[j+1]:
+                # Linear interpolation
+                t = (target_distances[i] - path_length_list[j]) / (path_length_list[j+1] - path_length_list[j])
+                current_geom = geometry_list[j] + t * (geometry_list[j+1] - geometry_list[j])
+                segment_start = j
+                break
+        
+        # Iteratively adjust to achieve exact distance from previous node
+        prev_geom = new_geometry_list[-1]
+        iteration = 0
+        
+        # Track current position parameters along the path
+        current_segment = segment_start
+        current_t = t
+        
+        while iteration < max_iterations:
+            # Calculate distance to previous node
+            prev_geom_centered = prev_geom - np.mean(prev_geom, axis=0)
+            current_geom_centered = current_geom - np.mean(current_geom, axis=0)
+            current_distance = np.linalg.norm(current_geom_centered - prev_geom_centered)
+            
+            # Exit if distance is accurate enough
+            if abs(current_distance - target_distance) < tolerance:
+                break
+            
+            # Move forward along the path if distance is too small, backward if too large
+            if current_distance < target_distance:
+                # Move forward
+                delta_dist = target_distance - current_distance
+                
+                # Advance within current segment
+                segment_length = path_length_list[current_segment+1] - path_length_list[current_segment]
+                remaining_segment = (1 - current_t) * segment_length
+                
+                if delta_dist <= remaining_segment:
+                    # Adjustment possible within current segment
+                    move_t = delta_dist / segment_length
+                    current_t += move_t
+                    current_geom = geometry_list[current_segment] + current_t * (geometry_list[current_segment+1] - geometry_list[current_segment])
+                else:
+                    # Need to move to next segment
+                    remaining_dist = delta_dist - remaining_segment
+                    
+                    # If next segment exists
+                    if current_segment + 1 < len(geometry_list) - 1:
+                        current_segment += 1
+                        next_segment_length = path_length_list[current_segment+1] - path_length_list[current_segment]
+                        current_t = min(remaining_dist / next_segment_length, 0.99)  # Prevent boundary crossing
+                        current_geom = geometry_list[current_segment] + current_t * (geometry_list[current_segment+1] - geometry_list[current_segment])
+                    else:
+                        # If we can't move further, stay at end of last segment
+                        current_t = 0.99
+                        current_geom = geometry_list[current_segment] + current_t * (geometry_list[current_segment+1] - geometry_list[current_segment])
+            else:
+                # Move backward
+                delta_dist = current_distance - target_distance
+                
+                # Retreat within current segment
+                segment_length = path_length_list[current_segment+1] - path_length_list[current_segment]
+                traversed_segment = current_t * segment_length
+                
+                if delta_dist <= traversed_segment:
+                    # Adjustment possible within current segment
+                    move_t = delta_dist / segment_length
+                    current_t -= move_t
+                    current_geom = geometry_list[current_segment] + current_t * (geometry_list[current_segment+1] - geometry_list[current_segment])
+                else:
+                    # Need to move to previous segment
+                    remaining_dist = delta_dist - traversed_segment
+                    
+                    # If previous segment exists
+                    if current_segment > 0:
+                        current_segment -= 1
+                        prev_segment_length = path_length_list[current_segment+1] - path_length_list[current_segment]
+                        current_t = max(1 - (remaining_dist / prev_segment_length), 0.01)  # Prevent boundary crossing
+                        current_geom = geometry_list[current_segment] + current_t * (geometry_list[current_segment+1] - geometry_list[current_segment])
+                    else:
+                        # If we can't move further, stay at start of first segment
+                        current_t = 0.01
+                        current_geom = geometry_list[current_segment] + current_t * (geometry_list[current_segment+1] - geometry_list[current_segment])
+            
+            iteration += 1
+        
+        new_geometry_list.append(current_geom)
+    
+    # Use last node from input
+    new_geometry_list.append(geometry_list[-1])
+
+    new_path_length_list = calc_path_length_list(new_geometry_list)
+    #print("Path length list (after the process) : ", new_path_length_list)
+    #print("Distances between nodes:")
+    #for x in range(len(new_path_length_list)-1):
+    #    print(new_path_length_list[x+1]-new_path_length_list[x])
+
+    return new_geometry_list