1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8 MODULE m_mixer 9 private 10 public :: mixer 11 12 CONTAINS 13 14 subroutine mixer( iscf ) 15 16 use precision, only: dp 17 use siesta_options, only: mixH, mix_scf_first 18 use siesta_options, only: mullipop, muldeb, orbmoms 19 ! Spin mixing options 20 use m_mixing_scf, only: MIX_SPIN_ALL, MIX_SPIN_SPINOR 21 use m_mixing_scf, only: MIX_SPIN_SUM, MIX_SPIN_SUM_DIFF 22 use m_mixing_scf, only: mix_spin 23 24 use sparse_matrices, only: Dold, Dscf, Hold, H, S 25 use sparse_matrices, only: maxnh, numh, listhptr, listh 26 use siesta_geom, only: na_u, isa 27 use atomlist, only: iaorb, iphorb, lasto, no_u, no_l 28 use atomlist, only: indxuo 29 30 use m_mixing_scf, only: scf_mix 31 use m_mixing, only: mixing 32 33 use m_spin, only: spin 34 35 use parallel, only: IONode 36 37 implicit none 38 39 real(dp), pointer :: Xin(:,:), Xout(:,:) 40 real(dp), allocatable :: F(:,:) 41 42 integer, intent(in) :: iscf 43 44 real(dp) :: dtmp 45 integer :: iiscf, i, is 46 integer :: nspin_mix 47 48 external :: mulliken 49!-------------------------------------------------------------------- BEGIN 50 51 call timer( 'MIXER', 1 ) 52 53 if ( mixH ) then 54 ! Mix Hamiltonian 55 Xin => Hold 56 Xout => H 57 nspin_mix = spin%H 58 else 59 ! Mix density matrix 60 Xin => Dold 61 Xout => Dscf 62 nspin_mix = spin%DM 63 end if 64 65 66 ! Create residual function to minimize 67 allocate(F(maxnh,nspin_mix)) 68 69!$OMP parallel default(shared), private(dtmp,i,is) 70 71 ! prepare input... 72 select case ( mix_spin ) 73 case ( MIX_SPIN_ALL , MIX_SPIN_SPINOR ) 74 ! Xin and Xout are unchanged... 75 case ( MIX_SPIN_SUM , MIX_SPIN_SUM_DIFF ) 76 77 ! transfer spin density to 78 ! spin-sum and spin-difference 79!$OMP do 80 do i = 1 , size(Xin,1) 81 dtmp = (Xin(i,1) + Xin(i,2)) * 0.5_dp 82 Xin(i,2) = (Xin(i,1) - Xin(i,2)) * 0.5_dp 83 Xin(i,1) = dtmp 84 end do 85!$OMP end do nowait 86 87!$OMP do 88 do i = 1 , size(Xout,1) 89 dtmp = (Xout(i,1) + Xout(i,2)) * 0.5_dp 90 Xout(i,2) = (Xout(i,1) - Xout(i,2)) * 0.5_dp 91 Xout(i,1) = dtmp 92 end do 93!$OMP end do 94 95 end select 96 97 ! Calculate the residual 98!$OMP do 99 do is = 1 , nspin_mix 100 do i = 1 , size(Xin, 1) 101 F(i,is) = Xout(i,is) - Xin(i,is) 102 end do 103 end do 104!$OMP end do nowait 105 106!$OMP end parallel 107 108 ! Call mixing routine 109 ! Xin contains the input element, F contains 110 ! F = Xout - Xin 111 ! Upon exit Xout contains the mixed quantity 112 select case ( mix_spin ) 113 case ( MIX_SPIN_ALL ) 114 call mixing( scf_mix, maxnh, nspin_mix, Xin, F, Xout) 115 case ( MIX_SPIN_SPINOR, MIX_SPIN_SUM_DIFF ) 116 call mixing( scf_mix, maxnh, nspin_mix, Xin, F, Xout, 117 & nsub=2) 118 case ( MIX_SPIN_SUM ) 119 call mixing( scf_mix, maxnh, nspin_mix, Xin, F, Xout, 120 & nsub=1) 121 end select 122 123!$OMP parallel default(shared), private(dtmp,i,is) 124 125 ! correct output... 126 select case ( mix_spin ) 127 case ( MIX_SPIN_ALL , MIX_SPIN_SPINOR ) 128 ! Xin and Xout are unchanged... 129 case ( MIX_SPIN_SUM , MIX_SPIN_SUM_DIFF ) 130 131 ! transfer spin-sum and difference to 132 ! spin-up and spin-down 133!$OMP do 134 do i = 1 , size(Xin,1) 135 dtmp = Xin(i,1) + Xin(i,2) 136 Xin(i,2) = Xin(i,1) - Xin(i,2) 137 Xin(i,1) = dtmp 138 end do 139!$OMP end do nowait 140 141!$OMP do 142 do i = 1 , size(Xout,1) 143 dtmp = Xout(i,1) + Xout(i,2) 144 Xout(i,2) = Xout(i,1) - Xout(i,2) 145 Xout(i,1) = dtmp 146 end do 147!$OMP end do 148 149 end select 150 151 152 ! Correctly handle mixed quantity 153 ! move over mixed quantity to "input" 154 if ( mix_scf_first ) then 155 156 ! always allow mixing 157 ! I.e. we do nothing 158 159 else if ( iscf == 1 ) then 160 161 ! We are not allowed to mix the first SCF 162 ! Remember that Xout contains the MIXED 163 ! quantity and F is the current residual 164 ! hence for no mixing we re-create Xout: 165 ! Xout = Xin + F 166!$OMP do 167 do is = 1 , nspin_mix 168 do i = 1 , size(Xin, 1) 169 Xout(i,is) = Xin(i,is) + F(i,is) 170 end do 171 end do 172!$OMP end do nowait 173 174 end if 175 176!$OMP end parallel 177 178 deallocate(F) 179 180 ! Print populations at each SCF step, if requested 181 ! Note that this is after mixing, which is not 182 ! entirely correct. It should be moved to the top, 183 ! or done somewhere else. 184 185 if (muldeb .and. (.not. MixH) ) then 186 if (IONode) write (6,"(/a)") 'Using DM after mixing:' 187 if ( spin%SO .and. orbmoms) then 188 call moments( 1, na_u, no_u, maxnh, numh, listhptr, 189 . listh, S, Dscf, isa, lasto, iaorb, iphorb, 190 . indxuo ) 191 endif 192 ! Call this unconditionally 193 call mulliken( mullipop, na_u, no_u, maxnh, 194 & numh, listhptr, listh, S, Dscf, isa, 195 & lasto, iaorb, iphorb ) 196 endif 197 198 call timer( 'MIXER', 2 ) 199 200!-------------------------------------------------------- END 201 END subroutine mixer 202 203 End MODULE m_mixer 204 205 206