Initial add

.gitignore

@@ -0,0 +1,30 @@
+
+\.idea/
+
+__pycache__/
+
+\.DS_Store
+
+dist/
+
+cgbind\.egg-info/
+
+build/
+
+tests/cage_L-1\.xyz
+
+tests/gen_cage/cage_L-1\.xyz
+
+tests/gen_cage/cage_L-1_orca_sp\.inp
+
+tests/gen_cage_subst/cage_L-1\.xyz
+
+tests/gen_cage_subst/cage_L-1_acetone\.xyz
+
+tests/gen_cage_subst/cage_L-1_bq\.xyz
+
+tests/gen_cage_subst/cage_L-1_methane\.xyz
+
+tests/gen_cage/cage_L_m4l6\.xyz
+
+tests/gen_cage/tmp_cage\.xyz

LICENSE.md

@@ -0,0 +1,22 @@
+
+The MIT License (MIT)
+
+Copyright (c) 2019 
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -1 +1,35 @@
-cgbind
+![alt text](cgbind/common/llogo.png)
+
+### *C*a*g*e*Bind*ing: A tool to for automated binding affinity calculations
+
+Constructing metallocage structures 
+
+
+## Requirements
+0. Python v. 3.7
+1. [rdkit](https://github.com/rdkit/rdkit)
+2. matplotlib
+3. ORCA v. 4.1 (optional)
+4. xTB (optional)
+5. MOPAC (optional)
+
+## Installation
+
+If the requirements are already satisfied
+```
+python setup.py install
+```
+
+Alternately for a one line fresh install with miniconda on Mac OSX
+```
+wget https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh && bash Miniconda3-latest-MacOSX-x86_64.sh && source ~/.bash_profile && conda install -c rdkit rdkit && git clone https://[email protected]/duartegroup/cgbind.git && cd cgbind && python setup.py install
+```
+
+
+## Usage
+To load as a module:
+```
+import cgbind
+```
+
+See _examples/_ for how to use _cgbind_

cgbind/MOPACio.py

@@ -0,0 +1,120 @@
+import os
+from .log import logger
+from .config import Config
+from .constants import Constants
+from subprocess import Popen
+
+
+def singlepoint(xyzs, name, charge=0, mult=1):
+    """
+    Run a single point calculation using MOPAC
+    :param xyzs: list of xyzs
+    :param name: name of the species, will be appended with '_mopac_sp' to form the filename
+    :param charge: total charge of the species
+    :param mult: multiplicity of the species
+    :return:
+    """
+    logger.info('Running MOPAC single point')
+
+    energy = 0
+    mop_filename = name + '_mopac_sp.mop'
+    mop_out_filename = mop_filename.replace('.mop', '.out')
+
+    gen_mopac_mop(mop_filename, xyzs, charge, mult, opt=False)
+    if not Config.suppress_print:
+        print("{:<30s}{:<50s}{:>10s}".format('Single point calculation of', name, 'Running'))
+    mopac_out_lines = run_mopac(mop_filename, mop_out_filename)
+    if not Config.suppress_print:
+        print("{:<30s}{:<50s}{:>10s}".format('', name, 'Done'))
+
+    for line in mopac_out_lines:
+        if 'TOTAL ENERGY' in line:
+            energy_ev = float(line.split()[-2])
+            energy = energy_ev * Constants.eV2ha
+
+    if energy == 0:
+        logger.warning('Didn\'t find an energy in MOPAC output')
+
+    return energy
+
+
+def gen_mopac_mop(mop_filename, xyzs, charge, mult, opt=True, opt_atom_ids=None):
+    logger.info('Generating a MOPAC input file')
+
+    keywords = ['AUX', 'THREADS=1', 'CHARGE=' + str(charge)]
+    if mult == 1:
+        keywords.append('SINGLET')
+    if mult == 2:
+        keywords.append('DOUBLET')
+    if mult == 3:
+        keywords.append('TRIPLET')
+    if Config.opt_method == 'pm6':
+        keywords.append('PM6')
+    elif Config.opt_method == 'pm7':
+        keywords.append('PM7')
+    else:
+        keywords.append('PM7')
+    if not opt:
+        keywords.append('1SCF')
+
+    with open(mop_filename, 'w') as mopac_mop:
+        print(*keywords, '\nTitle \n', file=mopac_mop)
+        if opt_atom_ids:
+            for i in range(len(xyzs)):
+                line = xyzs[i]
+                if i in opt_atom_ids:
+                    print(line[0], line[1], '1', line[2], '1', line[3], '1', sep='\t', file=mopac_mop)
+                else:
+                    print(line[0], line[1], '0', line[2], '0', line[3], '0', sep='\t', file=mopac_mop)
+        else:
+            [print(line[0], line[1], '1', line[2], '1', line[3], '1', sep='\t', file=mopac_mop) for line in xyzs]
+
+    return 0
+
+
+def run_mopac(mop_filename, out_filename):
+
+    dev_null = open(os.devnull, 'w')
+    mopac_run = Popen([Config.path_to_mopac, mop_filename],
+                      stderr=dev_null,
+                      env={'MOPAC_LICENSE': Config.path_to_mopac_licence})
+    mopac_run.wait()
+
+    return [line for line in open(out_filename, 'r')]
+
+
+def get_mopac_xyzs_energy(out_lines):
+
+    opt_converged, geom_section = False, False
+    opt_xyzs = []
+    energy = 0.0
+    n_blank_lines = 0
+
+    for line in out_lines:
+
+        if 'SCF FIELD WAS ACHIEVED' in line:
+            opt_converged = True
+
+        if 'CARTESIAN COORDINATES' in line and opt_converged:
+            geom_section = True
+
+        if geom_section and len(line.split()) == 0:
+            n_blank_lines += 1
+
+        if geom_section and n_blank_lines == 2:
+            geom_section = False
+
+        if geom_section and len(line.split()) == 5:
+            atom_label, x, y, z = line.split()[1:]
+            opt_xyzs.append([atom_label, float(x), float(y), float(z)])
+
+        if 'TOTAL ENERGY' in line:
+            energy_ev = float(line.split()[-2])
+            energy = Constants.eV2ha * energy_ev
+
+    if len(opt_xyzs) == 0:
+        logger.warning('Could not get optimised xyzs from MOPAC output')
+    if energy == 0:
+        logger.warning('Could not get optimised energy from MOPAC output')
+
+    return opt_xyzs, energy

cgbind/ORCAio.py

@@ -0,0 +1,153 @@
+import os
+from .log import logger
+from .config import Config
+from subprocess import Popen
+
+
+def singlepointenergy(xyzs, name, keywords, solvent=None, charge=0, mult=1, n_cores=1):
+    logger.info('Running an ORCA single point of {}'.format(name))
+
+    if not Config.code == 'orca':
+        logger.warning('Code was set at "{}", but only ORCA single points are supported, '
+                       'Using ORCA anyway'.format(Config.code))
+
+    energy = 0
+    inp_filename = name + '_orca_sp.inp'
+    out_filename = inp_filename.replace('.inp', '.out')
+
+    if ('PAL' + str(n_cores)) not in keywords and n_cores > 1:
+        if n_cores > 8:
+            logger.warning("Number of ORCA cores requested was {}. Reducing to 8".format(n_cores))
+            keywords.append('PAL8')
+        else:
+            keywords.append('PAL' + str(n_cores))
+
+    if solvent and 'CPCM' not in keywords:
+        keywords.append('CPCM')
+
+    xys2inp(xyzs, keywords, inp_filename, charge, mult, solvent)
+    if Config.path_to_orca is None:
+        logger.error('path_to_orca needs to be set for an ORCA single point. Skipping')
+        return
+
+    if not Config.suppress_print:
+        print("{:<30s}{:<50s}{:>10s}".format('Single point calculation of', name, 'Running'))
+    orca_out_lines = run_orca(inp_filename, out_filename)
+    if not Config.suppress_print:
+        print("{:<30s}{:<50s}{:>10s}".format(' ', name, 'Done'))
+    for line in orca_out_lines:
+        if 'FINAL SINGLE POINT ENERGY' in line:
+            energy = float(line.split()[4])
+
+    return energy
+
+
+def xys2inp(xyzs, keywords, filename, charge=0, mult=1, solvent=None, opt_atom_ids=None):
+
+    if filename.endswith('.inp'):
+        with open(filename, 'w') as inp_file:
+            print('!', *keywords, file=inp_file)
+            if solvent:
+                print('%cpcm\n smd true\n SMDsolvent \"' + solvent + '\"\n end', file=inp_file)
+            print('%scf \nmaxiter 250 \nend', file=inp_file)
+            print('% maxcore 4000', file=inp_file)
+            if opt_atom_ids:
+                print('%geom Constraints', file=inp_file)
+                for i in opt_atom_ids:
+                    print('{ C', i, 'C }', file=inp_file)
+                print('end\n', 'invertConstraints true\n', 'end', sep='', file=inp_file)
+            print('*xyz', charge, mult, file=inp_file)
+            [print(*line, sep='\t', file=inp_file) for line in xyzs]
+            print('*', file=inp_file)
+
+    return 0
+
+
+def gen_orca_inp(inp_filename, xyzs, charge, mult, opt=True, opt_atom_ids=None, n_cores=1):
+    logger.info('Generating {}'.format(inp_filename))
+
+    if Config.opt_method == 'xtb':
+        keywords = ['LooseOpt', 'XTB']
+    elif Config.opt_method == 'min_dft':
+        keywords = ['LooseOpt', 'PBE', 'D3BJ', 'RI', 'def2-SVP', 'SlowConv']
+    elif Config.opt_method == 'min_dft_normal':
+        keywords = ['Opt', 'PBE', 'RI', 'D3BJ', 'def2-SVP', 'SlowConv']
+    elif Config.opt_method == 'dft':
+        keywords = ['Opt', 'PBE0', 'RIJCOSX', 'D3BJ', 'def2-SVP']
+    else:
+        logger.warning('No opt method set. Defaulting to PBE-D3BJ/def2-SVP/LooseOpt')
+        keywords = ['LooseOpt', 'PBE', 'D3BJ', 'RI', 'def2-SVP', 'SlowConv']
+
+    if not opt:
+        if 'LooseOpt' in keywords:
+            keywords.remove('LooseOpt')
+        if 'Opt' in keywords:
+            keywords.remove('Opt')
+
+    if 1 < n_cores <= 8:
+        keywords.append('PAL' + str(n_cores))
+    if n_cores > 8:                                 # ORCA calculations are capped at using 8 cores
+        logger.warning("Number of ORCA cores requested was {}. Reducing to 8".format(n_cores))
+        keywords.append('PAL8')
+
+    xys2inp(xyzs, keywords, inp_filename, charge, mult, solvent=None, opt_atom_ids=opt_atom_ids)
+
+    return 0
+
+
+def run_orca(inp_filename, out_filename):
+    logger.info('Running ORCA calculation from {}'.format(inp_filename))
+
+    orca_done, orca_terminated_normally = False, False
+
+    if os.path.exists(out_filename):
+        logger.info('Out file already exists. Checking to see whether it completed')
+        orca_done = True
+
+    if orca_done:
+        out_lines = [line for line in open(out_filename, 'r', encoding="utf-8")]
+        for line in reversed(out_lines):
+            if 'ORCA TERMINATED NORMALLY' in line:
+                orca_terminated_normally = True
+                if not Config.suppress_print:
+                    print("{:<30s}{:<50s}{:>10s}".format('Found ORCA run done for', inp_filename, 'Skipping'))
+                break
+
+    if not orca_terminated_normally:
+        with open(out_filename, 'w') as orca_out:
+            orca_run = Popen([Config.path_to_orca, inp_filename], stdout=orca_out)
+        orca_run.wait()
+
+    return [line for line in open(out_filename, 'r', encoding="utf-8")]
+
+
+def get_orca_xyzs_energy(out_lines):
+    logger.info('Getting xyzs and final energy from an ORCA output file')
+
+    opt_converged, geom_section = False, False
+    opt_xyzs = []
+    energy = 0.0
+
+    for line in out_lines:
+
+        if 'THE OPTIMIZATION HAS CONVERGED' in line:
+            opt_converged = True
+        if 'CARTESIAN COORDINATES' in line and opt_converged:
+            geom_section = True
+
+        if geom_section and len(line.split()) == 0:
+            geom_section = False
+
+        if geom_section and len(line.split()) == 4:
+            atom_label, x, y, z = line.split()
+            opt_xyzs.append([atom_label, float(x), float(y), float(z)])
+
+        if 'FINAL SINGLE POINT ENERGY' in line:
+            energy = float(line.split()[4])             # e.g. line = 'FINAL SINGLE POINT ENERGY     -4143.815610365798'
+
+    if len(opt_xyzs) == 0:
+        logger.warning('Could not find any xyzs in ORCA output')
+    if energy == 0:
+        logger.warning('Could not find energy in ORCA output')
+
+    return opt_xyzs, energy