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