LCOV - code coverage report
Current view: top level - eigen - tlmplm_cholesky.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 66 89 74.2 %
Date: 2019-09-08 04:53:50 Functions: 2 3 66.7 %

          Line data    Source code
       1             : MODULE m_tlmplm_cholesky
       2             :   use m_judft
       3             :   IMPLICIT NONE
       4             :   !*********************************************************************
       5             :   !     sets up the local Hamiltonian, i.e. the Hamiltonian in the
       6             :   !     l',m',l,m,u- basis which is independent from k!
       7             :   !     shifts this local Hamiltonian to make it positive definite
       8             :   !     and does a cholesky decomposition
       9             :   !*********************************************************************
      10             : CONTAINS
      11         608 :   SUBROUTINE tlmplm_cholesky(sphhar,atoms,noco,enpara,&
      12             :        jspin,jsp,mpi,v,input,td,ud)
      13             :     USE m_tlmplm
      14             :     USE m_types
      15             :     USE m_gaunt, ONLY: gaunt2
      16             :     IMPLICIT NONE
      17             :     TYPE(t_mpi),INTENT(IN)      :: mpi
      18             :     TYPE(t_noco),INTENT(IN)     :: noco
      19             :     TYPE(t_input),INTENT(IN)    :: input
      20             :     TYPE(t_sphhar),INTENT(IN)   :: sphhar
      21             :     TYPE(t_atoms),INTENT(IN)    :: atoms
      22             :     TYPE(t_enpara),INTENT(IN)   :: enpara
      23             :     !     ..
      24             :     !     .. Scalar Arguments ..
      25             :     INTEGER, INTENT (IN) :: jspin,jsp !physical spin&spin index for data
      26             :     !     ..
      27             :     TYPE(t_potden),INTENT(IN)    :: v
      28             :     TYPE(t_tlmplm),INTENT(INOUT) :: td
      29             :     TYPE(t_usdus),INTENT(INOUT)  :: ud
      30             : 
      31             :     !     ..
      32             :     !     .. Local Scalars ..
      33             :     REAL temp
      34             :     INTEGER i,l,lm,lmin,lmin0,lmp,lmplm,lp,info,in
      35             :     INTEGER lpl ,mp,n,m,s,i_u
      36             :     LOGICAL OK
      37             :     !     ..
      38             :     !     .. Local Arrays ..
      39             : 
      40             : 
      41             :     REAL,PARAMETER:: e_shift_min=0.5
      42             :     REAL,PARAMETER:: e_shift_max=65.0
      43             : 
      44             : 
      45             :     !     ..e_shift
      46         608 :     td%e_shift(:,jsp)=0.0
      47         608 :     IF (jsp<3) td%e_shift(:,jsp)=e_shift_min
      48             : 
      49             : 
      50         608 :     td%tdulo(:,:,:,jsp) = CMPLX(0.0,0.0)
      51         608 :     td%tuulo(:,:,:,jsp) = CMPLX(0.0,0.0)
      52         608 :     td%tuloulo(:,:,:,jsp) = CMPLX(0.0,0.0)
      53             : 
      54             : 
      55             : 
      56         608 :     td%h_off=0.0
      57         608 :     !$    call gaunt2(atoms%lmaxd)
      58             :     !$OMP PARALLEL DO DEFAULT(NONE)&
      59             :     !$OMP PRIVATE(temp,i,l,lm,lmin,lmin0,lmp)&
      60             :     !$OMP PRIVATE(lmplm,lp,m,mp,n)&
      61             :     !$OMP PRIVATE(OK,s,in,info)&
      62         608 :     !$OMP SHARED(atoms,jspin,jsp,sphhar,enpara,td,ud,v,mpi,input)
      63             :     DO  n = 1,atoms%ntype
      64        1164 :        CALL tlmplm(n,sphhar,atoms,enpara,jspin,jsp,mpi,v,input,td,ud)
      65        1164 :        OK=.FALSE.
      66             :        cholesky_loop:DO WHILE(.NOT.OK)
      67      202660 :           td%h_loc(:,:,n,jsp)=0.0
      68        1164 :           OK=.TRUE.
      69             :           !
      70             :           !--->    generate the wavefunctions for each l
      71             :           !
      72        1164 :           s=atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+1
      73             :           !Setup local hamiltonian
      74       59608 :           DO lmp=0,atoms%lnonsph(n)*(atoms%lnonsph(n)+2)
      75       58444 :              lp=FLOOR(SQRT(1.0*lmp))
      76       58444 :              mp=lmp-lp*(lp+1)
      77       58444 :              IF (lp>atoms%lmax(n).OR.ABS(mp)>lp) STOP "BUG"
      78             :              !--->             loop over l,m
      79      894368 :              DO l = 0,atoms%lnonsph(n)
      80     3475132 :                 DO m = -l,l
      81     2999308 :                    lm = l* (l+1) + m
      82     2999308 :                    in = td%ind(lmp,lm,n,jsp)
      83     3416688 :                    IF (in/=-9999) THEN
      84      733960 :                       IF (in>=0) THEN
      85      395512 :                          td%h_loc(lm,lmp,n,jsp)    = CONJG(td%tuu(in,n,jsp))
      86      395512 :                          td%h_loc(lm+s,lmp,n,jsp)  = CONJG(td%tud(in,n,jsp))
      87      395512 :                          td%h_loc(lm,lmp+s,n,jsp)  = CONJG(td%tdu(in,n,jsp))
      88      395512 :                          td%h_loc(lm+s,lmp+s,n,jsp)= CONJG(td%tdd(in,n,jsp))
      89             :                       ELSE
      90      338448 :                          td%h_loc(lm,lmp,n,jsp)    = td%tuu(-in,n,jsp)
      91      338448 :                          td%h_loc(lm+s,lmp,n,jsp)  = td%tdu(-in,n,jsp)
      92      338448 :                          td%h_loc(lm,lmp+s,n,jsp)  = td%tud(-in,n,jsp)
      93      338448 :                          td%h_loc(lm+s,lmp+s,n,jsp)= td%tdd(-in,n,jsp)
      94             :                       END IF
      95             :                    END IF
      96             :                 END DO
      97             :              END DO
      98             :           ENDDO
      99             :           !Include contribution from LDA+U
     100        1404 :           DO i_u=1,atoms%n_u
     101         240 :              IF (n.NE.atoms%lda_u(i_u)%atomtype) CYCLE
     102             :              !Found a "U" for this atom type
     103         120 :              l=atoms%lda_u(i_u)%l
     104         120 :              lp=atoms%lda_u(i_u)%l
     105        1764 :              DO m = -l,l
     106         480 :                 lm = l* (l+1) + m
     107        2760 :                 DO mp = -lp,lp
     108        2040 :                    lmp = lp* (lp+1) + mp
     109        2040 :                    td%h_loc(lm,lmp,n,jsp)     =td%h_loc(lm,lmp,n,jsp) + v%mmpMat(m,mp,i_u,jsp)
     110        2520 :                    td%h_loc(lm+s,lmp+s,n,jsp) =td%h_loc(lm+s,lmp+s,n,jsp)+ v%mmpMat(m,mp,i_u,jsp)*ud%ddn(lp,n,jsp)
     111             :                 ENDDO
     112             :              ENDDO
     113             :           END DO
     114             :           
     115        1164 :           IF (jsp<3) THEN
     116             :              !Create Cholesky decomposition of local hamiltonian
     117             :              
     118             :              !--->    Add diagonal terms to make matrix positive definite
     119       17620 :              DO lp = 0,atoms%lnonsph(n)
     120       67836 :                 DO mp = -lp,lp
     121       58444 :                    lmp = lp* (lp+1) + mp
     122       58444 :                    td%h_loc(lmp,lmp,n,jsp)=td%e_shift(n,jsp)+td%h_loc(lmp,lmp,n,jsp)
     123       66672 :                    td%h_loc(lmp+s,lmp+s,n,jsp)=td%e_shift(n,jsp)*ud%ddn(lp,n,jsp)+td%h_loc(lmp+s,lmp+s,n,jsp)
     124             :                 END DO
     125             :              END DO
     126        1164 :              IF (lmp+1.NE.s) CALL judft_error("BUG in tlmpln_cholesky")
     127             :              !Perform cholesky decomposition
     128        1164 :              info=0
     129        1164 :              CALL zpotrf("L",2*s,td%h_loc(:,:,n,jsp),SIZE(td%h_loc,1),info)
     130             : 
     131             :              !Upper part to zero
     132      118052 :              DO l=0,2*s-1
     133    11998396 :                 DO lp=0,l-1
     134     6057060 :                    td%h_loc(lp,l,n,jsp)=0.0
     135             :                 ENDDO
     136             :              ENDDO
     137             : 
     138        1164 :              IF (info.NE.0) THEN
     139           0 :                 td%e_shift(n,jsp)=td%e_shift(n,jsp)*2.0
     140           0 :                 PRINT *,"Potential shift for atom type ",n," is too small. Increasing the value to:",td%e_shift(n,jsp)
     141           0 :                 IF (td%e_shift(n,jsp)>e_shift_max) THEN
     142           0 :                    CALL judft_error("Potential shift at maximum")
     143             :                 ENDIF
     144             :                 OK=.FALSE.
     145             :              ENDIF
     146             :           ENDIF
     147             :        ENDDO cholesky_loop
     148             :        !Now add diagonal contribution to matrices
     149        2380 :        IF (jsp<3) THEN
     150       11808 :           DO l = 0,atoms%lmax(n)
     151       99288 :              DO  m = -l,l
     152       98124 :                 lm = l* (l+1) + m
     153       98124 :                 lmplm = (lm* (lm+3))/2
     154       98124 :                 td%tuu(lmplm,n,jsp)=td%tuu(lmplm,n,jsp) + enpara%el0(l,n,jsp)
     155       98124 :                 td%tdd(lmplm,n,jsp)=td%tdd(lmplm,n,jsp) + enpara%el0(l,n,jsp)*ud%ddn(l,n,jsp)
     156       98124 :                 td%tud(lmplm,n,jsp)=td%tud(lmplm,n,jsp) + 0.5
     157      108768 :                 td%tdu(lmplm,n,jsp)=td%tdu(lmplm,n,jsp) + 0.5
     158             :              ENDDO
     159             :           ENDDO
     160             :        ENDIF
     161             :     ENDDO
     162             :     !$OMP END PARALLEL DO
     163         608 :     IF (noco%l_constr) CALL tlmplm_constrained(atoms,v,enpara,input,ud,noco,td)
     164             : 
     165             : 
     166             : 
     167         608 :   END SUBROUTINE tlmplm_cholesky
     168             : 
     169             : 
     170             :    
     171             : 
     172             : 
     173           0 :   SUBROUTINE tlmplm_constrained(atoms,v,enpara,input,ud,noco,td)
     174             :     USE m_radovlp
     175             :     USE m_types
     176             :     IMPLICIT NONE
     177             :     TYPE(t_input),INTENT(IN)    :: input
     178             :     TYPE(t_atoms),INTENT(IN)    :: atoms
     179             :     TYPE(t_enpara),INTENT(IN)   :: enpara
     180             :     TYPE(t_potden),INTENT(IN)   :: v
     181             :     TYPE(t_tlmplm),INTENT(INOUT):: td
     182             :     TYPE(t_usdus),INTENT(INOUT) :: ud
     183             :     TYPE(t_noco),INTENT(IN)     :: noco
     184             :     
     185           0 :     REAL, ALLOCATABLE :: uun21(:,:),udn21(:,:),dun21(:,:),ddn21(:,:)
     186             :     COMPLEX :: c
     187             :     INTEGER :: n,l,s  
     188             :       
     189             :     
     190             :     ALLOCATE(uun21(0:atoms%lmaxd,atoms%ntype),udn21(0:atoms%lmaxd,atoms%ntype),&
     191           0 :          dun21(0:atoms%lmaxd,atoms%ntype),ddn21(0:atoms%lmaxd,atoms%ntype) )
     192           0 :     CALL rad_ovlp(atoms,ud,input,v%mt,enpara%el0, uun21,udn21,dun21,ddn21)
     193             :     
     194           0 :     DO  n = 1,atoms%ntype
     195             :        !If we do  a constraint calculation, we have to calculate the
     196             :        !local spin off-diagonal contributions
     197           0 :        s=atoms%lnonsph(n)+1
     198             :        !first ispin=2,jspin=1 case
     199           0 :        DO l=0,atoms%lnonsph(n)
     200           0 :           c=(-0.5)*CMPLX(noco%b_con(1,n),noco%b_con(2,n))
     201           0 :           td%h_off(l  ,l  ,n,1)     =td%h_off(l  ,l  ,n,1) + uun21(l,n)*c
     202           0 :           td%h_loc(l  ,l+s,n,1)     =td%h_off(l  ,l+s,n,1) + udn21(l,n)*c
     203           0 :           td%h_loc(l+s,l  ,n,1)     =td%h_off(l+s,l  ,n,1) + dun21(l,n)*c
     204           0 :           td%h_loc(l+s,l+s,n,1)     =td%h_off(l+s,l+s,n,1) + ddn21(l,n)*c
     205             :        ENDDO
     206             :        
     207             :        
     208             :        !then ispin=2,jspin=1 case
     209           0 :        DO l=0,atoms%lnonsph(n)
     210           0 :           c=(-0.5)*CMPLX(noco%b_con(1,n),-noco%b_con(2,n))
     211           0 :           td%h_off(l  ,l  ,n,2)     =td%h_off(l  ,l  ,n,2) + uun21(l,n)*c
     212           0 :           td%h_loc(l  ,l+s,n,2)     =td%h_off(l  ,l+s,n,2) + udn21(l,n)*c
     213           0 :           td%h_loc(l+s,l  ,n,2)     =td%h_off(l+s,l  ,n,2) + dun21(l,n)*c
     214           0 :           td%h_loc(l+s,l+s,n,2)     =td%h_off(l+s,l+s,n,2) + ddn21(l,n)*c
     215             :        ENDDO
     216             :     END DO
     217           0 :   END SUBROUTINE tlmplm_constrained
     218             :   
     219             :   END MODULE m_tlmplm_cholesky
     220             :   

Generated by: LCOV version 1.13