TABLE OF CONTENTS


ABINIT/extrapwf [ Functions ]

[ Top ] [ Functions ]

NAME

 extrapwf

FUNCTION

 Extrapolate wavefunctions for new ionic positions
 from values of wavefunctions of previous SCF cycle.
 Use algorithm proposed by T. A.  Arias et al. in PRB 45, 1538 (1992)

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (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

  atindx(natom)=index table for atoms
  atindx1(natom)=index table for atoms, inverse of atindx
  dtset <type(dataset_type)>=all input variables in this dataset
  istep=number of call the routine
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mgfft=maximum size of 1D FFTs
  mpi_enreg=information about MPI parallelization
  nattyp(ntypat)=number of atoms of each type in cell.
  ngfft(18)=contain all needed information about 3D FFT
  npwarr(nkpt)=number of planewaves in basis at this k point
  ntypat=number of types of atoms in cell
  pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
  psps<type(pseudopotential_type)>=variables related to pseudopotentials
  rprimd(3,3)=dimensional primitive translation vectors (bohr)
  usepaw= 0 for non paw calculation; =1 for paw calculation
  xred_old(3,natom)=old reduced coordinates for atoms in unit cell
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point

SIDE EFFECTS

  cg(2,mcg)= plane wave wavefunction coefficient
                          Value from previous SCF cycle is input
                          Extrapolated value is output
  scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles

NOTES

  THIS ROUTINE IS NOT USEABLE AT PRESENT.
  SHOULD BE CAREFULY TESTED AND DEBUGGED (ESPECIALLY WITHIN PAW).

PARENTS

      extraprho

CHILDREN

      ctocprj,dotprod_g,getph,hermit,metric,pawcprj_alloc,pawcprj_copy
      pawcprj_free,pawcprj_get,pawcprj_getdim,pawcprj_lincom,pawcprj_put
      pawcprj_zaxpby,xmpi_allgather,xmpi_alltoallv,zgemm,zhpev

SOURCE

 59 #if defined HAVE_CONFIG_H
 60 #include "config.h"
 61 #endif
 62 
 63 #include "abi_common.h"
 64 
 65 subroutine extrapwf(atindx,atindx1,cg,dtset,istep,kg,mcg,mgfft,mpi_enreg,&
 66 & nattyp,ngfft,npwarr,ntypat,pawtab,psps,rprimd,scf_history,usepaw,xred_old,ylm)
 67 
 68  use defs_basis
 69  use defs_abitypes
 70  use m_scf_history
 71  use m_xmpi
 72  use m_profiling_abi
 73  use m_errors
 74  use m_cgtools
 75 
 76  use m_numeric_tools,   only : hermit
 77  use m_geometry, only : metric
 78  use m_kg, only : getph 
 79  use defs_datatypes, only : pseudopotential_type
 80  use m_pawtab, only : pawtab_type
 81  use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_get, pawcprj_lincom, &
 82 &                      pawcprj_free, pawcprj_zaxpby, pawcprj_put, pawcprj_getdim
 83 
 84 !This section has been created automatically by the script Abilint (TD).
 85 !Do not modify the following lines by hand.
 86 #undef ABI_FUNC
 87 #define ABI_FUNC 'extrapwf'
 88  use interfaces_32_util
 89  use interfaces_66_nonlocal
 90 !End of the abilint section
 91 
 92  implicit none
 93 
 94 !Arguments ------------------------------------
 95 !scalars
 96  integer,intent(in) :: istep,mcg,mgfft,ntypat,usepaw
 97  type(MPI_type),intent(in) :: mpi_enreg
 98  type(dataset_type),intent(in) :: dtset
 99  type(scf_history_type),intent(inout) :: scf_history
100  type(pseudopotential_type),intent(in) :: psps
101 !arrays
102  integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem),nattyp(ntypat),ngfft(18)
103  integer,intent(in) :: npwarr(dtset%nkpt)
104  real(dp),intent(in) :: rprimd(3,3)
105  real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
106  real(dp), intent(inout) :: cg(2,mcg)
107  real(dp),intent(in) :: xred_old(3,dtset%natom)
108  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
109 
110 !Local variables-------------------------------
111 !scalars
112  integer :: ia,iat,iatom,iband_max,iband_max1,iband_min,iband_min1,ibd,ibg,iblockbd,iblockbd1,icg,icgb,icgb1,icgb2
113  integer :: ierr,ig,ii,ikpt,ilmn1,ilmn2,inc,ind1,ind2,iorder_cprj
114  integer :: isize,isppol,istep1,istwf_k,itypat,klmn,me_distrb,my_nspinor
115  integer :: nband_k,nblockbd,nprocband,npw_k,npw_nk,spaceComm_band
116  real(dp) :: dotr,dotr1,doti,doti1,eigval
117  !character(len=500) :: message
118 !arrays
119  real(dp) :: alpha(2),beta(2),gmet(3,3),gprimd(3,3),rmet(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom),ucvol
120  integer,allocatable :: bufsize(:),bufsize_wf(:),bufdisp(:),bufdisp_wf(:),dimcprj(:),npw_block(:),npw_disp(:)
121  real(dp),allocatable :: al(:,:),anm(:),cwavef(:,:),cwavef1(:,:),cwavef_tmp(:,:),deltawf1(:,:),deltawf2(:,:)
122  real(dp),allocatable :: eig(:),evec(:,:)
123  real(dp),allocatable :: unm(:,:,:)
124  real(dp),allocatable :: work(:,:),work1(:,:),wf1(:,:),ylmgr_k(:,:,:),zhpev1(:,:),zhpev2(:)
125  complex(dpc),allocatable :: unm_tmp(:,:),anm_tmp(:,:)
126  type(pawcprj_type),allocatable :: cprj(:,:),cprj_k(:,:),cprj_k1(:,:),cprj_k2(:,:),cprj_k3(:,:),cprj_k4(:,:)
127 !complex(dpc) :: aa
128 
129 ! *************************************************************************
130 
131  if (istep==0) return
132 
133 !Useful array
134  if (usepaw==1) then
135    ABI_ALLOCATE(dimcprj,(dtset%natom))
136    call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O')
137  end if
138 
139 !Metric
140  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
141 
142 !History indexes
143  ind1=1;ind2=2
144 
145 !First step
146  if (istep==1) then
147    scf_history%cg(:,:,ind1)=cg(:,:)
148 !  scf_history%cg(:,:,ind2)=zero
149    scf_history%cg(:,:,ind2)= cg(:,:)
150    if(usepaw==1) then
151 !    WARNING: THIS SECTION IS USELESS; NOW crpj CAN BE READ FROM SCFCV
152      call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred_old)
153      iatom=0 ; iorder_cprj=0
154      call pawcprj_alloc(scf_history%cprj(:,:,ind1),0,dimcprj)
155      call pawcprj_alloc(scf_history%cprj(:,:,ind2),0,dimcprj)
156      ABI_ALLOCATE(ylmgr_k,(dtset%mpw,3,0))
157      call ctocprj(atindx,cg,1,scf_history%cprj(:,:,ind1),gmet,gprimd,&
158 &     iatom,0,iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,scf_history%mcprj,&
159 &     dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,&
160 &     dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,&
161 &     dtset%nloalg,npwarr,dtset%nspinor,dtset%nsppol,dtset%ntypat,&
162 &     dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,0,&
163 &     xred_old,ylm,ylmgr_k)
164      ABI_DEALLOCATE(ylmgr_k)
165 !    call pawcprj_set_zero(scf_history%cprj(:,:,ind2))
166      call pawcprj_copy(scf_history%cprj(:,:,ind1),scf_history%cprj(:,:,ind2))
167    end if
168  else
169 
170 !From 2nd step
171 
172 !  Init parallelism
173    me_distrb=mpi_enreg%me_kpt
174    if (mpi_enreg%paral_kgb==1.or.mpi_enreg%paralbd==1) then
175      spaceComm_band=mpi_enreg%comm_band
176      nprocband=mpi_enreg%nproc_band
177    else
178      spaceComm_band=xmpi_comm_self
179      nprocband=1
180    end if
181 
182 !  For the moment sequential part only
183    nprocband=1
184 
185 !  Additional statements if band-fft parallelism
186    if (nprocband>1) then
187      ABI_ALLOCATE(npw_block,(nprocband))
188      ABI_ALLOCATE(npw_disp,(nprocband))
189      ABI_ALLOCATE(bufsize,(nprocband))
190      ABI_ALLOCATE(bufdisp,(nprocband))
191      ABI_ALLOCATE(bufsize_wf,(nprocband))
192      ABI_ALLOCATE(bufdisp_wf,(nprocband))
193    end if
194 
195    icg=0
196    ibg=0
197 
198    if(usepaw==1) then
199 !    WARNING: THIS SECTION IS USELESS; NOW cprj CAN BE READ FROM SCFCV
200      call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred_old)
201      ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,scf_history%mcprj))
202      call pawcprj_alloc(cprj,0,dimcprj)
203      iatom=0 ; iorder_cprj=0
204      ABI_ALLOCATE(ylmgr_k,(dtset%mpw,3,0))
205      call ctocprj(atindx,cg,1,cprj,gmet,gprimd,iatom,0,iorder_cprj,&
206 &     dtset%istwfk,kg,dtset%kptns,mcg,scf_history%mcprj,dtset%mgfft,&
207 &     dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,&
208 &     nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,dtset%nloalg,&
209 &     npwarr,dtset%nspinor,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,&
210 &     ph1d,psps,rmet,dtset%typat,ucvol,0,xred_old,&
211 &     ylm,ylmgr_k)
212      ABI_DEALLOCATE(ylmgr_k)
213    end if  ! end usepaw=1
214 
215 !  LOOP OVER SPINS
216    my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
217    do isppol=1,dtset%nsppol
218 
219 !    BIG FAT k POINT LOOP
220      do ikpt=1,dtset%nkpt
221 
222 !      Select k point to be treated by this proc
223        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
224        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
225 
226        istwf_k=dtset%istwfk(ikpt)
227 
228 !      Retrieve number of plane waves
229        npw_k=npwarr(ikpt)
230        if (nprocband>1) then
231 !        Special treatment for band-fft //
232          call xmpi_allgather(npw_k,npw_block,spaceComm_band,ierr)
233          npw_nk=sum(npw_block);npw_disp(1)=0
234          do ii=2,nprocband
235            npw_disp(ii)=npw_disp(ii-1)+npw_block(ii-1)
236          end do
237        else
238          npw_nk=npw_k
239        end if
240 
241 !      Allocate arrays for a wave-function (or a block of WFs)
242        ABI_ALLOCATE(cwavef,(2,npw_nk*my_nspinor))
243        ABI_ALLOCATE(cwavef1,(2,npw_nk*my_nspinor))
244        if (nprocband>1) then
245          isize=2*my_nspinor;bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:)
246          isize=2*my_nspinor*npw_k;bufsize_wf(:)=isize
247          do ii=1,nprocband
248            bufdisp_wf(ii)=(ii-1)*isize
249          end do
250        end if
251 
252 !      Subspace alignment
253 
254 !      Loop over bands or blocks of bands
255        nblockbd=nband_k/nprocband
256        icgb=icg
257 
258        if(usepaw==1) then
259          ABI_DATATYPE_ALLOCATE( cprj_k,(dtset%natom,my_nspinor*nblockbd))
260          call pawcprj_alloc(cprj_k,cprj(1,1)%ncpgr,dimcprj)
261          call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,dtset%mband,&
262 &         dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,&
263 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
264          ABI_DATATYPE_ALLOCATE( cprj_k1,(dtset%natom,my_nspinor*nblockbd))
265          call pawcprj_alloc(cprj_k1,scf_history%cprj(1,1,ind1)%ncpgr,dimcprj)
266          call pawcprj_get(atindx1,cprj_k1,scf_history%cprj(:,:,ind1),dtset%natom,1,ibg,ikpt,1,isppol,&
267 &         dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,&
268 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
269          ABI_DATATYPE_ALLOCATE( cprj_k2,(dtset%natom,my_nspinor*nblockbd))
270          call pawcprj_alloc(cprj_k2,scf_history%cprj(1,1,ind2)%ncpgr,dimcprj)
271          call pawcprj_get(atindx1,cprj_k2,scf_history%cprj(:,:,ind2),dtset%natom,1,ibg,ikpt,1,isppol,&
272 &         dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,&
273 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
274        end if  !end usepaw=1
275 
276        ABI_ALLOCATE(unm,(2,nblockbd,nblockbd))
277        unm=zero
278        icgb2=0
279 
280        do iblockbd=1,nblockbd
281          iband_min=1+(iblockbd-1)*nprocband
282          iband_max=iblockbd*nprocband
283 
284          if(xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then
285            if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min,iband_max,isppol,me_distrb)) cycle
286          end if
287 
288 !        Extract wavefunction information
289          if (nprocband>1) then
290 !          Special treatment for band-fft //
291            ABI_ALLOCATE(cwavef_tmp,(2,npw_k*my_nspinor*nprocband))
292            do ig=1,npw_k*my_nspinor*nprocband
293              cwavef_tmp(1,ig)=cg(1,ig+icgb)
294              cwavef_tmp(2,ig)=cg(2,ig+icgb)
295            end do
296            call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef,bufsize,bufdisp,spaceComm_band,ierr)
297            ABI_DEALLOCATE(cwavef_tmp)
298          else
299            do ig=1,npw_k*my_nspinor
300              cwavef(1,ig)=cg(1,ig+icgb)
301              cwavef(2,ig)=cg(2,ig+icgb)
302            end do
303          end if
304 
305          icgb1=icg
306 
307          do iblockbd1=1,nblockbd
308            iband_min1=1+(iblockbd1-1)*nprocband
309            iband_max1=iblockbd1*nprocband
310 
311            if(xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then
312              if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min1,iband_max1,isppol,me_distrb)) cycle
313            end if
314 
315 !          Extract wavefunction information
316 
317            if (nprocband>1) then
318 !            Special treatment for band-fft //
319              ABI_ALLOCATE(cwavef_tmp,(2,npw_k*my_nspinor*nprocband))
320              do ig=1,npw_k*my_nspinor*nprocband
321                cwavef_tmp(1,ig)=scf_history%cg(1,ig+icgb1,ind1)
322                cwavef_tmp(2,ig)=scf_history%cg(2,ig+icgb1,ind1)
323              end do
324              call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef1,bufsize,bufdisp,spaceComm_band,ierr)
325              ABI_DEALLOCATE(cwavef_tmp)
326            else
327              do ig=1,npw_k*my_nspinor
328                cwavef1(1,ig)=scf_history%cg(1,ig+icgb1,ind1)
329                cwavef1(2,ig)=scf_history%cg(2,ig+icgb1,ind1)
330              end do
331            end if
332 
333 !          Calculate Unm=<psi_nk(t)|S|psi_mk(t-dt)>
334            call dotprod_g(dotr,doti,istwf_k,npw_k*my_nspinor,2,cwavef,cwavef1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
335            if(usepaw==1) then
336              ia =0
337              do itypat=1,ntypat
338                do iat=1+ia,nattyp(itypat)+ia
339                  do ilmn1=1,pawtab(itypat)%lmn_size
340                    do ilmn2=1,ilmn1
341                      klmn=((ilmn1-1)*ilmn1)/2+ilmn2
342                      dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)+&
343 &                     cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2))
344                      doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)-&
345 &                     cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2))
346                    end do
347                    do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
348                      klmn=((ilmn2-1)*ilmn2)/2+ilmn1
349                      dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)+&
350 &                     cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2))
351                      doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)-&
352 &                     cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2))
353                    end do
354                  end do
355                end do
356                ia=ia+nattyp(itypat)
357              end do
358            end if
359 !          unm(1,iblockbd,iblockbd1)=dotr
360 !          unm(2,iblockbd,iblockbd1)=doti
361            unm(1,iblockbd1,iblockbd)=dotr
362            unm(2,iblockbd1,iblockbd)=doti
363 !          End loop over bands iblockbd1
364            icgb1=icgb1+npw_k*my_nspinor*nprocband
365 
366          end do
367 
368 !        End loop over bands iblockbd
369          icgb2=icgb2+npw_k*my_nspinor*nprocband
370          icgb=icgb+npw_k*my_nspinor*nprocband
371        end do
372 
373 !      write(std_out,*) 'UNM'
374 !      do iblockbd=1,nblockbd
375 !      write(std_out,11) (unm(1,iblockbd,iblockbd1),unm(2,iblockbd,iblockbd1),iblockbd1=1,nblockbd)
376 !      end do
377 !      11 format(12(1x,f9.5),a)
378 !      Compute A=tU^*U
379        ABI_ALLOCATE(unm_tmp,(nblockbd,nblockbd))
380        ABI_ALLOCATE(anm_tmp,(nblockbd,nblockbd))
381        ABI_ALLOCATE(anm,(nblockbd*(nblockbd+1)))
382        unm_tmp(:,:)=cmplx(unm(1,:,:),unm(2,:,:),kind=dp)
383        call zgemm('C','N',nblockbd,nblockbd,nblockbd,dcmplx(1._dp), unm_tmp,nblockbd, &
384 &       unm_tmp,nblockbd,dcmplx(0._dp),anm_tmp,nblockbd)
385        do iblockbd=1,nblockbd
386          do iblockbd1=iblockbd,nblockbd
387            ii=iblockbd1*(iblockbd1-1)+2*(iblockbd-1)+1
388            anm(ii)=real(anm_tmp(iblockbd,iblockbd1))
389            anm(ii+1)=aimag(anm_tmp(iblockbd,iblockbd1))
390          end do
391        end do
392        call hermit(anm,anm,ierr,nblockbd)
393 !      aa=dcmplx(0._dp)
394 !      do iblockbd=1,nblockbd
395 !      aa=aa+conjg(unm_tmp(iblockbd,1))*unm_tmp(iblockbd,1)
396 !      end do
397 !      write(std_out,*) 'tU*U', aa
398 !      write(std_out,*) 'ANM_tmp'
399 !      do iblockbd=1,nblockbd
400 !      write(std_out,11) (anm_tmp(iblockbd,iblockbd1),iblockbd1=1,nblockbd)
401 !      end do
402 !      write(std_out,*) 'ANM'
403 !      do iblockbd=1,nblockbd*(nblockbd+1)
404 !      write(std_out,11) anm(iblockbd)
405 !      end do
406 
407 !      Diagonalize A
408        ABI_ALLOCATE(eig,(nblockbd))
409        ABI_ALLOCATE(evec,(2*nblockbd,nblockbd))
410        ABI_ALLOCATE(zhpev1,(2,2*nblockbd-1))
411        ABI_ALLOCATE(zhpev2,(3*nblockbd-2))
412        call zhpev('V','U',nblockbd,anm,eig,evec,nblockbd,zhpev1,&
413 &       zhpev2,ierr)
414        ABI_DEALLOCATE(anm)
415        ABI_DEALLOCATE(zhpev1)
416        ABI_DEALLOCATE(zhpev2)
417 !      aa=dcmplx(0._dp)
418 !      do iblockbd=1,nblockbd
419 !      aa=aa+anm_tmp(1,iblockbd)*cmplx(evec((2*iblockbd-1),1),evec(2*iblockbd,1),kind=dp)
420 !      end do
421 !      write(std_out,*) 'EIG', aa, eig(1)*evec(1,1),eig(1)*evec(2,1)
422 
423 !      Compute A'=evec*tU^/sqrt(eig)
424        call zgemm('C','C',nblockbd,nblockbd,nblockbd,dcmplx(1._dp),evec,nblockbd, &
425 &       unm_tmp,nblockbd,dcmplx(0._dp),anm_tmp,nblockbd)
426        do iblockbd=1,nblockbd
427          eigval=dsqrt(eig(iblockbd))
428          do iblockbd1=1,nblockbd
429            anm_tmp(iblockbd,iblockbd1)=anm_tmp(iblockbd,iblockbd1)/eigval
430          end do
431        end do
432 
433 !      Compute tA^A'to come back to the initial subspace for the cg's
434 
435        call zgemm('N','N',nblockbd,nblockbd,nblockbd,dcmplx(1._dp),evec,nblockbd, &
436 &       anm_tmp,nblockbd,dcmplx(0._dp),unm_tmp,nblockbd)
437        anm_tmp=unm_tmp
438 !      write(std_out,*) 'ANM_tmp'
439 !      do iblockbd=1,nblockbd
440 !      write(std_out,11) (anm_tmp(iblockbd,iblockbd1),iblockbd1=1,nblockbd)
441 !      end do
442 
443 !      Wavefunction alignment (istwfk=1 ?)
444        ABI_ALLOCATE(work,(2,npw_nk*my_nspinor*nblockbd))
445        ABI_ALLOCATE(work1,(2,my_nspinor*nblockbd*npw_nk))
446        work1(:,:)=scf_history%cg(:,icg+1:icg+my_nspinor*nblockbd*npw_nk,ind1)
447        call zgemm('N','N',npw_nk*my_nspinor,nblockbd,nblockbd,dcmplx(1._dp), &
448 &       work1,npw_nk*my_nspinor, &
449 &       anm_tmp,nblockbd,dcmplx(0._dp),work,npw_nk*my_nspinor)
450        scf_history%cg(:,1+icg:npw_nk*my_nspinor*nblockbd+icg,ind1)=work(:,:)
451 
452        work1(:,:)=scf_history%cg(:,icg+1:icg+my_nspinor*nblockbd*npw_nk,ind2)
453        call zgemm('N','N',npw_nk*my_nspinor,nblockbd,nblockbd,dcmplx(1._dp), &
454 &       work1,npw_nk*my_nspinor, &
455 &       anm_tmp,nblockbd,dcmplx(0._dp),work,npw_nk*my_nspinor)
456        scf_history%cg(:,1+icg:npw_nk*my_nspinor*nblockbd+icg,ind2)=work(:,:)
457        ABI_DEALLOCATE(work1)
458 !      If paw, must also align cprj:
459        if (usepaw==1) then
460 !        New version (MT):
461          ABI_DATATYPE_ALLOCATE(cprj_k3,(dtset%natom,my_nspinor))
462          ABI_DATATYPE_ALLOCATE(cprj_k4,(dtset%natom,my_nspinor))
463          call pawcprj_alloc(cprj_k3,cprj_k1(1,1)%ncpgr,dimcprj)
464          call pawcprj_alloc(cprj_k4,cprj_k2(1,1)%ncpgr,dimcprj)
465          ABI_ALLOCATE(al,(2,nblockbd))
466          do iblockbd=1,nblockbd
467            ii=(iblockbd-1)*my_nspinor
468            do iblockbd1=1,nblockbd
469              al(1,iblockbd1)=real (anm_tmp(iblockbd,iblockbd1))
470              al(2,iblockbd1)=aimag(anm_tmp(iblockbd,iblockbd1))
471            end do
472            call pawcprj_lincom(al,cprj_k1,cprj_k3,nblockbd)
473            call pawcprj_lincom(al,cprj_k2,cprj_k4,nblockbd)
474            call pawcprj_copy(cprj_k3,cprj_k1(:,ii+1:ii+my_nspinor))
475            call pawcprj_copy(cprj_k4,cprj_k2(:,ii+1:ii+my_nspinor))
476          end do
477          ABI_DEALLOCATE(al)
478 !        Old version (FJ):
479 !        allocate( cprj_k3(dtset%natom,my_nspinor*nblockbd))
480 !        call pawcprj_alloc(cprj_k3,cprj_k1(1,1)%ncpgr,dimcprj)
481 !        allocate( cprj_k4(dtset%natom,my_nspinor*nblockbd))
482 !        call pawcprj_alloc(cprj_k4,cprj_k2(1,1)%ncpgr,dimcprj)
483 !        beta(1)=one;beta(2)=zero
484 !        do iblockbd=1,nblockbd*my_nspinor
485 !        do iblockbd1=1,nblockbd*my_nspinor
486 !        alpha(1)=real(anm_tmp(iblockbd,iblockbd1));alpha(2)=aimag(anm_tmp(iblockbd,iblockbd1))
487 !        call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd1:iblockbd1),cprj_k3(:,iblockbd:iblockbd))
488 !        call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd1:iblockbd1),cprj_k4(:,iblockbd:iblockbd))
489 !        end do
490 !        end do
491 !        call pawcprj_copy(cprj_k3,cprj_k1)
492 !        call pawcprj_copy(cprj_k4,cprj_k2)
493 
494          call pawcprj_free(cprj_k3)
495          call pawcprj_free(cprj_k4)
496          ABI_DATATYPE_DEALLOCATE(cprj_k3)
497          ABI_DATATYPE_DEALLOCATE(cprj_k4)
498        end if
499        ABI_DEALLOCATE(anm_tmp)
500        ABI_DEALLOCATE(unm_tmp)
501        ABI_DEALLOCATE(work)
502 
503 !      Wavefunction extrapolation
504        ibd=0
505        inc=npw_nk*my_nspinor
506        ABI_ALLOCATE(deltawf2,(2,npw_nk*my_nspinor))
507        ABI_ALLOCATE(wf1,(2,npw_nk*my_nspinor))
508        ABI_ALLOCATE(deltawf1,(2,npw_nk*my_nspinor))
509        do iblockbd=1,nblockbd
510          deltawf2(:,:)=scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind2)
511          wf1(:,:)=scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind1)
512 !        wf1(2,1)=zero;deltawf2(2,1)=zero
513 
514          call dotprod_g(dotr,doti,istwf_k,npw_nk*my_nspinor,2,cg(:,icg+1+ibd:ibd+icg+inc),cg(:,icg+1+ibd:ibd+icg+inc),&
515 &         mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
516          call dotprod_g(dotr1,doti1,istwf_k,npw_nk*my_nspinor,2,cg(:,icg+1+ibd:ibd+icg+inc),wf1,&
517 &         mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
518          if(usepaw==1) then
519            ia =0
520            do itypat=1,ntypat
521              do iat=1+ia,nattyp(itypat)+ia
522                do ilmn1=1,pawtab(itypat)%lmn_size
523                  do ilmn2=1,ilmn1
524                    klmn=((ilmn1-1)*ilmn1)/2+ilmn2
525                    dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)+&
526 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2))
527                    doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)-&
528 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2))
529                    dotr1=dotr1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)+&
530 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2))
531                    doti1=doti1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)-&
532 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2))
533                  end do
534                  do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
535                    klmn=((ilmn2-1)*ilmn2)/2+ilmn1
536                    dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)+&
537 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2))
538                    doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)-&
539 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2))
540                    dotr1=dotr1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)+&
541 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2))
542                    doti1=doti1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)-&
543 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2))
544                  end do
545                end do
546              end do
547              ia=ia+nattyp(itypat)
548            end do
549          end if
550          dotr=sqrt(dotr**2+doti**2)
551          dotr1=sqrt(dotr1**2+doti1**2)
552          write(std_out,*)'DOTR, DOTR1',dotr,dotr1
553          dotr=dotr1/dotr
554          write(std_out,*)'DOTR',dotr
555          deltawf1=zero
556          if(dotr>=0.9d0) then
557            deltawf1(:,:)=cg(:,icg+1+ibd:ibd+icg+inc)-wf1(:,:)
558            if(usepaw==1) then
559              alpha(1)=one;alpha(2)=zero
560              beta(1)=-one;beta(2)=zero
561              ia =0
562              call pawcprj_zaxpby(alpha,beta,cprj_k(:,iblockbd:iblockbd),cprj_k1(:,iblockbd:iblockbd))
563            end if
564            istep1=istep
565          else
566            istep1=1
567          end if
568          scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind1)=cg(:,icg+1+ibd:ibd+icg+inc)
569          scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind2)=deltawf1(:,:)
570          if(usepaw==1) then
571            call pawcprj_put(atindx1,cprj_k,scf_history%cprj(:,:,ind1),dtset%natom,1,ibg,ikpt,1,isppol,&
572 &           dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
573 &           mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
574            call pawcprj_put(atindx1,cprj_k1,scf_history%cprj(:,:,ind2),dtset%natom,1,ibg,ikpt,1,isppol,&
575 &           dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
576 &           mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
577          end if
578 
579 !        if(istep1>=3) then
580          cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+scf_history%alpha*deltawf1(:,:) &
581 &         +scf_history%beta *deltawf2(:,:)
582 
583 !        to be used later
584 !        if(usepaw==1) then
585 !        alpha(2)=zero
586 !        beta(1)=one;beta(2)=zero
587 !        alpha(1)=scf_history%alpha
588 !        call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
589 !        alpha(1)=scf_history%beta
590 !        call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
591 !        call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,&
592 !        &    dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
593 !        &    mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
594 !        end if
595 !        else if (istep1==2) then
596 !          cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+scf_history%alpha*deltawf1(:,:)+scf_history%beta*wf1(:,:)
597 !       !     cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+deltawf1(:,:)
598 !        if(usepaw==1) then
599 !        alpha(2)=zero
600 !        beta(1)=one;beta(2)=zero
601 !        alpha(1)=scf_history%alpha
602 !        call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
603 !        alpha(1)=scf_history%beta
604 !        call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
605 !        call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,&
606 !        &    dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
607 !        &    mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
608 !        end if
609 !        end if
610          ibd=ibd+inc
611        end do ! end loop on iblockbd
612 
613        ABI_DEALLOCATE(deltawf1)
614        ABI_DEALLOCATE(deltawf2)
615        ABI_DEALLOCATE(wf1)
616        ABI_DEALLOCATE(cwavef)
617        ABI_DEALLOCATE(cwavef1)
618        ABI_DEALLOCATE(eig)
619        ABI_DEALLOCATE(evec)
620        ABI_DEALLOCATE(unm)
621        if(usepaw==1) then
622          call pawcprj_free(cprj_k)
623          ABI_DATATYPE_DEALLOCATE(cprj_k)
624          call pawcprj_free(cprj_k1)
625          ABI_DATATYPE_DEALLOCATE(cprj_k1)
626          call pawcprj_free(cprj_k2)
627          ABI_DATATYPE_DEALLOCATE(cprj_k2)
628        end if
629 
630        ibg=ibg+my_nspinor*nband_k
631        icg=icg+my_nspinor*nband_k*npw_k
632 
633 !      End big k point loop
634      end do
635 !    End loop over spins
636    end do
637 
638    if(usepaw==1) then
639      call pawcprj_free(cprj)
640      ABI_DATATYPE_DEALLOCATE(cprj)
641    end if
642    if (nprocband>1) then
643      ABI_DEALLOCATE(npw_block)
644      ABI_DEALLOCATE(npw_disp)
645      ABI_DEALLOCATE(bufsize)
646      ABI_DEALLOCATE(bufdisp)
647      ABI_DEALLOCATE(bufsize_wf)
648      ABI_DEALLOCATE(bufdisp_wf)
649    end if
650 
651  end if ! istep>=2
652 
653  if (usepaw==1) then
654    ABI_DEALLOCATE(dimcprj)
655  end if
656 
657 end subroutine extrapwf