rot_den_mat.f 3.53 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13 14 15 16 17
      MODULE m_rotdenmat
      use m_juDFT
      CONTAINS
      SUBROUTINE rot_den_mat(
     >                       alph,beta,
     X                       rho11,rho22,rho21)
c***********************************************************************
c This subroutine rotates the direction of the magnetization of the 
c density matrix by multiplying with the unitary 2x2 spin rotation
c matrix. --> U*rho*U^dagger
c Philipp Kurz 2000-02-03
18
c new method for improved stability (l_new=t) gb'19
19 20
c***********************************************************************

21
      use m_constants
22 23 24 25 26 27 28 29 30 31
      IMPLICIT NONE

C     .. Scalar Arguments ..
      REAL, INTENT    (IN) :: alph,beta
      REAL, INTENT    (INOUT) :: rho11
      REAL, INTENT    (INOUT) :: rho22
      COMPLEX, INTENT (INOUT) :: rho21
C     ..
C     .. Local Scalars ..
      INTEGER ispin
32 33 34
      REAL eps,r11n,r22n
      COMPLEX r21n
      LOGICAL l_new
35 36 37 38
C     ..
C     .. Local Arrays ..
      COMPLEX u2(2,2),rho(2,2),rhoh(2,2)
C     ..
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
      l_new = .true.

      IF (l_new) THEN

      r11n = 0.5*(1.0+cos(beta))*rho11 - sin(beta)*real(rho21) +
     +       0.5*(1.0-cos(beta))*rho22
      r22n = 0.5*(1.0-cos(beta))*rho11 + sin(beta)*real(rho21) +
     +       0.5*(1.0+cos(beta))*rho22
      r21n = CMPLX(cos(alph),-sin(alph))*(sin(beta)*(rho11-rho22) +
     +       2.0*(cos(beta)*real(rho21)-cmplx(0.0,aimag(rho21))))*0.5

      rho11 = r11n
      rho22 = r22n
      rho21 = r21n

      ELSE
55 56 57 58
      
      eps = 1.0e-10

c---> set up the unitary 2x2 spin rotation matrix U^(2)
59 60 61 62
      u2(1,1) =  exp(-ImagUnit*alph/2)*cos(beta/2)
      u2(1,2) = -exp(-ImagUnit*alph/2)*sin(beta/2)
      u2(2,1) =  exp( ImagUnit*alph/2)*sin(beta/2)
      u2(2,2) =  exp( ImagUnit*alph/2)*cos(beta/2)
63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
      
      rho(1,1) = cmplx(rho11,0.0)
      rho(2,2) = cmplx(rho22,0.0)
      rho(2,1) =       rho21
      rho(1,2) = conjg(rho21)

c---> first calculate U*rho
      rhoh(1,1) = u2(1,1)*rho(1,1) + u2(1,2)*rho(2,1)
      rhoh(1,2) = u2(1,1)*rho(1,2) + u2(1,2)*rho(2,2)
      rhoh(2,1) = u2(2,1)*rho(1,1) + u2(2,2)*rho(2,1)
      rhoh(2,2) = u2(2,1)*rho(1,2) + u2(2,2)*rho(2,2)
c---> now calculate (U*rho)*U^dagger
      rho(1,1) = rhoh(1,1)*conjg(u2(1,1))
     +         + rhoh(1,2)*conjg(u2(1,2))
      rho(1,2) = rhoh(1,1)*conjg(u2(2,1))
     +         + rhoh(1,2)*conjg(u2(2,2))
      rho(2,1) = rhoh(2,1)*conjg(u2(1,1))
     +         + rhoh(2,2)*conjg(u2(1,2))
      rho(2,2) = rhoh(2,1)*conjg(u2(2,1))
     +         + rhoh(2,2)*conjg(u2(2,2))

c---> check wether the diagonal elements of the rotated density
c---> are real.
      DO ispin = 1,2
         IF (aimag(rho(ispin,ispin)).GT.eps) THEN
            CALL juDFT_error("rotation of mag. failed",calledby
Daniel Wortmann's avatar
Daniel Wortmann committed
89
     +           ="rot_den_mat",hint=
90 91
     +        'After the rotation of the density matrix in the '//
     +       'muffin-tin sphere one diagonal element of the '//
92
     +       '(hermitian) density matrix is not real. That means '//
93
     +       'that the density matrix was probably damaged.')
94 95 96
         ENDIF
      ENDDO
 
97 98 99 100
      rho11 = real(rho(1,1))
      rho22 = real(rho(2,2))
      rho21 =      rho(2,1)

101 102
      ENDIF
      
103 104 105
      END SUBROUTINE rot_den_mat
      END MODULE m_rotdenmat