updates mainly to runjob and the Molecule object to reduce dependenci…

…es and run steps

qforce/initialize.py

@@ -82,7 +82,6 @@ def _get_job_info(filename):
     job['dir'] = f'{path}{job["name"]}_qforce'
     pathways = Pathways(job['dir'], name=job['name'])
     job['pathways'] = pathways
-    job['frag_dir'] = str(pathways["fragments"])
     #
     if job['coord_file'] is False:
         init = pathways['init.xyz']

qforce/main.py

@@ -12,76 +12,77 @@
 from .logger import LoggerExit
 #
 from .schemes import Computations, HessianCreator, CrestCreator
-from .schemes.bondorder import BondOrderCreator
-from .schemes import NoFragmentationDihedralsCreator
+from .schemes import DihedralCreator
 
 
-def save_molecule(filename, molecule, comment=None):
-    write(filename, molecule, plain=True, comment=comment)
+def runjob(config, job, ext_q=None, ext_lj=None):
 
+    qm_interface = QM(job, config.qm)
+    ff_interface = ForceField.implemented_md_software[config.ff.output_software]
 
-def runjob_(config, job, ext_q=None, ext_lj=None):
-    # setup the qm interface
-    qm = QM(job, config.qm)
-    # read initial structure, ase molecule structure
-    molecule = read(job.coord_file)
-    # run preopt
-    folder = job.pathways.getdir('preopt', create=True)
-    crest = CrestCreator(folder, molecule)
-    crest.run(qm)
-    molecule.set_positions(crest.most_stable())
-    # setup hessian calculation
-    save_molecule(qm.pathways.get('init.xyz'), molecule,
-                  comment=f'{job.name} - input geometry for hessian')
-    # get the bondorder to setup molecule
-    bo = BondOrderCreator(molecule)
-    bo.run(qm)
-    bondorder = bo.bondorder()
-    #
-    ffcls = ForceField.implemented_md_software.get(config.ff.output_software, None)
-    if ffcls is None:
-        raise ValueError(f"Forcefield '{config.ff.output_software}' unknown!")
+    mol = Molecule(job, config)
 
-    mol = Molecule(config, job, bondorder, ffcls, ext_q, ext_lj)
-    #
-    folder = job.pathways.getdir('addstruct')
-    structs = Computations(config.addstructs, folder)
-    structs.register('dihedrals', NoFragmentationDihedralsCreator(mol, job, config))
-    structs.register('hessian', HessianCreator(molecule))
-    #
-    structs.activate('fromfile')
-    structs.activate('xtbmd', crest.structures())
-    # do all additional calculations
-    structs.run(qm)
-    #  register hessian after structs were run!
-    # structs.register('hessian', hessian)
-    #
-    mol.qm_minimum_energy, mol.qm_minimum_coords = structs.normalize()
-    #
-    md_hessian = multi_fit(job.logger, config.terms, mol, structs)
-    #
-    return bondorder, md_hessian, mol, structs
+    do_crest(job, qm_interface, mol)
 
+    # This is ideally temporary if we can fix the global optimization
+    hessian_out = do_hessian(qm_interface, mol)
 
-def runjob(config, job, ext_q=None, ext_lj=None):
-    main_hessian, md_hessian, mol, structs = runjob_(config, job, ext_q=ext_q, ext_lj=ext_lj)
+    mol.setup(config, job, ff_interface, hessian_out, ext_q, ext_lj)
+
+    structs = do_all_structs(job, config, qm_interface, mol)
 
-    calc_qm_vs_md_frequencies(job, main_hessian, md_hessian)
+    md_hessian = multi_fit(job.logger, config.terms, mol, structs)
+    calc_qm_vs_md_frequencies(job, hessian_out, md_hessian)
 
-    if (main_hessian.dipole_deriv is not None
+    if (hessian_out.dipole_deriv is not None
        and 'charge_flux' in mol.terms
        and len(mol.terms['charge_flux']) > 0):
         raise NotImplementedError("Charge flux is not updated to new syntax")
         # fit_charge_flux(main_hessian, qm_energy_out, qm_gradient_out, mol)
 
     ff = ForceField(config.ff.output_software, job, config, mol, mol.topo.neighbors)
-    ff.software.write(job.dir, main_hessian.coords)
+    ff.software.write(job.dir, mol.coords)
 
     print_outcome(job.logger, job.dir, config.ff.output_software)
 
     return mol
 
 
+def do_hessian(qm_interface, mol):
+    hessian = HessianCreator(mol)
+    hessian.run(qm_interface)
+    main_hessian = hessian.main_hessian()
+    mol.update_coords(main_hessian.coords, 'Structure optimized for Hessian calculation')
+    return main_hessian
+
+
+def do_crest(job, qm_interface, mol):
+    folder = job.pathways.getdir('preopt', create=True)
+    crest = CrestCreator(folder, mol)
+    crest.run(qm_interface)
+    mol.update_coords(crest.get_most_stable(), 'CREST lowest energy structure')
+    mol.all_coords = crest.get_structures()
+    mol.bond_orders = crest.get_bond_orders()
+
+
+def do_all_structs(job, config, qm_interface, mol):
+    folder = job.pathways.jobdir
+
+    structs = Computations(config.addstructs, folder)
+    structs.register('dihedrals', DihedralCreator(mol, job, config))
+    structs.register('hessian', HessianCreator(mol))
+    #
+    structs.activate('fromfile')
+    structs.activate('xtbmd', mol.all_coords)
+    # do all additional calculations
+    structs.run(qm_interface)
+    #  register hessian after structs were run!
+    # structs.register('hessian', hessian)
+    #
+    mol.qm_minimum_energy, mol.qm_minimum_coords = structs.normalize()
+    return structs
+
+
 def load_keeper(job):
     file = job.pathways['calculations.json']
     if file.exists():

qforce/molecule/molecule.py

@@ -1,19 +1,48 @@
+from ase import Atoms
+from ase.io import read, write
+import numpy as np
+#
 from .topology import Topology
 from .terms import Terms
 from .non_bonded import NonBonded
 
 
 class Molecule(object):
 
-    def __init__(self, config, job, qm_out, ff, ext_q=None, ext_lj=None):
+    def __init__(self, job, config):
         self.name = job.name
-        self.atomids = qm_out.atomids
-        self.coords = qm_out.coords
-        self.charge = qm_out.charge
-        self.multiplicity = qm_out.multiplicity
+        self.job_dir = job.dir
+        self.charge = config.qm.charge
+        self.multiplicity = config.qm.multiplicity
+
+        self.coords = None
+        self.bond_orders = None
+        self.point_charges = None
+
+        coords, self.atomids = self._read_input_coords(job.coord_file)
+        self.update_coords(coords, 'Input Coordinates')
+        self.all_coords = [self.coords]
         self.n_atoms = len(self.atomids)
-        self.topo = Topology(config.ff, qm_out)
-        self.non_bonded = NonBonded.from_topology(config.ff, job, qm_out, self.topo, ext_q, ext_lj)
-        self.terms = Terms.from_topology(config.terms, self.topo, self.non_bonded, ff)
+
+        self.topo = None
+        self.non_bonded = None
+        self.terms = None
+
         self.qm_minimum_energy = None
         self.qm_minimum_coords = None
+
+    def _read_input_coords(self, file):
+        ase_molecule = read(file)
+        init_coords = ase_molecule.get_positions()
+        atomids = ase_molecule.get_atomic_numbers()
+        return init_coords, atomids
+
+    def update_coords(self, coords, comment=''):
+        self.coords = np.array(coords)
+        atoms = Atoms(positions=coords, numbers=self.atomids)
+        write(self.job_dir+'/coords.xyz', atoms, plain=True, comment=comment)
+
+    def setup(self, config, job, ff_interface, hessian_out, ext_q, ext_lj):
+        self.topo = Topology(config.ff, self)
+        self.non_bonded = NonBonded.from_topology(config.ff, job, hessian_out, self.topo, ext_q, ext_lj)
+        self.terms = Terms.from_topology(config.terms, self.topo, self.non_bonded, ff_interface)

qforce/molecule/topology.py

@@ -10,12 +10,12 @@ class Topology(object):
     Contains all bonding etc. information of the system
     """
 
-    def __init__(self, config, qm_out):
+    def __init__(self, config, mol):
         self.n_equiv = config.n_equiv
-        self.atomids = qm_out.atomids
+        self.atomids = mol.atomids
         self.n_atoms = len(self.atomids)
-        self.coords = qm_out.coords
-        self.b_order_matrix = qm_out.b_orders
+        self.coords = mol.coords
+        self.bond_order_matrix = mol.bond_orders
         #
         self.n_types = 0
         self.n_terms = 0
@@ -29,25 +29,24 @@ def __init__(self, config, qm_out):
         self.all_rigid = config.all_rigid
         self.ba_couple_1_shared = config.ba_couple_1_shared
         #
-        self._setup(qm_out)
+        self._setup()
 
-    def _setup(self, qm_out):
-        self._find_bonds_and_rings(qm_out)
+    def _setup(self):
+        self._find_bonds_and_rings()
         self._find_atom_types()
         self._find_neighbors()
         self._find_bonds_angles_dihedrals()
 
-    def _find_bonds_and_rings(self, qm_out):
+    def _find_bonds_and_rings(self):
         """Setup networkx graph """
         self.graph = nx.Graph()
         for i_idx, i_elem in enumerate(self.atomids):
-            self.graph.add_node(i_idx, idx=i_idx, elem=i_elem, n_bonds=qm_out.n_bonds[i_idx],
-                                q=qm_out.point_charges[i_idx], coords=self.coords[i_idx],
+            self.graph.add_node(i_idx, idx=i_idx, elem=i_elem, coords=self.coords[i_idx],
                                 neighs=[], unique_neighs=[], nonrepeat_neighs=[], hybrid=None,
                                 type=None, n_neighs=0, n_unique_neighs=0, n_nonrepeat_neighs=0)
             # add bonds
             for j_idx, j_elem in enumerate(self.atomids):
-                b_order = qm_out.b_orders[i_idx, j_idx]
+                b_order = self.bond_order_matrix[i_idx, j_idx]
                 if b_order > 0.3:
                     id1, id2 = sorted([i_elem, j_elem])
                     b_order_half_rounded = np.round(b_order*2)/2
@@ -58,10 +57,12 @@ def _find_bonds_and_rings(self, qm_out):
                     self.graph.add_edge(i_idx, j_idx, vector=vec, length=dist, order=b_order, vers=None,
                                         type=f'{id1}({b_order_half_rounded}){id2}', n_rings=0)
 
-            if qm_out.n_bonds[i_idx] > ELE_MAXB[i_elem]:
+            n_bonds = self.bond_order_matrix[i_idx].sum().round()
+            self.node(i_idx)['n_bonds'] = n_bonds
+            if n_bonds > ELE_MAXB[i_elem]:
                 print(f'WARNING: Atom {i_idx+1} ({ATOM_SYM[i_elem]}) has too many',
-                      f' ({qm_out.n_bonds[i_idx]}) bonds?')
-            elif qm_out.n_bonds[i_idx] == 0:
+                      f' ({n_bonds}) bonds?')
+            elif n_bonds == 0:
                 print(f'WARNING: Atom {i_idx+1} ({ATOM_SYM[i_elem]}) has no bonds')
 
             if self.node(i_idx)['n_neighs'] == 1:

qforce/pathkeeper.py

@@ -109,14 +109,14 @@ class Pathways:
         # folders
         'preopt': 'crest',
         'hessian': 'hessian',
-        'addstruct': '5_additional',
+        'addstruct': 'additional',
         'hessian_new': 'hessian_new',
         'hessian_charge': ['hessian', 'charge'],
         'hessian_energy': ['hessian', '${idx}_en_conformer'],
         'hessian_gradient': ['hessian', '${idx}_grad_conformer'],
         'hessian_step': ['hessian', '${idx}_conformer'],
-        'fragments': '2_fragments',
-        'frag': ['fragments', '${frag_id}'],
+        'dihedrals': 'dihedrals',
+        'frag': ['dihedrals', '${frag_id}'],
         'frag_charge': ['frag', 'charge'],
         'frag_step': ['frag', 'step_${idx}'],
         'frag_mm': ['frag', 'mm'],

qforce/qm/crest.py

@@ -1,7 +1,8 @@
 import os.path
-
 from ase.io import read, write
 from ase import Atoms
+import numpy as np
+#
 from .qm_base import WriteABC, ReadABC, QMInterface, Calculator
 
 
@@ -39,16 +40,19 @@ class ReadCrest(ReadABC):
 
     hessian_files = {}
 
-    opt_files = {'coord_file': ['crest_conformers.xyz'], }
+    opt_files = {'coord_file': ['crest_conformers.xyz'], 'wbo_file': ['wbo']}
 
     sp_files = {}
 
     charge_files = {}
 
     scan_files = {}
 
-    def opt(self, config, coord_file):
-        return self._read_xyzs(coord_file)
+    def opt(self, config, coord_file, wbo_file):
+        coords = self._read_xyzs(coord_file)
+        n_atoms = len(coords[0])
+        bond_orders = self._read_crest_wbo_analysis(wbo_file, n_atoms)
+        return coords, np.array(bond_orders)
 
     def sp(self, config, sp_file):
         raise NotImplementedError
@@ -146,6 +150,34 @@ def _read_xyzs(coord_file):
         mols = read(coord_file, index=':')
         return [mol.get_positions() for mol in mols]
 
+    @staticmethod
+    def _read_crest_wbo_analysis(out_file, n_atoms):
+        """ Read the wbo analysis from CREST.
+
+        Parameters
+        ----------
+        out_file : string
+            The name of the orca output file.
+        n_atoms : int
+            The number of atoms in the molecule.
+
+        Returns
+        -------
+        b_orders : list
+            A list (length: n_atoms) of list (length: n_atoms) of float.
+            representing the bond order between each atom pair.
+
+        """
+
+        b_orders = [[0 for _ in range(n_atoms)] for _ in range(n_atoms)]
+
+        file = np.loadtxt(out_file)
+        for x, y, bo in file:
+            b_orders[int(x) - 1][int(y) - 1] = bo
+            b_orders[int(y) - 1][int(x) - 1] = bo
+
+        return b_orders
+
 
 class WriteCrest(WriteABC):
 

qforce/qm/qm.py

@@ -161,10 +161,10 @@ def Calculation(self, filename, required_files, *, folder=None, software=None):
 
     def setup_hessian_calculation(self, folder, coords, atnums, preopt=False):
         """Setup hessian calculation"""
-        if preopt is True:
-            software = self.softwares['preopt']
-        else:
-            software = self.softwares['software']
+        # if preopt is True:
+        #     software = self.softwares['preopt']
+        # else:
+        software = self.softwares['software']
 
         calculation = self.Calculation(self.hessian_name(software),
                                        software.read.hessian_files,
@@ -360,7 +360,7 @@ def read_scan_data(self, files):
         software = self.softwares['scan_software']
 
         if self.config.dihedral_scanner == 'relaxed_scan':
-            out = software.read.scan_new(self.config, **files)
+            out = software.read.scan(self.config, **files)
         elif self.config.dihedral_scanner == 'torsiondrive':
             out = software.read.scan_torsiondrive(self.config, **files)