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
!---> 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
137
138
               WRITE(6,'(/," Error in atomic positions: atoms",i5,
     &             " and",i5," (at least) are the same")') n,i
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
                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