1 package ffx.algorithms.optimize;
2
3 import com.google.common.collect.MinMaxPriorityQueue;
4 import ffx.potential.AssemblyState;
5 import ffx.potential.ForceFieldEnergy;
6 import ffx.potential.MolecularAssembly;
7 import ffx.potential.bonded.Atom;
8 import ffx.potential.bonded.Bond;
9 import ffx.potential.bonded.Molecule;
10 import ffx.potential.bonded.RestrainDistance;
11 import ffx.potential.parameters.BondType;
12 import ffx.potential.parsers.XYZFilter;
13
14 import java.io.File;
15 import java.util.AbstractQueue;
16 import java.util.ArrayList;
17 import java.util.concurrent.ArrayBlockingQueue;
18 import java.util.logging.Logger;
19
20 import static ffx.potential.utils.Superpose.applyRotation;
21 import static java.lang.String.format;
22 import static org.apache.commons.math3.util.FastMath.abs;
23 import static org.apache.commons.math3.util.FastMath.acos;
24 import static org.apache.commons.math3.util.FastMath.cos;
25 import static org.apache.commons.math3.util.FastMath.min;
26 import static org.apache.commons.math3.util.FastMath.sin;
27 import static org.apache.commons.math3.util.FastMath.sqrt;
28
29
30
31
32
33
34
35
36
37
38
39
40
41 public class ConformationScan {
42
43 private static final Logger logger = Logger.getLogger(ConformationScan.class.getName());
44
45 private final MolecularAssembly mola;
46 private final ForceFieldEnergy forceFieldEnergy;
47 private final AssemblyState initState;
48 private double[] x;
49
50 private final Molecule[] s1;
51 private final Molecule[] s2;
52 private final Atom[] s1Atoms;
53 private final Atom[] s2Atoms;
54 private final ArrayList<Atom> s1TargetAtoms;
55 private final ArrayList<Atom> s2TargetAtoms;
56
57 private AbstractQueue<StateContainer> statesQueue;
58
59 private final double eps;
60 private final int maxIter;
61 private double hBondDist;
62 private double flatBottomRadius;
63 private final boolean tScan;
64 private final boolean excludeH;
65 private final boolean minimize;
66
67 private double m1MinEnergy;
68 private double m2MinEnergy;
69 private double totalMonomerMinimizedEnergy;
70
71 private double minimumEnergy;
72 private double averageEnergy;
73 private double averageEnergyNoOutlier;
74 private double stdOfEnergies;
75 private double stdOfEnergiesNoOutlier;
76
77 public ConformationScan(
78 MolecularAssembly molecularAssembly,
79 Molecule[] systemOne,
80 Molecule[] systemTwo,
81 double minimizeEps,
82 int minimizeMaxIterations,
83 boolean staticTorsionScan,
84 boolean excludeExtraHydrogen,
85 boolean minimize
86 )
87 {
88 mola = molecularAssembly;
89 forceFieldEnergy = molecularAssembly.getPotentialEnergy();
90 initState = new AssemblyState(mola);
91 s1 = systemOne;
92 s2 = systemTwo;
93 eps = minimizeEps;
94 maxIter = minimizeMaxIterations;
95 tScan = staticTorsionScan;
96 excludeH = excludeExtraHydrogen;
97 this.minimize = minimize;
98 s1Atoms = getAtomsFromMoleculeArray(s1);
99 s2Atoms = getAtomsFromMoleculeArray(s2);
100 s1TargetAtoms = new ArrayList<>();
101 s2TargetAtoms = new ArrayList<>();
102 setTargetAtoms(mola.getAtomArray());
103 statesQueue = MinMaxPriorityQueue.maximumSize(s1TargetAtoms.size() * s2TargetAtoms.size())
104 .create();
105 if(!minimize) {
106 statesQueue = new ArrayBlockingQueue<>(s1TargetAtoms.size() * s2TargetAtoms.size() + 1);
107 }
108 }
109
110 public void scan(){
111 minimizeEachMolecule(minimize);
112 double[] zAxis = new double[]{0,0,1};
113
114 int loopCounter = 1;
115 for(Atom a: s1TargetAtoms){
116 for(Atom b: s2TargetAtoms){
117 zAxis[2] = 1;
118 initState.revertState();
119 alignSystemCOMtoAtomVecWithAxis(a, zAxis, s1Atoms);
120 zAxis[2] = -1;
121 logger.info("\n ----- Trial " + loopCounter + " out of " +
122 (s2TargetAtoms.size() * s1TargetAtoms.size()) + " -----");
123 loopCounter++;
124 alignSystemCOMtoAtomVecWithAxis(b, zAxis, s2Atoms);
125
126 hBondDist = 2.0;
127 double[] hBondVector = new double[]{0, 0, a.getZ() - b.getZ() + hBondDist};
128 logger.info(" Initial H-bond distance: " + hBondVector[2]);
129
130 hBondVector[2] += minimizeVector(hBondVector, -15, 15)[2];
131 logger.info(" Best H-bond distance: " + hBondVector[2]);
132 hBondDist = hBondVector[2] - (a.getZ() - b.getZ());
133 flatBottomRadius = abs(hBondDist / 2.0);
134 logger.info(" Flat bottom radius: " + flatBottomRadius);
135
136
137 try {
138 mola.update();
139 forceFieldEnergy.getCoordinates(x);
140 int status = 0;
141 if(minimize) {
142 status = minimizeSystem(a, b);
143 }
144 if(status != -1) {
145 forceFieldEnergy.getCoordinates(x);
146 double e = forceFieldEnergy.energy(x, true) - totalMonomerMinimizedEnergy;
147 logger.info("\n Binding energy of trial " + (loopCounter - 1) + ": " + e);
148 statesQueue.add(new StateContainer(new AssemblyState(mola), e));
149 } else {
150 logger.warning(" Minimization failed. No state will be saved.");
151
152 }
153 } catch (Exception ignored) {
154 logger.warning(" Minimization failed. No state will be saved.");
155
156
157 }
158 }
159 }
160 logger.info("\n ------------------------- End of Trials -------------------------");
161 if(statesQueue.peek() != null) {
162 minimumEnergy = statesQueue.peek().getEnergy();
163 } else{
164 logger.warning(" No states were saved. No minimum energy found.");
165 minimumEnergy = Double.NaN;
166 return;
167 }
168 calculateMeanStd();
169 }
170
171 public void logAllEnergyInformation(){
172 logger.info(format(" Minimum energy of monomer 1: %12.5f", m1MinEnergy));
173 logger.info(format(" Minimum energy of monomer 2: %12.5f", m2MinEnergy));
174 logger.info(format(" Sum of minimum monomers: %12.5f", totalMonomerMinimizedEnergy));
175 logger.info(format(" Minimum binding energy of system: %12.5f", minimumEnergy));
176 logger.info(format(" Average binding energy: %12.5f +- %6.3f", averageEnergy, stdOfEnergies));
177 logger.info(format(" Average binding energy (no outlier): %12.5f +- %6.3f", averageEnergyNoOutlier, stdOfEnergiesNoOutlier));
178 }
179 public boolean writeStructuresToXYZ(File outputFile){
180 XYZFilter xyzFilter = new XYZFilter(outputFile, mola, mola.getForceField(), mola.getForceField().getProperties());
181 logger.info("\n Writing structures to " + outputFile.getAbsolutePath());
182 int count = 1;
183 StateContainer[] states = new StateContainer[statesQueue.size()];
184 if(minimize){
185
186 int index = 0;
187 while(!statesQueue.isEmpty() && index < states.length){
188 states[index] = statesQueue.poll();
189 index++;
190 }
191 if(index < states.length){
192 logger.warning(" Not all states will be saved. ");
193 return false;
194 }
195 }
196 else{
197 statesQueue.toArray(states);
198 }
199 if(states.length == 0){
200 logger.warning(" No states were saved. No structures will be written.");
201 return false;
202 }
203 for(StateContainer s: states){
204 AssemblyState state = s.getState();
205 double energy = s.getEnergy();
206 logger.info(" Writing configuration #" + count + " with energy: " + energy);
207 state.revertState();
208 xyzFilter.writeFile(outputFile, count != 1);
209 count++;
210 }
211 return true;
212 }
213
214 public static ArrayList<Double> logBindingEnergyCalculation(ConformationScan m1, ConformationScan m2, ConformationScan dimer){
215 ArrayList<Double> bindingEnergies = new ArrayList<>();
216 logger.info("\n ------------------------- Binding Energy -------------------------");
217 logger.info(" ddE (of binding) = dE (of coformer dimer binding) - dE (average of both monomer dimer bindings summed)");
218 logger.info("\n Calculation using minimized energies: ");
219 double averageMonomerEnergy = (m1.getMinimumEnergy() + m2.getMinimumEnergy()) / 2;
220 double bindingE = dimer.getMinimumEnergy() - averageMonomerEnergy;
221 logger.info(format(" Minimal Structure ddE = %8.3f - %8.3f = %12.5f",
222 dimer.getMinimumEnergy(),
223 averageMonomerEnergy,
224 bindingE));
225 bindingEnergies.add(bindingE);
226 logger.info("\n Calculation using average energies: ");
227 averageMonomerEnergy = (m1.getAverageEnergy() + m2.getAverageEnergy()) / 2;
228 bindingE = dimer.getAverageEnergy() - averageMonomerEnergy;
229 logger.info(format(" Average ddE = %8.3f - %8.3f = %12.5f",
230 dimer.getAverageEnergy(),
231 averageMonomerEnergy,
232 bindingE));
233 bindingEnergies.add(bindingE);
234 logger.info("\n Calculation using average energies (no outliers): ");
235 averageMonomerEnergy = (m1.getAverageEnergyNoOutlier() + m2.getAverageEnergyNoOutlier()) / 2;
236 bindingE = dimer.getAverageEnergyNoOutlier() - averageMonomerEnergy;
237 logger.info(format(" Average (no outlier) ddE = %8.3f - %8.3f = %12.5f",
238 dimer.getAverageEnergyNoOutlier(),
239 averageMonomerEnergy,
240 bindingE));
241 bindingEnergies.add(bindingE);
242 logger.info("\n N.B. -- Monomer energies are calculated using the minimized structures in all cases.");
243 logger.info(" ------------------------- End of Binding Energy -------------------------\n");
244 return bindingEnergies;
245 }
246
247 public ArrayList<AssemblyState> getStates(){
248 ArrayList<AssemblyState> states = new ArrayList<>();
249 for(StateContainer s: statesQueue){
250 states.add(s.getState());
251 }
252 return states;
253 }
254 public ArrayList<AssemblyState> getStatesWithinEnergy(double energy){
255 ArrayList<AssemblyState> states = new ArrayList<>();
256 StateContainer minState = statesQueue.peek();
257 if(minState != null){
258 double minEnergy = minState.getEnergy();
259 for(StateContainer s: statesQueue){
260 if(s.getEnergy() - minEnergy < energy){
261 states.add(s.getState());
262 }
263 }
264 }
265 return states;
266 }
267 public ArrayList<AssemblyState> getStatesFilteredByRMSD(double rmsdCutoff){return null;}
268 public ArrayList<AssemblyState> getStatesAroundAverage(double averageE, double radialCutoff){return null;}
269 public ArrayList<Double> getEnergies(){
270 ArrayList<Double> energies = new ArrayList<>();
271 for(StateContainer s: statesQueue){
272 energies.add(s.getEnergy());
273 }
274 return energies;
275 }
276 public ArrayList<Double> getEnergiesWithinEnergy(double energy){
277 ArrayList<Double> energies = new ArrayList<>();
278 StateContainer minState = statesQueue.peek();
279 if(minState != null){
280 double minEnergy = minState.getEnergy();
281 for(StateContainer s: statesQueue){
282 if(s.getEnergy() - minEnergy < energy){
283 energies.add(s.getEnergy());
284 }
285 }
286 }
287 return energies;
288 }
289 public double getM1MinEnergy() {return m1MinEnergy;}
290 public double getM2MinEnergy() {return m2MinEnergy;}
291 public double getTotalMonomerMinimizedEnergy() {return totalMonomerMinimizedEnergy;}
292 public double getMinimumEnergy() {return minimumEnergy;}
293 public double getAverageEnergy() {return averageEnergy;}
294 public double getAverageEnergyNoOutlier() {return averageEnergyNoOutlier;}
295 public double getStdOfEnergies() {return stdOfEnergies;}
296 public double getStdOfEnergiesNoOutlier() {return stdOfEnergiesNoOutlier;}
297 private double[] minimizeVector(double[] hBondVector, int lowBound, int highBound) {
298 highBound += 1;
299
300 double[] coarsePotentialSurface = new double[abs(highBound - lowBound)];
301 double[] zSearched = new double[abs(highBound - lowBound)];
302 double[] coarseVector = new double[3];
303 int minIndex = -1;
304 double minE = Double.MAX_VALUE;
305 coarseVector[2] = hBondVector[2] + lowBound;
306 for(int i = 0; i < abs(highBound - lowBound); i++){
307 zSearched[i] = i == 0 ? lowBound : zSearched[i - 1] + 1;
308
309 for(Atom a: s2Atoms){
310 a.move(coarseVector);
311 }
312
313 forceFieldEnergy.getCoordinates(x);
314 coarsePotentialSurface[i] = forceFieldEnergy.energy(x, false);
315 if(coarsePotentialSurface[i] < minE){
316 minE = coarsePotentialSurface[i];
317 minIndex = i;
318 }
319 coarseVector[2] = 1;
320 }
321
322
323 for(int i = 0; i < s2Atoms.length; i++){
324 s2Atoms[i].move(new double[]{0, 0, hBondVector[2] - zSearched[zSearched.length - 1]});
325 }
326
327
328
329
330 double[] refinedVector = new double[3];
331 double a;
332 double aPotential;
333 if (minIndex == 0) {
334 a = zSearched[minIndex + 1];
335 aPotential = coarsePotentialSurface[minIndex + 1];
336 } else if (minIndex == coarsePotentialSurface.length - 1) {
337 a = zSearched[minIndex - 1];
338 aPotential = coarsePotentialSurface[minIndex - 1];
339 } else {
340 a = coarsePotentialSurface[minIndex - 1] < coarsePotentialSurface[minIndex + 1] ?
341 zSearched[minIndex - 1] : zSearched[minIndex + 1];
342 aPotential = min(coarsePotentialSurface[minIndex - 1], coarsePotentialSurface[minIndex + 1]);
343 }
344 double b = zSearched[minIndex];
345 double bPotential = coarsePotentialSurface[minIndex];
346 double c = 0;
347 double convergence = 1e-5;
348 while(abs(aPotential - bPotential) > convergence){
349 refinedVector[2] = (a + b) / 2 - c;
350
351 for(Atom a1: s2Atoms){
352 a1.move(refinedVector);
353 }
354
355 forceFieldEnergy.getCoordinates(x);
356 double e = forceFieldEnergy.energy(x, false);
357
358 minE = min(e, minE);
359 if(aPotential > bPotential && e < aPotential){
360 a = (a+b)/2;
361 aPotential = e;
362 c = a;
363 } else if(e < bPotential){
364 b = (a+b)/2;
365 bPotential = e;
366 c = b;
367 } else{
368
369 c = (a + b) / 2;
370 if(aPotential < bPotential){
371 refinedVector[2] = a - c;
372 for(Atom a1: s2Atoms){
373 a1.move(refinedVector);
374 }
375 } else{
376 refinedVector[2] = b - c;
377 for(Atom a1: s2Atoms){
378 a1.move(refinedVector);
379 }
380 }
381 break;
382 }
383 }
384 refinedVector[2] = aPotential < bPotential ? a : b;
385
386
387 return refinedVector;
388 }
389 private int minimizeEachMolecule(boolean minimize){
390
391 int statOne = 0;
392 Minimize monomerMinEngine;
393 for(Atom a: s2Atoms){ a.setUse(false); }
394 if(minimize) {
395 logger.info("\n --------- Minimize System 1 --------- ");
396 if(tScan) {
397 logger.info("\n --------- System 1 Static Torsion Scan --------- ");
398 for(Molecule m: s1) {
399
400 TorsionSearch m1TorsionSearch = new TorsionSearch(mola, m, 32, 1);
401 m1TorsionSearch.staticAnalysis(0, 100);
402 if (!m1TorsionSearch.getStates().isEmpty()) {
403 AssemblyState minState = m1TorsionSearch.getStates().get(0);
404 minState.revertState();
405 }
406 }
407 }
408 monomerMinEngine = new Minimize(mola, forceFieldEnergy, null);
409 monomerMinEngine.minimize(eps, maxIter).getCoordinates(x);
410 statOne = monomerMinEngine.getStatus();
411 }
412 logger.info("\n --------- System 1 Energy Breakdown --------- ");
413 double monomerEnergy = forceFieldEnergy.energy(x, true);
414 for(Atom a: s2Atoms){ a.setUse(true); }
415
416
417 int statTwo = 0;
418 for(Atom a: s1Atoms){ a.setUse(false); }
419 if(minimize) {
420 logger.info("\n --------- Minimize System 2 --------- ");
421 if(tScan) {
422 logger.info("\n --------- System 2 Static Torsion Scan --------- ");
423 for(Molecule m: s2) {
424
425 TorsionSearch m2TorsionSearch = new TorsionSearch(mola, m, 32, 1);
426 m2TorsionSearch.staticAnalysis(0, 100);
427 if (!m2TorsionSearch.getStates().isEmpty()) {
428 AssemblyState minState = m2TorsionSearch.getStates().get(0);
429 minState.revertState();
430 }
431 }
432 }
433 monomerMinEngine = new Minimize(mola, forceFieldEnergy, null);
434 monomerMinEngine.minimize(eps, maxIter).getCoordinates(x);
435 statTwo = monomerMinEngine.getStatus();
436 }
437 for(Atom a: s1Atoms){ a.setUse(false); }
438 logger.info("\n --------- System 2 Energy Breakdown --------- ");
439 double monomerEnergy2 = forceFieldEnergy.energy(x, true);
440 for(Atom a: s1Atoms){ a.setUse(true); }
441
442
443 logger.info(format("\n %-29s%12.7f kcal/mol", "Monomer energy 1:", monomerEnergy));
444 logger.info(format(" %-29s%12.7f kcal/mol", "Monomer energy 2:", monomerEnergy2));
445
446 m1MinEnergy = monomerEnergy;
447 m2MinEnergy = monomerEnergy2;
448 totalMonomerMinimizedEnergy = monomerEnergy + monomerEnergy2;
449
450 if(statOne != -1 && statTwo != -1){
451 logger.info("\n --------- Monomer Minimization Converged --------- ");
452 return 0;
453 } else{
454 logger.warning("\n --------- Monomer Minimization Did Not Converge --------- ");
455
456 logger.info(" Saving state for reference. ");
457 statesQueue.add(new StateContainer(new AssemblyState(mola), totalMonomerMinimizedEnergy));
458 return -1;
459 }
460 }
461 private int minimizeSystem(Atom a, Atom b) throws Exception{
462 if(tScan){
463
464 logger.info("\n --------- System 1 Static Torsion Scan --------- ");
465 forceFieldEnergy.getCoordinates(x);
466 double tscanE = forceFieldEnergy.energy(x, false);
467 for(Molecule m: s1) {
468 TorsionSearch m1TorsionSearch = new TorsionSearch(mola, m, 32, 1);
469 m1TorsionSearch.staticAnalysis(0, 100);
470 if (!m1TorsionSearch.getStates().isEmpty()) {
471 AssemblyState minState = m1TorsionSearch.getStates().get(0);
472 minState.revertState();
473 }
474 }
475 forceFieldEnergy.getCoordinates(x);
476 double tscanEAfter = forceFieldEnergy.energy(x, false);
477 logger.info("\n Energy before static torsion scan of system 1: " + tscanE);
478 logger.info(" Energy after static torsion scan of system 1: " + tscanEAfter);
479
480 logger.info("\n --------- System 2 Static Torsion Scan --------- ");
481 forceFieldEnergy.getCoordinates(x);
482 tscanE = forceFieldEnergy.energy(x, false);
483 for(Molecule m: s2) {
484 TorsionSearch m2TorsionSearch = new TorsionSearch(mola, m, 32, 1);
485 m2TorsionSearch.staticAnalysis(0, 100);
486 if (!m2TorsionSearch.getStates().isEmpty()) {
487 AssemblyState minState = m2TorsionSearch.getStates().get(0);
488 minState.revertState();
489 }
490 }
491 forceFieldEnergy.getCoordinates(x);
492 tscanEAfter = forceFieldEnergy.energy(x, false);
493 logger.info("\n Energy before static torsion scan of system 2: " + tscanE);
494 logger.info(" Energy after static torsion scan of system 2: " + tscanEAfter);
495 }
496 forceFieldEnergy.getCoordinates(x);
497 double e = forceFieldEnergy.energy(x, true);
498 RestrainDistance restrainDistance = getRestraintBond(a, b, e);
499 Minimize minEngine = new Minimize(mola, forceFieldEnergy, null);
500 try {
501 minEngine.minimize(this.eps, this.maxIter);
502 } catch (Exception ex){
503
504 a.getBonds().remove(restrainDistance);
505 b.getBonds().remove(restrainDistance);
506 a.update();
507 b.update();
508 mola.getBondList().remove(restrainDistance);
509 mola.update();
510 return -1;
511 }
512
513 a.getBonds().remove(restrainDistance);
514 b.getBonds().remove(restrainDistance);
515 a.update();
516 b.update();
517 mola.getBondList().remove(restrainDistance);
518 mola.update();
519 return minEngine.getStatus();
520 }
521
522 private RestrainDistance getRestraintBond(Atom a, Atom b, double e) throws Exception {
523 if (e > 1000000){
524 throw new Exception(" Energy too high to minimize.");
525 }
526
527 BondType restraint = new BondType(new int[]{a.getAtomicNumber(), b.getAtomicNumber()},
528 100.0,
529 this.hBondDist,
530 BondType.BondFunction.FLAT_BOTTOM_QUARTIC,
531 this.flatBottomRadius);
532 RestrainDistance restrainDistance = new RestrainDistance(a, b,
533 null,
534 false,
535 0.0, 0.0,
536 null);
537 restrainDistance.setBondType(restraint);
538 return restrainDistance;
539 }
540
541 private void setTargetAtoms(Atom[] atoms){
542 for(Atom a: atoms){
543 if(a.getAtomType().atomicNumber == 7|| a.getAtomType().atomicNumber == 8
544 || a.getAtomType().atomicNumber == 9 || a.getAtomType().atomicNumber == 15
545 || a.getAtomType().atomicNumber == 16 || a.getAtomType().atomicNumber == 17){
546 if(a.getMoleculeNumber() == mola.getMoleculeNumbers()[0]) {
547 s1TargetAtoms.add(a);
548 } else{
549 s2TargetAtoms.add(a);
550 }
551
552
553 if(excludeH) {
554 for (Bond b : a.getBonds()) {
555 int num = b.get1_2(a).getAtomType().atomicNumber;
556 if (num == 1) {
557 if (a.getMoleculeNumber() == mola.getMoleculeNumbers()[0]) {
558 s1TargetAtoms.add(a);
559 } else {
560 s2TargetAtoms.add(a);
561 }
562 }
563 }
564 }
565 }
566 else if(a.getAtomType().atomicNumber == 1 && !excludeH){
567 if(a.getMoleculeNumber() == mola.getMoleculeNumbers()[0]) {
568 s1TargetAtoms.add(a);
569 } else{
570 s2TargetAtoms.add(a);
571 }
572 }
573 }
574
575 }
576 private void calculateMeanStd(){
577 ArrayList<Double> energies = this.getEnergies();
578 double highOutlierCutoff = Double.MAX_VALUE;
579 double lowOutlierCutoff = Double.MIN_VALUE;
580
581 if(energies.size() > 4) {
582 double iqr = energies.get(energies.size() / 4) - energies.get(3 * energies.size() / 4);
583 double q3 = energies.get(3 * energies.size() / 4);
584 double q1 = energies.get(energies.size() / 4);
585 highOutlierCutoff = q3 + 1.5 * iqr;
586 lowOutlierCutoff = q1 - 1.5 * iqr;
587 }
588
589
590 int count = 0;
591 double sum = 0.0;
592 double sumNoOutlier = 0.0;
593 for(double e: energies){
594 if(e < highOutlierCutoff && e > lowOutlierCutoff){
595 sumNoOutlier += e;
596 count++;
597 }
598 sum += e;
599 }
600 averageEnergy = sum / energies.size();
601 averageEnergyNoOutlier = sumNoOutlier / count;
602
603
604 double sumOfSquares = 0.0;
605 double sumOfSquaresNoOutlier = 0.0;
606 for(double e: energies){
607 sumOfSquares += (e - averageEnergy) * (e - averageEnergy);
608 if(e < highOutlierCutoff && e > lowOutlierCutoff){
609 sumOfSquaresNoOutlier += (e - averageEnergyNoOutlier) * (e - averageEnergyNoOutlier);
610 }
611 }
612 stdOfEnergies = sqrt(sumOfSquares / energies.size());
613 stdOfEnergiesNoOutlier = sqrt(sumOfSquaresNoOutlier / count);
614 }
615
616 private static void alignSystemCOMtoAtomVecWithAxis(Atom a, double[] axis, Atom[] mAtoms){
617
618 double[] moleculeOneCOM = getCOM(mAtoms);
619 for(int i = 0; i < mAtoms.length; i++){
620 mAtoms[i].move(new double[] {-moleculeOneCOM[0], -moleculeOneCOM[1], -moleculeOneCOM[2]});
621 }
622
623 double[] aCoords = a.getXYZ().copy().get();
624
625 double[][] rotation = getRotationBetween(aCoords, axis);
626
627 double[] moleculeAtomicPositions = new double[mAtoms.length * 3];
628 for(int i = 0; i < mAtoms.length; i++){
629 moleculeAtomicPositions[i*3] = mAtoms[i].getX();
630 moleculeAtomicPositions[i*3 + 1] = mAtoms[i].getY();
631 moleculeAtomicPositions[i*3 + 2] = mAtoms[i].getZ();
632 }
633 applyRotation(moleculeAtomicPositions, rotation);
634
635 for(int i = 0; i < mAtoms.length; i++){
636 mAtoms[i].setXYZ(new double[] {moleculeAtomicPositions[i*3],
637 moleculeAtomicPositions[i*3 + 1],
638 moleculeAtomicPositions[i*3 + 2]});
639 }
640 }
641
642
643
644
645
646
647 private static double[] getCOM(Atom[] atoms){
648
649 double[] COM = new double[3];
650 double totalMass = 0.0;
651 for(Atom s: atoms){
652 double[] pos = s.getXYZ().get();
653 COM[0] += pos[0] * s.getMass();
654 COM[1] += pos[1] * s.getMass();
655 COM[2] += pos[2] * s.getMass();
656 totalMass += s.getMass();
657 }
658 totalMass = 1 / totalMass;
659 COM[0] *= totalMass;
660 COM[1] *= totalMass;
661 COM[2] *= totalMass;
662
663 return COM;
664 }
665
666
667
668
669
670
671
672 private static double[][] getRotationBetween(double[] v1, double[] v2){
673
674 double v1Norm = 1/ sqrt(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]);
675 double v2Norm = 1 / sqrt(v2[0] * v2[0] + v2[1] * v2[1] + v2[2] * v2[2]);
676 for(int i = 0; i < 3; i++){
677 v1[i] *= v1Norm;
678 v2[i] *= v2Norm;
679 }
680
681 double[] crossProduct = new double[3];
682 crossProduct[0] = v1[1] * v2[2] - v1[2] * v2[1];
683 crossProduct[1] = v1[2] * v2[0] - v1[0] * v2[2];
684 crossProduct[2] = v1[0] * v2[1] - v1[1] * v2[0];
685
686 double crossProductNorm = 1 / sqrt(crossProduct[0] * crossProduct[0] + crossProduct[1] * crossProduct[1] + crossProduct[2] * crossProduct[2]);
687 for(int i = 0; i < 3; i++){
688 crossProduct[i] *= crossProductNorm;
689 }
690
691 double dotProduct = v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
692
693 double theta = acos(dotProduct);
694 double[] u = crossProduct;
695
696 double[] quaternion = new double[4];
697 quaternion[0] = cos(theta/2);
698 quaternion[1] = u[0] * sin(theta/2);
699 quaternion[2] = u[1] * sin(theta/2);
700 quaternion[3] = u[2] * sin(theta/2);
701
702 double quaternionNorm = 1 / sqrt(quaternion[0] * quaternion[0] + quaternion[1] * quaternion[1]
703 + quaternion[2] * quaternion[2] + quaternion[3] * quaternion[3]);
704 for(int i = 0; i < 4; i++){ quaternion[i] *= quaternionNorm; }
705
706 double q1q1 = quaternion[1] * quaternion[1];
707 double q2q2 = quaternion[2] * quaternion[2];
708 double q3q3 = quaternion[3] * quaternion[3];
709 double q0q1 = quaternion[0] * quaternion[1];
710 double q0q2 = quaternion[0] * quaternion[2];
711 double q0q3 = quaternion[0] * quaternion[3];
712 double q1q2 = quaternion[1] * quaternion[2];
713 double q1q3 = quaternion[1] * quaternion[3];
714 double q2q3 = quaternion[2] * quaternion[3];
715
716 double[][] rotation = new double[3][3];
717 rotation[0][0] = 1 - 2 * (q2q2 + q3q3);
718 rotation[0][1] = 2 * (q1q2 - q0q3);
719 rotation[0][2] = 2 * (q1q3 + q0q2);
720 rotation[1][0] = 2 * (q1q2 + q0q3);
721 rotation[1][1] = 1 - 2 * (q1q1 + q3q3);
722 rotation[1][2] = 2 * (q2q3 - q0q1);
723 rotation[2][0] = 2 * (q1q3 - q0q2);
724 rotation[2][1] = 2 * (q2q3 + q0q1);
725 rotation[2][2] = 1 - 2 * (q1q1 + q2q2);
726
727 return rotation;
728 }
729
730 private static Atom[] getAtomsFromMoleculeArray(Molecule[] system){
731 ArrayList<Atom> atoms = new ArrayList<>();
732 for(Molecule m: system){
733 atoms.addAll(m.getAtomList());
734 }
735 return atoms.toArray(new Atom[0]);
736 }
737
738
739
740
741
742 private record StateContainer(AssemblyState state, double e) implements Comparable<StateContainer>
743 {
744 AssemblyState getState() {
745 return state;
746 }
747
748 double getEnergy() {
749 return e;
750 }
751
752 @Override
753 public int compareTo(StateContainer o) {
754 return Double.compare(e, o.getEnergy());
755 }
756 }
757 }