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.bonded;
39
40 import static ffx.potential.bonded.BondedUtils.buildBond;
41 import static ffx.potential.bonded.BondedUtils.intxyz;
42 import static ffx.utilities.Constants.kB;
43 import static java.lang.String.format;
44 import static org.apache.commons.math3.util.FastMath.sqrt;
45
46 import ffx.potential.ForceFieldEnergy;
47 import ffx.potential.MolecularAssembly;
48 import ffx.potential.parameters.AngleType;
49 import ffx.potential.parameters.BondType;
50 import ffx.potential.parameters.ForceField;
51 import ffx.potential.parameters.TorsionType;
52 import java.util.ArrayList;
53 import java.util.List;
54 import java.util.concurrent.ThreadLocalRandom;
55 import java.util.logging.Level;
56 import java.util.logging.Logger;
57 import javax.swing.tree.MutableTreeNode;
58
59
60
61
62
63
64
65
66 public class MultiTerminus extends Residue {
67
68 private static final Logger logger = Logger.getLogger(MultiResidue.class.getName());
69 public final END end;
70
71 private final ForceField forceField;
72
73 private final ForceFieldEnergy forceFieldEnergy;
74
75 private final MolecularAssembly molecularAssembly;
76
77 public boolean isCharged;
78
79 private Atom uberH3;
80 private Atom uberHO;
81 private Bond bondH3;
82 private Bond bondHO;
83 private List<ROLS> rolsH3 = new ArrayList<>();
84 private List<ROLS> rolsHO = new ArrayList<>();
85
86
87
88
89
90
91
92
93
94 public MultiTerminus(
95 Residue residue,
96 ForceField forceField,
97 ForceFieldEnergy forceFieldEnergy,
98 MolecularAssembly molecularAssembly) {
99 super(
100 residue.getName(),
101 residue.getResidueNumber(),
102 residue.residueType,
103 residue.getChainID(),
104 residue.getChainID().toString());
105 try {
106 String ffString = forceField.getString("FORCEFIELD");
107 if (ForceField.ForceFieldName.valueOf(ffString)
108 != ForceField.ForceFieldName.AMOEBA_PROTEIN_2013) {
109 logger.severe("MultiTerminus supported only under AMOEBA_PROTEIN_2013.");
110 }
111 } catch (Exception ex) {
112 Logger.getLogger(MultiTerminus.class.getName()).log(Level.SEVERE, null, ex);
113 }
114 this.forceField = forceField;
115 this.forceFieldEnergy = forceFieldEnergy;
116 this.molecularAssembly = molecularAssembly;
117 if (residue.getNextResidue() == null) {
118 end = END.CTERM;
119 isCharged = (residue.getAtomNode("HO") == null);
120 } else if (residue.getPreviousResidue() == null) {
121 end = END.NTERM;
122 isCharged = (residue.getAtomNode("H3") != null);
123 } else {
124 end = null;
125 logger.severe("MultiTerminus constructed for non-terminal residue.");
126 }
127
128 }
129
130
131 @Override
132 public void add(MutableTreeNode mtn) {
133 super.add(mtn);
134 }
135
136
137 @Override
138 public boolean equals(Object object) {
139 if (this == object) {
140 return true;
141 } else if (object == null || getClass() != object.getClass()) {
142 return false;
143 }
144 MultiTerminus other = (MultiTerminus) object;
145 if (this.getResidueNumber() == other.getResidueNumber()
146 && this.isCharged == other.isCharged()) {
147 return true;
148 } else {
149 return false;
150 }
151 }
152
153
154
155
156
157
158 public boolean isCharged() {
159 return isCharged;
160 }
161
162
163
164
165
166 public void titrateTerminus_v0() {
167 logger.info(
168 format(" Titrating residue %s (currently %d).", this.toString(), (isCharged ? 1 : 0)));
169 StringBuilder sb = new StringBuilder();
170 sb.append(" Contents of children: ");
171 for (Atom atom : this.getAtomList()) {
172 sb.append(format("%s, ", atom.getName()));
173 }
174 logger.info(sb.toString());
175
176 Atom N = getBBAtom("N");
177 Atom CA = getBBAtom("CA");
178 Atom C = getBBAtom("C");
179 Atom O = getBBAtom("O");
180 Atom OXT = getBBAtom("OXT");
181 Atom H1 = getBBAtom("H1");
182 Atom H2 = getBBAtom("H2");
183 Atom H3 = getBBAtom("H3");
184 Atom HA = getBBAtom("HA");
185 Atom OH = getBBAtom("OH");
186 Atom HO = getBBAtom("HO");
187 String resName = C.getResidueName();
188 int resSeq = C.getResidueNumber();
189 Character chainID = C.getChainID();
190
191 if (end == END.NTERM) {
192 if (isCharged) {
193 N.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.N.neutType)));
194 H1.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.H1.neutType)));
195 H2.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.H2.neutType)));
196 intxyz(H1, N, 1.02, CA, 109.5, C, 180.0, 0);
197 intxyz(H2, N, 1.02, CA, 109.5, C, 0.0, 0);
198 Bond bondH3 = N.getBond(H3);
199 bondH3.removeFromParent();
200 H3.removeFromParent();
201 H3.setParent(null);
202
203
204
205
206 updateGeometry();
207 logger.info(
208 format(
209 " Finished titration. H3 status: %b %b %b",
210 N.getBond(H3) == null,
211 this.getAtomNode().contains(H3) == null,
212 H3.getParent() == null));
213 } else {
214 if (H3 != null) {
215 logger.warning("N-terminal use state toggled.");
216 }
217 N.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.N.chrgType)));
218 H1.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.H1.chrgType)));
219 H2.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.H2.chrgType)));
220 H3 =
221 new Atom(
222 molecularAssembly.getAtomArray().length,
223 "H3",
224 N.getAltLoc(),
225 new double[3],
226 resName,
227 resSeq,
228 chainID,
229 N.getOccupancy(),
230 N.getTempFactor(),
231 N.getSegID(),
232 true);
233 H3.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.H3.chrgType)));
234 intxyz(H1, N, 1.02, CA, 109.5, C, 180.0, 0);
235 intxyz(H2, N, 1.02, CA, 109.5, C, 60.0, 0);
236 intxyz(H3, N, 1.02, CA, 109.5, C, -60.0, 0);
237 double[] vel = new double[3];
238 N.getVelocity(vel);
239 H3.setVelocity(vel);
240 Bond bondH3 = new Bond(N, H3);
241 bondH3.setBondType(forceField.getBondType(N.getAtomType(), H3.getAtomType()));
242 this.getAtomNode().add(H3);
243 this.getBondList().add(bondH3);
244 this.add(bondH3);
245 updateGeometry();
246 logger.info(
247 format(
248 " Finished titration. H3 statuses: "
249 + "(They have each other: %b %b) (I have them: %b %b) (I have bond: %b)",
250 N.getBond(H3) != null,
251 H3.getBond(N) != null,
252 this.getAtomNode().contains(H3) != null,
253 H3.getParent() == this.getAtomNode(),
254 this.getBondList().contains(bondH3)));
255 logger.info(
256 format(
257 " Bonds from H3: %s %s",
258 H3.getBonds().get(0).get1_2(H3).describe(Atom.Descriptions.XyzIndex_Name),
259 H3.getBonds()
260 .get(0)
261 .get1_2(H3)
262 .getBonds()
263 .get(0)
264 .get1_2(H3.getBonds().get(0).get1_2(H3))
265 .describe(Atom.Descriptions.XyzIndex_Name)));
266 }
267 } else if (end == END.CTERM) {
268 if (isCharged) {
269 if (HO != null) {
270 logger.warning("C-terminal in unusual charge state.");
271 }
272 C.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.C.neutType)));
273 O.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.O.neutType)));
274 OXT.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.OH.neutType)));
275 OXT.setName("OH");
276 OH = OXT;
277 HO =
278 new Atom(
279 molecularAssembly.getAtomArray().length,
280 "HO",
281 OXT.getAltLoc(),
282 new double[3],
283 resName,
284 resSeq,
285 chainID,
286 OXT.getOccupancy(),
287 OXT.getTempFactor(),
288 OXT.getSegID(),
289 true);
290 HO.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.HO.neutType)));
291 intxyz(HO, OXT, 1.02, C, 109.5, CA, -1.7, 0);
292 double vel[] = new double[3];
293 OH.getVelocity(vel);
294 HO.setVelocity(vel);
295 buildBond(OH, HO, forceField, null);
296 this.getAtomNode().add(HO);
297 } else {
298 C.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.C.chrgType)));
299 O.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.O.chrgType)));
300 OH.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.OH.chrgType)));
301 OH.setName("OXT");
302 Bond bondHO = OH.getBond(HO);
303 bondHO.removeFromParent();
304 HO.removeFromParent();
305 HO.setParent(null);
306 this.getBondList().remove(bondHO);
307 this.getAtomNode().remove(HO);
308 updateGeometry();
309 }
310 }
311 isCharged = !isCharged;
312 forceFieldEnergy.reInit();
313 }
314
315
316 private Atom getBBAtom(String name) {
317 List<Atom> list = this.getAtomList();
318 for (Atom atom : list) {
319 if (atom.getName().equals(name)) {
320
321 return atom;
322
323
324
325
326
327
328
329
330
331
332 }
333 }
334 return null;
335 }
336
337 private void updateBondedTerms() {
338 StringBuilder sb = new StringBuilder();
339 sb.append("Updating bonded terms: \n");
340 for (Bond bond : getBondList()) {
341 BondType oldType = bond.bondType;
342 BondType newType = forceField.getBondType(bond.atoms[0].getAtomType(), bond.atoms[1].getAtomType());
343 if (oldType != newType) {
344 sb.append(format(" Bond: %s --> %s \n", bond.bondType, newType));
345 bond.setBondType(newType);
346 if (newType.distance < 0.9 * oldType.distance
347 || newType.distance > 1.1 * oldType.distance) {
348 logger.info(
349 format(
350 " Large bond distance change: %s %s, %.2f --> %.2f ",
351 bond.atoms[0].describe(Atom.Descriptions.XyzIndex_Name),
352 bond.atoms[1].describe(Atom.Descriptions.XyzIndex_Name),
353 oldType.distance,
354 newType.distance));
355 }
356 }
357 }
358 for (Angle angle : getAngleList()) {
359 AngleType oldType = angle.angleType;
360 Angle dummy = Angle.angleFactory(angle.bonds[0], angle.bonds[1], forceField);
361 AngleType newType = dummy.angleType;
362 if (oldType != newType) {
363 sb.append(format(" Angle: %s --> %s \n", angle.angleType, dummy.angleType));
364 angle.setAngleType(dummy.angleType);
365 if (newType.angle[0] < 0.9 * oldType.angle[0]
366 || newType.angle[0] > 1.1 * oldType.angle[0]) {
367 logger.info(
368 format(
369 " Large angle change: %s %s %s, %.2f --> %.2f ",
370 angle.atoms[0].describe(Atom.Descriptions.XyzIndex_Name),
371 angle.atoms[1].describe(Atom.Descriptions.XyzIndex_Name),
372 angle.atoms[2].describe(Atom.Descriptions.XyzIndex_Name),
373 oldType.angle[0],
374 newType.angle[0]));
375 }
376 }
377 }
378 for (Torsion tors : getTorsionList()) {
379 TorsionType oldType = tors.torsionType;
380 Torsion dummy =
381 Torsion.torsionFactory(tors.bonds[0], tors.bonds[1], tors.bonds[2], forceField);
382 TorsionType newType = dummy.torsionType;
383 if (oldType != newType) {
384 sb.append(format(" Torsion: %s --> %s \n", tors.torsionType, dummy.torsionType));
385 tors.torsionType = dummy.torsionType;
386 }
387 }
388 }
389
390
391 private void titrateTerminusByRebuilding() {
392 if (true) {
393 throw new UnsupportedOperationException();
394 }
395
396 Atom CA = (Atom) this.getAtomNode("CA");
397 Atom C = (Atom) this.getAtomNode("C");
398 Atom HA = (Atom) this.getAtomNode("HA");
399 Atom N = (Atom) this.getAtomNode("N");
400 Atom O = (Atom) this.getAtomNode("O");
401 Atom OXT = (Atom) this.getAtomNode("OXT");
402 Atom NH2 = (Atom) this.getAtomNode("NH2");
403
404 if (getNextResidue() == null) {
405
406 } else if (getPreviousResidue() == null) {
407 String resName = C.getResidueName();
408 int resSeq = C.getResidueNumber();
409 Character chainID = C.getChainID();
410 double Cxyz[] = new double[3];
411 double Oxyz[] = new double[3];
412 double OXTxyz[] = new double[3];
413 C.getXYZ(Cxyz);
414 O.getXYZ(Oxyz);
415 OXT.getXYZ(OXTxyz);
416
417 int protCkey = 235;
418 int protOkey = 236;
419 int protOHkey = 237;
420 int protHOkey = 238;
421 Atom protC =
422 new Atom(
423 0,
424 "C",
425 C.getAltLoc(),
426 Cxyz,
427 resName,
428 resSeq,
429 chainID,
430 C.getOccupancy(),
431 C.getTempFactor(),
432 C.getSegID(),
433 true);
434 Atom protO =
435 new Atom(
436 0,
437 "O",
438 O.getAltLoc(),
439 Oxyz,
440 resName,
441 resSeq,
442 chainID,
443 O.getOccupancy(),
444 O.getTempFactor(),
445 O.getSegID(),
446 true);
447 Atom protOH =
448 new Atom(
449 0,
450 "OH",
451 OXT.getAltLoc(),
452 OXTxyz,
453 resName,
454 resSeq,
455 chainID,
456 OXT.getOccupancy(),
457 OXT.getTempFactor(),
458 OXT.getSegID(),
459 true);
460 Atom protHO =
461 new Atom(
462 0,
463 "HO",
464 OXT.getAltLoc(),
465 new double[3],
466 resName,
467 resSeq,
468 chainID,
469 OXT.getOccupancy(),
470 OXT.getTempFactor(),
471 OXT.getSegID(),
472 true);
473 intxyz(protHO, protOH, 1.02, protC, 109.5, CA, -1.7, 0);
474 protC.setAtomType(forceField.getAtomType(Integer.toString(protCkey)));
475 protO.setAtomType(forceField.getAtomType(Integer.toString(protOkey)));
476 protOH.setAtomType(forceField.getAtomType(Integer.toString(protOHkey)));
477 protHO.setAtomType(forceField.getAtomType(Integer.toString(protHOkey)));
478
479 buildBond(CA, protO, forceField, null);
480 buildBond(protC, protO, forceField, null);
481 buildBond(protC, protOH, forceField, null);
482 buildBond(protOH, protHO, forceField, null);
483 }
484 }
485
486
487 private void updateGeometry() {
488
489 List<Atom> atoms = this.getAtomList();
490 List<Bond> bonds = this.getBondList();
491 List<Angle> angles = this.getAngleList();
492 List<Torsion> torsions = this.getTorsionList();
493
494 for (Atom atom : atoms) {
495 atom.clearGeometry();
496 }
497 for (Atom atom : atoms) {
498 for (Bond b : bonds) {
499 if (b.containsAtom(atom)) {
500 atom.setBond(b);
501 }
502 }
503 }
504 for (Atom atom : atoms) {
505 for (Angle a : angles) {
506 if (a.containsAtom(atom)) {
507 atom.setAngle(a);
508 }
509 }
510 }
511 for (Atom atom : atoms) {
512 for (Torsion t : torsions) {
513 if (t.containsAtom(atom)) {
514 atom.setTorsion(t);
515 }
516 }
517 }
518 }
519
520 private void maxwellMe(Atom atom, double temperature) {
521 double[] vv = new double[3];
522 for (int i = 0; i < 3; i++) {
523 vv[i] = ThreadLocalRandom.current().nextGaussian() * sqrt(kB * temperature / atom.getMass());
524 }
525 atom.setVelocity(vv);
526 }
527
528 public enum END {
529 NTERM,
530 CTERM;
531 }
532
533 public enum BB_TYPE {
534
535
536
537
538
539
540
541
542 N(231, 41, 225, 1),
543 C(233, 30, 235, 32),
544 O(234, 31, 236, 33),
545 OXT(234, 31, 236, 33),
546 OH(234, 31, 237, 34),
547 HO(-1, -1, 238, 35),
548 H1(232, 41, 226, 4),
549 H2(232, 41, 226, 4),
550 H3(232, 41, 226, 4);
551 public int chrgType, chrgClass;
552 public int neutType, neutClass;
553
554 BB_TYPE(int chrgType, int chrgClass, int neutType, int neutClass) {
555 this.chrgType = chrgType;
556 this.chrgClass = chrgClass;
557 this.neutType = neutType;
558 this.neutClass = neutClass;
559 }
560 }
561 }