Showing 1 of 1 files from the diff.

@@ -0,0 +1,1706 @@
Loading
1 +
            # !/usr/bin/env python
2 +
3 +
# =============================================================================================
4 +
# MODULE DOCSTRING
5 +
# =============================================================================================
6 +
"""
7 +
This module should convert a mol2 file int a Molecule object and expose all necessary information
8 +
.. todo::
9 +
   * Load in the molecule
10 +
   * Determine the equivalent atoms in the molecule
11 +
   * Determine the bonds based on a smirnoff input file
12 +
"""
13 +
14 +
# =============================================================================================
15 +
# GLOBAL IMPORTS
16 +
# =============================================================================================
17 +
18 +
import logging as log
19 +
import numpy as np
20 +
from openeye import oechem
21 +
import openforcefield.topology as openff
22 +
from openforcefield.typing.engines.smirnoff import ForceField
23 +
from scipy.spatial import distance
24 +
import os
25 +
# Initialize units
26 +
from pint import UnitRegistry
27 +
28 +
ureg = UnitRegistry()
29 +
Q_ = ureg.Quantity
30 +
ureg.define('bohr = 0.52917721067 * angstrom')
31 +
32 +
# =============================================================================================
33 +
# GLOBAL PARAMETERS
34 +
# =============================================================================================
35 +
biological_elements = [1, 3, 5, 6, 7, 8, 9, 11, 12, 14, 15, 16, 17, 19, 20, 34, 35, 53]
36 +
37 +
ROOT_DIR_PATH = os.path.join(os.path.dirname(os.path.realpath(__file__)), '..')
38 +
39 +
40 +
# =============================================================================================
41 +
# PRIVATE SUBROUTINES
42 +
# =============================================================================================
43 +
def find_eq_atoms(mol1):
44 +
    """
45 +
    Finds all equivalent atoms in a molecule.
46 +
    :return
47 +
    Array of pairs of equivalent atoms.
48 +
49 +
    :parameter
50 +
    mol1: Openeye molecule object
51 +
52 +
    TODO Include rdkit support for this function
53 +
    """
54 +
55 +
    qmol = oechem.OEQMol()
56 +
57 +
    # build OEQMol from OEGRAPHMOLECULE
58 +
    oechem.OEBuildMDLQueryExpressions(qmol, mol1)
59 +
    ss2 = oechem.OESubSearch(qmol)
60 +
    oechem.OEPrepareSearch(mol1, ss2)
61 +
62 +
    # Store the equivalent atoms
63 +
    eq_atoms = [[] for i in range(mol1.NumAtoms())]
64 +
65 +
    # Goes through all matches and compares the atom indices.
66 +
    for count, match in enumerate(ss2.Match(mol1)):
67 +
        for ma in match.GetAtoms():
68 +
            # if 2 atoms are not the same atoms
69 +
            if ma.pattern.GetIdx() != ma.target.GetIdx():
70 +
                # and the pair is not stored yet
71 +
                if ma.target.GetIdx() not in eq_atoms[ma.pattern.GetIdx()]:
72 +
                    # save it to the array
73 +
                    eq_atoms[ma.pattern.GetIdx()].append(ma.target.GetIdx())
74 +
75 +
    # goes through the array and returns pairs of equivalent atoms
76 +
    sorted_eq_atoms = []
77 +
    for i, ele in enumerate(eq_atoms):
78 +
        for element in ele:
79 +
            # Making sure we have not already covered this pair of atoms
80 +
            if [element, i] not in sorted_eq_atoms:
81 +
                sorted_eq_atoms.append([i, element])
82 +
    tmp = []
83 +
    for k, ele1 in enumerate(sorted_eq_atoms):
84 +
        exclude = 0
85 +
        for ele2 in sorted_eq_atoms[:k]:
86 +
            if ele1[0] == ele2[0]:
87 +
                exclude = 1
88 +
        if exclude == 0:
89 +
            tmp.append(ele1)
90 +
    sorted_eq_atoms = tmp
91 +
92 +
    return sorted_eq_atoms
93 +
94 +
95 +
96 +
# =============================================================================================
97 +
# TrainingSet
98 +
# =============================================================================================
99 +
100 +
class TrainingSet():
101 +
    """
102 +
    The training set class is the top level class of the resspol program.
103 +
104 +
    It consist of multiple molecule instances and combines the optimization matrices and vectors across multiple molecules.
105 +
    """
106 +
107 +
    def __init__(self, mode='q', scaleparameters=None, scf_scaleparameters=None, SCF=False, thole=False):
108 +
        """
109 +
        Initilize the class and sets the following parameters:
110 +
111 +
        :param mode:
112 +
        :param scaleparameters:
113 +
        :param scf_scaleparameters:
114 +
        :param SCF:
115 +
        :param thole:
116 +
        """
117 +
        self.molecules = list()
118 +
        self.B = np.zeros(0)
119 +
        self.A = np.zeros((0, 0))
120 +
        self.q = 0.0
121 +
        self.mode = mode
122 +
        self.scf_scaleparameters = scf_scaleparameters
123 +
        self.scaleparameters = scaleparameters
124 +
        self._SCF = SCF
125 +
        self._thole = thole
126 +
        self._mode = mode
127 +
        self.step = 0
128 +
        self.number_of_lines_in_X = 0
129 +
        self.X_BCC = 0.0
130 +
        self.Y_BCC = 0.0
131 +
        # Label the atoms and bonds using a offxml file
132 +
        forcefield = ForceField(os.path.join(ROOT_DIR_PATH, 'resppol/data/test_data/BCCPOL.offxml'))
133 +
134 +
        self._nalpha = len(forcefield.get_parameter_handler('vdW').parameters)
135 +
        self._nbccs = len(forcefield.get_parameter_handler('Bonds').parameters)
136 +
        self._bccs ={}
137 +
        for i,element in enumerate(forcefield.get_parameter_handler('Bonds').parameters):
138 +
            self._bccs[element.id] = i
139 +
140 +
141 +
        self._alpha = {}
142 +
        for i,element in enumerate(forcefield.get_parameter_handler('vdW').parameters):
143 +
            self._alpha[element.id] = i
144 +
145 +
        self.bccs_old = np.zeros(self._nbccs)
146 +
        self.alphas_old = np.zeros(self._nalpha)
147 +
148 +
    def load_from_file(self,txtfile):
149 +
        """
150 +
        Allows to build a TrainingSet instance from an text file.
151 +
        File format as followed:
152 +
153 +
        :return:
154 +
        """
155 +
        f = open(datei)
156 +
        lines = f.readlines()
157 +
        f.close()
158 +
        for i, line in enumerate(lines):
159 +
            mol2file = line.split()[1]
160 +
            # noinspection PyTypeChecker
161 +
            self.add_molecule(Molecule(mol2file))
162 +
163 +
    def add_molecule(self, datei):
164 +
        """
165 +
        Adds a molecule.
166 +
        :param datei: Mol2 file of a molecule
167 +
        :return:
168 +
        """
169 +
        self.molecules.append(Molecule(datei, position=self.number_of_lines_in_X, trainingset=self))
170 +
        self.number_of_lines_in_X += self.molecules[-1]._lines_in_X
171 +
172 +
173 +
174 +
175 +
176 +
177 +
178 +
179 +
    def build_matrix_A(self):
180 +
        """
181 +
        Builds the matrix A of the underlying molecules and combines them.
182 +
183 +
        This function is only used for charge optimizatoin RESP
184 +
        """
185 +
        for molecule in self.molecules:
186 +
            molecule.build_matrix_A()
187 +
            X12 = np.zeros((len(self.A), len(molecule.A)))
188 +
            self.A = np.concatenate(
189 +
                (np.concatenate((self.A, X12), axis=1), np.concatenate((X12.transpose(), molecule.A), axis=1)), axis=0)
190 +
191 +
    def build_vector_B(self):
192 +
        for molecule in self.molecules:
193 +
            molecule.build_vector_B()
194 +
            self.B = np.concatenate((self.B, molecule.B))
195 +
196 +
    def build_matrix_X(self):
197 +
        self.X = np.zeros((0,0))
198 +
        """
199 +
        Builds the matrix X of the underlying molecules and combines them.
200 +
201 +
        This function is only used for charge optimizatoin RESP
202 +
        """
203 +
        for molecule in self.molecules:
204 +
            molecule.build_matrix_X()
205 +
            X12 = np.zeros((len(self.X), len(molecule.X)))
206 +
            self.X = np.concatenate(
207 +
                (np.concatenate((self.X, X12), axis=1), np.concatenate((X12.transpose(), molecule.X), axis=1)), axis=0)
208 +
209 +
        X12 = np.zeros((len(self.X),len(self.intramolecular_polarization_rst)))
210 +
        X22 =  np.zeros((len(self.intramolecular_polarization_rst),len(self.intramolecular_polarization_rst)))
211 +
212 +
        self.X = np.concatenate((np.concatenate((self.X, X12), axis=1), np.concatenate((X12.transpose(),X22), axis =1)), axis=0)
213 +
        for i,atoms in enumerate(self.intramolecular_polarization_rst):
214 +
            self.X[self.number_of_lines_in_X + i][atoms[0]] = self.X[atoms[0]][self.number_of_lines_in_X + i] = 1
215 +
            self.X[self.number_of_lines_in_X + i][atoms[1]] = self.X[atoms[1]][self.number_of_lines_in_X + i] = -1
216 +
217 +
218 +
    def build_matrix_X_BCC(self):
219 +
        """
220 +
        Builds the matrix X of the underlying molecules and combines them.
221 +
222 +
        This function is only used for charge optimizatoin RESP
223 +
        """
224 +
225 +
        for molecule in self.molecules:
226 +
            molecule.build_matrix_X_BCC()
227 +
            self.X_BCC += molecule.X
228 +
229 +
230 +
    def build_vector_Y_BCC(self):
231 +
        """
232 +
        Builds the matrix X of the underlying molecules and combines them.
233 +
234 +
        This function is only used for charge optimizatoin RESP
235 +
        """
236 +
237 +
        for molecule in self.molecules:
238 +
            molecule.build_vector_Y_BCC()
239 +
            self.Y_BCC += molecule.Y
240 +
241 +
    def build_vector_Y(self):
242 +
        self.Y = np.zeros(0)
243 +
        for molecule in self.molecules:
244 +
            molecule.build_vector_Y()
245 +
            self.Y = np.concatenate((self.Y, molecule.Y))
246 +
247 +
        Y2 =  np.zeros(len(self.intramolecular_polarization_rst))
248 +
        self.Y = np.concatenate((self.Y, Y2))
249 +
    # def get_intramolecular_charge_rst()
250 +
251 +
    @property
252 +
    def intramolecular_polarization_rst(self):
253 +
254 +
        intramolecular_polarization_rst = []
255 +
        first_occurrence_of_parameter = {}
256 +
        for molecule in self.molecules:
257 +
            first_occurrence_of_parameter_in_molecule = {}
258 +
            for atom in molecule._atoms:
259 +
                if atom._parameter_id not in first_occurrence_of_parameter.keys():
260 +
                    first_occurrence_of_parameter[atom._parameter_id] = molecule._position_in_A + atom._id + molecule._lines_in_A
261 +
                elif atom._parameter_id not in first_occurrence_of_parameter_in_molecule.keys():
262 +
                    intramolecular_polarization_rst.append(
263 +
                        [first_occurrence_of_parameter[atom._parameter_id], molecule._position_in_A + atom._id + molecule._lines_in_A])
264 +
        return intramolecular_polarization_rst
265 +
266 +
    def optimize_charges_alpha(self,criteria = 10E-5):
267 +
        converged = False
268 +
        while converged == False:
269 +
            self.optimize_charges_alpha_step()
270 +
            converged = True
271 +
            for molecule in self.molecules:
272 +
                molecule.step +=1
273 +
                if not all(abs(molecule.q - molecule.q_old) <criteria) or not all(abs(molecule.alpha - molecule.alpha_old) <criteria) :
274 +
                    converged = False
275 +
                molecule.q_old = molecule.q
276 +
                print(molecule.q)
277 +
                print(molecule.alpha)
278 +
                molecule.alpha_old =molecule.alpha
279 +
280 +
281 +
282 +
    def optimize_charges_alpha_step(self):
283 +
        self.build_matrix_X()
284 +
        self.build_vector_Y()
285 +
        self.q_alpha = np.linalg.solve(self.X, self.Y)
286 +
        # Update the charges of the molecules below
287 +
        q_alpha_tmp = self.q_alpha
288 +
        for molecule in self.molecules:
289 +
            molecule.q_alpha = q_alpha_tmp[:len(molecule.X)]
290 +
            q_alpha_tmp = q_alpha_tmp[len(molecule.X):]
291 +
            molecule.update_q_alpha()
292 +
293 +
    def optimize_bcc_alpha(self,criteria = 10E-5):
294 +
        converged = False
295 +
        while converged == False:
296 +
            self.optimize_bcc_alpha_step()
297 +
            for molecule in self.molecules:
298 +
                molecule.step +=1
299 +
                print(molecule.q)
300 +
                print(molecule.alpha)
301 +
            if all(abs(self.bccs - self.bccs_old) < criteria) and all(abs(self.alphas - self.alphas_old)< criteria):
302 +
                converged = True
303 +
            else:
304 +
                self.alphas_old = self.alphas
305 +
                self.bccs_old = self.bccs
306 +
307 +
308 +
    def optimize_bcc_alpha_step(self):
309 +
        self.build_matrix_X_BCC()
310 +
        self.build_vector_Y_BCC()
311 +
312 +
        # Check if bccs or alphas are not in the training set and set them to zero
313 +
        for i,row in enumerate(self.X_BCC):
314 +
            if all(row == 0.0):
315 +
                self.X_BCC[i][i]=1
316 +
317 +
        self.bcc_alpha = np.linalg.solve(self.X_BCC, self.Y_BCC)
318 +
        # Update the charges of the molecules below
319 +
        self.bccs = self.bcc_alpha[:self._nbccs]
320 +
        self.alphas = self.bcc_alpha[self._nbccs:]
321 +
        for molecule in self.molecules:
322 +
            molecule.update_bcc(self.bccs, self.alphas)
323 +
324 +
325 +
326 +
    def optimize_charges(self):
327 +
        self.build_matrix_A()
328 +
        self.build_vector_B()
329 +
        self.q = Q_(np.linalg.solve(self.A, self.B), 'elementary_charge')
330 +
331 +
        # Update the charges of the molecules below
332 +
        q_tmp = self.q
333 +
        for molecule in self.molecules:
334 +
            molecule.q = q_tmp[:len(molecule.A)]
335 +
            q_tmp = q_tmp[len(molecule.A):]
336 +
            molecule.update_q()
337 +
338 +
339 +
# =============================================================================================
340 +
# Molecule
341 +
# =============================================================================================
342 +
343 +
344 +
class Molecule:
345 +
    """
346 +
        This class loads in a mol2 file and expose all relevant information. It combines the functionality of
347 +
        an openeye molecule with the functionality of the OFF molecule. The OE part is only needed to determine
348 +
        the chemical equivalent atoms. When this feature is implemented in OFF toolkit the OE molecule is not
349 +
        necessary anymore
350 +
351 +
        :param datei: mol2 file
352 +
353 +
        :return:
354 +
        """
355 +
356 +
    def __init__(self, datei, position=0, trainingset=None):
357 +
358 +
        # Get Molecle name from filename
359 +
        self.B = 0.0
360 +
        self.A = 0.0
361 +
        self.X = 0.0
362 +
        self.Y = 0.0
363 +
        self._name = datei.split('/')[-1].strip(".mol2")
364 +
        self._trainingset = trainingset
365 +
366 +
        # Number of optimization steps
367 +
        self.step = 0
368 +
369 +
        # Postion of this molecule in optimization matrix A
370 +
        self._position_in_A = position
371 +
372 +
        # Copy (scf) scaleparameters from the trainingset definition
373 +
        if trainingset is None:
374 +
            self.scf_scaleparameters = None
375 +
            self.scaleparameters = None
376 +
            self._thole = False
377 +
            self._SCF = False
378 +
            self._mode = 'q'
379 +
        else:
380 +
            self.scf_scaleparameters = self._trainingset.scf_scaleparameters
381 +
            self.scaleparameters = self._trainingset.scaleparameters
382 +
            self._thole = self._trainingset._thole
383 +
            self._SCF = self._trainingset._SCF
384 +
            self._mode = self._trainingset._mode
385 +
            # Initialize OE Molecule
386 +
        self.oemol = oechem.OEMol()
387 +
388 +
        # Open Input File Stream
389 +
        ifs = oechem.oemolistream(datei)
390 +
391 +
        # Read from IFS to molecule
392 +
        oechem.OEReadMol2File(ifs, self.oemol)
393 +
394 +
        # Check if molecule is a 3 dimensional object
395 +
        if self.oemol.GetDimension() != 3:
396 +
            log.error('Molecule either not found or does not have a 3D structure')
397 +
            raise Exception('Molecule either not found or does not have a 3D structure')
398 +
399 +
        # Check for strange atoms
400 +
        for atom in self.oemol.GetAtoms():
401 +
            if atom.GetAtomicNum() not in biological_elements:
402 +
                log.warning("""I detect an atom with atomic number: {}
403 +
                Are you sure you want to include such an atom in your dataset?
404 +
                Maybe you are using a mol2 file with amber atomtypes""".format(atom.GetAtomicNum()))
405 +
                raise Warning("""I detect an atom with atomic number: {}
406 +
                Are you sure you want to include such an atom in your dataset?
407 +
                Maybe you are using a mol2 file with amber atomtypes""".format(atom.GetAtomicNum()))
408 +
409 +
        # Generate the OFF molecule from the openeye molecule
410 +
        self.offmol = openff.Molecule.from_openeye(self.oemol)
411 +
        self.offtop = openff.Topology.from_molecules([self.offmol])
412 +
413 +
        # Label the atoms and bonds using a offxml file
414 +
        forcefield = ForceField(os.path.join(ROOT_DIR_PATH, 'resppol/data/test_data/BCCPOL.offxml'))
415 +
        if trainingset is None:
416 +
            self._trainingset = TrainingSet()
417 +
418 +
419 +
        # Run the parameter labeling
420 +
        molecule_parameter_list = forcefield.label_molecules(self.offtop)
421 +
422 +
        # set molecule charge
423 +
        # self._charge=openff.Molecule.total_charge
424 +
        self._charge = 0
425 +
426 +
        # Initialize the bonds, atoms
427 +
        self._bonds = list()
428 +
        self._atoms = list()
429 +
430 +
        # Define all BCC bonds
431 +
        for i, properties in enumerate(molecule_parameter_list[0]['Bonds'].items()):
432 +
            atom_indices, parameter = properties
433 +
            self.add_bond(i, atom_indices, parameter.id)
434 +
435 +
        self._nbonds = len(self._bonds)
436 +
437 +
        # Define all atomtypes for polarization
438 +
        for i, properties in enumerate(molecule_parameter_list[0]['vdW'].items()):
439 +
            atom_index, parameter = properties
440 +
            self.add_atom(i, atom_index, parameter.id)
441 +
442 +
        self._natoms = len(self._atoms)
443 +
444 +
        # Initialize and fill scaling matrix
445 +
        self.scale = np.ones((self._natoms, self._natoms))
446 +
        self.scale_scf = np.ones((self._natoms, self._natoms))
447 +
        self.scaling(scf_scaleparameters=self.scf_scaleparameters, scaleparameters=self.scaleparameters)
448 +
        self.create_BCCmatrix_T()
449 +
        self.create_POLmatrix_R()
450 +
451 +
        # Initialize conformers
452 +
        self.conformers = list()
453 +
454 +
        # Number of lines for matrix X
455 +
        if self._mode == 'q':
456 +
            self._lines_in_X = self._natoms + len(self.chemical_eq_atoms) + 1
457 +
            self.q_old = np.zeros(self._natoms)
458 +
        if self._mode == 'q_alpha':
459 +
            self._lines_in_X = self._natoms + len(
460 +
                self.chemical_eq_atoms) + 1 + 3 * self._natoms + 2 * self._natoms
461 +
            self.q_old = np.zeros(self._natoms)
462 +
            self.alpha_old = np.zeros(3*self._natoms)
463 +
        self._lines_in_A =  self._natoms + len(self.chemical_eq_atoms) + 1
464 +
        # Initiliaxe charges
465 +
        self.q_alpha = np.zeros(self._lines_in_X)
466 +
467 +
    def add_bond(self, index, atom_indices, parameter_id):
468 +
        """
469 +
        Adds a bond object ot the molecule
470 +
471 +
        :param index:
472 +
        :param atom_indices:
473 +
        :param parameter_id:
474 +
        :return:
475 +
        """
476 +
        atom_index1 = atom_indices[0]
477 +
        atom_index2 = atom_indices[1]
478 +
        self._bonds.append(Bond(index, atom_index1, atom_index2, parameter_id))
479 +
480 +
    def add_atom(self, index, atom_index, parameter_id):
481 +
        """
482 +
        Adds a atom object to the molecule
483 +
484 +
        :param index:
485 +
        :param atom_index:
486 +
        :param parameter_id:
487 +
        :return:
488 +
        """
489 +
        self._atoms.append(Atom(index, atom_index[0], parameter_id))
490 +
491 +
    def add_conformer_from_mol2(self, mol2file):
492 +
        """
493 +
        Adds a conformer from a mol2 file
494 +
        Automatically checks if the conformer corresponds to the molecule
495 +
496 +
        :param mol2file:
497 +
        :return:
498 +
        """
499 +
        conf = openff.Molecule.from_file(mol2file)
500 +
501 +
        # Check if molecule is a 3 dimensional object and has the correct dimensions
502 +
        # Checks if this conformer is from this molecule based on atom names
503 +
        if self.offmol.n_atoms != conf.n_atoms or \
504 +
                self.offmol.n_bonds != conf.n_bonds or \
505 +
                self.offmol.n_angles != conf.n_angles:
506 +
            log.error('Conformer and Molecule does not match')
507 +
            raise Exception('Conformer and Molecule does not match')
508 +
        else:
509 +
            self.conformers.append(Conformer(conf, molecule=self))
510 +
        # Checks if this conformer is from this moleule based on atom names
511 +
512 +
    @property
513 +
    def chemical_eq_atoms(self):
514 +
        """
515 +
        # Note: Something similar in compare_charges, which can be copied and slightly modified
516 +
        :return:
517 +
        """
518 +
        return find_eq_atoms(self.oemol)
519 +
520 +
    @property
521 +
    def same_polarization_atoms(self):
522 +
        array_of_same_pol_atoms = []
523 +
        first_occurrence = {}
524 +
        for atom in self._atoms:
525 +
            if atom._parameter_id not in first_occurrence.keys():
526 +
                first_occurrence[atom._parameter_id] = atom._id
527 +
            else:
528 +
                array_of_same_pol_atoms.append([first_occurrence[atom._parameter_id], atom._id ])
529 +
        return (array_of_same_pol_atoms)
530 +
531 +
532 +
    def create_BCCmatrix_T(self):
533 +
        """
534 +
        Creates the bond matrix T for a molecule.
535 +
536 +
        Parameters:
537 +
        ----------
538 +
            bondtyps: list of [int,int]
539 +
                Bonds in the set between different BCC groups
540 +
541 +
        Attribute:
542 +
        ----------
543 +
            bondmatrix T: 2 dim array
544 +
                1 dim atom
545 +
                2 dim BCC
546 +
        See also AM1-BCC paper:
547 +
        https://onlinelibrary.wiley.com/doi/epdf/10.1002/%28SICI%291096-987X%2820000130%2921%3A2%3C132%3A%3AAID-JCC5%3E3.0.CO%3B2-P
548 +
        """
549 +
        bondsexcluded = []
550 +
        self.T = np.zeros((self._natoms, self._trainingset._nbccs))
551 +
        for bond in self._bonds:
552 +
            if bond._parameter_id not in bondsexcluded:
553 +
                    self.T[bond._atom1][self._trainingset._bccs[bond._parameter_id]] += 1
554 +
                    self.T[bond._atom2][self._trainingset._bccs[bond._parameter_id]] -= 1
555 +
556 +
    def create_POLmatrix_R(self):
557 +
        """
558 +
        Create a polarization matrix. Defines which atom has which polarizability parameter.
559 +
560 +
        Parameters:
561 +
        ----------
562 +
        groups: dictionary atomtyp -> int
563 +
            Stores the polarization group  of each atom
564 +
565 +
        Defines which atom belongs to the which pol group.
566 +
        Takes the list of atomtypes and assign the corresponding group value.
567 +
        Here I am using a slightly different approach in contrast to the atomtyps.
568 +
        as the number of different pol types is maximal the number of atom types. No pre-screening for the involved
569 +
        atom-types is needed
570 +
        """
571 +
        self.R = np.zeros((self._natoms, self._trainingset._nalpha))
572 +
        for atom in self._atoms:
573 +
            self.R[atom._id][self._trainingset._alpha[atom._parameter_id]] += 1
574 +
575 +
    def update_bcc(self, bccs, alphas):
576 +
        """
577 +
        Converts the optimized bccs and alphas to a qd object.
578 +
579 +
        :param mol2: molecule class object
580 +
        :return:
581 +
        """
582 +
        for i in range(self._natoms):
583 +
            self.q_alpha[i] = np.dot(self.T[i], bccs)
584 +
        self.alpha = [np.dot(self.R[i], alphas) for i in range(len(self.R))]
585 +
        self.q_alpha[self._lines_in_A:self._lines_in_A + 3* len(self.alpha)] = np.concatenate(
586 +
            (np.concatenate((self.alpha, self.alpha)), self.alpha))
587 +
        self.q = self.q_alpha[:self._natoms]
588 +
        for conformer in self.conformers:
589 +
            conformer.q_alpha = self.q_alpha
590 +
591 +
592 +
    def build_matrix_A(self):
593 +
        for conformer in self.conformers:
594 +
            conformer.build_matrix_A()
595 +
            self.A += conformer.A
596 +
597 +
    def build_vector_B(self):
598 +
        for conformer in self.conformers:
599 +
            conformer.build_vector_B()
600 +
            self.B += conformer.B
601 +
602 +
    def build_matrix_X(self):
603 +
        for conformer in self.conformers:
604 +
            conformer.build_matrix_X()
605 +
            self.X += conformer.X
606 +
607 +
    def build_matrix_X_BCC(self):
608 +
        for conformer in self.conformers:
609 +
            conformer.build_matrix_X_BCC()
610 +
            self.X += conformer.X
611 +
612 +
    def build_vector_Y(self):
613 +
        for conformer in self.conformers:
614 +
            conformer.build_vector_Y()
615 +
            self.Y += conformer.Y
616 +
617 +
    def build_vector_Y_BCC(self):
618 +
        for conformer in self.conformers:
619 +
            conformer.build_vector_Y_BCC()
620 +
            self.Y += conformer.Y
621 +
622 +
    def optimize_charges(self):
623 +
        self.build_matrix_A()
624 +
        self.build_vector_B()
625 +
        self.q = Q_(np.linalg.solve(self.A, self.B), 'elementary_charge')
626 +
627 +
    def optimize_charges_alpha(self):
628 +
        self.build_matrix_X()
629 +
        self.build_vector_Y()
630 +
        self.q_alpha = Q_(np.linalg.solve(self.X, self.Y), 'elementary_charge')
631 +
632 +
    def update_q(self):
633 +
        for conformer in self.conformers:
634 +
            conformer.q = self.q
635 +
636 +
    def update_q_alpha(self):
637 +
        self.q = self.q_alpha[:self._natoms]
638 +
        self.alpha = self.q_alpha[self._natoms+1+len(self.chemical_eq_atoms):4*self._natoms+1+len(self.chemical_eq_atoms):]
639 +
        for conformer in self.conformers:
640 +
            conformer.q_alpha = self.q_alpha
641 +
642 +
    def scaling(self, scaleparameters=None, scf_scaleparameters=None):
643 +
        """
644 +
        Takes the bond information from a molecule instances and converts it to an scaling matrix.
645 +
646 +
        Parameters:
647 +
        ---------
648 +
        bonds:  list of [int,int]
649 +
            list of atoms connected with each other
650 +
        scaleparameters: [float,float,float]
651 +
            1-2 scaling, 1-3 scaling, 1-4 scaling parameter
652 +
653 +
        Attributes:
654 +
        ---------
655 +
        scale: matrix
656 +
            scaling matrix how atoms interact with each other
657 +
658 +
        """
659 +
660 +
        # Initializing
661 +
        bound12 = np.zeros((self._natoms, self._natoms))
662 +
        bound13 = np.zeros((self._natoms, self._natoms))
663 +
        bound14 = np.zeros((self._natoms, self._natoms))
664 +
        if scaleparameters is None:
665 +
            scaleparameters = [0.0, 0.0, 0.8333333333]
666 +
667 +
        # Building connection matrix
668 +
        for bond in self._bonds:
669 +
            bound12[bond._atom1][bond._atom2] = 1.0
670 +
            bound12[bond._atom2][bond._atom1] = 1.0
671 +
672 +
        for i in range(len(bound12)):
673 +
            b12 = np.where(bound12[i] == 1.0)[0]
674 +
            for j in range(len(b12)):
675 +
                b12t = np.where(bound12[b12[j]] == 1.0)[0]
676 +
                for k in range(len(b12t)):
677 +
                    if i != b12t[k]:
678 +
                        bound13[b12t[k]][i] = 1.0
679 +
                        bound13[i][b12t[k]] = 1.0
680 +
681 +
        for i in range(self._natoms):
682 +
            b13 = np.where(bound13[i] == 1.0)[0]
683 +
            for j in range(len(b13)):
684 +
                b13t = np.where(bound12[b13[j]] == 1.0)[0]
685 +
                for k in range(len(b13t)):
686 +
                    if bound12[b13t[k]][i] == 0.0:
687 +
                        bound14[b13t[k]][i] = 1.0
688 +
                        bound14[i][b13t[k]] = 1.0
689 +
690 +
        for i in range(self._natoms):
691 +
            self.scale[i][i] = 0.0
692 +
        # find values in matrix with value 1.0
693 +
        b12 = np.array(np.where(bound12 == 1.0)).transpose()
694 +
        b13 = np.array(np.where(bound13 == 1.0)).transpose()
695 +
        b14 = np.array(np.where(bound14 == 1.0)).transpose()
696 +
697 +
        # Fill scaling matrix with values
698 +
        for i in range(len(b12)):
699 +
            self.scale[b12[i][0]][b12[i][1]] = scaleparameters[
700 +
                0]  # Value for 1-2 interaction 0 means interactions are neglected
701 +
        for i in range(len(b13)):
702 +
            self.scale[b13[i][0]][b13[i][1]] = scaleparameters[
703 +
                1]  # Value for 1-3 interaction 0 means interactions are neglected
704 +
        for i in range(len(b14)):
705 +
            self.scale[b14[i][0]][b14[i][1]] = scaleparameters[2]  # Value for the 1-4 scaling
706 +
707 +
        # Different Scaling parameter for SCF
708 +
        if scf_scaleparameters != None:
709 +
            self.scale_scf = np.ones((self._natoms, self._natoms))
710 +
            for i in range(self._natoms):
711 +
                self.scale_scf[i][i] = 0.0
712 +
            for i in range(len(b12)):
713 +
                self.scale_scf[b12[i][0]][b12[i][1]] = scf_scaleparameters[
714 +
                    0]  # Value for 1-2 interaction 0 means interactions are neglected
715 +
            for i in range(len(b13)):
716 +
                self.scale_scf[b13[i][0]][b13[i][1]] = scf_scaleparameters[
717 +
                    1]  # Value for 1-3 interaction 0 means interactions are neglected
718 +
            for i in range(len(b14)):
719 +
                self.scale_scf[b14[i][0]][b14[i][1]] = scf_scaleparameters[2]  # Value for the 1-4 scaling
720 +
721 +
722 +
# =============================================================================================
723 +
# Atom
724 +
# =============================================================================================
725 +
726 +
class Atom:
727 +
    """
728 +
729 +
    """
730 +
731 +
    def __init__(self, index, atom_index, parameter_id):
732 +
        self._id = index
733 +
        self._atom = atom_index
734 +
        self._parameter_id = parameter_id
735 +
736 +
737 +
# =============================================================================================
738 +
# Bond
739 +
# =============================================================================================
740 +
741 +
class Bond:
742 +
743 +
    def __init__(self, index, atom_index1, atom_index2, parameter_id):
744 +
        self._id = index
745 +
        self._atom1 = atom_index1
746 +
        self._atom2 = atom_index2
747 +
        self._parameter_id = parameter_id
748 +
749 +
750 +
# =============================================================================================
751 +
# Conformer
752 +
# =============================================================================================
753 +
754 +
class Conformer:
755 +
    """
756 +
757 +
    """
758 +
759 +
    def __init__(self, conf, molecule=None):
760 +
        """
761 +
762 +
        :param conf: openff molecule file
763 +
        """
764 +
        # noinspection PyProtectedMember
765 +
        self.atom_positions = Q_(np.array(conf.conformers[0]._value), 'angstrom')
766 +
        self.atom_positions_angstrom = self.atom_positions.to('angstrom').magnitude
767 +
        self.natoms = len(self.atom_positions.magnitude)
768 +
        self.baseESP = None
769 +
        self.polESPs = list()
770 +
        self._molecule = molecule
771 +
        self.q_alpha = self._molecule.q_alpha
772 +
        self._lines_in_A = self._molecule._lines_in_A
773 +
774 +
        # Initiliaze Electric field vectors
775 +
        self.e_field_at_atom = np.zeros((3, self.natoms))
776 +
777 +
    def get_grid_coord(self):
778 +
        self.grid_coord_angstrom = self.baseESP.positions.to('angstrom').magnitude
779 +
        self.npoints = len(self.baseESP.positions.magnitude)
780 +
781 +
    def add_baseESP(self, *args, ):
782 +
        """
783 +
        Adds the unpolarized molecule to this conformation.
784 +
785 +
        :param args:  ESPF file
786 +
        1.) GESP file form g09
787 +
        2.) grid.dat file and esp.dat file generated via respyte and psi4
788 +
789 +
        :return:
790 +
        """
791 +
        self.baseESP = BCCUnpolESP(*args, conformer=self)
792 +
793 +
        # Check if atomic coordinates match
794 +
        for atom1 in self.baseESP.atom_positions:
795 +
            atom_is_present = 0
796 +
            for atom2 in self.atom_positions:
797 +
                if np.linalg.norm(atom1.to('angstrom').magnitude - atom2.to('angstrom').magnitude) < 0.01:
798 +
                    atom_is_present = 1
799 +
            if atom_is_present == 0:
800 +
                raise Exception("ESP grid does not belong to the conformer")
801 +
802 +
    def add_polESP(self, *args, e_field=Q_([0.0, 0.0, 0.0], 'bohr')):
803 +
        """
804 +
        Adds the unpolarized molecule to this conformation.
805 +
806 +
        :param args:  ESPF file
807 +
        1.) GESP file form g09
808 +
        2.) grid.dat file and esp.dat file generated via respyte and psi4
809 +
        :param:e_field: Pint formatted electrif field
810 +
        e.g e_field=Q_([0.0, 0.0, 0.0], 'bohr'))
811 +
812 +
        :return:
813 +
        """
814 +
        self.polESPs.append(BCCPolESP(*args, conformer=self, e_field=e_field))
815 +
816 +
    # Build the matrix A for the charge optimization
817 +
    def build_matrix_A(self):
818 +
819 +
        """
820 +
        Fast method for only optimizing charges.
821 +
822 +
        :return:
823 +
        """
824 +
        # Determine size of matrix for this molecule
825 +
        # every atom is one line
826 +
        # one line to restrain the overall charge of the molecule
827 +
        # one line for every pair of equivalent atoms
828 +
        self.get_distances()
829 +
        self.A = np.zeros((self._lines_in_A, self._lines_in_A))
830 +
        for j in range(self.natoms):
831 +
            for k in range(j + 1):
832 +
                self.A[j][k] = np.dot(self.dist[j], self.dist[k])
833 +
834 +
        # Symmetric matrix -> copy diagonal elements, add total charge restrain
835 +
        for j in range(self.natoms):
836 +
            for k in range(j):
837 +
                self.A[k][j] = self.A[j][k]
838 +
            self.A[self.natoms][j] = 1.0
839 +
            self.A[j][self.natoms] = 1.0
840 +
841 +
        # Add charge restraints for equivalent atoms
842 +
        for k, eq_atoms in enumerate(self._molecule.chemical_eq_atoms):
843 +
            if eq_atoms[1] > 0:
844 +
                self.A[self.natoms + 1 + k][eq_atoms[0]] = 1
845 +
                self.A[self.natoms + 1 + k][eq_atoms[1]] = -1
846 +
                self.A[eq_atoms[0]][self.natoms + 1 + k] = 1
847 +
                self.A[eq_atoms[1]][self.natoms + 1 + k] = -1
848 +
849 +
    def build_matrix_D(self):
850 +
        """
851 +
        Method for building the polarization matrix D. Only used for fitting polarizations withotut charges or BCCs.
852 +
853 +
        :return:
854 +
        """
855 +
        # 1 dipole vector of length 3 per atom
856 +
        # restrains for isotropic polarization
857 +
        # Restaint atoms with same polarization parameters
858 +
        self.Dlines = (3 * self.natoms + 2 * self.natoms) #+ len(self._molecule.same_polarization_atoms))
859 +
        self.D = np.zeros((self.Dlines, self.Dlines))
860 +
861 +
        for j in range(self.natoms):
862 +
            for k in range(self.natoms):
863 +
                self.D[k][j] = np.multiply(np.dot(self.dist_x[k], self.dist_x[j]),
864 +
                                           self.e_field_at_atom[0][k] * self.e_field_at_atom[0][j])
865 +
                self.D[j + self.natoms][k] = self.D[k][j + self.natoms] = np.multiply(
866 +
                    np.dot(self.dist_x[k], self.dist_y[j]), self.e_field_at_atom[0][k] * self.e_field_at_atom[1][j])
867 +
                self.D[j + 2 * self.natoms][k] = self.D[k][j + 2 * self.natoms] = np.multiply(
868 +
                    np.dot(self.dist_x[k], self.dist_z[j]), self.e_field_at_atom[0][k] * self.e_field_at_atom[2][j])
869 +
                self.D[k + self.natoms][j + self.natoms] = np.multiply(np.dot(self.dist_y[k], self.dist_y[j]),
870 +
                                                                       self.e_field_at_atom[1][k] *
871 +
                                                                       self.e_field_at_atom[1][j])
872 +
                self.D[j + 2 * self.natoms][k + self.natoms] = self.D[k + self.natoms][
873 +
                    j + 2 * self.natoms] = np.multiply(np.dot(self.dist_y[k], self.dist_z[j]),
874 +
                                                       self.e_field_at_atom[1][k] * self.e_field_at_atom[2][j])
875 +
                self.D[k + 2 * self.natoms][j + 2 * self.natoms] = np.multiply(
876 +
                    np.dot(self.dist_z[k], self.dist_z[j]), self.e_field_at_atom[2][k] * self.e_field_at_atom[2][j])
877 +
        # Add dipole restraints for equivalent atoms /only works for isotropic suff now
878 +
        for polESP in self.polESPs:
879 +
            for j in range(self.natoms):
880 +
                for k in range(self.natoms):
881 +
                    self.D[k][j] += np.multiply(np.dot(self.dist_x[k], self.dist_x[j]),
882 +
                                                polESP.e_field_at_atom[0][k] * polESP.e_field_at_atom[0][j])
883 +
                    self.D[j + self.natoms][k] = self.D[k][j + self.natoms] = self.D[k][j + self.natoms] + np.multiply(
884 +
                        np.dot(self.dist_x[k], self.dist_y[j]),
885 +
                        polESP.e_field_at_atom[0][k] * polESP.e_field_at_atom[1][j])
886 +
                    self.D[j + 2 * self.natoms][k] = self.D[k][j + 2 * self.natoms] = self.D[k][
887 +
                                                                                          j + 2 * self.natoms] + np.multiply(
888 +
                        np.dot(self.dist_x[k], self.dist_z[j]),
889 +
                        polESP.e_field_at_atom[0][k] * polESP.e_field_at_atom[2][j])
890 +
                    self.D[k + self.natoms][j + self.natoms] += np.multiply(np.dot(self.dist_y[k], self.dist_y[j]),
891 +
                                                                            polESP.e_field_at_atom[1][k] *
892 +
                                                                            polESP.e_field_at_atom[1][j])
893 +
                    self.D[j + 2 * self.natoms][k + self.natoms] = self.D[k + self.natoms][j + 2 * self.natoms] = \
894 +
                    self.D[k + self.natoms][j + 2 * self.natoms] + np.multiply(np.dot(self.dist_y[k], self.dist_z[j]),
895 +
                                                                               polESP.e_field_at_atom[1][k] *
896 +
                                                                               polESP.e_field_at_atom[2][j])
897 +
                    self.D[k + 2 * self.natoms][j + 2 * self.natoms] += np.multiply(
898 +
                        np.dot(self.dist_z[k], self.dist_z[j]),
899 +
                        polESP.e_field_at_atom[2][k] * polESP.e_field_at_atom[2][j])
900 +
901 +
        self.D = self.D / (len(self.polESPs) + 1)
902 +
903 +
        #for j, atoms in enumerate(self._molecule.same_polarization_atoms):
904 +
        #    self.D[5 * self.natoms + j][atoms[0]] = 1
905 +
        #    self.D[5 * self.natoms + j][atoms[1]] = -1
906 +
        #    self.D[atoms[0]][5 * self.natoms + j] = 1
907 +
        #    self.D[atoms[1]][5 * self.natoms + j] = -1
908 +
        # Code to keep polarization parameters at their initial value. Implmentation will change
909 +
        # elif self.eqdipoles[j][1] < 0:
910 +
        #    self.D[self.ndipoles + self.aniso_lines + j][self.eqdipoles[j][0]] = 1
911 +
        #    self.D[self.eqdipoles[j][0]][self.ndipoles + self.aniso_lines + j] = 1
912 +
913 +
        # Add restraints for polarization isotropy
914 +
        for j in range(self.natoms):
915 +
            self.D[3 * self.natoms + j][j] = self.D[j][3 * self.natoms + j] = 1.0
916 +
            self.D[3 * self.natoms + j][j + self.natoms] = self.D[j + self.natoms][3 * self.natoms + j] = -1.0
917 +
            self.D[4 * self.natoms + j][j] = self.D[j][4 * self.natoms + j] = 1.0
918 +
            self.D[4 * self.natoms + j][j + 2 * self.natoms] = self.D[j + 2 * self.natoms][4 * self.natoms + j] = -1.0
919 +
920 +
    def build_matrix_X(self):
921 +
        """
922 +
        Creates Matrix X for the RESP-POl method.
923 +
924 +
        RESP and Polarization with the large matrix.
925 +
        Probably worth changing it to the new model.
926 +
927 +
        Again the math is shown in the manuscript.
928 +
        """
929 +
        self.build_matrix_A()
930 +
        self.get_electric_field()
931 +
        self.build_matrix_D()
932 +
        self.B = np.zeros((self._lines_in_A, self.Dlines))
933 +
        self.C = np.zeros((self.Dlines, self._lines_in_A))
934 +
935 +
        # Matrix element B see notes
936 +
        for k in range(self.natoms):
937 +
            for j in range(self.natoms):
938 +
                self.B[k][j] = np.multiply(np.dot(self.dist[k], self.dist_x[j]), self.e_field_at_atom[0][j])  # B1
939 +
                self.B[k][self.natoms + j] = np.multiply(np.dot(self.dist[k], self.dist_y[j]),
940 +
                                                         self.e_field_at_atom[1][j])  # B2
941 +
                self.B[k][2 * self.natoms + j] = np.multiply(np.dot(self.dist[k], self.dist_z[j]),
942 +
                                                             self.e_field_at_atom[2][j])  # B3
943 +
944 +
        for polESP in self.polESPs:
945 +
            for k in range(self.natoms):
946 +
                for j in range(self.natoms):
947 +
                    self.B[k][j] += np.multiply(np.dot(self.dist[k], self.dist_x[j]),
948 +
                                                polESP.e_field_at_atom[0][j])  # B1
949 +
                    self.B[k][self.natoms + j] += np.multiply(np.dot(self.dist[k], self.dist_y[j]),
950 +
                                                              polESP.e_field_at_atom[1][j])  # B2
951 +
                    self.B[k][2 * self.natoms + j] += np.multiply(np.dot(self.dist[k], self.dist_z[j]),
952 +
                                                                  polESP.e_field_at_atom[2][j])  # B3
953 +
        # matrix element C see notes
954 +
        # matrix element C
955 +
        for j in range(self.natoms):
956 +
            for k in range(self.natoms):
957 +
                self.C[k][j] = np.multiply(np.dot(self.dist[j], self.dist_x[k]), self.e_field_at_atom[0][k])
958 +
                self.C[self.natoms + k][j] = np.multiply(np.dot(self.dist[j], self.dist_y[k]),
959 +
                                                         self.e_field_at_atom[1][k])
960 +
                self.C[2 * self.natoms + k][j] = np.multiply(np.dot(self.dist[j], self.dist_z[k]),
961 +
                                                             self.e_field_at_atom[2][k])
962 +
963 +
        for polESP in self.polESPs:
964 +
            for j in range(self.natoms):
965 +
                for k in range(self.natoms):
966 +
                    self.C[k][j] += np.multiply(np.dot(self.dist[j], self.dist_x[k]), polESP.e_field_at_atom[0][k])
967 +
                    self.C[self.natoms + k][j] += np.multiply(np.dot(self.dist[j], self.dist_y[k]),
968 +
                                                              polESP.e_field_at_atom[1][k])
969 +
                    self.C[2 * self.natoms + k][j] += np.multiply(np.dot(self.dist[j], self.dist_z[k]),
970 +
                                                                  polESP.e_field_at_atom[2][k])
971 +
        # Normalize B and C
972 +
        self.B = self.B / (len(self.polESPs) + 1)
973 +
        self.C = self.C / (len(self.polESPs) + 1)
974 +
975 +
        # Combine all matrices
976 +
        self.X = np.concatenate(
977 +
            (np.concatenate((self.A, self.B), axis=1), np.concatenate((self.C, self.D), axis=1)), axis=0)
978 +
979 +
980 +
    # Combines BCC and polarizabilities
981 +
    def build_matrix_X_BCC(self):
982 +
        """
983 +
        Creates the Matrix A for the BCC-POL approach:
984 +
985 +
        :param mol2: molecule class object
986 +
987 +
        :return:
988 +
989 +
        The math and the optimization method is explained in the current paper draft.
990 +
        """
991 +
        self.get_distances()
992 +
        nbcc = self._molecule._trainingset._nbccs
993 +
        nalpha = self._molecule._trainingset._nalpha
994 +
        T = self._molecule.T
995 +
        R = self._molecule.R
996 +
        self.get_electric_field()
997 +
        # BCC part
998 +
        self.A = np.zeros((nbcc, nbcc))
999 +
        if self._molecule._mode == 'alpha':  # Do not optimize BCCs in that case
1000 +
            for alpha in range(nbcc):
1001 +
                self.A[alpha][alpha] = 1
1002 +
        else:
1003 +
            for ESP in [self.baseESP] + self.polESPs:
1004 +
                for j in range(self.natoms):
1005 +
                    for k in range(self.natoms):
1006 +
                        for alpha in range(nbcc):
1007 +
                            for beta in range(nbcc):
1008 +
                                self.A[alpha][beta] += T[j][alpha] * T[k][beta] * np.dot(self.dist[j],
1009 +
                                                                                               self.dist[k])
1010 +
        # Polarizabilities part
1011 +
        self.D = np.zeros((nalpha, nalpha))
1012 +
        if self._molecule._mode != 'bcconly':  # Mode is not optimizing polarizabilies
1013 +
            for ESP in [self.baseESP] + self.polESPs:
1014 +
                for j in range(self.natoms):
1015 +
                    for k in range(self.natoms):
1016 +
                        for alpha in range(nalpha):
1017 +
                            for beta in range(nalpha):
1018 +
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1019 +
                                    np.dot(self.dist_x[j], self.dist_x[k]), ESP.e_field_at_atom[0][j] * ESP.e_field_at_atom[0][k])
1020 +
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1021 +
                                    np.dot(self.dist_x[j], self.dist_y[k]), ESP.e_field_at_atom[0][j] * ESP.e_field_at_atom[1][k])
1022 +
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1023 +
                                    np.dot(self.dist_x[j], self.dist_z[k]), ESP.e_field_at_atom[0][j] * ESP.e_field_at_atom[2][k])
1024 +
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1025 +
                                    np.dot(self.dist_y[j], self.dist_x[k]), ESP.e_field_at_atom[1][j] * ESP.e_field_at_atom[0][k])
1026 +
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1027 +
                                    np.dot(self.dist_y[j], self.dist_y[k]), ESP.e_field_at_atom[1][j] * ESP.e_field_at_atom[1][k])
1028 +
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1029 +
                                    np.dot(self.dist_y[j], self.dist_z[k]), ESP.e_field_at_atom[1][j] * ESP.e_field_at_atom[2][k])
1030 +
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1031 +
                                    np.dot(self.dist_z[j], self.dist_x[k]), ESP.e_field_at_atom[2][j] * ESP.e_field_at_atom[0][k])
1032 +
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1033 +
                                    np.dot(self.dist_z[j], self.dist_y[k]), ESP.e_field_at_atom[2][j] * ESP.e_field_at_atom[1][k])
1034 +
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1035 +
                                    np.dot(self.dist_z[j], self.dist_z[k]), ESP.e_field_at_atom[2][j] * ESP.e_field_at_atom[2][k])
1036 +
        else:
1037 +
            for alpha in range(nalpha):
1038 +
                self.D[alpha][alpha] = 1
1039 +
1040 +
        # Cross interaction between BCC charges and polarizations
1041 +
        self.B = np.zeros((nbcc, nalpha))
1042 +
        self.C = np.zeros((nalpha, nbcc))
1043 +
        if self._molecule._mode != 'alpha':
1044 +
            for ESP in [self.baseESP] + self.polESPs:
1045 +
                for j in range(self.natoms):
1046 +
                    for k in range(self.natoms):
1047 +
                        for alpha in range(nbcc):
1048 +
                            for beta in range(nalpha):
1049 +
                                self.B[alpha][beta] += T[j][alpha] * R[k][beta] * np.multiply(
1050 +
                                    np.dot(self.dist[j], self.dist_x[k]), ESP.e_field_at_atom[0][k])
1051 +
                                self.B[alpha][beta] += T[j][alpha] * R[k][beta] * np.multiply(
1052 +
                                    np.dot(self.dist[j], self.dist_y[k]), ESP.e_field_at_atom[1][k])
1053 +
                                self.B[alpha][beta] += T[j][alpha] * R[k][beta] * np.multiply(
1054 +
                                    np.dot(self.dist[j], self.dist_z[k]), ESP.e_field_at_atom[2][k])
1055 +
        if self._molecule._mode != 'bcconly':
1056 +
            for ESP in [self.baseESP] + self.polESPs:
1057 +
                for j in range(self.natoms):
1058 +
                    for k in range(self.natoms):
1059 +
                        for alpha in range(nalpha):
1060 +
                            for beta in range(nbcc):
1061 +
                                self.C[alpha][beta] += R[j][alpha] * T[k][beta] * np.multiply(
1062 +
                                    np.dot(self.dist[k], self.dist_x[j]), ESP.e_field_at_atom[0][j])
1063 +
                                self.C[alpha][beta] += R[j][alpha] * T[k][beta] * np.multiply(
1064 +
                                    np.dot(self.dist[k], self.dist_y[j]), ESP.e_field_at_atom[1][j])
1065 +
                                self.C[alpha][beta] += R[j][alpha] * T[k][beta] * np.multiply(
1066 +
                                    np.dot(self.dist[k], self.dist_z[j]), ESP.e_field_at_atom[2][j])
1067 +
1068 +
        # Restraints for polarizaton
1069 +
        #if hasattr(self, 'wrst1'):
1070 +
        #    for alpha in range(nalpha):
1071 +
        #        self.D[alpha][alpha] += mol2.group_pop[alpha] * self.wrst1 * self.dipscale / np.sqrt(
1072 +
        #           np.square(self.qd[nbcc + alpha]) + 0.01)
1073 +
1074 +
        # No implementation of bcc restraints.
1075 +
1076 +
1077 +
        #Combine all matrices
1078 +
        self.X = np.concatenate((np.concatenate((self.A, self.B), axis=1), np.concatenate((self.C, self.D), axis=1)),
1079 +
                                axis=0)
1080 +
1081 +
1082 +
        """
1083 +
        X12= np.zeros((len(self.X), len(self._molecule.same_polarization_atoms)))
1084 +
        X22 = np.zeros((len(self._molecule.same_polarization_atoms),len(self._molecule.same_polarization_atoms)))
1085 +
        self.X = np.concatenate((np.concatenate((self.X, X12), axis=1), np.concatenate((X12.transpose, X22), axis=1)),axis=0)
1086 +
        
1087 +
        
1088 +
        for j, atoms in enumerate(self._molecule.same_polarization_atoms):
1089 +
            self.X[5 * self.natoms + j] = 0.0
1090 +
1091 +
        """
1092 +
    def build_vector_Y_BCC(self):
1093 +
        """
1094 +
        Creates vector Y for the BCC Pol method
1095 +
1096 +
        :param mol2: molecule object
1097 +
1098 +
        :return:
1099 +
        """
1100 +
1101 +
        nbcc = self._molecule._trainingset._nbccs
1102 +
        nalpha = self._molecule._trainingset._nalpha
1103 +
        T = self._molecule.T
1104 +
        R = self._molecule.R
1105 +
1106 +
        # Vector belonging to the BCCs
1107 +
        self.Y1 = np.zeros(nbcc)
1108 +
        if self._molecule._mode != 'alpha':
1109 +
            for ESP in [self.baseESP] + self.polESPs:
1110 +
                esp_values = ESP.esp_values.to('elementary_charge / angstrom').magnitude
1111 +
                for beta in range(nbcc):
1112 +
                    for k in range(self.natoms):
1113 +
                        self.Y1[beta] += T[k][beta] * np.dot(esp_values, self.dist[k])
1114 +
        else:
1115 +
            self.Y1 = self.qbond
1116 +
1117 +
        # Vector belonging to the polarizabilities
1118 +
        self.Y2 = np.zeros(nalpha)
1119 +
        if self._molecule._mode != 'bcconly':
1120 +
            for ESP in [self.baseESP] + self.polESPs:
1121 +
                esp_values = ESP.esp_values.to('elementary_charge / angstrom').magnitude
1122 +
                for beta in range(nalpha):
1123 +
                    for k in range(self.natoms):
1124 +
                        self.Y2[beta] += R[k][beta] * np.multiply(np.dot(esp_values, self.dist_x[k]), ESP.e_field_at_atom[0][k])
1125 +
                        self.Y2[beta] += R[k][beta] * np.multiply(np.dot(esp_values, self.dist_y[k]), ESP.e_field_at_atom[1][k])
1126 +
                        self.Y2[beta] += R[k][beta] * np.multiply(np.dot(esp_values, self.dist_z[k]), ESP.e_field_at_atom[2][k])
1127 +
1128 +
1129 +
        self.Y = np.concatenate((self.Y1, self.Y2))
1130 +
1131 +
1132 +
    # noinspection PyProtectedMember
1133 +
    def build_vector_B(self):
1134 +
        """
1135 +
        Creates the Vector B for the charge fitting.
1136 +
        :return:
1137 +
        """
1138 +
        # Determine size of the vector for this molecule
1139 +
        # every atom is one line
1140 +
        # one line to restrain the overall charge of the molecule
1141 +
        # one line for every pair of equivalent atoms
1142 +
        self._lines_in_A = self.natoms + 1 + len(self._molecule.chemical_eq_atoms)
1143 +
        self.B = np.zeros(self._lines_in_A)
1144 +
        self.esp_values = self.baseESP.esp_values.to('elementary_charge / angstrom').magnitude
1145 +
1146 +
        for k in range(self.natoms):
1147 +
            self.B[k] = np.dot(self.esp_values, self.dist[k])
1148 +
        self.B[self.natoms] = self._molecule._charge
1149 +
        for polESP in self.polESPs:
1150 +
            esp_values = polESP.esp_values.to('elementary_charge / angstrom').magnitude
1151 +
            for k in range(self.natoms):
1152 +
                self.B[k] += np.dot(esp_values, self.dist[k])
1153 +
            self.B[self.natoms] = self._molecule._charge
1154 +
1155 +
        self.B = self.B / (len(self.polESPs) + 1)
1156 +
        for k, eq_atoms in enumerate(self._molecule.chemical_eq_atoms):
1157 +
            if eq_atoms[1] > 0:
1158 +
                self.B[self.natoms + 1 + k] = 0.0
1159 +
1160 +
    def build_vector_C(self):
1161 +
        """
1162 +
        Vector C corresponds to matrix D and is only for pur polarization fitting.
1163 +
1164 +
        :return:
1165 +
        """
1166 +
        self.C = np.zeros(self.Dlines)
1167 +
        self.esp_values = self.baseESP.esp_values.to('elementary_charge / angstrom').magnitude
1168 +
        for k in range(self.natoms):
1169 +
            self.C[k] = np.multiply(np.dot(self.esp_values, self.dist_x[k]), self.e_field_at_atom[0][k])
1170 +
            self.C[k + self.natoms] = np.multiply(np.dot(self.esp_values, self.dist_y[k]), self.e_field_at_atom[1][k])
1171 +
            self.C[k + self.natoms * 2] = np.multiply(np.dot(self.esp_values, self.dist_z[k]),
1172 +
                                                      self.e_field_at_atom[2][k])
1173 +
1174 +
        for polESP in self.polESPs:
1175 +
            esp_values = polESP.esp_values.to('elementary_charge / angstrom').magnitude
1176 +
            for k in range(self.natoms):
1177 +
                self.C[k] += np.multiply(np.dot(esp_values, self.dist_x[k]), polESP.e_field_at_atom[0][k])
1178 +
                self.C[k + self.natoms] += np.multiply(np.dot(esp_values, self.dist_y[k]), polESP.e_field_at_atom[1][k])
1179 +
                self.C[k + self.natoms * 2] += np.multiply(np.dot(esp_values, self.dist_z[k]),
1180 +
                                                           polESP.e_field_at_atom[2][k])
1181 +
1182 +
        self.C = self.C / (len(self.polESPs) + 1)
1183 +
1184 +
        #for j, atoms in enumerate(self._molecule.same_polarization_atoms):
1185 +
        #    self.C[5 * self.natoms + j] = 0.0
1186 +
1187 +
    def build_vector_Y(self):
1188 +
        """
1189 +
        Creates the Vector Y for the RESP-Pol Method.
1190 +
        :return:
1191 +
        """
1192 +
        self.build_vector_B()
1193 +
        self.build_vector_C()
1194 +
        self.Y = np.concatenate((self.B, self.C))
1195 +
1196 +
    def get_electric_field(self):
1197 +
        self.baseESP.get_electric_field()
1198 +
        self.e_field_at_atom = self.baseESP.e_field_at_atom
1199 +
        for polESP in self.polESPs:
1200 +
            polESP.get_electric_field()
1201 +
1202 +
    def optimize_charges_alpha(self):
1203 +
        """
1204 +
        Builds the necessary matrix and vector and performs a charge and polarizabilities optimization for this 1 conformation.
1205 +
        :return:
1206 +
        """
1207 +
        self.build_matrix_X()
1208 +
        self.build_vector_Y()
1209 +
        self.q_alpha = np.linalg.solve(self.X, self.Y)
1210 +
1211 +
    def optimize_charges_alpha_bcc(self):
1212 +
        """
1213 +
        Builds the necessary matrix and vector and performs a charge and polarizabilities optimization for this 1 conformation.
1214 +
        :return:
1215 +
        """
1216 +
        self.build_matrix_X_BCC()
1217 +
        self.build_vector_Y_BCC()
1218 +
        # Check if bccs or alphas are not in the training set and set them to zero
1219 +
        for i,row in enumerate(self.X):
1220 +
            if all(row == 0.0):
1221 +
                self.X[i][i]=1
1222 +
1223 +
1224 +
        self.q_alpha = np.linalg.solve(self.X, self.Y)
1225 +
1226 +
1227 +
    def get_distances(self):
1228 +
        self.get_grid_coord()
1229 +
1230 +
        # Distances between atoms and ESP points
1231 +
        self.dist = np.zeros((self.natoms, self.npoints))
1232 +
        self.dist_3 = np.zeros((self.natoms, self.npoints))
1233 +
        self.dist_x = np.zeros((self.natoms, self.npoints))
1234 +
        self.dist_y = np.zeros((self.natoms, self.npoints))
1235 +
        self.dist_z = np.zeros((self.natoms, self.npoints))
1236 +
1237 +
        self.dist = 1. / distance.cdist(self.atom_positions_angstrom, self.grid_coord_angstrom)
1238 +
        self.dist_3 = np.power(self.dist, 3)  # maybe free afterwards
1239 +
        self.dist_x = -np.multiply(
1240 +
            np.subtract.outer(np.transpose(self.atom_positions_angstrom)[0], np.transpose(self.grid_coord_angstrom)[0]),
1241 +
            self.dist_3)
1242 +
        # self.dist_x2=np.multiply(np.transpose(np.subtract.outer(np.transpose(self.grid_coord)[0],np.transpose(self.atom_positions)[0])),self.dist_3)
1243 +
        self.dist_y = -np.multiply(
1244 +
            np.subtract.outer(np.transpose(self.atom_positions_angstrom)[1], np.transpose(self.grid_coord_angstrom)[1]),
1245 +
            self.dist_3)
1246 +
        self.dist_z = -np.multiply(
1247 +
            np.subtract.outer(np.transpose(self.atom_positions_angstrom)[2], np.transpose(self.grid_coord_angstrom)[2]),
1248 +
            self.dist_3)
1249 +
        del self.dist_3
1250 +
1251 +
        # Distances between atoms and atoms
1252 +
        self.diatomic_dist = np.zeros((self.natoms, self.natoms))
1253 +
        self.diatomic_dist_3 = np.zeros((self.natoms, self.natoms))
1254 +
        self.diatomic_dist_5 = np.zeros((self.natoms, self.natoms))
1255 +
        self.diatomic_dist_x = np.zeros((self.natoms, self.natoms))
1256 +
        self.diatomic_dist_y = np.zeros((self.natoms, self.natoms))
1257 +
        self.diatomic_dist_z = np.zeros((self.natoms, self.natoms))
1258 +
        self.diatomic_distb_x = np.zeros((self.natoms, self.natoms))
1259 +
        self.diatomic_distb_y = np.zeros((self.natoms, self.natoms))
1260 +
        self.diatomic_distb_z = np.zeros((self.natoms, self.natoms))
1261 +
1262 +
        self.diatomic_dist = distance.cdist(self.atom_positions_angstrom, self.atom_positions_angstrom)
1263 +
        di = np.diag_indices(self.natoms)
1264 +
        self.diatomic_dist[di] = 1.0E10
1265 +
        # self.adist=np.fill_diagonal(self.adist,1.0)
1266 +
        self.diatomic_dist = 1. / self.diatomic_dist
1267 +
        self.diatomic_dist_3 = np.power(self.diatomic_dist, 3)
1268 +
        self.diatomic_dist_5 = np.power(self.diatomic_dist, 5)
1269 +
        self.diatomic_dist[di] = 0.0
1270 +
        self.diatomic_dist_x = np.multiply(np.subtract.outer(np.transpose(self.atom_positions_angstrom)[0],
1271 +
                                                             np.transpose(self.atom_positions_angstrom)[0]),
1272 +
                                           self.diatomic_dist_3)  # X distance between two atoms divided by the dist^3
1273 +
        self.diatomic_dist_y = np.multiply(np.subtract.outer(np.transpose(self.atom_positions_angstrom)[1],
1274 +
                                                             np.transpose(self.atom_positions_angstrom)[1]),
1275 +
                                           self.diatomic_dist_3)
1276 +
        self.diatomic_dist_z = np.multiply(np.subtract.outer(np.transpose(self.atom_positions_angstrom)[2],
1277 +
                                                             np.transpose(self.atom_positions_angstrom)[2]),
1278 +
                                           self.diatomic_dist_3)
1279 +
        self.diatomic_distb_x = np.subtract.outer(np.transpose(self.atom_positions_angstrom)[0],
1280 +
                                                  np.transpose(self.atom_positions_angstrom)[
1281 +
                                                      0])  # X distances between two atoms
1282 +
        self.diatomic_distb_y = np.subtract.outer(np.transpose(self.atom_positions_angstrom)[1],
1283 +
                                                  np.transpose(self.atom_positions_angstrom)[1])
1284 +
        self.diatomic_distb_z = np.subtract.outer(np.transpose(self.atom_positions_angstrom)[2],
1285 +
                                                  np.transpose(self.atom_positions_angstrom)[2])
1286 +
1287 +
    def delete_distances(self):
1288 +
        """Deletes the all calculated distances to free memory."""
1289 +
        del self.dist
1290 +
        del self.dist_x
1291 +
        del self.dist_y
1292 +
        del self.dist_z
1293 +
1294 +
        del self.diatomic_dist
1295 +
        del self.diatomic_dist_3
1296 +
        del self.diatomic_dist_5
1297 +
        del self.diatomic_dist_x
1298 +
        del self.diatomic_dist_y
1299 +
        del self.diatomic_dist_z
1300 +
        del self.diatomic_distb_x
1301 +
        del self.diatomic_distb_y
1302 +
        del self.diatomic_distb_z
1303 +
1304 +
1305 +
# =============================================================================================
1306 +
# ESPGRID
1307 +
# =============================================================================================
1308 +
1309 +
# noinspection PyTypeChecker
1310 +
class ESPGRID:
1311 +
    """
1312 +
1313 +
    """
1314 +
1315 +
    def define_grid(self, *args, **kwargs):
1316 +
        for ele in args:
1317 +
            if 'gesp' in ele:
1318 +
                self.gridtype = 'gesp'
1319 +
            if 'espf' in ele:
1320 +
                self.gridtype = 'respyte'
1321 +
            if 'grid.dat' in ele:
1322 +
                self.gridtype = 'psi4'
1323 +
1324 +
    def set_ext_e_field(self, vector):
1325 +
        self._ext_e_field = vector
1326 +
1327 +
    def load_grid(self, *args):
1328 +
        if self.gridtype == 'gesp':
1329 +
            self.name = args[0].rstrip('.gesp')
1330 +
            f = open(args[0], 'r')
1331 +
            lines = f.readlines()
1332 +
            f.close()
1333 +
            for i, line in enumerate(lines):
1334 +
                if 'ATOMIC' in line and 'COORDINATES' in line:
1335 +
                    self.natoms = int(line.strip('\n').split()[-1])
1336 +
                    for j in range(self.natoms):
1337 +
                        entry = lines[i + 1 + j].replace('D', 'E').split()
1338 +
                        self.atoms.append(entry[0])
1339 +
                        self.atom_positions.append(Q_([float(entry[k]) for k in range(1, 4, 1)], 'bohr'))
1340 +
                if 'GRID' in line:
1341 +
                    self.ngrid = int(line.strip('\n').split()[-1])
1342 +
                    grid = lines[i + 1:i + 1 + self.ngrid]
1343 +
                    break
1344 +
            # noinspection PyUnboundLocalVariable
1345 +
            for i, line in enumerate(grid):
1346 +
                grid[i] = [float(ele) for ele in line.replace('D', 'E').split()]
1347 +
1348 +
            self.positions = Q_(np.array(grid)[:, 1:4], 'bohr')
1349 +
            self.esp_values = Q_(np.array(grid)[:, 0], 'elementary_charge / bohr')
1350 +
        elif self.gridtype == 'respyte':
1351 +
            self.name = args[0].rstrip('.espf')
1352 +
            f = open(args[0], 'r')
1353 +
            lines = f.readlines()
1354 +
            f.close
1355 +
            ndata = int(len(lines) / 2) if len(lines) % 2 == 0 else int((len(lines) - 1) / 2)
1356 +
            grid = np.zeros((ndata, 4))
1357 +
            for i in range(ndata):
1358 +
                grid[i] = [float(ele) for ele in lines[2 * i].split()]
1359 +
            self.positions = Q_(np.array(grid)[:, 0:3], 'angstrom')
1360 +
            self.esp_values = Q_(np.array(grid)[:, 3], 'elementary_charge / bohr')
1361 +
        elif self.gridtype == 'psi4':
1362 +
            for ele in args:
1363 +
                if "grid.dat" in ele:
1364 +
                    gridfile = ele
1365 +
                elif 'esp.dat' in ele:
1366 +
                    espfile = ele
1367 +
            self.name = 'esp'
1368 +
            np.loadtxt(espfile)
1369 +
            self.positions = Q_(np.loadtxt(gridfile), 'angstrom')
1370 +
            self.esp_values = Q_(np.loadtxt(espfile), 'elementary_charge / bohr')
1371 +
1372 +
    def get_electric_field(self, ):
1373 +
        """
1374 +
        Calculates the electric field at every atomic positions.
1375 +
        :return:
1376 +
        """
1377 +
        alpha = self._conformer.q_alpha[self._conformer._lines_in_A:self._conformer._lines_in_A + 3*self._conformer.natoms]
1378 +
        alpha[np.where(alpha == 0.0)] += 10E-10
1379 +
1380 +
        # Load permanent charges for BCC method
1381 +
        # For all other methods this is set to 0.0 STILL HAVE to implment
1382 +
        try:
1383 +
            self._conformer.q_am1
1384 +
        except Exception:
1385 +
            log.warning('I do not have AM1-type charges')
1386 +
            self._conformer.q_am1 = np.zeros(self._conformer.natoms)
1387 +
1388 +
        for j in range(self._conformer.natoms):
1389 +
            self.e_field_at_atom[0][j] = np.dot(
1390 +
                np.multiply(self._conformer.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1391 +
                self._conformer.diatomic_dist_x[j]) + np.dot(
1392 +
                np.multiply(self._conformer.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1393 +
                self._conformer.diatomic_dist_x[j]) + self._ext_e_field[0].to(
1394 +
                'elementary_charge / angstrom / angstrom').magnitude
1395 +
            self.e_field_at_atom[1][j] = np.dot(
1396 +
                np.multiply(self._conformer.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1397 +
                self._conformer.diatomic_dist_y[j]) + np.dot(
1398 +
                np.multiply(self._conformer.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1399 +
                self._conformer.diatomic_dist_y[j]) + self._ext_e_field[1].to(
1400 +
                'elementary_charge / angstrom / angstrom').magnitude
1401 +
            self.e_field_at_atom[2][j] = np.dot(
1402 +
                np.multiply(self._conformer.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1403 +
                self._conformer.diatomic_dist_z[j]) + np.dot(
1404 +
                np.multiply(self._conformer.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1405 +
                self._conformer.diatomic_dist_z[j]) + self._ext_e_field[2].to(
1406 +
                'elementary_charge / angstrom / angstrom').magnitude
1407 +
1408 +
        self.e = self.e_field_at_atom.flatten()
1409 +
        #print(self.e)
1410 +
        if self._conformer._molecule._SCF and self._conformer._molecule.step > 0:
1411 +
            if not hasattr(self, 'Bdip') or self._conformer._molecule._thole:
1412 +
                if self._conformer._molecule._thole:
1413 +
                    thole_param=1.368711 # Change back to 1.368711
1414 +
                    #thole_param = 0.390
1415 +
                    dipole_tmp = np.where(alpha < 0.0, -alpha, alpha)
1416 +
                    thole_v = np.multiply(self._conformer.diatomic_dist, np.float_power(
1417 +
                        np.multiply(dipole_tmp[:self._conformer.natoms, None], dipole_tmp[:self._conformer.natoms]),
1418 +
                        1. / 6))
1419 +
                    di = np.diag_indices(self._conformer.natoms)
1420 +
                    thole_v[di] = 1.0
1421 +
                    thole_v = 1. / thole_v
1422 +
                    thole_v[di] = 0.0
1423 +
1424 +
                    # Exponential thole
1425 +
                    thole_fe = np.ones((self._conformer.natoms, self._conformer.natoms))
1426 +
                    thole_ft = np.ones((self._conformer.natoms, self._conformer.natoms))
1427 +
                    thole_fe -= np.exp(np.multiply(thole_param, np.power(-thole_v, 3)))
1428 +
                    thole_ft -= np.multiply(np.multiply(thole_param, np.power(thole_v, 3)) + 1.,
1429 +
                                            np.exp(np.multiply(thole_param, np.power(-thole_v, 3))))
1430 +
                    # 1.5 was found in the OpenMM code. Not sure whuy it is there
1431 +
1432 +
                    # In original thole these lines should not be here
1433 +
                    # thole_ft = np.multiply(thole_ft, self._conformer.scale)
1434 +
                    # thole_fe = np.multiply(thole_fe, self._conformer.scale)
1435 +
                    # Linear thole
1436 +
                    # thole_fe=np.zeros((self._conformer.natoms,self._conformer.natoms))
1437 +
                    # thole_fe=np.zeros((self._conformer.natoms,self._conformer.natoms))
1438 +
                    # thole_fe=np.where(thole_v>1.0,1.0,4*np.power(thole_v,3)-3*np.power(thole_v,4))
1439 +
                    # thole_ft=np.where(thole_v>1.0,1.0,np.power(thole_v,4))
1440 +
                else:
1441 +
                    try:
1442 +
                        thole_ft = self._conformer._molecule.scale_scf
1443 +
                        thole_fe = self._conformer._molecule.scale_scf
1444 +
                    except Exception:
1445 +
1446 +
                        thole_ft = self._conformer._molecule.scale
1447 +
                        thole_fe = self._conformer._molecule.scale
1448 +
                    else:
1449 +
                        print('Using different set of scaling for SCF interactions')
1450 +
                        log.info('Using different set of scaling for SCF interactions')
1451 +
                Bdip11 = np.add(np.multiply(thole_fe, self._conformer.diatomic_dist_3), np.multiply(thole_ft,
1452 +
                                                                                                    -3 * np.multiply(
1453 +
                                                                                                        np.multiply(
1454 +
                                                                                                            self._conformer.diatomic_distb_x,
1455 +
                                                                                                            self._conformer.diatomic_distb_x),
1456 +
                                                                                                        self._conformer.diatomic_dist_5)))
1457 +
                Bdip22 = np.add(np.multiply(thole_fe, self._conformer.diatomic_dist_3), np.multiply(thole_ft,
1458 +
                                                                                                    -3 * np.multiply(
1459 +
                                                                                                        np.multiply(
1460 +
                                                                                                            self._conformer.diatomic_distb_y,
1461 +
                                                                                                            self._conformer.diatomic_distb_y),
1462 +
                                                                                                        self._conformer.diatomic_dist_5)))
1463 +
                Bdip33 = np.add(np.multiply(thole_fe, self._conformer.diatomic_dist_3), np.multiply(thole_ft,
1464 +
                                                                                                    -3 * np.multiply(
1465 +
                                                                                                        np.multiply(
1466 +
                                                                                                            self._conformer.diatomic_distb_z,
1467 +
                                                                                                            self._conformer.diatomic_distb_z),
1468 +
                                                                                                        self._conformer.diatomic_dist_5)))
1469 +
                Bdip12 = np.multiply(thole_ft,
1470 +
                                     -3 * np.multiply(np.multiply(self._conformer.diatomic_distb_x,
1471 +
                                                                  self._conformer.diatomic_distb_y),
1472 +
                                                      self._conformer.diatomic_dist_5))
1473 +
                Bdip13 = np.multiply(thole_ft,
1474 +
                                     -3 * np.multiply(np.multiply(self._conformer.diatomic_distb_x,
1475 +
                                                                  self._conformer.diatomic_distb_z),
1476 +
                                                      self._conformer.diatomic_dist_5))
1477 +
                Bdip23 = np.multiply(thole_ft,
1478 +
                                     -3 * np.multiply(np.multiply(self._conformer.diatomic_distb_y,
1479 +
                                                                  self._conformer.diatomic_distb_z),
1480 +
                                                      self._conformer.diatomic_dist_5))
1481 +
                Bdip = np.concatenate((np.concatenate((Bdip11, Bdip12, Bdip13), axis=1),
1482 +
                                       np.concatenate((Bdip12, Bdip22, Bdip23), axis=1),
1483 +
                                       np.concatenate((Bdip13, Bdip23, Bdip33), axis=1)),
1484 +
                                      axis=0)
1485 +
1486 +
            for j in range(self._conformer.natoms):
1487 +
                for k in range(3):
1488 +
                    for l in range(3):
1489 +
                        Bdip[k * self._conformer.natoms + j][l * self._conformer.natoms + j] = 0.0
1490 +
1491 +
            for j in range(3 * self._conformer.natoms):
1492 +
                Bdip[j][j] = 1. / alpha[j]
1493 +
            dipole_scf = np.linalg.solve(Bdip, self.e)
1494 +
            self.e = np.divide(dipole_scf, alpha)
1495 +
            print(self.e)
1496 +
        self.e_field_at_atom[0] = 1.0 * self.e[:self._conformer.natoms]
1497 +
        self.e_field_at_atom[1] = 1.0 * self.e[self._conformer.natoms:2 * self._conformer.natoms]
1498 +
        self.e_field_at_atom[2] = 1.0 * self.e[2 * self._conformer.natoms:3 * self._conformer.natoms]
1499 +
1500 +
    # Analysation stuff
1501 +
    def calc_esp_q_alpha(self, q_alpha):
1502 +
        """
1503 +
        Calculates the ESP for a given set of point charges and polarizabilities.
1504 +
1505 +
        :param qd: list of float
1506 +
            Point charges and polarizabilities
1507 +
        :return:
1508 +
1509 +
        Calculates the ESP on every point of the initially read GESP file
1510 +
        """
1511 +
        self.q_pot = np.zeros(len(self.esp_values))
1512 +
        for i in range(len(self.esp_values)):
1513 +
            if self._conformer._molecule._mode == 'q':
1514 +
                self.q_pot[i] = np.dot(q_alpha[:self._conformer.natoms], np.transpose(self._conformer.dist)[i])
1515 +
            else:
1516 +
                dipole = q_alpha[self._conformer._lines_in_A:]
1517 +
                # self.dipole[self.dipole<0.0]=0.0
1518 +
                try:
1519 +
                    # e_dip=np.dot(self.dipole_scf[:self.natoms],np.transpose(self.dist_x)[i])+np.dot(self.dipole_scf[self.natoms:2*self.natoms],np.transpose(self.dist_y)[i])+np.dot(self.dipole_scf[2*self.natoms:3*self.natoms],np.transpose(self.dist_z)[i])
1520 +
                    e_dip = np.dot(np.multiply(dipole[:self.natoms], self.e_field_at_atom[0]),
1521 +
                                   np.transpose(self._conformer.dist_x)[i]) + np.dot(
1522 +
                        np.multiply(dipole[self.natoms:2 * self.natoms], self.e_field_at_atom[1]),
1523 +
                        np.transpose(self._conformer.dist_y)[i]) + np.dot(
1524 +
                        np.multiply(dipole[2 * self.natoms:3 * self.natoms], self.e_field_at_atom[2]),
1525 +
                        np.transpose(self._conformer.dist_z)[i])
1526 +
                    self.q_pot[i] = np.dot(q_alpha[:self.natoms], np.transpose(self._conformer.dist)[i]) + e_dip
1527 +
                except:
1528 +
                    self.q_pot[i] = np.dot(q_alpha[:self.natoms], np.transpose(self._conformer.dist)[i])
1529 +
                    # print('DEBUG no electric field ')
1530 +
            # if self.mode=='d':
1531 +
            #    self.dipole=qd
1532 +
            #    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])
1533 +
            #    self.q_pot[i]=e_dip
1534 +
1535 +
    def sub_esp_q_alpha(self, q_alpha):
1536 +
        """
1537 +
        Subtracts the ESP create by a set of point charges and polarizabilities.
1538 +
1539 +
        :param qd: list of float
1540 +
            Point charges and polarizabilities
1541 +
        :return:
1542 +
1543 +
        """
1544 +
        if self._conformer._molecule._mode != 'alpha':  # Change that in bcc that this line is unnecessary
1545 +
            self.calc_esp_q_alpha(q_alpha)
1546 +
            self.esp_values = np.subtract(self.esp_values.to('elementary_charge / angstrom').magnitude, self.q_pot)
1547 +
1548 +
    def calc_sse(self, q_alpha):
1549 +
        """
1550 +
        Calculate the Squared Sum of Errors between the stored ESP and a ESP created from qd.
1551 +
        :param qd: list of float
1552 +
            Point charges and polarizabilities
1553 +
        :return:
1554 +
        """
1555 +
        self.calc_esp_q_alpha(q_alpha)
1556 +
        self.sse = np.square(self.esp_values.to('elementary_charge / angstrom').magnitude - self.q_pot).sum()
1557 +
1558 +
    def write_res_esp(self, q_alpha):
1559 +
        """
1560 +
        NOT FINISHED YET!!!
1561 +
1562 +
        Writes the residual ESP to a file.
1563 +
1564 +
        :param qd: list of float
1565 +
            Point charges and polarizabilities
1566 +
        :return:
1567 +
        """
1568 +
        self.calc_esp_q_alpha(q_alpha)
1569 +
        res_pot = self.esp_values - self.q_pot
1570 +
        f = open(self.name + '.rgesp', 'w')
1571 +
        f.write(' ESP FILE - ATOMIC UNITS\n')
1572 +
        f.write(' CHARGE =  {0} - MULTIPLICITY =   1\n'.format(self._conformer._molecule._charge))
1573 +
        f.write(' ATOMIC COORDINATES AND ESP CHARGES. #ATOMS =     {} \n'.format(np.sum(self.natoms)))
1574 +
        for i in range(self.conformer._natoms):
1575 +
            f.write(
1576 +
                ' {} {} {} {} {}\n'.format(self.atoms[i], self.atomcrd[i][0], self.atomcrd[i][1], self.atomcrd[i][2],
1577 +
                                           self._conformer._moleccule.q_alpha[i]))
1578 +
        f.write(' ESP VALUES AND GRID POINT COORDINATES. #POINTS =   {}\n'.format(np.sum(self.npoints)))
1579 +
        for i in range(len(self.esp_values)):
1580 +
            f.write(' {} {} {} {}\n'.format(res_pot[i], self.crd[i][0], self.crd[i][1], self.crd[i][2]))
1581 +
        f.close()
1582 +
1583 +
# =============================================================================================
1584 +
# BCCUnpolESP
1585 +
# =============================================================================================
1586 +
1587 +
class BCCUnpolESP(ESPGRID):
1588 +
    """
1589 +
1590 +
    """
1591 +
1592 +
    def __init__(self, *args, conformer=None):
1593 +
        # Decide if we have a Gaussian grid or a psi4 grid
1594 +
        self.gridtype = None
1595 +
        self.natoms = -1
1596 +
        self.atoms = []
1597 +
        self.atom_positions = []
1598 +
        self.define_grid(*args)
1599 +
        self.esp_values = None
1600 +
        self.positions = None
1601 +
        self._conformer = conformer
1602 +
1603 +
        # External e-field is 0 in all directions
1604 +
        vector = Q_([0, 0, 0], 'elementary_charge / bohr / bohr')
1605 +
        self.set_ext_e_field(vector)
1606 +
1607 +
        self.load_grid(*args)
1608 +
1609 +
        self.e_field_at_atom = np.zeros((3, self._conformer.natoms))
1610 +
1611 +
1612 +
# =============================================================================================
1613 +
# BCCPolESP
1614 +
# =============================================================================================
1615 +
1616 +
class BCCPolESP(ESPGRID):
1617 +
    """
1618 +
1619 +
    """
1620 +
1621 +
    def __init__(self, *args, conformer=None, e_field=Q_([0.0, 0.0, 0.0], 'elementary_charge / bohr / bohr')):
1622 +
        # Decide if we have a Gaussian grid or a psi4 grid
1623 +
        self.gridtype = None
1624 +
        self.natoms = -1
1625 +
        self.atoms = []
1626 +
        self.atom_positions = []
1627 +
        self.define_grid(*args)
1628 +
        self.esp_values = None
1629 +
        self.positions = None
1630 +
        self._conformer = conformer
1631 +
1632 +
        # Set external e-field
1633 +
        self.set_ext_e_field(e_field)
1634 +
1635 +
        self.load_grid(*args)
1636 +
1637 +
        self.e_field_at_atom = np.zeros((3, self._conformer.natoms))
1638 +
1639 +
1640 +
if __name__ == '__main__':
1641 +
    """
1642 +
    datei = os.path.join(ROOT_DIR_PATH, 'resppol/tmp/phenol/conf0/mp2_0.mol2')
1643 +
    test = TrainingSet(scf_scaleparameters=[0.0, 0.0, 0.5])
1644 +
    test.add_molecule(datei)
1645 +
    test.molecules[0].add_conformer_from_mol2(datei)
1646 +
    print(test.molecules[0].same_polarization_atoms)
1647 +
    print(test.molecules[0].scale)
1648 +
    print(test.molecules[0].scale_scf)
1649 +
    espfile = '/home/michael/resppol/resppol/tmp/phenol/conf0/molecule0.gesp'
1650 +
    test.molecules[0].conformers[0].add_baseESP(espfile)
1651 +
    datei = os.path.join(ROOT_DIR_PATH, 'resppol/tmp/butanol/conf0/mp2_0.mol2')
1652 +
    test.add_molecule(datei)
1653 +
    test.molecules[1].add_conformer_from_mol2(datei)
1654 +
    espfile = '/home/michael/resppol/resppol/tmp/butanol/conf0/molecule0.gesp'
1655 +
    test.molecules[1].conformers[0].add_baseESP(espfile)
1656 +
    test.optimize_charges()
1657 +
    test.molecules[0].conformers[0].build_matrix_X()
1658 +
    for molecule in test.molecules:
1659 +
        print(molecule.q)
1660 +
1661 +
    print('FINISH')
1662 +
    
1663 +
    datei = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test2.mol2')
1664 +
    test = TrainingSet(mode='q_alpha',SCF= True,scf_scaleparameters=[1,1,1], scaleparameters=[1,1,1])
1665 +
    #test = TrainingSet(mode='q_alpha')
1666 +
    test.add_molecule(datei)
1667 +
    test.molecules[0].add_conformer_from_mol2(datei)
1668 +
    #test.molecules[0].add_conformer_from_mol2(datei)
1669 +
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3.gesp')
1670 +
    test.molecules[0].conformers[0].add_baseESP(espfile)
1671 +
    #test.molecules[0].conformers[1].add_baseESP(espfile)
1672 +
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3_Z+.gesp')
1673 +
    test.molecules[0].conformers[0].add_polESP(espfile, e_field=Q_([0.0, 0.0, 1.], 'elementary_charge / bohr / bohr'))
1674 +
    #test.molecules[0].conformers[1].add_polESP(espfile, e_field=Q_([0.0, 0.0, 0.8],'elementary_charge / bohr / bohr'))
1675 +
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3_Z-.gesp')
1676 +
    test.molecules[0].conformers[0].add_polESP(espfile, e_field=Q_([0.0, 0.0, -1.], 'elementary_charge / bohr / bohr') )
1677 +
    #test.molecules[0].conformers[1].add_polESP(espfile, e_field=Q_([0.0, 0.0, -0.8], 'elementary_charge / bohr / bohr') )
1678 +
    #test.build_matrix_X()
1679 +
    #test.build_vector_Y()
1680 +
    test.build_matrix_X_BCC()
1681 +
    test.build_vector_Y_BCC()
1682 +
    test.optimize_bcc_alpha()
1683 +
    print(test.molecules[0].conformers[0].q_alpha)
1684 +
    #test.optimize_charges_alpha()
1685 +
    #print(test.q_alpha)
1686 +
    print(test.molecules[0].conformers[0].baseESP.e_field_at_atom)
1687 +
    print(test.molecules[0].conformers[0].polESPs[0].e_field_at_atom)
1688 +
    """
1689 +
1690 +
    datei = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test2.mol2')
1691 +
    test = TrainingSet(mode='q_alpha',SCF= True, thole = True)
1692 +
    test.add_molecule(datei)
1693 +
    test.molecules[0].add_conformer_from_mol2(datei)
1694 +
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3.gesp')
1695 +
    test.molecules[0].conformers[0].add_baseESP(espfile)
1696 +
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3_Z+.gesp')
1697 +
    test.molecules[0].conformers[0].add_polESP(espfile, e_field=Q_([0.0, 0.0, 1.0], 'elementary_charge / bohr / bohr'))
1698 +
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3_Z-.gesp')
1699 +
    test.molecules[0].conformers[0].add_polESP(espfile, e_field=Q_([0.0, 0.0, -1.0], 'elementary_charge / bohr / bohr') )
1700 +
    test.optimize_charges_alpha()
1701 +
    test.molecules[0].conformers[0].baseESP.calc_esp_q_alpha(test.q_alpha)
1702 +
    test.molecules[0].conformers[0].baseESP.calc_sse(test.q_alpha)
1703 +
    test.molecules[0].conformers[0].baseESP.sub_esp_q_alpha(test.q_alpha)
1704 +
    print(test.molecules[0].conformers[0].baseESP.q_pot)
1705 +
1706 +
    print(test.molecules[0].conformers[0].q_alpha)
Files Coverage
resppol 49.91%
Project Totals (3 files) 49.91%
254.2
TRAVIS_OS_NAME=osx
254.1
TRAVIS_OS_NAME=osx
254.4
TRAVIS_PYTHON_VERSION=3.7
TRAVIS_OS_NAME=linux
254.3
TRAVIS_PYTHON_VERSION=3.6
TRAVIS_OS_NAME=linux
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