Merge pull request #10 from duartegroup/tom

Switch gradients to be consistent with unit convention

.gitignore

@@ -46,3 +46,9 @@ tests/.autode_calculations
 tests/data/complex_conf*.xyz
 
 tests/data/*.com
+
+tests/data/conformers/.autode_calculations
+
+tests/data/orca/.autode_calculations
+
+tests/data/xtb/.autode_calculations

autode/calculation.py

@@ -276,7 +276,7 @@ def get_gradients(self):
         calculation
 
         Returns:
-            (np.ndarray): Gradient vectors for each atom (Å)
+            (np.ndarray): Gradient vectors for each atom (Ha Å^-1)
                           gradients.shape = (n_atoms, 3)
         """
         logger.info(f'Getting gradients from {self.output.filename}')

autode/constants.py

@@ -1,7 +1,8 @@
-class Constants(object):
+class Constants:
 
-    ha2kcalmol = 627.509  # Hartee^-1 kcal mol^-1
-    ha2kJmol = 2625.50    # Hartree^-1 kJ mol^-1
-    eV2ha = 0.0367493     # Hartree ev^-1
-    a02ang = 0.529177     # Å bohr^-1
-    kcal2kJ = 4.184       # kJ kcal^-1
+    ha2kcalmol = 627.509   # Hartee^-1 kcal mol^-1
+    ha2kJmol = 2625.50     # Hartree^-1 kJ mol^-1
+    eV2ha = 0.0367493      # Hartree ev^-1
+    a02ang = 0.529177      # Å bohr^-1
+    ang2a0 = 1.0 / a02ang  # bohr Å^-1
+    kcal2kJ = 4.184        # kJ kcal^-1

autode/smiles/smiles.py

@@ -12,7 +12,9 @@
 
 
 def calc_multiplicity(molecule, n_radical_electrons):
-    """Calculate the spin multiplicity 2S + 1 where S is the number of unpaired electrons
+    """Calculate the spin multiplicity 2S + 1 where S is the number of
+    unpaired electrons
+
     Arguments:
         molecule (autode.molecule.Molecule):
         n_radical_electrons (int):
@@ -24,20 +26,23 @@ def calc_multiplicity(molecule, n_radical_electrons):
         return 1
 
     if molecule.mult == 1 and n_radical_electrons == 1:
-        # Cannot have multiplicity = 1 and 1 radical electrons – override default multiplicity
+        # Cannot have multiplicity = 1 and 1 radical electrons – override
+        # default multiplicity
         return 2
 
     if molecule.mult == 1 and n_radical_electrons > 1:
-        logger.warning('Diradicals by default singlets. Set mol.mult if it\'s any different')
+        logger.warning('Diradicals by default singlets. Set mol.mult if it\'s '
+                       'any different')
         return 1
 
     return molecule.mult
 
 
 def init_organic_smiles(molecule, smiles):
     """
-    Initialise a molecule from a SMILES string, set the charge, multiplicity (if it's not already specified) and the 3D
-    geometry using RDKit
+    Initialise a molecule from a SMILES string, set the charge, multiplicity (
+    if it's not already specified) and the 3D geometry using RDKit
+
     Arguments:
         molecule (autode.molecule.Molecule):
         smiles (str): SMILES string
@@ -87,7 +92,6 @@ def init_smiles(molecule, smiles):
         molecule (autode.molecule.Molecule):
         smiles (str): SMILES string
     """
-
     # Assume that the RDKit conformer generation algorithm is not okay for metals
     molecule.rdkit_conf_gen_is_fine = False
 
@@ -115,7 +119,8 @@ def init_smiles(molecule, smiles):
 
 def check_bonds(molecule, bonds):
     """
-    Ensure the SMILES string and the 3D structure have the same bonds, but don't override
+    Ensure the SMILES string and the 3D structure have the same bonds,
+    but don't override
 
     Arguments:
         molecule (autode.molecule.Molecule):

autode/transition_states/transition_state.py

@@ -36,7 +36,6 @@ def _run_opt_ts_calc(self, method, name_ext):
                                       other_input_block=method.keywords.optts_block)
         self.optts_calc.run()
 
-        # Was this intentionally removed?
         if not self.optts_calc.optimisation_converged():
             if self.optts_calc.optimisation_nearly_converged():
                 logger.info('Optimisation nearly converged')

autode/transition_states/ts_guess.py

@@ -57,12 +57,9 @@ def get_ts_guess_constrained_opt(reactant, method, keywords, name, distance_cons
     # corresponding energy
     try:
         mol_with_constraints.optimise(method=method, calc=hl_const_opt)
+
     except AtomsNotFound:
-        # Retrun with the low level
-        try:
-            mol_with_constraints.run_const_opt(ll_const_opt)
-        except AtomsNotFound:
-            return None
+        pass
 
     return get_ts_guess(species=mol_with_constraints, reactant=reactant,
                         product=product, name=f'ts_guess_{name}')

autode/wrappers/G09.py

@@ -1,5 +1,6 @@
 from copy import deepcopy
 import numpy as np
+from autode.constants import Constants
 from autode.wrappers.base import ElectronicStructureMethod
 from autode.utils import run_external
 from autode.atoms import Atom
@@ -440,16 +441,20 @@ def get_gradients(self, calc):
                     gradients_section = False
 
             if gradients_section and len(line.split()) == 5:
-                _, _, x, y, z = line.split()
+                _, _, fx, fy, fz = line.split()
                 try:
-                    gradients.append([float(x), float(y), float(z)])
+                    # Ha / a0
+                    force = np.array([float(fx), float(fy), float(fz)])
+
+                    grad = -force / Constants.a02ang
+                    gradients.append(grad)
                 except ValueError:
                     pass
         for line in gradients:
             for i in range(3):
                 line[i] *= -1
 
-        return gradients
+        return np.array(gradients)
 
     def __init__(self):
         super().__init__(name='g09', path=Config.G09.path,

autode/wrappers/MOPAC.py

@@ -293,7 +293,7 @@ def get_gradients(self, calc):
                 _, _, _, _, _, _, value, _ = line.split()
                 gradients.append(value)
         grad_array = np.asarray(gradients)
-        grad_array *= Constants.a02ang/Constants.ha2kcalmol
+        grad_array /= Constants.ha2kcalmol  # From kcal mol-1 Å^-1 to Ha Å^-1
         grad_array.reshape((calc.molecule.n_atoms, 3))
 
         return grad_array.tolist()

autode/wrappers/NWChem.py

@@ -356,9 +356,11 @@ def get_gradients(self, calc):
 
             if gradients_section and len(line.split()) == 8:
                 x, y, z = line.split()[5:]
-                gradients.append([float(x), float(y), float(z)])
+                gradients.append(np.array([float(x), float(y), float(z)]))
 
-        return gradients
+        # Convert from Ha a0^-1 to Ha A-1
+        gradients = [grad / Constants.a02ang for grad in gradients]
+        return np.array(gradients)
 
     def __init__(self):
         super().__init__('nwchem', path=Config.NWChem.path,

autode/wrappers/ORCA.py

@@ -1,5 +1,6 @@
 import numpy as np
 import os
+from autode.constants import Constants
 from autode.utils import run_external
 from autode.wrappers.base import ElectronicStructureMethod
 from autode.atoms import Atom
@@ -423,8 +424,12 @@ def get_gradients(self, calc):
                 first, last = i + 3, i + 3 + calc.molecule.n_atoms
                 for grad_line in calc.output.file_lines[first:last]:
                     dadx, dady, dadz = grad_line.split()[-3:]
-                    gradients.append([float(dadx), float(dady), float(dadz)])
+                    vec = [float(dadx), float(dady), float(dadz)]
 
+                    gradients.append(np.array(vec))
+
+        # Convert from Ha a0^-1 to Ha A-1
+        gradients = [grad / Constants.a02ang for grad in gradients]
         return np.array(gradients)
 
     def __init__(self):

autode/wrappers/XTB.py

@@ -130,7 +130,7 @@ def execute(self, calc):
             flags += ['--input', calc.input.additional_filenames[-1]]
 
         @work_in_tmp_dir(filenames_to_copy=calc.input.get_input_filenames(),
-                         kept_file_exts=('.xyz', '.out', '.pc'))
+                         kept_file_exts=('.xyz', '.out', '.pc', '.grad', 'gradient'))
         def execute_xtb():
             logger.info(f'Setting the number of OMP threads to {calc.n_cores}')
             os.environ['OMP_NUM_THREADS'] = str(calc.n_cores)
@@ -286,22 +286,32 @@ def get_atomic_charges(self, calc):
 
     def get_gradients(self, calc):
         gradients = []
+
         if os.path.exists(f'{calc.name}_xtb.grad'):
             grad_file_name = f'{calc.name}_xtb.grad'
             with open(grad_file_name, 'r') as grad_file:
                 for line in grad_file:
                     x, y, z = line.split()
-                    gradients.append([float(x), float(y), float(z)])
-        else:
+                    gradients.append(np.array([float(x), float(y), float(z)]))
+
+        elif os.path.exists('gradient'):
             with open('gradient', 'r') as grad_file:
                 for line_no, line in enumerate(grad_file):
                     if line_no > 1 and len(line.split()) == 3:
                         x, y, z = line.split()
-                        gradients.append([float(x.replace('D', 'E')), float(y.replace('D', 'E')), float(z.replace('D', 'E'))])
+                        vec = [float(x.replace('D', 'E')),
+                               float(y.replace('D', 'E')),
+                               float(z.replace('D', 'E'))]
+
+                        gradients.append(np.array(vec))
+
             with open(f'{calc.name}_xtb.grad', 'w') as new_grad_file:
-                [print('{:^12.8f} {:^12.8f} {:^12.8f}'.format(*line), file=new_grad_file) for line in gradients]
+                [print('{:^12.8f} {:^12.8f} {:^12.8f}'.format(*line),
+                       file=new_grad_file) for line in gradients]
             os.remove('gradient')
 
+        # Convert from Ha a0^-1 to Ha A-1
+        gradients = [grad / Constants.a02ang for grad in gradients]
         return np.array(gradients)
 
     def __init__(self):

doc/examples/conformers.rst

@@ -60,7 +60,7 @@ Metal Complex
 
 Arbitrary distance constraints can be added in a RR conformer generation. For
 example, to generate conformers of
-`Vaska's complex <https://en.wikipedia.org/wiki/Vaska%27s_complex/>`_
+`Vaska's complex <https://en.wikipedia.org/wiki/Vaska%27s_complex>`_
 while retaining the square planar geometry
 
 

doc/examples/molecules.rst

@@ -104,7 +104,8 @@ Calculations
 
 **autodE** provides wrappers around common electronic structure theory packages
 (ORCA, XTB, NWChem, MOPAC, Gaussian09) so geometries may be optimised and
-energies calculated.
+energies calculated. Energies are in atomic Hartrees and gradients in
+Ha / Å.
 
 For example, to optimise the geometry at the XTB level and then perform a
 single point energy evaluation with ORCA