Skip to content

GitLab

Commit 8c3878b2 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Mods

parent 116a99ba
Hide whitespace changes
Inline Side-by-side
Showing with 974 additions and 1250 deletions
platforms/brook/src/gpu/kObcBaseD2.br → platforms/brook/src/gpu/kgbsa.br
View file @ 8c3878b2
/****************************************************************
* This file is part of the gpu acceleration library for gromacs.
* Author: Mark Friedrichs
*
* This kernel was developed in collaboration with
*
* Copyright (C) Pande Group, Stanford, 2006
*****************************************************************/
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs, Chris Bruns *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
//124 FLOPS [??? total]
kernel void bornSumInternal( float3 d1, float3 d2, float3 d3, float3 d4, float4 jAtomicRadiiScaled,
float iAtomicRadii, out float4 bornSum<> ){
// ---------------------------------------------------------------------------------------
float4 r, rInverse, r2, r2Inverse;
float4 l_ij, l_ij2;
float4 u_ij, u_ij2;
float4 t3;
float4 logRatio, scaledRadiusJ2, rPlusScaledJ, rMinusScaledJ, rMinusScaledJAbs;
float4 iAtomicRadii4;
float4 diff_u2_l2;
const float bigValue = 1000000.0f;
const float4 bigValue4 = float4( bigValue, bigValue, bigValue, bigValue );
const float smallValue = 0.000001f;
const float4 smallValue4 = float4( smallValue, smallValue, smallValue, smallValue );
const float4 zero4 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
const float4 one4 = float4( 1.0f, 1.0f, 1.0f, 1.0f );
float4 mask;
// ---------------------------------------------------------------------------------------
iAtomicRadii4 = float4( iAtomicRadii, iAtomicRadii, iAtomicRadii, iAtomicRadii );
r2 = float4( dot(d1, d1), dot( d2, d2 ), dot( d3, d3 ), dot( d4, d4 ) ); //20
// ---------------------------------------------------------------------------------------
// not needed for force calculation since deltaX = deltaY = deltaZ = 0, if r = 0
// r2 not masked -- assume for i = j, that 0 * infinity i 0; appears to work
mask = r2 < smallValue4 ? zero4 : one4;
// ---------------------------------------------------------------------------------------
r = sqrt( r2 ); //4
rInverse = rsqrt( r2 ); //8
rPlusScaledJ = r + jAtomicRadiiScaled; //4
rMinusScaledJ = r - jAtomicRadiiScaled; //4
rMinusScaledJAbs = abs( rMinusScaledJ );
//40
// ---------------------------------------------------------------------------------------
mask = iAtomicRadii4 >= rPlusScaledJ ? zero4 : mask;
l_ij = iAtomicRadii4 > rMinusScaledJAbs ? iAtomicRadii4 : rMinusScaledJAbs;
// ---------------------------------------------------------------------------------------
l_ij = 1.0f/l_ij; //4
l_ij2 = l_ij*l_ij; //4
u_ij = 1.0f/rPlusScaledJ; //4
u_ij2 = u_ij*u_ij; //4
diff_u2_l2 = u_ij2 - l_ij2; //4
logRatio = log(u_ij/l_ij); //4
//24
// ---------------------------------------------------------------------------------------
// terms associated w/ dL/dr & dU/dr are zero and not included
// Born radius term
bornSum = l_ij - u_ij + 0.5f*rInverse*logRatio + 0.25f*( r - jAtomicRadiiScaled*jAtomicRadiiScaled*rInverse)*diff_u2_l2; //10*4
t3 = ( iAtomicRadii4 < (jAtomicRadiiScaled - r) ) ? 2.0f*( (1.0f/iAtomicRadii) - l_ij) : zero4; //12 for true, but don't know how to average so saying 0
bornSum += t3; //4
bornSum *= mask; //4
//60
}
//???? FLOPS
kernel void kCalculateBornRadii( float numberOfAtoms, float roundedUpAtoms, float duplicationFactor,
float streamWidth, float fstreamWidth, float3 posq[][],
float2 scaledAtomicRadii[][],
out float4 bornForce<> ){
// ---------------------------------------------------------------------------------------
// atomic radii: radius - dielectric offset
float i1AtomicRadii, i2AtomicRadii, i3AtomicRadii, i4AtomicRadii;
// scaled atomic radii: scaleFactor*( radius - dielectric offset )
float j1ScaledAtomicRadii, j2ScaledAtomicRadii, j3ScaledAtomicRadii, j4ScaledAtomicRadii;
float4 jScaledAtomicRadii;
// (Born radii**2)*BornForce*ObcChainForce
// i,j coordinates
float3 i1pos, i2pos, i3pos, i4pos;
float3 j1pos, j2pos, j3pos, j4pos;
// delta coordinates
float3 d1, d2, d3, d4;
// indices
float2 iAtom;
float forceIndex;
//This is forceIndex mod numberOfAtoms, the true i index
float4 de, bornSum;
float iAtomLinearIndex, jLinind;
float2 jAtom;
float jEnd, jStart, jBlock;
float whichRep;
float tmp;
float tmp2;
// ---------------------------------------------------------------------------------------
// constants
const float I_Unroll = 4.0f;
// ---------------------------------------------------------------------------------------
// set linear index
iAtom = indexof( bornForce );
forceIndex = I_Unroll*( iAtom.x + iAtom.y*fstreamWidth );
iAtomLinearIndex = fmod( forceIndex, roundedUpAtoms );
// ---------------------------------------------------------------------------------------
// set gather coordinates
iAtom.x = fmod( iAtomLinearIndex, streamWidth );
iAtom.y = round( (iAtomLinearIndex - iAtom.x)/streamWidth );
// ---------------------------------------------------------------------------------------
// fetch position, (atomic radius - dielectric offset )
i1pos = posq[ iAtom ];
i1AtomicRadii = scaledAtomicRadii[ iAtom ].y;
iAtom.x += 1;
i2pos = posq[ iAtom ];
i2AtomicRadii = scaledAtomicRadii[ iAtom ].y;
iAtom.x += 1;
i3pos = posq[ iAtom ];
i3AtomicRadii = scaledAtomicRadii[ iAtom ].y;
iAtom.x += 1;
i4pos = posq[ iAtom ];
i4AtomicRadii = scaledAtomicRadii[ iAtom ].y;
bornForce = float4( 0.0f, 0.0f, 0.0f, 0.0f );
// inner loop indices setup
//changed the following instruction for rounding issues on some ASICs
//whichRep = floor( forceIndex / roundedUpAtoms );
tmp = fmod(forceIndex, roundedUpAtoms);
whichRep = round((forceIndex - tmp)/roundedUpAtoms);
jBlock = 1 + round( (numberOfAtoms-fmod( numberOfAtoms, (duplicationFactor*streamWidth) )) / (duplicationFactor*streamWidth) );
jStart = whichRep*jBlock;
jEnd = ( whichRep > duplicationFactor - 1.5f ) ? 999999.0f : (jStart + jBlock);
jAtom.y = jStart;
jLinind = jAtom.y*streamWidth;
// ---------------------------------------------------------------------------------------
while ( jAtom.y < jEnd && ( numberOfAtoms - jLinind ) > 0.9f ){
jAtom.x = 0.0f;
while ( jAtom.x < streamWidth && ( numberOfAtoms - jLinind ) > 0.9f ){
// ---------------------------------------------------------------------------------------
// gather required values
j1ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j1pos = posq[ jAtom ];
jAtom.x += 1.0f;
j2ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j2pos = posq[ jAtom ];
jAtom.x += 1.0f;
j3ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j3pos = posq[ jAtom ];
jAtom.x += 1.0f;
j4ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j4pos = posq[ jAtom ];
jAtom.x += 1.0f;
jScaledAtomicRadii = float4( j1ScaledAtomicRadii, j2ScaledAtomicRadii, j3ScaledAtomicRadii, j4ScaledAtomicRadii );
// ---------------------------------------------------------------------------------------
// First set of 4 -- i == 1
d1 = j1pos - i1pos; //3
d2 = j2pos - i1pos; //3
d3 = j3pos - i1pos; //3
d4 = j4pos - i1pos; //3
// i = 1
bornSumInternal( d1, d2, d3, d4, jScaledAtomicRadii, i1AtomicRadii, bornSum ); //112 : 368(?)
bornForce.x += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// ---------------------------------------------------------------------------------------
// i = 2
d1 = j1pos - i2pos; //3
d2 = j2pos - i2pos; //3
d3 = j3pos - i2pos; //3
d4 = j4pos - i2pos; //3
bornSumInternal( d1, d2, d3, d4, jScaledAtomicRadii, i2AtomicRadii, bornSum );//112 : 368(?)
bornForce.y += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// ---------------------------------------------------------------------------------------
// i = 3
d1 = j1pos - i3pos; //3
d2 = j2pos - i3pos; //3
d3 = j3pos - i3pos; //3
d4 = j4pos - i3pos; //3
bornSumInternal( d1, d2, d3, d4, jScaledAtomicRadii, i3AtomicRadii, bornSum );//112 : 368(?)
bornForce.z += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// ---------------------------------------------------------------------------------------
// i = 4
d1 = j1pos - i4pos; //3
d2 = j2pos - i4pos; //3
d3 = j3pos - i4pos; //3
d4 = j4pos - i4pos; //3
bornSumInternal( d1, d2, d3, d4, jScaledAtomicRadii, i4AtomicRadii, bornSum );//112 : 368(?)
bornForce.w += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// ---------------------------------------------------------------------------------------
// increment linear index by 4
jLinind += 4.0f;
}
jAtom.y += 1.0f;
}
}
kernel void kPostCalculateBornRadii_nobranch(
float repfac,
float atomStreamWidth,
float pStreamWidth,
float natoms,
float roundNatoms,
float iUnroll,
float conversion,
float mergeNonObcForces,
float4 pstream1[][],
float atomicRadii<>,
out float bornRadii<>,
out float obcChain<> ){
// ---------------------------------------------------------------------------------------
float atomIndex, forceIndex, qIndex, qOff;
float2 pindex;
float i;
float sum2, sum3, bornSum, tanhSum, atomicRadiiOffset, obcIntermediate;
float4 o1;
float tmp;
float2 iAtom;
float4 forces;
float expPlus, expMinus;
// ---------------------------------------------------------------------------------------
// constants -- OBC Type II
const float alphaObc = 1.0f;
const float betaObc = 0.8f;
const float gammaObc = 4.85f;
const float dielectricOffset = 0.009f;
// ---------------------------------------------------------------------------------------
// given atom index find force indices and streams
pindex = indexof( obcChain );
atomIndex = pindex.x + pindex.y*atomStreamWidth;
forceIndex = atomIndex;
// add current forces in inStream to forces stored in pstreams
// the .w entry is Born sum values; it will be used to calculate the
// Born radii and obcChain term
bornSum = 0.0f;
// sum over j-loop 'duplications' by gathering from pstreams
for( i = 0.0f; i < repfac; i += 1.0f ){
qIndex = round( (forceIndex - fmod( forceIndex, iUnroll))/iUnroll );
qOff = forceIndex - iUnroll*qIndex;
pindex.y = round( (qIndex - fmod( qIndex, pStreamWidth ))/pStreamWidth );
pindex.x = qIndex - pindex.y*pStreamWidth;
o1 = pstream1[ pindex ];
tmp = qOff < 0.5f ? o1.x : o1.y;
tmp = qOff < 1.5f ? tmp : o1.z;
tmp = qOff < 2.5f ? tmp : o1.w;
bornSum += tmp;
forceIndex += roundNatoms;
}
// compute Born radii and ObcChain
atomicRadiiOffset = atomicRadii - dielectricOffset;
bornSum *= 0.5f*atomicRadiiOffset;
sum2 = bornSum*bornSum;
sum3 = bornSum*sum2;
// Tanh does not exist?
// calculate [ exp(x) - exp(-x) ]/[ exp(x) + exp(-x) ]
// tanhSum = tanh( bornSum - betaObc*sum2 + gammaObc*sum3 );
tanhSum = bornSum - betaObc*sum2 + gammaObc*sum3;
expPlus = exp( tanhSum );
expMinus = 1.0f/expPlus;
tanhSum = ( expPlus - expMinus )/( expPlus + expMinus );
bornRadii = 1.0f/( (1.0f/(atomicRadiiOffset)) - tanhSum/atomicRadii );
obcIntermediate = atomicRadiiOffset*( alphaObc - 2.0f*betaObc*bornSum + 3.0f*gammaObc*sum2 );
obcChain = (1.0f - tanhSum*tanhSum)*obcIntermediate/atomicRadii;
if( atomIndex >= natoms ){
bornRadii = 0.0f;
obcChain = 0.0f;
}
}
kernel void loop1Internal( float3 d1, float3 d2, float3 d3, float3 d4, float4 jBornR,
float4 jQ, float iBornR, float iQ, out float4 dGpol_dr<>,
out float4 dGpol_dalpha2_ij<> ){
// ---------------------------------------------------------------------------------------
float4 r2, alpha2_ij, D_ij, expTerm, denominator2, denominator, Gpol;
// ---------------------------------------------------------------------------------------
r2 = float4( dot(d1, d1), dot( d2, d2 ), dot( d3, d3 ), dot( d4, d4 ) );
alpha2_ij = jBornR*iBornR;
D_ij = r2/(4.0f*alpha2_ij);
expTerm = exp( -D_ij );
denominator2 = r2 + alpha2_ij*expTerm;
denominator = sqrt( denominator2 );
Gpol = jQ/denominator;
Gpol *= iQ;
dGpol_dr = -Gpol*( 1.0f - 0.25f*expTerm )/denominator2;
dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*( 1.0f + D_ij )*jBornR/denominator2;
}
/* ---------------------------------------------------------------------------------------
Calculate nonpolar ACE term (Simbios)
bornRadius: Born radius
vdwRadius: Vdw radius
duplicationFactor: duplication factor
aceForce: ACE term
--------------------------------------------------------------------------------------- */
kernel void kAceNonPolarLoop1( float iBornRadius, float iVdwRadius, float duplicationFactor,
out float aceForce<> ){
// ---------------------------------------------------------------------------------------
// nonpolar term
float iSurface;
float iAceTerm;
// ---------------------------------------------------------------------------------------
// constants
// solvent radius
const float probeRadius = 0.14f;
// PI*4*6*0.0054*1000 (0.0054=asolv from Tinker)
//const float PI_24_aI = -0.3694512961;
const float PI_24_aI = -407.1504079f;
// ---------------------------------------------------------------------------------------
// etch i position and partial charge
// e = ai * term * (ri+probe)**2 * (ri/rb)**6
// (drbi) = drb(i) - 6.0fd0*e/rb
// (rI+probe)**2
iSurface = (iVdwRadius+probeRadius);
iSurface = iSurface*iSurface;
// (rI/rB)**6
iAceTerm = iVdwRadius/iBornRadius;
iAceTerm = iAceTerm*iAceTerm*iAceTerm;
iAceTerm = iAceTerm*iAceTerm;
aceForce = iSurface*iAceTerm*PI_24_aI/(duplicationFactor*iBornRadius);
}
/* ---------------------------------------------------------------------------------------
Calculate first loop force terms (Simbios)
numberOfAtoms: no. of atoms
roundedUpAtoms: rounded up number of atoms -- accounts for unrolling
duplicationFactor: number of threads for inner loop
streamWidth: atom stream width
fstreamWidth: force stream width (output -- i-unroll)
soluteDielectric: solute dielectric
solventDielectric: solvent dielectric
includeAce: include ACE term
posq: atom positions and charge
bornRadii: Born radii
nonpolarForce: nonpolar force (0 if nonpolar not included, else
ACE value)
bornForce1: i-unroll first force component, including dBornR/dr in .w
bornForce2: i-unroll second force component, including dBornR/dr in .w
bornForce3: i-unroll first force component, including dBornR/dr in .w
bornForce4: i-unroll second force component, including dBornR/dr in .w
--------------------------------------------------------------------------------------- */
kernel void kObcLoop1( float numberOfAtoms, float roundedUpAtoms, float duplicationFactor,
float streamWidth, float fstreamWidth, float soluteDielectric,
float solventDielectric, float includeAce,
float3 posq[][], float bornRadii[][], float2 atomicRadii[][],
out float4 bornForce1<>, out float4 bornForce2<>,
out float4 bornForce3<>, out float4 bornForce4<> ){
// ---------------------------------------------------------------------------------------
// Born radii
float i1BornR, i2BornR, i3BornR, i4BornR;
float j1BornR, j2BornR, j3BornR, j4BornR;
float4 jBornR;
// atomic radii
float i1AtomicR, i2AtomicR, i3AtomicR, i4AtomicR;
// i,j coordinates
float3 i1Pos, i2Pos, i3Pos, i4Pos;
float3 j1Pos, j2Pos, j3Pos, j4Pos;
float4 j1PosQ, j2PosQ, j3PosQ, j4PosQ;
// i, j partial charges
float i1Q, i2Q, i3Q, i4Q;
float j1Q, j2Q, j3Q, j4Q;
float4 jQ;
float aceForce;
// delta coordinates
float3 d1, d2, d3, d4;
// intermediate terms
float4 dGpol_dr, dGpol_dalpha2_ij;
// indices
float2 iAtom;
float forceIndex;
// This is forceIndex mod numberOfAtoms, the true i index
float iAtomLinearIndex, jLinind;
float2 jAtom;
float jEnd, jStart, jBlock;
float whichRep;
float tmp;
// ---------------------------------------------------------------------------------------
// electricConstant = -166.0f2691;
// preFactor = 2.0f*electricConstant*(1.0f - (1.0f/waterDielectric))
float preFactor = -332.05382f;
const float I_Unroll = 4.0f;
const float3 zero3 = float3( 0.0f, 0.0f, 0.0f );
// ---------------------------------------------------------------------------------------
preFactor *= ( (1.0f/soluteDielectric) - (1.0f/solventDielectric) );
iAtom = indexof( bornForce1 );
forceIndex = I_Unroll*( iAtom.x + iAtom.y*fstreamWidth );
iAtomLinearIndex = fmod( forceIndex, roundedUpAtoms );
// ---------------------------------------------------------------------------------------
// set gather index
iAtom.x = fmod( iAtomLinearIndex, streamWidth );
iAtom.y = round( (iAtomLinearIndex - fmod(iAtomLinearIndex, streamWidth ))/streamWidth );
// ---------------------------------------------------------------------------------------
// etch i1 position and partial charge
jQ = posq[ iAtom ];
i1Pos = jQ.xyz;
i1Q = atomicRadii[ iAtom ].y;
i1Q *= preFactor;
i1BornR = bornRadii[ iAtom ];
i1AtomicR = atomicRadii[ iAtom ].x;
kAceNonPolarLoop1( i1BornR, i1AtomicR, duplicationFactor, aceForce );
bornForce1.xyz = zero3;
bornForce1.w = includeAce > 0.5f ? aceForce : 0.0f;
// ---------------------------------------------------------------------------------------
// etch i2 position and partial charge
iAtom.x += 1;
jQ = posq[ iAtom ];
i2Pos = jQ.xyz;
i2Q = atomicRadii[ iAtom ].y;
i2Q *= preFactor;
i2BornR = bornRadii[ iAtom ];
i2AtomicR = atomicRadii[ iAtom ].x;
kAceNonPolarLoop1( i2BornR, i2AtomicR, duplicationFactor, aceForce );
bornForce2.xyz = zero3;
bornForce2.w = includeAce > 0.5f ? aceForce : 0.0f;
// ---------------------------------------------------------------------------------------
// etch i3 position and partial charge
iAtom.x += 1;
jQ = posq[ iAtom ];
i3Pos = jQ.xyz;
i3Q = atomicRadii[ iAtom ].y;
i3Q *= preFactor;
i3BornR = bornRadii[ iAtom ];
i3AtomicR = atomicRadii[ iAtom ].x;
kAceNonPolarLoop1( i3BornR, i3AtomicR, duplicationFactor, aceForce );
bornForce3.xyz = zero3;
bornForce3.w = includeAce > 0.5f ? aceForce : 0.0f;
// ---------------------------------------------------------------------------------------
// etch i4 position and partial charge
iAtom.x += 1;
jQ = posq[ iAtom ];
i4Pos = jQ.xyz;
i4Q = atomicRadii[ iAtom ].y;
i4Q *= preFactor;
i4BornR = bornRadii[ iAtom ];
i4AtomicR = atomicRadii[ iAtom ].x;
kAceNonPolarLoop1( i4BornR, i4AtomicR, duplicationFactor, aceForce );
bornForce4.xyz = zero3;
bornForce4.w = includeAce > 0.5f ? aceForce : 0.0f;
// ---------------------------------------------------------------------------------------
// inner loop setup
// if dupFac == 4, I_UnRoll =2, then breaking inner loop into two segments
// to increase number of threads in flight
// forceStreamSz = N*RepFac/I_UnRoll
// forceIndex = I_UnRoll*( a.x + a.y*forceStreamSz )
// whichRep = 0 or 1
// jBlock = 1 + floor[ N/(duplicationFactor*streamWidth) ]
//changed the following instruction for rounding issues on some ASICs
//whichRep = floor( forceIndex / roundedUpAtoms );
tmp = fmod(forceIndex, roundedUpAtoms);
whichRep = round((forceIndex - tmp)/roundedUpAtoms);
jBlock = 1 + floor( numberOfAtoms/(duplicationFactor*streamWidth ) );
jStart = whichRep*jBlock;
jEnd = ( whichRep > duplicationFactor - 1.5f ) ? 999999.0f : (jStart + jBlock);
jAtom.y = jStart;
jLinind = jAtom.y*streamWidth;
// ---------------------------------------------------------------------------------------
while ( jAtom.y < jEnd && ( numberOfAtoms - jLinind ) > 0.9f ){
jAtom.x = 0.0f;
while ( jAtom.x < streamWidth && ( numberOfAtoms - jLinind ) > 0.9f ) {
// ---------------------------------------------------------------------------------------
// gather required values
j1Pos = posq[ jAtom ];
j1Q = atomicRadii[ jAtom ].y;
j1BornR = bornRadii[ jAtom ];
jAtom.x += 1.0f;
j2Pos = posq[ jAtom ];
j2Q = atomicRadii[ jAtom ].y;
j2BornR = bornRadii[ jAtom ];
jAtom.x += 1.0f;
j3Pos = posq[ jAtom ];
j3Q = atomicRadii[ jAtom ].y;
j3BornR = bornRadii[ jAtom ];
jAtom.x += 1.0f;
j4Pos = posq[ jAtom ];
j4Q = atomicRadii[ jAtom ].y;
j4BornR = bornRadii[ jAtom ];
jAtom.x += 1.0f;
jBornR = float4( j1BornR, j2BornR, j3BornR, j4BornR );
jQ = float4( j1Q, j2Q, j3Q, j4Q );
// ---------------------------------------------------------------------------------------
// i == 1
d1 = i1Pos - j1Pos;
d2 = i1Pos - j2Pos;
d3 = i1Pos - j3Pos;
d4 = i1Pos - j4Pos;
loop1Internal( d1, d2, d3, d4, jBornR, jQ, i1BornR, i1Q, dGpol_dr, dGpol_dalpha2_ij );
bornForce1.xyz += dGpol_dr.x*d1;
bornForce1.xyz += dGpol_dr.y*d2;
bornForce1.xyz += dGpol_dr.z*d3;
bornForce1.xyz += dGpol_dr.w*d4;
bornForce1.w += dGpol_dalpha2_ij.x + dGpol_dalpha2_ij.y + dGpol_dalpha2_ij.z + dGpol_dalpha2_ij.w;
// ---------------------------------------------------------------------------------------
// i == 2
d1 = i2Pos - j1Pos;
d2 = i2Pos - j2Pos;
d3 = i2Pos - j3Pos;
d4 = i2Pos - j4Pos;
loop1Internal( d1, d2, d3, d4, jBornR, jQ, i2BornR, i2Q, dGpol_dr, dGpol_dalpha2_ij );
bornForce2.xyz += dGpol_dr.x*d1;
bornForce2.xyz += dGpol_dr.y*d2;
bornForce2.xyz += dGpol_dr.z*d3;
bornForce2.xyz += dGpol_dr.w*d4;
bornForce2.w += dGpol_dalpha2_ij.x + dGpol_dalpha2_ij.y + dGpol_dalpha2_ij.z + dGpol_dalpha2_ij.w;
// ---------------------------------------------------------------------------------------
// i == 3
d1 = i3Pos - j1Pos;
d2 = i3Pos - j2Pos;
d3 = i3Pos - j3Pos;
d4 = i3Pos - j4Pos;
loop1Internal( d1, d2, d3, d4, jBornR, jQ, i3BornR, i3Q, dGpol_dr, dGpol_dalpha2_ij );
bornForce3.xyz += dGpol_dr.x*d1;
bornForce3.xyz += dGpol_dr.y*d2;
bornForce3.xyz += dGpol_dr.z*d3;
bornForce3.xyz += dGpol_dr.w*d4;
bornForce3.w += dGpol_dalpha2_ij.x + dGpol_dalpha2_ij.y + dGpol_dalpha2_ij.z + dGpol_dalpha2_ij.w;
// ---------------------------------------------------------------------------------------
// i == 4
d1 = i4Pos - j1Pos;
d2 = i4Pos - j2Pos;
d3 = i4Pos - j3Pos;
d4 = i4Pos - j4Pos;
loop1Internal( d1, d2, d3, d4, jBornR, jQ, i4BornR, i4Q, dGpol_dr, dGpol_dalpha2_ij );
bornForce4.xyz += dGpol_dr.x*d1;
bornForce4.xyz += dGpol_dr.y*d2;
bornForce4.xyz += dGpol_dr.z*d3;
bornForce4.xyz += dGpol_dr.w*d4;
bornForce4.w += dGpol_dalpha2_ij.x + dGpol_dalpha2_ij.y + dGpol_dalpha2_ij.z + dGpol_dalpha2_ij.w;
// ---------------------------------------------------------------------------------------
jLinind += 4.0f;
}
jAtom.y += 1.0f;
}
}
// The inner loop by definition creates divergent paths. Chances are
// fair that we will take all sides of the branch anyway, so this
// verion uses manual predication
kernel void kPostObcLoop1_nobranch(
float repfac,
float atomStreamWidth,
float pStreamWidth,
float natoms,
float roundNatoms,
float iUnroll,
float4 pstream1[][],
float4 pstream2[][],
float4 pstream3[][],
float4 pstream4[][],
float obcChain<>,
float bornRadii<>,
out float4 outstream<>,
out float bornRadii2Force<> )
{
// ---------------------------------------------------------------------------------------
float atomIndex, forceIndex, qIndex, qOff;
float2 pindex;
float i;
float4 o1,o2,o3,o4;
float4 tmp;
// given atom index find force indices and streams
pindex = indexof( outstream );
atomIndex = pindex.x + pindex.y*atomStreamWidth;
forceIndex = atomIndex;
outstream = float4( 0.0f, 0.0f, 0.0f, 0.0f );
for( i = 0.0f; i < repfac; i += 1.0f ){
// qIndex = floor( forceIndex/iUnroll );
qIndex = round( (forceIndex - fmod( forceIndex, iUnroll))/iUnroll );
qOff = forceIndex - iUnroll*qIndex;
// pindex.y = floor( qIndex/ pStreamWidth );
pindex.y = round( (qIndex - fmod( qIndex, pStreamWidth ))/pStreamWidth );
// pindex.x = qIndex - pindex.y*pStreamWidth + qOff;
pindex.x = qIndex - pindex.y*pStreamWidth;
o1 = pstream1[ pindex ];
o2 = pstream2[ pindex ];
o3 = pstream3[ pindex ];
o4 = pstream4[ pindex ];
tmp = qOff < 0.5f ? o1 : o2;
tmp = qOff < 1.5f ? tmp : o3;
tmp = qOff < 2.5f ? tmp : o4;
outstream += tmp;
forceIndex += roundNatoms;
}
bornRadii2Force = obcChain*bornRadii*bornRadii*outstream.w;
}
//156 FLOPS [172 total]
kernel void loop2Internal( float3 d1, float3 d2, float3 d3, float3 d4, float4 jAtomicRadiiScaled,
......@@ -260,90 +1109,10 @@ kernel void loop2InternalNoSum( float3 d1, float3 d2, float3 d3, float3 d4, floa
// terms associated w/ dL/dr & dU/dr are zero and not included
t3 = 0.125f*(1.0f + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + 0.25f*logRatio*r2Inverse; //8*4
de = mask*bornForce*t3*rInverse; //12
}
//124 FLOPS [??? total]
kernel void bornSumInternal( float3 d1, float3 d2, float3 d3, float3 d4, float4 jAtomicRadiiScaled,
float iAtomicRadii, out float4 bornSum<> ){
// ---------------------------------------------------------------------------------------
float4 r, rInverse, r2, r2Inverse;
float4 l_ij, l_ij2;
float4 u_ij, u_ij2;
float4 t3;
float4 logRatio, scaledRadiusJ2, rPlusScaledJ, rMinusScaledJ, rMinusScaledJAbs;
float4 iAtomicRadii4;
float4 diff_u2_l2;
const float bigValue = 1000000.0f;
const float4 bigValue4 = float4( bigValue, bigValue, bigValue, bigValue );
const float smallValue = 0.000001f;
const float4 smallValue4 = float4( smallValue, smallValue, smallValue, smallValue );
const float4 zero4 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
const float4 one4 = float4( 1.0f, 1.0f, 1.0f, 1.0f );
float4 mask;
// ---------------------------------------------------------------------------------------
iAtomicRadii4 = float4( iAtomicRadii, iAtomicRadii, iAtomicRadii, iAtomicRadii );
r2 = float4( dot(d1, d1), dot( d2, d2 ), dot( d3, d3 ), dot( d4, d4 ) ); //20
// ---------------------------------------------------------------------------------------
// not needed for force calculation since deltaX = deltaY = deltaZ = 0, if r = 0
// r2 not masked -- assume for i = j, that 0 * infinity i 0; appears to work
mask = r2 < smallValue4 ? zero4 : one4;
// ---------------------------------------------------------------------------------------
r = sqrt( r2 ); //4
rInverse = rsqrt( r2 ); //8
rPlusScaledJ = r + jAtomicRadiiScaled; //4
rMinusScaledJ = r - jAtomicRadiiScaled; //4
rMinusScaledJAbs = abs( rMinusScaledJ );
//40
// ---------------------------------------------------------------------------------------
mask = iAtomicRadii4 >= rPlusScaledJ ? zero4 : mask;
l_ij = iAtomicRadii4 > rMinusScaledJAbs ? iAtomicRadii4 : rMinusScaledJAbs;
// ---------------------------------------------------------------------------------------
l_ij = 1.0f/l_ij; //4
l_ij2 = l_ij*l_ij; //4
u_ij = 1.0f/rPlusScaledJ; //4
u_ij2 = u_ij*u_ij; //4
diff_u2_l2 = u_ij2 - l_ij2; //4
logRatio = log(u_ij/l_ij); //4
//24
// ---------------------------------------------------------------------------------------
// terms associated w/ dL/dr & dU/dr are zero and not included
// Born radius term
bornSum = l_ij - u_ij + 0.5f*rInverse*logRatio + 0.25f*( r - jAtomicRadiiScaled*jAtomicRadiiScaled*rInverse)*diff_u2_l2; //10*4
t3 = 0.125f*(1.0f + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + 0.25f*logRatio*r2Inverse; //8*4
t3 = ( iAtomicRadii4 < (jAtomicRadiiScaled - r) ) ? 2.0f*( (1.0f/iAtomicRadii) - l_ij) : zero4; //12 for true, but don't know how to average so saying 0
bornSum += t3; //4
bornSum *= mask; //4
de = mask*bornForce*t3*rInverse; //12
//60
}
/* ---------------------------------------------------------------------------------------
......@@ -370,9 +1139,8 @@ kernel void bornSumInternal( float3 d1, float3 d2, float3 d3, float3 d4, float4
//1376 FLOPS
kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicationFactor,
float streamWidth, float fstreamWidth, float4 posq[][],
float atomicRadii[][], float scaledAtomicRadii[][],
float bornForceFactor[][],
float streamWidth, float fstreamWidth, float3 posq[][],
float scaledAtomicRadii[][], float bornForceFactor[][],
out float4 bornForce1<>, out float4 bornForce2<>,
out float4 bornForce3<>, out float4 bornForce4<> ){
......@@ -451,37 +1219,41 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
// etch i1 position, Born prefactor, atomic radius
i1pos = posq[ iAtom ].xyz;
i1pos = posq[ iAtom ];
i1BornForceFactor = bornForceFactor[ iAtom ];
i1AtomicRadii = atomicRadii[ iAtom ];
i1ScaledAtomicRadii = scaledAtomicRadii[ iAtom ];
//i1AtomicRadii = atomicRadii[ iAtom ];
i1AtomicRadii = scaledAtomicRadii[ iAtom ].y;
i1ScaledAtomicRadii = scaledAtomicRadii[ iAtom ].x;
bornForce1 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
// etch i2 position, Born prefactor, atomic radius
iAtom.x += 1;
i2pos = posq[ iAtom ].xyz;
i2pos = posq[ iAtom ];
i2BornForceFactor = bornForceFactor[ iAtom ];
i2AtomicRadii = atomicRadii[ iAtom ];
i2ScaledAtomicRadii = scaledAtomicRadii[ iAtom ];
//i2AtomicRadii = atomicRadii[ iAtom ];
i2AtomicRadii = scaledAtomicRadii[ iAtom ].y;
i2ScaledAtomicRadii = scaledAtomicRadii[ iAtom ].x;
bornForce2 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
// etch i3 position, Born prefactor, atomic radius
iAtom.x += 1;
i3pos = posq[ iAtom ].xyz;
i3pos = posq[ iAtom ];
i3BornForceFactor = bornForceFactor[ iAtom ];
i3AtomicRadii = atomicRadii[ iAtom ];
i3ScaledAtomicRadii = scaledAtomicRadii[ iAtom ];
//i3AtomicRadii = atomicRadii[ iAtom ];
i3AtomicRadii = scaledAtomicRadii[ iAtom ].y;
i3ScaledAtomicRadii = scaledAtomicRadii[ iAtom ].x;
bornForce3 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
// etch i4 position, Born prefactor, atomic radius
iAtom.x += 1;
i4pos = posq[ iAtom ].xyz;
i4pos = posq[ iAtom ];
i4BornForceFactor = bornForceFactor[ iAtom ];
i4AtomicRadii = atomicRadii[ iAtom ];
i4ScaledAtomicRadii = scaledAtomicRadii[ iAtom ];
//i4AtomicRadii = atomicRadii[ iAtom ];
i4AtomicRadii = scaledAtomicRadii[ iAtom ].y;
i4ScaledAtomicRadii = scaledAtomicRadii[ iAtom ].x;
bornForce4 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
// create float4 for main vars
......@@ -493,11 +1265,11 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
// inner loop indices setup
//changed the following instruction for rounding issues on some ASICs
//whichRep = floor( forceIndex / roundedUpAtoms );
tmp = fmod(forceIndex, roundedUpAtoms);
whichRep = round((forceIndex - tmp)/roundedUpAtoms);
tmp = fmod(forceIndex, roundedUpAtoms);
whichRep = round((forceIndex - tmp)/roundedUpAtoms);
//jBlock = 1 + floor( numberOfAtoms/(duplicationFactor*streamWidth ) );
jBlock = 1 + round( (numberOfAtoms-fmod( numberOfAtoms, (duplicationFactor*streamWidth) )) / (duplicationFactor*streamWidth) );
jBlock = 1 + round( (numberOfAtoms-fmod( numberOfAtoms, (duplicationFactor*streamWidth) )) / (duplicationFactor*streamWidth) );
jStart = whichRep*jBlock;
......@@ -518,27 +1290,27 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
j1BornForceFactor = bornForceFactor[ jAtom ];
j1AtomicRadii = atomicRadii[ jAtom ];
j1ScaledAtomicRadii = scaledAtomicRadii[ jAtom ];
j1pos = posq[ jAtom ].xyz;
j1AtomicRadii = scaledAtomicRadii[ jAtom ].y;
j1ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j1pos = posq[ jAtom ];
jAtom.x += 1.0f;
j2BornForceFactor = bornForceFactor[ jAtom ];
j2AtomicRadii = atomicRadii[ jAtom ];
j2ScaledAtomicRadii = scaledAtomicRadii[ jAtom ];
j2pos = posq[ jAtom ].xyz;
j2AtomicRadii = scaledAtomicRadii[ jAtom ].y;
j2ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j2pos = posq[ jAtom ];
jAtom.x += 1.0f;
j3BornForceFactor = bornForceFactor[ jAtom ];
j3AtomicRadii = atomicRadii[ jAtom ];
j3ScaledAtomicRadii = scaledAtomicRadii[ jAtom ];
j3pos = posq[ jAtom ].xyz;
j3AtomicRadii = scaledAtomicRadii[ jAtom ].y;
j3ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j3pos = posq[ jAtom ];
jAtom.x += 1.0f;
j4BornForceFactor = bornForceFactor[ jAtom ];
j4AtomicRadii = atomicRadii[ jAtom ];
j4ScaledAtomicRadii = scaledAtomicRadii[ jAtom ];
j4pos = posq[ jAtom ].xyz;
j4AtomicRadii = scaledAtomicRadii[ jAtom ].y;
j4ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j4pos = posq[ jAtom ];
jAtom.x += 1.0f;
......@@ -557,14 +1329,12 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
// i = 1
//loop2Internal( d1, d2, d3, d4, jScaledAtomicRadii, i1BornForceFactor, i1AtomicRadii, de, bornSum ); //156 : 368
loop2Internal( d1, d2, d3, d4, jScaledAtomicRadii, i1BornForceFactor, i1AtomicRadii, de ); //156 : 368
bornForce1.xyz += de.x*d1; //3*2=6
bornForce1.xyz += de.y*d2; //3*2=6
bornForce1.xyz += de.z*d3; //3*2=6
bornForce1.xyz += de.w*d4; //3*2=6
//bornForce1.w += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// ---------------------------------------------------------------------------------------
......@@ -575,14 +1345,12 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
d3 = j3pos - i2pos; //3
d4 = j4pos - i2pos; //3
//loop2Internal( d1, d2, d3, d4, jScaledAtomicRadii, i2BornForceFactor, i2AtomicRadii, de, bornSum );//156 : 368
loop2Internal( d1, d2, d3, d4, jScaledAtomicRadii, i2BornForceFactor, i2AtomicRadii, de );//156 : 368
bornForce2.xyz += de.x*d1; //3*2=6
bornForce2.xyz += de.y*d2; //3*2=6
bornForce2.xyz += de.z*d3; //3*2=6
bornForce2.xyz += de.w*d4; //3*2=6
//bornForce2.w += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// ---------------------------------------------------------------------------------------
......@@ -599,7 +1367,6 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
bornForce3.xyz += de.y*d2; //3*2=6
bornForce3.xyz += de.z*d3; //3*2=6
bornForce3.xyz += de.w*d4; //3*2=6
//bornForce3.w += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// ---------------------------------------------------------------------------------------
......@@ -616,7 +1383,6 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
bornForce4.xyz += de.y*d2; //3*2=6
bornForce4.xyz += de.z*d3; //3*2=6
bornForce4.xyz += de.w*d4; //3*2=6
//bornForce4.w += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// ---------------------------------------------------------------------------------------
......@@ -627,7 +1393,6 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
d3 = j1pos - i3pos; //3
d4 = j1pos - i4pos; //3
// loop2Internal( d1, d2, d3, d4, iScaledAtomicRadii, j1BornForceFactor, j1AtomicRadii, de, bornSum );
loop2InternalNoSum( d1, d2, d3, d4, iScaledAtomicRadii, j1BornForceFactor, j1AtomicRadii, de ); //112 : 304
bornForce1.xyz += de.x*d1; //3*2=6
......@@ -644,7 +1409,6 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
d3 = j2pos - i3pos; //3
d4 = j2pos - i4pos; //3
// loop2Internal( d1, d2, d3, d4, iScaledAtomicRadii, j2BornForceFactor, j2AtomicRadii, de, bornSum );
loop2InternalNoSum( d1, d2, d3, d4, iScaledAtomicRadii, j2BornForceFactor, j2AtomicRadii, de ); //112 : 304
bornForce1.xyz += de.x*d1; //3*2=6
......@@ -661,7 +1425,6 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
d3 = j3pos - i3pos; //3
d4 = j3pos - i4pos; //3
// loop2Internal( d1, d2, d3, d4, iScaledAtomicRadii, j3BornForceFactor, j3AtomicRadii, de, bornSum );
loop2InternalNoSum( d1, d2, d3, d4, iScaledAtomicRadii, j3BornForceFactor, j3AtomicRadii, de ); //112 : 304
bornForce1.xyz += de.x*d1; //3*2=6
......@@ -678,7 +1441,6 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
d3 = j4pos - i3pos; //3
d4 = j4pos - i4pos; //3
// loop2Internal( d1, d2, d3, d4, iScaledAtomicRadii, j4BornForceFactor, j4AtomicRadii, de, bornSum );
loop2InternalNoSum( d1, d2, d3, d4, iScaledAtomicRadii, j4BornForceFactor, j4AtomicRadii, de ); //112 : 304
bornForce1.xyz += de.x*d1; //3*2=6
......@@ -697,27 +1459,37 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
}
}
/* ---------------------------------------------------------------------------------------
Calculate Born radii and Obc chain derivative term
numberOfAtoms: no. of atoms
streamWidth: atom stream width
forces input forces
atomicRadii atomic radii
bornRadii output Born radii
obcChain output Obc chain derivative
--------------------------------------------------------------------------------------- */
kernel void kObcBornRadii( float numberOfAtoms, float streamWidth, float4 forces<>,
float atomicRadii<>, out float bornRadii<>, out float obcChain<> ){
kernel void kPostObcLoop2_nobranch(
float repfac,
float atomStreamWidth,
float pStreamWidth,
float natoms,
float roundNatoms,
float iUnroll,
float conversion,
float mergeNonObcForces,
float4 inObcForces<>,
float3 nonObcForces<>,
float4 pstream1[][],
float4 pstream2[][],
float4 pstream3[][],
float4 pstream4[][],
float atomicRadii<>,
out float bornRadii<>,
out float obcChain<>,
out float3 outForces<> ){
// ---------------------------------------------------------------------------------------
float atomIndex, forceIndex, qIndex, qOff;
float2 pindex;
float i;
float sum2, sum3, bornSum, tanhSum, atomicRadiiOffset, obcIntermediate;
float4 o1,o2,o3,o4;
float4 tmp;
float2 iAtom;
float4 forces;
float expPlus, expMinus;
// ---------------------------------------------------------------------------------------
......@@ -728,8 +1500,54 @@ kernel void kObcBornRadii( float numberOfAtoms, float streamWidth, float4 forces
const float betaObc = 0.8f;
const float gammaObc = 4.85f;
const float dielectricOffset = 0.009f;
// ---------------------------------------------------------------------------------------
// given atom index find force indices and streams
pindex = indexof( outForces );
atomIndex = pindex.x + pindex.y*atomStreamWidth;
forceIndex = atomIndex;
// add current forces in inStream to forces stored in pstreams
// the .w entry is Born sum values; it will be used to calculate the
// Born radii and obcChain term
forces = inObcForces;
forces.w = 0.0f;
//forces = float4( 0.0f, 0.0f, 0.0f, 0.0f );
// sum over j-loop 'duplications' by gathering from pstreams
for( i = 0.0f; i < repfac; i += 1.0f ){
// qIndex = floor( forceIndex/iUnroll );
qIndex = round( (forceIndex - fmod( forceIndex, iUnroll))/iUnroll );
qOff = forceIndex - iUnroll*qIndex;
// pindex.y = floor( qIndex/ pStreamWidth );
pindex.y = round( (qIndex - fmod( qIndex, pStreamWidth ))/pStreamWidth );
// pindex.x = qIndex - pindex.y*pStreamWidth + qOff;
pindex.x = qIndex - pindex.y*pStreamWidth;
o1 = pstream1[ pindex ];
o2 = pstream2[ pindex ];
o3 = pstream3[ pindex ];
o4 = pstream4[ pindex ];
tmp = qOff < 0.5f ? o1 : o2;
tmp = qOff < 1.5f ? tmp : o3;
tmp = qOff < 2.5f ? tmp : o4;
forces += tmp;
forceIndex += roundNatoms;
}
// compute Born radii and ObcChain
atomicRadiiOffset = atomicRadii - dielectricOffset;
......@@ -748,223 +1566,21 @@ kernel void kObcBornRadii( float numberOfAtoms, float streamWidth, float4 forces
expPlus = exp( tanhSum );
expMinus = 1.0f/expPlus;
tanhSum = ( expPlus - expMinus )/( expPlus + expMinus );
bornRadii = 1.0f/( (1.0f/(atomicRadiiOffset)) - tanhSum/atomicRadii );
obcIntermediate = atomicRadiiOffset*( alphaObc - 2.0f*betaObc*bornSum + 3.0f*gammaObc*sum2 );
obcChain = (1.0f - tanhSum*tanhSum)*obcIntermediate/atomicRadii;
if( atomIndex >= natoms ){
bornRadii = 0.0f;
obcChain = 0.0f;
}
}
kernel float4 scalar_force_CDLJMerge( float4 qq, float epsfac, float4 sig, float4 eps, float4 r2 ) {
float4 invr, invrsig2, invrsig6;
float4 f;
invr = rsqrt( r2 );
invrsig2 = invr * sig;
invrsig2 = invrsig2 * invrsig2;
invrsig6 = invrsig2 * invrsig2 * invrsig2;
f = eps * ( 12.0f*invrsig6 - 6.0f ) * invrsig6;
f += epsfac*qq*invr;
f *= invr*invr*2.39f;
return f;
}
//???? FLOPS
kernel void kCalculateBornRadii( float numberOfAtoms, float roundedUpAtoms, float duplicationFactor,
float streamWidth, float fstreamWidth, float3 posq[][],
float2 scaledAtomicRadii[][],
out float4 bornForce<> ){
// ---------------------------------------------------------------------------------------
// atomic radii: radius - dielectric offset
float i1AtomicRadii, i2AtomicRadii, i3AtomicRadii, i4AtomicRadii;
// scaled atomic radii: scaleFactor*( radius - dielectric offset )
float j1ScaledAtomicRadii, j2ScaledAtomicRadii, j3ScaledAtomicRadii, j4ScaledAtomicRadii;
float4 jScaledAtomicRadii;
// (Born radii**2)*BornForce*ObcChainForce
// i,j coordinates
float3 i1pos, i2pos, i3pos, i4pos;
float3 j1pos, j2pos, j3pos, j4pos;
// delta coordinates
float3 d1, d2, d3, d4;
// indices
float2 iAtom;
float forceIndex;
//This is forceIndex mod numberOfAtoms, the true i index
float4 de, bornSum;
float iAtomLinearIndex, jLinind;
float2 jAtom;
float jEnd, jStart, jBlock;
float whichRep;
float tmp;
float tmp2;
// ---------------------------------------------------------------------------------------
// constants
const float I_Unroll = 4.0f;
// ---------------------------------------------------------------------------------------
// set linear index
iAtom = indexof( bornForce );
forceIndex = I_Unroll*( iAtom.x + iAtom.y*fstreamWidth );
iAtomLinearIndex = fmod( forceIndex, roundedUpAtoms );
// ---------------------------------------------------------------------------------------
// set gather coordinates
iAtom.x = fmod( iAtomLinearIndex, streamWidth );
iAtom.y = round( (iAtomLinearIndex - iAtom.x)/streamWidth );
// ---------------------------------------------------------------------------------------
// fetch position, (atomic radius - dielectric offset )
i1pos = posq[ iAtom ];
i1AtomicRadii = scaledAtomicRadii[ iAtom ].y;
iAtom.x += 1;
i2pos = posq[ iAtom ];
i2AtomicRadii = scaledAtomicRadii[ iAtom ].y;
iAtom.x += 1;
i3pos = posq[ iAtom ];
i3AtomicRadii = scaledAtomicRadii[ iAtom ].y;
iAtom.x += 1;
i4pos = posq[ iAtom ];
i4AtomicRadii = scaledAtomicRadii[ iAtom ].y;
bornForce = float4( 0.0f, 0.0f, 0.0f, 0.0f );
// inner loop indices setup
//changed the following instruction for rounding issues on some ASICs
//whichRep = floor( forceIndex / roundedUpAtoms );
tmp = fmod(forceIndex, roundedUpAtoms);
whichRep = round((forceIndex - tmp)/roundedUpAtoms);
jBlock = 1 + round( (numberOfAtoms-fmod( numberOfAtoms, (duplicationFactor*streamWidth) )) / (duplicationFactor*streamWidth) );
jStart = whichRep*jBlock;
jEnd = ( whichRep > duplicationFactor - 1.5f ) ? 999999.0f : (jStart + jBlock);
jAtom.y = jStart;
jLinind = jAtom.y*streamWidth;
// ---------------------------------------------------------------------------------------
while ( jAtom.y < jEnd && ( numberOfAtoms - jLinind ) > 0.9f ){
jAtom.x = 0.0f;
while ( jAtom.x < streamWidth && ( numberOfAtoms - jLinind ) > 0.9f ){
// ---------------------------------------------------------------------------------------
// gather required values
j1ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j1pos = posq[ jAtom ];
jAtom.x += 1.0f;
j2ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j2pos = posq[ jAtom ];
jAtom.x += 1.0f;
j3ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j3pos = posq[ jAtom ];
jAtom.x += 1.0f;
j4ScaledAtomicRadii = scaledAtomicRadii[ jAtom ].x;
j4pos = posq[ jAtom ];
jAtom.x += 1.0f;
jScaledAtomicRadii = float4( j1ScaledAtomicRadii, j2ScaledAtomicRadii, j3ScaledAtomicRadii, j4ScaledAtomicRadii );
// ---------------------------------------------------------------------------------------
// First set of 4 -- i == 1
d1 = j1pos - i1pos; //3
d2 = j2pos - i1pos; //3
d3 = j3pos - i1pos; //3
d4 = j4pos - i1pos; //3
// i = 1
bornSumInternal( d1, d2, d3, d4, jScaledAtomicRadii, i1AtomicRadii, bornSum ); //112 : 368(?)
bornForce.x += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// ---------------------------------------------------------------------------------------
// i = 2
d1 = j1pos - i2pos; //3
d2 = j2pos - i2pos; //3
d3 = j3pos - i2pos; //3
d4 = j4pos - i2pos; //3
bornSumInternal( d1, d2, d3, d4, jScaledAtomicRadii, i2AtomicRadii, bornSum );//112 : 368(?)
bornForce.y += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// ---------------------------------------------------------------------------------------
// i = 3
d1 = j1pos - i3pos; //3
d2 = j2pos - i3pos; //3
d3 = j3pos - i3pos; //3
d4 = j4pos - i3pos; //3
bornSumInternal( d1, d2, d3, d4, jScaledAtomicRadii, i3AtomicRadii, bornSum );//112 : 368(?)
bornForce.z += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// ---------------------------------------------------------------------------------------
// i = 4
d1 = j1pos - i4pos; //3
d2 = j2pos - i4pos; //3
d3 = j3pos - i4pos; //3
d4 = j4pos - i4pos; //3
bornSumInternal( d1, d2, d3, d4, jScaledAtomicRadii, i4AtomicRadii, bornSum );//112 : 368(?)
bornForce.w += bornSum.x + bornSum.y + bornSum.z + bornSum.w; //4
// add converted new forces to non-Obc forces
// ---------------------------------------------------------------------------------------
// increment linear index by 4
jLinind += 4.0f;
}
jAtom.y += 1.0f;
outForces = conversion*forces.xyz;
if( mergeNonObcForces > 0.1f ){
outForces += nonObcForces;
}
}
platforms/brook/src/gpu/kgbsa.h
View file @ 8c3878b2
......@@ -486,7 +486,6 @@ void kObcLoop2(
float streamWidth,
float fstreamWidth,
::brook::stream posq,
::brook::stream atomicRadii,
::brook::stream scaledAtomicRadii,
::brook::stream bornForceFactor,
::brook::stream bornForce1,
......
platforms/brook/src/gpu/kgbsa1.br deleted 100644 → 0
View file @ 116a99ba
/****************************************************************
* This file is part of the gpu acceleration library for gromacs.
* Author: Mark Friedrichs
*
* This kernel was developed in collaboration with
*
* Copyright (C) Pande Group, Stanford, 2006
*****************************************************************/
//96 : 228 Flops
kernel void loop1Internal( float3 d1, float3 d2, float3 d3, float3 d4, float4 jBornR,
float4 jQ, float iBornR, float iQ, out float4 dGpol_dr<>,
out float4 dGpol_dalpha2_ij<> ){
// ---------------------------------------------------------------------------------------
float4 r2, alpha2_ij, D_ij, expTerm, denominator2, denominator, Gpol;
// ---------------------------------------------------------------------------------------
r2 = float4( dot(d1, d1), dot( d2, d2 ), dot( d3, d3 ), dot( d4, d4 ) ); //20
alpha2_ij = jBornR*iBornR; //4
D_ij = r2/(4.0f*alpha2_ij); //8
expTerm = exp( -D_ij ); //4 : 80
denominator2 = r2 + alpha2_ij*expTerm; //8
denominator = sqrt( denominator2 ); //4 : 60
Gpol = jQ/denominator; //4
Gpol *= iQ;//4
dGpol_dr = -Gpol*( 1.0f - 0.25f*expTerm )/denominator2; //16
dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*( 1.0f + D_ij )*jBornR/denominator2; //24
}
/* ---------------------------------------------------------------------------------------
Calculate nonpolar ACE term (Simbios)
bornRadius: Born radius
vdwRadius: Vdw radius
duplicationFactor: duplication factor
aceForce: ACE term
--------------------------------------------------------------------------------------- */
kernel void kAceNonPolarLoop1( float iBornRadius, float iVdwRadius, float duplicationFactor,
out float aceForce<> ){
// ---------------------------------------------------------------------------------------
// nonpolar term
float iSurface;
float iAceTerm;
// ---------------------------------------------------------------------------------------
// constants
// solvent radius
const float probeRadius = 0.14f;
// PI*24*0.0f049 (0.0f054=asolv from Tinker)
//const float PI_24_aI = -0.3694512961;
const float PI_24_aI = -400.71504079f;
// ---------------------------------------------------------------------------------------
// etch i position and partial charge
// e = ai * term * (ri+probe)**2 * (ri/rb)**6
// (drbi) = drb(i) - 6.0fd0*e/rb
// (rI+probe)**2
iSurface = (iVdwRadius+probeRadius);
iSurface = iSurface*iSurface;
// (rI/rB)**6
iAceTerm = iVdwRadius/iBornRadius;
iAceTerm = iAceTerm*iAceTerm*iAceTerm;
iAceTerm = iAceTerm*iAceTerm;
aceForce = iSurface*iAceTerm*PI_24_aI/(duplicationFactor*iBornRadius);
}
/* ---------------------------------------------------------------------------------------
Calculate first loop force terms (Simbios)
numberOfAtoms: no. of atoms
roundedUpAtoms: rounded up number of atoms -- accounts for unrolling
duplicationFactor: number of threads for inner loop
streamWidth: atom stream width
fstreamWidth: force stream width (output -- i-unroll)
soluteDielectric: solute dielectric
solventDielectric: solvent dielectric
includeAce: include ACE term
posq: atom positions and charge
bornRadii: Born radii
nonpolarForce: nonpolar force (0 if nonpolar not included, else
ACE value)
bornForce1: i-unroll first force component, including dBornR/dr in .w
bornForce2: i-unroll second force component, including dBornR/dr in .w
bornForce3: i-unroll first force component, including dBornR/dr in .w
bornForce4: i-unroll second force component, including dBornR/dr in .w
--------------------------------------------------------------------------------------- */
//FLOPS = 544
kernel void kObcLoop1( float numberOfAtoms, float roundedUpAtoms, float duplicationFactor,
float streamWidth, float fstreamWidth, float soluteDielectric,
float solventDielectric, float includeAce,
float4 posq[][], float bornRadii[][], float atomicRadii[][],
out float4 bornForce1<>, out float4 bornForce2<>,
out float4 bornForce3<>, out float4 bornForce4<> ){
// ---------------------------------------------------------------------------------------
// Born radii
float i1BornR, i2BornR, i3BornR, i4BornR;
float j1BornR, j2BornR, j3BornR, j4BornR;
float4 jBornR;
// atomic radii
float i1AtomicR, i2AtomicR, i3AtomicR, i4AtomicR;
// i,j coordinates
float3 i1Pos, i2Pos, i3Pos, i4Pos;
float3 j1Pos, j2Pos, j3Pos, j4Pos;
float4 j1PosQ, j2PosQ, j3PosQ, j4PosQ;
// i, j partial charges
float i1Q, i2Q, i3Q, i4Q;
float j1Q, j2Q, j3Q, j4Q;
float4 jQ;
float aceForce;
// delta coordinates
float3 d1, d2, d3, d4;
// intermediate terms
float4 dGpol_dr, dGpol_dalpha2_ij;
// indices
float2 iAtom;
float forceIndex;
// This is forceIndex mod numberOfAtoms, the true i index
float iAtomLinearIndex, jLinind;
float2 jAtom;
float jEnd, jStart, jBlock;
float whichRep;
float tmp;
// ---------------------------------------------------------------------------------------
// electricConstant = -166.0f2691;
// preFactor = 2.0f*electricConstant*(1.0f - (1.0f/waterDielectric))
float preFactor = -332.05382f;
const float I_Unroll = 4.0f;
const float3 zero3 = float3( 0.0f, 0.0f, 0.0f );
// ---------------------------------------------------------------------------------------
preFactor *= ( (1.0f/soluteDielectric) - (1.0f/solventDielectric) );
iAtom = indexof( bornForce1 );
forceIndex = I_Unroll*( iAtom.x + iAtom.y*fstreamWidth );
iAtomLinearIndex = fmod( forceIndex, roundedUpAtoms );
// ---------------------------------------------------------------------------------------
// set gather index
iAtom.x = fmod( iAtomLinearIndex, streamWidth );
iAtom.y = round( (iAtomLinearIndex - fmod(iAtomLinearIndex, streamWidth ))/streamWidth );
// ---------------------------------------------------------------------------------------
// etch i1 position and partial charge
jQ = posq[ iAtom ];
i1BornR = bornRadii[ iAtom ];
i1AtomicR = atomicRadii[ iAtom ];
i1Pos = jQ.xyz;
i1Q = jQ.w;
i1Q *= preFactor;
kAceNonPolarLoop1( i1BornR, i1AtomicR, duplicationFactor, aceForce );
bornForce1.xyz = zero3;
bornForce1.w = includeAce > 0.5f ? aceForce : 0.0f;
// ---------------------------------------------------------------------------------------
// etch i2 position and partial charge
iAtom.x += 1;
jQ = posq[ iAtom ];
i2BornR = bornRadii[ iAtom ];
i2AtomicR = atomicRadii[ iAtom ];
i2Pos = jQ.xyz;
i2Q = jQ.w;
i2Q *= preFactor;
kAceNonPolarLoop1( i2BornR, i2AtomicR, duplicationFactor, aceForce );
bornForce2.xyz = zero3;
bornForce2.w = includeAce > 0.5f ? aceForce : 0.0f;
// ---------------------------------------------------------------------------------------
// etch i3 position and partial charge
iAtom.x += 1;
jQ = posq[ iAtom ];
i3BornR = bornRadii[ iAtom ];
i3AtomicR = atomicRadii[ iAtom ];
i3Pos = jQ.xyz;
i3Q = jQ.w;
i3Q *= preFactor;
kAceNonPolarLoop1( i3BornR, i3AtomicR, duplicationFactor, aceForce );
bornForce3.xyz = zero3;
bornForce3.w = includeAce > 0.5f ? aceForce : 0.0f;
// ---------------------------------------------------------------------------------------
// etch i4 position and partial charge
iAtom.x += 1;
jQ = posq[ iAtom ];
i4BornR = bornRadii[ iAtom ];
i4AtomicR = atomicRadii[ iAtom ];
i4Pos = jQ.xyz;
i4Q = jQ.w;
i4Q *= preFactor;
kAceNonPolarLoop1( i4BornR, i4AtomicR, duplicationFactor, aceForce );
bornForce4.xyz = zero3;
bornForce4.w = includeAce > 0.5f ? aceForce : 0.0f;
// ---------------------------------------------------------------------------------------
// inner loop setup
// if dupFac == 4, I_UnRoll =2, then breaking inner loop into two segments
// to increase number of threads in flight
// forceStreamSz = N*RepFac/I_UnRoll
// forceIndex = I_UnRoll*( a.x + a.y*forceStreamSz )
// whichRep = 0 or 1
// jBlock = 1 + floor[ N/(duplicationFactor*streamWidth) ]
//changed the following instruction for rounding issues on some ASICs
//whichRep = floor( forceIndex / roundedUpAtoms );
tmp = fmod(forceIndex, roundedUpAtoms);
whichRep = round((forceIndex - tmp)/roundedUpAtoms);
//jBlock = 1 + floor( numberOfAtoms/(duplicationFactor*streamWidth ) );
tmp = fmod( numberOfAtoms, (duplicationFactor*streamWidth ) );
jBlock = 1 + round( (numberOfAtoms-tmp)/(duplicationFactor*streamWidth ));
jStart = whichRep*jBlock;
jEnd = ( whichRep > duplicationFactor - 1.5f ) ? 999999.0f : (jStart + jBlock);
jAtom.y = jStart;
jLinind = jAtom.y*streamWidth;
// ---------------------------------------------------------------------------------------
while ( jAtom.y < jEnd && ( numberOfAtoms - jLinind ) > 0.9f ){
jAtom.x = 0.0f;
while ( jAtom.x < streamWidth && ( numberOfAtoms - jLinind ) > 0.9f ) {
// ---------------------------------------------------------------------------------------
// gather required values
j1PosQ = posq[ jAtom ];
j1Pos = j1PosQ.xyz;
j1Q = j1PosQ.w;
j1BornR = bornRadii[ jAtom ];
jAtom.x += 1.0f;
j2PosQ = posq[ jAtom ];
j2Pos = j2PosQ.xyz;
j2Q = j2PosQ.w;
j2BornR = bornRadii[ jAtom ];
jAtom.x += 1.0f;
j3PosQ = posq[ jAtom ];
j3Pos = j3PosQ.xyz;
j3Q = j3PosQ.w;
j3BornR = bornRadii[ jAtom ];
jAtom.x += 1.0f;
j4PosQ = posq[ jAtom ];
j4Pos = j4PosQ.xyz;
j4Q = j4PosQ.w;
j4BornR = bornRadii[ jAtom ];
jAtom.x += 1.0f;
jBornR = float4( j1BornR, j2BornR, j3BornR, j4BornR );
jQ = float4( j1Q, j2Q, j3Q, j4Q );
// ---------------------------------------------------------------------------------------
// i == 1
d1 = i1Pos - j1Pos; //3
d2 = i1Pos - j2Pos; //3
d3 = i1Pos - j3Pos; //3
d4 = i1Pos - j4Pos; //3
loop1Internal( d1, d2, d3, d4, jBornR, jQ, i1BornR, i1Q, dGpol_dr, dGpol_dalpha2_ij ); //96 : 228
bornForce1.xyz += dGpol_dr.x*d1; //6
bornForce1.xyz += dGpol_dr.y*d2; //6
bornForce1.xyz += dGpol_dr.z*d3; //6
bornForce1.xyz += dGpol_dr.w*d4; //6
bornForce1.w += dGpol_dalpha2_ij.x + dGpol_dalpha2_ij.y + dGpol_dalpha2_ij.z + dGpol_dalpha2_ij.w; //4
// ---------------------------------------------------------------------------------------
// i == 2
d1 = i2Pos - j1Pos; //3
d2 = i2Pos - j2Pos; //3
d3 = i2Pos - j3Pos; //3
d4 = i2Pos - j4Pos; //3
loop1Internal( d1, d2, d3, d4, jBornR, jQ, i2BornR, i2Q, dGpol_dr, dGpol_dalpha2_ij ); //96 : 228
bornForce2.xyz += dGpol_dr.x*d1; //6
bornForce2.xyz += dGpol_dr.y*d2; //6
bornForce2.xyz += dGpol_dr.z*d3; //6
bornForce2.xyz += dGpol_dr.w*d4; //6
bornForce2.w += dGpol_dalpha2_ij.x + dGpol_dalpha2_ij.y + dGpol_dalpha2_ij.z + dGpol_dalpha2_ij.w; //4
// ---------------------------------------------------------------------------------------
// i == 3
d1 = i3Pos - j1Pos; //3
d2 = i3Pos - j2Pos; //3
d3 = i3Pos - j3Pos; //3
d4 = i3Pos - j4Pos; //3
loop1Internal( d1, d2, d3, d4, jBornR, jQ, i3BornR, i3Q, dGpol_dr, dGpol_dalpha2_ij ); //96 : 228
bornForce3.xyz += dGpol_dr.x*d1; //6
bornForce3.xyz += dGpol_dr.y*d2; //6
bornForce3.xyz += dGpol_dr.z*d3; //6
bornForce3.xyz += dGpol_dr.w*d4; //6
bornForce3.w += dGpol_dalpha2_ij.x + dGpol_dalpha2_ij.y + dGpol_dalpha2_ij.z + dGpol_dalpha2_ij.w; //4
// ---------------------------------------------------------------------------------------
// i == 4
d1 = i4Pos - j1Pos; //3
d2 = i4Pos - j2Pos; //3
d3 = i4Pos - j3Pos; //3
d4 = i4Pos - j4Pos; //3
loop1Internal( d1, d2, d3, d4, jBornR, jQ, i4BornR, i4Q, dGpol_dr, dGpol_dalpha2_ij ); //96 : 228
bornForce4.xyz += dGpol_dr.x*d1; //6
bornForce4.xyz += dGpol_dr.y*d2; //6
bornForce4.xyz += dGpol_dr.z*d3; //6
bornForce4.xyz += dGpol_dr.w*d4; //6
bornForce4.w += dGpol_dalpha2_ij.x + dGpol_dalpha2_ij.y + dGpol_dalpha2_ij.z + dGpol_dalpha2_ij.w; //4
// ---------------------------------------------------------------------------------------
jLinind += 4.0f;
}
jAtom.y += 1.0f;
}
}
platforms/brook/src/gpu/kmerge.br
View file @ 8c3878b2
......@@ -225,154 +225,6 @@ kernel void kMergeFloat4_4(
}
}
kernel void kPostObcLoop1(
float repfac,
float atomStreamWidth,
float pStreamWidth,
float natoms,
float roundNatoms,
float iUnroll,
float4 pstream1[][],
float4 pstream2[][],
float4 pstream3[][],
float4 pstream4[][],
float obcChain<>,
float bornRadii<>,
out float4 outstream<>,
out float bornRadii2Force<> )
{
// ---------------------------------------------------------------------------------------
float atomIndex, forceIndex, qIndex, qOff;
float2 pindex;
float i;
float4 o1,o2,o3,o4;
float4 tmp;
float4 zero4;
// given atom index find force indices and streams
pindex = indexof( outstream );
atomIndex = pindex.x + pindex.y*atomStreamWidth;
forceIndex = atomIndex;
zero4 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
outstream = zero4;
for( i = 0.0f; i < repfac; i += 1.0f ){
// qIndex = floor( forceIndex/iUnroll );
qIndex = round( (forceIndex - fmod( forceIndex, iUnroll))/iUnroll );
qOff = forceIndex - iUnroll*qIndex;
// pindex.y = floor( qIndex/ pStreamWidth );
pindex.y = round( (qIndex - fmod( qIndex, pStreamWidth ))/pStreamWidth );
// pindex.x = qIndex - pindex.y*pStreamWidth + qOff;
pindex.x = qIndex - pindex.y*pStreamWidth;
// outstream += float4( forceIndex, qIndex, pindex.x, pindex.y );
// This is going to cause really divergent code and we are
// going to end up doing all the fetches anyway...
if ( qOff < 0.5f ){
outstream += pstream1[ pindex ];
} else if( qOff < 1.5f ){
outstream += pstream2[ pindex ];
} else if( qOff < 2.5f ){
outstream += pstream3[ pindex ];
} else {
outstream += pstream4[ pindex ];
}
// o1 = pstream1[ pindex ];
// o2 = pstream2[ pindex ];
// o3 = pstream3[ pindex ];
// o4 = pstream4[ pindex ];
// tmp = qOff < 0.5f ? o1 : o2;
// tmp = qOff < 1.5f ? tmp : o3;
// tmp = qOff < 2.5f ? tmp : o4;
// outstream += tmp;
forceIndex += roundNatoms;
}
bornRadii2Force = obcChain*bornRadii*bornRadii*outstream.w;
}
// The inner loop by definition creates divergent paths. Chances are
// fair that we will take all sides of the branch anyway, so this
// verion uses manual predication
kernel void kPostObcLoop1_nobranch(
float repfac,
float atomStreamWidth,
float pStreamWidth,
float natoms,
float roundNatoms,
float iUnroll,
float4 pstream1[][],
float4 pstream2[][],
float4 pstream3[][],
float4 pstream4[][],
float obcChain<>,
float bornRadii<>,
out float4 outstream<>,
out float bornRadii2Force<> )
{
// ---------------------------------------------------------------------------------------
float atomIndex, forceIndex, qIndex, qOff;
float2 pindex;
float i;
float4 o1,o2,o3,o4;
float4 tmp;
// given atom index find force indices and streams
pindex = indexof( outstream );
atomIndex = pindex.x + pindex.y*atomStreamWidth;
forceIndex = atomIndex;
outstream = float4( 0.0f, 0.0f, 0.0f, 0.0f );
for( i = 0.0f; i < repfac; i += 1.0f ){
// qIndex = floor( forceIndex/iUnroll );
qIndex = round( (forceIndex - fmod( forceIndex, iUnroll))/iUnroll );
qOff = forceIndex - iUnroll*qIndex;
// pindex.y = floor( qIndex/ pStreamWidth );
pindex.y = round( (qIndex - fmod( qIndex, pStreamWidth ))/pStreamWidth );
// pindex.x = qIndex - pindex.y*pStreamWidth + qOff;
pindex.x = qIndex - pindex.y*pStreamWidth;
o1 = pstream1[ pindex ];
o2 = pstream2[ pindex ];
o3 = pstream3[ pindex ];
o4 = pstream4[ pindex ];
tmp = qOff < 0.5f ? o1 : o2;
tmp = qOff < 1.5f ? tmp : o3;
tmp = qOff < 2.5f ? tmp : o4;
outstream += tmp;
forceIndex += roundNatoms;
}
bornRadii2Force = obcChain*bornRadii*bornRadii*outstream.w;
}
kernel void kSetValue4( float value, out float4 outstream<> ){
outstream = float4( value, value, value, value );
}
......@@ -522,354 +374,6 @@ kernel void kAddAndMergeFloat4_4(
}
}
kernel void kPostObcLoop2(
float repfac,
float atomStreamWidth,
float pStreamWidth,
float natoms,
float roundNatoms,
float iUnroll,
float conversion,
float mergeNonObcForces,
float4 inObcForces<>,
float3 nonObcForces<>,
float4 pstream1[][],
float4 pstream2[][],
float4 pstream3[][],
float4 pstream4[][],
float atomicRadii<>,
out float bornRadii<>,
out float obcChain<>,
out float3 outForces<> ){
// ---------------------------------------------------------------------------------------
float atomIndex, forceIndex, qIndex, qOff;
float2 pindex;
float i;
float sum2, sum3, bornSum, tanhSum, atomicRadiiOffset, obcIntermediate;
float2 iAtom;
float4 forces;
float expPlus, expMinus;
// ---------------------------------------------------------------------------------------
// constants -- OBC Type II
const float alphaObc = 1.0f;
const float betaObc = 0.8f;
const float gammaObc = 4.85f;
const float dielectricOffset = 0.009f;
// ---------------------------------------------------------------------------------------
// given atom index find force indices and streams
pindex = indexof( outForces );
atomIndex = pindex.x + pindex.y*atomStreamWidth;
forceIndex = atomIndex;
// add current forces in inStream to forces stored in pstreams
// the .w entry is Born sum values; it will be used to calculate the
// Born radii and obcChain term
forces = inObcForces;
forces.w = 0.0f;
//forces = float4( 0.0f, 0.0f, 0.0f, 0.0f );
// sum over j-loop 'duplications' by gathering from pstreams
for( i = 0.0f; i < repfac; i += 1.0f ){
// qIndex = floor( forceIndex/iUnroll );
qIndex = round( (forceIndex - fmod( forceIndex, iUnroll))/iUnroll );
qOff = forceIndex - iUnroll*qIndex;
// pindex.y = floor( qIndex/ pStreamWidth );
pindex.y = round( (qIndex - fmod( qIndex, pStreamWidth ))/pStreamWidth );
// pindex.x = qIndex - pindex.y*pStreamWidth + qOff;
pindex.x = qIndex - pindex.y*pStreamWidth;
if( qOff < 0.5f ){
forces += pstream1[ pindex ];
} else if( qOff < 1.5f ){
forces += pstream2[ pindex ];
} else if( qOff < 2.5f ){
forces += pstream3[ pindex ];
} else {
forces += pstream4[ pindex ];
}
forceIndex += roundNatoms;
}
// compute Born radii and ObcChain
atomicRadiiOffset = atomicRadii - dielectricOffset;
bornSum = forces.w;
bornSum *= 0.5f*atomicRadiiOffset;
sum2 = bornSum*bornSum;
sum3 = bornSum*sum2;
// Tanh does not exist?
// calculate [ exp(x) - exp(-x) ]/[ exp(x) + exp(-x) ]
// tanhSum = tanh( bornSum - betaObc*sum2 + gammaObc*sum3 );
tanhSum = bornSum - betaObc*sum2 + gammaObc*sum3;
expPlus = exp( tanhSum );
expMinus = 1.0f/expPlus;
tanhSum = ( expPlus - expMinus )/( expPlus + expMinus );
bornRadii = 1.0f/( (1.0f/(atomicRadiiOffset)) - tanhSum/atomicRadii );
obcIntermediate = atomicRadiiOffset*( alphaObc - 2.0f*betaObc*bornSum + 3.0f*gammaObc*sum2 );
obcChain = (1.0f - tanhSum*tanhSum)*obcIntermediate/atomicRadii;
if( atomIndex >= natoms ){
bornRadii = 0.0f;
obcChain = 0.0f;
}
// add converted new forces to non-Obc forces
outForces = conversion*forces.xyz;
if( mergeNonObcForces > 0.1f ){
outForces += nonObcForces;
}
}
kernel void kPostObcLoop2_nobranch(
float repfac,
float atomStreamWidth,
float pStreamWidth,
float natoms,
float roundNatoms,
float iUnroll,
float conversion,
float mergeNonObcForces,
float4 inObcForces<>,
float3 nonObcForces<>,
float4 pstream1[][],
float4 pstream2[][],
float4 pstream3[][],
float4 pstream4[][],
float atomicRadii<>,
out float bornRadii<>,
out float obcChain<>,
out float3 outForces<> ){
// ---------------------------------------------------------------------------------------
float atomIndex, forceIndex, qIndex, qOff;
float2 pindex;
float i;
float sum2, sum3, bornSum, tanhSum, atomicRadiiOffset, obcIntermediate;
float4 o1,o2,o3,o4;
float4 tmp;
float2 iAtom;
float4 forces;
float expPlus, expMinus;
// ---------------------------------------------------------------------------------------
// constants -- OBC Type II
const float alphaObc = 1.0f;
const float betaObc = 0.8f;
const float gammaObc = 4.85f;
const float dielectricOffset = 0.009f;
// ---------------------------------------------------------------------------------------
// given atom index find force indices and streams
pindex = indexof( outForces );
atomIndex = pindex.x + pindex.y*atomStreamWidth;
forceIndex = atomIndex;
// add current forces in inStream to forces stored in pstreams
// the .w entry is Born sum values; it will be used to calculate the
// Born radii and obcChain term
forces = inObcForces;
forces.w = 0.0f;
//forces = float4( 0.0f, 0.0f, 0.0f, 0.0f );
// sum over j-loop 'duplications' by gathering from pstreams
for( i = 0.0f; i < repfac; i += 1.0f ){
// qIndex = floor( forceIndex/iUnroll );
qIndex = round( (forceIndex - fmod( forceIndex, iUnroll))/iUnroll );
qOff = forceIndex - iUnroll*qIndex;
// pindex.y = floor( qIndex/ pStreamWidth );
pindex.y = round( (qIndex - fmod( qIndex, pStreamWidth ))/pStreamWidth );
// pindex.x = qIndex - pindex.y*pStreamWidth + qOff;
pindex.x = qIndex - pindex.y*pStreamWidth;
o1 = pstream1[ pindex ];
o2 = pstream2[ pindex ];
o3 = pstream3[ pindex ];
o4 = pstream4[ pindex ];
tmp = qOff < 0.5f ? o1 : o2;
tmp = qOff < 1.5f ? tmp : o3;
tmp = qOff < 2.5f ? tmp : o4;
forces += tmp;
forceIndex += roundNatoms;
}
// compute Born radii and ObcChain
atomicRadiiOffset = atomicRadii - dielectricOffset;
bornSum = forces.w;
bornSum *= 0.5f*atomicRadiiOffset;
sum2 = bornSum*bornSum;
sum3 = bornSum*sum2;
// Tanh does not exist?
// calculate [ exp(x) - exp(-x) ]/[ exp(x) + exp(-x) ]
// tanhSum = tanh( bornSum - betaObc*sum2 + gammaObc*sum3 );
tanhSum = bornSum - betaObc*sum2 + gammaObc*sum3;
expPlus = exp( tanhSum );
expMinus = 1.0f/expPlus;
tanhSum = ( expPlus - expMinus )/( expPlus + expMinus );
bornRadii = 1.0f/( (1.0f/(atomicRadiiOffset)) - tanhSum/atomicRadii );
obcIntermediate = atomicRadiiOffset*( alphaObc - 2.0f*betaObc*bornSum + 3.0f*gammaObc*sum2 );
obcChain = (1.0f - tanhSum*tanhSum)*obcIntermediate/atomicRadii;
if( atomIndex >= natoms ){
bornRadii = 0.0f;
obcChain = 0.0f;
}
// add converted new forces to non-Obc forces
outForces = conversion*forces.xyz;
if( mergeNonObcForces > 0.1f ){
outForces += nonObcForces;
}
}
kernel void kPostCalculateBornRadii_nobranch(
float repfac,
float atomStreamWidth,
float pStreamWidth,
float natoms,
float roundNatoms,
float iUnroll,
float conversion,
float mergeNonObcForces,
float4 pstream1[][],
float atomicRadii<>,
out float bornRadii<>,
out float obcChain<> ){
// ---------------------------------------------------------------------------------------
float atomIndex, forceIndex, qIndex, qOff;
float2 pindex;
float i;
float sum2, sum3, bornSum, tanhSum, atomicRadiiOffset, obcIntermediate;
float4 o1;
float tmp;
float2 iAtom;
float4 forces;
float expPlus, expMinus;
// ---------------------------------------------------------------------------------------
// constants -- OBC Type II
const float alphaObc = 1.0f;
const float betaObc = 0.8f;
const float gammaObc = 4.85f;
const float dielectricOffset = 0.009f;
// ---------------------------------------------------------------------------------------
// given atom index find force indices and streams
pindex = indexof( obcChain );
atomIndex = pindex.x + pindex.y*atomStreamWidth;
forceIndex = atomIndex;
// add current forces in inStream to forces stored in pstreams
// the .w entry is Born sum values; it will be used to calculate the
// Born radii and obcChain term
bornSum = 0.0f;
// sum over j-loop 'duplications' by gathering from pstreams
for( i = 0.0f; i < repfac; i += 1.0f ){
qIndex = round( (forceIndex - fmod( forceIndex, iUnroll))/iUnroll );
qOff = forceIndex - iUnroll*qIndex;
pindex.y = round( (qIndex - fmod( qIndex, pStreamWidth ))/pStreamWidth );
pindex.x = qIndex - pindex.y*pStreamWidth;
o1 = pstream1[ pindex ];
tmp = qOff < 0.5f ? o1.x : o1.y;
tmp = qOff < 1.5f ? tmp : o1.z;
tmp = qOff < 2.5f ? tmp : o1.w;
bornSum += tmp;
forceIndex += roundNatoms;
}
// compute Born radii and ObcChain
atomicRadiiOffset = atomicRadii - dielectricOffset;
bornSum *= 0.5f*atomicRadiiOffset;
sum2 = bornSum*bornSum;
sum3 = bornSum*sum2;
// Tanh does not exist?
// calculate [ exp(x) - exp(-x) ]/[ exp(x) + exp(-x) ]
// tanhSum = tanh( bornSum - betaObc*sum2 + gammaObc*sum3 );
tanhSum = bornSum - betaObc*sum2 + gammaObc*sum3;
expPlus = exp( tanhSum );
expMinus = 1.0f/expPlus;
tanhSum = ( expPlus - expMinus )/( expPlus + expMinus );
bornRadii = 1.0f/( (1.0f/(atomicRadiiOffset)) - tanhSum/atomicRadii );
obcIntermediate = atomicRadiiOffset*( alphaObc - 2.0f*betaObc*bornSum + 3.0f*gammaObc*sum2 );
obcChain = (1.0f - tanhSum*tanhSum)*obcIntermediate/atomicRadii;
if( atomIndex >= natoms ){
bornRadii = 0.0f;
obcChain = 0.0f;
}
}
/* Add forces from two streams */
kernel void kAddForces3_4( float conversion, float3 force1<>, float4 force2<>, out float3 outForce<> ){
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment