|
24 | 24 | #include "source_io/ctrl_runner_lcao.h" // use ctrl_runner_lcao() |
25 | 25 | #include "source_io/ctrl_iter_lcao.h" // use ctrl_iter_lcao() |
26 | 26 |
|
| 27 | +// tmp |
| 28 | +#include "source_estate/module_dm/cal_dm_psi.h" |
| 29 | + |
27 | 30 | namespace ModuleESolver |
28 | 31 | { |
29 | 32 |
|
@@ -389,12 +392,49 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const |
389 | 392 | = GlobalC::exx_info.info_ri.real_number ? this->exd->two_level_step : this->exc->two_level_step; |
390 | 393 | } |
391 | 394 | #endif |
392 | | - elecstate::setup_dm<TK>(ucell, estate, this->psi, this->chr, iter, exx_two_level_step); |
| 395 | +// elecstate::setup_dm<TK>(ucell, estate, this->psi, this->chr, iter, exx_two_level_step); |
| 396 | + |
| 397 | + |
| 398 | + if (iter == 1 && exx_two_level_step == 0) |
| 399 | + { |
| 400 | + std::cout << " WAVEFUN -> CHARGE " << std::endl; |
| 401 | + |
| 402 | + // calculate the density matrix using read in wave functions |
| 403 | + // and then calculate the charge density on grid. |
| 404 | + |
| 405 | + estate->skip_weights = true; |
| 406 | + elecstate::calculate_weights(estate->ekb, |
| 407 | + estate->wg, |
| 408 | + estate->klist, |
| 409 | + estate->eferm, |
| 410 | + estate->f_en, |
| 411 | + estate->nelec_spin, |
| 412 | + estate->skip_weights); |
| 413 | + |
| 414 | + elecstate::calEBand(estate->ekb, estate->wg, estate->f_en); |
| 415 | + elecstate::cal_dm_psi(estate->DM->get_paraV_pointer(), estate->wg, *this->psi, *(estate->DM)); |
| 416 | + estate->DM->cal_DMR(); |
| 417 | + |
| 418 | + estate->psiToRho(*this->psi); |
| 419 | + estate->skip_weights = false; |
| 420 | + |
| 421 | + elecstate::cal_ux(ucell); |
| 422 | + |
| 423 | + //! update the potentials by using new electron charge density |
| 424 | + estate->pot->update_from_charge(&this->chr, &ucell); |
| 425 | + |
| 426 | + //! compute the correction energy for metals |
| 427 | + estate->f_en.descf = estate->cal_delta_escf(); |
| 428 | + } |
| 429 | + |
| 430 | + |
| 431 | + |
393 | 432 | } |
394 | 433 |
|
395 | 434 | #ifdef __EXX |
396 | 435 | // calculate exact-exchange |
397 | 436 | if (PARAM.inp.calculation != "nscf") |
| 437 | +q |
398 | 438 | { |
399 | 439 | if (GlobalC::exx_info.info_ri.real_number) |
400 | 440 | { |
|
0 commit comments