Loading...
Searching...
No Matches
eonc::Lanczos Class Reference

#include <Lanczos.h>

Inheritance diagram for eonc::Lanczos:

Public Member Functions

 Lanczos (std::shared_ptr< Matter > matter, const Parameters &params, std::shared_ptr< Potential > pot)
 ~Lanczos ()=default
void compute (std::shared_ptr< Matter > matter, AtomMatrix initialDirection)
void compute (std::shared_ptr< Matter > matter, AtomMatrix initialDirection, const VectorXi &mobileAtoms)
 Same as compute(matter, dir) but with an explicit mobile atom list (intersected with free flags).
double getEigenvalue ()
AtomMatrix getEigenvector ()
Public Member Functions inherited from eonc::LowestEigenmode
 LowestEigenmode (std::shared_ptr< Potential > potPassed, const Parameters &parameters)
 ~LowestEigenmode ()=default

Private Attributes

AtomMatrix lowestEv
double lowestEw
eonc::log::Scoped log

Additional Inherited Members

Public Attributes inherited from eonc::LowestEigenmode
long totalForceCalls {0}
double statsTorque {0.0}
double statsCurvature {0.0}
double statsAngle {0.0}
long statsRotations {0}
long totalIterations {0}
Static Public Attributes inherited from eonc::LowestEigenmode
static const char MINMODE_DIMER [] = "dimer"
static const char MINMODE_GPRDIMER [] = "gprdimer"
static const char MINMODE_LANCZOS [] = "lanczos"
static const char MINMODE_DAVIDSON [] = "davidson"
Protected Attributes inherited from eonc::LowestEigenmode
std::shared_ptr< Potential > pot
const Parameters & params

Detailed Description

Definition at line 26 of file Lanczos.h.

Constructor & Destructor Documentation

◆ Lanczos()

Lanczos::Lanczos ( std::shared_ptr< Matter > matter,
const Parameters & params,
std::shared_ptr< Potential > pot )

Definition at line 32 of file Lanczos.cpp.

34 : LowestEigenmode(pot, params) {
35 lowestEv.resize(matter->numberOfAtoms(), 3);
36 lowestEv.setZero();
37 lowestEw = 0.0;
38}
eonc::Lanczos::lowestEv
AtomMatrix lowestEv
Definition Lanczos.h:41
eonc::Lanczos::lowestEw
double lowestEw
Definition Lanczos.h:42
eonc::LowestEigenmode::params
const Parameters & params
Definition LowestEigenmode.h:28
eonc::LowestEigenmode::pot
std::shared_ptr< Potential > pot
Definition LowestEigenmode.h:27
eonc::LowestEigenmode::LowestEigenmode
LowestEigenmode(std::shared_ptr< Potential > potPassed, const Parameters &parameters)
Definition LowestEigenmode.h:42

◆ ~Lanczos()

eonc::Lanczos::~Lanczos ( )
default

Member Function Documentation

◆ compute() [1/2]

void Lanczos::compute ( std::shared_ptr< Matter > matter,
AtomMatrix initialDirection )

Definition at line 40 of file Lanczos.cpp.

40 {
41 const VectorXi mobile =
42 resolveMobileAtoms(matter.get(), params.lanczos_options.phva_atoms);
43 compute(std::move(matter), std::move(direction), mobile);
44}
eonc::Lanczos::compute
void compute(std::shared_ptr< Matter > matter, AtomMatrix initialDirection)
Definition Lanczos.cpp:40
eonc::resolveMobileAtoms
VectorXi resolveMobileAtoms(const Matter *matter, const std::string &atomList)
PHVA-class mobile set for FD Hessian and matrix-free Krylov (Lanczos / Davidson).
Definition MobileAtoms.cpp:64

◆ compute() [2/2]

void Lanczos::compute ( std::shared_ptr< Matter > matter,
AtomMatrix initialDirection,
const VectorXi & mobileAtoms )

Same as compute(matter, dir) but with an explicit mobile atom list (intersected with free flags).

Eigenvector rows outside the list are 0.

Definition at line 46 of file Lanczos.cpp.

47 {
48 totalForceCalls = 0;
49 statsRotations = 0;
50 lowestEv.resize(matter->numberOfAtoms(), 3);
51 lowestEv.setZero();
52
53 const VectorXi mobile = resolveMobileAtoms(matter.get(), mobileIn);
54 const int size = 3 * static_cast<int>(mobile.size());
55 if (size == 0) {
56 lowestEw = 0.0;
57 return;
58 }
59
60 const long maxIters = params.lanczos_options.max_iterations;
61 MatrixXd T(size, maxIters), Q(size, maxIters);
62 T.setZero();
63 VectorXd u(size), r = packMobileRows(direction, mobile);
64
65 double alpha, beta = r.norm();
66 if (beta < eonc::safemath::eps) {
67 lowestEw = 0.0;
68 return;
69 }
70 double ew = 0, ewOld = 0, ewAbsRelErr;
71 const double dr = params.main_options.finiteDifference;
72 VectorXd evEst, evT, evOldEst;
73
74 auto tmpMatter = std::make_unique<Matter>(*matter);
75 const long forceCallsStart = tmpMatter->getForceCalls();
76 const AtomMatrix pos0 = matter->getPositions();
77 const VectorXd force0 = mobileForces(tmpMatter.get(), mobile);
78
79 auto applyH = [&](const VectorXd &v) -> VectorXd {
80 AtomMatrix pos = pos0;
81 unpackMobileRows(packMobileRows(pos0, mobile) + dr * v, mobile, pos);
82 tmpMatter->setPositions(pos);
83 return -(mobileForces(tmpMatter.get(), mobile) - force0) / dr;
84 };
85
86 for (int i = 0; i < size; i++) {
87 statsRotations = i;
88 Q.col(i) = r / beta;
89
90 u = applyH(Q.col(i));
91
92 if (i == 0) {
93 r = u;
94 } else {
95 r = u - beta * Q.col(i - 1);
96 }
97 alpha = Q.col(i).dot(r);
98 r = r - alpha * Q.col(i);
99
100 T(i, i) = alpha;
101 if (i > 0) {
102 T(i - 1, i) = beta;
103 T(i, i - 1) = beta;
104 }
105
106 beta = r.norm();
107
108 if (beta <= 1e-10 * std::fabs(alpha)) {
109 if (i == 0) {
110 ew = alpha;
111 evEst = Q.col(0);
112 }
113 QUILL_LOG_ERROR(log, "[ILanczos] ERROR: linear dependence");
114 break;
115 }
116 if (i >= 1) {
117 Eigen::SelfAdjointEigenSolver<MatrixXd> es(T.block(0, 0, i + 1, i + 1));
118 ew = es.eigenvalues()(0);
119 evT = es.eigenvectors().col(0);
120 ewAbsRelErr = eonc::safemath::safe_div(std::fabs(ew - ewOld),
121 std::fabs(ewOld), 1.0);
122 ewOld = ew;
123
124 evEst = Q.block(0, 0, size, i + 1) * evT;
125 evEst.normalize();
126 statsAngle = eonc::safemath::safe_acos(std::fabs(evEst.dot(evOldEst))) *
127 (180 / eonc::helpers::pi);
128 statsTorque = ewAbsRelErr;
129 evOldEst = evEst;
130 QUILL_LOG_INFO(log,
131 "[ILanczos] {:9s} {:9s} {:10s} {:14s} {:9.4f} "
132 "{:10.6f} {:7.3f} {:5} n_mobile={}",
133 "----", "----", "----", "----", ew, ewAbsRelErr,
134 statsAngle, i, mobile.size());
135 if (ewAbsRelErr < params.lanczos_options.tolerance) {
136 QUILL_LOG_INFO(log, "[ILanczos] Tolerance reached: {}",
137 params.lanczos_options.tolerance);
138 break;
139 }
140 } else {
141 ew = alpha;
142 ewOld = ew;
143 evEst = Q.col(0);
144 evOldEst = Q.col(0);
145 if (lowestEw != 0.0 && params.lanczos_options.quit_early) {
146 double Cprev = lowestEw;
147 double Cnew = u.dot(Q.col(i));
148 ewAbsRelErr = eonc::safemath::safe_div(std::fabs(Cnew - Cprev),
149 std::fabs(Cprev), 1.0);
150 if (ewAbsRelErr <= params.lanczos_options.tolerance) {
151 statsAngle = 0.0;
152 statsTorque = ewAbsRelErr;
153 QUILL_LOG_INFO(log, "[ILanczos] Tolerance reached: {}",
154 params.lanczos_options.tolerance);
155 break;
156 }
157 }
158 }
159
160 if (i >= params.lanczos_options.max_iterations - 1) {
161 QUILL_LOG_ERROR(log, "[ILanczos] Max iterations");
162 break;
163 }
164 }
165
166 lowestEw = ew;
167 totalForceCalls = tmpMatter->getForceCalls() - forceCallsStart;
168
169 lowestEv.setZero();
170 if (evEst.size() == size) {
171 unpackMobileRows(evEst, mobile, lowestEv);
172 }
173}
MatrixXd
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, eOnStorageOrder > MatrixXd
Definition Eigen.h:33
AtomMatrix
Eigen::Matrix< double, Eigen::Dynamic, 3, eOnStorageOrder > AtomMatrix
Definition Eigen.h:37
eonc::Lanczos::log
eonc::log::Scoped log
Definition Lanczos.h:43
eonc::helpers::pi
constexpr double pi
Definition HelperFunctions.h:28
eonc::safemath::safe_div
constexpr double safe_div(double num, double denom, double fallback=0.0)
Definition SafeMath.h:21
eonc::safemath::safe_acos
double safe_acos(double x)
Definition SafeMath.h:34
eonc::safemath::eps
constexpr double eps
Definition SafeMath.h:19
eonc::unpackMobileRows
void unpackMobileRows(const VectorXd &packed, const VectorXi &mobile, AtomMatrix &full)
Write a 3*n_mobile vector into full AtomMatrix rows (other rows unchanged).
Definition MobileAtoms.cpp:125
eonc::packMobileRows
VectorXd packMobileRows(const AtomMatrix &full, const VectorXi &mobile)
Pack full (n_atoms,3) rows of mobile atoms into a 3*n_mobile vector.
Definition MobileAtoms.cpp:117
eonc::mobileForces
VectorXd mobileForces(Matter *matter, const VectorXi &mobile)
Force components on mobile atoms after Matter has a valid force cache (calls getForces under the hood...
Definition MobileAtoms.cpp:132

◆ getEigenvalue()

double Lanczos::getEigenvalue ( )

Definition at line 175 of file Lanczos.cpp.

175{ return lowestEw; }

◆ getEigenvector()

AtomMatrix Lanczos::getEigenvector ( )

Definition at line 177 of file Lanczos.cpp.

177{ return lowestEv; }

Member Data Documentation

◆ log

eonc::log::Scoped eonc::Lanczos::log
private

Definition at line 43 of file Lanczos.h.

◆ lowestEv

AtomMatrix eonc::Lanczos::lowestEv
private

Definition at line 41 of file Lanczos.h.

◆ lowestEw

double eonc::Lanczos::lowestEw
private

Definition at line 42 of file Lanczos.h.

The documentation for this class was generated from the following files: