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