Line data Source code
1 : !!****m* ABINIT/m_extraprho
2 : !! NAME
3 : !! m_extraprho
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1998-2026 ABINIT group (MT, FJ)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
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 m_extfpmd, only : extfpmd_type
32 : use defs_datatypes, only : pseudopotential_type
33 : use defs_abitypes, only : MPI_type
34 : use m_atomdata, only : atom_length
35 : use m_numeric_tools, only : hermit
36 : use m_geometry, only : metric
37 : use m_kg, only : getph
38 : use m_jellium, only : jellium
39 : use m_atm2fft, only : atm2fft
40 : use m_pawtab, only : pawtab_type
41 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_inquire_dim, pawrhoij_filter
42 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_get, pawcprj_lincom, &
43 : pawcprj_free, pawcprj_zaxpby,pawcprj_axpby, pawcprj_put, pawcprj_getdim
44 : use m_mpinfo, only : proc_distrb_cycle
45 : use m_cgprj, only : ctocprj
46 :
47 : implicit none
48 :
49 : private
50 : !!***
51 :
52 : public :: extraprho
53 : !!***
54 :
55 : contains
56 : !!***
57 :
58 : !!****f* ABINIT/extraprho
59 : !!
60 : !! NAME
61 : !! extraprho
62 : !!
63 : !! FUNCTION
64 : !! Extrapolate electronic density for new ionic positions
65 : !! from values of density of previous SCF cycle.
66 : !! Use algorithm proposed by D. Alfe in Comp. Phys. Comm. 118 (1999), 31-33 [[cite:Alfe1999]]
67 : !!
68 : !! INPUTS
69 : !! atindx
70 : !! atindx1(natom)=index table for atoms, inverse of atindx
71 : !! cg(2,mcg)= plane wave wavefunction coefficient
72 : !! cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn>
73 : !! dtset <type(dataset_type)>=all input variables in this dataset
74 : !! | densty(ntypat,4)=parameters for initialisation of the gaussian density
75 : !! | jellslab,slabzbeg,slabzend,slabwsrad=parameters for jellium slab
76 : !! | natom=number of atoms in cell.
77 : !! | nspden=number of spin-density components
78 : !! gmet(3,3)=reciprocal space metric
79 : !! gprimd(3,3)=reciprocal space dimensional primitive translations
80 : !! gsqcut=cutoff value on G**2 for sphere inside fft box
81 : !! istep=number of call the routine
82 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
83 : !! mgfft=maximum size of 1D FFTs
84 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
85 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
86 : !! mpi_enreg=information about MPI parallelization
87 : !! mqgrid=number of grid pts in q array for f(q) spline.
88 : !! my_natom=number of atoms treated by current processor
89 : !! nattyp(ntypat)=number of atoms of each type in cell.
90 : !! nfft=(effective) number of FFT grid points (for this processor)
91 : !! ngfft(18)=contain all needed information about 3D FFT
92 : !! npwarr=(nkpt)=number of planewaves in basis at this k point
93 : !! ntypat=number of types of atoms in cell
94 : !! pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
95 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
96 : !! psps<type(pseudopotential_type)>=variables related to pseudopotentials
97 : !! qgrid(mqgrid)=q grid for spline from 0 to qmax
98 : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
99 : !! ucvol=unit cell volume (bohr**3).
100 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
101 : !! xred_new(3,natom)=new reduced coordinates for atoms in unit cell
102 : !! xred_old(3,natom)=old reduced coordinates for atoms in unit cell
103 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
104 : !! zion(ntypat)=charge on each type of atom
105 : !! znucl(ntypat)=atomic number of each atom type
106 : !!
107 : !! SIDE EFFECTS
108 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= PAW rhoij occupancies and related data
109 : !! Value from previous SCF cycle is input
110 : !! Extrapolated value is output
111 : !! rhor(nfft,nspden)=the density from previous SCF cycle is input
112 : !! the extrapolated density is output
113 : !! scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles
114 : !!
115 : !! SOURCE
116 :
117 684 : subroutine extraprho(atindx,atindx1,cg,cprj,dtset,gmet,gprimd,gsqcut,istep,&
118 342 : & kg,mcg,mcprj,mgfft,mpi_enreg,mqgrid,my_natom,nattyp,nfft,ngfft,npwarr,ntypat,pawrhoij,&
119 342 : & pawtab,ph1d,psps,qgrid,rhor,rprimd,scf_history,ucvol,usepaw,&
120 342 : & xred_new,xred_old,ylm,zion,znucl,extfpmd)
121 :
122 : !Arguments ------------------------------------
123 : !scalars
124 : integer,intent(in) :: istep,mcg,mcprj,mgfft,my_natom,mqgrid,nfft,ntypat,usepaw
125 : real(dp),intent(in) :: gsqcut,ucvol
126 : type(MPI_type),intent(in) :: mpi_enreg
127 : type(dataset_type),intent(in) :: dtset
128 : type(scf_history_type),intent(inout) :: scf_history
129 : type(pseudopotential_type),intent(in) :: psps
130 : type(extfpmd_type),intent(in),optional,pointer :: extfpmd
131 : !arrays
132 : integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem)
133 : integer,intent(in) :: nattyp(ntypat),ngfft(18),npwarr(dtset%nkpt)
134 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom)
135 : real(dp),intent(in) :: qgrid(mqgrid)
136 : real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
137 : real(dp),intent(in) :: zion(ntypat),znucl(ntypat)
138 : real(dp), intent(inout) :: cg(2,mcg)
139 : real(dp),intent(inout) :: rhor(nfft,dtset%nspden),rprimd(3,3),xred_new(3,dtset%natom)
140 : real(dp),intent(in) :: xred_old(3,dtset%natom)
141 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw)
142 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
143 : type(pawcprj_type),intent(inout) :: cprj(:,:)
144 :
145 : !Local variables-------------------------------
146 : !scalars
147 : integer :: cplex_rhoij,dplex,iatom,ii,ind1,ind1new,ind2,ind2new,iq,iq0,irhoij,ispden,itypat,jrhoij,klmn
148 : integer :: lmn2_size,nselect,nspden_rhoij,optatm,optdyfr,opteltfr,optgr,option,optn,optn2
149 : integer :: optstr,optv,qphase_rhoij
150 : real(dp) :: a11,a12,a22,a33,alpha,b1,b2,beta,detA,fact,ratio1,ratio2
151 : logical :: hasmoved,usegauss
152 : character(len=500) :: message
153 : !arrays
154 : integer :: dummy3(3)
155 : real(dp) :: diff_t(3),diff_tmdt(3),diff_tpdt(3),dummy2(2)
156 : real(dp) :: dummy_in(0)
157 : real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0)
158 : real(dp) :: strn_dummy6(6),strv_dummy6(6)
159 342 : real(dp),allocatable :: deltarho(:),gauss(:,:),rhoijtmp(:,:),work1(:)
160 342 : real(dp),allocatable :: work2(:,:),work3(:,:),xred_tpdt(:,:)
161 :
162 : ! *************************************************************************
163 :
164 : !---------------------------------------------------------------
165 : !----------- Inits
166 : !---------------------------------------------------------------
167 :
168 : !History indexes
169 342 : ind1=scf_history%hindex(1)
170 342 : ind2=scf_history%hindex(2)
171 :
172 : !Compatibility tests
173 342 : if (ind1==0.and.ind2>0)then
174 0 : ABI_BUG(' Incompatible history indexes !')
175 : end if
176 :
177 : !Rotated values of history indexes
178 342 : if (ind1>0.and.ind2>0) then
179 : ind1new=ind2;ind2new=ind1
180 191 : else if (ind1>0.and.ind2==0) then
181 93 : ind1new=3-ind1;ind2new=ind1
182 98 : else if (ind1==0.and.ind2==0) then
183 98 : ind1new=1;ind2new=0
184 : end if
185 :
186 : !Compute ionic positions at t+dt in red. coordinates
187 : !Has to take the boundary conditions into account
188 1026 : ABI_MALLOC(xred_tpdt,(3,dtset%natom))
189 4127 : do iatom=1,dtset%natom
190 3785 : xred_tpdt(1,iatom)=xred_old(1,iatom)+mod(xred_new(1,iatom)-xred_old(1,iatom)+1.5_dp,one)-half
191 3785 : xred_tpdt(2,iatom)=xred_old(2,iatom)+mod(xred_new(2,iatom)-xred_old(2,iatom)+1.5_dp,one)-half
192 4127 : xred_tpdt(3,iatom)=xred_old(3,iatom)+mod(xred_new(3,iatom)-xred_old(3,iatom)+1.5_dp,one)-half
193 : end do
194 :
195 : !---------------------------------------------------------------
196 : !----------- Compute Alpha and Beta
197 : !----------- see (4) in Comp. Phys. Comm. 118 (1999), 31-33 [[cite:Alfe1999]]
198 : !---------------------------------------------------------------
199 :
200 : !Compute a_ij matrix
201 342 : a11=zero;a12=zero;a22=zero;a33=zero;b1=zero;b2=zero
202 342 : diff_t=zero;diff_tmdt=zero;diff_tpdt=zero
203 4127 : do iatom=1,dtset%natom
204 :
205 15140 : diff_tpdt(1:3)=xred_tpdt(1:3,iatom)-xred_old(1:3,iatom)
206 3785 : if (ind1>0) then
207 13976 : diff_t(1:3)=scf_history%xreddiff(1:3,iatom,ind1)
208 13202 : if (ind2>0) diff_tmdt(1:3)=scf_history%xreddiff(1:3,iatom,ind2)
209 : end if
210 15140 : do ii=1,3
211 11355 : a11=a11+diff_t(ii)**2
212 11355 : a22=a22+diff_tmdt(ii)**2
213 11355 : a33=a33+diff_tpdt(ii)**2
214 11355 : a12=a12+diff_t(ii) *diff_tmdt(ii)
215 11355 : b1 =b1 +diff_t(ii) *diff_tpdt(ii)
216 15140 : b2 =b2 +diff_tmdt(ii)*diff_tpdt(ii)
217 : end do
218 :
219 : ! Store reduced coordinates diffs in SCF history
220 15482 : scf_history%xreddiff(1:3,iatom,ind1new)=diff_tpdt(1:3)
221 :
222 : end do
223 342 : ABI_FREE(xred_tpdt)
224 342 : hasmoved=(a11>=tol10.or.a22>=tol10.or.a33>=tol10)
225 :
226 : !Compute alpha and beta
227 342 : alpha=zero;beta=zero
228 342 : if (hasmoved.and.ind1>0) then
229 244 : ratio1=one;if (abs(a33)>=tol10) ratio1=(a11+a33-two*b1)/a33
230 244 : ratio2=one;if (abs(a33)>=tol10) ratio2=(a11+a33-two*b2)/a33
231 244 : detA=a11*a22-a12**2
232 244 : if (abs(a11)>=tol10.and.(abs(a22)<tol10.or.abs(detA)<tol10)) then
233 147 : alpha=b1/a11
234 97 : else if (abs(a22)>=tol10.and.(abs(a11)<tol10.or.abs(detA)<tol10)) then
235 0 : beta=b2/a22
236 97 : else if (abs(ratio1)+abs(ratio2)<tol6) then
237 0 : if (ind2>0) then
238 0 : alpha=two;beta=-one
239 : else
240 0 : alpha=one
241 : end if
242 : write(message,'(3a,f4.1,a,f4.1)')&
243 0 : & 'Ionic positions lead to a collinear system !',ch10,&
244 0 : & 'Mixing coeffs have been set to: alpha=',alpha,' beta=',beta
245 0 : ABI_WARNING(message)
246 97 : else if (abs(a11)>=tol10.and.abs(a22)>=tol10) then
247 97 : alpha=(b1*a22-b2*a12)/detA
248 97 : beta =(b2*a11-b1*a12)/detA
249 : end if
250 : end if
251 :
252 :
253 : !---------------------------------------------------------------
254 : !----------- Contribution from delta_rho(t), delta_rho(t-dt)
255 : !----------- and delta_rho(t-2dt) to predicted rho(t+dt)
256 : !---------------------------------------------------------------
257 :
258 : !deltarho(t+dt) <- deltarho(t) + alpha.[deltarho(t)-deltarho(t-dt)]
259 : !+ beta .[deltarho(t-dt)-deltarho(t-2dt)]
260 : !Note: scf_history%deltarhor is updated at the same time
261 :
262 1026 : ABI_MALLOC(deltarho,(nfft))
263 685 : do ispden=1,dtset%nspden
264 :
265 343 : if (ispden==1) then
266 14606454 : deltarho(:)=rhor(:,ispden)-scf_history%atmrho_last(:)
267 1 : else if (ispden==2.and.dtset%nspden==2) then
268 13825 : deltarho(:)=rhor(:,ispden)-half*scf_history%atmrho_last(:)
269 : end if
270 :
271 :
272 : ! rho(t+dt) <- deltarho(t) + alpha.deltarho(t)
273 343 : if (dtset%nspden/=4.or.ispden==1) then
274 14620279 : rhor(:,ispden)=(one+alpha)*deltarho(:)
275 : else
276 0 : rhor(:,ispden)=(one+alpha)*rhor(:,ispden)
277 : end if
278 :
279 343 : if (hasmoved) then
280 :
281 : ! rho(t+dt) <- -alpha.deltarho(t-dt) + beta.deltarho(t-dt)
282 341 : if (abs(beta-alpha)>tol14.and.ind1>0) then
283 10148692 : rhor(:,ispden)=rhor(:,ispden)+(beta-alpha)*scf_history%deltarhor(:,ispden,ind1)
284 : end if
285 :
286 : ! rho(t+dt) <- -beta.deltarho(t-2dt)
287 341 : if (abs(beta)>tol14.and.ind2>0) then
288 793569 : rhor(:,ispden)=rhor(:,ispden)-beta*scf_history%deltarhor(:,ispden,ind2)
289 : end if
290 :
291 : end if
292 :
293 : ! Store deltarho(t) in history
294 685 : if (dtset%nspden/=4.or.ispden==1) then
295 14620279 : scf_history%deltarhor(:,ispden,ind1new)=deltarho(:)
296 : else
297 0 : scf_history%deltarhor(:,ispden,ind1new)=rhor(:,ispden)
298 : end if
299 :
300 : end do
301 :
302 342 : ABI_FREE(deltarho)
303 :
304 : !---------------------------------------------------------------
305 : !----------- Contribution from rho_at(t+dt) to predicted rho(t+dt)
306 : !---------------------------------------------------------------
307 :
308 : !Determine whether a gaussian atomic density has to be used or not
309 : !MG: Note that there's a small inconsistency between initro and extraprho because in initrho
310 : ! we use `use_gaussian(ntypat)`.
311 342 : usegauss=.true.
312 441 : if (usepaw==0) usegauss = any(.not.psps%nctab(1:ntypat)%has_tvale)
313 535 : if (usepaw==1) usegauss=(minval(pawtab(1:ntypat)%has_tvale)==0)
314 342 : if (usegauss) then
315 216 : optn2=3
316 648 : ABI_MALLOC(gauss,(2,ntypat))
317 432 : do itypat=1,ntypat
318 216 : gauss(1,itypat)=zion(itypat)
319 432 : gauss(2,itypat) = atom_length(dtset%densty(itypat,1),zion(itypat),znucl(itypat))
320 : end do
321 216 : call wrtout(std_out," Extrapolating rho(t+dt) using gaussian functions as atomic densities", "COLL")
322 : else
323 126 : optn2=2
324 126 : ABI_MALLOC(gauss,(2,0))
325 126 : call wrtout(std_out," Extrapolating rho(t+dt) using atomic densities taken from pseudos", "COLL")
326 : end if
327 :
328 : !Compute rho_at(t+dt) as sum of atomic densities
329 : !Note: scf_history%atmrho_last is updated at the same time
330 342 : optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optv=0;optn=1
331 : call atm2fft(atindx1,scf_history%atmrho_last,dummy_out1,dummy_out2,dummy_out3,&
332 : & dummy_out4,gauss,gmet,gprimd,dummy_out5,dummy_out6,gsqcut,mgfft,mqgrid,dtset%natom,nattyp,&
333 : & nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,qgrid,&
334 : & dummy3,dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,usepaw,dummy_in,dummy_in,dummy_in,dummy2,dummy_in,&
335 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
336 342 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
337 342 : ABI_FREE(gauss)
338 :
339 : !Take eventually into account jellium slab
340 342 : if (dtset%jellslab/=0) then
341 0 : option=2
342 0 : ABI_MALLOC(work1,(nfft))
343 0 : ABI_MALLOC(work2,(nfft,1))
344 0 : ABI_MALLOC(work3,(2,nfft))
345 0 : work2(:,1)=scf_history%atmrho_last(:)
346 : call jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,1,option,&
347 0 : & dtset%slabwsrad,work3,work2,rprimd,work1,dtset%slabzbeg,dtset%slabzend)
348 0 : scf_history%atmrho_last(:)=work2(:,1)
349 0 : ABI_FREE(work1)
350 0 : ABI_FREE(work2)
351 0 : ABI_FREE(work3)
352 : end if
353 :
354 : !Add rho_at(t+dt) to rho(t+dt)
355 14606454 : rhor(:,1)=rhor(:,1)+scf_history%atmrho_last(:)
356 14166 : if (dtset%nspden==2) rhor(:,2)=rhor(:,2)+half*scf_history%atmrho_last(:)
357 :
358 : !---------------------------------------------------------------
359 : !----------- Extrapolation of PAW rhoij occupancy matrixes
360 : !---------------------------------------------------------------
361 :
362 342 : if (usepaw==1) then
363 :
364 147 : if (ind2==0) then
365 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
366 122 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cpxocc=dtset%pawcpxocc)
367 : call pawrhoij_alloc(scf_history%pawrhoij(:,ind1new),cplex_rhoij,nspden_rhoij,&
368 : & dtset%nspinor,dtset%nsppol,dtset%typat,pawtab=pawtab,&
369 122 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
370 : end if
371 :
372 467 : do iatom=1,my_natom
373 :
374 320 : nspden_rhoij=pawrhoij(iatom)%nspden
375 320 : lmn2_size=pawrhoij(iatom)%lmn2_size
376 320 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij;dplex=cplex_rhoij-1
377 320 : qphase_rhoij=pawrhoij(iatom)%qphase
378 :
379 320 : if (hasmoved) then
380 1280 : ABI_MALLOC(rhoijtmp,(cplex_rhoij*qphase_rhoij*lmn2_size,nspden_rhoij))
381 9490 : rhoijtmp=zero
382 :
383 640 : do ispden=1,nspden_rhoij
384 960 : do iq=1,qphase_rhoij
385 320 : iq0=merge(0,cplex_rhoij*lmn2_size,iq==1)
386 :
387 : ! rhoij(t+dt) <- rhoij(t) + alpha.rhoij(t)
388 320 : fact=one+alpha
389 320 : jrhoij=1+iq0
390 9170 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
391 8850 : klmn=cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex+iq0
392 : rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) &
393 17700 : & +fact*pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
394 9170 : jrhoij=jrhoij+cplex_rhoij
395 : end do
396 :
397 : ! rhoij(t+dt) <- -alpha.rhoij(t-dt) + beta.rhoij(t-dt)
398 320 : if (abs(beta-alpha)>tol14.and.ind1>0) then
399 234 : fact=beta-alpha
400 234 : jrhoij=1+iq0
401 4248 : do irhoij=1,scf_history%pawrhoij(iatom,ind1)%nrhoijsel
402 4014 : klmn=cplex_rhoij*scf_history%pawrhoij(iatom,ind1)%rhoijselect(irhoij)-dplex+iq0
403 : rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) &
404 8028 : & +fact*scf_history%pawrhoij(iatom,ind1)%rhoijp(jrhoij:jrhoij+dplex,ispden)
405 4248 : jrhoij=jrhoij+cplex_rhoij
406 : end do
407 : end if
408 :
409 : ! rho(t+dt) <- -beta.rhoij(t-2dt)
410 640 : if (abs(beta)>tol14.and.ind2>0) then
411 6 : fact=-beta
412 6 : jrhoij=1+iq0
413 222 : do irhoij=1,scf_history%pawrhoij(iatom,ind2)%nrhoijsel
414 216 : klmn=cplex_rhoij*scf_history%pawrhoij(iatom,ind2)%rhoijselect(irhoij)-dplex+iq0
415 : rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) &
416 432 : & +fact*scf_history%pawrhoij(iatom,ind2)%rhoijp(jrhoij:jrhoij+dplex,ispden)
417 222 : jrhoij=jrhoij+cplex_rhoij
418 : end do
419 : end if
420 :
421 : end do ! iq
422 : end do !ispden
423 : end if !hasmoved
424 :
425 : ! Store rhoij(t) in history
426 : ! (cannot use pawrhoij_copy here because update for single atom)
427 320 : nselect=pawrhoij(iatom)%nrhoijsel
428 320 : scf_history%pawrhoij(iatom,ind1new)%nrhoijsel=nselect
429 9170 : scf_history%pawrhoij(iatom,ind1new)%rhoijselect(:)=0
430 9170 : scf_history%pawrhoij(iatom,ind1new)%rhoijselect(1:nselect)=pawrhoij(iatom)%rhoijselect(1:nselect)
431 : scf_history%pawrhoij(iatom,ind1new)%rhoijp(1:cplex_rhoij*nselect,1:nspden_rhoij)= &
432 9490 : & pawrhoij(iatom)%rhoijp(1:cplex_rhoij*nselect,1:nspden_rhoij)
433 :
434 : ! Select non-zero values of rhoij(t+dt)
435 467 : if (hasmoved) then
436 : call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,pawrhoij(iatom)%nrhoijsel,&
437 : & cplex_rhoij,qphase_rhoij,lmn2_size,nspden_rhoij,&
438 320 : & rhoij_input=rhoijtmp)
439 320 : ABI_FREE(rhoijtmp)
440 : end if
441 :
442 : end do !iatom
443 : end if !usepaw
444 :
445 342 : if(present(extfpmd)) then
446 342 : if(associated(extfpmd)) then
447 0 : scf_history%nextfpmd(ind1new)=extfpmd%nelect
448 0 : if(hasmoved) then
449 0 : extfpmd%nelect=(one+alpha)*extfpmd%nelect
450 0 : if (abs(beta-alpha)>tol14.and.ind1>0) then
451 0 : extfpmd%nelect=extfpmd%nelect+(beta-alpha)*scf_history%nextfpmd(ind1)
452 : endif
453 0 : if (abs(beta)>tol14.and.ind2>0) then
454 0 : extfpmd%nelect=extfpmd%nelect-beta*scf_history%nextfpmd(ind2)
455 : endif
456 : endif
457 : endif
458 : endif
459 :
460 342 : scf_history%alpha=alpha
461 342 : scf_history%beta=beta
462 :
463 :
464 :
465 : !---------------------------------------------------------------
466 : !----------- End
467 : !---------------------------------------------------------------
468 :
469 342 : if(scf_history%usecg==1) then
470 8 : if (hasmoved) then
471 8 : if (dtset%extrapwf==1) then
472 : call extrapwf(atindx,atindx1,cg,dtset,istep,kg,mcg,mgfft,mpi_enreg,nattyp,&
473 8 : & ngfft,npwarr,ntypat,pawtab,psps,rprimd,scf_history,usepaw,xred_old,ylm)
474 0 : elseif(dtset%extrapwf==2) then
475 0 : scf_history%hindex(1)=ind1
476 0 : scf_history%hindex(2)=ind2
477 0 : scf_history%hindex(3)=ind1new
478 : call extrapwf_biortho(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,&
479 0 : & nattyp,npwarr,pawtab,scf_history)
480 : end if
481 : else
482 0 : scf_history%cg(:,:,2)=zero
483 : end if
484 :
485 : end if
486 : !Rotate history indexes
487 342 : scf_history%hindex(1)=ind1new
488 342 : scf_history%hindex(2)=ind2new
489 :
490 :
491 342 : end subroutine extraprho
492 : !!***
493 :
494 : !!****f* ABINIT/extrapwf
495 : !!
496 : !! NAME
497 : !! extrapwf
498 : !!
499 : !! FUNCTION
500 : !! Extrapolate wavefunctions for new ionic positions
501 : !! from values of wavefunctions of previous SCF cycle.
502 : !! Use algorithm proposed by T. A. Arias et al. in PRB 45, 1538 (1992) [[cite:Arias1992]]
503 : !!
504 : !! INPUTS
505 : !! atindx(natom)=index table for atoms
506 : !! atindx1(natom)=index table for atoms, inverse of atindx
507 : !! dtset <type(dataset_type)>=all input variables in this dataset
508 : !! istep=number of call the routine
509 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
510 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
511 : !! mgfft=maximum size of 1D FFTs
512 : !! mpi_enreg=information about MPI parallelization
513 : !! nattyp(ntypat)=number of atoms of each type in cell.
514 : !! ngfft(18)=contain all needed information about 3D FFT
515 : !! npwarr(nkpt)=number of planewaves in basis at this k point
516 : !! ntypat=number of types of atoms in cell
517 : !! pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
518 : !! psps<type(pseudopotential_type)>=variables related to pseudopotentials
519 : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
520 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
521 : !! xred_old(3,natom)=old reduced coordinates for atoms in unit cell
522 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
523 : !!
524 : !! SIDE EFFECTS
525 : !! cg(2,mcg)= plane wave wavefunction coefficient
526 : !! Value from previous SCF cycle is input
527 : !! Extrapolated value is output
528 : !! scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles
529 : !!
530 : !! NOTES
531 : !! THIS ROUTINE IS NOT USEABLE AT PRESENT.
532 : !! SHOULD BE CAREFULY TESTED AND DEBUGGED (ESPECIALLY WITHIN PAW).
533 : !!
534 : !! SOURCE
535 :
536 8 : subroutine extrapwf(atindx,atindx1,cg,dtset,istep,kg,mcg,mgfft,mpi_enreg,&
537 8 : & nattyp,ngfft,npwarr,ntypat,pawtab,psps,rprimd,scf_history,usepaw,xred_old,ylm)
538 :
539 : !Arguments ------------------------------------
540 : !scalars
541 : integer,intent(in) :: istep,mcg,mgfft,ntypat,usepaw
542 : type(MPI_type),intent(in) :: mpi_enreg
543 : type(dataset_type),intent(in) :: dtset
544 : type(scf_history_type),intent(inout) :: scf_history
545 : type(pseudopotential_type),intent(in) :: psps
546 : !arrays
547 : integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem),nattyp(ntypat),ngfft(18)
548 : integer,intent(in) :: npwarr(dtset%nkpt)
549 : real(dp),intent(in) :: rprimd(3,3)
550 : real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
551 : real(dp), intent(inout) :: cg(2,mcg)
552 : real(dp),intent(in) :: xred_old(3,dtset%natom)
553 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
554 :
555 :
556 : !Local variables-------------------------------
557 : !scalars
558 : integer :: ia,iat,iatom,iband_max,iband_max1,iband_min,iband_min1,ibd,ibg,iblockbd,iblockbd1,icg,icgb,icgb1,icgb2
559 : integer :: ierr,ig,ii,ikpt,ilmn1,ilmn2,inc,ind1,ind2,iorder_cprj
560 : integer :: isize,isppol,istep1,istwf_k,itypat,klmn,me_distrb,my_nspinor
561 : integer :: nband_k,nblockbd,nprocband,npw_k,npw_nk,spaceComm_band
562 : real(dp) :: dotr,dotr1,doti,doti1,eigval
563 : !character(len=500) :: message
564 : !arrays
565 16 : real(dp) :: alpha(2),beta(2),gmet(3,3),gprimd(3,3),rmet(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom),ucvol
566 8 : integer,allocatable :: bufsize(:),bufsize_wf(:),bufdisp(:),bufdisp_wf(:),dimcprj(:),npw_block(:),npw_disp(:)
567 8 : real(dp),allocatable :: al(:,:),anm(:),cwavef(:,:),cwavef1(:,:),cwavef_tmp(:,:),deltawf1(:,:),deltawf2(:,:)
568 8 : real(dp),allocatable :: eig(:),evec(:,:)
569 8 : real(dp),allocatable :: unm(:,:,:)
570 8 : real(dp),allocatable :: work(:,:),work1(:,:),wf1(:,:),ylmgr_k(:,:,:),zhpev1(:,:),zhpev2(:)
571 8 : complex(dp),allocatable :: unm_tmp(:,:),anm_tmp(:,:)
572 8 : type(pawcprj_type),allocatable :: cprj(:,:),cprj_k(:,:),cprj_k1(:,:),cprj_k2(:,:),cprj_k3(:,:),cprj_k4(:,:)
573 : !complex(dp) :: aa
574 : ! *************************************************************************
575 :
576 8 : if (istep==0) return
577 :
578 : !Useful array
579 8 : if (usepaw==1) then
580 24 : ABI_MALLOC(dimcprj,(dtset%natom))
581 8 : call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O')
582 : end if
583 :
584 : !Metric
585 8 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
586 :
587 : !History indexes
588 8 : ind1=1;ind2=2
589 :
590 : !First step
591 8 : if (istep==1) then
592 60481 : scf_history%cg(:,:,ind1)=cg(:,:)
593 : ! scf_history%cg(:,:,ind2)=zero
594 60481 : scf_history%cg(:,:,ind2)= cg(:,:)
595 1 : if(usepaw==1) then
596 : ! WARNING: THIS SECTION IS USELESS; NOW crpj CAN BE READ FROM SCFCV
597 1 : call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred_old)
598 1 : iatom=0 ; iorder_cprj=0
599 1 : call pawcprj_alloc(scf_history%cprj(:,:,ind1),0,dimcprj)
600 1 : call pawcprj_alloc(scf_history%cprj(:,:,ind2),0,dimcprj)
601 2 : ABI_MALLOC(ylmgr_k,(dtset%mpw,3,0))
602 : call ctocprj(atindx,cg,1,scf_history%cprj(:,:,ind1),gmet,gprimd,&
603 : & iatom,0,iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,scf_history%mcprj,&
604 : & dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,&
605 : & dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,&
606 : & dtset%nloalg,npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,dtset%ntypat,&
607 : & dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,0,&
608 1 : & xred_old,ylm,ylmgr_k)
609 1 : ABI_FREE(ylmgr_k)
610 : ! call pawcprj_set_zero(scf_history%cprj(:,:,ind2))
611 1 : call pawcprj_copy(scf_history%cprj(:,:,ind1),scf_history%cprj(:,:,ind2))
612 : end if
613 : else
614 :
615 : !From 2nd step
616 :
617 : ! Init parallelism
618 7 : me_distrb=mpi_enreg%me_kpt
619 7 : if (mpi_enreg%paral_kgb==1.or.mpi_enreg%paralbd==1) then
620 0 : spaceComm_band=mpi_enreg%comm_band
621 0 : nprocband=mpi_enreg%nproc_band
622 : else
623 7 : spaceComm_band=xmpi_comm_self
624 7 : nprocband=1
625 : end if
626 :
627 : ! For the moment sequential part only
628 7 : nprocband=1
629 :
630 : ! Additional statements if band-fft parallelism
631 : if (nprocband>1) then
632 : ABI_MALLOC(npw_block,(nprocband))
633 : ABI_MALLOC(npw_disp,(nprocband))
634 : ABI_MALLOC(bufsize,(nprocband))
635 : ABI_MALLOC(bufdisp,(nprocband))
636 : ABI_MALLOC(bufsize_wf,(nprocband))
637 : ABI_MALLOC(bufdisp_wf,(nprocband))
638 : end if
639 :
640 7 : icg=0
641 7 : ibg=0
642 :
643 7 : if(usepaw==1) then
644 : ! WARNING: THIS SECTION IS USELESS; NOW cprj CAN BE READ FROM SCFCV
645 7 : call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred_old)
646 1036 : ABI_MALLOC(cprj,(dtset%natom,scf_history%mcprj))
647 7 : call pawcprj_alloc(cprj,0,dimcprj)
648 7 : iatom=0 ; iorder_cprj=0
649 14 : ABI_MALLOC(ylmgr_k,(dtset%mpw,3,0))
650 : call ctocprj(atindx,cg,1,cprj,gmet,gprimd,iatom,0,iorder_cprj,&
651 : & dtset%istwfk,kg,dtset%kptns,mcg,scf_history%mcprj,dtset%mgfft,&
652 : & dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,&
653 : & nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,dtset%nloalg,&
654 : & npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,&
655 : & ph1d,psps,rmet,dtset%typat,ucvol,0,xred_old,&
656 7 : & ylm,ylmgr_k)
657 7 : ABI_FREE(ylmgr_k)
658 : end if ! end usepaw=1
659 :
660 : ! LOOP OVER SPINS
661 7 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
662 14 : do isppol=1,dtset%nsppol
663 :
664 : ! BIG FAT k POINT LOOP
665 42 : do ikpt=1,dtset%nkpt
666 :
667 : ! Select k point to be treated by this proc
668 28 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
669 28 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
670 :
671 28 : istwf_k=dtset%istwfk(ikpt)
672 :
673 : ! Retrieve number of plane waves
674 28 : npw_k=npwarr(ikpt)
675 : if (nprocband>1) then
676 : ! Special treatment for band-fft //
677 : call xmpi_allgather(npw_k,npw_block,spaceComm_band,ierr)
678 : npw_nk=sum(npw_block);npw_disp(1)=0
679 : do ii=2,nprocband
680 : npw_disp(ii)=npw_disp(ii-1)+npw_block(ii-1)
681 : end do
682 : else
683 : npw_nk=npw_k
684 : end if
685 :
686 : ! Allocate arrays for a wave-function (or a block of WFs)
687 84 : ABI_MALLOC(cwavef,(2,npw_nk*my_nspinor))
688 56 : ABI_MALLOC(cwavef1,(2,npw_nk*my_nspinor))
689 : if (nprocband>1) then
690 : isize=2*my_nspinor;bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:)
691 : isize=2*my_nspinor*npw_k;bufsize_wf(:)=isize
692 : do ii=1,nprocband
693 : bufdisp_wf(ii)=(ii-1)*isize
694 : end do
695 : end if
696 :
697 : ! Subspace alignment
698 :
699 : ! Loop over bands or blocks of bands
700 28 : nblockbd=nband_k/nprocband
701 28 : icgb=icg
702 :
703 28 : if(usepaw==1) then
704 1120 : ABI_MALLOC( cprj_k,(dtset%natom,my_nspinor*nblockbd))
705 28 : call pawcprj_alloc(cprj_k,cprj(1,1)%ncpgr,dimcprj)
706 : call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,dtset%mband,&
707 : & dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,&
708 28 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
709 1092 : ABI_MALLOC( cprj_k1,(dtset%natom,my_nspinor*nblockbd))
710 28 : call pawcprj_alloc(cprj_k1,scf_history%cprj(1,1,ind1)%ncpgr,dimcprj)
711 : call pawcprj_get(atindx1,cprj_k1,scf_history%cprj(:,:,ind1),dtset%natom,1,ibg,ikpt,1,isppol,&
712 : & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,&
713 28 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
714 1092 : ABI_MALLOC( cprj_k2,(dtset%natom,my_nspinor*nblockbd))
715 28 : call pawcprj_alloc(cprj_k2,scf_history%cprj(1,1,ind2)%ncpgr,dimcprj)
716 : call pawcprj_get(atindx1,cprj_k2,scf_history%cprj(:,:,ind2),dtset%natom,1,ibg,ikpt,1,isppol,&
717 : & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,&
718 28 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
719 : end if !end usepaw=1
720 :
721 112 : ABI_MALLOC(unm,(2,nblockbd,nblockbd))
722 12460 : unm=zero
723 : icgb2=0
724 :
725 364 : do iblockbd=1,nblockbd
726 336 : iband_min=1+(iblockbd-1)*nprocband
727 336 : iband_max=iblockbd*nprocband
728 :
729 336 : if(xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then
730 336 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min,iband_max,isppol,me_distrb)) cycle
731 : end if
732 :
733 : ! Extract wavefunction information
734 : if (nprocband>1) then
735 : ! Special treatment for band-fft //
736 : ABI_MALLOC(cwavef_tmp,(2,npw_k*my_nspinor*nprocband))
737 : do ig=1,npw_k*my_nspinor*nprocband
738 : cwavef_tmp(1,ig)=cg(1,ig+icgb)
739 : cwavef_tmp(2,ig)=cg(2,ig+icgb)
740 : end do
741 : call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef,bufsize,bufdisp,spaceComm_band,ierr)
742 : ABI_FREE(cwavef_tmp)
743 : else
744 141204 : do ig=1,npw_k*my_nspinor
745 140868 : cwavef(1,ig)=cg(1,ig+icgb)
746 141204 : cwavef(2,ig)=cg(2,ig+icgb)
747 : end do
748 : end if
749 :
750 : icgb1=icg
751 :
752 4368 : do iblockbd1=1,nblockbd
753 4032 : iband_min1=1+(iblockbd1-1)*nprocband
754 4032 : iband_max1=iblockbd1*nprocband
755 :
756 4032 : if(xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then
757 4032 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min1,iband_max1,isppol,me_distrb)) cycle
758 : end if
759 :
760 : ! Extract wavefunction information
761 :
762 : if (nprocband>1) then
763 : ! Special treatment for band-fft //
764 : ABI_MALLOC(cwavef_tmp,(2,npw_k*my_nspinor*nprocband))
765 : do ig=1,npw_k*my_nspinor*nprocband
766 : cwavef_tmp(1,ig)=scf_history%cg(1,ig+icgb1,ind1)
767 : cwavef_tmp(2,ig)=scf_history%cg(2,ig+icgb1,ind1)
768 : end do
769 : call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef1,bufsize,bufdisp,spaceComm_band,ierr)
770 : ABI_FREE(cwavef_tmp)
771 : else
772 1694448 : do ig=1,npw_k*my_nspinor
773 1690416 : cwavef1(1,ig)=scf_history%cg(1,ig+icgb1,ind1)
774 1694448 : cwavef1(2,ig)=scf_history%cg(2,ig+icgb1,ind1)
775 : end do
776 : end if
777 :
778 : ! Calculate Unm=<psi_nk(t)|S|psi_mk(t-dt)>
779 4032 : call dotprod_g(dotr,doti,istwf_k,npw_k*my_nspinor,2,cwavef,cwavef1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
780 4032 : if(usepaw==1) then
781 : ia =0
782 8064 : do itypat=1,ntypat
783 12096 : do iat=1+ia,nattyp(itypat)+ia
784 76608 : do ilmn1=1,pawtab(itypat)%lmn_size
785 354816 : do ilmn2=1,ilmn1
786 290304 : klmn=((ilmn1-1)*ilmn1)/2+ilmn2
787 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)+&
788 290304 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2))
789 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)-&
790 354816 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2))
791 : end do
792 298368 : do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
793 225792 : klmn=((ilmn2-1)*ilmn2)/2+ilmn1
794 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)+&
795 225792 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2))
796 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)-&
797 290304 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2))
798 : end do
799 : end do
800 : end do
801 8064 : ia=ia+nattyp(itypat)
802 : end do
803 : end if
804 : ! unm(1,iblockbd,iblockbd1)=dotr
805 : ! unm(2,iblockbd,iblockbd1)=doti
806 4032 : unm(1,iblockbd1,iblockbd)=dotr
807 4032 : unm(2,iblockbd1,iblockbd)=doti
808 : ! End loop over bands iblockbd1
809 4368 : icgb1=icgb1+npw_k*my_nspinor*nprocband
810 :
811 : end do
812 :
813 : ! End loop over bands iblockbd
814 : icgb2=icgb2+npw_k*my_nspinor*nprocband
815 364 : icgb=icgb+npw_k*my_nspinor*nprocband
816 : end do
817 :
818 : ! write(std_out,*) 'UNM'
819 : ! do iblockbd=1,nblockbd
820 : ! write(std_out,11) (unm(1,iblockbd,iblockbd1),unm(2,iblockbd,iblockbd1),iblockbd1=1,nblockbd)
821 : ! end do
822 : ! 11 format(12(1x,f9.5),a)
823 : ! Compute A=tU^*U
824 112 : ABI_MALLOC(unm_tmp,(nblockbd,nblockbd))
825 84 : ABI_MALLOC(anm_tmp,(nblockbd,nblockbd))
826 84 : ABI_MALLOC(anm,(nblockbd*(nblockbd+1)))
827 4396 : unm_tmp(:,:)=cmplx(unm(1,:,:),unm(2,:,:),kind=dp)
828 : call zgemm('C','N',nblockbd,nblockbd,nblockbd,dcmplx(1._dp), unm_tmp,nblockbd, &
829 28 : & unm_tmp,nblockbd,dcmplx(0._dp),anm_tmp,nblockbd)
830 364 : do iblockbd=1,nblockbd
831 2548 : do iblockbd1=iblockbd,nblockbd
832 2184 : ii=iblockbd1*(iblockbd1-1)+2*(iblockbd-1)+1
833 2184 : anm(ii)=real(anm_tmp(iblockbd,iblockbd1))
834 2520 : anm(ii+1)=aimag(anm_tmp(iblockbd,iblockbd1))
835 : end do
836 : end do
837 28 : call hermit(anm,anm,ierr,nblockbd)
838 : ! aa=dcmplx(0._dp)
839 : ! do iblockbd=1,nblockbd
840 : ! aa=aa+conjg(unm_tmp(iblockbd,1))*unm_tmp(iblockbd,1)
841 : ! end do
842 : ! write(std_out,*) 'tU*U', aa
843 : ! write(std_out,*) 'ANM_tmp'
844 : ! do iblockbd=1,nblockbd
845 : ! write(std_out,11) (anm_tmp(iblockbd,iblockbd1),iblockbd1=1,nblockbd)
846 : ! end do
847 : ! write(std_out,*) 'ANM'
848 : ! do iblockbd=1,nblockbd*(nblockbd+1)
849 : ! write(std_out,11) anm(iblockbd)
850 : ! end do
851 :
852 : ! Diagonalize A
853 84 : ABI_MALLOC(eig,(nblockbd))
854 112 : ABI_MALLOC(evec,(2*nblockbd,nblockbd))
855 84 : ABI_MALLOC(zhpev1,(2,2*nblockbd-1))
856 84 : ABI_MALLOC(zhpev2,(3*nblockbd-2))
857 : call zhpev('V','U',nblockbd,anm,eig,evec,nblockbd,zhpev1,&
858 28 : & zhpev2,ierr)
859 28 : ABI_FREE(anm)
860 28 : ABI_FREE(zhpev1)
861 28 : ABI_FREE(zhpev2)
862 : ! aa=dcmplx(0._dp)
863 : ! do iblockbd=1,nblockbd
864 : ! aa=aa+anm_tmp(1,iblockbd)*cmplx(evec((2*iblockbd-1),1),evec(2*iblockbd,1),kind=dp)
865 : ! end do
866 : ! write(std_out,*) 'EIG', aa, eig(1)*evec(1,1),eig(1)*evec(2,1)
867 :
868 : ! Compute A'=evec*tU^/sqrt(eig)
869 : call zgemm('C','C',nblockbd,nblockbd,nblockbd,dcmplx(1._dp),evec,nblockbd, &
870 28 : & unm_tmp,nblockbd,dcmplx(0._dp),anm_tmp,nblockbd)
871 364 : do iblockbd=1,nblockbd
872 336 : eigval=dsqrt(eig(iblockbd))
873 4396 : do iblockbd1=1,nblockbd
874 4368 : anm_tmp(iblockbd,iblockbd1)=anm_tmp(iblockbd,iblockbd1)/eigval
875 : end do
876 : end do
877 :
878 : ! Compute tA^A'to come back to the initial subspace for the cg's
879 :
880 : call zgemm('N','N',nblockbd,nblockbd,nblockbd,dcmplx(1._dp),evec,nblockbd, &
881 28 : & anm_tmp,nblockbd,dcmplx(0._dp),unm_tmp,nblockbd)
882 4424 : anm_tmp=unm_tmp
883 : ! write(std_out,*) 'ANM_tmp'
884 : ! do iblockbd=1,nblockbd
885 : ! write(std_out,11) (anm_tmp(iblockbd,iblockbd1),iblockbd1=1,nblockbd)
886 : ! end do
887 :
888 : ! Wavefunction alignment (istwfk=1 ?)
889 84 : ABI_MALLOC(work,(2,npw_nk*my_nspinor*nblockbd))
890 84 : ABI_MALLOC(work1,(2,my_nspinor*nblockbd*npw_nk))
891 422632 : work1(:,:)=scf_history%cg(:,icg+1:icg+my_nspinor*nblockbd*npw_nk,ind1)
892 : call zgemm('N','N',npw_nk*my_nspinor,nblockbd,nblockbd,dcmplx(1._dp), &
893 : & work1,npw_nk*my_nspinor, &
894 28 : & anm_tmp,nblockbd,dcmplx(0._dp),work,npw_nk*my_nspinor)
895 422632 : scf_history%cg(:,1+icg:npw_nk*my_nspinor*nblockbd+icg,ind1)=work(:,:)
896 :
897 422632 : work1(:,:)=scf_history%cg(:,icg+1:icg+my_nspinor*nblockbd*npw_nk,ind2)
898 : call zgemm('N','N',npw_nk*my_nspinor,nblockbd,nblockbd,dcmplx(1._dp), &
899 : & work1,npw_nk*my_nspinor, &
900 28 : & anm_tmp,nblockbd,dcmplx(0._dp),work,npw_nk*my_nspinor)
901 422632 : scf_history%cg(:,1+icg:npw_nk*my_nspinor*nblockbd+icg,ind2)=work(:,:)
902 28 : ABI_FREE(work1)
903 : ! If paw, must also align cprj:
904 28 : if (usepaw==1) then
905 : ! New version (MT):
906 196 : ABI_MALLOC(cprj_k3,(dtset%natom,my_nspinor))
907 168 : ABI_MALLOC(cprj_k4,(dtset%natom,my_nspinor))
908 28 : call pawcprj_alloc(cprj_k3,cprj_k1(1,1)%ncpgr,dimcprj)
909 28 : call pawcprj_alloc(cprj_k4,cprj_k2(1,1)%ncpgr,dimcprj)
910 84 : ABI_MALLOC(al,(2,nblockbd))
911 364 : do iblockbd=1,nblockbd
912 336 : ii=(iblockbd-1)*my_nspinor
913 4368 : do iblockbd1=1,nblockbd
914 4032 : al(1,iblockbd1)=real (anm_tmp(iblockbd,iblockbd1))
915 4368 : al(2,iblockbd1)=aimag(anm_tmp(iblockbd,iblockbd1))
916 : end do
917 336 : call pawcprj_lincom(al,cprj_k1,cprj_k3,nblockbd)
918 336 : call pawcprj_lincom(al,cprj_k2,cprj_k4,nblockbd)
919 336 : call pawcprj_copy(cprj_k3,cprj_k1(:,ii+1:ii+my_nspinor))
920 364 : call pawcprj_copy(cprj_k4,cprj_k2(:,ii+1:ii+my_nspinor))
921 : end do
922 28 : ABI_FREE(al)
923 : ! Old version (FJ):
924 : ! allocate( cprj_k3(dtset%natom,my_nspinor*nblockbd))
925 : ! call pawcprj_alloc(cprj_k3,cprj_k1(1,1)%ncpgr,dimcprj)
926 : ! allocate( cprj_k4(dtset%natom,my_nspinor*nblockbd))
927 : ! call pawcprj_alloc(cprj_k4,cprj_k2(1,1)%ncpgr,dimcprj)
928 : ! beta(1)=one;beta(2)=zero
929 : ! do iblockbd=1,nblockbd*my_nspinor
930 : ! do iblockbd1=1,nblockbd*my_nspinor
931 : ! alpha(1)=real(anm_tmp(iblockbd,iblockbd1));alpha(2)=aimag(anm_tmp(iblockbd,iblockbd1))
932 : ! call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd1:iblockbd1),cprj_k3(:,iblockbd:iblockbd))
933 : ! call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd1:iblockbd1),cprj_k4(:,iblockbd:iblockbd))
934 : ! end do
935 : ! end do
936 : ! call pawcprj_copy(cprj_k3,cprj_k1)
937 : ! call pawcprj_copy(cprj_k4,cprj_k2)
938 :
939 28 : call pawcprj_free(cprj_k3)
940 28 : call pawcprj_free(cprj_k4)
941 84 : ABI_FREE(cprj_k3)
942 84 : ABI_FREE(cprj_k4)
943 : end if
944 28 : ABI_FREE(anm_tmp)
945 28 : ABI_FREE(unm_tmp)
946 28 : ABI_FREE(work)
947 :
948 : ! Wavefunction extrapolation
949 28 : ibd=0
950 28 : inc=npw_nk*my_nspinor
951 56 : ABI_MALLOC(deltawf2,(2,npw_nk*my_nspinor))
952 56 : ABI_MALLOC(wf1,(2,npw_nk*my_nspinor))
953 56 : ABI_MALLOC(deltawf1,(2,npw_nk*my_nspinor))
954 364 : do iblockbd=1,nblockbd
955 422940 : deltawf2(:,:)=scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind2)
956 422940 : wf1(:,:)=scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind1)
957 : ! wf1(2,1)=zero;deltawf2(2,1)=zero
958 :
959 : 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),&
960 336 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
961 : call dotprod_g(dotr1,doti1,istwf_k,npw_nk*my_nspinor,2,cg(:,icg+1+ibd:ibd+icg+inc),wf1,&
962 336 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
963 336 : if(usepaw==1) then
964 : ia =0
965 672 : do itypat=1,ntypat
966 1008 : do iat=1+ia,nattyp(itypat)+ia
967 6384 : do ilmn1=1,pawtab(itypat)%lmn_size
968 29568 : do ilmn2=1,ilmn1
969 24192 : klmn=((ilmn1-1)*ilmn1)/2+ilmn2
970 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)+&
971 24192 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2))
972 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)-&
973 24192 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2))
974 : dotr1=dotr1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)+&
975 24192 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2))
976 : doti1=doti1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)-&
977 29568 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2))
978 : end do
979 24864 : do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
980 18816 : klmn=((ilmn2-1)*ilmn2)/2+ilmn1
981 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)+&
982 18816 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2))
983 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)-&
984 18816 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2))
985 : dotr1=dotr1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)+&
986 18816 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2))
987 : doti1=doti1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)-&
988 24192 : & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2))
989 : end do
990 : end do
991 : end do
992 672 : ia=ia+nattyp(itypat)
993 : end do
994 : end if
995 336 : dotr=sqrt(dotr**2+doti**2)
996 336 : dotr1=sqrt(dotr1**2+doti1**2)
997 336 : write(std_out,*)'DOTR, DOTR1',dotr,dotr1
998 336 : dotr=dotr1/dotr
999 336 : write(std_out,*)'DOTR',dotr
1000 422940 : deltawf1=zero
1001 336 : if(dotr>=0.9d0) then
1002 387704 : deltawf1(:,:)=cg(:,icg+1+ibd:ibd+icg+inc)-wf1(:,:)
1003 308 : if(usepaw==1) then
1004 308 : alpha(1)=one;alpha(2)=zero
1005 308 : beta(1)=-one;beta(2)=zero
1006 308 : ia =0
1007 308 : call pawcprj_zaxpby(alpha,beta,cprj_k(:,iblockbd:iblockbd),cprj_k1(:,iblockbd:iblockbd))
1008 : end if
1009 : istep1=istep
1010 : else
1011 : istep1=1
1012 : end if
1013 422940 : scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind1)=cg(:,icg+1+ibd:ibd+icg+inc)
1014 422940 : scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind2)=deltawf1(:,:)
1015 336 : if(usepaw==1) then
1016 : call pawcprj_put(atindx1,cprj_k,scf_history%cprj(:,:,ind1),dtset%natom,1,ibg,ikpt,1,isppol,&
1017 : & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
1018 336 : & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1019 : call pawcprj_put(atindx1,cprj_k1,scf_history%cprj(:,:,ind2),dtset%natom,1,ibg,ikpt,1,isppol,&
1020 : & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
1021 336 : & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1022 : end if
1023 :
1024 : ! if(istep1>=3) then
1025 : cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+scf_history%alpha*deltawf1(:,:) &
1026 422940 : & +scf_history%beta *deltawf2(:,:)
1027 :
1028 : ! to be used later
1029 : ! if(usepaw==1) then
1030 : ! alpha(2)=zero
1031 : ! beta(1)=one;beta(2)=zero
1032 : ! alpha(1)=scf_history%alpha
1033 : ! call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
1034 : ! alpha(1)=scf_history%beta
1035 : ! call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
1036 : ! call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,&
1037 : ! & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
1038 : ! & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1039 : ! end if
1040 : ! else if (istep1==2) then
1041 : ! cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+scf_history%alpha*deltawf1(:,:)+scf_history%beta*wf1(:,:)
1042 : ! ! cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+deltawf1(:,:)
1043 : ! if(usepaw==1) then
1044 : ! alpha(2)=zero
1045 : ! beta(1)=one;beta(2)=zero
1046 : ! alpha(1)=scf_history%alpha
1047 : ! call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
1048 : ! alpha(1)=scf_history%beta
1049 : ! call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd))
1050 : ! call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,&
1051 : ! & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,&
1052 : ! & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1053 : ! end if
1054 : ! end if
1055 364 : ibd=ibd+inc
1056 : end do ! end loop on iblockbd
1057 :
1058 28 : ABI_FREE(deltawf1)
1059 28 : ABI_FREE(deltawf2)
1060 28 : ABI_FREE(wf1)
1061 28 : ABI_FREE(cwavef)
1062 28 : ABI_FREE(cwavef1)
1063 28 : ABI_FREE(eig)
1064 28 : ABI_FREE(evec)
1065 28 : ABI_FREE(unm)
1066 28 : if(usepaw==1) then
1067 28 : call pawcprj_free(cprj_k)
1068 700 : ABI_FREE(cprj_k)
1069 28 : call pawcprj_free(cprj_k1)
1070 700 : ABI_FREE(cprj_k1)
1071 28 : call pawcprj_free(cprj_k2)
1072 700 : ABI_FREE(cprj_k2)
1073 : end if
1074 :
1075 28 : ibg=ibg+my_nspinor*nband_k
1076 63 : icg=icg+my_nspinor*nband_k*npw_k
1077 :
1078 : ! End big k point loop
1079 : end do
1080 : ! End loop over spins
1081 : end do
1082 :
1083 7 : if(usepaw==1) then
1084 7 : call pawcprj_free(cprj)
1085 679 : ABI_FREE(cprj)
1086 : end if
1087 : if (nprocband>1) then
1088 : ABI_FREE(npw_block)
1089 : ABI_FREE(npw_disp)
1090 : ABI_FREE(bufsize)
1091 : ABI_FREE(bufdisp)
1092 : ABI_FREE(bufsize_wf)
1093 : ABI_FREE(bufdisp_wf)
1094 : end if
1095 :
1096 : end if ! istep>=2
1097 :
1098 8 : if (usepaw==1) then
1099 8 : ABI_FREE(dimcprj)
1100 : end if
1101 :
1102 8 : end subroutine extrapwf
1103 : !!***
1104 :
1105 :
1106 : !!****f* ABINIT/extrapwf_biortho
1107 : !!
1108 : !! NAME
1109 : !! extrapwf_biortho
1110 : !!
1111 : !! FUNCTION
1112 : !! Extrapolate wavefunctions for new ionic positions
1113 : !! from values of wavefunctions of previous SCF cycle.
1114 : !! Use biorthogonal algorithm proposed XG
1115 : !!
1116 : !! INPUTS
1117 : !! atindx1(dtset%natom)=index table for atoms, inverse of atindx
1118 : !! dtset <type(dataset_type)>=all input variables in this dataset
1119 : !! istep=number of call the routine
1120 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
1121 : !! mcprj=size of cprj array
1122 : !! mpi_enreg=information about MPI parallelization
1123 : !! nattyp(dtset%ntypat)=number of atoms of each type in cell.
1124 : !! npwarr(nkpt)=number of planewaves in basis at this k point
1125 : !! pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
1126 : !!
1127 : !! SIDE EFFECTS
1128 : !! cg(2,mcg)= plane wave wavefunction coefficient
1129 : !! Value from previous SCF cycle is input and stored in some form
1130 : !! Extrapolated value is output
1131 : !! cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors
1132 : !! Value from previous SCF cycle is input and stored in some form
1133 : !! Extrapolated value is output
1134 : !! scf_history_wf <type(scf_history_type)>=arrays obtained from previous SCF cycles
1135 : !!
1136 : !! SOURCE
1137 :
1138 0 : subroutine extrapwf_biortho(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,&
1139 0 : & nattyp,npwarr,pawtab,scf_history_wf)
1140 :
1141 : !use m_scf_history
1142 : use m_cgcprj, only : dotprod_set_cgcprj,cgcprj_cholesky,lincom_cgcprj
1143 :
1144 : !Arguments ------------------------------------
1145 : !scalars
1146 : integer,intent(in) :: istep,mcg,mcprj
1147 : type(MPI_type),intent(in) :: mpi_enreg
1148 : type(dataset_type),intent(in) :: dtset
1149 : type(scf_history_type),intent(inout) :: scf_history_wf
1150 : !arrays
1151 : integer,intent(in) :: atindx1(dtset%natom),nattyp(dtset%ntypat)
1152 : integer,intent(in) :: npwarr(dtset%nkpt)
1153 : real(dp), intent(inout) :: cg(2,mcg)
1154 : type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj)
1155 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
1156 :
1157 : !Local variables-------------------------------
1158 : !scalars
1159 : integer :: hermitian
1160 : integer :: ibdmix,ibg,ibg_hist,icg,icg_hist !,iband
1161 : integer :: ierr,ikpt,indh,ind1,ind2,ind1new,inplace
1162 : integer :: isppol,istwf_k,kk,me_distrb,mband,my_nspinor,mcprj_k
1163 : integer :: nband_k,nbdmix,nbdmax,npw_k,ntypat
1164 : integer :: spaceComm_band,usepaw
1165 : real(dp) :: alpha,beta !,dotr,doti
1166 :
1167 : !arrays
1168 0 : integer,allocatable :: ipiv(:),dimcprj(:)
1169 0 : real(dp),allocatable ::psi_ortho(:,:),mmn(:,:,:)
1170 0 : real(dp),allocatable :: smn(:,:,:)
1171 0 : type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kh(:,:)
1172 :
1173 : ! *************************************************************************
1174 :
1175 0 : if (istep==0) return
1176 :
1177 0 : ntypat=dtset%ntypat
1178 0 : usepaw=dtset%usepaw
1179 0 : mband=dtset%mband
1180 0 : nbdmax=dtset%mband
1181 0 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
1182 0 : me_distrb=mpi_enreg%me_kpt
1183 0 : spaceComm_band=xmpi_comm_self
1184 :
1185 : !scf_history_wf%alpha contains dtset%wfmix
1186 0 : alpha=scf_history_wf%alpha
1187 0 : beta=scf_history_wf%beta
1188 0 : ind1=scf_history_wf%hindex(1)
1189 0 : ind2=scf_history_wf%hindex(2)
1190 0 : ind1new=scf_history_wf%hindex(3)
1191 0 : icg=0
1192 0 : icg_hist=0
1193 0 : ibg=0
1194 0 : ibg_hist=0
1195 :
1196 : !Useful array
1197 0 : ABI_MALLOC(dimcprj,(dtset%natom))
1198 0 : if (usepaw==1) then
1199 0 : call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O')
1200 : end if
1201 :
1202 0 : if(istep==1)then
1203 0 : do indh=1,scf_history_wf%history_size
1204 0 : call pawcprj_alloc(scf_history_wf%cprj(:,:,indh),0,dimcprj)
1205 : end do
1206 : end if
1207 :
1208 0 : mcprj_k=my_nspinor*nbdmax
1209 0 : ABI_MALLOC(cprj_k,(dtset%natom,mcprj_k))
1210 0 : ABI_MALLOC(cprj_kh,(dtset%natom,mcprj_k))
1211 :
1212 0 : if(usepaw==1) then
1213 0 : call pawcprj_alloc(cprj_k,0,dimcprj)
1214 0 : call pawcprj_alloc(cprj_kh,0,dimcprj)
1215 : end if
1216 0 : ABI_MALLOC(smn,(2,nbdmax,nbdmax))
1217 0 : ABI_MALLOC(mmn,(2,nbdmax,nbdmax))
1218 :
1219 : !Explanation for the index for the wavefunction stored in scf_history_wf
1220 : !The reference is the cg+cprj output after the wf optimization at istep 1.
1221 : !For wavefunction mixing for molecular dynamics, we use the same mixing as for the density in extraprho. To keep the same indexes,
1222 : ! we choose to take indh=3 for the reference.
1223 :
1224 : !First step
1225 0 : if (istep==1) then
1226 :
1227 0 : indh=3 ! This input wavefunction is the reference
1228 :
1229 : ! LOOP OVER SPINS
1230 0 : do isppol=1,dtset%nsppol
1231 :
1232 : ! BIG FAT k POINT LOOP
1233 0 : do ikpt=1,dtset%nkpt
1234 :
1235 : ! Select k point to be treated by this proc
1236 0 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1237 0 : nbdmix=min(nband_k,nbdmax)
1238 :
1239 0 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
1240 :
1241 0 : npw_k=npwarr(ikpt)
1242 :
1243 0 : scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
1244 :
1245 0 : if(usepaw==1) then
1246 : ! scf_history_wf%cprj(:,ibg_hist+1:ibg_hist+my_nspinor*nbdmix,1)=cprj(:,ibg+1:ibg+my_nspinor*nbdmix)
1247 : call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,0,isppol,mband,&
1248 : & dtset%mkmem,dtset%natom,nbdmax,nbdmix,my_nspinor,dtset%nsppol,0,&
1249 0 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1250 : call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,0,isppol,&
1251 : & mband,dtset%mkmem,dtset%natom,nbdmax,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
1252 0 : & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1253 : end if
1254 :
1255 : ! Update the counters
1256 0 : ibg=ibg+my_nspinor*nband_k
1257 0 : ibg_hist=ibg_hist+my_nspinor*nbdmix
1258 0 : icg=icg+my_nspinor*nband_k*npw_k
1259 0 : icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
1260 :
1261 : end do
1262 : end do
1263 :
1264 : else
1265 : ! From istep==2
1266 0 : if (istep==2) ind1=3
1267 0 : if (istep==3) ind2=3
1268 : ! biorthogonalization
1269 0 : indh=3 ! This input wavefunction is the reference
1270 :
1271 : ! LOOP OVER SPINS
1272 0 : do isppol=1,dtset%nsppol
1273 :
1274 : ! BIG FAT k POINT LOOP
1275 0 : do ikpt=1,dtset%nkpt
1276 :
1277 : ! Select k point to be treated by this proc
1278 0 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1279 0 : nbdmix=min(nband_k,nbdmax)
1280 0 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
1281 0 : istwf_k=dtset%istwfk(ikpt)
1282 0 : npw_k=npwarr(ikpt)
1283 0 : ABI_MALLOC(psi_ortho,(2,npw_k*my_nspinor*nbdmix))
1284 0 : psi_ortho=zero
1285 : ! Biorthogonalization
1286 :
1287 0 : if(usepaw==1) then
1288 : call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,0,isppol,mband,&
1289 : & dtset%mkmem,dtset%natom,nbdmax,nbdmix,my_nspinor,dtset%nsppol,0,&
1290 0 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1291 : call pawcprj_get(atindx1,cprj_kh,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,0,isppol,&
1292 : & mband,dtset%mkmem,dtset%natom,nbdmax,nbdmix,my_nspinor,dtset%nsppol,0,&
1293 0 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1294 : end if !end usepaw=1
1295 :
1296 0 : hermitian=0
1297 :
1298 : call dotprod_set_cgcprj(atindx1,scf_history_wf%cg(:,:,indh),cg,cprj_kh,cprj_k,dimcprj,hermitian,&
1299 : & 0,0,icg_hist,icg,ikpt,isppol,istwf_k,nbdmax,mcg,mcg,mcprj_k,mcprj_k,dtset%mkmem,&
1300 0 : & mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn(:,1:nbdmix,1:nbdmix),usepaw)
1301 :
1302 : ! Invert S matrix, that is NOT hermitian.
1303 : ! Calculate M=S^-1
1304 0 : mmn=zero
1305 0 : do kk=1,nbdmix
1306 0 : mmn(1,kk,kk)=one
1307 : end do
1308 :
1309 0 : ABI_MALLOC(ipiv,(nbdmix))
1310 : ! The smn is destroyed by the following inverse call
1311 0 : call zgesv(nbdmix,nbdmix,smn,nbdmax,ipiv,mmn,nbdmax,ierr)
1312 0 : ABI_FREE(ipiv)
1313 : !DEBUG
1314 0 : if(ierr/=0)then
1315 0 : ABI_ERROR(' The call to cgesv general inversion routine failed')
1316 : end if
1317 : !ENDDEBUG
1318 :
1319 : ! The M matrix is used to compute the biorthogonalized set of wavefunctions, and to store it at the proper place
1320 0 : inplace=0
1321 : call lincom_cgcprj(mmn(:,1:nbdmix,1:nbdmix),cg,cprj_k,dimcprj,&
1322 : & icg,inplace,mcg,mcprj_k,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw,&
1323 0 : & cgout=psi_ortho,cprjout=cprj_kh,icgout=0)
1324 :
1325 : !!!TEST
1326 : ! if (usepaw==0) then
1327 : ! do iband=1,nband_k
1328 : ! 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),&
1329 : !& psi_ortho(:,1+(iband-1)*my_nspinor*npw_k:iband*my_nspinor*npw_k),mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1330 : ! write(80+mpi_enreg%me,*) dotr,doti
1331 : ! flush(80+mpi_enreg%me)
1332 : ! end do
1333 : ! else
1334 : ! hermitian=0
1335 : ! call dotprod_set_cgcprj(atindx1,scf_history_wf%cg(:,:,indh),psi_ortho,scf_history_wf%cprj(:,:,indh)&
1336 : !& ,cprj_kh,dimcprj,hermitian,0,0,icg_hist,icg_hist,ikpt,isppol,istwf_k,nbdmax,mcg,mcg,mcprj_k,mcprj_k,dtset%mkmem,&
1337 : !& mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn(:,1:nbdmix,1:nbdmix),usepaw)
1338 : ! write(90+mpi_enreg%me,*) smn
1339 : ! flush(90+mpi_enreg%me)
1340 : ! end if
1341 : !!!TEST
1342 : ! The biorthogonalised set of wavefunctions is now stored at the proper place
1343 :
1344 : ! cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=zero
1345 :
1346 : ! psi(t+dt) <- psi(t) + alpha.psi(t)
1347 0 : cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=(one+alpha)*psi_ortho(:,1:my_nspinor*npw_k*nbdmix)
1348 0 : if(usepaw==1) then
1349 0 : do ibdmix=1,nbdmix
1350 0 : call pawcprj_axpby(one+alpha,zero,cprj_kh(:,ibdmix:ibdmix),cprj_k(:,ibdmix:ibdmix))
1351 : end do ! end loop on ibdmix
1352 : end if
1353 : ! psi(t+dt) <- -alpha.psi(t-dt) + beta.psi(t-dt)
1354 0 : if (abs(beta-alpha)>tol14.and.ind1>0) then
1355 : cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)+&
1356 0 : & (beta-alpha)*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind1)
1357 0 : if(usepaw==1) then
1358 0 : do ibdmix=1,nbdmix
1359 : call pawcprj_axpby(beta-alpha,one,scf_history_wf%cprj(:,ibg_hist+ibdmix:ibg_hist+ibdmix,ind1),&
1360 0 : & cprj_k(:,ibdmix:ibdmix))
1361 : end do ! end loop on ibdmix
1362 : end if
1363 : end if
1364 :
1365 : ! psi(t+dt) <- -beta.psi(t-2dt)
1366 0 : if (abs(beta)>tol14.and.ind2>0) then
1367 : cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)&
1368 0 : & -beta*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind2)
1369 0 : if(usepaw==1) then
1370 0 : do ibdmix=1,nbdmix
1371 0 : call pawcprj_axpby(-beta,one,scf_history_wf%cprj(:,ibg_hist+ibdmix:ibg_hist+ibdmix,ind2),cprj_k(:,ibdmix:ibdmix))
1372 : end do ! end loop on ibdmix
1373 : end if
1374 : end if
1375 :
1376 : ! Store psi(t) in history
1377 0 : scf_history_wf%cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix,ind1new)=psi_ortho(:,1:my_nspinor*npw_k*nbdmix)
1378 0 : if(usepaw==1) then
1379 : call pawcprj_put(atindx1,cprj_kh,scf_history_wf%cprj(:,:,ind1new),dtset%natom,1,ibg_hist,ikpt,0,isppol,&
1380 : & nbdmix,dtset%mkmem,dtset%natom,nbdmax,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
1381 0 : & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1382 : end if
1383 :
1384 : ! Back to usual orthonormalization for the cg and cprj_k
1385 : call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,mcprj_k,dtset%mkmem,&
1386 0 : & mpi_enreg,dtset%natom,nattyp,nbdmax,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw)
1387 :
1388 : ! Need to transfer cprj_k to cprj
1389 0 : if(usepaw==1) then
1390 : call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,0,isppol,&
1391 : & mband,dtset%mkmem,dtset%natom,nbdmax,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
1392 0 : & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
1393 : end if
1394 :
1395 0 : ibg=ibg+my_nspinor*nband_k
1396 0 : ibg_hist=ibg_hist+my_nspinor*nbdmix
1397 0 : icg=icg+my_nspinor*nband_k*npw_k
1398 0 : icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
1399 0 : ABI_FREE(psi_ortho)
1400 : ! End big k point loop
1401 : end do
1402 : ! End loop over spins
1403 : end do
1404 :
1405 :
1406 : end if !istep>1
1407 :
1408 :
1409 0 : if(usepaw==1) then
1410 0 : call pawcprj_free(cprj_k)
1411 0 : call pawcprj_free(cprj_kh)
1412 : end if
1413 0 : ABI_FREE(cprj_k)
1414 0 : ABI_FREE(cprj_kh)
1415 0 : ABI_FREE(dimcprj)
1416 0 : ABI_FREE(mmn)
1417 0 : ABI_FREE(smn)
1418 :
1419 :
1420 :
1421 0 : end subroutine extrapwf_biortho
1422 : !!***
1423 : end module m_extraprho
1424 : !!***
|
