@@ -75,6 +75,9 @@
Loading
75 75
class MolproEngineError(EngineError):
76 76
    pass
77 77
78 +
class GaussianEngineError(EngineError):
79 +
    pass
80 +
78 81
class QCEngineAPIEngineError(EngineError):
79 82
    pass
80 83

@@ -38,10 +38,11 @@
Loading
38 38
import os
39 39
import itertools
40 40
import numpy as np
41 +
import shutil
41 42
42 43
import os
43 44
from .internal import Distance, Angle, Dihedral, CartesianX, CartesianY, CartesianZ, TranslationX, TranslationY, TranslationZ, RotationA, RotationB, RotationC
44 -
from .engine import set_tcenv, load_tcin, TeraChem, ConicalIntersection, Psi4, QChem, Gromacs, Molpro, OpenMM, QCEngineAPI
45 +
from .engine import set_tcenv, load_tcin, TeraChem, ConicalIntersection, Psi4, QChem, Gromacs, Molpro, OpenMM, QCEngineAPI, Gaussian
45 46
from .rotate import calc_fac_dfac
46 47
from .molecule import Molecule, Elements
47 48
from .nifty import logger, isint, uncommadash, bohr2ang, ang2bohr
@@ -119,8 +120,8 @@
Loading
119 120
    if engine_str:
120 121
        engine_str = engine_str.lower()
121 122
        if engine_str[:4] == 'tera' : engine_str = 'tera'
122 -
        if engine_str not in ['tera', 'qchem', 'psi4', 'gmx', 'molpro', 'openmm', 'qcengine']:
123 -
            raise RuntimeError("Valid values of engine are: tera, qchem, psi4, gmx, molpro, openmm, qcengine")
123 +
        if engine_str not in ['tera', 'qchem', 'psi4', 'gmx', 'molpro', 'openmm', 'qcengine', "gaussian"]:
124 +
            raise RuntimeError("Valid values of engine are: tera, qchem, psi4, gmx, molpro, openmm, qcengine, gaussian")
124 125
        if customengine:
125 126
            raise RuntimeError("engine and customengine cannot simultaneously be set")
126 127
        if engine_str == 'tera':
@@ -218,6 +219,21 @@
Loading
218 219
            if molproexe is not None:
219 220
                engine.set_molproexe(molproexe)
220 221
            threads_enabled = True
222 +
        elif engine_str == "gaussian":
223 +
            logger.info("Gaussian engine selected. Expecting Gaussian input for gradient calculation. \n")
224 +
            M = Molecule(inputf, radii=radii, fragment=frag)
225 +
            # now work out which gaussian version we have
226 +
            if shutil.which("g16") is not None:
227 +
                exe = "g16"
228 +
            elif shutil.which("g09") is not None:
229 +
                exe = "g09"
230 +
            else:
231 +
                raise ValueError("Neither g16 or g09 was found, please check the environment.")
232 +
            engine = Gaussian(molecule=M, exe=exe, threads=threads)
233 +
            threads_enabled = True
234 +
            logger.info("The gaussian engine exe is set as %s" % engine.gaussian_exe)
235 +
            # load the template into the engine
236 +
            engine.load_gaussian_input(inputf)
221 237
        elif engine_str == 'qcengine':
222 238
            logger.info("QCEngine selected.\n")
223 239
            schema = kwargs.get('qcschema', False)

@@ -240,7 +240,8 @@
Loading
240 240
    grp_univ.add_argument('--engine', type=str, help='Specify engine for computing energies and gradients.\n'
241 241
                          '"tera" = TeraChem (default)         "qchem" = Q-Chem\n'
242 242
                          '"psi4" = Psi4                       "openmm" = OpenMM (pass a force field or XML input file)\n'
243 -
                          '"molpro" = Molpro                   "gmx" = Gromacs (pass conf.gro; requires topol.top and shot.mdp\n ')
243 +
                          '"molpro" = Molpro                   "gmx" = Gromacs (pass conf.gro; requires topol.top and shot.mdp\n '
244 +
                          '"gaussian" = Gaussian09/16\n ')
244 245
    grp_univ.add_argument('--nt', type=int, help='Specify number of threads for running in parallel\n(for TeraChem this should be number of GPUs)')
245 246
    
246 247
    grp_jobtype = parser.add_argument_group('jobtype', 'Control the type of optimization job')

@@ -3182,16 +3182,18 @@
Loading
3182 3182
            line = line.strip().expandtabs()
3183 3183
            # Everything after exclamation point is a comment
3184 3184
            sline = line.split('!')[0].split()
3185 -
            if len(sline) == 2:
3185 +
            if re.match(r"^ *[A-Z][a-z]?(.*[-+]?([0-9]*\.)?[0-9]+){3}$", line) is not None:
3186 +
                inxyz = 1
3187 +
                if sline[0].capitalize() in PeriodicTable and isfloat(sline[1]) and isfloat(sline[2]) and isfloat(
3188 +
                        sline[3]):
3189 +
                    elem.append(sline[0])
3190 +
                    xyz.append(np.array([float(sline[1]), float(sline[2]), float(sline[3])]))
3191 +
3192 +
            elif len(sline) == 2:
3186 3193
                if isint(sline[0]) and isint(sline[1]):
3187 3194
                    charge = int(sline[0])
3188 3195
                    mult = int(sline[1])
3189 3196
                    title_ln = ln - 2
3190 -
            elif len(sline) == 4:
3191 -
                inxyz = 1
3192 -
                if sline[0].capitalize() in PeriodicTable and isfloat(sline[1]) and isfloat(sline[2]) and isfloat(sline[3]):
3193 -
                    elem.append(sline[0])
3194 -
                    xyz.append(np.array([float(sline[1]),float(sline[2]),float(sline[3])]))
3195 3197
            elif inxyz:
3196 3198
                break
3197 3199
            ln += 1

@@ -49,7 +49,7 @@
Loading
49 49
from .molecule import Molecule
50 50
from .nifty import bak, eqcgmx, fqcgmx, bohr2ang, logger, getWorkQueue, queue_up_src_dest, rootdir, splitall, copy_tree_over
51 51
from .errors import EngineError, CheckCoordError, Psi4EngineError, QChemEngineError, TeraChemEngineError, \
52 -
    ConicalIntersectionEngineError, OpenMMEngineError, GromacsEngineError, MolproEngineError, QCEngineAPIEngineError
52 +
    ConicalIntersectionEngineError, OpenMMEngineError, GromacsEngineError, MolproEngineError, QCEngineAPIEngineError, GaussianEngineError
53 53
54 54
#=============================#
55 55
#| Useful TeraChem functions |#
@@ -686,6 +686,147 @@
Loading
686 686
    def copy_scratch(self, src, dest):
687 687
        return
688 688
689 +
class Gaussian(Engine):
690 +
    """
691 +
    Run a Gaussian energy and gradient calculation.
692 +
    """
693 +
    def __init__(self, molecule, exe=None, threads=None):
694 +
        super(Gaussian, self).__init__(molecule)
695 +
        self.threads = threads
696 +
        if exe.lower() in ("g16", "g09"):
697 +
            self.gaussian_exe = exe.lower()
698 +
        else:
699 +
            raise ValueError("Only g16 and g09 are supported.")
700 +
701 +
    def load_gaussian_input(self, gaussian_input):
702 +
        """
703 +
        We can read the .com files using molecule but we can not write them so use the template method.
704 +
705 +
        Note only Gaussian cartesian coordinates are supported.
706 +
        We also edit the checkpoint file name and change the processors flag if not present.
707 +
708 +
        Example input file:
709 +
710 +
        %Mem=6GB
711 +
        %NProcShared=2
712 +
        %Chk=lig
713 +
        # B3LYP/6-31G(d) Force=NoStep
714 +
715 +
        water energy
716 +
717 +
        0   1
718 +
        O  -0.464   0.177   0.0
719 +
        H  -0.464   1.137   0.0
720 +
        H   0.441  -0.143   0.0
721 +
722 +
723 +
        """
724 +
        reading_molecule, found_geo, found_force = False, False, False
725 +
        gauss_temp = []  # store a template of the input file for generating new ones
726 +
        with open(gaussian_input) as gauss_in:
727 +
            for line in gauss_in:
728 +
                match = re.search(r"^ *[A-Z][a-z]?(.*[-+]?([0-9]*\.)?[0-9]+){3}$", line)
729 +
                if match is not None:
730 +
                    reading_molecule = True
731 +
                    if not found_geo:
732 +
                        found_geo = True
733 +
                        gauss_temp.append("$!geometry@here")
734 +
735 +
                elif reading_molecule:
736 +
                    if line.strip() == '':
737 +
                        reading_molecule = False
738 +
                        gauss_temp.append(line)
739 +
740 +
                elif "%nprocshared" in line.lower():
741 +
                    # we should replace the line with our threads value
742 +
                    if self.threads is not None:
743 +
                        gauss_temp.append("%NProcShared=" + str(self.threads) + "\n")
744 +
                    else:
745 +
                        gauss_temp.append(line)
746 +
747 +
                elif "%chk" in line.lower():
748 +
                    # we should replace it with what we want
749 +
                    gauss_temp.append("%Chk=ligand\n")
750 +
751 +
                else:
752 +
                    gauss_temp.append(line)
753 +
754 +
                if "force=nostep" in line.lower():
755 +
                    found_force = True
756 +
        if not found_force:
757 +
            raise RuntimeError("Gaussian inputfile %s should have force=nostep command." % gaussian_input)
758 +
759 +
        # now we need to make sure the chk point file and threads were set
760 +
        if not any("%chk" in command.lower() for command in gauss_temp):
761 +
            # insert at the top
762 +
            gauss_temp.insert(0, "%Chk=ligand\n")
763 +
        if not any("%nprocshared" in command.lower() for command in gauss_temp):
764 +
            # if threads is not none set it else use 1 thread
765 +
            thread_str = "%NProcShared=" + str(self.threads or 1)
766 +
            gauss_temp.insert(0, thread_str + "\n")
767 +
768 +
        self.gauss_temp = gauss_temp
769 +
770 +
    def calc_new(self, coords, dirname):
771 +
        """
772 +
        Run the gaussian single point calculation using the given exe.
773 +
        """
774 +
        if not os.path.exists(dirname): os.makedirs(dirname)
775 +
        # Convert coordinates back to the xyz file
776 +
        self.M.xyzs[0] = coords.reshape(-1, 3) * bohr2ang
777 +
        # Write Gaussian com file
778 +
        with open(os.path.join(dirname, 'gaussian.com'), 'w') as outfile:
779 +
            for line in self.gauss_temp:
780 +
                if line == '$!geometry@here':
781 +
                    for i, (e, c) in enumerate(zip(self.M.elem, self.M.xyzs[0])):
782 +
                        outfile.write("%-7s %13.7f %13.7f %13.7f\n" % (e, c[0], c[1], c[2]))
783 +
                else:
784 +
                    outfile.write(line)
785 +
        try:
786 +
            # Run Gaussian
787 +
            subprocess.check_call('%s < gaussian.com > gaussian.log && formchk ligand.chk ligand.fchk > form_log.txt' % self.gaussian_exe, cwd=dirname, shell=True)
788 +
            # Read energy and gradients from Gaussian output
789 +
            result = self.read_result(dirname)
790 +
        except (OSError, IOError, RuntimeError, subprocess.CalledProcessError):
791 +
            raise GaussianEngineError
792 +
        return result
793 +
794 +
    def read_result(self, dirname, check_coord=None):
795 +
        """
796 +
        Read the result of the output file to get the gradient and energy.
797 +
        """
798 +
        if check_coord is not None:
799 +
            raise CheckCoordError("Coordinate checking not implemented")
800 +
        energy, gradient = None, None
801 +
        # first get the energy from the formatted checkpoint file, works for all methods
802 +
        fchk_out = os.path.join(dirname, "ligand.fchk")
803 +
        with open(fchk_out) as fchk:
804 +
            for line in fchk:
805 +
                if "Total Energy" in line:
806 +
                    energy = float(line.split()[-1])
807 +
                    break
808 +
        # now we get the gradient from the output in Hartrees/Bohr
809 +
        gaussian_out = os.path.join(dirname, 'gaussian.log')
810 +
        with open(gaussian_out) as outfile:
811 +
            found_grad = False
812 +
            for line in outfile:
813 +
                line_strip = line.strip()
814 +
                if " Forces (Hartrees/Bohr)" in line_strip:
815 +
                    found_grad = True
816 +
                    gradient = []
817 +
                elif found_grad:
818 +
                    ls = line_strip.split()
819 +
                    if len(ls) == 5 and ls[0].isdigit() and ls[1].isdigit():
820 +
                        gradient.append([-float(g) for g in ls[2:]])
821 +
                    elif "Cartesian Forces:  Max" in line:
822 +
                        found_grad = False
823 +
        if energy is None:
824 +
            raise RuntimeError("Gaussian energy is not found in %s, please check." % fchk_out)
825 +
        if gradient is None:
826 +
            raise RuntimeError("Gaussian gradient is not found in %s, please check." % gaussian_out)
827 +
        gradient = np.array(gradient, dtype=np.float64).ravel()
828 +
        return {'energy':energy, 'gradient':gradient}
829 +
689 830
class Psi4(Engine):
690 831
    """
691 832
    Run a Psi4 energy and gradient calculation.
Files Coverage
geometric 69.90%
Project Totals (15 files) 69.90%
551.2
TRAVIS_PYTHON_VERSION=3.5
TRAVIS_OS_NAME=linux
551.1
TRAVIS_PYTHON_VERSION=2.7
TRAVIS_OS_NAME=linux
551.3
TRAVIS_PYTHON_VERSION=3.6
TRAVIS_OS_NAME=linux
551.4
TRAVIS_PYTHON_VERSION=3.7
TRAVIS_OS_NAME=linux
1
coverage:
2
  ignore:
3
    - geometric/tests/.*
4
    - geometric/PDB.py
5
    - geometric/_version.py
6
    - docs/*
7
    - setup.py
8
  status:
9
    patch: false
10
comment:
11
  layout: "header"
12
  require_changes: false
13
  branches: null
14
  behavior: default
15
  flags: null
16
  paths: null
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading