diff --git a/src/multicolvar/AlphaBeta.cpp b/src/multicolvar/AlphaBeta.cpp index 8a7ecfdb61b9cdbd4f403e8bb957a92d3ac592fc..dc28588958f76f632f0569f529c64f0e218e67b9 100644 --- a/src/multicolvar/AlphaBeta.cpp +++ b/src/multicolvar/AlphaBeta.cpp @@ -31,15 +31,51 @@ using namespace std; namespace PLMD{ namespace multicolvar{ -//+PLUMEDOC MCOLVAR ALPHABETA +//+PLUMEDOC COLVAR ALPHABETA /* +Measures a distance including pbc between the instantaneous values of a set of torsional angles and set of reference values. + +This colvar calculates the following quantity. + +\f[ +s = \frac{1}{2} \sum_i \left[ 1 + \cos( \phi_i - \phi_i^{\textrm{Ref}} ) \right] +\f] + +where the \f$\phi_i\f$ values are the instantaneous values for the \ref TORSION angles of interest. +The \f$\phi_i^{\textrm{Ref}}\f$ values are the user-specified reference values for the torsional angles. + +\par Examples + +The following provides an example of the input for an alpha beta similarity. + +\verbatim +ALPHABETA ... +ATOMS1=168,170,172,188 REFERENCE1=3.14 +ATOMS2=170,172,188,190 REFERENCE2=3.14 +ATOMS3=188,190,192,230 REFERENCE3=3.14 +LABEL=ab +... ALPHABETA +PRINT ARG=ab FILE=colvar STRIDE=10 +\endverbatim + +Because all the reference values are the same we can calculate the same quantity using + +\verbatim +ALPHABETA ... +ATOMS1=168,170,172,188 REFERENCE=3.14 +ATOMS2=170,172,188,190 +ATOMS3=188,190,192,230 +LABEL=ab +... ALPHABETA +PRINT ARG=ab FILE=colvar STRIDE=10 +\endverbatim */ //+ENDPLUMEDOC class AlphaBeta : public MultiColvar { private: - double target; + std::vector<double> target; public: static void registerKeywords( Keywords& keys ); AlphaBeta(const ActionOptions&); @@ -53,8 +89,9 @@ PLUMED_REGISTER_ACTION(AlphaBeta,"ALPHABETA") void AlphaBeta::registerKeywords( Keywords& keys ){ MultiColvar::registerKeywords( keys ); keys.use("ATOMS"); - keys.add("optional","REFERENCE","a single reference value for all the dihedrals"); - //keys.add("optional","REFERENCE1","specific reference value for each dihedral"); + keys.add("numbered","REFERENCE","the reference values for each of the torsional angles. If you use a single REFERENCE value the " + "same reference value is used for all torsions"); + keys.reset_style("REFERENCE","compulsory"); } AlphaBeta::AlphaBeta(const ActionOptions&ao): @@ -62,11 +99,30 @@ PLUMED_MULTICOLVAR_INIT(ao) { // Read in the atoms int natoms=4; readAtoms( natoms ); - parse("REFERENCE",target); + // Resize target + target.resize( taskList.fullSize() ); + + // Read in reference values + unsigned ntarget=0; + for(unsigned i=0;i<target.size();++i){ + if( !parseNumbered( "REFERENCE", i+1, target[i] ) ) break; + ntarget++; + } + if( ntarget==0 ){ + parse("REFERENCE",target[0]); + for(unsigned i=1;i<target.size();++i) target[i]=target[0]; + } else if( ntarget!=target.size() ){ + error("found wrong number of REFERENCE values"); + } + // And setup the ActionWithVessel - std::string fake_input; - addVessel( "SUM", fake_input, -1 ); // -1 here means that this value will be named getLabel() readVesselKeywords(); + if( getNumberOfVessels()==0 ){ + std::string fake_input; + addVessel( "SUM", fake_input, -1 ); // -1 here means that this value will be named getLabel() + readVesselKeywords(); // This makes sure resizing is done + } + // And check everything has been read in correctly checkRead(); } @@ -80,8 +136,8 @@ double AlphaBeta::compute( const unsigned& j ){ Vector dd0,dd1,dd2; PLMD::Torsion t; double value = t.compute(d0,d1,d2,dd0,dd1,dd2); - double svalue = -0.5*sin(value-target); - double cvalue = 1.+cos(value-target); + double svalue = -0.5*sin(value-target[current]); + double cvalue = 1.+cos(value-target[current]); dd0 *= svalue; dd1 *= svalue; @@ -99,11 +155,9 @@ double AlphaBeta::compute( const unsigned& j ){ } Vector AlphaBeta::getCentralAtom(){ - addCentralAtomDerivatives( 0, 0.25*Tensor::identity() ); - addCentralAtomDerivatives( 1, 0.25*Tensor::identity() ); - addCentralAtomDerivatives( 2, 0.25*Tensor::identity() ); - addCentralAtomDerivatives( 3, 0.25*Tensor::identity() ); - return 0.25*( getPosition(0) + getPosition(1) + getPosition(2) + getPosition(3) ); + addCentralAtomDerivatives( 1, 0.5*Tensor::identity() ); + addCentralAtomDerivatives( 2, 0.5*Tensor::identity() ); + return 0.5*( getPosition(1) + getPosition(2) ); } }