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

PARENTS

      scfcv

CHILDREN

      atm2fft,extrapwf,jellium,pawrhoij_alloc,wrtout

SOURCE

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

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

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

  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

1175 subroutine extrapwf_biortho(atindx,atindx1,cg,dtset,istep,kg,mcg,mgfft,mpi_enreg,&
1176 & nattyp,ngfft,npwarr,pawtab,psps,rprimd,scf_history,xred_old,ylm)
1177 
1178 
1179 !This section has been created automatically by the script Abilint (TD).
1180 !Do not modify the following lines by hand.
1181 #undef ABI_FUNC
1182 #define ABI_FUNC 'extrapwf_biortho'
1183 !End of the abilint section
1184 
1185  implicit none
1186 
1187 !Arguments ------------------------------------
1188 !scalars
1189  integer,intent(in) :: istep,mcg,mgfft
1190  type(MPI_type),intent(in) :: mpi_enreg
1191  type(dataset_type),intent(in) :: dtset
1192  type(scf_history_type),intent(inout) :: scf_history
1193  type(pseudopotential_type),intent(in) :: psps
1194 !arrays
1195  integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem),nattyp(dtset%ntypat),ngfft(18)
1196  integer,intent(in) :: npwarr(dtset%nkpt)
1197  real(dp),intent(in) :: rprimd(3,3)
1198  real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
1199  real(dp), intent(inout) :: cg(2,mcg)
1200  real(dp),intent(in) :: xred_old(3,dtset%natom)
1201  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
1202 
1203 !Local variables-------------------------------
1204 !scalars
1205  integer :: ia,iat,iatom,iband_max,iband_max1,iband_min,iband_min1,ibd,ibg,iblockbd,iblockbd1,icg,icgb,icgb1,icgb2
1206  integer :: ierr,ig,ii,ikpt,ilmn1,ilmn2,inc,ind1,ind2,iorder_cprj
1207  integer :: isize,isppol,istep1,istwf_k,itypat,klmn,me_distrb,my_nspinor
1208  integer :: nband_k,nblockbd,nprocband,npw_k,npw_nk,ntypat,option,spaceComm_band, usepaw
1209  real(dp) :: dotr,dotr1,doti,doti1,eigval
1210  !character(len=500) :: message
1211 !arrays
1212  real(dp) :: alpha(2),beta(2),gmet(3,3),gprimd(3,3),rmet(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom),ucvol
1213  integer,allocatable :: bufsize(:),bufsize_wf(:),bufdisp(:),bufdisp_wf(:),dimcprj(:),npw_block(:),npw_disp(:)
1214  real(dp),allocatable :: al(:,:),anm(:),cwavef(:,:),cwavef1(:,:),cwavef_tmp(:,:),deltawf1(:,:),deltawf2(:,:)
1215  real(dp),allocatable :: eig(:),evec(:,:)
1216  real(dp),allocatable :: unm(:,:,:)
1217  real(dp),allocatable :: work(:,:),work1(:,:),wf1(:,:),ylmgr_k(:,:,:),zhpev1(:,:),zhpev2(:)
1218  complex(dpc),allocatable :: unm_tmp(:,:),anm_tmp(:,:)
1219  type(pawcprj_type),allocatable :: cprj(:,:),cprj_k(:,:),cprj_k1(:,:),cprj_k2(:,:),cprj_k3(:,:),cprj_k4(:,:)
1220 !complex(dpc) :: aa
1221 
1222 ! *************************************************************************
1223 
1224  if (istep==0) return
1225  usepaw=dtset%usepaw
1226  ntypat=dtset%ntypat
1227  option =1
1228 !Useful array
1229  if (usepaw==1) then
1230    ABI_ALLOCATE(dimcprj,(dtset%natom))
1231    call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O')
1232  end if
1233 
1234 !Metric
1235  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1236 
1237 !History indexes
1238  ind1=1;ind2=2
1239 
1240 !First step
1241  if (istep==1) then
1242    scf_history%cg(:,:,ind1)=cg(:,:)
1243 !  scf_history%cg(:,:,ind2)=zero
1244    scf_history%cg(:,:,ind2)= cg(:,:)
1245    if(usepaw==1) then
1246 !    WARNING: THIS SECTION IS USELESS; NOW crpj CAN BE READ FROM SCFCV
1247      call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred_old)
1248      iatom=0 ; iorder_cprj=0
1249      call pawcprj_alloc(scf_history%cprj(:,:,ind1),0,dimcprj)
1250      call pawcprj_alloc(scf_history%cprj(:,:,ind2),0,dimcprj)
1251      ABI_ALLOCATE(ylmgr_k,(dtset%mpw,3,0))
1252      call ctocprj(atindx,cg,1,scf_history%cprj(:,:,ind1),gmet,gprimd,&
1253 &     iatom,0,iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,scf_history%mcprj,&
1254 &     dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,&
1255 &     dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,&
1256 &     dtset%nloalg,npwarr,dtset%nspinor,dtset%nsppol,dtset%ntypat,&
1257 &     dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,0,&
1258 &     xred_old,ylm,ylmgr_k)
1259      ABI_DEALLOCATE(ylmgr_k)
1260 !    call pawcprj_set_zero(scf_history%cprj(:,:,ind2))
1261      call pawcprj_copy(scf_history%cprj(:,:,ind1),scf_history%cprj(:,:,ind2))
1262    end if
1263  end if
1264 
1265 
1266 ! call wf_mixing(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,&
1267 !&         nattyp,npwarr,pawtab,scf_history,option)
1268 
1269 
1270 end subroutine extrapwf_biortho

ABINIT/m_extraprho [ Modules ]

[ Top ] [ Modules ]

NAME

  m_extraprho

FUNCTION

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 .

PARENTS

CHILDREN

SOURCE

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