Scine::Readuct  6.0.0
This is the SCINE module ReaDuct.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Pages
AfirOptimizationTask.h
Go to the documentation of this file.
1 
7 #ifndef READUCT_AFIROPTIMIZATIONTASK_H_
8 #define READUCT_AFIROPTIMIZATIONTASK_H_
9 
10 /* Readuct */
11 #include "Tasks/Task.h"
12 /* Scine */
20 /* External */
21 #include <boost/exception/diagnostic_information.hpp>
22 #include <boost/filesystem.hpp>
23 /* std c++ */
24 #include <cstdio>
25 #include <fstream>
26 #include <iostream>
27 
28 namespace Scine {
29 namespace Readuct {
30 
31 class AfirOptimizationTask : public Task {
32  public:
39  AfirOptimizationTask(std::vector<std::string> input, std::vector<std::string> output, std::shared_ptr<Core::Log> logger = nullptr)
40  : Task(std::move(input), std::move(output), std::move(logger)) {
41  }
42 
43  std::string name() const override {
44  return "AFIR Optimization";
45  }
46 
47  bool run(SystemsMap& systems, Utils::UniversalSettings::ValueCollection taskSettings, bool testRunOnly = false,
48  std::vector<std::function<void(const int&, const Utils::AtomCollection&, const Utils::Results&, const std::string&)>>
49  observers = {}) const final {
52 
53  bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true);
54  std::shared_ptr<Core::Calculator> calc;
55  if (!testRunOnly) { // leave out in case of task chaining --> attention calc is NULL
56  // Note: _input is guaranteed not to be empty by Task constructor
57  calc = copyCalculator(systems, _input.front(), name());
58  Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator);
59 
60  // Check system size
61  if (calc->getStructure()->size() == 1) {
62  throw std::runtime_error("Cannot calculate AFIR for monoatomic systems.");
63  }
64  }
65 
66  // Generate optimizer
67  auto optimizertype = taskSettings.extract("optimizer", std::string{"BFGS"});
68  std::transform(optimizertype.begin(), optimizertype.end(), optimizertype.begin(), ::toupper);
69  std::shared_ptr<Utils::AfirOptimizerBase> optimizer;
70  if (optimizertype == "LBFGS") {
71  auto tmp = std::make_shared<Utils::AfirOptimizer<Utils::Lbfgs>>(*calc);
72  tmp->optimizer.useTrustRadius = true;
73  tmp->optimizer.trustRadius = 0.1;
74  optimizer = std::move(tmp);
75  }
76  else if (optimizertype == "BFGS") {
77  auto tmp = std::make_shared<Utils::AfirOptimizer<Utils::Bfgs>>(*calc);
78  tmp->optimizer.useTrustRadius = true;
79  tmp->optimizer.trustRadius = 0.1;
80  optimizer = std::move(tmp);
81  }
82  else if (optimizertype == "SD" || optimizertype == "STEEPESTDESCENT") {
83  auto tmp = std::make_shared<Utils::AfirOptimizer<Utils::SteepestDescent>>(*calc);
84  optimizer = std::move(tmp);
85  }
86  else {
87  throw std::runtime_error("Unknown Optimizer requested for AFIR, available are: SD, BFGS and LBFGS!");
88  }
89 
90  // Read and delete special settings
91  bool stopOnError = stopOnErrorExtraction(taskSettings);
92 
93  // Apply user settings
94  auto settings = optimizer->getSettings();
95  settings.merge(taskSettings);
96  if (!settings.valid()) {
97  settings.throwIncorrectSettings();
98  }
99  optimizer->setSettings(settings);
100  if (!testRunOnly && !Utils::GeometryOptimization::settingsMakeSense(*optimizer)) {
101  throw std::logic_error("The given calculator settings are too inaccurate for the given convergence criteria of "
102  "this optimization Task");
103  }
104 
105  // If no errors encountered until here, the basic settings should be alright
106  if (testRunOnly) {
107  return true;
108  }
109 
110  // Add observer
111  // Trajectory stream
112  using Writer = Utils::XyzStreamHandler;
113  const std::string& outputSystem = (!_output.empty() ? _output.front() : _input.front());
114  boost::filesystem::path dir(outputSystem);
115  boost::filesystem::create_directory(dir);
116  boost::filesystem::path trjfile(outputSystem + ".afir.trj.xyz");
117  std::ofstream trajectory((dir / trjfile).string(), std::ofstream::out);
118  double oldEnergy = 0.0;
119  Eigen::VectorXd oldParams;
120  auto cout = _logger->output;
121  auto func = [&](const int& cycle, const double& energy, const Eigen::VectorXd& params) {
122  if (oldParams.size() != params.size()) {
123  oldParams = params;
124  }
125  if (cycle == 1) {
126  cout.printf("%7s %16s %16s %16s %16s\n", "Cycle", "Energy", "Energy Diff.", "Step RMS", "Max. Step");
127  }
128  auto diff = (params - oldParams).eval();
129  cout.printf("%7d %+16.9f %+16.9f %+16.9f %+16.9f\n", cycle, energy, energy - oldEnergy,
130  sqrt(diff.squaredNorm() / diff.size()), diff.cwiseAbs().maxCoeff());
131  oldEnergy = energy;
132  oldParams = params;
133  auto structure = calc->getStructure();
134  Writer::write(trajectory, *structure, std::to_string(energy));
135  };
136  optimizer->addObserver(func);
137 
138  // Add custom observers
139  auto customObservers = [&calc, &observers](const int& cycle, const double& /*energy*/, const Eigen::VectorXd& /*params*/) {
140  for (auto& observer : observers) {
141  auto atoms = calc->getStructure();
142  Utils::Results& results = calc->results();
143  observer(cycle, *atoms, results, "afir_scan");
144  }
145  };
146  optimizer->addObserver(customObservers);
147 
148  // Run optimization
149  auto structure = calc->getStructure();
150  int cycles = 0;
151  try {
152  cycles = optimizer->optimize(*structure, *_logger);
153  }
154  catch (...) {
155  Writer::write(trajectory, *calc->getStructure());
156  trajectory.close();
157  _logger->error << "AFIR Optimization failed with error!" << Core::Log::endl;
158  if (stopOnError) {
159  throw;
160  }
161  _logger->error << boost::current_exception_diagnostic_information() << Core::Log::endl;
162  return false;
163  }
164  trajectory.close();
165 
166  int maxiter = settings.getInt("convergence_max_iterations");
167  if (maxiter > cycles) {
168  cout << Core::Log::endl
169  << " Converged after " << cycles << " iterations." << Core::Log::endl
170  << Core::Log::endl;
171  }
172  else {
173  cout << Core::Log::endl
174  << " Stopped after " << maxiter << " iterations." << Core::Log::endl
175  << Core::Log::endl;
176  if (stopOnError) {
177  throw std::runtime_error("Problem: AFIR optimization did not converge.");
178  }
179  }
180 
181  // Print/Store results
182  boost::filesystem::path xyzfile(outputSystem + ".xyz");
183  std::ofstream xyz((dir / xyzfile).string(), std::ofstream::out);
184  Writer::write(xyz, *(calc->getStructure()));
185  xyz.close();
186  systems[outputSystem] = std::move(calc);
187 
188  return cycles < maxiter;
189  }
190 };
191 
192 } // namespace Readuct
193 } // namespace Scine
194 
195 #endif // READUCT_AFIROPTIMIZATIONTASK_H_
std::string name() const override
Getter for the tasks name.
Definition: AfirOptimizationTask.h:43
void warningIfMultipleOutputsGiven() const
Warn if more than one output system was specified.
Definition: Task.h:107
AfirOptimizationTask(std::vector< std::string > input, std::vector< std::string > output, std::shared_ptr< Core::Log > logger=nullptr)
Construct a new AfirOptimizationTask.
Definition: AfirOptimizationTask.h:39
static std::ostream & endl(std::ostream &os)
const std::vector< std::string > & input() const
Getter for the expected names of the input systems.
Definition: Task.h:85
bool run(SystemsMap &systems, Utils::UniversalSettings::ValueCollection taskSettings, bool testRunOnly=false, std::vector< std::function< void(const int &, const Utils::AtomCollection &, const Utils::Results &, const std::string &)>> observers={}) const final
Executes the actual task represented by this class.
Definition: AfirOptimizationTask.h:47
Definition: AfirOptimizationTask.h:31
void warningIfMultipleInputsGiven() const
Warn if more than one input system was specified.
Definition: Task.h:98
const std::vector< std::string > & output() const
Getter for the names of the output systems generated by this task.
Definition: Task.h:92
The base class for all tasks in Readuct.
Definition: Task.h:29