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.IntegerForLoop;
41 import edu.rit.pj.IntegerSchedule;
42 import edu.rit.pj.ParallelRegion;
43 import edu.rit.pj.ParallelTeam;
44 import edu.rit.pj.reduction.SharedDouble;
45 import ffx.crystal.Crystal;
46 import ffx.crystal.CrystalPotential;
47 import ffx.crystal.LatticeSystem;
48 import ffx.crystal.NCSCrystal;
49 import ffx.crystal.ReplicatesCrystal;
50 import ffx.crystal.SpaceGroup;
51 import ffx.crystal.SpaceGroupDefinitions;
52 import ffx.crystal.SymOp;
53 import ffx.numerics.Constraint;
54 import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
55 import ffx.numerics.atomic.AtomicDoubleArray3D;
56 import ffx.numerics.math.Double3;
57 import ffx.numerics.switching.ConstantSwitch;
58 import ffx.numerics.switching.UnivariateFunctionFactory;
59 import ffx.numerics.switching.UnivariateSwitchingFunction;
60 import ffx.potential.bonded.Angle;
61 import ffx.potential.bonded.AngleTorsion;
62 import ffx.potential.bonded.Atom;
63 import ffx.potential.bonded.Bond;
64 import ffx.potential.bonded.BondedTerm;
65 import ffx.potential.bonded.ImproperTorsion;
66 import ffx.potential.bonded.LambdaInterface;
67 import ffx.potential.bonded.MSNode;
68 import ffx.potential.bonded.MultiResidue;
69 import ffx.potential.bonded.OutOfPlaneBend;
70 import ffx.potential.bonded.PiOrbitalTorsion;
71 import ffx.potential.bonded.RelativeSolvation;
72 import ffx.potential.bonded.RelativeSolvation.SolvationLibrary;
73 import ffx.potential.bonded.Residue;
74 import ffx.potential.bonded.RestrainDistance;
75 import ffx.potential.bonded.RestraintTorsion;
76 import ffx.potential.bonded.StretchBend;
77 import ffx.potential.bonded.StretchTorsion;
78 import ffx.potential.bonded.Torsion;
79 import ffx.potential.bonded.TorsionTorsion;
80 import ffx.potential.bonded.UreyBradley;
81 import ffx.potential.constraint.CcmaConstraint;
82 import ffx.potential.constraint.SettleConstraint;
83 import ffx.potential.extended.ExtendedSystem;
84 import ffx.potential.nonbonded.COMRestraint;
85 import ffx.potential.nonbonded.GeneralizedKirkwood;
86 import ffx.potential.nonbonded.NCSRestraint;
87 import ffx.potential.nonbonded.ParticleMeshEwald;
88 import ffx.potential.nonbonded.RestrainGroups;
89 import ffx.potential.nonbonded.RestrainPosition;
90 import ffx.potential.nonbonded.VanDerWaals;
91 import ffx.potential.nonbonded.VanDerWaalsTornado;
92 import ffx.potential.openmm.OpenMMEnergy;
93 import ffx.potential.parameters.AngleType;
94 import ffx.potential.parameters.AngleType.AngleMode;
95 import ffx.potential.parameters.BondType;
96 import ffx.potential.parameters.ForceField;
97 import ffx.potential.parameters.ForceField.ELEC_FORM;
98 import ffx.potential.parameters.OutOfPlaneBendType;
99 import ffx.potential.parameters.TorsionType;
100 import ffx.potential.utils.ConvexHullOps;
101 import ffx.potential.utils.EnergyException;
102 import ffx.potential.utils.PotentialsFunctions;
103 import ffx.potential.utils.PotentialsUtils;
104 import ffx.utilities.FFXProperty;
105 import org.apache.commons.configuration2.CompositeConfiguration;
106
107 import javax.annotation.Nullable;
108 import java.io.File;
109 import java.time.LocalDateTime;
110 import java.time.format.DateTimeFormatter;
111 import java.util.ArrayList;
112 import java.util.Arrays;
113 import java.util.Collections;
114 import java.util.HashSet;
115 import java.util.List;
116 import java.util.Set;
117 import java.util.logging.Level;
118 import java.util.logging.Logger;
119 import java.util.stream.Collectors;
120 import java.util.stream.Stream;
121
122 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
123 import static ffx.potential.bonded.BondedTerm.removeNeuralNetworkTerms;
124 import static ffx.potential.nonbonded.RestrainPosition.parseRestrainPositions;
125 import static ffx.potential.nonbonded.VanDerWaalsForm.DEFAULT_VDW_CUTOFF;
126 import static ffx.potential.nonbonded.pme.EwaldParameters.DEFAULT_EWALD_CUTOFF;
127 import static ffx.potential.parameters.ForceField.toEnumForm;
128 import static ffx.potential.parsers.XYZFileFilter.isXYZ;
129 import static ffx.utilities.PropertyGroup.NonBondedCutoff;
130 import static ffx.utilities.PropertyGroup.PotentialFunctionSelection;
131 import static java.lang.Double.isInfinite;
132 import static java.lang.Double.isNaN;
133 import static java.lang.String.format;
134 import static java.lang.System.arraycopy;
135 import static java.util.Arrays.sort;
136 import static org.apache.commons.io.FilenameUtils.removeExtension;
137 import static org.apache.commons.math3.util.FastMath.PI;
138 import static org.apache.commons.math3.util.FastMath.max;
139 import static org.apache.commons.math3.util.FastMath.min;
140 import static org.apache.commons.math3.util.FastMath.sqrt;
141
142
143
144
145
146
147
148 public class ForceFieldEnergy implements CrystalPotential, LambdaInterface {
149
150
151
152
153 public static final double DEFAULT_CONSTRAINT_TOLERANCE = 1E-4;
154
155
156
157 private static final Logger logger = Logger.getLogger(ForceFieldEnergy.class.getName());
158
159
160
161 private static final double toSeconds = 1.0e-9;
162
163
164
165 protected final MolecularAssembly molecularAssembly;
166
167
168
169
170 public final double maxDebugGradient;
171
172
173
174 private final ParallelTeam parallelTeam;
175
176
177
178 private final NCSRestraint ncsRestraint;
179
180
181
182 private final COMRestraint comRestraint;
183
184
185
186 private final VanDerWaals vanderWaals;
187
188
189
190 private final ANIEnergy aniEnergy;
191
192 private final List<Constraint> constraints;
193
194 @FFXProperty(name = "vdw-cutoff", propertyGroup = NonBondedCutoff, defaultValue = "12.0", description = """
195 Sets the cutoff distance value in Angstroms for van der Waals potential energy interactions.
196 The energy for any pair of van der Waals sites beyond the cutoff distance will be set to zero.
197 Other properties can be used to define the smoothing scheme near the cutoff distance.
198 The default cutoff distance in the absence of the vdw-cutoff keyword is infinite for nonperiodic
199 systems and 12.0 for periodic systems.
200 """)
201 private double vdwCutoff;
202
203 @FFXProperty(name = "ewald-cutoff", propertyGroup = NonBondedCutoff, defaultValue = "7.0", description = """
204 Sets the value in Angstroms of the real-space distance cutoff for use during Ewald summation.
205 By default, in the absence of the ewald-cutoff property, a value of 7.0 is used.
206 """)
207 private double ewaldCutoff;
208
209 @FFXProperty(name = "gk-cutoff", propertyGroup = NonBondedCutoff, defaultValue = "12.0", description = """
210 Sets the value in Angstroms of the generalized Kirkwood distance cutoff for use
211 during implicit solvent simulations. By default, in the absence of the gk-cutoff property,
212 no cutoff is used under aperiodic boundary conditions and the vdw-cutoff is used under PBC.
213 """)
214 private double gkCutoff;
215
216 private static final double DEFAULT_LIST_BUFFER = 2.0;
217
218 @FFXProperty(name = "list-buffer", propertyGroup = NonBondedCutoff, defaultValue = "2.0", description = """
219 Sets the size of the neighbor list buffer in Angstroms for potential energy functions.
220 This value is added to the actual cutoff distance to determine which pairs will be kept on the neighbor list.
221 This buffer value is used for all potential function neighbor lists.
222 The default value in the absence of the list-buffer keyword is 2.0 Angstroms.
223 """)
224 private double listBuffer;
225
226
227
228
229 private final double cutOff2;
230
231
232
233 private final double cutoffPlusBuffer;
234
235
236
237 private final ParticleMeshEwald particleMeshEwald;
238
239
240
241 private final ELEC_FORM elecForm;
242
243
244
245 private final boolean nnTermOrig;
246
247
248
249 private final boolean bondTermOrig;
250
251
252
253 private final boolean angleTermOrig;
254
255
256
257 private final boolean stretchBendTermOrig;
258
259
260
261 private final boolean ureyBradleyTermOrig;
262
263
264
265 private final boolean outOfPlaneBendTermOrig;
266
267
268
269 private final boolean torsionTermOrig;
270
271
272
273 private final boolean stretchTorsionTermOrig;
274
275
276
277 private final boolean angleTorsionTermOrig;
278
279
280
281 private final boolean improperTorsionTermOrig;
282
283
284
285 private final boolean piOrbitalTorsionTermOrig;
286
287
288
289 private final boolean torsionTorsionTermOrig;
290
291
292
293 private final boolean restrainPositionTermOrig;
294
295
296
297 private final boolean restrainGroupTermOrig;
298
299
300
301 private final boolean restrainDistanceTermOrig;
302
303
304
305 private final boolean restrainTorsionTermOrig;
306
307
308
309
310 private final boolean vanderWaalsTermOrig;
311
312
313
314 private final boolean multipoleTermOrig;
315
316
317
318 private final boolean polarizationTermOrig;
319
320
321
322 private final boolean generalizedKirkwoodTermOrig;
323
324 private boolean esvTermOrig;
325
326
327
328 private final boolean rigidHydrogens;
329
330
331
332 private final boolean lambdaTorsions;
333
334
335
336 private final RelativeSolvation relativeSolvation;
337 private final boolean relativeSolvationTerm;
338 private final Platform platform = Platform.FFX;
339
340
341
342
343 @FFXProperty(name = "lambdaterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
344 defaultValue = "false", description = "Specifies use of the Lambda state variable.")
345 protected boolean lambdaTerm;
346
347
348
349 private double lambda = 1.0;
350
351
352
353 protected double[] optimizationScaling = null;
354
355
356
357 public boolean lambdaBondedTerms = false;
358
359
360
361 boolean lambdaAllBondedTerms = false;
362
363
364
365
366 boolean destroyed = false;
367
368
369
370
371 private STATE state = STATE.BOTH;
372
373
374
375 private Atom[] atoms;
376
377
378
379 private Bond[] bonds;
380
381
382
383 private Angle[] angles;
384
385
386
387 private StretchBend[] stretchBends;
388
389
390
391 private UreyBradley[] ureyBradleys;
392
393
394
395 private OutOfPlaneBend[] outOfPlaneBends;
396
397
398
399 private Torsion[] torsions;
400
401
402
403 private ImproperTorsion[] improperTorsions;
404
405
406
407 private PiOrbitalTorsion[] piOrbitalTorsions;
408
409
410
411 private TorsionTorsion[] torsionTorsions;
412
413
414
415 private StretchTorsion[] stretchTorsions;
416
417
418
419 private AngleTorsion[] angleTorsions;
420
421
422
423 private RestrainDistance[] restrainDistances;
424
425
426
427 private RestraintTorsion[] restrainTorsions;
428
429
430
431 private RestrainPosition[] restrainPositions;
432
433
434
435 private final RestrainGroups restrainGroups;
436
437
438
439
440 private int nAtoms;
441
442
443
444 private int nBonds;
445
446
447
448 private int nAngles;
449
450
451
452 private int nStretchBends;
453
454
455
456 private int nUreyBradleys;
457
458
459
460 private int nOutOfPlaneBends;
461
462
463
464 private int nTorsions;
465
466
467
468 private int nImproperTorsions;
469
470
471
472 private int nPiOrbitalTorsions;
473
474
475
476 private int nTorsionTorsions;
477
478
479
480 private int nAngleTorsions;
481
482
483
484 private int nStretchTorsions;
485
486
487
488 private int nRestrainDistances;
489
490
491
492 private int nRestrainTorsions;
493
494
495
496 private int nRestrainPositions;
497
498
499
500 private int nRestrainGroups;
501
502
503
504
505 private int nVanDerWaalInteractions;
506
507
508
509 private int nPermanentInteractions;
510
511
512
513 private int nGKInteractions;
514
515
516
517
518 private Crystal crystal;
519
520
521
522 private BondedRegion bondedRegion;
523
524
525
526
527 private boolean nnTerm;
528
529
530
531
532 @FFXProperty(name = "bondterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
533 defaultValue = "true", description = "Specifies use of the bond stretch potential.")
534 private boolean bondTerm;
535
536
537
538
539 @FFXProperty(name = "angleterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
540 defaultValue = "true", description = "Specifies use of the angle bend potential.")
541 private boolean angleTerm;
542
543
544
545
546 @FFXProperty(name = "strbndterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
547 defaultValue = "true", description = "Specifies use of the stretch-bend potential.")
548 private boolean stretchBendTerm;
549
550
551
552
553 @FFXProperty(name = "ureyterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
554 defaultValue = "true", description = "Specifies use of the Urey-Bradley potential.")
555 private boolean ureyBradleyTerm;
556
557
558
559
560 @FFXProperty(name = "opbendterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
561 defaultValue = "true", description = "Specifies use of the out-of-plane potential.")
562 private boolean outOfPlaneBendTerm;
563
564
565
566
567 @FFXProperty(name = "torsionterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
568 defaultValue = "true", description = "Specifies use of the torsional potential.")
569 private boolean torsionTerm;
570
571
572
573
574 @FFXProperty(name = "strtorterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
575 defaultValue = "true", description = "Specifies use of the stretch-torsion potential.")
576 private boolean stretchTorsionTerm;
577
578
579
580
581 @FFXProperty(name = "angtorterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
582 defaultValue = "true", description = "Specifies use of the angle-torsion potential.")
583 private boolean angleTorsionTerm;
584
585
586
587
588 @FFXProperty(name = "imptorterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
589 defaultValue = "true", description = "Specifies use of the improper torsion potential.")
590 private boolean improperTorsionTerm;
591
592
593
594
595 @FFXProperty(name = "pitorsterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
596 defaultValue = "true", description = "Specifies use of the pi-system torsion potential.")
597 private boolean piOrbitalTorsionTerm;
598
599
600
601
602 @FFXProperty(name = "tortorterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
603 defaultValue = "true", description = "Specifies use of the pi-system torsion potential.")
604 private boolean torsionTorsionTerm;
605
606
607
608
609 @FFXProperty(name = "vdwterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
610 defaultValue = "true", description = """
611 Specifies use of the vdw der Waals potential.
612 If set to false, all non-bonded terms are turned off.
613 """)
614 private boolean vanderWaalsTerm;
615
616
617
618
619 @FFXProperty(name = "mpoleterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
620 defaultValue = "true", description = """
621 Specifies use of the fixed charge electrostatic potential.
622 Setting mpoleterm to false also turns off polarization and generalized Kirkwood,
623 overriding the polarizeterm and gkterm properties.
624 """)
625 private boolean multipoleTerm;
626
627
628
629
630 @FFXProperty(name = "polarizeterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
631 defaultValue = "true", description = """
632 Specifies use of the polarizable electrostatic potential.
633 Setting polarizeterm to false overrides the polarization property.
634 """)
635 private boolean polarizationTerm;
636
637
638
639
640 @FFXProperty(name = "gkterm", clazz = Boolean.class, propertyGroup = PotentialFunctionSelection,
641 defaultValue = "false", description = "Specifies use of generalized Kirkwood electrostatics.")
642 private boolean generalizedKirkwoodTerm;
643
644
645
646
647 private boolean comTerm;
648
649
650
651 private boolean comTermOrig;
652
653
654
655 private boolean ncsTerm;
656
657
658
659 private boolean ncsTermOrig;
660
661
662
663 private boolean restrainDistanceTerm;
664
665
666
667 private boolean restrainTorsionTerm;
668
669
670
671 private boolean restrainPositionTerm;
672
673
674
675 private boolean restrainGroupTerm;
676
677
678
679
680 private double rigidScale;
681
682
683
684
685 private double nnEnergy;
686
687
688
689 private double bondEnergy;
690
691
692
693 private double bondRMSD;
694
695
696
697 private double angleEnergy;
698
699
700
701 private double angleRMSD;
702
703
704
705 private double stretchBendEnergy;
706
707
708
709 private double ureyBradleyEnergy;
710
711
712
713 private double outOfPlaneBendEnergy;
714
715
716
717 private double torsionEnergy;
718
719
720
721 private double stretchTorsionEnergy;
722
723
724
725 private double angleTorsionEnergy;
726
727
728
729 private double improperTorsionEnergy;
730
731
732
733 private double piOrbitalTorsionEnergy;
734
735
736
737 private double torsionTorsionEnergy;
738
739
740
741
742 private double restrainDistanceEnergy;
743
744
745
746 private double restrainPositionEnergy;
747
748
749
750 private double restrainGroupEnergy;
751
752
753
754 private double restrainTorsionEnergy;
755
756
757
758
759 private double totalBondedEnergy;
760
761
762
763 private double vanDerWaalsEnergy;
764
765
766
767 private double permanentMultipoleEnergy;
768
769
770
771 private double permanentRealSpaceEnergy;
772
773
774
775 private double polarizationEnergy;
776
777
778
779 private double totalMultipoleEnergy;
780
781
782
783 private double totalNonBondedEnergy;
784
785
786
787 private double solvationEnergy;
788
789
790
791 private double ncsEnergy;
792
793
794
795
796 private double comRestraintEnergy;
797
798
799
800 private double totalEnergy;
801
802
803
804
805 private long nnTime;
806
807
808
809 private long bondTime;
810
811
812
813 private long angleTime;
814
815
816
817 private long stretchBendTime;
818
819
820
821 private long ureyBradleyTime;
822
823
824
825 private long outOfPlaneBendTime;
826
827
828
829 private long torsionTime;
830
831
832
833 private long angleTorsionTime;
834
835
836
837 private long stretchTorsionTime;
838
839
840
841 private long piOrbitalTorsionTime;
842
843
844
845 private long improperTorsionTime;
846
847
848
849 private long torsionTorsionTime;
850
851
852
853 private long vanDerWaalsTime;
854
855
856
857 private long electrostaticTime;
858
859
860
861 private long restrainDistanceTime;
862
863
864
865 private long restrainTorsionTime;
866
867
868
869
870 private long comRestraintTime;
871
872
873
874 private long restrainGroupTime;
875
876
877
878 private long totalTime;
879
880
881
882 private double[] xyz;
883
884
885
886 private ExtendedSystem esvSystem = null;
887
888 private boolean esvTerm;
889 private double esvBias;
890
891
892
893 private long ncsTime;
894 private int nRelativeSolvations;
895
896
897
898 private long restrainPositionTime;
899
900 private double relativeSolvationEnergy;
901
902
903
904 private boolean printOnFailure;
905
906
907
908
909
910
911 protected ForceFieldEnergy(MolecularAssembly molecularAssembly) {
912 this(molecularAssembly, ParallelTeam.getDefaultThreadCount());
913 }
914
915
916
917
918
919
920
921 protected ForceFieldEnergy(MolecularAssembly molecularAssembly, int numThreads) {
922
923 this.molecularAssembly = molecularAssembly;
924 atoms = molecularAssembly.getAtomArray();
925 nAtoms = atoms.length;
926 xyz = new double[nAtoms * 3];
927
928
929 for (int i = 0; i < nAtoms; i++) {
930 int index = atoms[i].getXyzIndex() - 1;
931 assert (i == index);
932 }
933
934
935 int nThreads = min(nAtoms, numThreads);
936 parallelTeam = new ParallelTeam(nThreads);
937
938 ForceField forceField = molecularAssembly.getForceField();
939 CompositeConfiguration properties = forceField.getProperties();
940 String name = forceField.toString().toUpperCase();
941
942 logger.info(format(" Constructing Force Field %s", name));
943 logger.info(format("\n SMP threads: %10d", nThreads));
944
945
946 boolean standardizeAtomNames = properties.getBoolean("standardizeAtomNames", true);
947 if (standardizeAtomNames) {
948 molecularAssembly.renameWaterProtons();
949 }
950
951 nnTerm = forceField.getBoolean("NNTERM", false);
952
953 if (nnTerm) {
954 for (Atom atom : atoms) {
955 atom.setNeuralNetwork(true);
956 }
957 }
958 bondTerm = forceField.getBoolean("BONDTERM", true);
959 angleTerm = forceField.getBoolean("ANGLETERM", true);
960 stretchBendTerm = forceField.getBoolean("STRBNDTERM", true);
961 ureyBradleyTerm = forceField.getBoolean("UREYTERM", true);
962 outOfPlaneBendTerm = forceField.getBoolean("OPBENDTERM", true);
963 torsionTerm = forceField.getBoolean("TORSIONTERM", true);
964 stretchTorsionTerm = forceField.getBoolean("STRTORTERM", true);
965 angleTorsionTerm = forceField.getBoolean("ANGTORTERM", true);
966 piOrbitalTorsionTerm = forceField.getBoolean("PITORSTERM", true);
967 torsionTorsionTerm = forceField.getBoolean("TORTORTERM", true);
968 improperTorsionTerm = forceField.getBoolean("IMPROPERTERM", true);
969 vanderWaalsTerm = forceField.getBoolean("VDWTERM", true);
970
971 boolean tornadoVM = forceField.getBoolean("tornado", false);
972 if (vanderWaalsTerm && !tornadoVM) {
973 multipoleTerm = forceField.getBoolean("MPOLETERM", true);
974 if (multipoleTerm) {
975 String polarizeString = forceField.getString("POLARIZATION", "NONE");
976 boolean defaultPolarizeTerm = !polarizeString.equalsIgnoreCase("NONE");
977 polarizationTerm = forceField.getBoolean("POLARIZETERM", defaultPolarizeTerm);
978 generalizedKirkwoodTerm = forceField.getBoolean("GKTERM", false);
979 } else {
980
981 polarizationTerm = false;
982 generalizedKirkwoodTerm = false;
983 }
984 } else {
985
986 multipoleTerm = false;
987 polarizationTerm = false;
988 generalizedKirkwoodTerm = false;
989 }
990
991
992 lambdaTerm = forceField.getBoolean("LAMBDATERM", false);
993 comTerm = forceField.getBoolean("COMRESTRAINTERM", false);
994 lambdaTorsions = forceField.getBoolean("TORSION_LAMBDATERM", false);
995 printOnFailure = forceField.getBoolean("PRINT_ON_FAILURE", false);
996
997
998 if (properties.containsKey("restrain-distance")) {
999 restrainDistanceTerm = true;
1000 } else {
1001 restrainDistanceTerm = false;
1002 }
1003
1004
1005 if (properties.containsKey("restrain-torsion-cos")) {
1006 restrainTorsionTerm = true;
1007 } else {
1008 restrainTorsionTerm = false;
1009 }
1010
1011
1012 if (properties.containsKey("restrain-groups")) {
1013 restrainGroupTerm = true;
1014 } else {
1015 restrainGroupTerm = false;
1016 }
1017
1018
1019 if (properties.containsKey("restrain-position") || properties.containsKey("restrain-position-lambda")) {
1020 restrainPositionTerm = true;
1021 } else {
1022 restrainPositionTerm = false;
1023 }
1024
1025
1026
1027 nnTermOrig = nnTerm;
1028 bondTermOrig = bondTerm;
1029 angleTermOrig = angleTerm;
1030 stretchBendTermOrig = stretchBendTerm;
1031 ureyBradleyTermOrig = ureyBradleyTerm;
1032 outOfPlaneBendTermOrig = outOfPlaneBendTerm;
1033 torsionTermOrig = torsionTerm;
1034 angleTorsionTermOrig = angleTorsionTerm;
1035 stretchTorsionTermOrig = stretchTorsionTerm;
1036 improperTorsionTermOrig = improperTorsionTerm;
1037 piOrbitalTorsionTermOrig = piOrbitalTorsionTerm;
1038 torsionTorsionTermOrig = torsionTorsionTerm;
1039 vanderWaalsTermOrig = vanderWaalsTerm;
1040 multipoleTermOrig = multipoleTerm;
1041 polarizationTermOrig = polarizationTerm;
1042 generalizedKirkwoodTermOrig = generalizedKirkwoodTerm;
1043 ncsTermOrig = ncsTerm;
1044 comTermOrig = comTerm;
1045 restrainDistanceTermOrig = restrainDistanceTerm;
1046 restrainTorsionTermOrig = restrainTorsionTerm;
1047 restrainPositionTermOrig = restrainPositionTerm;
1048 restrainGroupTermOrig = restrainGroupTerm;
1049
1050
1051 String spacegroup;
1052 double a, b, c, alpha, beta, gamma;
1053 boolean aperiodic;
1054 try {
1055 spacegroup = forceField.getString("SPACEGROUP", "P 1");
1056 SpaceGroup sg = SpaceGroupDefinitions.spaceGroupFactory(spacegroup);
1057 LatticeSystem latticeSystem = sg.latticeSystem;
1058 a = forceField.getDouble("A_AXIS");
1059 aperiodic = false;
1060 b = forceField.getDouble("B_AXIS", latticeSystem.getDefaultBAxis(a));
1061 c = forceField.getDouble("C_AXIS", latticeSystem.getDefaultCAxis(a, b));
1062 alpha = forceField.getDouble("ALPHA", latticeSystem.getDefaultAlpha());
1063 beta = forceField.getDouble("BETA", latticeSystem.getDefaultBeta());
1064 gamma = forceField.getDouble("GAMMA", latticeSystem.getDefaultGamma());
1065 if (!sg.latticeSystem.validParameters(a, b, c, alpha, beta, gamma)) {
1066 logger.severe(" Check lattice parameters.");
1067 }
1068 if (a == 1.0 && b == 1.0 && c == 1.0) {
1069 String message = " A-, B-, and C-axis values equal to 1.0.";
1070 logger.info(message);
1071 throw new Exception(message);
1072 }
1073 if (a <= 0.0 || b <= 0.0 || c <= 0.0 || alpha <= 0.0 || beta <= 0.0 || gamma <= 0.0) {
1074
1075 String message = " Crystal parameters are not valid due to negative or zero value.";
1076 logger.warning(message);
1077 throw new Exception(message);
1078 }
1079 } catch (Exception e) {
1080 aperiodic = true;
1081
1082
1083 double maxR = 0.0;
1084 if (nAtoms < 10) {
1085 for (int i = 0; i < nAtoms - 1; i++) {
1086 Double3 xi = atoms[i].getXYZ();
1087 for (int j = 1; j < nAtoms; j++) {
1088 double r = atoms[j].getXYZ().dist(xi);
1089 maxR = max(r, maxR);
1090 }
1091 }
1092 } else {
1093 try {
1094 maxR = ConvexHullOps.maxDist(ConvexHullOps.constructHull(atoms));
1095 } catch (Exception ignored) {
1096
1097 for (int i = 0; i < nAtoms - 1; i++) {
1098 Double3 xi = atoms[i].getXYZ();
1099 for (int j = 1; j < nAtoms; j++) {
1100 double r = atoms[j].getXYZ().dist(xi);
1101 maxR = max(r, maxR);
1102 }
1103 }
1104 }
1105 }
1106 maxR = max(10.0, maxR);
1107
1108 logger.info(format(" The system will be treated as aperiodic (max separation = %6.1f A).", maxR));
1109
1110
1111 forceField.addProperty("EWALD_ALPHA", "0.0");
1112
1113
1114 spacegroup = "P1";
1115 a = 2.0 * maxR;
1116 b = 2.0 * maxR;
1117 c = 2.0 * maxR;
1118 alpha = 90.0;
1119 beta = 90.0;
1120 gamma = 90.0;
1121 }
1122 Crystal unitCell = new Crystal(a, b, c, alpha, beta, gamma, spacegroup);
1123 unitCell.setAperiodic(aperiodic);
1124
1125
1126 vdwCutoff = DEFAULT_VDW_CUTOFF;
1127 ewaldCutoff = DEFAULT_EWALD_CUTOFF;
1128 gkCutoff = vdwCutoff;
1129 if (unitCell.aperiodic()) {
1130
1131 vdwCutoff = Double.POSITIVE_INFINITY;
1132 ewaldCutoff = Double.POSITIVE_INFINITY;
1133 gkCutoff = Double.POSITIVE_INFINITY;
1134 }
1135
1136 vdwCutoff = forceField.getDouble("VDW_CUTOFF", vdwCutoff);
1137 ewaldCutoff = forceField.getDouble("EWALD_CUTOFF", ewaldCutoff);
1138 gkCutoff = forceField.getDouble("GK_CUTOFF", gkCutoff);
1139
1140
1141 double nlistCutoff = vdwCutoff;
1142 if (multipoleTerm) {
1143 nlistCutoff = max(vdwCutoff, ewaldCutoff);
1144 }
1145 if (generalizedKirkwoodTerm) {
1146 nlistCutoff = max(nlistCutoff, gkCutoff);
1147 }
1148
1149
1150 boolean disabledNeighborUpdates = forceField.getBoolean("DISABLE_NEIGHBOR_UPDATES", false);
1151 if (disabledNeighborUpdates) {
1152 logger.info(format(" Neighbor list updates disabled; interactions will "
1153 + "only be calculated between atoms that started the simulation "
1154 + "within a radius of %9.3g Angstroms of each other.", nlistCutoff));
1155
1156 vdwCutoff = Double.POSITIVE_INFINITY;
1157 ewaldCutoff = Double.POSITIVE_INFINITY;
1158 gkCutoff = Double.POSITIVE_INFINITY;
1159 }
1160
1161 listBuffer = forceField.getDouble("LIST_BUFFER", DEFAULT_LIST_BUFFER);
1162 cutoffPlusBuffer = nlistCutoff + listBuffer;
1163 cutOff2 = 2.0 * cutoffPlusBuffer;
1164 unitCell = configureNCS(forceField, unitCell);
1165
1166
1167 if (!aperiodic) {
1168 this.crystal = ReplicatesCrystal.replicatesCrystalFactory(unitCell, cutOff2);
1169 logger.info(format("\n Density: %6.3f (g/cc)",
1170 crystal.getDensity(molecularAssembly.getMass())));
1171 logger.info(crystal.toString());
1172 } else {
1173 this.crystal = unitCell;
1174 }
1175
1176 if (!unitCell.aperiodic() && unitCell.spaceGroup.number == 1) {
1177 ncsTerm = forceField.getBoolean("NCSTERM", false);
1178 ncsTermOrig = ncsTerm;
1179 } else {
1180 ncsTerm = false;
1181 ncsTermOrig = false;
1182 }
1183
1184
1185 int nSpecial = checkForSpecialPositions(forceField);
1186
1187 rigidHydrogens = forceField.getBoolean("RIGID_HYDROGENS", false);
1188 rigidScale = forceField.getDouble("RIGID_SCALE", 10.0);
1189
1190 nRelativeSolvations = 0;
1191 String relSolvLibrary = forceField.getString("RELATIVE_SOLVATION", "NONE").toUpperCase();
1192 SolvationLibrary library = SolvationLibrary.valueOf(relSolvLibrary);
1193 if (library == SolvationLibrary.AUTO) {
1194 if (generalizedKirkwoodTerm && name.toUpperCase().contains("OPLS")) {
1195 library = SolvationLibrary.MACCALLUM_TIP4P;
1196
1197 } else if (generalizedKirkwoodTerm) {
1198 library = SolvationLibrary.GK;
1199 } else {
1200 library = SolvationLibrary.NONE;
1201 }
1202 }
1203 if (library != SolvationLibrary.NONE) {
1204 relativeSolvationTerm = true;
1205 relativeSolvation = new RelativeSolvation(library, forceField);
1206 } else {
1207 relativeSolvationTerm = false;
1208 relativeSolvation = null;
1209 }
1210
1211 boolean checkAllNodeCharges = forceField.getBoolean("CHECK_ALL_NODE_CHARGES", false);
1212
1213 if (rigidScale <= 1.0) {
1214 rigidScale = 1.0;
1215 }
1216
1217 logger.info("\n Bonded Terms");
1218 if (rigidHydrogens && rigidScale > 1.0) {
1219 logger.info(format(" Rigid hydrogens: %10.2f", rigidScale));
1220 }
1221
1222
1223 if (bondTerm) {
1224 List<Bond> bondList = molecularAssembly.getBondList();
1225 if (nnTerm) {
1226 removeNeuralNetworkTerms(bondList);
1227 }
1228 nBonds = bondList.size();
1229 bonds = bondList.toArray(new Bond[0]);
1230 sort(bonds);
1231 if (nBonds > 0) {
1232 logger.info(format(" Bonds: %10d", nBonds));
1233 }
1234 } else {
1235 nBonds = 0;
1236 bonds = null;
1237 }
1238
1239
1240 if (angleTerm) {
1241 List<Angle> angleList = molecularAssembly.getAngleList();
1242 if (nnTerm) {
1243 removeNeuralNetworkTerms(angleList);
1244 }
1245 nAngles = angleList.size();
1246 angles = angleList.toArray(new Angle[0]);
1247 sort(angles);
1248 if (nAngles > 0) {
1249 logger.info(format(" Angles: %10d", nAngles));
1250 }
1251 } else {
1252 nAngles = 0;
1253 angles = null;
1254 }
1255
1256
1257 if (stretchBendTerm) {
1258 List<StretchBend> stretchBendList = molecularAssembly.getStretchBendList();
1259 if (nnTerm) {
1260 removeNeuralNetworkTerms(stretchBendList);
1261 }
1262 nStretchBends = stretchBendList.size();
1263 stretchBends = stretchBendList.toArray(new StretchBend[0]);
1264 sort(stretchBends);
1265 if (nStretchBends > 0) {
1266 logger.info(format(" Stretch-Bends: %10d", nStretchBends));
1267 }
1268 } else {
1269 nStretchBends = 0;
1270 stretchBends = null;
1271 }
1272
1273
1274 if (ureyBradleyTerm) {
1275 List<UreyBradley> ureyBradleyList = molecularAssembly.getUreyBradleyList();
1276 if (nnTerm) {
1277 removeNeuralNetworkTerms(ureyBradleyList);
1278 }
1279 nUreyBradleys = ureyBradleyList.size();
1280 ureyBradleys = ureyBradleyList.toArray(new UreyBradley[0]);
1281 sort(ureyBradleys);
1282 if (nUreyBradleys > 0) {
1283 logger.info(format(" Urey-Bradleys: %10d", nUreyBradleys));
1284 }
1285 } else {
1286 nUreyBradleys = 0;
1287 ureyBradleys = null;
1288 }
1289
1290
1291 if (rigidHydrogens) {
1292 if (bonds != null) {
1293 for (Bond bond : bonds) {
1294 if (bond.containsHydrogen()) {
1295 bond.setRigidScale(rigidScale);
1296 }
1297 }
1298 }
1299 if (angles != null) {
1300 for (Angle angle : angles) {
1301 if (angle.containsHydrogen()) {
1302 angle.setRigidScale(rigidScale);
1303 }
1304 }
1305 }
1306 if (stretchBends != null) {
1307 for (StretchBend stretchBend : stretchBends) {
1308 if (stretchBend.containsHydrogen()) {
1309 stretchBend.setRigidScale(rigidScale);
1310 }
1311 }
1312 }
1313 if (ureyBradleys != null) {
1314 for (UreyBradley ureyBradley : ureyBradleys) {
1315 if (ureyBradley.containsHydrogen()) {
1316 ureyBradley.setRigidScale(rigidScale);
1317 }
1318 }
1319 }
1320 }
1321
1322
1323 if (outOfPlaneBendTerm) {
1324 List<OutOfPlaneBend> outOfPlaneBendList = molecularAssembly.getOutOfPlaneBendList();
1325 if (nnTerm) {
1326 removeNeuralNetworkTerms(outOfPlaneBendList);
1327 }
1328 nOutOfPlaneBends = outOfPlaneBendList.size();
1329 outOfPlaneBends = outOfPlaneBendList.toArray(new OutOfPlaneBend[0]);
1330 sort(outOfPlaneBends);
1331 if (nOutOfPlaneBends > 0) {
1332 logger.info(format(" Out-of-Plane Bends: %10d", nOutOfPlaneBends));
1333 }
1334 } else {
1335 nOutOfPlaneBends = 0;
1336 outOfPlaneBends = null;
1337 }
1338
1339
1340 if (torsionTerm) {
1341 List<Torsion> torsionList = molecularAssembly.getTorsionList();
1342 if (nnTerm) {
1343 removeNeuralNetworkTerms(torsionList);
1344 }
1345 nTorsions = torsionList.size();
1346 torsions = torsionList.toArray(new Torsion[0]);
1347 if (nTorsions > 0) {
1348 logger.info(format(" Torsions: %10d", nTorsions));
1349 }
1350 } else {
1351 nTorsions = 0;
1352 torsions = null;
1353 }
1354
1355
1356 if (stretchTorsionTerm) {
1357 List<StretchTorsion> stretchTorsionList = molecularAssembly.getStretchTorsionList();
1358 if (nnTerm) {
1359 removeNeuralNetworkTerms(stretchTorsionList);
1360 }
1361 nStretchTorsions = stretchTorsionList.size();
1362 stretchTorsions = stretchTorsionList.toArray(new StretchTorsion[0]);
1363 if (nStretchTorsions > 0) {
1364 logger.info(format(" Stretch-Torsions: %10d", nStretchTorsions));
1365 }
1366 } else {
1367 nStretchTorsions = 0;
1368 stretchTorsions = null;
1369 }
1370
1371
1372 if (angleTorsionTerm) {
1373 List<AngleTorsion> angleTorsionList = molecularAssembly.getAngleTorsionList();
1374 if (nnTerm) {
1375 removeNeuralNetworkTerms(angleTorsionList);
1376 }
1377 nAngleTorsions = angleTorsionList.size();
1378 angleTorsions = angleTorsionList.toArray(new AngleTorsion[0]);
1379 if (nAngleTorsions > 0) {
1380 logger.info(format(" Angle-Torsions: %10d", nAngleTorsions));
1381 }
1382 } else {
1383 nAngleTorsions = 0;
1384 angleTorsions = null;
1385 }
1386
1387
1388 if (piOrbitalTorsionTerm) {
1389 List<PiOrbitalTorsion> piOrbitalTorsionList = molecularAssembly.getPiOrbitalTorsionList();
1390 if (nnTerm) {
1391 removeNeuralNetworkTerms(piOrbitalTorsionList);
1392 }
1393 nPiOrbitalTorsions = piOrbitalTorsionList.size();
1394 piOrbitalTorsions = piOrbitalTorsionList.toArray(new PiOrbitalTorsion[0]);
1395 if (nPiOrbitalTorsions > 0) {
1396 logger.info(format(" Pi-Orbital Torsions: %10d", nPiOrbitalTorsions));
1397 }
1398 } else {
1399 nPiOrbitalTorsions = 0;
1400 piOrbitalTorsions = null;
1401 }
1402
1403
1404 if (torsionTorsionTerm) {
1405 List<TorsionTorsion> torsionTorsionList = molecularAssembly.getTorsionTorsionList();
1406 if (nnTerm) {
1407 removeNeuralNetworkTerms(torsionTorsionList);
1408 }
1409 nTorsionTorsions = torsionTorsionList.size();
1410 torsionTorsions = torsionTorsionList.toArray(new TorsionTorsion[0]);
1411 if (nTorsionTorsions > 0) {
1412 logger.info(format(" Torsion-Torsions: %10d", nTorsionTorsions));
1413 }
1414 } else {
1415 nTorsionTorsions = 0;
1416 torsionTorsions = null;
1417 }
1418
1419
1420 if (improperTorsionTerm) {
1421 List<ImproperTorsion> improperTorsionList = molecularAssembly.getImproperTorsionList();
1422 if (nnTerm) {
1423 removeNeuralNetworkTerms(improperTorsionList);
1424 }
1425 nImproperTorsions = improperTorsionList.size();
1426 improperTorsions = improperTorsionList.toArray(new ImproperTorsion[0]);
1427 if (nImproperTorsions > 0) {
1428 logger.info(format(" Improper Torsions: %10d", nImproperTorsions));
1429 }
1430 } else {
1431 nImproperTorsions = 0;
1432 improperTorsions = null;
1433 }
1434
1435 int[] molecule = molecularAssembly.getMoleculeNumbers();
1436 if (vanderWaalsTerm) {
1437 logger.info("\n Non-Bonded Terms");
1438 boolean[] nn = null;
1439 if (nnTerm) {
1440 nn = molecularAssembly.getNeuralNetworkIdentity();
1441 } else {
1442 nn = new boolean[nAtoms];
1443 Arrays.fill(nn, false);
1444 }
1445
1446 if (!tornadoVM) {
1447 vanderWaals = new VanDerWaals(atoms, molecule, nn, crystal, forceField, parallelTeam,
1448 vdwCutoff, nlistCutoff);
1449 } else {
1450 vanderWaals = new VanDerWaalsTornado(atoms, crystal, forceField, vdwCutoff);
1451 }
1452 } else {
1453 vanderWaals = null;
1454 }
1455
1456 if (nnTerm) {
1457 aniEnergy = new ANIEnergy(molecularAssembly);
1458 } else {
1459 aniEnergy = null;
1460 }
1461
1462 if (multipoleTerm) {
1463 if (name.contains("OPLS") || name.contains("AMBER") || name.contains("CHARMM")) {
1464 elecForm = ELEC_FORM.FIXED_CHARGE;
1465 } else {
1466 elecForm = ELEC_FORM.PAM;
1467 }
1468 particleMeshEwald = new ParticleMeshEwald(atoms, molecule, forceField, crystal,
1469 vanderWaals.getNeighborList(), elecForm, ewaldCutoff, gkCutoff, parallelTeam);
1470 double charge = molecularAssembly.getCharge(checkAllNodeCharges);
1471 logger.info(format("\n Overall system charge: %10.3f", charge));
1472 } else {
1473 elecForm = null;
1474 particleMeshEwald = null;
1475 }
1476
1477 if (ncsTerm) {
1478 String sg = forceField.getString("NCSGROUP", "P 1");
1479 Crystal ncsCrystal = new Crystal(a, b, c, alpha, beta, gamma, sg);
1480 ncsRestraint = new NCSRestraint(atoms, forceField, ncsCrystal);
1481 } else {
1482 ncsRestraint = null;
1483 }
1484
1485 if (restrainPositionTerm) {
1486 restrainPositions = parseRestrainPositions(molecularAssembly);
1487 if (restrainPositions != null) {
1488 nRestrainPositions = restrainPositions.length;
1489 } else {
1490 nRestrainPositions = 0;
1491 }
1492 } else {
1493 restrainPositions = null;
1494 }
1495
1496 if (restrainGroupTerm) {
1497 restrainGroups = new RestrainGroups(molecularAssembly);
1498 nRestrainGroups = restrainGroups.getNumberOfRestraints();
1499 } else {
1500 restrainGroups = null;
1501 }
1502
1503 if (comTerm) {
1504 comRestraint = new COMRestraint(molecularAssembly);
1505 } else {
1506 comRestraint = null;
1507 }
1508
1509 if (restrainDistanceTerm) {
1510
1511 configureRestraintBonds(properties);
1512 }
1513
1514 if (restrainTorsionTerm) {
1515 nRestrainTorsions = configureRestraintTorsions(properties, forceField);
1516 }
1517
1518
1519 bondedRegion = new BondedRegion();
1520
1521 maxDebugGradient = forceField.getDouble("MAX_DEBUG_GRADIENT", Double.POSITIVE_INFINITY);
1522
1523 molecularAssembly.setPotential(this);
1524
1525
1526 constraints = configureConstraints(forceField);
1527
1528 if (lambdaTerm) {
1529 this.setLambda(1.0);
1530 if (nSpecial == 0) {
1531
1532
1533 crystal.setSpecialPositionCutoff(0.0);
1534 } else {
1535
1536 logger.info(" Special positions checking will be performed during a lambda simulation.");
1537 }
1538 }
1539 }
1540
1541 private int checkForSpecialPositions(ForceField forceField) {
1542
1543 boolean specialPositionsInactive = forceField.getBoolean("SPECIAL_POSITIONS_INACTIVE", true);
1544 int nSpecial = 0;
1545 if (specialPositionsInactive) {
1546 int nSymm = crystal.getNumSymOps();
1547 if (nSymm > 1) {
1548 SpaceGroup spaceGroup = crystal.spaceGroup;
1549 double sp2 = crystal.getSpecialPositionCutoff2();
1550 double[] mate = new double[3];
1551 StringBuilder sb = new StringBuilder("\n Atoms at Special Positions set to Inactive:\n");
1552 for (int i = 0; i < nAtoms; i++) {
1553 Atom atom = atoms[i];
1554 double[] xyz = atom.getXYZ(null);
1555 for (int iSymm = 1; iSymm < nSymm; iSymm++) {
1556 SymOp symOp = spaceGroup.getSymOp(iSymm);
1557 crystal.applySymOp(xyz, mate, symOp);
1558 double dr2 = crystal.image(xyz, mate);
1559 if (dr2 < sp2) {
1560 sb.append(
1561 format(" %s separation with SymOp %d at %8.6f A.\n", atoms[i], iSymm, sqrt(dr2)));
1562 atom.setActive(false);
1563 nSpecial++;
1564 break;
1565 }
1566 }
1567 }
1568 if (nSpecial > 0) {
1569 logger.info(sb.toString());
1570 }
1571 }
1572 }
1573 return nSpecial;
1574 }
1575
1576
1577
1578
1579
1580
1581
1582
1583 private int configureRestraintTorsions(CompositeConfiguration properties, ForceField forceField) {
1584 String[] restrainTorsionCos = properties.getStringArray("restrain-torsion-cos");
1585 double torsUnits = forceField.getDouble("TORSIONUNIT", TorsionType.DEFAULT_TORSION_UNIT);
1586 List<RestraintTorsion> rTorsList = new ArrayList<>(restrainTorsionCos.length);
1587 for (String rtc : restrainTorsionCos) {
1588 String[] toks = rtc.split("\\s+");
1589 int nTok = toks.length;
1590 if (nTok < 7) {
1591 throw new IllegalArgumentException(format(
1592 "restrain-torsion-cos record %s must have 4 atom indices, amplitude, phase shift and periodicity!",
1593 rtc));
1594 }
1595 int ai1 = Integer.parseInt(toks[0]) - 1;
1596 int ai2 = Integer.parseInt(toks[1]) - 1;
1597 int ai3 = Integer.parseInt(toks[2]) - 1;
1598 int ai4 = Integer.parseInt(toks[3]) - 1;
1599 Atom a1 = atoms[ai1];
1600 Atom a2 = atoms[ai2];
1601 Atom a3 = atoms[ai3];
1602 Atom a4 = atoms[ai4];
1603 int[] atomClasses = new int[]{-1, -1, -1, -1};
1604
1605 int startTerms = 4;
1606 boolean lamEnabled = false;
1607 boolean revLam = false;
1608 if (toks[4].equalsIgnoreCase("lambda")) {
1609 ++startTerms;
1610 lamEnabled = true;
1611 if (toks[5].toLowerCase().startsWith("reverse")) {
1612 ++startTerms;
1613 revLam = true;
1614 }
1615 }
1616 if (!lambdaTerm) {
1617 lamEnabled = false;
1618 }
1619
1620 int nTerms = (nTok - startTerms) / 3;
1621 assert (nTok - startTerms) % 3 == 0;
1622 double[] amp = new double[nTerms];
1623 double[] phase = new double[nTerms];
1624 int[] period = new int[nTerms];
1625
1626 for (int i = 0; i < nTerms; i++) {
1627 int i0 = startTerms + (3 * i);
1628 amp[i] = Double.parseDouble(toks[i0]);
1629 phase[i] = Double.parseDouble(toks[i0 + 1]);
1630 period[i] = Integer.parseInt(toks[i0 + 2]);
1631 }
1632
1633
1634 assert !lamEnabled || lambdaTerm;
1635 TorsionType tType = new TorsionType(atomClasses, amp, phase, period);
1636 rTorsList.add(new RestraintTorsion(a1, a2, a3, a4, tType, lamEnabled, revLam, torsUnits));
1637 }
1638
1639 if (!rTorsList.isEmpty()) {
1640 logger.info(format(" Adding %4d cosine-based torsion restraints.", nRestrainTorsions));
1641 restrainTorsionTerm = true;
1642 restrainTorsions = rTorsList.toArray(new RestraintTorsion[0]);
1643 restrainTorsionEnergy = 0;
1644 restrainTorsionTime = 0;
1645 return rTorsList.size();
1646 }
1647 return 0;
1648 }
1649
1650
1651
1652
1653
1654
1655 private void configureRestraintBonds(CompositeConfiguration properties) {
1656
1657 String[] bondRestraints = properties.getStringArray("restrain-distance");
1658 for (String bondRest : bondRestraints) {
1659 try {
1660 String[] toks = bondRest.split("\\s+");
1661 if (toks.length < 2) {
1662 throw new IllegalArgumentException(
1663 format(" restrain-distance value %s could not be parsed!", bondRest));
1664 }
1665
1666
1667 int at1 = Integer.parseInt(toks[0]) - 1;
1668 int at2 = Integer.parseInt(toks[1]) - 1;
1669
1670 double forceConst = 100.0;
1671 double flatBottomRadius = 0;
1672 Atom a1 = atoms[at1];
1673 Atom a2 = atoms[at2];
1674
1675 if (toks.length > 2) {
1676 forceConst = Double.parseDouble(toks[2]);
1677 }
1678 double dist;
1679 switch (toks.length) {
1680 case 2:
1681 case 3:
1682 double[] xyz1 = new double[3];
1683 xyz1 = a1.getXYZ(xyz1);
1684 double[] xyz2 = new double[3];
1685 xyz2 = a2.getXYZ(xyz2);
1686
1687 dist = crystal.minDistOverSymOps(xyz1, xyz2);
1688 break;
1689 case 4:
1690 dist = Double.parseDouble(toks[3]);
1691 break;
1692 case 5:
1693 default:
1694 double minDist = Double.parseDouble(toks[3]);
1695 double maxDist = Double.parseDouble(toks[4]);
1696 dist = 0.5 * (minDist + maxDist);
1697 flatBottomRadius = 0.5 * Math.abs(maxDist - minDist);
1698 break;
1699 }
1700
1701 UnivariateSwitchingFunction switchF;
1702 double lamStart = RestrainDistance.DEFAULT_RB_LAM_START;
1703 double lamEnd = RestrainDistance.DEFAULT_RB_LAM_END;
1704 if (toks.length > 5) {
1705 int offset = 5;
1706 if (toks[5].matches("^[01](?:\\.[0-9]*)?")) {
1707 offset = 6;
1708 lamStart = Double.parseDouble(toks[5]);
1709 if (toks[6].matches("^[01](?:\\.[0-9]*)?")) {
1710 offset = 7;
1711 lamEnd = Double.parseDouble(toks[6]);
1712 }
1713 }
1714 switchF = UnivariateFunctionFactory.parseUSF(toks, offset);
1715 } else {
1716 switchF = new ConstantSwitch();
1717 }
1718
1719 setRestrainDistance(a1, a2, dist, forceConst, flatBottomRadius, lamStart, lamEnd, switchF);
1720 } catch (Exception ex) {
1721 logger.info(format(" Exception in parsing restrain-distance: %s", ex));
1722 }
1723 }
1724 }
1725
1726
1727
1728
1729
1730
1731
1732 private List<Constraint> configureConstraints(ForceField forceField) {
1733 String constraintStrings = forceField.getString("CONSTRAIN", forceField.getString("RATTLE", null));
1734 if (constraintStrings == null) {
1735 return Collections.emptyList();
1736 }
1737
1738 ArrayList<Constraint> constraints = new ArrayList<>();
1739 logger.info(format(" Experimental: parsing constraints option %s", constraintStrings));
1740
1741 Set<Bond> numericBonds = new HashSet<>(1);
1742 Set<Angle> numericAngles = new HashSet<>(1);
1743
1744
1745 if (constraintStrings.isEmpty() || constraintStrings.matches("^\\s*$")) {
1746
1747 logger.info(" Constraining X-H bonds.");
1748 numericBonds = Arrays.stream(bonds).filter(
1749 (Bond bond) -> bond.getAtom(0).getAtomicNumber() == 1
1750 || bond.getAtom(1).getAtomicNumber() == 1).collect(Collectors.toSet());
1751 } else {
1752 String[] constraintToks = constraintStrings.split("\\s+");
1753
1754
1755 for (String tok : constraintToks) {
1756 if (tok.equalsIgnoreCase("WATER")) {
1757 logger.info(" Constraining waters to be rigid based on angle & bonds.");
1758
1759
1760 Stream<MSNode> settleStream = molecularAssembly.getMolecules().stream()
1761 .filter((MSNode m) -> m.getAtomList().size() == 3).filter((MSNode m) -> {
1762 List<Atom> atoms = m.getAtomList();
1763 Atom O = null;
1764 List<Atom> H = new ArrayList<>(2);
1765 for (Atom at : atoms) {
1766 int atN = at.getAtomicNumber();
1767 if (atN == 8) {
1768 O = at;
1769 } else if (atN == 1) {
1770 H.add(at);
1771 }
1772 }
1773 return O != null && H.size() == 2;
1774 });
1775
1776 settleStream = Stream.concat(settleStream, molecularAssembly.getWater().stream());
1777
1778 List<SettleConstraint> settleConstraints = settleStream.map(
1779 (MSNode m) -> m.getAngleList().getFirst()).map(SettleConstraint::settleFactory).toList();
1780 constraints.addAll(settleConstraints);
1781
1782 } else if (tok.equalsIgnoreCase("DIATOMIC")) {
1783 logger.severe(" Diatomic distance constraints not yet implemented properly.");
1784 } else if (tok.equalsIgnoreCase("TRIATOMIC")) {
1785 logger.severe(
1786 " Triatomic SETTLE constraints for non-water molecules not yet implemented properly.");
1787 }
1788 }
1789
1790
1791 for (String tok : constraintToks) {
1792 if (tok.equalsIgnoreCase("BONDS")) {
1793 numericBonds = new HashSet<>(Arrays.asList(bonds));
1794 } else if (tok.equalsIgnoreCase("ANGLES")) {
1795 numericAngles = new HashSet<>(Arrays.asList(angles));
1796 }
1797 }
1798 }
1799
1800
1801 for (Angle angle : numericAngles) {
1802 angle.getBondList().forEach(numericBonds::remove);
1803 }
1804
1805
1806 List<Angle> ccmaAngles = numericAngles.stream().filter((Angle ang) -> !ang.isConstrained())
1807 .collect(Collectors.toList());
1808 List<Bond> ccmaBonds = numericBonds.stream().filter((Bond bond) -> !bond.isConstrained())
1809 .collect(Collectors.toList());
1810
1811 CcmaConstraint ccmaConstraint = CcmaConstraint.ccmaFactory(ccmaBonds, ccmaAngles, atoms,
1812 getMass(), CcmaConstraint.DEFAULT_CCMA_NONZERO_CUTOFF);
1813 constraints.add(ccmaConstraint);
1814 logger.info(format(" Added %d constraints.", constraints.size()));
1815
1816 return constraints;
1817 }
1818
1819
1820
1821
1822
1823
1824
1825 public static ForceFieldEnergy energyFactory(MolecularAssembly assembly) {
1826 return energyFactory(assembly, ParallelTeam.getDefaultThreadCount());
1827 }
1828
1829
1830
1831
1832
1833
1834
1835
1836 @SuppressWarnings("fallthrough")
1837 public static ForceFieldEnergy energyFactory(MolecularAssembly assembly, int numThreads) {
1838 ForceField forceField = assembly.getForceField();
1839 String platformString = toEnumForm(forceField.getString("PLATFORM", "FFX"));
1840 try {
1841 Platform platform = Platform.valueOf(platformString);
1842 switch (platform) {
1843 case OMM, OMM_REF, OMM_CUDA:
1844 try {
1845 return new OpenMMEnergy(assembly, platform, numThreads);
1846 } catch (Exception ex) {
1847 logger.warning(format(" Exception creating OpenMMEnergy: %s", ex));
1848
1849 ForceFieldEnergy ffxEnergy = assembly.getPotentialEnergy();
1850 if (ffxEnergy == null) {
1851 ffxEnergy = new ForceFieldEnergy(assembly, numThreads);
1852 assembly.setPotential(ffxEnergy);
1853 }
1854 return ffxEnergy;
1855 }
1856 case OMM_OPENCL, OMM_OPTCPU:
1857 logger.warning(format(" Platform %s not supported; defaulting to FFX", platform));
1858 default:
1859 ForceFieldEnergy ffxEnergy = assembly.getPotentialEnergy();
1860 if (ffxEnergy == null) {
1861 ffxEnergy = new ForceFieldEnergy(assembly, numThreads);
1862 assembly.setPotential(ffxEnergy);
1863 }
1864 return ffxEnergy;
1865 }
1866 } catch (IllegalArgumentException | NullPointerException ex) {
1867 logger.warning(
1868 format(" String %s did not match a known energy implementation", platformString));
1869 ForceFieldEnergy ffxEnergy = assembly.getPotentialEnergy();
1870 if (ffxEnergy == null) {
1871 ffxEnergy = new ForceFieldEnergy(assembly, numThreads);
1872 assembly.setPotential(ffxEnergy);
1873 }
1874 return ffxEnergy;
1875 }
1876 }
1877
1878
1879
1880
1881
1882
1883
1884 public void applyAllConstraintPositions(double[] xPrior, double[] xNew) {
1885 applyAllConstraintPositions(xPrior, xNew, DEFAULT_CONSTRAINT_TOLERANCE);
1886 }
1887
1888
1889
1890
1891
1892
1893
1894
1895 public void applyAllConstraintPositions(double[] xPrior, double[] xNew, double tol) {
1896 if (xPrior == null) {
1897 xPrior = Arrays.copyOf(xNew, xNew.length);
1898 }
1899 for (Constraint constraint : constraints) {
1900 constraint.applyConstraintToStep(xPrior, xNew, getMass(), tol);
1901 }
1902 }
1903
1904
1905
1906
1907
1908
1909
1910 public void attachExtendedSystem(ExtendedSystem system) {
1911 if (system == null) {
1912 throw new IllegalArgumentException();
1913 }
1914 esvTerm = true;
1915 esvTermOrig = esvTerm;
1916 esvSystem = system;
1917 if (vanderWaalsTerm) {
1918 if (vanderWaals == null) {
1919 logger.warning("Null VdW during ESV setup.");
1920 } else {
1921 vanderWaals.attachExtendedSystem(system);
1922 }
1923 }
1924 if (multipoleTerm) {
1925 if (particleMeshEwald == null) {
1926 logger.warning("Null PME during ESV setup.");
1927 } else {
1928 particleMeshEwald.attachExtendedSystem(system);
1929 }
1930 }
1931 }
1932
1933
1934
1935
1936
1937
1938 @Override
1939 public boolean dEdLZeroAtEnds() {
1940
1941
1942 return !lambdaTerm;
1943 }
1944
1945
1946
1947
1948
1949
1950 public boolean destroy() {
1951 if (destroyed) {
1952
1953
1954 logger.fine(format(" This ForceFieldEnergy is already destroyed: %s", this));
1955 return true;
1956 } else {
1957 try {
1958 if (parallelTeam != null) {
1959 parallelTeam.shutdown();
1960 }
1961 if (vanderWaals != null) {
1962 vanderWaals.destroy();
1963 }
1964 if (particleMeshEwald != null) {
1965 particleMeshEwald.destroy();
1966 }
1967 molecularAssembly.finishDestruction();
1968 destroyed = true;
1969 return true;
1970 } catch (Exception ex) {
1971 logger.warning(format(" Exception in shutting down a ForceFieldEnergy: %s", ex));
1972 logger.info(Utilities.stackTraceToString(ex));
1973 return false;
1974 }
1975 }
1976 }
1977
1978
1979
1980
1981
1982
1983 public double energy() {
1984 return energy(false, false);
1985 }
1986
1987
1988
1989
1990
1991
1992 public double energy(boolean gradient, boolean print) {
1993 try {
1994 totalTime = System.nanoTime();
1995 nnTime = 0;
1996 bondTime = 0;
1997 angleTime = 0;
1998 stretchBendTime = 0;
1999 ureyBradleyTime = 0;
2000 outOfPlaneBendTime = 0;
2001 torsionTime = 0;
2002 stretchTorsionTime = 0;
2003 angleTorsionTime = 0;
2004 piOrbitalTorsionTime = 0;
2005 torsionTorsionTime = 0;
2006 improperTorsionTime = 0;
2007 vanDerWaalsTime = 0;
2008 electrostaticTime = 0;
2009 restrainDistanceTime = 0;
2010 restrainTorsionTime = 0;
2011 restrainPositionTime = 0;
2012 restrainGroupTime = 0;
2013 ncsTime = 0;
2014
2015
2016 nnEnergy = 0.0;
2017 bondEnergy = 0.0;
2018 angleEnergy = 0.0;
2019 stretchBendEnergy = 0.0;
2020 ureyBradleyEnergy = 0.0;
2021 outOfPlaneBendEnergy = 0.0;
2022 torsionEnergy = 0.0;
2023 angleTorsionEnergy = 0.0;
2024 stretchTorsionEnergy = 0.0;
2025 piOrbitalTorsionEnergy = 0.0;
2026 torsionTorsionEnergy = 0.0;
2027 improperTorsionEnergy = 0.0;
2028 totalBondedEnergy = 0.0;
2029
2030
2031 restrainDistanceEnergy = 0.0;
2032 restrainTorsionEnergy = 0.0;
2033 restrainPositionEnergy = 0.0;
2034 restrainGroupEnergy = 0.0;
2035 ncsEnergy = 0.0;
2036
2037
2038 bondRMSD = 0.0;
2039 angleRMSD = 0.0;
2040
2041
2042 vanDerWaalsEnergy = 0.0;
2043 permanentMultipoleEnergy = 0.0;
2044 permanentRealSpaceEnergy = 0.0;
2045 polarizationEnergy = 0.0;
2046 totalMultipoleEnergy = 0.0;
2047 totalNonBondedEnergy = 0.0;
2048
2049
2050 solvationEnergy = 0.0;
2051
2052
2053 relativeSolvationEnergy = 0.0;
2054 nRelativeSolvations = 0;
2055
2056 esvBias = 0.0;
2057
2058
2059 totalEnergy = 0.0;
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069
2070 try {
2071 bondedRegion.setGradient(gradient);
2072 parallelTeam.execute(bondedRegion);
2073 } catch (RuntimeException ex) {
2074 logger.warning("Runtime exception during bonded term calculation.");
2075 throw ex;
2076 } catch (Exception ex) {
2077 logger.info(Utilities.stackTraceToString(ex));
2078 logger.severe(ex.toString());
2079 }
2080
2081 if (!lambdaBondedTerms) {
2082
2083 if (ncsTerm) {
2084 ncsTime = -System.nanoTime();
2085 ncsEnergy = ncsRestraint.residual(gradient, print);
2086 ncsTime += System.nanoTime();
2087 }
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099 if (restrainGroupTerm) {
2100 restrainGroupTime = -System.nanoTime();
2101 restrainGroupEnergy = restrainGroups.energy(gradient);
2102 restrainGroupTime += System.nanoTime();
2103 }
2104
2105 if (comTerm) {
2106 comRestraintTime = -System.nanoTime();
2107 comRestraintEnergy = comRestraint.residual(gradient, print);
2108 comRestraintTime += System.nanoTime();
2109 }
2110
2111
2112
2113 if (nnTerm) {
2114 nnTime = -System.nanoTime();
2115 nnEnergy = aniEnergy.energy(gradient, print);
2116 nnTime += System.nanoTime();
2117 }
2118
2119
2120 if (vanderWaalsTerm) {
2121 vanDerWaalsTime = -System.nanoTime();
2122 vanDerWaalsEnergy = vanderWaals.energy(gradient, print);
2123 nVanDerWaalInteractions = this.vanderWaals.getInteractions();
2124 vanDerWaalsTime += System.nanoTime();
2125 }
2126
2127 if (multipoleTerm) {
2128 electrostaticTime = -System.nanoTime();
2129 totalMultipoleEnergy = particleMeshEwald.energy(gradient, print);
2130 permanentMultipoleEnergy = particleMeshEwald.getPermanentEnergy();
2131 permanentRealSpaceEnergy = particleMeshEwald.getPermRealEnergy();
2132 polarizationEnergy = particleMeshEwald.getPolarizationEnergy();
2133 nPermanentInteractions = particleMeshEwald.getInteractions();
2134 solvationEnergy = particleMeshEwald.getSolvationEnergy();
2135 nGKInteractions = particleMeshEwald.getGKInteractions();
2136 electrostaticTime += System.nanoTime();
2137 }
2138 }
2139
2140 if (relativeSolvationTerm) {
2141 List<Residue> residuesList = molecularAssembly.getResidueList();
2142 for (Residue residue : residuesList) {
2143 if (residue instanceof MultiResidue) {
2144 Atom refAtom = residue.getSideChainAtoms().get(0);
2145 if (refAtom != null && refAtom.getUse()) {
2146
2147
2148 double thisSolvation = relativeSolvation.getSolvationEnergy(residue, false);
2149 relativeSolvationEnergy -= thisSolvation;
2150 if (thisSolvation != 0) {
2151 nRelativeSolvations++;
2152 }
2153 }
2154 }
2155 }
2156 }
2157
2158 totalTime = System.nanoTime() - totalTime;
2159
2160 totalBondedEnergy = nnEnergy + bondEnergy + angleEnergy + stretchBendEnergy + ureyBradleyEnergy
2161 + outOfPlaneBendEnergy + torsionEnergy + angleTorsionEnergy + stretchTorsionEnergy
2162 + piOrbitalTorsionEnergy + improperTorsionEnergy + torsionTorsionEnergy
2163 + restrainDistanceEnergy + restrainPositionEnergy + restrainGroupEnergy + restrainTorsionEnergy + ncsEnergy;
2164
2165 totalNonBondedEnergy = vanDerWaalsEnergy + totalMultipoleEnergy + relativeSolvationEnergy;
2166 totalEnergy = totalBondedEnergy + totalNonBondedEnergy + solvationEnergy;
2167 if (esvTerm) {
2168 esvBias = esvSystem.getBiasEnergy();
2169 totalEnergy += esvBias;
2170 }
2171
2172 if (print) {
2173 StringBuilder sb = new StringBuilder();
2174 if (gradient) {
2175 sb.append("\n Computed Potential Energy and Atomic Coordinate Gradients\n");
2176 } else {
2177 sb.append("\n Computed Potential Energy\n");
2178 }
2179 sb.append(this);
2180 logger.info(sb.toString());
2181 }
2182 return totalEnergy;
2183 } catch (EnergyException ex) {
2184 if (printOnFailure) {
2185 printFailure();
2186 }
2187 if (ex.doCauseSevere()) {
2188 logger.info(Utilities.stackTraceToString(ex));
2189 logger.log(Level.SEVERE, " Error in calculating energies or gradients", ex);
2190 } else {
2191 logger.log(Level.INFO, format(" Exception in energy calculation:\n %s", ex));
2192 }
2193 throw ex;
2194 }
2195 }
2196
2197
2198
2199
2200 @Override
2201 public double energy(double[] x) {
2202 return energy(x, false);
2203 }
2204
2205
2206
2207
2208 @Override
2209 public double energy(double[] x, boolean verbose) {
2210 assert Arrays.stream(x).allMatch(Double::isFinite);
2211
2212
2213 unscaleCoordinates(x);
2214
2215
2216 setCoordinates(x);
2217
2218 double e = this.energy(false, verbose);
2219
2220
2221 scaleCoordinates(x);
2222
2223 return e;
2224 }
2225
2226
2227
2228
2229 @Override
2230 public double energyAndGradient(double[] x, double[] g) {
2231 return energyAndGradient(x, g, false);
2232 }
2233
2234
2235
2236
2237 @Override
2238 public double energyAndGradient(double[] x, double[] g, boolean verbose) {
2239 assert Arrays.stream(x).allMatch(Double::isFinite);
2240
2241
2242 unscaleCoordinates(x);
2243
2244
2245 setCoordinates(x);
2246 double e = energy(true, verbose);
2247
2248
2249
2250 try {
2251 fillGradient(g);
2252
2253
2254 scaleCoordinatesAndGradient(x, g);
2255
2256 if (maxDebugGradient < Double.MAX_VALUE) {
2257 boolean extremeGrad = Arrays.stream(g)
2258 .anyMatch((double gi) -> (gi > maxDebugGradient || gi < -maxDebugGradient));
2259 if (extremeGrad) {
2260 File origFile = molecularAssembly.getFile();
2261 String timeString = LocalDateTime.now()
2262 .format(DateTimeFormatter.ofPattern("yyyy_MM_dd-HH_mm_ss"));
2263
2264 String filename = format("%s-LARGEGRAD-%s.pdb",
2265 removeExtension(molecularAssembly.getFile().getName()), timeString);
2266 PotentialsFunctions ef = new PotentialsUtils();
2267 filename = ef.versionFile(filename);
2268
2269 logger.warning(format(" Excessively large gradient detected; printing snapshot to file %s",
2270 filename));
2271 ef.saveAsPDB(molecularAssembly, new File(filename));
2272 molecularAssembly.setFile(origFile);
2273 }
2274 }
2275 return e;
2276 } catch (EnergyException ex) {
2277 if (printOnFailure) {
2278 printFailure();
2279 }
2280 if (ex.doCauseSevere()) {
2281 logger.info(Utilities.stackTraceToString(ex));
2282 logger.log(Level.SEVERE, " Error in calculating energies or gradients", ex);
2283 } else {
2284 logger.log(Level.INFO, format(" Exception in energy calculation:\n %s", ex));
2285 }
2286 throw ex;
2287 }
2288 }
2289
2290
2291
2292
2293 private void printFailure() {
2294 String timeString = LocalDateTime.now()
2295 .format(DateTimeFormatter.ofPattern("yyyy_MM_dd-HH_mm_ss"));
2296 File file = molecularAssembly.getFile();
2297 String ext = "pdb";
2298 if (isXYZ(file)) {
2299 ext = "xyz";
2300 }
2301 String filename = format("%s-ERROR-%s.%s", removeExtension(file.getName()), timeString, ext);
2302 PotentialsFunctions ef = new PotentialsUtils();
2303 filename = ef.versionFile(filename);
2304 logger.info(format(" Writing on-error snapshot to file %s", filename));
2305 ef.save(molecularAssembly, new File(filename));
2306 }
2307
2308
2309
2310
2311
2312
2313 @Override
2314 public double[] getAcceleration(double[] acceleration) {
2315 int n = getNumberOfVariables();
2316 if (acceleration == null || acceleration.length < n) {
2317 acceleration = new double[n];
2318 }
2319 int index = 0;
2320 double[] a = new double[3];
2321 for (int i = 0; i < nAtoms; i++) {
2322 if (atoms[i].isActive()) {
2323 atoms[i].getAcceleration(a);
2324 acceleration[index++] = a[0];
2325 acceleration[index++] = a[1];
2326 acceleration[index++] = a[2];
2327 }
2328 }
2329 return acceleration;
2330 }
2331
2332
2333
2334
2335
2336
2337 public double getAngleEnergy() {
2338 return angleEnergy;
2339 }
2340
2341
2342
2343
2344
2345
2346 public double getAngleTorsionEnergy() {
2347 return angleTorsionEnergy;
2348 }
2349
2350
2351
2352
2353
2354
2355 public AngleTorsion[] getAngleTorsions() {
2356 return angleTorsions;
2357 }
2358
2359
2360
2361
2362
2363
2364 public Angle[] getAngles() {
2365 return angles;
2366 }
2367
2368 public String getAngleEnergyString() {
2369 AngleType angleType = angles[0].angleType;
2370 String energy;
2371 if (angleType.angleFunction == AngleType.AngleFunction.SEXTIC) {
2372 energy = format("""
2373 k*(d^2 + %.15g*d^3 + %.15g*d^4 + %.15g*d^5 + %.15g*d^6);
2374 d=%.15g*theta-theta0;
2375 """,
2376 angleType.cubic, angleType.quartic, angleType.pentic, angleType.sextic, 180.0 / PI);
2377 } else {
2378 energy = format("""
2379 k*(d^2);
2380 d=%.15g*theta-theta0;
2381 """,
2382 180.0 / PI);
2383 }
2384 return energy;
2385 }
2386
2387 public String getInPlaneAngleEnergyString() {
2388 AngleType angleType = angles[0].angleType;
2389 String energy = format("""
2390 k*(d^2 + %.15g*d^3 + %.15g*d^4 + %.15g*d^5 + %.15g*d^6);
2391 d=theta-theta0;
2392 theta = %.15g*pointangle(x1, y1, z1, projx, projy, projz, x3, y3, z3);
2393 projx = x2-nx*dot;
2394 projy = y2-ny*dot;
2395 projz = z2-nz*dot;
2396 dot = nx*(x2-x3) + ny*(y2-y3) + nz*(z2-z3);
2397 nx = px/norm;
2398 ny = py/norm;
2399 nz = pz/norm;
2400 norm = sqrt(px*px + py*py + pz*pz);
2401 px = (d1y*d2z-d1z*d2y);
2402 py = (d1z*d2x-d1x*d2z);
2403 pz = (d1x*d2y-d1y*d2x);
2404 d1x = x1-x4;
2405 d1y = y1-y4;
2406 d1z = z1-z4;
2407 d2x = x3-x4;
2408 d2y = y3-y4;
2409 d2z = z3-z4;
2410 """,
2411 angleType.cubic, angleType.quartic, angleType.pentic, angleType.sextic, 180.0 / PI);
2412 return energy;
2413 }
2414
2415
2416
2417
2418
2419
2420
2421 public Angle[] getAngles(AngleMode angleMode) {
2422 if (angles == null || angles.length < 1) {
2423 return null;
2424 }
2425 int nAngles = angles.length;
2426 List<Angle> angleList = new ArrayList<>();
2427
2428 for (int i = 0; i < nAngles; i++) {
2429 if (angles[i].getAngleMode() == angleMode) {
2430 angleList.add(angles[i]);
2431 }
2432 }
2433 nAngles = angleList.size();
2434 if (nAngles < 1) {
2435 return null;
2436 }
2437 return angleList.toArray(new Angle[0]);
2438 }
2439
2440
2441
2442
2443
2444
2445
2446 public double getNeutralNetworkEnergy() {
2447 return nnEnergy;
2448 }
2449
2450
2451
2452
2453
2454
2455 public double getBondEnergy() {
2456 return bondEnergy;
2457 }
2458
2459
2460
2461
2462
2463
2464 public Bond[] getBonds() {
2465 return bonds;
2466 }
2467
2468 public String getBondEnergyString() {
2469 BondType bondType = bonds[0].getBondType();
2470 String energy;
2471 if (bondType.bondFunction == BondType.BondFunction.QUARTIC) {
2472 energy = format("""
2473 k*(d^2 + %.15g*d^3 + %.15g*d^4);
2474 d=r-r0;
2475 """,
2476 bondType.cubic / OpenMM_NmPerAngstrom,
2477 bondType.quartic / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom));
2478 } else {
2479 energy = """
2480 k*(d^2);
2481 d=r-r0;
2482 """;
2483 }
2484 return energy;
2485 }
2486
2487
2488
2489
2490
2491
2492 public double getCavitationEnergy() {
2493 return particleMeshEwald.getCavitationEnergy();
2494 }
2495
2496
2497
2498
2499
2500
2501 @Override
2502 public List<Constraint> getConstraints() {
2503 return constraints.isEmpty() ? Collections.emptyList() : new ArrayList<>(constraints);
2504 }
2505
2506
2507
2508
2509
2510
2511 public List<RestrainPosition> getRestrainPositions() {
2512 if (restrainPositions == null) {
2513 return Collections.emptyList();
2514 }
2515 return List.of(restrainPositions);
2516 }
2517
2518
2519
2520
2521 @Override
2522 public double[] getCoordinates(double[] x) {
2523 int n = getNumberOfVariables();
2524 if (x == null || x.length < n) {
2525 x = new double[n];
2526 }
2527 int index = 0;
2528 for (int i = 0; i < nAtoms; i++) {
2529 Atom a = atoms[i];
2530 if (a.isActive()) {
2531 x[index++] = a.getX();
2532 x[index++] = a.getY();
2533 x[index++] = a.getZ();
2534 }
2535 }
2536 return x;
2537 }
2538
2539
2540
2541
2542
2543
2544 @Override
2545 public Crystal getCrystal() {
2546 return crystal;
2547 }
2548
2549
2550
2551
2552
2553
2554 @Override
2555 public void setCrystal(Crystal crystal) {
2556 setCrystal(crystal, false);
2557 }
2558
2559
2560
2561
2562
2563
2564 public double getCutoffPlusBuffer() {
2565 return cutoffPlusBuffer;
2566 }
2567
2568
2569
2570
2571
2572
2573 public double getDispersionEnergy() {
2574 return particleMeshEwald.getDispersionEnergy();
2575 }
2576
2577
2578
2579
2580
2581
2582 public double getElectrostaticEnergy() {
2583 return totalMultipoleEnergy;
2584 }
2585
2586
2587
2588
2589 @Override
2590 public STATE getEnergyTermState() {
2591 return state;
2592 }
2593
2594
2595
2596
2597
2598
2599 @Override
2600 public void setEnergyTermState(STATE state) {
2601 this.state = state;
2602 switch (state) {
2603 case FAST:
2604 nnTerm = nnTermOrig;
2605 bondTerm = bondTermOrig;
2606 angleTerm = angleTermOrig;
2607 stretchBendTerm = stretchBendTermOrig;
2608 ureyBradleyTerm = ureyBradleyTermOrig;
2609 outOfPlaneBendTerm = outOfPlaneBendTermOrig;
2610 torsionTerm = torsionTermOrig;
2611 stretchTorsionTerm = stretchTorsionTermOrig;
2612 angleTorsionTerm = angleTorsionTermOrig;
2613 piOrbitalTorsionTerm = piOrbitalTorsionTermOrig;
2614 torsionTorsionTerm = torsionTorsionTermOrig;
2615 improperTorsionTerm = improperTorsionTermOrig;
2616 restrainDistanceTerm = restrainDistanceTermOrig;
2617 restrainTorsionTerm = restrainTorsionTermOrig;
2618 restrainPositionTerm = restrainPositionTermOrig;
2619 restrainGroupTerm = restrainGroupTermOrig;
2620 ncsTerm = ncsTermOrig;
2621 comTerm = comTermOrig;
2622 vanderWaalsTerm = false;
2623 multipoleTerm = false;
2624 polarizationTerm = false;
2625 generalizedKirkwoodTerm = false;
2626 esvTerm = false;
2627 break;
2628 case SLOW:
2629 vanderWaalsTerm = vanderWaalsTermOrig;
2630 multipoleTerm = multipoleTermOrig;
2631 polarizationTerm = polarizationTermOrig;
2632 generalizedKirkwoodTerm = generalizedKirkwoodTermOrig;
2633 esvTerm = esvTermOrig;
2634 nnTerm = false;
2635 bondTerm = false;
2636 angleTerm = false;
2637 stretchBendTerm = false;
2638 ureyBradleyTerm = false;
2639 outOfPlaneBendTerm = false;
2640 torsionTerm = false;
2641 stretchTorsionTerm = false;
2642 angleTorsionTerm = false;
2643 piOrbitalTorsionTerm = false;
2644 torsionTorsionTerm = false;
2645 improperTorsionTerm = false;
2646 restrainDistanceTerm = false;
2647 restrainPositionTerm = false;
2648 restrainGroupTerm = false;
2649 ncsTerm = false;
2650 comTerm = false;
2651 break;
2652 default:
2653 nnTerm = nnTermOrig;
2654 bondTerm = bondTermOrig;
2655 angleTerm = angleTermOrig;
2656 stretchBendTerm = stretchBendTermOrig;
2657 ureyBradleyTerm = ureyBradleyTermOrig;
2658 outOfPlaneBendTerm = outOfPlaneBendTermOrig;
2659 torsionTerm = torsionTermOrig;
2660 stretchTorsionTerm = stretchTorsionTermOrig;
2661 angleTorsionTerm = angleTorsionTermOrig;
2662 piOrbitalTorsionTerm = piOrbitalTorsionTermOrig;
2663 torsionTorsionTerm = torsionTorsionTermOrig;
2664 improperTorsionTerm = improperTorsionTermOrig;
2665 restrainDistanceTerm = restrainDistanceTermOrig;
2666 restrainTorsionTerm = restrainTorsionTermOrig;
2667 restrainPositionTerm = restrainPositionTermOrig;
2668 restrainGroupTerm = restrainGroupTermOrig;
2669 ncsTerm = ncsTermOrig;
2670 comTermOrig = comTerm;
2671 vanderWaalsTerm = vanderWaalsTermOrig;
2672 multipoleTerm = multipoleTermOrig;
2673 polarizationTerm = polarizationTermOrig;
2674 generalizedKirkwoodTerm = generalizedKirkwoodTermOrig;
2675 }
2676 }
2677
2678
2679
2680
2681
2682
2683 public double getEsvBiasEnergy() {
2684 return esvBias;
2685 }
2686
2687
2688
2689
2690
2691
2692 public ExtendedSystem getExtendedSystem() {
2693 return esvSystem;
2694 }
2695
2696
2697
2698
2699
2700
2701 public GeneralizedKirkwood getGK() {
2702 if (particleMeshEwald != null) {
2703 return particleMeshEwald.getGK();
2704 } else {
2705 return null;
2706 }
2707 }
2708
2709
2710
2711
2712
2713
2714 public double getGKEnergy() {
2715 return particleMeshEwald.getGKEnergy();
2716 }
2717
2718
2719
2720
2721
2722
2723
2724 public double[] getGradient(double[] g) {
2725 return fillGradient(g);
2726 }
2727
2728
2729
2730
2731
2732
2733 public double getImproperTorsionEnergy() {
2734 return improperTorsionEnergy;
2735 }
2736
2737
2738
2739
2740
2741
2742 public ImproperTorsion[] getImproperTorsions() {
2743 return improperTorsions;
2744 }
2745
2746
2747
2748
2749 @Override
2750 public double getLambda() {
2751 return lambda;
2752 }
2753
2754
2755
2756
2757 @Override
2758 public void setLambda(double lambda) {
2759 if (lambdaTerm) {
2760 if (lambda <= 1.0 && lambda >= 0.0) {
2761 this.lambda = lambda;
2762 if (vanderWaalsTerm) {
2763 vanderWaals.setLambda(lambda);
2764 }
2765 if (multipoleTerm) {
2766 particleMeshEwald.setLambda(lambda);
2767 }
2768
2769
2770
2771
2772
2773
2774
2775
2776 if (restrainDistanceTerm && nRestrainDistances > 0) {
2777 for (RestrainDistance restrainDistance : restrainDistances) {
2778 restrainDistance.setLambda(lambda);
2779 }
2780 }
2781 if (restrainTorsionTerm && nRestrainTorsions > 0) {
2782 for (RestraintTorsion restraintTorsion : restrainTorsions) {
2783 restraintTorsion.setLambda(lambda);
2784 }
2785 }
2786
2787 if (ncsTerm && ncsRestraint != null) {
2788 ncsRestraint.setLambda(lambda);
2789 }
2790 if (comTerm && comRestraint != null) {
2791 comRestraint.setLambda(lambda);
2792 }
2793 if (lambdaTorsions) {
2794 for (int i = 0; i < nTorsions; i++) {
2795 torsions[i].setLambda(lambda);
2796 }
2797 for (int i = 0; i < nPiOrbitalTorsions; i++) {
2798 piOrbitalTorsions[i].setLambda(lambda);
2799 }
2800 for (int i = 0; i < nTorsionTorsions; i++) {
2801 torsionTorsions[i].setLambda(lambda);
2802 }
2803 }
2804 } else {
2805 String message = format("Lambda value %8.3f is not in the range [0..1].", lambda);
2806 logger.warning(message);
2807 }
2808 } else {
2809 logger.fine(" Attempting to set a lambda value on a ForceFieldEnergy with lambdaterm false.");
2810 }
2811 }
2812
2813
2814
2815
2816 @Override
2817 public double[] getMass() {
2818 int n = getNumberOfVariables();
2819 double[] mass = new double[n];
2820 int index = 0;
2821 for (int i = 0; i < nAtoms; i++) {
2822 Atom a = atoms[i];
2823 if (a.isActive()) {
2824 double m = a.getMass();
2825 mass[index++] = m;
2826 mass[index++] = m;
2827 mass[index++] = m;
2828 }
2829 }
2830 return mass;
2831 }
2832
2833
2834
2835
2836 @Override
2837 public int getNumberOfVariables() {
2838 int nActive = 0;
2839 for (int i = 0; i < nAtoms; i++) {
2840 if (atoms[i].isActive()) {
2841 nActive++;
2842 }
2843 }
2844 return nActive * 3;
2845 }
2846
2847
2848
2849
2850
2851
2852 public int getNumberofAngleTorsions() {
2853 return nAngleTorsions;
2854 }
2855
2856
2857
2858
2859
2860
2861 public int getNumberofAngles() {
2862 return nAngles;
2863 }
2864
2865
2866
2867
2868
2869
2870 public int getNumberofBonds() {
2871 return nBonds;
2872 }
2873
2874
2875
2876
2877
2878
2879 public int getNumberofImproperTorsions() {
2880 return nImproperTorsions;
2881 }
2882
2883
2884
2885
2886
2887
2888 public int getNumberofOutOfPlaneBends() {
2889 return nOutOfPlaneBends;
2890 }
2891
2892
2893
2894
2895
2896
2897 public int getNumberofPiOrbitalTorsions() {
2898 return nPiOrbitalTorsions;
2899 }
2900
2901
2902
2903
2904
2905
2906 public int getNumberofStretchBends() {
2907 return nStretchBends;
2908 }
2909
2910
2911
2912
2913
2914
2915 public int getNumberofStretchTorsions() {
2916 return nStretchTorsions;
2917 }
2918
2919
2920
2921
2922
2923
2924 public int getNumberofTorsionTorsions() {
2925 return nTorsionTorsions;
2926 }
2927
2928
2929
2930
2931
2932
2933 public int getNumberofTorsions() {
2934 return nTorsions;
2935 }
2936
2937
2938
2939
2940
2941
2942 public int getNumberofUreyBradleys() {
2943 return nUreyBradleys;
2944 }
2945
2946
2947
2948
2949
2950
2951 public double getOutOfPlaneBendEnergy() {
2952 return outOfPlaneBendEnergy;
2953 }
2954
2955
2956
2957
2958
2959
2960 public OutOfPlaneBend[] getOutOfPlaneBends() {
2961 return outOfPlaneBends;
2962 }
2963
2964 public String getOutOfPlaneEnergyString() {
2965 OutOfPlaneBendType outOfPlaneBendType = outOfPlaneBends[0].outOfPlaneBendType;
2966 String energy = format("""
2967 k*(theta^2 + %.15g*theta^3 + %.15g*theta^4 + %.15g*theta^5 + %.15g*theta^6);
2968 theta = %.15g*pointangle(x2, y2, z2, x4, y4, z4, projx, projy, projz);
2969 projx = x2-nx*dot;
2970 projy = y2-ny*dot;
2971 projz = z2-nz*dot;
2972 dot = nx*(x2-x3) + ny*(y2-y3) + nz*(z2-z3);
2973 nx = px/norm;
2974 ny = py/norm;
2975 nz = pz/norm;
2976 norm = sqrt(px*px + py*py + pz*pz);
2977 px = (d1y*d2z-d1z*d2y);
2978 py = (d1z*d2x-d1x*d2z);
2979 pz = (d1x*d2y-d1y*d2x);
2980 d1x = x1-x4;
2981 d1y = y1-y4;
2982 d1z = z1-z4;
2983 d2x = x3-x4;
2984 d2y = y3-y4;
2985 d2z = z3-z4
2986 """,
2987 outOfPlaneBendType.cubic, outOfPlaneBendType.quartic,
2988 outOfPlaneBendType.pentic, outOfPlaneBendType.sextic, 180.0 / PI);
2989 return energy;
2990 }
2991
2992
2993
2994
2995
2996
2997 public String getPDBHeaderString() {
2998 energy(false, false);
2999 StringBuilder sb = new StringBuilder();
3000 sb.append("REMARK 3 CALCULATED POTENTIAL ENERGY\n");
3001 if (nnTerm) {
3002 sb.append(
3003 format("REMARK 3 %s %g (%d)\n", "NEUTRAL NETWORK : ", nnEnergy, nAtoms));
3004 }
3005 if (bondTerm) {
3006 sb.append(
3007 format("REMARK 3 %s %g (%d)\n", "BOND STRETCHING : ", bondEnergy, nBonds));
3008 sb.append(format("REMARK 3 %s %g\n", "BOND RMSD : ", bondRMSD));
3009 }
3010 if (angleTerm) {
3011 sb.append(format("REMARK 3 %s %g (%d)\n", "ANGLE BENDING : ", angleEnergy,
3012 nAngles));
3013 sb.append(format("REMARK 3 %s %g\n", "ANGLE RMSD : ", angleRMSD));
3014 }
3015 if (stretchBendTerm) {
3016 sb.append(
3017 format("REMARK 3 %s %g (%d)\n", "STRETCH-BEND : ", stretchBendEnergy,
3018 nStretchBends));
3019 }
3020 if (ureyBradleyTerm) {
3021 sb.append(
3022 format("REMARK 3 %s %g (%d)\n", "UREY-BRADLEY : ", ureyBradleyEnergy,
3023 nUreyBradleys));
3024 }
3025 if (outOfPlaneBendTerm) {
3026 sb.append(
3027 format("REMARK 3 %s %g (%d)\n", "OUT-OF-PLANE BEND : ", outOfPlaneBendEnergy,
3028 nOutOfPlaneBends));
3029 }
3030 if (torsionTerm) {
3031 sb.append(format("REMARK 3 %s %g (%d)\n", "TORSIONAL ANGLE : ", torsionEnergy,
3032 nTorsions));
3033 }
3034 if (piOrbitalTorsionTerm) {
3035 sb.append(format("REMARK 3 %s %g (%d)\n", "PI-ORBITAL TORSION : ",
3036 piOrbitalTorsionEnergy, nPiOrbitalTorsions));
3037 }
3038 if (torsionTorsionTerm) {
3039 sb.append(
3040 format("REMARK 3 %s %g (%d)\n", "TORSION-TORSION : ", torsionTorsionEnergy,
3041 nTorsionTorsions));
3042 }
3043 if (improperTorsionTerm) {
3044 sb.append(
3045 format("REMARK 3 %s %g (%d)\n", "IMPROPER TORSION : ", improperTorsionEnergy,
3046 nImproperTorsions));
3047 }
3048 if (restrainDistanceTerm) {
3049 sb.append(format("REMARK 3 %s %g (%d)\n", "RESTRAINT BOND STRETCHING : ",
3050 restrainDistanceEnergy, nRestrainDistances));
3051 }
3052
3053 if (ncsTerm) {
3054 sb.append(
3055 format("REMARK 3 %s %g (%d)\n", "NCS RESTRAINT : ", ncsEnergy, nAtoms));
3056 }
3057
3058 if (restrainPositionTerm && nRestrainPositions > 0) {
3059 int nRests = 0;
3060 for (RestrainPosition restrainPosition : restrainPositions) {
3061 nRests += restrainPosition.getNumAtoms();
3062 }
3063 sb.append(format("REMARK 3 %s %g (%d)\n", "COORDINATE RESTRAINTS : ", restrainPositionEnergy,
3064 nRests));
3065 }
3066
3067 if (comTerm) {
3068 sb.append(
3069 format("REMARK 3 %s %g (%d)\n", "COM RESTRAINT : ", comRestraintEnergy,
3070 nAtoms));
3071 }
3072
3073 if (vanderWaalsTerm) {
3074 sb.append(
3075 format("REMARK 3 %s %g (%d)\n", "VAN DER WAALS : ", vanDerWaalsEnergy,
3076 nVanDerWaalInteractions));
3077 }
3078 if (multipoleTerm) {
3079 sb.append(format("REMARK 3 %s %g (%d)\n", "ATOMIC MULTIPOLES : ",
3080 permanentMultipoleEnergy, nPermanentInteractions));
3081 }
3082 if (polarizationTerm) {
3083 sb.append(format("REMARK 3 %s %g (%d)\n", "POLARIZATION : ",
3084 polarizationEnergy, nPermanentInteractions));
3085 }
3086 sb.append(format("REMARK 3 %s %g\n", "TOTAL POTENTIAL (KCAL/MOL) : ", totalEnergy));
3087 int nsymm = crystal.getUnitCell().spaceGroup.getNumberOfSymOps();
3088 if (nsymm > 1) {
3089 sb.append(format("REMARK 3 %s %g\n", "UNIT CELL POTENTIAL : ", totalEnergy * nsymm));
3090 }
3091 if (crystal.getUnitCell() != crystal) {
3092 nsymm = crystal.spaceGroup.getNumberOfSymOps();
3093 if (nsymm > 1) {
3094 sb.append(format("REMARK 3 %s %g\n", "REPLICATES CELL POTENTIAL : ", totalEnergy * nsymm));
3095 }
3096 }
3097 sb.append("REMARK 3\n");
3098
3099 return sb.toString();
3100 }
3101
3102
3103
3104
3105
3106
3107 public ParallelTeam getParallelTeam() {
3108 return parallelTeam;
3109 }
3110
3111
3112
3113
3114
3115
3116 public int getPermanentInteractions() {
3117 return nPermanentInteractions;
3118 }
3119
3120
3121
3122
3123
3124
3125 public double getPermanentMultipoleEnergy() {
3126 return permanentMultipoleEnergy;
3127 }
3128
3129
3130
3131
3132
3133
3134 public double getPermanentRealSpaceEnergy() {
3135 return permanentRealSpaceEnergy;
3136 }
3137
3138
3139
3140
3141
3142
3143 public double getPermanentReciprocalMpoleEnergy() {
3144 return particleMeshEwald.getPermRecipEnergy();
3145 }
3146
3147
3148
3149
3150
3151
3152 public double getPermanentReciprocalSelfEnergy() {
3153 return particleMeshEwald.getPermSelfEnergy();
3154 }
3155
3156
3157
3158
3159
3160
3161 public double getPiOrbitalTorsionEnergy() {
3162 return piOrbitalTorsionEnergy;
3163 }
3164
3165
3166
3167
3168
3169
3170 public PiOrbitalTorsion[] getPiOrbitalTorsions() {
3171 return piOrbitalTorsions;
3172 }
3173
3174 public String getPiOrbitalTorsionEnergyString() {
3175 String energy = """
3176 2*k*sin(phi)^2;
3177 phi = pointdihedral(x3+c1x, y3+c1y, z3+c1z, x3, y3, z3, x4, y4, z4, x4+c2x, y4+c2y, z4+c2z);
3178 c1x = (d14y*d24z-d14z*d24y);
3179 c1y = (d14z*d24x-d14x*d24z);
3180 c1z = (d14x*d24y-d14y*d24x);
3181 c2x = (d53y*d63z-d53z*d63y);
3182 c2y = (d53z*d63x-d53x*d63z);
3183 c2z = (d53x*d63y-d53y*d63x);
3184 d14x = x1-x4;
3185 d14y = y1-y4;
3186 d14z = z1-z4;
3187 d24x = x2-x4;
3188 d24y = y2-y4;
3189 d24z = z2-z4;
3190 d53x = x5-x3;
3191 d53y = y5-y3;
3192 d53z = z5-z3;
3193 d63x = x6-x3;
3194 d63y = y6-y3;
3195 d63z = z6-z3;
3196 """;
3197 return energy;
3198 }
3199
3200
3201
3202
3203
3204
3205
3206 public Platform getPlatform() {
3207 return platform;
3208 }
3209
3210
3211
3212
3213
3214
3215 public ParticleMeshEwald getPmeNode() {
3216 return particleMeshEwald;
3217 }
3218
3219
3220
3221
3222
3223
3224 public double getPolarizationEnergy() {
3225 return polarizationEnergy;
3226 }
3227
3228
3229
3230
3231
3232
3233
3234 @Override
3235 public double[] getPreviousAcceleration(double[] previousAcceleration) {
3236 int n = getNumberOfVariables();
3237 if (previousAcceleration == null || previousAcceleration.length < n) {
3238 previousAcceleration = new double[n];
3239 }
3240 int index = 0;
3241 double[] a = new double[3];
3242 for (int i = 0; i < nAtoms; i++) {
3243 if (atoms[i].isActive()) {
3244 atoms[i].getPreviousAcceleration(a);
3245 previousAcceleration[index++] = a[0];
3246 previousAcceleration[index++] = a[1];
3247 previousAcceleration[index++] = a[2];
3248 }
3249 }
3250 return previousAcceleration;
3251 }
3252
3253
3254
3255
3256
3257
3258 public double getRelativeSolvationEnergy() {
3259 return relativeSolvationEnergy;
3260 }
3261
3262
3263
3264
3265 @Override
3266 public double[] getScaling() {
3267 return optimizationScaling;
3268 }
3269
3270
3271
3272
3273 @Override
3274 public void setScaling(double[] scaling) {
3275 optimizationScaling = scaling;
3276 }
3277
3278
3279
3280
3281
3282
3283 public double getSolvationEnergy() {
3284 return solvationEnergy;
3285 }
3286
3287
3288
3289
3290
3291
3292 public int getSolvationInteractions() {
3293 return nGKInteractions;
3294 }
3295
3296
3297
3298
3299
3300
3301 public double getStrenchBendEnergy() {
3302 return stretchBendEnergy;
3303 }
3304
3305
3306
3307
3308
3309
3310 public StretchBend[] getStretchBends() {
3311 return stretchBends;
3312 }
3313
3314 public String getStretchBendEnergyString() {
3315 String energy = format("(k1*(distance(p1,p2)-r12) + k2*(distance(p2,p3)-r23))*(%.15g*(angle(p1,p2,p3)-theta0))", 180.0 / PI);
3316 return energy;
3317 }
3318
3319
3320
3321
3322
3323
3324 public double getStretchTorsionEnergy() {
3325 return stretchTorsionEnergy;
3326 }
3327
3328
3329
3330
3331
3332
3333 public StretchTorsion[] getStretchTorsions() {
3334 return stretchTorsions;
3335 }
3336
3337
3338
3339
3340
3341
3342 public double getTorsionEnergy() {
3343 return torsionEnergy;
3344 }
3345
3346
3347
3348
3349
3350
3351 public double getTorsionTorsionEnergy() {
3352 return torsionTorsionEnergy;
3353 }
3354
3355
3356
3357
3358
3359
3360 public TorsionTorsion[] getTorsionTorsions() {
3361 return torsionTorsions;
3362 }
3363
3364
3365
3366
3367
3368
3369 public Torsion[] getTorsions() {
3370 return torsions;
3371 }
3372
3373
3374
3375
3376
3377
3378 public double getTotalElectrostaticEnergy() {
3379 return totalMultipoleEnergy + solvationEnergy;
3380 }
3381
3382
3383
3384
3385 @Override
3386 public double getTotalEnergy() {
3387 return totalEnergy;
3388 }
3389
3390
3391
3392
3393
3394
3395 public double getUreyBradleyEnergy() {
3396 return ureyBradleyEnergy;
3397 }
3398
3399
3400
3401
3402
3403
3404 public UreyBradley[] getUreyBradleys() {
3405 return ureyBradleys;
3406 }
3407
3408
3409
3410
3411
3412
3413 public double getVanDerWaalsEnergy() {
3414 return vanDerWaalsEnergy;
3415 }
3416
3417
3418
3419
3420
3421
3422 public int getVanDerWaalsInteractions() {
3423 return nVanDerWaalInteractions;
3424 }
3425
3426
3427
3428
3429
3430
3431 @Override
3432 public VARIABLE_TYPE[] getVariableTypes() {
3433 int n = getNumberOfVariables();
3434 VARIABLE_TYPE[] type = new VARIABLE_TYPE[n];
3435 int i = 0;
3436 for (int j = 0; j < nAtoms; j++) {
3437 if (atoms[j].isActive()) {
3438 type[i++] = VARIABLE_TYPE.X;
3439 type[i++] = VARIABLE_TYPE.Y;
3440 type[i++] = VARIABLE_TYPE.Z;
3441 }
3442 }
3443 return type;
3444 }
3445
3446
3447
3448
3449
3450
3451 public VanDerWaals getVdwNode() {
3452 return vanderWaals;
3453 }
3454
3455
3456
3457
3458
3459
3460 @Override
3461 public double[] getVelocity(double[] velocity) {
3462 int n = getNumberOfVariables();
3463 if (velocity == null || velocity.length < n) {
3464 velocity = new double[n];
3465 }
3466 int index = 0;
3467 double[] v = new double[3];
3468 for (int i = 0; i < nAtoms; i++) {
3469 Atom a = atoms[i];
3470 if (a.isActive()) {
3471 a.getVelocity(v);
3472 velocity[index++] = v[0];
3473 velocity[index++] = v[1];
3474 velocity[index++] = v[2];
3475 }
3476 }
3477 return velocity;
3478 }
3479
3480
3481
3482
3483 @Override
3484 public double getd2EdL2() {
3485 double d2EdLambda2 = 0.0;
3486 if (!lambdaBondedTerms) {
3487 if (vanderWaalsTerm) {
3488 d2EdLambda2 = vanderWaals.getd2EdL2();
3489 }
3490 if (multipoleTerm) {
3491 d2EdLambda2 += particleMeshEwald.getd2EdL2();
3492 }
3493
3494
3495
3496
3497
3498
3499
3500 if (restrainDistanceTerm && nRestrainDistances > 0) {
3501 for (int i = 0; i < nRestrainDistances; i++) {
3502 d2EdLambda2 += restrainDistances[i].getd2EdL2();
3503 }
3504 }
3505 if (ncsTerm && ncsRestraint != null) {
3506 d2EdLambda2 += ncsRestraint.getd2EdL2();
3507 }
3508 if (comTerm && comRestraint != null) {
3509 d2EdLambda2 += comRestraint.getd2EdL2();
3510 }
3511 if (lambdaTorsions) {
3512 for (int i = 0; i < nTorsions; i++) {
3513 d2EdLambda2 += torsions[i].getd2EdL2();
3514 }
3515 for (int i = 0; i < nPiOrbitalTorsions; i++) {
3516 d2EdLambda2 += piOrbitalTorsions[i].getd2EdL2();
3517 }
3518 for (int i = 0; i < nTorsionTorsions; i++) {
3519 d2EdLambda2 += torsionTorsions[i].getd2EdL2();
3520 }
3521 }
3522 }
3523 return d2EdLambda2;
3524 }
3525
3526
3527
3528
3529 @Override
3530 public double getdEdL() {
3531 double dEdLambda = 0.0;
3532 if (!lambdaBondedTerms) {
3533 if (vanderWaalsTerm) {
3534 dEdLambda = vanderWaals.getdEdL();
3535 }
3536 if (multipoleTerm) {
3537 dEdLambda += particleMeshEwald.getdEdL();
3538 }
3539 if (restrainDistanceTerm && nRestrainDistances > 0) {
3540 for (int i = 0; i < nRestrainDistances; i++) {
3541 dEdLambda += restrainDistances[i].getdEdL();
3542 }
3543 }
3544
3545
3546
3547
3548
3549
3550
3551 if (ncsTerm && ncsRestraint != null) {
3552 dEdLambda += ncsRestraint.getdEdL();
3553 }
3554 if (comTerm && comRestraint != null) {
3555 dEdLambda += comRestraint.getdEdL();
3556 }
3557 if (lambdaTorsions) {
3558 for (int i = 0; i < nTorsions; i++) {
3559 dEdLambda += torsions[i].getdEdL();
3560 }
3561 for (int i = 0; i < nPiOrbitalTorsions; i++) {
3562 dEdLambda += piOrbitalTorsions[i].getdEdL();
3563 }
3564 for (int i = 0; i < nTorsionTorsions; i++) {
3565 dEdLambda += torsionTorsions[i].getdEdL();
3566 }
3567 }
3568 }
3569 return dEdLambda;
3570 }
3571
3572
3573
3574
3575 @Override
3576 public void getdEdXdL(double[] gradients) {
3577 if (!lambdaBondedTerms) {
3578 if (vanderWaalsTerm) {
3579 vanderWaals.getdEdXdL(gradients);
3580 }
3581 if (multipoleTerm) {
3582 particleMeshEwald.getdEdXdL(gradients);
3583 }
3584 if (restrainDistanceTerm && nRestrainDistances > 0) {
3585 for (int i = 0; i < nRestrainDistances; i++) {
3586 restrainDistances[i].getdEdXdL(gradients);
3587 }
3588 }
3589
3590
3591
3592
3593
3594
3595
3596 if (ncsTerm && ncsRestraint != null) {
3597 ncsRestraint.getdEdXdL(gradients);
3598 }
3599 if (comTerm && comRestraint != null) {
3600 comRestraint.getdEdXdL(gradients);
3601 }
3602 if (lambdaTorsions) {
3603 double[] grad = new double[3];
3604 int index = 0;
3605 for (int i = 0; i < nAtoms; i++) {
3606 Atom a = atoms[i];
3607 if (a.isActive()) {
3608 a.getLambdaXYZGradient(grad);
3609 gradients[index++] += grad[0];
3610 gradients[index++] += grad[1];
3611 gradients[index++] += grad[2];
3612 }
3613 }
3614 }
3615 }
3616 }
3617
3618
3619
3620
3621
3622
3623 @Override
3624 public void setAcceleration(double[] acceleration) {
3625 if (acceleration == null) {
3626 return;
3627 }
3628 int index = 0;
3629 double[] accel = new double[3];
3630 for (int i = 0; i < nAtoms; i++) {
3631 if (atoms[i].isActive()) {
3632 accel[0] = acceleration[index++];
3633 accel[1] = acceleration[index++];
3634 accel[2] = acceleration[index++];
3635 atoms[i].setAcceleration(accel);
3636 }
3637 }
3638 }
3639
3640
3641
3642
3643
3644
3645 public void setCoordinates(double[] coords) {
3646 if (coords == null) {
3647 return;
3648 }
3649 int index = 0;
3650 for (int i = 0; i < nAtoms; i++) {
3651 Atom a = atoms[i];
3652 if (a.isActive()) {
3653 double x = coords[index++];
3654 double y = coords[index++];
3655 double z = coords[index++];
3656 a.moveTo(x, y, z);
3657 }
3658 }
3659 }
3660
3661
3662
3663
3664
3665
3666
3667 public void setCrystal(Crystal crystal, boolean checkReplicatesCell) {
3668 if (checkReplicatesCell) {
3669 this.crystal = ReplicatesCrystal.replicatesCrystalFactory(crystal.getUnitCell(), cutOff2);
3670 } else {
3671 this.crystal = crystal;
3672 }
3673
3674
3675
3676
3677 if (vanderWaalsTerm) {
3678 vanderWaals.setCrystal(this.crystal);
3679 }
3680 if (multipoleTerm) {
3681 particleMeshEwald.setCrystal(this.crystal);
3682 }
3683 }
3684
3685
3686
3687
3688
3689
3690
3691 @Override
3692 public void setPreviousAcceleration(double[] previousAcceleration) {
3693 if (previousAcceleration == null) {
3694 return;
3695 }
3696 int index = 0;
3697 double[] prev = new double[3];
3698 for (int i = 0; i < nAtoms; i++) {
3699 if (atoms[i].isActive()) {
3700 prev[0] = previousAcceleration[index++];
3701 prev[1] = previousAcceleration[index++];
3702 prev[2] = previousAcceleration[index++];
3703 atoms[i].setPreviousAcceleration(prev);
3704 }
3705 }
3706 }
3707
3708
3709
3710
3711
3712
3713
3714
3715
3716
3717
3718 public void setPrintOnFailure(boolean onFail, boolean override) {
3719 if (override) {
3720
3721 printOnFailure = onFail;
3722 } else {
3723 try {
3724 molecularAssembly.getForceField().getBoolean("PRINT_ON_FAILURE");
3725
3726
3727
3728
3729
3730 } catch (Exception ex) {
3731 printOnFailure = onFail;
3732 }
3733 }
3734 }
3735
3736
3737
3738
3739
3740
3741 @Override
3742 public void setVelocity(double[] velocity) {
3743 if (velocity == null) {
3744 return;
3745 }
3746 int index = 0;
3747 double[] vel = new double[3];
3748 for (int i = 0; i < nAtoms; i++) {
3749 if (atoms[i].isActive()) {
3750 vel[0] = velocity[index++];
3751 vel[1] = velocity[index++];
3752 vel[2] = velocity[index++];
3753 atoms[i].setVelocity(vel);
3754 }
3755 }
3756 }
3757
3758
3759
3760
3761 @Override
3762 public String toString() {
3763 StringBuilder sb = new StringBuilder();
3764 if (bondTerm && nBonds > 0) {
3765 sb.append(format(" %s %20.8f %12d %12.3f (%8.5f)\n", "Bond Stretching ", bondEnergy, nBonds,
3766 bondTime * toSeconds, bondRMSD));
3767 }
3768 if (angleTerm && nAngles > 0) {
3769 sb.append(
3770 format(" %s %20.8f %12d %12.3f (%8.5f)\n", "Angle Bending ", angleEnergy, nAngles,
3771 angleTime * toSeconds, angleRMSD));
3772 }
3773 if (stretchBendTerm && nStretchBends > 0) {
3774 sb.append(
3775 format(" %s %20.8f %12d %12.3f\n", "Stretch-Bend ", stretchBendEnergy, nStretchBends,
3776 stretchBendTime * toSeconds));
3777 }
3778 if (ureyBradleyTerm && nUreyBradleys > 0) {
3779 sb.append(
3780 format(" %s %20.8f %12d %12.3f\n", "Urey-Bradley ", ureyBradleyEnergy, nUreyBradleys,
3781 ureyBradleyTime * toSeconds));
3782 }
3783 if (outOfPlaneBendTerm && nOutOfPlaneBends > 0) {
3784 sb.append(format(" %s %20.8f %12d %12.3f\n", "Out-of-Plane Bend ", outOfPlaneBendEnergy,
3785 nOutOfPlaneBends, outOfPlaneBendTime * toSeconds));
3786 }
3787 if (torsionTerm && nTorsions > 0) {
3788 sb.append(format(" %s %20.8f %12d %12.3f\n", "Torsional Angle ", torsionEnergy, nTorsions,
3789 torsionTime * toSeconds));
3790 }
3791 if (piOrbitalTorsionTerm && nPiOrbitalTorsions > 0) {
3792 sb.append(format(" %s %20.8f %12d %12.3f\n", "Pi-Orbital Torsion", piOrbitalTorsionEnergy,
3793 nPiOrbitalTorsions, piOrbitalTorsionTime * toSeconds));
3794 }
3795 if (stretchTorsionTerm && nStretchTorsions > 0) {
3796 sb.append(format(" %s %20.8f %12d %12.3f\n", "Stretch-Torsion ", stretchTorsionEnergy,
3797 nStretchTorsions, stretchTorsionTime * toSeconds));
3798 }
3799 if (angleTorsionTerm && nAngleTorsions > 0) {
3800 sb.append(format(" %s %20.8f %12d %12.3f\n", "Angle-Torsion ", angleTorsionEnergy,
3801 nAngleTorsions, angleTorsionTime * toSeconds));
3802 }
3803 if (torsionTorsionTerm && nTorsionTorsions > 0) {
3804 sb.append(format(" %s %20.8f %12d %12.3f\n", "Torsion-Torsion ", torsionTorsionEnergy,
3805 nTorsionTorsions, torsionTorsionTime * toSeconds));
3806 }
3807 if (improperTorsionTerm && nImproperTorsions > 0) {
3808 sb.append(format(" %s %20.8f %12d %12.3f\n", "Improper Torsion ", improperTorsionEnergy,
3809 nImproperTorsions, improperTorsionTime * toSeconds));
3810 }
3811 if (restrainDistanceTerm && nRestrainDistances > 0) {
3812 sb.append(format(" %s %20.8f %12d %12.3f\n", "Restrain Distance ", restrainDistanceEnergy,
3813 nRestrainDistances, restrainDistanceTime * toSeconds));
3814 }
3815 if (restrainTorsionTerm && nRestrainTorsions > 0) {
3816 sb.append(format(" %s %20.8f %12d %12.3f\n", "Restrain Torsion ", restrainTorsionEnergy, nRestrainTorsions,
3817 restrainTorsionTime * toSeconds));
3818 }
3819 if (restrainPositionTerm && nRestrainPositions > 0) {
3820 int nRestrainPositionAtoms = 0;
3821 for (RestrainPosition restrainPosition : restrainPositions) {
3822 nRestrainPositionAtoms += restrainPosition.getNumAtoms();
3823 }
3824 sb.append(format(" %s %20.8f %12d %12.3f\n", "Restrain Position ", restrainPositionEnergy, nRestrainPositionAtoms,
3825 restrainPositionTime * toSeconds));
3826 }
3827 if (restrainGroupTerm && nRestrainGroups > 0) {
3828 sb.append(format(" %s %20.8f %12d %12.3f\n", "Restrain Groups ", restrainGroupEnergy,
3829 nRestrainGroups, restrainGroupTime * toSeconds));
3830 }
3831 if (ncsTerm) {
3832 sb.append(format(" %s %20.8f %12d %12.3f\n", "NCS Restraint ", ncsEnergy, nAtoms,
3833 ncsTime * toSeconds));
3834 }
3835 if (comTerm) {
3836 sb.append(format(" %s %20.8f %12d %12.3f\n", "COM Restraint ", comRestraintEnergy, nAtoms,
3837 comRestraintTime * toSeconds));
3838 }
3839 if (vanderWaalsTerm && nVanDerWaalInteractions > 0) {
3840 sb.append(format(" %s %20.8f %12d %12.3f\n", "Van der Waals ", vanDerWaalsEnergy,
3841 nVanDerWaalInteractions, vanDerWaalsTime * toSeconds));
3842 }
3843 if (multipoleTerm && nPermanentInteractions > 0) {
3844 if (polarizationTerm) {
3845 sb.append(format(" %s %20.8f %12d\n", "Atomic Multipoles ", permanentMultipoleEnergy,
3846 nPermanentInteractions));
3847 } else {
3848 if (elecForm == ELEC_FORM.FIXED_CHARGE) {
3849 sb.append(format(" %s %20.8f %12d %12.3f\n", "Atomic Charges ", permanentMultipoleEnergy,
3850 nPermanentInteractions, electrostaticTime * toSeconds));
3851 } else
3852 sb.append(format(" %s %20.8f %12d %12.3f\n", "Atomic Multipoles ", permanentMultipoleEnergy,
3853 nPermanentInteractions, electrostaticTime * toSeconds));
3854 }
3855 }
3856 if (polarizationTerm && nPermanentInteractions > 0) {
3857 sb.append(format(" %s %20.8f %12d %12.3f\n", "Polarization ", polarizationEnergy,
3858 nPermanentInteractions, electrostaticTime * toSeconds));
3859 }
3860 if (generalizedKirkwoodTerm && nGKInteractions > 0) {
3861 sb.append(
3862 format(" %s %20.8f %12d\n", "Solvation ", solvationEnergy, nGKInteractions));
3863 }
3864 if (relativeSolvationTerm) {
3865 sb.append(format(" %s %20.8f %12d\n", "Relative Solvation", relativeSolvationEnergy,
3866 nRelativeSolvations));
3867 }
3868 if (esvTerm) {
3869 sb.append(
3870 format(" %s %20.8f %s\n", "ExtendedSystemBias", esvBias, esvSystem.getLambdaList()));
3871 sb.append(esvSystem.getBiasDecomposition());
3872 }
3873 if (nnTerm) {
3874 sb.append(format(" %s %20.8f %12d %12.3f\n", "Neural Network ", nnEnergy, nAtoms,
3875 nnTime * toSeconds));
3876 }
3877 sb.append(
3878 format(" %s %20.8f %s %12.3f (sec)", "Total Potential ", totalEnergy, "(Kcal/mole)",
3879 totalTime * toSeconds));
3880 int nsymm = crystal.getUnitCell().spaceGroup.getNumberOfSymOps();
3881 if (nsymm > 1) {
3882 sb.append(format("\n %s %20.8f", "Unit Cell ", totalEnergy * nsymm));
3883 }
3884 if (crystal.getUnitCell() != crystal) {
3885 nsymm = crystal.spaceGroup.getNumberOfSymOps();
3886 sb.append(format("\n %s %20.8f", "Replicates Cell ", totalEnergy * nsymm));
3887 }
3888
3889 return sb.toString();
3890 }
3891
3892
3893
3894
3895
3896 public void logBondedTerms() {
3897 if (bondTerm && nBonds > 0) {
3898 logger.info("\n Bond Stretching Interactions:");
3899 Bond[] bonds = getBonds();
3900 for (Bond bond : bonds) {
3901 logger.info(" Bond \t" + bond.toString());
3902 }
3903 }
3904 if (angleTerm && nAngles > 0) {
3905 logger.info("\n Angle Bending Interactions:");
3906 Angle[] angles = getAngles();
3907 for (Angle angle : angles) {
3908 logger.info(" Angle \t" + angle.toString());
3909 }
3910 }
3911 if (stretchBendTerm && nStretchBends > 0) {
3912 logger.info("\n Stretch-Bend Interactions:");
3913 StretchBend[] stretchBends = getStretchBends();
3914 for (StretchBend stretchBend : stretchBends) {
3915 logger.info(" Stretch-Bend \t" + stretchBend.toString());
3916 }
3917 }
3918 if (ureyBradleyTerm && nUreyBradleys > 0) {
3919 logger.info("\n Urey-Bradley Interactions:");
3920 UreyBradley[] ureyBradleys = getUreyBradleys();
3921 for (UreyBradley ureyBradley : ureyBradleys) {
3922 logger.info("Urey-Bradley \t" + ureyBradley.toString());
3923 }
3924 }
3925 if (outOfPlaneBendTerm && nOutOfPlaneBends > 0) {
3926 logger.info("\n Out-of-Plane Bend Interactions:");
3927 OutOfPlaneBend[] outOfPlaneBends = getOutOfPlaneBends();
3928 for (OutOfPlaneBend outOfPlaneBend : outOfPlaneBends) {
3929 logger.info(" Out-of-Plane Bend \t" + outOfPlaneBend.toString());
3930 }
3931 }
3932 if (torsionTerm && nTorsions > 0) {
3933 logger.info("\n Torsion Angle Interactions:");
3934 Torsion[] torsions = getTorsions();
3935 for (Torsion torsion : torsions) {
3936 logger.info(" Torsion \t" + torsion.toString());
3937 }
3938 }
3939 if (piOrbitalTorsionTerm && nPiOrbitalTorsions > 0) {
3940 logger.info("\n Pi-Orbital Torsion Interactions:");
3941 PiOrbitalTorsion[] piOrbitalTorsions = getPiOrbitalTorsions();
3942 for (PiOrbitalTorsion piOrbitalTorsion : piOrbitalTorsions) {
3943 logger.info(" Pi-Torsion \t" + piOrbitalTorsion.toString());
3944 }
3945 }
3946 if (stretchTorsionTerm && nStretchTorsions > 0) {
3947 logger.info("\n Stretch-Torsion Interactions:");
3948 StretchTorsion[] stretchTorsions = getStretchTorsions();
3949 for (StretchTorsion stretchTorsion : stretchTorsions) {
3950 logger.info(" Stretch-Torsion \t" + stretchTorsion.toString());
3951 }
3952 }
3953 if (angleTorsionTerm && nAngleTorsions > 0) {
3954 logger.info("\n Angle-Torsion Interactions:");
3955 AngleTorsion[] angleTorsions = getAngleTorsions();
3956 for (AngleTorsion angleTorsion : angleTorsions) {
3957 logger.info(" Angle-Torsion \t" + angleTorsion.toString());
3958 }
3959 }
3960 if (torsionTorsionTerm && nTorsionTorsions > 0) {
3961 logger.info("\n Torsion-Torsion Interactions:");
3962 TorsionTorsion[] torsionTorsions = getTorsionTorsions();
3963 for (TorsionTorsion torsionTorsion : torsionTorsions) {
3964 logger.info(" Torsion-Torsion \t" + torsionTorsion.toString());
3965 }
3966 }
3967 if (improperTorsionTerm && nImproperTorsions > 0) {
3968 logger.info("\n Improper Interactions:");
3969 ImproperTorsion[] improperTorsions = getImproperTorsions();
3970 for (ImproperTorsion improperTorsion : improperTorsions) {
3971 logger.info(" Improper \t" + improperTorsion.toString());
3972 }
3973 }
3974 if (restrainDistanceTerm && nRestrainDistances > 0) {
3975 logger.info("\n Restrain Distance Interactions:");
3976 List<RestrainDistance> restrainDistances = getRestrainDistances(null);
3977 for (RestrainDistance restrainDistance : restrainDistances) {
3978 logger.info(" Restrain Distance \t" + restrainDistance.toString());
3979 }
3980 }
3981 if (restrainTorsionTerm && nRestrainTorsions > 0) {
3982 logger.info("\n Restrain Torsion Interactions:");
3983 for (RestraintTorsion restraintTorsion : restrainTorsions) {
3984 logger.info(" Restrain Torsion \t" + restraintTorsion.toString());
3985 }
3986 }
3987 if (restrainPositionTerm && nRestrainPositions > 0) {
3988 logger.info("\n Restrain Position Interactions:");
3989 for (RestrainPosition restrainPosition : restrainPositions) {
3990 logger.info(" Restrain Position \t" + restrainPosition.toString());
3991 }
3992 }
3993 if (restrainGroupTerm && nRestrainGroups > 0) {
3994 logger.info("\n Restrain Group Interactions:");
3995 logger.info(restrainGroups.toString());
3996 }
3997 }
3998
3999
4000
4001
4002
4003
4004
4005 void setLambdaBondedTerms(boolean lambdaBondedTerms, boolean lambdaAllBondedTerms) {
4006 this.lambdaBondedTerms = lambdaBondedTerms;
4007 this.lambdaAllBondedTerms = lambdaAllBondedTerms;
4008 }
4009
4010
4011
4012
4013
4014
4015
4016 private double getNonbondedEnergy(boolean includeSolv) {
4017 return (includeSolv ? (totalNonBondedEnergy + solvationEnergy) : totalNonBondedEnergy);
4018 }
4019
4020
4021
4022
4023
4024
4025
4026
4027 private double[] fillGradient(double[] g) {
4028 assert (g != null);
4029 double[] grad = new double[3];
4030 int n = getNumberOfVariables();
4031 if (g == null || g.length < n) {
4032 g = new double[n];
4033 }
4034 int index = 0;
4035 for (int i = 0; i < nAtoms; i++) {
4036 Atom a = atoms[i];
4037 if (a.isActive()) {
4038 a.getXYZGradient(grad);
4039 double gx = grad[0];
4040 double gy = grad[1];
4041 double gz = grad[2];
4042 if (isNaN(gx) || isInfinite(gx) || isNaN(gy) || isInfinite(gy) || isNaN(gz) || isInfinite(
4043 gz)) {
4044 StringBuilder sb = new StringBuilder(
4045 format("The gradient of atom %s is (%8.3f,%8.3f,%8.3f).", a, gx, gy, gz));
4046 double[] vals = new double[3];
4047 a.getVelocity(vals);
4048 sb.append(format("\n Velocities: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
4049 a.getAcceleration(vals);
4050 sb.append(format("\n Accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
4051 a.getPreviousAcceleration(vals);
4052 sb.append(
4053 format("\n Previous accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
4054
4055 throw new EnergyException(sb.toString());
4056 }
4057 g[index++] = gx;
4058 g[index++] = gy;
4059 g[index++] = gz;
4060 }
4061 }
4062 return g;
4063 }
4064
4065
4066
4067
4068
4069
4070
4071
4072
4073
4074
4075
4076
4077 private void setRestrainDistance(Atom a1, Atom a2, double distance, double forceConstant,
4078 double flatBottom, double lamStart, double lamEnd,
4079 UnivariateSwitchingFunction switchingFunction) {
4080 restrainDistanceTerm = true;
4081 boolean rbLambda = !(switchingFunction instanceof ConstantSwitch) && lambdaTerm;
4082 RestrainDistance restrainDistance = new RestrainDistance(a1, a2, crystal, rbLambda, lamStart, lamEnd, switchingFunction);
4083 int[] classes = {a1.getAtomType().atomClass, a2.getAtomType().atomClass};
4084 if (flatBottom != 0) {
4085 BondType bondType = new BondType(classes, forceConstant, distance,
4086 BondType.BondFunction.FLAT_BOTTOM_HARMONIC, flatBottom);
4087 restrainDistance.setBondType(bondType);
4088 } else {
4089 BondType bondType = new BondType(classes, forceConstant, distance,
4090 BondType.BondFunction.HARMONIC);
4091 restrainDistance.setBondType(bondType);
4092 }
4093
4094 if (restrainDistances == null) {
4095 nRestrainDistances = 0;
4096 }
4097
4098 RestrainDistance[] restrainDistances1 = new RestrainDistance[++nRestrainDistances];
4099 if (restrainDistances != null && restrainDistances.length != 0) {
4100 arraycopy(restrainDistances, 0, restrainDistances1, 0, (nRestrainDistances - 1));
4101 }
4102 restrainDistances1[nRestrainDistances - 1] = restrainDistance;
4103 restrainDistances = restrainDistances1;
4104 restrainDistance.energy(false);
4105 restrainDistance.log();
4106 }
4107
4108 private Crystal configureNCS(ForceField forceField, Crystal unitCell) {
4109
4110
4111 if (forceField.getProperties().containsKey("MTRIX1")
4112 && forceField.getProperties().containsKey("MTRIX2") && forceField.getProperties()
4113 .containsKey("MTRIX3")) {
4114 Crystal unitCell2 = new Crystal(unitCell.a, unitCell.b, unitCell.c, unitCell.alpha,
4115 unitCell.beta, unitCell.gamma, unitCell.spaceGroup.pdbName);
4116 SpaceGroup spaceGroup = unitCell2.spaceGroup;
4117
4118 CompositeConfiguration properties = forceField.getProperties();
4119 String[] MTRX1List = properties.getStringArray("MTRIX1");
4120 String[] MTRX2List = properties.getStringArray("MTRIX2");
4121 String[] MTRX3List = properties.getStringArray("MTRIX3");
4122 spaceGroup.symOps.clear();
4123 double number1;
4124 double number2;
4125 double number3;
4126 for (int i = 0; i < MTRX1List.length; i++) {
4127 double[][] Rot_MTRX = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
4128 double[] Tr_MTRX = {0, 0, 0};
4129 String[] tokens1 = MTRX1List[i].trim().split(" +");
4130 String[] tokens2 = MTRX2List[i].trim().split(" +");
4131 String[] tokens3 = MTRX3List[i].trim().split(" +");
4132 for (int k = 0; k < 4; k++) {
4133 number1 = Double.parseDouble(tokens1[k]);
4134 number2 = Double.parseDouble(tokens2[k]);
4135 number3 = Double.parseDouble(tokens3[k]);
4136 if (k != 3) {
4137 Rot_MTRX[0][k] = number1;
4138 Rot_MTRX[1][k] = number2;
4139 Rot_MTRX[2][k] = number3;
4140 } else {
4141 Tr_MTRX[0] = number1;
4142 Tr_MTRX[1] = number2;
4143 Tr_MTRX[2] = number3;
4144 }
4145 }
4146 SymOp symOp = new SymOp(Rot_MTRX, Tr_MTRX);
4147 if (logger.isLoggable(Level.FINEST)) {
4148 logger.info(
4149 format(" MTRIXn SymOp: %d of %d\n" + symOp, i + 1, MTRX1List.length));
4150 }
4151 spaceGroup.symOps.add(symOp);
4152 }
4153 unitCell = NCSCrystal.NCSCrystalFactory(unitCell, spaceGroup.symOps);
4154 unitCell.updateCrystal();
4155 }
4156 return unitCell;
4157 }
4158
4159
4160
4161
4162
4163
4164
4165 public List<RestrainDistance> getRestrainDistances(@Nullable BondType.BondFunction bondFunction) {
4166 List<RestrainDistance> list = new ArrayList<>();
4167 if (restrainDistances != null && restrainDistances.length > 0) {
4168
4169 if (bondFunction == null) {
4170 return Arrays.asList(restrainDistances);
4171 }
4172
4173 for (RestrainDistance restrainDistance : restrainDistances) {
4174 if (restrainDistance.getBondType().bondFunction == bondFunction) {
4175 list.add(restrainDistance);
4176 }
4177 }
4178 if (!list.isEmpty()) {
4179 return list;
4180 }
4181 }
4182 return null;
4183 }
4184
4185
4186
4187
4188
4189
4190 public List<RestraintTorsion> getRestrainTorsions() {
4191 if (restrainTorsions != null && restrainTorsions.length > 0) {
4192 return Arrays.asList(restrainTorsions);
4193 } else {
4194 return null;
4195 }
4196 }
4197
4198
4199
4200
4201
4202
4203 public RestrainGroups getRestrainGroups() {
4204 return restrainGroups;
4205 }
4206
4207 private class BondedRegion extends ParallelRegion {
4208
4209
4210 private final SharedDouble sharedBondRMSD;
4211 private final SharedDouble sharedAngleRMSD;
4212
4213 private final SharedDouble sharedBondEnergy;
4214 private final SharedDouble sharedAngleEnergy;
4215 private final SharedDouble sharedOutOfPlaneBendEnergy;
4216 private final SharedDouble sharedPiOrbitalTorsionEnergy;
4217 private final SharedDouble sharedStretchBendEnergy;
4218 private final SharedDouble sharedUreyBradleyEnergy;
4219 private final SharedDouble sharedImproperTorsionEnergy;
4220 private final SharedDouble sharedTorsionEnergy;
4221 private final SharedDouble sharedStretchTorsionEnergy;
4222 private final SharedDouble sharedAngleTorsionEnergy;
4223 private final SharedDouble sharedTorsionTorsionEnergy;
4224 private final SharedDouble sharedRestrainPositionEnergy;
4225 private final SharedDouble sharedRestrainDistanceEnergy;
4226 private final SharedDouble sharedRestrainTorsionEnergy;
4227
4228 private final int nThreads;
4229
4230 private final GradInitLoop[] gradInitLoops;
4231 private final GradReduceLoop[] gradReduceLoops;
4232
4233 private final BondedTermLoop[] bondLoops;
4234 private final BondedTermLoop[] angleLoops;
4235 private final BondedTermLoop[] outOfPlaneBendLoops;
4236 private final BondedTermLoop[] improperTorsionLoops;
4237 private final BondedTermLoop[] piOrbitalTorsionLoops;
4238 private final BondedTermLoop[] stretchBendLoops;
4239 private final BondedTermLoop[] torsionLoops;
4240 private final BondedTermLoop[] stretchTorsionLoops;
4241 private final BondedTermLoop[] angleTorsionLoops;
4242 private final BondedTermLoop[] torsionTorsionLoops;
4243 private final BondedTermLoop[] ureyBradleyLoops;
4244
4245 private final BondedTermLoop[] restrainPositionLoops;
4246 private final BondedTermLoop[] restrainDistanceLoops;
4247 private final BondedTermLoop[] restrainTorsionLoops;
4248 private final AtomicDoubleArray3D grad;
4249
4250 private boolean gradient = false;
4251 private AtomicDoubleArrayImpl atomicDoubleArrayImpl;
4252 private AtomicDoubleArray3D lambdaGrad;
4253
4254 BondedRegion() {
4255
4256
4257 sharedAngleRMSD = new SharedDouble();
4258 sharedBondRMSD = new SharedDouble();
4259
4260
4261 sharedBondEnergy = new SharedDouble();
4262 sharedAngleEnergy = new SharedDouble();
4263 sharedImproperTorsionEnergy = new SharedDouble();
4264 sharedOutOfPlaneBendEnergy = new SharedDouble();
4265 sharedPiOrbitalTorsionEnergy = new SharedDouble();
4266 sharedStretchBendEnergy = new SharedDouble();
4267 sharedTorsionEnergy = new SharedDouble();
4268 sharedStretchTorsionEnergy = new SharedDouble();
4269 sharedAngleTorsionEnergy = new SharedDouble();
4270 sharedTorsionTorsionEnergy = new SharedDouble();
4271 sharedUreyBradleyEnergy = new SharedDouble();
4272
4273
4274 sharedRestrainPositionEnergy = new SharedDouble();
4275 sharedRestrainDistanceEnergy = new SharedDouble();
4276 sharedRestrainTorsionEnergy = new SharedDouble();
4277
4278 nThreads = parallelTeam.getThreadCount();
4279
4280
4281 gradInitLoops = new GradInitLoop[nThreads];
4282 gradReduceLoops = new GradReduceLoop[nThreads];
4283
4284
4285 angleLoops = new BondedTermLoop[nThreads];
4286 bondLoops = new BondedTermLoop[nThreads];
4287 improperTorsionLoops = new BondedTermLoop[nThreads];
4288 outOfPlaneBendLoops = new BondedTermLoop[nThreads];
4289 piOrbitalTorsionLoops = new BondedTermLoop[nThreads];
4290 stretchBendLoops = new BondedTermLoop[nThreads];
4291 torsionLoops = new BondedTermLoop[nThreads];
4292 stretchTorsionLoops = new BondedTermLoop[nThreads];
4293 angleTorsionLoops = new BondedTermLoop[nThreads];
4294 torsionTorsionLoops = new BondedTermLoop[nThreads];
4295 ureyBradleyLoops = new BondedTermLoop[nThreads];
4296 restrainPositionLoops = new BondedTermLoop[nThreads];
4297 restrainTorsionLoops = new BondedTermLoop[nThreads];
4298 restrainDistanceLoops = new BondedTermLoop[nThreads];
4299
4300
4301 atomicDoubleArrayImpl = AtomicDoubleArrayImpl.MULTI;
4302 ForceField forceField = molecularAssembly.getForceField();
4303
4304 String value = forceField.getString("ARRAY_REDUCTION", "MULTI");
4305 try {
4306 atomicDoubleArrayImpl = AtomicDoubleArrayImpl.valueOf(toEnumForm(value));
4307 } catch (Exception e) {
4308 logger.info(format(" Unrecognized ARRAY-REDUCTION %s; defaulting to %s", value,
4309 atomicDoubleArrayImpl));
4310 }
4311 logger.fine(format(" Bonded using %s arrays.", atomicDoubleArrayImpl.toString()));
4312
4313 grad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, nThreads);
4314 if (lambdaTerm) {
4315 lambdaGrad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, nThreads);
4316 }
4317 }
4318
4319 @Override
4320 public void finish() {
4321
4322 if (bondTerm) {
4323 bondRMSD = sqrt(sharedBondRMSD.get() / bonds.length);
4324 }
4325 if (angleTerm) {
4326 angleRMSD = sqrt(sharedAngleRMSD.get() / angles.length);
4327 }
4328
4329
4330 angleEnergy = sharedAngleEnergy.get();
4331 bondEnergy = sharedBondEnergy.get();
4332 improperTorsionEnergy = sharedImproperTorsionEnergy.get();
4333 outOfPlaneBendEnergy = sharedOutOfPlaneBendEnergy.get();
4334 piOrbitalTorsionEnergy = sharedPiOrbitalTorsionEnergy.get();
4335 stretchBendEnergy = sharedStretchBendEnergy.get();
4336 ureyBradleyEnergy = sharedUreyBradleyEnergy.get();
4337 torsionEnergy = sharedTorsionEnergy.get();
4338 stretchTorsionEnergy = sharedStretchTorsionEnergy.get();
4339 angleTorsionEnergy = sharedAngleTorsionEnergy.get();
4340 torsionTorsionEnergy = sharedTorsionTorsionEnergy.get();
4341 ureyBradleyEnergy = sharedUreyBradleyEnergy.get();
4342
4343
4344 restrainPositionEnergy = sharedRestrainPositionEnergy.get();
4345 restrainDistanceEnergy = sharedRestrainDistanceEnergy.get();
4346 restrainTorsionEnergy = sharedRestrainTorsionEnergy.get();
4347 }
4348
4349 @Override
4350 public void run() throws Exception {
4351 int threadID = getThreadIndex();
4352
4353
4354 if (gradient) {
4355 if (gradInitLoops[threadID] == null) {
4356 gradInitLoops[threadID] = new GradInitLoop();
4357 }
4358 execute(0, nAtoms - 1, gradInitLoops[threadID]);
4359 }
4360
4361
4362 if (angleTerm) {
4363 if (angleLoops[threadID] == null) {
4364 angleLoops[threadID] = new BondedTermLoop(angles, sharedAngleEnergy, sharedAngleRMSD);
4365 }
4366 if (threadID == 0) {
4367 angleTime = -System.nanoTime();
4368 }
4369 execute(0, nAngles - 1, angleLoops[threadID]);
4370 if (threadID == 0) {
4371 angleTime += System.nanoTime();
4372 }
4373 }
4374
4375 if (bondTerm) {
4376 if (bondLoops[threadID] == null) {
4377 bondLoops[threadID] = new BondedTermLoop(bonds, sharedBondEnergy, sharedBondRMSD);
4378 }
4379 if (threadID == 0) {
4380 bondTime = -System.nanoTime();
4381 }
4382 execute(0, nBonds - 1, bondLoops[threadID]);
4383 if (threadID == 0) {
4384 bondTime += System.nanoTime();
4385 }
4386 }
4387
4388 if (improperTorsionTerm) {
4389 if (improperTorsionLoops[threadID] == null) {
4390 improperTorsionLoops[threadID] = new BondedTermLoop(improperTorsions,
4391 sharedImproperTorsionEnergy);
4392 }
4393 if (threadID == 0) {
4394 improperTorsionTime = -System.nanoTime();
4395 }
4396 execute(0, nImproperTorsions - 1, improperTorsionLoops[threadID]);
4397 if (threadID == 0) {
4398 improperTorsionTime += System.nanoTime();
4399 }
4400 }
4401
4402 if (outOfPlaneBendTerm) {
4403 if (outOfPlaneBendLoops[threadID] == null) {
4404 outOfPlaneBendLoops[threadID] = new BondedTermLoop(outOfPlaneBends,
4405 sharedOutOfPlaneBendEnergy);
4406 }
4407 if (threadID == 0) {
4408 outOfPlaneBendTime = -System.nanoTime();
4409 }
4410 execute(0, nOutOfPlaneBends - 1, outOfPlaneBendLoops[threadID]);
4411 if (threadID == 0) {
4412 outOfPlaneBendTime += System.nanoTime();
4413 }
4414 }
4415
4416 if (piOrbitalTorsionTerm) {
4417 if (piOrbitalTorsionLoops[threadID] == null) {
4418 piOrbitalTorsionLoops[threadID] = new BondedTermLoop(piOrbitalTorsions,
4419 sharedPiOrbitalTorsionEnergy);
4420 }
4421 if (threadID == 0) {
4422 piOrbitalTorsionTime = -System.nanoTime();
4423 }
4424 execute(0, nPiOrbitalTorsions - 1, piOrbitalTorsionLoops[threadID]);
4425 if (threadID == 0) {
4426 piOrbitalTorsionTime += System.nanoTime();
4427 }
4428 }
4429
4430 if (stretchBendTerm) {
4431 if (stretchBendLoops[threadID] == null) {
4432 stretchBendLoops[threadID] = new BondedTermLoop(stretchBends, sharedStretchBendEnergy);
4433 }
4434 if (threadID == 0) {
4435 stretchBendTime = -System.nanoTime();
4436 }
4437 execute(0, nStretchBends - 1, stretchBendLoops[threadID]);
4438 if (threadID == 0) {
4439 stretchBendTime += System.nanoTime();
4440 }
4441 }
4442
4443 if (torsionTerm) {
4444 if (torsionLoops[threadID] == null) {
4445 torsionLoops[threadID] = new BondedTermLoop(torsions, sharedTorsionEnergy);
4446 }
4447 if (threadID == 0) {
4448 torsionTime = -System.nanoTime();
4449 }
4450 execute(0, nTorsions - 1, torsionLoops[threadID]);
4451 if (threadID == 0) {
4452 torsionTime += System.nanoTime();
4453 }
4454 }
4455
4456 if (stretchTorsionTerm) {
4457 if (stretchTorsionLoops[threadID] == null) {
4458 stretchTorsionLoops[threadID] = new BondedTermLoop(stretchTorsions,
4459 sharedStretchTorsionEnergy);
4460 }
4461 if (threadID == 0) {
4462 stretchTorsionTime = -System.nanoTime();
4463 }
4464 execute(0, nStretchTorsions - 1, stretchTorsionLoops[threadID]);
4465 if (threadID == 0) {
4466 stretchTorsionTime += System.nanoTime();
4467 }
4468 }
4469
4470 if (angleTorsionTerm) {
4471 if (angleTorsionLoops[threadID] == null) {
4472 angleTorsionLoops[threadID] = new BondedTermLoop(angleTorsions, sharedAngleTorsionEnergy);
4473 }
4474 if (threadID == 0) {
4475 angleTorsionTime = -System.nanoTime();
4476 }
4477 execute(0, nAngleTorsions - 1, angleTorsionLoops[threadID]);
4478 if (threadID == 0) {
4479 angleTorsionTime += System.nanoTime();
4480 }
4481 }
4482
4483 if (torsionTorsionTerm) {
4484 if (torsionTorsionLoops[threadID] == null) {
4485 torsionTorsionLoops[threadID] = new BondedTermLoop(torsionTorsions,
4486 sharedTorsionTorsionEnergy);
4487 }
4488 if (threadID == 0) {
4489 torsionTorsionTime = -System.nanoTime();
4490 }
4491 execute(0, nTorsionTorsions - 1, torsionTorsionLoops[threadID]);
4492 if (threadID == 0) {
4493 torsionTorsionTime += System.nanoTime();
4494 }
4495 }
4496
4497 if (ureyBradleyTerm) {
4498 if (ureyBradleyLoops[threadID] == null) {
4499 ureyBradleyLoops[threadID] = new BondedTermLoop(ureyBradleys, sharedUreyBradleyEnergy);
4500 }
4501 if (threadID == 0) {
4502 ureyBradleyTime = -System.nanoTime();
4503 }
4504 execute(0, nUreyBradleys - 1, ureyBradleyLoops[threadID]);
4505 if (threadID == 0) {
4506 ureyBradleyTime += System.nanoTime();
4507 }
4508 }
4509
4510
4511 if (restrainPositionTerm && nRestrainPositions > 0) {
4512 if (restrainPositionLoops[threadID] == null) {
4513 restrainPositionLoops[threadID] = new BondedTermLoop(restrainPositions, sharedRestrainPositionEnergy);
4514 }
4515 if (threadID == 0) {
4516 restrainPositionTime = -System.nanoTime();
4517 }
4518 execute(0, nRestrainPositions - 1, restrainPositionLoops[threadID]);
4519 if (threadID == 0) {
4520 restrainPositionTime += System.nanoTime();
4521 }
4522 }
4523
4524
4525 if (restrainDistanceTerm && nRestrainDistances > 0) {
4526 if (restrainDistanceLoops[threadID] == null) {
4527 restrainDistanceLoops[threadID] = new BondedTermLoop(restrainDistances, sharedRestrainDistanceEnergy);
4528 }
4529 if (threadID == 0) {
4530 restrainDistanceTime = -System.nanoTime();
4531 }
4532 execute(0, nRestrainDistances - 1, restrainDistanceLoops[threadID]);
4533 if (threadID == 0) {
4534 restrainDistanceTime += System.nanoTime();
4535 }
4536 }
4537
4538
4539 if (restrainTorsionTerm && nRestrainTorsions > 0) {
4540 if (restrainTorsionLoops[threadID] == null) {
4541 restrainTorsionLoops[threadID] = new BondedTermLoop(restrainTorsions, sharedRestrainTorsionEnergy);
4542 }
4543 if (threadID == 0) {
4544 restrainTorsionTime = -System.nanoTime();
4545 }
4546 execute(0, nRestrainTorsions - 1, restrainTorsionLoops[threadID]);
4547 if (threadID == 0) {
4548 restrainTorsionTime += System.nanoTime();
4549 }
4550 }
4551
4552
4553 if (gradient) {
4554 if (gradReduceLoops[threadID] == null) {
4555 gradReduceLoops[threadID] = new GradReduceLoop();
4556 }
4557 execute(0, nAtoms - 1, gradReduceLoops[threadID]);
4558 }
4559 }
4560
4561 public void setGradient(boolean gradient) {
4562 this.gradient = gradient;
4563 }
4564
4565 @Override
4566 public void start() {
4567
4568 sharedAngleRMSD.set(0.0);
4569 sharedBondRMSD.set(0.0);
4570
4571
4572 sharedAngleEnergy.set(0.0);
4573 sharedBondEnergy.set(0.0);
4574 sharedImproperTorsionEnergy.set(0.0);
4575 sharedOutOfPlaneBendEnergy.set(0.0);
4576 sharedPiOrbitalTorsionEnergy.set(0.0);
4577 sharedStretchBendEnergy.set(0.0);
4578 sharedTorsionEnergy.set(0.0);
4579 sharedStretchTorsionEnergy.set(0.0);
4580 sharedAngleTorsionEnergy.set(0.0);
4581 sharedTorsionTorsionEnergy.set(0.0);
4582 sharedUreyBradleyEnergy.set(0.0);
4583
4584
4585 sharedRestrainPositionEnergy.set(0.0);
4586 sharedRestrainDistanceEnergy.set(0.0);
4587 sharedRestrainTorsionEnergy.set(0.0);
4588
4589
4590 if (gradient) {
4591 grad.alloc(nAtoms);
4592 }
4593 if (lambdaTerm) {
4594 lambdaGrad.alloc(nAtoms);
4595 }
4596 }
4597
4598 private class GradInitLoop extends IntegerForLoop {
4599
4600 @Override
4601 public void run(int first, int last) {
4602 int threadID = getThreadIndex();
4603 if (gradient) {
4604 grad.reset(threadID, first, last);
4605 for (int i = first; i <= last; i++) {
4606 atoms[i].setXYZGradient(0.0, 0.0, 0.0);
4607 }
4608 }
4609 if (lambdaTerm) {
4610 lambdaGrad.reset(threadID, first, last);
4611 for (int i = first; i <= last; i++) {
4612 atoms[i].setLambdaXYZGradient(0.0, 0.0, 0.0);
4613 }
4614 }
4615 }
4616
4617 @Override
4618 public IntegerSchedule schedule() {
4619 return IntegerSchedule.fixed();
4620 }
4621 }
4622
4623 private class GradReduceLoop extends IntegerForLoop {
4624
4625 @Override
4626 public void run(int first, int last) {
4627 if (gradient) {
4628 grad.reduce(first, last);
4629 for (int i = first; i <= last; i++) {
4630 Atom a = atoms[i];
4631 a.setXYZGradient(grad.getX(i), grad.getY(i), grad.getZ(i));
4632 }
4633 }
4634 if (lambdaTerm) {
4635 lambdaGrad.reduce(first, last);
4636 for (int i = first; i <= last; i++) {
4637 Atom a = atoms[i];
4638 a.setLambdaXYZGradient(lambdaGrad.getX(i), lambdaGrad.getY(i), lambdaGrad.getZ(i));
4639 }
4640 }
4641 }
4642
4643 @Override
4644 public IntegerSchedule schedule() {
4645 return IntegerSchedule.fixed();
4646 }
4647 }
4648
4649 private class BondedTermLoop extends IntegerForLoop {
4650
4651 private final BondedTerm[] terms;
4652 private final SharedDouble sharedEnergy;
4653 private final SharedDouble sharedRMSD;
4654 private final boolean computeRMSD;
4655 private double localEnergy;
4656 private double localRMSD;
4657 private int threadID;
4658
4659 BondedTermLoop(BondedTerm[] terms, SharedDouble sharedEnergy) {
4660 this(terms, sharedEnergy, null);
4661 }
4662
4663 BondedTermLoop(BondedTerm[] terms, SharedDouble sharedEnergy, SharedDouble sharedRMSD) {
4664 this.terms = terms;
4665 this.sharedEnergy = sharedEnergy;
4666 this.sharedRMSD = sharedRMSD;
4667 computeRMSD = (sharedRMSD != null);
4668 }
4669
4670 @Override
4671 public void finish() {
4672 sharedEnergy.addAndGet(localEnergy);
4673 if (computeRMSD) {
4674 sharedRMSD.addAndGet(localRMSD);
4675 }
4676 }
4677
4678 @Override
4679 public void run(int first, int last) {
4680 for (int i = first; i <= last; i++) {
4681 BondedTerm term = terms[i];
4682
4683
4684
4685
4686
4687
4688
4689
4690
4691
4692
4693
4694
4695
4696
4697
4698 boolean used = !lambdaBondedTerms || lambdaAllBondedTerms || (term.applyLambda() && !term.isLambdaScaled());
4699 if (used) {
4700 localEnergy += term.energy(gradient, threadID, grad, lambdaGrad);
4701 if (computeRMSD) {
4702 double value = term.getValue();
4703 localRMSD += value * value;
4704 }
4705 }
4706 }
4707 }
4708
4709 @Override
4710 public void start() {
4711 localEnergy = 0.0;
4712 localRMSD = 0.0;
4713 threadID = getThreadIndex();
4714 }
4715 }
4716 }
4717 }