MSchauperl / resppol

Compare f9a1533 ... +4 ... 9c01119

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 +

@@ -31,6 +31,7 @@
Loading
31 31
import os
32 32
# Initialize units
33 33
from pint import UnitRegistry
34 +
import resppol.utilities as util
34 35
35 36
ureg = UnitRegistry()
36 37
Q_ = ureg.Quantity
@@ -53,6 +54,7 @@
Loading
53 54
# =============================================================================================
54 55
# PRIVATE SUBROUTINES
55 56
# =============================================================================================
57 +
56 58
def find_eq_atoms(mol1):
57 59
    """
58 60
    Finds all equivalent atoms in a molecule.
@@ -113,7 +115,6 @@
Loading
113 115
    return sorted_eq_atoms
114 116
115 117
116 -
117 118
# =============================================================================================
118 119
# TrainingSet
119 120
# =============================================================================================
@@ -220,7 +221,7 @@
Loading
220 221
            # noinspection PyTypeChecker
221 222
            self.add_molecule(Molecule(mol2file))
222 223
223 -
    def add_molecule(self, datei):
224 +
    def add_molecule(self, datei,am1_mol2_file=None):
224 225
        """
225 226
        Adds a molecule to the TrainingSet object.
226 227
        :param datei: Mol2 file of a molecule
@@ -229,7 +230,7 @@
Loading
229 230
        #ToDo Check matrix composition
230 231
231 232
        # Adds the molecule
232 -
        self.molecules.append(Molecule(datei, position=self.number_of_lines_in_X, trainingset=self,id=len(self.molecules)))
233 +
        self.molecules.append(Molecule(datei, position=self.number_of_lines_in_X, trainingset=self,id=len(self.molecules),am1_mol2_file=am1_mol2_file))
233 234
        log.info('Added molecule number {} from file {}.'.format(len(self.molecules)-1,datei))
234 235
        # Defines the position of the molecule in the overall matrix.
235 236
        #   A   0       CT          =   B
@@ -390,9 +391,11 @@
Loading
390 391
        log.info('Updating charges and polarization parameters for all molecules')
391 392
392 393
393 -
    def optimize_bcc_alpha(self,criteria = 10E-5):
394 +
    def optimize_bcc_alpha(self,criteria = 10E-3):
394 395
        converged = False
395 396
        self.counter =0
397 +
        # Remove potential from AM1 charges
398 +
        self.substract_am1_potential()
396 399
        while converged == False and self.counter <100:
397 400
            log.warning('Optimization Step {}'.format(self.counter))
398 401
            self.optimize_bcc_alpha_step()
@@ -444,6 +447,9 @@
Loading
444 447
            molecule.update_q()
445 448
446 449
450 +
    def substract_am1_potential(self):
451 +
        for molecule in self.molecules:
452 +
            molecule.substract_am1_potential()
447 453
# =============================================================================================
448 454
# Molecule
449 455
# =============================================================================================
@@ -463,7 +469,7 @@
Loading
463 469
        :return:
464 470
        """
465 471
466 -
    def __init__(self, datei, position=0, trainingset=None,id=None):
472 +
    def __init__(self, datei, position=0, trainingset=None,id=None,am1_mol2_file=None):
467 473
468 474
        # Get Molecle name from filename
469 475
        self.id = id
@@ -573,6 +579,20 @@
Loading
573 579
        # Initialize conformers
574 580
        self.conformers = list()
575 581
582 +
        # Load in charges from mol2 file:
583 +
        self.q_am1 = None
584 +
        if am1_mol2_file != None:
585 +
            self.q_am1 = []
586 +
            self.am1mol = oechem.OEMol()
587 +
            # Open Input File Stream
588 +
            ifs = oechem.oemolistream(am1_mol2_file)
589 +
590 +
            # Read from IFS to molecule
591 +
            # Create OE molecule
592 +
            oechem.OEReadMol2File(ifs, self.am1mol)
593 +
            for atom in self.am1mol.GetAtoms():
594 +
                self.q_am1.append(atom.GetPartialCharge())
595 +
576 596
        # Number of lines for matrix X
577 597
        if self._mode == 'q':
578 598
            self._lines_in_X = self._natoms + len(self.chemical_eq_atoms) + 1
@@ -708,7 +728,7 @@
Loading
708 728
        :return:
709 729
        """
710 730
        for i in range(self._natoms):
711 -
            self.q_alpha[i] = np.dot(self.T[i], bccs)
731 +
            self.q_alpha[i] = self.q_am1[i] + np.dot(self.T[i], bccs)
712 732
        self.alpha = [np.dot(self.R[i], alphas) for i in range(len(self.R))]
713 733
        self.q_alpha[self._lines_in_A:self._lines_in_A + 3* len(self.alpha)] = np.concatenate(
714 734
            (np.concatenate((self.alpha, self.alpha)), self.alpha))
@@ -769,6 +789,15 @@
Loading
769 789
        for conformer in self.conformers:
770 790
            conformer.q_alpha = self.q_alpha
771 791
792 +
    def substract_am1_potential(self):
793 +
        if self.q_am1 is not None:
794 +
            for conformer in self.conformers:
795 +
                conformer.substract_am1_potential()
796 +
        else:
797 +
            self.q_am1 = np.zeros(self._natoms)
798 +
            log.warning('NO AM1 charges are definied')
799 +
800 +
772 801
    def scaling(self, scaleparameters=None, scf_scaleparameters=None):
773 802
        """
774 803
        Takes the bond information from a molecule instances and converts it to an scaling matrix.
@@ -902,6 +931,7 @@
Loading
902 931
        self._molecule = molecule
903 932
        self.q_alpha = self._molecule.q_alpha
904 933
        self._lines_in_A = self._molecule._lines_in_A
934 +
        self.q_am1 = self._molecule.q_am1
905 935
906 936
        # Initiliaze Electric field vectors
907 937
        self.e_field_at_atom = np.zeros((3, self.natoms))
@@ -1332,6 +1362,12 @@
Loading
1332 1362
        for polESP in self.polESPs:
1333 1363
            polESP.get_electric_field()
1334 1364
1365 +
    def substract_am1_potential(self):
1366 +
        self.get_distances()
1367 +
        for ESPGRID in [self.baseESP] + self.polESPs:
1368 +
            ESPGRID.substract_am1_potential()
1369 +
        self.delete_distances()
1370 +
1335 1371
    def optimize_charges_alpha(self):
1336 1372
        """
1337 1373
        Builds the necessary matrix and vector and performs a charge and polarizabilities optimization for this 1 conformation.
@@ -1546,11 +1582,12 @@
Loading
1546 1582
1547 1583
        # Load permanent charges for BCC method
1548 1584
        # For all other methods this is set to 0.0 STILL HAVE to implment
1549 -
        try:
1550 -
            self._conformer.q_am1
1551 -
        except Exception:
1585 +
        if self._conformer.q_am1 is None:
1552 1586
            log.warning('I do not have AM1-type charges')
1553 1587
            self._conformer.q_am1 = np.zeros(self._conformer.natoms)
1588 +
        #else:
1589 +
        #    log.info('Substracting AM1 charge potential from the ESP values')
1590 +
        # ToDo Add Code to substract am1 charge potential
1554 1591
1555 1592
        for j in range(self._conformer.natoms):
1556 1593
            self.e_field_at_atom[0][j] = np.dot(
@@ -1663,7 +1700,7 @@
Loading
1663 1700
        self.e_field_at_atom[2] = 1.0 * self.e[2 * self._conformer.natoms:3 * self._conformer.natoms]
1664 1701
1665 1702
    # Analysation stuff
1666 -
    def calc_esp_q_alpha(self, q_alpha):
1703 +
    def calc_esp_q_alpha(self, q_alpha, mode='q'):
1667 1704
        """
1668 1705
        Calculates the ESP for a given set of point charges and polarizabilities.
1669 1706
@@ -1675,7 +1712,7 @@
Loading
1675 1712
        """
1676 1713
        self.q_pot = np.zeros(len(self.esp_values))
1677 1714
        for i in range(len(self.esp_values)):
1678 -
            if self._conformer._molecule._mode == 'q':
1715 +
            if mode == 'q':
1679 1716
                self.q_pot[i] = np.dot(q_alpha[:self._conformer.natoms], np.transpose(self._conformer.dist)[i])
1680 1717
            else:
1681 1718
                dipole = q_alpha[self._conformer._lines_in_A:]
@@ -1696,6 +1733,11 @@
Loading
1696 1733
            #    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])
1697 1734
            #    self.q_pot[i]=e_dip
1698 1735
1736 +
    def substract_am1_potential(self):
1737 +
        self.calc_esp_q_alpha(self._conformer.q_am1, mode='q')
1738 +
        self.esp_values = Q_(np.subtract(self.esp_values.to('elementary_charge / angstrom').magnitude, self.q_pot),
1739 +
                         'elementary_charge / angstrom')
1740 +
1699 1741
    def sub_esp_q_alpha(self, q_alpha):
1700 1742
        """
1701 1743
        Subtracts the ESP create by a set of point charges and polarizabilities.

Learn more Showing 2 files with coverage changes found.

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