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.WcaDispersionForce;
42 import ffx.potential.bonded.Atom;
43 import ffx.potential.nonbonded.GeneralizedKirkwood;
44 import ffx.potential.nonbonded.VanDerWaals;
45 import ffx.potential.nonbonded.VanDerWaalsForm;
46 import ffx.potential.nonbonded.implicit.DispersionRegion;
47 import ffx.potential.parameters.VDWType;
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 java.lang.String.format;
55
56 public class AmoebaWcaDispersionForce extends WcaDispersionForce {
57
58 private static final Logger logger = Logger.getLogger(AmoebaGeneralizedKirkwoodForce.class.getName());
59
60
61
62
63
64
65 public AmoebaWcaDispersionForce(OpenMMEnergy openMMEnergy) {
66 GeneralizedKirkwood gk = openMMEnergy.getGK();
67 if (gk == null) {
68 destroy();
69 return;
70 }
71 DispersionRegion dispersionRegion = gk.getDispersionRegion();
72 if (dispersionRegion == null) {
73 return;
74 }
75
76 double epso = 0.1100;
77 double epsh = 0.0135;
78 double rmino = 1.7025;
79 double rminh = 1.3275;
80 double awater = 0.033428;
81 double slevy = 1.0;
82 double dispoff = dispersionRegion.getDispersionOffset();
83 double shctd = dispersionRegion.getDispersionOverlapFactor();
84
85 VanDerWaals vdW = openMMEnergy.getVdwNode();
86 VanDerWaalsForm vdwForm = vdW.getVDWForm();
87 double radScale = 1.0;
88 if (vdwForm.radiusSize == VanDerWaalsForm.RADIUS_SIZE.DIAMETER) {
89 radScale = 0.5;
90 }
91
92 Atom[] atoms = openMMEnergy.getMolecularAssembly().getAtomArray();
93 for (Atom atom : atoms) {
94 VDWType vdwType = atom.getVDWType();
95 double radius = vdwType.radius;
96 double eps = vdwType.wellDepth;
97 addParticle(OpenMM_NmPerAngstrom * radius * radScale, OpenMM_KJPerKcal * eps);
98 }
99
100 setEpso(epso * OpenMM_KJPerKcal);
101 setEpsh(epsh * OpenMM_KJPerKcal);
102 setRmino(rmino * OpenMM_NmPerAngstrom);
103 setRminh(rminh * OpenMM_NmPerAngstrom);
104 setDispoff(dispoff * OpenMM_NmPerAngstrom);
105 setAwater(awater / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom));
106 setSlevy(slevy);
107 setShctd(shctd);
108
109 int forceGroup = openMMEnergy.getMolecularAssembly().getForceField().getInteger("GK_FORCE_GROUP", 1);
110 setForceGroup(forceGroup);
111 logger.log(Level.INFO, format(" WCA dispersion force \t\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 AmoebaWcaDispersionForce(openMMEnergy);
130 }
131
132
133
134
135
136
137
138 public void updateForce(Atom[] atoms, OpenMMEnergy openMMEnergy) {
139 VanDerWaals vdW = openMMEnergy.getVdwNode();
140 VanDerWaalsForm vdwForm = vdW.getVDWForm();
141 double radScale = 1.0;
142 if (vdwForm.radiusSize == VanDerWaalsForm.RADIUS_SIZE.DIAMETER) {
143 radScale = 0.5;
144 }
145
146 double lambdaElec = openMMEnergy.getSystem().getLambdaElec();
147
148 for (Atom atom : atoms) {
149 int index = atom.getXyzIndex() - 1;
150 double useFactor = 1.0;
151 if (!atom.getUse()) {
152 useFactor = 0.0;
153 }
154
155
156
157 double lambdaScale = lambdaElec * lambdaElec;
158 if (!atom.applyLambda()) {
159 lambdaScale = 1.0;
160 }
161 useFactor *= lambdaScale;
162
163 VDWType vdwType = atom.getVDWType();
164 double radius = vdwType.radius;
165 double eps = vdwType.wellDepth;
166 setParticleParameters(index, OpenMM_NmPerAngstrom * radius * radScale, OpenMM_KJPerKcal * eps * useFactor);
167 }
168
169 updateParametersInContext(openMMEnergy.getContext());
170 }
171 }