C++ EnergyDrift::current方法代码示例

//.........这里部分代码省略.........
        //case 5: // volume move
          //sys+=iso.move();
         // break;
        }

      if ( slp_global.runtest(0.00001) ) {
        xtc.setbox( nonbonded->geometry.len );
        xtc.save(textio::prefix+"traj.xtc", spc);
      }
      //if ( slp_global.runtest(0.1) ) 
        //(*mfc).sample(spc.p, )
#ifdef SLIT
      if ( slp_global.runtest(0.1) ) {
         for (auto &g : pol) {
           double d = gouy->dist2surf(g.cm);
           dst_map["Q"]( d )+=g.charge(spc.p);
           Point p=shape.vectorgyrationRadiusSquared(g, spc);
           dst_map["Rg2"]( d )+=p.x()+p.y()+p.z();
           dst_map["Rg2x"]( d )+=p.x();
           dst_map["Rg2z"]( d )+=p.z();
           dst_map["Ree2"]( d )+=spc.geo->sqdist( spc.p[g.front()], spc.p[g.back()] );
           //dst_map["<Energy>"]( d )+=sys.current();
           for (int i=g.front(); i<=g.back(); i++){
            ostringstream o;
            o << "qres" << i ;
            dst_map[o.str()]( d )+=spc.p[i].charge;
           }
         } 
      }
#endif

    } // end of micro loop

    temper.setCurrentEnergy( sys.current() );
    sys+=temper.move();

    sys.checkDrift( Energy::systemEnergy(spc,pot,spc.p) );

    mpi.cout << loop.timing();
    cout << loop.timing();

    rdf.save(textio::prefix+"rdf_p2p.dat");

#ifdef SLIT
    surfdist.save(textio::prefix+"surfdist.dat");
    surfmapall.save(textio::prefix+"surfall.dat");
    dst_map["Q"].save(textio::prefix+"netq.dat");
    //dst_map["<Energy>"].save(textio::prefix+"aveenergy.dat");

    for (auto &g : pol){   //Saving averages 
      //surfmap[g.name].save(textio::prefix+g.name+"surfdist.dat");
      std::ofstream f(textio::prefix+"res-surfdist.dat");
      std::ofstream f1(textio::prefix+"Rg-comp.dat");
      std::ofstream f2(textio::prefix+"qres.dat");
      f.precision(7);
      f1.precision(7);
      f2.precision(7);
      if (f && f1 && f2) {
        //f << "#g(z) of each residue from 1 to " << g.size() << endl;
        for (double d=0; d<=nonbonded->geometry.len.z(); d+=0.25){
          f1 << d; 
          if (dst_map["Rg2"](d).cnt > 0) f1 << "\t" << dst_map["Rg2"](d); //If counter is empty, then it prints shit lot of warning
          else f1 << "\t0";

          if (dst_map["Rg2x"](d).cnt > 0)f1 << "\t" << dst_map["Rg2x"](d); 
          else f1 << "\t0";

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

#ifdef TEMPER
  #define cout mpi.cout
  Faunus::MPI::MPIController mpi;
#endif

  cout << textio::splash();

  InputMap mcp(textio::prefix+"gouychapman.input");
  MCLoop loop(mcp);                    // class for handling mc loops
  FormatPQR pqr;                       // PQR structure file I/O
  FormatAAM aam;                       // AAM structure file I/O
  FormatXTC xtc(1000);                 // XTC gromacs trajectory format for only cuboid geomtries
  EnergyDrift sys;                     // class for tracking system energy drifts

  Tspace spc(mcp);

  auto pot = Energy::Nonbonded<Tspace, Tpairpot>(mcp) + 
    Energy::ExternalPotential<Tspace,Potential::GouyChapman<> >(mcp);
  pot.second.setSurfPositionZ( &spc.geo.len_half.z() );
  pot.setSpace(spc);

  // Add rigid molecules
  vector<GroupMolecular> pol( mcp.get("polymer_N",0));
  string polyfile = mcp.get<string>("polymer_file", "");
  for (auto &g : pol) {                    // load molecules
    aam.load(polyfile);
    Geometry::FindSpace f;
    f.find(*spc.geo, spc.p, aam.p);        // find empty spot in particle vector
    g = spc.insert( aam.p );               // insert into space
    g.name="Molecule";
    spc.enroll(g);
  }

  Move::TranslateRotate<Tspace> gmv(mcp,pot,spc);
#ifdef TEMPER
  Move::ParallelTempering<Tspace> temper(mcp,pot,spc,mpi);
#endif

  Analysis::RadialDistribution<float,int> rdf(0.25);
  Analysis::LineDistribution<float,int> surfdist(0.25);

  spc.load(textio::prefix+"state");
  pqr.save(textio::prefix+"initial.pqr", spc.p);

  sys.init( Energy::systemEnergy(spc,pot,spc.p) );

  cout << atom.info() << spc.info() << pot.info() << textio::header("MC Simulation Begins!");

  while ( loop.macroCnt() ) {                      // Markov chain 
    while ( loop.microCnt() ) {
      int k,i= slump.rand() % 1;
      switch (i) {
        case 0:                                    // translate and rotate molecules
          k=pol.size();
          while (k-->0) {
            gmv.setGroup( pol[ slump.rand() % pol.size() ] );
            sys+=gmv.move();
          }

          for (auto &g : pol)
            surfdist( gouy->dist2surf(g.cm) )++;   // molecule mass center to GC surface histogram

          for (auto i=pol.begin(); i!=pol.end()-1; i++)
            for (auto j=i+1; j!=pol.end(); j++)
              rdf( spc.geo->dist(i->cm,j->cm) )++; // molecule-molecule rdf

          break;
      }

      if ( slump.runtest(0.0001) ) {
        xtc.setbox( spc.geo.len );
        xtc.save(textio::prefix+"traj.xtc", spc);  // gromacs xtc file output
      }
    } // end of micro loop
#ifdef TEMPER
    temper.setCurrentEnergy( sys.current() );
    sys+=temper.move();
#endif

    sys.checkDrift( Energy::systemEnergy(spc,pot,spc.p) );

    cout << loop.timing();

    rdf.save(textio::prefix+"rdf_p2p.dat");
    surfdist.save(textio::prefix+"surfdist.dat");
    pqr.save(textio::prefix+"confout.pqr", spc.p);
    spc.save(textio::prefix+"state");

  } // end of macro loop

  cout << loop.info() << pot.first.info() << sys.info() << gmv.info();   
#ifdef TEMPER
  cout << temper.info();
#endif
}