From 47de8880a5a81f370e0a4ac47d4df1a41c90c149 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 6 Feb 2025 23:04:26 -0500 Subject: [PATCH 1/5] Rename 'inertia' to 'identity' This was originally called 'I' and during the PEP8 conversion around the py2 to py3 time, they were all renamed inertia, even though some of them represent the identity matrix not the moments of inertia --- arkane/statmech.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/arkane/statmech.py b/arkane/statmech.py index af82a045bb8..802f7efe5be 100644 --- a/arkane/statmech.py +++ b/arkane/statmech.py @@ -1078,13 +1078,10 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_ d[3 * i + 2, 5] = (p[i, 0] * inertia_xyz[2, 1] - p[i, 1] * inertia_xyz[2, 0]) * amass[i] # Make sure projection matrix is orthonormal - - inertia = np.identity(n_atoms * 3, float) - + identity = np.identity(n_atoms * 3, float) p = np.zeros((n_atoms * 3, 3 * n_atoms + external), float) - p[:, 0:external] = d[:, 0:external] - p[:, external:external + 3 * n_atoms] = inertia[:, 0:3 * n_atoms] + p[:, external:external + 3 * n_atoms] = identity[:, 0:3 * n_atoms] for i in range(3 * n_atoms + external): norm = 0.0 @@ -1234,8 +1231,8 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_ # Do the projection d_int_proj = np.dot(vmw.T, d_int) proj = np.dot(d_int, d_int.T) - inertia = np.identity(n_atoms * 3, float) - proj = inertia - proj + identity = np.identity(n_atoms * 3, float) + proj = identity - proj fm = np.dot(proj, np.dot(fm, proj)) # Get eigenvalues of mass-weighted force constant matrix eig, v = np.linalg.eigh(fm) From 014019d67d83a7fb2267b547312f59a76b12c0be Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 6 Feb 2025 23:05:31 -0500 Subject: [PATCH 2/5] Refactor some column sorting. I think the previous might have left the final column in an undetermined state. And this might be quicker (fewer loops) and cleaner Python. --- arkane/statmech.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/arkane/statmech.py b/arkane/statmech.py index 802f7efe5be..8ad7e77e678 100644 --- a/arkane/statmech.py +++ b/arkane/statmech.py @@ -1101,14 +1101,10 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_ # Order p, there will be vectors that are 0.0 i = 0 - while i < 3 * n_atoms: - norm = 0.0 - for j in range(3 * n_atoms): - norm += p[j, i] * p[j, i] - if norm < 0.5: - p[:, i:3 * n_atoms + external - 1] = p[:, i + 1:3 * n_atoms + external] - else: - i += 1 + is_zero_column = (p*p).sum(axis=0) < 0.5 + # put the columns that are zeros at the end + temporary = np.hstack((p[:, ~is_zero_column], p[:, is_zero_column])) + p[:, :] = temporary # T is the transformation vector from cartesian to internal coordinates T = np.zeros((n_atoms * 3, 3 * n_atoms - external), float) From 515faa44dc7d102158341423432662339f577331 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 6 Feb 2025 23:23:33 -0500 Subject: [PATCH 3/5] Added /edited comments --- arkane/statmech.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/arkane/statmech.py b/arkane/statmech.py index 8ad7e77e678..184c6cab8a1 100644 --- a/arkane/statmech.py +++ b/arkane/statmech.py @@ -1083,6 +1083,7 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_ p[:, 0:external] = d[:, 0:external] p[:, external:external + 3 * n_atoms] = identity[:, 0:3 * n_atoms] + # modified Gram–Schmidt orthonormalization for i in range(3 * n_atoms + external): norm = 0.0 for j in range(3 * n_atoms): @@ -1091,7 +1092,7 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_ if norm > 1E-15: p[j, i] /= np.sqrt(norm) else: - p[j, i] = 0.0 + p[j, i] = 0.0 # zeroing out vectors that are nearly zero or dependent, could lose a basis for j in range(i + 1, 3 * n_atoms + external): proj = 0.0 for k in range(3 * n_atoms): @@ -1099,7 +1100,7 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_ for k in range(3 * n_atoms): p[k, j] -= proj * p[k, i] - # Order p, there will be vectors that are 0.0 + # Order p, since there will be vectors that are 0.0 i = 0 is_zero_column = (p*p).sum(axis=0) < 0.5 # put the columns that are zeros at the end From 075dcebf30af25219cf009673757de1f18c82dad Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 7 Feb 2025 00:09:14 -0500 Subject: [PATCH 4/5] Change orthonormalization when projecting out modes in statmech. We were doing a manually coded Gram-Schmidt orthonormalization. This is now replaced with a QR decomposition built in to Numpy, which should be more robust (and probably faster). Added in some checks that print warnings. --- arkane/statmech.py | 48 ++++++++++++++++++++-------------------------- 1 file changed, 21 insertions(+), 27 deletions(-) diff --git a/arkane/statmech.py b/arkane/statmech.py index 184c6cab8a1..d874b976509 100644 --- a/arkane/statmech.py +++ b/arkane/statmech.py @@ -1078,39 +1078,33 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_ d[3 * i + 2, 5] = (p[i, 0] * inertia_xyz[2, 1] - p[i, 1] * inertia_xyz[2, 0]) * amass[i] # Make sure projection matrix is orthonormal - identity = np.identity(n_atoms * 3, float) - p = np.zeros((n_atoms * 3, 3 * n_atoms + external), float) - p[:, 0:external] = d[:, 0:external] - p[:, external:external + 3 * n_atoms] = identity[:, 0:3 * n_atoms] + p = np.hstack((d, np.identity(n_atoms * 3, float))) - # modified Gram–Schmidt orthonormalization - for i in range(3 * n_atoms + external): - norm = 0.0 - for j in range(3 * n_atoms): - norm += p[j, i] * p[j, i] - for j in range(3 * n_atoms): - if norm > 1E-15: - p[j, i] /= np.sqrt(norm) - else: - p[j, i] = 0.0 # zeroing out vectors that are nearly zero or dependent, could lose a basis - for j in range(i + 1, 3 * n_atoms + external): - proj = 0.0 - for k in range(3 * n_atoms): - proj += p[k, i] * p[k, j] - for k in range(3 * n_atoms): - p[k, j] -= proj * p[k, i] + # Compute the QR decomposition in reduced mode + Q, R = np.linalg.qr(p, mode='reduced') + + # Verify that Q has orthonormal columns: + if np.isclose(np.dot(Q.T, Q), np.eye(Q.shape[1])).all(): + logging.debug("Q has orthonormal columns") + else: + logging.warning("Q does not have orthonormal columns") + # Verify the reconstruction + reconstructed_p = Q @ R + if np.isclose(p, reconstructed_p).all(): + logging.debug("Reconstruction of projection matrix is correct") + else: + logging.warning("Reconstruction of projection matrix is incorrect") + + # Check for nearly zero columns in the QR decomposition + for i, rkk in enumerate(R.diagonal()): + if abs(rkk) < 1E-15: + logging.warning(f'Column {i} of the QR decomposition is nearly zero, could lose a basis') - # Order p, since there will be vectors that are 0.0 - i = 0 - is_zero_column = (p*p).sum(axis=0) < 0.5 - # put the columns that are zeros at the end - temporary = np.hstack((p[:, ~is_zero_column], p[:, is_zero_column])) - p[:, :] = temporary # T is the transformation vector from cartesian to internal coordinates T = np.zeros((n_atoms * 3, 3 * n_atoms - external), float) - T[:, 0:3 * n_atoms - external] = p[:, external:3 * n_atoms] + T[:, :] = Q[:, external:3 * n_atoms] # Generate mass-weighted force constant matrix # This converts the axes to mass-weighted Cartesian axes From 3fded4d388c5042b1c3a332315b501cb856691cc Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 7 Feb 2025 00:40:35 -0500 Subject: [PATCH 5/5] Replace another gram-schmidt with QR decomposition Not sure it helps at all. --- arkane/statmech.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/arkane/statmech.py b/arkane/statmech.py index d874b976509..19d8a5c3e0e 100644 --- a/arkane/statmech.py +++ b/arkane/statmech.py @@ -1198,24 +1198,20 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_ d_int[j, i] += d_int_proj[k, i] * vmw[j, k] # Ortho normalize - for i in range(n_rotors): - norm = 0.0 - for j in range(3 * n_atoms): - norm += d_int[j, i] * d_int[j, i] - for j in range(3 * n_atoms): - d_int[j, i] /= np.sqrt(norm) - for j in range(i + 1, n_rotors): - proj = 0.0 - for k in range(3 * n_atoms): - proj += d_int[k, i] * d_int[k, j] - for k in range(3 * n_atoms): - d_int[k, j] -= proj * d_int[k, i] + # Here, Q will have orthonormal columns, and R is the triangular factor. + Q, R = np.linalg.qr(d_int, mode='reduced') + # replace d_int with its orthonormalized version: + d_int = Q # calculate the frequencies corresponding to the internal rotors int_proj = np.dot(fm, d_int) kmus = np.array([np.linalg.norm(int_proj[:, i]) for i in range(int_proj.shape[1])]) int_rotor_freqs = np.sqrt(kmus) / (2.0 * math.pi * constants.c * 100.0) + logging.debug('Frequencies from internal rotors:') + for i in range(n_rotors): + logging.debug(' rotor %d: %.6f cm^-1', i, int_rotor_freqs[i]) + if get_projected_out_freqs: return int_rotor_freqs