TABLE OF CONTENTS
- ABINIT/dfpt_ciflexoout
- ABINIT/dfpt_ddmdqout
- ABINIT/dfpt_flexo
- ABINIT/dfpt_isdqout
- ABINIT/dfpt_qdrpole
- ABINIT/dfpt_qdrpout
- ABINIT/m_dfpt_lw
ABINIT/dfpt_ciflexoout [ Functions ]
NAME
dfpt_ciflexoout
FUNCTION
This subroutine gathers the different terms entering the electrocic flexoelectric tensor, perfofms the transformation from reduced to cartesian coordinates and writes out the tensor in output files.
INPUTS
elqgradhart(2,3,3,3,3)=electronic electrostatic contribution from the q-gradient of the Hartree potential elflexoflg(3,3,3,3)=array that indicates which elements of the electronic contribution to the flexoelectric tensor have been calculated elflexowf(2,3,3,3,3)=total wave function contributions to the electronic flexoelectric tensor (except t4) elflexowf_t1(2,3,3,3,3)=term 1 of the wave function contribution elflexowf_t2(2,3,3,3,3)=term 2 of the wave function contribution elflexowf_t3(2,3,3,3,3)=term 3 of the wave function contribution elflexowf_t4(2,3,3,3,3)=term 4 of the wave function contribution (type-I) elflexowf_t5(2,3,3,3,3)=term 5 of the wave function contribution gprimd(3,3)=reciprocal space dimensional primitive translations kptopt=2 time reversal symmetry is enforced, 3 trs is not enforced (for debugging purposes) matom=number of atoms mpert=maximum number of perturbations nefipert=number of electric field perturbations nstrpert=number of strain perturbations nq1grad=number of q1 (q_{\gamma}) gradients pert_efield(3,nefipert)=array with the info for the electric field perturbations pert_strain(6,nstrpert)=array with the info for the strain perturbations prtvol=volume of information to be printed. 1-> The different contributions to the quadrupole are printed. q1grad(3,nq1grad)=array with the info for the q1 (q_{\gamma}) gradients rprimd(3,3)=dimensional primitive translations in real space (bohr) ucvol=unit cell volume in bohr**3.
OUTPUT
d3etot(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE
SIDE EFFECTS
NOTES
SOURCE
3131 subroutine dfpt_ciflexoout(d3etot,elflexoflg,elflexowf,elflexowf_t1,elflexowf_t2, & 3132 & elflexowf_t3,elflexowf_t4,elflexowf_t5, & 3133 & elqgradhart,gprimd,kptopt,matom,mpert,nefipert,& 3134 & nstrpert,nq1grad,pert_efield,pert_strain,prtvol,q1grad,rprimd,ucvol) 3135 3136 !Arguments ------------------------------------ 3137 !scalars 3138 integer,intent(in) :: kptopt,matom,mpert,nefipert,nstrpert,nq1grad,prtvol 3139 real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert) 3140 real(dp),intent(in) :: ucvol 3141 3142 !arrays 3143 integer,intent(in) :: elflexoflg(3,3,3,3) 3144 integer,intent(in) :: pert_efield(3,nefipert) 3145 integer,intent(in) :: pert_strain(6,nstrpert) 3146 integer,intent(in) :: q1grad(3,nq1grad) 3147 real(dp),intent(in) :: elflexowf(2,3,3,3,3) 3148 real(dp),intent(inout) :: elflexowf_t1(2,3,3,3,3) 3149 real(dp),intent(inout) :: elflexowf_t2(2,3,3,3,3) 3150 real(dp),intent(inout) :: elflexowf_t3(2,3,3,3,3) 3151 real(dp),intent(inout) :: elflexowf_t4(2,3,3,3,3) 3152 real(dp),intent(inout) :: elflexowf_t5(2,3,3,3,3) 3153 real(dp),intent(inout) :: elqgradhart(2,3,3,3,3) 3154 real(dp),intent(in) :: gprimd(3,3) 3155 real(dp),intent(in) :: rprimd(3,3) 3156 3157 !Local variables------------------------------- 3158 !scalars 3159 integer :: alpha,beta,delta,efipert,gamma 3160 integer :: ibuf,iefidir,iefipert,ii,iq1dir,iq1grad,istr1dir,istr2dir,istrpert 3161 integer :: q1pert,strcomp,strpert 3162 integer, parameter :: re=1,im=2 3163 real(dp) :: fac,tmpim,tmpre,ucvolinv 3164 character(len=500) :: msg 3165 3166 !arrays 3167 integer,allocatable :: cartflg_t4(:,:,:,:),redflg(:,:,:,:) 3168 integer :: flg1(3),flg2(3) 3169 real(dp) :: vec1(3),vec2(3) 3170 real(dp),allocatable :: elec_flexotens_cart(:,:,:,:,:),elec_flexotens_red(:,:,:,:,:) 3171 real(dp),allocatable :: elflexowf_buffer_cart(:,:,:,:,:,:),elflexowf_t4_cart(:,:,:,:,:) 3172 3173 ! ************************************************************************* 3174 3175 DBG_ENTER("COLL") 3176 3177 !Gather the different terms in the electronic contribution to the flexoelectric tensor 3178 ABI_MALLOC(elec_flexotens_red,(2,3,3,3,3)) 3179 ! elec_flexotens_red=zero 3180 ucvolinv= 1.0_dp/ucvol 3181 3182 if (kptopt==3) then 3183 3184 !Compute real and 'true' imaginary parts of flexoelectric tensor and independent terms 3185 !T4 term needs further treatment and it will be lately added to the cartesian coordinates 3186 !version of the flexoelectric tensor 3187 do istrpert=1,nstrpert 3188 istr1dir=pert_strain(3,istrpert) 3189 istr2dir=pert_strain(4,istrpert) 3190 do iq1grad=1,nq1grad 3191 iq1dir=q1grad(2,iq1grad) 3192 do iefipert=1,nefipert 3193 iefidir=pert_efield(2,iefipert) 3194 3195 if (elflexoflg(iefidir,iq1dir,istr1dir,istr2dir)==1) then 3196 3197 elec_flexotens_red(re,iefidir,iq1dir,istr1dir,istr2dir)=2.0_dp*ucvolinv* & 3198 & ( elqgradhart(re,iefidir,iq1dir,istr1dir,istr2dir) + & 3199 & elflexowf(re,iefidir,iq1dir,istr1dir,istr2dir) ) 3200 elec_flexotens_red(im,iefidir,iq1dir,istr1dir,istr2dir)=2.0_dp*ucvolinv* & 3201 & ( elqgradhart(im,iefidir,iq1dir,istr1dir,istr2dir) + & 3202 & elflexowf(im,iefidir,iq1dir,istr1dir,istr2dir) ) 3203 3204 !Multiply by the imaginary unit that has been factorized out 3205 tmpre=elec_flexotens_red(re,iefidir,iq1dir,istr1dir,istr2dir) 3206 tmpim=elec_flexotens_red(im,iefidir,iq1dir,istr1dir,istr2dir) 3207 elec_flexotens_red(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim 3208 elec_flexotens_red(im,iefidir,iq1dir,istr1dir,istr2dir)=tmpre 3209 3210 !Divide by the imaginary unit the T4 term 3211 !(this is because a -i factor has been factorized out in all the individual contributions to this therm) 3212 tmpre=elflexowf_t4(re,iefidir,istr1dir,istr2dir,iq1dir) 3213 tmpim=elflexowf_t4(im,iefidir,istr1dir,istr2dir,iq1dir) 3214 elflexowf_t4(re,iefidir,istr1dir,istr2dir,iq1dir)=tmpim*2.0_dp*ucvolinv 3215 elflexowf_t4(im,iefidir,istr1dir,istr2dir,iq1dir)=-tmpre*2.0_dp*ucvolinv 3216 3217 !Compute and save individual terms in mixed coordinates 3218 if (prtvol>=10) then 3219 3220 tmpre=elqgradhart(re,iefidir,iq1dir,istr1dir,istr2dir) 3221 tmpim=elqgradhart(im,iefidir,iq1dir,istr1dir,istr2dir) 3222 elqgradhart(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim*2.0_dp*ucvolinv 3223 elqgradhart(im,iefidir,iq1dir,istr1dir,istr2dir)=tmpre*2.0_dp*ucvolinv 3224 3225 tmpre=elflexowf_t1(re,iefidir,iq1dir,istr1dir,istr2dir) 3226 tmpim=elflexowf_t1(im,iefidir,iq1dir,istr1dir,istr2dir) 3227 elflexowf_t1(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim*2.0_dp*ucvolinv 3228 elflexowf_t1(im,iefidir,iq1dir,istr1dir,istr2dir)=tmpre*2.0_dp*ucvolinv 3229 3230 tmpre=elflexowf_t2(re,iefidir,iq1dir,istr1dir,istr2dir) 3231 tmpim=elflexowf_t2(im,iefidir,iq1dir,istr1dir,istr2dir) 3232 elflexowf_t2(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim*2.0_dp*ucvolinv 3233 elflexowf_t2(im,iefidir,iq1dir,istr1dir,istr2dir)=tmpre*2.0_dp*ucvolinv 3234 3235 tmpre=elflexowf_t3(re,iefidir,iq1dir,istr1dir,istr2dir) 3236 tmpim=elflexowf_t3(im,iefidir,iq1dir,istr1dir,istr2dir) 3237 elflexowf_t3(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim*2.0_dp*ucvolinv 3238 elflexowf_t3(im,iefidir,iq1dir,istr1dir,istr2dir)=tmpre*2.0_dp*ucvolinv 3239 3240 tmpre=elflexowf_t5(re,iefidir,iq1dir,istr1dir,istr2dir) 3241 tmpim=elflexowf_t5(im,iefidir,iq1dir,istr1dir,istr2dir) 3242 elflexowf_t5(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim*2.0_dp*ucvolinv 3243 elflexowf_t5(im,iefidir,iq1dir,istr1dir,istr2dir)=tmpre*2.0_dp*ucvolinv 3244 3245 end if 3246 3247 end if 3248 3249 end do 3250 end do 3251 end do 3252 3253 else if (kptopt==2) then 3254 3255 !Compute real part of flexoelectric tensor and independent terms 3256 !T4 term needs further treatment and it will be lately added to the cartesian coordinates 3257 !version of the flexoelectric tensor 3258 do istrpert=1,nstrpert 3259 istr1dir=pert_strain(3,istrpert) 3260 istr2dir=pert_strain(4,istrpert) 3261 do iq1grad=1,nq1grad 3262 iq1dir=q1grad(2,iq1grad) 3263 do iefipert=1,nefipert 3264 iefidir=pert_efield(2,iefipert) 3265 3266 if (elflexoflg(iefidir,iq1dir,istr1dir,istr2dir)==1) then 3267 3268 elec_flexotens_red(re,iefidir,iq1dir,istr1dir,istr2dir)=2.0_dp*ucvolinv* & 3269 & ( elqgradhart(re,iefidir,iq1dir,istr1dir,istr2dir) + & 3270 & elflexowf(re,iefidir,iq1dir,istr1dir,istr2dir) ) 3271 elec_flexotens_red(im,iefidir,iq1dir,istr1dir,istr2dir)=2.0_dp*ucvolinv* & 3272 & ( elqgradhart(im,iefidir,iq1dir,istr1dir,istr2dir) + & 3273 & elflexowf(im,iefidir,iq1dir,istr1dir,istr2dir) ) 3274 3275 !Multiply by the imaginary unit that has been factorized out 3276 tmpim=elec_flexotens_red(im,iefidir,iq1dir,istr1dir,istr2dir) 3277 elec_flexotens_red(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim 3278 elec_flexotens_red(im,iefidir,iq1dir,istr1dir,istr2dir)=0.0_dp 3279 3280 !Divide by the imaginary unit the T4 term 3281 !(this is because a -i factor has been factorized out in all the individual contributions to this therm) 3282 tmpim=elflexowf_t4(im,iefidir,istr1dir,istr2dir,iq1dir) 3283 elflexowf_t4(re,iefidir,istr1dir,istr2dir,iq1dir)=2.0_dp*tmpim*ucvolinv 3284 elflexowf_t4(im,iefidir,istr1dir,istr2dir,iq1dir)=0.0_dp 3285 3286 !Compute and save individual terms in mixed coordinates 3287 if (prtvol>=10) then 3288 3289 tmpim=elqgradhart(im,iefidir,iq1dir,istr1dir,istr2dir) 3290 elqgradhart(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim*2.0_dp*ucvolinv 3291 elqgradhart(im,iefidir,iq1dir,istr1dir,istr2dir)=0.0_dp 3292 3293 tmpim=elflexowf_t1(im,iefidir,iq1dir,istr1dir,istr2dir) 3294 elflexowf_t1(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim*2.0_dp*ucvolinv 3295 elflexowf_t1(im,iefidir,iq1dir,istr1dir,istr2dir)=0.0_dp 3296 3297 tmpim=elflexowf_t2(im,iefidir,iq1dir,istr1dir,istr2dir) 3298 elflexowf_t2(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim*2.0_dp*ucvolinv 3299 elflexowf_t2(im,iefidir,iq1dir,istr1dir,istr2dir)=0.0_dp 3300 3301 tmpim=elflexowf_t3(im,iefidir,iq1dir,istr1dir,istr2dir) 3302 elflexowf_t3(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim*2.0_dp*ucvolinv 3303 elflexowf_t3(im,iefidir,iq1dir,istr1dir,istr2dir)=0.0_dp 3304 3305 tmpim=elflexowf_t5(im,iefidir,iq1dir,istr1dir,istr2dir) 3306 elflexowf_t5(re,iefidir,iq1dir,istr1dir,istr2dir)=-tmpim*2.0_dp*ucvolinv 3307 elflexowf_t5(im,iefidir,iq1dir,istr1dir,istr2dir)=0.0_dp 3308 3309 end if 3310 3311 end if 3312 3313 end do 3314 end do 3315 end do 3316 3317 else 3318 write(msg,"(1a)") 'kptopt must be 2 or 3 for long-wave DFPT calculations' 3319 ABI_BUG(msg) 3320 end if 3321 3322 !Transormation to complete cartesian coordinates the flexoelectric tensor 3323 !and separately the T4 term 3324 ABI_MALLOC(elec_flexotens_cart,(2,3,3,3,3)) 3325 ABI_MALLOC(elflexowf_t4_cart,(2,3,3,3,3)) 3326 ABI_MALLOC(cartflg_t4,(3,3,3,3)) 3327 elec_flexotens_cart=elec_flexotens_red 3328 elflexowf_t4_cart=elflexowf_t4 3329 cartflg_t4=0 3330 3331 ! ABI_FREE(elec_flexotens_red) 3332 3333 if (prtvol>=10) then 3334 ABI_MALLOC(elflexowf_buffer_cart,(5,2,3,3,3,3)) 3335 elflexowf_buffer_cart(1,:,:,:,:,:)=elflexowf_t1(:,:,:,:,:) 3336 elflexowf_buffer_cart(2,:,:,:,:,:)=elflexowf_t2(:,:,:,:,:) 3337 elflexowf_buffer_cart(3,:,:,:,:,:)=elflexowf_t3(:,:,:,:,:) 3338 elflexowf_buffer_cart(4,:,:,:,:,:)=elqgradhart(:,:,:,:,:) 3339 elflexowf_buffer_cart(5,:,:,:,:,:)=elflexowf_t5(:,:,:,:,:) 3340 end if 3341 3342 !1st transform coordinates of the electric field derivative of the flexoelectric tensor 3343 do istr2dir=1,3 3344 do istr1dir=1,3 3345 do iq1dir=1,3 3346 do ii=1,2 3347 do iefidir=1,3 3348 vec1(iefidir)=elec_flexotens_cart(ii,iefidir,iq1dir,istr1dir,istr2dir) 3349 flg1(iefidir)=elflexoflg(iefidir,iq1dir,istr1dir,istr2dir) 3350 end do 3351 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 3352 do iefidir=1,3 3353 elec_flexotens_cart(ii,iefidir,iq1dir,istr1dir,istr2dir)=vec2(iefidir) 3354 end do 3355 end do 3356 end do 3357 end do 3358 end do 3359 3360 !Do now the transformation of the electric field derivative for the T4 term 3361 do istr2dir=1,3 3362 do istr1dir=1,3 3363 do iq1dir=1,3 3364 do ii=1,2 3365 do iefidir=1,3 3366 vec1(iefidir)=elflexowf_t4_cart(ii,iefidir,istr1dir,istr2dir,iq1dir) 3367 flg1(iefidir)=elflexoflg(iefidir,iq1dir,istr1dir,istr2dir) 3368 end do 3369 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 3370 do iefidir=1,3 3371 elflexowf_t4_cart(ii,iefidir,istr1dir,istr2dir,iq1dir)=vec2(iefidir) 3372 cartflg_t4(iefidir,istr1dir,istr2dir,iq1dir)=flg2(iefidir) 3373 end do 3374 end do 3375 end do 3376 end do 3377 end do 3378 3379 !Do now the transformation of the electric field derivative for the other therms 3380 if (prtvol>=10) then 3381 do ibuf=1,5 3382 do istr2dir=1,3 3383 do istr1dir=1,3 3384 do iq1dir=1,3 3385 do ii=1,2 3386 do iefidir=1,3 3387 vec1(iefidir)=elflexowf_buffer_cart(ibuf,ii,iefidir,iq1dir,istr1dir,istr2dir) 3388 flg1(iefidir)=elflexoflg(iefidir,iq1dir,istr1dir,istr2dir) 3389 end do 3390 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 3391 do iefidir=1,3 3392 elflexowf_buffer_cart(ibuf,ii,iefidir,iq1dir,istr1dir,istr2dir)=vec2(iefidir) 3393 end do 3394 end do 3395 end do 3396 end do 3397 end do 3398 end do 3399 end if 3400 3401 !2nd transform coordinates of the q-gradient (treat it as electric field) 3402 !of the flexoelectric tensor 3403 do istr2dir=1,3 3404 do istr1dir=1,3 3405 do iefidir=1,3 3406 do ii=1,2 3407 do iq1dir=1,3 3408 vec1(iq1dir)=elec_flexotens_cart(ii,iefidir,iq1dir,istr1dir,istr2dir) 3409 flg1(iq1dir)=elflexoflg(iefidir,iq1dir,istr1dir,istr2dir) 3410 end do 3411 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 3412 do iq1dir=1,3 3413 elec_flexotens_cart(ii,iefidir,iq1dir,istr1dir,istr2dir)=vec2(iq1dir) 3414 end do 3415 end do 3416 end do 3417 end do 3418 end do 3419 3420 !2nd transform coordinates of the q-gradient (treat it as electric field) 3421 !of the individual terms 3422 if (prtvol>=10) then 3423 do ibuf=1,5 3424 do istr2dir=1,3 3425 do istr1dir=1,3 3426 do iefidir=1,3 3427 do ii=1,2 3428 do iq1dir=1,3 3429 vec1(iq1dir)=elflexowf_buffer_cart(ibuf,ii,iefidir,iq1dir,istr1dir,istr2dir) 3430 flg1(iq1dir)=elflexoflg(iefidir,iq1dir,istr1dir,istr2dir) 3431 end do 3432 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 3433 do iq1dir=1,3 3434 elflexowf_buffer_cart(ibuf,ii,iefidir,iq1dir,istr1dir,istr2dir)=vec2(iq1dir) 3435 end do 3436 end do 3437 end do 3438 end do 3439 end do 3440 end do 3441 end if 3442 3443 !Write the flexoelectric tensor in cartesian coordinates 3444 !Open output files 3445 if (prtvol>=10) then 3446 open(unit=71,file='elec_flexo_wf_t1.out',status='unknown',form='formatted',action='write') 3447 open(unit=72,file='elec_flexo_wf_t2.out',status='unknown',form='formatted',action='write') 3448 open(unit=73,file='elec_flexo_wf_t3.out',status='unknown',form='formatted',action='write') 3449 open(unit=74,file='elec_flexo_wf_t4.out',status='unknown',form='formatted',action='write') 3450 open(unit=75,file='elec_flexo_wf_t5.out',status='unknown',form='formatted',action='write') 3451 open(unit=76,file='elec_flexo_elecstic.out',status='unknown',form='formatted',action='write') 3452 end if 3453 3454 write(ab_out,'(a)')' ' 3455 write(ab_out,'(a)')' Clamped-ion flexoelectric tensor, in cartesian coordinates,' 3456 write(ab_out,'(a)')' efidir qgrdir strdir1 strdir2 real part imaginary part' 3457 do istr2dir=1,3 3458 delta=istr2dir 3459 do istr1dir=1,3 3460 beta=istr1dir 3461 do iq1dir=1,3 3462 gamma=iq1dir 3463 do iefidir=1,3 3464 alpha=iefidir 3465 3466 if (cartflg_t4(alpha,beta,delta,gamma)==1 .and. cartflg_t4(alpha,delta,gamma,beta)==1 & 3467 & .and. cartflg_t4(alpha,gamma,beta,delta)==1) then 3468 3469 !Converts the T4 term to type-II form 3470 tmpre= elflexowf_t4_cart(re,alpha,beta,delta,gamma) + & 3471 & elflexowf_t4_cart(re,alpha,delta,gamma,beta) - & 3472 & elflexowf_t4_cart(re,alpha,gamma,beta,delta) 3473 3474 tmpim= elflexowf_t4_cart(im,alpha,beta,delta,gamma) + & 3475 & elflexowf_t4_cart(im,alpha,delta,gamma,beta) - & 3476 & elflexowf_t4_cart(im,alpha,gamma,beta,delta) 3477 3478 !Add the T4 term after conversion to type-II form 3479 elec_flexotens_cart(re,alpha,gamma,beta,delta)= & 3480 & elec_flexotens_cart(re,alpha,gamma,beta,delta) + tmpre 3481 elec_flexotens_cart(im,alpha,gamma,beta,delta)= & 3482 & elec_flexotens_cart(im,alpha,gamma,beta,delta) + tmpim 3483 3484 !Writes the complete flexoelectric tensor 3485 write(ab_out,'(4(i5,3x),2(1x,f20.10))') alpha,gamma,beta,delta, & 3486 & elec_flexotens_cart(re,alpha,gamma,beta,delta), & 3487 & elec_flexotens_cart(im,alpha,gamma,beta,delta) 3488 3489 if (prtvol>=10) then 3490 write(71,'(4(i5,3x),2(1x,f20.10))') alpha,gamma,beta,delta, & 3491 & elflexowf_buffer_cart(1,re,alpha,gamma,beta,delta), & 3492 & elflexowf_buffer_cart(1,im,alpha,gamma,beta,delta) 3493 3494 write(72,'(4(i5,3x),2(1x,f20.10))') alpha,gamma,beta,delta, & 3495 & elflexowf_buffer_cart(2,re,alpha,gamma,beta,delta), & 3496 & elflexowf_buffer_cart(2,im,alpha,gamma,beta,delta) 3497 3498 write(73,'(4(i5,3x),2(1x,f20.10))') alpha,gamma,beta,delta, & 3499 & elflexowf_buffer_cart(3,re,alpha,gamma,beta,delta), & 3500 & elflexowf_buffer_cart(3,im,alpha,gamma,beta,delta) 3501 3502 write(74,'(4(i5,3x),2(1x,f20.10))') alpha,gamma,beta,delta, & 3503 & tmpre, tmpim 3504 3505 write(75,'(4(i5,3x),2(1x,f20.10))') alpha,gamma,beta,delta, & 3506 & elflexowf_buffer_cart(5,re,alpha,gamma,beta,delta), & 3507 & elflexowf_buffer_cart(5,im,alpha,gamma,beta,delta) 3508 3509 write(76,'(4(i5,3x),2(1x,f20.10))') alpha,gamma,beta,delta, & 3510 & elflexowf_buffer_cart(4,re,alpha,gamma,beta,delta), & 3511 & elflexowf_buffer_cart(4,im,alpha,gamma,beta,delta) 3512 end if 3513 3514 end if 3515 3516 end do 3517 end do 3518 write(ab_out,'(a)')' ' 3519 if (prtvol>=10) then 3520 write(71,*)' ' 3521 write(72,*)' ' 3522 write(73,*)' ' 3523 write(74,*)' ' 3524 write(75,*)' ' 3525 write(76,*)' ' 3526 end if 3527 end do 3528 end do 3529 3530 if (prtvol>=10) then 3531 close(71) 3532 close(72) 3533 close(73) 3534 close(74) 3535 close(75) 3536 close(76) 3537 ABI_FREE(elflexowf_buffer_cart) 3538 end if 3539 3540 !Calculate the contribution to the d3etot in mixed (reduced/cartesian) coordinates 3541 elec_flexotens_red=elec_flexotens_cart 3542 ABI_FREE(elec_flexotens_cart) 3543 ABI_FREE(elflexowf_t4_cart) 3544 ABI_FREE(cartflg_t4) 3545 ABI_MALLOC(redflg,(3,3,3,3)) 3546 redflg=0 3547 3548 !1st transform back coordinates of the electric field derivative of the flexoelectric tensor 3549 fac=two_pi ** 2 3550 do istr2dir=1,3 3551 do istr1dir=1,3 3552 do iq1dir=1,3 3553 do ii=1,2 3554 do iefidir=1,3 3555 vec1(iefidir)=elec_flexotens_red(ii,iefidir,iq1dir,istr1dir,istr2dir) 3556 flg1(iefidir)=elflexoflg(iefidir,iq1dir,istr1dir,istr2dir) 3557 end do 3558 call cart39(flg1,flg2,transpose(rprimd),matom+2,matom,transpose(gprimd),vec1,vec2) 3559 do iefidir=1,3 3560 elec_flexotens_red(ii,iefidir,iq1dir,istr1dir,istr2dir)=vec2(iefidir)*fac 3561 redflg(iefidir,iq1dir,istr1dir,istr2dir)=flg2(iefidir) 3562 end do 3563 end do 3564 end do 3565 end do 3566 end do 3567 3568 !2nd transform back coordinates of the q-gradient (treat it as electric field) 3569 !of the flexoelectric tensor 3570 do istr2dir=1,3 3571 do istr1dir=1,3 3572 do iefidir=1,3 3573 do ii=1,2 3574 do iq1dir=1,3 3575 vec1(iq1dir)=elec_flexotens_red(ii,iefidir,iq1dir,istr1dir,istr2dir) 3576 flg1(iq1dir)=elflexoflg(iefidir,iq1dir,istr1dir,istr2dir) 3577 end do 3578 call cart39(flg1,flg2,transpose(rprimd),matom+2,matom,transpose(gprimd),vec1,vec2) 3579 do iq1dir=1,3 3580 elec_flexotens_red(ii,iefidir,iq1dir,istr1dir,istr2dir)=vec2(iq1dir)*fac 3581 end do 3582 end do 3583 end do 3584 end do 3585 end do 3586 3587 !Add contributions to d3etot 3588 efipert=matom+2 3589 q1pert=matom+8 3590 fac=ucvol/two 3591 do istrpert=1,nstrpert 3592 strpert=pert_strain(1,istrpert) 3593 strcomp=pert_strain(2,istrpert) 3594 istr1dir=pert_strain(3,istrpert) 3595 istr2dir=pert_strain(4,istrpert) 3596 do iq1grad=1,nq1grad 3597 iq1dir=q1grad(2,iq1grad) 3598 do iefipert=1,nefipert 3599 iefidir=pert_efield(2,iefipert) 3600 3601 if (redflg(iefidir,iq1dir,istr1dir,istr2dir)==1) then 3602 d3etot(re,iefidir,efipert,strcomp,strpert,iq1dir,q1pert)= & 3603 & elec_flexotens_red(im,iefidir,iq1dir,istr1dir,istr2dir)*fac 3604 d3etot(im,iefidir,efipert,strcomp,strpert,iq1dir,q1pert)= & 3605 & -elec_flexotens_red(re,iefidir,iq1dir,istr1dir,istr2dir)*fac 3606 end if 3607 3608 end do 3609 end do 3610 end do 3611 3612 ABI_FREE(elec_flexotens_red) 3613 ABI_FREE(redflg) 3614 3615 3616 DBG_EXIT("COLL") 3617 3618 end subroutine dfpt_ciflexoout
ABINIT/dfpt_ddmdqout [ Functions ]
NAME
dfpt_ddmdqout
FUNCTION
This subroutine gathers the different terms entering the first q derivative of the dynamical matrix, perfofms the transformation from reduced to cartesian coordinates and writes out the tensor in output files.
INPUTS
ddmdq_flg(natpert,natpert,nq1grad)=array that indicates which elements of the first q derivative
of dynamical matrix have been calculated
ddmdq_qgradhart(2,natpert,natpert,nq1grad)=electronic electrostatic contribution from the
q-gradient of the Hartree potential
ddmdqwf(2,natpert,natpert,nq1grad)=total wave function contribution to the first q derivative of dynamical matrix
ddmdqwf_t1(2,natpert,natpert,nq1grad)=term 1 of the wave function contribution
ddmdqwf_t2(2,natpert,natpert,nq1grad)=term 2 of the wave function contribution
ddmdqwf_t3(2,natpert,natpert,nq1grad)=term 3 of the wave function contribution
dyewdq(2,3,natom,3,natom,3)= First q-gradient of Ewald part of the dynamical matrix
gprimd(3,3)=reciprocal space dimensional primitive translations
kptopt=2 time reversal symmetry is enforced, 3 trs is not enforced (for debugging purposes)
matom=number of atoms
mpert=maximum number of perturbations
natpert=number of atomic displacement perturbations
nq1grad=number of q1 (q_{\gamma}) gradients
pert_atdis(3,natpert)=array with the info for the electric field perturbations
prtvol=volume of information to be printed. 1-> The different contributions to the quadrupole are printed.
q1grad(3,nq1grad)=array with the info for the q1 (q_{\gamma}) gradients
rprimd(3,3)=dimensional primitive translations in real space (bohr)
OUTPUT
d3etot(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE
SIDE EFFECTS
NOTES
SOURCE
3661 subroutine dfpt_ddmdqout(ddmdq_flg,ddmdq_qgradhart,ddmdqwf,ddmdqwf_t1,ddmdqwf_t2,ddmdqwf_t3,d3etot, & 3662 & dyewdq,gprimd,kptopt,matom,mpert,natpert,nq1grad,pert_atdis,prtvol,q1grad,rprimd) 3663 3664 !Arguments ------------------------------------ 3665 !scalars 3666 integer,intent(in) :: kptopt,matom,mpert,natpert,nq1grad,prtvol 3667 real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert) 3668 3669 !arrays 3670 integer,intent(in) :: ddmdq_flg(natpert,natpert,nq1grad) 3671 integer,intent(in) :: pert_atdis(3,natpert) 3672 integer,intent(in) :: q1grad(3,nq1grad) 3673 real(dp),intent(inout) :: ddmdq_qgradhart(2,natpert,natpert,nq1grad) 3674 real(dp),intent(in) :: ddmdqwf(2,natpert,natpert,nq1grad) 3675 real(dp),intent(inout) :: ddmdqwf_t1(2,natpert,natpert,nq1grad) 3676 real(dp),intent(inout) :: ddmdqwf_t2(2,natpert,natpert,nq1grad) 3677 real(dp),intent(inout) :: ddmdqwf_t3(2,natpert,natpert,nq1grad) 3678 real(dp),intent(in) :: dyewdq(2,3,matom,3,matom,3) 3679 real(dp),intent(in) :: gprimd(3,3) 3680 real(dp),intent(in) :: rprimd(3,3) 3681 3682 !Local variables------------------------------- 3683 !scalars 3684 integer :: iatdir,iatom,iatpert,ii,iq1dir,iq1grad,iq1pert,jatdir,jatom,jatpert 3685 integer, parameter :: im=2,re=1 3686 real(dp) :: piezofrim,piezofrre,tmpim,tmpre 3687 character(len=500) :: msg 3688 3689 !arrays 3690 integer :: flg1(3),flg2(3) 3691 integer, allocatable :: cartflg(:,:,:,:,:),ddmdq_cartflg(:,:,:,:,:) 3692 real(dp) :: vec1(3),vec2(3) 3693 real(dp), allocatable :: ddmdq_cart(:,:,:,:,:,:),ddmdq_red(:,:,:,:) 3694 3695 !**************************************************************************** 3696 3697 DBG_ENTER("COLL") 3698 3699 !Open output files 3700 if (prtvol>=10) then 3701 open(unit=71,file='ddmdq_wf_t1.out',status='unknown',form='formatted',action='write') 3702 open(unit=72,file='ddmdq_wf_t2.out',status='unknown',form='formatted',action='write') 3703 open(unit=73,file='ddmdq_wf_t3.out',status='unknown',form='formatted',action='write') 3704 open(unit=74,file='ddmdq_ewald.out',status='unknown',form='formatted',action='write') 3705 open(unit=76,file='ddmdq_elecstic.out',status='unknown',form='formatted',action='write') 3706 end if 3707 3708 !Gather the different terms in the tensors and print the result 3709 ABI_MALLOC(ddmdq_red,(2,natpert,natpert,nq1grad)) 3710 3711 if (kptopt==3) then 3712 3713 iq1pert=matom+8 3714 do iq1grad=1,nq1grad 3715 iq1dir=q1grad(2,iq1grad) 3716 do jatpert=1,natpert 3717 jatom=pert_atdis(1,jatpert) 3718 jatdir=pert_atdis(2,jatpert) 3719 do iatpert=1,natpert 3720 iatom=pert_atdis(1,iatpert) 3721 iatdir=pert_atdis(2,iatpert) 3722 3723 if (ddmdq_flg(iatpert,jatpert,iq1grad)==1) then 3724 3725 !Calculate and save the third order energy derivative 3726 tmpre=ddmdq_qgradhart(re,iatpert,jatpert,iq1grad)+ddmdqwf(re,iatpert,jatpert,iq1grad)+& 3727 & half*dyewdq(re,iatdir,iatom,jatdir,jatom,iq1dir) 3728 tmpim=ddmdq_qgradhart(im,iatpert,jatpert,iq1grad)+ddmdqwf(im,iatpert,jatpert,iq1grad)+& 3729 & half*dyewdq(im,iatdir,iatom,jatdir,jatom,iq1dir) 3730 d3etot(re,iatdir,iatom,jatdir,jatom,iq1dir,iq1pert)=tmpre 3731 d3etot(im,iatdir,iatom,jatdir,jatom,iq1dir,iq1pert)=tmpim 3732 3733 !Calculate and write the q-gradient of the dynamical matrix (twice 3734 !the Energy derivative, see Gonze and Lee 1997) in the form of first 3735 !real-space moment of IFCs 3736 ddmdq_red(re,iatpert,jatpert,iq1grad)=-two*tmpim 3737 ddmdq_red(im,iatpert,jatpert,iq1grad)=two*tmpre 3738 3739 if (prtvol>=10) then 3740 !Write individual contributions 3741 write(71,'(5(i5,4x),2(1x,f20.10))') iatom, iatdir, jatom, jatdir, iq1dir, & 3742 & -two*ddmdqwf_t1(:,iatpert,jatpert,iq1grad) 3743 3744 write(72,'(5(i5,4x),2(1x,f20.10))') iatom, iatdir, jatom, jatdir, iq1dir, & 3745 & -two*ddmdqwf_t2(:,iatpert,jatpert,iq1grad) 3746 3747 write(73,'(5(i5,4x),2(1x,f20.10))') iatom, iatdir, jatom, jatdir, iq1dir, & 3748 & -two*ddmdqwf_t3(:,iatpert,jatpert,iq1grad) 3749 3750 write(74,'(5(i5,4x),2(1x,f20.10))') iatom, iatdir, jatom, jatdir, iq1dir, & 3751 & -dyewdq(:,iatdir,iatom,jatdir,jatom,iq1dir) 3752 3753 write(76,'(5(i5,4x),2(1x,f20.10))') iatom, iatdir, jatom, jatdir, iq1dir, & 3754 & -two*ddmdq_qgradhart(:,iatpert,jatpert,iq1grad) 3755 end if 3756 3757 end if 3758 3759 end do 3760 end do 3761 end do 3762 3763 else if (kptopt==2) then 3764 3765 iq1pert=matom+8 3766 do iq1grad=1,nq1grad 3767 iq1dir=q1grad(2,iq1grad) 3768 do jatpert=1,natpert 3769 jatom=pert_atdis(1,jatpert) 3770 jatdir=pert_atdis(2,jatpert) 3771 do iatpert=1,natpert 3772 iatom=pert_atdis(1,iatpert) 3773 iatdir=pert_atdis(2,iatpert) 3774 3775 if (ddmdq_flg(iatpert,jatpert,iq1grad)==1) then 3776 3777 !Calculate and save the third order energy derivative 3778 tmpim=ddmdq_qgradhart(im,iatpert,jatpert,iq1grad)+ddmdqwf(im,iatpert,jatpert,iq1grad)+& 3779 & half*dyewdq(im,iatdir,iatom,jatdir,jatom,iq1dir) 3780 d3etot(re,iatdir,iatom,jatdir,jatom,iq1dir,iq1pert)=zero 3781 d3etot(im,iatdir,iatom,jatdir,jatom,iq1dir,iq1pert)=tmpim 3782 3783 !Calculate and write the q-gradient of the dynamical matrix (twice 3784 !the Energy derivative, see Gonze and Lee 1997) in the form of first 3785 !real-space moment of IFCs 3786 ddmdq_red(re,iatpert,jatpert,iq1grad)=-two*tmpim 3787 ddmdq_red(im,iatpert,jatpert,iq1grad)=zero 3788 3789 if (prtvol>=10) then 3790 !Write individual contributions 3791 write(71,'(5(i5,4x),2(1x,f20.10))') iatom, iatdir, jatom, jatdir, iq1dir, & 3792 & -two*ddmdqwf_t1(im,iatpert,jatpert,iq1grad) 3793 3794 write(72,'(5(i5,4x),2(1x,f20.10))') iatom, iatdir, jatom, jatdir, iq1dir, & 3795 & -two*ddmdqwf_t2(im,iatpert,jatpert,iq1grad) 3796 3797 write(73,'(5(i5,4x),2(1x,f20.10))') iatom, iatdir, jatom, jatdir, iq1dir, & 3798 & -two*ddmdqwf_t3(im,iatpert,jatpert,iq1grad) 3799 3800 write(74,'(5(i5,4x),2(1x,f20.10))') iatom, iatdir, jatom, jatdir, iq1dir, & 3801 & -dyewdq(im,iatdir,iatom,jatdir,jatom,iq1dir) 3802 3803 write(76,'(5(i5,4x),2(1x,f20.10))') iatom, iatdir, jatom, jatdir, iq1dir, & 3804 & -two*ddmdq_qgradhart(im,iatpert,jatpert,iq1grad) 3805 end if 3806 3807 end if 3808 3809 end do 3810 end do 3811 end do 3812 3813 else 3814 3815 write(msg,"(1a)") 'kptopt must be 2 or 3 for the ddmdq calculation' 3816 ABI_BUG(msg) 3817 3818 end if 3819 3820 if (prtvol>=10) then 3821 close(71) 3822 close(72) 3823 close(73) 3824 close(74) 3825 close(76) 3826 end if 3827 3828 !Transformation to cartesian coordinates of the ddmdq 3829 ABI_MALLOC(ddmdq_cart,(2,matom,3,matom,3,3)) 3830 ABI_MALLOC(ddmdq_cartflg,(matom,3,matom,3,3)) 3831 ABI_MALLOC(cartflg,(matom,3,matom,3,3)) 3832 cartflg=0 3833 do iq1grad=1,nq1grad 3834 iq1dir=q1grad(2,iq1grad) 3835 do jatpert=1,natpert 3836 jatom=pert_atdis(1,jatpert) 3837 jatdir=pert_atdis(2,jatpert) 3838 do iatpert=1,natpert 3839 iatom=pert_atdis(1,iatpert) 3840 iatdir=pert_atdis(2,iatpert) 3841 ddmdq_cartflg(iatom,iatdir,jatom,jatdir,iq1dir)=ddmdq_flg(iatpert,jatpert,iq1grad) 3842 ddmdq_cart(:,iatom,iatdir,jatom,jatdir,iq1dir)=ddmdq_red(:,iatpert,jatpert,iq1grad) 3843 end do 3844 end do 3845 end do 3846 ABI_FREE(ddmdq_red) 3847 3848 !1st transform coordenates of the first atomic displacement derivative 3849 do iq1dir=1,3 3850 do jatdir=1,3 3851 do jatom=1,matom 3852 do ii=1,2 3853 do iatom=1,matom 3854 3855 do iatdir=1,3 3856 vec1(iatdir)=ddmdq_cart(ii,iatom,iatdir,jatom,jatdir,iq1dir) 3857 flg1(iatdir)=ddmdq_cartflg(iatom,iatdir,jatom,jatdir,iq1dir) 3858 end do 3859 call cart39(flg1,flg2,gprimd,iatom,matom,rprimd,vec1,vec2) 3860 do iatdir=1,3 3861 ddmdq_cart(ii,iatom,iatdir,jatom,jatdir,iq1dir)=vec2(iatdir) 3862 cartflg(iatom,iatdir,jatom,jatdir,iq1dir)=flg2(iatdir) 3863 end do 3864 3865 end do 3866 end do 3867 end do 3868 end do 3869 end do 3870 3871 !2nd transform coordenates of the second atomic displacement derivative 3872 do iq1dir=1,3 3873 do iatdir=1,3 3874 do iatom=1,matom 3875 do ii=1,2 3876 do jatom=1,matom 3877 3878 do jatdir=1,3 3879 vec1(jatdir)=ddmdq_cart(ii,iatom,iatdir,jatom,jatdir,iq1dir) 3880 flg1(jatdir)=ddmdq_cartflg(iatom,iatdir,jatom,jatdir,iq1dir) 3881 end do 3882 call cart39(flg1,flg2,gprimd,jatom,matom,rprimd,vec1,vec2) 3883 do jatdir=1,3 3884 ddmdq_cart(ii,iatom,iatdir,jatom,jatdir,iq1dir)=vec2(jatdir) 3885 end do 3886 3887 end do 3888 end do 3889 end do 3890 end do 3891 end do 3892 3893 !3rd transform coordinates of the q-gradient (treat it as electric field) 3894 do jatdir=1,3 3895 do jatom=1,matom 3896 do iatdir=1,3 3897 do iatom=1,matom 3898 do ii=1,2 3899 3900 do iq1dir=1,3 3901 vec1(iq1dir)=ddmdq_cart(ii,iatom,iatdir,jatom,jatdir,iq1dir) 3902 flg1(iq1dir)=ddmdq_cartflg(iatom,iatdir,jatom,jatdir,iq1dir) 3903 end do 3904 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 3905 do iq1dir=1,3 3906 ddmdq_cart(ii,iatom,iatdir,jatom,jatdir,iq1dir)=vec2(iq1dir) 3907 end do 3908 3909 end do 3910 end do 3911 end do 3912 end do 3913 end do 3914 3915 !Write the tensor in cartesian coordinates 3916 write(ab_out,'(a)')' ' 3917 write(ab_out,'(a)')' First real-space moment of IFCs, in cartesian coordinates,' 3918 write(ab_out,'(a)')' iatom iatdir jatom jatddir qgrdir real part imaginary part' 3919 do iq1dir=1,3 3920 do jatdir=1,3 3921 do jatom=1,matom 3922 do iatdir=1,3 3923 do iatom=1,matom 3924 3925 if (cartflg(iatom,iatdir,jatom,jatdir,iq1dir)==1) then 3926 write(ab_out,'(5(i5,4x),2(1x,f20.10))') iatom, iatdir, jatom, jatdir, iq1dir, & 3927 & ddmdq_cart(:,iatom,iatdir,jatom,jatdir,iq1dir) 3928 end if 3929 3930 3931 end do 3932 end do 3933 end do 3934 end do 3935 write(ab_out,'(a)')' ' 3936 end do 3937 write(ab_out,'(80a)')('=',ii=1,80) 3938 3939 !Write the piezoelectric force-response tensor 3940 write(ab_out,'(a)')' ' 3941 write(ab_out,'(a)')' Piezoelectric force-response tensor, in cartesian coordinates ' 3942 write(ab_out,'(a)')' (from q-gradient of dynamical matrix),' 3943 write(ab_out,'(a)')' (for non-vanishing forces in the cell it lacks an improper contribution),' 3944 write(ab_out,'(a)')' iatom iatddir jatddir qgrdir real part imaginary part' 3945 do iq1dir=1,3 3946 do iatdir=1,3 3947 do iatom=1,matom 3948 do jatdir=1,3 3949 piezofrre=zero 3950 piezofrim=zero 3951 do jatom=1,matom 3952 3953 if (cartflg(iatom,iatdir,jatom,jatdir,iq1dir)==1) then 3954 piezofrre=piezofrre+ddmdq_cart(1,iatom,iatdir,jatom,jatdir,iq1dir) 3955 piezofrim=piezofrim+ddmdq_cart(2,iatom,iatdir,jatom,jatdir,iq1dir) 3956 end if 3957 3958 end do 3959 3960 write(ab_out,'(4(i5,4x),2(1x,f20.10))') iatom, iatdir, jatdir, iq1dir, & 3961 & piezofrre, piezofrim 3962 3963 end do 3964 end do 3965 end do 3966 write(ab_out,'(a)')' ' 3967 end do 3968 write(ab_out,'(80a)')('=',ii=1,80) 3969 3970 ABI_FREE(ddmdq_cart) 3971 ABI_FREE(ddmdq_cartflg) 3972 ABI_FREE(cartflg) 3973 3974 3975 DBG_EXIT("COLL") 3976 end subroutine dfpt_ddmdqout
ABINIT/dfpt_flexo [ Functions ]
NAME
dfpt_flexo
FUNCTION
This routine computes the elements of the flexoelectric tensor as:
--> Electronic contribution: second q-gradient of the second mixed derivative
of the energy w.r.t an electric field and a metric perturbation.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) codvsn=code version doccde(mband*nkpt*nsppol)=derivative of occupancies wrt the energy dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset dyewdq(2,3,natom,3,natom,3)= First q-gradient of Ewald part of the dynamical matrix dyewdqdq(2,3,natom,3,natom,3,3)= Second q-gradient of Ewald part of the dynamical matrix sumed over the second sublattice gmet(3,3)=reciprocal space metric tensor in bohr**-2. gprimd(3,3)=reciprocal space dimensional primitive translations kxc(nfft,nkxc)=exchange and correlation kernel mpert=maximum number of perturbations for output processing mpi_enreg=information about MPI parallelization nattyp(ntypat)= # atoms of each type. nfft=(effective) number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other nkpt=number of k points in the full BZ nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed. nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS (DUMMY) pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data (DUMMY) pertsy(3,natom+6)=set of perturbations that form a basis for all other perturbations psps <type(pseudopotential_type)>=variables related to pseudopotentials rmet(3,3)=real space metric (bohr**2) rprimd(3,3)=dimensional primitive translations in real space (bohr) rhog(2,nfftf)=array for Fourier transform of GS electron density rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3. timrev=1 if time-reversal preserves the q wavevector; 0 otherwise. ucvol=unit cell volume in bohr**3. xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
blkflg(3,mpert,3,mpert,3,mpert)= ( 1 if the element of the 3dte has been calculated ; 0 otherwise ) d3etot(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE
SIDE EFFECTS
NOTES
SOURCE
1615 subroutine dfpt_flexo(atindx,blkflg,codvsn,d3etot,doccde,dtfil,dtset,dyewdq,dyewdqdq,& 1616 & gmet,gprimd,kxc,mpert,& 1617 & mpi_enreg,nattyp,nfft,ngfft,nkpt,nkxc,& 1618 & nspden,nsppol,occ,pawrhoij,pawtab,pertsy,psps,rmet,rprimd,rhog,rhor,& 1619 & timrev,ucvol,xred) 1620 1621 !Arguments ------------------------------------ 1622 !scalars 1623 integer,intent(in) :: mpert,nfft,nkpt,nkxc,nspden,nsppol,timrev 1624 real(dp),intent(in) :: ucvol 1625 character(len=8), intent(in) :: codvsn 1626 type(MPI_type),intent(inout) :: mpi_enreg 1627 type(datafiles_type),intent(in) :: dtfil 1628 type(dataset_type),intent(in) :: dtset 1629 type(hdr_type) :: hdr_den 1630 type(pseudopotential_type),intent(in) :: psps 1631 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 1632 type(pawrhoij_type),intent(inout) :: pawrhoij(:) 1633 !arrays 1634 integer,intent(in) :: atindx(dtset%natom) 1635 integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) 1636 integer,intent(in) :: nattyp(dtset%ntypat),ngfft(18) 1637 integer,intent(in) :: pertsy(3,dtset%natom+6) 1638 real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert) 1639 real(dp),intent(in) :: doccde(dtset%mband*nkpt*dtset%nsppol) 1640 real(dp),intent(in) :: dyewdq(2,3,dtset%natom,3,dtset%natom,3) 1641 real(dp),intent(inout) :: dyewdqdq(2,3,dtset%natom,3,3,3) 1642 real(dp),intent(in) :: gmet(3,3), gprimd(3,3) 1643 real(dp),intent(in) :: kxc(nfft,nkxc) 1644 real(dp),intent(in) :: occ(dtset%mband*nkpt*dtset%nsppol) 1645 real(dp),intent(in) :: rmet(3,3),rprimd(3,3) 1646 real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,nspden) 1647 real(dp),intent(in) :: xred(3,dtset%natom) 1648 1649 !Local variables------------------------------- 1650 !scalars 1651 integer :: ask_accurate,bantot_rbz,bdtot_index,bdtot1_index 1652 integer :: cplex,formeig,forunit,gscase,icg 1653 integer :: ia1,iatdir,iatom,iatpert,iatpert_cnt,iatpol 1654 integer :: ii,iefipert,iefipert_cnt,ierr,ikg,ikpt,ikpt1,ilm 1655 integer :: iq1dir,iq2dir,iq1grad,iq1grad_cnt,iq1q2grad,iq1q2grad_var 1656 integer :: ireadwf0,isppol,istrdir,istrpert,istrtype,istrpert_cnt,istwf_k,itypat,jatpert,jj,ka,kb 1657 integer :: lw_flexo,master,matom,matpert,mcg,me,mgfft,mkmem_rbz,mpw,my_nkpt_rbz 1658 integer :: natpert,nband_k,nefipert,nfftot,nhat1grdim,nkpt_rbz 1659 integer :: npw_k,npw1_k,nq1grad,nq1q2grad,nstrpert,nsym1,n3xccc 1660 integer :: nylmgr,optene,option,optorth,optres 1661 integer :: pawread,pertcase,qcar,qdir,spaceworld 1662 integer :: usexcnhat,useylmgr 1663 integer :: mband_mem 1664 integer,parameter :: formeig1=1 1665 integer,parameter :: re=1,im=2 1666 real(dp) :: boxcut,delad,delag,delbd,delbg 1667 real(dp) :: doti,dotr,dum_scl,ecut_eff,ecut,etotal,fac,fermie,gsqcut,residm 1668 real(dp) :: vres2,wtk_k 1669 logical :: non_magnetic_xc 1670 character(len=500) :: msg 1671 character(len=fnlen) :: fi1o,fiwfatdis,fiwfstrain,fiwfefield,fiwfddk,fiwfdkdk 1672 type(ebands_t) :: bs_rbz 1673 type(hdr_type) :: hdr0 1674 type(wvl_data) :: wvl 1675 type(wffile_type) :: wffgs,wfftgs 1676 type(gs_hamiltonian_type) :: gs_hamkq 1677 1678 !arrays 1679 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 1680 integer,allocatable :: bz2ibz_smap(:,:) 1681 integer,allocatable :: ddmdq_flg(:,:,:,:,:),elflexoflg(:,:,:,:) 1682 integer,allocatable :: indkpt1(:), indkpt1_tmp(:) 1683 integer,allocatable :: indsy1(:,:,:),irrzon1(:,:,:) 1684 integer,allocatable :: istwfk_rbz(:),isdq_flg(:,:,:,:,:) 1685 integer,allocatable :: kg(:,:),kg_k(:,:) 1686 integer,allocatable :: nband_rbz(:),npwarr(:),npwtot(:) 1687 integer,allocatable :: pert_atdis(:,:),pert_atdis_tmp(:,:) 1688 integer,allocatable :: pert_efield(:,:),pert_efield_tmp(:,:) 1689 integer,allocatable :: pert_strain(:,:),pert_strain_tmp(:,:) 1690 integer,allocatable :: q1grad(:,:),q1grad_tmp(:,:),q1q2grad(:,:) 1691 integer,allocatable :: symaf1(:),symrc1(:,:,:),symrl1(:,:,:) 1692 real(dp) :: kpoint(3) 1693 real(dp),allocatable :: cg(:,:),doccde_rbz(:) 1694 real(dp),allocatable :: ddmdq_qgradhart(:,:,:,:) 1695 real(dp),allocatable :: ddmdqwf(:,:,:,:),ddmdqwf_k(:,:,:,:) 1696 real(dp),allocatable :: ddmdqwf_t1(:,:,:,:),ddmdqwf_t1_k(:,:,:,:) 1697 real(dp),allocatable :: ddmdqwf_t2(:,:,:,:),ddmdqwf_t2_k(:,:,:,:) 1698 real(dp),allocatable :: ddmdqwf_t3(:,:,:,:),ddmdqwf_t3_k(:,:,:,:) 1699 real(dp),allocatable :: eigen0(:) 1700 real(dp),allocatable :: elflexowf(:,:,:,:,:),elflexowf_k(:,:,:,:,:) 1701 real(dp),allocatable :: elflexowf_t1(:,:,:,:,:),elflexowf_t1_k(:,:,:,:,:) 1702 real(dp),allocatable :: elflexowf_t2(:,:,:,:,:),elflexowf_t2_k(:,:,:,:,:) 1703 real(dp),allocatable :: elflexowf_t3(:,:,:,:,:),elflexowf_t3_k(:,:,:,:,:) 1704 real(dp),allocatable :: elflexowf_t4(:,:,:,:,:),elflexowf_t4_k(:,:,:,:,:) 1705 real(dp),allocatable :: elflexowf_t5(:,:,:,:,:),elflexowf_t5_k(:,:,:,:,:) 1706 real(dp),allocatable :: elqgradhart(:,:,:,:,:) 1707 real(dp),allocatable :: frwfdq(:,:,:,:,:,:),frwfdq_k(:,:,:,:,:,:) 1708 real(dp),allocatable :: isdq_qgradhart(:,:,:,:,:,:) 1709 real(dp),allocatable :: isdqwf(:,:,:,:,:,:),isdqwf_k(:,:,:,:,:,:) 1710 real(dp),allocatable :: isdqwf_t1(:,:,:,:,:,:),isdqwf_t1_k(:,:,:,:,:,:) 1711 real(dp),allocatable :: isdqwf_t2(:,:,:,:,:,:),isdqwf_t2_k(:,:,:,:,:,:) 1712 real(dp),allocatable :: isdqwf_t3(:,:,:,:,:,:),isdqwf_t3_k(:,:,:,:,:,:) 1713 real(dp),allocatable :: isdqwf_t4(:,:,:,:,:,:),isdqwf_t4_k(:,:,:,:,:,:) 1714 real(dp),allocatable :: isdqwf_t5(:,:,:,:,:,:),isdqwf_t5_k(:,:,:,:,:,:) 1715 real(dp),allocatable :: kpt_rbz(:,:) 1716 real(dp),allocatable :: nhat(:,:),nhat1(:,:),nhat1gr(:,:,:) 1717 real(dp),allocatable :: occ_k(:),occ_rbz(:) 1718 real(dp),allocatable :: ph1d(:,:),phnons1(:,:,:) 1719 real(dp),allocatable :: rhog1_tmp(:,:),rhog1_atdis(:,:,:) 1720 real(dp),allocatable :: rhog1_efield(:,:,:),rhor1_atdis(:,:,:),rhor1_efield(:,:,:) 1721 real(dp),allocatable :: rhor1_cplx(:,:),rhor1_tmp(:,:),rhor1_real(:,:) 1722 real(dp),allocatable :: rhor1_strain(:,:,:) 1723 real(dp),allocatable :: vhartr1(:),vhxc1_atdis(:,:),vhxc1_efield(:,:),vhxc1_strain(:,:) 1724 real(dp),allocatable :: vpsp1(:),vqgradhart(:),vresid1(:,:) 1725 real(dp),allocatable :: vxc1(:,:),vxc1dq(:,:),vxc1dqc(:,:,:,:) 1726 real(dp),allocatable :: dum_vxc(:,:) 1727 real(dp),allocatable :: wtk_folded(:), wtk_rbz(:) 1728 real(dp),allocatable,target :: vtrial1(:,:) 1729 real(dp),allocatable :: tnons1(:,:) 1730 real(dp),allocatable :: xccc3d1(:) 1731 real(dp),allocatable :: ylm(:,:),ylm_k(:,:),ylmgr(:,:,:),ylmgr_k(:,:,:) 1732 type(pawrhoij_type),allocatable :: pawrhoij_read(:) 1733 type(wfk_t),allocatable :: wfk_t_atdis(:),wfk_t_efield(:),wfk_t_ddk(:) 1734 type(wfk_t),allocatable :: wfk_t_dkdk(:),wfk_t_strain(:,:) 1735 ! ************************************************************************* 1736 1737 DBG_ENTER("COLL") 1738 1739 !Anounce start of flexoelectric tensor calculation 1740 write(msg, '(a,80a,a,a,a)' ) ch10,('=',ii=1,80),ch10,& 1741 & ' ==> Compute Flexoelectric Tensor Related Magnitudes <== ',ch10 1742 call wrtout(std_out,msg,'COLL') 1743 call wrtout(ab_out,msg,'COLL') 1744 1745 !Init parallelism 1746 spaceworld=mpi_enreg%comm_cell 1747 me=mpi_enreg%me_kpt 1748 master=0 1749 1750 !Get FFT grid(s) sizes (be careful !) See NOTES in the comments at the beginning of respfn.F90 1751 mgfft=dtset%mgfft 1752 ecut=dtset%ecut 1753 ecut_eff=ecut*(dtset%dilatmx)**2 1754 1755 !Compute large sphere cut-off gsqcut 1756 call getcut(boxcut,ecut,gmet,gsqcut,dtset%iboxcut,std_out,dtset%qptn,ngfft) 1757 1758 !Various initializations 1759 lw_flexo=dtset%lw_flexo 1760 cplex=2-timrev 1761 matom=dtset%natom 1762 usexcnhat=0 1763 n3xccc=0 1764 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 1765 non_magnetic_xc=.true. 1766 1767 !Generate the 1-dimensional phases 1768 ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom)) 1769 call getph(atindx,dtset%natom,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),ph1d,xred) 1770 1771 !################# PERTURBATIONS AND q-GRADIENTS LABELLING ############################ 1772 1773 !Determine which atomic displacement, electric field, strain and q-gradient directions 1774 !have to be evaluated taking into account the perturbation symmetries 1775 matpert=dtset%natom*3 1776 1777 !Atomic displacement 1778 if (lw_flexo==1.or.lw_flexo==3.or.lw_flexo==4) then 1779 ABI_MALLOC(pert_atdis_tmp,(3,matpert)) 1780 pert_atdis_tmp=0 1781 iatpert_cnt=0 1782 do iatpol=1,matom 1783 do iatdir=1,3 1784 if (pertsy(iatdir,iatpol)==1) then 1785 iatpert_cnt=iatpert_cnt+1 1786 pert_atdis_tmp(1,iatpert_cnt)=iatpol !atom displaced 1787 pert_atdis_tmp(2,iatpert_cnt)=iatdir !direction of displacement 1788 pert_atdis_tmp(3,iatpert_cnt)=(iatpol-1)*3+iatdir !like pertcase in dfpt_loopert.f90 1789 end if 1790 end do 1791 end do 1792 natpert=iatpert_cnt 1793 ABI_MALLOC(pert_atdis,(3,natpert)) 1794 do iatpert=1,natpert 1795 pert_atdis(:,iatpert)=pert_atdis_tmp(:,iatpert) 1796 end do 1797 ABI_FREE(pert_atdis_tmp) 1798 end if 1799 1800 !Electric field 1801 if (lw_flexo==1.or.lw_flexo==2) then 1802 ABI_MALLOC(pert_efield_tmp,(3,3)) 1803 pert_efield_tmp=0 1804 iefipert_cnt=0 1805 do iefipert=1,3 1806 if (pertsy(iefipert,dtset%natom+2)==1) then 1807 iefipert_cnt=iefipert_cnt+1 1808 pert_efield_tmp(1,iefipert_cnt)=dtset%natom+2 !electric field pert 1809 pert_efield_tmp(2,iefipert_cnt)=iefipert !electric field direction 1810 pert_efield_tmp(3,iefipert_cnt)=matpert+3+iefipert !like pertcase in dfpt_loopert.f90 1811 end if 1812 end do 1813 nefipert=iefipert_cnt 1814 ABI_MALLOC(pert_efield,(3,nefipert)) 1815 do iefipert=1,nefipert 1816 pert_efield(:,iefipert)=pert_efield_tmp(:,iefipert) 1817 end do 1818 ABI_FREE(pert_efield_tmp) 1819 end if 1820 1821 !ddk 1822 !The q1grad is related with the response to the ddk 1823 ABI_MALLOC(q1grad_tmp,(3,3)) 1824 q1grad_tmp=0 1825 iq1grad_cnt=0 1826 do iq1grad=1,3 1827 if (pertsy(iq1grad,dtset%natom+1)==1) then 1828 iq1grad_cnt=iq1grad_cnt+1 1829 q1grad_tmp(1,iq1grad_cnt)=dtset%natom+1 !ddk perturbation 1830 q1grad_tmp(2,iq1grad_cnt)=iq1grad !ddk direction 1831 q1grad_tmp(3,iq1grad_cnt)=matpert+iq1grad !like pertcase in dfpt_loopert.f90 1832 end if 1833 end do 1834 nq1grad=iq1grad_cnt 1835 ABI_MALLOC(q1grad,(3,nq1grad)) 1836 do iq1grad=1,nq1grad 1837 q1grad(:,iq1grad)=q1grad_tmp(:,iq1grad) 1838 end do 1839 ABI_FREE(q1grad_tmp) 1840 1841 !d2_dkdk 1842 !For the evaluation of the 2nd order q-gradient, the 9 directios are activated because 1843 !currently the program calculates by defect all the components of the d2_dkdk perturbation. 1844 !TODO: This will have to be modified in the future when ABINIT enables to calculate specific 1845 !components of the d2_dkdk 1846 if (lw_flexo==1.or.lw_flexo==2) then 1847 nq1q2grad=9 1848 ABI_MALLOC(q1q2grad,(4,nq1q2grad)) 1849 do iq1q2grad=1,nq1q2grad 1850 call rf2_getidirs(iq1q2grad,iq1dir,iq2dir) 1851 if (iq1dir==iq2dir) then 1852 q1q2grad(1,iq1q2grad)=dtset%natom+10 !d2_dkdk perturbation diagonal elements 1853 else 1854 q1q2grad(1,iq1q2grad)=dtset%natom+11 !d2_dkdk perturbation off-diagonal elements 1855 end if 1856 q1q2grad(2,iq1q2grad)=iq1dir !dq1 direction 1857 q1q2grad(3,iq1q2grad)=iq2dir !dq2 direction 1858 iq1q2grad_var=iq1q2grad 1859 if (iq1q2grad>6) iq1q2grad_var=iq1q2grad-3 !Lower=Upper diagonal triangle matrix 1860 q1q2grad(4,iq1q2grad)=iq1q2grad_var+(dtset%natom+6)*3 !like pertcase in dfpt_loopert.f90 1861 end do 1862 end if 1863 1864 !Strain perturbation 1865 if (lw_flexo==1.or.lw_flexo==2.or.lw_flexo==4) then 1866 ABI_MALLOC(pert_strain_tmp,(6,9)) 1867 pert_strain_tmp=0 1868 !tmp uniaxial components 1869 istrpert_cnt=0 1870 do istrdir=1,3 1871 if (pertsy(istrdir,dtset%natom+3)==1) then 1872 istrpert_cnt=istrpert_cnt+1 1873 ka=idx(2*istrdir-1);kb=idx(2*istrdir) 1874 pert_strain_tmp(1,istrpert_cnt)=dtset%natom+3 !Uniaxial strain perturbation 1875 pert_strain_tmp(2,istrpert_cnt)=istrdir !Uniaxial strain case 1876 pert_strain_tmp(3,istrpert_cnt)=ka !Strain direction 1 1877 pert_strain_tmp(4,istrpert_cnt)=kb !Strain direction 2 1878 pert_strain_tmp(5,istrpert_cnt)=matpert+6+istrdir !like pertcase in dfpt_loopert.f90 1879 pert_strain_tmp(6,istrpert_cnt)=istrdir !Indexing for the second q-gradient of the metric Hamiltonian 1880 end if 1881 end do 1882 !tmp shear components 1883 do istrdir=1,3 1884 if (pertsy(istrdir,dtset%natom+4)==1) then 1885 istrpert_cnt=istrpert_cnt+1 1886 ka=idx(2*(istrdir+3)-1);kb=idx(2*(istrdir+3)) 1887 pert_strain_tmp(1,istrpert_cnt)=dtset%natom+4 !Shear strain perturbation 1888 pert_strain_tmp(2,istrpert_cnt)=istrdir !Shear strain case 1889 pert_strain_tmp(3,istrpert_cnt)=ka !Strain direction 1 1890 pert_strain_tmp(4,istrpert_cnt)=kb !Strain direction 2 1891 pert_strain_tmp(5,istrpert_cnt)=matpert+9+istrdir !like pertcase in dfpt_loopert.f90 1892 pert_strain_tmp(6,istrpert_cnt)=3+istrdir !Indexing for the second q-gradient of the metric Hamiltonian 1893 end if 1894 end do 1895 do istrdir=1,3 1896 if (pertsy(istrdir,dtset%natom+4)==1) then 1897 istrpert_cnt=istrpert_cnt+1 1898 ka=idx(2*(istrdir+3)-1);kb=idx(2*(istrdir+3)) 1899 pert_strain_tmp(1,istrpert_cnt)=dtset%natom+4 !Shear strain perturbation 1900 pert_strain_tmp(2,istrpert_cnt)=istrdir !Shear strain case 1901 pert_strain_tmp(3,istrpert_cnt)=kb !Strain direction 1 1902 pert_strain_tmp(4,istrpert_cnt)=ka !Strain direction 2 1903 pert_strain_tmp(5,istrpert_cnt)=matpert+9+istrdir !like pertcase in dfpt_loopert.f90 1904 pert_strain_tmp(6,istrpert_cnt)=6+istrdir !Indexing for the second q-gradient of the metric Hamiltonian 1905 end if 1906 end do 1907 nstrpert=istrpert_cnt 1908 ABI_MALLOC(pert_strain,(6,nstrpert)) 1909 do istrpert=1,nstrpert 1910 pert_strain(:,istrpert)=pert_strain_tmp(:,istrpert) 1911 end do 1912 ABI_FREE(pert_strain_tmp) 1913 end if 1914 1915 !################# ELECTROSTATIC CONTRIBUTIONS ####################################### 1916 1917 !This is necessary to deactivate paw options in the dfpt_rhotov routine 1918 ABI_MALLOC(pawrhoij_read,(0)) 1919 pawread=0 1920 nhat1grdim=0 1921 ABI_MALLOC(nhat1gr,(0,0,0)) 1922 ABI_MALLOC(nhat,(nfft,nspden)) 1923 nhat=zero 1924 ABI_MALLOC(nhat1,(cplex*nfft,nspden)) 1925 nhat1=zero 1926 1927 !Read the first order densities response from a disk file, calculates the FFT 1928 !(rhog1_tmp) and the first order Hartree and xc potentials(vhxc1_pert). 1929 !TODO: In the call to read_rhor there is a security option that compares with the header 1930 !hdr. Not activated at this moment. 1931 ABI_MALLOC(rhog1_tmp,(2,nfft)) 1932 ABI_MALLOC(rhor1_tmp,(cplex*nfft,nspden)) 1933 ABI_MALLOC(rhor1_real,(1*nfft,nspden)) 1934 ABI_MALLOC(vhartr1,(cplex*nfft)) 1935 ABI_MALLOC(vpsp1,(cplex*nfft)) 1936 ABI_MALLOC(vtrial1,(cplex*nfft,nspden)) 1937 ABI_MALLOC(vresid1,(cplex*nfft,nspden)) 1938 ABI_MALLOC(vxc1,(cplex*nfft,nspden)) 1939 ABI_MALLOC(dum_vxc,(nfft,nspden)) 1940 ABI_MALLOC(xccc3d1,(cplex*n3xccc)) 1941 vpsp1=zero; dum_vxc=zero 1942 optene=0; optres=1 1943 1944 !Atomic displacement 1945 if (lw_flexo==1.or.lw_flexo==3.or.lw_flexo==4) then 1946 ABI_MALLOC(rhor1_atdis,(natpert,cplex*nfft,nspden)) 1947 ABI_MALLOC(rhog1_atdis,(natpert,2,nfft)) 1948 ABI_MALLOC(vhxc1_atdis,(natpert,cplex*nfft)) 1949 vtrial1=zero 1950 do iatpert= 1, natpert 1951 iatpol=pert_atdis(1,iatpert) 1952 iatdir=pert_atdis(2,iatpert) 1953 pertcase=pert_atdis(3,iatpert) 1954 1955 !Reads a real first order density 1956 call appdig(pertcase,dtfil%fildens1in,fi1o) 1957 call read_rhor(fi1o, 1, nspden, nfft, ngfft, pawread, mpi_enreg, rhor1_real, & 1958 & hdr_den, pawrhoij_read, spaceworld) 1959 call hdr_den%free() 1960 1961 !Perform FFT rhor1 to rhog1 1962 call fourdp(cplex,rhog1_tmp,rhor1_real,-1,mpi_enreg,nfft,1,ngfft,0) 1963 1964 !Accumulate density in meaningful complex arrays 1965 if (timrev==0) then 1966 do ii=1,nfft 1967 jj=ii*2 1968 rhor1_tmp(jj-1,:)=rhor1_real(ii,:) 1969 rhor1_tmp(jj,:)=zero 1970 end do 1971 else if (timrev==1) then 1972 rhor1_tmp(:,:)=rhor1_real(:,:) 1973 end if 1974 rhog1_atdis(iatpert,:,:)=rhog1_tmp(:,:) 1975 rhor1_atdis(iatpert,:,:)=rhor1_tmp(:,:) 1976 1977 !Calculate first order Hartree and xc potentials 1978 call dfpt_rhotov(cplex,dum_scl,dum_scl,dum_scl,dum_scl,dum_scl, & 1979 & gsqcut,iatdir,iatpol,& 1980 & dtset%ixc,kxc,mpi_enreg,dtset%natom,nfft,ngfft,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,& 1981 & nspden,n3xccc,non_magnetic_xc,optene,optres,dtset%qptn,rhog,rhog1_tmp,rhor,rhor1_tmp,& 1982 & rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1,vpsp1,vresid1,vres2,vtrial1,dum_vxc,vxc1,& 1983 & xccc3d1,dtset%ixcrot) 1984 1985 !Accumulate the potential in meaningful arrays 1986 vhxc1_atdis(iatpert,:)=vtrial1(:,nspden) 1987 1988 end do 1989 end if 1990 1991 !Electric field 1992 if (lw_flexo==1.or.lw_flexo==2) then 1993 ABI_MALLOC(rhog1_efield,(nefipert,2,nfft)) 1994 ABI_MALLOC(rhor1_efield,(nefipert,cplex*nfft,nspden)) 1995 ABI_MALLOC(vhxc1_efield,(nefipert,cplex*nfft)) 1996 vtrial1=zero 1997 do iefipert=1,nefipert 1998 pertcase=pert_efield(3,iefipert) 1999 2000 !Reads a real first order density 2001 call appdig(pertcase,dtfil%fildens1in,fi1o) 2002 call read_rhor(fi1o, 1, nspden, nfft, ngfft, pawread, mpi_enreg, rhor1_real, & 2003 & hdr_den, pawrhoij_read, spaceworld) 2004 call hdr_den%free() 2005 2006 !Perform FFT rhor1 to rhog1 2007 call fourdp(cplex,rhog1_tmp,rhor1_real,-1,mpi_enreg,nfft,1,ngfft,0) 2008 2009 !Accumulate density in meaningful complex arrays 2010 if (timrev==0) then 2011 do ii=1,nfft 2012 jj=ii*2 2013 rhor1_tmp(jj-1,:)=rhor1_real(ii,:) 2014 rhor1_tmp(jj,:)=zero 2015 end do 2016 else if (timrev==1) then 2017 rhor1_tmp(:,:)=rhor1_real(:,:) 2018 end if 2019 rhog1_efield(iefipert,:,:)=rhog1_tmp(:,:) 2020 rhor1_efield(iefipert,:,:)=rhor1_tmp(:,:) 2021 2022 !Calculate first order Hartree and xc potentials 2023 call dfpt_rhotov(cplex,dum_scl,dum_scl,dum_scl,dum_scl,dum_scl, & 2024 & gsqcut,pert_efield(2,iefipert),dtset%natom+2,& 2025 & dtset%ixc,kxc,mpi_enreg,dtset%natom,nfft,ngfft,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,& 2026 & nspden,n3xccc,non_magnetic_xc,optene,optres,dtset%qptn,rhog,rhog1_tmp,rhor,rhor1_tmp,& 2027 & rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1,vpsp1,vresid1,vres2,vtrial1,dum_vxc,vxc1,& 2028 & xccc3d1,dtset%ixcrot) 2029 2030 !Accumulate the potential in meaningful arrays 2031 vhxc1_efield(iefipert,:)=vtrial1(:,nspden) 2032 2033 end do 2034 endif 2035 2036 !Strain 2037 if (lw_flexo==1.or.lw_flexo==2.or.lw_flexo==4) then 2038 ABI_MALLOC(rhor1_strain,(nstrpert,cplex*nfft,nspden)) 2039 ABI_MALLOC(vhxc1_strain,(nstrpert,cplex*nfft)) 2040 vtrial1=zero 2041 do istrpert= 1, nstrpert 2042 istrtype=pert_strain(1,istrpert) 2043 istrdir=pert_strain(2,istrpert) 2044 pertcase=pert_strain(5,istrpert) 2045 2046 !Reads a real first order density 2047 call appdig(pertcase,dtfil%fildens1in,fi1o) 2048 call read_rhor(fi1o, 1, nspden, nfft, ngfft, pawread, mpi_enreg, rhor1_real, & 2049 & hdr_den, pawrhoij_read, spaceworld) 2050 call hdr_den%free() 2051 2052 !Perform FFT rhor1 to rhog1 2053 call fourdp(cplex,rhog1_tmp,rhor1_real,-1,mpi_enreg,nfft,1,ngfft,0) 2054 2055 !Accumulate density in meaningful complex arrays 2056 if (timrev==0) then 2057 do ii=1,nfft 2058 jj=ii*2 2059 rhor1_tmp(jj-1,:)=rhor1_real(ii,:) 2060 rhor1_tmp(jj,:)=zero 2061 end do 2062 else if (timrev==1) then 2063 rhor1_tmp(:,:)=rhor1_real(:,:) 2064 end if 2065 rhor1_strain(istrpert,:,:)=rhor1_tmp(:,:) 2066 2067 !Calculate first order Hartree and xc potentials 2068 call dfpt_rhotov(cplex,dum_scl,dum_scl,dum_scl,dum_scl,dum_scl, & 2069 & gsqcut,istrdir,istrtype,& 2070 & dtset%ixc,kxc,mpi_enreg,dtset%natom,nfft,ngfft,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,& 2071 & nspden,n3xccc,non_magnetic_xc,optene,optres,dtset%qptn,rhog,rhog1_tmp,rhor,rhor1_tmp,& 2072 & rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1,vpsp1,vresid1,vres2,vtrial1,dum_vxc,vxc1,& 2073 & xccc3d1,dtset%ixcrot) 2074 2075 !Accumulate the potential in meaningful arrays 2076 vhxc1_strain(istrpert,:)=vtrial1(:,nspden) 2077 2078 end do 2079 end if 2080 2081 !These arrays will not be used anymore (for the moment) 2082 ABI_FREE(rhor1_real) 2083 ABI_FREE(rhor1_tmp) 2084 ABI_FREE(vhartr1) 2085 ABI_FREE(vpsp1) 2086 ABI_FREE(vtrial1) 2087 ABI_FREE(vresid1) 2088 ABI_FREE(vxc1) 2089 ABI_FREE(dum_vxc) 2090 ABI_FREE(xccc3d1) 2091 2092 ABI_FREE(pawrhoij_read) 2093 ABI_FREE(nhat1gr) 2094 ABI_FREE(nhat) 2095 ABI_FREE(nhat1) 2096 2097 !!Calculate the electrostatic term from the q-gradient of the Hxc kernel 2098 ABI_MALLOC(vqgradhart,(2*nfft)) 2099 ABI_MALLOC(rhor1_cplx,(2*nfft,nspden)) 2100 rhor1_cplx=zero 2101 if (nkxc == 7) then 2102 ABI_MALLOC(vxc1dq,(2*nfft,nspden)) 2103 ABI_MALLOC(rhor1_tmp,(cplex*nfft,nspden)) 2104 end if 2105 2106 !Electronic contribution 2107 if (lw_flexo==1.or.lw_flexo==2) then 2108 !If GGA xc first calculate the Cartesian q-gradient of the xc kernel 2109 if (nkxc == 7) then 2110 ABI_MALLOC(vxc1dqc,(2*nfft,nspden,nefipert,3)) 2111 do qcar=1,3 2112 do iefipert=1,nefipert 2113 rhor1_tmp(:,:)=rhor1_efield(iefipert,:,:) 2114 call dfpt_mkvxcggadq(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,qcar,rhor1_tmp,vxc1dq) 2115 vxc1dqc(:,:,iefipert,qcar)=vxc1dq(:,:) 2116 end do 2117 end do 2118 end if 2119 2120 ABI_MALLOC(elqgradhart,(2,3,3,3,3)) 2121 ABI_MALLOC(elflexoflg,(3,3,3,3)) 2122 elflexoflg=0 2123 do iq1grad=1,nq1grad 2124 qdir=q1grad(2,iq1grad) 2125 do iefipert=1,nefipert 2126 2127 !Calculate the gradient of the potential generated by the first order electric field density 2128 rhog1_tmp(:,:)=rhog1_efield(iefipert,:,:) 2129 call hartredq(2,gmet,gsqcut,mpi_enreg,nfft,ngfft,qdir,rhog1_tmp,vqgradhart) 2130 2131 !To ckeck 2132 !call appdig(pert_efield(3,iefipert)+q1grad(3,iq1grad),"Gradient_Hartree_potential",fi1o) 2133 !call fftdatar_write_from_hdr("first_order_potential",fi1o,dtset%iomode,hdr_den,& 2134 ! & ngfft,cplex,nfft,nspden,vqgradhart,mpi_enreg) 2135 2136 !If GGA convert the gradient of xc kernel to reduced coordinates and incorporate it to the Hartree part 2137 if (nkxc == 7) then 2138 vxc1dq=zero 2139 do qcar=1,3 2140 vxc1dq(:,:)=vxc1dq(:,:) + gprimd(qcar,qdir) * & 2141 & vxc1dqc(:,:,iefipert,qcar) 2142 end do 2143 vqgradhart(:)=vqgradhart(:)+vxc1dq(:,1) 2144 end if 2145 2146 do istrpert=1,nstrpert 2147 2148 !Calculate the electrostatic energy term with the first order strain density 2149 !I need a cplex density for the dotprod_vn 2150 do ii=1,nfft 2151 jj=ii*2 2152 rhor1_cplx(jj-1,:)=rhor1_strain(istrpert,ii,:) 2153 end do 2154 2155 call dotprod_vn(2,rhor1_cplx,dotr,doti,nfft,nfftot,nspden,2,vqgradhart,ucvol) 2156 elqgradhart(re,pert_efield(2,iefipert),q1grad(2,iq1grad),pert_strain(3,istrpert),pert_strain(4,istrpert))=dotr*half 2157 elqgradhart(im,pert_efield(2,iefipert),q1grad(2,iq1grad),pert_strain(3,istrpert),pert_strain(4,istrpert))=doti*half 2158 elflexoflg(pert_efield(2,iefipert),q1grad(2,iq1grad),pert_strain(3,istrpert),pert_strain(4,istrpert))=1 2159 2160 blkflg(pert_efield(2,iefipert),pert_efield(1,iefipert), & 2161 & pert_strain(2,istrpert),pert_strain(1,istrpert), & 2162 & q1grad(2,iq1grad),matom+8)=1 2163 2164 end do 2165 end do 2166 end do 2167 if (nkxc == 7) then 2168 ABI_FREE(vxc1dqc) 2169 end if 2170 end if 2171 2172 !Calculate here the Cartesian q-gradient of the GGA xc kernel which is the same 2173 !for the other two spatial-dispersion properties 2174 if (lw_flexo==1.or.lw_flexo==3.or.lw_flexo==4) then 2175 if (nkxc == 7) then 2176 ABI_MALLOC(vxc1dqc,(2*nfft,nspden,natpert,3)) 2177 do qcar=1,3 2178 do iatpert=1,natpert 2179 rhor1_tmp(:,:)=rhor1_atdis(iatpert,:,:) 2180 call dfpt_mkvxcggadq(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,qcar,rhor1_tmp,vxc1dq) 2181 vxc1dqc(:,:,iatpert,qcar)=vxc1dq(:,:) 2182 end do 2183 end do 2184 ABI_FREE(rhor1_tmp) 2185 end if 2186 end if 2187 2188 !1st q-gradient of DM contribution 2189 if (lw_flexo==1.or.lw_flexo==3) then 2190 ABI_MALLOC(ddmdq_qgradhart,(2,natpert,natpert,nq1grad)) 2191 ABI_MALLOC(ddmdq_flg,(matom,3,matom,3,3)) 2192 ddmdq_flg=0 2193 do iq1grad=1,nq1grad 2194 qdir=q1grad(2,iq1grad) 2195 do iatpert=1,natpert 2196 2197 rhog1_tmp(:,:)=rhog1_atdis(iatpert,:,:) 2198 call hartredq(2,gmet,gsqcut,mpi_enreg,nfft,ngfft,qdir,rhog1_tmp,vqgradhart) 2199 2200 !If GGA convert the gradient of xc kernel to reduced coordinates and incorporate it to the Hartree part 2201 if (nkxc == 7) then 2202 vxc1dq=zero 2203 do qcar=1,3 2204 vxc1dq(:,:)=vxc1dq(:,:) + gprimd(qcar,qdir) * & 2205 & vxc1dqc(:,:,iatpert,qcar) 2206 end do 2207 vqgradhart(:)=vqgradhart(:)+vxc1dq(:,1) 2208 end if 2209 2210 !TODO:Maybe it is only necessary to compute half of these elements by symmetry 2211 do jatpert=1,natpert 2212 2213 !Calculate the electrostatic energy term with the first order electric field density 2214 !I need a cplex density for the dotprod_vn 2215 do ii=1,nfft 2216 jj=ii*2 2217 rhor1_cplx(jj-1,:)=rhor1_atdis(jatpert,ii,:) 2218 end do 2219 2220 call dotprod_vn(2,rhor1_cplx,dotr,doti,nfft,nfftot,nspden,2,vqgradhart,ucvol) 2221 ddmdq_qgradhart(re,iatpert,jatpert,iq1grad)=dotr*half 2222 ddmdq_qgradhart(im,iatpert,jatpert,iq1grad)=doti*half 2223 ddmdq_flg(pert_atdis(1,iatpert),pert_atdis(2,iatpert),& 2224 & pert_atdis(1,jatpert),pert_atdis(2,jatpert),q1grad(2,iq1grad))=1 2225 2226 blkflg(pert_atdis(2,iatpert),pert_atdis(1,iatpert), & 2227 & pert_atdis(2,jatpert),pert_atdis(1,jatpert), & 2228 & q1grad(2,iq1grad),matom+8)=1 2229 2230 end do 2231 end do 2232 end do 2233 end if 2234 2235 !1st g-gradient of internal strain tensor contribution 2236 if (lw_flexo==1.or.lw_flexo==4) then 2237 ABI_MALLOC(isdq_qgradhart,(2,matom,3,3,3,3)) 2238 ABI_MALLOC(isdq_flg,(matom,3,3,3,3)) 2239 isdq_flg=0 2240 do iq1grad=1,nq1grad 2241 qdir=q1grad(2,iq1grad) 2242 do iatpert=1,natpert 2243 2244 rhog1_tmp(:,:)=rhog1_atdis(iatpert,:,:) 2245 call hartredq(2,gmet,gsqcut,mpi_enreg,nfft,ngfft,qdir,rhog1_tmp,vqgradhart) 2246 2247 !If GGA convert the gradient of xc kernel to reduced coordinates and incorporate it to the Hartree part 2248 if (nkxc == 7) then 2249 vxc1dq=zero 2250 do qcar=1,3 2251 vxc1dq(:,:)=vxc1dq(:,:) + gprimd(qcar,qdir) * & 2252 & vxc1dqc(:,:,iatpert,qcar) 2253 end do 2254 vqgradhart(:)=vqgradhart(:)+vxc1dq(:,1) 2255 end if 2256 2257 do istrpert=1,nstrpert 2258 2259 !Calculate the electrostatic energy term with the first order strain density 2260 !I need a cplex density for the dotprod_vn 2261 do ii=1,nfft 2262 jj=ii*2 2263 rhor1_cplx(jj-1,:)=rhor1_strain(istrpert,ii,:) 2264 end do 2265 2266 call dotprod_vn(2,rhor1_cplx,dotr,doti,nfft,nfftot,nspden,2,vqgradhart,ucvol) 2267 2268 isdq_qgradhart(re,pert_atdis(1,iatpert),pert_atdis(2,iatpert),q1grad(2,iq1grad), & 2269 & pert_strain(3,istrpert),pert_strain(4,istrpert))=dotr*half 2270 isdq_qgradhart(im,pert_atdis(1,iatpert),pert_atdis(2,iatpert),q1grad(2,iq1grad), & 2271 & pert_strain(3,istrpert),pert_strain(4,istrpert))=doti*half 2272 isdq_flg(pert_atdis(1,iatpert),pert_atdis(2,iatpert),q1grad(2,iq1grad), & 2273 & pert_strain(3,istrpert),pert_strain(4,istrpert))=1 2274 2275 blkflg(pert_atdis(2,iatpert),pert_atdis(1,iatpert), & 2276 & pert_strain(2,istrpert),pert_strain(1,istrpert), & 2277 & q1grad(2,iq1grad),matom+8)=1 2278 2279 end do 2280 end do 2281 end do 2282 end if 2283 2284 ABI_FREE(rhor1_cplx) 2285 ABI_FREE(rhog1_tmp) 2286 ABI_FREE(vqgradhart) 2287 if (nkxc == 7) then 2288 ABI_FREE(vxc1dq) 2289 ABI_SFREE(vxc1dqc) 2290 ABI_SFREE(rhor1_tmp) 2291 end if 2292 if (lw_flexo==1.or.lw_flexo==3.or.lw_flexo==4) then 2293 ABI_FREE(rhog1_atdis) 2294 ABI_FREE(rhor1_atdis) 2295 end if 2296 if (lw_flexo==1.or.lw_flexo==2) then 2297 ABI_FREE(rhog1_efield) 2298 ABI_FREE(rhor1_efield) 2299 end if 2300 if (lw_flexo==1.or.lw_flexo==2.or.lw_flexo==4) then 2301 ABI_FREE(rhor1_strain) 2302 end if 2303 2304 !################# WAVE FUNCTION CONTRIBUTIONS ####################################### 2305 2306 !Determine the subset of symmetry operations (nsym1 operations) 2307 !that leaves the perturbation invariant, and initialize corresponding arrays 2308 !symaf1, symrl1, tnons1 (and pawang1%zarot, if PAW).. 2309 !MR TODO: For the moment only the identiy symmetry is activated (nsym1=1) 2310 ! In a future I will try to activate perturbation dependent symmetries 2311 ! with littlegroup_pert.F90. 2312 nsym1 = 1 2313 ABI_MALLOC(indsy1,(4,nsym1,dtset%natom)) 2314 ABI_MALLOC(symrc1,(3,3,nsym1)) 2315 ABI_MALLOC(symaf1,(nsym1)) 2316 ABI_MALLOC(symrl1,(3,3,nsym1)) 2317 ABI_MALLOC(tnons1,(3,nsym1)) 2318 symaf1(1:nsym1)= 1 2319 symrl1(:,:,nsym1)= dtset%symrel(:,:,1) 2320 tnons1(:,nsym1)= 0_dp 2321 2322 !Set up corresponding symmetry data 2323 ABI_MALLOC(irrzon1,(dtset%nfft**(1-1/nsym1),2,(nspden/dtset%nsppol)-3*(nspden/4))) 2324 ABI_MALLOC(phnons1,(2,dtset%nfft**(1-1/nsym1),(nspden/dtset%nsppol)-3*(nspden/4))) 2325 call setsym(indsy1,irrzon1,1,dtset%natom,dtset%nfft,dtset%ngfft,nspden,dtset%nsppol,& 2326 &nsym1,phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred) 2327 2328 ABI_FREE(indsy1) 2329 ABI_FREE(symaf1) 2330 ABI_FREE(symrl1) 2331 ABI_FREE(tnons1) 2332 2333 !Determine the subset of k-points needed in the "reduced Brillouin zone", 2334 !and initialize other quantities 2335 ABI_MALLOC(indkpt1_tmp,(nkpt)) 2336 ABI_MALLOC(wtk_folded,(nkpt)) 2337 ABI_MALLOC(bz2ibz_smap,(6, nkpt)) 2338 indkpt1_tmp(:)=0 2339 2340 if (dtset%kptopt==2) then 2341 call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,& 2342 & nsym1,symrc1,timrev,dtset%wtk,wtk_folded,bz2ibz_smap,xmpi_comm_self) 2343 else if (dtset%kptopt==3) then 2344 call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,& 2345 & nsym1,symrc1,0,dtset%wtk,wtk_folded,bz2ibz_smap,xmpi_comm_self) 2346 else 2347 write(msg,"(1a)") 'kptopt must be 2 or 3 for the quadrupole calculation' 2348 ABI_BUG(msg) 2349 end if 2350 ABI_FREE(bz2ibz_smap) 2351 2352 ABI_MALLOC(doccde_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol)) 2353 ABI_MALLOC(indkpt1,(nkpt_rbz)) 2354 ABI_MALLOC(istwfk_rbz,(nkpt_rbz)) 2355 ABI_MALLOC(kpt_rbz,(3,nkpt_rbz)) 2356 ABI_MALLOC(nband_rbz,(nkpt_rbz*dtset%nsppol)) 2357 ABI_MALLOC(occ_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol)) 2358 ABI_MALLOC(wtk_rbz,(nkpt_rbz)) 2359 indkpt1(:)=indkpt1_tmp(1:nkpt_rbz) 2360 do ikpt=1,nkpt_rbz 2361 istwfk_rbz(ikpt)=dtset%istwfk(indkpt1(ikpt)) 2362 kpt_rbz(:,ikpt)=dtset%kptns(:,indkpt1(ikpt)) 2363 wtk_rbz(ikpt)=wtk_folded(indkpt1(ikpt)) 2364 end do 2365 ABI_FREE(indkpt1_tmp) 2366 ABI_FREE(wtk_folded) 2367 2368 !Transfer occ to occ_rbz 2369 !NOTE : this takes into account that indkpt1 is ordered 2370 !MG: What about using occ(band,kpt,spin) ??? 2371 bdtot_index=0;bdtot1_index=0 2372 do isppol=1,dtset%nsppol 2373 ikpt1=1 2374 do ikpt=1,nkpt 2375 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 2376 ! Must test against ikpt1/=nkpt_rbz+1, before evaluate indkpt1(ikpt1) 2377 if(ikpt1/=nkpt_rbz+1)then 2378 if(ikpt==indkpt1(ikpt1))then 2379 nband_rbz(ikpt1+(isppol-1)*nkpt_rbz)=nband_k 2380 occ_rbz(1+bdtot1_index:nband_k+bdtot1_index)=occ(1+bdtot_index:nband_k+bdtot_index) 2381 doccde_rbz(1+bdtot1_index:nband_k+bdtot1_index)=doccde(1+bdtot_index:nband_k+bdtot_index) 2382 ikpt1=ikpt1+1 2383 bdtot1_index=bdtot1_index+nband_k 2384 end if 2385 end if 2386 bdtot_index=bdtot_index+nband_k 2387 end do 2388 end do 2389 2390 !Compute maximum number of planewaves at k 2391 call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpt_rbz,mpi_enreg,mpw,nkpt_rbz) 2392 2393 !Allocate some k-dependent arrays at k 2394 ABI_MALLOC(kg,(3,mpw*nkpt_rbz)) 2395 ABI_MALLOC(kg_k,(3,mpw)) 2396 ABI_MALLOC(npwarr,(nkpt_rbz)) 2397 ABI_MALLOC(npwtot,(nkpt_rbz)) 2398 2399 !Determine distribution of k-points/bands over MPI processes 2400 if (allocated(mpi_enreg%my_kpttab)) then 2401 ABI_FREE(mpi_enreg%my_kpttab) 2402 end if 2403 ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt_rbz)) 2404 if(xmpi_paral==1) then 2405 ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt_rbz,dtset%mband,dtset%nsppol)) 2406 call distrb2(dtset%mband,mband_mem,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,dtset%nsppol,mpi_enreg) 2407 else 2408 mpi_enreg%my_kpttab(:)=(/(ii,ii=1,nkpt_rbz)/) 2409 end if 2410 my_nkpt_rbz=maxval(mpi_enreg%my_kpttab) 2411 mkmem_rbz =my_nkpt_rbz 2412 call initmpi_band(mkmem_rbz,mpi_enreg,nband_rbz,nkpt_rbz,dtset%nsppol) 2413 2414 !Set up the basis sphere of planewaves at k 2415 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg,& 2416 & kpt_rbz,mkmem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw,npwarr,npwtot,dtset%nsppol) 2417 ABI_FREE(npwtot) 2418 2419 !Set up the spherical harmonics (Ylm) and 1st gradients at k 2420 useylmgr=1; option=2 ; nylmgr=9 2421 ABI_MALLOC(ylm,(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm)) 2422 ABI_MALLOC(ylmgr,(mpw*mkmem_rbz,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 2423 if (psps%useylm==1) then 2424 call initylmg(gprimd,kg,kpt_rbz,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,nband_rbz,nkpt_rbz,& 2425 & npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr) 2426 end if 2427 2428 !Initialize band structure datatype at k 2429 bantot_rbz=sum(nband_rbz(1:nkpt_rbz*dtset%nsppol)) 2430 ABI_MALLOC(eigen0,(bantot_rbz)) 2431 eigen0(:)=zero 2432 ! CP modified 2433 ! call ebands_init(bantot_rbz,bs_rbz,dtset%nelect,doccde_rbz,eigen0,istwfk_rbz,kpt_rbz,& 2434 !& nband_rbz,nkpt_rbz,npwarr,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,& 2435 !& dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 2436 !& dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 2437 call ebands_init(bantot_rbz,bs_rbz,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,& 2438 & doccde_rbz,eigen0,istwfk_rbz,kpt_rbz,nband_rbz,nkpt_rbz,npwarr,dtset%nsppol,dtset%nspinor,& 2439 & dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz, dtset%cellcharge(1), dtset%kptopt, & 2440 & dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 2441 ! End CP modified 2442 ABI_FREE(eigen0) 2443 ABI_FREE(doccde_rbz) 2444 2445 !Initialize header, update it with evolving variables 2446 gscase=0 ! A GS WF file is read 2447 call hdr_init(bs_rbz,codvsn,dtset,hdr0,pawtab,gscase,psps,wvl%descr,& 2448 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 2449 2450 ! CP modified 2451 ! call hdr0%update(bantot_rbz,etotal,fermie,& 2452 !& residm,rprimd,occ_rbz,pawrhoij,xred,dtset%amu_orig(:,1),& 2453 !& comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 2454 call hdr0%update(bantot_rbz,etotal,fermie,fermie,& ! CP modified call to function; we assume occopt 9 does not work with DFPT 2455 & residm,rprimd,occ_rbz,pawrhoij,xred,dtset%amu_orig(:,1),& 2456 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 2457 ! End CP modified 2458 2459 !Clean band structure datatype (should use it more in the future !) 2460 call ebands_free(bs_rbz) 2461 2462 !Initialize GS wavefunctions at k 2463 ireadwf0=1; formeig=0 ; ask_accurate=1 ; optorth=0 2464 mcg=mpw*dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol 2465 if (one*mpw*dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol > huge(1)) then 2466 write (msg,'(4a, 5(a,i0), 2a)')& 2467 & "Default integer is not wide enough to store the size of the GS wavefunction array (WF0, mcg).",ch10,& 2468 & "Action: increase the number of processors. Consider also OpenMP threads.",ch10,& 2469 & "nspinor: ",dtset%nspinor, "mpw: ",mpw, "mband: ",dtset%mband, "mkmem_rbz: ",& 2470 & mkmem_rbz, "nsppol: ",dtset%nsppol,ch10,& 2471 & 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode' 2472 ABI_ERROR(msg) 2473 end if 2474 ABI_STAT_MALLOC(cg,(2,mcg), ierr) 2475 ABI_CHECK(ierr==0, "out-of-memory in cg") 2476 2477 ABI_MALLOC(eigen0,(dtset%mband*nkpt_rbz*dtset%nsppol)) 2478 call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,& 2479 & formeig,hdr0,ireadwf0,istwfk_rbz,kg,& 2480 & kpt_rbz,dtset%localrdwf,dtset%mband,mcg,& 2481 & mkmem_rbz,mpi_enreg,mpw,nband_rbz,dtset%ngfft,nkpt_rbz,npwarr,& 2482 & dtset%nsppol,dtset%nsym,occ_rbz,optorth,dtset%symafm,& 2483 & dtset%symrel,dtset%tnons,dtfil%unkg,wffgs,wfftgs,& 2484 & dtfil%unwffgs,dtfil%fnamewffk,wvl) 2485 ABI_FREE(eigen0) 2486 !Close wffgs%unwff, if it was ever opened (in inwffil) 2487 if (ireadwf0==1) then 2488 call WffClose(wffgs,ierr) 2489 end if 2490 2491 !==== Initialize most of the Hamiltonian ==== 2492 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 2493 !2) Perform the setup needed for the non-local factors: 2494 !3) Constant kleimann-Bylander energies are copied from psps to gs_hamkq. 2495 call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,nsppol,nspden,dtset%natom,& 2496 & dtset%typat,xred,nfft,dtset%mgfft,ngfft,rprimd,dtset%nloalg,ph1d=ph1d,& 2497 & use_gpu_cuda=dtset%use_gpu_cuda) 2498 2499 2500 !==== Initialize response functions files and handlers ==== 2501 !Atomic displacement 2502 if (lw_flexo==1.or.lw_flexo==3.or.lw_flexo==4) then 2503 ABI_MALLOC(wfk_t_atdis,(natpert)) 2504 do iatpert=1,natpert 2505 pertcase=pert_atdis(3,iatpert) 2506 call appdig(pertcase,dtfil%fnamewff1,fiwfatdis) 2507 2508 !The value 20 is taken arbitrarily I would say 2509 forunit=20+pertcase 2510 2511 !Check that atdis file exists and open it 2512 if (.not. file_exists(fiwfatdis)) then 2513 ! Trick needed to run Abinit test suite in netcdf mode. 2514 if (file_exists(nctk_ncify(fiwfatdis))) then 2515 write(msg,"(3a)")"- File: ",trim(fiwfatdis),& 2516 " does not exist but found netcdf file with similar name." 2517 call wrtout(std_out,msg,'COLL') 2518 fiwfatdis = nctk_ncify(fiwfatdis) 2519 end if 2520 if (.not. file_exists(fiwfatdis)) then 2521 ABI_ERROR('Missing file: '//TRIM(fiwfatdis)) 2522 end if 2523 end if 2524 write(msg,'(a,a)')'-open atomic displacement wf1 file :',trim(fiwfatdis) 2525 call wrtout(std_out,msg,'COLL') 2526 call wrtout(ab_out,msg,'COLL') 2527 call wfk_open_read(wfk_t_atdis(iatpert),fiwfatdis,formeig1,dtset%iomode,forunit,spaceworld) 2528 end do 2529 end if 2530 2531 !ddk files 2532 ABI_MALLOC(wfk_t_ddk,(nq1grad)) 2533 do iq1grad=1,nq1grad 2534 pertcase=q1grad(3,iq1grad) 2535 call appdig(pertcase,dtfil%fnamewffddk,fiwfddk) 2536 2537 !The value 20 is taken arbitrarily I would say 2538 forunit=20+pertcase 2539 2540 !Check that ddk file exists and open it 2541 if (.not. file_exists(fiwfddk)) then 2542 ! Trick needed to run Abinit test suite in netcdf mode. 2543 if (file_exists(nctk_ncify(fiwfddk))) then 2544 write(msg,"(3a)")"- File: ",trim(fiwfddk),& 2545 " does not exist but found netcdf file with similar name." 2546 call wrtout(std_out,msg,'COLL') 2547 fiwfddk = nctk_ncify(fiwfddk) 2548 end if 2549 if (.not. file_exists(fiwfddk)) then 2550 ABI_ERROR('Missing file: '//TRIM(fiwfddk)) 2551 end if 2552 end if 2553 write(msg, '(a,a)') '-open ddk wf1 file :',trim(fiwfddk) 2554 call wrtout(std_out,msg,'COLL') 2555 call wrtout(ab_out,msg,'COLL') 2556 call wfk_open_read(wfk_t_ddk(iq1grad),fiwfddk,formeig1,dtset%iomode,forunit,spaceworld) 2557 end do 2558 2559 !Electric field files 2560 if (lw_flexo==1.or.lw_flexo==2) then 2561 ABI_MALLOC(wfk_t_efield,(nefipert)) 2562 do iefipert=1,nefipert 2563 pertcase=pert_efield(3,iefipert) 2564 call appdig(pertcase,dtfil%fnamewff1,fiwfefield) 2565 2566 !The value 20 is taken arbitrarily I would say 2567 forunit=20+pertcase 2568 2569 !Check that efield file exists and open it 2570 if (.not. file_exists(fiwfefield)) then 2571 ! Trick needed to run Abinit test suite in netcdf mode. 2572 if (file_exists(nctk_ncify(fiwfefield))) then 2573 write(msg,"(3a)")"- File: ",trim(fiwfefield),& 2574 " does not exist but found netcdf file with similar name." 2575 call wrtout(std_out,msg,'COLL') 2576 fiwfefield = nctk_ncify(fiwfefield) 2577 end if 2578 if (.not. file_exists(fiwfefield)) then 2579 ABI_ERROR('Missing file: '//TRIM(fiwfefield)) 2580 end if 2581 end if 2582 write(msg, '(a,a)') '-open electric field wf1 file :',trim(fiwfefield) 2583 call wrtout(std_out,msg,'COLL') 2584 call wrtout(ab_out,msg,'COLL') 2585 call wfk_open_read(wfk_t_efield(iefipert),fiwfefield,formeig1,dtset%iomode,forunit,spaceworld) 2586 end do 2587 end if 2588 2589 !Strain files 2590 if (lw_flexo==1.or.lw_flexo==2.or.lw_flexo==4) then 2591 ABI_MALLOC(wfk_t_strain,(3,3)) 2592 do istrpert=1,nstrpert 2593 pertcase=pert_strain(5,istrpert) 2594 ka=pert_strain(3,istrpert) 2595 kb=pert_strain(4,istrpert) 2596 call appdig(pertcase,dtfil%fnamewff1,fiwfstrain) 2597 2598 !The value 20 is taken arbitrarily I would say 2599 forunit=20+pertcase 2600 2601 !Check that strain file exists and open it 2602 if (.not. file_exists(fiwfstrain)) then 2603 ! Trick needed to run Abinit test suite in netcdf mode. 2604 if (file_exists(nctk_ncify(fiwfstrain))) then 2605 write(msg,"(3a)")"- File: ",trim(fiwfstrain),& 2606 " does not exist but found netcdf file with similar name." 2607 call wrtout(std_out,msg,'COLL') 2608 fiwfstrain = nctk_ncify(fiwfstrain) 2609 end if 2610 if (.not. file_exists(fiwfstrain)) then 2611 ABI_ERROR('Missing file: '//TRIM(fiwfstrain)) 2612 end if 2613 end if 2614 if (ka>=kb) then 2615 write(msg, '(a,a)') '-open strain wf1 file :',trim(fiwfstrain) 2616 call wrtout(std_out,msg,'COLL') 2617 call wrtout(ab_out,msg,'COLL') 2618 call wfk_open_read(wfk_t_strain(ka,kb),fiwfstrain,formeig1,dtset%iomode,forunit,spaceworld) 2619 else 2620 wfk_t_strain(ka,kb)=wfk_t_strain(kb,ka) 2621 end if 2622 2623 end do 2624 end if 2625 2626 !d2_dkdk 2627 if (lw_flexo==1.or.lw_flexo==2) then 2628 ABI_MALLOC(wfk_t_dkdk,(nq1q2grad)) 2629 do iq1q2grad=1,nq1q2grad 2630 2631 pertcase=q1q2grad(4,iq1q2grad) 2632 call appdig(pertcase,dtfil%fnamewffdkdk,fiwfdkdk) 2633 2634 !The value 20 is taken arbitrarily I would say 2635 forunit=20+pertcase 2636 2637 !Check that d2_ddk file exists and open it 2638 if (.not. file_exists(fiwfdkdk)) then 2639 ! Trick needed to run Abinit test suite in netcdf mode. 2640 if (file_exists(nctk_ncify(fiwfdkdk))) then 2641 write(msg,"(3a)")"- File: ",trim(fiwfdkdk),& 2642 " does not exist but found netcdf file with similar name." 2643 call wrtout(std_out,msg,'COLL') 2644 fiwfdkdk = nctk_ncify(fiwfdkdk) 2645 end if 2646 if (.not. file_exists(fiwfdkdk)) then 2647 ABI_ERROR('Missing file: '//TRIM(fiwfdkdk)) 2648 end if 2649 end if 2650 if (iq1q2grad <= 6) then 2651 write(msg, '(a,a)') '-open d2_dkdk wf2 file :',trim(fiwfdkdk) 2652 call wrtout(std_out,msg,'COLL') 2653 call wrtout(ab_out,msg,'COLL') 2654 call wfk_open_read(wfk_t_dkdk(iq1q2grad),fiwfdkdk,formeig1,dtset%iomode,forunit,spaceworld) 2655 else 2656 wfk_t_dkdk(iq1q2grad)=wfk_t_dkdk(iq1q2grad-3) 2657 end if 2658 end do 2659 end if 2660 2661 !Allocate the electronic flexoelectric tensor part depending on the wave functions 2662 if (lw_flexo==1.or.lw_flexo==2) then 2663 ABI_MALLOC(elflexowf,(2,3,3,3,3)) 2664 ABI_MALLOC(elflexowf_k,(2,3,3,3,3)) 2665 ABI_MALLOC(elflexowf_t1,(2,3,3,3,3)) 2666 ABI_MALLOC(elflexowf_t1_k,(2,3,3,3,3)) 2667 ABI_MALLOC(elflexowf_t2,(2,3,3,3,3)) 2668 ABI_MALLOC(elflexowf_t2_k,(2,3,3,3,3)) 2669 ABI_MALLOC(elflexowf_t3,(2,3,3,3,3)) 2670 ABI_MALLOC(elflexowf_t3_k,(2,3,3,3,3)) 2671 ABI_MALLOC(elflexowf_t4,(2,3,3,3,3)) 2672 ABI_MALLOC(elflexowf_t4_k,(2,3,3,3,3)) 2673 ABI_MALLOC(elflexowf_t5,(2,3,3,3,3)) 2674 ABI_MALLOC(elflexowf_t5_k,(2,3,3,3,3)) 2675 elflexowf=zero 2676 elflexowf_t1=zero 2677 elflexowf_t2=zero 2678 elflexowf_t3=zero 2679 elflexowf_t4=zero 2680 elflexowf_t5=zero 2681 end if 2682 2683 !Allocate arrays for wf contributions to the first q-gradient of the dynamical matrix 2684 if (lw_flexo==1.or.lw_flexo==3) then 2685 ABI_MALLOC(ddmdqwf,(2,natpert,natpert,nq1grad)) 2686 ABI_MALLOC(ddmdqwf_k,(2,natpert,natpert,nq1grad)) 2687 ABI_MALLOC(ddmdqwf_t1,(2,natpert,natpert,nq1grad)) 2688 ABI_MALLOC(ddmdqwf_t1_k,(2,natpert,natpert,nq1grad)) 2689 ABI_MALLOC(ddmdqwf_t2,(2,natpert,natpert,nq1grad)) 2690 ABI_MALLOC(ddmdqwf_t2_k,(2,natpert,natpert,nq1grad)) 2691 ABI_MALLOC(ddmdqwf_t3,(2,natpert,natpert,nq1grad)) 2692 ABI_MALLOC(ddmdqwf_t3_k,(2,natpert,natpert,nq1grad)) 2693 ddmdqwf=zero 2694 ddmdqwf_t1=zero 2695 ddmdqwf_t2=zero 2696 ddmdqwf_t3=zero 2697 end if 2698 2699 !Allocate arrays for wf contributions to the first q-gradient of the internal strain tensor 2700 if (lw_flexo==1.or.lw_flexo==4) then 2701 ABI_MALLOC(frwfdq,(2,matom,3,3,3,nq1grad)) 2702 ABI_MALLOC(frwfdq_k,(2,matom,3,3,3,nq1grad)) 2703 ABI_MALLOC(isdqwf,(2,matom,3,nq1grad,3,3)) 2704 ABI_MALLOC(isdqwf_k,(2,matom,3,nq1grad,3,3)) 2705 ABI_MALLOC(isdqwf_t1,(2,matom,3,nq1grad,3,3)) 2706 ABI_MALLOC(isdqwf_t1_k,(2,matom,3,nq1grad,3,3)) 2707 ABI_MALLOC(isdqwf_t2,(2,matom,3,nq1grad,3,3)) 2708 ABI_MALLOC(isdqwf_t2_k,(2,matom,3,nq1grad,3,3)) 2709 ABI_MALLOC(isdqwf_t3,(2,matom,3,nq1grad,3,3)) 2710 ABI_MALLOC(isdqwf_t3_k,(2,matom,3,nq1grad,3,3)) 2711 ABI_MALLOC(isdqwf_t4,(2,matom,3,3,3,nq1grad)) 2712 ABI_MALLOC(isdqwf_t4_k,(2,matom,3,3,3,nq1grad)) 2713 ABI_MALLOC(isdqwf_t5,(2,matom,3,nq1grad,3,3)) 2714 ABI_MALLOC(isdqwf_t5_k,(2,matom,3,nq1grad,3,3)) 2715 frwfdq=zero 2716 isdqwf=zero 2717 isdqwf_t1=zero 2718 isdqwf_t2=zero 2719 isdqwf_t3=zero 2720 isdqwf_t4=zero 2721 isdqwf_t5=zero 2722 end if 2723 2724 !LOOP OVER SPINS 2725 bdtot_index=0 2726 icg=0 2727 do isppol=1,nsppol 2728 ikg=0 2729 2730 ! Continue to initialize the Hamiltonian 2731 call gs_hamkq%load_spin(isppol,with_nonlocal=.true.) 2732 2733 ! BIG FAT k POINT LOOP 2734 do ikpt=1,nkpt_rbz 2735 2736 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 2737 istwf_k=istwfk_rbz(ikpt) 2738 npw_k=npwarr(ikpt) 2739 npw1_k=npw_k 2740 2741 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 2742 bdtot_index=bdtot_index+nband_k 2743 2744 cycle ! Skip the rest of the k-point loop 2745 end if 2746 2747 ABI_MALLOC(occ_k,(nband_k)) 2748 ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 2749 ABI_MALLOC(ylmgr_k,(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 2750 occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index) 2751 kpoint(:)=kpt_rbz(:,ikpt) 2752 wtk_k=wtk_rbz(ikpt) 2753 2754 ! Get plane-wave vectors and related data at k 2755 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 2756 if (psps%useylm==1) then 2757 do ilm=1,psps%mpsang*psps%mpsang 2758 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 2759 end do 2760 if (useylmgr==1) then 2761 do ilm=1,psps%mpsang*psps%mpsang 2762 do ii=1,nylmgr 2763 ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm) 2764 end do 2765 end do 2766 end if 2767 end if 2768 2769 ! Compute the wf contributions to the electronic flexoelectric tensor 2770 if (lw_flexo==1.or.lw_flexo==2) then 2771 call dfpt_ciflexowf(cg,cplex,dtset,elflexowf_k,elflexowf_t1_k,elflexowf_t2_k, & 2772 & elflexowf_t3_k,elflexowf_t4_k,elflexowf_t5_k, & 2773 & gs_hamkq,gsqcut,icg,ikpt,indkpt1,isppol,istwf_k, & 2774 & kg_k,kpoint,kxc,mkmem_rbz, & 2775 & mpi_enreg,mpw,nattyp,nband_k,nefipert,nfft,ngfft,nkpt_rbz,nkxc, & 2776 & npw_k,nq1grad,nq1q2grad,nspden,nsppol,nstrpert,nylmgr,occ_k, & 2777 & pert_efield,pert_strain,ph1d,psps,q1grad,q1q2grad,rhog,rhor,rmet,ucvol,useylmgr, & 2778 & vhxc1_efield,vhxc1_strain,wfk_t_efield,wfk_t_ddk, & 2779 & wfk_t_dkdk,wfk_t_strain,wtk_k,ylm_k,ylmgr_k) 2780 2781 ! Add the contribution from each k-point 2782 elflexowf=elflexowf + elflexowf_k 2783 elflexowf_t1=elflexowf_t1 + elflexowf_t1_k 2784 elflexowf_t2=elflexowf_t2 + elflexowf_t2_k 2785 elflexowf_t3=elflexowf_t3 + elflexowf_t3_k 2786 elflexowf_t4=elflexowf_t4 + elflexowf_t4_k 2787 elflexowf_t5=elflexowf_t5 + elflexowf_t5_k 2788 end if 2789 2790 ! Compute the wf contributions to the first q-gradient of the dynamical matrix 2791 if (lw_flexo==1.or.lw_flexo==3) then 2792 call dfpt_ddmdqwf(atindx,cg,cplex,ddmdqwf_k,ddmdqwf_t1_k,ddmdqwf_t2_k,ddmdqwf_t3_k,& 2793 & dtset,gs_hamkq,gsqcut,icg,ikpt,indkpt1,isppol,istwf_k,kg_k,kpoint,mkmem_rbz, & 2794 & mpi_enreg,mpw,natpert,nattyp,nband_k,nfft,ngfft,nkpt_rbz,npw_k,nq1grad,nspden, & 2795 & nsppol,nylmgr,occ_k,pert_atdis,ph1d,psps,q1grad,rmet,ucvol,useylmgr, & 2796 & vhxc1_atdis,wfk_t_atdis,wfk_t_ddk,wtk_k,xred,ylm_k,ylmgr_k) 2797 2798 ! Add the contribution from each k-point 2799 ddmdqwf=ddmdqwf + ddmdqwf_k 2800 ddmdqwf_t1=ddmdqwf_t1 + ddmdqwf_t1_k 2801 ddmdqwf_t2=ddmdqwf_t2 + ddmdqwf_t2_k 2802 ddmdqwf_t3=ddmdqwf_t3 + ddmdqwf_t3_k 2803 end if 2804 2805 ! Compute the wf contributions to the first q-gradient of the internal strain tensor 2806 if (lw_flexo==1.or.lw_flexo==4) then 2807 2808 ! First calculate the frozen wf contribution (notice the type-I indexing of 2809 ! this term) 2810 call dfpt_isdqfr(atindx,cg,dtset,frwfdq_k,gs_hamkq,gsqcut,icg,ikpt,& 2811 & isppol,istwf_k,kg_k,kpoint,mkmem_rbz,mpi_enreg,matom,mpw,natpert,nattyp,nband_k,nfft,& 2812 & ngfft,npw_k,nq1grad,nspden,nsppol,nstrpert,nylmgr,occ_k,pert_atdis, & 2813 & pert_strain,ph1d,psps,rmet,ucvol,useylmgr,wtk_k,ylm_k,ylmgr_k) 2814 2815 ! Add the contribution from each k-point 2816 frwfdq=frwfdq + frwfdq_k 2817 2818 ! Now comute the 1st order wf contributions 2819 call dfpt_isdqwf(atindx,cg,cplex,dtset,gs_hamkq,gsqcut,icg,ikpt,indkpt1,isdqwf_k, & 2820 & isdqwf_t1_k,isdqwf_t2_k,isdqwf_t3_k,isdqwf_t4_k,isdqwf_t5_k,isppol,istwf_k, & 2821 & kg_k,kpoint,kxc,matom,mkmem_rbz,mpi_enreg,mpw,natpert,nattyp,nband_k,nfft,ngfft,nkpt_rbz,nkxc, & 2822 & npw_k,nq1grad,nspden,nsppol,nstrpert,nylmgr,occ_k, & 2823 & pert_atdis,pert_strain,ph1d,psps,q1grad,rhog,rhor,rmet,ucvol,useylmgr, & 2824 & vhxc1_atdis,vhxc1_strain,wfk_t_atdis,wfk_t_ddk, & 2825 & wfk_t_strain,wtk_k,xred,ylm_k,ylmgr_k) 2826 2827 ! Add the contribution from each k-point 2828 isdqwf=isdqwf + isdqwf_k 2829 isdqwf_t1=isdqwf_t1 + isdqwf_t1_k 2830 isdqwf_t2=isdqwf_t2 + isdqwf_t2_k 2831 isdqwf_t3=isdqwf_t3 + isdqwf_t3_k 2832 isdqwf_t4=isdqwf_t4 + isdqwf_t4_k 2833 isdqwf_t5=isdqwf_t5 + isdqwf_t5_k 2834 2835 end if 2836 2837 ! Keep track of total number of bands 2838 bdtot_index=bdtot_index+nband_k 2839 2840 ! Shift arrays memory 2841 if (mkmem_rbz/=0) then 2842 icg=icg+npw_k*dtset%nspinor*nband_k 2843 ikg=ikg+npw_k 2844 end if 2845 2846 ABI_FREE(occ_k) 2847 ABI_FREE(ylm_k) 2848 ABI_FREE(ylmgr_k) 2849 2850 end do 2851 ! END BIG FAT k POINT LOOP 2852 end do 2853 !END LOOP OVER SPINS 2854 2855 !Memory cleaning 2856 call gs_hamkq%free() 2857 2858 !Close response function files 2859 if (lw_flexo==1.or.lw_flexo==3.or.lw_flexo==4) then 2860 do iatpert=1,natpert 2861 call wfk_t_atdis(iatpert)%close() 2862 end do 2863 end if 2864 if (lw_flexo==1.or.lw_flexo==2.or.lw_flexo==4) then 2865 do istrpert=1,nstrpert 2866 if (istrpert <= 6) then 2867 ka=pert_strain(3,istrpert) 2868 kb=pert_strain(4,istrpert) 2869 call wfk_t_strain(ka,kb)%close() 2870 end if 2871 end do 2872 end if 2873 do iq1grad=1,nq1grad 2874 call wfk_t_ddk(iq1grad)%close() 2875 end do 2876 if (lw_flexo==1.or.lw_flexo==2) then 2877 do iefipert=1,nefipert 2878 call wfk_t_efield(iefipert)%close() 2879 end do 2880 do iq1q2grad=1,nq1q2grad 2881 if (iq1q2grad <= 6) call wfk_t_dkdk(iq1q2grad)%close() 2882 end do 2883 end if 2884 2885 !=== MPI communications ================== 2886 if (xmpi_paral==1) then 2887 2888 if (lw_flexo==1.or.lw_flexo==2) then 2889 call xmpi_sum(elflexowf,spaceworld,ierr) 2890 call xmpi_sum(elflexowf_t1,spaceworld,ierr) 2891 call xmpi_sum(elflexowf_t2,spaceworld,ierr) 2892 call xmpi_sum(elflexowf_t3,spaceworld,ierr) 2893 call xmpi_sum(elflexowf_t4,spaceworld,ierr) 2894 call xmpi_sum(elflexowf_t5,spaceworld,ierr) 2895 end if 2896 2897 if (lw_flexo==1.or.lw_flexo==3) then 2898 call xmpi_sum(ddmdqwf,spaceworld,ierr) 2899 call xmpi_sum(ddmdqwf_t1,spaceworld,ierr) 2900 call xmpi_sum(ddmdqwf_t2,spaceworld,ierr) 2901 call xmpi_sum(ddmdqwf_t3,spaceworld,ierr) 2902 end if 2903 2904 if (lw_flexo==1.or.lw_flexo==4) then 2905 call xmpi_sum(frwfdq,spaceworld,ierr) 2906 call xmpi_sum(isdqwf,spaceworld,ierr) 2907 call xmpi_sum(isdqwf_t1,spaceworld,ierr) 2908 call xmpi_sum(isdqwf_t2,spaceworld,ierr) 2909 call xmpi_sum(isdqwf_t3,spaceworld,ierr) 2910 call xmpi_sum(isdqwf_t4,spaceworld,ierr) 2911 call xmpi_sum(isdqwf_t5,spaceworld,ierr) 2912 end if 2913 2914 end if 2915 2916 !Sum the G=0 contribution to the geometric term of the first 2917 !q-gradient of the internal !strain tensor 2918 if (lw_flexo==1.or.lw_flexo==4) then 2919 fac=pi*pi 2920 2921 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATIONS 2922 do iatpert=1,natpert 2923 iatom=pert_atdis(1,iatpert) 2924 iatdir=pert_atdis(2,iatpert) 2925 2926 !Determination of the atom type 2927 ia1=0 2928 itypat=0 2929 do ii=1,dtset%ntypat 2930 ia1=ia1+nattyp(ii) 2931 if(atindx(iatom)<=ia1.and.itypat==0)itypat=ii 2932 end do 2933 2934 !LOOP OVER STRAIN PERTURBATIONS 2935 do istrpert= 1, nstrpert 2936 ka=pert_strain(3,istrpert) 2937 kb=pert_strain(4,istrpert) 2938 delad=zero ; if (iatdir==kb) delad=one 2939 delbd=zero ; if (ka==kb) delbd=one 2940 2941 !LOOP OVER Q1-GRADIENT 2942 do iq1grad=1,3 2943 delag=zero ; if(iatdir==iq1grad) delag=one 2944 delbg=zero ; if(ka==iq1grad) delbg=one 2945 2946 frwfdq(1,iatom,iatdir,ka,kb,iq1grad)=frwfdq(1,iatom,iatdir,ka,kb,iq1grad)+& 2947 & fac*rhog(1,1)*psps%vlspl(1,2,itypat)*(delag*delbd+delad*delbg) 2948 2949 end do 2950 end do 2951 end do 2952 end if 2953 2954 !Anounce finalization of calculations 2955 if (lw_flexo==1.or.lw_flexo==2) then 2956 write(msg, '(a,a,a)' ) ch10, & 2957 ' Clamped-ion flexoelectric tensor calculation completed ',ch10 2958 call wrtout(std_out,msg,'COLL') 2959 call wrtout(ab_out,msg,'COLL') 2960 end if 2961 if (lw_flexo==1.or.lw_flexo==3) then 2962 write(msg, '(a,a,a)' ) ch10, & 2963 ' Dynamical matrix 1st q-gradient calculation completed ',ch10 2964 call wrtout(std_out,msg,'COLL') 2965 call wrtout(ab_out,msg,'COLL') 2966 end if 2967 if (lw_flexo==1.or.lw_flexo==4) then 2968 write(msg, '(a,a,a)' ) ch10, & 2969 ' Piezoelectric force response 1st q-gradient calculation completed ',ch10 2970 call wrtout(std_out,msg,'COLL') 2971 call wrtout(ab_out,msg,'COLL') 2972 end if 2973 2974 !Gather the different terms in the flexoelectric tensor and print them out 2975 if (me==0) then 2976 if (lw_flexo==1.or.lw_flexo==2) then 2977 call dfpt_ciflexoout(d3etot,elflexoflg,elflexowf,elflexowf_t1,elflexowf_t2, & 2978 & elflexowf_t3,elflexowf_t4,elflexowf_t5, & 2979 & elqgradhart,gprimd,dtset%kptopt,matom,mpert,nefipert, & 2980 & nstrpert,nq1grad,pert_efield,pert_strain,dtset%prtvol,q1grad,rprimd,ucvol) 2981 end if 2982 if (lw_flexo==1.or.lw_flexo==3) then 2983 call dfpt_ddmdqout(ddmdq_flg,ddmdq_qgradhart,ddmdqwf,ddmdqwf_t1,ddmdqwf_t2,ddmdqwf_t3,d3etot,& 2984 & dyewdq,gprimd,dtset%kptopt,matom,mpert,natpert,nq1grad,pert_atdis,dtset%prtvol,q1grad,rprimd) 2985 end if 2986 if (lw_flexo==1.or.lw_flexo==4) then 2987 call dfpt_isdqout(d3etot,dyewdqdq,frwfdq,gprimd,isdq_flg,isdq_qgradhart,isdqwf,isdqwf_t1,isdqwf_t2,& 2988 & isdqwf_t3,isdqwf_t4,isdqwf_t5,dtset%kptopt,matom,mpert,natpert, & 2989 & nstrpert,nq1grad,pert_atdis,pert_strain,dtset%prtvol,q1grad,rprimd,ucvol) 2990 end if 2991 end if 2992 2993 !Deallocattions 2994 if (lw_flexo==1.or.lw_flexo==3.or.lw_flexo==4) then 2995 ABI_FREE(pert_atdis) 2996 ABI_FREE(vhxc1_atdis) 2997 ABI_FREE(wfk_t_atdis) 2998 end if 2999 if (lw_flexo==1.or.lw_flexo==2) then 3000 ABI_FREE(pert_efield) 3001 ABI_FREE(q1q2grad) 3002 ABI_FREE(vhxc1_efield) 3003 ABI_FREE(elqgradhart) 3004 ABI_FREE(elflexoflg) 3005 ABI_FREE(wfk_t_efield) 3006 ABI_FREE(wfk_t_dkdk) 3007 ABI_FREE(elflexowf) 3008 ABI_FREE(elflexowf_k) 3009 ABI_FREE(elflexowf_t1) 3010 ABI_FREE(elflexowf_t1_k) 3011 ABI_FREE(elflexowf_t2) 3012 ABI_FREE(elflexowf_t2_k) 3013 ABI_FREE(elflexowf_t3) 3014 ABI_FREE(elflexowf_t3_k) 3015 ABI_FREE(elflexowf_t4) 3016 ABI_FREE(elflexowf_t4_k) 3017 ABI_FREE(elflexowf_t5) 3018 ABI_FREE(elflexowf_t5_k) 3019 end if 3020 if (lw_flexo==1.or.lw_flexo==3) then 3021 ABI_FREE(ddmdq_qgradhart) 3022 ABI_FREE(ddmdq_flg) 3023 ABI_FREE(ddmdqwf) 3024 ABI_FREE(ddmdqwf_k) 3025 ABI_FREE(ddmdqwf_t1) 3026 ABI_FREE(ddmdqwf_t1_k) 3027 ABI_FREE(ddmdqwf_t2) 3028 ABI_FREE(ddmdqwf_t2_k) 3029 ABI_FREE(ddmdqwf_t3) 3030 ABI_FREE(ddmdqwf_t3_k) 3031 end if 3032 if (lw_flexo==1.or.lw_flexo==2.or.lw_flexo==4) then 3033 ABI_FREE(pert_strain) 3034 ABI_FREE(vhxc1_strain) 3035 ABI_FREE(wfk_t_strain) 3036 end if 3037 if (lw_flexo==1.or.lw_flexo==4) then 3038 ABI_FREE(frwfdq) 3039 ABI_FREE(frwfdq_k) 3040 ABI_FREE(isdqwf) 3041 ABI_FREE(isdqwf_k) 3042 ABI_FREE(isdqwf_t1) 3043 ABI_FREE(isdqwf_t1_k) 3044 ABI_FREE(isdqwf_t2) 3045 ABI_FREE(isdqwf_t2_k) 3046 ABI_FREE(isdqwf_t3) 3047 ABI_FREE(isdqwf_t3_k) 3048 ABI_FREE(isdqwf_t4) 3049 ABI_FREE(isdqwf_t4_k) 3050 ABI_FREE(isdqwf_t5) 3051 ABI_FREE(isdqwf_t5_k) 3052 ABI_FREE(isdq_qgradhart) 3053 ABI_FREE(isdq_flg) 3054 end if 3055 ABI_FREE(cg) 3056 ABI_FREE(q1grad) 3057 ABI_FREE(ph1d) 3058 ABI_FREE(indkpt1) 3059 ABI_FREE(istwfk_rbz) 3060 ABI_FREE(kpt_rbz) 3061 ABI_FREE(nband_rbz) 3062 ABI_FREE(occ_rbz) 3063 ABI_FREE(wtk_rbz) 3064 ABI_FREE(kg) 3065 ABI_FREE(kg_k) 3066 ABI_FREE(npwarr) 3067 !ABI_FREE(mpi_enreg%my_kpttab) 3068 !ABI_FREE(mpi_enreg%proc_distrb) 3069 ABI_FREE(irrzon1) 3070 ABI_FREE(phnons1) 3071 ABI_FREE(symrc1) 3072 ABI_FREE(ylm) 3073 ABI_FREE(ylmgr) 3074 ABI_FREE(wfk_t_ddk) 3075 if(xmpi_paral==1) then 3076 ABI_FREE(mpi_enreg%proc_distrb) 3077 end if 3078 3079 ! Clean the header 3080 call hdr0%free() 3081 3082 DBG_EXIT("COLL") 3083 3084 end subroutine dfpt_flexo
ABINIT/dfpt_isdqout [ Functions ]
NAME
dfpt_isdqout
FUNCTION
This subroutine gathers the different terms entering the q-gradient of the internal strain tensor, perfofms the transformation from reduced to cartesian coordinates and writes out the tensor in output files in type-II formulation.
INPUTS
dyewdqdq(2,3,natom,3,3,3)= Second q-gradient of Ewald part of the dynamical matrix
sumed over the second atomic sublattice
frwfdq(2,natpert,3,3,nq1grad)=frozen wf contribution (type-I)
gprimd(3,3)=reciprocal space dimensional primitive translations
isdq_flg(matom,3,nq1grad,3,3)=array that indicates which elements of the q-gradient of
internal strain tensor have been calculated
isdq_qgradhart(2,matom,3,nq1grad,3,3)=electronic electrostatic contribution from the
q-gradient of the Hartree potential
isdqwf(2,matom,3,nq1grad,3,3)=total wave function contributions to the q-gradient of
internal strain tensor (except t4)
isdqwf_t1(2,matom,3,nq1grad,3,3)=term 1 of the wave function contribution
isdqwf_t2(2,matom,3,nq1grad,3,3)=term 2 of the wave function contribution
isdqwf_t3(2,matom,3,nq1grad,3,3)=term 3 of the wave function contribution
isdqwf_t4(2,matom,3,3,3,nq1grad)=term 4 of the wave function contribution (type-I)
isdqwf_t5(2,matom,3,nq1grad,3,3)=term 5 of the wave function contribution
kptopt=2 time reversal symmetry is enforced, 3 trs is not enforced (for debugging purposes)
matom=number of atoms
mpert=maximum number of perturbations
natpert=number of atomic displacement perturbations
nstrpert=number of strain perturbations
nq1grad=number of q1 (q_{\gamma}) gradients
pert_atdis(3,natpert)=array with the info for the atomic displacement perturbations
pert_strain(6,nstrpert)=array with the info for the strain perturbations
prtvol=volume of information to be printed. 1-> The different contributions to the quadrupole are printed.
q1grad(3,nq1grad)=array with the info for the q1 (q_{\gamma}) gradients
rprimd(3,3)=dimensional primitive translations in real space (bohr)
ucvol=Unit cell volume
OUTPUT
d3etot(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE
SIDE EFFECTS
NOTES
SOURCE
4026 subroutine dfpt_isdqout(d3etot,dyewdqdq,frwfdq,gprimd,isdq_flg,isdq_qgradhart,isdqwf,isdqwf_t1,isdqwf_t2,& 4027 & isdqwf_t3,isdqwf_t4,isdqwf_t5,kptopt,matom,mpert,natpert, & 4028 & nstrpert,nq1grad,pert_atdis,pert_strain,prtvol,q1grad,rprimd,ucvol) 4029 4030 !Arguments ------------------------------------ 4031 !scalars 4032 integer,intent(in) :: kptopt,matom,mpert,natpert,nstrpert,nq1grad,prtvol 4033 real(dp),intent(in) :: ucvol 4034 4035 !arrays 4036 integer,intent(in) :: isdq_flg(matom,3,nq1grad,3,3) 4037 integer,intent(in) :: pert_atdis(3,natpert) 4038 integer,intent(in) :: pert_strain(6,nstrpert) 4039 integer,intent(in) :: q1grad(3,nq1grad) 4040 real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert) 4041 real(dp),intent(inout) :: dyewdqdq(2,3,matom,3,3,3) 4042 real(dp),intent(in) :: isdqwf(2,matom,3,nq1grad,3,3) 4043 real(dp),intent(inout) :: frwfdq(2,matom,3,3,3,nq1grad) 4044 real(dp),intent(inout) :: isdq_qgradhart(2,matom,3,nq1grad,3,3) 4045 real(dp),intent(inout) :: isdqwf_t1(2,matom,3,nq1grad,3,3) 4046 real(dp),intent(inout) :: isdqwf_t2(2,matom,3,nq1grad,3,3) 4047 real(dp),intent(inout) :: isdqwf_t3(2,matom,3,nq1grad,3,3) 4048 real(dp),intent(inout) :: isdqwf_t4(2,matom,3,3,3,nq1grad) 4049 real(dp),intent(inout) :: isdqwf_t5(2,matom,3,nq1grad,3,3) 4050 real(dp),intent(in) :: gprimd(3,3) 4051 real(dp),intent(in) :: rprimd(3,3) 4052 4053 !Local variables------------------------------- 4054 !scalars 4055 integer :: alpha,beta,delta,gamma 4056 integer :: iatdir,iatom,iatpert,ibuf,ii,iq1dir,iq1grad,istr1dir,istr2dir,istrpert 4057 integer :: q1pert,strcomp,strpert 4058 integer, parameter :: re=1,im=2 4059 real(dp) :: celastre,celastim,fac,tewim,tewre,tfrim,tfrre,t4im,tmpim,tmpre,t4re 4060 character(len=500) :: msg 4061 4062 !arrays 4063 integer :: flg1(3),flg2(3) 4064 integer, allocatable :: typeI_cartflag(:,:,:,:,:),redflg(:,:,:,:,:) 4065 real(dp) :: vec1(3),vec2(3) 4066 real(dp),allocatable :: dyewdqdq_cart(:,:,:,:,:,:),frwfdq_cart(:,:,:,:,:,:) 4067 real(dp),allocatable :: isdqtens_cart(:,:,:,:,:,:),isdqtens_red(:,:,:,:,:,:) 4068 real(dp),allocatable :: isdqwf_t4_cart(:,:,:,:,:,:) 4069 real(dp),allocatable :: isdqtens_buffer_cart(:,:,:,:,:,:,:) 4070 4071 ! ************************************************************************* 4072 4073 DBG_ENTER("COLL") 4074 4075 !Gather the different terms in the q-gradient of the internal strain tensor 4076 ABI_MALLOC(isdqtens_red,(2,matom,3,nq1grad,3,3)) 4077 ! isdqtens_red=zero 4078 4079 if (kptopt==3) then 4080 4081 !Compute real and 'true' imaginary parts of isdq tensor and independent 4082 !terms. The T4 term and the frozen wf contributions need further treatment 4083 !and they will be lately added to the cartesian coordinates 4084 !version of the isdq tensor 4085 do istrpert=1,nstrpert 4086 istr1dir=pert_strain(3,istrpert) 4087 istr2dir=pert_strain(4,istrpert) 4088 do iq1grad=1,nq1grad 4089 iq1dir=q1grad(2,iq1grad) 4090 do iatpert=1,natpert 4091 iatom=pert_atdis(1,iatpert) 4092 iatdir=pert_atdis(2,iatpert) 4093 4094 if (isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir)==1) then 4095 4096 !The interesting magnitude is 2 times the gradient of the second order Energy 4097 isdqtens_red(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=two* & 4098 & ( isdq_qgradhart(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) + & 4099 & isdqwf(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) ) 4100 isdqtens_red(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=two* & 4101 & ( isdq_qgradhart(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) + & 4102 & isdqwf(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) ) 4103 4104 !Multiply by the imaginary unit that has been factorized out 4105 tmpre=isdqtens_red(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4106 tmpim=isdqtens_red(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4107 isdqtens_red(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-tmpim 4108 isdqtens_red(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=tmpre 4109 4110 !Do the smae for the T4 term 4111 !(This is different from the flexoout routine because the factorized -i in the T4 4112 ! has to be multiplied by -2 as we did for the rest of contributions) 4113 tmpre=isdqwf_t4(re,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4114 tmpim=isdqwf_t4(im,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4115 isdqwf_t4(re,iatom,iatdir,istr1dir,istr2dir,iq1dir)=two*tmpim 4116 isdqwf_t4(im,iatom,iatdir,istr1dir,istr2dir,iq1dir)=-two*tmpre 4117 4118 !Multiply by 2 the frozen wf contribution 4119 tmpre=frwfdq(re,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4120 tmpim=frwfdq(im,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4121 frwfdq(re,iatom,iatdir,istr1dir,istr2dir,iq1dir)=two*tmpre 4122 frwfdq(im,iatom,iatdir,istr1dir,istr2dir,iq1dir)=two*tmpim 4123 4124 !Multiply by 2 the Ewald contribution 4125 tmpre=dyewdqdq(re,iatdir,iatom,istr1dir,istr2dir,iq1dir) 4126 tmpim=dyewdqdq(im,iatdir,iatom,istr1dir,istr2dir,iq1dir) 4127 dyewdqdq(re,iatdir,iatom,istr1dir,istr2dir,iq1dir)=two*tmpre 4128 dyewdqdq(im,iatdir,iatom,istr1dir,istr2dir,iq1dir)=two*tmpim 4129 4130 !Compute and save individual terms in mixed coordinates 4131 if (prtvol>=10) then 4132 4133 tmpre=isdq_qgradhart(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4134 tmpim=isdq_qgradhart(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4135 isdq_qgradhart(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-two*tmpim 4136 isdq_qgradhart(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=two*tmpre 4137 4138 tmpre=isdqwf_t1(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4139 tmpim=isdqwf_t1(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4140 isdqwf_t1(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-two*tmpim 4141 isdqwf_t1(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=two*tmpre 4142 4143 tmpre=isdqwf_t2(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4144 tmpim=isdqwf_t2(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4145 isdqwf_t2(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-two*tmpim 4146 isdqwf_t2(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=two*tmpre 4147 4148 tmpre=isdqwf_t3(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4149 tmpim=isdqwf_t3(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4150 isdqwf_t3(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-two*tmpim 4151 isdqwf_t3(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=two*tmpre 4152 4153 tmpre=isdqwf_t5(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4154 tmpim=isdqwf_t5(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4155 isdqwf_t5(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-two*tmpim 4156 isdqwf_t5(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=two*tmpre 4157 4158 end if 4159 4160 end if 4161 4162 end do 4163 end do 4164 end do 4165 4166 else if (kptopt==2) then 4167 4168 !Compute real part of isdq tensor and independent 4169 !terms. The T4 term and the frozen wf contributions need further treatment 4170 !and they will be lately added to the cartesian coordinates 4171 !version of the isdq tensor 4172 do istrpert=1,nstrpert 4173 istr1dir=pert_strain(3,istrpert) 4174 istr2dir=pert_strain(4,istrpert) 4175 do iq1grad=1,nq1grad 4176 iq1dir=q1grad(2,iq1grad) 4177 do iatpert=1,natpert 4178 iatom=pert_atdis(1,iatpert) 4179 iatdir=pert_atdis(2,iatpert) 4180 4181 if (isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir)==1) then 4182 4183 !The interesting magnitude is 2 times the gradient of the second order Energy 4184 isdqtens_red(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=two* & 4185 & ( isdq_qgradhart(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) + & 4186 & isdqwf(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) ) 4187 isdqtens_red(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=two* & 4188 & ( isdq_qgradhart(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) + & 4189 & isdqwf(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) ) 4190 4191 !Multiply by the imaginary unit that has been factorized out 4192 tmpre=isdqtens_red(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4193 tmpim=isdqtens_red(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4194 isdqtens_red(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-tmpim 4195 isdqtens_red(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=zero 4196 4197 !Do the smae for the T4 term 4198 !(This is different from the flexoout routine because the factorized -i in the T4 4199 ! has to be multiplied by 2 as we did for the rest of contributions) 4200 tmpre=isdqwf_t4(re,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4201 tmpim=isdqwf_t4(im,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4202 isdqwf_t4(re,iatom,iatdir,istr1dir,istr2dir,iq1dir)=two*tmpim 4203 isdqwf_t4(im,iatom,iatdir,istr1dir,istr2dir,iq1dir)=zero 4204 4205 !Multiply by 2 the frozen wf contribution 4206 tmpre=frwfdq(re,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4207 frwfdq(re,iatom,iatdir,istr1dir,istr2dir,iq1dir)=two*tmpre 4208 frwfdq(im,iatom,iatdir,istr1dir,istr2dir,iq1dir)=zero 4209 4210 !Multiply by 2 the Ewald contribution 4211 tmpre=dyewdqdq(re,iatdir,iatom,istr1dir,istr2dir,iq1dir) 4212 dyewdqdq(re,iatdir,iatom,istr1dir,istr2dir,iq1dir)=two*tmpre 4213 dyewdqdq(im,iatdir,iatom,istr1dir,istr2dir,iq1dir)=zero 4214 4215 !Compute and save individual terms in mixed coordinates 4216 if (prtvol>=10) then 4217 4218 tmpim=isdq_qgradhart(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4219 isdq_qgradhart(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-two*tmpim 4220 isdq_qgradhart(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=zero 4221 4222 tmpim=isdqwf_t1(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4223 isdqwf_t1(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-two*tmpim 4224 isdqwf_t1(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=zero 4225 4226 tmpim=isdqwf_t2(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4227 isdqwf_t2(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-two*tmpim 4228 isdqwf_t2(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=zero 4229 4230 tmpim=isdqwf_t3(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4231 isdqwf_t3(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-two*tmpim 4232 isdqwf_t3(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=zero 4233 4234 tmpim=isdqwf_t5(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4235 isdqwf_t5(re,iatom,iatdir,iq1dir,istr1dir,istr2dir)=-two*tmpim 4236 isdqwf_t5(im,iatom,iatdir,iq1dir,istr1dir,istr2dir)=zero 4237 4238 end if 4239 4240 end if 4241 4242 end do 4243 end do 4244 end do 4245 4246 else 4247 write(msg,"(1a)") 'kptopt must be 2 or 3 for long-wave DFPT calculations' 4248 ABI_BUG(msg) 4249 end if 4250 4251 !Transform to complete cartesian coordinates all the contributions 4252 ABI_MALLOC(isdqtens_cart,(2,matom,3,nq1grad,3,3)) 4253 ABI_MALLOC(isdqwf_t4_cart,(2,matom,3,3,3,nq1grad)) 4254 ABI_MALLOC(frwfdq_cart,(2,matom,3,3,3,nq1grad)) 4255 ABI_MALLOC(dyewdqdq_cart,(2,3,matom,3,3,nq1grad)) 4256 ABI_MALLOC(typeI_cartflag,(matom,3,3,3,nq1grad)) 4257 isdqtens_cart=isdqtens_red 4258 isdqwf_t4_cart=isdqwf_t4 4259 frwfdq_cart=frwfdq 4260 dyewdqdq_cart=dyewdqdq 4261 typeI_cartflag=0 4262 4263 if (prtvol>=10) then 4264 ABI_MALLOC(isdqtens_buffer_cart,(5,2,matom,3,nq1grad,3,3)) 4265 isdqtens_buffer_cart(1,:,:,:,:,:,:)=isdqwf_t1(:,:,:,:,:,:) 4266 isdqtens_buffer_cart(2,:,:,:,:,:,:)=isdqwf_t2(:,:,:,:,:,:) 4267 isdqtens_buffer_cart(3,:,:,:,:,:,:)=isdqwf_t3(:,:,:,:,:,:) 4268 isdqtens_buffer_cart(4,:,:,:,:,:,:)=isdq_qgradhart(:,:,:,:,:,:) 4269 isdqtens_buffer_cart(5,:,:,:,:,:,:)=isdqwf_t5(:,:,:,:,:,:) 4270 end if 4271 4272 !##1st transform coordinates of the atomic displacement derivative contributions 4273 do istr2dir=1,3 4274 do istr1dir=1,3 4275 do iq1dir=1,3 4276 do ii=1,2 4277 do iatom=1,matom 4278 4279 do iatdir=1,3 4280 vec1(iatdir)=isdqtens_cart(ii,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4281 flg1(iatdir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4282 end do 4283 call cart39(flg1,flg2,gprimd,iatom,matom,rprimd,vec1,vec2) 4284 do iatdir=1,3 4285 isdqtens_cart(ii,iatom,iatdir,iq1dir,istr1dir,istr2dir)=vec2(iatdir) 4286 end do 4287 4288 do iatdir=1,3 4289 vec1(iatdir)=isdqwf_t4_cart(ii,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4290 flg1(iatdir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4291 end do 4292 call cart39(flg1,flg2,gprimd,iatom,matom,rprimd,vec1,vec2) 4293 do iatdir=1,3 4294 isdqwf_t4_cart(ii,iatom,iatdir,istr1dir,istr2dir,iq1dir)=vec2(iatdir) 4295 typeI_cartflag(iatom,iatdir,istr1dir,istr2dir,iq1dir)=flg2(iatdir) 4296 end do 4297 4298 do iatdir=1,3 4299 vec1(iatdir)=frwfdq_cart(ii,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4300 flg1(iatdir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4301 end do 4302 call cart39(flg1,flg2,gprimd,iatom,matom,rprimd,vec1,vec2) 4303 do iatdir=1,3 4304 frwfdq_cart(ii,iatom,iatdir,istr1dir,istr2dir,iq1dir)=vec2(iatdir) 4305 end do 4306 4307 do iatdir=1,3 4308 vec1(iatdir)=dyewdqdq_cart(ii,iatdir,iatom,istr1dir,istr2dir,iq1dir) 4309 flg1(iatdir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4310 end do 4311 call cart39(flg1,flg2,gprimd,iatom,matom,rprimd,vec1,vec2) 4312 do iatdir=1,3 4313 dyewdqdq_cart(ii,iatdir,iatom,istr1dir,istr2dir,iq1dir)=vec2(iatdir) 4314 end do 4315 4316 end do 4317 end do 4318 end do 4319 end do 4320 end do 4321 4322 !For debugging, transform also all the individual contributions 4323 if (prtvol>=10) then 4324 do ibuf=1,5 4325 do istr2dir=1,3 4326 do istr1dir=1,3 4327 do iq1dir=1,3 4328 do ii=1,2 4329 do iatom=1,matom 4330 4331 do iatdir=1,3 4332 vec1(iatdir)=isdqtens_buffer_cart(ibuf,ii,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4333 flg1(iatdir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4334 end do 4335 call cart39(flg1,flg2,gprimd,iatom,matom,rprimd,vec1,vec2) 4336 do iatdir=1,3 4337 isdqtens_buffer_cart(ibuf,ii,iatom,iatdir,iq1dir,istr1dir,istr2dir)=vec2(iatdir) 4338 end do 4339 4340 end do 4341 end do 4342 end do 4343 end do 4344 end do 4345 end do 4346 end if 4347 4348 !##2nd transform coordinates of the q-gradient (treat it as electric field) 4349 do istr2dir=1,3 4350 do istr1dir=1,3 4351 do iatom=1,matom 4352 do iatdir=1,3 4353 do ii=1,2 4354 4355 do iq1dir=1,3 4356 vec1(iq1dir)=isdqtens_cart(ii,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4357 flg1(iq1dir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4358 end do 4359 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 4360 do iq1dir=1,3 4361 isdqtens_cart(ii,iatom,iatdir,iq1dir,istr1dir,istr2dir)=vec2(iq1dir) 4362 end do 4363 4364 do iq1dir=1,3 4365 vec1(iq1dir)=frwfdq_cart(ii,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4366 flg1(iq1dir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4367 end do 4368 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 4369 do iq1dir=1,3 4370 frwfdq_cart(ii,iatom,iatdir,istr1dir,istr2dir,iq1dir)=vec2(iq1dir) 4371 end do 4372 4373 do iq1dir=1,3 4374 vec1(iq1dir)=dyewdqdq_cart(ii,iatdir,iatom,istr1dir,istr2dir,iq1dir) 4375 flg1(iq1dir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4376 end do 4377 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 4378 do iq1dir=1,3 4379 dyewdqdq_cart(ii,iatdir,iatom,istr1dir,istr2dir,iq1dir)=vec2(iq1dir) 4380 end do 4381 4382 end do 4383 end do 4384 end do 4385 end do 4386 end do 4387 4388 !For debugging, transform also all the individual contributions 4389 if (prtvol>=10) then 4390 do ibuf=1,5 4391 do istr2dir=1,3 4392 do istr1dir=1,3 4393 do iatom=1,matom 4394 do iatdir=1,3 4395 do ii=1,2 4396 4397 do iq1dir=1,3 4398 vec1(iq1dir)=isdqtens_buffer_cart(ibuf,ii,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4399 flg1(iq1dir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4400 end do 4401 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 4402 do iq1dir=1,3 4403 isdqtens_buffer_cart(ibuf,ii,iatom,iatdir,iq1dir,istr1dir,istr2dir)=vec2(iq1dir) 4404 end do 4405 4406 end do 4407 end do 4408 end do 4409 end do 4410 end do 4411 end do 4412 endif 4413 4414 !3rd## Transform the metric perturbation direction to cartesian coordinates in 4415 !the frozen wf contributions (treat it as an atomic displacement) 4416 do istr2dir=1,3 4417 do iq1dir=1,3 4418 do iatom=1,matom 4419 do iatdir=1,3 4420 do ii=1,2 4421 4422 do istr1dir=1,3 4423 vec1(istr1dir)=frwfdq_cart(ii,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4424 flg1(istr1dir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4425 end do 4426 call cart39(flg1,flg2,gprimd,iatom,matom,rprimd,vec1,vec2) 4427 do istr1dir=1,3 4428 frwfdq_cart(ii,iatom,iatdir,istr1dir,istr2dir,iq1dir)=vec2(istr1dir) 4429 end do 4430 4431 do istr1dir=1,3 4432 vec1(istr1dir)=dyewdqdq_cart(ii,iatdir,iatom,istr1dir,istr2dir,iq1dir) 4433 flg1(istr1dir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4434 end do 4435 call cart39(flg1,flg2,gprimd,iatom,matom,rprimd,vec1,vec2) 4436 do istr1dir=1,3 4437 dyewdqdq_cart(ii,iatdir,iatom,istr1dir,istr2dir,iq1dir)=vec2(istr1dir) 4438 end do 4439 4440 end do 4441 end do 4442 end do 4443 end do 4444 end do 4445 4446 !4th## Transform the second q-gradient direction to cartesian coordinates in 4447 !the frozen wf contributions (treat it as an electric field) 4448 do istr1dir=1,3 4449 do iq1dir=1,3 4450 do iatom=1,matom 4451 do iatdir=1,3 4452 do ii=1,2 4453 4454 do istr2dir=1,3 4455 vec1(istr2dir)=frwfdq_cart(ii,iatom,iatdir,istr1dir,istr2dir,iq1dir) 4456 flg1(istr2dir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4457 end do 4458 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 4459 do istr2dir=1,3 4460 frwfdq_cart(ii,iatom,iatdir,istr1dir,istr2dir,iq1dir)=vec2(istr2dir) 4461 end do 4462 4463 do istr2dir=1,3 4464 vec1(istr2dir)=dyewdqdq_cart(ii,iatdir,iatom,istr1dir,istr2dir,iq1dir) 4465 flg1(istr2dir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4466 end do 4467 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 4468 do istr2dir=1,3 4469 dyewdqdq_cart(ii,iatdir,iatom,istr1dir,istr2dir,iq1dir)=vec2(istr2dir) 4470 end do 4471 4472 end do 4473 end do 4474 end do 4475 end do 4476 end do 4477 4478 !Write the q-gradient fo the internal strain tensor in cartesian coordinates 4479 !Open output files 4480 if (prtvol>=10) then 4481 open(unit=71,file='isdq_wf_t1.out',status='unknown',form='formatted',action='write') 4482 open(unit=72,file='isdq_wf_t2.out',status='unknown',form='formatted',action='write') 4483 open(unit=73,file='isdq_wf_t3.out',status='unknown',form='formatted',action='write') 4484 open(unit=74,file='isdq_wf_t4.out',status='unknown',form='formatted',action='write') 4485 open(unit=75,file='isdq_wf_t5.out',status='unknown',form='formatted',action='write') 4486 open(unit=76,file='isdq_elecstic.out',status='unknown',form='formatted',action='write') 4487 open(unit=77,file='isdq_frwf.out',status='unknown',form='formatted',action='write') 4488 open(unit=78,file='isdq_ewdqdq.out',status='unknown',form='formatted',action='write') 4489 end if 4490 4491 write(ab_out,'(a)')' ' 4492 write(ab_out,'(a)')' First real-space moment of piezoelectric force-response tensor, in cartesian coordinates,' 4493 write(ab_out,'(a)')'atom atdir qgrdir strdir1 strdir2 real part imaginary part' 4494 do istr2dir=1,3 4495 delta=istr2dir 4496 do istr1dir=1,3 4497 beta=istr1dir 4498 do iq1dir=1,3 4499 gamma=iq1dir 4500 do iatom=1,matom 4501 do iatdir=1,3 4502 alpha=iatdir 4503 4504 if (typeI_cartflag(iatom,alpha,beta,delta,gamma)==1 .and. & 4505 & typeI_cartflag(iatom,alpha,delta,gamma,beta)==1 .and. & 4506 & typeI_cartflag(iatom,alpha,gamma,beta,delta)==1 ) then 4507 4508 !Converts the T4 term to type-II form 4509 t4re= isdqwf_t4_cart(re,iatom,alpha,beta,delta,gamma) + & 4510 & isdqwf_t4_cart(re,iatom,alpha,delta,gamma,beta) - & 4511 & isdqwf_t4_cart(re,iatom,alpha,gamma,beta,delta) 4512 t4im= isdqwf_t4_cart(im,iatom,alpha,beta,delta,gamma) + & 4513 & isdqwf_t4_cart(im,iatom,alpha,delta,gamma,beta) - & 4514 & isdqwf_t4_cart(im,iatom,alpha,gamma,beta,delta) 4515 4516 !Converts the frozen wf term to type-II form 4517 tfrre= frwfdq_cart(re,iatom,alpha,beta,delta,gamma) + & 4518 & frwfdq_cart(re,iatom,alpha,delta,gamma,beta) - & 4519 & frwfdq_cart(re,iatom,alpha,gamma,beta,delta) 4520 tfrim= frwfdq_cart(im,iatom,alpha,beta,delta,gamma) + & 4521 & frwfdq_cart(im,iatom,alpha,delta,gamma,beta) - & 4522 & frwfdq_cart(im,iatom,alpha,gamma,beta,delta) 4523 4524 !Converts the Ewald term to type-II form 4525 tewre= dyewdqdq_cart(re,alpha,iatom,beta,delta,gamma) + & 4526 & dyewdqdq_cart(re,alpha,iatom,delta,gamma,beta) - & 4527 & dyewdqdq_cart(re,alpha,iatom,gamma,beta,delta) 4528 4529 tewim= dyewdqdq_cart(im,alpha,iatom,beta,delta,gamma) + & 4530 & dyewdqdq_cart(im,alpha,iatom,delta,gamma,beta) - & 4531 & dyewdqdq_cart(im,alpha,iatom,gamma,beta,delta) 4532 4533 !Add the type-II T4, frozen wf and Ewald contributions 4534 isdqtens_cart(re,iatom,alpha,gamma,beta,delta)= & 4535 & isdqtens_cart(re,iatom,alpha,gamma,beta,delta) + t4re + & 4536 & half*tfrre + quarter*tewre 4537 isdqtens_cart(im,iatom,alpha,gamma,beta,delta)= & 4538 & isdqtens_cart(im,iatom,alpha,gamma,beta,delta) + t4im + & 4539 & half*tfrim + quarter*tewim 4540 4541 !Writes the complete q-gradient of internal strain tensor 4542 write(ab_out,'(5(i5,3x),2(1x,f20.10))') iatom,alpha,gamma,beta,delta, & 4543 & isdqtens_cart(re,iatom,alpha,gamma,beta,delta), & 4544 & isdqtens_cart(im,iatom,alpha,gamma,beta,delta) 4545 4546 if (prtvol>=10) then 4547 write(71,'(5(i5,3x),2(1x,f20.10))') iatom,alpha,gamma,beta,delta, & 4548 & isdqtens_buffer_cart(1,re,iatom,alpha,gamma,beta,delta), & 4549 & isdqtens_buffer_cart(1,im,iatom,alpha,gamma,beta,delta) 4550 4551 write(72,'(5(i5,3x),2(1x,f20.10))') iatom,alpha,gamma,beta,delta, & 4552 & isdqtens_buffer_cart(2,re,iatom,alpha,gamma,beta,delta), & 4553 & isdqtens_buffer_cart(2,im,iatom,alpha,gamma,beta,delta) 4554 4555 write(73,'(5(i5,3x),2(1x,f20.10))') iatom,alpha,gamma,beta,delta, & 4556 & isdqtens_buffer_cart(3,re,iatom,alpha,gamma,beta,delta), & 4557 & isdqtens_buffer_cart(3,im,iatom,alpha,gamma,beta,delta) 4558 4559 write(74,'(5(i5,3x),2(1x,f20.10))') iatom,alpha,gamma,beta,delta,t4re,t4im 4560 4561 write(75,'(5(i5,3x),2(1x,f20.10))') iatom,alpha,gamma,beta,delta, & 4562 & isdqtens_buffer_cart(5,re,iatom,alpha,gamma,beta,delta), & 4563 & isdqtens_buffer_cart(5,im,iatom,alpha,gamma,beta,delta) 4564 4565 write(76,'(5(i5,3x),2(1x,f20.10))') iatom,alpha,gamma,beta,delta, & 4566 & isdqtens_buffer_cart(4,re,iatom,alpha,gamma,beta,delta), & 4567 & isdqtens_buffer_cart(4,im,iatom,alpha,gamma,beta,delta) 4568 4569 write(77,'(5(i5,3x),2(1x,f20.10))') iatom,alpha,gamma,beta,delta, & 4570 & half*tfrre, half*tfrim 4571 4572 write(78,'(5(i5,3x),2(1x,f20.10))') iatom,alpha,gamma,beta,delta, & 4573 & quarter*tewre, quarter*tewim 4574 4575 end if 4576 4577 end if 4578 4579 end do 4580 end do 4581 write(ab_out,'(a)')' ' 4582 if (prtvol>=10) then 4583 write(71,'(a)')' ' 4584 write(72,'(a)')' ' 4585 write(73,'(a)')' ' 4586 write(74,'(a)')' ' 4587 write(75,'(a)')' ' 4588 write(76,'(a)')' ' 4589 write(77,'(a)')' ' 4590 write(78,'(a)')' ' 4591 end if 4592 end do 4593 end do 4594 end do 4595 4596 if (prtvol>=10) then 4597 close(71) 4598 close(72) 4599 close(73) 4600 close(74) 4601 close(75) 4602 close(76) 4603 close(77) 4604 close(78) 4605 end if 4606 4607 !Write the elastic tensor 4608 write(ab_out,'(a)')' ' 4609 write(ab_out,'(a)')' Clamped-ion elastic tensor, in cartesian coordinates ' 4610 write(ab_out,'(a)')' (from q-gradient of piezofrt),' 4611 write(ab_out,'(a)')' (for stressed cells it lacks an improper contribution),' 4612 write(ab_out,'(a)')' atdir qgrdir strdir1 strdir2 real part imaginary part' 4613 do istr2dir=1,3 4614 delta=istr2dir 4615 do istr1dir=1,3 4616 beta=istr1dir 4617 do iq1dir=1,3 4618 gamma=iq1dir 4619 do iatdir=1,3 4620 alpha=iatdir 4621 celastre=zero 4622 celastim=zero 4623 do iatom=1,matom 4624 4625 if (isdq_flg(iatom,alpha,gamma,beta,delta)==1) then 4626 celastre=celastre+isdqtens_cart(re,iatom,alpha,gamma,beta,delta) 4627 celastim=celastim+isdqtens_cart(im,iatom,alpha,gamma,beta,delta) 4628 end if 4629 4630 end do 4631 4632 write(ab_out,'(4(i5,4x),2(1x,f20.10))') alpha,gamma,beta,delta, & 4633 & celastre/ucvol, celastim/ucvol 4634 4635 end do 4636 end do 4637 write(ab_out,'(a)')' ' 4638 end do 4639 end do 4640 write(ab_out,'(80a)')('=',ii=1,80) 4641 4642 !Calculate the contribution to the d3etot in mixed (reduced/cartesian) coordinates 4643 isdqtens_red=isdqtens_cart 4644 ABI_FREE(isdqtens_cart) 4645 ABI_FREE(isdqwf_t4_cart) 4646 ABI_FREE(frwfdq_cart) 4647 ABI_FREE(dyewdqdq_cart) 4648 ABI_FREE(typeI_cartflag) 4649 ABI_MALLOC(redflg,(matom,3,3,3,3)) 4650 4651 !1st transform back coordinates of the atomic displacement derivative 4652 do istr2dir=1,3 4653 do istr1dir=1,3 4654 do iq1dir=1,3 4655 do ii=1,2 4656 do iatom=1,matom 4657 do iatdir=1,3 4658 vec1(iatdir)=isdqtens_red(ii,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4659 flg1(iatdir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4660 end do 4661 call cart39(flg1,flg2,transpose(rprimd),iatom,matom,transpose(gprimd),vec1,vec2) 4662 do iatdir=1,3 4663 isdqtens_red(ii,iatom,iatdir,iq1dir,istr1dir,istr2dir)=vec2(iatdir) 4664 redflg(iatom,iatdir,iq1dir,istr1dir,istr2dir)=flg2(iatdir) 4665 end do 4666 end do 4667 end do 4668 end do 4669 end do 4670 end do 4671 4672 !2nd transform back coordinates of the q-gradient (treat it as electric field) 4673 fac=two_pi ** 2 4674 do istr2dir=1,3 4675 do istr1dir=1,3 4676 do iatom=1,matom 4677 do iatdir=1,3 4678 do ii=1,2 4679 do iq1dir=1,3 4680 vec1(iq1dir)=isdqtens_red(ii,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4681 flg1(iq1dir)=isdq_flg(iatom,iatdir,iq1dir,istr1dir,istr2dir) 4682 end do 4683 call cart39(flg1,flg2,transpose(rprimd),matom+2,matom,transpose(gprimd),vec1,vec2) 4684 do iq1dir=1,3 4685 isdqtens_red(ii,iatom,iatdir,iq1dir,istr1dir,istr2dir)=vec2(iq1dir)*fac 4686 redflg(iatom,iatdir,iq1dir,istr1dir,istr2dir)=flg2(iq1dir) 4687 end do 4688 end do 4689 end do 4690 end do 4691 end do 4692 end do 4693 4694 !Add contributions to d3etot (remove the -2 factor) 4695 q1pert=matom+8 4696 do istrpert=1,nstrpert 4697 strpert=pert_strain(1,istrpert) 4698 strcomp=pert_strain(2,istrpert) 4699 istr1dir=pert_strain(3,istrpert) 4700 istr2dir=pert_strain(4,istrpert) 4701 do iq1grad=1,nq1grad 4702 iq1dir=q1grad(2,iq1grad) 4703 do iatom=1,matom 4704 do iatdir=1,3 4705 4706 if (redflg(iatom,iatdir,iq1dir,istr1dir,istr2dir)==1) then 4707 d3etot(re,iatdir,iatom,strcomp,strpert,iq1dir,q1pert)= & 4708 & -half*isdqtens_red(re,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4709 d3etot(im,iatdir,iatom,strcomp,strpert,iq1dir,q1pert)= & 4710 & -half*isdqtens_red(im,iatom,iatdir,iq1dir,istr1dir,istr2dir) 4711 end if 4712 4713 end do 4714 end do 4715 end do 4716 end do 4717 4718 ABI_FREE(isdqtens_red) 4719 ABI_FREE(redflg) 4720 4721 DBG_EXIT("COLL") 4722 end subroutine dfpt_isdqout
ABINIT/dfpt_qdrpole [ Functions ]
NAME
dfpt_qdrpole
FUNCTION
This routine computes the elements of the quadrupole and the P^(1) tensor as the second order q-gradient (d2/dq1dq2) of the charge response to an atomic displacement (see, e.g., M.Royo and M. Stengel to be published). The atoms and atomic displacements directions that are evaluated are fixed by the rfatpol and rfdir variables in the input, rfdir also fixes the directions for the dq2 derivative, whereas dq1 is evaluated in all three directions.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) codvsn=code version doccde(mband*nkpt*nsppol)=derivative of occupancies wrt the energy dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset gmet(3,3)=reciprocal space metric tensor in bohr**-2. gprimd(3,3)=reciprocal space dimensional primitive translations kxc(nfft,nkxc)=exchange and correlation kernel mpert=maximum number of perturbations for output processing mpi_enreg=information about MPI parallelization nattyp(ntypat)= # atoms of each type. nfft=(effective) number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other nkpt=number of k points in the full BZ nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed. nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS (DUMMY) pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data (DUMMY) pertsy(3,natom+6)=set of perturbations that form a basis for all other perturbations psps <type(pseudopotential_type)>=variables related to pseudopotentials rmet(3,3)=real space metric (bohr**2) rprimd(3,3)=dimensional primitive translations in real space (bohr) rhog(2,nfftf)=array for Fourier transform of GS electron density rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3. timrev=1 if time-reversal preserves the q wavevector; 0 otherwise. ucvol=unit cell volume in bohr**3. xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
blkflg(3,mpert,3,mpert,3,mpert)= ( 1 if the element of the 3dte has been calculated ; 0 otherwise ) d3etot(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE
SIDE EFFECTS
NOTES
For historical reasons the Quadrupole and P^(1) tensor has been obtained from
the q-gradient of E^{\tau_{\kappa\beta}* \epsilon_{alpha}}. This is
different from the equations appearing in PRX 9, 021050 (2019) which
stand for the q-gradient of E^{\epsilon_{alpha}* \tau_{\kappa\beta}}.
By this reason there is an additional -1 factor here when computing
the physical observables.
SOURCE
138 subroutine dfpt_qdrpole(atindx,blkflg,codvsn,d3etot,doccde,dtfil,dtset,& 139 & gmet,gprimd,kxc,mpert,& 140 & mpi_enreg,nattyp,nfft,ngfft,nkpt,nkxc,& 141 & nspden,nsppol,occ,pawrhoij,pawtab,pertsy,psps,rmet,rprimd,rhog,rhor,& 142 & timrev,ucvol,xred) 143 144 !Arguments ------------------------------------ 145 !scalars 146 integer,intent(in) :: mpert,nfft,nkpt,nkxc,nspden,nsppol,timrev 147 integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) 148 real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert) 149 real(dp),intent(in) :: ucvol 150 character(len=8), intent(in) :: codvsn 151 type(MPI_type),intent(inout) :: mpi_enreg 152 type(datafiles_type),intent(in) :: dtfil 153 type(dataset_type),intent(in) :: dtset 154 type(pseudopotential_type),intent(in) :: psps 155 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 156 type(pawrhoij_type),intent(inout) :: pawrhoij(:) 157 !arrays 158 integer,intent(in) :: atindx(dtset%natom) 159 integer,intent(in) :: nattyp(dtset%ntypat),ngfft(18) 160 integer,intent(in) :: pertsy(3,dtset%natom+6) 161 real(dp),intent(in) :: doccde(dtset%mband*nkpt*dtset%nsppol) 162 real(dp),intent(in) :: gmet(3,3), gprimd(3,3) 163 real(dp),intent(in) :: kxc(nfft,nkxc) 164 real(dp),intent(in) :: occ(dtset%mband*nkpt*dtset%nsppol) 165 real(dp),intent(in) :: rmet(3,3),rprimd(3,3) 166 real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,nspden) 167 real(dp),intent(in) :: xred(3,dtset%natom) 168 169 !Local variables------------------------------- 170 !scalars 171 integer :: ask_accurate,bdtot_index,bdtot1_index,bantot_rbz 172 integer :: cplex,formeig,forunit,gscase,iatpert,iatpert_cnt,iatpol 173 integer :: iatdir,icg,ierr,ii,ikg,ikpt,ikpt1,ilm,qcar,iq1dir,iq1grad,iq1grad_cnt 174 integer :: iq1q2grad,iq1q2grad_var,iq2dir,iq2grad,iq2grad_cnt,ireadwf0,isppol,istwf_k 175 integer :: jj,master,matom,matpert,mcg,me,mgfft 176 integer :: mkmem_rbz,mpw,my_nkpt_rbz 177 integer :: natpert,nband_k,nfftot,nhat1grdim,nkpt_rbz,npw_k,npw1_k 178 integer :: nq1grad,nq1q2grad,nq2grad,nsym1,nylmgr,n3xccc 179 integer :: optorth,optene,option,optres 180 integer :: pawread,pertcase,qdir,spaceworld 181 integer :: usexcnhat,useylmgr 182 integer :: mband_mem 183 integer,parameter :: formeig1=1 184 integer,parameter :: re=1,im=2 185 real(dp) :: boxcut,doti,dotr,dum_scl,ecut_eff,ecut,etotal,fermie,gsqcut,residm 186 real(dp) :: vres2, wtk_k 187 logical :: non_magnetic_xc 188 character(len=500) :: msg 189 character(len=fnlen) :: fi1o,fiwfatdis,fiwfefield,fiwfddk,fiwfdkdk 190 type(hdr_type) :: hdr_den 191 type(ebands_t) :: bs_rbz 192 type(hdr_type) :: hdr0 193 type(wvl_data) :: wvl 194 type(wffile_type) :: wffgs,wfftgs 195 type(gs_hamiltonian_type) :: gs_hamkq 196 197 !arrays 198 integer,allocatable :: bz2ibz_smap(:,:) 199 integer,allocatable :: indkpt1(:), indkpt1_tmp(:),indsy1(:,:,:),irrzon1(:,:,:) 200 integer,allocatable :: istwfk_rbz(:),kg(:,:),kg_k(:,:) 201 integer,allocatable :: nband_rbz(:),npwarr(:),npwtot(:) 202 integer,allocatable :: pert_atdis(:,:), pert_atdis_tmp(:,:) 203 integer,allocatable :: q1grad(:,:),q1grad_tmp(:,:),q1q2grad(:,:),q2grad(:,:),q2grad_tmp(:,:) 204 integer,allocatable :: qdrflg(:,:,:,:) 205 integer,allocatable :: symaf1(:),symrc1(:,:,:),symrl1(:,:,:) 206 real(dp) :: kpoint(3) 207 real(dp),allocatable :: cg(:,:),doccde_rbz(:) 208 real(dp),allocatable :: eigen0(:),eqgradhart(:,:,:,:) 209 real(dp),allocatable :: kpt_rbz(:,:) 210 real(dp),allocatable :: nhat(:,:),nhat1(:,:),nhat1gr(:,:,:) 211 real(dp),allocatable :: occ_k(:),occ_rbz(:) 212 real(dp),allocatable :: ph1d(:,:),phnons1(:,:,:) 213 real(dp),allocatable :: qdrpwf(:,:,:,:),qdrpwf_k(:,:,:,:) 214 real(dp),allocatable :: qdrpwf_t1(:,:,:,:),qdrpwf_t1_k(:,:,:,:) 215 real(dp),allocatable :: qdrpwf_t2(:,:,:,:),qdrpwf_t2_k(:,:,:,:) 216 real(dp),allocatable :: qdrpwf_t3(:,:,:,:),qdrpwf_t3_k(:,:,:,:) 217 real(dp),allocatable :: qdrpwf_t4(:,:,:,:),qdrpwf_t4_k(:,:,:,:) 218 real(dp),allocatable :: qdrpwf_t5(:,:,:,:),qdrpwf_t5_k(:,:,:,:) 219 real(dp),allocatable :: rhog1_atdis(:,:,:) 220 real(dp),allocatable :: rhog1_tmp(:,:) 221 real(dp),allocatable :: rhor1_atdis(:,:,:),rhor1_efield(:,:,:) 222 real(dp),allocatable :: rhor1_cplx(:,:),rhor1_tmp(:,:),rhor1_real(:,:) 223 real(dp),allocatable :: tnons1(:,:) 224 real(dp),allocatable :: vhartr1(:),vhxc1_atdis(:,:),vhxc1_efield(:,:) 225 real(dp),allocatable :: vpsp1(:),vqgradhart(:),vresid1(:,:),vxc1(:,:) 226 real(dp),allocatable :: dum_vxc(:,:),vxc1dq(:,:),vxc1dqc(:,:,:,:) 227 real(dp),allocatable :: wtk_folded(:), wtk_rbz(:) 228 real(dp),allocatable,target :: vtrial1(:,:) 229 real(dp),allocatable :: xccc3d1(:) 230 real(dp),allocatable :: ylm(:,:),ylm_k(:,:),ylmgr(:,:,:),ylmgr_k(:,:,:) 231 type(pawrhoij_type),allocatable :: pawrhoij_read(:) 232 type(wfk_t),allocatable :: wfk_t_atdis(:),wfk_t_efield(:),wfk_t_ddk(:),wfk_t_dkdk(:) 233 234 ! ************************************************************************* 235 236 DBG_ENTER("COLL") 237 238 !Anounce start of quadrupole tensor calculation 239 write(msg, '(a,80a,a,a,a)' ) ch10,('=',ii=1,80),ch10,& 240 & ' ==> Compute Quadrupole Tensor <== ',ch10 241 call wrtout(std_out,msg,'COLL') 242 call wrtout(ab_out,msg,'COLL') 243 244 !Keep track of total time spent in dfpt_qdrpole 245 !!To be implemented if required 246 247 !Init parallelism 248 spaceworld=mpi_enreg%comm_cell 249 me=mpi_enreg%me_kpt 250 master=0 251 252 !Get FFT grid(s) sizes (be careful !) See NOTES in the comments at the beginning of respfn.F90 253 mgfft=dtset%mgfft 254 ecut=dtset%ecut 255 ecut_eff=ecut*(dtset%dilatmx)**2 256 257 !Compute large sphere cut-off gsqcut 258 call getcut(boxcut,ecut,gmet,gsqcut,dtset%iboxcut,std_out,dtset%qptn,ngfft) 259 260 !Various initializations 261 cplex=2-timrev 262 matom=dtset%natom 263 usexcnhat=0 264 n3xccc=0 265 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 266 non_magnetic_xc=.true. 267 268 269 !Generate the 1-dimensional phases 270 ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom)) 271 call getph(atindx,dtset%natom,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),ph1d,xred) 272 273 !################# PERTURBATIONS AND q-GRADIENTS LABELLING ############################ 274 275 !Determine which atomic displacements and q-gradient directions have to be evaluated 276 !taking into account the perturbation symmetries 277 matpert=dtset%natom*3 278 ABI_MALLOC(pert_atdis_tmp,(3,matpert)) 279 pert_atdis_tmp=0 280 iatpert_cnt=0 281 do iatpol=1,matom 282 do iatdir=1,3 283 if (pertsy(iatdir,iatpol)==1) then 284 iatpert_cnt=iatpert_cnt+1 285 pert_atdis_tmp(1,iatpert_cnt)=iatpol !atom displaced 286 pert_atdis_tmp(2,iatpert_cnt)=iatdir !direction of displacement 287 pert_atdis_tmp(3,iatpert_cnt)=(iatpol-1)*3+iatdir !like pertcase in dfpt_loopert.f90 288 end if 289 end do 290 end do 291 natpert=iatpert_cnt 292 ABI_MALLOC(pert_atdis,(3,natpert)) 293 do iatpert=1,natpert 294 pert_atdis(:,iatpert)=pert_atdis_tmp(:,iatpert) 295 end do 296 ABI_FREE(pert_atdis_tmp) 297 298 !The q2grad is related with the response to the electric field 299 ABI_MALLOC(q2grad_tmp,(3,3)) 300 q2grad_tmp=0 301 iq2grad_cnt=0 302 do iq2grad=1,3 303 if (pertsy(iq2grad,dtset%natom+2)==1) then 304 iq2grad_cnt=iq2grad_cnt+1 305 q2grad_tmp(1,iq2grad_cnt)=dtset%natom+2 !electric field pert 306 q2grad_tmp(2,iq2grad_cnt)=iq2grad !electric field direction 307 q2grad_tmp(3,iq2grad_cnt)=matpert+3+iq2grad !like pertcase in dfpt_loopert.f90 308 end if 309 end do 310 nq2grad=iq2grad_cnt 311 ABI_MALLOC(q2grad,(3,nq2grad)) 312 do iq2grad=1,nq2grad 313 q2grad(:,iq2grad)=q2grad_tmp(:,iq2grad) 314 end do 315 ABI_FREE(q2grad_tmp) 316 317 !The q1grad is related with the response to the ddk 318 ABI_MALLOC(q1grad_tmp,(3,3)) 319 q1grad_tmp=0 320 iq1grad_cnt=0 321 do iq1grad=1,3 322 if (pertsy(iq1grad,dtset%natom+1)==1) then 323 iq1grad_cnt=iq1grad_cnt+1 324 q1grad_tmp(1,iq1grad_cnt)=dtset%natom+1 !ddk perturbation 325 q1grad_tmp(2,iq1grad_cnt)=iq1grad !ddk or ddq direction 326 q1grad_tmp(3,iq1grad_cnt)=matpert+iq1grad !like pertcase in dfpt_loopert.f90 327 end if 328 end do 329 nq1grad=iq1grad_cnt 330 ABI_MALLOC(q1grad,(3,nq1grad)) 331 do iq1grad=1,nq1grad 332 q1grad(:,iq1grad)=q1grad_tmp(:,iq1grad) 333 end do 334 ABI_FREE(q1grad_tmp) 335 336 !For the evaluation of the 2nd order q-gradient, the 9 directios are activated because 337 !currently the program calculates by defect all the components of the d2_dkdk perturbation. 338 !TODO: This will have to be modified in the future when ABINIT enables to calculate specific 339 !components of the d2_dkdk 340 nq1q2grad=9 341 ABI_MALLOC(q1q2grad,(4,nq1q2grad)) 342 do iq1q2grad=1,nq1q2grad 343 call rf2_getidirs(iq1q2grad,iq1dir,iq2dir) 344 if (iq1dir==iq2dir) then 345 q1q2grad(1,iq1q2grad)=dtset%natom+10 !d2_dkdk perturbation diagonal elements 346 else 347 q1q2grad(1,iq1q2grad)=dtset%natom+11 !d2_dkdk perturbation off-diagonal elements 348 end if 349 q1q2grad(2,iq1q2grad)=iq1dir !dq1 direction 350 q1q2grad(3,iq1q2grad)=iq2dir !dq2 direction 351 iq1q2grad_var=iq1q2grad 352 if (iq1q2grad>6) iq1q2grad_var=iq1q2grad-3 !Lower=Upper diagonal triangle matrix 353 q1q2grad(4,iq1q2grad)=iq1q2grad_var+(dtset%natom+6)*3 !like pertcase in dfpt_loopert.f90 354 end do 355 356 !################# ELECTROSTATIC CONTRIBUTIONS ####################################### 357 358 !This is necessary to deactivate paw options in the dfpt_rhotov routine 359 ABI_MALLOC(pawrhoij_read,(0)) 360 pawread=0 361 nhat1grdim=0 362 ABI_MALLOC(nhat1gr,(0,0,0)) 363 ABI_MALLOC(nhat,(nfft,nspden)) 364 nhat=zero 365 ABI_MALLOC(nhat1,(cplex*nfft,nspden)) 366 nhat1=zero 367 368 !Read the electric field density response from a disk file(rhor1_efield), calculates the FFT 369 !(rhog1_tmp) and the first order Hartree and xc potentials(vhxc1_efield). 370 !TODO: In the call to read_rhor there is a security option that compares with the header 371 !hdr. Not activated at this moment. 372 ABI_MALLOC(rhog1_tmp,(2,nfft)) 373 ABI_MALLOC(rhor1_efield,(nq2grad,cplex*nfft,nspden)) 374 ABI_MALLOC(rhor1_tmp,(cplex*nfft,nspden)) 375 ABI_MALLOC(rhor1_real,(1*nfft,nspden)) 376 ABI_MALLOC(vhartr1,(cplex*nfft)) 377 ABI_MALLOC(vhxc1_efield,(nq2grad,cplex*nfft)) 378 ABI_MALLOC(vpsp1,(cplex*nfft)) 379 ABI_MALLOC(vtrial1,(cplex*nfft,nspden)) 380 ABI_MALLOC(vresid1,(cplex*nfft,nspden)) 381 ABI_MALLOC(vxc1,(cplex*nfft,nspden)) 382 ABI_MALLOC(dum_vxc,(nfft,nspden)) 383 ABI_MALLOC(xccc3d1,(cplex*n3xccc)) 384 vpsp1=zero; vtrial1=zero; dum_vxc=zero 385 optene=0; optres=1 386 do iq2grad=1,nq2grad 387 pertcase=q2grad(3,iq2grad) 388 389 !Reads a real first order density 390 call appdig(pertcase,dtfil%fildens1in,fi1o) 391 call read_rhor(fi1o, 1, nspden, nfft, ngfft, pawread, mpi_enreg, rhor1_real, & 392 & hdr_den, pawrhoij_read, spaceworld) 393 call hdr_den%free() 394 395 !Perform FFT rhor1 to rhog1 396 call fourdp(cplex,rhog1_tmp,rhor1_real,-1,mpi_enreg,nfft,1,ngfft,0) 397 398 !Accumulate density in meaningful complex arrays 399 if (timrev==0) then 400 do ii=1,nfft 401 jj=ii*2 402 rhor1_tmp(jj-1,:)=rhor1_real(ii,:) 403 rhor1_tmp(jj,:)=zero 404 end do 405 else if (timrev==1) then 406 rhor1_tmp(:,:)=rhor1_real(:,:) 407 end if 408 rhor1_efield(iq2grad,:,:)=rhor1_tmp(:,:) 409 410 !Calculate first order Hartree and xc potentials 411 call dfpt_rhotov(cplex,dum_scl,dum_scl,dum_scl,dum_scl,dum_scl, & 412 & gsqcut,q2grad(2,iq2grad),dtset%natom+2,& 413 & dtset%ixc,kxc,mpi_enreg,dtset%natom,nfft,ngfft,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,& 414 & nspden,n3xccc,non_magnetic_xc,optene,optres,dtset%qptn,rhog,rhog1_tmp,rhor,rhor1_tmp,& 415 & rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1,vpsp1,vresid1,vres2,vtrial1,dum_vxc,vxc1,& 416 & xccc3d1,dtset%ixcrot) 417 418 419 !Accumulate the potential in meaningful arrays 420 vhxc1_efield(iq2grad,:)=vtrial1(:,nspden) 421 422 end do 423 424 !Read the atomic displacement density response from a disk file, calculate the FFT 425 !(rhog1_atdis) and the first order Hartree and xc potentials(vhxc1_atdis). 426 ABI_MALLOC(rhog1_atdis,(natpert,2,nfft)) 427 ABI_MALLOC(rhor1_atdis,(natpert,cplex*nfft,nspden)) 428 ABI_MALLOC(vhxc1_atdis,(natpert,cplex*nfft)) 429 vtrial1=zero 430 do iatpert= 1, natpert 431 iatpol=pert_atdis(1,iatpert) 432 iatdir=pert_atdis(2,iatpert) 433 pertcase=pert_atdis(3,iatpert) 434 435 !Reads a real first order density 436 call appdig(pertcase,dtfil%fildens1in,fi1o) 437 call read_rhor(fi1o, 1, nspden, nfft, ngfft, pawread, mpi_enreg, rhor1_real, & 438 & hdr_den, pawrhoij_read, spaceworld) 439 call hdr_den%free() 440 441 !Perform FFT rhor1 to rhog1 442 call fourdp(cplex,rhog1_tmp,rhor1_real,-1,mpi_enreg,nfft,1,ngfft,0) 443 444 !Accumulate density in meaningful complex arrays 445 if (timrev==0) then 446 do ii=1,nfft 447 jj=ii*2 448 rhor1_tmp(jj-1,:)=rhor1_real(ii,:) 449 end do 450 else if (timrev==1) then 451 rhor1_tmp(:,:)=rhor1_real(:,:) 452 end if 453 rhor1_atdis(iatpert,:,:)=rhor1_tmp(:,:) 454 rhog1_atdis(iatpert,:,:)=rhog1_tmp(:,:) 455 456 !Calculate first order Hartree and xc potentials 457 call dfpt_rhotov(cplex,dum_scl,dum_scl,dum_scl,dum_scl,dum_scl, & 458 & gsqcut,iatdir,iatpol,& 459 & dtset%ixc,kxc,mpi_enreg,dtset%natom,nfft,ngfft,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,& 460 & nspden,n3xccc,non_magnetic_xc,optene,optres,dtset%qptn,rhog,rhog1_tmp,rhor,rhor1_tmp,& 461 & rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1,vpsp1,vresid1,vres2,vtrial1,dum_vxc,vxc1,& 462 & xccc3d1,dtset%ixcrot) 463 464 !Accumulate the potential in meaningful arrays 465 vhxc1_atdis(iatpert,:)=vtrial1(:,nspden) 466 467 end do 468 469 !These arrays will not be used anymore (for the moment) 470 ABI_FREE(rhor1_real) 471 ABI_FREE(rhor1_tmp) 472 ABI_FREE(vhartr1) 473 ABI_FREE(vpsp1) 474 ABI_FREE(vtrial1) 475 ABI_FREE(vresid1) 476 ABI_FREE(vxc1) 477 ABI_FREE(dum_vxc) 478 ABI_FREE(xccc3d1) 479 480 ABI_FREE(pawrhoij_read) 481 ABI_FREE(nhat1gr) 482 ABI_FREE(nhat) 483 ABI_FREE(nhat1) 484 485 !Calculate the electrostatic contribution from the q-gradient of the Hxc kernel 486 487 !If GGA xc first calculate the Cartesian q-gradient of the xc kernel 488 if (nkxc == 7) then 489 ABI_MALLOC(vxc1dq,(2*nfft,nspden)) 490 ABI_MALLOC(vxc1dqc,(2*nfft,nspden,natpert,3)) 491 ABI_MALLOC(rhor1_tmp,(cplex*nfft,nspden)) 492 do qcar=1,3 493 do iatpert=1,natpert 494 rhor1_tmp(:,:)=rhor1_atdis(iatpert,:,:) 495 call dfpt_mkvxcggadq(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,qcar,rhor1_tmp,vxc1dq) 496 vxc1dqc(:,:,iatpert,qcar)=vxc1dq(:,:) 497 end do 498 end do 499 ABI_FREE(rhor1_tmp) 500 end if 501 502 ABI_MALLOC(rhor1_cplx,(2*nfft,nspden)) 503 ABI_MALLOC(vqgradhart,(2*nfft)) 504 ABI_MALLOC(eqgradhart,(2,natpert,nq2grad,nq1grad)) 505 ABI_MALLOC(qdrflg,(matom,3,3,3)) 506 qdrflg=0 507 rhor1_cplx=zero 508 do iq1grad=1,nq1grad 509 qdir=q1grad(2,iq1grad) 510 do iatpert=1,natpert 511 512 !Calculate the gradient of the Hartree potential generated by the first order atomic displacement density 513 rhog1_tmp(:,:)=rhog1_atdis(iatpert,:,:) 514 call hartredq(2,gmet,gsqcut,mpi_enreg,nfft,ngfft,qdir,rhog1_tmp,vqgradhart) 515 516 !To ckeck 517 ! call appdig(pert_atdis(3,iatpert)+q1grad(3,iq1grad),"Hartree_",fi1o) 518 ! call fftdatar_write_from_hdr("first_order_potential",fi1o,dtset%iomode,hdr_den,& 519 ! & ngfft,2,nfft,nspden,vqgradhart,mpi_enreg) 520 521 !If GGA convert the gradient of xc kernel to reduced coordinates and incorporate it to the Hartree part 522 if (nkxc == 7) then 523 vxc1dq=zero 524 do qcar=1,3 525 vxc1dq(:,:)=vxc1dq(:,:) + gprimd(qcar,qdir) * & 526 & vxc1dqc(:,:,iatpert,qcar) 527 end do 528 vqgradhart(:)=vqgradhart(:)+vxc1dq(:,1) 529 end if 530 531 do iq2grad=1,nq2grad 532 533 !Calculate the electrostatic energy term with the first order electric field density 534 !I need a cplex density for the dotprod_vn 535 do ii=1,nfft 536 jj=ii*2 537 rhor1_cplx(jj-1,:)=rhor1_efield(iq2grad,ii,:) 538 end do 539 540 call dotprod_vn(2,rhor1_cplx,dotr,doti,nfft,nfftot,nspden,2,vqgradhart,ucvol) 541 eqgradhart(re,iatpert,iq2grad,iq1grad)=dotr*half 542 eqgradhart(im,iatpert,iq2grad,iq1grad)=doti*half 543 544 qdrflg(pert_atdis(1,iatpert),pert_atdis(2,iatpert),q2grad(2,iq2grad),q1grad(2,iq1grad))=1 545 546 blkflg(q2grad(2,iq2grad),q2grad(1,iq2grad),pert_atdis(2,iatpert),pert_atdis(1,iatpert),& 547 & q1grad(2,iq1grad),matom+8)=1 548 549 end do 550 end do 551 end do 552 ABI_FREE(rhor1_cplx) 553 ABI_FREE(rhog1_tmp) 554 ABI_FREE(rhog1_atdis) 555 ABI_FREE(rhor1_atdis) 556 ABI_FREE(rhor1_efield) 557 ABI_FREE(vqgradhart) 558 if (nkxc == 7) then 559 ABI_FREE(vxc1dq) 560 ABI_FREE(vxc1dqc) 561 end if 562 563 !################# WAVE FUNCTION CONTRIBUTIONS ####################################### 564 565 !Determine the subset of symmetry operations (nsym1 operations) 566 !that leaves the perturbation invariant, and initialize corresponding arrays 567 !symaf1, symrl1, tnons1 (and pawang1%zarot, if PAW).. 568 !MR TODO: For the moment only the identiy symmetry is activated (nsym1=1) 569 ! In a future I will try to activate perturbation dependent symmetries 570 ! with littlegroup_pert.F90. 571 nsym1 = 1 572 ABI_MALLOC(indsy1,(4,nsym1,dtset%natom)) 573 ABI_MALLOC(symrc1,(3,3,nsym1)) 574 ABI_MALLOC(symaf1,(nsym1)) 575 ABI_MALLOC(symrl1,(3,3,nsym1)) 576 ABI_MALLOC(tnons1,(3,nsym1)) 577 symaf1(1:nsym1)= 1 578 symrl1(:,:,nsym1)= dtset%symrel(:,:,1) 579 tnons1(:,nsym1)= 0_dp 580 581 !Set up corresponding symmetry data 582 ABI_MALLOC(irrzon1,(dtset%nfft**(1-1/nsym1),2,(nspden/dtset%nsppol)-3*(nspden/4))) 583 ABI_MALLOC(phnons1,(2,dtset%nfft**(1-1/nsym1),(nspden/dtset%nsppol)-3*(nspden/4))) 584 call setsym(indsy1,irrzon1,1,dtset%natom,dtset%nfft,dtset%ngfft,nspden,dtset%nsppol,& 585 &nsym1,phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred) 586 587 ABI_FREE(indsy1) 588 ABI_FREE(symaf1) 589 ABI_FREE(symrl1) 590 ABI_FREE(tnons1) 591 592 !Determine the subset of k-points needed in the "reduced Brillouin zone", 593 !and initialize other quantities 594 ABI_MALLOC(indkpt1_tmp,(nkpt)) 595 ABI_MALLOC(wtk_folded,(nkpt)) 596 ABI_MALLOC(bz2ibz_smap,(6, nkpt)) 597 indkpt1_tmp(:)=0 598 599 if (dtset%kptopt==2) then 600 call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,& 601 & nsym1,symrc1,timrev,dtset%wtk,wtk_folded,bz2ibz_smap,xmpi_comm_self) 602 else if (dtset%kptopt==3) then 603 call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,& 604 & nsym1,symrc1,0,dtset%wtk,wtk_folded,bz2ibz_smap,xmpi_comm_self) 605 else 606 write(msg,"(1a)") 'kptopt must be 2 or 3 for the quadrupole calculation' 607 ABI_BUG(msg) 608 end if 609 ABI_FREE(bz2ibz_smap) 610 611 ABI_MALLOC(doccde_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol)) 612 ABI_MALLOC(indkpt1,(nkpt_rbz)) 613 ABI_MALLOC(istwfk_rbz,(nkpt_rbz)) 614 ABI_MALLOC(kpt_rbz,(3,nkpt_rbz)) 615 ABI_MALLOC(nband_rbz,(nkpt_rbz*dtset%nsppol)) 616 ABI_MALLOC(occ_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol)) 617 ABI_MALLOC(wtk_rbz,(nkpt_rbz)) 618 indkpt1(:)=indkpt1_tmp(1:nkpt_rbz) 619 do ikpt=1,nkpt_rbz 620 istwfk_rbz(ikpt)=dtset%istwfk(indkpt1(ikpt)) 621 kpt_rbz(:,ikpt)=dtset%kptns(:,indkpt1(ikpt)) 622 wtk_rbz(ikpt)=wtk_folded(indkpt1(ikpt)) 623 end do 624 ABI_FREE(indkpt1_tmp) 625 ABI_FREE(wtk_folded) 626 627 !Transfer occ to occ_rbz 628 !NOTE : this takes into account that indkpt1 is ordered 629 !MG: What about using occ(band,kpt,spin) ??? 630 bdtot_index=0;bdtot1_index=0 631 do isppol=1,dtset%nsppol 632 ikpt1=1 633 do ikpt=1,nkpt 634 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 635 ! Must test against ikpt1/=nkpt_rbz+1, before evaluate indkpt1(ikpt1) 636 if(ikpt1/=nkpt_rbz+1)then 637 if(ikpt==indkpt1(ikpt1))then 638 nband_rbz(ikpt1+(isppol-1)*nkpt_rbz)=nband_k 639 occ_rbz(1+bdtot1_index:nband_k+bdtot1_index)=occ(1+bdtot_index:nband_k+bdtot_index) 640 doccde_rbz(1+bdtot1_index:nband_k+bdtot1_index)=doccde(1+bdtot_index:nband_k+bdtot_index) 641 ikpt1=ikpt1+1 642 bdtot1_index=bdtot1_index+nband_k 643 end if 644 end if 645 bdtot_index=bdtot_index+nband_k 646 end do 647 end do 648 649 !Compute maximum number of planewaves at k 650 call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpt_rbz,mpi_enreg,mpw,nkpt_rbz) 651 652 !Allocate some k-dependent arrays at k 653 ABI_MALLOC(kg,(3,mpw*nkpt_rbz)) 654 ABI_MALLOC(kg_k,(3,mpw)) 655 ABI_MALLOC(npwarr,(nkpt_rbz)) 656 ABI_MALLOC(npwtot,(nkpt_rbz)) 657 658 !Determine distribution of k-points/bands over MPI processes 659 if (allocated(mpi_enreg%my_kpttab)) then 660 ABI_FREE(mpi_enreg%my_kpttab) 661 end if 662 ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt_rbz)) 663 if(xmpi_paral==1) then 664 ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt_rbz,dtset%mband,dtset%nsppol)) 665 call distrb2(dtset%mband,mband_mem,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,dtset%nsppol,mpi_enreg) 666 else 667 mpi_enreg%my_kpttab(:)=(/(ii,ii=1,nkpt_rbz)/) 668 end if 669 my_nkpt_rbz=maxval(mpi_enreg%my_kpttab) 670 mkmem_rbz =my_nkpt_rbz 671 call initmpi_band(mkmem_rbz,mpi_enreg,nband_rbz,nkpt_rbz,dtset%nsppol) 672 673 !Set up the basis sphere of planewaves at k 674 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg,& 675 & kpt_rbz,mkmem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw,npwarr,npwtot,dtset%nsppol) 676 ABI_FREE(npwtot) 677 678 !Set up the spherical harmonics (Ylm) and 1st gradients at k 679 useylmgr=1; option=1 ; nylmgr=3 680 ABI_MALLOC(ylm,(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm)) 681 ABI_MALLOC(ylmgr,(mpw*mkmem_rbz,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 682 if (psps%useylm==1) then 683 call initylmg(gprimd,kg,kpt_rbz,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,nband_rbz,nkpt_rbz,& 684 & npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr) 685 end if 686 687 !Initialize band structure datatype at k 688 bantot_rbz=sum(nband_rbz(1:nkpt_rbz*dtset%nsppol)) 689 ABI_MALLOC(eigen0,(bantot_rbz)) 690 eigen0(:)=zero 691 ! CP modified 692 ! call ebands_init(bantot_rbz,bs_rbz,dtset%nelect,doccde_rbz,eigen0,istwfk_rbz,kpt_rbz,& 693 !& nband_rbz,nkpt_rbz,npwarr,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,& 694 !& dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 695 !& dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 696 call ebands_init(bantot_rbz,bs_rbz,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,& 697 & doccde_rbz,eigen0,istwfk_rbz,kpt_rbz,& 698 & nband_rbz,nkpt_rbz,npwarr,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,& 699 & dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 700 & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 701 ! End CP modified 702 ABI_FREE(eigen0) 703 704 ABI_FREE(doccde_rbz) 705 706 !Initialize header, update it with evolving variables 707 gscase=0 ! A GS WF file is read 708 call hdr_init(bs_rbz,codvsn,dtset,hdr0,pawtab,gscase,psps,wvl%descr,& 709 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 710 711 ! CP modified 712 ! call hdr0%update(bantot_rbz,etotal,fermie,& 713 !& residm,rprimd,occ_rbz,pawrhoij,xred,dtset%amu_orig(:,1),& 714 !& comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 715 call hdr0%update(bantot_rbz,etotal,fermie,fermie,& 716 & residm,rprimd,occ_rbz,pawrhoij,xred,dtset%amu_orig(:,1),& 717 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 718 ! End CP modified 719 720 !Clean band structure datatype (should use it more in the future !) 721 call ebands_free(bs_rbz) 722 723 !!Initialize GS wavefunctions at k 724 ireadwf0=1; formeig=0 ; ask_accurate=1 ; optorth=0 725 mcg=mpw*dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol 726 if (one*mpw*dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol > huge(1)) then 727 write (msg,'(4a, 5(a,i0), 2a)')& 728 & "Default integer is not wide enough to store the size of the GS wavefunction array (WF0, mcg).",ch10,& 729 & "Action: increase the number of processors. Consider also OpenMP threads.",ch10,& 730 & "nspinor: ",dtset%nspinor, "mpw: ",mpw, "mband: ",dtset%mband, "mkmem_rbz: ",& 731 & mkmem_rbz, "nsppol: ",dtset%nsppol,ch10,& 732 & 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode' 733 ABI_ERROR(msg) 734 end if 735 ABI_STAT_MALLOC(cg,(2,mcg), ierr) 736 ABI_CHECK(ierr==0, "out-of-memory in cg") 737 738 ABI_MALLOC(eigen0,(dtset%mband*nkpt_rbz*dtset%nsppol)) 739 call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,& 740 & formeig,hdr0,ireadwf0,istwfk_rbz,kg,& 741 & kpt_rbz,dtset%localrdwf,dtset%mband,mcg,& 742 & mkmem_rbz,mpi_enreg,mpw,nband_rbz,dtset%ngfft,nkpt_rbz,npwarr,& 743 & dtset%nsppol,dtset%nsym,occ_rbz,optorth,dtset%symafm,& 744 & dtset%symrel,dtset%tnons,dtfil%unkg,wffgs,wfftgs,& 745 & dtfil%unwffgs,dtfil%fnamewffk,wvl) 746 ABI_FREE(eigen0) 747 748 !Close wffgs%unwff, if it was ever opened (in inwffil) 749 if (ireadwf0==1) then 750 call WffClose(wffgs,ierr) 751 end if 752 753 !==== Initialize most of the Hamiltonian ==== 754 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 755 !2) Perform the setup needed for the non-local factors: 756 !3) Constant kleimann-Bylander energies are copied from psps to gs_hamkq. 757 call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,nsppol,nspden,dtset%natom,& 758 & dtset%typat,xred,nfft,dtset%mgfft,ngfft,rprimd,dtset%nloalg,ph1d=ph1d,& 759 & use_gpu_cuda=dtset%use_gpu_cuda) 760 761 762 !==== Initialize response functions files and handlers ==== 763 !Atomic displacement files 764 ABI_MALLOC(wfk_t_atdis,(natpert)) 765 766 do iatpert=1,natpert 767 768 pertcase=pert_atdis(3,iatpert) 769 call appdig(pertcase,dtfil%fnamewff1,fiwfatdis) 770 771 !The value 20 is taken arbitrarily I would say 772 forunit=20+pertcase 773 774 !Check that atdis file exists and open it 775 if (.not. file_exists(fiwfatdis)) then 776 ! Trick needed to run Abinit test suite in netcdf mode. 777 if (file_exists(nctk_ncify(fiwfatdis))) then 778 write(msg,"(3a)")"- File: ",trim(fiwfatdis),& 779 " does not exist but found netcdf file with similar name." 780 call wrtout(std_out,msg,'COLL') 781 fiwfatdis = nctk_ncify(fiwfatdis) 782 end if 783 if (.not. file_exists(fiwfatdis)) then 784 ABI_ERROR('Missing file: '//TRIM(fiwfatdis)) 785 end if 786 end if 787 write(msg,'(a,a)')'-open atomic displacement wf1 file :',trim(fiwfatdis) 788 call wrtout([std_out, ab_out], msg) 789 call wfk_open_read(wfk_t_atdis(iatpert),fiwfatdis,formeig1,dtset%iomode,forunit, xmpi_comm_self) 790 end do 791 792 !ddk files 793 ABI_MALLOC(wfk_t_ddk,(nq1grad)) 794 do iq1grad=1,nq1grad 795 796 pertcase=q1grad(3,iq1grad) 797 call appdig(pertcase,dtfil%fnamewffddk,fiwfddk) 798 799 !The value 20 is taken arbitrarily I would say 800 forunit=20+pertcase 801 802 !Check that ddk file exists and open it 803 if (.not. file_exists(fiwfddk)) then 804 ! Trick needed to run Abinit test suite in netcdf mode. 805 if (file_exists(nctk_ncify(fiwfddk))) then 806 write(msg,"(3a)")"- File: ",trim(fiwfddk),& 807 " does not exist but found netcdf file with similar name." 808 call wrtout(std_out,msg,'COLL') 809 fiwfddk = nctk_ncify(fiwfddk) 810 end if 811 if (.not. file_exists(fiwfddk)) then 812 ABI_ERROR('Missing file: '//TRIM(fiwfddk)) 813 end if 814 end if 815 write(msg, '(a,a)') '-open ddk wf1 file :',trim(fiwfddk) 816 call wrtout(std_out,msg,'COLL') 817 call wrtout(ab_out,msg,'COLL') 818 call wfk_open_read(wfk_t_ddk(iq1grad),fiwfddk,formeig1,dtset%iomode,forunit,spaceworld) 819 820 end do 821 822 !Electric field files 823 ABI_MALLOC(wfk_t_efield,(nq2grad)) 824 do iq2grad=1,nq2grad 825 826 pertcase=q2grad(3,iq2grad) 827 call appdig(pertcase,dtfil%fnamewff1,fiwfefield) 828 829 !The value 20 is taken arbitrarily I would say 830 forunit=20+pertcase 831 832 !Check that efield file exists and open it 833 if (.not. file_exists(fiwfefield)) then 834 ! Trick needed to run Abinit test suite in netcdf mode. 835 if (file_exists(nctk_ncify(fiwfefield))) then 836 write(msg,"(3a)")"- File: ",trim(fiwfefield),& 837 " does not exist but found netcdf file with similar name." 838 call wrtout(std_out,msg,'COLL') 839 fiwfefield = nctk_ncify(fiwfefield) 840 end if 841 if (.not. file_exists(fiwfefield)) then 842 ABI_ERROR('Missing file: '//TRIM(fiwfefield)) 843 end if 844 end if 845 write(msg, '(a,a)') '-open electric field wf1 file :',trim(fiwfefield) 846 call wrtout(std_out,msg,'COLL') 847 call wrtout(ab_out,msg,'COLL') 848 call wfk_open_read(wfk_t_efield(iq2grad),fiwfefield,formeig1,dtset%iomode,forunit,spaceworld) 849 850 end do 851 852 !d2_dkdk 853 ABI_MALLOC(wfk_t_dkdk,(nq1q2grad)) 854 do iq1q2grad=1,nq1q2grad 855 856 pertcase=q1q2grad(4,iq1q2grad) 857 call appdig(pertcase,dtfil%fnamewffdkdk,fiwfdkdk) 858 859 !The value 20 is taken arbitrarily I would say 860 forunit=20+pertcase 861 862 !Check that d2_ddk file exists and open it 863 if (.not. file_exists(fiwfdkdk)) then 864 ! Trick needed to run Abinit test suite in netcdf mode. 865 if (file_exists(nctk_ncify(fiwfdkdk))) then 866 write(msg,"(3a)")"- File: ",trim(fiwfdkdk),& 867 " does not exist but found netcdf file with similar name." 868 call wrtout(std_out,msg,'COLL') 869 fiwfdkdk = nctk_ncify(fiwfdkdk) 870 end if 871 if (.not. file_exists(fiwfdkdk)) then 872 ABI_ERROR('Missing file: '//TRIM(fiwfdkdk)) 873 end if 874 end if 875 if (iq1q2grad <= 6) then 876 write(msg, '(a,a)') '-open d2_dkdk wf2 file :',trim(fiwfdkdk) 877 call wrtout(std_out,msg,'COLL') 878 call wrtout(ab_out,msg,'COLL') 879 call wfk_open_read(wfk_t_dkdk(iq1q2grad),fiwfdkdk,formeig1,dtset%iomode,forunit,spaceworld) 880 else 881 wfk_t_dkdk(iq1q2grad)=wfk_t_dkdk(iq1q2grad-3) 882 end if 883 884 end do 885 886 !Allocate the quadrupole tensor part depending on the wave functions 887 ABI_MALLOC(qdrpwf,(2,natpert,nq2grad,nq1grad)) 888 ABI_MALLOC(qdrpwf_k,(2,natpert,nq2grad,nq1grad)) 889 ABI_MALLOC(qdrpwf_t1,(2,natpert,nq2grad,nq1grad)) 890 ABI_MALLOC(qdrpwf_t1_k,(2,natpert,nq2grad,nq1grad)) 891 ABI_MALLOC(qdrpwf_t2,(2,natpert,nq2grad,nq1grad)) 892 ABI_MALLOC(qdrpwf_t2_k,(2,natpert,nq2grad,nq1grad)) 893 ABI_MALLOC(qdrpwf_t3,(2,natpert,nq2grad,nq1grad)) 894 ABI_MALLOC(qdrpwf_t3_k,(2,natpert,nq2grad,nq1grad)) 895 ABI_MALLOC(qdrpwf_t4,(2,natpert,nq2grad,nq1grad)) 896 ABI_MALLOC(qdrpwf_t4_k,(2,natpert,nq2grad,nq1grad)) 897 ABI_MALLOC(qdrpwf_t5,(2,natpert,nq2grad,nq1grad)) 898 ABI_MALLOC(qdrpwf_t5_k,(2,natpert,nq2grad,nq1grad)) 899 qdrpwf=zero 900 qdrpwf_t1=zero 901 qdrpwf_t2=zero 902 qdrpwf_t3=zero 903 qdrpwf_t4=zero 904 qdrpwf_t5=zero 905 906 !LOOP OVER SPINS 907 bdtot_index=0 908 icg=0 909 do isppol=1,nsppol 910 ikg=0 911 912 ! Continue to initialize the Hamiltonian 913 call gs_hamkq%load_spin(isppol,with_nonlocal=.true.) 914 915 ! BIG FAT k POINT LOOP 916 do ikpt=1,nkpt_rbz 917 918 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 919 istwf_k=istwfk_rbz(ikpt) 920 npw_k=npwarr(ikpt) 921 npw1_k=npw_k 922 923 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 924 bdtot_index=bdtot_index+nband_k 925 926 cycle ! Skip the rest of the k-point loop 927 end if 928 929 ABI_MALLOC(occ_k,(nband_k)) 930 ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 931 ABI_MALLOC(ylmgr_k,(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 932 occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index) 933 kpoint(:)=kpt_rbz(:,ikpt) 934 wtk_k=wtk_rbz(ikpt) 935 936 ! Get plane-wave vectors and related data at k 937 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 938 if (psps%useylm==1) then 939 do ilm=1,psps%mpsang*psps%mpsang 940 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 941 end do 942 if (useylmgr==1) then 943 do ilm=1,psps%mpsang*psps%mpsang 944 do ii=1,nylmgr 945 ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm) 946 end do 947 end do 948 end if 949 end if 950 951 call dfpt_qdrpwf(atindx,cg,cplex,dtset,gs_hamkq,gsqcut, & 952 & icg,ikpt,indkpt1,isppol,istwf_k, & 953 & kg_k,kpoint,mkmem_rbz, & 954 & mpi_enreg,mpw,natpert,nattyp,nband_k,nfft,ngfft,nkpt_rbz, & 955 & npw_k,nq1grad, & 956 & nq2grad,nq1q2grad,nspden,nsppol,nylmgr,occ_k, & 957 & pert_atdis,ph1d,psps,qdrpwf_k,qdrpwf_t1_k,qdrpwf_t2_k,qdrpwf_t3_k, & 958 & qdrpwf_t4_k,qdrpwf_t5_k,q1grad,q2grad,q1q2grad,rmet,ucvol,useylmgr, & 959 & vhxc1_atdis,vhxc1_efield,wfk_t_atdis,wfk_t_efield,wfk_t_ddk, & 960 & wfk_t_dkdk,wtk_k,xred,ylm_k,ylmgr_k) 961 962 963 ! Add the contribution from each k-point 964 qdrpwf=qdrpwf + qdrpwf_k 965 qdrpwf_t1=qdrpwf_t1 + qdrpwf_t1_k 966 qdrpwf_t2=qdrpwf_t2 + qdrpwf_t2_k 967 qdrpwf_t3=qdrpwf_t3 + qdrpwf_t3_k 968 qdrpwf_t4=qdrpwf_t4 + qdrpwf_t4_k 969 qdrpwf_t5=qdrpwf_t5 + qdrpwf_t5_k 970 971 ! Keep track of total number of bands 972 bdtot_index=bdtot_index+nband_k 973 974 ! Shift arrays memory 975 if (mkmem_rbz/=0) then 976 icg=icg+npw_k*dtset%nspinor*nband_k 977 ikg=ikg+npw_k 978 end if 979 980 ABI_FREE(occ_k) 981 ABI_FREE(ylm_k) 982 ABI_FREE(ylmgr_k) 983 984 end do 985 ! END BIG FAT k POINT LOOP 986 end do 987 !END LOOP OVER SPINS 988 989 !Memory cleaning 990 call gs_hamkq%free() 991 992 !Close response function files 993 do iatpert=1,natpert 994 call wfk_t_atdis(iatpert)%close() 995 end do 996 do iq1grad=1,nq1grad 997 call wfk_t_ddk(iq1grad)%close() 998 end do 999 do iq2grad=1,nq2grad 1000 call wfk_t_efield(iq2grad)%close() 1001 end do 1002 do iq1q2grad=1,nq1q2grad 1003 if (iq1q2grad <= 6) call wfk_t_dkdk(iq1q2grad)%close() 1004 end do 1005 1006 !=== MPI communications ================== 1007 if (xmpi_paral==1) then 1008 1009 call xmpi_sum(qdrpwf,spaceworld,ierr) 1010 call xmpi_sum(qdrpwf_t1,spaceworld,ierr) 1011 call xmpi_sum(qdrpwf_t2,spaceworld,ierr) 1012 call xmpi_sum(qdrpwf_t3,spaceworld,ierr) 1013 call xmpi_sum(qdrpwf_t4,spaceworld,ierr) 1014 call xmpi_sum(qdrpwf_t5,spaceworld,ierr) 1015 1016 end if 1017 1018 !Anounce finalization of quadrupole tensor calculation 1019 write(msg, '(a,a,a)' ) ch10, & 1020 ' Quadrupole tensor calculation completed ',ch10 1021 call wrtout(std_out,msg,'COLL') 1022 call wrtout(ab_out,msg,'COLL') 1023 1024 !Gather the different terms in the quadrupole tensor and print them out 1025 if (me==0) then 1026 call dfpt_qdrpout(d3etot,eqgradhart,gprimd,dtset%kptopt,matom,mpert,natpert,& 1027 & nq1grad,nq2grad,pert_atdis,dtset%prtvol,q1grad,q2grad,qdrflg,qdrpwf,qdrpwf_t1,qdrpwf_t2, & 1028 & qdrpwf_t3,qdrpwf_t4,qdrpwf_t5,rprimd,ucvol) 1029 end if 1030 1031 !Deallocations 1032 ABI_FREE(cg) 1033 ABI_FREE(eqgradhart) 1034 ABI_FREE(indkpt1) 1035 ABI_FREE(istwfk_rbz) 1036 ABI_FREE(kg) 1037 ABI_FREE(kg_k) 1038 ABI_FREE(kpt_rbz) 1039 ! ABI_FREE(mpi_enreg%my_kpttab) 1040 ! ABI_FREE(mpi_enreg%proc_distrb) 1041 ABI_FREE(nband_rbz) 1042 ABI_FREE(npwarr) 1043 ABI_FREE(occ_rbz) 1044 ABI_FREE(ph1d) 1045 ABI_FREE(pert_atdis) 1046 ABI_FREE(irrzon1) 1047 ABI_FREE(phnons1) 1048 ABI_FREE(qdrflg) 1049 ABI_FREE(qdrpwf) 1050 ABI_FREE(qdrpwf_k) 1051 ABI_FREE(qdrpwf_t1) 1052 ABI_FREE(qdrpwf_t2) 1053 ABI_FREE(qdrpwf_t3) 1054 ABI_FREE(qdrpwf_t4) 1055 ABI_FREE(qdrpwf_t5) 1056 ABI_FREE(qdrpwf_t1_k) 1057 ABI_FREE(qdrpwf_t2_k) 1058 ABI_FREE(qdrpwf_t3_k) 1059 ABI_FREE(qdrpwf_t4_k) 1060 ABI_FREE(qdrpwf_t5_k) 1061 ABI_FREE(q1grad) 1062 ABI_FREE(q1q2grad) 1063 ABI_FREE(q2grad) 1064 ABI_FREE(symrc1) 1065 ABI_FREE(vhxc1_atdis) 1066 ABI_FREE(vhxc1_efield) 1067 ABI_FREE(wfk_t_atdis) 1068 ABI_FREE(wfk_t_ddk) 1069 ABI_FREE(wfk_t_dkdk) 1070 ABI_FREE(wfk_t_efield) 1071 ABI_FREE(wtk_rbz) 1072 ABI_FREE(ylm) 1073 ABI_FREE(ylmgr) 1074 if(xmpi_paral==1) then 1075 ABI_FREE(mpi_enreg%proc_distrb) 1076 end if 1077 1078 ! Clean the header 1079 call hdr0%free() 1080 1081 DBG_EXIT("COLL") 1082 1083 end subroutine dfpt_qdrpole
ABINIT/dfpt_qdrpout [ Functions ]
NAME
dfpt_qdrpout
FUNCTION
This subroutine gathers the different terms entering the quadrupole tensor, perfofms the transformation from reduced to cartesian coordinates and writes out the quadrupole tensor in external files. It also writes the firts q-gradient of the polarization response to an atomic displacement.
INPUTS
eqgradhart(2,natpert,nq2grad,nq1grad)=electrostatic contribution from the q-gradient of the Hartree potential gprimd(3,3)=reciprocal space dimensional primitive translations kptopt=2 time reversal symmetry is enforced, 3 trs is not enforced (for debugging purposes) matom=maximum number of atoms to displace natpert=number of atomic displacement perturbations nq1grad=number of q1 (q_{\gamma}) gradients nq2grad=number of q2 (q_{\delta}) gradients pert_atdis(3,natpert)=array with the info for the atomic displacement perturbations prtvol=volume of information to be printed. 1-> The different contributions to the quadrupole are printed. q1grad(3,nq1grad)=array with the info for the q1 (q_{\gamma}) gradients q2grad(3,nq2grad)=array with the info for the q2 (q_{\gamma}) gradients qdrflg(matom,3,3,3)=array that indicates which quadrupole tensor elements have been calculated qdrpwf_k(2,natpert,nq2grad,nq1grad)= wave function dependent part of the quadrupole tensor for the k-point kpt qdrpwf_t1_k(2,natpert,nq2grad,nq1grad)= t1 term (see notes) of qdrpwf_k qdrpwf_t2_k(2,natpert,nq2grad,nq1grad)= t2 term (see notes) of qdrpwf_k qdrpwf_t3_k(2,natpert,nq2grad,nq1grad)= t3 term (see notes) of qdrpwf_k qdrpwf_t4_k(2,natpert,nq2grad,nq1grad)= t4 term (see notes) of qdrpwf_k qdrpwf_t5_k(2,natpert,nq2grad,nq1grad)= t5 term (see notes) of qdrpwf_k rprimd(3,3)=dimensional primitive translations in real space (bohr) ucvol=unit cell volume in bohr**3.
OUTPUT
d3etot(2,3,mpert,6,mpert,3,mpert)= matrix of the 3DTE
SIDE EFFECTS
NOTES
SOURCE
1129 subroutine dfpt_qdrpout(d3etot,eqgradhart,gprimd,kptopt,matom,mpert,natpert, & 1130 & nq1grad,nq2grad,pert_atdis,prtvol,q1grad,q2grad,qdrflg,qdrpwf,qdrpwf_t1,qdrpwf_t2, & 1131 & qdrpwf_t3,qdrpwf_t4,qdrpwf_t5,rprimd,ucvol) 1132 1133 !Arguments ------------------------------------ 1134 !scalars 1135 integer,intent(in) :: kptopt,matom,mpert,natpert,nq1grad,nq2grad,prtvol 1136 real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert) 1137 real(dp),intent(in) :: ucvol 1138 1139 !arrays 1140 integer,intent(in) :: pert_atdis(3,natpert) 1141 integer,intent(in) :: qdrflg(matom,3,3,3) 1142 integer,intent(in) :: q1grad(3,nq1grad),q2grad(3,nq2grad) 1143 real(dp),intent(in) :: eqgradhart(2,natpert,nq2grad,nq1grad) 1144 real(dp),intent(in) :: gprimd(3,3) 1145 real(dp),intent(in) :: qdrpwf(2,natpert,nq2grad,nq1grad) 1146 real(dp),intent(in) :: qdrpwf_t1(2,natpert,nq2grad,nq1grad) 1147 real(dp),intent(in) :: qdrpwf_t2(2,natpert,nq2grad,nq1grad) 1148 real(dp),intent(in) :: qdrpwf_t3(2,natpert,nq2grad,nq1grad) 1149 real(dp),intent(in) :: qdrpwf_t4(2,natpert,nq2grad,nq1grad) 1150 real(dp),intent(in) :: qdrpwf_t5(2,natpert,nq2grad,nq1grad) 1151 real(dp),intent(in) :: rprimd(3,3) 1152 1153 !Local variables------------------------------- 1154 !scalars 1155 integer :: alpha,beta,gamma 1156 integer :: iatpert,iatdir,iatom,ii,iiq1grad,iiq2grad 1157 integer :: iq1dir,iq1grad,iq2dir,iq2grad 1158 integer :: iq1pert,iq2pert 1159 integer, parameter :: re=1,im=2 1160 real(dp) :: piezore,piezoim,tmpre, tmpim 1161 character(len=500) :: msg 1162 !arrays 1163 integer :: flg1(3),flg2(3) 1164 integer,allocatable :: cartflg(:,:,:,:) 1165 real(dp) :: vec1(3),vec2(3) 1166 real(dp),allocatable :: qdrptens_cart(:,:,:,:,:),qdrptens_red(:,:,:,:,:) 1167 real(dp),allocatable :: dqpol_cart(:,:,:,:,:),dqpol_red(:,:,:,:,:) 1168 1169 ! ************************************************************************* 1170 1171 DBG_ENTER("COLL") 1172 1173 !Open output files 1174 if (prtvol>=10) then 1175 open(unit=71,file='qdrpl_wf_t1.out',status='unknown',form='formatted',action='write') 1176 open(unit=72,file='qdrpl_wf_t2.out',status='unknown',form='formatted',action='write') 1177 open(unit=73,file='qdrpl_wf_t3.out',status='unknown',form='formatted',action='write') 1178 open(unit=74,file='qdrpl_wf_t4.out',status='unknown',form='formatted',action='write') 1179 open(unit=75,file='qdrpl_wf_t5.out',status='unknown',form='formatted',action='write') 1180 open(unit=76,file='qdrpl_elecstic.out',status='unknown',form='formatted',action='write') 1181 end if 1182 1183 !Gather the different terms in the tensors and print the result 1184 ABI_MALLOC(qdrptens_red,(2,matom,3,3,3)) 1185 ABI_MALLOC(dqpol_red,(2,matom,3,3,3)) 1186 1187 if (kptopt==3) then 1188 1189 !Write real and 'true' imaginary parts of quadrupole tensor and the 1190 !q-gradient of the polarization response 1191 iq1pert=matom+8 1192 iq2pert=matom+2 1193 do iq1grad=1,nq1grad 1194 iq1dir=q1grad(2,iq1grad) 1195 do iq2grad=1,nq2grad 1196 iq2dir=q2grad(2,iq2grad) 1197 do iatpert=1,natpert 1198 iatom=pert_atdis(1,iatpert) 1199 iatdir=pert_atdis(2,iatpert) 1200 1201 if (qdrflg(iatom,iatdir,iq2dir,iq1dir)==1) then 1202 1203 !Calculate and save the third order energy derivative 1204 tmpre=eqgradhart(re,iatpert,iq2grad,iq1grad)+qdrpwf(re,iatpert,iq2grad,iq1grad) 1205 tmpim=eqgradhart(im,iatpert,iq2grad,iq1grad)+qdrpwf(im,iatpert,iq2grad,iq1grad) 1206 d3etot(1,iq2dir,iq2pert,iatdir,iatom,iq1dir,iq1pert)=tmpre 1207 d3etot(2,iq2dir,iq2pert,iatdir,iatom,iq1dir,iq1pert)=tmpim 1208 1209 !Calculate and write the q-gradient of the polarization response 1210 !(without the inverse volume factor) 1211 dqpol_red(1,iatom,iatdir,iq2dir,iq1dir)=-two*tmpim 1212 dqpol_red(2,iatom,iatdir,iq2dir,iq1dir)=two*tmpre 1213 1214 if (qdrflg(iatom,iatdir,iq2dir,iq1dir)==1 .and. qdrflg(iatom,iatdir,iq1dir,iq2dir)==1 ) then 1215 1216 !Avoid double counting when summing up unsymmetrized contributions 1217 do iiq1grad=1,3 1218 if (q2grad(2,iq2grad)==q1grad(2,iiq1grad)) exit 1219 end do 1220 1221 do iiq2grad=1,3 1222 if (q1grad(2,iq1grad)==q2grad(2,iiq2grad)) exit 1223 end do 1224 1225 qdrptens_red(re,iatom,iatdir,iq2dir,iq1dir)=-two* & 1226 & ( eqgradhart(re,iatpert,iq2grad,iq1grad)+eqgradhart(re,iatpert,iiq2grad,iiq1grad) & 1227 & + qdrpwf(re,iatpert,iq2grad,iq1grad)+qdrpwf(re,iatpert,iiq2grad,iiq1grad) ) 1228 qdrptens_red(im,iatom,iatdir,iq2dir,iq1dir)=-two* & 1229 & ( eqgradhart(im,iatpert,iq2grad,iq1grad)+eqgradhart(im,iatpert,iiq2grad,iiq1grad) & 1230 & + qdrpwf(im,iatpert,iq2grad,iq1grad)+qdrpwf(im,iatpert,iiq2grad,iiq1grad) ) 1231 1232 !Divide by the imaginary unit to get the proper observable 1233 !(See M.Royo and M.Stengel paper) 1234 tmpre=qdrptens_red(re,iatom,iatdir,iq2dir,iq1dir) 1235 tmpim=qdrptens_red(im,iatom,iatdir,iq2dir,iq1dir) 1236 qdrptens_red(re,iatom,iatdir,iq2dir,iq1dir)=tmpim 1237 qdrptens_red(im,iatom,iatdir,iq2dir,iq1dir)=-tmpre 1238 1239 if (prtvol>=10) then 1240 !Write individual contributions 1241 write(71,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1242 & qdrpwf_t1(im,iatpert,iq2grad,iq1grad)+qdrpwf_t1(im,iatpert,iiq2grad,iiq1grad), & 1243 & -(qdrpwf_t1(re,iatpert,iq2grad,iq1grad)+qdrpwf_t1(re,iatpert,iiq2grad,iiq1grad)) 1244 1245 write(72,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1246 & qdrpwf_t2(im,iatpert,iq2grad,iq1grad)+qdrpwf_t2(im,iatpert,iiq2grad,iiq1grad), & 1247 & -(qdrpwf_t2(re,iatpert,iq2grad,iq1grad)+qdrpwf_t2(re,iatpert,iiq2grad,iiq1grad)) 1248 1249 write(73,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1250 & qdrpwf_t3(im,iatpert,iq2grad,iq1grad)+qdrpwf_t3(im,iatpert,iiq2grad,iiq1grad), & 1251 & -(qdrpwf_t3(re,iatpert,iq2grad,iq1grad)+qdrpwf_t3(re,iatpert,iiq2grad,iiq1grad)) 1252 1253 write(74,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1254 & qdrpwf_t4(im,iatpert,iq2grad,iq1grad)+qdrpwf_t4(im,iatpert,iiq2grad,iiq1grad), & 1255 & -(qdrpwf_t4(re,iatpert,iq2grad,iq1grad)+qdrpwf_t4(re,iatpert,iiq2grad,iiq1grad)) 1256 1257 write(75,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1258 & qdrpwf_t5(im,iatpert,iq2grad,iq1grad)+qdrpwf_t5(im,iatpert,iiq2grad,iiq1grad), & 1259 & -(qdrpwf_t5(re,iatpert,iq2grad,iq1grad)+qdrpwf_t5(re,iatpert,iiq2grad,iiq1grad)) 1260 1261 write(76,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1262 & eqgradhart(im,iatpert,iq2grad,iq1grad)+eqgradhart(im,iatpert,iiq2grad,iiq1grad), & 1263 & -(eqgradhart(re,iatpert,iq2grad,iq1grad)+eqgradhart(re,iatpert,iiq2grad,iiq1grad)) 1264 end if 1265 1266 end if 1267 end if 1268 end do 1269 end do 1270 end do 1271 1272 else if (kptopt==2) then 1273 1274 !Write real and zero imaginary parts of quadrupole tensor and the 1275 !q-gradient of the polarization response 1276 iq1pert=matom+8 1277 iq2pert=matom+2 1278 do iq1grad=1,nq1grad 1279 iq1dir=q1grad(2,iq1grad) 1280 do iq2grad=1,nq2grad 1281 iq2dir=q2grad(2,iq2grad) 1282 do iatpert=1,natpert 1283 iatom=pert_atdis(1,iatpert) 1284 iatdir=pert_atdis(2,iatpert) 1285 1286 if (qdrflg(iatom,iatdir,iq2dir,iq1dir)==1) then 1287 1288 !Calculate and save the third order energy derivative 1289 tmpim=eqgradhart(im,iatpert,iq2grad,iq1grad)+qdrpwf(im,iatpert,iq2grad,iq1grad) 1290 d3etot(1,iq2dir,iq2pert,iatdir,iatom,iq1dir,iq1pert)=zero 1291 d3etot(2,iq2dir,iq2pert,iatdir,iatom,iq1dir,iq1pert)=tmpim 1292 1293 !Calculate and write the q-gradient of the polarization response 1294 !(without the inverse volume factor) 1295 dqpol_red(1,iatom,iatdir,iq2dir,iq1dir)=-two*tmpim 1296 dqpol_red(2,iatom,iatdir,iq2dir,iq1dir)=0.0_dp 1297 1298 if (qdrflg(iatom,iatdir,iq2dir,iq1dir)==1 .and. qdrflg(iatom,iatdir,iq1dir,iq2dir)==1 ) then 1299 1300 do iiq1grad=1,3 1301 if (q2grad(2,iq2grad)==q1grad(2,iiq1grad)) exit 1302 end do 1303 1304 do iiq2grad=1,3 1305 if (q1grad(2,iq1grad)==q2grad(2,iiq2grad)) exit 1306 end do 1307 1308 qdrptens_red(re,iatom,iatdir,iq2dir,iq1dir)=-two* & 1309 & ( eqgradhart(re,iatpert,iq2grad,iq1grad)+eqgradhart(re,iatpert,iiq2grad,iiq1grad) & 1310 & + qdrpwf(re,iatpert,iq2grad,iq1grad)+qdrpwf(re,iatpert,iiq2grad,iiq1grad) ) 1311 qdrptens_red(im,iatom,iatdir,iq2dir,iq1dir)=-two* & 1312 & ( eqgradhart(im,iatpert,iq2grad,iq1grad)+eqgradhart(im,iatpert,iiq2grad,iiq1grad) & 1313 & + qdrpwf(im,iatpert,iq2grad,iq1grad)+qdrpwf(im,iatpert,iiq2grad,iiq1grad) ) 1314 1315 !Divide by the imaginary unit to get the proper observable 1316 !(See M.Royo and M.Stengel paper) 1317 tmpim=qdrptens_red(im,iatom,iatdir,iq2dir,iq1dir) 1318 qdrptens_red(re,iatom,iatdir,iq2dir,iq1dir)=tmpim 1319 qdrptens_red(im,iatom,iatdir,iq2dir,iq1dir)=0.0_dp 1320 1321 if (prtvol>=10) then 1322 !Write individual contributions 1323 write(71,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1324 & qdrpwf_t1(im,iatpert,iq2grad,iq1grad)+qdrpwf_t1(im,iatpert,iiq2grad,iiq1grad), & 1325 & 0.0_dp 1326 1327 write(72,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1328 & qdrpwf_t2(im,iatpert,iq2grad,iq1grad)+qdrpwf_t2(im,iatpert,iiq2grad,iiq1grad), & 1329 & 0.0_dp 1330 1331 write(73,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1332 & qdrpwf_t3(im,iatpert,iq2grad,iq1grad)+qdrpwf_t3(im,iatpert,iiq2grad,iiq1grad), & 1333 & 0.0_dp 1334 1335 write(74,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1336 & qdrpwf_t4(im,iatpert,iq2grad,iq1grad)+qdrpwf_t4(im,iatpert,iiq2grad,iiq1grad), & 1337 & 0.0_dp 1338 1339 write(75,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1340 & qdrpwf_t5(im,iatpert,iq2grad,iq1grad)+qdrpwf_t5(im,iatpert,iiq2grad,iiq1grad), & 1341 & 0.0_dp 1342 1343 write(76,'(4(i2,2x),2(f18.10,2x))') iq2dir,iatom,iatdir,iq1dir, & 1344 & eqgradhart(im,iatpert,iq2grad,iq1grad)+eqgradhart(im,iatpert,iiq2grad,iiq1grad), & 1345 & 0.0_dp 1346 end if 1347 1348 end if 1349 end if 1350 end do 1351 end do 1352 end do 1353 1354 else 1355 write(msg,"(1a)") 'kptopt must be 2 or 3 for the quadrupole calculation' 1356 ABI_BUG(msg) 1357 end if 1358 1359 if (prtvol>=10) then 1360 close(71) 1361 close(72) 1362 close(73) 1363 close(74) 1364 close(75) 1365 close(76) 1366 end if 1367 1368 !Transformation to cartesian coordinates of the quadrupole tensor 1369 ABI_MALLOC(qdrptens_cart,(2,matom,3,3,3)) 1370 ABI_MALLOC(dqpol_cart,(2,matom,3,3,3)) 1371 ABI_MALLOC(cartflg,(matom,3,3,3)) 1372 qdrptens_cart(:,:,:,:,:)=qdrptens_red(:,:,:,:,:) 1373 dqpol_cart(:,:,:,:,:)=dqpol_red(:,:,:,:,:) 1374 cartflg=0 1375 1376 ABI_FREE(qdrptens_red) 1377 ABI_FREE(dqpol_red) 1378 1379 !1st transform coordinates of the atomic displacement derivative 1380 do iq1dir=1,3 1381 do iq2dir=1,3 1382 do ii=1,2 1383 do iatom=1,matom 1384 1385 do iatdir=1,3 1386 vec1(iatdir)=qdrptens_cart(ii,iatom,iatdir,iq2dir,iq1dir) 1387 flg1(iatdir)=qdrflg(iatom,iatdir,iq2dir,iq1dir) 1388 end do 1389 call cart39(flg1,flg2,gprimd,iatom,matom,rprimd,vec1,vec2) 1390 do iatdir=1,3 1391 qdrptens_cart(ii,iatom,iatdir,iq2dir,iq1dir)=vec2(iatdir) 1392 cartflg(iatom,iatdir,iq2dir,iq1dir)=flg2(iatdir) 1393 end do 1394 1395 do iatdir=1,3 1396 vec1(iatdir)=dqpol_cart(ii,iatom,iatdir,iq2dir,iq1dir) 1397 flg1(iatdir)=qdrflg(iatom,iatdir,iq2dir,iq1dir) 1398 end do 1399 call cart39(flg1,flg2,gprimd,iatom,matom,rprimd,vec1,vec2) 1400 do iatdir=1,3 1401 dqpol_cart(ii,iatom,iatdir,iq2dir,iq1dir)=vec2(iatdir) 1402 end do 1403 1404 end do 1405 end do 1406 end do 1407 end do 1408 1409 !2nd transform coordinates of the electric field derivative 1410 do iq1dir=1,3 1411 do iatdir=1,3 1412 do iatom=1,matom 1413 do ii=1,2 1414 do iq2dir=1,3 1415 vec1(iq2dir)=qdrptens_cart(ii,iatom,iatdir,iq2dir,iq1dir) 1416 flg1(iq2dir)=qdrflg(iatom,iatdir,iq2dir,iq1dir) 1417 end do 1418 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 1419 do iq2dir=1,3 1420 qdrptens_cart(ii,iatom,iatdir,iq2dir,iq1dir)=vec2(iq2dir) 1421 cartflg(iatom,iatdir,iq2dir,iq1dir)=flg2(iq2dir) 1422 end do 1423 1424 do iq2dir=1,3 1425 vec1(iq2dir)=dqpol_cart(ii,iatom,iatdir,iq2dir,iq1dir) 1426 flg1(iq2dir)=qdrflg(iatom,iatdir,iq2dir,iq1dir) 1427 end do 1428 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 1429 do iq2dir=1,3 1430 dqpol_cart(ii,iatom,iatdir,iq2dir,iq1dir)=vec2(iq2dir) 1431 end do 1432 1433 end do 1434 end do 1435 end do 1436 end do 1437 1438 !3rd transform coordinates of the q-gradient (treat it as electric field) 1439 do iq2dir=1,3 1440 do iatdir=1,3 1441 do iatom=1,matom 1442 do ii=1,2 1443 do iq1dir=1,3 1444 vec1(iq1dir)=qdrptens_cart(ii,iatom,iatdir,iq2dir,iq1dir) 1445 flg1(iq1dir)=qdrflg(iatom,iatdir,iq2dir,iq1dir) 1446 end do 1447 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 1448 do iq1dir=1,3 1449 qdrptens_cart(ii,iatom,iatdir,iq2dir,iq1dir)=vec2(iq1dir) 1450 cartflg(iatom,iatdir,iq2dir,iq1dir)=flg2(iq1dir) 1451 end do 1452 1453 do iq1dir=1,3 1454 vec1(iq1dir)=dqpol_cart(ii,iatom,iatdir,iq2dir,iq1dir) 1455 flg1(iq1dir)=qdrflg(iatom,iatdir,iq2dir,iq1dir) 1456 end do 1457 call cart39(flg1,flg2,gprimd,matom+2,matom,rprimd,vec1,vec2) 1458 do iq1dir=1,3 1459 dqpol_cart(ii,iatom,iatdir,iq2dir,iq1dir)=vec2(iq1dir) 1460 end do 1461 1462 end do 1463 end do 1464 end do 1465 end do 1466 1467 1468 !Write the Quadrupole tensor in cartesian coordinates 1469 write(ab_out,'(a)')' ' 1470 write(ab_out,'(a)')' Quadrupole tensor, in cartesian coordinates,' 1471 write(ab_out,'(a)')' efidir atom atddir qgrdir real part imaginary part' 1472 do iq1dir=1,3 1473 do iq2dir=1,3 1474 do iatdir=1,3 1475 do iatom=1,matom 1476 1477 if (cartflg(iatom,iatdir,iq2dir,iq1dir)==1) then 1478 1479 write(ab_out,'(4(i5,3x),2(1x,f20.10))') iq2dir,iatom,iatdir,iq1dir, & 1480 & qdrptens_cart(re,iatom,iatdir,iq2dir,iq1dir),qdrptens_cart(im,iatom,iatdir,iq2dir,iq1dir) 1481 1482 end if 1483 1484 end do 1485 end do 1486 end do 1487 write(ab_out,'(a)')' ' 1488 end do 1489 1490 !Write the q-gradient of the Polarization response 1491 write(ab_out,*)' First real-space moment of the polarization response ' 1492 write(ab_out,*)' to an atomic displacementatom, in cartesian coordinates,' 1493 write(ab_out,*)' (1/ucvol factor not included),' 1494 write(ab_out,*)' efidir atom atddir qgrdir real part imaginary part' 1495 do iq1dir=1,3 1496 do iq2dir=1,3 1497 do iatdir=1,3 1498 do iatom=1,matom 1499 1500 if (cartflg(iatom,iatdir,iq2dir,iq1dir)==1) then 1501 1502 write(ab_out,'(4(i5,3x),2(1x,f20.10))') iq2dir,iatom,iatdir,iq1dir, & 1503 & dqpol_cart(re,iatom,iatdir,iq2dir,iq1dir),dqpol_cart(im,iatom,iatdir,iq2dir,iq1dir) 1504 1505 end if 1506 1507 end do 1508 end do 1509 end do 1510 write(ab_out,*)' ' 1511 end do 1512 1513 !Write the electronic (frozen-ion) contribution to the piezoelectric tensor 1514 !(R.M. Martin, PRB 5, 1607 (1972)) 1515 write(ab_out,'(a)')' ' 1516 write(ab_out,'(a)')' Electronic (clamped-ion) contribution to the piezoelectric tensor,' 1517 write(ab_out,'(a)')' in cartesian coordinates, (from quadrupole calculation)' 1518 write(ab_out,'(a)')' atddir qgrdir efidir real part imaginary part' 1519 do iq1dir=1,3 1520 gamma=iq1dir 1521 do iq2dir=1,3 1522 alpha=iq2dir 1523 do iatdir=1,3 1524 beta=iatdir 1525 1526 piezore=zero 1527 piezoim=zero 1528 do iatom=1,matom 1529 1530 if (cartflg(iatom,iatdir,iq2dir,iq1dir)==1) then 1531 1532 piezore=piezore+( qdrptens_cart(re,iatom,beta,alpha,gamma)-qdrptens_cart(re,iatom,alpha,gamma,beta) & 1533 & + qdrptens_cart(re,iatom,gamma,beta,alpha) ) 1534 1535 piezoim=piezoim+( qdrptens_cart(im,iatom,beta,alpha,gamma)-qdrptens_cart(im,iatom,alpha,gamma,beta) & 1536 & + qdrptens_cart(im,iatom,gamma,beta,alpha) ) 1537 end if 1538 1539 end do 1540 1541 piezore=-piezore*half/ucvol 1542 piezoim=-piezoim*half/ucvol 1543 1544 write(ab_out,'(3(i5,3x),2(1x,f20.10))') beta,gamma,alpha,piezore,piezoim 1545 1546 end do 1547 end do 1548 write(ab_out,'(a)')' ' 1549 end do 1550 write(ab_out,'(80a)')('=',ii=1,80) 1551 1552 ABI_FREE(qdrptens_cart) 1553 ABI_FREE(dqpol_cart) 1554 ABI_FREE(cartflg) 1555 1556 DBG_EXIT("COLL") 1557 1558 end subroutine dfpt_qdrpout
ABINIT/m_dfpt_lw [ Modules ]
NAME
m_dfpt_lw
FUNCTION
Calculation of spatial dispertion magnitudes (quadrupole and flexoelectric tensor) from the DFPT long-wave approach.
COPYRIGHT
Copyright (C) 2018-2022 ABINIT group (MR,MS) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
SOURCE
20 #if defined HAVE_CONFIG_H 21 #include "config.h" 22 #endif 23 24 #include "abi_common.h" 25 26 module m_dfpt_lw 27 28 use defs_basis 29 use defs_datatypes 30 use defs_abitypes 31 use defs_wvltypes 32 use m_dtset 33 use m_dtfil 34 use m_profiling_abi 35 use m_xmpi 36 use m_errors 37 use m_wfk 38 use m_nctk 39 use m_hamiltonian 40 use m_ebands 41 use m_pawrhoij 42 use m_pawtab 43 use m_wffile 44 45 use m_hdr 46 use m_io_tools, only : file_exists 47 use m_ioarr, only : read_rhor, fftdatar_write_from_hdr 48 use m_rf2, only : rf2_getidirs 49 use m_kg, only : getcut, getph, getmpw, kpgio 50 use m_abicore, only : appdig 51 use m_fft, only : fourdp 52 use m_dfpt_rhotov, only : dfpt_rhotov 53 use m_dfpt_mkvxc, only : dfpt_mkvxcggadq 54 use m_spacepar, only : hartredq, setsym 55 use m_cgtools, only : dotprod_vn 56 use m_symkpt, only : symkpt 57 use m_mpinfo, only : distrb2, initmpi_band, proc_distrb_cycle 58 use m_initylmg, only : initylmg 59 use m_inwffil, only : inwffil 60 use m_dfpt_lwwf 61 use m_dynmat, only : cart39 62 63 implicit none 64 65 private