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.IntArray;
43 import ffx.openmm.CentroidBondForce;
44 import ffx.potential.bonded.Atom;
45 import ffx.potential.nonbonded.RestrainGroups;
46
47 import java.util.logging.Level;
48 import java.util.logging.Logger;
49
50 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
51 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
52 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_False;
53 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_True;
54 import static java.lang.String.format;
55
56
57
58
59 public class RestrainGroupsForce extends CentroidBondForce {
60
61 private static final Logger logger = Logger.getLogger(RestrainGroupsForce.class.getName());
62
63 private static final String energy = "step(distance(g1,g2)-u)*k*(distance(g1,g2)-u)^2+step(l-distance(g1,g2))*k*(distance(g1,g2)-l)^2";
64
65
66
67
68
69
70 public RestrainGroupsForce(OpenMMEnergy openMMEnergy) {
71 super(2, energy);
72 RestrainGroups restrainGroups = openMMEnergy.getRestrainGroups();
73 if (restrainGroups == null) {
74 destroy();
75 return;
76 }
77
78 addPerBondParameter("k");
79 addPerBondParameter("l");
80 addPerBondParameter("u");
81
82 Atom[] atoms = openMMEnergy.getMolecularAssembly().getAtomArray();
83
84
85 int nGroups = restrainGroups.getNumberOfGroups();
86 IntArray group = new IntArray(0);
87 DoubleArray weight = new DoubleArray(0);
88 for (int j = 0; j < nGroups; j++) {
89 int[] groupMembers = restrainGroups.getGroupMembers(j);
90 for (int i : groupMembers) {
91 group.append(i);
92 weight.append(atoms[i].getMass());
93 }
94 addGroup(group, weight);
95 group.resize(0);
96 weight.resize(0);
97 }
98 group.destroy();
99 weight.destroy();
100
101
102 double convert = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
103 int nRestraints = restrainGroups.getNumberOfRestraints();
104 int[] group1 = restrainGroups.getGroup1();
105 int[] group2 = restrainGroups.getGroup2();
106 double[] forceConstants = restrainGroups.getForceConstants();
107 double[] smallerDistance = restrainGroups.getSmallerDistance();
108 double[] largerDistance = restrainGroups.getLargerDistance();
109 group = new IntArray(0);
110 DoubleArray parameters = new DoubleArray(0);
111 for (int i = 0; i < nRestraints; i++) {
112 group.append(group1[i]);
113 group.append(group2[i]);
114 parameters.append(forceConstants[i] * convert);
115 parameters.append(smallerDistance[i] * OpenMM_NmPerAngstrom);
116 parameters.append(largerDistance[i] * OpenMM_NmPerAngstrom);
117 addBond(group, parameters);
118 group.resize(0);
119 parameters.resize(0);
120 }
121 group.destroy();
122 parameters.destroy();
123
124
125 int forceGroup = openMMEnergy.getMolecularAssembly().getForceField().getInteger("RESTRAIN_GROUPS_FORCE_GROUP", 0);
126 setForceGroup(forceGroup);
127
128 if (openMMEnergy.getCrystal().aperiodic()) {
129 setUsesPeriodicBoundaryConditions(OpenMM_False);
130 } else {
131 setUsesPeriodicBoundaryConditions(OpenMM_True);
132 }
133 logger.log(Level.INFO, format(" Restrain Groups \t%6d\t\t%1d", nRestraints, forceGroup));
134 }
135
136
137
138
139 public static Force constructForce(OpenMMEnergy openMMEnergy) {
140 RestrainGroups restrainGroups = openMMEnergy.getRestrainGroups();
141 if (restrainGroups == null) {
142 return null;
143 }
144 return new RestrainGroupsForce(openMMEnergy);
145 }
146 }