Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-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 "Value.h"
      23             : #include "ActionWithValue.h"
      24             : #include "ActionAtomistic.h"
      25             : #include "ActionWithArguments.h"
      26             : #include "ActionWithVector.h"
      27             : #include "ActionWithVirtualAtom.h"
      28             : #include "tools/Exception.h"
      29             : #include "tools/OpenMP.h"
      30             : #include "tools/OFile.h"
      31             : #include "PlumedMain.h"
      32             : 
      33             : namespace PLMD {
      34             : 
      35         164 : Value::Value():
      36         164 :   action(NULL),
      37         164 :   value_set(false),
      38         164 :   hasForce(false),
      39         164 :   storedata(false),
      40             :   shape(std::vector<unsigned>()),
      41         164 :   hasDeriv(true),
      42         164 :   bufstart(0),
      43         164 :   streampos(0),
      44         164 :   matpos(0),
      45         164 :   ngrid_der(0),
      46         164 :   ncols(0),
      47         164 :   book_start(0),
      48         164 :   symmetric(false),
      49         164 :   periodicity(unset),
      50         164 :   min(0.0),
      51         164 :   max(0.0),
      52         164 :   max_minus_min(0.0),
      53         164 :   inv_max_minus_min(0.0),
      54         164 :   derivativeIsZeroWhenValueIsZero(false)
      55             : {
      56         164 :   data.resize(1); inputForce.resize(1);
      57         164 : }
      58             : 
      59          18 : Value::Value(const std::string& name):
      60          18 :   action(NULL),
      61          18 :   value_set(false),
      62          18 :   hasForce(false),
      63          18 :   name(name),
      64          18 :   storedata(false),
      65             :   shape(std::vector<unsigned>()),
      66          18 :   hasDeriv(true),
      67          18 :   bufstart(0),
      68          18 :   streampos(0),
      69          18 :   ngrid_der(0),
      70          18 :   matpos(0),
      71          18 :   ncols(0),
      72          18 :   book_start(0),
      73          18 :   symmetric(false),
      74          18 :   periodicity(unset),
      75          18 :   min(0.0),
      76          18 :   max(0.0),
      77          18 :   max_minus_min(0.0),
      78          18 :   inv_max_minus_min(0.0),
      79          18 :   derivativeIsZeroWhenValueIsZero(false)
      80             : {
      81          18 :   data.resize(1); inputForce.resize(1);
      82          18 :   data[0]=inputForce[0]=0;
      83          18 : }
      84             : 
      85      118722 : Value::Value(ActionWithValue* av, const std::string& name, const bool withderiv, const std::vector<unsigned>&ss):
      86      118722 :   action(av),
      87      118722 :   value_set(false),
      88      118722 :   hasForce(false),
      89      118722 :   name(name),
      90      118722 :   storedata(false),
      91      118722 :   hasDeriv(withderiv),
      92      118722 :   bufstart(0),
      93      118722 :   streampos(0),
      94      118722 :   ngrid_der(0),
      95      118722 :   matpos(0),
      96      118722 :   ncols(0),
      97      118722 :   book_start(0),
      98      118722 :   symmetric(false),
      99      118722 :   periodicity(unset),
     100      118722 :   min(0.0),
     101      118722 :   max(0.0),
     102      118722 :   max_minus_min(0.0),
     103      118722 :   inv_max_minus_min(0.0),
     104      118722 :   derivativeIsZeroWhenValueIsZero(false)
     105             : {
     106      118722 :   if( action ) {
     107      118266 :     if( action->getName()=="ACCUMULATE" || action->getName()=="COLLECT" ) valtype=average;
     108             :   }
     109      229748 :   if( action ) storedata=action->getName()=="PUT" || valtype==average;
     110      118722 :   if( ss.size() && withderiv ) storedata=true;
     111      118722 :   setShape( ss );
     112      118722 : }
     113             : 
     114         455 : void Value::setValType( const std::string& vtype ) {
     115         455 :   if( vtype=="normal" ) valtype=normal;
     116         455 :   else if( vtype=="constant" ) valtype=constant;
     117         455 :   else if( vtype=="average" ) valtype=average;
     118         455 :   else if( vtype=="calcFromAverage" ) valtype=calcFromAverage;
     119           0 :   else plumed_merror("invalid valtype " + vtype );
     120         455 : }
     121             : 
     122      131045 : void Value::setShape( const std::vector<unsigned>&ss ) {
     123      131045 :   std::size_t tot=1; shape.resize( ss.size() );
     124      161851 :   for(unsigned i=0; i<shape.size(); ++i) { tot = tot*ss[i]; shape[i]=ss[i]; }
     125             : 
     126      131045 :   if( shape.size()>0 && hasDeriv ) {
     127             :     // This is for grids
     128        2005 :     ngrid_der = shape.size();
     129        2005 :     if( action ) ngrid_der = action->getNumberOfDerivatives();
     130        2005 :     std::size_t ndata = tot*(1+ngrid_der);
     131        2005 :     data.resize( ndata ); inputForce.resize( tot );
     132      129040 :   } else if( shape.size()==0 ) {
     133             :     // This is for scalars
     134      106518 :     data.resize(1); inputForce.resize(1);
     135       22522 :   } else if( storedata && shape.size()<2 ) {
     136             :     // This is for vectors (matrices have special version because we have sparse storage)
     137       12419 :     data.resize( tot ); inputForce.resize( tot );
     138             :   }
     139      131045 : }
     140             : 
     141      109806 : void Value::setupPeriodicity() {
     142      109806 :   if( min==0 && max==0 ) {
     143      102039 :     periodicity=notperiodic;
     144             :   } else {
     145        7767 :     periodicity=periodic;
     146        7767 :     max_minus_min=max-min;
     147        7767 :     plumed_massert(max_minus_min>0, "your function has a very strange domain?");
     148        7767 :     inv_max_minus_min=1.0/max_minus_min;
     149             :   }
     150      109806 : }
     151             : 
     152     1685475 : bool Value::isPeriodic()const {
     153     1685475 :   plumed_massert(periodicity!=unset,"periodicity should be set");
     154     1685475 :   return periodicity==periodic;
     155             : }
     156             : 
     157      409249 : bool Value::applyForce(std::vector<double>& forces ) const {
     158      409249 :   if( !hasForce || valtype!=normal ) return false;
     159             :   plumed_dbg_massert( data.size()-1==forces.size()," forces array has wrong size" );
     160        9165 :   const unsigned N=data.size()-1;
     161    41672481 :   for(unsigned i=0; i<N; ++i) forces[i]=inputForce[0]*data[1+i];
     162             :   return true;
     163             : }
     164             : 
     165     1102380 : void Value::setNotPeriodic() {
     166     1102380 :   min=0; max=0; periodicity=notperiodic;
     167     1102380 : }
     168             : 
     169        7767 : void Value::setDomain(const std::string& pmin,const std::string& pmax) {
     170        7767 :   str_min=pmin;
     171        7767 :   if( !Tools::convertNoexcept(str_min,min) ) action->error("could not convert period string " + str_min + " to real");
     172        7767 :   str_max=pmax;
     173        7767 :   if( !Tools::convertNoexcept(str_max,max) ) action->error("could not convert period string " + str_max + " to read");
     174        7767 :   setupPeriodicity();
     175        7767 : }
     176             : 
     177       16759 : void Value::getDomain(std::string&minout,std::string&maxout) const {
     178       16759 :   plumed_massert(periodicity==periodic,"function should be periodic");
     179       16759 :   minout=str_min;
     180       16759 :   maxout=str_max;
     181       16759 : }
     182             : 
     183          34 : void Value::getDomain(double&minout,double&maxout) const {
     184          34 :   plumed_massert(periodicity==periodic,"function should be periodic");
     185          34 :   minout=min;
     186          34 :   maxout=max;
     187          34 : }
     188             : 
     189         190 : void Value::setGradients( ActionAtomistic* aa, unsigned& start ) {
     190             :   // Can't do gradients if we don't have derivatives
     191         190 :   if( !hasDeriv ) return;
     192         154 :   plumed_assert( shape.size()==0 );
     193        8362 :   for(unsigned j=0; j<aa->getNumberOfAtoms(); ++j) {
     194        8208 :     Vector der(data[1+start+3*j],data[1+start+3*j+1],data[1+start+3*j+2]);
     195        8208 :     aa->getGradient( j, der, gradients );
     196             :   }
     197         154 :   start += aa->getNumberOfAtoms();
     198             : }
     199             : 
     200         283 : void Value::passGradients( const double& der, std::map<AtomNumber,Vector>& g ) const {
     201        8921 :   for(const auto & p : gradients) { AtomNumber iatom=p.first; g[iatom] += p.second*der; }
     202         283 : }
     203             : 
     204         261 : double Value::projection(const Value& v1,const Value&v2) {
     205             :   double proj=0.0;
     206             :   const std::map<AtomNumber,Vector> & grad1(v1.gradients);
     207             :   const std::map<AtomNumber,Vector> & grad2(v2.gradients);
     208       16167 :   for(const auto & p1 : grad1) {
     209       15906 :     AtomNumber a=p1.first;
     210             :     const auto p2=grad2.find(a);
     211       15906 :     if(p2!=grad2.end()) {
     212        8224 :       proj+=dotProduct(p1.second,(*p2).second);
     213             :     }
     214             :   }
     215         261 :   return proj;
     216             : }
     217             : 
     218    18718290 : ActionWithValue* Value::getPntrToAction() {
     219    18718290 :   plumed_assert( action!=NULL );
     220    18718290 :   return action;
     221             : }
     222             : 
     223     4948688 : void Value::set(const std::size_t& n, const double& v ) {
     224     4948688 :   value_set=true;
     225     4948688 :   if( getRank()==0 ) { plumed_assert( n==0 ); data[n]=v; applyPeriodicity(n); }
     226     4923013 :   else if( !hasDeriv ) { plumed_dbg_massert( n<data.size(), "failing in " + getName() ); data[n]=v; applyPeriodicity(n); }
     227      317316 :   else { data[n*(1+ngrid_der)] = v; }
     228     4948688 : }
     229             : 
     230       74915 : void Value::push_back( const double& v ) {
     231       74915 :   value_set=true;
     232       74915 :   if( shape.size()==1 ) {
     233       36356 :     data.push_back(v); shape[0]++;
     234       38559 :   } else if( shape.size()==2 ) {
     235       38559 :     data.push_back(v);
     236       38559 :     shape[0] = std::ceil( data.size() / shape[1] );
     237             :   }
     238       74915 : }
     239             : 
     240     3223564 : std::size_t Value::getIndexInStore( const std::size_t& ival ) const {
     241     3223564 :   if( shape.size()==2 && ncols<shape[1] ) {
     242     2712800 :     unsigned irow = std::floor( ival / shape[1] ), jcol = ival%shape[1];
     243    11801200 :     for(unsigned i=0; i<getRowLength(irow); ++i) {
     244    11801200 :       if( getRowIndex(irow,i)==jcol ) return irow*ncols+i;
     245             :     }
     246           0 :     plumed_merror("cannot get store index");
     247             :   }
     248      510764 :   return ival;
     249             : }
     250             : 
     251   557639715 : double Value::get(const std::size_t& ival, const bool trueind) const {
     252   557639715 :   if( hasDeriv ) return data[ival*(1+ngrid_der)];
     253             : #ifdef DNDEBUG
     254             :   if( action ) plumed_dbg_massert( ival<getNumberOfValues(), "could not get value from " + name );
     255             : #endif
     256   542726522 :   if( shape.size()==2 && ncols<shape[1] && trueind ) {
     257     2000962 :     unsigned irow = std::floor( ival / shape[1] ), jcol = ival%shape[1];
     258             :     // This is a special treatment for the lower triangular matrices that are used when
     259             :     // we do ITRE with COLLECT_FRAMES
     260     2000962 :     if( ncols==0 ) {
     261           0 :       if( jcol<=irow ) return data[0.5*irow*(irow+1) + jcol];
     262             :       return 0;
     263             :     }
     264    20655470 :     for(unsigned i=0; i<getRowLength(irow); ++i) {
     265    19057356 :       if( getRowIndex(irow,i)==jcol ) return data[irow*ncols+i];
     266             :     }
     267             :     return 0.0;
     268             :   }
     269   540725560 :   plumed_massert( ival<data.size(), "cannot get value from " + name );
     270   540725560 :   return data[ival];
     271             : }
     272             : 
     273    39305298 : void Value::addForce(const std::size_t& iforce, double f, const bool trueind) {
     274    39305298 :   hasForce=true;
     275    39305298 :   if( shape.size()==2 && !hasDeriv && ncols<shape[1] && trueind ) {
     276           0 :     unsigned irow = std::floor( iforce / shape[0] ), jcol = iforce%shape[0];
     277           0 :     for(unsigned i=0; i<getRowLength(irow); ++i) {
     278           0 :       if( getRowIndex(irow,i)==jcol ) { inputForce[irow*ncols+i]+=f; return; }
     279             :     }
     280           0 :     plumed_assert( fabs(f)<epsilon ); return;
     281             :   }
     282    39305298 :   plumed_massert( iforce<inputForce.size(), "can't add force to " + name );
     283    39305298 :   inputForce[iforce]+=f;
     284             : }
     285             : 
     286             : 
     287       12065 : void Value::buildDataStore( const bool forprint ) {
     288       12065 :   if( getRank()==0 ) return;
     289        5289 :   storedata=true; setShape( shape );
     290        5289 :   if( !forprint ) return ;
     291         192 :   ActionWithVector* av=dynamic_cast<ActionWithVector*>( action );
     292         192 :   if( av ) (av->getFirstActionInChain())->never_reduce_tasks=true;
     293             : }
     294             : 
     295        7082 : void Value::reshapeMatrixStore( const unsigned& n ) {
     296             :   plumed_dbg_assert( shape.size()==2 && !hasDeriv );
     297        7082 :   if( !storedata ) return ;
     298        7082 :   ncols=n; if( ncols>shape[1] ) ncols=shape[1];
     299        7082 :   unsigned size=shape[0]*ncols;
     300        7082 :   if( matrix_bookeeping.size()!=(size+shape[0]) ) {
     301        2562 :     data.resize( size ); inputForce.resize( size );
     302        2562 :     matrix_bookeeping.resize( size + shape[0], 0 );
     303        2562 :     if( ncols>=shape[1] ) {
     304      248548 :       for(unsigned i=0; i<shape[0]; ++i) {
     305      246049 :         matrix_bookeeping[(1+ncols)*i] = shape[1];
     306    23653127 :         for(unsigned j=0; j<shape[1]; ++j) matrix_bookeeping[(1+ncols)*i+1+j]=j;
     307             :       }
     308             :     }
     309             :   }
     310        7082 :   if( ncols<shape[1] ) std::fill(matrix_bookeeping.begin(), matrix_bookeeping.end(), 0);
     311             : }
     312             : 
     313       35389 : void Value::setPositionInMatrixStash( const unsigned& p ) {
     314             :   plumed_dbg_assert( shape.size()==2 && !hasDeriv );
     315       35389 :   matpos=p;
     316       35389 : }
     317             : 
     318    15992076 : bool Value::ignoreStoredValue(const std::string& c) const {
     319    15992076 :   if( !storedata && shape.size()>0 ) return true;
     320     6832518 :   ActionWithVector* av=dynamic_cast<ActionWithVector*>(action);
     321     6832518 :   if( av ) return (av->getFirstActionInChain())->getLabel()==c;
     322             :   return false;
     323             : }
     324             : 
     325        5510 : void Value::setConstant() {
     326        5510 :   valtype=constant; storedata=true; setShape( shape );
     327        5510 :   if( getRank()==2 && !hasDeriv ) reshapeMatrixStore( shape[1] );
     328        5510 : }
     329             : 
     330         456 : void Value::writeBinary(std::ostream&o) const {
     331         456 :   o.write(reinterpret_cast<const char*>(&data[0]),data.size()*sizeof(double));
     332         456 : }
     333             : 
     334        1274 : void Value::setSymmetric( const bool& sym ) {
     335        1274 :   plumed_assert( shape.size()==2 && !hasDeriv );
     336        1274 :   if( sym && shape[0]!=shape[1] ) plumed_merror("non-square matrix cannot be symmetric");
     337        1274 :   symmetric=sym;
     338        1274 : }
     339             : 
     340        2220 : void Value::retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& active, std::vector<double>& elems ) {
     341        2220 :   nedge=0; plumed_dbg_assert( shape.size()==2 && !hasDeriv );
     342             :   // Check we have enough space to store the edge list
     343        2220 :   if( elems.size()<shape[0]*ncols ) { elems.resize( shape[0]*ncols ); active.resize( shape[0]*ncols ); }
     344             : 
     345      269579 :   for(unsigned i=0; i<shape[0]; ++i) {
     346             :     unsigned ncol = getRowLength(i);
     347    15002116 :     for(unsigned j=0; j<ncol; ++j) {
     348    14734757 :       if( fabs(get(i*ncols+j,false))<epsilon ) continue;
     349     3265037 :       if( symmetric && getRowIndex(i,j)>i ) continue;
     350     3213497 :       active[nedge].first = i; active[nedge].second = getRowIndex(i,j);
     351     3213497 :       elems[nedge] = get(i*ncols+j,false); nedge++;
     352             :     }
     353             :   }
     354        2220 : }
     355             : 
     356         456 : void Value::readBinary(std::istream&i) {
     357         456 :   i.read(reinterpret_cast<char*>(&data[0]),data.size()*sizeof(double));
     358         456 : }
     359             : 
     360      775666 : void Value::convertIndexToindices(const std::size_t& index, std::vector<unsigned>& indices ) const {
     361      775666 :   if( hasDeriv || getRank()==1 ) {
     362       12302 :     std::size_t kk=index; indices[0]=index%shape[0];
     363       12302 :     for(unsigned i=1; i<shape.size()-1; ++i) {
     364           0 :       kk=(kk-indices[i-1])/shape[i-1];
     365           0 :       indices[i]=kk%shape[i];
     366             :     }
     367       12302 :     if(shape.size()>=2) indices[shape.size()-1]=(kk-indices[shape.size()-2])/shape[shape.size()-2];
     368      763364 :   } else if( getRank()==2 ) {
     369      763364 :     indices[0]=std::floor( index/shape[1] ); indices[1] = index%shape[1];
     370             :   }
     371      775666 : }
     372             : 
     373      548216 : void Value::print( OFile& ofile ) const {
     374      571290 :   if( isPeriodic() ) { ofile.printField( "min_" + name, str_min ); ofile.printField("max_" + name, str_max ); }
     375      548216 :   if( shape.size()==0 || getNumberOfValues()==1 ) {
     376      546897 :     ofile.printField( name, get(0) );
     377             :   } else {
     378        1319 :     std::vector<unsigned> indices( shape.size() );
     379      776985 :     for(unsigned i=0; i<getNumberOfValues(); ++i) {
     380      775666 :       convertIndexToindices( i, indices ); std::string num, fname = name;
     381     2314696 :       for(unsigned i=0; i<shape.size(); ++i) { Tools::convert( indices[i]+1, num ); fname += "." + num; }
     382      775666 :       ofile.printField( fname, get(i) );
     383             :     }
     384             :   }
     385      548216 : }
     386             : 
     387      232611 : unsigned Value::getGoodNumThreads( const unsigned& j, const unsigned& k ) const {
     388      232611 :   return OpenMP::getGoodNumThreads( &data[j], (k-j) );
     389             : }
     390             : 
     391     7526197 : bool Value::calculateOnUpdate() const {
     392     7526197 :   return (valtype==average || valtype==calcFromAverage);
     393             : }
     394             : 
     395         129 : std::string Value::getValueType() const {
     396         129 :   if( getRank()==0 ) return "scalar";
     397          96 :   if( getRank()>0 && hasDerivatives() ) return "grid";
     398          91 :   if( getRank()==1 ) return "vector";
     399          21 :   if( getRank()==2 ) return "matrix";
     400           0 :   plumed_merror("unknown type for value " + getName() );
     401             :   return "";
     402             : }
     403             : 
     404             : }