dimen7.F90 13.6 KB
Newer Older
1 2 3 4 5 6
      MODULE m_dimen7
      use m_juDFT
      CONTAINS
      SUBROUTINE dimen7(&
     &                  input,sym,stars,&
     &                  atoms,sphhar,dimension,vacuum,&
7
     &                  kpts,oneD,hybrid,cell)
8 9 10 11 12 13 14 15 16 17 18 19 20 21

!
! This program reads the input files of the flapw-programm (inp & kpts)
! and creates a file 'fl7para' that contains dimensions 
! for the main flapw-programm.
!

      USE m_localsym
      USE m_socsym
      USE m_sssym
      USE m_spg2set
      USE m_constants
      USE m_rwinp
      USE m_inpnoco
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
!      USE m_julia
!      USE m_od_kptsgen
      USE m_types_input
      USE m_types_sym
      USE m_types_stars
      USE m_types_atoms
      USE m_types_sphhar
      USE m_types_dimension
      USE m_types_vacuum
      USE m_types_kpts
      USE m_types_oneD
      USE m_types_hybrid
      USE m_types_cell
      USE m_types_noco
      USE m_types_banddos
      USE m_types_sliceplot
      USE m_types_xcpot_inbuild_nofunction


41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
      USE m_firstglance
      USE m_inv3
      USE m_rwsymfile
      USE m_strgndim
      USE m_convndim
      USE m_inpeigdim
      USE m_ylm
      IMPLICIT NONE
!
! dimension-parameters for flapw:
!
      TYPE(t_input),INTENT(INOUT)   :: input
      TYPE(t_sym),INTENT(INOUT)     :: sym
      TYPE(t_stars),INTENT(INOUT)   :: stars 
      TYPE(t_atoms),INTENT(INOUT)   :: atoms
      TYPE(t_sphhar),INTENT(INOUT)  :: sphhar
      TYPE(t_dimension),INTENT(INOUT) :: dimension
      TYPE(t_vacuum),INTENT(INOUT)   :: vacuum
      TYPE(t_kpts),INTENT(INOUT)     :: kpts
      TYPE(t_oneD),INTENT(INOUT)     :: oneD
      TYPE(t_hybrid),INTENT(INOUT)   :: hybrid
      TYPE(t_cell),INTENT(INOUT)     :: cell
 
      TYPE(t_noco)      :: noco
      TYPE(t_sliceplot) :: sliceplot
      TYPE(t_banddos)   :: banddos
67
      TYPE(t_xcpot_inbuild_nf)     :: xcpot
68 69 70 71 72 73

!
!
!-------------------------------------------------------------------
! ..  Local Scalars ..
      REAL   :: thetad,xa,epsdisp,epsforce ,rmtmax,arltv1,arltv2,arltv3   
74
      REAL   :: s,r,d ,idsprs
75
      INTEGER :: ok,ilo,n,nstate,i,j,na,n1,n2,jrc,nopd,symfh
76
      INTEGER :: nmopq(3)
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
      
      CHARACTER(len=1) :: rw
      CHARACTER(len=4) :: namex 
      CHARACTER(len=7) :: symfn
      CHARACTER(len=12):: relcor
      LOGICAL  ::l_kpts,l_qpts,l_inpexist,l_tmp(2)
! ..
      REAL    :: a1(3),a2(3),a3(3)  
      REAL    :: q(3)

      CHARACTER(len=3), ALLOCATABLE :: noel(:)
      LOGICAL, ALLOCATABLE :: error(:) 
     
      INTEGER ntp1,ii
      INTEGER, ALLOCATABLE :: lmx1(:), nq1(:), nlhtp1(:)

!     added for HF and hybrid functionals
      LOGICAL               ::  l_gamma=.false.
95
      character(len=4) :: latnam
96

97
      EXTERNAL prp_xcfft_box!,parawrite
98 99 100 101 102 103 104
!     ..
      
    
!---> First, check whether an inp-file exists
!
      INQUIRE (file='inp',exist=l_inpexist)
      IF (.not.l_inpexist) THEN
105
         CALL juDFT_error("no inp- or input-file found!",calledby ="dimen7")
106 107 108 109
      ENDIF
!
!---> determine ntype,nop,natd,nwdd,nlod and layerd
!
Daniel Wortmann's avatar
Daniel Wortmann committed
110
      CALL first_glance(atoms%ntype,sym%nop,atoms%nat,atoms%nlod,vacuum%layerd,&
111
                        input%itmax,l_kpts,l_qpts,l_gamma,kpts%nkpt,kpts%nkpt3,nmopq)
Daniel Wortmann's avatar
Daniel Wortmann committed
112
      atoms%ntype=atoms%ntype
113 114 115
      atoms%nlod = max(atoms%nlod,1)

      ALLOCATE (&
Daniel Wortmann's avatar
Daniel Wortmann committed
116
     & atoms%lmax(atoms%ntype),atoms%ntypsy(atoms%nat),atoms%neq(atoms%ntype),atoms%nlhtyp(atoms%ntype),&
117 118
     & atoms%rmt(atoms%ntype),atoms%zatom(atoms%ntype),atoms%jri(atoms%ntype),atoms%dx(atoms%ntype), &
     & atoms%nlo(atoms%ntype),atoms%llo(atoms%nlod,atoms%ntype),atoms%nflip(atoms%ntype),atoms%bmu(atoms%ntype),&
119
     & noel(atoms%ntype),vacuum%izlay(vacuum%layerd,2),atoms%econf(atoms%ntype),atoms%lnonsph(atoms%ntype),&
Daniel Wortmann's avatar
Daniel Wortmann committed
120
     & atoms%taual(3,atoms%nat),atoms%pos(3,atoms%nat),&
121
     & atoms%nz(atoms%ntype),atoms%relax(3,atoms%ntype),&
122
     & atoms%l_geo(atoms%ntype),noco%alph(atoms%ntype),noco%beta(atoms%ntype),&
123 124
     & atoms%lda_u(atoms%ntype),noco%l_relax(atoms%ntype),&
     & noco%b_con(2,atoms%ntype),&
125
     & sphhar%clnu(1,1,1),sphhar%nlh(1),sphhar%llh(1,1),sphhar%nmem(1,1),sphhar%mlh(1,1,1),&
Daniel Wortmann's avatar
Daniel Wortmann committed
126
     & hybrid%select1(4,atoms%ntype),hybrid%lcutm1(atoms%ntype),&
127 128 129 130 131
     & hybrid%lcutwf(atoms%ntype), STAT=ok)
!
!---> read complete input and calculate nvacd,llod,lmaxd,jmtd,neigd and 
!
      CALL rw_inp('r',&
132
     &            atoms,vacuum,input,stars,sliceplot,banddos,&
133
     &                  cell,sym,xcpot,noco,oneD,hybrid,kpts,&
134
     &                  noel,namex,relcor,a1,a2,a3,latnam)
135 136 137 138 139 140

!---> pk non-collinear
!---> read the angle and spin-spiral information from nocoinp
      noco%qss = 0.0
      noco%l_ss = .false.
      IF (noco%l_noco) THEN 
141
         CALL inpnoco(atoms,input,vacuum,noco)
142 143 144 145 146 147 148 149 150
      ENDIF

      vacuum%nvacd = 2
      IF (sym%zrfs .OR. sym%invs .OR. oneD%odd%d1) vacuum%nvacd = 1
      atoms%llod  = 0
      atoms%lmaxd = 0
      atoms%jmtd  = 0
      rmtmax      = 0.0
      dimension%neigd = 0
151
      dimension%nstd  = maxval(atoms%econf%num_core_states)
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176
      atoms%lmaxd = maxval(atoms%lmax)
      atoms%jmtd  = maxval(atoms%jri)
      rmtmax      = maxval(atoms%rmt)
      DO n = 1,atoms%ntype
        DO ilo = 1,atoms%nlo(n)
!+apw
          IF (atoms%llo(ilo,n).LT.0) THEN
             atoms%llo(ilo,n) = -atoms%llo(ilo,n) - 1
          ELSE
             dimension%neigd = dimension%neigd + atoms%neq(n)*(2*abs(atoms%llo(ilo,n)) +1)
          ENDIF
!-apw
          atoms%llod = max(abs(atoms%llo(ilo,n)),atoms%llod)
        ENDDO
        nstate = 4
        IF ((atoms%nz(n).GE.21.AND.atoms%nz(n).LE.29) .OR. &
     &      (atoms%nz(n).GE.39.AND.atoms%nz(n).LE.47) .OR.&
     &      (atoms%nz(n).GE.57.AND.atoms%nz(n).LE.79)) nstate = 9
        IF ((atoms%nz(n).GE.58.AND.atoms%nz(n).LE.71) .OR.&
     &      (atoms%nz(n).GE.90.AND.atoms%nz(n).LE.103)) nstate = 16
        dimension%neigd = dimension%neigd + nstate*atoms%neq(n)
!
      ENDDO
      CALL ylmnorm_init(atoms%lmaxd)
!      IF (mod(lmaxd,2).NE.0) lmaxd = lmaxd + 1
177
      IF (2*DIMENSION%neigd.LT.MAX(5.0,input%zelec)) THEN
178
        WRITE(6,*) dimension%neigd,' states estimated in dimen7 ...'
179
        DIMENSION%neigd = MAX(5,NINT(0.75*input%zelec))
180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
        WRITE(6,*) 'changed dimension%neigd to ',dimension%neigd
      ENDIF
      IF (noco%l_soc .and. (.not. noco%l_noco)) dimension%neigd=2*dimension%neigd 
      IF (noco%l_soc .and. noco%l_ss) dimension%neigd=(3*dimension%neigd)/2  
       ! not as accurate, but saves much time

      rmtmax = rmtmax*stars%gmax
      CALL convn_dim(rmtmax,dimension%ncvd)
!
! determine core mesh
!
      dimension%msh = 0
      DO n = 1,atoms%ntype
         r = atoms%rmt(n)
         d = exp(atoms%dx(n))
         jrc = atoms%jri(n)
         DO WHILE (r < atoms%rmt(n) + 20.0)
            jrc = jrc + 1
            r = r*d
         ENDDO
         dimension%msh = max( dimension%msh, jrc ) 
      ENDDO
!
! ---> now, set the lattice harmonics, determine nlhd
!
205 206 207
      cell%amat(:,1) = a1(:)*input%scaleCell
      cell%amat(:,2) = a2(:)*input%scaleCell
      cell%amat(:,3) = a3(:)*input%scaleCell
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
      CALL inv3(cell%amat,cell%bmat,cell%omtil)
      IF (input%film) cell%omtil = cell%omtil/cell%amat(3,3)*vacuum%dvac
!-odim
      IF (oneD%odd%d1) cell%omtil = cell%amat(3,3)*pimach()*(vacuum%dvac**2)/4.
!+odim
      cell%bmat=tpi_const*cell%bmat
    
      na = 0
      DO n = 1,atoms%ntype
        DO n1 = 1,atoms%neq(n)
            na = na + 1
            IF (input%film) atoms%taual(3,na) = atoms%taual(3,na)/a3(3)
            atoms%pos(:,na) = matmul(cell%amat,atoms%taual(:,na))
        ENDDO
        atoms%zatom(n) = real( atoms%nz(n) )
      ENDDO
      ALLOCATE (sym%mrot(3,3,sym%nop),sym%tau(3,sym%nop))
      IF (sym%namgrp.EQ.'any ') THEN
         nopd = sym%nop ; rw = 'R'
         symfh = 94 ; symfn = 'sym.out'
         CALL rw_symfile(rw,symfh,symfn,nopd,cell%bmat,sym%mrot,sym%tau,sym%nop,sym%nop2,sym%symor)
      ELSE
230
         CALL spg2set(sym%nop,sym%zrfs,sym%invs,sym%namgrp,latnam,sym%mrot,sym%tau,sym%nop2,sym%symor)
231 232 233 234
      ENDIF
      sphhar%ntypsd = 0
      IF (.NOT.oneD%odd%d1) THEN
        CALL local_sym(atoms%lmaxd,atoms%lmax,sym%nop,sym%mrot,sym%tau,&
Daniel Wortmann's avatar
Daniel Wortmann committed
235
                       atoms%nat,atoms%ntype,atoms%neq,cell%amat,cell%bmat,&
236 237 238 239 240
                       atoms%taual,sphhar%nlhd,sphhar%memd,sphhar%ntypsd,.true.,&
                       atoms%nlhtyp,atoms%ntypsy,sphhar%nlh,sphhar%llh,&
                       sphhar%nmem,sphhar%mlh,sphhar%clnu)
!-odim
      ELSEIF (oneD%odd%d1) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
241
        ntp1 = atoms%nat
242 243 244 245 246 247 248 249 250 251
        ALLOCATE (nq1(ntp1),lmx1(ntp1),nlhtp1(ntp1))
        ii = 1
        nq1=1
        DO i = 1,atoms%ntype
          DO j = 1,atoms%neq(i)
            lmx1(ii) = atoms%lmax(i)
            ii = ii + 1
          END DO
        END DO
        CALL local_sym(atoms%lmaxd,lmx1,sym%nop,sym%mrot,sym%tau,&
Daniel Wortmann's avatar
Daniel Wortmann committed
252
              atoms%nat,ntp1,nq1,cell%amat,cell%bmat,atoms%taual,&
253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
              sphhar%nlhd,sphhar%memd,sphhar%ntypsd,.true.,nlhtp1,&
              atoms%ntypsy,sphhar%nlh,sphhar%llh,sphhar%nmem,&
              sphhar%mlh,sphhar%clnu)        
        ii = 1
        DO i = 1,atoms%ntype
          atoms%nlhtyp(i) = nlhtp1(ii)
          ii = ii + atoms%neq(i)
        END DO
        DEALLOCATE (nq1,lmx1,nlhtp1)
      END IF
!+odim
!
! Check if symmetry is compatible with SOC or SSDW
!
      IF (noco%l_soc .and. (.not.noco%l_noco)) THEN  
        ! test symmetry for spin-orbit coupling
        ALLOCATE ( error(sym%nop) )
270
        CALL soc_sym(sym%nop,sym%mrot,noco%theta,noco%phi,cell%amat,error)
271 272
        IF ( ANY(error(:)) ) THEN
          WRITE(*,fmt='(1x)')
273 274
          WRITE(*,fmt='(A)') 'Symmetry incompatible with SOC spin-quantization axis ,'  
          WRITE(*,fmt='(A)') 'do not perform self-consistent calculations !'    
275 276 277 278 279
          WRITE(*,fmt='(1x)')
          IF ( input%eonly .or. (noco%l_soc.and.noco%l_ss) .or. input%gw.ne.0 ) THEN  ! .or. .
            CONTINUE 
          ELSE 
            IF (input%itmax>1) THEN
280
               CALL juDFT_error("symmetry & SOC",calledby ="dimen7")
281 282 283 284 285 286 287 288
            ENDIF 
          ENDIF 
        ENDIF           
        DEALLOCATE ( error )
      ENDIF
      IF (noco%l_ss) THEN  ! test symmetry for spin-spiral
        ALLOCATE ( error(sym%nop) )
        CALL ss_sym(sym%nop,sym%mrot,noco%qss,error)
289
        IF ( ANY(error(:)) )  CALL juDFT_error("symmetry & SSDW", calledby="dimen7")
290 291 292 293 294 295 296
        DEALLOCATE ( error )
      ENDIF
!
! Dimensioning of the stars
!
      IF (input%film.OR.(sym%namgrp.ne.'any ')) THEN
         CALL strgn1_dim(stars%gmax,cell%bmat,sym%invs,sym%zrfs,sym%mrot,&
297 298
                    sym%tau,sym%nop,sym%nop2,stars%mx1,stars%mx2,stars%mx3,&
                    stars%ng3,stars%ng2,oneD%odd)
299 300 301

      ELSE
         CALL strgn2_dim(stars%gmax,cell%bmat,sym%invs,sym%zrfs,sym%mrot,&
302 303 304 305
                    sym%tau,sym%nop,stars%mx1,stars%mx2,stars%mx3,&
                    stars%ng3,stars%ng2)
         oneD%odd%n2d = stars%ng2
         oneD%odd%nq2 = stars%ng2
306 307 308 309
         oneD%odd%nop = sym%nop
      ENDIF

      IF ( xcpot%gmaxxc .le. 10.0**(-6) ) THEN
310 311
         WRITE (6,'(" xcpot%gmaxxc=0 : xcpot%gmaxxc=stars%gmax choosen as default value")')
         WRITE (6,'(" concerning memory, you may want to choose a smaller value for stars%gmax")')
312 313 314
         xcpot%gmaxxc=stars%gmax
      END IF

315
      CALL prp_xcfft_box(xcpot%gmaxxc,cell%bmat,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft)
316 317 318 319 320
!
! k-point generator provides kpts-file, if it's missing:
!
      IF (.not.l_kpts) THEN
       IF (.NOT.oneD%odd%d1) THEN
321
          IF(l_gamma .AND. banddos%ndir .EQ. 0) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
322
         call judft_error("gamma swtich not supported in old inp file anymore",calledby="dimen7")
323
         ELSE
324
!         CALL julia(sym,cell,input,noco,banddos,kpts,.false.,.FALSE.)
325 326
         ENDIF
       ELSE
327
!        CALL od_kptsgen (kpts%nkpt)
328 329 330 331 332
       ENDIF
      ELSE
        IF(input%gw.eq.2) THEN
          INQUIRE(file='QGpsi',exist=l_kpts) ! Use QGpsi if it exists ot
          IF(l_kpts) THEN
333 334
            WRITE(6,*) 'QGpsi exists and will be used to generate kpts-file'
            OPEN (15,file='QGpsi',form='unformatted',status='old',action='read')
335 336 337 338 339 340
            OPEN (41,file='kpts',form='formatted',status='unknown')
            REWIND(41)
            READ (15) kpts%nkpt
            WRITE (41,'(i5,f20.10)') kpts%nkpt,1.0
            DO n = 1, kpts%nkpt
              READ (15) q
341
              WRITE (41,'(4f10.5)') MATMUL(TRANSPOSE(cell%amat),q)/input%scaleCell,1.0
342 343 344 345 346 347 348 349 350 351 352 353 354
              READ (15)
            ENDDO
            CLOSE (15)
            CLOSE (41)
          ENDIF
        ENDIF
      ENDIF
      
      dimension%neigd = max(dimension%neigd,input%gw_neigd)

!
! Using the k-point generator also for creation of q-points for the
! J-constants calculation:
355 356 357 358 359 360 361 362
!      IF(.not.l_qpts)THEN
!        kpts%nkpt3=nmopq
!        l_tmp=(/noco%l_ss,noco%l_soc/)
!        noco%l_ss=.false.
!        noco%l_soc=.false.
!        CALL julia(sym,cell,input,noco,banddos,kpts,.true.,.FALSE.)
!        noco%l_ss=l_tmp(1); noco%l_soc=l_tmp(2)
!      ENDIF
363 364 365 366

!
! now proceed as usual
!
367
      CALL inpeig_dim(input,cell,noco,oneD,kpts,dimension,stars,latnam)
368 369
      vacuum%layerd = max(vacuum%layerd,1)
      dimension%nstd = max(dimension%nstd,30)
Daniel Wortmann's avatar
Daniel Wortmann committed
370
      atoms%ntype = atoms%ntype
371 372 373
      IF (noco%l_noco) dimension%neigd = 2*dimension%neigd

      atoms%nlod = max(atoms%nlod,2) ! for chkmt
374
      input%jspins=input%jspins
375
      !CALL parawrite(sym,stars,atoms,sphhar,DIMENSION,vacuum,kpts,oneD,input)
376

377
      DEALLOCATE( sym%mrot,sym%tau,&
378
     & atoms%lmax,atoms%ntypsy,atoms%neq,atoms%nlhtyp,atoms%rmt,atoms%zatom,atoms%jri,atoms%dx,atoms%nlo,atoms%llo,atoms%nflip,atoms%bmu,noel,&
379
     & vacuum%izlay,atoms%econf,atoms%lnonsph,atoms%taual,atoms%pos,atoms%nz,atoms%relax,&
380 381
     & atoms%l_geo,noco%alph,noco%beta,atoms%lda_u,noco%l_relax,noco%b_con,sphhar%clnu,sphhar%nlh,&
     & sphhar%llh,sphhar%nmem,sphhar%mlh,hybrid%select1,hybrid%lcutm1,&
382
     & hybrid%lcutwf)
383

384 385 386
      RETURN
      END SUBROUTINE dimen7
      END MODULE m_dimen7