diff --git a/arkane/ess/gaussian.py b/arkane/ess/gaussian.py
index 156823b987..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):
@@ -166,7 +179,7 @@ def load_geometry(self):
         last is returned.
         """
         number, coord, mass = [], [], []
-
+        
         with open(self.path, 'r') as f:
             line = f.readline()
             while line != '':