From cba999effdd368cd922050819750dfab4cf0437f 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 1a134b7d1675c3df671fd5bca945a9d82961c73d 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 51e359a6fb67b60f1e77ab3215b6c9f1870bb16e 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 786ebe635f9265920d2e0a1dec75fad8f1f8adb6 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 03c26251e5b4f232399fe8026abf44160e17d9dc 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