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.Force;
41 import ffx.openmm.amoeba.GKCavitationForce;
42 import ffx.potential.bonded.Atom;
43 import ffx.potential.nonbonded.GeneralizedKirkwood;
44 import ffx.potential.nonbonded.ParticleMeshEwald;
45 import ffx.potential.nonbonded.implicit.ChandlerCavitation;
46 import ffx.potential.nonbonded.implicit.DispersionRegion;
47 import ffx.potential.nonbonded.implicit.GaussVol;
48
49 import java.util.logging.Level;
50 import java.util.logging.Logger;
51
52
53 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
54 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
55 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_False;
56 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_True;
57 import static java.lang.String.format;
58
59
60
61
62 public class AmoebaGKCavitationForce extends GKCavitationForce {
63
64 private static final Logger logger = Logger.getLogger(AmoebaGKCavitationForce.class.getName());
65
66
67
68
69
70
71 public AmoebaGKCavitationForce(OpenMMEnergy openMMEnergy) {
72 logger.severe(" The AmoebaGKCavitationForce is not currently supported.");
73
74
75 GeneralizedKirkwood generalizedKirkwood = openMMEnergy.getGK();
76 if (generalizedKirkwood == null) {
77 destroy();
78 return;
79 }
80 ChandlerCavitation chandlerCavitation = generalizedKirkwood.getChandlerCavitation();
81 if (chandlerCavitation == null) {
82 destroy();
83 return;
84 }
85 GaussVol gaussVol = chandlerCavitation.getGaussVol();
86 if (gaussVol == null) {
87 destroy();
88 return;
89 }
90
91 double surfaceTension = chandlerCavitation.getSurfaceTension()
92 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom / OpenMM_NmPerAngstrom;
93 double[] rad = gaussVol.getRadii();
94
95 int index = 0;
96 Atom[] atoms = openMMEnergy.getMolecularAssembly().getAtomArray();
97 for (Atom atom : atoms) {
98 int isHydrogen = OpenMM_False;
99 double radius = rad[index++];
100 if (atom.isHydrogen()) {
101 isHydrogen = OpenMM_True;
102 radius = 0.0;
103 }
104 addParticle(radius * OpenMM_NmPerAngstrom, surfaceTension, isHydrogen);
105 }
106
107
108
109
110 int forceGroup = openMMEnergy.getMolecularAssembly().getForceField().getInteger("GK_FORCE_GROUP", 2);
111 setForceGroup(forceGroup);
112 logger.log(Level.INFO, format(" GaussVol cavitation force \t\t%d", forceGroup));
113 }
114
115
116
117
118
119
120
121 public static Force constructForce(OpenMMEnergy openMMEnergy) {
122 GeneralizedKirkwood gk = openMMEnergy.getGK();
123 if (gk == null) {
124 return null;
125 }
126 DispersionRegion dispersionRegion = gk.getDispersionRegion();
127 if (dispersionRegion == null) {
128 return null;
129 }
130 return new AmoebaGKCavitationForce(openMMEnergy);
131 }
132
133
134
135
136
137
138
139 public void updateForce(Atom[] atoms, OpenMMEnergy openMMEnergy) {
140 GeneralizedKirkwood generalizedKirkwood = openMMEnergy.getGK();
141 if (generalizedKirkwood == null) {
142 return;
143 }
144 ChandlerCavitation chandlerCavitation = generalizedKirkwood.getChandlerCavitation();
145 if (chandlerCavitation == null) {
146 return;
147 }
148 GaussVol gaussVol = chandlerCavitation.getGaussVol();
149 if (gaussVol == null) {
150 return;
151 }
152
153 double surfaceTension = chandlerCavitation.getSurfaceTension()
154 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom / OpenMM_NmPerAngstrom;
155
156 ParticleMeshEwald pme = openMMEnergy.getPmeNode();
157 double lambdaElec = pme.getAlchemicalParameters().permLambda;
158
159
160
161
162
163 double[] rad = gaussVol.getRadii();
164
165 for (Atom atom : atoms) {
166 int index = atom.getXyzIndex() - 1;
167 double useFactor = 1.0;
168 if (!atom.getUse()) {
169 useFactor = 0.0;
170 }
171
172
173 double lambdaScale = lambdaElec * lambdaElec;
174 if (!atom.applyLambda()) {
175 lambdaScale = 1.0;
176 }
177 useFactor *= lambdaScale;
178
179 double radius = rad[index];
180 int isHydrogen = OpenMM_False;
181 if (atom.isHydrogen()) {
182 isHydrogen = OpenMM_True;
183 radius = 0.0;
184 }
185
186 setParticleParameters(index, radius * OpenMM_NmPerAngstrom, surfaceTension * useFactor, isHydrogen);
187 }
188 updateParametersInContext(openMMEnergy.getContext());
189 }
190 }