1!! Copyright (C) 2002-2014 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira 2!! 3!! This program is free software; you can redistribute it and/or modify 4!! it under the terms of the GNU General Public License as published by 5!! the Free Software Foundation; either version 2, or (at your option) 6!! any later version. 7!! 8!! This program is distributed in the hope that it will be useful, 9!! but WITHOUT ANY WARRANTY; without even the implied warranty of 10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11!! GNU General Public License for more details. 12!! 13!! You should have received a copy of the GNU General Public License 14!! along with this program; if not, write to the Free Software 15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 16!! 02110-1301, USA. 17!! 18 19#include "global.h" 20 21module scf_oct_m 22 use batch_oct_m 23 use batch_ops_oct_m 24 use berry_oct_m 25 use density_oct_m 26 use eigensolver_oct_m 27 use energy_calc_oct_m 28 use forces_oct_m 29 use geometry_oct_m 30 use global_oct_m 31 use grid_oct_m 32 use hamiltonian_elec_oct_m 33 use io_oct_m 34 use kpoints_oct_m 35 use lcao_oct_m 36 use lda_u_oct_m 37 use lda_u_io_oct_m 38 use lda_u_mixer_oct_m 39 use loct_oct_m 40 use magnetic_oct_m 41 use mesh_oct_m 42 use mesh_function_oct_m 43 use messages_oct_m 44 use mix_oct_m 45 use modelmb_exchange_syms_oct_m 46 use mpi_oct_m 47 use multicomm_oct_m 48 use namespace_oct_m 49 use output_oct_m 50 use parser_oct_m 51 use partial_charges_oct_m 52 use preconditioners_oct_m 53 use profiling_oct_m 54 use restart_oct_m 55 use scdm_oct_m 56 use simul_box_oct_m 57 use smear_oct_m 58 use species_oct_m 59 use states_abst_oct_m 60 use states_elec_oct_m 61 use states_elec_dim_oct_m 62 use states_elec_group_oct_m 63 use states_elec_io_oct_m 64 use states_elec_restart_oct_m 65 use stress_oct_m 66 use symmetries_oct_m 67 use types_oct_m 68 use unit_oct_m 69 use unit_system_oct_m 70 use utils_oct_m 71 use v_ks_oct_m 72 use varinfo_oct_m 73 use vdw_ts_oct_m 74! use xc_functl_oct_m 75 use walltimer_oct_m 76 use wfs_elec_oct_m 77 use XC_F90(lib_m) 78 use xc_oep_oct_m 79 80 implicit none 81 82 private 83 public :: & 84 scf_t, & 85 scf_init, & 86 scf_mix_clear, & 87 scf_run, & 88 scf_end, & 89 scf_state_info, & 90 scf_print_mem_use 91 92 integer, public, parameter :: & 93 VERB_NO = 0, & 94 VERB_COMPACT = 1, & 95 VERB_FULL = 3 96 97 !> some variables used for the SCF cycle 98 type scf_t 99 private 100 integer, public :: max_iter !< maximum number of SCF iterations 101 integer :: max_iter_berry !< max number of electronic iterations before updating density, for Berry potential 102 103 FLOAT, public :: lmm_r 104 105 ! several convergence criteria 106 FLOAT :: conv_abs_dens, conv_rel_dens, conv_abs_ev, conv_rel_ev, conv_abs_force 107 FLOAT :: abs_dens, rel_dens, abs_ev, rel_ev, abs_force 108 FLOAT :: conv_energy_diff 109 FLOAT :: energy_diff 110 logical :: conv_eigen_error 111 logical :: check_conv 112 113 integer :: mix_field 114 logical :: lcao_restricted 115 logical :: calc_force 116 logical :: calc_stress 117 logical :: calc_dipole 118 logical :: calc_partial_charges 119 type(mix_t) :: smix 120 type(mixfield_t), pointer :: mixfield 121 type(eigensolver_t) :: eigens 122 integer :: mixdim1 123 logical :: forced_finish !< remember if 'touch stop' was triggered earlier. 124 type(lda_u_mixer_t) :: lda_u_mix 125 type(grid_t), pointer :: gr 126 end type scf_t 127 128contains 129 130 ! --------------------------------------------------------- 131 subroutine scf_init(scf, namespace, gr, geo, st, mc, hm, ks, conv_force) 132 type(scf_t), intent(inout) :: scf 133 type(namespace_t), intent(in) :: namespace 134 type(grid_t), target, intent(inout) :: gr 135 type(geometry_t), intent(in) :: geo 136 type(states_elec_t), intent(in) :: st 137 type(multicomm_t), intent(in) :: mc 138 type(hamiltonian_elec_t), intent(inout) :: hm 139 type(v_ks_t), intent(in) :: ks 140 FLOAT, optional, intent(in) :: conv_force 141 142 FLOAT :: rmin 143 integer :: mixdefault, ierr 144 type(type_t) :: mix_type 145 146 PUSH_SUB(scf_init) 147 148 !%Variable MaximumIter 149 !%Type integer 150 !%Default 200 151 !%Section SCF::Convergence 152 !%Description 153 !% Maximum number of SCF iterations. The code will stop even if convergence 154 !% has not been achieved. -1 means unlimited. 155 !% 0 means just do LCAO (or read from restart), compute the eigenvalues and energy, 156 !% and stop, without updating the wavefunctions or density. 157 !% 158 !% If convergence criteria are set, the SCF loop will only stop once the criteria 159 !% are fulfilled for two consecutive iterations. 160 !%End 161 call parse_variable(namespace, 'MaximumIter', 200, scf%max_iter) 162 163 !%Variable MaximumIterBerry 164 !%Type integer 165 !%Default 10 166 !%Section SCF::Convergence 167 !%Description 168 !% Maximum number of iterations for the Berry potential, within each SCF iteration. 169 !% Only applies if a <tt>StaticElectricField</tt> is applied in a periodic direction. 170 !% The code will move on to the next SCF iteration even if convergence 171 !% has not been achieved. -1 means unlimited. 172 !%End 173 if(associated(hm%vberry)) then 174 call parse_variable(namespace, 'MaximumIterBerry', 10, scf%max_iter_berry) 175 if(scf%max_iter_berry < 0) scf%max_iter_berry = huge(scf%max_iter_berry) 176 end if 177 178 !%Variable ConvEnergy 179 !%Type float 180 !%Default 0.0 181 !%Section SCF::Convergence 182 !%Description 183 !% Stop the SCF when the magnitude of change in energy during at 184 !% one SCF iteration is smaller than this value. 185 !% 186 !%A zero value (the default) means do not use this criterion. 187 !% 188 !% If this criterion is used, the SCF loop will only stop once it is 189 !% fulfilled for two consecutive iterations. 190 !%End 191 call parse_variable(namespace, 'ConvEnergy', M_ZERO, scf%conv_energy_diff, unit = units_inp%energy) 192 193 !%Variable ConvAbsDens 194 !%Type float 195 !%Default 0.0 196 !%Section SCF::Convergence 197 !%Description 198 !% Absolute convergence of the density: 199 !% 200 !% <math>\varepsilon = \int {\rm d}^3r \left| \rho^{out}(\bf r) -\rho^{inp}(\bf r) \right|</math>. 201 !% 202 !% A zero value (the default) means do not use this criterion. 203 !% 204 !% If this criterion is used, the SCF loop will only stop once it is 205 !% fulfilled for two consecutive iterations. 206 !%End 207 call parse_variable(namespace, 'ConvAbsDens', M_ZERO, scf%conv_abs_dens) 208 209 !%Variable ConvRelDens 210 !%Type float 211 !%Default 1e-6 212 !%Section SCF::Convergence 213 !%Description 214 !% Relative convergence of the density: 215 !% 216 !% <math>\varepsilon = \frac{1}{N} \mathrm{ConvAbsDens}</math>. 217 !% 218 !% <i>N</i> is the total number of electrons in the problem. A 219 !% zero value means do not use this criterion. 220 !% 221 !% If you reduce this value, you should also reduce 222 !% <tt>EigensolverTolerance</tt> to a value of roughly 1/10 of 223 !% <tt>ConvRelDens</tt> to avoid convergence problems. 224 !% 225 !% If this criterion is used, the SCF loop will only stop once it is 226 !% fulfilled for two consecutive iterations. 227 !%End 228 call parse_variable(namespace, 'ConvRelDens', CNST(1e-6), scf%conv_rel_dens) 229 230 !%Variable ConvAbsEv 231 !%Type float 232 !%Default 0.0 233 !%Section SCF::Convergence 234 !%Description 235 !% Absolute convergence of the sum of the eigenvalues: 236 !% 237 !% <math> \varepsilon = \left| \sum_{j=1}^{N_{occ}} \varepsilon_j^{out} - 238 !% \sum_{j=1}^{N_{occ}} \varepsilon_j^{inp} \right| </math> 239 !% 240 !% A zero value (the default) means do not use this criterion. 241 !% 242 !% If this criterion is used, the SCF loop will only stop once it is 243 !% fulfilled for two consecutive iterations. 244 !%End 245 call parse_variable(namespace, 'ConvAbsEv', M_ZERO, scf%conv_abs_ev, unit = units_inp%energy) 246 247 !%Variable ConvRelEv 248 !%Type float 249 !%Default 0.0 250 !%Section SCF::Convergence 251 !%Description 252 !% Relative convergence of the sum of the eigenvalues: 253 !% 254 !% <math>\varepsilon = \frac{ \left| \sum_{j=1}^{N_{occ}} ( \varepsilon_j^{out} - \varepsilon_j^{inp} ) \right|} 255 !% {\left| \sum_{j=1}^{N_{occ}} \varepsilon_j^{out} \right|} </math> 256 !% 257 !%A zero value (the default) means do not use this criterion. 258 !% 259 !% If this criterion is used, the SCF loop will only stop once it is 260 !% fulfilled for two consecutive iterations. 261 !%End 262 call parse_variable(namespace, 'ConvRelEv', M_ZERO, scf%conv_rel_ev) 263 264 call messages_obsolete_variable(namespace, 'ConvAbsForce', 'ConvForce') 265 call messages_obsolete_variable(namespace, 'ConvRelForce', 'ConvForce') 266 267 !%Variable ConvForce 268 !%Type float 269 !%Section SCF::Convergence 270 !%Description 271 !% Absolute convergence of the forces: maximum variation of any 272 !% component of the ionic forces in consecutive iterations. A 273 !% zero value means do not use this criterion. The default is 274 !% zero, except for geometry optimization, which sets a default of 275 !% 1e-8 H/b. 276 !% 277 !% If this criterion is used, the SCF loop will only stop once it is 278 !% fulfilled for two consecutive iterations. 279 !%End 280 call parse_variable(namespace, 'ConvForce', optional_default(conv_force, M_ZERO), scf%conv_abs_force, unit = units_inp%force) 281 282 scf%check_conv = & 283 scf%conv_energy_diff > M_ZERO .or. & 284 scf%conv_abs_dens > M_ZERO .or. scf%conv_rel_dens > M_ZERO .or. & 285 scf%conv_abs_ev > M_ZERO .or. scf%conv_rel_ev > M_ZERO .or. & 286 scf%conv_abs_force > M_ZERO 287 288 if(.not. scf%check_conv .and. scf%max_iter < 0) then 289 call messages_write("All convergence criteria are disabled. Octopus is cowardly refusing") 290 call messages_new_line() 291 call messages_write("to enter an infinite loop.") 292 call messages_new_line() 293 call messages_new_line() 294 call messages_write("Please set one of the following variables to a positive value:") 295 call messages_new_line() 296 call messages_new_line() 297 call messages_write(" | MaximumIter | ConvEnergy | ConvAbsDens | ConvRelDens |") 298 call messages_new_line() 299 call messages_write(" | ConvAbsEv | ConvRelEv | ConvForce |") 300 call messages_new_line() 301 call messages_fatal() 302 end if 303 304 !%Variable ConvEigenError 305 !%Type logical 306 !%Default false 307 !%Section SCF::Convergence 308 !%Description 309 !% If true, the calculation will not be considered converged unless all states have 310 !% individual errors less than <tt>EigensolverTolerance</tt>. 311 !% 312 !% If this criterion is used, the SCF loop will only stop once it is 313 !% fulfilled for two consecutive iterations. 314 !%End 315 call parse_variable(namespace, 'ConvEigenError', .false., scf%conv_eigen_error) 316 317 if(scf%max_iter < 0) scf%max_iter = huge(scf%max_iter) 318 319 call messages_obsolete_variable(namespace, 'What2Mix', 'MixField') 320 321 !%Variable MixField 322 !%Type integer 323 !%Section SCF::Mixing 324 !%Description 325 !% Selects what should be mixed during the SCF cycle. Note that 326 !% currently the exact-exchange part of hybrid functionals is not 327 !% mixed at all, which would require wavefunction-mixing, not yet 328 !% implemented. This may lead to instabilities in the SCF cycle, 329 !% so starting from a converged LDA/GGA calculation is recommended 330 !% for hybrid functionals. The default depends on the <tt>TheoryLevel</tt> 331 !% and the exchange-correlation potential used. 332 !%Option none 0 333 !% No mixing is done. This is the default for independent 334 !% particles. 335 !%Option potential 1 336 !% The Kohn-Sham potential is mixed. This is the default for other cases. 337 !%Option density 2 338 !% Mix the density. 339 !%Option states 3 340 !% (Experimental) Mix the states. In this case, the mixing is always linear. 341 !%End 342 343 mixdefault = OPTION__MIXFIELD__POTENTIAL 344 if(hm%theory_level == INDEPENDENT_PARTICLES) mixdefault = OPTION__MIXFIELD__NONE 345 346 call parse_variable(namespace, 'MixField', mixdefault, scf%mix_field) 347 if(.not.varinfo_valid_option('MixField', scf%mix_field)) call messages_input_error(namespace, 'MixField') 348 call messages_print_var_option(stdout, 'MixField', scf%mix_field, "what to mix during SCF cycles") 349 350 if (scf%mix_field == OPTION__MIXFIELD__POTENTIAL .and. hm%theory_level == INDEPENDENT_PARTICLES) then 351 call messages_write('Input: Cannot mix the potential for non-interacting particles.') 352 call messages_fatal() 353 end if 354 355 if (scf%mix_field == OPTION__MIXFIELD__POTENTIAL .and. hm%pcm%run_pcm) then 356 call messages_write('Input: You have selected to mix the potential.', new_line = .true.) 357 call messages_write(' This might produce convergence problems for solvated systems.', new_line = .true.) 358 call messages_write(' Mix the Density instead.') 359 call messages_warning() 360 end if 361 362 if(scf%mix_field == OPTION__MIXFIELD__DENSITY & 363 .and. bitand(hm%xc%family, XC_FAMILY_OEP + XC_FAMILY_MGGA + XC_FAMILY_HYB_MGGA) /= 0) then 364 365 call messages_write('Input: You have selected to mix the density with OEP or MGGA XC functionals.', new_line = .true.) 366 call messages_write(' This might produce convergence problems. Mix the potential instead.') 367 call messages_warning() 368 end if 369 370 if(scf%mix_field == OPTION__MIXFIELD__STATES) then 371 call messages_experimental('MixField = states') 372 end if 373 374 ! Handle mixing now... 375 select case(scf%mix_field) 376 case(OPTION__MIXFIELD__POTENTIAL) 377 scf%mixdim1 = gr%mesh%np 378 case(OPTION__MIXFIELD__DENSITY) 379 scf%mixdim1 = gr%fine%mesh%np 380 case(OPTION__MIXFIELD__STATES) 381 ! we do not really need the mixer, except for the value of the mixing coefficient 382 scf%mixdim1 = 1 383 end select 384 385 mix_type = TYPE_FLOAT 386 387 if(scf%mix_field == OPTION__MIXFIELD__DENSITY) then 388 call mix_init(scf%smix, namespace, gr%fine%der, scf%mixdim1, 1, st%d%nspin, func_type_ = mix_type) 389 else if(scf%mix_field /= OPTION__MIXFIELD__NONE) then 390 call mix_init(scf%smix, namespace, gr%der, scf%mixdim1, 1, st%d%nspin, func_type_ = mix_type) 391 end if 392 393 !If we use LDA+U, we also have do mix it 394 if(scf%mix_field /= OPTION__MIXFIELD__STATES) then 395 call lda_u_mixer_init(hm%lda_u, scf%lda_u_mix, st) 396 call lda_u_mixer_init_auxmixer(hm%lda_u, scf%lda_u_mix, scf%smix, st) 397 end if 398 call mix_get_field(scf%smix, scf%mixfield) 399 400 if(hm%lda_u_level /= DFT_U_NONE .and. hm%lda_u%basisfromstates) then 401 call lda_u_loadbasis(hm%lda_u, namespace, st, gr%mesh, mc, ierr) 402 if (ierr /= 0) then 403 message(1) = "Unable to load LDA+U basis from selected states." 404 call messages_fatal(1) 405 end if 406 call lda_u_periodic_coulomb_integrals(hm%lda_u, namespace, st, gr%der, mc, associated(hm%hm_base%phase)) 407 end if 408 409 ! now the eigensolver stuff 410 call eigensolver_init(scf%eigens, namespace, gr, st, geo, mc) 411 412 !The evolution operator is a very specific propagation that requires a specific 413 !setting to work in the current framework 414 if(scf%eigens%es_type == RS_EVO) then 415 if(scf%mix_field /= OPTION__MIXFIELD__DENSITY) then 416 message(1) = "Evolution eigensolver is only compatible with MixField = density." 417 call messages_fatal(1) 418 end if 419 if(mix_coefficient(scf%smix) /= M_ONE) then 420 message(1) = "Evolution eigensolver is only compatible with Mixing = 1." 421 call messages_fatal(1) 422 end if 423 if(mix_scheme(scf%smix) /= OPTION__MIXINGSCHEME__LINEAR) then 424 message(1) = "Evolution eigensolver is only compatible with MixingScheme = linear." 425 call messages_fatal(1) 426 end if 427 end if 428 429 !%Variable SCFinLCAO 430 !%Type logical 431 !%Default no 432 !%Section SCF 433 !%Description 434 !% Performs the SCF cycle with the calculation restricted to the LCAO subspace. 435 !% This may be useful for systems with convergence problems (first do a 436 !% calculation within the LCAO subspace, then restart from that point for 437 !% an unrestricted calculation). 438 !%End 439 call parse_variable(namespace, 'SCFinLCAO', .false., scf%lcao_restricted) 440 if(scf%lcao_restricted) then 441 call messages_experimental('SCFinLCAO') 442 message(1) = 'Info: SCF restricted to LCAO subspace.' 443 call messages_info(1) 444 445 if(scf%conv_eigen_error) then 446 message(1) = "ConvEigenError cannot be used with SCFinLCAO, since error is unknown." 447 call messages_fatal(1) 448 end if 449 end if 450 451 452 !%Variable SCFCalculateForces 453 !%Type logical 454 !%Section SCF 455 !%Description 456 !% This variable controls whether the forces on the ions are 457 !% calculated at the end of a self-consistent iteration. The 458 !% default is yes, unless the system only has user-defined 459 !% species. 460 !%End 461 call parse_variable(namespace, 'SCFCalculateForces', .not. geo%only_user_def, scf%calc_force) 462 463 if(scf%calc_force .and. gr%der%boundaries%spiralBC) then 464 message(1) = 'Forces cannot be calculated when using spiral boundary conditions.' 465 write(message(2),'(a)') 'Please use SCFCalculateForces = no.' 466 call messages_fatal(2, namespace=namespace) 467 end if 468 if(scf%calc_force) then 469 if(associated(hm%ep%B_field) .or. associated(hm%ep%A_static)) then 470 write(message(1),'(a)') 'The forces are currently not properly calculated if static' 471 write(message(2),'(a)') 'magnetic fields or static vector potentials are present.' 472 write(message(3),'(a)') 'Please use SCFCalculateForces = no.' 473 call messages_fatal(3, namespace=namespace) 474 end if 475 end if 476 477 !%Variable SCFCalculateStress 478 !%Type logical 479 !%Section SCF 480 !%Description 481 !% This variable controls whether the stress on the lattice is 482 !% calculated at the end of a self-consistent iteration. The 483 !% default is no. 484 !%End 485 call parse_variable(namespace, 'SCFCalculateStress', .false. , scf%calc_stress) 486 if(scf%calc_stress) then 487 if(ks%theory_level == HARTREE .or. ks%theory_level == HARTREE_FOCK .or. & 488 (ks%theory_level == KOHN_SHAM_DFT .and. bitand(hm%xc%family, XC_FAMILY_LDA) == 0)) then 489 write(message(1),'(a)') 'The stress tensor is currently only properly computed at the DFT-LDA level' 490 write(message(2),'(a)') 'Please use SCFCalculateStress = no.' 491 call messages_fatal(2, namespace=namespace) 492 end if 493 if(ks%vdw_correction /= OPTION__VDWCORRECTION__NONE) then 494 write(message(1),'(a)') 'The stress tensor is currently not properly computed with vdW corrections' 495 write(message(2),'(a)') 'Please use SCFCalculateStress = no.' 496 call messages_fatal(2, namespace=namespace) 497 end if 498 end if 499 500 !%Variable SCFCalculateDipole 501 !%Type logical 502 !%Section SCF 503 !%Description 504 !% This variable controls whether the dipole is calculated at the 505 !% end of a self-consistent iteration. For finite systems the 506 !% default is yes. For periodic systems the default is no, unless 507 !% an electric field is being applied in a periodic direction. 508 !% The single-point Berry`s phase approximation is used for 509 !% periodic directions. Ref: 510 !% E Yaschenko, L Fu, L Resca, and R Resta, <i>Phys. Rev. B</i> <b>58</b>, 1222-1229 (1998). 511 !%End 512 call parse_variable(namespace, 'SCFCalculateDipole', .not. simul_box_is_periodic(gr%sb), scf%calc_dipole) 513 if(associated(hm%vberry)) scf%calc_dipole = .true. 514 515 !%Variable SCFCalculatePartialCharges 516 !%Type logical 517 !%Default no 518 !%Section SCF 519 !%Description 520 !% (Experimental) This variable controls whether partial charges 521 !% are calculated at the end of a self-consistent iteration. 522 !%End 523 call parse_variable(namespace, 'SCFCalculatePartialCharges', .false., scf%calc_partial_charges) 524 if(scf%calc_partial_charges) call messages_experimental('SCFCalculatePartialCharges') 525 526 rmin = geometry_min_distance(geo) 527 if(geo%natoms == 1) then 528 if(simul_box_is_periodic(gr%sb)) then 529 rmin = minval(gr%sb%lsize(1:gr%sb%periodic_dim)) 530 else 531 rmin = CNST(100.0) 532 end if 533 end if 534 535 !%Variable LocalMagneticMomentsSphereRadius 536 !%Type float 537 !%Section Output 538 !%Description 539 !% The local magnetic moments are calculated by integrating the 540 !% magnetization density in spheres centered around each atom. 541 !% This variable controls the radius of the spheres. 542 !% The default is half the minimum distance between two atoms 543 !% in the input coordinates, or 100 a.u. if there is only one atom (for isolated systems). 544 !%End 545 call parse_variable(namespace, 'LocalMagneticMomentsSphereRadius', rmin*M_HALF, scf%lmm_r, unit = units_inp%length) 546 ! this variable is also used in td/td_write.F90 547 548 scf%forced_finish = .false. 549 550 scf%gr => gr 551 552 POP_SUB(scf_init) 553 end subroutine scf_init 554 555 556 ! --------------------------------------------------------- 557 subroutine scf_end(scf) 558 type(scf_t), intent(inout) :: scf 559 560 PUSH_SUB(scf_end) 561 562 call eigensolver_end(scf%eigens, scf%gr) 563 564 if(scf%mix_field /= OPTION__MIXFIELD__NONE) call mix_end(scf%smix) 565 566 nullify(scf%mixfield) 567 568 if(scf%mix_field /= OPTION__MIXFIELD__STATES) call lda_u_mixer_end(scf%lda_u_mix, scf%smix) 569 570 571 POP_SUB(scf_end) 572 end subroutine scf_end 573 574 575 ! --------------------------------------------------------- 576 subroutine scf_mix_clear(scf) 577 type(scf_t), intent(inout) :: scf 578 579 PUSH_SUB(scf_mix_clear) 580 581 call mix_clear(scf%smix) 582 583 if(scf%mix_field /= OPTION__MIXFIELD__STATES) call lda_u_mixer_clear(scf%lda_u_mix, scf%smix) 584 585 POP_SUB(scf_mix_clear) 586 end subroutine scf_mix_clear 587 588 589 ! --------------------------------------------------------- 590 subroutine scf_run(scf, namespace, mc, gr, geo, st, ks, hm, outp, gs_run, verbosity, iters_done, & 591 restart_load, restart_dump) 592 type(scf_t), intent(inout) :: scf !< self consistent cycle 593 type(namespace_t), intent(in) :: namespace 594 type(multicomm_t), intent(in) :: mc 595 type(grid_t), intent(inout) :: gr !< grid 596 type(geometry_t), intent(inout) :: geo !< geometry 597 type(states_elec_t), intent(inout) :: st !< States 598 type(v_ks_t), intent(inout) :: ks !< Kohn-Sham 599 type(hamiltonian_elec_t), intent(inout) :: hm !< Hamiltonian 600 type(output_t), intent(in) :: outp 601 logical, optional, intent(in) :: gs_run 602 integer, optional, intent(in) :: verbosity 603 integer, optional, intent(out) :: iters_done 604 type(restart_t), optional, intent(in) :: restart_load 605 type(restart_t), optional, intent(in) :: restart_dump 606 607 logical :: finish, converged_current, converged_last, gs_run_, berry_conv 608 integer :: iter, is, iatom, nspin, ierr, iberry, idir, verbosity_, ib, iqn 609 FLOAT :: evsum_out, evsum_in, forcetmp, dipole(MAX_DIM), dipole_prev(MAX_DIM) 610 FLOAT :: etime, itime 611 character(len=MAX_PATH_LEN) :: dirname 612 type(lcao_t) :: lcao !< Linear combination of atomic orbitals 613 type(profile_t), save :: prof 614 FLOAT, allocatable :: rhoout(:,:,:), rhoin(:,:,:) 615 FLOAT, allocatable :: vhxc_old(:,:) 616 FLOAT, allocatable :: forceout(:,:), forcein(:,:), forcediff(:), tmp(:) 617 class(wfs_elec_t), allocatable :: psioutb(:, :) 618#ifdef HAVE_MPI 619 logical :: forced_finish_tmp 620#endif 621 622 PUSH_SUB(scf_run) 623 624 if(scf%forced_finish) then 625 message(1) = "Previous clean stop, not doing SCF and quitting." 626 call messages_fatal(1, only_root_writes = .true.) 627 end if 628 629 gs_run_ = .true. 630 if(present(gs_run)) gs_run_ = gs_run 631 632 verbosity_ = VERB_FULL 633 if(present(verbosity)) verbosity_ = verbosity 634 635 if(scf%lcao_restricted) then 636 call lcao_init(lcao, namespace, gr, geo, st) 637 if(.not. lcao_is_available(lcao)) then 638 message(1) = 'LCAO is not available. Cannot do SCF in LCAO.' 639 call messages_fatal(1) 640 end if 641 end if 642 643 nspin = st%d%nspin 644 645 if (present(restart_load)) then 646 if (restart_has_flag(restart_load, RESTART_FLAG_RHO)) then 647 ! Load density and used it to recalculated the KS potential. 648 call states_elec_load_rho(restart_load, st, gr, ierr) 649 if (ierr /= 0) then 650 message(1) = 'Unable to read density. Density will be calculated from states.' 651 call messages_warning(1) 652 else 653 if(bitand(ks%xc_family, XC_FAMILY_OEP) == 0) then 654 call v_ks_calc(ks, namespace, hm, st, geo) 655 else 656 if (.not. restart_has_flag(restart_load, RESTART_FLAG_VHXC) .and. ks%oep%level /= XC_OEP_FULL) then 657 call v_ks_calc(ks, namespace, hm, st, geo) 658 end if 659 end if 660 end if 661 end if 662 663 if (restart_has_flag(restart_load, RESTART_FLAG_VHXC)) then 664 call hamiltonian_elec_load_vhxc(restart_load, hm, gr%mesh, ierr) 665 if (ierr /= 0) then 666 message(1) = 'Unable to read Vhxc. Vhxc will be calculated from states.' 667 call messages_warning(1) 668 else 669 call hamiltonian_elec_update(hm, gr%mesh, namespace) 670 if(bitand(ks%xc_family, XC_FAMILY_OEP) /= 0) then 671 if (ks%oep%level == XC_OEP_FULL) then 672 do is = 1, st%d%nspin 673 ks%oep%vxc(1:gr%mesh%np, is) = hm%vhxc(1:gr%mesh%np, is) - hm%vhartree(1:gr%mesh%np) 674 end do 675 call v_ks_calc(ks, namespace, hm, st, geo) 676 end if 677 end if 678 end if 679 end if 680 681 if (restart_has_flag(restart_load, RESTART_FLAG_MIX)) then 682 select case (scf%mix_field) 683 case (OPTION__MIXFIELD__DENSITY) 684 call mix_load(restart_load, scf%smix, gr%fine%mesh, ierr) 685 case (OPTION__MIXFIELD__POTENTIAL) 686 call mix_load(restart_load, scf%smix, gr%mesh, ierr) 687 end select 688 if (ierr /= 0) then 689 message(1) = "Unable to read mixing information. Mixing will start from scratch." 690 call messages_warning(1) 691 end if 692 end if 693 694 if(hm%lda_u_level /= DFT_U_NONE) then 695 call lda_u_load(restart_load, hm%lda_u, st, hm%energy%dft_u, ierr) 696 if (ierr /= 0) then 697 message(1) = "Unable to read LDA+U information. LDA+U data will be calculated from states." 698 call messages_warning(1) 699 end if 700 end if 701 end if 702 703 SAFE_ALLOCATE(rhoout(1:gr%fine%mesh%np, 1:1, 1:nspin)) 704 SAFE_ALLOCATE(rhoin (1:gr%fine%mesh%np, 1:1, 1:nspin)) 705 706 rhoin(1:gr%fine%mesh%np, 1, 1:nspin) = st%rho(1:gr%fine%mesh%np, 1:nspin) 707 rhoout = M_ZERO 708 709 !We store the Hxc potential for the contribution to the forces 710 if(scf%calc_force .or. scf%conv_abs_force > M_ZERO & 711 .or. (outp%duringscf .and. bitand(outp%what, OPTION__OUTPUT__FORCES) /= 0)) then 712 SAFE_ALLOCATE(vhxc_old(1:gr%mesh%np, 1:nspin)) 713 vhxc_old(1:gr%mesh%np, 1:nspin) = hm%vhxc(1:gr%mesh%np, 1:nspin) 714 end if 715 716 717 select case(scf%mix_field) 718 case(OPTION__MIXFIELD__POTENTIAL) 719 call mixfield_set_vin(scf%mixfield, hm%vhxc) 720 case(OPTION__MIXFIELD__DENSITY) 721 call mixfield_set_vin(scf%mixfield, rhoin) 722 723 case(OPTION__MIXFIELD__STATES) 724 725 SAFE_ALLOCATE_TYPE_ARRAY(wfs_elec_t, psioutb, (st%group%block_start:st%group%block_end, st%d%kpt%start:st%d%kpt%end)) 726 727 do iqn = st%d%kpt%start, st%d%kpt%end 728 do ib = st%group%block_start, st%group%block_end 729 call st%group%psib(ib, iqn)%copy_to(psioutb(ib, iqn)) 730 end do 731 end do 732 733 end select 734 735 call lda_u_update_occ_matrices(hm%lda_u, namespace, gr%mesh, st, hm%hm_base, hm%energy) 736 !If we use LDA+U, we also have do mix it 737 if(scf%mix_field /= OPTION__MIXFIELD__STATES) call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix) 738 739 evsum_in = states_elec_eigenvalues_sum(st) 740 741 ! allocate and compute forces only if they are used as convergence criteria 742 if (scf%conv_abs_force > M_ZERO) then 743 SAFE_ALLOCATE( forcein(1:geo%natoms, 1:gr%sb%dim)) 744 SAFE_ALLOCATE( forceout(1:geo%natoms, 1:gr%sb%dim)) 745 SAFE_ALLOCATE(forcediff(1:gr%sb%dim)) 746 call forces_calculate(gr, namespace, geo, hm, st, ks) 747 do iatom = 1, geo%natoms 748 forcein(iatom, 1:gr%sb%dim) = geo%atom(iatom)%f(1:gr%sb%dim) 749 end do 750 end if 751 752 call create_convergence_file(STATIC_DIR, "convergence") 753 754 if ( verbosity_ /= VERB_NO ) then 755 if(scf%max_iter > 0) then 756 write(message(1),'(a)') 'Info: Starting SCF iteration.' 757 else 758 write(message(1),'(a)') 'Info: No SCF iterations will be done.' 759 ! we cannot tell whether it is converged. 760 finish = .false. 761 end if 762 call messages_info(1) 763 end if 764 converged_current = .false. 765 766 ! SCF cycle 767 itime = loct_clock() 768 769 do iter = 1, scf%max_iter 770 call profiling_in(prof, "SCF_CYCLE") 771 772 ! reset scdm flag 773 scdm_is_local = .false. 774 775 ! this initialization seems redundant but avoids improper optimization at -O3 by PGI 7 on chum, 776 ! which would cause a failure of testsuite/linear_response/04-vib_modes.03-vib_modes_fd.inp 777 scf%eigens%converged = 0 778 779 scf%energy_diff = hm%energy%total 780 781 !Used for the contribution to the forces 782 if(scf%calc_force .or. scf%conv_abs_force > M_ZERO .or. & 783 (outp%duringscf .and. bitand(outp%what, OPTION__OUTPUT__FORCES) /= 0)) & 784 vhxc_old(1:gr%mesh%np, 1:nspin) = hm%vhxc(1:gr%mesh%np, 1:nspin) 785 786 if(scf%lcao_restricted) then 787 call lcao_init_orbitals(lcao, st, gr, geo) 788 call lcao_wf(lcao, st, gr, geo, hm, namespace) 789 else 790 if(associated(hm%vberry)) then 791 ks%frozen_hxc = .true. 792 do iberry = 1, scf%max_iter_berry 793 scf%eigens%converged = 0 794 call eigensolver_run(scf%eigens, namespace, gr, st, hm, iter) 795 796 call v_ks_calc(ks, namespace, hm, st, geo, calc_current=outp%duringscf) 797 798 dipole_prev = dipole 799 call calc_dipole(dipole) 800 write(message(1),'(a,9f12.6)') 'Dipole = ', dipole(1:gr%sb%dim) 801 call messages_info(1) 802 803 berry_conv = .true. 804 do idir = 1, gr%sb%periodic_dim 805 berry_conv = berry_conv .and. & 806 (abs((dipole(idir) - dipole_prev(idir)) / dipole_prev(idir)) < CNST(1e-5) & 807 .or. abs(dipole(idir) - dipole_prev(idir)) < CNST(1e-5)) 808 end do 809 if(berry_conv) exit 810 end do 811 ks%frozen_hxc = .false. 812 else 813 scf%eigens%converged = 0 814 call eigensolver_run(scf%eigens, namespace, gr, st, hm, iter) 815 end if 816 end if 817 818 ! occupations 819 call states_elec_fermi(st, namespace, gr%mesh) 820 call lda_u_update_occ_matrices(hm%lda_u, namespace, gr%mesh, st, hm%hm_base, hm%energy) 821 822 ! compute output density, potential (if needed) and eigenvalues sum 823 call density_calc(st, gr, st%rho) 824 825 rhoout(1:gr%fine%mesh%np, 1, 1:nspin) = st%rho(1:gr%fine%mesh%np, 1:nspin) 826 827 select case(scf%mix_field) 828 case(OPTION__MIXFIELD__POTENTIAL) 829 call v_ks_calc(ks, namespace, hm, st, geo, calc_current=outp%duringscf) 830 call mixfield_set_vout(scf%mixfield, hm%vhxc) 831 case (OPTION__MIXFIELD__DENSITY) 832 call mixfield_set_vout(scf%mixfield, rhoout) 833 case(OPTION__MIXFIELD__STATES) 834 835 do iqn = st%d%kpt%start, st%d%kpt%end 836 do ib = st%group%block_start, st%group%block_end 837 call st%group%psib(ib, iqn)%copy_data_to(gr%mesh%np, psioutb(ib, iqn)) 838 end do 839 end do 840 end select 841 842 if(scf%mix_field /= OPTION__MIXFIELD__STATES) call lda_u_mixer_set_vout(hm%lda_u, scf%lda_u_mix) 843 844 evsum_out = states_elec_eigenvalues_sum(st) 845 846 ! recalculate total energy 847 call energy_calc_total(namespace, hm, gr, st, iunit = 0) 848 849 ! compute convergence criteria 850 scf%energy_diff = hm%energy%total - scf%energy_diff 851 scf%abs_dens = M_ZERO 852 SAFE_ALLOCATE(tmp(1:gr%fine%mesh%np)) 853 do is = 1, nspin 854 tmp = abs(rhoin(1:gr%fine%mesh%np, 1, is) - rhoout(1:gr%fine%mesh%np, 1, is)) 855 scf%abs_dens = scf%abs_dens + dmf_integrate(gr%fine%mesh, tmp) 856 end do 857 SAFE_DEALLOCATE_A(tmp) 858 859 ! compute forces only if they are used as convergence criterion 860 if (scf%conv_abs_force > M_ZERO) then 861 call forces_calculate(gr, namespace, geo, hm, st, ks, vhxc_old=vhxc_old) 862 scf%abs_force = M_ZERO 863 do iatom = 1, geo%natoms 864 forceout(iatom,1:gr%sb%dim) = geo%atom(iatom)%f(1:gr%sb%dim) 865 forcediff(1:gr%sb%dim) = abs( forceout(iatom,1:gr%sb%dim) - forcein(iatom,1:gr%sb%dim) ) 866 forcetmp = maxval( forcediff ) 867 if ( forcetmp > scf%abs_force ) then 868 scf%abs_force = forcetmp 869 end if 870 end do 871 else 872 if(outp%duringscf .and. bitand(outp%what, OPTION__OUTPUT__FORCES) /= 0 & 873 .and. outp%output_interval /= 0 & 874 .and. gs_run_ .and. mod(iter, outp%output_interval) == 0) & 875 call forces_calculate(gr, namespace, geo, hm, st, ks, vhxc_old=vhxc_old) 876 end if 877 878 if(abs(st%qtot) <= M_EPSILON) then 879 scf%rel_dens = M_HUGE 880 else 881 scf%rel_dens = scf%abs_dens / st%qtot 882 end if 883 884 scf%abs_ev = abs(evsum_out - evsum_in) 885 if(abs(evsum_out) <= M_EPSILON) then 886 scf%rel_ev = M_HUGE 887 else 888 scf%rel_ev = scf%abs_ev / abs(evsum_out) 889 end if 890 891 scf%eigens%current_rel_dens_error = scf%rel_dens 892 893 ! are we finished? 894 converged_last = converged_current 895 converged_current = scf%check_conv .and. & 896 (scf%conv_abs_dens <= M_ZERO .or. scf%abs_dens <= scf%conv_abs_dens) .and. & 897 (scf%conv_rel_dens <= M_ZERO .or. scf%rel_dens <= scf%conv_rel_dens) .and. & 898 (scf%conv_abs_force <= M_ZERO .or. scf%abs_force <= scf%conv_abs_force) .and. & 899 (scf%conv_abs_ev <= M_ZERO .or. scf%abs_ev <= scf%conv_abs_ev) .and. & 900 (scf%conv_rel_ev <= M_ZERO .or. scf%rel_ev <= scf%conv_rel_ev) .and. & 901 (scf%conv_energy_diff <= M_ZERO .or. abs(scf%energy_diff) <= scf%conv_energy_diff) .and. & 902 (.not. scf%conv_eigen_error .or. all(scf%eigens%converged == st%nst)) 903 ! only finish if the convergence criterion is fulfilled in two 904 ! consecutive iterations 905 finish = converged_last .and. converged_current 906 907 etime = loct_clock() - itime 908 itime = etime + itime 909 call scf_write_iter() 910 911 ! mixing 912 select case (scf%mix_field) 913 case (OPTION__MIXFIELD__DENSITY) 914 ! mix input and output densities and compute new potential 915 call mixing(scf%smix) 916 call mixfield_get_vnew(scf%mixfield, st%rho) 917 ! for spinors, having components 3 or 4 be negative is not unphysical 918 if(minval(st%rho(1:gr%fine%mesh%np, 1:st%d%spin_channels)) < -CNST(1e-6)) then 919 write(message(1),*) 'Negative density after mixing. Minimum value = ', & 920 minval(st%rho(1:gr%fine%mesh%np, 1:st%d%spin_channels)) 921 call messages_warning(1) 922 end if 923 call lda_u_mixer_get_vnew(hm%lda_u, scf%lda_u_mix, st) 924 call v_ks_calc(ks, namespace, hm, st, geo, calc_current=outp%duringscf) 925 case (OPTION__MIXFIELD__POTENTIAL) 926 ! mix input and output potentials 927 call mixing(scf%smix) 928 call mixfield_get_vnew(scf%mixfield, hm%vhxc) 929 call lda_u_mixer_get_vnew(hm%lda_u, scf%lda_u_mix, st) 930 call hamiltonian_elec_update_pot(hm, gr%mesh, accel_copy=.true.) 931 932 case(OPTION__MIXFIELD__STATES) 933 934 do iqn = st%d%kpt%start, st%d%kpt%end 935 do ib = st%group%block_start, st%group%block_end 936 call batch_scal(gr%mesh%np, CNST(1.0) - mix_coefficient(scf%smix), st%group%psib(ib, iqn)) 937 call batch_axpy(gr%mesh%np, mix_coefficient(scf%smix), psioutb(ib, iqn), st%group%psib(ib, iqn)) 938 end do 939 end do 940 941 call density_calc(st, gr, st%rho) 942 call v_ks_calc(ks, namespace, hm, st, geo, calc_current=outp%duringscf) 943 944 case(OPTION__MIXFIELD__NONE) 945 call v_ks_calc(ks, namespace, hm, st, geo, calc_current=outp%duringscf) 946 end select 947 948 949 ! Are we asked to stop? (Whenever Fortran is ready for signals, this should go away) 950 scf%forced_finish = clean_stop(mc%master_comm) .or. walltimer_alarm() 951 952#ifdef HAVE_MPI 953 call MPI_Allreduce(scf%forced_finish, forced_finish_tmp, 1, MPI_LOGICAL, MPI_LOR, mc%master_comm, mpi_err) 954 scf%forced_finish = forced_finish_tmp 955#endif 956 957 if (finish .and. st%modelmbparticles%nparticle > 0) then 958 call modelmb_sym_all_states (gr, st) 959 end if 960 961 if (gs_run_ .and. present(restart_dump)) then 962 ! save restart information 963 964 if ( (finish .or. (modulo(iter, outp%restart_write_interval) == 0) & 965 .or. iter == scf%max_iter .or. scf%forced_finish) ) then 966 967 call states_elec_dump(restart_dump, st, gr, ierr, iter=iter) 968 if (ierr /= 0) then 969 message(1) = 'Unable to write states wavefunctions.' 970 call messages_warning(1) 971 end if 972 973 call states_elec_dump_rho(restart_dump, st, gr, ierr, iter=iter) 974 if (ierr /= 0) then 975 message(1) = 'Unable to write density.' 976 call messages_warning(1) 977 end if 978 979 if(hm%lda_u_level /= DFT_U_NONE) then 980 call lda_u_dump(restart_dump, hm%lda_u, st, ierr) 981 if (ierr /= 0) then 982 message(1) = 'Unable to write LDA+U information.' 983 call messages_warning(1) 984 end if 985 end if 986 987 select case (scf%mix_field) 988 case (OPTION__MIXFIELD__DENSITY) 989 call mix_dump(restart_dump, scf%smix, gr%fine%mesh, ierr) 990 if (ierr /= 0) then 991 message(1) = 'Unable to write mixing information.' 992 call messages_warning(1) 993 end if 994 case (OPTION__MIXFIELD__POTENTIAL) 995 call hamiltonian_elec_dump_vhxc(restart_dump, hm, gr%mesh, ierr) 996 if (ierr /= 0) then 997 message(1) = 'Unable to write Vhxc.' 998 call messages_warning(1) 999 end if 1000 1001 call mix_dump(restart_dump, scf%smix, gr%mesh, ierr) 1002 if (ierr /= 0) then 1003 message(1) = 'Unable to write mixing information.' 1004 call messages_warning(1) 1005 end if 1006 end select 1007 end if 1008 end if 1009 1010 call write_convergence_file(STATIC_DIR, "convergence") 1011 1012 if(finish) then 1013 if(present(iters_done)) iters_done = iter 1014 if(verbosity_ >= VERB_COMPACT) then 1015 write(message(1), '(a, i4, a)') 'Info: SCF converged in ', iter, ' iterations' 1016 write(message(2), '(a)') '' 1017 call messages_info(2) 1018 end if 1019 call profiling_out(prof) 1020 exit 1021 end if 1022 1023 if((outp%what+outp%what_lda_u+outp%whatBZ)/=0 .and. outp%duringscf .and. outp%output_interval /= 0 & 1024 .and. gs_run_ .and. mod(iter, outp%output_interval) == 0) then 1025 write(dirname,'(a,a,i4.4)') trim(outp%iter_dir),"scf.",iter 1026 call output_all(outp, namespace, dirname, gr, geo, st, hm, ks) 1027 end if 1028 1029 ! save information for the next iteration 1030 rhoin(1:gr%fine%mesh%np, 1, 1:nspin) = st%rho(1:gr%fine%mesh%np, 1:nspin) 1031 1032 select case(scf%mix_field) 1033 case(OPTION__MIXFIELD__POTENTIAL) 1034 call mixfield_set_vin(scf%mixfield, hm%vhxc(1:gr%mesh%np, 1:nspin)) 1035 case (OPTION__MIXFIELD__DENSITY) 1036 call mixfield_set_vin(scf%mixfield, rhoin) 1037 end select 1038 if(scf%mix_field /= OPTION__MIXFIELD__STATES) call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix) 1039 1040 evsum_in = evsum_out 1041 if (scf%conv_abs_force > M_ZERO) then 1042 forcein(1:geo%natoms, 1:gr%sb%dim) = forceout(1:geo%natoms, 1:gr%sb%dim) 1043 end if 1044 1045 1046 if(scf%forced_finish) then 1047 call profiling_out(prof) 1048 exit 1049 end if 1050 1051 ! check if debug mode should be enabled or disabled on the fly 1052 call io_debug_on_the_fly(namespace) 1053 1054 call profiling_out(prof) 1055 end do !iter 1056 1057 if(scf%lcao_restricted) call lcao_end(lcao) 1058 1059 if((scf%max_iter > 0 .and. scf%mix_field == OPTION__MIXFIELD__POTENTIAL) .or. output_needs_current(outp, states_are_real(st))) then 1060 call v_ks_calc(ks, namespace, hm, st, geo) 1061 end if 1062 1063 select case(scf%mix_field) 1064 case(OPTION__MIXFIELD__STATES) 1065 1066 do iqn = st%d%kpt%start, st%d%kpt%end 1067 do ib = st%group%block_start, st%group%block_end 1068 call psioutb(ib, iqn)%end() 1069 end do 1070 end do 1071 1072 SAFE_DEALLOCATE_A(psioutb) 1073 end select 1074 1075 SAFE_DEALLOCATE_A(rhoout) 1076 SAFE_DEALLOCATE_A(rhoin) 1077 1078 if(scf%max_iter > 0 .and. any(scf%eigens%converged < st%nst) .and. .not. scf%lcao_restricted) then 1079 write(message(1),'(a)') 'Some of the states are not fully converged!' 1080 call messages_warning(1) 1081 end if 1082 1083 if(.not.finish) then 1084 write(message(1), '(a,i4,a)') 'SCF *not* converged after ', iter - 1, ' iterations.' 1085 call messages_warning(1) 1086 end if 1087 1088 ! calculate forces 1089 if(scf%calc_force) then 1090 call forces_calculate(gr, namespace, geo, hm, st, ks, vhxc_old=vhxc_old) 1091 end if 1092 1093 ! calculate stress 1094 if(scf%calc_stress) call stress_calculate(namespace, gr, hm, st, geo, ks) 1095 1096 if(scf%max_iter == 0) then 1097 call energy_calc_eigenvalues(namespace, hm, gr%der, st) 1098 call states_elec_fermi(st, namespace, gr%mesh) 1099 call states_elec_write_eigenvalues(stdout, st%nst, st, gr%sb) 1100 end if 1101 1102 if(gs_run_) then 1103 ! output final information 1104 call scf_write_static(STATIC_DIR, "info") 1105 call output_all(outp, namespace, STATIC_DIR, gr, geo, st, hm, ks) 1106 end if 1107 1108 if(simul_box_is_periodic(gr%sb) .and. st%d%nik > st%d%nspin) then 1109 if(bitand(gr%sb%kpoints%method, KPOINTS_PATH) /= 0) then 1110 call states_elec_write_bandstructure(STATIC_DIR, namespace, st%nst, st, gr%sb, & 1111 geo, gr%mesh, hm%hm_base%phase, vec_pot = hm%hm_base%uniform_vector_potential, & 1112 vec_pot_var = hm%hm_base%vector_potential) 1113 end if 1114 end if 1115 1116 if( ks%vdw_correction == OPTION__VDWCORRECTION__VDW_TS) then 1117 call vdw_ts_write_c6ab(ks%vdw_ts, geo, STATIC_DIR, 'c6ab_eff', namespace) 1118 end if 1119 1120 SAFE_DEALLOCATE_A(vhxc_old) 1121 1122 POP_SUB(scf_run) 1123 1124 contains 1125 1126 1127 ! --------------------------------------------------------- 1128 subroutine scf_write_iter() 1129 character(len=50) :: str 1130 1131 PUSH_SUB(scf_run.scf_write_iter) 1132 1133 if ( verbosity_ == VERB_FULL ) then 1134 1135 write(str, '(a,i5)') 'SCF CYCLE ITER #' ,iter 1136 call messages_print_stress(stdout, trim(str)) 1137 write(message(1),'(a,es15.8,2(a,es9.2))') ' etot = ', units_from_atomic(units_out%energy, hm%energy%total), & 1138 ' abs_ev = ', units_from_atomic(units_out%energy, scf%abs_ev), ' rel_ev = ', scf%rel_ev 1139 write(message(2),'(a,es15.2,2(a,es9.2))') & 1140 ' ediff = ', scf%energy_diff, ' abs_dens = ', scf%abs_dens, ' rel_dens = ', scf%rel_dens 1141 ! write info about forces only if they are used as convergence criteria 1142 if (scf%conv_abs_force > M_ZERO) then 1143 write(message(3),'(23x,a,es9.2)') & 1144 ' force = ', units_from_atomic(units_out%force, scf%abs_force) 1145 call messages_info(3) 1146 else 1147 call messages_info(2) 1148 end if 1149 1150 if(.not.scf%lcao_restricted) then 1151 write(message(1),'(a,i6)') 'Matrix vector products: ', scf%eigens%matvec 1152 write(message(2),'(a,i6)') 'Converged eigenvectors: ', sum(scf%eigens%converged(1:st%d%nik)) 1153 call messages_info(2) 1154 call states_elec_write_eigenvalues(stdout, st%nst, st, gr%sb, scf%eigens%diff, compact = .true.) 1155 else 1156 call states_elec_write_eigenvalues(stdout, st%nst, st, gr%sb, compact = .true.) 1157 end if 1158 1159 if(associated(hm%vberry)) then 1160 call calc_dipole(dipole) 1161 call write_dipole(stdout, dipole) 1162 end if 1163 1164 if(st%d%ispin > UNPOLARIZED) then 1165 call write_magnetic_moments(stdout, gr%mesh, st, geo, gr%der%boundaries, scf%lmm_r) 1166 end if 1167 1168 if(hm%lda_u_level == DFT_U_ACBN0) then 1169 call lda_u_write_U(hm%lda_u, stdout) 1170 call lda_u_write_V(hm%lda_u, stdout) 1171 end if 1172 1173 write(message(1),'(a)') '' 1174 write(message(2),'(a,i5,a,f14.2)') 'Elapsed time for SCF step ', iter,':', etime 1175 call messages_info(2) 1176 1177 call scf_print_mem_use() 1178 1179 call messages_print_stress(stdout) 1180 1181 end if 1182 1183 if ( verbosity_ == VERB_COMPACT ) then 1184 ! write info about forces only if they are used as convergence criteria 1185 if (scf%conv_abs_force > M_ZERO) then 1186 write(message(1),'(a,i4,a,es15.8, 2(a,es9.2), a, f7.1, a)') & 1187 'iter ', iter, & 1188 ' : etot ', units_from_atomic(units_out%energy, hm%energy%total), & 1189 ' : abs_dens', scf%abs_dens, & 1190 ' : force ', units_from_atomic(units_out%force, scf%abs_force), & 1191 ' : etime ', etime, 's' 1192 else 1193 write(message(1),'(a,i4,a,es15.8, a,es9.2, a, f7.1, a)') & 1194 'iter ', iter, & 1195 ' : etot ', units_from_atomic(units_out%energy, hm%energy%total), & 1196 ' : abs_dens', scf%abs_dens, & 1197 ' : etime ', etime, 's' 1198 end if 1199 call messages_info(1) 1200 end if 1201 1202 POP_SUB(scf_run.scf_write_iter) 1203 end subroutine scf_write_iter 1204 1205 1206 ! --------------------------------------------------------- 1207 subroutine scf_write_static(dir, fname) 1208 character(len=*), intent(in) :: dir, fname 1209 1210 integer :: iunit, iatom 1211 FLOAT, allocatable :: hirshfeld_charges(:) 1212 1213 PUSH_SUB(scf_run.scf_write_static) 1214 1215 if(mpi_grp_is_root(mpi_world)) then ! this the absolute master writes 1216 call io_mkdir(dir, namespace) 1217 iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write') 1218 1219 call grid_write_info(gr, geo, iunit) 1220 1221 call symmetries_write_info(gr%mesh%sb%symm, namespace, gr%sb%dim, gr%sb%periodic_dim, iunit) 1222 1223 if(simul_box_is_periodic(gr%sb)) then 1224 call kpoints_write_info(gr%mesh%sb%kpoints, namespace, iunit) 1225 write(iunit,'(1x)') 1226 end if 1227 1228 call v_ks_write_info(ks, iunit, namespace) 1229 1230 ! scf information 1231 if(finish) then 1232 write(iunit, '(a, i4, a)')'SCF converged in ', iter, ' iterations' 1233 else 1234 write(iunit, '(a)') 'SCF *not* converged!' 1235 end if 1236 write(iunit, '(1x)') 1237 1238 if(any(scf%eigens%converged < st%nst) .and. .not. scf%lcao_restricted) then 1239 write(iunit,'(a)') 'Some of the states are not fully converged!' 1240 end if 1241 1242 call states_elec_write_eigenvalues(iunit, st%nst, st, gr%sb) 1243 write(iunit, '(1x)') 1244 1245 if(simul_box_is_periodic(gr%sb)) then 1246 call states_elec_write_gaps(iunit, st, gr%sb) 1247 write(iunit, '(1x)') 1248 end if 1249 1250 write(iunit, '(3a)') 'Energy [', trim(units_abbrev(units_out%energy)), ']:' 1251 else 1252 iunit = 0 1253 end if 1254 1255 call energy_calc_total(namespace, hm, gr, st, iunit, full = .true.) 1256 1257 if(mpi_grp_is_root(mpi_world)) write(iunit, '(1x)') 1258 if(st%d%ispin > UNPOLARIZED) then 1259 call write_magnetic_moments(iunit, gr%mesh, st, geo, gr%der%boundaries, scf%lmm_r) 1260 if(mpi_grp_is_root(mpi_world)) write(iunit, '(1x)') 1261 end if 1262 1263 if(hm%lda_u_level == DFT_U_ACBN0) then 1264 call lda_u_write_U(hm%lda_u, iunit) 1265 call lda_u_write_V(hm%lda_u, iunit) 1266 if(mpi_grp_is_root(mpi_world)) write(iunit, '(1x)') 1267 end if 1268 1269 if(scf%calc_dipole) then 1270 call calc_dipole(dipole) 1271 call write_dipole(iunit, dipole) 1272 end if 1273 1274 if(mpi_grp_is_root(mpi_world)) then 1275 if(scf%max_iter > 0) then 1276 write(iunit, '(a)') 'Convergence:' 1277 write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'abs_dens = ', scf%abs_dens, & 1278 ' (', scf%conv_abs_dens, ')' 1279 write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'rel_dens = ', scf%rel_dens, & 1280 ' (', scf%conv_rel_dens, ')' 1281 write(iunit, '(6x, a, es15.8,a,es15.8,4a)') 'abs_ev = ', & 1282 units_from_atomic(units_out%energy, scf%abs_ev), & 1283 ' (', units_from_atomic(units_out%energy, scf%conv_abs_ev), ')', & 1284 ' [', trim(units_abbrev(units_out%energy)), ']' 1285 write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'rel_ev = ', scf%rel_ev, & 1286 ' (', scf%conv_rel_ev, ')' 1287 write(iunit,'(1x)') 1288 end if 1289 ! otherwise, these values are uninitialized, and unknown. 1290 1291 if (bitand(ks%xc_family, XC_FAMILY_OEP) /= 0 .and. ks%theory_level /= HARTREE_FOCK) then 1292 if (((ks%oep%level == XC_OEP_FULL) .or. (ks%oep%level == XC_OEP_KLI)) .and. ks%oep%has_photons) then 1293 write(iunit, '(a)') 'Photon observables:' 1294 write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'Photon number = ', ks%oep%pt%number(1) 1295 write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'Photon ex. = ', ks%oep%pt%ex 1296 write(iunit,'(1x)') 1297 end if 1298 end if 1299 1300 if(scf%calc_force) call forces_write_info(iunit, geo, gr%sb, dir, namespace) 1301 1302 if(scf%calc_stress) then 1303 write(iunit,'(a)') "Stress tensor [H/b^3]" 1304 write(iunit,'(a9,2x,3a18)')"T_{ij}","x","y","z" 1305 write(iunit,'(a9,2x,3es18.6)')"x", gr%sb%stress_tensor(1,1:3) 1306 write(iunit,'(a9,2x,3es18.6)')"y", gr%sb%stress_tensor(2,1:3) 1307 write(iunit,'(a9,2x,3es18.6)')"z", gr%sb%stress_tensor(3,1:3) 1308 end if 1309 1310 end if 1311 1312 if(scf%calc_partial_charges) then 1313 SAFE_ALLOCATE(hirshfeld_charges(1:geo%natoms)) 1314 1315 call partial_charges_calculate(namespace, gr%fine%mesh, st, geo, hirshfeld_charges = hirshfeld_charges) 1316 1317 if(mpi_grp_is_root(mpi_world)) then 1318 1319 write(iunit,'(a)') 'Partial ionic charges' 1320 write(iunit,'(a)') ' Ion Hirshfeld' 1321 1322 do iatom = 1, geo%natoms 1323 write(iunit,'(i4,a10,f16.3)') iatom, trim(species_label(geo%atom(iatom)%species)), hirshfeld_charges(iatom) 1324 1325 end do 1326 1327 end if 1328 1329 SAFE_DEALLOCATE_A(hirshfeld_charges) 1330 1331 end if 1332 1333 if(mpi_grp_is_root(mpi_world)) then 1334 call io_close(iunit) 1335 end if 1336 1337 POP_SUB(scf_run.scf_write_static) 1338 end subroutine scf_write_static 1339 1340 1341 ! --------------------------------------------------------- 1342 subroutine calc_dipole(dipole) 1343 FLOAT, intent(out) :: dipole(:) 1344 1345 integer :: ispin, idir 1346 FLOAT :: e_dip(MAX_DIM + 1, st%d%nspin), n_dip(MAX_DIM), nquantumpol 1347 1348 PUSH_SUB(scf_run.calc_dipole) 1349 1350 dipole(1:MAX_DIM) = M_ZERO 1351 1352 do ispin = 1, st%d%nspin 1353 call dmf_multipoles(gr%fine%mesh, st%rho(:, ispin), 1, e_dip(:, ispin)) 1354 end do 1355 1356 call geometry_dipole(geo, n_dip) 1357 1358 do idir = 1, gr%sb%dim 1359 ! in periodic directions use single-point Berry`s phase calculation 1360 if(idir <= gr%sb%periodic_dim) then 1361 dipole(idir) = -n_dip(idir) - berry_dipole(st, gr%mesh, idir) 1362 1363 ! use quantum of polarization to reduce to smallest possible magnitude 1364 nquantumpol = NINT(dipole(idir)/(CNST(2.0)*gr%sb%lsize(idir))) 1365 dipole(idir) = dipole(idir) - nquantumpol * (CNST(2.0) * gr%sb%lsize(idir)) 1366 1367 ! in aperiodic directions use normal dipole formula 1368 else 1369 e_dip(idir + 1, 1) = sum(e_dip(idir + 1, :)) 1370 dipole(idir) = -n_dip(idir) - e_dip(idir + 1, 1) 1371 end if 1372 end do 1373 1374 POP_SUB(scf_run.calc_dipole) 1375 end subroutine calc_dipole 1376 1377 1378 ! --------------------------------------------------------- 1379 subroutine write_dipole(iunit, dipole) 1380 integer, intent(in) :: iunit 1381 FLOAT, intent(in) :: dipole(:) 1382 1383 PUSH_SUB(scf_run.write_dipole) 1384 1385 if(mpi_grp_is_root(mpi_world)) then 1386 call output_dipole(iunit, dipole, gr%mesh%sb%dim) 1387 1388 if (simul_box_is_periodic(gr%sb)) then 1389 write(iunit, '(a)') "Defined only up to quantum of polarization (e * lattice vector)." 1390 write(iunit, '(a)') "Single-point Berry's phase method only accurate for large supercells." 1391 1392 if (gr%sb%kpoints%full%npoints > 1) then 1393 write(iunit, '(a)') & 1394 "WARNING: Single-point Berry's phase method for dipole should not be used when there is more than one k-point." 1395 write(iunit, '(a)') & 1396 "Instead, finite differences on k-points (not yet implemented) are needed." 1397 end if 1398 1399 if(.not. smear_is_semiconducting(st%smear)) then 1400 write(iunit, '(a)') "Single-point Berry's phase dipole calculation not correct without integer occupations." 1401 end if 1402 end if 1403 1404 write(iunit, *) 1405 end if 1406 1407 POP_SUB(scf_run.write_dipole) 1408 end subroutine write_dipole 1409 1410 ! ----------------------------------------------------- 1411 1412 subroutine create_convergence_file(dir, fname) 1413 character(len=*), intent(in) :: dir 1414 character(len=*), intent(in) :: fname 1415 1416 integer :: iunit 1417 character(len=12) :: label 1418 if(mpi_grp_is_root(mpi_world)) then ! this the absolute master writes 1419 call io_mkdir(dir, namespace) 1420 iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write') 1421 write(iunit, '(a)', advance = 'no') '#iter energy ' 1422 label = 'energy_diff' 1423 write(iunit, '(1x,a)', advance = 'no') label 1424 label = 'abs_dens' 1425 write(iunit, '(1x,a)', advance = 'no') label 1426 label = 'rel_dens' 1427 write(iunit, '(1x,a)', advance = 'no') label 1428 label = 'abs_ev' 1429 write(iunit, '(1x,a)', advance = 'no') label 1430 label = 'rel_ev' 1431 write(iunit, '(1x,a)', advance = 'no') label 1432 if (scf%conv_abs_force > M_ZERO) then 1433 label = 'force_diff' 1434 write(iunit, '(1x,a)', advance = 'no') label 1435 end if 1436 if (bitand(ks%xc_family, XC_FAMILY_OEP) /= 0 .and. ks%theory_level /= HARTREE_FOCK) then 1437 if (ks%oep%level == XC_OEP_FULL) then 1438 label = 'OEP norm2ss' 1439 write(iunit, '(1x,a)', advance = 'no') label 1440 end if 1441 end if 1442 write(iunit,'(a)') '' 1443 call io_close(iunit) 1444 end if 1445 1446 end subroutine create_convergence_file 1447 1448 ! ----------------------------------------------------- 1449 1450 subroutine write_convergence_file(dir, fname) 1451 character(len=*), intent(in) :: dir 1452 character(len=*), intent(in) :: fname 1453 1454 integer :: iunit 1455 1456 if(mpi_grp_is_root(mpi_world)) then ! this the absolute master writes 1457 call io_mkdir(dir, namespace) 1458 iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write', position='append') 1459 write(iunit, '(i5,es18.8)', advance = 'no') iter, units_from_atomic(units_out%energy, hm%energy%total) 1460 write(iunit, '(es13.5)', advance = 'no') units_from_atomic(units_out%energy, scf%energy_diff) 1461 write(iunit, '(es13.5)', advance = 'no') scf%abs_dens 1462 write(iunit, '(es13.5)', advance = 'no') scf%rel_dens 1463 write(iunit, '(es13.5)', advance = 'no') units_from_atomic(units_out%energy, scf%abs_ev) 1464 write(iunit, '(es13.5)', advance = 'no') scf%rel_ev 1465 if (scf%conv_abs_force > M_ZERO) then 1466 write(iunit, '(es13.5)', advance = 'no') units_from_atomic(units_out%force, scf%abs_force) 1467 end if 1468 if (bitand(ks%xc_family, XC_FAMILY_OEP) /= 0 .and. ks%theory_level /= HARTREE_FOCK) then 1469 if (ks%oep%level == XC_OEP_FULL) & 1470 write(iunit, '(es13.5)', advance = 'no') ks%oep%norm2ss 1471 end if 1472 write(iunit,'(a)') '' 1473 call io_close(iunit) 1474 end if 1475 1476 end subroutine write_convergence_file 1477 1478 end subroutine scf_run 1479 1480 ! --------------------------------------------------------- 1481 subroutine scf_state_info(st) 1482 class(states_abst_t), intent(in) :: st 1483 1484 PUSH_SUB(scf_state_info) 1485 1486 if (states_are_real(st)) then 1487 call messages_write('Info: SCF using real wavefunctions.') 1488 else 1489 call messages_write('Info: SCF using complex wavefunctions.') 1490 end if 1491 call messages_info() 1492 1493 POP_SUB(scf_state_info) 1494 1495 end subroutine scf_state_info 1496 1497 ! --------------------------------------------------------- 1498 subroutine scf_print_mem_use() 1499 FLOAT :: mem 1500#ifdef HAVE_MPI 1501 FLOAT :: mem_tmp 1502#endif 1503 1504 PUSH_SUB(scf_print_mem_use) 1505 1506 if(conf%report_memory) then 1507 mem = loct_get_memory_usage()/(CNST(1024.0)**2) 1508#ifdef HAVE_MPI 1509 call MPI_Allreduce(mem, mem_tmp, 1, MPI_FLOAT, MPI_SUM, mpi_world%comm, mpi_err) 1510 mem = mem_tmp 1511#endif 1512 write(message(1),'(a,f14.2)') 'Memory usage [Mbytes] :', mem 1513 call messages_info(1) 1514 end if 1515 1516 POP_SUB(scf_print_mem_use) 1517 end subroutine scf_print_mem_use 1518 1519end module scf_oct_m 1520 1521 1522!! Local Variables: 1523!! mode: f90 1524!! coding: utf-8 1525!! End: 1526