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.CustomExternalForce;
41 import ffx.openmm.DoubleArray;
42 import ffx.openmm.Force;
43 import ffx.potential.bonded.Atom;
44 import ffx.potential.bonded.RestrainPosition;
45 import ffx.potential.terms.RestrainPositionPotentialEnergy;
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 java.lang.String.format;
53
54
55
56
57 public class RestrainPositionsForce extends CustomExternalForce {
58
59 private static final Logger logger = Logger.getLogger(RestrainPositionsForce.class.getName());
60
61
62
63
64
65
66 public RestrainPositionsForce(RestrainPositionPotentialEnergy restrainPositionPotentialEnergy) {
67 super(RestrainPositionPotentialEnergy.getRestrainPositionEnergyString());
68 RestrainPosition[] restrainPositions = restrainPositionPotentialEnergy.getRestrainPositionArray();
69
70 addPerParticleParameter("k0");
71 addPerParticleParameter("x0");
72 addPerParticleParameter("y0");
73 addPerParticleParameter("z0");
74
75 int nRestraints = restrainPositions.length;
76 double convert = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
77
78 for (RestrainPosition restrainPosition : restrainPositions) {
79 double forceConstant = restrainPosition.getForceConstant() * convert;
80 Atom[] restrainPositionAtoms = restrainPosition.getAtoms();
81 int numAtoms = restrainPosition.getNumAtoms();
82 double[][] equilibriumCoordinates = restrainPosition.getEquilibriumCoordinates();
83 for (int i = 0; i < numAtoms; i++) {
84 equilibriumCoordinates[i][0] *= OpenMM_NmPerAngstrom;
85 equilibriumCoordinates[i][1] *= OpenMM_NmPerAngstrom;
86 equilibriumCoordinates[i][2] *= OpenMM_NmPerAngstrom;
87 }
88
89 DoubleArray parameters = new DoubleArray(4);
90 for (int i = 0; i < numAtoms; i++) {
91 int index = restrainPositionAtoms[i].getXyzIndex() - 1;
92 parameters.set(0, forceConstant);
93 for (int j = 0; j < 3; j++) {
94 parameters.set(j + 1, equilibriumCoordinates[i][j]);
95 }
96 addParticle(index, parameters);
97 }
98 parameters.destroy();
99 }
100
101 int forceGroup = restrainPositionPotentialEnergy.getForceGroup();
102 setForceGroup(forceGroup);
103 logger.log(Level.INFO, format(" Restrain Positions\t%6d\t\t%d", nRestraints, forceGroup));
104 }
105
106
107
108
109 public static Force constructForce(OpenMMEnergy openMMEnergy) {
110 RestrainPositionPotentialEnergy restrainPositionPotentialEnergy =
111 openMMEnergy.getRestrainPositionPotentialEnergy();
112 if (restrainPositionPotentialEnergy == null) {
113 return null;
114 }
115 return new RestrainPositionsForce(restrainPositionPotentialEnergy);
116 }
117
118 }