AmoebaGeneralizedKirkwoodForce.h 13.5 KB
Newer Older
Mark Friedrichs's avatar
Mark Friedrichs committed
1
2
3
4
5
6
7
8
9
10
11
#ifndef AMOEBA_OPENMM_GK_FORCE_FIELD_H_
#define AMOEBA_OPENMM_GK_FORCE_FIELD_H_

/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
12
13
 * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
 * Authors: Mark Friedrichs, Peter Eastman                                    *
Mark Friedrichs's avatar
Mark Friedrichs committed
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

#include "openmm/Force.h"
36
#include "internal/windowsExportAmoeba.h"
Mark Friedrichs's avatar
Cleanup  
Mark Friedrichs committed
37
#include <vector>
Mark Friedrichs's avatar
Mark Friedrichs committed
38
39
40
41

namespace OpenMM {

/**
42
 * This class implements an implicit solvation force using the generalized Kirkwood/Grycuk model.
Mark Friedrichs's avatar
Mark Friedrichs committed
43
 * <p>
44
45
 * To use this class, create an AmoebaGeneralizedKirkwoodForce object, then call addParticle() once for each particle in the
 * System to define its parameters.  The number of particles for which you define parameters must
Mark Friedrichs's avatar
Cleanup  
Mark Friedrichs committed
46
 * be equal to the number of particles in the System, or else an exception will be thrown when you
Mark Friedrichs's avatar
Mark Friedrichs committed
47
 * try to create a Context.  After a particle has been added, you can modify its force field parameters
48
49
 * by calling setParticleParameters().  This will have no effect on Contexts that already exist unless you
 * call updateParametersInContext().
Mark Friedrichs's avatar
Mark Friedrichs committed
50
51
 */

52
class OPENMM_EXPORT_AMOEBA AmoebaGeneralizedKirkwoodForce : public Force {
Mark Friedrichs's avatar
Cleanup  
Mark Friedrichs committed
53

Mark Friedrichs's avatar
Mark Friedrichs committed
54
55
56
public:

    /*
57
     * Create an AmoebaGeneralizedKirkwoodForce.
Mark Friedrichs's avatar
Mark Friedrichs committed
58
59
60
61
62
63
64
65
66
67
68
     */
    AmoebaGeneralizedKirkwoodForce();

    /**
     * Get the number of particles in the system.
     */
    int getNumParticles() const {
        return particles.size();
    }

    /**
69
     * Add the  parameters for a particle.  This should be called once for each particle
Mark Friedrichs's avatar
Mark Friedrichs committed
70
71
     * in the System.  When it is called for the i'th time, it specifies the parameters for the i'th particle.
     *
72
73
74
75
     * This method is provided for backwards compatibility. Compared to the alternative five parameter addParticle
     * method, the descreenRadius parameter is set to base radius value and the neckFactor is set to zero
     * (no neck descreening).
     *
Mark Friedrichs's avatar
Mark Friedrichs committed
76
     * @param charge         the charge of the particle, measured in units of the proton charge
77
78
     * @param radius         the atomic radius of the particle, measured in nm
     * @param scalingFactor  the scaling factor for the particle
Mark Friedrichs's avatar
Mark Friedrichs committed
79
80
81
82
     * @return the index of the particle that was added
     */
    int addParticle(double charge, double radius, double scalingFactor);

83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
    /**
     * Add the  parameters for a particle.  This should be called once for each particle
     * in the System.  When it is called for the i'th time, it specifies the parameters for the i'th particle.
     * <p>
     * For generalized Born / generalized Kirkwood methods, the radius of each atom has two roles. The first
     * is to define the base radius of the atom when computing its effective radius. This base radius is usually
     * parameterized against solvation free energy differences. The second role is to describe how much continuum
     * water is displaced when the atom descreens water for the calculation of the Born radii of other atoms.
     * Separation of the two roles into the "radius" and "descreenRadius" parameters gives model developers more
     * control over these separate roles.
     * <p>
     * For example, the fitting of base "radius" values will usually result in deviation from the force field's
     * van der Waals definition of Rmin (or sigma). The descreenRadius can be defined separately using force field
     * van der Waals Rmin values, which maintains consistency of atomic sizes during the HCT pairwise
     * descreening integral. The "scalingFactor" is applied to the descreenRadius during the HCT pairwise
     * descreening integral, while the neckFactor (if greater than zero) includes neck contributions to descreening.
     *
     * @param charge         the charge of the particle, measured in units of the proton charge
     * @param radius         the atomic radius of the particle, measured in nm
     * @param scalingFactor  the scaling factor for the particle (unitless)
     * @param descreenRadius the atomic radius of the particle for descreening, measure in nm
     * @param neckFactor     the scaling factor for interstitial neck descreening (unitless)
     * @return the index of the particle that was added
    */
    int addParticle(double charge, double radius, double scalingFactor, double descreenRadius, double neckFactor);

Mark Friedrichs's avatar
Mark Friedrichs committed
109
110
    /**
     * Get the force field parameters for a particle.
111
     *
Robert McGibbon's avatar
Robert McGibbon committed
112
113
114
115
     * @param      index          the index of the particle for which to get parameters
     * @param[out] charge         the charge of the particle, measured in units of the proton charge
     * @param[out] radius         the atomic radius of the particle, measured in nm
     * @param[out] scalingFactor  the scaling factor for the particle
116
117
118
119
     * @param[out] descreenRadius the atomic radius of the particle for descreening, measure in nm
     * @param[out] neckFactor     the scaling factor for interstitial neck descreening (unitless)
    */
    void getParticleParameters(int index, double& charge, double& radius, double& scalingFactor, double& descreenRadius, double& neckFactor) const;
Mark Friedrichs's avatar
Mark Friedrichs committed
120
121
122

    /**
     * Set the force field parameters for a particle.
123
124
125
126
127
128
129
130
131
    *
    * @param index          the index of the particle for which to set parameters
    * @param charge         the charge of the particle, measured in units of the proton charge
    * @param radius         the atomic radius of the particle, measured in nm
    * @param scalingFactor  the scaling factor for the particle
    * @param descreenRadius the atomic radius of the particle for descreening, measure in nm
    * @param neckFactor     the scaling factor for interstitial neck descreening (unitless)
    */
    void setParticleParameters(int index, double charge, double radius, double scalingFactor, double descreenRadius, double neckFactor);
Mark Friedrichs's avatar
Mark Friedrichs committed
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

    /**
     * Get the dielectric constant for the solvent.
     */
    double getSolventDielectric() const {
        return solventDielectric;
    }

    /**
     * Set the dielectric constant for the solvent.
     */
    void setSolventDielectric(double dielectric) {
        solventDielectric = dielectric;
    }

    /**
     * Get the dielectric constant for the solute.
     */
    double getSoluteDielectric() const {
        return soluteDielectric;
    }

    /**
     * Set the dielectric constant for the solute.
156
157
     *
     * @param dielectric The solute dielectric constant.
Mark Friedrichs's avatar
Mark Friedrichs committed
158
159
160
161
162
     */
    void setSoluteDielectric(double dielectric) {
        soluteDielectric = dielectric;
    }

163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
    /**
     * Get the flag signaling whether the solute descreening integral is rescaled by a Tanh function
     * to account for interstitial spaces. If True then tanh rescaling is used.
    */
    bool getTanhRescaling() const {
        return tanhRescaling;
    }

    /**
     * Set the flag signaling whether the solute descreening integral is rescaled by a Tanh function
     * to account for interstitial spaces.
     *
     * @param tanhRescale False to turn off Tanh rescaling; true to turn on.
    */
    void setTanhRescaling(bool tanhRescale) {
        this->tanhRescaling = tanhRescale;
    }

    /**
     * Get Tanh function parameters b0, b1 and b2.
     *
     * @param b0 The first tanh parameter.
     * @param b1 The second tanh parameter.
     * @param b2 The third tanh parameter.
     */
    void getTanhParameters(double& b0, double& b1, double& b2) const {
        b0 = this->beta0;
        b1 = this->beta1;
        b2 = this->beta2;
    }

    /**
     * Set the the Tanh function parameters to account for interstitial spaces.
     *
     * @param b0 The first tanh parameter.
     * @param b1 The second tanh parameter.
     * @param b2 The third tanh parameter.
     */
    void setTanhParameters(double b0, double b1, double b2) {
        this->beta0 = b0;
        this->beta1 = b1;
        this->beta2 = b2;
    }

    /**
    * Get the offset added to the atomic radius of each atom that sets the beginning of the
    * descreening integral when calculating effective Born radii.
    */
    double getDescreenOffset() const;

    /**
     * Get the offset added to the atomic radius of each atom that sets the beginning of the
     * descreening integral when calculating effective Born radii.
     *
     * @param descreenOffet The descreening offset (nm).
     */
    void setDescreenOffset(double descreenOffet);

Mark Friedrichs's avatar
Mark Friedrichs committed
221
    /**
222
     * Get the flag signaling whether the cavity term should be included
Mark Friedrichs's avatar
Mark Friedrichs committed
223
     */
peastman's avatar
peastman committed
224
    int getIncludeCavityTerm() const;
Mark Friedrichs's avatar
Mark Friedrichs committed
225
226

    /**
227
228
229
     * Set the flag signaling whether the cavity term should be included.
     *
     * @param includeCavityTerm Zero to turn off the cavity term; one to turn on.
Mark Friedrichs's avatar
Mark Friedrichs committed
230
231
232
233
     */
    void setIncludeCavityTerm(int includeCavityTerm);

    /**
234
     * Get the probe radius (nm) used for the cavity contribution.
Mark Friedrichs's avatar
Mark Friedrichs committed
235
236
237
238
     */
    double getProbeRadius() const;

    /**
239
240
241
     * Set the probe radius (nm) used for the cavity contribution.
     *
     * @param probeRadius The probeRadius for the cavity term.
Mark Friedrichs's avatar
Mark Friedrichs committed
242
243
244
     */
    void setProbeRadius(double probeRadius);

245
246
247
248
249
250
251
252
253
254
255
256
    /**
     * Get the dielectric offset (nm) used for cavity contribution.
     */
    double getDielectricOffset() const;

    /**
     * Set the dielectric offset (nm) used for cavity contribution.
     *
     * @param dielectricOffset The dielectric offset (nm).
     */
    void setDielectricOffset(double dielectricOffset);

Mark Friedrichs's avatar
Mark Friedrichs committed
257
258
259
260
261
262
    /**
     * Get the surface area factor kJ/(nm*nm) used in SASA contribution
     */
    double getSurfaceAreaFactor() const;

    /**
263
264
265
     * Set the surface area factor kJ/(nm*nm) used in SASA contribution.
     *
     * @param surfaceAreaFactor The surface area factor in kJ/(nm*nm).
Mark Friedrichs's avatar
Mark Friedrichs committed
266
     */
267
    void setSurfaceAreaFactor(double surfaceAreaFactor);
268
269
270
    /**
     * Update the per-particle parameters in a Context to match those stored in this Force object.  This method provides
     * an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
peastman's avatar
peastman committed
271
     * Simply call setParticleParameters() to modify this object's parameters, then call updateParametersInContext()
272
273
274
275
276
277
     * to copy them over to the Context.
     * 
     * The only information this method updates is the values of per-particle parameters.  All other aspects of the Force
     * (the probe radius, the surface area factor, etc.) are unaffected and can only be changed by reinitializing the Context.
     */
    void updateParametersInContext(Context& context);
278
279
280
281
282
283
    /**
     * Returns whether or not this force makes use of periodic boundary
     * conditions.
     *
     * @returns true if nonbondedMethod uses PBC and false otherwise
     */
284
285
286
    bool usesPeriodicBoundaryConditions() const {
        return false;
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
287
protected:
288
    ForceImpl* createImpl() const;
Mark Friedrichs's avatar
Mark Friedrichs committed
289
290
291
private:
    class ParticleInfo;
    int includeCavityTerm;
292
    bool tanhRescaling;
293
    double solventDielectric, soluteDielectric, dielectricOffset,
Mark Friedrichs's avatar
Mark Friedrichs committed
294
           probeRadius, surfaceAreaFactor;
295
    double beta0, beta1, beta2, descreenOffset;
Mark Friedrichs's avatar
Mark Friedrichs committed
296
297
298
    std::vector<ParticleInfo> particles;
};

299
300
301
302
/**
 * This is an internal class used to record information about a particle.
 * @private
 */
Mark Friedrichs's avatar
Mark Friedrichs committed
303
304
class AmoebaGeneralizedKirkwoodForce::ParticleInfo {
public:
305
    double charge, radius, scalingFactor, descreenRadius, neckFactor;
Mark Friedrichs's avatar
Mark Friedrichs committed
306
    ParticleInfo() {
307
        charge = radius = scalingFactor = descreenRadius = neckFactor = 0.0;
Mark Friedrichs's avatar
Mark Friedrichs committed
308
    }
309
310
    ParticleInfo(double charge, double radius, double scalingFactor, double descreenRadius, double neckFactor) :
        charge(charge), radius(radius), scalingFactor(scalingFactor), descreenRadius(descreenRadius), neckFactor(neckFactor) {
Mark Friedrichs's avatar
Mark Friedrichs committed
311
312
313
314
315
    }
};

} // namespace OpenMM

316
#endif /*AMOEBA_OPENMM_GK_FORCE_FIELD_H_*/