<commit_before>/*
Highly Optimized Object-Oriented Molecular Dynamics (HOOMD) Open
Source Software License
Copyright (c) 2008 Ames Laboratory Iowa State University
All rights reserved.

Redistribution and use of HOOMD, in source and binary forms, with or
without modification, are permitted, provided that the following
conditions are met:

* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.

* Neither the name of the copyright holder nor the names HOOMD's
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.

Disclaimer

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND
CONTRIBUTORS ``AS IS''  AND ANY EXPRESS OR IMPLIED WARRANTIES,
INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 

IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS  BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
*/

// $Id$
// $URL$

#include "gpu_forces.h"
#include "gpu_pdata.h"

#ifdef WIN32
#include <cassert>
#else
#include <assert.h>
#endif

/*! \file ljforcesum_kernel.cu
	\brief Contains code for the LJ force sum kernel on the GPU
	\details Functions in this file are NOT to be called by anyone not knowing 
		exactly what they are doing. They are designed to be used solely by 
		LJForceComputeGPU.
*/

///////////////////////////////////////// LJ params

//! The last used id generated by gpu_alloc_ljparam_data()
static unsigned int lj_current_id = 0;
//! The id of the last gpu_ljparam_data selected by gpu_select_ljparam_data()
static unsigned int lj_last_selected_id = INT_MAX;

//! Storage of the lj params in constant memory on the GPU
__constant__ gpu_ljparam_data c_ljparams;

//! Allocate ljparams
/*! \post Memory is allocated for the parameters and all values are set to 0.0f;
	\return Pointer to the allocated structure
	\note NOT THREAD SAFE! Only one set of parameters should be allocated on the host, though
		so this shouldn't matter too much.
*/
gpu_ljparam_data *gpu_alloc_ljparam_data()
	{
	// allocate memory
	gpu_ljparam_data *data = (gpu_ljparam_data *)malloc(sizeof(gpu_ljparam_data));
	assert(data);
 
	// clear
	for (int i = 0; i < MAX_NTYPE_PAIRS; i++)
		{
		data->lj1[i] = 0.0f;
		data->lj2[i] = 0.0f;
		}

	data->id = lj_current_id++;
	return data;
	}

/*! \post Allocated memory is freed on the host
	\param ljparams Parameter set to free
*/
void gpu_free_ljparam_data(gpu_ljparam_data *ljparams)
	{
	assert(ljparams);
	free(ljparams);
	}

/*! \post The parameters are copied to the constant memory on the device. If this was the 
	last structure selected (determined by id), then don't copy unless \a force is specified
	\returns error code passed on from cudaMemcpyToSymbol
	
	\param ljparams Parameter set to make active on the device
	\param force Force an update of the constant memory on the gpu.
*/
cudaError_t gpu_select_ljparam_data(gpu_ljparam_data *ljparams, bool force)
	{
	assert(ljparams);
	cudaError_t retval = cudaSuccess;
	if (ljparams->id != lj_last_selected_id || force)
		{
		retval = cudaMemcpyToSymbol(c_ljparams, ljparams, sizeof(gpu_ljparam_data));
		lj_last_selected_id = ljparams->id;
		}
	return retval;
	}

//! Texture for reading particle positions
texture<float4, 1, cudaReadModeElementType> pdata_pos_tex;

//! Kernel for calculating lj forces
/*! This kerenel is called to calculate the lennard-jones forces on all N particles

	\param d_forces Device memory array to write calculated forces to
	\param pdata Particle data on the GPU to calculate forces on
	\param nlist Neigbhor list data on the GPU to use to calculate the forces
	\param r_cutsq Precalculated r_cut*r_cut, where r_cut is the radius beyond which forces are
		set to 0
	\param box Box dimensions used to implement periodic boundary conditions
	
	Developer information:
	Each block will calculate the forces on a block of particles.
	Each thread will calculate the total force on one particle.
	The neighborlist is arranged in columns so that reads are fully coalesced when doing this.
*/
extern "C" __global__ void calcLJForces_kernel(float4 *d_forces, gpu_pdata_arrays pdata, gpu_nlist_array nlist, float r_cutsq, gpu_boxsize box)
	{
	// start by identifying which particle we are to handle
	int pidx = blockIdx.x * blockDim.x + threadIdx.x;
	
	if (pidx >= pdata.N)
		return;
	
	// load in the length of the list
	int n_neigh = nlist.n_neigh[pidx];

	// read in the position of our particle. Texture reads of float4's are faster than global reads on compute 1.0 hardware
	float4 pos = tex1Dfetch(pdata_pos_tex, pidx);
	
	// initialize the force to 0
	float4 force = make_float4(0.0f, 0.0f, 0.0f, 0.0f);

	// loop over neighbors
	for (int neigh_idx = 0; neigh_idx < n_neigh; neigh_idx++)
		{
		int cur_neigh = nlist.list[nlist.pitch*neigh_idx + pidx];
		
		// get the neighbor's position
		float4 neigh_pos = tex1Dfetch(pdata_pos_tex, cur_neigh);
	
		// calculate dr (with periodic boundary conditions)
		float dx = pos.x - neigh_pos.x;
		float dy = pos.y - neigh_pos.y;
		float dz = pos.z - neigh_pos.z;
			
		dx -= box.Lx * rintf(dx * box.Lxinv);
		dy -= box.Ly * rintf(dy * box.Lyinv);
		dz -= box.Lz * rintf(dz * box.Lzinv);
			
		// calculate r^2
		float rsq = dx*dx + dy*dy + dz*dz;
		
		// calculate 1/r^2
		float r2inv;
		if (rsq >= r_cutsq)
			r2inv = 0.0f;
		else
			r2inv = 1.0f / rsq;

		// lookup the coefficients between this combination of particle types
		int typ_pair = __float_as_int(neigh_pos.w) * MAX_NTYPES + __float_as_int(pos.w);
		float lj1 = c_ljparams.lj1[typ_pair];
		float lj2 = c_ljparams.lj2[typ_pair];
	
		// calculate 1/r^6
		float r6inv = r2inv*r2inv*r2inv;
		// calculate the force magnitude / r
		float fforce = r2inv * r6inv * (12.0f * lj1  * r6inv - 6.0f * lj2);
			
		// add up the force vector components
		force.x += dx * fforce;
		force.y += dy * fforce;
		force.z += dz * fforce;
		force.w += r6inv * (lj1 * r6inv - lj2);
		}

	// potential energy per particle must be halved
	force.w *= 0.5f;
	// now that the force calculation is complete, write out the result
	d_forces[pidx] = force;
	}


/*! \param d_forces Device memory to write forces to
	\param pdata Particle data on the GPU to perform the calculation on
	\param box Box dimensions (in GPU format) to use for periodic boundary conditions
	\param nlist Neighbor list stored on the gpu
	\param r_cutsq Precomputed r_cut*r_cut, where r_cut is the radius beyond which the 
		force is set to 0
	\param M Block size to execute
	
	\returns Any error code resulting from the kernel launch
	\note Always returns cudaSuccess in release builds to avoid the cudaThreadSynchronize()
*/
cudaError_t gpu_ljforce_sum(float4 *d_forces, gpu_pdata_arrays *pdata, gpu_boxsize *box, gpu_nlist_array *nlist, float r_cutsq, int M)
	{
	assert(pdata);
	assert(nlist);

    // setup the grid to run the kernel
    dim3 grid( (int)ceil((double)pdata->N/ (double)M), 1, 1);
    dim3 threads(M, 1, 1);

	// bind the texture
	cudaError_t error = cudaBindTexture(0, pdata_pos_tex, pdata->pos, sizeof(float4) * pdata->N);
	if (error != cudaSuccess)
		return error;

    // run the kernel
    calcLJForces_kernel<<< grid, threads >>>(d_forces, *pdata, *nlist, r_cutsq, *box);

	#ifdef NDEBUG
	return cudaSuccess;
	#else
	cudaThreadSynchronize();
	return cudaGetLastError();
	#endif
	}

// vim:syntax=cpp
<commit_msg>This change hopefully fixes the damnd ULF bug. Again.<commit_after>/*
Highly Optimized Object-Oriented Molecular Dynamics (HOOMD) Open
Source Software License
Copyright (c) 2008 Ames Laboratory Iowa State University
All rights reserved.

Redistribution and use of HOOMD, in source and binary forms, with or
without modification, are permitted, provided that the following
conditions are met:

* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.

* Neither the name of the copyright holder nor the names HOOMD's
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.

Disclaimer

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND
CONTRIBUTORS ``AS IS''  AND ANY EXPRESS OR IMPLIED WARRANTIES,
INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 

IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS  BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
*/

// $Id$
// $URL$

#include "gpu_forces.h"
#include "gpu_pdata.h"

#ifdef WIN32
#include <cassert>
#else
#include <assert.h>
#endif

/*! \file ljforcesum_kernel.cu
	\brief Contains code for the LJ force sum kernel on the GPU
	\details Functions in this file are NOT to be called by anyone not knowing 
		exactly what they are doing. They are designed to be used solely by 
		LJForceComputeGPU.
*/

///////////////////////////////////////// LJ params

//! The last used id generated by gpu_alloc_ljparam_data()
static unsigned int lj_current_id = 0;
//! The id of the last gpu_ljparam_data selected by gpu_select_ljparam_data()
static unsigned int lj_last_selected_id = INT_MAX;

//! Storage of the lj params in constant memory on the GPU
__constant__ gpu_ljparam_data c_ljparams;

//! Allocate ljparams
/*! \post Memory is allocated for the parameters and all values are set to 0.0f;
	\return Pointer to the allocated structure
	\note NOT THREAD SAFE! Only one set of parameters should be allocated on the host, though
		so this shouldn't matter too much.
*/
gpu_ljparam_data *gpu_alloc_ljparam_data()
	{
	// allocate memory
	gpu_ljparam_data *data = (gpu_ljparam_data *)malloc(sizeof(gpu_ljparam_data));
	assert(data);
 
	// clear
	for (int i = 0; i < MAX_NTYPE_PAIRS; i++)
		{
		data->lj1[i] = 0.0f;
		data->lj2[i] = 0.0f;
		}

	data->id = lj_current_id++;
	return data;
	}

/*! \post Allocated memory is freed on the host
	\param ljparams Parameter set to free
*/
void gpu_free_ljparam_data(gpu_ljparam_data *ljparams)
	{
	assert(ljparams);
	free(ljparams);
	}

/*! \post The parameters are copied to the constant memory on the device. If this was the 
	last structure selected (determined by id), then don't copy unless \a force is specified
	\returns error code passed on from cudaMemcpyToSymbol
	
	\param ljparams Parameter set to make active on the device
	\param force Force an update of the constant memory on the gpu.
*/
cudaError_t gpu_select_ljparam_data(gpu_ljparam_data *ljparams, bool force)
	{
	assert(ljparams);
	cudaError_t retval = cudaSuccess;
	if (ljparams->id != lj_last_selected_id || force)
		{
		retval = cudaMemcpyToSymbol(c_ljparams, ljparams, sizeof(gpu_ljparam_data));
		lj_last_selected_id = ljparams->id;
		}
	return retval;
	}

//! Texture for reading particle positions
texture<float4, 1, cudaReadModeElementType> pdata_pos_tex;

//! Kernel for calculating lj forces
/*! This kerenel is called to calculate the lennard-jones forces on all N particles

	\param d_forces Device memory array to write calculated forces to
	\param pdata Particle data on the GPU to calculate forces on
	\param nlist Neigbhor list data on the GPU to use to calculate the forces
	\param r_cutsq Precalculated r_cut*r_cut, where r_cut is the radius beyond which forces are
		set to 0
	\param box Box dimensions used to implement periodic boundary conditions
	
	Developer information:
	Each block will calculate the forces on a block of particles.
	Each thread will calculate the total force on one particle.
	The neighborlist is arranged in columns so that reads are fully coalesced when doing this.
*/
extern "C" __global__ void calcLJForces_kernel(float4 *d_forces, gpu_pdata_arrays pdata, gpu_nlist_array nlist, float r_cutsq, gpu_boxsize box)
	{
	// start by identifying which particle we are to handle
	int pidx = blockIdx.x * blockDim.x + threadIdx.x;
	
	if (pidx >= pdata.N)
		return;
	
	// load in the length of the list
	int n_neigh = nlist.n_neigh[pidx];

	// read in the position of our particle. Texture reads of float4's are faster than global reads on compute 1.0 hardware
	float4 pos = tex1Dfetch(pdata_pos_tex, pidx);
	
	// initialize the force to 0
	float4 force = make_float4(0.0f, 0.0f, 0.0f, 0.0f);

	// loop over neighbors
	for (int neigh_idx = 0; neigh_idx < nlist.height; neigh_idx++)
		{
		if (neigh_idx < n_neigh)
		{
		int cur_neigh = nlist.list[nlist.pitch*neigh_idx + pidx];
		
		// get the neighbor's position
		float4 neigh_pos = tex1Dfetch(pdata_pos_tex, cur_neigh);
	
		// calculate dr (with periodic boundary conditions)
		float dx = pos.x - neigh_pos.x;
		float dy = pos.y - neigh_pos.y;
		float dz = pos.z - neigh_pos.z;
			
		dx -= box.Lx * rintf(dx * box.Lxinv);
		dy -= box.Ly * rintf(dy * box.Lyinv);
		dz -= box.Lz * rintf(dz * box.Lzinv);
			
		// calculate r^2
		float rsq = dx*dx + dy*dy + dz*dz;
		
		// calculate 1/r^2
		float r2inv;
		if (rsq >= r_cutsq)
			r2inv = 0.0f;
		else
			r2inv = 1.0f / rsq;

		// lookup the coefficients between this combination of particle types
		int typ_pair = __float_as_int(neigh_pos.w) * MAX_NTYPES + __float_as_int(pos.w);
		float lj1 = c_ljparams.lj1[typ_pair];
		float lj2 = c_ljparams.lj2[typ_pair];
	
		// calculate 1/r^6
		float r6inv = r2inv*r2inv*r2inv;
		// calculate the force magnitude / r
		float fforce = r2inv * r6inv * (12.0f * lj1  * r6inv - 6.0f * lj2);
			
		// add up the force vector components
		force.x += dx * fforce;
		force.y += dy * fforce;
		force.z += dz * fforce;
		force.w += r6inv * (lj1 * r6inv - lj2);
		}
		}

	// potential energy per particle must be halved
	force.w *= 0.5f;
	// now that the force calculation is complete, write out the result
	d_forces[pidx] = force;
	}


/*! \param d_forces Device memory to write forces to
	\param pdata Particle data on the GPU to perform the calculation on
	\param box Box dimensions (in GPU format) to use for periodic boundary conditions
	\param nlist Neighbor list stored on the gpu
	\param r_cutsq Precomputed r_cut*r_cut, where r_cut is the radius beyond which the 
		force is set to 0
	\param M Block size to execute
	
	\returns Any error code resulting from the kernel launch
	\note Always returns cudaSuccess in release builds to avoid the cudaThreadSynchronize()
*/
cudaError_t gpu_ljforce_sum(float4 *d_forces, gpu_pdata_arrays *pdata, gpu_boxsize *box, gpu_nlist_array *nlist, float r_cutsq, int M)
	{
	assert(pdata);
	assert(nlist);

    // setup the grid to run the kernel
    dim3 grid( (int)ceil((double)pdata->N/ (double)M), 1, 1);
    dim3 threads(M, 1, 1);

	// bind the texture
	cudaError_t error = cudaBindTexture(0, pdata_pos_tex, pdata->pos, sizeof(float4) * pdata->N);
	if (error != cudaSuccess)
		return error;

    // run the kernel
    calcLJForces_kernel<<< grid, threads >>>(d_forces, *pdata, *nlist, r_cutsq, *box);

	#ifdef NDEBUG
	return cudaSuccess;
	#else
	cudaThreadSynchronize();
	return cudaGetLastError();
	#endif
	}

// vim:syntax=cpp

{'fraction_non_alphanumeric': 0.0837864557587741, 'fraction_numerical': 0.01050420168067227, 'mean_word_length': 3.8986077481840193, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7531939', 'n_tokens_mistral': 5581, 'n_tokens_neox': 5163, 'n_words': 2033}