SafeHyperJob.cpp

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** SPDX-License-Identifier: BSD-3-Clause
**
** Copyright (c) 2010--present, eOn Development Team
** All rights reserved.
**
** Repo:
** https://github.com/TheochemUI/eOn
*/
#include "eon/SafeHyperJob.h"
#include "eon/BondBoost.h"
#include "eon/Dynamics.h"
#include "eon/ForceCallTimer.h"
 
#include <cmath>
 
namespace eonc {
 
void SafeHyperJob::reportResults() {
  if (newStateFlag) {
    QUILL_LOG_DEBUG(log, "Transition time: {:.2e} s",
                    minCorrectedTime * 1.0e-15 * params.constants().timeUnit);
  } else {
    QUILL_LOG_DEBUG(log,
                    "No new state was found in {} dynamics steps ({:.3e} s)",
                    params.dynamics_options().steps,
                    time * 1.0e-15 * params.constants().timeUnit);
  }
}
 
int SafeHyperJob::dynamics() {
  bool transitionFlag = false, recordFlag = true, stopFlag = false,
       firstTransitFlag = false;
  long nFreeCoord = reactant->numberOfFreeAtoms() * 3;
  long mdBufferLength;
  long step = 0, refineStep, newStateStep = 0;
  long nCheck = 0, nRecord = 0, nBoost = 0, nState = 0;
  long StateCheckInterval, RecordInterval;
  double kinE, kinT, avgT, varT;
  double kB = params.constants().kB;
  double correctedTime = 0.0, sumCorrectedTime = 0.0, firstTransitionTime = 0.0;
  double Temp = 0.0, sumT = 0.0, sumT2 = 0.0;
  double sumboost = 0.0, boost = 1.0, boostPotential = 0.0;
  double transitionTime_current = 0.0, transitionTime_pre = 0.0;
  AtomMatrix velocity;
 
  minCorrectedTime = 1.0e200;
  const auto clock = prdClock();
  StateCheckInterval = clock.state_check;
  RecordInterval = clock.record;
  Temp = params.main_options().temperature;
  newStateFlag = metaStateFlag = false;
 
  mdBufferLength = clock.buffer;
  std::vector<std::shared_ptr<Matter>> mdBuffer(mdBufferLength);
  for (long i = 0; i < mdBufferLength; i++) {
    mdBuffer[i] = std::make_shared<Matter>(pot, params);
  }
  timeBuffer.resize(mdBufferLength);
  biasBuffer.resize(mdBufferLength);
 
  Dynamics safeHyper(current.get(), params);
  BondBoost bondBoost(current.get(), params);
 
  if (params.hyperdynamics_options().bias_potential ==
      Hyperdynamics::BOND_BOOST) {
    bondBoost.initialize();
    current->setBiasPotential(&bondBoost);
  }
 
  safeHyper.setThermalVelocity();
 
  {
    eonc::ForceCallTimer timer(dephaseFCalls);
    dephase();
  }
 
  QUILL_LOG_DEBUG(
      log,
      "Starting MD run\nTemperature: {:.2f} Kelvin\n"
      "Total Simulation Time: {:.2f} fs\nTime Step: {:.2f} fs\nTotal Steps: {}",
      Temp,
      params.dynamics_options().steps * params.dynamics_options().time_step *
          params.constants().timeUnit,
      params.dynamics_options().time_step * params.constants().timeUnit,
      params.dynamics_options().steps);
  QUILL_LOG_DEBUG(log, "MD buffer length: {}", mdBufferLength);
 
  long tenthSteps = params.dynamics_options().steps / 10;
  if (tenthSteps == 0) {
    tenthSteps = params.dynamics_options().steps;
  }
 
  while (!stopFlag) {
    boost = 1.0;
    boostPotential = 0.0;
    if ((params.hyperdynamics_options().bias_potential ==
         Hyperdynamics::BOND_BOOST) &&
        !newStateFlag) {
      bondBoost.advance();
      boostPotential = bondBoost.boost();
      QUILL_LOG_TRACE_L1(log, "step= {} , boost = {:.5f}", step,
                         boostPotential);
      if (Temp > 0.0 && kB > 0.0) {
        boost = std::exp(boostPotential / kB / Temp);
      } else {
        boost = 1.0;
      }
      if (boost > 1.0) {
        sumboost += boost;
        nBoost++;
      }
    }
    time += params.dynamics_options().time_step * boost;
 
    kinE = current->getKineticEnergy();
    kinT = (2.0 * kinE / nFreeCoord / kB);
    sumT += kinT;
    sumT2 += kinT * kinT;
    QUILL_LOG_TRACE_L1(log, "steps = {:10} temp = {:10.5f}", step, kinT);
 
    safeHyper.oneStep();
    mdFCalls++;
 
    nCheck++;
    step++;
    QUILL_LOG_TRACE_L1(log, "step = {:4}, time = {:10.4f}", step, time);
 
    if (params.parallel_replica_options().refine_transition && recordFlag &&
        !newStateFlag) {
      if (nCheck % RecordInterval == 0) {
        *mdBuffer[nRecord] = *current;
        timeBuffer[nRecord] = time;
        biasBuffer[nRecord] = boostPotential;
        nRecord++;
      }
    }
 
    if ((nCheck == StateCheckInterval) && !newStateFlag) {
      nCheck = 0;
      nRecord = 0;
      {
        eonc::ForceCallTimer timer(minimizeFCalls);
        transitionFlag = checkState(current.get(), reactant.get());
      }
      if (transitionFlag) {
        nState++;
        QUILL_LOG_DEBUG(log, "New State {}: ", nState);
        *finalStateTmp = *current;
        transitionTime = time;
        newStateStep = step;
        transitionStep = newStateStep;
        firstTransitFlag = 1;
      }
    }
 
    if (transitionFlag) {
      QUILL_LOG_TRACE_L1(log, "Refining transition time.");
      const bool can_refine =
          params.parallel_replica_options().refine_transition && nRecord >= 2;
      if (can_refine) {
        eonc::ForceCallTimer timer(refineFCalls);
        refineStep = refine(mdBuffer, reactant.get());
        transitionStep =
            newStateStep - StateCheckInterval + refineStep * RecordInterval;
        transitionTime_current = timeBuffer[static_cast<size_t>(refineStep)];
        transitionPot = biasBuffer[static_cast<size_t>(refineStep)];
        const long prev = refineStep > 0 ? refineStep - 1 : 0;
        *current = *mdBuffer[static_cast<size_t>(prev)];
      } else {
        refineStep = 0;
        transitionTime_current = time;
        transitionPot = boostPotential;
      }
      transitionTime = transitionTime_current - transitionTime_pre;
      transitionTime_pre = transitionTime_current;
      correctedTime =
          (Temp > 0.0 && kB > 0.0)
              ? transitionTime * std::exp((-1) * transitionPot / kB / Temp)
              : transitionTime;
      sumCorrectedTime += correctedTime;
      if (nState == 1) {
        firstTransitionTime = transitionTime;
      }
      velocity = current->getVelocities();
      velocity = velocity * (-1);
      current->setVelocities(velocity);
 
      if (correctedTime < minCorrectedTime) {
        minCorrectedTime = correctedTime;
        if (can_refine) {
          *saddle = *mdBuffer[static_cast<size_t>(refineStep)];
        } else {
          *saddle = *current;
        }
        *finalState = *finalStateTmp;
      }
      QUILL_LOG_DEBUG(log,
                      "tranisitonTime= {:.3e} s, biasPot= {:.3f} eV, "
                      "correctedTime= {:.3e} s, "
                      "sumCorrectedTime= {:.3e} s, minCorTime= {:.3e} s",
                      transitionTime * 1e-15 * params.constants().timeUnit,
                      transitionPot,
                      correctedTime * 1e-15 * params.constants().timeUnit,
                      sumCorrectedTime * 1e-15 * params.constants().timeUnit,
                      minCorrectedTime * 1.0e-15 * params.constants().timeUnit);
 
      transitionFlag = false;
    }
 
    if (firstTransitFlag && sumCorrectedTime > firstTransitionTime) {
      stopFlag = true;
      newStateFlag = true;
    }
 
    if ((step % tenthSteps == 0) || (step == params.dynamics_options().steps)) {
      double maxAtomDistance = current->perAtomNorm(*reactant);
      QUILL_LOG_DEBUG(
          log, "progress: {:.0f}%, max displacement: {:.3f}, step {} / {}",
          static_cast<double>(100.0 * step / params.dynamics_options().steps),
          maxAtomDistance, step, params.dynamics_options().steps);
    }
 
    // Honor dynamics step budget (parity with TADJob); without this the
    // loop only exits on a confidence-gated transition and can run forever.
    if (step >= params.dynamics_options().steps) {
      stopFlag = true;
    }
  }
 
  avgT = sumT / step;
  varT = sumT2 / step - avgT * avgT;
 
  if (nBoost > 0) {
    QUILL_LOG_DEBUG(log,
                    "Temperature : Average = {:.6f} ; Stddev = {:.6f} ; "
                    "Factor = {:.6f}; Boost = {:.6f}",
                    avgT, std::sqrt(varT), varT / avgT / avgT * nFreeCoord / 2,
                    sumboost / nBoost);
  } else {
    QUILL_LOG_DEBUG(
        log,
        "Temperature : Average = {:.6f} ; Stddev = {:.6f} ; Factor = {:.6f}",
        avgT, std::sqrt(varT), varT / avgT / avgT * nFreeCoord / 2);
  }
  if (std::isfinite(avgT) == 0) {
    QUILL_LOG_DEBUG(log, "Infinite average temperature, something went wrong!");
    newStateFlag = false;
  }
 
  // finalState is only filled on transition; keep product valid otherwise.
  current->setBiasPotential(nullptr);
 
  if (newStateFlag && finalState) {
    *product = *finalState;
  } else if (current) {
    *product = *current;
  }
 
  if (newStateFlag) {
    return 1;
  } else {
    return 0;
  }
}
 
} // namespace eonc