From caf19462bf18583f4afa19f2e79e68153158f19d Mon Sep 17 00:00:00 2001 From: Zhi Wang Date: Sat, 25 Feb 2023 00:28:57 -0600 Subject: [PATCH 1/4] start working on extfield --- include/ff/elec.h | 3 + src/acc/amoeba/empolenonewald.cpp | 66 +++++++++++-- src/acc/echarge.cpp | 47 +++++++++- src/amoeba/empole.cpp | 1 + src/elec.cpp | 148 +++++++++++++++++++++--------- test/cmakesrc.txt | 1 + test/extfield.cpp | 121 ++++++++++++++++++++++++ test/file/extfield/water4.key | 14 +++ test/file/extfield/water4.xyz | 13 +++ test/ref/extfield.1.txt | 25 +++++ test/ref/extfield.2.txt | 0 test/ref/extfield.3.txt | 0 12 files changed, 388 insertions(+), 51 deletions(-) create mode 100644 test/extfield.cpp create mode 100644 test/file/extfield/water4.key create mode 100644 test/file/extfield/water4.xyz create mode 100644 test/ref/extfield.1.txt create mode 100644 test/ref/extfield.2.txt create mode 100644 test/ref/extfield.3.txt diff --git a/include/ff/elec.h b/include/ff/elec.h index c0543cb01..61b887021 100644 --- a/include/ff/elec.h +++ b/include/ff/elec.h @@ -66,6 +66,9 @@ 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. +); //====================================================================// // // diff --git a/src/acc/amoeba/empolenonewald.cpp b/src/acc/amoeba/empolenonewald.cpp index 265af8dfa..00fc96447 100644 --- a/src/acc/amoeba/empolenonewald.cpp +++ b/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 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( // 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( // 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,57 @@ void empoleNonEwald_acc(int vers) else if (vers == calc::v6) empoleNonEwald_acc1(); } + +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; + real vxz = zi * frx + xi * frz; + real vyz = zi * fry + yi * frz; + atomic_add(vxx, vxy, vxz, vyy, vyz, vzz, vir_em, offset); + } + } + } +} } diff --git a/src/acc/echarge.cpp b/src/acc/echarge.cpp index 630fd675e..6d6678b35 100644 --- a/src/acc/echarge.cpp +++ b/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 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(); } + +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; + real vxz = zi * frx + xi * frz; + real vyz = zi * fry + yi * frz; + atomic_add(vxx, vxy, vxz, vyy, vyz, vzz, vir_ec, offset); + } + } + } +} } diff --git a/src/amoeba/empole.cpp b/src/amoeba/empole.cpp index 6d73d1c5e..59ce9b095 100644 --- a/src/amoeba/empole.cpp +++ b/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; diff --git a/src/elec.cpp b/src/elec.cpp index 5fb8b1ca6..efbf48b65 100644 --- a/src/elec.cpp +++ b/src/elec.cpp @@ -1,7 +1,7 @@ #include "ff/elec.h" -#include "ff/modamoeba.h" #include "ff/echarge.h" #include "ff/energy.h" +#include "ff/modamoeba.h" #include "ff/modhippo.h" #include "ff/potent.h" #include "tool/iofortstr.h" @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -20,11 +21,16 @@ namespace tinker { static void pchgData(RcOp op) { - if (not use(Potent::CHARGE)) return; + if (not use(Potent::CHARGE)) + return; - if (op & RcOp::DEALLOC) { darray::deallocate(pchg); } + if (op & RcOp::DEALLOC) { + darray::deallocate(pchg); + } - if (op & RcOp::ALLOC) { darray::allocate(n, &pchg); } + if (op & RcOp::ALLOC) { + darray::allocate(n, &pchg); + } if (op & RcOp::INIT) { std::vector pchgbuf(n); @@ -46,8 +52,7 @@ static void pchgData(RcOp op) TINKER_FVOID2(cpp0, cu1, mpoleDataBinding, RcOp); static void mpoleData(RcOp op) { - if (not use(Potent::MPOLE) and not use(Potent::POLAR) - and not use(Potent::REPULS)) + if (not use(Potent::MPOLE) and not use(Potent::POLAR) and not use(Potent::REPULS)) return; if (op & RcOp::DEALLOC) { @@ -140,7 +145,8 @@ static void mpoleData(RcOp op) static void mdpuscaleData(RcOp op) { - if (not use(Potent::MPOLE) and not use(Potent::POLAR)) return; + if (not use(Potent::MPOLE) and not use(Potent::POLAR)) + return; if (op & RcOp::DEALLOC) { nmexclude = 0; @@ -156,8 +162,7 @@ static void mdpuscaleData(RcOp op) real m; }; std::map ikm; - auto insert_m = [](std::map& a, int i, int k, real val, - char ch) { + auto insert_m = [](std::map& a, int i, int k, real val, char ch) { key_t key; key.first = i; key.second = k; @@ -165,10 +170,12 @@ static void mdpuscaleData(RcOp op) if (it == a.end()) { m x; x.m = 1; - if (ch == 'm') x.m = val; + if (ch == 'm') + x.m = val; a[key] = x; } else { - if (ch == 'm') it->second.m = val; + if (ch == 'm') + it->second.m = val; } }; struct mdpu @@ -176,8 +183,7 @@ static void mdpuscaleData(RcOp op) real m, d, p, u; }; std::map ik_scale; - auto insert_mdpu = [](std::map& a, int i, int k, real val, - char ch) { + auto insert_mdpu = [](std::map& a, int i, int k, real val, char ch) { key_t key; key.first = i; key.second = k; @@ -310,10 +316,13 @@ static void mdpuscaleData(RcOp op) int k = couple_i12[bask + j]; real val = p2scale; for (int jj = 0; jj < polgrp::np11[i]; ++jj) { - if (k == polgrp::ip11[i * maxp11 + jj]) val = p2iscale; + if (k == polgrp::ip11[i * maxp11 + jj]) + val = p2iscale; } k -= 1; - if (k > i) { insert_mdpu(ik_scale, i, k, val, 'p'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, val, 'p'); + } } } @@ -324,10 +333,13 @@ static void mdpuscaleData(RcOp op) int k = couple_i13[bask + j]; real val = p3scale; for (int jj = 0; jj < polgrp::np11[i]; ++jj) { - if (k == polgrp::ip11[i * maxp11 + jj]) val = p3iscale; + if (k == polgrp::ip11[i * maxp11 + jj]) + val = p3iscale; } k -= 1; - if (k > i) { insert_mdpu(ik_scale, i, k, val, 'p'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, val, 'p'); + } } } @@ -338,10 +350,13 @@ static void mdpuscaleData(RcOp op) int k = couple_i14[bask + j]; real val = p4scale; for (int jj = 0; jj < polgrp::np11[i]; ++jj) { - if (k == polgrp::ip11[i * maxp11 + jj]) val = p4iscale; + if (k == polgrp::ip11[i * maxp11 + jj]) + val = p4iscale; } k -= 1; - if (k > i) { insert_mdpu(ik_scale, i, k, val, 'p'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, val, 'p'); + } } } @@ -352,10 +367,13 @@ static void mdpuscaleData(RcOp op) int k = couple_i15[bask + j]; real val = p5scale; for (int jj = 0; jj < polgrp::np11[i]; ++jj) { - if (k == polgrp::ip11[i * maxp11 + jj]) val = p5iscale; + if (k == polgrp::ip11[i * maxp11 + jj]) + val = p5iscale; } k -= 1; - if (k > i) { insert_mdpu(ik_scale, i, k, val, 'p'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, val, 'p'); + } } } @@ -365,7 +383,9 @@ static void mdpuscaleData(RcOp op) bask = i * maxp11; for (int j = 0; j < nn; ++j) { int k = polgrp::ip11[bask + j] - 1; - if (k > i) { insert_mdpu(ik_scale, i, k, d1scale, 'd'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, d1scale, 'd'); + } } } @@ -374,7 +394,9 @@ static void mdpuscaleData(RcOp op) bask = i * maxp12; for (int j = 0; j < nn; ++j) { int k = polgrp::ip12[bask + j] - 1; - if (k > i) { insert_mdpu(ik_scale, i, k, d2scale, 'd'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, d2scale, 'd'); + } } } @@ -383,7 +405,9 @@ static void mdpuscaleData(RcOp op) bask = i * maxp13; for (int j = 0; j < nn; ++j) { int k = polgrp::ip13[bask + j] - 1; - if (k > i) { insert_mdpu(ik_scale, i, k, d3scale, 'd'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, d3scale, 'd'); + } } } @@ -392,7 +416,9 @@ static void mdpuscaleData(RcOp op) bask = i * maxp14; for (int j = 0; j < nn; ++j) { int k = polgrp::ip14[bask + j] - 1; - if (k > i) { insert_mdpu(ik_scale, i, k, d4scale, 'd'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, d4scale, 'd'); + } } } @@ -402,7 +428,9 @@ static void mdpuscaleData(RcOp op) bask = i * maxp11; for (int j = 0; j < nn; ++j) { int k = polgrp::ip11[bask + j] - 1; - if (k > i) { insert_mdpu(ik_scale, i, k, u1scale, 'u'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, u1scale, 'u'); + } } } @@ -411,7 +439,9 @@ static void mdpuscaleData(RcOp op) bask = i * maxp12; for (int j = 0; j < nn; ++j) { int k = polgrp::ip12[bask + j] - 1; - if (k > i) { insert_mdpu(ik_scale, i, k, u2scale, 'u'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, u2scale, 'u'); + } } } @@ -420,7 +450,9 @@ static void mdpuscaleData(RcOp op) bask = i * maxp13; for (int j = 0; j < nn; ++j) { int k = polgrp::ip13[bask + j] - 1; - if (k > i) { insert_mdpu(ik_scale, i, k, u3scale, 'u'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, u3scale, 'u'); + } } } @@ -429,7 +461,9 @@ static void mdpuscaleData(RcOp op) bask = i * maxp14; for (int j = 0; j < nn; ++j) { int k = polgrp::ip14[bask + j] - 1; - if (k > i) { insert_mdpu(ik_scale, i, k, u4scale, 'u'); } + if (k > i) { + insert_mdpu(ik_scale, i, k, u4scale, 'u'); + } } } } @@ -497,8 +531,7 @@ static void chgpenData(RcOp op) }; // mdw excl list - auto insert_mdw = [](std::map, mdw>& a, int i, int k, - real val, char ch) { + auto insert_mdw = [](std::map, mdw>& a, int i, int k, real val, char ch) { std::pair key; key.first = i; key.second = k; @@ -540,7 +573,8 @@ static void chgpenData(RcOp op) nn = couple::n12[i]; for (int j = 0; j < nn; ++j) { int k = couple::i12[i][j] - 1; - if (k > i) insert_mdw(ik_mdw, i, k, m2scale, 'm'); + if (k > i) + insert_mdw(ik_mdw, i, k, m2scale, 'm'); } } @@ -549,7 +583,8 @@ static void chgpenData(RcOp op) bask = i * maxn13; for (int j = 0; j < nn; ++j) { int k = couple::i13[bask + j] - 1; - if (k > i) insert_mdw(ik_mdw, i, k, m3scale, 'm'); + if (k > i) + insert_mdw(ik_mdw, i, k, m3scale, 'm'); } } @@ -558,7 +593,8 @@ static void chgpenData(RcOp op) bask = i * maxn14; for (int j = 0; j < nn; ++j) { int k = couple::i14[bask + j] - 1; - if (k > i) insert_mdw(ik_mdw, i, k, m4scale, 'm'); + if (k > i) + insert_mdw(ik_mdw, i, k, m4scale, 'm'); } } @@ -567,7 +603,8 @@ static void chgpenData(RcOp op) bask = i * maxn15; for (int j = 0; j < nn; ++j) { int k = couple::i15[bask + j] - 1; - if (k > i) insert_mdw(ik_mdw, i, k, m5scale, 'm'); + if (k > i) + insert_mdw(ik_mdw, i, k, m5scale, 'm'); } } } @@ -591,10 +628,13 @@ static void chgpenData(RcOp op) int k = couple_i12[bask + j]; real val = p2scale; for (int jj = 0; jj < polgrp::np11[i]; ++jj) { - if (k == polgrp::ip11[i * maxp11 + jj]) val = p2iscale; + if (k == polgrp::ip11[i * maxp11 + jj]) + val = p2iscale; } k -= 1; - if (k > i) { insert_mdw(ik_mdw, i, k, val, 'd'); } + if (k > i) { + insert_mdw(ik_mdw, i, k, val, 'd'); + } } } @@ -605,10 +645,13 @@ static void chgpenData(RcOp op) int k = couple_i13[bask + j]; real val = p3scale; for (int jj = 0; jj < polgrp::np11[i]; ++jj) { - if (k == polgrp::ip11[i * maxp11 + jj]) val = p3iscale; + if (k == polgrp::ip11[i * maxp11 + jj]) + val = p3iscale; } k -= 1; - if (k > i) { insert_mdw(ik_mdw, i, k, val, 'd'); } + if (k > i) { + insert_mdw(ik_mdw, i, k, val, 'd'); + } } } @@ -619,10 +662,13 @@ static void chgpenData(RcOp op) int k = couple_i14[bask + j]; real val = p4scale; for (int jj = 0; jj < polgrp::np11[i]; ++jj) { - if (k == polgrp::ip11[i * maxp11 + jj]) val = p4iscale; + if (k == polgrp::ip11[i * maxp11 + jj]) + val = p4iscale; } k -= 1; - if (k > i) { insert_mdw(ik_mdw, i, k, val, 'd'); } + if (k > i) { + insert_mdw(ik_mdw, i, k, val, 'd'); + } } } @@ -633,10 +679,13 @@ static void chgpenData(RcOp op) int k = couple_i15[bask + j]; real val = p5scale; for (int jj = 0; jj < polgrp::np11[i]; ++jj) { - if (k == polgrp::ip11[i * maxp11 + jj]) val = p5iscale; + if (k == polgrp::ip11[i * maxp11 + jj]) + val = p5iscale; } k -= 1; - if (k > i) { insert_mdw(ik_mdw, i, k, val, 'd'); } + if (k > i) { + insert_mdw(ik_mdw, i, k, val, 'd'); + } } } } @@ -765,4 +814,17 @@ void elecData(RcOp op) RcMan mdpuscale42{mdpuscaleData, op}; RcMan chgpen42{chgpenData, op}; } + +TINKER_FVOID2(acc1, cu1, exfieldCharge, int); +TINKER_FVOID2(acc1, cu1, exfieldDipole, int); +void exfield(int vers, int useDipole) +{ + if (not extfld::use_exfld) + return; + + if (useDipole) + TINKER_FCALL2(acc1, cu1, exfieldDipole, vers); + else + TINKER_FCALL2(acc1, cu1, exfieldCharge, vers); +} } diff --git a/test/cmakesrc.txt b/test/cmakesrc.txt index d3351bf63..cc18a782d 100644 --- a/test/cmakesrc.txt +++ b/test/cmakesrc.txt @@ -14,6 +14,7 @@ disp.cpp emhippo.cpp ephippo.cpp expol.cpp +extfield.cpp geom.cpp improp.cpp imptor.cpp diff --git a/test/extfield.cpp b/test/extfield.cpp new file mode 100644 index 000000000..43ff548f5 --- /dev/null +++ b/test/extfield.cpp @@ -0,0 +1,121 @@ +#include "ff/modamoeba.h" + +#include "test.h" +#include "testrt.h" + +using namespace tinker; + +static const std::string externalField = "external-field 150 -300 450\n"; +static const std::string mpoleOnly = "multipoleterm only\n"; +static const std::string polarOnly = "polarizeterm only\n"; + +TEST_CASE("External-Fields-MPole-Analyze", "[ff][extfield]") +{ + TestFile fx1(TINKER9_DIRSTR "/test/file/extfield/water4.xyz"); + TestFile fk1(TINKER9_DIRSTR "/test/file/extfield/water4.key", "", externalField + mpoleOnly); + TestFile fp1(TINKER9_DIRSTR "/test/file/commit_6fe8e913/water03.prm"); + const char* xn = "water4.xyz"; + const char* argv[] = {"dummy", xn}; + int argc = 2; + + const double eps_e = testGetEps(0.0001, 0.0001); + const double eps_g = testGetEps(0.0001, 0.0001); + const double eps_v = testGetEps(0.001, 0.001); + + TestReference r(TINKER9_DIRSTR "/test/ref/extfield.1.txt"); + auto ref_c = r.getCount(); + auto ref_e = r.getEnergy(); + auto ref_v = r.getVirial(); + auto ref_g = r.getGradient(); + + rc_flag = calc::xyz | calc::energy | calc::grad | calc::virial | calc::analyz; + testBeginWithArgs(argc, argv); + initialize(); + + energy(calc::v0); + COMPARE_REALS(esum, ref_e, eps_e); + + energy(calc::v1); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + energy(calc::v3); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_INTS(countReduce(nem), ref_c); + + energy(calc::v4); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v5); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v6); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + finish(); + testEnd(); +} + +TEST_CASE("External-Fields-2-Analyze", "[ff][extfield]") {} + +TEST_CASE("External-Fields-3-Analyze", "[ff][extfield]") {} + +TEST_CASE("External-Fields-MPole", "[ff][extfield]") +{ + TestFile fx1(TINKER9_DIRSTR "/test/file/extfield/water4.xyz"); + TestFile fk1(TINKER9_DIRSTR "/test/file/extfield/water4.key", "", externalField + mpoleOnly); + TestFile fp1(TINKER9_DIRSTR "/test/file/commit_6fe8e913/water03.prm"); + const char* xn = "water4.xyz"; + const char* argv[] = {"dummy", xn}; + int argc = 2; + + const double eps_e = testGetEps(0.0001, 0.0001); + const double eps_g = testGetEps(0.0001, 0.0001); + const double eps_v = testGetEps(0.001, 0.001); + + TestReference r(TINKER9_DIRSTR "/test/ref/extfield.1.txt"); + auto ref_e = r.getEnergy(); + auto ref_v = r.getVirial(); + auto ref_g = r.getGradient(); + + rc_flag = calc::xyz | calc::energy | calc::grad | calc::virial; + testBeginWithArgs(argc, argv); + initialize(); + + energy(calc::v0); + COMPARE_REALS(esum, ref_e, eps_e); + + energy(calc::v1); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + energy(calc::v4); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v5); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v6); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + finish(); + testEnd(); +} + +TEST_CASE("External-Fields-2", "[ff][extfield]") {} + +TEST_CASE("External-Fields-3", "[ff][extfield]") {} diff --git a/test/file/extfield/water4.key b/test/file/extfield/water4.key new file mode 100644 index 000000000..534ac00ec --- /dev/null +++ b/test/file/extfield/water4.key @@ -0,0 +1,14 @@ +parameters water03.prm +a-axis 18.643 +neighbor-list +list-buffer 0.2 +vdw-cutoff 9.0 +ewald +ewald-cutoff 7.0 +polarization MUTUAL +polar-eps 0.00001 + +bondterm none +angleterm none +ureybradterm none +vdwterm none diff --git a/test/file/extfield/water4.xyz b/test/file/extfield/water4.xyz new file mode 100644 index 000000000..d68696253 --- /dev/null +++ b/test/file/extfield/water4.xyz @@ -0,0 +1,13 @@ + 12 Water Cubic Box (18.643 Ang, 4 AMOEBA) + 1 O 8.679662 7.087692 -0.696862 1 2 3 + 2 H 7.809455 6.755792 -0.382259 2 1 + 3 H 8.722232 6.814243 -1.617561 2 1 + 4 O -0.117313 8.244447 6.837616 1 5 6 + 5 H 0.216892 7.895445 6.050027 2 4 + 6 H 0.444268 7.826013 7.530196 2 4 + 7 O 8.379057 -0.092611 6.814631 1 8 9 + 8 H 9.340423 0.098069 6.734062 2 7 + 9 H 7.939619 0.573676 6.269838 2 7 + 10 O 6.589952 1.844323 -6.923167 1 11 12 + 11 H 5.885429 2.402305 -6.717934 2 10 + 12 H 6.181533 1.062747 -7.273678 2 10 diff --git a/test/ref/extfield.1.txt b/test/ref/extfield.1.txt new file mode 100644 index 000000000..0354ac169 --- /dev/null +++ b/test/ref/extfield.1.txt @@ -0,0 +1,25 @@ + extfield.cpp + + Energy Component Breakdown : Kcal/mole Interactions + Atomic Multipoles 39.7407 33 + + Internal Virial Tensor : 8.034 -1.807 32.838 + -1.807 -20.769 -9.504 + 32.838 -9.504 43.938 + + Cartesian Gradient Breakdown over Individual Atoms : + + Type Atom dE/dX dE/dY dE/dZ Norm + + Anlyt 1 23.8380 -46.8453 61.6292 80.9993 + Anlyt 2 -12.7966 22.2736 -32.6743 41.5629 + Anlyt 3 -11.0559 24.5926 -28.9025 39.5270 + Anlyt 4 18.2390 -42.3533 59.5494 75.3167 + Anlyt 5 -12.0590 23.7001 -30.8210 40.7068 + Anlyt 6 -6.1755 18.6394 -28.7569 34.8213 + Anlyt 7 19.5687 -37.7986 62.2519 75.4120 + Anlyt 8 -10.3750 21.5653 -33.2617 40.9761 + Anlyt 9 -9.1377 16.2887 -29.0281 34.5174 + Anlyt 10 19.7605 -44.4330 68.7687 84.2253 + Anlyt 11 -8.6191 21.7018 -34.8026 41.9104 + Anlyt 12 -11.2502 22.6736 -33.9490 42.3461 diff --git a/test/ref/extfield.2.txt b/test/ref/extfield.2.txt new file mode 100644 index 000000000..e69de29bb diff --git a/test/ref/extfield.3.txt b/test/ref/extfield.3.txt new file mode 100644 index 000000000..e69de29bb From a2b2b331a4e933062b6414b2220075f94e45bd1d Mon Sep 17 00:00:00 2001 From: Zhi Wang Date: Wed, 8 Mar 2023 01:48:06 -0800 Subject: [PATCH 2/4] add cpu code for extfield+epolar --- include/ff/amoeba/epolar.h | 1 + include/ff/elec.h | 3 + src/acc/amoeba/empolenonewald.cpp | 26 ++++ src/acc/amoeba/epolarnonewald.cpp | 18 +++ src/amoeba/epolar.cpp | 9 ++ src/amoeba/field.cpp | 1 + src/elec.cpp | 9 ++ src/hippo/epolar.cpp | 1 + src/hippo/field.cpp | 2 + test/extfield.cpp | 204 +++++++++++++++++++++++++++++- test/ref/extfield.1.txt | 2 +- test/ref/extfield.2.txt | 25 ++++ test/ref/extfield.3.txt | 24 ++++ 13 files changed, 320 insertions(+), 5 deletions(-) diff --git a/include/ff/amoeba/epolar.h b/include/ff/amoeba/epolar.h index 7ed8e5d59..44b16b0f0 100644 --- a/include/ff/amoeba/epolar.h +++ b/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]); /// \} } diff --git a/include/ff/elec.h b/include/ff/elec.h index 61b887021..4c402d59e 100644 --- a/include/ff/elec.h +++ b/include/ff/elec.h @@ -69,6 +69,9 @@ 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. +); //====================================================================// // // diff --git a/src/acc/amoeba/empolenonewald.cpp b/src/acc/amoeba/empolenonewald.cpp index 00fc96447..6c1c5305f 100644 --- a/src/acc/amoeba/empolenonewald.cpp +++ b/src/acc/amoeba/empolenonewald.cpp @@ -261,4 +261,30 @@ void exfieldDipole_acc(int vers) } } } + +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; + } + } +} } diff --git a/src/acc/amoeba/epolarnonewald.cpp b/src/acc/amoeba/epolarnonewald.cpp index d73202568..07320b701 100644 --- a/src/acc/amoeba/epolarnonewald.cpp +++ b/src/acc/amoeba/epolarnonewald.cpp @@ -5,6 +5,7 @@ #include "ff/switch.h" #include "seq/pair_polar.h" #include "tool/gpucard.h" +#include #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 { diff --git a/src/amoeba/epolar.cpp b/src/amoeba/epolar.cpp index de08ed6dd..ebe817ec9 100644 --- a/src/amoeba/epolar.cpp +++ b/src/amoeba/epolar.cpp @@ -13,6 +13,7 @@ #include "tool/iofortstr.h" #include "tool/ioprint.h" #include +#include #include #include #include @@ -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); + } +} } diff --git a/src/amoeba/field.cpp b/src/amoeba/field.cpp index 2bcd7940c..8382bf971 100644 --- a/src/amoeba/field.cpp +++ b/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); } } diff --git a/src/elec.cpp b/src/elec.cpp index efbf48b65..a5bf56094 100644 --- a/src/elec.cpp +++ b/src/elec.cpp @@ -827,4 +827,13 @@ void exfield(int vers, int useDipole) else TINKER_FCALL2(acc1, cu1, exfieldCharge, vers); } + +TINKER_FVOID2(acc1, cu1, extfieldModifyDField, real (*)[3], real (*)[3]); +void extfieldModifyDField(real (*field)[3], real (*fieldp)[3]) +{ + if (not extfld::use_exfld) + return; + + TINKER_FCALL2(acc1, cu1, extfieldModifyDField, field, fieldp); +} } diff --git a/src/hippo/epolar.cpp b/src/hippo/epolar.cpp index d3e0eadb2..7b6a0ba5f 100644 --- a/src/hippo/epolar.cpp +++ b/src/hippo/epolar.cpp @@ -228,6 +228,7 @@ void epolarChgpen(int vers) epolarChgpenEwald(vers, use_cfgrad); else epolarChgpenNonEwald(vers, use_cfgrad); + epolarPairwiseExtfield(vers, uind); torque(vers, depx, depy, depz); if (use_cfgrad) dcflux(vers, depx, depy, depz, vir_ep); diff --git a/src/hippo/field.cpp b/src/hippo/field.cpp index 44c84371c..f132ab4f2 100644 --- a/src/hippo/field.cpp +++ b/src/hippo/field.cpp @@ -31,6 +31,7 @@ void dfieldChgpen(real (*field)[3]) dfieldChgpenEwald(field); else dfieldChgpenNonEwald(field); + extfieldModifyDField(field, nullptr); } } @@ -108,6 +109,7 @@ void dfieldAplus(real (*field)[3]) dfieldAplusEwald(field); else dfieldAplusNonEwald(field); + extfieldModifyDField(field, nullptr); } } diff --git a/test/extfield.cpp b/test/extfield.cpp index 43ff548f5..721f1a6a9 100644 --- a/test/extfield.cpp +++ b/test/extfield.cpp @@ -63,9 +63,111 @@ TEST_CASE("External-Fields-MPole-Analyze", "[ff][extfield]") testEnd(); } -TEST_CASE("External-Fields-2-Analyze", "[ff][extfield]") {} +TEST_CASE("External-Fields-Polarize-Analyze", "[ff][extfield]") +{ + TestFile fx1(TINKER9_DIRSTR "/test/file/extfield/water4.xyz"); + TestFile fk1(TINKER9_DIRSTR "/test/file/extfield/water4.key", "", externalField + polarOnly); + TestFile fp1(TINKER9_DIRSTR "/test/file/commit_6fe8e913/water03.prm"); + const char* xn = "water4.xyz"; + const char* argv[] = {"dummy", xn}; + int argc = 2; + + const double eps_e = testGetEps(0.0001, 0.0001); + const double eps_g = testGetEps(0.0001, 0.0001); + const double eps_v = testGetEps(0.001, 0.001); + + TestReference r(TINKER9_DIRSTR "/test/ref/extfield.2.txt"); + auto ref_c = r.getCount(); + auto ref_e = r.getEnergy(); + auto ref_v = r.getVirial(); + auto ref_g = r.getGradient(); + + rc_flag = calc::xyz | calc::energy | calc::grad | calc::virial | calc::analyz; + testBeginWithArgs(argc, argv); + initialize(); + + energy(calc::v0); + COMPARE_REALS(esum, ref_e, eps_e); + + energy(calc::v1); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + energy(calc::v3); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_INTS(countReduce(nep), ref_c); + + energy(calc::v4); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v5); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v6); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + finish(); + testEnd(); +} + +TEST_CASE("External-Fields-MPolar-Analyze", "[ff][extfield]") +{ + TestFile fx1(TINKER9_DIRSTR "/test/file/extfield/water4.xyz"); + TestFile fk1(TINKER9_DIRSTR "/test/file/extfield/water4.key", "", externalField); + TestFile fp1(TINKER9_DIRSTR "/test/file/commit_6fe8e913/water03.prm"); + const char* xn = "water4.xyz"; + const char* argv[] = {"dummy", xn}; + int argc = 2; + + const double eps_e = testGetEps(0.0001, 0.0001); + const double eps_g = testGetEps(0.0001, 0.0001); + const double eps_v = testGetEps(0.001, 0.001); + + TestReference r(TINKER9_DIRSTR "/test/ref/extfield.3.txt"); + auto ref_e = r.getEnergy(); + auto ref_v = r.getVirial(); + auto ref_g = r.getGradient(); + + rc_flag = calc::xyz | calc::energy | calc::grad | calc::virial | calc::analyz; + testBeginWithArgs(argc, argv); + initialize(); -TEST_CASE("External-Fields-3-Analyze", "[ff][extfield]") {} + energy(calc::v0); + COMPARE_REALS(esum, ref_e, eps_e); + + energy(calc::v1); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + energy(calc::v3); + COMPARE_REALS(esum, ref_e, eps_e); + + energy(calc::v4); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v5); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v6); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + finish(); + testEnd(); +} TEST_CASE("External-Fields-MPole", "[ff][extfield]") { @@ -116,6 +218,100 @@ TEST_CASE("External-Fields-MPole", "[ff][extfield]") testEnd(); } -TEST_CASE("External-Fields-2", "[ff][extfield]") {} +TEST_CASE("External-Fields-Polarize", "[ff][extfield]") +{ + TestFile fx1(TINKER9_DIRSTR "/test/file/extfield/water4.xyz"); + TestFile fk1(TINKER9_DIRSTR "/test/file/extfield/water4.key", "", externalField + polarOnly); + TestFile fp1(TINKER9_DIRSTR "/test/file/commit_6fe8e913/water03.prm"); + const char* xn = "water4.xyz"; + const char* argv[] = {"dummy", xn}; + int argc = 2; + + const double eps_e = testGetEps(0.0001, 0.0001); + const double eps_g = testGetEps(0.0001, 0.0001); + const double eps_v = testGetEps(0.001, 0.001); + + TestReference r(TINKER9_DIRSTR "/test/ref/extfield.2.txt"); + auto ref_e = r.getEnergy(); + auto ref_v = r.getVirial(); + auto ref_g = r.getGradient(); + + rc_flag = calc::xyz | calc::energy | calc::grad | calc::virial; + testBeginWithArgs(argc, argv); + initialize(); + + energy(calc::v0); + COMPARE_REALS(esum, ref_e, eps_e); -TEST_CASE("External-Fields-3", "[ff][extfield]") {} + energy(calc::v1); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + energy(calc::v4); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v5); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v6); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + finish(); + testEnd(); +} + +TEST_CASE("External-Fields-MPolar", "[ff][extfield]") +{ + TestFile fx1(TINKER9_DIRSTR "/test/file/extfield/water4.xyz"); + TestFile fk1(TINKER9_DIRSTR "/test/file/extfield/water4.key", "", externalField); + TestFile fp1(TINKER9_DIRSTR "/test/file/commit_6fe8e913/water03.prm"); + const char* xn = "water4.xyz"; + const char* argv[] = {"dummy", xn}; + int argc = 2; + + const double eps_e = testGetEps(0.0001, 0.0001); + const double eps_g = testGetEps(0.0001, 0.0001); + const double eps_v = testGetEps(0.001, 0.001); + + TestReference r(TINKER9_DIRSTR "/test/ref/extfield.3.txt"); + auto ref_e = r.getEnergy(); + auto ref_v = r.getVirial(); + auto ref_g = r.getGradient(); + + rc_flag = calc::xyz | calc::energy | calc::grad | calc::virial; + testBeginWithArgs(argc, argv); + initialize(); + + energy(calc::v0); + COMPARE_REALS(esum, ref_e, eps_e); + + energy(calc::v1); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + energy(calc::v4); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v5); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v6); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + finish(); + testEnd(); +} diff --git a/test/ref/extfield.1.txt b/test/ref/extfield.1.txt index 0354ac169..7943cde2f 100644 --- a/test/ref/extfield.1.txt +++ b/test/ref/extfield.1.txt @@ -1,4 +1,4 @@ - extfield.cpp + extfield.cpp -- multipole term Energy Component Breakdown : Kcal/mole Interactions Atomic Multipoles 39.7407 33 diff --git a/test/ref/extfield.2.txt b/test/ref/extfield.2.txt index e69de29bb..209b60b44 100644 --- a/test/ref/extfield.2.txt +++ b/test/ref/extfield.2.txt @@ -0,0 +1,25 @@ + extfield.cpp -- polarizeterm term + + Energy Component Breakdown : Kcal/mole Interactions + Polarization -142.5085 21 + + Internal Virial Tensor : -15.107 7.941 -9.516 + 7.941 -18.664 14.278 + -9.516 14.278 -35.502 + + Cartesian Gradient Breakdown over Individual Atoms : + + Type Atom dE/dX dE/dY dE/dZ Norm + + Anlyt 1 -4.0805 -1.1001 -5.8071 7.1821 + Anlyt 2 1.8315 4.4798 -6.1608 7.8344 + Anlyt 3 2.2631 -3.4294 11.8448 12.5372 + Anlyt 4 7.4558 -7.7717 5.6988 12.1846 + Anlyt 5 0.4048 -3.4424 8.1070 8.8168 + Anlyt 6 -7.9318 11.2394 -13.8113 19.4933 + Anlyt 7 -1.5917 9.5515 -9.4165 13.5068 + Anlyt 8 -3.9535 1.8939 -4.1680 6.0489 + Anlyt 9 5.5135 -11.5339 13.7019 18.7396 + Anlyt 10 -5.0904 1.2074 -1.7957 5.5312 + Anlyt 11 4.4017 -6.5569 2.0021 8.1471 + Anlyt 12 0.7963 5.3972 -0.2778 5.4627 diff --git a/test/ref/extfield.3.txt b/test/ref/extfield.3.txt index e69de29bb..23ffe37c8 100644 --- a/test/ref/extfield.3.txt +++ b/test/ref/extfield.3.txt @@ -0,0 +1,24 @@ + extfield.cpp -- multipole and polarization terms + + Total Potential Energy : -102.7678 Kcal/mole + + Internal Virial Tensor : -7.074 6.134 23.322 + 6.134 -39.433 4.774 + 23.322 4.774 8.435 + + Cartesian Gradient Breakdown over Individual Atoms : + + Type Atom dE/dX dE/dY dE/dZ Norm + + Anlyt 1 19.7575 -47.9454 55.8221 76.1920 + Anlyt 2 -10.9650 26.7534 -38.8351 48.4163 + Anlyt 3 -8.7928 21.1632 -17.0576 28.5685 + Anlyt 4 25.6949 -50.1250 65.2482 86.1979 + Anlyt 5 -11.6542 20.2577 -22.7140 32.5902 + Anlyt 6 -14.1073 29.8788 -42.5682 53.8870 + Anlyt 7 17.9770 -28.2471 52.8354 62.5512 + Anlyt 8 -14.3286 23.4592 -37.4297 46.4394 + Anlyt 9 -3.6242 4.7548 -15.3262 16.4510 + Anlyt 10 14.6701 -43.2256 66.9730 81.0497 + Anlyt 11 -4.2175 15.1449 -32.8006 36.3735 + Anlyt 12 -10.4539 28.0708 -34.2267 45.4832 From cabea238673bfb856199b2fbff05af044f22ac3a Mon Sep 17 00:00:00 2001 From: Zhi Wang Date: Wed, 8 Mar 2023 02:41:08 -0800 Subject: [PATCH 3/4] add cuda code for extfield --- src/amoeba/emplar.cpp | 5 +- src/cu/amoeba/empole.cu | 100 ++++++++++++++++++++++++++++++++++++++-- src/cu/amoeba/epolar.cu | 25 +++++++++- src/cu/echarge.cu | 50 ++++++++++++++++++++ 4 files changed, 172 insertions(+), 8 deletions(-) diff --git a/src/amoeba/emplar.cpp b/src/amoeba/emplar.cpp index 328a4ccba..e04594ef6 100644 --- a/src/amoeba/emplar.cpp +++ b/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; diff --git a/src/cu/amoeba/empole.cu b/src/cu/amoeba/empole.cu index eed092671..857892f8c 100644 --- a/src/cu/amoeba/empole.cu +++ b/src/cu/amoeba/empole.cu @@ -1,5 +1,5 @@ -#include "ff/modamoeba.h" #include "ff/image.h" +#include "ff/modamoeba.h" #include "ff/pme.h" #include "ff/spatial.h" #include "ff/switch.h" @@ -7,6 +7,7 @@ #include "seq/launch.h" #include "seq/pair_mpole.h" #include "seq/triangle.h" +#include namespace tinker { #include "empole_cu1.cc" @@ -36,10 +37,9 @@ static void empole_cu() } } int ngrid = gpuGridSize(BLOCK_DIM); - empole_cu1<<>>(st.n, TINKER_IMAGE_ARGS, nem, em, vir_em, - demx, demy, demz, off, st.si1.bit0, nmdpuexclude, mdpuexclude, mdpuexclude_scale, st.x, st.y, - st.z, st.sorted, st.nakpl, st.iakpl, st.niak, st.iak, st.lst, trqx, trqy, trqz, rpole, f, - aewald); + empole_cu1<<>>(st.n, TINKER_IMAGE_ARGS, nem, em, vir_em, demx, demy, demz, + off, st.si1.bit0, nmdpuexclude, mdpuexclude, mdpuexclude_scale, st.x, st.y, st.z, st.sorted, st.nakpl, st.iakpl, + st.niak, st.iak, st.lst, trqx, trqy, trqz, rpole, f, aewald); } void empoleNonEwald_cu(int vers) @@ -75,4 +75,94 @@ void empoleEwaldRealSelf_cu(int vers) empole_cu(); } } + +__global__ +static void exfieldDipole_cu1(CountBuffer restrict nem, EnergyBuffer restrict em, VirialBuffer vir_em, + grad_prec* restrict demx, grad_prec* restrict demy, grad_prec* restrict demz, real* restrict trqx, + real* restrict trqy, real* restrict trqz, int vers, int n, real f, real ef1, real ef2, real ef3, + const real (*restrict rpole)[10], const real* restrict x, const real* restrict y, const real* restrict z) +{ + bool do_e = vers & calc::energy; + bool do_a = vers & calc::analyz; + bool do_g = vers & calc::grad; + bool do_v = vers & calc::virial; + + int ithread = ITHREAD; + for (int ii = ithread; ii < n; ii += STRIDE) { + 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, ithread); + if (do_a) + atomic_add(1, nem, ithread); + } + 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; + real vxz = zi * frx + xi * frz; + real vyz = zi * fry + yi * frz; + atomic_add(vxx, vxy, vxz, vyy, vyz, vzz, vir_em, ithread); + } + } + } +} + +void exfieldDipole_cu(int vers) +{ + real f = electric / dielec; + real ef1 = extfld::exfld[0], ef2 = extfld::exfld[1], ef3 = extfld::exfld[2]; + launch_k1s(g::s0, n, exfieldDipole_cu1, nem, em, vir_em, demx, demy, demz, trqx, trqy, trqz, vers, n, f, ef1, ef2, + ef3, rpole, x, y, z); +} + +__global__ +static void extfieldModifyDField_cu1(real (*restrict field)[3], real (*restrict fieldp)[3], int n, real ex1, real ex2, + real ex3) +{ + int ithread = ITHREAD; + if (fieldp) { + for (int i = ithread; i < n; i += STRIDE) { + 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 { + for (int i = ithread; i < n; i += STRIDE) { + field[i][0] += ex1; + field[i][1] += ex2; + field[i][2] += ex3; + } + } +} + +void extfieldModifyDField_cu(real (*field)[3], real (*fieldp)[3]) +{ + real ex1 = extfld::exfld[0]; + real ex2 = extfld::exfld[1]; + real ex3 = extfld::exfld[2]; + launch_k1b(g::s0, n, extfieldModifyDField_cu1, field, fieldp, n, ex1, ex2, ex3); +} } diff --git a/src/cu/amoeba/epolar.cu b/src/cu/amoeba/epolar.cu index de9558c92..7c3a30590 100644 --- a/src/cu/amoeba/epolar.cu +++ b/src/cu/amoeba/epolar.cu @@ -1,6 +1,6 @@ -#include "ff/modamoeba.h" #include "ff/cumodamoeba.h" #include "ff/image.h" +#include "ff/modamoeba.h" #include "ff/pme.h" #include "ff/spatial.h" #include "ff/switch.h" @@ -8,10 +8,11 @@ #include "seq/launch.h" #include "seq/pair_polar.h" #include "seq/triangle.h" +#include namespace tinker { __global__ -void epolar0DotProd_cu1(int n, real f, EnergyBuffer restrict ep, const real (*restrict gpu_uind)[3], +static void epolar0DotProd_cu1(int n, real f, EnergyBuffer restrict ep, const real (*restrict gpu_uind)[3], const real (*restrict gpu_udirp)[3], const real* restrict polarity_inv) { int ithread = ITHREAD; @@ -27,6 +28,26 @@ void epolar0DotProd_cu(const real (*gpu_uind)[3], const real (*gpu_udirp)[3]) const real f = -0.5 * electric / dielec; launch_k1b(g::s0, n, epolar0DotProd_cu1, n, f, ep, gpu_uind, gpu_udirp, polarity_inv); } + +__global__ +static void epolarPairwiseExtfield_cu1(EnergyBuffer restrict ep, const real (*uind)[3], int n, real f, real ex1, + real ex2, real ex3) +{ + int ithread = ITHREAD; + for (int i = ithread; i < n; i += STRIDE) { + real e = uind[i][0] * ex1 + uind[i][1] * ex2 + uind[i][2] * ex3; + atomic_add(f * e, ep, ithread); + } +} + +void epolarPairwiseExtfield_cu(const real (*uind)[3]) +{ + const real f = -0.5 * electric / dielec; + real ex1 = extfld::exfld[0]; + real ex2 = extfld::exfld[1]; + real ex3 = extfld::exfld[2]; + launch_k1b(g::s0, n, epolarPairwiseExtfield_cu1, ep, uind, n, f, ex1, ex2, ex3); +} } namespace tinker { diff --git a/src/cu/echarge.cu b/src/cu/echarge.cu index f5f724a52..f2572ab06 100644 --- a/src/cu/echarge.cu +++ b/src/cu/echarge.cu @@ -8,6 +8,7 @@ #include "seq/launch.h" #include "seq/pair_charge.h" #include "seq/triangle.h" +#include namespace tinker { #include "echarge_cu1.cc" @@ -218,4 +219,53 @@ void echargeEwaldFphiSelf_cu(int vers) else if (vers == calc::v6) echargeFphiSelf_cu(); } + +__global__ +static void exfieldCharge_cu1(CountBuffer restrict nec, EnergyBuffer restrict ec, VirialBuffer restrict vir_ec, + grad_prec* restrict decx, grad_prec* restrict decy, grad_prec* restrict decz, int vers, int n, real f, real ef1, + real ef2, real ef3, const real* restrict pchg, const real* restrict x, const real* restrict y, + const real* restrict z) +{ + bool do_e = vers & calc::energy; + bool do_a = vers & calc::analyz; + bool do_g = vers & calc::grad; + bool do_v = vers & calc::virial; + + int ithread = ITHREAD; + for (int ii = ithread; ii < n; ii += STRIDE) { + 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, ithread); + if (do_a) + atomic_add(1, nec, ithread); + } + 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; + real vxz = zi * frx + xi * frz; + real vyz = zi * fry + yi * frz; + atomic_add(vxx, vxy, vxz, vyy, vyz, vzz, vir_ec, ithread); + } + } + } +} + +void exfieldCharge_cu(int vers) +{ + real f = electric / dielec; + real ef1 = extfld::exfld[0], ef2 = extfld::exfld[1], ef3 = extfld::exfld[2]; + launch_k1s(g::s0, n, exfieldCharge_cu1, nec, ec, vir_ec, decx, decy, decz, vers, n, f, ef1, ef2, ef3, pchg, x, y, z); +} } From 869c4e3802c8135b8b336a14721b962200327ce9 Mon Sep 17 00:00:00 2001 From: Zhi Wang Date: Wed, 8 Mar 2023 21:36:18 -0800 Subject: [PATCH 4/4] pchg extfield code; divide off diag virial by 2 --- src/acc/amoeba/empolenonewald.cpp | 6 +- src/acc/echarge.cpp | 6 +- src/cu/amoeba/empole.cu | 6 +- src/cu/echarge.cu | 6 +- src/echarge.cpp | 4 +- src/echglj.cpp | 1 + test/extfield.cpp | 101 ++++++++++++++++++++++++++++++ test/file/extfield/amber.key | 10 +++ test/file/extfield/amber.xyz | 13 ++++ test/ref/extfield.1.txt | 6 +- test/ref/extfield.3.txt | 6 +- test/ref/extfield.4.txt | 24 +++++++ 12 files changed, 170 insertions(+), 19 deletions(-) create mode 100644 test/file/extfield/amber.key create mode 100644 test/file/extfield/amber.xyz create mode 100644 test/ref/extfield.4.txt diff --git a/src/acc/amoeba/empolenonewald.cpp b/src/acc/amoeba/empolenonewald.cpp index 6c1c5305f..bdef159d4 100644 --- a/src/acc/amoeba/empolenonewald.cpp +++ b/src/acc/amoeba/empolenonewald.cpp @@ -253,9 +253,9 @@ void exfieldDipole_acc(int vers) real vxx = xi * frx; real vyy = yi * fry; real vzz = zi * frz; - real vxy = yi * frx + xi * fry; - real vxz = zi * frx + xi * frz; - real vyz = zi * fry + yi * 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); } } diff --git a/src/acc/echarge.cpp b/src/acc/echarge.cpp index 6d6678b35..7023ed7e0 100644 --- a/src/acc/echarge.cpp +++ b/src/acc/echarge.cpp @@ -404,9 +404,9 @@ void exfieldCharge_acc(int vers) real vxx = xi * frx; real vyy = yi * fry; real vzz = zi * frz; - real vxy = yi * frx + xi * fry; - real vxz = zi * frx + xi * frz; - real vyz = zi * fry + yi * 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); } } diff --git a/src/cu/amoeba/empole.cu b/src/cu/amoeba/empole.cu index 857892f8c..1b77ff5f6 100644 --- a/src/cu/amoeba/empole.cu +++ b/src/cu/amoeba/empole.cu @@ -118,9 +118,9 @@ static void exfieldDipole_cu1(CountBuffer restrict nem, EnergyBuffer restrict em real vxx = xi * frx; real vyy = yi * fry; real vzz = zi * frz; - real vxy = yi * frx + xi * fry; - real vxz = zi * frx + xi * frz; - real vyz = zi * fry + yi * 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, ithread); } } diff --git a/src/cu/echarge.cu b/src/cu/echarge.cu index f2572ab06..18104bd82 100644 --- a/src/cu/echarge.cu +++ b/src/cu/echarge.cu @@ -253,9 +253,9 @@ static void exfieldCharge_cu1(CountBuffer restrict nec, EnergyBuffer restrict ec real vxx = xi * frx; real vyy = yi * fry; real vzz = zi * frz; - real vxy = yi * frx + xi * fry; - real vxz = zi * frx + xi * frz; - real vyz = zi * fry + yi * 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, ithread); } } diff --git a/src/echarge.cpp b/src/echarge.cpp index 99b934602..5f9ae1d24 100644 --- a/src/echarge.cpp +++ b/src/echarge.cpp @@ -196,8 +196,10 @@ void echarge(int vers) echargeEwaldRecipSelf(vers); TINKER_FCALL2(acc1, cu1, echargeEwaldReal, vers); pmeStreamFinishWait(use_pme_stream and static_cast(vers & calc::analyz)); - } else + } else { echargeNonEwald(vers); + } + exfield(vers, 0); if (rc_a) { if (do_e) { diff --git a/src/echglj.cpp b/src/echglj.cpp index adab908e5..f337d35df 100644 --- a/src/echglj.cpp +++ b/src/echglj.cpp @@ -288,6 +288,7 @@ void echglj(int vers) } else { TINKER_FCALL2(acc0, cu1, echgljRadArithEpsGeomNonEwald, vers); } + exfield(vers, 0); if (do_e) { if (elrc_vol != 0) { diff --git a/test/extfield.cpp b/test/extfield.cpp index 721f1a6a9..ba41312d0 100644 --- a/test/extfield.cpp +++ b/test/extfield.cpp @@ -169,6 +169,58 @@ TEST_CASE("External-Fields-MPolar-Analyze", "[ff][extfield]") testEnd(); } +TEST_CASE("External-Fields-VdwPchg-Analyze", "[ff][extfield]") +{ + TestFile fx1(TINKER9_DIRSTR "/test/file/extfield/amber.xyz"); + TestFile fk1(TINKER9_DIRSTR "/test/file/extfield/amber.key"); + TestFile fp1(TINKER9_DIRSTR "/test/file/commit_350df099/amber99sb.prm"); + const char* xn = "amber.xyz"; + const char* argv[] = {"dummy", xn}; + int argc = 2; + + const double eps_e = testGetEps(0.0001, 0.0001); + const double eps_g = testGetEps(0.0001, 0.0001); + const double eps_v = testGetEps(0.001, 0.001); + + TestReference r(TINKER9_DIRSTR "/test/ref/extfield.4.txt"); + auto ref_e = r.getEnergy(); + auto ref_v = r.getVirial(); + auto ref_g = r.getGradient(); + + rc_flag = calc::xyz | calc::energy | calc::grad | calc::virial | calc::analyz; + testBeginWithArgs(argc, argv); + initialize(); + + energy(calc::v0); + COMPARE_REALS(esum, ref_e, eps_e); + + energy(calc::v1); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + energy(calc::v3); + COMPARE_REALS(esum, ref_e, eps_e); + + energy(calc::v4); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v5); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v6); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + finish(); + testEnd(); +} + TEST_CASE("External-Fields-MPole", "[ff][extfield]") { TestFile fx1(TINKER9_DIRSTR "/test/file/extfield/water4.xyz"); @@ -315,3 +367,52 @@ TEST_CASE("External-Fields-MPolar", "[ff][extfield]") finish(); testEnd(); } + +TEST_CASE("External-Fields-VdwPchg", "[ff][extfield]") +{ + TestFile fx1(TINKER9_DIRSTR "/test/file/extfield/amber.xyz"); + TestFile fk1(TINKER9_DIRSTR "/test/file/extfield/amber.key"); + TestFile fp1(TINKER9_DIRSTR "/test/file/commit_350df099/amber99sb.prm"); + const char* xn = "amber.xyz"; + const char* argv[] = {"dummy", xn}; + int argc = 2; + + const double eps_e = testGetEps(0.0001, 0.0001); + const double eps_g = testGetEps(0.0001, 0.0001); + const double eps_v = testGetEps(0.001, 0.001); + + TestReference r(TINKER9_DIRSTR "/test/ref/extfield.4.txt"); + auto ref_e = r.getEnergy(); + auto ref_v = r.getVirial(); + auto ref_g = r.getGradient(); + + rc_flag = calc::xyz | calc::energy | calc::grad | calc::virial; + testBeginWithArgs(argc, argv); + initialize(); + + energy(calc::v0); + COMPARE_REALS(esum, ref_e, eps_e); + + energy(calc::v1); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + energy(calc::v4); + COMPARE_REALS(esum, ref_e, eps_e); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v5); + COMPARE_GRADIENT(ref_g, eps_g); + + energy(calc::v6); + COMPARE_GRADIENT(ref_g, eps_g); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v); + + finish(); + testEnd(); +} diff --git a/test/file/extfield/amber.key b/test/file/extfield/amber.key new file mode 100644 index 000000000..a64a8511b --- /dev/null +++ b/test/file/extfield/amber.key @@ -0,0 +1,10 @@ +parameters amber99sb.prm +a-axis 18.643 +neighbor-list +list-buffer 0.2 +cutoff 9.0 + +bondterm none +angleterm none + +external-field 150 -300 450 diff --git a/test/file/extfield/amber.xyz b/test/file/extfield/amber.xyz new file mode 100644 index 000000000..3bd4ea3fc --- /dev/null +++ b/test/file/extfield/amber.xyz @@ -0,0 +1,13 @@ + 12 Water Cubic Box (18.643 Ang) + 1 O 8.679662 7.087692 -0.696862 2001 2 3 + 2 H 7.809455 6.755792 -0.382259 2002 1 + 3 H 8.722232 6.814243 -1.617561 2002 1 + 4 O -0.117313 8.244447 6.837616 2001 5 6 + 5 H 0.216892 7.895445 6.050027 2002 4 + 6 H 0.444268 7.826013 7.530196 2002 4 + 7 O 8.379057 -0.092611 6.814631 2001 8 9 + 8 H 9.340423 0.098069 6.734062 2002 7 + 9 H 7.939619 0.573676 6.269838 2002 7 + 10 O 6.589952 1.844323 -6.923167 2001 11 12 + 11 H 5.885429 2.402305 -6.717934 2002 10 + 12 H 6.181533 1.062747 -7.273678 2002 10 diff --git a/test/ref/extfield.1.txt b/test/ref/extfield.1.txt index 7943cde2f..368fb9408 100644 --- a/test/ref/extfield.1.txt +++ b/test/ref/extfield.1.txt @@ -3,9 +3,9 @@ Energy Component Breakdown : Kcal/mole Interactions Atomic Multipoles 39.7407 33 - Internal Virial Tensor : 8.034 -1.807 32.838 - -1.807 -20.769 -9.504 - 32.838 -9.504 43.938 + Internal Virial Tensor : 8.034 -0.430 19.176 + -0.430 -20.769 -6.245 + 19.176 -6.245 43.938 Cartesian Gradient Breakdown over Individual Atoms : diff --git a/test/ref/extfield.3.txt b/test/ref/extfield.3.txt index 23ffe37c8..2780c6e6b 100644 --- a/test/ref/extfield.3.txt +++ b/test/ref/extfield.3.txt @@ -2,9 +2,9 @@ Total Potential Energy : -102.7678 Kcal/mole - Internal Virial Tensor : -7.074 6.134 23.322 - 6.134 -39.433 4.774 - 23.322 4.774 8.435 + Internal Virial Tensor : -7.074 7.510 9.659 + 7.510 -39.433 8.033 + 9.659 8.033 8.435 Cartesian Gradient Breakdown over Individual Atoms : diff --git a/test/ref/extfield.4.txt b/test/ref/extfield.4.txt new file mode 100644 index 000000000..155bccac1 --- /dev/null +++ b/test/ref/extfield.4.txt @@ -0,0 +1,24 @@ + extfield.cpp -- van der waals and charge terms + + Total Potential Energy : 49.6917 Kcal/mole + + Internal Virial Tensor : 8.093 -2.613 21.038 + -2.613 -20.550 -4.143 + 21.038 -4.143 65.691 + + Cartesian Gradient Breakdown over Individual Atoms : + + Type Atom dE/dX dE/dY dE/dZ Norm + + Anlyt 1 28.8952 -57.6770 86.5850 107.9747 + Anlyt 2 -14.4549 28.8391 -43.3016 53.9969 + Anlyt 3 -14.4257 29.0381 -43.0993 53.9339 + Anlyt 4 28.8477 -57.6955 86.5432 107.9384 + Anlyt 5 -14.4239 28.8477 -43.2716 53.9692 + Anlyt 6 -14.4221 28.8460 -43.2702 53.9666 + Anlyt 7 29.4518 -57.2409 86.9148 108.1578 + Anlyt 8 -14.5032 28.5911 -43.6526 54.1604 + Anlyt 9 -14.7044 28.6645 -43.4550 54.0945 + Anlyt 10 28.6892 -59.1323 85.3912 107.7560 + Anlyt 11 -14.3542 29.5544 -42.6159 53.8110 + Anlyt 12 -14.5956 29.3647 -42.7681 53.8927