TABLE OF CONTENTS
ABINIT/m_paw_onsite [ Modules ]
NAME
m_paw_onsite
FUNCTION
This module contains a set of routines to compute various PAW on-site quantities i.e. quantities expressed with <Phi_i|...|Phi_j> and/or <tild_Phi_i|...|tild_Phi_j>.
COPYRIGHT
Copyright (C) 2013-2022 ABINIT group (MT,FJ) 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
FOR DEVELOPPERS: in order to preserve the portability of libPAW library, please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
SOURCE
21 #include "libpaw.h" 22 23 MODULE m_paw_onsite 24 25 USE_DEFS 26 USE_MSG_HANDLING 27 USE_MEMORY_PROFILING 28 29 use m_pawrad, only : pawrad_type, pawrad_deducer0, simp_gen, nderiv_gen 30 use m_pawtab, only : pawtab_type 31 use m_paw_sphharm, only : setnabla_ylm 32 33 implicit none 34 35 private 36 37 !public procedures. 38 public :: pawnabla_init ! Evaluate valence-valence on-site contribs of the nabla operator in cart. coord. 39 public :: pawnabla_core_init ! Evaluate core-valence on-site contribs of the nabla operator in cart. coord.
m_paw_onsite/pawnabla_core_init [ Functions ]
[ Top ] [ m_paw_onsite ] [ Functions ]
NAME
pawnabla_core_init
FUNCTION
Evaluate core-valence onsite contributions of the nabla operator in cartesian coordinates, i.e. <Phi_i|Nabla|Phi_core_j>-<tPhi_i|Nabla|tPhi_core_j>. Core wave-functions are only given for one atom type.
INPUTS
mpsang=1+maximum angular momentum ntypat=Number of types of atoms in cell Pawrad(ntypat)<Pawrad_type>=PAW radial mesh and related data: %rad(mesh_size)=The coordinates of all the points of the radial mesh Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data: %mesh_size=Dimension of radial mesh %lmn_size=Number of (l,m,n) elements for the PAW basis phi_cor(mesh_size,nphicor)=core wave-functions for the current type of atoms. indlmn_cor(6,nlmn_core)= array giving l,m,n,lm,ln,s for i=lmn, for the core wave functions.
OUTPUT
See side effects
SIDE EFFECTS
Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data: %has_nabla=set to 1 in matrix elements are calculated and stored %nabla_ij(3,lmn_size,lmn_size)= <phi_i|nabla|phi_core_j>-<tphi_i|nabla|tphi_core_j>
NOTES
MG extracted this piece of code from optics_paw.F90 in order to have something more reusable! Note however the storage mode of nabla_ij differs from optics_paw (here Cartesian coordinates run faster). Besides nabla_ij contains the matrix elements of \nabla instead of the elements of the momentum operator p.
SOURCE
261 subroutine pawnabla_core_init(mpsang,ntypat,pawrad,pawtab,phi_cor,indlmn_cor,diracflag) 262 263 !Arguments ------------------------------------ 264 !scalars 265 integer,intent(in) :: mpsang,ntypat 266 integer,intent(in),optional :: diracflag 267 !arrays 268 integer,intent(in) :: indlmn_cor(:,:) 269 real(dp),intent(in) :: phi_cor(:,:) 270 type(pawtab_type),target,intent(inout) :: pawtab(ntypat) 271 type(pawrad_type),intent(in) :: pawrad(ntypat) 272 273 !Local variables------------------------------- 274 !scalars 275 integer :: nln,nln_cor,ilm,ilmn,iln,itypat,sgnkappa 276 integer :: jl,jm,jlm,jlmn,jln,js,jlm_re,jlm_im,jm_re,jm_im 277 integer :: lmn_size,lmncmax,lcmax,ltmax,mesh_size,mesh_size_cor 278 real(dp) :: intg,jmj,cgc 279 logical :: dirac 280 character(len=500) :: msg 281 !arrays 282 integer, LIBPAW_CONTIGUOUS pointer :: indlmn(:,:) 283 real(dp) , allocatable:: ang_phipphj(:,:,:) 284 real(dp),allocatable :: dphi(:),ff(:),int1(:,:),int2(:,:),rad(:) 285 286 ! ************************************************************************* 287 288 mesh_size_cor=size(phi_cor,1) 289 nln_cor=size(phi_cor,2) 290 291 dirac=.false. 292 if(present(diracflag)) then 293 dirac=(diracflag==1) 294 if(diracflag==1.and.size(indlmn_cor,1)<8) then 295 msg='Wrong first dimension of indlmn_cor in pawnabla_core_init for diracrelativism!' 296 LIBPAW_BUG(msg) 297 endif 298 endif 299 300 !To be checked 301 lmncmax=size(indlmn_cor,2) !Includes spinors if diracrel core wf are used 302 lcmax=maxval(indlmn_cor(1,:)) 303 ltmax=max(lcmax+1,mpsang) 304 LIBPAW_ALLOCATE(ang_phipphj,(ltmax**2,ltmax**2,8)) 305 306 if (ltmax>4)then 307 write(msg,'(3a)')& 308 & 'Not designed for angular momentum greater than 3!',ch10,& 309 & 'Modification in the table defined in routine setnabla_ylm is required.' 310 LIBPAW_BUG(msg) 311 end if 312 313 !if (mesh_size_cor/=pawrad(1)%mesh_size) then 314 ! write(msg,'(a)') 'Wrong mesh_size_cor value (1)!' 315 ! LIBPAW_BUG(msg) 316 !end if 317 !if (any(mesh_size_cor/=pawtab(:)%mesh_size)) then 318 ! write(msg,'(3a)') 'Wrong mesh_size_cor value (2)!',ch10,& 319 !& 'Should have only one type of atom.' 320 ! LIBPAW_ERROR(msg) 321 !end if 322 323 !Integration of the angular part: all angular integrals have been computed 324 !outside Abinit and tabulated for each (l,m) value up to l=3 325 call setnabla_ylm(ang_phipphj,ltmax) 326 327 do itypat=1,ntypat 328 329 ! COMPUTE nabla_ij := <phi_i|nabla|phi_cor> for this type 330 mesh_size=min(pawtab(itypat)%partialwave_mesh_size,pawrad(itypat)%mesh_size) 331 mesh_size=min(mesh_size_cor,mesh_size) 332 lmn_size=pawtab(itypat)%lmn_size 333 nln=pawtab(itypat)%basis_size 334 335 if (mesh_size_cor<mesh_size) then 336 msg='mesh_size and mesh_sier_cor not compatible!' 337 LIBPAW_BUG(msg) 338 endif 339 340 if (allocated(pawtab(itypat)%nabla_ij)) then 341 LIBPAW_DEALLOCATE(pawtab(itypat)%nabla_ij) 342 end if 343 LIBPAW_ALLOCATE(pawtab(itypat)%nabla_ij,(3,lmn_size,lmncmax)) 344 345 if (dirac) then 346 if (allocated(pawtab(itypat)%nabla_im_ij)) then 347 LIBPAW_DEALLOCATE(pawtab(itypat)%nabla_im_ij) 348 end if 349 LIBPAW_ALLOCATE(pawtab(itypat)%nabla_im_ij,(3,lmn_size,lmncmax)) 350 end if 351 352 pawtab(itypat)%has_nabla=1 353 354 LIBPAW_ALLOCATE(ff,(mesh_size)) 355 LIBPAW_ALLOCATE(rad,(mesh_size)) 356 LIBPAW_ALLOCATE(int1,(lmn_size,lmncmax)) 357 LIBPAW_ALLOCATE(int2,(lmn_size,lmncmax)) 358 LIBPAW_ALLOCATE(dphi,(mesh_size)) 359 360 indlmn => pawtab(itypat)%indlmn 361 rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size) 362 363 ! int1= \int Phi d/dr(Phi_core) r^2 dr 364 ! = \int (phi d/dr(phi_core) - phi phj_core/r) dr 365 ! with Phi=phi/r and Phi_core=phi_core/r 366 do jln=1,nln_cor 367 ff(1:mesh_size)=phi_cor(1:mesh_size,jln) 368 call nderiv_gen(dphi,ff,pawrad(itypat)) 369 do iln=1,nln 370 ff(2:mesh_size)=pawtab(itypat)%phi(2:mesh_size,iln)*dphi(2:mesh_size) & 371 & -pawtab(itypat)%phi(2:mesh_size,iln)*phi_cor(2:mesh_size,jln)/rad(2:mesh_size) 372 call pawrad_deducer0(ff,mesh_size,pawrad(itypat)) 373 call simp_gen(intg,ff,pawrad(itypat)) 374 int1(iln,jln)=intg 375 end do 376 end do 377 378 ! int2= \int Phi Phi_core /r r^2 dr = \int phi phi_core /r dr 379 ! with Phi=phi/r and Phi_core=phi_core/r 380 do jln=1,nln_cor 381 do iln=1,nln 382 ff(2:mesh_size)=(pawtab(itypat)%phi(2:mesh_size,iln)*phi_cor(2:mesh_size,jln))/rad(2:mesh_size) 383 call pawrad_deducer0(ff,mesh_size,pawrad(itypat)) 384 call simp_gen(intg,ff,pawrad(itypat)) 385 int2(iln,jln)=intg 386 end do 387 end do 388 389 ! ===== FULLY-RELATIVISTIC ===== 390 if(dirac) then 391 392 ! Integration of the radial part, Note unpacked loop 393 do jlmn=1,lmncmax 394 jl=indlmn_cor(1,jlmn) 395 jm=indlmn_cor(2,jlmn) 396 397 sgnkappa=indlmn_cor(3,jlmn) 398 jmj=half*indlmn_cor(8,jlmn) ! 2mj is stored in indlmn_cor 399 js=indlmn_cor(6,jlmn) ! 1 is up, 2 is down 400 401 ! Calculate spinor dependend coeffs 402 ! (Clebsch-Gordan, I guess) 403 cgc=one ! so nothing changes without core spinors 404 if (sgnkappa==1) then 405 if(js==1) then 406 cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1)) 407 else 408 cgc=-sqrt((dble(jl)+jmj+half)/dble(2*jl+1)) 409 endif 410 else 411 if(js==1) then 412 cgc= sqrt((dble(jl)+jmj+half)/dble(2*jl+1)) 413 else 414 cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1)) 415 endif 416 endif 417 418 jlm=indlmn_cor(4,jlmn) 419 jln=indlmn_cor(5,jlmn) 420 421 do ilmn=1,lmn_size 422 ilm=indlmn(4,ilmn) 423 iln=indlmn(5,ilmn) 424 425 ! jl was set as a flag for invalid combinations 426 ! i.e. m=-(l+1) or m=(l+1) 427 ! In these cases, cgc=0 ; so nabla_ij=0 428 if(jl==-1) then 429 pawtab(itypat)%nabla_ij(1:3,ilmn,jlmn)= zero 430 pawtab(itypat)%nabla_im_ij(1:3,ilmn,jlmn) = zero 431 432 else 433 434 ! if jm<>0, need to convert from complex 435 ! to real spherical harmonics 436 if(jm<0) then 437 jm_re=abs(jm) 438 jm_im=-abs(jm) 439 jlm_re=jl*(jl+1)+jm_re+1 440 jlm_im=jl*(jl+1)+jm_im+1 441 pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=half_sqrt2*cgc*( & 442 & int1(iln,jln)* ang_phipphj(ilm,jlm_re,1) & 443 & +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,2)+ang_phipphj(ilm,jlm_re,3))) 444 pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=half_sqrt2*cgc*( & 445 & int1(iln,jln)* ang_phipphj(ilm,jlm_re,4) & 446 & +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,5)+ang_phipphj(ilm,jlm_re,6))) 447 pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=half_sqrt2*cgc*( & 448 & int1(iln,jln)* ang_phipphj(ilm,jlm_re,7) & 449 & +int2(iln,jln)* ang_phipphj(ilm,jlm_re,8)) 450 pawtab(itypat)%nabla_im_ij(1,ilmn,jlmn)=-half_sqrt2*cgc*( & 451 & int1(iln,jln)* ang_phipphj(ilm,jlm_im,1) & 452 & +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,2)+ang_phipphj(ilm,jlm_im,3))) 453 pawtab(itypat)%nabla_im_ij(2,ilmn,jlmn)=-half_sqrt2*cgc*( & 454 & int1(iln,jln)* ang_phipphj(ilm,jlm_im,4) & 455 & +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,5)+ang_phipphj(ilm,jlm_im,6))) 456 pawtab(itypat)%nabla_im_ij(3,ilmn,jlmn)=-half_sqrt2*cgc*( & 457 & int1(iln,jln)* ang_phipphj(ilm,jlm_im,7) & 458 & +int2(iln,jln)* ang_phipphj(ilm,jlm_im,8)) 459 460 else if (jm>0) then 461 jm_re=abs(jm) 462 jm_im=-abs(jm) 463 jlm_re=jl*(jl+1)+jm_re+1 464 jlm_im=jl*(jl+1)+jm_im+1 465 pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( & 466 & int1(iln,jln)* ang_phipphj(ilm,jlm_re,1) & 467 & +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,2)+ang_phipphj(ilm,jlm_re,3))) 468 pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( & 469 & int1(iln,jln)* ang_phipphj(ilm,jlm_re,4) & 470 & +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,5)+ang_phipphj(ilm,jlm_re,6))) 471 pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( & 472 & int1(iln,jln)* ang_phipphj(ilm,jlm_re,7) & 473 & +int2(iln,jln)* ang_phipphj(ilm,jlm_re,8)) 474 pawtab(itypat)%nabla_im_ij(1,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( & 475 & int1(iln,jln)* ang_phipphj(ilm,jlm_im,1) & 476 & +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,2)+ang_phipphj(ilm,jlm_im,3))) 477 pawtab(itypat)%nabla_im_ij(2,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( & 478 & int1(iln,jln)* ang_phipphj(ilm,jlm_im,4) & 479 & +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,5)+ang_phipphj(ilm,jlm_im,6))) 480 pawtab(itypat)%nabla_im_ij(3,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( & 481 & int1(iln,jln)* ang_phipphj(ilm,jlm_im,7) & 482 & +int2(iln,jln)* ang_phipphj(ilm,jlm_im,8)) 483 484 else ! jm=0 : no conversion necessary if m=0 485 pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=cgc*( & 486 & int1(iln,jln)* ang_phipphj(ilm,jlm,1) & 487 & +int2(iln,jln)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3))) 488 pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=cgc*( & 489 & int1(iln,jln)* ang_phipphj(ilm,jlm,4) & 490 & +int2(iln,jln)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6))) 491 pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=cgc*( & 492 & int1(iln,jln)* ang_phipphj(ilm,jlm,7) & 493 & +int2(iln,jln)* ang_phipphj(ilm,jlm,8)) 494 pawtab(itypat)%nabla_im_ij(1,ilmn,jlmn)= zero 495 pawtab(itypat)%nabla_im_ij(2,ilmn,jlmn)= zero 496 pawtab(itypat)%nabla_im_ij(3,ilmn,jlmn)= zero 497 end if 498 499 endif ! jl==-1? 500 501 end do !ilmn 502 end do !jlmn 503 504 pawtab(itypat)%has_nabla=4 505 506 ! ===== NON-RELATIVISTIC OR SCALAR-RELATICISTIC ===== 507 else 508 509 ! Integration of the radial part, Note unpacked loop 510 do jlmn=1,lmncmax 511 jl=indlmn_cor(1,jlmn) 512 jlm=indlmn_cor(4,jlmn) 513 jln =indlmn_cor(5,jlmn) 514 do ilmn=1,lmn_size 515 ilm=indlmn(4,ilmn) 516 iln =indlmn(5,ilmn) 517 pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=( & 518 & int1(iln,jln)* ang_phipphj(ilm,jlm,1) & 519 & +int2(iln,jln)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3))) 520 521 pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=( & 522 & int1(iln,jln)* ang_phipphj(ilm,jlm,4) & 523 & +int2(iln,jln)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6))) 524 525 pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=( & 526 & int1(iln,jln)* ang_phipphj(ilm,jlm,7) & 527 & +int2(iln,jln)* ang_phipphj(ilm,jlm,8)) 528 end do !ilmn 529 end do !jlmn 530 531 pawtab(itypat)%has_nabla=3 532 533 end if ! Relativistic? 534 535 LIBPAW_DEALLOCATE(ff) 536 LIBPAW_DEALLOCATE(rad) 537 LIBPAW_DEALLOCATE(int1) 538 LIBPAW_DEALLOCATE(int2) 539 LIBPAW_DEALLOCATE(dphi) 540 541 end do !itypat 542 543 LIBPAW_DEALLOCATE(ang_phipphj) 544 545 end subroutine pawnabla_core_init
m_paw_onsite/pawnabla_init [ Functions ]
[ Top ] [ m_paw_onsite ] [ Functions ]
NAME
pawnabla_init
FUNCTION
Evaluate all valence-valence onsite contributions of the nabla operator in cartesian coordinates, i.e. <Phi_i|Nabla|Phi_j>-<tPhi_i|Nabla|tPhi_j>.
INPUTS
mpsang=1+maximum angular momentum ntypat=Number of types of atoms in cell Pawrad(ntypat)<Pawrad_type>=PAW radial mesh and related data: %rad(mesh_size)=The coordinates of all the points of the radial mesh Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data: %mesh_size=Dimension of radial mesh %lmn_size=Number of (l,m,n) elements for the PAW basis
OUTPUT
See side effects
SIDE EFFECTS
Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data: %has_nabla=set to 1 in matrix elements are calculated and stored %nabla_ij(3,lmn_size,lmn_size)= <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j>
NOTES
MG extracted this piece of code from optics_paw.F90 in order to have something more reusable! Note however the storage mode of nabla_ij differs from optics_paw (here Cartesian coordinates run faster). Besides nabla_ij contains the matrix elements of \nabla instead of the elements of the momentum operator p.
SOURCE
83 subroutine pawnabla_init(mpsang,ntypat,pawrad,pawtab) 84 85 !Arguments ------------------------------------ 86 !scalars 87 integer,intent(in) :: mpsang,ntypat 88 !arrays 89 type(pawtab_type),target,intent(inout) :: pawtab(ntypat) 90 type(pawrad_type),intent(in) :: pawrad(ntypat) 91 92 !Local variables------------------------------- 93 !scalars 94 integer :: ii,nln,il,ilm,ilmn,iln,itypat 95 integer :: jl,jlm,jlmn,jln,lmn_size,mesh_size 96 real(dp) :: avg,intg 97 character(len=500) :: msg 98 !arrays 99 integer, LIBPAW_CONTIGUOUS pointer :: indlmn(:,:) 100 real(dp) :: ang_phipphj(mpsang**2,mpsang**2,8) 101 real(dp),allocatable :: dphi(:),dtphi(:),ff(:),int1(:,:),int2(:,:),rad(:) 102 103 ! ************************************************************************* 104 105 if (mpsang>4)then 106 write(msg,'(3a)')& 107 & 'Not designed for angular momentum greater than 3 ',ch10,& 108 & 'Modification in the table defined in routine setnabla_ylm is required.' 109 LIBPAW_BUG(msg) 110 end if 111 112 !Integration of the angular part: all angular integrals have been computed 113 !outside Abinit and tabulated for each (l,m) value up to l=3 114 call setnabla_ylm(ang_phipphj,mpsang) 115 116 do itypat=1,ntypat 117 118 ! COMPUTE nabla_ij := <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for this type 119 mesh_size=pawtab(itypat)%mesh_size 120 lmn_size=pawtab(itypat)%lmn_size 121 nln=pawtab(itypat)%basis_size 122 123 if (allocated(pawtab(itypat)%nabla_ij)) then 124 LIBPAW_DEALLOCATE(pawtab(itypat)%nabla_ij) 125 end if 126 LIBPAW_ALLOCATE(pawtab(itypat)%nabla_ij,(3,lmn_size,lmn_size)) 127 pawtab(itypat)%has_nabla=1 128 129 LIBPAW_ALLOCATE(ff,(mesh_size)) 130 LIBPAW_ALLOCATE(rad,(mesh_size)) 131 LIBPAW_ALLOCATE(int1,(lmn_size,lmn_size)) 132 LIBPAW_ALLOCATE(int2,(lmn_size,lmn_size)) 133 LIBPAW_ALLOCATE(dphi,(mesh_size)) 134 LIBPAW_ALLOCATE(dtphi,(mesh_size)) 135 136 indlmn => pawtab(itypat)%indlmn 137 rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size) 138 139 ! int1= \int [ Phi d/dr(Phj) - tPhi d/dr(tPhj) ] r^2 dr 140 ! = \int [ (phi d/dr(phj) - phi phj/r) - (tphi d/dr(tphj) - tphi tphj/r) ] dr 141 ! with Phi=phi/r and tPhi=phi/r 142 do jln=1,nln 143 ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,jln) 144 call nderiv_gen(dphi,ff,pawrad(itypat)) 145 ff(1:mesh_size)=pawtab(itypat)%tphi(1:mesh_size,jln) 146 call nderiv_gen(dtphi,ff,pawrad(itypat)) 147 do iln=1,nln 148 ff(2:mesh_size)= & 149 & pawtab(itypat)%phi (2:mesh_size,iln)*dphi (2:mesh_size) & 150 & -pawtab(itypat)%phi (2:mesh_size,iln)*pawtab(itypat)%phi (2:mesh_size,jln)/rad(2:mesh_size) & 151 & -( pawtab(itypat)%tphi(2:mesh_size,iln)*dtphi(2:mesh_size) & 152 & -pawtab(itypat)%tphi(2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln)/rad(2:mesh_size) ) 153 call pawrad_deducer0(ff,mesh_size,pawrad(itypat)) 154 call simp_gen(intg,ff,pawrad(itypat)) 155 int1(iln,jln)=intg 156 end do 157 end do 158 159 ! int2= \int [ Phi Phj /r - \int tPhi tPhj /r ] r^2 dr 160 ! = \int [ phi phj /r - \int tphi tphj /r ] dr 161 ! with Phi=phi/r and tPhi=phi/r 162 do jln=1,nln 163 do iln=1,nln 164 ff(2:mesh_size)= ( & 165 & pawtab(itypat)%phi (2:mesh_size,iln)*pawtab(itypat)%phi (2:mesh_size,jln) & 166 & -pawtab(itypat)%tphi(2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln) ) /rad(2:mesh_size) 167 call pawrad_deducer0(ff,mesh_size,pawrad(itypat)) 168 call simp_gen(intg,ff,pawrad(itypat)) 169 int2(iln,jln)=intg 170 end do 171 end do 172 173 ! Integration of the radial part, Note unpacked loop 174 do jlmn=1,lmn_size 175 jlm=indlmn(4,jlmn) 176 jl =indlmn(5,jlmn) 177 do ilmn=1,lmn_size 178 ilm=indlmn(4,ilmn) 179 il =indlmn(5,ilmn) 180 181 pawtab(itypat)%nabla_ij(1,ilmn,jlmn)= & 182 & int1(il,jl)* ang_phipphj(ilm,jlm,1) & 183 & +int2(il,jl)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3)) 184 185 pawtab(itypat)%nabla_ij(2,ilmn,jlmn)= & 186 & int1(il,jl)* ang_phipphj(ilm,jlm,4) & 187 & +int2(il,jl)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6)) 188 189 pawtab(itypat)%nabla_ij(3,ilmn,jlmn)= & 190 & int1(il,jl)* ang_phipphj(ilm,jlm,7) & 191 & +int2(il,jl)* ang_phipphj(ilm,jlm,8) 192 193 end do !ilmn 194 end do !jlmn 195 196 ! Symetrization 197 if (lmn_size>1) then 198 do jlmn=2,lmn_size 199 do ilmn=1,jlmn-1 200 do ii=1,3 201 avg=half*(pawtab(itypat)%nabla_ij(ii,ilmn,jlmn)-pawtab(itypat)%nabla_ij(ii,jlmn,ilmn)) 202 pawtab(itypat)%nabla_ij(ii,ilmn,jlmn)= avg 203 pawtab(itypat)%nabla_ij(ii,jlmn,ilmn)=-avg 204 end do 205 end do 206 end do 207 end if 208 209 ! End 210 pawtab(itypat)%has_nabla=2 211 LIBPAW_DEALLOCATE(ff) 212 LIBPAW_DEALLOCATE(rad) 213 LIBPAW_DEALLOCATE(int2) 214 LIBPAW_DEALLOCATE(int1) 215 LIBPAW_DEALLOCATE(dphi) 216 LIBPAW_DEALLOCATE(dtphi) 217 218 end do !itypat 219 220 end subroutine pawnabla_init