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