TABLE OF CONTENTS


ABINIT/wf_mixing [ Functions ]

[ Top ] [ Functions ]

NAME

 wf_mixing

FUNCTION

 Mixing of wavefunctions in the outer loop of a double loop SCF approach.
 Different algorithms are implemented, depending on the value of wfmixalg.

COPYRIGHT

 Copyright (C) 2017-2018 ABINIT group (XG,MT,FJ)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  atindx1(dtset%natom)=index table for atoms, inverse of atindx
  dtset <type(dataset_type)>=all input variables in this dataset
  istep=number of call the routine (usually the outer loop in the SCF double loop)
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of cprj array 
  mpi_enreg=information about MPI parallelization
  nattyp(dtset%ntypat)=number of atoms of each type in cell.
  npwarr(nkpt)=number of planewaves in basis at this k point
  pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data

SIDE EFFECTS

  cg(2,mcg)= plane wave wavefunction coefficient
                          Value from previous SCF cycle is input and stored in some form
                          Extrapolated value is output
  cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors
                          Value from previous SCF cycle is input and stored in some form
                          Extrapolated value is output
  scf_history_wf <type(scf_history_type)>=arrays obtained from previous SCF cycles

PARENTS

      scfcv

CHILDREN

      cgcprj_cholesky,dotprod_set_cgcprj,dotprodm_sumdiag_cgcprj
      lincom_cgcprj,pawcprj_alloc,pawcprj_axpby,pawcprj_free,pawcprj_get
      pawcprj_getdim,pawcprj_lincom,pawcprj_put,timab,xmpi_sum,zgesv

SOURCE

 48 #if defined HAVE_CONFIG_H
 49 #include "config.h"
 50 #endif
 51 
 52 #include "abi_common.h"
 53 
 54 subroutine wf_mixing(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,&
 55 & nattyp,npwarr,pawtab,scf_history_wf)
 56 
 57  use defs_basis
 58  use defs_abitypes
 59  use m_scf_history
 60  use m_xmpi
 61  use m_profiling_abi
 62  use m_errors
 63  use m_cgtools
 64 
 65  use m_pawtab, only : pawtab_type
 66  use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_get, pawcprj_lincom, &
 67 &                      pawcprj_free, pawcprj_axpby, pawcprj_put, pawcprj_getdim
 68 
 69 !This section has been created automatically by the script Abilint (TD).
 70 !Do not modify the following lines by hand.
 71 #undef ABI_FUNC
 72 #define ABI_FUNC 'wf_mixing'
 73  use interfaces_18_timing
 74  use interfaces_32_util
 75  use interfaces_66_wfs
 76 !End of the abilint section
 77 
 78  implicit none
 79 
 80 !Arguments ------------------------------------
 81 !scalars
 82  integer,intent(in) :: istep,mcg,mcprj
 83  type(MPI_type),intent(in) :: mpi_enreg
 84  type(dataset_type),intent(in) :: dtset
 85  type(scf_history_type),intent(inout) :: scf_history_wf
 86 !arrays
 87  integer,intent(in) :: atindx1(dtset%natom),nattyp(dtset%ntypat)
 88  integer,intent(in) :: npwarr(dtset%nkpt)
 89  real(dp), intent(inout) :: cg(2,mcg)
 90  type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj)
 91  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
 92 
 93 !Local variables-------------------------------
 94 !scalars
 95  integer :: hermitian
 96  integer :: ibdmix,ibdsp,ibg,ibg_hist,icg,icg_hist
 97  integer :: ierr,ikpt,indh,ind_biorthog,ind_biorthog_eff,ind_newwf,ind_residual,inplace
 98  integer :: iset2,isppol,istep_cycle,istep_new,istwf_k,kk,me_distrb,my_nspinor
 99  integer :: nband_k,nbdmix,npw_k,nset1,nset2,ntypat
100  integer :: shift_set1,shift_set2,spaceComm_band,spare_mem,usepaw,wfmixalg
101  real(dp) :: alpha,beta
102  complex(dpc) :: sum_coeffs
103 !arrays
104  integer,allocatable :: ipiv(:),dimcprj(:)
105  real(dp) :: tsec(2)
106  real(dp),allocatable :: al(:,:),mmn(:,:,:)
107  real(dp),allocatable :: dotprod_res(:,:,:),dotprod_res_k(:,:,:),res_mn(:,:,:),smn(:,:,:)
108  complex(dpc),allocatable :: coeffs(:)
109  type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kh(:,:)
110 
111 ! *************************************************************************
112 
113 !DEBUG
114 !write(std_out,*)
115 !write(std_out,*)' wf_mixing : enter, istep= ',istep
116 !call flush(std_out)
117 !write(std_out,*)' istep,scf_history_wf%alpha=',istep,scf_history_wf%alpha
118 !write(std_out,*)' cg(1,1)=',cg(1,1)
119 !write(std_out,*)' scf_history_wf%cg(1,1,1:5)=',scf_history_wf%cg(1,1,1:5)
120 !ABI_ALLOCATE(cg_ref,(2,mcg))
121 !cg_ref(:,:)=cg(:,:)
122 !ABI_DATATYPE_ALLOCATE(cprj_ref,(dtset%natom,mcprj))
123 !cprj_ref(:,:)=cprj(:,:)
124 !      write(std_out,*)' scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)=',&
125 !&       scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)
126 !       call flush(std_out)
127 !ENDDEBUG
128 
129  if (istep==0) return
130 
131  ntypat=dtset%ntypat
132  usepaw=dtset%usepaw
133  wfmixalg=scf_history_wf%wfmixalg
134  nbdmix=dtset%nbandhf
135  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
136  me_distrb=mpi_enreg%me_kpt
137  spaceComm_band=xmpi_comm_self
138 
139  spare_mem=0
140  if(scf_history_wf%history_size==wfmixalg-1)spare_mem=1
141 
142 !scf_history_wf%alpha contains dtset%wfmix
143  alpha=scf_history_wf%alpha
144  beta=one-scf_history_wf%alpha
145 
146  icg=0
147  icg_hist=0
148  ibg=0
149  ibg_hist=0
150 
151 !Useful array
152  ABI_ALLOCATE(dimcprj,(dtset%natom))
153  if (usepaw==1) then
154    call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O')
155  end if
156 
157  if(istep==1)then
158    do indh=1,scf_history_wf%history_size
159      call pawcprj_alloc(scf_history_wf%cprj(:,:,indh),0,dimcprj)
160    end do
161  end if
162 
163  ABI_DATATYPE_ALLOCATE(cprj_k,(dtset%natom,my_nspinor*nbdmix))
164  ABI_DATATYPE_ALLOCATE(cprj_kh,(dtset%natom,my_nspinor*nbdmix))
165  if(usepaw==1) then
166    call pawcprj_alloc(cprj_k,0,dimcprj)
167    call pawcprj_alloc(cprj_kh,0,dimcprj)
168  end if
169  ABI_ALLOCATE(smn,(2,nbdmix,nbdmix))
170  ABI_ALLOCATE(mmn,(2,nbdmix,nbdmix))
171 
172  if(wfmixalg>2)then
173    nset1=1 
174    nset2=min(istep-1,wfmixalg-1) 
175    ABI_ALLOCATE(dotprod_res_k,(2,1,nset2))
176    ABI_ALLOCATE(dotprod_res,(2,1,nset2))
177    ABI_ALLOCATE(res_mn,(2,wfmixalg-1,wfmixalg-1)) 
178    dotprod_res=zero
179    if(istep==1)then
180      scf_history_wf%dotprod_sumdiag_cgcprj_ij=zero
181    end if
182  end if
183 
184 !Explanation for the index for the wavefunction stored in scf_history_wf
185 !The reference is the cg+cprj output after the wf optimization at istep 1. 
186 !It comes as input to the present routine as cgcprj input at step 2, and is usually found at indh=1.
187 
188 !In the simple mixing case (wfmixalg==1), the reference is never stored, because it is used "on-the-fly" to biothogonalize the
189 !previous input (that was stored in indh=1), then generate the next input, which is stored again in indh=1
190 
191 !When the storage is not spared: 
192 !- the values of indh from 2 to wfmixalg store the (computed here) biorthogonalized input cgcprj, then the residual
193 !- the values of indh from wfmixalg+1 to 2*wfmixalg-1 store the biorthogonalized output cgcprj (coming as argument)
194 
195 !First step
196  if (istep==1 .or. (wfmixalg==2 .and. abs(scf_history_wf%alpha-one)<tol8) ) then
197 
198    indh=2   ! This input wavefunction is NOT the reference
199    if(wfmixalg==2)indh=1 ! But this does not matter in the simple mixing case that has history_size=1
200 
201 !  Simply store the wavefunctions and cprj. However, nband_k might be different from nbandhf...
202 !  LOOP OVER SPINS
203    do isppol=1,dtset%nsppol
204 
205 !    BIG FAT k POINT LOOP
206      do ikpt=1,dtset%nkpt
207 
208 !      Select k point to be treated by this proc
209        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
210        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
211 
212        npw_k=npwarr(ikpt)
213 
214        scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
215        if(usepaw==1) then
216 !        scf_history_wf%cprj(:,ibg_hist+1:ibg_hist+my_nspinor*nbdmix,1)=cprj(:,ibg+1:ibg+my_nspinor*nbdmix)
217          call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,dtset%mband,&
218 &         dtset%mkmem,dtset%natom,nbdmix,nband_k,my_nspinor,dtset%nsppol,0,&
219 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
220          call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,1,isppol,&
221 &         nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
222 &         mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
223        end if
224 
225 !      Update the counters
226        ibg=ibg+my_nspinor*nband_k
227        ibg_hist=ibg_hist+my_nspinor*nbdmix
228        icg=icg+my_nspinor*nband_k*npw_k
229        icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
230 
231      end do
232    end do
233 
234  else
235 !  From istep==2
236 
237 !  First part of the computation : biorthogonalization, and computation of the residual (possibly, prediction of the next input in the case of simple mixing)
238 !  Index for the wavefunctions stored in scf_history_wf whose scalar products with the argument cgcprj will have to be computed.
239    indh=1   ! This input wavefunction is the reference
240    if(wfmixalg/=2 .and. istep==2)indh=2 ! except for istep=2 in the rmm-diis
241 
242    if(wfmixalg>2)then
243 !    istep inside the cycle defined by wfmixalg, and next index. Then, indices of the wavefunction sets.
244      istep_cycle=mod((istep-2),wfmixalg-1)
245      istep_new=mod((istep-1),wfmixalg-1)
246      ind_biorthog=1+wfmixalg+istep_cycle
247      ind_residual=2+istep_cycle
248      ind_newwf=2+istep_new
249      shift_set1=ind_residual-1
250      shift_set2=1
251    end if
252 
253 !  LOOP OVER SPINS
254    do isppol=1,dtset%nsppol
255 
256 !    BIG FAT k POINT LOOP
257      do ikpt=1,dtset%nkpt
258 
259 !      Select k point to be treated by this proc
260        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
261        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
262 
263        istwf_k=dtset%istwfk(ikpt)
264        npw_k=npwarr(ikpt)
265 
266 !      Biorthogonalization
267 
268        if(usepaw==1) then
269          call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,dtset%mband,&
270 &         dtset%mkmem,dtset%natom,nbdmix,nband_k,my_nspinor,dtset%nsppol,0,&
271 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
272          call pawcprj_get(atindx1,cprj_kh,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,1,isppol,&
273 &         nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,my_nspinor,dtset%nsppol,0,&
274 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
275        end if  !end usepaw=1
276 
277        hermitian=0
278        if(wfmixalg==2 .or. istep==2)then
279          call dotprod_set_cgcprj(atindx1,cg,scf_history_wf%cg(:,:,indh),cprj_k,cprj_kh,dimcprj,hermitian,&
280 &         0,0,icg,icg_hist,ikpt,isppol,istwf_k,nbdmix,mcg,mcg,mcprj,mcprj,dtset%mkmem,&
281 &         mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn,usepaw)
282        else
283          call dotprod_set_cgcprj(atindx1,scf_history_wf%cg(:,:,indh),cg,cprj_kh,cprj_k,dimcprj,hermitian,&
284 &         0,0,icg,icg_hist,ikpt,isppol,istwf_k,nbdmix,mcg,mcg,mcprj,mcprj,dtset%mkmem,&
285 &         mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn,usepaw)
286        end if
287 
288 !      Invert S matrix, that is NOT hermitian. 
289 !      Calculate M=S^-1
290        mmn=zero
291        do kk=1,nbdmix
292          mmn(1,kk,kk)=one
293        end do
294 
295        ABI_ALLOCATE(ipiv,(nbdmix))
296 !      The smn is destroyed by the following inverse call
297        call zgesv(nbdmix,nbdmix,smn,nbdmix,ipiv,mmn,nbdmix,ierr)
298        ABI_DEALLOCATE(ipiv)
299 !DEBUG
300        if(ierr/=0)then
301          MSG_ERROR(' The call to cgesv general inversion routine failed')
302        end if
303 !ENDDEBUG
304 
305 !      The M matrix is used to compute the biorthogonalized set of wavefunctions, and to store it at the proper place
306        if(wfmixalg==2 .or. istep==2)then
307          inplace=1
308          call lincom_cgcprj(mmn,scf_history_wf%cg(:,:,indh),cprj_kh,dimcprj,&
309 &         icg_hist,inplace,mcg,my_nspinor*nbdmix,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw)
310        else
311          inplace=0
312          call lincom_cgcprj(mmn,cg,cprj_k,dimcprj,&
313 &         icg,inplace,mcg,my_nspinor*nbdmix,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw,&
314 &         cgout=scf_history_wf%cg(:,:,ind_biorthog),cprjout=scf_history_wf%cprj(:,:,ind_biorthog),icgout=icg_hist)
315        end if
316 
317 !      The biorthogonalised set of wavefunctions is now stored at the proper place
318 
319 !      Finalize this first part of the computation, depending on the algorithm and the step.
320 
321        if(wfmixalg==2)then
322 
323 !        Wavefunction extrapolation, simple mixing case
324 !        alpha contains dtset%wfmix, beta contains one-alpha
325          cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=&
326 &         alpha*cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)&
327 &         +beta*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh) 
328          if(usepaw==1) then
329            do ibdmix=1,nbdmix
330              call pawcprj_axpby(beta,alpha,cprj_kh(:,ibdmix:ibdmix),cprj_k(:,ibdmix:ibdmix))
331            end do ! end loop on ibdmix
332          end if
333 
334 !        Back to usual orthonormalization 
335          call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,my_nspinor*nband_k,dtset%mkmem,&
336 &         mpi_enreg,dtset%natom,nattyp,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw)
337 
338 !        Store the newly extrapolated wavefunctions, orthonormalized, in scf_history_wf
339          scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
340          if(usepaw==1) then
341            do ibdmix=1,nbdmix
342              call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,1,isppol,&
343 &             nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
344 &             mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
345            end do ! end loop on ibdmix
346          end if
347 
348        else  !  wfmixalg/=2
349 !        RMM-DIIS
350 
351          if (istep==2)then
352 !          Store the argument wf as the reference for all future steps, in scf_history_wf with index 1. 
353            scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,1)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
354            if(usepaw==1) then
355              do ibdmix=1,nbdmix
356                call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,1),dtset%natom,1,ibg_hist,ikpt,1,isppol,&
357 &               nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
358 &               mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
359              end do ! end loop on ibdmix
360            end if
361          end if
362 
363          ind_biorthog_eff=ind_biorthog
364          if(istep==2)ind_biorthog_eff=1 ! The argument wf has not been stored in ind_biorthog
365 
366 !        Compute the residual of the wavefunctions for this istep, 
367 !        that replaces the previously stored set of (biorthogonalized) input wavefunctions
368          scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual)=&
369 &         scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog_eff)&
370 &         -scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual)
371          if(usepaw==1) then
372            do ibdmix=1,nbdmix
373              call pawcprj_axpby(one,-one,scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog_eff),&
374 &             scf_history_wf%cprj(:,ibdmix:ibdmix,ind_residual))
375            end do ! end loop on ibdmix
376          end if
377          
378 !        Compute the new scalar products to fill the res_mn matrix
379          call dotprodm_sumdiag_cgcprj(atindx1,scf_history_wf%cg,scf_history_wf%cprj,dimcprj,&
380 &         ibg,icg,ikpt,isppol,istwf_k,nbdmix,mcg,mcprj,dtset%mkmem,&
381 &         mpi_enreg,scf_history_wf%history_size,dtset%natom,nattyp,nbdmix,npw_k,nset1,nset2,my_nspinor,dtset%nsppol,ntypat,&
382 &         shift_set1,shift_set2,pawtab,dotprod_res_k,usepaw)
383 
384          dotprod_res=dotprod_res+dotprod_res_k
385 
386 !        scf_history_wf for index ind_biorthog will contain the extrapolated wavefunctions (and no more the output of the SCF loop).
387          scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog)=&
388 &         scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog_eff)+&
389 &         (alpha-one)*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual)
390          if(usepaw==1) then
391            do ibdmix=1,nbdmix
392              if(ind_biorthog/=ind_biorthog_eff)then
393                scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog)=scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog_eff) 
394              end if
395              call pawcprj_axpby((alpha-one),one,scf_history_wf%cprj(:,ibdmix:ibdmix,ind_residual),&
396 &             scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog))
397            end do ! end loop on ibdmix
398          end if
399 
400        end if
401 
402        ibg=ibg+my_nspinor*nband_k
403        ibg_hist=ibg_hist+my_nspinor*nbdmix
404        icg=icg+my_nspinor*nband_k*npw_k
405        icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
406 
407 !      End big k point loop
408      end do
409 !    End loop over spins
410    end do
411 
412  end if ! istep>=2
413 
414  if(wfmixalg>2 .and. istep>1)then
415 
416 !DEBUG
417 !      write(std_out,*)' '
418 !      write(std_out,*)' Entering the residual minimisation part '
419 !      write(std_out,*)' '
420 !      call flush(std_out)
421 !ENDDEBUG
422 
423    call timab(48,1,tsec)
424    call xmpi_sum(dotprod_res,mpi_enreg%comm_kpt,ierr)
425    call timab(48,2,tsec)
426 
427    scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,1+shift_set1,1+shift_set2:nset2+shift_set2)=dotprod_res(:,1,1:nset2)
428    scf_history_wf%dotprod_sumdiag_cgcprj_ij(1,1+shift_set2:nset2+shift_set2,1+shift_set1)=dotprod_res(1,1,1:nset2)
429    scf_history_wf%dotprod_sumdiag_cgcprj_ij(2,1+shift_set2:nset2+shift_set2,1+shift_set1)=-dotprod_res(2,1,1:nset2)
430 
431  end if ! wfmixalg>2 and istep>1
432 
433  if(wfmixalg>2 .and. istep>2)then
434 
435 !  Extract the relevant matrix R_mn
436    res_mn(:,1:nset2,1:nset2)=&
437 &   scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,1+shift_set2:nset2+shift_set2,1+shift_set2:nset2+shift_set2)
438 
439 !DEBUG
440 !      write(std_out,*)' The matrix res_mn(:,1:nset2,1:nset2) is :'
441 !      write(std_out,*)res_mn(:,1:nset2,1:nset2)
442 !      call flush(std_out)
443 !ENDDEBUG
444 
445 !  Solve R_mn \alpha_n = 1_m
446    ABI_ALLOCATE(ipiv,(nset2))
447    ABI_ALLOCATE(coeffs,(nset2))
448    coeffs(:)=cone
449 !  The res_mn is destroyed by the following inverse call
450    call zgesv(nset2,1,res_mn,wfmixalg-1,ipiv,coeffs,nset2,ierr)
451    ABI_DEALLOCATE(ipiv)
452 !  The coefficients must sum to one
453    sum_coeffs=sum(coeffs)
454    coeffs=coeffs/sum_coeffs
455 
456 !DEBUG
457 !      write(std_out,*)' The coefficients that minimize the residual have been found'
458 !      write(std_out,*)' coeffs =',coeffs
459 !      call flush(std_out)
460 !ENDDEBUG
461  end if ! wfmixalg>2 and istep>2
462 
463  if(wfmixalg>2 .and. istep>1)then
464 
465 !  Find the new "input" wavefunction, bi-orthogonalized, and store it replacing the adequate "old" input wavefunction.
466 
467    icg=0
468    icg_hist=0
469    ibg=0
470    ibg_hist=0
471    ABI_ALLOCATE(al,(2,nset2))
472    if(istep>2)then
473      do iset2=1,nset2
474        al(1,iset2)=real(coeffs(iset2)) ; al(2,iset2)=aimag(coeffs(iset2))
475      end do
476    else
477      al(1,1)=one ; al(2,1)=zero
478    end if
479 
480 !DEBUG
481 !      write(std_out,*)' Overload the coefficients, in order to simulate a simple mixing with wfmix '
482 !      write(std_out,*)' Set al(1,ind_biorthog-3)=one, for ind_biorthog=',ind_biorthog
483 !      write(std_out,*)' This will feed scf_history for set ind_biorthog-3+wfmixalg=',ind_biorthog-3+wfmixalg
484 !      al(:,:)=zero
485 !      al(1,ind_biorthog-3)=one
486 !      call flush(std_out)
487 !ENDDEBUG
488 
489 !  LOOP OVER SPINS
490    do isppol=1,dtset%nsppol
491 
492 !    BIG FAT k POINT LOOP
493      do ikpt=1,dtset%nkpt
494 
495 !      Select k point to be treated by this proc
496        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
497        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
498 
499        istwf_k=dtset%istwfk(ikpt)
500        npw_k=npwarr(ikpt)
501 
502        if(istep>2)then
503 !        Make the appropriate linear combination (from the extrapolated wfs)
504          cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=zero
505          do iset2=1,nset2
506            cg(1,icg+1:icg+my_nspinor*npw_k*nband_k)=cg(1,icg+1:icg+my_nspinor*npw_k*nband_k)&
507 &           +al(1,iset2)*scf_history_wf%cg(1,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)&
508 &           -al(2,iset2)*scf_history_wf%cg(2,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)
509            cg(2,icg+1:icg+my_nspinor*npw_k*nband_k)=cg(2,icg+1:icg+my_nspinor*npw_k*nband_k)&
510 &           +al(1,iset2)*scf_history_wf%cg(2,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)&        
511 &           +al(2,iset2)*scf_history_wf%cg(1,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)
512          end do
513        else ! One needs a simple copy from the extrapolated wavefunctions
514          cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,1+wfmixalg)
515        end if
516 !      Note the storage in cprj_k. By the way, a simple copy might also be used in case istep=2.
517        if(usepaw==1) then
518          do ibdsp=1,my_nspinor*nbdmix
519            call pawcprj_lincom(al,scf_history_wf%cprj(:,ibdsp,1+wfmixalg:nset2+wfmixalg),cprj_k(:,ibdsp:ibdsp),nset2)
520          end do
521        end if
522 
523 !      Store the newly extrapolated wavefunctions for this k point, still bi-orthonormalized, in scf_history_wf
524        scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_newwf)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
525        if(usepaw==1) then
526          call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,ind_newwf),dtset%natom,1,ibg_hist,ikpt,1,isppol,&
527 &         nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
528 &         mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
529        end if
530 
531 !      Back to usual orthonormalization for the cg and cprj_k
532        call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,my_nspinor*nband_k,dtset%mkmem,&
533 &       mpi_enreg,dtset%natom,nattyp,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw)
534 
535 !      Need to transfer cprj_k to cprj
536        if(usepaw==1) then
537          call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg_hist,ikpt,1,isppol,&
538 &         nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
539 &         mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
540        end if
541 
542        ibg=ibg+my_nspinor*nband_k
543        ibg_hist=ibg_hist+my_nspinor*nbdmix
544        icg=icg+my_nspinor*nband_k*npw_k
545        icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
546 
547 !      End big k point loop
548      end do
549 !    End loop over spins
550    end do
551 
552    if(istep>2)then
553      ABI_DEALLOCATE(coeffs)
554    end if
555    ABI_DEALLOCATE(al)
556 
557  end if ! wfmixalg>2 and istep>1
558 
559 !DEBUG
560 ! write(std_out,*)' wf_mixing : exit '
561 !      write(std_out,*)' scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)=',&
562 !&       scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)
563 ! write(std_out,*)' cg(1:2,1:2)=',cg(1:2,1:2)
564 ! write(std_out,*)' scf_history_wf%cg(1:2,1:2,1)=',scf_history_wf%cg(1:2,1:2,1)
565 ! ABI_DEALLOCATE(cg_ref)
566 ! ABI_DATATYPE_DEALLOCATE(cprj_ref)
567 !ENDDEBUG
568 
569  if(usepaw==1) then
570    call pawcprj_free(cprj_k)
571    call pawcprj_free(cprj_kh)
572  end if
573  ABI_DATATYPE_DEALLOCATE(cprj_k)
574  ABI_DATATYPE_DEALLOCATE(cprj_kh)
575  ABI_DEALLOCATE(dimcprj)
576  ABI_DEALLOCATE(mmn)
577  ABI_DEALLOCATE(smn)
578  if(wfmixalg>2)then
579    ABI_DEALLOCATE(dotprod_res_k)
580    ABI_DEALLOCATE(dotprod_res)
581    ABI_DEALLOCATE(res_mn)
582  end if
583 end subroutine wf_mixing