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 package ffx.potential;
39
40 import edu.rit.pj.ParallelRegion;
41 import edu.rit.pj.ParallelSection;
42 import edu.rit.pj.ParallelTeam;
43 import ffx.crystal.Crystal;
44 import ffx.crystal.CrystalPotential;
45 import ffx.crystal.SymOp;
46 import ffx.numerics.Potential;
47 import ffx.numerics.switching.UnivariateSwitchingFunction;
48 import ffx.potential.bonded.Atom;
49 import ffx.potential.bonded.LambdaInterface;
50 import ffx.potential.parameters.ForceField;
51 import ffx.potential.utils.EnergyException;
52
53 import java.util.ArrayList;
54 import java.util.Arrays;
55 import java.util.List;
56 import java.util.logging.Level;
57 import java.util.logging.Logger;
58
59 import static ffx.crystal.SymOp.applyCartesianSymOp;
60 import static ffx.crystal.SymOp.applyCartesianSymRot;
61 import static ffx.crystal.SymOp.invertSymOp;
62 import static ffx.potential.utils.Superpose.rmsd;
63 import static ffx.utilities.StringUtils.parseAtomRanges;
64 import static java.lang.Double.parseDouble;
65 import static java.lang.String.format;
66 import static java.util.Arrays.fill;
67 import static org.apache.commons.math3.util.FastMath.sqrt;
68
69
70
71
72
73
74
75 public class DualTopologyEnergy implements CrystalPotential, LambdaInterface {
76
77
78
79
80 private static final Logger logger = Logger.getLogger(DualTopologyEnergy.class.getName());
81
82
83
84 private final int nShared;
85
86
87
88 private final int nVariables;
89
90
91
92 private final EnergyRegion energyRegion;
93
94
95
96 private final double[] mass;
97
98
99
100 private final VARIABLE_TYPE[] variableTypes;
101
102
103
104 private final double[] x1;
105
106
107
108 private final double[] x2;
109
110
111
112 private final double[] g1;
113
114
115
116 private final double[] g2;
117
118
119
120 private final double[] rg1;
121
122
123
124 private final double[] rg2;
125
126
127
128 private final double[] gl1;
129
130
131
132 private final double[] gl2;
133
134
135
136
137 private final double[] rgl1;
138
139
140
141
142 private final double[] rgl2;
143
144
145
146 private final CrystalPotential potential1;
147
148
149
150 private final CrystalPotential potential2;
151
152
153
154 private final LambdaInterface lambdaInterface1;
155
156
157
158 private final LambdaInterface lambdaInterface2;
159
160
161
162 private final ForceFieldEnergy forceFieldEnergy1;
163
164
165
166 private final ForceFieldEnergy forceFieldEnergy2;
167
168
169
170 private final int nActive1;
171
172
173
174 private final int nActive2;
175
176
177
178 private final Atom[] activeAtoms1;
179
180
181
182 private final Atom[] activeAtoms2;
183
184
185
186 private final boolean[] sharedAtoms1;
187
188
189
190 private final boolean[] sharedAtoms2;
191
192
193
194 private final UnivariateSwitchingFunction switchFunction;
195
196
197
198 private double energy1 = 0;
199
200
201
202 private double energy2 = 0;
203
204
205
206 private ParallelTeam parallelTeam;
207
208
209
210 private final boolean doValenceRestraint1;
211
212
213
214 private final boolean doValenceRestraint2;
215
216
217
218 private final boolean useFirstSystemBondedEnergy;
219
220
221
222 private double restraintEnergy1 = 0;
223
224
225
226 private double restraintEnergy2 = 0;
227
228
229
230 private double dEdL_1 = 0;
231
232
233
234 private double dEdL_2 = 0;
235
236
237
238 private double restraintdEdL_1 = 0;
239
240
241
242 private double restraintdEdL_2 = 0;
243
244
245
246 private double d2EdL2_1 = 0;
247
248
249
250 private double d2EdL2_2 = 0;
251
252
253
254 private double restraintd2EdL2_1 = 0;
255
256
257
258 private double restraintd2EdL2_2 = 0;
259
260
261
262 private double totalEnergy = 0;
263
264
265
266 private double lambda = 1.0;
267
268
269
270 private double f1L = 1.0;
271
272
273
274 private double f2L = 0.0;
275
276
277
278 private double dF1dL = 0.0;
279
280
281
282 private double dF2dL = 0.0;
283
284
285
286 private double d2F1dL2 = 0.0;
287
288
289
290 private double d2F2dL2 = 0.0;
291
292
293
294 private double[] scaling = null;
295
296
297
298 private STATE state = STATE.BOTH;
299
300
301
302 private SymOp[] symOp;
303
304
305
306 private SymOp[] inverse;
307
308
309
310 private final ArrayList<List<Integer>> symOpAtoms = new ArrayList<>();
311
312
313
314 private final int[] molecule2;
315
316
317
318 private boolean useSymOp = false;
319
320
321
322
323
324
325
326
327 public DualTopologyEnergy(
328 MolecularAssembly topology1,
329 MolecularAssembly topology2,
330 UnivariateSwitchingFunction switchFunction) {
331 forceFieldEnergy1 = topology1.getPotentialEnergy();
332 forceFieldEnergy2 = topology2.getPotentialEnergy();
333 potential1 = forceFieldEnergy1;
334 potential2 = forceFieldEnergy2;
335 lambdaInterface1 = forceFieldEnergy1;
336 lambdaInterface2 = forceFieldEnergy2;
337
338
339 Atom[] atoms1 = topology1.getAtomArray();
340
341 Atom[] atoms2 = topology2.getAtomArray();
342
343 molecule2 = topology2.getMoleculeNumbers();
344
345 ForceField forceField1 = topology1.getForceField();
346 doValenceRestraint1 = forceField1.getBoolean("LAMBDA_VALENCE_RESTRAINTS", true);
347 ForceField forceField2 = topology2.getForceField();
348 doValenceRestraint2 = forceField2.getBoolean("LAMBDA_VALENCE_RESTRAINTS", true);
349
350 useFirstSystemBondedEnergy = forceField2.getBoolean("USE_FIRST_SYSTEM_BONDED_ENERGY", false);
351
352
353 int shared1 = 0;
354 int shared2 = 0;
355 int activeCount1 = 0;
356 int activeCount2 = 0;
357 for (Atom a1 : atoms1) {
358 if (a1.isActive()) {
359 activeCount1++;
360 if (!a1.applyLambda()) {
361 shared1++;
362 }
363 }
364 }
365 for (Atom a2 : atoms2) {
366 if (a2.isActive()) {
367 activeCount2++;
368 if (!a2.applyLambda()) {
369 shared2++;
370 }
371 }
372 }
373 nActive1 = activeCount1;
374 nActive2 = activeCount2;
375 activeAtoms1 = new Atom[activeCount1];
376 activeAtoms2 = new Atom[activeCount2];
377 sharedAtoms1 = new boolean[activeCount1];
378 sharedAtoms2 = new boolean[activeCount2];
379 Arrays.fill(sharedAtoms1, true);
380 Arrays.fill(sharedAtoms2, true);
381 int index = 0;
382 for (Atom a1 : atoms1) {
383 if (a1.isActive()) {
384 activeAtoms1[index] = a1;
385 if (a1.applyLambda()) {
386
387 sharedAtoms1[index] = false;
388 }
389 index++;
390 }
391 }
392 index = 0;
393 for (Atom a2 : atoms2) {
394 if (a2.isActive()) {
395 activeAtoms2[index] = a2;
396 if (a2.applyLambda()) {
397
398 sharedAtoms2[index] = false;
399 }
400 index++;
401 }
402 }
403
404 assert (shared1 == shared2);
405 nShared = shared1;
406
407 int nSoftCore1 = nActive1 - nShared;
408
409 int nSoftCore2 = nActive2 - nShared;
410
411 int nTotal = nShared + nSoftCore1 + nSoftCore2;
412 nVariables = 3 * nTotal;
413
414
415 x1 = new double[nActive1 * 3];
416 x2 = new double[nActive2 * 3];
417 g1 = new double[nActive1 * 3];
418 g2 = new double[nActive2 * 3];
419 rg1 = new double[nActive1 * 3];
420 rg2 = new double[nActive2 * 3];
421 gl1 = new double[nActive1 * 3];
422 gl2 = new double[nActive2 * 3];
423 rgl1 = new double[nActive1 * 3];
424 rgl2 = new double[nActive2 * 3];
425
426
427 int i1 = 0;
428 int i2 = 0;
429 logger.info(format(" Shared atoms: %d (1: %d of %d; 2: %d of %d)\n ", nShared, shared1, nActive1, shared2, nActive2));
430 for (int i = 0; i < nShared; i++) {
431 Atom a1 = atoms1[i1++];
432 while (a1.applyLambda()) {
433 a1 = atoms1[i1++];
434 }
435 Atom a2 = atoms2[i2++];
436 while (a2.applyLambda()) {
437 a2 = atoms2[i2++];
438 }
439
440
441
442
443
444 }
445
446
447 index = 0;
448 variableTypes = new VARIABLE_TYPE[nVariables];
449 for (int i = 0; i < nTotal; i++) {
450 variableTypes[index++] = VARIABLE_TYPE.X;
451 variableTypes[index++] = VARIABLE_TYPE.Y;
452 variableTypes[index++] = VARIABLE_TYPE.Z;
453 }
454
455
456 int commonIndex = 0;
457 int softcoreIndex = 3 * nShared;
458 mass = new double[nVariables];
459 for (int i = 0; i < nActive1; i++) {
460 Atom a = activeAtoms1[i];
461 double m = a.getMass();
462 if (sharedAtoms1[i]) {
463 mass[commonIndex++] = m;
464 mass[commonIndex++] = m;
465 mass[commonIndex++] = m;
466 } else {
467 mass[softcoreIndex++] = m;
468 mass[softcoreIndex++] = m;
469 mass[softcoreIndex++] = m;
470 }
471 }
472 commonIndex = 0;
473 for (int i = 0; i < nActive2; i++) {
474 Atom a = activeAtoms2[i];
475 double m = a.getMass();
476 if (sharedAtoms2[i]) {
477 for (int j = 0; j < 3; j++) {
478 double massCommon = mass[commonIndex];
479 massCommon = Math.max(massCommon, m);
480 mass[commonIndex++] = massCommon;
481 }
482 } else {
483 mass[softcoreIndex++] = m;
484 mass[softcoreIndex++] = m;
485 mass[softcoreIndex++] = m;
486 }
487 }
488 energyRegion = new EnergyRegion();
489 parallelTeam = new ParallelTeam(1);
490
491 String symOpString = forceField2.getString("symop", null);
492
493 if (symOpString != null) {
494
495 readSymOp(symOpString);
496 } else {
497
498 potential2.setCrystal(potential1.getCrystal());
499 useSymOp = false;
500 symOp = null;
501 inverse = null;
502 }
503
504 this.switchFunction = switchFunction;
505 logger.info(format("\n Dual topology using switching function:\n %s", switchFunction));
506 }
507
508
509
510
511
512
513 private void readSymOp(String symOpString) {
514 int numSymOps = 0;
515 try {
516 String[] tokens = symOpString.split(" +");
517 int numTokens = tokens.length;
518
519 if (numTokens == 12) {
520 numSymOps = 1;
521 symOp = new SymOp[numSymOps];
522 inverse = new SymOp[numSymOps];
523
524 symOpAtoms.add(parseAtomRanges("symmetryAtoms", "1-N", nActive2));
525
526 symOp[0] = new SymOp(new double[][]{
527 {parseDouble(tokens[0]), parseDouble(tokens[1]), parseDouble(tokens[2])},
528 {parseDouble(tokens[3]), parseDouble(tokens[4]), parseDouble(tokens[5])},
529 {parseDouble(tokens[6]), parseDouble(tokens[7]), parseDouble(tokens[8])}},
530 new double[]{parseDouble(tokens[9]), parseDouble(tokens[10]), parseDouble(tokens[11])});
531 } else if (numTokens % 14 == 0) {
532 numSymOps = numTokens / 14;
533 symOp = new SymOp[numSymOps];
534 inverse = new SymOp[numSymOps];
535 if (logger.isLoggable(Level.FINER)) {
536 logger.finer(format(" Number of Tokens: %3d Number of Sym Ops: %2d", numTokens, numSymOps));
537 }
538 for (int i = 0; i < numSymOps; i++) {
539 int j = i * 14;
540
541 symOpAtoms.add(parseAtomRanges("symmetryAtoms", tokens[j], nActive2));
542 symOp[i] = new SymOp(new double[][]{
543 {parseDouble(tokens[j + 2]), parseDouble(tokens[j + 3]), parseDouble(tokens[j + 4])},
544 {parseDouble(tokens[j + 5]), parseDouble(tokens[j + 6]), parseDouble(tokens[j + 7])},
545 {parseDouble(tokens[j + 8]), parseDouble(tokens[j + 9]), parseDouble(tokens[j + 10])}},
546 new double[]{parseDouble(tokens[j + 11]), parseDouble(tokens[j + 12]), parseDouble(tokens[j + 13])});
547 }
548 } else {
549 logger.warning(" Sym Op is formatted incorrectly.");
550 return;
551 }
552 useSymOp = true;
553
554
555 potential1.getCoordinates(x1);
556 potential2.getCoordinates(x2);
557
558
559 double[] x1Shared = new double[nShared * 3];
560 double[] x2Shared = new double[nShared * 3];
561 int sharedIndex = 0;
562 for (int i = 0; i < nActive1; i++) {
563 if (sharedAtoms1[i]) {
564 int index = i * 3;
565 x1Shared[sharedIndex++] = x1[index];
566 x1Shared[sharedIndex++] = x1[index + 1];
567 x1Shared[sharedIndex++] = x1[index + 2];
568 }
569 }
570 sharedIndex = 0;
571 for (int i = 0; i < nActive2; i++) {
572 if (sharedAtoms2[i]) {
573 int index = i * 3;
574 x2Shared[sharedIndex++] = x2[index];
575 x2Shared[sharedIndex++] = x2[index + 1];
576 x2Shared[sharedIndex++] = x2[index + 2];
577 }
578 }
579
580 double[] m = new double[nShared];
581 Arrays.fill(m, 1.0);
582 double origRMSD = rmsd(x1Shared, x2Shared, m);
583 logger.info("\n Topology 2 Coordinates from Input File");
584 logger.info(format(" RMSD: %12.3f", origRMSD));
585
586
587 double[] origX1 = Arrays.copyOf(x1Shared, nShared * 3);
588 double[] origX2 = Arrays.copyOf(x2Shared, nShared * 3);
589
590
591
592 logger.info(" Number Sym Ops: " + numSymOps);
593 for (int i = 0; i < numSymOps; i++) {
594 List<Integer> symAtoms = symOpAtoms.get(i);
595 boolean[] mask = new boolean[nShared];
596 int index = 0;
597 for (int j = 0; j < nActive2; j++) {
598 if (sharedAtoms2[j]) {
599 if (symAtoms.contains(j)) {
600 mask[index++] = true;
601 } else {
602 index++;
603 }
604 }
605 }
606
607 inverse[i] = invertSymOp(symOp[i]);
608
609 applyCartesianSymOp(origX1, origX1, symOp[i], mask);
610
611 logger.info("\n SymOp between topologies:\n" + symOp[i]);
612 logger.info("\n Inverse SymOp:\n" + inverse[i]);
613 }
614 double symOpRMSD = rmsd(origX1, origX2, m);
615 logger.info("\n Topology 2 Coordinates from Loaded SymOp");
616 logger.info(format(" RMSD: %12.3f", symOpRMSD));
617
618 if (logger.isLoggable(Level.FINE)) {
619 for (int i = 0; i < nShared; i++) {
620 int i3 = i * 3;
621 double rmsd = rmsd(new double[]{origX1[i3 + 0], origX1[i3 + 1], origX1[i3 + 2]}, new double[]{origX2[i3 + 0], origX2[i3 + 1], origX2[i3 + 2]}, new double[]{1.0});
622 if (rmsd > 0.5) {
623 logger.fine(format(" Atom %3d has an RMSD of %8.3f", i + 1, rmsd));
624 }
625 }
626 }
627 } catch (Exception ex) {
628 logger.warning(ex.toString());
629 logger.severe(" Error parsing SymOp for Dual Topology:\n (" + symOpString + ")");
630 }
631 }
632
633
634
635
636
637
638
639 @Override
640 public boolean dEdLZeroAtEnds() {
641 if (!forceFieldEnergy1.dEdLZeroAtEnds() || !forceFieldEnergy2.dEdLZeroAtEnds()) {
642 return false;
643 }
644 return (switchFunction.highestOrderZeroDerivativeAtZeroBound() > 0);
645 }
646
647
648
649
650 @Override
651 public boolean destroy() {
652 boolean ffe1Destroy = forceFieldEnergy1.destroy();
653 boolean ffe2Destroy = forceFieldEnergy2.destroy();
654 try {
655 if (parallelTeam != null) {
656 parallelTeam.shutdown();
657 }
658 return ffe1Destroy && ffe2Destroy;
659 } catch (Exception ex) {
660 logger.warning(format(" Exception in shutting down DualTopologyEnergy: %s", ex));
661 logger.info(Utilities.stackTraceToString(ex));
662 return false;
663 }
664 }
665
666
667
668
669 @Override
670 public double energy(double[] x) {
671 return energy(x, false);
672 }
673
674
675
676
677 @Override
678 public double energy(double[] x, boolean verbose) {
679 try {
680 energyRegion.setX(x);
681 energyRegion.setVerbose(verbose);
682 parallelTeam.execute(energyRegion);
683 } catch (Exception ex) {
684 throw new EnergyException(format(" Exception in calculating dual-topology energy: %s", ex));
685 }
686 return totalEnergy;
687 }
688
689
690
691
692
693
694 @Override
695 public double energyAndGradient(double[] x, double[] g) {
696 return energyAndGradient(x, g, false);
697 }
698
699
700
701
702
703
704 @Override
705 public double energyAndGradient(double[] x, double[] g, boolean verbose) {
706 assert Arrays.stream(x).allMatch(Double::isFinite);
707 try {
708 energyRegion.setX(x);
709 energyRegion.setG(g);
710 energyRegion.setVerbose(verbose);
711 parallelTeam.execute(energyRegion);
712 } catch (Exception ex) {
713 throw new EnergyException(format(" Exception in calculating dual-topology energy: %s", ex));
714 }
715 return totalEnergy;
716 }
717
718
719
720
721 @Override
722 public double[] getAcceleration(double[] acceleration) {
723 if (acceleration == null || acceleration.length < nVariables) {
724 acceleration = new double[nVariables];
725 }
726 int indexCommon = 0;
727 int indexUnique = nShared * 3;
728 double[] accel = new double[3];
729 for (int i = 0; i < nActive1; i++) {
730 Atom atom = activeAtoms1[i];
731 atom.getAcceleration(accel);
732 if (sharedAtoms1[i]) {
733 acceleration[indexCommon++] = accel[0];
734 acceleration[indexCommon++] = accel[1];
735 acceleration[indexCommon++] = accel[2];
736 } else {
737 acceleration[indexUnique++] = accel[0];
738 acceleration[indexUnique++] = accel[1];
739 acceleration[indexUnique++] = accel[2];
740 }
741 }
742 for (int i = 0; i < nActive2; i++) {
743 if (!sharedAtoms2[i]) {
744 Atom atom = activeAtoms2[i];
745 atom.getAcceleration(accel);
746 acceleration[indexUnique++] = accel[0];
747 acceleration[indexUnique++] = accel[1];
748 acceleration[indexUnique++] = accel[2];
749 }
750 }
751
752 return acceleration;
753 }
754
755
756
757
758 @Override
759 public double[] getCoordinates(double[] x) {
760 if (x == null) {
761 x = new double[nVariables];
762 }
763 int indexCommon = 0;
764 int indexUnique = nShared * 3;
765 for (int i = 0; i < nActive1; i++) {
766 Atom a = activeAtoms1[i];
767 if (sharedAtoms1[i]) {
768 x[indexCommon++] = a.getX();
769 x[indexCommon++] = a.getY();
770 x[indexCommon++] = a.getZ();
771 } else {
772 x[indexUnique++] = a.getX();
773 x[indexUnique++] = a.getY();
774 x[indexUnique++] = a.getZ();
775 }
776 }
777 for (int i = 0; i < nActive2; i++) {
778 if (!sharedAtoms2[i]) {
779 Atom a = activeAtoms2[i];
780 x[indexUnique++] = a.getX();
781 x[indexUnique++] = a.getY();
782 x[indexUnique++] = a.getZ();
783 }
784 }
785 return x;
786 }
787
788
789
790
791 @Override
792 public Crystal getCrystal() {
793
794 if (useSymOp) {
795 logger.warning(" Get method only returned first crystal.");
796 }
797 return potential1.getCrystal();
798 }
799
800
801
802
803 @Override
804 public void setCrystal(Crystal crystal) {
805
806 if (useSymOp) {
807 logger.warning(" Both systems set to the same crystal.");
808 }
809 potential1.setCrystal(crystal);
810 potential2.setCrystal(crystal);
811 }
812
813
814
815
816 @Override
817 public STATE getEnergyTermState() {
818 return state;
819 }
820
821
822
823
824 @Override
825 public void setEnergyTermState(STATE state) {
826 this.state = state;
827 potential1.setEnergyTermState(state);
828 potential2.setEnergyTermState(state);
829 }
830
831
832
833
834 @Override
835 public double getLambda() {
836 return lambda;
837 }
838
839
840
841
842 @Override
843 public void setLambda(double lambda) {
844 if (lambda <= 1.0 && lambda >= 0.0) {
845 this.lambda = lambda;
846 double oneMinusLambda = 1.0 - lambda;
847 lambdaInterface1.setLambda(lambda);
848 lambdaInterface2.setLambda(oneMinusLambda);
849
850 f1L = switchFunction.valueAt(lambda);
851 dF1dL = switchFunction.firstDerivative(lambda);
852 d2F1dL2 = switchFunction.secondDerivative(lambda);
853
854 f2L = switchFunction.valueAt(oneMinusLambda);
855 dF2dL = -1.0 * switchFunction.firstDerivative(oneMinusLambda);
856 d2F2dL2 = switchFunction.secondDerivative(oneMinusLambda);
857 } else {
858 String message = format("Lambda value %8.3f is not in the range [0..1].", lambda);
859 throw new IllegalArgumentException(message);
860 }
861 }
862
863
864
865
866 @Override
867 public double[] getMass() {
868 return mass;
869 }
870
871
872
873
874
875
876 public int getNumSharedVariables() {
877 return 3 * nShared;
878 }
879
880
881
882
883 @Override
884 public int getNumberOfVariables() {
885 return nVariables;
886 }
887
888
889
890
891 @Override
892 public double[] getPreviousAcceleration(double[] previousAcceleration) {
893 if (previousAcceleration == null || previousAcceleration.length < nVariables) {
894 previousAcceleration = new double[nVariables];
895 }
896 int indexCommon = 0;
897 int indexUnique = nShared * 3;
898 double[] prev = new double[3];
899 for (int i = 0; i < nActive1; i++) {
900 Atom atom = activeAtoms1[i];
901 atom.getPreviousAcceleration(prev);
902 if (sharedAtoms1[i]) {
903 previousAcceleration[indexCommon++] = prev[0];
904 previousAcceleration[indexCommon++] = prev[1];
905 previousAcceleration[indexCommon++] = prev[2];
906 } else {
907 previousAcceleration[indexUnique++] = prev[0];
908 previousAcceleration[indexUnique++] = prev[1];
909 previousAcceleration[indexUnique++] = prev[2];
910 }
911 }
912 for (int i = 0; i < nActive2; i++) {
913 Atom atom = activeAtoms2[i];
914 if (!sharedAtoms2[i]) {
915 atom.getPreviousAcceleration(prev);
916 previousAcceleration[indexUnique++] = prev[0];
917 previousAcceleration[indexUnique++] = prev[1];
918 previousAcceleration[indexUnique++] = prev[2];
919 }
920 }
921
922 return previousAcceleration;
923 }
924
925
926
927
928 @Override
929 public double[] getScaling() {
930 return scaling;
931 }
932
933
934
935
936 @Override
937 public void setScaling(double[] scaling) {
938 this.scaling = scaling;
939 }
940
941
942
943
944
945
946
947 public UnivariateSwitchingFunction getSwitch() {
948 return switchFunction;
949 }
950
951
952
953
954 @Override
955 public double getTotalEnergy() {
956 return totalEnergy;
957 }
958
959 @Override
960 public List<Potential> getUnderlyingPotentials() {
961 List<Potential> under = new ArrayList<>(2);
962 under.add(forceFieldEnergy1);
963 under.add(forceFieldEnergy2);
964 under.addAll(forceFieldEnergy1.getUnderlyingPotentials());
965 under.addAll(forceFieldEnergy2.getUnderlyingPotentials());
966 return under;
967 }
968
969
970
971
972 @Override
973 public VARIABLE_TYPE[] getVariableTypes() {
974 return variableTypes;
975 }
976
977
978
979
980 @Override
981 public double[] getVelocity(double[] velocity) {
982 if (velocity == null || velocity.length < nVariables) {
983 velocity = new double[nVariables];
984 }
985 int indexCommon = 0;
986 int indexUnique = nShared * 3;
987 double[] vel = new double[3];
988 for (int i = 0; i < nActive1; i++) {
989 Atom atom = activeAtoms1[i];
990 atom.getVelocity(vel);
991 if (sharedAtoms1[i]) {
992 velocity[indexCommon++] = vel[0];
993 velocity[indexCommon++] = vel[1];
994 velocity[indexCommon++] = vel[2];
995 } else {
996 velocity[indexUnique++] = vel[0];
997 velocity[indexUnique++] = vel[1];
998 velocity[indexUnique++] = vel[2];
999 }
1000 }
1001 for (int i = 0; i < nActive2; i++) {
1002 if (!sharedAtoms2[i]) {
1003 Atom atom = activeAtoms2[i];
1004 atom.getVelocity(vel);
1005 velocity[indexUnique++] = vel[0];
1006 velocity[indexUnique++] = vel[1];
1007 velocity[indexUnique++] = vel[2];
1008 }
1009 }
1010
1011 return velocity;
1012 }
1013
1014
1015
1016
1017 @Override
1018 public double getd2EdL2() {
1019 double e1 =
1020 f1L * d2EdL2_1
1021 + 2.0 * dF1dL * dEdL_1
1022 + d2F1dL2 * energy1
1023 + f2L * restraintd2EdL2_1
1024 + 2.0 * dF2dL * restraintdEdL_1
1025 + d2F2dL2 * restraintEnergy1;
1026 double e2;
1027 if (!useFirstSystemBondedEnergy) {
1028 e2 = f2L * d2EdL2_2
1029 + 2.0 * dF2dL * dEdL_2
1030 + d2F2dL2 * energy2
1031 + f1L * restraintd2EdL2_2
1032 + 2.0 * dF1dL * restraintdEdL_2
1033 + d2F1dL2 * restraintEnergy2;
1034 } else {
1035 e2 = f2L * d2EdL2_2
1036 + 2.0 * dF2dL * dEdL_2
1037 + d2F2dL2 * energy2
1038 - f2L * restraintd2EdL2_2
1039 - 2.0 * dF2dL * restraintdEdL_2
1040 - d2F2dL2 * restraintEnergy2;
1041 }
1042 return e1 + e2;
1043 }
1044
1045
1046
1047
1048 @Override
1049 public double getdEdL() {
1050
1051
1052 double e1 = f1L * dEdL_1 + dF1dL * energy1 + f2L * restraintdEdL_1 + dF2dL * restraintEnergy1;
1053 double e2;
1054 if (!useFirstSystemBondedEnergy) {
1055 e2 = f2L * dEdL_2 + dF2dL * energy2 + f1L * restraintdEdL_2 + dF1dL * restraintEnergy2;
1056 } else {
1057 e2 = f2L * dEdL_2 + dF2dL * energy2 - f2L * restraintdEdL_2 - dF2dL * restraintEnergy2;
1058 }
1059 return e1 + e2;
1060 }
1061
1062
1063
1064
1065 @Override
1066 public void getdEdXdL(double[] g) {
1067 if (g == null) {
1068 g = new double[nVariables];
1069 }
1070
1071 int index = 0;
1072 int indexCommon = 0;
1073 int indexUnique = nShared * 3;
1074
1075 for (int i = 0; i < nActive1; i++) {
1076 if (sharedAtoms1[i]) {
1077 g[indexCommon++] =
1078 f1L * gl1[index] + dF1dL * g1[index] + f2L * rgl1[index] + dF2dL * rg1[index++];
1079 g[indexCommon++] =
1080 f1L * gl1[index] + dF1dL * g1[index] + f2L * rgl1[index] + dF2dL * rg1[index++];
1081 g[indexCommon++] =
1082 f1L * gl1[index] + dF1dL * g1[index] + f2L * rgl1[index] + dF2dL * rg1[index++];
1083 } else {
1084 g[indexUnique++] =
1085 f1L * gl1[index] + dF1dL * g1[index] + f2L * rgl1[index] + dF2dL * rg1[index++];
1086 g[indexUnique++] =
1087 f1L * gl1[index] + dF1dL * g1[index] + f2L * rgl1[index] + dF2dL * rg1[index++];
1088 g[indexUnique++] =
1089 f1L * gl1[index] + dF1dL * g1[index] + f2L * rgl1[index] + dF2dL * rg1[index++];
1090 }
1091 }
1092
1093
1094 if (!useFirstSystemBondedEnergy) {
1095 index = 0;
1096 indexCommon = 0;
1097 for (int i = 0; i < nActive2; i++) {
1098 if (sharedAtoms2[i]) {
1099 g[indexCommon++] += (-f2L * gl2[index] + dF2dL * g2[index] - f1L * rgl2[index]
1100 + dF1dL * rg2[index++]);
1101 g[indexCommon++] += (-f2L * gl2[index] + dF2dL * g2[index] - f1L * rgl2[index]
1102 + dF1dL * rg2[index++]);
1103 g[indexCommon++] += (-f2L * gl2[index] + dF2dL * g2[index] - f1L * rgl2[index]
1104 + dF1dL * rg2[index++]);
1105 } else {
1106 g[indexUnique++] = (-f2L * gl2[index] + dF2dL * g2[index] - f1L * rgl2[index]
1107 + dF1dL * rg2[index++]);
1108 g[indexUnique++] = (-f2L * gl2[index] + dF2dL * g2[index] - f1L * rgl2[index]
1109 + dF1dL * rg2[index++]);
1110 g[indexUnique++] = (-f2L * gl2[index] + dF2dL * g2[index] - f1L * rgl2[index]
1111 + dF1dL * rg2[index++]);
1112 }
1113 }
1114 } else {
1115 index = 0;
1116 indexCommon = 0;
1117 for (int i = 0; i < nActive2; i++) {
1118 if (sharedAtoms2[i]) {
1119 g[indexCommon++] += (-f2L * gl2[index] + dF2dL * g2[index] + f2L * rgl2[index]
1120 - dF2dL * rg2[index++]);
1121 g[indexCommon++] += (-f2L * gl2[index] + dF2dL * g2[index] + f2L * rgl2[index]
1122 - dF2dL * rg2[index++]);
1123 g[indexCommon++] += (-f2L * gl2[index] + dF2dL * g2[index] + f2L * rgl2[index]
1124 - dF2dL * rg2[index++]);
1125 } else {
1126 g[indexUnique++] = (-f2L * gl2[index] + dF2dL * g2[index] + f2L * rgl2[index]
1127 - dF2dL * rg2[index++]);
1128 g[indexUnique++] = (-f2L * gl2[index] + dF2dL * g2[index] + f2L * rgl2[index]
1129 - dF2dL * rg2[index++]);
1130 g[indexUnique++] = (-f2L * gl2[index] + dF2dL * g2[index] + f2L * rgl2[index]
1131 - dF2dL * rg2[index++]);
1132 }
1133 }
1134
1135 }
1136 }
1137
1138
1139
1140
1141 @Override
1142 public void setAcceleration(double[] acceleration) {
1143 double[] accel = new double[3];
1144 int indexCommon = 0;
1145 int indexUnique = 3 * nShared;
1146 for (int i = 0; i < nActive1; i++) {
1147 Atom atom = activeAtoms1[i];
1148 if (sharedAtoms1[i]) {
1149 accel[0] = acceleration[indexCommon++];
1150 accel[1] = acceleration[indexCommon++];
1151 accel[2] = acceleration[indexCommon++];
1152 } else {
1153 accel[0] = acceleration[indexUnique++];
1154 accel[1] = acceleration[indexUnique++];
1155 accel[2] = acceleration[indexUnique++];
1156 }
1157 atom.setAcceleration(accel);
1158 }
1159 indexCommon = 0;
1160 for (int i = 0; i < nActive2; i++) {
1161 Atom atom = activeAtoms2[i];
1162 if (sharedAtoms2[i]) {
1163 accel[0] = acceleration[indexCommon++];
1164 accel[1] = acceleration[indexCommon++];
1165 accel[2] = acceleration[indexCommon++];
1166 } else {
1167 accel[0] = acceleration[indexUnique++];
1168 accel[1] = acceleration[indexUnique++];
1169 accel[2] = acceleration[indexUnique++];
1170 }
1171 atom.setAcceleration(accel);
1172 }
1173 }
1174
1175
1176
1177
1178
1179
1180 public void setParallel(boolean parallel) {
1181 if (parallelTeam != null) {
1182 try {
1183 parallelTeam.shutdown();
1184 } catch (Exception e) {
1185 logger.severe(format(" Exception in shutting down old ParallelTeam for DualTopologyEnergy: %s", e));
1186 }
1187 }
1188 parallelTeam = parallel ? new ParallelTeam(2) : new ParallelTeam(1);
1189 }
1190
1191
1192
1193
1194 @Override
1195 public void setPreviousAcceleration(double[] previousAcceleration) {
1196 double[] prev = new double[3];
1197 int indexCommon = 0;
1198 int indexUnique = 3 * nShared;
1199 for (int i = 0; i < nActive1; i++) {
1200 Atom atom = activeAtoms1[i];
1201 if (sharedAtoms1[i]) {
1202 prev[0] = previousAcceleration[indexCommon++];
1203 prev[1] = previousAcceleration[indexCommon++];
1204 prev[2] = previousAcceleration[indexCommon++];
1205 } else {
1206 prev[0] = previousAcceleration[indexUnique++];
1207 prev[1] = previousAcceleration[indexUnique++];
1208 prev[2] = previousAcceleration[indexUnique++];
1209 }
1210 atom.setPreviousAcceleration(prev);
1211 }
1212 indexCommon = 0;
1213 for (int i = 0; i < nActive2; i++) {
1214 Atom atom = activeAtoms2[i];
1215 if (sharedAtoms2[i]) {
1216 prev[0] = previousAcceleration[indexCommon++];
1217 prev[1] = previousAcceleration[indexCommon++];
1218 prev[2] = previousAcceleration[indexCommon++];
1219 } else {
1220 prev[0] = previousAcceleration[indexUnique++];
1221 prev[1] = previousAcceleration[indexUnique++];
1222 prev[2] = previousAcceleration[indexUnique++];
1223 }
1224 atom.setPreviousAcceleration(prev);
1225 }
1226 }
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238 public void setPrintOnFailure(boolean onFail, boolean override) {
1239 forceFieldEnergy1.setPrintOnFailure(onFail, override);
1240 forceFieldEnergy2.setPrintOnFailure(onFail, override);
1241 }
1242
1243
1244
1245
1246 @Override
1247 public void setVelocity(double[] velocity) {
1248 double[] vel = new double[3];
1249 int indexCommon = 0;
1250 int indexUnique = 3 * nShared;
1251 for (int i = 0; i < nActive1; i++) {
1252 Atom atom = activeAtoms1[i];
1253 if (sharedAtoms1[i]) {
1254 vel[0] = velocity[indexCommon++];
1255 vel[1] = velocity[indexCommon++];
1256 vel[2] = velocity[indexCommon++];
1257 } else {
1258 vel[0] = velocity[indexUnique++];
1259 vel[1] = velocity[indexUnique++];
1260 vel[2] = velocity[indexUnique++];
1261 }
1262 atom.setVelocity(vel);
1263 }
1264
1265 indexCommon = 0;
1266 for (int i = 0; i < nActive2; i++) {
1267 Atom atom = activeAtoms2[i];
1268 if (sharedAtoms2[i]) {
1269 vel[0] = velocity[indexCommon++];
1270 vel[1] = velocity[indexCommon++];
1271 vel[2] = velocity[indexCommon++];
1272 } else {
1273 vel[0] = velocity[indexUnique++];
1274 vel[1] = velocity[indexUnique++];
1275 vel[2] = velocity[indexUnique++];
1276 }
1277 atom.setVelocity(vel);
1278 }
1279 }
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289 private void reconcileAtoms(Atom a1, Atom a2, Level warnlev) {
1290 double dist = 0;
1291 double[] xyz1 = a1.getXYZ(null);
1292 double[] xyz2 = a2.getXYZ(null);
1293 double[] xyzAv = new double[3];
1294 for (int i = 0; i < 3; i++) {
1295 double dx = xyz1[i] - xyz2[i];
1296 dist += (dx * dx);
1297 xyzAv[i] = xyz1[i] + (0.5 * dx);
1298 }
1299
1300
1301 double maxDist2 = 0.09;
1302
1303
1304
1305
1306
1307 double minDistWarn2 = 0.00001;
1308
1309 if (dist > maxDist2) {
1310 logger.log(
1311 Level.SEVERE,
1312 format(
1313 " Distance between atoms %s " + "and %s is %7.4f >> maximum allowed %7.4f",
1314 a1, a2, sqrt(dist), sqrt(maxDist2)));
1315 } else if (dist > minDistWarn2) {
1316 logger.log(
1317 warnlev,
1318 format(
1319 " Distance between atoms %s " + "and %s is %7.4f; moving atoms together.",
1320 a1, a2, sqrt(dist)));
1321 a1.setXYZ(xyzAv);
1322 a2.setXYZ(xyzAv);
1323 } else if (dist > 0) {
1324
1325 a1.setXYZ(xyzAv);
1326 a2.setXYZ(xyzAv);
1327 }
1328 }
1329
1330 private void packGradient(double[] x, double[] g) {
1331 if (g == null) {
1332 g = new double[nVariables];
1333 }
1334 int indexCommon = 0;
1335 int indexUnique = nShared * 3;
1336
1337
1338 int index = 0;
1339 for (int i = 0; i < nActive1; i++) {
1340 if (sharedAtoms1[i]) {
1341 g[indexCommon++] = f1L * g1[index] + f2L * rg1[index++];
1342 g[indexCommon++] = f1L * g1[index] + f2L * rg1[index++];
1343 g[indexCommon++] = f1L * g1[index] + f2L * rg1[index++];
1344 } else {
1345 g[indexUnique++] = f1L * g1[index] + f2L * rg1[index++];
1346 g[indexUnique++] = f1L * g1[index] + f2L * rg1[index++];
1347 g[indexUnique++] = f1L * g1[index] + f2L * rg1[index++];
1348 }
1349 }
1350
1351 if (!useFirstSystemBondedEnergy) {
1352
1353 indexCommon = 0;
1354 index = 0;
1355 for (int i = 0; i < nActive2; i++) {
1356 if (sharedAtoms2[i]) {
1357 g[indexCommon++] += f2L * g2[index] + f1L * rg2[index++];
1358 g[indexCommon++] += f2L * g2[index] + f1L * rg2[index++];
1359 g[indexCommon++] += f2L * g2[index] + f1L * rg2[index++];
1360 } else {
1361 g[indexUnique++] = f2L * g2[index] + f1L * rg2[index++];
1362 g[indexUnique++] = f2L * g2[index] + f1L * rg2[index++];
1363 g[indexUnique++] = f2L * g2[index] + f1L * rg2[index++];
1364 }
1365 }
1366 } else {
1367
1368 indexCommon = 0;
1369 index = 0;
1370 for (int i = 0; i < nActive2; i++) {
1371 if (sharedAtoms2[i]) {
1372 g[indexCommon++] += f2L * g2[index] - f2L * rg2[index++];
1373 g[indexCommon++] += f2L * g2[index] - f2L * rg2[index++];
1374 g[indexCommon++] += f2L * g2[index] - f2L * rg2[index++];
1375 } else {
1376 g[indexUnique++] = f2L * g2[index] - f2L * rg2[index++];
1377 g[indexUnique++] = f2L * g2[index] - f2L * rg2[index++];
1378 g[indexUnique++] = f2L * g2[index] - f2L * rg2[index++];
1379 }
1380 }
1381 }
1382
1383 scaleCoordinatesAndGradient(x, g);
1384 }
1385
1386 private void unpackCoordinates(double[] x) {
1387
1388
1389 unscaleCoordinates(x);
1390
1391 int index = 0;
1392 int indexCommon = 0;
1393 int indexUnique = 3 * nShared;
1394 for (int i = 0; i < nActive1; i++) {
1395 if (sharedAtoms1[i]) {
1396 x1[index++] = x[indexCommon++];
1397 x1[index++] = x[indexCommon++];
1398 x1[index++] = x[indexCommon++];
1399 } else {
1400 x1[index++] = x[indexUnique++];
1401 x1[index++] = x[indexUnique++];
1402 x1[index++] = x[indexUnique++];
1403 }
1404 }
1405
1406 index = 0;
1407 indexCommon = 0;
1408 for (int i = 0; i < nActive2; i++) {
1409 if (sharedAtoms2[i]) {
1410 x2[index++] = x[indexCommon++];
1411 x2[index++] = x[indexCommon++];
1412 x2[index++] = x[indexCommon++];
1413 } else {
1414 x2[index++] = x[indexUnique++];
1415 x2[index++] = x[indexUnique++];
1416 x2[index++] = x[indexUnique++];
1417 }
1418 }
1419 }
1420
1421
1422
1423
1424
1425
1426
1427
1428 void reloadCommonMasses(boolean secondTopology) {
1429 int commonIndex = 0;
1430 if (secondTopology) {
1431 for (int i = 0; i < nActive2; i++) {
1432 Atom a = activeAtoms2[i];
1433 double m = a.getMass();
1434 if (sharedAtoms2[i]) {
1435 mass[commonIndex++] = m;
1436 mass[commonIndex++] = m;
1437 mass[commonIndex++] = m;
1438 }
1439 }
1440 } else {
1441 for (int i = 0; i < nActive1; i++) {
1442 Atom a = activeAtoms1[i];
1443 double m = a.getMass();
1444 if (sharedAtoms1[i]) {
1445 mass[commonIndex++] = m;
1446 mass[commonIndex++] = m;
1447 mass[commonIndex++] = m;
1448 }
1449 }
1450 }
1451 }
1452
1453 private class EnergyRegion extends ParallelRegion {
1454
1455 private final Energy1Section e1sect;
1456 private final Energy2Section e2sect;
1457 private double[] x;
1458 private double[] g;
1459 private boolean gradient = false;
1460 private boolean verbose = false;
1461
1462 EnergyRegion() {
1463 e1sect = new Energy1Section();
1464 e2sect = new Energy2Section();
1465 }
1466
1467 @Override
1468 public void finish() {
1469
1470 if (!useFirstSystemBondedEnergy) {
1471 totalEnergy =
1472 f1L * energy1 + f2L * restraintEnergy1 + f2L * energy2 + f1L * restraintEnergy2;
1473 } else {
1474 totalEnergy =
1475 f1L * energy1 + f2L * restraintEnergy1 + f2L * energy2 - f2L * restraintEnergy2;
1476 }
1477
1478 if (gradient) {
1479 packGradient(x, g);
1480 } else {
1481 scaleCoordinates(x);
1482 }
1483
1484 if (verbose) {
1485 logger.info(format(" Total dual-topology energy: %12.4f", totalEnergy));
1486 }
1487 setVerbose(false);
1488 setGradient(false);
1489 }
1490
1491 @Override
1492 public void run() throws Exception {
1493 execute(e1sect, e2sect);
1494 }
1495
1496 public void setG(double[] g) {
1497 this.g = g;
1498 setGradient(true);
1499 }
1500
1501 public void setGradient(boolean grad) {
1502 this.gradient = grad;
1503 e1sect.setGradient(grad);
1504 e2sect.setGradient(grad);
1505 }
1506
1507 public void setVerbose(boolean verb) {
1508 this.verbose = verb;
1509 e1sect.setVerbose(verb);
1510 e2sect.setVerbose(verb);
1511 }
1512
1513 public void setX(double[] x) {
1514 this.x = x;
1515 }
1516
1517 @Override
1518 public void start() {
1519 unpackCoordinates(x);
1520 }
1521 }
1522
1523 private class Energy1Section extends ParallelSection {
1524
1525 private boolean gradient = false;
1526 private boolean verbose = false;
1527
1528 @Override
1529 public void run() {
1530 if (gradient) {
1531 fill(gl1, 0.0);
1532 fill(rgl1, 0.0);
1533 if (f1L > 0.0) {
1534 energy1 = potential1.energyAndGradient(x1, g1, verbose);
1535 dEdL_1 = lambdaInterface1.getdEdL();
1536 d2EdL2_1 = lambdaInterface1.getd2EdL2();
1537 lambdaInterface1.getdEdXdL(gl1);
1538 } else {
1539
1540 energy1 = 0.0;
1541 dEdL_1 = 0.0;
1542 d2EdL2_1 = 0.0;
1543 }
1544 if (doValenceRestraint1 && (f2L > 0.0) && potential1 instanceof ForceFieldEnergy) {
1545 forceFieldEnergy1.setLambdaBondedTerms(true, useFirstSystemBondedEnergy);
1546 if (verbose) {
1547 logger.info(" Calculating lambda bonded terms for topology 1");
1548 }
1549 restraintEnergy1 = forceFieldEnergy1.energyAndGradient(x1, rg1, verbose);
1550 restraintdEdL_1 = forceFieldEnergy1.getdEdL();
1551 restraintd2EdL2_1 = forceFieldEnergy1.getd2EdL2();
1552 forceFieldEnergy1.getdEdXdL(rgl1);
1553 forceFieldEnergy1.setLambdaBondedTerms(false, false);
1554 } else {
1555 restraintEnergy1 = 0.0;
1556 restraintdEdL_1 = 0.0;
1557 restraintd2EdL2_1 = 0.0;
1558 }
1559 if (logger.isLoggable(Level.FINE)) {
1560 logger.fine(format(" Topology 1 Energy & Restraints: %15.8f %15.8f\n", f1L * energy1, f2L * restraintEnergy1));
1561 logger.fine(format(" Topology 1: %15.8f * (%.2f)", energy1, f1L));
1562 logger.fine(format(" T1 Restraints: %15.8f * (%.2f)", restraintEnergy1, f2L));
1563 }
1564 } else {
1565 if (f1L > 0.0) {
1566 energy1 = potential1.energy(x1, verbose);
1567 } else {
1568 energy1 = 0.0;
1569 }
1570 if (doValenceRestraint1 && (f2L > 0.0) && potential1 instanceof ForceFieldEnergy) {
1571 forceFieldEnergy1.setLambdaBondedTerms(true, useFirstSystemBondedEnergy);
1572 if (verbose) {
1573 logger.info(" Calculating lambda bonded terms for topology 1");
1574 }
1575 restraintEnergy1 = potential1.energy(x1, verbose);
1576 forceFieldEnergy1.setLambdaBondedTerms(false, false);
1577 } else {
1578 restraintEnergy1 = 0.0;
1579 }
1580 if (logger.isLoggable(Level.FINE)) {
1581 logger.fine(format(" Topology 1 Energy & Restraints: %15.8f %15.8f\n", f1L * energy1, f2L * restraintEnergy1));
1582 }
1583 }
1584 }
1585
1586 public void setGradient(boolean grad) {
1587 this.gradient = grad;
1588 }
1589
1590 public void setVerbose(boolean verb) {
1591 this.verbose = verb;
1592 }
1593 }
1594
1595 private class Energy2Section extends ParallelSection {
1596
1597 private boolean gradient = false;
1598 private boolean verbose = false;
1599
1600 @Override
1601 public void run() {
1602 int numSymOps = 0;
1603 if (useSymOp) {
1604 numSymOps = symOpAtoms.size();
1605 for (int i = 0; i < numSymOps; i++) {
1606 List<Integer> symAtoms = symOpAtoms.get(i);
1607 boolean[] mask = new boolean[nActive2];
1608 for (int j = 0; j < nActive2; j++) {
1609 if (sharedAtoms2[j]) {
1610 if (symAtoms.contains(j)) {
1611 mask[j] = true;
1612 }
1613 }
1614 }
1615 applyCartesianSymOp(x2, x2, symOp[i], mask);
1616 }
1617 }
1618
1619 if (gradient) {
1620 fill(gl2, 0.0);
1621 fill(rgl2, 0.0);
1622
1623
1624 if (f2L > 0.0) {
1625 energy2 = potential2.energyAndGradient(x2, g2, verbose);
1626 dEdL_2 = -lambdaInterface2.getdEdL();
1627 d2EdL2_2 = lambdaInterface2.getd2EdL2();
1628 lambdaInterface2.getdEdXdL(gl2);
1629 if (useSymOp) {
1630
1631 for (int i = 0; i < numSymOps; i++) {
1632 List<Integer> symAtoms = symOpAtoms.get(i);
1633 boolean[] mask = new boolean[nActive2];
1634 for (int j = 0; j < nActive2; j++) {
1635 if (sharedAtoms2[j]) {
1636 if (symAtoms.contains(j)) {
1637 mask[j] = true;
1638 }
1639 }
1640 }
1641 applyCartesianSymRot(g2, g2, inverse[i], mask);
1642 applyCartesianSymRot(gl2, gl2, inverse[i], mask);
1643 }
1644 }
1645 } else {
1646
1647 energy2 = 0.0;
1648 dEdL_2 = 0.0;
1649 d2EdL2_2 = 0.0;
1650 }
1651
1652 if (doValenceRestraint2 && f1L > 0.0) {
1653 forceFieldEnergy2.setLambdaBondedTerms(true, useFirstSystemBondedEnergy);
1654 if (verbose) {
1655 logger.info(" Calculating lambda bonded terms for topology 2");
1656 }
1657 restraintEnergy2 = forceFieldEnergy2.energyAndGradient(x2, rg2, verbose);
1658 restraintdEdL_2 = -forceFieldEnergy2.getdEdL();
1659 restraintd2EdL2_2 = forceFieldEnergy2.getd2EdL2();
1660 forceFieldEnergy2.getdEdXdL(rgl2);
1661 if (useSymOp) {
1662
1663 for (int i = 0; i < numSymOps; i++) {
1664 List<Integer> symAtoms = symOpAtoms.get(i);
1665 boolean[] mask = new boolean[nActive2];
1666 for (int j = 0; j < nActive2; j++) {
1667 if (sharedAtoms2[j]) {
1668 if (symAtoms.contains(j)) {
1669 mask[j] = true;
1670 }
1671 }
1672 }
1673 applyCartesianSymRot(rg2, rg2, inverse[i], mask);
1674 applyCartesianSymRot(rgl2, rgl2, inverse[i], mask);
1675 }
1676 }
1677 forceFieldEnergy2.setLambdaBondedTerms(false, false);
1678 } else {
1679 restraintEnergy2 = 0.0;
1680 restraintdEdL_2 = 0.0;
1681 restraintd2EdL2_2 = 0.0;
1682 }
1683 if (logger.isLoggable(Level.FINE)) {
1684 logger.fine(format(" Topology 2 Energy & Restraints: %15.8f %15.8f", f2L * energy2, f1L * restraintEnergy2));
1685 logger.fine(format(" Topology 2: %15.8f * (%.2f)", energy2, f2L));
1686 logger.fine(format(" T2 Restraints: %15.8f * (%.2f)", restraintEnergy2, f1L));
1687 }
1688 } else {
1689 if (f2L > 0.0) {
1690 energy2 = potential2.energy(x2, verbose);
1691 } else {
1692 energy2 = 0.0;
1693 }
1694 if (doValenceRestraint2 && (f1L > 0.0) && potential2 instanceof ForceFieldEnergy) {
1695 forceFieldEnergy2.setLambdaBondedTerms(true, useFirstSystemBondedEnergy);
1696 if (verbose) {
1697 logger.info(" Calculating lambda bonded terms for topology 2");
1698 }
1699 restraintEnergy2 = potential2.energy(x2, verbose);
1700 forceFieldEnergy2.setLambdaBondedTerms(false, false);
1701 } else {
1702 restraintEnergy2 = 0.0;
1703 }
1704 if (logger.isLoggable(Level.FINE)) {
1705 logger.fine(format(" Topology 2 Energy & Restraints: %15.8f %15.8f", f2L * energy2, f1L * restraintEnergy2));
1706 }
1707 }
1708 }
1709
1710 public void setGradient(boolean grad) {
1711 this.gradient = grad;
1712 }
1713
1714 public void setVerbose(boolean verb) {
1715 this.verbose = verb;
1716 }
1717 }
1718 }