Line data Source code
1 : !!****m* ABINIT/m_rf2
2 : !! NAME
3 : !! m_rf2
4 : !!
5 : !! FUNCTION
6 : !! This module defines structures and provides procedures used to compute the 2nd order Sternheimer
7 : !! equation.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2015-2026 ABINIT group (LB,MT)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_rf2
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_errors
28 : use m_hamiltonian
29 : use m_cgtools
30 :
31 : use defs_abitypes, only : MPI_type
32 : use m_getgh1c, only : getgh1c
33 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_free, pawcprj_output
34 : use m_getghc, only : getghc
35 : use m_nonlop, only : nonlop
36 : use m_getgh2c, only : getgh2c
37 :
38 : implicit none
39 :
40 : private
41 : !!***
42 :
43 : !----------------------------------------------------------------------
44 :
45 : !!****t* m_rf2/rf2_t
46 : !! NAME
47 : !! rf2_t
48 : !!
49 : !! FUNCTION
50 : !! Datatype gathering all needed data for the computation of the 2nd order Sternheimer equation.
51 : !!
52 : !! SOURCE
53 :
54 : type,public :: rf2_t
55 :
56 : !scalars
57 : integer :: ndir ! number of directions to consider
58 : integer :: nband_k ! number of bands
59 : integer :: size_wf ! number of components in a wavefunction
60 : integer :: size_cprj ! number of components of a cprj variable (i.e. <p_i|WF>)
61 :
62 : !arrays
63 : integer :: iperts(2) ! perturbations to compute
64 : ! ipert = natom + 10 :
65 : ! iperts(1) = iperts(2) = natom+1 (wavevector k)
66 : !
67 : ! ipert = natom + 11 :
68 : ! iperts(1) = natom+1 (wavevector k)
69 : ! iperts(2) = natom+2 (electric field)
70 :
71 : integer :: idirs(2) ! directions of the perturbations (ndir=1 : idirs(1)=idirs(2) , ndir=2 : idirs(1)/=idirs(2))
72 :
73 : real(dp), ABI_CONTIGUOUS pointer :: RHS_Stern(:,:) => null() ! need a pointer because is a ptr target
74 : ! Right-hand side of the 2nd order Sternheimer equation, for every bands.
75 : ! Namely, for a band "n" :
76 : ! |(RHS_Stern)_n> = (H^(2)-epsilon_n S^(2)) |u^(0)_n> + 2(H^(1)-epsilon_n S^(1))|u^(1)_n>
77 : ! - 2 sum_m ( lambda_mn^(1) ( S^(1) |u^(0)_m> + S^(0) |u^(1)_m> )
78 : ! ( /!\ : in the sum, the index m runs over occupied bands )
79 : ! where :
80 : ! - epsilon_n = eigenvalue of the GS Hamiltonian
81 : ! - lambda_mn^(1) = <u^(0)_m| (H^(1)-(epsilon_n+epsilon_m)/2 S^(1) |u^(0)_n> (1st order Lagrange multiplier)
82 : ! - X^(2) = d^2X/(dk_dir1 dlambda_dir2)
83 : ! - 2X^(1)Y^(1) = dX/dk_dir1 dY/dlambda_dir2 + dX/dlambda_dir2 dY/dk_dir1
84 : ! lambda_dir2 = k_dir2 (ipert==natom+10) , E_dir2 (ipert==natom+11)
85 : ! Note : P_c^* will be apply to |(RHS_Stern)_n> in dfpt_cgwf.F90
86 : ! **
87 : ! Computed in "rf2_init"
88 :
89 : real(dp),allocatable :: amn(:,:)
90 : ! Scalar needed for the "orhtonormalization condition", see "dcwavef(:,:)" below
91 : ! Namely :
92 : ! A_mn = <u^(0)_m| S^(2) |u^(0)_n> + 2 <u^(1)_m| S^(0) |u^(1)_n>
93 : ! + 2 <u^(1)_m| S^(1) |u^(0)_n> + 2 <u^(0)_m| S^(1) |u^(1)_n>
94 : ! **
95 : ! Computed in "rf2_init", stored only for testing purposes
96 :
97 : real(dp),allocatable :: dcwavef(:,:)
98 : ! Vector needed to enforce the "orhtonormalization condition" on 2nd order wave functions.
99 : ! Namely :
100 : ! |dcwavef_n> = -1/2 sum_m A_mn |u^(0)_m>
101 : ! **
102 : ! Computed in "rf2_init"
103 :
104 : real(dp),allocatable :: lambda_mn(:,:)
105 : ! 2nd order Lagrange multiplier.
106 : ! Namely :
107 : ! lambda_mn = <u^(0)_m| H^(2) |u^(0)_n> + 2 <u^(1)_m| H^(0) |u^(1)_n>
108 : ! + 2 <u^(1)_m| H^(1) |u^(0)_n> + 2 <u^(0)_m| H^(1) |u^(1)_n>
109 : ! - A_mn * (epsilon_m + epsilon_n) / 2
110 : ! **
111 : ! Computed in "rf2_init"
112 :
113 : end type rf2_t
114 :
115 : public :: rf2_getidir
116 : public :: rf2_getidirs
117 : public :: rf2_accumulate_bands
118 : public :: rf2_apply_hamiltonian
119 : public :: rf2_destroy
120 :
121 : !!***
122 :
123 : !----------------------------------------------------------------------
124 :
125 : CONTAINS !==============================================================================
126 : !!***
127 :
128 : !----------------------------------------------------------------------
129 :
130 : !!****f* m_rf2/rf2_getidir
131 : !! NAME
132 : !! rf2_getidir
133 : !!
134 : !! FUNCTION
135 : !! Get the direction of the 2nd order perturbation according to inputs idir1 and idir2
136 : !!
137 : !! INPUTS
138 : !! idir1 : index of the 1st direction (1<=idir1<=3)
139 : !! idir2 : index of the 2nd direction (1<=idir2<=3)
140 : !!
141 : !! OUTPUT
142 : !! idir : integer between 1 and 9
143 : !!
144 : !! SOURCE
145 :
146 8584 : subroutine rf2_getidir(idir1,idir2,idir)
147 :
148 : !Arguments ---------------------------------------------
149 : !scalars
150 : integer,intent(in) :: idir1,idir2
151 : integer,intent(out) :: idir
152 : integer :: dir_from_dirs(9) = (/1,6,5,9,2,4,8,7,3/)
153 : ! *************************************************************************
154 :
155 8584 : idir = dir_from_dirs(3*(idir1-1)+idir2)
156 :
157 8584 : end subroutine rf2_getidir
158 : !!***
159 :
160 : !----------------------------------------------------------------------
161 :
162 : !!****f* m_rf2/rf2_getidirs
163 : !! NAME
164 : !! rf2_getidirs
165 : !!
166 : !! FUNCTION
167 : !! Get directions of the 1st and 2nd perturbations according to the input idir
168 : !!
169 : !! INPUTS
170 : !! idir: integer between 1 and 9
171 : !!
172 : !! OUTPUT
173 : !! idir1: index of the 1st direction (1<=idir1<=3)
174 : !! idir2: index of the 2nd direction (1<=idir2<=3)
175 : !!
176 : !! SOURCE
177 :
178 231905 : subroutine rf2_getidirs(idir,idir1,idir2)
179 :
180 : !Arguments ---------------------------------------------
181 : !scalars
182 : integer,intent(in) :: idir
183 : integer,intent(out) :: idir1,idir2
184 : ! *************************************************************************
185 :
186 262190 : select case(idir)
187 : ! Diagonal terms :
188 : case (1)
189 30285 : idir1 = 1
190 30285 : idir2 = 1
191 : case (2)
192 30285 : idir1 = 2
193 30285 : idir2 = 2
194 : case (3)
195 20529 : idir1 = 3
196 20529 : idir2 = 3
197 : ! Upper triangular terms :
198 : case (4)
199 31829 : idir1 = 2
200 31829 : idir2 = 3
201 : case (5)
202 21545 : idir1 = 1
203 21545 : idir2 = 3
204 : case (6)
205 31829 : idir1 = 1
206 31829 : idir2 = 2
207 : ! Lower triangular terms :
208 : case (7)
209 22302 : idir1 = 3
210 22302 : idir2 = 2
211 : case (8)
212 20999 : idir1 = 3
213 20999 : idir2 = 1
214 : case (9)
215 22302 : idir1 = 2
216 231905 : idir2 = 1
217 : end select
218 :
219 231905 : end subroutine rf2_getidirs
220 : !!***
221 :
222 : !----------------------------------------------------------------------
223 :
224 : !!****f* m_rf2/rf2_accumulate_bands
225 : !! NAME
226 : !! rf2_init
227 : !!
228 : !! FUNCTION
229 : !! Compute the scalar product < vi | v1j > and add it to rf2%lambda_mn
230 : !! If necessary, compute < vi | v2j > and add it to rf2%amn.
231 : !!
232 : !! INPUTS
233 : !! rf2 : rf2_t object containing all rf2 data
234 : !! choice : 4 possible values, see "select case (choice)" below
235 : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
236 : !! mpi_enreg=information about MPI parallelization
237 : !! iband : band index for vi
238 : !! idir1 (used only if debug_mode/=0) : direction of the 1st perturbation
239 : !! idir2 (used only if debug_mode/=0) : direction of the 2nd perturbation
240 : !! ipert1 (used only if debug_mode/=0) : 1st perturbation
241 : !! ipert2 (used only if debug_mode/=0) : 2nd perturbation
242 : !! jband : band index for v1j and v2j
243 : !! debug_mode : if /=0 : all < vi | v1j > and < vi | v2j > are printed in std_out
244 : !! vi,v1j,v2j : input vectors
245 : !!
246 : !! OUTPUT
247 : !!
248 : !! SOURCE
249 :
250 1480320 : subroutine rf2_accumulate_bands(rf2,choice,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,&
251 1480320 : jband,debug_mode,vi,v1j,v2j,factor_in)
252 :
253 : !Arguments ---------------------------------------------
254 : !scalars
255 : class(rf2_t),intent(inout) :: rf2
256 : integer,intent(in) :: choice,iband,idir1,idir2,ipert1,ipert2,jband,debug_mode
257 : type(gs_hamiltonian_type),intent(in) :: gs_hamkq
258 : type(MPI_type),intent(in) :: mpi_enreg
259 : real(dp), intent(in) :: factor_in
260 : !arrays
261 : real(dp),intent(in) :: vi(2,rf2%size_wf),v1j(2,rf2%size_wf),v2j(2,rf2%size_wf)
262 :
263 : !Local variables ---------------------------------------
264 : !scalars
265 : integer :: nband_k,size_wf
266 : real(dp) :: dotr,dot2r,doti,dot2i
267 : character(len=500) :: msg
268 : character(len=15) :: bra_i,ket_j,op1,op2
269 : character(len=2) :: pert1,pert2
270 : ! *************************************************************************
271 :
272 : DBG_ENTER("COLL")
273 :
274 1480320 : nband_k = rf2%nband_k
275 1480320 : size_wf = rf2%size_wf
276 :
277 1480320 : call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,vi,v1j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
278 :
279 1480320 : if(debug_mode/=0) then
280 24576 : if (ipert1 == gs_hamkq%natom+1) then
281 17664 : pert1 = "dk"
282 : else
283 6912 : pert1 = "dE"
284 : end if
285 24576 : if (ipert2 == gs_hamkq%natom+1) then
286 15360 : pert2 = "dk"
287 : else
288 9216 : pert2 = "dE"
289 : end if
290 31488 : select case (choice)
291 : case (1) ! < u^(pert2) | H^(0) | u^(pert1) >
292 6912 : write(bra_i,'(2a,2(i1,a))') ' < du/',pert2,idir2,'(',iband,') | '
293 6912 : write(ket_j,'(2a,2(i1,a))') ' | du/',pert1,idir1,'(',jband,') > '
294 6912 : write(op1,'(a)') ' H^(0) '
295 6912 : write(op2,'(a)') ' S^(0) '
296 : case (2) ! < u^(pert2) | H^(pert1) | u^(0) >
297 6912 : write(bra_i,'(2a,2(i1,a))') ' < du/',pert2,idir2,'(',iband,') | '
298 6912 : write(ket_j,'(a,i1,a)') ' | u^(0) (',jband,') > '
299 6912 : write(op1,'(2a,i1,a)') ' dH/',pert1,idir1,' '
300 6912 : write(op2,'(2a,i1,a)') ' dS/',pert1,idir1,' '
301 : case (3) ! < u^(0) | H^(pert1pert2) | u^(0) >
302 3840 : write(bra_i,'(a,i1,a)') ' < u^(0) (',iband,') | '
303 3840 : write(ket_j,'(a,i1,a)') ' | u^(0) (',jband,') > '
304 3840 : write(op1,'(2(2a,i1),a)') 'd^2H/(',pert1,idir1,' ',pert2,idir2,')'
305 3840 : write(op2,'(2(2a,i1),a)') 'd^2S/(',pert1,idir1,' ',pert2,idir2,')'
306 : case (4) ! < u^(0) | H^(pert2) | u^(pert1) >
307 6912 : write(bra_i,'(a,i1,a)') ' < u^(0) (',iband,') | '
308 6912 : write(ket_j,'(2a,2(i1,a))') ' | du/',pert1,idir1,'(',jband,') > '
309 6912 : write(op1,'(2a,i1,a)') ' dH/',pert2,idir2,' '
310 31488 : write(op2,'(2a,i1,a)') ' dS/',pert2,idir2,' '
311 : end select
312 24576 : dot2r = dotr ; if (abs(dot2r)<tol9) dot2r = zero ! in order to hide the numerical noise
313 24576 : dot2i = doti ; if (abs(dot2i)<tol9) dot2i = zero ! in order to hide the numerical noise
314 24576 : write(msg,'(3a,2(a,es17.8E3))') bra_i,op1,ket_j,' = ',dot2r,',',dot2i
315 24576 : call wrtout(std_out,msg)
316 : end if
317 :
318 4440960 : rf2%lambda_mn(:,iband+(jband-1)*nband_k) = factor_in*(/dotr,doti/) + rf2%lambda_mn(:,iband+(jband-1)*nband_k)
319 :
320 1480320 : if (choice == 1 .or. gs_hamkq%usepaw==1) then
321 529344 : call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,vi,v2j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
322 :
323 529344 : if(debug_mode/=0) then
324 24576 : dot2r = dotr ; if (abs(dot2r)<tol9) dot2r = zero ! in order to hide the numerical noise
325 24576 : dot2i = doti ; if (abs(dot2i)<tol9) dot2i = zero ! in order to hide the numerical noise
326 24576 : write(msg,'(3a,2(a,es17.8E3))') bra_i,op2,ket_j,' = ',dot2r,',',dot2i
327 24576 : call wrtout(std_out,msg)
328 : end if
329 :
330 2539008 : rf2%amn(:,iband+(jband-1)*nband_k) = factor_in*(/dotr,doti/) + rf2%amn(:,iband+(jband-1)*nband_k)
331 :
332 : end if ! end choice
333 :
334 : DBG_EXIT("COLL")
335 :
336 1480320 : end subroutine rf2_accumulate_bands
337 : !!***
338 :
339 : !----------------------------------------------------------------------
340 :
341 : !!****f* m_rf2/rf2_apply_hamiltonian
342 : !! NAME
343 : !! rf2_apply_hamiltonian
344 : !!
345 : !! FUNCTION
346 : !! Apply the KS Hamiltonian (or derivative) to an input wave function.
347 : !! If asked, it also does some checks.
348 : !!
349 : !! INPUTS
350 : !! rf2 : rf2_t object containing all rf2 data
351 : !! cg_jband (used only if debug_mode/=0) : array containing |u^(0)(jband)> for all bands
352 : !! cprj_jband(natom,nspinor*usecprj)= u^(0) wave functions for all bands
353 : !! projected with non-local projectors: cprj_jband=<p_i|u^(0)(jband)>
354 : !! cwave(2,size_wf) : input wave function |u>
355 : !! cwaveprj(natom,nspinor) : input wave function |u>
356 : !! projected with non-local projectors <p_i|u>
357 : !! eig0 : 0-order eigenvalue for the present wavefunction at k
358 : !! [enl]=optional (if not present, use hamk%ekb); non-local coeffs connecting projectors
359 : !! eig1_k_jband : first-order lagrange multipliers for the band j (Lambda^(1)_ji for all i)
360 : !! [ffnl1]=nonlocal form factors (needed for tests of ipert>=natom+10)
361 : !! [ffnl1_test]=nonlocal form factors used for tests (needed for tests of ipert>=natom+10)
362 : !! jband : band index of the input wavefunction
363 : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
364 : !! gvnlx1(2,npw1*nspinor)= part of <G|K^(1)+Vnl^(1)|C> not depending on VHxc^(1) (sij_opt/=-1)
365 : !! or part of <G|K^(1)+Vnl^(1)-lambda.S^(1)|C> not depending on VHxc^(1) (sij_opt==-1)
366 : !! idir=direction of the perturbation
367 : !! ipert=type of the perturbation of the Hamiltonian :
368 : !! ipert = 0 : GS calculation, call of getghc
369 : !! ipert = natom+1,2 : 1st order calculation, call of getgh1c
370 : !! ipert = natom+10,11 : 2nd order calculation, call of getgh2c
371 : !! ikpt=number of the k-point
372 : !! isppol=1 index of current spin component
373 : !! mkmem =number of k points trated by this node (GS data).
374 : !! mpi_enreg=information about MPI parallelization
375 : !! nband_k=number of bands at this k point for that spin polarization
376 : !! nsppol=1 for unpolarized, 2 for spin-polarized
377 : !! debug_mode : if /=0 : some tests are done (see NOTES below). Wrong results are printed in std_out
378 : !! prtvol=control print volume and debugging output (for getghc)
379 : !! rf_hamk_idir <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q (here q=0)
380 : !! size_cprj=size of a cprj array (=gs_hamkq%nspinor)
381 : !! size_wf=size of a wavefunction array (=gs_hamkq%npw_k*gs_hamkq%nspinor)
382 : !!
383 : !! OUTPUT
384 : !! h_cwave(2,size_wf) : array containing H^(ipert)|cwave>
385 : !! s_cwave(2,size_wf) : array containing S^(ipert)|cwave>
386 : !!
387 : !! NOTES
388 : !! * Tests are done if debug_mode/=0. In that case : cg_jband(:,jband,1) = |u^(0)(jband)>
389 : !! For ipert==natom+2 (electric field) : cg_jband(:,jband,2) = i * |du/dk_idir(jband)>
390 : !! According to ipert, we check that :
391 : !! -- ipert = 0 :
392 : !! < u^(0)(jband)| H^(0) | u^(0)(jband)> = eig0(jband)
393 : !! -- ipert = natom+1 or 2 :
394 : !! < u^(0)(jband)| ( H^(1) - (eig0(jband)+eig0(iband))/2 * S^(1) ) | u^(0)(iband)> = Lambda^(1)(jband,iband)
395 : !! --ipert >= natom+10 :
396 : !! < u^(0) | H^(2) | u^(0) > from getgh2c is equal to nonlop with signs=1
397 : !! * Use of cprj :
398 : !! The cprj array is computed in dfpt_looppert.F90. In the context of rf2 calculation, it contains
399 : !! <Proj_i^(0)|u^(0)> and <Proj_i^(1)|u^(0)> (dir=1,2 and 3) for all wavefunctions u^(0) or
400 : !! <Proj_i^(0)|u^(1)> and <Proj_i^(1)|u^(1)> (dir=1,2 and 3) for all 1st-order wavefunctions u^(1).
401 : !! Note that <Proj_i^(2)|u^(0)> is always computed on the fly.
402 : !!
403 : !! SOURCE
404 :
405 393312 : subroutine rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave,cwaveprj,h_cwave,s_cwave,eig0,eig1_k_jband,&
406 393312 : & jband,gs_hamkq,gvnlx1,idir,ipert,ikpt,isppol,mkmem,mpi_enreg,nband_k,nsppol,&
407 : & debug_mode,prtvol,rf_hamk_idir,size_cprj,size_wf,&
408 409440 : & conj,enl,ffnl1,ffnl1_test) ! optional
409 :
410 : !Arguments ---------------------------------------------
411 : !scalars
412 : logical,intent(in),optional :: conj
413 : integer,intent(in) :: idir,ipert,ikpt,isppol,jband,mkmem,nband_k,nsppol,debug_mode,prtvol,size_wf,size_cprj
414 : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
415 : ! type(rf2_t),intent(in) :: rf2
416 : type(rf_hamiltonian_type),intent(inout),target :: rf_hamk_idir
417 : type(MPI_type),intent(in) :: mpi_enreg
418 :
419 : !arrays
420 : real(dp),intent(in),target :: cg_jband(2,size_wf*debug_mode*nband_k,2)
421 : real(dp),intent(in),optional,target :: enl(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,gs_hamkq%dimekbq)
422 : real(dp),intent(in),optional :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:)
423 : real(dp),intent(in) :: eig0(nband_k),eig1_k_jband(2*nband_k)
424 : real(dp),intent(inout) :: gvnlx1(2,size_wf)
425 : real(dp),intent(inout) :: cwave(2,size_wf),h_cwave(2,size_wf),s_cwave(2,size_wf)
426 : type(pawcprj_type),intent(inout),target :: cprj_jband(:,:),cwaveprj(:,:)
427 :
428 : !Local variables ---------------------------------------
429 : !scalars
430 : integer,parameter :: berryopt=0,tim_getghc=0,tim_getgh1c=0,tim_getgh2c=0,tim_nonlop=0
431 : integer,parameter :: alpha(9)=(/1,2,3,2,1,1,3,3,2/),beta(9)=(/1,2,3,3,3,2,2,1,1/)
432 : integer :: choice,cpopt,iatom,iband,idirc,idir1,idir2,idir_dum,natom,nnlout,paw_opt,signs,sij_opt
433 : integer :: opt_gvnlx1,opt_gvnl2,optlocal,optnl,usevnl
434 : logical :: compute_conjugate,has_cprj_jband,has_cwaveprj,pert_phon_elfd
435 786624 : real(dp) :: dotr,doti,dotr2,doti2,enlout(18*gs_hamkq%natom),tol_test
436 393312 : real(dp), pointer :: enl_ptr(:,:,:,:)
437 393312 : real(dp),allocatable,target :: enl_temp(:,:,:,:)
438 : character(len=500) :: msg
439 :
440 : !arrays
441 : real(dp),target :: cwave_empty(0,0),svectout_dum(0,0)
442 393312 : real(dp),allocatable :: gvnlxc(:,:),iddk(:,:)
443 393312 : real(dp), ABI_CONTIGUOUS pointer :: cwave_i(:,:),cwave_j(:,:)
444 393312 : type(pawcprj_type),target :: cprj_empty(0,0)
445 393312 : type(pawcprj_type),pointer :: cprj_j(:,:)
446 : ! *********************************************************************
447 :
448 : DBG_ENTER("COLL")
449 :
450 393312 : compute_conjugate = .false.
451 0 : if(present(conj)) compute_conjugate = conj
452 :
453 : ABI_UNUSED(ikpt)
454 : ABI_UNUSED(isppol)
455 : ABI_UNUSED(mkmem)
456 : ABI_UNUSED(nsppol)
457 :
458 : !Check sizes
459 1179936 : if (size(cprj_jband)/=0) then
460 18432 : if (size(cprj_jband)/=nband_k*gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj) then
461 : write(msg,'(2(a,i10))') &
462 0 : 'Wrong cprj size! actual size = ',size(cprj_jband),&
463 0 : ' good size = ',nband_k*gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj
464 0 : ABI_BUG(msg)
465 : end if
466 : end if
467 1179936 : if (size(cwaveprj)/=0) then
468 49536 : if (size(cwaveprj)/=gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj) then
469 : write(msg,'(2(a,i10))') &
470 0 : 'Wrong cwaveprj size! actual size = ',size(cwaveprj),&
471 0 : ' good size = ',gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj
472 0 : ABI_BUG(msg)
473 : end if
474 : end if
475 :
476 393312 : natom = gs_hamkq%natom
477 :
478 393312 : pert_phon_elfd = .false.
479 393312 : if (ipert>natom+11.and.ipert<=2*natom+11) pert_phon_elfd = .true.
480 393312 : usevnl = 0
481 393312 : opt_gvnlx1 = 0
482 393312 : opt_gvnl2 = 0
483 393312 : optnl = 2
484 393312 : sij_opt=1;if (gs_hamkq%usepaw==0) sij_opt=0
485 393312 : if (ipert==natom+2.or.(ipert==natom+11.and.gs_hamkq%usepaw==1).or.pert_phon_elfd) then
486 78528 : usevnl = 1
487 78528 : opt_gvnlx1 = 2
488 78528 : opt_gvnl2 = 1
489 : end if
490 314784 : optlocal = 1
491 78528 : if (pert_phon_elfd) then
492 16128 : optlocal = 0
493 16128 : sij_opt = 0
494 16128 : optnl = 1
495 : end if
496 393312 : tol_test=tol8
497 :
498 : ! In the PAW case : manage cprj_jband, cwaveprj
499 1179936 : has_cprj_jband=(gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1.and.size(cprj_jband)/=0)
500 1179936 : has_cwaveprj=(gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1.and.size(cwaveprj)/=0)
501 393312 : cprj_j => cprj_empty
502 :
503 393312 : if (debug_mode/=0) then
504 19968 : if (ipert/=0) then
505 18240 : write(msg,'(5(a,i4))') 'RF2 TEST rf2_apply_hamiltonian ipert = ',ipert,' idir =',idir,' isppol =',isppol,&
506 36480 : & ' ikpt = ',ikpt,' jband = ',jband
507 : else
508 1728 : write(msg,'(3(a,i4))') 'RF2 TEST rf2_apply_hamiltonian ipert = 0 idir = 0 isppol =',isppol,&
509 3456 : & ' ikpt = ',ikpt,' jband = ',jband
510 : end if
511 19968 : call wrtout(std_out,msg)
512 : end if
513 :
514 : ! *******************************************************************************************
515 : ! apply H^(0)
516 : ! *******************************************************************************************
517 393312 : if (ipert == 0) then
518 :
519 308016 : ABI_MALLOC(gvnlxc,(2,size_wf))
520 37976508 : gvnlxc(:,:) = zero
521 :
522 : ! Test if < u^(0) | H^(0) | u^(0) > = eig0(jband)
523 102672 : if(debug_mode/=0) then
524 1728 : cwave_j => cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,1)
525 1728 : if (has_cprj_jband) cprj_j => cprj_jband(:,1+(jband-1)*size_cprj:jband*size_cprj)
526 1728 : cpopt = -1+3*gs_hamkq%usecprj*gs_hamkq%usepaw
527 : call getghc(cpopt,cwave_j,cprj_j,h_cwave,s_cwave,gs_hamkq,gvnlxc,zero,mpi_enreg,1,prtvol,&
528 1728 : sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
529 :
530 1728 : call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_j,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
531 1728 : dotr = dotr - eig0(jband)
532 1728 : dotr = sqrt(dotr**2+doti**2)
533 1728 : if (dotr > tol_test) then
534 0 : write(msg,'(a,es17.8E3)') 'RF2 TEST GETGHC : NOT PASSED dotr = ',dotr
535 0 : call wrtout(ab_out,msg)
536 0 : call wrtout(std_out,msg)
537 : end if
538 : end if ! end tests
539 :
540 102672 : cpopt=-1;if (has_cwaveprj) cpopt=2
541 : call getghc(cpopt,cwave,cwaveprj,h_cwave,s_cwave,gs_hamkq,gvnlxc,zero,mpi_enreg,1,prtvol,&
542 102672 : sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
543 102672 : ABI_FREE(gvnlxc)
544 :
545 : ! *******************************************************************************************
546 : ! apply H^(1)
547 : ! *******************************************************************************************
548 290640 : else if (ipert<=natom+2) then
549 :
550 : ! Test if < u^(0) | ( H^(1) - eps^(0) S^(1) ) | u^(0) > = eig^(1)
551 226416 : if(debug_mode/=0) then
552 31104 : ABI_MALLOC(iddk,(2,size_wf))
553 10368 : cwave_j => cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,1)
554 10368 : if (has_cprj_jband) cprj_j => cprj_jband(:,1+(jband-1)*size_cprj:jband*size_cprj)
555 4570560 : iddk(:,:) = zero;if (ipert==natom+2) iddk(:,:)=cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,2)
556 : call getgh1c(berryopt,cwave_j,cprj_j,h_cwave,cwave_empty,s_cwave,gs_hamkq,iddk,idir,ipert,(/zero/),&
557 10368 : mpi_enreg,1,optlocal,optnl,opt_gvnlx1,rf_hamk_idir,sij_opt,tim_getgh1c,usevnl,conj=compute_conjugate)
558 134784 : do iband=1,nband_k
559 124416 : cwave_i => cg_jband(:,1+(iband-1)*size_wf:iband*size_wf,1)
560 124416 : call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_i,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
561 124416 : if (gs_hamkq%usepaw==1.and.ipert/=natom+2) then ! S^(1) is zero for ipert=natom+2
562 55296 : call dotprod_g(dotr2,doti2,gs_hamkq%istwf_k,size_wf,2,cwave_i,s_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
563 55296 : dotr = dotr - (eig0(iband)+eig0(jband))*dotr2/two
564 55296 : doti = doti - (eig0(iband)+eig0(jband))*doti2/two
565 : end if
566 124416 : dotr = dotr - eig1_k_jband(1+2*(iband-1))
567 124416 : doti = doti - eig1_k_jband(2+2*(iband-1))
568 124416 : dotr = sqrt(dotr**2+doti**2)
569 134784 : if (dotr > tol_test) then
570 0 : write(msg,'(4(a,i2),a,es17.8E3)') 'RF2 TEST GETGH1 : ipert=',ipert,' idir=',idir,&
571 0 : ' jband=',jband,' iband=',iband,' NOT PASSED dotr = ',dotr
572 0 : call wrtout(ab_out,msg)
573 0 : call wrtout(std_out,msg)
574 : end if
575 : end do ! end iband
576 10368 : ABI_FREE(iddk)
577 : end if ! end tests
578 :
579 : call getgh1c(berryopt,cwave,cwaveprj,h_cwave,cwave_empty,s_cwave,gs_hamkq,gvnlx1,idir,ipert,(/zero/),&
580 226416 : mpi_enreg,1,optlocal,optnl,opt_gvnlx1,rf_hamk_idir,sij_opt,tim_getgh1c,usevnl,conj=compute_conjugate)
581 :
582 : ! *******************************************************************************************
583 : ! apply H^(2)
584 : ! *******************************************************************************************
585 64224 : else if (ipert==natom+10.or.ipert==natom+11.or.pert_phon_elfd) then
586 :
587 : ! *******************************************************************************************
588 : ! Test if < u^(0) | H^(2) | u^(0) > from getgh2c is equal to nonlop with signs=1
589 64224 : if(debug_mode/=0.and.present(ffnl1).and.present(ffnl1_test)) then
590 :
591 7872 : cwave_j => cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,1)
592 :
593 23616 : ABI_MALLOC(iddk,(2,size_wf))
594 2233680 : iddk(:,:) = cwave_j(:,:)
595 :
596 : call getgh2c(cwave_j,cprj_empty,h_cwave,s_cwave,gs_hamkq,iddk,idir,ipert,zero,&
597 : mpi_enreg,optlocal,optnl,opt_gvnl2,rf_hamk_idir,sij_opt,tim_getgh2c,usevnl,&
598 7872 : conj=compute_conjugate,optkin=0,enl=enl)
599 7872 : ABI_FREE(iddk)
600 :
601 7872 : call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_j,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
602 :
603 7872 : idir1=alpha(idir);idir2=beta(idir)
604 7872 : idirc=3*(idir1-1)+idir2
605 :
606 291264 : enlout = zero
607 :
608 : ! Change the pointer ffnl_k to ffnl1_test (idir_ffnl=0, for nonlop with signs=1)
609 7872 : call gs_hamkq%load_k(ffnl_k=ffnl1_test)
610 :
611 7872 : signs=1
612 7872 : dotr2 = 0
613 7872 : doti2 = 0
614 :
615 : ! ************************
616 : ! IPERT == NATOM+10
617 : ! ************************
618 7872 : if (ipert==natom+10) then
619 :
620 384 : cpopt=-1; choice=8; paw_opt=1; if (gs_hamkq%usepaw==0) paw_opt=0 ; nnlout = 6
621 :
622 : call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idirc,(/zero/),mpi_enreg,1,nnlout,&
623 384 : & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum)
624 384 : dotr2 = enlout(idir)
625 :
626 : end if ! IPERT == NATOM+10
627 :
628 : ! ************************
629 : ! IPERT == NATOM+11
630 : ! ************************
631 7872 : if (ipert==natom+11) then
632 :
633 576 : if (present(enl)) then
634 0 : enl_ptr => enl
635 576 : else if (associated(rf_hamk_idir%e1kbfr).and.associated(rf_hamk_idir%e1kbsc).and.optnl==2) then
636 3456 : ABI_MALLOC(enl_temp,(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,rf_hamk_idir%cplex))
637 44352 : enl_temp(:,:,:,:) = rf_hamk_idir%e1kbfr(:,:,:,:) + rf_hamk_idir%e1kbsc(:,:,:,:)
638 576 : enl_ptr => enl_temp
639 0 : else if (associated(rf_hamk_idir%e1kbfr)) then
640 0 : enl_ptr => rf_hamk_idir%e1kbfr
641 : end if
642 :
643 : ! Compute application of dS/dk1 to cwavef
644 : ! sum_{i,j} s_ij < psi^(0) | d(|p_i><p_j|)/dk(idir1) | psi^(0) >
645 576 : cpopt=-1 ; choice=5 ; paw_opt=3 ; nnlout=3
646 : call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
647 576 : & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum)
648 576 : dotr2 = enlout(idir1)
649 :
650 : ! Compute part of H^(2) due to derivative of projectors (idir1) and derivative of Dij (idir2)
651 : ! sum_{i,j} chi_ij(idir2) < psi^(0) | d(|p_i><p_j|)/dk(idir1) | psi^(0) >
652 576 : cpopt=-1 ; choice=5 ; paw_opt=1 ; nnlout=3
653 : call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
654 576 : & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum,enl=enl_ptr)
655 576 : dotr2 = dotr2 + enlout(idir1)
656 :
657 : ! Compute derivatives due to projectors |d^2[p_i]/dk1dk2>,|d[p_i]/dk1>,|d[p_i]/dk2>
658 : ! i * sum_{i,j} < psi^(0) | (d(|p_i><dp_j/dk(idir2)|)/dk(idir1) | psi^(0) >
659 576 : cpopt=-1 ; choice=81 ; paw_opt=3 ; nnlout=18
660 : call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
661 576 : & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum)
662 576 : dotr2 = dotr2 - enlout(2*idirc )
663 576 : doti2 = doti2 + enlout(2*idirc-1)
664 :
665 : end if ! IPERT == NATOM+11
666 :
667 : ! ************************
668 : ! PERT_PHON_ELFD
669 : ! ************************
670 7872 : if (pert_phon_elfd) then
671 :
672 6912 : if (present(enl)) then
673 6912 : enl_ptr => enl
674 0 : else if (associated(rf_hamk_idir%e1kbfr).and.associated(rf_hamk_idir%e1kbsc).and.optnl==2) then
675 0 : ABI_MALLOC(enl_temp,(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,rf_hamk_idir%cplex))
676 0 : enl_temp(:,:,:,:) = rf_hamk_idir%e1kbfr(:,:,:,:) + rf_hamk_idir%e1kbsc(:,:,:,:)
677 0 : enl_ptr => enl_temp
678 0 : else if (associated(rf_hamk_idir%e1kbfr)) then
679 0 : enl_ptr => rf_hamk_idir%e1kbfr
680 : end if
681 :
682 6912 : iatom = gs_hamkq%atindx(ipert - natom - 11) ! Atoms are type-sorted in enlout()
683 :
684 : ! Compute application of dS/dtau1 to i*d[cwavef]/dk2
685 : ! sum_{i,j} s_ij < psi^(0) | d(|p_i><p_j|)/dtau(idir1) | i*psi^(k(idir2)) >
686 6912 : cpopt=-1 ; choice=2 ; paw_opt=3 ; nnlout=3*natom
687 : call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
688 6912 : & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum)
689 6912 : dotr2 = enlout(3*(iatom-1)+idir1)
690 :
691 : ! Compute part of H^(2) due to derivative of projectors (idir1) and derivative of Dij (idir2)
692 : ! sum_{i,j} chi_ij(idir2) < psi^(0) | d(|p_i><p_j|)/dtau(idir1) | psi^(0) >
693 6912 : cpopt=-1 ; choice=2 ; paw_opt=1 ; nnlout=3*natom
694 : call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
695 6912 : & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum,enl=enl_ptr)
696 6912 : dotr2 = dotr2 + enlout(3*(iatom-1)+idir1)
697 :
698 : ! Compute derivatives due to projectors |d^2[p_i]/dtau1dk2>,|d[p_i]/dtau1>,|d[p_i]/dk2>
699 : ! i * sum_{i,j} < psi^(0) | (d(|p_i><dp_j/dk(idir2)|)/dtau(idir1) | psi^(0) >
700 6912 : cpopt=-1 ; choice=54 ; paw_opt=3 ; nnlout=18*natom
701 : call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
702 6912 : & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum)
703 6912 : dotr2 = dotr2 - enlout(18*(iatom-1)+2*idirc )
704 6912 : doti2 = doti2 + enlout(18*(iatom-1)+2*idirc-1)
705 :
706 : end if ! PERT_PHON_ELFD
707 :
708 7872 : dotr = dotr - dotr2
709 7872 : doti = doti - doti2
710 7872 : dotr = sqrt(dotr**2+doti**2)
711 7872 : if (dotr > 10*tol_test) then
712 0 : write(msg,'(a,es17.8E3)') 'RF2 TEST GETGH2C : NOT PASSED dotr = ',dotr
713 0 : call wrtout(ab_out,msg)
714 0 : call wrtout(std_out,msg)
715 : end if
716 :
717 : ! Change the pointer ffnl_k back to ffnl1 (idir_ffnl=4, for nonlop with signs=2)
718 7872 : call gs_hamkq%load_k(ffnl_k=ffnl1)
719 :
720 7872 : if (associated(enl_ptr)) then
721 7871 : nullify(enl_ptr)
722 : end if
723 7872 : ABI_SFREE(enl_temp)
724 :
725 : end if ! end tests
726 : ! *******************************************************************************************
727 :
728 : call getgh2c(cwave,cwaveprj,h_cwave,s_cwave,gs_hamkq,gvnlx1,idir,ipert,zero,&
729 : mpi_enreg,optlocal,optnl,opt_gvnl2,rf_hamk_idir,sij_opt,tim_getgh2c,usevnl,&
730 64224 : conj=compute_conjugate,enl=enl)
731 :
732 : else
733 :
734 0 : h_cwave = zero
735 0 : ABI_ERROR(" rf2_apply_hamiltonian can be used only for: 0<=ipert<=natom+2 and natom+10<=ipert<=2*natom+11.")
736 0 : return
737 :
738 : end if
739 :
740 : DBG_EXIT("COLL")
741 :
742 802752 : end subroutine rf2_apply_hamiltonian
743 : !!***
744 :
745 : !----------------------------------------------------------------------
746 :
747 : !!****f* m_rf2/rf2_destroy
748 : !! NAME
749 : !! rf2_init
750 : !!
751 : !! FUNCTION
752 : !! Free all allocated arrays in a rf2_t object
753 : !!
754 : !! SOURCE
755 :
756 14616 : subroutine rf2_destroy(rf2)
757 :
758 : !Arguments ---------------------------------------------
759 : !scalars
760 : type(rf2_t),intent(inout) :: rf2
761 : ! *************************************************************************
762 :
763 14616 : if (associated(rf2%RHS_Stern)) then
764 14616 : ABI_FREE(rf2%RHS_Stern)
765 : end if
766 14616 : ABI_SFREE(rf2%dcwavef)
767 14616 : ABI_SFREE(rf2%amn)
768 14616 : ABI_SFREE(rf2%lambda_mn)
769 :
770 14616 : end subroutine rf2_destroy
771 : !!***
772 :
773 0 : END MODULE m_rf2
|
