LCOV - code coverage report
Current view: top level - inpgen - spg_gen.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 192 265 72.5 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_spggen
       2             :       use m_juDFT
       3             : !********************************************************************
       4             : !     calculates the space group operations given the lattice vectors
       5             : !     and the atomic positions.                              mw 12-99
       6             : !     Modified by GM (2016)
       7             : !********************************************************************
       8             :       CONTAINS
       9           3 :       SUBROUTINE spg_gen(
      10             :      >                   dispfh,outfh,errfh,dispfn,
      11             :      >                   cartesian,symor,as,bs,scale,
      12           3 :      >                   atomid,atompos,natin,nop48,neig12,
      13             :      <                   ntype,nat,nops,mrot,tau,
      14             :      <                   neq,ntyrep,zatom,natype,natrep,natmap,pos)
      15             : 
      16             :       USE m_bravaissymm
      17             :       USE m_closure, ONLY : close_pt, closure
      18             :       USE m_supercheck
      19             : 
      20             :       IMPLICIT NONE
      21             : 
      22             :       INTEGER, INTENT (IN)    :: dispfh,outfh,errfh
      23             :       INTEGER, INTENT (IN)    :: natin
      24             :       INTEGER, INTENT (IN)    :: nop48,neig12
      25             :       LOGICAL, INTENT (IN)    :: cartesian,symor
      26             :       REAL,    INTENT (IN)    :: as(3,3),bs(3,3),scale(3)
      27             :       REAL,    INTENT (INOUT) :: atompos(3,natin)
      28             :       REAL,    INTENT (IN)    :: atomid(natin)
      29             :       INTEGER, INTENT (OUT)   :: ntype,nat,nops
      30             :       CHARACTER(len=4), INTENT (IN) :: dispfn
      31             : !--> actually, intent out:
      32             :       INTEGER, ALLOCATABLE :: neq(:), ntyrep(:)    ! here, these variables are allocated with
      33             :       REAL,    ALLOCATABLE :: zatom(:)             ! the dimension 'ntype'
      34             :       INTEGER, ALLOCATABLE :: natype(:),natrep(:),natmap(:)  ! or  'nat'
      35             :       REAL,    ALLOCATABLE :: pos(:,:)                       ! or  '3,nat'
      36             :       INTEGER, ALLOCATABLE :: mrot(:,:,:)                    ! or  '3,3,nop'
      37             :       REAL,    ALLOCATABLE :: tau(:,:)                       ! or  '3,nop'
      38             : 
      39             : 
      40             : !   point group operations from bravais lattice
      41           6 :       INTEGER mops,mmrot(3,3,nop48),m_inv(nop48)
      42           6 :       INTEGER ncyl(nop48)
      43          12 :       INTEGER index_op(nop48),num_tr(nop48)
      44             : 
      45             :       INTEGER mtmp(3,3),mp(3,3)
      46             :       INTEGER u,v,w
      47             :       INTEGER i,j,k,n,mop,nc,new,ns,nt,ntypm,ind(1),iTr,maxTrVecs
      48          12 :       INTEGER ity(natin),npaint,ipaint(natin)
      49             :       INTEGER ios,istep0, binSize, maxBinSize
      50           6 :       INTEGER locBinDim(3), secondAtom(natin)
      51             :       INTEGER binDim(3), iBin(3)
      52           6 :       INTEGER trIndices(natin)
      53             :       CHARACTER(len=30) :: filen
      54             : 
      55          12 :       REAL    posr(3,natin),rtau(3),tr(3),disp(3,natin)
      56          12 :       REAL    ttau(3,nop48),trs(3,natin)
      57             :       REAL    eps7
      58           6 :       REAL    trVecs(3,natin)
      59             : 
      60             :       LOGICAL lnew,lclose, foundOne, boundary(3)
      61             :       LOGICAL l_exist
      62             : 
      63           6 :       INTEGER, ALLOCATABLE :: mtable(:,:), binSizes(:,:,:)
      64           3 :       INTEGER, ALLOCATABLE :: atomIndexBins(:,:,:,:)
      65             :       REAL,    ALLOCATABLE :: inipos(:,:)
      66             : 
      67           3 :       eps7= 1.0e-7 ; istep0 = 0
      68             : 
      69             : !---> determine the point group of the Bravais lattice
      70             :       CALL bravais_symm(
      71             :      >                  as,bs,scale,nop48,neig12,
      72           3 :      <                  mops,mmrot)
      73             : 
      74           3 :       ALLOCATE ( mtable(mops,mops) )
      75             :       CALL close_pt(                        ! check closure and get
      76             :      >              mops,mmrot,             ! multiplication table
      77           3 :      <              mtable)
      78         147 :       DO mop = 1, mops                      ! inverses
      79        3531 :          DO i = 1, mops
      80        3528 :             IF ( mtable(i,mop) == 1 ) THEN
      81         144 :                m_inv(mop) = i
      82         144 :                EXIT
      83             :             ENDIF
      84             :          ENDDO
      85             :       ENDDO
      86             : 
      87             : !---> read in atomic positions and shift to (-1/2,1/2] in lattice
      88             : !---> coords. also read in identification (atomic) number (atomid)
      89             : !---> to distinguish different atom types (need not be atomic number)
      90           9 :       DO n=1,natin
      91           6 :          IF (cartesian) THEN  ! convert to lattice coords. if necessary
      92           0 :             atompos(:,n) = matmul( bs, atompos(:,n) )
      93             :          ENDIF
      94           9 :          pos(:,n) = atompos(:,n) - anint( atompos(:,n) - eps7 )
      95             :       ENDDO
      96             : 
      97             : !---> store the positions (in lattice coord.s) given in the input file
      98             : !     ALLOCATE ( inipos(3,natin) )
      99             : !     DO n=1,natin
     100             : !        inipos(:,n) = pos(:,n)
     101             : !     enddo
     102             : !     l_inipos = .true.
     103             : 
     104             : !--->  check whether a displacement file exists (dispfh)
     105             : !---> these displacements should be in the same units as input file,
     106             : !---> i.e., lattice  or scaled cartesian
     107             : !
     108           3 :       l_exist = .false.
     109           3 :       INQUIRE (file=dispfn,exist=l_exist)
     110           3 :       IF (l_exist) then
     111           0 :         OPEN(dispfh,file=dispfn,status='old',form='formatted')
     112           0 :         DO n = 1, natin
     113           0 :           READ( dispfh,*,iostat=ios) tr(1:3)
     114           0 :           IF (ios .ne. 0) EXIT
     115           0 :           disp(1:3,n) = tr(1:3)
     116           0 :           IF (cartesian) THEN    ! convert to lattice coords. if necessary
     117           0 :             tr = matmul( bs, tr )
     118             :           ENDIF
     119           0 :           pos(:,n) = pos(:,n) + tr(:)
     120           0 :           pos(:,n) = pos(:,n) - anint( pos(:,n)- eps7 )
     121             :         ENDDO
     122           0 :         CLOSE (dispfh)
     123           0 :         IF ( ios==0 ) THEN
     124           0 :           WRITE (filen,'(a,"_",i2.2)') trim(adjustl(dispfn)),istep0+1
     125           0 :           OPEN (dispfh,file=filen,status='unknown',form='formatted')
     126           0 :           DO n = 1, natin
     127           0 :             WRITE (dispfh,'(3f16.8)') disp(1:3,n)
     128             :           ENDDO
     129           0 :           CLOSE (dispfh)
     130             :         ENDIF
     131             :       ENDIF ! l_exist
     132             : 
     133             : !---> sanity check: atoms must be distinct
     134           9 :       DO n = 1, natin
     135          12 :          DO i = n+1, natin
     136           9 :             IF ( all( abs( pos(:,n) - pos(:,i) ) < eps7 ) ) then
     137             :                WRITE(6,'(/," Error in atomic positions: atoms",i5,
     138           0 :      &             " and",i5," (at least) are the same")') n,i
     139           0 :                 CALL juDFT_error("atom positions",calledby="spg_gen")
     140             :             ENDIF
     141             :          ENDDO
     142             :       ENDDO
     143             : 
     144             : !---> determine the number of distinct atoms based on atomic number,
     145             : !---> etc. (not necessarily symmetry inequivalent)
     146             : 
     147           3 :       ntypm = 1
     148           3 :       ity(1) = 1
     149           6 :       DO n=2, natin
     150           3 :          lnew = .true.
     151           3 :          DO i=1,n-1
     152           3 :             IF ( abs( atomid(i)-atomid(n) ) < eps7 ) THEN
     153           3 :                ity(n) = ity(i)
     154           3 :                lnew = .false.
     155             :                EXIT
     156             :             ENDIF
     157             :          ENDDO
     158           3 :          IF (lnew) then
     159           0 :             ntypm = ntypm + 1
     160           0 :             ity(n) = ntypm
     161             :          ENDIF
     162             :       ENDDO
     163             : 
     164             : !---> determine the order of cyclic (sub)group for each operation
     165             : 
     166           3 :       ncyl(1) = 1
     167         147 :       DO mop=1,mops
     168         144 :          ncyl(mop) = 2
     169         144 :          i = mop
     170         195 :          DO              ! multiply on left until identity
     171         336 :            i = mtable(mop, i)
     172         336 :            IF ( i == 1 ) EXIT
     173         336 :            ncyl(mop) = ncyl(mop)+1
     174             :          ENDDO
     175             :       ENDDO
     176             : 
     177             : !---> check if this is a supercell
     178             : 
     179           3 :       nat = natin
     180             :       CALL super_check(
     181             :      >                 nat,pos,ity,ntypm,
     182           3 :      <                 ns,trs)
     183             : 
     184             : !---> determine the independent atoms in the supercell and
     185             : !---> mark off the equivalent atoms
     186             : 
     187           3 :       IF ( ns > 1 ) then
     188           0 :          ipaint(1:nat) = 0
     189           0 :          DO n = 1, nat
     190           0 :             IF ( ipaint(n) < 0 ) CYCLE
     191           0 :             ipaint(n) = 0
     192           0 :             DO i = 2, ns        ! want to mark the equivalent atoms
     193           0 :                tr(1:3) = pos(1:3,n) + trs(1:3,i)
     194           0 :                tr(:) = tr(:) - anint( tr(:) - eps7 )
     195           0 :                DO j = n+1, nat
     196           0 :                   IF ( all( abs( tr(:) - pos(:,j) ) < eps7 ) ) THEN
     197           0 :                      ipaint(j) = -1
     198           0 :                      EXIT
     199             :                   ENDIF
     200             :                ENDDO
     201             :             ENDDO
     202             :          ENDDO
     203             :       ENDIF
     204             : 
     205             : 
     206             : !---> loop over operations (first operation is identity) and
     207             : !---> determine the possible non-primitive translations
     208             : !---> for supercells, it is possible to have more than 1 translation
     209             : !---> that also satisfies the cyclic condition; to generate subgroup
     210             : !---> consistent with the overall cell (and other conventions), we
     211             : !---> break the supercell symmetry if necessary; in this case, only
     212             : !---> a single non-primitive tranlation per rotation is allowed.
     213             : !---> (if one wants the full symmetry of the supercell, can easily
     214             : !---> modify the code here to calculate all translations.)
     215             : 
     216           3 :       npaint = 0
     217             :   100 CONTINUE
     218             : 
     219           3 :       num_tr(1) = 1
     220          12 :       ttau(1:3,1) = (/ 0.00 , 0.00 , 0.00 /)
     221             : 
     222           3 :       num_tr(2:mops) = 0
     223         144 :       DO mop = 2,  nop48            ! DO mop = 2, mops seems to be equivalent,
     224         144 :        IF (mop <= mops ) THEN       ! but not for ifc (7.0) with -O3
     225             : 
     226         141 :          IF ( num_tr(mop) < 0 ) cycle  ! this operation removed previously
     227             : 
     228             : !--->    rotate all atoms:
     229         423 :          DO n=1,nat
     230         282 :             posr(:,n) = matmul( mmrot(:,:,mop) , pos(:,n) )
     231         423 :             posr(:,n) = posr(:,n) - anint(posr(:,n) - eps7)
     232             :          ENDDO
     233             : 
     234             : !        Start code section A (replacing the commented part following the section)
     235             : 
     236             : !        Determine possible translation vectors. Each translation vector has
     237             : !        to work for all atoms.
     238             : 
     239         423 :          trVecs = 0.0
     240             :          maxTrVecs = 0
     241         423 :          secondAtom = 0
     242             : 
     243             : !!       1. Determine all possible translation vectors for the first atom
     244             : 
     245         423 :          trIndices = -1
     246         141 :          j = 1
     247         423 :          DO i = 1, nat
     248         282 :             IF (ity(i).NE.ity(j)) CYCLE
     249         282 :             tr(1:3) = pos(1:3,i) - posr(1:3,j)
     250        1128 :             tr(1:3) = tr(1:3) - anint(tr(1:3) - eps7)
     251         282 :             maxTrVecs = maxTrVecs + 1
     252         282 :             trVecs(:,maxTrVecs) = tr(:)
     253         282 :             trIndices(maxTrVecs) = maxTrVecs
     254         423 :             secondAtom(maxTrVecs) = i
     255             :          END DO
     256             : 
     257             : !!       2. Sort atoms into "location bins"
     258             : !!          (position vectors are in the region -0.5 to 0.5)
     259             : 
     260         141 :          binDim(:) = CEILING(natin**(1.0/3.0)*0.5)
     261             :          !TODO: The dimension binDim should better be adjusted to the unit cell shape.
     262             :          !      This simple setting might be bad for very elongated unit cells.
     263             : 
     264             :          ALLOCATE(binSizes(-binDim(1)-1:binDim(1)+1,
     265             :      &                     -binDim(2)-1:binDim(2)+1,
     266         141 :      &                     -binDim(3)-1:binDim(3)+1))
     267         846 :          binSizes = 0
     268         423 :          DO n = 1, nat
     269         282 :             iBin(:) = NINT(binDim(:) * pos(:,n) / 0.501)
     270             : 
     271        1128 :             boundary(:) = (ABS(pos(:,n))-0.5).LE.eps7
     272         987 :             DO i = -1, 1, 2
     273         564 :                IF((i.EQ.-1).AND.(.NOT.boundary(1))) CYCLE
     274        1410 :                DO j = -1, 1, 2
     275        1128 :                   IF((j.EQ.-1).AND.(.NOT.boundary(2))) CYCLE
     276        2820 :                   DO k = -1, 1, 2
     277        2256 :                      IF((k.EQ.-1).AND.(.NOT.boundary(3))) CYCLE
     278        2256 :                      binSize = binSizes(i*iBin(1),j*iBin(2),k*iBin(3))
     279        2256 :                      binSize = binSize + 1
     280        3384 :                      binSizes(i*iBin(1),j*iBin(2),k*iBin(3)) = binSize
     281             :                   END DO
     282             :                END DO
     283             :             END DO
     284             : 
     285             :          END DO
     286             : 
     287         141 :          maxBinSize = 0
     288         846 :          DO i = -binDim(1)-1, binDim(1)+1
     289        7896 :             DO j = -binDim(2)-1, binDim(2)+1
     290       39480 :                DO k = -binDim(3)-1, binDim(3)+1
     291       21150 :                   IF (binSizes(i,j,k).GT.maxBinSize) THEN
     292         141 :                      maxBinSize = binSizes(i,j,k)
     293             :                   END IF
     294             :                END DO
     295             :             END DO
     296             :          END DO
     297             : 
     298             :          ALLOCATE(atomIndexBins(maxBinSize,-binDim(1)-1:binDim(1)+1,
     299             :      &                                     -binDim(2)-1:binDim(2)+1,
     300         141 :      &                                     -binDim(3)-1:binDim(3)+1))
     301             : 
     302         846 :          binSizes = 0
     303         423 :          DO n = 1, nat
     304         282 :             iBin(:) = NINT(binDim(:) * pos(:,n) / 0.501)
     305             : 
     306        1128 :             boundary(:) = (ABS(pos(:,n))-0.5).LE.eps7
     307         987 :             DO i = -1, 1, 2
     308         564 :                IF((i.EQ.-1).AND.(.NOT.boundary(1))) CYCLE
     309        1410 :                DO j = -1, 1, 2
     310        1128 :                   IF((j.EQ.-1).AND.(.NOT.boundary(2))) CYCLE
     311        2820 :                   DO k = -1, 1, 2
     312        2256 :                      IF((k.EQ.-1).AND.(.NOT.boundary(3))) CYCLE
     313        2256 :                      binSize = binSizes(i*iBin(1),j*iBin(2),k*iBin(3))
     314        2256 :                      binSize = binSize + 1
     315        2256 :                      binSizes(i*iBin(1),j*iBin(2),k*iBin(3)) = binSize
     316             :                      atomIndexBins(binSize,i*iBin(1),j*iBin(2),
     317        3384 :      &                                     k*iBin(3)) = n
     318             :                   END DO
     319             :                END DO
     320             :             END DO
     321             :          END DO
     322             : 
     323             : !!       3. Check for every other atom which of the first atom's translation
     324             : !!          vectors are applicable.
     325             : 
     326         423 :          DO j = 2, nat
     327             :             iTr = 1
     328         141 :             DO WHILE (iTr.LE.maxTrVecs)
     329         282 :                tr(1:3) = posr(1:3,j) + trVecs(1:3,trIndices(iTr))
     330        1128 :                tr(1:3) = tr(1:3) - anint(tr(1:3) - eps7)
     331             : 
     332        1128 :                iBin(:) = NINT(binDim(:) * tr(:) / 0.501)
     333             : 
     334         282 :                foundOne = .FALSE.
     335             : 
     336             :                !Search for atoms in the nearest bin
     337         834 :                DO k = 1,binSizes(iBin(1),iBin(2),iBin(3))
     338         693 :                   i = atomIndexBins(k,iBin(1),iBin(2),iBin(3))
     339         693 :                   rtau(:) = tr(:)-pos(:,i)
     340        2772 :                   rtau(:) = rtau(:) - anint(rtau(:) - eps7)
     341        1257 :                   IF(ALL(ABS(rtau(:)).LE.eps7)) THEN
     342         141 :                      IF (ity(i).NE.ity(j)) EXIT
     343             :                      foundOne = .TRUE.
     344             :                      EXIT
     345             :                   END IF
     346             :                END DO
     347             : 
     348             :                IF(.NOT.foundOne) THEN
     349             :                !Search for atoms in the surrounding bins
     350         987 :       binLoop: DO u = -1, 1
     351        3102 :                   DO v = -1, 1
     352        9306 :                      DO w = -1, 1
     353        3807 :                         IF((u.EQ.0).AND.(v.EQ.0).AND.(w.EQ.0)) CYCLE
     354        7191 :                         DO k = 1,binSizes(iBin(1)+u,iBin(2)+v,iBin(3)+w)
     355             :                            i = atomIndexBins(k,iBin(1)+u,iBin(2)+v,
     356        2256 :      &                                         iBin(3)+w)
     357        2256 :                            rtau(:) = tr(:)-pos(:,i)
     358        9024 :                            rtau(:) = rtau(:) - anint(rtau(:) - eps7)
     359        6063 :                            IF(ALL(ABS(rtau(:)).LE.eps7)) THEN
     360           0 :                               IF (ity(i).NE.ity(j)) EXIT binLoop
     361             :                               foundOne = .TRUE.
     362             :                               EXIT binLoop
     363             :                            END IF
     364             :                         END DO
     365             :                      END DO
     366             :                   END DO
     367             :                END DO binLoop
     368             :                END IF
     369             : 
     370         423 :                IF (foundOne) THEN
     371         141 :                   iTr = iTr + 1
     372             :                ELSE
     373         141 :                   trIndices(iTr) = trIndices(maxTrVecs)
     374         141 :                   maxTrVecs = maxTrVecs - 1
     375             :                END IF
     376             :             END DO
     377             :          END DO
     378             : 
     379             : !        Check which translation vectors are consistent with the cyclic
     380             : !        part of the group
     381             : 
     382         141 :          DO iTr = 1, maxTrVecs
     383         141 :             j = 1
     384         141 :             i = secondAtom(trIndices(iTr))
     385         141 :             tr(:) = trVecs(:,trIndices(iTr))
     386         141 :             rtau(1:3) = tr(1:3)
     387             : 
     388             : !--->       check that this is consistent with cyclic part of the group
     389         474 :             DO nc = 1,ncyl(mop)-1
     390         474 :                rtau(:) = matmul( mmrot(:,:,mop) , rtau(:) ) + tr(:)
     391             :             END DO
     392         564 :             rtau(1:3) = rtau(1:3) - anint( rtau(1:3) - eps7 )
     393         564 :             IF ( any( abs(rtau(:)) > eps7 ) ) CYCLE  ! skip if not 0
     394             : 
     395         141 :             num_tr(mop) = num_tr(mop) + 1
     396         141 :             ttau(1:3,mop) = tr(1:3)
     397           0 :             EXIT                  ! have one, go to next operation
     398             :          END DO
     399             : 
     400         141 :          DEALLOCATE(atomIndexBins,binSizes)
     401             : 
     402             : !        End code section A (replacing the commented part following the section)
     403             : 
     404             : !--->    generate possible non-symmorphic pieces
     405             :  !        DO j=1,nat
     406             :  !           DO i=1,nat
     407             :  !              IF ( ity(i) .ne. ity(j) ) CYCLE
     408             :  !
     409             :  !              tr(1:3) = pos(1:3,j) - posr(1:3,i)
     410             :  !              tr(1:3) = tr(1:3) - anint( tr(1:3) - eps7 )
     411             :  !              rtau(1:3) = tr(1:3)
     412             : 
     413             : !--->          check that this is consistent with cyclic part of the group
     414             :  !              DO nc = 1,ncyl(mop)-1
     415             :  !                 rtau(:) = matmul( mmrot(:,:,mop) , rtau(:) ) + tr(:)
     416             :  !              ENDDO
     417             :  !              rtau(1:3) = rtau(1:3) - anint( rtau(1:3) - eps7 )
     418             :  !              IF ( any( abs(rtau(:)) > eps7 ) ) CYCLE  ! skip if not 0
     419             : 
     420             : !--->          test if this vector brings the system into registry
     421             :  !              IF ( l_rotmatch(tr) ) THEN   ! new translation
     422             :  !                 num_tr(mop) = num_tr(mop) + 1
     423             :  !                 ttau(1:3,mop) = tr(1:3)
     424             :  !                 EXIT                  ! have one, go to next operation
     425             :  !              ENDIF
     426             :  !
     427             :  !           ENDDO
     428             :  !           IF ( num_tr(mop) > 0 ) EXIT  ! have one, go to next operation
     429             :  !        ENDDO
     430             : 
     431         141 :          IF ( num_tr(mop) < 1 ) num_tr(m_inv(mop)) = -1 ! remove inverse also
     432             :        ENDIF
     433             :       ENDDO
     434             : 
     435           3 :       nops = 0
     436         147 :       DO mop = 1, mops
     437         144 :          IF ( num_tr(mop) < 1 ) CYCLE
     438         144 :          nops = nops + 1
     439         147 :          index_op(nops) = mop
     440             :       ENDDO
     441             : 
     442             : 
     443             : !---> check closure of group
     444             :       CALL closure(
     445             :      >             mops,mmrot,ttau,nops,index_op,
     446           3 :      <             lclose)
     447             : 
     448           3 :       IF ( ( ns==1 ) .AND. ( .not. lclose ) ) THEN
     449             :          WRITE (6,'(/," Congratulations, you have found a system (not",
     450           0 :      &    " a supercell) that breaks the algorithms. Sorry...")')
     451           0 :           CALL juDFT_error("Program failed :(",calledby="spg_gen")
     452             :       ENDIF
     453             : 
     454             : !---> supercells: if we have determined a (sub)group directly, great;
     455             : !---> however, may not have found a subgroup consistent with the
     456             : !---> supercell, so need to search harder. to do this, we will
     457             : !---> "paint" the different atoms in the primitive cell and see
     458             : !---> which choice gives the largest number of operations. because
     459             : !---> of this requirement, there is a "go to" construct. there is
     460             : !---> extra work done for the case of supercells, but seems to
     461             : !---> be necessary
     462             : 
     463           3 :       IF ( ns > 1 ) THEN
     464             : 
     465           0 :          IF ( ( .not. lclose ) .or. ( npaint > 0 ) ) THEN
     466           0 :             IF ( npaint == 0 ) THEN
     467           0 :                WRITE (6,'(/," Generating subgroups of supercell...")')
     468             :             ELSE
     469           0 :                IF (lclose) THEN
     470           0 :                   ipaint(npaint) = nops
     471             :                   WRITE (6,'(5x,"painting atom",i5,":",i3,
     472           0 :      &                  " operations")') npaint,nops
     473             :                ELSE
     474           0 :                   ipaint(npaint) = -1
     475             :                   WRITE (6,'(5x,"painting atom",i5,": not a group")')
     476           0 :      &                 npaint
     477             :                ENDIF
     478             :             ENDIF
     479           0 :             ity(1:nat) = abs( ity(1:nat) )
     480           0 :             npaint = -1
     481           0 :             DO i = 1, nat
     482           0 :                IF ( ipaint(i) .ne. 0 ) CYCLE
     483           0 :                npaint = i
     484           0 :                ity(i) = -ity(i)
     485           0 :                exit
     486             :             ENDDO
     487           0 :             IF ( npaint > 0 ) go to 100
     488             :          ENDIF
     489             : 
     490           0 :          IF ( npaint == -1 ) THEN  ! (re)calculate group best one
     491           0 :             ind = maxloc( ipaint(1:3) )  ! ind MUST be an array
     492           0 :             IF ( ipaint(ind(1)) < 1 ) THEN
     493           0 :                WRITE (6,'(/," Algorithm didnot work. Sorry...")')
     494             :                CALL juDFT_error("Supercell failure ;(",
     495           0 :      &                   calledby ="spg_gen")
     496             :             ENDIF
     497           0 :             ity(ind(1)) = -abs(ity(ind(1)))
     498           0 :             npaint = -2
     499           0 :             GO TO 100
     500             :          ENDIF
     501             : 
     502             :       ENDIF
     503             : 
     504           3 :       ity(1:nat) = abs( ity(1:nat) )
     505             : 
     506             : !---> set up various mapping functions, etc., in mod_crystal
     507             : !---> allocate arrays for space group information (mod_spgsym)
     508           3 :       DEALLOCATE ( mtable)          ! deallocate to avoid fragmenting heap
     509             : !      if( nopd < nops )then
     510             : !        nopd = nops
     511             : !      endif
     512           3 :       ALLOCATE ( mrot(3,3,nops),tau(3,nops) )
     513             : 
     514         147 :       DO n=1,nops
     515         144 :          mrot(:,:,n) = mmrot(:,:,index_op(n))
     516         147 :          tau(:,n) = ttau(:,index_op(n))
     517             :       ENDDO
     518             : 
     519           3 :       IF ( symor ) THEN
     520             : !--->   reduce symmetry to the largest symmorphic subgroup
     521           0 :         j = 1
     522           0 :         DO i = 1, nops
     523           0 :           IF ( all ( abs( tau(:,i) ) < eps7 ) ) THEN
     524           0 :             IF ( j<i ) THEN
     525           0 :               mrot(:,:,j) = mrot(:,:,i)
     526             :             ENDIF
     527           0 :             j = j + 1
     528             :           ENDIF
     529             :         ENDDO
     530           0 :         tau = 0.00
     531           0 :         IF ( nops > j-1 ) THEN
     532           0 :           WRITE (outfh,*) 'System has non-symmorphic symmetry with',
     533           0 :      &                   nops, 'operations.'
     534           0 :           nops = j - 1
     535           0 :           WRITE(outfh,*) 'Symmetry reduced to symmorphic subgroup with',
     536           0 :      &   nops, 'operations.'
     537             :         ENDIF
     538             : 
     539             :       ENDIF ! symor
     540             : 
     541             : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     542             : !--->    at this point, the symmetry is correct (assumed here)
     543             : !--->    now fill in arrays in mod_crystal
     544             : 
     545           3 :       natype(1:nat) = 0
     546           3 :       ntype = 0
     547           9 :       DO i =1,nat
     548           6 :          IF ( natype(i) .ne. 0 ) cycle
     549           3 :          ntype = ntype + 1   ! new atom type
     550           3 :          natype(i) = ntype   ! atom type
     551           3 :          natrep(i) = i       ! this is the representative
     552             : !--->    rotate representative and get symmetry equavalent atoms
     553         150 :          DO n=1,nops
     554         144 :             tr(:) = matmul( mrot(:,:,n) , pos(:,i) ) + tau(:,n)
     555         576 :             tr(:) = tr(:) - anint( tr(:) -eps7 )
     556             : !--->       this (rotated) atom done already? (correct symmetry assumed)
     557         291 :             DO j=i+1,nat
     558         144 :                IF ( natype(j) .ne. 0 ) CYCLE
     559           9 :                IF ( ity(j) .ne. ity(i) ) CYCLE
     560           9 :                IF ( any( abs( tr(:) - pos(:,j) ) > eps7 ) ) CYCLE
     561           3 :                natrep(j) = i      ! representative atom
     562           3 :                natype(j) = ntype  ! atom type
     563         150 :                EXIT
     564             :             ENDDO
     565             :          ENDDO
     566             :       ENDDO
     567             : 
     568             : !      if( ntypd < ntype )then
     569             : !        ntypd = ntype
     570             : !      endif
     571           3 :       ALLOCATE( neq(ntype),ntyrep(ntype),zatom(ntype) )
     572             : 
     573           6 :       neq(1:ntype) = 0
     574           6 :       ntyrep(1:ntype) = 0
     575           9 :       DO n=1,nat
     576           6 :          neq( natype(n) ) = neq( natype(n) ) + 1
     577           6 :          zatom( natype(n) ) = atomid(n)
     578           9 :          IF ( ntyrep( natype(n) ) == 0 ) ntyrep( natype(n) ) = n
     579             :       ENDDO
     580             : 
     581           3 :       natmap(1:nat) = 0
     582             :       j = 0
     583          12 :       DO nt = 1,ntype
     584          12 :          DO n=1,nat
     585           9 :             IF ( natype(n) == nt ) THEN
     586           6 :                j = j+ 1
     587           6 :                natmap(j) = n
     588             :             ENDIF
     589             :          ENDDO
     590             :       ENDDO
     591             : 
     592             : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     593             : 
     594             :       CONTAINS ! INTERNAL SUBROUTINES
     595             : 
     596             :       LOGICAL FUNCTION l_rotmatch(tr)
     597             : !********************************************************************
     598             : !     determines whether the vector tr will bring the rotated
     599             : !     positions posr into regestry with the original
     600             : !********************************************************************
     601             :       IMPLICIT NONE
     602             : 
     603             :       REAL, INTENT (IN) :: tr(3)
     604             : 
     605             :       INTEGER i,j,in
     606             :       REAL    rp(3)
     607             : 
     608             :       l_rotmatch = .false.
     609             : 
     610             :       DO i=1,nat
     611             : !--->       rotated and shifted atom, reduced to (-1/2,1/2]
     612             :          rp(:) = posr(:,i) + tr(:) - anint( posr(:,i) + tr(:) - eps7 )
     613             : !--->       find which atom, if any, this matches
     614             :          in = 0
     615             :          DO j=1,nat
     616             :             IF ( ity(i).ne.ity(j) ) CYCLE
     617             :             IF ( all( abs(pos(:,j)-rp(:)) < eps7 ) ) THEN
     618             :                in = j
     619             :                EXIT
     620             :             ENDIF
     621             :          ENDDO
     622             :          IF (in == 0 ) RETURN
     623             :       ENDDO
     624             : 
     625             : !--->   here only if everything matches
     626             :       l_rotmatch = .true.
     627             : 
     628             :       END FUNCTION l_rotmatch
     629             : 
     630             :       END SUBROUTINE spg_gen
     631             :       END MODULE m_spggen
     632             : 

Generated by: LCOV version 1.13