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.log(Level.INFO, format(" Angles \t\t%6d\t\t%1d", nAngles, forceGroup));
113 }
114 }
115
116
117
118
119
120
121
122 public static Force constructForce(OpenMMEnergy openMMEnergy) {
123 Angle[] angles = openMMEnergy.getAngles();
124 if (angles == null || angles.length < 1) {
125 return null;
126 }
127 AngleForce angleForce = new AngleForce(openMMEnergy);
128 if (angleForce.nAngles > 0) {
129 return angleForce;
130 }
131 return null;
132 }
133
134
135
136
137
138
139 public void updateForce(OpenMMEnergy openMMEnergy) {
140 Angle[] angles = openMMEnergy.getAngles();
141 if (angles == null || angles.length < 1) {
142 return;
143 }
144
145 DoubleArray parameters = new DoubleArray(0);
146 int index = 0;
147 for (Angle angle : angles) {
148 AngleType.AngleMode angleMode = angle.angleType.angleMode;
149 if (!manyBodyTitration && angleMode == AngleType.AngleMode.IN_PLANE) {
150
151 }
152
153 else if (!rigidHydrogenAngles || !isHydrogenAngle(angle)) {
154 int i1 = angle.getAtom(0).getXyzIndex() - 1;
155 int i2 = angle.getAtom(1).getXyzIndex() - 1;
156 int i3 = angle.getAtom(2).getXyzIndex() - 1;
157 double theta0 = angle.angleType.angle[angle.nh];
158 double k = OpenMM_KJPerKcal * angle.angleType.angleUnit * angle.angleType.forceConstant;
159 if (angleMode == AngleType.AngleMode.IN_PLANE) {
160
161 k = 0.0;
162 }
163 parameters.append(theta0);
164 parameters.append(k);
165 setAngleParameters(index++, i1, i2, i3, parameters);
166 parameters.resize(0);
167 }
168 }
169 parameters.destroy();
170 updateParametersInContext(openMMEnergy.getContext());
171 }
172
173
174
175
176
177
178
179
180 private boolean isHydrogenAngle(Angle angle) {
181 if (angle.containsHydrogen()) {
182
183 double angleVal = angle.angleType.angle[angle.nh];
184
185
186
187 if (angleVal < 160.0) {
188 Atom atom1 = angle.getAtom(0);
189 Atom atom2 = angle.getAtom(1);
190 Atom atom3 = angle.getAtom(2);
191
192
193 return atom1.isHydrogen() && atom3.isHydrogen() && !atom2.isHydrogen();
194 }
195 }
196 return false;
197 }
198
199 }