src/gromacs/applied_forces/qmmm/qmmmoptions.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, 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
  37  * Implements QMMMOptions class
  38  *
  39  * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
  40  * \author Christian Blau <blau@kth.se>
  41  * \ingroup module_applied_forces
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "qmmmoptions.h"
  46
  47 #include <map>
  48
  49 #include "gromacs/fileio/warninp.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/options/basicoptions.h"
  52 #include "gromacs/options/optionsection.h"
  53 #include "gromacs/selection/indexutil.h"
  54 #include "gromacs/topology/mtop_lookup.h"
  55 #include "gromacs/topology/mtop_util.h"
  56 #include "gromacs/topology/topology.h"
  57 #include "gromacs/utility/exceptions.h"
  58 #include "gromacs/utility/keyvaluetreebuilder.h"
  59 #include "gromacs/utility/keyvaluetreetransform.h"
  60 #include "gromacs/utility/logger.h"
  61 #include "gromacs/utility/mdmodulesnotifiers.h"
  62 #include "gromacs/utility/path.h"
  63 #include "gromacs/utility/strconvert.h"
  64 #include "gromacs/utility/stringutil.h"
  65 #include "gromacs/utility/textreader.h"
  66
  67 #include "qmmminputgenerator.h"
  68 #include "qmmmtopologypreprocessor.h"
  69
  70 namespace gmx
  71 {
  72
  73 namespace
  74 {
  75
  76 /*! \brief Helper to declare mdp transform rules.
  77  *
  78  * Enforces uniform mdp options that are always prepended with the correct
  79  * string for the QMMM mdp options.
  80  *
  81  * \tparam ToType type to be transformed to
  82  * \tparam TransformWithFunctionType type of transformation function to be used
  83  *
  84  * \param[in] rules KVT transformation rules
  85  * \param[in] transformationFunction the function to transform the flat kvt tree
  86  * \param[in] optionTag string tag that describes the mdp option, appended to the
  87  *                      default string for the QMMM simulation
  88  */
  89 template<class ToType, class TransformWithFunctionType>
  90 void QMMMMdpTransformFromString(IKeyValueTreeTransformRules* rules,
  91                                 TransformWithFunctionType    transformationFunction,
  92                                 const std::string&           optionTag)
  93 {
  94     rules->addRule()
  95             .from<std::string>("/" + c_qmmmCP2KModuleName + "-" + optionTag)
  96             .to<ToType>("/" + c_qmmmCP2KModuleName + "/" + optionTag)
  97             .transformWith(transformationFunction);
  98 }
  99
 100 /*! \brief Helper to declare mdp output.
 101  *
 102  * Enforces uniform mdp options output strings that are always prepended with the
 103  * correct string for the QMMM mdp options and are consistent with the
 104  * options name and transformation type.
 105  *
 106  * \tparam OptionType the type of the mdp option
 107  * \param[in] builder the KVT builder to generate the output
 108  * \param[in] option the mdp option
 109  * \param[in] optionTag string tag that describes the mdp option, appended to the
 110  *                      default string for the QMMM simulation
 111  */
 112 template<class OptionType>
 113 void addQMMMMdpOutputValue(KeyValueTreeObjectBuilder* builder, const OptionType& option, const std::string& optionTag)
 114 {
 115     builder->addValue<OptionType>(c_qmmmCP2KModuleName + "-" + optionTag, option);
 116 }
 117
 118 /*! \brief Helper to declare mdp output comments.
 119  *
 120  * Enforces uniform mdp options comment output strings that are always prepended
 121  * with the correct string for the QMMM mdp options and are consistent
 122  * with the options name and transformation type.
 123  *
 124  * \param[in] builder the KVT builder to generate the output
 125  * \param[in] comment on the mdp option
 126  * \param[in] optionTag string tag that describes the mdp option
 127  */
 128 void addQMMMMdpOutputValueComment(KeyValueTreeObjectBuilder* builder,
 129                                   const std::string&         comment,
 130                                   const std::string&         optionTag)
 131 {
 132     builder->addValue<std::string>("comment-" + c_qmmmCP2KModuleName + "-" + optionTag, comment);
 133 }
 134
 135 } // namespace
 136
 137 void QMMMOptions::initMdpTransform(IKeyValueTreeTransformRules* rules)
 138 {
 139     const auto& stringIdentityTransform = [](std::string s) { return s; };
 140     QMMMMdpTransformFromString<bool>(rules, &fromStdString<bool>, c_activeTag_);
 141     QMMMMdpTransformFromString<std::string>(rules, stringIdentityTransform, c_qmGroupTag_);
 142     QMMMMdpTransformFromString<std::string>(rules, stringIdentityTransform, c_qmMethodTag_);
 143     QMMMMdpTransformFromString<int>(rules, &fromStdString<int>, c_qmChargeTag_);
 144     QMMMMdpTransformFromString<int>(rules, &fromStdString<int>, c_qmMultTag_);
 145     QMMMMdpTransformFromString<std::string>(rules, stringIdentityTransform, c_qmUserInputFileNameTag_);
 146 }
 147
 148 void QMMMOptions::buildMdpOutput(KeyValueTreeObjectBuilder* builder) const
 149 {
 150
 151     addQMMMMdpOutputValueComment(builder, "", "empty-line");
 152
 153     // Active flag
 154     addQMMMMdpOutputValueComment(builder, "; QM/MM with CP2K", "module");
 155     addQMMMMdpOutputValue(builder, parameters_.active_, c_activeTag_);
 156
 157     if (parameters_.active_)
 158     {
 159         // Index group for QM atoms, default System
 160         addQMMMMdpOutputValueComment(builder, "; Index group with QM atoms", c_qmGroupTag_);
 161         addQMMMMdpOutputValue(builder, groupString_, c_qmGroupTag_);
 162
 163         // QM method (DFT functional), default PBE
 164         addQMMMMdpOutputValueComment(builder, "; DFT functional for QM calculations", c_qmMethodTag_);
 165         addQMMMMdpOutputValue<std::string>(
 166                 builder, c_qmmmQMMethodNames[parameters_.qmMethod_], c_qmMethodTag_);
 167
 168         // QM charge, default 0
 169         addQMMMMdpOutputValueComment(builder, "; QM charge", c_qmChargeTag_);
 170         addQMMMMdpOutputValue(builder, parameters_.qmCharge_, c_qmChargeTag_);
 171
 172         // QM mutiplicity, default 1
 173         addQMMMMdpOutputValueComment(builder, "; QM multiplicity", c_qmMultTag_);
 174         addQMMMMdpOutputValue(builder, parameters_.qmMultiplicity_, c_qmMultTag_);
 175
 176         // QM input filename, default empty (will be deduced from *.tpr name during mdrun)
 177         addQMMMMdpOutputValueComment(
 178                 builder, "; Names of CP2K files during simulation", c_qmUserInputFileNameTag_);
 179         addQMMMMdpOutputValue(builder, parameters_.qmFileNameBase_, c_qmUserInputFileNameTag_);
 180     }
 181 }
 182
 183 void QMMMOptions::initMdpOptions(IOptionsContainerWithSections* options)
 184 {
 185     auto section = options->addSection(OptionSection(c_qmmmCP2KModuleName.c_str()));
 186
 187     section.addOption(BooleanOption(c_activeTag_.c_str()).store(&parameters_.active_));
 188     section.addOption(StringOption(c_qmGroupTag_.c_str()).store(&groupString_));
 189     section.addOption(EnumOption<QMMMQMMethod>(c_qmMethodTag_.c_str())
 190                               .enumValue(c_qmmmQMMethodNames)
 191                               .store(&parameters_.qmMethod_));
 192     section.addOption(StringOption(c_qmUserInputFileNameTag_.c_str()).store(&parameters_.qmFileNameBase_));
 193     section.addOption(IntegerOption(c_qmChargeTag_.c_str()).store(&parameters_.qmCharge_));
 194     section.addOption(IntegerOption(c_qmMultTag_.c_str()).store(&parameters_.qmMultiplicity_));
 195 }
 196
 197 bool QMMMOptions::active() const
 198 {
 199     return parameters_.active_;
 200 }
 201
 202 const QMMMParameters& QMMMOptions::parameters()
 203 {
 204     return parameters_;
 205 }
 206
 207 void QMMMOptions::setLogger(const MDLogger& logger)
 208 {
 209     // Exit if QMMM module is not active
 210     if (!parameters_.active_)
 211     {
 212         return;
 213     }
 214
 215     logger_ = &logger;
 216 }
 217
 218 void QMMMOptions::setWarninp(warninp* wi)
 219 {
 220     // Exit if QMMM module is not active
 221     if (!parameters_.active_)
 222     {
 223         return;
 224     }
 225
 226     wi_ = wi;
 227 }
 228
 229 void QMMMOptions::appendLog(const std::string& msg)
 230 {
 231     if (logger_)
 232     {
 233         GMX_LOG(logger_->info).asParagraph().appendText(msg);
 234     }
 235 }
 236
 237 void QMMMOptions::appendWarning(const std::string& msg)
 238 {
 239     if (wi_)
 240     {
 241         warning(wi_, msg);
 242     }
 243 }
 244
 245 void QMMMOptions::setQMMMGroupIndices(const IndexGroupsAndNames& indexGroupsAndNames)
 246 {
 247     // Exit if QMMM module is not active
 248     if (!parameters_.active_)
 249     {
 250         return;
 251     }
 252
 253     // Create QM index
 254     parameters_.qmIndices_ = indexGroupsAndNames.indices(groupString_);
 255
 256     // Check that group is not empty
 257     if (parameters_.qmIndices_.empty())
 258     {
 259         GMX_THROW(InconsistentInputError(formatString(
 260                 "Group %s defining QM atoms should not be empty.", groupString_.c_str())));
 261     }
 262
 263     // Create temporary index for the whole System
 264     auto systemIndices = indexGroupsAndNames.indices(std::string("System"));
 265
 266     // Sort qmindices_ and sindices_
 267     std::sort(parameters_.qmIndices_.begin(), parameters_.qmIndices_.end());
 268     std::sort(systemIndices.begin(), systemIndices.end());
 269
 270     // Create MM index
 271     parameters_.mmIndices_.reserve(systemIndices.size());
 272
 273     // Position in qmindicies_
 274     size_t j = 0;
 275     // Now loop over sindices_ and write to mmindices_ only the atoms which does not belong to qmindices_
 276     for (size_t i = 0; i < systemIndices.size(); i++)
 277     {
 278         if (systemIndices[i] != parameters_.qmIndices_[j])
 279         {
 280             parameters_.mmIndices_.push_back(systemIndices[i]);
 281         }
 282         else
 283         {
 284             if (j < parameters_.qmIndices_.size() - 1)
 285             {
 286                 j++;
 287             }
 288         }
 289     }
 290 }
 291
 292 void QMMMOptions::processTprFilename(const MdRunInputFilename& tprFilename)
 293 {
 294     // Exit if QMMM module is not active
 295     if (!parameters_.active_)
 296     {
 297         return;
 298     }
 299
 300     // Provided name should not be empty
 301     GMX_RELEASE_ASSERT(!tprFilename.mdRunFilename_.empty(),
 302                        "Filename of the *.tpr simulation file is empty");
 303
 304     // Exit if parameters_.qmFileNameBase_ has been provided
 305     if (!parameters_.qmFileNameBase_.empty())
 306     {
 307         return;
 308     }
 309
 310     parameters_.qmFileNameBase_ =
 311             gmx::Path::stripExtension(gmx::Path::getFilename(tprFilename.mdRunFilename_)) + "_cp2k";
 312 }
 313
 314 void QMMMOptions::processExternalInputFile()
 315 {
 316     // First check if we could read qmExternalInputFileName_
 317     TextReader fInp(qmExternalInputFileName_);
 318
 319     // Then we need to build a map with all CP2K parameters found in file
 320     std::map<std::string, std::string> cp2kParams;
 321
 322     // Key-value pair of the parameters
 323     std::pair<std::string, std::string> kv;
 324
 325     // Loop over all lines in the file
 326     std::string line;
 327     while (fInp.readLine(&line))
 328     {
 329         // Split line into words
 330         auto words = splitString(line);
 331
 332         // If we have 2 or more words in the line then build key-value pair
 333         if (words.size() >= 2)
 334         {
 335             kv.first  = words[0];
 336             kv.second = words[1];
 337
 338             // Convert kv.first to the upper case
 339             kv.first = toUpperCase(kv.first);
 340
 341             // Put into the map
 342             cp2kParams.insert(kv);
 343         }
 344     }
 345     fInp.close();
 346
 347     // Check if CHARGE found in file
 348     if (cp2kParams.count("CHARGE") == 0)
 349     {
 350         GMX_THROW(InconsistentInputError(
 351                 formatString("Parameter CHARGE not found in the external CP2K input file %s",
 352                              qmExternalInputFileName_.c_str())));
 353     }
 354
 355     // Check if MULTIPLICITY found in file
 356     if (cp2kParams.count("MULTIPLICITY") == 0)
 357     {
 358         GMX_THROW(InconsistentInputError(
 359                 formatString("Parameter MULTIPLICITY not found in the external CP2K input file %s",
 360                              qmExternalInputFileName_.c_str())));
 361     }
 362
 363     // Check if COORD_FILE_NAME found in file
 364     if (cp2kParams.count("COORD_FILE_NAME") == 0)
 365     {
 366         GMX_THROW(InconsistentInputError(formatString(
 367                 "Parameter COORD_FILE_NAME not found in the external CP2K input file %s",
 368                 qmExternalInputFileName_.c_str())));
 369     }
 370
 371     // Check if RUN_TYPE in the file present and is ENERGY_FORCE
 372     if (cp2kParams.count("RUN_TYPE") == 0)
 373     {
 374         GMX_THROW(InconsistentInputError(
 375                 formatString("Parameter RUN_TYPE not found in the external CP2K input file %s",
 376                              qmExternalInputFileName_.c_str())));
 377     }
 378
 379     // Check if RUN_TYPE in the file is equal to ENERGY_FORCE
 380     if (toUpperCase(cp2kParams["RUN_TYPE"]) != "ENERGY_FORCE")
 381     {
 382         GMX_THROW(InconsistentInputError(formatString(
 383                 "Parameter RUN_TYPE should be ENERGY_FORCE in the external CP2K input file %s",
 384                 qmExternalInputFileName_.c_str())));
 385     }
 386
 387     // Set parameters_.qmCharge_ and qmMult_
 388     parameters_.qmCharge_       = fromStdString<int>(cp2kParams["CHARGE"]);
 389     parameters_.qmMultiplicity_ = fromStdString<int>(cp2kParams["MULTIPLICITY"]);
 390
 391     // Read the whole input as one string and replace COORD_FILE_NAME value with %s placeholder
 392     std::string str      = TextReader::readFileToString(qmExternalInputFileName_);
 393     parameters_.qmInput_ = replaceAllWords(str, cp2kParams["COORD_FILE_NAME"], "%s");
 394 }
 395
 396 void QMMMOptions::setQMExternalInputFile(const QMInputFileName& qmExternalInputFileName)
 397 {
 398     // Exit if QMMM module is not active
 399     if (!parameters_.active_)
 400     {
 401         return;
 402     }
 403
 404     if (parameters_.qmMethod_ != QMMMQMMethod::INPUT)
 405     {
 406         if (qmExternalInputFileName.hasQMInputFileName_)
 407         {
 408             // If parameters_.qmMethod_ != INPUT then user should not provide external input file
 409             GMX_THROW(InconsistentInputError(
 410                     "External CP2K input file has been provided with -qmi option, but "
 411                     + c_qmmmCP2KModuleName + "-" + c_qmMethodTag_ + " is not INPUT"));
 412         }
 413
 414         // Exit if we dont need to process external input file
 415         return;
 416     }
 417
 418     // Case where user should provide external input file with -qmi option
 419     if (parameters_.qmMethod_ == QMMMQMMethod::INPUT && !qmExternalInputFileName.hasQMInputFileName_)
 420     {
 421         GMX_THROW(InconsistentInputError(c_qmmmCP2KModuleName + "-" + c_qmMethodTag_
 422                                          + " = INPUT requested, but external CP2K "
 423                                            "input file is not provided with -qmi option"));
 424     }
 425
 426     // If external input is provided by the user then we should process it and save into the parameters_
 427     qmExternalInputFileName_ = qmExternalInputFileName.qmInputFileName_;
 428     processExternalInputFile();
 429 }
 430
 431 void QMMMOptions::processCoordinates(const CoordinatesAndBoxPreprocessed& coord)
 432 {
 433     // Exit if QMMM module is not active
 434     if (!parameters_.active_)
 435     {
 436         return;
 437     }
 438
 439     QMMMInputGenerator inpGen(
 440             parameters_, coord.pbc_, coord.box_, atomCharges_, coord.coordinates_.unpaddedConstArrayRef());
 441
 442     // Generate pdb file with point charges for CP2K
 443     parameters_.qmPdb_ = inpGen.generateCP2KPdb();
 444
 445     // In case parameters_.qmMethod_ != INPUT we should generate CP2K Input, QM box and translation
 446     if (parameters_.qmMethod_ != QMMMQMMethod::INPUT)
 447     {
 448         /* Check if some of the box vectors dimension lower that 1 nm.
 449          * For SCF stability box should be big enough.
 450          */
 451         matrix box;
 452         copy_mat(coord.box_, box);
 453         if (norm(box[0]) < 1.0 || norm(box[1]) < 1.0 || norm(box[2]) < 1.0)
 454         {
 455             GMX_THROW(InconsistentInputError(
 456                     "One of the box vectors is shorter than 1 nm.\n"
 457                     "For stable CP2K SCF convergence all simulation box vectors should be "
 458                     ">= 1 nm. Please consider to increase simulation box or provide custom CP2K "
 459                     "input using "
 460                     + c_qmmmCP2KModuleName + "-" + c_qmMethodTag_ + " = INPUT"));
 461         }
 462
 463         parameters_.qmInput_ = inpGen.generateCP2KInput();
 464         copy_mat(inpGen.qmBox(), parameters_.qmBox_);
 465         parameters_.qmTrans_ = inpGen.qmTrans();
 466     }
 467 }
 468
 469 void QMMMOptions::modifyQMMMTopology(gmx_mtop_t* mtop)
 470 {
 471     // Exit if QMMM module is not active
 472     if (!parameters_.active_)
 473     {
 474         return;
 475     }
 476
 477     // Process topology
 478     QMMMTopologyPreprocessor topPrep(parameters_.qmIndices_);
 479     topPrep.preprocess(mtop);
 480
 481     // Get atom numbers
 482     parameters_.atomNumbers_ = copyOf(topPrep.atomNumbers());
 483
 484     // Get atom point charges
 485     atomCharges_ = copyOf(topPrep.atomCharges());
 486
 487     // Get Link Frontier
 488     parameters_.link_ = copyOf(topPrep.linkFrontier());
 489
 490     // Get info about modifications
 491     QMMMTopologyInfo topInfo = topPrep.topInfo();
 492
 493     // Cast int qmCharge_ to real as further calculations use floating point
 494     real qmC = static_cast<real>(parameters_.qmCharge_);
 495
 496     // Print message to the log about performed modifications
 497     std::string msg = "\nQMMM Interface with CP2K is active, topology was modified!\n";
 498
 499     msg += formatString(
 500             "Number of QM atoms: %d\nNumber of MM atoms: %d\n", topInfo.numQMAtoms, topInfo.numMMAtoms);
 501
 502     msg += formatString("Total charge of the classical system (before modifications): %.5f\n",
 503                         topInfo.remainingMMCharge + topInfo.totalClassicalChargeOfQMAtoms);
 504
 505     msg += formatString("Classical charge removed from QM atoms: %.5f\n",
 506                         topInfo.totalClassicalChargeOfQMAtoms);
 507
 508     if (topInfo.numVirtualSitesModified > 0)
 509     {
 510         msg += formatString(
 511                 "Note: There are %d virtual sites found, which are built from QM atoms only. "
 512                 "Classical charges on them have been removed as well.\n",
 513                 topInfo.numVirtualSitesModified);
 514     }
 515
 516     msg += formatString("Total charge of QMMM system (after modifications): %.5f\n",
 517                         qmC + topInfo.remainingMMCharge);
 518
 519     if (topInfo.numBondsRemoved > 0)
 520     {
 521         msg += formatString("Bonds removed: %d\n", topInfo.numBondsRemoved);
 522     }
 523
 524     if (topInfo.numAnglesRemoved > 0)
 525     {
 526         msg += formatString("Angles removed: %d\n", topInfo.numAnglesRemoved);
 527     }
 528
 529     if (topInfo.numDihedralsRemoved > 0)
 530     {
 531         msg += formatString("Dihedrals removed: %d\n", topInfo.numDihedralsRemoved);
 532     }
 533
 534     if (topInfo.numSettleRemoved > 0)
 535     {
 536         msg += formatString("Settles removed: %d\n", topInfo.numSettleRemoved);
 537     }
 538
 539     if (topInfo.numConnBondsAdded > 0)
 540     {
 541         msg += formatString("F_CONNBONDS (type 5 bonds) added: %d\n", topInfo.numConnBondsAdded);
 542     }
 543
 544     if (topInfo.numLinkBonds > 0)
 545     {
 546         msg += formatString("QM-MM broken bonds found: %d\n", topInfo.numLinkBonds);
 547     }
 548
 549     appendLog(msg + "\n");
 550
 551     /* We should warn the user if there is inconsistence between removed classical charges
 552      * on QM atoms and total QM charge
 553      */
 554     if (std::abs(topInfo.totalClassicalChargeOfQMAtoms - qmC) > 1E-5)
 555     {
 556         msg = formatString(
 557                 "Total charge of your QMMM system differs from classical system! "
 558                 "Consider manually spreading %.5lf charge over MM atoms nearby to the QM "
 559                 "region\n",
 560                 topInfo.totalClassicalChargeOfQMAtoms - qmC);
 561         appendWarning(msg);
 562     }
 563
 564     // If there are many constrained bonds in QM system then we should also warn the user
 565     if (topInfo.numConstrainedBondsInQMSubsystem > 2)
 566     {
 567         msg = "Your QM subsystem has a lot of constrained bonds. They probably have been "
 568               "generated automatically. That could produce an artifacts in the simulation. "
 569               "Consider constraints = none in the mdp file.\n";
 570         appendWarning(msg);
 571     }
 572 }
 573
 574 void QMMMOptions::writeInternalParametersToKvt(KeyValueTreeObjectBuilder treeBuilder)
 575 {
 576     // Write QM atoms index
 577     auto GroupIndexAdder =
 578             treeBuilder.addUniformArray<std::int64_t>(c_qmmmCP2KModuleName + "-" + c_qmGroupTag_);
 579     for (const auto& indexValue : parameters_.qmIndices_)
 580     {
 581         GroupIndexAdder.addValue(indexValue);
 582     }
 583
 584     // Write MM atoms index
 585     GroupIndexAdder =
 586             treeBuilder.addUniformArray<std::int64_t>(c_qmmmCP2KModuleName + "-" + c_mmGroupTag_);
 587     for (const auto& indexValue : parameters_.mmIndices_)
 588     {
 589         GroupIndexAdder.addValue(indexValue);
 590     }
 591
 592     // Write atoms numbers
 593     GroupIndexAdder =
 594             treeBuilder.addUniformArray<std::int64_t>(c_qmmmCP2KModuleName + "-" + c_atomNumbersTag_);
 595     for (const auto& indexValue : parameters_.atomNumbers_)
 596     {
 597         GroupIndexAdder.addValue(indexValue);
 598     }
 599
 600     // Write link
 601     GroupIndexAdder = treeBuilder.addUniformArray<std::int64_t>(c_qmmmCP2KModuleName + "-" + c_qmLinkTag_);
 602     for (const auto& indexValue : parameters_.link_)
 603     {
 604         GroupIndexAdder.addValue(indexValue.qm);
 605     }
 606     GroupIndexAdder = treeBuilder.addUniformArray<std::int64_t>(c_qmmmCP2KModuleName + "-" + c_mmLinkTag_);
 607     for (const auto& indexValue : parameters_.link_)
 608     {
 609         GroupIndexAdder.addValue(indexValue.mm);
 610     }
 611
 612     // Write CP2K input file content
 613     treeBuilder.addValue<std::string>(c_qmmmCP2KModuleName + "-" + c_qmInputTag_, parameters_.qmInput_);
 614
 615     // Write CP2K pdb file content
 616     treeBuilder.addValue<std::string>(c_qmmmCP2KModuleName + "-" + c_qmPdbTag_, parameters_.qmPdb_);
 617
 618     // Write QM box matrix
 619     auto DoubleArrayAdder = treeBuilder.addUniformArray<double>(c_qmmmCP2KModuleName + "-" + c_qmBoxTag_);
 620     for (int i = 0; i < DIM; i++)
 621     {
 622         for (int j = 0; j < DIM; j++)
 623         {
 624             DoubleArrayAdder.addValue(static_cast<double>(parameters_.qmBox_[i][j]));
 625         }
 626     }
 627
 628     // Write QM Translation vector
 629     DoubleArrayAdder = treeBuilder.addUniformArray<double>(c_qmmmCP2KModuleName + "-" + c_qmTransTag_);
 630     for (int i = 0; i < DIM; i++)
 631     {
 632         DoubleArrayAdder.addValue(static_cast<double>(parameters_.qmTrans_[i]));
 633     }
 634 }
 635
 636 void QMMMOptions::readInternalParametersFromKvt(const KeyValueTreeObject& tree)
 637 {
 638     // Check if active
 639     if (!parameters_.active_)
 640     {
 641         return;
 642     }
 643
 644     // Try to read QM atoms index
 645     if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmGroupTag_))
 646     {
 647         GMX_THROW(InconsistentInputError(
 648                 "Cannot find QM atoms index vector required for QM/MM simulation.\nThis could be "
 649                 "caused by incompatible or corrupted tpr input file."));
 650     }
 651     auto kvtIndexArray = tree[c_qmmmCP2KModuleName + "-" + c_qmGroupTag_].asArray().values();
 652     parameters_.qmIndices_.resize(kvtIndexArray.size());
 653     std::transform(std::begin(kvtIndexArray),
 654                    std::end(kvtIndexArray),
 655                    std::begin(parameters_.qmIndices_),
 656                    [](const KeyValueTreeValue& val) { return val.cast<std::int64_t>(); });
 657
 658     // Try to read MM atoms index
 659     if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_mmGroupTag_))
 660     {
 661         GMX_THROW(InconsistentInputError(
 662                 "Cannot find MM atoms index vector required for QM/MM simulation.\nThis could be "
 663                 "caused by incompatible or corrupted tpr input file."));
 664     }
 665     kvtIndexArray = tree[c_qmmmCP2KModuleName + "-" + c_mmGroupTag_].asArray().values();
 666     parameters_.mmIndices_.resize(kvtIndexArray.size());
 667     std::transform(std::begin(kvtIndexArray),
 668                    std::end(kvtIndexArray),
 669                    std::begin(parameters_.mmIndices_),
 670                    [](const KeyValueTreeValue& val) { return val.cast<std::int64_t>(); });
 671
 672     // Try to read atoms numbers
 673     if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_atomNumbersTag_))
 674     {
 675         GMX_THROW(InconsistentInputError(
 676                 "Cannot find Atom Numbers vector required for QM/MM simulation.\nThis could be "
 677                 "caused by incompatible or corrupted tpr input file."));
 678     }
 679     kvtIndexArray = tree[c_qmmmCP2KModuleName + "-" + c_atomNumbersTag_].asArray().values();
 680     parameters_.atomNumbers_.resize(kvtIndexArray.size());
 681     std::transform(std::begin(kvtIndexArray),
 682                    std::end(kvtIndexArray),
 683                    std::begin(parameters_.atomNumbers_),
 684                    [](const KeyValueTreeValue& val) { return val.cast<std::int64_t>(); });
 685
 686     // Try to read Link Frontier (two separate vectors and then combine)
 687     std::vector<index> qmLink;
 688     std::vector<index> mmLink;
 689
 690     if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmLinkTag_))
 691     {
 692         GMX_THROW(InconsistentInputError(
 693                 "Cannot find QM Link Frontier vector required for QM/MM simulation.\nThis could be "
 694                 "caused by incompatible or corrupted tpr input file."));
 695     }
 696     kvtIndexArray = tree[c_qmmmCP2KModuleName + "-" + c_qmLinkTag_].asArray().values();
 697     qmLink.resize(kvtIndexArray.size());
 698     std::transform(std::begin(kvtIndexArray),
 699                    std::end(kvtIndexArray),
 700                    std::begin(qmLink),
 701                    [](const KeyValueTreeValue& val) { return val.cast<std::int64_t>(); });
 702
 703     if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_mmLinkTag_))
 704     {
 705         GMX_THROW(InconsistentInputError(
 706                 "Cannot find MM Link Frontier vector required for QM/MM simulation.\nThis could be "
 707                 "caused by incompatible or corrupted tpr input file."));
 708     }
 709     kvtIndexArray = tree[c_qmmmCP2KModuleName + "-" + c_mmLinkTag_].asArray().values();
 710     mmLink.resize(kvtIndexArray.size());
 711     std::transform(std::begin(kvtIndexArray),
 712                    std::end(kvtIndexArray),
 713                    std::begin(mmLink),
 714                    [](const KeyValueTreeValue& val) { return val.cast<std::int64_t>(); });
 715
 716     parameters_.link_.resize(qmLink.size());
 717     for (size_t i = 0; i < qmLink.size(); i++)
 718     {
 719         parameters_.link_[i].qm = qmLink[i];
 720         parameters_.link_[i].mm = mmLink[i];
 721     }
 722
 723     // Try to read CP2K input and pdb strings from *.tpr
 724     if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmInputTag_))
 725     {
 726         GMX_THROW(InconsistentInputError(
 727                 "Cannot find CP2K input string required for QM/MM simulation.\nThis could be "
 728                 "caused by incompatible or corrupted tpr input file."));
 729     }
 730     parameters_.qmInput_ = tree[c_qmmmCP2KModuleName + "-" + c_qmInputTag_].cast<std::string>();
 731
 732     if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmPdbTag_))
 733     {
 734         GMX_THROW(InconsistentInputError(
 735                 "Cannot find CP2K pdb string required for QM/MM simulation.\nThis could be "
 736                 "caused by incompatible or corrupted tpr input file."));
 737     }
 738     parameters_.qmPdb_ = tree[c_qmmmCP2KModuleName + "-" + c_qmPdbTag_].cast<std::string>();
 739
 740     // Try to read QM box
 741     if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmBoxTag_))
 742     {
 743         GMX_THROW(InconsistentInputError(
 744                 "Cannot find QM box matrix required for QM/MM simulation.\nThis could be "
 745                 "caused by incompatible or corrupted tpr input file."));
 746     }
 747     auto kvtDoubleArray = tree[c_qmmmCP2KModuleName + "-" + c_qmBoxTag_].asArray().values();
 748     for (int i = 0; i < DIM; i++)
 749     {
 750         for (int j = 0; j < DIM; j++)
 751         {
 752             parameters_.qmBox_[i][j] = static_cast<real>(kvtDoubleArray[i * 3 + j].cast<double>());
 753         }
 754     }
 755
 756     // Try to read QM translation vector
 757     if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmTransTag_))
 758     {
 759         GMX_THROW(InconsistentInputError(
 760                 "Cannot find QM subsystem centering information for QM/MM simulation.\nThis could "
 761                 "be caused by incompatible or corrupted tpr input file."));
 762     }
 763     kvtDoubleArray = tree[c_qmmmCP2KModuleName + "-" + c_qmTransTag_].asArray().values();
 764     for (int i = 0; i < DIM; i++)
 765     {
 766         parameters_.qmTrans_[i] = static_cast<real>(kvtDoubleArray[i].cast<double>());
 767     }
 768 }
 769
 770 } // namespace gmx