LCOV - code coverage report
Current view: top level - init - strgn.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 338 371 91.1 %
Date: 2019-09-08 04:53:50 Functions: 3 3 100.0 %

          Line data    Source code
       1             : MODULE m_strgn
       2             :   USE m_juDFT
       3             :   !
       4             :   !     *********************************************************
       5             :   !     generate two- and three-dimensional stars
       6             :   !     for slab geometry
       7             :   !     e. wimmer   nov.1984    c.l.fu  1987
       8             :   !     implementation of new box-dimension: to treat nonorthogonal
       9             :   !     lattice systems
      10             :   !     S. Bl"ugel, IFF, 17.Nov.97
      11             :   !
      12             :   !     OpenMP paralleliation added
      13             :   !     U.Alekseeva          Jan.2019
      14             :   !     *********************************************************
      15             : CONTAINS
      16          18 :   SUBROUTINE strgn1(&
      17             :        &                  stars,sym,atoms,&
      18             :        &                  vacuum,sphhar,input,cell,xcpot)
      19             :     !
      20             :     USE m_spgrot
      21             :     USE m_angle
      22             :     USE m_types
      23             :     USE m_boxdim
      24             :     USE m_sort
      25             :     USE m_cdn_io
      26             :     IMPLICIT NONE
      27             :     TYPE(t_stars),INTENT(INOUT)  :: stars
      28             :     TYPE(t_sym),INTENT(IN)       :: sym
      29             :     TYPE(t_atoms),INTENT(IN)     :: atoms
      30             :     TYPE(t_vacuum),INTENT(IN)    :: vacuum
      31             :     TYPE(t_sphhar),INTENT(IN)    :: sphhar
      32             :     TYPE(t_input),INTENT(IN)     :: input
      33             :     TYPE(t_cell),INTENT(IN)      :: cell
      34             :     CLASS(t_xcpot),INTENT(IN)    :: xcpot
      35             : 
      36             :     !     ..
      37             :     !     ..
      38             :     !     .. Local Scalars ..
      39             :     REAL arltv1,arltv2,arltv3,s
      40             :     REAL gmi,gla,eps 
      41             :     REAL gfx,gfy,pon,pon2
      42             :     INTEGER j,k,k1,k2,k3,m0,mxx1,mxx2,n
      43             :     INTEGER ned1,nint,kdone,i,i_sym,n_sym
      44             :     LOGICAL NEW,l_cdn1,l_xcExtended, l_error
      45             :     INTEGER kfx,kfy,kfz,kidx,nfftx,nffty,nfftz,kfft
      46             :     INTEGER nfftxy,norm,n1,kidx2,k2i
      47             :     !     ..
      48             :     !     .. Local Arrays ..
      49             :     REAL,    ALLOCATABLE :: gsk3(:)
      50          36 :     INTEGER, ALLOCATABLE :: ig2p(:),INDEX(:),index3(:),kv3rev(:,:)
      51          36 :     REAL g(3),phi3(stars%ng2),phi
      52          36 :     COMPLEX phas(sym%nop)
      53          36 :     INTEGER kr(3,sym%nop),kv(3)
      54             :     INTEGER index2(stars%ng2)
      55             :     !     ..
      56             :     !     ..
      57             :     !
      58             :     !
      59          18 :     nfftx = 3*stars%mx1
      60          18 :     nffty = 3*stars%mx2
      61          18 :     nfftz = 3*stars%mx3
      62          18 :     nfftxy= 9*stars%mx1*stars%mx2
      63             : 
      64          18 :     ALLOCATE (gsk3(stars%ng3),INDEX(stars%ng3),index3(stars%ng3),kv3rev(stars%ng3,3))
      65             : 
      66             :     !
      67             :     !WRITE (*,*) ' stars are always ordered '
      68             : 
      69          18 :     l_xcExtended = xcpot%needs_grad()
      70             :     !--->    read in information if exists
      71          18 :     CALL readStars(stars,l_xcExtended,.TRUE.,l_error)
      72          18 :     IF(.NOT.l_error) THEN
      73             :        GOTO 270
      74             :     END IF
      75             :  
      76          12 :     IF (input%film.AND.sym%invs.AND.(.not.sym%zrfs).AND.(.not.sym%symor)) THEN
      77             :       n_sym = 2 ! needs reordering of 2d-stars
      78             :     ELSE
      79             :       n_sym = 1 ! as before...
      80             :     ENDIF
      81             : 
      82          12 :     mxx1 = 0
      83          12 :     mxx2 = 0
      84          12 :     stars%ng2 = 0
      85          12 :     kv(3) = 0
      86             : 
      87          12 :     stars%kv2 = 0
      88         308 :     DO  k1 = stars%mx1,-stars%mx1,-1
      89         296 :        kv(1) = k1
      90        8280 :        k2_loop:DO  k2 = stars%mx2,-stars%mx2,-1
      91        7972 :           kv(2) = k2
      92             : 
      93       13597 :           DO i_sym = 1, n_sym
      94        7972 :             IF (i_sym == 2) THEN
      95           0 :                kv(1) = - kv(1) ; kv(2) = - kv(2)
      96             :             ENDIF
      97             : 
      98       39860 :             DO j = 1,2
      99       23916 :                g(j) = kv(1)*cell%bmat(1,j) + kv(2)*cell%bmat(2,j)
     100             :             ENDDO
     101        7972 :             s = SQRT(g(1)**2+g(2)**2)
     102             :             !--->   determine the angle of the G_{||} needed in odim calculations
     103             :             !+odim YM cute little 'angle' is written by me to determine the phi2
     104        7972 :             phi = angle(g(1),g(2))
     105             :             !-odim
     106             :             !
     107             :             !--->   check if generated vector belongs to a star already
     108             :             !--->   stars should be within the g_max-sphere !   (Oct.97) sbluegel
     109             :             !
     110       13301 :             IF (s.LT.stars%gmax) THEN      
     111             :                CALL spgrot(&
     112             :                     &                     sym%nop,sym%symor,sym%mrot,sym%tau,sym%invtab,&
     113             :                     &                     kv,&
     114        5958 :                     &                     kr)
     115       32976 :                DO n = 1,sym%nop2
     116       27018 :                   IF (mxx1.LT.kr(1,n)) mxx1 = kr(1,n)
     117       32976 :                   IF (mxx2.LT.kr(2,n)) mxx2 = kr(2,n)
     118             :                ENDDO
     119     1992405 :                DO k = 1,stars%ng2 
     120     7163113 :                   DO n = 1,sym%nop2
     121     2587997 :                      IF (kr(1,n).EQ.stars%kv2(1,k) .AND.&
     122     1986447 :                           &                   kr(2,n).EQ.stars%kv2(2,k)) CYCLE k2_loop
     123             :                   ENDDO
     124             :                ENDDO
     125             :                !--->    new representative found
     126        3315 :                stars%ng2 = stars%ng2 + 1
     127             :                IF (stars%ng2.GT.stars%ng2) THEN
     128             :                   WRITE (6,8070) stars%ng2,stars%ng2
     129             :                   CALL juDFT_error("nq2.GT.n2d",calledby="strgn")
     130             :                ENDIF
     131        9945 :                DO j = 1,2
     132        9945 :                   stars%kv2(j,stars%ng2) = kv(j)
     133             :                ENDDO
     134        3315 :                stars%sk2(stars%ng2) = s
     135        3315 :                stars%phi2(stars%ng2) = phi
     136             :             ENDIF
     137             :             ENDDO ! i_sym
     138             : 
     139             :        ENDDO k2_loop
     140             :     ENDDO
     141             : 8070 FORMAT ('nq2 = ',i5,' > n2d =',i5)
     142             : 
     143          12 :     IF( mxx1.GT.stars%mx1) THEN
     144           0 :        WRITE (6,'(/'' mxx1.gt.k1d. mxx1,k1d='',2i4)')  mxx1,stars%mx1
     145           0 :        CALL juDFT_error("mxx1.gt.k1d",calledby="strgn")
     146             :     ENDIF
     147             : 
     148          12 :     IF (mxx2.GT.stars%mx2) THEN
     149           0 :        WRITE (6,'(/'' mxx2.gt.k2d. mxx2,k2d='',2i4)')  mxx2,stars%mx2
     150           0 :        CALL juDFT_error("mxx2.gt.k2d",calledby="strgn")
     151             :     ENDIF
     152             : 
     153             :     !--->    sort for increasing length sk2
     154        3327 :     DO  k = 1,stars%ng2
     155        3327 :        gsk3(k) = (stars%mx1+stars%kv2(1,k)) + (stars%mx2+stars%kv2(2,k))*(2*stars%mx1+1)
     156             :     ENDDO
     157          12 :     CALL sort(INDEX(:stars%ng2),stars%sk2(:stars%ng2),gsk3(:stars%ng2))
     158        3327 :     DO k = 1,stars%ng2
     159        3315 :        kv3rev(k,1) = stars%kv2(1,INDEX(k))
     160        3315 :        kv3rev(k,2) = stars%kv2(2,INDEX(k))
     161        3315 :        gsk3(k) = stars%sk2(INDEX(k))
     162        3327 :        phi3(k) = stars%phi2(INDEX(k))
     163             :     ENDDO
     164        3327 :     DO k = 1,stars%ng2
     165        3315 :        stars%kv2(1,k) = kv3rev(k,1)
     166        3315 :        stars%kv2(2,k) = kv3rev(k,2)
     167        3315 :        stars%sk2(k) = gsk3(k)
     168        3327 :        stars%phi2(k) = phi3(k)
     169             :     ENDDO
     170             : 
     171             : 
     172             :     WRITE (6,'(/'' nq2='',i4/'' k,kv2(1,2), sk2, phi2''&
     173             :          &     /(3i4,f10.5,f10.5))')&
     174          12 :          &  stars%ng2,(k,stars%kv2(1,k),stars%kv2(2,k),stars%sk2(k),stars%phi2(k),k=1,stars%ng2)
     175             :     !
     176             :     !     three dimensional stars
     177             :     !
     178          12 :     stars%ng3 = 0
     179          12 :     stars%ig = 0
     180         482 :     DO k3 = -stars%mx3,stars%mx3
     181       12476 :        DO k2 = -stars%mx2,stars%mx2
     182      351230 :           DO k1 = -stars%mx1,stars%mx1
     183      338766 :              stars%ig(k1,k2,k3) = 0
     184      350760 :              stars%rgphs(k1,k2,k3) = cmplx(0.0,0.0)
     185             :           ENDDO
     186             :        ENDDO
     187             :     ENDDO
     188             :     !+gu
     189          12 :     stars%igfft2(:,:)=0
     190          12 :     stars%igfft(:,:)=0
     191          12 :     stars%pgfft=0.0
     192          12 :     stars%pgfft2=0.0
     193             :     !-gu
     194          12 :     WRITE (6,'(/'' bmat(3,3),mx3='',f10.5,i5)') cell%bmat(3,3),stars%mx3
     195             : 
     196             :     IF (stars%mx3.GT.stars%mx3) THEN
     197             :        WRITE ( 6,FMT=8000) stars%mx3,stars%mx3
     198             :        CALL juDFT_error("mx3.gt.k3d",calledby="strgn")
     199             :     ENDIF
     200             : 8000 FORMAT('   mx3.gt.k3d:',2i6)
     201             : 
     202          12 :     m0 = -stars%mx3
     203             :     !     zrfs,invs: z-reflection, inversion.
     204          12 :     IF (sym%zrfs .OR. sym%invs) m0 = 0
     205             : 
     206          12 :     stars%ig2 = 0
     207          12 :     stars%sk3 = 0.0
     208          12 :     stars%kv3 = 0
     209        3327 :     DO  k2 = 1,stars%ng2
     210      250689 :        DO  k3 = m0,stars%mx3
     211      123681 :           s = SQRT(stars%sk2(k2)**2+ (k3*cell%bmat(3,3))**2)
     212             :           !
     213             :           !--->   stars should be within the g_max-sphere !   (Oct.97) sbluegel
     214      126996 :           IF (s.LT.stars%gmax) THEN      
     215             :              !
     216       78966 :              stars%ng3 = stars%ng3 + 1
     217             :              IF (stars%ng3.GT.stars%ng3) THEN
     218             :                 WRITE (6,*) stars%ng3,stars%ng3
     219             :                 CALL juDFT_error("nq3.GT.n3d",calledby="strgn")
     220             :              ENDIF
     221      236898 :              DO j = 1,2
     222      236898 :                 stars%kv3(j,stars%ng3) = stars%kv2(j,k2)
     223             :              ENDDO
     224       78966 :              stars%kv3(3,stars%ng3) = k3
     225       78966 :              stars%ig2(stars%ng3) = k2
     226       78966 :              stars%sk3(stars%ng3) = s
     227             :           ENDIF
     228             :        ENDDO
     229             :     ENDDO
     230             : 
     231             :     !--->    sort for increasing length sk3
     232             :     ! secondary key for equal length stars
     233       78978 :     DO  k = 1,stars%ng3
     234             :        gsk3(k) = (stars%mx1+stars%kv3(1,k)) +&
     235             :             &           (stars%mx2+stars%kv3(2,k))*(2*stars%mx1+1) +&
     236       78978 :             &           (stars%mx3+stars%kv3(3,k))*(2*stars%mx1+1)*(2*stars%mx2+1)
     237             :     ENDDO
     238          12 :     CALL sort(index(:stars%ng3),stars%sk3,gsk3)
     239             : 
     240          12 :     ALLOCATE (ig2p(stars%ng3))
     241             : 
     242       78978 :     DO k = 1,stars%ng3
     243       78966 :        kv3rev(k,1) = stars%kv3(1,INDEX(k))
     244       78966 :        kv3rev(k,2) = stars%kv3(2,INDEX(k))
     245       78966 :        kv3rev(k,3) = stars%kv3(3,INDEX(k))
     246       78966 :        gsk3(k) = stars%sk3(INDEX(k))
     247       78978 :        ig2p(k) = stars%ig2(INDEX(k))
     248             :     ENDDO
     249       78978 :     DO k = 1,stars%ng3
     250       78966 :        stars%kv3(1,k) = kv3rev(k,1)
     251       78966 :        stars%kv3(2,k) = kv3rev(k,2)
     252       78966 :        stars%kv3(3,k) = kv3rev(k,3)
     253       78966 :        stars%sk3(k) = gsk3(k)
     254       78978 :        stars%ig2(k) = ig2p(k)
     255             :     ENDDO
     256             :   
     257             :     !
     258             :     !--->  determine true gmax and change old gmax to new gmax
     259             :     !
     260          12 :     WRITE (6,8060) stars%gmax, stars%sk3(stars%ng3)
     261          12 :     stars%gmax = stars%sk3(stars%ng3)
     262             : 8060 FORMAT (/,1x,'old stars%gmax    =',f10.5, '(a.u.)**(-1) ==>  new stars%gmax  '&
     263             :          &       ,'  =',f10.5,'(a.u.)**(-1) ',/,t38,'==>  new E_cut   =',&
     264             :          &            f10.5,' Ry')
     265             : 
     266             :     !--->    generate all star members
     267             :     !+gu
     268          12 :     kidx=0
     269          12 :     kidx2=0
     270             :     !-gu
     271          12 :     stars%rgphs(:,:,:) = cmplx(0.0,0.0)
     272          12 :     stars%ft2_gfx = 0.0
     273          12 :     stars%ft2_gfy = 0.0
     274       78978 :     DO  k = 1,stars%ng3
     275             : 
     276             :        CALL spgrot(&
     277             :             &               sym%nop,sym%symor,sym%mrot,sym%tau,sym%invtab,&
     278             :             &               stars%kv3(:,k),&
     279       78966 :             &               kr,phas)
     280             : 
     281       78978 :        IF (stars%kv3(3,k).EQ.0) THEN
     282             : 
     283        9897 :           DO  n = 1,sym%nop2
     284             :              !+gu
     285             :              ! -->       set up the igfft(*,3) array as (1d) fft-pointer:
     286             :              !
     287             :              !           star ------------> g-vector ------------> fft-grid & phase
     288             :              !                igfft(*,1)             igfft(*,2)           igfft(*,3)
     289             :              !
     290             :              !           size of fft-grid is chosen to be ( 3*k1d x 3*k2d x 3*k3d )
     291             :              !
     292             :              NEW=.TRUE.
     293             : 
     294       31868 :              DO n1=1,n-1
     295             :                 norm=(kr(1,n)-kr(1,n1))**2 +&
     296             :                      &                (kr(2,n)-kr(2,n1))**2 +&
     297       12643 :                      &                (kr(3,n)-kr(3,n1))**2
     298       19225 :                 IF (norm.EQ.0) NEW=.FALSE.
     299             :              ENDDO
     300             : 
     301        6582 :              IF (NEW) THEN
     302             : 
     303        5958 :                 kfx = kr(1,n)
     304        5958 :                 kfy = kr(2,n)
     305        5958 :                 kfz = kr(3,n)
     306             :                 !+guta
     307        5958 :                 gfx = cell%bmat(1,1)*kfx+cell%bmat(2,1)*kfy+cell%bmat(3,1)*kfz 
     308        5958 :                 gfy = cell%bmat(1,2)*kfx+cell%bmat(2,2)*kfy+cell%bmat(3,2)*kfz
     309             :                 !-guta
     310        5958 :                 IF (kfx.LT.0) kfx = kfx+nfftx
     311        5958 :                 IF (kfy.LT.0) kfy = kfy+nffty
     312        5958 :                 IF (kfz.LT.0) kfz = kfz+nfftz
     313        5958 :                 kfft = kfx + kfy*nfftx + kfz*nfftxy
     314             :                 !
     315             :                 ! -->            store the number of the star, its position 
     316             :                 !c                 on fft-grid and phase
     317             :                 !
     318        5958 :                 stars%igfft(kidx,1) = k
     319        5958 :                 stars%igfft(kidx,2) = kfft
     320        5958 :                 stars%pgfft(kidx)   = phas(n)
     321        5958 :                 kidx          = kidx+1
     322             :                 !
     323             :                 ! -->            now for 2d - stars
     324             :                 !
     325        5958 :                 kfft=kfx + kfy*nfftx
     326     3997019 :                 DO k2 = 1,stars%ng2
     327     3991061 :                    IF ((stars%kv3(1,k).EQ.stars%kv2(1,k2)).AND.&
     328       11916 :                         &                 (stars%kv3(2,k).EQ.stars%kv2(2,k2))) k2i = k2
     329             :                 ENDDO
     330        5958 :                 stars%igfft2(kidx2,1) = k2i
     331        5958 :                 stars%igfft2(kidx2,2) = kfft
     332        5958 :                 stars%pgfft2(kidx2)   = phas(n)
     333             :                 !+guta
     334        5958 :                 IF (xcpot%needs_grad()) THEN
     335             :                    !!                   pgft2x: exp(i*(gfx,gfy,gfz)*tau)*gfx.
     336             :                    !!                        y                             y.
     337             :                    !!                   pgft2xx: exp(i*(gfx,gfy,gfz)*tau)*gfx*gfx.
     338             :                    !!                        yy                             y   y
     339             :                    !!                        xy                             x   y
     340             : 
     341        3263 :                    stars%ft2_gfx(kidx2)  = gfx
     342        3263 :                    stars%ft2_gfy(kidx2)  = gfy
     343             :                    !                    stars%pgft2xx(kidx2) = phas(n)*gfx*gfx
     344             :                    !                    stars%pgft2yy(kidx2) = phas(n)*gfy*gfy
     345             :                    !                    stars%pgft2xy(kidx2) = phas(n)*gfx*gfy
     346             :                 ENDIF
     347             :                 !-guta
     348        5958 :                 kidx2=kidx2+1
     349             : 
     350             :              ENDIF
     351             :              !-gu
     352        6582 :              stars%ig(kr(1,n),kr(2,n),kr(3,n)) = k
     353             :              stars%rgphs(kr(1,n),kr(2,n),kr(3,n)) = &
     354        9897 :                   &         stars%rgphs(kr(1,n),kr(2,n),kr(3,n)) + phas(n)
     355             : 
     356             :           ENDDO
     357             : 
     358             :        ELSE
     359             :           !        here: kv3(3,k) =/= 0 
     360             : 
     361      249709 :           DO  n = 1,sym%nop
     362             :              !+gu
     363             :              NEW=.TRUE.
     364     1501532 :              DO n1 = 1,n-1
     365             :                 norm=(kr(1,n)-kr(1,n1))**2 +&
     366             :                      &                (kr(2,n)-kr(2,n1))**2 +&
     367      663737 :                      &                (kr(3,n)-kr(3,n1))**2
     368      837795 :                 IF (norm.EQ.0) NEW = .FALSE.
     369             :              ENDDO
     370             : 
     371      174058 :              IF (NEW) THEN
     372             : 
     373      156754 :                 kfx = kr(1,n)
     374      156754 :                 kfy = kr(2,n)
     375      156754 :                 kfz = kr(3,n)
     376      156754 :                 IF (kfx.LT.0) kfx = kfx+nfftx
     377      156754 :                 IF (kfy.LT.0) kfy = kfy+nffty
     378      156754 :                 IF (kfz.LT.0) kfz = kfz+nfftz
     379             : 
     380      156754 :                 kfft=kfx + kfy*nfftx + kfz*nfftxy
     381      156754 :                 stars%igfft(kidx,1)=k
     382      156754 :                 stars%igfft(kidx,2)=kfft
     383      156754 :                 stars%pgfft(kidx)=phas(n)
     384      156754 :                 kidx=kidx+1
     385             : 
     386             :              ENDIF
     387             :              !-gu
     388      174058 :              stars%ig(kr(1,n),kr(2,n),kr(3,n)) = k
     389             :              stars%rgphs(kr(1,n),kr(2,n),kr(3,n)) = &
     390      249709 :                   &         stars%rgphs(kr(1,n),kr(2,n),kr(3,n)) + phas(n)
     391             : 
     392             :           ENDDO
     393             : 
     394             :        ENDIF
     395             : 
     396             :     ENDDO
     397             :     !    
     398          12 :     stars%kimax=kidx-1
     399          12 :     stars%kimax2=kidx2-1
     400             :     !
     401             :     !     count number of members for each star
     402             :     !     nstr2 ... members of 2-dim stars
     403             :     !
     404             : 
     405          12 :     stars%nstr2(:) = 0
     406          12 :     stars%nstr(:) = 0
     407             : 
     408         482 :     DO k3 = -stars%mx3,stars%mx3
     409       12476 :        DO k2 = -mxx2,mxx2
     410      351230 :           DO k1 = -mxx1,mxx1
     411      338766 :              k = stars%ig(k1,k2,k3)
     412      350760 :              IF ( k .NE. 0 ) THEN
     413      162712 :                 stars%nstr(k) = stars%nstr(k) + 1
     414      162712 :                 stars%nstr2(stars%ig2(k)) = stars%nstr2(stars%ig2(k)) +(1-min0(iabs(k3),1))
     415             :              ENDIF
     416             :           ENDDO
     417             :        ENDDO
     418             :     ENDDO
     419             :     !
     420             :     ! normalize phases:
     421             :     !
     422          12 :     IF (sym%symor) THEN
     423          12 :        stars%rgphs(:,:,:) = cmplx(1.0,0.0)
     424             :     ELSE
     425           0 :        pon = 1.0 / sym%nop
     426           0 :        pon2 = 1.0 / sym%nop2
     427           0 :        DO k3 = -stars%mx3,stars%mx3
     428           0 :           DO k2 = -mxx2,mxx2
     429           0 :              DO k1 = -mxx1,mxx1
     430             : 
     431           0 :                 kfx = k1 ; kfy = k2; kfz = k3
     432           0 :                 k = stars%ig(k1,k2,k3)
     433             : 
     434           0 :                 IF (kfx.LT.0) kfx = kfx+nfftx
     435           0 :                 IF (kfy.LT.0) kfy = kfy+nffty
     436           0 :                 IF (kfz.LT.0) kfz = kfz+nfftz
     437           0 :                 kfft=kfx + kfy*nfftx + kfz*nfftxy
     438           0 :                 IF (k.GT.0) THEN
     439           0 :                    IF (stars%kv3(3,k).EQ.0) THEN
     440           0 :                       stars%rgphs(k1,k2,k3) = stars%rgphs(k1,k2,k3) * stars%nstr(k)*pon2
     441             :                    ELSE
     442           0 :                       stars%rgphs(k1,k2,k3) = stars%rgphs(k1,k2,k3) * stars%nstr(k)*pon
     443             :                    ENDIF
     444           0 :                    kidx = -90
     445           0 :                    DO i = 1, stars%kimax
     446           0 :                       IF ( stars%igfft(i,2) == kfft ) kidx = i
     447             :                    ENDDO
     448           0 :                    IF ( kidx > 0 ) stars%pgfft(kidx)=stars%rgphs(k1,k2,k3)
     449             :                 ENDIF
     450             :              ENDDO
     451             :           ENDDO
     452             :        ENDDO
     453             :     ENDIF
     454             : 
     455             :     
     456          12 :     IF (stars%mx1< mxx1 .or.  stars%mx2< mxx2) then
     457           0 :           print *,stars%mx1, mxx1, stars%mx2, mxx2       
     458           0 :           CALL judft_error("BUG in strgn")
     459             :     endif
     460             :     !stars%mx1=mxx1
     461             :     !stars%mx2=mxx2
     462             : 
     463             :     !--->    write /str0/ and /str1/ to file
     464          18 :     CALL writeStars(stars,l_xcExtended,.TRUE.)
     465             : 
     466             : 270 CONTINUE
     467             :     !
     468             :     !-->  listing
     469             :     !
     470             : 8010 FORMAT (' gmax=',f10.6,/,' nq3=  ',i5,/,' nq2=  ',i5,/)
     471             : 8020 FORMAT (' mx1= ',i5,/,' mx2= ',i5,/)
     472          18 :     WRITE (6,FMT=8030)
     473             : 8030 FORMAT (/,/,/,'   s t a r   l i s t',/)
     474             : 
     475          18 :     WRITE (6,FMT=8010) stars%gmax,stars%ng3,stars%ng2
     476          18 :     WRITE (6,'('' mx1,mx2,mx3='',3i3)') stars%mx1,stars%mx2,stars%mx3
     477          18 :     WRITE (6,'('' kimax2,kimax='',2i7,'', (start from 0)'')') stars%kimax2,&
     478          36 :          &  stars%kimax
     479             : 
     480          18 :     WRITE (6,FMT=8040)
     481             : 8040 FORMAT(/4x,'no.',5x,'kv3',9x,'sk3',9x,'sk2',5x,&
     482             :          &  'ig2',1x,'nstr',2x,'nstr2'/)
     483             : 
     484          18 :     ned1=9
     485          18 :     nint=30
     486         180 :     DO k = 1,ned1
     487         810 :        WRITE (6,FMT=8050) k,(stars%kv3(j,k),j=1,3),stars%sk3(k),&
     488         162 :             &                     stars%sk2(stars%ig2(k)),&
     489         990 :             &                     stars%ig2(k),stars%nstr(k),stars%nstr2(stars%ig2(k))
     490             :     ENDDO
     491             : 8050 FORMAT (1x,i5,3i4,2f12.6,i4,2i6)
     492             : 
     493        2941 :     DO k = ned1+1,stars%ng3,nint
     494       14615 :        WRITE (6,FMT=8050) k,(stars%kv3(j,k),j=1,3),stars%sk3(k),&
     495        2923 :             &                     stars%sk2(stars%ig2(k)),&
     496       17538 :             &                     stars%ig2(k),stars%nstr(k),stars%nstr2(stars%ig2(k))
     497        2941 :        kdone = k
     498             :     ENDDO
     499             : 
     500          18 :     IF (kdone.LT.stars%ng3) THEN
     501          90 :        WRITE (6,FMT=8050) stars%ng3,(stars%kv3(j,stars%ng3),j=1,3),stars%sk3(stars%ng3),&
     502          18 :             &                     stars%sk2(stars%ig2(stars%ng3)),&
     503         108 :             &                     stars%ig2(stars%ng3),stars%nstr(stars%ng3),stars%nstr2(stars%ig2(stars%ng3))
     504             :     ENDIF
     505             : 
     506          18 :     DEALLOCATE (gsk3,index,index3,kv3rev)
     507             : 
     508          18 :   END SUBROUTINE strgn1
     509             :   !----------------------------------------------------------------
     510          20 :   SUBROUTINE strgn2(&
     511             :        &                  stars,sym,atoms,&
     512             :        &                  vacuum,sphhar,input,cell,xcpot)
     513             :     USE m_boxdim
     514             :     USE m_sort
     515             :     USE m_spgrot
     516             :     USE m_constants
     517             :     USE m_types
     518             :     USE m_cdn_io
     519             :     IMPLICIT NONE
     520             :     TYPE(t_stars),INTENT(INOUT)  :: stars
     521             :     TYPE(t_sym),INTENT(IN)       :: sym
     522             :     TYPE(t_atoms),INTENT(IN)     :: atoms
     523             :     TYPE(t_vacuum),INTENT(IN)    :: vacuum
     524             :     TYPE(t_sphhar),INTENT(IN)    :: sphhar
     525             :     TYPE(t_input),INTENT(IN)     :: input
     526             :     TYPE(t_cell),INTENT(IN)      :: cell
     527             :     CLASS(t_xcpot),INTENT(IN)    :: xcpot
     528             :     !     ..
     529             :     !     .. Local Scalars ..
     530             :     REAL arltv1,arltv2,arltv3,s
     531             :     REAL gmi,gla,eps,gmax2,pon
     532             :     REAL gfx,gfy
     533             :     INTEGER j,k,k1,k2,k3,m0,mxx1,mxx2,mxx3,n
     534             :     INTEGER ned1,nint,kdone,i
     535             :     LOGICAL NEW,l_cdn1,l_error,l_xcExtended
     536             :     INTEGER kfx,kfy,kfz,kidx,nfftx,nffty,nfftz,kfft
     537             :     INTEGER nfftxy,norm,n1,kidx2,k2i
     538             :     !     ..
     539             :     !     .. Local Arrays ..
     540             :     REAL, ALLOCATABLE :: gsk3(:)
     541          40 :     INTEGER, ALLOCATABLE :: INDEX(:),index3(:),kv3rev(:,:)
     542             :     REAL g(3)
     543          40 :     COMPLEX phas(sym%nop)
     544          40 :     INTEGER kr(3,sym%nop),kv(3)
     545             :     INTEGER index2(stars%ng2)
     546             : 
     547             :     !     ..
     548             :     !
     549          20 :     nfftx = 3*stars%mx1
     550          20 :     nffty = 3*stars%mx2
     551          20 :     nfftz = 3*stars%mx3
     552          20 :     nfftxy= 9*stars%mx1*stars%mx2
     553             : 
     554          20 :     ALLOCATE (gsk3(stars%ng3),INDEX(stars%ng3),index3(stars%ng3),kv3rev(stars%ng3,3))
     555             : 
     556             :     !
     557          20 :     WRITE (*,*) ' stars are always ordered '
     558             : 
     559          20 :     l_xcExtended = xcpot%needs_grad()
     560             :     !--->    read in information if exists
     561          20 :     CALL readStars(stars,l_xcExtended,.FALSE.,l_error)
     562          20 :     IF(.NOT.l_error) THEN
     563             :        GOTO 270
     564             :     END IF
     565             : 
     566          15 :     mxx1 = 0
     567          15 :     mxx2 = 0
     568          15 :     mxx3 = 0
     569          15 :     stars%ng3 = 0
     570          15 :     stars%ig = 0
     571          15 :     gmax2 = stars%gmax * stars%gmax
     572             : 
     573         320 :     x_dim: DO k1 = stars%mx1,-stars%mx1,-1
     574         305 :        kv(1) = k1
     575        7095 :        y_dim: DO k2 = stars%mx2,-stars%mx2,-1
     576        6775 :           kv(2) = k2
     577      215515 :           z_dim: DO k3 = stars%mx3,-stars%mx3,-1
     578      208435 :              IF ( stars%ig(k1,k2,k3) .NE. 0 ) CYCLE z_dim
     579      120833 :              kv(3) = k3
     580             : 
     581      483332 :              DO j = 1,3
     582             :                 g(j) = kv(1)*cell%bmat(1,j) + kv(2)*cell%bmat(2,j) +&
     583      483332 :                      &                kv(3)*cell%bmat(3,j)
     584             :              ENDDO
     585      120833 :              s = g(1)**2 + g(2)**2 + g(3)**2
     586             :              !
     587             :              !--->   check if generated vector belongs to a star already
     588      127608 :              IF (s.LT.gmax2) THEN      
     589             :                 !
     590             :                 !--->    new representative found
     591             :                 CALL spgrot(&
     592             :                      &                     sym%nop,sym%symor,sym%mrot,sym%tau,sym%invtab,&
     593             :                      &                     kv,&
     594       13931 :                      &                     kr)
     595       13931 :                 stars%ng3 = stars%ng3 + 1
     596             :                 IF (stars%ng3.GT.stars%ng3) THEN
     597             :                    WRITE (6,'("nq3 = ",i5," > n3d =",i5)') stars%ng3,stars%ng3
     598             :                    CALL juDFT_error("nq3>n3d",calledby ="strgn")
     599             :                 ENDIF
     600       55724 :                 DO j = 1,3
     601       55724 :                    stars%kv3(j,stars%ng3) = kv(j)
     602             :                 ENDDO
     603       13931 :                 stars%sk3(stars%ng3) = SQRT(s)
     604      142474 :                 DO n = 1,sym%nop
     605      128543 :                    IF (mxx1.LT.kr(1,n)) mxx1 = kr(1,n)
     606      128543 :                    IF (mxx2.LT.kr(2,n)) mxx2 = kr(2,n)
     607      128543 :                    IF (mxx3.LT.kr(3,n)) mxx3 = kr(3,n)
     608      142474 :                    stars%ig(kr(1,n),kr(2,n),kr(3,n)) = stars%ng3
     609             :                 ENDDO
     610             :              ENDIF
     611             :           ENDDO z_dim
     612             :        ENDDO y_dim
     613             :     ENDDO x_dim
     614             : 
     615          15 :     IF( mxx1.GT.stars%mx1) THEN
     616           0 :        WRITE (6,'(/'' mxx1.gt.k1d. mxx1,k1d='',2i4)')  mxx1,stars%mx1
     617           0 :        CALL juDFT_error("mxx1>k1d",calledby ="strgn")
     618             :     ENDIF
     619             : 
     620          15 :     IF (mxx2.GT.stars%mx2) THEN
     621           0 :        WRITE (6,'(/'' mxx2.gt.k2d. mxx2,k2d='',2i4)')  mxx2,stars%mx2
     622           0 :        CALL juDFT_error("mxx2>k2d",calledby ="strgn")
     623             :     ENDIF
     624             : 
     625          15 :     IF (mxx3.GT.stars%mx3) THEN
     626           0 :        WRITE (6,'(/'' mxx3.gt.k3d. mxx3,k3d='',2i4)')  mxx3,stars%mx3
     627           0 :        CALL juDFT_error("mxx3>k3d",calledby ="strgn")
     628             :     ENDIF
     629             : 
     630             :     !--->    sort for increasing length sk3
     631             :     ! secondary key for equal length stars
     632       13946 :     DO  k = 1,stars%ng3
     633             :        gsk3(k) = (stars%mx1+stars%kv3(1,k)) +&
     634             :             &           (stars%mx2+stars%kv3(2,k))*(2*stars%mx1+1) +&
     635       13946 :             &           (stars%mx3+stars%kv3(3,k))*(2*stars%mx1+1)*(2*stars%mx2+1)
     636             :     ENDDO
     637          15 :     CALL sort(index(:stars%ng3),stars%sk3,gsk3)
     638       13946 :     DO k = 1,stars%ng3
     639       13931 :        kv3rev(k,1) = stars%kv3(1,INDEX(k))
     640       13931 :        kv3rev(k,2) = stars%kv3(2,INDEX(k))
     641       13931 :        kv3rev(k,3) = stars%kv3(3,INDEX(k))
     642       13946 :        gsk3(k) = stars%sk3(INDEX(k))
     643             :     ENDDO
     644       13946 :     DO k = 1,stars%ng3
     645       13931 :        stars%kv3(1,k) = kv3rev(k,1)
     646       13931 :        stars%kv3(2,k) = kv3rev(k,2)
     647       13931 :        stars%kv3(3,k) = kv3rev(k,3)
     648       13946 :        stars%sk3(k) = gsk3(k)
     649             :     ENDDO
     650             :     !
     651             :     !--->  determine true gmax and change old gmax to new gmax
     652             :     !
     653          15 :     WRITE (6,8060) stars%gmax, stars%sk3(stars%ng3)
     654          15 :     stars%gmax = stars%sk3(stars%ng3)
     655             : 8060 FORMAT (/,1x,'gmax    =',f10.5, '(a.u.)**(-1) ==>  new gmax  '&
     656             :          &       ,'  =',f10.5,'(a.u.)**(-1) ',/,t38,'==>  new E_cut   =',&
     657             :          &            f10.5,' Ry')
     658             : 
     659             :     !--->    generate all star members
     660             :     !+gu
     661          15 :     kidx=0
     662          15 :     kidx2=0
     663         458 :     DO k3 = -mxx3,mxx3
     664        9703 :        DO k2 = -mxx2,mxx2
     665      218123 :           DO k1 = -mxx1,mxx1
     666      217680 :              stars%ig(k1,k2,k3) = 0
     667             :           ENDDO
     668             :        ENDDO
     669             :     ENDDO
     670             : 
     671             :     !-gu
     672             :     !
     673             :     ! sum over phases
     674             :     !
     675          15 :     stars%rgphs(:,:,:) = cmplx(0.0,0.0)
     676          15 :     stars%igfft = 0
     677          15 :     stars%pgfft = cmplx(0.0,0.0)
     678       13946 :     starloop: DO k = 1,stars%ng3
     679             : 
     680             :        CALL spgrot(&
     681             :             &               sym%nop,sym%symor,sym%mrot,sym%tau,sym%invtab,&
     682             :             &               stars%kv3(:,k),&
     683       13931 :             &               kr,phas)
     684             :        !
     685             :        ! -->    set up the igfft(*,3) array as (1d) fft-pointer:
     686             :        !
     687             :        !        star ------------> g-vector ------------> fft-grid & phase
     688             :        !             igfft(*,1)             igfft(*,2)           igfft(*,3)
     689             :        !
     690             :        !        size of fft-grid is chosen to be ( 3*k1d x 3*k2d x 3*k3d )
     691             :        !
     692      142489 :        ops: DO n = 1,sym%nop
     693             : 
     694             :           NEW=.TRUE.
     695     1742699 :           DO n1 = 1,n-1
     696             :              norm=(kr(1,n)-kr(1,n1))**2 +&
     697             :                   &             (kr(2,n)-kr(2,n1))**2 +&
     698      807078 :                   &             (kr(3,n)-kr(3,n1))**2
     699      935621 :              IF (norm.EQ.0) NEW = .FALSE.
     700             :           ENDDO
     701             : 
     702      128543 :           IF (NEW) THEN
     703             : 
     704      101533 :              kfx = kr(1,n)
     705      101533 :              kfy = kr(2,n)
     706      101533 :              kfz = kr(3,n)
     707      101533 :              IF (kfx.LT.0) kfx = kfx+nfftx
     708      101533 :              IF (kfy.LT.0) kfy = kfy+nffty
     709      101533 :              IF (kfz.LT.0) kfz = kfz+nfftz
     710             : 
     711      101533 :              kfft=kfx + kfy*nfftx + kfz*nfftxy
     712      101533 :              stars%igfft(kidx,1)=k
     713      101533 :              stars%igfft(kidx,2)=kfft
     714      101533 :              stars%pgfft(kidx)=phas(n)
     715      101533 :              kidx=kidx+1
     716             : 
     717             :           ENDIF
     718      128543 :           stars%ig(kr(1,n),kr(2,n),kr(3,n)) = k
     719             :           stars%rgphs(kr(1,n),kr(2,n),kr(3,n)) = &
     720      142474 :                &      stars%rgphs(kr(1,n),kr(2,n),kr(3,n)) + phas(n)
     721             : 
     722             :        ENDDO ops
     723             :     ENDDO starloop
     724             :     !    
     725          15 :     stars%kimax=kidx-1
     726             :     !
     727             :     !     count number of members for each star
     728             :     !
     729       13946 :     DO k = 1,stars%ng3
     730       13946 :        stars%nstr(k) = 0
     731             :     ENDDO
     732             : 
     733         901 :     DO k3 = -mxx3,mxx3
     734        9703 :        DO k2 = -mxx2,mxx2
     735      218123 :           DO k1 = -mxx1,mxx1
     736      208435 :              k = stars%ig(k1,k2,k3)
     737      217680 :              IF ( k .NE. 0 ) THEN
     738      101533 :                 stars%nstr(k) = stars%nstr(k) + 1
     739             :              ENDIF
     740             :              !               DO k = 1,nq3
     741             :              !                  IF (ig(k1,k2,k3).eq.k) THEN
     742             :              !                     nstr(k) = nstr(k) + 1
     743             :              !                  ENDIF
     744             :              !               ENDDO
     745             :           ENDDO
     746             :        ENDDO
     747             :     ENDDO
     748             :     !
     749             :     ! normalize phases:
     750             :     !
     751          15 :     IF (sym%symor) THEN
     752          13 :        stars%rgphs(:,:,:) = cmplx(1.0,0.0)
     753             :     ELSE
     754           2 :        pon = 1.0 / sym%nop
     755             :        !$OMP PARALLEL DO &
     756             :        !$OMP DEFAULT(none) &
     757             :        !$OMP SHARED(mxx1,mxx2,mxx3,stars,nfftx,nffty,nfftz,nfftxy,pon) &
     758             :        !$OMP PRIVATE(k1,k2,k3,k,kfx,kfy,kfz,kfft,kidx,i)
     759             :        DO k3 = -mxx3,mxx3
     760        1004 :           DO k2 = -mxx2,mxx2
     761       21300 :              DO k1 = -mxx1,mxx1
     762       21250 :                 k = stars%ig(k1,k2,k3)
     763       22250 :                 IF (k.GT.0) THEN
     764        9094 :                    stars%rgphs(k1,k2,k3) = stars%rgphs(k1,k2,k3) * stars%nstr(k)*pon
     765        9094 :                    kfx = k1 ; kfy = k2; kfz = k3
     766       13321 :                    IF (kfx.LT.0) kfx = kfx+nfftx
     767       13321 :                    IF (kfy.LT.0) kfy = kfy+nffty
     768       13376 :                    IF (kfz.LT.0) kfz = kfz+nfftz
     769        9094 :                    kfft=kfx + kfy*nfftx + kfz*nfftxy
     770        9094 :                    kidx = -90
     771    47022130 :                    DO i = 1, stars%kimax
     772    47031222 :                       IF ( stars%igfft(i,2) == kfft ) kidx = i
     773             :                    ENDDO
     774        9094 :                    IF (kidx > 0 ) stars%pgfft(kidx)=stars%rgphs(k1,k2,k3)
     775             :                 ENDIF
     776             :              ENDDO
     777             :           ENDDO
     778             :        ENDDO
     779             :        !OMP END PARALLLEL DO
     780             :     ENDIF
     781          15 :     if ( stars%mx1 < mxx1 .or. stars%mx2 < mxx2 .or. stars%mx3 < mxx3 ) call &
     782           0 :          judft_error("BUG 1 in strgen") 
     783          15 :     stars%ng2 = 2 ; stars%kv2 = 0 ; stars%ig2 = 0 ; stars%kimax2= 0 ; stars%igfft2 = 0
     784          15 :     stars%sk2 = 0.0 ; stars%pgfft2 = 0.0  ; stars%nstr2 = 0
     785          15 :     IF (xcpot%needs_grad()) THEN
     786          15 :        stars%ft2_gfx = 0.0 ; stars%ft2_gfy = 0.0 
     787             :     ENDIF
     788             : 
     789             :     !--->    write /str0/ and /str1/ to file
     790          15 :     CALL timestart("writeStars") 
     791          15 :     CALL writeStars(stars,l_xcExtended,.FALSE.)
     792          20 :     CALL timestop("writeStars") 
     793             : 
     794             : 270 CONTINUE
     795             : 
     796             :     !
     797             :     !-->  listing
     798             : 8010 FORMAT (' gmax=',f10.6,/,' nq3=  ',i7,/)
     799             : 8020 FORMAT (' mx1= ',i5,/,' mx2= ',i5,' mx3= ',i5,/)
     800          20 :     WRITE (6,FMT=8030)
     801             : 8030 FORMAT (/,/,/,'   s t a r   l i s t',/)
     802             : 
     803             : 
     804          20 :     WRITE (6,FMT=8010) stars%gmax,stars%ng3
     805          20 :     WRITE (6,'('' mx1,mx2,mx3='',3i3)') stars%mx1,stars%mx2,stars%mx3
     806          20 :     WRITE (6,'('' kimax2,kimax='',2i7,'', (start from 0)'')') stars%kimax2,&
     807          40 :          &  stars%kimax
     808             : 
     809          20 :     WRITE (6,FMT=8040)
     810             : 8040 FORMAT(/6x,'no.',5x,'kv3',9x,'sk3',7x,'nstr'/)
     811             : 
     812          20 :     ned1=9
     813          20 :     nint=30
     814         200 :     DO k = 1,ned1
     815         200 :        WRITE (6,FMT=8050) k,(stars%kv3(j,k),j=1,3),stars%sk3(k),stars%nstr(k)
     816             :     ENDDO
     817             : 8050 FORMAT (1x,i7,3i4,f12.6,i6)
     818             : 
     819         665 :     DO k = ned1+1,stars%ng3,nint
     820         645 :        WRITE (6,FMT=8050) k,(stars%kv3(j,k),j=1,3),stars%sk3(k),stars%nstr(k)
     821         665 :        kdone = k
     822             :     ENDDO
     823             : 
     824          20 :     IF (kdone.LT.stars%ng3) THEN
     825         100 :        WRITE (6,FMT=8050) stars%ng3,(stars%kv3(j,stars%ng3),j=1,3),&
     826         120 :             &                     stars%sk3(stars%ng3),stars%nstr(stars%ng3)
     827             :     ENDIF
     828             : 
     829             : 
     830          20 :     DEALLOCATE (gsk3,index,index3,kv3rev)
     831             : 
     832          20 :   END SUBROUTINE strgn2
     833             : END MODULE m_strgn

Generated by: LCOV version 1.13