Skip to content

Commit 7ed63bc

Browse files
author
dyzheng
committed
Refactor: new real space projection DeltaSpin method
1 parent edfa4c1 commit 7ed63bc

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

48 files changed

+2317
-1463
lines changed

‎source/Makefile.Objects‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -691,7 +691,6 @@ OBJS_DFTU=dftu.o\
691691
dftu_hamilt.o
692692

693693
OBJS_DELTASPIN=basic_funcs.o\
694-
cal_h_lambda.o\
695694
cal_mw_from_lambda.o\
696695
cal_mw_helper.o\
697696
cal_mw.o\

‎source/module_cell/atom_spec.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,8 @@ class Atom
3939
ModuleBase::Vector3<double> *taud;// Direct coordinates of each atom in this type.
4040
ModuleBase::Vector3<double> *vel;// velocities of each atom in this type.
4141
ModuleBase::Vector3<double> *force; // force acting on each atom in this type.
42+
ModuleBase::Vector3<double> *lambda; // Lagrange multiplier for each atom in this type. used in deltaspin
43+
ModuleBase::Vector3<int> *constrain; // constrain for each atom in this type. used in deltaspin
4244
std::string label_orb; // atomic Element symbol in the orbital file of lcao
4345

4446
double* mag;

‎source/module_cell/read_atoms.cpp‎

Lines changed: 56 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -515,16 +515,20 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
515515
delete[] atoms[it].angle1;
516516
delete[] atoms[it].angle2;
517517
delete[] atoms[it].m_loc_;
518-
atoms[it].tau = new ModuleBase::Vector3<double>[na];
518+
delete[] atoms[it].lambda;
519+
delete[] atoms[it].constrain;
520+
atoms[it].tau = new ModuleBase::Vector3<double>[na];
519521
atoms[it].dis = new ModuleBase::Vector3<double>[na];
520-
atoms[it].taud = new ModuleBase::Vector3<double>[na];
522+
atoms[it].taud = new ModuleBase::Vector3<double>[na];
521523
atoms[it].vel = new ModuleBase::Vector3<double>[na];
522-
atoms[it].mbl = new ModuleBase::Vector3<int>[na];
524+
atoms[it].mbl = new ModuleBase::Vector3<int>[na];
523525
atoms[it].mag = new double[na];
524526
atoms[it].angle1 = new double[na];
525527
atoms[it].angle2 = new double[na];
526528
atoms[it].m_loc_ = new ModuleBase::Vector3<double>[na];
527-
atoms[it].mass = this->atom_mass[it]; //mohan add 2011-11-07
529+
atoms[it].mass = this->atom_mass[it]; //mohan add 2011-11-07
530+
atoms[it].lambda = new ModuleBase::Vector3<double>[na];
531+
atoms[it].constrain = new ModuleBase::Vector3<int>[na];
528532
ModuleBase::GlobalFunc::ZEROS(atoms[it].mag,na);
529533
for (int ia = 0;ia < na; ia++)
530534
{
@@ -538,6 +542,8 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
538542
atoms[it].angle1[ia]=0;
539543
atoms[it].angle2[ia]=0;
540544
atoms[it].m_loc_[ia].set(0,0,0);
545+
atoms[it].lambda[ia].set(0,0,0);
546+
atoms[it].constrain[ia].set(0,0,0);
541547

542548
std::string tmpid;
543549
tmpid = ifpos.get();
@@ -616,7 +622,52 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
616622
atoms[it].angle2[ia]=atoms[it].angle2[ia]/180 *ModuleBase::PI;
617623
input_angle_mag=true;
618624
set_element_mag_zero = true;
619-
}
625+
}
626+
else if ( tmpid == "lambda")
627+
{
628+
double tmplam=0;
629+
ifpos >> tmplam;
630+
tmp=ifpos.get();
631+
while (tmp==' ')
632+
{
633+
tmp=ifpos.get();
634+
}
635+
if((tmp >= 48 && tmp <= 57) or tmp=='-')
636+
{
637+
ifpos.putback(tmp);
638+
ifpos >> atoms[it].lambda[ia].y>>atoms[it].lambda[ia].z;
639+
atoms[it].lambda[ia].x=tmplam;
640+
}
641+
else
642+
{
643+
ifpos.putback(tmp);
644+
atoms[it].lambda[ia].z=tmplam;
645+
}
646+
atoms[it].lambda[ia].x /= ModuleBase::Ry_to_eV;
647+
atoms[it].lambda[ia].y /= ModuleBase::Ry_to_eV;
648+
atoms[it].lambda[ia].z /= ModuleBase::Ry_to_eV;
649+
}
650+
else if ( tmpid == "sc")
651+
{
652+
double tmplam=0;
653+
ifpos >> tmplam;
654+
tmp=ifpos.get();
655+
while (tmp==' ')
656+
{
657+
tmp=ifpos.get();
658+
}
659+
if((tmp >= 48 && tmp <= 57) or tmp=='-')
660+
{
661+
ifpos.putback(tmp);
662+
ifpos >> atoms[it].constrain[ia].y>>atoms[it].constrain[ia].z;
663+
atoms[it].constrain[ia].x=tmplam;
664+
}
665+
else
666+
{
667+
ifpos.putback(tmp);
668+
atoms[it].constrain[ia].z=tmplam;
669+
}
670+
}
620671
}
621672
// move to next line
622673
while ( (tmpid != "\n") && (ifpos.good()) )

‎source/module_cell/unitcell.cpp‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -302,6 +302,48 @@ std::vector<std::vector<int>> UnitCell::get_lnchiCounts() const {
302302
return lnchiCounts;
303303
}
304304

305+
std::vector<ModuleBase::Vector3<double>> UnitCell::get_target_mag() const
306+
{
307+
std::vector<ModuleBase::Vector3<double>> target_mag(this->nat);
308+
for (int it = 0; it < this->ntype; it++)
309+
{
310+
for (int ia = 0; ia < this->atoms[it].na; ia++)
311+
{
312+
int iat = itia2iat(it, ia);
313+
target_mag[iat] = this->atoms[it].m_loc_[ia];
314+
}
315+
}
316+
return target_mag;
317+
}
318+
319+
std::vector<ModuleBase::Vector3<double>> UnitCell::get_lambda() const
320+
{
321+
std::vector<ModuleBase::Vector3<double>> lambda(this->nat);
322+
for (int it = 0; it < this->ntype; it++)
323+
{
324+
for (int ia = 0; ia < this->atoms[it].na; ia++)
325+
{
326+
int iat = itia2iat(it, ia);
327+
lambda[iat] = this->atoms[it].lambda[ia];
328+
}
329+
}
330+
return lambda;
331+
}
332+
333+
std::vector<ModuleBase::Vector3<int>> UnitCell::get_constrain() const
334+
{
335+
std::vector<ModuleBase::Vector3<int>> constrain(this->nat);
336+
for (int it = 0; it < this->ntype; it++)
337+
{
338+
for (int ia = 0; ia < this->atoms[it].na; ia++)
339+
{
340+
int iat = itia2iat(it, ia);
341+
constrain[iat] = this->atoms[it].constrain[ia];
342+
}
343+
}
344+
return constrain;
345+
}
346+
305347
void UnitCell::update_pos_tau(const double* pos) {
306348
int iat = 0;
307349
for (int it = 0; it < this->ntype; it++) {

‎source/module_cell/unitcell.h‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -328,6 +328,12 @@ class UnitCell {
328328
/// @brief get lnchiCounts, which is a vector of element type with the
329329
/// l:nchi vector
330330
std::vector<std::vector<int>> get_lnchiCounts() const;
331+
/// @brief get target magnetic moment for deltaspin
332+
std::vector<ModuleBase::Vector3<double>> get_target_mag() const;
333+
/// @brief get lagrange multiplier for deltaspin
334+
std::vector<ModuleBase::Vector3<double>> get_lambda() const;
335+
/// @brief get constrain for deltaspin
336+
std::vector<ModuleBase::Vector3<int>> get_constrain() const;
331337
};
332338

333339
/**

‎source/module_elecstate/elecstate_lcao.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -138,14 +138,14 @@ void ElecStateLCAO<TK>::init_DM(const K_Vectors* kv, const Parallel_Orbitals* pa
138138
template <>
139139
double ElecStateLCAO<double>::get_spin_constrain_energy()
140140
{
141-
SpinConstrain<double, base_device::DEVICE_CPU>& sc = SpinConstrain<double>::getScInstance();
141+
SpinConstrain<double>& sc = SpinConstrain<double>::getScInstance();
142142
return sc.cal_escon();
143143
}
144144

145145
template <>
146146
double ElecStateLCAO<std::complex<double>>::get_spin_constrain_energy()
147147
{
148-
SpinConstrain<std::complex<double>, base_device::DEVICE_CPU>& sc
148+
SpinConstrain<std::complex<double>>& sc
149149
= SpinConstrain<std::complex<double>>::getScInstance();
150150
return sc.cal_escon();
151151
}

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 12 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -673,7 +673,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
673673
// method
674674
if (PARAM.inp.sc_mag_switch && iter > PARAM.inp.sc_scf_nmin)
675675
{
676-
SpinConstrain<TK, base_device::DEVICE_CPU>& sc = SpinConstrain<TK, base_device::DEVICE_CPU>::getScInstance();
676+
SpinConstrain<TK>& sc = SpinConstrain<TK>::getScInstance();
677677
sc.run_lambda_loop(iter - 1);
678678
}
679679
}
@@ -779,30 +779,23 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2density(int istep, int iter, double ethr)
779779
}
780780
#endif
781781

782-
// 8) for delta spin
783-
if (PARAM.inp.sc_mag_switch)
784-
{
785-
SpinConstrain<TK, base_device::DEVICE_CPU>& sc = SpinConstrain<TK, base_device::DEVICE_CPU>::getScInstance();
786-
sc.cal_MW(iter, this->p_hamilt);
787-
}
788-
789-
// 9) use new charge density to calculate energy
782+
// 8) use new charge density to calculate energy
790783
this->pelec->cal_energies(1);
791784

792-
// 10) symmetrize the charge density
785+
// 9) symmetrize the charge density
793786
Symmetry_rho srho;
794787
for (int is = 0; is < PARAM.inp.nspin; is++)
795788
{
796789
srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::ucell.symm);
797790
}
798791

799-
// 11) compute magnetization, only for spin==2
792+
// 10) compute magnetization, only for nspin != 1
800793
GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx,
801794
this->pelec->charge->nxyz,
802795
this->pelec->charge->rho,
803796
this->pelec->nelec_spin.data());
804797

805-
// 12) calculate delta energy
798+
// 11) calculate delta energy for hartree and exchange-correlation double countings
806799
this->pelec->f_en.deband = this->pelec->cal_delta_eband();
807800
}
808801

@@ -911,7 +904,6 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
911904
//! 2) output charge density
912905
//! 3) output exx matrix
913906
//! 4) output charge density and density matrix
914-
//! 5) cal_MW? (why put it here?)
915907
//------------------------------------------------------------------------------
916908
template <typename TK, typename TR>
917909
void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
@@ -1054,15 +1046,6 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
10541046
}
10551047
}
10561048

1057-
// 5) cal_MW?
1058-
// escon: energy of spin constraint depends on Mi, so cal_energies should be
1059-
// called after cal_MW
1060-
if (PARAM.inp.sc_mag_switch)
1061-
{
1062-
SpinConstrain<TK, base_device::DEVICE_CPU>& sc = SpinConstrain<TK, base_device::DEVICE_CPU>::getScInstance();
1063-
sc.cal_MW(iter, this->p_hamilt);
1064-
}
1065-
10661049
// 6) use the converged occupation matrix for next MD/Relax SCF calculation
10671050
if (PARAM.inp.dft_plus_u && this->conv_esolver)
10681051
{
@@ -1223,13 +1206,14 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
12231206
}
12241207
}
12251208

1226-
// 15) write spin constrian MW?
1227-
// spin constrain calculations, added by Tianqi Zhao.
1228-
if (PARAM.inp.sc_mag_switch)
1209+
// 15) write atomic magnetization only when spin_constraint is on
1210+
// spin constrain calculations.
1211+
if (PARAM.inp.sc_mag_switch)
12291212
{
1230-
SpinConstrain<TK, base_device::DEVICE_CPU>& sc = SpinConstrain<TK, base_device::DEVICE_CPU>::getScInstance();
1231-
sc.cal_MW(istep, true);
1232-
sc.print_Mag_Force();
1213+
SpinConstrain<TK>& sc = SpinConstrain<TK>::getScInstance();
1214+
sc.cal_mi_lcao(istep);
1215+
sc.print_Mi(GlobalV::ofs_running);
1216+
sc.print_Mag_Force(GlobalV::ofs_running);
12331217
}
12341218

12351219
// 16) delete grid

‎source/module_esolver/lcao_before_scf.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -143,7 +143,7 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
143143
#endif
144144
if (PARAM.inp.sc_mag_switch)
145145
{
146-
SpinConstrain<TK, base_device::DEVICE_CPU>& sc = SpinConstrain<TK, base_device::DEVICE_CPU>::getScInstance();
146+
SpinConstrain<TK>& sc = SpinConstrain<TK>::getScInstance();
147147
sc.init_sc(PARAM.inp.sc_thr,
148148
PARAM.inp.nsc,
149149
PARAM.inp.nsc_min,

‎source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ if(ENABLE_LCAO)
1313
operator_lcao/nonlocal_new.cpp
1414
operator_lcao/td_ekinetic_lcao.cpp
1515
operator_lcao/td_nonlocal_lcao.cpp
16-
operator_lcao/sc_lambda_lcao.cpp
16+
operator_lcao/dspin_lcao.cpp
1717
operator_lcao/dftu_lcao.cpp
1818
pulay_force_stress_center2.cpp
1919
FORCE_STRESS.cpp

‎source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp‎

Lines changed: 13 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,8 @@
3333
#include "operator_lcao/op_dftu_lcao.h"
3434
#include "operator_lcao/op_exx_lcao.h"
3535
#include "operator_lcao/overlap_new.h"
36-
#include "operator_lcao/sc_lambda_lcao.h"
36+
#include "operator_lcao/dspin_lcao.h"
37+
#include "module_hamilt_lcao/module_deltaspin/spin_constrain.h"
3738
#include "operator_lcao/td_ekinetic_lcao.h"
3839
#include "operator_lcao/td_nonlocal_lcao.h"
3940
#include "operator_lcao/veff_lcao.h"
@@ -377,11 +378,18 @@ HamiltLCAO<TK, TR>::HamiltLCAO(Gint_Gamma* GG_in,
377378
}
378379
if (PARAM.inp.sc_mag_switch)
379380
{
380-
Operator<TK>* sc_lambda = new OperatorScLambda<OperatorLCAO<TK, TR>>(this->hsk,
381-
kv->kvec_d,
382-
this->hR, // no explicit call yet
383-
this->kv->isk);
381+
Operator<TK>* sc_lambda =
382+
new DeltaSpin<OperatorLCAO<TK, TR>>(
383+
this->hsk,
384+
kv->kvec_d,
385+
this->hR,
386+
GlobalC::ucell,
387+
&GlobalC::GridD,
388+
two_center_bundle.overlap_orb_onsite.get()
389+
);
384390
this->getOperator()->add(sc_lambda);
391+
SpinConstrain<TK>& sc = SpinConstrain<TK>::getScInstance();
392+
sc.set_operator(sc_lambda);
385393
}
386394
}
387395

0 commit comments

Comments
 (0)