From 838a1938bb2b3073eae4e45ab27f037caa8beb7a Mon Sep 17 00:00:00 2001
From: Trung Nguyen <ndactrung@gmail.com>
Date: Sun, 20 May 2018 17:03:12 -0500
Subject: [PATCH] Added pair styles lj/expand/coul/long to USER-MISC and its
gpu version to GPU package, added lj/cut/dipole/long/gpu; and added an
example Makfile to build GPU package for multi-arch
---
lib/gpu/Makefile.linux_multi | 53 ++
lib/gpu/Nvidia.makefile | 32 +-
lib/gpu/Nvidia.makefile_multi | 854 ++++++++++++++++++
lib/gpu/lal_dipole_long_lj.cpp | 173 ++++
lib/gpu/lal_dipole_long_lj.cu | 640 +++++++++++++
lib/gpu/lal_dipole_long_lj.h | 85 ++
lib/gpu/lal_dipole_long_lj_ext.cpp | 129 +++
lib/gpu/lal_lj_expand_coul_long.cpp | 183 ++++
lib/gpu/lal_lj_expand_coul_long.cu | 265 ++++++
lib/gpu/lal_lj_expand_coul_long.h | 88 ++
lib/gpu/lal_lj_expand_coul_long_ext.cpp | 152 ++++
src/GPU/Install.sh | 2 +
src/GPU/pair_lj_cut_dipole_long_gpu.cpp | 434 +++++++++
src/GPU/pair_lj_cut_dipole_long_gpu.h | 66 ++
src/GPU/pair_lj_expand_coul_long_gpu.cpp | 325 +++++++
src/GPU/pair_lj_expand_coul_long_gpu.h | 67 ++
src/USER-MISC/pair_lj_expand_coul_long.cpp | 990 +++++++++++++++++++++
src/USER-MISC/pair_lj_expand_coul_long.h | 93 ++
18 files changed, 4629 insertions(+), 2 deletions(-)
create mode 100644 lib/gpu/Makefile.linux_multi
create mode 100644 lib/gpu/Nvidia.makefile_multi
create mode 100644 lib/gpu/lal_dipole_long_lj.cpp
create mode 100644 lib/gpu/lal_dipole_long_lj.cu
create mode 100644 lib/gpu/lal_dipole_long_lj.h
create mode 100644 lib/gpu/lal_dipole_long_lj_ext.cpp
create mode 100644 lib/gpu/lal_lj_expand_coul_long.cpp
create mode 100644 lib/gpu/lal_lj_expand_coul_long.cu
create mode 100644 lib/gpu/lal_lj_expand_coul_long.h
create mode 100644 lib/gpu/lal_lj_expand_coul_long_ext.cpp
create mode 100644 src/GPU/pair_lj_cut_dipole_long_gpu.cpp
create mode 100644 src/GPU/pair_lj_cut_dipole_long_gpu.h
create mode 100644 src/GPU/pair_lj_expand_coul_long_gpu.cpp
create mode 100644 src/GPU/pair_lj_expand_coul_long_gpu.h
create mode 100644 src/USER-MISC/pair_lj_expand_coul_long.cpp
create mode 100644 src/USER-MISC/pair_lj_expand_coul_long.h
diff --git a/lib/gpu/Makefile.linux_multi b/lib/gpu/Makefile.linux_multi
new file mode 100644
index 0000000000..02b405ac35
--- /dev/null
+++ b/lib/gpu/Makefile.linux_multi
@@ -0,0 +1,53 @@
+# /* ----------------------------------------------------------------------
+# Generic Linux Makefile for CUDA
+# - Change CUDA_ARCH for your GPU
+# ------------------------------------------------------------------------- */
+
+# which file will be copied to Makefile.lammps
+
+EXTRAMAKE = Makefile.lammps.standard
+
+ifeq ($(CUDA_HOME),)
+CUDA_HOME = /usr/local/cuda
+endif
+
+NVCC = nvcc
+
+# Kepler CUDA
+#CUDA_ARCH = -arch=sm_35
+# newer CUDA
+#CUDA_ARCH = -arch=sm_13
+# older CUDA
+#CUDA_ARCH = -arch=sm_10 -DCUDA_PRE_THREE
+
+CUDA_ARCH = -gencode arch=compute_60,code=sm_60 -gencode arch=compute_61,code=sm_61
+
+# this setting should match LAMMPS Makefile
+# one of LAMMPS_SMALLBIG (default), LAMMPS_BIGBIG and LAMMPS_SMALLSMALL
+
+LMP_INC = -DLAMMPS_SMALLBIG
+
+# precision for GPU calculations
+# -D_SINGLE_SINGLE # Single precision for all calculations
+# -D_DOUBLE_DOUBLE # Double precision for all calculations
+# -D_SINGLE_DOUBLE # Accumulation of forces, etc. in double
+
+CUDA_PRECISION = -D_SINGLE_DOUBLE
+
+CUDA_INCLUDE = -I$(CUDA_HOME)/include
+CUDA_LIB = -L$(CUDA_HOME)/lib64
+CUDA_OPTS = -DUNIX -O3 -Xptxas -v --use_fast_math $(LMP_INC) -Xcompiler "-fPIC -std=c++98"
+
+CUDR_CPP = mpicxx -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK -DOMPI_SKIP_MPICXX=1 -fPIC
+CUDR_OPTS = -O2 $(LMP_INC) # -xHost -no-prec-div -ansi-alias
+
+BIN_DIR = ./
+OBJ_DIR = ./
+LIB_DIR = ./
+AR = ar
+BSH = /bin/sh
+
+CUDPP_OPT = -DUSE_CUDPP -Icudpp_mini
+
+include Nvidia.makefile_multi
+
diff --git a/lib/gpu/Nvidia.makefile b/lib/gpu/Nvidia.makefile
index 5f692cf66c..21663ad85c 100644
--- a/lib/gpu/Nvidia.makefile
+++ b/lib/gpu/Nvidia.makefile
@@ -77,7 +77,9 @@ OBJS = $(OBJ_DIR)/lal_atom.o $(OBJ_DIR)/lal_ans.o \
$(OBJ_DIR)/lal_coul_debye.o $(OBJ_DIR)/lal_coul_debye_ext.o \
$(OBJ_DIR)/lal_zbl.o $(OBJ_DIR)/lal_zbl_ext.o \
$(OBJ_DIR)/lal_lj_cubic.o $(OBJ_DIR)/lal_lj_cubic_ext.o \
- $(OBJ_DIR)/lal_ufm.o $(OBJ_DIR)/lal_ufm_ext.o
+ $(OBJ_DIR)/lal_ufm.o $(OBJ_DIR)/lal_ufm_ext.o \
+ $(OBJ_DIR)/lal_dipole_long_lj.o $(OBJ_DIR)/lal_dipole_long_lj_ext.o \
+ $(OBJ_DIR)/lal_lj_expand_coul_long.o $(OBJ_DIR)/lal_lj_expand_coul_long_ext.o
CBNS = $(OBJ_DIR)/device.cubin $(OBJ_DIR)/device_cubin.h \
$(OBJ_DIR)/atom.cubin $(OBJ_DIR)/atom_cubin.h \
@@ -133,7 +135,9 @@ CBNS = $(OBJ_DIR)/device.cubin $(OBJ_DIR)/device_cubin.h \
$(OBJ_DIR)/coul_debye.cubin $(OBJ_DIR)/coul_debye_cubin.h \
$(OBJ_DIR)/zbl.cubin $(OBJ_DIR)/zbl_cubin.h \
$(OBJ_DIR)/lj_cubic.cubin $(OBJ_DIR)/lj_cubic_cubin.h \
- $(OBJ_DIR)/ufm.cubin $(OBJ_DIR)/ufm_cubin.h
+ $(OBJ_DIR)/ufm.cubin $(OBJ_DIR)/ufm_cubin.h \
+ $(OBJ_DIR)/dipole_long_lj.cubin $(OBJ_DIR)/dipole_long_lj_cubin.h \
+ $(OBJ_DIR)/lj_expand_coul_long.cubin $(OBJ_DIR)/lj_expand_coul_long_cubin.h
all: $(OBJ_DIR) $(GPU_LIB) $(EXECS)
@@ -809,6 +813,30 @@ $(OBJ_DIR)/lal_lj_cubic.o: $(ALL_H) lal_lj_cubic.h lal_lj_cubic.cpp $(OBJ_DIR)/l
$(OBJ_DIR)/lal_lj_cubic_ext.o: $(ALL_H) lal_lj_cubic.h lal_lj_cubic_ext.cpp lal_base_atomic.h
$(CUDR) -o $@ -c lal_lj_cubic_ext.cpp -I$(OBJ_DIR)
+$(OBJ_DIR)/dipole_long_lj.cubin: lal_dipole_long_lj.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --cubin -DNV_KERNEL -o $@ lal_dipole_long_lj.cu
+
+$(OBJ_DIR)/dipole_long_lj_cubin.h: $(OBJ_DIR)/dipole_long_lj.cubin $(OBJ_DIR)/dipole_long_lj.cubin
+ $(BIN2C) -c -n dipole_long_lj $(OBJ_DIR)/dipole_long_lj.cubin > $(OBJ_DIR)/dipole_long_lj_cubin.h
+
+$(OBJ_DIR)/lal_dipole_long_lj.o: $(ALL_H) lal_dipole_long_lj.h lal_dipole_long_lj.cpp $(OBJ_DIR)/dipole_long_lj_cubin.h $(OBJ_DIR)/lal_base_dipole.o
+ $(CUDR) -o $@ -c lal_dipole_long_lj.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_dipole_long_lj_ext.o: $(ALL_H) lal_dipole_long_lj.h lal_dipole_long_lj_ext.cpp lal_base_dipole.h
+ $(CUDR) -o $@ -c lal_dipole_long_lj_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_expand_coul_long.cubin: lal_lj_expand_coul_long.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --cubin -DNV_KERNEL -o $@ lal_lj_expand_coul_long.cu
+
+$(OBJ_DIR)/lj_expand_coul_long_cubin.h: $(OBJ_DIR)/lj_expand_coul_long.cubin $(OBJ_DIR)/lj_expand_coul_long.cubin
+ $(BIN2C) -c -n lj_expand_coul_long $(OBJ_DIR)/lj_expand_coul_long.cubin > $(OBJ_DIR)/lj_expand_coul_long_cubin.h
+
+$(OBJ_DIR)/lal_lj_expand_coul_long.o: $(ALL_H) lal_lj_expand_coul_long.h lal_lj_expand_coul_long.cpp $(OBJ_DIR)/lj_expand_coul_long_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_lj_expand_coul_long.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_expand_coul_long_ext.o: $(ALL_H) lal_lj_expand_coul_long.h lal_lj_expand_coul_long_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_lj_expand_coul_long_ext.cpp -I$(OBJ_DIR)
+
$(BIN_DIR)/nvc_get_devices: ./geryon/ucl_get_devices.cpp $(NVD_H)
$(CUDR) -o $@ ./geryon/ucl_get_devices.cpp -DUCL_CUDADR $(CUDA_LIB) -lcuda
diff --git a/lib/gpu/Nvidia.makefile_multi b/lib/gpu/Nvidia.makefile_multi
new file mode 100644
index 0000000000..5fb35cce3c
--- /dev/null
+++ b/lib/gpu/Nvidia.makefile_multi
@@ -0,0 +1,854 @@
+CUDA = $(NVCC) $(CUDA_INCLUDE) $(CUDA_OPTS) -Icudpp_mini $(CUDA_ARCH) \
+ $(CUDA_PRECISION)
+CUDR = $(CUDR_CPP) $(CUDR_OPTS) $(CUDA_PRECISION) $(CUDA_INCLUDE) \
+ $(CUDPP_OPT)
+CUDA_LINK = $(CUDA_LIB) -lcudart
+BIN2C = $(CUDA_HOME)/bin/bin2c
+
+GPU_LIB = $(LIB_DIR)/libgpu.a
+
+# Headers for Geryon
+UCL_H = $(wildcard ./geryon/ucl*.h)
+NVC_H = $(wildcard ./geryon/nvc*.h) $(UCL_H)
+NVD_H = $(wildcard ./geryon/nvd*.h) $(UCL_H) lal_preprocessor.h
+# Headers for Pair Stuff
+PAIR_H = lal_atom.h lal_answer.h lal_neighbor_shared.h \
+ lal_neighbor.h lal_precision.h lal_device.h \
+ lal_balance.h lal_pppm.h
+
+ALL_H = $(NVD_H) $(PAIR_H)
+
+EXECS = $(BIN_DIR)/nvc_get_devices
+ifdef CUDPP_OPT
+CUDPP = $(OBJ_DIR)/cudpp.o $(OBJ_DIR)/cudpp_plan.o \
+ $(OBJ_DIR)/cudpp_maximal_launch.o $(OBJ_DIR)/cudpp_plan_manager.o \
+ $(OBJ_DIR)/radixsort_app.cu_o $(OBJ_DIR)/scan_app.cu_o
+endif
+OBJS = $(OBJ_DIR)/lal_atom.o $(OBJ_DIR)/lal_ans.o \
+ $(OBJ_DIR)/lal_neighbor.o $(OBJ_DIR)/lal_neighbor_shared.o \
+ $(OBJ_DIR)/lal_device.o $(OBJ_DIR)/lal_base_atomic.o \
+ $(OBJ_DIR)/lal_base_charge.o $(OBJ_DIR)/lal_base_ellipsoid.o \
+ $(OBJ_DIR)/lal_base_dipole.o $(OBJ_DIR)/lal_base_three.o \
+ $(OBJ_DIR)/lal_base_dpd.o \
+ $(OBJ_DIR)/lal_pppm.o $(OBJ_DIR)/lal_pppm_ext.o \
+ $(OBJ_DIR)/lal_gayberne.o $(OBJ_DIR)/lal_gayberne_ext.o \
+ $(OBJ_DIR)/lal_re_squared.o $(OBJ_DIR)/lal_re_squared_ext.o \
+ $(OBJ_DIR)/lal_lj.o $(OBJ_DIR)/lal_lj_ext.o \
+ $(OBJ_DIR)/lal_lj96.o $(OBJ_DIR)/lal_lj96_ext.o \
+ $(OBJ_DIR)/lal_lj_expand.o $(OBJ_DIR)/lal_lj_expand_ext.o \
+ $(OBJ_DIR)/lal_lj_coul.o $(OBJ_DIR)/lal_lj_coul_ext.o \
+ $(OBJ_DIR)/lal_lj_coul_long.o $(OBJ_DIR)/lal_lj_coul_long_ext.o \
+ $(OBJ_DIR)/lal_lj_dsf.o $(OBJ_DIR)/lal_lj_dsf_ext.o \
+ $(OBJ_DIR)/lal_lj_class2_long.o $(OBJ_DIR)/lal_lj_class2_long_ext.o \
+ $(OBJ_DIR)/lal_coul_long.o $(OBJ_DIR)/lal_coul_long_ext.o \
+ $(OBJ_DIR)/lal_morse.o $(OBJ_DIR)/lal_morse_ext.o \
+ $(OBJ_DIR)/lal_charmm_long.o $(OBJ_DIR)/lal_charmm_long_ext.o \
+ $(OBJ_DIR)/lal_lj_sdk.o $(OBJ_DIR)/lal_lj_sdk_ext.o \
+ $(OBJ_DIR)/lal_lj_sdk_long.o $(OBJ_DIR)/lal_lj_sdk_long_ext.o \
+ $(OBJ_DIR)/lal_eam.o $(OBJ_DIR)/lal_eam_ext.o \
+ $(OBJ_DIR)/lal_eam_fs_ext.o $(OBJ_DIR)/lal_eam_alloy_ext.o \
+ $(OBJ_DIR)/lal_buck.o $(OBJ_DIR)/lal_buck_ext.o \
+ $(OBJ_DIR)/lal_buck_coul.o $(OBJ_DIR)/lal_buck_coul_ext.o \
+ $(OBJ_DIR)/lal_buck_coul_long.o $(OBJ_DIR)/lal_buck_coul_long_ext.o \
+ $(OBJ_DIR)/lal_table.o $(OBJ_DIR)/lal_table_ext.o \
+ $(OBJ_DIR)/lal_yukawa.o $(OBJ_DIR)/lal_yukawa_ext.o \
+ $(OBJ_DIR)/lal_born.o $(OBJ_DIR)/lal_born_ext.o \
+ $(OBJ_DIR)/lal_born_coul_wolf.o $(OBJ_DIR)/lal_born_coul_wolf_ext.o \
+ $(OBJ_DIR)/lal_born_coul_long.o $(OBJ_DIR)/lal_born_coul_long_ext.o \
+ $(OBJ_DIR)/lal_dipole_lj.o $(OBJ_DIR)/lal_dipole_lj_ext.o \
+ $(OBJ_DIR)/lal_dipole_lj_sf.o $(OBJ_DIR)/lal_dipole_lj_sf_ext.o \
+ $(OBJ_DIR)/lal_colloid.o $(OBJ_DIR)/lal_colloid_ext.o \
+ $(OBJ_DIR)/lal_gauss.o $(OBJ_DIR)/lal_gauss_ext.o \
+ $(OBJ_DIR)/lal_yukawa_colloid.o $(OBJ_DIR)/lal_yukawa_colloid_ext.o \
+ $(OBJ_DIR)/lal_lj_coul_debye.o $(OBJ_DIR)/lal_lj_coul_debye_ext.o \
+ $(OBJ_DIR)/lal_coul_dsf.o $(OBJ_DIR)/lal_coul_dsf_ext.o \
+ $(OBJ_DIR)/lal_sw.o $(OBJ_DIR)/lal_sw_ext.o \
+ $(OBJ_DIR)/lal_vashishta.o $(OBJ_DIR)/lal_vashishta_ext.o \
+ $(OBJ_DIR)/lal_beck.o $(OBJ_DIR)/lal_beck_ext.o \
+ $(OBJ_DIR)/lal_mie.o $(OBJ_DIR)/lal_mie_ext.o \
+ $(OBJ_DIR)/lal_soft.o $(OBJ_DIR)/lal_soft_ext.o \
+ $(OBJ_DIR)/lal_lj_coul_msm.o $(OBJ_DIR)/lal_lj_coul_msm_ext.o \
+ $(OBJ_DIR)/lal_lj_gromacs.o $(OBJ_DIR)/lal_lj_gromacs_ext.o \
+ $(OBJ_DIR)/lal_dpd.o $(OBJ_DIR)/lal_dpd_ext.o \
+ $(OBJ_DIR)/lal_tersoff.o $(OBJ_DIR)/lal_tersoff_ext.o \
+ $(OBJ_DIR)/lal_tersoff_zbl.o $(OBJ_DIR)/lal_tersoff_zbl_ext.o \
+ $(OBJ_DIR)/lal_tersoff_mod.o $(OBJ_DIR)/lal_tersoff_mod_ext.o \
+ $(OBJ_DIR)/lal_coul.o $(OBJ_DIR)/lal_coul_ext.o \
+ $(OBJ_DIR)/lal_coul_debye.o $(OBJ_DIR)/lal_coul_debye_ext.o \
+ $(OBJ_DIR)/lal_zbl.o $(OBJ_DIR)/lal_zbl_ext.o \
+ $(OBJ_DIR)/lal_lj_cubic.o $(OBJ_DIR)/lal_lj_cubic_ext.o \
+ $(OBJ_DIR)/lal_ufm.o $(OBJ_DIR)/lal_ufm_ext.o \
+ $(OBJ_DIR)/lal_dipole_long_lj.o $(OBJ_DIR)/lal_dipole_long_lj_ext.o \
+ $(OBJ_DIR)/lal_lj_expand_coul_long.o $(OBJ_DIR)/lal_lj_expand_coul_long_ext.o
+
+CBNS = $(OBJ_DIR)/device.cubin $(OBJ_DIR)/device_cubin.h \
+ $(OBJ_DIR)/atom.cubin $(OBJ_DIR)/atom_cubin.h \
+ $(OBJ_DIR)/neighbor_cpu.cubin $(OBJ_DIR)/neighbor_cpu_cubin.h \
+ $(OBJ_DIR)/neighbor_gpu.cubin $(OBJ_DIR)/neighbor_gpu_cubin.h \
+ $(OBJ_DIR)/pppm_f.cubin $(OBJ_DIR)/pppm_f_cubin.h \
+ $(OBJ_DIR)/pppm_d.cubin $(OBJ_DIR)/pppm_d_cubin.h \
+ $(OBJ_DIR)/ellipsoid_nbor.cubin $(OBJ_DIR)/ellipsoid_nbor_cubin.h \
+ $(OBJ_DIR)/gayberne.cubin $(OBJ_DIR)/gayberne_lj.cubin \
+ $(OBJ_DIR)/gayberne_cubin.h $(OBJ_DIR)/gayberne_lj_cubin.h \
+ $(OBJ_DIR)/re_squared.cubin $(OBJ_DIR)/re_squared_lj.cubin \
+ $(OBJ_DIR)/re_squared_cubin.h $(OBJ_DIR)/re_squared_lj_cubin.h \
+ $(OBJ_DIR)/lj.cubin $(OBJ_DIR)/lj_cubin.h \
+ $(OBJ_DIR)/lj96.cubin $(OBJ_DIR)/lj96_cubin.h \
+ $(OBJ_DIR)/lj_expand.cubin $(OBJ_DIR)/lj_expand_cubin.h \
+ $(OBJ_DIR)/lj_coul.cubin $(OBJ_DIR)/lj_coul_cubin.h \
+ $(OBJ_DIR)/lj_coul_long.cubin $(OBJ_DIR)/lj_coul_long_cubin.h \
+ $(OBJ_DIR)/lj_dsf.cubin $(OBJ_DIR)/lj_dsf_cubin.h \
+ $(OBJ_DIR)/lj_class2_long.cubin $(OBJ_DIR)/lj_class2_long_cubin.h \
+ $(OBJ_DIR)/coul_long.cubin $(OBJ_DIR)/coul_long_cubin.h \
+ $(OBJ_DIR)/morse.cubin $(OBJ_DIR)/morse_cubin.h \
+ $(OBJ_DIR)/charmm_long.cubin $(OBJ_DIR)/charmm_long_cubin.h \
+ $(OBJ_DIR)/lj_sdk.cubin $(OBJ_DIR)/lj_sdk_cubin.h \
+ $(OBJ_DIR)/lj_sdk_long.cubin $(OBJ_DIR)/lj_sdk_long_cubin.h \
+ $(OBJ_DIR)/eam.cubin $(OBJ_DIR)/eam_cubin.h \
+ $(OBJ_DIR)/buck.cubin $(OBJ_DIR)/buck_cubin.h \
+ $(OBJ_DIR)/buck_coul_long.cubin $(OBJ_DIR)/buck_coul_long_cubin.h \
+ $(OBJ_DIR)/buck_coul.cubin $(OBJ_DIR)/buck_coul_cubin.h \
+ $(OBJ_DIR)/table.cubin $(OBJ_DIR)/table_cubin.h \
+ $(OBJ_DIR)/yukawa.cubin $(OBJ_DIR)/yukawa_cubin.h \
+ $(OBJ_DIR)/born.cubin $(OBJ_DIR)/born_cubin.h \
+ $(OBJ_DIR)/born_coul_wolf.cubin $(OBJ_DIR)/born_coul_wolf_cubin.h \
+ $(OBJ_DIR)/born_coul_long.cubin $(OBJ_DIR)/born_coul_long_cubin.h \
+ $(OBJ_DIR)/dipole_lj.cubin $(OBJ_DIR)/dipole_lj_cubin.h \
+ $(OBJ_DIR)/dipole_lj_sf.cubin $(OBJ_DIR)/dipole_lj_sf_cubin.h \
+ $(OBJ_DIR)/colloid.cubin $(OBJ_DIR)/colloid_cubin.h \
+ $(OBJ_DIR)/gauss.cubin $(OBJ_DIR)/gauss_cubin.h \
+ $(OBJ_DIR)/yukawa_colloid.cubin $(OBJ_DIR)/yukawa_colloid_cubin.h \
+ $(OBJ_DIR)/lj_coul_debye.cubin $(OBJ_DIR)/lj_coul_debye_cubin.h \
+ $(OBJ_DIR)/coul_dsf.cubin $(OBJ_DIR)/coul_dsf_cubin.h \
+ $(OBJ_DIR)/sw.cubin $(OBJ_DIR)/sw_cubin.h \
+ $(OBJ_DIR)/vashishta.cubin $(OBJ_DIR)/vashishta_cubin.h \
+ $(OBJ_DIR)/beck.cubin $(OBJ_DIR)/beck_cubin.h \
+ $(OBJ_DIR)/mie.cubin $(OBJ_DIR)/mie_cubin.h \
+ $(OBJ_DIR)/soft.cubin $(OBJ_DIR)/soft_cubin.h \
+ $(OBJ_DIR)/lj_coul_msm.cubin $(OBJ_DIR)/lj_coul_msm_cubin.h \
+ $(OBJ_DIR)/lj_gromacs.cubin $(OBJ_DIR)/lj_gromacs_cubin.h \
+ $(OBJ_DIR)/dpd.cubin $(OBJ_DIR)/dpd_cubin.h \
+ $(OBJ_DIR)/tersoff.cubin $(OBJ_DIR)/tersoff_cubin.h \
+ $(OBJ_DIR)/tersoff_zbl.cubin $(OBJ_DIR)/tersoff_zbl_cubin.h \
+ $(OBJ_DIR)/tersoff_mod.cubin $(OBJ_DIR)/tersoff_mod_cubin.h \
+ $(OBJ_DIR)/coul.cubin $(OBJ_DIR)/coul_cubin.h \
+ $(OBJ_DIR)/coul_debye.cubin $(OBJ_DIR)/coul_debye_cubin.h \
+ $(OBJ_DIR)/zbl.cubin $(OBJ_DIR)/zbl_cubin.h \
+ $(OBJ_DIR)/lj_cubic.cubin $(OBJ_DIR)/lj_cubic_cubin.h \
+ $(OBJ_DIR)/ufm.cubin $(OBJ_DIR)/ufm_cubin.h \
+ $(OBJ_DIR)/dipole_long_lj.cubin $(OBJ_DIR)/dipole_long_lj_cubin.h \
+ $(OBJ_DIR)/lj_expand_coul_long.cubin $(OBJ_DIR)/lj_expand_coul_long_cubin.h
+
+all: $(OBJ_DIR) $(GPU_LIB) $(EXECS)
+
+$(OBJ_DIR):
+ mkdir -p $@
+
+$(OBJ_DIR)/cudpp.o: cudpp_mini/cudpp.cpp
+ $(CUDR) -o $@ -c cudpp_mini/cudpp.cpp -Icudpp_mini
+
+$(OBJ_DIR)/cudpp_plan.o: cudpp_mini/cudpp_plan.cpp
+ $(CUDR) -o $@ -c cudpp_mini/cudpp_plan.cpp -Icudpp_mini
+
+$(OBJ_DIR)/cudpp_maximal_launch.o: cudpp_mini/cudpp_maximal_launch.cpp
+ $(CUDR) -o $@ -c cudpp_mini/cudpp_maximal_launch.cpp -Icudpp_mini
+
+$(OBJ_DIR)/cudpp_plan_manager.o: cudpp_mini/cudpp_plan_manager.cpp
+ $(CUDR) -o $@ -c cudpp_mini/cudpp_plan_manager.cpp -Icudpp_mini
+
+$(OBJ_DIR)/radixsort_app.cu_o: cudpp_mini/radixsort_app.cu
+ $(CUDA) -o $@ -c cudpp_mini/radixsort_app.cu
+
+$(OBJ_DIR)/scan_app.cu_o: cudpp_mini/scan_app.cu
+ $(CUDA) -o $@ -c cudpp_mini/scan_app.cu
+
+$(OBJ_DIR)/atom.cubin: lal_atom.cu lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_atom.cu
+
+$(OBJ_DIR)/atom_cubin.h: $(OBJ_DIR)/atom.cubin
+ $(BIN2C) -c -n atom $(OBJ_DIR)/atom.cubin > $(OBJ_DIR)/atom_cubin.h
+
+$(OBJ_DIR)/lal_atom.o: lal_atom.cpp lal_atom.h $(NVD_H) $(OBJ_DIR)/atom_cubin.h
+ $(CUDR) -o $@ -c lal_atom.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_ans.o: lal_answer.cpp lal_answer.h $(NVD_H)
+ $(CUDR) -o $@ -c lal_answer.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/neighbor_cpu.cubin: lal_neighbor_cpu.cu lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_neighbor_cpu.cu
+
+$(OBJ_DIR)/neighbor_cpu_cubin.h: $(OBJ_DIR)/neighbor_cpu.cubin
+ $(BIN2C) -c -n neighbor_cpu $(OBJ_DIR)/neighbor_cpu.cubin > $(OBJ_DIR)/neighbor_cpu_cubin.h
+
+$(OBJ_DIR)/neighbor_gpu.cubin: lal_neighbor_gpu.cu lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_neighbor_gpu.cu
+
+$(OBJ_DIR)/neighbor_gpu_cubin.h: $(OBJ_DIR)/neighbor_gpu.cubin
+ $(BIN2C) -c -n neighbor_gpu $(OBJ_DIR)/neighbor_gpu.cubin > $(OBJ_DIR)/neighbor_gpu_cubin.h
+
+$(OBJ_DIR)/lal_neighbor_shared.o: lal_neighbor_shared.cpp lal_neighbor_shared.h $(OBJ_DIR)/neighbor_cpu_cubin.h $(OBJ_DIR)/neighbor_gpu_cubin.h $(NVD_H)
+ $(CUDR) -o $@ -c lal_neighbor_shared.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_neighbor.o: lal_neighbor.cpp lal_neighbor.h lal_neighbor_shared.h $(NVD_H)
+ $(CUDR) -o $@ -c lal_neighbor.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/device.cubin: lal_device.cu lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_device.cu
+
+$(OBJ_DIR)/device_cubin.h: $(OBJ_DIR)/device.cubin
+ $(BIN2C) -c -n device $(OBJ_DIR)/device.cubin > $(OBJ_DIR)/device_cubin.h
+
+$(OBJ_DIR)/lal_device.o: lal_device.cpp lal_device.h $(ALL_H) $(OBJ_DIR)/device_cubin.h
+ $(CUDR) -o $@ -c lal_device.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_base_atomic.o: $(ALL_H) lal_base_atomic.h lal_base_atomic.cpp
+ $(CUDR) -o $@ -c lal_base_atomic.cpp
+
+$(OBJ_DIR)/lal_base_charge.o: $(ALL_H) lal_base_charge.h lal_base_charge.cpp
+ $(CUDR) -o $@ -c lal_base_charge.cpp
+
+$(OBJ_DIR)/lal_base_ellipsoid.o: $(ALL_H) lal_base_ellipsoid.h lal_base_ellipsoid.cpp $(OBJ_DIR)/ellipsoid_nbor_cubin.h
+ $(CUDR) -o $@ -c lal_base_ellipsoid.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_base_dipole.o: $(ALL_H) lal_base_dipole.h lal_base_dipole.cpp
+ $(CUDR) -o $@ -c lal_base_dipole.cpp
+
+$(OBJ_DIR)/lal_base_three.o: $(ALL_H) lal_base_three.h lal_base_three.cpp
+ $(CUDR) -o $@ -c lal_base_three.cpp
+
+$(OBJ_DIR)/lal_base_dpd.o: $(ALL_H) lal_base_dpd.h lal_base_dpd.cpp
+ $(CUDR) -o $@ -c lal_base_dpd.cpp
+
+$(OBJ_DIR)/pppm_f.cubin: lal_pppm.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -Dgrdtyp=float -Dgrdtyp4=float4 -o $@ lal_pppm.cu
+
+$(OBJ_DIR)/pppm_f_cubin.h: $(OBJ_DIR)/pppm_f.cubin
+ $(BIN2C) -c -n pppm_f $(OBJ_DIR)/pppm_f.cubin > $(OBJ_DIR)/pppm_f_cubin.h
+
+$(OBJ_DIR)/pppm_d.cubin: lal_pppm.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -Dgrdtyp=double -Dgrdtyp4=double4 -o $@ lal_pppm.cu
+
+$(OBJ_DIR)/pppm_d_cubin.h: $(OBJ_DIR)/pppm_d.cubin
+ $(BIN2C) -c -n pppm_d $(OBJ_DIR)/pppm_d.cubin > $(OBJ_DIR)/pppm_d_cubin.h
+
+$(OBJ_DIR)/lal_pppm.o: $(ALL_H) lal_pppm.h lal_pppm.cpp $(OBJ_DIR)/pppm_f_cubin.h $(OBJ_DIR)/pppm_d_cubin.h
+ $(CUDR) -o $@ -c lal_pppm.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_pppm_ext.o: $(ALL_H) lal_pppm.h lal_pppm_ext.cpp
+ $(CUDR) -o $@ -c lal_pppm_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/ellipsoid_nbor.cubin: lal_ellipsoid_nbor.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_ellipsoid_nbor.cu
+
+$(OBJ_DIR)/ellipsoid_nbor_cubin.h: $(OBJ_DIR)/ellipsoid_nbor.cubin
+ $(BIN2C) -c -n ellipsoid_nbor $(OBJ_DIR)/ellipsoid_nbor.cubin > $(OBJ_DIR)/ellipsoid_nbor_cubin.h
+
+$(OBJ_DIR)/gayberne.cubin: lal_gayberne.cu lal_precision.h lal_ellipsoid_extra.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_gayberne.cu
+
+$(OBJ_DIR)/gayberne_lj.cubin: lal_gayberne_lj.cu lal_precision.h lal_ellipsoid_extra.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_gayberne_lj.cu
+
+$(OBJ_DIR)/gayberne_cubin.h: $(OBJ_DIR)/gayberne.cubin
+ $(BIN2C) -c -n gayberne $(OBJ_DIR)/gayberne.cubin > $(OBJ_DIR)/gayberne_cubin.h
+
+$(OBJ_DIR)/gayberne_lj_cubin.h: $(OBJ_DIR)/gayberne_lj.cubin
+ $(BIN2C) -c -n gayberne_lj $(OBJ_DIR)/gayberne_lj.cubin > $(OBJ_DIR)/gayberne_lj_cubin.h
+
+$(OBJ_DIR)/lal_gayberne.o: $(ALL_H) lal_gayberne.h lal_gayberne.cpp $(OBJ_DIR)/gayberne_cubin.h $(OBJ_DIR)/gayberne_lj_cubin.h $(OBJ_DIR)/lal_base_ellipsoid.o
+ $(CUDR) -o $@ -c lal_gayberne.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_gayberne_ext.o: $(ALL_H) $(OBJ_DIR)/lal_gayberne.o lal_gayberne_ext.cpp
+ $(CUDR) -o $@ -c lal_gayberne_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/re_squared.cubin: lal_re_squared.cu lal_precision.h lal_ellipsoid_extra.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_re_squared.cu
+
+$(OBJ_DIR)/re_squared_lj.cubin: lal_re_squared_lj.cu lal_precision.h lal_ellipsoid_extra.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_re_squared_lj.cu
+
+$(OBJ_DIR)/re_squared_cubin.h: $(OBJ_DIR)/re_squared.cubin
+ $(BIN2C) -c -n re_squared $(OBJ_DIR)/re_squared.cubin > $(OBJ_DIR)/re_squared_cubin.h
+
+$(OBJ_DIR)/re_squared_lj_cubin.h: $(OBJ_DIR)/re_squared_lj.cubin
+ $(BIN2C) -c -n re_squared_lj $(OBJ_DIR)/re_squared_lj.cubin > $(OBJ_DIR)/re_squared_lj_cubin.h
+
+$(OBJ_DIR)/lal_re_squared.o: $(ALL_H) lal_re_squared.h lal_re_squared.cpp $(OBJ_DIR)/re_squared_cubin.h $(OBJ_DIR)/re_squared_lj_cubin.h $(OBJ_DIR)/lal_base_ellipsoid.o
+ $(CUDR) -o $@ -c lal_re_squared.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_re_squared_ext.o: $(ALL_H) $(OBJ_DIR)/lal_re_squared.o lal_re_squared_ext.cpp
+ $(CUDR) -o $@ -c lal_re_squared_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj.cubin: lal_lj.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj.cu
+
+$(OBJ_DIR)/lj_cubin.h: $(OBJ_DIR)/lj.cubin $(OBJ_DIR)/lj.cubin
+ $(BIN2C) -c -n lj $(OBJ_DIR)/lj.cubin > $(OBJ_DIR)/lj_cubin.h
+
+$(OBJ_DIR)/lal_lj.o: $(ALL_H) lal_lj.h lal_lj.cpp $(OBJ_DIR)/lj_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_lj.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_ext.o: $(ALL_H) lal_lj.h lal_lj_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_lj_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_coul.cubin: lal_lj_coul.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_coul.cu
+
+$(OBJ_DIR)/lj_coul_cubin.h: $(OBJ_DIR)/lj_coul.cubin $(OBJ_DIR)/lj_coul.cubin
+ $(BIN2C) -c -n lj_coul $(OBJ_DIR)/lj_coul.cubin > $(OBJ_DIR)/lj_coul_cubin.h
+
+$(OBJ_DIR)/lal_lj_coul.o: $(ALL_H) lal_lj_coul.h lal_lj_coul.cpp $(OBJ_DIR)/lj_coul_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_lj_coul.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_coul_ext.o: $(ALL_H) lal_lj_coul.h lal_lj_coul_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_lj_coul_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_class2_long.cubin: lal_lj_class2_long.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_class2_long.cu
+
+$(OBJ_DIR)/lj_class2_long_cubin.h: $(OBJ_DIR)/lj_class2_long.cubin $(OBJ_DIR)/lj_class2_long.cubin
+ $(BIN2C) -c -n lj_class2_long $(OBJ_DIR)/lj_class2_long.cubin > $(OBJ_DIR)/lj_class2_long_cubin.h
+
+$(OBJ_DIR)/lal_lj_class2_long.o: $(ALL_H) lal_lj_class2_long.h lal_lj_class2_long.cpp $(OBJ_DIR)/lj_class2_long_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_lj_class2_long.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_class2_long_ext.o: $(ALL_H) lal_lj_class2_long.h lal_lj_class2_long_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_lj_class2_long_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/coul_long.cubin: lal_coul_long.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_coul_long.cu
+
+$(OBJ_DIR)/coul_long_cubin.h: $(OBJ_DIR)/coul_long.cubin $(OBJ_DIR)/coul_long.cubin
+ $(BIN2C) -c -n coul_long $(OBJ_DIR)/coul_long.cubin > $(OBJ_DIR)/coul_long_cubin.h
+
+$(OBJ_DIR)/lal_coul_long.o: $(ALL_H) lal_coul_long.h lal_coul_long.cpp $(OBJ_DIR)/coul_long_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_coul_long.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_coul_long_ext.o: $(ALL_H) lal_coul_long.h lal_coul_long_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_coul_long_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_coul_long.cubin: lal_lj_coul_long.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_coul_long.cu
+
+$(OBJ_DIR)/lj_coul_long_cubin.h: $(OBJ_DIR)/lj_coul_long.cubin $(OBJ_DIR)/lj_coul_long.cubin
+ $(BIN2C) -c -n lj_coul_long $(OBJ_DIR)/lj_coul_long.cubin > $(OBJ_DIR)/lj_coul_long_cubin.h
+
+$(OBJ_DIR)/lal_lj_coul_long.o: $(ALL_H) lal_lj_coul_long.h lal_lj_coul_long.cpp $(OBJ_DIR)/lj_coul_long_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_lj_coul_long.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_coul_long_ext.o: $(ALL_H) lal_lj_coul_long.h lal_lj_coul_long_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_lj_coul_long_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_dsf.cubin: lal_lj_dsf.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_dsf.cu
+
+$(OBJ_DIR)/lj_dsf_cubin.h: $(OBJ_DIR)/lj_dsf.cubin $(OBJ_DIR)/lj_dsf.cubin
+ $(BIN2C) -c -n lj_dsf $(OBJ_DIR)/lj_dsf.cubin > $(OBJ_DIR)/lj_dsf_cubin.h
+
+$(OBJ_DIR)/lal_lj_dsf.o: $(ALL_H) lal_lj_dsf.h lal_lj_dsf.cpp $(OBJ_DIR)/lj_dsf_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_lj_dsf.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_dsf_ext.o: $(ALL_H) lal_lj_dsf.h lal_lj_dsf_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_lj_dsf_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/morse.cubin: lal_morse.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_morse.cu
+
+$(OBJ_DIR)/morse_cubin.h: $(OBJ_DIR)/morse.cubin $(OBJ_DIR)/morse.cubin
+ $(BIN2C) -c -n morse $(OBJ_DIR)/morse.cubin > $(OBJ_DIR)/morse_cubin.h
+
+$(OBJ_DIR)/lal_morse.o: $(ALL_H) lal_morse.h lal_morse.cpp $(OBJ_DIR)/morse_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_morse.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_morse_ext.o: $(ALL_H) lal_morse.h lal_morse_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_morse_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/charmm_long.cubin: lal_charmm_long.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_charmm_long.cu
+
+$(OBJ_DIR)/charmm_long_cubin.h: $(OBJ_DIR)/charmm_long.cubin $(OBJ_DIR)/charmm_long.cubin
+ $(BIN2C) -c -n charmm_long $(OBJ_DIR)/charmm_long.cubin > $(OBJ_DIR)/charmm_long_cubin.h
+
+$(OBJ_DIR)/lal_charmm_long.o: $(ALL_H) lal_charmm_long.h lal_charmm_long.cpp $(OBJ_DIR)/charmm_long_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_charmm_long.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_charmm_long_ext.o: $(ALL_H) lal_charmm_long.h lal_charmm_long_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_charmm_long_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj96.cubin: lal_lj96.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj96.cu
+
+$(OBJ_DIR)/lj96_cubin.h: $(OBJ_DIR)/lj96.cubin $(OBJ_DIR)/lj96.cubin
+ $(BIN2C) -c -n lj96 $(OBJ_DIR)/lj96.cubin > $(OBJ_DIR)/lj96_cubin.h
+
+$(OBJ_DIR)/lal_lj96.o: $(ALL_H) lal_lj96.h lal_lj96.cpp $(OBJ_DIR)/lj96_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_lj96.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj96_ext.o: $(ALL_H) lal_lj96.h lal_lj96_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_lj96_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_expand.cubin: lal_lj_expand.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_expand.cu
+
+$(OBJ_DIR)/lj_expand_cubin.h: $(OBJ_DIR)/lj_expand.cubin $(OBJ_DIR)/lj_expand.cubin
+ $(BIN2C) -c -n lj_expand $(OBJ_DIR)/lj_expand.cubin > $(OBJ_DIR)/lj_expand_cubin.h
+
+$(OBJ_DIR)/lal_lj_expand.o: $(ALL_H) lal_lj_expand.h lal_lj_expand.cpp $(OBJ_DIR)/lj_expand_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_lj_expand.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_expand_ext.o: $(ALL_H) lal_lj_expand.h lal_lj_expand_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_lj_expand_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_sdk.cubin: lal_lj_sdk.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_sdk.cu
+
+$(OBJ_DIR)/lj_sdk_cubin.h: $(OBJ_DIR)/lj_sdk.cubin $(OBJ_DIR)/lj_sdk.cubin
+ $(BIN2C) -c -n lj_sdk $(OBJ_DIR)/lj_sdk.cubin > $(OBJ_DIR)/lj_sdk_cubin.h
+
+$(OBJ_DIR)/lal_lj_sdk.o: $(ALL_H) lal_lj_sdk.h lal_lj_sdk.cpp $(OBJ_DIR)/lj_sdk_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_lj_sdk.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_sdk_ext.o: $(ALL_H) lal_lj_sdk.h lal_lj_sdk_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_lj_sdk_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_sdk_long.cubin: lal_lj_sdk_long.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_sdk_long.cu
+
+$(OBJ_DIR)/lj_sdk_long_cubin.h: $(OBJ_DIR)/lj_sdk_long.cubin $(OBJ_DIR)/lj_sdk_long.cubin
+ $(BIN2C) -c -n lj_sdk_long $(OBJ_DIR)/lj_sdk_long.cubin > $(OBJ_DIR)/lj_sdk_long_cubin.h
+
+$(OBJ_DIR)/lal_lj_sdk_long.o: $(ALL_H) lal_lj_sdk_long.h lal_lj_sdk_long.cpp $(OBJ_DIR)/lj_sdk_long_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_lj_sdk_long.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_sdk_long_ext.o: $(ALL_H) lal_lj_sdk_long.h lal_lj_sdk_long_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_lj_sdk_long_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/eam.cubin: lal_eam.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_eam.cu
+
+$(OBJ_DIR)/eam_cubin.h: $(OBJ_DIR)/eam.cubin $(OBJ_DIR)/eam.cubin
+ $(BIN2C) -c -n eam $(OBJ_DIR)/eam.cubin > $(OBJ_DIR)/eam_cubin.h
+
+$(OBJ_DIR)/lal_eam.o: $(ALL_H) lal_eam.h lal_eam.cpp $(OBJ_DIR)/eam_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_eam.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_eam_ext.o: $(ALL_H) lal_eam.h lal_eam_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_eam_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_eam_fs_ext.o: $(ALL_H) lal_eam.h lal_eam_fs_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_eam_fs_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_eam_alloy_ext.o: $(ALL_H) lal_eam.h lal_eam_alloy_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_eam_alloy_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/buck.cubin: lal_buck.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_buck.cu
+
+$(OBJ_DIR)/buck_cubin.h: $(OBJ_DIR)/buck.cubin $(OBJ_DIR)/buck.cubin
+ $(BIN2C) -c -n buck $(OBJ_DIR)/buck.cubin > $(OBJ_DIR)/buck_cubin.h
+
+$(OBJ_DIR)/lal_buck.o: $(ALL_H) lal_buck.h lal_buck.cpp $(OBJ_DIR)/buck_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_buck.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_buck_ext.o: $(ALL_H) lal_buck.h lal_buck_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_buck_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/buck_coul.cubin: lal_buck_coul.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_buck_coul.cu
+
+$(OBJ_DIR)/buck_coul_cubin.h: $(OBJ_DIR)/buck_coul.cubin $(OBJ_DIR)/buck_coul.cubin
+ $(BIN2C) -c -n buck_coul $(OBJ_DIR)/buck_coul.cubin > $(OBJ_DIR)/buck_coul_cubin.h
+
+$(OBJ_DIR)/lal_buck_coul.o: $(ALL_H) lal_buck_coul.h lal_buck_coul.cpp $(OBJ_DIR)/buck_coul_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_buck_coul.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_buck_coul_ext.o: $(ALL_H) lal_buck_coul.h lal_buck_coul_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_buck_coul_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/buck_coul_long.cubin: lal_buck_coul_long.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_buck_coul_long.cu
+
+$(OBJ_DIR)/buck_coul_long_cubin.h: $(OBJ_DIR)/buck_coul_long.cubin $(OBJ_DIR)/buck_coul_long.cubin
+ $(BIN2C) -c -n buck_coul_long $(OBJ_DIR)/buck_coul_long.cubin > $(OBJ_DIR)/buck_coul_long_cubin.h
+
+$(OBJ_DIR)/lal_buck_coul_long.o: $(ALL_H) lal_buck_coul_long.h lal_buck_coul_long.cpp $(OBJ_DIR)/buck_coul_long_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_buck_coul_long.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_buck_coul_long_ext.o: $(ALL_H) lal_buck_coul_long.h lal_buck_coul_long_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_buck_coul_long_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/table.cubin: lal_table.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_table.cu
+
+$(OBJ_DIR)/table_cubin.h: $(OBJ_DIR)/table.cubin $(OBJ_DIR)/table.cubin
+ $(BIN2C) -c -n table $(OBJ_DIR)/table.cubin > $(OBJ_DIR)/table_cubin.h
+
+$(OBJ_DIR)/lal_table.o: $(ALL_H) lal_table.h lal_table.cpp $(OBJ_DIR)/table_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_table.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_table_ext.o: $(ALL_H) lal_table.h lal_table_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_table_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/yukawa.cubin: lal_yukawa.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_yukawa.cu
+
+$(OBJ_DIR)/yukawa_cubin.h: $(OBJ_DIR)/yukawa.cubin $(OBJ_DIR)/yukawa.cubin
+ $(BIN2C) -c -n yukawa $(OBJ_DIR)/yukawa.cubin > $(OBJ_DIR)/yukawa_cubin.h
+
+$(OBJ_DIR)/lal_yukawa.o: $(ALL_H) lal_yukawa.h lal_yukawa.cpp $(OBJ_DIR)/yukawa_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_yukawa.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_yukawa_ext.o: $(ALL_H) lal_yukawa.h lal_yukawa_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_yukawa_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/born.cubin: lal_born.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_born.cu
+
+$(OBJ_DIR)/born_cubin.h: $(OBJ_DIR)/born.cubin $(OBJ_DIR)/born.cubin
+ $(BIN2C) -c -n born $(OBJ_DIR)/born.cubin > $(OBJ_DIR)/born_cubin.h
+
+$(OBJ_DIR)/lal_born.o: $(ALL_H) lal_born.h lal_born.cpp $(OBJ_DIR)/born_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_born.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_born_ext.o: $(ALL_H) lal_born.h lal_born_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_born_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/born_coul_wolf.cubin: lal_born_coul_wolf.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_born_coul_wolf.cu
+
+$(OBJ_DIR)/born_coul_wolf_cubin.h: $(OBJ_DIR)/born_coul_wolf.cubin $(OBJ_DIR)/born_coul_wolf.cubin
+ $(BIN2C) -c -n born_coul_wolf $(OBJ_DIR)/born_coul_wolf.cubin > $(OBJ_DIR)/born_coul_wolf_cubin.h
+
+$(OBJ_DIR)/lal_born_coul_wolf.o: $(ALL_H) lal_born_coul_wolf.h lal_born_coul_wolf.cpp $(OBJ_DIR)/born_coul_wolf_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_born_coul_wolf.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_born_coul_wolf_ext.o: $(ALL_H) lal_born_coul_wolf.h lal_born_coul_wolf_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_born_coul_wolf_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/born_coul_long.cubin: lal_born_coul_long.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_born_coul_long.cu
+
+$(OBJ_DIR)/born_coul_long_cubin.h: $(OBJ_DIR)/born_coul_long.cubin $(OBJ_DIR)/born_coul_long.cubin
+ $(BIN2C) -c -n born_coul_long $(OBJ_DIR)/born_coul_long.cubin > $(OBJ_DIR)/born_coul_long_cubin.h
+
+$(OBJ_DIR)/lal_born_coul_long.o: $(ALL_H) lal_born_coul_long.h lal_born_coul_long.cpp $(OBJ_DIR)/born_coul_long_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_born_coul_long.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_born_coul_long_ext.o: $(ALL_H) lal_born_coul_long.h lal_born_coul_long_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_born_coul_long_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/dipole_lj.cubin: lal_dipole_lj.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_dipole_lj.cu
+
+$(OBJ_DIR)/dipole_lj_cubin.h: $(OBJ_DIR)/dipole_lj.cubin $(OBJ_DIR)/dipole_lj.cubin
+ $(BIN2C) -c -n dipole_lj $(OBJ_DIR)/dipole_lj.cubin > $(OBJ_DIR)/dipole_lj_cubin.h
+
+$(OBJ_DIR)/lal_dipole_lj.o: $(ALL_H) lal_dipole_lj.h lal_dipole_lj.cpp $(OBJ_DIR)/dipole_lj_cubin.h $(OBJ_DIR)/lal_base_dipole.o
+ $(CUDR) -o $@ -c lal_dipole_lj.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_dipole_lj_ext.o: $(ALL_H) lal_dipole_lj.h lal_dipole_lj_ext.cpp lal_base_dipole.h
+ $(CUDR) -o $@ -c lal_dipole_lj_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/dipole_lj_sf.cubin: lal_dipole_lj_sf.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_dipole_lj_sf.cu
+
+$(OBJ_DIR)/dipole_lj_sf_cubin.h: $(OBJ_DIR)/dipole_lj_sf.cubin $(OBJ_DIR)/dipole_lj_sf.cubin
+ $(BIN2C) -c -n dipole_lj_sf $(OBJ_DIR)/dipole_lj_sf.cubin > $(OBJ_DIR)/dipole_lj_sf_cubin.h
+
+$(OBJ_DIR)/lal_dipole_lj_sf.o: $(ALL_H) lal_dipole_lj_sf.h lal_dipole_lj_sf.cpp $(OBJ_DIR)/dipole_lj_sf_cubin.h $(OBJ_DIR)/lal_base_dipole.o
+ $(CUDR) -o $@ -c lal_dipole_lj_sf.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_dipole_lj_sf_ext.o: $(ALL_H) lal_dipole_lj_sf.h lal_dipole_lj_sf_ext.cpp lal_base_dipole.h
+ $(CUDR) -o $@ -c lal_dipole_lj_sf_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/colloid.cubin: lal_colloid.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_colloid.cu
+
+$(OBJ_DIR)/colloid_cubin.h: $(OBJ_DIR)/colloid.cubin $(OBJ_DIR)/colloid.cubin
+ $(BIN2C) -c -n colloid $(OBJ_DIR)/colloid.cubin > $(OBJ_DIR)/colloid_cubin.h
+
+$(OBJ_DIR)/lal_colloid.o: $(ALL_H) lal_colloid.h lal_colloid.cpp $(OBJ_DIR)/colloid_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_colloid.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_colloid_ext.o: $(ALL_H) lal_colloid.h lal_colloid_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_colloid_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/gauss.cubin: lal_gauss.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_gauss.cu
+
+$(OBJ_DIR)/gauss_cubin.h: $(OBJ_DIR)/gauss.cubin $(OBJ_DIR)/gauss.cubin
+ $(BIN2C) -c -n gauss $(OBJ_DIR)/gauss.cubin > $(OBJ_DIR)/gauss_cubin.h
+
+$(OBJ_DIR)/lal_gauss.o: $(ALL_H) lal_gauss.h lal_gauss.cpp $(OBJ_DIR)/gauss_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_gauss.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_gauss_ext.o: $(ALL_H) lal_gauss.h lal_gauss_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_gauss_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/yukawa_colloid.cubin: lal_yukawa_colloid.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_yukawa_colloid.cu
+
+$(OBJ_DIR)/yukawa_colloid_cubin.h: $(OBJ_DIR)/yukawa_colloid.cubin $(OBJ_DIR)/yukawa_colloid.cubin
+ $(BIN2C) -c -n yukawa_colloid $(OBJ_DIR)/yukawa_colloid.cubin > $(OBJ_DIR)/yukawa_colloid_cubin.h
+
+$(OBJ_DIR)/lal_yukawa_colloid.o: $(ALL_H) lal_yukawa_colloid.h lal_yukawa_colloid.cpp $(OBJ_DIR)/yukawa_colloid_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_yukawa_colloid.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_yukawa_colloid_ext.o: $(ALL_H) lal_yukawa_colloid.h lal_yukawa_colloid_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_yukawa_colloid_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_coul_debye.cubin: lal_lj_coul_debye.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_coul_debye.cu
+
+$(OBJ_DIR)/lj_coul_debye_cubin.h: $(OBJ_DIR)/lj_coul_debye.cubin $(OBJ_DIR)/lj_coul_debye.cubin
+ $(BIN2C) -c -n lj_coul_debye $(OBJ_DIR)/lj_coul_debye.cubin > $(OBJ_DIR)/lj_coul_debye_cubin.h
+
+$(OBJ_DIR)/lal_lj_coul_debye.o: $(ALL_H) lal_lj_coul_debye.h lal_lj_coul_debye.cpp $(OBJ_DIR)/lj_coul_debye_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_lj_coul_debye.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_coul_debye_ext.o: $(ALL_H) lal_lj_coul_debye.h lal_lj_coul_debye_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_lj_coul_debye_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/coul_dsf.cubin: lal_coul_dsf.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_coul_dsf.cu
+
+$(OBJ_DIR)/coul_dsf_cubin.h: $(OBJ_DIR)/coul_dsf.cubin $(OBJ_DIR)/coul_dsf.cubin
+ $(BIN2C) -c -n coul_dsf $(OBJ_DIR)/coul_dsf.cubin > $(OBJ_DIR)/coul_dsf_cubin.h
+
+$(OBJ_DIR)/lal_coul_dsf.o: $(ALL_H) lal_coul_dsf.h lal_coul_dsf.cpp $(OBJ_DIR)/coul_dsf_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_coul_dsf.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_coul_dsf_ext.o: $(ALL_H) lal_coul_dsf.h lal_coul_dsf_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_coul_dsf_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/vashishta.cubin: lal_vashishta.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_vashishta.cu
+
+$(OBJ_DIR)/vashishta_cubin.h: $(OBJ_DIR)/vashishta.cubin $(OBJ_DIR)/vashishta.cubin
+ $(BIN2C) -c -n vashishta $(OBJ_DIR)/vashishta.cubin > $(OBJ_DIR)/vashishta_cubin.h
+
+$(OBJ_DIR)/lal_vashishta.o: $(ALL_H) lal_vashishta.h lal_vashishta.cpp $(OBJ_DIR)/vashishta_cubin.h $(OBJ_DIR)/lal_base_three.o
+ $(CUDR) -o $@ -c lal_vashishta.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_vashishta_ext.o: $(ALL_H) lal_vashishta.h lal_vashishta_ext.cpp lal_base_three.h
+ $(CUDR) -o $@ -c lal_vashishta_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/sw.cubin: lal_sw.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_sw.cu
+
+$(OBJ_DIR)/sw_cubin.h: $(OBJ_DIR)/sw.cubin $(OBJ_DIR)/sw.cubin
+ $(BIN2C) -c -n sw $(OBJ_DIR)/sw.cubin > $(OBJ_DIR)/sw_cubin.h
+
+$(OBJ_DIR)/lal_sw.o: $(ALL_H) lal_sw.h lal_sw.cpp $(OBJ_DIR)/sw_cubin.h $(OBJ_DIR)/lal_base_three.o
+ $(CUDR) -o $@ -c lal_sw.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_sw_ext.o: $(ALL_H) lal_sw.h lal_sw_ext.cpp lal_base_three.h
+ $(CUDR) -o $@ -c lal_sw_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/beck.cubin: lal_beck.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_beck.cu
+
+$(OBJ_DIR)/beck_cubin.h: $(OBJ_DIR)/beck.cubin $(OBJ_DIR)/beck.cubin
+ $(BIN2C) -c -n beck $(OBJ_DIR)/beck.cubin > $(OBJ_DIR)/beck_cubin.h
+
+$(OBJ_DIR)/lal_beck.o: $(ALL_H) lal_beck.h lal_beck.cpp $(OBJ_DIR)/beck_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_beck.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_beck_ext.o: $(ALL_H) lal_beck.h lal_beck_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_beck_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/mie.cubin: lal_mie.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_mie.cu
+
+$(OBJ_DIR)/mie_cubin.h: $(OBJ_DIR)/mie.cubin $(OBJ_DIR)/mie.cubin
+ $(BIN2C) -c -n mie $(OBJ_DIR)/mie.cubin > $(OBJ_DIR)/mie_cubin.h
+
+$(OBJ_DIR)/lal_mie.o: $(ALL_H) lal_mie.h lal_mie.cpp $(OBJ_DIR)/mie_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_mie.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_mie_ext.o: $(ALL_H) lal_mie.h lal_mie_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_mie_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/soft.cubin: lal_soft.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_soft.cu
+
+$(OBJ_DIR)/soft_cubin.h: $(OBJ_DIR)/soft.cubin $(OBJ_DIR)/soft.cubin
+ $(BIN2C) -c -n soft $(OBJ_DIR)/soft.cubin > $(OBJ_DIR)/soft_cubin.h
+
+$(OBJ_DIR)/lal_soft.o: $(ALL_H) lal_soft.h lal_soft.cpp $(OBJ_DIR)/soft_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_soft.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_soft_ext.o: $(ALL_H) lal_soft.h lal_soft_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_soft_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_coul_msm.cubin: lal_lj_coul_msm.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_coul_msm.cu
+
+$(OBJ_DIR)/lj_coul_msm_cubin.h: $(OBJ_DIR)/lj_coul_msm.cubin $(OBJ_DIR)/lj_coul_msm.cubin
+ $(BIN2C) -c -n lj_coul_msm $(OBJ_DIR)/lj_coul_msm.cubin > $(OBJ_DIR)/lj_coul_msm_cubin.h
+
+$(OBJ_DIR)/lal_lj_coul_msm.o: $(ALL_H) lal_lj_coul_msm.h lal_lj_coul_msm.cpp $(OBJ_DIR)/lj_coul_msm_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_lj_coul_msm.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_coul_msm_ext.o: $(ALL_H) lal_lj_coul_msm.h lal_lj_coul_msm_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_lj_coul_msm_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_gromacs.cubin: lal_lj_gromacs.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_gromacs.cu
+
+$(OBJ_DIR)/lj_gromacs_cubin.h: $(OBJ_DIR)/lj_gromacs.cubin $(OBJ_DIR)/lj_gromacs.cubin
+ $(BIN2C) -c -n lj_gromacs $(OBJ_DIR)/lj_gromacs.cubin > $(OBJ_DIR)/lj_gromacs_cubin.h
+
+$(OBJ_DIR)/lal_lj_gromacs.o: $(ALL_H) lal_lj_gromacs.h lal_lj_gromacs.cpp $(OBJ_DIR)/lj_gromacs_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_lj_gromacs.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_gromacs_ext.o: $(ALL_H) lal_lj_gromacs.h lal_lj_gromacs_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_lj_gromacs_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/dpd.cubin: lal_dpd.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_dpd.cu
+
+$(OBJ_DIR)/dpd_cubin.h: $(OBJ_DIR)/dpd.cubin $(OBJ_DIR)/dpd.cubin
+ $(BIN2C) -c -n dpd $(OBJ_DIR)/dpd.cubin > $(OBJ_DIR)/dpd_cubin.h
+
+$(OBJ_DIR)/ufm.cubin: lal_ufm.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_ufm.cu
+
+$(OBJ_DIR)/ufm_cubin.h: $(OBJ_DIR)/ufm.cubin $(OBJ_DIR)/ufm.cubin
+ $(BIN2C) -c -n ufm $(OBJ_DIR)/ufm.cubin > $(OBJ_DIR)/ufm_cubin.h
+
+$(OBJ_DIR)/lal_ufm.o: $(ALL_H) lal_ufm.h lal_ufm.cpp $(OBJ_DIR)/ufm_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_ufm.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_ufm_ext.o: $(ALL_H) lal_ufm.h lal_ufm_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_ufm_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_dpd.o: $(ALL_H) lal_dpd.h lal_dpd.cpp $(OBJ_DIR)/dpd_cubin.h $(OBJ_DIR)/lal_base_dpd.o
+ $(CUDR) -o $@ -c lal_dpd.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_dpd_ext.o: $(ALL_H) lal_dpd.h lal_dpd_ext.cpp lal_base_dpd.h
+ $(CUDR) -o $@ -c lal_dpd_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/tersoff.cubin: lal_tersoff.cu lal_precision.h lal_tersoff_extra.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_tersoff.cu
+
+$(OBJ_DIR)/tersoff_cubin.h: $(OBJ_DIR)/tersoff.cubin $(OBJ_DIR)/tersoff.cubin
+ $(BIN2C) -c -n tersoff $(OBJ_DIR)/tersoff.cubin > $(OBJ_DIR)/tersoff_cubin.h
+
+$(OBJ_DIR)/lal_tersoff.o: $(ALL_H) lal_tersoff.h lal_tersoff.cpp $(OBJ_DIR)/tersoff_cubin.h $(OBJ_DIR)/lal_base_three.o
+ $(CUDR) -o $@ -c lal_tersoff.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_tersoff_ext.o: $(ALL_H) lal_tersoff.h lal_tersoff_ext.cpp lal_base_three.h
+ $(CUDR) -o $@ -c lal_tersoff_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/tersoff_zbl.cubin: lal_tersoff_zbl.cu lal_precision.h lal_tersoff_zbl_extra.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_tersoff_zbl.cu
+
+$(OBJ_DIR)/tersoff_zbl_cubin.h: $(OBJ_DIR)/tersoff_zbl.cubin $(OBJ_DIR)/tersoff_zbl.cubin
+ $(BIN2C) -c -n tersoff_zbl $(OBJ_DIR)/tersoff_zbl.cubin > $(OBJ_DIR)/tersoff_zbl_cubin.h
+
+$(OBJ_DIR)/lal_tersoff_zbl.o: $(ALL_H) lal_tersoff_zbl.h lal_tersoff_zbl.cpp $(OBJ_DIR)/tersoff_zbl_cubin.h $(OBJ_DIR)/lal_base_three.o
+ $(CUDR) -o $@ -c lal_tersoff_zbl.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_tersoff_zbl_ext.o: $(ALL_H) lal_tersoff_zbl.h lal_tersoff_zbl_ext.cpp lal_base_three.h
+ $(CUDR) -o $@ -c lal_tersoff_zbl_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/tersoff_mod.cubin: lal_tersoff_mod.cu lal_precision.h lal_tersoff_mod_extra.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_tersoff_mod.cu
+
+$(OBJ_DIR)/tersoff_mod_cubin.h: $(OBJ_DIR)/tersoff_mod.cubin $(OBJ_DIR)/tersoff_mod.cubin
+ $(BIN2C) -c -n tersoff_mod $(OBJ_DIR)/tersoff_mod.cubin > $(OBJ_DIR)/tersoff_mod_cubin.h
+
+$(OBJ_DIR)/lal_tersoff_mod.o: $(ALL_H) lal_tersoff_mod.h lal_tersoff_mod.cpp $(OBJ_DIR)/tersoff_mod_cubin.h $(OBJ_DIR)/lal_base_three.o
+ $(CUDR) -o $@ -c lal_tersoff_mod.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_tersoff_mod_ext.o: $(ALL_H) lal_tersoff_mod.h lal_tersoff_mod_ext.cpp lal_base_three.h
+ $(CUDR) -o $@ -c lal_tersoff_mod_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/coul.cubin: lal_coul.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_coul.cu
+
+$(OBJ_DIR)/coul_cubin.h: $(OBJ_DIR)/coul.cubin $(OBJ_DIR)/coul.cubin
+ $(BIN2C) -c -n coul $(OBJ_DIR)/coul.cubin > $(OBJ_DIR)/coul_cubin.h
+
+$(OBJ_DIR)/lal_coul.o: $(ALL_H) lal_coul.h lal_coul.cpp $(OBJ_DIR)/coul_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_coul.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_coul_ext.o: $(ALL_H) lal_coul.h lal_coul_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_coul_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/coul_debye.cubin: lal_coul_debye.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_coul_debye.cu
+
+$(OBJ_DIR)/coul_debye_cubin.h: $(OBJ_DIR)/coul_debye.cubin $(OBJ_DIR)/coul_debye.cubin
+ $(BIN2C) -c -n coul_debye $(OBJ_DIR)/coul_debye.cubin > $(OBJ_DIR)/coul_debye_cubin.h
+
+$(OBJ_DIR)/lal_coul_debye.o: $(ALL_H) lal_coul_debye.h lal_coul_debye.cpp $(OBJ_DIR)/coul_debye_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_coul_debye.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_coul_debye_ext.o: $(ALL_H) lal_coul_debye.h lal_coul_debye_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_coul_debye_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/zbl.cubin: lal_zbl.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_zbl.cu
+
+$(OBJ_DIR)/zbl_cubin.h: $(OBJ_DIR)/zbl.cubin $(OBJ_DIR)/zbl.cubin
+ $(BIN2C) -c -n zbl $(OBJ_DIR)/zbl.cubin > $(OBJ_DIR)/zbl_cubin.h
+
+$(OBJ_DIR)/lal_zbl.o: $(ALL_H) lal_zbl.h lal_zbl.cpp $(OBJ_DIR)/zbl_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_zbl.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_zbl_ext.o: $(ALL_H) lal_zbl.h lal_zbl_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_zbl_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_cubic.cubin: lal_lj_cubic.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_cubic.cu
+
+$(OBJ_DIR)/lj_cubic_cubin.h: $(OBJ_DIR)/lj_cubic.cubin $(OBJ_DIR)/lj_cubic.cubin
+ $(BIN2C) -c -n lj_cubic $(OBJ_DIR)/lj_cubic.cubin > $(OBJ_DIR)/lj_cubic_cubin.h
+
+$(OBJ_DIR)/lal_lj_cubic.o: $(ALL_H) lal_lj_cubic.h lal_lj_cubic.cpp $(OBJ_DIR)/lj_cubic_cubin.h $(OBJ_DIR)/lal_base_atomic.o
+ $(CUDR) -o $@ -c lal_lj_cubic.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_cubic_ext.o: $(ALL_H) lal_lj_cubic.h lal_lj_cubic_ext.cpp lal_base_atomic.h
+ $(CUDR) -o $@ -c lal_lj_cubic_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/dipole_long_lj.cubin: lal_dipole_long_lj.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_dipole_long_lj.cu
+
+$(OBJ_DIR)/dipole_long_lj_cubin.h: $(OBJ_DIR)/dipole_long_lj.cubin $(OBJ_DIR)/dipole_long_lj.cubin
+ $(BIN2C) -c -n dipole_long_lj $(OBJ_DIR)/dipole_long_lj.cubin > $(OBJ_DIR)/dipole_long_lj_cubin.h
+
+$(OBJ_DIR)/lal_dipole_long_lj.o: $(ALL_H) lal_dipole_long_lj.h lal_dipole_long_lj.cpp $(OBJ_DIR)/dipole_long_lj_cubin.h $(OBJ_DIR)/lal_base_dipole.o
+ $(CUDR) -o $@ -c lal_dipole_long_lj.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_dipole_long_lj_ext.o: $(ALL_H) lal_dipole_long_lj.h lal_dipole_long_lj_ext.cpp lal_base_dipole.h
+ $(CUDR) -o $@ -c lal_dipole_long_lj_ext.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lj_expand_coul_long.cubin: lal_lj_expand_coul_long.cu lal_precision.h lal_preprocessor.h
+ $(CUDA) --fatbin -DNV_KERNEL -o $@ lal_lj_expand_coul_long.cu
+
+$(OBJ_DIR)/lj_expand_coul_long_cubin.h: $(OBJ_DIR)/lj_expand_coul_long.cubin $(OBJ_DIR)/lj_expand_coul_long.cubin
+ $(BIN2C) -c -n lj_expand_coul_long $(OBJ_DIR)/lj_expand_coul_long.cubin > $(OBJ_DIR)/lj_expand_coul_long_cubin.h
+
+$(OBJ_DIR)/lal_lj_expand_coul_long.o: $(ALL_H) lal_lj_expand_coul_long.h lal_lj_expand_coul_long.cpp $(OBJ_DIR)/lj_expand_coul_long_cubin.h $(OBJ_DIR)/lal_base_charge.o
+ $(CUDR) -o $@ -c lal_lj_expand_coul_long.cpp -I$(OBJ_DIR)
+
+$(OBJ_DIR)/lal_lj_expand_coul_long_ext.o: $(ALL_H) lal_lj_expand_coul_long.h lal_lj_expand_coul_long_ext.cpp lal_base_charge.h
+ $(CUDR) -o $@ -c lal_lj_expand_coul_long_ext.cpp -I$(OBJ_DIR)
+
+$(BIN_DIR)/nvc_get_devices: ./geryon/ucl_get_devices.cpp $(NVD_H)
+ $(CUDR) -o $@ ./geryon/ucl_get_devices.cpp -DUCL_CUDADR $(CUDA_LIB) -lcuda
+
+$(GPU_LIB): $(OBJS) $(CUDPP)
+ $(AR) -crusv $(GPU_LIB) $(OBJS) $(CUDPP)
+ @cp $(EXTRAMAKE) Makefile.lammps
+
+clean:
+ -rm -f $(EXECS) $(GPU_LIB) $(OBJS) $(CUDPP) $(CBNS) *.linkinfo
+
+veryclean: clean
+ -rm -rf *~ *.linkinfo
+
+cleanlib:
+ -rm -f $(EXECS) $(GPU_LIB) $(OBJS) $(CBNS) *.linkinfo
diff --git a/lib/gpu/lal_dipole_long_lj.cpp b/lib/gpu/lal_dipole_long_lj.cpp
new file mode 100644
index 0000000000..9cad926a2d
--- /dev/null
+++ b/lib/gpu/lal_dipole_long_lj.cpp
@@ -0,0 +1,173 @@
+/***************************************************************************
+ dipole_lj.cpp
+ -------------------
+ Trung Dac Nguyen (ORNL)
+
+ Class for acceleration of the dipole/cut pair style.
+
+ __________________________________________________________________________
+ This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
+ __________________________________________________________________________
+
+ begin :
+ email : nguyentd@ornl.gov
+ ***************************************************************************/
+
+#ifdef USE_OPENCL
+#include "dipole_long_lj_cl.h"
+#elif defined(USE_CUDART)
+const char *dipole_long_lj=0;
+#else
+#include "dipole_long_lj_cubin.h"
+#endif
+
+#include "lal_dipole_long_lj.h"
+#include <cassert>
+using namespace LAMMPS_AL;
+#define DipoleLongLJT DipoleLongLJ<numtyp, acctyp>
+
+extern Device<PRECISION,ACC_PRECISION> device;
+
+template <class numtyp, class acctyp>
+DipoleLongLJT::DipoleLongLJ() : BaseDipole<numtyp,acctyp>(),
+ _allocated(false) {
+}
+
+template <class numtyp, class acctyp>
+DipoleLongLJT::~DipoleLongLJ() {
+ clear();
+}
+
+template <class numtyp, class acctyp>
+int DipoleLongLJT::bytes_per_atom(const int max_nbors) const {
+ return this->bytes_per_atom_atomic(max_nbors);
+}
+
+template <class numtyp, class acctyp>
+int DipoleLongLJT::init(const int ntypes,
+ double **host_cutsq, 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 double g_ewald) {
+ int success;
+ success=this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,gpu_split,
+ _screen,dipole_long_lj,"k_dipole_long_lj");
+ if (success!=0)
+ return success;
+
+ // If atom type constants fit in shared memory use fast kernel
+ int lj_types=ntypes;
+ shared_types=false;
+ int max_shared_types=this->device->max_shared_types();
+ 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<numtyp> host_write(lj_types*lj_types*32,*(this->ucl_device),
+ UCL_WRITE_ONLY);
+
+ for (int i=0; i<lj_types*lj_types; i++)
+ host_write[i]=0.0;
+
+ lj1.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
+ this->atom->type_pack4(ntypes,lj_types,lj1,host_write,host_lj1,host_lj2,
+ host_cut_ljsq);
+
+ lj3.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
+ this->atom->type_pack4(ntypes,lj_types,lj3,host_write,host_lj3,host_lj4,
+ host_offset);
+
+ cutsq.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
+ this->atom->type_pack1(ntypes,lj_types,cutsq,host_write,host_cutsq);
+
+ 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;
+ _g_ewald=g_ewald;
+
+ _allocated=true;
+ this->_max_bytes=lj1.row_bytes()+lj3.row_bytes()+cutsq.row_bytes()+
+ sp_lj.row_bytes();
+ return 0;
+}
+
+template <class numtyp, class acctyp>
+void DipoleLongLJT::clear() {
+ if (!_allocated)
+ return;
+ _allocated=false;
+
+ lj1.clear();
+ lj3.clear();
+ cutsq.clear();
+ sp_lj.clear();
+ this->clear_atomic();
+}
+
+template <class numtyp, class acctyp>
+double DipoleLongLJT::host_memory_usage() const {
+ return this->host_memory_usage_atomic()+sizeof(DipoleLongLJ<numtyp,acctyp>);
+}
+
+// ---------------------------------------------------------------------------
+// Calculate energies, forces, and torques
+// ---------------------------------------------------------------------------
+template <class numtyp, class acctyp>
+void DipoleLongLJT::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<int>(ceil(static_cast<double>(this->ans->inum())/
+ (BX/this->_threads_per_atom)));
+
+ int ainum=this->ans->inum();
+ 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->x, &lj1, &lj3, &sp_lj,
+ &this->nbor->dev_nbor,
+ &this->_nbor_data->begin(),
+ &this->ans->force, &this->ans->engv, &eflag, &vflag,
+ &ainum, &nbor_pitch, &this->atom->q,
+ &this->atom->quat, &cutsq, &_cut_coulsq,
+ &_qqrd2e, &this->_threads_per_atom);
+ } else {
+ this->k_pair.set_size(GX,BX);
+ this->k_pair.run(&this->atom->x, &lj1, &lj3,
+ &_lj_types, &sp_lj, &this->nbor->dev_nbor,
+ &this->_nbor_data->begin(), &this->ans->force,
+ &this->ans->engv, &eflag, &vflag, &ainum,
+ &nbor_pitch, &this->atom->q,
+ &this->atom->quat, &cutsq, &_cut_coulsq,
+ &_qqrd2e, &this->_threads_per_atom);
+ }
+ this->time_pair.stop();
+}
+
+template class DipoleLongLJ<PRECISION,ACC_PRECISION>;
diff --git a/lib/gpu/lal_dipole_long_lj.cu b/lib/gpu/lal_dipole_long_lj.cu
new file mode 100644
index 0000000000..4c1cf9352d
--- /dev/null
+++ b/lib/gpu/lal_dipole_long_lj.cu
@@ -0,0 +1,640 @@
+// **************************************************************************
+// dipole_lj.cu
+// -------------------
+// Trung Dac Nguyen (ORNL)
+//
+// Device code for acceleration of the dipole/cut pair style
+//
+// __________________________________________________________________________
+// This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
+// __________________________________________________________________________
+//
+// begin :
+// email : nguyentd@ornl.gov
+// ***************************************************************************/
+
+#ifdef NV_KERNEL
+#include "lal_aux_fun1.h"
+#ifndef _DOUBLE_DOUBLE
+texture<float4> pos_tex;
+texture<float> q_tex;
+texture<float4> mu_tex;
+#else
+texture<int4,1> pos_tex;
+texture<int2> q_tex;
+texture<int4,1> mu_tex;
+#endif
+
+#else
+#define pos_tex x_
+#define q_tex q_
+#define mu_tex mu_
+#endif
+
+#if (ARCH < 300)
+
+#define store_answers_tq(f, tor, energy, ecoul, virial, ii, inum, tid, \
+ t_per_atom, offset, eflag, vflag, ans, engv) \
+ if (t_per_atom>1) { \
+ __local acctyp red_acc[8][BLOCK_PAIR]; \
+ red_acc[0][tid]=f.x; \
+ red_acc[1][tid]=f.y; \
+ red_acc[2][tid]=f.z; \
+ red_acc[3][tid]=tor.x; \
+ red_acc[4][tid]=tor.y; \
+ red_acc[5][tid]=tor.z; \
+ for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \
+ if (offset < s) { \
+ for (int r=0; r<6; r++) \
+ red_acc[r][tid] += red_acc[r][tid+s]; \
+ } \
+ } \
+ f.x=red_acc[0][tid]; \
+ f.y=red_acc[1][tid]; \
+ f.z=red_acc[2][tid]; \
+ tor.x=red_acc[3][tid]; \
+ tor.y=red_acc[4][tid]; \
+ tor.z=red_acc[5][tid]; \
+ if (eflag>0 || vflag>0) { \
+ for (int r=0; r<6; r++) \
+ red_acc[r][tid]=virial[r]; \
+ red_acc[6][tid]=energy; \
+ red_acc[7][tid]=ecoul; \
+ for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \
+ if (offset < s) { \
+ for (int r=0; r<8; r++) \
+ red_acc[r][tid] += red_acc[r][tid+s]; \
+ } \
+ } \
+ for (int r=0; r<6; r++) \
+ virial[r]=red_acc[r][tid]; \
+ energy=red_acc[6][tid]; \
+ ecoul=red_acc[7][tid]; \
+ } \
+ } \
+ if (offset==0) { \
+ int ei=ii; \
+ if (eflag>0) { \
+ engv[ei]=energy*(acctyp)0.5; \
+ ei+=inum; \
+ engv[ei]=e_coul*(acctyp)0.5; \
+ ei+=inum; \
+ } \
+ if (vflag>0) { \
+ for (int i=0; i<6; i++) { \
+ engv[ei]=virial[i]*(acctyp)0.5; \
+ ei+=inum; \
+ } \
+ } \
+ ans[ii]=f; \
+ ans[ii+inum]=tor; \
+ }
+
+#else
+
+#define store_answers_tq(f, tor, energy, e_coul, virial, ii, inum, tid, \
+ t_per_atom, offset, eflag, vflag, ans, engv) \
+ if (t_per_atom>1) { \
+ for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \
+ f.x += shfl_xor(f.x, s, t_per_atom); \
+ f.y += shfl_xor(f.y, s, t_per_atom); \
+ f.z += shfl_xor(f.z, s, t_per_atom); \
+ tor.x += shfl_xor(tor.x, s, t_per_atom); \
+ tor.y += shfl_xor(tor.y, s, t_per_atom); \
+ tor.z += shfl_xor(tor.z, s, t_per_atom); \
+ energy += shfl_xor(energy, s, t_per_atom); \
+ e_coul += shfl_xor(e_coul, s, t_per_atom); \
+ } \
+ if (vflag>0) { \
+ for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \
+ for (int r=0; r<6; r++) \
+ virial[r] += shfl_xor(virial[r], s, t_per_atom); \
+ } \
+ } \
+ } \
+ if (offset==0) { \
+ int ei=ii; \
+ if (eflag>0) { \
+ engv[ei]=energy*(acctyp)0.5; \
+ ei+=inum; \
+ engv[ei]=e_coul*(acctyp)0.5; \
+ ei+=inum; \
+ } \
+ if (vflag>0) { \
+ for (int i=0; i<6; i++) { \
+ engv[ei]=virial[i]*(acctyp)0.5; \
+ ei+=inum; \
+ } \
+ } \
+ ans[ii]=f; \
+ ans[ii+inum]=tor; \
+ }
+
+#endif
+
+#define MY_PIS (acctyp)1.77245385090551602729
+
+__kernel void k_dipole_long_lj(const __global numtyp4 *restrict x_,
+ const __global numtyp4 *restrict lj1,
+ const __global numtyp4 *restrict lj3,
+ const int lj_types,
+ const __global numtyp *restrict sp_lj_in,
+ const __global int *dev_nbor,
+ const __global int *dev_packed,
+ __global acctyp4 *restrict ans,
+ __global acctyp *restrict engv,
+ const int eflag, const int vflag, const int inum,
+ const int nbor_pitch,
+ const __global numtyp *restrict q_,
+ const __global numtyp4 *restrict mu_,
+ const __global numtyp *restrict cutsq,
+ const numtyp cut_coulsq, const numtyp qqrd2e,
+ const numtyp g_ewald, const int t_per_atom) {
+ int tid, ii, offset;
+ atom_info(t_per_atom,ii,tid,offset);
+
+ __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];
+
+ acctyp energy=(acctyp)0;
+ acctyp e_coul=(acctyp)0;
+ acctyp4 f;
+ f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
+ acctyp4 tor;
+ tor.x=(acctyp)0;
+ tor.y=(acctyp)0;
+ tor.z=(acctyp)0;
+ acctyp virial[6];
+ for (int i=0; i<6; i++)
+ virial[i]=(acctyp)0;
+
+ numtyp pre1 = numtyp(2.0) * g_ewald / MY_PIS;
+ numtyp pre2 = numtyp(4.0) * (g_ewald*g_ewald*g_ewald) / MY_PIS;
+ numtyp pre3 = numtyp(8.0) * (g_ewald*g_ewald*g_ewald*g_ewald*g_ewald) / MY_PIS;
+
+ if (ii<inum) {
+ int nbor, nbor_end;
+ int i, numj;
+ __local int n_stride;
+ nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
+ n_stride,nbor_end,nbor);
+
+ numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
+ numtyp qtmp; fetch(qtmp,i,q_tex);
+ numtyp4 mui; fetch4(mui,i,mu_tex); //mu_[i];
+ int itype=ix.w;
+
+ for ( ; nbor<nbor_end; nbor+=n_stride) {
+ int j=dev_packed[nbor];
+
+ numtyp factor_lj, factor_coul;
+ factor_lj = sp_lj[sbmask(j)];
+ factor_coul = sp_lj[sbmask(j)+4];
+ j &= NEIGHMASK;
+
+ numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
+ numtyp qj; fetch(qj,j,q_tex);
+ numtyp4 muj; fetch4(muj,j,mu_tex); //mu_[j];
+ int jtype=jx.w;
+
+ // Compute r12
+ numtyp delx = ix.x-jx.x;
+ numtyp dely = ix.y-jx.y;
+ numtyp delz = ix.z-jx.z;
+ numtyp rsq = delx*delx+dely*dely+delz*delz;
+
+ int mtype=itype*lj_types+jtype;
+ if (rsq<cutsq[mtype]) {
+ numtyp r2inv=ucl_recip(rsq);
+ numtyp force_lj,rinv,r6inv;
+ numtyp pdotp, pidotr, pjdotr, _erfc;
+ numtyp g0,g1,g2,b0,b1,b2,b3,d0,d1,d2,d3;
+ numtyp zdix,zdiy,zdiz,zdjx,zdjy,zdjz,zaix,zaiy,zaiz,zajx,zajy,zajz;
+ numtyp g0b1_g1b2_g2b3,g0d1_g1d2_g2d3,facm1;
+ numtyp fdx,fdy,fdz,fax,fay,faz;
+ acctyp4 forcecoul, ticoul;
+ acctyp4 force;
+
+ forcecoul.x = forcecoul.y = forcecoul.z = (acctyp)0;
+ ticoul.x = ticoul.y = ticoul.z = (acctyp)0;
+
+ if (rsq < lj1[mtype].z) {
+ r6inv = r2inv*r2inv*r2inv;
+ force_lj = factor_lj*r6inv*(lj1[mtype].x*r6inv-lj1[mtype].y)*r2inv;
+ } else force_lj = (numtyp)0.0;
+
+ if (rsq < cut_coulsq) {
+ rinv = ucl_rsqrt(rsq);
+ numtyp r = ucl_rsqrt(r2inv);
+ numtyp grij = g_ewald * r;
+ numtyp expm2 = ucl_exp(-grij*grij);
+ numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij);
+ _erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+
+ pdotp = mui.x*muj.x + mui.y*muj.y + mui.z*muj.z;
+ pidotr = mui.x*delx + mui.y*dely + mui.z*delz;
+ pjdotr = muj.x*delx + muj.y*dely + muj.z*delz;
+
+ g0 = qtmp*qj;
+ g1 = qtmp*pjdotr - qj*pidotr + pdotp;
+ g2 = -pidotr*pjdotr;
+
+ if (factor_coul > (numtyp)0.0) {
+ b0 = _erfc * rinv;
+ b1 = (b0 + pre1*expm2) * r2inv;
+ b2 = ((numtyp)3.0*b1 + pre2*expm2) * r2inv;
+ b3 = ((numtyp)5.0*b2 + pre3*expm2) * r2inv;
+
+ g0b1_g1b2_g2b3 = g0*b1 + g1*b2 + g2*b3;
+ fdx = delx * g0b1_g1b2_g2b3 -
+ b1 * (qtmp*muj.x - qj*mui.x) +
+ b2 * (pjdotr*mui.x + pidotr*muj.x);
+ fdy = dely * g0b1_g1b2_g2b3 -
+ b1 * (qtmp*muj.y - qj*mui.y) +
+ b2 * (pjdotr*mui.y + pidotr*muj.y);
+ fdz = delz * g0b1_g1b2_g2b3 -
+ b1 * (qtmp*muj.z - qj*mui.z) +
+ b2 * (pjdotr*mui.z + pidotr*muj.z);
+
+ zdix = delx * (qj*b1 + b2*pjdotr) - b1*muj.x;
+ zdiy = dely * (qj*b1 + b2*pjdotr) - b1*muj.y;
+ zdiz = delz * (qj*b1 + b2*pjdotr) - b1*muj.z;
+ zdjx = delx * (-qtmp*b1 + b2*pidotr) - b1*mui.x;
+ zdjy = dely * (-qtmp*b1 + b2*pidotr) - b1*mui.y;
+ zdjz = delz * (-qtmp*b1 + b2*pidotr) - b1*mui.z;
+
+ if (factor_coul < (numtyp)1.0) {
+ fdx *= factor_coul;
+ fdy *= factor_coul;
+ fdz *= factor_coul;
+ zdix *= factor_coul;
+ zdiy *= factor_coul;
+ zdiz *= factor_coul;
+ zdjx *= factor_coul;
+ zdjy *= factor_coul;
+ zdjz *= factor_coul;
+ }
+ } else {
+ fdx = fdy = fdz = (numtyp)0.0;
+ zdix = zdiy = zdiz = (numtyp)0.0;
+ zdjx = zdjy = zdjz = (numtyp)0.0;
+ }
+
+ if (factor_coul < (numtyp)1.0) {
+ d0 = (_erfc - (numtyp)1.0) * rinv;
+ d1 = (d0 + pre1*expm2) * r2inv;
+ d2 = ((numtyp)3.0*d1 + pre2*expm2) * r2inv;
+ d3 = ((numtyp)5.0*d2 + pre3*expm2) * r2inv;
+
+ g0d1_g1d2_g2d3 = g0*d1 + g1*d2 + g2*d3;
+ fax = delx * g0d1_g1d2_g2d3 -
+ d1 * (qtmp*muj.x - qj*mui.x) +
+ d2 * (pjdotr*mui.x + pidotr*muj.x);
+ fay = dely * g0d1_g1d2_g2d3 -
+ d1 * (qtmp*muj.y - qj*mui.y) +
+ d2 * (pjdotr*mui.y + pidotr*muj.y);
+ faz = delz * g0d1_g1d2_g2d3 -
+ d1 * (qtmp*muj.z - qj*mui.z) +
+ d2 * (pjdotr*mui.z + pidotr*muj.z);
+
+ zaix = delx * (qj*d1 + d2*pjdotr) - d1*muj.x;
+ zaiy = dely * (qj*d1 + d2*pjdotr) - d1*muj.y;
+ zaiz = delz * (qj*d1 + d2*pjdotr) - d1*muj.z;
+ zajx = delx * (-qtmp*d1 + d2*pidotr) - d1*mui.x;
+ zajy = dely * (-qtmp*d1 + d2*pidotr) - d1*mui.y;
+ zajz = delz * (-qtmp*d1 + d2*pidotr) - d1*mui.z;
+
+ if (factor_coul > (numtyp)0.0) {
+ facm1 = (numtyp)1.0 - factor_coul;
+ fax *= facm1;
+ fay *= facm1;
+ faz *= facm1;
+ zaix *= facm1;
+ zaiy *= facm1;
+ zaiz *= facm1;
+ zajx *= facm1;
+ zajy *= facm1;
+ zajz *= facm1;
+ }
+ } else {
+ fax = fay = faz = (numtyp)0.0;
+ zaix = zaiy = zaiz = (numtyp)0.0;
+ zajx = zajy = zajz = (numtyp)0.0;
+ }
+
+ forcecoul.x = fdx + fax;
+ forcecoul.y = fdy + fay;
+ forcecoul.z = fdz + faz;
+
+ ticoul.x = mui.y*(zdiz + zaiz) - mui.z*(zdiy + zaiy);
+ ticoul.y = mui.z*(zdix + zaix) - mui.x*(zdiz + zaiz);
+ ticoul.z = mui.x*(zdiy + zaiy) - mui.y*(zdix + zaix);
+
+ } else {
+ forcecoul.x = forcecoul.y = forcecoul.z = (numtyp)0.0;
+ ticoul.x = ticoul.y = ticoul.z = (numtyp)0.0;
+ }
+
+ force.x = qqrd2e*forcecoul.x + delx*force_lj;
+ force.y = qqrd2e*forcecoul.y + dely*force_lj;
+ force.z = qqrd2e*forcecoul.z + delz*force_lj;
+ f.x+=force.x;
+ f.y+=force.y;
+ f.z+=force.z;
+ tor.x+=qqrd2e*ticoul.x;
+ tor.y+=qqrd2e*ticoul.y;
+ tor.z+=qqrd2e*ticoul.z;
+
+ if (eflag>0) {
+ acctyp e = (acctyp)0.0;
+ if (rsq < lj1[mtype].w && factor_coul > (numtyp)0.0) {
+ e = qqrd2e*(b0*g0 + b1*g1 + b2*g2);
+ if (factor_coul < (numtyp)1.0) {
+ e_coul *= factor_coul;
+ e_coul += ((numtyp)1.0-factor_coul) * qqrd2e * (d0*g0 + d1*g1 + d2*g2);
+ }
+ } else e = (acctyp)0.0;
+ e_coul += e;
+
+ if (rsq < lj1[mtype].z) {
+ e=r6inv*(lj3[mtype].x*r6inv-lj3[mtype].y);
+ energy+=factor_lj*(e-lj3[mtype].z);
+ }
+ }
+ if (vflag>0) {
+ virial[0] += delx*force.x;
+ virial[1] += dely*force.y;
+ virial[2] += delz*force.z;
+ virial[3] += delx*force.y;
+ virial[4] += delx*force.z;
+ virial[5] += dely*force.z;
+ }
+ }
+
+ } // for nbor
+ store_answers_tq(f,tor,energy,e_coul,virial,ii,inum,tid,t_per_atom,offset,eflag,
+ vflag,ans,engv);
+ } // if ii
+}
+
+__kernel void k_dipole_long_lj_fast(const __global numtyp4 *restrict x_,
+ const __global numtyp4 *restrict lj1_in,
+ const __global numtyp4 *restrict lj3_in,
+ const __global numtyp *restrict sp_lj_in,
+ const __global int *dev_nbor,
+ const __global int *dev_packed,
+ __global acctyp4 *restrict ans,
+ __global acctyp *restrict engv,
+ const int eflag, const int vflag, const int inum,
+ const int nbor_pitch,
+ const __global numtyp *restrict q_,
+ const __global numtyp4 *restrict mu_,
+ const __global numtyp *restrict _cutsq,
+ const numtyp cut_coulsq, const numtyp qqrd2e,
+ const numtyp g_ewald, const int t_per_atom) {
+ int tid, ii, offset;
+ atom_info(t_per_atom,ii,tid,offset);
+
+ __local numtyp4 lj1[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
+ __local numtyp4 lj3[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
+ __local numtyp cutsq[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
+ __local numtyp sp_lj[8];
+ if (tid<8)
+ sp_lj[tid]=sp_lj_in[tid];
+ if (tid<MAX_SHARED_TYPES*MAX_SHARED_TYPES) {
+ lj1[tid]=lj1_in[tid];
+ cutsq[tid]=_cutsq[tid];
+ if (eflag>0)
+ lj3[tid]=lj3_in[tid];
+ }
+
+ acctyp energy=(acctyp)0;
+ acctyp e_coul=(acctyp)0;
+ acctyp4 f;
+ f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
+ acctyp4 tor;
+ tor.x=(acctyp)0;
+ tor.y=(acctyp)0;
+ tor.z=(acctyp)0;
+ acctyp virial[6];
+ for (int i=0; i<6; i++)
+ virial[i]=(acctyp)0;
+
+ __syncthreads();
+
+ numtyp pre1 = numtyp(2.0) * g_ewald / MY_PIS;
+ numtyp pre2 = numtyp(4.0) * (g_ewald*g_ewald*g_ewald) / MY_PIS;
+ numtyp pre3 = numtyp(8.0) * (g_ewald*g_ewald*g_ewald*g_ewald*g_ewald) / MY_PIS;
+
+ if (ii<inum) {
+ int nbor, nbor_end;
+ int i, numj;
+ __local int n_stride;
+ nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
+ n_stride,nbor_end,nbor);
+
+ numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
+ numtyp qtmp; fetch(qtmp,i,q_tex);
+ numtyp4 mui; fetch4(mui,i,mu_tex); //mu_[i];
+ int iw=ix.w;
+ int itype=fast_mul((int)MAX_SHARED_TYPES,iw);
+
+ for ( ; nbor<nbor_end; nbor+=n_stride) {
+ int j=dev_packed[nbor];
+
+ numtyp factor_lj, factor_coul;
+ factor_lj = sp_lj[sbmask(j)];
+ factor_coul = sp_lj[sbmask(j)+4];
+ j &= NEIGHMASK;
+
+ numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
+ numtyp qj; fetch(qj,j,q_tex);
+ numtyp4 muj; fetch4(muj,j,mu_tex); //mu_[j];
+ int mtype=itype+jx.w;
+
+ // Compute r12
+ numtyp delx = ix.x-jx.x;
+ numtyp dely = ix.y-jx.y;
+ numtyp delz = ix.z-jx.z;
+ numtyp rsq = delx*delx+dely*dely+delz*delz;
+
+ if (rsq<cutsq[mtype]) {
+ numtyp r2inv=ucl_recip(rsq);
+ numtyp force_lj,rinv,r6inv;
+ numtyp pdotp, pidotr, pjdotr, _erfc;
+ numtyp g0,g1,g2,b0,b1,b2,b3,d0,d1,d2,d3;
+ numtyp zdix,zdiy,zdiz,zdjx,zdjy,zdjz,zaix,zaiy,zaiz,zajx,zajy,zajz;
+ numtyp g0b1_g1b2_g2b3,g0d1_g1d2_g2d3,facm1;
+ numtyp fdx,fdy,fdz,fax,fay,faz;
+ acctyp4 forcecoul, ticoul;
+ acctyp4 force;
+
+ forcecoul.x = forcecoul.y = forcecoul.z = (acctyp)0;
+ ticoul.x = ticoul.y = ticoul.z = (acctyp)0;
+
+ if (rsq < lj1[mtype].z) {
+ r6inv = r2inv*r2inv*r2inv;
+ force_lj = factor_lj*r6inv*(lj1[mtype].x*r6inv-lj1[mtype].y)*r2inv;
+ } else force_lj = (numtyp)0.0;
+
+ if (rsq < cut_coulsq) {
+ rinv = ucl_rsqrt(rsq);
+ numtyp r = ucl_rsqrt(r2inv);
+ numtyp grij = g_ewald * r;
+ numtyp expm2 = ucl_exp(-grij*grij);
+ numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij);
+ _erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+
+ pdotp = mui.x*muj.x + mui.y*muj.y + mui.z*muj.z;
+ pidotr = mui.x*delx + mui.y*dely + mui.z*delz;
+ pjdotr = muj.x*delx + muj.y*dely + muj.z*delz;
+
+ g0 = qtmp*qj;
+ g1 = qtmp*pjdotr - qj*pidotr + pdotp;
+ g2 = -pidotr*pjdotr;
+
+ if (factor_coul > (numtyp)0.0) {
+ b0 = _erfc * rinv;
+ b1 = (b0 + pre1*expm2) * r2inv;
+ b2 = ((numtyp)3.0*b1 + pre2*expm2) * r2inv;
+ b3 = ((numtyp)5.0*b2 + pre3*expm2) * r2inv;
+
+ g0b1_g1b2_g2b3 = g0*b1 + g1*b2 + g2*b3;
+ fdx = delx * g0b1_g1b2_g2b3 -
+ b1 * (qtmp*muj.x - qj*mui.x) +
+ b2 * (pjdotr*mui.x + pidotr*muj.x);
+ fdy = dely * g0b1_g1b2_g2b3 -
+ b1 * (qtmp*muj.y - qj*mui.y) +
+ b2 * (pjdotr*mui.y + pidotr*muj.y);
+ fdz = delz * g0b1_g1b2_g2b3 -
+ b1 * (qtmp*muj.z - qj*mui.z) +
+ b2 * (pjdotr*mui.z + pidotr*muj.z);
+
+ zdix = delx * (qj*b1 + b2*pjdotr) - b1*muj.x;
+ zdiy = dely * (qj*b1 + b2*pjdotr) - b1*muj.y;
+ zdiz = delz * (qj*b1 + b2*pjdotr) - b1*muj.z;
+ zdjx = delx * (-qtmp*b1 + b2*pidotr) - b1*mui.x;
+ zdjy = dely * (-qtmp*b1 + b2*pidotr) - b1*mui.y;
+ zdjz = delz * (-qtmp*b1 + b2*pidotr) - b1*mui.z;
+
+ if (factor_coul < (numtyp)1.0) {
+ fdx *= factor_coul;
+ fdy *= factor_coul;
+ fdz *= factor_coul;
+ zdix *= factor_coul;
+ zdiy *= factor_coul;
+ zdiz *= factor_coul;
+ zdjx *= factor_coul;
+ zdjy *= factor_coul;
+ zdjz *= factor_coul;
+ }
+ } else {
+ fdx = fdy = fdz = (numtyp)0.0;
+ zdix = zdiy = zdiz = (numtyp)0.0;
+ zdjx = zdjy = zdjz = (numtyp)0.0;
+ }
+
+ if (factor_coul < (numtyp)1.0) {
+ d0 = (_erfc - (numtyp)1.0) * rinv;
+ d1 = (d0 + pre1*expm2) * r2inv;
+ d2 = ((numtyp)3.0*d1 + pre2*expm2) * r2inv;
+ d3 = ((numtyp)5.0*d2 + pre3*expm2) * r2inv;
+
+ g0d1_g1d2_g2d3 = g0*d1 + g1*d2 + g2*d3;
+ fax = delx * g0d1_g1d2_g2d3 -
+ d1 * (qtmp*muj.x - qj*mui.x) +
+ d2 * (pjdotr*mui.x + pidotr*muj.x);
+ fay = dely * g0d1_g1d2_g2d3 -
+ d1 * (qtmp*muj.y - qj*mui.y) +
+ d2 * (pjdotr*mui.y + pidotr*muj.y);
+ faz = delz * g0d1_g1d2_g2d3 -
+ d1 * (qtmp*muj.z - qj*mui.z) +
+ d2 * (pjdotr*mui.z + pidotr*muj.z);
+
+ zaix = delx * (qj*d1 + d2*pjdotr) - d1*muj.x;
+ zaiy = dely * (qj*d1 + d2*pjdotr) - d1*muj.y;
+ zaiz = delz * (qj*d1 + d2*pjdotr) - d1*muj.z;
+ zajx = delx * (-qtmp*d1 + d2*pidotr) - d1*mui.x;
+ zajy = dely * (-qtmp*d1 + d2*pidotr) - d1*mui.y;
+ zajz = delz * (-qtmp*d1 + d2*pidotr) - d1*mui.z;
+
+ if (factor_coul > (numtyp)0.0) {
+ facm1 = (numtyp)1.0 - factor_coul;
+ fax *= facm1;
+ fay *= facm1;
+ faz *= facm1;
+ zaix *= facm1;
+ zaiy *= facm1;
+ zaiz *= facm1;
+ zajx *= facm1;
+ zajy *= facm1;
+ zajz *= facm1;
+ }
+ } else {
+ fax = fay = faz = (numtyp)0.0;
+ zaix = zaiy = zaiz = (numtyp)0.0;
+ zajx = zajy = zajz = (numtyp)0.0;
+ }
+
+ forcecoul.x = fdx + fax;
+ forcecoul.y = fdy + fay;
+ forcecoul.z = fdz + faz;
+
+ ticoul.x = mui.y*(zdiz + zaiz) - mui.z*(zdiy + zaiy);
+ ticoul.y = mui.z*(zdix + zaix) - mui.x*(zdiz + zaiz);
+ ticoul.z = mui.x*(zdiy + zaiy) - mui.y*(zdix + zaix);
+
+ } else {
+ forcecoul.x = forcecoul.y = forcecoul.z = (numtyp)0.0;
+ ticoul.x = ticoul.y = ticoul.z = (numtyp)0.0;
+ }
+
+ force.x = qqrd2e*forcecoul.x + delx*force_lj;
+ force.y = qqrd2e*forcecoul.y + dely*force_lj;
+ force.z = qqrd2e*forcecoul.z + delz*force_lj;
+ f.x+=force.x;
+ f.y+=force.y;
+ f.z+=force.z;
+ tor.x+=qqrd2e*ticoul.x;
+ tor.y+=qqrd2e*ticoul.y;
+ tor.z+=qqrd2e*ticoul.z;
+
+ if (eflag>0) {
+ acctyp e = (acctyp)0.0;
+ if (rsq < lj1[mtype].w && factor_coul > (numtyp)0.0) {
+ e = qqrd2e*(b0*g0 + b1*g1 + b2*g2);
+ if (factor_coul < (numtyp)1.0) {
+ e_coul *= factor_coul;
+ e_coul += ((numtyp)1.0-factor_coul) * qqrd2e * (d0*g0 + d1*g1 + d2*g2);
+ }
+ } else e = (acctyp)0.0;
+ e_coul += e;
+
+ if (rsq < lj1[mtype].z) {
+ e=r6inv*(lj3[mtype].x*r6inv-lj3[mtype].y);
+ energy+=factor_lj*(e-lj3[mtype].z);
+ }
+ }
+ if (vflag>0) {
+ virial[0] += delx*force.x;
+ virial[1] += dely*force.y;
+ virial[2] += delz*force.z;
+ virial[3] += delx*force.y;
+ virial[4] += delx*force.z;
+ virial[5] += dely*force.z;
+ }
+ }
+
+ } // for nbor
+ store_answers_tq(f,tor,energy,e_coul,virial,ii,inum,tid,t_per_atom,offset,eflag,
+ vflag,ans,engv);
+ } // if ii
+}
+
diff --git a/lib/gpu/lal_dipole_long_lj.h b/lib/gpu/lal_dipole_long_lj.h
new file mode 100644
index 0000000000..1381e24326
--- /dev/null
+++ b/lib/gpu/lal_dipole_long_lj.h
@@ -0,0 +1,85 @@
+/***************************************************************************
+ dipole_long_lj.h
+ -------------------
+ Trung Dac Nguyen (Northwestern)
+
+ Class for acceleration of the lj/cut/dipole/long pair style.
+
+ __________________________________________________________________________
+ This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
+ __________________________________________________________________________
+
+ begin :
+ email : ndactrung@gmail.com
+ ***************************************************************************/
+
+#ifndef LAL_DIPOLE_LONG_LJ_H
+#define LAL_DIPOLE_LONG_LJ_H
+
+#include "lal_base_dipole.h"
+
+namespace LAMMPS_AL {
+
+template <class numtyp, class acctyp>
+class DipoleLongLJ : public BaseDipole<numtyp, acctyp> {
+ public:
+ DipoleLongLJ();
+ ~DipoleLongLJ();
+
+ /// 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
+ *
+ * Returns:
+ * - 0 if successfull
+ * - -1 if fix gpu not found
+ * - -3 if there is an out of memory error
+ * - -4 if the GPU library was not compiled for GPU
+ * - -5 Double precision is not supported on card **/
+ int init(const int ntypes, double **host_cutsq, 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 double g_ewald);
+
+ /// 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 = lj1, lj1.y = lj2, lj1.z = cutsq_vdw, lj1.w = cutsq_coul
+ UCL_D_Vec<numtyp4> lj1;
+ /// lj3.x = lj3, lj3.y = lj4, lj3.z = offset
+ UCL_D_Vec<numtyp4> lj3;
+ /// cutsq
+ UCL_D_Vec<numtyp> cutsq;
+ /// Special LJ values [0-3] and Special Coul values [4-7]
+ UCL_D_Vec<numtyp> 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, _g_ewald;
+
+ private:
+ bool _allocated;
+ void loop(const bool _eflag, const bool _vflag);
+};
+
+}
+
+#endif
diff --git a/lib/gpu/lal_dipole_long_lj_ext.cpp b/lib/gpu/lal_dipole_long_lj_ext.cpp
new file mode 100644
index 0000000000..b58bcea5d0
--- /dev/null
+++ b/lib/gpu/lal_dipole_long_lj_ext.cpp
@@ -0,0 +1,129 @@
+/***************************************************************************
+ dipole_long_lj_ext.cpp
+ -------------------
+ Trung Dac Nguyen (ORNL)
+
+ Functions for LAMMPS access to dipole/cut acceleration routines.
+
+ __________________________________________________________________________
+ This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
+ __________________________________________________________________________
+
+ begin :
+ email : nguyentd@ornl.gov
+ ***************************************************************************/
+
+#include <iostream>
+#include <cassert>
+#include <math.h>
+
+#include "lal_dipole_long_lj.h"
+
+using namespace std;
+using namespace LAMMPS_AL;
+
+static DipoleLongLJ<PRECISION,ACC_PRECISION> DPLJMF;
+
+// ---------------------------------------------------------------------------
+// Allocate memory on host and device and copy constants to device
+// ---------------------------------------------------------------------------
+int dplj_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 inum,
+ const int nall, const int max_nbors, const int maxspecial,
+ const double cell_size, int &gpu_mode, FILE *screen,
+ double **host_cut_ljsq, const double host_cut_coulsq,
+ double *host_special_coul, const double qqrd2e,
+ const double g_ewald) {
+ DPLJMF.clear();
+ gpu_mode=DPLJMF.device->gpu_mode();
+ double gpu_split=DPLJMF.device->particle_split();
+ int first_gpu=DPLJMF.device->first_device();
+ int last_gpu=DPLJMF.device->last_device();
+ int world_me=DPLJMF.device->world_me();
+ int gpu_rank=DPLJMF.device->gpu_rank();
+ int procs_per_gpu=DPLJMF.device->procs_per_gpu();
+
+ DPLJMF.device->init_message(screen,"lj/cut/dipole/long",first_gpu,last_gpu);
+
+ bool message=false;
+ if (DPLJMF.device->replica_me()==0 && screen)
+ message=true;
+
+ if (message) {
+ fprintf(screen,"Initializing Device and compiling on process 0...");
+ fflush(screen);
+ }
+
+ int init_ok=0;
+ if (world_me==0)
+ init_ok=DPLJMF.init(ntypes, cutsq, 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, g_ewald);
+
+ DPLJMF.device->world_barrier();
+ if (message)
+ fprintf(screen,"Done.\n");
+
+ for (int i=0; i<procs_per_gpu; i++) {
+ if (message) {
+ if (last_gpu-first_gpu==0)
+ fprintf(screen,"Initializing Device %d on core %d...",first_gpu,i);
+ else
+ fprintf(screen,"Initializing Devices %d-%d on core %d...",first_gpu,
+ last_gpu,i);
+ fflush(screen);
+ }
+ if (gpu_rank==i && world_me!=0)
+ init_ok=DPLJMF.init(ntypes, cutsq, 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, g_ewald);
+
+ DPLJMF.device->gpu_barrier();
+ if (message)
+ fprintf(screen,"Done.\n");
+ }
+ if (message)
+ fprintf(screen,"\n");
+
+ if (init_ok==0)
+ DPLJMF.estimate_gpu_overhead();
+ return init_ok;
+}
+
+void dplj_gpu_clear() {
+ DPLJMF.clear();
+}
+
+int** dplj_gpu_compute_n(const int ago, const int inum_full,
+ const int nall, double **host_x, int *host_type,
+ double *sublo, double *subhi, tagint *tag, int **nspecial,
+ tagint **special, const bool eflag, const bool vflag,
+ const bool eatom, const bool vatom, int &host_start,
+ int **ilist, int **jnum, const double cpu_time,
+ bool &success, double *host_q, double **host_mu,
+ double *boxlo, double *prd) {
+ return DPLJMF.compute(ago, inum_full, nall, host_x, host_type, sublo,
+ subhi, tag, nspecial, special, eflag, vflag, eatom,
+ vatom, host_start, ilist, jnum, cpu_time, success,
+ host_q, host_mu, boxlo, prd);
+}
+
+void dplj_gpu_compute(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,
+ double **host_mu, const int nlocal, double *boxlo, double *prd) {
+ DPLJMF.compute(ago,inum_full,nall,host_x,host_type,ilist,numj,firstneigh,eflag,
+ vflag,eatom,vatom,host_start,cpu_time,success,host_q,host_mu,
+ nlocal,boxlo,prd);
+}
+
+double dplj_gpu_bytes() {
+ return DPLJMF.host_memory_usage();
+}
+
+
diff --git a/lib/gpu/lal_lj_expand_coul_long.cpp b/lib/gpu/lal_lj_expand_coul_long.cpp
new file mode 100644
index 0000000000..6a3da9b686
--- /dev/null
+++ b/lib/gpu/lal_lj_expand_coul_long.cpp
@@ -0,0 +1,183 @@
+/***************************************************************************
+ lj_expand_coul_long.cpp
+ --------------------------
+ Trung Nguyen (Northwestern)
+
+ Class for acceleration of the lj/expand/coul/long pair style.
+
+ __________________________________________________________________________
+ This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
+ __________________________________________________________________________
+
+ begin :
+ email : trung.nguyen@northwestern.edu
+ ***************************************************************************/
+
+#if defined(USE_OPENCL)
+#include "lj_expand_coul_long_cl.h"
+#elif defined(USE_CUDART)
+const char *lj_expand_coul_long=0;
+#else
+#include "lj_expand_coul_long_cubin.h"
+#endif
+
+#include "lal_lj_expand_coul_long.h"
+#include <cassert>
+using namespace LAMMPS_AL;
+#define LJExpandCoulLongT LJExpandCoulLong<numtyp, acctyp>
+
+extern Device<PRECISION,ACC_PRECISION> device;
+
+template <class numtyp, class acctyp>
+LJExpandCoulLongT::LJExpandCoulLong() : BaseCharge<numtyp,acctyp>(),
+ _allocated(false) {
+}
+
+template <class numtyp, class acctyp>
+LJExpandCoulLongT::~LJExpandCoulLong() {
+ clear();
+}
+
+template <class numtyp, class acctyp>
+int LJExpandCoulLongT::bytes_per_atom(const int max_nbors) const {
+ return this->bytes_per_atom_atomic(max_nbors);
+}
+
+template <class numtyp, class acctyp>
+int LJExpandCoulLongT::init(const int ntypes,
+ double **host_cutsq, double **host_lj1,
+ double **host_lj2, double **host_lj3,
+ double **host_lj4, double **host_offset, double **host_shift,
+ 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 double g_ewald) {
+ int success;
+ success=this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,gpu_split,
+ _screen,lj_expand_coul_long,"k_lj_expand_coul_long");
+ if (success!=0)
+ return success;
+
+ // If atom type constants fit in shared memory use fast kernel
+ int lj_types=ntypes;
+ shared_types=false;
+ int max_shared_types=this->device->max_shared_types();
+ 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<numtyp> host_write(lj_types*lj_types*32,*(this->ucl_device),
+ UCL_WRITE_ONLY);
+
+ for (int i=0; i<lj_types*lj_types; i++)
+ host_write[i]=0.0;
+
+ lj1.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
+ this->atom->type_pack4(ntypes,lj_types,lj1,host_write,host_lj1,host_lj2,
+ host_cutsq, host_cut_ljsq);
+
+ lj3.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
+ this->atom->type_pack4(ntypes,lj_types,lj3,host_write,host_lj3,host_lj4,
+ host_offset,host_shift);
+
+ 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;
+ _g_ewald=g_ewald;
+
+ _allocated=true;
+ this->_max_bytes=lj1.row_bytes()+lj3.row_bytes()+sp_lj.row_bytes();
+ return 0;
+}
+
+template <class numtyp, class acctyp>
+void LJExpandCoulLongT::reinit(const int ntypes, double **host_cutsq, double **host_lj1,
+ double **host_lj2, double **host_lj3, double **host_lj4,
+ double **host_offset, double **host_shift, double **host_cut_ljsq) {
+ // Allocate a host write buffer for data initialization
+ UCL_H_Vec<numtyp> host_write(_lj_types*_lj_types*32,*(this->ucl_device),
+ UCL_WRITE_ONLY);
+
+ for (int i=0; i<_lj_types*_lj_types; i++)
+ host_write[i]=0.0;
+
+ this->atom->type_pack4(ntypes,_lj_types,lj1,host_write,host_lj1,host_lj2,
+ host_cutsq, host_cut_ljsq);
+ this->atom->type_pack4(ntypes,_lj_types,lj3,host_write,host_lj3,host_lj4,
+ host_offset,host_shift);
+}
+
+template <class numtyp, class acctyp>
+void LJExpandCoulLongT::clear() {
+ if (!_allocated)
+ return;
+ _allocated=false;
+
+ lj1.clear();
+ lj3.clear();
+ sp_lj.clear();
+ this->clear_atomic();
+}
+
+template <class numtyp, class acctyp>
+double LJExpandCoulLongT::host_memory_usage() const {
+ return this->host_memory_usage_atomic()+sizeof(LJExpandCoulLong<numtyp,acctyp>);
+}
+
+// ---------------------------------------------------------------------------
+// Calculate energies, forces, and torques
+// ---------------------------------------------------------------------------
+template <class numtyp, class acctyp>
+void LJExpandCoulLongT::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<int>(ceil(static_cast<double>(this->ans->inum())/
+ (BX/this->_threads_per_atom)));
+
+ int ainum=this->ans->inum();
+ 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->x, &lj1, &lj3, &sp_lj,
+ &this->nbor->dev_nbor, &this->_nbor_data->begin(),
+ &this->ans->force, &this->ans->engv, &eflag,
+ &vflag, &ainum, &nbor_pitch, &this->atom->q,
+ &_cut_coulsq, &_qqrd2e, &_g_ewald,
+ &this->_threads_per_atom);
+ } else {
+ this->k_pair.set_size(GX,BX);
+ this->k_pair.run(&this->atom->x, &lj1, &lj3,
+ &_lj_types, &sp_lj, &this->nbor->dev_nbor,
+ &this->_nbor_data->begin(), &this->ans->force,
+ &this->ans->engv, &eflag, &vflag, &ainum,
+ &nbor_pitch, &this->atom->q, &_cut_coulsq,
+ &_qqrd2e, &_g_ewald, &this->_threads_per_atom);
+ }
+ this->time_pair.stop();
+}
+
+template class LJExpandCoulLong<PRECISION,ACC_PRECISION>;
diff --git a/lib/gpu/lal_lj_expand_coul_long.cu b/lib/gpu/lal_lj_expand_coul_long.cu
new file mode 100644
index 0000000000..eaea5b0e12
--- /dev/null
+++ b/lib/gpu/lal_lj_expand_coul_long.cu
@@ -0,0 +1,265 @@
+// **************************************************************************
+// lj_coul_long.cu
+// -------------------
+// Trung Nguyen (Northwestern)
+//
+// Device code for acceleration of the lj/expand/coul/long pair style
+//
+// __________________________________________________________________________
+// This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
+// __________________________________________________________________________
+//
+// begin :
+// email : trung.nguyen@northwestern.edu
+// ***************************************************************************/
+
+#ifdef NV_KERNEL
+
+#include "lal_aux_fun1.h"
+#ifndef _DOUBLE_DOUBLE
+texture<float4> pos_tex;
+texture<float> q_tex;
+#else
+texture<int4,1> pos_tex;
+texture<int2> q_tex;
+#endif
+
+#else
+#define pos_tex x_
+#define q_tex q_
+#endif
+
+__kernel void k_lj_expand_coul_long(const __global numtyp4 *restrict x_,
+ const __global numtyp4 *restrict lj1,
+ const __global numtyp4 *restrict lj3,
+ const int lj_types,
+ const __global numtyp *restrict sp_lj_in,
+ const __global int *dev_nbor,
+ const __global int *dev_packed,
+ __global acctyp4 *restrict ans,
+ __global acctyp *restrict engv,
+ const int eflag, const int vflag, const int inum,
+ const int nbor_pitch,
+ const __global numtyp *restrict q_,
+ const numtyp cut_coulsq, const numtyp qqrd2e,
+ const numtyp g_ewald, const int t_per_atom) {
+ int tid, ii, offset;
+ atom_info(t_per_atom,ii,tid,offset);
+
+ __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];
+
+ acctyp energy=(acctyp)0;
+ acctyp e_coul=(acctyp)0;
+ acctyp4 f;
+ f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
+ acctyp virial[6];
+ for (int i=0; i<6; i++)
+ virial[i]=(acctyp)0;
+
+ if (ii<inum) {
+ int nbor, nbor_end;
+ int i, numj;
+ __local int n_stride;
+ nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
+ n_stride,nbor_end,nbor);
+
+ numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
+ numtyp qtmp; fetch(qtmp,i,q_tex);
+ int itype=ix.w;
+
+ for ( ; nbor<nbor_end; nbor+=n_stride) {
+ int j=dev_packed[nbor];
+
+ numtyp factor_lj, factor_coul;
+ factor_lj = sp_lj[sbmask(j)];
+ factor_coul = (numtyp)1.0-sp_lj[sbmask(j)+4];
+ j &= NEIGHMASK;
+
+ numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
+ int jtype=jx.w;
+
+ // Compute r12
+ numtyp delx = ix.x-jx.x;
+ numtyp dely = ix.y-jx.y;
+ numtyp delz = ix.z-jx.z;
+ numtyp rsq = delx*delx+dely*dely+delz*delz;
+
+ int mtype=itype*lj_types+jtype;
+ if (rsq<lj1[mtype].z) {
+ numtyp r2inv=ucl_recip(rsq);
+ numtyp forcecoul, force_lj, force, r6inv, prefactor, _erfc;
+
+ if (rsq < lj1[mtype].w) {
+ r6inv = r2inv*r2inv*r2inv;
+ force_lj = factor_lj*r6inv*(lj1[mtype].x*r6inv-lj1[mtype].y);
+ } else
+ force_lj = (numtyp)0.0;
+
+ if (rsq < cut_coulsq) {
+ numtyp r = ucl_rsqrt(r2inv);
+ numtyp grij = g_ewald * r;
+ numtyp expm2 = ucl_exp(-grij*grij);
+ numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij);
+ _erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+ fetch(prefactor,j,q_tex);
+ prefactor *= qqrd2e * qtmp/r;
+ forcecoul = prefactor * (_erfc + EWALD_F*grij*expm2-factor_coul);
+ } else
+ forcecoul = (numtyp)0.0;
+
+ force = (force_lj + forcecoul) * r2inv;
+
+ f.x+=delx*force;
+ f.y+=dely*force;
+ f.z+=delz*force;
+
+ if (eflag>0) {
+ if (rsq < cut_coulsq)
+ e_coul += prefactor*(_erfc-factor_coul);
+ if (rsq < lj1[mtype].w) {
+ numtyp e=r6inv*(lj3[mtype].x*r6inv-lj3[mtype].y);
+ energy+=factor_lj*(e-lj3[mtype].z);
+ }
+ }
+ 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_q(f,energy,e_coul,virial,ii,inum,tid,t_per_atom,offset,eflag,
+ vflag,ans,engv);
+ } // if ii
+}
+
+__kernel void k_lj_expand_coul_long_fast(const __global numtyp4 *restrict x_,
+ const __global numtyp4 *restrict lj1_in,
+ const __global numtyp4 *restrict lj3_in,
+ const __global numtyp *restrict sp_lj_in,
+ const __global int *dev_nbor,
+ const __global int *dev_packed,
+ __global acctyp4 *restrict ans,
+ __global acctyp *restrict engv,
+ const int eflag, const int vflag,
+ const int inum, const int nbor_pitch,
+ const __global numtyp *restrict q_,
+ const numtyp cut_coulsq, const numtyp qqrd2e,
+ const numtyp g_ewald, const int t_per_atom) {
+ int tid, ii, offset;
+ atom_info(t_per_atom,ii,tid,offset);
+
+ __local numtyp4 lj1[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
+ __local numtyp4 lj3[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
+ __local numtyp sp_lj[8];
+ if (tid<8)
+ sp_lj[tid]=sp_lj_in[tid];
+ if (tid<MAX_SHARED_TYPES*MAX_SHARED_TYPES) {
+ lj1[tid]=lj1_in[tid];
+ if (eflag>0)
+ lj3[tid]=lj3_in[tid];
+ }
+
+ acctyp energy=(acctyp)0;
+ acctyp e_coul=(acctyp)0;
+ acctyp4 f;
+ f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
+ acctyp virial[6];
+ for (int i=0; i<6; i++)
+ virial[i]=(acctyp)0;
+
+ __syncthreads();
+
+ if (ii<inum) {
+ int nbor, nbor_end;
+ int i, numj;
+ __local int n_stride;
+ nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
+ n_stride,nbor_end,nbor);
+
+ numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
+ numtyp qtmp; fetch(qtmp,i,q_tex);
+ int iw=ix.w;
+ int itype=fast_mul((int)MAX_SHARED_TYPES,iw);
+
+ for ( ; nbor<nbor_end; nbor+=n_stride) {
+ int j=dev_packed[nbor];
+
+ numtyp factor_lj, factor_coul;
+ factor_lj = sp_lj[sbmask(j)];
+ factor_coul = (numtyp)1.0-sp_lj[sbmask(j)+4];
+ j &= NEIGHMASK;
+
+ numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
+ int mtype=itype+jx.w;
+
+ // Compute r12
+ numtyp delx = ix.x-jx.x;
+ numtyp dely = ix.y-jx.y;
+ numtyp delz = ix.z-jx.z;
+ numtyp rsq = delx*delx+dely*dely+delz*delz;
+
+ if (rsq<lj1[mtype].z) {
+ numtyp r2inv=ucl_recip(rsq);
+ numtyp forcecoul, force_lj, force, r6inv, prefactor, _erfc;
+
+ if (rsq < lj1[mtype].w) {
+ r6inv = r2inv*r2inv*r2inv;
+ force_lj = factor_lj*r6inv*(lj1[mtype].x*r6inv-lj1[mtype].y);
+ } else
+ force_lj = (numtyp)0.0;
+
+ if (rsq < cut_coulsq) {
+ numtyp r = ucl_rsqrt(r2inv);
+ numtyp grij = g_ewald * r;
+ numtyp expm2 = ucl_exp(-grij*grij);
+ numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij);
+ _erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+ fetch(prefactor,j,q_tex);
+ prefactor *= qqrd2e * qtmp/r;
+ forcecoul = prefactor * (_erfc + EWALD_F*grij*expm2-factor_coul);
+ } else
+ forcecoul = (numtyp)0.0;
+
+ force = (force_lj + forcecoul) * r2inv;
+
+ f.x+=delx*force;
+ f.y+=dely*force;
+ f.z+=delz*force;
+
+ if (eflag>0) {
+ if (rsq < cut_coulsq)
+ e_coul += prefactor*(_erfc-factor_coul);
+ if (rsq < lj1[mtype].w) {
+ numtyp e=r6inv*(lj3[mtype].x*r6inv-lj3[mtype].y);
+ energy+=factor_lj*(e-lj3[mtype].z);
+ }
+ }
+ 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_q(f,energy,e_coul,virial,ii,inum,tid,t_per_atom,offset,eflag,
+ vflag,ans,engv);
+ } // if ii
+}
+
diff --git a/lib/gpu/lal_lj_expand_coul_long.h b/lib/gpu/lal_lj_expand_coul_long.h
new file mode 100644
index 0000000000..c5c7a468a9
--- /dev/null
+++ b/lib/gpu/lal_lj_expand_coul_long.h
@@ -0,0 +1,88 @@
+/***************************************************************************
+ lj_expand_coul_long.h
+ -------------------
+ Trung Nguyen (Northwestern)
+
+ Class for acceleration of the lj/expand/coul/long pair style.
+
+ __________________________________________________________________________
+ This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
+ __________________________________________________________________________
+
+ begin :
+ email : trung.nguyen@northwestern.edu
+ ***************************************************************************/
+
+#ifndef LAL_LJ_EXPAND_COUL_LONG_H
+#define LAL_LJ_EXPAND_COUL_LONG_H
+
+#include "lal_base_charge.h"
+
+namespace LAMMPS_AL {
+
+template <class numtyp, class acctyp>
+class LJExpandCoulLong : public BaseCharge<numtyp, acctyp> {
+ public:
+ LJExpandCoulLong();
+ ~LJExpandCoulLong();
+
+ /// 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
+ *
+ * Returns:
+ * - 0 if successfull
+ * - -1 if fix gpu not found
+ * - -3 if there is an out of memory error
+ * - -4 if the GPU library was not compiled for GPU
+ * - -5 Double precision is not supported on card **/
+ int init(const int ntypes, double **host_cutsq,
+ double **host_lj1, double **host_lj2, double **host_lj3,
+ double **host_lj4, double **host_offset, double **host_shift, 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 double g_ewald);
+
+ /// Send updated coeffs from host to device (to be compatible with fix adapt)
+ void reinit(const int ntypes, double **host_cutsq,
+ double **host_lj1, double **host_lj2, double **host_lj3,
+ double **host_lj4, double **host_offset, double **host_shift, double **host_cut_ljsq);
+
+ /// 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 = lj1, lj1.y = lj2, lj1.z = cutsq, lj1.w = cutsq_vdw
+ UCL_D_Vec<numtyp4> lj1;
+ /// lj3.x = lj3, lj3.y = lj4, lj3.z = offset, lj3.w = shift
+ UCL_D_Vec<numtyp4> lj3;
+ /// Special LJ values [0-3] and Special Coul values [4-7]
+ UCL_D_Vec<numtyp> 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, _g_ewald;
+
+ private:
+ bool _allocated;
+ void loop(const bool _eflag, const bool _vflag);
+};
+
+}
+
+#endif
diff --git a/lib/gpu/lal_lj_expand_coul_long_ext.cpp b/lib/gpu/lal_lj_expand_coul_long_ext.cpp
new file mode 100644
index 0000000000..e669316e4e
--- /dev/null
+++ b/lib/gpu/lal_lj_expand_coul_long_ext.cpp
@@ -0,0 +1,152 @@
+/***************************************************************************
+ lj_coul_long_ext.cpp
+ ------------------------
+ Trung Nguyen (Northwestern)
+
+ Functions for LAMMPS access to lj/expand/coul/long acceleration routines.
+
+ __________________________________________________________________________
+ This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
+ __________________________________________________________________________
+
+ begin :
+ email : trung.nguyen@northwestern.edu
+ ***************************************************************************/
+
+#include <iostream>
+#include <cassert>
+#include <math.h>
+
+#include "lal_lj_expand_coul_long.h"
+
+using namespace std;
+using namespace LAMMPS_AL;
+
+static LJExpandCoulLong<PRECISION,ACC_PRECISION> LJECLMF;
+
+// ---------------------------------------------------------------------------
+// Allocate memory on host and device and copy constants to device
+// ---------------------------------------------------------------------------
+int ljecl_gpu_init(const int ntypes, double **cutsq, double **host_lj1,
+ double **host_lj2, double **host_lj3, double **host_lj4,
+ double **offset, double **shift, 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 double g_ewald) {
+ LJECLMF.clear();
+ gpu_mode=LJECLMF.device->gpu_mode();
+ double gpu_split=LJECLMF.device->particle_split();
+ int first_gpu=LJECLMF.device->first_device();
+ int last_gpu=LJECLMF.device->last_device();
+ int world_me=LJECLMF.device->world_me();
+ int gpu_rank=LJECLMF.device->gpu_rank();
+ int procs_per_gpu=LJECLMF.device->procs_per_gpu();
+
+ LJECLMF.device->init_message(screen,"lj/expand/coul/long",first_gpu,last_gpu);
+
+ bool message=false;
+ if (LJECLMF.device->replica_me()==0 && screen)
+ message=true;
+
+ if (message) {
+ fprintf(screen,"Initializing Device and compiling on process 0...");
+ fflush(screen);
+ }
+
+ int init_ok=0;
+ if (world_me==0)
+ init_ok=LJECLMF.init(ntypes, cutsq, host_lj1, host_lj2, host_lj3, host_lj4,
+ offset, shift, special_lj, inum, nall, 300, maxspecial,
+ cell_size, gpu_split, screen, host_cut_ljsq,
+ host_cut_coulsq, host_special_coul, qqrd2e, g_ewald);
+
+ LJECLMF.device->world_barrier();
+ if (message)
+ fprintf(screen,"Done.\n");
+
+ for (int i=0; i<procs_per_gpu; i++) {
+ if (message) {
+ if (last_gpu-first_gpu==0)
+ fprintf(screen,"Initializing Device %d on core %d...",first_gpu,i);
+ else
+ fprintf(screen,"Initializing Devices %d-%d on core %d...",first_gpu,
+ last_gpu,i);
+ fflush(screen);
+ }
+ if (gpu_rank==i && world_me!=0)
+ init_ok=LJECLMF.init(ntypes, cutsq, host_lj1, host_lj2, host_lj3, host_lj4,
+ offset, shift, special_lj, inum, nall, 300, maxspecial,
+ cell_size, gpu_split, screen, host_cut_ljsq,
+ host_cut_coulsq, host_special_coul, qqrd2e, g_ewald);
+
+ LJECLMF.device->gpu_barrier();
+ if (message)
+ fprintf(screen,"Done.\n");
+ }
+ if (message)
+ fprintf(screen,"\n");
+
+ if (init_ok==0)
+ LJECLMF.estimate_gpu_overhead();
+ return init_ok;
+}
+
+// ---------------------------------------------------------------------------
+// Copy updated coeffs from host to device
+// ---------------------------------------------------------------------------
+void ljecl_gpu_reinit(const int ntypes, double **cutsq, double **host_lj1,
+ double **host_lj2, double **host_lj3, double **host_lj4,
+ double **offset, double **shift, double **host_cut_ljsq) {
+ int world_me=LJECLMF.device->world_me();
+ int gpu_rank=LJECLMF.device->gpu_rank();
+ int procs_per_gpu=LJECLMF.device->procs_per_gpu();
+
+ if (world_me==0)
+ LJECLMF.reinit(ntypes, cutsq, host_lj1, host_lj2, host_lj3, host_lj4,
+ offset, shift, host_cut_ljsq);
+ LJECLMF.device->world_barrier();
+
+ for (int i=0; i<procs_per_gpu; i++) {
+ if (gpu_rank==i && world_me!=0)
+ LJECLMF.reinit(ntypes, cutsq, host_lj1, host_lj2, host_lj3, host_lj4,
+ offset, shift, host_cut_ljsq);
+ LJECLMF.device->gpu_barrier();
+ }
+}
+
+void ljecl_gpu_clear() {
+ LJECLMF.clear();
+}
+
+int** ljecl_gpu_compute_n(const int ago, const int inum_full,
+ const int nall, double **host_x, int *host_type,
+ double *sublo, double *subhi, tagint *tag, int **nspecial,
+ tagint **special, const bool eflag, const bool vflag,
+ const bool eatom, const bool vatom, int &host_start,
+ int **ilist, int **jnum, const double cpu_time,
+ bool &success, double *host_q, double *boxlo,
+ double *prd) {
+ return LJECLMF.compute(ago, inum_full, nall, host_x, host_type, sublo,
+ subhi, tag, nspecial, special, eflag, vflag, eatom,
+ vatom, host_start, ilist, jnum, cpu_time, success,
+ host_q, boxlo, prd);
+}
+
+void ljecl_gpu_compute(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,
+ const int nlocal, double *boxlo, double *prd) {
+ LJECLMF.compute(ago,inum_full,nall,host_x,host_type,ilist,numj,
+ firstneigh,eflag,vflag,eatom,vatom,host_start,cpu_time,success,
+ host_q,nlocal,boxlo,prd);
+}
+
+double ljecl_gpu_bytes() {
+ return LJECLMF.host_memory_usage();
+}
+
+
diff --git a/src/GPU/Install.sh b/src/GPU/Install.sh
index 88f47a3dc4..a850a9b2b3 100755
--- a/src/GPU/Install.sh
+++ b/src/GPU/Install.sh
@@ -133,6 +133,8 @@ action pppm_gpu.cpp pppm.cpp
action pppm_gpu.h pppm.cpp
action pair_ufm_gpu.cpp
action pair_ufm_gpu.h
+action pair_lj_cut_dipole_long_gpu.cpp pair_lj_cut_dipole_long.cpp
+action pair_lj_cut_dipole_long_gpu.h pair_lj_cut_dipole_long.cpp
# edit 2 Makefile.package files to include/exclude package info
diff --git a/src/GPU/pair_lj_cut_dipole_long_gpu.cpp b/src/GPU/pair_lj_cut_dipole_long_gpu.cpp
new file mode 100644
index 0000000000..045b26b7ce
--- /dev/null
+++ b/src/GPU/pair_lj_cut_dipole_long_gpu.cpp
@@ -0,0 +1,434 @@
+/* ----------------------------------------------------------------------
+ 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: Trung Dac Nguyen (Northwestern)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "pair_lj_cut_dipole_long_gpu.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "comm.h"
+#include "force.h"
+#include "kspace.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "integrate.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+#include "neigh_request.h"
+#include "universe.h"
+#include "update.h"
+#include "domain.h"
+#include "gpu_extra.h"
+
+#define EWALD_F 1.12837917
+#define EWALD_P 0.3275911
+#define A1 0.254829592
+#define A2 -0.284496736
+#define A3 1.421413741
+#define A4 -1.453152027
+#define A5 1.061405429
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+// External functions from cuda library for atom decomposition
+
+int dplj_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,
+ double **host_cut_ljsq, const double host_cut_coulsq,
+ double *host_special_coul, const double qqrd2e, const double g_ewald);
+void dplj_gpu_clear();
+int ** dplj_gpu_compute_n(const int ago, const int inum,
+ const int nall, double **host_x, int *host_type,
+ double *sublo, double *subhi, tagint *tag,
+ int **nspecial, tagint **special, const bool eflag,
+ const bool vflag, const bool eatom, const bool vatom,
+ int &host_start, int **ilist, int **jnum,
+ const double cpu_time, bool &success,
+ double *host_q, double **host_mu,
+ double *boxlo, double *prd);
+void dplj_gpu_compute(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 **host_mu,
+ const int nlocal, double *boxlo, double *prd);
+double dplj_gpu_bytes();
+
+/* ---------------------------------------------------------------------- */
+
+PairLJCutDipoleLongGPU::PairLJCutDipoleLongGPU(LAMMPS *lmp) : PairLJCutDipoleLong(lmp),
+ gpu_mode(GPU_FORCE)
+{
+ respa_enable = 0;
+ reinitflag = 0;
+ cpu_time = 0.0;
+ GPU_EXTRA::gpu_ready(lmp->modify, lmp->error);
+}
+
+/* ----------------------------------------------------------------------
+ free all arrays
+------------------------------------------------------------------------- */
+
+PairLJCutDipoleLongGPU::~PairLJCutDipoleLongGPU()
+{
+ dplj_gpu_clear();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJCutDipoleLongGPU::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;
+ int *ilist, *numneigh, **firstneigh;
+ if (gpu_mode != GPU_FORCE) {
+ inum = atom->nlocal;
+ firstneigh = dplj_gpu_compute_n(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, &ilist, &numneigh, cpu_time,
+ success, atom->q, atom->mu, domain->boxlo,
+ domain->prd);
+ } else {
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+ dplj_gpu_compute(neighbor->ago, inum, nall, atom->x, atom->type,
+ ilist, numneigh, firstneigh, eflag, vflag, eflag_atom,
+ vflag_atom, host_start, cpu_time, success, atom->q,
+ atom->mu, atom->nlocal, domain->boxlo, domain->prd);
+ }
+ if (!success)
+ error->one(FLERR,"Insufficient memory on accelerator");
+
+ if (host_start<inum) {
+ cpu_time = MPI_Wtime();
+ cpu_compute(host_start, inum, eflag, vflag, ilist, numneigh, firstneigh);
+ cpu_time = MPI_Wtime() - cpu_time;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairLJCutDipoleLongGPU::init_style()
+{
+ if (!atom->q_flag || !atom->mu_flag || !atom->torque_flag)
+ error->all(FLERR,"Pair dipole/cut/gpu requires atom attributes q, mu, torque");
+
+ if (force->newton_pair)
+ error->all(FLERR,"Cannot use newton pair with dipole/cut/gpu pair style");
+
+ if (strcmp(update->unit_style,"electron") == 0)
+ error->all(FLERR,"Cannot (yet) use 'electron' units with dipoles");
+
+ // 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;
+
+ cut_coulsq = cut_coul * cut_coul;
+
+ // insure use of KSpace long-range solver, set g_ewald
+
+ if (force->kspace == NULL)
+ error->all(FLERR,"Pair style requires a KSpace style");
+ g_ewald = force->kspace->g_ewald;
+
+ // setup force tables
+
+ if (ncoultablebits) init_tables(cut_coul,NULL);
+
+ int maxspecial=0;
+ if (atom->molecular)
+ maxspecial=atom->maxspecial;
+ int success = dplj_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, cut_ljsq, cut_coulsq,
+ force->special_coul, force->qqrd2e, g_ewald);
+ GPU_EXTRA::check_flag(success,error,world);
+
+ if (gpu_mode == GPU_FORCE) {
+ int irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->full = 1;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLJCutDipoleLongGPU::memory_usage()
+{
+ double bytes = Pair::memory_usage();
+ return bytes + dplj_gpu_bytes();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJCutDipoleLongGPU::cpu_compute(int start, int inum, int eflag, int vflag,
+ int *ilist, int *numneigh,
+ int **firstneigh)
+{
+ int i,j,ii,jj,jnum,itype,jtype;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz;
+ double rsq,r,rinv,r2inv,r6inv;
+ double forcecoulx,forcecouly,forcecoulz,fforce;
+ double tixcoul,tiycoul,tizcoul,tjxcoul,tjycoul,tjzcoul;
+ double fx,fy,fz,fdx,fdy,fdz,fax,fay,faz;
+ double pdotp,pidotr,pjdotr,pre1,pre2,pre3;
+ double grij,expm2,t,erfc;
+ double g0,g1,g2,b0,b1,b2,b3,d0,d1,d2,d3;
+ double zdix,zdiy,zdiz,zdjx,zdjy,zdjz,zaix,zaiy,zaiz,zajx,zajy,zajz;
+ double g0b1_g1b2_g2b3,g0d1_g1d2_g2d3;
+ double forcelj,factor_coul,factor_lj,facm1;
+ double evdwl,ecoul;
+ int *jlist;
+
+ evdwl = ecoul = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ double **mu = atom->mu;
+ double **torque = atom->torque;
+ int *type = atom->type;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ double qqrd2e = force->qqrd2e;
+
+ pre1 = 2.0 * g_ewald / MY_PIS;
+ pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS;
+ pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS;
+
+ // 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];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ 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;
+ rinv = sqrt(r2inv);
+
+ if (rsq < cut_coulsq) {
+ r = sqrt(rsq);
+ grij = g_ewald * r;
+ expm2 = exp(-grij*grij);
+ t = 1.0 / (1.0 + EWALD_P*grij);
+ erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+
+ pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2];
+ pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
+ pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
+
+ g0 = qtmp*q[j];
+ g1 = qtmp*pjdotr - q[j]*pidotr + pdotp;
+ g2 = -pidotr*pjdotr;
+
+ if (factor_coul > 0.0) {
+ b0 = erfc * rinv;
+ b1 = (b0 + pre1*expm2) * r2inv;
+ b2 = (3.0*b1 + pre2*expm2) * r2inv;
+ b3 = (5.0*b2 + pre3*expm2) * r2inv;
+
+ g0b1_g1b2_g2b3 = g0*b1 + g1*b2 + g2*b3;
+ fdx = delx * g0b1_g1b2_g2b3 -
+ b1 * (qtmp*mu[j][0] - q[j]*mu[i][0]) +
+ b2 * (pjdotr*mu[i][0] + pidotr*mu[j][0]);
+ fdy = dely * g0b1_g1b2_g2b3 -
+ b1 * (qtmp*mu[j][1] - q[j]*mu[i][1]) +
+ b2 * (pjdotr*mu[i][1] + pidotr*mu[j][1]);
+ fdz = delz * g0b1_g1b2_g2b3 -
+ b1 * (qtmp*mu[j][2] - q[j]*mu[i][2]) +
+ b2 * (pjdotr*mu[i][2] + pidotr*mu[j][2]);
+
+ zdix = delx * (q[j]*b1 + b2*pjdotr) - b1*mu[j][0];
+ zdiy = dely * (q[j]*b1 + b2*pjdotr) - b1*mu[j][1];
+ zdiz = delz * (q[j]*b1 + b2*pjdotr) - b1*mu[j][2];
+ zdjx = delx * (-qtmp*b1 + b2*pidotr) - b1*mu[i][0];
+ zdjy = dely * (-qtmp*b1 + b2*pidotr) - b1*mu[i][1];
+ zdjz = delz * (-qtmp*b1 + b2*pidotr) - b1*mu[i][2];
+
+ if (factor_coul < 1.0) {
+ fdx *= factor_coul;
+ fdy *= factor_coul;
+ fdz *= factor_coul;
+ zdix *= factor_coul;
+ zdiy *= factor_coul;
+ zdiz *= factor_coul;
+ zdjx *= factor_coul;
+ zdjy *= factor_coul;
+ zdjz *= factor_coul;
+ }
+ } else {
+ fdx = fdy = fdz = 0.0;
+ zdix = zdiy = zdiz = 0.0;
+ zdjx = zdjy = zdjz = 0.0;
+ }
+
+ if (factor_coul < 1.0) {
+ d0 = (erfc - 1.0) * rinv;
+ d1 = (d0 + pre1*expm2) * r2inv;
+ d2 = (3.0*d1 + pre2*expm2) * r2inv;
+ d3 = (5.0*d2 + pre3*expm2) * r2inv;
+
+ g0d1_g1d2_g2d3 = g0*d1 + g1*d2 + g2*d3;
+ fax = delx * g0d1_g1d2_g2d3 -
+ d1 * (qtmp*mu[j][0] - q[j]*mu[i][0]) +
+ d2 * (pjdotr*mu[i][0] + pidotr*mu[j][0]);
+ fay = dely * g0d1_g1d2_g2d3 -
+ d1 * (qtmp*mu[j][1] - q[j]*mu[i][1]) +
+ d2 * (pjdotr*mu[i][1] + pidotr*mu[j][1]);
+ faz = delz * g0d1_g1d2_g2d3 -
+ d1 * (qtmp*mu[j][2] - q[j]*mu[i][2]) +
+ d2 * (pjdotr*mu[i][2] + pidotr*mu[j][2]);
+
+ zaix = delx * (q[j]*d1 + d2*pjdotr) - d1*mu[j][0];
+ zaiy = dely * (q[j]*d1 + d2*pjdotr) - d1*mu[j][1];
+ zaiz = delz * (q[j]*d1 + d2*pjdotr) - d1*mu[j][2];
+ zajx = delx * (-qtmp*d1 + d2*pidotr) - d1*mu[i][0];
+ zajy = dely * (-qtmp*d1 + d2*pidotr) - d1*mu[i][1];
+ zajz = delz * (-qtmp*d1 + d2*pidotr) - d1*mu[i][2];
+
+ if (factor_coul > 0.0) {
+ facm1 = 1.0 - factor_coul;
+ fax *= facm1;
+ fay *= facm1;
+ faz *= facm1;
+ zaix *= facm1;
+ zaiy *= facm1;
+ zaiz *= facm1;
+ zajx *= facm1;
+ zajy *= facm1;
+ zajz *= facm1;
+ }
+ } else {
+ fax = fay = faz = 0.0;
+ zaix = zaiy = zaiz = 0.0;
+ zajx = zajy = zajz = 0.0;
+ }
+
+ forcecoulx = fdx + fax;
+ forcecouly = fdy + fay;
+ forcecoulz = fdz + faz;
+
+ tixcoul = mu[i][1]*(zdiz + zaiz) - mu[i][2]*(zdiy + zaiy);
+ tiycoul = mu[i][2]*(zdix + zaix) - mu[i][0]*(zdiz + zaiz);
+ tizcoul = mu[i][0]*(zdiy + zaiy) - mu[i][1]*(zdix + zaix);
+ tjxcoul = mu[j][1]*(zdjz + zajz) - mu[j][2]*(zdjy + zajy);
+ tjycoul = mu[j][2]*(zdjx + zajx) - mu[j][0]*(zdjz + zajz);
+ tjzcoul = mu[j][0]*(zdjy + zajy) - mu[j][1]*(zdjx + zajx);
+
+ } else {
+ forcecoulx = forcecouly = forcecoulz = 0.0;
+ tixcoul = tiycoul = tizcoul = 0.0;
+ tjxcoul = tjycoul = tjzcoul = 0.0;
+ }
+
+ // LJ interaction
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ fforce = factor_lj * forcelj*r2inv;
+ } else fforce = 0.0;
+
+ // total force
+
+ fx = qqrd2e*forcecoulx + delx*fforce;
+ fy = qqrd2e*forcecouly + dely*fforce;
+ fz = qqrd2e*forcecoulz + delz*fforce;
+
+ // force & torque accumulation
+
+ f[i][0] += fx;
+ f[i][1] += fy;
+ f[i][2] += fz;
+ torque[i][0] += qqrd2e*tixcoul;
+ torque[i][1] += qqrd2e*tiycoul;
+ torque[i][2] += qqrd2e*tizcoul;
+
+ if (eflag) {
+ if (rsq < cut_coulsq && factor_coul > 0.0) {
+ ecoul = qqrd2e*(b0*g0 + b1*g1 + b2*g2);
+ if (factor_coul < 1.0) {
+ ecoul *= factor_coul;
+ ecoul += (1-factor_coul) * qqrd2e * (d0*g0 + d1*g1 + d2*g2);
+ }
+ } else ecoul = 0.0;
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
+ offset[itype][jtype];
+ evdwl *= factor_lj;
+ } else evdwl = 0.0;
+ }
+
+ if (evflag) ev_tally_xyz_full(i,evdwl,ecoul,fx,fy,fz,delx,dely,delz);
+ }
+ }
+ }
+}
diff --git a/src/GPU/pair_lj_cut_dipole_long_gpu.h b/src/GPU/pair_lj_cut_dipole_long_gpu.h
new file mode 100644
index 0000000000..dc72921583
--- /dev/null
+++ b/src/GPU/pair_lj_cut_dipole_long_gpu.h
@@ -0,0 +1,66 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(lj/cut/dipole/long/gpu,PairLJCutDipoleLongGPU)
+
+#else
+
+#ifndef LMP_PAIR_LJ_CUT_DIPOLE_LONG_GPU_H
+#define LMP_PAIR_LJ_CUT_DIPOLE_LONG_GPU_H
+
+#include "pair_lj_cut_dipole_long.h"
+
+namespace LAMMPS_NS {
+
+class PairLJCutDipoleLongGPU : public PairLJCutDipoleLong {
+ public:
+ PairLJCutDipoleLongGPU(LAMMPS *lmp);
+ ~PairLJCutDipoleLongGPU();
+ void cpu_compute(int, int, int, int, int *, int *, int **);
+ void compute(int, int);
+ void init_style();
+ double memory_usage();
+
+ enum { GPU_FORCE, GPU_NEIGH, GPU_HYB_NEIGH };
+
+ private:
+ int gpu_mode;
+ double cpu_time;
+};
+
+}
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Insufficient memory on accelerator
+
+There is insufficient memory on one of the devices specified for the gpu
+package
+
+E: Pair lj/cutdipole/long/gpu requires atom attributes q, mu, torque
+
+The atom style defined does not have this attribute.
+
+E: Cannot use newton pair with lj/cutdipole/long/gpu pair style
+
+Self-explanatory.
+
+E: Cannot (yet) use 'electron' units with dipoles
+
+This feature is not yet supported.
+
+*/
diff --git a/src/GPU/pair_lj_expand_coul_long_gpu.cpp b/src/GPU/pair_lj_expand_coul_long_gpu.cpp
new file mode 100644
index 0000000000..c60c0f847b
--- /dev/null
+++ b/src/GPU/pair_lj_expand_coul_long_gpu.cpp
@@ -0,0 +1,325 @@
+/* ----------------------------------------------------------------------
+ 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: Trung Nguyen (Northwestern)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "pair_lj_expand_coul_long_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"
+#include "gpu_extra.h"
+
+#define EWALD_F 1.12837917
+#define EWALD_P 0.3275911
+#define A1 0.254829592
+#define A2 -0.284496736
+#define A3 1.421413741
+#define A4 -1.453152027
+#define A5 1.061405429
+
+using namespace LAMMPS_NS;
+
+// External functions from cuda library for atom decomposition
+
+int ljecl_gpu_init(const int ntypes, double **cutsq, double **host_lj1,
+ double **host_lj2, double **host_lj3, double **host_lj4,
+ double **offset, double **shift, 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 double g_ewald);
+int ljecl_gpu_reinit(const int ntypes, double **cutsq, double **host_lj1,
+ double **host_lj2, double **host_lj3, double **host_lj4,
+ double **offset, double **shift, double **host_lj_cutsq);
+void ljecl_gpu_clear();
+int ** ljecl_gpu_compute_n(const int ago, const int inum,
+ const int nall, double **host_x, int *host_type,
+ double *sublo, double *subhi, tagint *tag,
+ int **nspecial, tagint **special, const bool eflag,
+ const bool vflag, const bool eatom, const bool vatom,
+ int &host_start, int **ilist, int **jnum,
+ const double cpu_time, bool &success, double *host_q,
+ double *boxlo, double *prd);
+void ljecl_gpu_compute(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,
+ const int nlocal, double *boxlo, double *prd);
+double ljecl_gpu_bytes();
+
+/* ---------------------------------------------------------------------- */
+
+PairLJExpandCoulLongGPU::PairLJExpandCoulLongGPU(LAMMPS *lmp) :
+ PairLJExpandCoulLong(lmp), gpu_mode(GPU_FORCE)
+{
+ respa_enable = 0;
+ cpu_time = 0.0;
+ GPU_EXTRA::gpu_ready(lmp->modify, lmp->error);
+}
+
+/* ----------------------------------------------------------------------
+ free all arrays
+------------------------------------------------------------------------- */
+
+PairLJExpandCoulLongGPU::~PairLJExpandCoulLongGPU()
+{
+ ljecl_gpu_clear();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJExpandCoulLongGPU::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;
+ int *ilist, *numneigh, **firstneigh;
+ if (gpu_mode != GPU_FORCE) {
+ inum = atom->nlocal;
+ firstneigh = ljecl_gpu_compute_n(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, &ilist, &numneigh, cpu_time,
+ success, atom->q, domain->boxlo,
+ domain->prd);
+ } else {
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+ ljecl_gpu_compute(neighbor->ago, inum, nall, atom->x, atom->type,
+ ilist, numneigh, firstneigh, eflag, vflag, eflag_atom,
+ vflag_atom, host_start, cpu_time, success, atom->q,
+ atom->nlocal, domain->boxlo, domain->prd);
+ }
+ if (!success)
+ error->one(FLERR,"Insufficient memory on accelerator");
+
+ if (host_start<inum) {
+ cpu_time = MPI_Wtime();
+ cpu_compute(host_start, inum, eflag, vflag, ilist, numneigh, firstneigh);
+ cpu_time = MPI_Wtime() - cpu_time;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairLJExpandCoulLongGPU::init_style()
+{
+ cut_respa = NULL;
+
+ if (!atom->q_flag)
+ error->all(FLERR,"Pair style lj/cut/coul/long/gpu requires atom attribute q");
+ if (force->newton_pair)
+ error->all(FLERR,"Cannot use newton pair with lj/cut/coul/long/gpu pair style");
+
+ // 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;
+
+ cut_coulsq = cut_coul * cut_coul;
+
+ // insure use of KSpace long-range solver, set g_ewald
+
+ if (force->kspace == NULL)
+ error->all(FLERR,"Pair style requires a KSpace style");
+ g_ewald = force->kspace->g_ewald;
+
+ // setup force tables
+
+ if (ncoultablebits) init_tables(cut_coul,cut_respa);
+
+ int maxspecial=0;
+ if (atom->molecular)
+ maxspecial=atom->maxspecial;
+ int success = ljecl_gpu_init(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4,
+ offset, shift, force->special_lj, atom->nlocal,
+ atom->nlocal+atom->nghost, 300, maxspecial,
+ cell_size, gpu_mode, screen, cut_ljsq, cut_coulsq,
+ force->special_coul, force->qqrd2e, g_ewald);
+ GPU_EXTRA::check_flag(success,error,world);
+
+ if (gpu_mode == GPU_FORCE) {
+ int irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->full = 1;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJExpandCoulLongGPU::reinit()
+{
+ Pair::reinit();
+
+ ljecl_gpu_reinit(atom->ntypes+1, cutsq, lj1, lj2, lj3, lj4, offset, shift, cut_ljsq);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLJExpandCoulLongGPU::memory_usage()
+{
+ double bytes = Pair::memory_usage();
+ return bytes + ljecl_gpu_bytes();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJExpandCoulLongGPU::cpu_compute(int start, int inum, int eflag,
+ int vflag, int *ilist, int *numneigh,
+ int **firstneigh)
+{
+ int i,j,ii,jj,jnum,itype,jtype,itable;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
+ double fraction,table;
+ double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
+ double grij,expm2,prefactor,t,erfc;
+ int *jlist;
+ double rsq;
+
+ evdwl = ecoul = 0.0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *type = atom->type;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ double qqrd2e = force->qqrd2e;
+
+ // 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];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ 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;
+
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ r = sqrt(rsq);
+ grij = g_ewald * r;
+ expm2 = exp(-grij*grij);
+ t = 1.0 / (1.0 + EWALD_P*grij);
+ erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+ prefactor = qqrd2e * qtmp*q[j]/r;
+ forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+ } else {
+ union_int_float_t rsq_lookup;
+ rsq_lookup.f = rsq;
+ itable = rsq_lookup.i & ncoulmask;
+ itable >>= ncoulshiftbits;
+ fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
+ table = ftable[itable] + fraction*dftable[itable];
+ forcecoul = qtmp*q[j] * table;
+ if (factor_coul < 1.0) {
+ table = ctable[itable] + fraction*dctable[itable];
+ prefactor = qtmp*q[j] * table;
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else forcecoul = 0.0;
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ } else forcelj = 0.0;
+
+ fpair = (forcecoul + factor_lj*forcelj) * r2inv;
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+
+ if (eflag) {
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq)
+ ecoul = prefactor*erfc;
+ else {
+ table = etable[itable] + fraction*detable[itable];
+ ecoul = qtmp*q[j] * table;
+ }
+ if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+ } else ecoul = 0.0;
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
+ offset[itype][jtype];
+ evdwl *= factor_lj;
+ } else evdwl = 0.0;
+ }
+
+ if (evflag) ev_tally_full(i,evdwl,ecoul,fpair,delx,dely,delz);
+ }
+ }
+ }
+}
diff --git a/src/GPU/pair_lj_expand_coul_long_gpu.h b/src/GPU/pair_lj_expand_coul_long_gpu.h
new file mode 100644
index 0000000000..4b39d8679b
--- /dev/null
+++ b/src/GPU/pair_lj_expand_coul_long_gpu.h
@@ -0,0 +1,67 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(lj/expand/coul/long/gpu,PairLJExpandCoulLongGPU)
+
+#else
+
+#ifndef LMP_PAIR_LJ_EXPAND_COUL_LONG_GPU_H
+#define LMP_PAIR_LJ_EXPAND_COUL_LONG_GPU_H
+
+#include "pair_lj_expand_coul_long.h"
+
+namespace LAMMPS_NS {
+
+class PairLJExpandCoulLongGPU : public PairLJExpandCoulLong {
+ public:
+ PairLJExpandCoulLongGPU(LAMMPS *lmp);
+ ~PairLJExpandCoulLongGPU();
+ void cpu_compute(int, int, int, int, int *, int *, int **);
+ void compute(int, int);
+ void init_style();
+ void reinit();
+ double memory_usage();
+
+ enum { GPU_FORCE, GPU_NEIGH, GPU_HYB_NEIGH };
+
+ private:
+ int gpu_mode;
+ double cpu_time;
+};
+
+}
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Insufficient memory on accelerator
+
+There is insufficient memory on one of the devices specified for the gpu
+package
+
+E: Pair style lj/cut/coul/long/gpu requires atom attribute q
+
+The atom style defined does not have this attribute.
+
+E: Cannot use newton pair with lj/cut/coul/long/gpu pair style
+
+Self-explanatory.
+
+E: Pair style requires a KSpace style
+
+No kspace style is defined.
+
+*/
diff --git a/src/USER-MISC/pair_lj_expand_coul_long.cpp b/src/USER-MISC/pair_lj_expand_coul_long.cpp
new file mode 100644
index 0000000000..c73581df15
--- /dev/null
+++ b/src/USER-MISC/pair_lj_expand_coul_long.cpp
@@ -0,0 +1,990 @@
+/* ----------------------------------------------------------------------
+ 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: Trung Nguyen (Northwestern)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "pair_lj_expand_coul_long.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "kspace.h"
+#include "update.h"
+#include "integrate.h"
+#include "respa.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+#define EWALD_F 1.12837917
+#define EWALD_P 0.3275911
+#define A1 0.254829592
+#define A2 -0.284496736
+#define A3 1.421413741
+#define A4 -1.453152027
+#define A5 1.061405429
+
+/* ---------------------------------------------------------------------- */
+
+PairLJExpandCoulLong::PairLJExpandCoulLong(LAMMPS *lmp) : Pair(lmp)
+{
+ ewaldflag = pppmflag = 1;
+ respa_enable = 1;
+ writedata = 1;
+ ftable = NULL;
+ qdist = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairLJExpandCoulLong::~PairLJExpandCoulLong()
+{
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+
+ memory->destroy(cut_lj);
+ memory->destroy(cut_ljsq);
+ memory->destroy(epsilon);
+ memory->destroy(sigma);
+ memory->destroy(lj1);
+ memory->destroy(lj2);
+ memory->destroy(lj3);
+ memory->destroy(lj4);
+ memory->destroy(offset);
+ memory->destroy(shift);
+ }
+ if (ftable) free_tables();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::compute(int eflag, int vflag)
+{
+ int i,ii,j,jj,inum,jnum,itype,jtype,itable;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
+ 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,rshift,rshiftsq;
+
+ evdwl = ecoul = 0.0;
+ 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;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+ double qqrd2e = force->qqrd2e;
+
+ 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];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ 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;
+
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ r = sqrt(rsq);
+ grij = g_ewald * r;
+ expm2 = exp(-grij*grij);
+ t = 1.0 / (1.0 + EWALD_P*grij);
+ erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+ prefactor = qqrd2e * qtmp*q[j]/r;
+ forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+ } else {
+ union_int_float_t rsq_lookup;
+ rsq_lookup.f = rsq;
+ itable = rsq_lookup.i & ncoulmask;
+ itable >>= ncoulshiftbits;
+ fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
+ table = ftable[itable] + fraction*dftable[itable];
+ forcecoul = qtmp*q[j] * table;
+ if (factor_coul < 1.0) {
+ table = ctable[itable] + fraction*dctable[itable];
+ prefactor = qtmp*q[j] * table;
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else forcecoul = 0.0;
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ r = sqrt(rsq);
+ rshift = r - shift[itype][jtype];
+ rshiftsq = rshift*rshift;
+ r2inv = 1.0/rshiftsq;
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ forcelj = factor_lj*forcelj/rshift/r;
+ } else forcelj = 0.0;
+
+ fpair = forcecoul*r2inv + forcelj;
+
+ 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 (eflag) {
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq)
+ ecoul = prefactor*erfc;
+ else {
+ table = etable[itable] + fraction*detable[itable];
+ ecoul = qtmp*q[j] * table;
+ }
+ if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+ } else ecoul = 0.0;
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
+ offset[itype][jtype];
+ evdwl *= factor_lj;
+ } else evdwl = 0.0;
+ }
+
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,ecoul,fpair,delx,dely,delz);
+ }
+ }
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::compute_inner()
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
+ double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
+ double rsw,r,rshift,rshiftsq;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+ double qqrd2e = force->qqrd2e;
+
+ inum = list->inum_inner;
+ ilist = list->ilist_inner;
+ numneigh = list->numneigh_inner;
+ firstneigh = list->firstneigh_inner;
+
+ double cut_out_on = cut_respa[0];
+ double cut_out_off = cut_respa[1];
+
+ double cut_out_diff = cut_out_off - cut_out_on;
+ double cut_out_on_sq = cut_out_on*cut_out_on;
+ double cut_out_off_sq = cut_out_off*cut_out_off;
+
+ // 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];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cut_out_off_sq) {
+ r2inv = 1.0/rsq;
+ forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
+
+ jtype = type[j];
+ if (rsq < cut_ljsq[itype][jtype]) {
+ r = sqrt(rsq);
+ rshift = r - shift[itype][jtype];
+ rshiftsq = rshift*rshift;
+ r2inv = 1.0/rshiftsq;
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ forcelj = factor_lj*forcelj/rshift/r;
+ } else forcelj = 0.0;
+
+ fpair = forcecoul*r2inv + forcelj;
+ if (rsq > cut_out_on_sq) {
+ rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+ fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
+ }
+
+ 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;
+ }
+ }
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::compute_middle()
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
+ double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
+ double rsw,r,rshift,rshiftsq;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+ double qqrd2e = force->qqrd2e;
+
+ inum = list->inum_middle;
+ ilist = list->ilist_middle;
+ numneigh = list->numneigh_middle;
+ firstneigh = list->firstneigh_middle;
+
+ double cut_in_off = cut_respa[0];
+ double cut_in_on = cut_respa[1];
+ double cut_out_on = cut_respa[2];
+ double cut_out_off = cut_respa[3];
+
+ double cut_in_diff = cut_in_on - cut_in_off;
+ double cut_out_diff = cut_out_off - cut_out_on;
+ double cut_in_off_sq = cut_in_off*cut_in_off;
+ double cut_in_on_sq = cut_in_on*cut_in_on;
+ double cut_out_on_sq = cut_out_on*cut_out_on;
+ double cut_out_off_sq = cut_out_off*cut_out_off;
+
+ // 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];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
+ r2inv = 1.0/rsq;
+ forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
+
+ jtype = type[j];
+ if (rsq < cut_ljsq[itype][jtype]) {
+ r = sqrt(rsq);
+ rshift = r - shift[itype][jtype];
+ rshiftsq = rshift*rshift;
+ r2inv = 1.0/rshiftsq;
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ forcelj = factor_lj*forcelj/rshift/r;
+ } else forcelj = 0.0;
+
+ fpair = forcecoul*r2inv + forcelj;
+ if (rsq < cut_in_on_sq) {
+ rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+ fpair *= rsw*rsw*(3.0 - 2.0*rsw);
+ }
+ if (rsq > cut_out_on_sq) {
+ rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+ fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
+ }
+
+ 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;
+ }
+ }
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::compute_outer(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype,itable;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
+ double fraction,table;
+ double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
+ double grij,expm2,prefactor,t,erfc;
+ double rsw,rshift,rshiftsq;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ double rsq;
+
+ evdwl = ecoul = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+ double qqrd2e = force->qqrd2e;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ double cut_in_off = cut_respa[2];
+ double cut_in_on = cut_respa[3];
+
+ double cut_in_diff = cut_in_on - cut_in_off;
+ double cut_in_off_sq = cut_in_off*cut_in_off;
+ double cut_in_on_sq = cut_in_on*cut_in_on;
+
+ // 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];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ 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;
+
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ r = sqrt(rsq);
+ grij = g_ewald * r;
+ expm2 = exp(-grij*grij);
+ t = 1.0 / (1.0 + EWALD_P*grij);
+ erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+ prefactor = qqrd2e * qtmp*q[j]/r;
+ forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0);
+ if (rsq > cut_in_off_sq) {
+ if (rsq < cut_in_on_sq) {
+ rsw = (r - cut_in_off)/cut_in_diff;
+ forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw);
+ if (factor_coul < 1.0)
+ forcecoul -=
+ (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw);
+ } else {
+ forcecoul += prefactor;
+ if (factor_coul < 1.0)
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else {
+ union_int_float_t rsq_lookup;
+ rsq_lookup.f = rsq;
+ itable = rsq_lookup.i & ncoulmask;
+ itable >>= ncoulshiftbits;
+ fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
+ table = ftable[itable] + fraction*dftable[itable];
+ forcecoul = qtmp*q[j] * table;
+ if (factor_coul < 1.0) {
+ table = ctable[itable] + fraction*dctable[itable];
+ prefactor = qtmp*q[j] * table;
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else forcecoul = 0.0;
+
+ if (rsq < cut_ljsq[itype][jtype] && rsq > cut_in_off_sq) {
+ r = sqrt(rsq);
+ rshift = r - shift[itype][jtype];
+ rshiftsq = rshift*rshift;
+ r2inv = 1.0/rshiftsq;
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ forcelj = factor_lj*forcelj/rshift/r;
+ if (rsq < cut_in_on_sq) {
+ rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+ forcelj *= rsw*rsw*(3.0 - 2.0*rsw);
+ }
+ } else forcelj = 0.0;
+
+ fpair = forcecoul*r2inv + forcelj;
+
+ 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 (eflag) {
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ ecoul = prefactor*erfc;
+ if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+ } else {
+ table = etable[itable] + fraction*detable[itable];
+ ecoul = qtmp*q[j] * table;
+ if (factor_coul < 1.0) {
+ table = ptable[itable] + fraction*dptable[itable];
+ prefactor = qtmp*q[j] * table;
+ ecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else ecoul = 0.0;
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ r6inv = r2inv*r2inv*r2inv;
+ evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
+ offset[itype][jtype];
+ evdwl *= factor_lj;
+ } else evdwl = 0.0;
+ }
+
+ if (vflag) {
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+ } else {
+ table = vtable[itable] + fraction*dvtable[itable];
+ forcecoul = qtmp*q[j] * table;
+ if (factor_coul < 1.0) {
+ table = ptable[itable] + fraction*dptable[itable];
+ prefactor = qtmp*q[j] * table;
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else forcecoul = 0.0;
+
+ if (rsq <= cut_in_off_sq) {
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ } else if (rsq <= cut_in_on_sq) {
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ }
+ fpair = (forcecoul + factor_lj*forcelj) * r2inv;
+ }
+
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,ecoul,fpair,delx,dely,delz);
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ memory->create(setflag,n+1,n+1,"pair:setflag");
+ for (int i = 1; i <= n; i++)
+ for (int j = i; j <= n; j++)
+ setflag[i][j] = 0;
+
+ memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+ memory->create(cut_lj,n+1,n+1,"pair:cut_lj");
+ memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq");
+ memory->create(epsilon,n+1,n+1,"pair:epsilon");
+ memory->create(sigma,n+1,n+1,"pair:sigma");
+ memory->create(lj1,n+1,n+1,"pair:lj1");
+ memory->create(lj2,n+1,n+1,"pair:lj2");
+ memory->create(lj3,n+1,n+1,"pair:lj3");
+ memory->create(lj4,n+1,n+1,"pair:lj4");
+ memory->create(offset,n+1,n+1,"pair:offset");
+ memory->create(shift,n+1,n+1,"pair:shift");
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::settings(int narg, char **arg)
+{
+ if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command");
+
+ cut_lj_global = force->numeric(FLERR,arg[0]);
+ if (narg == 1) cut_coul = cut_lj_global;
+ else cut_coul = force->numeric(FLERR,arg[1]);
+
+ // reset cutoffs that have been explicitly set
+
+ if (allocated) {
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++)
+ if (setflag[i][j]) cut_lj[i][j] = cut_lj_global;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::coeff(int narg, char **arg)
+{
+ if (narg < 5 || narg > 6)
+ error->all(FLERR,"Incorrect args for pair coefficients");
+ if (!allocated) allocate();
+
+ int ilo,ihi,jlo,jhi;
+ force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
+ force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
+
+ double epsilon_one = force->numeric(FLERR,arg[2]);
+ double sigma_one = force->numeric(FLERR,arg[3]);
+ double shift_one = force->numeric(FLERR,arg[4]);
+
+ double cut_lj_one = cut_lj_global;
+ if (narg == 6) cut_lj_one = force->numeric(FLERR,arg[5]);
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ epsilon[i][j] = epsilon_one;
+ sigma[i][j] = sigma_one;
+ shift[i][j] = shift_one;
+ cut_lj[i][j] = cut_lj_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::init_style()
+{
+ if (!atom->q_flag)
+ error->all(FLERR,"Pair style lj/cut/coul/long requires atom attribute q");
+
+ // request regular or rRESPA neighbor lists
+
+ int irequest;
+ int respa = 0;
+
+ if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
+ if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
+ if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
+ }
+
+ irequest = neighbor->request(this,instance_me);
+
+ if (respa >= 1) {
+ neighbor->requests[irequest]->respaouter = 1;
+ neighbor->requests[irequest]->respainner = 1;
+ }
+ if (respa == 2) neighbor->requests[irequest]->respamiddle = 1;
+
+ cut_coulsq = cut_coul * cut_coul;
+
+ // set rRESPA cutoffs
+
+ if (strstr(update->integrate_style,"respa") &&
+ ((Respa *) update->integrate)->level_inner >= 0)
+ cut_respa = ((Respa *) update->integrate)->cutoff;
+ else cut_respa = NULL;
+
+ // insure use of KSpace long-range solver, set g_ewald
+
+ if (force->kspace == NULL)
+ error->all(FLERR,"Pair style requires a KSpace style");
+ g_ewald = force->kspace->g_ewald;
+
+ // setup force tables
+
+ if (ncoultablebits) init_tables(cut_coul,cut_respa);
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairLJExpandCoulLong::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) {
+ epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
+ sigma[i][i],sigma[j][j]);
+ sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
+ cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]);
+ }
+
+ // include TIP4P qdist in full cutoff, qdist = 0.0 if not TIP4P
+
+ double cut = MAX(cut_lj[i][j],cut_coul+2.0*qdist);
+ cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
+
+ lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+ lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+ lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+ lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+
+ if (offset_flag && (cut_lj[i][j] > 0.0)) {
+ double ratio = sigma[i][j] / cut_lj[i][j];
+ offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
+ } else offset[i][j] = 0.0;
+
+ cut_ljsq[j][i] = cut_ljsq[i][j];
+ lj1[j][i] = lj1[i][j];
+ lj2[j][i] = lj2[i][j];
+ lj3[j][i] = lj3[i][j];
+ lj4[j][i] = lj4[i][j];
+ shift[j][i] = shift[i][j];
+ offset[j][i] = offset[i][j];
+
+ // check interior rRESPA cutoff
+
+ if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3])
+ error->all(FLERR,"Pair cutoff < Respa interior cutoff");
+
+ // compute I,J contribution to long-range tail correction
+ // count total # of atoms of type I and J via Allreduce
+
+ if (tail_flag) {
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+
+ double count[2],all[2];
+ count[0] = count[1] = 0.0;
+ for (int k = 0; k < nlocal; k++) {
+ if (type[k] == i) count[0] += 1.0;
+ if (type[k] == j) count[1] += 1.0;
+ }
+ MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
+
+ double sig2 = sigma[i][j]*sigma[i][j];
+ double sig6 = sig2*sig2*sig2;
+ double rc1 = cut_lj[i][j];
+ double rc2 = rc1*rc1;
+ double rc3 = rc2*rc1;
+ double rc6 = rc3*rc3;
+ double rc9 = rc3*rc6;
+ double shift1 = shift[i][j];
+ double shift2 = shift1*shift1;
+ double shift3 = shift2*shift1;
+ etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * sig6 *
+ ((1.0/9.0 + 2.0*shift1/(10.0*rc1) + shift2/(11.0*rc2))*sig6/rc9 -
+ (1.0/3.0 + 2.0*shift1/(4.0*rc1) + shift2/(5.0*rc2))/rc3);
+ ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] * sig6 *
+ ((1.0/9.0 + 3.0*shift1/(10.0*rc1) +
+ 3.0*shift2/(11.0*rc2) + shift3/(12.0*rc3))*2.0*sig6/rc9 -
+ (1.0/3.0 + 3.0*shift1/(4.0*rc1) +
+ 3.0*shift2/(5.0*rc2) + shift3/(6.0*rc3))/rc3);
+ }
+
+ return cut;
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::write_restart(FILE *fp)
+{
+ write_restart_settings(fp);
+
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ fwrite(&setflag[i][j],sizeof(int),1,fp);
+ if (setflag[i][j]) {
+ fwrite(&epsilon[i][j],sizeof(double),1,fp);
+ fwrite(&sigma[i][j],sizeof(double),1,fp);
+ fwrite(&shift[i][j],sizeof(double),1,fp);
+ fwrite(&cut_lj[i][j],sizeof(double),1,fp);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::read_restart(FILE *fp)
+{
+ read_restart_settings(fp);
+
+ allocate();
+
+ int i,j;
+ int me = comm->me;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+ MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+ if (setflag[i][j]) {
+ if (me == 0) {
+ fread(&epsilon[i][j],sizeof(double),1,fp);
+ fread(&sigma[i][j],sizeof(double),1,fp);
+ fread(&shift[i][j],sizeof(double),1,fp);
+ fread(&cut_lj[i][j],sizeof(double),1,fp);
+ }
+ MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&shift[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::write_restart_settings(FILE *fp)
+{
+ fwrite(&cut_lj_global,sizeof(double),1,fp);
+ fwrite(&cut_coul,sizeof(double),1,fp);
+ fwrite(&offset_flag,sizeof(int),1,fp);
+ fwrite(&mix_flag,sizeof(int),1,fp);
+ fwrite(&tail_flag,sizeof(int),1,fp);
+ fwrite(&ncoultablebits,sizeof(int),1,fp);
+ fwrite(&tabinner,sizeof(double),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::read_restart_settings(FILE *fp)
+{
+ if (comm->me == 0) {
+ fread(&cut_lj_global,sizeof(double),1,fp);
+ fread(&cut_coul,sizeof(double),1,fp);
+ fread(&offset_flag,sizeof(int),1,fp);
+ fread(&mix_flag,sizeof(int),1,fp);
+ fread(&tail_flag,sizeof(int),1,fp);
+ fread(&ncoultablebits,sizeof(int),1,fp);
+ fread(&tabinner,sizeof(double),1,fp);
+ }
+ MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
+ MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+ MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
+ MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world);
+ MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world);
+}
+
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::write_data(FILE *fp)
+{
+ for (int i = 1; i <= atom->ntypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],shift[i][i]);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes all pairs to data file
+------------------------------------------------------------------------- */
+
+void PairLJExpandCoulLong::write_data_all(FILE *fp)
+{
+ for (int i = 1; i <= atom->ntypes; i++)
+ for (int j = i; j <= atom->ntypes; j++)
+ fprintf(fp,"%d %d %g %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],shift[i][j],cut_lj[i][j]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLJExpandCoulLong::single(int i, int j, int itype, int jtype,
+ double rsq,
+ double factor_coul, double factor_lj,
+ double &fforce)
+{
+ double r2inv,r6inv,r,grij,expm2,t,erfc,prefactor;
+ double fraction,table,forcecoul,forcelj,phicoul,philj;
+ int itable;
+
+ r2inv = 1.0/rsq;
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ r = sqrt(rsq);
+ grij = g_ewald * r;
+ expm2 = exp(-grij*grij);
+ t = 1.0 / (1.0 + EWALD_P*grij);
+ erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+ prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
+ forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+ } else {
+ union_int_float_t rsq_lookup_single;
+ rsq_lookup_single.f = rsq;
+ itable = rsq_lookup_single.i & ncoulmask;
+ itable >>= ncoulshiftbits;
+ fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
+ table = ftable[itable] + fraction*dftable[itable];
+ forcecoul = atom->q[i]*atom->q[j] * table;
+ if (factor_coul < 1.0) {
+ table = ctable[itable] + fraction*dctable[itable];
+ prefactor = atom->q[i]*atom->q[j] * table;
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else forcecoul = 0.0;
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ } else forcelj = 0.0;
+
+ fforce = (forcecoul + factor_lj*forcelj) * r2inv;
+
+ double eng = 0.0;
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq)
+ phicoul = prefactor*erfc;
+ else {
+ table = etable[itable] + fraction*detable[itable];
+ phicoul = atom->q[i]*atom->q[j] * table;
+ }
+ if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
+ eng += phicoul;
+ }
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
+ offset[itype][jtype];
+ eng += factor_lj*philj;
+ }
+
+ return eng;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void *PairLJExpandCoulLong::extract(const char *str, int &dim)
+{
+ dim = 0;
+ if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
+ dim = 2;
+ if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
+ if (strcmp(str,"sigma") == 0) return (void *) sigma;
+ if (strcmp(str,"delta") == 0) return (void *) shift;
+ return NULL;
+}
diff --git a/src/USER-MISC/pair_lj_expand_coul_long.h b/src/USER-MISC/pair_lj_expand_coul_long.h
new file mode 100644
index 0000000000..d88acb99b7
--- /dev/null
+++ b/src/USER-MISC/pair_lj_expand_coul_long.h
@@ -0,0 +1,93 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(lj/expand/coul/long,PairLJExpandCoulLong)
+
+#else
+
+#ifndef LMP_PAIR_LJ_EXPAND_COUL_LONG_H
+#define LMP_PAIR_LJ_EXPAND_COUL_LONG_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairLJExpandCoulLong : public Pair {
+
+ public:
+ PairLJExpandCoulLong(class LAMMPS *);
+ virtual ~PairLJExpandCoulLong();
+ virtual void compute(int, int);
+ virtual void settings(int, char **);
+ void coeff(int, char **);
+ virtual void init_style();
+ virtual double init_one(int, int);
+ void write_restart(FILE *);
+ void read_restart(FILE *);
+ virtual void write_restart_settings(FILE *);
+ virtual void read_restart_settings(FILE *);
+ void write_data(FILE *);
+ void write_data_all(FILE *);
+ virtual double single(int, int, int, int, double, double, double, double &);
+
+ void compute_inner();
+ void compute_middle();
+ virtual void compute_outer(int, int);
+ virtual void *extract(const char *, int &);
+
+ protected:
+ double cut_lj_global;
+ double **cut_lj,**cut_ljsq;
+ double cut_coul,cut_coulsq;
+ double **epsilon,**sigma;
+ double **lj1,**lj2,**lj3,**lj4,**offset,**shift;
+ double *cut_respa;
+ double qdist; // TIP4P distance from O site to negative charge
+ double g_ewald;
+
+ virtual void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: Pair style lj/cut/coul/long requires atom attribute q
+
+The atom style defined does not have this attribute.
+
+E: Pair style requires a KSpace style
+
+No kspace style is defined.
+
+E: Pair cutoff < Respa interior cutoff
+
+One or more pairwise cutoffs are too short to use with the specified
+rRESPA cutoffs.
+
+*/
--
GitLab