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