TABLE OF CONTENTS


ABINIT/extraprho [ Functions ]

[ Top ] [ Functions ]

NAME

 extraprho

FUNCTION

 Extrapolate electronic density for new ionic positions
 from values of density of previous SCF cycle.
 Use algorithm proposed by D. Alfe in Comp. Phys. Comm. 118 (1999), 31-33 [[cite:Alfe1999]]

INPUTS

  atindx
  atindx1(natom)=index table for atoms, inverse of atindx
  cg(2,mcg)= plane wave wavefunction coefficient
  cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn>
  dtset <type(dataset_type)>=all input variables in this dataset
   | densty(ntypat,4)=parameters for initialisation of the gaussian density
   | jellslab,slabzbeg,slabzend,slabwsrad=parameters for jellium slab
   | natom=number of atoms in cell.
   | nspden=number of spin-density components
  gmet(3,3)=reciprocal space metric
  gprimd(3,3)=reciprocal space dimensional primitive translations
  gsqcut=cutoff value on G**2 for sphere inside fft box
  istep=number of call the routine
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  mgfft=maximum size of 1D FFTs
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mpi_enreg=information about MPI parallelization
  mqgrid=number of grid pts in q array for f(q) spline.
  my_natom=number of atoms treated by current processor
  nattyp(ntypat)=number of atoms of each type in cell.
  nfft=(effective) number of FFT grid points (for this processor)
  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
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
  psps<type(pseudopotential_type)>=variables related to pseudopotentials
  qgrid(mqgrid)=q grid for spline from 0 to qmax
  rprimd(3,3)=dimensional primitive translation vectors (bohr)
  ucvol=unit cell volume (bohr**3).
  usepaw= 0 for non paw calculation; =1 for paw calculation
  xred_new(3,natom)=new reduced coordinates for atoms in unit cell
  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
  zion(ntypat)=charge on each type of atom
  znucl(ntypat)=atomic number of each atom type

SIDE EFFECTS

  pawrhoij(my_natom) <type(pawrhoij_type)>= PAW rhoij occupancies and related data
                                            Value from previous SCF cycle is input
                                            Extrapolated value is output
  rhor(nfft,nspden)=the density from previous SCF cycle is input
                    the extrapolated density is output
  scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles

SOURCE

116 subroutine extraprho(atindx,atindx1,cg,cprj,dtset,gmet,gprimd,gsqcut,istep,&
117 & kg,mcg,mcprj,mgfft,mpi_enreg,mqgrid,my_natom,nattyp,nfft,ngfft,npwarr,ntypat,pawrhoij,&
118 & pawtab,ph1d,psps,qgrid,rhor,rprimd,scf_history,ucvol,usepaw,&
119 & xred_new,xred_old,ylm,zion,znucl)
120 
121 !Arguments ------------------------------------
122 !scalars
123  integer,intent(in) :: istep,mcg,mcprj,mgfft,my_natom,mqgrid,nfft,ntypat,usepaw
124  real(dp),intent(in) :: gsqcut,ucvol
125  type(MPI_type),intent(in) :: mpi_enreg
126  type(dataset_type),intent(in) :: dtset
127  type(scf_history_type),intent(inout) :: scf_history
128  type(pseudopotential_type),intent(in) :: psps
129 !arrays
130  integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem)
131  integer,intent(in) :: nattyp(ntypat),ngfft(18),npwarr(dtset%nkpt)
132  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom)
133  real(dp),intent(in) :: qgrid(mqgrid)
134  real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
135  real(dp),intent(in) ::  zion(ntypat),znucl(ntypat)
136  real(dp), intent(inout) :: cg(2,mcg)
137  real(dp),intent(inout) :: rhor(nfft,dtset%nspden),rprimd(3,3),xred_new(3,dtset%natom)
138  real(dp),intent(in) :: xred_old(3,dtset%natom)
139  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw)
140  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
141  type(pawcprj_type),intent(inout) :: cprj(:,:)
142 
143 !Local variables-------------------------------
144 !scalars
145  integer :: cplex_rhoij,dplex,iatom,ii,ind1,ind1new,ind2,ind2new,iq,iq0,irhoij,ispden,itypat,jrhoij,klmn
146  integer :: lmn2_size,nselect,nspden_rhoij,optatm,optdyfr,opteltfr,optgr,option,optn,optn2
147  integer :: optstr,optv,qphase_rhoij
148  real(dp) :: a11,a12,a22,a33,alpha,b1,b2,beta,detA,fact,ratio1,ratio2
149  logical :: hasmoved,usegauss
150  character(len=500) :: message
151 !arrays
152  integer :: dummy3(3)
153  real(dp) :: diff_t(3),diff_tmdt(3),diff_tpdt(3),dummy2(2)
154  real(dp) :: dummy_in(0)
155  real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0)
156  real(dp) :: strn_dummy6(6),strv_dummy6(6)
157  real(dp),allocatable :: deltarho(:),gauss(:,:),rhoijtmp(:,:),work1(:)
158  real(dp),allocatable :: work2(:,:),work3(:,:),xred_tpdt(:,:)
159 
160 ! *************************************************************************
161 
162 !---------------------------------------------------------------
163 !----------- Inits
164 !---------------------------------------------------------------
165 
166 !History indexes
167  ind1=scf_history%hindex(1)
168  ind2=scf_history%hindex(2)
169 
170 !Compatibility tests
171  if (ind1==0.and.ind2>0)then
172    ABI_BUG(' Incompatible history indexes !')
173  end if
174 
175 !Rotated values of history indexes
176  if (ind1>0.and.ind2>0) then
177    ind1new=ind2;ind2new=ind1
178  else if (ind1>0.and.ind2==0) then
179    ind1new=3-ind1;ind2new=ind1
180  else if (ind1==0.and.ind2==0) then
181    ind1new=1;ind2new=0
182  end if
183 
184 !Compute ionic positions at t+dt in red. coordinates
185 !Has to take the boundary conditions into account
186  ABI_MALLOC(xred_tpdt,(3,dtset%natom))
187  do iatom=1,dtset%natom
188    xred_tpdt(1,iatom)=xred_old(1,iatom)+mod(xred_new(1,iatom)-xred_old(1,iatom)+1.5_dp,one)-half
189    xred_tpdt(2,iatom)=xred_old(2,iatom)+mod(xred_new(2,iatom)-xred_old(2,iatom)+1.5_dp,one)-half
190    xred_tpdt(3,iatom)=xred_old(3,iatom)+mod(xred_new(3,iatom)-xred_old(3,iatom)+1.5_dp,one)-half
191  end do
192 
193 !---------------------------------------------------------------
194 !----------- Compute Alpha and Beta
195 !----------- see (4) in Comp. Phys. Comm. 118 (1999), 31-33 [[cite:Alfe1999]]
196 !---------------------------------------------------------------
197 
198 !Compute a_ij matrix
199  a11=zero;a12=zero;a22=zero;a33=zero;b1=zero;b2=zero
200  diff_t=zero;diff_tmdt=zero;diff_tpdt=zero
201  do iatom=1,dtset%natom
202 
203    diff_tpdt(1:3)=xred_tpdt(1:3,iatom)-xred_old(1:3,iatom)
204    if (ind1>0) then
205      diff_t(1:3)=scf_history%xreddiff(1:3,iatom,ind1)
206      if (ind2>0) diff_tmdt(1:3)=scf_history%xreddiff(1:3,iatom,ind2)
207    end if
208    do ii=1,3
209      a11=a11+diff_t(ii)**2
210      a22=a22+diff_tmdt(ii)**2
211      a33=a33+diff_tpdt(ii)**2
212      a12=a12+diff_t(ii)   *diff_tmdt(ii)
213      b1 =b1 +diff_t(ii)   *diff_tpdt(ii)
214      b2 =b2 +diff_tmdt(ii)*diff_tpdt(ii)
215    end do
216 
217 !  Store reduced coordinates diffs in SCF history
218    scf_history%xreddiff(1:3,iatom,ind1new)=diff_tpdt(1:3)
219 
220  end do
221  ABI_FREE(xred_tpdt)
222  hasmoved=(a11>=tol10.or.a22>=tol10.or.a33>=tol10)
223 
224 !Compute alpha and beta
225  alpha=zero;beta=zero
226  if (hasmoved.and.ind1>0) then
227    ratio1=one;if (abs(a33)>=tol10) ratio1=(a11+a33-two*b1)/a33
228    ratio2=one;if (abs(a33)>=tol10) ratio2=(a11+a33-two*b2)/a33
229    detA=a11*a22-a12**2
230    if (abs(a11)>=tol10.and.(abs(a22)<tol10.or.abs(detA)<tol10)) then
231      alpha=b1/a11
232    else if (abs(a22)>=tol10.and.(abs(a11)<tol10.or.abs(detA)<tol10)) then
233      beta=b2/a22
234    else if (abs(ratio1)+abs(ratio2)<tol6) then
235      if (ind2>0) then
236        alpha=two;beta=-one
237      else
238        alpha=one
239      end if
240      write(message,'(3a,f4.1,a,f4.1)')&
241 &     'Ionic positions lead to a collinear system !',ch10,&
242 &     'Mixing coeffs have been set to: alpha=',alpha,' beta=',beta
243      ABI_WARNING(message)
244    else if (abs(a11)>=tol10.and.abs(a22)>=tol10) then
245      alpha=(b1*a22-b2*a12)/detA
246      beta =(b2*a11-b1*a12)/detA
247    end if
248  end if
249 
250 
251 !---------------------------------------------------------------
252 !----------- Contribution from delta_rho(t), delta_rho(t-dt)
253 !----------- and delta_rho(t-2dt) to predicted rho(t+dt)
254 !---------------------------------------------------------------
255 
256 !deltarho(t+dt) <- deltarho(t) + alpha.[deltarho(t)-deltarho(t-dt)]
257 !+ beta .[deltarho(t-dt)-deltarho(t-2dt)]
258 !Note: scf_history%deltarhor is updated at the same time
259 
260  ABI_MALLOC(deltarho,(nfft))
261  do ispden=1,dtset%nspden
262 
263    if (ispden==1) then
264      deltarho(:)=rhor(:,ispden)-scf_history%atmrho_last(:)
265    else if (ispden==2.and.dtset%nspden==2) then
266      deltarho(:)=rhor(:,ispden)-half*scf_history%atmrho_last(:)
267    end if
268 
269 
270 !  rho(t+dt) <- deltarho(t) + alpha.deltarho(t)
271    if (dtset%nspden/=4.or.ispden==1) then
272      rhor(:,ispden)=(one+alpha)*deltarho(:)
273    else
274      rhor(:,ispden)=(one+alpha)*rhor(:,ispden)
275    end if
276 
277    if (hasmoved) then
278 
279 !    rho(t+dt) <- -alpha.deltarho(t-dt) + beta.deltarho(t-dt)
280      if (abs(beta-alpha)>tol14.and.ind1>0) then
281        rhor(:,ispden)=rhor(:,ispden)+(beta-alpha)*scf_history%deltarhor(:,ispden,ind1)
282      end if
283 
284 !    rho(t+dt) <- -beta.deltarho(t-2dt)
285      if (abs(beta)>tol14.and.ind2>0) then
286        rhor(:,ispden)=rhor(:,ispden)-beta*scf_history%deltarhor(:,ispden,ind2)
287      end if
288 
289    end if
290 
291 !  Store deltarho(t) in history
292    if (dtset%nspden/=4.or.ispden==1) then
293      scf_history%deltarhor(:,ispden,ind1new)=deltarho(:)
294    else
295      scf_history%deltarhor(:,ispden,ind1new)=rhor(:,ispden)
296    end if
297 
298  end do
299 
300  ABI_FREE(deltarho)
301 
302 !---------------------------------------------------------------
303 !----------- Contribution from rho_at(t+dt) to predicted rho(t+dt)
304 !---------------------------------------------------------------
305 
306 !Determine whether a gaussian atomic density has to be used or not
307 !MG: Note that there's a small inconsistency between initro and extraprho because in initrho
308 ! we use `use_gaussian(ntypat)`.
309  usegauss=.true.
310  if (usepaw==0) usegauss = any(.not.psps%nctab(1:ntypat)%has_tvale)
311  if (usepaw==1) usegauss=(minval(pawtab(1:ntypat)%has_tvale)==0)
312  if (usegauss) then
313    optn2=3
314    ABI_MALLOC(gauss,(2,ntypat))
315    do itypat=1,ntypat
316      gauss(1,itypat)=zion(itypat)
317      gauss(2,itypat) = atom_length(dtset%densty(itypat,1),zion(itypat),znucl(itypat))
318    end do
319    call wrtout(std_out," Extrapolating rho(t+dt) using gaussian functions as atomic densities", "COLL")
320  else
321    optn2=2
322    ABI_MALLOC(gauss,(2,0))
323    call wrtout(std_out," Extrapolating rho(t+dt) using atomic densities taken from pseudos", "COLL")
324  end if
325 
326 !Compute rho_at(t+dt) as sum of atomic densities
327 !Note: scf_history%atmrho_last is updated at the same time
328  optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optv=0;optn=1
329  call atm2fft(atindx1,scf_history%atmrho_last,dummy_out1,dummy_out2,dummy_out3,&
330 & dummy_out4,gauss,gmet,gprimd,dummy_out5,dummy_out6,gsqcut,mgfft,mqgrid,dtset%natom,nattyp,&
331 & nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,qgrid,&
332 & dummy3,dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,usepaw,dummy_in,dummy_in,dummy_in,dummy2,dummy_in,&
333 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
334 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
335  ABI_FREE(gauss)
336 
337 !Take eventually into account jellium slab
338  if (dtset%jellslab/=0) then
339    option=2
340    ABI_MALLOC(work1,(nfft))
341    ABI_MALLOC(work2,(nfft,1))
342    ABI_MALLOC(work3,(2,nfft))
343    work2(:,1)=scf_history%atmrho_last(:)
344    call jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,1,option,&
345 &   dtset%slabwsrad,work3,work2,rprimd,work1,dtset%slabzbeg,dtset%slabzend)
346    scf_history%atmrho_last(:)=work2(:,1)
347    ABI_FREE(work1)
348    ABI_FREE(work2)
349    ABI_FREE(work3)
350  end if
351 
352 !Add rho_at(t+dt) to rho(t+dt)
353  rhor(:,1)=rhor(:,1)+scf_history%atmrho_last(:)
354  if (dtset%nspden==2) rhor(:,2)=rhor(:,2)+half*scf_history%atmrho_last(:)
355 
356 !---------------------------------------------------------------
357 !----------- Extrapolation of PAW rhoij occupancy matrixes
358 !---------------------------------------------------------------
359 
360  if (usepaw==1) then
361 
362    if (ind2==0) then
363      call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
364 &                nspden=dtset%nspden,spnorb=dtset%pawspnorb,cpxocc=dtset%pawcpxocc)
365      call pawrhoij_alloc(scf_history%pawrhoij(:,ind1new),cplex_rhoij,nspden_rhoij,&
366 &     dtset%nspinor,dtset%nsppol,dtset%typat,pawtab=pawtab,&
367 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
368    end if
369 
370    do iatom=1,my_natom
371 
372      nspden_rhoij=pawrhoij(iatom)%nspden
373      lmn2_size=pawrhoij(iatom)%lmn2_size
374      cplex_rhoij=pawrhoij(iatom)%cplex_rhoij;dplex=cplex_rhoij-1
375      qphase_rhoij=pawrhoij(iatom)%qphase
376 
377      if (hasmoved) then
378        ABI_MALLOC(rhoijtmp,(cplex_rhoij*qphase_rhoij*lmn2_size,nspden_rhoij))
379        rhoijtmp=zero
380 
381        do ispden=1,nspden_rhoij
382          do iq=1,qphase_rhoij
383            iq0=merge(0,cplex_rhoij*lmn2_size,iq==1)
384 
385 !          rhoij(t+dt) <- rhoij(t) + alpha.rhoij(t)
386            fact=one+alpha
387            jrhoij=1+iq0
388            do irhoij=1,pawrhoij(iatom)%nrhoijsel
389              klmn=cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex+iq0
390              rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) &
391 &             +fact*pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
392              jrhoij=jrhoij+cplex_rhoij
393            end do
394 
395 !          rhoij(t+dt) <- -alpha.rhoij(t-dt) + beta.rhoij(t-dt)
396            if (abs(beta-alpha)>tol14.and.ind1>0) then
397              fact=beta-alpha
398              jrhoij=1+iq0
399              do irhoij=1,scf_history%pawrhoij(iatom,ind1)%nrhoijsel
400                klmn=cplex_rhoij*scf_history%pawrhoij(iatom,ind1)%rhoijselect(irhoij)-dplex+iq0
401                rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) &
402 &               +fact*scf_history%pawrhoij(iatom,ind1)%rhoijp(jrhoij:jrhoij+dplex,ispden)
403                jrhoij=jrhoij+cplex_rhoij
404              end do
405            end if
406 
407 !          rho(t+dt) <- -beta.rhoij(t-2dt)
408            if (abs(beta)>tol14.and.ind2>0) then
409              fact=-beta
410              jrhoij=1+iq0
411              do irhoij=1,scf_history%pawrhoij(iatom,ind2)%nrhoijsel
412                klmn=cplex_rhoij*scf_history%pawrhoij(iatom,ind2)%rhoijselect(irhoij)-dplex+iq0
413                rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) &
414 &               +fact*scf_history%pawrhoij(iatom,ind2)%rhoijp(jrhoij:jrhoij+dplex,ispden)
415                jrhoij=jrhoij+cplex_rhoij
416              end do
417            end if
418 
419          end do ! iq
420        end do !ispden
421      end if !hasmoved
422 
423 !    Store rhoij(t) in history
424 !    (cannot use pawrhoij_copy here because update for single atom)
425      nselect=pawrhoij(iatom)%nrhoijsel
426      scf_history%pawrhoij(iatom,ind1new)%nrhoijsel=nselect
427      scf_history%pawrhoij(iatom,ind1new)%rhoijselect(:)=0
428      scf_history%pawrhoij(iatom,ind1new)%rhoijselect(1:nselect)=pawrhoij(iatom)%rhoijselect(1:nselect)
429      scf_history%pawrhoij(iatom,ind1new)%rhoijp(1:cplex_rhoij*nselect,1:nspden_rhoij)= &
430 &     pawrhoij(iatom)%rhoijp(1:cplex_rhoij*nselect,1:nspden_rhoij)
431 
432 !    Select non-zero values of rhoij(t+dt)
433      if (hasmoved) then
434        call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,pawrhoij(iatom)%nrhoijsel,&
435 &                           cplex_rhoij,qphase_rhoij,lmn2_size,nspden_rhoij,&
436 &                           rhoij_input=rhoijtmp)
437        ABI_FREE(rhoijtmp)
438      end if
439 
440    end do !iatom
441  end if !usepaw
442 
443 
444  scf_history%alpha=alpha
445  scf_history%beta=beta
446 
447 
448 
449 !---------------------------------------------------------------
450 !----------- End
451 !---------------------------------------------------------------
452 
453  if(scf_history%usecg==1) then
454    if (hasmoved) then
455      if (dtset%extrapwf==1) then
456        call extrapwf(atindx,atindx1,cg,dtset,istep,kg,mcg,mgfft,mpi_enreg,nattyp,&
457 &       ngfft,npwarr,ntypat,pawtab,psps,rprimd,scf_history,usepaw,xred_old,ylm)
458      elseif(dtset%extrapwf==2) then
459        scf_history%hindex(1)=ind1
460        scf_history%hindex(2)=ind2
461        scf_history%hindex(3)=ind1new
462        call extrapwf_biortho(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,&
463 &         nattyp,npwarr,pawtab,scf_history)
464      end if
465    else
466      scf_history%cg(:,:,2)=zero
467    end if
468 
469  end if
470 !Rotate history indexes
471  scf_history%hindex(1)=ind1new
472  scf_history%hindex(2)=ind2new
473 
474 
475 end subroutine extraprho

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) [[cite:Arias1992]]

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).

SOURCE

 520 subroutine extrapwf(atindx,atindx1,cg,dtset,istep,kg,mcg,mgfft,mpi_enreg,&
 521 & nattyp,ngfft,npwarr,ntypat,pawtab,psps,rprimd,scf_history,usepaw,xred_old,ylm)
 522 
 523 !Arguments ------------------------------------
 524 !scalars
 525  integer,intent(in) :: istep,mcg,mgfft,ntypat,usepaw
 526  type(MPI_type),intent(in) :: mpi_enreg
 527  type(dataset_type),intent(in) :: dtset
 528  type(scf_history_type),intent(inout) :: scf_history
 529  type(pseudopotential_type),intent(in) :: psps
 530 !arrays
 531  integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem),nattyp(ntypat),ngfft(18)
 532  integer,intent(in) :: npwarr(dtset%nkpt)
 533  real(dp),intent(in) :: rprimd(3,3)
 534  real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
 535  real(dp), intent(inout) :: cg(2,mcg)
 536  real(dp),intent(in) :: xred_old(3,dtset%natom)
 537  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
 538 
 539 
 540 !Local variables-------------------------------
 541 !scalars
 542  integer :: ia,iat,iatom,iband_max,iband_max1,iband_min,iband_min1,ibd,ibg,iblockbd,iblockbd1,icg,icgb,icgb1,icgb2
 543  integer :: ierr,ig,ii,ikpt,ilmn1,ilmn2,inc,ind1,ind2,iorder_cprj
 544  integer :: isize,isppol,istep1,istwf_k,itypat,klmn,me_distrb,my_nspinor
 545  integer :: nband_k,nblockbd,nprocband,npw_k,npw_nk,spaceComm_band
 546  real(dp) :: dotr,dotr1,doti,doti1,eigval
 547  !character(len=500) :: message
 548 !arrays
 549  real(dp) :: alpha(2),beta(2),gmet(3,3),gprimd(3,3),rmet(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom),ucvol
 550  integer,allocatable :: bufsize(:),bufsize_wf(:),bufdisp(:),bufdisp_wf(:),dimcprj(:),npw_block(:),npw_disp(:)
 551  real(dp),allocatable :: al(:,:),anm(:),cwavef(:,:),cwavef1(:,:),cwavef_tmp(:,:),deltawf1(:,:),deltawf2(:,:)
 552  real(dp),allocatable :: eig(:),evec(:,:)
 553  real(dp),allocatable :: unm(:,:,:)
 554  real(dp),allocatable :: work(:,:),work1(:,:),wf1(:,:),ylmgr_k(:,:,:),zhpev1(:,:),zhpev2(:)
 555  complex(dpc),allocatable :: unm_tmp(:,:),anm_tmp(:,:)
 556  type(pawcprj_type),allocatable :: cprj(:,:),cprj_k(:,:),cprj_k1(:,:),cprj_k2(:,:),cprj_k3(:,:),cprj_k4(:,:)
 557 !complex(dpc) :: aa
 558 
 559 ! *************************************************************************
 560 
 561  if (istep==0) return
 562 
 563 !Useful array
 564  if (usepaw==1) then
 565    ABI_MALLOC(dimcprj,(dtset%natom))
 566    call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O')
 567  end if
 568 
 569 !Metric
 570  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 571 
 572 !History indexes
 573  ind1=1;ind2=2
 574 
 575 !First step
 576  if (istep==1) then
 577    scf_history%cg(:,:,ind1)=cg(:,:)
 578 !  scf_history%cg(:,:,ind2)=zero
 579    scf_history%cg(:,:,ind2)= cg(:,:)
 580    if(usepaw==1) then
 581 !    WARNING: THIS SECTION IS USELESS; NOW crpj CAN BE READ FROM SCFCV
 582      call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred_old)
 583      iatom=0 ; iorder_cprj=0
 584      call pawcprj_alloc(scf_history%cprj(:,:,ind1),0,dimcprj)
 585      call pawcprj_alloc(scf_history%cprj(:,:,ind2),0,dimcprj)
 586      ABI_MALLOC(ylmgr_k,(dtset%mpw,3,0))
 587      call ctocprj(atindx,cg,1,scf_history%cprj(:,:,ind1),gmet,gprimd,&
 588 &     iatom,0,iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,scf_history%mcprj,&
 589 &     dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,&
 590 &     dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,&
 591 &     dtset%nloalg,npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,dtset%ntypat,&
 592 &     dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,0,&
 593 &     xred_old,ylm,ylmgr_k)
 594      ABI_FREE(ylmgr_k)
 595 !    call pawcprj_set_zero(scf_history%cprj(:,:,ind2))
 596      call pawcprj_copy(scf_history%cprj(:,:,ind1),scf_history%cprj(:,:,ind2))
 597    end if
 598  else
 599 
 600 !From 2nd step
 601 
 602 !  Init parallelism
 603    me_distrb=mpi_enreg%me_kpt
 604    if (mpi_enreg%paral_kgb==1.or.mpi_enreg%paralbd==1) then
 605      spaceComm_band=mpi_enreg%comm_band
 606      nprocband=mpi_enreg%nproc_band
 607    else
 608      spaceComm_band=xmpi_comm_self
 609      nprocband=1
 610    end if
 611 
 612 !  For the moment sequential part only
 613    nprocband=1
 614 
 615 !  Additional statements if band-fft parallelism
 616    if (nprocband>1) then
 617      ABI_MALLOC(npw_block,(nprocband))
 618      ABI_MALLOC(npw_disp,(nprocband))
 619      ABI_MALLOC(bufsize,(nprocband))
 620      ABI_MALLOC(bufdisp,(nprocband))
 621      ABI_MALLOC(bufsize_wf,(nprocband))
 622      ABI_MALLOC(bufdisp_wf,(nprocband))
 623    end if
 624 
 625    icg=0
 626    ibg=0
 627 
 628    if(usepaw==1) then
 629 !    WARNING: THIS SECTION IS USELESS; NOW cprj CAN BE READ FROM SCFCV
 630      call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred_old)
 631      ABI_MALLOC(cprj,(dtset%natom,scf_history%mcprj))
 632      call pawcprj_alloc(cprj,0,dimcprj)
 633      iatom=0 ; iorder_cprj=0
 634      ABI_MALLOC(ylmgr_k,(dtset%mpw,3,0))
 635      call ctocprj(atindx,cg,1,cprj,gmet,gprimd,iatom,0,iorder_cprj,&
 636 &     dtset%istwfk,kg,dtset%kptns,mcg,scf_history%mcprj,dtset%mgfft,&
 637 &     dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,&
 638 &     nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,dtset%nloalg,&
 639 &     npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,&
 640 &     ph1d,psps,rmet,dtset%typat,ucvol,0,xred_old,&
 641 &     ylm,ylmgr_k)
 642      ABI_FREE(ylmgr_k)
 643    end if  ! end usepaw=1
 644 
 645 !  LOOP OVER SPINS
 646    my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
 647    do isppol=1,dtset%nsppol
 648 
 649 !    BIG FAT k POINT LOOP
 650      do ikpt=1,dtset%nkpt
 651 
 652 !      Select k point to be treated by this proc
 653        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
 654        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
 655 
 656        istwf_k=dtset%istwfk(ikpt)
 657 
 658 !      Retrieve number of plane waves
 659        npw_k=npwarr(ikpt)
 660        if (nprocband>1) then
 661 !        Special treatment for band-fft //
 662          call xmpi_allgather(npw_k,npw_block,spaceComm_band,ierr)
 663          npw_nk=sum(npw_block);npw_disp(1)=0
 664          do ii=2,nprocband
 665            npw_disp(ii)=npw_disp(ii-1)+npw_block(ii-1)
 666          end do
 667        else
 668          npw_nk=npw_k
 669        end if
 670 
 671 !      Allocate arrays for a wave-function (or a block of WFs)
 672        ABI_MALLOC(cwavef,(2,npw_nk*my_nspinor))
 673        ABI_MALLOC(cwavef1,(2,npw_nk*my_nspinor))
 674        if (nprocband>1) then
 675          isize=2*my_nspinor;bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:)
 676          isize=2*my_nspinor*npw_k;bufsize_wf(:)=isize
 677          do ii=1,nprocband
 678            bufdisp_wf(ii)=(ii-1)*isize
 679          end do
 680        end if
 681 
 682 !      Subspace alignment
 683 
 684 !      Loop over bands or blocks of bands
 685        nblockbd=nband_k/nprocband
 686        icgb=icg
 687 
 688        if(usepaw==1) then
 689          ABI_MALLOC( cprj_k,(dtset%natom,my_nspinor*nblockbd))
 690          call pawcprj_alloc(cprj_k,cprj(1,1)%ncpgr,dimcprj)
 691          call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,dtset%mband,&
 692 &         dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,&
 693 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
 694          ABI_MALLOC( cprj_k1,(dtset%natom,my_nspinor*nblockbd))
 695          call pawcprj_alloc(cprj_k1,scf_history%cprj(1,1,ind1)%ncpgr,dimcprj)
 696          call pawcprj_get(atindx1,cprj_k1,scf_history%cprj(:,:,ind1),dtset%natom,1,ibg,ikpt,1,isppol,&
 697 &         dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,&
 698 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
 699          ABI_MALLOC( cprj_k2,(dtset%natom,my_nspinor*nblockbd))
 700          call pawcprj_alloc(cprj_k2,scf_history%cprj(1,1,ind2)%ncpgr,dimcprj)
 701          call pawcprj_get(atindx1,cprj_k2,scf_history%cprj(:,:,ind2),dtset%natom,1,ibg,ikpt,1,isppol,&
 702 &         dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,&
 703 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
 704        end if  !end usepaw=1
 705 
 706        ABI_MALLOC(unm,(2,nblockbd,nblockbd))
 707        unm=zero
 708        icgb2=0
 709 
 710        do iblockbd=1,nblockbd
 711          iband_min=1+(iblockbd-1)*nprocband
 712          iband_max=iblockbd*nprocband
 713 
 714          if(xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then
 715            if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min,iband_max,isppol,me_distrb)) cycle
 716          end if
 717 
 718 !        Extract wavefunction information
 719          if (nprocband>1) then
 720 !          Special treatment for band-fft //
 721            ABI_MALLOC(cwavef_tmp,(2,npw_k*my_nspinor*nprocband))
 722            do ig=1,npw_k*my_nspinor*nprocband
 723              cwavef_tmp(1,ig)=cg(1,ig+icgb)
 724              cwavef_tmp(2,ig)=cg(2,ig+icgb)
 725            end do
 726            call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef,bufsize,bufdisp,spaceComm_band,ierr)
 727            ABI_FREE(cwavef_tmp)
 728          else
 729            do ig=1,npw_k*my_nspinor
 730              cwavef(1,ig)=cg(1,ig+icgb)
 731              cwavef(2,ig)=cg(2,ig+icgb)
 732            end do
 733          end if
 734 
 735          icgb1=icg
 736 
 737          do iblockbd1=1,nblockbd
 738            iband_min1=1+(iblockbd1-1)*nprocband
 739            iband_max1=iblockbd1*nprocband
 740 
 741            if(xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then
 742              if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min1,iband_max1,isppol,me_distrb)) cycle
 743            end if
 744 
 745 !          Extract wavefunction information
 746 
 747            if (nprocband>1) then
 748 !            Special treatment for band-fft //
 749              ABI_MALLOC(cwavef_tmp,(2,npw_k*my_nspinor*nprocband))
 750              do ig=1,npw_k*my_nspinor*nprocband
 751                cwavef_tmp(1,ig)=scf_history%cg(1,ig+icgb1,ind1)
 752                cwavef_tmp(2,ig)=scf_history%cg(2,ig+icgb1,ind1)
 753              end do
 754              call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef1,bufsize,bufdisp,spaceComm_band,ierr)
 755              ABI_FREE(cwavef_tmp)
 756            else
 757              do ig=1,npw_k*my_nspinor
 758                cwavef1(1,ig)=scf_history%cg(1,ig+icgb1,ind1)
 759                cwavef1(2,ig)=scf_history%cg(2,ig+icgb1,ind1)
 760              end do
 761            end if
 762 
 763 !          Calculate Unm=<psi_nk(t)|S|psi_mk(t-dt)>
 764            call dotprod_g(dotr,doti,istwf_k,npw_k*my_nspinor,2,cwavef,cwavef1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
 765            if(usepaw==1) then
 766              ia =0
 767              do itypat=1,ntypat
 768                do iat=1+ia,nattyp(itypat)+ia
 769                  do ilmn1=1,pawtab(itypat)%lmn_size
 770                    do ilmn2=1,ilmn1
 771                      klmn=((ilmn1-1)*ilmn1)/2+ilmn2
 772                      dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)+&
 773 &                     cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2))
 774                      doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)-&
 775 &                     cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2))
 776                    end do
 777                    do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
 778                      klmn=((ilmn2-1)*ilmn2)/2+ilmn1
 779                      dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)+&
 780 &                     cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2))
 781                      doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)-&
 782 &                     cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2))
 783                    end do
 784                  end do
 785                end do
 786                ia=ia+nattyp(itypat)
 787              end do
 788            end if
 789 !          unm(1,iblockbd,iblockbd1)=dotr
 790 !          unm(2,iblockbd,iblockbd1)=doti
 791            unm(1,iblockbd1,iblockbd)=dotr
 792            unm(2,iblockbd1,iblockbd)=doti
 793 !          End loop over bands iblockbd1
 794            icgb1=icgb1+npw_k*my_nspinor*nprocband
 795 
 796          end do
 797 
 798 !        End loop over bands iblockbd
 799          icgb2=icgb2+npw_k*my_nspinor*nprocband
 800          icgb=icgb+npw_k*my_nspinor*nprocband
 801        end do
 802 
 803 !      write(std_out,*) 'UNM'
 804 !      do iblockbd=1,nblockbd
 805 !      write(std_out,11) (unm(1,iblockbd,iblockbd1),unm(2,iblockbd,iblockbd1),iblockbd1=1,nblockbd)
 806 !      end do
 807 !      11 format(12(1x,f9.5),a)
 808 !      Compute A=tU^*U
 809        ABI_MALLOC(unm_tmp,(nblockbd,nblockbd))
 810        ABI_MALLOC(anm_tmp,(nblockbd,nblockbd))
 811        ABI_MALLOC(anm,(nblockbd*(nblockbd+1)))
 812        unm_tmp(:,:)=cmplx(unm(1,:,:),unm(2,:,:),kind=dp)
 813        call zgemm('C','N',nblockbd,nblockbd,nblockbd,dcmplx(1._dp), unm_tmp,nblockbd, &
 814 &       unm_tmp,nblockbd,dcmplx(0._dp),anm_tmp,nblockbd)
 815        do iblockbd=1,nblockbd
 816          do iblockbd1=iblockbd,nblockbd
 817            ii=iblockbd1*(iblockbd1-1)+2*(iblockbd-1)+1
 818            anm(ii)=real(anm_tmp(iblockbd,iblockbd1))
 819            anm(ii+1)=aimag(anm_tmp(iblockbd,iblockbd1))
 820          end do
 821        end do
 822        call hermit(anm,anm,ierr,nblockbd)
 823 !      aa=dcmplx(0._dp)
 824 !      do iblockbd=1,nblockbd
 825 !      aa=aa+conjg(unm_tmp(iblockbd,1))*unm_tmp(iblockbd,1)
 826 !      end do
 827 !      write(std_out,*) 'tU*U', aa
 828 !      write(std_out,*) 'ANM_tmp'
 829 !      do iblockbd=1,nblockbd
 830 !      write(std_out,11) (anm_tmp(iblockbd,iblockbd1),iblockbd1=1,nblockbd)
 831 !      end do
 832 !      write(std_out,*) 'ANM'
 833 !      do iblockbd=1,nblockbd*(nblockbd+1)
 834 !      write(std_out,11) anm(iblockbd)
 835 !      end do
 836 
 837 !      Diagonalize A
 838        ABI_MALLOC(eig,(nblockbd))
 839        ABI_MALLOC(evec,(2*nblockbd,nblockbd))
 840        ABI_MALLOC(zhpev1,(2,2*nblockbd-1))
 841        ABI_MALLOC(zhpev2,(3*nblockbd-2))
 842        call zhpev('V','U',nblockbd,anm,eig,evec,nblockbd,zhpev1,&
 843 &       zhpev2,ierr)
 844        ABI_FREE(anm)
 845        ABI_FREE(zhpev1)
 846        ABI_FREE(zhpev2)
 847 !      aa=dcmplx(0._dp)
 848 !      do iblockbd=1,nblockbd
 849 !      aa=aa+anm_tmp(1,iblockbd)*cmplx(evec((2*iblockbd-1),1),evec(2*iblockbd,1),kind=dp)
 850 !      end do
 851 !      write(std_out,*) 'EIG', aa, eig(1)*evec(1,1),eig(1)*evec(2,1)
 852 
 853 !      Compute A'=evec*tU^/sqrt(eig)
 854        call zgemm('C','C',nblockbd,nblockbd,nblockbd,dcmplx(1._dp),evec,nblockbd, &
 855 &       unm_tmp,nblockbd,dcmplx(0._dp),anm_tmp,nblockbd)
 856        do iblockbd=1,nblockbd
 857          eigval=dsqrt(eig(iblockbd))
 858          do iblockbd1=1,nblockbd
 859            anm_tmp(iblockbd,iblockbd1)=anm_tmp(iblockbd,iblockbd1)/eigval
 860          end do
 861        end do
 862 
 863 !      Compute tA^A'to come back to the initial subspace for the cg's
 864 
 865        call zgemm('N','N',nblockbd,nblockbd,nblockbd,dcmplx(1._dp),evec,nblockbd, &
 866 &       anm_tmp,nblockbd,dcmplx(0._dp),unm_tmp,nblockbd)
 867        anm_tmp=unm_tmp
 868 !      write(std_out,*) 'ANM_tmp'
 869 !      do iblockbd=1,nblockbd
 870 !      write(std_out,11) (anm_tmp(iblockbd,iblockbd1),iblockbd1=1,nblockbd)
 871 !      end do
 872 
 873 !      Wavefunction alignment (istwfk=1 ?)
 874        ABI_MALLOC(work,(2,npw_nk*my_nspinor*nblockbd))
 875        ABI_MALLOC(work1,(2,my_nspinor*nblockbd*npw_nk))
 876        work1(:,:)=scf_history%cg(:,icg+1:icg+my_nspinor*nblockbd*npw_nk,ind1)
 877        call zgemm('N','N',npw_nk*my_nspinor,nblockbd,nblockbd,dcmplx(1._dp), &
 878 &       work1,npw_nk*my_nspinor, &
 879 &       anm_tmp,nblockbd,dcmplx(0._dp),work,npw_nk*my_nspinor)
 880        scf_history%cg(:,1+icg:npw_nk*my_nspinor*nblockbd+icg,ind1)=work(:,:)
 881 
 882        work1(:,:)=scf_history%cg(:,icg+1:icg+my_nspinor*nblockbd*npw_nk,ind2)
 883        call zgemm('N','N',npw_nk*my_nspinor,nblockbd,nblockbd,dcmplx(1._dp), &
 884 &       work1,npw_nk*my_nspinor, &
 885 &       anm_tmp,nblockbd,dcmplx(0._dp),work,npw_nk*my_nspinor)
 886        scf_history%cg(:,1+icg:npw_nk*my_nspinor*nblockbd+icg,ind2)=work(:,:)
 887        ABI_FREE(work1)
 888 !      If paw, must also align cprj:
 889        if (usepaw==1) then
 890 !        New version (MT):
 891          ABI_MALLOC(cprj_k3,(dtset%natom,my_nspinor))
 892          ABI_MALLOC(cprj_k4,(dtset%natom,my_nspinor))
 893          call pawcprj_alloc(cprj_k3,cprj_k1(1,1)%ncpgr,dimcprj)
 894          call pawcprj_alloc(cprj_k4,cprj_k2(1,1)%ncpgr,dimcprj)
 895          ABI_MALLOC(al,(2,nblockbd))
 896          do iblockbd=1,nblockbd
 897            ii=(iblockbd-1)*my_nspinor
 898            do iblockbd1=1,nblockbd
 899              al(1,iblockbd1)=real (anm_tmp(iblockbd,iblockbd1))
 900              al(2,iblockbd1)=aimag(anm_tmp(iblockbd,iblockbd1))
 901            end do
 902            call pawcprj_lincom(al,cprj_k1,cprj_k3,nblockbd)
 903            call pawcprj_lincom(al,cprj_k2,cprj_k4,nblockbd)
 904            call pawcprj_copy(cprj_k3,cprj_k1(:,ii+1:ii+my_nspinor))
 905            call pawcprj_copy(cprj_k4,cprj_k2(:,ii+1:ii+my_nspinor))
 906          end do
 907          ABI_FREE(al)
 908 !        Old version (FJ):
 909 !        allocate( cprj_k3(dtset%natom,my_nspinor*nblockbd))
 910 !        call pawcprj_alloc(cprj_k3,cprj_k1(1,1)%ncpgr,dimcprj)
 911 !        allocate( cprj_k4(dtset%natom,my_nspinor*nblockbd))
 912 !        call pawcprj_alloc(cprj_k4,cprj_k2(1,1)%ncpgr,dimcprj)
 913 !        beta(1)=one;beta(2)=zero
 914 !        do iblockbd=1,nblockbd*my_nspinor
 915 !        do iblockbd1=1,nblockbd*my_nspinor
 916 !        alpha(1)=real(anm_tmp(iblockbd,iblockbd1));alpha(2)=aimag(anm_tmp(iblockbd,iblockbd1))
 917 !        call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd1:iblockbd1),cprj_k3(:,iblockbd:iblockbd))
 918 !        call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd1:iblockbd1),cprj_k4(:,iblockbd:iblockbd))
 919 !        end do
 920 !        end do
 921 !        call pawcprj_copy(cprj_k3,cprj_k1)
 922 !        call pawcprj_copy(cprj_k4,cprj_k2)
 923 
 924          call pawcprj_free(cprj_k3)
 925          call pawcprj_free(cprj_k4)
 926          ABI_FREE(cprj_k3)
 927          ABI_FREE(cprj_k4)
 928        end if
 929        ABI_FREE(anm_tmp)
 930        ABI_FREE(unm_tmp)
 931        ABI_FREE(work)
 932 
 933 !      Wavefunction extrapolation
 934        ibd=0
 935        inc=npw_nk*my_nspinor
 936        ABI_MALLOC(deltawf2,(2,npw_nk*my_nspinor))
 937        ABI_MALLOC(wf1,(2,npw_nk*my_nspinor))
 938        ABI_MALLOC(deltawf1,(2,npw_nk*my_nspinor))
 939        do iblockbd=1,nblockbd
 940          deltawf2(:,:)=scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind2)
 941          wf1(:,:)=scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind1)
 942 !        wf1(2,1)=zero;deltawf2(2,1)=zero
 943 
 944          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),&
 945 &         mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
 946          call dotprod_g(dotr1,doti1,istwf_k,npw_nk*my_nspinor,2,cg(:,icg+1+ibd:ibd+icg+inc),wf1,&
 947 &         mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
 948          if(usepaw==1) then
 949            ia =0
 950            do itypat=1,ntypat
 951              do iat=1+ia,nattyp(itypat)+ia
 952                do ilmn1=1,pawtab(itypat)%lmn_size
 953                  do ilmn2=1,ilmn1
 954                    klmn=((ilmn1-1)*ilmn1)/2+ilmn2
 955                    dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)+&
 956 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2))
 957                    doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)-&
 958 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2))
 959                    dotr1=dotr1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)+&
 960 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2))
 961                    doti1=doti1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)-&
 962 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2))
 963                  end do
 964                  do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
 965                    klmn=((ilmn2-1)*ilmn2)/2+ilmn1
 966                    dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)+&
 967 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2))
 968                    doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)-&
 969 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2))
 970                    dotr1=dotr1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)+&
 971 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2))
 972                    doti1=doti1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)-&
 973 &                   cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2))
 974                  end do
 975                end do
 976              end do
 977              ia=ia+nattyp(itypat)
 978            end do
 979          end if
 980          dotr=sqrt(dotr**2+doti**2)
 981          dotr1=sqrt(dotr1**2+doti1**2)
 982          write(std_out,*)'DOTR, DOTR1',dotr,dotr1
 983          dotr=dotr1/dotr
 984          write(std_out,*)'DOTR',dotr
 985          deltawf1=zero
 986          if(dotr>=0.9d0) then
 987            deltawf1(:,:)=cg(:,icg+1+ibd:ibd+icg+inc)-wf1(:,:)
 988            if(usepaw==1) then
 989              alpha(1)=one;alpha(2)=zero
 990              beta(1)=-one;beta(2)=zero
 991              ia =0
 992              call pawcprj_zaxpby(alpha,beta,cprj_k(:,iblockbd:iblockbd),cprj_k1(:,iblockbd:iblockbd))
 993            end if
 994            istep1=istep
 995          else
 996            istep1=1
 997          end if
 998          scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind1)=cg(:,icg+1+ibd:ibd+icg+inc)
 999          scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind2)=deltawf1(:,:)
1000          if(usepaw==1) then
1001            call pawcprj_put(atindx1,cprj_k,scf_history%cprj(:,:,ind1),dtset%natom,1,ibg,ikpt,1,isppol,&
1002 &           dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
1003 &           mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1004            call pawcprj_put(atindx1,cprj_k1,scf_history%cprj(:,:,ind2),dtset%natom,1,ibg,ikpt,1,isppol,&
1005 &           dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
1006 &           mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1007          end if
1008 
1009 !        if(istep1>=3) then
1010          cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+scf_history%alpha*deltawf1(:,:) &
1011 &         +scf_history%beta *deltawf2(:,:)
1012 
1013 !        to be used later
1014 !        if(usepaw==1) then
1015 !        alpha(2)=zero
1016 !        beta(1)=one;beta(2)=zero
1017 !        alpha(1)=scf_history%alpha
1018 !        call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
1019 !        alpha(1)=scf_history%beta
1020 !        call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
1021 !        call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,&
1022 !        &    dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
1023 !        &    mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1024 !        end if
1025 !        else if (istep1==2) then
1026 !          cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+scf_history%alpha*deltawf1(:,:)+scf_history%beta*wf1(:,:)
1027 !       !     cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+deltawf1(:,:)
1028 !        if(usepaw==1) then
1029 !        alpha(2)=zero
1030 !        beta(1)=one;beta(2)=zero
1031 !        alpha(1)=scf_history%alpha
1032 !        call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
1033 !        alpha(1)=scf_history%beta
1034 !        call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
1035 !        call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,&
1036 !        &    dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
1037 !        &    mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1038 !        end if
1039 !        end if
1040          ibd=ibd+inc
1041        end do ! end loop on iblockbd
1042 
1043        ABI_FREE(deltawf1)
1044        ABI_FREE(deltawf2)
1045        ABI_FREE(wf1)
1046        ABI_FREE(cwavef)
1047        ABI_FREE(cwavef1)
1048        ABI_FREE(eig)
1049        ABI_FREE(evec)
1050        ABI_FREE(unm)
1051        if(usepaw==1) then
1052          call pawcprj_free(cprj_k)
1053          ABI_FREE(cprj_k)
1054          call pawcprj_free(cprj_k1)
1055          ABI_FREE(cprj_k1)
1056          call pawcprj_free(cprj_k2)
1057          ABI_FREE(cprj_k2)
1058        end if
1059 
1060        ibg=ibg+my_nspinor*nband_k
1061        icg=icg+my_nspinor*nband_k*npw_k
1062 
1063 !      End big k point loop
1064      end do
1065 !    End loop over spins
1066    end do
1067 
1068    if(usepaw==1) then
1069      call pawcprj_free(cprj)
1070      ABI_FREE(cprj)
1071    end if
1072    if (nprocband>1) then
1073      ABI_FREE(npw_block)
1074      ABI_FREE(npw_disp)
1075      ABI_FREE(bufsize)
1076      ABI_FREE(bufdisp)
1077      ABI_FREE(bufsize_wf)
1078      ABI_FREE(bufdisp_wf)
1079    end if
1080 
1081  end if ! istep>=2
1082 
1083  if (usepaw==1) then
1084    ABI_FREE(dimcprj)
1085  end if
1086 
1087 end subroutine extrapwf

ABINIT/extrapwf_biortho [ Functions ]

[ Top ] [ Functions ]

NAME

 extrapwf_biortho

FUNCTION

 Extrapolate wavefunctions for new ionic positions
 from values of wavefunctions of previous SCF cycle.
 Use biorthogonal algorithm proposed XG

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
  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

SOURCE

1123  subroutine extrapwf_biortho(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,&
1124 & nattyp,npwarr,pawtab,scf_history_wf)
1125 
1126  !use m_scf_history
1127  use m_cgcprj,  only : dotprod_set_cgcprj,cgcprj_cholesky,lincom_cgcprj
1128 
1129 !Arguments ------------------------------------
1130 !scalars
1131  integer,intent(in) :: istep,mcg,mcprj
1132  type(MPI_type),intent(in) :: mpi_enreg
1133  type(dataset_type),intent(in) :: dtset
1134  type(scf_history_type),intent(inout) :: scf_history_wf
1135 !arrays
1136  integer,intent(in) :: atindx1(dtset%natom),nattyp(dtset%ntypat)
1137  integer,intent(in) :: npwarr(dtset%nkpt)
1138  real(dp), intent(inout) :: cg(2,mcg)
1139  type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj)
1140  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
1141 
1142 !Local variables-------------------------------
1143 !scalars
1144  integer :: hermitian
1145  integer :: ibdmix,ibg,ibg_hist,icg,icg_hist !,iband
1146  integer :: ierr,ikpt,indh,ind1,ind2,ind1new,inplace
1147  integer :: isppol,istwf_k,kk,me_distrb,mband,my_nspinor,mcprj_k
1148  integer :: nband_k,nbdmix,nbdmax,npw_k,ntypat
1149  integer :: spaceComm_band,usepaw
1150  real(dp) :: alpha,beta !,dotr,doti
1151 
1152 !arrays
1153  integer,allocatable :: ipiv(:),dimcprj(:)
1154  real(dp),allocatable ::psi_ortho(:,:),mmn(:,:,:)
1155  real(dp),allocatable :: smn(:,:,:)
1156  type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kh(:,:)
1157 
1158 ! *************************************************************************
1159 
1160  if (istep==0) return
1161 
1162  ntypat=dtset%ntypat
1163  usepaw=dtset%usepaw
1164  mband=dtset%mband
1165  nbdmax=dtset%mband
1166  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
1167  me_distrb=mpi_enreg%me_kpt
1168  spaceComm_band=xmpi_comm_self
1169 
1170 !scf_history_wf%alpha contains dtset%wfmix
1171  alpha=scf_history_wf%alpha
1172  beta=scf_history_wf%beta
1173  ind1=scf_history_wf%hindex(1)
1174  ind2=scf_history_wf%hindex(2)
1175  ind1new=scf_history_wf%hindex(3)
1176  icg=0
1177  icg_hist=0
1178  ibg=0
1179  ibg_hist=0
1180 
1181 !Useful array
1182  ABI_MALLOC(dimcprj,(dtset%natom))
1183  if (usepaw==1) then
1184    call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O')
1185  end if
1186 
1187  if(istep==1)then
1188    do indh=1,scf_history_wf%history_size
1189      call pawcprj_alloc(scf_history_wf%cprj(:,:,indh),0,dimcprj)
1190    end do
1191  end if
1192 
1193  mcprj_k=my_nspinor*nbdmax
1194  ABI_MALLOC(cprj_k,(dtset%natom,mcprj_k))
1195  ABI_MALLOC(cprj_kh,(dtset%natom,mcprj_k))
1196 
1197  if(usepaw==1) then
1198    call pawcprj_alloc(cprj_k,0,dimcprj)
1199    call pawcprj_alloc(cprj_kh,0,dimcprj)
1200  end if
1201  ABI_MALLOC(smn,(2,nbdmax,nbdmax))
1202  ABI_MALLOC(mmn,(2,nbdmax,nbdmax))
1203 
1204 !Explanation for the index for the wavefunction stored in scf_history_wf
1205 !The reference is the cg+cprj output after the wf optimization at istep 1.
1206 !For wavefunction mixing for molecular dynamics, we use the same mixing as for the density in extraprho. To keep the same indexes,
1207 ! we choose to take indh=3 for the reference.
1208 
1209 !First step
1210  if (istep==1) then
1211 
1212    indh=3   ! This input wavefunction is the reference
1213 
1214 !  LOOP OVER SPINS
1215    do isppol=1,dtset%nsppol
1216 
1217 !    BIG FAT k POINT LOOP
1218      do ikpt=1,dtset%nkpt
1219 
1220 !      Select k point to be treated by this proc
1221        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1222        nbdmix=min(nband_k,nbdmax)
1223 
1224        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
1225 
1226        npw_k=npwarr(ikpt)
1227 
1228        scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
1229 
1230        if(usepaw==1) then
1231 !        scf_history_wf%cprj(:,ibg_hist+1:ibg_hist+my_nspinor*nbdmix,1)=cprj(:,ibg+1:ibg+my_nspinor*nbdmix)
1232          call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,0,isppol,mband,&
1233 &         dtset%mkmem,dtset%natom,nbdmax,nbdmix,my_nspinor,dtset%nsppol,0,&
1234 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1235          call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,0,isppol,&
1236 &         mband,dtset%mkmem,dtset%natom,nbdmax,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
1237 &         mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1238        end if
1239 
1240 !      Update the counters
1241        ibg=ibg+my_nspinor*nband_k
1242        ibg_hist=ibg_hist+my_nspinor*nbdmix
1243        icg=icg+my_nspinor*nband_k*npw_k
1244        icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
1245 
1246      end do
1247    end do
1248 
1249  else
1250 !  From istep==2
1251    if (istep==2) ind1=3
1252    if (istep==3) ind2=3
1253 !  biorthogonalization
1254    indh=3   ! This input wavefunction is the reference
1255 
1256 !  LOOP OVER SPINS
1257    do isppol=1,dtset%nsppol
1258 
1259 !    BIG FAT k POINT LOOP
1260      do ikpt=1,dtset%nkpt
1261 
1262 !      Select k point to be treated by this proc
1263        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1264        nbdmix=min(nband_k,nbdmax)
1265        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
1266        istwf_k=dtset%istwfk(ikpt)
1267        npw_k=npwarr(ikpt)
1268        ABI_MALLOC(psi_ortho,(2,npw_k*my_nspinor*nbdmix))
1269        psi_ortho=zero
1270 !      Biorthogonalization
1271 
1272        if(usepaw==1) then
1273          call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,0,isppol,mband,&
1274 &         dtset%mkmem,dtset%natom,nbdmax,nbdmix,my_nspinor,dtset%nsppol,0,&
1275 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1276          call pawcprj_get(atindx1,cprj_kh,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,0,isppol,&
1277 &         mband,dtset%mkmem,dtset%natom,nbdmax,nbdmix,my_nspinor,dtset%nsppol,0,&
1278 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1279        end if  !end usepaw=1
1280 
1281        hermitian=0
1282 
1283        call dotprod_set_cgcprj(atindx1,scf_history_wf%cg(:,:,indh),cg,cprj_kh,cprj_k,dimcprj,hermitian,&
1284 &         0,0,icg_hist,icg,ikpt,isppol,istwf_k,nbdmax,mcg,mcg,mcprj_k,mcprj_k,dtset%mkmem,&
1285 &         mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn(:,1:nbdmix,1:nbdmix),usepaw)
1286 
1287 !      Invert S matrix, that is NOT hermitian.
1288 !      Calculate M=S^-1
1289        mmn=zero
1290        do kk=1,nbdmix
1291          mmn(1,kk,kk)=one
1292        end do
1293 
1294        ABI_MALLOC(ipiv,(nbdmix))
1295 !      The smn is destroyed by the following inverse call
1296        call zgesv(nbdmix,nbdmix,smn,nbdmax,ipiv,mmn,nbdmax,ierr)
1297        ABI_FREE(ipiv)
1298 !DEBUG
1299        if(ierr/=0)then
1300          ABI_ERROR(' The call to cgesv general inversion routine failed')
1301        end if
1302 !ENDDEBUG
1303 
1304 !      The M matrix is used to compute the biorthogonalized set of wavefunctions, and to store it at the proper place
1305        inplace=0
1306        call lincom_cgcprj(mmn(:,1:nbdmix,1:nbdmix),cg,cprj_k,dimcprj,&
1307 &         icg,inplace,mcg,mcprj_k,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw,&
1308 &         cgout=psi_ortho,cprjout=cprj_kh,icgout=0)
1309 
1310 !!!TEST
1311 !      if (usepaw==0) then
1312 !        do iband=1,nband_k
1313 !          call dotprod_g(dotr,doti,istwf_k,npw_k,2,scf_history_wf%cg(:,icg+1+my_nspinor*npw_k:icg+2*my_nspinor*npw_k,indh),&
1314 !&            psi_ortho(:,1+(iband-1)*my_nspinor*npw_k:iband*my_nspinor*npw_k),mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1315 !          write(80+mpi_enreg%me,*) dotr,doti
1316 !          flush(80+mpi_enreg%me)
1317 !        end do
1318 !      else
1319 !        hermitian=0
1320 !        call dotprod_set_cgcprj(atindx1,scf_history_wf%cg(:,:,indh),psi_ortho,scf_history_wf%cprj(:,:,indh)&
1321 !&         ,cprj_kh,dimcprj,hermitian,0,0,icg_hist,icg_hist,ikpt,isppol,istwf_k,nbdmax,mcg,mcg,mcprj_k,mcprj_k,dtset%mkmem,&
1322 !&         mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn(:,1:nbdmix,1:nbdmix),usepaw)
1323 !        write(90+mpi_enreg%me,*) smn
1324 !        flush(90+mpi_enreg%me)
1325 !      end if
1326 !!!TEST
1327 !      The biorthogonalised set of wavefunctions is now stored at the proper place
1328 
1329 !       cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=zero
1330 
1331 !      psi(t+dt) <- psi(t) + alpha.psi(t)
1332        cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=(one+alpha)*psi_ortho(:,1:my_nspinor*npw_k*nbdmix)
1333        if(usepaw==1) then
1334          do ibdmix=1,nbdmix
1335            call pawcprj_axpby(one+alpha,zero,cprj_kh(:,ibdmix:ibdmix),cprj_k(:,ibdmix:ibdmix))
1336          end do ! end loop on ibdmix
1337        end if
1338 !      psi(t+dt) <- -alpha.psi(t-dt) + beta.psi(t-dt)
1339        if (abs(beta-alpha)>tol14.and.ind1>0) then
1340          cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)+&
1341 &             (beta-alpha)*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind1)
1342          if(usepaw==1) then
1343            do ibdmix=1,nbdmix
1344              call pawcprj_axpby(beta-alpha,one,scf_history_wf%cprj(:,ibg_hist+ibdmix:ibg_hist+ibdmix,ind1),&
1345 &                                 cprj_k(:,ibdmix:ibdmix))
1346            end do ! end loop on ibdmix
1347          end if
1348        end if
1349 
1350 !      psi(t+dt) <- -beta.psi(t-2dt)
1351        if (abs(beta)>tol14.and.ind2>0) then
1352          cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)&
1353 &                               -beta*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind2)
1354          if(usepaw==1) then
1355            do ibdmix=1,nbdmix
1356              call pawcprj_axpby(-beta,one,scf_history_wf%cprj(:,ibg_hist+ibdmix:ibg_hist+ibdmix,ind2),cprj_k(:,ibdmix:ibdmix))
1357            end do ! end loop on ibdmix
1358          end if
1359        end if
1360 
1361 !      Store psi(t) in history
1362        scf_history_wf%cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix,ind1new)=psi_ortho(:,1:my_nspinor*npw_k*nbdmix)
1363        if(usepaw==1) then
1364          call pawcprj_put(atindx1,cprj_kh,scf_history_wf%cprj(:,:,ind1new),dtset%natom,1,ibg_hist,ikpt,0,isppol,&
1365 &         nbdmix,dtset%mkmem,dtset%natom,nbdmax,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
1366 &         mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1367        end if
1368 
1369 !      Back to usual orthonormalization for the cg and cprj_k
1370        call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,mcprj_k,dtset%mkmem,&
1371 &       mpi_enreg,dtset%natom,nattyp,nbdmax,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw)
1372 
1373 !      Need to transfer cprj_k to cprj
1374        if(usepaw==1) then
1375          call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,0,isppol,&
1376 &         mband,dtset%mkmem,dtset%natom,nbdmax,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
1377 &         mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1378        end if
1379 
1380        ibg=ibg+my_nspinor*nband_k
1381        ibg_hist=ibg_hist+my_nspinor*nbdmix
1382        icg=icg+my_nspinor*nband_k*npw_k
1383        icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
1384        ABI_FREE(psi_ortho)
1385 !      End big k point loop
1386      end do
1387 !    End loop over spins
1388    end do
1389 
1390 
1391  end if !istep>1
1392 
1393 
1394  if(usepaw==1) then
1395    call pawcprj_free(cprj_k)
1396    call pawcprj_free(cprj_kh)
1397  end if
1398  ABI_FREE(cprj_k)
1399  ABI_FREE(cprj_kh)
1400  ABI_FREE(dimcprj)
1401  ABI_FREE(mmn)
1402  ABI_FREE(smn)
1403 
1404 
1405 
1406 end subroutine extrapwf_biortho

ABINIT/m_extraprho [ Modules ]

[ Top ] [ Modules ]

NAME

  m_extraprho

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2022 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 .

SOURCE

15 #if defined HAVE_CONFIG_H
16 #include "config.h"
17 #endif
18 
19 #include "abi_common.h"
20 
21 module m_extraprho
22 
23  use defs_basis
24  use m_abicore
25  use m_scf_history
26  use m_errors
27  use m_xmpi
28  use m_cgtools
29  use m_dtset
30 
31  use defs_datatypes, only : pseudopotential_type
32  use defs_abitypes, only : MPI_type
33  use m_atomdata, only : atom_length
34  use m_numeric_tools,   only : hermit
35  use m_geometry, only : metric
36  use m_kg, only : getph
37  use m_jellium,  only : jellium
38  use m_atm2fft,  only : atm2fft
39  use m_pawtab,   only : pawtab_type
40  use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_inquire_dim, pawrhoij_filter
41  use m_pawcprj,  only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_get, pawcprj_lincom, &
42                         pawcprj_free, pawcprj_zaxpby,pawcprj_axpby, pawcprj_put, pawcprj_getdim
43  use m_mpinfo,   only : proc_distrb_cycle
44  use m_cgprj,    only : ctocprj
45 
46  implicit none
47 
48  private