Branch data 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), 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), 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
# # # # #
# # # # #
# # # # ]
|
