CpuGBVI.cpp 38.8 KB
Newer Older
Mark Friedrichs's avatar
Mark Friedrichs committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
 * Contributors: Pande Group
 *
 * 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 <string.h>
Mark Friedrichs's avatar
Mark Friedrichs committed
25
#include <math.h>
Mark Friedrichs's avatar
Mark Friedrichs committed
26
#include <sstream>
Mark Friedrichs's avatar
Mark Friedrichs committed
27
#include <stdio.h>
Mark Friedrichs's avatar
Mark Friedrichs committed
28

29
30
#include "SimTKOpenMMCommon.h"
#include "ReferenceForce.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
31
#include "CpuGBVI.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
32

33
34
35
using namespace std;
using namespace OpenMM;

Mark Friedrichs's avatar
Mark Friedrichs committed
36
37
/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
38
    CpuGBVI constructor
Mark Friedrichs's avatar
Mark Friedrichs committed
39

Mark Friedrichs's avatar
Mark Friedrichs committed
40
41
42
    gbviParameters      gbviParameters object
    
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
43

44
45
CpuGBVI::CpuGBVI(GBVIParameters* gbviParameters) : _gbviParameters(gbviParameters) {
    _switchDeriviative.resize(gbviParameters->getNumberOfAtoms());
Mark Friedrichs's avatar
Mark Friedrichs committed
46
47
48
49
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
50
    CpuGBVI destructor
Mark Friedrichs's avatar
Mark Friedrichs committed
51

Mark Friedrichs's avatar
Mark Friedrichs committed
52
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
53

54
CpuGBVI::~CpuGBVI() {
Mark Friedrichs's avatar
Mark Friedrichs committed
55
56
57
58
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
59
    Get GBVIParameters reference
Mark Friedrichs's avatar
Mark Friedrichs committed
60

Mark Friedrichs's avatar
Mark Friedrichs committed
61
    @return GBVIParameters reference
Mark Friedrichs's avatar
Mark Friedrichs committed
62

Mark Friedrichs's avatar
Mark Friedrichs committed
63
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
64

65
GBVIParameters* CpuGBVI::getGBVIParameters() const {
Mark Friedrichs's avatar
Mark Friedrichs committed
66
    return _gbviParameters;
Mark Friedrichs's avatar
Mark Friedrichs committed
67
68
69
70
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
71
    Set GBVIParameters reference
Mark Friedrichs's avatar
Mark Friedrichs committed
72

Mark Friedrichs's avatar
Mark Friedrichs committed
73
    @param GBVIParameters reference
Mark Friedrichs's avatar
Mark Friedrichs committed
74

Mark Friedrichs's avatar
Mark Friedrichs committed
75
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
76

77
void CpuGBVI::setGBVIParameters(GBVIParameters* gbviParameters) {
Mark Friedrichs's avatar
Mark Friedrichs committed
78
    _gbviParameters = gbviParameters;
79
80
81
82
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
83
    Return OBC chain derivative: size = _obcParameters->getNumberOfAtoms()
84

Mark Friedrichs's avatar
Mark Friedrichs committed
85
    @return array
86

Mark Friedrichs's avatar
Mark Friedrichs committed
87
    --------------------------------------------------------------------------------------- */
88

89
RealOpenMMVector& CpuGBVI::getSwitchDeriviative() {
Mark Friedrichs's avatar
Mark Friedrichs committed
90
    return _switchDeriviative;
91
92
93
94
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
95
    Compute quintic spline value and associated derviative
96

Mark Friedrichs's avatar
Mark Friedrichs committed
97
98
99
100
101
    @param x                   value to compute spline at
    @param rl                  lower cutoff value
    @param ru                  upper cutoff value
    @param outValue            value of spline at x
    @param outDerivative       value of derivative of spline at x
102

Mark Friedrichs's avatar
Mark Friedrichs committed
103
    --------------------------------------------------------------------------------------- */
104

105
106
void CpuGBVI::quinticSpline(RealOpenMM x, RealOpenMM rl, RealOpenMM ru,
                            RealOpenMM* outValue, RealOpenMM* outDerivative) {
107

Mark Friedrichs's avatar
Mark Friedrichs committed
108
    // ---------------------------------------------------------------------------------------
109

110
111
112
113
114
115
    static const RealOpenMM one           = static_cast<RealOpenMM>(  1.0);
    static const RealOpenMM minusSix      = static_cast<RealOpenMM>( -6.0);
    static const RealOpenMM minusTen      = static_cast<RealOpenMM>(-10.0);
    static const RealOpenMM minusThirty   = static_cast<RealOpenMM>(-30.0);
    static const RealOpenMM fifteen       = static_cast<RealOpenMM>( 15.0);
    static const RealOpenMM sixty         = static_cast<RealOpenMM>( 60.0);
116

Mark Friedrichs's avatar
Mark Friedrichs committed
117
    // ---------------------------------------------------------------------------------------
118

Mark Friedrichs's avatar
Mark Friedrichs committed
119
120
121
122
123
    RealOpenMM numerator    = x  - rl;
    RealOpenMM denominator  = ru - rl;
    RealOpenMM ratio        = numerator/denominator;
    RealOpenMM ratio2       = ratio*ratio;
    RealOpenMM ratio3       = ratio2*ratio;
124

Mark Friedrichs's avatar
Mark Friedrichs committed
125
126
    *outValue               = one + ratio3*(minusTen + fifteen*ratio + minusSix*ratio2);
    *outDerivative          = ratio2*(minusThirty + sixty*ratio + minusThirty*ratio2)/denominator;
127
128
129
130
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
131
132
    Compute Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
    and quintic splice switching function
133

Mark Friedrichs's avatar
Mark Friedrichs committed
134
135
136
137
138
139
    @param atomicRadius3       atomic radius cubed
    @param bornSum             Born sum (volume integral)
    @param gbviParameters      Gbvi parameters (parameters used in spline
                               QuinticLowerLimitFactor & QuinticUpperBornRadiusLimit)
    @param bornRadius          output Born radius
    @param switchDeriviative   output switching function deriviative
140

Mark Friedrichs's avatar
Mark Friedrichs committed
141
    --------------------------------------------------------------------------------------- */
142

143
144
145
void CpuGBVI::computeBornRadiiUsingQuinticSpline(RealOpenMM atomicRadius3, RealOpenMM bornSum,
                                                 GBVIParameters* gbviParameters, 
                                                 RealOpenMM* bornRadius, RealOpenMM* switchDeriviative) {
146

Mark Friedrichs's avatar
Mark Friedrichs committed
147
    // ---------------------------------------------------------------------------------------
148

149
150
151
152
153
154
155
    static const RealOpenMM zero          = static_cast<RealOpenMM>( 0.0);
    static const RealOpenMM one           = static_cast<RealOpenMM>( 1.0);
    static const RealOpenMM minusOne      = static_cast<RealOpenMM>(-1.0);
    static const RealOpenMM minusThree    = static_cast<RealOpenMM>(-3.0);
    static const RealOpenMM oneEighth     = static_cast<RealOpenMM>( 0.125);
    static const RealOpenMM minusOneThird = static_cast<RealOpenMM>((-1.0/3.0));
    static const RealOpenMM three         = static_cast<RealOpenMM>( 3.0);
156

Mark Friedrichs's avatar
Mark Friedrichs committed
157
    // ---------------------------------------------------------------------------------------
158

Mark Friedrichs's avatar
Mark Friedrichs committed
159
    // R                = [ S(V)*(A - V) ]**(-1/3)
160

Mark Friedrichs's avatar
Mark Friedrichs committed
161
162
163
    // S(V)             = 1                                 V < L
    // S(V)             = qSpline + U/(A-V)                 L < V < A
    // S(V)             = U/(A-V)                           U < V 
164

Mark Friedrichs's avatar
Mark Friedrichs committed
165
    // dR/dr            = (-1/3)*[ S(V)*(A - V) ]**(-4/3)*[ d{ S(V)*(A-V) }/dr
166

Mark Friedrichs's avatar
Mark Friedrichs committed
167
    // d{ S(V)*(A-V) }/dr   = (dV/dr)*[ (A-V)*dS/dV - S(V) ]
168

Mark Friedrichs's avatar
Mark Friedrichs committed
169
    //  (A - V)*dS/dV - S(V)  = 0 - 1                             V < L
170

Mark Friedrichs's avatar
Mark Friedrichs committed
171
    //  (A - V)*dS/dV - S(V)  = (A-V)*d(qSpline) + (A-V)*U/(A-V)**2 - qSpline - U/(A-V) 
172

173
    //                        = (A-V)*d(qSpline) - qSpline        L < V < A**(-3)
174

Mark Friedrichs's avatar
Mark Friedrichs committed
175
176
177
178
    //  (A - V)*dS/dV - S(V)  = (A-V)*U*/(A-V)**2 - U/(A-V) = 0   U < V

    RealOpenMM splineL          = gbviParameters->getQuinticLowerLimitFactor()*atomicRadius3;
    RealOpenMM sum;
179
180
    if (bornSum > splineL) {
        if (bornSum < atomicRadius3) {
Mark Friedrichs's avatar
Mark Friedrichs committed
181
            RealOpenMM splineValue, splineDerivative;
182
            quinticSpline(bornSum, splineL, atomicRadius3, &splineValue, &splineDerivative); 
Mark Friedrichs's avatar
Mark Friedrichs committed
183
184
185
186
187
188
189
190
191
192
            sum                 = (atomicRadius3 - bornSum)*splineValue + gbviParameters->getQuinticUpperBornRadiusLimit();
            *switchDeriviative  = splineValue - (atomicRadius3 - bornSum)*splineDerivative;
        } else {   
            sum                = gbviParameters->getQuinticUpperBornRadiusLimit();
            *switchDeriviative = zero;
        }
    } else {
        sum                = atomicRadius3 - bornSum; 
        *switchDeriviative = one;
    }
193
    *bornRadius = POW(sum, minusOneThird);
Mark Friedrichs's avatar
Mark Friedrichs committed
194
195
196
197
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
198
    Get Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
Mark Friedrichs's avatar
Mark Friedrichs committed
199

Mark Friedrichs's avatar
Mark Friedrichs committed
200
201
202
    @param atomCoordinates     atomic coordinates
    @param bornRadii           output array of Born radii
    @param chain               not used here
Mark Friedrichs's avatar
Mark Friedrichs committed
203

Mark Friedrichs's avatar
Mark Friedrichs committed
204
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
205

206
void CpuGBVI::computeBornRadii(const vector<RealVec>& atomCoordinates, RealOpenMMVector& bornRadii) {
Mark Friedrichs's avatar
Mark Friedrichs committed
207

Mark Friedrichs's avatar
Mark Friedrichs committed
208
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
209

210
211
212
213
214
215
    static const RealOpenMM zero          = static_cast<RealOpenMM>(0.0);
    static const RealOpenMM one           = static_cast<RealOpenMM>(1.0);
    static const RealOpenMM minusThree    = static_cast<RealOpenMM>(-3.0);
    static const RealOpenMM oneEighth     = static_cast<RealOpenMM>(0.125);
    static const RealOpenMM minusOneThird = static_cast<RealOpenMM>((-1.0/3.0));
    static const RealOpenMM three         = static_cast<RealOpenMM>(3.0);
Mark Friedrichs's avatar
Mark Friedrichs committed
216

Mark Friedrichs's avatar
Mark Friedrichs committed
217
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
218

Mark Friedrichs's avatar
Mark Friedrichs committed
219
220
221
222
    GBVIParameters* gbviParameters        = getGBVIParameters();
    int numberOfAtoms                     = gbviParameters->getNumberOfAtoms();
    const RealOpenMMVector& atomicRadii   = gbviParameters->getAtomicRadii();
    const RealOpenMMVector& scaledRadii   = gbviParameters->getScaledRadii();
Mark Friedrichs's avatar
Mark Friedrichs committed
223

Mark Friedrichs's avatar
Mark Friedrichs committed
224
    RealOpenMMVector& switchDeriviatives  = getSwitchDeriviative();
Mark Friedrichs's avatar
Mark Friedrichs committed
225

Mark Friedrichs's avatar
Mark Friedrichs committed
226
    // ---------------------------------------------------------------------------------------
227

Mark Friedrichs's avatar
Mark Friedrichs committed
228
    // calculate Born radii
Mark Friedrichs's avatar
Mark Friedrichs committed
229

230
    for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
231
232
233
234
235
      
        RealOpenMM radiusI         = atomicRadii[atomI];
        RealOpenMM sum             = zero;
 
        // sum over volumes
Mark Friedrichs's avatar
Mark Friedrichs committed
236

237
        for (int atomJ = 0; atomJ < numberOfAtoms; atomJ++) {
238

239
            if (atomJ != atomI) {
Mark Friedrichs's avatar
Mark Friedrichs committed
240
241
242
  
                RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
                if (_gbviParameters->getPeriodic())
243
                    ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR);
Mark Friedrichs's avatar
Mark Friedrichs committed
244
                else
245
                    ReferenceForce::getDeltaR(atomCoordinates[atomI], atomCoordinates[atomJ], deltaR);
Mark Friedrichs's avatar
Mark Friedrichs committed
246
247
248
249
250
251
   
                RealOpenMM r               = deltaR[ReferenceForce::RIndex];
   
                if (_gbviParameters->getUseCutoff() && r > _gbviParameters->getCutoffDistance())
                    continue;
   
252
                sum  += CpuGBVI::getVolume(r, radiusI, scaledRadii[atomJ]);
Mark Friedrichs's avatar
Mark Friedrichs committed
253
254
255
256
   
            }
        }

257
258
        RealOpenMM atomicRadius3 = POW(radiusI, minusThree);
        if (_gbviParameters->getBornRadiusScalingMethod() != GBVIParameters::QuinticSpline) {
Mark Friedrichs's avatar
Mark Friedrichs committed
259
           sum                        = atomicRadius3 - sum; 
260
           bornRadii[atomI]           = POW(sum, minusOneThird);
Mark Friedrichs's avatar
Mark Friedrichs committed
261
262
263
           switchDeriviatives[atomI]  = one; 
        } else {
           RealOpenMM bornRadius, switchDeriviative;
264
265
           computeBornRadiiUsingQuinticSpline(atomicRadius3, sum, gbviParameters, 
                                              &bornRadius, &switchDeriviative);
Mark Friedrichs's avatar
Mark Friedrichs committed
266
267
268
269
           bornRadii[atomI]           = bornRadius;
           switchDeriviatives[atomI]  = switchDeriviative;
        }    
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
270
271
272
273
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
274
    Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
Mark Friedrichs's avatar
Mark Friedrichs committed
275

Mark Friedrichs's avatar
Mark Friedrichs committed
276
277
278
    @param r                   distance between atoms i & j
    @param R                   atomic radius
    @param S                   scaled atomic radius
Mark Friedrichs's avatar
Mark Friedrichs committed
279

Mark Friedrichs's avatar
Mark Friedrichs committed
280
    @return volume
Mark Friedrichs's avatar
Mark Friedrichs committed
281

Mark Friedrichs's avatar
Mark Friedrichs committed
282
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
283

284
RealOpenMM CpuGBVI::getVolume(RealOpenMM r, RealOpenMM R, RealOpenMM S) {
Mark Friedrichs's avatar
Mark Friedrichs committed
285

Mark Friedrichs's avatar
Mark Friedrichs committed
286
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
287

288
289
    static const RealOpenMM zero         = static_cast<RealOpenMM>( 0.0);
    static const RealOpenMM minusThree   = static_cast<RealOpenMM>(-3.0);
Mark Friedrichs's avatar
Mark Friedrichs committed
290

Mark Friedrichs's avatar
Mark Friedrichs committed
291
    RealOpenMM              diff         = (S - R);
292
    if (FABS(diff) < r) {
Mark Friedrichs's avatar
Mark Friedrichs committed
293

Mark Friedrichs's avatar
Mark Friedrichs committed
294
        RealOpenMM lowerBound = (R > (r - S)) ? R : (r - S);
295
296
        return (CpuGBVI::getL(r, (r + S),    S) -
                CpuGBVI::getL(r, lowerBound, S));
Mark Friedrichs's avatar
Mark Friedrichs committed
297
 
298
    } else if (r <= diff) {
Mark Friedrichs's avatar
Mark Friedrichs committed
299

300
301
302
        return CpuGBVI::getL(r, (r + S), S) -
               CpuGBVI::getL(r, (r - S), S) + 
               POW(R, minusThree);
Mark Friedrichs's avatar
Mark Friedrichs committed
303

Mark Friedrichs's avatar
Mark Friedrichs committed
304
305
306
    } else {
        return zero;
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
307
308
309
310
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
311
    Get L (used in analytical solution for volume integrals) 
Mark Friedrichs's avatar
Mark Friedrichs committed
312

Mark Friedrichs's avatar
Mark Friedrichs committed
313
314
315
    @param r                   distance between atoms i & j
    @param R                   atomic radius
    @param S                   scaled atomic radius
Mark Friedrichs's avatar
Mark Friedrichs committed
316

Mark Friedrichs's avatar
Mark Friedrichs committed
317
    @return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
Mark Friedrichs's avatar
Mark Friedrichs committed
318

Mark Friedrichs's avatar
Mark Friedrichs committed
319
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
320

321
RealOpenMM CpuGBVI::getL(RealOpenMM r, RealOpenMM x, RealOpenMM S) {
Mark Friedrichs's avatar
Mark Friedrichs committed
322

Mark Friedrichs's avatar
Mark Friedrichs committed
323
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
324

325
326
327
328
329
    static const RealOpenMM one           = static_cast<RealOpenMM>(1.0);
    static const RealOpenMM threeHalves   = static_cast<RealOpenMM>(1.5);
    static const RealOpenMM third         = static_cast<RealOpenMM>((1.0/3.0));
    static const RealOpenMM fourth        = static_cast<RealOpenMM>(0.25);
    static const RealOpenMM eighth        = static_cast<RealOpenMM>(0.125);
Mark Friedrichs's avatar
Mark Friedrichs committed
330

Mark Friedrichs's avatar
Mark Friedrichs committed
331
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
332

Mark Friedrichs's avatar
Mark Friedrichs committed
333
    RealOpenMM rInv   = one/r;
Mark Friedrichs's avatar
Mark Friedrichs committed
334

Mark Friedrichs's avatar
Mark Friedrichs committed
335
336
337
    RealOpenMM xInv   = one/x;
    RealOpenMM xInv2  = xInv*xInv;
    RealOpenMM xInv3  = xInv2*xInv;
Mark Friedrichs's avatar
Mark Friedrichs committed
338

Mark Friedrichs's avatar
Mark Friedrichs committed
339
    RealOpenMM diff2  = (r + S)*(r - S);
Mark Friedrichs's avatar
Mark Friedrichs committed
340

341
    return (threeHalves*xInv)*((xInv*fourth*rInv) - (xInv2*third) + (diff2*xInv3*eighth*rInv));
Mark Friedrichs's avatar
Mark Friedrichs committed
342
343
344
345
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
346
    Get partial derivative of L wrt r
Mark Friedrichs's avatar
Mark Friedrichs committed
347

Mark Friedrichs's avatar
Mark Friedrichs committed
348
349
350
    @param r                   distance between atoms i & j
    @param R                   atomic radius
    @param S                   scaled atomic radius
Mark Friedrichs's avatar
Mark Friedrichs committed
351

Mark Friedrichs's avatar
Mark Friedrichs committed
352
    @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
Mark Friedrichs's avatar
Mark Friedrichs committed
353

Mark Friedrichs's avatar
Mark Friedrichs committed
354
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
355

356
RealOpenMM CpuGBVI::dL_dr(RealOpenMM r, RealOpenMM x, RealOpenMM S) {
Mark Friedrichs's avatar
Mark Friedrichs committed
357

Mark Friedrichs's avatar
Mark Friedrichs committed
358
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
359

360
361
362
363
364
365
    static const RealOpenMM one           = static_cast<RealOpenMM>(1.0);
    static const RealOpenMM threeHalves   = static_cast<RealOpenMM>(1.5);
    static const RealOpenMM threeEights   = static_cast<RealOpenMM>(0.375);
    static const RealOpenMM third         = static_cast<RealOpenMM>((1.0/3.0));
    static const RealOpenMM fourth        = static_cast<RealOpenMM>(0.25);
    static const RealOpenMM eighth        = static_cast<RealOpenMM>(0.125);
Mark Friedrichs's avatar
Mark Friedrichs committed
366

Mark Friedrichs's avatar
Mark Friedrichs committed
367
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
368

Mark Friedrichs's avatar
Mark Friedrichs committed
369
370
    RealOpenMM rInv   = one/r;
    RealOpenMM rInv2  = rInv*rInv;
Mark Friedrichs's avatar
Mark Friedrichs committed
371

Mark Friedrichs's avatar
Mark Friedrichs committed
372
373
374
    RealOpenMM xInv   = one/x;
    RealOpenMM xInv2  = xInv*xInv;
    RealOpenMM xInv3  = xInv2*xInv;
Mark Friedrichs's avatar
Mark Friedrichs committed
375

Mark Friedrichs's avatar
Mark Friedrichs committed
376
    RealOpenMM diff2  = (r + S)*(r - S);
Mark Friedrichs's avatar
Mark Friedrichs committed
377

378
    return ((-threeHalves*xInv2*rInv2)*(fourth + eighth*diff2*xInv2) + threeEights*xInv3*xInv);
Mark Friedrichs's avatar
Mark Friedrichs committed
379
380
381
382
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
383
    Get partial derivative of L wrt x
Mark Friedrichs's avatar
Mark Friedrichs committed
384

Mark Friedrichs's avatar
Mark Friedrichs committed
385
386
387
    @param r                   distance between atoms i & j
    @param R                   atomic radius
    @param S                   scaled atomic radius
Mark Friedrichs's avatar
Mark Friedrichs committed
388

Mark Friedrichs's avatar
Mark Friedrichs committed
389
    @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
Mark Friedrichs's avatar
Mark Friedrichs committed
390

Mark Friedrichs's avatar
Mark Friedrichs committed
391
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
392

393
RealOpenMM CpuGBVI::dL_dx(RealOpenMM r, RealOpenMM x, RealOpenMM S) {
Mark Friedrichs's avatar
Mark Friedrichs committed
394

Mark Friedrichs's avatar
Mark Friedrichs committed
395
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
396

397
398
399
400
    static const RealOpenMM one           = static_cast<RealOpenMM>( 1.0);
    static const RealOpenMM half          = static_cast<RealOpenMM>( 0.5);
    static const RealOpenMM threeHalvesM  = static_cast<RealOpenMM>(-1.5);
    static const RealOpenMM third         = static_cast<RealOpenMM>( (1.0/3.0));
Mark Friedrichs's avatar
Mark Friedrichs committed
401

Mark Friedrichs's avatar
Mark Friedrichs committed
402
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
403

Mark Friedrichs's avatar
Mark Friedrichs committed
404
    RealOpenMM rInv   = one/r;
Mark Friedrichs's avatar
Mark Friedrichs committed
405

Mark Friedrichs's avatar
Mark Friedrichs committed
406
407
408
    RealOpenMM xInv   = one/x;
    RealOpenMM xInv2  = xInv*xInv;
    RealOpenMM xInv3  = xInv2*xInv;
Mark Friedrichs's avatar
Mark Friedrichs committed
409

Mark Friedrichs's avatar
Mark Friedrichs committed
410
    RealOpenMM diff   = (r + S)*(r - S);
Mark Friedrichs's avatar
Mark Friedrichs committed
411

412
    return (threeHalvesM*xInv3)*((half*rInv) - xInv + (half*diff*xInv2*rInv));
Mark Friedrichs's avatar
Mark Friedrichs committed
413
414
415
416
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
417
    Sgb function
Mark Friedrichs's avatar
Mark Friedrichs committed
418

Mark Friedrichs's avatar
Mark Friedrichs committed
419
    @param t                   r*r*G_i*G_j
Mark Friedrichs's avatar
Mark Friedrichs committed
420

Mark Friedrichs's avatar
Mark Friedrichs committed
421
    @return Sgb (top of p. 1694 of Labute paper [JCC 29 p. 1693-1698 2008])
Mark Friedrichs's avatar
Mark Friedrichs committed
422

Mark Friedrichs's avatar
Mark Friedrichs committed
423
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
424

425
RealOpenMM CpuGBVI::Sgb(RealOpenMM t) {
Mark Friedrichs's avatar
Mark Friedrichs committed
426

Mark Friedrichs's avatar
Mark Friedrichs committed
427
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
428

Mark Friedrichs's avatar
Mark Friedrichs committed
429
    // static const char* methodName = "CpuGBVI::Sgb";
Mark Friedrichs's avatar
Mark Friedrichs committed
430

431
432
433
    static const RealOpenMM zero    = static_cast<RealOpenMM>(0.0);
    static const RealOpenMM one     = static_cast<RealOpenMM>(1.0);
    static const RealOpenMM fourth  = static_cast<RealOpenMM>(0.25);
Mark Friedrichs's avatar
Mark Friedrichs committed
434

Mark Friedrichs's avatar
Mark Friedrichs committed
435
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
436

437
    return ((t != zero) ? one/SQRT((one + (fourth*EXP(-t))/t)) : zero);
Mark Friedrichs's avatar
Mark Friedrichs committed
438
439
440
441
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
442
    Get GB/VI energy
Mark Friedrichs's avatar
Mark Friedrichs committed
443

Mark Friedrichs's avatar
Mark Friedrichs committed
444
445
    @param atomCoordinates     atomic coordinates
    @param partialCharges      partial charges
Mark Friedrichs's avatar
Mark Friedrichs committed
446

Mark Friedrichs's avatar
Mark Friedrichs committed
447
    @return energy
Mark Friedrichs's avatar
Mark Friedrichs committed
448

Mark Friedrichs's avatar
Mark Friedrichs committed
449
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
450

451
RealOpenMM CpuGBVI::computeBornEnergy(const vector<RealVec>& atomCoordinates, const RealOpenMMVector& partialCharges) {
Mark Friedrichs's avatar
Mark Friedrichs committed
452

Mark Friedrichs's avatar
Mark Friedrichs committed
453
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
454

455
456
457
458
459
460
461
462
    static const RealOpenMM zero          = static_cast<RealOpenMM>(0.0);
    static const RealOpenMM one           = static_cast<RealOpenMM>(1.0);
    static const RealOpenMM two           = static_cast<RealOpenMM>(2.0);
    static const RealOpenMM three         = static_cast<RealOpenMM>(3.0);
    static const RealOpenMM four          = static_cast<RealOpenMM>(4.0);
    static const RealOpenMM half          = static_cast<RealOpenMM>(0.5);
    static const RealOpenMM fourth        = static_cast<RealOpenMM>(0.25);
    static const RealOpenMM eighth        = static_cast<RealOpenMM>(0.125);
Mark Friedrichs's avatar
Mark Friedrichs committed
463

Mark Friedrichs's avatar
Mark Friedrichs committed
464
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
465

Mark Friedrichs's avatar
Mark Friedrichs committed
466
467
468
469
470
    const GBVIParameters* gbviParameters       = getGBVIParameters();
    const RealOpenMM preFactor                 = gbviParameters->getElectricConstant();
    const int numberOfAtoms                    = gbviParameters->getNumberOfAtoms();
    const RealOpenMMVector& atomicRadii        = gbviParameters->getAtomicRadii();
    const RealOpenMMVector& gammaParameters    = gbviParameters->getGammaParameters();
Mark Friedrichs's avatar
Mark Friedrichs committed
471

Mark Friedrichs's avatar
Mark Friedrichs committed
472
    // compute Born radii
Mark Friedrichs's avatar
Mark Friedrichs committed
473

474
475
    RealOpenMMVector bornRadii(numberOfAtoms);
    computeBornRadii(atomCoordinates, bornRadii);
Mark Friedrichs's avatar
Mark Friedrichs committed
476

Mark Friedrichs's avatar
Mark Friedrichs committed
477
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
478

Mark Friedrichs's avatar
Mark Friedrichs committed
479
480
481
    // Eq.2 of Labute paper [JCC 29 p. 1693-1698 2008]
    // to minimze roundoff error sum cavityEnergy separately since in general much
    // smaller than other contributions
Mark Friedrichs's avatar
Mark Friedrichs committed
482

Mark Friedrichs's avatar
Mark Friedrichs committed
483
484
    RealOpenMM energy                 = zero;
    RealOpenMM cavityEnergy           = zero;
Mark Friedrichs's avatar
Mark Friedrichs committed
485

486
    for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
487
 
Mark Friedrichs's avatar
Mark Friedrichs committed
488
        RealOpenMM partialChargeI   = partialCharges[atomI];
Mark Friedrichs's avatar
Mark Friedrichs committed
489
 
Mark Friedrichs's avatar
Mark Friedrichs committed
490
491
492
493
494
495
496
497
498
        // self-energy term
 
        RealOpenMM  atomIEnergy     = half*partialChargeI/bornRadii[atomI];
 
        // cavity term
 
        RealOpenMM ratio            = (atomicRadii[atomI]/bornRadii[atomI]);
        cavityEnergy               += gammaParameters[atomI]*ratio*ratio*ratio;
 
499
        for (int atomJ = atomI + 1; atomJ < numberOfAtoms; atomJ++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
500
501
502
 
            RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
            if (_gbviParameters->getPeriodic())
503
                ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR);
Mark Friedrichs's avatar
Mark Friedrichs committed
504
            else
505
                ReferenceForce::getDeltaR(atomCoordinates[atomI], atomCoordinates[atomJ], deltaR);
Mark Friedrichs's avatar
Mark Friedrichs committed
506
507
508
509
510
            if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
                continue;
  
            RealOpenMM r2           = deltaR[ReferenceForce::R2Index];
            RealOpenMM t            = fourth*r2/(bornRadii[atomI]*bornRadii[atomJ]);         
511
            atomIEnergy            += partialCharges[atomJ]*Sgb(t)/deltaR[ReferenceForce::RIndex];
Mark Friedrichs's avatar
Mark Friedrichs committed
512
513
514
515
516
517
        }
 
        energy += two*partialChargeI*atomIEnergy;
    }
    energy               *= preFactor;
    energy               -= cavityEnergy;
Mark Friedrichs's avatar
Mark Friedrichs committed
518

Mark Friedrichs's avatar
Mark Friedrichs committed
519
520
521
522
    RealOpenMM conversion = static_cast<RealOpenMM>(gbviParameters->getTau());  
    return (conversion*energy);
 
}
Mark Friedrichs's avatar
Mark Friedrichs committed
523

Mark Friedrichs's avatar
Mark Friedrichs committed
524
525
/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
526
    Get GB/VI forces
Mark Friedrichs's avatar
Mark Friedrichs committed
527

Mark Friedrichs's avatar
Mark Friedrichs committed
528
529
530
    @param atomCoordinates     atomic coordinates
    @param partialCharges      partial charges
    @param forces              forces
Mark Friedrichs's avatar
Mark Friedrichs committed
531

Mark Friedrichs's avatar
Mark Friedrichs committed
532
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
533
534


535
536
void CpuGBVI::computeBornForces(std::vector<RealVec>& atomCoordinates, const RealOpenMMVector& partialCharges,
                                 std::vector<OpenMM::RealVec>& inputForces) {
Mark Friedrichs's avatar
Mark Friedrichs committed
537

Mark Friedrichs's avatar
Mark Friedrichs committed
538
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
539

540
541
542
543
544
545
546
547
548
    static const RealOpenMM zero               = static_cast<RealOpenMM>(0.0);
    static const RealOpenMM one                = static_cast<RealOpenMM>(1.0);
    static const RealOpenMM two                = static_cast<RealOpenMM>(2.0);
    static const RealOpenMM three              = static_cast<RealOpenMM>(3.0);
    static const RealOpenMM four               = static_cast<RealOpenMM>(4.0);
    static const RealOpenMM half               = static_cast<RealOpenMM>(0.5);
    static const RealOpenMM oneThird           = static_cast<RealOpenMM>((1.0/3.0));
    static const RealOpenMM fourth             = static_cast<RealOpenMM>(0.25);
    static const RealOpenMM eighth             = static_cast<RealOpenMM>(0.125);
Mark Friedrichs's avatar
Mark Friedrichs committed
549

Mark Friedrichs's avatar
Mark Friedrichs committed
550
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
551

Mark Friedrichs's avatar
Mark Friedrichs committed
552
553
554
555
    const GBVIParameters* gbviParameters       = getGBVIParameters();
    const int numberOfAtoms                    = gbviParameters->getNumberOfAtoms();
    const RealOpenMMVector& atomicRadii        = gbviParameters->getAtomicRadii();
    const RealOpenMMVector& gammaParameters    = gbviParameters->getGammaParameters();
Mark Friedrichs's avatar
Mark Friedrichs committed
556

Mark Friedrichs's avatar
Mark Friedrichs committed
557
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
558

Mark Friedrichs's avatar
Mark Friedrichs committed
559
    // constants
Mark Friedrichs's avatar
Mark Friedrichs committed
560

Mark Friedrichs's avatar
Mark Friedrichs committed
561
    const RealOpenMM preFactor                 = two*gbviParameters->getElectricConstant();
Mark Friedrichs's avatar
Mark Friedrichs committed
562

Mark Friedrichs's avatar
Mark Friedrichs committed
563
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
564

Mark Friedrichs's avatar
Mark Friedrichs committed
565
    // compute Born radii
Mark Friedrichs's avatar
Mark Friedrichs committed
566

567
568
    RealOpenMMVector bornRadii(numberOfAtoms);
    computeBornRadii(atomCoordinates, bornRadii);
Mark Friedrichs's avatar
Mark Friedrichs committed
569

Mark Friedrichs's avatar
Mark Friedrichs committed
570
    // set energy/forces to zero
Mark Friedrichs's avatar
Mark Friedrichs committed
571

572
573
    std::vector<OpenMM::RealVec> forces(numberOfAtoms);
    for (int ii = 0; ii < numberOfAtoms; ii++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
574
575
576
577
        forces[ii][0] = zero;
        forces[ii][1] = zero;
        forces[ii][2] = zero;
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
578

579
    RealOpenMMVector bornForces(numberOfAtoms, 0.0);
Mark Friedrichs's avatar
Mark Friedrichs committed
580

Mark Friedrichs's avatar
Mark Friedrichs committed
581
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
582

Mark Friedrichs's avatar
Mark Friedrichs committed
583
    // first main loop
Mark Friedrichs's avatar
Mark Friedrichs committed
584

585
    for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
586
 
Mark Friedrichs's avatar
Mark Friedrichs committed
587
588
589
590
        // partial of polar term wrt Born radius
        // and (dGpol/dr)(dr/dx)
 
        RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
591
        for (int atomJ = atomI; atomJ < numberOfAtoms; atomJ++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
592
593
594
 
            RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
            if (_gbviParameters->getPeriodic())
595
                ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR);
Mark Friedrichs's avatar
Mark Friedrichs committed
596
            else
597
                ReferenceForce::getDeltaR(atomCoordinates[atomI], atomCoordinates[atomJ], deltaR);
Mark Friedrichs's avatar
Mark Friedrichs committed
598
599
            if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
                continue;
Mark Friedrichs's avatar
Mark Friedrichs committed
600
  
Mark Friedrichs's avatar
Mark Friedrichs committed
601
            RealOpenMM r2                 = deltaR[ReferenceForce::R2Index];
Mark Friedrichs's avatar
Mark Friedrichs committed
602
603
604
605
606
607
608
609
  
            RealOpenMM deltaX             = deltaR[ReferenceForce::XIndex];
            RealOpenMM deltaY             = deltaR[ReferenceForce::YIndex];
            RealOpenMM deltaZ             = deltaR[ReferenceForce::ZIndex];
  
            RealOpenMM alpha2_ij          = bornRadii[atomI]*bornRadii[atomJ];
            RealOpenMM D_ij               = r2/(four*alpha2_ij);
  
610
            RealOpenMM expTerm            = EXP(-D_ij);
Mark Friedrichs's avatar
Mark Friedrichs committed
611
            RealOpenMM denominator2       = r2 + alpha2_ij*expTerm; 
612
            RealOpenMM denominator        = SQRT(denominator2); 
Mark Friedrichs's avatar
Mark Friedrichs committed
613
614
            
            RealOpenMM Gpol               = (partialChargeI*partialCharges[atomJ])/denominator; 
615
            RealOpenMM dGpol_dr           = -Gpol*(one - fourth*expTerm)/denominator2;  
Mark Friedrichs's avatar
Mark Friedrichs committed
616
  
617
            RealOpenMM dGpol_dalpha2_ij   = -half*Gpol*expTerm*(one + D_ij)/denominator2;
Mark Friedrichs's avatar
Mark Friedrichs committed
618
  
619
            if (atomI != atomJ) {
Mark Friedrichs's avatar
Mark Friedrichs committed
620
621
622
623
624
625
  
                bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
   
                deltaX            *= dGpol_dr;
                deltaY            *= dGpol_dr;
                deltaZ            *= dGpol_dr;
Mark Friedrichs's avatar
Mark Friedrichs committed
626
   
Mark Friedrichs's avatar
Mark Friedrichs committed
627
628
629
                forces[atomI][0]  += deltaX;
                forces[atomI][1]  += deltaY;
                forces[atomI][2]  += deltaZ;
Mark Friedrichs's avatar
Mark Friedrichs committed
630
  
Mark Friedrichs's avatar
Mark Friedrichs committed
631
632
633
                forces[atomJ][0]  -= deltaX;
                forces[atomJ][1]  -= deltaY;
                forces[atomJ][2]  -= deltaZ;
Mark Friedrichs's avatar
Mark Friedrichs committed
634
   
Mark Friedrichs's avatar
Mark Friedrichs committed
635
636
637
638
639
            }
            bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
 
        }
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
640

Mark Friedrichs's avatar
Mark Friedrichs committed
641
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
642

Mark Friedrichs's avatar
Mark Friedrichs committed
643
    // second main loop: (dGpol/dBornRadius)(dBornRadius/dr)(dr/dx)
Mark Friedrichs's avatar
Mark Friedrichs committed
644

Mark Friedrichs's avatar
Mark Friedrichs committed
645
646
    // dGpol/dBornRadius) = bornForces[]
    // dBornRadius/dr     = (1/3)*(bR**4)*(dV/dr)
Mark Friedrichs's avatar
Mark Friedrichs committed
647

Mark Friedrichs's avatar
Mark Friedrichs committed
648
649
    const RealOpenMMVector& scaledRadii           = gbviParameters->getScaledRadii();
    const RealOpenMMVector& switchDeriviative     = getSwitchDeriviative();
650
    for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
651
652
653
654
655
656
657
658
659
660
661
 
        RealOpenMM R        = atomicRadii[atomI];
 
        // partial of cavity term wrt Born radius
       
        RealOpenMM  ratio   = (atomicRadii[atomI]/bornRadii[atomI]);
        bornForces[atomI]  += (three*gammaParameters[atomI]*ratio*ratio*ratio)/bornRadii[atomI]; 
 
        RealOpenMM b2       = bornRadii[atomI]*bornRadii[atomI];
        bornForces[atomI]  *= switchDeriviative[atomI]*oneThird*b2*b2;
 
662
        for (int atomJ = 0; atomJ < numberOfAtoms; atomJ++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
663
 
664
            if (atomJ != atomI) {
Mark Friedrichs's avatar
Mark Friedrichs committed
665
666
667
668
669
670
671
  
                RealOpenMM deltaX             = atomCoordinates[atomJ][0] - atomCoordinates[atomI][0];
                RealOpenMM deltaY             = atomCoordinates[atomJ][1] - atomCoordinates[atomI][1];
                RealOpenMM deltaZ             = atomCoordinates[atomJ][2] - atomCoordinates[atomI][2];
        
                RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
                if (_gbviParameters->getPeriodic())
672
                    ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR);
Mark Friedrichs's avatar
Mark Friedrichs committed
673
                else
674
                    ReferenceForce::getDeltaR(atomCoordinates[atomI], atomCoordinates[atomJ], deltaR);
Mark Friedrichs's avatar
Mark Friedrichs committed
675
676
677
678
679
680
681
                if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
                    continue;
       
                RealOpenMM r2                 = deltaR[ReferenceForce::R2Index];
                           deltaX             = deltaR[ReferenceForce::XIndex];
                           deltaY             = deltaR[ReferenceForce::YIndex];
                           deltaZ             = deltaR[ReferenceForce::ZIndex];
Mark Friedrichs's avatar
Mark Friedrichs committed
682
   
683
                RealOpenMM r                  = SQRT(r2);
Mark Friedrichs's avatar
Mark Friedrichs committed
684
685
686
687
688
689
690
691
    
                RealOpenMM S                  = scaledRadii[atomJ];
                RealOpenMM diff               = (S - R);
   
                RealOpenMM de                 = zero;
   
                // find dRb/dr, where Rb is the Born radius
   
692
693
694
695
                if (FABS(diff) < r) {
                    de = CpuGBVI::dL_dr(r, r+S, S) + CpuGBVI::dL_dx(r, r+S, S);   
                    if (R > (r - S)) {
                       de -= CpuGBVI::dL_dr(r, R, S);  
Mark Friedrichs's avatar
Mark Friedrichs committed
696
                    } else {
697
                       de -= (CpuGBVI::dL_dr(r, (r-S), S) + CpuGBVI::dL_dx(r, (r-S), S));
Mark Friedrichs's avatar
Mark Friedrichs committed
698
                    }
699
700
701
                } else if (r < (S - R)) {
                    de  = CpuGBVI::dL_dr(r, r+S, S) + CpuGBVI::dL_dx(r, r+S, S);   
                    de -= (CpuGBVI::dL_dr(r, r-S, S) + CpuGBVI::dL_dx(r, r-S, S));   
Mark Friedrichs's avatar
Mark Friedrichs committed
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
                }
   
                 // de = (dG/dRb)(dRb/dr)
   
                de                      *= bornForces[atomI]/r;
   
                deltaX                  *= de;
                deltaY                  *= de;
                deltaZ                  *= de;
       
                forces[atomI][0]        += deltaX;
                forces[atomI][1]        += deltaY;
                forces[atomI][2]        += deltaZ;
     
                forces[atomJ][0]        -= deltaX;
                forces[atomJ][1]        -= deltaY;
                forces[atomJ][2]        -= deltaZ;
   
            }
        }
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
723

724
    //printGbvi(atomCoordinates, partialCharges, bornRadii, bornForces, forces, "GBVI: Post loop2", stderr);
Mark Friedrichs's avatar
Mark Friedrichs committed
725

Mark Friedrichs's avatar
Mark Friedrichs committed
726
    // convert from cal to Joule & apply prefactor tau = (1/diel_solute - 1/diel_solvent)
Mark Friedrichs's avatar
Mark Friedrichs committed
727

Mark Friedrichs's avatar
Mark Friedrichs committed
728
    RealOpenMM conversion = static_cast<RealOpenMM>(gbviParameters->getTau());  
729
    for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
730
731
732
733
       inputForces[atomI][0] += conversion*forces[atomI][0];
       inputForces[atomI][1] += conversion*forces[atomI][1];
       inputForces[atomI][2] += conversion*forces[atomI][2];
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
734
735
736

}

Mark Friedrichs's avatar
Mark Friedrichs committed
737
738
739
740
741
742
743
744
745
746
747
748
749
750
/**---------------------------------------------------------------------------------------

    Print GB/VI parameters, radii, forces, ...

    @param atomCoordinates     atomic coordinates
    @param partialCharges      partial charges
    @param bornRadii           Born radii (may be empty)
    @param bornForces          Born forces (may be empty)
    @param forces              forces (may be empty)
    @param idString            id string (who is calling)
    @param log                 log file

    --------------------------------------------------------------------------------------- */

751
752
753
754
755
void CpuGBVI::printGbvi(const std::vector<OpenMM::RealVec>& atomCoordinates, const RealOpenMMVector& partialCharges,
                        const RealOpenMMVector& bornRadii,
                        const RealOpenMMVector& bornForces,
                        const std::vector<OpenMM::RealVec>& forces,
                        const std::string& idString, FILE* log) {
Mark Friedrichs's avatar
Mark Friedrichs committed
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775

    // ---------------------------------------------------------------------------------------

    const GBVIParameters* gbviParameters       = getGBVIParameters();
    const int numberOfAtoms                    = gbviParameters->getNumberOfAtoms();
    const RealOpenMMVector& atomicRadii        = gbviParameters->getAtomicRadii();
    const RealOpenMMVector& gammaParameters    = gbviParameters->getGammaParameters();

    // ---------------------------------------------------------------------------------------

    // constants

    const RealOpenMM preFactor                 = 2.0*gbviParameters->getElectricConstant();

    // ---------------------------------------------------------------------------------------

    const RealOpenMMVector& scaledRadii       = gbviParameters->getScaledRadii();
    const RealOpenMMVector& switchDeriviative = getSwitchDeriviative();
    RealOpenMM tau                            = static_cast<RealOpenMM>(gbviParameters->getTau());  

776
777
    int useComparisonFormat                   = 1;

778
779
780
781
782
    (void) fprintf(log, "Reference Gbvi     %s atoms=%d\n", idString.c_str(), numberOfAtoms);
    (void) fprintf(log, "    tau            %15.7e\n", tau); 
    (void) fprintf(log, "    scaleMethod    %d (QuinticEnum=%d)\n", 
                    _gbviParameters->getBornRadiusScalingMethod(), GBVIParameters::QuinticSpline);
    (void) fprintf(log, "    preFactor      %15.7e)\n", preFactor);
Mark Friedrichs's avatar
Mark Friedrichs committed
783
 
784
785
786
787
788
789
    if (useComparisonFormat) {
        (void) fprintf(log, "  br bF swd r scR tau*gamma q)\n");
        for (unsigned int atomI = 0; atomI < static_cast<unsigned int>(numberOfAtoms); atomI++) {
            (void) fprintf(log, "%6d ", atomI);
            if (bornRadii.size() > atomI) {
                 (void) fprintf(log, "%15.7e ", bornRadii[atomI]);
790
            }
791
792
            if (bornForces.size() > atomI) {
                 (void) fprintf(log, "%15.7e ", tau*bornForces[atomI]);    
793
            }
794
            (void) fprintf(log, " %15.7e %15.7e %15.7e %15.7e %15.7e",
795
796
797
798
                            switchDeriviative[atomI],    
                            atomicRadii[atomI],    
                            scaledRadii[atomI],    
                            tau*gammaParameters[atomI],    
799
800
                            partialCharges[atomI]);
            (void) fprintf(log, "\n");
801
802
        }   
    } else {
803
804
        for (unsigned int atomI = 0; atomI < static_cast<unsigned int>(numberOfAtoms); atomI++) {
            (void) fprintf(log, "%6d r=%15.7e rSc=%15.7e swd=%15.7e tau*gam=%15.7e q=%15.7e", atomI,
805
806
807
808
                            atomicRadii[atomI],    
                            scaledRadii[atomI],    
                            switchDeriviative[atomI],    
                            tau*gammaParameters[atomI],    
809
810
811
                            partialCharges[atomI]);
            if (bornRadii.size() > atomI) {
                 (void) fprintf(log, " bR=%15.7e", bornRadii[atomI]);
812
            }
813
814
            if (bornForces.size() > atomI) {
                 (void) fprintf(log, " tau*bF=%15.7e", tau*bornForces[atomI]);    
815
            }
816
            (void) fprintf(log, "\n");
817
        }   
Mark Friedrichs's avatar
Mark Friedrichs committed
818
819
820
821
822
823
    }   

    return;

}

Mark Friedrichs's avatar
Mark Friedrichs committed
824
825
/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
826
    Use double precision 
Mark Friedrichs's avatar
Mark Friedrichs committed
827

Mark Friedrichs's avatar
Mark Friedrichs committed
828
    Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
Mark Friedrichs's avatar
Mark Friedrichs committed
829

Mark Friedrichs's avatar
Mark Friedrichs committed
830
831
832
    @param r                   distance between atoms i & j
    @param R                   atomic radius
    @param S                   scaled atomic radius
Mark Friedrichs's avatar
Mark Friedrichs committed
833

Mark Friedrichs's avatar
Mark Friedrichs committed
834
    @return volume
Mark Friedrichs's avatar
Mark Friedrichs committed
835

Mark Friedrichs's avatar
Mark Friedrichs committed
836
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
837

838
double CpuGBVI::getVolumeD(double r, double R, double S) {
Mark Friedrichs's avatar
Mark Friedrichs committed
839

Mark Friedrichs's avatar
Mark Friedrichs committed
840
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
841

Mark Friedrichs's avatar
Mark Friedrichs committed
842
843
    static const double zero         = 0.0;
    static const double minusThree   = -3.0;
Mark Friedrichs's avatar
Mark Friedrichs committed
844

Mark Friedrichs's avatar
Mark Friedrichs committed
845
    double              diff    = (S - R);
846
    if (fabs(diff) < r) {
Mark Friedrichs's avatar
Mark Friedrichs committed
847

Mark Friedrichs's avatar
Mark Friedrichs committed
848
       double lowerBound = (R > (r - S)) ? R : (r - S);
Mark Friedrichs's avatar
Mark Friedrichs committed
849

850
851
       return (CpuGBVI::getLD(r, (r + S),    S) -
               CpuGBVI::getLD(r, lowerBound, S));
Mark Friedrichs's avatar
Mark Friedrichs committed
852

853
    } else if (r < diff) {
Mark Friedrichs's avatar
Mark Friedrichs committed
854

855
856
857
       return CpuGBVI::getLD(r, (r + S), S) -
              CpuGBVI::getLD(r, (r - S), S) + 
              pow(R, minusThree);
Mark Friedrichs's avatar
Mark Friedrichs committed
858

Mark Friedrichs's avatar
Mark Friedrichs committed
859
860
861
    } else {
       return zero;
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
862
863
864
865
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
866
    Use double precision 
Mark Friedrichs's avatar
Mark Friedrichs committed
867

Mark Friedrichs's avatar
Mark Friedrichs committed
868
    Get L (used in analytical solution for volume integrals) 
Mark Friedrichs's avatar
Mark Friedrichs committed
869

Mark Friedrichs's avatar
Mark Friedrichs committed
870
871
872
    @param r                   distance between atoms i & j
    @param R                   atomic radius
    @param S                   scaled atomic radius
Mark Friedrichs's avatar
Mark Friedrichs committed
873

Mark Friedrichs's avatar
Mark Friedrichs committed
874
    @return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
Mark Friedrichs's avatar
Mark Friedrichs committed
875

Mark Friedrichs's avatar
Mark Friedrichs committed
876
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
877

878
double CpuGBVI::getLD(double r, double x, double S) {
Mark Friedrichs's avatar
Mark Friedrichs committed
879

Mark Friedrichs's avatar
Mark Friedrichs committed
880
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
881

Mark Friedrichs's avatar
Mark Friedrichs committed
882
883
884
885
886
    static const double one           =  1.0;
    static const double threeHalves   =  1.5;
    static const double third         =  1.0/3.0;
    static const double fourth        =  0.25;
    static const double eighth        =  0.125;
Mark Friedrichs's avatar
Mark Friedrichs committed
887

Mark Friedrichs's avatar
Mark Friedrichs committed
888
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
889

Mark Friedrichs's avatar
Mark Friedrichs committed
890
    double rInv   = one/r;
Mark Friedrichs's avatar
Mark Friedrichs committed
891

Mark Friedrichs's avatar
Mark Friedrichs committed
892
893
894
    double xInv   = one/x;
    double xInv2  = xInv*xInv;
    double xInv3  = xInv2*xInv;
Mark Friedrichs's avatar
Mark Friedrichs committed
895

Mark Friedrichs's avatar
Mark Friedrichs committed
896
    double diff2  = (r + S)*(r - S);
Mark Friedrichs's avatar
Mark Friedrichs committed
897

898
    return (threeHalves*xInv)*((xInv*fourth*rInv) - (xInv2*third) + (diff2*xInv3*eighth*rInv));
Mark Friedrichs's avatar
Mark Friedrichs committed
899
900
901
902
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
903
    Use double precision 
Mark Friedrichs's avatar
Mark Friedrichs committed
904

Mark Friedrichs's avatar
Mark Friedrichs committed
905
    Get partial derivative of L wrt r
Mark Friedrichs's avatar
Mark Friedrichs committed
906

Mark Friedrichs's avatar
Mark Friedrichs committed
907
908
909
    @param r                   distance between atoms i & j
    @param R                   atomic radius
    @param S                   scaled atomic radius
Mark Friedrichs's avatar
Mark Friedrichs committed
910

Mark Friedrichs's avatar
Mark Friedrichs committed
911
    @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
Mark Friedrichs's avatar
Mark Friedrichs committed
912

Mark Friedrichs's avatar
Mark Friedrichs committed
913
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
914

915
double CpuGBVI::dL_drD(double r, double x, double S) {
Mark Friedrichs's avatar
Mark Friedrichs committed
916

Mark Friedrichs's avatar
Mark Friedrichs committed
917
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
918

Mark Friedrichs's avatar
Mark Friedrichs committed
919
920
921
922
923
924
    static const double one           =  1.0;
    static const double threeHalves   =  1.5;
    static const double threeEights   =  0.375;
    static const double third         =  1.0/3.0;
    static const double fourth        =  0.25;
    static const double eighth        =  0.125;
Mark Friedrichs's avatar
Mark Friedrichs committed
925

Mark Friedrichs's avatar
Mark Friedrichs committed
926
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
927

Mark Friedrichs's avatar
Mark Friedrichs committed
928
929
    double rInv   = one/r;
    double rInv2  = rInv*rInv;
Mark Friedrichs's avatar
Mark Friedrichs committed
930

Mark Friedrichs's avatar
Mark Friedrichs committed
931
932
933
    double xInv   = one/x;
    double xInv2  = xInv*xInv;
    double xInv3  = xInv2*xInv;
Mark Friedrichs's avatar
Mark Friedrichs committed
934

Mark Friedrichs's avatar
Mark Friedrichs committed
935
    double diff2  = (r + S)*(r - S);
Mark Friedrichs's avatar
Mark Friedrichs committed
936

937
    return ((-threeHalves*xInv2*rInv2)*(fourth + eighth*diff2*xInv2) + threeEights*xInv3*xInv);
Mark Friedrichs's avatar
Mark Friedrichs committed
938
939
940
941
}

/**---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
942
    Use double precision 
Mark Friedrichs's avatar
Mark Friedrichs committed
943

Mark Friedrichs's avatar
Mark Friedrichs committed
944
    Get partial derivative of L wrt x
Mark Friedrichs's avatar
Mark Friedrichs committed
945

Mark Friedrichs's avatar
Mark Friedrichs committed
946
947
948
    @param r                   distance between atoms i & j
    @param R                   atomic radius
    @param S                   scaled atomic radius
Mark Friedrichs's avatar
Mark Friedrichs committed
949

Mark Friedrichs's avatar
Mark Friedrichs committed
950
    @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
Mark Friedrichs's avatar
Mark Friedrichs committed
951

Mark Friedrichs's avatar
Mark Friedrichs committed
952
    --------------------------------------------------------------------------------------- */
Mark Friedrichs's avatar
Mark Friedrichs committed
953

954
double CpuGBVI::dL_dxD(double r, double x, double S) {
Mark Friedrichs's avatar
Mark Friedrichs committed
955

Mark Friedrichs's avatar
Mark Friedrichs committed
956
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
957

Mark Friedrichs's avatar
Mark Friedrichs committed
958
959
960
961
    static const double one           =   1.0;
    static const double half          =   0.5;
    static const double threeHalvesM  =  -1.5;
    static const double third         =   1.0/3.0;
Mark Friedrichs's avatar
Mark Friedrichs committed
962

Mark Friedrichs's avatar
Mark Friedrichs committed
963
    // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
964

Mark Friedrichs's avatar
Mark Friedrichs committed
965
    double rInv   = one/r;
Mark Friedrichs's avatar
Mark Friedrichs committed
966

Mark Friedrichs's avatar
Mark Friedrichs committed
967
968
969
    double xInv   = one/x;
    double xInv2  = xInv*xInv;
    double xInv3  = xInv2*xInv;
Mark Friedrichs's avatar
Mark Friedrichs committed
970

Mark Friedrichs's avatar
Mark Friedrichs committed
971
    double diff   = (r + S)*(r - S);
Mark Friedrichs's avatar
Mark Friedrichs committed
972

973
    return (threeHalvesM*xInv3)*((half*rInv) - xInv + (half*diff*xInv2*rInv));
Mark Friedrichs's avatar
Mark Friedrichs committed
974
}