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