@@ -132,15 +132,23 @@ def load_force_constant_matrix(self):
132132 only the last is returned. The units of the returned force constants
133133 are J/m^2. If no force constant matrix can be found in the log file,
134134 ``None`` is returned.
135+ Also checks that the force constant matrix was computed using the correct
136+ (input orientation Cartesian) coordinates.
137+ IOP(2/9=2000) must be specified for large cases (14+ atoms)
135138 """
136139 force = None
137140
141+ iop2_9_equals_2000 = False
142+
138143 n_atoms = self .get_number_of_atoms ()
139144 n_rows = n_atoms * 3
140145
141146 with open (self .path , 'r' ) as f :
142147 line = f .readline ()
143148 while line != '' :
149+ if '2/9=2000' in line :
150+ iop2_9_equals_2000 = True
151+
144152 # Read force constant matrix
145153 if 'Force constants in Cartesian coordinates:' in line :
146154 force = np .zeros ((n_rows , n_rows ), float )
@@ -157,25 +165,24 @@ def load_force_constant_matrix(self):
157165 force *= 4.35974417e-18 / 5.291772108e-11 ** 2
158166 line = f .readline ()
159167
168+ if n_atoms > 13 and not iop2_9_equals_2000 :
169+ raise LogError (f'Gaussian output file { self .path } contains more than 13 atoms. '
170+ f'Please add the `iop(2/9=2000)` keyword to your input file '
171+ f'so Gaussian will compute force matrix using the input orientation Cartesians.' )
172+
160173 return force
161174
162175 def load_geometry (self ):
163176 """
164177 Return the optimum geometry of the molecular configuration from the
165178 Gaussian log file. If multiple such geometries are identified, only the
166179 last is returned.
167- Also checks that the Cartesian coordinates are printed in the input orientation:
168- IOP(2/9=2000) must be specified for large cases (14+ atoms)
169180 """
170181 number , coord , mass = [], [], []
171-
172- iop2_9_equals_2000 = False
182+
173183 with open (self .path , 'r' ) as f :
174184 line = f .readline ()
175185 while line != '' :
176- if '2/9=2000' in line :
177- iop2_9_equals_2000 = True
178-
179186 # Automatically determine the number of atoms
180187 if 'Input orientation:' in line :
181188 number , coord = [], []
@@ -202,11 +209,6 @@ def load_geometry(self):
202209 '50 or more atoms, you will need to add the `iop(2/9=2000)` keyword to your '
203210 'input file so Gaussian will print the input orientation geometry.' .format (self .path ))
204211
205- if len (number ) > 13 and not iop2_9_equals_2000 :
206- raise LogError (f'Gaussian output file { self .path } contains more than 13 atoms. '
207- f'Please add the `iop(2/9=2000)` keyword to your input file '
208- f'so Gaussian will print the geometry using the input orientation.' )
209-
210212 return coord , number , mass
211213
212214 def load_conformer (self , symmetry = None , spin_multiplicity = 0 , optical_isomers = None , label = '' ):
0 commit comments