spg_gen.f 21.2 KB
Newer Older
1 2 3 4 5
      MODULE m_spggen
      use m_juDFT
!********************************************************************
!     calculates the space group operations given the lattice vectors
!     and the atomic positions.                              mw 12-99
6
!     Modified by GM (2016)
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
!********************************************************************
      CONTAINS
      SUBROUTINE spg_gen(
     >                   dispfh,outfh,errfh,dispfn,
     >                   cartesian,symor,as,bs,scale,
     >                   atomid,atompos,natin,nop48,neig12,
     <                   ntype,nat,nops,mrot,tau,
     <                   neq,ntyrep,zatom,natype,natrep,natmap,pos)

      USE m_bravaissymm
      USE m_closure, ONLY : close_pt, closure
      USE m_supercheck

      IMPLICIT NONE

      INTEGER, INTENT (IN)    :: dispfh,outfh,errfh
      INTEGER, INTENT (IN)    :: natin
      INTEGER, INTENT (IN)    :: nop48,neig12
      LOGICAL, INTENT (IN)    :: cartesian,symor
      REAL,    INTENT (IN)    :: as(3,3),bs(3,3),scale(3)
      REAL,    INTENT (INOUT) :: atompos(3,natin)
      REAL,    INTENT (IN)    :: atomid(natin)
      INTEGER, INTENT (OUT)   :: ntype,nat,nops
      CHARACTER(len=4), INTENT (IN) :: dispfn
!--> actually, intent out:
      INTEGER, ALLOCATABLE :: neq(:), ntyrep(:)    ! here, these variables are allocated with
      REAL,    ALLOCATABLE :: zatom(:)             ! the dimension 'ntype'
      INTEGER, ALLOCATABLE :: natype(:),natrep(:),natmap(:)  ! or  'nat'
      REAL,    ALLOCATABLE :: pos(:,:)                       ! or  '3,nat'
      INTEGER, ALLOCATABLE :: mrot(:,:,:)                    ! or  '3,3,nop'
      REAL,    ALLOCATABLE :: tau(:,:)                       ! or  '3,nop'


!   point group operations from bravais lattice
      INTEGER mops,mmrot(3,3,nop48),m_inv(nop48)
      INTEGER ncyl(nop48)
      INTEGER index_op(nop48),num_tr(nop48)

      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
      INTEGER ity(natin),npaint,ipaint(natin)
49
      INTEGER ios,istep0, binSize, maxBinSize
Gregor Michalicek's avatar
Gregor Michalicek committed
50
      INTEGER locBinDim(3), secondAtom(natin)
51
      INTEGER binDim(3), iBin(3)
Daniel Wortmann's avatar
Daniel Wortmann committed
52
      INTEGER trIndices(natin)
53 54 55 56 57
      CHARACTER(len=30) :: filen

      REAL    posr(3,natin),rtau(3),tr(3),disp(3,natin)
      REAL    ttau(3,nop48),trs(3,natin)
      REAL    eps7
Daniel Wortmann's avatar
Daniel Wortmann committed
58
      REAL    trVecs(3,natin)
59

60
      LOGICAL lnew,lclose, foundOne, boundary(3)
61 62
      LOGICAL l_exist

63
      INTEGER, ALLOCATABLE :: mtable(:,:), binSizes(:,:,:)
Gregor Michalicek's avatar
Gregor Michalicek committed
64
      INTEGER, ALLOCATABLE :: atomIndexBins(:,:,:,:)
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
      REAL,    ALLOCATABLE :: inipos(:,:)

      eps7= 1.0e-7 ; istep0 = 0

!---> determine the point group of the Bravais lattice
      CALL bravais_symm(
     >                  as,bs,scale,nop48,neig12,
     <                  mops,mmrot)

      ALLOCATE ( mtable(mops,mops) )
      CALL close_pt(                        ! check closure and get
     >              mops,mmrot,             ! multiplication table
     <              mtable)
      DO mop = 1, mops                      ! inverses
         DO i = 1, mops
            IF ( mtable(i,mop) == 1 ) THEN
               m_inv(mop) = i
               EXIT
            ENDIF
         ENDDO
      ENDDO

87 88 89 90 91 92 93 94 95 96
!---> read in atomic positions and shift to (-1/2,1/2] in lattice
!---> coords. also read in identification (atomic) number (atomid)
!---> to distinguish different atom types (need not be atomic number)
      DO n=1,natin
         IF (cartesian) THEN  ! convert to lattice coords. if necessary
            atompos(:,n) = matmul( bs, atompos(:,n) )
         ENDIF
         pos(:,n) = atompos(:,n) - anint( atompos(:,n) - eps7 )
      ENDDO

97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230
!---> store the positions (in lattice coord.s) given in the input file
!     ALLOCATE ( inipos(3,natin) )
!     DO n=1,natin
!        inipos(:,n) = pos(:,n)
!     enddo
!     l_inipos = .true.

!--->  check whether a displacement file exists (dispfh)
!---> these displacements should be in the same units as input file,
!---> i.e., lattice  or scaled cartesian
!
      l_exist = .false.
      INQUIRE (file=dispfn,exist=l_exist)
      IF (l_exist) then
        OPEN(dispfh,file=dispfn,status='old',form='formatted')
        DO n = 1, natin
          READ( dispfh,*,iostat=ios) tr(1:3)
          IF (ios .ne. 0) EXIT
          disp(1:3,n) = tr(1:3)
          IF (cartesian) THEN    ! convert to lattice coords. if necessary
            tr = matmul( bs, tr )
          ENDIF
          pos(:,n) = pos(:,n) + tr(:)
          pos(:,n) = pos(:,n) - anint( pos(:,n)- eps7 )
        ENDDO
        CLOSE (dispfh)
        IF ( ios==0 ) THEN
          WRITE (filen,'(a,"_",i2.2)') trim(adjustl(dispfn)),istep0+1
          OPEN (dispfh,file=filen,status='unknown',form='formatted')
          DO n = 1, natin
            WRITE (dispfh,'(3f16.8)') disp(1:3,n)
          ENDDO
          CLOSE (dispfh)
        ENDIF
      ENDIF ! l_exist

!---> sanity check: atoms must be distinct
      DO n = 1, natin
         DO i = n+1, natin
            IF ( all( abs( pos(:,n) - pos(:,i) ) < eps7 ) ) then
               WRITE(6,'(/," Error in atomic positions: atoms",i3,
     &             " and",i3," (at least) are the same")') n,i
                CALL juDFT_error("atom positions",calledby="spg_gen")
            ENDIF
         ENDDO
      ENDDO

!---> determine the number of distinct atoms based on atomic number,
!---> etc. (not necessarily symmetry inequivalent)

      ntypm = 1
      ity(1) = 1
      DO n=2, natin
         lnew = .true.
         DO i=1,n-1
            IF ( abs( atomid(i)-atomid(n) ) < eps7 ) THEN
               ity(n) = ity(i)
               lnew = .false.
               EXIT
            ENDIF
         ENDDO
         IF (lnew) then
            ntypm = ntypm + 1
            ity(n) = ntypm
         ENDIF
      ENDDO

!---> determine the order of cyclic (sub)group for each operation

      ncyl(1) = 1
      DO mop=1,mops
         ncyl(mop) = 2
         i = mop
         DO              ! multiply on left until identity
           i = mtable(mop, i)
           IF ( i == 1 ) EXIT
           ncyl(mop) = ncyl(mop)+1
         ENDDO
      ENDDO

!---> check if this is a supercell

      nat = natin
      CALL super_check(
     >                 nat,pos,ity,ntypm,
     <                 ns,trs)

!---> determine the independent atoms in the supercell and
!---> mark off the equivalent atoms

      IF ( ns > 1 ) then
         ipaint(1:nat) = 0
         DO n = 1, nat
            IF ( ipaint(n) < 0 ) CYCLE
            ipaint(n) = 0
            DO i = 2, ns        ! want to mark the equivalent atoms
               tr(1:3) = pos(1:3,n) + trs(1:3,i)
               tr(:) = tr(:) - anint( tr(:) - eps7 )
               DO j = n+1, nat
                  IF ( all( abs( tr(:) - pos(:,j) ) < eps7 ) ) THEN
                     ipaint(j) = -1
                     EXIT
                  ENDIF
               ENDDO
            ENDDO
         ENDDO
      ENDIF


!---> loop over operations (first operation is identity) and
!---> determine the possible non-primitive translations
!---> for supercells, it is possible to have more than 1 translation
!---> that also satisfies the cyclic condition; to generate subgroup
!---> consistent with the overall cell (and other conventions), we
!---> break the supercell symmetry if necessary; in this case, only
!---> a single non-primitive tranlation per rotation is allowed.
!---> (if one wants the full symmetry of the supercell, can easily
!---> modify the code here to calculate all translations.)

      npaint = 0
  100 CONTINUE

      num_tr(1) = 1
      ttau(1:3,1) = (/ 0.00 , 0.00 , 0.00 /)

      num_tr(2:mops) = 0
      DO mop = 2,  nop48            ! DO mop = 2, mops seems to be equivalent,
       IF (mop <= mops ) THEN       ! but not for ifc (7.0) with -O3

         IF ( num_tr(mop) < 0 ) cycle  ! this operation removed previously

!--->    rotate all atoms:
         DO n=1,nat
            posr(:,n) = matmul( mmrot(:,:,mop) , pos(:,n) )
Gregor Michalicek's avatar
Gregor Michalicek committed
231
            posr(:,n) = posr(:,n) - anint(posr(:,n) - eps7)
232 233
         ENDDO

234 235 236 237 238 239 240 241 242 243 244
!        Start code section A (replacing the commented part following the section)

!        Determine possible translation vectors. Each translation vector has
!        to work for all atoms.

         trVecs = 0.0
         maxTrVecs = 0
         secondAtom = 0

!!       1. Determine all possible translation vectors for the first atom

Gregor Michalicek's avatar
Gregor Michalicek committed
245
         trIndices = -1
246 247 248
         j = 1
         DO i = 1, nat
            IF (ity(i).NE.ity(j)) CYCLE
249
            tr(1:3) = pos(1:3,i) - posr(1:3,j)
250 251 252
            tr(1:3) = tr(1:3) - anint(tr(1:3) - eps7)
            maxTrVecs = maxTrVecs + 1
            trVecs(:,maxTrVecs) = tr(:)
Gregor Michalicek's avatar
Gregor Michalicek committed
253
            trIndices(maxTrVecs) = maxTrVecs
254 255 256
            secondAtom(maxTrVecs) = i
         END DO

257
!!       2. Sort atoms into "location bins"
Gregor Michalicek's avatar
Gregor Michalicek committed
258
!!          (position vectors are in the region -0.5 to 0.5)
259 260 261 262 263 264 265 266 267

         binDim(:) = CEILING(natin**(1.0/3.0)*0.5)
         !TODO: The dimension binDim should better be adjusted to the unit cell shape.
         !      This simple setting might be bad for very elongated unit cells.

         ALLOCATE(binSizes(-binDim(1)-1:binDim(1)+1,
     &                     -binDim(2)-1:binDim(2)+1,
     &                     -binDim(3)-1:binDim(3)+1))
         binSizes = 0
Gregor Michalicek's avatar
Gregor Michalicek committed
268
         DO n = 1, nat
269
            iBin(:) = NINT(binDim(:) * pos(:,n) / 0.501)
270

271
            boundary(:) = (ABS(pos(:,n))-0.5).LE.eps7
272 273 274 275 276 277 278
            DO i = -1, 1, 2
               IF((i.EQ.-1).AND.(.NOT.boundary(1))) CYCLE
               DO j = -1, 1, 2
                  IF((j.EQ.-1).AND.(.NOT.boundary(2))) CYCLE
                  DO k = -1, 1, 2
                     IF((k.EQ.-1).AND.(.NOT.boundary(3))) CYCLE
                     binSize = binSizes(i*iBin(1),j*iBin(2),k*iBin(3))
279
                     binSize = binSize + 1
280
                     binSizes(i*iBin(1),j*iBin(2),k*iBin(3)) = binSize
281 282 283
                  END DO
               END DO
            END DO
284

285 286 287 288 289 290 291 292 293 294 295 296 297
         END DO

         maxBinSize = 0
         DO i = -binDim(1)-1, binDim(1)+1
            DO j = -binDim(2)-1, binDim(2)+1
               DO k = -binDim(3)-1, binDim(3)+1
                  IF (binSizes(i,j,k).GT.maxBinSize) THEN
                     maxBinSize = binSizes(i,j,k)
                  END IF
               END DO
            END DO
         END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
298 299 300
         ALLOCATE(atomIndexBins(maxBinSize,-binDim(1)-1:binDim(1)+1,
     &                                     -binDim(2)-1:binDim(2)+1,
     &                                     -binDim(3)-1:binDim(3)+1))
301 302

         binSizes = 0
Gregor Michalicek's avatar
Gregor Michalicek committed
303
         DO n = 1, nat
304
            iBin(:) = NINT(binDim(:) * pos(:,n) / 0.501)
305

306
            boundary(:) = (ABS(pos(:,n))-0.5).LE.eps7
307 308 309 310 311 312 313
            DO i = -1, 1, 2
               IF((i.EQ.-1).AND.(.NOT.boundary(1))) CYCLE
               DO j = -1, 1, 2
                  IF((j.EQ.-1).AND.(.NOT.boundary(2))) CYCLE
                  DO k = -1, 1, 2
                     IF((k.EQ.-1).AND.(.NOT.boundary(3))) CYCLE
                     binSize = binSizes(i*iBin(1),j*iBin(2),k*iBin(3))
314
                     binSize = binSize + 1
315 316 317
                     binSizes(i*iBin(1),j*iBin(2),k*iBin(3)) = binSize
                     atomIndexBins(binSize,i*iBin(1),j*iBin(2),
     &                                     k*iBin(3)) = n
318 319 320 321 322
                  END DO
               END DO
            END DO
         END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
323 324
!!       3. Check for every other atom which of the first atom's translation
!!          vectors are applicable.
325

326
         DO j = 2, nat
Gregor Michalicek's avatar
Gregor Michalicek committed
327 328
            iTr = 1
            DO WHILE (iTr.LE.maxTrVecs)
329
               tr(1:3) = posr(1:3,j) + trVecs(1:3,trIndices(iTr))
330 331
               tr(1:3) = tr(1:3) - anint(tr(1:3) - eps7)

Gregor Michalicek's avatar
Gregor Michalicek committed
332 333 334
               iBin(:) = NINT(binDim(:) * tr(:) / 0.501)

               foundOne = .FALSE.
335 336 337

               !Search for atoms in the nearest bin
               DO k = 1,binSizes(iBin(1),iBin(2),iBin(3))
Gregor Michalicek's avatar
Gregor Michalicek committed
338
                  i = atomIndexBins(k,iBin(1),iBin(2),iBin(3))
339
                  rtau(:) = tr(:)-pos(:,i)
340 341
                  rtau(:) = rtau(:) - anint(rtau(:) - eps7)
                  IF(ALL(ABS(rtau(:)).LE.eps7)) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
342 343
                     IF (ity(i).NE.ity(j)) EXIT
                     foundOne = .TRUE.
344 345 346
                     EXIT
                  END IF
               END DO
347 348 349 350 351 352 353 354 355 356

               IF(.NOT.foundOne) THEN
               !Search for atoms in the surrounding bins
      binLoop: DO u = -1, 1
                  DO v = -1, 1
                     DO w = -1, 1
                        IF((u.EQ.0).AND.(v.EQ.0).AND.(w.EQ.0)) CYCLE
                        DO k = 1,binSizes(iBin(1)+u,iBin(2)+v,iBin(3)+w)
                           i = atomIndexBins(k,iBin(1)+u,iBin(2)+v,
     &                                         iBin(3)+w)
357
                           rtau(:) = tr(:)-pos(:,i)
358 359 360 361 362 363 364 365 366 367 368 369
                           rtau(:) = rtau(:) - anint(rtau(:) - eps7)
                           IF(ALL(ABS(rtau(:)).LE.eps7)) THEN
                              IF (ity(i).NE.ity(j)) EXIT binLoop
                              foundOne = .TRUE.
                              EXIT binLoop
                           END IF
                        END DO
                     END DO
                  END DO
               END DO binLoop
               END IF

Gregor Michalicek's avatar
Gregor Michalicek committed
370 371 372 373 374 375
               IF (foundOne) THEN
                  iTr = iTr + 1
               ELSE
                  trIndices(iTr) = trIndices(maxTrVecs)
                  maxTrVecs = maxTrVecs - 1
               END IF
376 377 378 379 380 381 382 383
            END DO
         END DO

!        Check which translation vectors are consistent with the cyclic
!        part of the group

         DO iTr = 1, maxTrVecs
            j = 1
Gregor Michalicek's avatar
Gregor Michalicek committed
384 385
            i = secondAtom(trIndices(iTr))
            tr(:) = trVecs(:,trIndices(iTr))
386 387 388 389 390 391 392 393 394 395 396 397 398 399
            rtau(1:3) = tr(1:3)

!--->       check that this is consistent with cyclic part of the group
            DO nc = 1,ncyl(mop)-1
               rtau(:) = matmul( mmrot(:,:,mop) , rtau(:) ) + tr(:)
            END DO
            rtau(1:3) = rtau(1:3) - anint( rtau(1:3) - eps7 )
            IF ( any( abs(rtau(:)) > eps7 ) ) CYCLE  ! skip if not 0

            num_tr(mop) = num_tr(mop) + 1
            ttau(1:3,mop) = tr(1:3)
            EXIT                  ! have one, go to next operation
         END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
400
         DEALLOCATE(atomIndexBins,binSizes)
401 402 403 404 405 406 407 408 409 410 411

!        End code section A (replacing the commented part following the section)

!--->    generate possible non-symmorphic pieces
 !        DO j=1,nat
 !           DO i=1,nat
 !              IF ( ity(i) .ne. ity(j) ) CYCLE
 !
 !              tr(1:3) = pos(1:3,j) - posr(1:3,i)
 !              tr(1:3) = tr(1:3) - anint( tr(1:3) - eps7 )
 !              rtau(1:3) = tr(1:3)
412 413

!--->          check that this is consistent with cyclic part of the group
414 415 416 417 418
 !              DO nc = 1,ncyl(mop)-1
 !                 rtau(:) = matmul( mmrot(:,:,mop) , rtau(:) ) + tr(:)
 !              ENDDO
 !              rtau(1:3) = rtau(1:3) - anint( rtau(1:3) - eps7 )
 !              IF ( any( abs(rtau(:)) > eps7 ) ) CYCLE  ! skip if not 0
419 420

!--->          test if this vector brings the system into registry
421 422 423 424 425 426 427 428 429
 !              IF ( l_rotmatch(tr) ) THEN   ! new translation
 !                 num_tr(mop) = num_tr(mop) + 1
 !                 ttau(1:3,mop) = tr(1:3)
 !                 EXIT                  ! have one, go to next operation
 !              ENDIF
 !
 !           ENDDO
 !           IF ( num_tr(mop) > 0 ) EXIT  ! have one, go to next operation
 !        ENDDO
430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632

         IF ( num_tr(mop) < 1 ) num_tr(m_inv(mop)) = -1 ! remove inverse also
       ENDIF
      ENDDO

      nops = 0
      DO mop = 1, mops
         IF ( num_tr(mop) < 1 ) CYCLE
         nops = nops + 1
         index_op(nops) = mop
      ENDDO


!---> check closure of group
      CALL closure(
     >             mops,mmrot,ttau,nops,index_op,
     <             lclose)

      IF ( ( ns==1 ) .AND. ( .not. lclose ) ) THEN
         WRITE (6,'(/," Congratulations, you have found a system (not",
     &    " a supercell) that breaks the algorithms. Sorry...")')
          CALL juDFT_error("Program failed :(",calledby="spg_gen")
      ENDIF

!---> supercells: if we have determined a (sub)group directly, great;
!---> however, may not have found a subgroup consistent with the
!---> supercell, so need to search harder. to do this, we will
!---> "paint" the different atoms in the primitive cell and see
!---> which choice gives the largest number of operations. because
!---> of this requirement, there is a "go to" construct. there is
!---> extra work done for the case of supercells, but seems to
!---> be necessary

      IF ( ns > 1 ) THEN

         IF ( ( .not. lclose ) .or. ( npaint > 0 ) ) THEN
            IF ( npaint == 0 ) THEN
               WRITE (6,'(/," Generating subgroups of supercell...")')
            ELSE
               IF (lclose) THEN
                  ipaint(npaint) = nops
                  WRITE (6,'(5x,"painting atom",i5,":",i3,
     &                  " operations")') npaint,nops
               ELSE
                  ipaint(npaint) = -1
                  WRITE (6,'(5x,"painting atom",i5,": not a group")')
     &                 npaint
               ENDIF
            ENDIF
            ity(1:nat) = abs( ity(1:nat) )
            npaint = -1
            DO i = 1, nat
               IF ( ipaint(i) .ne. 0 ) CYCLE
               npaint = i
               ity(i) = -ity(i)
               exit
            ENDDO
            IF ( npaint > 0 ) go to 100
         ENDIF

         IF ( npaint == -1 ) THEN  ! (re)calculate group best one
            ind = maxloc( ipaint(1:3) )  ! ind MUST be an array
            IF ( ipaint(ind(1)) < 1 ) THEN
               WRITE (6,'(/," Algorithm didnot work. Sorry...")')
               CALL juDFT_error("Supercell failure ;(",
     &                   calledby ="spg_gen")
            ENDIF
            ity(ind(1)) = -abs(ity(ind(1)))
            npaint = -2
            GO TO 100
         ENDIF

      ENDIF

      ity(1:nat) = abs( ity(1:nat) )

!---> set up various mapping functions, etc., in mod_crystal
!---> allocate arrays for space group information (mod_spgsym)
      DEALLOCATE ( mtable)          ! deallocate to avoid fragmenting heap
!      if( nopd < nops )then
!        nopd = nops
!      endif
      ALLOCATE ( mrot(3,3,nops),tau(3,nops) )

      DO n=1,nops
         mrot(:,:,n) = mmrot(:,:,index_op(n))
         tau(:,n) = ttau(:,index_op(n))
      ENDDO

      IF ( symor ) THEN
!--->   reduce symmetry to the largest symmorphic subgroup
        j = 1
        DO i = 1, nops
          IF ( all ( abs( tau(:,i) ) < eps7 ) ) THEN
            IF ( j<i ) THEN
              mrot(:,:,j) = mrot(:,:,i)
            ENDIF
            j = j + 1
          ENDIF
        ENDDO
        tau = 0.00
        IF ( nops > j-1 ) THEN
          WRITE (outfh,*) 'System has non-symmorphic symmetry with',
     &                   nops, 'operations.'
          nops = j - 1
          WRITE(outfh,*) 'Symmetry reduced to symmorphic subgroup with',
     &   nops, 'operations.'
        ENDIF

      ENDIF ! symor

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!--->    at this point, the symmetry is correct (assumed here)
!--->    now fill in arrays in mod_crystal

      natype(1:nat) = 0
      ntype = 0
      DO i =1,nat
         IF ( natype(i) .ne. 0 ) cycle
         ntype = ntype + 1   ! new atom type
         natype(i) = ntype   ! atom type
         natrep(i) = i       ! this is the representative
!--->    rotate representative and get symmetry equavalent atoms
         DO n=1,nops
            tr(:) = matmul( mrot(:,:,n) , pos(:,i) ) + tau(:,n)
            tr(:) = tr(:) - anint( tr(:) -eps7 )
!--->       this (rotated) atom done already? (correct symmetry assumed)
            DO j=i+1,nat
               IF ( natype(j) .ne. 0 ) CYCLE
               IF ( ity(j) .ne. ity(i) ) CYCLE
               IF ( any( abs( tr(:) - pos(:,j) ) > eps7 ) ) CYCLE
               natrep(j) = i      ! representative atom
               natype(j) = ntype  ! atom type
               EXIT
            ENDDO
         ENDDO
      ENDDO

!      if( ntypd < ntype )then
!        ntypd = ntype
!      endif
      ALLOCATE( neq(ntype),ntyrep(ntype),zatom(ntype) )

      neq(1:ntype) = 0
      ntyrep(1:ntype) = 0
      DO n=1,nat
         neq( natype(n) ) = neq( natype(n) ) + 1
         zatom( natype(n) ) = atomid(n)
         IF ( ntyrep( natype(n) ) == 0 ) ntyrep( natype(n) ) = n
      ENDDO

      natmap(1:nat) = 0
      j = 0
      DO nt = 1,ntype
         DO n=1,nat
            IF ( natype(n) == nt ) THEN
               j = j+ 1
               natmap(j) = n
            ENDIF
         ENDDO
      ENDDO

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

      CONTAINS ! INTERNAL SUBROUTINES

      LOGICAL FUNCTION l_rotmatch(tr)
!********************************************************************
!     determines whether the vector tr will bring the rotated
!     positions posr into regestry with the original
!********************************************************************
      IMPLICIT NONE

      REAL, INTENT (IN) :: tr(3)

      INTEGER i,j,in
      REAL    rp(3)

      l_rotmatch = .false.

      DO i=1,nat
!--->       rotated and shifted atom, reduced to (-1/2,1/2]
         rp(:) = posr(:,i) + tr(:) - anint( posr(:,i) + tr(:) - eps7 )
!--->       find which atom, if any, this matches
         in = 0
         DO j=1,nat
            IF ( ity(i).ne.ity(j) ) CYCLE
            IF ( all( abs(pos(:,j)-rp(:)) < eps7 ) ) THEN
               in = j
               EXIT
            ENDIF
         ENDDO
         IF (in == 0 ) RETURN
      ENDDO

!--->   here only if everything matches
      l_rotmatch = .true.

      END FUNCTION l_rotmatch

      END SUBROUTINE spg_gen
      END MODULE m_spggen