LCOV - code coverage report
Current view: top level - kpoints - kptmop.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 108 279 38.7 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_kptmop
       2             :       use m_juDFT
       3             : !-----------------------------------------------------------------------
       4             : ! ---> This program generates k-points
       5             : !           in irreducible wedge of BZ  (for nreg=0)
       6             : !           in total BZ                 (for nreg=1)
       7             : !      (BZ = 1. Brillouin-zone) for all canonical Bravais lattices
       8             : !      in 2 and 3 dimensions,
       9             : !      using the basis vectors of the reciprocal lattice
      10             : !      and the bordering planes of the irreducible wedge.
      11             : !
      12             : !      The k-points are generated by the Monkhorst--Pack method.
      13             : !      The information on the bordering planes of the irr wedge
      14             : !      is taken from BRZONE.
      15             : !
      16             : !      The program checks the compatibility of the dimension and
      17             : !      symmetry and of the provided Monkhorst-Pack-parameters.
      18             : !-----------------------------------------------------------------------
      19             :       CONTAINS
      20          25 :       SUBROUTINE kptmop(
      21             :      >                  iofile,iokpt,kpri,ktest,
      22             :      >                  idsyst,idtype,nmop,ikzero,kzero,
      23             :      >                  rltv,bltv,nreg,nfulst,nbound,idimens,
      24          50 :      >                  xvec,fnorm,fdist,ncorn,nface,nedge,cpoint,
      25             :      >                  nsym,ccr,rlsymr,talfa,mkpt,mface,mdir,
      26          25 :      <                  nkpt,divis,vkxyz,nkstar,wghtkp)
      27             : c
      28             : c    Meaning of variables:
      29             : c    INPUT:
      30             : c
      31             : c    Symmetry of lattice:
      32             : c    idsyst   : crystal system identification in MDDFT programs
      33             : c    idtype   : lattice type identification in MDDFT programs
      34             : c    bltv     : cartesian coordinates of basis vectors for
      35             : c               Bravais lattice: bltv(ix,jn), ix=1,3; jn=1,3
      36             : c    rltv     : cartesian coordinates of basis vectors for
      37             : c               reciprocal lattice: rltv(ix,jn), ix=1,3; jn=1,3
      38             : c    nsym     : number of symmetry elements of points group
      39             : c    ccr     : rotation matrix for symmetry element
      40             : c                   in cartesian representation
      41             : c    rlsymr   : rotation matrix for symmetry element
      42             : c                   in reciprocal lattice basis representation
      43             : c    talfa    : translation vector associated with (non-symmorphic)
      44             : c               symmetry elements in Bravais lattice representation
      45             : c
      46             : c    representation of the irreducible part of the BZ:
      47             : c    fnorm    : normal vector of the planes bordering the irrBZ
      48             : c    fdist    : distance vector of the planes bordering the irrBZ
      49             : c    ncorn    : number of corners of the irrBZ
      50             : c    nface    : number of faces of the irrBZ
      51             : c    nedge    : number of edges of the irrBZ
      52             : c    xvec     : arbitrary vector lying in the irrBZ (FOR SURE!!)
      53             : c               components are:
      54             : c
      55             : c    characterization of the Monkhorst-Pack k-point set:
      56             : c    idimens  : number of dimensions for k-point set (2 or 3)
      57             : c    nreg     : 1 kpoints in full BZ; 0 kpoints in irrBZ
      58             : c    nfulst   : 1 kpoints ordered in full stars 
      59             : c                  (meaningful only for nreg =1; full BZ)
      60             : c    nbound   : 0 no primary points on BZ boundary; 
      61             : c               1 with boundary points (not for BZ integration!!!)
      62             : c    ikzero   : 0 no shift of k-points; 
      63             : c               1 shift of k-points for better use of sym in irrBZ
      64             : c    kzero    : shifting vector to bring one k-point to (0,0,0)(for even nmop)
      65             : c                                             away from (0,0,0)(for odd nmop)  
      66             : c    nmop     : integer number triple: nmop(i), i=1,3; nmop(i)
      67             : c               determines number of k-points in direction of rltv(ix,i)
      68             : c
      69             : c    OUTPUT: k-point set
      70             : c    nkpt     : number of k-points generated in set
      71             : c    vkxyz    : vector of kpoint generated; in cartesian representation
      72             : c    wghtkp   : weight associated with k-points for BZ integration
      73             : c    divis    : integer triple divis(i); i=1,4.
      74             : c               Used to find more accurate representation of k-points
      75             : c               vklmn(i,kpt)/divis(i) and weights as wght(kpt)/divis(4)
      76             : c    nkstar   : number of stars for k-points generated in full stars
      77             : c
      78             : c-----------------------------------------------------------------------
      79             :       USE m_constants, ONLY : pimach
      80             :       USE m_ordstar
      81             :       USE m_fulstar
      82             :       IMPLICIT NONE
      83             : C
      84             : C-----> PARAMETER STATEMENTS
      85             : C
      86             :       INTEGER, INTENT (IN) :: mkpt,mface,mdir
      87             : c
      88             : c ---> file number for read and write
      89             : c
      90             :       INTEGER, INTENT (IN) :: iofile,iokpt
      91             : c
      92             : c ---> running mode parameter
      93             : c
      94             :       INTEGER, INTENT (IN) :: kpri,ktest
      95             : C
      96             : C----->  Symmetry information
      97             : C
      98             :       INTEGER, INTENT (IN) :: nsym,idsyst,idtype
      99             :       REAL,    INTENT (IN) :: ccr(3,3,48)
     100             :       REAL,    INTENT (IN) :: rlsymr(3,3,48),talfa(3,48)
     101             : C
     102             : C----->  BRAVAIS LATTICE INFORMATION
     103             : C
     104             :       REAL,    INTENT (IN) :: bltv(3,3),cpoint(3,mface)
     105             : C
     106             : C----->  RECIPROCAL LATTICE INFORMATION
     107             : C
     108             :       INTEGER, INTENT (IN) :: ncorn,nface,nedge
     109             :       REAL,    INTENT (IN) :: xvec(3),rltv(3,3)
     110             :       REAL,    INTENT (IN) :: fnorm(3,mface),fdist(mface)
     111             : C
     112             : C----->  BRILLOUINE ZONE INTEGRATION
     113             : C
     114             :       INTEGER, INTENT (IN) :: nreg,ikzero,nfulst,nbound,idimens
     115             :       INTEGER, INTENT (INOUT) :: nmop(3)
     116             :       REAL,    INTENT (INOUT) :: kzero(3)
     117             :       INTEGER, INTENT (OUT):: nkpt,nkstar
     118             :       REAL,    INTENT (OUT):: vkxyz(3,mkpt),wghtkp(mkpt),divis(4)
     119             : C
     120             : C --->  local variables
     121             : c
     122             :       CHARACTER*80 blank
     123             :       INTEGER  i,idim,i1,i2,i3,ii,ij,ik,is,isym,ifac, iik,iiik
     124             :       INTEGER  ikc, i1red,nred,isumkpt,nleft,nirrbz
     125             :       INTEGER  dirmin,dirmax,ndir1,ndir2,idir
     126             :       INTEGER  kpl,kpm,kpn,nstar(mdir),nstnew
     127             :       INTEGER  iplus,iminus,nc2d,n
     128             :       REAL     invtpi,zero,one,half,eps,eps1,orient,sum,denom,aivnkpt
     129             : 
     130             :       INTEGER  nfract(3),lim(3),isi(3)
     131          50 :       REAL     cp2d(3,mface)
     132             :       REAL     vktra(3),vkstar(3,48),ainvnmop(3),fsig(2),vktes(3)
     133          75 :       INTEGER, ALLOCATABLE :: ikpn(:,:),irrkpn(:),nkrep(:)
     134          50 :       INTEGER, ALLOCATABLE :: iside(:),iostar(:)
     135          50 :       REAL,    ALLOCATABLE :: fract(:,:), vkrep(:,:)
     136             : C
     137             : C --->  intrinsic functions
     138             : c
     139             :       INTRINSIC   real,abs
     140             : C
     141             : C --->  save and data statements
     142             : c
     143             :       SAVE     one,zero,half,eps,eps1,iplus,iminus
     144             :       DATA     zero/0.00/,one/1.00/,half/0.50/,
     145             :      +         eps/1.0e-8/,eps1/1.0e-6/,iplus/1/,iminus/-1/
     146             : c
     147             : c-----------------------------------------------------------------------
     148             : c
     149          25 :       ALLOCATE (fract(mkpt,3),vkrep(3,mkpt),ikpn(48,mkpt),irrkpn(mkpt))
     150          25 :       ALLOCATE (nkrep(mkpt),iostar(mkpt),iside(mface))
     151          25 :       if (kpri .ge. 1) then
     152             : c       write(iofile,'(/)')
     153           0 :         write(iofile,'(3x,'' *<* kptmop *>* '')')
     154           0 :         write(iofile,'(3x,'' generate k-vectors'')')
     155           0 :         write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~'')')
     156           0 :         write(iofile,'(3x,'' by Monkhorst-Pack-method'')')
     157           0 :         write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~'')')
     158           0 :         if (idimens .eq. 2) then
     159           0 :             write(iofile,'(3x,'' in 2 dimensions'')')
     160           0 :             write(iofile,'(3x,'' ~~~~~~~~~~~~~~~'')')
     161             :             write(iokpt,'(''#'',i4,20x,'' idimens: '',
     162           0 :      >      ''k-points for 2-dimensional lattice'')') idimens
     163           0 :         else if (idimens .lt. 2 .or. idimens.gt.3) then
     164             :             write(iofile,'(1x,i4,'' idimens wrong, choose 2 or 3'')')
     165           0 :      >                                                       idimens
     166           0 :              CALL juDFT_error("idimens",calledby="kptmop")
     167             :         end if
     168           0 :         if (nreg .eq. 1) then
     169           0 :            if (nfulst .eq. 1) then
     170           0 :             write(iofile,'(3x,'' in 1. Brillouin zone'')')
     171           0 :             write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~'')')
     172           0 :             write(iofile,'(3x,'' full stars generated'')')
     173           0 :             write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~'')')
     174             :             write(iokpt,'(''#'',2(i4,1x),14x,'' nreg,nfulst: '',
     175           0 :      >      ''k-points in totBZ, ordered in full stars'')') nreg,nfulst
     176           0 :            else if (nfulst .eq. 0) then
     177           0 :             write(iofile,'(3x,'' in parallelepiped'')')
     178           0 :             write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~'')')
     179           0 :             write(iofile,'(3x,'' full stars not generated'')')
     180           0 :             write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~'')')
     181             :             write(iokpt,'(''#'',2(i4,1x),14x,'' nreg,nfulst: '',
     182           0 :      >      ''k-points in par-epiped, not in full stars'')') nreg,nfulst
     183             :            else
     184             :             write(iofile,'(2(1x,i4),15x,'' nreg,nfulst: wrong choice '',
     185             :      >       ''of parameters; allowed combinations: (1,1); (1,0)'')' )
     186           0 :      >                                                      nreg,nfulst
     187           0 :               CALL juDFT_error("nfulst",calledby="kptmop")
     188             :            end if
     189           0 :         else if (nreg .eq. 0) then
     190           0 :         write(iofile,'(3x,'' in irred wedge of 1. Brillouin zone'')')
     191           0 :         write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')')
     192           0 :         write(iokpt,'(''#'',i4,21x,''nreg: k-points in irrBZ'')') nreg
     193             :  
     194             :         else
     195           0 :         write(iofile,'(3x,'' wrong choice of nreg: '', i4)') nreg
     196             :  
     197           0 :         write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')')
     198           0 :          CALL juDFT_error("nreg",calledby="kptmop")
     199             :         end if
     200             :         write(iokpt,'(''#'',3(i4,1x),10x,''nmop(1), nmop(2), nmop(3)'',
     201           0 :      +     '' as read in'')')                      (nmop(i), i=1,3)
     202             :         write(iokpt,'(''#'',2(i4,1x),14x,'' ikzero, nbound'')')
     203           0 :      +                                               ikzero, nbound
     204           0 :         write(iofile,'(/)')
     205           0 :         if (nbound .eq. 1) then
     206           0 :         write(iofile,'(3x,'' k-points on boundary included'')')
     207           0 :         write(iofile,'(3x,'' irregular Monkhorst-Pack-ratios'')')
     208           0 :         write(iofile,'(3x,'' cannot be used for BZ-integration'')')
     209           0 :         else if (nbound.eq. 0) then
     210           0 :         write(iofile,'(3x,'' no k-points on boundary of BZ'')')
     211           0 :         write(iofile,'(3x,'' regular Monkhorst-Pack-ratios'')')
     212           0 :         write(iofile,'(3x,'' can be used for BZ-integration'')')
     213             :         else
     214           0 :         write(iofile,'(3x,'' wrong choice of nbound: '', i4)') nbound
     215           0 :         write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')')
     216           0 :          CALL juDFT_error("nbound",calledby="kptmop")
     217             :         end if
     218             : c
     219           0 :         write(iofile,'(1x,i4,10x,''iofile'')') iofile
     220           0 :         write(iofile,'(1x,i4,10x,''iokpt'')')  iokpt
     221           0 :         write(iofile,'(1x,i4,10x,''kpri'')')  kpri
     222           0 :         write(iofile,'(1x,i4,10x,''ktest'')')  ktest
     223           0 :         write(iofile,'(1x,3(f10.7,1x),10x,''xvec'')') (xvec(ii),ii=1,3)
     224           0 :         write(iofile,'(1x,i4,10x,''ncorn'')')  ncorn
     225           0 :         write(iofile,'(1x,i4,10x,''nedge'')')  nedge
     226           0 :         write(iofile,'(1x,i4,10x,''nface'')')  nface
     227           0 :         do  5 ifac = 1,nface
     228             :         write(iofile,'(1x,i4,1x,3(f10.7,1x),10x,''fnorm'')')
     229           0 :      +                 ifac,(fnorm(ii,ifac),ii=1,3)
     230             :         write(iofile,'(1x,i4,1x,f10.7,1x,10x,''fdist'')')
     231           0 :      +                 ifac,fdist(ifac)
     232           0 :   5     continue
     233             : c
     234             :       end if
     235             : !
     236             : ! --->   for 2 dimensions only the following Bravais lattices exist:
     237             : !          TYPE                    EQUIVALENT 3-DIM        idsyst/idtype
     238             : !         square               = p-tetragonal ( 1+2 axis )      2/1
     239             : !         rectangular          = p-orthorhomb ( 1+2 axis )      3/1
     240             : !         centered rectangular = c-face-orthorhomb( 1+2 axis)   3/6
     241             : !         hexagonal            = p-hexagonal  ( 1+2 axis )      4/1
     242             : !         oblique              = p-monoclinic ( 1+2 axis )      6/1
     243             : !
     244          25 :       IF (idimens .EQ. 2) THEN
     245             : !
     246             : ! --->   identify the allowed symmetries
     247             : !        and check the consistency of the Monkhorst-Pack-parameters
     248             : !
     249           6 :         IF (idsyst.EQ.2 .OR. idsyst.EQ.4) THEN
     250           6 :           IF (idtype.EQ.1) THEN
     251           6 :              IF (nmop(3).EQ.1) nmop(3) = 0
     252           6 :              IF (nmop(1).NE.nmop(2) .OR. nmop(3).NE.0) THEN
     253           0 :                 nmop(2) = nmop(1)
     254           0 :                 nmop(3) = 0
     255             :                 WRITE (iofile,'(1x,''WARNING!!!!!!!'',/,
     256             :      +            ''nmop-Parameters not in accordance with symmetry'',/,
     257             :      +            2(1x,i4),/,
     258             :      +            '' we have set nmop(2) = nmop(1)'',/,
     259           0 :      +            '' and/or nmop(3) = 0'')') idsyst, idtype
     260             :                 WRITE (iofile,'(3(1x,i4),'' new val for nmop: '')')
     261           0 :      +                                              (nmop(i),i=1,3)
     262             :                 CALL juDFT_warn(
     263             :      +             "k point mesh not compatible with symmetry (1)",
     264           0 :      +             calledby='kptmop')
     265             :              ELSE
     266           6 :                 WRITE (iofile,'('' values accepted unchanged'')')
     267             :                 WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     268          30 :      +                                            (nmop(i),i=1,3)
     269             :              ENDIF
     270             :           ENDIF
     271           0 :         ELSEIF (idsyst.EQ.3) THEN
     272           0 :           IF (idtype.EQ.1 .OR. idtype.EQ.6) THEN
     273           0 :              IF (nmop(3).EQ.1) nmop(3) = 0
     274           0 :              IF (nmop(3).NE.0) THEN
     275           0 :                 nmop(3) = 0
     276             :                 WRITE (iofile,'(1x,''WARNING!!!!!!!'',/,
     277             :      +            ''nmop-Parameters not in accordance with symmetry'',/,
     278             :      +            2(1x,i4),/,
     279           0 :      +            '' we have set nmop(3) = 0'')') idsyst, idtype
     280             :                 WRITE (iofile,'(3(1x,i4),'' new val for nmop: '')')
     281           0 :      +                                              (nmop(i),i=1,3)
     282             :                 CALL juDFT_warn(
     283             :      +             "k point mesh not compatible with symmetry (2)",
     284           0 :      +             calledby='kptmop')
     285             :              ELSE
     286           0 :                 WRITE (iofile,'('' values accepted unchanged'')')
     287             :                 WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     288           0 :      +                                            (nmop(i),i=1,3)
     289             :              ENDIF
     290             :           ENDIF
     291           0 :         ELSEIF (idsyst.EQ.6) THEN
     292           0 :           IF (idtype.EQ.1) THEN
     293           0 :              IF (nmop(3).EQ.1) nmop(3) = 0
     294           0 :              IF (nmop(3).NE.0) THEN
     295           0 :                 nmop(3) = 0
     296             :                 WRITE (iofile,'(1x,''WARNING!!!!!!!'',/,
     297             :      +            ''nmop-Parameters not in accordance with symmetry'',/,
     298             :      +            2(1x,i4),/,
     299           0 :      +            '' we have set nmop(3) = 0'')') idsyst, idtype
     300             :                 WRITE (iofile,'(3(1x,i4),'' new val for nmop: '')')
     301           0 :      +                                              (nmop(i),i=1,3)
     302             :                 CALL juDFT_warn(
     303             :      +             "k point mesh not compatible with symmetry (3)",
     304           0 :      +             calledby='kptmop')
     305             :              ELSE
     306           0 :                 WRITE (iofile,'('' values accepted unchanged'')')
     307             :                 WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     308           0 :      +                                             (nmop(i),i=1,3)
     309             :              ENDIF
     310             :           ENDIF
     311             :         ELSE
     312             : !
     313             : ! --->   in all other cases:
     314             : !
     315             :           WRITE (iofile,'(3(1x,i4),20x,'' idimens,idsyst,idtype: '',
     316             :      >      ''wrong choice for 2-dimensional crystal structure'')')
     317           0 :      >                                 idimens,idsyst,idtype
     318           0 :            CALL juDFT_error("2-dim crystal",calledby="kptmop")
     319             :         ENDIF 
     320             : !
     321             : ! --->   check consistency of nmop-parameters with crystal symmetry
     322             : !
     323          19 :       ELSEIF (idimens.EQ.3) THEN
     324          19 :          IF (idsyst.EQ.1 .OR. idsyst.EQ.5) THEN
     325             :              IF (nmop(1).NE.nmop(2) .OR. nmop(1).NE.nmop(3)
     326           6 :      +                              .OR. nmop(2).NE.nmop(3)) THEN
     327           0 :                 nmop(3) = nmop(1)
     328           0 :                 nmop(2) = nmop(1)
     329             :                 WRITE (iofile,'(1x,''WARNING!!!!!!!'',/,
     330             :      +          ''nmop-Parameters not in accordance with symmetry'',/,
     331             :      +          2(1x,i4),/,
     332           0 :      +          '' we have set all nmop(i) = nmop(1)'')') idsyst, idtype
     333             :                 WRITE (iofile,'(3(1x,i4),'' new val for nmop(i): '')')
     334           0 :      +                                                 (nmop(i),i=1,3)
     335             :                 CALL juDFT_warn(
     336             :      +             "k point mesh not compatible with symmetry (4)",
     337           0 :      +             calledby='kptmop')
     338             :              ELSE
     339           6 :                 WRITE (iofile,'('' values accepted unchanged'')')
     340             :                 WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     341          30 :      +                                                 (nmop(i),i=1,3)
     342             :              ENDIF
     343          13 :          ELSEIF (idsyst.EQ.2 .OR. idsyst.eq.4) THEN
     344          13 :              if((nmop(3).eq.nmop(2)).and.idsyst.eq.2)then
     345           4 :                 WRITE (iofile,'('' values accepted unchanged'')')
     346             :                 WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     347           4 :      +                                            (nmop(i),i=1,3)
     348           9 :              elseif(nmop(1).NE.nmop(2)) THEN
     349           0 :                 nmop(2) = nmop(1)
     350             :                 WRITE (iofile,'(1x,''WARNING!!!!!!!'',/,
     351             :      +            ''nmop-Parameters not in accordance with symmetry'',/,
     352             :      +            2(1x,i4),/,
     353           0 :      +            '' we have set nmop(2) = nmop(1)'')') idsyst, idtype
     354             :                 WRITE (iofile,'(3(1x,i4),'' new val for nmop: '')')
     355           0 :      +                                              (nmop(i),i=1,3)
     356             :                 CALL juDFT_warn(
     357             :      +             "k point mesh not compatible with symmetry (5)",
     358           0 :      +             calledby='kptmop')
     359             :              ELSE
     360           9 :                 WRITE (iofile,'('' values accepted unchanged'')')
     361             :                 WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')')
     362           9 :      +                                            (nmop(i),i=1,3)
     363             :              ENDIF
     364           0 :          ELSEIF (idsyst.LT.1 .OR. idsyst.GT.7) THEN
     365             :              WRITE (iofile,'(1x,''wrong choice of symmetry'',/,
     366           0 :      +                               2(1x,i4))') idsyst, idtype
     367           0 :              WRITE (iofile,'(''only values 1.le.idsyst.le.7 allowed'')')
     368           0 :               CALL juDFT_error("wrong idsyst",calledby="kptmop")
     369             :          ELSE
     370           0 :              WRITE (iofile,'('' values accepted unchanged'')')
     371             :              WRITE (iofile,'(3(1x,i4),11x,''nmop(i),i=1,3'')')
     372           0 :      +                                          (nmop(i),i=1,3)
     373             :          ENDIF
     374             :       ELSE
     375           0 :           CALL juDFT_error("idimens =/= 2,3 ",calledby="kptmop")
     376             :       ENDIF
     377             : !
     378             : ! --->   start calculation
     379             : ! =====================================================================
     380             : !
     381             : ! ---> set sign constants
     382          25 :        isi(1) = 0
     383          25 :        isi(2) = iminus
     384          25 :        isi(3) = iplus
     385             : !
     386             : ! ---> calc orientation of boundary faces of irr wedge of BZ
     387             : !      characterized by
     388             : !          iside(i)= sign( (xvec,fnorm(i))-fdist(i) ) ;(i=1,nface )
     389             : !
     390          25 :       WRITE (iofile,'(1x,''orientation of boundary faces'')')
     391         170 :       DO ifac = 1, nface
     392         145 :          orient = zero
     393         145 :          iside(ifac) = iplus
     394         580 :          DO ii = 1, 3
     395         580 :             orient = orient + xvec(ii)*fnorm(ii,ifac)
     396             :          ENDDO
     397         145 :          orient = orient - fdist(ifac)
     398         145 :          IF (orient .LT. 0) iside(ifac) = iminus
     399             :          WRITE (iofile,'(1x,2(i4,2x),f10.7,10x,''ifac,iside,orient'',
     400         170 :      +                       '' for xvec'')') ifac,iside(ifac),orient
     401             :       ENDDO
     402             : 
     403          25 :       invtpi = one / ( 2.0 * pimach() )
     404             : 
     405          25 :       WRITE (iofile,'(''Bravais lattice vectors'')' )
     406         100 :       DO ii = 1, 3
     407         100 :          WRITE (iofile,'(43x,3(1x,f11.6))') (bltv(ii,ikc), ikc=1,3)
     408             :       ENDDO
     409          25 :       WRITE (iofile,'(''reciprocal lattice vectors'')' )
     410         100 :       DO ii = 1, 3
     411         100 :          WRITE (iofile,'(43x,3(1x,f11.6))' ) (rltv(ii,ikc), ikc=1,3)
     412             :       ENDDO
     413             : !
     414             : ! ---> nmop(i) are Monkhorst-Pack parameters; they determine the
     415             : !                           number of k-points in i-direction
     416             : !      if basis vector lengths are not related by symmetry,
     417             : !      we can use independent fractions for each direction
     418             : !
     419             :       WRITE (iofile,'(3(1x,i4),10x,'' Monkhorst-Pack-parameters'')')
     420          25 :      +                                             (nmop(i1),i1=1,3)
     421             : 
     422          94 :       DO idim = 1, idimens
     423          94 :          IF (nmop(idim).GT.0) THEN
     424          69 :             ainvnmop(idim) = one/ real(nmop(idim))
     425             :          ELSE
     426             :             WRITE (iofile,'('' nmop('',i4,'') ='',i4,
     427           0 :      +                     '' not allowed'')') idim, nmop(idim)
     428           0 :              CALL juDFT_error("nmop wrong",calledby="kptmop")
     429             :          ENDIF
     430             :       ENDDO
     431             : !
     432             : ! --->nreg determines region of BZ in which k-points are generated
     433             : !
     434          25 :       IF ( nreg .EQ. 1) THEN
     435           0 :         IF (nfulst .EQ. 1) THEN
     436             :           WRITE (iofile,'(2(1x,i4),15x,'' nreg,nfulst: '',
     437           0 :      >      ''k-points in totBZ, ordered in full stars'')') nreg,nfulst
     438           0 :         ELSEIF (nfulst .EQ. 0) then
     439             :           WRITE (iofile,'(2(1x,i4),15x,'' nreg,nfulst: '',
     440           0 :      >     ''k-points in par-epiped, not in full stars'')') nreg,nfulst
     441             :         ENDIF
     442          25 :       ELSEIF ( nreg .EQ. 0) THEN
     443             :         WRITE (iofile,'(1x,i4,10x,''nreg; k-points in'',
     444          25 :      +                '' irreducible wedge of BZ'')' ) nreg
     445             :       ENDIF
     446             : 
     447          25 :       WRITE (iofile,'(1x,''Monkhorst-Pack-fractions'')' )
     448             : !
     449             : ! ---> nbound=1: k-points are generated on boundary of BZ
     450             : !        include  fract(1) =       -1/2
     451             : !            and  fract(2*nmop+1) = 1/2     for surface points of BZ
     452             : !
     453          25 :       IF ( nbound .EQ. 1) THEN
     454             :         WRITE (iofile,'(1x,i4,10x,''nbound; k-points on boundary'',
     455           0 :      +   '' of BZ included'')' ) nbound
     456             : !
     457             : ! ---> irregular Monkhorst--Pack--fractions
     458             : !                                   fract(r) = r / (2*nmop)
     459             : !
     460           0 :         DO idim = 1,idimens
     461           0 :            denom = half*ainvnmop(idim)
     462           0 :            divis(idim) = one / denom
     463             : 
     464           0 :            DO kpn = -nmop(idim),nmop(idim)
     465           0 :              fract(kpn+nmop(idim)+1,idim) = denom * real (kpn)
     466           0 :              WRITE (iofile,'(10x,f10.7)' ) fract(kpn+nmop(idim)+1,idim)
     467             :            ENDDO
     468           0 :            nfract(idim) = 2*nmop(idim) + 1
     469             :         ENDDO
     470           0 :         IF (idimens .eq. 2) THEN
     471           0 :            nfract(3) = 1
     472           0 :            fract(1,3) = 0
     473           0 :            divis(3) = one
     474             :         END IF
     475             : !
     476             : ! ---> nbound=0: k-points are NOT generated on boundary of BZ
     477             : !                This is the regular Monkhorst-Pack-method
     478             : !
     479          25 :       ELSEIF ( nbound .eq. 0) then
     480             :          WRITE (iofile,'(1x,i4,10x,''nbound; no k-points '',
     481          25 :      +                 '' on boundary of BZ'')' ) nbound
     482             : !
     483             : ! --->   regular Monkhorst--Pack--fractions
     484             : !                                   fract(r) =(2*r-nmop-1) / (2*nmop)
     485             : !
     486          94 :          DO idim = 1,idimens
     487          69 :             denom = half*ainvnmop(idim)
     488          69 :             divis(idim) = one / denom
     489          69 :             WRITE(iofile,'(5x,i4,5x,''idim'')' ) idim
     490         333 :             DO  kpn = 1,nmop(idim)
     491         264 :               fract(kpn,idim) = denom * real (2*kpn -nmop(idim)-1)
     492         333 :               write(iofile,'(10x,f10.7)' ) fract(kpn,idim)
     493             :             ENDDO
     494          94 :             nfract(idim) = nmop(idim)
     495             :          ENDDO
     496             : 
     497          25 :          IF (idimens .EQ. 2) THEN
     498           6 :             nfract(3) = 1
     499           6 :             fract(1,3) = 0
     500           6 :             divis(3) = one
     501             :          ENDIF
     502             : 
     503             :       ELSE
     504           0 :          WRITE (iofile,'(3x,'' wrong choice of nbound:'', i4)') nbound
     505           0 :          WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')')
     506           0 :           CALL juDFT_error("nbound",calledby="kptmop")
     507             :       ENDIF
     508             : !
     509             : ! ---> set kzero
     510             : !
     511          25 :       kzero(1) = zero
     512          25 :       kzero(2) = zero
     513          25 :       kzero(3) = zero
     514             : 
     515          25 :       IF (ikzero.NE.0) THEN
     516             : !
     517             : ! --->   set kzero .ne. zero
     518             : !        for nmop even and non-orthogonal  rec lattice vectors
     519             : !          shift all k-points to bring one to (0,0,0) 
     520             : !        for nmop odd and non-orthogonal  rec lattice vectors
     521             : !          shift k-points to avoid (0,0,0) 
     522             : !            implemented for hexagon   (idsyst =4, idtype arbitrary)
     523             : !                        and ort(2-dim)(idsyst =3, idtype =1)
     524             : !                        and fcc       (idsyst =1, idtype =3)
     525             : !                        and p-monocli (idsyst =6, idtype =1)
     526             : !
     527           0 :         denom = half*ainvnmop(1)
     528             : 
     529           0 :         IF( idsyst .EQ. 4) THEN    !  hexagonal Bravais lattice
     530             : 
     531           0 :           DO ii = 1, 3
     532           0 :             kzero(ii) = denom * (rltv(ii,1)+rltv(ii,2))
     533             :           ENDDO
     534             : 
     535           0 :         ELSEIF ( idsyst .EQ. 3 .and. idtype .EQ. 1) THEN ! orthorhombic (2D)
     536             : 
     537           0 :           IF (nmop(1).EQ.1) THEN
     538           0 :               denom = half*ainvnmop(3)
     539           0 :             DO ii = 1, 3
     540           0 :               kzero(ii) = denom * (rltv(ii,2)+rltv(ii,3))
     541             :             ENDDO
     542             :           ELSE
     543             :             WRITE (iofile,'(''ikzero.ne.0 is NOT applied '',
     544           0 :      >                     ''for present choice of parameters'')')
     545             :             WRITE (iofile,'(5(1x,i4),
     546             :      >               ''; idsyst,idtype,nmop(1),nmop(2),nmop(3)'')')
     547           0 :      >                       idsyst,idtype,nmop(1),nmop(2),nmop(3)
     548             :           ENDIF
     549             : 
     550           0 :         ELSEIF ( idsyst .eq. 1 .and. idtype .eq. 3) then !  face centered cubic
     551             : 
     552           0 :           DO ii = 1, 3
     553           0 :             kzero(ii) = denom * (rltv(ii,1)+rltv(ii,2)+rltv(ii,3))
     554             :           ENDDO
     555             : 
     556           0 :         ELSEIF ( idsyst .eq. 6 .and. idtype .eq. 1) then ! p-monoclinic 
     557             : 
     558           0 :           DO ii = 1, 3
     559             :             kzero(ii) = half*ainvnmop(1) * rltv(ii,1)
     560           0 :      +                + half*ainvnmop(2) * rltv(ii,2)
     561             :           ENDDO
     562             : 
     563             :         ELSE
     564             :           WRITE (iofile,'(''ikzero.ne.0 is NOT applied '',
     565           0 :      >                   ''for this system'')')
     566             :           WRITE (iofile,'(5(1x,i4),
     567             :      >             ''; idsyst,idtype,nmop(1),nmop(2),nmop(3)'')')
     568           0 :      >                     idsyst,idtype,nmop(1),nmop(2),nmop(3)
     569             :         ENDIF
     570             :       ENDIF
     571             : !
     572             : !
     573             : ! --->   initialize k-points = zero and weights = 1.0
     574             : !
     575        6061 :       DO  kpn = 1,mkpt
     576        6036 :          vkxyz(1,kpn) = zero
     577        6036 :          vkxyz(2,kpn) = zero
     578        6036 :          vkxyz(3,kpn) = zero
     579        6061 :          wghtkp(kpn) = one
     580             :       ENDDO
     581             : !
     582             : ! ---> generate equidistant k-vectors in cartesian coordinates
     583             : !                                   with off-set kzero
     584             : !
     585          25 :       nkpt = 0
     586          95 :       DO i3 = 1,nfract(3)
     587         597 :          DO i2 = 1,nfract(2)
     588        6608 :             DO i1 = 1,nfract(1)
     589        6036 :               nkpt = nkpt + 1
     590        6036 :               IF (nkpt>mkpt)  CALL juDFT_error("nkpt > mkpt",calledby
     591           0 :      +             ="kptmop")
     592             :               vkxyz(1,nkpt) = kzero(1) + rltv(1,1)*fract(i1,1)
     593             :      +                                 + rltv(1,2)*fract(i2,2)
     594        6036 :      +                                 + rltv(1,3)*fract(i3,3)
     595             :               vkxyz(2,nkpt) = kzero(2) + rltv(2,1)*fract(i1,1)
     596             :      +                                 + rltv(2,2)*fract(i2,2)
     597        6036 :      +                                 + rltv(2,3)*fract(i3,3)
     598             :               vkxyz(3,nkpt) = kzero(3) + rltv(3,1)*fract(i1,1)
     599             :      +                                 + rltv(3,2)*fract(i2,2)
     600        6538 :      +                                 + rltv(3,3)*fract(i3,3)
     601             :             ENDDO
     602             :          ENDDO
     603             :       ENDDO
     604             : !
     605             : ! --->   calculate weights of k-points and print out k-points
     606             : !         wghtkp = 1/nkpt
     607             : !              ( = 1/(nmop(1)*nmop(2)*nmop(3)) for reg Monk-Pack-method)
     608             : !
     609          25 :       divis(4) = real(nkpt)
     610          25 :       aivnkpt  = one/real(nkpt)
     611             : 
     612        6061 :       DO  kpn= 1,nkpt
     613        6061 :         wghtkp(kpn) = wghtkp(kpn)*aivnkpt
     614             :       ENDDO
     615             : 
     616          25 :       IF (kpri .ge. 1) THEN
     617           0 :       WRITE (iofile,'(''generated k-points vkxyz'')')
     618           0 :       DO  kpn= 1,nkpt
     619             :         WRITE (iofile,'(1x,i4,1x,4(f13.10,1x),10x,''vkxyz, wghtkp'')')
     620           0 :      +               kpn,(vkxyz(ii,kpn),ii=1,3),wghtkp(kpn)
     621             :       ENDDO
     622             :       ENDIF
     623             : 
     624             : !
     625             : ! ---> for nreg=1 and nfulst=0
     626             : !                 (k-points in total zone of arbitrary shape)
     627             : !       k-point generation finished
     628             : !       (k-points are actually generated in parallel-epiped
     629             : !       (r1,r2,r3) which is an equivalent unit cell)
     630             : !       The generated set does not contain all points of the
     631             : !       symmetry related stars.
     632             : !
     633          25 :       IF (nreg.EQ.1 .AND. nfulst.EQ.0) GOTO 1000
     634             : !
     635             : ! ====================================================================
     636             : !
     637             : ! --->   order generated k-points in stars by applying symmetry:
     638             : !        - determine number of different stars nkstar .le. nkpt
     639             : !        - determine order of star iostar(kpn) .le. nsym
     640             : !        - assign pointer ikpn(i,ik); i=1,iostar(ik); ik=1,nkstar
     641             : !        - determine representative vector in irrBZ for each star:
     642             : !                                   vkrep(ix,ik); ix=1,3; ik=1,nkstar
     643             : !
     644             :         CALL ordstar(
     645             :      >               iokpt,kpri,ktest,
     646             :      >               fnorm,fdist,nface,iside,
     647             :      >               nsym,ccr,rltv,mkpt,mface,mdir,
     648             :      =               nkpt,vkxyz,
     649          25 :      <               nkstar,iostar,ikpn,vkrep,nkrep)
     650             : !
     651             : ! ---> for nreg=0 (k-points in irr part of BZ)
     652             : !
     653             : !      (a) calculate weights for k-points in irrBZ
     654             : !           - wghtkp(ik)=iostar(ik)/nkpt_old ; ik=1,nkstar
     655             : !
     656         347 :         DO ik = 1, nkstar
     657         347 :              wghtkp(ik) = wghtkp(ik)*iostar(ik)
     658             :         ENDDO
     659             : !
     660             : !      (b) final preparation of k-points for transfer to file
     661             : !           - assign nkpt= nkstar
     662             : !           - assign vkxyz(ix,ik) = vkrep(ix,ik); ix=1,3; ik=1,nkstar
     663             : !
     664          25 :         IF (nreg.EQ.0) THEN
     665             : 
     666         175 :           DO i1 = 1,3
     667        2032 :             DO ik = 1,nkstar
     668        1041 :               vkxyz(i1,ik) = vkrep(i1,ik)
     669             :             ENDDO
     670             :           ENDDO
     671          25 :           nkpt = nkstar
     672          25 :           GOTO 1000
     673             :         ENDIF
     674             : c
     675             : c ---> for nreg.eq.0: k-point generation finished
     676             : c
     677             : c ================================================================
     678             : c
     679             : c
     680           0 :        IF (nreg.EQ.1 .AND. nfulst.EQ.1) THEN
     681             : c
     682             : c --->   generate full stars for all representative k-points
     683             : c        - for nreg=1 and nfulst=1:
     684             : c              - determine order of full star ifstar(kpn).le.nsym
     685             : c              - assign nkpt= sum {ifstar(ik)} (ik=1,nkstar)
     686             : c              - assign vkxyz(ix,kpn) = vkstar(ix,ikpn(is,ik));
     687             : c                     ix=1,3; kpn=1,nkpt; ik=1,nstar; is=1,ifstar(ik)
     688             : c              - calculate wghtkp(kpn)=wghtkp_old(ik)/(ifstar(ik)
     689             : c                                kpn=1,nkpt; ik=1,nstar
     690             : c
     691             : c
     692             :             CALL fulstar(
     693             :      >                   iofile,iokpt,kpri,ktest,
     694             :      >                   ccr,nsym,
     695             :      >                   vkrep,nkstar,mkpt,mface,mdir,
     696           0 :      =                   nkpt,vkxyz,wghtkp)
     697             : 
     698           0 :             divis(4) = divis(4) * nsym
     699             : 
     700             :       ENDIF
     701             : 
     702             :  1000 CONTINUE
     703             : !      
     704             : ! --> check for corner points, include them into k-point set:
     705             : !
     706          25 :       IF (nbound.EQ.1) THEN
     707           0 :         n = 1
     708           0 :         nc2d = 1               ! determine 2D corner points
     709           0 :         cp2d(:,nc2d) = cpoint(:,n)
     710           0 :         corn: DO n = 2, ncorn
     711           0 :           DO i = 1, n-1
     712           0 :             IF ((abs(cpoint(1,n)-cpoint(1,i)).LT.0.0001).AND.
     713           0 :      +          (abs(cpoint(2,n)-cpoint(2,i)).LT.0.0001)) CYCLE corn
     714             :           ENDDO
     715           0 :           nc2d = nc2d + 1
     716           0 :           cp2d(:,nc2d) = cpoint(:,n)
     717             :         ENDDO corn
     718           0 :         WRITE (iofile,'(''2D corner points in internal units'')')
     719           0 :         corn2d: DO n = 1, nc2d 
     720           0 :           WRITE (iofile,'(i3,3x,2(f10.7,1x))') n,cp2d(1,n),cp2d(2,n)
     721           0 :           DO i = 1, nkpt
     722           0 :             IF ((abs(cp2d(1,n)-vkxyz(1,i)).LT.0.0001).AND.
     723           0 :      +          (abs(cp2d(2,n)-vkxyz(2,i)).LT.0.0001)) CYCLE corn2d
     724             :           ENDDO
     725           0 :           nkpt = nkpt + 1
     726           0 :           vkxyz(:,nkpt) = cp2d(:,n)
     727             :         ENDDO corn2d
     728             :       ENDIF 
     729             : !
     730             : ! --->   print out k-points and weights
     731             : !
     732          25 :       IF (ktest .GE. 1) THEN
     733           0 :         WRITE (iofile,'(''generated k-points vkxyz'')')
     734             :         WRITE (iofile,'(1x,i4,20x,
     735           0 :      +               ''nkpt,number of generated k-points'')') nkpt
     736             :  
     737           0 :        DO kpn = 1, nkpt
     738             :          WRITE (iofile,'(3(f10.7,1x),f12.10,1x,i4,3x,
     739           0 :      +   ''vkxyz, wghtkp'')') (vkxyz(ii,kpn),ii=1,3),wghtkp(kpn), kpn
     740             :        ENDDO
     741             : 
     742             :       ENDIF
     743          25 :       DEALLOCATE (fract,vkrep,ikpn,irrkpn,nkrep,iostar,iside)
     744             : 
     745          25 :       RETURN
     746             :       END SUBROUTINE kptmop
     747             :       END MODULE m_kptmop

Generated by: LCOV version 1.13