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
^^^^^^^^^^^^^^^^^^^^^
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");
}
// 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
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<int>(cp2kParams["CHARGE"]);
parameters_.qmMultiplicity_ = fromStdString<int>(cp2kParams["MULTIPLICITY"]);
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
// 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_
"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();
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
// 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);
}
}
; Density guided simulation
density-guided-simulation-active = false
+
+; QM/MM with CP2K
+qmmm-cp2k-active = false
</String>
</ReferenceData>
; Density guided simulation
density-guided-simulation-active = false
+
+; QM/MM with CP2K
+qmmm-cp2k-active = false
</String>
</ReferenceData>
; Density guided simulation
density-guided-simulation-active = false
+
+; QM/MM with CP2K
+qmmm-cp2k-active = false
</String>
</ReferenceData>
; Density guided simulation
density-guided-simulation-active = false
+
+; QM/MM with CP2K
+qmmm-cp2k-active = false
</String>
</ReferenceData>
; Density guided simulation
density-guided-simulation-active = false
+
+; QM/MM with CP2K
+qmmm-cp2k-active = false
</String>
</ReferenceData>
; Density guided simulation
density-guided-simulation-active = false
+
+; QM/MM with CP2K
+qmmm-cp2k-active = false
</String>
</ReferenceData>
; Density guided simulation
density-guided-simulation-active = false
+
+; QM/MM with CP2K
+qmmm-cp2k-active = false
</String>
</ReferenceData>
; Density guided simulation
density-guided-simulation-active = false
+
+; QM/MM with CP2K
+qmmm-cp2k-active = false
</String>
</ReferenceData>
; Density guided simulation
density-guided-simulation-active = false
+
+; QM/MM with CP2K
+qmmm-cp2k-active = false
</String>
</ReferenceData>
#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"
densityFitting_(DensityFittingModuleInfo::create()),
field_(createElectricFieldModule()),
imd_(createInteractiveMolecularDynamicsModule()),
+ qmmm_(QMMMModuleInfo::create()),
swapCoordinates_(createSwapCoordinatesModule())
{
}
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.
}
std::unique_ptr<IMDModule> field_;
std::unique_ptr<ForceProviders> forceProviders_;
std::unique_ptr<IMDModule> imd_;
+ std::unique_ptr<IMDModule> qmmm_;
std::unique_ptr<IMDModule> swapCoordinates_;
/*! \brief List of registered MDModules
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)
impl_->forceProviders_ = std::make_unique<ForceProviders>();
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());
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<gmx::IMDModule> module)
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff