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.CustomExternalForce;
43 import ffx.potential.bonded.Atom;
44 import ffx.potential.nonbonded.RestrainPosition;
45
46 import java.util.List;
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(OpenMMEnergy openMMEnergy) {
67 super("k0*periodicdistance(x,y,z,x0,y0,z0)^2");
68
69 List<RestrainPosition> restrainPositionList = openMMEnergy.getRestrainPositions();
70 if (restrainPositionList == null || restrainPositionList.isEmpty()) {
71 destroy();
72 return;
73 }
74
75
76 addPerParticleParameter("k0");
77 addPerParticleParameter("x0");
78 addPerParticleParameter("y0");
79 addPerParticleParameter("z0");
80
81 int nRestraints = restrainPositionList.size();
82 double convert = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
83
84 for (RestrainPosition restrainPosition : openMMEnergy.getRestrainPositions()) {
85 double forceConstant = restrainPosition.getForceConstant() * convert;
86 Atom[] restrainPositionAtoms = restrainPosition.getAtoms();
87 int numAtoms = restrainPosition.getNumAtoms();
88 double[][] equilibriumCoordinates = restrainPosition.getEquilibriumCoordinates();
89 for (int i = 0; i < numAtoms; i++) {
90 equilibriumCoordinates[i][0] *= OpenMM_NmPerAngstrom;
91 equilibriumCoordinates[i][1] *= OpenMM_NmPerAngstrom;
92 equilibriumCoordinates[i][2] *= OpenMM_NmPerAngstrom;
93 }
94
95 DoubleArray parameters = new DoubleArray(4);
96 for (int i = 0; i < numAtoms; i++) {
97 int index = restrainPositionAtoms[i].getXyzIndex() - 1;
98 parameters.set(0, forceConstant);
99 for (int j = 0; j < 3; j++) {
100 parameters.set(j + 1, equilibriumCoordinates[i][j]);
101 }
102 addParticle(index, parameters);
103 }
104 parameters.destroy();
105 }
106
107 int forceGroup = openMMEnergy.getMolecularAssembly().getForceField().getInteger("RESTRAIN_POSITION_FORCE_GROUP", 0);
108 setForceGroup(forceGroup);
109 logger.log(Level.INFO, format(" Restrain Positions\t%6d\t\t%d", nRestraints, forceGroup));
110 }
111
112
113
114
115 public static Force constructForce(OpenMMEnergy openMMEnergy) {
116 List<RestrainPosition> restrainPositionList = openMMEnergy.getRestrainPositions();
117 if (restrainPositionList == null || restrainPositionList.isEmpty()) {
118 return null;
119 }
120 return new RestrainPositionsForce(openMMEnergy);
121 }
122
123 }