LCOV - code coverage report
Current view: top level - optional - plotdop.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 126 216 58.3 %
Date: 2019-09-08 04:53:50 Functions: 2 2 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_plotdop
       8             : 
       9             : use m_juDFT
      10             : use m_types
      11             : 
      12             : !+++++++++++++++++++++++++++++++++++++++++++++++++
      13             : !  plot the charge density for fleur  code output
      14             : !     
      15             : !  if twodim = .false. a 3-D plot with nplot layers in z-direction
      16             : !  is constructed; the 3x3 matrix gives the 3 vectors of the cell ..
      17             : !  .gustav
      18             : !
      19             : !  Changed the input/output for easier use. 
      20             : !  This subroutine uses the file plot_inp for input. 
      21             : !  The old plotin-file is still supported by the old subroutine at the 
      22             : !  end of the module
      23             : !                     Juelich, 21.1.06 DW
      24             : !
      25             : !     +++++++++++++++++++++++++++++++++++++++++++++++++
      26             : 
      27             : CONTAINS
      28             : 
      29           2 : SUBROUTINE plotdop(oneD,dimension,stars,vacuum,sphhar,atoms,&
      30             :                    input,sym,cell,sliceplot,noco,cdnfname)
      31             : 
      32             :    USE m_outcdn
      33             :    USE m_loddop
      34             :    USE m_xsf_io
      35             :    USE m_cdn_io
      36             :    USE m_constants
      37             : 
      38             :    IMPLICIT NONE
      39             : 
      40             :    TYPE(t_oneD),                INTENT(IN)    :: oneD
      41             :    TYPE(t_dimension),           INTENT(IN)    :: dimension
      42             :    TYPE(t_stars),               INTENT(IN)    :: stars
      43             :    TYPE(t_vacuum),              INTENT(IN)    :: vacuum
      44             :    TYPE(t_sphhar),              INTENT(IN)    :: sphhar
      45             :    TYPE(t_atoms),               INTENT(IN)    :: atoms
      46             :    TYPE(t_input),               INTENT(IN)    :: input
      47             :    TYPE(t_sym),                 INTENT(IN)    :: sym
      48             :    TYPE(t_cell),                INTENT(IN)    :: cell
      49             :    TYPE(t_sliceplot),           INTENT(IN)    :: sliceplot
      50             :    TYPE(t_noco),                INTENT(IN)    :: noco
      51             :    CHARACTER(len=10), OPTIONAL, INTENT(IN)    :: cdnfname
      52             : 
      53             : !  .. Local Scalars ..
      54             :    REAL          :: tec,qint,fermiEnergyTemp,phi0,angss
      55             :    INTEGER       :: i,j,ix,iy,iz,jsp,na,nplo,iv,iflag,nfile
      56             :    INTEGER       :: nplot,nt,jm,jspin,numInFiles,numOutFiles
      57             :    LOGICAL       :: twodim,oldform,newform,l_qfix
      58             :    LOGICAL       :: cartesian,xsf,unwind,polar
      59             : 
      60             : !  .. Local Arrays ..
      61           2 :    TYPE(t_potden), ALLOCATABLE :: den(:)
      62           2 :    REAL, ALLOCATABLE    :: xdnout(:)
      63             :    REAL    :: pt(3),vec1(3),vec2(3),vec3(3),zero(3),help(3),qssc(3)
      64             :    INTEGER :: grid(3)
      65           4 :    REAL    :: rhocc(atoms%jmtd)
      66             :    REAL    :: point(3)
      67           2 :    CHARACTER (len=10), ALLOCATABLE :: cdnFilenames(:)
      68           2 :    CHARACTER (len=15), ALLOCATABLE :: outFilenames(:)
      69             :    CHARACTER (len=30)              :: filename
      70             :    CHARACTER (len=7)               :: textline
      71             : 
      72             :    REAL, PARAMETER :: eps = 1.0e-15
      73             : 
      74             :    NAMELIST /plot/twodim,cartesian,unwind,vec1,vec2,vec3,grid,zero,phi0,filename
      75             : 
      76           2 :    oldform = .false.
      77           2 :    INQUIRE(file ="plotin",exist = oldform) 
      78           2 :    IF (oldform) THEN 
      79           0 :       CALL juDFT_error("Use of plotin file no longer supported",calledby = "plotdop")
      80             :    END IF
      81             : 
      82           2 :    INQUIRE(file ="plot_inp",exist= newform)
      83           2 :    IF (.NOT.newform) THEN !no input file exists, create a template and exit
      84           0 :       OPEN(20,file ="plot_inp")
      85           0 :       WRITE(20,'(i2,a5,l1)') 2,",xsf=",.true.
      86           0 :       WRITE(20,*) "&PLOT twodim=t,cartesian=t"
      87           0 :       WRITE(20,*) "  vec1(1)=10.0 vec2(2)=10.0"
      88           0 :       WRITE(20,*) "  filename='plot1' /"
      89           0 :       WRITE(20,*) "&PLOT twodim=f,cartesian=f"
      90           0 :       WRITE(20,*) "  vec1(1)=1.0 vec1(2)=0.0 vec1(3)=0.0 "
      91           0 :       WRITE(20,*) "  vec2(1)=0.0 vec2(2)=1.0 vec2(3)=0.0 "
      92           0 :       WRITE(20,*) "  vec3(1)=0.0 vec3(2)=0.0 vec3(3)=1.0 "
      93           0 :       WRITE(20,*) "  grid(1)=30  grid(2)=30  grid(3)=30  "
      94           0 :       WRITE(20,*) "  zero(1)=0.0 zero(2)=0.0 zero(3)=0.5 "
      95           0 :       WRITE(20,*) "  filename ='plot2' /"
      96           0 :       CLOSE(20)
      97           0 :       WRITE(*,*) "No plot_inp file found. Created a template"
      98           0 :       CALL juDFT_error("Missing input for plot; modify plot_inp",calledby ="plotdop")
      99             :    END IF
     100             : 
     101           2 :    nfile = 120
     102           2 :    numInFiles = 0
     103           2 :    numOutFiles = 0
     104           2 :    IF(PRESENT(cdnfname)) THEN
     105             :       numInFiles = 1
     106             :       numOutFiles = 1
     107             :    ELSE
     108           2 :       IF(noco%l_noco) THEN
     109             :          numInFiles = 4
     110             :          numOutFiles = 4
     111             :       ELSE
     112           2 :          numInFiles = 1
     113           2 :          numOutFiles = 1
     114             :       END IF
     115             :    END IF
     116           2 :    ALLOCATE(den(numInFiles))
     117           2 :    ALLOCATE(cdnFilenames(numInFiles))
     118           2 :    IF(PRESENT(cdnfname)) THEN
     119           0 :       cdnFilenames(1) = cdnfname
     120             :    ELSE
     121           2 :       IF(noco%l_noco) THEN
     122           0 :         cdnFilenames(1)='cdn'
     123           0 :         cdnFilenames(2)='mdnx'
     124           0 :         cdnFilenames(3)='mdny'
     125           0 :         cdnFilenames(4)='mdnz'
     126             :       ELSE
     127           2 :          IF (sliceplot%slice) THEN
     128           1 :             cdnFilenames(1)='cdn_slice'
     129             :          ELSE
     130           1 :             cdnFilenames(1)='cdn1'
     131             :          END IF
     132             :       END IF
     133             :    END IF
     134             : 
     135             :    ! Read in charge/potential
     136           6 :    DO i = 1, numInFiles
     137           2 :       CALL den(i)%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
     138           2 :       IF(TRIM(ADJUSTL(cdnFilenames(i))).EQ.'cdn1') THEN
     139             :          CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,&
     140           1 :                           CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,den(i))
     141           1 :       ELSE IF(TRIM(ADJUSTL(cdnFilenames(i))).EQ.'cdn') THEN
     142             :          CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,&
     143           0 :                           CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,den(i))
     144             :       ELSE
     145             :          CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,&
     146           1 :                           CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,den(i),TRIM(ADJUSTL(cdnFilenames(i))))
     147             :       END IF
     148             : 
     149             :       ! Subtract core charge if input%score is set
     150           4 :       IF ((.NOT.noco%l_noco).AND.(input%score)) THEN
     151           0 :          OPEN (17,file='cdnc',form='unformatted',status='old')
     152           0 :          REWIND 17
     153           0 :          DO jspin = 1, input%jspins
     154           0 :             DO nt = 1, atoms%ntype
     155           0 :                jm = atoms%jri(nt)
     156           0 :                READ (17) (rhocc(j),j=1,jm)
     157           0 :                DO j = 1, atoms%jri(nt)
     158           0 :                   den(i)%mt(j,0,nt,jspin) = den(i)%mt(j,0,nt,jspin) - rhocc(j)/2.0/SQRT(pi_const)
     159             :                END DO
     160           0 :                READ (17) tec
     161             :             END DO
     162           0 :             READ (17) qint
     163           0 :             den(i)%pw(1,jspin) = den(i)%pw(1,jspin) - qint/cell%volint
     164             :          END DO
     165           0 :          CLOSE (17)
     166           2 :       ELSE IF (input%score) THEN
     167           0 :          CALL juDFT_error('Subtracting core charge in noco calculations not supported', calledby = 'plotdop')
     168             :       END IF
     169             :    END DO
     170             : 
     171           2 :    IF (noco%l_ss) THEN 
     172           0 :       qssc = MATMUL(TRANSPOSE(cell%bmat),noco%qss) 
     173             :    END IF 
     174             : 
     175             :    ! Open the plot_inp file for input
     176           2 :    OPEN (18,file='plot_inp')
     177           2 :    READ(18,'(i2,5x,l1,1x,a)') nplot,xsf,textline
     178           2 :    polar = .FALSE.
     179           2 :    IF ((noco%l_noco).AND.(numInFiles.EQ.4)) THEN
     180           0 :       polar = (textline(1:7)=='polar=T').OR.(textline(1:7)=='polar=t')
     181             :       IF (polar) THEN
     182             :          numOutFiles = 7
     183             :       END IF
     184             :    END IF
     185           2 :    ALLOCATE(outFilenames(numOutFiles))
     186           2 :    ALLOCATE(xdnout(numOutFiles))
     187           2 :    IF(numOutFiles.EQ.1) THEN
     188           2 :       outFilenames(1) = 'plot'
     189             :    ELSE
     190           0 :       DO i = 1, numInFiles
     191           0 :          outFilenames(i) = TRIM(ADJUSTL(cdnFilenames(i)))//'_pl'
     192             :       END DO
     193           0 :       IF (polar) THEN
     194           0 :          outFilenames(5) = 'mabs_pl'
     195           0 :          outFilenames(6) = 'mtha_pl'
     196           0 :          outFilenames(7) = 'mphi_pl'
     197             :       END IF
     198             :    END IF
     199             : 
     200             :    ! If xsf is specified we create input files for xcrysden
     201           2 :    IF (xsf) THEN
     202           6 :       DO i = 1, numOutFiles
     203           2 :          OPEN(nfile+i,file=TRIM(ADJUSTL(outFilenames(i)))//'.xsf',form='formatted')
     204           4 :          CALL xsf_WRITE_atoms(nfile+i,atoms,input%film,oneD%odi%d1,cell%amat)
     205             :       END DO
     206             :    END IF
     207             : 
     208             :    ! Loop over all plots
     209           4 :    DO nplo = 1, nplot
     210             : 
     211             :       ! the defaults
     212           2 :       twodim = .TRUE.
     213           2 :       cartesian = .TRUE.
     214           2 :       grid = (/100,100,100/)
     215           2 :       vec1 = (/0.,0.,0./)
     216           2 :       vec2 = (/0.,0.,0./)
     217           2 :       vec3 = (/0.,0.,0./)
     218           2 :       zero = (/0.,0.,0./)
     219           2 :       filename = "default"
     220           2 :       READ(18,plot)
     221           6 :       IF (twodim.AND.ANY(grid(1:2)<1)) &
     222           0 :          CALL juDFT_error("Illegal grid size in plot",calledby="plotdop")
     223           8 :       IF (.NOT.twodim.AND.ANY(grid<1)) &
     224           0 :          CALL juDFT_error("Illegal grid size in plot",calledby="plotdop")
     225           2 :       IF (twodim) grid(3) = 1
     226             : 
     227             :       !calculate cartesian coordinates if needed
     228           2 :       IF (.NOT.cartesian) THEN
     229           2 :          vec1=matmul(cell%amat,vec1)
     230           2 :          vec2=matmul(cell%amat,vec2)
     231           2 :          vec3=matmul(cell%amat,vec3)
     232           2 :          zero=matmul(cell%amat,zero)
     233             :       END IF
     234             : 
     235             :       !Open the file
     236           2 :       IF (filename =="default") WRITE(filename,'(a,i2)') "plot",nplo
     237           6 :       DO i = 1, numOutFiles
     238           4 :          IF (xsf) THEN
     239           2 :             CALL xsf_WRITE_header(nfile+i,twodim,filename,vec1,vec2,vec3,zero,grid)
     240             :          ELSE
     241           0 :             IF (numOutFiles.NE.1) THEN
     242           0 :                OPEN (nfile+i,file = filename//outFilenames(i),form='formatted')
     243             :             ELSE
     244           0 :                OPEN (nfile+i,file = filename,form='formatted')
     245             :             END IF
     246             :          END IF
     247             :       END DO
     248             : 
     249             :       !loop over spins
     250           4 :       DO jsp = 1, input%jspins
     251             :          !loop over all points
     252          62 :          DO iz = 0, grid(3)-1
     253         662 :             DO iy = 0, grid(2)-1
     254        6660 :                DO ix = 0, grid(1)-1
     255             : 
     256             :                   point = zero + vec1*REAL(ix)/(grid(1)-1) +&
     257        6000 :                                  vec2*REAL(iy)/(grid(2)-1)
     258        6000 :                   IF (.NOT.twodim) point = point + vec3*REAL(iz)/(grid(3)-1)
     259             : 
     260             :                   ! Set region specific parameters for point
     261             :                   
     262             :                   ! Get MT sphere for point if point is in MT sphere
     263        6000 :                   CALL getMTSphere(input,cell,atoms,oneD,point,nt,na,pt)
     264        6000 :                   IF (na.NE.0) THEN
     265             :                      ! In MT sphere
     266        1312 :                      iv = 0
     267        1312 :                      iflag = 1
     268        4688 :                   ELSE IF (input%film.AND..NOT.oneD%odi%d1.AND.ABS(point(3))>=cell%z1) THEN
     269             :                      ! In vacuum in 2D system
     270        1200 :                      iv = 1
     271        1200 :                      iflag = 0
     272        1200 :                      pt(:) = point(:)
     273        3488 :                   ELSE IF ((oneD%odi%d1).AND.(SQRT((point(1))**2 + (point(2))**2)>=cell%z1)) THEN
     274             :                      ! In vacuum in 1D system
     275           0 :                      iv = 1
     276           0 :                      iflag = 0
     277           0 :                      pt(:) = point(:)
     278             :                   ELSE
     279             :                      ! In interstitial region
     280        3488 :                      iv = 0
     281        3488 :                      iflag = 2
     282        3488 :                      pt(:) = point(:)
     283             :                   END IF
     284             : 
     285       18000 :                   DO i = 1, numInFiles
     286             :                      CALL outcdn(pt,nt,na,iv,iflag,jsp,sliceplot,stars,&
     287             :                                  vacuum,sphhar,atoms,sym,cell,oneD,&
     288             :                                  den(i)%pw,den(i)%vacxy,den(i)%mt,&
     289       12000 :                                  den(i)%vacz,xdnout(i))
     290             :                   END DO
     291             : 
     292        6000 :                   IF (na.NE.0) THEN
     293        1312 :                      IF (noco%l_ss) THEN 
     294             :                         ! rotate magnetization "backward"
     295           0 :                         angss = DOT_PRODUCT(qssc,pt-atoms%pos(:,na))
     296           0 :                         help(1) = xdnout(2)
     297           0 :                         help(2) = xdnout(3)
     298           0 :                         xdnout(2) = +help(1)*COS(angss)+help(2)*SIN(angss) 
     299           0 :                         xdnout(3) = -help(1)*SIN(angss)+help(2)*COS(angss) 
     300             :                         ! xdnout(2)=0. ; xdnout(3)=0. ; xdnout(4)=0. 
     301             :                      END IF
     302             :                   END IF
     303             : 
     304        6000 :                   IF (noco%l_ss .AND. (.NOT. unwind)) THEN
     305             :                      ! rotate magnetization
     306           0 :                      angss = DOT_PRODUCT(qssc,point)
     307           0 :                      help(1) = xdnout(2)
     308           0 :                      help(2) = xdnout(3)
     309           0 :                      xdnout(2) = +help(1)*COS(angss) -help(2)*SIN(angss)
     310           0 :                      xdnout(3) = +help(1)*SIN(angss) +help(2)*COS(angss)
     311             :                   END IF
     312             : 
     313        6000 :                   IF (polar) THEN
     314           0 :                      xdnout(5) = SQRT(ABS(xdnout(2)**2+xdnout(3)**2+xdnout(4)**2))
     315           0 :                      IF (xdnout(5)<eps) THEN
     316           0 :                         xdnout(5)= 0.0
     317           0 :                         xdnout(6)= -tpi_const
     318           0 :                         xdnout(7)= -tpi_const
     319             :                      ELSE
     320           0 :                         DO j = 1, 3
     321           0 :                            help(j) = xdnout(1+j)/xdnout(5) 
     322             :                         END DO
     323           0 :                         IF (help(3)<0.5) THEN
     324           0 :                            xdnout(6)= ACOS(help(3))
     325             :                         ELSE
     326           0 :                            xdnout(6)= pi_const/2.0-ASIN(help(3))
     327             :                         END IF
     328           0 :                         IF (SQRT(ABS(help(1)**2+help(2)**2)) < eps) THEN
     329           0 :                            xdnout(7)= -tpi_const
     330             :                         ELSE
     331           0 :                            IF ( ABS(help(1)) > ABS(help(2)) ) THEN
     332           0 :                               xdnout(7)= ABS(ATAN(help(2)/help(1)))
     333             :                            ELSE
     334           0 :                               xdnout(7)= pi_const/2.0-ABS(ATAN(help(1)/help(2)))
     335             :                            END IF
     336           0 :                            IF (help(2)<0.0) THEN
     337           0 :                               xdnout(7)= -xdnout(7)
     338             :                            END IF
     339           0 :                            IF (help(1)<0.0) THEN
     340           0 :                               xdnout(7)= pi_const-xdnout(7)
     341             :                            END IF
     342           0 :                            DO WHILE (xdnout(7)-pi_const*phi0 > +pi_const)
     343           0 :                               xdnout(7)= xdnout(7)-tpi_const
     344             :                            END DO
     345           0 :                            DO WHILE (xdnout(7)-pi_const*phi0 < -pi_const)
     346           0 :                               xdnout(7)= xdnout(7)+tpi_const
     347             :                            END DO
     348             :                         END IF
     349             :                      END IF
     350           0 :                      xdnout(6)= xdnout(6)/pi_const
     351           0 :                      xdnout(7)= xdnout(7)/pi_const
     352             :                   END IF ! (polar)
     353             : 
     354       18600 :                   DO i = 1, numOutFiles
     355       12000 :                      IF (xsf) THEN
     356        6000 :                         WRITE(nfile+i,*) xdnout(i)
     357             :                      ELSE
     358           0 :                         WRITE(nfile+i,'(4e15.7)') point ,xdnout(i)
     359             :                      END IF
     360             :                   END DO
     361             : 
     362             :                END DO
     363             :             END DO
     364             :          END DO !z-loop
     365           8 :          DO i = 1, numOutFiles
     366           2 :             IF (xsf.AND.jsp /= input%jspins) &
     367           2 :                CALL xsf_WRITE_newblock(nfile+i,twodim,vec1,vec2,vec3,zero,grid)
     368             :          END DO
     369             :       END DO !Spin-loop
     370             : 
     371           8 :       DO i = 1, numOutFiles
     372           4 :          IF (xsf) THEN
     373           2 :             CALL xsf_WRITE_endblock(nfile+i,twodim)
     374             :          ELSE
     375           0 :             CLOSE(nfile+i)
     376             :          END IF
     377             :       END DO
     378             :    END DO !nplot  
     379             :     
     380           2 :    CLOSE(18)
     381           2 :    IF (xsf) THEN
     382           6 :       DO i = 1, numOutFiles
     383           4 :          CLOSE(nfile+i)
     384             :       END DO
     385             :    END IF
     386             : 
     387           2 :    DEALLOCATE(xdnout, cdnFilenames, outFilenames)
     388             : 
     389           2 :    END SUBROUTINE plotdop
     390             : 
     391             : 
     392             : 
     393        6000 :    SUBROUTINE getMTSphere(input,cell,atoms,oneD,point,iType,iAtom,pt)
     394             : 
     395             :    IMPLICIT NONE
     396             : 
     397             :    TYPE(t_input), INTENT(IN)    :: input
     398             :    TYPE(t_cell),  INTENT(IN)    :: cell
     399             :    TYPE(t_atoms), INTENT(IN)    :: atoms
     400             :    TYPE(t_oneD),  INTENT(IN)    :: oneD
     401             : 
     402             :    INTEGER,       INTENT(OUT)   :: iType, iAtom
     403             :    REAL,          INTENT(OUT)   :: pt(3)
     404             :    REAL,          INTENT(IN)    :: point(3)
     405             : 
     406             :    INTEGER                      :: ii1, ii2, ii3, i1, i2, i3, nq
     407             :    REAL                         :: s
     408             : 
     409        6000 :    ii1 = 3
     410        6000 :    ii2 = 3
     411        6000 :    ii3 = 3
     412        6000 :    IF (input%film .AND. .NOT.oneD%odi%d1) ii3 = 0
     413        6000 :    IF (oneD%odi%d1) THEN
     414           0 :       ii1 = 0
     415           0 :       ii2 = 0
     416             :    END IF
     417             : 
     418       42096 :    DO i1 = -ii1, ii1
     419      554000 :       DO i2 = -ii2, ii2
     420      549312 :          DO i3 = -ii3, ii3
     421      257264 :             pt = point+MATMUL(cell%amat,(/i1,i2,i3/))
     422      257264 :             iAtom = 0
     423     1282152 :             DO iType = 1, atoms%ntype
     424     2564728 :                DO nq = 1, atoms%neq(iType)
     425     1539840 :                   iAtom = iAtom + 1
     426     1539840 :                   s = SQRT(DOT_PRODUCT(atoms%pos(:,iAtom)-pt,atoms%pos(:,iAtom)-pt))
     427     2308776 :                   IF (s<atoms%rmsh(atoms%jri(iType),iType)) THEN
     428             :                      ! Return with the current iType, iAtom, pt
     429             :                      RETURN
     430             :                   END IF
     431             :                END DO
     432             :             END DO
     433             :          END DO
     434             :       END DO
     435             :    END DO !i1
     436             : 
     437             :    ! If no MT sphere encloses the point return 0 for iType, iAtom
     438        4688 :    iType = 0
     439        4688 :    iAtom = 0
     440        4688 :    pt(:) = point(:)
     441             : 
     442             :    END SUBROUTINE getMTSphere
     443             : 
     444             : END MODULE m_plotdop

Generated by: LCOV version 1.13