C++ PDB类代码示例

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 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 RMSD::set(const PDB&pdb, string mytype ){

	setReference(pdb.getPositions());
	setAlign(pdb.getOccupancy());
	setDisplace(pdb.getBeta());
        setType(mytype);
}

int main()
{
	char szError[cbErrMax];
	EC ec;
	PDB* ppdb = PDB::Open("001.PDB", pdbWrite, 0, &ec, szError);

	assert(ppdb);
	DBI* pdbi = ppdb->OpenDBI(pdbWrite, "<target>");
	assert(pdbi);
	TS ts = pdbi->OpenTpi(pdbWrite);
	assert(ts);

	struct {
		unsigned short	len;
		unsigned short	leaf;
		unsigned long	signature;
		unsigned long	age;
		unsigned char	name[9];
	} typeUseTypeServer = {
		20,	
		LF_TYPESERVER,
		123456789,
		1,
		"\x07""002.PDB"
	};

	struct {
		unsigned short	reclen;	
		unsigned short	rectyp;	
		CV_uoff32_t 	off;
		unsigned short	seg;
		CV_typ_t		typind;	
		unsigned char	name[5];
	} aSymbol = {
		16,
		S_GDATA32,
		0,
		0,
		0,
		"\x03""ABC"
	};

	Mod* pmod = pdbi->OpenMod(ts, 1, "<mod>");
	assert(pmod->AddTypes((PB)&typeUseTypeServer, 20));
	for (TI ti = 0x1000; ti < 0x12e0; ti++) {
		aSymbol.typind = ti;
		assert(pmod->AddSymbols((PB)&aSymbol, 16));
	}
	pmod->Commit();
	pmod->Close();

	ts->Commit();
	ts->Close();
	ppdb->Commit();
	pdbi->Close();
	ppdb->Commit();
	ppdb->Close();

	return 0;
}

int main(int argc, char *argv[]) {
  string hdr = invocationHeader(argc, argv);
  opts::BasicOptions* bopts = new opts::BasicOptions;
  opts::ModelWithCoords* mopts = new opts::ModelWithCoords;
  ToolOptions* topts = new ToolOptions;
  opts::AggregateOptions options;
  options.add(bopts).add(mopts).add(topts);

  if (!options.parse(argc, argv))
    exit(-1);

  AtomicGroup model = mopts->model;
  
  if (topts->reimage) {
    if (!model.isPeriodic()) {
      cerr << "WARNING- Reimaging requested, but the model has no periodic box information\n";
    } else {
      if (!topts->bonds_name.empty()) {
        AtomicGroup bonds = createSystem(topts->bonds_name);
        copyBonds(model, bonds);
      }

      if (!model.hasBonds()) {
        cerr << "WARNING- The model has no connectivity.  Assigning bonds based on distance.\n";
        model.findBonds();
      }
    }
  }

  AtomicGroup center_mol = selectAtoms(model, topts->center_sel);
  GCoord center = center_mol.centroid();

  AtomicGroup apply_mol = selectAtoms(model, topts->apply_sel);
  GCoord offset = topts->translate - center;
  if (topts->center_xy)
    offset.z() = 0.0;
  
  for (AtomicGroup::iterator atom = apply_mol.begin(); atom != apply_mol.end(); ++atom)
    (*atom)->coords() += offset;

  if (topts->reimage && model.isPeriodic()) {
    vGroup molecules = model.splitByMolecule();
    vGroup segments = model.splitByUniqueSegid();
      
    for (vGroup::iterator seg = segments.begin(); seg != segments.end(); ++seg)
      seg->reimage();

    for (vGroup::iterator mol =molecules.begin(); mol != molecules.end(); ++mol)
      mol->reimage();
  }

  AtomicGroup write_mol = selectAtoms(model, topts->write_sel);
  PDB pdb = PDB::fromAtomicGroup(write_mol);
  pdb.remarks().add(hdr);
  cout << pdb;
}

void AdaptivePath::update() {
  double weight2 = -1.*mypathv->dx;
  double weight1 = 1.0 + mypathv->dx;
  if( weight1>1.0 ) {
    weight1=1.0; weight2=0.0;
  } else if( weight2>1.0 ) {
    weight1=0.0; weight2=1.0;
  }
  // Add projections to dispalcement accumulators
  ReferenceConfiguration* myref = getReferenceConfiguration( mypathv->iclose1 );
  myref->extractDisplacementVector( getPositions(), getArguments(), mypathv->cargs, false, displacement );
  getReferenceConfiguration( mypathv->iclose2 )->extractDisplacementVector( myref->getReferencePositions(), getArguments(), myref->getReferenceArguments(), false, displacement2 );
  displacement.addDirection( -mypathv->dx, displacement2 );
  pdisplacements[mypathv->iclose1].addDirection( weight1, displacement );
  pdisplacements[mypathv->iclose2].addDirection( weight2, displacement );
  // Update weight accumulators
  wsum[mypathv->iclose1] *= fadefact;
  wsum[mypathv->iclose2] *= fadefact;
  wsum[mypathv->iclose1] += weight1;
  wsum[mypathv->iclose2] += weight2;

  // This does the update of the path if it is time to
  if( (getStep()>0) && (getStep()%update_str==0) ) {
    wsum[fixedn[0]]=wsum[fixedn[1]]=0.;
    for(unsigned inode=0; inode<getNumberOfReferencePoints(); ++inode) {
      if( wsum[inode]>0 ) {
        // First displace the node by the weighted direction
        getReferenceConfiguration( inode )->displaceReferenceConfiguration( 1./wsum[inode], pdisplacements[inode] );
        // Reset the displacement
        pdisplacements[inode].zeroDirection();
      }
    }
    // Now ensure all the nodes of the path are equally spaced
    PathReparameterization myspacings( getPbc(), getArguments(), getAllReferenceConfigurations() );
    myspacings.reparameterize( fixedn[0], fixedn[1], tolerance );
  }
  if( (getStep()>0) && (getStep()%wstride==0) ) {
    pathfile.printf("# PATH AT STEP %d TIME %f \n", getStep(), getTime() );
    std::vector<std::unique_ptr<ReferenceConfiguration>>& myconfs=getAllReferenceConfigurations();
    std::vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
    if( moldat.size()>1 ) error("you should only have one MOLINFO action in your input file");
    SetupMolInfo* mymoldat=NULL; if( moldat.size()==1 ) mymoldat=moldat[0];
    std::vector<std::string> argument_names( getNumberOfArguments() );
    for(unsigned i=0; i<getNumberOfArguments(); ++i) argument_names[i] = getPntrToArgument(i)->getName();
    PDB mypdb; mypdb.setArgumentNames( argument_names );
    for(unsigned i=0; i<myconfs.size(); ++i) {
      pathfile.printf("REMARK TYPE=%s\n", myconfs[i]->getName().c_str() );
      mypdb.setAtomPositions( myconfs[i]->getReferencePositions() );
      for(unsigned j=0; j<getNumberOfArguments(); ++j) mypdb.setArgumentValue( getPntrToArgument(j)->getName(), myconfs[i]->getReferenceArgument(j) );
      mypdb.print( atoms.getUnits().getLength()/0.1, mymoldat, pathfile, ofmt );
    }
    pathfile.flush();
  }
}

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

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 CCPDB::parseFrame(void) {
    if (ifs->eof())
      return(false);

    // We cheat here...  Maybe it would be better to have a PDB::clear()
    // or AtomicGroup::clear() member function???
    // Note:  For some reason, PDB newframe(*(ifs)) doesn't parse correctly...
    PDB newframe;
    newframe.read(*(ifs));
    frame = newframe;
    if (frame.size() == 0)
      return(false);

    return(true);
  }

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

TrigonometricPathVessel::TrigonometricPathVessel( const vesselbase::VesselOptions& da ):
  StoreDataVessel(da),
  projdir(ReferenceConfigurationOptions("DIRECTION")),
  mydpack1( 1, getAction()->getNumberOfDerivatives() ),
  mydpack2( 1, getAction()->getNumberOfDerivatives() ),
  mydpack3( 1, getAction()->getNumberOfDerivatives() ),
  mypack1( 0, 0, mydpack1 ),
  mypack2( 0, 0, mydpack2 ),
  mypack3( 0, 0, mydpack3 )
{
  mymap=dynamic_cast<Mapping*>( getAction() );
  plumed_massert( mymap, "Trigonometric path vessel can only be used with mappings");
  // Retrieve the index of the property in the underlying mapping
  if( mymap->getNumberOfProperties()!=1 ) error("cannot use trigonometric paths when there are multiple properties");

  for(unsigned i=0; i<mymap->getFullNumberOfTasks(); ++i) {
    if( mymap->getTaskCode(i)!=mymap->getPositionInFullTaskList(i) ) error("mismatched tasks and codes");
  }
  mymap->addComponentWithDerivatives("gspath"); mymap->componentIsNotPeriodic("gspath");
  sp=mymap->copyOutput( mymap->getNumberOfComponents()-1 ); sp->resizeDerivatives( mymap->getNumberOfDerivatives() );
  mymap->addComponentWithDerivatives("gzpath"); mymap->componentIsNotPeriodic("gzpath");
  zp=mymap->copyOutput( mymap->getNumberOfComponents()-1 ); zp->resizeDerivatives( mymap->getNumberOfDerivatives() );

  // Check we have PCA
  ReferenceConfiguration* ref0=mymap->getReferenceConfiguration(0);
  for(unsigned i=0; i<mymap->getFullNumberOfTasks(); ++i) {
    if( !(mymap->getReferenceConfiguration(i))->pcaIsEnabledForThisReference() ) error("pca must be implemented in order to use trigometric path");
    if( ref0->getName()!=(mymap->getReferenceConfiguration(i))->getName() ) error("cannot use mixed metrics");
    if( mymap->getNumberOfAtoms()!=(mymap->getReferenceConfiguration(i))->getNumberOfReferencePositions() ) error("all frames must use the same set of atoms");
    if( mymap->getNumberOfArguments()!=(mymap->getReferenceConfiguration(i))->getNumberOfReferenceArguments() ) error("all frames must use the same set of arguments");
  }

  cargs.resize( mymap->getNumberOfArguments() ); std::vector<std::string> argument_names( mymap->getNumberOfArguments() );
  for(unsigned i=0; i<mymap->getNumberOfArguments(); ++i) argument_names[i] = (mymap->getPntrToArgument(i))->getName();
  PDB mypdb; mypdb.setAtomNumbers( mymap->getAbsoluteIndexes() ); mypdb.addBlockEnd( mymap->getAbsoluteIndexes().size() );
  if( argument_names.size()>0 ) mypdb.setArgumentNames( argument_names );
  projdir.read( mypdb );
  mypack1.resize( mymap->getNumberOfArguments(), mymap->getNumberOfAtoms() ); ref0->setupPCAStorage( mypack1 );
  mypack2.resize( mymap->getNumberOfArguments(), mymap->getNumberOfAtoms() ); ref0->setupPCAStorage( mypack2 );
  mypack3.resize( mymap->getNumberOfArguments(), mymap->getNumberOfAtoms() );
  for(unsigned i=0; i<mymap->getNumberOfAtoms(); ++i) { mypack1.setAtomIndex(i,i); mypack2.setAtomIndex(i,i); mypack3.setAtomIndex(i,i); }
  mypack1_stashd_atoms.resize( mymap->getNumberOfAtoms() ); mypack1_stashd_args.resize( mymap->getNumberOfArguments() );
}

int colvarproxy_namd::load_atoms(char const *pdb_filename,
                                 cvm::atom_group &atoms,
                                 std::string const &pdb_field_str,
                                 double const pdb_field_value)
{
  if (pdb_field_str.size() == 0)
    cvm::error("Error: must define which PDB field to use "
               "in order to define atoms from a PDB file.\n", INPUT_ERROR);

  PDB *pdb = new PDB(pdb_filename);
  size_t const pdb_natoms = pdb->num_atoms();

  e_pdb_field pdb_field_index = pdb_field_str2enum(pdb_field_str);

  for (size_t ipdb = 0; ipdb < pdb_natoms; ipdb++) {

    double atom_pdb_field_value = 0.0;

    switch (pdb_field_index) {
    case e_pdb_occ:
      atom_pdb_field_value = (pdb->atom(ipdb))->occupancy();
      break;
    case e_pdb_beta:
      atom_pdb_field_value = (pdb->atom(ipdb))->temperaturefactor();
      break;
    case e_pdb_x:
      atom_pdb_field_value = (pdb->atom(ipdb))->xcoor();
      break;
    case e_pdb_y:
      atom_pdb_field_value = (pdb->atom(ipdb))->ycoor();
      break;
    case e_pdb_z:
      atom_pdb_field_value = (pdb->atom(ipdb))->zcoor();
      break;
    default:
      break;
    }

    if ( (pdb_field_value) &&
         (atom_pdb_field_value != pdb_field_value) ) {
      continue;
    } else if (atom_pdb_field_value == 0.0) {
      continue;
    }

    if (atoms.is_enabled(colvardeps::f_ag_scalable)) {
      atoms.add_atom_id(ipdb);
    } else {
      atoms.add_atom(cvm::atom(ipdb+1));
    }
  }

  delete pdb;
  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}

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

int colvarproxy_namd::load_coords(char const *pdb_filename,
                                  std::vector<cvm::atom_pos> &pos,
                                  const std::vector<int> &indices,
                                  std::string const &pdb_field_str,
                                  double const pdb_field_value)
{
  if (pdb_field_str.size() == 0 && indices.size() == 0) {
    cvm::error("Bug alert: either PDB field should be defined or list of "
               "atom IDs should be available when loading atom coordinates!\n", BUG_ERROR);
  }

  e_pdb_field pdb_field_index;
  bool const use_pdb_field = (pdb_field_str.size() > 0);
  if (use_pdb_field) {
    pdb_field_index = pdb_field_str2enum(pdb_field_str);
  }

  // next index to be looked up in PDB file (if list is supplied)
  std::vector<int>::const_iterator current_index = indices.begin();

  PDB *pdb = new PDB(pdb_filename);
  size_t const pdb_natoms = pdb->num_atoms();

  if (pos.size() != pdb_natoms) {

    bool const pos_allocated = (pos.size() > 0);

    size_t ipos = 0, ipdb = 0;
    for ( ; ipdb < pdb_natoms; ipdb++) {

      if (use_pdb_field) {
        // PDB field mode: skip atoms with wrong value in PDB field
        double atom_pdb_field_value = 0.0;

        switch (pdb_field_index) {
        case e_pdb_occ:
          atom_pdb_field_value = (pdb->atom(ipdb))->occupancy();
          break;
        case e_pdb_beta:
          atom_pdb_field_value = (pdb->atom(ipdb))->temperaturefactor();
          break;
        case e_pdb_x:
          atom_pdb_field_value = (pdb->atom(ipdb))->xcoor();
          break;
        case e_pdb_y:
          atom_pdb_field_value = (pdb->atom(ipdb))->ycoor();
          break;
        case e_pdb_z:
          atom_pdb_field_value = (pdb->atom(ipdb))->zcoor();
          break;
        default:
          break;
        }

        if ( (pdb_field_value) &&
             (atom_pdb_field_value != pdb_field_value) ) {
          continue;
        } else if (atom_pdb_field_value == 0.0) {
          continue;
        }

      } else {
        // Atom ID mode: use predefined atom IDs from the atom group
        if (((int) ipdb) != *current_index) {
          // Skip atoms not in the list
          continue;
        } else {
          current_index++;
        }
      }

      if (!pos_allocated) {
        pos.push_back(cvm::atom_pos(0.0, 0.0, 0.0));
      } else if (ipos >= pos.size()) {
        cvm::error("Error: the PDB file \""+
                   std::string(pdb_filename)+
                   "\" contains coordinates for "
                   "more atoms than needed.\n", BUG_ERROR);
      }

      pos[ipos] = cvm::atom_pos((pdb->atom(ipdb))->xcoor(),
                                (pdb->atom(ipdb))->ycoor(),
                                (pdb->atom(ipdb))->zcoor());
      ipos++;
      if (!use_pdb_field && current_index == indices.end())
        break;
    }

    if (ipos < pos.size() || (!use_pdb_field && current_index != indices.end())) {
      size_t n_requested = use_pdb_field ? pos.size() : indices.size();
      cvm::error("Error: number of matching records in the PDB file \""+
                 std::string(pdb_filename)+"\" ("+cvm::to_str(ipos)+
                 ") does not match the number of requested coordinates ("+
                 cvm::to_str(n_requested)+").\n", INPUT_ERROR);
      return COLVARS_ERROR;
    }
  } else {

    // when the PDB contains exactly the number of atoms of the array,
    // ignore the fields and just read coordinates
//.........这里部分代码省略.........