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
|
from ase.calculators.emt import EMT
from ase.constraints import FixInternals
from ase.optimize.bfgs import BFGS
from ase.build import molecule
import copy
import pytest
# Convenience functions to compute linear combinations of internal coordinates
def get_bondcombo(atoms, bondcombo_def):
return sum([defin[2] * atoms.get_distance(*defin[0:2]) for
defin in bondcombo_def])
def get_anglecombo(atoms, anglecombo_def):
return sum([defin[3] * atoms.get_angle(*defin[0:3]) for
defin in anglecombo_def])
def get_dihedralcombo(atoms, dihedralcombo_def):
return sum([defin[4] * atoms.get_dihedral(*defin[0:4]) for
defin in dihedralcombo_def])
def setup_atoms():
atoms = molecule('CH3CH2OH', vacuum=5.0)
atoms.rattle(stdev=0.3)
return atoms
def setup_fixinternals():
atoms = setup_atoms()
# Angles, Bonds, Dihedrals are built up with pairs of constraint
# value and indices defining the constraint
# Linear combinations of bond lengths are built up similarly with the
# coefficients appended to the indices defining the constraint
# Fix bond between atoms 1 and 2 to 1.4
bond_def = [1, 2]
target_bond = 1.4
# Fix angle to whatever it was from the start
angle_def = [6, 0, 1]
target_angle = atoms.get_angle(*angle_def)
# Fix this dihedral angle to whatever it was from the start
dihedral_def = [6, 0, 1, 2]
target_dihedral = atoms.get_dihedral(*dihedral_def)
# Initialize constraint
constr = FixInternals(bonds=[(target_bond, bond_def)],
angles_deg=[(target_angle, angle_def)],
dihedrals_deg=[(target_dihedral, dihedral_def)],
epsilon=1e-10)
print(constr)
return (atoms, constr, bond_def, target_bond, angle_def, target_angle,
dihedral_def, target_dihedral)
def test_fixinternals():
(atoms, constr, bond_def, target_bond, angle_def, target_angle,
dihedral_def, target_dihedral) = setup_fixinternals()
calc = EMT()
opt = BFGS(atoms)
previous_angle = atoms.get_angle(*angle_def)
previous_dihedral = atoms.get_dihedral(*dihedral_def)
print('angle before', previous_angle)
print('dihedral before', previous_dihedral)
print('bond length before', atoms.get_distance(*bond_def))
print('target bondlength', target_bond)
atoms.calc = calc
atoms.set_constraint(constr)
print('-----Optimization-----')
opt.run(fmax=0.01)
new_angle = atoms.get_angle(*angle_def)
new_dihedral = atoms.get_dihedral(*dihedral_def)
new_bondlength = atoms.get_distance(*bond_def)
print('angle after', new_angle)
print('dihedral after', new_dihedral)
print('bondlength after', new_bondlength)
err1 = new_angle - previous_angle
err2 = new_dihedral - previous_dihedral
err3 = new_bondlength - target_bond
print('error in angle', repr(err1))
print('error in dihedral', repr(err2))
print('error in bondlength', repr(err3))
assert err1 < 1e-11
assert err2 < 1e-12
assert err3 < 1e-12
def setup_combos():
atoms = setup_atoms()
# Fix linear combination of two bond lengths with atom indices 0-8 and
# 0-6 with weighting coefficients 1.0 and -1.0 to the current value.
# In other words, fulfil the following constraint:
# 1.0 * atoms.get_distance(2, 1) + -1.0 * atoms.get_distance(2, 3) = const.
bondcombo_def = [[2, 1, 1.0], [2, 3, -1.0]]
target_bondcombo = get_bondcombo(atoms, bondcombo_def)
# Fix linear combination of two angles
# 1. * atoms.get_angle(7, 0, 8) + 1. * atoms.get_angle(7, 0, 6) = const.
anglecombo_def = [[7, 0, 8, 1.], [7, 0, 6, 1]]
target_anglecombo = get_anglecombo(atoms, anglecombo_def)
# Fix linear combination of two dihedrals
dihedralcombo_def = [[3, 2, 1, 4, 1.0], [2, 1, 0, 7, 1.0]]
target_dihedralcombo = get_dihedralcombo(atoms, dihedralcombo_def)
# Initialize constraint
constr = FixInternals(bondcombos=[(target_bondcombo, bondcombo_def)],
anglecombos=[(target_anglecombo, anglecombo_def)],
dihedralcombos=[(target_dihedralcombo,
dihedralcombo_def)], epsilon=1e-10)
print(constr)
return (atoms, constr, bondcombo_def, target_bondcombo, anglecombo_def,
target_anglecombo, dihedralcombo_def, target_dihedralcombo)
def test_combos():
(atoms, constr, bondcombo_def, target_bondcombo, anglecombo_def,
target_anglecombo, dihedralcombo_def,
target_dihedralcombo) = setup_combos()
ref_bondcombo = get_bondcombo(atoms, bondcombo_def)
ref_anglecombo = get_anglecombo(atoms, anglecombo_def)
ref_dihedralcombo = get_dihedralcombo(atoms, dihedralcombo_def)
atoms.calc = EMT()
atoms.set_constraint(constr)
opt = BFGS(atoms)
opt.run(fmax=0.01)
new_bondcombo = get_bondcombo(atoms, bondcombo_def)
new_anglecombo = get_anglecombo(atoms, anglecombo_def)
new_dihedralcombo = get_dihedralcombo(atoms, dihedralcombo_def)
err_bondcombo = new_bondcombo - ref_bondcombo
err_anglecombo = new_anglecombo - ref_anglecombo
err_dihedralcombo = new_dihedralcombo - ref_dihedralcombo
print('error in bondcombo:', repr(err_bondcombo))
print('error in anglecombo:', repr(err_anglecombo))
print('error in dihedralcombo:', repr(err_dihedralcombo))
for err in [err_bondcombo, err_anglecombo, err_dihedralcombo]:
assert err < 2e-12
def test_index_shuffle():
(atoms, constr, bond_def, target_bond, angle_def, target_angle,
dihedral_def, target_dihedral) = setup_fixinternals()
constr2 = copy.deepcopy(constr)
# test no change, test constr.get_indices()
assert all(a == b for a, b in zip(constr.get_indices(), (0, 1, 2, 6, 8)))
constr.index_shuffle(atoms, range(len(atoms)))
assert all(a == b for a, b in zip(constr.get_indices(), (0, 1, 2, 6, 8)))
# test full constraint is not part of new slice
with pytest.raises(IndexError):
constr.index_shuffle(atoms, [0])
# test correct shuffling
constr2.index_shuffle(atoms, [1, 2, 0, 6])
assert constr2.bonds[0][1] == [0, 1]
assert constr2.angles[0][1] == [3, 2, 0]
assert constr2.dihedrals[0][1] == [3, 2, 0, 1]
def test_combo_index_shuffle():
(atoms, constr, bondcombo_def, target_bondcombo, anglecombo_def,
target_anglecombo, dihedralcombo_def,
target_dihedralcombo) = setup_combos()
# test no change, test constr.get_indices()
answer = (0, 1, 2, 3, 4, 6, 7, 8)
assert all(a == b for a, b in zip(constr.get_indices(), answer))
constr.index_shuffle(atoms, range(len(atoms)))
assert all(a == b for a, b in zip(constr.get_indices(), answer))
# test anglecombo not part of slice
constr.index_shuffle(atoms, [1, 2, 3, 4, 0, 7])
assert constr.bondcombos[0][1] == [[1, 0, 1.0], [1, 2, -1.0]]
assert constr.dihedralcombos[0][1] == [[2, 1, 0, 3, 1.0], [1, 0, 4, 5, 1.0]]
|
