From 5748002cd3b6c8f473e2abd51b759a30d1e6e084 Mon Sep 17 00:00:00 2001
From: Carlo Camilloni <carlo.camilloni@gmail.com>
Date: Mon, 14 Jul 2014 09:39:50 +0200
Subject: [PATCH] Movingrestraint: add kappa as a component
this should close #99
@davidebr can you double check?
---
src/bias/MovingRestraint.cpp | 20 ++++++++++++++++----
1 file changed, 16 insertions(+), 4 deletions(-)
diff --git a/src/bias/MovingRestraint.cpp b/src/bias/MovingRestraint.cpp
index 884815b57..1f3b8f355 100644
--- a/src/bias/MovingRestraint.cpp
+++ b/src/bias/MovingRestraint.cpp
@@ -121,12 +121,18 @@ PLUMED_REGISTER_ACTION(MovingRestraint,"MOVINGRESTRAINT")
void MovingRestraint::registerKeywords( Keywords& keys ){
Bias::registerKeywords(keys);
keys.use("ARG");
- keys.add("compulsory","VERSE","B","Tells plumed whether the restraint is only acting for CV larger (U) or smaller (L) than the restraint or whether it is acting on both sides (B)");
- keys.add("numbered","STEP","This keyword appears multiple times as STEPx with x=0,1,2,...,n. Each value given represents the MD step at which the restraint parameters take the values KAPPAx and ATx.");
+ keys.add("compulsory","VERSE","B","Tells plumed whether the restraint is only acting for CV larger (U) or smaller (L) than "
+ "the restraint or whether it is acting on both sides (B)");
+ keys.add("numbered","STEP","This keyword appears multiple times as STEPx with x=0,1,2,...,n. Each value given represents "
+ "the MD step at which the restraint parameters take the values KAPPAx and ATx.");
keys.reset_style("STEP","compulsory");
- keys.add("numbered","AT","ATx is equal to the position of the restraint at time STEPx. For intermediate times this parameter is linearly interpolated. If no ATx is specified for STEPx then the values of AT are kept constant during the interval of time between STEPx-1 and STEPx.");
+ keys.add("numbered","AT","ATx is equal to the position of the restraint at time STEPx. For intermediate times this parameter "
+ "is linearly interpolated. If no ATx is specified for STEPx then the values of AT are kept constant "
+ "during the interval of time between STEPx-1 and STEPx.");
keys.reset_style("AT","compulsory");
- keys.add("numbered","KAPPA","KAPPAx is equal to the value of the force constants at time STEPx. For intermediate times this parameter is linearly interpolated. If no KAPPAx is specified for STEPx then the values of KAPPAx are kept constant during the interval of time between STEPx-1 and STEPx.");
+ keys.add("numbered","KAPPA","KAPPAx is equal to the value of the force constants at time STEPx. For intermediate times this "
+ "parameter is linearly interpolated. If no KAPPAx is specified for STEPx then the values of KAPPAx "
+ "are kept constant during the interval of time between STEPx-1 and STEPx.");
keys.reset_style("KAPPA","compulsory");
componentsAreNotOptional(keys);
keys.addOutputComponent("bias","default","the instantaneous value of the bias potential");
@@ -139,6 +145,9 @@ void MovingRestraint::registerKeywords( Keywords& keys ){
"These quantities will named with the arguments of the bias followed by "
"the character string _work. These quantities tell the user how much work has "
"been done by the potential in dragging the system along the various colvar axis.");
+ keys.addOutputComponent("_kappa","default","one or multiple instances of this quantity will be refereceable elsewhere in the input file. "
+ "These quantities will named with the arguments of the bias followed by "
+ "the character string _kappa. These quantities tell the user the time dependent value of kappa.");
}
MovingRestraint::MovingRestraint(const ActionOptions&ao):
@@ -187,6 +196,8 @@ verse(getNumberOfArguments())
addComponent(comp); componentIsNotPeriodic(comp);
comp=getPntrToArgument(i)->getName()+"_work"; // each spring has its own work
addComponent(comp); componentIsNotPeriodic(comp);
+ comp=getPntrToArgument(i)->getName()+"_kappa"; // each spring has its own kappa
+ addComponent(comp); componentIsNotPeriodic(comp);
work.push_back(0.); // initialize the work value
}
@@ -231,6 +242,7 @@ void MovingRestraint::calculate(){
dpotdk[i]=0.5*cv*cv;
if(oldaa.size()==aa.size() && oldf.size()==f.size()) work[i]+=0.5*(oldf[i]+f[i])*(aa[i]-oldaa[i]) + 0.5*( dpotdk[i]+olddpotdk[i] )*(kk[i]-oldk[i]);
getPntrToComponent(getPntrToArgument(i)->getName()+"_work")->set(work[i]);
+ getPntrToComponent(getPntrToArgument(i)->getName()+"_kappa")->set(kk[i]);
ene+=0.5*k*cv*cv;
setOutputForce(i,f[i]);
totf2+=f[i]*f[i];
--
GitLab