From: Dmitry Morozov Date: Wed, 29 Sep 2021 17:49:01 +0000 (+0000) Subject: Activate QMMM MDModule X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=34dadf240683e112bd305e765ddf08fb2d50bfdf;p=alexxy%2Fgromacs.git Activate QMMM MDModule --- diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst index 95b3f226e3..d88973e007 100644 --- a/docs/user-guide/mdp-options.rst +++ b/docs/user-guide/mdp-options.rst @@ -3312,6 +3312,52 @@ electron-microscopy experiments. (See the `reference manual`_ for details) z-axis by :math:`\theta` degress by using following input: :math:`(\cos \theta , -\sin \theta , 0 , \sin \theta , \cos \theta , 0 , 0 , 0 , 1)` . +QM/MM simulations with CP2K Interface +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +These options enable and control the calculation and application of additional +QM/MM forces that are computed by the CP2K package if it is linked into |Gromacs|. +For further details about QM/MM interface implementation follow :ref:`qmmm`. + +.. mdp:: qmmm-cp2k-active + + (false) Activate QM/MM simulations. Requires CP2K to be linked with |Gromacs| + +.. mdp:: qmmm-cp2k-qmgroup + + (System) Index group with atoms that are treated with QM. + +.. mdp:: qmmm-cp2k-qmmethod + + (PBE) Method used to describe the QM part of the system. + + .. mdp-value:: PBE + + DFT using PBE functional and DZVP-MOLOPT basis set. + + .. mdp-value:: BLYP + + DFT using BLYP functional and DZVP-MOLOPT basis set. + + .. mdp-value:: INPUT + + Provide an external input file for CP2K when running :ref:`gmx grompp` with the ``-qmi`` command-line option. + External input files are subject to the limitations that are described in :ref:`qmmm`. + +.. mdp:: qmmm-cp2k-qmcharge + + (0) Total charge of the QM part. + +.. mdp:: qmmm-cp2k-qmmultiplicity + + (1) Multiplicity or spin-state of QM part. Default value 1 means singlet state. + +.. mdp:: qmmm-cp2k-qmfilenames + + () Names of the CP2K files that will be generated during the simulation. + When using the default, empty, value the name of the simulation input file will be used + with an additional ``_cp2k`` suffix. + User defined thingies ^^^^^^^^^^^^^^^^^^^^^ diff --git a/src/gromacs/applied_forces/qmmm/qmmmforceprovider.cpp b/src/gromacs/applied_forces/qmmm/qmmmforceprovider.cpp index f685efdcf7..451aae269d 100644 --- a/src/gromacs/applied_forces/qmmm/qmmmforceprovider.cpp +++ b/src/gromacs/applied_forces/qmmm/qmmmforceprovider.cpp @@ -103,7 +103,7 @@ void QMMMForceProvider::appendLog(const std::string& msg) void QMMMForceProvider::initCP2KForceEnvironment(const t_commrec& cr) { // Check that we have filename either defined in KVT or deduced from *.tpr name - GMX_RELEASE_ASSERT(parameters_.qmFileNameBase_.empty(), + GMX_RELEASE_ASSERT(!parameters_.qmFileNameBase_.empty(), "Names of CP2K input/output files is not defined in QMMMParameters"); diff --git a/src/gromacs/applied_forces/qmmm/qmmmoptions.cpp b/src/gromacs/applied_forces/qmmm/qmmmoptions.cpp index 7c314f72a1..6f8810854d 100644 --- a/src/gromacs/applied_forces/qmmm/qmmmoptions.cpp +++ b/src/gromacs/applied_forces/qmmm/qmmmoptions.cpp @@ -298,7 +298,7 @@ void QMMMOptions::processTprFilename(const MdRunInputFilename& tprFilename) } // Provided name should not be empty - GMX_RELEASE_ASSERT(tprFilename.mdRunFilename_.empty(), + GMX_RELEASE_ASSERT(!tprFilename.mdRunFilename_.empty(), "Filename of the *.tpr simulation file is empty"); // Exit if parameters_.qmFileNameBase_ has been provided @@ -368,6 +368,22 @@ void QMMMOptions::processExternalInputFile() qmExternalInputFileName_.c_str()))); } + // Check if RUN_TYPE in the file present and is ENERGY_FORCE + if (cp2kParams.count("RUN_TYPE") == 0) + { + GMX_THROW(InconsistentInputError( + formatString("Parameter RUN_TYPE not found in the external CP2K input file %s", + qmExternalInputFileName_.c_str()))); + } + + // Check if RUN_TYPE in the file is equal to ENERGY_FORCE + if (toUpperCase(cp2kParams["RUN_TYPE"]) != "ENERGY_FORCE") + { + GMX_THROW(InconsistentInputError(formatString( + "Parameter RUN_TYPE should be ENERGY_FORCE in the external CP2K input file %s", + qmExternalInputFileName_.c_str()))); + } + // Set parameters_.qmCharge_ and qmMult_ parameters_.qmCharge_ = fromStdString(cp2kParams["CHARGE"]); parameters_.qmMultiplicity_ = fromStdString(cp2kParams["MULTIPLICITY"]); @@ -390,9 +406,9 @@ void QMMMOptions::setQMExternalInputFile(const QMInputFileName& qmExternalInputF if (qmExternalInputFileName.hasQMInputFileName_) { // If parameters_.qmMethod_ != INPUT then user should not provide external input file - GMX_THROW( - InconsistentInputError("External CP2K Input file has been provided with -qmi " - "option, but qmmm-qmmethod is not INPUT")); + GMX_THROW(InconsistentInputError( + "External CP2K input file has been provided with -qmi option, but " + + c_qmmmCP2KModuleName + "-" + c_qmMethodTag_ + " is not INPUT")); } // Exit if we dont need to process external input file @@ -402,9 +418,9 @@ void QMMMOptions::setQMExternalInputFile(const QMInputFileName& qmExternalInputF // Case where user should provide external input file with -qmi option if (parameters_.qmMethod_ == QMMMQMMethod::INPUT && !qmExternalInputFileName.hasQMInputFileName_) { - GMX_THROW( - InconsistentInputError("qmmm-qmmethod = INPUT requested, but external CP2K Input " - "file is not provided with -qmi option")); + GMX_THROW(InconsistentInputError(c_qmmmCP2KModuleName + "-" + c_qmMethodTag_ + + " = INPUT requested, but external CP2K " + "input file is not provided with -qmi option")); } // If external input is provided by the user then we should process it and save into the parameters_ @@ -440,7 +456,8 @@ void QMMMOptions::processCoordinates(const CoordinatesAndBoxPreprocessed& coord) "One of the box vectors is shorter than 1 nm.\n" "For stable CP2K SCF convergence all simulation box vectors should be " ">= 1 nm. Please consider to increase simulation box or provide custom CP2K " - "input using qmmm-qmmethod = INPUT")); + "input using " + + c_qmmmCP2KModuleName + "-" + c_qmMethodTag_ + " = INPUT")); } parameters_.qmInput_ = inpGen.generateCP2KInput(); @@ -529,7 +546,7 @@ void QMMMOptions::modifyQMMMTopology(gmx_mtop_t* mtop) msg += formatString("QM-MM broken bonds found: %d\n", topInfo.numLinkBonds); } - appendLog(msg); + appendLog(msg + "\n"); /* We should warn the user if there is inconsistence between removed classical charges * on QM atoms and total QM charge @@ -547,9 +564,9 @@ void QMMMOptions::modifyQMMMTopology(gmx_mtop_t* mtop) // If there are many constrained bonds in QM system then we should also warn the user if (topInfo.numConstrainedBondsInQMSubsystem > 2) { - msg += "Your QM subsystem has a lot of constrained bonds. They probably have been " - "generated automatically. That could produce an artifacts in the simulation. " - "Consider constraints = none in the mdp file.\n"; + msg = "Your QM subsystem has a lot of constrained bonds. They probably have been " + "generated automatically. That could produce an artifacts in the simulation. " + "Consider constraints = none in the mdp file.\n"; appendWarning(msg); } } diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml index 58d37a7ef6..c88c2aa792 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsDefineParametersWithValuesIncludingAssignment.xml @@ -298,5 +298,8 @@ electric-field-z = 0 0 0 0 ; Density guided simulation density-guided-simulation-active = false + +; QM/MM with CP2K +qmmm-cp2k-active = false diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml index 8dd0919803..7b6968ea7c 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricField.xml @@ -298,5 +298,8 @@ electric-field-z = 0 0 0 0 ; Density guided simulation density-guided-simulation-active = false + +; QM/MM with CP2K +qmmm-cp2k-active = false diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml index bbefcc9748..4afc56f511 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldOscillating.xml @@ -298,5 +298,8 @@ electric-field-z = 3.7 7.5 0 0 ; Density guided simulation density-guided-simulation-active = false + +; QM/MM with CP2K +qmmm-cp2k-active = false diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml index 8b21cdf5ef..4aac61bfa2 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsElectricFieldPulsed.xml @@ -298,5 +298,8 @@ electric-field-z = 0 0 0 0 ; Density guided simulation density-guided-simulation-active = false + +; QM/MM with CP2K +qmmm-cp2k-active = false diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml index 72c39c0c0e..2f3bebf7c1 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsEmptyLines.xml @@ -298,5 +298,8 @@ electric-field-z = 0 0 0 0 ; Density guided simulation density-guided-simulation-active = false + +; QM/MM with CP2K +qmmm-cp2k-active = false diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml index 72c39c0c0e..2f3bebf7c1 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsImplicitSolventNo.xml @@ -298,5 +298,8 @@ electric-field-z = 0 0 0 0 ; Density guided simulation density-guided-simulation-active = false + +; QM/MM with CP2K +qmmm-cp2k-active = false diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml index 72c39c0c0e..2f3bebf7c1 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsKeyWithoutValue.xml @@ -298,5 +298,8 @@ electric-field-z = 0 0 0 0 ; Density guided simulation density-guided-simulation-active = false + +; QM/MM with CP2K +qmmm-cp2k-active = false diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml index 40e229819f..54545d90e3 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_AcceptsMimic.xml @@ -298,5 +298,8 @@ electric-field-z = 0 0 0 0 ; Density guided simulation density-guided-simulation-active = false + +; QM/MM with CP2K +qmmm-cp2k-active = false diff --git a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml index ed1446672a..136d8617e7 100644 --- a/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml +++ b/src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_HandlesDifferentKindsOfMdpLines.xml @@ -298,5 +298,8 @@ electric-field-z = 0 0 0 0 ; Density guided simulation density-guided-simulation-active = false + +; QM/MM with CP2K +qmmm-cp2k-active = false diff --git a/src/gromacs/mdrun/mdmodules.cpp b/src/gromacs/mdrun/mdmodules.cpp index 9022b26e80..ac1c57bd02 100644 --- a/src/gromacs/mdrun/mdmodules.cpp +++ b/src/gromacs/mdrun/mdmodules.cpp @@ -40,6 +40,7 @@ #include "gromacs/applied_forces/electricfield.h" #include "gromacs/applied_forces/densityfitting/densityfitting.h" +#include "gromacs/applied_forces/qmmm/qmmm.h" #include "gromacs/imd/imd.h" #include "gromacs/mdtypes/iforceprovider.h" #include "gromacs/mdtypes/imdmodule.h" @@ -66,6 +67,7 @@ public: densityFitting_(DensityFittingModuleInfo::create()), field_(createElectricFieldModule()), imd_(createInteractiveMolecularDynamicsModule()), + qmmm_(QMMMModuleInfo::create()), swapCoordinates_(createSwapCoordinatesModule()) { } @@ -76,6 +78,7 @@ public: auto appliedForcesOptions = options->addSection(OptionSection("applied-forces")); field_->mdpOptionProvider()->initMdpOptions(&appliedForcesOptions); densityFitting_->mdpOptionProvider()->initMdpOptions(&appliedForcesOptions); + qmmm_->mdpOptionProvider()->initMdpOptions(&appliedForcesOptions); // In future, other sections would also go here. } @@ -102,6 +105,7 @@ public: std::unique_ptr field_; std::unique_ptr forceProviders_; std::unique_ptr imd_; + std::unique_ptr qmmm_; std::unique_ptr swapCoordinates_; /*! \brief List of registered MDModules @@ -125,12 +129,14 @@ void MDModules::initMdpTransform(IKeyValueTreeTransformRules* rules) auto appliedForcesScope = rules->scopedTransform("/applied-forces"); impl_->field_->mdpOptionProvider()->initMdpTransform(appliedForcesScope.rules()); impl_->densityFitting_->mdpOptionProvider()->initMdpTransform(appliedForcesScope.rules()); + impl_->qmmm_->mdpOptionProvider()->initMdpTransform(appliedForcesScope.rules()); } void MDModules::buildMdpOutput(KeyValueTreeObjectBuilder* builder) { impl_->field_->mdpOptionProvider()->buildMdpOutput(builder); impl_->densityFitting_->mdpOptionProvider()->buildMdpOutput(builder); + impl_->qmmm_->mdpOptionProvider()->buildMdpOutput(builder); } void MDModules::assignOptionsToModules(const KeyValueTreeObject& params, IKeyValueTreeErrorHandler* errorHandler) @@ -167,6 +173,7 @@ ForceProviders* MDModules::initForceProviders() impl_->forceProviders_ = std::make_unique(); impl_->field_->initForceProviders(impl_->forceProviders_.get()); impl_->densityFitting_->initForceProviders(impl_->forceProviders_.get()); + impl_->qmmm_->initForceProviders(impl_->forceProviders_.get()); for (auto&& module : impl_->modules_) { module->initForceProviders(impl_->forceProviders_.get()); @@ -177,11 +184,13 @@ ForceProviders* MDModules::initForceProviders() void MDModules::subscribeToPreProcessingNotifications() { impl_->densityFitting_->subscribeToPreProcessingNotifications(&impl_->notifiers_); + impl_->qmmm_->subscribeToPreProcessingNotifications(&impl_->notifiers_); } void MDModules::subscribeToSimulationSetupNotifications() { impl_->densityFitting_->subscribeToSimulationSetupNotifications(&impl_->notifiers_); + impl_->qmmm_->subscribeToSimulationSetupNotifications(&impl_->notifiers_); } void MDModules::add(std::shared_ptr module)