1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2    Copyright (c) 2014-2021 The plumed team
3    (see the PEOPLE file at the root of the distribution for a list of names)
4 
5    See http://www.plumed.org for more information.
6 
7    This file is part of plumed, version 2.
8 
9    plumed is free software: you can redistribute it and/or modify
10    it under the terms of the GNU Lesser General Public License as published by
11    the Free Software Foundation, either version 3 of the License, or
12    (at your option) any later version.
13 
14    plumed 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
17    GNU Lesser General Public License for more details.
18 
19    You should have received a copy of the GNU Lesser General Public License
20    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
21 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 #include "BridgedMultiColvarFunction.h"
23 #include "core/PlumedMain.h"
24 #include "core/ActionSet.h"
25 #include "CatomPack.h"
26 
27 namespace PLMD {
28 namespace multicolvar {
29 
registerKeywords(Keywords & keys)30 void BridgedMultiColvarFunction::registerKeywords( Keywords& keys ) {
31   MultiColvarBase::registerKeywords( keys );
32   keys.add("compulsory","DATA","The multicolvar that calculates the set of base quantities that we are interested in");
33 }
34 
BridgedMultiColvarFunction(const ActionOptions & ao)35 BridgedMultiColvarFunction::BridgedMultiColvarFunction(const ActionOptions&ao):
36   Action(ao),
37   MultiColvarBase(ao)
38 {
39   std::string mlab; parse("DATA",mlab);
40   mycolv = plumed.getActionSet().selectWithLabel<MultiColvarBase*>(mlab);
41   if(!mycolv) error("action labeled " + mlab + " does not exist or is not a multicolvar");
42 
43   // When using numerical derivatives here we must use numerical derivatives
44   // in base multicolvar
45   if( checkNumericalDerivatives() ) mycolv->useNumericalDerivatives();
46 
47   myBridgeVessel = mycolv->addBridgingVessel( this ); addDependency(mycolv);
48   weightHasDerivatives=true; usespecies=mycolv->usespecies;
49   // Number of tasks is the same as the number in the underlying MultiColvar
50   for(unsigned i=0; i<mycolv->getFullNumberOfTasks(); ++i) addTaskToList( mycolv->getTaskCode(i) );
51 }
52 
turnOnDerivatives()53 void BridgedMultiColvarFunction::turnOnDerivatives() {
54   BridgedMultiColvarFunction* check = dynamic_cast<BridgedMultiColvarFunction*>( mycolv );
55   if( check ) {
56     if( check->getNumberOfAtoms()>0 ) error("cannot calculate required derivatives of this quantity");
57   }
58   MultiColvarBase::turnOnDerivatives();
59 }
60 
transformBridgedDerivatives(const unsigned & current,MultiValue & invals,MultiValue & outvals) const61 void BridgedMultiColvarFunction::transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) const {
62   completeTask( current, invals, outvals );
63 
64   // Now update the outvals derivatives lists
65   if( derivativesAreRequired() ) {
66     outvals.emptyActiveMembers();
67     if( mycolv->isDensity() ) {
68       for(unsigned j=0; j<3; ++j) outvals.putIndexInActiveArray( 3*current+j );
69       for(unsigned j=invals.getNumberOfDerivatives()-9; j<invals.getNumberOfDerivatives(); ++j) outvals.putIndexInActiveArray(j);
70     } else {
71       for(unsigned j=0; j<invals.getNumberActive(); ++j) outvals.putIndexInActiveArray( invals.getActiveIndex(j) );
72     }
73     for(unsigned j=invals.getNumberOfDerivatives(); j<outvals.getNumberOfDerivatives(); ++j) outvals.putIndexInActiveArray( j );
74     outvals.completeUpdate();
75   }
76 }
77 
performTask(const unsigned & taskIndex,const unsigned & current,MultiValue & myvals) const78 void BridgedMultiColvarFunction::performTask( const unsigned& taskIndex, const unsigned& current, MultiValue& myvals ) const {
79   // This allows us to speed up the code as we don't need to reallocate memory on every call of perform task
80   MultiValue& invals=myBridgeVessel->getTemporyMultiValue();
81   if( invals.getNumberOfValues()!=mycolv->getNumberOfQuantities() ||
82       invals.getNumberOfDerivatives()!=mycolv->getNumberOfDerivatives() ) {
83     invals.resize( mycolv->getNumberOfQuantities(), mycolv->getNumberOfDerivatives() );
84   }
85   invals.clearAll(); mycolv->performTask( taskIndex, current, invals );
86   transformBridgedDerivatives( taskIndex, invals, myvals );
87 }
88 
calculateNumericalDerivatives(ActionWithValue * a)89 void BridgedMultiColvarFunction::calculateNumericalDerivatives( ActionWithValue* a ) {
90   if(!a) {
91     a=dynamic_cast<ActionWithValue*>(this);
92     plumed_massert(a,"cannot compute numerical derivatives for an action without values");
93   }
94   if( myBridgeVessel ) {
95     myBridgeVessel->completeNumericalDerivatives();
96   } else {
97     error("numerical derivatives are not implemented");
98   }
99 }
100 
applyBridgeForces(const std::vector<double> & bb)101 void BridgedMultiColvarFunction::applyBridgeForces( const std::vector<double>& bb ) {
102   if( getNumberOfAtoms()==0 ) return ;
103 
104   std::vector<Vector>& f( modifyForces() );
105   for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
106     f[i][0]+=bb[3*i+0]; f[i][1]+=bb[3*i+1]; f[i][2]+=bb[3*i+2];
107   }
108   applyForces();
109 }
110 
isPeriodic()111 bool BridgedMultiColvarFunction::isPeriodic() {
112   return mycolv->isPeriodic();
113 }
114 
deactivate_task(const unsigned & taskno)115 void BridgedMultiColvarFunction::deactivate_task( const unsigned& taskno ) {
116   plumed_merror("This should never be called");
117 }
118 
getCentralAtomPack(const unsigned & basn,const unsigned & curr,CatomPack & mypack)119 void BridgedMultiColvarFunction::getCentralAtomPack( const unsigned& basn, const unsigned& curr, CatomPack& mypack ) {
120   return mycolv->getCentralAtomPack( basn, curr, mypack );
121 }
122 
123 }
124 }
125