1
# !/usr/bin/env python
2

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

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

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

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

45

46
# =============================================================================================
47
# PRIVATE SUBROUTINES
48
# =============================================================================================
49 4
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 4
    qmol = oechem.OEQMol()
62

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

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

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

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

98 4
    return sorted_eq_atoms
99

100

101

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

106 4
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 4
    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 4
        self.molecules = list()
124 4
        self.B = np.zeros(0)
125 4
        self.A = np.zeros((0, 0))
126 4
        self.q = 0.0
127 4
        self.mode = mode
128 4
        self.scf_scaleparameters = scf_scaleparameters
129 4
        self.scaleparameters = scaleparameters
130 4
        self._SCF = SCF
131 4
        self._thole = thole
132 4
        self._mode = mode
133 4
        self.step = 0
134 4
        self.number_of_lines_in_X = 0
135 4
        self.X_BCC = 0.0
136 4
        self.Y_BCC = 0.0
137 4
        self._FF = FF
138
        # Label the atoms and bonds using a offxml file
139 4
        forcefield = ForceField(os.path.join(ROOT_DIR_PATH, FF))
140

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

147

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

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

155 4
    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 4
    def add_molecule(self, datei):
171
        """
172
        Adds a molecule.
173
        :param datei: Mol2 file of a molecule
174
        :return:
175
        """
176 4
        self.molecules.append(Molecule(datei, position=self.number_of_lines_in_X, trainingset=self))
177 4
        self.number_of_lines_in_X += self.molecules[-1]._lines_in_X
178

179

180

181

182

183

184

185

186 4
    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 4
        for molecule in self.molecules:
193 4
            molecule.build_matrix_A()
194 4
            X12 = np.zeros((len(self.A), len(molecule.A)))
195 4
            self.A = np.concatenate(
196
                (np.concatenate((self.A, X12), axis=1), np.concatenate((X12.transpose(), molecule.A), axis=1)), axis=0)
197

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

203 4
    def build_matrix_X(self):
204 4
        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 4
        for molecule in self.molecules:
211 4
            molecule.build_matrix_X()
212 4
            X12 = np.zeros((len(self.X), len(molecule.X)))
213 4
            self.X = np.concatenate(
214
                (np.concatenate((self.X, X12), axis=1), np.concatenate((X12.transpose(), molecule.X), axis=1)), axis=0)
215

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

219 4
        self.X = np.concatenate((np.concatenate((self.X, X12), axis=1), np.concatenate((X12.transpose(),X22), axis =1)), axis=0)
220 4
        for i,atoms in enumerate(self.intramolecular_polarization_rst):
221 4
            self.X[self.number_of_lines_in_X + i][atoms[0]] = self.X[atoms[0]][self.number_of_lines_in_X + i] = 1
222 4
            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 4
    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 4
        for molecule in self.molecules:
233 4
            molecule.build_matrix_X_BCC()
234 4
            self.X_BCC += molecule.X
235

236

237 4
    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 4
        for molecule in self.molecules:
245 4
            molecule.build_vector_Y_BCC()
246 4
            self.Y_BCC += molecule.Y
247

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

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

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

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

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

287

288

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

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

314

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

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

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

331

332

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

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

345

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

350

351 4
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 4
    def __init__(self, datei, position=0, trainingset=None):
364

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

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

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

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

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

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

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

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

411 4
            if atom.GetAtomicNum() not in biological_elements:
412 4
                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 4
                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 4
        self.offmol = openff.Molecule.from_openeye(self.oemol)
421 4
        self.offtop = openff.Topology.from_molecules([self.offmol])
422

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

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

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

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

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

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

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

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

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

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

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

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

477 4
    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 4
        atom_index1 = atom_indices[0]
487 4
        atom_index2 = atom_indices[1]
488 4
        self._bonds.append(Bond(index, atom_index1, atom_index2, parameter_id))
489

490 4
    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 4
        self._atoms.append(Atom(index, atom_index, parameter_id, atomic_number= atomic_number))
500

501 4
    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 4
        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 4
        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 4
            log.error('Conformer and Molecule does not match')
517 4
            raise Exception('Conformer and Molecule does not match')
518
        else:
519 4
            self.conformers.append(Conformer(conf, molecule=self))
520
        # Checks if this conformer is from this moleule based on atom names
521

522 4
    @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 4
        return find_eq_atoms(self.oemol)
529

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

541

542 4
    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 4
        bondsexcluded = []
560 4
        self.T = np.zeros((self._natoms, self._trainingset._nbccs))
561 4
        for bond in self._bonds:
562 4
            if bond._parameter_id not in bondsexcluded:
563 4
                    self.T[bond._atom1][self._trainingset._bccs[bond._parameter_id]] += 1
564 4
                    self.T[bond._atom2][self._trainingset._bccs[bond._parameter_id]] -= 1
565

566 4
    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 4
        self.R = np.zeros((self._natoms, self._trainingset._nalpha))
582 4
        for atom in self._atoms:
583 4
            self.R[atom._id][self._trainingset._alpha[atom._parameter_id]] += 1
584

585 4
    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 4
        for i in range(self._natoms):
593 4
            self.q_alpha[i] = np.dot(self.T[i], bccs)
594 4
        self.alpha = [np.dot(self.R[i], alphas) for i in range(len(self.R))]
595 4
        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 4
        self.q = self.q_alpha[:self._natoms]
598 4
        for conformer in self.conformers:
599 4
            conformer.q_alpha = self.q_alpha
600

601

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

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

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

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

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

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

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

637 4
    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 4
    def update_q(self):
643 4
        for conformer in self.conformers:
644 4
            conformer.q = self.q
645

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

652 4
    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 4
        bound12 = np.zeros((self._natoms, self._natoms))
672 4
        bound13 = np.zeros((self._natoms, self._natoms))
673 4
        bound14 = np.zeros((self._natoms, self._natoms))
674 4
        if scaleparameters is None:
675 4
            scaleparameters = [0.0, 0.0, 0.8333333333]
676

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

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

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

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

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

717
        # Different Scaling parameter for SCF
718 4
        if scf_scaleparameters != None:
719 4
            self.scale_scf = np.ones((self._natoms, self._natoms))
720 4
            for i in range(self._natoms):
721 4
                self.scale_scf[i][i] = 0.0
722 4
            for i in range(len(b12)):
723 4
                self.scale_scf[b12[i][0]][b12[i][1]] = scf_scaleparameters[
724
                    0]  # Value for 1-2 interaction 0 means interactions are neglected
725 4
            for i in range(len(b13)):
726 4
                self.scale_scf[b13[i][0]][b13[i][1]] = scf_scaleparameters[
727
                    1]  # Value for 1-3 interaction 0 means interactions are neglected
728 4
            for i in range(len(b14)):
729 4
                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 4
class Atom:
737
    """
738

739
    """
740

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

747

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

752 4
class Bond:
753

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

760

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

765 4
class Conformer:
766
    """
767

768
    """
769

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

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

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

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

792 4
    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 4
        self.baseESP = BCCUnpolESP(*args, conformer=self)
803

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

813 4
    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 4
        self.polESPs.append(BCCPolESP(*args, conformer=self, e_field=e_field))
826

827
    # Build the matrix A for the charge optimization
828 4
    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 line to restrain the overall charge of the molecule
838
        # one line for every pair of equivalent atoms
839 4
        self.get_distances()
840 4
        self.A = np.zeros((self._lines_in_A, self._lines_in_A))
841 4
        for j in range(self.natoms):
842 4
            for k in range(j + 1):
843 4
                self.A[j][k] = np.dot(self.dist[j], self.dist[k])
844

845
        # Symmetric matrix -> copy diagonal elements, add total charge restrain
846 4
        for j in range(self.natoms):
847 4
            for k in range(j):
848 4
                self.A[k][j] = self.A[j][k]
849 4
            self.A[self.natoms][j] = 1.0
850 4
            self.A[j][self.natoms] = 1.0
851

852
        # Add charge restraints for equivalent atoms
853 4
        for k, eq_atoms in enumerate(self._molecule.chemical_eq_atoms):
854 4
            if eq_atoms[1] > 0:
855 4
                self.A[self.natoms + 1 + k][eq_atoms[0]] = 1
856 4
                self.A[self.natoms + 1 + k][eq_atoms[1]] = -1
857 4
                self.A[eq_atoms[0]][self.natoms + 1 + k] = 1
858 4
                self.A[eq_atoms[1]][self.natoms + 1 + k] = -1
859

860 4
    def build_matrix_D(self):
861
        """
862
        Method for building the polarization matrix D. Only used for fitting polarizations withotut charges or BCCs.
863

864
        :return:
865
        """
866
        # 1 dipole vector of length 3 per atom
867
        # restrains for isotropic polarization
868
        # Restaint atoms with same polarization parameters
869 4
        self.Dlines = (3 * self.natoms + 2 * self.natoms) #+ len(self._molecule.same_polarization_atoms))
870 4
        self.D = np.zeros((self.Dlines, self.Dlines))
871

872 4
        for j in range(self.natoms):
873 4
            for k in range(self.natoms):
874 4
                self.D[k][j] = np.multiply(np.dot(self.dist_x[k], self.dist_x[j]),
875
                                           self.e_field_at_atom[0][k] * self.e_field_at_atom[0][j])
876 4
                self.D[j + self.natoms][k] = self.D[k][j + self.natoms] = np.multiply(
877
                    np.dot(self.dist_x[k], self.dist_y[j]), self.e_field_at_atom[0][k] * self.e_field_at_atom[1][j])
878 4
                self.D[j + 2 * self.natoms][k] = self.D[k][j + 2 * self.natoms] = np.multiply(
879
                    np.dot(self.dist_x[k], self.dist_z[j]), self.e_field_at_atom[0][k] * self.e_field_at_atom[2][j])
880 4
                self.D[k + self.natoms][j + self.natoms] = np.multiply(np.dot(self.dist_y[k], self.dist_y[j]),
881
                                                                       self.e_field_at_atom[1][k] *
882
                                                                       self.e_field_at_atom[1][j])
883 4
                self.D[j + 2 * self.natoms][k + self.natoms] = self.D[k + self.natoms][
884
                    j + 2 * self.natoms] = np.multiply(np.dot(self.dist_y[k], self.dist_z[j]),
885
                                                       self.e_field_at_atom[1][k] * self.e_field_at_atom[2][j])
886 4
                self.D[k + 2 * self.natoms][j + 2 * self.natoms] = np.multiply(
887
                    np.dot(self.dist_z[k], self.dist_z[j]), self.e_field_at_atom[2][k] * self.e_field_at_atom[2][j])
888
        # Add dipole restraints for equivalent atoms /only works for isotropic suff now
889 4
        for polESP in self.polESPs:
890 4
            for j in range(self.natoms):
891 4
                for k in range(self.natoms):
892 4
                    self.D[k][j] += np.multiply(np.dot(self.dist_x[k], self.dist_x[j]),
893
                                                polESP.e_field_at_atom[0][k] * polESP.e_field_at_atom[0][j])
894 4
                    self.D[j + self.natoms][k] = self.D[k][j + self.natoms] = self.D[k][j + self.natoms] + np.multiply(
895
                        np.dot(self.dist_x[k], self.dist_y[j]),
896
                        polESP.e_field_at_atom[0][k] * polESP.e_field_at_atom[1][j])
897 4
                    self.D[j + 2 * self.natoms][k] = self.D[k][j + 2 * self.natoms] = self.D[k][
898
                                                                                          j + 2 * self.natoms] + np.multiply(
899
                        np.dot(self.dist_x[k], self.dist_z[j]),
900
                        polESP.e_field_at_atom[0][k] * polESP.e_field_at_atom[2][j])
901 4
                    self.D[k + self.natoms][j + self.natoms] += np.multiply(np.dot(self.dist_y[k], self.dist_y[j]),
902
                                                                            polESP.e_field_at_atom[1][k] *
903
                                                                            polESP.e_field_at_atom[1][j])
904 4
                    self.D[j + 2 * self.natoms][k + self.natoms] = self.D[k + self.natoms][j + 2 * self.natoms] = \
905
                    self.D[k + self.natoms][j + 2 * self.natoms] + np.multiply(np.dot(self.dist_y[k], self.dist_z[j]),
906
                                                                               polESP.e_field_at_atom[1][k] *
907
                                                                               polESP.e_field_at_atom[2][j])
908 4
                    self.D[k + 2 * self.natoms][j + 2 * self.natoms] += np.multiply(
909
                        np.dot(self.dist_z[k], self.dist_z[j]),
910
                        polESP.e_field_at_atom[2][k] * polESP.e_field_at_atom[2][j])
911

912 4
        self.D = self.D / (len(self.polESPs) + 1)
913

914
        #for j, atoms in enumerate(self._molecule.same_polarization_atoms):
915
        #    self.D[5 * self.natoms + j][atoms[0]] = 1
916
        #    self.D[5 * self.natoms + j][atoms[1]] = -1
917
        #    self.D[atoms[0]][5 * self.natoms + j] = 1
918
        #    self.D[atoms[1]][5 * self.natoms + j] = -1
919
        # Code to keep polarization parameters at their initial value. Implmentation will change
920
        # elif self.eqdipoles[j][1] < 0:
921
        #    self.D[self.ndipoles + self.aniso_lines + j][self.eqdipoles[j][0]] = 1
922
        #    self.D[self.eqdipoles[j][0]][self.ndipoles + self.aniso_lines + j] = 1
923

924
        # Add restraints for polarization isotropy
925 4
        for j in range(self.natoms):
926 4
            self.D[3 * self.natoms + j][j] = self.D[j][3 * self.natoms + j] = 1.0
927 4
            self.D[3 * self.natoms + j][j + self.natoms] = self.D[j + self.natoms][3 * self.natoms + j] = -1.0
928 4
            self.D[4 * self.natoms + j][j] = self.D[j][4 * self.natoms + j] = 1.0
929 4
            self.D[4 * self.natoms + j][j + 2 * self.natoms] = self.D[j + 2 * self.natoms][4 * self.natoms + j] = -1.0
930

931 4
    def build_matrix_X(self):
932
        """
933
        Creates Matrix X for the RESP-POl method.
934

935
        RESP and Polarization with the large matrix.
936
        Probably worth changing it to the new model.
937

938
        Again the math is shown in the manuscript.
939
        """
940 4
        self.build_matrix_A()
941 4
        self.get_electric_field()
942 4
        self.build_matrix_D()
943 4
        self.B = np.zeros((self._lines_in_A, self.Dlines))
944 4
        self.C = np.zeros((self.Dlines, self._lines_in_A))
945

946
        # Matrix element B see notes
947 4
        for k in range(self.natoms):
948 4
            for j in range(self.natoms):
949 4
                self.B[k][j] = np.multiply(np.dot(self.dist[k], self.dist_x[j]), self.e_field_at_atom[0][j])  # B1
950 4
                self.B[k][self.natoms + j] = np.multiply(np.dot(self.dist[k], self.dist_y[j]),
951
                                                         self.e_field_at_atom[1][j])  # B2
952 4
                self.B[k][2 * self.natoms + j] = np.multiply(np.dot(self.dist[k], self.dist_z[j]),
953
                                                             self.e_field_at_atom[2][j])  # B3
954

955 4
        for polESP in self.polESPs:
956 4
            for k in range(self.natoms):
957 4
                for j in range(self.natoms):
958 4
                    self.B[k][j] += np.multiply(np.dot(self.dist[k], self.dist_x[j]),
959
                                                polESP.e_field_at_atom[0][j])  # B1
960 4
                    self.B[k][self.natoms + j] += np.multiply(np.dot(self.dist[k], self.dist_y[j]),
961
                                                              polESP.e_field_at_atom[1][j])  # B2
962 4
                    self.B[k][2 * self.natoms + j] += np.multiply(np.dot(self.dist[k], self.dist_z[j]),
963
                                                                  polESP.e_field_at_atom[2][j])  # B3
964
        # matrix element C see notes
965
        # matrix element C
966 4
        for j in range(self.natoms):
967 4
            for k in range(self.natoms):
968 4
                self.C[k][j] = np.multiply(np.dot(self.dist[j], self.dist_x[k]), self.e_field_at_atom[0][k])
969 4
                self.C[self.natoms + k][j] = np.multiply(np.dot(self.dist[j], self.dist_y[k]),
970
                                                         self.e_field_at_atom[1][k])
971 4
                self.C[2 * self.natoms + k][j] = np.multiply(np.dot(self.dist[j], self.dist_z[k]),
972
                                                             self.e_field_at_atom[2][k])
973

974 4
        for polESP in self.polESPs:
975 4
            for j in range(self.natoms):
976 4
                for k in range(self.natoms):
977 4
                    self.C[k][j] += np.multiply(np.dot(self.dist[j], self.dist_x[k]), polESP.e_field_at_atom[0][k])
978 4
                    self.C[self.natoms + k][j] += np.multiply(np.dot(self.dist[j], self.dist_y[k]),
979
                                                              polESP.e_field_at_atom[1][k])
980 4
                    self.C[2 * self.natoms + k][j] += np.multiply(np.dot(self.dist[j], self.dist_z[k]),
981
                                                                  polESP.e_field_at_atom[2][k])
982
        # Normalize B and C
983 4
        self.B = self.B / (len(self.polESPs) + 1)
984 4
        self.C = self.C / (len(self.polESPs) + 1)
985

986
        # Combine all matrices
987 4
        self.X = np.concatenate(
988
            (np.concatenate((self.A, self.B), axis=1), np.concatenate((self.C, self.D), axis=1)), axis=0)
989

990

991
    # Combines BCC and polarizabilities
992 4
    def build_matrix_X_BCC(self):
993
        """
994
        Creates the Matrix A for the BCC-POL approach:
995

996
        :param mol2: molecule class object
997

998
        :return:
999

1000
        The math and the optimization method is explained in the current paper draft.
1001
        """
1002 4
        self.get_distances()
1003 4
        nbcc = self._molecule._trainingset._nbccs
1004 4
        nalpha = self._molecule._trainingset._nalpha
1005 4
        T = self._molecule.T
1006 4
        R = self._molecule.R
1007 4
        self.get_electric_field()
1008
        # BCC part
1009 4
        self.A = np.zeros((nbcc, nbcc))
1010 4
        if self._molecule._mode == 'alpha':  # Do not optimize BCCs in that case
1011 0
            for alpha in range(nbcc):
1012 0
                self.A[alpha][alpha] = 1
1013
        else:
1014
            #for ESP in [self.baseESP] + self.polESPs: #lgtm [py/unused-loop-variable]
1015 4
            for j in range(self.natoms):
1016 4
                    for k in range(self.natoms):
1017 4
                        for alpha in range(nbcc):
1018 4
                            for beta in range(nbcc):
1019 4
                                self.A[alpha][beta] += T[j][alpha] * T[k][beta] * np.dot(self.dist[j],
1020
                                                                                               self.dist[k])
1021 4
            self.A = self.A *(len(self.polESPs)+1)
1022
        # Polarizabilities part
1023 4
        self.D = np.zeros((nalpha, nalpha))
1024 4
        if self._molecule._mode != 'bcconly':  # Mode is not optimizing polarizabilies
1025 4
            for ESP in [self.baseESP] + self.polESPs:
1026 4
                for j in range(self.natoms):
1027 4
                    for k in range(self.natoms):
1028 4
                        for alpha in range(nalpha):
1029 4
                            for beta in range(nalpha):
1030 4
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1031
                                    np.dot(self.dist_x[j], self.dist_x[k]), ESP.e_field_at_atom[0][j] * ESP.e_field_at_atom[0][k])
1032 4
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1033
                                    np.dot(self.dist_x[j], self.dist_y[k]), ESP.e_field_at_atom[0][j] * ESP.e_field_at_atom[1][k])
1034 4
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1035
                                    np.dot(self.dist_x[j], self.dist_z[k]), ESP.e_field_at_atom[0][j] * ESP.e_field_at_atom[2][k])
1036 4
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1037
                                    np.dot(self.dist_y[j], self.dist_x[k]), ESP.e_field_at_atom[1][j] * ESP.e_field_at_atom[0][k])
1038 4
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1039
                                    np.dot(self.dist_y[j], self.dist_y[k]), ESP.e_field_at_atom[1][j] * ESP.e_field_at_atom[1][k])
1040 4
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1041
                                    np.dot(self.dist_y[j], self.dist_z[k]), ESP.e_field_at_atom[1][j] * ESP.e_field_at_atom[2][k])
1042 4
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1043
                                    np.dot(self.dist_z[j], self.dist_x[k]), ESP.e_field_at_atom[2][j] * ESP.e_field_at_atom[0][k])
1044 4
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1045
                                    np.dot(self.dist_z[j], self.dist_y[k]), ESP.e_field_at_atom[2][j] * ESP.e_field_at_atom[1][k])
1046 4
                                self.D[alpha][beta] += R[j][alpha] * R[k][beta] * np.multiply(
1047
                                    np.dot(self.dist_z[j], self.dist_z[k]), ESP.e_field_at_atom[2][j] * ESP.e_field_at_atom[2][k])
1048
        else:
1049 0
            for alpha in range(nalpha):
1050 0
                self.D[alpha][alpha] = 1
1051

1052
        # Cross interaction between BCC charges and polarizations
1053 4
        self.B = np.zeros((nbcc, nalpha))
1054 4
        self.C = np.zeros((nalpha, nbcc))
1055 4
        if self._molecule._mode != 'alpha':
1056 4
            for ESP in [self.baseESP] + self.polESPs:
1057 4
                for j in range(self.natoms):
1058 4
                    for k in range(self.natoms):
1059 4
                        for alpha in range(nbcc):
1060 4
                            for beta in range(nalpha):
1061 4
                                self.B[alpha][beta] += T[j][alpha] * R[k][beta] * np.multiply(
1062
                                    np.dot(self.dist[j], self.dist_x[k]), ESP.e_field_at_atom[0][k])
1063 4
                                self.B[alpha][beta] += T[j][alpha] * R[k][beta] * np.multiply(
1064
                                    np.dot(self.dist[j], self.dist_y[k]), ESP.e_field_at_atom[1][k])
1065 4
                                self.B[alpha][beta] += T[j][alpha] * R[k][beta] * np.multiply(
1066
                                    np.dot(self.dist[j], self.dist_z[k]), ESP.e_field_at_atom[2][k])
1067 4
        if self._molecule._mode != 'bcconly':
1068 4
            for ESP in [self.baseESP] + self.polESPs:
1069 4
                for j in range(self.natoms):
1070 4
                    for k in range(self.natoms):
1071 4
                        for alpha in range(nalpha):
1072 4
                            for beta in range(nbcc):
1073 4
                                self.C[alpha][beta] += R[j][alpha] * T[k][beta] * np.multiply(
1074
                                    np.dot(self.dist[k], self.dist_x[j]), ESP.e_field_at_atom[0][j])
1075 4
                                self.C[alpha][beta] += R[j][alpha] * T[k][beta] * np.multiply(
1076
                                    np.dot(self.dist[k], self.dist_y[j]), ESP.e_field_at_atom[1][j])
1077 4
                                self.C[alpha][beta] += R[j][alpha] * T[k][beta] * np.multiply(
1078
                                    np.dot(self.dist[k], self.dist_z[j]), ESP.e_field_at_atom[2][j])
1079

1080
        # Restraints for polarizaton
1081
        #if hasattr(self, 'wrst1'):
1082
        #    for alpha in range(nalpha):
1083
        #        self.D[alpha][alpha] += mol2.group_pop[alpha] * self.wrst1 * self.dipscale / np.sqrt(
1084
        #           np.square(self.qd[nbcc + alpha]) + 0.01)
1085

1086
        # No implementation of bcc restraints.
1087

1088

1089
        #Combine all matrices
1090 4
        self.X = np.concatenate((np.concatenate((self.A, self.B), axis=1), np.concatenate((self.C, self.D), axis=1)),
1091
                                axis=0)
1092

1093

1094 4
        """
1095
        X12= np.zeros((len(self.X), len(self._molecule.same_polarization_atoms)))
1096
        X22 = np.zeros((len(self._molecule.same_polarization_atoms),len(self._molecule.same_polarization_atoms)))
1097
        self.X = np.concatenate((np.concatenate((self.X, X12), axis=1), np.concatenate((X12.transpose, X22), axis=1)),axis=0)
1098
        
1099
        
1100
        for j, atoms in enumerate(self._molecule.same_polarization_atoms):
1101
            self.X[5 * self.natoms + j] = 0.0
1102

1103
        """
1104 4
    def build_vector_Y_BCC(self):
1105
        """
1106
        Creates vector Y for the BCC Pol method
1107

1108
        :param mol2: molecule object
1109

1110
        :return:
1111
        """
1112

1113 4
        nbcc = self._molecule._trainingset._nbccs
1114 4
        nalpha = self._molecule._trainingset._nalpha
1115 4
        T = self._molecule.T
1116 4
        R = self._molecule.R
1117

1118
        # Vector belonging to the BCCs
1119 4
        self.Y1 = np.zeros(nbcc)
1120 4
        if self._molecule._mode != 'alpha':
1121 4
            for ESP in [self.baseESP] + self.polESPs:
1122 4
                esp_values = ESP.esp_values.to('elementary_charge / angstrom').magnitude
1123 4
                for beta in range(nbcc):
1124 4
                    for k in range(self.natoms):
1125 4
                        self.Y1[beta] += T[k][beta] * np.dot(esp_values, self.dist[k])
1126
        else:
1127 0
            self.Y1 = self.qbond
1128

1129
        # Vector belonging to the polarizabilities
1130 4
        self.Y2 = np.zeros(nalpha)
1131 4
        if self._molecule._mode != 'bcconly':
1132 4
            for ESP in [self.baseESP] + self.polESPs:
1133 4
                esp_values = ESP.esp_values.to('elementary_charge / angstrom').magnitude
1134 4
                for beta in range(nalpha):
1135 4
                    for k in range(self.natoms):
1136 4
                        self.Y2[beta] += R[k][beta] * np.multiply(np.dot(esp_values, self.dist_x[k]), ESP.e_field_at_atom[0][k])
1137 4
                        self.Y2[beta] += R[k][beta] * np.multiply(np.dot(esp_values, self.dist_y[k]), ESP.e_field_at_atom[1][k])
1138 4
                        self.Y2[beta] += R[k][beta] * np.multiply(np.dot(esp_values, self.dist_z[k]), ESP.e_field_at_atom[2][k])
1139

1140

1141 4
        self.Y = np.concatenate((self.Y1, self.Y2))
1142

1143

1144
    # noinspection PyProtectedMember
1145 4
    def build_vector_B(self):
1146
        """
1147
        Creates the Vector B for the charge fitting.
1148
        :return:
1149
        """
1150
        # Determine size of the vector for this molecule
1151
        # every atom is one line
1152
        # one line to restrain the overall charge of the molecule
1153
        # one line for every pair of equivalent atoms
1154 4
        self._lines_in_A = self.natoms + 1 + len(self._molecule.chemical_eq_atoms)
1155 4
        self.B = np.zeros(self._lines_in_A)
1156 4
        self.esp_values = self.baseESP.esp_values.to('elementary_charge / angstrom').magnitude
1157

1158 4
        for k in range(self.natoms):
1159 4
            self.B[k] = np.dot(self.esp_values, self.dist[k])
1160 4
        self.B[self.natoms] = self._molecule._charge
1161 4
        for polESP in self.polESPs:
1162 4
            esp_values = polESP.esp_values.to('elementary_charge / angstrom').magnitude
1163 4
            for k in range(self.natoms):
1164 4
                self.B[k] += np.dot(esp_values, self.dist[k])
1165 4
            self.B[self.natoms] = self._molecule._charge
1166

1167 4
        self.B = self.B / (len(self.polESPs) + 1)
1168 4
        for k, eq_atoms in enumerate(self._molecule.chemical_eq_atoms):
1169 4
            if eq_atoms[1] > 0:
1170 4
                self.B[self.natoms + 1 + k] = 0.0
1171

1172 4
    def build_vector_C(self):
1173
        """
1174
        Vector C corresponds to matrix D and is only for pur polarization fitting.
1175

1176
        :return:
1177
        """
1178 4
        self.C = np.zeros(self.Dlines)
1179 4
        self.esp_values = self.baseESP.esp_values.to('elementary_charge / angstrom').magnitude
1180 4
        for k in range(self.natoms):
1181 4
            self.C[k] = np.multiply(np.dot(self.esp_values, self.dist_x[k]), self.e_field_at_atom[0][k])
1182 4
            self.C[k + self.natoms] = np.multiply(np.dot(self.esp_values, self.dist_y[k]), self.e_field_at_atom[1][k])
1183 4
            self.C[k + self.natoms * 2] = np.multiply(np.dot(self.esp_values, self.dist_z[k]),
1184
                                                      self.e_field_at_atom[2][k])
1185

1186 4
        for polESP in self.polESPs:
1187 4
            esp_values = polESP.esp_values.to('elementary_charge / angstrom').magnitude
1188 4
            for k in range(self.natoms):
1189 4
                self.C[k] += np.multiply(np.dot(esp_values, self.dist_x[k]), polESP.e_field_at_atom[0][k])
1190 4
                self.C[k + self.natoms] += np.multiply(np.dot(esp_values, self.dist_y[k]), polESP.e_field_at_atom[1][k])
1191 4
                self.C[k + self.natoms * 2] += np.multiply(np.dot(esp_values, self.dist_z[k]),
1192
                                                           polESP.e_field_at_atom[2][k])
1193

1194 4
        self.C = self.C / (len(self.polESPs) + 1)
1195

1196
        #for j, atoms in enumerate(self._molecule.same_polarization_atoms):
1197
        #    self.C[5 * self.natoms + j] = 0.0
1198

1199 4
    def build_vector_Y(self):
1200
        """
1201
        Creates the Vector Y for the RESP-Pol Method.
1202
        :return:
1203
        """
1204 4
        self.build_vector_B()
1205 4
        self.build_vector_C()
1206 4
        self.Y = np.concatenate((self.B, self.C))
1207

1208 4
    def get_electric_field(self):
1209 4
        self.baseESP.get_electric_field()
1210 4
        self.e_field_at_atom = self.baseESP.e_field_at_atom
1211 4
        for polESP in self.polESPs:
1212 4
            polESP.get_electric_field()
1213

1214 4
    def optimize_charges_alpha(self):
1215
        """
1216
        Builds the necessary matrix and vector and performs a charge and polarizabilities optimization for this 1 conformation.
1217
        :return:
1218
        """
1219 0
        self.build_matrix_X()
1220 0
        self.build_vector_Y()
1221 0
        self.q_alpha = np.linalg.solve(self.X, self.Y)
1222

1223 4
    def optimize_charges_alpha_bcc(self):
1224
        """
1225
        Builds the necessary matrix and vector and performs a charge and polarizabilities optimization for this 1 conformation.
1226
        :return:
1227
        """
1228 0
        self.build_matrix_X_BCC()
1229 0
        self.build_vector_Y_BCC()
1230
        # Check if bccs or alphas are not in the training set and set them to zero
1231 0
        for i,row in enumerate(self.X):
1232 0
            if all(row == 0.0):
1233 0
                self.X[i][i]=1
1234

1235

1236 0
        self.q_alpha = np.linalg.solve(self.X, self.Y)
1237

1238

1239 4
    def get_distances(self):
1240 4
        self.get_grid_coord()
1241

1242
        # Distances between atoms and ESP points
1243 4
        self.dist = np.zeros((self.natoms, self.npoints))
1244 4
        self.dist_3 = np.zeros((self.natoms, self.npoints))
1245 4
        self.dist_x = np.zeros((self.natoms, self.npoints))
1246 4
        self.dist_y = np.zeros((self.natoms, self.npoints))
1247 4
        self.dist_z = np.zeros((self.natoms, self.npoints))
1248

1249 4
        self.dist = 1. / distance.cdist(self.atom_positions_angstrom, self.grid_coord_angstrom)
1250 4
        self.dist_3 = np.power(self.dist, 3)  # maybe free afterwards
1251 4
        self.dist_x = -np.multiply(
1252
            np.subtract.outer(np.transpose(self.atom_positions_angstrom)[0], np.transpose(self.grid_coord_angstrom)[0]),
1253
            self.dist_3)
1254
        # 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)
1255 4
        self.dist_y = -np.multiply(
1256
            np.subtract.outer(np.transpose(self.atom_positions_angstrom)[1], np.transpose(self.grid_coord_angstrom)[1]),
1257
            self.dist_3)
1258 4
        self.dist_z = -np.multiply(
1259
            np.subtract.outer(np.transpose(self.atom_positions_angstrom)[2], np.transpose(self.grid_coord_angstrom)[2]),
1260
            self.dist_3)
1261 4
        del self.dist_3
1262

1263
        # Distances between atoms and atoms
1264 4
        self.diatomic_dist = np.zeros((self.natoms, self.natoms))
1265 4
        self.diatomic_dist_3 = np.zeros((self.natoms, self.natoms))
1266 4
        self.diatomic_dist_5 = np.zeros((self.natoms, self.natoms))
1267 4
        self.diatomic_dist_x = np.zeros((self.natoms, self.natoms))
1268 4
        self.diatomic_dist_y = np.zeros((self.natoms, self.natoms))
1269 4
        self.diatomic_dist_z = np.zeros((self.natoms, self.natoms))
1270 4
        self.diatomic_distb_x = np.zeros((self.natoms, self.natoms))
1271 4
        self.diatomic_distb_y = np.zeros((self.natoms, self.natoms))
1272 4
        self.diatomic_distb_z = np.zeros((self.natoms, self.natoms))
1273

1274 4
        self.diatomic_dist = distance.cdist(self.atom_positions_angstrom, self.atom_positions_angstrom)
1275 4
        di = np.diag_indices(self.natoms)
1276 4
        self.diatomic_dist[di] = 1.0E10
1277
        # self.adist=np.fill_diagonal(self.adist,1.0)
1278 4
        self.diatomic_dist = 1. / self.diatomic_dist
1279 4
        self.diatomic_dist_3 = np.power(self.diatomic_dist, 3)
1280 4
        self.diatomic_dist_5 = np.power(self.diatomic_dist, 5)
1281 4
        self.diatomic_dist[di] = 0.0
1282 4
        self.diatomic_dist_x = np.multiply(np.subtract.outer(np.transpose(self.atom_positions_angstrom)[0],
1283
                                                             np.transpose(self.atom_positions_angstrom)[0]),
1284
                                           self.diatomic_dist_3)  # X distance between two atoms divided by the dist^3
1285 4
        self.diatomic_dist_y = np.multiply(np.subtract.outer(np.transpose(self.atom_positions_angstrom)[1],
1286
                                                             np.transpose(self.atom_positions_angstrom)[1]),
1287
                                           self.diatomic_dist_3)
1288 4
        self.diatomic_dist_z = np.multiply(np.subtract.outer(np.transpose(self.atom_positions_angstrom)[2],
1289
                                                             np.transpose(self.atom_positions_angstrom)[2]),
1290
                                           self.diatomic_dist_3)
1291 4
        self.diatomic_distb_x = np.subtract.outer(np.transpose(self.atom_positions_angstrom)[0],
1292
                                                  np.transpose(self.atom_positions_angstrom)[
1293
                                                      0])  # X distances between two atoms
1294 4
        self.diatomic_distb_y = np.subtract.outer(np.transpose(self.atom_positions_angstrom)[1],
1295
                                                  np.transpose(self.atom_positions_angstrom)[1])
1296 4
        self.diatomic_distb_z = np.subtract.outer(np.transpose(self.atom_positions_angstrom)[2],
1297
                                                  np.transpose(self.atom_positions_angstrom)[2])
1298

1299 4
    def delete_distances(self):
1300
        """Deletes the all calculated distances to free memory."""
1301 4
        del self.dist
1302 4
        del self.dist_x
1303 4
        del self.dist_y
1304 4
        del self.dist_z
1305

1306 4
        del self.diatomic_dist
1307 4
        del self.diatomic_dist_3
1308 4
        del self.diatomic_dist_5
1309 4
        del self.diatomic_dist_x
1310 4
        del self.diatomic_dist_y
1311 4
        del self.diatomic_dist_z
1312 4
        del self.diatomic_distb_x
1313 4
        del self.diatomic_distb_y
1314 4
        del self.diatomic_distb_z
1315

1316 4
    def write_res_esp(self, q_alpha= None ):
1317
        """
1318
        NOT FINISHED YET!!!
1319

1320
        Writes the residual ESP to a file.
1321

1322
        :param qd: list of float
1323
            Point charges and polarizabilities
1324
        :return:
1325
        """
1326 4
        if q_alpha is None:
1327 4
            q_alpha = self.q_alpha
1328 4
        atoms = []
1329 4
        for i,atom in enumerate(self.atom_positions.to('bohr').magnitude):
1330 4
            atoms.append([self._molecule._atoms[i]._atomic_number,atom])
1331

1332 4
        for ESP in [self.baseESP] + self.polESPs:
1333 4
            ESP.write_res_esp(q_alpha, atoms= atoms)
1334

1335
# =============================================================================================
1336
# ESPGRID
1337
# =============================================================================================
1338

1339
# noinspection PyTypeChecker
1340 4
class ESPGRID:
1341
    """
1342

1343
    """
1344

1345 4
    def define_grid(self, *args, **kwargs):
1346 4
        for ele in args:
1347 4
            if 'gesp' in ele:
1348 4
                self.gridtype = 'gesp'
1349 4
            if 'espf' in ele:
1350 0
                self.gridtype = 'respyte'
1351 4
            if 'grid.dat' in ele:
1352 0
                self.gridtype = 'psi4'
1353

1354 4
    def set_ext_e_field(self, vector):
1355 4
        self._ext_e_field = vector
1356

1357 4
    def load_grid(self, *args):
1358 4
        if self.gridtype == 'gesp':
1359 4
            self.name = args[0].rstrip('.gesp')
1360 4
            f = open(args[0], 'r')
1361 4
            lines = f.readlines()
1362 4
            f.close()
1363 4
            for i, line in enumerate(lines):
1364 4
                if 'ATOMIC' in line and 'COORDINATES' in line:
1365 4
                    self.natoms = int(line.strip('\n').split()[-1])
1366 4
                    for j in range(self.natoms):
1367 4
                        entry = lines[i + 1 + j].replace('D', 'E').split()
1368 4
                        self.atoms.append(entry[0])
1369 4
                        self.atom_positions.append(Q_([float(entry[k]) for k in range(1, 4, 1)], 'bohr'))
1370 4
                if 'GRID' in line:
1371 4
                    self.ngrid = int(line.strip('\n').split()[-1])
1372 4
                    grid = lines[i + 1:i + 1 + self.ngrid]
1373 4
                    break
1374
            # noinspection PyUnboundLocalVariable
1375 4
            for i, line in enumerate(grid):
1376 4
                grid[i] = [float(ele) for ele in line.replace('D', 'E').split()]
1377

1378 4
            self.positions = Q_(np.array(grid)[:, 1:4], 'bohr')
1379 4
            self.esp_values = Q_(np.array(grid)[:, 0], 'elementary_charge / bohr')
1380 0
        elif self.gridtype == 'respyte':
1381 0
            self.name = args[0].rstrip('.espf')
1382 0
            f = open(args[0], 'r')
1383 0
            lines = f.readlines()
1384 0
            f.close
1385 0
            ndata = int(len(lines) / 2) if len(lines) % 2 == 0 else int((len(lines) - 1) / 2)
1386 0
            grid = np.zeros((ndata, 4))
1387 0
            for i in range(ndata):
1388 0
                grid[i] = [float(ele) for ele in lines[2 * i].split()]
1389 0
            self.positions = Q_(np.array(grid)[:, 0:3], 'angstrom')
1390 0
            self.esp_values = Q_(np.array(grid)[:, 3], 'elementary_charge / bohr')
1391 0
        elif self.gridtype == 'psi4':
1392 0
            for ele in args:
1393 0
                if "grid.dat" in ele:
1394 0
                    gridfile = ele
1395 0
                elif 'esp.dat' in ele:
1396 0
                    espfile = ele
1397 0
            self.name = 'esp'
1398 0
            np.loadtxt(espfile)
1399 0
            self.positions = Q_(np.loadtxt(gridfile), 'angstrom')
1400 0
            self.esp_values = Q_(np.loadtxt(espfile), 'elementary_charge / bohr')
1401

1402 4
    def get_electric_field(self, ):
1403
        """
1404
        Calculates the electric field at every atomic positions.
1405
        :return:
1406
        """
1407 4
        alpha = self._conformer.q_alpha[self._conformer._lines_in_A:self._conformer._lines_in_A + 3*self._conformer.natoms]
1408 4
        alpha[np.where(alpha == 0.0)] += 10E-10
1409

1410
        # Load permanent charges for BCC method
1411
        # For all other methods this is set to 0.0 STILL HAVE to implment
1412 4
        try:
1413 4
            self._conformer.q_am1
1414 4
        except Exception:
1415 4
            log.warning('I do not have AM1-type charges')
1416 4
            self._conformer.q_am1 = np.zeros(self._conformer.natoms)
1417

1418 4
        for j in range(self._conformer.natoms):
1419 4
            self.e_field_at_atom[0][j] = np.dot(
1420
                np.multiply(self._conformer.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1421
                self._conformer.diatomic_dist_x[j]) + np.dot(
1422
                np.multiply(self._conformer.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1423
                self._conformer.diatomic_dist_x[j]) + self._ext_e_field[0].to(
1424
                'elementary_charge / angstrom / angstrom').magnitude
1425 4
            self.e_field_at_atom[1][j] = np.dot(
1426
                np.multiply(self._conformer.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1427
                self._conformer.diatomic_dist_y[j]) + np.dot(
1428
                np.multiply(self._conformer.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1429
                self._conformer.diatomic_dist_y[j]) + self._ext_e_field[1].to(
1430
                'elementary_charge / angstrom / angstrom').magnitude
1431 4
            self.e_field_at_atom[2][j] = np.dot(
1432
                np.multiply(self._conformer.q_alpha[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1433
                self._conformer.diatomic_dist_z[j]) + np.dot(
1434
                np.multiply(self._conformer.q_am1[:self._conformer.natoms], self._conformer._molecule.scale[j]),
1435
                self._conformer.diatomic_dist_z[j]) + self._ext_e_field[2].to(
1436
                'elementary_charge / angstrom / angstrom').magnitude
1437

1438 4
        self.e = self.e_field_at_atom.flatten()
1439
        #print(self.e)
1440 4
        if self._conformer._molecule._SCF and self._conformer._molecule.step > 0:
1441 4
            if not hasattr(self, 'Bdip') or self._conformer._molecule._thole:
1442 4
                if self._conformer._molecule._thole:
1443 4
                    thole_param=1.368711 # Change back to 1.368711
1444
                    #thole_param = 0.390
1445 4
                    dipole_tmp = np.where(alpha < 0.0, -alpha, alpha)
1446 4
                    thole_v = np.multiply(self._conformer.diatomic_dist, np.float_power(
1447
                        np.multiply(dipole_tmp[:self._conformer.natoms, None], dipole_tmp[:self._conformer.natoms]),
1448
                        1. / 6))
1449 4
                    di = np.diag_indices(self._conformer.natoms)
1450 4
                    thole_v[di] = 1.0
1451 4
                    thole_v = 1. / thole_v
1452 4
                    thole_v[di] = 0.0
1453

1454
                    # Exponential thole
1455 4
                    thole_fe = np.ones((self._conformer.natoms, self._conformer.natoms))
1456 4
                    thole_ft = np.ones((self._conformer.natoms, self._conformer.natoms))
1457 4
                    thole_fe -= np.exp(np.multiply(thole_param, np.power(-thole_v, 3)))
1458 4
                    thole_ft -= np.multiply(np.multiply(thole_param, np.power(thole_v, 3)) + 1.,
1459
                                            np.exp(np.multiply(thole_param, np.power(-thole_v, 3))))
1460
                    # 1.5 was found in the OpenMM code. Not sure whuy it is there
1461

1462
                    # In original thole these lines should not be here
1463
                    # thole_ft = np.multiply(thole_ft, self._conformer.scale)
1464
                    # thole_fe = np.multiply(thole_fe, self._conformer.scale)
1465
                    # Linear thole
1466
                    # thole_fe=np.zeros((self._conformer.natoms,self._conformer.natoms))
1467
                    # thole_fe=np.zeros((self._conformer.natoms,self._conformer.natoms))
1468
                    # thole_fe=np.where(thole_v>1.0,1.0,4*np.power(thole_v,3)-3*np.power(thole_v,4))
1469
                    # thole_ft=np.where(thole_v>1.0,1.0,np.power(thole_v,4))
1470
                else:
1471 4
                    try:
1472 4
                        thole_ft = self._conformer._molecule.scale_scf
1473 4
                        thole_fe = self._conformer._molecule.scale_scf
1474 0
                    except Exception:
1475

1476 0
                        thole_ft = self._conformer._molecule.scale
1477 0
                        thole_fe = self._conformer._molecule.scale
1478
                    else:
1479 4
                        print('Using different set of scaling for SCF interactions')
1480 4
                        log.info('Using different set of scaling for SCF interactions')
1481 4
                Bdip11 = np.add(np.multiply(thole_fe, self._conformer.diatomic_dist_3), np.multiply(thole_ft,
1482
                                                                                                    -3 * np.multiply(
1483
                                                                                                        np.multiply(
1484
                                                                                                            self._conformer.diatomic_distb_x,
1485
                                                                                                            self._conformer.diatomic_distb_x),
1486
                                                                                                        self._conformer.diatomic_dist_5)))
1487 4
                Bdip22 = np.add(np.multiply(thole_fe, self._conformer.diatomic_dist_3), np.multiply(thole_ft,
1488
                                                                                                    -3 * np.multiply(
1489
                                                                                                        np.multiply(
1490
                                                                                                            self._conformer.diatomic_distb_y,
1491
                                                                                                            self._conformer.diatomic_distb_y),
1492
                                                                                                        self._conformer.diatomic_dist_5)))
1493 4
                Bdip33 = np.add(np.multiply(thole_fe, self._conformer.diatomic_dist_3), np.multiply(thole_ft,
1494
                                                                                                    -3 * np.multiply(
1495
                                                                                                        np.multiply(
1496
                                                                                                            self._conformer.diatomic_distb_z,
1497
                                                                                                            self._conformer.diatomic_distb_z),
1498
                                                                                                        self._conformer.diatomic_dist_5)))
1499 4
                Bdip12 = np.multiply(thole_ft,
1500
                                     -3 * np.multiply(np.multiply(self._conformer.diatomic_distb_x,
1501
                                                                  self._conformer.diatomic_distb_y),
1502
                                                      self._conformer.diatomic_dist_5))
1503 4
                Bdip13 = np.multiply(thole_ft,
1504
                                     -3 * np.multiply(np.multiply(self._conformer.diatomic_distb_x,
1505
                                                                  self._conformer.diatomic_distb_z),
1506
                                                      self._conformer.diatomic_dist_5))
1507 4
                Bdip23 = np.multiply(thole_ft,
1508
                                     -3 * np.multiply(np.multiply(self._conformer.diatomic_distb_y,
1509
                                                                  self._conformer.diatomic_distb_z),
1510
                                                      self._conformer.diatomic_dist_5))
1511 4
                Bdip = np.concatenate((np.concatenate((Bdip11, Bdip12, Bdip13), axis=1),
1512
                                       np.concatenate((Bdip12, Bdip22, Bdip23), axis=1),
1513
                                       np.concatenate((Bdip13, Bdip23, Bdip33), axis=1)),
1514
                                      axis=0)
1515

1516 4
            for j in range(self._conformer.natoms):
1517 4
                for k in range(3):
1518 4
                    for l in range(3):
1519 4
                        Bdip[k * self._conformer.natoms + j][l * self._conformer.natoms + j] = 0.0
1520

1521 4
            for j in range(3 * self._conformer.natoms):
1522 4
                Bdip[j][j] = 1. / alpha[j]
1523 4
            dipole_scf = np.linalg.solve(Bdip, self.e)
1524 4
            self.e = np.divide(dipole_scf, alpha)
1525 4
            print(self.e)
1526 4
        self.e_field_at_atom[0] = 1.0 * self.e[:self._conformer.natoms]
1527 4
        self.e_field_at_atom[1] = 1.0 * self.e[self._conformer.natoms:2 * self._conformer.natoms]
1528 4
        self.e_field_at_atom[2] = 1.0 * self.e[2 * self._conformer.natoms:3 * self._conformer.natoms]
1529

1530
    # Analysation stuff
1531 4
    def calc_esp_q_alpha(self, q_alpha):
1532
        """
1533
        Calculates the ESP for a given set of point charges and polarizabilities.
1534

1535
        :param qd: list of float
1536
            Point charges and polarizabilities
1537
        :return:
1538

1539
        Calculates the ESP on every point of the initially read GESP file
1540
        """
1541 4
        self.q_pot = np.zeros(len(self.esp_values))
1542 4
        for i in range(len(self.esp_values)):
1543 4
            if self._conformer._molecule._mode == 'q':
1544 0
                self.q_pot[i] = np.dot(q_alpha[:self._conformer.natoms], np.transpose(self._conformer.dist)[i])
1545
            else:
1546 4
                dipole = q_alpha[self._conformer._lines_in_A:]
1547
                # self.dipole[self.dipole<0.0]=0.0
1548 4
                try:
1549
                    # 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])
1550 4
                    e_dip = np.dot(np.multiply(dipole[:self.natoms], self.e_field_at_atom[0]),
1551
                                   np.transpose(self._conformer.dist_x)[i]) + np.dot(
1552
                        np.multiply(dipole[self.natoms:2 * self.natoms], self.e_field_at_atom[1]),
1553
                        np.transpose(self._conformer.dist_y)[i]) + np.dot(
1554
                        np.multiply(dipole[2 * self.natoms:3 * self.natoms], self.e_field_at_atom[2]),
1555
                        np.transpose(self._conformer.dist_z)[i])
1556 4
                    self.q_pot[i] = np.dot(q_alpha[:self.natoms], np.transpose(self._conformer.dist)[i]) + e_dip
1557 0
                except:
1558 0
                    self.q_pot[i] = np.dot(q_alpha[:self.natoms], np.transpose(self._conformer.dist)[i])
1559
                    # print('DEBUG no electric field ')
1560
            # if self.mode=='d':
1561
            #    self.dipole=qd
1562
            #    e_dip=np.dot(np.multiply(self.dipole[:self.natoms],self.e_x),np.transpose(self.dist_x)[i])+np.dot(np.multiply(self.dipole[self.natoms:2*self.natoms],self.e_y),np.transpose(self.dist_y)[i])+np.dot(np.multiply(self.dipole[2*self.natoms:3*self.natoms],self.e_z),np.transpose(self.dist_z)[i])
1563
            #    self.q_pot[i]=e_dip
1564

1565 4
    def sub_esp_q_alpha(self, q_alpha):
1566
        """
1567
        Subtracts the ESP create by a set of point charges and polarizabilities.
1568

1569
        :param qd: list of float
1570
            Point charges and polarizabilities
1571
        :return:
1572

1573
        """
1574 4
        if self._conformer._molecule._mode != 'alpha':  # Change that in bcc that this line is unnecessary
1575 4
            self.calc_esp_q_alpha(q_alpha)
1576 4
            self.esp_values = Q_(np.subtract(self.esp_values.to('elementary_charge / angstrom').magnitude, self.q_pot),'elementary_charge / angstrom')
1577

1578 4
    def calc_sse(self, q_alpha):
1579
        """
1580
        Calculate the Squared Sum of Errors between the stored ESP and a ESP created from qd.
1581
        :param qd: list of float
1582
            Point charges and polarizabilities
1583
        :return:
1584
        """
1585 4
        self.calc_esp_q_alpha(q_alpha)
1586 4
        self.sse = np.square(self.esp_values.to('elementary_charge / angstrom').magnitude - self.q_pot).sum()
1587

1588 4
    def write_res_esp(self, q_alpha, atoms = []):
1589
        """
1590
        NOT FINISHED YET!!!
1591

1592
        Writes the residual ESP to a file.
1593

1594
        :param qd: list of float
1595
            Point charges and polarizabilities
1596
        :return:
1597
        """
1598 4
        self.calc_esp_q_alpha(q_alpha)
1599 4
        res_pot = np.subtract(self.esp_values.to('elementary_charge / angstrom').magnitude, Q_(self.q_pot, 'elementary_charge / angstrom'))
1600
        #res_pot = (self.esp_values - self.q_pot)#.to('elementary_charge / bohr').magnitude
1601 4
        f = open(self.name + '.rgesp', 'w')
1602 4
        f.write(' ESP FILE - ATOMIC UNITS\n')
1603 4
        f.write(' CHARGE =  {0} - MULTIPLICITY =   1\n'.format(self._conformer._molecule._charge))
1604 4
        f.write(' ATOMIC COORDINATES AND ESP CHARGES. #ATOMS =     {} \n'.format(np.sum(self.natoms)))
1605 4
        for i in range(self._conformer.natoms):
1606 4
            f.write(
1607
                ' {} {} {} {} {}\n'.format(atoms[i][0], atoms[i][1][0], atoms[i][1][1], atoms[i][1][2],
1608
                                           self._conformer._molecule.q_alpha[i]))
1609 4
        f.write(' ESP VALUES AND GRID POINT COORDINATES. #POINTS =   {}\n'.format(len(self.esp_values)))
1610 4
        for i in range(len(self.esp_values)):
1611 4
            try:
1612 4
                f.write(' {} {} {} {}\n'.format(res_pot[i].magnitude, self.positions[i][0].magnitude, self.positions[i][1].magnitude, self.positions[i][2].magnitude))
1613 2
            except: # Difference between python 3.6 and 3.7
1614 2
                f.write(' {} {} {} {}\n'.format(res_pot[i], self.positions[i][0].magnitude, self.positions[i][1].magnitude, self.positions[i][2].magnitude))
1615 4
        f.close()
1616

1617
# =============================================================================================
1618
# BCCUnpolESP
1619
# =============================================================================================
1620

1621 4
class BCCUnpolESP(ESPGRID):
1622
    """
1623

1624
    """
1625

1626 4
    def __init__(self, *args, conformer=None):
1627
        # Decide if we have a Gaussian grid or a psi4 grid
1628 4
        self.gridtype = None
1629 4
        self.natoms = -1
1630 4
        self.atoms = []
1631 4
        self.atom_positions = []
1632 4
        self.define_grid(*args)
1633 4
        self.esp_values = None
1634 4
        self.positions = None
1635 4
        self._conformer = conformer
1636

1637
        # External e-field is 0 in all directions
1638 4
        vector = Q_([0, 0, 0], 'elementary_charge / bohr / bohr')
1639 4
        self.set_ext_e_field(vector)
1640

1641 4
        self.load_grid(*args)
1642

1643 4
        self.e_field_at_atom = np.zeros((3, self._conformer.natoms))
1644

1645

1646
# =============================================================================================
1647
# BCCPolESP
1648
# =============================================================================================
1649

1650 4
class BCCPolESP(ESPGRID):
1651
    """
1652

1653
    """
1654

1655 4
    def __init__(self, *args, conformer=None, e_field=Q_([0.0, 0.0, 0.0], 'elementary_charge / bohr / bohr')):
1656
        # Decide if we have a Gaussian grid or a psi4 grid
1657 4
        self.gridtype = None
1658 4
        self.natoms = -1
1659 4
        self.atoms = []
1660 4
        self.atom_positions = []
1661 4
        self.define_grid(*args)
1662 4
        self.esp_values = None
1663 4
        self.positions = None
1664 4
        self._conformer = conformer
1665

1666
        # Set external e-field
1667 4
        self.set_ext_e_field(e_field)
1668

1669 4
        self.load_grid(*args)
1670

1671 4
        self.e_field_at_atom = np.zeros((3, self._conformer.natoms))
1672

1673

1674 4
if __name__ == '__main__':
1675
    pass
1676 4
    """
1677
    datei = os.path.join(ROOT_DIR_PATH, 'resppol/tmp/phenol/conf0/mp2_0.mol2')
1678
    test = TrainingSet(scf_scaleparameters=[0.0, 0.0, 0.5])
1679
    test.add_molecule(datei)
1680
    test.molecules[0].add_conformer_from_mol2(datei)
1681
    print(test.molecules[0].same_polarization_atoms)
1682
    print(test.molecules[0].scale)
1683
    print(test.molecules[0].scale_scf)
1684
    espfile = '/home/michael/resppol/resppol/tmp/phenol/conf0/molecule0.gesp'
1685
    test.molecules[0].conformers[0].add_baseESP(espfile)
1686
    datei = os.path.join(ROOT_DIR_PATH, 'resppol/tmp/butanol/conf0/mp2_0.mol2')
1687
    test.add_molecule(datei)
1688
    test.molecules[1].add_conformer_from_mol2(datei)
1689
    espfile = '/home/michael/resppol/resppol/tmp/butanol/conf0/molecule0.gesp'
1690
    test.molecules[1].conformers[0].add_baseESP(espfile)
1691
    test.optimize_charges()
1692
    test.molecules[0].conformers[0].build_matrix_X()
1693
    for molecule in test.molecules:
1694
        print(molecule.q)
1695

1696
    print('FINISH')
1697
    
1698
    datei = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test2.mol2')
1699
    test = TrainingSet(mode='q_alpha',SCF= True,scf_scaleparameters=[1,1,1], scaleparameters=[1,1,1])
1700
    #test = TrainingSet(mode='q_alpha')
1701
    test.add_molecule(datei)
1702
    test.molecules[0].add_conformer_from_mol2(datei)
1703
    #test.molecules[0].add_conformer_from_mol2(datei)
1704
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3.gesp')
1705
    test.molecules[0].conformers[0].add_baseESP(espfile)
1706
    #test.molecules[0].conformers[1].add_baseESP(espfile)
1707
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3_Z+.gesp')
1708
    test.molecules[0].conformers[0].add_polESP(espfile, e_field=Q_([0.0, 0.0, 1.], 'elementary_charge / bohr / bohr'))
1709
    #test.molecules[0].conformers[1].add_polESP(espfile, e_field=Q_([0.0, 0.0, 0.8],'elementary_charge / bohr / bohr'))
1710
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3_Z-.gesp')
1711
    test.molecules[0].conformers[0].add_polESP(espfile, e_field=Q_([0.0, 0.0, -1.], 'elementary_charge / bohr / bohr') )
1712
    #test.molecules[0].conformers[1].add_polESP(espfile, e_field=Q_([0.0, 0.0, -0.8], 'elementary_charge / bohr / bohr') )
1713
    #test.build_matrix_X()
1714
    #test.build_vector_Y()
1715
    test.build_matrix_X_BCC()
1716
    test.build_vector_Y_BCC()
1717
    test.optimize_bcc_alpha()
1718
    print(test.molecules[0].conformers[0].q_alpha)
1719
    #test.optimize_charges_alpha()
1720
    #print(test.q_alpha)
1721
    print(test.molecules[0].conformers[0].baseESP.e_field_at_atom)
1722
    print(test.molecules[0].conformers[0].polESPs[0].e_field_at_atom)
1723
    
1724

1725
    datei = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test2.mol2')
1726
    test = TrainingSet(mode='q_alpha',SCF= True, thole = True)
1727
    test.add_molecule(datei)
1728
    test.molecules[0].add_conformer_from_mol2(datei)
1729
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3.gesp')
1730
    test.molecules[0].conformers[0].add_baseESP(espfile)
1731
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3_Z+.gesp')
1732
    test.molecules[0].conformers[0].add_polESP(espfile, e_field=Q_([0.0, 0.0, 1.0], 'elementary_charge / bohr / bohr'))
1733
    espfile = os.path.join(ROOT_DIR_PATH, 'resppol/data/fast_test_data/test3_Z-.gesp')
1734
    test.molecules[0].conformers[0].add_polESP(espfile, e_field=Q_([0.0, 0.0, -1.0], 'elementary_charge / bohr / bohr') )
1735
    test.optimize_charges_alpha()
1736
    test.molecules[0].conformers[0].baseESP.calc_esp_q_alpha(test.q_alpha)
1737
    test.molecules[0].conformers[0].baseESP.calc_sse(test.q_alpha)
1738
    test.molecules[0].conformers[0].baseESP.sub_esp_q_alpha(test.q_alpha)
1739
    test.molecules[0].conformers[0].write_res_esp()
1740
    print(test.molecules[0].conformers[0].baseESP.q_pot)
1741

1742
    print(test.molecules[0].conformers[0].q_alpha)
1743
    """

Read our documentation on viewing source code .

Loading