@@ -1198,24 +1198,20 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_
11981198 d_int [j , i ] += d_int_proj [k , i ] * vmw [j , k ]
11991199
12001200 # Ortho normalize
1201- for i in range (n_rotors ):
1202- norm = 0.0
1203- for j in range (3 * n_atoms ):
1204- norm += d_int [j , i ] * d_int [j , i ]
1205- for j in range (3 * n_atoms ):
1206- d_int [j , i ] /= np .sqrt (norm )
1207- for j in range (i + 1 , n_rotors ):
1208- proj = 0.0
1209- for k in range (3 * n_atoms ):
1210- proj += d_int [k , i ] * d_int [k , j ]
1211- for k in range (3 * n_atoms ):
1212- d_int [k , j ] -= proj * d_int [k , i ]
1201+ # Here, Q will have orthonormal columns, and R is the triangular factor.
1202+ Q , R = np .linalg .qr (d_int , mode = 'reduced' )
1203+ # replace d_int with its orthonormalized version:
1204+ d_int = Q
12131205
12141206 # calculate the frequencies corresponding to the internal rotors
12151207 int_proj = np .dot (fm , d_int )
12161208 kmus = np .array ([np .linalg .norm (int_proj [:, i ]) for i in range (int_proj .shape [1 ])])
12171209 int_rotor_freqs = np .sqrt (kmus ) / (2.0 * math .pi * constants .c * 100.0 )
12181210
1211+ logging .debug ('Frequencies from internal rotors:' )
1212+ for i in range (n_rotors ):
1213+ logging .debug (' rotor %d: %.6f cm^-1' , i , int_rotor_freqs [i ])
1214+
12191215 if get_projected_out_freqs :
12201216 return int_rotor_freqs
12211217
0 commit comments