Add virtual site with one constructing atom
authorBerk Hess <hess@kth.se>
Tue, 23 Jun 2020 12:54:12 +0000 (12:54 +0000)
committerPaul Bauer <paul.bauer.q@gmail.com>
Tue, 23 Jun 2020 12:54:12 +0000 (12:54 +0000)
Also fixed the order of the function type update table in tpxio.cpp.

docs/reference-manual/functions/interaction-methods.rst
docs/release-notes/2021/major/features.rst
src/gromacs/fileio/tpxio.cpp
src/gromacs/gmxlib/nrnb.cpp
src/gromacs/gmxlib/nrnb.h
src/gromacs/gmxpreprocess/convparm.cpp
src/gromacs/mdlib/vsite.cpp
src/gromacs/mdlib/vsite.h
src/gromacs/topology/ifunc.cpp
src/gromacs/topology/ifunc.h

index 25820ca0d1395f21cf015efdf859c50191915122..4d88abfb0722eaaa8827f6e5ea96c9f5a8f01104 100644 (file)
@@ -197,6 +197,20 @@ below. Constructing atoms in virtual sites can be virtual sites
 themselves, but only if they are higher in the list, i.e. virtual sites
 can be constructed from “particles” that are simpler virtual sites.
 
+-  On top of an atom. This allows giving an atom multiple atom types and
+   with that also assigned multiple, different bonded interactions. This
+   can espically be of use in free-energy calculations.
+
+-  The coordinates of the virtual site equal that of the constructing atom:
+
+   .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i
+             :label: eqnvsite1
+
+-  The force is moved to the constructing atom:
+
+   .. math:: \mathbf{F}_i ~=~ \mathbf{F}_{s}
+             :label: eqnvsite1force
+
 -  As a linear combination of two atoms
    (:numref:`Fig. %s <fig-vsites>` 2):
 
index 82a35ba5d15d17f57408bf2c4b1c61809d7d5bd4..7a5626f7cac69e84eb7fccd55b52ce93013b9a84 100644 (file)
@@ -7,6 +7,12 @@ New and improved features
    Also, please use the syntax :issue:`number` to reference issues on GitLab, without the
    a space between the colon and number!
 
+Virtual site with single constructing atom
+""""""""""""""""""""""""""""""""""""""""""
+
+Added a virtual site that is constructed on top if its single constructing
+atom. This can be useful for free-energy calculations.
+
 Lower energy drift due to SETTLE
 """"""""""""""""""""""""""""""""
 
index 68b4f3e95c6a4a9797f30bf71a30113aff4b52cd..e20769039028d42235c0c149972c50252aef6ff9 100644 (file)
@@ -131,6 +131,7 @@ enum tpxv
     tpxv_VSite2FD,                  /**< Added 2FD type virtual site */
     tpxv_AddSizeField, /**< Added field with information about the size of the serialized tpr file in bytes, excluding the header */
     tpxv_StoreNonBondedInteractionExclusionGroup, /**< Store the non bonded interaction exclusion group in the topology */
+    tpxv_VSite1,                                  /**< Added 1 type virtual site */
     tpxv_Count                                    /**< the total number of tpxv versions */
 };
 
@@ -205,10 +206,13 @@ static const t_ftupd ftupd[] = {
     { 72, F_GBPOL_NOLONGERUSED },
     { 72, F_NPSOLVATION_NOLONGERUSED },
     { 93, F_LJ_RECIP },
+    { 76, F_ANHARM_POL },
     { 90, F_FBPOSRES },
+    { tpxv_VSite1, F_VSITE1 },
+    { tpxv_VSite2FD, F_VSITE2FD },
+    { tpxv_GenericInternalParameters, F_DENSITYFITTING },
     { 69, F_VTEMP_NOLONGERUSED },
     { 66, F_PDISPCORR },
-    { 76, F_ANHARM_POL },
     { 79, F_DVDL_COUL },
     {
             79,
@@ -220,8 +224,6 @@ static const t_ftupd ftupd[] = {
     },
     { 79, F_DVDL_RESTRAINT },
     { 79, F_DVDL_TEMPERATURE },
-    { tpxv_GenericInternalParameters, F_DENSITYFITTING },
-    { tpxv_VSite2FD, F_VSITE2FD },
 };
 #define NFTUPD asize(ftupd)
 
@@ -1940,6 +1942,7 @@ static void do_iparams(gmx::ISerializer* serializer, t_functype ftype, t_iparams
             serializer->doReal(&iparams->settle.doh);
             serializer->doReal(&iparams->settle.dhh);
             break;
+        case F_VSITE1: break; // VSite1 has 0 parameters
         case F_VSITE2:
         case F_VSITE2FD: serializer->doReal(&iparams->vsite.a); break;
         case F_VSITE3:
index 1fca45d7cddb063fb7e15a2a27810514ac6dcfa5..1366d351e51c93f6abbec57f2234fd6300bea1e5 100644 (file)
@@ -175,6 +175,7 @@ static const t_nrnb_data nbdata[eNRNB] = {
     { "Shake-Init", 10 },
     { "Constraint-Vir", 24 },
     { "Settle", 370 },
+    { "Virtual Site 1", 1 },
     { "Virtual Site 2", 23 },
     { "Virtual Site 2fd", 63 },
     { "Virtual Site 3", 37 },
index e5129fcdde983e0385c287ac007deaa398aae1a0..847c6cdb9d345ca3d168c95e7cbf42d3c7cdc0c4 100644 (file)
@@ -155,6 +155,7 @@ enum
     eNR_SHAKE_RIJ,
     eNR_CONSTR_VIR,
     eNR_SETTLE,
+    eNR_VSITE1,
     eNR_VSITE2,
     eNR_VSITE2FD,
     eNR_VSITE3,
index 9b4da4cce4f5f8aaf4bbdc7cd9715e56f76a9c2c..c9efb733cbe83608742057356d7cdd837eff6c5d 100644 (file)
@@ -402,6 +402,7 @@ static int assign_param(t_functype ftype, t_iparams* newparam, gmx::ArrayRef<con
             newparam->settle.doh = old[0];
             newparam->settle.dhh = old[1];
             break;
+        case F_VSITE1:
         case F_VSITE2:
         case F_VSITE2FD:
         case F_VSITE3:
index f2c8b0ef666388eff2c6ea76dc6c5e5be4e588ab..aab8dffd747d508f0d7af2505b03f5d1f5d158ae 100644 (file)
@@ -334,6 +334,13 @@ static inline real inverseNorm(const rvec x)
 #ifndef DOXYGEN
 /* Vsite construction routines */
 
+static void constr_vsite1(const rvec xi, rvec x)
+{
+    copy_rvec(xi, x);
+
+    /* TOTAL: 0 flops */
+}
+
 static void constr_vsite2(const rvec xi, const rvec xj, rvec x, real a, const t_pbc* pbc)
 {
     real b = 1 - a;
@@ -675,6 +682,7 @@ static void construct_vsites_thread(ArrayRef<RVec>                  x,
                 real b1, c1;
                 switch (ftype)
                 {
+                    case F_VSITE1: constr_vsite1(x[ai], x[avsite]); break;
                     case F_VSITE2:
                         aj = ia[3];
                         constr_vsite2(x[ai], x[aj], x[avsite], a1, pbc_null2);
@@ -857,6 +865,14 @@ void constructVirtualSites(ArrayRef<RVec> x, ArrayRef<const t_iparams> ip, Array
 #ifndef DOXYGEN
 /* Force spreading routines */
 
+static void spread_vsite1(const t_iatom ia[], ArrayRef<RVec> f)
+{
+    const int av = ia[1];
+    const int ai = ia[2];
+
+    f[av] += f[ai];
+}
+
 template<VirialHandling virialHandling>
 static void spread_vsite2(const t_iatom        ia[],
                           real                 a,
@@ -1761,6 +1777,7 @@ static void spreadForceForThread(ArrayRef<const RVec>            x,
                 /* Construct the vsite depending on type */
                 switch (ftype)
                 {
+                    case F_VSITE1: spread_vsite1(ia, f); break;
                     case F_VSITE2:
                         spread_vsite2<virialHandling>(ia, a1, x, f, fshift, pbc_null2);
                         break;
@@ -2034,6 +2051,7 @@ void VirtualSitesHandler::Impl::spreadForces(ArrayRef<const RVec> x,
         dd_move_f_vsites(*domainInfo_.domdec_, f, fshift);
     }
 
+    inc_nrnb(nrnb, eNR_VSITE1, vsite_count(ilists_, F_VSITE1));
     inc_nrnb(nrnb, eNR_VSITE2, vsite_count(ilists_, F_VSITE2));
     inc_nrnb(nrnb, eNR_VSITE2FD, vsite_count(ilists_, F_VSITE2FD));
     inc_nrnb(nrnb, eNR_VSITE3, vsite_count(ilists_, F_VSITE3));
index be1d8c9ebc0f46c6fd4ba3738ac8cc9d1151b9e1..ae7c932b5e5caf390dccc20acd6cc96d107eeccf 100644 (file)
@@ -74,7 +74,7 @@ class RangePartitioning;
  * This is used to avoid loops over all ftypes just to get the vsite entries.
  * (We should replace the fixed ilist array by only the used entries.)
  */
-static constexpr int c_ftypeVsiteStart = F_VSITE2;
+static constexpr int c_ftypeVsiteStart = F_VSITE1;
 //! The start and end value of the vsite indices in the ftype enum
 static constexpr int c_ftypeVsiteEnd = F_VSITEN + 1;
 
index 469d2bfbb1663eb137474c899969dd7125bdcb7f..e6b4028814a2a0069bb0ce823c60b66c0bd3ad55 100644 (file)
@@ -141,9 +141,10 @@ const t_interaction_function interaction_function[F_NRE] = {
     def_bonded("ANGRESZ", "Angle Rest. Z", 2, 3, 3), def_bonded("DIHRES", "Dih. Rest.", 4, 3, 3),
     def_nofc("DIHRESVIOL", "Dih. Rest. Viol."), /* obsolete */
     def_shkcb("CONSTR", "Constraint", 2, 1, 1), def_shk("CONSTRNC", "Constr. No Conn.", 2, 1, 1),
-    def_shkcb("SETTLE", "Settle", 3, 2, 0), def_vsite("VSITE2", "Virtual site 2", 3, 1),
-    def_vsite("VSITE2FD", "Virtual site 2fd", 3, 1), def_vsite("VSITE3", "Virtual site 3", 4, 2),
-    def_vsite("VSITE3FD", "Virtual site 3fd", 4, 2), def_vsite("VSITE3FAD", "Virtual site 3fad", 4, 2),
+    def_shkcb("SETTLE", "Settle", 3, 2, 0), def_vsite("VSITE1", "Virtual site 1", 2, 0),
+    def_vsite("VSITE2", "Virtual site 2", 3, 1), def_vsite("VSITE2FD", "Virtual site 2fd", 3, 1),
+    def_vsite("VSITE3", "Virtual site 3", 4, 2), def_vsite("VSITE3FD", "Virtual site 3fd", 4, 2),
+    def_vsite("VSITE3FAD", "Virtual site 3fad", 4, 2),
     def_vsite("VSITE3OUT", "Virtual site 3out", 4, 3), def_vsite("VSITE4FD", "Virtual site 4fd", 5, 3),
     def_vsite("VSITE4FDN", "Virtual site 4fdn", 5, 3), def_vsite("VSITEN", "Virtual site N", 2, 2),
     def_nofc("COM_PULL", "COM Pull En."), def_nofc("DENSITYFIT", "Density fitting"),
index 64cce554628be3c19c5d59249e6526112b94b6da..584c57f0b76fb6c34276478de5798b2b9363a74e 100644 (file)
@@ -188,6 +188,7 @@ enum
     F_CONSTR,
     F_CONSTRNC,
     F_SETTLE,
+    F_VSITE1,
     F_VSITE2,
     F_VSITE2FD,
     F_VSITE3,