← Back to branch summary

~siesta-pseudos-bases/siesta/trunk-psml

« back to all changes in this revision

Viewing changes to Src/m_mixing.F90

  • Committer: Alberto Garcia
  • Date: 2019-09-02 14:09:43 UTC
  • mfrom: (427.6.323 trunk)
  • Revision ID: albertog@icmab.es-20190902140943-mzmbe1jacgefpgxw
Sync to trunk-776 (notably nc/soc wavefunction support)


Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/m_mixing.F90
423
423
 
424
424
         ! Do options so that a pulay option may refer to
425
425
         ! the actual names of the constants
 
426
         
426
427
         if ( m%m == MIX_PULAY ) then
427
 
 
428
 
            ! The linear mixing weight
429
 
            if ( leqi(opt,'weight.linear') &
430
 
                 .or. leqi(opt,'w.linear') ) then
431
 
               
432
 
               m%rv(1) = fdf_breals(pline,1)
433
 
 
434
 
            else if ( leqi(opt,'svd.cond') ) then
435
 
 
436
 
               ! This is only applicable to the Pulay
437
 
               ! mixing scheme...
438
 
               
439
 
               m%rv(I_SVD_COND) = fdf_bvalues(pline,1)
440
 
               
441
 
            end if
442
 
 
 
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
           
443
465
         end if
444
 
 
445
466
         
446
467
         ! Generic options for all advanced methods...
447
468
         if ( leqi(opt,'next.p') ) then
473
494
    type(tMixer), intent(inout) :: mix
474
495
    integer :: n
475
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
 
476
501
    ! Correct amount of history in the mixing.
477
502
    if ( 0 < mix%restart .and. &
478
503
         mix%restart < mix%n_hist ) then
518
543
       
519
544
       ! allocate temporary array
520
545
       mix%n_hist = max(2, mix%n_hist)
521
 
       n = 1 + mix%n_hist
 
546
       n = 2 + mix%n_hist
522
547
       allocate(mix%rv(I_SVD_COND:n))
523
548
       mix%rv(1:n) = mix%w
524
549
 
702
727
    type(tMixer), intent(in), target :: mixers(:)
703
728
 
704
729
    type(tMixer), pointer :: m
705
 
    character(len=50) :: fmt
 
730
    character(len=64) :: fmt
706
731
 
707
732
    logical :: bool
708
733
    integer :: i
795
820
          write(*,'(2a,t50,''= '',i0)') trim(fmt), &
796
821
               '    History steps',m%n_hist
797
822
          write(*,'(2a,t50,''= '',f12.6)') trim(fmt), &
798
 
               '    Jacobian weight',m%w
799
 
          write(*,'(2a,t50,''= '',f12.6)') trim(fmt), &
800
 
               '    Weight prime',m%rv(1)
 
823
               '    Linear mixing weight',m%rv(1)
 
824
          write(*,'(2a,t50,''= '',f12.6)') trim(fmt), &
 
825
               '    Inverse Jacobian weight',m%w
 
826
          write(*,'(2a,t50,''= '',f12.6)') trim(fmt), &
 
827
               '    Weight prime',m%rv(2)
 
828
          write(*,'(2a,t50,''= '',e10.4)') trim(fmt), &
 
829
               '    SVD condition',m%rv(I_SVD_COND)
801
830
          if ( m%rv(I_P_NEXT) > 0._dp ) then
802
831
             write(*,'(2a,t50,''= '',f6.4)') trim(fmt), &
803
832
                  '    Step mixer parameter',m%rv(I_P_NEXT)
921
950
       ! For Broyden this is the inverse Jacobian
922
951
       write(*,'(t3,a,f6.4)') 'weight ', m%w
923
952
       select case ( m%m )
924
 
       case ( MIX_PULAY, MIX_BROYDEN )
925
 
          write(*,'(t3,a,f6.4)') 'weight.linear ', m%rv(1)
 
953
       case ( MIX_PULAY )
 
954
         write(*,'(t3,a,f6.4)') 'weight.linear ', m%rv(1)
 
955
       case ( MIX_BROYDEN )
 
956
         write(*,'(t3,a,f6.4)') 'weight.linear ', m%rv(1)
 
957
         write(*,'(t3,a,f6.4)') 'weight.prime ', m%rv(2)
926
958
       end select
927
959
          
928
960
       if ( m%n_hist > 0 ) then
1685
1717
      ! This is the modified Broyden algorithm...
1686
1718
 
1687
1719
      ! Retrieve the previous weights
1688
 
      w => mix%rv(2:1+nh)
 
1720
      w => mix%rv(3:2+nh)
1689
1721
      select case ( mix%v )
1690
1722
      case ( 2 ) ! Unity Broyden
1691
1723
         
1779
1811
 
1780
1812
      ! Add the diagonal term
1781
1813
      ! This should also prevent it from being
1782
 
      ! singular (unless mix%w == 0)
 
1814
      ! singular (unless mix%rv(2) == 0)
1783
1815
      do i = 1 , nh
1784
 
         A(i,i) = mix%rv(1) ** 2 + A(i,i)
 
1816
         A(i,i) = mix%rv(2) ** 2 + A(i,i)
1785
1817
      end do
1786
1818
 
1787
1819
      ! Calculate the inverse
2006
2038
         
2007
2039
      end if
2008
2040
 
2009
 
      ! Get the linear mixing term...
 
2041
      ! Get the inverse Jacobian term...
2010
2042
      G = mix%w
2011
2043
 
2012
2044
      ! if debugging print out the different variables
2077
2109
 
2078
2110
 
2079
2111
    ! Fix the action to finalize it..
2080
 
    if ( mix%restart > 0 .and. &
2081
 
         mod(current_itt(mix),mix%restart) == 0 ) then
2082
 
 
2083
 
       mix%action = IOR(mix%action, ACTION_RESTART)
 
2112
    if ( mix%restart > 0 ) then
 
2113
      if ( mod(current_itt(mix),mix%restart) == 0 ) then
 
2114
        mix%action = IOR(mix%action, ACTION_RESTART)
 
2115
      end if
2084
2116
 
2085
2117
    end if
2086
2118
 
2242
2274
 
2243
2275
      ! Update weights (if necessary)
2244
2276
      if ( nh > 1 ) then
2245
 
         do i = 2 , nh
 
2277
         do i = 3 , nh + 1
2246
2278
            mix%rv(i) = mix%rv(i+1)
2247
2279
         end do
2248
2280
      end if