size_t OFile::llwrite(const char*ptr,size_t s) {
  size_t r;
  if(linked) return linked->llwrite(ptr,s);
  if(! (comm && comm->Get_rank()>0)) {
    if(!fp) plumed_merror("writing on uninitilized File");
    if(gzfp) {
#ifdef __PLUMED_HAS_ZLIB
      r=gzwrite(gzFile(gzfp),ptr,s);
#else
      plumed_merror("file " + getPath() + ": trying to use a gz file without zlib being linked");
#endif
    } else {
      r=fwrite(ptr,1,s,fp);
    }
  }
//  This barrier is apparently useless since it comes
//  just before a Bcast.
//
//  Anyway, it looks like it is solving an issue that appeared on
//  TRAVIS (at least on my laptop) so I add it here.
//  GB
  if(comm) comm->Barrier();


  if(comm) comm->Bcast(r,0);
  return r;
}

void PlumedMain::update() {
  if(!active)return;

  stopwatch.start("6 Update");
// update step (for statistics, etc)
  updateFlags.push(true);
  for(const auto & p : actionSet) {
    p->beforeUpdate();
    if(p->isActive() && p->checkUpdate() && updateFlagsTop()) p->update();
  }
  while(!updateFlags.empty()) updateFlags.pop();
  if(!updateFlags.empty()) plumed_merror("non matching changes in the update flags");
// Check that no action has told the calculation to stop
  if(stopNow) {
    if(stopFlag) (*stopFlag)=1;
    else plumed_merror("your md code cannot handle plumed stop events - add a call to plumed.comm(stopFlag,stopCondition)");
  }

// flush by default every 10000 steps
// hopefully will not affect performance
// also if receive checkpointing signal
  if(step%10000==0||doCheckPoint) {
    fflush();
    log.flush();
    for(const auto & p : actionSet) p->fflush();
  }
  stopwatch.stop("6 Update");
}

OFile& OFile::open(const std::string&path) {
  plumed_assert(!cloned);
  eof=false;
  err=false;
  fp=NULL;
  gzfp=NULL;
  this->path=path;
  this->path=appendSuffix(path,getSuffix());
  if(checkRestart()) {
    fp=std::fopen(const_cast<char*>(this->path.c_str()),"a");
    mode="a";
    if(Tools::extension(this->path)=="gz") {
#ifdef __PLUMED_HAS_ZLIB
      gzfp=(void*)gzopen(const_cast<char*>(this->path.c_str()),"a9");
#else
      plumed_merror("file " + getPath() + ": trying to use a gz file without zlib being linked");
#endif
    }
  } else {
    backupFile( backstring, this->path );
    if(comm)comm->Barrier();
    fp=std::fopen(const_cast<char*>(this->path.c_str()),"w");
    mode="w";
    if(Tools::extension(this->path)=="gz") {
#ifdef __PLUMED_HAS_ZLIB
      gzfp=(void*)gzopen(const_cast<char*>(this->path.c_str()),"w9");
#else
      plumed_merror("file " + getPath() + ": trying to use a gz file without zlib being linked");
#endif
    }
  }
  if(plumed) plumed->insertFile(*this);
  return *this;
}

void PlumedMain::update(){
  if(!active)return;

  stopwatch.start("6 Update");
// update step (for statistics, etc)
  updateFlags.push(true);
  for(ActionSet::iterator p=actionSet.begin();p!=actionSet.end();++p){
    (*p)->beforeUpdate();
    if((*p)->isActive() && (*p)->checkUpdate() && updateFlagsTop()) (*p)->update();
  }
  while(!updateFlags.empty()) updateFlags.pop();
  if(!updateFlags.empty()) plumed_merror("non matching changes in the update flags");
// Check that no action has told the calculation to stop
  if(stopNow){
     if(stopFlag) (*stopFlag)=1;
     else plumed_merror("your md code cannot handle plumed stop events - add a call to plumed.comm(stopFlag,stopCondition)");
  }  

// flush by default every 10000 steps
// hopefully will not affect performance
  if(step%10000==0){
    fflush();
    log.flush();
    for(ActionSet::const_iterator p=actionSet.begin();p!=actionSet.end();++p) (*p)->fflush();
  }
  stopwatch.stop("6 Update");
}

void CoeffsBase::checkCoeffsInfo(const std::string& msg_header, const std::string& coeffs_type_f, const unsigned int ndimensions_f, const size_t ncoeffs_total_f, const std::vector<unsigned int>& indices_shape_f) {

  if(coeffs_type_f != getTypeStr()) {
    std::string msg = msg_header + " coeffs type " + coeffs_type_f + " from file doesn't match the defined value " + getTypeStr();
    plumed_merror(msg);
  }
  if(ndimensions_f != numberOfDimensions() ) {
    std::string s1; Tools::convert(ndimensions_f,s1);
    std::string s2; Tools::convert(numberOfDimensions(),s2);
    std::string msg = msg_header + " the number of dimensions " + s1 + " in file doesn't match the defined value " + s2;
    plumed_merror(msg);
  }
  if(ncoeffs_total_f != numberOfCoeffs() ) {
    std::string s1; Tools::convert(ncoeffs_total_f,s1);
    std::string s2; Tools::convert(numberOfCoeffs(),s2);
    std::string msg = msg_header + " the number of coeffs " + s1 + " in file doesn't match the defined value " + s2;
    plumed_merror(msg);
  }
  for(unsigned int k=0; k<numberOfDimensions(); k++) {
    if(indices_shape_f[k] != shapeOfIndices(k) ) {
      std::string s1; Tools::convert(indices_shape_f[k],s1);
      std::string s2; Tools::convert(shapeOfIndices(k),s2);
      std::string msg = msg_header + " for dimension labeled " + getDimensionLabel(k) + " the shape of indices " + s1 + " in file doesn't match defined value " + s2;
      plumed_merror(msg);
    }
  }
}

void PlumedMain::load(const std::string& ss) {
  if(DLLoader::installed()) {
    string s=ss;
    size_t n=s.find_last_of(".");
    string extension="";
    string base=s;
    if(n!=std::string::npos && n<s.length()-1) extension=s.substr(n+1);
    if(n!=std::string::npos && n<s.length())   base=s.substr(0,n);
    if(extension=="cpp") {
// full path command, including environment setup
// this will work even if plumed is not in the execution path or if it has been
// installed with a name different from "plumed"
      string cmd=config::getEnvCommand()+" \""+config::getPlumedRoot()+"\"/scripts/mklib.sh "+s;
      log<<"Executing: "<<cmd;
      if(comm.Get_size()>0) log<<" (only on master node)";
      log<<"\n";
      if(comm.Get_rank()==0) system(cmd.c_str());
      comm.Barrier();
      base="./"+base;
    }
    s=base+"."+config::getSoExt();
    void *p=dlloader.load(s);
    if(!p) {
      const std::string error_msg="I cannot load library " + ss + " " + dlloader.error();
      log<<"ERROR\n";
      log<<error_msg<<"\n";
      plumed_merror(error_msg);
    }
    log<<"Loading shared library "<<s.c_str()<<"\n";
    log<<"Here is the new list of available actions\n";
    log<<actionRegister();
  } else plumed_merror("loading not enabled, please recompile with -D__PLUMED_HAS_DLOPEN");
}

void OFile::backupFile( const std::string& bstring, const std::string& fname ) {
  if(fname=="/dev/null") return;
  int maxbackup=100;
  if(std::getenv("PLUMED_MAXBACKUP")) Tools::convert(std::getenv("PLUMED_MAXBACKUP"),maxbackup);
  if(maxbackup>0 && (!comm || comm->Get_rank()==0)) {
    FILE* ff=std::fopen(const_cast<char*>(fname.c_str()),"r");
    if(ff) {
      std::fclose(ff);
      std::string backup;
      size_t found=fname.find_last_of("/\\");
      std::string directory=fname.substr(0,found+1);
      std::string file=fname.substr(found+1);
      for(int i=0;; i++) {
        std::string num;
        Tools::convert(i,num);
        if(i>maxbackup) plumed_merror("cannot backup file "+file+" maximum number of backup is "+num+"\n");
        backup=directory+bstring +"."+num+"."+file;
        FILE* fff=std::fopen(backup.c_str(),"r");
        if(!fff) break;
        else std::fclose(fff);
      }
      int check=rename(fname.c_str(),backup.c_str());
      plumed_massert(check==0,"renaming "+fname+" into "+backup+" failed for reason: "+strerror(errno));
    }
  }
}

Value* ActionWithValue::copyOutput( const std::string& name ) const {
  for(unsigned i=0; i<values.size(); ++i) {
    if (values[i]->name==name) return values[i];
  }
  plumed_merror("there is no pointer with name " + name);
  return NULL;
}

unsigned PointWiseMapping::getPropertyIndex( const std::string& name ) const {
  for(unsigned i=0;i<property.size();++i){
     if( name==property[i] ) return i;
  }
  plumed_merror("no property with name " + name + " found");
  return 0;
}

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;
}

void MolDataClass::getBackboneForResidue( const std::string& type, const unsigned& residuenum, const PDB& mypdb, std::vector<AtomNumber>& atoms ){
  std::string residuename=mypdb.getResidueName( residuenum );
  plumed_massert( MolDataClass::allowedResidue( type, residuename ), "residue " + residuename + " unrecognized for molecule type " + type );
  if( type=="protein" ){
     if( residuename=="GLY"){
         atoms.resize(5);
         atoms[0]=mypdb.getNamedAtomFromResidue("N",residuenum);
         atoms[1]=mypdb.getNamedAtomFromResidue("CA",residuenum);
         atoms[2]=mypdb.getNamedAtomFromResidue("HA1",residuenum);
         atoms[3]=mypdb.getNamedAtomFromResidue("C",residuenum);
         atoms[4]=mypdb.getNamedAtomFromResidue("O",residuenum);
     } else if( residuename=="ACE"){
         atoms.resize(1); 
         atoms[0]=mypdb.getNamedAtomFromResidue("C",residuenum);
     } else if( residuename=="NME"){
         atoms.resize(1); 
         atoms[0]=mypdb.getNamedAtomFromResidue("N",residuenum);
     } else {
         atoms.resize(5);
         atoms[0]=mypdb.getNamedAtomFromResidue("N",residuenum);
         atoms[1]=mypdb.getNamedAtomFromResidue("CA",residuenum); 
         atoms[2]=mypdb.getNamedAtomFromResidue("CB",residuenum);
         atoms[3]=mypdb.getNamedAtomFromResidue("C",residuenum);
         atoms[4]=mypdb.getNamedAtomFromResidue("O",residuenum);
     }
  } else {
     plumed_merror(type + " is not a valid molecule type");
  }
}

void PlumedMain::load(const std::string& ss){
  if(DLLoader::installed()){
     string s=ss;
     size_t n=s.find_last_of(".");
     string extension="";
     string base=s;
     if(n!=std::string::npos && n<s.length()-1) extension=s.substr(n+1);
     if(n!=std::string::npos && n<s.length())   base=s.substr(0,n);
     if(extension=="cpp"){
       string cmd="plumed mklib "+s;
       log<<"Executing: "<<cmd;
       if(comm.Get_size()>0) log<<" (only on master node)";
       log<<"\n";
       if(comm.Get_rank()==0) system(cmd.c_str());
       comm.Barrier();
       base="./"+base;
     }
     s=base+"."+config::getSoExt();
     void *p=dlloader.load(s);
     if(!p){
       log<<"ERROR\n";
       log<<"I cannot load library "<<ss<<"\n";
       log<<dlloader.error();
       log<<"\n";
       this->exit(1);
     }
     log<<"Loading shared library "<<s.c_str()<<"\n";
     log<<"Here is the new list of available actions\n";
     log<<actionRegister();
  } else plumed_merror("loading not enabled, please recompile with -D__PLUMED_HAS_DLOPEN");
}

std::string SwitchingFunction::description() const {
  std::ostringstream ostr;
  ostr<<1./invr0<<".  Using ";
  if(type==rational){
     ostr<<"rational";
  } else if(type==exponential){
     ostr<<"exponential";
  } else if(type==gaussian){
     ostr<<"gaussian";
  } else if(type==smap){
     ostr<<"smap";
  } else if(type==cubic){
     ostr<<"cubic";
  } else{
     plumed_merror("Unknown switching function type");
  }
  ostr<<" swiching function with parameters d0="<<d0;
  if(type==rational){
    ostr<<" nn="<<nn<<" mm="<<mm;
  } else if(type==smap){
    ostr<<" a="<<a<<" b="<<b;
  } else if(type==cubic){
    ostr<<" dmax="<<dmax;
  }
  return ostr.str(); 
}

void LandmarkSelectionBase::checkRead() const {
  if(!input.empty()) {
    std::string msg="cannot understand the following words from landmark selection input : ";
    for(unsigned i=0; i<input.size(); ++i) msg = msg + input[i] + ", ";
    plumed_merror(msg);
  }
}

void PlumedMain::readInputWords(const std::vector<std::string> & words) {
  plumed_assert(initialized);
  if(words.empty())return;
  else if(words[0]=="ENDPLUMED") return;
  else if(words[0]=="_SET_SUFFIX") {
    plumed_assert(words.size()==2);
    setSuffix(words[1]);
  } else {
    std::vector<std::string> interpreted(words);
    Tools::interpretLabel(interpreted);
    Action* action=actionRegister().create(ActionOptions(*this,interpreted));
    if(!action) {
      log<<"ERROR\n";
      log<<"I cannot understand line:";
      for(unsigned i=0; i<interpreted.size(); ++i) log<<" "<<interpreted[i];
      log<<"\n";
      log.flush();
      plumed_merror("I cannot understand line " + interpreted[0] + " " + interpreted[1]);
    };
    action->checkRead();
    actionSet.push_back(action);
  };

  pilots=actionSet.select<ActionPilot*>();
}

AtomNumber PDB::getNamedAtomFromResidue( const std::string& aname, const unsigned& resnum ) const {
  for(unsigned i=0;i<size();++i){
     if( residue[i]==resnum && atomsymb[i]==aname ) return numbers[i];
  }
  std::string num; Tools::convert( resnum, num );
  plumed_merror("residue " + num + " does not contain an atom named " + aname );
  return numbers[0]; // This is to stop compiler errors
}

int ActionWithValue::getComponent( const std::string& name ) const {
  plumed_massert( !exists( getLabel() ), "You should not be calling this routine if you are using a value");
  std::string thename; thename=getLabel() + "." + name;
  for(unsigned i=0; i<values.size(); ++i) {
    if (values[i]->name==thename) return i;
  }
  plumed_merror("there is no component with name " + name);
  return -1;
}

size_t OFile::llwrite(const char*ptr,size_t s){
  size_t r;
  if(linked) return linked->llwrite(ptr,s);
  if(! (comm && comm->Get_rank()>0)){
    if(!fp) plumed_merror("writing on uninitilized File");
    if(gzfp){
#ifdef __PLUMED_HAS_ZLIB
      r=gzwrite(gzFile(gzfp),ptr,s);
#else
      plumed_merror("trying to use a gz file without zlib being linked");
#endif
    } else {
      r=fwrite(ptr,1,s,fp);
    }
  }
  if(comm) comm->Bcast(r,0);
  return r;
}

void BridgeVessel::completeNumericalDerivatives(){
  unsigned nextra = myOutputAction->getNumberOfDerivatives() - getAction()->getNumberOfDerivatives();
  Matrix<double> tmpder( myOutputValues->getNumberOfComponents(), nextra );
  ActionWithVessel* vval=dynamic_cast<ActionWithVessel*>( myOutputAction );
  for(unsigned i=0;i<nextra;++i){
      vval->bridgeVariable=i; getAction()->calculate();
      for(int j=0;j<myOutputValues->getNumberOfComponents();++j) tmpder(j,i) = myOutputValues->getOutputQuantity(j);
  }
  vval->bridgeVariable=nextra; getAction()->calculate(); 
  plumed_assert( inum==mynumerical_values.size() ); inum=0;  // Reset inum now that we have finished calling calculate
  std::vector<double> base( myOutputValues->getNumberOfComponents() );
  for(int j=0;j<myOutputValues->getNumberOfComponents();++j) base[j] = myOutputValues->getOutputQuantity(j);

  const double delta=sqrt(epsilon);
  ActionAtomistic* aa=dynamic_cast<ActionAtomistic*>( getAction() );
  unsigned nvals=myOutputValues->getNumberOfComponents();
  for(unsigned j=0;j<nvals;++j) ( myOutputValues->copyOutput(j) )->clearDerivatives();   

  if( aa ){
      ActionWithArguments* aarg=dynamic_cast<ActionWithArguments*>( getAction() );
      plumed_assert( !aarg ); Tensor box=aa->getBox(); 
      unsigned natoms=aa->getNumberOfAtoms();
      for(unsigned j=0;j<nvals;++j){
          double ref=( myOutputValues->copyOutput(j) )->get();
          if( ( myOutputValues->copyOutput(j) )->getNumberOfDerivatives()>0 ){
              for(unsigned i=0;i<3*natoms;++i){
                  double d=( mynumerical_values[i*nvals+j] - ref)/delta;
                  ( myOutputValues->copyOutput(j) )->addDerivative(i,d);
              }
              Tensor virial;
              for(int i=0;i<3;i++) for(int k=0;k<3;k++){
                 virial(i,k)=( mynumerical_values[ nvals*(3*natoms + 3*i + k) + j ]-ref)/delta;
              }
              virial=-matmul(box.transpose(),virial);
              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));
          }
      }
  } else {
      plumed_merror("not implemented or tested yet");
//      unsigned nder=myOutputAction->getNumberOfDerivatives();
//      for(unsigned j=0;j<nvals;++j){
//          double ref=( myOutputValues->copyOutput(j) )->get();
//              for(unsigned i=0;i<nder;++i){
//                  double d=( mynumerical_values[i*nvals+j] - ref)/delta;
//                  ( myOutputValues->copyOutput(j) )->addDerivative(i,d);
//              }
//          }
//      }
  }
  // Add the derivatives wrt to the local quantities we are working with
  for(unsigned j=0;j<nvals;++j){
     unsigned k=0;
     for(unsigned i=getAction()->getNumberOfDerivatives();i<myOutputAction->getNumberOfDerivatives();++i){
        ( myOutputValues->copyOutput(j) )->addDerivative( i, (tmpder(j,k)-base[j])/sqrt(epsilon) ); k++;
     }
  }
}

void Communicator::Set_comm(void*val){
#ifdef __PLUMED_MPI
 plumed_massert(initialized(),"you are trying to use an MPI function, but MPI is not initialized");
 if(val) Set_comm(*(MPI_Comm*)val);
#else
 (void) val;
 plumed_merror("you are trying to use an MPI function, but PLUMED has been compiled without MPI support");
#endif
}