C++ PDB::size方法代码示例

void SingleDomainRMSD::readReference( const PDB& pdb ){
  readAtomsFromPDB( pdb );
  double wa=0, wd=0;
  for(unsigned i=0;i<pdb.size();++i){ wa+=align[i]; wd+=displace[i]; }

  Vector center;
  for(unsigned i=0;i<pdb.size();++i){
     align[i]=align[i] / wa; displace[i]=displace[i] / wd;
     center+=reference_atoms[i]*align[i];
  }
  for(unsigned i=0;i<pdb.size();++i) reference_atoms[i]-=center;
} 

void ReferenceAtoms::readAtomsFromPDB( const PDB& pdb, const bool allowblocks  ){
  if( !allowblocks && pdb.getNumberOfAtomBlocks()!=1 ) error("found multi-atom-block pdb format but expecting only one block of atoms");  

  for(unsigned i=0;i<pdb.size();++i){
     indices.push_back( pdb.getAtomNumbers()[i] ); reference_atoms.push_back( pdb.getPositions()[i] );
     align.push_back( pdb.getOccupancy()[i] ); displace.push_back( pdb.getBeta()[i] );
  }
  atom_der_index.resize( reference_atoms.size() );
}

void MultiDomainRMSD::read( const PDB& pdb ) {
  unsigned nblocks =  pdb.getNumberOfAtomBlocks();
  if( nblocks<2 ) error("multidomain RMSD only has one block of atoms");

  std::vector<Vector> positions; std::vector<double> align, displace;
  std::string num; blocks.resize( nblocks+1 ); blocks[0]=0;
  for(unsigned i=0; i<nblocks; ++i) blocks[i+1]=pdb.getAtomBlockEnds()[i];

  double tmp, lower=0.0, upper=std::numeric_limits<double>::max( );
  if( pdb.getArgumentValue("LOWER_CUTOFF",tmp) ) lower=tmp;
  if( pdb.getArgumentValue("UPPER_CUTOFF",tmp) ) upper=tmp;
  bool nopbc=pdb.hasFlag("NOPBC");

  domains.resize(0); weights.resize(0);
  for(unsigned i=1; i<=nblocks; ++i) {
    Tools::convert(i,num);
    if( ftype=="RMSD" ) {
      // parse("TYPE"+num, ftype );
      lower=0.0; upper=std::numeric_limits<double>::max( );
      if( pdb.getArgumentValue("LOWER_CUTOFF"+num,tmp) ) lower=tmp;
      if( pdb.getArgumentValue("UPPER_CUTOFF"+num,tmp) ) upper=tmp;
      nopbc=pdb.hasFlag("NOPBC");
    }
    domains.emplace_back( metricRegister().create<SingleDomainRMSD>( ftype ) );
    positions.resize( blocks[i] - blocks[i-1] );
    align.resize( blocks[i] - blocks[i-1] );
    displace.resize( blocks[i] - blocks[i-1] );
    unsigned n=0;
    for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) {
      positions[n]=pdb.getPositions()[j];
      align[n]=pdb.getOccupancy()[j];
      displace[n]=pdb.getBeta()[j];
      n++;
    }
    domains[i-1]->setBoundsOnDistances( !nopbc, lower, upper );
    domains[i-1]->setReferenceAtoms( positions, align, displace );
    domains[i-1]->setupRMSDObject();

    double ww=0;
    if( !pdb.getArgumentValue("WEIGHT"+num,ww) ) weights.push_back( 1.0 );
    else weights.push_back( ww );
  }
  // And set the atom numbers for this object
  indices.resize(0); atom_der_index.resize(0);
  for(unsigned i=0; i<pdb.size(); ++i) { indices.push_back( pdb.getAtomNumbers()[i] ); atom_der_index.push_back(i); }
  // setAtomNumbers( pdb.getAtomNumbers() );
}

int Driver<real>::main(FILE* in,FILE*out,Communicator& pc){

  Units units;
  PDB pdb;

// Parse everything
  bool printhelpdebug; parseFlag("--help-debug",printhelpdebug);
  if( printhelpdebug ){
      fprintf(out,"%s",
         "Additional options for debug (only to be used in regtest):\n"
         "  [--debug-float]         : turns on the single precision version (to check float interface)\n"
         "  [--debug-dd]            : use a fake domain decomposition\n"
         "  [--debug-pd]            : use a fake particle decomposition\n"
      );
      return 0;
  }
  // Are we reading trajectory data
  bool noatoms; parseFlag("--noatoms",noatoms);

  std::string fakein; 
  bool debugfloat=parse("--debug-float",fakein);
  if(debugfloat && sizeof(real)!=sizeof(float)){
      CLTool* cl=cltoolRegister().create(CLToolOptions("driver-float"));    //new Driver<float>(*this);
      cl->setInputData(this->getInputData());
      int ret=cl->main(in,out,pc);
      delete cl;
      return ret;
  }

  bool debug_pd=parse("--debug-pd",fakein);
  bool debug_dd=parse("--debug-dd",fakein);
  if( debug_pd || debug_dd ){
    if(noatoms) error("cannot debug without atoms");
  }

// set up for multi replica driver:
  int multi=0;
  parse("--multi",multi);
  Communicator intracomm;
  Communicator intercomm;
  if(multi){
    int ntot=pc.Get_size();
    int nintra=ntot/multi;
    if(multi*nintra!=ntot) error("invalid number of processes for multi environment");
    pc.Split(pc.Get_rank()/nintra,pc.Get_rank(),intracomm);
    pc.Split(pc.Get_rank()%nintra,pc.Get_rank(),intercomm);
  } else {
    intracomm.Set_comm(pc.Get_comm());
  }

// set up for debug replica exchange:
  bool debug_grex=parse("--debug-grex",fakein);
  int  grex_stride=0;
  FILE*grex_log=NULL;
  if(debug_grex){
    if(noatoms) error("must have atoms to debug_grex");
    if(multi<2)  error("--debug_grex needs --multi with at least two replicas");
    Tools::convert(fakein,grex_stride);
    string n; Tools::convert(intercomm.Get_rank(),n);
    string file;
    parse("--debug-grex-log",file);
    if(file.length()>0){
      file+="."+n;
      grex_log=fopen(file.c_str(),"w");
    }
  }

// Read the plumed input file name  
  string plumedFile; parse("--plumed",plumedFile);
// the timestep
  double t; parse("--timestep",t);
  real timestep=real(t);
// the stride
  unsigned stride; parse("--trajectory-stride",stride);
// are we writing forces
  string dumpforces(""), dumpforcesFmt("%f");; 
  if(!noatoms) parse("--dump-forces",dumpforces);
  if(dumpforces!="") parse("--dump-forces-fmt",dumpforcesFmt);

  string trajectory_fmt;

// Read in an xyz file
  string trajectoryFile(""), pdbfile("");
  bool pbc_cli_given=false; vector<double> pbc_cli_box(9,0.0);
  if(!noatoms){
     std::string traj_xyz; parse("--ixyz",traj_xyz);
     std::string traj_gro; parse("--igro",traj_gro);
     if(traj_xyz.length()>0 && traj_gro.length()>0){
       fprintf(stderr,"ERROR: cannot provide more than one trajectory file\n");
       if(grex_log)fclose(grex_log);
       return 1;
     }
     if(traj_xyz.length()>0 && trajectoryFile.length()==0){
       trajectoryFile=traj_xyz;
       trajectory_fmt="xyz";
     }
     if(traj_gro.length()>0 && trajectoryFile.length()==0){
       trajectoryFile=traj_gro;
       trajectory_fmt="gro";
     }
//.........这里部分代码省略.........