/home/julian/svn/gpucpumd/mdgpu2/gpu_forces.cu

Go to the documentation of this file.

#include "gpu_forces.cuh"




__device__ void gpuAddLJForce(FLOAT_ARRAY_TYPE *dForce, FLOAT_ARRAY_TYPE *dA, FLOAT_ARRAY_TYPE *dB, float *boxsize)
{
    float invR2, invR6, fd;
    FLOAT_ARRAY_TYPE d;
    gpuDistPBC( &d, dA, dB, boxsize);
    invR2 = 1.0f / ( d.x * d.x + d.y * d.y + d.z * d.z );
    invR6 = powf( invR2, 3.0f );
    fd = 48.0f * invR2 * invR6 * ( invR6 - 0.5f );
    (*dForce).x += fd * d.x;
    (*dForce).y += fd * d.y;
    (*dForce).z += fd * d.z;
}

__device__ void gpuAddLJForceCut(FLOAT_ARRAY_TYPE *dForce, FLOAT_ARRAY_TYPE *dA, FLOAT_ARRAY_TYPE *dB, float *boxsize, float *rCutSq)
{
    float invR2, invR6, fd;
    FLOAT_ARRAY_TYPE d;
    gpuDistPBC( &d, dA, dB, boxsize);
    invR2 = d.x * d.x + d.y * d.y + d.z * d.z;
    if ( (*rCutSq) > invR2 )
    {
        invR2 = 1.0f / invR2;
        invR6 = powf( invR2, 3.0f );
        fd = 48.0f * invR2 * invR6 * ( invR6 - 0.5f );
        (*dForce).x += fd * d.x;
        (*dForce).y += fd * d.y;
        (*dForce).z += fd * d.z;
    }
}


__device__ void gpuAddHarmonicBondForce(FLOAT_ARRAY_TYPE *dForce, FLOAT_ARRAY_TYPE *dA, FLOAT_ARRAY_TYPE *dB, float *dK, float *boxsize)
{
    FLOAT_ARRAY_TYPE d;
    gpuDistPBC(&d, dA, dB, boxsize);
    (*dForce).x -= (*dK) * d.x;
    (*dForce).y -= (*dK) * d.y;
    (*dForce).z -= (*dK) * d.z;
}


__global__ void gpuLJForces(FLOAT_ARRAY_TYPE *dPos, FLOAT_ARRAY_TYPE *dForce, int N, float boxsize)
{
    int tid = threadIdx.x;
    int idx = blockDim.x * blockIdx.x + tid;
    if (idx < N)
    {
        FLOAT_ARRAY_TYPE tPos = dPos[idx];
        FLOAT_ARRAY_TYPE f;
        f.x = 0.0f;
        f.y = 0.0f;
        f.z = 0.0f;
        for (int i = 0; i < N; ++i)
        {
            if ( i != idx) //exclude self interaction
            {
                gpuAddLJForce( &f, &tPos, &dPos[i], &boxsize);
            }
        }
        gpuCumulateVec(&dForce[idx], &f);
    }
}


__global__ void gpuLJForcesSmem(FLOAT_ARRAY_TYPE *dPos, FLOAT_ARRAY_TYPE *dForce, int N, float boxsize)
{
    //Reduced units are used -> epsilon=sigma=1
    int tid = threadIdx.x;
    int offset = blockDim.x * blockIdx.x;
    int rowOffset = offset + blockDim.x;
    int idx = offset + tid;
    //load particle positions into shared memory
    extern __shared__ FLOAT_ARRAY_TYPE sPos[];
    sPos[tid] = dPos[idx];

    //sync such that everything is in shared memory before going on
    __syncthreads();

    //calculate forces for particles which are already in shared memory
    FLOAT_ARRAY_TYPE f;

    f.x = 0.0f;
    f.y = 0.0f;
    f.z = 0.0f;
    //block before thread block in device memory
    for (int i = 0; i < offset; ++i)
    {
        gpuAddLJForce( &f, &sPos[tid], &dPos[i], &boxsize);
    }
    //block in shared memory
    for (int i = 0; i < blockDim.x; ++i)
    {
        if ( i != tid) //excude self interaction
        {
            gpuAddLJForce( &f, &sPos[tid], &sPos[i], &boxsize);
        }
    }
    //block after threadblock in device memory
    for (int i = rowOffset; i < N; ++i)
    {
        gpuAddLJForce( &f, &sPos[tid], &dPos[i], &boxsize);
    }
    gpuCumulateVec(&dForce[idx], &f);
}


__global__ void gpuLJForcesSmemSmem(FLOAT_ARRAY_TYPE *dPos, FLOAT_ARRAY_TYPE *dForce, int N, float boxsize)
{
    //Reduced units are used -> epsilon=sigma=1
    int tid = threadIdx.x;
    int bid = blockIdx.x;
    int bOffset;//current block offset
    int bdim = blockDim.x;
    int offset = blockDim.x * blockIdx.x;
    int idx = offset + tid;
    //load particle positions into shared memory
    extern __shared__ FLOAT_ARRAY_TYPE sPos[];
    sPos[tid] = dPos[idx];

    //sync such that everything is in shared memory before going on
    __syncthreads();

    //calculate forces for particles which are already in shared memory
    FLOAT_ARRAY_TYPE f;
    f.x = 0.0f;
    f.y = 0.0f;
    f.z = 0.0f;

    for (int j = 0; j < bid; ++j)
    {
        bOffset = j * bdim;
        //load block of positions into shared memory
        sPos[blockDim.x+tid] = dPos[bOffset+tid];
        __syncthreads();
        //do force calculation for this block
        for (int i = 0; i < bdim; ++i)
        {
            gpuAddLJForce( &f, &sPos[tid], &sPos[bdim+i], &boxsize);
        }
        __syncthreads();
    }

    //do force calculation inside thread block
    for (int i = 0; i < blockDim.x; ++i)
    {
        if ( i != tid )
        {
            gpuAddLJForce( &f, &sPos[tid], &sPos[i], &boxsize);
        }
    }
    __syncthreads();

    //do force calculation for blocks after thread block
    for (int j = bid + 1; j < gridDim.x; ++j)
    {
        bOffset = j * bdim;
        //load block of positions into shared memory
        sPos[bdim+tid] = dPos[bOffset+tid];
        __syncthreads();
        for (int i = 0; i < bdim; ++i)
        {
            gpuAddLJForce( &f, &sPos[tid], &sPos[bdim+i], &boxsize);
        }
        __syncthreads();
    }
    gpuCumulateVec(&dForce[idx], &f);
}


__global__ void gpuLJForcesCut(FLOAT_ARRAY_TYPE *dPos, FLOAT_ARRAY_TYPE *dForce, int N, float boxsize, float rCutSq)
{
    //Reduced units are used -> epsilon=sigma=1
    int tid = threadIdx.x;
    int idx = blockDim.x * blockIdx.x + tid;
    if (idx < N)
    {
        FLOAT_ARRAY_TYPE f;
        f.x = 0.0f;
        f.y = 0.0f;
        f.z = 0.0f;
        for (int i = 0; i < N; ++i)
        {
            if ( i != idx) //excude self interaction
            {
                gpuAddLJForceCut( &f, &dPos[idx], &dPos[i], &boxsize, &rCutSq);
            }
        }
        gpuCumulateVec(&dForce[idx], &f);
    }
}


__global__ void gpuLJForcesCutVar1(FLOAT_ARRAY_TYPE *dPos, FLOAT_ARRAY_TYPE *dForce, int N, float boxsize, float rCutSq)
{
    //Reduced units are used -> epsilon=sigma=1
    int tid = threadIdx.x;
    int idx = blockDim.x * blockIdx.x + tid;
    if (idx < N)
    {
        FLOAT_ARRAY_TYPE tPos = dPos[idx]; //reference particle stored in register
        FLOAT_ARRAY_TYPE f;
        f.x = 0.0f;
        f.y = 0.0f;
        f.z = 0.0f;
        for (int i = 0; i < N; ++i)
        {
            if ( i != idx) //excude self interaction
            {
                gpuAddLJForceCut( &f, &tPos, &dPos[i], &boxsize, &rCutSq);
            }
        }
        gpuCumulateVec(&dForce[idx], &f);
    }
}


__global__ void gpuLJForcesCutRSmem(FLOAT_ARRAY_TYPE *dPos, FLOAT_ARRAY_TYPE *dForce, int N, float boxsize, float rCutSq)
{
    //Reduced units are used -> epsilon=sigma=1
    int tid = threadIdx.x;
    int offset = blockDim.x * blockIdx.x;
    int rowOffset = offset + blockDim.x;
    int idx = offset + tid;
    if(idx < N)
    {
        //load particle positions into shared memory
        extern __shared__ FLOAT_ARRAY_TYPE sPos[];
        sPos[tid] = dPos[idx];

        //sync such that everything is in shared memory before going on
        __syncthreads();

        //calculate forces for particles which are already in shared memory
        FLOAT_ARRAY_TYPE f;
        f.x = 0.;
        f.y = 0.;
        f.z = 0.;

        //block before thread block in device memory
        for (int i = 0; i < offset; ++i)
        {
            gpuAddLJForceCut( &f, &sPos[tid], &dPos[i], &boxsize, &rCutSq);
        }
        //block in shared memory
        for (int i = 0; i < blockDim.x; ++i)
        {
            if ( i != tid) //excude self interaction
            {
                gpuAddLJForceCut( &f, &sPos[tid], &sPos[i], &boxsize, &rCutSq);
            }
        }
        //block after threadblock in device memory
        for (int i = rowOffset; i < N; ++i)
        {
            gpuAddLJForceCut( &f, &sPos[tid], &dPos[i], &boxsize, &rCutSq);
        }
        gpuCumulateVec(&dForce[idx], &f);
    }
}



__global__ void gpuLJForcesCutSmem(FLOAT_ARRAY_TYPE *dPos, FLOAT_ARRAY_TYPE *dForce, int N, float boxsize, float rCutSq)
{
    //Reduced units are used -> epsilon=sigma=1
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + tid;
    extern __shared__ FLOAT_ARRAY_TYPE sPos[];
    if(idx < N)
    {
        FLOAT_ARRAY_TYPE tPos = dPos[idx];//save reference particle position in register
        FLOAT_ARRAY_TYPE f;
        f.x = 0.;
        f.y = 0.;
        f.z = 0.;
        for(int j=0; j<gridDim.x; ++j)
        {
            int bOffset = j * blockDim.x;
            sPos[tid] = dPos[bOffset+tid]; //fetch block of particle positions and save in shared memory
            __syncthreads();
            for (int i=0; i<blockDim.x; ++i)
            {
                int id = bOffset + i;
                if ( id != idx && id < N)
                {
                    gpuAddLJForceCut( &f, &tPos, &sPos[i], &boxsize, &rCutSq);
                }
            }
            __syncthreads();
        }
        gpuCumulateVec(&dForce[idx], &f);
    }
}



__global__ void gpuLJForcesCutSmemSmem(FLOAT_ARRAY_TYPE *dPos, FLOAT_ARRAY_TYPE *dForce, int N, float boxsize, float rCutSq)
{
    //Reduced units are used -> epsilon=sigma=1
    int tid = threadIdx.x;
    int bid = blockIdx.x;
    int bOffset;//current block offset
    int bdim = blockDim.x;
    int offset = blockDim.x * blockIdx.x;
    int idx = offset + tid;
    if(idx < N)
    {
        //load particle positions into shared memory
        extern __shared__ FLOAT_ARRAY_TYPE sPos[];
        sPos[tid] = dPos[idx];

        //sync such that everything is in shared memory before going on
        __syncthreads();

        //calculate forces for particles which are already in shared memory
        FLOAT_ARRAY_TYPE f;//temporary force storage
        f.x = 0.;
        f.y = 0.;
        f.z = 0.;

        for (int j = 0; j < bid; ++j)
        {
            bOffset = j * bdim;
            //load block of positions into shared memory
            sPos[blockDim.x+tid] = dPos[bOffset+tid];
            __syncthreads();
            for (int i = 0; i < bdim; ++i)
            {
                gpuAddLJForceCut( &f, &sPos[tid], &sPos[bdim+i], &boxsize, &rCutSq);
            }
            __syncthreads();
        }

        //now do force calculation inside thread block
        for (int i = 0; i < blockDim.x; ++i)
        {
            if ( i != tid )
            {
                gpuAddLJForceCut( &f, &sPos[tid], &sPos[i], &boxsize, &rCutSq);
            }
        }
        __syncthreads();

        //do force calculation for blocks after thread block
        for (int j = bid + 1; j < gridDim.x; ++j)
        {
            bOffset = j * bdim;
            //load block of positions into shared memory
            sPos[bdim+tid] = dPos[bOffset+tid];
            __syncthreads();
            //do force calculation for this block
            for (int i = 0; i < bdim; ++i)
            {
                gpuAddLJForceCut( &f, &sPos[tid], &sPos[bdim+i], &boxsize, &rCutSq);
            }
            __syncthreads();
        }
        gpuCumulateVec(&dForce[idx], &f);
    }
}

__global__ void gpuResetForces(FLOAT_ARRAY_TYPE *dForces, int N)
{
    int idx = blockDim.x * blockIdx.x + threadIdx.x;
    if ( idx < N )
    {
        dForces[idx].x = 0.;
        dForces[idx].y = 0.;
        dForces[idx].z = 0.;
    }
}


__global__ void gpuLJForcesVerletList(FLOAT_ARRAY_TYPE *dPos, FLOAT_ARRAY_TYPE *dForce, int N, int *dVlists, size_t vlists_pitch, int *dVcount, float boxsize, float rCutSq)
{
    int idx = blockDim.x * blockIdx.x + threadIdx.x;

    if ( idx < N )
    {
        FLOAT_ARRAY_TYPE tPos = dPos[idx];
        FLOAT_ARRAY_TYPE f;
        f.x = 0.;
        f.y = 0.;
        f.z = 0.;
        int* dVlist;
        for (int i = 0; i < dVcount[idx]; ++i)
        {
            dVlist = (int*)((char*)dVlists + i * vlists_pitch);
            gpuAddLJForceCut( &f, &tPos, &dPos[dVlist[idx]], &boxsize, &rCutSq);
        }
        gpuCumulateVec(&dForce[idx], &f);
    }
}


__global__ void gpuHarmonicBondForces(FLOAT_ARRAY_TYPE *dPos, FLOAT_ARRAY_TYPE *dForce, int *cnn, size_t pitch, int N, int numNN, float k, float boxsize)
{
    int idx = blockDim.x * blockIdx.x + threadIdx.x;
    if (idx < N)
    {
        FLOAT_ARRAY_TYPE tPos = dPos[idx];
        int pid;
        int *row;
        for (int c = 0; c < numNN; ++c)
        {
            row = (int*)((char*)cnn + c * pitch);
            pid = row[idx];
            if (pid != -1)
            {
                gpuAddHarmonicBondForce(&dForce[idx], &tPos, &dPos[pid], &k, &boxsize);
            }
        }
    }
}


__global__ void gpuHarmonicBondForcesSmem(FLOAT_ARRAY_TYPE *dPos, FLOAT_ARRAY_TYPE *dForce, int *cnn, size_t pitch, int N, int numNN, float k, float boxsize)
{
    extern __shared__ FLOAT_ARRAY_TYPE sPos[];
    int tid = threadIdx.x;
    int idx = blockDim.x * blockIdx.x + tid;
    if (idx < N)
    {
        sPos[tid] = dPos[idx];
        int pid;
        int *row;
        for (int c = 0; c < numNN; ++c)
        {
            row = (int*)((char*)cnn + c * pitch);
            pid = row[idx];
            if (pid != -1)
            {
                gpuAddHarmonicBondForce(&dForce[idx], &sPos[tid], &dPos[pid], &k, &boxsize);
            }
        }
    }
}

---