Branch data 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 : : !!***
|
