add check for iop(2/9=2000) in gaussian load geometry

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=''):