1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
|
import numpy as np
from ase import units
from ase.calculators.calculator import Calculator, all_changes
from ase.calculators.tip3p import TIP3P, angleHOH, rOH
__all__ = ['rOH', 'angleHOH', 'TIP4P', 'sigma0', 'epsilon0']
# Electrostatic constant and parameters:
k_c = units.Hartree * units.Bohr
qH = 0.52
A = 600e3 * units.kcal / units.mol
B = 610 * units.kcal / units.mol
sigma0 = (A / B)**(1 / 6.)
epsilon0 = B**2 / (4 * A)
# https://doi.org/10.1063/1.445869
class TIP4P(TIP3P):
def __init__(self, rc=7.0, width=1.0):
""" TIP4P potential for water.
:doi:`10.1063/1.445869`
Requires an atoms object of OHH,OHH, ... sequence
Correct TIP4P charges and LJ parameters set automatically.
Virtual interaction sites implemented in the following scheme:
Original atoms object has no virtual sites.
When energy/forces are requested:
* virtual sites added to temporary xatoms object
* energy / forces calculated
* forces redistributed from virtual sites to actual atoms object
This means you do not get into trouble when propagating your system
with MD while having to skip / account for massless virtual sites.
This also means that if using for QM/MM MD with GPAW, the EmbedTIP4P
class must be used.
"""
TIP3P.__init__(self, rc, width)
self.atoms_per_mol = 3
self.sites_per_mol = 4
self.energy = None
self.forces = None
def calculate(self, atoms=None,
properties=['energy', 'forces'],
system_changes=all_changes):
Calculator.calculate(self, atoms, properties, system_changes)
assert (atoms.numbers[::3] == 8).all()
assert (atoms.numbers[1::3] == 1).all()
assert (atoms.numbers[2::3] == 1).all()
xpos = self.add_virtual_sites(atoms.positions)
xcharges = self.get_virtual_charges(atoms)
cell = atoms.cell
pbc = atoms.pbc
natoms = len(atoms)
nmol = natoms // 3
self.energy = 0.0
self.forces = np.zeros((4 * natoms // 3, 3))
C = cell.diagonal()
assert (cell == np.diag(C)).all(), 'not orthorhombic'
assert ((C >= 2 * self.rc) | ~pbc).all(), 'cutoff too large'
# Get dx,dy,dz from first atom of each mol to same atom of all other
# and find min. distance. Everything moves according to this analysis.
for a in range(nmol - 1):
D = xpos[(a + 1) * 4::4] - xpos[a * 4]
shift = np.zeros_like(D)
for i, periodic in enumerate(pbc):
if periodic:
shift[:, i] = np.rint(D[:, i] / C[i]) * C[i]
q_v = xcharges[(a + 1) * 4:]
# Min. img. position list as seen for molecule !a!
position_list = np.zeros(((nmol - 1 - a) * 4, 3))
for j in range(4):
position_list[j::4] += xpos[(a + 1) * 4 + j::4] - shift
# Make the smooth cutoff:
pbcRoo = position_list[::4] - xpos[a * 4]
pbcDoo = np.sum(np.abs(pbcRoo)**2, axis=-1)**(1 / 2)
x1 = pbcDoo > self.rc - self.width
x2 = pbcDoo < self.rc
x12 = np.logical_and(x1, x2)
y = (pbcDoo[x12] - self.rc + self.width) / self.width
t = np.zeros(len(pbcDoo))
t[x2] = 1.0
t[x12] -= y**2 * (3.0 - 2.0 * y)
dtdd = np.zeros(len(pbcDoo))
dtdd[x12] -= 6.0 / self.width * y * (1.0 - y)
self.energy_and_forces(a, xpos, position_list, q_v, nmol, t, dtdd)
if self.pcpot:
e, f = self.pcpot.calculate(xcharges, xpos)
self.energy += e
self.forces += f
f = self.redistribute_forces(self.forces)
self.results['energy'] = self.energy
self.results['forces'] = f
def energy_and_forces(self, a, xpos, position_list, q_v, nmol, t, dtdd):
""" energy and forces on molecule a from all other molecules.
cutoff is based on O-O Distance. """
# LJ part - only O-O interactions
epsil = np.tile([epsilon0], nmol - 1 - a)
sigma = np.tile([sigma0], nmol - 1 - a)
DOO = position_list[::4] - xpos[a * 4]
d2 = (DOO**2).sum(1)
d = np.sqrt(d2)
e_lj = 4 * epsil * (sigma**12 / d**12 - sigma**6 / d**6)
f_lj = (4 * epsil * (12 * sigma**12 / d**13 -
6 * sigma**6 / d**7) * t -
e_lj * dtdd)[:, np.newaxis] * DOO / d[:, np.newaxis]
self.forces[a * 4] -= f_lj.sum(0)
self.forces[(a + 1) * 4::4] += f_lj
# Electrostatics
e_elec = 0
all_cut = np.repeat(t, 4)
for i in range(4):
D = position_list - xpos[a * 4 + i]
d2_all = (D**2).sum(axis=1)
d_all = np.sqrt(d2_all)
e = k_c * q_v[i] * q_v / d_all
e_elec += np.dot(all_cut, e).sum()
e_f = e.reshape(nmol - a - 1, 4).sum(1)
F = (e / d_all * all_cut)[:, np.newaxis] * D / d_all[:, np.newaxis]
FOO = -(e_f * dtdd)[:, np.newaxis] * DOO / d[:, np.newaxis]
self.forces[(a + 1) * 4 + 0::4] += FOO
self.forces[a * 4] -= FOO.sum(0)
self.forces[(a + 1) * 4:] += F
self.forces[a * 4 + i] -= F.sum(0)
self.energy += np.dot(e_lj, t) + e_elec
def add_virtual_sites(self, pos):
# Order: OHHM,OHHM,...
# DOI: 10.1002/(SICI)1096-987X(199906)20:8
b = 0.15
xatomspos = np.zeros((4 * len(pos) // 3, 3))
for w in range(0, len(pos), 3):
r_i = pos[w] # O pos
r_j = pos[w + 1] # H1 pos
r_k = pos[w + 2] # H2 pos
n = (r_j + r_k) / 2 - r_i
n /= np.linalg.norm(n)
r_d = r_i + b * n
x = 4 * w // 3
xatomspos[x + 0] = r_i
xatomspos[x + 1] = r_j
xatomspos[x + 2] = r_k
xatomspos[x + 3] = r_d
return xatomspos
def get_virtual_charges(self, atoms):
charges = np.empty(len(atoms) * 4 // 3)
charges[0::4] = 0.00 # O
charges[1::4] = qH # H1
charges[2::4] = qH # H2
charges[3::4] = - 2 * qH # X1
return charges
def redistribute_forces(self, forces):
f = forces
b = 0.15
a = 0.5
pos = self.atoms.positions
for w in range(0, len(pos), 3):
r_i = pos[w] # O pos
r_j = pos[w + 1] # H1 pos
r_k = pos[w + 2] # H2 pos
r_ij = r_j - r_i
r_jk = r_k - r_j
r_d = r_i + b * (r_ij + a * r_jk) / np.linalg.norm(r_ij + a * r_jk)
r_id = r_d - r_i
gamma = b / np.linalg.norm(r_ij + a * r_jk)
x = w * 4 // 3
Fd = f[x + 3] # force on M
F1 = (np.dot(r_id, Fd) / np.dot(r_id, r_id)) * r_id
Fi = Fd - gamma * (Fd - F1) # Force from M on O
Fj = (1 - a) * gamma * (Fd - F1) # Force from M on H1
Fk = a * gamma * (Fd - F1) # Force from M on H2
f[x] += Fi
f[x + 1] += Fj
f[x + 2] += Fk
# remove virtual sites from force array
f = np.delete(f, list(range(3, f.shape[0], 4)), axis=0)
return f
|
