Line data Source code
1 : !!****m* ABINIT/m_cgwf
2 : !! NAME
3 : !! m_cgwf
4 : !!
5 : !! FUNCTION
6 : !! Conjugate-gradient eigensolver.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (DCA, XG, GMR, MT, MVeithen, ISouza, JIniguez, MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_cgwf
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_xmpi
27 : use m_abicore
28 : use m_cgtools
29 : use m_efield
30 : use m_dtfil
31 : use m_distribfft
32 :
33 : use defs_abitypes, only : MPI_type
34 : use defs_datatypes, only : pseudopotential_type
35 : use m_fstrings, only : sjoin, itoa, ftoa, ktoa
36 : use m_dtset, only : dataset_type
37 : use m_hdr, only : hdr_type, fform_contains
38 : use m_time, only : timab
39 : use m_numeric_tools, only : rhophi
40 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_put, pawcprj_copy, &
41 : pawcprj_get, pawcprj_mpi_allgather, pawcprj_free, pawcprj_symkn
42 : use m_hamiltonian, only : gs_hamiltonian_type, gspot_transgrid_and_pack
43 : use m_crystal, only : crystal_t
44 : use m_fock, only : fock_set_ieigen, fock_set_getghc_call
45 : use m_getghc, only : getghc
46 : use m_berrytk, only : smatrix
47 : use m_nonlop, only : nonlop
48 : use m_paw_overlap, only : smatrix_k_paw
49 : use m_cgprj, only : getcprj
50 : use m_mpinfo, only : initmpi_seq, destroy_mpi_enreg
51 : use m_pawtab, only : pawtab_type
52 : use m_fftcore, only : get_kg, ngfft_seq
53 : use m_kg, only : mkkin, mkkpg, getph
54 : use m_ioarr, only : read_rhor
55 : use m_pawrhoij, only : pawrhoij_type
56 : use m_pawfgr, only : pawfgr_type
57 : use m_mkffnl, only : mkffnl_objs
58 : use m_initylmg, only : initylmg_k
59 : use m_abi_linalg, only : abi_linalg_init, abi_linalg_finalize
60 : use m_cgtk, only : cgtk_fixphase
61 : use m_gemm_nonlop_projectors, only : gemm_nonlop_use_gemm
62 :
63 : implicit none
64 :
65 : private
66 : !!***
67 :
68 : !!****t* m_cgwf/nscf_t
69 : !! NAME
70 : !! nscf_t
71 : !!
72 : !! FUNCTION
73 : !! Simplified interface to the cgwf routine (congjugate gradient) to perform
74 : !! NSCF calculations starting from the KS potential read from file.
75 : !!
76 : !! SOURCE
77 :
78 : type,public :: nscf_t
79 :
80 : integer :: paral_kgb0 = 0
81 :
82 : integer :: ngfft(18), ngfftf(18)
83 : ! FFT meshes (coarse and fine)
84 :
85 : type(mpi_type) :: mpi_enreg
86 : ! Info on parallelism.
87 :
88 : real(dp),allocatable :: vtrial(:,:)
89 : ! (nfftf, nspden)
90 : ! KS potential in real-space.
91 :
92 : contains
93 :
94 : procedure :: init => nscf_init
95 : ! Initialize the object.
96 :
97 : procedure :: setup_spin => nscf_setup_spin
98 :
99 : procedure :: setup_kpt => nscf_setup_kpt
100 : ! Prepare call to nscf_solve_kpt.
101 : ! Computes k-dependent terms, gs_ham_k and allocates wavefunction block for this k-point
102 :
103 : procedure :: solve_kpt => nscf_solve_kpt
104 : ! Solves the NSCF equation for given k-point and spin.
105 :
106 : procedure :: free => nscf_free
107 : ! Free dynamic memory.
108 :
109 : end type nscf_t
110 : !!***
111 :
112 : public :: cgwf ! Conjugate-gradient solver.
113 : !!***
114 :
115 : contains
116 : !!***
117 :
118 : !!****f* m_cgwf/cgwf
119 : !! NAME
120 : !! cgwf
121 : !!
122 : !! FUNCTION
123 : !! Update all wavefunction |C>, non self-consistently.
124 : !! also compute the corresponding H|C> and Vnl|C> (and S|C> if paw).
125 : !! Uses a conjugate-gradient algorithm.
126 : !! In case of PAW, resolves a generalized eigenproblem using an
127 : !! overlap matrix (not used for norm conserving psps).
128 : !!
129 : !! INPUTS
130 : !! berryopt == 4/14: electric field is on;
131 : ! 6/7/16/7: electric displacement field is on;
132 : !! all other values, no field is present
133 : !! chkexit= if non-zero, check whether the user wishes to exit
134 : !! cpus = CPU time limit
135 : !! filnam_ds1=name of input file (used for exit checking)
136 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
137 : !! icg=shift to be applied on the location of data in the array cg
138 : !! igsc=shift to be applied on the location of data in the array gsc
139 : !! ikpt=number of the k-point
140 : !! inonsc=index of non self-consistent loop
141 : !! isppol=spin polarization currently treated
142 : !! mband =maximum number of bands
143 : !! mcg=second dimension of the cg array
144 : !! mcgq=second dimension of the cgq array
145 : !! mgsc=second dimension of the gsc array
146 : !! mkgq = second dimension of pwnsfacq
147 : !! mpi_enreg=information about MPI parallelization
148 : !! mpw=maximum dimensioned size of npw
149 : !! nband=number of bands.
150 : !! nbdblock=number of bands in a block
151 : !! nkpt=number of k points
152 : !! nline=number of line minimizations per band.
153 : !! npw=number of planewaves in basis sphere at given k.
154 : !! npwarr(nkpt)=number of planewaves in basis at this k point
155 : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
156 : !! nsppol=number of spin polarizations
157 : !! ortalg=governs the choice of the algorithm for orthogonalisation.
158 : !! prtvol=control print volume and debugging output
159 : !! pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat between k-points (see initberry.f)
160 : !! pwind_alloc = first dimension of pwind
161 : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations (see initberry.f)
162 : !! pwnsfacq(2,mkgq) = phase factors for the nearest neighbours of the current k-point (electric field, MPI //)
163 : !! quit= if 1, proceeds to smooth ending of the job.
164 : !! tolrde=tolerance on the ratio of differences of energies (for the line minimisation)
165 : !! tolwfr=tolerance on largest wf residual
166 : !! use_subovl=1 if the overlap matrix is not identity in WFs subspace
167 : !! use_subvnlx=1 if subvnlx has to be computed
168 : !! wfoptalg=govern the choice of algorithm for wf optimisation
169 : !! (0, 1, 10 and 11 : in the present routine, usual CG algorithm;
170 : !! (2 and 3: use shifted square Hamiltonian)
171 : !! zshift(nband)=in case wfoptalg is 2 or 3, shift of the Hamiltonian
172 : !!
173 : !! OUTPUT
174 : !! dphase_k(3) = change in Zak phase for the current k-point in case berryopt = 4/14,6/16,7/17 (electric (displacement) field)
175 : !! resid(nband)=wf residual for new states=|(H-e)|C>|^2 (hartree^2)
176 : !! subham(nband*(nband+1))=Hamiltonian expressed in the WFs subspace
177 : !! subovl(nband*(nband+1)*use_subovl)=overlap matrix expressed in the WFs subspace
178 : !! subvnlx(nband*(nband+1)*use_subvnlx))=non-local Hamiltonian (if NCPP) plus Fock ACE operator (if usefock_ACE) expressed in the WFs subspace
179 : !!
180 : !! SIDE EFFECTS
181 : !! cg(2,mcg)
182 : !! at input =wavefunction <G|C band,k> coefficients for ALL bands
183 : !! at output same as input except that the current band, with number 'band' has been updated
184 : !! dtefield <type(efield_type)> = variables related to Berry phase calculations (see initberry.f)
185 : !!
186 : !! if(gs_hamk%usepaw==1)
187 : !! gsc(2,mgsc)=<G|S|C band,k> coefficients for ALL bands where S is the overlap matrix (used only for paw)
188 : !!
189 : !! NOTES
190 : !! 1) cg should not be filtered and normalized: it should already be OK at input !
191 : !! 2) Not sure that that the generalized eigenproblem (when gs_hamk%usepaw=1)
192 : !! is compatible with wfoptalg=2 or 3 (use of shifted square Hamiltonian) - to be verified
193 : !!
194 : !! SOURCE
195 :
196 442732 : subroutine cgwf(berryopt,cg,cgq,chkexit,cpus,dphase_k,dtefield,&
197 442732 : filnam_ds1,gsc,gs_hamk,icg,igsc,ikpt,inonsc,&
198 : isppol,mband,mcg,mcgq,mgsc,mkgq,mpi_enreg,&
199 442732 : mpw,nband,nbdblock,nkpt,nline,npw,npwarr,&
200 442732 : nspinor,nsppol,ortalg,prtvol,pwind,&
201 442732 : pwind_alloc,pwnsfac,pwnsfacq,quit,resid,subham,subovl,&
202 442732 : subvnlx,tolrde,tolwfr,use_subovl,use_subvnlx,wfoptalg,zshift)
203 :
204 : !Arguments ------------------------------------
205 : integer,intent(in) :: berryopt,chkexit,icg,igsc,ikpt,inonsc,isppol
206 : integer,intent(in) :: mband,mcg,mcgq,mgsc,mkgq,mpw,nband,nbdblock,nkpt,nline
207 : integer,intent(in) :: npw,nspinor,nsppol,ortalg,prtvol,pwind_alloc
208 : integer,intent(in) :: use_subovl,use_subvnlx,wfoptalg, quit
209 : real(dp),intent(in) :: cpus,tolrde,tolwfr
210 : character(len=*),intent(in) :: filnam_ds1
211 : type(MPI_type),intent(in) :: mpi_enreg
212 : type(efield_type),intent(inout) :: dtefield
213 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
214 : !arrays
215 : integer,intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3)
216 : real(dp),intent(in) :: cgq(2,mcgq)
217 : real(dp),intent(in) :: pwnsfac(2,pwind_alloc),pwnsfacq(2,mkgq),zshift(nband)
218 : real(dp),intent(inout) :: cg(2,mcg),gsc(2,mgsc)
219 : real(dp),intent(inout) :: dphase_k(3)
220 : real(dp),intent(out) :: subham(nband*(nband+1)),subovl(nband*(nband+1)*use_subovl)
221 : real(dp),intent(out) :: subvnlx(nband*(nband+1)*use_subvnlx)
222 : real(dp),intent(out) :: resid(nband)
223 :
224 : !Local variables-------------------------------
225 : integer,parameter :: level=113,tim_getghc=1,tim_projbd=1,type_calc=0, enough_warning = 3
226 : integer,save :: nskip=0
227 : integer :: choice,counter,cpopt,ddkflag,dimenlc1,dimenlr1,dimenl2,iat,iatom,itypat
228 : integer :: iband,ibandmin,ibandmax,me_g0
229 : integer :: ibdblock,iblock,icg1,icg_shift,icp1,icp2,idir,idum1,ierr,ifor,igs,igsc_shift,ii,ikgf
230 : integer :: ikpt2,ikpt2f,ikptf,iline,iproc,ipw,ispinor,istwf_k,isubh,isubo,itrs
231 : integer :: job,mcg_q,me_distrb,natom,ncpgr,nblock,nproc_distrb,npw_k2
232 : integer :: optekin,paw_opt,signs,shiftbd,sij_opt,spaceComm_distrb
233 : integer :: useoverlap,wfopta10,gpu_option_save, num_warning
234 : real(dp) :: chc,costh,deltae,deold,dhc,dhd,diff,dotgg,dotgp,doti,dotr
235 : real(dp) :: dphase_aux2,e0,e0_old,e1,e1_old,eval,gamma
236 : real(dp) :: lam0,lamold,root,sinth,sintn,swap,tan2th,theta,thetam, xnorm
237 : logical :: gen_eigenpb, finite_field
238 : character(len=500) :: msg
239 : integer :: hel(2,3)
240 442732 : integer,allocatable :: dimlmn(:),dimlmn_srt(:),ikptf_recv(:),pwind_k(:),sflag_k(:)
241 : real(dp) :: bcut(2,3),dphase_aux1(3),dtm_k(2),phase_end(3), phase_init(3),tsec(2), gsc_dummy(0,0)
242 885464 : real(dp),allocatable :: cg1_k(:,:),cgq_k(:,:),conjgr(:,:),cwavef(:,:)
243 442732 : real(dp),allocatable :: detovc(:,:,:),detovd(:,:,:),direc(:,:),direc_tmp(:,:)
244 442732 : real(dp),allocatable :: gh_direc(:,:),gh_direcws(:,:),ghc(:,:),ghc_all(:,:),ghcws(:,:)
245 442732 : real(dp),allocatable :: grad_berry(:,:),grad_total(:,:),gs_direc(:,:)
246 442732 : real(dp),allocatable :: gvnlxc(:,:),gvnlx_direc(:,:),gvnlx_dummy(:,:)
247 442732 : real(dp),allocatable :: pcon(:),pwnsfac_k(:,:),scprod(:,:),scwavef(:,:)
248 442732 : real(dp),allocatable :: smat_inv(:,:,:),smat_k(:,:,:),smat_k_paw(:,:,:),swork(:,:),vresid(:,:),work(:,:)
249 3099124 : type(pawcprj_type) :: cprj_dum(1,1)
250 442732 : type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kb(:,:)
251 442732 : type(pawcprj_type),allocatable :: cprj_direc(:,:),cprj_band_srt(:,:),cprj_gat(:,:), cprj_fkn(:,:),cprj_ikn(:,:)
252 : ! *********************************************************************
253 :
254 : DBG_ENTER("COLL")
255 :
256 : !Starting the routine
257 442732 : call timab(22,1,tsec)
258 :
259 : !Touching chkexit, cpus,filnam_ds to avoid warning for abirules. This is dirty...
260 442732 : if(chkexit<0)then
261 0 : ABI_BUG('chkexit should be positive!')
262 : end if
263 :
264 442732 : if(cpus<0 .and. filnam_ds1=='a')then
265 0 : ABI_BUG('cpus should be positive!')
266 : end if
267 :
268 : !======================================================================
269 : !========= LOCAL VARIABLES DEFINITIONS AND ALLOCATIONS ================
270 : !======================================================================
271 :
272 : !MPI data
273 442732 : spaceComm_distrb=mpi_enreg%comm_cell
274 442732 : nproc_distrb=xmpi_comm_size(spaceComm_distrb)
275 442732 : me_distrb=mpi_enreg%me_kpt
276 442732 : me_g0 = mpi_enreg%me_g0
277 442732 : gpu_option_save=gs_hamk%gpu_option
278 : !Using cgwf with ndat=1 in getghc calls isn't beneficial in most cases, hence we disable GPU here
279 442732 : gs_hamk%gpu_option=ABI_GPU_DISABLED
280 442732 : if(gpu_option_save==ABI_GPU_OPENMP) gemm_nonlop_use_gemm=.false.
281 :
282 : !if PAW, one has to solve a generalized eigenproblem (H|Psi>=Lambda.S|Psi>)
283 : !else, one has to solve a classical eigenproblem (H|Psi>=Lambda.|Psi>)
284 442732 : gen_eigenpb=(gs_hamk%usepaw==1)
285 442732 : useoverlap=0;if (gen_eigenpb) useoverlap=1
286 :
287 : !Initializations and allocations
288 442732 : isubh=1;isubo=1
289 442732 : nblock=(nband-1)/nbdblock+1
290 442732 : istwf_k=gs_hamk%istwf_k
291 442732 : wfopta10=mod(wfoptalg,10)
292 442732 : optekin=0;if (wfoptalg>=10) optekin=1
293 442732 : natom=gs_hamk%natom
294 442732 : cpopt=-1
295 442732 : num_warning = 0
296 :
297 1328196 : ABI_MALLOC(pcon,(npw))
298 1328196 : ABI_MALLOC(ghc,(2,npw*nspinor))
299 885464 : ABI_MALLOC(gvnlxc,(2,npw*nspinor))
300 885464 : ABI_MALLOC(conjgr,(2,npw*nspinor))
301 885464 : ABI_MALLOC(cwavef,(2,npw*nspinor))
302 885464 : ABI_MALLOC(direc,(2,npw*nspinor))
303 1328196 : ABI_MALLOC(scprod,(2,nband))
304 :
305 885464 : ABI_MALLOC(gh_direc,(2,npw*nspinor))
306 885464 : ABI_MALLOC(gvnlx_direc,(2,npw*nspinor))
307 885464 : ABI_MALLOC(vresid,(2,npw*nspinor))
308 :
309 442732 : if (gen_eigenpb) then
310 195190 : ABI_MALLOC(gs_direc,(2,npw*nspinor))
311 : else
312 345137 : ABI_MALLOC(gs_direc,(0,0))
313 : end if
314 :
315 442732 : if (gen_eigenpb) then
316 195190 : ABI_MALLOC(scwavef,(2,npw*nspinor))
317 195190 : ABI_MALLOC(direc_tmp,(2,npw*nspinor))
318 : end if
319 :
320 97595 : if (gen_eigenpb.and.(inonsc==1)) then
321 219549 : ABI_MALLOC_OR_DIE(ghc_all,(2,nband*npw*nspinor), ierr)
322 : end if
323 :
324 442732 : if (wfopta10==2.or.wfopta10==3) then
325 118 : ABI_MALLOC(work,(2,npw*nspinor))
326 : end if
327 :
328 442732 : if (gen_eigenpb.and.(wfopta10==2.or.wfopta10==3)) then
329 0 : ABI_MALLOC(swork,(2,npw*nspinor))
330 : else
331 442732 : ABI_MALLOC(swork,(0,0))
332 : end if
333 :
334 442732 : if (wfopta10==2 .or. wfopta10==3) then
335 118 : ABI_MALLOC(ghcws,(2,npw*nspinor))
336 118 : ABI_MALLOC(gh_direcws,(2,npw*nspinor))
337 118 : ABI_MALLOC(gvnlx_dummy,(2,npw*nspinor))
338 : else
339 442673 : ABI_MALLOC(gvnlx_dummy,(0,0))
340 : end if
341 :
342 : !if "generalized eigenproblem", not sure of wfoptalg=2,3 algorithms
343 442732 : if ((gen_eigenpb).and.(wfopta10==2.or.wfopta10==3)) then
344 : write(msg, '(a,a,a,a,a)' ) &
345 0 : & 'Conjugate gradient algorithm not tested with',ch10,&
346 0 : & 'wfoptalg=2 or 3 and usepaw==1 !',ch10,&
347 0 : & 'Program will continue at your own risk...'
348 0 : ABI_WARNING(msg)
349 : end if
350 :
351 : !Electric field: definition of local variables:
352 : !detovc(1:2,ifor,idir) determinant of the overlap matrix
353 : !S_{nm}(k,k+dk)=<u_{n,k}|u_{m,k+dk}>, with the states at
354 : !k as bras (neighbor is specified by ifor and idir)
355 : !detovd same as detovc but with <u_{n,k}| replaced by
356 : !<D| (search direction) in the band-th line
357 : !grad_berry(1:2,ipw) Berry phase term contribution to the gradient vector
358 : !hel(ifor,idir) helicity of the ellipse associated w/ (ifor,idir)
359 : !bcut(ifor,idir) branch cut of the ellipse associated w/ (ifor,idir)
360 : !theta_min optimal angle theta in line_minimization when electric
361 : !field is on
362 : !grad_total(1:2,ipw) total gradient (zero field term plus Berry phase term)
363 :
364 : finite_field = ( (berryopt == 4) .or. (berryopt == 6) .or. (berryopt == 7) .or. &
365 442732 : & (berryopt == 14) .or. (berryopt == 16) .or. (berryopt == 17) )
366 442732 : ncpgr = 0 ! do not think the cprj's here need gradients (no force computation in cgwf)
367 :
368 : if (finite_field) then
369 :
370 : ! ji : These could be a couple of new input variables (but it is OK to define them here)
371 8330 : ikptf = dtefield%i2fbz(ikpt)
372 8330 : ikgf = dtefield%fkgindex(ikptf) ! this is the shift for pwind
373 8330 : mcg_q = mpw*mband*nspinor
374 8330 : ABI_MALLOC(detovc,(2,2,3))
375 8330 : ABI_MALLOC(detovd,(2,2,3))
376 16660 : ABI_MALLOC(grad_berry,(2,npw*nspinor))
377 24990 : ABI_MALLOC(cg1_k,(2,mpw))
378 24990 : ABI_MALLOC(cgq_k,(2,mcg_q))
379 16660 : ABI_MALLOC(grad_total,(2,npw*nspinor))
380 24990 : ABI_MALLOC(sflag_k,(mband))
381 24990 : ABI_MALLOC(pwind_k,(mpw))
382 24990 : ABI_MALLOC(pwnsfac_k,(4,mpw))
383 33320 : ABI_MALLOC(smat_k,(2,mband,mband))
384 24990 : ABI_MALLOC(smat_inv,(2,mband,mband))
385 : ! now the special features if using PAW
386 8330 : if (gs_hamk%usepaw /= 0) then
387 22512 : ABI_MALLOC(smat_k_paw,(2,gs_hamk%usepaw*mband,gs_hamk%usepaw*mband))
388 298284 : smat_k_paw(:,:,:) = zero
389 : ! the following are arguments to nonlop used to apply the on-site dipole to direc vector
390 5628 : choice = 1 ! only apply projectors
391 5628 : paw_opt = 1 ! only apply Dij
392 5628 : signs = 2 ! apply nonlop to vector in k space
393 : ! following two items are the nonlocal potential strength dij due to the on-site dipoles
394 5628 : dimenlc1 = 2*gs_hamk%lmnmax*(gs_hamk%lmnmax+1)/2
395 5628 : dimenlr1 = gs_hamk%lmnmax*(gs_hamk%lmnmax+1)/2
396 5628 : dimenl2 = natom
397 : ! cprj structures for finite_field case
398 16884 : ABI_MALLOC(dimlmn,(natom))
399 16884 : do iatom = 1, natom
400 11256 : itypat = gs_hamk%typat(iatom)
401 16884 : dimlmn(iatom)=dtefield%lmn_size(itypat)
402 : end do
403 11256 : ABI_MALLOC(dimlmn_srt,(natom))
404 5628 : iatom = 0
405 16884 : do itypat = 1, gs_hamk%ntypat
406 28140 : do iat = 1, gs_hamk%nattyp(itypat)
407 11256 : iatom = iatom + 1
408 22512 : dimlmn_srt(iatom) = dtefield%lmn_size(itypat)
409 : end do
410 : end do
411 16884 : ABI_MALLOC(ikptf_recv,(nproc_distrb))
412 90048 : ABI_MALLOC(cprj_k,(natom,mband*nspinor))
413 84420 : ABI_MALLOC(cprj_kb,(natom,mband*nspinor))
414 84420 : ABI_MALLOC(cprj_direc,(natom,mband*nspinor))
415 39396 : ABI_MALLOC(cprj_band_srt,(natom,nspinor))
416 178176 : ABI_MALLOC(cprj_gat,(natom,mband*nspinor*nproc_distrb))
417 5628 : call pawcprj_alloc(cprj_k,ncpgr,dimlmn)
418 5628 : call pawcprj_alloc(cprj_kb,ncpgr,dimlmn)
419 5628 : call pawcprj_alloc(cprj_direc,ncpgr,dimlmn)
420 5628 : call pawcprj_alloc(cprj_band_srt,ncpgr,dimlmn_srt)
421 5628 : call pawcprj_alloc(cprj_gat,ncpgr,dimlmn)
422 5628 : if (nkpt /= dtefield%fnkpt) then
423 84420 : ABI_MALLOC(cprj_fkn,(natom,mband*nspinor))
424 84420 : ABI_MALLOC(cprj_ikn,(natom,mband*nspinor))
425 5628 : call pawcprj_alloc(cprj_fkn,ncpgr,dimlmn)
426 5628 : call pawcprj_alloc(cprj_ikn,ncpgr,dimlmn)
427 : end if
428 : else
429 2702 : ABI_MALLOC(smat_k_paw,(0,0,0))
430 2702 : ABI_MALLOC(dimlmn,(0))
431 2702 : ABI_MALLOC(dimlmn_srt,(0))
432 2702 : ABI_MALLOC(ikptf_recv,(0))
433 2702 : ABI_MALLOC(cprj_k,(0,0))
434 2702 : ABI_MALLOC(cprj_kb,(0,0))
435 2702 : ABI_MALLOC(cprj_direc,(0,0))
436 2702 : ABI_MALLOC(cprj_band_srt,(0,0))
437 2702 : ABI_MALLOC(cprj_gat,(0,0))
438 : end if
439 : end if ! finite_field
440 :
441 : ! ======================================================================
442 : ! If generalized eigenproblem: has to know <g|S|c> for all
443 : ! bands (for orthogonalization purpose); take benefit of this
444 : ! calculation to compute <g|H|c> at the same time, and cprj_k if finite_field
445 : ! ======================================================================
446 442732 : if (gen_eigenpb.and.inonsc==1) then
447 958813 : do iblock=1,nblock
448 885630 : ibandmin=1+(iblock-1)*nbdblock
449 885630 : ibandmax=min(iblock*nbdblock,nband)
450 1844443 : do iband=ibandmin,ibandmax
451 885630 : ibdblock=iband-(iblock-1)*nbdblock
452 885630 : icg_shift=npw*nspinor*(iband-1)+icg
453 885630 : igsc_shift=npw*nspinor*(iband-1)+igsc
454 :
455 885630 : call cg_zcopy(npw*nspinor,cg(1,1+icg_shift),cwavef)
456 :
457 : ! Compute <g|H|c>
458 : ! By setting ieigen to iband, Fock contrib. of this iband to the energy will be calculated
459 885630 : call fock_set_ieigen(gs_hamk%fockcommon,iband)
460 885630 : sij_opt=1
461 885630 : if (finite_field .and. gs_hamk%usepaw == 1) then
462 17536 : call getghc(0,cwavef,cprj_band_srt,ghc,scwavef,gs_hamk,gvnlxc,eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,type_calc)
463 : call pawcprj_put(gs_hamk%atindx,cprj_band_srt,cprj_k,natom,iband,0,ikpt,&
464 17536 : 1,isppol,mband,1,natom,1,mband,dimlmn,nspinor,nsppol,0,mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
465 : else
466 868094 : call getghc(cpopt,cwavef,cprj_dum,ghc,scwavef,gs_hamk,gvnlxc,eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,type_calc)
467 : end if
468 :
469 885630 : call cg_zcopy(npw*nspinor,ghc,ghc_all(1,1+icg_shift-icg))
470 1771260 : call cg_zcopy(npw*nspinor,scwavef,gsc(1,1+igsc_shift))
471 : end do
472 : end do
473 : end if
474 :
475 : ! Loop over blocks of bands. In the standard band-sequential algorithm, nblock=nband.
476 3592359 : do iblock=1,nblock
477 3149627 : counter=100*iblock*nbdblock+inonsc
478 :
479 : ! Loop over bands in a block
480 : ! This loop can be MPI-parallelized, over processors attached to the same k point
481 3149627 : ibandmin=1+(iblock-1)*nbdblock
482 3149627 : ibandmax=min(iblock*nbdblock,nband)
483 :
484 : ! Big iband loop
485 6299254 : do iband=ibandmin,ibandmax
486 3149627 : ibdblock=iband-(iblock-1)*nbdblock
487 3149627 : counter=100*iband+inonsc
488 3149627 : icg_shift=npw*nspinor*(iband-1)+icg
489 3149627 : igsc_shift=npw*nspinor*(iband-1)+igsc
490 :
491 : ! ======================================================================
492 : ! ========== INITIALISATION OF MINIMIZATION ITERATIONS =================
493 : ! ======================================================================
494 :
495 3149627 : if (prtvol>=10) then ! Tell us what is going on:
496 52505 : write(msg, '(a,i6,2x,a,i3,a)' )' --- cgwf is called for band',iband,'for',nline,' lines'
497 52505 : call wrtout(std_out,msg,'PERS')
498 : end if
499 :
500 3149627 : dotgp=one
501 3149627 : if (finite_field) then
502 2043360 : detovc(:,:,:) = zero ; detovd(:,:,:) = zero
503 47520 : phase_init(:) = zero
504 47520 : dphase_aux1(:) = zero
505 47520 : phase_end(:) = zero
506 47520 : bcut(:,:) = zero
507 47520 : hel(:,:) = 0
508 : end if
509 :
510 : ! Extraction of the vector that is iteratively updated
511 3149627 : call cg_zcopy(npw*nspinor,cg(1,1+icg_shift),cwavef)
512 :
513 : ! If generalized eigenproblem: extraction of the overlap information
514 3149627 : if (gen_eigenpb) then
515 1243378 : call cg_zcopy(npw*nspinor,gsc(1,1+igsc_shift),scwavef)
516 : end if
517 :
518 : ! Normalize incoming wf (and S.wf, if generalized eigenproblem):
519 : ! WARNING: It might be interesting to skip the following operation.
520 : ! The associated routines should be reexamined to see whether cwavef is not already normalized.
521 : if (gen_eigenpb) then
522 1243378 : call dotprod_g(dotr,doti,istwf_k,npw*nspinor,2,cwavef,scwavef,me_g0,mpi_enreg%comm_spinorfft)
523 1243378 : dotr=sqrt(dotr**2+doti**2); xnorm=one/sqrt(dotr)
524 3730134 : call cg_zscal(npw*nspinor,(/xnorm,zero/),cwavef)
525 3730134 : call cg_zscal(npw*nspinor,(/xnorm,zero/),scwavef)
526 : else
527 1906249 : call sqnorm_g(dotr,istwf_k,npw*nspinor,cwavef,me_g0,mpi_enreg%comm_fft)
528 1906249 : xnorm=one/sqrt(abs(dotr))
529 5718747 : call cg_zscal(npw*nspinor,(/xnorm,zero/),cwavef)
530 : end if
531 :
532 3149627 : if (prtvol==-level) then
533 0 : write(msg,'(a,f14.6)')' cgwf: xnorm = ',xnorm
534 0 : call wrtout(std_out,msg,'PERS')
535 : end if
536 :
537 : ! Compute (or extract) <g|H|c>
538 3149627 : if (gen_eigenpb.and.(inonsc==1)) then
539 :
540 : !$OMP PARALLEL DO PRIVATE(ipw)
541 324256285 : do ipw=1,npw*nspinor
542 323370655 : ghc(1,ipw)=xnorm*ghc_all(1,ipw+icg_shift-icg)
543 324256285 : ghc(2,ipw)=xnorm*ghc_all(2,ipw+icg_shift-icg)
544 : end do
545 :
546 : else
547 : ! By setting ieigen to iband, Fock contrib. of this iband to the energy will be calculated
548 2263997 : call fock_set_ieigen(gs_hamk%fockcommon,iband)
549 2263997 : sij_opt=0
550 2263997 : call getghc(cpopt,cwavef,cprj_dum,ghc,gsc_dummy,gs_hamk,gvnlxc, eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,type_calc)
551 : end if
552 :
553 : ! Minimisation of the residual: compute <G|(H-zshift)^2|C iband,k>
554 3149627 : if(wfopta10==2 .or. wfopta10==3) then
555 92036 : ghcws(:,:)=ghc(:,:)
556 236 : if (gen_eigenpb) then
557 0 : sij_opt=1
558 0 : work(:,:)=ghc(:,:)-zshift(iband)*scwavef(:,:)
559 : else
560 236 : sij_opt=0
561 92036 : work(:,:)=ghc(:,:)-zshift(iband)*cwavef(:,:)
562 : end if
563 236 : call getghc(cpopt,work,cprj_dum,ghc,swork,gs_hamk,gvnlx_dummy, eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,type_calc)
564 236 : if (gen_eigenpb) then
565 0 : ghc(:,:)=ghc(:,:)-zshift(iband)*swork(:,:)
566 : else
567 92036 : ghc(:,:)=ghc(:,:)-zshift(iband)*work(:,:)
568 : end if
569 : end if
570 :
571 : ! ======================================================================
572 : ! ====== BEGIN LOOP FOR A GIVEN BAND: MINIMIZATION ITERATIONS ==========
573 : ! ======================================================================
574 3149627 : if(nline/=0)then
575 12043660 : do iline=1,nline
576 :
577 : ! === COMPUTE THE RESIDUAL ===
578 :
579 : ! Compute lambda = <C|H|C> or <C|(H-zshift)**2|C>
580 10558896 : call dotprod_g(chc,doti,istwf_k,npw*nspinor,1,cwavef,ghc,me_g0,mpi_enreg%comm_spinorfft)
581 10558896 : lam0=chc
582 :
583 : ! Check that lam0 is decreasing on succeeding lines:
584 10558896 : if (.not.finite_field) then
585 10383816 : if (iline==1) then
586 3102107 : lamold=lam0
587 : else
588 7281709 : if (lam0 > lamold+tol12 .and. num_warning <= enough_warning) then
589 0 : num_warning = num_warning + 1
590 : write(msg, '(a,i0,a,e14.6,a,e14.6)')&
591 0 : 'New trial energy at line: ',iline,' = ',lam0,' is higher than former: ',lamold
592 0 : ABI_WARNING(msg)
593 : end if
594 7281709 : lamold=lam0
595 : end if
596 : end if
597 :
598 : ! Compute residual vector:
599 : ! Note that vresid is precomputed to guarantee cancellation of errors
600 : ! and allow residuals to reach values as small as 1.0d-24 or better.
601 :
602 10558896 : if (wfopta10<=1) then
603 10556675 : eval=chc
604 10556675 : if (gen_eigenpb) then
605 :
606 : !$OMP PARALLEL DO
607 1536626771 : do ipw=1,npw*nspinor
608 1532326805 : vresid(1,ipw)=ghc(1,ipw)-chc*scwavef(1,ipw)
609 1536626771 : vresid(2,ipw)=ghc(2,ipw)-chc*scwavef(2,ipw)
610 : end do
611 : else
612 : !$OMP PARALLEL DO
613 2118144519 : do ipw=1,npw*nspinor
614 2111887810 : vresid(1,ipw)=ghc(1,ipw)-chc*cwavef(1,ipw)
615 2118144519 : vresid(2,ipw)=ghc(2,ipw)-chc*cwavef(2,ipw)
616 : end do
617 : end if
618 : else
619 2221 : call dotprod_g(eval,doti,istwf_k,npw*nspinor,1,cwavef,ghcws,me_g0,mpi_enreg%comm_spinorfft)
620 2221 : if (gen_eigenpb) then
621 : !$OMP PARALLEL DO
622 0 : do ipw=1,npw*nspinor
623 0 : vresid(1,ipw)=ghcws(1,ipw)-eval*scwavef(1,ipw)
624 0 : vresid(2,ipw)=ghcws(2,ipw)-eval*scwavef(2,ipw)
625 : end do
626 : else
627 : !$OMP PARALLEL DO
628 288769 : do ipw=1,npw*nspinor
629 286548 : vresid(1,ipw)=ghcws(1,ipw)-eval*cwavef(1,ipw)
630 288769 : vresid(2,ipw)=ghcws(2,ipw)-eval*cwavef(2,ipw)
631 : end do
632 : end if
633 : end if
634 :
635 : ! Compute residual (squared) norm
636 10558896 : call sqnorm_g(resid(iband),istwf_k,npw*nspinor,vresid,me_g0,mpi_enreg%comm_fft)
637 :
638 10558896 : if (prtvol==-level) then
639 0 : write(msg,'(a,i0,2f14.6)')' cgwf: iline,eval,resid = ',iline,eval,resid(iband)
640 0 : call wrtout(std_out,msg,'PERS')
641 : end if
642 :
643 : ! ======================================================================
644 : ! ============== CHECK FOR CONVERGENCE CRITERIA ========================
645 : ! ======================================================================
646 :
647 : ! If residual sufficiently small stop line minimizations
648 10558896 : if (resid(iband)<tolwfr) then
649 322225 : if (prtvol>=10) then
650 : write(msg, '(a,i4,a,i2,a,es12.4)' ) &
651 2191 : ' cgwf: band ',iband,' converged after ',iline,' line minimizations: resid =',resid(iband)
652 2191 : call wrtout(std_out,msg,'PERS')
653 : end if
654 322225 : nskip=nskip+(nline-iline+1) ! Number of two-way 3D ffts skipped
655 322225 : exit ! Exit from the loop on iline
656 : end if
657 :
658 : ! If user require exiting the job, stop line minimisations
659 10236671 : if (quit==1) then
660 0 : write(msg, '(a,i0)' )' cgwf: user require exiting => skip update of band ',iband
661 0 : call wrtout(std_out,msg,'PERS')
662 :
663 0 : nskip=nskip+(nline-iline+1) ! Number of two-way 3D ffts skipped
664 0 : exit ! Exit from the loop on iline
665 : end if
666 :
667 : ! ======================================================================
668 : ! =========== COMPUTE THE STEEPEST DESCENT DIRECTION ===================
669 : ! ======================================================================
670 :
671 : ! Compute the steepest descent direction
672 10236671 : if (gen_eigenpb) then
673 4274655 : call cg_zcopy(npw*nspinor,vresid,direc) ! Store <G|H-lambda.S|C> in direc
674 : else
675 5962016 : call cg_zcopy(npw*nspinor,ghc,direc) ! Store <G|H|C> in direc
676 : end if
677 :
678 : ! Electric field: compute the gradient of the Berry phase part of the energy functional.
679 : ! See PRL 89, 117602 (2002) [[cite:Souza2002]], grad_berry(:,:) is the second term of Eq. (4)
680 10236671 : if (finite_field) then
681 :
682 : call make_grad_berry(cg,cgq,cprj_k,detovc,dimlmn,dimlmn_srt,direc,dtefield,grad_berry,&
683 : gs_hamk,iband,icg,ikpt,isppol,mband,mcg,mcgq,mkgq,mpi_enreg,mpw,natom,nkpt,npw,npwarr,&
684 175080 : nspinor,nsppol,pwind,pwind_alloc,pwnsfac,pwnsfacq)
685 :
686 : ! Add grad_berry to direc and store original gradient
687 72453204 : direc(:,:) = direc(:,:) + grad_berry(:,:)
688 72453204 : grad_total(:,:) = direc(:,:)
689 : ! DEBUG: check that grad_berry is orthogonal to the occupied manifold at k
690 : ! do jband = 1, dtefield%mband_occ
691 : ! dotr = zero ; doti = zero
692 : ! do ipw = 1, npw*nspinor
693 : ! if(.not.gen_eigenpb) then
694 : ! dotr = dotr + cg(1,icg + (jband-1)*npw*nspinor + ipw)*grad_berry(1,ipw) + &
695 : ! & cg(2,icg + (jband-1)*npw*nspinor + ipw)*grad_berry(2,ipw)
696 : ! doti = doti + cg(1,icg + (jband-1)*npw*nspinor + ipw)*grad_berry(2,ipw) - &
697 : ! & cg(2,icg + (jband-1)*npw*nspinor + ipw)*grad_berry(1,ipw)
698 : ! end if
699 : ! end do
700 : ! if ((abs(dotr) > tol12).or.(abs(doti) > tol12)) then
701 : ! write(std_out,'(a)')'cgwf-berry : ERROR (orthogonality)'
702 : ! write(std_out,'(3(2x,i3),2(5x,e16.9))')ikpt,iband,jband,dotr,doti
703 : ! stop
704 : ! end if
705 : ! end do
706 : ! ENDDEBUG
707 : end if ! finite_field
708 :
709 : ! =========== PROJECT THE STEEPEST DESCENT DIRECTION ===================
710 : ! ========= OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ================
711 :
712 : ! The following projection over the subspace orthogonal to occupied bands
713 : ! is optional. It is a bit more accurate, but doubles the number of N^3 ops.
714 : ! It is done only if ortalg>=0.
715 :
716 : ! Project the steepest descent direction:
717 : ! direc(2,npw)=<G|H|Cnk> - \sum_{(i<=n)} <G|H|Cik> , normalized.
718 :
719 10236671 : if(ortalg>=0)then
720 5963096 : if (gen_eigenpb) then
721 : call projbd(cg,direc,iband,icg,igsc,istwf_k,mcg,mgsc,nband,npw,nspinor,&
722 1642 : gsc,scprod,0,tim_projbd,useoverlap,me_g0,mpi_enreg%comm_fft)
723 : else
724 : call projbd(cg,direc,-1 ,icg,igsc,istwf_k,mcg,mgsc,nband,npw,nspinor,&
725 5961454 : gsc,scprod,0,tim_projbd,useoverlap,me_g0,mpi_enreg%comm_fft)
726 : end if
727 : else
728 : ! For negative ortalg must still project current band out of conjugate vector (unneeded if gen_eigenpb)
729 4273575 : if (.not.gen_eigenpb) then
730 562 : call dotprod_g(dotr,doti,istwf_k,npw*nspinor,3,cwavef,direc,me_g0,mpi_enreg%comm_spinorfft)
731 562 : if(istwf_k==1)then
732 1686 : call cg_zaxpy(npw*nspinor,-(/dotr,doti/),cwavef,direc)
733 : else
734 0 : call cg_zaxpy(npw*nspinor,(/-dotr,zero/),cwavef,direc)
735 : end if
736 : end if
737 : end if
738 :
739 : ! For a generalized eigenpb, store the steepest descent direction
740 4595630330 : if (gen_eigenpb) direc_tmp=direc
741 :
742 : ! ======================================================================
743 : ! ======== PRECONDITION THE STEEPEST DESCENT DIRECTION =================
744 : ! ======================================================================
745 :
746 : ! If wfoptalg>=10, the precondition matrix is kept constant during iteration ; otherwise it is recomputed
747 10236671 : if (wfoptalg<10.or.iline==1) then
748 7196099 : if (gs_hamk%use_gbt == 0) then
749 7086996 : call cg_precon(cwavef,zero,istwf_k,gs_hamk%kinpw_k,npw,nspinor,me_g0,optekin,pcon,direc,mpi_enreg%comm_fft)
750 : else
751 109103 : call cg_precon(cwavef,zero,istwf_k,gs_hamk%kinpw_k,npw,1,me_g0,optekin,pcon,direc,mpi_enreg%comm_fft)
752 : call cg_precon(cwavef(:,npw+1:),zero,istwf_k,gs_hamk%kinpw_kp,npw,1,me_g0,optekin,pcon,&
753 109103 : direc(:,npw+1:),mpi_enreg%comm_fft)
754 : end if
755 :
756 7196099 : if(wfopta10==2 .or. wfopta10==3)then
757 : ! Minimisation of the residual: must precondition twice
758 : ! (might make only one call, with modified precon routine - might also make a shift !!!)
759 2216 : call cg_precon(cwavef,zero,istwf_k,gs_hamk%kinpw_k,npw,nspinor,me_g0,optekin,pcon,direc,mpi_enreg%comm_fft)
760 2216 : if(iline==1)then
761 : !$OMP PARALLEL DO
762 30437 : do ipw=1,npw
763 30204 : pcon(ipw)=pcon(ipw)**2
764 30437 : pcon(ipw)=pcon(ipw)**2
765 : end do
766 : end if
767 : end if
768 : else
769 6523988 : do ispinor=1,nspinor
770 3483416 : igs=(ispinor-1)*npw
771 : !$OMP PARALLEL DO
772 1109398972 : do ipw=1+igs,npw+igs
773 1102874984 : direc(1,ipw)=direc(1,ipw)*pcon(ipw-igs)
774 1106358400 : direc(2,ipw)=direc(2,ipw)*pcon(ipw-igs)
775 : end do
776 : end do
777 : end if
778 :
779 : ! ======= PROJECT THE PRECOND. STEEPEST DESCENT DIRECTION ==============
780 : ! ========= OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ================
781 : ! Projecting again out all bands (not normalized).
782 : call projbd(cg,direc,-1,icg,igsc,istwf_k,mcg,mgsc,nband,npw,nspinor,&
783 10236671 : gsc,scprod,0,tim_projbd,useoverlap,me_g0,mpi_enreg%comm_fft)
784 :
785 : ! ======================================================================
786 : ! ================= COMPUTE THE CONJUGATE-GRADIENT =====================
787 : ! ======================================================================
788 :
789 10236671 : if (finite_field) then
790 175080 : call dotprod_g(dotgg,doti,istwf_k,npw*nspinor,1,direc,grad_total,me_g0,mpi_enreg%comm_spinorfft)
791 : !DEBUG (electric field)
792 : !check that the dotproduct is real
793 : !if (abs(doti) > tol8) then
794 : ! write(std_out,*) ' cgwf-berry: ERROR'; write(std_out,*) ' doti = ',doti
795 : ! stop
796 : !end if
797 : !ENDDEBUG
798 : else
799 10061591 : if (gen_eigenpb) then
800 4193801 : call dotprod_g(dotgg,doti,istwf_k,npw*nspinor,1,direc,direc_tmp,me_g0,mpi_enreg%comm_spinorfft)
801 : else
802 5867790 : call dotprod_g(dotgg,doti,istwf_k,npw*nspinor,1,direc,ghc,me_g0,mpi_enreg%comm_spinorfft)
803 : end if
804 : end if
805 :
806 : ! MJV: added 5 Feb 2012 - causes divide by 0 on next iteration of iline
807 10236671 : if (abs(dotgg) < TINY(0.0_dp)*1.e50_dp) dotgg = TINY(0.0_dp)*1.e50_dp
808 :
809 : ! At first iteration, gamma is set to zero
810 10236671 : if (iline==1) then
811 2901348 : gamma=zero
812 2901348 : dotgp=dotgg
813 2901348 : call cg_zcopy(npw*nspinor,direc,conjgr)
814 2901348 : if (prtvol==-level)then
815 0 : write(msg,'(a,es21.10e3)')' cgwf: dotgg = ',dotgg
816 0 : call wrtout(std_out,msg,'PERS')
817 : end if
818 :
819 : else
820 7335323 : gamma=dotgg/dotgp
821 7335323 : dotgp=dotgg
822 :
823 7335323 : if (prtvol==-level)then
824 0 : write(msg,'(a,2es16.6)')' cgwf: dotgg,gamma = ',dotgg,gamma
825 0 : call wrtout(std_out,msg,'PERS')
826 : end if
827 :
828 : ! Note: another way to compute gamma: Polak, Ribiere no real improvement ; to be more carefully tested
829 : ! call dotprod_g(dotgg,doti,istwf_k,mpi_enreg,npw*nspinor,1,direc,direc_tmp)
830 : ! !direcp must be set to zero at the beginning
831 : ! direcp=direc-direcp
832 : ! call dotprod_g(dotgmg,doti,istwf_k,mpi_enreg,npw*nspinor,1,direcp,direc_tmp)
833 : ! direcp=direc;gamma=dotgmg/dotgp;dotgp=dotgmg
834 :
835 : !$OMP PARALLEL DO
836 2568986867 : do ipw=1,npw*nspinor
837 2561651544 : conjgr(1,ipw)=direc(1,ipw)+gamma*conjgr(1,ipw)
838 2568986867 : conjgr(2,ipw)=direc(2,ipw)+gamma*conjgr(2,ipw)
839 : end do
840 : !call cg_zaxpby(npw*nspinor,cg_one,direc,(/gamma,zero/),conjgr)
841 : end if
842 :
843 : ! ======================================================================
844 : ! ============ PROJECTION OF THE CONJUGATED GRADIENT ===================
845 : ! ======================================================================
846 :
847 10236671 : if (gen_eigenpb) then
848 4274655 : call dotprod_g(dotr,doti,istwf_k,npw*nspinor,3,scwavef,conjgr,me_g0,mpi_enreg%comm_spinorfft)
849 : else
850 5962016 : call dotprod_g(dotr,doti,istwf_k,npw*nspinor,3,cwavef,conjgr,me_g0,mpi_enreg%comm_spinorfft)
851 : end if
852 :
853 : ! Project the conjugated gradient onto the current band
854 : ! MG: TODO: this is an hot spot that could be rewritten with BLAS! provided
855 : ! that direc --> conjgr
856 10236671 : if(istwf_k==1)then
857 :
858 : !$OMP PARALLEL DO
859 3032748788 : do ipw=1,npw*nspinor
860 3023273562 : direc(1,ipw)=conjgr(1,ipw)-(dotr*cwavef(1,ipw)-doti*cwavef(2,ipw))
861 3032748788 : direc(2,ipw)=conjgr(2,ipw)-(dotr*cwavef(2,ipw)+doti*cwavef(1,ipw))
862 : end do
863 : else
864 : !$OMP PARALLEL DO
865 510717479 : do ipw=1,npw*nspinor
866 509956034 : direc(1,ipw)=conjgr(1,ipw)-dotr*cwavef(1,ipw)
867 510717479 : direc(2,ipw)=conjgr(2,ipw)-dotr*cwavef(2,ipw)
868 : end do
869 : end if
870 :
871 : ! In case of generalized eigenproblem, normalization of direction vector
872 : ! cannot be done here (because S|D> is not known here).
873 10236671 : if (.not.gen_eigenpb) then
874 5962016 : call sqnorm_g(dotr,istwf_k,npw*nspinor,direc,me_g0,mpi_enreg%comm_fft)
875 5962016 : xnorm=one/sqrt(abs(dotr))
876 17886048 : call cg_zscal(npw*nspinor,(/xnorm,zero/),direc)
877 5962016 : xnorm=one
878 : end if
879 :
880 : ! ======================================================================
881 : ! ===== COMPUTE CONTRIBUTIONS TO 1ST AND 2ND DERIVATIVES OF ENERGY =====
882 : ! ======================================================================
883 :
884 : ! Compute gh_direc = <G|H|D> and eventually gs_direc = <G|S|D>
885 10236671 : sij_opt=0;if (gen_eigenpb) sij_opt=1
886 :
887 10236671 : call getghc(cpopt,direc,cprj_dum,gh_direc,gs_direc,gs_hamk,gvnlx_direc,eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,type_calc)
888 :
889 10236671 : if(wfopta10==2 .or. wfopta10==3)then
890 : ! Minimisation of the residual, so compute <G|(H-zshift)^2|D>
891 859832 : gh_direcws(:,:)=gh_direc(:,:)
892 2216 : if (gen_eigenpb) then
893 0 : sij_opt=1
894 0 : work(:,:)=gh_direc(:,:)-zshift(iband)*gs_direc(:,:)
895 : else
896 2216 : sij_opt=0
897 859832 : work(:,:)=gh_direc(:,:)-zshift(iband)*direc(:,:)
898 : end if
899 :
900 2216 : call getghc(cpopt,work,cprj_dum,gh_direc,swork,gs_hamk,gvnlx_dummy,eval,mpi_enreg,1,prtvol,0,tim_getghc,type_calc)
901 :
902 2216 : if (gen_eigenpb) then
903 0 : gh_direc(:,:)=gh_direc(:,:)-zshift(iband)*swork(:,:)
904 : else
905 859832 : gh_direc(:,:)=gh_direc(:,:)-zshift(iband)*work(:,:)
906 : end if
907 : end if
908 :
909 : ! In case of generalized eigenproblem, compute now the norm of the conjugated gradient
910 10236671 : if (gen_eigenpb) then
911 4274655 : call dotprod_g(dotr,doti,istwf_k,npw*nspinor,1,direc,gs_direc,me_g0,mpi_enreg%comm_spinorfft)
912 4274655 : xnorm=one/sqrt(abs(dotr))
913 : end if
914 :
915 : ! Compute dhc = Re{<D|H|C>}
916 10236671 : call dotprod_g(dhc,doti,istwf_k,npw*nspinor,1,direc,ghc,me_g0,mpi_enreg%comm_spinorfft)
917 10236671 : dhc=dhc*xnorm
918 :
919 : ! Compute <D|H|D> or <D|(H-zshift)^2|D>
920 10236671 : call dotprod_g(dhd,doti,istwf_k,npw*nspinor,1,direc,gh_direc,me_g0,mpi_enreg%comm_spinorfft)
921 10236671 : dhd=dhd*xnorm**2
922 :
923 10236671 : if(prtvol==-level)then
924 0 : write(msg,'(a,3f14.6)') 'cgwf: chc,dhc,dhd=',chc,dhc,dhd
925 0 : call wrtout(std_out,msg,'PERS')
926 : end if
927 :
928 : ! ======================================================================
929 : ! ======= COMPUTE MIXING FACTORS - CHECK FOR CONVERGENCE ===============
930 : ! ======================================================================
931 :
932 10236671 : if (.not.finite_field) then
933 : ! Compute tan(2 theta),sin(theta) and cos(theta)
934 10061591 : tan2th=2.0_dp*dhc/(chc-dhd)
935 :
936 10061591 : if (abs(tan2th)<1.d-05) then
937 3580694 : costh=1.0_dp-0.125_dp*tan2th**2
938 3580694 : sinth=0.5_dp*tan2th*(1.0_dp-0.375_dp*tan2th**2)
939 :
940 : ! Check that result is above machine precision
941 3580694 : if (abs(sinth)<epsilon(0._dp)) then
942 4404 : if (prtvol > 0) then
943 4187 : write(msg, '(a,es16.4)' ) ' cgwf: converged with tan2th=',tan2th
944 4187 : call wrtout(std_out,msg,'PERS')
945 : end if
946 : ! Number of one-way 3D ffts skipped
947 4404 : nskip=nskip+2*(nline-iline)
948 4404 : exit ! Exit from the loop on iline
949 : end if
950 :
951 : else
952 6480897 : root=sqrt(1.0_dp+tan2th**2)
953 6480897 : costh=sqrt(0.5_dp+0.5_dp/root)
954 6480897 : sinth=sign(sqrt(0.5_dp-0.5_dp/root),tan2th)
955 : end if
956 :
957 : ! Check for lower of two possible roots (same sign as curvature at theta where slope is zero)
958 10057187 : diff=(chc-dhd)
959 : ! Swap c and d if value of diff is positive
960 10057187 : if (diff>zero) then
961 31225 : swap=costh
962 31225 : costh=-sinth
963 31225 : sinth=swap
964 31225 : if(prtvol<0 .or. prtvol>=10)then
965 1775 : write(msg,*)' Note: swap roots, iline,diff=',iline,diff
966 1775 : call wrtout(std_out,msg,'PERS')
967 : end if
968 : end if
969 :
970 : else
971 : ! In case the electric field is on, the line minimization has to be done numerically
972 :
973 : ! Compute determinant of the overlap matrix where in the band-th line
974 : ! the wavefunction is replaced by the search direction
975 175080 : job = 10 ; shiftbd = 0
976 700320 : do idir = 1, 3
977 : ! do not do this for efield_dot(idir)=0
978 525240 : if (abs(dtefield%efield_dot(idir)) < tol12) cycle
979 1748952 : do ifor = 1, 2
980 1049248 : ikpt2f = dtefield%ikpt_dk(ikptf,ifor,idir)
981 1049248 : if (dtefield%indkk_f2ibz(ikpt2f,6) == 1) then
982 188684 : itrs = 10
983 : else
984 860564 : itrs = 0
985 : end if
986 1049248 : ikpt2 = dtefield%indkk_f2ibz(ikpt2f,1)
987 1049248 : npw_k2 = npwarr(ikpt2)
988 134708540 : pwind_k(1:npw) = pwind(ikgf+1:ikgf+npw,ifor,idir)
989 402027124 : pwnsfac_k(1:2,1:npw) = pwnsfac(1:2,ikgf+1:ikgf+npw)
990 15128228 : sflag_k(:) = dtefield%sflag(:,ikpt+(isppol-1)*nkpt,ifor,idir)
991 1171302620 : smat_k(:,:,:) = dtefield%smat(:,:,:,ikpt+(isppol-1)*nkpt,ifor,idir)
992 :
993 1049248 : if (mpi_enreg%nproc_cell > 1) then
994 529884 : icg1 = dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt)
995 : cgq_k(:,1:dtefield%mband_occ*nspinor*npw_k2) = &
996 4005752010 : & cgq(:,icg1+1:icg1+dtefield%mband_occ*nspinor*npw_k2)
997 529884 : idum1 = dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt)
998 184855848 : pwnsfac_k(3:4,1:npw_k2) = pwnsfacq(1:2,idum1+1:idum1+npw_k2)
999 : else
1000 519364 : icg1 = dtefield%cgindex(ikpt2,isppol)
1001 : cgq_k(:,1:dtefield%mband_occ*nspinor*npw_k2) = &
1002 2893814284 : & cg(:,icg1+1:icg1+dtefield%mband_occ*nspinor*npw_k2)
1003 519364 : idum1=dtefield%fkgindex(ikpt2f)
1004 217483654 : pwnsfac_k(3:4,1:npw_k2) = pwnsfac(1:2,idum1+1:idum1+npw_k2)
1005 : end if
1006 :
1007 1049248 : icg1 = 0 ; ddkflag = 0
1008 1049248 : if (gen_eigenpb) then
1009 : !$OMP PARALLEL DO
1010 39536280 : do ipw=1,npw*nspinor
1011 39051156 : direc_tmp(1,ipw)=direc(1,ipw)*xnorm
1012 39536280 : direc_tmp(2,ipw)=direc(2,ipw)*xnorm
1013 : end do
1014 : ! need cprj corresponding to direc_tmp in order to make smat_k_paw properly
1015 : call getcprj(1,0,direc_tmp,cprj_band_srt,&
1016 : & gs_hamk%ffnl_k,0,gs_hamk%indlmn,gs_hamk%istwf_k,gs_hamk%kg_k,&
1017 : & gs_hamk%kpg_k,gs_hamk%kpt_k,gs_hamk%lmnmax,gs_hamk%mgfft,&
1018 : & mpi_enreg,1,gs_hamk%natom,gs_hamk%nattyp,gs_hamk%ngfft,gs_hamk%nloalg,&
1019 : & gs_hamk%npw_k,gs_hamk%nspinor,gs_hamk%ntypat,gs_hamk%phkxred,gs_hamk%ph1d,&
1020 485124 : & gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm)
1021 :
1022 485124 : call pawcprj_copy(cprj_k,cprj_direc)
1023 : call pawcprj_put(gs_hamk%atindx,cprj_band_srt,cprj_direc,gs_hamk%natom,iband,0,&
1024 : & ikpt,1,isppol,mband,1,gs_hamk%natom,1,mband,dimlmn,gs_hamk%nspinor,nsppol,0,&
1025 485124 : & mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
1026 :
1027 : ! icp1=dtefield%mband_occ*(ikptf-1)
1028 485124 : icp2=mband*nspinor*(ikpt2-1)
1029 : call pawcprj_get(gs_hamk%atindx,cprj_kb,dtefield%cprj,natom,1,icp2,ikpt,0,isppol,&
1030 : & mband,dtefield%fnkpt,natom,mband,mband,nspinor,nsppol,0,&
1031 485124 : & mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
1032 :
1033 485124 : if (ikpt2 /= ikpt2f) then ! construct cprj_kb by symmetry
1034 482759 : call pawcprj_copy(cprj_kb,cprj_ikn)
1035 : call pawcprj_symkn(cprj_fkn,cprj_ikn,dtefield%atom_indsym,dimlmn,-1,gs_hamk%indlmn,&
1036 : & dtefield%indkk_f2ibz(ikpt2f,2),dtefield%indkk_f2ibz(ikpt2f,6),&
1037 : & dtefield%fkptns(:,dtefield%i2fbz(ikpt2)),&
1038 : & dtefield%lmax,dtefield%lmnmax,mband,natom,dtefield%mband_occ,nspinor,&
1039 482759 : & dtefield%nsym,gs_hamk%ntypat,gs_hamk%typat,dtefield%zarot)
1040 482759 : call pawcprj_copy(cprj_fkn,cprj_kb)
1041 : end if
1042 :
1043 485124 : call smatrix_k_paw(cprj_direc,cprj_kb,dtefield,idir,ifor,mband,natom,smat_k_paw,gs_hamk%typat)
1044 :
1045 : call smatrix(direc_tmp,cgq_k,cg1_k,ddkflag,dtm_k,icg1,icg1,&
1046 : & itrs,job,iband,npw*nspinor,mcg_q,mpw,iband,&
1047 : & mpw,dtefield%mband_occ,dtefield%nband_occ(isppol),&
1048 : & npw,npw_k2,nspinor,pwind_k,pwnsfac_k,sflag_k,&
1049 485124 : & shiftbd,smat_inv,smat_k,smat_k_paw,gs_hamk%usepaw)
1050 : else
1051 : call smatrix(direc,cgq_k,cg1_k,ddkflag,dtm_k,icg1,icg1,&
1052 : & itrs,job,iband,npw*nspinor,mcg_q,mpw,iband,&
1053 : & mpw,dtefield%mband_occ,dtefield%nband_occ(isppol),&
1054 : & npw,npw_k2,nspinor,pwind_k,pwnsfac_k,sflag_k,&
1055 564124 : & shiftbd,smat_inv,smat_k,smat_k_paw,gs_hamk%usepaw)
1056 : end if
1057 3672984 : detovd(:,ifor,idir) = dtm_k(:) ! Store the determinant of the overlap
1058 : ! matrix (required to compute theta_min)
1059 : ! DEBUG
1060 : ! write(std_out,*)'cgwf-berry: detovc and detovd'
1061 : ! write(std_out,*)detovc(:,ifor,idir)
1062 : ! write(std_out,*)detovd(:,ifor,idir)
1063 : ! write(std_out,*)'smat_k'
1064 : ! do jband = 1, 4
1065 : ! write(std_out,'(4(2x,e14.6))')smat_k(1,jband,:)
1066 : ! write(std_out,'(4(2x,e14.6))')smat_k(2,jband,:)
1067 : ! write(std_out,*)
1068 : ! end do
1069 : ! ENDDEBUG
1070 : end do ! ifor
1071 : end do ! idir
1072 :
1073 : call linemin(bcut,chc,costh,detovc,detovd,dhc,dhd,&
1074 : & dphase_aux1,dtefield%efield_dot,iline,&
1075 : & dtefield%fnkpt,dtefield%nstr,hel,phase_end,&
1076 175080 : & phase_init,dtefield%sdeg,sinth,thetam)
1077 : ! DEBUG
1078 : ! if (mpi_enreg%me == 1) then
1079 : ! write(std_out,*)'after linemin '
1080 : ! write(std_out,'(a,3(2x,f16.9))')'phase_init = ',phase_init(:)
1081 : ! write(std_out,'(a,3(2x,f16.9))')'phase_end = ',phase_end(:)
1082 : ! write(std_out,'(a,3(2x,f16.9))')'dphase_aux1 = ',dphase_aux1(:)
1083 : ! write(std_out,*) 'thetam',thetam
1084 : ! end if
1085 : ! ENDDEBUG
1086 : end if ! finite_field
1087 :
1088 : ! ======================================================================
1089 : ! =========== GENERATE NEW |wf>, H|wf>, Vnl|Wf>, S|Wf> ... =============
1090 : ! ======================================================================
1091 :
1092 10232267 : sintn=sinth*xnorm
1093 :
1094 : !$OMP PARALLEL DO
1095 3542332894 : do ipw=1,npw*nspinor
1096 3532100627 : cwavef(1,ipw)=cwavef(1,ipw)*costh+direc(1,ipw)*sintn
1097 3542332894 : cwavef(2,ipw)=cwavef(2,ipw)*costh+direc(2,ipw)*sintn
1098 : end do
1099 :
1100 : ! call cg_zaxpby(npw*nspinor,(/sintn,zero/),direc,(/costh,zero/),cwavef)
1101 10232267 : call cg_zcopy(npw*nspinor,cwavef,cg(1,1+icg_shift))
1102 :
1103 3542332894 : do ipw=1,npw*nspinor
1104 3532100627 : ghc(1,ipw) =ghc(1,ipw)*costh + gh_direc(1,ipw)*sintn
1105 3542332894 : ghc(2,ipw) =ghc(2,ipw)*costh + gh_direc(2,ipw)*sintn
1106 : end do
1107 :
1108 :
1109 10232267 : if (use_subvnlx==1) then
1110 : !$OMP PARALLEL DO
1111 1954943524 : do ipw=1,npw*nspinor
1112 1949021064 : gvnlxc(1,ipw)=gvnlxc(1,ipw)*costh + gvnlx_direc(1,ipw)*sintn
1113 1954943524 : gvnlxc(2,ipw)=gvnlxc(2,ipw)*costh + gvnlx_direc(2,ipw)*sintn
1114 : end do
1115 : ! call cg_zaxpby(npw*nspinor,(/sintn,zero/),gvnlx_direc,(/costh,zero/),gvnlxc)
1116 : end if
1117 :
1118 10232267 : if (gen_eigenpb) then
1119 : !$OMP PARALLEL DO
1120 1530223677 : do ipw=1,npw*nspinor
1121 1525953220 : scwavef(1,ipw)=scwavef(1,ipw)*costh+gs_direc(1,ipw)*sintn
1122 1530223677 : scwavef(2,ipw)=scwavef(2,ipw)*costh+gs_direc(2,ipw)*sintn
1123 : ! gsc(1,ipw+igsc_shift)=scwavef(1,ipw)
1124 : ! gsc(2,ipw+igsc_shift)=scwavef(2,ipw)
1125 : end do
1126 : ! call cg_zaxpby(npw*nspinor,(/sintn,zero/),gs_direc,(/costh,zero/),scwavef)
1127 4270457 : call cg_zcopy(npw*nspinor,scwavef,gsc(1,1+igsc_shift))
1128 :
1129 4270457 : if (finite_field) then ! must update cprj for the new wavefunction
1130 : call getcprj(1,0,cwavef,cprj_band_srt,&
1131 : & gs_hamk%ffnl_k,0,gs_hamk%indlmn,istwf_k,gs_hamk%kg_k,gs_hamk%kpg_k,gs_hamk%kpt_k,&
1132 : & gs_hamk%lmnmax,gs_hamk%mgfft,mpi_enreg,1,natom,gs_hamk%nattyp,&
1133 : & gs_hamk%ngfft,gs_hamk%nloalg,gs_hamk%npw_k,gs_hamk%nspinor,gs_hamk%ntypat,&
1134 80854 : & gs_hamk%phkxred,gs_hamk%ph1d,gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm)
1135 : call pawcprj_put(gs_hamk%atindx,cprj_band_srt,cprj_k,gs_hamk%natom,iband,0,ikpt,&
1136 : & 1,isppol,mband,1,gs_hamk%natom,1,mband,dimlmn,gs_hamk%nspinor,nsppol,0,&
1137 80854 : & mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
1138 : end if
1139 : end if
1140 :
1141 10232267 : if(wfopta10==2 .or. wfopta10==3)then
1142 : ! Need to keep track of ghcws, in order to avoid recomputing it
1143 : !$OMP PARALLEL DO
1144 288088 : do ipw=1,npw*nspinor
1145 285872 : ghcws(1,ipw)=ghcws(1,ipw)*costh + gh_direcws(1,ipw)*sintn
1146 288088 : ghcws(2,ipw)=ghcws(2,ipw)*costh + gh_direcws(2,ipw)*sintn
1147 : end do
1148 : ! call cg_zaxpby(npw*nspinor,(/sintn,zero/),gh_direcws,(/costh,zero/),ghcws)
1149 : end if
1150 :
1151 : ! ======================================================================
1152 : ! =========== CHECK CONVERGENCE AGAINST TRIAL ENERGY ===================
1153 : ! ======================================================================
1154 :
1155 : ! Compute delta(E)
1156 10232267 : if (.not.finite_field) then
1157 10057187 : deltae=chc*(costh**2-1._dp)+dhd*sinth**2+2._dp*costh*sinth*dhc
1158 : else
1159 : ! Compute deltae
1160 : call etheta(bcut,chc,detovc,detovd,dhc,dhd,dtefield%efield_dot,e0,e1,&
1161 175080 : & hel,dtefield%fnkpt,dtefield%nstr,dtefield%sdeg,thetam)
1162 175080 : theta = zero
1163 :
1164 : call etheta(bcut,chc,detovc,detovd,dhc,dhd,&
1165 : & dtefield%efield_dot,e0_old,e1_old,&
1166 175080 : & hel,dtefield%fnkpt,dtefield%nstr,dtefield%sdeg,theta)
1167 175080 : deltae = e0 - e0_old
1168 : ! write(std_out,*) 'e0, e0_old, deltae', e0, e0_old, deltae
1169 : ! Check that e0 is decreasing on succeeding lines:
1170 : ! if (deltae > zero) then
1171 175080 : if (deltae > tol12 .and. num_warning <= enough_warning) then ! exploring different checks for finit_field
1172 0 : num_warning = num_warning + 1
1173 : write(msg, '(3a,i8,a,1p,e14.6,a1,3x,a,1p,e14.6,a1)')&
1174 0 : ' (electric field)',ch10,&
1175 0 : ' New trial energy at line',iline,' = ',e0,ch10,&
1176 0 : ' is higher than former:',e0_old,ch10
1177 0 : ABI_WARNING(msg)
1178 : end if
1179 : end if ! finite_field
1180 :
1181 : ! Check convergence and eventually exit
1182 11717031 : if (iline==1) then
1183 2900980 : deold=deltae
1184 7331287 : else if (abs(deltae)<tolrde*abs(deold) .and. iline/=nline .and. wfopta10<2)then
1185 1338234 : if(prtvol>=10)then
1186 : write(msg, '(a,i4,1x,a,1p,e12.4,a,e12.4,a)' ) &
1187 18002 : ' cgwf: line',iline,' deltae=',deltae,' < tolrde*',deold,' =>skip lines'
1188 18002 : call wrtout(std_out,msg,'PERS')
1189 : end if
1190 1338234 : nskip=nskip+2*(nline-iline) ! Number of one-way 3D ffts skipped
1191 1338234 : exit ! Exit from the loop on iline
1192 : end if
1193 :
1194 : end do ! END LOOP FOR A GIVEN BAND Note that there are three "exit" instructions inside
1195 :
1196 : ! Additional computations in case of electric field
1197 3149627 : if (finite_field) then
1198 : ! Bring present contribution to dphasek(idir) into [-pi,pi]
1199 190080 : do idir = 1, 3
1200 142560 : dphase_aux2 = mod(phase_end(idir) - phase_init(idir) + 100*two_pi,two_pi)
1201 142560 : if (dphase_aux2 > pi) dphase_aux2 = dphase_aux2 - two_pi
1202 : ! DEBUG
1203 : ! dphase_aux1(idir)=mod(dphase_aux1(idir)+100*two_pi,two_pi)
1204 : ! if(dphase_aux1(idir)>pi) dphase_aux1(idir)=dphase_aux1(idir)-two_pi
1205 : ! diff = dphase_aux2 - dphase_aux1(idir)
1206 : ! if (abs(diff) > tol10) then
1207 : ! write(std_out,*)'cgwf-berry: ERROR'
1208 : ! write(std_out,'(a,3(2x,i3),f16.9)')'ikpt,iband,idir,diff',ikpt,iband,idir,diff
1209 : ! stop
1210 : ! end if
1211 : ! write(100,*) idir, dphase_aux2
1212 : ! ENDDEBUG
1213 190080 : dphase_k(idir) = dphase_k(idir) + dphase_aux2
1214 : ! DEBUG
1215 : ! write(std_out,*) 'idir,phase_init,phase_end,dphase_k'
1216 : ! write(std_out,*) idir,phase_init(idir),phase_end(idir),dphase_k(idir)
1217 : ! ENDDEBUG
1218 : end do
1219 : end if ! finite_field
1220 :
1221 : else ! nline==0 , needs to provide a residual
1222 0 : resid(iband)=-one
1223 : end if ! End nline==0 case
1224 :
1225 : ! ======================================================================
1226 : ! =============== END OF CURRENT BAND: CLEANING ========================
1227 : ! ======================================================================
1228 :
1229 : ! It was checked that getghc is NOT needed here : equivalent results with the copy below.
1230 3241427 : if(wfopta10==2 .or. wfopta10==3) ghc(:,:)=ghcws(:,:)
1231 :
1232 7364987 : if (finite_field) dtefield%sflag(:,ikpt + (isppol-1)*nkpt,:,:) = 0
1233 :
1234 : ! At the end of the treatment of a set of bands, write the number of one-way 3D ffts skipped
1235 3149627 : if (xmpi_paral==0 .and. mpi_enreg%paral_kgb==0 .and. iband==nband .and. prtvol/=0) then
1236 : write(msg,'(a,i0)')' cgwf: number of one-way 3D ffts skipped in cgwf until now =',nskip
1237 : call wrtout(std_out,msg,'PERS')
1238 : end if
1239 :
1240 : end do ! End big iband loop. iband in a block
1241 :
1242 : ! ======================================================================
1243 : ! ============= COMPUTE HAMILTONIAN IN WFs SUBSPACE ====================
1244 : ! ======================================================================
1245 : call mksubham(cg,ghc,gsc,gvnlxc,iblock,icg,igsc,istwf_k,&
1246 : isubh,isubo,mcg,mgsc,nband,nbdblock,npw,&
1247 3592359 : nspinor,subham,subovl,subvnlx,use_subovl,use_subvnlx,me_g0)
1248 :
1249 : end do ! iblock End loop over block of bands
1250 :
1251 442732 : ABI_SFREE(dimlmn_srt)
1252 :
1253 442732 : if (finite_field .and. gs_hamk%usepaw == 1) then ! store updated cprjs for this kpt
1254 : ! switch from ikptf to ikpt
1255 5628 : ikptf = ikpt
1256 5628 : call xmpi_allgather(ikptf,ikptf_recv,spaceComm_distrb,ierr)
1257 : call pawcprj_mpi_allgather(cprj_k,cprj_gat,natom,nspinor*mband,1,dimlmn,ncpgr,nproc_distrb,&
1258 5628 : & spaceComm_distrb,ierr,rank_ordered=.true.)
1259 18600 : do iproc = 1, nproc_distrb
1260 12972 : icp2=nspinor*mband*(iproc-1)
1261 : call pawcprj_get(gs_hamk%atindx1,cprj_k,cprj_gat,natom,1,icp2,ikpt,0,isppol,&
1262 : & mband,nproc_distrb,natom,mband,mband,nspinor,nsppol,0,&
1263 12972 : & mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
1264 : ! ikptf = ikptf_recv(iproc)
1265 12972 : icp1 = nspinor*mband*(ikptf_recv(iproc)-1)
1266 : call pawcprj_put(gs_hamk%atindx1,cprj_k,dtefield%cprj,natom,1,icp1,ikpt,0,isppol,&
1267 : & mband,nkpt,natom,mband,mband,dimlmn,nspinor,nsppol,0,&
1268 18600 : & mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
1269 : end do
1270 : end if
1271 :
1272 : ! Debugging outputs
1273 442732 : if(prtvol==-level)then
1274 0 : isubh=1
1275 0 : if (use_subvnlx==1) write(msg,'(a)') ' cgwf : isubh subham(isubh:isubh+1) subvnlx(isubh:isubh+1)'
1276 0 : if (use_subvnlx==0) write(msg,'(a)') ' cgwf : isubh subham(isubh:isubh+1)'
1277 0 : do iband=1,nband
1278 0 : do ii=1,iband
1279 0 : if (use_subvnlx==1) then
1280 0 : write(msg,'(i5,4es16.6)')isubh,subham(isubh:isubh+1),subvnlx(isubh:isubh+1)
1281 : else
1282 0 : write(msg,'(i5,2es16.6)')isubh,subham(isubh:isubh+1)
1283 : end if
1284 0 : call wrtout(std_out,msg,'PERS')
1285 0 : isubh=isubh+2
1286 : end do
1287 : end do
1288 : end if
1289 :
1290 : ! ===================
1291 : ! FINAL DEALLOCATIONS
1292 : ! ===================
1293 442732 : ABI_FREE(conjgr)
1294 442732 : ABI_FREE(cwavef)
1295 442732 : ABI_FREE(direc)
1296 442732 : ABI_FREE(pcon)
1297 442732 : ABI_FREE(scprod)
1298 442732 : ABI_FREE(ghc)
1299 442732 : ABI_FREE(gvnlxc)
1300 442732 : ABI_FREE(gh_direc)
1301 442732 : ABI_FREE(gvnlx_direc)
1302 442732 : ABI_FREE(vresid)
1303 442732 : ABI_FREE(gs_direc)
1304 442732 : ABI_FREE(gvnlx_dummy)
1305 442732 : ABI_FREE(swork)
1306 :
1307 442732 : ABI_SFREE(scwavef)
1308 442732 : ABI_SFREE(direc_tmp)
1309 442732 : ABI_SFREE(ghc_all)
1310 442732 : ABI_SFREE(ghcws)
1311 442732 : ABI_SFREE(gh_direcws)
1312 442732 : ABI_SFREE(work)
1313 :
1314 442732 : if (finite_field) then
1315 8330 : ABI_FREE(cg1_k)
1316 8330 : ABI_FREE(cgq_k)
1317 8330 : ABI_FREE(detovc)
1318 8330 : ABI_FREE(detovd)
1319 8330 : ABI_FREE(grad_berry)
1320 8330 : ABI_FREE(sflag_k)
1321 8330 : ABI_FREE(smat_inv)
1322 8330 : ABI_FREE(smat_k)
1323 8330 : ABI_FREE(pwind_k)
1324 8330 : ABI_FREE(pwnsfac_k)
1325 8330 : ABI_FREE(grad_total)
1326 8330 : if (gs_hamk%usepaw /= 0) then
1327 5628 : call pawcprj_free(cprj_k)
1328 5628 : call pawcprj_free(cprj_kb)
1329 5628 : call pawcprj_free(cprj_direc)
1330 5628 : call pawcprj_free(cprj_band_srt)
1331 5628 : call pawcprj_free(cprj_gat)
1332 5628 : if (nkpt /= dtefield%fnkpt) then
1333 5628 : call pawcprj_free(cprj_fkn)
1334 5628 : call pawcprj_free(cprj_ikn)
1335 50652 : ABI_FREE(cprj_fkn)
1336 50652 : ABI_FREE(cprj_ikn)
1337 : end if
1338 : end if
1339 8330 : ABI_FREE(smat_k_paw)
1340 8330 : ABI_FREE(dimlmn)
1341 53354 : ABI_FREE(cprj_k)
1342 53354 : ABI_FREE(cprj_kb)
1343 53354 : ABI_FREE(cprj_direc)
1344 112106 : ABI_FREE(cprj_gat)
1345 8330 : ABI_FREE(ikptf_recv)
1346 19586 : ABI_FREE(cprj_band_srt)
1347 : end if
1348 :
1349 442732 : gs_hamk%gpu_option=gpu_option_save
1350 442732 : if(gpu_option_save==ABI_GPU_OPENMP) gemm_nonlop_use_gemm=.true.
1351 :
1352 : ! Do not delete this line, needed to run with open MP
1353 442732 : write(unit=msg,fmt=*) resid(1)
1354 :
1355 442732 : call timab(22,2,tsec)
1356 :
1357 : DBG_EXIT("COLL")
1358 :
1359 885464 : end subroutine cgwf
1360 : !!***
1361 :
1362 : !!****f* m_cgwf/linemin
1363 : !! NAME
1364 : !! linemin
1365 : !!
1366 : !! FUNCTION
1367 : !! Performs the "line minimization" w.r.t. the angle theta on a unit circle
1368 : !! to update the wavefunction associated with the current k-point and
1369 : !! band label.
1370 : !! This routine is used only when the electric field is on (otherwise it could
1371 : !! in principle also be used, but there is a simpler procedure, as originally
1372 : !! coded in abinit).
1373 : !!
1374 : !! INPUTS
1375 : !! chc = <C|H_0|C> where |C> is the wavefunction of the current band
1376 : !! detovc = determinant of the overlap matrix S
1377 : !! detovd = determinant of the overlap matrix where for the band
1378 : !! that is being updated <C| is replaced by <D| (search direction)
1379 : !! dhc = Re[<D|H_0|C>]
1380 : !! dhd = <D|H_0|D>
1381 : !! efield_dot = reciprocal lattice coordinates of the electric field
1382 : !! iline = index of the current line minimization
1383 : !! nkpt = number of k-points
1384 : !! nstr(idir) = number of strings along the idir-th direction
1385 : !! sdeg = spin degeneracy
1386 : !!
1387 : !! OUTPUT
1388 : !! bcut(ifor,idir) = branch cut of the ellipse associated with (ifor,idir)
1389 : !! costh = cos(thetam)
1390 : !! hel(ifor,idir) = helicity of the ellipse associated with (ifor,idir)
1391 : !! phase_end = total change in Zak phase, must be equal to
1392 : !! dphase_aux1 + n*two_pi
1393 : !! sinth = sin(thetam)
1394 : !! thetam = optimal angle theta in line_minimization
1395 : !!
1396 : !! SIDE EFFECTS
1397 : !! Input/Output
1398 : !! dphase_aux1 = change in Zak phase accumulated during the loop over iline
1399 : !! (can be used for debugging in cgwf.f)
1400 : !! phase_init = initial Zak phase (before doing the first line minimization)
1401 : !!
1402 : !! NOTES
1403 : !! We are making the "frozen Hamiltonian approximation", i.e., the
1404 : !! Hamiltonian does not change with theta (we are neglecting the dependence
1405 : !! of the Hartree and exchange-correlation terms on theta; the original
1406 : !! abinit routine does the same)
1407 : !!
1408 : !! SOURCE
1409 :
1410 175080 : subroutine linemin(bcut,chc,costh,detovc,detovd,dhc,dhd,dphase_aux1,&
1411 : & efield_dot,iline,nkpt,nstr,hel,phase_end,phase_init,sdeg,sinth,thetam)
1412 :
1413 : !Arguments ------------------------------------
1414 : !scalars
1415 : integer,intent(in) :: iline,nkpt
1416 : real(dp),intent(in) :: chc,dhc,dhd,sdeg
1417 : real(dp),intent(out) :: costh,sinth,thetam
1418 : !arrays
1419 : integer,intent(in) :: nstr(3)
1420 : integer,intent(out) :: hel(2,3)
1421 : real(dp),intent(in) :: detovc(2,2,3),detovd(2,2,3),efield_dot(3)
1422 : real(dp),intent(inout) :: dphase_aux1(3),phase_init(3)
1423 : real(dp),intent(out) :: bcut(2,3),phase_end(3)
1424 :
1425 : !Local variables -------------------------
1426 : !scalars
1427 : integer :: idir,ifor,igrid,iter,maxiter,ngrid
1428 : real(dp) :: aa,angle,bb,big_axis,cc,cphi_0,delta_theta,e0,e1
1429 : real(dp) :: excentr,iab,phase0,phase_min,phi_0,rdum,sgn,small_axis,sphi_0
1430 : real(dp) :: theta,theta_0,val
1431 : logical :: flag_neg
1432 : character(len=500) :: msg
1433 : !arrays
1434 : real(dp) :: g_theta(2),theta_min(2),theta_try(2)
1435 : real(dp) :: esave(251),e1save(251) !!REC
1436 :
1437 : ! ***********************************************************************
1438 :
1439 : !Compute the helicity and the branch cut of the ellipse in the complex
1440 : !plane associated with the overlap between a k-point and one of its neighbours
1441 :
1442 700320 : do idir = 1, 3
1443 :
1444 525240 : if (abs(efield_dot(idir)) < tol12) cycle
1445 :
1446 1748952 : do ifor = 1, 2
1447 :
1448 : aa = half*(detovc(1,ifor,idir)*detovc(1,ifor,idir) + &
1449 : & detovc(2,ifor,idir)*detovc(2,ifor,idir) + &
1450 : & detovd(1,ifor,idir)*detovd(1,ifor,idir) + &
1451 1049248 : & detovd(2,ifor,idir)*detovd(2,ifor,idir))
1452 :
1453 : bb = half*(detovc(1,ifor,idir)*detovc(1,ifor,idir) + &
1454 : & detovc(2,ifor,idir)*detovc(2,ifor,idir) - &
1455 : & detovd(1,ifor,idir)*detovd(1,ifor,idir) - &
1456 1049248 : & detovd(2,ifor,idir)*detovd(2,ifor,idir))
1457 :
1458 : cc = detovc(1,ifor,idir)*detovd(1,ifor,idir) + &
1459 1049248 : & detovc(2,ifor,idir)*detovd(2,ifor,idir)
1460 :
1461 : iab = detovc(1,ifor,idir)*detovd(2,ifor,idir) - &
1462 1049248 : & detovc(2,ifor,idir)*detovd(1,ifor,idir)
1463 :
1464 1049248 : if (iab >= zero) then
1465 522698 : hel(ifor,idir) = 1
1466 : else
1467 526550 : hel(ifor,idir) = -1
1468 : end if
1469 :
1470 1049248 : if (abs(bb) > tol8) then
1471 1049248 : theta_0 = half*atan(cc/bb)
1472 : else
1473 : theta_0 = quarter*pi
1474 : end if
1475 :
1476 1049248 : if (bb < zero) theta_0 = theta_0 + pi*half
1477 :
1478 3147744 : g_theta(:) = cos(theta_0)*detovc(:,ifor,idir) + sin(theta_0)*detovd(:,ifor,idir)
1479 : ! write(std_out,*)'before rhophi, g_theta =',g_theta
1480 1049248 : call rhophi(g_theta,phi_0,rdum)
1481 : ! write(std_out,*)'after rhophi, phi_0 = ',phi_0
1482 :
1483 1049248 : cphi_0 = cos(phi_0)
1484 1049248 : sphi_0 = sin(phi_0)
1485 :
1486 1049248 : rdum = aa - sqrt(bb*bb + cc*cc)
1487 1049248 : if (rdum < zero) rdum = zero
1488 1049248 : small_axis = sqrt(rdum)
1489 1049248 : big_axis = sqrt(aa + sqrt(bb*bb + cc*cc))
1490 1049248 : excentr = hel(ifor,idir)*small_axis/big_axis
1491 :
1492 : ! Find angle for which phi = pi
1493 1049248 : if (abs(excentr) > tol8) then
1494 1049178 : angle = atan(tan(pi-phi_0)/excentr)
1495 : else
1496 70 : if (tan(pi-phi_0)*hel(ifor,idir) > zero) then
1497 : angle = half*pi
1498 : else
1499 26 : angle = -0.5_dp*pi
1500 : end if
1501 : end if
1502 1049248 : bcut(ifor,idir) = angle + theta_0
1503 :
1504 :
1505 : ! Compute the branch-cut angle
1506 1049248 : if (hel(ifor,idir) == 1) then
1507 522698 : if ((sphi_0 > 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) + pi
1508 522698 : if ((sphi_0 < 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) - pi
1509 : else
1510 526550 : if ((sphi_0 > 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) - pi
1511 526550 : if ((sphi_0 < 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) + pi
1512 : end if
1513 :
1514 1049248 : if (bcut(ifor,idir) > pi) bcut(ifor,idir) = bcut(ifor,idir) - two_pi
1515 2623736 : if (bcut(ifor,idir) < -1_dp*pi) bcut(ifor,idir) = bcut(ifor,idir) + two_pi
1516 :
1517 : ! DEBUG
1518 : ! write(std_out,'(a,2x,i3,2x,i3,5x,f16.9,5x,i2)')'linemin: ifor,idir,bcut,hel',&
1519 : ! & ifor,idir,bcut(ifor,idir),hel(ifor,idir)
1520 : ! write(std_out,'(a,5x,f16.9,5x,f16.9)')'linemin: big_axis,small_axis ',&
1521 : ! & big_axis,small_axis
1522 : ! ENDDEBUG
1523 :
1524 : end do ! ifor
1525 : end do ! idir
1526 :
1527 : !---------------------------------------------------------------------------
1528 :
1529 : !Perform the "line minimization" w.r.t. the angle theta on a unit circle
1530 : !to update the wavefunction associated with the current k-point and band label.
1531 :
1532 : ngrid = 250 ! initial number of subdivisions in [-pi/2,pi/2]
1533 : !for finding extrema
1534 : maxiter = 100
1535 525240 : delta_theta = pi/ngrid
1536 :
1537 : !DEBUG
1538 : !write(std_out,*)'linemin: theta, e0, e1, e1fdiff'
1539 : !ENDDEBUG
1540 :
1541 :
1542 : !Get the interval where the absolute minimum of E(theta) is located
1543 :
1544 525240 : val = huge(one) ! large number
1545 525240 : flag_neg=.false.
1546 525240 : theta_min(:) = ten
1547 44120160 : do igrid = 1, ngrid+1
1548 :
1549 43945080 : theta = (igrid - 1)*delta_theta - pi*half
1550 : call etheta(bcut,chc,detovc,detovd,dhc,dhd,efield_dot,e0,e1,&
1551 43945080 : & hel,nkpt,nstr,sdeg,theta)
1552 :
1553 43945080 : esave(igrid)=e0 !!REC
1554 43945080 : e1save(igrid)=e1 !!REC
1555 :
1556 : ! It is important to detect when the slope changes from negative to positive
1557 : ! Moreover, a slope being extremely close to zero must be ignored
1558 :
1559 : ! DEBUG
1560 : ! write(std_out,*)' igrid,e0,e1,val,theta_min(:)=',igrid,theta,e0,e1,val,theta_min(:)
1561 : ! ENDDEBUG
1562 :
1563 : ! Store e1 and theta if negative ...
1564 43945080 : if(e1 < -tol10)then
1565 21897104 : theta_try(1)=theta
1566 21897104 : flag_neg=.true.
1567 : end if
1568 : ! A change of sign is just happening
1569 88065240 : if(e1 > tol10 .and. flag_neg)then
1570 209837 : theta_try(2)=theta
1571 209837 : flag_neg=.false.
1572 : ! Still, must be better than the previous minimum in order to succeed
1573 209837 : if (e0 < val-tol10) then
1574 191589 : val=e0
1575 191589 : theta_min(:)=theta_try(:)
1576 : end if
1577 : end if
1578 : end do
1579 :
1580 : !In case the minimum was not found
1581 :
1582 175080 : if (abs(theta_min(1) - ten) < tol10) then
1583 : ! REC start
1584 0 : write(msg,'(a,a)')ch10,' linemin: ERROR- cannot find theta_min.'
1585 0 : call wrtout(std_out,msg,'COLL')
1586 0 : write(msg,'(a,a)')ch10,' igrid theta esave(igrid) e1save(igrid) '
1587 0 : call wrtout(std_out,msg,'COLL')
1588 0 : do igrid = 1, ngrid+1
1589 0 : theta = (igrid - 1)*delta_theta - pi*half
1590 0 : write(std_out,'(i6,3f16.9)')igrid,theta,esave(igrid),e1save(igrid)
1591 : !write(101,'(i6,3f16.9)')igrid,theta,esave(igrid),e1save(igrid)
1592 : end do
1593 0 : write(msg,'(6a)')ch10,&
1594 0 : ' linemin: ERROR - ',ch10,&
1595 0 : ' Cannot find theta_min. No minimum exists: the field is too strong ! ',ch10,&
1596 0 : ' Try decreasing difference between D and 4 Pi P by changing structure or D (only for fixed D calculation)'
1597 0 : call wrtout(std_out,msg,'COLL')
1598 0 : ABI_ERROR('linemin cannot find theta_min')
1599 : end if
1600 :
1601 : !Compute the mimum of E(theta)
1602 :
1603 :
1604 : iter = 0
1605 3853907 : do while ((delta_theta > tol8).and.(iter < maxiter))
1606 3678827 : delta_theta = half*(theta_min(2) - theta_min(1))
1607 3678827 : theta = theta_min(1) + delta_theta
1608 : call etheta(bcut,chc,detovc,detovd,dhc,dhd,efield_dot,e0,e1,&
1609 3678827 : & hel,nkpt,nstr,sdeg,theta)
1610 3678827 : if (e1 > zero) then
1611 1282472 : theta_min(2) = theta
1612 : else
1613 2396355 : theta_min(1) = theta
1614 : end if
1615 3678827 : iter = iter + 1
1616 :
1617 : ! DEBUG
1618 : ! write(std_out,'(a,2x,i3,2(2x,f16.9))')'iter,e0,e1 = ',iter,e0,e1
1619 : ! ENDDEBUG
1620 :
1621 : end do
1622 :
1623 175080 : costh = cos(theta)
1624 175080 : sinth = sin(theta)
1625 :
1626 175080 : thetam = theta
1627 :
1628 : !DEBUG
1629 : !write(std_out,*)'linemin : thetam = ',thetam
1630 : !ENDDEBUG
1631 :
1632 : !---------------------------------------------------------------------------
1633 :
1634 : !Compute and store the change in electronic polarization
1635 :
1636 175080 : sgn = one
1637 700320 : do idir = 1, 3
1638 :
1639 525240 : if (abs(efield_dot(idir)) < tol12) cycle
1640 :
1641 524624 : phase_end(idir) = zero
1642 1748952 : do ifor = 1, 2
1643 :
1644 3147744 : g_theta(:) = detovc(:,ifor,idir)
1645 : ! write(std_out,*)'before rhophi (2nd call), g_theta =',g_theta
1646 1049248 : call rhophi(g_theta,phase0,rdum)
1647 : ! write(std_out,*)'after rhophi, phase0 = ',phase0
1648 :
1649 1049248 : if(iline == 1) phase_init(idir) = phase_init(idir) + sgn*phase0
1650 :
1651 3147744 : g_theta(:) = costh*detovc(:,ifor,idir) + sinth*detovd(:,ifor,idir)
1652 1049248 : call rhophi(g_theta,phase_min,rdum)
1653 :
1654 1049248 : phase_end(idir) = phase_end(idir) + sgn*phase_min
1655 :
1656 : ! Correct for branch cuts (remove jumps)
1657 1049248 : if (bcut(ifor,idir) <= zero) phase0 = phase0 + hel(ifor,idir)*two_pi
1658 1049248 : if(thetam >= bcut(ifor,idir)) phase_min = phase_min + hel(ifor,idir)*two_pi
1659 :
1660 1049248 : dphase_aux1(idir) = dphase_aux1(idir) + sgn*(phase_min - phase0)
1661 :
1662 3672984 : sgn = -1_dp*sgn
1663 :
1664 : end do ! idir
1665 : end do ! ifor
1666 :
1667 : !DEBUG
1668 : !write(std_out,'(a,3(2x,f16.9))')'dphase_aux1 = ',(dphase_aux1(idir),idir = 1, 3)
1669 : !write(std_out,*)' linemin: debug, exit.'
1670 : !ENDDEBUG
1671 :
1672 175080 : end subroutine linemin
1673 : !!***
1674 :
1675 : !!****f* m_cgwf/etheta
1676 : !! NAME
1677 : !! etheta
1678 : !!
1679 : !! FUNCTION
1680 : !! Computes the energy per unit cell and its first derivative
1681 : !! for a given angle theta. More precisely, computes only the part of
1682 : !! the energy that changes with theta.
1683 : !!
1684 : !! INPUTS
1685 : !! bcut(ifor,idir) = branch cut of the ellipse associated with (ifor,idir)
1686 : !! chc = <C|H_0|C> where |C> is the wavefunction of the current band
1687 : !! detovc = determinant of the overlap matrix S
1688 : !! detovd = determinant of the overlap matrix where for the band
1689 : !! that is being updated <C| is replaced by <D| (search direction)
1690 : !! dhc = Re[<D|H_0|C>]
1691 : !! dhd = <D|H_0|D>
1692 : !! efield_dot = reciprocal lattice coordinates of the electric field
1693 : !! hel(ifor,idir) = helicity of the ellipse associated with (ifor,idir)
1694 : !! nkpt = number of k-points
1695 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
1696 : !! nstr(idir) = number of strings along the idir-th direction
1697 : !! sdeg = spin degeneracy
1698 : !! theta = value of the angle for which the energy (e0) and its
1699 : !! derivative (e1) are computed
1700 : !!
1701 : !! OUTPUT
1702 : !! e0 = energy for the given value of theta
1703 : !! e1 = derivative of the energy with respect to theta
1704 : !!
1705 : !! SOURCE
1706 :
1707 47974067 : subroutine etheta(bcut,chc,detovc,detovd,dhc,dhd,efield_dot,e0,e1,&
1708 : & hel,nkpt,nstr,sdeg,theta)
1709 :
1710 : !Arguments ------------------------------------
1711 : !scalars
1712 : integer,intent(in) :: nkpt
1713 : real(dp),intent(in) :: chc,dhc,dhd,sdeg,theta
1714 : real(dp),intent(out) :: e0,e1
1715 : !arrays
1716 : integer,intent(in) :: hel(2,3),nstr(3)
1717 : real(dp),intent(in) :: bcut(2,3),detovc(2,2,3),detovd(2,2,3),efield_dot(3)
1718 :
1719 : !Local variables -------------------------
1720 : !scalars
1721 : integer :: idir,ifor
1722 : real(dp) :: c2theta,ctheta,dphase,gnorm,phase,rho,s2theta,sgn,stheta
1723 : !arrays
1724 : real(dp) :: dg_theta(2),g_theta(2)
1725 :
1726 : ! ***********************************************************************
1727 :
1728 : e0 = zero ; e1 = zero
1729 :
1730 47974067 : ctheta = cos(theta)
1731 47974067 : stheta = sin(theta)
1732 47974067 : c2theta = ctheta*ctheta - stheta*stheta ! cos(2*theta)
1733 47974067 : s2theta = two*ctheta*stheta ! sin(2*theta)
1734 :
1735 47974067 : e0 = chc*ctheta*ctheta + dhd*stheta*stheta + dhc*s2theta
1736 47974067 : e0 = e0*sdeg/nkpt
1737 :
1738 : !DEBUG
1739 : !e0 = zero
1740 : !ENDDEBUG
1741 :
1742 47974067 : e1 = (dhd - chc)*s2theta + two*dhc*c2theta
1743 47974067 : e1 = e1*sdeg/nkpt
1744 :
1745 47974067 : sgn = -1_dp
1746 191896268 : do idir = 1, 3
1747 :
1748 143922201 : if (abs(efield_dot(idir)) < tol12) cycle
1749 :
1750 479234315 : do ifor = 1, 2
1751 :
1752 : g_theta(:) = ctheta*detovc(:,ifor,idir) + &
1753 862520496 : & stheta*detovd(:,ifor,idir)
1754 : dg_theta(:) = -1_dp*stheta*detovc(:,ifor,idir) + &
1755 862520496 : & ctheta*detovd(:,ifor,idir)
1756 :
1757 : ! Compute E(theta)
1758 :
1759 287506832 : call rhophi(g_theta,phase,rho)
1760 287506832 : if (theta >= bcut(ifor,idir)) phase = phase + hel(ifor,idir)*two_pi
1761 :
1762 : ! DEBUG
1763 : ! unit = 100 + 10*idir + ifor; write(unit,'(4(f16.9))')theta,g_theta(:),phase
1764 : ! ENDDEBUG
1765 :
1766 287506832 : e0 = e0 + sgn*sdeg*efield_dot(idir)*phase/(two_pi*nstr(idir))
1767 :
1768 :
1769 : ! Compute dE/dtheta
1770 :
1771 : ! imaginary part of the derivative of ln(g_theta)
1772 287506832 : gnorm = g_theta(1)*g_theta(1) + g_theta(2)*g_theta(2)
1773 287506832 : dphase = (dg_theta(2)*g_theta(1) - dg_theta(1)*g_theta(2))/gnorm
1774 :
1775 287506832 : e1 = e1 + sgn*sdeg*efield_dot(idir)*dphase/(two_pi*nstr(idir))
1776 :
1777 718935865 : sgn = -1_dp*sgn
1778 :
1779 : end do
1780 : end do
1781 :
1782 47974067 : end subroutine etheta
1783 : !!***
1784 :
1785 : !!****f* m_cgwf/mksubham
1786 : !! NAME
1787 : !! mksubham
1788 : !!
1789 : !! FUNCTION
1790 : !! Build the Hamiltonian matrix in the eigenfunctions subspace,
1791 : !! for one given band (or for one given block of bands)
1792 : !!
1793 : !! INPUTS
1794 : !! cg(2,mcg)=wavefunctions
1795 : !! gsc(2,mgsc)=<g|S|c> matrix elements (S=overlap)
1796 : !! iblock=index of block of bands
1797 : !! icg=shift to be applied on the location of data in the array cg
1798 : !! igsc=shift to be applied on the location of data in the array cg
1799 : !! istwf_k=input parameter that describes the storage of wfs
1800 : !! mcg=second dimension of the cg array
1801 : !! mgsc=second dimension of the gsc array
1802 : !! nband_k=number of bands at this k point for that spin polarization
1803 : !! nbdblock=number of bands in a block
1804 : !! npw_k=number of plane waves at this k point
1805 : !! nspinor=number of spinorial components of the wavefunctions
1806 : !! use_subovl=1 if the overlap matrix is not identity in WFs subspace
1807 : !! use_subvnlx= 1 if <C band,k|H|C band_prime,k> has to be computed
1808 : !! me_g0=1 if this processors has G=0, 0 otherwise
1809 : !!
1810 : !! OUTPUT
1811 : !!
1812 : !! SIDE EFFECTS
1813 : !! ghc(2,npw_k*nspinor)=<G|H|C band,k> for the current state
1814 : !! This is an input in non-blocked algorithm
1815 : !! an output in blocked algorithm
1816 : !! gvnlxc(2,npw_k*nspinor)=<G|Vnl|C band,k> for the current state
1817 : !! This is an input in non-blocked algorithm
1818 : !! an output in blocked algorithm
1819 : !! isubh=index of current state in array subham
1820 : !! isubo=index of current state in array subovl
1821 : !! subham(nband_k*(nband_k+1))=Hamiltonian expressed in the WFs subspace
1822 : !! subovl(nband_k*(nband_k+1)*use_subovl)=overlap matrix expressed in the WFs subspace
1823 : !! subvnlx(nband_k*(nband_k+1)*use_subvnlx)=non-local Hamiltonian (if NCPP) plus Fock ACE operator (if usefock_ACE)
1824 : !! expressed in the WFs subspace
1825 : !!
1826 : !! SOURCE
1827 :
1828 3149627 : subroutine mksubham(cg,ghc,gsc,gvnlxc,iblock,icg,igsc,istwf_k,&
1829 : & isubh,isubo,mcg,mgsc,nband_k,nbdblock,npw_k,&
1830 3149627 : & nspinor,subham,subovl,subvnlx,use_subovl,use_subvnlx,me_g0)
1831 :
1832 : !Arguments ------------------------------------
1833 : !scalars
1834 : integer,intent(in) :: iblock,icg,igsc,istwf_k,mcg,mgsc,nband_k
1835 : integer,intent(in) :: nbdblock,npw_k,nspinor,use_subovl,use_subvnlx,me_g0
1836 : integer,intent(inout) :: isubh,isubo
1837 : !arrays
1838 : real(dp),intent(in) :: cg(2,mcg)
1839 : real(dp),intent(in) :: gsc(2,mgsc)
1840 : real(dp),intent(inout) :: ghc(2,npw_k*nspinor),gvnlxc(2,npw_k*nspinor)
1841 : real(dp),intent(inout) :: subham(nband_k*(nband_k+1))
1842 : real(dp),intent(inout) :: subovl(nband_k*(nband_k+1)*use_subovl)
1843 : real(dp),intent(inout) :: subvnlx(nband_k*(nband_k+1)*use_subvnlx)
1844 :
1845 : !Local variables-------------------------------
1846 : !scalars
1847 : integer :: iband,ibdblock,ii,ipw,ipw1,isp,iwavef,jwavef
1848 : real(dp) :: cgimipw,cgreipw,chcim,chcre,cscim,cscre,cvcim,cvcre
1849 : !real(dp) :: chc(2),cvc(2),csc(2)
1850 :
1851 : ! *********************************************************************
1852 :
1853 : !Loop over bands in a block This loop can be parallelized
1854 6299254 : do iband=1+(iblock-1)*nbdblock,min(iblock*nbdblock,nband_k)
1855 3149627 : ibdblock=iband-(iblock-1)*nbdblock
1856 :
1857 : ! Compute elements of subspace Hamiltonian <C(i)|H|C(n)> and <C(i)|Vnl|C(n)>
1858 3149627 : if(istwf_k==1)then
1859 :
1860 29990204 : do ii=1,iband
1861 27077335 : iwavef=(ii-1)*npw_k*nspinor+icg
1862 27077335 : chcre=zero ; chcim=zero
1863 27077335 : if (use_subvnlx==0) then
1864 4205249040 : do ipw=1,npw_k*nspinor
1865 4193373153 : cgreipw=cg(1,ipw+iwavef)
1866 4193373153 : cgimipw=cg(2,ipw+iwavef)
1867 4193373153 : chcre=chcre+cgreipw*ghc(1,ipw)+cgimipw*ghc(2,ipw)
1868 4205249040 : chcim=chcim+cgreipw*ghc(2,ipw)-cgimipw*ghc(1,ipw)
1869 : end do
1870 : ! chc = cg_zdotc(npw_k*nspinor,cg(1,1+iwavef),ghc)
1871 : else
1872 : #if 1
1873 5851386581 : do ipw=1,npw_k*nspinor
1874 5836185133 : cgreipw=cg(1,ipw+iwavef)
1875 5836185133 : cgimipw=cg(2,ipw+iwavef)
1876 5836185133 : chcre=chcre+cgreipw*ghc(1,ipw)+cgimipw*ghc(2,ipw)
1877 5851386581 : chcim=chcim+cgreipw*ghc(2,ipw)-cgimipw*ghc(1,ipw)
1878 : end do
1879 : cvcre=zero ; cvcim=zero
1880 5851386581 : do ipw=1,npw_k*nspinor
1881 5836185133 : cgreipw=cg(1,ipw+iwavef)
1882 5836185133 : cgimipw=cg(2,ipw+iwavef)
1883 5836185133 : cvcre=cvcre+cgreipw*gvnlxc(1,ipw)+cgimipw*gvnlxc(2,ipw)
1884 5851386581 : cvcim=cvcim+cgreipw*gvnlxc(2,ipw)-cgimipw*gvnlxc(1,ipw)
1885 : end do
1886 15201448 : subvnlx(isubh )=cvcre
1887 15201448 : subvnlx(isubh+1)=cvcim
1888 : #else
1889 : ! New version with BLAS1, will require some update of the refs.
1890 : cvc = cg_zdotc(npw_k*nspinor,cg(1,1+iwavef),gvnlxc)
1891 : subvnlx(isubh )=cvc(1)
1892 : subvnlx(isubh+1)=cvc(2)
1893 : chc = cg_zdotc(npw_k*nspinor,cg(1,1+iwavef),ghc)
1894 : chcre = chc(1)
1895 : chcim = chc(2)
1896 : #endif
1897 : ! Store real and imag parts in Hermitian storage mode:
1898 : end if
1899 27077335 : subham(isubh )=chcre
1900 27077335 : subham(isubh+1)=chcim
1901 : ! subham(isubh )=chc(1)
1902 : ! subham(isubh+1)=chc(2)
1903 29990204 : isubh=isubh+2
1904 : end do
1905 :
1906 236758 : else if(istwf_k>=2)then
1907 3024374 : do ii=1,iband
1908 2787616 : iwavef=(ii-1)*npw_k+icg
1909 : ! Use the time-reversal symmetry, but should not double-count G=0
1910 2787616 : if(istwf_k==2 .and. me_g0==1) then
1911 1028770 : chcre = half*cg(1,1+iwavef)*ghc(1,1)
1912 1028770 : if (use_subvnlx==1) cvcre=half*cg(1,1+iwavef)*gvnlxc(1,1)
1913 : ipw1=2
1914 : else
1915 1758846 : chcre=zero; ipw1=1
1916 1758846 : if (use_subvnlx==1) cvcre=zero
1917 : end if
1918 2787616 : if (use_subvnlx==0) then
1919 2091306 : do isp=1,nspinor
1920 992545161 : do ipw=ipw1+(isp-1)*npw_k,npw_k*isp
1921 990453855 : cgreipw=cg(1,ipw+iwavef)
1922 990453855 : cgimipw=cg(2,ipw+iwavef)
1923 991499508 : chcre=chcre+cgreipw*ghc(1,ipw)+cgimipw*ghc(2,ipw)
1924 : end do
1925 : end do
1926 1045653 : chcre=two*chcre
1927 : else
1928 3483926 : do isp=1,nspinor
1929 983592409 : do ipw=ipw1+(isp-1)*npw_k,npw_k*isp
1930 980108483 : cgreipw=cg(1,ipw+iwavef)
1931 980108483 : cgimipw=cg(2,ipw+iwavef)
1932 980108483 : chcre=chcre+cgreipw*ghc(1,ipw)+cgimipw*ghc(2,ipw)
1933 981850446 : cvcre=cvcre+cgreipw*gvnlxc(1,ipw)+cgimipw*gvnlxc(2,ipw)
1934 : end do
1935 : end do
1936 1741963 : chcre=two*chcre
1937 1741963 : cvcre=two*cvcre
1938 : ! Store real and imag parts in Hermitian storage mode:
1939 1741963 : subvnlx(isubh )=cvcre
1940 1741963 : subvnlx(isubh+1)=zero
1941 : end if
1942 2787616 : subham(isubh )=chcre
1943 2787616 : subham(isubh+1)=zero
1944 3024374 : isubh=isubh+2
1945 : end do
1946 : end if
1947 :
1948 : ! Compute elements of subspace <C(i)|S|C(n)> (S=overlap matrix)
1949 : ! <C(i)|S|C(n)> should be closed to Identity.
1950 6299254 : if (use_subovl==1) then
1951 0 : jwavef=(iband-1)*npw_k*nspinor+igsc
1952 0 : if(istwf_k==1)then
1953 0 : do ii=1,iband
1954 0 : iwavef=(ii-1)*npw_k*nspinor+icg
1955 0 : cscre=zero ; cscim=zero
1956 0 : do ipw=1,npw_k*nspinor
1957 0 : cgreipw=cg(1,ipw+iwavef)
1958 0 : cgimipw=cg(2,ipw+iwavef)
1959 0 : cscre=cscre+cgreipw*gsc(1,ipw+jwavef)+cgimipw*gsc(2,ipw+jwavef)
1960 0 : cscim=cscim+cgreipw*gsc(2,ipw+jwavef)-cgimipw*gsc(1,ipw+jwavef)
1961 : end do
1962 : ! csc = cg_zdotc(npw_k*nspinor,cg(1,1+iwavef),gsc)
1963 : ! subovl(isubo )=csc(1)
1964 : ! subovl(isubo+1)=csc(2)
1965 : ! Store real and imag parts in Hermitian storage mode:
1966 0 : subovl(isubo )=cscre
1967 0 : subovl(isubo+1)=cscim
1968 0 : isubo=isubo+2
1969 : end do
1970 0 : else if(istwf_k>=2)then
1971 0 : do ii=1,iband
1972 0 : iwavef=(ii-1)*npw_k*nspinor+icg
1973 0 : if(istwf_k==2 .and. me_g0==1)then
1974 0 : cscre=half*cg(1,1+iwavef)*gsc(1,1+jwavef)
1975 0 : ipw1=2
1976 : else
1977 : cscre=zero; ipw1=1
1978 : end if
1979 0 : do isp=1,nspinor
1980 0 : do ipw=ipw1+(isp-1)*npw_k,npw_k*isp
1981 0 : cgreipw=cg(1,ipw+iwavef)
1982 0 : cgimipw=cg(2,ipw+iwavef)
1983 0 : cscre=cscre+cg(1,ipw+iwavef)*gsc(1,ipw+jwavef)+cg(2,ipw+iwavef)*gsc(2,ipw+jwavef)
1984 : end do
1985 : end do
1986 0 : cscre=two*cscre
1987 : ! Store real and imag parts in Hermitian storage mode:
1988 0 : subovl(isubo )=cscre
1989 0 : subovl(isubo+1)=zero
1990 0 : isubo=isubo+2
1991 : end do
1992 : end if
1993 : end if
1994 :
1995 : end do ! iband in a block
1996 :
1997 3149627 : end subroutine mksubham
1998 : !!***
1999 :
2000 : !!****f* ABINIT/make_grad_berry
2001 : !! NAME
2002 : !! make_grad_berry
2003 : !!
2004 : !! FUNCTION
2005 : !! compute gradient contribution from berry phase in finite
2006 : !! electric field case
2007 : !!
2008 : !! INPUTS
2009 : !! cg(2,mcg)=input wavefunctions
2010 : !! cgq(2,mcgq) = wavefunctions at neighboring k points
2011 : !! cprj_k(natom,nband_k*usepaw)=cprj at this k point
2012 : !! dimlmn(natom)=lmn_size for each atom in input order
2013 : !! dimlmn_srt(natom)=lmn_size for each atom sorted by type
2014 : !! direc(2,npw*nspinor)=gradient vector
2015 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
2016 : !! iband=index of band currently being treated
2017 : !! icg=shift to be applied on the location of data in the array cg
2018 : !! ikpt=number of the k-point currently being treated
2019 : !! isppol=spin polarization currently treated
2020 : !! natom=number of atoms in cell.
2021 : !! mband =maximum number of bands
2022 : !! mpw=maximum dimensioned size of npw
2023 : !! mcg=second dimension of the cg array
2024 : !! mcgq=second dimension of the cgq array
2025 : !! mkgq = second dimension of pwnsfacq
2026 : !! nkpt=number of k points
2027 : !! mpi_enreg=information about MPI parallelization
2028 : !! npw=number of planewaves in basis sphere at given k.
2029 : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
2030 : !! nsppol=number of spin polarizations
2031 : !! pwind(pwind_alloc,2,3) = array used to compute
2032 : !! the overlap matrix smat between k-points (see initberry.f)
2033 : !! pwind_alloc = first dimension of pwind
2034 : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
2035 : !! (see initberry.f)
2036 : !! pwnsfacq(2,mkgq) = phase factors for the nearest neighbours of the
2037 : !! current k-point (electric field, MPI //)
2038 : !!
2039 : !! OUTPUT
2040 : !! grad_berry(2,npw*nspinor) :: contribution to gradient in finite electric field case
2041 : !!
2042 : !! SIDE EFFECTS
2043 : !! dtefield <type(efield_type)> = variables related to Berry phase calculations (see initberry.f)
2044 : !!
2045 : !! NOTES
2046 : !!
2047 : !! SOURCE
2048 :
2049 175080 : subroutine make_grad_berry(cg,cgq,cprj_k,detovc,dimlmn,dimlmn_srt,direc,dtefield,grad_berry,&
2050 : & gs_hamk,iband,icg,ikpt,isppol,mband,mcg,mcgq,mkgq,mpi_enreg,mpw,natom,nkpt,&
2051 175080 : & npw,npwarr,nspinor,nsppol,pwind,pwind_alloc,pwnsfac,pwnsfacq)
2052 :
2053 : !Arguments ------------------------------------
2054 : !scalars
2055 : integer,intent(in) :: iband,icg,ikpt,isppol,mband,mcg,mcgq
2056 : integer,intent(in) :: mkgq,mpw,natom,nkpt,npw,nspinor,nsppol,pwind_alloc
2057 : type(gs_hamiltonian_type),intent(in) :: gs_hamk
2058 : type(efield_type),intent(inout) :: dtefield
2059 : type(MPI_type),intent(in) :: mpi_enreg
2060 :
2061 : !arrays
2062 : integer,intent(in) :: dimlmn(natom),dimlmn_srt(natom)
2063 : integer,intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3)
2064 : real(dp),intent(in) :: cg(2,mcg),cgq(2,mcgq)
2065 : real(dp),intent(inout) :: direc(2,npw*nspinor)
2066 : real(dp),intent(in) :: pwnsfac(2,pwind_alloc),pwnsfacq(2,mkgq)
2067 : real(dp),intent(out) :: detovc(2,2,3),grad_berry(2,npw*nspinor)
2068 : type(pawcprj_type),intent(in) :: cprj_k(natom,dtefield%mband_occ*gs_hamk%usepaw*dtefield%nspinor)
2069 :
2070 : !Local variables-------------------------------
2071 : !scalars
2072 : integer :: choice,cpopt,ddkflag,dimenlr1,iatom,icg1,icp2,idum1
2073 : integer :: idir,ifor,ikgf,ikptf,ikpt2,ikpt2f,ipw,i_paw_band,ispinor,itrs,itypat,job
2074 : integer :: klmn,mcg1_k,mcg_q,nbo,npw_k2,nspinortot,paw_opt,shiftbd,signs
2075 : real(dp) :: fac
2076 : character(len=500) :: msg
2077 : !arrays
2078 525240 : integer :: pwind_k(npw),sflag_k(dtefield%mband_occ)
2079 525240 : real(dp) :: cg1_k(2,npw*nspinor),dtm_k(2),pwnsfac_k(4,mpw)
2080 350160 : real(dp) :: smat_k(2,dtefield%mband_occ,dtefield%mband_occ)
2081 350160 : real(dp) :: smat_inv(2,dtefield%mband_occ,dtefield%mband_occ),svectout_dum(2,0)
2082 : real(dp) :: dummy_enlout(0)
2083 175080 : real(dp),allocatable :: cgq_k(:,:),enl_rij(:,:,:,:),grad_berry_ev(:,:)
2084 175080 : real(dp),allocatable :: qijbkk(:,:,:,:),smat_k_paw(:,:,:)
2085 : ! type(pawcprj_type) :: cprj_dum(1,1) ! was used in on-site dipole, now suppressed
2086 : ! 15 June 2012 J Zwanziger
2087 175080 : type(pawcprj_type),allocatable :: cprj_kb(:,:),cprj_band_srt(:,:)
2088 175080 : type(pawcprj_type),allocatable :: cprj_fkn(:,:),cprj_ikn(:,:)
2089 : ! *********************************************************************
2090 :
2091 175080 : nbo = dtefield%mband_occ
2092 :
2093 : !allocations
2094 :
2095 : !Electric field: compute the gradient of the Berry phase part of the energy functional.
2096 : !See PRL 89, 117602 (2002) [[cite:Souza2002]], grad_berry(:,:) is the second term of Eq. (4)
2097 72453204 : grad_berry(:,:) = zero
2098 175080 : job = 11 ; shiftbd = 1
2099 175080 : mcg_q = mpw*mband*nspinor
2100 175080 : mcg1_k = npw*nspinor
2101 :
2102 175080 : if (gs_hamk%usepaw /= 0) then
2103 80854 : dimenlr1 = gs_hamk%lmnmax*(gs_hamk%lmnmax+1)/2
2104 485124 : ABI_MALLOC(qijbkk,(dimenlr1,natom,nspinor**2,2))
2105 404270 : ABI_MALLOC(enl_rij,(nspinor*dimenlr1,natom,nspinor**2,1))
2106 323416 : ABI_MALLOC(smat_k_paw,(2,nbo,nbo))
2107 242562 : ABI_MALLOC(grad_berry_ev,(2,npw*nspinor))
2108 6225758 : enl_rij = zero
2109 12370662 : qijbkk = zero
2110 4285262 : smat_k_paw = zero
2111 1293664 : ABI_MALLOC(cprj_kb,(natom,nbo*nspinor))
2112 80854 : call pawcprj_alloc(cprj_kb,0,dimlmn)
2113 565978 : ABI_MALLOC(cprj_band_srt,(natom,nspinor))
2114 80854 : call pawcprj_alloc(cprj_band_srt,0,dimlmn_srt)
2115 80854 : if (nkpt /= dtefield%fnkpt) then
2116 1212810 : ABI_MALLOC(cprj_fkn,(natom,nbo*nspinor))
2117 1212810 : ABI_MALLOC(cprj_ikn,(natom,nbo*nspinor))
2118 80854 : call pawcprj_alloc(cprj_fkn,0,dimlmn)
2119 80854 : call pawcprj_alloc(cprj_ikn,0,dimlmn)
2120 : else
2121 0 : ABI_MALLOC(cprj_fkn,(0,0))
2122 0 : ABI_MALLOC(cprj_ikn,(0,0))
2123 : end if
2124 : else
2125 94226 : ABI_MALLOC(qijbkk,(0,0,0,0))
2126 94226 : ABI_MALLOC(enl_rij,(0,0,0,0))
2127 94226 : ABI_MALLOC(smat_k_paw,(0,0,0))
2128 94226 : ABI_MALLOC(grad_berry_ev,(0,0))
2129 94226 : ABI_MALLOC(cprj_kb,(0,0))
2130 94226 : ABI_MALLOC(cprj_band_srt,(0,0))
2131 94226 : ABI_MALLOC(cprj_fkn,(0,0))
2132 94226 : ABI_MALLOC(cprj_ikn,(0,0))
2133 : end if
2134 :
2135 175080 : ikptf = dtefield%i2fbz(ikpt)
2136 175080 : ikgf = dtefield%fkgindex(ikptf) ! this is the shift for pwind
2137 :
2138 700320 : do idir = 1, 3
2139 : ! skip idir values for which efield_dot(idir)=0
2140 525240 : if (abs(dtefield%efield_dot(idir)) < tol12) cycle
2141 : ! Implicitly, we use the gradient multiplied by the number of k points in the FBZ
2142 : fac = dtefield%efield_dot(idir)*dble(dtefield%fnkpt)/&
2143 524624 : & (dble(dtefield%nstr(idir))*four_pi)
2144 1748952 : do ifor = 1, 2
2145 : ! Handle dtefield%i2fbz properly and ask whether t.r.s. is used
2146 1049248 : ikpt2f = dtefield%ikpt_dk(ikptf,ifor,idir)
2147 1049248 : if (dtefield%indkk_f2ibz(ikpt2f,6) == 1) then
2148 188684 : itrs = 10
2149 : else
2150 860564 : itrs = 0
2151 : end if
2152 1049248 : ikpt2 = dtefield%indkk_f2ibz(ikpt2f,1)
2153 1049248 : npw_k2 = npwarr(ikpt2)
2154 3147744 : ABI_MALLOC(cgq_k,(2,nbo*nspinor*npw_k2))
2155 134708540 : pwind_k(1:npw) = pwind(ikgf+1:ikgf+npw,ifor,idir)
2156 402027124 : pwnsfac_k(1:2,1:npw) = pwnsfac(1:2,ikgf+1:ikgf+npw)
2157 15128228 : sflag_k(:) = dtefield%sflag(:,ikpt+(isppol-1)*nkpt,ifor,idir)
2158 1171302620 : smat_k(:,:,:) = dtefield%smat(:,:,:,ikpt+(isppol-1)*nkpt,ifor,idir)
2159 1049248 : if (mpi_enreg%nproc_cell > 1) then
2160 529884 : icg1 = dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt)
2161 : cgq_k(:,1:nbo*nspinor*npw_k2) = &
2162 4005752010 : & cgq(:,icg1+1:icg1+nbo*nspinor*npw_k2)
2163 529884 : idum1 = dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt)
2164 184855848 : pwnsfac_k(3:4,1:npw_k2) = pwnsfacq(1:2,idum1+1:idum1+npw_k2)
2165 : else
2166 519364 : icg1 = dtefield%cgindex(ikpt2,isppol)
2167 : cgq_k(:,1:nbo*nspinor*npw_k2) = &
2168 2893814284 : & cg(:,icg1+1:icg1+nbo*nspinor*npw_k2)
2169 519364 : idum1 = dtefield%fkgindex(ikpt2f)
2170 217483654 : pwnsfac_k(3:4,1:npw_k2) = pwnsfac(1:2,idum1+1:idum1+npw_k2)
2171 : end if
2172 1049248 : if (gs_hamk%usepaw == 1) then
2173 485124 : icp2=nbo*(ikpt2-1)*nspinor
2174 : call pawcprj_get(gs_hamk%atindx1,cprj_kb,dtefield%cprj,natom,1,icp2,ikpt,0,isppol,&
2175 : & nbo,dtefield%fnkpt,natom,nbo,nbo,nspinor,nsppol,0,&
2176 485124 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2177 485124 : if (ikpt2 /= ikpt2f) then ! construct cprj_kb by symmetry
2178 482759 : call pawcprj_copy(cprj_kb,cprj_ikn)
2179 : call pawcprj_symkn(cprj_fkn,cprj_ikn,dtefield%atom_indsym,dimlmn,-1,gs_hamk%indlmn,&
2180 : & dtefield%indkk_f2ibz(ikpt2f,2),dtefield%indkk_f2ibz(ikpt2f,6),&
2181 : & dtefield%fkptns(:,dtefield%i2fbz(ikpt2)),&
2182 : & dtefield%lmax,dtefield%lmnmax,mband,natom,nbo,nspinor,&
2183 482759 : & dtefield%nsym,gs_hamk%ntypat,gs_hamk%typat,dtefield%zarot)
2184 482759 : call pawcprj_copy(cprj_fkn,cprj_kb)
2185 : end if
2186 485124 : call smatrix_k_paw(cprj_k,cprj_kb,dtefield,idir,ifor,mband,natom,smat_k_paw,gs_hamk%typat)
2187 : end if
2188 :
2189 1049248 : icg1 = 0 ; ddkflag = 1
2190 : call smatrix(cg,cgq_k,cg1_k,ddkflag,dtm_k,icg,icg1,itrs,&
2191 : & job,iband,mcg,mcg_q,mcg1_k,iband,mpw,nbo,dtefield%nband_occ(isppol),&
2192 : & npw,npw_k2,nspinor,pwind_k,pwnsfac_k,sflag_k,&
2193 1049248 : & shiftbd,smat_inv,smat_k,smat_k_paw,gs_hamk%usepaw)
2194 1049248 : ABI_FREE(cgq_k)
2195 3147744 : detovc(:,ifor,idir) = dtm_k(:) !store the determinant of the overlap
2196 1049248 : if (sqrt(dtm_k(1)*dtm_k(1) + dtm_k(2)*dtm_k(2)) < tol12) then
2197 : write(msg,'(3a,i5,a,i3,a,a,a)') &
2198 0 : ' (electric field)',ch10,&
2199 0 : ' For k-point #',ikpt,' and band # ',iband,',',ch10,&
2200 0 : ' the determinant of the overlap matrix is found to be 0. Fixing...'
2201 : ! REC try this:
2202 0 : write(std_out,*)msg,dtm_k(1:2)
2203 0 : if(abs(dtm_k(1))<=1d-12)dtm_k(1)=1d-12
2204 0 : if(abs(dtm_k(2))<=1d-12)dtm_k(2)=1d-12
2205 0 : write(std_out,*)' Changing to:',dtm_k(1:2)
2206 : end if
2207 :
2208 1049248 : if (gs_hamk%usepaw == 1) then
2209 : ! this loop applies discretized derivative of projectors
2210 : ! note that qijb_kk is sorted by input atom order, but nonlop wants it sorted by type
2211 1455372 : do iatom = 1, natom
2212 970248 : itypat = gs_hamk%typat(gs_hamk%atindx1(iatom))
2213 36384300 : do klmn = 1, dtefield%lmn2_size(itypat)
2214 : ! note: D_ij-like terms have 4 spinor components: 11, 22, 12, and 21. Here the qijb is diagonal
2215 : ! in spin space so only the first two are nonzero and they are equal
2216 70828104 : do ispinor = 1, nspinor
2217 34928928 : qijbkk(klmn,iatom,ispinor,1) = dtefield%qijb_kk(1,klmn,gs_hamk%atindx1(iatom),idir)
2218 34928928 : qijbkk(klmn, iatom,ispinor,2) = dtefield%qijb_kk(2,klmn,gs_hamk%atindx1(iatom),idir)
2219 69857856 : if (ifor > 1) qijbkk(klmn,iatom,ispinor,2) = -qijbkk(klmn,iatom,ispinor,2)
2220 : end do
2221 : end do ! end loop over lmn2_size
2222 : end do ! end loop over natom
2223 :
2224 485124 : choice = 1
2225 485124 : signs = 2
2226 485124 : paw_opt = 1
2227 485124 : cpopt = 2 ! use cprj_kb in memory
2228 485124 : nspinortot=min(2,nspinor*(1+mpi_enreg%paral_spinor))
2229 2425620 : do i_paw_band = 1, nbo
2230 :
2231 : call pawcprj_get(gs_hamk%atindx,cprj_band_srt,cprj_kb,natom,i_paw_band,0,ikpt,1,&
2232 : & isppol,nbo,1,natom,1,nbo,nspinor,nsppol,0,&
2233 1940496 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2234 :
2235 : ! Pass dummy_enlout to avoid aliasing (enl, enlout)
2236 : call nonlop(choice,cpopt,cprj_band_srt,dummy_enlout,gs_hamk,idir,(/zero/),mpi_enreg,1,0,&
2237 1940496 : & paw_opt,signs,svectout_dum,0,direc,grad_berry_ev,enl=qijbkk)
2238 :
2239 : ! Add i*fac*smat_inv(i_paw_band,iband)*grad_berry_ev to the gradient
2240 158630244 : do ipw = 1, npw*nspinor
2241 :
2242 : grad_berry(1,ipw) = grad_berry(1,ipw) - &
2243 : & fac*(smat_inv(2,i_paw_band,iband)*grad_berry_ev(1,ipw) + &
2244 156204624 : & smat_inv(1,i_paw_band,iband)*grad_berry_ev(2,ipw))
2245 :
2246 : grad_berry(2,ipw) = grad_berry(2,ipw) + &
2247 : & fac*(smat_inv(1,i_paw_band,iband)*grad_berry_ev(1,ipw) - &
2248 158145120 : & smat_inv(2,i_paw_band,iband)*grad_berry_ev(2,ipw))
2249 :
2250 : end do
2251 : end do
2252 : end if ! end if PAW
2253 :
2254 : ! Add i*fac*cg1_k to the gradient
2255 145554332 : do ipw = 1, npw*nspinor
2256 144505084 : grad_berry(1,ipw) = grad_berry(1,ipw) - fac*cg1_k(2,ipw)
2257 145554332 : grad_berry(2,ipw) = grad_berry(2,ipw) + fac*cg1_k(1,ipw)
2258 : end do
2259 1049248 : fac = -1._dp*fac
2260 15128228 : dtefield%sflag(:,ikpt+(isppol-1)*nkpt,ifor,idir) = sflag_k(:)
2261 1049248 : dtefield%sflag(iband,ikpt+(isppol-1)*nkpt,ifor,idir) = 0
2262 1171827860 : dtefield%smat(:,:,:,ikpt+(isppol-1)*nkpt,ifor,idir) = smat_k(:,:,:)
2263 : end do ! ifor
2264 :
2265 : ! if (gs_hamk%usepaw == 1) then
2266 : ! ! call nonlop to apply on-site dipole <EV> part to direc
2267 : ! ! note that rij is sorted by input atom order, but nonlop wants it sorted by type
2268 : ! do iatom = 1, natom
2269 : ! itypat = gs_hamk%typat(gs_hamk%atindx1(iatom))
2270 : ! do klmn = 1, dtefield%lmn2_size(itypat)
2271 : ! ! note: D_ij-like terms have 4 spinor components: 11, 22, 12, and 21. Here the enl_rij is diagonal
2272 : ! ! in spin space so only the first two are nonzero and they are equal
2273 : ! do ispinor = 1, nspinor
2274 : ! if (nspinor == 1) then
2275 : ! enl_rij(klmn,iatom,ispinor) = dtefield%rij(klmn,itypat,idir)
2276 : ! else
2277 : ! enl_rij(2*klmn-1,iatom,ispinor) = dtefield%rij(klmn,itypat,idir)
2278 : ! end if
2279 : ! end do
2280 : ! end do ! end loop over lmn2_size
2281 : ! end do ! end loop over natom
2282 : ! cpopt = -1 ! compute cprj inside nonlop because we do not have them for direc
2283 : ! call nonlop(choice,cpopt,cprj_dum,dummy_enlout,gs_hamk,idir,zero,mpi_enreg,1,0,&
2284 : ! & paw_opt,signs,svectout_dum,0,direc,grad_berry_ev,enl=enl_rij)
2285 : ! grad_berry(:,:) = grad_berry(:,:) - dtefield%efield_dot(idir)*grad_berry_ev(:,:)/two_pi
2286 : ! end if
2287 :
2288 : end do ! idir
2289 :
2290 : !deallocations
2291 175080 : if(gs_hamk%usepaw /= 0) then
2292 80854 : call pawcprj_free(cprj_kb)
2293 80854 : call pawcprj_free(cprj_band_srt)
2294 80854 : if (nkpt /= dtefield%fnkpt) then
2295 80854 : call pawcprj_free(cprj_fkn)
2296 80854 : call pawcprj_free(cprj_ikn)
2297 : end if
2298 : end if
2299 175080 : ABI_FREE(grad_berry_ev)
2300 175080 : ABI_FREE(qijbkk)
2301 175080 : ABI_FREE(enl_rij)
2302 175080 : ABI_FREE(smat_k_paw)
2303 821912 : ABI_FREE(cprj_kb)
2304 336788 : ABI_FREE(cprj_band_srt)
2305 821912 : ABI_FREE(cprj_fkn)
2306 821912 : ABI_FREE(cprj_ikn)
2307 :
2308 175080 : end subroutine make_grad_berry
2309 : !!***
2310 :
2311 : !!****f* m_cgwf/nscf_init
2312 : !! NAME
2313 : !! nscf_init
2314 : !!
2315 : !! FUNCTION
2316 : !! Initialize the object.
2317 : !!
2318 : !! INPUT
2319 : !! dtset<dataset_type>=All input variables for this dataset.
2320 : !! dtfil<datafiles_type>=Variables related to files.
2321 : !! cryst=Crystalline structure
2322 : !! comm=MPI communicator.
2323 : !!
2324 : !! SOURCE
2325 :
2326 4 : subroutine nscf_init(nscf, dtset, dtfil, cryst, comm)
2327 :
2328 : !Arguments ------------------------------------
2329 : class(nscf_t),intent(out) :: nscf
2330 : type(dataset_type),intent(in) :: dtset
2331 : type(datafiles_type),intent(in) :: dtfil
2332 : type(crystal_t),intent(in) :: cryst
2333 : integer,intent(in) :: comm
2334 :
2335 : !Local variables ------------------------------
2336 : !scalars
2337 : integer,parameter :: cplex1 = 1, pawread0 = 0
2338 : integer :: fform, nfftf
2339 104 : type(crystal_t) :: pot_cryst
2340 2 : type(hdr_type) :: pot_hdr
2341 : character(len=500) :: msg
2342 : !arrays
2343 : integer :: units(2)
2344 2 : type(pawrhoij_type),allocatable :: pot_pawrhoij(:)
2345 : ! *************************************************************************
2346 :
2347 6 : units = [std_out, ab_out]
2348 2 : ABI_CHECK(dtset%usepaw == 0, "PAW not implemented!")
2349 2 : if (dtset%usekden /= 0) then
2350 0 : ABI_ERROR("nscf_init with MGGA not yet coded")
2351 : end if
2352 :
2353 2 : call wrtout(units, sjoin(" Reading KS GS potential from: ", dtfil%filpotin))
2354 2 : call pot_hdr%from_fname(dtfil%filpotin, fform, comm)
2355 2 : ABI_CHECK(fform /= 0, "hdr_read_from_fname returned fform 0")
2356 2 : ABI_CHECK(fform_contains(fform, "vtrial", msg), msg)
2357 :
2358 : ! Init FFT mesh from file as we don't want to interpolate the KS potential.
2359 2 : call ngfft_seq(nscf%ngfftf, pot_hdr%ngfft)
2360 2 : call ngfft_seq(nscf%ngfft, pot_hdr%ngfft)
2361 2 : call pot_hdr%free()
2362 :
2363 2 : if (dtset%usepaw == 0) then
2364 14 : ABI_CHECK(all(nscf%ngfft(1:6) == dtset%ngfft(1:6)), "different ngfft meshes")
2365 : else
2366 0 : ABI_CHECK(all(nscf%ngfftf(1:6) == dtset%ngfftdg(1:6)), "different ngfftf meshes")
2367 : end if
2368 :
2369 2 : call initmpi_seq(nscf%mpi_enreg)
2370 2 : call nscf%mpi_enreg%distribfft%init_seq('c', nscf%ngfft(2), nscf%ngfft(3), 'all')
2371 2 : call nscf%mpi_enreg%distribfft%init_seq('f', nscf%ngfftf(2), nscf%ngfftf(3), 'all')
2372 :
2373 : ! Read KS potential from file.
2374 8 : nfftf = product(nscf%ngfftf(1:3))
2375 8 : ABI_MALLOC(nscf%vtrial, (nfftf, dtset%nspden))
2376 :
2377 : call read_rhor(dtfil%filpotin, cplex1, dtset%nspden, nfftf, nscf%ngfftf, pawread0, nscf%mpi_enreg, &
2378 2 : nscf%vtrial, pot_hdr, pot_pawrhoij, comm, allow_interp=.False., want_varname="vtrial")
2379 :
2380 2 : pot_cryst = pot_hdr%get_crystal()
2381 2 : if (cryst%compare(pot_cryst, header=" Comparing input crystal with POT crystal") /= 0) then
2382 0 : ABI_ERROR("Crystal structure from WFK and POT files do not agree! Check messages above!")
2383 : end if
2384 2 : call pot_cryst%free(); call pot_hdr%free()
2385 :
2386 2 : end subroutine nscf_init
2387 : !!***
2388 :
2389 : !!****f* m_cgwf/nscf_setup_spin
2390 : !! NAME
2391 : !! nscf_setup_spin
2392 : !!
2393 : !! FUNCTION
2394 : !!
2395 : !! INPUT
2396 : !! isppol=Spin index.
2397 : !! dtset<dataset_type>=All input variables for this dataset.
2398 : !! pawfgr=<type(pawfgr_type)>=fine grid parameters and related data
2399 : !! gs_ham_k=<type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
2400 : !!
2401 : !! OUTPUT
2402 : !! vlocal: KS potential for this spin.
2403 : !!
2404 : !! SOURCE
2405 :
2406 0 : subroutine nscf_setup_spin(nscf, isppol, dtset, pawfgr, gs_ham_k, vlocal)
2407 :
2408 : !Arguments ------------------------------------
2409 : class(nscf_t),intent(inout) :: nscf
2410 : integer,intent(in) :: isppol
2411 : type(dataset_type),intent(in) :: dtset
2412 : type(pawfgr_type),intent(in) :: pawfgr
2413 : type(gs_hamiltonian_type),intent(inout) :: gs_ham_k
2414 : !arrays
2415 : real(dp),allocatable,intent(out) :: vlocal(:,:,:,:)
2416 :
2417 : !Local variables ------------------------------
2418 : integer :: nvloc, n1, n2, n3, n4, n5, n6, nfft, nfftf, mgfft, mgfftf
2419 : ! *************************************************************************
2420 :
2421 : ! See vtorho.F90 for the sequence of calls needed to initialize the GS Hamiltonian.
2422 : ! The Hamiltonian has references to the _k arrays allocated here and returned
2423 : associate (mpi_enreg => nscf%mpi_enreg)
2424 :
2425 : !==== Initialize most of the Hamiltonian ====
2426 : ! Allocate all arrays and initialize quantities that do not depend on k and spin.
2427 :
2428 : ! FFT meshes from input file, not necessary equal to the ones found in the external files.
2429 0 : nfftf = product(nscf%ngfftf(1:3)); mgfftf = maxval(nscf%ngfftf(1:3))
2430 0 : nfft = product(nscf%ngfft(1:3)) ; mgfft = maxval(nscf%ngfft(1:3))
2431 0 : n1 = nscf%ngfft(1); n2 = nscf%ngfft(2); n3 = nscf%ngfft(3); n4 = nscf%ngfft(4); n5 = nscf%ngfft(5); n6 = nscf%ngfft(6)
2432 :
2433 : ! Set up local potential vlocal on the coarse FFT mesh from vtrial taking into account the spin.
2434 : ! Also, continue to initialize the Hamiltonian.
2435 0 : nvloc = gs_ham_k%nvloc
2436 0 : ABI_CALLOC(vlocal, (n4, n5, n6, nvloc))
2437 :
2438 : call gspot_transgrid_and_pack(isppol, dtset%usepaw, nscf%paral_kgb0, nfft, nscf%ngfft, nfftf, &
2439 0 : dtset%nspden, gs_ham_k%nvloc, 1, pawfgr, mpi_enreg, nscf%vtrial, vlocal)
2440 :
2441 0 : call gs_ham_k%load_spin(isppol, vlocal=vlocal, with_nonlocal=.true.)
2442 :
2443 0 : if (dtset%usekden /= 0) then
2444 0 : ABI_ERROR("nscf_init with mgga not yet coded")
2445 : !call gspot_transgrid_and_pack(isppol, dtset%usepaw, paral_kgb0, dtset%nfft, dtset%ngfft, nfftf, &
2446 : ! dtset%nspden, gs_ham_k%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal)
2447 : !call gs_ham_k%load_spin(isppol, vxctaulocal=vxctaulocal)
2448 : end if
2449 :
2450 : end associate
2451 :
2452 0 : end subroutine nscf_setup_spin
2453 : !!***
2454 :
2455 : !!****f* m_cgwf/nscf_setup_kpt
2456 : !! NAME
2457 : !! nscf_setup_kpt
2458 : !!
2459 : !! FUNCTION
2460 : !! Prepare call to nscf_solve_kpt.
2461 : !! Compute k-dependent terms, gs_ham_k and allocate wavefunction block for this k-point.
2462 : !!
2463 : !! INPUT
2464 : !! isppol=Spin index
2465 : !! kpt(3)=K-point
2466 : !! istwf_k=wavefunction storage.
2467 : !! nband_k=Number of bands.
2468 : !! cryst=Crystalline structure.
2469 : !! dtset<dataset_type>=All input variables for this dataset.
2470 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
2471 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
2472 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data.
2473 : !!
2474 : !! OUTPUT
2475 : !! kg_k=
2476 : !! cg_k
2477 : !! gsc_k
2478 : !! eig_k
2479 : !! gs_ham_k=all data for the Hamiltonian at k-point kpt.
2480 : !!
2481 : !! SOURCE
2482 :
2483 127 : subroutine nscf_setup_kpt(nscf, isppol, kpt, istwf_k, nband_k, cryst, dtset, psps, pawtab, pawfgr, & ! in
2484 : npw_k, kg_k, kpg_k, ph3d_k, kinpw_k, ffnl_k, vlocal, cg_k, gsc_k, gs_ham_k) ! out
2485 :
2486 : !Arguments ------------------------------------
2487 : class(nscf_t),intent(inout) :: nscf
2488 : integer,intent(in) :: isppol, istwf_k, nband_k
2489 : real(dp),intent(in) :: kpt(3)
2490 : type(dataset_type),intent(in) :: dtset
2491 : type(crystal_t),intent(in) :: cryst
2492 : type(pseudopotential_type),intent(in) :: psps
2493 : type(pawtab_type),intent(in) :: pawtab(cryst%ntypat*psps%usepaw)
2494 : type(pawfgr_type),intent(in) :: pawfgr
2495 : type(gs_hamiltonian_type),intent(out) :: gs_ham_k
2496 : !arrays
2497 : integer,intent(out) :: npw_k
2498 : integer,allocatable,intent(out) :: kg_k(:,:)
2499 : real(dp),allocatable,intent(out) :: kpg_k(:,:), ph3d_k(:,:,:), kinpw_k(:), ffnl_k(:,:,:,:), vlocal(:,:,:,:)
2500 : real(dp),allocatable,intent(out) :: cg_k(:,:,:), gsc_k(:,:,:)
2501 :
2502 : !Local variables ------------------------------
2503 : !scalars
2504 : integer,parameter :: nkpt1 = 1, use_subovl0 = 0, ider0 = 0, idir0 = 0, mkmem1 = 1, useylmgr0 = 0, optder0=0
2505 : integer :: nvloc, nkpg, n1, n2, n3, n4, n5, n6, nfft, nfftf, mgfft, mgfftf, nspinor, ncomp
2506 : !character(len=500) :: msg
2507 : !arrays
2508 : real(dp) :: ylmgr_dum(1,1,1)
2509 127 : real(dp),allocatable :: ph1d(:,:), ylm_k(:,:)
2510 : ! *************************************************************************
2511 :
2512 : ABI_CHECK_IEQ(gs_ham_k%use_gbt, 0, "use_gbt /= 0 not coded")
2513 :
2514 : ! See vtorho.F90 for the sequence of calls needed to initialize the GS Hamiltonian.
2515 : ! The Hamiltonian has references to the _k arrays that allocated here and returned to the caller.
2516 : associate (mpi_enreg => nscf%mpi_enreg)
2517 :
2518 : !==== Initialize most of the Hamiltonian ====
2519 : ! Allocate all arrays and initialize quantities that do not depend on k and spin.
2520 :
2521 : ! FFT meshes from input file, not necessary equal to the ones found in the external files.
2522 508 : nfftf = product(nscf%ngfftf(1:3)); mgfftf = maxval(nscf%ngfftf(1:3))
2523 889 : nfft = product(nscf%ngfft(1:3)) ; mgfft = maxval(nscf%ngfft(1:3))
2524 127 : n1 = nscf%ngfft(1); n2 = nscf%ngfft(2); n3 = nscf%ngfft(3); n4 = nscf%ngfft(4); n5 = nscf%ngfft(5); n6 = nscf%ngfft(6)
2525 127 : nspinor = dtset%nspinor
2526 :
2527 : ! Compute g-sphere for this k-point from ecut
2528 127 : call get_kg(kpt, istwf_k, dtset%ecut, cryst%gmet, npw_k, kg_k)
2529 :
2530 : ! Compute kinetic energy for this k-point.
2531 381 : ABI_MALLOC(kinpw_k, (npw_k))
2532 127 : call mkkin(dtset%ecut, dtset%ecutsm, dtset%effmass_free, cryst%gmet, kg_k, kinpw_k, kpt, npw_k, 0, 0)
2533 :
2534 : ! Compute (k+G) vectors (only if useylm=1)
2535 127 : nkpg = 3 * dtset%nloalg(3)
2536 508 : ABI_MALLOC(kpg_k, (npw_k, nkpg))
2537 127 : if (nscf%paral_kgb0 /= 1 .and. nkpg > 0) call mkkpg(kg_k, kpg_k, kpt, nkpg, npw_k)
2538 :
2539 : ! Get one-dimensional structure factor information on the coarse grid.
2540 381 : ABI_MALLOC(ph1d, (2,3*(2*mgfft+1)*cryst%natom))
2541 127 : call getph(cryst%atindx, cryst%natom, n1, n2, n3, ph1d, cryst%xred)
2542 :
2543 : ! Initi GS Hamiltonian.
2544 : call gs_ham_k%init(psps, pawtab, nspinor, dtset%nsppol, dtset%nspden, cryst%natom, &
2545 : dtset%typat, cryst%xred, nfft, mgfft, nscf%ngfft, cryst%rprimd, dtset%nloalg, &
2546 : comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab, mpi_spintab=mpi_enreg%my_isppoltab, &
2547 127 : usecprj=dtset%usepaw, ph1d=ph1d, nucdipmom=dtset%nucdipmom, gpu_option=dtset%gpu_option)
2548 127 : ABI_FREE(ph1d)
2549 :
2550 : ! Set up local potential vlocal on the coarse FFT mesh from vtrial taking into account the spin.
2551 : ! Also, continue to initialize the Hamiltonian.
2552 127 : nvloc = gs_ham_k%nvloc
2553 4564761 : ABI_CALLOC(vlocal, (n4, n5, n6, nvloc))
2554 :
2555 : ! ncomp=Number of extra components in vtrial and vlocal (e.g. 1 if LDA/GGA pot, 4 for Meta-GGA, etc).
2556 127 : ncomp = 1
2557 : call gspot_transgrid_and_pack(isppol, psps%usepaw, nscf%paral_kgb0, nfft, nscf%ngfft, nfftf, &
2558 127 : dtset%nspden, gs_ham_k%nvloc, ncomp, pawfgr, mpi_enreg, nscf%vtrial, vlocal)
2559 :
2560 127 : call gs_ham_k%load_spin(isppol, vlocal=vlocal, with_nonlocal=.true.)
2561 :
2562 127 : if (dtset%usekden /= 0) then
2563 0 : ABI_ERROR("nscf_init with mgga not yet coded")
2564 : !call gspot_transgrid_and_pack(isppol, psps%usepaw, nscf%paral_kgb0, dtset%nfft, dtset%ngfft, nfftf, &
2565 : ! dtset%nspden, gs_ham_k%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal)
2566 : !call gs_ham_k%load_spin(isppol, vxctaulocal=vxctaulocal)
2567 : end if
2568 :
2569 : ! Set up the spherical harmonics (Ylm) at k.
2570 508 : ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2 * psps%useylm))
2571 127 : if (psps%useylm == 1) then
2572 84 : call initylmg_k(npw_k, psps%mpsang, optder0, cryst%rprimd, cryst%gprimd, kpt, kg_k, ylm_k, ylmgr_dum)
2573 : end if
2574 :
2575 : ! Compute nonlocal form factors ffnl_k at (k+G)
2576 635 : ABI_MALLOC(ffnl_k, (npw_k, 1, psps%lmnmax, psps%ntypat))
2577 127 : call mkffnl_objs(cryst, psps, 1, ffnl_k, ider0, idir0, kg_k, kpg_k, kpt, nkpg, npw_k, ylm_k, ylmgr_dum)
2578 127 : ABI_FREE(ylm_k)
2579 :
2580 : ! Load k-dependent part in the Hamiltonian datastructure
2581 : ! - Compute 3D phase factors
2582 : ! - Prepare various tabs in case of band-FFT parallelism
2583 : ! - Load k-dependent quantities in the Hamiltonian
2584 :
2585 508 : ABI_MALLOC(ph3d_k, (2, npw_k, gs_ham_k%matblk))
2586 : call gs_ham_k%load_k(kpt_k=kpt, istwf_k=istwf_k, npw_k=npw_k, &
2587 : kinpw_k=kinpw_k, kg_k=kg_k, kpg_k=kpg_k, ffnl_k=ffnl_k, ph3d_k=ph3d_k, &
2588 127 : compute_ph3d=(nscf%paral_kgb0/=1), compute_gbound=(nscf%paral_kgb0/=1))
2589 :
2590 : ! Allocate output buffers.
2591 508 : ABI_MALLOC(cg_k, (2, npw_k*nspinor, nband_k))
2592 635 : ABI_MALLOC(gsc_k, (2, npw_k*nspinor, nband_k*dtset%usepaw))
2593 : end associate
2594 :
2595 127 : end subroutine nscf_setup_kpt
2596 : !!***
2597 :
2598 : !!****f* m_cgwf/nscf_solve_kpt
2599 : !! NAME
2600 : !! nscf_solve_kpt
2601 : !!
2602 : !! FUNCTION
2603 : !! Solves the NSCF equation. Simplified wrapper around cgwf.
2604 : !!
2605 : !! INPUT
2606 : !! isppol=Spin index
2607 : !! kpt(3)=K-point
2608 : !! cryst=Crystalline structure
2609 : !! dtset<dataset_type>=All input variables for this dataset.
2610 : !! dtfil <type(datafiles_type)>=variables related to files
2611 : !! gs_ham_k <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
2612 : !! use_cg_k: True if input cg_k should be used to initialize the eigensolver.
2613 : !! npw_k=Number of planewaves
2614 : !!
2615 : !! OUTPUT
2616 : !! kg_k=g-vectors for this k-point
2617 : !! cg_k=Wavefunction block.
2618 : !! gsc_k=<g|S|c> for PAW
2619 : !! eig_k=Eigenvalues.
2620 : !! msg=Error message.
2621 : !! ierr=Exit status.
2622 : !!
2623 : !! SOURCE
2624 :
2625 127 : subroutine nscf_solve_kpt(nscf, isppol, kpt, istwf_k, nband_k, cryst, dtset, dtfil, gs_ham_k, use_cg_k, & ! in
2626 127 : npw_k, cg_k, gsc_k, eig_k, msg, ierr) ! out
2627 :
2628 : !Arguments ------------------------------------
2629 : class(nscf_t),intent(inout) :: nscf
2630 : logical,intent(in) :: use_cg_k
2631 : integer,intent(in) :: isppol, istwf_k, nband_k, npw_k
2632 : real(dp),intent(in) :: kpt(3)
2633 : type(dataset_type),intent(in) :: dtset
2634 : type(datafiles_type), intent(in) :: dtfil
2635 : type(crystal_t),intent(in) :: cryst
2636 : type(gs_hamiltonian_type),intent(inout) :: gs_ham_k
2637 : !arrays
2638 : real(dp),intent(inout) :: cg_k(2, npw_k * dtset%nspinor, nband_k)
2639 : real(dp),intent(inout) :: gsc_k(2, npw_k * dtset%nspinor, nband_k * dtset%usepaw)
2640 : real(dp),allocatable,intent(out) :: eig_k(:)
2641 : integer,intent(out) :: ierr
2642 : character(len=*),intent(out) :: msg
2643 :
2644 : !Local variables ------------------------------
2645 : !scalars
2646 : integer,parameter :: mcgq0 = 0, mkgq0 = 0, nkpt1 = 1, pwind_alloc0 = 0, use_subvnlx0 = 0, use_subovl0 = 0, ider0 = 0, idir0 = 0
2647 : integer,parameter :: icg0 = 0, igsc0 = 0, ikpt0 = 0, quit0 = 0, ortalgo_3 = 3, mkmem1 = 1, useylmgr0 = 0
2648 : integer :: mcg, mgsc, n1, n2, n3, n4, n5, n6, nfft, nfftf, mgfft, mgfftf, inonsc, npwsp, me_g0, linalg_max_size
2649 : integer :: nspinor, ii, iband
2650 : real(dp),parameter :: cpus0 = zero
2651 : real(dp) :: max_resid
2652 127 : type(efield_type) :: dtefield
2653 : !arrays
2654 : integer :: npwarr_k(1), pwind(pwind_alloc0,2,3)
2655 254 : real(dp) :: pwnsfac(2,pwind_alloc0), pwnsfacq(2,mkgq0), zshift(nband_k), cgq(2, mcgq0), dphase_k(3)
2656 254 : real(dp) :: subovl(nband_k*(nband_k+1)*use_subovl0), subvnlx(nband_k*(nband_k+1)*use_subvnlx0)
2657 127 : real(dp),allocatable :: subham(:), resid_k(:), evec(:,:)
2658 : ! *************************************************************************
2659 :
2660 : ! See vtorho.F90 for the sequence of calls required to initialize the GS Hamiltonian.
2661 : associate (mpi_enreg => nscf%mpi_enreg, kg_k => gs_ham_k%kg_k)
2662 :
2663 : ! FFT meshes from input file, not necessary equal to the ones found in the external files.
2664 : nfftf = product(nscf%ngfftf(1:3)); mgfftf = maxval(nscf%ngfftf(1:3))
2665 : nfft = product(nscf%ngfft(1:3)) ; mgfft = maxval(nscf%ngfft(1:3))
2666 127 : n1 = nscf%ngfft(1); n2 = nscf%ngfft(2); n3 = nscf%ngfft(3); n4 = nscf%ngfft(4); n5 = nscf%ngfft(5); n6 = nscf%ngfft(6)
2667 127 : nspinor = dtset%nspinor
2668 :
2669 254 : npwarr_k = npw_k; npwsp = npw_k * nspinor; me_g0 = 1; mcg = npw_k * nspinor * nband_k; mgsc = mcg * dtset%usepaw
2670 :
2671 381 : ABI_MALLOC(resid_k, (nband_k))
2672 254 : ABI_MALLOC(eig_k, (nband_k))
2673 :
2674 127 : if (.not. use_cg_k) then
2675 : ! Initialize the wavefunctions with random numbers.
2676 44 : call cg_randomize(istwf_k, npw_k, nspinor, nband_k, me_g0, cg_k)
2677 : ! Multiply with envelope function to reduce kinetic energy.
2678 44 : call cg_envlop(cg_k, dtset%ecut, cryst%gmet, icg0, kg_k, kpt, mcg, nband_k, npw_k, nspinor)
2679 : end if
2680 :
2681 : ! Ortoghonalize input trial states (this is important, even when cg_k is already initialized from a previous k-point.
2682 127 : call pw_orthon(icg0, igsc0, istwf_k, mcg, mgsc, npwsp, nband_k, ortalgo_3, gsc_k, dtset%usepaw, cg_k, me_g0, xmpi_comm_self)
2683 : !call cg_kfilter(npw_k, nspinor, nband_k, gs_ham_k%kinpw_k, cg_k)
2684 :
2685 : ! linalg initialisation (required by subdiago)
2686 9271 : linalg_max_size = maxval(dtset%nband(:))
2687 : call abi_linalg_init(linalg_max_size, RUNL_GSTATE, dtset%wfoptalg, nscf%paral_kgb0,&
2688 127 : dtset%gpu_option, dtset%use_slk, dtset%np_slk, nscf%mpi_enreg%comm_bandspinorfft)
2689 :
2690 381 : ABI_MALLOC(subham, (nband_k*(nband_k+1)))
2691 508 : ABI_MALLOC(evec, (2*nband_k, nband_k))
2692 :
2693 : ! NSCF iterations.
2694 127 : ierr = 1; msg = ""
2695 909 : do inonsc=1,dtset%nstep
2696 :
2697 : call cgwf(dtset%berryopt, cg_k, cgq, dtset%chkexit, cpus0, dphase_k, dtefield, dtfil%filnam_ds(1), &
2698 : gsc_k, gs_ham_k, icg0, igsc0, ikpt0, inonsc, isppol, nband_k, mcg, mcgq0, mgsc, mkgq0, &
2699 : mpi_enreg, npw_k, nband_k, dtset%nbdblock, nkpt1, dtset%nline, npw_k, npwarr_k, dtset%nspinor, &
2700 : dtset%nsppol, dtset%ortalg, dtset%prtvol, &
2701 : pwind, pwind_alloc0, pwnsfac, pwnsfacq, quit0, resid_k, &
2702 909 : subham, subovl, subvnlx, dtset%tolrde, dtset%tolwfr_diago, use_subovl0, use_subvnlx0, mod(dtset%wfoptalg, 100), zshift)
2703 :
2704 : ! subspace rotation (without this, cgwf will never converge!)
2705 : call subdiago(cg_k, eig_k, evec, gsc_k, icg0, igsc0, istwf_k, mcg, mgsc, nband_k, npw_k, dtset%nspinor, nscf%paral_kgb0, &
2706 909 : subham, subovl, use_subovl0, gs_ham_k%usepaw, me_g0)
2707 :
2708 : ! Check for convergence.
2709 909 : if (dtset%nbdbuf >= 0) then
2710 10908 : max_resid = maxval(resid_k(1:max(1,nband_k-dtset%nbdbuf)))
2711 : else
2712 0 : ABI_ERROR(sjoin('Bad value for nbdbuf:', itoa(dtset%nbdbuf)))
2713 : end if
2714 :
2715 : ! Exit loop over inonsc if converged
2716 909 : if (max_resid < dtset%tolwfr) then
2717 127 : ierr = 0
2718 127 : msg = sjoin(" NSCF for kpt:", ktoa(kpt), "spin:", itoa(isppol))
2719 127 : msg = sjoin(msg, ", completed in: ", itoa(inonsc), "steps. max_resid:", ftoa(max_resid))
2720 127 : call wrtout(std_out, msg)
2721 :
2722 : ! Print energies and residuals
2723 127 : if (dtset%prtvol > 10) then
2724 0 : do ii=0,(nband_k-1)/8
2725 0 : write(msg, '(a,8es10.2)' )' ene:',(eig_k(iband) * Ha_eV,iband=1+ii*8,min(nband_k,8+ii*8)); call wrtout(std_out, msg)
2726 0 : write(msg, '(a,8es10.2)' )' res:',(resid_k(iband), iband=1+ii*8,min(nband_k,8+ii*8)); call wrtout(std_out, msg)
2727 : end do
2728 : end if
2729 :
2730 : exit ! EXIT NSCF cycle
2731 : end if
2732 : end do ! inonsc (NON SELF-CONSISTENT LOOP)
2733 :
2734 127 : if (ierr /= 0) then
2735 0 : msg = sjoin(" NSCF run for kpt:", ktoa(kpt), "spin", itoa(isppol))
2736 0 : msg = sjoin(msg, " didn't converge after", itoa(dtset%nstep), " steps", ch10)
2737 0 : msg = sjoin(msg, "max_resid:", ftoa(max_resid), " >= tolwfr:", ftoa(dtset%tolwfr))
2738 : end if
2739 :
2740 : ! Fix the phase of the wavefunctions.
2741 127 : call cgtk_fixphase(cg_k, gsc_k, icg0, igsc0, istwf_k, mcg, mgsc, mpi_enreg, nband_k, npw_k, dtset%usepaw)
2742 :
2743 127 : call abi_linalg_finalize(dtset%gpu_option)
2744 :
2745 127 : ABI_FREE(subham)
2746 127 : ABI_FREE(evec)
2747 127 : ABI_FREE(resid_k)
2748 : end associate
2749 :
2750 254 : end subroutine nscf_solve_kpt
2751 : !!***
2752 :
2753 : !!****f* m_cgwf/nscf_free
2754 : !! NAME
2755 : !! nscf_free
2756 : !!
2757 : !! FUNCTION
2758 : !! Free dynamic memory
2759 : !!
2760 : !! SOURCE
2761 :
2762 2 : subroutine nscf_free(nscf)
2763 :
2764 : !Arguments ------------------------------------
2765 : class(nscf_t),intent(inout) :: nscf
2766 : !************************************************************************
2767 :
2768 2 : ABI_SFREE(nscf%vtrial)
2769 2 : call destroy_mpi_enreg(nscf%mpi_enreg)
2770 :
2771 2 : end subroutine nscf_free
2772 : !!***
2773 :
2774 4 : end module m_cgwf
2775 : !!***
|
