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