1
# !/usr/bin/env python
2

3
# =============================================================================================
4
# MODULE DOCSTRING
5
# =============================================================================================
6 100
"""
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 100
import logging as log
19 100
import numpy as np
20 100
try:
21 100
	from openeye import oechem
22 0
except:
23 0
	print('Could not import openeye')
24 100
try:
25 100
	import openforcefield.topology as openff
26 100
	from openforcefield.typing.engines.smirnoff import ForceField
27 0
except:
28 0
	print("Could not import openforcefield")
29 100
from scipy.spatial import distance
30 100
import os
31
# Initialize units
32 100
from pint import UnitRegistry
33

34 100
ureg = UnitRegistry()
35 100
Q_ = ureg.Quantity
36 100
ureg.define('bohr = 0.52917721067 * angstrom')
37

38
# =============================================================================================
39
# GLOBAL PARAMETERS
40
# =============================================================================================
41 100
biological_elements = [1, 3, 5, 6, 7, 8, 9, 11, 12, 14, 15, 16, 17, 19, 20, 34, 35, 53]
42

43 100
ROOT_DIR_PATH = os.path.join(os.path.dirname(os.path.realpath(__file__)), '..')
44

45

46
# =============================================================================================
47
# PRIVATE SUBROUTINES
48
# =============================================================================================
49 100
def find_eq_atoms(mol1):
50
    """
51
    Finds all equivalent atoms in a molecule.
52
    :return
53
    Array of pairs of equivalent atoms.
54

55
    :parameter
56
    mol1: Openeye molecule object
57

58
    TODO Include rdkit support for this function
59
    """
60

61 100
    qmol = oechem.OEQMol()
62

63
    # build OEQMol from OEGRAPHMOLECULE
64 100
    oechem.OEBuildMDLQueryExpressions(qmol, mol1)
65 100
    ss2 = oechem.OESubSearch(qmol)
66 100
    oechem.OEPrepareSearch(mol1, ss2)
67

68
    # Store the equivalent atoms
69 100
    eq_atoms = [[] for i in range(mol1.NumAtoms())]
70

71
    # Goes through all matches and compares the atom indices.
72 100
    for count, match in enumerate(ss2.Match(mol1)):
73 100
        for ma in match.GetAtoms():
74
            # if 2 atoms are not the same atoms
75 100
            if ma.pattern.GetIdx() != ma.target.GetIdx():
76
                # and the pair is not stored yet
77 100
                if ma.target.GetIdx() not in eq_atoms[ma.pattern.GetIdx()]:
78
                    # save it to the array
79 100
                    eq_atoms[ma.pattern.GetIdx()].append(ma.target.GetIdx())
80

81
    # goes through the array and returns pairs of equivalent atoms
82 100
    sorted_eq_atoms = []
83 100
    for i, ele in enumerate(eq_atoms):
84 100
        for element in ele:
85
            # Making sure we have not already covered this pair of atoms
86 100
            if [element, i] not in sorted_eq_atoms:
87 100
                sorted_eq_atoms.append([i, element])
88 100
    tmp = []
89 100
    for k, ele1 in enumerate(sorted_eq_atoms):
90 100
        exclude = 0
91 100
        for ele2 in sorted_eq_atoms[:k]:
92 100
            if ele1[0] == ele2[0]:
93 100
                exclude = 1
94 100
        if exclude == 0:
95 100
            tmp.append(ele1)
96 100
    sorted_eq_atoms = tmp
97

98 100
    return sorted_eq_atoms
99

100

101

102
# =============================================================================================
103
# TrainingSet
104
# =============================================================================================
105

106 100
class TrainingSet():
107
    """
108
    The training set class is the top level class of the resspol program.
109

110
    It consist of multiple molecule instances and combines the optimization matrices and vectors across multiple molecules.
111
    """
112

113 100
    def __init__(self, mode='q', scaleparameters=None, scf_scaleparameters=None, SCF=False, thole=False, FF='resppol/data/test_data/BCCPOL.offxml'):
114
        """
115
        Initilize the class and sets the following parameters:
116

117
        :param mode:
118
        :param scaleparameters:
119
        :param scf_scaleparameters:
120
        :param SCF:
121
        :param thole:
122
        """
123 100
        self.molecules = list()
124 100
        self.B = np.zeros(0)
125 100
        self.A = np.zeros((0, 0))
126 100
        self.q = 0.0
127 100
        self.mode = mode
128 100
        self.scf_scaleparameters = scf_scaleparameters
129 100
        self.scaleparameters = scaleparameters
130 100
        self._SCF = SCF
131 100
        self._thole = thole
132 100
        self._mode = mode
133 100
        self.step = 0
134 100
        self.number_of_lines_in_X = 0
135 100
        self.X_BCC = 0.0
136 100
        self.Y_BCC = 0.0
137 100
        self._FF = FF
138
        # Label the atoms and bonds using a offxml file
139 100
        forcefield = ForceField(os.path.join(ROOT_DIR_PATH, FF))
140

141 100
        self._nalpha = len(forcefield.get_parameter_handler('vdW').parameters)
142 100
        self._nbccs = len(forcefield.get_parameter_handler('Bonds').parameters)
143 100
        self._bccs ={}
144 100
        for i,element in enumerate(forcefield.get_parameter_handler('Bonds').parameters):
145 100
            self._bccs[element.id] = i
146

147

148 100
        self._alpha = {}
149 100
        for i,element in enumerate(forcefield.get_parameter_handler('vdW').parameters):
150 100
            self._alpha[element.id] = i
151

152 100
        self.bccs_old = np.zeros(self._nbccs)
153 100
        self.alphas_old = np.zeros(self._nalpha)
154

155 100
    def load_from_file(self,txtfile):
156
        """
157
        Allows to build a TrainingSet instance from an text file.
158
        File format as followed:
159

160
        :return:
161
        """
162 0
        f = open(txtfile)
163 0
        lines = f.readlines()
164 0
        f.close()
165 0
        for i, line in enumerate(lines):
166 0
            mol2file = line.split()[1]
167
            # noinspection PyTypeChecker
168 0
            self.add_molecule(Molecule(mol2file))
169

170 100
    def add_molecule(self, datei):
171
        """
172
        Adds a molecule.
173
        :param datei: Mol2 file of a molecule
174
        :return:
175
        """
176 100
        self.molecules.append(Molecule(datei, position=self.number_of_lines_in_X, trainingset=self))
177 100
        self.number_of_lines_in_X += self.molecules[-1]._lines_in_X
178

179

180

181

182

183

184

185

186 100
    def build_matrix_A(self):
187
        """
188
        Builds the matrix A of the underlying molecules and combines them.
189

190
        This function is only used for charge optimizatoin RESP
191
        """
192 100
        for molecule in self.molecules:
193 100
            molecule.build_matrix_A()
194 100
            X12 = np.zeros((len(self.A), len(molecule.A)))
195 100
            self.A = np.concatenate(
196
                (np.concatenate((self.A, X12), axis=1), np.concatenate((X12.transpose(), molecule.A), axis=1)), axis=0)
197

198 100
    def build_vector_B(self):
199 100
        for molecule in self.molecules:
200 100
            molecule.build_vector_B()
201 100
            self.B = np.concatenate((self.B, molecule.B))
202

203 100
    def build_matrix_X(self):
204 100
        self.X = np.zeros((0,0))
205
        """
206
        Builds the matrix X of the underlying molecules and combines them.
207

208
        This function is only used for charge optimizatoin RESP
209
        """
210 100
        for molecule in self.molecules:
211 100
            molecule.build_matrix_X()
212 100
            X12 = np.zeros((len(self.X), len(molecule.X)))
213 100
            self.X = np.concatenate(
214
                (np.concatenate((self.X, X12), axis=1), np.concatenate((X12.transpose(), molecule.X), axis=1)), axis=0)
215

216 100
        X12 = np.zeros((len(self.X),len(self.intramolecular_polarization_rst)))
217 100
        X22 =  np.zeros((len(self.intramolecular_polarization_rst),len(self.intramolecular_polarization_rst)))
218

219 100
        self.X = np.concatenate((np.concatenate((self.X, X12), axis=1), np.concatenate((X12.transpose(),X22), axis =1)), axis=0)
220 100
        for i,atoms in enumerate(self.intramolecular_polarization_rst):
221 100
            self.X[self.number_of_lines_in_X + i][atoms[0]] = self.X[atoms[0]][self.number_of_lines_in_X + i] = 1
222 100
            self.X[self.number_of_lines_in_X + i][atoms[1]] = self.X[atoms[1]][self.number_of_lines_in_X + i] = -1
223

224

225 100
    def build_matrix_X_BCC(self):
226
        """
227
        Builds the matrix X of the underlying molecules and combines them.
228

229
        This function is only used for charge optimizatoin RESP
230
        """
231

232 100
        for molecule in self.molecules:
233 100
            molecule.build_matrix_X_BCC()
234 100
            self.X_BCC += molecule.X
235

236

237 100
    def build_vector_Y_BCC(self):
238
        """
239
        Builds the matrix X of the underlying molecules and combines them.
240

241
        This function is only used for charge optimizatoin RESP
242
        """
243

244 100
        for molecule in self.molecules:
245 100
            molecule.build_vector_Y_BCC()
246 100
            self.Y_BCC += molecule.Y
247

248 100
    def build_vector_Y(self):
249 100
        self.Y = np.zeros(0)
250 100
        for molecule in self.molecules:
251 100
            molecule.build_vector_Y()
252 100
            self.Y = np.concatenate((self.Y, molecule.Y))
253

254 100
        Y2 =  np.zeros(len(self.intramolecular_polarization_rst))
255 100
        self.Y = np.concatenate((self.Y, Y2))
256
    # def get_intramolecular_charge_rst()
257

258 100
    @property
259
    def intramolecular_polarization_rst(self):
260

261 100
        intramolecular_polarization_rst = []
262 100
        first_occurrence_of_parameter = {}
263 100
        for molecule in self.molecules:
264 100
            first_occurrence_of_parameter_in_molecule = {}
265 100
            for atom in molecule._atoms:
266 100
                if atom._parameter_id not in first_occurrence_of_parameter.keys():
267 100
                    first_occurrence_of_parameter[atom._parameter_id] = molecule._position_in_A + atom._id + molecule._lines_in_A
268 100
                elif atom._parameter_id not in first_occurrence_of_parameter_in_molecule.keys():
269 100
                    intramolecular_polarization_rst.append(
270
                        [first_occurrence_of_parameter[atom._parameter_id], molecule._position_in_A + atom._id + molecule._lines_in_A])
271 100
        return intramolecular_polarization_rst
272

273 100
    def optimize_charges_alpha(self,criteria = 10E-5):
274 100
        converged = False
275 100
        while converged == False:
276 100
            self.optimize_charges_alpha_step()
277 100
            converged = True
278 100
            for molecule in self.molecules:
279 100
                molecule.step +=1
280 100
                if not all(abs(molecule.q - molecule.q_old) <criteria) or not all(abs(molecule.alpha - molecule.alpha_old) <criteria) :
281 100
                    converged = False
282 100
                molecule.q_old = molecule.q
283 100
                print(molecule.q)
284 100
                print(molecule.alpha)
285 100
                molecule.alpha_old =molecule.alpha
286

287

288

289 100
    def optimize_charges_alpha_step(self):
290 100
        self.build_matrix_X()
291 100
        self.build_vector_Y()
292 100
        self.q_alpha = np.linalg.solve(self.X, self.Y)
293
        # Update the charges of the molecules below
294 100
        q_alpha_tmp = self.q_alpha
295 100
        for molecule in self.molecules:
296 100
            molecule.q_alpha = q_alpha_tmp[:len(molecule.X)]
297 100
            q_alpha_tmp = q_alpha_tmp[len(molecule.X):]
298 100
            molecule.update_q_alpha()
299

300 100
    def optimize_bcc_alpha(self,criteria = 10E-5):
301 100
        converged = False
302 100
        while converged == False:
303 100
            self.optimize_bcc_alpha_step()
304 100
            for molecule in self.molecules:
305 100
                molecule.step +=1
306 100
                print(molecule.q)
307 100
                print(molecule.alpha)
308 100
            if all(abs(self.bccs - self.bccs_old) < criteria) and all(abs(self.alphas - self.alphas_old)< criteria):
309 100
                converged = True
310
            else:
311 100
                self.alphas_old = self.alphas
312 100
                self.bccs_old = self.bccs
313

314

315 100
    def optimize_bcc_alpha_step(self):
316 100
        self.build_matrix_X_BCC()
317 100
        self.build_vector_Y_BCC()
318

319
        # Check if bccs or alphas are not in the training set and set them to zero
320 100
        for i,row in enumerate(self.X_BCC):
321 100
            if all(row == 0.0):
322 100
                self.X_BCC[i][i]=1
323

324 100
        self.bcc_alpha = np.linalg.solve(self.X_BCC, self.Y_BCC)
325
        # Update the charges of the molecules below
326 100
        self.bccs = self.bcc_alpha[:self._nbccs]
327 100
        self.alphas = self.bcc_alpha[self._nbccs:]
328 100
        for molecule in self.molecules:
329 100
            molecule.update_bcc(self.bccs, self.alphas)
330

331

332

333 100
    def optimize_charges(self):
334 100
        self.build_matrix_A()
335 100
        self.build_vector_B()
336 100
        self.q = Q_(np.linalg.solve(self.A, self.B), 'elementary_charge')
337

338
        # Update the charges of the molecules below
339 100
        q_tmp = self.q
340 100
        for molecule in self.molecules:
341 100
            molecule.q = q_tmp[:len(molecule.A)]
342 100
            q_tmp = q_tmp[len(molecule.A):]
343 100
            molecule.update_q()
344

345

346
# =============================================================================================
347
# Molecule
348
# =============================================================================================
349

350

351 100
class Molecule:
352
    """
353
        This class loads in a mol2 file and expose all relevant information. It combines the functionality of
354
        an openeye molecule with the functionality of the OFF molecule. The OE part is only needed to determine
355
        the chemical equivalent atoms. When this feature is implemented in OFF toolkit the OE molecule is not
356
        necessary anymore
357

358
        :param datei: mol2 file
359

360
        :return:
361
        """
362

363 100
    def __init__(self, datei, position=0, trainingset=None):
364

365
        # Get Molecle name from filename
366 100
        self.B = 0.0
367 100
        self.A = 0.0
368 100
        self.X = 0.0
369 100
        self.Y = 0.0
370 100
        self._name = datei.split('/')[-1].strip(".mol2")
371 100
        self._trainingset = trainingset
372

373
        # Number of optimization steps
374 100
        self.step = 0
375

376
        # Postion of this molecule in optimization matrix A
377 100
        self._position_in_A = position
378

379
        # Copy (scf) scaleparameters from the trainingset definition
380 100
        if trainingset is None:
381 100
            self.scf_scaleparameters = None
382 100
            self.scaleparameters = None
383 100
            self._thole = False
384 100
            self._SCF = False
385 100
            self._mode = 'q'
386
        else:
387 100
            self.scf_scaleparameters = self._trainingset.scf_scaleparameters
388 100
            self.scaleparameters = self._trainingset.scaleparameters
389 100
            self._thole = self._trainingset._thole
390 100
            self._SCF = self._trainingset._SCF
391 100
            self._mode = self._trainingset._mode
392
            # Initialize OE Molecule
393 100
        self.oemol = oechem.OEMol()
394

395
        # Open Input File Stream
396 100
        ifs = oechem.oemolistream(datei)
397

398
        # Read from IFS to molecule
399 100
        oechem.OEReadMol2File(ifs, self.oemol)
400

401
        # Check if molecule is a 3 dimensional object
402 100
        if self.oemol.GetDimension() != 3:
403 100
            log.error('Molecule either not found or does not have a 3D structure')
404 100
            raise Exception('Molecule either not found or does not have a 3D structure')
405

406
        # Check for strange atoms
407 100
        AtomicNumbers = []
408 100
        for atom in self.oemol.GetAtoms():
409 100
            AtomicNumbers.append(atom.GetAtomicNum())
410

411 100
            if atom.GetAtomicNum() not in biological_elements:
412 100
                log.warning("""I detect an atom with atomic number: {}
413
                Are you sure you want to include such an atom in your dataset?
414
                Maybe you are using a mol2 file with amber atomtypes""".format(atom.GetAtomicNum()))
415 100
                raise Warning("""I detect an atom with atomic number: {}
416
                Are you sure you want to include such an atom in your dataset?
417
                Maybe you are using a mol2 file with amber atomtypes""".format(atom.GetAtomicNum()))
418

419
        # Generate the OFF molecule from the openeye molecule
420 100
        self.offmol = openff.Molecule.from_openeye(self.oemol)
421 100
        self.offtop = openff.Topology.from_molecules([self.offmol])
422

423
        # Label the atoms and bonds using a offxml file
424 100
        if trainingset is None:
425 100
            self._trainingset = TrainingSet()
426

427 100
        forcefield = ForceField(os.path.join(ROOT_DIR_PATH, self._trainingset._FF))
428

429
        # Run the parameter labeling
430 100
        molecule_parameter_list = forcefield.label_molecules(self.offtop)
431

432
        # set molecule charge
433
        # self._charge=openff.Molecule.total_charge
434 100
        self._charge = 0
435

436
        # Initialize the bonds, atoms
437 100
        self._bonds = list()
438 100
        self._atoms = list()
439

440
        # Define all BCC bonds
441 100
        for i, properties in enumerate(molecule_parameter_list[0]['Bonds'].items()):
442 100
            atom_indices, parameter = properties
443 100
            self.add_bond(i, atom_indices, parameter.id)
444

445 100
        self._nbonds = len(self._bonds)
446

447
        # Define all atomtypes for polarization
448 100
        for i, properties in enumerate(molecule_parameter_list[0]['vdW'].items()):
449 100
            atom_index, parameter = properties
450 100
            self.add_atom(i, atom_index[0], parameter.id, atomic_number = AtomicNumbers[atom_index[0]])
451

452 100
        self._natoms = len(self._atoms)
453

454
        # Initialize and fill scaling matrix
455 100
        self.scale = np.ones((self._natoms, self._natoms))
456 100
        self.scale_scf = np.ones((self._natoms, self._natoms))
457 100
        self.scaling(scf_scaleparameters=self.scf_scaleparameters, scaleparameters=self.scaleparameters)
458 100
        self.create_BCCmatrix_T()
459 100
        self.create_POLmatrix_R()
460

461
        # Initialize conformers
462 100
        self.conformers = list()
463

464
        # Number of lines for matrix X
465 100
        if self._mode == 'q':
466 100
            self._lines_in_X = self._natoms + len(self.chemical_eq_atoms) + 1
467 100
            self.q_old = np.zeros(self._natoms)
468 100
        if self._mode == 'q_alpha':
469 100
            self._lines_in_X = self._natoms + len(
470
                self.chemical_eq_atoms) + 1 + 3 * self._natoms + 2 * self._natoms
471 100
            self.q_old = np.zeros(self._natoms)
472 100
            self.alpha_old = np.zeros(3*self._natoms)
473 100
        self._lines_in_A =  self._natoms + len(self.chemical_eq_atoms) + 1
474
        # Initiliaxe charges
475 100
        self.q_alpha = np.zeros(self._lines_in_X)
476

477 100
    def add_bond(self, index, atom_indices, parameter_id):
478
        """
479
        Adds a bond object ot the molecule
480

481
        :param index:
482
        :param atom_indices:
483
        :param parameter_id:
484
        :return:
485
        """
486 100
        atom_index1 = atom_indices[0]
487 100
        atom_index2 = atom_indices[1]
488 100
        self._bonds.append(Bond(index, atom_index1, atom_index2, parameter_id))
489

490 100
    def add_atom(self, index, atom_index, parameter_id, atomic_number=0):
491
        """
492
        Adds a atom object to the molecule
493

494
        :param index:
495
        :param atom_index:
496
        :param parameter_id:
497
        :return:
498
        """
499 100
        self._atoms.append(Atom(index, atom_index, parameter_id, atomic_number= atomic_number))
500

501 100
    def add_conformer_from_mol2(self, mol2file):
502
        """
503
        Adds a conformer from a mol2 file
504
        Automatically checks if the conformer corresponds to the molecule
505

506
        :param mol2file:
507
        :return:
508
        """
509 100
        conf = openff.Molecule.from_file(mol2file)
510

511
        # Check if molecule is a 3 dimensional object and has the correct dimensions
512
        # Checks if this conformer is from this molecule based on atom names
513 100
        if self.offmol.n_atoms != conf.n_atoms or \
514
                self.offmol.n_bonds != conf.n_bonds or \
515
                self.offmol.n_angles != conf.n_angles:
516 100
            log.error('Conformer and Molecule does not match')
517 100
            raise Exception('Conformer and Molecule does not match')
518
        else:
519 100
            self.conformers.append(Conformer(conf, molecule=self))
520
        # Checks if this conformer is from this moleule based on atom names
521

522 100
    @property
523
    def chemical_eq_atoms(self):
524
        """
525
        # Note: Something similar in compare_charges, which can be copied and slightly modified
526
        :return:
527
        """
528 100
        return find_eq_atoms(self.oemol)
529

530 100
    @property
531
    def same_polarization_atoms(self):
532 100
        array_of_same_pol_atoms = []
533 100
        first_occurrence = {}
534 100
        for atom in self._atoms:
535 100
            if atom._parameter_id not in first_occurrence.keys():
536 100
                first_occurrence[atom._parameter_id] = atom._id
537
            else:
538 100
                array_of_same_pol_atoms.append([first_occurrence[atom._parameter_id], atom._id ])
539 100
        return (array_of_same_pol_atoms)
540

541

542 100
    def create_BCCmatrix_T(self):
543
        """
544
        Creates the bond matrix T for a molecule.
545

546
        Parameters:
547
        ----------
548
            bondtyps: list of [int,int]
549
                Bonds in the set between different BCC groups
550

551
        Attribute:
552
        ----------
553
            bondmatrix T: 2 dim array
554
                1 dim atom
555
                2 dim BCC
556
        See also AM1-BCC paper:
557
        https://onlinelibrary.wiley.com/doi/epdf/10.1002/%28SICI%291096-987X%2820000130%2921%3A2%3C132%3A%3AAID-JCC5%3E3.0.CO%3B2-P
558
        """
559 100
        bondsexcluded = []
560 100
        self.T = np.zeros((self._natoms, self._trainingset._nbccs))
561 100
        for bond in self._bonds:
562 100
            if bond._parameter_id not in bondsexcluded:
563 100
                    self.T[bond._atom1][self._trainingset._bccs[bond._parameter_id]] += 1
564 100
                    self.T[bond._atom2][self._trainingset._bccs[bond._parameter_id]] -= 1
565

566 100
    def create_POLmatrix_R(self):
567
        """
568
        Create a polarization matrix. Defines which atom has which polarizability parameter.
569

570
        Parameters:
571
        ----------
572
        groups: dictionary atomtyp -> int
573
            Stores the polarization group  of each atom
574

575
        Defines which atom belongs to the which pol group.
576
        Takes the list of atomtypes and assign the corresponding group value.
577
        Here I am using a slightly different approach in contrast to the atomtyps.
578
        as the number of different pol types is maximal the number of atom types. No pre-screening for the involved
579
        atom-types is needed
580
        """
581 100
        self.R = np.zeros((self._natoms, self._trainingset._nalpha))
582 100
        for atom in self._atoms:
583 100
            self.R[atom._id][self._trainingset._alpha[atom._parameter_id]] += 1
584

585 100
    def update_bcc(self, bccs, alphas):
586
        """
587
        Converts the optimized bccs and alphas to a qd object.
588

589
        :param mol2: molecule class object
590
        :return:
591
        """
592 100
        for i in range(self._natoms):
593 100
            self.q_alpha[i] = np.dot(self.T[i], bccs)
594 100
        self.alpha = [np.dot(self.R[i], alphas) for i in range(len(self.R))]
595 100
        self.q_alpha[self._lines_in_A:self._lines_in_A + 3* len(self.alpha)] = np.concatenate(
596
            (np.concatenate((self.alpha, self.alpha)), self.alpha))
597 100
        self.q = self.q_alpha[:self._natoms]
598 100
        for conformer in self.conformers:
599 100
            conformer.q_alpha = self.q_alpha
600

601

602 100
    def build_matrix_A(self):
603 100
        for conformer in self.conformers:
604 100
            conformer.build_matrix_A()
605 100
            self.A += conformer.A
606

607 100
    def build_vector_B(self):
608 100
        for conformer in self.conformers:
609 100
            conformer.build_vector_B()
610 100
            self.B += conformer.B
611

612 100
    def build_matrix_X(self):
613 100
        for conformer in self.conformers:
614 100
            conformer.build_matrix_X()
615 100
            self.X += conformer.X
616

617 100
    def build_matrix_X_BCC(self):
618 100
        for conformer in self.conformers:
619 100
            conformer.build_matrix_X_BCC()
620 100
            self.X += conformer.X
621

622 100
    def build_vector_Y(self):
623 100
        for conformer in self.conformers:
624 100
            conformer.build_vector_Y()
625 100
            self.Y += conformer.Y
626

627 100
    def build_vector_Y_BCC(self):
628 100
        for conformer in self.conformers:
629 100
            conformer.build_vector_Y_BCC()
630 100
            self.Y += conformer.Y
631

632 100
    def optimize_charges(self):
633 100
        self.build_matrix_A()
634 100
        self.build_vector_B()
635 100
        self.q = Q_(np.linalg.solve(self.A, self.B), 'elementary_charge')
636

637 100
    def optimize_charges_alpha(self):
638 0
        self.build_matrix_X()
639 0
        self.build_vector_Y()
640 0
        self.q_alpha = Q_(np.linalg.solve(self.X, self.Y), 'elementary_charge')
641

642 100
    def update_q(self):
643 100
        for conformer in self.conformers:
644 100
            conformer.q = self.q
645

646 100
    def update_q_alpha(self):
647 100
        self.q = self.q_alpha[:self._natoms]
648 100
        self.alpha = self.q_alpha[self._natoms+1+len(self.chemical_eq_atoms):4*self._natoms+1+len(self.chemical_eq_atoms):]
649 100
        for conformer in self.conformers:
650 100
            conformer.q_alpha = self.q_alpha
651

652 100
    def scaling(self, scaleparameters=None, scf_scaleparameters=None):
653
        """
654
        Takes the bond information from a molecule instances and converts it to an scaling matrix.
655

656
        Parameters:
657
        ---------
658
        bonds:  list of [int,int]
659
            list of atoms connected with each other
660
        scaleparameters: [float,float,float]
661
            1-2 scaling, 1-3 scaling, 1-4 scaling parameter
662

663
        Attributes:
664
        ---------
665
        scale: matrix
666
            scaling matrix how atoms interact with each other
667

668
        """
669

670
        # Initializing
671 100
        bound12 = np.zeros((self._natoms, self._natoms))
672 100
        bound13 = np.zeros((self._natoms, self._natoms))
673 100
        bound14 = np.zeros((self._natoms, self._natoms))
674 100
        if scaleparameters is None:
675 100
            scaleparameters = [0.0, 0.0, 0.8333333333]
676

677
        # Building connection matrix
678 100
        for bond in self._bonds:
679 100
            bound12[bond._atom1][bond._atom2] = 1.0
680 100
            bound12[bond._atom2][bond._atom1] = 1.0
681

682 100
        for i in range(len(bound12)):
683 100
            b12 = np.where(bound12[i] == 1.0)[0]
684 100
            for j in range(len(b12)):
685 100
                b12t = np.where(bound12[b12[j]] == 1.0)[0]
686 100
                for k in range(len(b12t)):
687 100
                    if i != b12t[k]:
688 100
                        bound13[b12t[k]][i] = 1.0
689 100
                        bound13[i][b12t[k]] = 1.0
690

691 100
        for i in range(self._natoms):
692 100
            b13 = np.where(bound13[i] == 1.0)[0]
693 100
            for j in range(len(b13)):
694 100
                b13t = np.where(bound12[b13[j]] == 1.0)[0]
695 100
                for k in range(len(b13t)):
696 100
                    if bound12[b13t[k]][i] == 0.0:
697 100
                        bound14[b13t[k]][i] = 1.0
698 100
                        bound14[i][b13t[k]] = 1.0
699

700 100
        for i in range(self._natoms):
701 100
            self.scale[i][i] = 0.0
702
        # find values in matrix with value 1.0
703 100
        b12 = np.array(np.where(bound12 == 1.0)).transpose()
704 100
        b13 = np.array(np.where(bound13 == 1.0)).transpose()
705 100
        b14 = np.array(np.where(bound14 == 1.0)).transpose()
706

707
        # Fill scaling matrix with values
708 100
        for i in range(len(b12)):
709 100
            self.scale[b12[i][0]][b12[i][1]] = scaleparameters[
710
                0]  # Value for 1-2 interaction 0 means interactions are neglected
711 100
        for i in range(len(b13)):
712 100
            self.scale[b13[i][0]][b13[i][1]] = scaleparameters[
713
                1]  # Value for 1-3 interaction 0 means interactions are neglected
714 100
        for i in range(len(b14)):
715 100
            self.scale[b14[i][0]][b14[i][1]] = scaleparameters[2]  # Value for the 1-4 scaling
716

717
        # Different Scaling parameter for SCF
718 100
        if scf_scaleparameters != None:
719 100
            self.scale_scf = np.ones((self._natoms, self._natoms))
720 100
            for i in range(self._natoms):
721 100
                self.scale_scf[i][i] = 0.0
722 100
            for i in range(len(b12)):
723 100
                self.scale_scf[b12[i][0]][b12[i][1]] = scf_scaleparameters[
724
                    0]  # Value for 1-2 interaction 0 means interactions are neglected
725 100
            for i in range(len(b13)):
726 100
                self.scale_scf[b13[i][0]][b13[i][1]] = scf_scaleparameters[
727
                    1]  # Value for 1-3 interaction 0 means interactions are neglected
728 100
            for i in range(len(b14)):
729 100
                self.scale_scf[b14[i][0]][b14[i][1]] = scf_scaleparameters[2]  # Value for the 1-4 scaling
730

731

732
# =============================================================================================
733
# Atom
734
# =============================================================================================
735

736 100
class Atom:
737
    """
738

739
    """
740

741 100
    def __init__(self, index, atom_index, parameter_id, atomic_number = 0):
742 100
        self._id = index
743 100
        self._atom = atom_index
744 100
        self._parameter_id = parameter_id
745 100
        self._atomic_number = atomic_number
746

747

748
# =============================================================================================
749
# Bond
750
# =============================================================================================
751

752 100
class Bond:
753

754 100
    def __init__(self, index, atom_index1, atom_index2, parameter_id):
755 100
        self._id = index
756 100
        self._atom1 = atom_index1
757 100
        self._atom2 = atom_index2
758 100
        self._parameter_id = parameter_id
759

760

761
# =============================================================================================
762
# Conformer
763
# =============================================================================================
764

765 100
class Conformer:
766
    """
767

768
    """
769

770 100
    def __init__(self, conf, molecule=None):
771
        """
772

773
        :param conf: openff molecule file
774
        """
775
        # noinspection PyProtectedMember
776 100
        self.atom_positions = Q_(np.array(conf.conformers[0]._value), 'angstrom')
777 100
        self.atom_positions_angstrom = self.atom_positions.to('angstrom').magnitude
778 100
        self.natoms = len(self.atom_positions.magnitude)
779 100
        self.baseESP = None
780 100
        self.polESPs = list()
781 100
        self._molecule = molecule
782 100
        self.q_alpha = self._molecule.q_alpha
783 100
        self._lines_in_A = self._molecule._lines_in_A
784

785
        # Initiliaze Electric field vectors
786 100
        self.e_field_at_atom = np.zeros((3, self.natoms))
787

788 100
    def get_grid_coord(self):
789 100
        self.grid_coord_angstrom = self.baseESP.positions.to('angstrom').magnitude
790 100
        self.npoints = len(self.baseESP.positions.magnitude)
791

792 100
    def add_baseESP(self, *args, ):
793
        """
794
        Adds the unpolarized molecule to this conformation.
795

796
        :param args:  ESPF file
797
        1.) GESP file form g09
798
        2.) grid.dat file and esp.dat file generated via respyte and psi4
799

800
        :return:
801
        """
802 100
        self.baseESP = BCCUnpolESP(*args, conformer=self)
803

804
        # Check if atomic coordinates match
805 100
        for atom1 in self.baseESP.atom_positions:
806 100
            atom_is_present = 0
807 100
            for atom2 in self.atom_positions:
808 100
                if np.linalg.norm(atom1.to('angstrom').magnitude - atom2.to('angstrom').magnitude) < 0.01:
809 100
                    atom_is_present = 1
810 100
            if atom_is_present == 0:
811 0
                raise Exception("ESP grid does not belong to the conformer")
812

813 100
    def add_polESP(self, *args, e_field=Q_([0.0, 0.0, 0.0], 'bohr')):
814
        """
815
        Adds the unpolarized molecule to this conformation.
816

817
        :param args:  ESPF file
818
        1.) GESP file form g09
819
        2.) grid.dat file and esp.dat file generated via respyte and psi4
820
        :param:e_field: Pint formatted electrif field
821
        e.g e_field=Q_([0.0, 0.0, 0.0], 'bohr'))
822

823
        :return:
824
        """
825 100
        self.polESPs.append(BCCPolESP(*args, conformer=self, e_field=e_field))
826

827
    # Build the matrix A for the charge optimization
828 100
    def build_matrix_A(self):
829

830
        """
831
        Fast method for only optimizing charges.
832

833
        :return:
834
        """
835
        # Determine size of matrix for this molecule
836
        # every atom is one line
837
        # one li