Line data    Source code

       1             : MODULE m_dfpt_dynmat_sym
       2             :    USE m_juDFT
       3             :    USE m_types
       4             :    USE m_constants
       5             : 
       6             :    IMPLICIT NONE
       7             : CONTAINS
       8           0 :    SUBROUTINE ft_dyn(atoms, qpts, sym, amat, dyn_mat_q, dyn_mat_r, dyn_mat_q_full)
       9             :       !! Transforms the dynamical matrices for a set of q vectors in the
      10             :       !! irreducible Brillouin zone onto the full set of q vector in the BZ
      11             :       !! and subsequently transforms it to real space (lattice vector grid), to
      12             :       !! calculate the mass-normalized Force Constant Matrix.
      13             :       type(t_atoms), INTENT(IN)   :: atoms
      14             :       type(t_kpts),  INTENT(IN)   :: qpts
      15             :       type(t_sym),   INTENT(IN)   :: sym
      16             :       REAL,          INTENT(IN)   :: amat(3,3)
      17             :       COMPLEX,       INTENT(IN)  :: dyn_mat_q(:,:,:) ! (nqpt,dyn_dim,dyn_dim)
      18             :       COMPLEX,       INTENT(OUT) :: dyn_mat_r(:,:,:) ! (nqptf,dyn_dim,dyn_dim)
      19             : 
      20             :       COMPLEX, ALLOCATABLE, INTENT(OUT) :: dyn_mat_q_full(:,:,:)
      21             : 
      22             :       INTEGER :: iq, dyn_dim, iqfull
      23             :       INTEGER :: isym, r_lim(2,3)
      24             :       INTEGER :: j1, j2, j3
      25             :       REAL    :: q_full(3), q_full_BZ(3)
      26             : 
      27           0 :       INTEGER, ALLOCATABLE :: qvec_to_index(:,:,:)
      28             :       COMPLEX, ALLOCATABLE :: dyn_mat_qsym(:,:)
      29             : 
      30           0 :       dyn_dim = 3*atoms%nat
      31             : 
      32           0 :       allocate(qvec_to_index(0:qpts%nkpt3(1)-1,0:qpts%nkpt3(2)-1,0:qpts%nkpt3(3)-1))
      33           0 :       allocate(dyn_mat_qsym(dyn_dim,dyn_dim))
      34           0 :       allocate(dyn_mat_q_full(qpts%nkptf,dyn_dim,dyn_dim))
      35             : 
      36           0 :       qvec_to_index = 0
      37             : 
      38             :       ! Create an array that maps the q coordinates to the index of their q vector
      39             :       ! in the full BZ
      40           0 :       DO j1 = 0, qpts%nkpt3(1)-1
      41           0 :          q_full(1) = j1*1.0/qpts%nkpt3(1)
      42           0 :          DO j2 = 0, qpts%nkpt3(2)-1
      43           0 :             q_full(2) = j2*1.0/qpts%nkpt3(2)
      44           0 :             DO j3 = 0, qpts%nkpt3(3)-1
      45           0 :                q_full(3) = j3*1.0/qpts%nkpt3(3)
      46           0 :                DO iq = 1, qpts%nkptf
      47           0 :                   IF (norm2(q_full-qpts%bkf(:,iq))<1e-8) qvec_to_index(j1, j2, j3) = iq
      48             :                END DO
      49             :             END DO
      50             :          END DO
      51             :       END DO
      52             : 
      53           0 :       r_lim(2,:) = qpts%nkpt3(:)/2
      54           0 :       r_lim(1,:) = r_lim(2,:) - qpts%nkpt3(:) + 1
      55             : 
      56           0 :       dyn_mat_r(:,:,:) = cmplx(0.0,0.0)
      57             : 
      58           0 :       DO j1 = 0, qpts%nkpt3(1)-1
      59           0 :          q_full(1) = j1*1.0/qpts%nkpt3(1)
      60           0 :          DO j2 = 0, qpts%nkpt3(2)-1
      61           0 :             q_full(2) = j2*1.0/qpts%nkpt3(2)
      62           0 :             DO j3 = 0, qpts%nkpt3(3)-1
      63           0 :                q_full(3) = j3*1.0/qpts%nkpt3(3)
      64             :                ! Get q vector index and that of its representative in the irreducible wedge
      65           0 :                iqfull = qvec_to_index(j1, j2, j3)
      66           0 :                iq = qpts%bkp(iqfull)
      67             :                ! Fold vector back to 1st BZ if necessary
      68           0 :                q_full_BZ(:)=q_full(:)
      69           0 :                q_full_BZ = qpts%to_first_bz(q_full)
      70             : 
      71           0 :                isym = sym%invtab(qpts%bksym(iqfull))
      72           0 :                dyn_mat_qsym(:,:) = cmplx(0.0,0.0)
      73           0 :                CALL rotate_dynmat(atoms,sym,isym,amat,qpts%bk(:,iq),dyn_mat_q(iq,:,:),dyn_mat_qsym)
      74             :                dyn_mat_qsym(:,:) = dyn_mat_qsym(:,:)
      75           0 :                dyn_mat_q_full(iqfull,:,:) = dyn_mat_qsym
      76             : 
      77             :                ! Perform the actual FT onto the lattice vector grid
      78           0 :                CALL ft_dyn_direct(r_lim,1,q_full,dyn_mat_qsym,dyn_mat_r)
      79             :             END DO
      80             :          END DO
      81             :       END DO
      82           0 :       dyn_mat_r(:,:,:)=dyn_mat_r(:,:,:)/qpts%nkptf
      83           0 :    END SUBROUTINE
      84             : 
      85           0 :    SUBROUTINE ft_dyn_direct(ft_lim,isn,bqpt,dyn_mat_q,dyn_mat_r)
      86             :       INTEGER, INTENT(IN) :: ft_lim(2,3), isn
      87             :       REAL,    INTENT(IN) :: bqpt(3)
      88             : 
      89             :       COMPLEX :: dyn_mat_q(:,:)
      90             :       COMPLEX :: dyn_mat_r(:,:,:)
      91             : 
      92             :       INTEGER :: iGrid, ix, iy, iz, iout
      93             :       REAL    :: phas
      94             :       COMPLEX :: phase_fac
      95           0 :       iGrid=0
      96           0 :       DO iz=ft_lim(1,3),ft_lim(2,3)
      97           0 :          DO iy=ft_lim(1,2),ft_lim(2,2)
      98           0 :             DO ix=ft_lim(1,1),ft_lim(2,1)
      99           0 :                iGrid = iGrid+1
     100           0 :                phas=isn*tpi_const*(bqpt(1)*ix+bqpt(2)*iy+bqpt(3)*iz)
     101           0 :                phase_fac=cmplx(cos(phas),sin(phas))
     102           0 :                IF (isn==1) THEN
     103           0 :                   dyn_mat_r(iGrid,:,:) = dyn_mat_r(iGrid,:,:) + phase_fac*dyn_mat_q(:,:)
     104           0 :                ELSE IF (isn==-1) THEN
     105           0 :                   dyn_mat_q(:,:)    = dyn_mat_q(:,:)    + phase_fac*dyn_mat_r(iGrid,:,:)
     106             :                END IF
     107             :             END DO
     108             :          END DO
     109             :       END DO
     110           0 :    END SUBROUTINE
     111             : 
     112           0 :    SUBROUTINE rotate_dynmat(atoms,sym,isym,amat,bqpt,dyn,dyn_mat_qsym)
     113             :       !! Applies a symmetry operation to the dynamical matrix of an IBZ q vector
     114             :       !! to find the matrix of its mapped q vector in the full BZ. This is done
     115             :       !! by using the symmetry relation of the FCM when a symmetry operation
     116             :       !! \(\underline{B}\) maps atoms \(\beta',\alpha'\) onto \(\beta,\alpha\)
     117             :       !! in Cartesian coordinates 
     118             :       !! (\(\underline{B}\boldsymbol{\tau}_{\beta'}=\boldsymbol{\tau}_{\beta}\))
     119             :       !! $$\underline{\Phi}_{\alpha'+\boldsymbol{R},\beta'}=\underline{B}\underline{\Phi}_{\alpha+\underline{B}\boldsymbol{R},\beta}\underline{B}^{-1}$$
     120             :       !! Resulting (with the definition of the DM as the mass-scaled Fourier Transform
     121             :       !! of the FCM) in the corresponding relation:
     122             :       !! $$\underline{D}_{\alpha',\beta'}(\boldsymbol{q})=p(\alpha,\beta)*\underline{B}\underline{D}_{\alpha,\beta}(\boldsymbol{q}_{\mathrm{rep}})\underline{B}^{-1}$$
     123             :       !! with a phase factor
     124             :       !! $$f(\alpha,\beta)=exp(ix),x=\boldsymbol{q}\cdot(\boldsymbol{\tau}_{\beta'}-\boldsymbol{\tau}_{\alpha'})-\boldsymbol{q}_{\mathrm{red}}\cdot(\boldsymbol{\tau}_{\beta}-\boldsymbol{\tau}_{\alpha})$$
     125             :       !! which can be written as a product
     126             :       !! $$f(\alpha,\beta) = f(\alpha)f^{*}(\beta), f(\alpha)=exp(i(\boldsymbol{q}_{\mathrm{red}}\cdot\boldsymbol{\tau}_{\alpha}-\boldsymbol{q}\cdot\boldsymbol{\tau}_{\alpha'})).$$
     127             :       !! The real space rotation is related to the rotation matrix of the symmetry operation by the Bravais matrix of the system
     128             :       !! $$\underline{B}=\underline{A}\underline{S}^{-1}\underline{A}^{-1}$$
     129             :       USE m_inv3
     130             :       type(t_atoms), INTENT(IN)    :: atoms
     131             :       type(t_sym),   INTENT(IN)    :: sym
     132             :       INTEGER,       INTENT(IN)    :: isym
     133             :       REAL,          INTENT(IN)    :: amat(3,3)
     134             :       REAL,          INTENT(IN)    :: bqpt(3)
     135             :       COMPLEX,       INTENT(IN)    :: dyn(:,:)
     136             :       COMPLEX,       INTENT(INOUT) :: dyn_mat_qsym(:,:)
     137             : 
     138             :       INTEGER :: iAtom, jAtom
     139             :       INTEGER :: iAlpha, iBeta, iAlpha_map, iBeta_map
     140             :       REAL    :: mrot(3,3), invmrot(3,3), invamat(3,3), q_full(3), phas, det
     141             :       COMPLEX :: brot(3,3), temp_mat_1(3,3), temp_mat_2(3,3)
     142             :       COMPLEX :: phase_fac
     143             : 
     144           0 :       INTEGER :: map(atoms%nat)
     145           0 :       COMPLEX :: phase_map(atoms%nat)
     146             : 
     147           0 :       mrot = sym%mrot(:,:,isym)
     148           0 :       invmrot = sym%mrot(:,:,sym%invtab(isym))
     149             : 
     150           0 :       CALL inv3(amat,invamat,det)
     151           0 :       temp_mat_1 = MATMUL(invmrot,invamat)
     152           0 :       brot = MATMUL(amat,temp_mat_1)
     153             :       ! Find the q vector in the full BZ
     154           0 :       q_full = MATMUL(mrot,bqpt)
     155             : 
     156             :       ! Calculate the array of phases
     157           0 :       DO iAtom = 1, atoms%nat
     158           0 :          jAtom = sym%mapped_atom(isym,iAtom)
     159           0 :          map(iAtom)=jAtom
     160             :          ! Calculate the phase factor f(alpha)
     161           0 :          phas=tpi_const*(dot_product(bqpt(:),atoms%taual(:,iAtom))-dot_product(q_full(:),atoms%taual(:,jAtom)))
     162           0 :          phase_fac=cmplx(cos(phas),sin(phas))
     163           0 :          phase_map(iAtom)=phase_fac
     164             :       END DO
     165             :       ! Transform the dynamical matrix from the representative atom and q vector to the unfolded ones
     166           0 :       DO iAtom=1, atoms%nat
     167           0 :          iAlpha = 3*(iAtom-1)
     168           0 :          iAlpha_map = 3*(map(iAtom)-1)
     169           0 :          DO jAtom=1, atoms%nat
     170           0 :             iBeta = 3*(jAtom-1)
     171           0 :             iBeta_map = 3*(map(jAtom)-1)
     172           0 :             temp_mat_1 = dyn(iAlpha+1:iAlpha+3,iBeta+1:iBeta+3)
     173           0 :             temp_mat_2 = MATMUL(brot,temp_mat_1)
     174           0 :             temp_mat_1 = MATMUL(temp_mat_2,TRANSPOSE(brot))
     175           0 :             phase_fac=phase_map(iAtom)*conjg(phase_map(jAtom))
     176             : 
     177             :             dyn_mat_qsym(iAlpha_map+1:iAlpha_map+3,iBeta_map+1:iBeta_map+3) &
     178           0 :           = dyn_mat_qsym(iAlpha_map+1:iAlpha_map+3,iBeta_map+1:iBeta_map+3) + phase_fac*temp_mat_1
     179             :          END DO
     180             :       END DO
     181           0 :    END SUBROUTINE
     182             : 
     183           0 :    SUBROUTINE ift_dyn(atoms,qpts,sym,amat,bqpt,dyn_mat_r,dyn_mat_q)
     184             :       !! Transforms the dynamical matrix on a real space lattice vector grid
     185             :       !! (--> mass-normalized FCM) back onto a specific q vector provided as
     186             :       !! input (bqpt) by the inverse Fourier Transformation as compared to
     187             :       !! SUBROUTINE ft_dyn.
     188             :       type(t_atoms), INTENT(IN) :: atoms
     189             :       type(t_kpts),  INTENT(IN) :: qpts
     190             :       type(t_sym),   INTENT(IN) :: sym
     191             :       REAL,          INTENT(IN) :: amat(3,3)
     192             :       REAL,          INTENT(IN) :: bqpt(3)
     193             :       COMPLEX,       INTENT(IN) :: dyn_mat_r(:,:,:)
     194             : 
     195             :       COMPLEX, ALLOCATABLE, INTENT(OUT) :: dyn_mat_q(:,:)
     196             : 
     197             :       INTEGER :: dyn_dim, q_lim(2,3)
     198             : 
     199           0 :       dyn_dim = 3*atoms%nat
     200             : 
     201           0 :       allocate(dyn_mat_q(dyn_dim,dyn_dim))
     202           0 :       dyn_mat_q(:,:) = cmplx(0.0,0.0)
     203             : 
     204           0 :       q_lim(2,:) = qpts%nkpt3(:)/2
     205           0 :       q_lim(1,:) = q_lim(2,:) - qpts%nkpt3(:) + 1
     206             : 
     207           0 :       CALL ft_dyn_direct(q_lim,-1,bqpt,dyn_mat_q,dyn_mat_r)
     208           0 :    END SUBROUTINE
     209             : 
     210           0 :    SUBROUTINE make_sym_list(sym, bqpt, sym_count, sym_list)
     211             :       TYPE(t_sym), INTENT(IN)    :: sym
     212             :       REAL,        INTENT(IN)    :: bqpt(3)
     213             :       INTEGER,     INTENT(OUT)   :: sym_count
     214             :       INTEGER,     INTENT(INOUT) :: sym_list(:)
     215             :       
     216             :       INTEGER :: iSym
     217             :       
     218           0 :       sym_count = 0
     219           0 :       sym_list = 0
     220           0 :       DO iSym = 1, sym%nop
     221           0 :          IF (norm2(bqpt-MATMUL(bqpt,sym%mrot(:,:,iSym)))<1e-8) THEN
     222           0 :             sym_count = sym_count + 1
     223           0 :             sym_list(sym_count) = iSym
     224             :          END IF
     225             :       END DO
     226           0 :    END SUBROUTINE
     227             : 
     228           0 :    SUBROUTINE make_sym_dynvec(atoms, sym, amat, bqpt, iDtype, iDir, sym_count, sym_list, dynvec, sym_dynvec)
     229             :       USE m_inv3
     230             : 
     231             :       TYPE(t_atoms), INTENT(IN)    :: atoms
     232             :       TYPE(t_sym),   INTENT(IN)    :: sym
     233             :       REAL,          INTENT(IN)    :: amat(3,3), bqpt(3)
     234             :       INTEGER,       INTENT(IN)    :: iDtype, iDir, sym_count
     235             :       INTEGER,       INTENT(IN)    :: sym_list(sym_count)
     236             :       COMPLEX,       INTENT(IN)    :: dynvec(:)
     237             :       COMPLEX,       INTENT(INOUT) :: sym_dynvec(:,:,:)
     238             :      
     239             :       INTEGER :: iSym, iAtom, iRow, iCol
     240             :       REAL    :: phas, det, mrot(3,3), invmrot(3,3), invamat(3,3)
     241             :       COMPLEX :: brot(3,3), temp_mat_1(3,3), temp_mat_2(3,3)
     242             :       COMPLEX :: phase_fac, rotvec(3)
     243             : 
     244           0 :       iRow = 3 * (iDtype-1) + iDir
     245             : 
     246           0 :       DO iSym = 1, sym_count
     247           0 :          mrot = sym%mrot(:,:,sym_list(iSym))
     248           0 :          invmrot = sym%mrot(:,:,sym%invtab(sym_list(iSym)))
     249           0 :          CALL inv3(amat,invamat,det)
     250           0 :          temp_mat_1 = MATMUL(invmrot,invamat)
     251           0 :          brot = MATMUL(amat,temp_mat_1)
     252           0 :          DO iAtom = 1, atoms%nat
     253           0 :             iCol = 3 * (iAtom-1)
     254             : 
     255           0 :             phas = -tpi_const*(dot_product(bqpt(:),atoms%taual(:,iAtom)-atoms%taual(:,iDtype)))
     256           0 :             phase_fac = cmplx(cos(phas),sin(phas))
     257             : 
     258           0 :             rotvec = MATMUL(brot,dynvec(iCol+1:iCol+3))
     259           0 :             sym_dynvec(iCol+1:iCol+3,iRow,iSym) = phase_fac * rotvec
     260             :          END DO
     261             :       END DO 
     262           0 :    END SUBROUTINE
     263             : 
     264           0 :    SUBROUTINE cheat_dynmat(atoms, sym, amat, bqpt, iBetaPr, jDirPr, sym_count, sym_list, sym_dynvec, dynmat, sym_dynmat, l_cheated)
     265             :       USE m_inv3
     266             : 
     267             :       TYPE(t_atoms), INTENT(IN)    :: atoms
     268             :       TYPE(t_sym),   INTENT(IN)    :: sym
     269             :       REAL,          INTENT(IN)    :: amat(3,3), bqpt(3)
     270             :       INTEGER,       INTENT(IN)    :: iBetaPr, jDirPr, sym_count
     271             :       INTEGER,       INTENT(IN)    :: sym_list(sym_count)
     272             :       COMPLEX,       INTENT(IN)    :: sym_dynvec(:,:,:), dynmat(:,:)
     273             :       COMPLEX,       INTENT(INOUT) :: sym_dynmat(:,:)
     274             :       LOGICAL,       INTENT(OUT)   :: l_cheated
     275             : 
     276             :       INTEGER :: iSym, iDatom, iAtom, iRowPr, iColPr, iRow, iCol, iDir, jDir, iDirPr, kDir, iBeta, iAlpha, iAlphaPr, symsumcount, aAlpha(atoms%nat), iDone, iNeed, kNeed
     277             :       REAL    :: phas, det, mrot(3,3), invmrot(3,3), invamat(3,3)
     278             :       COMPLEX :: brot(3,3), temp_mat_1(3,3), temp_mat_2(3,3), symsum(3,3), symvec(3)
     279           0 :       COMPLEX :: phase_fac, rotvec(3), rhs_vecs(3,3 * (iBetaPr-1) + jDirPr - 1), rhs_vec_full(3)
     280             :       LOGICAL :: l_mapped, l_done(3), l_need(3)
     281             : 
     282           0 :       l_cheated = .FALSE.
     283           0 :       iRowPr = 3 * (iBetaPr-1) + jDirPr
     284           0 :       alphaPrloop: DO iAlphaPr = 1, atoms%nat
     285           0 :          iColPr = 3 * (iAlphaPr-1)
     286           0 :          DO iBeta = 1, iBetaPr
     287           0 :             iRow = 3 * (iBeta-1)
     288           0 :             alphaloop: DO iAlpha = 1, atoms%nat
     289           0 :                iCol = 3 * (iAlpha-1)
     290           0 :                symloop: DO iSym = 1, sym_count
     291           0 :                   IF (.NOT.(iBetaPr==sym%mapped_atom(sym_list(iSym),iBeta))) CYCLE
     292           0 :                   IF (.NOT.(iAlphaPr==sym%mapped_atom(sym_list(iSym),iAlpha))) CYCLE
     293             :                   ! Get all rhs vectors that are possible
     294           0 :                   rhs_vecs = sym_dynvec(iCol+1:iCol+3,:iRowPr-1,iSym)
     295           0 :                   phas = tpi_const*(dot_product(bqpt(:),atoms%taual(:,iAlphaPr)-atoms%taual(:,iBetaPr)))
     296           0 :                   phase_fac = cmplx(cos(phas),sin(phas))
     297             :                   
     298           0 :                   invmrot = sym%mrot(:,:,sym%invtab(sym_list(iSym)))
     299           0 :                   CALL inv3(amat,invamat,det)
     300           0 :                   temp_mat_1 = MATMUL(invmrot,invamat)
     301           0 :                   brot = MATMUL(amat,temp_mat_1)
     302             :                   
     303           0 :                   DO jDir = 1, 3
     304           0 :                      IF (iRow+jDir>iRowPr) CYCLE symloop
     305           0 :                      rhs_vec_full = phase_fac*rhs_vecs(:, iRow + jDir)
     306             :                      !write(*,*) "---------------"
     307             :                      !write(*,*) rhs_vec_full
     308           0 :                      symvec = brot(:, jDir)
     309             :                      !write(*,*) symvec
     310             :                   
     311             :                      l_done = .FALSE.
     312           0 :                      iDone = 0
     313           0 :                      l_need = .TRUE.
     314             :                      iNeed = 3
     315             :                      kNeed = 0
     316           0 :                      DO kDir = 1, 3
     317             :                         !IF (ABS(symvec(kDir))>1e-8.AND.ANY((ABS(dynmat(3 * (iBetaPr-1) + kDir,:))<1e-15)))  
     318           0 :                         IF (ALL(ABS(dynmat(3 * (iBetaPr-1) + kDir,iColPr+1:iColPr+3))>1e-15)) THEN 
     319           0 :                            l_done(kDir) = .TRUE.
     320           0 :                            iDone = iDone + 1
     321           0 :                            l_need(kDir) = .FALSE.
     322           0 :                            iNeed = iNeed - 1
     323           0 :                            CYCLE
     324             :                         END IF
     325           0 :                         IF (ABS(symvec(kDir))<1e-8) THEN
     326           0 :                            l_need(kDir) = .FALSE.
     327           0 :                             iNeed = iNeed - 1
     328           0 :                             CYCLE
     329             :                         END IF
     330           0 :                         kNeed = kDir
     331             :                      END DO
     332             :                      !write(*,*) iDone, l_done
     333             :                      !write(*,*) iNeed, l_need, kNeed
     334           0 :                      IF (iDone==0.OR.iDone==3) CYCLE
     335           0 :                      IF (iNeed/=1) CYCLE
     336           0 :                      DO kDir = 1, 3
     337           0 :                         IF (.NOT.l_need(kDir)) rhs_vec_full = rhs_vec_full - symvec(kDir)*dynmat(3 * (iBetaPr-1) + kDir,iColPr+1:iColPr+3)
     338             :                      END DO
     339           0 :                      IF (SQRT(REAL(DOT_PRODUCT(rhs_vec_full,rhs_vec_full)))<1e-15) CYCLE
     340             :                      !write(*,*) "Newrhs:", rhs_vec_full
     341             :                      !write(*,*) "Thesym:", symvec(kNeed)
     342           0 :                      IF (ABS(symvec(kNeed))>1e-8) THEN 
     343           0 :                         sym_dynmat(iRowPr-jDirPr+kNeed,iColPr+1:iColPr+3) = rhs_vec_full/symvec(kNeed)
     344           0 :                         l_cheated = .TRUE.
     345             :                      ELSE
     346             :                         CYCLE
     347             :                      END IF
     348             :                      !write(*,*) "Cheat: ", sym_dynmat(iRowPr-jDirPr+kNeed,:)
     349           0 :                      IF (ALL(ABS(sym_dynmat(iRowPr-jDirPr+kNeed,iColPr+1:iColPr+3))>1e-15)) CYCLE alphaPrloop
     350             :                   END DO
     351             :                END DO symloop
     352             :             END DO alphaloop
     353             :          END DO
     354             :       END DO alphaPrloop
     355           0 :    END SUBROUTINE
     356             : 
     357           0 : END MODULE m_dfpt_dynmat_sym