Skip to content

Commit 5ff084e

Browse files
author
abacus_fixer
committed
deltaspin: lift 6 lambda-loop helpers out of SpinConstrain class
Phase 1 of god-class slimming: lift six lambda-loop helper methods (check_rms_stop, check_restriction, check_gradient_decay, cal_alpha_opt, print_header, print_termination) out of the SpinConstrain<TK> class template and into free function templates in the spinconstrain namespace. Motivation ---------- spin_constrain.h is a god class with 7+ mixed responsibilities. The lambda-loop workflow helpers are pure algorithms / I/O routines that do not need to be part of the class interface; lifting them shrinks the public API surface and clarifies the boundary between state (SpinConstrain) and algorithm (free functions). Changes ------- - New header lambda_loop_helper.h declares 6 free function templates taking 'const SpinConstrain<TK>&' as first parameter. - lambda_loop_helper.cpp: 6 explicit specializations on SpinConstrain<std::complex<double>> rewritten as free function templates with explicit instantiation for std::complex<double>. - template_helpers.cpp: 6 corresponding no-op stubs for SpinConstrain<double> re-expressed as free function explicit specializations (preserves linker surface for the nspin=2 stub). - spin_constrain.h: removed 6 member-function declarations; added public getters get_Mi() and get_current_sc_thr() so helpers can read internal state without friendship. - spin_constrain.cpp: get_nat/get_ntype/get_decay_grad(itype)/ get_decay_grad() promoted to const (required for const ref parameter); get_decay_grad(itype) uses map::find instead of operator[] to remain const-correct. - lambda_loop.cpp: 7 call sites switched from 'this->helper(...)' to 'helper(*this, ...)'. - test/template_helpers_test.cpp: call sites updated. - test/CMakeLists.txt: MODULE_LCAO_deltaspin_template_helpers now links lambda_loop_helper.cpp and basic_funcs.cpp (transitive deps of the new free function implementations). Verification ------------ - 'make -j 30' in build_max_para_test: clean build, abacus + abacus_basic_para + all MODULE_LCAO_deltaspin_* targets linked. - 'ctest -R deltaspin --output-on-failure': 5/5 tests passed (basic_func, spin_constrain, template_helpers, pw, core). Out of scope ------------ PW-specific methods (cal_mi_pw, update_psi_charge_pw_*, calculate_delta_hcc) and LCAO-specific helpers (cal_mi_lcao, convert, calculate_MW, collect_MW) remain as member functions for now; they will be lifted in subsequent phases.
1 parent ac49785 commit 5ff084e

8 files changed

Lines changed: 336 additions & 192 deletions

File tree

‎source/source_lcao/module_deltaspin/lambda_loop.cpp‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
11
#include "spin_constrain.h"
2-
32
#include <iostream>
43
#include <cmath>
54
#include <chrono>
65
#include <fstream>
76
#include <iomanip>
87

98
#include "basic_funcs.h"
9+
#include "lambda_loop_helper.h"
1010
#include "source_base/constants.h"
1111
#include "source_io/module_parameter/parameter.h"
1212

@@ -106,7 +106,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
106106

107107
double inner_loop_duration = 0.0;
108108

109-
this->print_header(ofs_running);
109+
print_header(*this, ofs_running);
110110

111111
// =============================================================
112112
// MAIN OPTIMIZATION LOOP
@@ -170,7 +170,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
170170
new_spin = this->Mi_;
171171

172172
// Check if gradient dM/dlambda has decayed below threshold
173-
bool GradLessThanBound = this->check_gradient_decay(new_spin, spin, delta_lambda, dnu_last_step, false, ofs_running);
173+
bool GradLessThanBound = check_gradient_decay(*this, new_spin, spin, delta_lambda, dnu_last_step, false, ofs_running);
174174
if (i_step >= this->nsc_min_ && GradLessThanBound)
175175
{
176176
// Gradient has decayed: further optimization yields diminishing returns
@@ -186,7 +186,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
186186
#endif
187187
inner_loop_duration += duration;
188188
ofs_running << " Total TIME(s) = " << inner_loop_duration << std::endl;
189-
this->print_termination(ofs_running);
189+
print_termination(*this, ofs_running);
190190
break;
191191
}
192192
spin = new_spin;
@@ -259,7 +259,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
259259
- iterstart)).count() / static_cast<double>(1e6);
260260
#endif
261261
inner_loop_duration += duration;
262-
if (this->check_rms_stop(outer_step, i_step, rms_error, duration, inner_loop_duration, ofs_running))
262+
if (check_rms_stop(*this, outer_step, i_step, rms_error, duration, inner_loop_duration, ofs_running))
263263
{
264264
// Save RMS for ESolver to display in the SCF iteration table.
265265
this->last_rms_error_ = rms_error;
@@ -316,7 +316,7 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
316316
}
317317

318318
// Cap step size to prevent overshooting
319-
this->check_restriction(search, alpha_trial, ofs_running);
319+
check_restriction(*this, search, alpha_trial, ofs_running);
320320

321321
// =============================================================
322322
// CUMULATIVE STEP UPDATE
@@ -354,8 +354,8 @@ void spinconstrain::SpinConstrain<std::complex<double>>::run_lambda_loop(int out
354354
spin_plus = this->Mi_;
355355

356356
// Find optimal step size via linear interpolation
357-
alpha_opt = this->cal_alpha_opt(spin, spin_plus, alpha_trial);
358-
this->check_restriction(search, alpha_opt, ofs_running);
357+
alpha_opt = cal_alpha_opt(*this, spin, spin_plus, alpha_trial);
358+
check_restriction(*this, search, alpha_opt, ofs_running);
359359

360360
// Correct dnu: dnu += (alpha_opt - alpha_trial) * search
361361
alpha_plus = alpha_opt - alpha_trial;

‎source/source_lcao/module_deltaspin/lambda_loop_helper.cpp‎

Lines changed: 112 additions & 102 deletions
Large diffs are not rendered by default.

‎source/source_lcao/module_deltaspin/lambda_loop_helper.h‎

Lines changed: 139 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,139 @@
1+
#ifndef LAMBDA_LOOP_HELPER_H
2+
#define LAMBDA_LOOP_HELPER_H
3+
4+
#include <ostream>
5+
#include <vector>
6+
7+
#include "source_base/vector3.h"
8+
#include "spin_constrain.h"
9+
10+
/**
11+
* @file lambda_loop_helper.h
12+
* @brief Free-function helpers for the DeltaSpin lambda optimization loop.
13+
*
14+
* @par Background
15+
* Originally these routines were member functions of SpinConstrain<TK>.
16+
* They have been lifted to free functions in the spinconstrain namespace to
17+
* shrink the SpinConstrain god class. Each helper takes the SpinConstrain
18+
* instance as its first parameter and accesses internal state through the
19+
* public getters (get_Mi, get_sc_lambda, get_constrain, ...).
20+
*
21+
* @par Template parameter TK
22+
* - std::complex<double>: full implementation (nspin=2 and nspin=4)
23+
* - double: stub specialization (no-ops) provided elsewhere
24+
*/
25+
26+
namespace spinconstrain
27+
{
28+
29+
/**
30+
* @brief Print final spin and lambda values when the lambda loop terminates.
31+
*
32+
* @param sc SpinConstrain instance (read Mi_ and lambda_)
33+
* @param ofs_running Log output stream
34+
*/
35+
template <typename TK>
36+
void print_termination(const SpinConstrain<TK>& sc, std::ostream& ofs_running);
37+
38+
/**
39+
* @brief Check whether RMS error is below the convergence threshold or the
40+
* maximum number of inner steps has been reached.
41+
*
42+
* @param sc SpinConstrain instance
43+
* @param outer_step Current SCF outer iteration
44+
* @param i_step Current inner lambda step
45+
* @param rms_error Current RMS error of Mi - M_target
46+
* @param duration Wall time for this step (s)
47+
* @param total_duration Cumulative wall time for the inner loop (s)
48+
* @param ofs_running Log output stream
49+
* @return true if the inner loop should terminate, false otherwise
50+
*/
51+
template <typename TK>
52+
bool check_rms_stop(const SpinConstrain<TK>& sc,
53+
int outer_step,
54+
int i_step,
55+
double rms_error,
56+
double duration,
57+
double total_duration,
58+
std::ostream& ofs_running);
59+
60+
/**
61+
* @brief Print header at the start of the lambda optimization loop.
62+
*
63+
* @param sc SpinConstrain instance
64+
* @param ofs_running Log output stream
65+
*/
66+
template <typename TK>
67+
void print_header(const SpinConstrain<TK>& sc, std::ostream& ofs_running);
68+
69+
/**
70+
* @brief Cap the step size to prevent the optimizer from overshooting.
71+
*
72+
* @details If |alpha_trial * max(search)| exceeds restrict_current_, the
73+
* trial step is reduced so that the maximum lambda change per step is
74+
* bounded. alpha_trial is modified in place.
75+
*
76+
* @param sc SpinConstrain instance
77+
* @param search Current search direction (per atom, 3 components)
78+
* @param alpha_trial Trial step size, modified in place if capped
79+
* @param ofs_running Log output stream
80+
*/
81+
template <typename TK>
82+
void check_restriction(const SpinConstrain<TK>& sc,
83+
const std::vector<ModuleBase::Vector3<double>>& search,
84+
double& alpha_trial,
85+
std::ostream& ofs_running);
86+
87+
/**
88+
* @brief Compute the optimal step size via two-point linear interpolation.
89+
*
90+
* @par Algorithm
91+
* alpha_opt = sum_k / sum_k2 * alpha_trial
92+
* where
93+
* sum_k = sum((target - spin) . (spin_plus - spin)) over constrained components
94+
* sum_k2 = sum(|spin - spin_plus|^2) over constrained components
95+
*
96+
* @param sc SpinConstrain instance
97+
* @param spin Mi at current lambda
98+
* @param spin_plus Mi at trial lambda (current + alpha_trial * search)
99+
* @param alpha_trial Current trial step size
100+
* @return Optimal step size; falls back to alpha_trial when sum_k2 ~ 0
101+
*/
102+
template <typename TK>
103+
double cal_alpha_opt(const SpinConstrain<TK>& sc,
104+
std::vector<ModuleBase::Vector3<double>> spin,
105+
std::vector<ModuleBase::Vector3<double>> spin_plus,
106+
const double alpha_trial);
107+
108+
/**
109+
* @brief Check whether the magnetic susceptibility gradient dM/dlambda has
110+
* decayed below the per-atom-type threshold.
111+
*
112+
* @par Algorithm
113+
* 1. Compute spin_change = new_spin - spin
114+
* 2. Compute nu_change = delta_lambda - dnu_last_step
115+
* 3. Build full gradient matrix dM[ia][ic]/dlambda[ja][jc]
116+
* 4. Extract diagonal; pick max abs per atom type
117+
* 5. Return true if max(|diag|) < decay_grad[itype] for any type
118+
*
119+
* @param sc SpinConstrain instance
120+
* @param new_spin Mi at current lambda
121+
* @param spin Mi at previous lambda
122+
* @param delta_lambda Current lambda change
123+
* @param dnu_last_step Previous cumulative step
124+
* @param print Whether to print detailed gradient info
125+
* @param ofs_running Log output stream
126+
* @return true if gradient decayed below threshold for any atom type
127+
*/
128+
template <typename TK>
129+
bool check_gradient_decay(const SpinConstrain<TK>& sc,
130+
std::vector<ModuleBase::Vector3<double>> new_spin,
131+
std::vector<ModuleBase::Vector3<double>> spin,
132+
std::vector<ModuleBase::Vector3<double>> delta_lambda,
133+
std::vector<ModuleBase::Vector3<double>> dnu_last_step,
134+
bool print,
135+
std::ostream& ofs_running);
136+
137+
} // namespace spinconstrain
138+
139+
#endif // LAMBDA_LOOP_HELPER_H

‎source/source_lcao/module_deltaspin/spin_constrain.cpp‎

Lines changed: 21 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -257,10 +257,10 @@ int SpinConstrain<TK>::get_iwt(int itype, int iat, int orbital_index) const
257257

258258
/// @brief Get total number of atoms across all element types
259259
template <typename TK>
260-
int SpinConstrain<TK>::get_nat()
260+
int SpinConstrain<TK>::get_nat() const
261261
{
262262
int nat = 0;
263-
for (std::map<int, int>::iterator it = this->atomCounts.begin(); it != this->atomCounts.end(); ++it)
263+
for (std::map<int, int>::const_iterator it = this->atomCounts.begin(); it != this->atomCounts.end(); ++it)
264264
{
265265
nat += it->second;
266266
}
@@ -269,7 +269,7 @@ int SpinConstrain<TK>::get_nat()
269269

270270
/// @brief Get number of element types
271271
template <typename TK>
272-
int SpinConstrain<TK>::get_ntype()
272+
int SpinConstrain<TK>::get_ntype() const
273273
{
274274
return this->atomCounts.size();
275275
}
@@ -618,9 +618,10 @@ void SpinConstrain<TK>::zero_Mi()
618618
/// this function can only be called by the root process because only
619619
/// root process reads the ScDecayGrad from json file
620620
template <typename TK>
621-
double SpinConstrain<TK>::get_decay_grad(int itype)
621+
double SpinConstrain<TK>::get_decay_grad(int itype) const
622622
{
623-
return this->ScDecayGrad[itype];
623+
std::map<int, double>::const_iterator it = this->ScDecayGrad.find(itype);
624+
return it != this->ScDecayGrad.end() ? it->second : 0.0;
624625
}
625626

626627
/// set grad_decy
@@ -638,7 +639,7 @@ void SpinConstrain<TK>::set_decay_grad()
638639

639640
/// get decay_grad
640641
template <typename TK>
641-
const std::vector<double>& SpinConstrain<TK>::get_decay_grad()
642+
const std::vector<double>& SpinConstrain<TK>::get_decay_grad() const
642643
{
643644
return this->decay_grad_;
644645
}
@@ -684,6 +685,20 @@ double SpinConstrain<TK>::get_sc_thr() const
684685
return this->sc_thr_;
685686
}
686687

688+
/// get current adaptive sc threshold
689+
template <typename TK>
690+
double SpinConstrain<TK>::get_current_sc_thr() const
691+
{
692+
return this->current_sc_thr_;
693+
}
694+
695+
/// get computed magnetic moments Mi per atom
696+
template <typename TK>
697+
const std::vector<ModuleBase::Vector3<double>>& SpinConstrain<TK>::get_Mi() const
698+
{
699+
return this->Mi_;
700+
}
701+
687702
/// get nsc
688703
template <typename TK>
689704
int SpinConstrain<TK>::get_nsc() const

‎source/source_lcao/module_deltaspin/spin_constrain.h‎

Lines changed: 12 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -339,35 +339,10 @@ class SpinConstrain
339339
int isk);
340340
#endif
341341

342-
/// Lambda loop helper: check if RMS error below threshold or max steps reached
343-
bool check_rms_stop(int outer_step, int i_step, double rms_error, double duration, double total_duration, std::ostream& ofs_running);
344-
345-
/// Lambda loop helper: cap step size via restrict_current_ to prevent overshooting
346-
void check_restriction(const std::vector<ModuleBase::Vector3<double>>& search, double& alpha_trial, std::ostream& ofs_running);
347-
348-
/**
349-
* @brief Lambda loop helper: check if dM/dlambda gradient has decayed below threshold.
350-
*
351-
* @details Computes the diagonal of the susceptibility matrix dM/dlambda for each
352-
* atom type. If max gradient < decay_grad[itype], the lambda optimization has
353-
* reached diminishing returns and should stop.
354-
*
355-
* @return true if gradient decayed below threshold, false otherwise
356-
*/
357-
bool check_gradient_decay(std::vector<ModuleBase::Vector3<double>> new_spin,
358-
std::vector<ModuleBase::Vector3<double>> old_spin,
359-
std::vector<ModuleBase::Vector3<double>> new_delta_lambda,
360-
std::vector<ModuleBase::Vector3<double>> old_delta_lambda,
361-
bool print,
362-
std::ostream& ofs_running);
363-
/// @brief Lambda loop helper: calculate optimal step size via linear interpolation
364-
double cal_alpha_opt(std::vector<ModuleBase::Vector3<double>> spin,
365-
std::vector<ModuleBase::Vector3<double>> spin_plus,
366-
const double alpha_trial);
367-
/// Print header at start of lambda loop
368-
void print_header(std::ostream& ofs_running);
369-
/// Print termination message with final spin and lambda values
370-
void print_termination(std::ostream& ofs_running);
342+
/// Lambda loop helpers (print_rms_stop, check_restriction, check_gradient_decay,
343+
/// cal_alpha_opt, print_header, print_termination) have been lifted to free
344+
/// functions in lambda_loop_helper.h. The class now only carries state and
345+
/// the core lambda-loop driver (run_lambda_loop).
371346

372347
/// Print magnetic moments to output stream
373348
void print_Mi(std::ofstream& ofs_running);
@@ -450,9 +425,9 @@ class SpinConstrain
450425
/// get constrain
451426
const std::vector<ModuleBase::Vector3<int>>& get_constrain() const;
452427
/// get nat
453-
int get_nat();
428+
int get_nat() const;
454429
/// get ntype
455-
int get_ntype();
430+
int get_ntype() const;
456431
/// check atomCounts
457432
void check_atomCounts();
458433
/// get iat
@@ -464,11 +439,11 @@ class SpinConstrain
464439
/// zero atomic magnetic moment
465440
void zero_Mi();
466441
/// get decay_grad
467-
double get_decay_grad(int itype);
442+
double get_decay_grad(int itype) const;
468443
/// set decay_grad
469444
void set_decay_grad();
470445
/// get decay_grad
471-
const std::vector<double>& get_decay_grad();
446+
const std::vector<double>& get_decay_grad() const;
472447
/// set decay_grad from variable
473448
void set_decay_grad(const double* decay_grad_in, int ntype_in);
474449
/// set decay grad switch
@@ -482,6 +457,8 @@ class SpinConstrain
482457
double sc_drop_thr_in);
483458
/// get sc_thr
484459
double get_sc_thr() const;
460+
/// get current adaptive sc threshold (max(initial_rms * sc_drop_thr_, sc_thr_))
461+
double get_current_sc_thr() const;
485462
/// get nsc
486463
int get_nsc() const;
487464
/// get nsc_min
@@ -492,6 +469,8 @@ class SpinConstrain
492469
double get_sccut() const;
493470
/// get sc_drop_thr
494471
double get_sc_drop_thr() const;
472+
/// get computed magnetic moments Mi per atom
473+
const std::vector<ModuleBase::Vector3<double>>& get_Mi() const;
495474
/// @brief set orbital parallel info
496475
void set_ParaV(Parallel_Orbitals* ParaV_in);
497476
/// @brief set parameters for solver

0 commit comments

Comments
 (0)