abacus-develop/source/source_esolver/esolver_ks.cpp at develop · deepmodeling/abacus-develop · GitHub

History

335 lines (267 loc) · 11.6 KB

Raw

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

#include "esolver_ks.h"

#include "source_base/timer_wrapper.h"

// for jason output information

#include "source_io/module_json/output_info.h"

#include "source_estate/update_pot.h" // mohan add 20251016

#include "source_estate/module_charge/chgmixing.h" // mohan add 20251018

#include "source_pw/module_pwdft/setup_pwwfc.h" // mohan add 20251018

#include "source_hsolver/hsolver.h"

#include "source_io/module_energy/write_eig_occ.h"

#include "source_io/module_energy/write_bands.h"

#include "source_hamilt/module_xc/xc_functional.h"

#include "source_io/module_output/output_log.h" // use write_head

#include "source_estate/elecstate_print.h" // print_etot

#include "source_pw/module_pwdft/dftu_base.h" // Plus_U_Base::u_converged in iter_finish

#include "source_hamilt/module_xc/general_exx_info.h" // for init_general_exx_info

namespace ModuleESolver

{

ESolver_KS::ESolver_KS() {}

ESolver_KS::~ESolver_KS()

{

//****************************************************

// do not add any codes in this deconstructor funcion

//****************************************************

delete this->p_hamilt;

delete this->p_chgmix;

this->ppcell.release_memory();

// mohan add 2025-10-18, should be put int clean() function

pw::teardown_pwwfc(this->pw_wfc);

}

void ESolver_KS::before_all_runners(BaseCell& basecell, const Input_para& inp)

{

basecell.require_kind(BaseCell::Kind::unitcell, __FUNCTION__);

UnitCell& ucell = static_cast<UnitCell&>(basecell);

ModuleBase::TITLE("ESolver_KS", "before_all_runners");

//! 1) setup "before_all_runniers" in ESolver_FP

ESolver_FP::before_all_runners(ucell, inp);

//! 2) initialize General_Exx_Info from input parameters

init_general_exx_info(general_exx_info_, inp);

//! 3) setup some parameters

classname = "ESolver_KS";

basisname = "";

this->scf_thr = inp.scf_thr;

this->scf_ene_thr = inp.scf_ene_thr;

this->maxniter = inp.scf_nmax;

this->niter = maxniter;

this->drho = 0.0;

// cell_factor

this->ppcell.cell_factor = inp.cell_factor;

//! 3) setup charge mixing

p_chgmix = new Charge_Mixing();

p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod);

p_chgmix->set_mixing(inp.mixing_mode, inp.mixing_beta, inp.mixing_ndim,

inp.mixing_gg0, inp.mixing_tau, inp.mixing_beta_mag, inp.mixing_gg0_mag,

inp.mixing_gg0_min, inp.mixing_angle, inp.mixing_dmr, ucell.omega, ucell.tpiba);

p_chgmix->init_mixing();

//! 4) setup plane wave for electronic wave functions

pw::setup_pwwfc(inp, ucell, *this->pw_rho, this->kv, this->pw_wfc);

//! 5) read in charge density, mohan add 2025-11-28

//! Inititlize the charge density.

this->chr.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv, this->pw_wfc);

this->chr.check_rho(); // check the rho

}

void ESolver_KS::hamilt2rho_single(UnitCell& ucell, const int istep, const int iter, const double ethr)

{}

void ESolver_KS::hamilt2rho(UnitCell& ucell, const int istep, const int iter, const double ethr)

{

// 1) use Hamiltonian to obtain charge density

this->hamilt2rho_single(ucell, istep, iter, diag_ethr);

// 2) for MPI: STOGROUP? need to rewrite

//<Temporary> It may be changed when more clever parallel algorithm is put forward.

// When parallel algorithm for bands are adopted. Density will only be treated in the first group.

//(Different ranks should have abtained the same, but small differences always exist in practice.)

// Maybe in the future, density and wavefunctions should use different

// parallel algorithms, in which they do not occupy all processors, for

// example wavefunctions uses 20 processors while density uses 10.

if (PARAM.globalv.ks_run)

{

drho = p_chgmix->get_drho(&this->chr, this->inp_->nelec);

hsolver_error = 0.0;

if (iter == 1 && this->inp_->calculation != "nscf")

{

hsolver_error

= hsolver::cal_hsolve_error(this->inp_->basis_type, this->inp_->esolver_type, diag_ethr, this->inp_->nelec);

// The error of HSolver is larger than drho,

// so a more precise HSolver should be executed.

if (hsolver_error > drho)

{

diag_ethr = hsolver::reset_diag_ethr(GlobalV::ofs_running, this->inp_->basis_type,

this->inp_->esolver_type, this->inp_->precision, hsolver_error,

drho, diag_ethr, this->inp_->nelec);

this->hamilt2rho_single(ucell, istep, iter, diag_ethr);

drho = p_chgmix->get_drho(&this->chr, this->inp_->nelec);

hsolver_error = hsolver::cal_hsolve_error(this->inp_->basis_type,

this->inp_->esolver_type, diag_ethr, this->inp_->nelec);

}

void ESolver_KS::runner(BaseCell& basecell, const int istep)

{

basecell.require_kind(BaseCell::Kind::unitcell, __FUNCTION__);

UnitCell& ucell = static_cast<UnitCell&>(basecell);

ModuleBase::TITLE("ESolver_KS", "runner");

ModuleBase::timer::start(this->classname, "runner");

// 1) before_scf (electronic iteration loops)

this->before_scf(ucell, istep);

ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");

// 2) SCF iterations

bool conv_esolver = false;

this->niter = this->maxniter;

this->diag_ethr = this->inp_->pw_diag_thr;

this->scf_nmax_flag = false; // mohan add 2025-09-21

for (int iter = 1; iter <= this->maxniter; ++iter)

{

if(iter == this->maxniter)

{

this->scf_nmax_flag=true;

}

// 3) initialization of SCF iterations

this->iter_init(ucell, istep, iter);

// 4) use Hamiltonian to obtain charge density

this->hamilt2rho(ucell, istep, iter, diag_ethr);

// 5) finish scf iterations

this->iter_finish(ucell, istep, iter, conv_esolver);

// 6) check convergence

if (conv_esolver || this->oscillate_esolver)

{

this->niter = iter;

if (this->oscillate_esolver)

{

std::cout << " !! Density oscillation is found, STOP HERE !!" << std::endl;

}

break;

}

} // end scf iterations

// 7) after scf

this->after_scf(ucell, istep, conv_esolver);

this->conv_esolver = conv_esolver;

ModuleBase::timer::end(this->classname, "runner");

return;

};

void ESolver_KS::before_scf(UnitCell& ucell, const int istep)

{

ModuleBase::TITLE("ESolver_KS", "before_scf");

ESolver_FP::before_scf(ucell, istep);

}

void ESolver_KS::iter_init(UnitCell& ucell, const int istep, const int iter)

{

if(this->inp_->esolver_type != "tddft")

{

ModuleIO::write_head(GlobalV::ofs_running, istep, iter, this->basisname);

}

iter_time = ModuleBase::get_time();

// Reset DeltaSpin RMS for this SCF iteration. If DeltaSpin is enabled, use 0

// (meaning "lambda loop not yet run this iteration"); otherwise -1 (no RMS column).

this->ds_rms_ = this->inp_->sc_mag_switch ? 0.0 : -1.0;

if (this->inp_->esolver_type == "ksdft")

{

diag_ethr = hsolver::set_diagethr_ks(this->inp_->basis_type, this->inp_->esolver_type,

this->inp_->calculation, this->inp_->init_chg, this->inp_->precision, istep, iter,

drho, this->inp_->pw_diag_thr, diag_ethr, this->inp_->nelec, this->inp_->scf_thr);

}

else if (this->inp_->esolver_type == "sdft")

{

diag_ethr = hsolver::set_diagethr_sdft(this->inp_->basis_type, this->inp_->esolver_type,

this->inp_->calculation, this->inp_->init_chg, istep, iter, drho,

this->inp_->pw_diag_thr, diag_ethr, this->inp_->nbands, esolver_KS_ne,

this->inp_->nelec, this->inp_->scf_thr);

}

// save input charge density (rho)

this->chr.save_rho_before_sum_band();

}

void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver)

{

// 1.1) print out band gap

if (!PARAM.globalv.two_fermi)

{

this->pelec->cal_bandgap();

}

else

{

this->pelec->cal_bandgap_updw();

}

// 1.2) print out eigenvalues and occupations

if (this->inp_->out_band[0])

{

if (iter % this->inp_->out_freq_elec == 0 || iter == this->inp_->scf_nmax || conv_esolver)

{

ModuleIO::write_eig_iter(this->pelec->ekb,this->pelec->wg,*this->pelec->klist);

}

// 2.1) compute magnetization, only for spin==2

ucell.magnet.compute_mag(ucell.omega, this->chr.nrxx, this->chr.nxyz, this->chr.rho,

this->inp_->nspin, PARAM.globalv.two_fermi, this->inp_->nelec,

this->pelec->nelec_spin.data());

// 2.2) charge mixing

// SCF will continue if U is not converged for uramping calculation

bool converged_u = true;

// to avoid unnecessary dependence on dft+u, refactor is needed

#ifdef __LCAO

if (this->inp_->dft_plus_u)

{

converged_u = this->dftu_->u_converged();

}

#endif

module_charge::chgmixing_ks(iter, ucell, this->pelec, this->chr, this->p_chgmix,

this->pw_rhod->nrxx, this->drho, this->oscillate_esolver, conv_esolver, hsolver_error,

this->scf_thr, this->scf_ene_thr, converged_u, *this->inp_);

// 2.3) Update potentials (should be done every SF iter)

elecstate::update_pot(ucell, this->pelec, this->chr, conv_esolver);

// 3.1) calculate energies

this->pelec->cal_energies(1); // Harris-Foulkes functional

this->pelec->cal_energies(2); // Kohn-Sham functional

if (iter == 1)

{

this->pelec->f_en.etot_old = this->pelec->f_en.etot;

}

this->pelec->f_en.etot_delta = this->pelec->f_en.etot - this->pelec->f_en.etot_old;

this->pelec->f_en.etot_old = this->pelec->f_en.etot;

// 4) get meta-GGA related parameters

double dkin = 0.0; // for meta-GGA

if (XC_Functional::get_ked_flag())

{

dkin = p_chgmix->get_dkin(&this->chr, this->inp_->nelec);

}

// Iter finish

ESolver_FP::iter_finish(ucell, istep, iter, conv_esolver);

// the end, print time

double duration = ModuleBase::get_duration(iter_time, ModuleBase::get_time());

// print energies

elecstate::print_etot(ucell.magnet, *pelec, conv_esolver, iter, drho,

dkin, duration, diag_ethr, 0, true, this->ds_rms_);

#ifdef __RAPIDJSON

// add Json of scf mag

Json::add_output_scf_mag(ucell.magnet.tot_mag, ucell.magnet.abs_mag,

this->pelec->f_en.etot * ModuleBase::Ry_to_eV,

this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV,

drho, duration);

#endif //__RAPIDJSON

}

//! Something to do after SCF iterations when SCF is converged or comes to the max iter step.

void ESolver_KS::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver)

{

ModuleBase::TITLE("ESolver_KS", "after_scf");

// 1) calculate the kinetic energy density tau

if (this->inp_->out_elf[0] > 0)

{

assert(this->psi != nullptr);

this->pelec->cal_tau(*(this->psi));

}

// 2) call after_scf() of ESolver_FP

ESolver_FP::after_scf(ucell, istep, conv_esolver);

// 3) write eigenvalues and occupations to eig_occ.txt

ModuleIO::write_eig_file(this->pelec->ekb, this->pelec->wg, this->kv, istep);

// 4) write band information to band.txt

ModuleIO::write_bands(*this->inp_, this->pelec->ekb, this->kv);

}

void ESolver_KS::after_all_runners(BaseCell& basecell)

{

basecell.require_kind(BaseCell::Kind::unitcell, __FUNCTION__);

UnitCell& ucell = static_cast<UnitCell&>(basecell);

// 1) write Etot information

ESolver_FP::after_all_runners(ucell);

}

} // namespace ModuleESolver

abacus-develop/source/source_esolver/esolver_ks.cpp at develop · deepmodeling/abacus-develop · GitHub

Uh oh!

Expand file tree

Breadcrumbs

esolver_ks.cpp

Latest commit

History

Breadcrumbs

esolver_ks.cpp

File metadata and controls