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

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

PathMSDBase::PathMSDBase(const ActionOptions&ao):
PLUMED_COLVAR_INIT(ao),
neigh_size(-1),
neigh_stride(-1),
nframes(0)
{
  parse("LAMBDA",lambda);
  parse("NEIGH_SIZE",neigh_size);
  parse("NEIGH_STRIDE",neigh_stride);
  parse("REFERENCE",reference);

  // open the file
  FILE* fp=fopen(reference.c_str(),"r");
  std::vector<AtomNumber> aaa;
  if (fp!=NULL)
  {
    log<<"Opening reference file "<<reference.c_str()<<"\n";
    bool do_read=true;
    while (do_read){
         PDB mypdb; 
         RMSD mymsd(log); 
         do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
         if(do_read){
            unsigned nat=0;
            nframes++;
            if(mypdb.getAtomNumbers().size()==0) error("number of atoms in a frame should be more than zero");
            if(nat==0) nat=mypdb.getAtomNumbers().size();
            if(nat!=mypdb.getAtomNumbers().size()) error("frames should have the same number of atoms");
            if(aaa.empty()) aaa=mypdb.getAtomNumbers();
            if(aaa!=mypdb.getAtomNumbers()) error("frames should contain same atoms in same order");
            log<<"Found PDB: "<<nframes<<" containing  "<<mypdb.getAtomNumbers().size()<<" atoms\n"; 
	    pdbv.push_back(mypdb); 
//            requestAtoms(mypdb.getAtomNumbers()); // is done in non base classes 
            derivs_s.resize(mypdb.getAtomNumbers().size());
            derivs_z.resize(mypdb.getAtomNumbers().size());
            mymsd.set(mypdb,"OPTIMAL");
            msdv.push_back(mymsd); // the vector that stores the frames
            //log<<mypdb; 
         }else{break ;}
    }
    fclose (fp);
    log<<"Found TOTAL "<<nframes<< " PDB in the file "<<reference.c_str()<<" \n"; 
    if(nframes==0) error("at least one frame expected");
  } 
  if(neigh_stride>0 || neigh_size>0){
           if(neigh_size>int(nframes)){
           	log.printf(" List size required ( %d ) is too large: resizing to the maximum number of frames required: %u  \n",neigh_size,nframes);
 		neigh_size=nframes;
           }
           log.printf("  Neighbor list enabled: \n");
           log.printf("                size   :  %d elements\n",neigh_size);
           log.printf("                stride :  %d timesteps \n",neigh_stride);
  }else{
           log.printf("  Neighbor list NOT enabled \n");
  }

}

bool DataCollectionObject::transferDataToPDB( PDB& mypdb ) {
  // Check if PDB contains argument names
  std::vector<std::string> pdb_args( mypdb.getArgumentNames() );
  // Now set the argument values
  std::map<std::string,double>::iterator it;
  for(unsigned i=0; i<pdb_args.size(); ++i) {
    it=args.find( pdb_args[i] );
    if( it==args.end() ) return false;
    mypdb.setArgumentValue( pdb_args[i], it->second );
  }
  // Now set the atomic positions
  std::vector<AtomNumber> pdb_pos( mypdb.getAtomNumbers() );
  if( pdb_pos.size()==positions.size() ) mypdb.setAtomPositions( positions );
  else if( pdb_pos.size()>0 ) plumed_merror("This feature is currently not ready");
  return true;
}

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

void MultiDomainRMSD::read( const PDB& pdb ){
   unsigned nblocks =  pdb.getNumberOfAtomBlocks();
   if( nblocks<2 ) error("multidomain RMSD only has one block of atoms");
  
   std::vector<AtomNumber> atomnumbers;
   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 lower=0.0, upper=std::numeric_limits<double>::max( );
   parse("LOWER_CUTOFF",lower,true); 
   parse("UPPER_CUTOFF",upper,true);

   for(unsigned i=1;i<=nblocks;++i){
       Tools::convert(i,num);
       if( ftype=="RMSD" ){
          parse("TYPE"+num, ftype );
          parse("LOWER_CUTOFF"+num,lower,true); 
          parse("UPPER_CUTOFF"+num,upper,true); 
       }
       domains.push_back( metricRegister().create<SingleDomainRMSD>( ftype ) );
       positions.resize( blocks[i] - blocks[i-1] + 1 );
       align.resize( blocks[i] - blocks[i-1] + 1 );
       displace.resize( blocks[i] - blocks[i-1] + 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( true, lower, upper );  // Currently no option for nopbc
       domains[i-1]->setReferenceAtoms( positions, align, displace );
       domains[i-1]->setNumberOfAtoms( positions.size() );
       
       double ww=0; parse("WEIGHT"+num, ww, true );
       if( ww==0 ) weights.push_back( 1.0 );
       else weights.push_back( ww );
   }   
   // And set the atom numbers for this object
   setAtomNumbers( pdb.getAtomNumbers() );
}

//.........这里部分代码省略.........
      	sw2.start();	
        for(unsigned int j=0;j<ntest;j++)rmsd->calc_DDistDRef( r_run, derivatives,DDistDRef, squared); 
      	sw2.stop();	
      	cout<<"WITH REFERENCE FRAME: \n"<<sw2<<endl;

      	Tensor Rotation;
      	Matrix<std::vector<Vector> > DRotDPos(3,3);	
      	Stopwatch sw3;
      	sw3.start();	
        for(unsigned int j=0;j<ntest;j++)rmsd->calc_DDistDRef_Rot_DRotDPos( r_run, derivatives,DDistDRef, Rotation , DRotDPos, squared); 
      	sw3.stop();	
      	cout<<"WITH ROTATION MATRIX DERIVATIVE: \n"<<sw3<<endl;

        Matrix<std::vector<Vector> > DRotDRef(3,3);
      	Stopwatch sw4;
      	sw4.start();	
        for(unsigned int j=0;j<ntest;j++)rmsd->calc_DDistDRef_Rot_DRotDPos_DRotDRef( r_run, derivatives ,DDistDRef, Rotation , DRotDPos, DRotDRef, squared); 
      	sw4.stop();	
      	cout<<"WITH ROTATION MATRIX DERIVATIVE OF REEFERENCE: \n"<<sw4<<endl;
      	// duplicate the atoms
      	unsigned s=r_run.size();
      	for (unsigned int i=0;i<s;i++){r_run.push_back(r_run[i]);r_ref.push_back(r_ref[i]);r_al.push_back(r_al[i]);r_disp.push_back(r_disp[i]);}
      
      }
  } 
  // Task 9: check the rotation
  if(std::find(task.begin(), task.end(), 9)!=task.end()){
      cout<<"Task 9: dump some pdbs so to check if all makes sense when using inverse transform. In particular positions_aligned.pdb should overlap with reference_centered.pdb "<<endl;
	rmsd->clear();
        rmsd->set(align, displace, ref,type,remove_com,normalize_weights  );	
	// dump the reference
	ofstream myfile;
	myfile.open ("reference.pdb");		
	std::vector<AtomNumber> at=pdb.getAtomNumbers();
	std::vector<Vector>   pos=pdb.getPositions();
	unsigned k=0;
	for(std::vector<AtomNumber>::iterator i=at.begin(); i!=at.end(); i++){
		myfile<<"ATOM";
                myfile.width(7);myfile<<std::right<<(*i).serial()<<" "; 
                myfile.width(4);myfile<<std::left<<pdb.getAtomName(*i); 
                myfile.width(4);myfile<<std::right<<pdb.getResidueName(*i)<<" A"; 
                myfile.width(4);myfile<<std::right<<pdb.getResidueNumber(*i)<<"    "; 
		myfile.setf( std::ios::fixed, std:: ios::floatfield );
                myfile.width(8);myfile.precision(3); myfile<<std::right<<pos[k][0]*10; 
                myfile.width(8);myfile.precision(3); myfile<<std::right<<pos[k][1]*10; 
                myfile.width(8);myfile.precision(3); myfile<<std::right<<pos[k][2]*10<<"  1.00  1.00\n"; 
		k++;	
	}	
	myfile.close();			
	// dump the position
	myfile.open ("positions.pdb");		
	at=pdbrun.getAtomNumbers();
	std::vector<Vector>   runpos=pdbrun.getPositions();
	k=0;
	for(std::vector<AtomNumber>::iterator i=at.begin(); i!=at.end(); i++){
		myfile<<"ATOM";
                myfile.width(7);myfile<<std::right<<(*i).serial()<<" "; 
                myfile.width(4);myfile<<std::left<<pdbrun.getAtomName(*i); 
                myfile.width(4);myfile<<std::right<<pdbrun.getResidueName(*i)<<" A"; 
                myfile.width(4);myfile<<std::right<<pdbrun.getResidueNumber(*i)<<"    "; 
		myfile.setf( std::ios::fixed, std:: ios::floatfield );
                myfile.width(8);myfile.precision(3); myfile<<std::right<<runpos[k][0]*10; 
                myfile.width(8);myfile.precision(3); myfile<<std::right<<runpos[k][1]*10; 
                myfile.width(8);myfile.precision(3); myfile<<std::right<<runpos[k][2]*10<<"  1.00  1.00\n"; 
		k++;	
	}	

//.........这里部分代码省略.........
// variables to test particle decomposition
  int pd_nlocal;
  int pd_start;
// variables to test random decomposition (=domain decomposition)
  std::vector<int>  dd_gatindex;
  std::vector<int>  dd_g2l;
  std::vector<real> dd_masses;
  std::vector<real> dd_charges;
  std::vector<real> dd_forces;
  std::vector<real> dd_coordinates;
  int dd_nlocal;
// random stream to choose decompositions
  Random rnd;

  while(true){
    if(!noatoms){
       if(!Tools::getline(fp,line)) break;
    }

    int natoms;
    bool first_step=false;
    if(!noatoms){
      if(trajectory_fmt=="gro") if(!Tools::getline(fp,line)) error("premature end of trajectory file");
      sscanf(line.c_str(),"%100d",&natoms);
    }
    if(checknatoms<0 && !noatoms){
      pd_nlocal=natoms;
      pd_start=0;
      first_step=true;
      masses.assign(natoms,real(1.0));
      charges.assign(natoms,real(0.0));
      if(pdbfile.length()>0){
        for(unsigned i=0;i<pdb.size();++i){
          AtomNumber an=pdb.getAtomNumbers()[i];
          unsigned index=an.index();
          if( index>=unsigned(natoms) ) error("atom index in pdb exceeds the number of atoms in trajectory");
          masses[index]=pdb.getOccupancy()[i];
          charges[index]=pdb.getBeta()[i];
        }
      }
    } else if( checknatoms<0 && noatoms ){ 
      natoms=0; 
    }
    if( checknatoms<0 ){
      checknatoms=natoms;
      p.cmd("setNatoms",&natoms);
      p.cmd("init");
    }
    if(checknatoms!=natoms){
       std::string stepstr; Tools::convert(step,stepstr);
       error("number of atoms in frame " + stepstr + " does not match number of atoms in first frame");
    }

    coordinates.assign(3*natoms,real(0.0));
    forces.assign(3*natoms,real(0.0));
    cell.assign(9,real(0.0));
    virial.assign(9,real(0.0));

    if( first_step || rnd.U01()>0.5){
      if(debug_pd){
        int npe=intracomm.Get_size();
        vector<int> loc(npe,0);
        vector<int> start(npe,0);
        for(int i=0;i<npe-1;i++){
          int cc=(natoms*2*rnd.U01())/npe;
          if(start[i]+cc>natoms) cc=natoms-start[i];