search

hhatch / feasst

FEASST: Free Energy and Advanced Sampling Simulation Toolkit (prototype of https://pages.nist.gov/feasst)
Other
2 stars 1 forks source link

.num_values() for accumulator of moments does not work #5

Open mahynski opened 3 years ago

mahynski commented 3 years ago

This was compiled using the hwh/acc branch:

commit ba51c74c66cf162e6d49cf7c5cf4ac20e60a321f (HEAD -> hwh/acc, origin/hwh/acc) Author: Harold Hatch harold.hatch@nist.gov Date: Tue Dec 1 15:38:29 2020 -0500

The run.cc file

#include <assert.h>
#include <fstream>
#include "feasst.h"

feasst::ArgumentParse args("A grand canonical ensemble flat histogram Monte Carlo simulation of a bulk fluid", {
  {"--task", "SLURM job array index", "0"},
  {"--num_procs", "number of processors", "1"},
  {"--num_hours", "number of hours before restart", feasst::str(3*24)},
  {"--dccb_cutoff", "cutoff for dual cut configurational bias", "4"},
  {"--dccb_begin", "number of molecules before using DCCB", "300"},
  {"--lx", "box length in x", "10.0"},
  {"--ly", "box length in y", "10.0"},
  {"--lz", "box length in z", "10.0"},
  {"--min_particles", "minimum number of particles", "0"},
  {"--max_particles", "maximum number of particles", "1"},
  {"--temperature", "temperature", "263.15"},
  {"--particle0", "particle type 0"},
  {"--particle1", "particle type 1, not used if empty (the default)"},
  {"--collect_flatness", "number of WL flatness to begin collection", "18"},
  {"--min_flatness", "number of WL flatness to switch to TM", "22"},
  {"--beta_mu", "baseline chemical potential of each species", "0."},
  {"--delta_beta_mu1", "beta_mu1 = beta_mu0 + delta_beta_mu1", "0."},
  {"--cyl_radius", "radius of cylindrical confinement, not used if < 0", "-1"},
  {"--cyl_cutoff", "square well interaction distance from cylinder to point", "10"},
  {"--cyl_epsilon", "square well interaction strength", "500"},
  {"--file_xyz", "file name of xyz positions of initial configuration.\n"
    "Assume all particle0, and do not use if empty (the default)"}});

std::shared_ptr<feasst::MonteCarlo> mc(const int thread, const int mn, const int mx) {
  const std::string steps_per = feasst::str(1e5);
  auto mc = feasst::MakeMonteCarlo();
  feasst::argtype domain_args = {
    {"side_length0", args.get("--lx")},
    {"side_length1", args.get("--ly")},
    {"side_length2", args.get("--lz")}};
  const int dccb_begin = args.get_int("--dccb_begin");
  std::string ref("-1");
  std::string num("1");
  const int dccb_cutoff = args.get_int("--dccb_cutoff");
  if (mx > dccb_begin) {
    domain_args.insert({"init_cells", feasst::str(dccb_cutoff)});
    ref = "0";
    num = "4";
  }
  const double beta = 1./args.get_double("--temperature");
  const double beta_mu = args.get_double("--beta_mu");
  feasst::argtype config_args = {{"particle_type0", args.get("--particle0")}};
  feasst::argtype criteria_args = {{"beta", feasst::str(beta)},
    {"chemical_potential0", feasst::str(beta_mu/beta)}};
  if (args.option_given("--particle1")) {
    config_args.insert({"particle_type1", args.get("--particle1")});
    criteria_args.insert({"chemical_potential1", feasst::str((beta_mu + args.get_double("--delta_beta_mu1"))/beta)});
  }
  if (args.option_given("--file_xyz")) {
    assert(thread == 0); // cannot initialize multiple threads to same configuration
    feasst::Configuration config(feasst::MakeDomain(domain_args), config_args);
    feasst::FileXYZ().load(args.get("--file_xyz"), &config);
    mc->add(config);
  } else {
    mc->add(feasst::Configuration(feasst::MakeDomain(domain_args), config_args));
  }
  mc->add(feasst::Potential(feasst::MakeLennardJones()));
  mc->add(feasst::Potential(feasst::MakeLongRangeCorrections()));
  const double cyl_radius = args.get_double("--cyl_radius");
  if (cyl_radius > 0) {
    feasst::Potential cylinder(feasst::MakeModelSquareWellShape(feasst::MakeCylinder(
      {{"radius", feasst::str(cyl_radius)}},
      feasst::Position({{"x", "0"}, {"y", "0"}, {"z", "0"}}),
      feasst::Position({{"x", "0"}, {"y", "0"}, {"z", "1"}}))));
    cylinder.set_model_params(mc->configuration());
    for (int site_type = 0; site_type < mc->configuration().num_site_types(); ++site_type) {
      cylinder.set_model_param("cutoff", site_type, args.get_double("--cyl_cutoff"));
      cylinder.set_model_param("epsilon", site_type, args.get_double("--cyl_epsilon"));
    }
    mc->add(cylinder);
  }
  if (mx > dccb_begin) {
    feasst::Potential reference(feasst::MakeLennardJones());
    if (mc->configuration().domain().num_cells() > 0) {
      reference = feasst::Potential(feasst::MakeLennardJones(), feasst::MakeVisitModelCell());
    }
    reference.set_model_params(mc->configuration());
    for (int site_type = 0; site_type < mc->configuration().num_site_types(); ++site_type) {
      reference.set_model_param("cutoff", site_type, dccb_cutoff);
    }
    mc->add_to_reference(reference);
  }
  mc->set(feasst::MakeThermoParams(criteria_args));
  mc->set(feasst::MakeFlatHistogram(
    feasst::MakeMacrostateNumParticles(
      feasst::Histogram({{"width", "1"}, {"max", feasst::str(mx)}, {"min", feasst::str(mn)}})),
    feasst::MakeWLTM({
      {"collect_flatness", args.get("--collect_flatness")},
      {"min_flatness", args.get("--min_flatness")},
      {"min_sweeps", "1000000"}}),
    criteria_args));
  std::cout << "Initial energy: " << MAX_PRECISION << mc->criteria().current_energy() << std::endl;
  for (int particle_type = 0; particle_type < mc->configuration().num_particle_types(); ++particle_type) {
    mc->add(feasst::MakeTrialTranslate({
      {"particle_type", feasst::str(particle_type)},
      {"weight", "10."},
      {"tunable_param", "1."},
      {"reference_index", ref},
      {"num_steps", num}}));
    mc->add(feasst::MakeTrialRotate({
      {"particle_type", feasst::str(particle_type)},
      {"weight", "4."},
      {"tunable_param", "1."},
      {"reference_index", ref},
      {"num_steps", num}}));
    mc->add(feasst::MakeTrialTransfer({
      {"particle_type", feasst::str(particle_type)},
      {"weight", "1"},
      {"reference_index", ref},
      {"num_steps", num}}));
  }
  mc->add(feasst::MakeCheckEnergy({{"steps_per", steps_per}, {"tolerance", "0.0001"}}));
  mc->add(feasst::MakeTuner({{"steps_per", steps_per}, {"stop_after_phase", "0"}}));
  //mc->add(feasst::MakeLogAndMovie({{"steps_per", steps_per},
  //                                 {"file_name", "clones" + feasst::str(thread)},
  //                                 {"file_name_append_phase", "True"}}));
  mc->add(feasst::MakeEnergy({
    {"file_name", "en" + feasst::str(thread) + ".txt"},
    {"file_name_append_phase", "true"},
    {"start_after_phase", "0"},
    {"steps_per_write", steps_per},
    {"steps_per_update", "100"},
    {"multistate", "True"}}));
  if (mc->configuration().num_particle_types() > 1) {
    mc->add(feasst::MakeNumParticles({
      {"particle_type", "0"},
      {"file_name", "num" + feasst::str(thread) + ".txt"},
      {"file_name_append_phase", "true"},
      {"start_after_phase", "0"},
      {"steps_per_write", steps_per},
      {"steps_per_update", "100"},
      {"multistate", "True"}}));
  }
  mc->add(feasst::MakeExtensiveMoments({
    {"steps_per_write", steps_per},
    {"steps_per_update", "100"},
    {"file_name", "extmom"+feasst::str(thread) + ".txt"},
    {"file_name_append_phase", "true"},
    {"start_after_phase", "0"}, // do not update until equilibration complete (phase 1)
    {"max_order", "3"},
    {"multistate", "True"}}));
  mc->add(feasst::MakeCriteriaUpdater({{"steps_per", steps_per}}));
  mc->add(feasst::MakeCriteriaWriter({{"steps_per", steps_per},
                                     {"file_name", "clones" + feasst::str(thread) + "_crit.txt"},
                                     {"file_name_append_phase", "True"}}));
  mc->set(feasst::MakeCheckpoint({{"file_name", "checkpoint" + feasst::str(thread) + ".fst"},
                                 {"num_hours_terminate", feasst::str(0.99*args.get_int("--num_procs")*args.get_double("--num_hours"))}}));
  mc->add(feasst::MakeCheckRigidBonds({{"steps_per", steps_per}}));
  return mc;
}

int main(int argc, char ** argv) {
    std::cout << feasst::ThreadOMP().is_enabled() << std::endl;

  std::cout << feasst::version() << std::endl;
  std::cout << "args: " << args.parse(argc, argv) << std::endl;
  const int max_particles = args.get_int("--max_particles");
  const int min_particles = args.get_int("--min_particles");
  const int num_procs = args.get_int("--num_procs");
  auto windows = feasst::WindowExponential({
    {"alpha", "1.75"},
    {"num", feasst::str(num_procs)},
    {"minimum", feasst::str(min_particles)},
    {"maximum", feasst::str(max_particles)},
    {"extra_overlap", "2"}}).boundaries();
  std::cout << feasst::feasst_str(windows) << std::endl;
  auto clones = feasst::MakeClones();
  if (args.get_int("--task") == 0) {
    for (int proc = 0; proc < static_cast<int>(windows.size()); ++proc) {
      clones->add(mc(proc, windows[proc][0], windows[proc][1]));
    }
    clones->set(feasst::MakeCheckpoint({{"file_name", "checkpoint.fst"}}));
  } else {
    clones = feasst::MakeClones("checkpoint", num_procs);
  }
  clones->initialize_and_run_until_complete({{"ln_prob_file", "ln_prob.txt"}});
  std::cout << feasst::feasst_str(clones->ln_prob().values()) << std::endl;
  std::ofstream file("clones.fst");
  file << clones->serialize();
  file.close();
}

The CMake file

cmake_minimum_required(VERSION 3.5)
set(CMAKE_CXX_STANDARD 11)

# Tell cmake where to find the feasst library
set(CMAKE_PREFIX_PATH "$ENV{HOME}/feasst/build/")

find_package(feasst REQUIRED)

set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -pedantic -g -pthread")

# set the C++11
if (CMAKE_VERSION VERSION_LESS "3.1")
  if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
    set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++0x")
  endif ()
  if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
    add_definitions(-std=c++0x)
  endif ()
  if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
    add_definitions(-std=c++11)
  endif ()
else ()
  set (CMAKE_CXX_STANDARD 11)
endif ()

# OMP
find_package(OpenMP)
if(OPENMP_FOUND)
  set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
  set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
endif()

include_directories("${CMAKE_PREFIX_PATH}/include")
add_executable (run run.cc)
target_link_libraries (run LINK_PUBLIC feasst)

The execution

$ run --dccb_cutoff 5.0 --dccb_begin 400 --lx 42.0 --ly 42.0 --lz 42.0 --temperature 236.0 --particle0 ../data.ethane --max_particles 800 --beta_mu 0.0 --cyl_radius -1 --num_procs 18 --task 0 --num_hours 24

The particle data (data.ethane)

# LAMMPS-inspired data file

2 sites
1 bonds

1 site types
1 bond types

Units

length Angstrom
energy K

Site Labels

0 C

Site Properties

0 sigma 3.75 epsilon 98 cutoff 14

Sites

0 0  0                    0                    0
1 0  1.54                 0                    0

Bond Properties

0 length 1.54 delta 0.00001

Bonds

0 0 0 1

The problem

In python, the .average() for a moment returns "inf" because the .num_values() ~ 1e-310. For example

>>> import feasst as fst
>>> import pyfeasst
>>> clones = fst.MakeClones('T=236/checkpoint', 18)
>>> extmom_index = fst.SeekAnalyze().index("ExtensiveMoments", clones.clone(0))[0]
>>> def extensive_moment(window, state):
              return fst.ExtensiveMoments(clones.clone(window).analyze(extmom_index).analyze(state))
>>> extensive_moment(0,0).moments(0,0,0,0,0).num_values()
4.66104911131207e-310

The num_values() should be identical to the .sum_dble() when called on moments(0,0,0,0,0).
>>> extensive_moment(0,0).moments(0,0,0,0,0).sum_dble() 
416627.0
hhatch commented 3 years ago

Thank you for this well organized example. Indeed there is something really wrong here and I haven't figured it out yet. One thing I did find is that a cc implementation of the above python analysis returns a reasonable Accumulator::num_values(). I will continue searching for a fix for the python interface.

#include <iostream>                                                               
#include "feasst.h"                                                               

// Manually add CheckRigidBonds to deserialize_map. Shouldn't have to...          
feasst::CheckRigidBonds check;                                                    

int main() {                                                                      
  std::cout << feasst::version() << std::endl;                                    
  auto clones = feasst::MakeClones("checkpoint", 1);                              
  const int extmom_index =                                                        
    feasst::SeekAnalyze().index("ExtensiveMoments", clones->clone(0))[0];         
  feasst::ExtensiveMoments extensive_moment(clones->clone(0).analyze(extmom_index).analyze(0));
  std::cout << extensive_moment.moments(0,0,0,0,0).num_values() << std::endl;     
  std::cout << extensive_moment.moments(0,0,0,0,0).sum_dble() << std::endl;
}

Note: this brings up a secondary issue. For some reason, Analyze::deserialize_map() doesn't populate every Analyze and I had to manually initialize feasst::CheckRigidBonds check to deserialize the checkpoint file.