customHbondForce.cc 9.75 KB
Newer Older
1
/**
Peter Eastman's avatar
Peter Eastman committed
2
 * Compute the difference between two vectors, optionally taking periodic boundary conditions into account
3
4
 * and setting the fourth component to the squared magnitude.
 */
5
6
inline DEVICE real4 delta(real4 vec1, real4 vec2, real4 periodicBoxSize, real4 invPeriodicBoxSize, real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ) {
    real4 result = make_real4(vec1.x-vec2.x, vec1.y-vec2.y, vec1.z-vec2.z, 0);
7
#ifdef USE_PERIODIC
8
    APPLY_PERIODIC_TO_DELTA(result)
9
10
11
12
13
14
15
16
#endif
    result.w = result.x*result.x + result.y*result.y + result.z*result.z;
    return result;
}

/**
 * Compute the angle between two vectors.  The w component of each vector should contain the squared magnitude.
 */
17
inline DEVICE real computeAngle(real4 vec1, real4 vec2) {
18
19
20
    real dotProduct = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
    real cosine = dotProduct*RSQRT(vec1.w*vec2.w);
    real angle;
21
22
23
    if (cosine > 0.99f || cosine < -0.99f) {
        // We're close to the singularity in acos(), so take the cross product and use asin() instead.

24
        real3 crossProduct = cross(trimTo3(vec1), trimTo3(vec2));
25
        real scale = vec1.w*vec2.w;
26
27
28
        angle = ASIN(SQRT(dot(crossProduct, crossProduct)/scale));
        if (cosine < 0)
            angle = M_PI-angle;
29
30
    }
    else
31
       angle = ACOS(cosine);
32
33
34
35
    return angle;
}

/**
36
 * Compute the cross product of two vectors, setting the fourth component to the squared magnitude.
37
 */
38
39
40
inline DEVICE real4 computeCross(real4 vec1, real4 vec2) {
    real3 cp = cross(trimTo3(vec1), trimTo3(vec2));
    return make_real4(cp.x, cp.y, cp.z, cp.x*cp.x+cp.y*cp.y+cp.z*cp.z);
41
}
42

43
44
45
/**
 * Compute forces on donors.
 */
46
47
48
49
KERNEL void computeDonorForces(
	GLOBAL mm_ulong* RESTRICT force,
	GLOBAL mixed* RESTRICT energyBuffer, GLOBAL const real4* RESTRICT posq, GLOBAL const int4* RESTRICT exclusions,
        GLOBAL const int4* RESTRICT donorAtoms, GLOBAL const int4* RESTRICT acceptorAtoms, real4 periodicBoxSize, real4 invPeriodicBoxSize,
50
        real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ
51
        PARAMETER_ARGUMENTS) {
52
    LOCAL real4 posBuffer[3*THREAD_BLOCK_SIZE];
53
    mixed energy = 0;
54
55
56
57
    real3 f1 = make_real3(0);
    real3 f2 = make_real3(0);
    real3 f3 = make_real3(0);
    for (int donorStart = 0; donorStart < NUM_DONORS; donorStart += GLOBAL_SIZE) {
58
59
        // Load information about the donor this thread will compute forces on.

60
        int donorIndex = donorStart+GLOBAL_ID;
61
        int4 atoms, exclusionIndices;
62
        real4 d1, d2, d3;
63
64
        if (donorIndex < NUM_DONORS) {
            atoms = donorAtoms[donorIndex];
65
66
67
            d1 = (atoms.x > -1 ? posq[atoms.x] : make_real4(0));
            d2 = (atoms.y > -1 ? posq[atoms.y] : make_real4(0));
            d3 = (atoms.z > -1 ? posq[atoms.z] : make_real4(0));
68
69
70
#ifdef USE_EXCLUSIONS
            exclusionIndices = exclusions[donorIndex];
#endif
71
72
        }
        else
73
74
            atoms = make_int4(-1, -1, -1, -1);
        for (int acceptorStart = 0; acceptorStart < NUM_ACCEPTORS; acceptorStart += LOCAL_SIZE) {
75
76
            // Load the next block of acceptors into local memory.

77
78
79
80
81
82
83
84
85
            SYNC_THREADS;
            int blockSize = min((int) LOCAL_SIZE, NUM_ACCEPTORS-acceptorStart);
            if (LOCAL_ID < blockSize) {
                int4 atoms2 = acceptorAtoms[acceptorStart+LOCAL_ID];
                posBuffer[3*LOCAL_ID] = (atoms2.x > -1 ? posq[atoms2.x] : make_real4(0));
                posBuffer[3*LOCAL_ID+1] = (atoms2.y > -1 ? posq[atoms2.y] : make_real4(0));
                posBuffer[3*LOCAL_ID+2] = (atoms2.z > -1 ? posq[atoms2.z] : make_real4(0));
            }
            SYNC_THREADS;
86
87
            if (donorIndex < NUM_DONORS) {
                for (int index = 0; index < blockSize; index++) {
88
                    int acceptorIndex = acceptorStart+index;
89
#ifdef USE_EXCLUSIONS
90
91
92
                    if (acceptorIndex == exclusionIndices.x || acceptorIndex == exclusionIndices.y || acceptorIndex == exclusionIndices.z || acceptorIndex == exclusionIndices.w)
                        continue;
#endif
93
94
                    // Compute the interaction between a donor and an acceptor.

95
96
97
                    real4 a1 = posBuffer[3*index];
                    real4 a2 = posBuffer[3*index+1];
                    real4 a3 = posBuffer[3*index+2];
Peter Eastman's avatar
Peter Eastman committed
98
                    real4 deltaD1A1 = delta(d1, a1, periodicBoxSize, invPeriodicBoxSize, periodicBoxVecX, periodicBoxVecY, periodicBoxVecZ);
99
#ifdef USE_CUTOFF
100
                    if (deltaD1A1.w < CUTOFF_SQUARED) {
101
#endif
102
103
104
                        COMPUTE_DONOR_FORCE
#ifdef USE_CUTOFF
                    }
105
#endif
106
107
108
                }
            }
        }
109

110
        // Write results
111

112
        if (donorIndex < NUM_DONORS) {
113
            if (atoms.x > -1) {
114
115
116
                ATOMIC_ADD(&force[atoms.x], (mm_ulong) realToFixedPoint(f1.x));
                ATOMIC_ADD(&force[atoms.x+PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f1.y));
                ATOMIC_ADD(&force[atoms.x+2*PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f1.z));
117
118
119
                MEM_FENCE;
            }
            if (atoms.y > -1) {
120
121
122
                ATOMIC_ADD(&force[atoms.y], (mm_ulong) realToFixedPoint(f2.x));
                ATOMIC_ADD(&force[atoms.y+PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f2.y));
                ATOMIC_ADD(&force[atoms.y+2*PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f2.z));
123
124
125
                MEM_FENCE;
            }
            if (atoms.z > -1) {
126
127
128
                ATOMIC_ADD(&force[atoms.z], (mm_ulong) realToFixedPoint(f3.x));
                ATOMIC_ADD(&force[atoms.z+PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f3.y));
                ATOMIC_ADD(&force[atoms.z+2*PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f3.z));
129
130
                MEM_FENCE;
            }
131
132
        }
    }
133
    energyBuffer[GLOBAL_ID] += energy;
134
135
136
137
}
/**
 * Compute forces on acceptors.
 */
138
139
140
141
KERNEL void computeAcceptorForces(
	GLOBAL mm_ulong* RESTRICT force,
        GLOBAL mixed* RESTRICT energyBuffer, GLOBAL const real4* RESTRICT posq, GLOBAL const int4* RESTRICT exclusions,
        GLOBAL const int4* RESTRICT donorAtoms, GLOBAL const int4* RESTRICT acceptorAtoms, real4 periodicBoxSize, real4 invPeriodicBoxSize,
142
        real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ
143
        PARAMETER_ARGUMENTS) {
144
145
146
147
148
    LOCAL real4 posBuffer[3*THREAD_BLOCK_SIZE];
    real3 f1 = make_real3(0);
    real3 f2 = make_real3(0);
    real3 f3 = make_real3(0);
    for (int acceptorStart = 0; acceptorStart < NUM_ACCEPTORS; acceptorStart += GLOBAL_SIZE) {
149
150
        // Load information about the acceptor this thread will compute forces on.

151
        int acceptorIndex = acceptorStart+GLOBAL_ID;
152
        int4 atoms, exclusionIndices;
153
        real4 a1, a2, a3;
154
155
        if (acceptorIndex < NUM_ACCEPTORS) {
            atoms = acceptorAtoms[acceptorIndex];
156
157
158
            a1 = (atoms.x > -1 ? posq[atoms.x] : make_real4(0));
            a2 = (atoms.y > -1 ? posq[atoms.y] : make_real4(0));
            a3 = (atoms.z > -1 ? posq[atoms.z] : make_real4(0));
159
160
161
#ifdef USE_EXCLUSIONS
            exclusionIndices = exclusions[acceptorIndex];
#endif
162
163
        }
        else
164
165
            atoms = make_int4(-1, -1, -1, -1);
        for (int donorStart = 0; donorStart < NUM_DONORS; donorStart += LOCAL_SIZE) {
166
167
            // Load the next block of donors into local memory.

168
169
170
171
172
173
174
175
176
            SYNC_THREADS;
            int blockSize = min((int) LOCAL_SIZE, NUM_DONORS-donorStart);
            if (LOCAL_ID < blockSize) {
                int4 atoms2 = donorAtoms[donorStart+LOCAL_ID];
                posBuffer[3*LOCAL_ID] = (atoms2.x > -1 ? posq[atoms2.x] : make_real4(0));
                posBuffer[3*LOCAL_ID+1] = (atoms2.y > -1 ? posq[atoms2.y] : make_real4(0));
                posBuffer[3*LOCAL_ID+2] = (atoms2.z > -1 ? posq[atoms2.z] : make_real4(0));
            }
            SYNC_THREADS;
177
178
            if (acceptorIndex < NUM_ACCEPTORS) {
                for (int index = 0; index < blockSize; index++) {
179
                    int donorIndex = donorStart+index;
180
#ifdef USE_EXCLUSIONS
181
182
183
                    if (donorIndex == exclusionIndices.x || donorIndex == exclusionIndices.y || donorIndex == exclusionIndices.z || donorIndex == exclusionIndices.w)
                        continue;
#endif
184
185
                    // Compute the interaction between a donor and an acceptor.

186
187
188
                    real4 d1 = posBuffer[3*index];
                    real4 d2 = posBuffer[3*index+1];
                    real4 d3 = posBuffer[3*index+2];
Peter Eastman's avatar
Peter Eastman committed
189
                    real4 deltaD1A1 = delta(d1, a1, periodicBoxSize, invPeriodicBoxSize, periodicBoxVecX, periodicBoxVecY, periodicBoxVecZ);
190
191
#ifdef USE_CUTOFF
                    if (deltaD1A1.w < CUTOFF_SQUARED) {
192
#endif
193
                        COMPUTE_ACCEPTOR_FORCE
194
#ifdef USE_CUTOFF
195
                    }
196
#endif
197
                }
198
199
200
201
202
            }
        }

        // Write results

203
        if (acceptorIndex < NUM_ACCEPTORS) {
204
            if (atoms.x > -1) {
205
206
207
                ATOMIC_ADD(&force[atoms.x], (mm_ulong) realToFixedPoint(f1.x));
                ATOMIC_ADD(&force[atoms.x+PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f1.y));
                ATOMIC_ADD(&force[atoms.x+2*PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f1.z));
208
209
210
                MEM_FENCE;
            }
            if (atoms.y > -1) {
211
212
213
                ATOMIC_ADD(&force[atoms.y], (mm_ulong) realToFixedPoint(f2.x));
                ATOMIC_ADD(&force[atoms.y+PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f2.y));
                ATOMIC_ADD(&force[atoms.y+2*PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f2.z));
214
215
216
                MEM_FENCE;
            }
            if (atoms.z > -1) {
217
218
219
                ATOMIC_ADD(&force[atoms.z], (mm_ulong) realToFixedPoint(f3.x));
                ATOMIC_ADD(&force[atoms.z+PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f3.y));
                ATOMIC_ADD(&force[atoms.z+2*PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(f3.z));
220
221
                MEM_FENCE;
            }
222
        }
223
224
    }
}