src/gromacs/mdlib/tests/constr.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2018,2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief SHAKE and LINCS tests.
  37  *
  38  * \author Artem Zhmurov <zhmurov@gmail.com>
  39  * \ingroup module_mdlib
  40  */
  41
  42 #include "gmxpre.h"
  43
  44 #include "gromacs/mdlib/constr.h"
  45
  46 #include <assert.h>
  47
  48 #include <cmath>
  49
  50 #include <algorithm>
  51 #include <unordered_map>
  52 #include <vector>
  53
  54 #include <gtest/gtest.h>
  55
  56 #include "gromacs/fileio/gmxfio.h"
  57 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
  58 #include "gromacs/math/paddedvector.h"
  59 #include "gromacs/math/vec.h"
  60 #include "gromacs/math/vectypes.h"
  61 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  62 #include "gromacs/mdlib/lincs.h"
  63 #include "gromacs/mdlib/shake.h"
  64 #include "gromacs/mdtypes/commrec.h"
  65 #include "gromacs/mdtypes/inputrec.h"
  66 #include "gromacs/mdtypes/mdatom.h"
  67 #include "gromacs/pbcutil/pbc.h"
  68 #include "gromacs/topology/block.h"
  69 #include "gromacs/topology/idef.h"
  70 #include "gromacs/topology/ifunc.h"
  71 #include "gromacs/topology/topology.h"
  72 #include "gromacs/utility/smalloc.h"
  73 #include "gromacs/utility/stringutil.h"
  74 #include "gromacs/utility/unique_cptr.h"
  75
  76 #include "testutils/refdata.h"
  77 #include "testutils/testasserts.h"
  78
  79 namespace gmx
  80 {
  81 namespace test
  82 {
  83
  84 /*! \brief
  85  * Constraints test data structure.
  86  *
  87  * Structure to collect all the necessary data, including system coordinates and topology,
  88  * constraints information, etc. The structure can be reset and reused.
  89  */
  90 struct ConstraintsTestData
  91 {
  92     public:
  93         std::string           title_;                   //!< Human-friendly name for a system
  94         int                   nAtom_;                   //!< Number of atoms
  95         gmx_mtop_t            mtop_;                    //!< Topology
  96         std::vector<real>     masses_;                  //!< Masses
  97         std::vector<real>     invmass_;                 //!< Inverse masses
  98         t_commrec             cr_;                      //!< Communication record
  99         t_inputrec            ir_;                      //!< Input record (info that usually in .mdp file)
 100         t_idef                idef_;                    //!< Local topology
 101         t_mdatoms             md_;                      //!< MD atoms
 102         gmx_multisim_t        ms_;                      //!< Multisim data
 103         t_nrnb                nrnb_;                    //!< Computational time array (normally used in benchmarks)
 104
 105         real                  invdt_;                   //!< Inverse timestep
 106         int                   nflexcon_   = 0;          //!< Number of flexible constraints
 107         bool                  computeVirial_;           //!< Whether the virial should be computed
 108         tensor                virialScaled_;            //!< Scaled virial
 109         tensor                virialScaledRef_;         //!< Scaled virial (reference values)
 110         bool                  compute_dHdLambda_;       //!< If the free energy is computed
 111         real                  dHdLambda_;               //!< For free energy computation
 112         real                  dHdLambdaRef_;            //!< For free energy computation (reference value)
 113
 114         PaddedVector<RVec>    x_;                       //!< Coordinates before the timestep
 115         PaddedVector<RVec>    xPrime_;                  //!< Coordinates after timestep, output for the constraints
 116         PaddedVector<RVec>    xPrime0_;                 //!< Backup for coordinates (for reset)
 117         PaddedVector<RVec>    xPrime2_;                 //!< Intermediate set of coordinates used by LINCS and
 118                                                         //!< SHAKE for different purposes
 119         PaddedVector<RVec>    v_;                       //!< Velocities
 120         PaddedVector<RVec>    v0_;                      //!< Backup for velocities (for reset)
 121
 122         // Fields to store constraints data for testing
 123         std::vector<int>      constraints_;             //!< Constraints data (type1-i1-j1-type2-i2-j2-...)
 124         std::vector<real>     constraintsR0_;           //!< Target lengths for all constraint types
 125
 126         /*! \brief
 127          * Constructor for the object with all parameters and variables needed by constraints algorithms.
 128          *
 129          * This constructor assembles stubs for all the data structures, required to initialize
 130          * and apply LINCS and SHAKE constraints. The coordinates and velocities before constraining
 131          * are saved to allow for reset. The constraints data are stored for testing after constraints
 132          * were applied.
 133          *
 134          * \param[in]  title              Human-friendly name of the system.
 135          * \param[in]  nAtom              Number of atoms in the system.
 136          * \param[in]  masses             Atom masses. Size of this vector should be equal to nAtom.
 137          * \param[in]  constraints        List of constraints, organized in triples of integers.
 138          *                                First integer is the index of type for a constraint, second
 139          *                                and third are the indices of constrained atoms. The types
 140          *                                of constraints should be sequential but not necessarily
 141          *                                start from zero (which is the way they normally are in
 142          *                                GROMACS).
 143          * \param[in]  constraintsR0      Target values for bond lengths for bonds of each type. The
 144          *                                size of this vector should be equal to the total number of
 145          *                                unique types in constraints vector.
 146          * \param[in]  computeVirial      Whether the virial should be computed.
 147          * \param[in]  virialScaledRef    Reference values for scaled virial tensor.
 148          * \param[in]  compute_dHdLambda  Whether free energy should be computed.
 149          * \param[in]  dHdLambdaRef       Reference value for dHdLambda.
 150          * \param[in]  initialTime        Initial time.
 151          * \param[in]  timestep           Timestep.
 152          * \param[in]  x                  Coordinates before integration step.
 153          * \param[in]  xPrime             Coordinates after integration step, but before constraining.
 154          * \param[in]  v                  Velocities before constraining.
 155          * \param[in]  shakeTolerance     Target tolerance for SHAKE.
 156          * \param[in]  shakeUseSOR        Use successive over-relaxation method for SHAKE iterations.
 157          *                                The general formula is:
 158          *                                   x_n+1 = (1-omega)*x_n + omega*f(x_n),
 159          *                                where omega = 1 if SOR is off and may be < 1 if SOR is on.
 160          * \param[in]  nLincsIter         Number of iterations used to compute the inverse matrix.
 161          * \param[in]  nProjOrder         The order for algorithm that adjusts the direction of the
 162          *                                bond after constraints are applied.
 163          * \param[in]  lincsWarnAngle     The threshold value for the change in bond angle. When
 164          *                                exceeded the program will issue a warning.
 165          *
 166          */
 167         ConstraintsTestData(const std::string &title,
 168                             int nAtom, std::vector<real> masses,
 169                             std::vector<int> constraints, std::vector<real> constraintsR0,
 170                             bool computeVirial, tensor virialScaledRef,
 171                             bool compute_dHdLambda, float dHdLambdaRef,
 172                             real initialTime, real timestep,
 173                             const std::vector<RVec> &x, const std::vector<RVec> &xPrime, const std::vector<RVec> &v,
 174                             real shakeTolerance, gmx_bool shakeUseSOR,
 175                             int nLincsIter, int nProjOrder, real lincsWarnAngle)
 176         {
 177             title_ = title;   // Human-friendly name of the system
 178             nAtom_ = nAtom;   // Number of atoms
 179
 180             // Masses of atoms
 181             masses_ = masses;
 182             invmass_.resize(nAtom); // Vector of inverse masses
 183
 184             for (int i = 0; i < nAtom; i++)
 185             {
 186                 invmass_[i] = 1.0/masses.at(i);
 187             }
 188
 189             // Saving constraints to check if they are satisfied after algorithm was applied
 190             constraints_   = constraints;   // Constraints indexes (in type-i-j format)
 191             constraintsR0_ = constraintsR0; // Equilibrium distances for each type of constraint
 192
 193             invdt_  = 1.0/timestep;         // Inverse timestep
 194
 195             // Communication record
 196             cr_.nnodes = 1;
 197             cr_.dd     = nullptr;
 198
 199             // Multisim data
 200             ms_.sim  = 0;
 201             ms_.nsim = 1;
 202
 203             // Input record - data that usually comes from configuration file (.mdp)
 204             ir_.efep           = 0;
 205             ir_.init_t         = initialTime;
 206             ir_.delta_t        = timestep;
 207             ir_.eI             = 0;
 208
 209             // MD atoms data
 210             md_.nMassPerturbed = 0;
 211             md_.lambda         = 0.0;
 212             md_.invmass        = invmass_.data();
 213             md_.nr             = nAtom;
 214             md_.homenr         = nAtom;
 215
 216             // Virial evaluation
 217             computeVirial_ = computeVirial;
 218             if (computeVirial)
 219             {
 220                 for (int i = 0; i < DIM; i++)
 221                 {
 222                     for (int j = 0; j < DIM; j++)
 223                     {
 224                         virialScaled_[i][j]    = 0;
 225                         virialScaledRef_[i][j] = virialScaledRef[i][j];
 226                     }
 227                 }
 228             }
 229
 230
 231             // Free energy evaluation
 232             compute_dHdLambda_ = compute_dHdLambda;
 233             dHdLambda_         = 0;
 234             if (compute_dHdLambda_)
 235             {
 236                 ir_.efep                    = efepYES;
 237                 dHdLambdaRef_               = dHdLambdaRef;
 238             }
 239             else
 240             {
 241                 ir_.efep                    = efepNO;
 242                 dHdLambdaRef_               = 0;
 243             }
 244
 245             // Constraints and their parameters (local topology)
 246             for (int i = 0; i < F_NRE; i++)
 247             {
 248                 idef_.il[i].nr = 0;
 249             }
 250             idef_.il[F_CONSTR].nr   = constraints.size();
 251
 252             snew(idef_.il[F_CONSTR].iatoms, constraints.size());
 253             int maxType = 0;
 254             for (unsigned i = 0; i < constraints.size(); i++)
 255             {
 256                 if (i % 3 == 0)
 257                 {
 258                     if (maxType < constraints.at(i))
 259                     {
 260                         maxType = constraints.at(i);
 261                     }
 262                 }
 263                 idef_.il[F_CONSTR].iatoms[i] = constraints.at(i);
 264             }
 265             snew(idef_.iparams, maxType + 1);
 266             for (unsigned i = 0; i < constraints.size()/3; i++)
 267             {
 268                 idef_.iparams[constraints.at(3*i)].constr.dA = constraintsR0.at(constraints.at(3*i));
 269                 idef_.iparams[constraints.at(3*i)].constr.dB = constraintsR0.at(constraints.at(3*i));
 270             }
 271
 272             // Constraints and their parameters (global topology)
 273             InteractionList interactionList;
 274             interactionList.iatoms.resize(constraints.size());
 275             for (unsigned i = 0; i < constraints.size(); i++)
 276             {
 277                 interactionList.iatoms.at(i) = constraints.at(i);
 278             }
 279             InteractionList interactionListEmpty;
 280             interactionListEmpty.iatoms.resize(0);
 281
 282             gmx_moltype_t molType;
 283             molType.atoms.nr             = nAtom;
 284             molType.ilist.at(F_CONSTR)   = interactionList;
 285             molType.ilist.at(F_CONSTRNC) = interactionListEmpty;
 286             mtop_.moltype.push_back(molType);
 287
 288             gmx_molblock_t molBlock;
 289             molBlock.type = 0;
 290             molBlock.nmol = 1;
 291             mtop_.molblock.push_back(molBlock);
 292
 293             mtop_.natoms = nAtom;
 294             mtop_.ffparams.iparams.resize(maxType + 1);
 295             for (int i = 0; i <= maxType; i++)
 296             {
 297                 mtop_.ffparams.iparams.at(i) = idef_.iparams[i];
 298             }
 299             mtop_.bIntermolecularInteractions = false;
 300
 301             // Coordinates and velocities
 302             x_.resizeWithPadding(nAtom_);
 303             xPrime_.resizeWithPadding(nAtom_);
 304             xPrime0_.resizeWithPadding(nAtom_);
 305             xPrime2_.resizeWithPadding(nAtom_);
 306
 307             v_.resizeWithPadding(nAtom_);
 308             v0_.resizeWithPadding(nAtom_);
 309
 310             std::copy(x.begin(), x.end(), x_.begin());
 311             std::copy(xPrime.begin(), xPrime.end(), xPrime_.begin());
 312             std::copy(xPrime.begin(), xPrime.end(), xPrime0_.begin());
 313             std::copy(xPrime.begin(), xPrime.end(), xPrime2_.begin());
 314
 315             std::copy(v.begin(), v.end(), v_.begin());
 316             std::copy(v.begin(), v.end(), v0_.begin());
 317
 318             // SHAKE-specific parameters
 319             ir_.shake_tol           = shakeTolerance;
 320             ir_.bShakeSOR           = shakeUseSOR;
 321
 322             // LINCS-specific parameters
 323             ir_.nLincsIter          = nLincsIter;
 324             ir_.nProjOrder          = nProjOrder;
 325             ir_.LincsWarnAngle      = lincsWarnAngle;
 326         }
 327
 328         /*! \brief
 329          * Reset the data structure so it can be reused.
 330          *
 331          * Set the coordinates and velocities back to their values before
 332          * constraining. The scaled virial tensor and dHdLambda are zeroed.
 333          *
 334          */
 335         void reset()
 336         {
 337             xPrime_  = xPrime0_;
 338             xPrime2_ = xPrime0_;
 339             v_       = v0_;
 340
 341             if (computeVirial_)
 342             {
 343                 for (int i = 0; i < DIM; i++)
 344                 {
 345                     for (int j = 0; j < DIM; j++)
 346                     {
 347                         virialScaled_[i][j] = 0;
 348                     }
 349                 }
 350             }
 351             dHdLambda_         = 0;
 352         }
 353
 354         /*! \brief
 355          * Cleaning up the memory.
 356          */
 357         ~ConstraintsTestData()
 358         {
 359             sfree(idef_.il[F_CONSTR].iatoms);
 360             sfree(idef_.iparams);
 361         }
 362
 363 };
 364
 365 /*! \brief The two-dimensional parameter space for test.
 366  *
 367  * The test will run for all possible combinations of accessible
 368  * values of the:
 369  * 1. PBC setup ("PBCNONE" or "PBCXYZ")
 370  * 2. The algorithm ("SHAKE" or "LINCS").
 371  */
 372 typedef std::tuple<std::string, std::string> ConstraintsTestParameters;
 373
 374 /*! \brief Test fixture for constraints.
 375  *
 376  * The fixture uses following test systems:
 377  * 1. Two atoms, connected with one constraint (e.g. NH).
 378  * 2. Three atoms, connected consequently with two constraints (e.g. CH2).
 379  * 3. Three atoms, constrained to the fourth atom (e.g. CH3).
 380  * 4. Four atoms, connected by two independent constraints.
 381  * 5. Three atoms, connected by three constraints in a triangle
 382  *      (e.g. H2O with constrained H-O-H angle).
 383  * 6. Four atoms, connected by three consequential constraints.
 384  *
 385  * For all systems, the final lengths of the constraints are tested against the
 386  * reference values, the direction of each constraint is checked.
 387  * Test also verifies that the center of mass has not been
 388  * shifted by the constraints and that its velocity has not changed.
 389  * For some systems, the value for scaled virial tensor is checked against
 390  * pre-computed data.
 391  */
 392 class ConstraintsTest : public ::testing::TestWithParam<ConstraintsTestParameters>
 393 {
 394     public:
 395         //! PBC setups
 396         std::unordered_map <std::string, t_pbc>                                             pbcs_;
 397         //! Algorithms (SHAKE and LINCS)
 398         std::unordered_map <std::string, void(*)(ConstraintsTestData *testData, t_pbc pbc)> algorithms_;
 399
 400         /*! \brief Test setup function.
 401          *
 402          * Setting up the pbcs and algorithms. Note, that corresponding string keywords
 403          * have to be explicitly added at the end of this file when the tests are called.
 404          *
 405          */
 406         void SetUp() override
 407         {
 408
 409             //
 410             // PBC initialization
 411             //
 412             t_pbc pbc;
 413
 414             // Infinitely small box
 415             matrix boxNone = { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} };
 416             set_pbc(&pbc, epbcNONE, boxNone);
 417             pbcs_["PBCNone"] = pbc;
 418
 419             // Rectangular box
 420             matrix boxXyz = { {10.0, 0.0, 0.0}, {0.0, 20.0, 0.0}, {0.0, 0.0, 15.0} };
 421             set_pbc(&pbc, epbcXYZ, boxXyz);
 422             pbcs_["PBCXYZ"] = pbc;
 423
 424             //
 425             // Algorithms
 426             //
 427             // SHAKE
 428             algorithms_["SHAKE"] = applyShake;
 429             // LINCS
 430             algorithms_["LINCS"] = applyLincs;
 431
 432         }
 433
 434         /*! \brief
 435          * Initialize and apply SHAKE constraints.
 436          *
 437          * \param[in] testData        Test data structure.
 438          * \param[in] pbc             Periodic boundary data (not used in SHAKE).
 439          */
 440         static void applyShake(ConstraintsTestData *testData, t_pbc gmx_unused pbc)
 441         {
 442             shakedata* shaked = shake_init();
 443             make_shake_sblock_serial(shaked, &testData->idef_, testData->md_);
 444             bool       success = constrain_shake(
 445                         nullptr,
 446                         shaked,
 447                         testData->invmass_.data(),
 448                         testData->idef_,
 449                         testData->ir_,
 450                         as_rvec_array(testData->x_.data()),
 451                         as_rvec_array(testData->xPrime_.data()),
 452                         as_rvec_array(testData->xPrime2_.data()),
 453                         &testData->nrnb_,
 454                         testData->md_.lambda,
 455                         &testData->dHdLambda_,
 456                         testData->invdt_,
 457                         as_rvec_array(testData->v_.data()),
 458                         testData->computeVirial_,
 459                         testData->virialScaled_,
 460                         false,
 461                         gmx::ConstraintVariable::Positions);
 462             EXPECT_TRUE(success) << "Test failed with a false return value in SHAKE.";
 463             done_shake(shaked);
 464         }
 465
 466         /*! \brief
 467          * Initialize and apply LINCS constraints.
 468          *
 469          * \param[in] testData        Test data structure.
 470          * \param[in] pbc             Periodic boundary data.
 471          */
 472         static void applyLincs(ConstraintsTestData *testData, t_pbc pbc)
 473         {
 474
 475             Lincs                *lincsd;
 476             int                   maxwarn         = 100;
 477             int                   warncount_lincs = 0;
 478             gmx_omp_nthreads_set(emntLINCS, 1);
 479
 480             // Make blocka structure for faster LINCS setup
 481             std::vector<t_blocka> at2con_mt;
 482             at2con_mt.reserve(testData->mtop_.moltype.size());
 483             for (const gmx_moltype_t &moltype : testData->mtop_.moltype)
 484             {
 485                 // This function is in constr.cpp
 486                 at2con_mt.push_back(make_at2con(moltype,
 487                                                 testData->mtop_.ffparams.iparams,
 488                                                 flexibleConstraintTreatment(EI_DYNAMICS(testData->ir_.eI))));
 489             }
 490             // Initialize LINCS
 491             lincsd = init_lincs(nullptr, testData->mtop_,
 492                                 testData->nflexcon_, at2con_mt,
 493                                 false,
 494                                 testData->ir_.nLincsIter, testData->ir_.nProjOrder);
 495             set_lincs(testData->idef_, testData->md_, EI_DYNAMICS(testData->ir_.eI), &testData->cr_, lincsd);
 496
 497             // Evaluate constraints
 498             bool sucess = constrain_lincs(false, testData->ir_, 0, lincsd, testData->md_,
 499                                           &testData->cr_,
 500                                           &testData->ms_,
 501                                           as_rvec_array(testData->x_.data()),
 502                                           as_rvec_array(testData->xPrime_.data()),
 503                                           as_rvec_array(testData->xPrime2_.data()),
 504                                           pbc.box, &pbc,
 505                                           testData->md_.lambda, &testData->dHdLambda_,
 506                                           testData->invdt_,
 507                                           as_rvec_array(testData->v_.data()),
 508                                           testData->computeVirial_, testData->virialScaled_,
 509                                           gmx::ConstraintVariable::Positions,
 510                                           &testData->nrnb_,
 511                                           maxwarn, &warncount_lincs);
 512             EXPECT_TRUE(sucess) << "Test failed with a false return value in LINCS.";
 513             EXPECT_EQ(warncount_lincs, 0) << "There were warnings in LINCS.";
 514             for (unsigned int i = 0; i < testData->mtop_.moltype.size(); i++)
 515             {
 516                 sfree(at2con_mt.at(i).index);
 517                 sfree(at2con_mt.at(i).a);
 518             }
 519             done_lincs(lincsd);
 520         }
 521
 522         /*! \brief
 523          * The test on the final length of constrained bonds.
 524          *
 525          * Goes through all the constraints and checks if the final length of all the constraints is equal
 526          * to the target length with provided tolerance.
 527          *
 528          * \param[in] tolerance       Allowed tolerance in final lengths.
 529          * \param[in] testData        Test data structure.
 530          * \param[in] pbc             Periodic boundary data.
 531          */
 532         void checkConstrainsLength(FloatingPointTolerance tolerance, const ConstraintsTestData &testData, t_pbc pbc)
 533         {
 534
 535             // Test if all the constraints are satisfied
 536             for (unsigned c = 0; c < testData.constraints_.size()/3; c++)
 537             {
 538                 real r0 = testData.constraintsR0_.at(testData.constraints_.at(3*c));
 539                 int  i  = testData.constraints_.at(3*c + 1);
 540                 int  j  = testData.constraints_.at(3*c + 2);
 541                 RVec xij0, xij1;
 542                 real d0, d1;
 543                 if (pbc.ePBC == epbcXYZ)
 544                 {
 545                     pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
 546                     pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
 547                 }
 548                 else
 549                 {
 550                     rvec_sub(testData.x_[i], testData.x_[j], xij0);
 551                     rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
 552                 }
 553                 d0 = norm(xij0);
 554                 d1 = norm(xij1);
 555                 EXPECT_REAL_EQ_TOL(r0, d1, tolerance) << gmx::formatString(
 556                         "rij = %f, which is not equal to r0 = %f for constraint #%u, between atoms %d and %d"
 557                         " (before constraining rij was %f).", d1, r0, c, i, j, d0);
 558             }
 559         }
 560
 561         /*! \brief
 562          * The test on the final length of constrained bonds.
 563          *
 564          * Goes through all the constraints and checks if the direction of constraint has not changed
 565          * by the algorithm (i.e. the constraints algorithm arrived to the solution that is closest
 566          * to the initial system conformation).
 567          *
 568          * \param[in] testData        Test data structure.
 569          * \param[in] pbc             Periodic boundary data.
 570          */
 571         void checkConstrainsDirection(const ConstraintsTestData &testData, t_pbc pbc)
 572         {
 573
 574             for (unsigned c = 0; c < testData.constraints_.size()/3; c++)
 575             {
 576                 int  i  = testData.constraints_.at(3*c + 1);
 577                 int  j  = testData.constraints_.at(3*c + 2);
 578                 RVec xij0, xij1;
 579                 if (pbc.ePBC == epbcXYZ)
 580                 {
 581                     pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
 582                     pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
 583                 }
 584                 else
 585                 {
 586                     rvec_sub(testData.x_[i], testData.x_[j], xij0);
 587                     rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
 588                 }
 589
 590                 real dot = xij0.dot(xij1);
 591
 592                 EXPECT_GE(dot, 0.0) << gmx::formatString(
 593                         "The constraint %u changed direction. Constraining algorithm might have returned the wrong root "
 594                         "of the constraints equation.", c);
 595
 596             }
 597         }
 598
 599         /*! \brief
 600          * The test on the coordinates of the center of the mass (COM) of the system.
 601          *
 602          * Checks if the center of mass has not been shifted by the constraints. Note,
 603          * that this test does not take into account the periodic boundary conditions.
 604          * Hence it will not work should the constraints decide to move atoms across
 605          * PBC borders.
 606          *
 607          * \param[in] tolerance       Allowed tolerance in COM coordinates.
 608          * \param[in] testData        Test data structure.
 609          */
 610         void checkCOMCoordinates(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
 611         {
 612
 613             RVec comPrime0({0.0, 0.0, 0.0});
 614             RVec comPrime({0.0, 0.0, 0.0});
 615             for (int i = 0; i < testData.nAtom_; i++)
 616             {
 617                 comPrime0 += testData.masses_[i]*testData.xPrime0_[i];
 618                 comPrime  += testData.masses_[i]*testData.xPrime_[i];
 619             }
 620
 621             comPrime0 /= testData.nAtom_;
 622             comPrime  /= testData.nAtom_;
 623
 624             EXPECT_REAL_EQ_TOL(comPrime[XX], comPrime0[XX], tolerance)
 625             << "Center of mass was shifted by constraints in x-direction.";
 626             EXPECT_REAL_EQ_TOL(comPrime[YY], comPrime0[YY], tolerance)
 627             << "Center of mass was shifted by constraints in y-direction.";
 628             EXPECT_REAL_EQ_TOL(comPrime[ZZ], comPrime0[ZZ], tolerance)
 629             << "Center of mass was shifted by constraints in z-direction.";
 630
 631         }
 632
 633         /*! \brief
 634          * The test on the velocity of the center of the mass (COM) of the system.
 635          *
 636          * Checks if the velocity of the center of mass has not changed.
 637          *
 638          * \param[in] tolerance       Allowed tolerance in COM velocity components.
 639          * \param[in] testData        Test data structure.
 640          */
 641         void checkCOMVelocity(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
 642         {
 643
 644             RVec comV0({0.0, 0.0, 0.0});
 645             RVec comV({0.0, 0.0, 0.0});
 646             for (int i = 0; i < testData.nAtom_; i++)
 647             {
 648                 comV0 += testData.masses_[i]*testData.v0_[i];
 649                 comV  += testData.masses_[i]*testData.v_[i];
 650             }
 651             comV0 /= testData.nAtom_;
 652             comV  /= testData.nAtom_;
 653
 654             EXPECT_REAL_EQ_TOL(comV[XX], comV0[XX], tolerance)
 655             << "Velocity of the center of mass in x-direction has been changed by constraints.";
 656             EXPECT_REAL_EQ_TOL(comV[YY], comV0[YY], tolerance)
 657             << "Velocity of the center of mass in y-direction has been changed by constraints.";
 658             EXPECT_REAL_EQ_TOL(comV[ZZ], comV0[ZZ], tolerance)
 659             << "Velocity of the center of mass in z-direction has been changed by constraints.";
 660         }
 661
 662         /*! \brief
 663          * The test on the final coordinates (not used).
 664          *
 665          * Goes through all atoms and checks if the final positions correspond to the
 666          * provided reference set of coordinates.
 667          *
 668          * \param[in] xPrimeRef       The reference set of coordinates.
 669          * \param[in] tolerance       Tolerance for the coordinates test.
 670          * \param[in] testData        Test data structure.
 671          */
 672         void checkFinalCoordinates(std::vector<RVec> xPrimeRef, FloatingPointTolerance tolerance,
 673                                    const ConstraintsTestData &testData)
 674         {
 675             for (int i = 0; i < testData.nAtom_; i++)
 676             {
 677                 for (int d = 0; d < DIM; d++)
 678                 {
 679                     EXPECT_REAL_EQ_TOL(xPrimeRef.at(i)[d], testData.xPrime_[i][d], tolerance) << gmx::formatString(
 680                             "Coordinates after constrains were applied differ from these in the reference set for atom #%d.", i);
 681                 }
 682             }
 683         }
 684
 685         /*! \brief
 686          * The test of virial tensor.
 687          *
 688          * Checks if the values in the scaled virial tensor are equal to pre-computed values.
 689          *
 690          * \param[in] tolerance       Tolerance for the tensor values.
 691          * \param[in] testData        Test data structure.
 692          */
 693         void checkVirialTensor(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
 694         {
 695             for (int i = 0; i < DIM; i++)
 696             {
 697                 for (int j = 0; j < DIM; j++)
 698                 {
 699                     EXPECT_REAL_EQ_TOL(testData.virialScaledRef_[i][j], testData.virialScaled_[i][j],
 700                                        tolerance) << gmx::formatString(
 701                             "Values in virial tensor at [%d][%d] are not within the tolerance from reference value.", i, j);
 702                 }
 703             }
 704         }
 705
 706         /*! \brief
 707          * The test for FEP (not used).
 708          *
 709          * Checks if the value of dH/dLambda is equal to the reference value.
 710          * \todo Add tests for dHdLambda values.
 711          *
 712          * \param[in] dHdLambdaRef    Reference value.
 713          * \param[in] tolerance       Tolerance.
 714          * \param[in] testData        Test data structure.
 715          */
 716         void checkFEP(const real dHdLambdaRef, FloatingPointTolerance tolerance,
 717                       const ConstraintsTestData &testData)
 718         {
 719             EXPECT_REAL_EQ_TOL(dHdLambdaRef, testData.dHdLambda_, tolerance) <<
 720             "Computed value for dV/dLambda is not equal to the reference value. ";
 721         }
 722
 723
 724
 725 };
 726
 727 TEST_P(ConstraintsTest, SingleConstraint){
 728     std::string       title = "one constraint (e.g. OH)";
 729     int               nAtom = 2;
 730
 731     std::vector<real> masses        = {1.0, 12.0};
 732     std::vector<int>  constraints   = {0, 0, 1};
 733     std::vector<real> constraintsR0 = {0.1};
 734
 735     real              oneTenthOverSqrtTwo    = 0.1_real / std::sqrt(2.0_real);
 736
 737     std::vector<RVec> x = {{ 0.0, oneTenthOverSqrtTwo, 0.0 },
 738                            { oneTenthOverSqrtTwo, 0.0, 0.0 }};
 739     std::vector<RVec> xPrime = {{ 0.01,  0.08,  0.01 },
 740                                 { 0.06,  0.01, -0.01 }};
 741     std::vector<RVec> v = {{ 1.0, 2.0, 3.0 },
 742                            { 3.0, 2.0, 1.0 }};
 743
 744     tensor            virialScaledRef = {{-5.58e-04,  5.58e-04, 0.00e+00 },
 745                                          { 5.58e-04, -5.58e-04, 0.00e+00 },
 746                                          { 0.00e+00,  0.00e+00, 0.00e+00 } };
 747
 748     real              shakeTolerance         = 0.0001;
 749     gmx_bool          shakeUseSOR            = false;
 750
 751     int               lincsNIter             = 1;
 752     int               lincsNProjOrder        = 4;
 753     real              lincsWarnAngle         = 30.0;
 754
 755     std::unique_ptr<ConstraintsTestData>   testData = std::make_unique<ConstraintsTestData>
 756             (title, nAtom, masses,
 757             constraints, constraintsR0,
 758             true, virialScaledRef,
 759             false, 0,
 760             real(0.0), real(0.001),
 761             x, xPrime, v,
 762             shakeTolerance, shakeUseSOR,
 763             lincsNIter, lincsNProjOrder, lincsWarnAngle);
 764     std::string pbcName;
 765     std::string algorithmName;
 766     std::tie(pbcName, algorithmName) = GetParam();
 767     t_pbc                                   pbc       = pbcs_.at(pbcName);
 768
 769     // Apply constraints
 770     algorithms_.at(algorithmName)(testData.get(), pbc);
 771
 772     checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
 773     checkConstrainsDirection(*testData, pbc);
 774     checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
 775     checkCOMVelocity(absoluteTolerance(0.0001), *testData);
 776
 777     checkVirialTensor(absoluteTolerance(0.0001), *testData);
 778
 779 }
 780
 781 TEST_P(ConstraintsTest, TwoDisjointConstraints){
 782
 783     std::string       title         = "two disjoint constraints";
 784     int               nAtom         = 4;
 785     std::vector<real> masses        = {0.5, 1.0/3.0, 0.25, 1.0};
 786     std::vector<int>  constraints   = {0, 0, 1, 1, 2, 3};
 787     std::vector<real> constraintsR0 = {2.0, 1.0};
 788
 789
 790     std::vector<RVec> x = {{  2.50, -3.10, 15.70 },
 791                            {  0.51, -3.02, 15.55 },
 792                            { -0.50, -3.00, 15.20 },
 793                            { -1.51, -2.95, 15.05 }};
 794
 795     std::vector<RVec> xPrime = {{  2.50, -3.10, 15.70 },
 796                                 {  0.51, -3.02, 15.55 },
 797                                 { -0.50, -3.00, 15.20 },
 798                                 { -1.51, -2.95, 15.05 }};
 799
 800     std::vector<RVec> v = {{ 0.0, 1.0, 0.0 },
 801                            { 1.0, 0.0, 0.0 },
 802                            { 0.0, 0.0, 1.0 },
 803                            { 0.0, 0.0, 0.0 }};
 804
 805     tensor            virialScaledRef = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
 806
 807     real              shakeTolerance         = 0.0001;
 808     gmx_bool          shakeUseSOR            = false;
 809
 810     int               lincsNIter             = 1;
 811     int               lincsNProjOrder        = 4;
 812     real              lincsWarnAngle         = 30.0;
 813
 814     std::unique_ptr<ConstraintsTestData>   testData = std::make_unique<ConstraintsTestData>
 815             (title, nAtom, masses,
 816             constraints, constraintsR0,
 817             false, virialScaledRef,
 818             false, 0,
 819             real(0.0), real(0.001),
 820             x, xPrime, v,
 821             shakeTolerance, shakeUseSOR,
 822             lincsNIter, lincsNProjOrder, lincsWarnAngle);
 823
 824     std::string pbcName;
 825     std::string algorithmName;
 826     std::tie(pbcName, algorithmName) = GetParam();
 827     t_pbc                                   pbc       = pbcs_.at(pbcName);
 828
 829     // Apply constraints
 830     algorithms_.at(algorithmName)(testData.get(), pbc);
 831
 832     checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
 833     checkConstrainsDirection(*testData, pbc);
 834     checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
 835     checkCOMVelocity(absoluteTolerance(0.0001), *testData);
 836
 837 }
 838
 839 TEST_P(ConstraintsTest, ThreeSequentialConstraints){
 840
 841     std::string       title         = "three atoms, connected longitudinally (e.g. CH2)";
 842     int               nAtom         = 3;
 843     std::vector<real> masses        = {1.0, 12.0, 16.0 };
 844     std::vector<int>  constraints   = {0, 0, 1, 1, 1, 2};
 845     std::vector<real> constraintsR0 = {0.1, 0.2};
 846
 847     real              oneTenthOverSqrtTwo    = 0.1_real / std::sqrt(2.0_real);
 848     real              twoTenthsOverSqrtThree = 0.2_real / std::sqrt(3.0_real);
 849
 850     std::vector<RVec> x = {{ oneTenthOverSqrtTwo, oneTenthOverSqrtTwo, 0.0 },
 851                            { 0.0, 0.0, 0.0 },
 852                            { twoTenthsOverSqrtThree, twoTenthsOverSqrtThree, twoTenthsOverSqrtThree }};
 853
 854     std::vector<RVec> xPrime = {{  0.08, 0.07,  0.01 },
 855                                 { -0.02, 0.01, -0.02 },
 856                                 {  0.10, 0.12,  0.11 }};
 857
 858     std::vector<RVec> v = {{ 1.0, 0.0, 0.0 },
 859                            { 0.0, 1.0, 0.0 },
 860                            { 0.0, 0.0, 1.0 }};
 861
 862     tensor            virialScaledRef = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
 863
 864     real              shakeTolerance         = 0.0001;
 865     gmx_bool          shakeUseSOR            = false;
 866
 867     int               lincsNIter             = 1;
 868     int               lincsNProjOrder        = 4;
 869     real              lincsWarnAngle         = 30.0;
 870
 871     std::unique_ptr<ConstraintsTestData>   testData = std::make_unique<ConstraintsTestData>
 872             (title, nAtom, masses,
 873             constraints, constraintsR0,
 874             false, virialScaledRef,
 875             false, 0,
 876             real(0.0), real(0.001),
 877             x, xPrime, v,
 878             shakeTolerance, shakeUseSOR,
 879             lincsNIter, lincsNProjOrder, lincsWarnAngle);
 880
 881     std::string pbcName;
 882     std::string algorithmName;
 883     std::tie(pbcName, algorithmName) = GetParam();
 884     t_pbc                                   pbc       = pbcs_.at(pbcName);
 885
 886     // Apply constraints
 887     algorithms_.at(algorithmName)(testData.get(), pbc);
 888
 889     checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
 890     checkConstrainsDirection(*testData, pbc);
 891     checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
 892     checkCOMVelocity(absoluteTolerance(0.0001), *testData);
 893
 894 }
 895
 896 TEST_P(ConstraintsTest, ThreeConstraintsWithCentralAtom){
 897
 898     std::string       title         = "three atoms, connected to the central atom (e.g. CH3)";
 899     int               nAtom         = 4;
 900     std::vector<real> masses        = {12.0, 1.0, 1.0, 1.0 };
 901     std::vector<int>  constraints   = {0, 0, 1, 0, 0, 2, 0, 0, 3};
 902     std::vector<real> constraintsR0 = {0.1};
 903
 904
 905     std::vector<RVec> x = {{ 0.00,  0.00,  0.00 },
 906                            { 0.10,  0.00,  0.00 },
 907                            { 0.00, -0.10,  0.00 },
 908                            { 0.00,  0.00,  0.10 }};
 909
 910     std::vector<RVec> xPrime = {{ 0.004,  0.009, -0.010 },
 911                                 { 0.110, -0.006,  0.003 },
 912                                 {-0.007, -0.102, -0.007 },
 913                                 {-0.005,  0.011,  0.102 }};
 914
 915     std::vector<RVec> v = {{ 1.0, 0.0, 0.0 },
 916                            { 1.0, 0.0, 0.0 },
 917                            { 1.0, 0.0, 0.0 },
 918                            { 1.0, 0.0, 0.0 }};
 919
 920     tensor            virialScaledRef = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
 921
 922     real              shakeTolerance         = 0.0001;
 923     gmx_bool          shakeUseSOR            = false;
 924
 925     int               lincsNIter             = 1;
 926     int               lincsNProjOrder        = 4;
 927     real              lincsWarnAngle         = 30.0;
 928
 929     std::unique_ptr<ConstraintsTestData>   testData = std::make_unique<ConstraintsTestData>
 930             (title, nAtom, masses,
 931             constraints, constraintsR0,
 932             false, virialScaledRef,
 933             false, 0,
 934             real(0.0), real(0.001),
 935             x, xPrime, v,
 936             shakeTolerance, shakeUseSOR,
 937             lincsNIter, lincsNProjOrder, lincsWarnAngle);
 938
 939     std::string pbcName;
 940     std::string algorithmName;
 941     std::tie(pbcName, algorithmName) = GetParam();
 942     t_pbc                                   pbc       = pbcs_.at(pbcName);
 943
 944     // Apply constraints
 945     algorithms_.at(algorithmName)(testData.get(), pbc);
 946
 947     checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
 948     checkConstrainsDirection(*testData, pbc);
 949     checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
 950     checkCOMVelocity(absoluteTolerance(0.0001), *testData);
 951
 952 }
 953
 954 TEST_P(ConstraintsTest, FourSequentialConstraints){
 955
 956     std::string       title         = "four atoms, connected longitudinally";
 957     int               nAtom         = 4;
 958     std::vector<real> masses        = {0.5, 1.0/3.0, 0.25, 1.0};
 959     std::vector<int>  constraints   = {0, 0, 1, 1, 1, 2, 2, 2, 3};
 960     std::vector<real> constraintsR0 = {2.0, 1.0, 1.0};
 961
 962
 963     std::vector<RVec> x = {{  2.50, -3.10, 15.70 },
 964                            {  0.51, -3.02, 15.55 },
 965                            { -0.50, -3.00, 15.20 },
 966                            { -1.51, -2.95, 15.05 }};
 967
 968     std::vector<RVec> xPrime = {{  2.50, -3.10, 15.70 },
 969                                 {  0.51, -3.02, 15.55 },
 970                                 { -0.50, -3.00, 15.20 },
 971                                 { -1.51, -2.95, 15.05 }};
 972
 973     std::vector<RVec> v = {{ 0.0, 0.0,  2.0 },
 974                            { 0.0, 0.0,  3.0 },
 975                            { 0.0, 0.0, -4.0 },
 976                            { 0.0, 0.0, -1.0 }};
 977
 978     tensor            virialScaledRef = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
 979
 980     real              shakeTolerance         = 0.0001;
 981     gmx_bool          shakeUseSOR            = false;
 982
 983     int               lincsNIter             = 4;
 984     int               lincsNProjOrder        = 8;
 985     real              lincsWarnAngle         = 30.0;
 986
 987     std::unique_ptr<ConstraintsTestData>   testData = std::make_unique<ConstraintsTestData>
 988             (title, nAtom, masses,
 989             constraints, constraintsR0,
 990             false, virialScaledRef,
 991             false, 0,
 992             real(0.0), real(0.001),
 993             x, xPrime, v,
 994             shakeTolerance, shakeUseSOR,
 995             lincsNIter, lincsNProjOrder, lincsWarnAngle);
 996
 997     std::string pbcName;
 998     std::string algorithmName;
 999     std::tie(pbcName, algorithmName) = GetParam();
1000     t_pbc                                   pbc       = pbcs_.at(pbcName);
1001
1002     // Apply constraints
1003     algorithms_.at(algorithmName)(testData.get(), pbc);
1004
1005     checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
1006     checkConstrainsDirection(*testData, pbc);
1007     checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
1008     checkCOMVelocity(absoluteTolerance(0.0001), *testData);
1009
1010 }
1011
1012 TEST_P(ConstraintsTest, TriangleOfConstraints){
1013
1014     std::string       title         = "basic triangle (tree atoms, connected to each other)";
1015     int               nAtom         = 3;
1016     std::vector<real> masses        = {1.0, 1.0, 1.0};
1017     std::vector<int>  constraints   = {0, 0, 1, 2, 0, 2, 1, 1, 2};
1018     std::vector<real> constraintsR0 = {0.1, 0.1, 0.1};
1019
1020     real              oneTenthOverSqrtTwo    = 0.1_real / std::sqrt(2.0_real);
1021
1022     std::vector<RVec> x = {{ oneTenthOverSqrtTwo, 0.0, 0.0 },
1023                            { 0.0, oneTenthOverSqrtTwo, 0.0 },
1024                            { 0.0, 0.0, oneTenthOverSqrtTwo }};
1025
1026     std::vector<RVec> xPrime = {{  0.09, -0.02,  0.01 },
1027                                 { -0.02,  0.10, -0.02 },
1028                                 {  0.03, -0.01,  0.07 }};
1029
1030     std::vector<RVec> v = {{  1.0,  1.0,  1.0 },
1031                            { -2.0, -2.0, -2.0 },
1032                            {  1.0,  1.0,  1.0 }};
1033
1034     tensor            virialScaledRef = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
1035
1036     real              shakeTolerance         = 0.0001;
1037     gmx_bool          shakeUseSOR            = false;
1038
1039     int               lincsNIter             = 1;
1040     int               lincsNProjOrder        = 4;
1041     real              lincsWarnAngle         = 30.0;
1042
1043     std::unique_ptr<ConstraintsTestData>   testData = std::make_unique<ConstraintsTestData>
1044             (title, nAtom, masses,
1045             constraints, constraintsR0,
1046             false, virialScaledRef,
1047             false, 0,
1048             real(0.0), real(0.001),
1049             x, xPrime, v,
1050             shakeTolerance, shakeUseSOR,
1051             lincsNIter, lincsNProjOrder, lincsWarnAngle);
1052
1053     std::string pbcName;
1054     std::string algorithmName;
1055     std::tie(pbcName, algorithmName) = GetParam();
1056     t_pbc                                   pbc       = pbcs_.at(pbcName);
1057
1058     // Apply constraints
1059     algorithms_.at(algorithmName)(testData.get(), pbc);
1060
1061     checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
1062     checkConstrainsDirection(*testData, pbc);
1063     checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
1064     checkCOMVelocity(absoluteTolerance(0.0001), *testData);
1065
1066 }
1067
1068
1069 INSTANTIATE_TEST_CASE_P(WithParameters, ConstraintsTest,
1070                             ::testing::Combine(::testing::Values("PBCNone", "PBCXYZ"),
1071                                                    ::testing::Values("SHAKE", "LINCS")));
1072
1073 } // namespace test
1074 } // namespace gmx