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 package ffx.potential.openmm;
39
40 import ffx.openmm.DoubleArray;
41 import ffx.openmm.Force;
42 import ffx.openmm.CustomAngleForce;
43 import ffx.potential.bonded.Angle;
44 import ffx.potential.bonded.Atom;
45 import ffx.potential.parameters.AngleType;
46 import ffx.potential.parameters.ForceField;
47
48 import java.util.logging.Level;
49 import java.util.logging.Logger;
50
51 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
52 import static java.lang.String.format;
53
54
55
56
57 public class AngleForce extends CustomAngleForce {
58
59 private static final Logger logger = Logger.getLogger(AngleForce.class.getName());
60
61 private int nAngles = 0;
62 private final boolean manyBodyTitration;
63 private final boolean rigidHydrogenAngles;
64
65
66
67
68
69
70 public AngleForce(OpenMMEnergy openMMEnergy) {
71 super(openMMEnergy.getAngleEnergyString());
72 ForceField forceField = openMMEnergy.getMolecularAssembly().getForceField();
73 manyBodyTitration = forceField.getBoolean("MANYBODY_TITRATION", false);
74 rigidHydrogenAngles = forceField.getBoolean("RIGID_HYDROGEN_ANGLES", false);
75 addPerAngleParameter("theta0");
76 addPerAngleParameter("k");
77 setName("Angle");
78
79 DoubleArray parameters = new DoubleArray(0);
80 Angle[] angles = openMMEnergy.getAngles();
81 for (Angle angle : angles) {
82 AngleType angleType = angle.getAngleType();
83 AngleType.AngleMode angleMode = angleType.angleMode;
84 if (!manyBodyTitration && angleMode == AngleType.AngleMode.IN_PLANE) {
85
86 } else if (isHydrogenAngle(angle) && rigidHydrogenAngles) {
87 logger.log(Level.INFO, " Constrained angle %s was not added the AngleForce.", angle);
88 } else {
89 int i1 = angle.getAtom(0).getXyzIndex() - 1;
90 int i2 = angle.getAtom(1).getXyzIndex() - 1;
91 int i3 = angle.getAtom(2).getXyzIndex() - 1;
92
93 double theta0 = angleType.angle[angle.nh];
94 double k = OpenMM_KJPerKcal * angleType.angleUnit * angleType.forceConstant;
95 if (angleMode == AngleType.AngleMode.IN_PLANE) {
96
97
98 k = 0.0;
99 }
100 parameters.append(theta0);
101 parameters.append(k);
102 addAngle(i1, i2, i3, parameters);
103 nAngles++;
104 parameters.resize(0);
105 }
106 }
107 parameters.destroy();
108
109 if (nAngles > 0) {
110 int forceGroup = forceField.getInteger("ANGLE_FORCE_GROUP", 0);
111 setForceGroup(forceGroup);
112 logger.info(format(" Angles: %10d", nAngles));
113 logger.fine(format(" Force Group: %10d", forceGroup));
114 }
115 }
116
117
118
119
120
121
122
123 public static Force constructForce(OpenMMEnergy openMMEnergy) {
124 Angle[] angles = openMMEnergy.getAngles();
125 if (angles == null || angles.length < 1) {
126 return null;
127 }
128 AngleForce angleForce = new AngleForce(openMMEnergy);
129 if (angleForce.nAngles > 0) {
130 return angleForce;
131 }
132 return null;
133 }
134
135
136
137
138
139
140 public void updateForce(OpenMMEnergy openMMEnergy) {
141 Angle[] angles = openMMEnergy.getAngles();
142 if (angles == null || angles.length < 1) {
143 return;
144 }
145
146 DoubleArray parameters = new DoubleArray(0);
147 int index = 0;
148 for (Angle angle : angles) {
149 AngleType.AngleMode angleMode = angle.angleType.angleMode;
150 if (!manyBodyTitration && angleMode == AngleType.AngleMode.IN_PLANE) {
151
152 }
153
154 else if (!rigidHydrogenAngles || !isHydrogenAngle(angle)) {
155 int i1 = angle.getAtom(0).getXyzIndex() - 1;
156 int i2 = angle.getAtom(1).getXyzIndex() - 1;
157 int i3 = angle.getAtom(2).getXyzIndex() - 1;
158 double theta0 = angle.angleType.angle[angle.nh];
159 double k = OpenMM_KJPerKcal * angle.angleType.angleUnit * angle.angleType.forceConstant;
160 if (angleMode == AngleType.AngleMode.IN_PLANE) {
161
162 k = 0.0;
163 }
164 parameters.append(theta0);
165 parameters.append(k);
166 setAngleParameters(index++, i1, i2, i3, parameters);
167 parameters.resize(0);
168 }
169 }
170 parameters.destroy();
171 updateParametersInContext(openMMEnergy.getContext());
172 }
173
174
175
176
177
178
179
180
181 private boolean isHydrogenAngle(Angle angle) {
182 if (angle.containsHydrogen()) {
183
184 double angleVal = angle.angleType.angle[angle.nh];
185
186
187
188 if (angleVal < 160.0) {
189 Atom atom1 = angle.getAtom(0);
190 Atom atom2 = angle.getAtom(1);
191 Atom atom3 = angle.getAtom(2);
192
193
194 return atom1.isHydrogen() && atom3.isHydrogen() && !atom2.isHydrogen();
195 }
196 }
197 return false;
198 }
199
200 }