ReferenceLJCoulomb14.cpp 4.2 KB
Newer Older
1

2
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
 * 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>

28
#include "SimTKOpenMMUtilities.h"
29
#include "ReferenceLJCoulomb14.h"
30
31
#include "ReferenceForce.h"

32
using std::vector;
33
using namespace OpenMM;
34

35
36
/**---------------------------------------------------------------------------------------

37
   ReferenceLJCoulomb14 constructor
38
39
40

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

41
ReferenceLJCoulomb14::ReferenceLJCoulomb14() : periodic(false) {
42
43
44
45
}

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

46
   ReferenceLJCoulomb14 destructor
47
48
49

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

50
ReferenceLJCoulomb14::~ReferenceLJCoulomb14() {
51
52
}

53
54
55
56
57
58
59
void ReferenceLJCoulomb14::setPeriodic(OpenMM::Vec3* vectors) {
    periodic = true;
    periodicBoxVectors[0] = vectors[0];
    periodicBoxVectors[1] = vectors[1];
    periodicBoxVectors[2] = vectors[2];
}

60
61
62
63
64
65
/**---------------------------------------------------------------------------------------

   Calculate LJ 1-4 ixn

   @param atomIndices      atom indices of 4 atoms in bond
   @param atomCoordinates  atom coordinates
66
67
   @param parameters       three parameters:
                                        parameters[0]= (c12/c6)**1/6  (sigma)
68
69
                                        parameters[1]= c6*c6/c12      (4*epsilon)
                                        parameters[2]= epsfac*q1*q2
70
   @param forces           force array (forces added to current values)
71
   @param totalEnergy      if not null, the energy will be added to this
72
73
74

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

75
76
void ReferenceLJCoulomb14::calculateBondIxn(vector<int>& atomIndices, vector<Vec3>& atomCoordinates,
                                     vector<double>& parameters, vector<Vec3>& forces,
peastman's avatar
peastman committed
77
                                     double* totalEnergy, double* energyParamDerivs) {
78
    double deltaR[2][ReferenceForce::LastDeltaRIndex];
79

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

82
83
84
85
86
87
    int atomAIndex = atomIndices[0];
    int atomBIndex = atomIndices[1];
    if (periodic)
        ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], periodicBoxVectors, deltaR[0]);
    else
        ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);  
88

89
90
91
92
    double inverseR  = 1.0/(deltaR[0][ReferenceForce::RIndex]);
    double sig2      = inverseR*parameters[0];
           sig2     *= sig2;
    double sig6      = sig2*sig2*sig2;
93

94
95
96
    double dEdR      = parameters[1]*(12.0*sig6 - 6.0)*sig6;
           dEdR     += ONE_4PI_EPS0*parameters[2]*inverseR;
           dEdR     *= inverseR*inverseR;
97

98
    // accumulate forces
99

100
101
102
103
104
    for (int ii = 0; ii < 3; ii++) {
        double force        = dEdR*deltaR[0][ii];
        forces[atomAIndex][ii] += force;
        forces[atomBIndex][ii] -= force;
    }
105

106
    // accumulate energies
107

108
109
    if (totalEnergy != NULL)
        *totalEnergy += parameters[1]*(sig6 - 1.0)*sig6 + (ONE_4PI_EPS0*parameters[2]*inverseR);
110
}