From 9028fbfcda08e7dc378ca3b505a89e5a12e58a1c Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Fri, 7 Feb 2025 11:27:02 -0500 Subject: [PATCH 1/2] add check for iop(2/9=2000) in gaussian load geometry --- arkane/ess/gaussian.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/arkane/ess/gaussian.py b/arkane/ess/gaussian.py index 156823b987..e9c55a1c9a 100644 --- a/arkane/ess/gaussian.py +++ b/arkane/ess/gaussian.py @@ -164,12 +164,18 @@ def load_geometry(self): Return the optimum geometry of the molecular configuration from the Gaussian log file. If multiple such geometries are identified, only the last is returned. + Also checks that the Cartesian coordinates are printed in the input orientation: + IOP(2/9=2000) must be specified for large cases (14+ atoms) """ number, coord, mass = [], [], [] + iop2_9_equals_2000 = False with open(self.path, 'r') as f: line = f.readline() while line != '': + if '2/9=2000' in line: + iop2_9_equals_2000 = True + # Automatically determine the number of atoms if 'Input orientation:' in line: number, coord = [], [] @@ -196,6 +202,11 @@ def load_geometry(self): '50 or more atoms, you will need to add the `iop(2/9=2000)` keyword to your ' 'input file so Gaussian will print the input orientation geometry.'.format(self.path)) + if len(number) > 13 and not iop2_9_equals_2000: + raise LogError(f'Gaussian output file {self.path} contains more than 13 atoms. ' + f'Please add the `iop(2/9=2000)` keyword to your input file ' + f'so Gaussian will print the geometry using the input orientation.') + return coord, number, mass def load_conformer(self, symmetry=None, spin_multiplicity=0, optical_isomers=None, label=''): From 1a26ad32b4fd48b903edfe5395b7af337b209906 Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Thu, 13 Feb 2025 15:16:35 -0500 Subject: [PATCH 2/2] move iop 2/9=2000 check from geometry to force constant matrix section --- arkane/ess/gaussian.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/arkane/ess/gaussian.py b/arkane/ess/gaussian.py index e9c55a1c9a..7e0e144840 100644 --- a/arkane/ess/gaussian.py +++ b/arkane/ess/gaussian.py @@ -132,15 +132,23 @@ def load_force_constant_matrix(self): only the last is returned. The units of the returned force constants are J/m^2. If no force constant matrix can be found in the log file, ``None`` is returned. + Also checks that the force constant matrix was computed using the correct + (input orientation Cartesian) coordinates. + IOP(2/9=2000) must be specified for large cases (14+ atoms) """ force = None + iop2_9_equals_2000 = False + n_atoms = self.get_number_of_atoms() n_rows = n_atoms * 3 with open(self.path, 'r') as f: line = f.readline() while line != '': + if '2/9=2000' in line: + iop2_9_equals_2000 = True + # Read force constant matrix if 'Force constants in Cartesian coordinates:' in line: force = np.zeros((n_rows, n_rows), float) @@ -157,6 +165,11 @@ def load_force_constant_matrix(self): force *= 4.35974417e-18 / 5.291772108e-11 ** 2 line = f.readline() + if n_atoms > 13 and not iop2_9_equals_2000: + raise LogError(f'Gaussian output file {self.path} contains more than 13 atoms. ' + f'Please add the `iop(2/9=2000)` keyword to your input file ' + f'so Gaussian will compute force matrix using the input orientation Cartesians.') + return force def load_geometry(self): @@ -164,18 +177,12 @@ def load_geometry(self): Return the optimum geometry of the molecular configuration from the Gaussian log file. If multiple such geometries are identified, only the last is returned. - Also checks that the Cartesian coordinates are printed in the input orientation: - IOP(2/9=2000) must be specified for large cases (14+ atoms) """ number, coord, mass = [], [], [] - - iop2_9_equals_2000 = False + with open(self.path, 'r') as f: line = f.readline() while line != '': - if '2/9=2000' in line: - iop2_9_equals_2000 = True - # Automatically determine the number of atoms if 'Input orientation:' in line: number, coord = [], [] @@ -202,11 +209,6 @@ def load_geometry(self): '50 or more atoms, you will need to add the `iop(2/9=2000)` keyword to your ' 'input file so Gaussian will print the input orientation geometry.'.format(self.path)) - if len(number) > 13 and not iop2_9_equals_2000: - raise LogError(f'Gaussian output file {self.path} contains more than 13 atoms. ' - f'Please add the `iop(2/9=2000)` keyword to your input file ' - f'so Gaussian will print the geometry using the input orientation.') - return coord, number, mass def load_conformer(self, symmetry=None, spin_multiplicity=0, optical_isomers=None, label=''):