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_lapwdim
       8             : CONTAINS
       9         160 :   SUBROUTINE lapw_dim(kpts,cell,input,noco,nococonv,forcetheo,atoms,nbasfcn,juPhon)
      10             :     !
      11             :     !*********************************************************************
      12             :     !     determines dimensions of the lapw basis set with |k+G|<rkmax.
      13             :     !  Generalization of the old apws_dim routine
      14             :     !*********************************************************************
      15             :     USE m_boxdim
      16             :     USE m_types_fleurinput
      17             :     USE m_types_forcetheo_extended
      18             :     IMPLICIT NONE
      19             :     TYPE(t_kpts),INTENT(IN)      :: kpts
      20             :     TYPE(t_cell),INTENT(IN)      :: cell
      21             :     TYPE(t_input),INTENT(IN)     :: input
      22             :     TYPE(t_nococonv),INTENT(IN)  :: nococonv
      23             :     TYPE(t_noco),INTENT(IN)      :: noco
      24             :      
      25             :     CLASS(t_forcetheo),INTENT(IN):: forcetheo
      26             :     TYPE(t_atoms),INTENT(IN)     :: atoms
      27             :     INTEGER, INTENT(OUT)         :: nbasfcn
      28             : 
      29             :     TYPE(t_juPhon), INTENT(IN)   :: juPhon
      30             : 
      31             :     !local variable for init
      32             :     INTEGER               :: nvd,nv2d,addx,addy,addz
      33         160 :     TYPE(t_lapw) :: lapw
      34             : 
      35             :     INTEGER j1,j2,j3,mk1,mk2,mk3,iofile,ksfft,q,nk,nv,nv2
      36             :     INTEGER ispin,nvh(2),nv2h(2)
      37             : 
      38             :     REAL arltv1,arltv2,arltv3,rkm,rk2,r2,s(3),gmaxp,qss(3)
      39         160 :     REAL,ALLOCATABLE:: q_vectors(:,:)
      40             :     REAL            :: bkpt(3)
      41             :     ! ..
      42             :     !
      43             :     !------->          ABBREVIATIONS
      44             :     !
      45             : 
      46             :     !   iofile      : device number for in and output
      47             :     !   gmax        : cut-off wavevector for charge density
      48             :     !   rkmax       : cut-off for |g+k|
      49             :     !   gmaxp       : gmaxp = gmax/rkmax, ideal: gmaxp=2
      50             :     !   arltv(i)    : length of reciprical lattice vector along
      51             :     !                 direction (i)
      52             :     !
      53             :     !---> Determine rkmax box of size mk1, mk2, mk3,
      54             :     !     for which |G(mk1,mk2,mk3) + (k1,k2,k3)| < rkmax
      55             :     !
      56         160 :     CALL boxdim(cell%bmat,arltv1,arltv2,arltv3)
      57             : 
      58             :     !     (add 1+1 due to integer rounding, strange k_vector in BZ)
      59         160 :     mk1 = int(input%rkmax/arltv1) + 2
      60         160 :     mk2 = int(input%rkmax/arltv2) + 2
      61         160 :     mk3 = int(input%rkmax/arltv3) + 2
      62             : 
      63         160 :     rkm = input%rkmax
      64         160 :     rk2 = rkm*rkm
      65             : 
      66             :     !Determine the q-vector(s) to use
      67         160 :     IF (juPhon%l_dfpt) THEN
      68           0 :        ALLOCATE(q_vectors(3,SIZE(juPhon%qvec,2)))
      69           0 :        q_vectors=juPhon%qvec
      70             :     ELSE
      71             :        SELECT TYPE(forcetheo)
      72             :        TYPE IS (t_forcetheo_ssdisp)
      73           0 :           ALLOCATE(q_vectors(3,SIZE(forcetheo%qvec,2)))
      74           0 :           q_vectors=forcetheo%qvec
      75             :        TYPE IS (t_forcetheo_dmi)
      76           0 :           ALLOCATE(q_vectors(3,SIZE(forcetheo%qvec,2)))
      77           0 :           q_vectors=forcetheo%qvec
      78             :        TYPE IS (t_forcetheo_jij)
      79           0 :           ALLOCATE(q_vectors(3,SIZE(forcetheo%qvec,2)))
      80           0 :           q_vectors=forcetheo%qvec
      81             :        CLASS IS (t_forcetheo) ! DEFAULT
      82         160 :           ALLOCATE(q_vectors(3,1))
      83         640 :           q_vectors(:,1)=nococonv%qss
      84             :        END SELECT
      85             :     END IF
      86             : 
      87         640 :     if (any(abs(nococonv%qss-q_vectors(:,1))>1E-4)) CALL judft_warn("q-vector for self-consistency should be first in list for force-theorem")
      88             : 
      89             : 
      90         160 :     nvd = 0 ; nv2d = 0
      91         320 :     DO q=1,SIZE(q_vectors,2)
      92         640 :        qss=q_vectors(:,q)
      93        1596 :        DO nk=1,kpts%nkpt
      94        5104 :           bkpt=kpts%bk(:,nk)
      95             :           !---> obtain vectors
      96             :           !---> in a spin-spiral calculation different basis sets are used for
      97             :           !---> the two spin directions, because the cutoff radius is defined
      98             :           !---> by |G + k +/- qss/2| < rkmax.
      99        3828 :           DO ispin = 1,2
     100        2552 :              addX = abs(NINT((bkpt(1) + (2*ispin - 3)/2.0*qss(1))/arltv1))
     101        2552 :              addY = abs(NINT((bkpt(2) + (2*ispin - 3)/2.0*qss(2))/arltv2))
     102        2552 :              addZ = abs(NINT((bkpt(3) + (2*ispin - 3)/2.0*qss(3))/arltv2))
     103        2552 :              nv = 0
     104        2552 :              nv2 = 0
     105       28656 :              DO j1 = -mk1-addX,mk1+addX
     106       26104 :                 s(1) = bkpt(1) + j1 + (2*ispin - 3)/2.0*qss(1)
     107      306072 :                 DO j2 = -mk2-addY,mk2+addY
     108      277416 :                    s(2) = bkpt(2) + j2 + (2*ispin - 3)/2.0*qss(2)
     109             :                    !--->          nv2 for films
     110      277416 :                    s(3) = 0.0
     111             :                    !r2 = dotirp(s,s,cell%bbmat)
     112     4438656 :                    r2 = dot_product(matmul(s,cell%bbmat),s)
     113      277416 :                    IF (r2.LE.rk2) nv2 = nv2 + 1
     114     3530848 :                    DO j3 = -mk3-addz,mk3+addz
     115     3227328 :                       s(3) = bkpt(3) + j3 + (2*ispin - 3)/2.0*qss(3)
     116             :                       !r2 = dotirp(s,s,cell%bbmat)
     117    51637248 :                       r2 = dot_product(matmul(s,cell%bbmat),s)
     118     3504744 :                       IF (r2.LE.rk2) THEN
     119      412068 :                          nv = nv + 1
     120             :                       END IF
     121             :                    END DO
     122             :                 END DO
     123             :              END DO
     124             :         
     125        2552 :              nvh(ispin)  = nv
     126        3828 :              nv2h(ispin) = nv2
     127             :           END DO
     128        1276 :           nvd=MAX(nvd,MAX(nvh(1),nvh(2)))
     129        1436 :           nv2d=MAX(nv2d,MAX(nv2h(1),nv2h(2)))
     130             : 
     131             :        ENDDO !k-loop
     132             :     ENDDO !q-loop
     133             : 
     134         160 :     nbasfcn = nvd + atoms%nat*atoms%nlod*(2*atoms%llod+1)
     135         160 :     IF (noco%l_noco) nbasfcn = 2*nbasfcn
     136         160 :     call lapw%init_dim(nvd,nv2d,nbasfcn)
     137             : 
     138         160 :   END SUBROUTINE lapw_dim
     139             : 
     140             :   
     141         160 : END MODULE m_lapwdim