LCOV - code coverage report
Current view: top level - orbdep - nabla.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 79 81 97.5 %
Date: 2019-09-08 04:53:50 Functions: 2 2 100.0 %

          Line data    Source code
       1             :       MODULE m_nabla
       2             :       use m_juDFT
       3             :         
       4             :       CONTAINS
       5             : 
       6          48 :       SUBROUTINE nabla(
       7             :      >                 ispecies,number_of_j1,grid_size,delta_x,
       8          48 :      >                 nstd,ntypd,j1,l1,lmax,ms,ri,psi,phi,dphi,
       9          48 :      <                 psi_phi)
      10             : !----------------------------------------------------------------
      11             : !
      12             : !     ispecies     ... number of the atom (itype)
      13             : !     number_of_j1 ... number of core wavefunction
      14             : !     grid_size    ... rumber of radial gridpoints (jri)
      15             : !     delta_x      ... logarithmic increment of grid (dx)
      16             : !     ri(grid_size)... radial mesh
      17             : !     nstd     ... number of corelevels (dimension)
      18             : !     ntypd    ... number of atoms (dimension)
      19             : !     j1,l1    ... quantum numbers of the core wavefunction
      20             : !     ms       ... -1/2 or +1/2 for spin 1 or 2
      21             : !     lmax     ... parameter, 3 (s,p,d,f)
      22             : !     psi(r)   ... core wavefunction
      23             : !     phi(r,l) ... valence wavefunction
      24             : !    dphi(r,l) ... radial derivative of valence wavefunction
      25             : !     
      26             : !----------------------------------------------------------------
      27             :        USE m_clebsch
      28             :        USE m_intgr, ONLY : intgr3
      29             :        IMPLICIT NONE  
      30             : 
      31             :        INTEGER, INTENT(IN) :: ispecies, number_of_j1, grid_size
      32             :        INTEGER, INTENT(IN) :: l1, lmax, nstd, ntypd
      33             :        REAL,    INTENT(IN) :: delta_x, j1, ms
      34             :        REAL,    INTENT(IN) :: psi(grid_size), ri(grid_size)
      35             :        REAL,    INTENT(IN) :: phi(grid_size,0:lmax)
      36             :        REAL,    INTENT(IN) ::dphi(grid_size,0:lmax)
      37             :        COMPLEX, INTENT(OUT):: psi_phi(nstd,(lmax+1)**2,3*ntypd)
      38             : 
      39             :        INTEGER :: m1, l2, m2, index, alloc_error, lmn1, lmn2
      40             :        REAL  :: result, result1, total_result, spin
      41             :        REAL  :: mu, cnst_one_over_sqrt_two, cnst_zero
      42             :        COMPLEX :: cnst_i
      43             :        REAL, DIMENSION(:), POINTER :: f
      44          48 :        spin = 0.50
      45          48 :        cnst_one_over_sqrt_two = 1.0/sqrt(2.0)
      46          48 :        cnst_i = cmplx(0.0,1.0)
      47          48 :        cnst_zero = 0.0
      48             : 
      49             :        NULLIFY(f)
      50             : 
      51             :        IF ( ASSOCIATED(f) ) THEN
      52             :         WRITE(6,*)'nabla: f association status:',ASSOCIATED(f)
      53             :         STOP
      54             :        ENDIF
      55             : 
      56          48 :        ALLOCATE ( f(grid_size),STAT=alloc_error )
      57          48 :        IF (alloc_error /= 0)  CALL juDFT_error("Couldn't allocate f",
      58           0 :      +      calledby ="nabla")
      59             :       
      60          48 :        lmn1 = 2 * (number_of_j1 - 1)  * l1
      61          48 :        mu = -j1
      62             : 
      63         304 :        DO WHILE (mu <= j1)
      64         128 :         lmn1 = lmn1 + 1
      65         128 :         m1 = INT(mu - ms) 
      66         128 :         lmn2 = 0
      67         640 :         DO l2 = 0, lmax
      68        2688 :          DO m2 = -l2, l2
      69        2048 :             lmn2 = lmn2 + 1
      70        2560 :             IF(l1 == l2 + 1)THEN  
      71          96 :              total_result = 0.00
      72          96 :              result  = 0.00
      73          96 :              result1 = 0.00
      74             : !     
      75             : ! (l+1)/srqt[(2l+1)(2l+3)] < phi_core | (d phi_valence / dr ) >
      76             : !
      77          96 :              f(:) = psi(:) * dphi(:,l2) ! assumed to be already multiplied with * ri(:) * ri(:)
      78             :       
      79          96 :              CALL intgr3(f,ri,delta_x,grid_size,result)
      80             :       
      81             :              result = result * (l2 + 1.00) /
      82          96 :      +                         sqrt((2.00*l2 +1.00)*(2.*l2+3.00))
      83             :       
      84             : !
      85             : !  - l(l+1)/srqt[(2l+1)(2l+3)] < phi_core | (1/r) phi_valence >
      86             : !
      87          96 :              f(:) = psi(:) * phi(:,l2) ! assumed to be already multiplied with 1G* ri(:)
      88          96 :              CALL intgr3(f,ri,delta_x,grid_size,result1)
      89             :              result1 = - result1 * l2 * (l2 + 1.0) /              
      90          96 :      +                   sqrt((2.0 * l2 + 1.00) * (2. * l2 + 3.00)) 
      91             : !
      92             : ! Sum up and decorate with Clebsch-Gordon coefficients
      93             : !
      94          96 :              result  = result + result1
      95          96 :              total_result = result*clebsch(real(l1),spin,mu-ms,ms,j1,mu)
      96             :       
      97          96 :              index = (ispecies - 1) * 3 + 1
      98             :              psi_phi(lmn1,lmn2,index)= cgc(l2,1,l1,m2,1,m1) *       ! left polarization
      99          96 :      +                                 total_result / cgc(l2,1,l1,0,0,0)
     100             :       
     101          96 :              index = index + 1
     102             :              psi_phi(lmn1,lmn2,index)= cgc(l2,1,l1,m2,-1,m1) *      ! right polarization
     103          96 :      +                                 total_result / cgc(l2,1,l1,0,0,0)
     104             :       
     105          96 :              index = index + 1
     106             :              psi_phi(lmn1,lmn2,index)= cgc(l2,1,l1,m2,0,m1) *        ! z-polarization
     107          96 :      +                                 total_result / cgc(l2,1,l1,0,0,0)
     108             :       
     109        1952 :             ELSEIF(l1== l2-1)THEN  
     110             :       
     111         576 :               result  =  cnst_zero 
     112         576 :               result1 =  cnst_zero
     113             : !
     114             : ! l/srqt[(2l-1)(2l+1)] < phi_core | (d phi_valence / dr ) >
     115             : !
     116         576 :               f(:) = psi(:)* dphi(:,l2)   * ri(:) * ri(:)
     117         576 :               CALL intgr3(f,ri,delta_x,grid_size,result)
     118         576 :               result = result * l2 / sqrt((2.0*l2 - 1.0)*(2.0*l2 + 1.0)) 
     119             : !
     120             : !   l(l+1)/srqt[(2l-1)(2l+1)] < phi_core | (1/r) phi_valence >
     121             : !
     122         576 :               f(:) = psi(:)* phi(:,l2) * ri(:)
     123         576 :               CALL intgr3(f,ri,delta_x,grid_size,result1)
     124             :               result1 = result1 * l2 * (l2 + 1.0) /
     125         576 :      +                  sqrt((2.00 * l2 - 1.00) * (2.00 * l2 + 1.00)) 
     126             : !
     127             : ! Sum up and decorate with Clebsch-Gordon coefficients
     128             : !
     129         576 :               result = result + result1
     130         576 :               total_result= result*clebsch(real(l1),spin,mu-ms,ms,j1,mu)
     131             :       
     132             :       
     133         576 :               index = (ispecies - 1) * 3 + 1
     134             :               psi_phi(lmn1,lmn2,index) = cgc(l2,1,l1,m2,1,m1) *          ! left polarization
     135         576 :      +                                 total_result / cgc(l2,1,l1,0,0,0)
     136         576 :               index = index + 1
     137             :               psi_phi(lmn1,lmn2,index) = cgc(l2,1,l1,m2,-1,m1) *         ! right polarization
     138         576 :      +                                 total_result / cgc(l2,1,l1,0,0,0)
     139         576 :               index = index + 1
     140             :               psi_phi(lmn1,lmn2,index) = cgc(l2,1,l1,m2,0,m1) *          ! z-polarization
     141         576 :      +                                 total_result / cgc(l2,1,l1,0,0,0)
     142             :             ENDIF         
     143             :           ENDDO
     144             :         ENDDO
     145         176 :         mu = mu + 1.00
     146             :        ENDDO
     147          48 :        DEALLOCATE(f,STAT=alloc_error)
     148             :        IF (alloc_error /= 0)  CALL juDFT_error("Couldn't deallocate f"
     149           0 :      +      ,calledby ="nabla")
     150          48 :       END SUBROUTINE nabla
     151             : 
     152        4032 :       FUNCTION cgc(l1,l2,l3,m1,m2,m3)
     153             : 
     154             :       IMPLICIT NONE  
     155             :       INTEGER :: l1, l2, l3, m1, m2, m3
     156             :       REAL  :: two_l1p1, two_l1p2, l1pm3, l1pm3p1, l1mm3p1, l1mm3, cgc 
     157             : 
     158        4032 :       IF (m3 /= m1 + m2) THEN
     159             :        cgc = 0.0
     160             :        RETURN
     161             :       END IF 
     162             : !     gb  m3 = m1 + m2
     163        2448 :       two_l1p1 = 2 * l1 + 1
     164        2448 :       two_l1p2 = 2 * l1 + 2
     165        2448 :       l1pm3 = l1 + m3
     166        2448 :       l1pm3p1 = l1 + m3 + 1
     167        2448 :       l1mm3p1 = l1 - m3 + 1
     168        2448 :       l1mm3 = l1 - m3 
     169        2448 :       cgc = 0.0 
     170        2448 :       IF (l3 == l1 + 1) THEN
     171         368 :           IF (m2 == 1) then
     172          24 :            cgc = sqrt( (l1pm3 * l1pm3p1) / (two_l1p1 * two_l1p2))   
     173         344 :           ELSEIF (m2 == 0) THEN
     174         320 :            cgc = sqrt( (l1mm3p1 * l1pm3p1) / (two_l1p1 * (l1 + 1)))
     175          24 :           ELSEIF (m2 == -1) THEN
     176          24 :            cgc = sqrt( (l1mm3 * l1mm3p1) / (two_l1p1 * two_l1p2))   
     177             :           END IF
     178        2080 :       ELSE IF(l3 == l1 -1) THEN
     179        2080 :           IF (m2 == 1) then
     180         112 :            cgc = sqrt( (l1mm3 * l1mm3p1) / (2.0 * l1 * two_l1p1))   
     181        1968 :           ELSEIF (m2 == 0) THEN
     182        1856 :            cgc = -sqrt( (l1mm3 * l1pm3) / (l1 * (two_l1p1)))   
     183         112 :           ELSEIF (m2 == -1) THEN
     184         112 :            cgc = sqrt( (l1pm3p1 * l1pm3) / (2.0 * l1 * two_l1p1))   
     185             :           END IF
     186             :       END IF
     187             :       END FUNCTION cgc
     188             : 
     189             : 
     190             :       END MODULE m_nabla

Generated by: LCOV version 1.13