void StoreDataVessel::chainRule( const unsigned& ival, const std::vector<double>& df ){
  plumed_dbg_assert( getAction()->derivativesAreRequired() && df.size()==vecsize );
  // Clear final derivatives array
  final_derivatives.assign( final_derivatives.size(), 0.0 );

  if(getAction()->lowmem){
     plumed_dbg_assert( ival<max_lowmem_stash );
     unsigned maxder = getAction()->getNumberOfDerivatives();
     for(unsigned ider=0;ider<active_der[ival];++ider){
        final_derivatives[ider]=0.0;
        unsigned ibuf=ival*(vecsize*maxder) + ider; 
        for(unsigned jcomp=0;jcomp<vecsize;++jcomp){
            final_derivatives[ider]+=df[jcomp]*local_derivatives[ibuf]; 
            ibuf+=maxder;
        }
     }
  } else {
     plumed_dbg_assert( ival<getAction()->getFullNumberOfTasks() );
     for(unsigned ider=0;ider<active_der[ival];++ider){
         final_derivatives[ider]=0.0;
         unsigned ibuf=ival*(vecsize*nspace) + 1 + ider; 
         for(unsigned jcomp=0;jcomp<vecsize;++jcomp){
             final_derivatives[ider]+=df[jcomp]*getBufferElement(ibuf);
             ibuf+=nspace;
         }
     }
  }
}

void StoreDataVessel::storeValues( const unsigned& myelem, MultiValue& myvals, std::vector<double>& buffer ) const {
  plumed_dbg_assert( vecsize>0 );
  unsigned jelem = getAction()->getPositionInCurrentTaskList( myelem ); plumed_dbg_assert( jelem<getNumberOfStoredValues() );
  unsigned ibuf = bufstart + jelem * vecsize * nspace; 
  for(unsigned icomp=0;icomp<vecsize;++icomp){
      buffer[ibuf] += myvals.get(icomp); ibuf+=nspace;
  }
}

void PointWiseMapping::readRestOfFrame(){
  if( getFrame( getNumberOfReferenceFrames() - 1 )->isDirection() ) plumed_merror("cannot use directions in mapping");
  plumed_dbg_assert( property.size()>0 );

  std::vector<double> labelvals;
  if( !ispath ){
      labelvals.resize( property.size() );
      for(unsigned i=0;i<property.size();++i) parse( property[i], labelvals[i] );
  } else {
      labelvals.resize(1);
      labelvals[0]=static_cast<double>( frames.size() ); 
  }
  low_dim.push_back( labelvals ); 
  plumed_dbg_assert( low_dim.size()==getNumberOfReferenceFrames() );
}

double ArgumentOnlyDistance::calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg,
                                   ReferenceValuePack& myder, const bool& squared ) const {
    plumed_dbg_assert( pos.size()==0 );
    double d=calculateArgumentDistance( vals, arg, myder, squared );
    if( !myder.updateComplete() ) myder.updateDynamicLists();
    return d;
}

double MultiColvarBase::getCentralAtomDerivative( const unsigned& iatom, const unsigned& jcomp, const Vector& df ){
  plumed_dbg_assert( atomsWithCatomDer.isActive(iatom) && jcomp<3 );
  unsigned nder = 3*getNumberOfAtoms() + 9;
  return df[0]*getElementDerivative( (getCentralAtomElementIndex()+0)*nder + 3*iatom + jcomp ) +
         df[1]*getElementDerivative( (getCentralAtomElementIndex()+1)*nder + 3*iatom + jcomp ) +
         df[2]*getElementDerivative( (getCentralAtomElementIndex()+2)*nder + 3*iatom + jcomp ); 
}

void VectorMultiColvar::addForcesOnAtoms( const std::vector<double>& inforces ){
  plumed_dbg_assert( inforces.size()==getNumberOfDerivatives() );
  std::vector<double> oldforces( getNumberOfDerivatives() ); 
  getForcesFromVessels( oldforces ); 
  for(unsigned i=0;i<getNumberOfDerivatives();++i) oldforces[i]+=inforces[i];
  setForcesOnAtoms( oldforces );
}

double OptimalRMSD::projectAtomicDisplacementOnVector( const bool& normalized, const std::vector<Vector>& vecs, ReferenceValuePack& mypack ) const {
  plumed_dbg_assert( mypack.calcUsingPCAOption() );

  double proj=0.0; mypack.clear();
  for(unsigned i=0; i<vecs.size(); ++i) {
    proj += dotProduct( mypack.getAtomsDisplacementVector()[i], vecs[i] );
  }
  for(unsigned a=0; a<3; a++) {
    for(unsigned b=0; b<3; b++) {
      for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
        double tmp1=0.;
        for(unsigned n=0; n<getNumberOfAtoms(); n++) tmp1+=mypack.centeredpos[n][b]*vecs[n][a];
        if( normalized ) mypack.addAtomDerivatives( iat, getDisplace()[iat]*mypack.DRotDPos[a][b][iat]*tmp1 );
        else mypack.addAtomDerivatives( iat, mypack.DRotDPos[a][b][iat]*tmp1 );
      }
    }
  }
  Tensor trot=mypack.rot[0].transpose();
  Vector v1; v1.zero(); double prefactor = 1. / static_cast<double>( getNumberOfAtoms() );
  for(unsigned n=0; n<getNumberOfAtoms(); n++) v1+=prefactor*matmul(trot,vecs[n]);
  if( normalized ) {
    for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) mypack.addAtomDerivatives( iat, getDisplace()[iat]*(matmul(trot,vecs[iat]) - v1) );
  } else {
    for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) mypack.addAtomDerivatives( iat, (matmul(trot,vecs[iat]) - v1) );
  }
  if( !mypack.updateComplete() ) mypack.updateDynamicLists();

  return proj;
}

void MultiColvarBase::activateIndexes( const unsigned& istart, const unsigned& number, const std::vector<unsigned>& indexes ){
  plumed_assert( number>0 );
  for(unsigned i=0;i<number-9;i+=3){
      plumed_dbg_assert( indexes[istart+i]%3==0 ); unsigned iatom=indexes[istart+i]/3; 
      atoms_with_derivatives.activate( iatom ); 
  }
}

void ActionWithVessel::chainRuleForElementDerivatives( const unsigned& iout, const unsigned& ider, const unsigned& stride, 
                                                       const unsigned& off, const double& df, Vessel* valout ){
  plumed_dbg_assert( off<stride );
  current_buffer_stride=stride;
  current_buffer_start=valout->bufstart + stride*(getNumberOfDerivatives()+1)*iout + stride + off;
  mergeDerivatives( ider, df );
}

void ReferenceValuePack::copyScaledDerivatives( const unsigned& from, const double& scalef, const MultiValue& tvals ){
  plumed_dbg_assert( tvals.getNumberOfDerivatives()==myvals.getNumberOfDerivatives() );
  for(unsigned i=0;i<tvals.getNumberActive();++i){
      unsigned ider=tvals.getActiveIndex(i);
      myvals.addDerivative( oind, ider, scalef*tvals.getDerivative( from, ider ) );
  }
}

void BridgedMultiColvarFunction::clearDerivativesAfterTask( const unsigned& ider ){
  unsigned vstart=getNumberOfDerivatives()*ider;
  if( derivativesAreRequired() ){
     // Clear atom derivatives
     for(unsigned i=0;i<atoms_with_derivatives.getNumberActive();++i){
        unsigned iatom=vstart+3*atoms_with_derivatives[i];
        setElementDerivative( iatom, 0.0 ); iatom++;
        setElementDerivative( iatom, 0.0 ); iatom++;
        setElementDerivative( iatom, 0.0 );
     }
     // Clear virial contribution
     unsigned nvir=vstart+3*mycolv->getNumberOfAtoms();
     for(unsigned j=0;j<9;++j){
        setElementDerivative( nvir, 0.0 ); nvir++;
     }
     // Clear derivatives of local atoms
     for(unsigned j=0;j<getNumberOfAtoms();++j){
        setElementDerivative( nvir, 0.0 ); nvir++;
        setElementDerivative( nvir, 0.0 ); nvir++;
        setElementDerivative( nvir, 0.0 ); nvir++;
     }
     plumed_dbg_assert( (nvir-vstart)==getNumberOfDerivatives() );
  }
  // Clear values
  thisval_wasset[ider]=false; setElementValue( ider, 0.0 ); thisval_wasset[ider]=false;
}

double HistogramBead::calculateWithCutoff( double x, double& df ) const {
  plumed_dbg_assert(init && periodicity!=unset );

  double lowB, upperB, f;
  lowB = difference( x, lowb ) / width ; upperB = difference( x, highb ) / width;
  if( upperB<=-cutoff || lowB>=cutoff ) { df=0; return 0; }

  if( type==gaussian ) {
    lowB /= sqrt(2.0); upperB /= sqrt(2.0);
    df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( sqrt(2*pi)*width );
    f = 0.5*( erf( upperB ) - erf( lowB ) );
  } else if( type==triangular ) {
    df=0;
    if( fabs(lowB)<1. ) df = (1 - fabs(lowB)) / width;
    if( fabs(upperB)<1. ) df -= (1 - fabs(upperB)) / width;
    if (upperB<=-1. || lowB >=1.) {
      f=0.;
    } else {
      double ia, ib;
      if( lowB>-1.0 ) { ia=lowB; } else { ia=-1.0; }
      if( upperB<1.0 ) { ib=upperB; } else { ib=1.0; }
      f = (ib*(2.-fabs(ib))-ia*(2.-fabs(ia)))*0.5;
    }
  } else {
    plumed_merror("function type does not exist");
  }
  return f;
}

unsigned CInterpolation::search1( const unsigned& kk, const double& x, const unsigned& jold ) const {
    int inc=stride[kk], jl=jold*stride[kk], ju=(jold+1)*stride[kk], jm; 
    if ( x>=splinepoints(jl,kk) && x<splinepoints( ju, kk ) ) return jl;
    else {
        if( x>=splinepoints(jl, kk ) ){
            while(true){
                ju=jl+inc;
                if( x<splinepoints( ju, kk ) ) break;
                else if( ju>=(np[kk]-1)*inc ){ju=(np[kk]-1)*inc; break; } 
                jl=ju; 
            }
        }
        else{
            ju=jl;
            while(true){
                jl=jl-inc;
                if( x>=splinepoints( jl, kk ) ) break;
                else if( jl<=0 ){ jl=0; break; }
                ju=jl; 
            }
        }
    }
    while( ju-jl>inc ){
      jm = (ju+jl) / (2*inc) ;
      if ( x>splinepoints(jm*inc,kk) ) jl=jm*inc; else ju=jm*inc;
    }
    plumed_dbg_assert( jl%stride[kk]==0 && ju==jl+stride[kk] );
    return jl;
}

void MultiColvarBase::setAtomsForCentralAtom( const std::vector<bool>& catom_ind ){
  unsigned nat=0; plumed_assert( catom_ind.size()==ablocks.size() );
  for(unsigned i=0;i<catom_ind.size();++i){
      use_for_central_atom[i]=catom_ind[i]; 
      if( use_for_central_atom[i] ) nat++;
  }
  plumed_dbg_assert( nat>0 ); numberForCentralAtom = 1.0 / static_cast<double>( nat );
}

void VolumeGradientBase::addBridgeForces( const std::vector<double>& bb ) {
  plumed_dbg_assert( bb.size()==tmpforces.size()-9 );
  // Forces on local atoms
  for(unsigned i=0; i<bb.size(); ++i) tmpforces[i]=bb[i];
  // Virial contribution is zero
  for(unsigned i=bb.size(); i<bb.size()+9; ++i) tmpforces[i]=0.0;
  setForcesOnAtoms( tmpforces, 0 );
}

void MultiDomainRMSD::setupPCAStorage( ReferenceValuePack& mypack ) {
  plumed_dbg_assert( pcaIsEnabledForThisReference() );
  mypack.switchOnPCAOption();
  mypack.displacement.resize( getNumberOfAtoms() );
  mypack.centeredpos.resize( getNumberOfAtoms() );
  mypack.DRotDPos.resize(3,3); mypack.rot.resize( domains.size() );
  for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j) mypack.DRotDPos(i,j).resize( getNumberOfAtoms() );
}

void MultiColvarBase::decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const {
  plumed_dbg_assert( atoms.size()==ablocks.size() && !usespecies && ablocks.size()<4 );
  unsigned scode = taskCode;
  for(unsigned i=0;i<ablocks.size();++i){
      unsigned ind=( scode / decoder[i] );
      atoms[i] = ablocks[i][ind];
      scode -= ind*decoder[i];
  }
}

bool CoeffsBase::indicesExist(const std::vector<unsigned int>& indices) const {
  plumed_dbg_assert(indices.size()==ndimensions_);
  for(unsigned int k=0; k<ndimensions_; k++) {
    if(indices[k]>=indices_shape_[k]) {
      return false;
    }
  }
  return true;
}

double InterpolateCubic::get_fdf( const std::vector<double>& pos ){
  plumed_dbg_assert( pos.size()==1 );
  
  unsigned mybox=findBox( pos );
  double d1=ub[0] - lb[0]; 
  double b=( pos[0] - lb[0] ) / d1, a=( ub[0] - pos[0] ) / d1;
  
  double *cbase=&clist[(mybox*4)+3], *c3=cbase-1, *c2=c3-1, *c1=c2-1;
  double f=a*(*c1) + b*(*c2) + (a*a*a-a)*(*c3) + (b*b*b-b)*(*cbase);
  return f;  
}

void StoreCentralAtomsVessel::getIndexList( const unsigned& ntotal, const unsigned& jstore, const unsigned& maxder, std::vector<unsigned>& aindexes ){
  plumed_dbg_assert( mycolv->derivativesAreRequired() );

  aindexes[jstore]=3*mycolv->atomsWithCatomDer.getNumberActive();
  if( aindexes[jstore]>maxder ) error("too many derivatives to store. Run with LOWMEM");
  unsigned kder = ntotal + jstore*maxder;
  for(unsigned jder=0;jder<mycolv->atomsWithCatomDer.getNumberActive();++jder){
     unsigned iatom = 3*mycolv->atomsWithCatomDer[jder];
     for(unsigned icomp=0;icomp<3;++icomp){ aindexes[ kder ] = iatom+icomp; kder++; }
  }
}