ReferenceLJCoulombIxn.cpp 24.8 KB
Newer Older
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

/* Portions copyright (c) 2006 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>
#include <sstream>
27
#include <complex>
28

29
30
31
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
32
33
#include "ReferenceLJCoulombIxn.h"
#include "ReferenceForce.h"
34
#include "PME.h"
35

36
37
// In case we're using some primitive version of Visual Studio this will
// make sure that erf() and erfc() are defined.
38
#include "openmm/internal/MSVC_erfc.h"
39

40
41
using std::vector;

42
43
44
45
46
47
/**---------------------------------------------------------------------------------------

   ReferenceLJCoulombIxn constructor

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

48
ReferenceLJCoulombIxn::ReferenceLJCoulombIxn( ) : cutoff(false), periodic(false), ewald(false), pme(false) {
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

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

   // static const char* methodName = "\nReferenceLJCoulombIxn::ReferenceLJCoulombIxn";

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

}

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

   ReferenceLJCoulombIxn destructor

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

ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){

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

   // static const char* methodName = "\nReferenceLJCoulombIxn::~ReferenceLJCoulombIxn";

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

}

74
75
76
77
78
79
80
81
82
83
84
85
86
  /**---------------------------------------------------------------------------------------

     Set the force to use a cutoff.

     @param distance            the cutoff distance
     @param neighbors           the neighbor list to use
     @param solventDielectric   the dielectric constant of the bulk solvent

     @return ReferenceForce::DefaultReturn

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

  int ReferenceLJCoulombIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric ) {
87

88
89
90
    cutoff = true;
    cutoffDistance = distance;
    neighborList = &neighbors;
91
92
    krf = pow(cutoffDistance, -3.0)*(solventDielectric-1.0)/(2.0*solventDielectric+1.0);
    crf = (1.0/cutoffDistance)*(3.0*solventDielectric)/(2.0*solventDielectric+1.0);
93

94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
    return ReferenceForce::DefaultReturn;
  }

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

     Set the force to use periodic boundary conditions.  This requires that a cutoff has
     also been set, and the smallest side of the periodic box is at least twice the cutoff
     distance.

     @param boxSize             the X, Y, and Z widths of the periodic box

     @return ReferenceForce::DefaultReturn

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

  int ReferenceLJCoulombIxn::setPeriodic( RealOpenMM* boxSize ) {

    assert(cutoff);
    assert(boxSize[0] >= 2.0*cutoffDistance);
    assert(boxSize[1] >= 2.0*cutoffDistance);
    assert(boxSize[2] >= 2.0*cutoffDistance);
    periodic = true;
    periodicBoxSize[0] = boxSize[0];
    periodicBoxSize[1] = boxSize[1];
    periodicBoxSize[2] = boxSize[2];
    return ReferenceForce::DefaultReturn;

  }

123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
  /**---------------------------------------------------------------------------------------

     Set the force to use Ewald summation.

     @param alpha  the Ewald separation parameter
     @param kmaxx  the largest wave vector in the x direction
     @param kmaxy  the largest wave vector in the y direction
     @param kmaxz  the largest wave vector in the z direction

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

  void ReferenceLJCoulombIxn::setUseEwald(RealOpenMM alpha, int kmaxx, int kmaxy, int kmaxz) {
      alphaEwald = alpha;
      numRx = kmaxx;
      numRy = kmaxy;
      numRz = kmaxz;
      ewald = true;
  }

142
143
144
145
146
  /**---------------------------------------------------------------------------------------

     Set the force to use Particle-Mesh Ewald (PME) summation.

     @param alpha  the Ewald separation parameter
147
     @param gridSize the dimensions of the mesh
148
149
150

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

151
  void ReferenceLJCoulombIxn::setUsePME(RealOpenMM alpha, int meshSize[3]) {
152
      alphaEwald = alpha;
153
154
155
      meshDim[0] = meshSize[0];
      meshDim[1] = meshSize[1];
      meshDim[2] = meshSize[2];
156
157
158
      pme = true;
  }

159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
/**---------------------------------------------------------------------------------------

   Calculate Ewald ixn

   @param numberOfAtoms    number of atoms
   @param atomCoordinates  atom coordinates
   @param atomParameters   atom parameters                             atomParameters[atomIndex][paramterIndex]
   @param exclusions       atom exclusion indices                      exclusions[atomIndex][atomToExcludeIndex]
                           exclusions[atomIndex][0] = number of exclusions
                           exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
                           interacting w/ atom atomIndex
   @param fixedParameters  non atom parameters (not currently used)
   @param forces           force array (forces added)
   @param energyByAtom     atom energy
   @param totalEnergy      total energy

   @return ReferenceForce::DefaultReturn
176

177
   --------------------------------------------------------------------------------------- */
178

179
180
181
int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
                                             RealOpenMM** atomParameters, int** exclusions,
                                             RealOpenMM* fixedParameters, RealOpenMM** forces,
182
                                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const {
183
184
    typedef std::complex<RealOpenMM> d_complex;

185
    static const RealOpenMM epsilon     =  1.0;
186
    static const RealOpenMM one         =  1.0;
187
188
    static const RealOpenMM six         =  6.0;
    static const RealOpenMM twelve      = 12.0;
189

190
    int kmax                            = (ewald ? std::max(numRx, std::max(numRy,numRz)) : 0);
191
    RealOpenMM  factorEwald             = -1 / (4*alphaEwald*alphaEwald);
Peter Eastman's avatar
Peter Eastman committed
192
193
194
    RealOpenMM SQRT_PI                  = sqrt(PI_M);
    RealOpenMM TWO_PI                   = 2.0 * PI_M;
    RealOpenMM recipCoeff               = (RealOpenMM)(ONE_4PI_EPS0*4*PI_M/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
195

196
    RealOpenMM totalSelfEwaldEnergy     = 0.0;
197
198
    RealOpenMM realSpaceEwaldEnergy     = 0.0;
    RealOpenMM recipEnergy              = 0.0;
199
    RealOpenMM totalRecipEnergy         = 0.0;
200
    RealOpenMM vdwEnergy                = 0.0;
201
202
203
204
205
206

// **************************************************************************************
// SELF ENERGY
// **************************************************************************************

    for( int atomID = 0; atomID < numberOfAtoms; atomID++ ){
Peter Eastman's avatar
Peter Eastman committed
207
        RealOpenMM selfEwaldEnergy       = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI);
208
        totalSelfEwaldEnergy            -= selfEwaldEnergy;
209
210

        if( energyByAtom ){
211
           energyByAtom[atomID]         -= selfEwaldEnergy;
212
        }
213
214
    }

215
216
217
218
    if( totalEnergy ){
        *totalEnergy += totalSelfEwaldEnergy;
    }

219
220
221
// **************************************************************************************
// RECIPROCAL SPACE EWALD ENERGY AND FORCES
// **************************************************************************************
222
    // PME
223

224
225
226
227
  if (pme) {
	pme_t          pmedata; /* abstract handle for PME data */
	RealOpenMM virial[3][3];

228
	pme_init(&pmedata,alphaEwald,numberOfAtoms,meshDim,5,1);
229
230
231
232
233
234
235
236
237
238

	pme_exec(pmedata,atomCoordinates,forces,atomParameters,periodicBoxSize,&recipEnergy,virial);

	if( totalEnergy )
       *totalEnergy += recipEnergy;

    if( energyByAtom )
        for(int n = 0; n < numberOfAtoms; n++)
            energyByAtom[n] += recipEnergy;

239
        pme_destroy(pmedata);
240
241
242
243
244
245
246
  }

    // Ewald method

  else if (ewald) {

    // setup reciprocal box
247
248
249
250

           RealOpenMM recipBoxSize[3] = { TWO_PI / periodicBoxSize[0], TWO_PI / periodicBoxSize[1], TWO_PI / periodicBoxSize[2]};


251
    // setup K-vectors
252

253
  #define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
254
255
256
  vector<d_complex> eir(kmax*numberOfAtoms*3);
  vector<d_complex> tab_xy(numberOfAtoms);
  vector<d_complex> tab_qxyz(numberOfAtoms);
257
258
259
260
261
262
263

  if (kmax < 1) {
      std::stringstream message;
      message << " kmax < 1 , Aborting" << std::endl;
      SimTKOpenMMLog::printError( message );
  }

264
265
  for(int i = 0; (i < numberOfAtoms); i++) {
    for(int m = 0; (m < 3); m++)
266
      EIR(0, i, m) = d_complex(1,0);
267

268
    for(int m=0; (m<3); m++)
269
270
      EIR(1, i, m) = d_complex(cos(atomCoordinates[i][m]*recipBoxSize[m]),
                               sin(atomCoordinates[i][m]*recipBoxSize[m]));
271

272
273
    for(int j=2; (j<kmax); j++)
      for(int m=0; (m<3); m++)
274
        EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
275
276
  }

277
278
    // calculate reciprocal space energy and forces

279
280
281
282
283
284
285
286
287
288
289
290
    int lowry = 0;
    int lowrz = 1;

    for(int rx = 0; rx < numRx; rx++) {

      RealOpenMM kx = rx * recipBoxSize[0];

      for(int ry = lowry; ry < numRy; ry++) {

        RealOpenMM ky = ry * recipBoxSize[1];

        if(ry >= 0) {
291
          for(int n = 0; n < numberOfAtoms; n++)
292
            tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
293
294
295
        }

        else {
296
          for(int n = 0; n < numberOfAtoms; n++)
297
            tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
298
299
300
301
302
        }

        for (int rz = lowrz; rz < numRz; rz++) {

          if( rz >= 0) {
303
304
           for( int n = 0; n < numberOfAtoms; n++)
             tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2));
305
306
307
          }

          else {
308
309
            for( int n = 0; n < numberOfAtoms; n++)
              tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2)));
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
          }

          RealOpenMM cs = 0.0f;
          RealOpenMM ss = 0.0f;

          for( int n = 0; n < numberOfAtoms; n++) {
            cs += tab_qxyz[n].real();
            ss += tab_qxyz[n].imag();
          }

          RealOpenMM kz = rz * recipBoxSize[2];
          RealOpenMM k2 = kx * kx + ky * ky + kz * kz;
          RealOpenMM ak = exp(k2*factorEwald) / k2;

          for(int n = 0; n < numberOfAtoms; n++) {
            RealOpenMM force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real());
326
327
328
            forces[n][0] += 2 * recipCoeff * force * kx ;
            forces[n][1] += 2 * recipCoeff * force * ky ;
            forces[n][2] += 2 * recipCoeff * force * kz ;
329
          }
330

331
332
          recipEnergy       = recipCoeff * ak * ( cs * cs + ss * ss);
          totalRecipEnergy += recipEnergy;
333
334
335
336
337
338
339
340

          if( totalEnergy )
             *totalEnergy += recipEnergy;

          if( energyByAtom )
             for(int n = 0; n < numberOfAtoms; n++)
               energyByAtom[n] += recipEnergy;

341
342
343
344
345
          lowrz = 1 - numRz;
        }
        lowry = 1 - numRy;
      }
    }
346
  }
347

348
349
350
351
352
  else {
      std::stringstream message;
      message << " Wrong method for Ewald summation, Aborting" << std::endl;
      SimTKOpenMMLog::printError( message );
  }
353

354

355
356
357
358
// **************************************************************************************
// SHORT-RANGE ENERGY AND FORCES
// **************************************************************************************

359
360
361
    RealOpenMM totalVdwEnergy            = 0.0f;
    RealOpenMM totalRealSpaceEwaldEnergy = 0.0f;

362
363
364
365
    for (int i = 0; i < (int) neighborList->size(); i++) {
       OpenMM::AtomPair pair = (*neighborList)[i];
       int ii = pair.first;
       int jj = pair.second;
366

367
368
       RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
       ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
369
370
371
       RealOpenMM r         = deltaR[0][ReferenceForce::RIndex];
       RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
       RealOpenMM alphaR    = alphaEwald * r;
372
373


Peter Eastman's avatar
Peter Eastman committed
374
375
       RealOpenMM dEdR      = (RealOpenMM) (ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR);
                  dEdR      = (RealOpenMM) (dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
376

377
378
379
380
381
       RealOpenMM sig       = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
       RealOpenMM sig2      = inverseR*sig;
                  sig2     *= sig2;
       RealOpenMM sig6      = sig2*sig2*sig2;
       RealOpenMM eps       = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
382
                  dEdR     += eps*( twelve*sig6 - six )*sig6*inverseR*inverseR;
383

384
       // accumulate forces
385

386
387
388
389
       for( int kk = 0; kk < 3; kk++ ){
          RealOpenMM force  = dEdR*deltaR[0][kk];
          forces[ii][kk]   += force;
          forces[jj][kk]   -= force;
390
391
       }

392
       // accumulate energies
393

394
       realSpaceEwaldEnergy        = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR));
395
       vdwEnergy                   = eps*(sig6-one)*sig6;
396

397
398
       totalVdwEnergy             += vdwEnergy;
       totalRealSpaceEwaldEnergy  += realSpaceEwaldEnergy;
399

400
401
402
403
        if( energyByAtom ){
           energyByAtom[ii] += realSpaceEwaldEnergy + vdwEnergy;
           energyByAtom[jj] += realSpaceEwaldEnergy + vdwEnergy;
        }
404

405
    }
406

407
408
409
    if( totalEnergy )
        *totalEnergy += totalRealSpaceEwaldEnergy + totalVdwEnergy;

Peter Eastman's avatar
Peter Eastman committed
410
411
    // Now subtract off the exclusions, since they were implicitly included in the reciprocal space sum.

412
    RealOpenMM totalExclusionEnergy = 0.0f;
Peter Eastman's avatar
Peter Eastman committed
413
414
415
416
417
418
419
420
421
    for (int i = 0; i < numberOfAtoms; i++)
        for (int j = 1; j <= exclusions[i][0]; j++)
            if (exclusions[i][j] > i) {
               int ii = i;
               int jj = exclusions[i][j];

               RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
               ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
               RealOpenMM r         = deltaR[0][ReferenceForce::RIndex];
422
423
               if (r >= cutoffDistance)
                   continue;
Peter Eastman's avatar
Peter Eastman committed
424
425
               RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
               RealOpenMM alphaR    = alphaEwald * r;
Peter Eastman's avatar
Peter Eastman committed
426
427
               RealOpenMM dEdR      = (RealOpenMM) (ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR);
                          dEdR      = (RealOpenMM) (dEdR * (erf(alphaR) - 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
Peter Eastman's avatar
Peter Eastman committed
428
429
430
431
432
433
434
435
436
437
438

               // accumulate forces

               for( int kk = 0; kk < 3; kk++ ){
                  RealOpenMM force  = dEdR*deltaR[0][kk];
                  forces[ii][kk]   -= force;
                  forces[jj][kk]   += force;
               }

               // accumulate energies

439
               realSpaceEwaldEnergy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erf(alphaR));
Peter Eastman's avatar
Peter Eastman committed
440

441
442
               totalExclusionEnergy += realSpaceEwaldEnergy;
               if( energyByAtom ){
Peter Eastman's avatar
Peter Eastman committed
443
444
                   energyByAtom[ii] -= realSpaceEwaldEnergy;
                   energyByAtom[jj] -= realSpaceEwaldEnergy;
445
               }
Peter Eastman's avatar
Peter Eastman committed
446
447
            }

448
449
    if( totalEnergy )
        *totalEnergy -= totalExclusionEnergy;
450

451
// ***********************************************************************
452
453
454
455
456

   return ReferenceForce::DefaultReturn;
}


457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
/**---------------------------------------------------------------------------------------

   Calculate LJ Coulomb pair ixn

   @param numberOfAtoms    number of atoms
   @param atomCoordinates  atom coordinates
   @param atomParameters   atom parameters                             atomParameters[atomIndex][paramterIndex]
   @param exclusions       atom exclusion indices                      exclusions[atomIndex][atomToExcludeIndex]
                           exclusions[atomIndex][0] = number of exclusions
                           exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
                           interacting w/ atom atomIndex
   @param fixedParameters  non atom parameters (not currently used)
   @param forces           force array (forces added)
   @param energyByAtom     atom energy
   @param totalEnergy      total energy

   @return ReferenceForce::DefaultReturn
474

475
   --------------------------------------------------------------------------------------- */
476

477
478
479
480
481
int ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
                                             RealOpenMM** atomParameters, int** exclusions,
                                             RealOpenMM* fixedParameters, RealOpenMM** forces,
                                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {

482
   if (ewald || pme)
483
        return calculateEwaldIxn(numberOfAtoms, atomCoordinates, atomParameters, exclusions, fixedParameters, forces, energyByAtom, totalEnergy);
484
   if (cutoff) {
485
       for (int i = 0; i < (int) neighborList->size(); i++) {
486
487
488
489
490
491
           OpenMM::AtomPair pair = (*neighborList)[i];
           calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
       }
   }
   else {
       // allocate and initialize exclusion array
492

493
494
495
496
       int* exclusionIndices = new int[numberOfAtoms];
       for( int ii = 0; ii < numberOfAtoms; ii++ ){
          exclusionIndices[ii] = -1;
       }
497

498
       for( int ii = 0; ii < numberOfAtoms; ii++ ){
499

500
          // set exclusions
501

502
503
504
          for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
             exclusionIndices[exclusions[ii][jj]] = ii;
          }
505

506
          // loop over atom pairs
507

508
          for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
509

510
511
512
513
514
             if( exclusionIndices[jj] != ii ){
                 calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
             }
          }
       }
515

516
517
       delete[] exclusionIndices;
   }
518

519
520
   return ReferenceForce::DefaultReturn;
}
521

522
  /**---------------------------------------------------------------------------------------
523

524
     Calculate LJ Coulomb pair ixn between two atoms
525

526
527
528
529
530
531
532
     @param ii               the index of the first atom
     @param jj               the index of the second atom
     @param atomCoordinates  atom coordinates
     @param atomParameters   atom parameters (charges, c6, c12, ...)     atomParameters[atomIndex][paramterIndex]
     @param forces           force array (forces added)
     @param energyByAtom     atom energy
     @param totalEnergy      total energy
533

534
     @return ReferenceForce::DefaultReturn
535

536
     --------------------------------------------------------------------------------------- */
537

Peter Eastman's avatar
Peter Eastman committed
538
int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoordinates,
539
540
                        RealOpenMM** atomParameters, RealOpenMM** forces,
                        RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
541

542
    // ---------------------------------------------------------------------------------------
543

544
    static const std::string methodName = "\nReferenceLJCoulombIxn::calculateOneIxn";
545

546
    // ---------------------------------------------------------------------------------------
547

548
    // constants -- reduce Visual Studio warnings regarding conversions between float & double
549

550
551
552
553
554
555
556
    static const RealOpenMM zero        =  0.0;
    static const RealOpenMM one         =  1.0;
    static const RealOpenMM two         =  2.0;
    static const RealOpenMM three       =  3.0;
    static const RealOpenMM six         =  6.0;
    static const RealOpenMM twelve      = 12.0;
    static const RealOpenMM oneM        = -1.0;
557

558
    static const int threeI             = 3;
559

560
561
562
563
564
565
566
567
568
569
570
    // debug flag

    static const int debug              = -1;

    static const int LastAtomIndex      = 2;

    RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];

    // get deltaR, R2, and R between 2 atoms

    if (periodic)
571
        ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
572
    else
573
        ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] );
574
575
576
577
578
579
580
581
582
583
584

    RealOpenMM r2        = deltaR[0][ReferenceForce::R2Index];
    RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
    RealOpenMM sig       = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
    RealOpenMM sig2      = inverseR*sig;
               sig2     *= sig2;
    RealOpenMM sig6      = sig2*sig2*sig2;

    RealOpenMM eps       = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
    RealOpenMM dEdR      = eps*( twelve*sig6 - six )*sig6;
               if (cutoff)
Peter Eastman's avatar
Peter Eastman committed
585
                   dEdR += (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2));
586
               else
Peter Eastman's avatar
Peter Eastman committed
587
                   dEdR += (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR);
588
               dEdR     *= inverseR*inverseR;
589

590
591
592
593
594
595
596
597
598
599
600
601
602
603
    // accumulate forces

    for( int kk = 0; kk < 3; kk++ ){
       RealOpenMM force  = dEdR*deltaR[0][kk];
       forces[ii][kk]   += force;
       forces[jj][kk]   -= force;
    }

    RealOpenMM energy = 0.0;

    // accumulate energies

    if( totalEnergy || energyByAtom ) {
        if (cutoff)
Peter Eastman's avatar
Peter Eastman committed
604
            energy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf));
605
        else
Peter Eastman's avatar
Peter Eastman committed
606
            energy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR);
607
608
609
610
611
612
613
614
615
        energy += eps*(sig6-one)*sig6;
        if( totalEnergy )
           *totalEnergy += energy;
        if( energyByAtom ){
           energyByAtom[ii] += energy;
           energyByAtom[jj] += energy;
        }
    }

616
    // debug
617
618
619
620
621
622
623

    if( debug == ii ){
       static bool printHeader = false;
       std::stringstream message;
       message << methodName;
       message << std::endl;
       int pairArray[2] = { ii, jj };
624
       if( !printHeader  ){
625
626
627
          printHeader = true;
          message << std::endl;
          message << methodName.c_str() << " a0 k [c q p s] r1 r2  angle dt rp p[] dot cosine angle dEdR*r F[]" << std::endl;
628
       }
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667

       message << std::endl;
       for( int kk = 0; kk < 2; kk++ ){
          message << " Atm " << pairArray[kk] << " [" << atomCoordinates[pairArray[kk]][0] << " " << atomCoordinates[pairArray[kk]][1] << " " << atomCoordinates[pairArray[kk]][2] << "] ";
       }
       message << std::endl << " Delta:";
       for( int kk = 0; kk < (LastAtomIndex - 1); kk++ ){
          message << " [";
          for( int jj = 0; jj < ReferenceForce::LastDeltaRIndex; jj++ ){
             message << deltaR[kk][jj] << " ";
          }
          message << "]";
       }
       message << std::endl;

       for( int kk = 0; kk < 2; kk++ ){
          message << " p" << pairArray[kk] << " [";
          message << atomParameters[pairArray[kk]][0] << " " << atomParameters[pairArray[kk]][1] << " " << atomParameters[pairArray[kk]][2];
          message << "]";
       }
      message << std::endl;

       message << " dEdR=" << dEdR;
       message << " E=" << energy << " force factors: ";
       message << "F=compute force; f=cumulative force";

       message << std::endl << "  ";
       message << " f" << ii << "[";
       SimTKOpenMMUtilities::formatRealStringStream( message, deltaR[0], threeI, dEdR );
       message << "]";

       for( int kk = 0; kk < 2; kk++ ){
          message << " F" <<  pairArray[kk] << " [";
          SimTKOpenMMUtilities::formatRealStringStream( message, forces[pairArray[kk]], threeI );
          message << "]";
       }

       SimTKOpenMMLog::printMessage( message );
    }
668
    return ReferenceForce::DefaultReturn;
669
  }
670