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

int main() {
  cout << textio::splash();
  InputMap mcp("cigars2fibrils.json");
  FormatMXYZ mxyz;
  EnergyDrift sys;

  // Energy functions and space
  Tspace spc(mcp);
  Energy::NonbondedVector<Tspace,Tpairpot> pot(mcp);

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

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

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

  MCLoop loop( mcp );
  mxyz.save( "cigars2fibrils-mov", spc.p, spc.geo.len, loop.innerCount() );

  while ( loop[0] ) {  // Markov chain 
    while ( loop[1] )
      sys += mv.move();

    sys.checkDrift(Energy::systemEnergy(spc,pot,spc.p)); // compare energy sum with current

    cout << loop.timing();

    mxyz.save( "cigars2fibrils-mov", spc.p, spc.geo.len, loop.innerCount() );
  } // end of macro loop

  // print information about simulationat the end
  cout << loop.info() + sys.info() + mv.info();
}

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() {
  cout << textio::splash();                 // show faunus banner and credits
  //InputMap mcp("cigars2fibrils.input");     // open user input file
  InputMap mcp("psc-nvt.input");
  MCLoop loop(mcp);                         // class for handling mc loops
  FormatMXYZ mxyz;                          // MXYZ structure file I/O
  EnergyDrift sys;                         // class for tracking system energy drifts

    
  // Energy functions and space
  Tspace spc(mcp);
  auto pot = Energy::NonbondedVector<Tspace,Tpairpot>(mcp);

  // Markov moves and analysis
  Move::AtomicTranslation<Tspace> mv(mcp, pot, spc);
  Move::AtomicRotation<Tspace> rot(mcp, pot, spc);

  cout <<"before adding cigars";

  // Add cigars
  Group cigars;
  cigars.addParticles(spc, mcp);
  cigars.name="PSC";
  for (auto i : cigars) {
    spc.p[i].dir.ranunit(slump);
    spc.p[i].patchdir.ranunit(slump);
    Geometry::cigar_initialize(spc.geo, spc.p[i]);
    spc.trial[i]=spc.p[i];
  }
  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!");

  mxyz.save("cigars2fibrils-mov", spc.p, spc.geo.len, loop.count());
  while ( loop.macroCnt() ) {  // Markov chain 
    while ( loop.microCnt() ) {
      int i= slump.rand() % 2;
      switch (i) {
        case 0:
          mv.setGroup(cigars);
          sys+=mv.move( cigars.size() );  // translate cigars
          break;
        case 1:
          rot.setGroup(cigars);
          sys+=rot.move( cigars.size() ); // translate cigars
          break;
      }
    } // end of micro loop
  sys.checkDrift(Energy::systemEnergy(spc,pot,spc.p)); // compare energy sum with current

  cout << loop.timing();

  mxyz.save("cigars2fibrils-mov", spc.p, spc.geo.len, loop.count());
  } // end of macro loop

  // print information
  cout << loop.info() << sys.info() << mv.info() << rot.info();
}

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<Tpair,Tgeo> pot(in);   // create Hamiltonian, non-bonded only
  EnergyDrift sys;                               // class for tracking system energy drifts
  Space spc( pot.getGeometry() );                // create simulation space, particles etc.
  GroupAtomic sol(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.05);       // particle-particle g(r)

  Move::AtomicTranslation trans(in, pot, spc);   // particle move class
  Move::AtomicRotation rot(in, pot, spc);        // particle move class
  //PolarizeMove<AtomicTranslation> trans(in,pot,spc);
  //PolarizeMove<AtomicRotation> 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_ST");


  double limit = 1e-6;
  Point ExternalField(0,0,0);

  sys.init( Energy::systemEnergy(spc,pot,spc.p)  );      // store initial total system 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()<0.5)
        for (auto i=sol.front(); i<sol.back(); i++) { // salt radial distribution function
          for (auto j=i+1; j<=sol.back(); j++) {
            double r =pot.geometry.dist(spc.p[i],spc.p[j]); 
            rdf(r)++;
            mucorr(r) += spc.p[i].mu.dot(spc.p[j].mu)/(spc.p[i].muscalar*spc.p[j].muscalar);
          }
        }
    }

    cout << "Eps: " << getDielectricConstant(pot,spc,in.get<double>("dipdip_cutoff",pc::infty)) << "\n";
    sys.checkDrift(Energy::systemEnergy(spc,pot,spc.p)); // compare energy sum with current
    cout << loop.timing();
  }

  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(); // final info
  spc.save("state_ST");
}

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() {
  cout << textio::splash();            // faunus splash info!
  ::atom.includefile("titrate.json");  // load atom properties
  InputMap mcp("titrate.input");
  EnergyDrift sys;                     // track system energy drift
  UnitTest test(mcp);

  // Space and hamiltonian
  Tspace spc(mcp);
  auto pot = Energy::Nonbonded<Tspace,Potential::DebyeHuckel>(mcp)
    + Energy::EquilibriumEnergy<Tspace>(mcp);

  // Add molecule to middle of simulation container
  string file = mcp.get<string>("molecule", string());
  Tspace::ParticleVector v;
  FormatAAM::load(file,v);
  Geometry::cm2origo(spc.geo,v); // center molecule
  Group g = spc.insert(v);       // insert into space
  g.name="peptide"; // babtise
  spc.enroll(g);    // enroll group in space

  spc.load("state");

  Move::SwapMove<Tspace> tit(mcp,pot,spc,pot.second);
  Analysis::ChargeMultipole mp;

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

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

  MCLoop loop(mcp);            // class for handling mc loops
  while ( loop[0] ) {  // Markov chain 
    while ( loop[1] ) {
      sys+=tit.move();
      mp.sample(g,spc);
    } // end of micro loop

    sys.checkDrift( Energy::systemEnergy(spc,pot,spc.p) );
    cout << loop.timing();

  } // end of macro loop

  cout << loop.info() + sys.info() + tit.info() + g.info() + mp.info();

  FormatPQR().save("confout.pqr", spc.p);
  spc.save("state");
}

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

//.........这里部分代码省略.........
            FormatPQR::load(file,v);
        else
            FormatAAM::load(file,v);

        Geometry::FindSpace fs;
        if (inPlane)
            fs.dir=Point(0,0,1);
        fs.find(spc.geo,spc.p,v);
        pol[i] = spc.insert(v);
        pol[i].name=file+std::to_string(i);
        spc.enroll( pol[i] );
    }
    Group allpol( pol.front().front(), pol.back().back() );

    // Add atomic species
    Group salt;
    salt.addParticles(spc, mcp);
    salt.name="Atomic Species";
    spc.enroll(salt);

    Analysis::LineDistribution<> rdf(0.5);
    Analysis::ChargeMultipole mpol;
    Analysis::MultipoleDistribution<Tspace> mpoldist(mcp);
    FormatQtraj qtrj("qtrj.dat");

    spc.load("state");

    Move::Isobaric<Tspace> iso(mcp,pot,spc);
    Move::TranslateRotate<Tspace> gmv(mcp,pot,spc);
    Move::AtomicTranslation<Tspace> mv(mcp, pot, spc);
    Move::SwapMove<Tspace> tit(mcp,pot,spc,*eqenergy);
    if (inPlane)
        for (auto &m : pol)
            gmv.directions[ m.name ]=Point(0,0,1);

    eqenergy->eq.findSites(spc.p);

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

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

    while ( loop.macroCnt() ) {  // Markov chain
        while ( loop.microCnt() ) {
            int k,i=rand() % 2;
            switch (i) {
            case 0:
                k=pol.size();
                while (k-->0) {
                    gmv.setGroup( pol[ rand() % pol.size() ] );
                    sys+=gmv.move();
                }
                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) )++;
                break;
            case 1:
                sys+=tit.move();
                break;
            case 2:
                sys+=iso.move();
                break;
            case 3:
                mv.setGroup(salt);
                sys+=mv.move();
                break;
            }

            double xi = slp_global(); // random number [0,1)

            // Sample multipolar moments and distribution
            if ( xi>0.95 ) {
                pol[0].setMassCenter(spc);
                pol[1].setMassCenter(spc);
                mpol.sample(pol,spc);
                mpoldist.sample(spc, pol[0], pol[1]);
            }

            if ( movie==true && xi>0.995 ) {
                xtc.setbox( 1000. );
                xtc.save("traj.xtc", spc.p);
                qtrj.save(spc.p);
            }
        } // end of micro loop

        sys.checkDrift( Energy::systemEnergy(spc,pot,spc.p) ); // detect energy drift
        cout << loop.timing();

        spc.save("state");
        mcp.save("mdout.mdp");
        rdf.save("rdf_p2p.dat");
        mpoldist.save("multipole.dat");
        FormatPQR::save("confout.pqr", spc.p);

    } // end of macro loop

    cout << loop.info() + sys.info() + gmv.info() + mv.info()
         + iso.info() + tit.info() + mpol.info();

}

//.........这里部分代码省略.........
  auto nonbonded = pot.create(nb);
#endif
#ifdef org
  Energy::Nonbonded<Potential::PotentialMap<Tpairpot>,Tgeometry> nb(mcp); //  Non-tabulation
  nb.pairpot.add(atom["Ca"].id,atom["Ca"].id,Tpairpot(mcp));
  nb.pairpot.add(atom["Ca"].id,atom["Cl"].id,Tpairpot(mcp));
  nb.pairpot.add(atom["Ca"].id,atom["Na"].id,Tpairpot(mcp));
  nb.pairpot.add(atom["Na"].id,atom["Na"].id,Tpairpot(mcp));
  nb.pairpot.add(atom["Na"].id,atom["Cl"].id,Tpairpot(mcp));
  nb.pairpot.add(atom["Cl"].id,atom["Cl"].id,Tpairpot(mcp));
  auto nonbonded = pot.create(nb);
#endif
#ifdef orgopt
  Energy::Nonbonded<Potential::PotentialVec<Tpairpot>,Tgeometry> nb(mcp); //  Non-tabulation
  nb.pairpot.add(atom["Ca"].id,atom["Ca"].id,Tpairpot(mcp));
  nb.pairpot.add(atom["Ca"].id,atom["Cl"].id,Tpairpot(mcp));
  nb.pairpot.add(atom["Ca"].id,atom["Na"].id,Tpairpot(mcp));
  nb.pairpot.add(atom["Na"].id,atom["Na"].id,Tpairpot(mcp));
  nb.pairpot.add(atom["Na"].id,atom["Cl"].id,Tpairpot(mcp));
  nb.pairpot.add(atom["Cl"].id,atom["Cl"].id,Tpairpot(mcp));
  auto nonbonded = pot.create(nb);
#endif
  

#ifdef tabsingle
  auto nonbonded = pot.create( Energy::Nonbonded<Potential::PotentialTabulate<Tpairpot>,Tgeometry>(mcp) );
#endif 
#ifdef taboptsingle
  auto nonbonded = pot.create( Energy::Nonbonded<Potential::PotentialTabulateVec<Tpairpot>,Tgeometry>(mcp) );
#endif 
#ifdef orgsingle
  auto nonbonded = pot.create( Energy::Nonbonded<Tpairpot,Tgeometry>(mcp) );
#endif 
  
  //nb.pairpot.print_tabulation(); // Print files to compare tabulation with real potential
  
  
  Space spc( pot.getGeometry() );

  // Add salt
  GroupAtomic salt(spc, mcp);
  salt.name="Salt";
  spc.enroll(salt);

  spc.load("state");

  Move::AtomicTranslation mv(mcp, pot, spc);
  mv.setGroup(salt);   // specify atomic particles to be moved

  short na = atom["Na"].id;
  short ca = atom["Ca"].id;
  short cl = atom["Cl"].id;
  
  Analysis::RadialDistribution<float,unsigned int> rdf1(0.2); // 0.2 Å resolution
  Analysis::RadialDistribution<float,unsigned int> rdf2(0.2); // 0.2 Å resolution

  
  
  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() ) {
      xtc.save("out.xtc", spc.p);
      int i=rand() % 1;
      switch (i) {
        case 0:
          mv.setGroup(salt);
          sys+=mv.move( salt.size()/2+1 );
          break;
      }
      if (slp_global() > 0.9) {
        rdf1.sample( spc, salt, na, cl );
        rdf2.sample( spc, salt, ca, cl );
      }
    } // end of micro loop

    sys.checkDrift( Energy::systemEnergy(spc,pot,spc.p) );
    spc.save("state");
    cout << loop.timing();
  } // end of macro loop
  
#ifdef tab
  rdf1.save("naclrdf_tab.dat");
  rdf2.save("caclrdf_tab.dat");
#endif
#ifdef tabopt
  rdf1.save("naclrdf_tabopt.dat");
  rdf2.save("caclrdf_tabopt.dat");
#endif
#ifdef org
  rdf1.save("naclrdf_notab.dat");
  rdf2.save("caclrdf_notab.dat");
#endif

  pqr.save("confout.pqr", spc.p);

  cout << loop.info() << spc.info() << sys.info() << mv.info();
}

int main(int argc, char** argv) {
  InputMap mcp("cyl.input");
  MCLoop loop(mcp);                    // class for handling mc loops
  FormatXTC xtc(1000);                 // XTC gromacs trajectory format
  EnergyDrift sys;                     // class for tracking system energy drifts
  UnitTest test(mcp);
  bool inPlane = mcp.get<bool>("molecule_plane");

  Tspace spc(mcp);
  auto pot = Energy::Nonbonded<Tspace,Tpairpot>(mcp)
    + Energy::MassCenterConstrain<Tspace>(mcp) 
    + Energy::EquilibriumEnergy<Tspace>(mcp);

  auto constrainenergy = &pot.first.second;
  auto eqenergy = &pot.second;

  // Add molecules
  int N1 = mcp.get("molecule1_N",0);
  int N2 = mcp.get("molecule2_N",0);
  string file;
  vector<Group> pol(N1+N2);
  for (int i=0; i<N1+N2; i++) {
    if (i>=N1)
      file = mcp.get<string>("molecule2");
    else
      file = mcp.get<string>("molecule1");
    Tspace::ParticleVector v;
    FormatAAM::load(file,v);
    Geometry::FindSpace fs;
    if (inPlane)
      fs.dir=Point(0,0,1);
    fs.find(spc.geo,spc.p,v);
    pol[i] = spc.insert(v);
    pol[i].name=file+std::to_string(i);
    spc.enroll( pol[i] );
  }
  Group allpol( pol.front().front(), pol.back().back() );

  constrainenergy->addPair(pol[0], pol[1]);

  // Add atomic species
  Group salt;
  salt.addParticles(spc, mcp);
  salt.name="Atomic Species";
  spc.enroll(salt);

  Analysis::LineDistribution<> rdf(0.5);
  Analysis::ChargeMultipole mpol;

  spc.load("state");

  Move::Isobaric<Tspace> iso(mcp,pot,spc);
  Move::TranslateRotate<Tspace> gmv(mcp,pot,spc);
  Move::AtomicTranslation<Tspace> mv(mcp, pot, spc);
  Move::SwapMove<Tspace> tit(mcp,pot,spc,*eqenergy);
  if (inPlane)
    for (auto &m : pol)
      gmv.directions[ m.name ]=Point(0,0,1);

  eqenergy->eq.findSites(spc.p);

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

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

  while ( loop.macroCnt() ) {  // Markov chain 
    while ( loop.microCnt() ) {
      int k,i=rand() % 2;
      switch (i) {
        case 0:
          k=pol.size();
          while (k-->0) {
            gmv.setGroup( pol[ rand() % pol.size() ] );
            sys+=gmv.move();
          }
          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) )++;
          break;
        case 1:
          sys+=tit.move();
          if( slp_global()>0.9 )
            mpol.sample(pol,spc);
          break;
        case 2:
          sys+=iso.move();
          break;
        case 3:
          mv.setGroup(salt);
          sys+=mv.move();
          break;
      }
      if ( slp_global()<-0.001 ) {
        xtc.setbox( 1000. );
        xtc.save("traj.xtc", spc);
      }
    } // end of micro loop

    sys.checkDrift( Energy::systemEnergy(spc,pot,spc.p) ); // detect energy drift
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            sys+=crank.move();          // crank shaft move
          }
          break;
        case 5:
          k=pol.size();
          while (k-->0) {
            pivot.setGroup( pol[ slp_global.rand() % pol.size() ] );
            sys+=pivot.move();          // pivot move
          }
          break;
      }

      // polymer-polymer mass center rdf
      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) )++;


      // SASA
      sasaCoords.clear();
      sasaWeights.clear();

/*
      for (auto i = pol.begin(); i != pol.end(); i++) {
		sasaCoords.push_back(i->cm);
		sasaWeights.push_back(5.0);
	  }
*/

      /*
      for (auto i = allpol.begin(); i != allpol.end(); i++) {
		auto monomer = spc.p.at(*i);
		sasaCoords.push_back(monomer);
		sasaWeights.push_back(10.0);
	  }
      clk_a = clock();
      auto ps = new POWERSASA::PowerSasa<Scalar,Coord>(sasaCoords, sasaWeights, 1, 1, 1, 1);
      ps->calc_sasa_all();
      clk_b = clock();
      clk_sum += clk_b - clk_a;
      
      {
		  Scalar volume = 0.0, surf = 0.0;
		  for (unsigned int i = 0; i < sasaCoords.size(); ++i)
		  {
			volume += (ps->getVol())[i];
			surf += (ps->getSasa())[i];
		  }
		  //printf("volume=%lf\n", volume);
		  //printf("sasa  =%lf\n", surf);
	  }
      delete ps;
*/
    } // end of micro loop

    // sample scattering
    debye.sample(spc.p);

    sys.checkDrift(Energy::systemEnergy(spc,pot,spc.p)); // compare energy sum with current
    cout << loop.timing();

  } // end of macro loop

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

  // perform unit tests
  iso.test(test);
  gmv.test(test);
  mv.test(test);
  sys.test(test);

  // print information
  cout << loop.info() << sys.info() << mv.info() << gmv.info() << crank.info()
    << pivot.info() << iso.info() << shape.info() << test.info();

  cout << endl << endl << "SASA" << endl;
  cout << "Time " << (float) clk_sum/CLOCKS_PER_SEC << " s = " << (float) clk_sum/CLOCKS_PER_SEC/3600 << " h" << endl;
  
  cout << sasaCoords.size() << endl;

  cout << pot.first.first.second.info() ;
  
  /*
  cout << allpol.size() << endl;
  cout << allpol.info() << endl;
  cout << spc.info() << endl;
  cout << typeid(spc.p.at(34)).name() << endl;
  cout << typeid(allpol.begin()).name() << endl;
  cout << spc.findGroup(34)->info() << endl;
*/	
  
  //for (auto i = allpol.begin(); i != allpol.end(); i++) {
  //}

  return test.numFailed();
}