|
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
|
|
"""
|
MSchauperl