1 2! Module for all mixing methods in a standard way 3 4! This module implements mixing of the Pulay and Broyden 5! type. 6 7! The Pulay method is implemented in the fast calculation 8! setup and in the stable method. 9! The stable method is executed if the inversion fails. 10! - Stable: G.Kresse and J.Furthmuller, Comp. Mat. Sci. 6, 15, 1996 11! - gr (guarenteed-reduction) : http://arxiv.org/pdf/cond-mat/0005521.pdf 12 13! All implemented methods employ a restart with variable 14! history saving. 15 16module m_mixing 17 18 use precision, only: dp 19 20#ifdef MPI 21 ! MPI stuff 22 use mpi_siesta 23#endif 24 25 ! Intrinsic classes for retaining history 26 use class_dData1D 27 use class_Fstack_dData1D 28 29 implicit none 30 31 private 32 33 save 34 35 integer, parameter :: MIX_LINEAR = 1 36 integer, parameter :: MIX_PULAY = 2 37 integer, parameter :: MIX_BROYDEN = 3 38 integer, parameter :: MIX_FIRE = 4 39 40 ! Action tokens (binary: 0, 1, 2, 4, 8, ...!) 41 integer, parameter :: ACTION_MIX = 0 42 integer, parameter :: ACTION_RESTART = 1 43 integer, parameter :: ACTION_NEXT = 2 44 45 type tMixer 46 47 ! Name of mixer 48 character(len=24) :: name 49 50 ! The different saved variables per iteration 51 ! and their respective stacks 52 type(Fstack_dData1D), allocatable :: stack(:) 53 54 ! The method of the mixer 55 integer :: m = MIX_PULAY 56 57 ! In case the mixing method has a variant 58 ! this denote the variant 59 ! This value is thus specific for each method 60 integer :: v = 0 61 62 ! The currently reached iteration 63 integer :: cur_itt = 0, start_itt = 0 64 65 ! Different mixers may have different histories 66 integer :: n_hist = 2 67 68 ! Number of iterations using this mixer 69 ! There are a couple of signals here 70 ! == 0 : 71 ! only use this mixer until convergence 72 ! > 0 : 73 ! after having runned n_itt step to "next" 74 integer :: n_itt = 0 75 76 ! When mod(cur_itt,restart_itt) == 0 the history will 77 ! be _reset_ 78 integer :: restart = 0 79 integer :: restart_save = 0 80 81 ! This is an action token specifying the current 82 ! action 83 integer :: action = ACTION_MIX 84 85 ! The next mixing method following this method 86 type(tMixer), pointer :: next => null() 87 88 ! The next mixing method following this method 89 ! Only used if mixing method achieved convergence 90 ! using this method 91 type(tMixer), pointer :: next_conv => null() 92 93 ! ** Parameters specific for the method: 94 95 ! The mixing parameter used for this mixer 96 real(dp) :: w = 0._dp 97 98 ! linear array of real variables used specifically 99 ! for this mixing type 100 real(dp), pointer :: rv(:) => null() 101 integer, pointer :: iv(:) => null() 102 103#ifdef MPI 104 ! In case we have MPI the mixing scheme 105 ! can implement a reduction scheme. 106 ! This can be MPI_Comm_Self to not employ any 107 ! reductions 108 integer :: Comm = MPI_Comm_Self 109#endif 110 111 end type tMixer 112 113 114 ! Indices for special constanst 115 integer, parameter :: I_PREVIOUS_RES = 0 116 integer, parameter :: I_P_RESTART = -1 117 integer, parameter :: I_P_NEXT = -2 118 ! This index should always be the lowest index 119 ! This is used to allocate the correct bounds for the 120 ! additional array of information 121 integer, parameter :: I_SVD_COND = -3 122 123 124 ! Debug mixing runs 125 logical :: debug_mix = .false. 126 ! In case of parallel mixing this also contains the node number 127 character(len=20) :: debug_msg = 'mix:' 128 129 public :: tMixer 130 131 ! Routines are divided in three sections 132 133 ! 1. Routines used to construct the mixers 134 ! Routines used to print information regarding 135 ! the mixers 136 public :: mixers_init, mixer_init 137 public :: mixers_print, mixers_print_block 138 public :: mixers_history_init 139 public :: mixers_reset 140 141 ! 2. Public functions for retrieving information 142 ! from external routines 143 public :: mix_method, mix_method_variant 144 145 ! 3. Actual mixing methods 146 public :: mixing 147 148 public :: MIX_LINEAR, MIX_FIRE, MIX_PULAY, MIX_BROYDEN 149 150 interface mixing 151 module procedure mixing_1d, mixing_2d 152 end interface mixing 153 154contains 155 156 !> Initialize a set of mixers by reading in fdf information. 157 !! @param[in] prefix the fdf-label prefixes 158 !! @param[pointer] mixers the mixers that are to be initialized 159 !! @param[in] Comm @opt optional MPI-communicator 160 subroutine mixers_init( prefix, mixers, Comm ) 161 162 use parallel, only: IONode, Node 163 use fdf 164 165 ! FDF-prefix for searching keywords 166 character(len=*), intent(in) :: prefix 167 ! The array of mixers (has to be nullified upon entry) 168 type(tMixer), pointer :: mixers(:) 169 integer, intent(in), optional :: Comm 170 171 ! Block constructs 172 type(block_fdf) :: bfdf 173 type(parsed_line), pointer :: pline 174 175 ! number of history steps saved 176 integer :: n_hist, n_restart, n_save 177 real(dp) :: w 178 179 integer :: nm, im, im2 180 character(len=10) :: lp 181 character(len=70) :: method, variant 182 183 ! Default mixing options... 184 if ( fdf_get('Mixer.Debug',.false.) ) then 185 debug_mix = IONode 186 debug_msg = 'mix:' 187 end if 188 if ( fdf_get('Mixer.Debug.MPI',.false.) ) then 189 debug_mix = .true. 190 write(debug_msg,'(a,i0,a)') 'mix (',Node,'):' 191 end if 192 193 lp = trim(prefix)//'.Mixer' 194 195 ! ensure nullification 196 call mixers_reset(mixers) 197 198 ! Return immediately if the user hasn't defined 199 ! an fdf-block for the mixing options... 200 if ( .not. fdf_block(trim(lp)//'s', bfdf) ) return 201 202 ! update mixing weight and kick mixing weight 203 w = fdf_get(trim(lp)//'.Weight',0.1_dp) 204 ! Get history length 205 n_hist = fdf_get(trim(lp)//'.History',6) 206 ! Restart after this number of iterations 207 n_restart = fdf_get(trim(lp)//'.Restart',0) 208 n_save = fdf_get(trim(lp)//'.Restart.Save',1) 209 ! negative savings are not allowed 210 n_save = max(0, n_save) 211 212 213 214 ! Read in the options regarding the mixing options 215 nm = 0 216 do while ( fdf_bline(bfdf,pline) ) 217 if ( fdf_bnnames(pline) == 0 ) cycle 218 nm = nm + 1 219 end do 220 if ( nm == 0 ) then 221 call die('mixing: No mixing schemes selected. & 222 &Please at least add one mixer.') 223 end if 224 225 226 ! Allocate all denoted mixers... 227 allocate(mixers(nm)) 228 mixers(:)%w = w 229 mixers(:)%n_hist = n_hist 230 mixers(:)%restart = n_restart 231 mixers(:)%restart_save = n_save 232 233 234 ! Rewind to grab names. 235 call fdf_brewind(bfdf) 236 nm = 0 237 do while ( fdf_bline(bfdf,pline) ) 238 if ( fdf_bnnames(pline) == 0 ) cycle 239 240 nm = nm + 1 241 mixers(nm)%name = fdf_bnames(pline,1) 242 243 end do 244 245 ! Now read all mixers for this segment and their options 246 do im = 1 , nm 247 248 call read_block( mixers(im) ) 249 250 end do 251 252 ! Create history stack and associate correct 253 ! stack pointers 254 call mixers_history_init(mixers) 255 256#ifdef MPI 257 if ( present(Comm) ) then 258 mixers(:)%Comm = Comm 259 else 260 mixers(:)%Comm = MPI_Comm_World 261 end if 262#endif 263 264 contains 265 266 subroutine read_block( m ) 267 type(tMixer), intent(inout), target :: m 268 269 character(len=64) :: opt 270 271 ! create block string 272 opt = trim(lp)//'.'//trim(m%name) 273 274 if ( .not. fdf_block(opt,bfdf) ) then 275 call die('Block: '//trim(opt)//' does not exist!') 276 end if 277 278 ! Default to the pulay method... 279 ! This enables NOT writing this in the block 280 method = 'pulay' 281 variant = ' ' 282 283 ! read method 284 do while ( fdf_bline(bfdf,pline) ) 285 if ( fdf_bnnames(pline) == 0 ) cycle 286 287 opt = fdf_bnames(pline,1) 288 289 if ( leqi(opt,'method') ) then 290 291 method = fdf_bnames(pline,2) 292 293 else if ( leqi(opt,'variant') ) then 294 295 variant = fdf_bnames(pline,2) 296 297 end if 298 299 end do 300 301 ! Retrieve the method and the variant 302 m%m = mix_method(method) 303 m%v = mix_method_variant(m%m, variant) 304 305 ! Define separate defaults which are 306 ! not part of the default input options 307 select case ( m%m ) 308 case ( MIX_LINEAR ) 309 m%n_hist = 0 310 end select 311 312 call fdf_brewind(bfdf) 313 314 ! read options 315 do while ( fdf_bline(bfdf,pline) ) 316 if ( fdf_bnnames(pline) == 0 ) cycle 317 318 opt = fdf_bnames(pline,1) 319 320 if ( leqi(opt,'iterations') & 321 .or. leqi(opt,'itt') ) then 322 323 m%n_itt = fdf_bintegers(pline,1) 324 325 else if ( leqi(opt,'history') ) then 326 327 m%n_hist = fdf_bintegers(pline,1) 328 329 else if ( leqi(opt,'weight') .or. leqi(opt, 'w') ) then 330 331 m%w = fdf_breals(pline,1) 332 333 else if ( leqi(opt,'restart') ) then 334 335 m%restart = fdf_bintegers(pline,1) 336 337 else if ( leqi(opt,'restart.save') ) then 338 339 m%restart_save = fdf_bintegers(pline,1) 340 m%restart_save = max(0,m%restart_save) 341 342 end if 343 344 end do 345 346 ! Initialize the mixer by setting the correct 347 ! standard options and allocate space in the mixers... 348 call mixer_init( m ) 349 350 ! Read the options for this mixer 351 call fdf_brewind(bfdf) 352 353 ! read options 354 do while ( fdf_bline(bfdf,pline) ) 355 if ( fdf_bnnames(pline) == 0 ) cycle 356 357 opt = fdf_bnames(pline,1) 358 359 if ( leqi(opt,'next') ) then 360 361 nullify(m%next) 362 363 opt = fdf_bnames(pline,2) 364 do im2 = 1 , nm 365 if ( leqi(opt,mixers(im2)%name) ) then 366 m%next => mixers(im2) 367 exit 368 end if 369 end do 370 371 if ( .not. associated(m%next) ) then 372 call die('mixing: Could not find next mixer. & 373 &Ensure all mixers exist and their names.') 374 end if 375 376 if ( associated(m%next, target=m) ) then 377 call die('mixing: Next *must* not be it-self. & 378 &Please change accordingly.') 379 end if 380 381 else if ( leqi(opt,'next.conv') ) then 382 383 nullify(m%next_conv) 384 385 opt = fdf_bnames(pline,2) 386 do im2 = 1 , nm 387 if ( leqi(opt,mixers(im2)%name) ) then 388 m%next_conv => mixers(im2) 389 exit 390 end if 391 end do 392 393 if ( .not. associated(m%next_conv) ) then 394 call die('mixing: Could not find next convergence mixer. & 395 &Ensure all mixers exist and their names.') 396 end if 397 398 if ( associated(m%next_conv,target=m) ) then 399 call die('mixing: next.conv *must* not be it-self. & 400 &Please change accordingly.') 401 end if 402 403 end if 404 405 end do 406 407 ! Ensure that if a next have not been specified 408 ! it will continue indefinitely. 409 if ( .not. associated(m%next) ) then 410 m%n_itt = 0 411 end if 412 413 414 ! Read the options for this mixer 415 call fdf_brewind(bfdf) 416 417 ! read options 418 do while ( fdf_bline(bfdf,pline) ) 419 ! skip lines without associated content 420 if ( fdf_bnnames(pline) == 0 ) cycle 421 422 opt = fdf_bnames(pline,1) 423 424 ! Do options so that a pulay option may refer to 425 ! the actual names of the constants 426 427 if ( m%m == MIX_PULAY ) then 428 if ( leqi(opt,'weight.linear') & 429 .or. leqi(opt,'w.linear') ) then 430 431 m%rv(1) = fdf_breals(pline,1) 432 if ( m%rv(1) <= 0._dp .or. 1._dp < m%rv(1) ) then 433 call die("m_mixing: Mixing weight should be 0 < weight <= 1") 434 end if 435 436 else if ( leqi(opt,'svd.cond') ) then 437 m%rv(I_SVD_COND) = fdf_bvalues(pline,1) 438 end if 439 440 else if ( m%m == MIX_BROYDEN ) then 441 442 ! The linear mixing weight 443 if ( leqi(opt,'weight.linear') & 444 .or. leqi(opt,'w.linear') ) then 445 446 m%rv(1) = fdf_breals(pline,1) 447 if ( m%rv(1) <= 0._dp .or. 1._dp < m%rv(1) ) then 448 call die("m_mixing: Linear mixing weight should be 0 < weight <= 1") 449 end if 450 451 else if ( leqi(opt,'weight.prime') & 452 .or. leqi(opt,'w.prime') ) then 453 454 m%rv(2) = fdf_breals(pline,1) 455 if ( m%rv(2) < 0._dp .or. 1._dp < m%rv(2) ) then 456 call die("m_mixing: Prime weight should be 0 <= weight <= 1") 457 end if 458 459 else if ( leqi(opt,'svd.cond') ) then 460 461 m%rv(I_SVD_COND) = fdf_bvalues(pline,1) 462 463 end if 464 465 end if 466 467 ! Generic options for all advanced methods... 468 if ( leqi(opt,'next.p') ) then 469 470 ! Only allow stepping to the next when 471 ! having a next associated 472 if ( associated(m%next) ) then 473 m%rv(I_P_NEXT) = fdf_bvalues(pline,1) 474 end if 475 476 else if ( leqi(opt,'restart.p') ) then 477 478 m%rv(I_P_RESTART) = fdf_bvalues(pline,1) 479 480 end if 481 482 end do 483 484 end subroutine read_block 485 486 end subroutine mixers_init 487 488 489 !> Initialize a single mixer depending on the preset 490 !! options. Useful for external correct setup. 491 !! 492 !! @param[inout] mix mixer to be initialized 493 subroutine mixer_init( mix ) 494 type(tMixer), intent(inout) :: mix 495 integer :: n 496 497 if ( mix%w <= 0._dp .or. 1._dp < mix%w ) then 498 call die("m_mixing: Mixing weight should be: 0 < weight <= 1") 499 end if 500 501 ! Correct amount of history in the mixing. 502 if ( 0 < mix%restart .and. & 503 mix%restart < mix%n_hist ) then 504 ! This is if we restart this scheme, 505 ! then it does not make sense to have a history 506 ! greater than the restart count 507 mix%n_hist = mix%restart 508 end if 509 if ( 0 < mix%n_itt .and. & 510 mix%n_itt < mix%n_hist ) then 511 ! If this only runs for n_itt itterations, 512 ! it makes no sense to have a history greater 513 ! than this. 514 mix%n_hist = mix%n_itt 515 end if 516 517 select case ( mix%m ) 518 case ( MIX_LINEAR ) 519 520 allocate(mix%rv(I_SVD_COND:0)) 521 ! Kill any history settings that do not apply to the 522 ! linear mixer. 523 mix%restart = 0 524 mix%restart_save = 0 525 526 case ( MIX_PULAY ) 527 528 allocate(mix%rv(I_SVD_COND:1)) 529 mix%rv(1) = mix%w 530 ! We allocate the double residual (n_hist-1) 531 mix%n_hist = max(2, mix%n_hist) 532 if ( mix%v == 1 .or. mix%v == 3 ) then 533 534 ! The GR method requires an even number 535 ! of restart steps 536 ! And then we ensure the history to be aligned 537 ! with a restart (restart has precedence) 538 if ( mix%restart /= 0 ) then 539 mix%restart = mix%restart + mod(mix%restart, 2) 540 end if 541 542 end if 543 544 case ( MIX_BROYDEN ) 545 546 ! allocate temporary array 547 mix%n_hist = max(2, mix%n_hist) 548 n = 2 + mix%n_hist 549 allocate(mix%rv(I_SVD_COND:n)) 550 mix%rv(1:n) = mix%w 551 552 end select 553 554 if ( mix%restart < 0 ) then 555 call die('mixing: restart count must be positive') 556 end if 557 558 mix%restart_save = min(mix%n_hist - 1, mix%restart_save) 559 mix%restart_save = max(0, mix%restart_save) 560 561 ! This is the restart parameter 562 ! I.e. if |f_k / f - 1| < rp 563 ! only works for positive rp 564 mix%rv(I_PREVIOUS_RES) = huge(1._dp) 565 mix%rv(I_P_RESTART) = -1._dp 566 mix%rv(I_P_NEXT) = -1._dp 567 mix%rv(I_SVD_COND) = 1.e-8_dp 568 569 end subroutine mixer_init 570 571 572 !> Initialize all history for the mixers 573 !! 574 !! Routine for clearing all history and setting up the 575 !! arrays so that they may be used subsequently. 576 !! 577 !! @param[inout] mixers the mixers to be initialized 578 subroutine mixers_history_init( mixers ) 579 type(tMixer), intent(inout), target :: mixers(:) 580 581 type(tMixer), pointer :: m 582 integer :: im, is, ns 583 logical :: is_GR 584 585 do im = 1 , size(mixers) 586 m => mixers(im) 587 if ( debug_mix .and. current_itt(m) >= 1 ) then 588 write(*,'(a,a)') trim(debug_msg), & 589 ' resetting history of all mixers' 590 exit 591 end if 592 end do 593 594 ! Clean up all arrays and reference counted 595 ! objects 596 do im = 1 , size(mixers) 597 598 m => mixers(im) 599 600 ! reset history track 601 m%start_itt = 0 602 m%cur_itt = 0 603 604 ! do not try and de-allocate something not 605 ! allocated 606 if ( allocated(m%stack) ) then 607 608 ns = size(m%stack) 609 do is = 1 , ns 610 call delete(m%stack(is)) 611 end do 612 613 ! clean-up 614 deallocate(m%stack) 615 616 end if 617 618 ! Re-populate 619 select case ( m%m ) 620 case ( MIX_LINEAR ) 621 622 ! do nothing 623 624 case ( MIX_PULAY ) 625 626 is_GR = (m%v == 1) .or. (m%v == 3) 627 628 if ( .not. is_GR ) then 629 allocate(m%stack(3)) 630 else 631 allocate(m%stack(2)) 632 end if 633 634 ! These arrays contains these informations 635 ! s1 = m%stack(1) 636 ! s2 = m%stack(2) 637 ! s3 = m%stack(3) 638 ! Here <> is input function, x[in], and 639 ! <>' is the corresponding output, x[out]. 640 ! First iteration: 641 ! s1 = { 1' - 1 } 642 ! s3 = { 1' } 643 ! Second iteration 644 ! s2 = { 2' - 2 - (1' - 1) } 645 ! s1 = { 2 - 1 , 2' - 2 } 646 ! s3 = { 2' } 647 ! Third iteration 648 ! s2 = { 2' - 2 - (1' - 1) , 3' - 3 - (2' - 2) } 649 ! s1 = { 2 - 1 , 3 - 2, 3' - 3 } 650 ! s3 = { 3' } 651 ! and so on 652 653 ! allocate x[i+1] - x[i] 654 call new(m%stack(1), m%n_hist) 655 ! allocate F[i+1] - F[i] 656 call new(m%stack(2), m%n_hist-1) 657 658 if ( .not. is_GR ) then 659 call new(m%stack(3), 1) 660 end if 661 662 case ( MIX_BROYDEN ) 663 664 ! Same as original Pulay 665 allocate(m%stack(3)) 666 call new(m%stack(1), m%n_hist) 667 call new(m%stack(2), m%n_hist-1) 668 call new(m%stack(3), 1) 669 670 end select 671 672 end do 673 674 end subroutine mixers_history_init 675 676 677 !> Reset the mixers, i.e. clean _everything_ 678 !! 679 !! Also deallocates (and nullifies) the input array! 680 !! 681 !! @param[inout] mixers array of mixers to be cleaned 682 subroutine mixers_reset( mixers ) 683 type(tMixer), pointer :: mixers(:) 684 685 type(tMixer), pointer :: m 686 687 integer :: im, is, ns 688 689 if ( .not. associated(mixers) ) return 690 691 do im = 1 , size(mixers) 692 m => mixers(im) 693 694 if ( allocated(m%stack) ) then 695 ns = size(m%stack) 696 do is = 1 , ns 697 call delete(m%stack(is)) 698 end do 699 deallocate(m%stack) 700 end if 701 702 if ( associated(m%rv) ) then 703 deallocate(m%rv) 704 nullify(m%rv) 705 end if 706 707 if ( associated(m%iv) ) then 708 deallocate(m%iv) 709 nullify(m%iv) 710 end if 711 712 end do 713 714 deallocate(mixers) 715 nullify(mixers) 716 717 end subroutine mixers_reset 718 719 720 !> Print (to std-out) information regarding the mixers 721 !! 722 !! @param[in] prefix the prefix (fdf) for the mixers 723 !! @param[in] mixers array of mixers allocated 724 subroutine mixers_print( prefix, mixers ) 725 726 use parallel, only: IONode 727 728 character(len=*), intent(in) :: prefix 729 type(tMixer), intent(in), target :: mixers(:) 730 731 type(tMixer), pointer :: m 732 character(len=64) :: fmt 733 734 logical :: bool 735 integer :: i 736 737 if ( .not. IONode ) return 738 739 fmt = 'mix.'//trim(prefix)//':' 740 741 if ( debug_mix ) then 742 write(*,'(2a,t50,''= '',l)') trim(fmt), & 743 ' Debug messages',debug_mix 744 end if 745 746 ! Print out options for all mixers 747 do i = 1 , size(mixers) 748 749 m => mixers(i) 750 751 select case ( m%m ) 752 case ( MIX_LINEAR ) 753 754 write(*,'(2a,t50,''= '',a)') trim(fmt), & 755 ' Linear mixing',trim(m%name) 756 write(*,'(2a,t50,''= '',f12.6)') trim(fmt), & 757 ' Mixing weight',m%w 758 759 if ( m%n_hist > 0 .and. (& 760 associated(m%next) & 761 .or. associated(m%next_conv)) ) then 762 write(*,'(2a,t50,''= '',i0)') trim(fmt), & 763 ' Carried history steps',m%n_hist 764 end if 765 766 case ( MIX_PULAY ) 767 768 write(*,'(2a,t50,''= '',a)') trim(fmt), & 769 ' Pulay mixing',trim(m%name) 770 771 select case ( m%v ) 772 case ( 0 ) 773 write(*,'(2a,t50,''= '',a)') trim(fmt), & 774 ' Variant','stable' 775 case ( 1 ) 776 write(*,'(2a,t50,''= '',a)') trim(fmt), & 777 ' Variant','GR' 778 case ( 2 ) 779 write(*,'(2a,t50,''= '',a)') trim(fmt), & 780 ' Variant','stable-SVD' 781 case ( 3 ) 782 write(*,'(2a,t50,''= '',a)') trim(fmt), & 783 ' Variant','GR-SVD' 784 end select 785 786 write(*,'(2a,t50,''= '',i0)') trim(fmt), & 787 ' History steps',m%n_hist 788 write(*,'(2a,t50,''= '',f12.6)') trim(fmt), & 789 ' Linear mixing weight',m%rv(1) 790 write(*,'(2a,t50,''= '',f12.6)') trim(fmt), & 791 ' Mixing weight',m%w 792 write(*,'(2a,t50,''= '',e10.4)') trim(fmt), & 793 ' SVD condition',m%rv(I_SVD_COND) 794 if ( m%rv(I_P_NEXT) > 0._dp ) then 795 write(*,'(2a,t50,''= '',f6.4)') trim(fmt), & 796 ' Step mixer parameter',m%rv(I_P_NEXT) 797 end if 798 bool = .false. 799 if ( m%rv(I_P_RESTART) > 0._dp ) then 800 write(*,'(2a,t50,''= '',f6.4)') trim(fmt), & 801 ' Restart parameter',m%rv(I_P_RESTART) 802 bool = .true. 803 end if 804 if ( m%restart > 0 ) then 805 write(*,'(2a,t50,''= '',i0)') trim(fmt), & 806 ' Restart steps',m%restart 807 bool = .true. 808 end if 809 if ( bool ) then 810 write(*,'(2a,t50,''= '',i0)') trim(fmt), & 811 ' Restart save steps',m%restart_save 812 end if 813 814 case ( MIX_BROYDEN ) 815 816 write(*,'(2a,t50,''= '',a)') trim(fmt), & 817 ' Broyden mixing',trim(m%name) 818 819 !write(*,'(2a,t50,''= '',a)') trim(fmt), & 820 ! ' Variant','original' 821 822 write(*,'(2a,t50,''= '',i0)') trim(fmt), & 823 ' History steps',m%n_hist 824 write(*,'(2a,t50,''= '',f12.6)') trim(fmt), & 825 ' Linear mixing weight',m%rv(1) 826 write(*,'(2a,t50,''= '',f12.6)') trim(fmt), & 827 ' Inverse Jacobian weight',m%w 828 write(*,'(2a,t50,''= '',f12.6)') trim(fmt), & 829 ' Weight prime',m%rv(2) 830 write(*,'(2a,t50,''= '',e10.4)') trim(fmt), & 831 ' SVD condition',m%rv(I_SVD_COND) 832 if ( m%rv(I_P_NEXT) > 0._dp ) then 833 write(*,'(2a,t50,''= '',f6.4)') trim(fmt), & 834 ' Step mixer parameter',m%rv(I_P_NEXT) 835 end if 836 bool = .false. 837 if ( m%rv(I_P_RESTART) > 0._dp ) then 838 write(*,'(2a,t50,''= '',f6.4)') trim(fmt), & 839 ' Restart parameter',m%rv(I_P_RESTART) 840 bool = .true. 841 end if 842 if ( m%restart > 0 ) then 843 write(*,'(2a,t50,''= '',i0)') trim(fmt), & 844 ' Restart steps',m%restart 845 bool = .true. 846 end if 847 if ( bool ) then 848 write(*,'(2a,t50,''= '',i0)') trim(fmt), & 849 ' Restart save steps',m%restart_save 850 end if 851 852 case ( MIX_FIRE ) 853 854 write(*,'(2a,t50,''= '',a)') trim(fmt), & 855 ' Fire mixing',trim(m%name) 856 857 end select 858 859 if ( m%n_itt > 0 ) then 860 write(*,'(2a,t50,''= '',i0)') trim(fmt), & 861 ' Number of mixing iterations',m%n_itt 862 if ( associated(m%next) ) then 863 write(*,'(2a,t50,''= '',a)') trim(fmt), & 864 ' Following mixer',trim(m%next%name) 865 else 866 call die('Something went wrong, if the mixer does not go & 867 &indefinitely it should have a following method.') 868 end if 869 end if 870 871 if ( associated(m%next_conv) ) then 872 write(*,'(2a,t50,''= '',a)') trim(fmt), & 873 ' Following mixer upon convergence',trim(m%next_conv%name) 874 end if 875 876 end do 877 878 end subroutine mixers_print 879 880 881 !> Print (to std-out) the fdf-blocks that recreate the mixer settings 882 !! 883 !! @param[in] prefix the fdf-prefix for reading the blocks 884 !! @param[in] mixers array of mixers that should be printed 885 !! their fdf-blocks 886 subroutine mixers_print_block( prefix, mixers ) 887 use parallel, only: IONode 888 889 character(len=*), intent(in) :: prefix 890 type(tMixer), intent(in), target :: mixers(:) 891 892 type(tMixer), pointer :: m 893 894 logical :: bool 895 integer :: i 896 897 if ( .not. IONode ) return 898 899 ! Write block of input 900 write(*,'(/3a)')'%block ',trim(prefix), '.Mixers' 901 do i = 1 , size(mixers) 902 m => mixers(i) 903 write(*,'(t3,a)') trim(m%name) 904 end do 905 write(*,'(3a)')'%endblock ',trim(prefix), '.Mixers' 906 907 908 ! Print out options for all mixers 909 do i = 1 , size(mixers) 910 911 m => mixers(i) 912 913 ! Write out this block 914 write(*,'(/4a)')'%block ',trim(prefix),'.Mixer.',trim(m%name) 915 916 write(*,'(t3,a)') '# Mixing method' 917 918 ! Write out method 919 select case ( m%m ) 920 case ( MIX_LINEAR ) 921 922 write(*,'(t2,2(tr1,a))') 'method','linear' 923 924 case ( MIX_PULAY ) 925 926 write(*,'(t2,2(tr1,a))') 'method','pulay' 927 select case ( m%v ) 928 case ( 0 ) 929 write(*,'(t2,2(tr1,a))') 'variant','stable' 930 case ( 1 ) 931 write(*,'(t2,2(tr1,a))') 'variant','GR' 932 case ( 2 ) 933 write(*,'(t2,2(tr1,a))') 'variant','stable+SVD' 934 case ( 3 ) 935 write(*,'(t2,2(tr1,a))') 'variant','GR+SVD' 936 end select 937 938 case ( MIX_BROYDEN ) 939 940 write(*,'(t2,2(tr1,a))') 'method','broyden' 941 942 ! currently no variants exists 943 944 end select 945 946 947 948 ! remark 949 write(*,'(/,t3,a)') '# Mixing options' 950 951 ! Weight 952 ! For Broyden this is the inverse Jacobian 953 write(*,'(t3,a,f6.4)') 'weight ', m%w 954 select case ( m%m ) 955 case ( MIX_PULAY ) 956 write(*,'(t3,a,f6.4)') 'weight.linear ', m%rv(1) 957 case ( MIX_BROYDEN ) 958 write(*,'(t3,a,f6.4)') 'weight.linear ', m%rv(1) 959 write(*,'(t3,a,f6.4)') 'weight.prime ', m%rv(2) 960 end select 961 962 if ( m%n_hist > 0 ) then 963 write(*,'(t3,a,i0)') 'history ', m%n_hist 964 end if 965 966 bool = .false. 967 if ( m%restart > 0 ) then 968 write(*,'(t3,a,i0)') 'restart ', m%restart 969 bool = .true. 970 end if 971 select case ( m%m ) 972 case ( MIX_PULAY, MIX_BROYDEN ) 973 if ( m%rv(I_P_RESTART) > 0._dp ) then 974 write(*,'(t3,a,e10.5)')'restart.p ', m%rv(I_P_RESTART) 975 bool = .true. 976 end if 977 end select 978 if ( bool ) then 979 write(*,'(t3,a,i0)')'restart.save ', m%restart_save 980 end if 981 982 983 984 ! remark 985 986 bool = .false. 987 if ( m%n_itt > 0 ) then 988 write(*,'(/,t3,a)') '# Continuation options' 989 write(*,'(t3,a,i0)')'iterations ', m%n_itt 990 bool = .true. 991 end if 992 select case ( m%m ) 993 case ( MIX_PULAY , MIX_BROYDEN ) 994 if ( m%rv(I_P_NEXT) > 0._dp ) then 995 if ( .not. bool ) & 996 write(*,'(/,t3,a)') '# Continuation options' 997 write(*,'(t3,a,f6.4)')'next.p ', m%rv(I_P_NEXT) 998 bool = .true. 999 end if 1000 end select 1001 if ( bool .and. associated(m%next) ) then 1002 write(*,'(t2,2(tr1,a))') 'next', trim(m%next%name) 1003 else if ( bool ) then 1004 call die('Something went wrong, if the mixer does not go & 1005 &indefinitely it should have a following method.') 1006 end if 1007 1008 if ( associated(m%next_conv) ) then 1009 if ( .not. bool ) & 1010 write(*,'(/,t3,a)') '# Continuation options' 1011 write(*,'(t2,2(tr1,a))') 'next.conv', trim(m%next_conv%name) 1012 end if 1013 1014 1015 write(*,'(4a)')'%endblock ',trim(prefix),'.Mixer.',trim(m%name) 1016 1017 end do 1018 1019 write(*,*) ! new-line 1020 1021 end subroutine mixers_print_block 1022 1023 1024 !> Return the integer specification of the mixing type 1025 !! 1026 !! @param[in] str the character representation of the mixing type 1027 !! @return the integer corresponding to the mixing type 1028 function mix_method(str) result(m) 1029 use fdf, only: leqi 1030 character(len=*), intent(in) :: str 1031 integer :: m 1032 1033 if ( leqi(str,'linear') ) then 1034 m = MIX_LINEAR 1035 else if ( leqi(str,'pulay') .or. & 1036 leqi(str,'diis') .or. & 1037 leqi(str, 'anderson') ) then 1038 m = MIX_PULAY 1039 else if ( leqi(str,'broyden') ) then 1040 m = MIX_BROYDEN 1041 else if ( leqi(str,'fire') ) then 1042 m = MIX_FIRE 1043 call die('mixing: FIRE currently not supported.') 1044 else 1045 call die('mixing: Unknown mixing variant.') 1046 end if 1047 1048 end function mix_method 1049 1050 1051 !> Return the variant of the mixing method 1052 !! 1053 !! @param[in] m the integer type of the mixing method 1054 !! @param[in] str the character specification of the mixing method variant 1055 !! @return the variant of the mixing method 1056 function mix_method_variant(m, str) result(v) 1057 use fdf, only: leqi 1058 integer, intent(in) :: m 1059 character(len=*), intent(in) :: str 1060 integer :: v 1061 1062 v = 0 1063 select case ( m ) 1064 case ( MIX_LINEAR ) 1065 ! no variants 1066 case ( MIX_PULAY ) 1067 v = 0 1068 ! We do not implement tho non-stable version 1069 ! There is no need to have an inferior Pulay mixer... 1070 if ( leqi(str,'original') .or. & 1071 leqi(str,'kresse') .or. leqi(str,'stable') ) then 1072 ! stable version, will nearly always succeed on inversion 1073 v = 0 1074 else if ( leqi(str,'original+svd') .or. & 1075 leqi(str,'kresse+svd') .or. leqi(str,'stable+svd') ) then 1076 ! stable version, will nearly always succeed on inversion 1077 v = 2 1078 else if ( leqi(str,'gr') .or. & 1079 leqi(str,'guarenteed-reduction') .or. & 1080 leqi(str,'bowler-gillan') ) then 1081 ! Guarenteed reduction version 1082 v = 1 1083 else if ( leqi(str,'gr+svd') .or. & 1084 leqi(str,'guarenteed-reduction+svd') .or. & 1085 leqi(str,'bowler-gillan+svd') ) then 1086 ! Guarenteed reduction version 1087 v = 3 1088 end if 1089 case ( MIX_BROYDEN ) 1090 ! Currently only one variant 1091 v = 0 1092 case ( MIX_FIRE ) 1093 ! no variants 1094 end select 1095 1096 end function mix_method_variant 1097 1098 1099 1100 1101 ! The basic mixing procedure is this: 1102 ! 1. Initialize the mixing algorithm 1103 ! This will typically mean that one needs to 1104 ! push input and output matrices 1105 ! 2. Calculate the mixing coefficients 1106 ! 3. Use coefficients to calculate the 1107 ! next (optimized) guess 1108 ! 4. Finalize the mixing method 1109 ! 1110 ! Having the routines split up in this manner 1111 ! allows one to skip step 2 and use coefficients 1112 ! from another set of input/output to retrieve the 1113 ! mixing coefficients. 1114 ! Say we may retrieve mixing coefficients from 1115 ! the Hamiltonian, but use them for the density-matrix 1116 1117 !> Initialize the mixing algorithm 1118 !! 1119 !! @param[pointer] mix the mixing method 1120 !! @param[in] n size of the arrays to be used in the algorithm 1121 !! @param[in] xin array of the input variables 1122 !! @param[in] xout array of the output variables 1123 subroutine mixing_init( mix, n, xin, F) 1124 1125 ! The current mixing method 1126 type(tMixer), pointer :: mix 1127 integer, intent(in) :: n 1128 ! In/out of the function 1129 real(dp), intent(in) :: xin(n), F(n) 1130 1131 real(dp), pointer :: res(:), rres(:) 1132 1133 integer :: i, ns 1134 real(dp) :: dnorm, dtmp 1135 logical :: p_next, p_restart 1136 1137 ! Initialize action for mixer 1138 mix%action = ACTION_MIX 1139 1140 ! Step iterator (so first mixing has cur_itt == 1) 1141 mix%cur_itt = mix%cur_itt + 1 1142 1143 ! If we are going to skip to next, we signal it 1144 ! before entering 1145 if ( mix%n_itt > 0 .and. & 1146 mix%n_itt <= current_itt(mix) ) then 1147 mix%action = IOR(mix%action, ACTION_NEXT) 1148 end if 1149 1150 1151 ! Check whether the residual norm is below a certain 1152 ! criteria 1153 p_next = mix%rv(I_P_NEXT) > 0._dp 1154 p_restart = mix%rv(I_P_RESTART) > 0._dp 1155 1156 ! Check whether a parameter next/restart is required 1157 if ( p_restart .or. p_next ) then 1158 1159 ! Calculate norm: ||f_k|| 1160 dnorm = norm(n, F, F) 1161#ifdef MPI 1162 dtmp = dnorm 1163 call MPI_AllReduce(dtmp,dnorm,1, & 1164 MPI_double_precision, MPI_Sum, & 1165 mix%Comm, i) 1166#endif 1167 1168 ! Calculate the relative difference 1169 dtmp = abs(dnorm / mix%rv(I_PREVIOUS_RES) - 1._dp) 1170 1171 ! We first check for next, that has precedence 1172 if ( p_next ) then 1173 1174 if ( dtmp < mix%rv(I_P_NEXT) ) then 1175 ! Signal stepping mixer 1176 mix%action = IOR(mix%action, ACTION_NEXT) 1177 end if 1178 1179 if ( debug_mix .and. current_itt(mix) > 1 ) & 1180 write(*,'(a,2(a,e8.3))') trim(debug_msg), & 1181 ' | ||f_k|| - ||f_k-1|| |/||f_k-1|| < np : ', & 1182 dtmp, ' < ', mix%rv(I_P_NEXT) 1183 1184 end if 1185 1186 if ( p_restart ) then 1187 1188 if ( dtmp < mix%rv(I_P_RESTART) ) then 1189 ! Signal restart 1190 mix%action = IOR(mix%action, ACTION_RESTART) 1191 end if 1192 1193 if ( debug_mix .and. current_itt(mix) > 1 ) & 1194 write(*,'(a,2(a,e8.3))') trim(debug_msg), & 1195 ' | ||f_k|| - ||f_k-1|| |/||f_k-1|| < rp : ', & 1196 dtmp, ' < ', mix%rv(I_P_RESTART) 1197 1198 end if 1199 1200 ! Store the new residual norm 1201 mix%rv(I_PREVIOUS_RES) = dnorm 1202 1203 end if 1204 1205 1206 ! Push information to the stack 1207 select case ( mix%m ) 1208 case ( MIX_LINEAR ) 1209 if ( debug_mix ) & 1210 write(*,'(2a)') trim(debug_msg),' linear' 1211 call init_linear() 1212 case ( MIX_PULAY ) 1213 if ( debug_mix ) then 1214 select case ( mix%v ) 1215 case ( 0 ) 1216 write(*,'(2a)') trim(debug_msg),' Pulay' 1217 case ( 1 ) 1218 write(*,'(2a)') trim(debug_msg),' Pulay, GR' 1219 case ( 2 ) 1220 write(*,'(2a)') trim(debug_msg),' Pulay-SVD' 1221 case ( 3 ) 1222 write(*,'(2a)') trim(debug_msg),' Pulay-SVD, GR' 1223 end select 1224 end if 1225 call init_pulay() 1226 case ( MIX_BROYDEN ) 1227 if ( debug_mix ) & 1228 write(*,'(2a)') trim(debug_msg),' Broyden' 1229 call init_broyden() 1230 end select 1231 1232 contains 1233 1234 subroutine init_linear() 1235 1236 ! information for this depends on the 1237 ! following method 1238 call fake_history_from_linear(mix%next) 1239 call fake_history_from_linear(mix%next_conv) 1240 1241 end subroutine init_linear 1242 1243 subroutine init_pulay() 1244 logical :: GR_linear 1245 1246 select case ( mix%v ) 1247 case ( 0 , 2 ) ! Stable Pulay 1248 1249 ! Add the residual to the stack 1250 call push_F(mix%stack(1), n, F) 1251 1252 ns = n_items(mix%stack(1)) 1253 1254 ! Add the residuals of the residuals if applicable 1255 if ( ns > 1 ) then 1256 1257 ! Create F[i+1] - F[i] 1258 call push_diff(mix%stack(2), mix%stack(1)) 1259 1260 ! Update the residual to reflect the input residual 1261 res => getstackval(mix, 1, ns-1) 1262 rres => getstackval(mix, 3) 1263 1264!$OMP parallel do default(shared), private(i) 1265 do i = 1 , n 1266 res(i) = res(i) - rres(i) + xin(i) 1267 end do 1268!$OMP end parallel do 1269 1270 end if 1271 1272 case ( 1 , 3 ) 1273 1274 ! Whether this is the linear cycle... 1275 GR_linear = mod(current_itt(mix), 2) == 1 1276 1277 ! Add the residual to the stack 1278 call push_F(mix%stack(1), n, F, mix%rv(1)) 1279 1280 ns = n_items(mix%stack(1)) 1281 1282 if ( GR_linear .and. current_itt(mix) > 1 .and. & 1283 ns > 1 ) then 1284 1285 res => getstackval(mix, 1) 1286 rres => getstackval(mix, 2) 1287!$OMP parallel do default(shared), private(i) 1288 do i = 1 , n 1289 rres(i) = rres(i) + res(i) 1290 end do 1291!$OMP end parallel do 1292 1293 else if ( ns > 1 .and. .not. GR_linear ) then 1294 1295 ! now we can calculate RRes[i] 1296 call push_diff(mix%stack(2),mix%stack(1)) 1297 1298 end if 1299 1300 end select 1301 1302 end subroutine init_pulay 1303 1304 subroutine init_broyden() 1305 1306 ! Add the residual to the stack 1307 call push_F(mix%stack(1), n, F) 1308 1309 ns = n_items(mix%stack(1)) 1310 1311 ! Add the residuals of the residuals if applicable 1312 if ( ns > 1 ) then 1313 1314 ! Create F[i+1] - F[i] 1315 call push_diff(mix%stack(2), mix%stack(1)) 1316 1317 ! Update the residual to reflect the input residual 1318 res => getstackval(mix, 1, ns-1) 1319 rres => getstackval(mix, 3) 1320 1321!$OMP parallel do default(shared), private(i) 1322 do i = 1 , n 1323 res(i) = res(i) - rres(i) + xin(i) 1324 end do 1325!$OMP end parallel do 1326 1327 else 1328 1329 ! Store F[x_in] (used to create the input residual) 1330 call push_stack_data(mix%stack(3), n) 1331 res => getstackval(mix, 3) 1332!$OMP parallel do default(shared), private(i) 1333 do i = 1 , n 1334 res(i) = xin(i) + F(i) 1335 end do 1336!$OMP end parallel do 1337 1338 end if 1339 1340 end subroutine init_broyden 1341 1342 subroutine fake_history_from_linear(next) 1343 type(tMixer), pointer :: next 1344 real(dp), pointer :: t1(:), t2(:) 1345 integer :: ns, nh, i, nhl 1346 1347 if ( .not. associated(next) ) return 1348 1349 ! Reduce to # history of linear 1350 nhl = mix%n_hist 1351 ! if the number of fake-history steps saved is 1352 ! zero we immediately return. 1353 ! Only if mix%n_hist > 0 will the below 1354 ! occur. 1355 if ( nhl == 0 ) return 1356 1357 ! Check for the type of following method 1358 select case ( next%m ) 1359 case ( MIX_PULAY ) 1360 1361 ! Here it depends on the variant 1362 select case ( next%v ) 1363 case ( 0 , 2 ) ! stable pulay mixing 1364 1365 ! Add the residual to the stack 1366 call push_F(next%stack(1), n, F) 1367 1368 ns = n_items(next%stack(1)) 1369 1370 ! Add the residuals of the residuals if applicable 1371 if ( ns > 1 ) then 1372 1373 ! Create F[i+1] - F[i] 1374 call push_diff(next%stack(2), next%stack(1)) 1375 1376 ! Update the residual to reflect the input residual 1377 t1 => getstackval(next, 1, ns-1) 1378 t2 => getstackval(next, 3) 1379 1380!$OMP parallel do default(shared), private(i) 1381 do i = 1 , n 1382 t1(i) = t1(i) - t2(i) + xin(i) 1383 t2(i) = xin(i) + F(i) 1384 end do 1385!$OMP end parallel do 1386 1387 else 1388 1389 call push_stack_data(next%stack(3), n) 1390 t1 => getstackval(next, 3) 1391!$OMP parallel do default(shared), private(i) 1392 do i = 1 , n 1393 t1(i) = xin(i) + F(i) 1394 end do 1395!$OMP end parallel do 1396 1397 end if 1398 1399 ! Clean up the data... 1400 if ( ns >= nhl ) then 1401 call reset(next%stack(1), 1) 1402 call reset(next%stack(2), 1) 1403 ns = ns - 1 1404 end if 1405 1406 nh = max_size(next%stack(1)) 1407 1408 if ( debug_mix ) & 1409 write(*,'(a,2(a,i0))') trim(debug_msg), & 1410 ' next%n_hist = ', ns, ' / ', nh 1411 1412 end select 1413 1414 case ( MIX_BROYDEN ) 1415 1416 ! Add the residual to the stack 1417 call push_F(next%stack(1), n, F) 1418 1419 ns = n_items(next%stack(1)) 1420 1421 ! Add the residuals of the residuals if applicable 1422 if ( ns >= 2 ) then 1423 1424 ! Create F[i+1] - F[i] 1425 call push_diff(next%stack(2), next%stack(1)) 1426 1427 ! Update the residual to reflect the input residual 1428 t1 => getstackval(next, 1, ns-1) 1429 t2 => getstackval(next, 3) 1430 1431!$OMP parallel do default(shared), private(i) 1432 do i = 1 , n 1433 t1(i) = t1(i) - t2(i) + xin(i) 1434 t2(i) = xin(i) + F(i) 1435 end do 1436!$OMP end parallel do 1437 1438 else 1439 1440 call push_stack_data(next%stack(3), n) 1441 t1 => getstackval(next, 3) 1442!$OMP parallel do default(shared), private(i) 1443 do i = 1 , n 1444 t1(i) = xin(i) + F(i) 1445 end do 1446!$OMP end parallel do 1447 1448 end if 1449 1450 ! Clean up the data... 1451 if ( ns >= nhl ) then 1452 call reset(next%stack(1), 1) 1453 call reset(next%stack(2), 1) 1454 ns = ns - 1 1455 end if 1456 1457 nh = max_size(next%stack(1)) 1458 1459 if ( debug_mix ) & 1460 write(*,'(a,2(a,i0))') trim(debug_msg), & 1461 ' next%n_hist = ', ns, ' / ', nh 1462 1463 end select 1464 1465 end subroutine fake_history_from_linear 1466 1467 end subroutine mixing_init 1468 1469 1470 !> Function to retrieve the number of coefficients 1471 !! calculated in this iteration. 1472 !! This is so external routines can query the size 1473 !! of the arrays used. 1474 !! 1475 !! @param[in] mix the used mixer 1476 function mixing_ncoeff( mix ) result(n) 1477 type(tMixer), intent(in) :: mix 1478 integer :: n 1479 1480 n = 0 1481 1482 select case ( mix%m ) 1483 case ( MIX_PULAY ) 1484 n = n_items(mix%stack(2)) 1485 case ( MIX_BROYDEN ) 1486 n = n_items(mix%stack(2)) 1487 end select 1488 1489 end function mixing_ncoeff 1490 1491 1492 !> Calculate the mixing coefficients for the 1493 !! current mixer 1494 !! 1495 !! @param[in] mix the current mixer 1496 !! @param[in] n the number of elements used to calculate 1497 !! the coefficients 1498 !! @param[in] xin the input value 1499 !! @param[in] F xout - xin, (residual) 1500 !! @param[out] coeff the coefficients 1501 subroutine mixing_coeff( mix, n, xin, F, coeff) 1502 1503 use parallel, only: IONode 1504 1505 type(tMixer), intent(inout) :: mix 1506 integer, intent(in) :: n 1507 real(dp), intent(in) :: xin(n), F(n) 1508 real(dp), intent(out) :: coeff(:) 1509 integer :: ncoeff 1510 1511 ncoeff = size(coeff) 1512 1513 if ( ncoeff < mixing_ncoeff(mix) ) then 1514 write(*,'(a)') 'mix: Error in calculating coefficients' 1515 ! Do not allow this... 1516 return 1517 end if 1518 1519 select case ( mix%m ) 1520 case ( MIX_LINEAR ) 1521 call linear_coeff() 1522 case ( MIX_PULAY ) 1523 call pulay_coeff() 1524 case ( MIX_BROYDEN ) 1525 call broyden_coeff() 1526 end select 1527 1528 contains 1529 1530 subroutine linear_coeff() 1531 integer :: i 1532 do i = 1 , ncoeff 1533 coeff(i) = 0._dp 1534 end do 1535 end subroutine linear_coeff 1536 1537 subroutine pulay_coeff() 1538 integer :: ns, nh, nmax 1539 integer :: i, j, info 1540 logical :: lreturn 1541 1542 ! Calculation quantities 1543 real(dp) :: dnorm, G 1544 real(dp), pointer :: res(:), rres(:), rres1(:), rres2(:) 1545 real(dp), allocatable :: A(:,:), Ainv(:,:) 1546 1547 ns = n_items(mix%stack(1)) 1548 nmax = max_size(mix%stack(1)) 1549 nh = n_items(mix%stack(2)) 1550 1551 lreturn = .false. 1552 1553 ! Easy check for initial step... 1554 select case ( mix%v ) 1555 case ( 0 , 2 ) ! Stable Pulay 1556 1557 lreturn = ns == 1 1558 1559 case ( 1 , 3 ) ! Guaranteed Pulay 1560 1561 lreturn = mod(current_itt(mix), 2) == 1 1562 1563 end select 1564 1565 ! In case we return we are actually doing 1566 ! linear mixing 1567 if ( lreturn ) return 1568 1569 1570 ! Print out number of currently used history steps 1571 if ( debug_mix ) & 1572 write(*,'(a,2(a,i0))') trim(debug_msg), & 1573 ' n_hist = ',ns, ' / ',nmax 1574 1575 ! Allocate arrays for calculating the 1576 ! coefficients 1577 allocate(A(nh,nh), Ainv(nh,nh)) 1578 1579 1580 ! Calculate A_ij coefficients for inversion 1581 do i = 1 , nh 1582 1583 ! Get RRes[i] array 1584 rres1 => getstackval(mix, 2, i) 1585 1586 do j = 1 , i - 1 1587 1588 ! Get RRes[j] array 1589 rres2 => getstackval(mix, 2, j) 1590 1591 ! A(i,j) = A(j,i) = norm(RRes[i],RRes[j]) 1592 A(i,j) = norm(n, rres1, rres2) 1593 A(j,i) = A(i,j) 1594 1595 end do 1596 1597 ! Diagonal 1598 A(i,i) = norm(n, rres1, rres1) 1599 1600 end do 1601 1602#ifdef MPI 1603 ! Global operations, but only for the non-extended entries 1604 call MPI_AllReduce(A(1,1),Ainv(1,1),nh*nh, & 1605 MPI_double_precision, MPI_Sum, & 1606 mix%Comm, i) 1607 ! copy over reduced arrays 1608 A(:,:) = Ainv(:,:) 1609#endif 1610 1611 ! Get inverse of matrix 1612 select case ( mix%v ) 1613 case ( 0 , 1 ) 1614 1615 call inverse(nh, A, Ainv, info) 1616 1617 if ( info /= 0 ) then 1618 1619 ! We will first try the SVD routine 1620 call svd(nh, A, Ainv, mix%rv(I_SVD_COND), j, info) 1621 1622 ! only inform if we should not use SVD per default 1623 if ( IONode ) & 1624 write(*,'(2a,i0,a,i0)') trim(debug_msg), & 1625 ' Pulay -- inversion failed, > SVD [rank/size] ', & 1626 j, ' / ', nh 1627 1628 end if 1629 1630 case ( 2 , 3 ) 1631 1632 ! We forcefully use the SVD routine 1633 call svd(nh, A, Ainv, mix%rv(I_SVD_COND), j, info) 1634 1635 end select 1636 1637 1638 ! NOTE, although mix%stack(1) contains 1639 ! the x[i] - x[i-1], the tip of the stack 1640 ! contains F[i]! 1641 ! res == F[i] 1642 res => getstackval(mix, 1) 1643 1644 ! Initialize the coefficients 1645 do i = 1 , nh 1646 coeff(i) = 0._dp 1647 end do 1648 1649 if ( info == 0 ) then 1650 1651 ! Calculate the coefficients on all processors 1652 do j = 1 , nh 1653 1654 ! res == F[i] 1655 ! rres == F[j+1] - F[j] 1656 rres => getstackval(mix, 2, j) 1657 dnorm = norm(n, rres, res) 1658 1659 do i = 1 , nh 1660 coeff(i) = coeff(i) - Ainv(i,j) * dnorm 1661 end do 1662 1663 end do 1664 1665#ifdef MPI 1666 ! Reduce the coefficients 1667 call MPI_AllReduce(coeff(1),A(1,1),nh, & 1668 MPI_double_precision, MPI_Sum, & 1669 mix%Comm, i) 1670 do i = 1 , nh 1671 coeff(i) = A(i,1) 1672 end do 1673#endif 1674 1675 else 1676 1677 info = 0 1678 1679 ! reset to linear mixing 1680 write(*,'(2a)') trim(debug_msg), & 1681 ' Pulay -- inversion failed, SVD failed, > linear' 1682 1683 end if 1684 1685 ! Clean up memory 1686 deallocate(A, Ainv) 1687 1688 end subroutine pulay_coeff 1689 1690 subroutine broyden_coeff() 1691 integer :: ns, nh, nmax 1692 integer :: i, j, k, info 1693 1694 ! Calculation quantities 1695 real(dp) :: dnorm, dtmp 1696 real(dp), pointer :: w(:), res(:), rres(:), rres1(:), rres2(:) 1697 real(dp), allocatable :: c(:), A(:,:), Ainv(:,:) 1698 1699 ns = n_items(mix%stack(1)) 1700 nmax = max_size(mix%stack(1)) 1701 nh = n_items(mix%stack(2)) 1702 1703 ! Easy check for initial step... 1704 if ( ns == 1 ) then 1705 1706 ! reset 1707 coeff = 0._dp 1708 1709 return 1710 1711 end if 1712 1713 1714 ! Print out number of currently used history steps 1715 if ( debug_mix ) & 1716 write(*,'(a,2(a,i0))') trim(debug_msg), & 1717 ' n_hist = ',ns, ' / ',nmax 1718 1719 1720 ! This is the modified Broyden algorithm... 1721 1722 ! Retrieve the previous weights 1723 w => mix%rv(3:2+nh) 1724 select case ( mix%v ) 1725 case ( 2 ) ! Unity Broyden 1726 1727 w(nh) = 1._dp 1728 1729 case ( 1 ) ! RMS Broyden 1730 1731 dnorm = norm(n, F, F) 1732#ifdef MPI 1733 call MPI_AllReduce(dnorm, dtmp, 1, & 1734 MPI_Double_Precision, MPI_Max, & 1735 mix%Comm, i) 1736 dnorm = dtmp 1737#endif 1738 w(:) = 1._dp / sqrt(dnorm) 1739 if ( debug_mix ) & 1740 write(*,'(2(a,e10.4))') & 1741 trim(debug_msg)//' weight = ',w(1), & 1742 ' , norm = ',dnorm 1743 1744 case ( 0 ) ! Varying weight 1745 1746 dnorm = 0._dp 1747!$OMP parallel do default(shared), private(i), & 1748!$OMP& reduction(max:dnorm) 1749 do i = 1 , n 1750 dnorm = max( dnorm, abs(F(i)) ) 1751 end do 1752!$OMP end parallel do 1753#ifdef MPI 1754 call MPI_AllReduce(dnorm, dtmp, 1, & 1755 MPI_Double_Precision, MPI_Max, & 1756 mix%Comm, i) 1757 dnorm = dtmp 1758#endif 1759 ! Problay 0.2 should be changed to user-defined 1760 w(nh) = exp( 1._dp / (dnorm + 0.2_dp) ) 1761 if ( debug_mix ) & 1762 write(*,'(2a,1000(tr1,e10.4))') & 1763 trim(debug_msg),' weights = ',w(1:nh) 1764 1765 end select 1766 1767 ! Allocate arrays used 1768 allocate(c(nh)) 1769 allocate(A(nh,nh), Ainv(nh,nh)) 1770 1771 1772 ! < RRes[i] | Res[n] > 1773 do i = 1 , nh 1774 rres => getstackval(mix, 2, i) 1775 c(i) = norm(n, rres, F) 1776 end do 1777#ifdef MPI 1778 call MPI_AllReduce(c(1), A(1,1), nh, & 1779 MPI_Double_Precision, MPI_Sum, & 1780 mix%Comm, i) 1781 do i = 1 , nh 1782 c(i) = A(i,1) 1783 end do 1784#endif 1785 1786 ! Create A_ij coefficients for inversion 1787 do i = 1 , nh 1788 1789 ! Get RRes[i] array 1790 rres1 => getstackval(mix, 2, i) 1791 1792 do j = 1 , i -1 1793 1794 ! Get RRes[j] array 1795 rres2 => getstackval(mix, 2, j) 1796 1797 ! A(i,j) = A(j,i) = dot_product(RRes[i],RRes[j]) 1798 A(i,j) = w(i) * w(j) * norm(n, rres1, rres2) 1799 A(j,i) = A(i,j) 1800 1801 end do 1802 1803 ! Do the diagonal term 1804 A(i,i) = w(i) * w(i) * norm(n, rres1, rres1) 1805 1806 end do 1807 1808#ifdef MPI 1809 call MPI_AllReduce(A(1,1), Ainv(1,1), nh*nh, & 1810 MPI_double_precision, MPI_Sum, & 1811 mix%Comm, i) 1812 A(:,:) = Ainv(:,:) 1813#endif 1814 1815 ! Add the diagonal term 1816 ! This should also prevent it from being 1817 ! singular (unless mix%rv(2) == 0) 1818 do i = 1 , nh 1819 A(i,i) = mix%rv(2) ** 2 + A(i,i) 1820 end do 1821 1822 ! Calculate the inverse 1823 call inverse(nh, A, Ainv, info) 1824 1825 if ( info /= 0 ) then 1826 1827 ! We will first try the SVD routine 1828 call svd(nh, A, Ainv, mix%rv(I_SVD_COND), j, info) 1829 1830 ! only inform if we should not use SVD per default 1831 if ( IONode ) & 1832 write(*,'(2a,i0,a,i0)') trim(debug_msg), & 1833 ' Broyden -- inversion failed, > SVD [rank/size] ', & 1834 j, ' / ', nh 1835 1836 end if 1837 1838 do i = 1 , nh 1839 coeff(i) = 0._dp 1840 end do 1841 1842 if ( info == 0 ) then 1843 1844 ! Calculate the coefficients... 1845 do i = 1 , nh 1846 1847 do j = 1 , nh 1848 ! Ainv should be symmetric (A is) 1849 coeff(i) = coeff(i) + w(j) * c(j) * Ainv(j,i) 1850 end do 1851 1852 ! Calculate correct weight... 1853 coeff(i) = - w(i) * coeff(i) 1854 1855 end do 1856 1857 else 1858 1859 ! reset to linear mixing 1860 write(*,'(2a)') trim(debug_msg), & 1861 ' Broyden -- inversion failed, SVD failed, > linear' 1862 1863 end if 1864 1865 deallocate(A, Ainv) 1866 1867 end subroutine broyden_coeff 1868 1869 end subroutine mixing_coeff 1870 1871 1872 1873 !> Calculate the guess for the next iteration 1874 !! 1875 !! Note this gets passed the coefficients. Hence, 1876 !! they may be calculated from another set of history 1877 !! steps. 1878 !! This may be useful in certain situations. 1879 !! 1880 !! @param[in] mix the current mixer 1881 !! @param[in] n the number of elements used to calculate 1882 !! the coefficients 1883 !! @param[in] xin the input value 1884 !! @param[in] F the xin residual 1885 !! @param[out] xnext the input for the following iteration 1886 !! @param[in] coeff the coefficients 1887 subroutine mixing_calc_next( mix, n, xin, F, xnext, coeff ) 1888 ! The current mixing method 1889 type(tMixer), pointer :: mix 1890 integer, intent(in) :: n 1891 real(dp), intent(in) :: xin(n) 1892 real(dp), intent(in) :: F(n) 1893 real(dp), intent(out) :: xnext(n) 1894 real(dp), intent(in) :: coeff(:) 1895 1896 select case ( mix%m ) 1897 case ( MIX_LINEAR ) 1898 call mixing_linear() 1899 case ( MIX_PULAY ) 1900 call mixing_pulay() 1901 case ( MIX_BROYDEN ) 1902 call mixing_broyden() 1903 end select 1904 1905 contains 1906 1907 subroutine mixing_linear() 1908 integer :: i 1909 real(dp) :: w 1910 w = mix%w 1911 1912 if ( debug_mix ) write(*,'(2a,e10.4)') & 1913 trim(debug_msg),' alpha = ', w 1914 1915!$OMP parallel do default(shared), private(i) 1916 do i = 1 , n 1917 xnext(i) = xin(i) + w * F(i) 1918 end do 1919!$OMP end parallel do 1920 1921 end subroutine mixing_linear 1922 1923 subroutine mixing_pulay() 1924 integer :: ns, nh 1925 integer :: i, j 1926 logical :: lreturn 1927 real(dp) :: G 1928 real(dp), pointer :: res(:), rres(:) 1929 1930 ns = n_items(mix%stack(1)) 1931 nh = size(coeff) 1932 if ( nh /= n_items(mix%stack(2)) ) then 1933 write(*,'(a)') 'mix: Error in mixing of Pulay' 1934 xnext = 0._dp 1935 return 1936 end if 1937 1938 ! Easy check for initial step... 1939 select case ( mix%v ) 1940 case ( 0 , 2 ) ! Stable Pulay 1941 1942 lreturn = ns == 1 1943 1944 if ( lreturn .and. debug_mix ) & 1945 write(*,'(2a,e10.4)') trim(debug_msg), & 1946 ' Pulay (initial), alpha = ', mix%rv(1) 1947 1948 case ( 1 , 3 ) ! Guaranteed Pulay 1949 1950 lreturn = mod(current_itt(mix), 2) == 1 1951 1952 if ( lreturn .and. debug_mix ) & 1953 write(*,'(2a,e10.4)') trim(debug_msg), & 1954 ' Direct mixing, alpha = ', mix%rv(1) 1955 1956 end select 1957 1958 ! In case we return we are actually doing 1959 ! linear mixing 1960 if ( lreturn ) then 1961 1962 ! We are doing a linear mixing 1963!$OMP parallel do default(shared), private(i) 1964 do i = 1 , n 1965 xnext(i) = xin(i) + F(i) * mix%rv(1) 1966 end do 1967!$OMP end parallel do 1968 1969 return 1970 1971 end if 1972 1973 ! Get the linear mixing term... 1974 G = mix%w 1975 1976 ! if debugging print out the different variables 1977 if ( debug_mix ) then 1978 write(*,'(2a,f10.6,a,e10.4,a,100(tr1,e10.4))') & 1979 trim(debug_msg),& 1980 ' G = ',G,', sum(alpha) = ',sum(coeff), & 1981 ', alpha = ',coeff 1982 end if 1983 1984 1985!$OMP parallel default(shared), private(i, j) 1986 1987 ! x^opt[i+1] = x[i] + G F[i] 1988!$OMP do 1989 do i = 1 , n 1990 xnext(i) = xin(i) + G * F(i) 1991 end do 1992!$OMP end do 1993 1994 do j = 1 , nh 1995 1996 ! res == x[j] - x[j-1] 1997 ! rres == F[j+1] - F[j] 1998!$OMP single 1999 res => getstackval(mix, 1, j) 2000 rres => getstackval(mix, 2, j) 2001!$OMP end single ! implicit barrier 2002 2003 ! x^opt[i+1] = x^opt[i+1] + 2004 ! alpha_j ( x[j] - x[j-1] + G (F[j+1] - F[j]) ) 2005!$OMP do 2006 do i = 1 , n 2007 xnext(i) = xnext(i) + coeff(j) * ( res(i) + G * rres(i) ) 2008 end do 2009!$OMP end do 2010 2011 end do 2012 2013!$OMP end parallel 2014 2015 end subroutine mixing_pulay 2016 2017 subroutine mixing_broyden() 2018 integer :: ns, nh 2019 integer :: i, j 2020 real(dp) :: G 2021 real(dp), pointer :: res(:), rres(:) 2022 2023 ns = n_items(mix%stack(1)) 2024 nh = size(coeff) 2025 if ( nh /= n_items(mix%stack(2)) ) then 2026 write(*,'(a)') 'mix: Error in mixing of Broyden' 2027 xnext = 0._dp 2028 return 2029 end if 2030 2031 if ( ns == 1 ) then 2032 if ( debug_mix ) & 2033 write(*,'(2a,e10.4)') trim(debug_msg), & 2034 ' Broyden (initial), alpha = ', mix%rv(1) 2035 2036 ! We are doing a linear mixing 2037!$OMP parallel do default(shared), private(i) 2038 do i = 1 , n 2039 xnext(i) = xin(i) + F(i) * mix%rv(1) 2040 end do 2041!$OMP end parallel do 2042 2043 return 2044 2045 end if 2046 2047 ! Get the inverse Jacobian term... 2048 G = mix%w 2049 2050 ! if debugging print out the different variables 2051 if ( debug_mix ) then 2052 write(*,'(2a,f10.6,a,e10.4,a,100(tr1,e10.4))') & 2053 trim(debug_msg), ' G = ', G, & 2054 ', sum(coeff) = ',sum(coeff), & 2055 ', coeff = ',coeff 2056 end if 2057 2058 2059!$OMP parallel default(shared), private(i, j) 2060 2061 ! x^opt[i+1] = x[i] + G F[i] 2062!$OMP do 2063 do i = 1 , n 2064 xnext(i) = xin(i) + G * F(i) 2065 end do 2066!$OMP end do 2067 2068 do j = 1 , nh 2069 2070 ! res == x[j] - x[j-1] 2071 ! rres == F[j+1] - F[j] 2072!$OMP single 2073 res => getstackval(mix, 1, j) 2074 rres => getstackval(mix, 2, j) 2075!$OMP end single ! implicit barrier 2076 2077 ! x^opt[i+1] = x^opt[i+1] + 2078 ! alpha_j ( x[j] - x[j-1] + G (F[j+1] - F[j]) ) 2079!$OMP do 2080 do i = 1 , n 2081 xnext(i) = xnext(i) + coeff(j) * ( res(i) + G * rres(i) ) 2082 end do 2083!$OMP end do 2084 2085 end do 2086 2087!$OMP end parallel 2088 2089 end subroutine mixing_broyden 2090 2091 end subroutine mixing_calc_next 2092 2093 2094 !> Finalize the mixing algorithm 2095 !! 2096 !! @param[inout] mix mixer to be finalized 2097 !! @param[in] n size of the input arrays 2098 !! @param[in] xin the input for this iteration 2099 !! @param[in] F the residual for this iteration 2100 !! @param[in] xnext the optimized input for the next iteration 2101 subroutine mixing_finalize(mix, n, xin, F, xnext) 2102 use parallel, only: IONode 2103 2104 type(tMixer), pointer :: mix 2105 integer, intent(in) :: n 2106 real(dp), intent(in) :: xin(n), F(n), xnext(n) 2107 integer :: rsave, is 2108 2109 select case ( mix%m ) 2110 case ( MIX_LINEAR ) 2111 call fin_linear() 2112 case ( MIX_PULAY ) 2113 call fin_pulay() 2114 case ( MIX_BROYDEN ) 2115 call fin_broyden() 2116 end select 2117 2118 2119 ! Fix the action to finalize it.. 2120 if ( mix%restart > 0 ) then 2121 if ( mod(current_itt(mix),mix%restart) == 0 ) then 2122 mix%action = IOR(mix%action, ACTION_RESTART) 2123 end if 2124 end if 2125 2126 2127 ! Check the actual finalization... 2128 ! First check whether we should restart history 2129 if ( IAND(mix%action, ACTION_RESTART) == ACTION_RESTART ) then 2130 2131 ! The user has requested to restart the 2132 ! mixing scheme now 2133 rsave = mix%restart_save 2134 2135 select case ( mix%m ) 2136 case ( MIX_PULAY ) 2137 2138 if ( IONode ) then 2139 write(*,'(a)')'mix: Pulay -- resetting history' 2140 end if 2141 2142 if ( rsave == 0 ) then 2143 do is = 1, size(mix%stack) 2144 call reset(mix%stack(is)) 2145 end do 2146 else 2147 call reset(mix%stack(1),-rsave) 2148 call reset(mix%stack(2),-rsave+1) 2149 end if 2150 2151 case ( MIX_BROYDEN ) 2152 2153 if ( IONode ) then 2154 write(*,'(a)')'mix: Broyden -- resetting history' 2155 end if 2156 2157 if ( rsave == 0 ) then 2158 call reset(mix%stack(1)) 2159 call reset(mix%stack(2)) 2160 call reset(mix%stack(3)) 2161 else 2162 call reset(mix%stack(1),-rsave) 2163 call reset(mix%stack(2),-rsave+1) 2164 end if 2165 2166 end select 2167 2168 if ( allocated(mix%stack) ) then 2169 if ( debug_mix ) & 2170 write(*,'(a,a,i0)') trim(debug_msg), & 2171 ' saved hist = ',n_items(mix%stack(1)) 2172 end if 2173 2174 end if 2175 2176 2177 ! check whether we should change the mixer 2178 if ( IAND(mix%action, ACTION_NEXT) == ACTION_NEXT ) then 2179 2180 call mixing_step( mix ) 2181 2182 end if 2183 2184 contains 2185 2186 subroutine fin_linear() 2187 2188 ! do nothing... 2189 2190 end subroutine fin_linear 2191 2192 subroutine fin_pulay() 2193 2194 integer :: ns 2195 integer :: i 2196 logical :: GR_linear 2197 real(dp), pointer :: res(:), rres(:) 2198 2199 ns = n_items(mix%stack(1)) 2200 2201 select case ( mix%v ) 2202 case ( 0 , 2 ) ! stable Pulay 2203 2204 if ( n_items(mix%stack(3)) == 0 ) then 2205 call push_stack_data(mix%stack(3), n) 2206 end if 2207 2208 res => getstackval(mix, 3) 2209 2210!$OMP parallel do default(shared), private(i) 2211 do i = 1 , n 2212 ! Input: 2213 ! res == x[i-1] + F[i-1] 2214 res(i) = xin(i) + F(i) 2215 ! Output: 2216 ! res == x[i] + F[i] 2217 end do 2218!$OMP end parallel do 2219 2220 case ( 1 , 3 ) ! GR Pulay 2221 2222 GR_linear = mod(current_itt(mix), 2) == 1 2223 2224 if ( n_items(mix%stack(2)) > 0 .and. & 2225 .not. GR_linear ) then 2226 2227 res => getstackval(mix, 1) 2228 rres => getstackval(mix, 2) 2229 2230!$OMP parallel do default(shared), private(i) 2231 do i = 1 , n 2232 ! Input: 2233 ! rres == F[i] - F[i-1] 2234 rres(i) = rres(i) - res(i) 2235 ! Output: 2236 ! rres == - F[i-1] 2237 end do 2238!$OMP end parallel do 2239 2240 call pop(mix%stack(1)) 2241 2242 ! Note that this is Res[i-1] = (F^i-1_out - F^i-1_in) 2243 res => getstackval(mix,1) 2244!$OMP parallel do default(shared), private(i) 2245 do i = 1 , n 2246 res(i) = res(i) - xin(i) + xnext(i) 2247 end do 2248!$OMP end parallel do 2249 2250 end if 2251 2252 end select 2253 2254 end subroutine fin_pulay 2255 2256 subroutine fin_broyden() 2257 2258 integer :: ns, nh 2259 integer :: i 2260 real(dp), pointer :: res(:), rres(:) 2261 2262 ns = current_itt(mix) 2263 nh = n_items(mix%stack(2)) 2264 2265 if ( ns >= 2 .and. n_items(mix%stack(3)) > 0 ) then 2266 2267 ! Update the residual to reflect the input residual 2268 res => getstackval(mix, 3) 2269 2270!$OMP parallel do default(shared), private(i) 2271 do i = 1 , n 2272 ! Input: 2273 ! res == x[i-1] + F[i-1] 2274 res(i) = xin(i) + F(i) 2275 ! Output: 2276 ! res == x[i] + F[i] 2277 end do 2278!$OMP end parallel do 2279 2280 end if 2281 2282 ! Update weights (if necessary) 2283 if ( nh > 1 ) then 2284 do i = 3 , nh + 1 2285 mix%rv(i) = mix%rv(i+1) 2286 end do 2287 end if 2288 2289 end subroutine fin_broyden 2290 2291 end subroutine mixing_finalize 2292 2293 2294 ! Perform the actual mixing... 2295 subroutine mixing_1d( mix, n, xin, F, xnext, nsub) 2296 ! The current mixing method 2297 type(tMixer), pointer :: mix 2298 ! The current step in the SCF and size of arrays 2299 integer, intent(in) :: n 2300 ! x1 == Input function, 2301 ! F1 == Residual from x1 2302 real(dp), intent(in) :: xin(n), F(n) 2303 ! x2 == Next input function 2304 real(dp), intent(inout) :: xnext(n) 2305 ! Number of elements used for calculating the mixing 2306 ! coefficients 2307 integer, intent(in), optional :: nsub 2308 2309 ! Coefficients 2310 integer :: ncoeff 2311 real(dp), allocatable :: coeff(:) 2312 2313 call mixing_init( mix, n, xin, F ) 2314 2315 ncoeff = mixing_ncoeff( mix ) 2316 allocate(coeff(ncoeff)) 2317 2318 ! Calculate coefficients 2319 if ( present(nsub) ) then 2320 call mixing_coeff( mix, nsub, xin, F, coeff) 2321 else 2322 call mixing_coeff( mix, n, xin, F, coeff) 2323 end if 2324 2325 ! Calculate the following output 2326 call mixing_calc_next( mix, n, xin, F, xnext, coeff) 2327 2328 ! Coefficients are not needed anymore... 2329 deallocate(coeff) 2330 2331 ! Finalize the mixer 2332 call mixing_finalize( mix, n, xin, F, xnext) 2333 2334 end subroutine mixing_1d 2335 2336 2337 subroutine mixing_2d( mix, n1, n2, xin, F, xnext, nsub) 2338 2339 type(tMixer), pointer :: mix 2340 integer, intent(in) :: n1, n2 2341 real(dp), intent(in) :: xin(n1,n2), F(n1,n2) 2342 real(dp), intent(inout) :: xnext(n1,n2) 2343 integer, intent(in), optional :: nsub 2344 2345 ! Simple wrapper for 1D 2346 if ( present(nsub) ) then 2347 call mixing_1d( mix, n1*n2 , xin(1,1), F(1,1), xnext(1,1) ,& 2348 nsub = n1 * nsub) 2349 else 2350 call mixing_1d( mix, n1*n2 , xin(1,1), F(1,1), xnext(1,1)) 2351 end if 2352 2353 end subroutine mixing_2d 2354 2355 2356 2357 2358 ! Step the mixing object and ensure that 2359 ! the old history is either copied over or freed 2360 subroutine mixing_step( mix ) 2361 2362 use parallel, only: IONode 2363 2364 type(tMixer), pointer :: mix 2365 2366 type(tMixer), pointer :: next => null() 2367 2368 type(dData1D), pointer :: d1D 2369 integer :: i, is, n, init_itt 2370 logical :: reset_stack, copy_stack 2371 2372 ! First try and 2373 next => mix%next 2374 if ( associated(next) ) then 2375 2376 ! Whether or not the two methods are allowed 2377 ! to share history 2378 copy_stack = mix%m == next%m 2379 ! Only Pulay original and Broyden methods are compatible. 2380 ! The Pulay-GR variant does not have the same data stored. 2381 select case ( mix%m ) 2382 case ( MIX_PULAY ) 2383 select case ( next%m ) 2384 case ( MIX_PULAY ) 2385 copy_stack = mod(mix%v, 2) == mod(next%v, 2) 2386 case ( MIX_BROYDEN ) 2387 copy_stack = mod(mix%v, 2) == 0 2388 end select 2389 case ( MIX_BROYDEN ) 2390 select case ( next%m ) 2391 case ( MIX_PULAY ) 2392 copy_stack = mod(mix%v, 2) == 0 2393 ! otherwise both are broyden 2394 end select 2395 end select 2396 2397 copy_stack = copy_stack .and. allocated(mix%stack) 2398 2399 ! If the two methods are similar 2400 if ( copy_stack ) then 2401 2402 ! They are similar, copy over the history stack 2403 do is = 1 , size(mix%stack) 2404 2405 ! Get maximum size of the current stack, 2406 n = n_items(mix%stack(is)) 2407 2408 ! Note that this will automatically take care of 2409 ! wrap-arounds and delete the unneccesry elements 2410 do i = 1 , n 2411 d1D => get_pointer(mix%stack(is),i) 2412 call push(next%stack(is),d1D) 2413 end do 2414 2415 ! nullify 2416 nullify(d1D) 2417 2418 end do 2419 2420 end if 2421 2422 end if 2423 2424 reset_stack = .true. 2425 if ( associated(next) ) then 2426 if ( associated(next%next, mix) .and. & 2427 next%n_itt > 0 ) then 2428 ! if this is a circular mixing routine 2429 ! we should not reset the history... 2430 reset_stack = .false. 2431 end if 2432 end if 2433 2434 if ( reset_stack ) then 2435 select case ( mix%m ) 2436 case ( MIX_PULAY , MIX_BROYDEN ) 2437 2438 n = size(mix%stack) 2439 do is = 1 , n 2440 call reset(mix%stack(is)) 2441 end do 2442 2443 end select 2444 end if 2445 2446 if ( associated(next) ) then 2447 init_itt = 0 2448 ! Set-up the next mixer 2449 select case ( next%m ) 2450 case ( MIX_PULAY , MIX_BROYDEN ) 2451 init_itt = n_items(next%stack(1)) 2452 end select 2453 next%start_itt = init_itt 2454 next%cur_itt = init_itt 2455 if ( IONode ) then 2456 if ( copy_stack ) then 2457 write(*,'(3a)') trim(debug_msg),' switching mixer (with history) --> ', & 2458 trim(next%name) 2459 else 2460 write(*,'(3a)') trim(debug_msg),' switching mixer --> ', & 2461 trim(next%name) 2462 end if 2463 end if 2464 mix => mix%next 2465 end if 2466 2467 end subroutine mixing_step 2468 2469 2470 2471 ! Calculate the inverse of a matrix 2472 subroutine inverse(n, A, B, info ) 2473 2474 integer, intent(in) :: n 2475 real(dp), intent(in) :: A(n,n) 2476 real(dp), intent(out) :: B(n,n) 2477 integer, intent(out) :: info 2478 2479 integer :: i, j 2480 ! Local arrays 2481 real(dp) :: pm(n,n), work(n*4), err 2482 ! Relative tolerance dependent on the magnitude 2483 ! For now we retain the old tolerance 2484 real(dp), parameter :: etol = 1.e-4_dp 2485 integer :: ipiv(n) 2486 2487 ! initialize info 2488 info = 0 2489 2490 ! simple check and fast return 2491 if ( n == 1 ) then 2492 2493 B(1,1) = 1._dp / A(1,1) 2494 2495 return 2496 2497 end if 2498 2499 call lapack_inv() 2500 if ( info /= 0 ) call simple_inv() 2501 2502 contains 2503 2504 subroutine lapack_inv() 2505 2506 B(:, :) = A(:, :) 2507 2508 call dgetrf(n,n,B,n,ipiv,info) 2509 if ( info /= 0 ) return 2510 2511 call dgetri(n,B,n,ipiv,work,n*4,info) 2512 if ( info /= 0 ) return 2513 2514 ! This sets info appropriately 2515 call check_inv() 2516 2517 end subroutine lapack_inv 2518 2519 subroutine simple_inv() 2520 2521 real(dp) :: x 2522 integer :: k 2523 2524 ! Copy over A 2525 B(:, :) = A(:, :) 2526 2527 do i = 1 , n 2528 if ( B(i,i) == 0._dp ) then 2529 info = -n 2530 return 2531 end if 2532 2533 x = 1._dp / B(i,i) 2534 B(i,i) = 1._dp 2535 do j = 1 , n 2536 B(j,i) = B(j,i) * x 2537 end do 2538 2539 do k = 1 , n 2540 if ( (k-i) /= 0 ) then 2541 x = B(i,k) 2542 B(i,k) = 0._dp 2543 do j = 1 , n 2544 B(j,k) = B(j,k) - B(j,i) * x 2545 end do 2546 end if 2547 end do 2548 end do 2549 2550 ! This sets info appropriately 2551 call check_inv() 2552 2553 end subroutine simple_inv 2554 2555 subroutine check_inv() 2556 2557 ! Check correcteness 2558 pm = matmul(A,B) 2559 2560 do j = 1 , n 2561 do i = 1 , n 2562 if ( i == j ) then 2563 err = pm(i,j) - 1._dp 2564 else 2565 err = pm(i,j) 2566 end if 2567 2568 ! This is pretty strict tolerance! 2569 if ( abs(err) > etol ) then 2570 ! Signal failure in inversion 2571 info = - n - 1 2572 return 2573 end if 2574 2575 end do 2576 end do 2577 2578 end subroutine check_inv 2579 2580 end subroutine inverse 2581 2582 2583 ! Calculate the svd of a matrix 2584 ! With ||(Ax - b)|| << 2585 subroutine svd(n, A, B, cond, rank, info) 2586 2587 integer, intent(in) :: n 2588 real(dp), intent(in) :: A(n,n) 2589 real(dp), intent(out) :: B(n,n) 2590 real(dp), intent(in) :: cond 2591 integer, intent(out) :: rank, info 2592 2593 ! Local arrays 2594 integer :: i 2595 character(len=64) :: fmt 2596 real(dp) :: AA(n,n), S(n), work(n*5) 2597 2598 if ( n == 1 ) then 2599 ! Simple scheme 2600 B(1,1) = 1._dp / A(1,1) 2601 return 2602 end if 2603 2604 ! Copy A matrix 2605 AA(:, :) = A(:, :) 2606 ! setup pseudo inverse solution for minimizing 2607 ! constraints 2608 B(:, :) = 0._dp 2609 do i = 1 , n 2610 B(i,i) = 1._dp 2611 end do 2612 call dgelss(n,n,n,AA,n,B,n,S,cond,rank,work,n*5,info) 2613 2614 ! if debugging print out the different variables 2615 if ( debug_mix ) then 2616 ! also mark the rank 2617 2618 if ( rank == n ) then 2619 ! complete rank 2620 write(*,'(2a,100(tr1,e10.4))') & 2621 trim(debug_msg),' SVD singular = ',S 2622 else 2623 ! this prints the location of the SVD rank, if not full 2624 write(fmt,'(i0,2a)') rank, '(tr1,e10.4),'' >'',100(tr1,e10.4)' 2625 write(*,'(2a,'//trim(fmt)//')') & 2626 trim(debug_msg),' SVD singular = ',S 2627 end if 2628 end if 2629 2630 ! One could potentially search a smaller part of the iterations 2631 ! for a complete satisfactory SVD solution. 2632 ! However, typically the singular values are because the 2633 ! two latest iterations which means that one will end up with 2634 ! the linear mixing scheme with only a coefficient for the latest 2635 ! iteration. 2636 ! This, is per see not a bad choice since it forces a linear mixing 2637 ! in the end. However, the SVD solution seems better if one tries to 2638 ! really push the tolerances. 2639 2640 ! One could do this (which requires the routine to be recursive) 2641 ! if ( rank < n ) then 2642 ! call svd(n-1,A(1:n-1,1:n-1), B(1:n-1, 1:n-1), cond, rank, info) 2643 ! ! zero out everything else 2644 ! B(n,:) = 0._dp 2645 ! B(:,n) = 0._dp 2646 ! end if 2647 2648 end subroutine svd 2649 2650 2651 ! ************************************************* 2652 ! * Helper routines * 2653 ! * LOCAL * 2654 ! ************************************************* 2655 2656 ! Returns the value array from the stack(:) 2657 ! Returns this array: 2658 ! mix%stack(sidx)(hidx) ! defaults to the last item 2659 function getstackval(mix,sidx,hidx) result(d1) 2660 type(tMixer), intent(in) :: mix 2661 integer, intent(in) :: sidx 2662 integer, intent(in), optional :: hidx 2663 2664 real(dp), pointer :: d1(:) 2665 2666 type(dData1D), pointer :: dD1 2667 2668 if ( present(hidx) ) then 2669 dD1 => get_pointer(mix%stack(sidx),hidx) 2670 else 2671 dD1 => get_pointer(mix%stack(sidx), & 2672 n_items(mix%stack(sidx))) 2673 end if 2674 2675 d1 => val(dD1) 2676 2677 end function getstackval 2678 2679 2680 ! Returns true if the following 2681 ! "advanced" mixer is 'method' 2682 function is_next(mix,method,next) result(bool) 2683 type(tMixer), intent(in), target :: mix 2684 integer, intent(in) :: method 2685 type(tMixer), pointer, optional :: next 2686 2687 logical :: bool 2688 2689 type(tMixer), pointer :: m 2690 2691 bool = .false. 2692 m => mix%next 2693 do while ( associated(m) ) 2694 2695 if ( m%m == MIX_LINEAR ) then 2696 m => m%next 2697 else if ( m%m == method ) then 2698 bool = .true. 2699 exit 2700 else 2701 ! Quit if it does not do anything 2702 exit 2703 end if 2704 2705 ! this will prevent cyclic combinations 2706 if ( associated(m,mix) ) exit 2707 2708 end do 2709 2710 if ( present(next) ) then 2711 next => m 2712 end if 2713 2714 end function is_next 2715 2716 2717 !> Get current iteration count 2718 !! 2719 !! This is abstracted because the initial iteration 2720 !! and the current iteration may be uniquely defined. 2721 function current_itt(mix) result(itt) 2722 type(tMixer), intent(in) :: mix 2723 integer :: itt 2724 2725 itt = mix%cur_itt - mix%start_itt 2726 2727 end function current_itt 2728 2729 2730 ! Stack handling routines 2731 function stack_check(stack,n) result(check) 2732 type(Fstack_dData1D), intent(inout) :: stack 2733 integer, intent(in) :: n 2734 logical :: check 2735 2736 ! Local arrays 2737 type(dData1D), pointer :: dD1 2738 2739 if ( n_items(stack) == 0 ) then 2740 check = .true. 2741 else 2742 2743 ! Check that the stack stored arrays are 2744 ! of same size... 2745 2746 dD1 => get_pointer(stack,1) 2747 check = n == size(dD1) 2748 2749 end if 2750 2751 end function stack_check 2752 2753 subroutine push_stack_data(s_F,n) 2754 type(Fstack_dData1D), intent(inout) :: s_F 2755 integer, intent(in) :: n 2756 2757 type(dData1D) :: dD1 2758 2759 if ( .not. stack_check(s_F,n) ) then 2760 call die('mixing: history has changed size...') 2761 end if 2762 2763 call newdData1D(dD1, n, '(F)') 2764 2765 ! Push the data to the stack 2766 call push(s_F,dD1) 2767 2768 ! Delete double reference 2769 call delete(dD1) 2770 2771 end subroutine push_stack_data 2772 2773 subroutine push_F(s_F,n,F,fact) 2774 type(Fstack_dData1D), intent(inout) :: s_F 2775 integer, intent(in) :: n 2776 real(dp), intent(in) :: F(n) 2777 real(dp), intent(in), optional :: fact 2778 2779 type(dData1D) :: dD1 2780 real(dp), pointer :: sF(:) 2781 integer :: in, ns 2782 integer :: i 2783 2784 if ( .not. stack_check(s_F,n) ) then 2785 call die('mixing: history has changed size...') 2786 end if 2787 2788 in = n_items(s_F) 2789 ns = max_size(s_F) 2790 2791 if ( in == ns ) then 2792 2793 ! we have to cycle the storage 2794 call get(s_F,1,dD1) 2795 2796 else 2797 2798 call newdData1D(dD1, n, '(F)') 2799 2800 end if 2801 2802 sF => val(dD1) 2803 2804 if ( present(fact) ) then 2805!$OMP parallel do default(shared), private(i) 2806 do i = 1 , n 2807 sF(i) = F(i) * fact 2808 end do 2809!$OMP end parallel do 2810 else 2811 call dcopy(n, F, 1, sF, 1) 2812 end if 2813 2814 ! Push the data to the stack 2815 call push(s_F,dD1) 2816 2817 ! Delete double reference 2818 call delete(dD1) 2819 2820 end subroutine push_F 2821 2822 subroutine update_F(s_F,n,F) 2823 type(Fstack_dData1D), intent(inout) :: s_F 2824 integer, intent(in) :: n 2825 real(dp), intent(in) :: F(n) 2826 2827 type(dData1D), pointer :: dD1 2828 real(dp), pointer :: FF(:) 2829 integer :: in 2830 2831 if ( .not. stack_check(s_F,n) ) then 2832 call die('mixing: history has changed size...') 2833 end if 2834 2835 in = n_items(s_F) 2836 2837 if ( in == 0 ) then 2838 2839 ! We need to add it as it does not exist 2840 call push_F(s_F,n,F) 2841 2842 else 2843 2844 ! we have an entry, update the latest 2845 dD1 => get_pointer(s_F,in) 2846 2847 FF => val(dD1) 2848 2849 call dcopy(n, F, 1, FF, 1) 2850 2851 end if 2852 2853 end subroutine update_F 2854 2855 subroutine push_diff(s_rres,s_res,alpha) 2856 type(Fstack_dData1D), intent(inout) :: s_rres 2857 type(Fstack_dData1D), intent(in) :: s_res 2858 real(dp), intent(in), optional :: alpha 2859 2860 type(dData1D) :: dD1 2861 type(dData1D), pointer :: pD1 2862 real(dp), pointer :: res1(:), res2(:), rres(:) 2863 integer :: in, ns, i, n 2864 2865 if ( n_items(s_res) < 2 ) then 2866 call die('mixing: Residual residuals cannot be calculated, & 2867 &inferior residual size.') 2868 end if 2869 2870 in = n_items(s_res) 2871 2872 ! First get the value of in 2873 pD1 => get_pointer(s_res,in-1) 2874 res1 => val(pD1) 2875 ! get the value of in 2876 pD1 => get_pointer(s_res,in) 2877 res2 => val(pD1) 2878 2879 in = n_items(s_rres) 2880 ns = max_size(s_rres) 2881 2882 if ( in == ns ) then 2883 2884 ! we have to cycle the storage 2885 call get(s_rres,1,dD1) 2886 2887 else 2888 2889 call newdData1D(dD1, size(res1), '(res)') 2890 2891 end if 2892 2893 ! Get the residual of the residual 2894 rres => val(dD1) 2895 n = size(rres) 2896 2897 if ( present(alpha) ) then 2898!$OMP parallel do default(shared), private(i) 2899 do i = 1 , n 2900 rres(i) = (res2(i) - res1(i)) * alpha 2901 end do 2902!$OMP end parallel do 2903 else 2904!$OMP parallel do default(shared), private(i) 2905 do i = 1 , n 2906 rres(i) = res2(i) - res1(i) 2907 end do 2908!$OMP end parallel do 2909 end if 2910 2911 ! Push the data to the stack 2912 call push(s_rres,dD1) 2913 2914 ! Delete double reference 2915 call delete(dD1) 2916 2917 end subroutine push_diff 2918 2919 2920 !> Calculate the norm of two arrays 2921 function norm(n, x1, x2) 2922 integer, intent(in) :: n 2923 real(dp), intent(in) :: x1(n), x2(n) 2924 real(dp) :: norm 2925 2926 ! Currently we use an external routine 2927 integer :: i 2928 2929 ! Calculate dot product 2930 2931 norm = 0._dp 2932!$OMP parallel do default(shared), private(i) & 2933!$OMP& reduction(+:norm) 2934 do i = 1 , n 2935 norm = norm + x1(i) * x2(i) 2936 end do 2937!$OMP end parallel do 2938 2939 end function norm 2940 2941 2942end module m_mixing 2943 2944