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

int main() {

  InputMap mcp( "gcmol.json" );

  Tspace spc( mcp );
  spc.load( "state", Tspace::RESIZE );

  auto pot = Energy::Nonbonded<Tspace,Tpairpot>(mcp);
  Move::Propagator<Tspace> mv( mcp, pot, spc );

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

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

  MCLoop loop( mcp );
  while ( loop[0] ) {
    while ( loop[1] )
      sys += mv.move();
    cout << loop.timing();
  }

  sys.checkDrift( Energy::systemEnergy( spc,pot,spc.p ) );
  spc.save("state");
  FormatPQR::save("confout.pqr", spc.p, spc.geo.len);

  UnitTest test( mcp );
  sys.test( test );
  mv.test( test );

  cout << loop.info() << sys.info() << spc.info() << mv.info() << test.info() << endl;

  return test.numFailed();
}

int main() {
    ::atom.includefile("stockmayer.json");         // load atom properties
    InputMap in("stockmayer.input");               // open parameter file for user input
    Energy::NonbondedVector<Tspace,Tpair> pot(in); // non-bonded only
    EnergyDrift sys;                               // class for tracking system energy drifts
    Tspace spc(in);                // create simulation space, particles etc.
    Group sol;
    sol.addParticles(spc, in);                     // group for particles
    MCLoop loop(in);                               // class for mc loop counting
    Analysis::RadialDistribution<> rdf(0.05);       // particle-particle g(r)
    Analysis::Table2D<double,Average<double> > mucorr(0.1);       // particle-particle g(r)
    TmoveTran trans(in,pot,spc);
    TmoveRot rot(in,pot,spc);
    trans.setGroup(sol);                                // tells move class to act on sol group
    rot.setGroup(sol);                                  // tells move class to act on sol group
    spc.load("state");
    spc.p = spc.trial;
    UnitTest test(in);               // class for unit testing
    FormatXTC xtc(spc.geo.len.norm());

    sys.init( Energy::systemEnergy(spc,pot,spc.p)  );   // initial energy
    while ( loop.macroCnt() ) {                         // Markov chain
        while ( loop.microCnt() ) {
            if (slp_global() > 0.5)
                sys+=trans.move( sol.size() );                // translate
            else
                sys+=rot.move( sol.size() );                  // rotate

            if (slp_global()<1.5)
                for (auto i=sol.front(); i<sol.back(); i++) { // salt rdf
                    for (auto j=i+1; j<=sol.back(); j++) {
                        double r =spc.geo.dist(spc.p[i],spc.p[j]);
                        rdf(r)++;
                        mucorr(r) += spc.p[i].mu.dot(spc.p[j].mu);
                    }
                }
            if (slp_global()>0.99)
                xtc.save(textio::prefix+"out.xtc", spc.p);
        }
        sys.checkDrift(Energy::systemEnergy(spc,pot,spc.p)); // compare energy sum with current
        cout << loop.timing();
    }

    // perform unit tests
    trans.test(test);
    rot.test(test);
    sys.test(test);

    FormatPQR().save("confout.pqr", spc.p);
    rdf.save("gofr.dat");                               // save g(r) to disk
    mucorr.save("mucorr.dat");                               // save g(r) to disk
    std::cout << spc.info() + pot.info() + trans.info()
              + rot.info() + sys.info() + test.info(); // final info
    spc.save("state");

    return test.numFailed();
}

int main() {
  InputMap mcp("grand.input");                  // open user input file
  Tspace spc(mcp);                              // simulation space
  Energy::Nonbonded<Tspace,CoulombHS> pot(mcp); // hamiltonian
  pot.setSpace(spc);                            // share space w. hamiltonian

  Group salt;                                   // group for salt particles
  salt.addParticles(spc, mcp);
  spc.load("state",Tspace::RESIZE);             // load old config. from disk (if any)

  // Two different Widom analysis methods
  double lB = Coulomb(mcp).bjerrumLength();     // get bjerrum length
  Analysis::Widom<PointParticle> widom1;        // widom analysis (I)
  Analysis::WidomScaled<Tspace> widom2(lB,1);   // ...and (II)
  widom1.add(spc.p);
  widom2.add(spc.p);

  Move::GrandCanonicalSalt<Tspace> gc(mcp,pot,spc,salt);
  Move::AtomicTranslation<Tspace> mv(mcp,pot,spc);
  mv.setGroup(salt);

  EnergyDrift sys;                              // class for tracking system energy drifts
  sys.init(Energy::systemEnergy(spc,pot,spc.p));// store initial total system energy

  cout << atom.info() + spc.info() + pot.info() + "\n";

  MCLoop loop(mcp);                             // class for handling mc loops
  while ( loop[0] ) {
    while ( loop[1] ) {
      if (slp_global() < 0.5)
        sys+=mv.move( salt.size() );            // translate salt
      else 
        sys+=gc.move();                         // grand canonical exchange
      widom1.sample(spc,pot,1);
      widom2.sample(spc.p,spc.geo);
    }                                           // end of micro loop
    sys.checkDrift(Energy::systemEnergy(spc,pot,spc.p)); // calc. energy drift
    cout << loop.timing();
  }                                             // end of macro loop

  FormatPQR::save("confout.pqr", spc.p);        // PQR snapshot for VMD etc.
  spc.save("state");                            // final simulation state

  UnitTest test(mcp);                           // class for unit testing
  gc.test(test);
  mv.test(test);
  sys.test(test);
  widom1.test(test);

  cout << loop.info() + sys.info() + mv.info() + gc.info() + test.info()
    + widom1.info() + widom2.info();

  return test.numFailed();
}

int main() {
  //::atom.includefile("stockmayer.json");         // load atom properties
  InputMap in("stockmayer.input");               // open parameter file for user input
  Energy::NonbondedVector<Tspace,Tpair> pot(in); // non-bonded only
  EnergyDrift sys;                               // class for tracking system energy drifts
  Tspace spc(in);                // create simulation space, particles etc.
  Group sol;
  sol.addParticles(spc, in);                     // group for particles
  MCLoop loop(in);                               // class for mc loop counting
  TmoveTran trans(in,pot,spc);
  TmoveRot rot(in,pot,spc);
  trans.setGroup(sol);                                // tells move class to act on sol group
  rot.setGroup(sol);                                  // tells move class to act on sol group
  spc.load("state");
  spc.trial = spc.p;
  UnitTest test(in);               // class for unit testing
  Analysis::DipoleAnalysis dian(spc,in);
  DipoleWRL sdp;
  FormatXTC xtc(spc.geo.len.norm());
  sys.init( Energy::systemEnergy(spc,pot,spc.p)  );   // initial energy
  while ( loop[0] ) {                         // Markov chain 
    while ( loop[1] ) {
      if (slp_global() > 0.5)
        sys+=trans.move( sol.size() );                // translate
      else
        sys+=rot.move( sol.size() );                  // rotate

      if (slp_global()<0.5)
        dian.sampleMuCorrelationAndKirkwood(spc);
      if (slp_global()>0.99)
        xtc.save(textio::prefix+"out.xtc", spc.p);  
      dian.sampleDP(spc);
    }    
    sys.checkDrift(Energy::systemEnergy(spc,pot,spc.p)); // compare energy sum with current
    cout << loop.timing() << std::flush;
  }

  // perform unit tests
  trans.test(test);
  rot.test(test);
  sys.test(test);
  sdp.saveDipoleWRL("stockmayer.wrl",spc,sol);
  FormatPQR().save("confout.pqr", spc.p);
  std::cout << spc.info() + pot.info() + trans.info()
    + rot.info() + sys.info() + test.info() + dian.info(); // final info
  dian.save();
  spc.save("state");

  return test.numFailed();
}

int main() {
  InputMap mcp("grand.json");                   // open user input file
  Tspace spc(mcp);                              // simulation space
  Energy::Nonbonded<Tspace,CoulombHS> pot(mcp); // hamiltonian
  pot.setSpace(spc);                            // share space w. hamiltonian

  spc.load("state",Tspace::RESIZE);             // load old config. from disk (if any)

  // Two different Widom analysis methods
  double lB = pot.pairpot.first.bjerrumLength();// get bjerrum length
  Analysis::Widom<PointParticle> widom1;        // widom analysis (I)
  Analysis::WidomScaled<Tspace> widom2(lB,1);   // ...and (II)
  widom1.add(spc.p);
  widom2.add(spc.p);

  Move::Propagator<Tspace> mv(mcp,pot,spc);

  EnergyDrift sys;                              // class for tracking system energy drifts
  sys.init(Energy::systemEnergy(spc,pot,spc.p));// store initial total system energy

  cout << atom.info() + spc.info() + pot.info() + "\n";

  MCLoop loop(mcp);                             // class for handling mc loops
  while ( loop[0] ) {
    while ( loop[1] ) {
      sys+=mv.move();                           // move!
      widom1.sample(spc,pot,1);
      widom2.sample(spc.p,spc.geo);
    }                                           // end of micro loop
    sys.checkDrift(Energy::systemEnergy(spc,pot,spc.p)); // calc. energy drift
    cout << loop.timing();
  }                                             // end of macro loop

  FormatPQR::save("confout.pqr", spc.p);        // PQR snapshot for VMD etc.
  spc.save("state");                            // final simulation state

  UnitTest test(mcp);                           // class for unit testing
  mv.test(test);
  sys.test(test);
  widom1.test(test);

  cout << loop.info() + sys.info() + mv.info() + test.info()
    + widom1.info() + widom2.info();

  return test.numFailed();
}

int main() {
  
  InputMap mcp("gcmol.input");     // open user input file
  MCLoop loop(mcp);                   // class for handling mc loops
  EnergyDrift sys;                    // class for tracking system energy drifts
  UnitTest test(mcp);

  // Create Space and a Hamiltonian with three terms
  Tspace spc(mcp);
  auto pot = Energy::Nonbonded<Tspace,Tpairpot>(mcp);

  //    ADD CONFIGURATIONS FOR POOL INSERTS
  string file = mcp.get<string>("polymer_file", "");
  p_vec conf;
  FormatAAM::load(file,conf);
  molecule.pushConfiguration("polymer", conf);   // p_vec is one molecule large
  molecule.pushConfiguration("polymer2",conf);

  // Markov moves and analysis
  Move::GCMolecular<Tspace> gc(mcp, pot, spc);

  sys.init( Energy::systemEnergy(spc,pot,spc.p)  );      // store initial total system energy

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

  while ( loop[0] ) {  // Markov chain
    while ( loop[1] ) {
      sys+=gc.move();
    }
    sys.checkDrift(Energy::systemEnergy(spc,pot,spc.p)); // compare energy sum with current
    cout << loop.timing();
  } // end of macro loop

  sys.test(test);
  gc.test(test);

  // print information
  cout << loop.info() << sys.info() << gc.info() << spc.info() <<test.info() << endl;

  // clean allocated memory
  spc.freeGroups();

  return test.numFailed();
}

int main() {

  cout << textio::splash();      // show faunus banner and credits
  InputMap mcp("membrane.json"); //read input file

  FormatXTC xtc(1000);
  EnergyDrift sys;               // class for tracking system energy drifts

  // Energy functions and space
  auto pot = Energy::NonbondedCutg2g<Tspace,Tpairpot>(mcp)
    + Energy::Bonded<Tspace>()
    + Energy::ExternalPressure<Tspace>(mcp)
    + Energy::EquilibriumEnergy<Tspace>(mcp);

  auto nonbonded = std::get<0>( pot.tuple() );
  auto bonded    = std::get<1>( pot.tuple() );

  nonbonded->noPairPotentialCutoff=true;
  Tspace spc(mcp);

  auto lipids = spc.findMolecules("lipid");

  Group allLipids( lipids.front()->front(), lipids.back()->back() );
  allLipids.setMolSize(3);
  MakeDesernoMembrane(allLipids, *bonded, nonbonded->pairpot, mcp);

  // Place all lipids in xy plane (z=0);
  for ( auto g : lipids ) {
    double dz = spc.p[ g->back() ].z();
    for (auto i : *g) {
      spc.p[i].z() -= dz;
      spc.geo.boundary( spc.p[i] );
    }
    g->setMassCenters( spc );
  }
  spc.trial=spc.p;   // sync. particle trial vector
  spc.load("state"); // load old config. from disk (if any)

  // Markov moves and analysis
  Move::Propagator<Tspace> mv(mcp, pot, spc);
  Analysis::BilayerStructure lipidstruct;
  Analysis::VirialPressure<Tspace> virial(mcp, pot, spc);

  sys.init( Energy::systemEnergy(spc,pot,spc.p)  ); // store total energy

  cout << atom.info() + spc.info() + pot.info();

  MCLoop loop(mcp);                      // class for handling mc loops
  while ( loop[0] ) {                    // Markov chain 
    while ( loop[1] ) {
      sys += mv.move();
      double ran = slump();
      if ( ran > 0.99 ) {
        xtc.setbox( spc.geo.len );
        xtc.save("traj.xtc", spc.p);
      }
      if ( ran > 0.90 ) {
        virial.sample();
        lipidstruct.sample(spc.geo, spc.p, allLipids);
      }

    } // end of micro loop

    sys.checkDrift(Energy::systemEnergy(spc,pot,spc.p)); // energy drift?

    // save to disk
    FormatPQR::save("confout.pqr", spc.p);
    spc.save("state");

    cout << loop.timing();
  } // end of macro loop

  // perform unit tests
  UnitTest test(mcp);
  mv.test(test);
  sys.test(test);
  lipidstruct.test(test);

  cout << loop.info() + mv.info()
    + lipidstruct.info() + sys.info() + virial.info() + test.info();

  return test.numFailed();
}

int main() {
  
  std::vector<Coord> sasaCoords;     //vector of sasa coordinates
  std::vector<Scalar> sasaWeights;   //vector of sasa weights
  
  clock_t clk_a, clk_b, clk_sum = 0;
  
  /*
  sasaCoords.push_back(Coord(1.0,2.0,3.0));
  sasaWeights.push_back(1.0);
  
  POWERSASA::PowerSasa<Scalar,Coord> *ps = 
    new POWERSASA::PowerSasa<Scalar,Coord>(sasaCoords, sasaWeights, 1, 1, 1, 1);
  ps->calc_sasa_all();
  
  
  Scalar volume = 0.0, surf = 0.0;
  for (unsigned int i = 0; i < sasaCoords.size(); ++i)
  {
	printf("%4d sasa=%7.3lf vol=%7.3lf\n", i, (ps->getSasa())[i], (ps->getVol())[i]);
	volume += (ps->getVol())[i];
	surf += (ps->getSasa())[i];
  }
  printf("volume=%lf\n", volume);
  printf("sasa  =%lf\n", surf);

  delete ps;
  */
  cout << textio::splash();           // show faunus banner and credits
  
  InputMap mcp("polymers.input");     // open user input file
  MCLoop loop(mcp);                   // class for handling mc loops
  EnergyDrift sys;                    // class for tracking system energy drifts
  UnitTest test(mcp);                 // class for unit testing

  // Create Space and a Hamiltonian with three terms
  Tspace spc(mcp);
  auto pot = Energy::Nonbonded<Tspace,Tpairpot>(mcp)
    + Energy::SASAEnergy<Tspace>(mcp)
    + Energy::ExternalPressure<Tspace>(mcp) + Energy::Bonded<Tspace>();

  auto bonded = &pot.second; // pointer to bond energy class

  // Add salt
  Group salt;
  salt.addParticles(spc, mcp);

  // Add polymers
  vector<Group> pol( mcp.get("polymer_N",0));            // vector of polymers
  string file = mcp.get<string>("polymer_file", "");
  double req = mcp.get<double>("polymer_eqdist", 0);
  double k   = mcp.get<double>("polymer_forceconst", 0);
  for (auto &g : pol) {                                  // load polymers
    Tspace::ParticleVector v;                            // temporary, empty particle vector
    FormatAAM::load(file,v);                             // load AAM structure into v
    Geometry::FindSpace().find(spc.geo,spc.p,v);         // find empty spot in particle vector
    g = spc.insert(v);                                   // insert into space
    g.name="Polymer";
    spc.enroll(g);
    for (int i=g.front(); i<g.back(); i++)
      bonded->add(i, i+1, Potential::Harmonic(k,req));   // add bonds
  }
  Group allpol( pol.front().front(), pol.back().back() );// make group w. all polymers

  // Markov moves and analysis
  Move::Isobaric<Tspace> iso(mcp,pot,spc);
  Move::TranslateRotate<Tspace> gmv(mcp,pot,spc);
  Move::AtomicTranslation<Tspace> mv(mcp, pot, spc);
  Move::CrankShaft<Tspace> crank(mcp, pot, spc);
  Move::Pivot<Tspace> pivot(mcp, pot, spc);
  Analysis::PolymerShape shape;
  Analysis::RadialDistribution<> rdf(0.2);
  Scatter::DebyeFormula<Scatter::FormFactorUnity<>> debye(mcp);

  spc.load("state");                                     // load old config. from disk (if any)
  sys.init( Energy::systemEnergy(spc,pot,spc.p)  );      // store initial total system energy

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

  while ( loop.macroCnt() ) {  // Markov chain 
    while ( loop.microCnt() ) {
      int k,i=slp_global.rand() % 6;
      switch (i) {
        case 0:
          mv.setGroup(salt);
          sys+=mv.move( salt.size() );  // translate salt
          break;
        case 1:
          mv.setGroup(allpol);
          sys+=mv.move( allpol.size() );// translate monomers
          for (auto &g : pol) {
            g.setMassCenter(spc);
            shape.sample(g,spc);        // sample gyration radii etc.
          }
          break;
        case 2:
          k=pol.size();
          while (k-->0) {
            gmv.setGroup( pol[ slp_global.rand() % pol.size() ] );
            sys+=gmv.move();            // translate/rotate polymers
//.........这里部分代码省略.........

int main() {
  cout << textio::splash();          // Spam

  InputMap mcp("slitpolymer.input"); // Open input parameter file
  MCLoop loop(mcp);                  // handle mc loops
  EnergyDrift sys;                   // track system energy drifts
  UnitTest test(mcp);                // unit testing
  Geometry::Cuboidslit geo(mcp);     // rectangular simulation box w. XY periodicity

  Energy::ExternalPotential< Potential::GouyChapman<> > pot(mcp);
  pot.setGeometry(geo);              // Pass on geometry to potential
  pot.expot.setSurfPositionZ( &geo.len_half.z() ); // Pass position of GC surface
  Space spc( pot.getGeometry() );    // Simulation space (all particles and group info)

  // Load and add polymer to Space
  FormatAAM aam;                                      // AAM structure file I/O
  string polyfile = mcp.get<string>("polymer_file", "");
  aam.load(polyfile);                                 // Load polymer structure into aam class
  Geometry::FindSpace().find(*spc.geo, spc.p, aam.particles()); // find empty spot
  GroupMolecular pol = spc.insert( aam.particles() ); // Insert into Space and return matching group
  pol.name="polymer";                                 // Give polymer arbitrary name
  spc.enroll(pol);                                    // Enroll polymer in Space

  // MC moves
  Move::TranslateRotate gmv(mcp,pot,spc);
  Move::CrankShaft crank(mcp,pot,spc);
  Move::Pivot pivot(mcp,pot,spc);
  Move::Reptation rep(mcp,pot,spc);

  Analysis::PolymerShape shape;                       // sample polymer shape
  Analysis::LineDistribution<> surfmapall;            // monomer-surface histogram
  spc.load("state");                                  // Load start configuration, if any
  sys.init( Energy::systemEnergy(spc,pot,spc.p) );    // Store total system energy

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

  while ( loop.macroCnt() ) {  // Markov chain 
    while ( loop.microCnt() ) {
      int i=slp_global.rand() % 4;
      switch (i) {
        case 0: // translate and rotate polymer
          gmv.setGroup(pol);
          sys+=gmv.move(); 
          break;
        case 1: // pivot
          pivot.setGroup(pol);
          sys+=pivot.move();
          break;
        case 2: // crankshaft
          crank.setGroup(pol);
          sys+=crank.move();
          break;
        case 3: // reptation
          rep.setGroup(pol);
          sys+=rep.move();
          break;
      }

      double rnd = slp_global(); // [0:1[
      if (rnd<0.05)
        shape.sample(pol,spc);   // sample polymer shape - gyration radius etc.
      if (rnd<0.05)
        for (auto i : pol)
          surfmapall( pot.expot.surfDist( spc.p[i] ) )++;  // sample monomer distribution

    } // end of micro loop

    sys.checkDrift( Energy::systemEnergy(spc,pot,spc.p) ); // re-calc system energy and detect drift
    cout << loop.timing();                                 // print timing and ETA information

  } // end of macro loop

  spc.save("state");               // save final state of simulation (positions etc)
  surfmapall.save("surfall.dat");  // save monomer-surface distribution
  FormatPQR().save("confout.pqr", spc.p);  // save PQR file

  // Perform unit tests (only for faunus integrity)
  gmv.test(test);
  sys.test(test);
  pivot.test(test);
  crank.test(test);
  rep.test(test);
  shape.test(test);

  cout << loop.info() + sys.info() + gmv.info() + crank.info()
    + pivot.info() + rep.info() + shape.info() + spc.info() + test.info();

  return test.numFailed();
}