1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8module trialorbitalclass 9 10! Includes: 11! Derived types: 12! trialorbital 13! Variables: 14! cutoff: cutoff radii 15! T: tolerances controling the cutoff radii 16! Procedures: 17! function gettrialwavefunction 18! real(dp) function gettrialrcut( ofwhat ) 19! integer function gettriallmax( ofwhat ) 20! function gettrialcenter( ofwhat ) 21! subroutine print_trialorb 22 23use precision, only:dp 24 25implicit none 26 27type trialorbital 28 real(dp),dimension(3) :: center ! Projection function center in 29 ! cristallographic coordinates relative 30 ! to the direct lattice vectors. 31 real(dp),dimension(3) :: zaxis ! Defines the axis from which the polar 32 ! angle theta in spherical polar coordinates 33 ! is measured. 34 ! Default: (0.0 0.0 1.0) 35 real(dp),dimension(3) :: xaxis ! Defines the axis from which the azimuthal 36 ! angle phi in spherical coordinates is 37 ! measured. 38 ! Must be orthogonal to z-axis. 39 real(dp),dimension(3) :: yaxis ! Angular momentum y-axis 40 real(dp) :: zovera ! z/a, diffusivity, spread. 41 ! Read from the nnkp file in Ang^-1 42 ! Transformed later to Ang^-1 43 integer :: r ! Radial quantum number 44 integer :: l ! Angular momentum 45 integer :: mr ! z-projection quantum number 46 real(dp) :: rcut ! Siesta's cut-off radius: Bohr 47 integer :: lmax ! Maximum total angular momentum 48end type 49 50! Cut-off radii in units of \alpha^-1 51real(dp),parameter,dimension(3) :: & 52 cutoffs = (/6.934_dp,18.87_dp,35.44_dp/) 53! Squared norm tolerance governing the cut-off radii 54real(dp),parameter :: T = 0.0001_dp 55 56CONTAINS 57 58! 59!<-----------------------WAVE FUNCTIONS-----------------------------> 60! 61real(dp) function gettrialwavefunction( orbital, atpoint ) 62! 63! Trial orbital as defined in Wannier90. 64! See pages 35-38 of the Wannier90 Users Guide, Version 1.2 65! Yields the function value at a given point relative to its center. 66! It contains the library of Wannier90 trial orbitals in the form 67! of statement functions. 68! 69! Argument: coordinates of the point in Bohrs 70! Output unit: Bohr^{-3/2} 71! 72 73 use units, only: pi ! Value of pi 74 75 implicit none 76! 77! Passed arguments 78! 79 80 real(dp), dimension(3), intent(in) :: atpoint 81 type(trialorbital), intent(in) :: orbital 82 83 real(dp),dimension(3) :: arg 84! Angular-dependent factor of the wave-function 85 real(dp) :: angular 86! 87! Statement functions: declaration 88! 89! x, y, and z coordinates 90 real(dp) :: x, y, z 91! Distance with respect to the origin 92 real(dp) :: rr 93! Spherical function 94 real(dp) :: sphere 95! Angular functions associated with particular values of l (trialorbital%l), 96! and m_r (trialorbital%mr) for l >= 0 97! s orbital 98 real(dp) :: s 99! p orbitals 100 real(dp) :: px, py, pz 101! d orbitals 102 real(dp) :: dz2, dxz, dyz, dxy, dx2y2 103! f orbitals 104 real(dp) :: fz3, fxz2, fyz2, fzx2y2, fxyz, fxx23y2, fy3x2y2 105! Angular functions associated with particular values of l (trialorbital%l), 106! and m_r (trialorbital%mr) for l < 0 (hybrid functions) 107! sp hybrids 108 real(dp) :: sp_1, sp_2 109! sp2 hybrids 110 real(dp) :: sp2_1, sp2_2, sp2_3 111! sp3 hybrids 112 real(dp) :: sp3_1, sp3_2, sp3_3, sp3_4 113! sp3d hybrids 114 real(dp) :: sp3d_1, sp3d_2, sp3d_3, sp3d_4, sp3d_5 115! sp3d2 hybrids 116 real(dp) :: sp3d2_1, sp3d2_2, sp3d2_3, sp3d2_4, sp3d2_5, sp3d2_6 117 integer :: rank 118! Radial functions associated with different values of the r value 119! defined in (trialorbital%r) 120 real(dp) :: R1, R2, R3 121 122! 123! Inverse square roots 124! 125 real(dp), parameter :: rs2 = 1.0_dp/dsqrt(2.0_dp) 126 real(dp), parameter :: rs3 = 1.0_dp/dsqrt(3.0_dp) 127 real(dp), parameter :: rs6 = 1.0_dp/dsqrt(6.0_dp) 128 real(dp), parameter :: rs12 = 1.0_dp/dsqrt(12.0_dp) 129! 130! Constants depending on l and power of z 131! See the prefactors before the angular functions in page 36 132! of the Wannier90 Users Guide, Version 1.2 133! 134 real(dp), parameter :: l0norm = 1.0_dp/dsqrt(4.0_dp*pi) 135 real(dp), parameter :: l1norm = dsqrt(3.0_dp)*l0norm 136 real(dp), parameter :: l2z2norm = dsqrt((5.0_dp/16.0_dp)/pi) 137 real(dp), parameter :: l2z1norm = dsqrt((15.0_dp/4.0_dp)/pi) 138 real(dp), parameter :: l2z0norm = l2z1norm*0.5_dp 139 real(dp), parameter :: l3z3norm = dsqrt(7.0_dp/pi)/4.0_dp 140 real(dp), parameter :: l3z2norm = dsqrt(21.0_dp/2.0_dp/pi)/4.0_dp 141 real(dp), parameter :: l3z1norm = dsqrt(105.0_dp/pi)/4.0_dp 142 real(dp), parameter :: l3z0norm = dsqrt(35.0_dp/2.0_dp/pi)/4.0_dp 143! 144! Real Spherical Harmonics 145! 146! To get a dimensionless spherical harmonics 147! 148 sphere(x,y,z,rank) = sqrt( x**2 + y**2 + z**2 )**rank 149! 150! s-orbital 151! 152 s(x,y,z) = l0norm 153! 154! p-orbitals 155! 156 px(x,y,z) = l1norm * x / sphere(x,y,z,1) 157 py(x,y,z) = l1norm * y / sphere(x,y,z,1) 158 pz(x,y,z) = l1norm * z / sphere(x,y,z,1) 159! 160! d-orbitals 161! 162 dz2(x,y,z) = l2z2norm * (2.0_dp * z**2 - x**2 - y**2 )/sphere(x,y,z,2) 163 dxz(x,y,z) = l2z1norm * z * x / sphere(x,y,z,2) 164 dyz(x,y,z) = l2z1norm * z * y / sphere(x,y,z,2) 165 dx2y2(x,y,z) = l2z0norm * (x**2 - y**2)/ sphere(x,y,z,2) 166 dxy(x,y,z) = l2z1norm * x * y / sphere(x,y,z,2) 167! 168! f-orbitals 169! 170 fz3(x,y,z) = l3z3norm * & 171 & (2.0_dp*z**2-3.0_dp*x**2-3.0_dp*y**2)*z/sphere(x,y,z,3) 172 fxz2(x,y,z) = l3z2norm * & 173 & (4.0_dp*z**2-x**2-y**2)*x/sphere(x,y,z,3) 174 fyz2(x,y,z) = l3z2norm * & 175 & (4.0_dp*z**2-x**2-y**2)*y/sphere(x,y,z,3) 176 fzx2y2(x,y,z) = l3z1norm * & 177 & z*(x**2-y**2)/sphere(x,y,z,3) 178 fxyz(x,y,z) = l3z1norm * & 179 & z*x*y*2.0_dp/sphere(x,y,z,3) 180 fxx23y2(x,y,z) = l3z0norm * & 181 & (x**2-3.0_dp*y**2)*x/sphere(x,y,z,3) 182 fy3x2y2(x,y,z) = l3z0norm * & 183 & (3.0_dp*x**2-y**2)*y/sphere(x,y,z,3) 184! 185! Hybrids. They follow the definitions given in page 37 186! of the Wannier90 Users Guide, Version 1.2 187! 188! sp hybrids 189! 190 sp_1(x,y,z) = rs2 * s(x,y,z) + rs2 * px(x,y,z) 191 sp_2(x,y,z) = rs2 * s(x,y,z) - rs2 * px(x,y,z) 192! 193! sp2 hybrids 194! 195 sp2_1(x,y,z) = rs3 * s(x,y,z) - rs6 * px(x,y,z) + rs2 * py(x,y,z) 196 sp2_2(x,y,z) = rs3 * s(x,y,z) - rs6 * px(x,y,z) - rs2 * py(x,y,z) 197 sp2_3(x,y,z) = rs3 * s(x,y,z) + rs6 * px(x,y,z) * 2.0_dp 198! 199! sp3 hybrids 200! 201 sp3_1(x,y,z) = 0.5_dp * ( s(x,y,z) + px(x,y,z) + py(x,y,z) + pz(x,y,z) ) 202 sp3_2(x,y,z) = 0.5_dp * ( s(x,y,z) + px(x,y,z) - py(x,y,z) - pz(x,y,z) ) 203 sp3_3(x,y,z) = 0.5_dp * ( s(x,y,z) - px(x,y,z) + py(x,y,z) - pz(x,y,z) ) 204 sp3_4(x,y,z) = 0.5_dp * ( s(x,y,z) - px(x,y,z) - py(x,y,z) + pz(x,y,z) ) 205! 206! sp3d hybrids 207! 208 sp3d_1(x,y,z) = rs3 * s(x,y,z) - rs6 * px(x,y,z) + rs2 * py(x,y,z) 209 sp3d_2(x,y,z) = rs3 * s(x,y,z) - rs6 * px(x,y,z) - rs2 * py(x,y,z) 210 sp3d_3(x,y,z) = rs3 * s(x,y,z) + rs6 * px(x,y,z) * 2.0_dp 211 sp3d_4(x,y,z) = rs2 * pz(x,y,z) + rs2 * dz2(x,y,z) 212 sp3d_5(x,y,z) = -rs2 * pz(x,y,z) + rs2 * dz2(x,y,z) 213! 214! sp3d2 hybrids 215! (Bug corrected by J. Junquera in the definition of sp3d2_3 and sp3d2_4. 216! In the original version by R. Korytar the px orbital was wrongly used 217! instead of the correct py orbital) 218! 219 sp3d2_1(x,y,z) = rs6 * s(x,y,z) - rs2 * px(x,y,z) - rs12 * dz2(x,y,z) + & 220 & dx2y2(x,y,z) / 2.0_dp 221 sp3d2_2(x,y,z) = rs6 * s(x,y,z) + rs2 * px(x,y,z) - rs12 * dz2(x,y,z) + & 222 & dx2y2(x,y,z) / 2.0_dp 223 sp3d2_3(x,y,z) = rs6 * s(x,y,z) - rs2 * py(x,y,z) - rs12 * dz2(x,y,z) - & 224 & dx2y2(x,y,z) / 2.0_dp 225 sp3d2_4(x,y,z) = rs6 * s(x,y,z) + rs2 * py(x,y,z) - rs12 * dz2(x,y,z) - & 226 & dx2y2(x,y,z) / 2.0_dp 227 sp3d2_5(x,y,z) = rs6 * s(x,y,z) - rs2 * pz(x,y,z) + rs3 * dz2(x,y,z) 228 sp3d2_6(x,y,z) = rs6 * s(x,y,z) + rs2 * pz(x,y,z) + rs3 * dz2(x,y,z) 229! 230! Radial part 231! They follow the definitions given in page 38 232! of the Wannier90 Users Guide, Version 1.2 233! 234 R1(rr) = 2.0_dp * orbital%zovera**(3.0_dp/2.0_dp) * & 235 & exp( -orbital%zovera * rr ) 236 R2(rr) = 0.5_dp/sqrt(2.0_dp) * orbital%zovera**(3.0_dp/2.0_dp) * & 237 & (2.0_dp - orbital%zovera * rr ) * exp( -orbital%zovera * rr / 2.0_dp) 238 R3(rr) = sqrt(4.0_dp/27.0_dp) * orbital%zovera**(3.0_dp/2.0_dp) * & 239 & ( 1.0_dp - 2.0_dp * orbital%zovera * rr/3.0_dp + & 240 & 2.0_dp * orbital%zovera**2 * rr**2 / 27.0_dp ) * & 241 & exp( -orbital%zovera * rr / 3.0_dp ) 242 243! 244! Executables !-----------> 245! 246! argument is the relative position of a given point with respect the center 247! of the trial function 248! arg = atpoint - orbital%center. 249! Recall that this is done in phiatm() 250 arg = atpoint 251 252! Compute the x, y, and z components of arg 253 x = dot_product( orbital%xaxis, arg ) 254 y = dot_product( orbital%yaxis, arg ) 255 z = dot_product( orbital%zaxis, arg ) 256! Compute the distance between the point and the center of the trial function 257 rr = sphere(x,y,z,1) 258 259! 260! If out of the rcut sphere then vanish 261! 262 if ( rr .gt. orbital%rcut ) then 263 gettrialwavefunction = 0.0_dp 264 return 265 endif 266 267! 268! Decipher arguments: 269! Compute the radial part of the trial function at the given point 270! 271 select case(orbital%r) 272 case(1) 273 gettrialwavefunction = R1(rr) 274 case(2) 275 gettrialwavefunction = R2(rr) 276 case(3) 277 gettrialwavefunction = R3(rr) 278 end select 279! Renormalize the wave function, since we cut it at R_c 280 gettrialwavefunction = gettrialwavefunction / dsqrt(1.0_dp-T) 281 282! 283! Decipher arguments: 284! Compute the angular part of the trial function at the given point 285! 286 if ( rr .eq. 0.0_dp ) then 287! Special treatment of the origin 288 select case(orbital%l) 289 case (0) 290 angular = l0norm 291! Hybrids are combinations of dz2 and s limits 292 case (-1) 293 angular = l0norm * rs2 294 case (-2) 295 angular = l0norm * rs3 296 case (-3) 297 angular = l0norm / 2.0_dp 298 case (-4) 299 if ( orbital%mr .lt. 4 ) then 300! angular = l0norm * rs3 301 else 302 angular = 0.0_dp 303 endif 304 case(-5) 305 angular = l0norm * rs6 306 case default 307 angular = 0.0_dp 308 end select 309 310!! Special treatment of the origin 311! select case(orbital%l) 312! case (0) 313! angular = l0norm 314! case (2) 315! if(orbital%mr.eq.1) then 316!! l=2,mr=1 doesn't vanish! 317!! angular = -l2z2norm ... this is the limit in the z=0 plane 318!! but other limiting value is +2*l2z2norm along the z axis 319!! so there is a cusp and we take it's average: zero 320! else 321! angular = 0.0_dp 322! endif 323!! Hybrids are combinations of dz2 and s limits 324! case (-1) 325! angular = l0norm/sqrt(2.0_dp) 326! case (-2) 327! angular = l0norm/sqrt(3.0_dp) 328! case (-3) 329! angular = l0norm/2.0_dp 330! case (-4) 331! if (orbital%mr.lt.4) then 332! angular = l0norm/sqrt(3.0_dp) 333! else 334! !angular = -l2z2norm/sqrt(2.0_dp) 335! endif 336! case(-5) 337! if (orbital%mr.lt.5) then 338! angular = l0norm/sqrt(6.0_dp)!+l2z2norm/sqrt(12.0_dp) 339! else 340! angular = l0norm/sqrt(6.0_dp)!-l2z2norm/sqrt(3.0_dp) 341! endif 342! case default 343! angular = 0.0_dp 344! end select 345 346 else 347! 348! rr.not equal.0 349! 350 select case(orbital%l) 351 case(0) 352 angular = s(x,y,z) 353 case(1) 354 select case(orbital%mr) 355 case(1) 356 angular = pz(x,y,z) 357 case(2) 358 angular = px(x,y,z) 359 case(3) 360 angular = py(x,y,z) 361 end select 362 case(2) 363 select case(orbital%mr) 364 case(1) 365 angular = dz2(x,y,z) 366 case(2) 367 angular = dxz(x,y,z) 368 case(3) 369 angular = dyz(x,y,z) 370 case(4) 371 angular = dx2y2(x,y,z) 372 case(5) 373 angular = dxy(x,y,z) 374 end select 375 case(3) 376 select case(orbital%mr) 377 case(1) 378 angular = fz3(x,y,z) 379 case(2) 380 angular = fxz2(x,y,z) 381 case(3) 382 angular = fyz2(x,y,z) 383 case(4) 384 angular = fzx2y2(x,y,z) 385 case(5) 386 angular = fxyz(x,y,z) 387 case(6) 388 angular = fxx23y2(x,y,z) 389 case(7) 390 angular = fy3x2y2(x,y,z) 391 end select 392! 393! Hybrids 394! 395 case(-1) 396 select case(orbital%mr) 397 case(1) 398 angular = sp_1(x,y,z) 399 case(2) 400 angular = sp_2(x,y,z) 401 end select 402 case(-2) 403 select case(orbital%mr) 404 case(1) 405 angular = sp2_1(x,y,z) 406 case(2) 407 angular = sp2_2(x,y,z) 408 case(3) 409 angular = sp2_3(x,y,z) 410 end select 411 case(-3) 412 select case(orbital%mr) 413 case(1) 414 angular = sp3_1(x,y,z) 415 case(2) 416 angular = sp3_2(x,y,z) 417 case(3) 418 angular = sp3_3(x,y,z) 419 case(4) 420 angular = sp3_4(x,y,z) 421 end select 422 case(-4) 423 select case(orbital%mr) 424 case(1) 425 angular = sp3d_1(x,y,z) 426 case(2) 427 angular = sp3d_2(x,y,z) 428 case(3) 429 angular = sp3d_3(x,y,z) 430 case(4) 431 angular = sp3d_4(x,y,z) 432 case(5) 433 angular = sp3d_5(x,y,z) 434 end select 435 case(-5) 436 select case(orbital%mr) 437 case(1) 438 angular = sp3d2_1(x,y,z) 439 case(2) 440 angular = sp3d2_2(x,y,z) 441 case(3) 442 angular = sp3d2_3(x,y,z) 443 case(4) 444 angular = sp3d2_4(x,y,z) 445 case(5) 446 angular = sp3d2_5(x,y,z) 447 case(6) 448 angular = sp3d2_6(x,y,z) 449 end select 450 end select 451 endif ! this ends the first if regarding rr = 0 or not 452 453 !print *,3.0_dp/2_dp,1.5_dp,3.0_dp/2.0_dp,3.0_dp/2.0_dp 454 gettrialwavefunction = gettrialwavefunction * angular 455end function gettrialwavefunction 456 457real(dp) function gettrialrcut( ofwhat ) 458 type(trialorbital),intent(in) :: ofwhat 459 gettrialrcut = ofwhat%rCut 460end function 461 462integer function gettriallmax( ofwhat ) 463 type(trialorbital),intent(in) :: ofwhat 464 gettriallmax = ofwhat%lMax 465end function 466 467function gettrialcenter( ofwhat ) 468 real(dp),dimension(3) :: gettrialcenter 469 type(trialorbital),intent(in) :: ofwhat 470 gettrialcenter = ofwhat%center 471end function gettrialcenter 472 473subroutine print_trialorb( what ) 474! 475! This subroutine prints the information about the trial projection function 476! 477 type(trialorbital),intent(in) :: what 478 479 write(*,fmt='(a,3f8.3,a)') "print_trialorb: center = ",what%center," Bohr" 480 write(*,fmt='(a,3f8.3)') "print_trialorb: zaxis = ",what%zaxis 481 write(*,fmt='(a,3f8.3)') "print_trialorb: xaxis = ",what%xaxis 482 write(*,fmt='(a,3f8.3)') "print_trialorb: yaxis = ",what%yaxis 483 write(*,fmt='(a,1f8.3,a)') "print_trialorb: zovera = ",what%zovera," Bohr**-1" 484 write(*,fmt='(a,i5)') "print_trialorb: r = ",what%r 485 write(*,fmt='(a,i5)') "print_trialorb: mr = ",what%mr 486 write(*,fmt='(a,i5)') "print_trialorb: l = ",what%l 487 write(*,'(a,1f8.3,a)') "print_trialorb: rcut = ",what%rcut," Bohr" 488 write(*,fmt='(a,i5)') "print_trialorb: lmax = ",what%lmax 489end subroutine print_trialorb 490 491endmodule trialorbitalclass 492