LCOV - code coverage report
Current view: top level - optional - pldngen.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 172 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.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_pldngen
       8             : 
       9             : !**********************************************************************
      10             : !     This subroutine generates the charge and magetization densities
      11             : !     (mx,my,mz) and writes them to the files cdn, mdnx, mdny, mdnz.
      12             : !     These files are needed to generate plots of the density.
      13             : !
      14             : !    i) The components of the hermitian density matrix (rho_11, rho_22,
      15             : !     rho_21) are reloaded from the file rhomat_inp.
      16             : !    ii) The density matrix in fouriertransformed to real space.
      17             : !    iii) The charge and magnetization density (n, mx, my, mz) are
      18             : !     calculated on the real space mesh.
      19             : !    iv) n, mx, my, and mz are Fouriertransformed and stored in terms
      20             : !     of stars.
      21             : !
      22             : !     Philipp Kurz 99/10/29
      23             : !**********************************************************************
      24             : 
      25             : CONTAINS
      26             : 
      27           0 : SUBROUTINE pldngen(mpi,sym,stars,atoms,sphhar,vacuum,&
      28             :                    cell,input,noco,oneD,sliceplot)
      29             : 
      30             :           !******** ABBREVIATIONS ***********************************************
      31             :           !     ifft3    : size of the 3d real space mesh
      32             :           !     ifft2    : size of the 2d real space mesh
      33             :           !     rpw      : first diagonal components of the interstitial density
      34             :           !                matrix
      35             :           !                later charge and mag. density (n, mx, my, mz)
      36             :           !                all stored in terms of 3d-stars
      37             :           !     ris      : first componets of the density matrix
      38             :           !                later charge and mag. density (n, mx, my, mz)
      39             :           !                all stored on real space mesh
      40             :           !**********************************************************************
      41             : 
      42             :    USE m_juDFT
      43             :    USE m_constants
      44             :    USE m_cdn_io
      45             :    USE m_loddop
      46             :    USE m_wrtdop
      47             :    USE m_qfix
      48             :    USE m_fft2d
      49             :    USE m_fft3d
      50             :    USE m_types
      51             :    USE m_rotdenmat 
      52             : 
      53             :    IMPLICIT NONE
      54             : 
      55             :    TYPE(t_mpi),INTENT(IN)    :: mpi
      56             :    TYPE(t_sym),INTENT(IN)    :: sym
      57             :    TYPE(t_stars),INTENT(IN)  :: stars
      58             :    TYPE(t_vacuum),INTENT(IN) :: vacuum
      59             :    TYPE(t_atoms),INTENT(IN)  :: atoms
      60             :    TYPE(t_sphhar),INTENT(IN) :: sphhar
      61             :    TYPE(t_input),INTENT(IN)  :: input
      62             :    TYPE(t_cell),INTENT(IN)   :: cell
      63             :    TYPE(t_oneD),INTENT(IN)    :: oneD
      64             :    TYPE(t_noco),INTENT(IN)   :: noco
      65             :    TYPE(t_sliceplot),INTENT(IN):: sliceplot
      66             : 
      67             :    ! Local type instances
      68             :    TYPE(t_input)  :: inp
      69           0 :    TYPE(t_potden) :: den
      70             : 
      71             :    ! Local Scalars
      72             :    INTEGER iden,ivac,ifft2,ifft3,archiveType
      73             :    INTEGER imz,ityp,iri,ilh,imesh,iter
      74             :    REAL cdnup,cdndown,chden,mgden,theta,phi,zero,rho_11,rziw,fermiEnergyTemp
      75             :    REAL rho_22,rho_21r,rho_21i,rhotot,mx,my,mz,fix,vz_r,vz_i
      76             :    COMPLEX czero
      77             : 
      78             :    ! Local Arrays
      79             :    !---> off-diagonal part of the density matrix
      80           0 :    COMPLEX, ALLOCATABLE :: cdom(:),cdomvz(:,:),cdomvxy(:,:,:)
      81           0 :    COMPLEX, ALLOCATABLE :: qpw(:,:),rhtxy(:,:,:,:)
      82           0 :    REAL,    ALLOCATABLE :: rht(:,:,:),rho(:,:,:,:)
      83           0 :    REAL,    ALLOCATABLE :: rvacxy(:,:,:,:),ris(:,:),fftwork(:)
      84             : 
      85             :    !---> for testing: output of offdiag. output density matrix. to plot the
      86             :    !---> offdiag. part of the output density matrix, that part has to be
      87             :    !---> written the file rhomt21 in cdnmt.
      88             :    LOGICAL :: l_qfix
      89             :    REAL    :: cdn11, cdn22
      90             :    COMPLEX :: cdn21
      91             :    !---> end of test part
      92             : 
      93           0 :    iter = 0 ! This is not clean!
      94             : 
      95           0 :    zero = 0.0 ; czero = CMPLX(0.0,0.0)
      96           0 :    ifft3 = 27*stars%mx1*stars%mx2*stars%mx3
      97           0 :    ifft2 = 9*stars%mx1*stars%mx2
      98             : 
      99             :    ALLOCATE (qpw(stars%ng3,4),rhtxy(vacuum%nmzxyd,stars%ng2-1,2,4),&
     100             :              cdom(stars%ng3),cdomvz(vacuum%nmzd,2),cdomvxy(vacuum%nmzxyd,stars%ng2-1,2),&
     101             :              ris(0:27*stars%mx1*stars%mx2*stars%mx3-1,4),fftwork(0:27*stars%mx1*stars%mx2*stars%mx3-1),&
     102             :              rvacxy(0:9*stars%mx1*stars%mx2-1,vacuum%nmzxyd,2,4),&
     103           0 :              rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,4),rht(vacuum%nmzd,2,4) )
     104             : 
     105             :    !---> initialize arrays for the density matrix
     106           0 :    rho(:,:,:,:) = zero ; qpw(:,:) = czero ; cdom(:) = czero
     107           0 :    IF (input%film) THEN
     108           0 :       cdomvz(:,:) = czero ;    rhtxy(:,:,:,:) = czero
     109           0 :       cdomvxy(:,:,:) = czero ; rht(:,:,:) = zero
     110             :    END IF
     111             : 
     112           0 :    IF (input%jspins .NE. 2) THEN
     113           0 :       WRITE (6,*) 'This is the non-collinear version of the flapw-'
     114           0 :       WRITE (6,*) 'program. It can only perform spin-polarized'
     115           0 :       WRITE (6,*) 'calculations.'
     116             :       CALL juDFT_error("jspins not equal 2",calledby = "pldngen",hint=&
     117             :                        "This is the non-collinear version of the flapw-"//&
     118             :                        ' PROGRAM. It can ONLY perform spin-polarized '//&
     119           0 :                        'calculations.')
     120             :    END IF
     121             : 
     122             :    !---> reload the density matrix from file rhomat_inp
     123           0 :    archiveType = CDN_ARCHIVE_TYPE_CDN1_const
     124           0 :    IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
     125           0 :    CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
     126           0 :    IF (.NOT.sliceplot%slice) THEN
     127             :       CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
     128           0 :                        0,fermiEnergyTemp,l_qfix,den)
     129             :    ELSE
     130             :       CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
     131           0 :                        0,fermiEnergyTemp,l_qfix,den,'cdn_slice')
     132             :    END IF
     133           0 :    rho(:,0:,1:,:input%jspins) = den%mt(:,0:,1:,:input%jspins)
     134           0 :    qpw(1:,:input%jspins) = den%pw(1:,:input%jspins)
     135           0 :    rht(1:,1:,:input%jspins) = den%vacz(1:,1:,:input%jspins)
     136           0 :    rhtxy(1:,1:,1:,:input%jspins) = den%vacxy(1:,1:,1:,:input%jspins)
     137           0 :    IF(noco%l_noco) THEN
     138           0 :       cdom = den%pw(:,3)
     139           0 :       cdomvz(:,:) = CMPLX(den%vacz(:,:,3),den%vacz(:,:,4))
     140           0 :       cdomvxy = den%vacxy(:,:,:,3)
     141             :    END IF
     142             : 
     143           0 :    IF (.NOT. sliceplot%slice) THEN
     144           0 :       CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
     145           0 :       den%iter = iter
     146           0 :       den%mt(:,0:,1:,:input%jspins) = rho(:,0:,1:,:input%jspins)
     147           0 :       den%pw(1:,:input%jspins) = qpw(1:,:input%jspins)
     148           0 :       den%vacz(1:,1:,:input%jspins) = rht(1:,1:,:input%jspins)
     149           0 :       den%vacxy(1:,1:,1:,:input%jspins) = rhtxy(1:,1:,1:,:input%jspins)
     150           0 :       IF(noco%l_noco) THEN
     151           0 :          den%pw(:,3) = cdom
     152           0 :          den%vacz(:,:,3) = REAL(cdomvz(:,:))
     153           0 :          den%vacz(:,:,4) = AIMAG(cdomvz(:,:))
     154           0 :          den%vacxy(:,:,:,3) = cdomvxy
     155             :       END IF
     156           0 :       CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,den,noco%l_noco,.FALSE.,.true.,fix)
     157           0 :       rho(:,0:,1:,:input%jspins) = den%mt(:,0:,1:,:input%jspins)
     158           0 :       qpw(1:,:input%jspins) = den%pw(1:,:input%jspins)
     159           0 :       rht(1:,1:,:input%jspins) = den%vacz(1:,1:,:input%jspins)
     160           0 :       rhtxy(1:,1:,1:,:input%jspins) = den%vacxy(1:,1:,1:,:input%jspins)
     161           0 :       IF(noco%l_noco) THEN
     162           0 :          cdom = den%pw(:,3)
     163           0 :          cdomvz(:,:) = CMPLX(den%vacz(:,:,3),den%vacz(:,:,4))
     164           0 :          cdomvxy = den%vacxy(:,:,:,3)
     165             :       END IF
     166             :    END IF
     167             : 
     168             :    !---> calculate the charge and magnetization density in the muffin tins
     169           0 :    DO ityp = 1,atoms%ntype
     170           0 :       DO ilh = 0,sphhar%nlh(atoms%ntypsy(ityp))
     171           0 :          DO iri = 1,atoms%jri(ityp)
     172           0 :             IF (SIZE(den%mt,4).LE.2) THEN 
     173           0 :                cdnup   = rho(iri,ilh,ityp,1)
     174           0 :                cdndown = rho(iri,ilh,ityp,2)
     175           0 :                theta = noco%beta(ityp)
     176           0 :                phi   = noco%alph(ityp)
     177           0 :                chden  = cdnup + cdndown
     178           0 :                mgden  = cdnup - cdndown
     179           0 :                rho(iri,ilh,ityp,1) = chden
     180           0 :                rho(iri,ilh,ityp,2) = mgden*COS(phi)*SIN(theta)
     181           0 :                rho(iri,ilh,ityp,3) = mgden*SIN(phi)*SIN(theta)
     182           0 :                rho(iri,ilh,ityp,4) = mgden*COS(theta)
     183             :             ELSE 
     184             :                !--->            for testing: output of offdiag. output density matrix
     185           0 :                cdn11 = rho(iri,ilh,ityp,1)
     186           0 :                cdn22 = rho(iri,ilh,ityp,2)
     187           0 :                cdn21 = CMPLX(den%mt(iri,ilh,ityp,3),den%mt(iri,ilh,ityp,4))
     188           0 :                CALL rot_den_mat(noco%alph(ityp),noco%beta(ityp),cdn11,cdn22,cdn21)
     189           0 :                rho(iri,ilh,ityp,1) = cdn11 + cdn22
     190           0 :                rho(iri,ilh,ityp,2) = 2.0*REAL(cdn21)
     191             :                ! Note: The minus sign in the following line is temporary to adjust for differences in the offdiagonal
     192             :                !       part of the density between this fleur version and ancient (v0.26) fleur.
     193           0 :                rho(iri,ilh,ityp,3) = -2.0*AIMAG(cdn21)
     194           0 :                rho(iri,ilh,ityp,4) = cdn11 - cdn22
     195             :                !--->            end of test part
     196             :             END IF
     197             :          END DO
     198             :       END DO
     199             :    END DO
     200             : 
     201             : 
     202             :    !---> fouriertransform the diagonal part of the density matrix
     203             :    !---> in the interstitial, qpw, to real space (ris)
     204           0 :    DO iden = 1,2
     205           0 :       CALL fft3d(ris(0,iden),fftwork,qpw(1,iden),stars,1)
     206             :    END DO
     207             :    !---> fouriertransform the off-diagonal part of the density matrix
     208           0 :    CALL fft3d(ris(0,3),ris(0,4),cdom(1),stars,+1)
     209             : 
     210             :    !---> calculate the charge and magnetization density on the
     211             :    !---> real space mesh
     212           0 :    DO imesh = 0,ifft3-1
     213           0 :       rho_11  = ris(imesh,1)
     214           0 :       rho_22  = ris(imesh,2)
     215           0 :       rho_21r = ris(imesh,3)
     216           0 :       rho_21i = ris(imesh,4)
     217           0 :       rhotot  = rho_11 + rho_22
     218           0 :       mx      =  2*rho_21r
     219           0 :       my      = -2*rho_21i
     220           0 :       mz      = (rho_11-rho_22)
     221             : 
     222           0 :       ris(imesh,1) = rhotot
     223           0 :       ris(imesh,2) = mx
     224           0 :       ris(imesh,3) = my
     225           0 :       ris(imesh,4) = mz
     226             :    END DO
     227             : 
     228             :    !---> Fouriertransform the density matrix back to reciprocal space
     229           0 :    DO iden = 1,4
     230           0 :       fftwork=zero
     231           0 :       CALL fft3d(ris(0,iden),fftwork,qpw(1,iden),stars,-1)
     232             :    END DO
     233             : 
     234             :    !---> fouriertransform the diagonal part of the density matrix
     235             :    !---> in the vacuum, rz & rxy, to real space (rvacxy)
     236           0 :    IF (input%film) THEN
     237           0 :       DO iden = 1,2
     238           0 :          DO ivac = 1,vacuum%nvac
     239           0 :             DO imz = 1,vacuum%nmzxyd
     240           0 :                rziw = 0.0
     241             :                CALL fft2d(stars,rvacxy(0,imz,ivac,iden),fftwork,rht(imz,ivac,iden),&
     242           0 :                           rziw,rhtxy(imz,1,ivac,iden),vacuum%nmzxyd,1)
     243             :             END DO
     244             :          END DO
     245             :       END DO
     246             :       !--->    fouriertransform the off-diagonal part of the density matrix
     247           0 :       DO ivac = 1,vacuum%nvac
     248           0 :          DO imz = 1,vacuum%nmzxyd
     249           0 :             rziw = 0.0
     250           0 :             vz_r = REAL(cdomvz(imz,ivac))
     251           0 :             vz_i = AIMAG(cdomvz(imz,ivac))
     252             :             CALL fft2d(stars,rvacxy(0,imz,ivac,3),rvacxy(0,imz,ivac,4),&
     253           0 :                        vz_r,vz_i,cdomvxy(imz,1,ivac),vacuum%nmzxyd,1)
     254             :          END DO
     255             :       END DO
     256             : 
     257             :       !--->    calculate the four components of the matrix potential on
     258             :       !--->    real space mesh
     259           0 :       DO ivac = 1,vacuum%nvac
     260           0 :          DO imz = 1,vacuum%nmzxyd
     261           0 :             DO imesh = 0,ifft2-1
     262           0 :                rho_11  = rvacxy(imesh,imz,ivac,1)
     263           0 :                rho_22  = rvacxy(imesh,imz,ivac,2)
     264           0 :                rho_21r = rvacxy(imesh,imz,ivac,3)
     265           0 :                rho_21i = rvacxy(imesh,imz,ivac,4)
     266           0 :                rhotot  = rho_11 + rho_22
     267           0 :                mx      =  2*rho_21r
     268           0 :                my      = -2*rho_21i
     269           0 :                mz      = (rho_11-rho_22)
     270             : 
     271           0 :                rvacxy(imesh,imz,ivac,1) = rhotot
     272           0 :                rvacxy(imesh,imz,ivac,2) = mx
     273           0 :                rvacxy(imesh,imz,ivac,3) = my
     274           0 :                rvacxy(imesh,imz,ivac,4) = mz
     275             :             END DO
     276             :          END DO
     277           0 :          DO imz = vacuum%nmzxyd+1,vacuum%nmzd
     278           0 :             rho_11  = rht(imz,ivac,1)
     279           0 :             rho_22  = rht(imz,ivac,2)
     280           0 :             rho_21r = REAL(cdomvz(imz,ivac))
     281           0 :             rho_21i = AIMAG(cdomvz(imz,ivac))
     282           0 :             rhotot  = rho_11 + rho_22
     283           0 :             mx      =  2*rho_21r
     284           0 :             my      = -2*rho_21i
     285           0 :             mz      = (rho_11-rho_22)
     286             : 
     287           0 :             rht(imz,ivac,1) = rhotot
     288           0 :             rht(imz,ivac,2) = mx
     289           0 :             rht(imz,ivac,3) = my
     290           0 :             rht(imz,ivac,4) = mz
     291             :          END DO
     292             :       END DO
     293             :       !--->    Fouriertransform the matrix potential back to reciprocal space
     294           0 :       DO iden = 1,4
     295           0 :          DO ivac = 1,vacuum%nvac
     296           0 :             DO imz = 1,vacuum%nmzxyd
     297           0 :                fftwork=zero
     298             :                CALL fft2d(stars,rvacxy(0,imz,ivac,iden),fftwork,rht(imz,ivac,iden),&
     299           0 :                           rziw,rhtxy(imz,1,ivac,iden),vacuum%nmzxyd,-1)
     300             :             END DO
     301             :          END DO
     302             :       END DO
     303             :    END IF
     304             : 
     305             :    !---> save charge density to file cdn
     306           0 :    inp=input
     307           0 :    inp%jspins=1
     308             : 
     309           0 :    CALL den%init(stars,atoms,sphhar,vacuum,noco,inp%jspins,POTDEN_TYPE_DEN)
     310           0 :    den%iter = iter
     311           0 :    den%mt(:,0:,1:,1:1) = rho(:,0:,1:,1:1)
     312           0 :    den%pw(1:,1:1) = qpw(1:,1:1)
     313           0 :    den%vacz(1:,1:,1:1) = rht(1:,1:,1:1)
     314           0 :    den%vacxy(1:,1:,1:,1:1) = rhtxy(1:,1:,1:,1:1)
     315           0 :    IF(noco%l_noco) THEN
     316           0 :       den%pw(:,3) = cdom
     317           0 :       den%vacz(:,:,3) = REAL(cdomvz(:,:))
     318           0 :       den%vacz(:,:,4) = AIMAG(cdomvz(:,:))
     319           0 :       den%vacxy(:,:,:,3) = cdomvxy
     320             :    END IF
     321             : 
     322             :    CALL writeDensity(stars,vacuum,atoms,cell,sphhar,inp,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
     323           0 :                      0,-1.0,0.0,.FALSE.,den)
     324             : 
     325             :    !---> save mx to file mdnx
     326           0 :    den%mt(:,0:,1:,1) = rho(:,0:,1:,2)
     327           0 :    den%pw(1:,1) = qpw(1:,2)
     328           0 :    den%vacz(1:,1:,1) = rht(1:,1:,2)
     329           0 :    den%vacxy(1:,1:,1:,1) = rhtxy(1:,1:,1:,2)
     330             :    CALL writeDensity(stars,vacuum,atoms,cell,sphhar,inp,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
     331           0 :                      0,-1.0,0.0,.FALSE.,den,'mdnx')
     332             : 
     333             :    !---> save my to file mdny
     334           0 :    den%mt(:,0:,1:,1) = rho(:,0:,1:,3)
     335           0 :    den%pw(1:,1) = qpw(1:,3)
     336           0 :    den%vacz(1:,1:,1) = rht(1:,1:,3)
     337           0 :    den%vacxy(1:,1:,1:,1) = rhtxy(1:,1:,1:,3)
     338             :    CALL writeDensity(stars,vacuum,atoms,cell,sphhar,inp,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
     339           0 :                      0,-1.0,0.0,.FALSE.,den,'mdny')
     340             : 
     341             :    !---> save mz to file mdnz
     342           0 :    den%mt(:,0:,1:,1) = rho(:,0:,1:,4)
     343           0 :    den%pw(1:,1) = qpw(1:,4)
     344           0 :    den%vacz(1:,1:,1) = rht(1:,1:,4)
     345           0 :    den%vacxy(1:,1:,1:,1) = rhtxy(1:,1:,1:,4)
     346             :    CALL writeDensity(stars,vacuum,atoms,cell,sphhar,inp,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
     347           0 :                      0,-1.0,0.0,.FALSE.,den,'mdnz')
     348             : 
     349           0 :    DEALLOCATE (qpw,rhtxy,cdom,cdomvz,cdomvxy,ris,fftwork,rvacxy,rho,rht)
     350             : 
     351           0 : END SUBROUTINE pldngen
     352             : 
     353             : END MODULE m_pldngen

Generated by: LCOV version 1.13