MSchauperl / resppol

Compare 3483a2a ... +7 ... 598132b

Coverage Reach
resppol.py utilities.py __init__.py

No flags found

Use flags to group coverage reports by test type, project and/or folders.
Then setup custom commit statuses and notifications for each flag.

e.g., #unittest #integration

#production #enterprise

#frontend #backend

Learn more about Codecov Flags here.


@@ -0,0 +1,74 @@
Loading
1 +
from openeye import oechem
2 +
import sys
3 +
4 +
5 +
atomic_number={'H':1,
6 +
               'C':6,
7 +
               'N':7,
8 +
               'O':8}
9 +
10 +
def read_mulliken_from_g09(g09_output):
11 +
    """
12 +
    Searches a gaussian output file for the Mulliken charge block and return the element and charge of every atom
13 +
14 +
    :param g09_output: Gaussian output file
15 +
    :return: List of [element,charge]
16 +
    """
17 +
    f=open(g09_output)
18 +
    lines=f.readlines()
19 +
    f.close()
20 +
    element_charge_array=[]
21 +
    for linenum,line in enumerate(lines):
22 +
        if 'Mulliken charges' in line:
23 +
            while True:
24 +
                if 'Sum' in lines[linenum+2]:
25 +
                    break
26 +
                else:
27 +
                    num,element,charge=lines[linenum+2].split()
28 +
                    element_charge_array.append([element,charge])
29 +
                    linenum+=1
30 +
            return(element_charge_array)
31 +
32 +
33 +
34 +
def load_g09_charges_to_oemol(oemol, g09_output, mol2_output=None):
35 +
    """
36 +
    This functions takes a oemol and a link to a gaussian output file and matches the charges to the mol2 file.
37 +
    Optional a new mol2 file with updated charges can be created with the optional keyword mol2_output.
38 +
39 +
    :param oemol: openeye molecule object
40 +
    :param g09_output: gaussian output file
41 +
    :param mol2_output: mol2 output file with updated charges (optional)
42 +
43 +
    :return: True
44 +
    """
45 +
46 +
    #1.) Read in Mulliken charge from a gaussian output file.
47 +
    element_charges=read_mulliken_from_g09(g09_output)
48 +
    #2.) Check if structure in gaussian output file is the same as in the mol2 file
49 +
    # This comparison is just based on elements are therefore not ideal
50 +
    for i,atom in enumerate(oemol.GetAtoms()):
51 +
        if oemol.NumAtoms() != len(element_charges):
52 +
            raise LookupError('Molecule and gaussian output file do not match')
53 +
        if atom.GetAtomicNum == atomic_number[element_charges[i][0]]:
54 +
            raise LookupError('Molecule and gaussian output file do not match')
55 +
56 +
    #3.) Assign Charges to oemol
57 +
    for i,atom in enumerate(oemol.GetAtoms()):
58 +
        atom.SetPartialCharge(float(element_charges[i][1]))
59 +
60 +
    #4.) Write mol2 output if requested.
61 +
    ofs = oechem.oemolostream()
62 +
63 +
    #Check if I can open the output file
64 +
    if mol2_output is not None:
65 +
        if ofs.open(mol2_output):
66 +
            #If yes I write the mol2 file
67 +
            ofs.SetFormat(oechem.OEFormat_MOL2)
68 +
            oechem.OEWriteMolecule(ofs, oemol)
69 +
70 +
        else:
71 +
            # If not, I throw an exception
72 +
            oechem.OEThrow.Fatal("Unable to create {}".format(mol2_output))
73 +
74 +

@@ -17,24 +17,32 @@
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:
23 -
	print('Could not import openeye')
24 +
	log.warning('Could not import openeye')
24 25
try:
25 26
	import openforcefield.topology as openff
26 27
	from openforcefield.typing.engines.smirnoff import ForceField
27 28
except:
28 -
	print("Could not import openforcefield")
29 +
	log.warning("Could not import openforcefield")
29 30
from scipy.spatial import distance
30 31
import os
31 32
# Initialize units
32 33
from pint import UnitRegistry
34 +
import resppol.utilities as util
33 35
34 36
ureg = UnitRegistry()
35 37
Q_ = ureg.Quantity
36 38
ureg.define('bohr = 0.52917721067 * angstrom')
37 39
40 +
##### LOGGER
41 +
42 +
FORMAT = '%(asctime)s  - %(levelname)s - %(message)s'
43 +
log.basicConfig(format=FORMAT, level=log.INFO)
44 +
log.getLogger().addHandler(log.StreamHandler(sys.stdout))
45 +
38 46
# =============================================================================================
39 47
# GLOBAL PARAMETERS
40 48
# =============================================================================================
@@ -46,6 +54,7 @@
Loading
46 54
# =============================================================================================
47 55
# PRIVATE SUBROUTINES
48 56
# =============================================================================================
57 +
49 58
def find_eq_atoms(mol1):
50 59
    """
51 60
    Finds all equivalent atoms in a molecule.
@@ -55,7 +64,8 @@
Loading
55 64
    :parameter
56 65
    mol1: Openeye molecule object
57 66
58 -
    TODO Include rdkit support for this function
67 +
    TODO:
68 +
     Include rdkit support for this function
59 69
    """
60 70
61 71
    qmol = oechem.OEQMol()
@@ -64,6 +74,12 @@
Loading
64 74
    oechem.OEBuildMDLQueryExpressions(qmol, mol1)
65 75
    ss2 = oechem.OESubSearch(qmol)
66 76
    oechem.OEPrepareSearch(mol1, ss2)
77 +
    qmol = oechem.OEQMol()
78 +
    oechem.OEBuildMDLQueryExpressions(qmol, mol1)
79 +
    ss2 = oechem.OESubSearch(qmol)
80 +
    oechem.OEPrepareSearch(mol1, ss2)
81 +
    # it is not an error that i actually write this twice.
82 +
    # For phenol it seems to be necessary for an undefined reason.
67 83
68 84
    # Store the equivalent atoms
69 85
    eq_atoms = [[] for i in range(mol1.NumAtoms())]
@@ -85,6 +101,7 @@
Loading
85 101
            # Making sure we have not already covered this pair of atoms
86 102
            if [element, i] not in sorted_eq_atoms:
87 103
                sorted_eq_atoms.append([i, element])
104 +
88 105
    tmp = []
89 106
    for k, ele1 in enumerate(sorted_eq_atoms):
90 107
        exclude = 0
@@ -98,7 +115,6 @@
Loading
98 115
    return sorted_eq_atoms
99 116
100 117
101 -
102 118
# =============================================================================================
103 119
# TrainingSet
104 120
# =============================================================================================
@@ -120,38 +136,76 @@
Loading
120 136
        :param SCF:
121 137
        :param thole:
122 138
        """
139 +
        #Reading input and setting default values
140 +
        # ToDo add logging messages to the definition
123 141
        self.molecules = list()
142 +
        # Initialize matrix B with no values
124 143
        self.B = np.zeros(0)
144 +
        # Initialize matrix A with 2 dimensions and no values
125 145
        self.A = np.zeros((0, 0))
146 +
        # Initialize vector q with charge 0.0
147 +
        #ToDo: Check if this makes sense
126 148
        self.q = 0.0
149 +
        # Mode defines if we are optimizing charges, BCCs or Polarizabilites
150 +
        # Simulatenous optimization of charges or BCCs with Polarizabilities is possible
127 151
        self.mode = mode
152 +
        # Set SCF scaling parameter for all molecules
153 +
        # ToDo maybe i want to add a default value here
128 154
        self.scf_scaleparameters = scf_scaleparameters
155 +
        # Scaling parameters for charge charge interaction and charge dipole interaction
156 +
        # ToDo maybe i want to add a default value here
129 157
        self.scaleparameters = scaleparameters
158 +
        log.info('Charge Dipole interactions are scaled using the following parameters {}. '.format(self.scaleparameters))
159 +
        log.info('If non are specified the default is 0.0, 0.0, 0.8333. '
160 +
                 'This means that 1-2 and 1-3 interactions are neglected and 1-4 are scaled by a factor 0.8333. ')
161 +
        # Define if we use SCF or not default is the direct approach
130 162
        self._SCF = SCF
163 +
        # If thole scaling is applied
131 164
        self._thole = thole
165 +
        if self._SCF == True and self._thole == False:
166 +
            log.info('''Dipoles are calculated using the self-consistent field approach
167 +
                        The interactions are scaled using the following scaling factors: {}'''.format(self.scf_scaleparameters))
168 +
        if self._SCF == True and self._thole == True:
169 +
            log.info('''Dipoles are calculated using the self-consistent field approach
170 +
                    The interactions are scaled using a thole distance based scheme''')
171 +
        # ToDo Delete due to repetion. Check if always _mode is used before doing that thought.
132 172
        self._mode = mode
173 +
        log.info('Running in Mode: {}'.format(self._mode))
174 +
        # counts the number of optimization steps required
133 175
        self.step = 0
176 +
        # Counts the number of the charge /BCC matrix part
134 177
        self.number_of_lines_in_X = 0
178 +
        # Not sure
179 +
        # ToDo Check what this is good for
135 180
        self.X_BCC = 0.0
136 181
        self.Y_BCC = 0.0
182 +
        # Set the force-field file which is used for all molecules
137 183
        self._FF = FF
138 -
        # Label the atoms and bonds using a offxml file
139 -
        forcefield = ForceField(os.path.join(ROOT_DIR_PATH, FF))
184 +
        # This block describes how to label the atoms and bonds using a offxml file
140 185
186 +
        # Read in the forcefield using the ForceField class of the openforcefield toolkit
187 +
        forcefield = ForceField(os.path.join(ROOT_DIR_PATH, FF))
188 +
        log.info('The forcefield used for this run is {} '.format(os.path.join(ROOT_DIR_PATH, FF)))
189 +
        # Number of different polarization types
141 190
        self._nalpha = len(forcefield.get_parameter_handler('vdW').parameters)
191 +
192 +
        #Number of different BCCs definied in the FF
142 193
        self._nbccs = len(forcefield.get_parameter_handler('Bonds').parameters)
194 +
        # Assigns a number for every BCC
143 195
        self._bccs ={}
144 196
        for i,element in enumerate(forcefield.get_parameter_handler('Bonds').parameters):
145 197
            self._bccs[element.id] = i
146 198
147 -
199 +
        # The same for polarizabilities
148 200
        self._alpha = {}
149 201
        for i,element in enumerate(forcefield.get_parameter_handler('vdW').parameters):
150 202
            self._alpha[element.id] = i
151 203
152 204
        self.bccs_old = np.zeros(self._nbccs)
153 205
        self.alphas_old = np.zeros(self._nalpha)
154 206
207 +
    # This function is not done yet
208 +
    # Maybe replace with json format reader and writer
155 209
    def load_from_file(self,txtfile):
156 210
        """
157 211
        Allows to build a TrainingSet instance from an text file.
@@ -167,13 +221,21 @@
Loading
167 221
            # noinspection PyTypeChecker
168 222
            self.add_molecule(Molecule(mol2file))
169 223
170 -
    def add_molecule(self, datei):
224 +
    def add_molecule(self, datei,am1_file=None):
171 225
        """
172 -
        Adds a molecule.
226 +
        Adds a molecule to the TrainingSet object.
173 227
        :param datei: Mol2 file of a molecule
174 228
        :return:
175 229
        """
176 -
        self.molecules.append(Molecule(datei, position=self.number_of_lines_in_X, trainingset=self))
230 +
        #ToDo Check matrix composition
231 +
232 +
        # Adds the molecule
233 +
        self.molecules.append(Molecule(datei, position=self.number_of_lines_in_X, trainingset=self,id=len(self.molecules),am1_file=am1_file))
234 +
        log.info('Added molecule number {} from file {}.'.format(len(self.molecules)-1,datei))
235 +
        # Defines the position of the molecule in the overall matrix.
236 +
        #   A   0       CT          =   B
237 +
        #   0    A(new) C(new)T     =   B(new)
238 +
        #   C    C(new) P           =   B(Pol)
177 239
        self.number_of_lines_in_X += self.molecules[-1]._lines_in_X
178 240
179 241
@@ -187,17 +249,25 @@
Loading
187 249
        """
188 250
        Builds the matrix A of the underlying molecules and combines them.
189 251
190 -
        This function is only used for charge optimizatoin RESP
252 +
        This function is only used for charge optimization with RESP
191 253
        """
192 -
        for molecule in self.molecules:
254 +
        #
255 +
        for i,molecule in enumerate(self.molecules):
256 +
            # Building matrix A on the molecule level
193 257
            molecule.build_matrix_A()
258 +
            log.info('Build matrix A (only charge optimization) for molecule {}'.format(i))
259 +
            #Defines the interaction matrix 0 for charge of different molecules
194 260
            X12 = np.zeros((len(self.A), len(molecule.A)))
261 +
            # Combines the matrix A of the TrainingSet and the matrix A of the molecule
195 262
            self.A = np.concatenate(
196 263
                (np.concatenate((self.A, X12), axis=1), np.concatenate((X12.transpose(), molecule.A), axis=1)), axis=0)
197 264
198 265
    def build_vector_B(self):
199 -
        for molecule in self.molecules:
266 +
        for i,molecule in enumerate(self.molecules):
267 +
            # Building the vector B of the molecule.
200 268
            molecule.build_vector_B()
269 +
            log.info('Build vector B (only charge optimization) for molecule {}'.format(i))
270 +
            # Adds the vector B of the molecuel to the existing vector B
201 271
            self.B = np.concatenate((self.B, molecule.B))
202 272
203 273
    def build_matrix_X(self):
@@ -207,8 +277,9 @@
Loading
207 277
208 278
        This function is only used for charge optimizatoin RESP
209 279
        """
210 -
        for molecule in self.molecules:
280 +
        for i,molecule in enumerate(self.molecules):
211 281
            molecule.build_matrix_X()
282 +
            log.info('Build matrix X for molecule {}'.format(i))
212 283
            X12 = np.zeros((len(self.X), len(molecule.X)))
213 284
            self.X = np.concatenate(
214 285
                (np.concatenate((self.X, X12), axis=1), np.concatenate((X12.transpose(), molecule.X), axis=1)), axis=0)
@@ -222,43 +293,52 @@
Loading
222 293
            self.X[self.number_of_lines_in_X + i][atoms[1]] = self.X[atoms[1]][self.number_of_lines_in_X + i] = -1
223 294
224 295
296 +
    def build_vector_Y(self):
297 +
        self.Y = np.zeros(0)
298 +
        for i,molecule in enumerate(self.molecules):
299 +
            molecule.build_vector_Y()
300 +
            log.info('Build vector Y for molecule {}'.format(i))
301 +
            self.Y = np.concatenate((self.Y, molecule.Y))
302 +
303 +
        Y2 =  np.zeros(len(self.intramolecular_polarization_rst))
304 +
        self.Y = np.concatenate((self.Y, Y2))
305 +
    # def get_intramolecular_charge_rst()
306 +
225 307
    def build_matrix_X_BCC(self):
226 308
        """
227 309
        Builds the matrix X of the underlying molecules and combines them.
228 310
229 -
        This function is only used for charge optimizatoin RESP
311 +
        This function is only used for optimization of BCCs
230 312
        """
231 313
232 -
        for molecule in self.molecules:
314 +
        for i,molecule in enumerate(self.molecules):
233 315
            molecule.build_matrix_X_BCC()
316 +
            log.info('Build matrix X BCC for molecule {}'.format(i))
234 317
            self.X_BCC += molecule.X
235 318
236 319
237 320
    def build_vector_Y_BCC(self):
238 321
        """
239 322
        Builds the matrix X of the underlying molecules and combines them.
240 323
241 -
        This function is only used for charge optimizatoin RESP
324 +
        This function is only used for optimization of BCCs
242 325
        """
243 326
244 -
        for molecule in self.molecules:
327 +
        for i,molecule in enumerate(self.molecules):
245 328
            molecule.build_vector_Y_BCC()
329 +
            log.info('Build vector Y BCC for molecule {}'.format(i))
246 330
            self.Y_BCC += molecule.Y
247 331
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 332
254 -
        Y2 =  np.zeros(len(self.intramolecular_polarization_rst))
255 -
        self.Y = np.concatenate((self.Y, Y2))
256 -
    # def get_intramolecular_charge_rst()
257 333
258 334
    @property
259 335
    def intramolecular_polarization_rst(self):
260 -
336 +
        # Defines a list of atompairs which are restrained to be equal between different molecules
337 +
        # ToDo Check I do not include intramolecular polarization restraints here,just intermolecular ones
338 +
        # Doublecheck this behviour
339 +
        # Uses Lagrange formalism
261 340
        intramolecular_polarization_rst = []
341 +
        # First occurance of an atom with a given polarization type
262 342
        first_occurrence_of_parameter = {}
263 343
        for molecule in self.molecules:
264 344
            first_occurrence_of_parameter_in_molecule = {}
@@ -268,48 +348,71 @@
Loading
268 348
                elif atom._parameter_id not in first_occurrence_of_parameter_in_molecule.keys():
269 349
                    intramolecular_polarization_rst.append(
270 350
                        [first_occurrence_of_parameter[atom._parameter_id], molecule._position_in_A + atom._id + molecule._lines_in_A])
351 +
        log.info('Apply the following intramolecular polarization restaraints {}'.format(intramolecular_polarization_rst))
271 352
        return intramolecular_polarization_rst
272 353
273 354
    def optimize_charges_alpha(self,criteria = 10E-5):
274 355
        converged = False
275 -
        while converged == False:
356 +
        self.counter = 0
357 +
        while converged == False and self.counter<100:
358 +
        #while self.counter < 10:
359 +
            log.warning('Optimization Step {}'.format(self.counter))
276 360
            self.optimize_charges_alpha_step()
277 361
            converged = True
278 -
            for molecule in self.molecules:
362 +
            self.counter += 1
363 +
            for num,molecule in enumerate(self.molecules):
279 364
                molecule.step +=1
280 365
                if not all(abs(molecule.q - molecule.q_old) <criteria) or not all(abs(molecule.alpha - molecule.alpha_old) <criteria) :
281 366
                    converged = False
367 +
                    log.info('Optimization step {}: Molekule {}: Charges'.format(self.counter,num))
368 +
                    log.info(molecule.q)
369 +
                    log.info('Optimization step {}: Molekule {}: Dipoles'.format(self.counter,num))
370 +
                    log.info(molecule.alpha)
282 371
                molecule.q_old = molecule.q
283 -
                print(molecule.q)
284 -
                print(molecule.alpha)
285 -
                molecule.alpha_old =molecule.alpha
372 +
                molecule.alpha_old = molecule.alpha
373 +
                log.debug(molecule.q)
374 +
                log.debug(molecule.alpha)
375 +
            if self.counter ==99:
376 +
                log.warning('Optimization did not converge. Stopped after 100 cycles.')
286 377
287 378
288 379
289 380
    def optimize_charges_alpha_step(self):
290 381
        self.build_matrix_X()
291 382
        self.build_vector_Y()
383 +
        #log.info(self.X)
292 384
        self.q_alpha = np.linalg.solve(self.X, self.Y)
293 385
        # Update the charges of the molecules below
294 386
        q_alpha_tmp = self.q_alpha
295 387
        for molecule in self.molecules:
296 388
            molecule.q_alpha = q_alpha_tmp[:len(molecule.X)]
297 389
            q_alpha_tmp = q_alpha_tmp[len(molecule.X):]
298 390
            molecule.update_q_alpha()
391 +
        log.info('Updating charges and polarization parameters for all molecules')
299 392
300 -
    def optimize_bcc_alpha(self,criteria = 10E-5):
393 +
394 +
    def optimize_bcc_alpha(self,criteria = 10E-3):
301 395
        converged = False
302 -
        while converged == False:
396 +
        self.counter =0
397 +
        # Remove potential from AM1 charges
398 +
        #self.substract_am1_potential()
399 +
        while converged == False and self.counter <100:
400 +
            log.warning('Optimization Step {}'.format(self.counter))
303 401
            self.optimize_bcc_alpha_step()
304 -
            for molecule in self.molecules:
402 +
            self.counter+=1
403 +
            for num,molecule in enumerate(self.molecules):
305 404
                molecule.step +=1
306 -
                print(molecule.q)
307 -
                print(molecule.alpha)
405 +
            log.info('Optimization step {}: BCCs'.format(self.counter))
406 +
            log.info(self.bccs)
407 +
            log.info('Optimization step {}: Dipoles'.format(self.counter))
408 +
            log.info(self.alphas)
308 409
            if all(abs(self.bccs - self.bccs_old) < criteria) and all(abs(self.alphas - self.alphas_old)< criteria):
309 410
                converged = True
310 411
            else:
311 412
                self.alphas_old = self.alphas
312 413
                self.bccs_old = self.bccs
414 +
            if self.counter ==99:
415 +
                log.warning('Optimization did not converge. Stopped after 100 cycles.')
313 416
314 417
315 418
    def optimize_bcc_alpha_step(self):
@@ -327,6 +430,7 @@
Loading
327 430
        self.alphas = self.bcc_alpha[self._nbccs:]
328 431
        for molecule in self.molecules:
329 432
            molecule.update_bcc(self.bccs, self.alphas)
433 +
        log.info('Updating Bccs for all molecules')
330 434
331 435
332 436
@@ -343,6 +447,9 @@
Loading
343 447
            molecule.update_q()
344 448
345 449
450 +
    def substract_am1_potential(self):
451 +
        for molecule in self.molecules:
452 +
            molecule.substract_am1_potential()
346 453
# =============================================================================================
347 454
# Molecule
348 455
# =============================================================================================
@@ -351,7 +458,9 @@
Loading
351 458
class Molecule:
352 459
    """
353 460
        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
461 +
        an openeye molecule with the functionality of the OFF molecule.
462 +
463 +
        The OE part is only needed to determine
355 464
        the chemical equivalent atoms. When this feature is implemented in OFF toolkit the OE molecule is not
356 465
        necessary anymore
357 466
@@ -360,9 +469,10 @@
Loading
360 469
        :return:
361 470
        """
362 471
363 -
    def __init__(self, datei, position=0, trainingset=None):
472 +
    def __init__(self, datei, position=0, trainingset=None,id=None,am1_file=None):
364 473
365 474
        # Get Molecle name from filename
475 +
        self.id = id
366 476
        self.B = 0.0
367 477
        self.A = 0.0
368 478
        self.X = 0.0
@@ -396,6 +506,7 @@
Loading
396 506
        ifs = oechem.oemolistream(datei)
397 507
398 508
        # Read from IFS to molecule
509 +
        # Create OE molecule
399 510
        oechem.OEReadMol2File(ifs, self.oemol)
400 511
401 512
        # Check if molecule is a 3 dimensional object
@@ -421,6 +532,7 @@
Loading
421 532
        self.offtop = openff.Topology.from_molecules([self.offmol])
422 533
423 534
        # Label the atoms and bonds using a offxml file
535 +
        # Not sure if that is the most sensible thing to do. Maybe revisit this part at a later stage and define differently
424 536
        if trainingset is None:
425 537
            self._trainingset = TrainingSet()
426 538
@@ -431,20 +543,23 @@
Loading
431 543
432 544
        # set molecule charge
433 545
        # self._charge=openff.Molecule.total_charge
546 +
        # Change this back at a later point. To make charged molecules possible
434 547
        self._charge = 0
435 548
436 549
        # Initialize the bonds, atoms
437 550
        self._bonds = list()
438 551
        self._atoms = list()
439 552
440 553
        # Define all BCC bonds
554 +
        # ToDo write test for this function
441 555
        for i, properties in enumerate(molecule_parameter_list[0]['Bonds'].items()):
442 556
            atom_indices, parameter = properties
443 557
            self.add_bond(i, atom_indices, parameter.id)
444 558
445 559
        self._nbonds = len(self._bonds)
446 560
447 561
        # Define all atomtypes for polarization
562 +
        # ToDo write test for this function
448 563
        for i, properties in enumerate(molecule_parameter_list[0]['vdW'].items()):
449 564
            atom_index, parameter = properties
450 565
            self.add_atom(i, atom_index[0], parameter.id, atomic_number = AtomicNumbers[atom_index[0]])
@@ -455,12 +570,30 @@
Loading
455 570
        self.scale = np.ones((self._natoms, self._natoms))
456 571
        self.scale_scf = np.ones((self._natoms, self._natoms))
457 572
        self.scaling(scf_scaleparameters=self.scf_scaleparameters, scaleparameters=self.scaleparameters)
573 +
574 +
        # Defines the bond matrix T:
458 575
        self.create_BCCmatrix_T()
576 +
        # Defines the polarization type matrix R
459 577
        self.create_POLmatrix_R()
460 578
461 579
        # Initialize conformers
462 580
        self.conformers = list()
463 581
582 +
        # Load in charges from mol2 file:
583 +
        self.q_am1 = None
584 +
        if am1_file != None:
585 +
            if '.mol2' in am1_file:
586 +
                self.q_am1 = []
587 +
                self.am1mol = oechem.OEMol()
588 +
                # Open Input File Stream
589 +
                ifs = oechem.oemolistream(am1_file)
590 +
591 +
                # Read from IFS to molecule
592 +
                # Create OE molecule
593 +
                oechem.OEReadMol2File(ifs, self.am1mol)
594 +
                for atom in self.am1mol.GetAtoms():
595 +
                    self.q_am1.append(atom.GetPartialCharge())
596 +
464 597
        # Number of lines for matrix X
465 598
        if self._mode == 'q':
466 599
            self._lines_in_X = self._natoms + len(self.chemical_eq_atoms) + 1
@@ -472,6 +605,8 @@
Loading
472 605
            self.alpha_old = np.zeros(3*self._natoms)
473 606
        self._lines_in_A =  self._natoms + len(self.chemical_eq_atoms) + 1
474 607
        # Initiliaxe charges
608 +
        # Maybe I want to change that later to read charges from the mol2 file
609 +
        # or at least give the option to do so
475 610
        self.q_alpha = np.zeros(self._lines_in_X)
476 611
477 612
    def add_bond(self, index, atom_indices, parameter_id):
@@ -486,6 +621,7 @@
Loading
486 621
        atom_index1 = atom_indices[0]
487 622
        atom_index2 = atom_indices[1]
488 623
        self._bonds.append(Bond(index, atom_index1, atom_index2, parameter_id))
624 +
        log.info('Add bond with type {} between atom {} and atom {}'.format(parameter_id,atom_index1,atom_index2))
489 625
490 626
    def add_atom(self, index, atom_index, parameter_id, atomic_number=0):
491 627
        """
@@ -497,6 +633,7 @@
Loading
497 633
        :return:
498 634
        """
499 635
        self._atoms.append(Atom(index, atom_index, parameter_id, atomic_number= atomic_number))
636 +
        log.info('Add atom type {} on atom {}'.format(parameter_id, atom_index))
500 637
501 638
    def add_conformer_from_mol2(self, mol2file):
502 639
        """
@@ -517,6 +654,7 @@
Loading
517 654
            raise Exception('Conformer and Molecule does not match')
518 655
        else:
519 656
            self.conformers.append(Conformer(conf, molecule=self))
657 +
            log.info('Added conformation {} to molecule {}'.format(len(self.conformers)-1,self.id))
520 658
        # Checks if this conformer is from this moleule based on atom names
521 659
522 660
    @property
@@ -536,6 +674,7 @@
Loading
536 674
                first_occurrence[atom._parameter_id] = atom._id
537 675
            else:
538 676
                array_of_same_pol_atoms.append([first_occurrence[atom._parameter_id], atom._id ])
677 +
        log.info('The following atoms have the same polariztion type: {}'.format(array_of_same_pol_atoms))
539 678
        return (array_of_same_pol_atoms)
540 679
541 680
@@ -590,10 +729,11 @@
Loading
590 729
        :return:
591 730
        """
592 731
        for i in range(self._natoms):
593 -
            self.q_alpha[i] = np.dot(self.T[i], bccs)
732 +
            self.q_alpha[i] = self.q_am1[i] + np.dot(self.T[i], bccs)
594 733
        self.alpha = [np.dot(self.R[i], alphas) for i in range(len(self.R))]
595 734
        self.q_alpha[self._lines_in_A:self._lines_in_A + 3* len(self.alpha)] = np.concatenate(
596 735
            (np.concatenate((self.alpha, self.alpha)), self.alpha))
736 +
597 737
        self.q = self.q_alpha[:self._natoms]
598 738
        for conformer in self.conformers:
599 739
            conformer.q_alpha = self.q_alpha
@@ -610,6 +750,7 @@
Loading
610 750
            self.B += conformer.B
611 751
612 752
    def build_matrix_X(self):
753 +
        self.X=0.0
613 754
        for conformer in self.conformers:
614 755
            conformer.build_matrix_X()
615 756
            self.X += conformer.X
@@ -620,6 +761,7 @@
Loading
620 761
            self.X += conformer.X
621 762
622 763
    def build_vector_Y(self):
764 +
        self.Y=0.0
623 765
        for conformer in self.conformers:
624 766
            conformer.build_vector_Y()
625 767
            self.Y += conformer.Y
@@ -649,6 +791,15 @@
Loading
649 791
        for conformer in self.conformers:
650 792
            conformer.q_alpha = self.q_alpha
651 793
794 +
    def substract_am1_potential(self):
795 +
        if self.q_am1 is not None:
796 +
            for conformer in self.conformers:
797 +
                conformer.substract_am1_potential()
798 +
        else:
799 +
            self.q_am1 = np.zeros(self._natoms)
800 +
            log.warning('NO AM1 charges are definied')
801 +
802 +
652 803
    def scaling(self, scaleparameters=None, scf_scaleparameters=None):
653 804
        """
654 805
        Takes the bond information from a molecule instances and converts it to an scaling matrix.
@@ -772,14 +923,14 @@
Loading
772 923
773 924
        :param conf: openff molecule file
774 925
        """
926 +
        # Read in all atomic positions and convert them to Angstrom
775 927
        # noinspection PyProtectedMember
776 928
        self.atom_positions = Q_(np.array(conf.conformers[0]._value), 'angstrom')
777 929
        self.atom_positions_angstrom = self.atom_positions.to('angstrom').magnitude
778 930
        self.natoms = len(self.atom_positions.magnitude)
779 931
        self.baseESP = None
780 932
        self.polESPs = list()
781 933
        self._molecule = molecule
782 -
        self.q_alpha = self._molecule.q_alpha
783 934
        self._lines_in_A = self._molecule._lines_in_A
784 935
785 936
        # Initiliaze Electric field vectors
@@ -810,6 +961,15 @@
Loading
810 961
            if atom_is_present == 0:
811 962
                raise Exception("ESP grid does not belong to the conformer")
812 963
964 +
        #Substract the AM1 Potential if AM1 charges are specified
965 +
        if self._molecule.q_am1 != None:
966 +
            # calculated diatomic distances between atoms
967 +
            self.get_distances()
968 +
            # substract potential
969 +
            self.baseESP.substract_am1_potential()
970 +
            # delete distances()
971 +
            self.delete_distances()
972 +
813 973
    def add_polESP(self, *args, e_field=Q_([0.0, 0.0, 0.0], 'bohr')):
814 974
        """
815 975
        Adds the unpolarized molecule to this conformation.
@@ -824,6 +984,15 @@
Loading
824 984
        """
825 985
        self.polESPs.append(BCCPolESP(*args, conformer=self, e_field=e_field))
826 986
987 +
        #Substract the AM1 Potential if AM1 charges are specified
988 +
        if self._molecule.q_am1 != None:
989 +
            # calculated diatomic distances between atoms
990 +
            self.get_distances()
991 +
            # substract potential
992 +
            self.polESPs[-1].substract_am1_potential()
993 +
            # delete distances()
994 +
            self.delete_distances()
995 +
827 996
    # Build the matrix A for the charge optimization
828 997
    def build_matrix_A(self):
829 998
@@ -1211,30 +1380,40 @@
Loading
1211 1380
        for polESP in self.polESPs:
1212 1381
            polESP.get_electric_field()
1213 1382
1214 -
    def optimize_charges_alpha(self):
1215 -
        """
1216 -
        Builds the necessary matrix and vector and performs a charge and polarizabilities optimization for this 1 conformation.
1217 -
        :return:
1218 -
        """
1219 -
        self.build_matrix_X()
1220 -
        self.build_vector_Y()
1221 -
        self.q_alpha = np.linalg.solve(self.X, self.Y)
1222 -
1223 -
    def optimize_charges_alpha_bcc(self):
1224 -
        """
1225 -
        Builds the necessary matrix and vector and performs a charge and polarizabilities optimization for this 1 conformation.
1226 -
        :return:
1227 -
        """
1228 -
        self.build_matrix_X_BCC()
1229 -
        self.build_vector_Y_BCC()
1230 -
        # Check if bccs or alphas are not in the training set and set them to zero
1231 -
        for i,row in enumerate(self.X):
1232 -
            if all(row == 0.0):
1233 -
                self.X[i][i]=1
1234 -
1235 -
1236 -
        self.q_alpha = np.linalg.solve(self.X, self.Y)
1237 -
1383 +
    # def substract_am1_potential(self):
1384 +
    #     self.get_distances()
1385 +
    #     for ESPGRID in [self.baseESP] + self.polESPs:
1386 +
    #         ESPGRID.substract_am1_potential()
1387 +
    #     self.delete_distances()
1388 +
1389 +
    # def optimize_charges_alpha(self):
1390 +
    #     """
1391 +
    #     Builds the necessary matrix and vector and performs a charge and polarizabilities optimization for this 1 conformation.
1392 +
    #     :return:
1393 +
    #     """
1394 +
    #     self.build_matrix_X()
1395 +
    #     self.build_vector_Y()
1396 +
    #     self.q_alpha = np.linalg.solve(self.X, self.Y)
1397 +
1398 +
    # def optimize_charges_alpha_bcc(self):
1399 +
    #     """
1400 +
    #     Builds the necessary matrix and vector and performs a charge and polarizabilities optimization for this 1 conformation.
1401 +
    #     :return:
1402 +
    #     """
1403 +
    #     self.build_matrix_X_BCC()
1404 +
    #     self.build_vector_Y_BCC()
1405 +
    #     # Check if bccs or alphas are not in the training set and set them to zero
1406 +
    #     for i,row in enumerate(self.X):
1407 +
    #         if all(row == 0.0):
1408 +
    #             self.X[i][i]=1
1409 +
    #
1410 +
    #
1411 +
    #     self.q_alpha = np.linalg.solve(self.X, self.Y)
1412 +
1413 +
1414 +
    # i constantly have to delete the distatnce matrixes as storing them is too expensive in terms of memory
1415 +
    # recomputing them every optimization cyclces has shown to be faster as long as i have not an extensive
1416 +
    # amount of RAM available.
1238 1417
1239 1418
    def get_distances(self):
1240 1419
        self.get_grid_coord()
@@ -1344,10 +1523,14 @@
Loading
1344 1523
1345 1524
    def define_grid(self, *args, **kwargs):
1346 1525
        for ele in args:
1526 +
            # Check if the ending of the file is gesp
1347 1527
            if 'gesp' in ele:
1348 1528
                self.gridtype = 'gesp'
1529 +
            # Check if one of the files end with espf
1530 +
            # This is the same as a psi4 file but respyte uses this name
1349 1531
            if 'espf' in ele:
1350 1532
                self.gridtype = 'respyte'
1533 +
            # psi4 type of grid
1351 1534
            if 'grid.dat' in ele:
1352 1535
                self.gridtype = 'psi4'
1353 1536
@@ -1364,13 +1547,19 @@
Loading
1364 1547
                if 'ATOMIC' in line and 'COORDINATES' in line:
1365 1548
                    self.natoms = int(line.strip('\n').split()[-1])
1366 1549
                    for j in range(self.natoms):
1550 +
                        # Convert the very unique D in gaussian number to E (the usual)
1367 1551
                        entry = lines[i + 1 + j].replace('D', 'E').split()
1368 1552
                        self.atoms.append(entry[0])
1553 +
                        # Units are in Bohr in Gaussian
1554 +
                        # pint should allow us to make the conversion easy
1369 1555
                        self.atom_positions.append(Q_([float(entry[k]) for k in range(1, 4, 1)], 'bohr'))
1370 1556
                if 'GRID' in line:
1557 +
                    # number of grid points
1371 1558
                    self.ngrid = int(line.strip('\n').split()[-1])
1559 +
                    # Stores the grid
1372 1560
                    grid = lines[i + 1:i + 1 + self.ngrid]
1373 1561
                    break
1562 +
            # ToDo combine this with above. Current implementation is correct but dose not make a lot of sense
1374 1563
            # noinspection PyUnboundLocalVariable
1375 1564
            for i, line in enumerate(grid):
1376 1565
                grid[i] = [float(ele) for ele in line.replace('D', 'E').split()]
@@ -1382,12 +1571,14 @@
Loading
1382 1571
            f = open(args[0], 'r')
1383 1572
            lines = f.readlines()
1384 1573
            f.close
1574 +
            # I only need every second line because the other line is the electric field at point i
1385 1575
            ndata = int(len(lines) / 2) if len(lines) % 2 == 0 else int((len(lines) - 1) / 2)
1386 1576
            grid = np.zeros((ndata, 4))
1387 1577
            for i in range(ndata):
1388 1578
                grid[i] = [float(ele) for ele in lines[2 * i].split()]
1389 1579
            self.positions = Q_(np.array(grid)[:, 0:3], 'angstrom')
1390 1580
            self.esp_values = Q_(np.array(grid)[:, 3], 'elementary_charge / bohr')
1581 +
        # Assuming that we only have the ESP in the esp.dat file
1391 1582
        elif self.gridtype == 'psi4':
1392 1583
            for ele in args:
1393 1584
                if "grid.dat" in ele:
@@ -1399,44 +1590,45 @@
Loading
1399 1590
            self.positions = Q_(np.loadtxt(gridfile), 'angstrom')
1400 1591
            self.esp_values = Q_(np.loadtxt(espfile), 'elementary_charge / bohr')
1401 1592
1593 +
1402 1594
    def get_electric_field(self, ):
1403 1595
        """
1404 1596
        Calculates the electric field at every atomic positions.
1405 1597
        :return:
1406 1598
        """
1407 -
        alpha = self._conformer.q_alpha[self._conformer._lines_in_A:self._conformer._lines_in_A + 3*self._conformer.natoms]
1599 +
        alpha = self._molecule.q_alpha[self._conformer._lines_in_A:self._conformer._lines_in_A + 3*self._conformer.natoms]
1408 1600
        alpha[np.where(alpha == 0.0)] += 10E-10
1409 1601
1410 1602
        # Load permanent charges for BCC method
1411 1603
        # For all other methods this is set to 0.0 STILL HAVE to implment
1412 -
        try:
1413 -
            self._conformer.q_am1
1414 -
        except Exception:
1604 +
        if self._molecule.q_am1 is None:
1415 1605
            log.warning('I do not have AM1-type charges')
1416 -
            self._conformer.q_am1 = np.zeros(self._conformer.natoms)
1606 +
            self._molecule.q_am1 = np.zeros(self._conformer.natoms)
1607 +
        #else:
1608 +
        #    log.info('Substracting AM1 charge potential from the ESP values')
1609 +
        # ToDo Add Code to substract am1 charge potential
1417 1610
1418 1611
        for j in range(self._conformer.natoms):
1419 1612
            self.e_field_at_atom[0][j] = np.dot(
1420 -
                np.multiply(self._conformer.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1613 +
                np.multiply(self._molecule.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1421 1614
                self._conformer.diatomic_dist_x[j]) + np.dot(
1422 -
                np.multiply(self._conformer.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1615 +
                np.multiply(self._molecule.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1423 1616
                self._conformer.diatomic_dist_x[j]) + self._ext_e_field[0].to(
1424 1617
                'elementary_charge / angstrom / angstrom').magnitude
1425 1618
            self.e_field_at_atom[1][j] = np.dot(
1426 -
                np.multiply(self._conformer.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1619 +
                np.multiply(self._molecule.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1427 1620
                self._conformer.diatomic_dist_y[j]) + np.dot(
1428 -
                np.multiply(self._conformer.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1621 +
                np.multiply(self._molecule.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1429 1622
                self._conformer.diatomic_dist_y[j]) + self._ext_e_field[1].to(
1430 1623
                'elementary_charge / angstrom / angstrom').magnitude
1431 1624
            self.e_field_at_atom[2][j] = np.dot(
1432 -
                np.multiply(self._conformer.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1625 +
                np.multiply(self._molecule.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1433 1626
                self._conformer.diatomic_dist_z[j]) + np.dot(
1434 -
                np.multiply(self._conformer.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1627 +
                np.multiply(self._molecule.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1435 1628
                self._conformer.diatomic_dist_z[j]) + self._ext_e_field[2].to(
1436 1629
                'elementary_charge / angstrom / angstrom').magnitude
1437 1630
1438 1631
        self.e = self.e_field_at_atom.flatten()
1439 -
        #print(self.e)
1440 1632
        if self._conformer._molecule._SCF and self._conformer._molecule.step > 0:
1441 1633
            if not hasattr(self, 'Bdip') or self._conformer._molecule._thole:
1442 1634
                if self._conformer._molecule._thole:
@@ -1476,7 +1668,6 @@
Loading
1476 1668
                        thole_ft = self._conformer._molecule.scale
1477 1669
                        thole_fe = self._conformer._molecule.scale
1478 1670
                    else:
1479 -
                        print('Using different set of scaling for SCF interactions')
1480 1671
                        log.info('Using different set of scaling for SCF interactions')
1481 1672
                Bdip11 = np.add(np.multiply(thole_fe, self._conformer.diatomic_dist_3), np.multiply(thole_ft,
1482 1673
                                                                                                    -3 * np.multiply(
@@ -1522,13 +1713,13 @@
Loading
1522 1713
                Bdip[j][j] = 1. / alpha[j]
1523 1714
            dipole_scf = np.linalg.solve(Bdip, self.e)
1524 1715
            self.e = np.divide(dipole_scf, alpha)
1525 -
            print(self.e)
1716 +
            log.debug(self.e)
1526 1717
        self.e_field_at_atom[0] = 1.0 * self.e[:self._conformer.natoms]
1527 1718
        self.e_field_at_atom[1] = 1.0 * self.e[self._conformer.natoms:2 * self._conformer.natoms]
1528 1719
        self.e_field_at_atom[2] = 1.0 * self.e[2 * self._conformer.natoms:3 * self._conformer.natoms]
1529 1720
1530 1721
    # Analysation stuff
1531 -
    def calc_esp_q_alpha(self, q_alpha):
1722 +
    def calc_esp_q_alpha(self, q_alpha, mode='q'):
1532 1723
        """
1533 1724
        Calculates the ESP for a given set of point charges and polarizabilities.
1534 1725
@@ -1540,7 +1731,7 @@
Loading
1540 1731
        """
1541 1732
        self.q_pot = np.zeros(len(self.esp_values))
1542 1733
        for i in range(len(self.esp_values)):
1543 -
            if self._conformer._molecule._mode == 'q':
1734 +
            if mode == 'q':
1544 1735
                self.q_pot[i] = np.dot(q_alpha[:self._conformer.natoms], np.transpose(self._conformer.dist)[i])
1545 1736
            else:
1546 1737
                dipole = q_alpha[self._conformer._lines_in_A:]
@@ -1556,12 +1747,16 @@
Loading
1556 1747
                    self.q_pot[i] = np.dot(q_alpha[:self.natoms], np.transpose(self._conformer.dist)[i]) + e_dip
1557 1748
                except:
1558 1749
                    self.q_pot[i] = np.dot(q_alpha[:self.natoms], np.transpose(self._conformer.dist)[i])
1559 -
                    # print('DEBUG no electric field ')
1560 1750
            # if self.mode=='d':
1561 1751
            #    self.dipole=qd
1562 1752
            #    e_dip=np.dot(np.multiply(self.dipole[:self.natoms],self.e_x),np.transpose(self.dist_x)[i])+np.dot(np.multiply(self.dipole[self.natoms:2*self.natoms],self.e_y),np.transpose(self.dist_y)[i])+np.dot(np.multiply(self.dipole[2*self.natoms:3*self.natoms],self.e_z),np.transpose(self.dist_z)[i])
1563 1753
            #    self.q_pot[i]=e_dip
1564 1754
1755 +
    def substract_am1_potential(self):
1756 +
        self.calc_esp_q_alpha(self._molecule.q_am1, mode='q')
1757 +
        self.esp_values = Q_(np.subtract(self.esp_values.to('elementary_charge / angstrom').magnitude, self.q_pot),
1758 +
                         'elementary_charge / angstrom')
1759 +
1565 1760
    def sub_esp_q_alpha(self, q_alpha):
1566 1761
        """
1567 1762
        Subtracts the ESP create by a set of point charges and polarizabilities.
@@ -1605,7 +1800,7 @@
Loading
1605 1800
        for i in range(self._conformer.natoms):
1606 1801
            f.write(
1607 1802
                ' {} {} {} {} {}\n'.format(atoms[i][0], atoms[i][1][0], atoms[i][1][1], atoms[i][1][2],
1608 -
                                           self._conformer._molecule.q_alpha[i]))
1803 +
                                           self._molecule.q_alpha[i]))
1609 1804
        f.write(' ESP VALUES AND GRID POINT COORDINATES. #POINTS =   {}\n'.format(len(self.esp_values)))
1610 1805
        for i in range(len(self.esp_values)):
1611 1806
            try:
@@ -1618,6 +1813,7 @@
Loading
1618 1813
# BCCUnpolESP
1619 1814
# =============================================================================================
1620 1815
1816 +
# Inheritace from the ESPGRID class
1621 1817
class BCCUnpolESP(ESPGRID):
1622 1818
    """
1623 1819
@@ -1626,27 +1822,41 @@
Loading
1626 1822
    def __init__(self, *args, conformer=None):
1627 1823
        # Decide if we have a Gaussian grid or a psi4 grid
1628 1824
        self.gridtype = None
1825 +
1826 +
        # Initilize values
1629 1827
        self.natoms = -1
1828 +
1829 +
        #Initialize array to store the atom position and the element from the ESP file
1830 +
        #ToDo Check if those are the same as in the Conformer file
1630 1831
        self.atoms = []
1631 1832
        self.atom_positions = []
1833 +
1834 +
        self._conformer = conformer
1835 +
        self._molecule = conformer._molecule
1836 +
        # Checks what grid type psi4 or gaussian was used
1632 1837
        self.define_grid(*args)
1838 +
1633 1839
        self.esp_values = None
1634 1840
        self.positions = None
1635 -
        self._conformer = conformer
1841 +
1636 1842
1637 1843
        # External e-field is 0 in all directions
1638 1844
        vector = Q_([0, 0, 0], 'elementary_charge / bohr / bohr')
1639 1845
        self.set_ext_e_field(vector)
1640 1846
1847 +
        # Load the grid to self.positions and self.esp_values
1641 1848
        self.load_grid(*args)
1642 1849
1850 +
        # As the electric field at the atom positions is different for every external polarization
1851 +
        # the electric field for each atom has to be stored on the ESPGrid level
1643 1852
        self.e_field_at_atom = np.zeros((3, self._conformer.natoms))
1853 +
        # No calculation of the e-fields at this stage only initializing
1644 1854
1645 1855
1646 1856
# =============================================================================================
1647 1857
# BCCPolESP
1648 1858
# =============================================================================================
1649 -
1859 +
# Inheritace from the ESPGRID class
1650 1860
class BCCPolESP(ESPGRID):
1651 1861
    """
1652 1862
@@ -1662,6 +1872,7 @@
Loading
1662 1872
        self.esp_values = None
1663 1873
        self.positions = None
1664 1874
        self._conformer = conformer
1875 +
        self._molecule = conformer._molecule
1665 1876
1666 1877
        # Set external e-field
1667 1878
        self.set_ext_e_field(e_field)
@@ -1672,7 +1883,23 @@
Loading
1672 1883
1673 1884
1674 1885
if __name__ == '__main__':
1675 -
    pass
1886 +
    ifs = oechem.oemolistream('/home/mschauperl/kirk/charge_method/medium_set/molecule2/conf0/mp2_0.mol2')
1887 +
    oemol=oechem.OEMol()
1888 +
    oechem.OEReadMol2File(ifs, oemol)
1889 +
    ifs2 = oechem.oemolistream('/home/mschauperl/kirk/charge_method/medium_set/molecule3/conf0/mp2_0.mol2')
1890 +
    oemol2=oechem.OEMol()
1891 +
    oechem.OEReadMol2File(ifs2, oemol2)
1892 +
    ifs3 = oechem.oemolistream('/home/mschauperl/programs/resppol/resppol/data/test_data/butanol_0.mol2')
1893 +
    oemol3=oechem.OEMol()
1894 +
    oechem.OEReadMol2File(ifs3, oemol3)
1895 +
    log.debug("This is a test")
1896 +
    print(find_eq_atoms(oemol))
1897 +
    print(find_eq_atoms(oemol))
1898 +
    print(find_eq_atoms(oemol2))
1899 +
    print(find_eq_atoms(oemol2))
1900 +
    print(find_eq_atoms(oemol3))
1901 +
    print(find_eq_atoms(oemol3))
1902 +
1676 1903
    """
1677 1904
    datei = os.path.join(ROOT_DIR_PATH, 'resppol/tmp/phenol/conf0/mp2_0.mol2')
1678 1905
    test = TrainingSet(scf_scaleparameters=[0.0, 0.0, 0.5])

Learn more Showing 2 files with coverage changes found.

Changes in resppol/resppol.py
-2
+2
Loading file...
New file resppol/utilities.py
New
Loading file...
Files Coverage
resppol -1.70% 91.61%
Project Totals (3 files) 91.61%
Loading