Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-2023 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 "BridgeVessel.h"
      23             : #include "tools/Matrix.h"
      24             : #include "core/ActionWithArguments.h"
      25             : #include "StoreDataVessel.h"
      26             : 
      27             : namespace PLMD {
      28             : namespace vesselbase {
      29             : 
      30          46 : BridgeVessel::BridgeVessel( const VesselOptions& da ):
      31             :   Vessel(da),
      32          46 :   inum(0),
      33             : // in_normal_calculate(false)
      34          46 :   myOutputAction(NULL),
      35          46 :   myOutputValues(NULL),
      36          46 :   my_tmp_val(0,0)
      37             : {
      38          46 : }
      39             : 
      40         513 : void BridgeVessel::resize() {
      41         513 :   if( myOutputAction->checkNumericalDerivatives() ) {
      42          10 :     mynumerical_values.resize( getAction()->getNumberOfDerivatives()*myOutputValues->getNumberOfComponents() );
      43          10 :     inum=0;
      44             :   }
      45             :   // This bit ensures that we can store data in a bridge function if needs be
      46             :   // Notice we need to resize der_list in the underlying action as this is called
      47             :   // from a bridge
      48         513 :   if( myOutputAction->mydata ) {
      49             :     unsigned dsize=(myOutputAction->mydata)->getSizeOfDerivativeList();
      50          12 :     if( getAction()->der_list.size()!=dsize ) getAction()->der_list.resize( dsize );
      51             :   }
      52         513 :   unsigned tmp=0; resizeBuffer( myOutputAction->getSizeOfBuffer( tmp ) );
      53         513 : }
      54             : 
      55          46 : void BridgeVessel::setOutputAction( ActionWithVessel* myact ) {
      56          46 :   ActionWithValue* checkme=dynamic_cast<ActionWithValue*>( getAction() );
      57          46 :   plumed_massert( checkme, "vessel in bridge must inherit from ActionWithValue");
      58             : 
      59          46 :   myOutputAction=myact;
      60          46 :   myOutputValues=dynamic_cast<ActionWithValue*>( myact );
      61          46 :   plumed_massert( myOutputValues, "bridging vessel must inherit from ActionWithValue");
      62          46 : }
      63             : 
      64           0 : std::string BridgeVessel::description() {
      65           0 :   plumed_merror("I shouldn't end up here");
      66             : }
      67             : 
      68        3025 : void BridgeVessel::prepare() {
      69        3025 :   myOutputAction->doJobsRequiredBeforeTaskList();
      70        3025 : }
      71             : 
      72        3025 : void BridgeVessel::setBufferStart( unsigned& start ) {
      73        3025 :   myOutputAction->getSizeOfBuffer( start );
      74        3025 : }
      75             : 
      76      102743 : MultiValue& BridgeVessel::transformDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) {
      77      205022 :   if( outvals.getNumberOfValues()!=myOutputAction->getNumberOfQuantities() ||
      78      102279 :       outvals.getNumberOfDerivatives()!=myOutputAction->getNumberOfDerivatives() ) {
      79        2287 :     outvals.resize( myOutputAction->getNumberOfQuantities(), myOutputAction->getNumberOfDerivatives() );
      80             :   }
      81      102743 :   myOutputAction->transformBridgedDerivatives( current, invals, outvals );
      82      102743 :   return outvals;
      83             : }
      84             : 
      85      102743 : void BridgeVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
      86             :   // in_normal_calculate=true;
      87      102743 :   if( myvals.get(0)<myOutputAction->getTolerance() ) return;
      88       62508 :   myOutputAction->calculateAllVessels( current, myvals, myvals, buffer, der_list );
      89       62508 :   return;
      90             : }
      91             : 
      92        3025 : void BridgeVessel::finish( const std::vector<double>& buffer ) {
      93        3025 :   myOutputAction->finishComputations( buffer );
      94        3025 :   if( myOutputAction->checkNumericalDerivatives() ) {
      95        1930 :     if ( inum<mynumerical_values.size() ) {
      96        2160 :       for(int i=0; i<myOutputValues->getNumberOfComponents(); ++i) {
      97        1080 :         mynumerical_values[inum]=myOutputValues->getOutputQuantity(i);
      98        1080 :         inum++;
      99             :       }
     100             :       plumed_dbg_assert( inum<=mynumerical_values.size() );
     101             :     } else {
     102         850 :       plumed_assert( inum==mynumerical_values.size() );
     103             :     }
     104             :   }
     105        3025 : }
     106             : 
     107          65 : void BridgeVessel::completeNumericalDerivatives() {
     108          65 :   unsigned nextra = myOutputAction->getNumberOfDerivatives() - getAction()->getNumberOfDerivatives();
     109          65 :   Matrix<double> tmpder( myOutputValues->getNumberOfComponents(), nextra );
     110          65 :   ActionWithVessel* vval=dynamic_cast<ActionWithVessel*>( myOutputAction );
     111         785 :   for(unsigned i=0; i<nextra; ++i) {
     112         720 :     vval->bridgeVariable=i; getAction()->calculate();
     113        1440 :     for(int j=0; j<myOutputValues->getNumberOfComponents(); ++j) tmpder(j,i) = myOutputValues->getOutputQuantity(j);
     114             :   }
     115          65 :   vval->bridgeVariable=nextra; getAction()->calculate();
     116          65 :   plumed_assert( inum==mynumerical_values.size() ); inum=0;  // Reset inum now that we have finished calling calculate
     117          65 :   std::vector<double> base( myOutputValues->getNumberOfComponents() );
     118         130 :   for(int j=0; j<myOutputValues->getNumberOfComponents(); ++j) base[j] = myOutputValues->getOutputQuantity(j);
     119             : 
     120             :   const double delta=std::sqrt(epsilon);
     121          65 :   ActionAtomistic* aa=dynamic_cast<ActionAtomistic*>( getAction() );
     122          65 :   unsigned nvals=myOutputValues->getNumberOfComponents();
     123         130 :   for(unsigned j=0; j<nvals; ++j) ( myOutputValues->copyOutput(j) )->clearDerivatives();
     124             : 
     125          65 :   if( aa ) {
     126          65 :     ActionWithArguments* aarg=dynamic_cast<ActionWithArguments*>( getAction() );
     127         130 :     plumed_assert( !aarg ); Tensor box=aa->getBox();
     128             :     unsigned natoms=aa->getNumberOfAtoms();
     129         130 :     for(unsigned j=0; j<nvals; ++j) {
     130          65 :       double ref=( myOutputValues->copyOutput(j) )->get();
     131          65 :       if( ( myOutputValues->copyOutput(j) )->getNumberOfDerivatives()>0 ) {
     132         560 :         for(unsigned i=0; i<3*natoms; ++i) {
     133         495 :           double d=( mynumerical_values[i*nvals+j] - ref)/delta;
     134         495 :           ( myOutputValues->copyOutput(j) )->addDerivative(i,d);
     135             :         }
     136          65 :         Tensor virial;
     137         845 :         for(int i=0; i<3; i++) for(int k=0; k<3; k++) {
     138         585 :             virial(i,k)=( mynumerical_values[ nvals*(3*natoms + 3*i + k) + j ]-ref)/delta;
     139             :           }
     140          65 :         virial=-matmul(box.transpose(),virial);
     141         845 :         for(int i=0; i<3; i++) for(int k=0; k<3; k++) ( myOutputValues->copyOutput(j) )->addDerivative(3*natoms+3*k+i,virial(k,i));
     142             :       }
     143             :     }
     144             :   } else {
     145           0 :     plumed_merror("not implemented or tested yet");
     146             : //      unsigned nder=myOutputAction->getNumberOfDerivatives();
     147             : //      for(unsigned j=0;j<nvals;++j){
     148             : //          double ref=( myOutputValues->copyOutput(j) )->get();
     149             : //              for(unsigned i=0;i<nder;++i){
     150             : //                  double d=( mynumerical_values[i*nvals+j] - ref)/delta;
     151             : //                  ( myOutputValues->copyOutput(j) )->addDerivative(i,d);
     152             : //              }
     153             : //          }
     154             : //      }
     155             :   }
     156             :   // Add the derivatives wrt to the local quantities we are working with
     157         130 :   for(unsigned j=0; j<nvals; ++j) {
     158             :     unsigned k=0;
     159         785 :     for(unsigned i=getAction()->getNumberOfDerivatives(); i<myOutputAction->getNumberOfDerivatives(); ++i) {
     160         720 :       ( myOutputValues->copyOutput(j) )->addDerivative( i, (tmpder(j,k)-base[j])/std::sqrt(epsilon) ); k++;
     161             :     }
     162             :   }
     163          65 : }
     164             : 
     165         795 : bool BridgeVessel::applyForce( std::vector<double>& outforces ) {
     166         795 :   bool hasforce=false; outforces.assign(outforces.size(),0.0);
     167         795 :   unsigned ndertot = myOutputAction->getNumberOfDerivatives();
     168         795 :   unsigned nextra = ndertot - getAction()->getNumberOfDerivatives();
     169         795 :   std::vector<double> forces( ndertot ), eforces( nextra, 0.0 );
     170        1730 :   for(unsigned i=0; i<myOutputAction->getNumberOfVessels(); ++i) {
     171         935 :     if( ( myOutputAction->getPntrToVessel(i) )->applyForce( forces ) ) {
     172             :       hasforce=true;
     173      260815 :       for(unsigned j=0; j<outforces.size(); ++j) outforces[j]+=forces[j];
     174          94 :       for(unsigned j=0; j<nextra; ++j) eforces[j]+=forces[ outforces.size()+j ];
     175             :     }
     176             :   }
     177         795 :   if(hasforce) myOutputAction->applyBridgeForces( eforces );
     178         795 :   return hasforce;
     179             : }
     180             : 
     181         795 : void BridgeVessel::copyTaskFlags() {
     182         795 :   myOutputAction->deactivateAllTasks();
     183       89536 :   for(unsigned i=0; i<getAction()->nactive_tasks; ++i) myOutputAction->taskFlags[ getAction()->indexOfTaskInFullList[i] ] = 1;
     184         795 :   myOutputAction->lockContributors();
     185         795 : }
     186             : 
     187       13042 : MultiValue& BridgeVessel::getTemporyMultiValue() {
     188       13042 :   return my_tmp_val;
     189             : }
     190             : 
     191             : }
     192             : }
     193             :