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