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