2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * Implements QMMMOptions class
39 * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
40 * \author Christian Blau <blau@kth.se>
41 * \ingroup module_applied_forces
45 #include "qmmmoptions.h"
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"
67 #include "qmmminputgenerator.h"
68 #include "qmmmtopologypreprocessor.h"
76 /*! \brief Helper to declare mdp transform rules.
78 * Enforces uniform mdp options that are always prepended with the correct
79 * string for the QMMM mdp options.
81 * \tparam ToType type to be transformed to
82 * \tparam TransformWithFunctionType type of transformation function to be used
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
89 template<class ToType, class TransformWithFunctionType>
90 void QMMMMdpTransformFromString(IKeyValueTreeTransformRules* rules,
91 TransformWithFunctionType transformationFunction,
92 const std::string& optionTag)
95 .from<std::string>("/" + c_qmmmCP2KModuleName + "-" + optionTag)
96 .to<ToType>("/" + c_qmmmCP2KModuleName + "/" + optionTag)
97 .transformWith(transformationFunction);
100 /*! \brief Helper to declare mdp output.
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.
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
112 template<class OptionType>
113 void addQMMMMdpOutputValue(KeyValueTreeObjectBuilder* builder, const OptionType& option, const std::string& optionTag)
115 builder->addValue<OptionType>(c_qmmmCP2KModuleName + "-" + optionTag, option);
118 /*! \brief Helper to declare mdp output comments.
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.
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
128 void addQMMMMdpOutputValueComment(KeyValueTreeObjectBuilder* builder,
129 const std::string& comment,
130 const std::string& optionTag)
132 builder->addValue<std::string>("comment-" + c_qmmmCP2KModuleName + "-" + optionTag, comment);
137 void QMMMOptions::initMdpTransform(IKeyValueTreeTransformRules* rules)
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_);
148 void QMMMOptions::buildMdpOutput(KeyValueTreeObjectBuilder* builder) const
151 addQMMMMdpOutputValueComment(builder, "", "empty-line");
154 addQMMMMdpOutputValueComment(builder, "; QM/MM with CP2K", "module");
155 addQMMMMdpOutputValue(builder, parameters_.active_, c_activeTag_);
157 if (parameters_.active_)
159 // Index group for QM atoms, default System
160 addQMMMMdpOutputValueComment(builder, "; Index group with QM atoms", c_qmGroupTag_);
161 addQMMMMdpOutputValue(builder, groupString_, c_qmGroupTag_);
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_);
168 // QM charge, default 0
169 addQMMMMdpOutputValueComment(builder, "; QM charge", c_qmChargeTag_);
170 addQMMMMdpOutputValue(builder, parameters_.qmCharge_, c_qmChargeTag_);
172 // QM mutiplicity, default 1
173 addQMMMMdpOutputValueComment(builder, "; QM multiplicity", c_qmMultTag_);
174 addQMMMMdpOutputValue(builder, parameters_.qmMultiplicity_, c_qmMultTag_);
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_);
183 void QMMMOptions::initMdpOptions(IOptionsContainerWithSections* options)
185 auto section = options->addSection(OptionSection(c_qmmmCP2KModuleName.c_str()));
187 section.addOption(BooleanOption(c_activeTag_.c_str()).store(¶meters_.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(¶meters_.qmMethod_));
192 section.addOption(StringOption(c_qmUserInputFileNameTag_.c_str()).store(¶meters_.qmFileNameBase_));
193 section.addOption(IntegerOption(c_qmChargeTag_.c_str()).store(¶meters_.qmCharge_));
194 section.addOption(IntegerOption(c_qmMultTag_.c_str()).store(¶meters_.qmMultiplicity_));
197 bool QMMMOptions::active() const
199 return parameters_.active_;
202 const QMMMParameters& QMMMOptions::parameters()
207 void QMMMOptions::setLogger(const MDLogger& logger)
209 // Exit if QMMM module is not active
210 if (!parameters_.active_)
218 void QMMMOptions::setWarninp(warninp* wi)
220 // Exit if QMMM module is not active
221 if (!parameters_.active_)
229 void QMMMOptions::appendLog(const std::string& msg)
233 GMX_LOG(logger_->info).asParagraph().appendText(msg);
237 void QMMMOptions::appendWarning(const std::string& msg)
245 void QMMMOptions::setQMMMGroupIndices(const IndexGroupsAndNames& indexGroupsAndNames)
247 // Exit if QMMM module is not active
248 if (!parameters_.active_)
254 parameters_.qmIndices_ = indexGroupsAndNames.indices(groupString_);
256 // Check that group is not empty
257 if (parameters_.qmIndices_.empty())
259 GMX_THROW(InconsistentInputError(formatString(
260 "Group %s defining QM atoms should not be empty.", groupString_.c_str())));
263 // Create temporary index for the whole System
264 auto systemIndices = indexGroupsAndNames.indices(std::string("System"));
266 // Sort qmindices_ and sindices_
267 std::sort(parameters_.qmIndices_.begin(), parameters_.qmIndices_.end());
268 std::sort(systemIndices.begin(), systemIndices.end());
271 parameters_.mmIndices_.reserve(systemIndices.size());
273 // Position in qmindicies_
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++)
278 if (systemIndices[i] != parameters_.qmIndices_[j])
280 parameters_.mmIndices_.push_back(systemIndices[i]);
284 if (j < parameters_.qmIndices_.size() - 1)
292 void QMMMOptions::processTprFilename(const MdRunInputFilename& tprFilename)
294 // Exit if QMMM module is not active
295 if (!parameters_.active_)
300 // Provided name should not be empty
301 GMX_RELEASE_ASSERT(tprFilename.mdRunFilename_.empty(),
302 "Filename of the *.tpr simulation file is empty");
304 // Exit if parameters_.qmFileNameBase_ has been provided
305 if (!parameters_.qmFileNameBase_.empty())
310 parameters_.qmFileNameBase_ =
311 gmx::Path::stripExtension(gmx::Path::getFilename(tprFilename.mdRunFilename_)) + "_cp2k";
314 void QMMMOptions::processExternalInputFile()
316 // First check if we could read qmExternalInputFileName_
317 TextReader fInp(qmExternalInputFileName_);
319 // Then we need to build a map with all CP2K parameters found in file
320 std::map<std::string, std::string> cp2kParams;
322 // Key-value pair of the parameters
323 std::pair<std::string, std::string> kv;
325 // Loop over all lines in the file
327 while (fInp.readLine(&line))
329 // Split line into words
330 auto words = splitString(line);
332 // If we have 2 or more words in the line then build key-value pair
333 if (words.size() >= 2)
336 kv.second = words[1];
338 // Convert kv.first to the upper case
339 kv.first = toUpperCase(kv.first);
342 cp2kParams.insert(kv);
347 // Check if CHARGE found in file
348 if (cp2kParams.count("CHARGE") == 0)
350 GMX_THROW(InconsistentInputError(
351 formatString("Parameter CHARGE not found in the external CP2K input file %s",
352 qmExternalInputFileName_.c_str())));
355 // Check if MULTIPLICITY found in file
356 if (cp2kParams.count("MULTIPLICITY") == 0)
358 GMX_THROW(InconsistentInputError(
359 formatString("Parameter MULTIPLICITY not found in the external CP2K input file %s",
360 qmExternalInputFileName_.c_str())));
363 // Check if COORD_FILE_NAME found in file
364 if (cp2kParams.count("COORD_FILE_NAME") == 0)
366 GMX_THROW(InconsistentInputError(formatString(
367 "Parameter COORD_FILE_NAME not found in the external CP2K input file %s",
368 qmExternalInputFileName_.c_str())));
371 // Set parameters_.qmCharge_ and qmMult_
372 parameters_.qmCharge_ = fromStdString<int>(cp2kParams["CHARGE"]);
373 parameters_.qmMultiplicity_ = fromStdString<int>(cp2kParams["MULTIPLICITY"]);
375 // Read the whole input as one string and replace COORD_FILE_NAME value with %s placeholder
376 std::string str = TextReader::readFileToString(qmExternalInputFileName_);
377 parameters_.qmInput_ = replaceAllWords(str, cp2kParams["COORD_FILE_NAME"], "%s");
380 void QMMMOptions::setQMExternalInputFile(const QMInputFileName& qmExternalInputFileName)
382 // Exit if QMMM module is not active
383 if (!parameters_.active_)
388 if (parameters_.qmMethod_ != QMMMQMMethod::INPUT)
390 if (qmExternalInputFileName.hasQMInputFileName_)
392 // If parameters_.qmMethod_ != INPUT then user should not provide external input file
394 InconsistentInputError("External CP2K Input file has been provided with -qmi "
395 "option, but qmmm-qmmethod is not INPUT"));
398 // Exit if we dont need to process external input file
402 // Case where user should provide external input file with -qmi option
403 if (parameters_.qmMethod_ == QMMMQMMethod::INPUT && !qmExternalInputFileName.hasQMInputFileName_)
406 InconsistentInputError("qmmm-qmmethod = INPUT requested, but external CP2K Input "
407 "file is not provided with -qmi option"));
410 // If external input is provided by the user then we should process it and save into the parameters_
411 qmExternalInputFileName_ = qmExternalInputFileName.qmInputFileName_;
412 processExternalInputFile();
415 void QMMMOptions::processCoordinates(const CoordinatesAndBoxPreprocessed& coord)
417 // Exit if QMMM module is not active
418 if (!parameters_.active_)
423 QMMMInputGenerator inpGen(
424 parameters_, coord.pbc_, coord.box_, atomCharges_, coord.coordinates_.unpaddedConstArrayRef());
426 // Generate pdb file with point charges for CP2K
427 parameters_.qmPdb_ = inpGen.generateCP2KPdb();
429 // In case parameters_.qmMethod_ != INPUT we should generate CP2K Input, QM box and translation
430 if (parameters_.qmMethod_ != QMMMQMMethod::INPUT)
432 /* Check if some of the box vectors dimension lower that 1 nm.
433 * For SCF stability box should be big enough.
436 copy_mat(coord.box_, box);
437 if (norm(box[0]) < 1.0 || norm(box[1]) < 1.0 || norm(box[2]) < 1.0)
439 GMX_THROW(InconsistentInputError(
440 "One of the box vectors is shorter than 1 nm.\n"
441 "For stable CP2K SCF convergence all simulation box vectors should be "
442 ">= 1 nm. Please consider to increase simulation box or provide custom CP2K "
443 "input using qmmm-qmmethod = INPUT"));
446 parameters_.qmInput_ = inpGen.generateCP2KInput();
447 copy_mat(inpGen.qmBox(), parameters_.qmBox_);
448 parameters_.qmTrans_ = inpGen.qmTrans();
452 void QMMMOptions::modifyQMMMTopology(gmx_mtop_t* mtop)
454 // Exit if QMMM module is not active
455 if (!parameters_.active_)
461 QMMMTopologyPreprocessor topPrep(parameters_.qmIndices_);
462 topPrep.preprocess(mtop);
465 parameters_.atomNumbers_ = copyOf(topPrep.atomNumbers());
467 // Get atom point charges
468 atomCharges_ = copyOf(topPrep.atomCharges());
471 parameters_.link_ = copyOf(topPrep.linkFrontier());
473 // Get info about modifications
474 QMMMTopologyInfo topInfo = topPrep.topInfo();
476 // Cast int qmCharge_ to real as further calculations use floating point
477 real qmC = static_cast<real>(parameters_.qmCharge_);
479 // Print message to the log about performed modifications
480 std::string msg = "\nQMMM Interface with CP2K is active, topology was modified!\n";
483 "Number of QM atoms: %d\nNumber of MM atoms: %d\n", topInfo.numQMAtoms, topInfo.numMMAtoms);
485 msg += formatString("Total charge of the classical system (before modifications): %.5f\n",
486 topInfo.remainingMMCharge + topInfo.totalClassicalChargeOfQMAtoms);
488 msg += formatString("Classical charge removed from QM atoms: %.5f\n",
489 topInfo.totalClassicalChargeOfQMAtoms);
491 if (topInfo.numVirtualSitesModified > 0)
494 "Note: There are %d virtual sites found, which are built from QM atoms only. "
495 "Classical charges on them have been removed as well.\n",
496 topInfo.numVirtualSitesModified);
499 msg += formatString("Total charge of QMMM system (after modifications): %.5f\n",
500 qmC + topInfo.remainingMMCharge);
502 if (topInfo.numBondsRemoved > 0)
504 msg += formatString("Bonds removed: %d\n", topInfo.numBondsRemoved);
507 if (topInfo.numAnglesRemoved > 0)
509 msg += formatString("Angles removed: %d\n", topInfo.numAnglesRemoved);
512 if (topInfo.numDihedralsRemoved > 0)
514 msg += formatString("Dihedrals removed: %d\n", topInfo.numDihedralsRemoved);
517 if (topInfo.numSettleRemoved > 0)
519 msg += formatString("Settles removed: %d\n", topInfo.numSettleRemoved);
522 if (topInfo.numConnBondsAdded > 0)
524 msg += formatString("F_CONNBONDS (type 5 bonds) added: %d\n", topInfo.numConnBondsAdded);
527 if (topInfo.numLinkBonds > 0)
529 msg += formatString("QM-MM broken bonds found: %d\n", topInfo.numLinkBonds);
534 /* We should warn the user if there is inconsistence between removed classical charges
535 * on QM atoms and total QM charge
537 if (std::abs(topInfo.totalClassicalChargeOfQMAtoms - qmC) > 1E-5)
540 "Total charge of your QMMM system differs from classical system! "
541 "Consider manually spreading %.5lf charge over MM atoms nearby to the QM "
543 topInfo.totalClassicalChargeOfQMAtoms - qmC);
547 // If there are many constrained bonds in QM system then we should also warn the user
548 if (topInfo.numConstrainedBondsInQMSubsystem > 2)
550 msg += "Your QM subsystem has a lot of constrained bonds. They probably have been "
551 "generated automatically. That could produce an artifacts in the simulation. "
552 "Consider constraints = none in the mdp file.\n";
557 void QMMMOptions::writeInternalParametersToKvt(KeyValueTreeObjectBuilder treeBuilder)
559 // Write QM atoms index
560 auto GroupIndexAdder =
561 treeBuilder.addUniformArray<std::int64_t>(c_qmmmCP2KModuleName + "-" + c_qmGroupTag_);
562 for (const auto& indexValue : parameters_.qmIndices_)
564 GroupIndexAdder.addValue(indexValue);
567 // Write MM atoms index
569 treeBuilder.addUniformArray<std::int64_t>(c_qmmmCP2KModuleName + "-" + c_mmGroupTag_);
570 for (const auto& indexValue : parameters_.mmIndices_)
572 GroupIndexAdder.addValue(indexValue);
575 // Write atoms numbers
577 treeBuilder.addUniformArray<std::int64_t>(c_qmmmCP2KModuleName + "-" + c_atomNumbersTag_);
578 for (const auto& indexValue : parameters_.atomNumbers_)
580 GroupIndexAdder.addValue(indexValue);
584 GroupIndexAdder = treeBuilder.addUniformArray<std::int64_t>(c_qmmmCP2KModuleName + "-" + c_qmLinkTag_);
585 for (const auto& indexValue : parameters_.link_)
587 GroupIndexAdder.addValue(indexValue.qm);
589 GroupIndexAdder = treeBuilder.addUniformArray<std::int64_t>(c_qmmmCP2KModuleName + "-" + c_mmLinkTag_);
590 for (const auto& indexValue : parameters_.link_)
592 GroupIndexAdder.addValue(indexValue.mm);
595 // Write CP2K input file content
596 treeBuilder.addValue<std::string>(c_qmmmCP2KModuleName + "-" + c_qmInputTag_, parameters_.qmInput_);
598 // Write CP2K pdb file content
599 treeBuilder.addValue<std::string>(c_qmmmCP2KModuleName + "-" + c_qmPdbTag_, parameters_.qmPdb_);
601 // Write QM box matrix
602 auto DoubleArrayAdder = treeBuilder.addUniformArray<double>(c_qmmmCP2KModuleName + "-" + c_qmBoxTag_);
603 for (int i = 0; i < DIM; i++)
605 for (int j = 0; j < DIM; j++)
607 DoubleArrayAdder.addValue(static_cast<double>(parameters_.qmBox_[i][j]));
611 // Write QM Translation vector
612 DoubleArrayAdder = treeBuilder.addUniformArray<double>(c_qmmmCP2KModuleName + "-" + c_qmTransTag_);
613 for (int i = 0; i < DIM; i++)
615 DoubleArrayAdder.addValue(static_cast<double>(parameters_.qmTrans_[i]));
619 void QMMMOptions::readInternalParametersFromKvt(const KeyValueTreeObject& tree)
622 if (!parameters_.active_)
627 // Try to read QM atoms index
628 if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmGroupTag_))
630 GMX_THROW(InconsistentInputError(
631 "Cannot find QM atoms index vector required for QM/MM simulation.\nThis could be "
632 "caused by incompatible or corrupted tpr input file."));
634 auto kvtIndexArray = tree[c_qmmmCP2KModuleName + "-" + c_qmGroupTag_].asArray().values();
635 parameters_.qmIndices_.resize(kvtIndexArray.size());
636 std::transform(std::begin(kvtIndexArray),
637 std::end(kvtIndexArray),
638 std::begin(parameters_.qmIndices_),
639 [](const KeyValueTreeValue& val) { return val.cast<std::int64_t>(); });
641 // Try to read MM atoms index
642 if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_mmGroupTag_))
644 GMX_THROW(InconsistentInputError(
645 "Cannot find MM atoms index vector required for QM/MM simulation.\nThis could be "
646 "caused by incompatible or corrupted tpr input file."));
648 kvtIndexArray = tree[c_qmmmCP2KModuleName + "-" + c_mmGroupTag_].asArray().values();
649 parameters_.mmIndices_.resize(kvtIndexArray.size());
650 std::transform(std::begin(kvtIndexArray),
651 std::end(kvtIndexArray),
652 std::begin(parameters_.mmIndices_),
653 [](const KeyValueTreeValue& val) { return val.cast<std::int64_t>(); });
655 // Try to read atoms numbers
656 if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_atomNumbersTag_))
658 GMX_THROW(InconsistentInputError(
659 "Cannot find Atom Numbers vector required for QM/MM simulation.\nThis could be "
660 "caused by incompatible or corrupted tpr input file."));
662 kvtIndexArray = tree[c_qmmmCP2KModuleName + "-" + c_atomNumbersTag_].asArray().values();
663 parameters_.atomNumbers_.resize(kvtIndexArray.size());
664 std::transform(std::begin(kvtIndexArray),
665 std::end(kvtIndexArray),
666 std::begin(parameters_.atomNumbers_),
667 [](const KeyValueTreeValue& val) { return val.cast<std::int64_t>(); });
669 // Try to read Link Frontier (two separate vectors and then combine)
670 std::vector<index> qmLink;
671 std::vector<index> mmLink;
673 if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmLinkTag_))
675 GMX_THROW(InconsistentInputError(
676 "Cannot find QM Link Frontier vector required for QM/MM simulation.\nThis could be "
677 "caused by incompatible or corrupted tpr input file."));
679 kvtIndexArray = tree[c_qmmmCP2KModuleName + "-" + c_qmLinkTag_].asArray().values();
680 qmLink.resize(kvtIndexArray.size());
681 std::transform(std::begin(kvtIndexArray),
682 std::end(kvtIndexArray),
684 [](const KeyValueTreeValue& val) { return val.cast<std::int64_t>(); });
686 if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_mmLinkTag_))
688 GMX_THROW(InconsistentInputError(
689 "Cannot find MM Link Frontier vector required for QM/MM simulation.\nThis could be "
690 "caused by incompatible or corrupted tpr input file."));
692 kvtIndexArray = tree[c_qmmmCP2KModuleName + "-" + c_mmLinkTag_].asArray().values();
693 mmLink.resize(kvtIndexArray.size());
694 std::transform(std::begin(kvtIndexArray),
695 std::end(kvtIndexArray),
697 [](const KeyValueTreeValue& val) { return val.cast<std::int64_t>(); });
699 parameters_.link_.resize(qmLink.size());
700 for (size_t i = 0; i < qmLink.size(); i++)
702 parameters_.link_[i].qm = qmLink[i];
703 parameters_.link_[i].mm = mmLink[i];
706 // Try to read CP2K input and pdb strings from *.tpr
707 if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmInputTag_))
709 GMX_THROW(InconsistentInputError(
710 "Cannot find CP2K input string required for QM/MM simulation.\nThis could be "
711 "caused by incompatible or corrupted tpr input file."));
713 parameters_.qmInput_ = tree[c_qmmmCP2KModuleName + "-" + c_qmInputTag_].cast<std::string>();
715 if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmPdbTag_))
717 GMX_THROW(InconsistentInputError(
718 "Cannot find CP2K pdb string required for QM/MM simulation.\nThis could be "
719 "caused by incompatible or corrupted tpr input file."));
721 parameters_.qmPdb_ = tree[c_qmmmCP2KModuleName + "-" + c_qmPdbTag_].cast<std::string>();
723 // Try to read QM box
724 if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmBoxTag_))
726 GMX_THROW(InconsistentInputError(
727 "Cannot find QM box matrix required for QM/MM simulation.\nThis could be "
728 "caused by incompatible or corrupted tpr input file."));
730 auto kvtDoubleArray = tree[c_qmmmCP2KModuleName + "-" + c_qmBoxTag_].asArray().values();
731 for (int i = 0; i < DIM; i++)
733 for (int j = 0; j < DIM; j++)
735 parameters_.qmBox_[i][j] = static_cast<real>(kvtDoubleArray[i * 3 + j].cast<double>());
739 // Try to read QM translation vector
740 if (!tree.keyExists(c_qmmmCP2KModuleName + "-" + c_qmTransTag_))
742 GMX_THROW(InconsistentInputError(
743 "Cannot find QM subsystem centering information for QM/MM simulation.\nThis could "
744 "be caused by incompatible or corrupted tpr input file."));
746 kvtDoubleArray = tree[c_qmmmCP2KModuleName + "-" + c_qmTransTag_].asArray().values();
747 for (int i = 0; i < DIM; i++)
749 parameters_.qmTrans_[i] = static_cast<real>(kvtDoubleArray[i].cast<double>());