Scine::Readuct  6.0.0
This is the SCINE module ReaDuct.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Pages
NtOptimization2Task.h
Go to the documentation of this file.
1 
7 #ifndef READUCT_NTOPTIMIZATION2TASK_H_
8 #define READUCT_NTOPTIMIZATION2TASK_H_
9 
10 /* Readuct */
11 #include "Tasks/Task.h"
12 /* Scine */
15 #include <Utils/IO/Yaml.h>
16 /* External */
17 #include <boost/exception/diagnostic_information.hpp>
18 #include <boost/filesystem.hpp>
19 /* std c++ */
20 #include <cstdio>
21 #include <fstream>
22 #include <iostream>
23 
24 namespace Scine {
25 namespace Readuct {
26 
27 class NtOptimization2Task : public Task {
28  public:
35  NtOptimization2Task(std::vector<std::string> input, std::vector<std::string> output, std::shared_ptr<Core::Log> logger = nullptr)
36  : Task(std::move(input), std::move(output), std::move(logger)) {
37  }
38 
39  std::string name() const override {
40  return "NT2 Optimization";
41  }
42 
43  bool run(SystemsMap& systems, Utils::UniversalSettings::ValueCollection taskSettings, bool testRunOnly = false,
44  std::vector<std::function<void(const int&, const Utils::AtomCollection&, const Utils::Results&, const std::string&)>>
45  observers = {}) const final {
48 
49  bool silentCalculator = taskSettings.extract("silent_stdout_calculator", true);
50  std::shared_ptr<Core::Calculator> calc;
51  if (!testRunOnly) { // leave out in case of task chaining --> attention calc is NULL
52  // Note: _input is guaranteed not to be empty by Task constructor
53  calc = copyCalculator(systems, _input.front(), name());
54  Utils::CalculationRoutines::setLog(*calc, true, true, !silentCalculator);
55 
56  // Check system size
57  if (calc->getStructure()->size() == 1) {
58  throw std::runtime_error("Cannot calculate NT2 optimization for monoatomic systems.");
59  }
60  }
61 
62  // Generate optimizer
63  auto optimizer = std::make_unique<Utils::NtOptimizer2>(*calc);
64  // Read and delete special settings
65  bool stopOnError = stopOnErrorExtraction(taskSettings);
66  // Apply user settings
67  auto settings = optimizer->getSettings();
68  settings.merge(taskSettings);
69  if (!settings.valid()) {
70  settings.throwIncorrectSettings();
71  }
72  optimizer->setSettings(settings);
73 
74  // If no errors encountered until here, the basic settings should be alright
75  if (testRunOnly) {
76  return true;
77  }
78 
79  // Add observer
80  // Trajectory stream
81  using Writer = Utils::XyzStreamHandler;
82  const std::string& outputSystem = ((!_output.empty()) ? _output[0] : _input[0]);
83  boost::filesystem::path dir(outputSystem);
84  boost::filesystem::create_directory(dir);
85  boost::filesystem::path trjfile(outputSystem + ".nt.trj.xyz");
86  std::ofstream trajectory((dir / trjfile).string(), std::ofstream::out);
87  double oldEnergy = 0.0;
88  Eigen::VectorXd oldParams;
89  auto func = [&](const int& cycle, const double& energy, const Eigen::VectorXd& params) {
90  if (oldParams.size() != params.size()) {
91  oldParams = params;
92  }
93  if (cycle == 1) {
94  printf("%7s %16s %16s %16s %16s\n", "Cycle", "Energy", "Energy Diff.", "Step RMS", "Max. Step");
95  }
96  auto diff = (params - oldParams).eval();
97  printf("%7d %+16.9f %+16.9f %+16.9f %+16.9f\n", cycle, energy, energy - oldEnergy,
98  sqrt(diff.squaredNorm() / diff.size()), diff.cwiseAbs().maxCoeff());
99  oldEnergy = energy;
100  oldParams = params;
101  auto structure = calc->getStructure();
102  Writer::write(trajectory, *structure, std::to_string(energy));
103  };
104  optimizer->addObserver(func);
105  auto cout = _logger->output;
106 
107  // Add custom observers
108  auto customObservers = [&calc, &observers](const int& cycle, const double& /*energy*/, const Eigen::VectorXd& /*params*/) {
109  for (auto& observer : observers) {
110  auto atoms = calc->getStructure();
111  Utils::Results& results = calc->results();
112  observer(cycle, *atoms, results, "nt2_scan");
113  }
114  };
115  optimizer->addObserver(customObservers);
116 
117  // Run optimization
118  auto structure = calc->getStructure();
119  int cycles = 0;
120  try {
121  cycles = optimizer->optimize(*structure, *_logger);
122  }
123  catch (...) {
124  trajectory.close();
125  if (stopOnError) {
126  boost::filesystem::path failedPath(outputSystem + ".nt.failed.xyz");
127  std::ofstream failedFile((dir / failedPath).string(), std::ofstream::out);
128  Writer::write(failedFile, *structure);
129  failedFile.close();
130  throw;
131  }
132  // check if thrown exception corresponds to no actual calculation failures but simply no guess found
133  // --> no thrown exception intended in this case if stopOnError equals false
134  std::string noGuess = "No transition state guess was found in Newton Trajectory scan";
135  size_t found = boost::current_exception_diagnostic_information().find(noGuess);
136  if (found == std::string::npos) { // did not find harmless exception -> throw
137  boost::filesystem::path failedPath(outputSystem + ".nt.failed.xyz");
138  std::ofstream failedFile((dir / failedPath).string(), std::ofstream::out);
139  Writer::write(failedFile, *structure);
140  failedFile.close();
141  throw;
142  }
143  cout << Core::Log::endl << " No TS guess found in NT2 scan." << Core::Log::endl << Core::Log::endl;
144  return false;
145  }
146 
147  trajectory.close();
148 
149  cout << Core::Log::endl
150  << " Found TS guess after " << cycles << " iterations." << Core::Log::endl
151  << Core::Log::endl;
152 
153  // Print/Store results
154  systems[outputSystem] = calc;
155  boost::filesystem::path xyzfile(outputSystem + ".xyz");
156  std::ofstream xyz((dir / xyzfile).string(), std::ofstream::out);
157  Writer::write(xyz, *(calc->getStructure()));
158  xyz.close();
159 
160  return true;
161  }
162 };
163 
164 } // namespace Readuct
165 } // namespace Scine
166 
167 #endif // READUCT_NTOPTIMIZATION2TASK_H_
void warningIfMultipleOutputsGiven() const
Warn if more than one output system was specified.
Definition: Task.h:107
Definition: NtOptimization2Task.h:27
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: NtOptimization2Task.h:43
NtOptimization2Task(std::vector< std::string > input, std::vector< std::string > output, std::shared_ptr< Core::Log > logger=nullptr)
Construct a new NtOptimization2Task.
Definition: NtOptimization2Task.h:35
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
std::string name() const override
Getter for the tasks name.
Definition: NtOptimization2Task.h:39