@@ -21,12 +21,12 @@
Loading
21 21
try:
22 22
	from openeye import oechem
23 23
except:
24 -
	print('Could not import openeye')
24 +
	log.warning('Could not import openeye')
25 25
try:
26 26
	import openforcefield.topology as openff
27 27
	from openforcefield.typing.engines.smirnoff import ForceField
28 28
except:
29 -
	print("Could not import openforcefield")
29 +
	log.warning("Could not import openforcefield")
30 30
from scipy.spatial import distance
31 31
import os
32 32
# Initialize units
@@ -39,7 +39,7 @@
Loading
39 39
##### LOGGER
40 40
41 41
FORMAT = '%(asctime)s  - %(levelname)s - %(message)s'
42 -
log.basicConfig(format=FORMAT, level=log.DEBUG)
42 +
log.basicConfig(format=FORMAT, level=log.INFO)
43 43
log.getLogger().addHandler(log.StreamHandler(sys.stdout))
44 44
45 45
# =============================================================================================
@@ -99,6 +99,7 @@
Loading
99 99
            # Making sure we have not already covered this pair of atoms
100 100
            if [element, i] not in sorted_eq_atoms:
101 101
                sorted_eq_atoms.append([i, element])
102 +
102 103
    tmp = []
103 104
    for k, ele1 in enumerate(sorted_eq_atoms):
104 105
        exclude = 0
@@ -153,13 +154,22 @@
Loading
153 154
        # Scaling parameters for charge charge interaction and charge dipole interaction
154 155
        # ToDo maybe i want to add a default value here
155 156
        self.scaleparameters = scaleparameters
157 +
        log.info('Charge Dipole interactions are scaled using the following parameters {}. '.format(self.scaleparameters))
158 +
        log.info('If non are specified the default is 0.0, 0.0, 0.8333. '
159 +
                 'This means that 1-2 and 1-3 interactions are neglected and 1-4 are scaled by a factor 0.8333. ')
156 160
        # Define if we use SCF or not default is the direct approach
157 161
        self._SCF = SCF
158 162
        # If thole scaling is applied
159 -
        # ToDo Check if in adding to scaling or instead
160 163
        self._thole = thole
164 +
        if self._SCF == True and self._thole == False:
165 +
            log.info('''Dipoles are calculated using the self-consistent field approach
166 +
                        The interactions are scaled using the following scaling factors: {}'''.format(self.scf_scaleparameters))
167 +
        if self._SCF == True and self._thole == True:
168 +
            log.info('''Dipoles are calculated using the self-consistent field approach
169 +
                    The interactions are scaled using a thole distance based scheme''')
161 170
        # ToDo Delete due to repetion. Check if always _mode is used before doing that thought.
162 171
        self._mode = mode
172 +
        log.info('Running in Mode: {}'.format(self._mode))
163 173
        # counts the number of optimization steps required
164 174
        self.step = 0
165 175
        # Counts the number of the charge /BCC matrix part
@@ -170,12 +180,11 @@
Loading
170 180
        self.Y_BCC = 0.0
171 181
        # Set the force-field file which is used for all molecules
172 182
        self._FF = FF
173 -
174 183
        # This block describes how to label the atoms and bonds using a offxml file
175 184
176 185
        # Read in the forcefield using the ForceField class of the openforcefield toolkit
177 186
        forcefield = ForceField(os.path.join(ROOT_DIR_PATH, FF))
178 -
187 +
        log.info('The forcefield used for this run is {} '.format(os.path.join(ROOT_DIR_PATH, FF)))
179 188
        # Number of different polarization types
180 189
        self._nalpha = len(forcefield.get_parameter_handler('vdW').parameters)
181 190
@@ -220,8 +229,8 @@
Loading
220 229
        #ToDo Check matrix composition
221 230
222 231
        # Adds the molecule
223 -
        self.molecules.append(Molecule(datei, position=self.number_of_lines_in_X, trainingset=self))
224 -
232 +
        self.molecules.append(Molecule(datei, position=self.number_of_lines_in_X, trainingset=self,id=len(self.molecules)))
233 +
        log.info('Added molecule number {} from file {}.'.format(len(self.molecules)-1,datei))
225 234
        # Defines the position of the molecule in the overall matrix.
226 235
        #   A   0       CT          =   B
227 236
        #   0    A(new) C(new)T     =   B(new)
@@ -242,9 +251,10 @@
Loading
242 251
        This function is only used for charge optimization with RESP
243 252
        """
244 253
        #
245 -
        for molecule in self.molecules:
254 +
        for i,molecule in enumerate(self.molecules):
246 255
            # Building matrix A on the molecule level
247 256
            molecule.build_matrix_A()
257 +
            log.info('Build matrix A (only charge optimization) for molecule {}'.format(i))
248 258
            #Defines the interaction matrix 0 for charge of different molecules
249 259
            X12 = np.zeros((len(self.A), len(molecule.A)))
250 260
            # Combines the matrix A of the TrainingSet and the matrix A of the molecule
@@ -252,9 +262,10 @@
Loading
252 262
                (np.concatenate((self.A, X12), axis=1), np.concatenate((X12.transpose(), molecule.A), axis=1)), axis=0)
253 263
254 264
    def build_vector_B(self):
255 -
        for molecule in self.molecules:
265 +
        for i,molecule in enumerate(self.molecules):
256 266
            # Building the vector B of the molecule.
257 267
            molecule.build_vector_B()
268 +
            log.info('Build vector B (only charge optimization) for molecule {}'.format(i))
258 269
            # Adds the vector B of the molecuel to the existing vector B
259 270
            self.B = np.concatenate((self.B, molecule.B))
260 271
@@ -265,8 +276,9 @@
Loading
265 276
266 277
        This function is only used for charge optimizatoin RESP
267 278
        """
268 -
        for molecule in self.molecules:
279 +
        for i,molecule in enumerate(self.molecules):
269 280
            molecule.build_matrix_X()
281 +
            log.info('Build matrix X for molecule {}'.format(i))
270 282
            X12 = np.zeros((len(self.X), len(molecule.X)))
271 283
            self.X = np.concatenate(
272 284
                (np.concatenate((self.X, X12), axis=1), np.concatenate((X12.transpose(), molecule.X), axis=1)), axis=0)
@@ -282,8 +294,9 @@
Loading
282 294
283 295
    def build_vector_Y(self):
284 296
        self.Y = np.zeros(0)
285 -
        for molecule in self.molecules:
297 +
        for i,molecule in enumerate(self.molecules):
286 298
            molecule.build_vector_Y()
299 +
            log.info('Build vector Y for molecule {}'.format(i))
287 300
            self.Y = np.concatenate((self.Y, molecule.Y))
288 301
289 302
        Y2 =  np.zeros(len(self.intramolecular_polarization_rst))
@@ -297,8 +310,9 @@
Loading
297 310
        This function is only used for optimization of BCCs
298 311
        """
299 312
300 -
        for molecule in self.molecules:
313 +
        for i,molecule in enumerate(self.molecules):
301 314
            molecule.build_matrix_X_BCC()
315 +
            log.info('Build matrix X BCC for molecule {}'.format(i))
302 316
            self.X_BCC += molecule.X
303 317
304 318
@@ -309,8 +323,9 @@
Loading
309 323
        This function is only used for optimization of BCCs
310 324
        """
311 325
312 -
        for molecule in self.molecules:
326 +
        for i,molecule in enumerate(self.molecules):
313 327
            molecule.build_vector_Y_BCC()
328 +
            log.info('Build vector Y BCC for molecule {}'.format(i))
314 329
            self.Y_BCC += molecule.Y
315 330
316 331
@@ -332,27 +347,39 @@
Loading
332 347
                elif atom._parameter_id not in first_occurrence_of_parameter_in_molecule.keys():
333 348
                    intramolecular_polarization_rst.append(
334 349
                        [first_occurrence_of_parameter[atom._parameter_id], molecule._position_in_A + atom._id + molecule._lines_in_A])
350 +
        log.info('Apply the following intramolecular polarization restaraints {}'.format(intramolecular_polarization_rst))
335 351
        return intramolecular_polarization_rst
336 352
337 353
    def optimize_charges_alpha(self,criteria = 10E-5):
338 354
        converged = False
339 -
        while converged == False:
355 +
        self.counter = 0
356 +
        while converged == False and self.counter<100:
357 +
        #while self.counter < 10:
358 +
            log.warning('Optimization Step {}'.format(self.counter))
340 359
            self.optimize_charges_alpha_step()
341 360
            converged = True
342 -
            for molecule in self.molecules:
361 +
            self.counter += 1
362 +
            for num,molecule in enumerate(self.molecules):
343 363
                molecule.step +=1
344 364
                if not all(abs(molecule.q - molecule.q_old) <criteria) or not all(abs(molecule.alpha - molecule.alpha_old) <criteria) :
345 365
                    converged = False
366 +
                    log.info('Optimization step {}: Molekule {}: Charges'.format(self.counter,num))
367 +
                    log.info(molecule.q)
368 +
                    log.info('Optimization step {}: Molekule {}: Dipoles'.format(self.counter,num))
369 +
                    log.info(molecule.alpha)
346 370
                molecule.q_old = molecule.q
347 -
                print(molecule.q)
348 -
                print(molecule.alpha)
349 -
                molecule.alpha_old =molecule.alpha
371 +
                molecule.alpha_old = molecule.alpha
372 +
                log.debug(molecule.q)
373 +
                log.debug(molecule.alpha)
374 +
            if self.counter ==99:
375 +
                log.warning('Optimization did not converge. Stopped after 100 cycles.')
350 376
351 377
352 378
353 379
    def optimize_charges_alpha_step(self):
354 380
        self.build_matrix_X()
355 381
        self.build_vector_Y()
382 +
        #log.info(self.X)
356 383
        self.q_alpha = np.linalg.solve(self.X, self.Y)
357 384
        # Update the charges of the molecules below
358 385
        q_alpha_tmp = self.q_alpha
@@ -360,20 +387,29 @@
Loading
360 387
            molecule.q_alpha = q_alpha_tmp[:len(molecule.X)]
361 388
            q_alpha_tmp = q_alpha_tmp[len(molecule.X):]
362 389
            molecule.update_q_alpha()
390 +
        log.info('Updating charges and polarization parameters for all molecules')
391 +
363 392
364 393
    def optimize_bcc_alpha(self,criteria = 10E-5):
365 394
        converged = False
366 -
        while converged == False:
395 +
        self.counter =0
396 +
        while converged == False and self.counter <100:
397 +
            log.warning('Optimization Step {}'.format(self.counter))
367 398
            self.optimize_bcc_alpha_step()
368 -
            for molecule in self.molecules:
399 +
            self.counter+=1
400 +
            for num,molecule in enumerate(self.molecules):
369 401
                molecule.step +=1
370 -
                print(molecule.q)
371 -
                print(molecule.alpha)
402 +
            log.info('Optimization step {}: BCCs'.format(self.counter))
403 +
            log.info(self.bccs)
404 +
            log.info('Optimization step {}: Dipoles'.format(self.counter))
405 +
            log.info(self.alphas)
372 406
            if all(abs(self.bccs - self.bccs_old) < criteria) and all(abs(self.alphas - self.alphas_old)< criteria):
373 407
                converged = True
374 408
            else:
375 409
                self.alphas_old = self.alphas
376 410
                self.bccs_old = self.bccs
411 +
            if self.counter ==99:
412 +
                log.warning('Optimization did not converge. Stopped after 100 cycles.')
377 413
378 414
379 415
    def optimize_bcc_alpha_step(self):
@@ -391,6 +427,7 @@
Loading
391 427
        self.alphas = self.bcc_alpha[self._nbccs:]
392 428
        for molecule in self.molecules:
393 429
            molecule.update_bcc(self.bccs, self.alphas)
430 +
        log.info('Updating Bccs for all molecules')
394 431
395 432
396 433
@@ -426,9 +463,10 @@
Loading
426 463
        :return:
427 464
        """
428 465
429 -
    def __init__(self, datei, position=0, trainingset=None):
466 +
    def __init__(self, datei, position=0, trainingset=None,id=None):
430 467
431 468
        # Get Molecle name from filename
469 +
        self.id = id
432 470
        self.B = 0.0
433 471
        self.A = 0.0
434 472
        self.X = 0.0
@@ -562,6 +600,7 @@
Loading
562 600
        atom_index1 = atom_indices[0]
563 601
        atom_index2 = atom_indices[1]
564 602
        self._bonds.append(Bond(index, atom_index1, atom_index2, parameter_id))
603 +
        log.info('Add bond with type {} between atom {} and atom {}'.format(parameter_id,atom_index1,atom_index2))
565 604
566 605
    def add_atom(self, index, atom_index, parameter_id, atomic_number=0):
567 606
        """
@@ -573,6 +612,7 @@
Loading
573 612
        :return:
574 613
        """
575 614
        self._atoms.append(Atom(index, atom_index, parameter_id, atomic_number= atomic_number))
615 +
        log.info('Add atom type {} on atom {}'.format(parameter_id, atom_index))
576 616
577 617
    def add_conformer_from_mol2(self, mol2file):
578 618
        """
@@ -593,6 +633,7 @@
Loading
593 633
            raise Exception('Conformer and Molecule does not match')
594 634
        else:
595 635
            self.conformers.append(Conformer(conf, molecule=self))
636 +
            log.info('Added conformation {} to molecule {}'.format(len(self.conformers)-1,self.id))
596 637
        # Checks if this conformer is from this moleule based on atom names
597 638
598 639
    @property
@@ -612,6 +653,7 @@
Loading
612 653
                first_occurrence[atom._parameter_id] = atom._id
613 654
            else:
614 655
                array_of_same_pol_atoms.append([first_occurrence[atom._parameter_id], atom._id ])
656 +
        log.info('The following atoms have the same polariztion type: {}'.format(array_of_same_pol_atoms))
615 657
        return (array_of_same_pol_atoms)
616 658
617 659
@@ -686,6 +728,7 @@
Loading
686 728
            self.B += conformer.B
687 729
688 730
    def build_matrix_X(self):
731 +
        self.X=0.0
689 732
        for conformer in self.conformers:
690 733
            conformer.build_matrix_X()
691 734
            self.X += conformer.X
@@ -696,6 +739,7 @@
Loading
696 739
            self.X += conformer.X
697 740
698 741
    def build_vector_Y(self):
742 +
        self.Y=0.0
699 743
        for conformer in self.conformers:
700 744
            conformer.build_vector_Y()
701 745
            self.Y += conformer.Y
@@ -1529,7 +1573,6 @@
Loading
1529 1573
                'elementary_charge / angstrom / angstrom').magnitude
1530 1574
1531 1575
        self.e = self.e_field_at_atom.flatten()
1532 -
        #print(self.e)
1533 1576
        if self._conformer._molecule._SCF and self._conformer._molecule.step > 0:
1534 1577
            if not hasattr(self, 'Bdip') or self._conformer._molecule._thole:
1535 1578
                if self._conformer._molecule._thole:
@@ -1569,7 +1612,6 @@
Loading
1569 1612
                        thole_ft = self._conformer._molecule.scale
1570 1613
                        thole_fe = self._conformer._molecule.scale
1571 1614
                    else:
1572 -
                        print('Using different set of scaling for SCF interactions')
1573 1615
                        log.info('Using different set of scaling for SCF interactions')
1574 1616
                Bdip11 = np.add(np.multiply(thole_fe, self._conformer.diatomic_dist_3), np.multiply(thole_ft,
1575 1617
                                                                                                    -3 * np.multiply(
@@ -1615,7 +1657,7 @@
Loading
1615 1657
                Bdip[j][j] = 1. / alpha[j]
1616 1658
            dipole_scf = np.linalg.solve(Bdip, self.e)
1617 1659
            self.e = np.divide(dipole_scf, alpha)
1618 -
            print(self.e)
1660 +
            log.debug(self.e)
1619 1661
        self.e_field_at_atom[0] = 1.0 * self.e[:self._conformer.natoms]
1620 1662
        self.e_field_at_atom[1] = 1.0 * self.e[self._conformer.natoms:2 * self._conformer.natoms]
1621 1663
        self.e_field_at_atom[2] = 1.0 * self.e[2 * self._conformer.natoms:3 * self._conformer.natoms]
@@ -1649,7 +1691,6 @@
Loading
1649 1691
                    self.q_pot[i] = np.dot(q_alpha[:self.natoms], np.transpose(self._conformer.dist)[i]) + e_dip
1650 1692
                except:
1651 1693
                    self.q_pot[i] = np.dot(q_alpha[:self.natoms], np.transpose(self._conformer.dist)[i])
1652 -
                    # print('DEBUG no electric field ')
1653 1694
            # if self.mode=='d':
1654 1695
            #    self.dipole=qd
1655 1696
            #    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])
Files Coverage
resppol 91.80%
Project Totals (2 files) 91.80%
452.4
TRAVIS_PYTHON_VERSION=3.7
TRAVIS_OS_NAME=linux
452.3
TRAVIS_PYTHON_VERSION=3.6
TRAVIS_OS_NAME=linux
452.1
TRAVIS_OS_NAME=osx
452.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