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
|
# fmt: off
from math import cos, pi, sin
import numpy as np
from ase import Atoms, units
from ase.calculators.combine_mm import CombineMM
from ase.calculators.counterions import AtomicCounterIon as ACI
from ase.calculators.fd import calculate_numerical_forces
from ase.calculators.qmmm import LJInteractionsGeneral
from ase.calculators.tip3p import TIP3P, angleHOH, rOH
from ase.calculators.tip3p import epsilon0 as eps3
from ase.calculators.tip3p import sigma0 as sig3
from ase.calculators.tip4p import TIP4P
from ase.calculators.tip4p import epsilon0 as eps4
from ase.calculators.tip4p import sigma0 as sig4
def test_combine_mm():
"""Test CombineMM forces by combining tip3p and tip4p with them selves, and
by combining tip3p with tip4p and testing against numerical forces.
Also test LJInterationsGeneral with CombineMM """
def make_atoms():
r = rOH
a = angleHOH * pi / 180
dimer = Atoms('H2OH2O',
[(r * cos(a), 0, r * sin(a)),
(r, 0, 0),
(0, 0, 0),
(r * cos(a / 2), r * sin(a / 2), 0),
(r * cos(a / 2), -r * sin(a / 2), 0),
(0, 0, 0)])
dimer = dimer[[2, 0, 1, 5, 3, 4]]
# put O-O distance in the cutoff range
dimer.positions[3:, 0] += 2.8
return dimer
dimer = make_atoms()
rc = 3.0
for (TIPnP, (eps, sig), nm) in zip([TIP3P, TIP4P],
((eps3, sig3), (eps4, sig4)),
[3, 3]):
dimer.calc = TIPnP(rc=rc, width=1.0)
F1 = dimer.get_forces()
sigma = np.array([sig, 0, 0])
epsilon = np.array([eps, 0, 0])
dimer.calc = CombineMM([0, 1, 2], nm, nm,
TIPnP(rc=rc, width=1.0),
TIPnP(rc=rc, width=1.0),
sigma, epsilon, sigma, epsilon,
rc=rc, width=1.0)
F2 = dimer.get_forces()
dF = F1 - F2
print(TIPnP)
print(dF)
assert abs(dF).max() < 1e-8
# Also check a TIP3P/TIP4P combination against numerical forces:
eps1 = np.array([eps3, 0, 0])
sig1 = np.array([sig3, 0, 0])
eps2 = np.array([eps4, 0, 0])
sig2 = np.array([sig4, 0, 0])
dimer.calc = CombineMM([0, 1, 2], 3, 3, TIP3P(rc, 1.0), TIP4P(rc, 1.0),
sig1, eps1, sig2, eps2, rc, 1.0)
F2 = dimer.get_forces()
Fn = calculate_numerical_forces(dimer, 1e-7)
dF = F2 - Fn
print('TIP3P/TIP4P')
print(dF)
assert abs(dF).max() < 1e-8
# LJInteractionsGeneral with CombineMM.
# As it is used within EIQMMM, but avoiding long calculations for tests.
# Total system is a unit test system comprised of:
# 2 Na ions playing the role of a 'QM' subsystem
# 2 Na ions as the 'Counterions'
# 2 Water molecules
dimer.calc = []
faux_qm = Atoms('2Na', positions=np.array([[1.4, -4, 0], [1.4, 4, 0]]))
ions = Atoms('2Na', positions=np.array([[1.4, 0, -4], [1.4, 0, 4]]))
mmatoms = ions + dimer
sigNa = 1.868 * (1.0 / 2.0)**(1.0 / 6.0) * 10
epsNa = 0.00277 * units.kcal / units.mol
# ACI for atoms 0 and 1 of the MM subsystem (2 and 3 for the total system)
# 1 atom 'per molecule'. The rest is TIP4P, 3 atoms per molecule:
calc = CombineMM([0, 1], 1, 3,
ACI(1, epsNa, sigNa), # calc 1
TIP4P(), # calc 2
[sigNa], [epsNa], # LJs for subsystem 1
sig2, eps2, # arrays for TIP4P from earlier
rc=7.5)
mmatoms.calc = calc
mmatoms.calc.initialize(mmatoms)
# LJ arrays for the 'QM' subsystem
sig_qm = np.array([sigNa, sigNa])
eps_qm = np.array([epsNa, epsNa])
# For the MM subsystem, tuple of arrays (counterion, water)
sig_mm = (np.array([sigNa]), sig2)
eps_mm = (np.array([epsNa]), eps2)
lj = LJInteractionsGeneral(sig_qm, eps_qm, sig_mm, eps_mm, 2)
ecomb, _fcomb1, _fcomb2 = lj.calculate(faux_qm,
mmatoms,
np.array([0, 0, 0]))
# This should give the same result as if not using CombineMM, on a sum
# of these systems:
# A: All the Na atoms in the 'QM' region
# B: LJInteractions between the 'QM' Na and the 'MM' Na
# C: -LJinteractions between the 'MM' Na and the 'MM' Water
# A:
mmatoms = dimer
mmatoms.calc = []
sig_qm = np.concatenate((sig_qm, sig_qm))
eps_qm = np.concatenate((eps_qm, eps_qm))
sig_mm = sig_mm[1]
eps_mm = eps_mm[1]
lj = LJInteractionsGeneral(sig_qm, eps_qm, sig_mm, eps_mm, 4)
ea, _fa1, _fa2 = lj.calculate(faux_qm + ions, mmatoms, np.array([0, 0, 0]))
# B:
lj = LJInteractionsGeneral(sig_qm[:2], eps_qm[:2],
sig_qm[:2], eps_qm[:2], 2, 2)
eb, _fb1, _fb2, = lj.calculate(faux_qm, ions, np.array([0, 0, 0]))
# C:
lj = LJInteractionsGeneral(sig_qm[:2], eps_qm[:2],
sig_mm, eps_mm,
2, 3)
ec, _fc1, _fc2, = lj.calculate(ions, dimer, np.array([0, 0, 0]))
assert ecomb - (ea + eb) + ec == 0
|
