External Field

include/ff/amoeba/epolar.h

@@ -10,5 +10,6 @@ void epolar(int vers);
 void epolarEwaldRecipSelf(int vers);
 // see also subroutine epolar0e in epolar.f
 void epolar0DotProd(const real (*uind)[3], const real (*udirp)[3]);
+void epolarPairwiseExtfield(int vers, const real (*uind)[3]);
 /// \}
 }

include/ff/elec.h

@@ -66,6 +66,12 @@ namespace tinker {
 
 bool useEwald();
 void elecData(RcOp);
+void exfield(int vers,     ///< Common integer flag for the energy components.
+             int useDipole ///< If 0, use partial charge; otherwise, also include dipole.
+);
+void extfieldModifyDField(real (*field)[3], ///< Permanent field.
+                          real (*fieldp)[3] ///< Set to \c nullptr if not using AMOEBA.
+);
 
 //====================================================================//
 //                                                                    //

src/acc/amoeba/empolenonewald.cpp

@@ -1,10 +1,11 @@
-#include "ff/modamoeba.h"
 #include "ff/energy.h"
 #include "ff/image.h"
+#include "ff/modamoeba.h"
 #include "ff/nblist.h"
 #include "ff/switch.h"
 #include "seq/pair_mpole.h"
 #include "tool/gpucard.h"
+#include <tinker/detail/extfld.hh>
 
 namespace tinker {
 #define DEVICE_PTRS x, y, z, demx, demy, demz, rpole, nem, em, vir_em, trqx, trqy, trqz
@@ -65,9 +66,8 @@ void empoleNonEwald_acc1()
             pair_mpole<do_e, do_g, NON_EWALD>(                        //
                r2, xr, yr, zr, 1,                                     //
                ci, dix, diy, diz, qixx, qixy, qixz, qiyy, qiyz, qizz, //
-               rpole[k][MPL_PME_0], rpole[k][MPL_PME_X], rpole[k][MPL_PME_Y], rpole[k][MPL_PME_Z],
-               rpole[k][MPL_PME_XX], rpole[k][MPL_PME_XY], rpole[k][MPL_PME_XZ],
-               rpole[k][MPL_PME_YY], rpole[k][MPL_PME_YZ],
+               rpole[k][MPL_PME_0], rpole[k][MPL_PME_X], rpole[k][MPL_PME_Y], rpole[k][MPL_PME_Z], rpole[k][MPL_PME_XX],
+               rpole[k][MPL_PME_XY], rpole[k][MPL_PME_XZ], rpole[k][MPL_PME_YY], rpole[k][MPL_PME_YZ],
                rpole[k][MPL_PME_ZZ], //
                f, 0, e, pgrad);
 
@@ -151,9 +151,8 @@ void empoleNonEwald_acc1()
          pair_mpole<do_e, do_g, NON_EWALD>(                        //
             r2, xr, yr, zr, mscale,                                //
             ci, dix, diy, diz, qixx, qixy, qixz, qiyy, qiyz, qizz, //
-            rpole[k][MPL_PME_0], rpole[k][MPL_PME_X], rpole[k][MPL_PME_Y], rpole[k][MPL_PME_Z],
-            rpole[k][MPL_PME_XX], rpole[k][MPL_PME_XY], rpole[k][MPL_PME_XZ], rpole[k][MPL_PME_YY],
-            rpole[k][MPL_PME_YZ],
+            rpole[k][MPL_PME_0], rpole[k][MPL_PME_X], rpole[k][MPL_PME_Y], rpole[k][MPL_PME_Z], rpole[k][MPL_PME_XX],
+            rpole[k][MPL_PME_XY], rpole[k][MPL_PME_XZ], rpole[k][MPL_PME_YY], rpole[k][MPL_PME_YZ],
             rpole[k][MPL_PME_ZZ], //
             f, 0, e, pgrad);
 
@@ -209,4 +208,83 @@ void empoleNonEwald_acc(int vers)
    else if (vers == calc::v6)
       empoleNonEwald_acc1<calc::V6>();
 }
+
+void exfieldDipole_acc(int vers)
+{
+   bool do_e = vers & calc::energy;
+   bool do_a = vers & calc::analyz;
+   bool do_g = vers & calc::grad;
+   bool do_v = vers & calc::virial;
+
+   auto bufsize = bufferSize();
+   real f = electric / dielec;
+   real ef1 = extfld::exfld[0], ef2 = extfld::exfld[1], ef3 = extfld::exfld[2];
+
+   #pragma acc parallel async deviceptr(rpole,nem,em,vir_em,x,y,z,demx,demy,demz,trqx,trqy,trqz)
+   #pragma acc loop independent
+   for (int ii = 0; ii < n; ++ii) {
+      int offset = ii & (bufsize - 1);
+      real xi = x[ii], yi = y[ii], zi = z[ii];
+      real ci = rpole[ii][0], dix = rpole[ii][1], diy = rpole[ii][2], diz = rpole[ii][3];
+
+      if (do_e) {
+         real phi = xi * ef1 + yi * ef2 + zi * ef3; // negative potential
+         real e = -f * (ci * phi + dix * ef1 + diy * ef2 + diz * ef3);
+         atomic_add(e, em, offset);
+         if (do_a)
+            atomic_add(1, nem, offset);
+      }
+      if (do_g) {
+         // torque due to the dipole
+         real tx = f * (diy * ef3 - diz * ef2);
+         real ty = f * (diz * ef1 - dix * ef3);
+         real tz = f * (dix * ef2 - diy * ef1);
+         atomic_add(tx, trqx, ii);
+         atomic_add(ty, trqy, ii);
+         atomic_add(tz, trqz, ii);
+         // gradient and virial due to the monopole
+         real frx = -f * ef1 * ci;
+         real fry = -f * ef2 * ci;
+         real frz = -f * ef3 * ci;
+         atomic_add(frx, demx, ii);
+         atomic_add(fry, demy, ii);
+         atomic_add(frz, demz, ii);
+         if (do_v) {
+            real vxx = xi * frx;
+            real vyy = yi * fry;
+            real vzz = zi * frz;
+            real vxy = (yi * frx + xi * fry) / 2;
+            real vxz = (zi * frx + xi * frz) / 2;
+            real vyz = (zi * fry + yi * frz) / 2;
+            atomic_add(vxx, vxy, vxz, vyy, vyz, vzz, vir_em, offset);
+         }
+      }
+   }
+}
+
+void extfieldModifyDField_acc(real (*field)[3], real (*fieldp)[3])
+{
+   real ex1 = extfld::exfld[0];
+   real ex2 = extfld::exfld[1];
+   real ex3 = extfld::exfld[2];
+
+   if (fieldp) {
+      #pragma acc parallel loop independent async deviceptr(field,fieldp)
+      for (int i = 0; i < n; ++i) {
+         field[i][0] += ex1;
+         field[i][1] += ex2;
+         field[i][2] += ex3;
+         fieldp[i][0] += ex1;
+         fieldp[i][1] += ex2;
+         fieldp[i][2] += ex3;
+      }
+   } else {
+      #pragma acc parallel loop independent async deviceptr(field)
+      for (int i = 0; i < n; ++i) {
+         field[i][0] += ex1;
+         field[i][1] += ex2;
+         field[i][2] += ex3;
+      }
+   }
+}
 }

src/acc/amoeba/epolarnonewald.cpp

@@ -5,6 +5,7 @@
 #include "ff/switch.h"
 #include "seq/pair_polar.h"
 #include "tool/gpucard.h"
+#include <tinker/detail/extfld.hh>
 
 #define TINKER9_POLPAIR 2
 
@@ -25,6 +26,23 @@ void epolar0DotProd_acc(const real (*gpu_uind)[3], const real (*gpu_udirp)[3])
       atomic_add(f * e, ep, offset);
    }
 }
+
+void epolarPairwiseExtfield_acc(const real (*uind)[3]) {
+   const real f = -0.5 * electric / dielec;
+
+   auto bufsize = bufferSize();
+
+   real ex1 = extfld::exfld[0];
+   real ex2 = extfld::exfld[1];
+   real ex3 = extfld::exfld[2];
+
+   #pragma acc parallel loop independent async deviceptr(ep,uind)
+   for (int i = 0; i < n; ++i) {
+      int offset = i & (bufsize - 1);
+      real e = uind[i][0] * ex1 + uind[i][1] * ex2 + uind[i][2] * ex3;
+      atomic_add(f * e, ep, offset);
+   }
+}
 }
 
 namespace tinker {

src/acc/echarge.cpp

@@ -1,5 +1,5 @@
-#include "ff/echarge.h"
 #include "ff/atom.h"
+#include "ff/echarge.h"
 #include "ff/image.h"
 #include "ff/nblist.h"
 #include "ff/pme.h"
@@ -8,6 +8,7 @@
 #include "seq/bsplgen.h"
 #include "seq/pair_charge.h"
 #include "tool/gpucard.h"
+#include <tinker/detail/extfld.hh>
 
 namespace tinker {
 #define DEVICE_PTRS x, y, z, decx, decy, decz, pchg, nec, ec, vir_ec
@@ -367,4 +368,48 @@ void echargeEwaldFphiSelf_acc(int vers)
    else if (vers == calc::v6)
       echarge_acc3<calc::V6, 5>();
 }
+
+void exfieldCharge_acc(int vers)
+{
+   bool do_e = vers & calc::energy;
+   bool do_a = vers & calc::analyz;
+   bool do_g = vers & calc::grad;
+   bool do_v = vers & calc::virial;
+
+   auto bufsize = bufferSize();
+   real f = electric / dielec;
+   real ef1 = extfld::exfld[0], ef2 = extfld::exfld[1], ef3 = extfld::exfld[2];
+
+   #pragma acc parallel async deviceptr(pchg,nec,ec,vir_ec,x,y,z,decx,decy,decz)
+   #pragma acc loop independent
+   for (int ii = 0; ii < n; ++ii) {
+      int offset = ii & (bufsize - 1);
+      real xi = x[ii], yi = y[ii], zi = z[ii], ci = pchg[ii];
+
+      if (do_e) {
+         real phi = xi * ef1 + yi * ef2 + zi * ef3; // negative potential
+         real e = -f * ci * phi;
+         atomic_add(e, ec, offset);
+         if (do_a)
+            atomic_add(1, nec, offset);
+      }
+      if (do_g) {
+         real frx = -f * ef1 * ci;
+         real fry = -f * ef2 * ci;
+         real frz = -f * ef3 * ci;
+         atomic_add(frx, decx, ii);
+         atomic_add(fry, decy, ii);
+         atomic_add(frz, decz, ii);
+         if (do_v) {
+            real vxx = xi * frx;
+            real vyy = yi * fry;
+            real vzz = zi * frz;
+            real vxy = (yi * frx + xi * fry) / 2;
+            real vxz = (zi * frx + xi * frz) / 2;
+            real vyz = (zi * fry + yi * frz) / 2;
+            atomic_add(vxx, vxy, vxz, vyy, vyz, vzz, vir_ec, offset);
+         }
+      }
+   }
+}
 }

src/amoeba/emplar.cpp

@@ -1,6 +1,7 @@
 #include "ff/amoeba/empole.h"
-#include "ff/modamoeba.h"
+#include "ff/elec.h"
 #include "ff/energy.h"
+#include "ff/modamoeba.h"
 #include "math/zero.h"
 #include "tool/externfunc.h"
 
@@ -14,6 +15,8 @@ void emplar(int vers)
 
    mpoleInit(vers);
    TINKER_FCALL2(acc0, cu1, emplar, vers);
+   exfield(vers, 1);
+   // epolarPairwiseExtfield(vers, uind); // emplar uses the dot product version
    torque(vers, demx, demy, demz);
    if (do_v) {
       VirialBuffer u2 = vir_trq;

src/amoeba/empole.cpp

@@ -104,6 +104,7 @@ void empole(int vers)
       empoleEwald(vers);
    else
       empoleNonEwald(vers);
+   exfield(vers, 1);
    torque(vers, demx, demy, demz);
    if (do_v) {
       VirialBuffer u2 = vir_trq;

src/amoeba/epolar.cpp

@@ -13,6 +13,7 @@
 #include "tool/iofortstr.h"
 #include "tool/ioprint.h"
 #include <tinker/detail/couple.hh>
+#include <tinker/detail/extfld.hh>
 #include <tinker/detail/mplpot.hh>
 #include <tinker/detail/polar.hh>
 #include <tinker/detail/polgrp.hh>
@@ -616,6 +617,7 @@ void epolar(int vers)
       else
          epolarNonEwald(vers);
    }
+   epolarPairwiseExtfield(vers, uind);
    torque(vers, depx, depy, depz);
    if (use_cfgrad) dcflux(vers, depx, depy, depz, vir_ep);
    if (do_v) {
@@ -653,4 +655,11 @@ void epolar0DotProd(const real (*uind)[3], const real (*udirp)[3])
 {
    TINKER_FCALL2(acc1, cu1, epolar0DotProd, uind, udirp);
 }
+
+TINKER_FVOID2(acc1, cu1, epolarPairwiseExtfield, const real (*)[3]);
+void epolarPairwiseExtfield(int vers, const real (*uind)[3]) {
+   if (extfld::use_exfld and (vers & calc::analyz)) {
+      TINKER_FCALL2(acc1, cu1, epolarPairwiseExtfield, uind);
+   }
+}
 }

src/amoeba/field.cpp

@@ -59,6 +59,7 @@ void dfield(real (*field)[3], real (*fieldp)[3])
       dfieldEwald(field, fieldp);
    else
       dfieldNonEwald(field, fieldp);
+   extfieldModifyDField(field, fieldp);
 }
 }