Add files via upload

Fixed an issue in Fischer's model Hessian generation where numerical instability occurred when four atoms were aligned linearly.

Author: ss0832

multioptpy/ModelHessian/fischer.py

@@ -131,51 +131,84 @@ def fischer_angle(self, coord, element_list):
                     self.cart_hess[3*k+n, 3*j+m] += force_const * b_vec[2][n] * b_vec[1][m]
 
     def fischer_dihedral(self, coord, element_list, bond_mat):
-        """Calculate Hessian components for dihedral torsions"""
+        """Calculate Hessian components for dihedral torsions with linearity check"""
         BC = BondConnectivity()
         b_c_mat = BC.bond_connect_matrix(element_list, coord)
         dihedral_indices = BC.dihedral_angle_connect_table(b_c_mat)
 
+        # Helper to calculate sin^2 of bond angle
+        def get_sin_sq_angle(idx1, idx2, idx3):
+            v1 = coord[idx1] - coord[idx2]
+            v2 = coord[idx3] - coord[idx2]
+            # Cross product magnitude squared
+            cross_prod = np.cross(v1, v2)
+            cross_sq = np.dot(cross_prod, cross_prod)
+            # Squared norms
+            n1_sq = np.dot(v1, v1)
+            n2_sq = np.dot(v2, v2)
+            if n1_sq * n2_sq < 1e-12:
+                return 0.0
+            return cross_sq / (n1_sq * n2_sq)
+
         for idx in dihedral_indices:
-            i, j, k, l = idx  # i-j-k-l dihedral
+            # i-j-k-l dihedral
+            i, j, k, l = idx 
+            
+            # --- Linearity Check ---
+            # Check bond angles I-J-K and J-K-L.
+            # If atoms are linear, the torsion definition is singular (B-matrix blows up).
+            # We use sin^2(theta) because it avoids acos and is efficient.
+            sin_sq_ijk = get_sin_sq_angle(i, j, k)
+            sin_sq_jkl = get_sin_sq_angle(j, k, l)
 
-            # Central bond in dihedral
+            # Threshold: if sin(theta) < 0.17 (approx 10 degrees from linear), dampen or skip.
+            # Using sin^2 < 0.03 as cutoff.
+            # The Wilson B-matrix scales as 1/sin(theta). 
+            # If we are too close to linear, we must skip to avoid numerical explosion.
+            if sin_sq_ijk < 1.0e-3 or sin_sq_jkl < 1.0e-3:
+                continue
+            # -----------------------
+
+            # Central bond in dihedral (j-k)
             r_jk = np.linalg.norm(coord[j] - coord[k])
             r_jk_cov = covalent_radii_lib(element_list[j]) + covalent_radii_lib(element_list[k])
 
             # Count bonds to central atoms
             bond_sum = self.count_bonds_for_dihedral(bond_mat, (j, k))
 
-            # Calculate force constant using Fischer's formula
+            # Calculate force constant
             force_const = self.calc_dihedral_force_const(r_jk, r_jk_cov, bond_sum)
 
+            # Optional: Additional damping can be applied here if desired, 
+            # but the skip above is the most robust fix for the singularity.
+            
             # Convert to Cartesian coordinates
             t_xyz = np.array([coord[i], coord[j], coord[k], coord[l]])
-            tau, b_vec = torsion2(t_xyz)
 
-            for n in range(3):
-                for m in range(3):
-                    self.cart_hess[3*i+n, 3*i+m] += force_const * b_vec[0][n] * b_vec[0][m]
-                    self.cart_hess[3*j+n, 3*j+m] += force_const * b_vec[1][n] * b_vec[1][m]
-                    self.cart_hess[3*k+n, 3*k+m] += force_const * b_vec[2][n] * b_vec[2][m]
-                    self.cart_hess[3*l+n, 3*l+m] += force_const * b_vec[3][n] * b_vec[3][m]
-                    
-                    self.cart_hess[3*i+n, 3*j+m] += force_const * b_vec[0][n] * b_vec[1][m]
-                    self.cart_hess[3*i+n, 3*k+m] += force_const * b_vec[0][n] * b_vec[2][m]
-                    self.cart_hess[3*i+n, 3*l+m] += force_const * b_vec[0][n] * b_vec[3][m]
-                    
-                    self.cart_hess[3*j+n, 3*i+m] += force_const * b_vec[1][n] * b_vec[0][m]
-                    self.cart_hess[3*j+n, 3*k+m] += force_const * b_vec[1][n] * b_vec[2][m]
-                    self.cart_hess[3*j+n, 3*l+m] += force_const * b_vec[1][n] * b_vec[3][m]
+            # Note: torsion2 might still be unstable if not skipped
+            try:
+                tau, b_vec = torsion2(t_xyz)
+            except Exception:
+                # Fallback if torsion2 fails numerically
+                continue
+            
+            # Iterate over all pairs of atoms in the dihedral
+            atom_indices = [i, j, k, l]
+            
+            for a in range(4):
+                for b in range(4):
+                    atom_a = atom_indices[a]
+                    atom_b = atom_indices[b]
 
-                    self.cart_hess[3*k+n, 3*i+m] += force_const * b_vec[2][n] * b_vec[0][m]
-                    self.cart_hess[3*k+n, 3*j+m] += force_const * b_vec[2][n] * b_vec[1][m]
-                    self.cart_hess[3*k+n, 3*l+m] += force_const * b_vec[2][n] * b_vec[3][m]
+                    vec_a = b_vec[a]
+                    vec_b = b_vec[b]
 
-                    self.cart_hess[3*l+n, 3*i+m] += force_const * b_vec[3][n] * b_vec[0][m]
-                    self.cart_hess[3*l+n, 3*j+m] += force_const * b_vec[3][n] * b_vec[1][m]
-                    self.cart_hess[3*l+n, 3*k+m] += force_const * b_vec[3][n] * b_vec[2][m]
-    
+                    # Update the 3x3 Hessian block
+                    for n in range(3):
+                        for m in range(3):
+                            val = force_const * vec_a[n] * vec_b[m]
+                            self.cart_hess[3*atom_a+n, 3*atom_b+m] += val
+                            
     def main(self, coord, element_list, cart_gradient):
         """Main function to generate Fischer's approximate Hessian"""
         print("Generating Fischer's approximate hessian...")

multioptpy/ModelHessian/fischerd3.py

@@ -211,8 +211,9 @@ def fischer_angle(self, coord, element_list):
                     H_mn_block = force_const * np.outer(b_vec[m_idx], b_vec[n_idx])
                     self.cart_hess[start_m:end_m, start_n:end_n] += H_mn_block
 
+
     def fischer_dihedral(self, coord, element_list, bond_mat):
-        """Calculate Hessian components for dihedral torsion (Optimized with slicing)"""
+        """Calculate Hessian components for dihedral torsion (Optimized with singularity damping)"""
         BC = BondConnectivity()
         b_c_mat = BC.bond_connect_matrix(element_list, coord)
         dihedral_indices = BC.dihedral_angle_connect_table(b_c_mat)
@@ -221,22 +222,71 @@ def fischer_dihedral(self, coord, element_list, bond_mat):
         tors_atom_bonds = {}
         for idx in dihedral_indices:
             i, j, k, l = idx  # i-j-k-l dihedral
-            # bond_sum is count_bonds_for_dihedral(bond_mat, (j, k))
             bond_sum = bond_mat[j].sum() + bond_mat[k].sum() - 2
             tors_atom_bonds[(j, k)] = bond_sum
 
         for idx in dihedral_indices:
             i, j, k, l = idx
 
-            r_jk = np.linalg.norm(coord[j] - coord[k])
+            # Vector calculations
+            vec_ji = coord[i] - coord[j]
+            vec_jk = coord[k] - coord[j]
+            vec_kl = coord[l] - coord[k]
+            
+            r_jk = np.linalg.norm(vec_jk)
             r_jk_cov = covalent_radii_lib(element_list[j]) + covalent_radii_lib(element_list[k])
 
             bond_sum = tors_atom_bonds.get((j, k), 0)
 
+            # Calculate base force constant
             force_const = self.calc_dihedral_force_const(r_jk, r_jk_cov, bond_sum)
 
+            # --- Singularity Handling (Damping for linear angles) ---
+            # Determine linearity of angles i-j-k and j-k-l
+            
+            # Angle 1: i-j-k 
+            n_ji = np.linalg.norm(vec_ji)
+            n_jk = r_jk # already calculated
+            
+            # Avoid division by zero if atoms overlap
+            if n_ji < 1e-8 or n_jk < 1e-8: continue
+            
+            cos_theta1 = np.dot(vec_ji, vec_jk) / (n_ji * n_jk)
+            
+            # Angle 2: j-k-l (Note: vec_kj is -vec_jk)
+            vec_kj = -vec_jk
+            n_kl = np.linalg.norm(vec_kl)
+            
+            if n_kl < 1e-8: continue
+            
+            cos_theta2 = np.dot(vec_kj, vec_kl) / (n_jk * n_kl)
+
+            # --- Damping Factor Calculation ---
+            # The Wilson B-matrix contains 1/sin(theta) terms which diverge at 180 deg.
+            # We scale the force constant by sin^2(theta) to cancel this divergence.
+            # sin^2(theta) = 1 - cos^2(theta)
+            
+            sin2_theta1 = 1.0 - min(cos_theta1**2, 1.0)
+            sin2_theta2 = 1.0 - min(cos_theta2**2, 1.0)
+            
+            # Hard cutoff: If geometry is extremely linear, skip to avoid NaN
+            if sin2_theta1 < 1e-4 or sin2_theta2 < 1e-4:
+                continue
+
+            # Apply scaling factor
+            # This ensures the force constant goes to 0 as the angle becomes linear
+            scaling_factor = sin2_theta1 * sin2_theta2
+            force_const *= scaling_factor
+            
+            # --------------------------------------------------------
+
             t_xyz = np.array([coord[i], coord[j], coord[k], coord[l]])
-            tau, b_vec = torsion2(t_xyz)
+            
+            try:
+                tau, b_vec = torsion2(t_xyz)
+            except (ValueError, ArithmeticError):
+                # Skip if numerical errors occur in torsion calculation
+                continue
 
             # Optimized: Use NumPy slicing and outer product
             atoms = [i, j, k, l]