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