Showing 1 of 5 files from the diff.

@@ -17,6 +17,7 @@
Loading
17 17
18 18
import logging as log
19 19
import numpy as np
20 +
import sys
20 21
try:
21 22
	from openeye import oechem
22 23
except:
@@ -35,6 +36,12 @@
Loading
35 36
Q_ = ureg.Quantity
36 37
ureg.define('bohr = 0.52917721067 * angstrom')
37 38
39 +
##### LOGGER
40 +
41 +
FORMAT = '%(asctime)s  - %(levelname)s - %(message)s'
42 +
log.basicConfig(format=FORMAT, level=log.DEBUG)
43 +
log.getLogger().addHandler(log.StreamHandler(sys.stdout))
44 +
38 45
# =============================================================================================
39 46
# GLOBAL PARAMETERS
40 47
# =============================================================================================
@@ -55,7 +62,8 @@
Loading
55 62
    :parameter
56 63
    mol1: Openeye molecule object
57 64
58 -
    TODO Include rdkit support for this function
65 +
    TODO:
66 +
     Include rdkit support for this function
59 67
    """
60 68
61 69
    qmol = oechem.OEQMol()
@@ -64,6 +72,12 @@
Loading
64 72
    oechem.OEBuildMDLQueryExpressions(qmol, mol1)
65 73
    ss2 = oechem.OESubSearch(qmol)
66 74
    oechem.OEPrepareSearch(mol1, ss2)
75 +
    qmol = oechem.OEQMol()
76 +
    oechem.OEBuildMDLQueryExpressions(qmol, mol1)
77 +
    ss2 = oechem.OESubSearch(qmol)
78 +
    oechem.OEPrepareSearch(mol1, ss2)
79 +
    # it is not an error that i actually write this twice.
80 +
    # For phenol it seems to be necessary for an undefined reason.
67 81
68 82
    # Store the equivalent atoms
69 83
    eq_atoms = [[] for i in range(mol1.NumAtoms())]
@@ -120,31 +134,59 @@
Loading
120 134
        :param SCF:
121 135
        :param thole:
122 136
        """
137 +
        #Reading input and setting default values
138 +
        # ToDo add logging messages to the definition
123 139
        self.molecules = list()
140 +
        # Initialize matrix B with no values
124 141
        self.B = np.zeros(0)
142 +
        # Initialize matrix A with 2 dimensions and no values
125 143
        self.A = np.zeros((0, 0))
144 +
        # Initialize vector q with charge 0.0
145 +
        #ToDo: Check if this makes sense
126 146
        self.q = 0.0
147 +
        # Mode defines if we are optimizing charges, BCCs or Polarizabilites
148 +
        # Simulatenous optimization of charges or BCCs with Polarizabilities is possible
127 149
        self.mode = mode
150 +
        # Set SCF scaling parameter for all molecules
151 +
        # ToDo maybe i want to add a default value here
128 152
        self.scf_scaleparameters = scf_scaleparameters
153 +
        # Scaling parameters for charge charge interaction and charge dipole interaction
154 +
        # ToDo maybe i want to add a default value here
129 155
        self.scaleparameters = scaleparameters
156 +
        # Define if we use SCF or not default is the direct approach
130 157
        self._SCF = SCF
158 +
        # If thole scaling is applied
159 +
        # ToDo Check if in adding to scaling or instead
131 160
        self._thole = thole
161 +
        # ToDo Delete due to repetion. Check if always _mode is used before doing that thought.
132 162
        self._mode = mode
163 +
        # counts the number of optimization steps required
133 164
        self.step = 0
165 +
        # Counts the number of the charge /BCC matrix part
134 166
        self.number_of_lines_in_X = 0
167 +
        # Not sure
168 +
        # ToDo Check what this is good for
135 169
        self.X_BCC = 0.0
136 170
        self.Y_BCC = 0.0
171 +
        # Set the force-field file which is used for all molecules
137 172
        self._FF = FF
138 -
        # Label the atoms and bonds using a offxml file
173 +
174 +
        # This block describes how to label the atoms and bonds using a offxml file
175 +
176 +
        # Read in the forcefield using the ForceField class of the openforcefield toolkit
139 177
        forcefield = ForceField(os.path.join(ROOT_DIR_PATH, FF))
140 178
179 +
        # Number of different polarization types
141 180
        self._nalpha = len(forcefield.get_parameter_handler('vdW').parameters)
181 +
182 +
        #Number of different BCCs definied in the FF
142 183
        self._nbccs = len(forcefield.get_parameter_handler('Bonds').parameters)
184 +
        # Assigns a number for every BCC
143 185
        self._bccs ={}
144 186
        for i,element in enumerate(forcefield.get_parameter_handler('Bonds').parameters):
145 187
            self._bccs[element.id] = i
146 188
147 -
189 +
        # The same for polarizabilities
148 190
        self._alpha = {}
149 191
        for i,element in enumerate(forcefield.get_parameter_handler('vdW').parameters):
150 192
            self._alpha[element.id] = i
@@ -152,6 +194,8 @@
Loading
152 194
        self.bccs_old = np.zeros(self._nbccs)
153 195
        self.alphas_old = np.zeros(self._nalpha)
154 196
197 +
    # This function is not done yet
198 +
    # Maybe replace with json format reader and writer
155 199
    def load_from_file(self,txtfile):
156 200
        """
157 201
        Allows to build a TrainingSet instance from an text file.
@@ -169,11 +213,19 @@
Loading
169 213
170 214
    def add_molecule(self, datei):
171 215
        """
172 -
        Adds a molecule.
216 +
        Adds a molecule to the TrainingSet object.
173 217
        :param datei: Mol2 file of a molecule
174 218
        :return:
175 219
        """
220 +
        #ToDo Check matrix composition
221 +
222 +
        # Adds the molecule
176 223
        self.molecules.append(Molecule(datei, position=self.number_of_lines_in_X, trainingset=self))
224 +
225 +
        # Defines the position of the molecule in the overall matrix.
226 +
        #   A   0       CT          =   B
227 +
        #   0    A(new) C(new)T     =   B(new)
228 +
        #   C    C(new) P           =   B(Pol)
177 229
        self.number_of_lines_in_X += self.molecules[-1]._lines_in_X
178 230
179 231
@@ -187,17 +239,23 @@
Loading
187 239
        """
188 240
        Builds the matrix A of the underlying molecules and combines them.
189 241
190 -
        This function is only used for charge optimizatoin RESP
242 +
        This function is only used for charge optimization with RESP
191 243
        """
244 +
        #
192 245
        for molecule in self.molecules:
246 +
            # Building matrix A on the molecule level
193 247
            molecule.build_matrix_A()
248 +
            #Defines the interaction matrix 0 for charge of different molecules
194 249
            X12 = np.zeros((len(self.A), len(molecule.A)))
250 +
            # Combines the matrix A of the TrainingSet and the matrix A of the molecule
195 251
            self.A = np.concatenate(
196 252
                (np.concatenate((self.A, X12), axis=1), np.concatenate((X12.transpose(), molecule.A), axis=1)), axis=0)
197 253
198 254
    def build_vector_B(self):
199 255
        for molecule in self.molecules:
256 +
            # Building the vector B of the molecule.
200 257
            molecule.build_vector_B()
258 +
            # Adds the vector B of the molecuel to the existing vector B
201 259
            self.B = np.concatenate((self.B, molecule.B))
202 260
203 261
    def build_matrix_X(self):
@@ -222,11 +280,21 @@
Loading
222 280
            self.X[self.number_of_lines_in_X + i][atoms[1]] = self.X[atoms[1]][self.number_of_lines_in_X + i] = -1
223 281
224 282
283 +
    def build_vector_Y(self):
284 +
        self.Y = np.zeros(0)
285 +
        for molecule in self.molecules:
286 +
            molecule.build_vector_Y()
287 +
            self.Y = np.concatenate((self.Y, molecule.Y))
288 +
289 +
        Y2 =  np.zeros(len(self.intramolecular_polarization_rst))
290 +
        self.Y = np.concatenate((self.Y, Y2))
291 +
    # def get_intramolecular_charge_rst()
292 +
225 293
    def build_matrix_X_BCC(self):
226 294
        """
227 295
        Builds the matrix X of the underlying molecules and combines them.
228 296
229 -
        This function is only used for charge optimizatoin RESP
297 +
        This function is only used for optimization of BCCs
230 298
        """
231 299
232 300
        for molecule in self.molecules:
@@ -238,27 +306,23 @@
Loading
238 306
        """
239 307
        Builds the matrix X of the underlying molecules and combines them.
240 308
241 -
        This function is only used for charge optimizatoin RESP
309 +
        This function is only used for optimization of BCCs
242 310
        """
243 311
244 312
        for molecule in self.molecules:
245 313
            molecule.build_vector_Y_BCC()
246 314
            self.Y_BCC += molecule.Y
247 315
248 -
    def build_vector_Y(self):
249 -
        self.Y = np.zeros(0)
250 -
        for molecule in self.molecules:
251 -
            molecule.build_vector_Y()
252 -
            self.Y = np.concatenate((self.Y, molecule.Y))
253 316
254 -
        Y2 =  np.zeros(len(self.intramolecular_polarization_rst))
255 -
        self.Y = np.concatenate((self.Y, Y2))
256 -
    # def get_intramolecular_charge_rst()
257 317
258 318
    @property
259 319
    def intramolecular_polarization_rst(self):
260 -
320 +
        # Defines a list of atompairs which are restrained to be equal between different molecules
321 +
        # ToDo Check I do not include intramolecular polarization restraints here,just intermolecular ones
322 +
        # Doublecheck this behviour
323 +
        # Uses Lagrange formalism
261 324
        intramolecular_polarization_rst = []
325 +
        # First occurance of an atom with a given polarization type
262 326
        first_occurrence_of_parameter = {}
263 327
        for molecule in self.molecules:
264 328
            first_occurrence_of_parameter_in_molecule = {}
@@ -351,7 +415,9 @@
Loading
351 415
class Molecule:
352 416
    """
353 417
        This class loads in a mol2 file and expose all relevant information. It combines the functionality of
354 -
        an openeye molecule with the functionality of the OFF molecule. The OE part is only needed to determine
418 +
        an openeye molecule with the functionality of the OFF molecule.
419 +
420 +
        The OE part is only needed to determine
355 421
        the chemical equivalent atoms. When this feature is implemented in OFF toolkit the OE molecule is not
356 422
        necessary anymore
357 423
@@ -396,6 +462,7 @@
Loading
396 462
        ifs = oechem.oemolistream(datei)
397 463
398 464
        # Read from IFS to molecule
465 +
        # Create OE molecule
399 466
        oechem.OEReadMol2File(ifs, self.oemol)
400 467
401 468
        # Check if molecule is a 3 dimensional object
@@ -421,6 +488,7 @@
Loading
421 488
        self.offtop = openff.Topology.from_molecules([self.offmol])
422 489
423 490
        # Label the atoms and bonds using a offxml file
491 +
        # Not sure if that is the most sensible thing to do. Maybe revisit this part at a later stage and define differently
424 492
        if trainingset is None:
425 493
            self._trainingset = TrainingSet()
426 494
@@ -431,6 +499,7 @@
Loading
431 499
432 500
        # set molecule charge
433 501
        # self._charge=openff.Molecule.total_charge
502 +
        # Change this back at a later point. To make charged molecules possible
434 503
        self._charge = 0
435 504
436 505
        # Initialize the bonds, atoms
@@ -438,6 +507,7 @@
Loading
438 507
        self._atoms = list()
439 508
440 509
        # Define all BCC bonds
510 +
        # ToDo write test for this function
441 511
        for i, properties in enumerate(molecule_parameter_list[0]['Bonds'].items()):
442 512
            atom_indices, parameter = properties
443 513
            self.add_bond(i, atom_indices, parameter.id)
@@ -445,6 +515,7 @@
Loading
445 515
        self._nbonds = len(self._bonds)
446 516
447 517
        # Define all atomtypes for polarization
518 +
        # ToDo write test for this function
448 519
        for i, properties in enumerate(molecule_parameter_list[0]['vdW'].items()):
449 520
            atom_index, parameter = properties
450 521
            self.add_atom(i, atom_index[0], parameter.id, atomic_number = AtomicNumbers[atom_index[0]])
@@ -455,7 +526,10 @@
Loading
455 526
        self.scale = np.ones((self._natoms, self._natoms))
456 527
        self.scale_scf = np.ones((self._natoms, self._natoms))
457 528
        self.scaling(scf_scaleparameters=self.scf_scaleparameters, scaleparameters=self.scaleparameters)
529 +
530 +
        # Defines the bond matrix T:
458 531
        self.create_BCCmatrix_T()
532 +
        # Defines the polarization type matrix R
459 533
        self.create_POLmatrix_R()
460 534
461 535
        # Initialize conformers
@@ -472,6 +546,8 @@
Loading
472 546
            self.alpha_old = np.zeros(3*self._natoms)
473 547
        self._lines_in_A =  self._natoms + len(self.chemical_eq_atoms) + 1
474 548
        # Initiliaxe charges
549 +
        # Maybe I want to change that later to read charges from the mol2 file
550 +
        # or at least give the option to do so
475 551
        self.q_alpha = np.zeros(self._lines_in_X)
476 552
477 553
    def add_bond(self, index, atom_indices, parameter_id):
@@ -772,6 +848,7 @@
Loading
772 848
773 849
        :param conf: openff molecule file
774 850
        """
851 +
        # Read in all atomic positions and convert them to Angstrom
775 852
        # noinspection PyProtectedMember
776 853
        self.atom_positions = Q_(np.array(conf.conformers[0]._value), 'angstrom')
777 854
        self.atom_positions_angstrom = self.atom_positions.to('angstrom').magnitude
@@ -1236,6 +1313,10 @@
Loading
1236 1313
        self.q_alpha = np.linalg.solve(self.X, self.Y)
1237 1314
1238 1315
1316 +
    # i constantly have to delete the distatnce matrixes as storing them is too expensive in terms of memory
1317 +
    # recomputing them every optimization cyclces has shown to be faster as long as i have not an extensive
1318 +
    # amount of RAM available.
1319 +
1239 1320
    def get_distances(self):
1240 1321
        self.get_grid_coord()
1241 1322
@@ -1344,10 +1425,14 @@
Loading
1344 1425
1345 1426
    def define_grid(self, *args, **kwargs):
1346 1427
        for ele in args:
1428 +
            # Check if the ending of the file is gesp
1347 1429
            if 'gesp' in ele:
1348 1430
                self.gridtype = 'gesp'
1431 +
            # Check if one of the files end with espf
1432 +
            # This is the same as a psi4 file but respyte uses this name
1349 1433
            if 'espf' in ele:
1350 1434
                self.gridtype = 'respyte'
1435 +
            # psi4 type of grid
1351 1436
            if 'grid.dat' in ele:
1352 1437
                self.gridtype = 'psi4'
1353 1438
@@ -1364,13 +1449,19 @@
Loading
1364 1449
                if 'ATOMIC' in line and 'COORDINATES' in line:
1365 1450
                    self.natoms = int(line.strip('\n').split()[-1])
1366 1451
                    for j in range(self.natoms):
1452 +
                        # Convert the very unique D in gaussian number to E (the usual)
1367 1453
                        entry = lines[i + 1 + j].replace('D', 'E').split()
1368 1454
                        self.atoms.append(entry[0])
1455 +
                        # Units are in Bohr in Gaussian
1456 +
                        # pint should allow us to make the conversion easy
1369 1457
                        self.atom_positions.append(Q_([float(entry[k]) for k in range(1, 4, 1)], 'bohr'))
1370 1458
                if 'GRID' in line:
1459 +
                    # number of grid points
1371 1460
                    self.ngrid = int(line.strip('\n').split()[-1])
1461 +
                    # Stores the grid
1372 1462
                    grid = lines[i + 1:i + 1 + self.ngrid]
1373 1463
                    break
1464 +
            # ToDo combine this with above. Current implementation is correct but dose not make a lot of sense
1374 1465
            # noinspection PyUnboundLocalVariable
1375 1466
            for i, line in enumerate(grid):
1376 1467
                grid[i] = [float(ele) for ele in line.replace('D', 'E').split()]
@@ -1382,12 +1473,14 @@
Loading
1382 1473
            f = open(args[0], 'r')
1383 1474
            lines = f.readlines()
1384 1475
            f.close
1476 +
            # I only need every second line because the other line is the electric field at point i
1385 1477
            ndata = int(len(lines) / 2) if len(lines) % 2 == 0 else int((len(lines) - 1) / 2)
1386 1478
            grid = np.zeros((ndata, 4))
1387 1479
            for i in range(ndata):
1388 1480
                grid[i] = [float(ele) for ele in lines[2 * i].split()]
1389 1481
            self.positions = Q_(np.array(grid)[:, 0:3], 'angstrom')
1390 1482
            self.esp_values = Q_(np.array(grid)[:, 3], 'elementary_charge / bohr')
1483 +
        # Assuming that we only have the ESP in the esp.dat file
1391 1484
        elif self.gridtype == 'psi4':
1392 1485
            for ele in args:
1393 1486
                if "grid.dat" in ele:
@@ -1618,6 +1711,7 @@
Loading
1618 1711
# BCCUnpolESP
1619 1712
# =============================================================================================
1620 1713
1714 +
# Inheritace from the ESPGRID class
1621 1715
class BCCUnpolESP(ESPGRID):
1622 1716
    """
1623 1717
@@ -1626,27 +1720,39 @@
Loading
1626 1720
    def __init__(self, *args, conformer=None):
1627 1721
        # Decide if we have a Gaussian grid or a psi4 grid
1628 1722
        self.gridtype = None
1723 +
1724 +
        # Initilize values
1629 1725
        self.natoms = -1
1726 +
1727 +
        #Initialize array to store the atom position and the element from the ESP file
1728 +
        #ToDo Check if those are the same as in the Conformer file
1630 1729
        self.atoms = []
1631 1730
        self.atom_positions = []
1731 +
1732 +
        # Checks what grid type psi4 or gaussian was used
1632 1733
        self.define_grid(*args)
1734 +
1633 1735
        self.esp_values = None
1634 1736
        self.positions = None
1737 +
1635 1738
        self._conformer = conformer
1636 1739
1637 1740
        # External e-field is 0 in all directions
1638 1741
        vector = Q_([0, 0, 0], 'elementary_charge / bohr / bohr')
1639 1742
        self.set_ext_e_field(vector)
1640 1743
1744 +
        # Load the grid to self.positions and self.esp_values
1641 1745
        self.load_grid(*args)
1642 1746
1747 +
        # As the electric field at the atom positions is different for every external polarization
1748 +
        # the electric field for each atom has to be stored on the ESPGrid level
1643 1749
        self.e_field_at_atom = np.zeros((3, self._conformer.natoms))
1644 -
1750 +
        # No calculation of the e-fields at this stage only initializing
1645 1751
1646 1752
# =============================================================================================
1647 1753
# BCCPolESP
1648 1754
# =============================================================================================
1649 -
1755 +
# Inheritace from the ESPGRID class
1650 1756
class BCCPolESP(ESPGRID):
1651 1757
    """
1652 1758
@@ -1672,7 +1778,23 @@
Loading
1672 1778
1673 1779
1674 1780
if __name__ == '__main__':
1675 -
    pass
1781 +
    ifs = oechem.oemolistream('/home/mschauperl/kirk/charge_method/medium_set/molecule2/conf0/mp2_0.mol2')
1782 +
    oemol=oechem.OEMol()
1783 +
    oechem.OEReadMol2File(ifs, oemol)
1784 +
    ifs2 = oechem.oemolistream('/home/mschauperl/kirk/charge_method/medium_set/molecule3/conf0/mp2_0.mol2')
1785 +
    oemol2=oechem.OEMol()
1786 +
    oechem.OEReadMol2File(ifs2, oemol2)
1787 +
    ifs3 = oechem.oemolistream('/home/mschauperl/programs/resppol/resppol/data/test_data/butanol_0.mol2')
1788 +
    oemol3=oechem.OEMol()
1789 +
    oechem.OEReadMol2File(ifs3, oemol3)
1790 +
    log.debug("This is a test")
1791 +
    print(find_eq_atoms(oemol))
1792 +
    print(find_eq_atoms(oemol))
1793 +
    print(find_eq_atoms(oemol2))
1794 +
    print(find_eq_atoms(oemol2))
1795 +
    print(find_eq_atoms(oemol3))
1796 +
    print(find_eq_atoms(oemol3))
1797 +
1676 1798
    """
1677 1799
    datei = os.path.join(ROOT_DIR_PATH, 'resppol/tmp/phenol/conf0/mp2_0.mol2')
1678 1800
    test = TrainingSet(scf_scaleparameters=[0.0, 0.0, 0.5])
Files Coverage
resppol 91.66%
Project Totals (2 files) 91.66%
444.1
TRAVIS_OS_NAME=osx
444.3
TRAVIS_PYTHON_VERSION=3.6
TRAVIS_OS_NAME=linux
444.4
TRAVIS_PYTHON_VERSION=3.7
TRAVIS_OS_NAME=linux
444.2
TRAVIS_OS_NAME=osx
1
# Codecov configuration to make it a bit less noisy
2
coverage:
3
  status:
4
    patch: false
5
    project:
6
      default:
7
        threshold: 50%
8
comment:
9
  layout: "header"
10
  require_changes: false
11
  branches: null
12
  behavior: default
13
  flags: null
14
  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