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

int main(int argc, char* argv[]) {

  // findiff step: eps
  double eps=1.e-6;

  // various setups: com needs to be reset
  bool remove_com=true;
  // normalize weights or not (to check if the rmsd is proportional to weights) 
  bool normalize_weights=true;
  // use msd instead of rmsd to check consistency
  bool squared=false;
  // enhance com to enphasize the issues with COM treatment
  bool enhance_com=false;

  // parse input instructions
  // default task, calculate the RMSD of one frame respect to a set of others 
  vector<int> task(1);task[0]=0;

  // this test wants to be only for OPTIMAL case
  string type; type.assign("OPTIMAL"); 

  // first parse the task: in this applications the tasks are simple integers that pass the action that is required
  for(int i = 1; i < argc; i++){ 
      task.push_back(atoi(argv[i]));
  }  
  if(std::find(task.begin(), task.end(), -1)!=task.end()){cout<<"squared=true (default false)"<<endl;squared=true;}
  if(std::find(task.begin(), task.end(), -2)!=task.end()){cout<<"normalize_weights=false (default true)"<<endl;normalize_weights=false;}
  if(std::find(task.begin(), task.end(), -3)!=task.end()){cout<<"remove_com=false (default true)"<<endl; remove_com=false;}
  if(std::find(task.begin(), task.end(), -4)!=task.end()){cout<<"OPTIMAL-FAST (default OPTIMAL)"<<endl; type.assign("OPTIMAL-FAST");}
  if(std::find(task.begin(), task.end(), -5)!=task.end()){cout<<"enhance_com=true (default false) include option -3 (no com removal) "<<endl;enhance_com=true ; }
  if(enhance_com)remove_com=false;


  cout<<"ARGUMENTS: \n";
  cout<<"OPTIONS that go on top of tasks:\n";
  cout<<" -1 : squared=true (default=false)\n";  
  cout<<" -2 : normalize_weights=false (default=true)\n";
  cout<<" -3 : remove_com=false (default=true) \n";
  cout<<" -4 : OPTIMAL-FAST (default=OPTIMAL) \n";
  cout<<" -5 : enhance_com=true (default false) automatically set option -3 (i.e. check the trouble with com when you do not remove it)\n";
  cout<<"TASKS (can choose more than one):\n";
  cout<<"  0 : normal rmsd/msd calculation  and derivative dumps (default: always done)\n";
  cout<<"  1 : findiff test for  d msd / d position  (inhomogenehous weights)\n";
  cout<<"  2 : findiff test for  d msd / d reference (inhomogenehous weights)\n";
  cout<<"  3 : findiff test for  d msd / d position  (homogenehous weights)\n";
  cout<<"  4 : findiff test for  d msd / d reference (homogenehous weights)\n";
  cout<<"  5 : findiff test for  d Rot / d position  (inhomogenehous weights)  (reference->position)\n";
  cout<<"  6 : findiff test for  d Rot / d reference  (inhomogenehous weights) (reference->position)\n";
  cout<<"  7 : consistency check for MSD proportionality (works with squared=true through option -1 )\n";
  cout<<"  8 : do some timings for all the above routines and for a growing number of atoms\n";
  cout<<"  9 : test the rotation order: print position.pdb reference.pdb aligned.pdb and check that it makes sense(should be reference aligned onto positions)\n";
  cout<<" 10 : findiff test for  d Rot / d position  ( position -> reference ) \n";
  cout<<" 11 : findiff test for  d Rot / d position  (homogenehous weights)  (reference->position)\n";
  cout<<" 12 : findiff test for  d Rot / d reference  (homogenehous weights) (reference->position)\n";
  cout<<" 13 : do timings only for the most common use (aligment +derivatives) and repeat for homogeneous weights\n";
  cout<<" 20 : a simple test calculating only derivative of positions which is needed to compare new version and old version (need to compile with -DOLDRMSD both plumed and the test)\n";

  PDB pdbref;

  // now create the object: does not do anything but set the typer to SIMPLE
  PLMD::RMSD* rmsd=new RMSD();

  // set the reference pdb 
  string reference; reference.assign("1GB1_mdl1_rototranslated.pdb");
  PDB pdb;
  if( !pdb.read(reference,false,0.1) ) 
        cout<<"missing input file 1GB1_mdl1_rototranslated.pdb "<<"\n";
  if (enhance_com){
	vector<Vector> v=pdb.getPositions();
	for(unsigned i=0;i<v.size();i++){v[i][0]+=10.;v[i][1]+=20.;v[i][2]+=30.;}
	pdb.setPositions(v);
  }

  
  cout<<"NOW CREATING THE RMSD OBJECT... with set() method";  
  // remember that "set" method parses the reference, align and displace from the PDB by calling the following methods 
  //      setReference(pdb.getPositions());  -> set the reference and put the weights as 1/n, remove the com according to such weight 
  //      setAlign(pdb.getOccupancy());    -> normalizes, set the alignment vector and remove the Com using such weights 
  //      setDisplace(pdb.getBeta());    -> normalizes and set the displacement vector 
  //      setType(mytype);
  rmsd->set(pdb,type,remove_com,normalize_weights);
  // store other vectors from further manipulation
  std::vector<Vector> ref ;  ref=pdb.getPositions() ;
  std::vector<double> align; align=pdb.getOccupancy(); // non-normalized !
  std::vector<double> displace;  displace=pdb.getBeta(); // non-normalized ! 

  cout<<"DONE!"<<endl;  

  // now take another conformation to compare with: the running frame
  PDB pdbrun; 
  // mimic gromacs: do it in nm
  if( !pdbrun.read("1GB1_mdl2.pdb",false,0.1) )
        cout<<"missing input file 1GB1_mdl2.pdb\n" ;
  std::vector<Vector> run ;  run=pdbrun.getPositions() ;
  std::vector<Vector> derivatives ; 
  if (enhance_com){
        for(unsigned i=0;i<run.size();i++){run[i][0]-=10.;run[i][1]-=20.;run[i][2]-=30.;}
  }


//.........这里部分代码省略.........