diff --git a/lib/gpu/Nvidia.makefile b/lib/gpu/Nvidia.makefile index 4a4166458b..5876017bef 100644 --- a/lib/gpu/Nvidia.makefile +++ b/lib/gpu/Nvidia.makefile @@ -50,6 +50,7 @@ OBJS = $(OBJ_DIR)/pair_gpu_atom.o $(OBJ_DIR)/pair_gpu_nbor.o \ $(OBJ_DIR)/crml_gpu_memory.o $(OBJ_DIR)/crml_gpu.o \ $(OBJ_DIR)/cmm_cut_gpu_memory.o $(OBJ_DIR)/cmm_cut_gpu.o \ $(OBJ_DIR)/cmmc_long_gpu_memory.o $(OBJ_DIR)/cmmc_long_gpu.o \ + $(OBJ_DIR)/cmmc_msm_gpu_memory.o $(OBJ_DIR)/cmmc_msm_gpu.o \ $(CUDPP) PTXS = $(OBJ_DIR)/pair_gpu_atom_kernel.ptx $(OBJ_DIR)/pair_gpu_atom_ptx.h \ $(OBJ_DIR)/pair_gpu_nbor_kernel.ptx $(OBJ_DIR)/pair_gpu_nbor_ptx.h \ @@ -62,7 +63,8 @@ PTXS = $(OBJ_DIR)/pair_gpu_atom_kernel.ptx $(OBJ_DIR)/pair_gpu_atom_ptx.h \ $(OBJ_DIR)/ljcl_cut_gpu_kernel.ptx $(OBJ_DIR)/ljcl_cut_gpu_ptx.h \ $(OBJ_DIR)/crml_gpu_kernel.ptx $(OBJ_DIR)/crml_gpu_ptx.h \ $(OBJ_DIR)/cmm_cut_gpu_kernel.ptx $(OBJ_DIR)/cmm_cut_gpu_ptx.h \ - $(OBJ_DIR)/cmmc_long_gpu_kernel.ptx $(OBJ_DIR)/cmmc_long_gpu_ptx.h + $(OBJ_DIR)/cmmc_long_gpu_kernel.ptx $(OBJ_DIR)/cmmc_long_gpu_ptx.h \ + $(OBJ_DIR)/cmmc_msm_gpu_kernel.ptx $(OBJ_DIR)/cmmc_msm_gpu_ptx.h all: $(GPU_LIB) $(EXECS) @@ -219,6 +221,18 @@ $(OBJ_DIR)/cmmc_long_gpu_memory.o: $(ALL_H) cmmc_long_gpu_memory.h cmmc_long_gpu $(OBJ_DIR)/cmmc_long_gpu.o: $(ALL_H) cmmc_long_gpu_memory.h cmmc_long_gpu.cpp $(CUDR) -o $@ -c cmmc_long_gpu.cpp -I$(OBJ_DIR) +$(OBJ_DIR)/cmmc_msm_gpu_kernel.ptx: cmmc_msm_gpu_kernel.cu pair_gpu_precision.h + $(CUDA) --ptx -DNV_KERNEL -o $@ cmmc_msm_gpu_kernel.cu + +$(OBJ_DIR)/cmmc_msm_gpu_ptx.h: $(OBJ_DIR)/cmmc_msm_gpu_kernel.ptx $(OBJ_DIR)/cmmc_msm_gpu_kernel.ptx + $(BSH) ./geryon/file_to_cstr.sh $(OBJ_DIR)/cmmc_msm_gpu_kernel.ptx $(OBJ_DIR)/cmmc_msm_gpu_ptx.h + +$(OBJ_DIR)/cmmc_msm_gpu_memory.o: $(ALL_H) cmmc_msm_gpu_memory.h cmmc_msm_gpu_memory.cpp $(OBJ_DIR)/cmmc_msm_gpu_ptx.h $(OBJ_DIR)/atomic_gpu_memory.o + $(CUDR) -o $@ -c cmmc_msm_gpu_memory.cpp -I$(OBJ_DIR) + +$(OBJ_DIR)/cmmc_msm_gpu.o: $(ALL_H) cmmc_msm_gpu_memory.h cmmc_msm_gpu.cpp + $(CUDR) -o $@ -c cmmc_msm_gpu.cpp -I$(OBJ_DIR) + $(BIN_DIR)/nvc_get_devices: ./geryon/ucl_get_devices.cpp $(NVC_H) $(CUDR) -o $@ ./geryon/ucl_get_devices.cpp -DUCL_CUDART $(CUDA_LINK) diff --git a/lib/gpu/cmmc_msm_gpu.cpp b/lib/gpu/cmmc_msm_gpu.cpp new file mode 100644 index 0000000000..c81f3a53dd --- /dev/null +++ b/lib/gpu/cmmc_msm_gpu.cpp @@ -0,0 +1,130 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Mike Brown (ORNL), brownw@ornl.gov +------------------------------------------------------------------------- */ + +#include +#include +#include + +#include "cmmc_msm_gpu_memory.h" + +using namespace std; + +static CMMM_GPU_Memory CMMMMF; + +// --------------------------------------------------------------------------- +// Allocate memory on host and device and copy constants to device +// --------------------------------------------------------------------------- +bool cmmm_gpu_init(const int ntypes, double **cutsq, int **cg_type, + double **host_lj1, double **host_lj2, double **host_lj3, + double **host_lj4, double **offset, double *special_lj, + const int inum, const int nall, const int max_nbors, + const int maxspecial, const double cell_size, int &gpu_mode, + FILE *screen, double **host_cut_ljsq, double host_cut_coulsq, + double *host_special_coul, const double qqrd2e, + const int smooth) { + CMMMMF.clear(); + gpu_mode=CMMMMF.device->gpu_mode(); + double gpu_split=CMMMMF.device->particle_split(); + int first_gpu=CMMMMF.device->first_device(); + int last_gpu=CMMMMF.device->last_device(); + int world_me=CMMMMF.device->world_me(); + int gpu_rank=CMMMMF.device->gpu_rank(); + int procs_per_gpu=CMMMMF.device->procs_per_gpu(); + + CMMMMF.device->init_message(screen,"cg/cmm/coul/msm",first_gpu,last_gpu); + + bool message=false; + if (CMMMMF.device->replica_me()==0 && screen) + message=true; + + if (message) { + fprintf(screen,"Initializing GPU and compiling on process 0..."); + fflush(screen); + } + + if (world_me==0) { + bool init_ok=CMMMMF.init(ntypes, cutsq, cg_type, host_lj1, host_lj2, + host_lj3, host_lj4, offset, special_lj, inum, + nall, 300, maxspecial, cell_size, gpu_split, + screen, host_cut_ljsq, host_cut_coulsq, + host_special_coul, qqrd2e,smooth); + if (!init_ok) + return false; + } + + CMMMMF.device->world_barrier(); + if (message) + fprintf(screen,"Done.\n"); + + for (int i=0; igpu_barrier(); + if (message) + fprintf(screen,"Done.\n"); + } + if (message) + fprintf(screen,"\n"); + return true; +} + +void cmmm_gpu_clear() { + CMMMMF.clear(); +} + +int * cmmm_gpu_compute_n(const int timestep, const int ago, const int inum_full, + const int nall, double **host_x, int *host_type, + double *boxlo, double *boxhi, int *tag, int **nspecial, + int **special, const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + const double cpu_time, bool &success, double *host_q) { + return CMMMMF.compute(timestep, ago, inum_full, nall, host_x, host_type, boxlo, + boxhi, tag, nspecial, special, eflag, vflag, eatom, + vatom, host_start, cpu_time, success, host_q); +} + +void cmmm_gpu_compute(const int timestep, const int ago, const int inum_full, + const int nall, double **host_x, int *host_type, + int *ilist, int *numj, int **firstneigh, + const bool eflag, const bool vflag, const bool eatom, + const bool vatom, int &host_start, const double cpu_time, + bool &success, double *host_q) { + CMMMMF.compute(timestep,ago,inum_full,nall,host_x,host_type,ilist,numj, + firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success, + host_q); +} + +double cmmm_gpu_bytes() { + return CMMMMF.host_memory_usage(); +} + + diff --git a/lib/gpu/cmmc_msm_gpu_kernel.cu b/lib/gpu/cmmc_msm_gpu_kernel.cu new file mode 100644 index 0000000000..316e1cfa9e --- /dev/null +++ b/lib/gpu/cmmc_msm_gpu_kernel.cu @@ -0,0 +1,422 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Mike Brown (ORNL), brownw@ornl.gov +------------------------------------------------------------------------- */ + +#ifndef CMMM_GPU_KERNEL +#define CMMM_GPU_KERNEL + +#define MAX_SHARED_TYPES 8 + +#ifdef _DOUBLE_DOUBLE +#define numtyp double +#define numtyp2 double2 +#define numtyp4 double4 +#define acctyp double +#define acctyp4 double4 +#endif + +#ifdef _SINGLE_DOUBLE +#define numtyp float +#define numtyp2 float2 +#define numtyp4 float4 +#define acctyp double +#define acctyp4 double4 +#endif + +#ifndef numtyp +#define numtyp float +#define numtyp2 float2 +#define numtyp4 float4 +#define acctyp float +#define acctyp4 float4 +#endif + +#ifdef NV_KERNEL + +#include "geryon/ucl_nv_kernel.h" +texture pos_tex; +texture q_tex; + +#ifdef _DOUBLE_DOUBLE +__inline double4 fetch_pos(const int& i, const double4 *pos) +{ + return pos[i]; +} +__inline double fetch_q(const int& i, const double *q) +{ + return q[i]; +} +#else +__inline float4 fetch_pos(const int& i, const float4 *pos) +{ + return tex1Dfetch(pos_tex, i); +} +__inline float fetch_q(const int& i, const float *q) +{ + return tex1Dfetch(q_tex, i); +} +#endif + +#else + +#pragma OPENCL EXTENSION cl_khr_fp64: enable +#define GLOBAL_ID_X get_global_id(0) +#define THREAD_ID_X get_local_id(0) +#define BLOCK_ID_X get_group_id(0) +#define BLOCK_SIZE_X get_local_size(0) +#define __syncthreads() barrier(CLK_LOCAL_MEM_FENCE) +#define __inline inline + +#define fetch_pos(i,y) x_[i] +#define fetch_q(i,y) q_[i] + +#endif + +__kernel void kernel_pair(__global numtyp4 *x_, __global numtyp4 *lj1, + __global numtyp4* lj3, const int lj_types, + __global numtyp *sp_lj_in, __global int *dev_nbor, + __global acctyp4 *ans, __global acctyp *engv, + const int eflag, const int vflag, const int inum, + const int nall, const int nbor_pitch, + __global numtyp *q_ , const numtyp cut_coulsq, + const numtyp qqrd2e, const int smooth) { + // ii indexes the two interacting particles in gi + int ii=GLOBAL_ID_X; + __local numtyp sp_lj[8]; + sp_lj[0]=sp_lj_in[0]; + sp_lj[1]=sp_lj_in[1]; + sp_lj[2]=sp_lj_in[2]; + sp_lj[3]=sp_lj_in[3]; + sp_lj[4]=sp_lj_in[4]; + sp_lj[5]=sp_lj_in[5]; + sp_lj[6]=sp_lj_in[6]; + sp_lj[7]=sp_lj_in[7]; + __local numtyp _ia; + __local numtyp _ia2; + __local numtyp _ia3; + + if (ii0) { + if (rsq < cut_coulsq) + if (smooth==0) + e_coul += prefactor*(ir+_ia*((numtyp)2.1875-(numtyp)2.1875*r2_ia2+ + (numtyp)1.3125*r4_ia4- + (numtyp)0.3125*r4_ia4*r2_ia2)- + factor_coul*ir); + else + e_coul += prefactor*(ir+_ia*((numtyp)2.4609375-(numtyp)3.28125* + r2_ia2+(numtyp)2.953125*r4_ia4- + (numtyp)1.40625*r6_ia6+ + (numtyp)0.2734375*r4_ia4*r4_ia4)); + + if (rsq < lj1[mtype].y) { + energy += factor_lj*inv1*(lj3[mtype].y*inv2-lj3[mtype].z)- + lj3[mtype].w; + } + } + if (vflag>0) { + virial[0] += delx*delx*force; + virial[1] += dely*dely*force; + virial[2] += delz*delz*force; + virial[3] += delx*dely*force; + virial[4] += delx*delz*force; + virial[5] += dely*delz*force; + } + } + + } // for nbor + + // Store answers + __global acctyp *ap1=engv+ii; + if (eflag>0) { + *ap1=energy; + ap1+=inum; + *ap1=e_coul; + ap1+=inum; + } + if (vflag>0) { + for (int i=0; i<6; i++) { + *ap1=virial[i]; + ap1+=inum; + } + } + ans[ii]=f; + } // if ii +} + +__kernel void kernel_pair_fast(__global numtyp4 *x_, __global numtyp4 *lj1_in, + __global numtyp4* lj3_in, + __global numtyp* sp_lj_in, __global int *dev_nbor, + __global acctyp4 *ans, __global acctyp *engv, + const int eflag, const int vflag, const int inum, + const int nall, const int nbor_pitch, + __global numtyp *q_ , const numtyp cut_coulsq, + const numtyp qqrd2e, const int smooth) { + // ii indexes the two interacting particles in gi + int ii=THREAD_ID_X; + __local numtyp4 lj1[MAX_SHARED_TYPES*MAX_SHARED_TYPES]; + __local numtyp4 lj3[MAX_SHARED_TYPES*MAX_SHARED_TYPES]; + __local numtyp sp_lj[8]; + if (ii<8) + sp_lj[ii]=sp_lj_in[ii]; + if (ii0) { + if (rsq < cut_coulsq) + if (smooth==0) + e_coul += prefactor*(ir+_ia*((numtyp)2.1875-(numtyp)2.1875*r2_ia2+ + (numtyp)1.3125*r4_ia4- + (numtyp)0.3125*r4_ia4*r2_ia2)- + factor_coul*ir); + else + e_coul += prefactor*(ir+_ia*((numtyp)2.4609375-(numtyp)3.28125* + r2_ia2+(numtyp)2.953125*r4_ia4- + (numtyp)1.40625*r6_ia6+ + (numtyp)0.2734375*r4_ia4*r4_ia4)); + if (rsq < lj1[mtype].y) { + energy += factor_lj*inv1*(lj3[mtype].y*inv2-lj3[mtype].z)- + lj3[mtype].w; + } + } + if (vflag>0) { + virial[0] += delx*delx*force; + virial[1] += dely*dely*force; + virial[2] += delz*delz*force; + virial[3] += delx*dely*force; + virial[4] += delx*delz*force; + virial[5] += dely*delz*force; + } + } + + } // for nbor + + // Store answers + __global acctyp *ap1=engv+ii; + if (eflag>0) { + *ap1=energy; + ap1+=inum; + *ap1=e_coul; + ap1+=inum; + } + if (vflag>0) { + for (int i=0; i<6; i++) { + *ap1=virial[i]; + ap1+=inum; + } + } + ans[ii]=f; + } // if ii*/ +} + +#endif + diff --git a/lib/gpu/cmmc_msm_gpu_memory.cpp b/lib/gpu/cmmc_msm_gpu_memory.cpp new file mode 100644 index 0000000000..6fced81d86 --- /dev/null +++ b/lib/gpu/cmmc_msm_gpu_memory.cpp @@ -0,0 +1,164 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Mike Brown (ORNL), brownw@ornl.gov +------------------------------------------------------------------------- */ + +#ifdef USE_OPENCL +#include "cmmc_msm_gpu_cl.h" +#else +#include "cmmc_msm_gpu_ptx.h" +#endif + +#include "cmmc_msm_gpu_memory.h" +#include +#define CMMM_GPU_MemoryT CMMM_GPU_Memory + +extern PairGPUDevice pair_gpu_device; + +template +CMMM_GPU_MemoryT::CMMM_GPU_Memory() : ChargeGPUMemory(), + _allocated(false) { +} + +template +CMMM_GPU_MemoryT::~CMMM_GPU_Memory() { + clear(); +} + +template +int CMMM_GPU_MemoryT::bytes_per_atom(const int max_nbors) const { + return this->bytes_per_atom_atomic(max_nbors); +} + +template +bool CMMM_GPU_MemoryT::init(const int ntypes, double **host_cutsq, + int **host_cg_type, double **host_lj1, + double **host_lj2, double **host_lj3, + double **host_lj4, double **host_offset, + double *host_special_lj, const int nlocal, + const int nall, const int max_nbors, + const int maxspecial, const double cell_size, + const double gpu_split, FILE *_screen, + double **host_cut_ljsq, + const double host_cut_coulsq, + double *host_special_coul, const double qqrd2e, + const int smooth) { + this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,gpu_split, + _screen,cmmc_msm_gpu_kernel); + + // If atom type constants fit in shared memory use fast kernel + int lj_types=ntypes; + shared_types=false; + if (lj_types<=MAX_SHARED_TYPES && this->_block_size>=MAX_SHARED_TYPES) { + lj_types=MAX_SHARED_TYPES; + shared_types=true; + } + _lj_types=lj_types; + + // Allocate a host write buffer for data initialization + UCL_H_Vec host_write(lj_types*lj_types*32,*(this->ucl_device), + UCL_WRITE_OPTIMIZED); + + for (int i=0; iucl_device),UCL_READ_ONLY); + this->atom->type_pack4(ntypes,lj_types,lj1,host_write,host_cutsq, + host_cut_ljsq,host_lj1,host_lj2); + + lj3.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY); + this->atom->type_pack4(ntypes,lj_types,lj3,host_write,host_cg_type,host_lj3, + host_lj4,host_offset); + + sp_lj.alloc(8,*(this->ucl_device),UCL_READ_ONLY); + for (int i=0; i<4; i++) { + host_write[i]=host_special_lj[i]; + host_write[i+4]=host_special_coul[i]; + } + ucl_copy(sp_lj,host_write,8,false); + + _cut_coulsq=host_cut_coulsq; + _qqrd2e=qqrd2e; + _smooth=smooth; + + _allocated=true; + this->_max_bytes=lj1.row_bytes()+lj3.row_bytes()+sp_lj.row_bytes(); + return true; +} + +template +void CMMM_GPU_MemoryT::clear() { + if (!_allocated) + return; + _allocated=false; + + lj1.clear(); + lj3.clear(); + sp_lj.clear(); + this->clear_atomic(); +} + +template +double CMMM_GPU_MemoryT::host_memory_usage() const { + return this->host_memory_usage_atomic()+sizeof(CMMM_GPU_Memory); +} + +// --------------------------------------------------------------------------- +// Calculate energies, forces, and torques +// --------------------------------------------------------------------------- +template +void CMMM_GPU_MemoryT::loop(const bool _eflag, const bool _vflag) { + // Compute the block size and grid size to keep all cores busy + const int BX=this->block_size(); + int eflag, vflag; + if (_eflag) + eflag=1; + else + eflag=0; + + if (_vflag) + vflag=1; + else + vflag=0; + + int GX=static_cast(ceil(static_cast(this->atom->inum())/BX)); + + int ainum=this->atom->inum(); + int anall=this->atom->nall(); + int nbor_pitch=this->nbor->nbor_pitch(); + this->time_pair.start(); + if (shared_types) { + this->k_pair_fast.set_size(GX,BX); + this->k_pair_fast.run(&this->atom->dev_x.begin(), &lj1.begin(), + &lj3.begin(), &sp_lj.begin(), + &this->nbor->dev_nbor.begin(), + &this->atom->dev_ans.begin(), + &this->atom->dev_engv.begin(), &eflag, &vflag, + &ainum, &anall, &nbor_pitch, + &this->atom->dev_q.begin(), &_cut_coulsq, + &_qqrd2e, &_smooth); + } else { + this->k_pair.set_size(GX,BX); + this->k_pair.run(&this->atom->dev_x.begin(), &lj1.begin(), &lj3.begin(), + &_lj_types, &sp_lj.begin(), &this->nbor->dev_nbor.begin(), + &this->atom->dev_ans.begin(), + &this->atom->dev_engv.begin(), &eflag, &vflag, &ainum, + &anall, &nbor_pitch, &this->atom->dev_q.begin(), + &_cut_coulsq, &_qqrd2e, &_smooth); + } + this->time_pair.stop(); +} + +template class CMMM_GPU_Memory; diff --git a/lib/gpu/cmmc_msm_gpu_memory.h b/lib/gpu/cmmc_msm_gpu_memory.h new file mode 100644 index 0000000000..ddaef5944d --- /dev/null +++ b/lib/gpu/cmmc_msm_gpu_memory.h @@ -0,0 +1,76 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Mike Brown (ORNL), brownw@ornl.gov +------------------------------------------------------------------------- */ + +#ifndef CMMM_GPU_MEMORY_H +#define CMMM_GPU_MEMORY_H + +#include "charge_gpu_memory.h" + +template +class CMMM_GPU_Memory : public ChargeGPUMemory { + public: + CMMM_GPU_Memory(); + ~CMMM_GPU_Memory(); + + /// Clear any previous data and set up for a new LAMMPS run + /** \param max_nbors initial number of rows in the neighbor matrix + * \param cell_size cutoff + skin + * \param gpu_split fraction of particles handled by device **/ + bool init(const int ntypes, double **host_cutsq, int ** cg_type, + double **host_lj1, double **host_lj2, double **host_lj3, + double **host_lj4, double **host_offset, double *host_special_lj, + const int nlocal, const int nall, const int max_nbors, + const int maxspecial, const double cell_size, + const double gpu_split, FILE *screen, double **host_cut_ljsq, + const double host_cut_coulsq, double *host_special_coul, + const double qqrd2e, const int smooth); + + /// Clear all host and device data + /** \note This is called at the beginning of the init() routine **/ + void clear(); + + /// Returns memory usage on device per atom + int bytes_per_atom(const int max_nbors) const; + + /// Total host memory used by library for pair style + double host_memory_usage() const; + + // --------------------------- TYPE DATA -------------------------- + + /// lj1.x = cutsq, lj1.y = cutsq_vdw, lj1.z = lj1, lj1.w = lj2, + UCL_D_Vec lj1; + /// lj3.x = cg_type, lj3.y = lj3, lj3.z = lj4, lj3.w = offset + UCL_D_Vec lj3; + /// Special LJ values [0-3] and Special Coul values [4-7] + UCL_D_Vec sp_lj; + + /// If atom type constants fit in shared memory, use fast kernels + bool shared_types; + + /// Number of atom types + int _lj_types; + + numtyp _cut_coulsq, _qqrd2e; + + private: + bool _allocated; + int _smooth; + void loop(const bool _eflag, const bool _vflag); +}; + +#endif + diff --git a/lib/gpu/pair_gpu_device.cpp b/lib/gpu/pair_gpu_device.cpp index c2d980cf99..d02618e22c 100644 --- a/lib/gpu/pair_gpu_device.cpp +++ b/lib/gpu/pair_gpu_device.cpp @@ -19,6 +19,9 @@ #include "pair_gpu_precision.h" #include #include +#ifdef _OPENMP +#include +#endif #define PairGPUDeviceT PairGPUDevice @@ -39,6 +42,9 @@ bool PairGPUDeviceT::init_device(MPI_Comm world, MPI_Comm replica, const int gpu_mode, const double p_split, const int nthreads) { _nthreads=nthreads; + #ifdef _OPENMP + omp_set_num_threads(nthreads); + #endif if (_device_init) return true; @@ -153,6 +159,9 @@ void PairGPUDeviceT::init_message(FILE *screen, const char *name, fprintf(screen,"-------------------------------------\n"); fprintf(screen,"- Using GPGPU acceleration for %s:\n",name); fprintf(screen,"- with %d procs per device.\n",_procs_per_gpu); + #ifdef _OPENMP + fprintf(screen,"- with %d threads per proc.\n",_nthreads); + #endif fprintf(screen,"-------------------------------------"); fprintf(screen,"-------------------------------------\n"); diff --git a/src/GPU/Install.sh b/src/GPU/Install.sh index fc36233e2a..a702c50e57 100644 --- a/src/GPU/Install.sh +++ b/src/GPU/Install.sh @@ -14,6 +14,11 @@ if (test $1 = 1) then sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(gpu_SYSLIB) |' ../Makefile.package fi + if (test -e ../pppm.cpp) then + cp pppm_gpu.cpp .. + cp pppm_gpu.h .. + fi + if (test -e ../pair_gayberne.cpp) then cp pair_gayberne_gpu.cpp .. cp pair_gayberne_gpu.h .. @@ -37,18 +42,27 @@ if (test $1 = 1) then if (test -e ../pair_cg_cmm_coul_long.cpp) then cp pair_cg_cmm_coul_long_gpu.cpp .. cp pair_cg_cmm_coul_long_gpu.h .. + cp pair_cg_cmm_coul_msm.cpp .. + cp pair_cg_cmm_coul_msm.h .. + cp pair_cg_cmm_coul_msm_gpu.cpp .. + cp pair_cg_cmm_coul_msm_gpu.h .. fi cp pair_lj_cut_gpu.cpp .. cp pair_lj96_cut_gpu.cpp .. cp pair_lj_cut_coul_cut_gpu.cpp .. + cp pair_lj_cut_tgpu.cpp .. cp pair_lj_cut_gpu.h .. cp pair_lj96_cut_gpu.h .. cp pair_lj_cut_coul_cut_gpu.h .. + cp pair_lj_cut_tgpu.h .. cp fix_gpu.cpp .. cp fix_gpu.h .. + cp pair_omp_gpu.cpp .. + cp pair_omp_gpu.h .. + elif (test $1 = 0) then if (test -e ../Makefile.package) then @@ -56,15 +70,20 @@ elif (test $1 = 0) then sed -i -e 's/[^ \t]*gpu_[^ \t]*) //' ../Makefile.package fi + rm ../pppm_gpu.cpp rm ../pair_gayberne_gpu.cpp rm ../pair_lj_cut_gpu.cpp rm ../pair_lj96_cut_gpu.cpp rm ../pair_lj_cut_coul_cut_gpu.cpp rm ../pair_lj_cut_coul_long_gpu.cpp rm ../pair_lj_charmm_coul_long_gpu.cpp + rm ../pair_lj_cut_tgpu.cpp rm ../pair_cg_cmm_gpu.cpp rm ../pair_cg_cmm_coul_long_gpu.cpp + rm ../pair_cg_cmm_coul_msm.cpp + rm ../pair_cg_cmm_coul_msm_gpu.cpp rm ../fix_gpu.cpp + rm ../pair_omp_gpu.cpp rm ../pppm_gpu.h rm ../pair_gayberne_gpu.h @@ -73,9 +92,13 @@ elif (test $1 = 0) then rm ../pair_lj_cut_coul_cut_gpu.h rm ../pair_lj_cut_coul_long_gpu.h rm ../pair_lj_charmm_coul_long_gpu.h + rm ../pair_lj_cut_tgpu.cpp rm ../pair_cg_cmm_gpu.h rm ../pair_cg_cmm_coul_long_gpu.h + rm ../pair_cg_cmm_coul_msm.h + rm ../pair_cg_cmm_coul_msm_gpu.h rm ../fix_gpu.h + rm ../pair_omp_gpu.h fi diff --git a/src/GPU/fix_gpu.cpp b/src/GPU/fix_gpu.cpp index 75ce1e83f3..4f148bbb5f 100644 --- a/src/GPU/fix_gpu.cpp +++ b/src/GPU/fix_gpu.cpp @@ -42,7 +42,7 @@ extern double lmp_gpu_forces(double **f, double **tor, double *eatom, FixGPU::FixGPU(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { - if (narg != 7) error->all("Illegal fix gpu command"); + if (narg < 7) error->all("Illegal fix gpu command"); if (strcmp(arg[1],"all") != 0) error->all("Illegal fix gpu command"); @@ -66,8 +66,25 @@ FixGPU::FixGPU(LAMMPS *lmp, int narg, char **arg) : if (particle_split==0 || particle_split>1) error->all("Illegal fix gpu command."); + int nthreads = 1; + if (narg == 9) { + if (strcmp(arg[7],"nthreads") == 0) + nthreads = atoi(arg[8]); + else + error->all("Illegal fix gpu command."); + } else if (narg != 7) + error->all("Illegal fix gpu command."); + + if (nthreads < 1) + error->all("Illegal fix gpu command."); + + #ifndef _OPENMP + if (nthreads > 1) + error->all("No OpenMP support compiled in."); + #endif + if (!lmp_init_device(universe->uworld,world,first_gpu,last_gpu,gpu_mode, - particle_split,1)) + particle_split,nthreads)) error->one("Could not find or initialize a specified accelerator device."); } diff --git a/src/GPU/pair_cg_cmm_coul_msm.cpp b/src/GPU/pair_cg_cmm_coul_msm.cpp new file mode 100644 index 0000000000..139f125513 --- /dev/null +++ b/src/GPU/pair_cg_cmm_coul_msm.cpp @@ -0,0 +1,296 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + CMM coarse grained MD potentials. Coulomb with MSM version. + Contributing author: Mike Brown +------------------------------------------------------------------------- */ + +#include "pair_cg_cmm_coul_msm.h" +#include "memory.h" +#include "atom.h" +#include "force.h" +#include "kspace.h" + +#include "string.h" + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +enum {C3=0,C4=1}; + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairCGCMMCoulMSM::PairCGCMMCoulMSM(LAMMPS *lmp) : PairCMMCommon(lmp) +{ + respa_enable = 0; + single_enable = 0; +} + +/* ---------------------------------------------------------------------- */ + +PairCGCMMCoulMSM::~PairCGCMMCoulMSM() +{ + if (allocated_coul) { + memory->destroy_2d_double_array(cut_lj); + memory->destroy_2d_double_array(cut_ljsq); + memory->destroy_2d_double_array(cut_coul); + memory->destroy_2d_double_array(cut_coulsq); + allocated_coul=0; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairCGCMMCoulMSM::allocate() +{ + PairCMMCommon::allocate(); + allocated_coul = 1; + + int n = atom->ntypes; + + cut_lj = memory->create_2d_double_array(n+1,n+1,"paircg:cut_lj"); + cut_ljsq = memory->create_2d_double_array(n+1,n+1,"paircg:cut_ljsq"); + cut_coul = memory->create_2d_double_array(n+1,n+1,"paircg:cut_coul"); + cut_coulsq = memory->create_2d_double_array(n+1,n+1,"paircg:cut_coulsq"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairCGCMMCoulMSM::settings(int narg, char **arg) +{ + // strip off smoothing type and send args to parent + + if (narg < 1) error->all("Illegal pair_style command"); + + if (strcmp(arg[0],"C3") == 0) + _smooth = C3; + else if (strcmp(arg[0],"C4") == 0) + _smooth = C4; + else error->all("Illegal pair_style command"); + + PairCMMCommon::settings(narg-1,&arg[1]); +} + +/* ---------------------------------------------------------------------- */ + +void PairCGCMMCoulMSM::init_style() +{ + if (!atom->q_flag) + error->all("Pair style cg/cut/coul/msm requires atom attribute q"); + + PairCMMCommon::init_style(); + _ia=-1.0/cut_coul_global; + _ia2=_ia*_ia; + _ia3=_ia2*_ia; + + cut_respa = NULL; +} + +/* ---------------------------------------------------------------------- */ + +double PairCGCMMCoulMSM::init_one(int i, int j) +{ + double mycut = PairCMMCommon::init_one(i,j); + + return mycut; +} + +/* ---------------------------------------------------------------------- */ + +void PairCGCMMCoulMSM::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype,itable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz; + double fraction,table; + double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double grij,expm2,prefactor,t,erfc; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + double qqrd2e = force->qqrd2e; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + + if (j < nall) factor_coul = factor_lj = 1.0; + else { + factor_coul = special_coul[j/nall]; + factor_lj = special_lj[j/nall]; + j %= nall; + } + + const double delx = xtmp - x[j][0]; + const double dely = ytmp - x[j][1]; + const double delz = ztmp - x[j][2]; + const double rsq = delx*delx + dely*dely + delz*delz; + const int jtype = type[j]; + + double evdwl = 0.0; + double ecoul = 0.0; + double fpair = 0.0; + + if (rsq < cutsq[itype][jtype]) { + const double r2inv = 1.0/rsq; + const int cgt=cg_type[itype][jtype]; + + double forcelj = 0.0; + double forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + forcelj=factor_lj; + if (eflag) evdwl=factor_lj; + + if (cgt == CG_LJ12_4) { + const double r4inv=r2inv*r2inv; + forcelj *= r4inv*(lj1[itype][jtype]*r4inv*r4inv + - lj2[itype][jtype]); + if (eflag) { + evdwl *= r4inv*(lj3[itype][jtype]*r4inv*r4inv + - lj4[itype][jtype]) - offset[itype][jtype]; + } + } else if (cgt == CG_LJ9_6) { + const double r3inv = r2inv*sqrt(r2inv); + const double r6inv = r3inv*r3inv; + forcelj *= r6inv*(lj1[itype][jtype]*r3inv + - lj2[itype][jtype]); + if (eflag) { + evdwl *= r6inv*(lj3[itype][jtype]*r3inv + - lj4[itype][jtype]) - offset[itype][jtype]; + } + } else { + const double r6inv = r2inv*r2inv*r2inv; + forcelj *= r6inv*(lj1[itype][jtype]*r6inv + - lj2[itype][jtype]); + if (eflag) { + evdwl *= r6inv*(lj3[itype][jtype]*r6inv + - lj4[itype][jtype]) - offset[itype][jtype]; + } + } + } + + if (rsq < cut_coulsq_global) { + const double ir = 1.0/sqrt(rsq); + const double prefactor = qqrd2e * qtmp*q[j]; + const double r2_ia2 = rsq*_ia2; + const double r4_ia4 = r2_ia2*r2_ia2; + if (_smooth==C3) { + forcecoul = prefactor*(_ia3*(-4.375+5.25*r2_ia2-1.875*r4_ia4)- + ir/rsq); + if (eflag) + ecoul = prefactor*(ir+_ia*(2.1875-2.1875*r2_ia2+ + 1.3125*r4_ia4-0.3125*r2_ia2*r4_ia4)); + } else { + const double r6_ia6 = r2_ia2*r4_ia4; + forcecoul = prefactor*(_ia3*(-6.5625+11.8125*r2_ia2-8.4375*r4_ia4+ + 2.1875*r6_ia6)-ir/rsq); + if (eflag) + ecoul = prefactor*(ir+_ia*(2.4609375-3.28125*r2_ia2+ + 2.953125*r4_ia4-1.40625*r6_ia6+ + 0.2734375*r4_ia4*r4_ia4)); + } + if (factor_coul < 1.0) { + forcecoul -= (1.0-factor_coul)*prefactor*ir; + if (eflag) ecoul -= (1.0-factor_coul)*prefactor*ir; + } + } + fpair = forcecoul + forcelj * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairCGCMMCoulMSM::write_restart(FILE *fp) +{ + write_restart_settings(fp); + PairCMMCommon::write_restart(fp); +} + +/* ---------------------------------------------------------------------- */ + +void PairCGCMMCoulMSM::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + PairCMMCommon::read_restart(fp); +} + +/* ---------------------------------------------------------------------- */ + +double PairCGCMMCoulMSM::memory_usage() +{ + double bytes=PairCMMCommon::memory_usage(); + + int n = atom->ntypes; + + // cut_coul/cut_coulsq/cut_ljsq + bytes += (n+1)*(n+1)*sizeof(double)*4; + + return bytes; +} + +/* ---------------------------------------------------------------------- */ + +void *PairCGCMMCoulMSM::extract(char *str) +{ + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul_global; + return NULL; +} + +/* ---------------------------------------------------------------------- */ diff --git a/src/GPU/pair_cg_cmm_coul_msm.h b/src/GPU/pair_cg_cmm_coul_msm.h new file mode 100644 index 0000000000..91fd9c17a4 --- /dev/null +++ b/src/GPU/pair_cg_cmm_coul_msm.h @@ -0,0 +1,53 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(cg/cmm/coul/msm,PairCGCMMCoulMSM) + +#else + +#ifndef LMP_PAIR_CG_CMM_COUL_MSM_H +#define LMP_PAIR_CG_CMM_COUL_MSM_H + +#include "pair_cmm_common.h" + +namespace LAMMPS_NS { + +class PairCGCMMCoulMSM : public PairCMMCommon { + public: + PairCGCMMCoulMSM(class LAMMPS *); + ~PairCGCMMCoulMSM(); + + void compute(int, int); + void settings(int, char **); + void init_style(); + double init_one(int, int); + + void write_restart(FILE *); + void read_restart(FILE *); + + double memory_usage(); + + void *extract(char *str); + + protected: + void allocate(); + double _ia, _ia2, _ia3; + int _smooth; +}; + +} + +#endif +#endif diff --git a/src/GPU/pair_cg_cmm_coul_msm_gpu.cpp b/src/GPU/pair_cg_cmm_coul_msm_gpu.cpp new file mode 100644 index 0000000000..9927445053 --- /dev/null +++ b/src/GPU/pair_cg_cmm_coul_msm_gpu.cpp @@ -0,0 +1,459 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Mike Brown (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "pair_cg_cmm_coul_msm_gpu.h" +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "integrate.h" +#include "memory.h" +#include "error.h" +#include "neigh_request.h" +#include "universe.h" +#include "update.h" +#include "domain.h" +#include "string.h" +#include "kspace.h" + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) +enum {C3=0,C4=1}; + +// External functions from cuda library for atom decomposition + +bool cmmm_gpu_init(const int ntypes, double **cutsq, int **cg_type, + double **host_lj1, double **host_lj2, double **host_lj3, + double **host_lj4, double **offset, double *special_lj, + const int nlocal, const int nall, const int max_nbors, + const int maxspecial, const double cell_size, int &gpu_mode, + FILE *screen, double **host_cut_ljsq, double host_cut_coulsq, + double *host_special_coul, const double qqrd2e, + const int smooth); +void cmmm_gpu_clear(); +int * cmmm_gpu_compute_n(const int timestep, const int ago, const int inum, + const int nall, double **host_x, int *host_type, + double *boxlo, double *boxhi, int *tag, int **nspecial, + int **special, const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + const double cpu_time, bool &success, double *host_q); +void cmmm_gpu_compute(const int timestep, const int ago, const int inum, + const int nall, double **host_x, int *host_type, + int *ilist, int *numj, int **firstneigh, + const bool eflag, const bool vflag, const bool eatom, + const bool vatom, int &host_start, const double cpu_time, + bool &success, double *host_q); +double cmmm_gpu_bytes(); + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairCGCMMCoulMSMGPU::PairCGCMMCoulMSMGPU(LAMMPS *lmp) : PairCGCMMCoulMSM(lmp), + gpu_mode(GPU_PAIR) +{ + respa_enable = 0; + cpu_time = 0.0; +} + +/* ---------------------------------------------------------------------- + free all arrays +------------------------------------------------------------------------- */ + +PairCGCMMCoulMSMGPU::~PairCGCMMCoulMSMGPU() +{ + cmmm_gpu_clear(); +} + +/* ---------------------------------------------------------------------- */ + +void PairCGCMMCoulMSMGPU::compute(int eflag, int vflag) +{ + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + int nall = atom->nlocal + atom->nghost; + int inum, host_start; + + bool success = true; + + if (gpu_mode == GPU_NEIGH) { + inum = atom->nlocal; + gpulist = cmmm_gpu_compute_n(update->ntimestep, neighbor->ago, inum, nall, + atom->x, atom->type, domain->sublo, + domain->subhi, atom->tag, atom->nspecial, + atom->special, eflag, vflag, eflag_atom, + vflag_atom, host_start, cpu_time, success, + atom->q); + } else { + inum = list->inum; + cmmm_gpu_compute(update->ntimestep, neighbor->ago, inum, nall, atom->x, + atom->type, list->ilist, list->numneigh, list->firstneigh, + eflag, vflag, eflag_atom, vflag_atom, host_start, cpu_time, + success, atom->q); + } + if (!success) + error->one("Out of memory on GPGPU"); + + if (host_startpair_match("gpu",0) == NULL) + error->all("Cannot use pair hybrid with multiple GPU pair styles"); + + // Repeat cutsq calculation because done after call to init_style + double maxcut = -1.0; + double cut; + for (int i = 1; i <= atom->ntypes; i++) { + for (int j = i; j <= atom->ntypes; j++) { + if (setflag[i][j] != 0 || (setflag[i][i] != 0 && setflag[j][j] != 0)) { + cut = init_one(i,j); + cut *= cut; + if (cut > maxcut) + maxcut = cut; + cutsq[i][j] = cutsq[j][i] = cut; + } else + cutsq[i][j] = cutsq[j][i] = 0.0; + } + } + double cell_size = sqrt(maxcut) + neighbor->skin; + + int maxspecial=0; + if (atom->molecular) + maxspecial=atom->maxspecial; + bool init_ok = cmmm_gpu_init(atom->ntypes+1, cutsq, cg_type, lj1, lj2, lj3, + lj4, offset, force->special_lj, atom->nlocal, + atom->nlocal+atom->nghost, 300, maxspecial, + cell_size, gpu_mode, screen, cut_ljsq, + cut_coulsq_global, force->special_coul, + force->qqrd2e,_smooth); + if (!init_ok) + error->one("Insufficient memory on accelerator (or no fix gpu).\n"); + + if (force->newton_pair) + error->all("Cannot use newton pair with GPU cg/cmm pair style"); + + if (gpu_mode != GPU_NEIGH) { + int irequest = neighbor->request(this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + } +} + +/* ---------------------------------------------------------------------- */ + +double PairCGCMMCoulMSMGPU::memory_usage() +{ + double bytes = Pair::memory_usage(); + return bytes + cmmm_gpu_bytes(); +} + +/* ---------------------------------------------------------------------- */ + +void PairCGCMMCoulMSMGPU::cpu_compute(int start, int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype,itable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz; + double fraction,table; + double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double grij,expm2,prefactor,t,erfc; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + double qqrd2e = force->qqrd2e; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = start; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + + if (j < nall) factor_coul = factor_lj = 1.0; + else { + factor_coul = special_coul[j/nall]; + factor_lj = special_lj[j/nall]; + j %= nall; + } + + const double delx = xtmp - x[j][0]; + const double dely = ytmp - x[j][1]; + const double delz = ztmp - x[j][2]; + const double rsq = delx*delx + dely*dely + delz*delz; + const int jtype = type[j]; + + double evdwl = 0.0; + double ecoul = 0.0; + double fpair = 0.0; + + if (rsq < cutsq[itype][jtype]) { + const double r2inv = 1.0/rsq; + const int cgt=cg_type[itype][jtype]; + + double forcelj = 0.0; + double forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + forcelj=factor_lj; + if (eflag) evdwl=factor_lj; + + if (cgt == CG_LJ12_4) { + const double r4inv=r2inv*r2inv; + forcelj *= r4inv*(lj1[itype][jtype]*r4inv*r4inv + - lj2[itype][jtype]); + if (eflag) { + evdwl *= r4inv*(lj3[itype][jtype]*r4inv*r4inv + - lj4[itype][jtype]) - offset[itype][jtype]; + } + } else if (cgt == CG_LJ9_6) { + const double r3inv = r2inv*sqrt(r2inv); + const double r6inv = r3inv*r3inv; + forcelj *= r6inv*(lj1[itype][jtype]*r3inv + - lj2[itype][jtype]); + if (eflag) { + evdwl *= r6inv*(lj3[itype][jtype]*r3inv + - lj4[itype][jtype]) - offset[itype][jtype]; + } + } else { + const double r6inv = r2inv*r2inv*r2inv; + forcelj *= r6inv*(lj1[itype][jtype]*r6inv + - lj2[itype][jtype]); + if (eflag) { + evdwl *= r6inv*(lj3[itype][jtype]*r6inv + - lj4[itype][jtype]) - offset[itype][jtype]; + } + } + } + + if (rsq < cut_coulsq_global) { + const double ir = 1.0/sqrt(rsq); + const double prefactor = qqrd2e * qtmp*q[j]; + const double r2_ia2 = rsq*_ia2; + const double r4_ia4 = r2_ia2*r2_ia2; + if (_smooth==C3) { + forcecoul = prefactor*(_ia3*(-4.375+5.25*r2_ia2-1.875*r4_ia4)- + ir/rsq); + if (eflag) + ecoul = prefactor*(ir+_ia*(2.1875-2.1875*r2_ia2+ + 1.3125*r4_ia4-0.3125*r2_ia2*r4_ia4)); + } else { + const double r6_ia6 = r2_ia2*r4_ia4; + forcecoul = prefactor*(_ia3*(-6.5625+11.8125*r2_ia2-8.4375*r4_ia4+ + 2.1875*r6_ia6)-ir/rsq); + if (eflag) + ecoul = prefactor*(ir+_ia*(2.4609375-3.28125*r2_ia2+ + 2.953125*r4_ia4-1.40625*r6_ia6+ + 0.2734375*r4_ia4*r4_ia4)); + } + if (factor_coul < 1.0) { + forcecoul -= (1.0-factor_coul)*prefactor*ir; + if (eflag) ecoul -= (1.0-factor_coul)*prefactor*ir; + } + } + fpair = forcecoul + forcelj * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (evflag) ev_tally_full(i,evdwl,ecoul,fpair,delx,dely,delz); + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void PairCGCMMCoulMSMGPU::cpu_compute(int *nbors, int start, int eflag, + int vflag) +{ + int i,j,jnum,itype,jtype,itable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz; + double fraction,table; + double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double grij,expm2,prefactor,t,erfc; + double rsq; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + int stride = nlocal-start; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + double qqrd2e = force->qqrd2e; + + // loop over neighbors of my atoms + + for (i = start; i < nlocal; i++) { + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + int *nbor = nbors + i - start; + jnum = *nbor; + nbor += stride; + int *nbor_end = nbor + stride * jnum; + + for (; nbor (b) ? (a) : (b)) + +// External functions from cuda library for atom decomposition + +bool ljl_gpu_init(const int ntypes, double **cutsq, double **host_lj1, + double **host_lj2, double **host_lj3, double **host_lj4, + double **offset, double *special_lj, const int nlocal, + const int nall, const int max_nbors, const int maxspecial, + const double cell_size, int &gpu_mode, FILE *screen); +void ljl_gpu_clear(); +int * ljl_gpu_compute_n(const int timestep, const int ago, const int inum, + const int nall, double **host_x, int *host_type, + double *boxlo, double *boxhi, int *tag, int **nspecial, + int **special, const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + const double cpu_time, bool &success); +void ljl_gpu_compute(const int timestep, const int ago, const int inum, + const int nall, double **host_x, int *host_type, + int *ilist, int *numj, int **firstneigh, + const bool eflag, const bool vflag, const bool eatom, + const bool vatom, int &host_start, const double cpu_time, + bool &success); +double ljl_gpu_bytes(); + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairLJCutTGPU::PairLJCutTGPU(LAMMPS *lmp) : PairLJCut(lmp), gpu_mode(GPU_PAIR) +{ + respa_enable = 0; + cpu_time = 0.0; + + omp = new PairOMPGPU(lmp); +} + +/* ---------------------------------------------------------------------- + free all arrays +------------------------------------------------------------------------- */ + +PairLJCutTGPU::~PairLJCutTGPU() +{ + ljl_gpu_clear(); + delete omp; +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutTGPU::compute(int eflag, int vflag) +{ + if (eflag || vflag) { + ev_setup(eflag,vflag); + omp->ev_setup_thr(eflag,vflag,eflag_either,eflag_global,eflag_atom, + vflag_either,vflag_global,vflag_atom); + } else evflag = vflag_fdotr = 0; + + int nall = atom->nlocal + atom->nghost; + int inum, host_start; + + bool success = true; + + if (gpu_mode == GPU_NEIGH) { + inum = atom->nlocal; + gpulist = ljl_gpu_compute_n(update->ntimestep, neighbor->ago, inum, nall, + atom->x, atom->type, domain->sublo, + domain->subhi, atom->tag, atom->nspecial, + atom->special, eflag, vflag, eflag_atom, + vflag_atom, host_start, cpu_time, success); + } else { + inum = list->inum; + ljl_gpu_compute(update->ntimestep, neighbor->ago, inum, nall, atom->x, + atom->type, list->ilist, list->numneigh, list->firstneigh, + eflag, vflag, eflag_atom, vflag_atom, host_start, cpu_time, + success); + } + if (!success) + error->one("Out of memory on GPGPU"); + + if (host_startev_reduce_thr(*this); + cpu_time = MPI_Wtime() - cpu_time; + } +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJCutTGPU::init_style() +{ + omp->init_style(); + cut_respa = NULL; + + if (force->pair_match("gpu",0) == NULL) + error->all("Cannot use pair hybrid with multiple GPU pair styles"); + + // Repeat cutsq calculation because done after call to init_style + double maxcut = -1.0; + double cut; + for (int i = 1; i <= atom->ntypes; i++) { + for (int j = i; j <= atom->ntypes; j++) { + if (setflag[i][j] != 0 || (setflag[i][i] != 0 && setflag[j][j] != 0)) { + cut = init_one(i,j); + cut *= cut; + if (cut > maxcut) + maxcut = cut; + cutsq[i][j] = cutsq[j][i] = cut; + } else + cutsq[i][j] = cutsq[j][i] = 0.0; + } + } + double cell_size = sqrt(maxcut) + neighbor->skin; + + int maxspecial=0; + if (atom->molecular) + maxspecial=atom->maxspecial; + bool init_ok = ljl_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4, + offset, force->special_lj, atom->nlocal, + atom->nlocal+atom->nghost, 300, maxspecial, + cell_size, gpu_mode, screen); + if (!init_ok) + error->one("Insufficient memory on accelerator (or no fix gpu).\n"); + + if (force->newton_pair) + error->all("Cannot use newton pair with GPU LJ pair style"); + + if (gpu_mode != GPU_NEIGH) { + int irequest = neighbor->request(this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + } +} + +/* ---------------------------------------------------------------------- */ + +double PairLJCutTGPU::memory_usage() +{ + double bytes = Pair::memory_usage(); + return bytes + ljl_gpu_bytes()+omp->memory_usage(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutTGPU::cpu_compute(int start, int eflag, int vflag) { + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq,r2inv,r6inv,forcelj,factor_lj; + int *ilist,*jlist,*numneigh,**firstneigh; + + double **x = atom->x; + int *type = atom->type; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + double *special_lj = force->special_lj; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + int iifrom, iito, tid; + double **f = omp->loop_setup_thr_full(atom->f,iifrom,iito,tid,start,inum, + nall); + for (ii = iifrom; ii < iito; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + + if (j < nall) factor_lj = 1.0; + else { + factor_lj = special_lj[j/nall]; + j %= nall; + } + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r2inv = 1.0/rsq; + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + fpair = factor_lj*forcelj*r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + + if (eflag) { + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + evdwl *= factor_lj; + } + + #ifdef _OPENMP + if (evflag) omp->ev_tally_full_thr(i,evdwl,0.0,fpair,delx,dely,delz, + tid); + #else + if (evflag) ev_tally_full(i,evdwl,0.0,fpair,delx,dely,delz); + #endif + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJCutTGPU::cpu_compute(int *nbors, int start, int eflag, int vflag) { + int i,j,itype,jtype; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + int stride = nlocal-start; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq,r2inv,r6inv,forcelj,factor_lj; + double *special_lj = force->special_lj; + + double **x = atom->x; + int *type = atom->type; + + // loop over neighbors of my atoms + + int iifrom, iito, tid; + double **f = omp->loop_setup_thr(atom->f,iifrom,iito,tid,start,nlocal,nlocal); + for (i = iifrom; i < iito; i++) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + int *nbor = nbors + i - start; + int jnum = *nbor; + nbor += stride; + int *nbor_end = nbor + stride * jnum; + + for (; nborev_tally_full_thr(i,evdwl,0.0,fpair,delx,dely, + delz,tid); + #else + if (evflag) ev_tally_full(i,evdwl,0.0,fpair,delx,dely,delz); + #endif + } else { + if (jev_tally_thr(i,j,nlocal,0,evdwl,0.0,fpair, + delx,dely,delz,tid); + #else + if (evflag) ev_tally(i,j,nlocal,0, + evdwl,0.0,fpair,delx,dely,delz); + #endif + } + } + } + } + omp->force_reduce_thr(atom->f, nlocal, tid); +} diff --git a/src/GPU/pair_lj_cut_tgpu.h b/src/GPU/pair_lj_cut_tgpu.h new file mode 100644 index 0000000000..5a8af718a7 --- /dev/null +++ b/src/GPU/pair_lj_cut_tgpu.h @@ -0,0 +1,51 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(lj/cut/tgpu,PairLJCutTGPU) + +#else + +#ifndef LMP_PAIR_LJ_LIGHT_TGPU_H +#define LMP_PAIR_LJ_LIGHT_TGPU_H + +#include "pair_lj_cut.h" +#include "pair_omp_gpu.h" + +namespace LAMMPS_NS { + +class PairLJCutTGPU : public PairLJCut { + public: + PairLJCutTGPU(LAMMPS *lmp); + ~PairLJCutTGPU(); + void cpu_compute(int, int, int); + void cpu_compute(int *, int, int, int); + void compute(int, int); + void init_style(); + double memory_usage(); + + enum { GPU_PAIR, GPU_NEIGH }; + + private: + int gpu_mode; + double cpu_time; + int *gpulist; + + PairOMPGPU *omp; +}; + +} +#endif +#endif + diff --git a/src/GPU/pair_omp_gpu.cpp b/src/GPU/pair_omp_gpu.cpp new file mode 100644 index 0000000000..90c0dddb66 --- /dev/null +++ b/src/GPU/pair_omp_gpu.cpp @@ -0,0 +1,667 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) + Modified by Mike for use with GPU library +------------------------------------------------------------------------- */ + +#if defined(_OPENMP) + +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "pair_omp_gpu.h" +#include "memory.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairOMPGPU::PairOMPGPU(LAMMPS *lmp) : Pointers(lmp) +{ + eng_vdwl_thr = NULL; + eng_coul_thr = NULL; + virial_thr = NULL; + eatom_thr = NULL; + vatom_thr = NULL; + f_thr = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairOMPGPU::~PairOMPGPU() +{ + mem_free(); +} + +/* ---------------------------------------------------------------------- + free any allocated memory +------------------------------------------------------------------------- */ + +void PairOMPGPU::mem_free() { + memory->sfree(eng_vdwl_thr); + memory->sfree(eng_coul_thr); + memory->destroy_2d_double_array(virial_thr); + memory->destroy_2d_double_array(eatom_thr); + memory->destroy_3d_double_array(vatom_thr); + memory->destroy_2d_double_array(f_thr); + eng_vdwl_thr = NULL; + eng_coul_thr = NULL; + virial_thr = NULL; + eatom_thr = NULL; + vatom_thr = NULL; + f_thr = NULL; + _nmax = 0; +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairOMPGPU::init_style() +{ + mem_free(); + + #pragma omp parallel + { + int th_id = omp_get_thread_num(); + #pragma omp barrier + if (th_id == 0) + _nthreads = omp_get_num_threads(); + } + + // for hybrid OpenMP/MPI we need multiple copies + // of some accumulators to avoid race conditions + + eng_vdwl_thr = (double *)memory->smalloc(_nthreads*sizeof(double), + "pair:eng_vdwl_thr"); + eng_coul_thr = (double *)memory->smalloc(_nthreads*sizeof(double), + "pair:eng_coul_thr"); + virial_thr = memory->create_2d_double_array(_nthreads,6,"pair:virial_thr"); + + maxeatom_thr = maxvatom_thr = 0; +} + +/* ---------------------------------------------------------------------- + setup for energy, virial computation. additional code for multi-threading + see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) +------------------------------------------------------------------------- */ + +void PairOMPGPU::ev_setup_thr(int eflag, int vflag, int _eflag_either, + int _eflag_global, int _eflag_atom, + int _vflag_either, int _vflag_global, + int _vflag_atom) +{ + eflag_either=_eflag_either; + eflag_global=_eflag_global; + eflag_atom=_eflag_atom; + vflag_either=_vflag_either; + vflag_global=_vflag_global; + vflag_atom=_vflag_atom; + int i,n,t; + + // reallocate per-atom arrays if necessary + if (eflag_atom && atom->nmax > maxeatom_thr) { + maxeatom_thr = atom->nmax; + memory->destroy_2d_double_array(eatom_thr); + eatom_thr = memory->create_2d_double_array(_nthreads, + maxeatom_thr,"pair:eatom_thr"); + } + if (vflag_atom && atom->nmax > maxvatom_thr) { + maxvatom_thr = atom->nmax; + memory->destroy_3d_double_array(vatom_thr); + vatom_thr = memory->create_3d_double_array(_nthreads, + maxvatom_thr,6,"pair:vatom_thr"); + } + + // zero per thread accumulators + // use force->newton instead of newton_pair + // b/c some bonds/dihedrals call pair::ev_tally with pairwise info + const int ntotal = (force->newton) ? + (atom->nlocal + atom->nghost) : atom->nlocal; + for (t = 0; t < _nthreads; ++t) { + if (eflag_global) eng_vdwl_thr[t] = eng_coul_thr[t] = 0.0; + if (vflag_global) for (i = 0; i < 6; ++i) virial_thr[t][i] = 0.0; + if (eflag_atom) { + for (i = 0; i < ntotal; ++i) eatom_thr[t][i] = 0.0; + } + if (vflag_atom) { + for (i = 0; i < ntotal; ++i) { + vatom_thr[t][i][0] = 0.0; + vatom_thr[t][i][1] = 0.0; + vatom_thr[t][i][2] = 0.0; + vatom_thr[t][i][3] = 0.0; + vatom_thr[t][i][4] = 0.0; + vatom_thr[t][i][5] = 0.0; + } + } + } +} + +/* ---------------------------------------------------------------------- + tally eng_vdwl and virial into per thread global and per-atom accumulators + need i < nlocal test since called by bond_quartic and dihedral_charmm +------------------------------------------------------------------------- */ + +void PairOMPGPU::ev_tally_thr(int i, int j, int nlocal, int newton_pair, + double evdwl, double ecoul, double fpair, + double delx, double dely, double delz, int tid) +{ + double evdwlhalf,ecoulhalf,epairhalf,v[6]; + + if (eflag_either) { + if (eflag_global) { + if (newton_pair) { + eng_vdwl_thr[tid] += evdwl; + eng_coul_thr[tid] += ecoul; + } else { + evdwlhalf = 0.5*evdwl; + ecoulhalf = 0.5*ecoul; + if (i < nlocal) { + eng_vdwl_thr[tid] += evdwlhalf; + eng_coul_thr[tid] += ecoulhalf; + } + if (j < nlocal) { + eng_vdwl_thr[tid] += evdwlhalf; + eng_coul_thr[tid] += ecoulhalf; + } + } + } + if (eflag_atom) { + epairhalf = 0.5 * (evdwl + ecoul); + if (newton_pair || i < nlocal) eatom_thr[tid][i] += epairhalf; + if (newton_pair || j < nlocal) eatom_thr[tid][j] += epairhalf; + } + } + + if (vflag_either) { + v[0] = delx*delx*fpair; + v[1] = dely*dely*fpair; + v[2] = delz*delz*fpair; + v[3] = delx*dely*fpair; + v[4] = delx*delz*fpair; + v[5] = dely*delz*fpair; + + if (vflag_global) { + if (newton_pair) { + virial_thr[tid][0] += v[0]; + virial_thr[tid][1] += v[1]; + virial_thr[tid][2] += v[2]; + virial_thr[tid][3] += v[3]; + virial_thr[tid][4] += v[4]; + virial_thr[tid][5] += v[5]; + } else { + if (i < nlocal) { + virial_thr[tid][0] += 0.5*v[0]; + virial_thr[tid][1] += 0.5*v[1]; + virial_thr[tid][2] += 0.5*v[2]; + virial_thr[tid][3] += 0.5*v[3]; + virial_thr[tid][4] += 0.5*v[4]; + virial_thr[tid][5] += 0.5*v[5]; + } + if (j < nlocal) { + virial_thr[tid][0] += 0.5*v[0]; + virial_thr[tid][1] += 0.5*v[1]; + virial_thr[tid][2] += 0.5*v[2]; + virial_thr[tid][3] += 0.5*v[3]; + virial_thr[tid][4] += 0.5*v[4]; + virial_thr[tid][5] += 0.5*v[5]; + } + } + } + + if (vflag_atom) { + if (newton_pair || i < nlocal) { + vatom_thr[tid][i][0] += 0.5*v[0]; + vatom_thr[tid][i][1] += 0.5*v[1]; + vatom_thr[tid][i][2] += 0.5*v[2]; + vatom_thr[tid][i][3] += 0.5*v[3]; + vatom_thr[tid][i][4] += 0.5*v[4]; + vatom_thr[tid][i][5] += 0.5*v[5]; + } + if (newton_pair || j < nlocal) { + vatom_thr[tid][j][0] += 0.5*v[0]; + vatom_thr[tid][j][1] += 0.5*v[1]; + vatom_thr[tid][j][2] += 0.5*v[2]; + vatom_thr[tid][j][3] += 0.5*v[3]; + vatom_thr[tid][j][4] += 0.5*v[4]; + vatom_thr[tid][j][5] += 0.5*v[5]; + } + } + } +} + +/* ---------------------------------------------------------------------- + tally eng_vdwl and virial into per thread global and per-atom accumulators + need i < nlocal test since called by bond_quartic and dihedral_charmm +------------------------------------------------------------------------- */ + +void PairOMPGPU::ev_tally_full_thr(int i, double evdwl, + double ecoul, double fpair, double delx, + double dely, double delz, int tid) +{ + double evdwlhalf,ecoulhalf,epairhalf,v[6]; + + if (eflag_either) { + if (eflag_global) { + evdwlhalf = 0.5*evdwl; + ecoulhalf = 0.5*ecoul; + eng_vdwl_thr[tid] += evdwlhalf; + eng_coul_thr[tid] += ecoulhalf; + } + if (eflag_atom) { + epairhalf = 0.5 * (evdwl + ecoul); + eatom_thr[tid][i] += epairhalf; + } + } + + if (vflag_either) { + v[0] = delx*delx*fpair; + v[1] = dely*dely*fpair; + v[2] = delz*delz*fpair; + v[3] = delx*dely*fpair; + v[4] = delx*delz*fpair; + v[5] = dely*delz*fpair; + + if (vflag_global) { + virial_thr[tid][0] += 0.5*v[0]; + virial_thr[tid][1] += 0.5*v[1]; + virial_thr[tid][2] += 0.5*v[2]; + virial_thr[tid][3] += 0.5*v[3]; + virial_thr[tid][4] += 0.5*v[4]; + virial_thr[tid][5] += 0.5*v[5]; + } + + if (vflag_atom) { + vatom_thr[tid][i][0] += 0.5*v[0]; + vatom_thr[tid][i][1] += 0.5*v[1]; + vatom_thr[tid][i][2] += 0.5*v[2]; + vatom_thr[tid][i][3] += 0.5*v[3]; + vatom_thr[tid][i][4] += 0.5*v[4]; + vatom_thr[tid][i][5] += 0.5*v[5]; + } + } +} + +/* ---------------------------------------------------------------------- + tally eng_vdwl and virial into global and per-atom accumulators + for virial, have delx,dely,delz and fx,fy,fz +------------------------------------------------------------------------- */ + +void PairOMPGPU::ev_tally_xyz_thr(int i, int j, int nlocal, int newton_pair, + double evdwl, double ecoul, + double fx, double fy, double fz, + double delx, double dely, double delz, int tid) +{ + double evdwlhalf,ecoulhalf,epairhalf,v[6]; + + if (eflag_either) { + if (eflag_global) { + if (newton_pair) { + eng_vdwl_thr[tid] += evdwl; + eng_coul_thr[tid] += ecoul; + } else { + evdwlhalf = 0.5*evdwl; + ecoulhalf = 0.5*ecoul; + if (i < nlocal) { + eng_vdwl_thr[tid] += evdwlhalf; + eng_coul_thr[tid] += ecoulhalf; + } + if (j < nlocal) { + eng_vdwl_thr[tid] += evdwlhalf; + eng_coul_thr[tid] += ecoulhalf; + } + } + } + if (eflag_atom) { + epairhalf = 0.5 * (evdwl + ecoul); + if (newton_pair || i < nlocal) eatom_thr[tid][i] += epairhalf; + if (newton_pair || j < nlocal) eatom_thr[tid][j] += epairhalf; + } + } + + if (vflag_either) { + v[0] = delx*fx; + v[1] = dely*fy; + v[2] = delz*fz; + v[3] = delx*fy; + v[4] = delx*fz; + v[5] = dely*fz; + + if (vflag_global) { + if (newton_pair) { + virial_thr[tid][0] += v[0]; + virial_thr[tid][1] += v[1]; + virial_thr[tid][2] += v[2]; + virial_thr[tid][3] += v[3]; + virial_thr[tid][4] += v[4]; + virial_thr[tid][5] += v[5]; + } else { + if (i < nlocal) { + virial_thr[tid][0] += 0.5*v[0]; + virial_thr[tid][1] += 0.5*v[1]; + virial_thr[tid][2] += 0.5*v[2]; + virial_thr[tid][3] += 0.5*v[3]; + virial_thr[tid][4] += 0.5*v[4]; + virial_thr[tid][5] += 0.5*v[5]; + } + if (j < nlocal) { + virial_thr[tid][0] += 0.5*v[0]; + virial_thr[tid][1] += 0.5*v[1]; + virial_thr[tid][2] += 0.5*v[2]; + virial_thr[tid][3] += 0.5*v[3]; + virial_thr[tid][4] += 0.5*v[4]; + virial_thr[tid][5] += 0.5*v[5]; + } + } + } + + if (vflag_atom) { + if (newton_pair || i < nlocal) { + vatom_thr[tid][i][0] += 0.5*v[0]; + vatom_thr[tid][i][1] += 0.5*v[1]; + vatom_thr[tid][i][2] += 0.5*v[2]; + vatom_thr[tid][i][3] += 0.5*v[3]; + vatom_thr[tid][i][4] += 0.5*v[4]; + vatom_thr[tid][i][5] += 0.5*v[5]; + } + if (newton_pair || j < nlocal) { + vatom_thr[tid][j][0] += 0.5*v[0]; + vatom_thr[tid][j][1] += 0.5*v[1]; + vatom_thr[tid][j][2] += 0.5*v[2]; + vatom_thr[tid][j][3] += 0.5*v[3]; + vatom_thr[tid][j][4] += 0.5*v[4]; + vatom_thr[tid][j][5] += 0.5*v[5]; + } + } + } +} + +/* ---------------------------------------------------------------------- + tally eng_vdwl and virial into global and per-atom accumulators + called by SW potential, newton_pair is always on + virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk + ------------------------------------------------------------------------- */ + +void PairOMPGPU::ev_tally3_thr(int i, int j, int k, double evdwl, double ecoul, + double *fj, double *fk, double *drji, double *drki, int tid, + double THIRD) +{ + double epairthird,v[6]; + + if (eflag_either) { + if (eflag_global) { + eng_vdwl_thr[tid] += evdwl; + eng_coul_thr[tid] += ecoul; + } + if (eflag_atom) { + epairthird = THIRD * (evdwl + ecoul); + eatom_thr[tid][i] += epairthird; + eatom_thr[tid][j] += epairthird; + eatom_thr[tid][k] += epairthird; + } + } + + if (vflag_atom) { + v[0] = THIRD * (drji[0]*fj[0] + drki[0]*fk[0]); + v[1] = THIRD * (drji[1]*fj[1] + drki[1]*fk[1]); + v[2] = THIRD * (drji[2]*fj[2] + drki[2]*fk[2]); + v[3] = THIRD * (drji[0]*fj[1] + drki[0]*fk[1]); + v[4] = THIRD * (drji[0]*fj[2] + drki[0]*fk[2]); + v[5] = THIRD * (drji[1]*fj[2] + drki[1]*fk[2]); + + vatom_thr[tid][i][0] += v[0]; vatom_thr[tid][i][1] += v[1]; + vatom_thr[tid][i][2] += v[2]; vatom_thr[tid][i][3] += v[3]; + vatom_thr[tid][i][4] += v[4]; vatom_thr[tid][i][5] += v[5]; + vatom_thr[tid][j][0] += v[0]; vatom_thr[tid][j][1] += v[1]; + vatom_thr[tid][j][2] += v[2]; vatom_thr[tid][j][3] += v[3]; + vatom_thr[tid][j][4] += v[4]; vatom_thr[tid][j][5] += v[5]; + vatom_thr[tid][k][0] += v[0]; vatom_thr[tid][k][1] += v[1]; + vatom_thr[tid][k][2] += v[2]; vatom_thr[tid][k][3] += v[3]; + vatom_thr[tid][k][4] += v[4]; vatom_thr[tid][k][5] += v[5]; + } +} + +/* ---------------------------------------------------------------------- + tally eng_vdwl and virial into global and per-atom accumulators + called by AIREBO potential, newton_pair is always on + ------------------------------------------------------------------------- */ + +void PairOMPGPU::ev_tally4_thr(int i, int j, int k, int m, double evdwl, + double *fi, double *fj, double *fk, + double *drim, double *drjm, double *drkm, int tid) +{ + double epairfourth,v[6]; + + if (eflag_either) { + if (eflag_global) eng_vdwl_thr[tid] += evdwl; + if (eflag_atom) { + epairfourth = 0.25 * evdwl; + eatom_thr[tid][i] += epairfourth; + eatom_thr[tid][j] += epairfourth; + eatom_thr[tid][k] += epairfourth; + eatom_thr[tid][m] += epairfourth; + } + } + + if (vflag_atom) { + v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); + v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); + v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); + v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); + v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); + v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); + + vatom_thr[tid][i][0] += v[0]; vatom_thr[tid][i][1] += v[1]; + vatom_thr[tid][i][2] += v[2]; vatom_thr[tid][i][3] += v[3]; + vatom_thr[tid][i][4] += v[4]; vatom_thr[tid][i][5] += v[5]; + vatom_thr[tid][j][0] += v[0]; vatom_thr[tid][j][1] += v[1]; + vatom_thr[tid][j][2] += v[2]; vatom_thr[tid][j][3] += v[3]; + vatom_thr[tid][j][4] += v[4]; vatom_thr[tid][j][5] += v[5]; + vatom_thr[tid][k][0] += v[0]; vatom_thr[tid][k][1] += v[1]; + vatom_thr[tid][k][2] += v[2]; vatom_thr[tid][k][3] += v[3]; + vatom_thr[tid][k][4] += v[4]; vatom_thr[tid][k][5] += v[5]; + vatom_thr[tid][m][0] += v[0]; vatom_thr[tid][m][1] += v[1]; + vatom_thr[tid][m][2] += v[2]; vatom_thr[tid][m][3] += v[3]; + vatom_thr[tid][m][4] += v[4]; vatom_thr[tid][m][5] += v[5]; + } +} + +/* ---------------------------------------------------------------------- + tally ecoul and virial into each of n atoms in list + called by TIP4P potential, newton_pair is always on + changes v values by dividing by n + ------------------------------------------------------------------------- */ + +void PairOMPGPU::ev_tally_list_thr(int n, int *list, double ecoul, double *v, int tid) +{ + int i,j; + + if (eflag_either) { + if (eflag_global) eng_coul_thr[tid] += ecoul; + if (eflag_atom) { + double epairatom = ecoul/n; + for (i = 0; i < n; i++) eatom_thr[tid][list[i]] += epairatom; + } + } + + if (vflag_either) { + if (vflag_global) { + virial_thr[tid][0] += v[0]; + virial_thr[tid][1] += v[1]; + virial_thr[tid][2] += v[2]; + virial_thr[tid][3] += v[3]; + virial_thr[tid][4] += v[4]; + virial_thr[tid][5] += v[5]; + } + + if (vflag_atom) { + v[0] /= n; + v[1] /= n; + v[2] /= n; + v[3] /= n; + v[4] /= n; + v[5] /= n; + for (i = 0; i < n; i++) { + j = list[i]; + vatom_thr[tid][j][0] += v[0]; + vatom_thr[tid][j][1] += v[1]; + vatom_thr[tid][j][2] += v[2]; + vatom_thr[tid][j][3] += v[3]; + vatom_thr[tid][j][4] += v[4]; + vatom_thr[tid][j][5] += v[5]; + } + } + } +} + +/* ---------------------------------------------------------------------- + tally virial into per-atom accumulators + called by AIREBO potential, newton_pair is always on + fpair is magnitude of force on atom I +------------------------------------------------------------------------- */ + +void PairOMPGPU::v_tally2_thr(int i, int j, double fpair, double *drij, int tid) +{ + double v[6]; + + v[0] = 0.5 * drij[0]*drij[0]*fpair; + v[1] = 0.5 * drij[1]*drij[1]*fpair; + v[2] = 0.5 * drij[2]*drij[2]*fpair; + v[3] = 0.5 * drij[0]*drij[1]*fpair; + v[4] = 0.5 * drij[0]*drij[2]*fpair; + v[5] = 0.5 * drij[1]*drij[2]*fpair; + + vatom_thr[tid][i][0] += v[0]; vatom_thr[tid][i][1] += v[1]; + vatom_thr[tid][i][2] += v[2]; vatom_thr[tid][i][3] += v[3]; + vatom_thr[tid][i][4] += v[4]; vatom_thr[tid][i][5] += v[5]; + vatom_thr[tid][j][0] += v[0]; vatom_thr[tid][j][1] += v[1]; + vatom_thr[tid][j][2] += v[2]; vatom_thr[tid][j][3] += v[3]; + vatom_thr[tid][j][4] += v[4]; vatom_thr[tid][j][5] += v[5]; +} + +/* ---------------------------------------------------------------------- + tally virial into per-atom accumulators + called by AIREBO and Tersoff potential, newton_pair is always on +------------------------------------------------------------------------- */ + +void PairOMPGPU::v_tally3_thr(int i, int j, int k, double *fi, double *fj, + double *drik, double *drjk, int tid, + double THIRD) +{ + double v[6]; + + v[0] = THIRD * (drik[0]*fi[0] + drjk[0]*fj[0]); + v[1] = THIRD * (drik[1]*fi[1] + drjk[1]*fj[1]); + v[2] = THIRD * (drik[2]*fi[2] + drjk[2]*fj[2]); + v[3] = THIRD * (drik[0]*fi[1] + drjk[0]*fj[1]); + v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]); + v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]); + + vatom_thr[tid][i][0] += v[0]; vatom_thr[tid][i][1] += v[1]; + vatom_thr[tid][i][2] += v[2]; vatom_thr[tid][i][3] += v[3]; + vatom_thr[tid][i][4] += v[4]; vatom_thr[tid][i][5] += v[5]; + vatom_thr[tid][j][0] += v[0]; vatom_thr[tid][j][1] += v[1]; + vatom_thr[tid][j][2] += v[2]; vatom_thr[tid][j][3] += v[3]; + vatom_thr[tid][j][4] += v[4]; vatom_thr[tid][j][5] += v[5]; + vatom_thr[tid][k][0] += v[0]; vatom_thr[tid][k][1] += v[1]; + vatom_thr[tid][k][2] += v[2]; vatom_thr[tid][k][3] += v[3]; + vatom_thr[tid][k][4] += v[4]; vatom_thr[tid][k][5] += v[5]; +} + +/* ---------------------------------------------------------------------- + tally virial into per-atom accumulators + called by AIREBO potential, newton_pair is always on +------------------------------------------------------------------------- */ + +void PairOMPGPU::v_tally4_thr(int i, int j, int k, int m, + double *fi, double *fj, double *fk, + double *drim, double *drjm, double *drkm, int tid) +{ + double v[6]; + + v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); + v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); + v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); + v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); + v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); + v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); + + vatom_thr[tid][i][0] += v[0]; vatom_thr[tid][i][1] += v[1]; + vatom_thr[tid][i][2] += v[2]; vatom_thr[tid][i][3] += v[3]; + vatom_thr[tid][i][4] += v[4]; vatom_thr[tid][i][5] += v[5]; + vatom_thr[tid][j][0] += v[0]; vatom_thr[tid][j][1] += v[1]; + vatom_thr[tid][j][2] += v[2]; vatom_thr[tid][j][3] += v[3]; + vatom_thr[tid][j][4] += v[4]; vatom_thr[tid][j][5] += v[5]; + vatom_thr[tid][k][0] += v[0]; vatom_thr[tid][k][1] += v[1]; + vatom_thr[tid][k][2] += v[2]; vatom_thr[tid][k][3] += v[3]; + vatom_thr[tid][k][4] += v[4]; vatom_thr[tid][k][5] += v[5]; + vatom_thr[tid][m][0] += v[0]; vatom_thr[tid][m][1] += v[1]; + vatom_thr[tid][m][2] += v[2]; vatom_thr[tid][m][3] += v[3]; + vatom_thr[tid][m][4] += v[4]; vatom_thr[tid][m][5] += v[5]; +} + +/* ---------------------------------------------------------------------- + reduce the per thread accumulated E/V data into the canonical accumulators. +------------------------------------------------------------------------- */ +void PairOMPGPU::ev_reduce_thr(Pair &p) +{ + const int ntotal = (force->newton) ? + (atom->nlocal + atom->nghost) : atom->nlocal; + + for (int n = 0; n < _nthreads; ++n) { + p.eng_vdwl += eng_vdwl_thr[n]; + p.eng_coul += eng_coul_thr[n]; + if (vflag_either) { + p.virial[0] += virial_thr[n][0]; + p.virial[1] += virial_thr[n][1]; + p.virial[2] += virial_thr[n][2]; + p.virial[3] += virial_thr[n][3]; + p.virial[4] += virial_thr[n][4]; + p.virial[5] += virial_thr[n][5]; + if (vflag_atom) { + for (int i = 0; i < ntotal; ++i) { + p.vatom[i][0] += vatom_thr[n][i][0]; + p.vatom[i][1] += vatom_thr[n][i][1]; + p.vatom[i][2] += vatom_thr[n][i][2]; + p.vatom[i][3] += vatom_thr[n][i][3]; + p.vatom[i][4] += vatom_thr[n][i][4]; + p.vatom[i][5] += vatom_thr[n][i][5]; + } + } + } + if (eflag_atom) { + for (int i = 0; i < ntotal; ++i) { + p.eatom[i] += eatom_thr[n][i]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +double PairOMPGPU::memory_usage() +{ + double bytes = 0.0; + + bytes += _nthreads * (2 + 7) * sizeof(double); + bytes += _nthreads * maxeatom_thr * sizeof(double); + bytes += _nthreads * maxvatom_thr * 6 * sizeof(double); + + if (f_thr != NULL) + bytes += _nthreads * _nmax * sizeof(double); + return bytes; +} + +#endif diff --git a/src/GPU/pair_omp_gpu.h b/src/GPU/pair_omp_gpu.h new file mode 100644 index 0000000000..df968a87a5 --- /dev/null +++ b/src/GPU/pair_omp_gpu.h @@ -0,0 +1,255 @@ +/* -*- c++ -*- ------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) + Modified by Mike for use with GPU library +------------------------------------------------------------------------- */ + +#ifndef LMP_PAIR_OMP_GPU_H +#define LMP_PAIR_OMP_GPU_H + +#include "pair.h" +#include + +#if defined(_OPENMP) +#include "atom.h" +#include "memory.h" +#include +#endif + +namespace LAMMPS_NS { + +#if defined(_OPENMP) + +class PairOMPGPU : protected Pointers { + + protected: + double *eng_vdwl_thr; // per thread accumulated vdw energy + double *eng_coul_thr; // per thread accumulated coulomb energies + double **virial_thr; // per thread virial + double **eatom_thr; // per thread per atom energy + double ***vatom_thr; // per thread per atom virial + double **f_thr; // per thread force (for thread id>0) + + int maxeatom_thr, maxvatom_thr; + int _nthreads, _nmax; + + int eflag_either,eflag_global,eflag_atom; + int vflag_either,vflag_global,vflag_atom; + + public: + PairOMPGPU(LAMMPS *); + ~PairOMPGPU(); + + void mem_free(); + void init_style(); + double memory_usage(); + + // set loop range for, thread id, and force array offset for threaded runs. + double **loop_setup_thr(double **f, int &ifrom, int &ito, int &tid, + const int start, const int inum, const int nall) + { + if (_nthreads > 1) { + tid = omp_get_thread_num(); + if (nall > _nmax) { +#pragma omp master + f_thr = memory->grow_2d_double_array(f_thr,atom->nmax*(_nthreads-1), + 3,"pair:f_thr"); +#pragma omp barrier + if (tid == 0) + _nmax = atom->nmax; + } + + // each thread works on a fixed chunk of atoms. + const int idelta = 1 + (inum - start)/_nthreads; + ifrom = tid * idelta + start; + ito = ifrom + idelta; + if (ito > inum) + ito = inum; + + // zero per thread force array + // keep thread memory access same as for force accumulation + if (tid > 0) { + double **f_zero = f_thr + nall * (tid - 1); + for (int i = 0; i < nall; i++) { + f_zero[i][0] = 0.0; + f_zero[i][1] = 0.0; + f_zero[i][2] = 0.0; + } + return f_zero; + } + return f; + } else { + tid = 0; + ifrom = start; + ito = inum; + return f; + } + }; + + // set loop range for, thread id, and force array offset for threaded runs. + double **loop_setup_thr(double **f, int &ifrom, int &ito, int &tid, + const int inum, const int nall) + { return loop_setup_thr(f,ifrom,ito,tid,0,inum,nall); } + + // loop setup with full neighbor list and nonzero starting index + double **loop_setup_thr_full(double **f, int &ifrom, int &ito, int &tid, + const int start, const int inum, const int nall) + { + if (_nthreads > 1) { + tid = omp_get_thread_num(); + + // each thread works on a fixed chunk of atoms. + const int idelta = 1 + (inum - start) / _nthreads; + ifrom = start + tid * idelta; + ito = ifrom + idelta; + if (ito > inum) + ito = inum; + + return f; + + } else { + tid = 0; + ifrom = start; + ito = inum; + return f; + } + }; + + // threading adapted versions of the ev_tally infrastructure. + void ev_setup_thr(int, int, int, int, int, int, int, int); + void ev_reduce_thr(Pair &); + void ev_tally_thr(int, int, int, int, double, double, double, + double, double, double, int); + void ev_tally_full_thr(int, double, double, double, + double, double, double, int); + void ev_tally_xyz_thr(int, int, int, int, double, double, + double, double, double, double, double, double, int); + void ev_tally3_thr(int, int, int, double, double, + double *, double *, double *, double *, int, double); + void ev_tally4_thr(int, int, int, int, double, + double *, double *, double *, double *, double *, + double *, int); + void ev_tally_list_thr(int, int *, double, double *, int); + void v_tally2_thr(int, int, double, double *, int); + void v_tally3_thr(int, int, int, double *, double *, double *, double *, int, + double); + void v_tally4_thr(int, int, int, int, double *, double *, double *, + double *, double *, double *, int); + + // reduce per thread forces into the first part of the force + // array that is used for the non-threaded parts and reset + // the temporary storage to 0.0. + // this is in the header to be inlined. + // need to post a barrier to wait until all threads are done + // with computing forces. the reduction can be threaded as well. + void force_reduce_thr(double **fall, const int nall, const int tid) + { + // NOOP in non-threaded execution. + if (_nthreads == 1) return; +#pragma omp barrier + { + double **f; + const int idelta = 1 + nall/_nthreads; + const int ifrom = tid*idelta; + const int ito = ((ifrom + idelta) > nall) ? nall : (ifrom + idelta); + for (int n = 1; n < _nthreads; ++n) { + const int toffs = (n-1)*nall; + f = f_thr + toffs; + for (int m = ifrom; m < ito; ++m) { + fall[m][0] += f[m][0]; + f[m][0] = 0.0; + fall[m][1] += f[m][1]; + f[m][1] = 0.0; + fall[m][2] += f[m][2]; + f[m][2] = 0.0; + } + } + } + }; + + // reduce per thread density into the first part of the rho + // array that is used for the non-threaded parts. for use with EAM. + // this is in the header to be inlined. + // we need to post a barrier to wait until all threads are done. + // the reduction can be threaded as well. + void rho_reduce_thr(double *rho, const int nmax, const int nrange, + const int tid) + { + // NOOP in non-threaded execution. + if (_nthreads == 1) return; +#pragma omp barrier + { + double *rho_thr; + const int idelta = 1 + nrange/_nthreads; + const int ifrom = tid*idelta; + const int ito = ((ifrom + idelta) > nrange) ? nrange : (ifrom + idelta); + for (int n = 1; n < _nthreads; ++n) { + const int toffs = n*nmax; + rho_thr = rho + toffs; + for (int m = ifrom; m < ito; ++m) + rho[m] += rho_thr[m]; + } + } + }; + +}; + +#else + +class PairOMPGPU { + + public: + inline PairOMPGPU(LAMMPS *) {}; + virtual inline ~PairOMPGPU() {}; + + inline void init_style() {} + inline double memory_usage() { return 0.0; } + + inline void ev_setup_thr(int, int, int, int, int, int, int, int) {} + inline void ev_reduce_thr(Pair &) {} + inline double **loop_setup_thr(double **f, int &ifrom, int &ito, int &tid, + const int start, const int inum, + const int nall) + { + ifrom = start; + ito = inum; + return f; + }; + + inline double **loop_setup_thr(double **f, int &ifrom, int &ito, int &tid, + const int inum, const int nall) + { return loop_setup_thr(f,ifrom,ito,tid,0,inum,nall); } + + // loop setup with full neighbor list and nonzero starting index + double **loop_setup_thr_full(double **f, int &ifrom, int &ito, int &tid, + const int start, const int inum, const int nall) + { + ifrom = start; + ito = inum; + return f; + }; + + inline void force_reduce_thr(double **fall, const int nall, const int tid) {}; + inline void rho_reduce_thr(double *rho, const int nmax, const int nrange, + const int tid) {}; + +}; + +#endif + +} + +#endif + diff --git a/src/GPU/pppm_gpu.cpp b/src/GPU/pppm_gpu.cpp new file mode 100644 index 0000000000..4cf9118e0c --- /dev/null +++ b/src/GPU/pppm_gpu.cpp @@ -0,0 +1,163 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Roy Pollock (LLNL), Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "string.h" +#include "stdio.h" +#include "stdlib.h" +#include "math.h" +#include "pppm_gpu.h" +#include "atom.h" +#include "comm.h" +#include "neighbor.h" +#include "force.h" +#include "pair.h" +#include "bond.h" +#include "angle.h" +#include "domain.h" +#include "fft3d_wrap.h" +#include "remap_wrap.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MAXORDER 7 +#define OFFSET 16384 +#define SMALL 0.00001 +#define LARGE 10000.0 +#define EPS_HOC 1.0e-7 + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PPPM_GPU::PPPM_GPU(LAMMPS *lmp, int narg, char **arg) : PPPM(lmp, narg, arg) +{ +} + +/* ---------------------------------------------------------------------- + free all memory +------------------------------------------------------------------------- */ + +PPPM_GPU::~PPPM_GPU() +{ +} + +/* ---------------------------------------------------------------------- + called once before run +------------------------------------------------------------------------- */ + +void PPPM_GPU::init() +{ + PPPM::init(); +} + +/* ---------------------------------------------------------------------- + adjust PPPM_GPU coeffs, called initially and whenever volume has changed +------------------------------------------------------------------------- */ + +void PPPM_GPU::setup() +{ + PPPM::setup(); +} + +/* ---------------------------------------------------------------------- + compute the PPPM_GPU long-range force, energy, virial +------------------------------------------------------------------------- */ + +void PPPM_GPU::compute(int eflag, int vflag) +{ + int i; + + // convert atoms from box to lamda coords + + if (triclinic == 0) boxlo = domain->boxlo; + else { + boxlo = domain->boxlo_lamda; + domain->x2lamda(atom->nlocal); + } + + // extend size of per-atom arrays if necessary + + if (atom->nlocal > nmax) { + memory->destroy_2d_int_array(part2grid); + nmax = atom->nmax; + part2grid = memory->create_2d_int_array(nmax,3,"pppm:part2grid"); + } + + energy = 0.0; + if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; + + // find grid points for all my particles + // map my particle charge onto my local 3d density grid + + particle_map(); + make_rho(); + + // all procs communicate density values from their ghost cells + // to fully sum contribution in their 3d bricks + // remap from 3d decomposition to FFT decomposition + + brick2fft(); + + // compute potential gradient on my FFT grid and + // portion of e_long on this proc's FFT grid + // return gradients (electric fields) in 3d brick decomposition + + poisson(eflag,vflag); + + // all procs communicate E-field values to fill ghost cells + // surrounding their 3d bricks + + fillbrick(); + + // calculate the force on my particles + + fieldforce(); + + // sum energy across procs and add in volume-dependent term + + if (eflag) { + double energy_all; + MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); + energy = energy_all; + + energy *= 0.5*volume; + energy -= g_ewald*qsqsum/1.772453851 + + 0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume); + energy *= qqrd2e*scale; + } + + // sum virial across procs + + if (vflag) { + double virial_all[6]; + MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world); + for (i = 0; i < 6; i++) virial[i] = 0.5*qqrd2e*scale*volume*virial_all[i]; + } + + // 2d slab correction + + if (slabflag) slabcorr(eflag); + + // convert atoms back from lamda to box coords + + if (triclinic) domain->lamda2x(atom->nlocal); +} + diff --git a/src/GPU/pppm_gpu.h b/src/GPU/pppm_gpu.h new file mode 100644 index 0000000000..4e38c1274d --- /dev/null +++ b/src/GPU/pppm_gpu.h @@ -0,0 +1,41 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef KSPACE_CLASS + +KSpaceStyle(pppm/gpu,PPPM_GPU) + +#else + +#ifndef LMP_PPPM_GPU_H +#define LMP_PPPM_GPU_H + +#include "pppm.h" + +namespace LAMMPS_NS { + +class PPPM_GPU : public PPPM { + public: + PPPM_GPU(class LAMMPS *, int, char **); + ~PPPM_GPU(); + void init(); + void setup(); + void compute(int, int); + + protected: +}; + +} + +#endif +#endif