LCOV - code coverage report
Current view: top level - vgen - rotate_int_den_to_local.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 46 136 33.8 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_rotate_int_den_to_local
       8             :   USE m_juDFT
       9             :   !**********************************************************************
      10             :   !     This subroutine calculates the spin-up and -down density, in the INT-region,
      11             :   !     e.i. it take the non-colinear density and rotates it locally into the
      12             :   !     spin-frame that make it spin-diagonal.
      13             :   !     The rotated density is needed to calculate the potential-energy integrals 
      14             :   !     in vgen_xcpot. For accuracy reasons, the magnetisation for the potential
      15             :   !     itself is regeneated from the unrotated densities.
      16             :   !     In addition this routine stores the angle used in the rotation.
      17             :   !     These angles are needed in vgen->vgen_finalize->vmatgen to rotate the up- and down-
      18             :   !     potentials back to the global frame.   DW 2018
      19             :   !     Based on rhodirgen by
      20             :   !     Philipp Kurz 99/11/01
      21             :   !**********************************************************************
      22             : CONTAINS
      23         248 :   SUBROUTINE rotate_int_den_to_local(DIMENSION,sym,stars,atoms,sphhar,vacuum,&
      24             :        cell,input,noco,oneD,den)
      25             :     !******** ABBREVIATIONS ***********************************************
      26             :     !     ifft3    : size of the 3d real space mesh
      27             :     !     ifft2    : size of the 2d real space mesh
      28             :     !     rpw      : diagonal components of the density matrix (rho_11 ,
      29             :     !                rho_22)
      30             :     !                later interstitial spin-up and -down density
      31             :     !                all stored in terms of 3d-stars
      32             :     !     ris      : first components of the density matrix
      33             :     !                later interstitial spin-up and -down density and
      34             :     !                direction of magnetic field (theta and phi)
      35             :     !                all stored on real space mesh
      36             :     !**********************************************************************
      37             : 
      38             :     USE m_constants
      39             :     USE m_fft2d
      40             :     USE m_fft3d
      41             :     USE m_types
      42             :     IMPLICIT NONE
      43             : 
      44             :     TYPE(t_dimension),INTENT(IN)   :: DIMENSION
      45             :     TYPE(t_noco),INTENT(IN)        :: noco
      46             :     TYPE(t_oneD),INTENT(IN)        :: oneD
      47             :     TYPE(t_input),INTENT(IN)       :: input
      48             :     TYPE(t_vacuum),INTENT(IN)      :: vacuum
      49             :     TYPE(t_sym),INTENT(IN)         :: sym
      50             :     TYPE(t_stars),INTENT(IN)       :: stars
      51             :     TYPE(t_cell),INTENT(IN)        :: cell
      52             :     TYPE(t_sphhar),INTENT(IN)      :: sphhar
      53             :     TYPE(t_atoms),INTENT(IN)       :: atoms
      54             :     TYPE(t_potden),INTENT(INOUT)   :: den
      55             : 
      56             :     !     .. Local Scalars ..
      57             :     INTEGER iden,jspin,ivac,ifft2,ifft3
      58             :     INTEGER imz,ityp,iri,ilh,imesh,iq2,iq3
      59             :     REAL   rho_11,rho_22,rho_21r,rho_21i,rhotot,magmom,phi
      60             :     REAL rho_up,rho_down,mx,my,mz,eps,vz_r,vz_i,rziw,theta
      61             :     !     ..
      62             :     !     .. Local Arrays ..
      63         248 :     REAL,    ALLOCATABLE :: rz(:,:,:)
      64         248 :     REAL,    ALLOCATABLE :: rvacxy(:,:,:,:),ris(:,:),fftwork(:)
      65             :     !     ..
      66         248 :     eps = 1.0e-20
      67             : 
      68             : 
      69             : 
      70             :     !
      71             :     !---> initialize arrays for the density matrix
      72             :     !
      73             : 
      74         248 :     ifft3 = 27*stars%mx1*stars%mx2*stars%mx3
      75         248 :     IF (input%film) THEN
      76           0 :        ifft2 = 9*stars%mx1*stars%mx2
      77           0 :        IF (oneD%odi%d1) ifft2 = 9*stars%mx3*oneD%odi%M
      78             :     ELSE
      79             :        ifft2=0
      80             :     END IF
      81             : 
      82         248 :     IF (ALLOCATED(den%phi_pw)) THEN
      83           0 :        DEALLOCATE(den%phi_pw,den%phi_vacz,den%phi_vacxy)
      84           0 :        DEALLOCATE(den%theta_pw,den%theta_vacz,den%theta_vacxy)
      85             :     ENDIF
      86         248 :     ALLOCATE(den%phi_pw(ifft3),den%theta_pw(ifft3))
      87         248 :     ALLOCATE(den%phi_vacz(vacuum%nmzd,2),den%theta_vacz(vacuum%nmzd,2))
      88         248 :     ALLOCATE(den%phi_vacxy(ifft2,vacuum%nmzxyd,2),den%theta_vacxy(ifft2,vacuum%nmzxyd,2))
      89             : 
      90             :      
      91         248 :     ALLOCATE (ris(ifft3,4),fftwork(ifft3))
      92             :     !---> fouriertransform the diagonal part of the density matrix
      93             :     !---> in the interstitial, den%pw, to real space (ris)
      94        1240 :     DO iden = 1,2
      95         744 :        CALL fft3d(ris(:,iden),fftwork,den%pw(:,iden),stars,+1)
      96             :     ENDDO
      97             :     !---> fouriertransform the off-diagonal part of the density matrix
      98         248 :     CALL fft3d(ris(:,3),ris(:,4),den%pw(:,3),stars,+1)
      99             : 
     100             :     !test
     101             :     !      DO iden=1,4
     102             :     !         write(*,*)'iden=',iden
     103             :     !         write(*,8500)(ris(imesh,iden),imesh=0,ifft3-1)
     104             :     !      enddo
     105             :     !test
     106             :     !---> calculate the charge and magnetization density on the
     107             :     !---> real space mesh
     108     3397712 :     DO imesh = 1,ifft3
     109     3397464 :        rho_11  = ris(imesh,1)
     110     3397464 :        rho_22  = ris(imesh,2)
     111     3397464 :        rho_21r = ris(imesh,3)
     112     3397464 :        rho_21i = ris(imesh,4)
     113     3397464 :        mx      =  2*rho_21r
     114     3397464 :        my      = -2*rho_21i
     115     3397464 :        mz      = (rho_11-rho_22)
     116     3397464 :        magmom  = SQRT(mx**2 + my**2 + mz**2)
     117     3397464 :        rhotot  = rho_11 + rho_22
     118     3397464 :        rho_up  = (rhotot + magmom)/2
     119     3397464 :        rho_down= (rhotot - magmom)/2
     120             : 
     121     3397464 :        IF (ABS(mz) .LE. eps) THEN
     122             :           theta = pi_const/2
     123     3397464 :        ELSEIF (mz .GE. 0.0) THEN
     124     1951282 :           theta = ATAN(SQRT(mx**2 + my**2)/mz)
     125             :        ELSE
     126     1446182 :           theta = ATAN(SQRT(mx**2 + my**2)/mz) + pi_const
     127             :        ENDIF
     128             : 
     129     3397464 :        IF (ABS(mx) .LE. eps) THEN
     130      193158 :           IF (ABS(my) .LE. eps) THEN
     131             :              phi = 0.0
     132           0 :           ELSEIF (my .GE. 0.0) THEN
     133             :              phi = pi_const/2
     134             :           ELSE
     135           0 :              phi = -pi_const/2
     136             :           ENDIF
     137     3204306 :        ELSEIF (mx .GE. 0.0) THEN
     138     1688082 :           phi = ATAN(my/mx)
     139             :        ELSE
     140     1516224 :           IF (my .GE. 0.0) THEN
     141      825740 :              phi = ATAN(my/mx) + pi_const
     142             :           ELSE
     143      690484 :              phi = ATAN(my/mx) - pi_const
     144             :           ENDIF
     145             :        ENDIF
     146             : 
     147             :        !         write(36,'(i4,2f12.6)') mod(imesh,33),rho_11,rho_22
     148     3397464 :        ris(imesh,1) = rho_up
     149     3397464 :        ris(imesh,2) = rho_down
     150     3397464 :        den%theta_pw(imesh) = theta
     151     3397712 :        den%phi_pw(imesh) = phi
     152             :     ENDDO
     153             : 
     154         744 :     DO jspin = 1,input%jspins
     155     6795424 :        fftwork=0.0
     156         744 :        CALL fft3d(ris(:,jspin),fftwork,den%pw(:,jspin),stars,-1)
     157             :     ENDDO
     158             : 
     159         248 :     IF (.NOT.input%film) RETURN
     160             : 
     161             : 
     162             :      !Now the vacuum part starts
     163             : 
     164             :    
     165           0 :     ALLOCATE(rvacxy(ifft2,vacuum%nmzxyd,2,4))
     166           0 :     ALLOCATE (rz(vacuum%nmzd,2,2))
     167             :     !---> fouriertransform the diagonal part of the density matrix
     168             :     !---> in the vacuum, rz & rxy, to real space (rvacxy)
     169           0 :     DO iden = 1,2
     170           0 :        DO ivac = 1,vacuum%nvac
     171           0 :           DO imz = 1,vacuum%nmzxyd
     172           0 :              rziw = 0.0
     173           0 :              IF (oneD%odi%d1) THEN
     174           0 :                 CALL judft_error("oneD not implemented",calledby="rhodirgen")
     175             :                 !CALL fft2d(oneD%k3,odi%M,odi%n2d,rvacxy(0,imz,ivac,iden),fftwork,&
     176             :                 !           den%vacz(imz,ivac,iden),rziw,den%vacxy(imz,1,ivac,iden),&
     177             :                 !           vacuum,odi%nq2,odi%kimax2,1,&
     178             :                 !     &                  %igf,odl%pgf,odi%nst2)
     179             :              ELSE
     180             :                 CALL fft2d(stars,rvacxy(:,imz,ivac,iden),fftwork,&
     181             :                      den%vacz(imz,ivac,iden),rziw,den%vacxy(imz,1,ivac,iden),&
     182           0 :                      vacuum%nmzxyd,1)
     183             :              ENDIF
     184             :           ENDDO
     185             :        ENDDO
     186             :     ENDDO
     187             :     !--->    fouriertransform the off-diagonal part of the density matrix
     188           0 :     DO ivac = 1,vacuum%nvac
     189           0 :        DO imz = 1,vacuum%nmzxyd
     190           0 :           rziw = 0.0
     191           0 :           vz_r = den%vacz(imz,ivac,3)
     192           0 :           vz_i = den%vacz(imz,ivac,4)
     193           0 :           IF (oneD%odi%d1) THEN
     194           0 :              CALL judft_error("oneD not implemented",calledby="rhodirgen")
     195             :              !CALL fft2d(oneD%k3,odi%M,odi%n2d,&
     196             :              !           rvacxy(0,imz,ivac,3),rvacxy(0,imz,ivac,4),&
     197             :              !           vz_r,vz_i,den%vacxy(imz,1,ivac,3),&
     198             :              !           vacuum,odi%nq2,odi%kimax2,1,&
     199             :              !     &               %igf,odl%pgf,odi%nst2)
     200             :           ELSE
     201             :              CALL fft2d(stars,rvacxy(:,imz,ivac,3),rvacxy(:,imz,ivac,4),&
     202           0 :                   vz_r,vz_i,den%vacxy(imz,1,ivac,3),vacuum%nmzxyd,1)
     203             :           ENDIF
     204             :        ENDDO
     205             :     ENDDO
     206             : 
     207             :     !--->    calculate the four components of the matrix potential on
     208             :     !--->    real space mesh
     209           0 :     DO ivac = 1,vacuum%nvac
     210           0 :        DO imz = 1,vacuum%nmzxyd
     211           0 :           DO imesh = 1,ifft2
     212           0 :              rho_11  = rvacxy(imesh,imz,ivac,1)
     213           0 :              rho_22  = rvacxy(imesh,imz,ivac,2)
     214           0 :              rho_21r = rvacxy(imesh,imz,ivac,3)
     215           0 :              rho_21i = rvacxy(imesh,imz,ivac,4)
     216           0 :              mx      =  2*rho_21r
     217           0 :              my      = -2*rho_21i
     218           0 :              mz      = (rho_11-rho_22)
     219           0 :              magmom  = SQRT(mx**2 + my**2 + mz**2)
     220           0 :              rhotot  = rho_11 + rho_22
     221           0 :              rho_up  = (rhotot + magmom)/2
     222           0 :              rho_down= (rhotot - magmom)/2
     223             : 
     224           0 :              IF (ABS(mz) .LE. eps) THEN
     225             :                 theta = pi_const/2
     226           0 :              ELSEIF (mz .GE. 0.0) THEN
     227           0 :                 theta = ATAN(SQRT(mx**2 + my**2)/mz)
     228             :              ELSE
     229           0 :                 theta = ATAN(SQRT(mx**2 + my**2)/mz) + pi_const
     230             :              ENDIF
     231             : 
     232           0 :              IF (ABS(mx) .LE. eps) THEN
     233           0 :                 IF (ABS(my) .LE. eps) THEN
     234             :                    phi = 0.0
     235           0 :                 ELSEIF (my .GE. 0.0) THEN
     236             :                    phi = pi_const/2
     237             :                 ELSE
     238           0 :                    phi = -pi_const/2
     239             :                 ENDIF
     240           0 :              ELSEIF (mx .GE. 0.0) THEN
     241           0 :                 phi = ATAN(my/mx)
     242             :              ELSE
     243           0 :                 IF (my .GE. 0.0) THEN
     244           0 :                    phi = ATAN(my/mx) + pi_const
     245             :                 ELSE
     246           0 :                    phi = ATAN(my/mx) - pi_const
     247             :                 ENDIF
     248             :              ENDIF
     249             : 
     250           0 :              rvacxy(imesh,imz,ivac,1) = rho_up
     251           0 :              rvacxy(imesh,imz,ivac,2) = rho_down
     252           0 :              den%theta_vacxy(imesh,imz,ivac) = theta
     253           0 :              den%phi_vacxy(imesh,imz,ivac) = phi
     254             :           ENDDO
     255             :        ENDDO
     256           0 :        DO imz = vacuum%nmzxyd+1,vacuum%nmzd
     257           0 :           rho_11  = den%vacz(imz,ivac,1)
     258           0 :           rho_22  = den%vacz(imz,ivac,2)
     259           0 :           rho_21r = den%vacz(imz,ivac,3)
     260           0 :           rho_21i = den%vacz(imz,ivac,4)
     261           0 :           mx      =  2*rho_21r
     262           0 :           my      = -2*rho_21i
     263           0 :           mz      = (rho_11-rho_22)
     264           0 :           magmom  = SQRT(mx**2 + my**2 + mz**2)
     265           0 :           rhotot  = rho_11 + rho_22
     266           0 :           rho_up  = (rhotot + magmom)/2
     267           0 :           rho_down= (rhotot - magmom)/2
     268             : 
     269           0 :           IF (ABS(mz) .LE. eps) THEN
     270             :              theta = pi_const/2
     271           0 :           ELSEIF (mz .GE. 0.0) THEN
     272           0 :              theta = ATAN(SQRT(mx**2 + my**2)/mz)
     273             :           ELSE
     274           0 :              theta = ATAN(SQRT(mx**2 + my**2)/mz) + pi_const
     275             :           ENDIF
     276             : 
     277           0 :           IF (ABS(mx) .LE. eps) THEN
     278           0 :              IF (ABS(my) .LE. eps) THEN
     279             :                 phi = 0.0
     280           0 :              ELSEIF (my .GE. 0.0) THEN
     281             :                 phi = pi_const/2
     282             :              ELSE
     283           0 :                 phi = -pi_const/2
     284             :              ENDIF
     285           0 :           ELSEIF (mx .GE. 0.0) THEN
     286           0 :              phi = ATAN(my/mx)
     287             :           ELSE
     288           0 :              IF (my .GE. 0.0) THEN
     289           0 :                 phi = ATAN(my/mx) + pi_const
     290             :              ELSE
     291           0 :                 phi = ATAN(my/mx) - pi_const
     292             :              ENDIF
     293             :           ENDIF
     294             : 
     295           0 :           den%vacz(imz,ivac,1) = rho_up
     296           0 :           den%vacz(imz,ivac,2) = rho_down
     297           0 :           den%theta_vacz(imz,ivac) = theta
     298           0 :           den%phi_vacz(imz,ivac) = phi
     299             :        ENDDO
     300             :     ENDDO
     301             :     !--->    Fouriertransform the matrix potential back to reciprocal space
     302           0 :     DO jspin = 1,input%jspins
     303           0 :        DO ivac = 1,vacuum%nvac
     304           0 :           DO imz = 1,vacuum%nmzxyd
     305           0 :              fftwork=0.0
     306           0 :              IF (oneD%odi%d1) THEN
     307           0 :                 call judft_error("oneD not implemented",calledby="rhodirgen")
     308             :                 !CALL fft2d(oneD%k3,odi%M,odi%n2d,&
     309             :                 !           rvacxy(0,imz,ivac,jspin),fftwork,&
     310             :                 !           den%vacz(imz,ivac,jspin),rziw,den%vacxy(imz,1,ivac,jspin),&
     311             :                 !           vacuum,odi%nq2,odi%kimax2,-1,&
     312             :                 !     &                  %igf,odl%pgf,odi%nst2)
     313             :              ELSE
     314             :                 CALL fft2d(stars,rvacxy(:,imz,ivac,jspin),fftwork,&
     315             :                      den%vacz(imz,ivac,jspin),rziw,den%vacxy(imz,1,ivac,jspin),&
     316           0 :                      vacuum%nmzxyd,-1)
     317             :              END IF
     318             :           ENDDO
     319             :        ENDDO
     320             :     ENDDO
     321             :     
     322             :     RETURN
     323         248 :   END SUBROUTINE rotate_int_den_to_local
     324             : END MODULE m_rotate_int_den_to_local

Generated by: LCOV version 1.13