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 package ffx.potential.extended;
40
41 import edu.rit.pj.reduction.SharedDouble;
42 import ffx.numerics.Constraint;
43 import ffx.numerics.Potential;
44 import ffx.potential.ForceFieldEnergy;
45 import ffx.potential.MolecularAssembly;
46 import ffx.potential.SystemState;
47 import ffx.potential.bonded.AminoAcidUtils;
48 import ffx.potential.bonded.AminoAcidUtils.AminoAcid3;
49 import ffx.potential.bonded.Atom;
50 import ffx.potential.bonded.Residue;
51 import ffx.potential.constraint.ShakeChargeConstraint;
52 import ffx.potential.parameters.ForceField;
53 import ffx.potential.parameters.PolarizeType;
54 import ffx.potential.parameters.TitrationUtils;
55 import ffx.potential.parsers.ESVFilter;
56 import ffx.utilities.Constants;
57 import ffx.utilities.FileUtils;
58 import org.apache.commons.configuration2.CompositeConfiguration;
59 import org.apache.commons.io.FilenameUtils;
60
61 import java.io.File;
62 import java.util.ArrayList;
63 import java.util.Arrays;
64 import java.util.Collections;
65 import java.util.List;
66 import java.util.Random;
67 import java.util.logging.Level;
68 import java.util.logging.Logger;
69
70 import static ffx.potential.bonded.BondedUtils.hasAttachedAtom;
71 import static ffx.utilities.Constants.kB;
72 import static java.lang.String.format;
73 import static org.apache.commons.math3.util.FastMath.log;
74 import static org.apache.commons.math3.util.FastMath.random;
75 import static org.apache.commons.math3.util.FastMath.sin;
76 import static org.apache.commons.math3.util.FastMath.sqrt;
77
78
79
80
81
82
83
84 public class ExtendedSystem implements Potential {
85
86 private static final double DISCR_BIAS = 1.0;
87 private static final double LOG10 = log(10.0);
88 private static final double THETA_FRICTION = 0.5;
89 private static final double THETA_MASS = 5.0;
90 private static final int dDiscr_dTautIndex = 6;
91 private static final int dDiscr_dTitrIndex = 3;
92 private static final int dModel_dTautIndex = 8;
93 private static final int dModel_dTitrIndex = 5;
94 private static final int dPh_dTautIndex = 7;
95 private static final int dPh_dTitrIndex = 4;
96 private static final int discrBiasIndex = 0;
97 private static final Logger logger = Logger.getLogger(ExtendedSystem.class.getName());
98 private static final int modelBiasIndex = 2;
99 private static final int pHBiasIndex = 1;
100 public final boolean guessTitrState;
101
102
103
104
105 public final AminoAcid3[] residueNames;
106
107
108
109
110
111
112 public final int[] tautomerDirections;
113 private final double[] ASH1toASH2 = new double[4];
114
115
116
117
118 private final double[] ASHFmod = new double[4];
119 private final double ASHrestraintConstant;
120
121
122
123 private final double ASHtautBiasMag;
124 private final double ASHtitrBiasMag;
125 private final double[] CYSFmod = new double[4];
126 private final double CYSrestraintConstant;
127
128
129
130 private final double CYStitrBiasMag;
131 private final double[] GLH1toGLH2 = new double[4];
132 private final double[] GLHFmod = new double[4];
133 private final double GLHrestraintConstant;
134
135
136
137 private final double GLHtautBiasMag;
138 private final double GLHtitrBiasMag;
139
140
141
142
143 private final double[] HIDFmod = new double[4];
144 private final double[] HIDtoHIEFmod = new double[4];
145 private final double[] HIEFmod = new double[4];
146 private final double HISrestraintConstant;
147
148
149
150 private final double HIStautBiasMag;
151 private final double HIStitrBiasMag;
152 private final double[] LYSFmod = new double[4];
153 private final double LYSrestraintConstant;
154
155
156
157 private final double LYStitrBiasMag;
158 private final List<Constraint> constraints;
159 public final boolean useTotalChargeCorrection;
160 private final boolean doBias;
161 private final boolean doElectrostatics;
162 private final boolean doPolarization;
163
164 private final boolean doVDW;
165
166
167
168 final private int[][][] esvHistogram;
169 private final SharedDouble[] esvIndElecDerivs;
170 private final SharedDouble[] esvPermElecDerivs;
171
172
173
174
175 private final SystemState esvState;
176
177
178
179
180 private final SharedDouble[] esvVdwDerivs;
181
182
183
184 private final Atom[] extendedAtoms;
185
186
187
188 private final double[] extendedLambdas;
189
190
191
192 private final int[] extendedMolecules;
193
194
195
196 private final List<Residue> extendedResidueList;
197
198
199
200 private final ForceFieldEnergy forceFieldEnergy;
201
202
203
204
205
206 private final boolean[] isTautomerizing;
207
208
209
210
211
212 private final boolean[] isTitrating;
213
214
215
216
217 private final boolean[] isTitratingHeavy;
218
219
220
221
222 private final boolean[] isTitratingHydrogen;
223
224
225
226
227
228
229 private final boolean lockStates;
230
231
232
233 private final int nAtoms;
234
235
236
237 private final int nESVs;
238
239
240
241 private final int nTitr;
242
243
244
245
246 private final ArrayList<Double> specialResiduePKAs;
247 private final ArrayList<Double> specialResidues;
248
249
250
251 private final ArrayList<Double> selectResidues;
252 private final int[] tautomerIndexMap;
253
254
255
256
257 private final double[] tautomerLambdas;
258
259
260
261 private final List<Residue> tautomerizingResidueList;
262
263
264
265 private final double thetaFriction;
266
267
268
269 private final double thetaMass;
270
271
272
273 private final List<Residue> titratingResidueList;
274
275
276
277
278 private final int[] titrationIndexMap;
279
280
281
282
283 private final double[] titrationLambdas;
284
285
286
287
288 private final TitrationUtils titrationUtils;
289 private final boolean useChargeConstraint;
290
291
292
293 File restartFile;
294
295
296
297
298 private boolean fixTautomerState;
299 private boolean fixTitrationState;
300
301
302
303 private ESVFilter esvFilter = null;
304
305
306
307 private double constantSystemPh = 7.4;
308
309
310
311 private double currentTemperature = Constants.ROOM_TEMPERATURE;
312
313
314
315
316
317
318 public ExtendedSystem(MolecularAssembly mola, double pH, final File esvFile) {
319 extendedAtoms = mola.getAtomArray();
320 extendedMolecules = mola.getMoleculeNumbers();
321 setConstantPh(pH);
322 ForceField forceField = mola.getForceField();
323 forceFieldEnergy = mola.getPotentialEnergy();
324 if (forceFieldEnergy == null) {
325 logger.severe("No potential energy found?");
326 }
327 CompositeConfiguration properties = mola.getProperties();
328 titrationUtils = new TitrationUtils(forceField);
329
330 doVDW = properties.getBoolean("esv.vdW", true);
331 doElectrostatics = properties.getBoolean("esv.elec", true);
332 doBias = properties.getBoolean("esv.bias", true);
333 doPolarization = properties.getBoolean("esv.polarization", true);
334 useTotalChargeCorrection = properties.getBoolean("esv.total.charge.correction", false);
335 thetaFriction = properties.getDouble("esv.friction", ExtendedSystem.THETA_FRICTION);
336 thetaMass = properties.getDouble("esv.mass", ExtendedSystem.THETA_MASS);
337
338 lockStates = properties.getBoolean("lock.esv.states", false);
339 double initialTitrationLambda = properties.getDouble("initial.titration.lambda", 0.5);
340 double initialTautomerLambda = properties.getDouble("initial.tautomer.lambda", 0.5);
341 guessTitrState = properties.getBoolean("guess.titration.state", false);
342 specialResidues = getPropertyList(properties, "esv.special.residues");
343 specialResiduePKAs = getPropertyList(properties, "esv.special.residues.pka");
344 selectResidues = getPropertyList(properties, "esv.select.residues");
345 initSpecialResidues(specialResidues, specialResiduePKAs, mola);
346
347 fixTitrationState = properties.getBoolean("fix.titration.lambda", false);
348 fixTautomerState = properties.getBoolean("fix.tautomer.lambda", false);
349 useChargeConstraint = properties.getBoolean("esv.charge.constraint", false);
350 int totalCharge = properties.getInt("esv.charge.constraint.value", 0);
351
352
353 ASHFmod[3] = properties.getDouble("ASH.cubic", TitrationUtils.Titration.ASHtoASP.cubic);
354 ASHFmod[2] = properties.getDouble("ASH.quadratic", TitrationUtils.Titration.ASHtoASP.quadratic);
355 ASHFmod[1] = properties.getDouble("ASH.linear", TitrationUtils.Titration.ASHtoASP.linear);
356 ASHFmod[0] = TitrationUtils.Titration.ASHtoASP.offset;
357 ASH1toASH2[3] = properties.getDouble("ASH1toASH2.cubic", TitrationUtils.Titration.ASH1toASH2.cubic);
358 ASH1toASH2[2] = properties.getDouble("ASH1toASH2.quadratic", TitrationUtils.Titration.ASH1toASH2.quadratic);
359 ASH1toASH2[1] = properties.getDouble("ASH1toASH2.linear", TitrationUtils.Titration.ASH1toASH2.linear);
360 ASH1toASH2[0] = TitrationUtils.Titration.ASH1toASH2.offset;
361 ASHtitrBiasMag = properties.getDouble("ASH.titration.bias.magnitude", DISCR_BIAS);
362 ASHtautBiasMag = properties.getDouble("ASH.tautomer.bias.magnitude", DISCR_BIAS);
363 ASHrestraintConstant = properties.getDouble("ASH.restraint.constant", 0.0);
364
365 GLHFmod[3] = properties.getDouble("GLH.cubic", TitrationUtils.Titration.GLHtoGLU.cubic);
366 GLHFmod[2] = properties.getDouble("GLH.quadratic", TitrationUtils.Titration.GLHtoGLU.quadratic);
367 GLHFmod[1] = properties.getDouble("GLH.linear", TitrationUtils.Titration.GLHtoGLU.linear);
368 GLHFmod[0] = TitrationUtils.Titration.GLHtoGLU.offset;
369 GLH1toGLH2[3] = properties.getDouble("GLH1toGLH2.cubic", TitrationUtils.Titration.GLH1toGLH2.cubic);
370 GLH1toGLH2[2] = properties.getDouble("GLH1toGLH2.quadratic", TitrationUtils.Titration.GLH1toGLH2.quadratic);
371 GLH1toGLH2[1] = properties.getDouble("GLH1toGLH2.linear", TitrationUtils.Titration.GLH1toGLH2.linear);
372 GLH1toGLH2[0] = TitrationUtils.Titration.GLH1toGLH2.offset;
373 GLHtitrBiasMag = properties.getDouble("GLH.titration.bias.magnitude", DISCR_BIAS);
374 GLHtautBiasMag = properties.getDouble("GLH.tautomer.bias.magnitude", DISCR_BIAS);
375 GLHrestraintConstant = properties.getDouble("GLH.restraint.constant", 0.0);
376
377 LYSFmod[3] = properties.getDouble("LYS.cubic", TitrationUtils.Titration.LYStoLYD.cubic);
378 LYSFmod[2] = properties.getDouble("LYS.quadratic", TitrationUtils.Titration.LYStoLYD.quadratic);
379 LYSFmod[1] = properties.getDouble("LYS.linear", TitrationUtils.Titration.LYStoLYD.linear);
380 LYSFmod[0] = TitrationUtils.Titration.LYStoLYD.offset;
381 LYStitrBiasMag = properties.getDouble("LYS.titration.bias.magnitude", DISCR_BIAS);
382 LYSrestraintConstant = properties.getDouble("LYS.restraint.constant", 0.0);
383
384 CYSFmod[3] = properties.getDouble("CYS.cubic", TitrationUtils.Titration.CYStoCYD.cubic);
385 CYSFmod[2] = properties.getDouble("CYS.quadratic", TitrationUtils.Titration.CYStoCYD.quadratic);
386 CYSFmod[1] = properties.getDouble("CYS.linear", TitrationUtils.Titration.CYStoCYD.linear);
387 CYSFmod[0] = TitrationUtils.Titration.CYStoCYD.offset;
388 CYStitrBiasMag = properties.getDouble("CYS.titration.bias.magnitude", DISCR_BIAS);
389 CYSrestraintConstant = properties.getDouble("CYS.restraint.constant", 0.0);
390
391 HIDFmod[3] = properties.getDouble("HID.cubic", TitrationUtils.Titration.HIStoHID.cubic);
392 HIDFmod[2] = properties.getDouble("HID.quadratic", TitrationUtils.Titration.HIStoHID.quadratic);
393 HIDFmod[1] = properties.getDouble("HID.linear", TitrationUtils.Titration.HIStoHID.linear);
394 HIDFmod[0] = TitrationUtils.Titration.HIStoHID.offset;
395 HIEFmod[3] = properties.getDouble("HIE.cubic", TitrationUtils.Titration.HIStoHIE.cubic);
396 HIEFmod[2] = properties.getDouble("HIE.quadratic", TitrationUtils.Titration.HIStoHIE.quadratic);
397 HIEFmod[1] = properties.getDouble("HIE.linear", TitrationUtils.Titration.HIStoHIE.linear);
398 HIEFmod[0] = TitrationUtils.Titration.HIStoHIE.offset;
399 HIDtoHIEFmod[3] = properties.getDouble("HIDtoHIE.cubic", TitrationUtils.Titration.HIDtoHIE.cubic);
400 HIDtoHIEFmod[2] = properties.getDouble("HIDtoHIE.quadratic", TitrationUtils.Titration.HIDtoHIE.quadratic);
401 HIDtoHIEFmod[1] = properties.getDouble("HIDtoHIE.linear", TitrationUtils.Titration.HIDtoHIE.linear);
402 HIDtoHIEFmod[0] = TitrationUtils.Titration.HIDtoHIE.offset;
403 HIStitrBiasMag = properties.getDouble("HIS.titration.bias.magnitude", DISCR_BIAS);
404 HIStautBiasMag = properties.getDouble("HIS.tautomer.bias.magnitude", DISCR_BIAS);
405 HISrestraintConstant = properties.getDouble("HIS.restraint.constant", 0.0);
406
407
408 logger.info("\n Titration bias magnitudes:");
409 logger.info(" Lysine titration bias magnitude: " + LYStitrBiasMag);
410 logger.info(" Cysteine titration bias magnitude: " + CYStitrBiasMag);
411 logger.info(" Histidine titration bias magnitude: " + HIStitrBiasMag);
412 logger.info(" Glutamic acid titration bias magnitude: " + GLHtitrBiasMag);
413 logger.info(" Aspartic acid titration bias magnitude: " + ASHtitrBiasMag);
414
415
416 logger.info("\n Tautomer bias magnitudes:");
417 logger.info(" Histidine tautomer bias magnitude: " + HIStautBiasMag);
418 logger.info(" Glutamic acid tautomer bias magnitude: " + GLHtautBiasMag);
419 logger.info(" Aspartic acid tautomer bias magnitude: " + ASHtautBiasMag);
420
421
422 logger.info("\n Titration model bias (Umod) terms:");
423 logger.info(" Lysine cubic term: " + LYSFmod[3]);
424 logger.info(" Lysine quadratic term: " + LYSFmod[2]);
425 logger.info(" Lysine linear term: " + LYSFmod[1]);
426 logger.info(" Cysteine cubic term: " + CYSFmod[3]);
427 logger.info(" Cysteine quadratic term: " + CYSFmod[2]);
428 logger.info(" Cysteine linear term: " + CYSFmod[1]);
429 logger.info(" Histidine cubic term: " + HIDFmod[3]);
430 logger.info(" HID quadratic term: " + HIDFmod[2]);
431 logger.info(" HID linear term: " + HIDFmod[1]);
432 logger.info(" HIE cubic term: " + HIEFmod[3]);
433 logger.info(" HIE quadratic term: " + HIEFmod[2]);
434 logger.info(" HIE linear term: " + HIEFmod[1]);
435 logger.info(" HID-HIE cubic term: " + HIDtoHIEFmod[3]);
436 logger.info(" HID-HIE quadratic term: " + HIDtoHIEFmod[2]);
437 logger.info(" HID-HIE linear term: " + HIDtoHIEFmod[1]);
438 logger.info(" Glutamic acid cubic term: " + GLHFmod[3]);
439 logger.info(" Glutamic acid quadratic term: " + GLHFmod[2]);
440 logger.info(" Glutamic acid linear term: " + GLHFmod[1]);
441 logger.info(" Aspartic acid cubic term: " + ASHFmod[3]);
442 logger.info(" Aspartic acid quadratic term: " + ASHFmod[2]);
443 logger.info(" Aspartic acid linear term: " + ASHFmod[1] + "\n");
444
445 titratingResidueList = new ArrayList<>();
446 tautomerizingResidueList = new ArrayList<>();
447 extendedResidueList = new ArrayList<>();
448
449 Atom[] atoms = mola.getAtomArray();
450 nAtoms = atoms.length;
451 isTitrating = new boolean[nAtoms];
452 isTitratingHydrogen = new boolean[nAtoms];
453 isTitratingHeavy = new boolean[nAtoms];
454 isTautomerizing = new boolean[nAtoms];
455 titrationLambdas = new double[nAtoms];
456 tautomerLambdas = new double[nAtoms];
457 titrationIndexMap = new int[nAtoms];
458 tautomerIndexMap = new int[nAtoms];
459 tautomerDirections = new int[nAtoms];
460 residueNames = new AminoAcid3[nAtoms];
461
462 Arrays.fill(isTitrating, false);
463 Arrays.fill(isTitratingHydrogen, false);
464 Arrays.fill(isTitratingHeavy, false);
465 Arrays.fill(isTautomerizing, false);
466 Arrays.fill(titrationLambdas, 1.0);
467 Arrays.fill(tautomerLambdas, 0.0);
468 Arrays.fill(titrationIndexMap, -1);
469 Arrays.fill(tautomerIndexMap, -1);
470 Arrays.fill(tautomerDirections, 0);
471 Arrays.fill(residueNames, AminoAcid3.UNK);
472
473
474
475
476
477
478 List<Residue> residueList = mola.getResidueList();
479 List<Residue> preprocessList = new ArrayList<>(residueList);
480 for (Residue residue : preprocessList) {
481 List<Atom> atomList = residue.getSideChainAtoms();
482 for (Atom atom : atomList) {
483
484 if (atom.getAtomicNumber() == 16) {
485 if (hasAttachedAtom(atom, 16)) {
486 residueList.remove(residue);
487 }
488 }
489 }
490 }
491
492 residueList.removeIf(residue -> (residue.getResidueType() == Residue.ResidueType.NA));
493 for (Residue residue : residueList) {
494 if (isTitratable(residue)) {
495 if(selectResidues.isEmpty() || selectResidues.contains((double) residue.getResidueNumber())){
496 titratingResidueList.add(residue);
497 List<Atom> atomList = residue.getSideChainAtoms();
498 for (Atom atom : atomList) {
499 int atomIndex = atom.getArrayIndex();
500 residueNames[atomIndex] = residue.getAminoAcid3();
501 isTitrating[atomIndex] = true;
502 titrationLambdas[atomIndex] = initialTitrationState(residue, initialTitrationLambda, guessTitrState);
503 int titrationIndex = titratingResidueList.indexOf(residue);
504 titrationIndexMap[atomIndex] = titrationIndex;
505 isTitratingHydrogen[atomIndex] = TitrationUtils.isTitratingHydrogen(residue.getAminoAcid3(), atom);
506 isTitratingHeavy[atomIndex] = TitrationUtils.isTitratingHeavy(residue.getAminoAcid3(), atom);
507
508
509
510 if (isTitratingHeavy(atomIndex)) {
511
512 if (atom.getPolarizeType() != null) {
513 double deprotPolar = titrationUtils.getPolarizability(atom, 0.0, 0.0, atom.getPolarizeType().polarizability);
514 double protPolar = titrationUtils.getPolarizability(atom, 1.0, 1.0, atom.getPolarizeType().polarizability);
515 double avgPolar = 0.5 * deprotPolar + 0.5 * protPolar;
516 PolarizeType esvPolarizingHeavyAtom = new PolarizeType(atom.getType(), avgPolar, atom.getPolarizeType().thole,
517 atom.getPolarizeType().ddp, atom.getPolarizeType().polarizationGroup);
518 atom.setPolarizeType(esvPolarizingHeavyAtom);
519 }
520 }
521 }
522
523 if (isTautomer(residue)) {
524 tautomerizingResidueList.add(residue);
525 for (Atom atom : atomList) {
526 int atomIndex = atom.getArrayIndex();
527 isTautomerizing[atomIndex] = true;
528 tautomerLambdas[atomIndex] = initialTautomerLambda;
529 int tautomerIndex = tautomerizingResidueList.indexOf(residue);
530 tautomerIndexMap[atomIndex] = tautomerIndex;
531 tautomerDirections[atomIndex] = TitrationUtils.getTitratingHydrogenDirection(residue.getAminoAcid3(), atom);
532 }
533 }
534
535 assert (titrationUtils.testResidueTypes(residue));
536 }
537 }
538 }
539
540
541 extendedResidueList.addAll(titratingResidueList);
542 extendedResidueList.addAll(tautomerizingResidueList);
543
544 nESVs = extendedResidueList.size();
545 nTitr = titratingResidueList.size();
546 extendedLambdas = new double[nESVs];
547 esvState = new SystemState(nESVs);
548
549
550
551
552 esvHistogram = new int[nTitr][10][10];
553 esvVdwDerivs = new SharedDouble[nESVs];
554 esvPermElecDerivs = new SharedDouble[nESVs];
555 esvIndElecDerivs = new SharedDouble[nESVs];
556 for (int i = 0; i < nESVs; i++) {
557 esvVdwDerivs[i] = new SharedDouble(0.0);
558 esvPermElecDerivs[i] = new SharedDouble(0.0);
559 esvIndElecDerivs[i] = new SharedDouble(0.0);
560 }
561 if (useChargeConstraint) {
562 constraints = new ArrayList<>();
563 ShakeChargeConstraint chargeConstraint = new ShakeChargeConstraint(nTitr, totalCharge, 0.001);
564 constraints.add(chargeConstraint);
565 } else {
566 constraints = Collections.emptyList();
567 }
568
569
570 Arrays.fill(esvState.getMass(), thetaMass);
571
572 for (int i = 0; i < nESVs; i++) {
573 if (i < nTitr) {
574 Residue residue = extendedResidueList.get(i);
575 double initialTitrLambda = initialTitrationState(residue, initialTitrationLambda, guessTitrState);
576 initializeThetaArrays(i, initialTitrLambda);
577 } else {
578 initializeThetaArrays(i, initialTautomerLambda);
579 }
580 }
581 if (esvFilter == null) {
582 esvFilter = new ESVFilter(mola.getName());
583 }
584 if (esvFile == null) {
585 String firstFileName = FilenameUtils.removeExtension(mola.getFile().getAbsolutePath());
586 restartFile = new File(firstFileName + ".esv");
587 } else {
588 double[] thetaPosition = esvState.x();
589 double[] thetaVelocity = esvState.v();
590 double[] thetaAccel = esvState.a();
591 if (!esvFilter.readESV(esvFile, thetaPosition, thetaVelocity, thetaAccel, esvHistogram)) {
592 String message = " Could not load the restart file - dynamics terminated.";
593 logger.log(Level.WARNING, message);
594 throw new IllegalStateException(message);
595 } else {
596 restartFile = esvFile;
597 updateLambdas();
598 }
599 }
600 logger.info("\n Extended System created for the following residues: " + titratingResidueList);
601 }
602
603
604
605
606 private void initSpecialResidues(ArrayList<Double> specialResidues, ArrayList<Double> specialResiduePKAs, MolecularAssembly mola) {
607 if (specialResidues.size() != specialResiduePKAs.size() && !specialResiduePKAs.isEmpty()) {
608 logger.severe("The number of special residues and their associated values do not match.");
609 } else if (!specialResidues.isEmpty()) {
610 logger.info("\n Special residues and their associated values:");
611 for (int i = 0; i < specialResidues.size(); i++) {
612 int resNum = (int) (double) specialResidues.get(i) - mola.getResidueList().get(0).getResidueNumber();
613 if (!specialResiduePKAs.isEmpty()) {
614 logger.info(" Residue: " + specialResidues.get(i) + "-" +
615 mola.getResidueList().get(resNum).getName()
616 + " Pka: " + specialResiduePKAs.get(i));
617 }
618 }
619 logger.info(" ");
620 }
621 logger.info(" Special residues: " + specialResidues);
622 logger.info(" Special residues pKa: " + specialResiduePKAs);
623 for (Residue res : mola.getResidueList()) {
624 if (!isTitratable(res) && specialResidues.contains((double) res.getResidueNumber())) {
625 logger.severe("Given special residue: " + res + " is not titratable.");
626 }
627 }
628 }
629
630
631
632
633 public boolean isTitratable(Residue residue) {
634 if (residue.getResidueType() == Residue.ResidueType.NA) {
635 return false;
636 }
637 AminoAcidUtils.AminoAcid3 AA3 = residue.getAminoAcid3();
638 return AA3.isConstantPhTitratable;
639 }
640
641
642 private ArrayList<Double> getPropertyList(CompositeConfiguration properties, String s) {
643 ArrayList<Double> list = new ArrayList<>();
644 String[] split = properties.getString(s, "").trim()
645 .replace("[", "")
646 .replace("]", "")
647 .replace(",", " ")
648 .split(" ");
649 for (String s1 : split) {
650 if (s1.isEmpty()) {
651 continue;
652 }
653 list.add(Double.parseDouble(s1));
654 }
655 return list;
656 }
657
658
659
660
661
662
663
664
665
666
667 private void initializeThetaArrays(int index, double lambda) {
668 extendedLambdas[index] = lambda;
669 double[] thetaPosition = esvState.x();
670 double[] thetaVelocity = esvState.v();
671 double[] thetaAccel = esvState.a();
672 thetaPosition[index] = Math.asin(Math.sqrt(lambda));
673
674 Random random = new Random();
675 thetaVelocity[index] = random.nextGaussian() * sqrt(kB * 298.15 / thetaMass);
676 double dUdL = getDerivatives()[index];
677 double dUdTheta = dUdL * sin(2 * thetaPosition[index]);
678 thetaAccel[index] = -Constants.KCAL_TO_GRAM_ANG2_PER_PS2 * dUdTheta / thetaMass;
679
680 }
681
682
683
684
685
686
687 public double[] getDerivatives() {
688 double[] esvDeriv = new double[nESVs];
689 double[] biasDerivComponents = new double[9];
690
691 for (Residue residue : titratingResidueList) {
692 int resTitrIndex = titratingResidueList.indexOf(residue);
693
694 if (doBias) {
695 getBiasTerms(residue, biasDerivComponents);
696
697 esvDeriv[resTitrIndex] += biasDerivComponents[dDiscr_dTitrIndex] + biasDerivComponents[dPh_dTitrIndex] + biasDerivComponents[dModel_dTitrIndex];
698
699 if (isTautomer(residue)) {
700 int resTautIndex = tautomerizingResidueList.indexOf(residue) + nTitr;
701 esvDeriv[resTautIndex] += biasDerivComponents[dDiscr_dTautIndex] + biasDerivComponents[dPh_dTautIndex] + biasDerivComponents[dModel_dTautIndex];
702 }
703 }
704
705 if (doVDW) {
706 esvDeriv[resTitrIndex] += getVdwDeriv(resTitrIndex);
707
708 if (isTautomer(residue)) {
709 int resTautIndex = tautomerizingResidueList.indexOf(residue) + nTitr;
710 esvDeriv[resTautIndex] += getVdwDeriv(resTautIndex);
711 }
712 }
713
714 if (doElectrostatics) {
715 esvDeriv[resTitrIndex] += getPermElecDeriv(resTitrIndex);
716
717 if (isTautomer(residue)) {
718 int resTautIndex = tautomerizingResidueList.indexOf(residue) + nTitr;
719 esvDeriv[resTautIndex] += getPermElecDeriv(resTautIndex);
720 }
721 if (doPolarization) {
722 esvDeriv[resTitrIndex] += getIndElecDeriv(resTitrIndex);
723 if (isTautomer(residue)) {
724 int resTautIndex = tautomerizingResidueList.indexOf(residue) + nTitr;
725 esvDeriv[resTautIndex] += getIndElecDeriv(resTautIndex);
726 }
727 }
728 }
729 }
730 return esvDeriv;
731 }
732
733
734
735
736 private void getBiasTerms(Residue residue, double[] biasEnergyAndDerivs) {
737 AminoAcidUtils.AminoAcid3 AA3 = residue.getAminoAcid3();
738 double titrationLambda = getTitrationLambda(residue);
739 double titrationLambdaSquared = titrationLambda * titrationLambda;
740 double titrationLambdaCubed = titrationLambdaSquared * titrationLambda;
741 double discrBias;
742 double pHBias;
743 double modelBias;
744 double dDiscr_dTitr;
745 double dDiscr_dTaut;
746 double dPh_dTitr;
747 double dPh_dTaut;
748 double dMod_dTitr;
749 double dMod_dTaut;
750
751 if (!doBias) {
752 AA3 = AminoAcid3.UNK;
753 }
754 switch (AA3) {
755 case ASD:
756 case ASH:
757 case ASP:
758
759 double tautomerLambda = getTautomerLambda(residue);
760 discrBias = -4.0 * ASHtitrBiasMag * (titrationLambda - 0.5) * (titrationLambda - 0.5);
761 discrBias += -4.0 * ASHtautBiasMag * (tautomerLambda - 0.5) * (tautomerLambda - 0.5);
762 dDiscr_dTitr = -8.0 * ASHtitrBiasMag * (titrationLambda - 0.5);
763 dDiscr_dTaut = -8.0 * ASHtautBiasMag * (tautomerLambda - 0.5);
764
765
766 double pKa1 = TitrationUtils.Titration.ASHtoASP.pKa;
767 double pKa2 = pKa1;
768 pHBias = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
769 * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
770 dPh_dTitr = LOG10 * Constants.R * currentTemperature * -1.0
771 * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
772 dPh_dTaut = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
773 * ((pKa1 - constantSystemPh) - (pKa2 - constantSystemPh));
774
775
776 double restraint = ASHrestraintConstant;
777 double pHSignedSquared = (pKa1 - constantSystemPh) * Math.abs(pKa1 - constantSystemPh);
778
779 pHBias += restraint * (1.0 - titrationLambda) * pHSignedSquared;
780 dPh_dTitr += restraint * -1.0 * pHSignedSquared;
781
782
783 double coeffA0 = ASH1toASH2[3];
784 double coeffA1 = ASH1toASH2[2];
785 double coeffA2 = ASH1toASH2[1];
786 double coeffB0 = ASHFmod[3];
787 double coeffB1 = ASHFmod[2];
788 double coeffB2 = ASHFmod[1];
789 double coeffC0 = ASHFmod[3];
790 double coeffC1 = ASHFmod[2];
791 double coeffC2 = ASHFmod[1];
792 double tautomerLambdaSquared = tautomerLambda * tautomerLambda;
793 double tautomerLambdaCubed = tautomerLambdaSquared * tautomerLambda;
794 double oneMinusTitrationLambda = (1.0 - titrationLambda);
795 double oneMinusTautomerLambda = (1.0 - tautomerLambda);
796 double coeffBSum = coeffB0 + coeffB1 + coeffB2;
797 double coeffCSum = coeffC0 + coeffC1 + coeffC2;
798 modelBias = titrationLambda * (coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
799 tautomerLambda * (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
800 oneMinusTautomerLambda * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
801
802 oneMinusTitrationLambda * (coeffCSum - coeffBSum) * tautomerLambda;
803 dMod_dTitr = (coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
804 tautomerLambda * (3.0 * coeffB0 * titrationLambdaSquared + 2.0 * coeffB1 * titrationLambda + coeffB2) +
805 oneMinusTautomerLambda * (3.0 * coeffC0 * titrationLambdaSquared + 2.0 * coeffC1 * titrationLambda + coeffC2) +
806 -tautomerLambda * (coeffCSum - coeffBSum);
807 dMod_dTaut = titrationLambda * (3.0 * coeffA0 * tautomerLambdaSquared + 2.0 * coeffA1 * tautomerLambda + coeffA2) +
808 (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
809 -1.0 * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
810 oneMinusTautomerLambda * (coeffCSum - coeffBSum);
811 break;
812 case GLD:
813 case GLH:
814 case GLU:
815
816 tautomerLambda = getTautomerLambda(residue);
817 discrBias = -4.0 * GLHtitrBiasMag * (titrationLambda - 0.5) * (titrationLambda - 0.5);
818 discrBias += -4.0 * GLHtautBiasMag * (tautomerLambda - 0.5) * (tautomerLambda - 0.5);
819 dDiscr_dTitr = -8.0 * GLHtitrBiasMag * (titrationLambda - 0.5);
820 dDiscr_dTaut = -8.0 * GLHtautBiasMag * (tautomerLambda - 0.5);
821
822
823 pKa1 = TitrationUtils.Titration.GLHtoGLU.pKa;
824 pKa2 = pKa1;
825 pHBias = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
826 * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
827 dPh_dTitr = LOG10 * Constants.R * currentTemperature * -1.0
828 * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
829 dPh_dTaut = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
830 * ((pKa1 - constantSystemPh) - (pKa2 - constantSystemPh));
831
832
833 restraint = GLHrestraintConstant;
834 pHSignedSquared = (pKa1 - constantSystemPh) * Math.abs(pKa1 - constantSystemPh);
835
836 pHBias += restraint * (1.0 - titrationLambda) * pHSignedSquared;
837 dPh_dTitr += restraint * -1.0 * pHSignedSquared;
838
839
840 coeffA0 = GLH1toGLH2[3];
841 coeffA1 = GLH1toGLH2[2];
842 coeffA2 = GLH1toGLH2[1];
843 coeffB0 = GLHFmod[3];
844 coeffB1 = GLHFmod[2];
845 coeffB2 = GLHFmod[1];
846 coeffC0 = GLHFmod[3];
847 coeffC1 = GLHFmod[2];
848 coeffC2 = GLHFmod[1];
849 tautomerLambdaSquared = tautomerLambda * tautomerLambda;
850 tautomerLambdaCubed = tautomerLambdaSquared * tautomerLambda;
851 oneMinusTitrationLambda = (1.0 - titrationLambda);
852 oneMinusTautomerLambda = (1.0 - tautomerLambda);
853 coeffBSum = coeffB0 + coeffB1 + coeffB2;
854 coeffCSum = coeffC0 + coeffC1 + coeffC2;
855 modelBias = titrationLambda * (coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
856 tautomerLambda * (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
857 oneMinusTautomerLambda * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
858
859 oneMinusTitrationLambda * (coeffCSum - coeffBSum) * tautomerLambda;
860 dMod_dTitr = (coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
861 tautomerLambda * (3.0 * coeffB0 * titrationLambdaSquared + 2.0 * coeffB1 * titrationLambda + coeffB2) +
862 oneMinusTautomerLambda * (3.0 * coeffC0 * titrationLambdaSquared + 2.0 * coeffC1 * titrationLambda + coeffC2) +
863 -tautomerLambda * (coeffCSum - coeffBSum);
864 dMod_dTaut = titrationLambda * (3.0 * coeffA0 * tautomerLambdaSquared + 2.0 * coeffA1 * tautomerLambda + coeffA2) +
865 (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
866 -1.0 * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
867 oneMinusTautomerLambda * (coeffCSum - coeffBSum);
868 break;
869 case HIS:
870 case HID:
871 case HIE:
872
873 tautomerLambda = getTautomerLambda(residue);
874
875 discrBias = -4.0 * HIStitrBiasMag * (titrationLambda - 0.5) * (titrationLambda - 0.5);
876 discrBias += -4.0 * HIStautBiasMag * (tautomerLambda - 0.5) * (tautomerLambda - 0.5);
877 dDiscr_dTitr = -8.0 * HIStitrBiasMag * (titrationLambda - 0.5);
878 dDiscr_dTaut = -8.0 * HIStautBiasMag * (tautomerLambda - 0.5);
879
880
881
882 pKa1 = TitrationUtils.Titration.HIStoHIE.pKa;
883
884 pKa2 = TitrationUtils.Titration.HIStoHID.pKa;
885 pHBias = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
886 * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
887 dPh_dTitr = LOG10 * Constants.R * currentTemperature * -1.0
888 * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
889 dPh_dTaut = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
890 * ((pKa1 - constantSystemPh) - (pKa2 - constantSystemPh));
891
892
893 restraint = HISrestraintConstant;
894 double pH1SignedSquared = (pKa1 - constantSystemPh) * Math.abs(pKa1 - constantSystemPh);
895 double pH2SignedSquared = (pKa2 - constantSystemPh) * Math.abs(pKa2 - constantSystemPh);
896
897 pHBias += restraint * (1.0 - titrationLambda)
898 * (tautomerLambda * pH1SignedSquared + (1.0 - tautomerLambda) * pH2SignedSquared);
899 dPh_dTitr += restraint * -1.0
900 * (tautomerLambda * pH1SignedSquared + (1.0 - tautomerLambda) * pH2SignedSquared);
901 dPh_dTaut += restraint * (1.0 - titrationLambda)
902 * (pH1SignedSquared - pH2SignedSquared);
903
904
905 coeffA0 = HIDtoHIEFmod[3];
906 coeffA1 = HIDtoHIEFmod[2];
907 coeffA2 = HIDtoHIEFmod[1];
908 coeffB0 = HIEFmod[3];
909 coeffB1 = HIEFmod[2];
910 coeffB2 = HIEFmod[1];
911 coeffC0 = HIDFmod[3];
912 coeffC1 = HIDFmod[2];
913 coeffC2 = HIDFmod[1];
914 tautomerLambdaSquared = tautomerLambda * tautomerLambda;
915 tautomerLambdaCubed = tautomerLambdaSquared * tautomerLambda;
916 oneMinusTitrationLambda = (1.0 - titrationLambda);
917 oneMinusTautomerLambda = (1.0 - tautomerLambda);
918 coeffBSum = coeffB0 + coeffB1 + coeffB2;
919 coeffCSum = coeffC0 + coeffC1 + coeffC2;
920 modelBias = oneMinusTitrationLambda * (coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
921 tautomerLambda * (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
922 oneMinusTautomerLambda * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
923
924 titrationLambda * (coeffCSum - coeffBSum) * tautomerLambda;
925 dMod_dTitr = -(coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
926 tautomerLambda * (3.0 * coeffB0 * titrationLambdaSquared + 2.0 * coeffB1 * titrationLambda + coeffB2) +
927 oneMinusTautomerLambda * (3.0 * coeffC0 * titrationLambdaSquared + 2.0 * coeffC1 * titrationLambda + coeffC2) +
928 tautomerLambda * (coeffCSum - coeffBSum);
929 dMod_dTaut = oneMinusTitrationLambda * (3.0 * coeffA0 * tautomerLambdaSquared + 2.0 * coeffA1 * tautomerLambda + coeffA2) +
930 (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
931 -1.0 * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
932 titrationLambda * (coeffCSum - coeffBSum);
933
934 break;
935 case LYS:
936 case LYD:
937
938 discrBias = -4.0 * LYStitrBiasMag * (titrationLambda - 0.5) * (titrationLambda - 0.5);
939 dDiscr_dTitr = -8.0 * LYStitrBiasMag * (titrationLambda - 0.5);
940 dDiscr_dTaut = 0.0;
941
942
943 pKa1 = TitrationUtils.Titration.LYStoLYD.pKa;
944 pHBias = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda) * (pKa1 - constantSystemPh);
945 dPh_dTitr = LOG10 * Constants.R * currentTemperature * -1.0 * (pKa1 - constantSystemPh);
946 dPh_dTaut = 0.0;
947
948
949 restraint = LYSrestraintConstant;
950 pHSignedSquared = (pKa1 - constantSystemPh) * Math.abs(pKa1 - constantSystemPh);
951
952 pHBias += restraint * (1.0 - titrationLambda) * pHSignedSquared;
953 dPh_dTitr += restraint * -1.0 * pHSignedSquared;
954
955
956 double cubic = LYSFmod[3];
957 double quadratic = LYSFmod[2];
958 double linear = LYSFmod[1];
959 modelBias = cubic * titrationLambdaCubed + quadratic * titrationLambdaSquared + linear * titrationLambda;
960 dMod_dTitr = 3 * cubic * titrationLambdaSquared + 2 * quadratic * titrationLambda + linear;
961 dMod_dTaut = 0.0;
962 break;
963 case CYS:
964 case CYD:
965
966 discrBias = -4.0 * CYStitrBiasMag * (titrationLambda - 0.5) * (titrationLambda - 0.5);
967 dDiscr_dTitr = -8.0 * CYStitrBiasMag * (titrationLambda - 0.5);
968 dDiscr_dTaut = 0.0;
969
970
971 pKa1 = TitrationUtils.Titration.CYStoCYD.pKa;
972 pHBias = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda) * (pKa1 - constantSystemPh);
973 dPh_dTitr = LOG10 * Constants.R * currentTemperature * -1.0 * (pKa1 - constantSystemPh);
974 dPh_dTaut = 0.0;
975
976
977 restraint = CYSrestraintConstant;
978 pHSignedSquared = (pKa1 - constantSystemPh) * Math.abs(pKa1 - constantSystemPh);
979
980 pHBias += restraint * (1.0 - titrationLambda) * pHSignedSquared;
981 dPh_dTitr += restraint * -1.0 * pHSignedSquared;
982
983
984
985 cubic = CYSFmod[3];
986 quadratic = CYSFmod[2];
987 linear = CYSFmod[1];
988 modelBias = cubic * titrationLambdaCubed + quadratic * titrationLambdaSquared + linear * titrationLambda;
989 dMod_dTitr = 3 * cubic * titrationLambdaSquared + 2 * quadratic * titrationLambda + linear;
990 dMod_dTaut = 0.0;
991 break;
992 default:
993 discrBias = 0.0;
994 pHBias = 0.0;
995 modelBias = 0.0;
996 dDiscr_dTitr = 0.0;
997 dDiscr_dTaut = 0.0;
998 dPh_dTitr = 0.0;
999 dPh_dTaut = 0.0;
1000 dMod_dTitr = 0.0;
1001 dMod_dTaut = 0.0;
1002 break;
1003 }
1004 biasEnergyAndDerivs[0] = discrBias;
1005 biasEnergyAndDerivs[1] = pHBias;
1006 biasEnergyAndDerivs[2] = -modelBias;
1007 biasEnergyAndDerivs[3] = dDiscr_dTitr;
1008 biasEnergyAndDerivs[4] = dPh_dTitr;
1009 biasEnergyAndDerivs[5] = -dMod_dTitr;
1010 biasEnergyAndDerivs[6] = dDiscr_dTaut;
1011 biasEnergyAndDerivs[7] = dPh_dTaut;
1012 biasEnergyAndDerivs[8] = -dMod_dTaut;
1013 }
1014
1015
1016
1017
1018
1019
1020
1021 public double getTitrationLambda(Residue residue) {
1022 if (titratingResidueList.contains(residue)) {
1023 int resIndex = titratingResidueList.indexOf(residue);
1024 return extendedLambdas[resIndex];
1025 } else {
1026 return 1.0;
1027 }
1028 }
1029
1030
1031
1032
1033
1034
1035
1036 public double getTautomerLambda(Residue residue) {
1037 if (tautomerizingResidueList.contains(residue)) {
1038 int resIndex = tautomerizingResidueList.indexOf(residue);
1039 return extendedLambdas[nTitr + resIndex];
1040 } else {
1041 return 1.0;
1042 }
1043 }
1044
1045
1046
1047
1048 private double getVdwDeriv(int esvID) {
1049 return esvVdwDerivs[esvID].get();
1050 }
1051
1052
1053
1054
1055 private double getPermElecDeriv(int esvID) {
1056 return esvPermElecDerivs[esvID].get();
1057 }
1058
1059
1060
1061
1062 private double getIndElecDeriv(int esvID) {
1063 return esvIndElecDerivs[esvID].get();
1064 }
1065
1066
1067
1068
1069 public boolean isTautomer(Residue residue) {
1070 if (residue.getResidueType() == Residue.ResidueType.NA) {
1071 return false;
1072 }
1073 AminoAcidUtils.AminoAcid3 AA3 = residue.getAminoAcid3();
1074 return AA3.isConstantPhTautomer;
1075 }
1076
1077
1078
1079
1080 private void updateLambdas() {
1081
1082 if (lockStates) {
1083 return;
1084 }
1085 double[] thetaPosition = esvState.x();
1086
1087 for (int i = 0; i < nESVs; i++) {
1088
1089 if ((!fixTitrationState && i < nTitr) || (!fixTautomerState && i >= nTitr)) {
1090 double sinTheta = Math.sin(thetaPosition[i]);
1091 extendedLambdas[i] = sinTheta * sinTheta;
1092 }
1093 }
1094 for (int i = 0; i < nAtoms; i++) {
1095 int mappedTitrationIndex = titrationIndexMap[i];
1096 int mappedTautomerIndex = tautomerIndexMap[i] + nTitr;
1097 if (isTitrating(i) && mappedTitrationIndex != -1) {
1098 titrationLambdas[i] = extendedLambdas[mappedTitrationIndex];
1099 }
1100 if (isTautomerizing(i) && mappedTautomerIndex >= nTitr) {
1101 tautomerLambdas[i] = extendedLambdas[mappedTautomerIndex];
1102 }
1103 }
1104 setESVHistogram();
1105 }
1106
1107
1108
1109
1110 public boolean isTitrating(int atomIndex) {
1111 return isTitrating[atomIndex];
1112 }
1113
1114
1115
1116
1117 public boolean isTautomerizing(int atomIndex) {
1118 return isTautomerizing[atomIndex];
1119 }
1120
1121
1122
1123
1124 private void setESVHistogram() {
1125 for (Residue residue : titratingResidueList) {
1126 int index = titratingResidueList.indexOf(residue);
1127 if (residue.getAminoAcid3().equals(AminoAcid3.LYS)) {
1128 double titrLambda = getTitrationLambda(residue);
1129 esvHistogram(index, titrLambda);
1130 } else {
1131 double titrLambda = getTitrationLambda(residue);
1132 double tautLambda = getTautomerLambda(residue);
1133 esvHistogram(index, titrLambda, tautLambda);
1134 }
1135 }
1136 }
1137
1138
1139
1140
1141
1142
1143
1144 private void esvHistogram(int esv, double lambda) {
1145 int value = (int) (lambda * 10.0);
1146
1147 if (value == 10) {
1148 value = 9;
1149 }
1150 esvHistogram[esv][value][0]++;
1151 }
1152
1153
1154
1155
1156
1157
1158
1159
1160 private void esvHistogram(int esv, double titrLambda, double tautLambda) {
1161 int titrValue = (int) (titrLambda * 10.0);
1162
1163 if (titrValue == 10) {
1164 titrValue = 9;
1165 }
1166 int tautValue = (int) (tautLambda * 10.0);
1167 if (tautValue == 10) {
1168 tautValue = 9;
1169 }
1170 esvHistogram[esv][titrValue][tautValue]++;
1171 }
1172
1173
1174
1175
1176 private double initialTitrationState(Residue residue, double initialLambda, boolean guessTitrState) {
1177 AminoAcid3 AA3 = residue.getAminoAcid3();
1178 double residueNumber = residue.getResidueNumber();
1179 double initialTitrationLambda = 0.0;
1180 if (specialResidues.contains(residueNumber)) {
1181 if (!specialResiduePKAs.isEmpty()) {
1182 initialTitrationLambda =
1183 (constantSystemPh < specialResiduePKAs.get(specialResidues.indexOf(residueNumber))) ? 1.0 : 0.0;
1184 }
1185 } else if (!guessTitrState) {
1186
1187 initialTitrationLambda = initialLambda;
1188 } else {
1189 initialTitrationLambda = switch (AA3) {
1190 case ASD -> (constantSystemPh < TitrationUtils.Titration.ASHtoASP.pKa) ? 1.0 : 0.0;
1191 case GLD -> (constantSystemPh < TitrationUtils.Titration.GLHtoGLU.pKa) ? 1.0 : 0.0;
1192 case HIS -> (constantSystemPh < TitrationUtils.Titration.HIStoHID.pKa) ? 1.0 : 0.0;
1193 case LYS -> (constantSystemPh < TitrationUtils.Titration.LYStoLYD.pKa) ? 1.0 : 0.0;
1194 case CYS -> (constantSystemPh < TitrationUtils.Titration.CYStoCYD.pKa) ? 1.0 : 0.0;
1195 default -> initialLambda;
1196 };
1197 }
1198 return initialTitrationLambda;
1199 }
1200
1201
1202
1203
1204
1205
1206
1207
1208 public boolean isTitratingHeavy(int atomIndex) {
1209 return isTitratingHeavy[atomIndex];
1210 }
1211
1212
1213 public ArrayList<Double> getSpecialResidueList() {
1214 return specialResidues;
1215 }
1216
1217
1218
1219
1220 public void setFixedTautomerState(boolean fixTautomerState) {
1221 this.fixTautomerState = fixTautomerState;
1222 }
1223
1224
1225
1226
1227 public void setFixedTitrationState(boolean fixTitrationState) {
1228 this.fixTitrationState = fixTitrationState;
1229 }
1230
1231
1232
1233
1234 public void reGuessLambdas() {
1235 logger.info(" Reinitializing lambdas to match pH");
1236 for (Residue residue : titratingResidueList) {
1237 double lambda = initialTitrationState(residue, 1.0, true);
1238 setTitrationLambda(residue, lambda);
1239 double tautomerLambda = (int) Math.round(random());
1240 setTautomerLambda(residue, tautomerLambda);
1241 }
1242 }
1243
1244
1245
1246
1247
1248
1249
1250 public void setTautomerLambda(Residue residue, double lambda) {
1251 setTautomerLambda(residue, lambda, true);
1252 }
1253
1254
1255
1256
1257
1258
1259
1260
1261 public void setTautomerLambda(Residue residue, double lambda, boolean changeThetas) {
1262 double[] thetaPosition = esvState.x();
1263 if (tautomerizingResidueList.contains(residue) && !lockStates) {
1264
1265
1266 int index = tautomerizingResidueList.indexOf(residue) + nTitr;
1267 extendedLambdas[index] = lambda;
1268 if (changeThetas) {
1269 thetaPosition[index] = Math.asin(Math.sqrt(lambda));
1270 }
1271 List<Atom> currentAtomList = residue.getSideChainAtoms();
1272 for (Atom atom : currentAtomList) {
1273 int atomIndex = atom.getArrayIndex();
1274 tautomerLambdas[atomIndex] = lambda;
1275 }
1276 }
1277
1278
1279 }
1280
1281
1282
1283
1284
1285
1286
1287 public void setTitrationLambda(Residue residue, double lambda) {
1288 setTitrationLambda(residue, lambda, true);
1289 }
1290
1291
1292
1293
1294
1295
1296
1297
1298 public void setTitrationLambda(Residue residue, double lambda, boolean changeThetas) {
1299 double[] thetaPosition = esvState.x();
1300 if (titratingResidueList.contains(residue) && !lockStates) {
1301 int index = titratingResidueList.indexOf(residue);
1302 extendedLambdas[index] = lambda;
1303 if (changeThetas) {
1304 thetaPosition[index] = Math.asin(Math.sqrt(lambda));
1305 }
1306 List<Atom> currentAtomList = residue.getSideChainAtoms();
1307 for (Atom atom : currentAtomList) {
1308 int atomIndex = atom.getArrayIndex();
1309 titrationLambdas[atomIndex] = lambda;
1310 }
1311 }
1312 }
1313
1314
1315
1316
1317
1318
1319
1320 public int[][] getESVHistogram(int[][] histogram) {
1321 for (int i = 0; i < titratingResidueList.size(); i++) {
1322 int h = 0;
1323 for (int j = 0; j < 10; j++) {
1324 for (int k = 0; k < 10; k++) {
1325 histogram[i][h++] = esvHistogram[i][j][k];
1326 }
1327 }
1328 }
1329 return histogram;
1330 }
1331
1332
1333
1334
1335
1336
1337 public void copyESVHistogramTo(int[][] histogram) {
1338 for (int i = 0; i < titratingResidueList.size(); i++) {
1339 int h = 0;
1340 for (int j = 0; j < 10; j++) {
1341 for (int k = 0; k < 10; k++) {
1342 esvHistogram[i][j][k] = histogram[i][h++];
1343 }
1344 }
1345 }
1346 }
1347
1348
1349
1350
1351 public void initEsvVdw() {
1352 for (int i = 0; i < nESVs; i++) {
1353 esvVdwDerivs[i].set(0.0);
1354 }
1355 }
1356
1357
1358
1359
1360 public void initEsvPermElec() {
1361 for (int i = 0; i < nESVs; i++) {
1362 esvPermElecDerivs[i].set(0.0);
1363 }
1364 }
1365
1366
1367
1368
1369 public void initEsvIndElec() {
1370 for (int i = 0; i < nESVs; i++) {
1371 esvIndElecDerivs[i].set(0.0);
1372 }
1373 }
1374
1375
1376
1377
1378
1379
1380
1381 public double getTitrationLambda(int atomIndex) {
1382 return titrationLambdas[atomIndex];
1383 }
1384
1385
1386
1387
1388
1389
1390
1391 public int getTitrationESVIndex(int i) {
1392 return titrationIndexMap[i];
1393 }
1394
1395
1396
1397
1398
1399
1400
1401 public double getTautomerLambda(int atomIndex) {
1402 return tautomerLambdas[atomIndex];
1403 }
1404
1405 public int getTautomerESVIndex(int i) {
1406 return tautomerIndexMap[i];
1407 }
1408
1409
1410
1411
1412
1413
1414 public List<Residue> getTitratingResidueList() {
1415 return titratingResidueList;
1416 }
1417
1418
1419
1420
1421
1422
1423 public List<Residue> getTautomerizingResidueList() {
1424 return tautomerizingResidueList;
1425 }
1426
1427
1428
1429
1430
1431
1432 public List<Residue> getExtendedResidueList() {
1433 return extendedResidueList;
1434 }
1435
1436 public double getThetaMass() {
1437 return thetaMass;
1438 }
1439
1440 public double getThetaFriction() {
1441 return thetaFriction;
1442 }
1443
1444
1445
1446
1447
1448
1449 public double[] getExtendedLambdas() {
1450 double[] lambdas = new double[nESVs];
1451 System.arraycopy(extendedLambdas, 0, lambdas, 0, lambdas.length);
1452 return lambdas;
1453 }
1454
1455
1456
1457
1458
1459
1460 public String getLambdaList() {
1461 if (nESVs < 1) {
1462 return "";
1463 }
1464 StringBuilder sb = new StringBuilder();
1465 for (int i = 0; i < nESVs; i++) {
1466 if (i == 0) {
1467 sb.append("\n Titration Lambdas: ");
1468 }
1469 if (i > 0) {
1470 sb.append(", ");
1471 }
1472 if (i == nTitr) {
1473 sb.append("\n Tautomer Lambdas: ");
1474 }
1475 sb.append(format("%6.4f", extendedLambdas[i]));
1476 }
1477 return sb.toString();
1478 }
1479
1480
1481
1482
1483
1484
1485 public Atom[] getExtendedAtoms() {
1486 return extendedAtoms;
1487 }
1488
1489
1490
1491
1492
1493
1494 public int[] getExtendedMolecule() {
1495 return extendedMolecules;
1496 }
1497
1498
1499
1500
1501
1502
1503 public double getBiasEnergy() {
1504 double totalBiasEnergy = 0.0;
1505 double[] biasEnergyComponents = new double[9];
1506
1507 for (Residue residue : titratingResidueList) {
1508 getBiasTerms(residue, biasEnergyComponents);
1509 double biasEnergy = biasEnergyComponents[discrBiasIndex] + biasEnergyComponents[pHBiasIndex] + biasEnergyComponents[modelBiasIndex];
1510 totalBiasEnergy += biasEnergy;
1511 }
1512 return totalBiasEnergy;
1513 }
1514
1515
1516
1517
1518
1519
1520 public String getBiasDecomposition() {
1521 double discrBias = 0.0;
1522 double phBias = 0.0;
1523 double modelBias = 0.0;
1524 double[] biasEnergyAndDerivs = new double[9];
1525 for (Residue residue : titratingResidueList) {
1526 getBiasTerms(residue, biasEnergyAndDerivs);
1527 discrBias += biasEnergyAndDerivs[discrBiasIndex];
1528 phBias += biasEnergyAndDerivs[pHBiasIndex];
1529
1530 modelBias += biasEnergyAndDerivs[modelBiasIndex];
1531 }
1532 return format(" %-16s %16.8f\n", "Discretizer", discrBias)
1533 + format(" %-16s %16.8f\n", "Acidostat", phBias)
1534 + format(" %-16s %16.8f\n", "Fmod", modelBias);
1535 }
1536
1537
1538
1539
1540 public void getVdwPrefactor(int atomIndex, double[] vdwPrefactorAndDerivs) {
1541 double prefactor = 1.0;
1542 double titrationDeriv = 0.0;
1543 double tautomerDeriv = 0.0;
1544 if (!isTitratingHydrogen(atomIndex) || !doVDW) {
1545 vdwPrefactorAndDerivs[0] = prefactor;
1546 vdwPrefactorAndDerivs[1] = titrationDeriv;
1547 vdwPrefactorAndDerivs[2] = tautomerDeriv;
1548 return;
1549 }
1550 AminoAcid3 AA3 = residueNames[atomIndex];
1551 switch (AA3) {
1552 case ASD:
1553 case GLD:
1554 if (tautomerDirections[atomIndex] == 1) {
1555 prefactor = titrationLambdas[atomIndex] * tautomerLambdas[atomIndex];
1556 titrationDeriv = tautomerLambdas[atomIndex];
1557 tautomerDeriv = titrationLambdas[atomIndex];
1558 } else if (tautomerDirections[atomIndex] == -1) {
1559 prefactor = titrationLambdas[atomIndex] * (1.0 - tautomerLambdas[atomIndex]);
1560 titrationDeriv = (1.0 - tautomerLambdas[atomIndex]);
1561 tautomerDeriv = -titrationLambdas[atomIndex];
1562 }
1563 break;
1564 case HIS:
1565 if (tautomerDirections[atomIndex] == 1) {
1566 prefactor = (1.0 - titrationLambdas[atomIndex]) * tautomerLambdas[atomIndex] + titrationLambdas[atomIndex];
1567 titrationDeriv = -tautomerLambdas[atomIndex] + 1.0;
1568 tautomerDeriv = (1 - titrationLambdas[atomIndex]);
1569 } else if (tautomerDirections[atomIndex] == -1) {
1570 prefactor = (1.0 - titrationLambdas[atomIndex]) * (1.0 - tautomerLambdas[atomIndex]) + titrationLambdas[atomIndex];
1571 titrationDeriv = tautomerLambdas[atomIndex];
1572 tautomerDeriv = -(1.0 - titrationLambdas[atomIndex]);
1573 }
1574 break;
1575 case LYS:
1576 case CYS:
1577 prefactor = titrationLambdas[atomIndex];
1578 titrationDeriv = 1.0;
1579 tautomerDeriv = 0.0;
1580 break;
1581 }
1582 vdwPrefactorAndDerivs[0] = prefactor;
1583 vdwPrefactorAndDerivs[1] = titrationDeriv;
1584 vdwPrefactorAndDerivs[2] = tautomerDeriv;
1585 }
1586
1587
1588
1589
1590
1591
1592
1593 public boolean isTitratingHydrogen(int atomIndex) {
1594 return isTitratingHydrogen[atomIndex];
1595 }
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605 public void addVdwDeriv(int atomI, double vdwEnergy, double[] vdwPrefactorAndDerivI, double vdwPrefactorJ) {
1606 if (!isTitratingHydrogen(atomI)) {
1607 return;
1608 }
1609
1610
1611
1612 int titrationEsvIndex = titrationIndexMap[atomI];
1613 int tautomerEsvIndex = tautomerIndexMap[atomI] + nTitr;
1614 double dTitr_dLambda;
1615 double dTaut_dLambda;
1616
1617 dTitr_dLambda = vdwPrefactorAndDerivI[1] * vdwPrefactorJ * vdwEnergy;
1618 dTaut_dLambda = vdwPrefactorAndDerivI[2] * vdwPrefactorJ * vdwEnergy;
1619
1620 esvVdwDerivs[titrationEsvIndex].addAndGet(dTitr_dLambda);
1621 if (tautomerEsvIndex >= nTitr) {
1622 esvVdwDerivs[tautomerEsvIndex].addAndGet(dTaut_dLambda);
1623 }
1624 }
1625
1626
1627
1628
1629
1630
1631
1632
1633 public void addPermElecDeriv(int atomI, double titrationEnergy, double tautomerEnergy) {
1634 int titrationEsvIndex = titrationIndexMap[atomI];
1635 int tautomerEsvIndex = tautomerIndexMap[atomI] + nTitr;
1636 esvPermElecDerivs[titrationEsvIndex].addAndGet(titrationEnergy);
1637 if (tautomerEsvIndex >= nTitr) {
1638 esvPermElecDerivs[tautomerEsvIndex].addAndGet(tautomerEnergy);
1639 }
1640 }
1641
1642
1643
1644
1645
1646
1647
1648
1649 public void addIndElecDeriv(int atomI, double titrationEnergy, double tautomerEnergy) {
1650 int titrationEsvIndex = titrationIndexMap[atomI];
1651 int tautomerEsvIndex = tautomerIndexMap[atomI] + nTitr;
1652 esvIndElecDerivs[titrationEsvIndex].addAndGet(titrationEnergy);
1653 if (tautomerEsvIndex >= nTitr) {
1654 esvIndElecDerivs[tautomerEsvIndex].addAndGet(tautomerEnergy);
1655 }
1656 }
1657
1658
1659
1660
1661
1662
1663 public double getConstantPh() {
1664 return constantSystemPh;
1665 }
1666
1667
1668
1669
1670
1671
1672 public void setConstantPh(double pH) {
1673 constantSystemPh = pH;
1674 }
1675
1676
1677
1678
1679 public TitrationUtils getTitrationUtils() {
1680 return titrationUtils;
1681 }
1682
1683
1684
1685
1686
1687
1688
1689 public void setRestartFile(File esvFile) {
1690 restartFile = esvFile;
1691 }
1692
1693
1694
1695
1696
1697
1698
1699 public boolean writeESVInfoTo(File esvFile) {
1700 logger.info("Writing pH Dynamics out to: " + esvFile.getParentFile().getName() + File.separator + esvFile.getName());
1701 double[] thetaPosition = esvState.x();
1702 double[] thetaVelocity = esvState.v();
1703 double[] thetaAccel = esvState.a();
1704 return esvFilter.writeESV(esvFile, thetaPosition, thetaVelocity, thetaAccel, titratingResidueList, esvHistogram, constantSystemPh);
1705 }
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715 public boolean readESVInfoFrom(File esvFile) {
1716 double[] thetaPosition = esvState.x();
1717 double[] thetaVelocity = esvState.v();
1718 double[] thetaAccel = esvState.a();
1719 return esvFilter.readESV(esvFile, thetaPosition, thetaVelocity, thetaAccel, esvHistogram);
1720 }
1721
1722
1723
1724
1725
1726
1727 public void setTemperature(double set) {
1728 currentTemperature = set;
1729 }
1730
1731
1732
1733
1734 public void preForce() {
1735 updateLambdas();
1736 }
1737
1738
1739
1740
1741 public void writeRestart() {
1742 String esvName = FileUtils.relativePathTo(restartFile).toString();
1743 double[] thetaPosition = esvState.x();
1744 double[] thetaVelocity = esvState.v();
1745 double[] thetaAccel = esvState.a();
1746 if (esvFilter.writeESV(restartFile, thetaPosition, thetaVelocity, thetaAccel, titratingResidueList, esvHistogram, constantSystemPh)) {
1747 logger.info(" Wrote PhDynamics restart file to " + esvName);
1748 } else {
1749 logger.info(" Writing PhDynamics restart file to " + esvName + " failed");
1750 }
1751 }
1752
1753
1754
1755
1756
1757 public void writeLambdaHistogram(boolean printHistograms) {
1758 printProtonationRatios();
1759 if (printHistograms) {
1760 logger.info(esvFilter.getLambdaHistogram(titratingResidueList, esvHistogram, constantSystemPh));
1761 }
1762 }
1763
1764
1765
1766
1767 public void printProtonationRatios() {
1768 for (int i = 0; i < esvHistogram.length; i++) {
1769 int[] rowSums = new int[esvHistogram[i].length];
1770 for (int j = 0; j < esvHistogram[i].length; j++) {
1771 for (int k = 0; k < esvHistogram[i][j].length; k++) {
1772 rowSums[j] += esvHistogram[i][j][k];
1773 }
1774 }
1775 int i1 = rowSums[0] + rowSums[rowSums.length - 1];
1776 double buf = i1 == 0 ? 0.0 : .001;
1777 logger.info(" " + extendedResidueList.get(i).toString() + " Deprotonation Fraction at pH " + constantSystemPh + ": " + (rowSums[0] / (i1 + buf)));
1778 if (buf == 0.0) {
1779 logger.info(" Buffer required to avoid division by 0");
1780 }
1781 }
1782 }
1783
1784 @Override
1785 public double[] getAcceleration(double[] acceleration) {
1786 return getThetaAccel();
1787 }
1788
1789 public double[] getThetaAccel() {
1790 return esvState.a();
1791 }
1792
1793 @Override
1794 public double energyAndGradient(double[] x, double[] g) {
1795 double[] thetaPosition = esvState.x();
1796 System.arraycopy(x, 0, thetaPosition, 0, thetaPosition.length);
1797 updateLambdas();
1798 fillESVgradient(g);
1799 return energy(x);
1800 }
1801
1802 @Override
1803 public List<Constraint> getConstraints() {
1804 return constraints.isEmpty() ? Collections.emptyList() : new ArrayList<>(constraints);
1805 }
1806
1807 @Override
1808 public STATE getEnergyTermState() {
1809 return null;
1810 }
1811
1812 private double[] fillESVgradient(double[] g) {
1813 double[] gradESV = postForce();
1814 System.arraycopy(gradESV, 0, g, 0, g.length);
1815 return g;
1816 }
1817
1818 @Override
1819 public void setEnergyTermState(STATE state) {
1820
1821 }
1822
1823 @Override
1824 public double[] getMass() {
1825 return getThetaMassArray();
1826 }
1827
1828
1829
1830
1831
1832
1833 public double[] postForce() {
1834 double[] thetaPosition = esvState.x();
1835 double[] dEdL = ExtendedSystem.this.getDerivatives();
1836 double[] dEdTheta = new double[dEdL.length];
1837 for (int i = 0; i < nESVs; i++) {
1838 dEdTheta[i] = dEdL[i] * sin(2 * thetaPosition[i]);
1839
1840 }
1841 return dEdTheta;
1842 }
1843
1844 public double[] getThetaMassArray() {
1845 return esvState.getMass();
1846 }
1847
1848 @Override
1849 public double[] getPreviousAcceleration(double[] previousAcceleration) {
1850 return new double[0];
1851 }
1852 @Override
1853 public double energy(double[] x) {
1854 double[] coords = new double[forceFieldEnergy.getNumberOfVariables()];
1855 forceFieldEnergy.getCoordinates(coords);
1856 return forceFieldEnergy.energy(coords);
1857 }
1858
1859 @Override
1860 public VARIABLE_TYPE[] getVariableTypes() {
1861 return new VARIABLE_TYPE[0];
1862 }
1863
1864 @Override
1865 public double[] getVelocity(double[] velocity) {
1866 return getThetaVelocity();
1867 }
1868
1869 public double[] getThetaVelocity() {
1870 return esvState.v();
1871 }
1872
1873 @Override
1874 public void setAcceleration(double[] acceleration) {
1875 }
1876
1877 @Override
1878 public double[] getCoordinates(double[] parameters) {
1879 return getThetaPosition();
1880 }
1881
1882 @Override
1883 public void setCoordinates(double[] parameters) {
1884 throw new UnsupportedOperationException("Not supported yet.");
1885 }
1886
1887 @Override
1888 public void setPreviousAcceleration(double[] previousAcceleration) {
1889
1890 }
1891
1892 @Override
1893 public void setVelocity(double[] velocity) {
1894
1895 }
1896
1897 public SystemState getState() {
1898 return esvState;
1899 }
1900
1901 public double[] getThetaPosition() {
1902 return esvState.x();
1903 }
1904
1905 @Override
1906 public int getNumberOfVariables() {
1907 return nESVs;
1908 }
1909
1910 @Override
1911 public double[] getScaling() {
1912 double[] scaling = new double[nESVs];
1913 Arrays.fill(scaling, 1.0);
1914 return scaling;
1915 }
1916
1917 @Override
1918 public void setScaling(double[] scaling) {
1919
1920 }
1921
1922 @Override
1923 public double getTotalEnergy() {
1924 return 0;
1925 }
1926
1927 }