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.implicit.ChandlerCavitation;
45 import ffx.potential.nonbonded.implicit.DispersionRegion;
46 import ffx.potential.nonbonded.implicit.GaussVol;
47
48 import java.util.logging.Level;
49 import java.util.logging.Logger;
50
51 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGKCavitationForce_NonbondedMethod.OpenMM_AmoebaGKCavitationForce_NoCutoff;
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 GeneralizedKirkwood generalizedKirkwood = openMMEnergy.getGK();
72 if (generalizedKirkwood == null) {
73 destroy();
74 return;
75 }
76 ChandlerCavitation chandlerCavitation = generalizedKirkwood.getChandlerCavitation();
77 if (chandlerCavitation == null) {
78 destroy();
79 return;
80 }
81 GaussVol gaussVol = chandlerCavitation.getGaussVol();
82 if (gaussVol == null) {
83 destroy();
84 return;
85 }
86
87 double surfaceTension = chandlerCavitation.getSurfaceTension()
88 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom / OpenMM_NmPerAngstrom;
89 double[] rad = gaussVol.getRadii();
90
91 int index = 0;
92 Atom[] atoms = openMMEnergy.getMolecularAssembly().getAtomArray();
93 for (Atom atom : atoms) {
94 int isHydrogen = OpenMM_False;
95 double radius = rad[index++];
96 if (atom.isHydrogen()) {
97 isHydrogen = OpenMM_True;
98 radius = 0.0;
99 }
100 addParticle(radius * OpenMM_NmPerAngstrom, surfaceTension, isHydrogen);
101 }
102
103 setNonbondedMethod(OpenMM_AmoebaGKCavitationForce_NoCutoff);
104
105 int forceGroup = openMMEnergy.getMolecularAssembly().getForceField().getInteger("GK_FORCE_GROUP", 2);
106 setForceGroup(forceGroup);
107 logger.log(Level.INFO, format(" GaussVol cavitation force \t\t%d", forceGroup));
108 }
109
110
111
112
113
114
115
116 public static Force constructForce(OpenMMEnergy openMMEnergy) {
117 GeneralizedKirkwood gk = openMMEnergy.getGK();
118 if (gk == null) {
119 return null;
120 }
121 DispersionRegion dispersionRegion = gk.getDispersionRegion();
122 if (dispersionRegion == null) {
123 return null;
124 }
125 return new AmoebaGKCavitationForce(openMMEnergy);
126 }
127
128
129
130
131
132
133
134 public void updateForce(Atom[] atoms, OpenMMEnergy openMMEnergy) {
135 GeneralizedKirkwood generalizedKirkwood = openMMEnergy.getGK();
136 if (generalizedKirkwood == null) {
137 return;
138 }
139 ChandlerCavitation chandlerCavitation = generalizedKirkwood.getChandlerCavitation();
140 if (chandlerCavitation == null) {
141 return;
142 }
143 GaussVol gaussVol = chandlerCavitation.getGaussVol();
144 if (gaussVol == null) {
145 return;
146 }
147
148 double surfaceTension = chandlerCavitation.getSurfaceTension()
149 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom / OpenMM_NmPerAngstrom;
150 double lambdaElec = openMMEnergy.getSystem().getLambdaElec();
151
152
153
154
155
156 double[] rad = gaussVol.getRadii();
157
158 for (Atom atom : atoms) {
159 int index = atom.getXyzIndex() - 1;
160 double useFactor = 1.0;
161 if (!atom.getUse()) {
162 useFactor = 0.0;
163 }
164
165
166 double lambdaScale = lambdaElec * lambdaElec;
167 if (!atom.applyLambda()) {
168 lambdaScale = 1.0;
169 }
170 useFactor *= lambdaScale;
171
172 double radius = rad[index];
173 int isHydrogen = OpenMM_False;
174 if (atom.isHydrogen()) {
175 isHydrogen = OpenMM_True;
176 radius = 0.0;
177 }
178
179 setParticleParameters(index, radius * OpenMM_NmPerAngstrom, surfaceTension * useFactor, isHydrogen);
180 }
181 updateParametersInContext(openMMEnergy.getContext());
182 }
183 }