Remove unused thole polarization rfac parameter
[alexxy/gromacs.git] / src / gromacs / topology / idef.h
index 3ebb53784cd4999fbd5dfa7e9d291fe98918cd50..60cf6bee0e122ce330e3cd199655bf9be1552800 100644 (file)
@@ -3,7 +3,8 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2018 by the GROMACS development team.
+ * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
 #ifndef GMX_TOPOLOGY_IDEF_H
 #define GMX_TOPOLOGY_IDEF_H
 
-#include <stdio.h>
+#include <cstdio>
 
 #include <array>
 #include <vector>
 
 #include "gromacs/math/vectypes.h"
 #include "gromacs/topology/ifunc.h"
-#include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/real.h"
 
+struct gmx_ffparams_t;
+
 typedef union t_iparams
 {
     /* Some parameters have A and B values for free energy calculations.
@@ -55,117 +57,159 @@ typedef union t_iparams
      * The harmonic type is used for all harmonic potentials:
      * bonds, angles and improper dihedrals
      */
-    struct {
+    struct
+    {
         real a, b, c;
     } bham;
-    struct {
+    struct
+    {
         real rA, krA, rB, krB;
     } harmonic;
-    struct {
+    struct
+    {
         real klinA, aA, klinB, aB;
     } linangle;
-    struct {
+    struct
+    {
         real lowA, up1A, up2A, kA, lowB, up1B, up2B, kB;
     } restraint;
     /* No free energy supported for cubic bonds, FENE, WPOL or cross terms */
-    struct {
+    struct
+    {
         real b0, kb, kcub;
     } cubic;
-    struct {
+    struct
+    {
         real bm, kb;
     } fene;
-    struct {
+    struct
+    {
         real r1e, r2e, krr;
     } cross_bb;
-    struct {
+    struct
+    {
         real r1e, r2e, r3e, krt;
     } cross_ba;
-    struct {
+    struct
+    {
         real thetaA, kthetaA, r13A, kUBA, thetaB, kthetaB, r13B, kUBB;
     } u_b;
-    struct {
+    struct
+    {
         real theta, c[5];
     } qangle;
-    struct {
+    struct
+    {
         real alpha;
     } polarize;
-    struct {
+    struct
+    {
         real alpha, drcut, khyp;
     } anharm_polarize;
-    struct {
+    struct
+    {
         real al_x, al_y, al_z, rOH, rHH, rOD;
     } wpol;
-    struct {
-        real a, alpha1, alpha2, rfac;
+    struct
+    {
+        real a, alpha1, alpha2;
     } thole;
-    struct {
+    struct
+    {
         real c6, c12;
     } lj;
-    struct {
+    struct
+    {
         real c6A, c12A, c6B, c12B;
     } lj14;
-    struct {
+    struct
+    {
         real fqq, qi, qj, c6, c12;
     } ljc14;
-    struct {
+    struct
+    {
         real qi, qj, c6, c12;
     } ljcnb;
     /* Proper dihedrals can not have different multiplicity when
      * doing free energy calculations, because the potential would not
      * be periodic anymore.
      */
-    struct {
-        real phiA, cpA; int mult; real phiB, cpB;
+    struct
+    {
+        real phiA, cpA;
+        int  mult;
+        real phiB, cpB;
     } pdihs;
-    struct {
+    struct
+    {
         real dA, dB;
     } constr;
     /* Settle can not be used for Free energy calculations of water bond geometry.
      * Use shake (or lincs) instead if you have to change the water bonds.
      */
-    struct {
+    struct
+    {
         real doh, dhh;
     } settle;
-    struct {
+    struct
+    {
         real b0A, cbA, betaA, b0B, cbB, betaB;
     } morse;
-    struct {
+    struct
+    {
         real pos0A[DIM], fcA[DIM], pos0B[DIM], fcB[DIM];
     } posres;
-    struct {
-        real pos0[DIM], r, k; int geom;
+    struct
+    {
+        real pos0[DIM], r, k;
+        int  geom;
     } fbposres;
-    struct {
+    struct
+    {
         real rbcA[NR_RBDIHS], rbcB[NR_RBDIHS];
     } rbdihs;
-    struct {
+    struct
+    {
         real cbtcA[NR_CBTDIHS], cbtcB[NR_CBTDIHS];
     } cbtdihs;
-    struct {
+    struct
+    {
         real a, b, c, d, e, f;
     } vsite;
-    struct {
-        int  n; real a;
+    struct
+    {
+        int  n;
+        real a;
     } vsiten;
     /* NOTE: npair is only set after reading the tpx file */
-    struct {
-        real low, up1, up2, kfac; int type, label, npair;
+    struct
+    {
+        real low, up1, up2, kfac;
+        int  type, label, npair;
     } disres;
-    struct {
+    struct
+    {
         real phiA, dphiA, kfacA, phiB, dphiB, kfacB;
     } dihres;
-    struct {
-        int  ex, power, label; real c, obs, kfac;
+    struct
+    {
+        int  ex, power, label;
+        real c, obs, kfac;
     } orires;
-    struct {
-        int  table; real kA; real kB;
+    struct
+    {
+        int  table;
+        real kA;
+        real kB;
     } tab;
-    struct {
+    struct
+    {
         int cmapA, cmapB;
     } cmap;
-    struct {
+    struct
+    {
         real buf[MAXFORCEPARAM];
-    } generic;                                               /* Conversion */
+    } generic; /* Conversion */
 } t_iparams;
 
 typedef int t_functype;
@@ -178,11 +222,45 @@ typedef int t_functype;
 struct InteractionList
 {
     /* Returns the total number of elements in iatoms */
-    int size() const
+    int size() const { return static_cast<int>(iatoms.size()); }
+
+    /* Returns whether the list is empty */
+    bool empty() const { return iatoms.empty(); }
+
+    /* Adds one interaction to the list */
+    template<std::size_t numAtoms>
+    void push_back(const int parameterType, const std::array<int, numAtoms>& atoms)
+    {
+        const std::size_t oldSize = iatoms.size();
+        iatoms.resize(iatoms.size() + 1 + numAtoms);
+        iatoms[oldSize] = parameterType;
+        for (std::size_t i = 0; i < numAtoms; i++)
+        {
+            iatoms[oldSize + 1 + i] = atoms[i];
+        }
+    }
+
+    /* Adds one interaction to the list */
+    void push_back(const int parameterType, const int numAtoms, const int* atoms)
+    {
+        const std::size_t oldSize = iatoms.size();
+        iatoms.resize(iatoms.size() + 1 + numAtoms);
+        iatoms[oldSize] = parameterType;
+        for (int i = 0; i < numAtoms; i++)
+        {
+            iatoms[oldSize + 1 + i] = atoms[i];
+        }
+    }
+
+    /* Appends \p ilist at the back of the list */
+    void append(const InteractionList& ilist)
     {
-        return iatoms.size();
+        iatoms.insert(iatoms.end(), ilist.iatoms.begin(), ilist.iatoms.end());
     }
 
+    /* Clears the list */
+    void clear() { iatoms.clear(); }
+
     /* List of interactions, see explanation further down */
     std::vector<int> iatoms;
 };
@@ -191,35 +269,23 @@ struct InteractionList
  *
  * TODO: Consider only including entries in use instead of all F_NRE
  */
-typedef std::array<InteractionList, F_NRE> InteractionLists;
+using InteractionLists = std::array<InteractionList, F_NRE>;
 
-/* Deprecated list of listed interactions.
- *
- * The nonperturbed/perturbed interactions are now separated (sorted) in the
- * ilist, such that the first 0..(nr_nonperturbed-1) ones are exactly that, and
- * the remaining ones from nr_nonperturbed..(nr-1) are perturbed bonded
- * interactions.
- */
+/* Deprecated list of listed interactions */
 struct t_ilist
 {
     /* Returns the total number of elements in iatoms */
-    int size() const
-    {
-        return nr;
-    }
+    int size() const { return nr; }
+
+    /* Returns whether the list is empty */
+    bool empty() const { return nr == 0; }
 
     int      nr;
-    int      nr_nonperturbed;
-    t_iatom *iatoms;
+    t_iatom* iatoms;
     int      nalloc;
 };
 
-/* TODO: Replace t_ilist in gmx_localtop_t by InteractionList.
- *       The nr_nonperturbed functionality needs to be ported.
- *       Remove t_topology.
- *       Remove t_ilist and remove templating on list type
- *       in mshift.cpp, constr.cpp, vsite.cpp and domdec_topology.cpp.
- */
+/* TODO: Remove t_ilist and remove templating on list type in mshift.cpp */
 
 /*
  * The structs InteractionList and t_ilist defines a list of atoms with their interactions.
@@ -240,49 +306,111 @@ struct t_ilist
  *      identifier is an index in a params[] and functype[] array.
  */
 
-typedef struct
+/*! \brief Type for returning a list of InteractionList references
+ *
+ * TODO: Remove when the function type is made part of InteractionList
+ */
+struct InteractionListHandle
+{
+    const int               functionType; //!< The function type
+    const std::vector<int>& iatoms;       //!< Reference to interaction list
+};
+
+/*! \brief Returns a list of all non-empty InteractionList entries with any of the interaction flags in \p flags set
+ *
+ * \param[in] ilists  Set of interaction lists
+ * \param[in] flags   Bit mask with one or more IF_... bits set
+ */
+static inline std::vector<InteractionListHandle> extractILists(const InteractionLists& ilists, int flags)
+{
+    std::vector<InteractionListHandle> handles;
+    for (size_t ftype = 0; ftype < ilists.size(); ftype++)
+    {
+        if ((interaction_function[ftype].flags & flags) && !ilists[ftype].empty())
+        {
+            handles.push_back({ static_cast<int>(ftype), ilists[ftype].iatoms });
+        }
+    }
+    return handles;
+}
+
+/*! \brief Returns the stride for the iatoms array in \p ilistHandle
+ *
+ * \param[in] ilistHandle  The ilist to return the stride for
+ */
+static inline int ilistStride(const InteractionListHandle& ilistHandle)
+{
+    return 1 + NRAL(ilistHandle.functionType);
+}
+
+struct gmx_cmapdata_t
 {
-    real *cmap; /* Has length 4*grid_spacing*grid_spacing, */
+    std::vector<real> cmap; /* Has length 4*grid_spacing*grid_spacing, */
     /* there are 4 entries for each cmap type (V,dVdx,dVdy,d2dVdxdy) */
-} gmx_cmapdata_t;
+};
 
-typedef struct gmx_cmap_t
+struct gmx_cmap_t
 {
-    int             ngrid;        /* Number of allocated cmap (cmapdata_t ) grids */
-    int             grid_spacing; /* Grid spacing */
-    gmx_cmapdata_t *cmapdata;     /* Pointer to grid with actual, pre-interpolated data */
-} gmx_cmap_t;
+    int                         grid_spacing = 0; /* Grid spacing */
+    std::vector<gmx_cmapdata_t> cmapdata; /* Lists of grids with actual, pre-interpolated data */
+};
 
 
-typedef struct gmx_ffparams_t
+enum
 {
-    int         ntypes;
-    int         atnr;
-    t_functype *functype;
-    t_iparams  *iparams;
-    double      reppow;    /* The repulsion power for VdW: C12*r^-reppow   */
-    real        fudgeQQ;   /* The scaling factor for Coulomb 1-4: f*q1*q2  */
-    gmx_cmap_t  cmap_grid; /* The dihedral correction maps                 */
-} gmx_ffparams_t;
-
-enum {
-    ilsortUNKNOWN, ilsortNO_FE, ilsortFE_UNSORTED, ilsortFE_SORTED
+    ilsortUNKNOWN,
+    ilsortNO_FE,
+    ilsortFE_SORTED
 };
 
-typedef struct t_idef
+/* Struct with list of interaction parameters and lists of interactions
+ *
+ * TODO: Convert to a proper class with private data members so we can
+ * ensure that the free-energy sorting and sorting setting is consistent.
+ */
+class InteractionDefinitions
+{
+public:
+    /* Constructor
+     *
+     * \param[in] ffparams  The interaction parameters, the lifetime of the created object should not exceed the lifetime of the passed parameters
+     */
+    InteractionDefinitions(const gmx_ffparams_t& ffparams);
+
+    // Clears data not read in from ffparams
+    void clear();
+
+    // The interaction parameters
+    const std::vector<t_iparams>& iparams;
+    // The function type per type
+    const std::vector<int>& functype;
+    // Position restraint interaction parameters
+    std::vector<t_iparams> iparams_posres;
+    // Flat-bottomed position restraint parameters
+    std::vector<t_iparams> iparams_fbposres;
+    // The list of interactions for each type. Note that some, such as LJ and COUL will have 0 entries.
+    std::array<InteractionList, F_NRE> il;
+    /* The number of non-perturbed interactions at the start of each entry in il */
+    std::array<int, F_NRE> numNonperturbedInteractions;
+    // The sorting state of interaction in il
+    int ilsort = ilsortUNKNOWN;
+    // The dihedral correction maps
+    gmx_cmap_t cmap_grid;
+};
+
+/* Deprecated interation definitions, used in t_topology */
+struct t_idef
 {
     int         ntypes;
     int         atnr;
-    t_functype *functype;
-    t_iparams  *iparams;
+    t_functypefunctype;
+    t_iparams*  iparams;
     real        fudgeQQ;
-    gmx_cmap_t  cmap_grid;
-    t_iparams  *iparams_posres, *iparams_fbposres;
-    int         iparams_posres_nalloc, iparams_fbposres_nalloc;
+    t_iparams * iparams_posres, *iparams_fbposres;
 
-    t_ilist     il[F_NRE];
-    int         ilsort;
-} t_idef;
+    t_ilist il[F_NRE];
+    int     ilsort;
+};
 
 /*
  * The struct t_idef defines all the interactions for the complete
@@ -318,14 +446,35 @@ typedef struct t_idef
  *      The state of the sorting of il, values are provided above.
  */
 
-void pr_iparams(FILE *fp, t_functype ftype, const t_iparams *iparams);
-void pr_ilist(FILE *fp, int indent, const char *title,
-              const t_functype *functype, const InteractionList &ilist,
-              gmx_bool bShowNumbers,
-              gmx_bool bShowParameters, const t_iparams *iparams);
-void pr_ffparams(FILE *fp, int indent, const char *title,
-                 const gmx_ffparams_t *ffparams, gmx_bool bShowNumbers);
-void pr_idef(FILE *fp, int indent, const char *title, const t_idef *idef,
-             gmx_bool bShowNumbers, gmx_bool bShowParameters);
+namespace gmx
+{
+class TextWriter;
+} // namespace gmx
+
+void printInteractionParameters(gmx::TextWriter* writer, t_functype ftype, const t_iparams& iparams);
+void pr_iparams(FILE* fp, t_functype ftype, const t_iparams& iparams);
+void pr_ilist(FILE*                  fp,
+              int                    indent,
+              const char*            title,
+              const t_functype*      functype,
+              const InteractionList& ilist,
+              bool                   bShowNumbers,
+              bool                   bShowParameters,
+              const t_iparams*       iparams);
+void pr_idef(FILE* fp, int indent, const char* title, const t_idef* idef, bool bShowNumbers, bool bShowParameters);
+
+/*! \brief
+ * Properly initialize idef struct.
+ *
+ * \param[in] idef Pointer to idef struct to initialize.
+ */
+void init_idef(t_idef* idef);
+
+/*! \brief
+ * Properly clean up idef struct.
+ *
+ * \param[in] idef Pointer to idef struct to clean up.
+ */
+void done_idef(t_idef* idef);
 
 #endif