Branch data 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 : 444775 : subroutine cgwf(berryopt,cg,cgq,chkexit,cpus,dphase_k,dtefield,&
197 : 444775 : filnam_ds1,gsc,gs_hamk,icg,igsc,ikpt,inonsc,&
198 : : isppol,mband,mcg,mcgq,mgsc,mkgq,mpi_enreg,&
199 : 444775 : mpw,nband,nbdblock,nkpt,nline,npw,npwarr,&
200 : 444775 : nspinor,nsppol,ortalg,prtvol,pwind,&
201 : 444775 : pwind_alloc,pwnsfac,pwnsfacq,quit,resid,subham,subovl,&
202 : 444775 : 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 : 444775 : 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 : 889550 : real(dp),allocatable :: cg1_k(:,:),cgq_k(:,:),conjgr(:,:),cwavef(:,:)
243 : 444775 : real(dp),allocatable :: detovc(:,:,:),detovd(:,:,:),direc(:,:),direc_tmp(:,:)
244 : 444775 : real(dp),allocatable :: gh_direc(:,:),gh_direcws(:,:),ghc(:,:),ghc_all(:,:),ghcws(:,:)
245 : 444775 : real(dp),allocatable :: grad_berry(:,:),grad_total(:,:),gs_direc(:,:)
246 : 444775 : real(dp),allocatable :: gvnlxc(:,:),gvnlx_direc(:,:),gvnlx_dummy(:,:)
247 : 444775 : real(dp),allocatable :: pcon(:),pwnsfac_k(:,:),scprod(:,:),scwavef(:,:)
248 : 444775 : real(dp),allocatable :: smat_inv(:,:,:),smat_k(:,:,:),smat_k_paw(:,:,:),swork(:,:),vresid(:,:),work(:,:)
249 [ + + + + : 3113425 : type(pawcprj_type) :: cprj_dum(1,1)
+ + + - +
- ]
250 : 444775 : type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kb(:,:)
251 : 444775 : 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 : 444775 : call timab(22,1,tsec)
258 : :
259 : : !Touching chkexit, cpus,filnam_ds to avoid warning for abirules. This is dirty...
260 [ - + ]: 444775 : if(chkexit<0)then
261 : 0 : ABI_BUG('chkexit should be positive!')
262 : : end if
263 : :
264 [ - + - - ]: 444775 : 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 : 444775 : spaceComm_distrb=mpi_enreg%comm_cell
274 : 444775 : nproc_distrb=xmpi_comm_size(spaceComm_distrb)
275 : 444775 : me_distrb=mpi_enreg%me_kpt
276 : 444775 : me_g0 = mpi_enreg%me_g0
277 : 444775 : 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 : 444775 : gs_hamk%gpu_option=ABI_GPU_DISABLED
280 [ - + ]: 444775 : 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 : 444775 : gen_eigenpb=(gs_hamk%usepaw==1)
285 [ + + ]: 444775 : useoverlap=0;if (gen_eigenpb) useoverlap=1
286 : :
287 : : !Initializations and allocations
288 : 444775 : isubh=1;isubo=1
289 : 444775 : nblock=(nband-1)/nbdblock+1
290 : 444775 : istwf_k=gs_hamk%istwf_k
291 : 444775 : wfopta10=mod(wfoptalg,10)
292 [ + + ]: 444775 : optekin=0;if (wfoptalg>=10) optekin=1
293 : 444775 : natom=gs_hamk%natom
294 : 444775 : cpopt=-1
295 : 444775 : num_warning = 0
296 : :
297 [ - + + - : 1334325 : ABI_MALLOC(pcon,(npw))
+ - + - -
+ - + ]
298 [ - + + - : 1334325 : ABI_MALLOC(ghc,(2,npw*nspinor))
+ - + - -
+ - + -
+ ]
299 [ - + + - : 889550 : ABI_MALLOC(gvnlxc,(2,npw*nspinor))
- + - + -
+ ]
300 [ - + + - : 889550 : ABI_MALLOC(conjgr,(2,npw*nspinor))
- + - + ]
301 [ - + + - : 889550 : ABI_MALLOC(cwavef,(2,npw*nspinor))
- + - + -
+ ]
302 [ - + + - : 889550 : ABI_MALLOC(direc,(2,npw*nspinor))
- + - + -
+ ]
303 [ - + + - : 1334325 : ABI_MALLOC(scprod,(2,nband))
+ - + - -
+ - + ]
304 : :
305 [ - + + - : 889550 : ABI_MALLOC(gh_direc,(2,npw*nspinor))
- + - + -
+ ]
306 [ - + + - : 889550 : ABI_MALLOC(gvnlx_direc,(2,npw*nspinor))
- + - + -
+ ]
307 [ - + + - : 889550 : ABI_MALLOC(vresid,(2,npw*nspinor))
- + - + ]
308 : :
309 [ + + ]: 444775 : if (gen_eigenpb) then
310 [ - + + - : 195976 : ABI_MALLOC(gs_direc,(2,npw*nspinor))
- + - + -
+ ]
311 : : else
312 [ - + - + ]: 346787 : ABI_MALLOC(gs_direc,(0,0))
313 : : end if
314 : :
315 [ + + ]: 444775 : if (gen_eigenpb) then
316 [ - + + - : 195976 : ABI_MALLOC(scwavef,(2,npw*nspinor))
- + - + -
+ ]
317 [ - + + - : 195976 : ABI_MALLOC(direc_tmp,(2,npw*nspinor))
- + - + ]
318 : : end if
319 : :
320 [ + + ]: 97988 : if (gen_eigenpb.and.(inonsc==1)) then
321 [ - + + - : 220662 : ABI_MALLOC_OR_DIE(ghc_all,(2,nband*npw*nspinor), ierr)
+ - + - +
- - + -
+ ]
322 : : end if
323 : :
324 [ + + ]: 444775 : if (wfopta10==2.or.wfopta10==3) then
325 [ - + + - : 118 : ABI_MALLOC(work,(2,npw*nspinor))
- + - + -
+ ]
326 : : end if
327 : :
328 [ + + - + ]: 444775 : if (gen_eigenpb.and.(wfopta10==2.or.wfopta10==3)) then
329 [ # # # # : 0 : ABI_MALLOC(swork,(2,npw*nspinor))
# # # # #
# ]
330 : : else
331 [ - + - + ]: 444775 : ABI_MALLOC(swork,(0,0))
332 : : end if
333 : :
334 [ + + ]: 444775 : 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 [ - + - + ]: 444716 : ABI_MALLOC(gvnlx_dummy,(0,0))
340 : : end if
341 : :
342 : : !if "generalized eigenproblem", not sure of wfoptalg=2,3 algorithms
343 [ + + - + ]: 444775 : 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 [ + - - - : 444775 : & (berryopt == 14) .or. (berryopt == 16) .or. (berryopt == 17) )
+ - ]
366 : 444775 : 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 [ + + + + ]: 444775 : if (gen_eigenpb.and.inonsc==1) then
447 [ + + ]: 969600 : do iblock=1,nblock
448 : 896046 : ibandmin=1+(iblock-1)*nbdblock
449 : 896046 : ibandmax=min(iblock*nbdblock,nband)
450 [ + + ]: 1865646 : do iband=ibandmin,ibandmax
451 : 896046 : ibdblock=iband-(iblock-1)*nbdblock
452 : 896046 : icg_shift=npw*nspinor*(iband-1)+icg
453 : 896046 : igsc_shift=npw*nspinor*(iband-1)+igsc
454 : :
455 : 896046 : 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 : 896046 : call fock_set_ieigen(gs_hamk%fockcommon,iband)
460 : 896046 : sij_opt=1
461 [ + + + - ]: 896046 : 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 : 878510 : 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 : 896046 : call cg_zcopy(npw*nspinor,ghc,ghc_all(1,1+icg_shift-icg))
470 : 1792092 : 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 [ + + ]: 3627827 : do iblock=1,nblock
477 : 3183052 : 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 : 3183052 : ibandmin=1+(iblock-1)*nbdblock
482 : 3183052 : ibandmax=min(iblock*nbdblock,nband)
483 : :
484 : : ! Big iband loop
485 [ + + ]: 6366104 : do iband=ibandmin,ibandmax
486 : 3183052 : ibdblock=iband-(iblock-1)*nbdblock
487 : 3183052 : counter=100*iband+inonsc
488 : 3183052 : icg_shift=npw*nspinor*(iband-1)+icg
489 : 3183052 : igsc_shift=npw*nspinor*(iband-1)+igsc
490 : :
491 : : ! ======================================================================
492 : : ! ========== INITIALISATION OF MINIMIZATION ITERATIONS =================
493 : : ! ======================================================================
494 : :
495 [ + + ]: 3183052 : 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 : 3183052 : dotgp=one
501 [ + + ]: 3183052 : 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 : 3183052 : call cg_zcopy(npw*nspinor,cg(1,1+icg_shift),cwavef)
512 : :
513 : : ! If generalized eigenproblem: extraction of the overlap information
514 [ + + ]: 3183052 : if (gen_eigenpb) then
515 : 1254370 : 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 : 1254370 : call dotprod_g(dotr,doti,istwf_k,npw*nspinor,2,cwavef,scwavef,me_g0,mpi_enreg%comm_spinorfft)
523 : 1254370 : dotr=sqrt(dotr**2+doti**2); xnorm=one/sqrt(dotr)
524 [ + + ]: 3763110 : call cg_zscal(npw*nspinor,(/xnorm,zero/),cwavef)
525 [ + + ]: 3763110 : call cg_zscal(npw*nspinor,(/xnorm,zero/),scwavef)
526 : : else
527 : 1928682 : call sqnorm_g(dotr,istwf_k,npw*nspinor,cwavef,me_g0,mpi_enreg%comm_fft)
528 : 1928682 : xnorm=one/sqrt(abs(dotr))
529 [ + + ]: 5786046 : call cg_zscal(npw*nspinor,(/xnorm,zero/),cwavef)
530 : : end if
531 : :
532 [ - + ]: 3183052 : 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 [ + + + + ]: 3183052 : if (gen_eigenpb.and.(inonsc==1)) then
539 : :
540 : : !$OMP PARALLEL DO PRIVATE(ipw)
541 [ + + ]: 326743741 : do ipw=1,npw*nspinor
542 : 325847695 : ghc(1,ipw)=xnorm*ghc_all(1,ipw+icg_shift-icg)
543 : 326743741 : 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 : 2287006 : call fock_set_ieigen(gs_hamk%fockcommon,iband)
549 : 2287006 : sij_opt=0
550 : 2287006 : 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 [ + + ]: 3183052 : 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 [ + - ]: 3183052 : if(nline/=0)then
575 [ + + ]: 12175728 : do iline=1,nline
576 : :
577 : : ! === COMPUTE THE RESIDUAL ===
578 : :
579 : : ! Compute lambda = <C|H|C> or <C|(H-zshift)**2|C>
580 : 10669929 : call dotprod_g(chc,doti,istwf_k,npw*nspinor,1,cwavef,ghc,me_g0,mpi_enreg%comm_spinorfft)
581 : 10669929 : lam0=chc
582 : :
583 : : ! Check that lam0 is decreasing on succeeding lines:
584 [ + + ]: 10669929 : if (.not.finite_field) then
585 [ + + ]: 10494849 : if (iline==1) then
586 : 3135532 : lamold=lam0
587 : : else
588 [ - + - - ]: 7359317 : 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 : 7359317 : 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 [ + + ]: 10669929 : if (wfopta10<=1) then
603 : 10667708 : eval=chc
604 [ + + ]: 10667708 : if (gen_eigenpb) then
605 : :
606 : : !$OMP PARALLEL DO
607 [ + + ]: 1545482720 : do ipw=1,npw*nspinor
608 : 1541145613 : vresid(1,ipw)=ghc(1,ipw)-chc*scwavef(1,ipw)
609 : 1545482720 : vresid(2,ipw)=ghc(2,ipw)-chc*scwavef(2,ipw)
610 : : end do
611 : : else
612 : : !$OMP PARALLEL DO
613 [ + + ]: 2189359997 : do ipw=1,npw*nspinor
614 : 2183029396 : vresid(1,ipw)=ghc(1,ipw)-chc*cwavef(1,ipw)
615 : 2189359997 : 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 : 10669929 : call sqnorm_g(resid(iband),istwf_k,npw*nspinor,vresid,me_g0,mpi_enreg%comm_fft)
637 : :
638 [ - + ]: 10669929 : 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 [ + + ]: 10669929 : if (resid(iband)<tolwfr) then
649 [ + + ]: 327037 : 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 : 327037 : nskip=nskip+(nline-iline+1) ! Number of two-way 3D ffts skipped
655 : 327037 : exit ! Exit from the loop on iline
656 : : end if
657 : :
658 : : ! If user require exiting the job, stop line minimisations
659 [ - + ]: 10342892 : 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 [ + + ]: 10342892 : if (gen_eigenpb) then
673 : 4311687 : call cg_zcopy(npw*nspinor,vresid,direc) ! Store <G|H-lambda.S|C> in direc
674 : : else
675 : 6031205 : 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 [ + + ]: 10342892 : 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 [ + + ]: 10342892 : if(ortalg>=0)then
720 [ + + ]: 6032285 : 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 : 6030643 : 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 [ + + ]: 4310607 : 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 [ + + + - : 4622176022 : 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 [ + + + + ]: 10342892 : if (wfoptalg<10.or.iline==1) then
748 [ + + ]: 7276245 : if (gs_hamk%use_gbt == 0) then
749 : 7115439 : 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 : 160806 : 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 : 160806 : direc(:,npw+1:),mpi_enreg%comm_fft)
754 : : end if
755 : :
756 [ + + ]: 7276245 : 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 [ + + ]: 6576138 : do ispinor=1,nspinor
770 : 3509491 : igs=(ispinor-1)*npw
771 : : !$OMP PARALLEL DO
772 [ + + ]: 1115647718 : do ipw=1+igs,npw+igs
773 : 1109071580 : direc(1,ipw)=direc(1,ipw)*pcon(ipw-igs)
774 : 1112581071 : 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 : 10342892 : gsc,scprod,0,tim_projbd,useoverlap,me_g0,mpi_enreg%comm_fft)
784 : :
785 : : ! ======================================================================
786 : : ! ================= COMPUTE THE CONJUGATE-GRADIENT =====================
787 : : ! ======================================================================
788 : :
789 [ + + ]: 10342892 : 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 [ + + ]: 10167812 : if (gen_eigenpb) then
800 : 4230833 : call dotprod_g(dotgg,doti,istwf_k,npw*nspinor,1,direc,direc_tmp,me_g0,mpi_enreg%comm_spinorfft)
801 : : else
802 : 5936979 : 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 [ - + ]: 10342892 : 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 [ + + ]: 10342892 : if (iline==1) then
811 : 2930489 : gamma=zero
812 : 2930489 : dotgp=dotgg
813 : 2930489 : call cg_zcopy(npw*nspinor,direc,conjgr)
814 [ - + ]: 2930489 : 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 : 7412403 : gamma=dotgg/dotgp
821 : 7412403 : dotgp=dotgg
822 : :
823 [ - + ]: 7412403 : 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 [ + + ]: 2626956370 : do ipw=1,npw*nspinor
837 : 2619543967 : conjgr(1,ipw)=direc(1,ipw)+gamma*conjgr(1,ipw)
838 : 2626956370 : 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 [ + + ]: 10342892 : if (gen_eigenpb) then
848 : 4311687 : call dotprod_g(dotr,doti,istwf_k,npw*nspinor,3,scwavef,conjgr,me_g0,mpi_enreg%comm_spinorfft)
849 : : else
850 : 6031205 : 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 [ + + ]: 10342892 : if(istwf_k==1)then
857 : :
858 : : !$OMP PARALLEL DO
859 [ + + ]: 3102362559 : do ipw=1,npw*nspinor
860 : 3092820576 : direc(1,ipw)=conjgr(1,ipw)-(dotr*cwavef(1,ipw)-doti*cwavef(2,ipw))
861 : 3102362559 : direc(2,ipw)=conjgr(2,ipw)-(dotr*cwavef(2,ipw)+doti*cwavef(1,ipw))
862 : : end do
863 : : else
864 : : !$OMP PARALLEL DO
865 [ + + ]: 519907084 : do ipw=1,npw*nspinor
866 : 519106175 : direc(1,ipw)=conjgr(1,ipw)-dotr*cwavef(1,ipw)
867 : 519907084 : 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 [ + + ]: 10342892 : if (.not.gen_eigenpb) then
874 : 6031205 : call sqnorm_g(dotr,istwf_k,npw*nspinor,direc,me_g0,mpi_enreg%comm_fft)
875 : 6031205 : xnorm=one/sqrt(abs(dotr))
876 [ + + ]: 18093615 : call cg_zscal(npw*nspinor,(/xnorm,zero/),direc)
877 : 6031205 : 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 : 10342892 : sij_opt=0;if (gen_eigenpb) sij_opt=1
886 : :
887 : 10342892 : 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 [ + + ]: 10342892 : 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 [ + + ]: 10342892 : if (gen_eigenpb) then
911 : 4311687 : call dotprod_g(dotr,doti,istwf_k,npw*nspinor,1,direc,gs_direc,me_g0,mpi_enreg%comm_spinorfft)
912 : 4311687 : xnorm=one/sqrt(abs(dotr))
913 : : end if
914 : :
915 : : ! Compute dhc = Re{<D|H|C>}
916 : 10342892 : call dotprod_g(dhc,doti,istwf_k,npw*nspinor,1,direc,ghc,me_g0,mpi_enreg%comm_spinorfft)
917 : 10342892 : dhc=dhc*xnorm
918 : :
919 : : ! Compute <D|H|D> or <D|(H-zshift)^2|D>
920 : 10342892 : call dotprod_g(dhd,doti,istwf_k,npw*nspinor,1,direc,gh_direc,me_g0,mpi_enreg%comm_spinorfft)
921 : 10342892 : dhd=dhd*xnorm**2
922 : :
923 [ - + ]: 10342892 : 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 [ + + ]: 10342892 : if (.not.finite_field) then
933 : : ! Compute tan(2 theta),sin(theta) and cos(theta)
934 : 10167812 : tan2th=2.0_dp*dhc/(chc-dhd)
935 : :
936 [ + + ]: 10167812 : if (abs(tan2th)<1.d-05) then
937 : 3616091 : costh=1.0_dp-0.125_dp*tan2th**2
938 : 3616091 : sinth=0.5_dp*tan2th*(1.0_dp-0.375_dp*tan2th**2)
939 : :
940 : : ! Check that result is above machine precision
941 [ + + ]: 3616091 : 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 : 6551721 : root=sqrt(1.0_dp+tan2th**2)
953 : 6551721 : costh=sqrt(0.5_dp+0.5_dp/root)
954 : 6551721 : 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 : 10163408 : diff=(chc-dhd)
959 : : ! Swap c and d if value of diff is positive
960 [ + + ]: 10163408 : if (diff>zero) then
961 : 31537 : swap=costh
962 : 31537 : costh=-sinth
963 : 31537 : sinth=swap
964 [ + + ]: 31537 : 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 : 10338488 : sintn=sinth*xnorm
1093 : :
1094 : : !$OMP PARALLEL DO
1095 [ + + ]: 3621136270 : do ipw=1,npw*nspinor
1096 : 3610797782 : cwavef(1,ipw)=cwavef(1,ipw)*costh+direc(1,ipw)*sintn
1097 : 3621136270 : 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 : 10338488 : call cg_zcopy(npw*nspinor,cwavef,cg(1,1+icg_shift))
1102 : :
1103 [ + + ]: 3621136270 : do ipw=1,npw*nspinor
1104 : 3610797782 : ghc(1,ipw) =ghc(1,ipw)*costh + gh_direc(1,ipw)*sintn
1105 : 3621136270 : ghc(2,ipw) =ghc(2,ipw)*costh + gh_direc(2,ipw)*sintn
1106 : : end do
1107 : :
1108 : :
1109 [ + + ]: 10338488 : if (use_subvnlx==1) then
1110 : : !$OMP PARALLEL DO
1111 [ + + ]: 2024909055 : do ipw=1,npw*nspinor
1112 : 2018917406 : gvnlxc(1,ipw)=gvnlxc(1,ipw)*costh + gvnlx_direc(1,ipw)*sintn
1113 : 2024909055 : 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 [ + + ]: 10338488 : if (gen_eigenpb) then
1119 : : !$OMP PARALLEL DO
1120 [ + + ]: 1539061522 : do ipw=1,npw*nspinor
1121 : 1534754033 : scwavef(1,ipw)=scwavef(1,ipw)*costh+gs_direc(1,ipw)*sintn
1122 : 1539061522 : 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 : 4307489 : call cg_zcopy(npw*nspinor,scwavef,gsc(1,1+igsc_shift))
1128 : :
1129 [ + + ]: 4307489 : 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 [ + + ]: 10338488 : 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 [ + + ]: 10338488 : if (.not.finite_field) then
1157 : 10163408 : 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 [ + + ]: 11844287 : if (iline==1) then
1183 : 2930121 : deold=deltae
1184 [ + + + + : 7408367 : else if (abs(deltae)<tolrde*abs(deold) .and. iline/=nline .and. wfopta10<2)then
+ + ]
1185 [ + + ]: 1345812 : 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 : 1345812 : nskip=nskip+2*(nline-iline) ! Number of one-way 3D ffts skipped
1191 : 1345812 : 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 [ + + ]: 3183052 : 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 [ + + + + : 3274852 : if(wfopta10==2 .or. wfopta10==3) ghc(:,:)=ghcws(:,:)
+ + ]
1231 : :
1232 [ + + + + : 7398412 : 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 : 3183052 : 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 : 3627827 : nspinor,subham,subovl,subvnlx,use_subovl,use_subvnlx,me_g0)
1248 : :
1249 : : end do ! iblock End loop over block of bands
1250 : :
1251 [ + + ]: 444775 : ABI_SFREE(dimlmn_srt)
1252 : :
1253 [ + + + + ]: 444775 : 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 [ - + ]: 444775 : 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 [ - + ]: 444775 : ABI_FREE(conjgr)
1294 [ - + ]: 444775 : ABI_FREE(cwavef)
1295 [ - + ]: 444775 : ABI_FREE(direc)
1296 [ - + ]: 444775 : ABI_FREE(pcon)
1297 [ - + ]: 444775 : ABI_FREE(scprod)
1298 [ - + ]: 444775 : ABI_FREE(ghc)
1299 [ - + ]: 444775 : ABI_FREE(gvnlxc)
1300 [ - + ]: 444775 : ABI_FREE(gh_direc)
1301 [ - + ]: 444775 : ABI_FREE(gvnlx_direc)
1302 [ - + ]: 444775 : ABI_FREE(vresid)
1303 [ - + ]: 444775 : ABI_FREE(gs_direc)
1304 [ - + ]: 444775 : ABI_FREE(gvnlx_dummy)
1305 [ - + ]: 444775 : ABI_FREE(swork)
1306 : :
1307 [ + + ]: 444775 : ABI_SFREE(scwavef)
1308 [ + + ]: 444775 : ABI_SFREE(direc_tmp)
1309 [ + + ]: 444775 : ABI_SFREE(ghc_all)
1310 [ + + ]: 444775 : ABI_SFREE(ghcws)
1311 [ + + ]: 444775 : ABI_SFREE(gh_direcws)
1312 [ + + ]: 444775 : ABI_SFREE(work)
1313 : :
1314 [ + + ]: 444775 : 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 : 444775 : gs_hamk%gpu_option=gpu_option_save
1350 [ - + ]: 444775 : 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 : 444775 : write(unit=msg,fmt=*) resid(1)
1354 : :
1355 : 444775 : call timab(22,2,tsec)
1356 : :
1357 : : DBG_EXIT("COLL")
1358 : :
1359 [ - + - + : 889550 : 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 : 3183052 : subroutine mksubham(cg,ghc,gsc,gvnlxc,iblock,icg,igsc,istwf_k,&
1829 : : & isubh,isubo,mcg,mgsc,nband_k,nbdblock,npw_k,&
1830 : 3183052 : & 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 [ + + ]: 6366104 : do iband=1+(iblock-1)*nbdblock,min(iblock*nbdblock,nband_k)
1855 : 3183052 : ibdblock=iband-(iblock-1)*nbdblock
1856 : :
1857 : : ! Compute elements of subspace Hamiltonian <C(i)|H|C(n)> and <C(i)|Vnl|C(n)>
1858 [ + + ]: 3183052 : if(istwf_k==1)then
1859 : :
1860 [ + + ]: 30215585 : do ii=1,iband
1861 : 27281659 : iwavef=(ii-1)*npw_k*nspinor+icg
1862 : 27281659 : chcre=zero ; chcim=zero
1863 [ + + ]: 27281659 : 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 [ + + ]: 6074015885 : do ipw=1,npw_k*nspinor
1874 : 6058610113 : cgreipw=cg(1,ipw+iwavef)
1875 : 6058610113 : cgimipw=cg(2,ipw+iwavef)
1876 : 6058610113 : chcre=chcre+cgreipw*ghc(1,ipw)+cgimipw*ghc(2,ipw)
1877 : 6074015885 : chcim=chcim+cgreipw*ghc(2,ipw)-cgimipw*ghc(1,ipw)
1878 : : end do
1879 : : cvcre=zero ; cvcim=zero
1880 [ + + ]: 6074015885 : do ipw=1,npw_k*nspinor
1881 : 6058610113 : cgreipw=cg(1,ipw+iwavef)
1882 : 6058610113 : cgimipw=cg(2,ipw+iwavef)
1883 : 6058610113 : cvcre=cvcre+cgreipw*gvnlxc(1,ipw)+cgimipw*gvnlxc(2,ipw)
1884 : 6074015885 : cvcim=cvcim+cgreipw*gvnlxc(2,ipw)-cgimipw*gvnlxc(1,ipw)
1885 : : end do
1886 : 15405772 : subvnlx(isubh )=cvcre
1887 : 15405772 : 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 : 27281659 : subham(isubh )=chcre
1900 : 27281659 : subham(isubh+1)=chcim
1901 : : ! subham(isubh )=chc(1)
1902 : : ! subham(isubh+1)=chc(2)
1903 : 30215585 : isubh=isubh+2
1904 : : end do
1905 : :
1906 [ + - ]: 249126 : else if(istwf_k>=2)then
1907 [ + + ]: 3206926 : do ii=1,iband
1908 : 2957800 : iwavef=(ii-1)*npw_k+icg
1909 : : ! Use the time-reversal symmetry, but should not double-count G=0
1910 [ + + + - ]: 2957800 : if(istwf_k==2 .and. me_g0==1) then
1911 : 1071602 : chcre = half*cg(1,1+iwavef)*ghc(1,1)
1912 [ + + ]: 1071602 : if (use_subvnlx==1) cvcre=half*cg(1,1+iwavef)*gvnlxc(1,1)
1913 : : ipw1=2
1914 : : else
1915 : 1886198 : chcre=zero; ipw1=1
1916 [ + + ]: 1886198 : if (use_subvnlx==1) cvcre=zero
1917 : : end if
1918 [ + + ]: 2957800 : if (use_subvnlx==0) then
1919 [ + + ]: 2419290 : do isp=1,nspinor
1920 [ + + ]: 1032233141 : do ipw=ipw1+(isp-1)*npw_k,npw_k*isp
1921 : 1029813851 : cgreipw=cg(1,ipw+iwavef)
1922 : 1029813851 : cgimipw=cg(2,ipw+iwavef)
1923 : 1031023496 : chcre=chcre+cgreipw*ghc(1,ipw)+cgimipw*ghc(2,ipw)
1924 : : end do
1925 : : end do
1926 : 1209645 : chcre=two*chcre
1927 : : else
1928 [ + + ]: 3496310 : do isp=1,nspinor
1929 [ + + ]: 984496585 : do ipw=ipw1+(isp-1)*npw_k,npw_k*isp
1930 : 981000275 : cgreipw=cg(1,ipw+iwavef)
1931 : 981000275 : cgimipw=cg(2,ipw+iwavef)
1932 : 981000275 : chcre=chcre+cgreipw*ghc(1,ipw)+cgimipw*ghc(2,ipw)
1933 : 982748430 : cvcre=cvcre+cgreipw*gvnlxc(1,ipw)+cgimipw*gvnlxc(2,ipw)
1934 : : end do
1935 : : end do
1936 : 1748155 : chcre=two*chcre
1937 : 1748155 : cvcre=two*cvcre
1938 : : ! Store real and imag parts in Hermitian storage mode:
1939 : 1748155 : subvnlx(isubh )=cvcre
1940 : 1748155 : subvnlx(isubh+1)=zero
1941 : : end if
1942 : 2957800 : subham(isubh )=chcre
1943 : 2957800 : subham(isubh+1)=zero
1944 : 3206926 : 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 [ - + ]: 6366104 : 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 : 3183052 : 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 : : !!***
|
