void ReferenceArguments::readArgumentsFromPDB( const PDB& pdb ) {
  ReferenceAtoms* aref=dynamic_cast<ReferenceAtoms*>( this );
  arg_names.resize( pdb.getArgumentNames().size() );
  for(unsigned i=0; i<arg_names.size(); ++i) arg_names[i]=pdb.getArgumentNames()[i];
  if( !aref && arg_names.size()==0 ) error("no arguments in input PDB file");

  reference_args.resize( arg_names.size() ); arg_der_index.resize( arg_names.size() );
  for(unsigned i=0; i<arg_names.size(); ++i) {
    if( !pdb.getArgumentValue(arg_names[i], reference_args[i]) ) error("argument " + arg_names[i] + " was not set in pdb input");
    arg_der_index[i]=i;
  }

  if( hasweights ) {
    plumed_massert( !hasmetric, "should not have weights if we are using metric");
    weights.resize( arg_names.size() ); sqrtweight.resize( arg_names.size() );
    for(unsigned i=0; i<reference_args.size(); ++i) {
      if( !pdb.getArgumentValue("sigma_" + arg_names[i], weights[i]) ) error("value sigma_" + arg_names[i] + " was not set in pdb input");
      sqrtweight[i] = sqrt( weights[i] );
    }
  } else if( hasmetric ) {
    plumed_massert( !hasweights, "should not have weights if we are using metric");
    double thissig; metric.resize( arg_names.size(), arg_names.size() );
    for(unsigned i=0; i<reference_args.size(); ++i) {
      for(unsigned j=i; j<reference_args.size(); ++j) {
        if( !pdb.getArgumentValue("sigma_" + arg_names[i] + "_" + arg_names[j], thissig) ) {
          error("value sigma_" + arg_names[i] + "_" + arg_names[j] + " was not set in pdb input");
        }
        metric(i,j)=metric(j,i)=thissig;
      }
    }
  } else {
    weights.resize( arg_names.size() ); sqrtweight.resize( arg_names.size() );
    for(unsigned i=0; i<weights.size(); ++i) sqrtweight[i]=weights[i]=1.0;
  }
}

void HistogramBead::generateBins( const std::string& params, std::vector<std::string>& bins ) {
  std::vector<std::string> data=Tools::getWords(params);
  plumed_massert(data.size()>=1,"There is no input for this keyword");

  std::string name=data[0];

  unsigned nbins; std::vector<double> range(2); std::string smear;
  bool found_nb=Tools::parse(data,"NBINS",nbins);
  plumed_massert(found_nb,"Number of bins in histogram not found");
  bool found_r=Tools::parse(data,"LOWER",range[0]);
  plumed_massert(found_r,"Lower bound for histogram not specified");
  found_r=Tools::parse(data,"UPPER",range[1]);
  plumed_massert(found_r,"Upper bound for histogram not specified");
  plumed_massert(range[0]<range[1],"Range specification is dubious");
  bool found_b=Tools::parse(data,"SMEAR",smear);
  if(!found_b) { Tools::convert(0.5,smear); }

  std::string lb,ub; double delr = ( range[1]-range[0] ) / static_cast<double>( nbins );
  for(unsigned i=0; i<nbins; ++i) {
    Tools::convert( range[0]+i*delr, lb );
    Tools::convert( range[0]+(i+1)*delr, ub );
    bins.push_back( name + " " +  "LOWER=" + lb + " " + "UPPER=" + ub + " " + "SMEAR=" + smear );
  }
  plumed_assert(bins.size()==nbins);
}

void BridgeVessel::setOutputAction( ActionWithVessel* myact ){
  ActionWithValue* checkme=dynamic_cast<ActionWithValue*>( getAction() );
  plumed_massert( checkme, "vessel in bridge must inherit from ActionWithValue");

  myOutputAction=myact;
  myOutputValues=dynamic_cast<ActionWithValue*>( myact );
  plumed_massert( myOutputValues, "bridging vessel must inherit from ActionWithValue");
}

void CLTool::parseFlag( const std::string&key, bool&t ){
  plumed_massert(keywords.exists(key),"keyword " + key + " has not been registered");
  plumed_massert(keywords.style(key,"flag"),"keyword " + key + " has not been registered as a flag");
  plumed_assert(inputData.count(key)>0); 
  if( inputData[key]=="true") t=true;
  else if( inputData[key]=="false") t=false;
  else plumed_error();
}

void RMSD::setReference(const vector<Vector> & reference){
  unsigned n=reference.size();
  this->reference=reference;
  plumed_massert(align.empty(),"you should first clear() an RMSD object, then set a new referece");
  plumed_massert(displace.empty(),"you should first clear() an RMSD object, then set a new referece");
  align.resize(n,1.0);
  displace.resize(n,1.0);
}

void  BiasRepresentation::addGrid( vector<string> gmin, vector<string> gmax, vector<unsigned> nbin ){
    plumed_massert(hills.size()==0,"you can set the grid before loading the hills");
    plumed_massert(hasgrid==false,"to build the grid you should not having the grid in this bias representation");
    string ss; ss="file.free"; 
    vector<Value*> vv;for(unsigned i=0;i<values.size();i++)vv.push_back(values[i]);
    //cerr<<" initializing grid "<<endl;
    BiasGrid_=new Grid(ss,vv,gmin,gmax,nbin,false,true);
    hasgrid=true;	
}

CoeffsBase::CoeffsBase(
  const std::string& label,
  std::vector<std::vector<Value*> >& multicoeffs_args,
  std::vector<std::vector<BasisFunctions*> >& multicoeffs_basisf,
  const bool use_iteration_counter,
  const std::string& multicoeffs_label):
  label_(label),
  data_label_(label),
  coeffs_type_(MultiCoeffs_LinearBasisSet),
  iteration_and_time_active_(use_iteration_counter),
  iteration_opt(0),
  time_md(-1.0),
  active(true),
  action_pntr_(NULL),
  vesbias_pntr_(NULL),
  ndimensions_(0),
  indices_shape_(0),
  ncoeffs_(0),
  coeffs_descriptions_(0),
  dimension_labels_(0),
  args_(0),
  basisf_(0),
  multicoeffs_(true),
  multicoeffs_args_(multicoeffs_args),
  multicoeffs_basisf_(multicoeffs_basisf),
  field_type_("type"),
  field_ndimensions_("ndimensions"),
  field_ncoeffs_total_("ncoeffs_total"),
  field_shape_prefix_("shape_"),
  field_time_("time"),
  field_iteration_("iteration"),
  output_fmt_("%30.16e")
{
  plumed_massert(multicoeffs_args.size()==multicoeffs_basisf.size(),"Multi Coeffs: number of arguments vectors does not match number of basis functions vectors");
  unsigned int num_args = multicoeffs_args[0].size();
  unsigned int dim = num_args+1;
  std::vector<std::string> dimension_labels(dim);
  std::vector<unsigned int> indices_shape(dim);
  for(unsigned int i=0; i<num_args; i++) {
    std::string ip;
    Tools::convert(i+1,ip);
    dimension_labels[i] = "bf" + ip;
    indices_shape[i] = multicoeffs_basisf[0][i]->getNumberOfBasisFunctions();
  }
  indices_shape[dim-1] = multicoeffs_args.size();
  dimension_labels[dim-1] = multicoeffs_label;
  for(unsigned int k=0; k<multicoeffs_args.size(); k++) {
    plumed_massert(multicoeffs_args[k].size()==num_args && multicoeffs_basisf[k].size()==num_args,"Multi Coeffs: arguments and basis functions vectors for each bias should be of the same size");
    for(unsigned int i=0; i<num_args; i++) {
      plumed_massert(indices_shape[i]==multicoeffs_basisf[k][i]->getNumberOfBasisFunctions(),"Multi Coeffs: the coeffs shape for each bias should be identical");
    }
  }
  initializeIndices(indices_shape,dimension_labels);
  setupBasisFunctionsInfo();
}

void CoeffsBase::reinitializeIndices(const std::vector<unsigned int>& indices_shape_new) {
  plumed_massert(indices_shape_.size()>0,"indices must have been previously initialized before using this function");
  plumed_massert(dimension_labels_.size()>0,"indices must have been previously initialized before using this function");
  plumed_massert(indices_shape_new.size()==numberOfDimensions(),"when resizeing Coeffs the dimension must be constant");
  indices_shape_=indices_shape_new;
  ncoeffs_=1;
  for(unsigned int i=0; i<ndimensions_; i++) {
    ncoeffs_*=indices_shape_[i];
  }
  coeffs_descriptions_.clear();
  coeffs_descriptions_.resize(ncoeffs_);
}

ZpathVessel::ZpathVessel( const vesselbase::VesselOptions& da ):
FunctionVessel(da)
{
  Mapping* mymap=dynamic_cast<Mapping*>( getAction() );
  plumed_massert( mymap, "ZpathVessel should only be used with mappings");
  invlambda = 1.0 / mymap->getLambda(); usetol=true;
}

OFile& OFile::rewind() {
// we use here "hard" rewind, which means close/reopen
// the reason is that normal rewind does not work when in append mode
// moreover, we can take a backup of the file
  plumed_assert(fp);
  clearFields();
  if(gzfp) {
#ifdef __PLUMED_HAS_ZLIB
    gzclose((gzFile)gzfp);
#endif
  } else fclose(fp);
  if(!comm || comm->Get_rank()==0) {
    std::string fname=this->path;
    size_t found=fname.find_last_of("/\\");
    std::string directory=fname.substr(0,found+1);
    std::string file=fname.substr(found+1);
    std::string backup=directory+backstring +".last."+file;
    int check=rename(fname.c_str(),backup.c_str());
    plumed_massert(check==0,"renaming "+fname+" into "+backup+" failed for reason: "+strerror(errno));
  }
  if(gzfp) {
#ifdef __PLUMED_HAS_ZLIB
    gzfp=(void*)gzopen(const_cast<char*>(this->path.c_str()),"w9");
#endif
  } else fp=std::fopen(const_cast<char*>(path.c_str()),"w");
  return *this;
}

void ActionWithValue::addComponentWithDerivatives( const std::string& name ) {
  if( !keywords.outputComponentExists(name,true) ) {
    warning("a description of component " + name + " has not been added to the manual. Components should be registered like keywords in "
            "registerKeywords as described in the developer doc.");
  }
  std::string thename; thename=getLabel() + "." + name;
  for(unsigned i=0; i<values.size(); ++i) {
    plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components");
    plumed_massert(values[i]->name!=thename&&name!="bias","Since PLUMED 2.3 the component 'bias' is automatically added to all biases by the general constructor!\n"
                   "Remove the line addComponentWithDerivatives(\"bias\") from your bias.");
    plumed_massert(values[i]->name!=thename,"there is already a value with this name");
  }
  values.push_back(new Value(this,thename, true ) );
  std::string msg="  added component to this action:  "+thename+" \n";
  log.printf(msg.c_str());
}

InterpolateBicubic::InterpolateBicubic( const std::vector<unsigned>& dd, const std::vector<double>& fmin, const std::vector<double>& fmax ) :
CInterpolation(dd,fmin,fmax)
{
  plumed_massert(np.size()==2,"should be two dimensional data");
  static int wt_d[16*16]=
     {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
     -3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1, 0, 0, 0, 0,
     2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0,
     0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
     0, 0, 0, 0,-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1,
     0, 0, 0, 0, 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1,
     -3, 3, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0,
     9,-9, 9,-9, 6, 3,-3,-6, 6,-6,-3, 3, 4, 2, 1, 2,
     -6, 6,-6, 6,-4,-2, 2, 4,-3, 3, 3,-3,-2,-1,-1,-2,
     2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0,
     -6, 6,-6, 6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1,
     4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1};

  // This is to set up the coefficient matrix
  unsigned l=0; wt.resize(16,16); t1.resize(16); t2.resize(16);
  for (unsigned i=0;i<16;i++) for (unsigned j=0;j<16;j++){ wt(i,j)=wt_d[l++]; }
  // Resize everything
  dcross.resize( np[0], np[1] ); clist.resize( np[0] * np[1] * 4 * 4 );
}

void Colvar::requestAtoms(const vector<AtomNumber> & a){
  plumed_massert(!isEnergy,"request atoms should not be called if this is energy");
// Tell actionAtomistic what atoms we are getting
  ActionAtomistic::requestAtoms(a);
// Resize the derivatives of all atoms
  for(int i=0;i<getNumberOfComponents();++i) getPntrToComponent(i)->resizeDerivatives(3*a.size()+9);
}

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 SecondaryStructureRMSD::setSecondaryStructure( std::vector<Vector>& structure, double bondlength, double units ){
  // If we are in natural units get conversion factor from nm into natural length units
  if( plumed.getAtoms().usingNaturalUnits() ){
      error("cannot use this collective variable when using natural units");
  }
  plumed_massert( !(align_strands && align_atom_1==0 && align_atom_2==0), "you must use setAtomsFromStrands with strands cutoff"); 

  // Convert into correct units
  for(unsigned i=0;i<structure.size();++i){
     structure[i][0]*=units; structure[i][1]*=units; structure[i][2]*=units;
  }

  if( references.size()==0 ){
     finishTaskListUpdate();

     readVesselKeywords();
     if( getNumberOfVessels()==0 ){
         double r0; parse("R_0",r0); double d0; parse("D_0",d0);
         int nn; parse("NN",nn); int mm; parse("MM",mm);
         std::ostringstream ostr;
         ostr<<"RATIONAL R_0="<<r0<<" D_0="<<d0<<" NN="<<nn<<" MM="<<mm;
         std::string input=ostr.str(); addVessel( "LESS_THAN", input, -1 ); // -1 here means that this value will be named getLabel()
         readVesselKeywords();  // This makes sure resizing is done
     } 
  }

  // Set the reference structure
  references.push_back( metricRegister().create<SingleDomainRMSD>( alignType ) ); 
  unsigned nn=references.size()-1;
  std::vector<double> align( structure.size(), 1.0 ), displace( structure.size(), 1.0 );
  references[nn]->setBoundsOnDistances( true , bondlength );  // We always use pbc
  references[nn]->setReferenceAtoms( structure, align, displace );
//  references[nn]->setNumberOfAtoms( structure.size() );
}

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

IFile& IFile::getline(std::string &str){
  char tmp=0;
  str="";
// I comment out this (see note below on commented fgetpos):
// fpos_t pos;
// fgetpos(fp,&pos);
  while(llread(&tmp,1)==1 && tmp && tmp!='\n' && tmp!='\r' && !eof && !err){
    str+=tmp;
  }
  if(tmp=='\r'){
    llread(&tmp,1);
    plumed_massert(tmp=='\n',"plumed only accepts \\n (unix) or \\r\\n (dos) new lines");
  }
  if(eof){
    if(str.length()>0) eof=false;
  } else if(err || tmp!='\n'){
    eof = true;
    str="";
// there was a fsetpos here that apparently is not necessary
//  fsetpos(fp,&pos);
// I think it was necessary to have rewind working correctly
// after end of file. Since rewind is not used now anywhere,
// it should be ok not to reset position.
// This is necessary so that eof works properly for emacs files
// with no endline at end of file.
  }
  return *this;
}

void Vessel::parseFlag(const std::string&key, bool & t){
  // Check keyword has been registered
  plumed_massert(keywords.exists(key), "keyword " + key + " has not been registered");
  // Check keyword is a flag
  if(!keywords.style(key,"nohtml")){
     plumed_massert(keywords.style(key,"flag"), "keyword " + key + " is not a flag");
  }

  // Read in the flag otherwise get the default value from the keywords object
  if(!Tools::parseFlag(line,key,t)){
     if( keywords.style(key,"nohtml") ){ 
        t=false; 
     } else if ( !keywords.getLogicalDefault(key,t) ){
        plumed_merror("there is a flag with no logical default in a vessel - weird");
     }
  }
}

void Action::parseFlag(const std::string&key,bool & t) {
  // Check keyword has been registered
  plumed_massert(keywords.exists(key), "keyword " + key + " has not been registered");
  // Check keyword is a flag
  if(!keywords.style(key,"nohtml")) {
    plumed_massert( keywords.style(key,"vessel") || keywords.style(key,"flag") || keywords.style(key,"hidden"), "keyword " + key + " is not a flag");
  }

  // Read in the flag otherwise get the default value from the keywords object
  if(!Tools::parseFlag(line,key,t)) {
    if( keywords.style(key,"nohtml") || keywords.style(key,"vessel") ) {
      t=false;
    } else if ( !keywords.getLogicalDefault(key,t) ) {
      log.printf("ERROR in action %s with label %s : flag %s has no default",name.c_str(),label.c_str(),key.c_str() );
      plumed_error();
    }
  }
}

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
}