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.openmm;
39
40 import ffx.openmm.CustomBondForce;
41 import ffx.openmm.DoubleArray;
42 import ffx.openmm.Force;
43 import ffx.potential.bonded.Atom;
44 import ffx.potential.bonded.RestrainDistance;
45 import ffx.potential.parameters.BondType;
46 import ffx.potential.terms.RestrainDistancePotentialEnergy;
47
48 import java.util.List;
49 import java.util.logging.Level;
50 import java.util.logging.Logger;
51
52 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
53 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
54 import static java.lang.String.format;
55
56
57
58
59 public class RestrainDistanceForce extends CustomBondForce {
60
61 private static final Logger logger = Logger.getLogger(RestrainDistanceForce.class.getName());
62
63
64
65
66
67
68
69 public RestrainDistanceForce(BondType.BondFunction bondFunction,
70 RestrainDistancePotentialEnergy restrainDistancePotentialEnergy) {
71 super(bondFunction.toMathematicalForm());
72 List<RestrainDistance> restrainDistances = restrainDistancePotentialEnergy.getRestrainDistances(bondFunction);
73 if (restrainDistances == null || restrainDistances.isEmpty()) {
74 destroy();
75 return;
76 }
77
78 addPerBondParameter("k");
79 addPerBondParameter("r0");
80 if (bondFunction.hasFlatBottom()) {
81 addPerBondParameter("fb");
82 }
83
84 BondType bondType = restrainDistances.getFirst().bondType;
85 switch (bondFunction) {
86 case QUARTIC, FLAT_BOTTOM_QUARTIC -> {
87 addGlobalParameter("cubic", bondType.cubic / OpenMM_NmPerAngstrom);
88 addGlobalParameter("quartic", bondType.quartic / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom));
89 }
90 }
91
92
93 double forceConvert = 2.0 * OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
94 DoubleArray parameters = new DoubleArray(0);
95 for (RestrainDistance restrainDistance : restrainDistances) {
96 bondType = restrainDistance.bondType;
97 double forceConstant = bondType.forceConstant * bondType.bondUnit * forceConvert;
98 double distance = bondType.distance * OpenMM_NmPerAngstrom;
99 Atom[] atoms = restrainDistance.getAtomArray();
100 int i1 = atoms[0].getXyzIndex() - 1;
101 int i2 = atoms[1].getXyzIndex() - 1;
102 parameters.append(forceConstant);
103 parameters.append(distance);
104 if (bondFunction.hasFlatBottom()) {
105 parameters.append(bondType.flatBottomRadius * OpenMM_NmPerAngstrom);
106 }
107 addBond(i1, i2, parameters);
108 parameters.destroy();
109 }
110 parameters.destroy();
111
112 int forceGroup = restrainDistancePotentialEnergy.getForceGroup();
113 setForceGroup(forceGroup);
114 logger.log(Level.INFO, format(" Restrain-Distance force \t%6d\t%d", restrainDistances.size(), forceGroup));
115 }
116
117
118
119
120
121
122
123 public static Force constructForce(BondType.BondFunction bondFunction, OpenMMEnergy openMMEnergy) {
124 RestrainDistancePotentialEnergy restrainDistancePotentialEnergy =
125 openMMEnergy.getRestrainDistancePotentialEnergy();
126 if (restrainDistancePotentialEnergy == null) {
127 return null;
128 }
129 List<RestrainDistance> restrainDistances = restrainDistancePotentialEnergy.getRestrainDistances(bondFunction);
130 if (restrainDistances == null || restrainDistances.isEmpty()) {
131 return null;
132 }
133 return new RestrainDistanceForce(bondFunction, restrainDistancePotentialEnergy);
134 }
135 }