prepare for issue #3

.gitignore

@@ -15,3 +15,6 @@ build*/*
 
 # vi
 *.swp
+
+# gdb
+*.out

example/dhfr2.key

@@ -0,0 +1,15 @@
+echo ################################################################
+echo Joint Amber-CHARMM Benchmark on Dihydrofolate Reductase in Water
+echo 23558 Atoms, 62.23 Ang Cube, 9 Ang Nonbond Cutoffs, 64x64x64 PME
+echo ################################################################
+
+parameters    ../params/amoebabio09
+neighbor-list
+a-axis                        62.23
+vdw-cutoff                     12.0
+ewald
+ewald-cutoff                    7.0
+pme-grid                   64 64 64
+pme-order                         5
+polarization                 MUTUAL
+polar-eps                   0.00001

include/acc_add.h

@@ -0,0 +1,38 @@
+#ifndef TINKER_ACC_ADD_H_
+#define TINKER_ACC_ADD_H_
+
+#include "macro.h"
+#include <type_traits>
+
+TINKER_NAMESPACE_BEGIN
+/**
+ * @brief
+ * add @c value to @c buffer[@c offset] atomically
+ */
+#pragma acc routine seq
+template <class T>
+inline void atomic_add_value(T value, T* buffer, int offset = 0) {
+  #pragma acc atomic update
+  buffer[offset] += value;
+}
+
+/**
+ * @brief
+ * add @c value to @c buffer[@c offset] atomically via fixed-point arithmetic
+ *
+ * @tparam T
+ * must be a floating point type
+ */
+#pragma acc routine seq
+template <
+    class T,
+    class = typename std::enable_if<std::is_floating_point<T>::value>::type>
+inline void atomic_add_value(T value, unsigned long long* buffer,
+                             int offset = 0) {
+  #pragma acc atomic update
+  buffer[offset] += static_cast<unsigned long long>(
+      static_cast<long long>(value * fixed_point));
+}
+TINKER_NAMESPACE_END
+
+#endif

include/acc_image.h

@@ -0,0 +1,26 @@
+#ifndef TINKER_ACC_IMAGE_H_
+#define TINKER_ACC_IMAGE_H_
+
+#include "box.h"
+
+TINKER_NAMESPACE_BEGIN
+/**
+ * @brief
+ * apply periodic boundary conditions to displacement (@c xr, @c yr, @c zr) and
+ * preserve the correct signs
+ */
+#pragma acc routine seq
+void image(real& __restrict__ xr, real& __restrict__ yr, real& __restrict__ zr,
+           const Box* __restrict__ pb);
+
+/**
+ * @brief
+ * apply periodic boundary conditions to displacement (@c xr, @c yr, @c zr) but
+ * only guarantee the lengths are correct
+ */
+#pragma acc routine seq
+void imagen(real& __restrict__ xr, real& __restrict__ yr, real& __restrict__ zr,
+            const Box* __restrict__ pb);
+TINKER_NAMESPACE_END
+
+#endif

include/acc_mathfunc.h

@@ -0,0 +1,38 @@
+#ifndef TINKER_ACC_MATHFUNC_H_
+#define TINKER_ACC_MATHFUNC_H_
+
+#include "mathfunc.h"
+
+#pragma acc routine(abs) seq
+
+#pragma acc routine(sqrt) seq
+#pragma acc routine(exp) seq
+#pragma acc routine(floor) seq
+#pragma acc routine(fabs) seq
+#pragma acc routine(pow) seq
+#pragma acc routine(cos) seq
+#pragma acc routine(sin) seq
+#pragma acc routine(acos) seq
+#pragma acc routine(asin) seq
+#pragma acc routine(erf) seq
+#pragma acc routine(erfc) seq
+#pragma acc routine(fmin) seq
+#pragma acc routine(fmax) seq
+#pragma acc routine(copysign) seq
+
+#pragma acc routine(sqrtf) seq
+#pragma acc routine(expf) seq
+#pragma acc routine(floorf) seq
+#pragma acc routine(fabsf) seq
+#pragma acc routine(powf) seq
+#pragma acc routine(cosf) seq
+#pragma acc routine(sinf) seq
+#pragma acc routine(acosf) seq
+#pragma acc routine(asinf) seq
+#pragma acc routine(erff) seq
+#pragma acc routine(erfcf) seq
+#pragma acc routine(fminf) seq
+#pragma acc routine(fmaxf) seq
+#pragma acc routine(copysignf) seq
+
+#endif

include/acc_switch.h

@@ -0,0 +1,49 @@
+#ifndef TINKER_ACC_SWITCH_H_
+#define TINKER_ACC_SWITCH_H_
+
+#include "macro.h"
+
+TINKER_NAMESPACE_BEGIN
+/**
+ * @brief
+ * second order smooth step function
+ * @f[
+ * f: [cut,off]\rightarrow[1,0]
+ * @f]
+ *
+ * deriving from
+ * @f[ S_2: [0,1]\rightarrow[0,1] @f]
+ * @f[ S_2(x) = 6 x^5 - 15 x^4 + 10 x^3 @f]
+ *
+ * @param[in] rik
+ * distance
+ *
+ * @param[in] cut
+ * distance at which switching of the potential begins
+ *
+ * @param[out] off
+ * distance at which the potential energy goes to zero
+ *
+ * @param[out] taper
+ * @f$ f @f$ value
+ *
+ * @param[out] dtaper
+ * @f$ df/dx @f$ value
+ *
+ * @tparam DO_DTAPER
+ * if @c false, @c dtaper will not be calculated
+ */
+#pragma acc routine seq
+template <int DO_DTAPER>
+inline void switch_taper5(real rik, real cut, real off,
+                          real& __restrict__ taper, real& __restrict__ dtaper) {
+  real _1_ab = REAL_RECIP(cut - off);
+  real x = (rik - off) * _1_ab;
+  real x2 = x * x;
+  real x3 = x2 * x;
+  taper = x3 * (6 * x2 - 15 * x + 10);
+  if_constexpr(DO_DTAPER) { dtaper = 30 * REAL_SQ(x * (1 - x)) * _1_ab; }
+}
+TINKER_NAMESPACE_END
+
+#endif

include/array.h

@@ -9,6 +9,7 @@ TINKER_NAMESPACE_BEGIN
 /// zero out the first n elements of an array on device
 /// @{
 void zero_array(int* dst, int nelem);
+void zero_array(unsigned long long* dst, int nelem);
 void zero_array(float* dst, int nelem);
 void zero_array(double* dst, int nelem);
 /// @}
@@ -24,6 +25,11 @@ void zero_array(double* dst, int nelem);
 void copyin_array(int* dst, const int* src, int nelem);
 void copyout_array(int* dst, const int* src, int nelem);
 
+void copyin_array(unsigned long long* dst, const unsigned long long* src,
+                  int nelem);
+void copyout_array(unsigned long long* dst, const unsigned long long* src,
+                   int nelem);
+
 void copyin_array(float* dst, const float* src, int nelem);
 void copyout_array(float* dst, const float* src, int nelem);
 
@@ -35,13 +41,8 @@ void copyout_array(double* dst, const double* src, int nelem);
 /// @}
 
 /// @brief
-/// copy the @c idx0-th of every @c ndim elements from @c src to @c dst
-///
-/// @param[in] idx0
-/// ranges from 0 to ndim-1
-///
-/// @param[in] nelem
-/// number of elements copied to @c dst
+/// pick the @c idx0-th of every @c ndim elements from @c src to @c dst of
+/// length @c nelem
 /// @{
 void copyin_array2(int idx0, int ndim, float* dst, const float* src, int nelem);
 void copyout_array2(int idx0, int ndim, float* dst, const float* src,

include/couple.h

@@ -23,7 +23,7 @@ struct Couple {
   ~Couple();
 };
 
-typedef GenericUnit<Couple, GenericUnitVersion::V1> CoupleUnit;
+typedef GenericUnit<Couple, GenericUnitVersion::EnableOnDevice> CoupleUnit;
 TINKER_EXTERN CoupleUnit couple_unit;
 
 void couple_data(rc_op);

include/e_angle.h

@@ -1,6 +1,7 @@
 #ifndef TINKER_E_ANGLE_H_
 #define TINKER_E_ANGLE_H_
 
+#include "energy_buffer.h"
 #include "rc_man.h"
 
 TINKER_NAMESPACE_BEGIN
@@ -16,8 +17,7 @@ TINKER_EXTERN real angunit;
 TINKER_EXTERN real cang, qang, pang, sang;
 TINKER_EXTERN eangle_t* angtyp;
 
-TINKER_EXTERN real* ea;
-TINKER_EXTERN real* vir_ea;
+TINKER_EXTERN BondedEnergy ea_handle;
 
 void eangle_data(rc_op op);
 

include/e_bond.h

@@ -1,6 +1,7 @@
 #ifndef TINKER_E_BOND_H_
 #define TINKER_E_BOND_H_
 
+#include "energy_buffer.h"
 #include "rc_man.h"
 
 TINKER_NAMESPACE_BEGIN
@@ -12,8 +13,7 @@ TINKER_EXTERN int nbond;
 TINKER_EXTERN int (*ibnd)[2];
 TINKER_EXTERN real *bl, *bk;
 
-TINKER_EXTERN real* eb;
-TINKER_EXTERN real* vir_eb;
+TINKER_EXTERN BondedEnergy eb_handle;
 
 void ebond_data(rc_op op);
 

include/e_mpole.h

@@ -2,15 +2,14 @@
 #define TINKER_E_MPOLE_H_
 
 #include "elec.h"
+#include "energy_buffer.h"
 
 TINKER_NAMESPACE_BEGIN
 TINKER_EXTERN elec_t empole_electyp;
 
 TINKER_EXTERN real m2scale, m3scale, m4scale, m5scale;
 
-TINKER_EXTERN real* em;
-TINKER_EXTERN int* nem;
-TINKER_EXTERN real* vir_em;
+TINKER_EXTERN NonbondedEnergy em_handle;
 
 void empole_data(rc_op op);
 

include/e_opbend.h

@@ -1,6 +1,7 @@
 #ifndef TINKER_E_OPBEND_H_
 #define TINKER_E_OPBEND_H_
 
+#include "energy_buffer.h"
 #include "rc_man.h"
 
 TINKER_NAMESPACE_BEGIN
@@ -13,8 +14,7 @@ TINKER_EXTERN real* opbk;
 TINKER_EXTERN real opbunit;
 TINKER_EXTERN real copb, qopb, popb, sopb;
 
-TINKER_EXTERN real* eopb;
-TINKER_EXTERN real* vir_eopb;
+TINKER_EXTERN BondedEnergy eopb_handle;
 
 void eopbend_data(rc_op op);