LCOV - code coverage report
Current view: top level - types - types_enpara.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 244 294 83.0 %
Date: 2019-09-08 04:53:50 Functions: 8 10 80.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_types_enpara
       8             :   USE m_judft
       9             :   IMPLICIT NONE
      10             :   PRIVATE
      11             :   TYPE t_enpara
      12             :      REAL, ALLOCATABLE CPP_MANAGED   :: el0(:,:,:)
      13             :      REAL                 :: evac0(2,2)
      14             :      REAL                 :: evac(2,2)
      15             :      REAL, ALLOCATABLE    :: ello0(:,:,:)
      16             :      REAL, ALLOCATABLE    :: el1(:,:,:)
      17             :      REAL                 :: evac1(2,2)
      18             :      REAL, ALLOCATABLE    :: ello1(:,:,:)
      19             :      REAL, ALLOCATABLE    :: enmix(:)
      20             :      INTEGER, ALLOCATABLE :: skiplo(:,:)
      21             :      LOGICAL, ALLOCATABLE :: lchange(:,:,:)
      22             :      LOGICAL, ALLOCATABLE :: lchg_v(:,:)
      23             :      LOGICAL, ALLOCATABLE :: llochg(:,:,:)
      24             :      REAL                 :: epara_min
      25             :      LOGICAL              :: ready ! are the enpara's ok for calculation?
      26             :      LOGICAL              :: floating !floating energy parameters are relative to potential
      27             :      INTEGER,ALLOCATABLE  :: qn_el(:,:,:)    !>if these are .ne.0 they are understood as
      28             :      INTEGER,ALLOCATABLE  :: qn_ello(:,:,:)  !>quantum numbers
      29             :    CONTAINS
      30             :      PROCEDURE :: init
      31             :      PROCEDURE :: update
      32             :      PROCEDURE :: read
      33             :      PROCEDURE :: write
      34             :      PROCEDURE :: mix
      35             :      PROCEDURE :: calcOutParams
      36             :   END TYPE t_enpara
      37             : 
      38             : 
      39             : 
      40             :   PUBLIC:: t_enpara
      41             : 
      42             : CONTAINS
      43          79 :   SUBROUTINE init(this,atoms,jspins,l_defaults)
      44             :     USE m_types_setup
      45             :     USE m_constants
      46             :     CLASS(t_enpara),INTENT(inout):: this
      47             :     TYPE(t_atoms),INTENT(IN)     :: atoms
      48             :     INTEGER,INTENT(IN)           :: jspins
      49             :     LOGICAL,INTENT(IN),OPTIONAL  :: l_defaults
      50             : 
      51             :     INTEGER :: n,i,jsp,l
      52             : 
      53          79 :     ALLOCATE(this%el0(0:atoms%lmaxd,atoms%ntype,jspins),this%el1(0:atoms%lmaxd,atoms%ntype,jspins))
      54          79 :     ALLOCATE(this%ello0(atoms%nlod,atoms%ntype,jspins),this%ello1(atoms%nlod,atoms%ntype,jspins))
      55         188 :     this%el0=-1E99
      56         188 :     this%ello0=-1E99
      57         237 :     this%evac0=-1E99
      58             : 
      59             : 
      60          79 :     ALLOCATE(this%llochg(atoms%nlod,atoms%ntype,jspins))
      61          79 :     ALLOCATE(this%lchg_v(2,jspins))
      62          79 :     ALLOCATE(this%skiplo(atoms%ntype,jspins))
      63          79 :     ALLOCATE(this%lchange(0:atoms%lmaxd,atoms%ntype,jspins))
      64         188 :     this%llochg=.FALSE.;this%lchg_v=.FALSE.;this%lchange=.FALSE.
      65         188 :     this%skiplo=0
      66          79 :     ALLOCATE(this%enmix(jspins))
      67         188 :     this%enmix=0.0
      68             : 
      69          79 :     this%ready=.FALSE.
      70          79 :     this%floating=.FALSE.
      71             : 
      72          79 :     ALLOCATE(this%qn_el(0:3,atoms%ntype,jspins))
      73          79 :     ALLOCATE(this%qn_ello(atoms%nlod,atoms%ntype,jspins))
      74             : 
      75          79 :     IF (PRESENT(l_defaults)) THEN
      76          38 :        IF (.NOT.l_defaults) RETURN
      77             :     ENDIF
      78             :     !Set most simple defaults
      79         153 :     DO jsp=1,jspins
      80         215 :        DO n = 1,atoms%ntype
      81         118 :           IF ( atoms%nz(n) < 3 ) THEN
      82           0 :              this%qn_el(0:3,n,jsp) =  (/1,2,3,4/) 
      83         118 :           ELSEIF ( atoms%nz(n) < 11 ) THEN
      84          10 :              this%qn_el(0:3,n,jsp) =  (/2,2,3,4/) 
      85         108 :           ELSEIF ( atoms%nz(n) < 19 ) THEN
      86          16 :              this%qn_el(0:3,n,jsp) =  (/3,3,3,4/) 
      87          92 :           ELSEIF ( atoms%nz(n) < 31 ) THEN
      88          81 :              this%qn_el(0:3,n,jsp) =  (/4,4,3,4/) 
      89          11 :           ELSEIF ( atoms%nz(n) < 37 ) THEN
      90           6 :              this%qn_el(0:3,n,jsp) =  (/4,4,4,4/) 
      91           5 :           ELSEIF ( atoms%nz(n) < 49 ) THEN
      92           0 :              this%qn_el(0:3,n,jsp) =  (/5,5,4,4/) 
      93           5 :           ELSEIF ( atoms%nz(n) < 55 ) THEN
      94           0 :              this%qn_el(0:3,n,jsp) =  (/5,5,5,4/) 
      95           5 :           ELSEIF ( atoms%nz(n) < 72 ) THEN
      96           0 :              this%qn_el(0:3,n,jsp) =  (/6,6,5,4/) 
      97           5 :           ELSEIF ( atoms%nz(n) < 81 ) THEN
      98           0 :              this%qn_el(0:3,n,jsp) =  (/6,6,5,5/) 
      99           5 :           ELSEIF ( atoms%nz(n) < 87 ) THEN
     100           5 :              this%qn_el(0:3,n,jsp) =  (/6,6,6,5/) 
     101             :           ELSE
     102           0 :              this%qn_el(0:3,n,jsp) =  (/7,7,6,5/) 
     103             :           ENDIF
     104             :           
     105         212 :           DO i = 1, atoms%nlo(n)
     106         156 :              IF (atoms%llo(i,n)<0) THEN
     107             :                 !llo might not be initialized
     108             :                 !in this case defaults broken
     109          30 :                 this%qn_ello(i,n,jsp) = 0
     110          30 :                 this%skiplo(n,jsp) = 0
     111             :              ELSE
     112           8 :                 this%qn_ello(i,n,jsp) = this%qn_el(atoms%llo(i,n),n,jsp) - 1
     113           8 :                 this%skiplo(n,jsp) = this%skiplo(n,jsp) + (2*atoms%llo(i,n)+1)
     114             :              ENDIF
     115             :           ENDDO
     116             :        ENDDO
     117             :     ENDDO
     118             : 
     119          41 :     this%evac0=eVac0Default_const
     120             : 
     121             :   END SUBROUTINE init
     122             : 
     123             :   !> This subroutine adjusts the energy parameters to the potential. In particular, it
     124             :   !! calculated them in case of qn_el>-1,qn_ello>-1
     125             :   !! Before this was done in lodpot.F
     126         340 :   SUBROUTINE update(enpara,mpi,atoms,vacuum,input,v)
     127             :     USE m_types_setup
     128             :     USE m_types_mpi
     129             :     USE m_xmlOutput
     130             :     USE m_types_potden
     131             :     USE m_find_enpara
     132             :     CLASS(t_enpara),INTENT(inout):: enpara
     133             :     TYPE(t_mpi),INTENT(IN)      :: mpi
     134             :     TYPE(t_atoms),INTENT(IN)    :: atoms
     135             :     TYPE(t_input),INTENT(IN)    :: input
     136             :     TYPE(t_vacuum),INTENT(IN)   :: vacuum
     137             :     TYPE(t_potden),INTENT(IN)   :: v
     138             : 
     139             : 
     140             :     LOGICAL ::  l_enpara
     141         680 :     LOGICAL ::  l_done(0:atoms%lmaxd,atoms%ntype,input%jspins)
     142         680 :     LOGICAL ::  lo_done(atoms%nlod,atoms%ntype,input%jspins)
     143             :     REAL    ::  vbar,vz0,rj
     144             :     INTEGER ::  n,jsp,l,ilo,j,ivac
     145             :     CHARACTER(LEN=20)    :: attributes(5)
     146         680 :     REAL ::  e_lo(0:3,atoms%ntype)!Store info on branches to do IO after OMP
     147         680 :     REAL ::  e_up(0:3,atoms%ntype)
     148         680 :     REAL ::  elo_lo(atoms%nlod,atoms%ntype)
     149         680 :     REAL ::  elo_up(atoms%nlod,atoms%ntype)
     150             : 
     151         340 :     IF (mpi%irank  == 0) CALL openXMLElement('energyParameters',(/'units'/),(/'Htr'/))
     152             : 
     153         948 :     l_done = .FALSE.;lo_done=.FALSE.
     154         948 :     DO jsp = 1,input%jspins
     155             :        !$OMP PARALLEL DO DEFAULT(none) &
     156             :        !$OMP SHARED(atoms,enpara,jsp,l_done,mpi,v,lo_done,e_lo,e_up,elo_lo,elo_up) &
     157             :        !$OMP PRIVATE(n,l,ilo)
     158             :        !! First calculate energy paramter from quantum numbers if these are given...
     159             :        !! l_done stores the index of those energy parameter updated
     160             :        DO n = 1, atoms%ntype
     161        5820 :           DO l = 0,3
     162        5820 :              IF( enpara%qn_el(l,n,jsp).ne.0)THEN 
     163        3808 :                 l_done(l,n,jsp) = .TRUE.
     164        3808 :                 enpara%el0(l,n,jsp)=find_enpara(.FALSE.,l,n,jsp,enpara%qn_el(l,n,jsp),atoms,mpi,v%mt(:,0,n,jsp),e_lo(l,n),e_up(l,n))
     165        3808 :                 IF( l .EQ. 3 ) THEN
     166        5816 :                    enpara%el0(4:,n,jsp) = enpara%el0(3,n,jsp)
     167         952 :                    l_done(4:,n,jsp) = .TRUE.
     168             :                 END IF
     169             :              ELSE 
     170         848 :                 l_done(l,n,jsp) = .FALSE.
     171             :              END IF
     172             :           ENDDO ! l
     173             :           ! Now for the lo's
     174        2576 :           DO ilo = 1, atoms%nlo(n)
     175         196 :              l = atoms%llo(ilo,n)
     176        1360 :              IF( enpara%qn_ello(ilo,n,jsp).NE.0) THEN
     177         188 :                 lo_done(ilo,n,jsp) = .TRUE.
     178         188 :                 enpara%ello0(ilo,n,jsp)=find_enpara(.TRUE.,l,n,jsp,enpara%qn_ello(ilo,n,jsp),atoms,mpi,v%mt(:,0,n,jsp),elo_lo(ilo,n),elo_up(ilo,n))
     179             :              ELSE
     180           8 :                 lo_done(ilo,n,jsp) = .FALSE.
     181             :              ENDIF
     182             :           ENDDO
     183             :        ENDDO ! n
     184             :        !$OMP END PARALLEL DO
     185         608 :        IF (mpi%irank==0) THEN
     186         304 :           WRITE(6,*)
     187         304 :           WRITE(6,*) "Updated energy parameters for spin:",jsp
     188             :           !Same loop for IO
     189         886 :           DO n = 1, atoms%ntype
     190        2910 :              DO l = 0,3
     191        2910 :                 IF( l_done(l,n,jsp)) CALL priv_write(.FALSE.,l,n,jsp,enpara%qn_el(l,n,jsp),e_lo(l,n),e_up(l,n),enpara%el0(l,n,jsp))
     192             :              ENDDO ! l
     193             :              ! Now for the lo's
     194         984 :              DO ilo = 1, atoms%nlo(n)
     195          98 :                 l = atoms%llo(ilo,n)
     196         680 :                 IF( lo_done(ilo,n,jsp)) CALL priv_write(.TRUE.,l,n,jsp,enpara%qn_ello(ilo,n,jsp),elo_lo(ilo,n),elo_up(ilo,n),enpara%ello0(ilo,n,jsp))
     197             :              END DO
     198             :           END DO
     199             :        ENDIF
     200             : 
     201             :        !!   Now check for floating energy parameters (not for those with l_done=T)
     202         608 :        IF (enpara%floating) THEN
     203           0 :           types_loop: DO n = 1,atoms%ntype 
     204             :              !
     205             :              !--->    determine energy parameters if lepr=1. the reference energy
     206             :              !--->    is the value of the l=0 potential at approximately rmt/4.
     207             :              !
     208           0 :              j = atoms%jri(n) - (LOG(4.0)/atoms%dx(n)+1.51)
     209           0 :              rj = atoms%rmt(n)*EXP(atoms%dx(n)* (j-atoms%jri(n)))
     210           0 :              vbar = v%mt(j,0,n,jsp)/rj
     211           0 :              IF (mpi%irank.EQ.0) THEN
     212           0 :                 WRITE ( 6,'('' spin'',i2,'', atom type'',i3,'' ='',f12.6,''   r='',f8.5)') jsp,n,vbar,rj
     213             :              ENDIF
     214           0 :              DO l = 0,atoms%lmax(n)
     215           0 :                 IF ( .NOT.l_done(l,n,jsp) ) THEN
     216           0 :                    enpara%el0(l,n,jsp) = vbar + enpara%el0(l,n,jsp)
     217             :                 END IF
     218             :              ENDDO
     219           0 :              IF (atoms%nlo(n).GE.1) THEN
     220           0 :                 DO ilo = 1,atoms%nlo(n)
     221           0 :                    IF ( .NOT. lo_done(ilo,n,jsp) ) THEN
     222           0 :                       enpara%ello0(ilo,n,jsp) = vbar + enpara%ello0(ilo,n,jsp)
     223             :                       !+apw+lo
     224           0 :                       IF (atoms%l_dulo(ilo,n)) THEN
     225           0 :                          enpara%ello0(ilo,n,jsp) = enpara%el0(atoms%llo(ilo,n),n,jsp)
     226             :                       ENDIF
     227             :                       !-apw+lo
     228             :                    END IF
     229             :                 END DO
     230             :              ENDIF
     231             :           END DO types_loop
     232             :        ENDIF
     233         948 :        IF (input%film) THEN
     234             : 
     235          16 :           INQUIRE (file ='enpara',exist= l_enpara)
     236          16 :           IF(l_enpara) enpara%evac0(:,jsp) = enpara%evac(:,jsp)
     237             : 
     238             :           !--->    vacuum energy parameters: for floating: relative to potential
     239             :           !--->    at vacuum-interstitial interface (better for electric field)
     240             : 
     241          32 :           DO ivac = 1,vacuum%nvac
     242          16 :              vz0 = 0.0
     243          16 :              IF (enpara%floating) THEN
     244           0 :                 vz0 = v%vacz(1,ivac,jsp)
     245           0 :                 IF (mpi%irank.EQ.0) THEN
     246           0 :                    WRITE ( 6,'('' spin'',i2,'', vacuum   '',i3,'' ='',f12.6)') jsp,ivac,vz0 
     247             :                 ENDIF
     248             :              ENDIF
     249          16 :              enpara%evac(ivac,jsp) = enpara%evac0(ivac,jsp) + vz0
     250          16 :              IF (.NOT.l_enpara) THEN
     251           8 :                 enpara%evac(ivac,jsp) = v%vacz(vacuum%nmz,ivac,jsp) + enpara%evac0(ivac,jsp)
     252             :              END IF
     253          32 :              IF (mpi%irank.EQ.0) THEN
     254          48 :                 attributes = ''
     255           8 :                 WRITE(attributes(1),'(i0)') ivac
     256           8 :                 WRITE(attributes(2),'(i0)') jsp
     257           8 :                 WRITE(attributes(3),'(f16.10)') v%vacz(1,ivac,jsp)
     258           8 :                 WRITE(attributes(4),'(f16.10)') v%vacz(vacuum%nmz,ivac,jsp)
     259           8 :                 WRITE(attributes(5),'(f16.10)') enpara%evac(ivac,jsp)
     260             :                 CALL writeXMLElementForm('vacuumEP',(/'vacuum','spin  ','vzIR  ','vzInf ','value '/),&
     261           8 :                      attributes(1:5),RESHAPE((/6+4,4,4,5,5+13,8,1,16,16,16/),(/5,2/)))
     262             :              END IF
     263             :           ENDDO
     264          16 :           IF (vacuum%nvac.EQ.1) THEN
     265          16 :              enpara%evac(2,jsp) = enpara%evac(1,jsp)
     266             :           END IF
     267             :        END IF
     268             :     END DO
     269             : 
     270             : !    enpara%ready=(ALL(enpara%el0>-1E99).AND.ALL(enpara%ello0>-1E99))
     271         340 :     enpara%epara_min=MIN(MINVAL(enpara%el0),MINVAL(enpara%ello0))
     272             :     
     273         340 :     IF (mpi%irank  == 0) CALL closeXMLElement('energyParameters')
     274         340 :   END SUBROUTINE update
     275             : 
     276          22 :   SUBROUTINE READ(enpara,atoms,jspins,film,l_required)
     277             :     USE m_types_setup
     278             :     IMPLICIT NONE
     279             :     CLASS(t_enpara),INTENT(INOUT):: enpara
     280             :     INTEGER, INTENT (IN)        :: jspins
     281             :     TYPE(t_atoms),INTENT(IN)    :: atoms
     282             :     LOGICAL,INTENT(IN)          :: film,l_required
     283             : 
     284             :     INTEGER :: n,l,lo,skip_t,io_err,jsp
     285             :     logical :: l_exist
     286             : 
     287          22 :     INQUIRE(file="enpara",exist=l_exist)
     288          22 :     IF (.NOT.l_exist.AND.l_required) CALL juDFT_error("No enpara file found")
     289          44 :     IF (.NOT.l_exist) RETURN
     290             : 
     291          14 :     OPEN(40,file="enpara",form="formatted",status="old")
     292             : 
     293          32 :     DO jsp=1,jspins
     294             : 
     295             :        !-->  first line contains mixing parameter!
     296             : 
     297          18 :        READ (40,FMT ='(48x,f10.6)',END=200) enpara%enmix(jsp)
     298          18 :        READ (40,*)                       ! skip next line
     299          18 :        IF (enpara%enmix(jsp).EQ.0.0) enpara%enmix(jsp) = 1.0
     300          18 :        WRITE (6,FMT=8001) jsp
     301          18 :        WRITE (6,FMT=8000)
     302          18 :        skip_t = 0
     303          54 :        DO n = 1,atoms%ntype
     304         216 :           READ (40,FMT=8040,END=200) (enpara%el0(l,n,jsp),l=0,3),&
     305         396 :                (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp)    
     306         216 :           WRITE (6,FMT=8140) n,(enpara%el0(l,n,jsp),l=0,3),&
     307         396 :                (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp)    
     308             :           !
     309             :           !--->    energy parameters for the local orbitals
     310             :           !
     311          36 :           IF (atoms%nlo(n).GE.1) THEN
     312           8 :              skip_t = skip_t + enpara%skiplo(n,jsp) * atoms%neq(n)
     313           8 :              READ (40,FMT=8039,END=200)  (enpara%ello0(lo,n,jsp),lo=1,atoms%nlo(n))
     314          24 :              READ (40,FMT=8038,END=200) (enpara%llochg(lo,n,jsp),lo=1,atoms%nlo(n))
     315          24 :              WRITE (6,FMT=8139)          (enpara%ello0(lo,n,jsp),lo=1,atoms%nlo(n))
     316           8 :              WRITE (6,FMT=8138)         (enpara%llochg(lo,n,jsp),lo=1,atoms%nlo(n))
     317          28 :           ELSEIF (enpara%skiplo(n,jsp).GT.0) THEN
     318           0 :              WRITE (6,*) "for atom",n," no LO's were specified"
     319           0 :              WRITE (6,*) 'but skiplo was set to',enpara%skiplo 
     320             :              CALL juDFT_error("No LO's but skiplo",calledby ="enpara",&
     321           0 :                   hint="If no LO's are set skiplo must be 0 in enpara")
     322             :           END IF
     323             :           !Integer values mean we have to set the qn-arrays
     324          36 :           enpara%qn_el(:,n,jsp)=0
     325         324 :           DO l=0,3
     326         180 :              IF (enpara%el0(l,n,jsp)==NINT(enpara%el0(l,n,jsp))) enpara%qn_el(l,n,jsp)=NINT(enpara%el0(l,n,jsp))
     327             :           ENDDO
     328          36 :           enpara%qn_ello(:,n,jsp)=0
     329          44 :           DO l=1,atoms%nlo(n)
     330          44 :              IF (enpara%ello0(l,n,jsp)==NINT(enpara%ello0(l,n,jsp))) enpara%qn_ello(l,n,jsp)=NINT(enpara%ello0(l,n,jsp))
     331             :           ENDDO
     332             :           !
     333             :           !--->    set the energy parameters with l>3 to the value of l=3
     334             :           !
     335          54 :           enpara%el0(4:,n,jsp) = enpara%el0(3,n,jsp)
     336             :        ENDDO   ! atoms%ntype
     337             : 
     338          18 :        IF (film) THEN
     339           4 :           enpara%lchg_v = .TRUE.
     340           4 :           READ (40,FMT=8050,END=200) enpara%evac0(1,jsp),enpara%lchg_v(1,jsp),enpara%evac0(2,jsp)
     341           4 :           WRITE (6,FMT=8150)         enpara%evac0(1,jsp),enpara%lchg_v(1,jsp),enpara%evac0(2,jsp)
     342             :        ENDIF
     343          32 :        IF (atoms%nlod.GE.1) THEN               
     344          18 :           WRITE (6,FMT=8090) jsp,skip_t
     345          18 :           WRITE (6,FMT=8091) 
     346             :        END IF
     347             :     END DO
     348             : 
     349          14 :     enpara%evac(:,:) = enpara%evac0(:,:)
     350             : 
     351          14 :     CLOSE(40)
     352             :     ! input formats
     353             : 
     354             : 8038 FORMAT (14x,60(l1,8x))
     355             : 8039 FORMAT (8x,60f9.5)
     356             : 8040 FORMAT (8x,4f9.5,9x,4l1,9x,i3)
     357             : 8050 FORMAT (19x,f9.5,9x,l1,15x,f9.5)
     358             : 
     359             :     ! output formats
     360             : 
     361             : 8138 FORMAT (' --> change   ',60(l1,8x))
     362             : 8139 FORMAT (' --> lo ',60f9.5)
     363             : 8140 FORMAT (' -->',i3,1x,4f9.5,' change: ',4l1,' skiplo: ',i3)
     364             : 8150 FORMAT ('  vacuum parameter=',f9.5,' change: ',l1, ' second vacuum=',f9.5)
     365             : 8001 FORMAT ('READING enpara for spin: ',i1)
     366             : 8000 FORMAT (/,' energy parameters:',/,t10,'s',t20, 'p',t30,'d',t37,'higher l - - -')
     367             : 8090 FORMAT ('Spin: ',i1,' -- ',i3,'eigenvalues')
     368             : 8091 FORMAT ('will be skipped for energyparameter computation')
     369             : 
     370          14 :     RETURN
     371             : 
     372           0 : 200 WRITE (6,*) 'the end of the file enpara has been reached while'
     373           0 :     WRITE (6,*) 'reading the energy-parameters.'
     374           0 :     WRITE (6,*) 'possible reason: energy parameters have not been'
     375           0 :     WRITE (6,*) 'specified for all atom types.'
     376           0 :     WRITE (6,FMT='(a,i4)') 'the actual number of atom-types is: ntype=',atoms%ntype
     377           0 :     CALL juDFT_error ("unexpected end of file enpara reached while reading")
     378             :   END SUBROUTINE read
     379             : 
     380             : 
     381          29 :   SUBROUTINE WRITE(enpara, atoms,jspins,film)
     382             : 
     383             :     ! write enpara-file
     384             :     !
     385             :     USE m_types_setup
     386             :     IMPLICIT NONE
     387             :     CLASS(t_enpara),INTENT(IN) :: enpara
     388             :     INTEGER, INTENT (IN) :: jspins
     389             :     LOGICAL,INTENT(IN)   :: film
     390             :     TYPE(t_atoms),INTENT(IN) :: atoms
     391             : 
     392             :     INTEGER n,l,lo,jspin
     393          58 :     REAL el0Temp(0:3), ello0Temp(1:atoms%nlod)
     394             : 
     395          29 :     OPEN(unit=40,file="enpara",form="formatted",status="replace")
     396             : 
     397          81 :     DO jspin=1,jspins
     398          52 :        WRITE (40,FMT=8035) jspin,enpara%enmix(jspin)
     399          52 :        WRITE (40,FMT=8036)
     400             : 8035   FORMAT (5x,'energy parameters          for spin ',i1,' mix=',f10.6)
     401             : 8036   FORMAT (t6,'atom',t15,'s',t24,'p',t33,'d',t42,'f')
     402         156 :        DO n = 1,atoms%ntype
     403         104 :           el0Temp(0:3) = enpara%el0(0:3,n,jspin)
     404         936 :           DO l = 0, 3
     405         520 :              IF (enpara%qn_el(l,n,jspin).NE.0) el0Temp(l) =  REAL(enpara%qn_el(l,n,jspin))
     406             :           END DO
     407         624 :           WRITE (6,FMT=8040)  n, (el0Temp(l),l=0,3),&
     408        1144 :                &                          (enpara%lchange(l,n,jspin),l=0,3),enpara%skiplo(n,jspin)
     409         624 :           WRITE (40,FMT=8040) n, (el0Temp(l),l=0,3),&
     410        1144 :                &                          (enpara%lchange(l,n,jspin),l=0,3),enpara%skiplo(n,jspin)
     411             :           !--->    energy parameters for the local orbitals
     412         156 :           IF (atoms%nlo(n).GE.1) THEN
     413           8 :              ello0Temp(1:atoms%nlo(n)) = enpara%ello0(1:atoms%nlo(n),n,jspin)
     414          24 :              DO lo = 1, atoms%nlo(n)
     415          16 :                 IF (enpara%qn_ello(lo,n,jspin).NE.0) ello0Temp(lo) = enpara%qn_ello(lo,n,jspin)
     416             :              END DO
     417           8 :              WRITE (6,FMT=8039) (ello0Temp(lo),lo=1,atoms%nlo(n))
     418           8 :              WRITE (6,FMT=8038) (enpara%llochg(lo,n,jspin),lo=1,atoms%nlo(n))
     419           8 :              WRITE (40,FMT=8039) (ello0Temp(lo),lo=1,atoms%nlo(n))
     420           8 :              WRITE (40,FMT=8038) (enpara%llochg(lo,n,jspin),lo=1,atoms%nlo(n))
     421             :           END IF
     422             : 
     423             :        ENDDO
     424             : 8038   FORMAT (' --> change   ',60(l1,8x))
     425             : 8039   FORMAT (' --> lo ',60f9.5)
     426             : 8040   FORMAT (' -->',i3,1x,4f9.5,' change: ',4l1,' skiplo: ',i3)
     427             : 
     428          81 :        IF (film) THEN
     429           4 :           WRITE (40,FMT=8050) enpara%evac(1,jspin),enpara%lchg_v(1,jspin),enpara%evac(2,jspin)
     430           4 :           WRITE (6,FMT=8050)  enpara%evac(1,jspin),enpara%lchg_v(1,jspin),enpara%evac(2,jspin)
     431             : 8050      FORMAT ('  vacuum parameter=',f9.5,' change: ',l1,&
     432             :                &           ' second vacuum=',f9.5)
     433             :        ENDIF
     434             :     ENDDO
     435          29 :     CLOSE(40)
     436          29 :     RETURN
     437             :   END SUBROUTINE WRITE
     438             : 
     439             : 
     440             : 
     441             : 
     442             : 
     443         322 :   SUBROUTINE mix(enpara,mpi,atoms,vacuum,input,vr,vz)
     444             :     !------------------------------------------------------------------
     445             :     USE m_types_setup
     446             :     USE m_types_mpi
     447             :     IMPLICIT NONE
     448             :     CLASS(t_enpara),INTENT(INOUT)  :: enpara
     449             :     TYPE(t_mpi),INTENT(IN)         :: mpi
     450             :     TYPE(t_atoms),INTENT(IN)       :: atoms
     451             :     TYPE(t_vacuum),INTENT(IN)      :: vacuum
     452             :     TYPE(t_input),INTENT(IN)       :: input
     453             : 
     454             :     REAL,    INTENT(IN) :: vr(:,:,:)
     455             :     REAL,    INTENT(IN) :: vz(vacuum%nmzd,2)
     456             : 
     457             :     INTEGER ityp,j,l,lo,jsp,n
     458             :     REAL    vbar,maxdist,maxdist2
     459         644 :     INTEGER same(atoms%nlod)
     460             :     LOGICAL l_enpara
     461             : #ifdef CPP_MPI
     462             :     INCLUDE 'mpif.h'
     463             :     INTEGER :: ierr
     464             : #endif    
     465             : 
     466         322 :     IF (mpi%irank==0) THEN
     467         161 :        maxdist2=0.0
     468         455 :        DO jsp=1,SIZE(enpara%el0,3)
     469         294 :           maxdist=0.0
     470         855 :           DO ityp = 1,atoms%ntype
     471             :              !        look for LO's energy parameters equal to the LAPW (and previous LO) ones
     472        1683 :              same = 0
     473         659 :              DO lo = 1,atoms%nlo(ityp)
     474          98 :                 IF(enpara%el0(atoms%llo(lo,ityp),ityp,jsp).EQ.enpara%ello0(lo,ityp,jsp)) same(lo)=-1
     475         695 :                 DO l = 1,lo-1
     476          36 :                    IF(atoms%llo(l,ityp).NE.atoms%llo(lo,ityp)) CYCLE
     477          98 :                    IF(enpara%ello0(l,ityp,jsp).EQ.enpara%ello0(lo,ityp,jsp).AND.same(lo).EQ.0) same(lo)=l
     478             :                 ENDDO
     479             :              ENDDO
     480             :              !
     481             :              !--->   change energy parameters
     482             :              !
     483        2805 :              DO l = 0,3
     484        2244 :                 WRITE(6,*) 'Type:',ityp,' l:',l
     485        2244 :                 WRITE(6,FMT=777) enpara%el0(l,ityp,jsp),enpara%el1(l,ityp,jsp),&
     486        4488 :                      ABS(enpara%el0(l,ityp,jsp)-enpara%el1(l,ityp,jsp))
     487        2244 :                 maxdist=MAX(maxdist,ABS(enpara%el0(l,ityp,jsp)-enpara%el1(l,ityp,jsp)))
     488        2805 :                 IF ( enpara%lchange(l,ityp,jsp) ) THEN
     489         360 :                    maxdist2=MAX(maxdist2,ABS(enpara%el0(l,ityp,jsp)-enpara%el1(l,ityp,jsp)))
     490             :                    enpara%el0(l,ityp,jsp) =(1.0-enpara%enmix(jsp))*enpara%el0(l,ityp,jsp) + &
     491         360 :                         enpara%enmix(jsp)*enpara%el1(l,ityp,jsp)
     492             :                 ENDIF
     493             :              ENDDO
     494         561 :              IF ( enpara%lchange(3,ityp,jsp) ) enpara%el0(4:,ityp,jsp) = enpara%el0(3,ityp,jsp)
     495             :              !
     496             :              !--->    determine and change local orbital energy parameters
     497             :              !
     498         659 :              DO lo = 1,atoms%nlo(ityp)
     499         659 :                 IF (atoms%l_dulo(lo,ityp)) THEN
     500           0 :                    enpara%ello0(lo,ityp,jsp) =enpara%el0(atoms%llo(lo,ityp),ityp,jsp)
     501             :                 ELSE
     502          98 :                    IF (enpara%llochg(lo,ityp,jsp) ) THEN
     503           0 :                       IF(same(lo).EQ.-1) THEN
     504           0 :                          enpara%ello0(lo,ityp,jsp) = enpara%el0(atoms%llo(lo,ityp),ityp,jsp)
     505           0 :                          CYCLE
     506           0 :                       ELSE IF(same(lo).GT.0) THEN
     507           0 :                          enpara%ello0(lo,ityp,jsp) = enpara%ello0(same(lo),ityp,jsp)
     508           0 :                          CYCLE
     509             :                       ENDIF
     510             :                    ENDIF
     511          98 :                    WRITE(6,*) 'Type:',ityp,' lo:',lo
     512          98 :                    WRITE(6,FMT=777) enpara%ello0(lo,ityp,jsp),enpara%ello1(lo,ityp,jsp),&
     513         196 :                         ABS(enpara%ello0(lo,ityp,jsp)-enpara%ello1(lo,ityp,jsp))
     514          98 :                    maxdist=MAX(maxdist,ABS(enpara%ello0(lo,ityp,jsp)-enpara%ello1(lo,ityp,jsp)))
     515          98 :                    IF (enpara%llochg(lo,ityp,jsp) ) THEN
     516           0 :                       maxdist2=MAX(maxdist2,ABS(enpara%ello0(lo,ityp,jsp)-enpara%ello1(lo,ityp,jsp)))
     517             :                       enpara%ello0(lo,ityp,jsp) =(1.0-enpara%enmix(jsp))*enpara%ello0(lo,ityp,jsp)+&
     518           0 :                            enpara%enmix(jsp)*enpara%ello1(lo,ityp,jsp)
     519             :                    ENDIF
     520             :                 END IF
     521             :              END DO
     522             :              !Shift if floating energy parameters are used
     523         855 :              IF (enpara%floating) THEN
     524           0 :                 j = atoms%jri(ityp) - (LOG(4.0)/atoms%dx(ityp)+1.51)
     525           0 :                 vbar = vr(j,ityp,jsp)/( atoms%rmt(ityp)*EXP(atoms%dx(ityp)*(j-atoms%jri(ityp))) )
     526           0 :                 enpara%el0(:,n,jsp)=enpara%el0(:,n,jsp)-vbar
     527             :              ENDIF
     528             :           END DO
     529             :           
     530             :           
     531         294 :           IF (input%film) THEN
     532           7 :              WRITE(6,*) 'Vacuum:'
     533          14 :              DO n=1,vacuum%nvac
     534           7 :                 WRITE(6,FMT=777) enpara%evac(n,jsp),enpara%evac1(n,jsp),ABS(enpara%evac(n,jsp)-enpara%evac1(n,jsp))
     535           7 :                 maxdist=MAX(maxdist,ABS(enpara%evac(n,jsp)-enpara%evac1(n,jsp)))
     536          14 :                 IF (enpara%lchg_v(n,jsp) ) THEN
     537           4 :                    maxdist2=MAX(maxdist2,ABS(enpara%evac(n,jsp)-enpara%evac1(n,jsp)))
     538             :                    enpara%evac(n,jsp) =(1.0-enpara%enmix(jsp))*enpara%evac(n,jsp)+&
     539           4 :                         enpara%enmix(jsp)*enpara%evac1(n,jsp)
     540             :                 END IF
     541             :              END DO
     542           7 :              IF (vacuum%nvac==1) enpara%evac(2,jsp) = enpara%evac(1,jsp)
     543           7 :              IF (enpara%floating) enpara%evac(:,jsp)=enpara%evac(:,jsp)-vz(:,jsp)
     544             :           ENDIF
     545         455 :           WRITE(6,'(a36,f12.6)') 'Max. mismatch of energy parameters:', maxdist
     546             :        END DO
     547         161 :        IF (maxdist2>1.0) CALL juDFT_warn&
     548             :             ("Energy parameter mismatch too large",hint&
     549             :             ="If any energy parameters calculated from the output "//&
     550             :             "differ from the input by more than 1Htr, chances are "//&
     551           0 :             "high that your initial setup was broken.")
     552             :     ENDIF
     553         322 :     INQUIRE(file='enpara',exist=l_enpara)
     554         322 :     IF (mpi%irank==0.AND.l_enpara) CALL enpara%WRITE(atoms,input%jspins,input%film)
     555             : #ifdef CPP_MPI
     556         322 :     CALL MPI_BCAST(enpara%el0,SIZE(enpara%el0),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
     557         322 :     CALL MPI_BCAST(enpara%ello0,SIZE(enpara%ello0),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
     558         322 :     CALL MPI_BCAST(enpara%evac,SIZE(enpara%evac),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
     559             : #endif    
     560         322 :     RETURN
     561             : 777 FORMAT('Old:',f8.5,' new:',f8.5,' diff:',f8.5)
     562             : 
     563             :   END SUBROUTINE mix
     564             : 
     565         324 :   SUBROUTINE calcOutParams(enpara,input,atoms,vacuum,regCharges)
     566             :     USE m_types_setup
     567             :     USE m_types_regionCharges
     568             :     IMPLICIT NONE
     569             :     CLASS(t_enpara),INTENT(INOUT)    :: enpara
     570             :     TYPE(t_input),INTENT(IN)         :: input
     571             :     TYPE(t_atoms),INTENT(IN)         :: atoms
     572             :     TYPE(t_vacuum),INTENT(IN)        :: vacuum
     573             :     TYPE(t_regionCharges),INTENT(IN) :: regCharges
     574             : 
     575             :     INTEGER :: ispin, n
     576             : 
     577         914 :     DO ispin = 1,input%jspins
     578        1716 :        DO n=1,atoms%ntype
     579        1126 :           enpara%el1(0:3,n,ispin)=regCharges%ener(0:3,n,ispin)/regCharges%sqal(0:3,n,ispin)
     580        1716 :           IF (atoms%nlo(n)>0) enpara%ello1(:atoms%nlo(n),n,ispin)=regCharges%enerlo(:atoms%nlo(n),n,ispin)/regCharges%sqlo(:atoms%nlo(n),n,ispin)
     581             :        END DO
     582         914 :        IF (input%film) enpara%evac1(:vacuum%nvac,ispin)=regCharges%pvac(:vacuum%nvac,ispin)/regCharges%svac(:vacuum%nvac,ispin)
     583             :     END DO
     584         324 :   END SUBROUTINE calcOutParams
     585             : 
     586        1998 :   SUBROUTINE priv_write(lo,l,n,jsp,nqn,e_lo,e_up,e)
     587             :     !subroutine to write energy parameters to output
     588             :     USE m_xmlOutput
     589             :     IMPLICIT NONE
     590             :     LOGICAL,INTENT(IN):: lo
     591             :     INTEGER,INTENT(IN):: l,n,jsp,nqn
     592             :     REAL,INTENT(IN)   :: e_lo,e_up,e
     593             : 
     594             :     CHARACTER(LEN=20)    :: attributes(6)
     595        1998 :     CHARACTER(len=:),ALLOCATABLE:: label
     596             :     CHARACTER(len=1),PARAMETER,DIMENSION(0:9):: ch=(/'s','p','d','f','g','h','i','j','k','l'/)
     597             : 
     598       13986 :     attributes = ''
     599        1998 :     WRITE(attributes(1),'(i0)') n
     600        1998 :     WRITE(attributes(2),'(i0)') jsp
     601        1998 :     WRITE(attributes(3),'(i0,a1)') abs(nqn), ch(l)
     602        1998 :     WRITE(attributes(4),'(f8.2)') e_lo
     603        1998 :     WRITE(attributes(5),'(f8.2)') e_up
     604        1998 :     WRITE(attributes(6),'(f16.10)') e
     605        1998 :     IF (nqn>0) THEN
     606        1998 :        IF (lo) THEN
     607          94 :           label='loAtomicEP'
     608             :        ELSE
     609        1904 :           label='atomicEP'
     610             :        ENDIF
     611             :     ELSE
     612           0 :        IF (lo) THEN
     613           0 :           label='heloAtomicEP'
     614             :        ELSE
     615           0 :           label='heAtomicEP'
     616             :        ENDIF
     617             :     END IF
     618             : 
     619             :     CALL writeXMLElementForm(label,(/'atomType     ','spin         ','branch       ',&
     620             :          'branchLowest ','branchHighest','value        '/),&
     621        1998 :          attributes,RESHAPE((/10,4,6,12,13,5,6,1,3,8,8,16/),(/6,2/)))
     622        1998 :     IF (lo) THEN
     623          94 :        WRITE(6,'(a6,i5,i2,a1,a12,f6.2,a3,f6.2,a13,f8.4)') '  Atom',n,nqn,ch(l),' branch from',&
     624         188 :          e_lo, ' to',e_up,' htr. ; e_l(lo) =',e
     625             :     ELSE
     626        1904 :        WRITE(6,'(a6,i5,i2,a1,a12,f6.2,a3,f6.2,a13,f8.4)') '  Atom',n,nqn,ch(l),' branch from',&
     627        3808 :          e_lo, ' to',e_up,' htr. ; e_l =',e
     628             :     END IF
     629        1998 :   END SUBROUTINE priv_write
     630             : 
     631             : 
     632           0 : END MODULE m_types_enpara

Generated by: LCOV version 1.13