LCOV - code coverage report
Current view: top level - inpgen - atom_sym.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 155 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             :       MODULE m_atomsym
       2             :       use m_juDFT
       3             :       CONTAINS
       4           0 :       SUBROUTINE atom_sym(
       5             :      >                    dispfh,outfh,errfh,dispfn,natmax,
       6           0 :      X                    natin,atomid,atompos,
       7           0 :      X                    ngen,mmrot,ttr,
       8             :      >                    cartesian,i_c,symor,as,bs,nop48,
       9             :      <                    ntype,nat,nops,mrot,tau,
      10             :      <                    neq,ntyrep,zatom,natype,natrep,natmap,pos)
      11             : 
      12             : !***********************************************************************
      13             : !        reads in atomic positions and space group information
      14             : !        for the case when natin < 0 or when cal_symm = .false.
      15             : !
      16             : !--->    atomic positions:
      17             : !--->
      18             : !--->    atomic positions input can be either in scaled cartesian
      19             : !--->    or lattice vector units, as determined by logical cartesian.
      20             : !--->    (for supercells, sometimes more natural to input positions
      21             : !--->    in scaled cartesian.) only represenative atoms need be
      22             : !--->    given, but will handle all (generates all of them and then
      23             : !--->    throws out duplicates).
      24             : !--->
      25             : !--->    for each atom, give identification (atomic) number and
      26             : !--->    position, e.g.,
      27             : !                 26    0.5  0.5  0.25
      28             : !
      29             : !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      30             : !--->    space group symmetry information:
      31             : !--->
      32             : !--->    to input symmetry, give the generators (except identity)
      33             : !--->    of the group as a matrix plus a non-primitive translation,
      34             : !--->    all in lattice coordinates.
      35             : !--->
      36             : !--->    for each generator, there will be 4 lines of input
      37             : !--->    corresponding to the matrix as one would write it:
      38             : !--->
      39             : !--->        m(1,1)  m(1,2) m(1,3)    ! integer defining
      40             : !--->        m(2,1)  m(2,2) m(2,3)    ! rotation matrix for
      41             : !--->        m(3,1)  m(3,2) m(3,3)    ! column vectors
      42             : !--->        tau(1)  tau(2) tau(3)    ! non-primitive translation
      43             : !                                                               mw 12-99
      44             : !***********************************************************************
      45             : 
      46             :       USE m_closure, ONLY : closure
      47             :       IMPLICIT NONE
      48             : 
      49             : !==> Arguments
      50             : 
      51             :       INTEGER, INTENT (IN)    :: dispfh,outfh,errfh
      52             :       LOGICAL, INTENT (IN)    :: cartesian         ! how atomic positions are given
      53             :       LOGICAL, INTENT (IN)    :: symor             ! whether to reduce to symmorphic subgroup
      54             :       INTEGER, INTENT (INOUT) :: natin             ! formerly 'ntype0'
      55             :       INTEGER, INTENT (IN)    :: ngen              ! Number of generators
      56             :       INTEGER, INTENT (IN)    :: i_c               ! centering of lattice
      57             :       INTEGER, INTENT (IN)    :: nop48, natmax     ! dimensioning
      58             :       INTEGER, INTENT (INOUT) :: mmrot(3,3,nop48)
      59             :       REAL,    INTENT (INOUT) :: ttr(3,nop48)
      60             :       REAL,    INTENT (IN)    :: as(3,3),bs(3,3)
      61             :       REAL,    INTENT (INOUT) :: atomid(natmax)    ! formerly 'zat0'
      62             :       REAL,    INTENT (INOUT) :: atompos(3,natmax) ! formerly 'postype0'
      63             :       INTEGER, INTENT (OUT)   :: ntype,nat,nops
      64             :       CHARACTER(len=4), INTENT (IN) :: dispfn
      65             : !--> actually, intent out:
      66             :       INTEGER, ALLOCATABLE :: neq(:), ntyrep(:)    ! here, these variables are allocated with
      67             :       REAL,    ALLOCATABLE :: zatom(:)             ! the dimension 'ntype'
      68             :       INTEGER, ALLOCATABLE :: natype(:),natrep(:),natmap(:)  ! or  'nat'
      69             :       REAL,    ALLOCATABLE :: pos(:,:)                       ! or  '3,nat'
      70             :       INTEGER, ALLOCATABLE :: mrot(:,:,:)                    ! or  '3,3,nop'
      71             :       REAL,    ALLOCATABLE :: tau(:,:)                       ! or  '3,nop'
      72             : 
      73             : !---->Automatic arrays
      74           0 :       INTEGER nneq(natmax),icount(nop48*natmax),imap(natmax)
      75           0 :       INTEGER ity(nop48*natmax),index_op(nop48)
      76           0 :       REAL    tpos(3,nop48*natmax),atomid2(natmax)
      77           0 :       REAL    tr(3),tt(3),disp(3,natmax)
      78             : 
      79             :       INTEGER mp(3,3),mtmp(3,3)
      80             :       REAL    ttau(3),orth(3,3),tc(3,3),td(3,3)
      81             : 
      82           0 :       INTEGER mmrot2(3,3,ngen)
      83           0 :       REAL    ttr2(3,ngen)
      84             : 
      85             :       LOGICAL l_exist,lclose,l_inipos
      86             :       INTEGER n,na,ng,ncyl,nc,no,nop0,nn,nt,i,j,mops
      87             :       INTEGER ios,istep0
      88             :       CHARACTER(len=30) :: filen
      89             :       REAL,    ALLOCATABLE :: inipos(:,:)
      90             : 
      91             :       REAL, PARAMETER :: eps = 1.0e-7, isqrt3 = 1.0/sqrt(3.0), 
      92             :      &                   thrd = 1.0/3.0, mtthrd = -2.0/3.0
      93             : 
      94             :       REAL,PARAMETER :: lmat(3,3,8)=reshape((/
      95             :      &              1.0,  0.0,  0.0,      ! 1: primitive     : P
      96             :      &              0.0,  1.0,  0.0,
      97             :      &              0.0,  0.0,  1.0,
      98             :      +             -1.0,  1.0,  1.0,      ! 2: Inverse (F) 
      99             :      &              1.0, -1.0,  1.0,
     100             :      &              1.0,  1.0, -1.0,
     101             :      +              0.0,  1.0,  1.0,      ! 3: Inverse (I)
     102             :      &              1.0,  0.0,  1.0,
     103             :      &              1.0,  1.0,  0.0,
     104             :      +              1.0,  1.0,  0.0,      ! 4: Inverse (hP)
     105             :      &         -isqrt3, isqrt3, 0.0,
     106             :      &              0.0,  0.0,  1.0,
     107             :      +         0.0, isqrt3, -isqrt3,      ! 5: Inverse (hR)
     108             :      &           mtthrd, thrd, thrd,
     109             :      &             thrd, thrd, thrd,
     110             :      +              1.0,  1.0,  0.0,      ! 6: Inverse ( S (C) )
     111             :      &             -1.0,  1.0,  0.0,
     112             :      &              0.0,  0.0,  1.0,
     113             :      +              1.0,  0.0,  1.0,      ! 7: Inverse (B)
     114             :      &              0.0,  1.0,  0.0,
     115             :      &             -1.0,  0.0,  1.0, 
     116             :      +              1.0,  0.0,  0.0,      ! 8: Inverse (A)
     117             :      &              0.0,  1.0, -1.0,
     118             :      &              0.0,  1.0,  1.0/),(/3,3,8/))
     119             : 
     120           0 :        istep0 = 0
     121             : !
     122             : !---> take atomic positions and shift to (-1/2,1/2] in lattice coords.
     123             : !
     124           0 :       natin = abs(natin)
     125           0 :       DO n=1,natin
     126           0 :          IF (cartesian) THEN  ! convert to lattice coords. if necessary
     127             : !            atompos(:,n) = matmul( bs, atompos(:,n) )
     128           0 :             atompos(:,n) = matmul( lmat(:,:,i_c), atompos(:,n) )
     129             :          ENDIF
     130           0 :          atompos(:,n) = atompos(:,n) - anint( atompos(:,n) - eps )
     131             :       ENDDO
     132             : 
     133             : !--->  store the positions (in lattice coord.s) given in the input file
     134             : !      ALLOCATE ( inipos(3,natin) )
     135             : !      DO n = 1, natin
     136             : !         inipos(:,n) = atompos(:,n)
     137             : !      ENDDO
     138             : !      l_inipos = .false.
     139             : 
     140             : !---> check whether a displacement file exists (dispfh)
     141             : !---> these displacement should be in the same units as input file,
     142             : !---> i.e., lattice  or scaled cartesian
     143             : 
     144           0 :       l_exist = .false.
     145           0 :       INQUIRE (file=dispfn,exist=l_exist)
     146           0 :       IF (l_exist) THEN
     147           0 :         OPEN (dispfh,file=dispfn,status='old',form='formatted')
     148             : 
     149           0 :         DO n = 1, natin
     150           0 :           READ (dispfh,*,iostat=ios) tr(1:3)
     151           0 :           IF (ios .NE. 0) EXIT
     152           0 :           disp(1:3,n) = tr(1:3)
     153           0 :           IF (cartesian) THEN    ! convert to lattice coords. if necessary
     154           0 :             tr = matmul( bs, tr )
     155             :           ENDIF
     156           0 :           atompos(:,n) = atompos(:,n) + tr(:)
     157           0 :           atompos(:,n) = atompos(:,n) - anint( atompos(:,n)- eps )
     158             :         ENDDO
     159           0 :         CLOSE (dispfh)
     160           0 :         IF ( ios==0 ) THEN
     161           0 :           WRITE(filen,'(a,"_",i2.2)') trim(adjustl(dispfn)),istep0+1
     162           0 :           OPEN(dispfh,file=filen,status='unknown',form='formatted')
     163           0 :           DO n=1,natin
     164           0 :             WRITE (dispfh,'(3f16.8)') disp(1:3,n)
     165             :           END DO
     166           0 :           CLOSE (dispfh)
     167             :         ENDIF
     168             :       ENDIF ! l_exist
     169             : 
     170           0 :       IF ( mmrot(1,1,1)==0 ) THEN  ! &gen was used
     171             : 
     172             : !--->   save generators
     173           0 :         IF (cartesian) THEN    ! convert to lattice coords. if necessary
     174           0 :           DO ng = 2, ngen+1
     175             : !            mmrot2(:,:,1) = matmul( bs, mmrot(:,:,ng) )
     176             : !            mmrot(:,:,ng) = matmul( mmrot2(:,:,1), as )
     177           0 :             tc = mmrot(:,:,ng)
     178           0 :             td = matmul( bs, tc ) 
     179           0 :             tc = matmul( td, as ) 
     180           0 :             mmrot(:,:,ng) = NINT(tc)
     181           0 :             ttr2(:,1) = matmul( lmat(:,:,i_c), ttr(:,ng) )
     182           0 :             ttr(:,ng) = ttr2(:,1)
     183             :           ENDDO
     184             :         ENDIF
     185           0 :         mmrot2(:,:,1:ngen) = mmrot(:,:,2:ngen+1) 
     186           0 :         ttr2(:,1:ngen)     = ttr(:,2:ngen+1) 
     187             : 
     188             : !--->   set-up identity
     189           0 :         nops = 1
     190           0 :         mmrot(:,:,1) = reshape( (/ 1,0,0, 0,1,0, 0,0,1 /), (/ 3,3 /) )
     191           0 :         ttr(:,1) = 0.00 
     192             : 
     193           0 :         DO ng = 1, ngen
     194             : 
     195           0 :           n = nops + 1
     196             : 
     197             : !--->    next generator
     198           0 :           mmrot(:,:,n) = mmrot2(:,:,ng) 
     199           0 :           ttr(:,n)     = ttr2(:,ng) 
     200             : 
     201             : !--->    checking for order of cyclic (sub)group
     202             :           ncyl=1
     203           0 :           mtmp = mmrot(:,:,n)
     204           0 :           ttau = ttr(:,n)
     205             : 
     206             :           DO                    ! multiply on left until identity
     207           0 :             mp = matmul( mmrot(:,:,n) , mtmp )
     208             :             
     209             :             IF (mp(1,1)==1 .and. mp(2,2)==1 .and. mp(3,3)==1 .and.
     210             :      &          mp(1,2)==0 .and. mp(2,1)==0 .and. mp(3,2)==0 .and.
     211           0 :      &          mp(3,2)==0 .and. mp(3,1)==0 .and. mp(1,3)==0)  EXIT
     212             : 
     213             : !--->    new cyclic operation
     214           0 :             ncyl = ncyl+1
     215           0 :             mmrot(:,:,nops+ncyl) = mp
     216           0 :             ttr(:,nops+ncyl) = matmul( mmrot(:,:,n) , ttau ) + ttr(:,n)
     217             : 
     218           0 :             mtmp = mp
     219           0 :             ttau = ttr(:,nops+ncyl)
     220             :           ENDDO
     221             : 
     222             : !--->    now multiply these new operations into previous ones
     223             : !--->    (not the identity since that already included above)
     224           0 :           nop0=nops
     225           0 :           nops=nops+ncyl
     226             : 
     227           0 :           DO nc=1,ncyl
     228           0 :             DO nn=2,nop0 ! don't multiply by identity, nn=1
     229           0 :               nops = nops+1
     230           0 :               IF (nops.gt.48) then
     231           0 :                 WRITE (6,'(" Error in generators: nops>48")')
     232             :                  CALL juDFT_error("atom_sym: nops>48",
     233           0 :      &                    calledby="atom_sym")
     234             :               ENDIF
     235             :               mmrot(:,:,nops) =
     236           0 :      &                matmul( mmrot(:,:,nop0+nc) , mmrot(:,:,nn) )
     237             :               ttr(:,nops) = ttr(:,nop0+nc) +
     238           0 :      &                matmul( mmrot(:,:,nop0+nc),ttr(:,nn) )
     239             :             ENDDO
     240             :           ENDDO
     241             : 
     242             :         ENDDO ! end loop on generators
     243             : 
     244             :       ELSE ! &sym was used
     245             : 
     246           0 :         nops = ngen + 1
     247             :         
     248             :       ENDIF
     249             : 
     250             : !---> rewrite all the non-primitive translations so in (-1/2,1/2]
     251           0 :       ttr(:,1:nops) = ttr(:,1:nops) - anint( ttr(:,1:nops)- eps )
     252             : 
     253             : !--->  allocate arrays for space group information (mod_spgsym)
     254             : !      if( nopd < nops )then
     255             : !        nopd = nops
     256             : !      endif
     257           0 :       ALLOCATE ( mrot(3,3,nops), tau(3,nops) )
     258             : 
     259           0 :       DO n = 1, nops
     260           0 :          mrot(:,:,n) = mmrot(:,:,n)
     261           0 :          tau(:,n) = ttr(:,n)
     262           0 :          index_op(n) = n
     263             :       ENDDO
     264             : 
     265             : !--->    check that the group is closed, etc.
     266           0 :       WRITE(6,*) 'nops=',nops
     267             : 
     268           0 :       mops = nops   ! different values appear in the call in spg_gen
     269             :       CALL closure(
     270             :      >             mops,mrot,tau,nops,index_op,
     271           0 :      <             lclose)
     272             : 
     273           0 :       IF ( .not. lclose ) then
     274           0 :          WRITE (6,'(/," Group did not close; check input")')
     275           0 :           CALL juDFT_error("atom_sym: not closed",calledby="atom_sym")
     276             :       ENDIF
     277             : 
     278           0 :       IF (symor) THEN
     279             : !--->   reduce symmetry to the largest symmorphic subgroup
     280           0 :         j = 1
     281           0 :         DO i = 1, nops
     282           0 :           IF ( all ( abs( tau(:,i) ) < eps ) ) THEN
     283           0 :             IF ( j<i ) then
     284           0 :               mrot(:,:,j) = mrot(:,:,i)
     285             :             ENDIF
     286           0 :             j = j + 1
     287             :           ENDIF
     288             :         END DO
     289           0 :         tau = 0.00
     290           0 :         IF ( nops > j-1 ) THEN
     291           0 :           WRITE (outfh,*) 'System has non-symmorphic symmetry with',
     292           0 :      &                   nops, 'operations.'
     293           0 :           nops = j - 1
     294           0 :           WRITE (outfh,*)'Symmetry reduced to symmorphic subgroup with',
     295           0 :      &   nops, 'operations.'
     296             :         ENDIF
     297             : 
     298             :       ENDIF ! symor
     299             : 
     300             : !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     301             : !---> generate other symmetry equivalent positions from representatives
     302           0 :       ntype = 0
     303           0 :       nat = 0
     304           0 :       icount(:) = 0
     305           0 :       imap(:) = 0
     306             : 
     307           0 :       repres_atoms: DO nt = 1, natin  ! loop over input atoms
     308             : 
     309             : !--->    check whether this atom exists already
     310           0 :          DO n = 1, nat
     311             :             tt = ( atompos(:,nt) - tpos(:,n) )
     312           0 :      &           - anint( atompos(:,nt) - tpos(:,n) )
     313           0 :             IF ( all( abs(tt) < eps ) ) THEN
     314           0 :                icount(n) = icount(n) + 1
     315           0 :                imap(nt) = n
     316           0 :                IF ( abs( atomid(nt)-atomid(ity(n)) ) < eps )
     317             :      &             CYCLE repres_atoms
     318             :                CALL juDFT_error("ERROR! mismatch between atoms."
     319           0 :      +              ,calledby ="atom_sym")
     320             :             ENDIF
     321             :          ENDDO
     322             : 
     323             : !--->    add in the representative
     324           0 :          ntype = ntype+1
     325           0 :          nneq(ntype) = 1
     326           0 :          tpos(:,nat+1) = atompos(:,nt)
     327           0 :          ity(nat+1) = nt
     328           0 :          icount(nat+1) = icount(nat+1) + 1
     329           0 :          imap(nt) = -(nat+1)
     330             : 
     331             : !--->    loop over operations
     332           0 :          opts: DO no = 2, nops
     333           0 :             tr = matmul( mrot(:,:,no) , atompos(:,nt) ) + tau(:,no)
     334           0 :             tr = tr - anint( tr - eps )
     335             : !--->       check whether this is a new atom
     336           0 :             DO n = 1, nneq(ntype)
     337           0 :                tt = ( tr-tpos(:,nat+n) ) - anint( tr-tpos(:,nat+n) )
     338           0 :                IF ( all( abs(tt) < eps ) ) THEN
     339             : 
     340           0 :                   nn = ity(nat+n)
     341           0 :                   IF ( abs( atomid(nt)-atomid(nn) ) < eps ) CYCLE opts
     342             :                   WRITE (6,'(" Mismatch between atoms and",
     343           0 :      &                       " symmetry input")')
     344             :                   CALL juDFT_error("atom_sym: mismatch rotated",calledby
     345           0 :      +                 ="atom_sym")
     346             :                ENDIF
     347             :             ENDDO
     348             : !--->       new position
     349           0 :             nneq(ntype) = nneq(ntype) + 1
     350           0 :             tpos(:,nat+nneq(ntype)) = tr
     351           0 :             ity(nat+nneq(ntype)) = nt
     352             : 
     353             :          ENDDO opts
     354             : 
     355           0 :          nat = nat + nneq(ntype)
     356             :       ENDDO repres_atoms
     357             : 
     358             : !--->    allocate arrays in mod_crystal and fill as needed
     359             : 
     360             : !      if( ntypd < ntype )then
     361             : !        ntypd = ntype
     362             : !      endif
     363             : !      if ( natd < nat ) THEN
     364             : !        natd = nat
     365             : !      endif
     366           0 :       ALLOCATE ( neq(ntype), ntyrep(ntype), zatom(ntype) )
     367           0 :       ALLOCATE ( natype(nat), natrep(nat), natmap(nat), pos(3,nat) )
     368             : 
     369           0 :       na=0
     370           0 :       DO n = 1, ntype
     371           0 :          neq(n) = nneq(n)
     372           0 :          zatom(n) = atomid( ity(na+1) )
     373           0 :          ntyrep(n) = na + 1
     374           0 :          DO nn=1,neq(n)
     375           0 :             natype(na+nn) = n
     376           0 :             natrep(na+nn) = na + 1
     377           0 :             natmap(na+nn) = na + nn
     378           0 :             pos(:,na+nn)  = tpos(:,na+nn)
     379           0 :             atomid2(na+nn) = atomid(n)
     380             :          ENDDO
     381           0 :          na = na + neq(n)
     382             :       ENDDO
     383           0 :       atomid = atomid2
     384             : 
     385             : 
     386             : !---> check, that we have  
     387             : !--->    either a list of representatives
     388             : !--->    or a list of exactly all the atoms in the unit cell.
     389             : !---> stop if input is inconsistent 
     390             : 
     391           0 :       IF ( natin > ntype ) THEN
     392           0 :         IF ( .not. all( icount(1:nat)==1 ) ) THEN
     393             : 
     394             :           WRITE (outfh,'("ERROR!: atom_sym",/,5x,
     395           0 :      &       "  atom   -   type   -   rep   - times in input")')
     396             :           WRITE (outfh,'(4i10)') 
     397           0 :      &       (na, natype(na),natrep(na),icount(na),na=1,nat)
     398             : 
     399           0 :           WRITE (outfh,'("input atom --> atom, is rep, duplicate")') 
     400           0 :           DO nt = 1, natin
     401           0 :             WRITE (errfh,'(2i8,7x,l1,8x,l1)') nt,abs(imap(nt)),
     402           0 :      &         (imap(nt)<0),(icount(abs(imap(nt)))>1)
     403             :           ENDDO
     404             : 
     405             :           CALL juDFT_error
     406             :      +         ("atom_sym:too many representative atoms given in input."
     407           0 :      +         ,calledby ="atom_sym")
     408             : 
     409             :         ENDIF
     410             : !       l_inipos = .true.
     411             :       ENDIF
     412             : !     DEALLOCATE ( inipos )
     413             : 
     414           0 :       END SUBROUTINE atom_sym
     415             :       END MODULE m_atomsym

Generated by: LCOV version 1.13