LCOV - code coverage report
Current view: top level - inpgen - crystal.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 122 154 79.2 %
Date: 2019-09-08 04:53:50 Functions: 2 2 100.0 %

          Line data    Source code
       1             :       MODULE m_crystal
       2             :       use m_juDFT
       3             : !********************************************************************
       4             : !      generate space group operations from lattice information
       5             : !********************************************************************
       6             :       CONTAINS
       7           3 :       SUBROUTINE crystal(
       8             :      >                   dbgfh,errfh,outfh,dispfh,dispfn,
       9             :      >                   cal_symm,cartesian,symor,oldfleur,
      10             :      >                   natin,natmax,nop48,
      11           3 :      >                   atomid,atompos,a1,a2,a3,aa,scale,i_c,
      12             :      <                   invs,zrfs,invs2,nop,nop2,
      13           3 :      <                   ngen,mmrot,ttr,ntype,nat,nops,
      14             :      <                   neq,ntyrep,zatom,natype,natrep,natmap,
      15             :      <                   mrot,tau,pos,amat,bmat,omtil)
      16             : 
      17             :       USE m_setab
      18             :       USE m_lattice, ONLY : angles
      19             :       USE m_atomsym
      20             :       USE m_spggen
      21             :       USE m_closure, ONLY : check_close
      22             :       USE m_symproperties
      23             :       USE m_generator
      24             : 
      25             :       IMPLICIT NONE
      26             : 
      27             : !===> Arguments
      28             :       LOGICAL, INTENT(IN)    :: cal_symm,cartesian,oldfleur
      29             :       INTEGER, INTENT(IN)    :: ngen,natmax,nop48,i_c
      30             :       INTEGER, INTENT(IN)    :: dbgfh,errfh,outfh,dispfh ! file handles, mainly 6
      31             :       REAL,    INTENT(IN)    :: aa
      32             :       LOGICAL, INTENT(INOUT) :: symor                    ! on input: if true, reduce symmetry if oldfleur
      33             :                                                          ! on output : if its symmorphic
      34             :       INTEGER, INTENT(INOUT) :: natin                    ! might change if atoms are  to be completed
      35             :                                                          ! by symmetry operations (if natin < 0 )
      36             :       REAL,    INTENT(INOUT) :: atomid(natmax)
      37             :       REAL,    INTENT(INOUT) :: atompos(3,natmax)
      38             :       REAL,    INTENT(IN)    :: a1(3),a2(3),a3(3)
      39             :       REAL,    INTENT(INOUT) :: scale(3)
      40             :       INTEGER, INTENT(INOUT) :: mmrot(3,3,nop48)         ! calculated here, if cal_symm is true
      41             :       REAL,    INTENT(INOUT) :: ttr(3,nop48)             ! or completed, or if only generators 
      42             :       INTEGER, INTENT (OUT)  :: ntype,nat,nops,nop,nop2
      43             :       LOGICAL, INTENT (OUT)  :: invs,zrfs,invs2
      44             :       REAL,    INTENT (OUT)  :: amat(3,3),bmat(3,3),omtil
      45             :       CHARACTER(len=4), INTENT (IN) :: dispfn
      46             : !--> actually, intent out:
      47             :       INTEGER, ALLOCATABLE :: neq(:), ntyrep(:)              ! these variables are allocated with
      48             :       REAL,    ALLOCATABLE :: zatom(:)                       ! dim 'ntype'
      49             :       INTEGER, ALLOCATABLE :: natype(:),natrep(:),natmap(:)  ! or  'nat'
      50             :       REAL,    ALLOCATABLE :: pos(:,:)                       ! or  '3,nat'
      51             :       INTEGER, ALLOCATABLE :: mrot(:,:,:)                    ! or  '3,3,nop'
      52             :       REAL,    ALLOCATABLE :: tau(:,:)                       ! or  '3,nop' here, or in atom_sym
      53             : 
      54             :                                                          ! are given (if mmrot(1,1,1) = 0 )
      55             : ! additional arguments are all variables in mod_lattice
      56             : 
      57             : !===> Parameters
      58             : !  dimensions for group operations, etc.; passed down to other routines
      59             : !                                         to force automatic storage
      60             :       INTEGER, PARAMETER  :: neig12=12
      61             : 
      62             : !===> Local Variables
      63             :       LOGICAL lerr,err_setup,invsym
      64             :       INTEGER i,j,k,n,m,na,nt,mdet,mtr,nop0,fh,inversionOp
      65             :       REAL    t,volume,eps7,eps12
      66             : 
      67           6 :       INTEGER optype(nop48)
      68             :       REAL    orth(3,3),ttau(3),rdummy(3,3),as(3,3),bs(3,3)
      69             :       REAL    amatinv(3,3),aamat(3,3),bbmat(3,3)
      70           3 :       REAL,   ALLOCATABLE :: poscc(:,:)
      71           3 :       INTEGER, ALLOCATABLE :: multtab(:,:),inv_op(:)
      72             : 
      73           3 :       eps7 = 1.0e-7 ; eps12 = 1.0e-12 
      74             : !
      75             : !---> set up the as and bs matrices needed to go between
      76             : !---> (scaled) cartesian and lattice coordinates
      77             : !
      78             :       CALL setab_scaled(
      79             :      >                  a1,a2,a3,
      80           3 :      <                  as,bs,volume)
      81             : 
      82           3 :       IF (volume < 0.00 ) THEN
      83             :         WRITE(6,'(/," Input coordinate system islefthanded;",
      84           0 :      &          " interchange a1 and a2 and try again.")')
      85           0 :          CALL juDFT_error("lefthanded system",calledby="crystal")
      86             :       ENDIF
      87             : !
      88             : !---> modify scale as necessary; note that scale(:) will
      89             : !---> be needed to convert to scaled cartesian coordinates
      90             : !
      91          12 :       DO i=1,3
      92          12 :          IF ( scale(i) .LT. 0.00 ) THEN
      93           0 :             scale(i) = sqrt( abs(scale(i)) )
      94           9 :          ELSEIF ( abs(scale(i)) .LT. 1.0e-10 ) THEN
      95           0 :             scale(i) = 1.00
      96             :          ENDIF
      97             :       ENDDO
      98             : !
      99             : !--->    generate lattice matrices (everything in mod_lattice)
     100             : !
     101             :       CALL setab(
     102             :      >           a1,a2,a3,aa,scale,
     103           3 :      <           amat,bmat,aamat,bbmat,amatinv,omtil)
     104             : 
     105             : 
     106             : !---> output: lattice
     107             : 
     108           3 :       WRITE (6,'(//," Lattice information:",/,1x,20("-"))')
     109             :       WRITE (6,'(/," overall lattice constant a0     =",
     110           3 :      &                f15.6," bohr")') aa
     111             :       WRITE (6,'(/," real-space primitive lattice vectors in units",
     112           3 :      &            " of a_{x,y,z}")')
     113           3 :       WRITE (6,'("      a_1:",3f12.6)') a1
     114           3 :       WRITE (6,'("      a_2:",3f12.6)') a2
     115           3 :       WRITE (6,'("      a_3:",3f12.6)') a3
     116             :       WRITE (6,'(/," lattice constants a_x, a_y, a_z =   ",3f12.6)')
     117           3 :      &     aa*scale(:)
     118             :       WRITE (6,'(" volume of unit cell (a.u.^3)    =",f15.6)')
     119           3 :      &      omtil
     120             : 
     121             : !---> lattice matrices
     122             : 
     123           3 :       WRITE (dbgfh,'(/,"dbg: lattice matrices")')
     124           3 :       WRITE (dbgfh,*) volume
     125           3 :       WRITE (dbgfh,'("dbg:   as      :",3(/,10x,3f12.6) )') as
     126           3 :       WRITE (dbgfh,'("dbg:   bs      :",3(/,10x,3f12.6) )') bs
     127           3 :       WRITE (dbgfh,'("dbg:   amat    :",3(/,10x,3f12.6) )') amat
     128           3 :       WRITE (dbgfh,'("dbg:   bmat    :",3(/,10x,3f12.6) )') bmat
     129           3 :       WRITE (dbgfh,'("dbg:   amatinv :",3(/,10x,3f12.6) )') amatinv
     130           3 :       WRITE (dbgfh,'("dbg:   aamat   :",3(/,10x,3f12.6) )') aamat
     131           3 :       WRITE (dbgfh,'("dbg:   bbmat   :",3(/,10x,3f12.6) )') bbmat
     132           3 :       WRITE (dbgfh,'(/,"dbg:   lattice vectors :")')
     133           3 :       CALL angles( amat )
     134           3 :       WRITE (dbgfh,'(/,"dbg:   reciprocal lattice vectors :")')
     135          12 :       rdummy = transpose( bmat )
     136           3 :       CALL angles( rdummy )
     137             : 
     138             : !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     139             : !---> atomic positions:
     140             : !--->
     141             : !---> atomic positions input can be either in scaled cartesian
     142             : !---> or lattice vector units, as determined by logical cartesian.
     143             : !---> (for supercells, sometimes more natural to input positions
     144             : !---> in scaled cartesian.)
     145             : !
     146             : !---> if natin < 0, then the representative atoms only are given;
     147             : !---> this requires that the space group symmetry be given as input.
     148             : !
     149             : !---> read in number of atoms or types
     150             : 
     151           3 :       IF ( .not.cal_symm ) THEN  ! generate atoms and 
     152             :                                  ! read in symmetry information
     153             :          CALL atom_sym(
     154             :      >                 dispfh,outfh,errfh,dispfn,natmax,
     155             :      X                 natin,atomid,atompos,
     156             :      X                 ngen,mmrot,ttr,
     157             :      >                 cartesian,i_c,symor,as,bs,nop48,
     158             :      <                 ntype,nat,nops,mrot,tau,
     159           0 :      <                 neq,ntyrep,zatom,natype,natrep,natmap,pos)
     160             : 
     161             :       ELSE
     162             : !----->  allocate arrays in mod_crystal
     163           3 :          nat = natin
     164           3 :          ALLOCATE( natype(nat),natrep(nat),natmap(nat),pos(3,nat) )
     165             : 
     166             : !--->    calculate space group symmetry
     167             :          CALL spg_gen(
     168             :      >                dispfh,outfh,errfh,dispfn,
     169             :      >                cartesian,symor,as,bs,scale,
     170             :      >                atomid,atompos,natin,nop48,neig12,
     171             :      <                ntype,nat,nops,mrot,tau,
     172           3 :      <                neq,ntyrep,zatom,natype,natrep,natmap,pos)
     173             :          ! Check whether there is an inversion center that is not at the
     174             :          ! origin and if one is found shift the crystal such that the
     175             :          ! inversion is with respect to the origin. Then recalculate
     176             :          ! symmetry operations.
     177           3 :          inversionOp = -1
     178          21 :          symOpLoop: DO k = 1, nops
     179          33 :             DO i = 1, 3
     180          93 :                DO j = 1, 3
     181          63 :                   IF (i.EQ.j) THEN
     182          30 :                      IF (mrot(i,j,k).NE.-1) CYCLE symOpLoop
     183             :                   ELSE
     184          33 :                      IF (mrot(i,j,k).NE.0) CYCLE symOpLoop
     185             :                   END IF
     186          57 :                   IF ((i.EQ.3).AND.(j.EQ.3)) THEN
     187             :                      inversionOp = k
     188             :                      EXIT symOpLoop
     189             :                   END IF
     190             :                END DO
     191             :             END DO
     192             :          END DO symOpLoop
     193           3 :          IF (inversionOp.GT.0) THEN
     194           3 :             IF(ANY(ABS(tau(:,inversionOp)).GT.eps7)) THEN
     195           0 :                WRITE(*,*) 'Found inversion center at finite position.'
     196           0 :                WRITE(*,*) 'Shifting crystal by:'
     197           0 :                WRITE(*,'(3f15.10)') -0.5*tau(:,inversionOp)            ! changed to minus
     198           0 :                WRITE(*,*) ''
     199           0 :                DO k = 1, ABS(natin)
     200           0 :                   atompos(:,k) = atompos(:,k) - 0.5*tau(:,inversionOp) ! GB`18
     201             :                END DO
     202           0 :                DEALLOCATE(neq,ntyrep,zatom,mrot,tau)
     203             :                CALL spg_gen(
     204             :      >                      dispfh,outfh,errfh,dispfn,
     205             :      >                      .FALSE.,symor,as,bs,scale,
     206             :      >                      atomid,atompos,natin,nop48,neig12,
     207             :      <                      ntype,nat,nops,mrot,tau,
     208           0 :      <                      neq,ntyrep,zatom,natype,natrep,natmap,pos)
     209             :             END IF
     210             :          END IF
     211             :       ENDIF
     212             : 
     213           3 :       WHERE ( abs( tau ) < eps7 ) tau = 0.00
     214             : 
     215             : !---> atom positions in cartesian coordinates
     216             : 
     217           3 :       ALLOCATE ( poscc(3,nat) )
     218           3 :       WHERE ( abs( pos ) < eps12 ) pos = 0.00
     219           9 :       DO n=1,nat
     220           9 :         poscc(:,n) = matmul( amat , pos(:,n) )
     221             :       ENDDO
     222           9 :       WHERE ( abs( poscc ) < eps12 ) poscc = 0.00
     223             : 
     224             : !---> check order of atoms
     225             : 
     226           3 :       lerr = .false.
     227           3 :       na = 0
     228           6 :       DO nt = 1, ntype
     229          12 :         DO n = 1, neq(nt)
     230           6 :           na = na + 1
     231           9 :           IF ( natmap(na).ne.na ) lerr = .true.
     232             :         ENDDO
     233             :       ENDDO
     234           3 :       IF ( lerr ) THEN
     235           0 :         WRITE (errfh,*)
     236           0 :         WRITE (errfh,*) '_err: crystal: ERROR. ',
     237           0 :      &  'Order of atoms is incompatible with fleur21 code.'
     238           0 :         WRITE (errfh,*) '_err: Change order of atoms in input.'
     239           0 :         err_setup = .true.
     240             : 
     241           0 :         WRITE (outfh,1030) '!===> suggested order of atoms'
     242           0 :         WRITE (outfh,1020) nat, ' ! number of atoms'
     243           0 :         WRITE (outfh,1010) '! atomid','x','y','z','type','oldidx'
     244           0 :         na = 0
     245           0 :         DO nt = 1, ntype
     246           0 :           DO n = 1, neq(nt)
     247           0 :             na = na + 1
     248           0 :             i = natmap(na)
     249           0 :             IF ( cartesian ) THEN
     250           0 :               WRITE (outfh,1000) atomid(i),matmul(as,pos(:,i)),nt,i
     251             :             ELSE
     252           0 :               WRITE (outfh,1000) atomid(i),pos(:,i),nt,i
     253             :             ENDIF
     254             :           ENDDO
     255             :         ENDDO
     256             :  1000   FORMAT (f9.2,3f19.12,' !',2i5)
     257             :  1010   FORMAT (a9,a12,a19,a19,a13,a7)
     258             :  1020   FORMAT (i10,5x,a)
     259             :  1030   FORMAT (a)
     260             :       ENDIF
     261             : 
     262             : !---> closure, multiplication table and some mapping functions
     263             : 
     264           3 :       ALLOCATE ( inv_op(nops),multtab(nops,nops) )
     265             :       CALL check_close( 
     266             :      >                 nops,mrot,tau,
     267           3 :      <                 multtab,inv_op,optype)
     268             : 
     269             : !---> determine properties of symmmetry operations, 
     270             : !---> rearrange mrot,tau as needed
     271             : 
     272             :       CALL symproperties(
     273             :      >                   nop48,optype,oldfleur,nops,multtab,amat,
     274             :      X                   symor,mrot,tau,
     275           3 :      <                   invsym,invs,zrfs,invs2,nop,nop2)
     276             : 
     277             : 
     278             : !---> redo to ensure proper mult. table and mapping functions
     279             :       CALL check_close(
     280             :      >                 nops,mrot,tau,
     281           3 :      <                 multtab,inv_op,optype)
     282             : 
     283             : !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     284             : 
     285             : !---> output: space group information
     286           3 :       WRITE (6,'(/," Space group information:",/,1x,24("-"))')
     287           3 :       WRITE (6,'(i5," operations")') nops
     288           3 :       IF ( symor ) THEN
     289           0 :         WRITE (6,'(6x,"space group is symmorphic")')
     290             :       ELSE
     291           3 :         WRITE (6,'(6x,"space group is nonsymmorphic")')
     292             :       ENDIF
     293           3 :       IF ( invs ) THEN
     294           3 :         WRITE (6,'(18x,"has inversion symmetry")')
     295           3 :         IF ( .not.invsym ) THEN
     296             :           WRITE (6,
     297           0 :      &          '(18x,"but inversion is NOT a symmetry operation")')
     298             :         ENDIF
     299             :       ELSE
     300           0 :         WRITE (6,'(18x,"does NOT have inversion symmetry")')
     301             :       ENDIF
     302           3 :       IF ( invs2 ) WRITE (6,'(18x,"has 2d inversion")')
     303           3 :       IF ( zrfs )  WRITE (6,'(18x,"has z-reflection")')
     304             : 
     305             :       WRITE (6,'(//," Operations: (in International notation)",/,
     306             :      &             " ---------------------------------------",/,3x,
     307           3 :      &    "lattice coordinates",17x,"(scaled) Cartesian coordinates")')
     308             : 
     309         147 :       DO n=1,nops
     310             : 
     311         144 :          IF ( optype(n) .lt. 0 ) THEN
     312          72 :              WRITE (6,'(16x,"_")')
     313             :          ELSE
     314          72 :              WRITE (6,*)
     315             :          ENDIF
     316             :          WRITE (6,
     317             :      &            '(" operation",i3,":  ",i1,"  (inverse =",i3,")")')
     318         144 :      &         n,abs(optype(n)),inv_op(n)
     319             : 
     320         144 :          orth = matmul( amat, matmul( mrot(:,:,n) , amatinv ) )
     321         144 :          ttau = matmul( amat, tau(:,n) )
     322         576 :          where( abs( ttau ) < 1.0e-13 ) ttau = 0.00
     323             :          WRITE (6,'("  (",3i3," )  (",f6.3," )", 7x,
     324             :      &             "  (",3f9.5," )  (",f6.3," )")')
     325        2160 :      &    ((mrot(j,i,n),i=1,3),tau(j,n),
     326        2307 :      &            (orth(j,i),i=1,3),ttau(j),j=1,3)
     327             :       ENDDO
     328             : 
     329           3 :       WRITE(outfh,'(/,"   Multiplcation table: {R_j|t_j}{R_i|t_i}")')
     330         147 :       DO j=1,nops
     331             :          WRITE(outfh,'(6x,"operation j=",i2," :",12i4,:/,(22x,12i4))')
     332         147 :      &         j,multtab(j,1:nops)
     333             :       ENDDO
     334             : 
     335             : !--->    determine a set of generators for this group
     336           3 :       CALL generator(nops,mrot,tau,outfh,errfh)
     337             : 
     338             : !---> output: the atomic positions, etc.
     339             : 
     340           3 :       WRITE (6,'(//," Atomic positions:",/,1x,17("-"))')
     341           3 :       WRITE (6,'(" atom types =",i5/,"      total =",i5)') ntype,nat
     342             :       WRITE (6,'(/,7x,"lattice coordinates",15x,
     343           3 :      &          "(scaled) Cartesian coordinates   atom")')
     344             : 
     345           3 :       na = 0
     346           6 :       DO nt=1,ntype
     347             :          WRITE (6,'(/," atom type",i4,":",2x,
     348             :      &            "atomic identification number =",f5.1,
     349             :      &            "     representative =",i4)')
     350           3 :      &             nt,zatom(nt),ntyrep(nt)
     351           9 :          DO n=1,neq(nt)
     352             :             WRITE (6,'(3f10.6,10x,3f10.6,i7)')
     353           9 :      &           pos(:,natmap(na+n)),poscc(:,natmap(na+n)),natmap(na+n)
     354             :          ENDDO
     355           6 :          na = na + neq(nt)
     356             :       ENDDO
     357             : 
     358           3 :       DEALLOCATE ( poscc,inv_op,multtab )
     359           3 :       RETURN
     360             : 
     361             :       CONTAINS   ! INTERNAL subroutines
     362             : 
     363           3 :       SUBROUTINE setab_scaled(
     364             :      >                        a1,a2,a3,
     365             :      <                        as,bs,volume)
     366             : 
     367             : !*****************************************************************
     368             : !     set up matrices needed to convert rotation matrices between
     369             : !     (scaled) cartesian and lattice coordinates
     370             : !*****************************************************************
     371             :       IMPLICIT NONE
     372             : 
     373             :       REAL, INTENT (IN)  :: a1(3),a2(3),a3(3)
     374             :       REAL, INTENT (OUT) :: as(3,3),bs(3,3),volume
     375             : 
     376           3 :       as(1,1) = a1(1)
     377           3 :       as(2,1) = a1(2)
     378           3 :       as(3,1) = a1(3)
     379           3 :       as(1,2) = a2(1)
     380           3 :       as(2,2) = a2(2)
     381           3 :       as(3,2) = a2(3)
     382           3 :       as(1,3) = a3(1)
     383           3 :       as(2,3) = a3(2)
     384           3 :       as(3,3) = a3(3)
     385             : 
     386             :       volume  = a1(1)*a2(2)*a3(3) + a2(1)*a3(2)*a1(3) +
     387             :      &          a3(1)*a1(2)*a2(3) - a1(3)*a2(2)*a3(1) -
     388           3 :      &          a2(3)*a3(2)*a1(1) - a3(3)*a1(2)*a2(1)
     389             : 
     390           3 :       bs(1,1) = (a2(2)*a3(3) - a2(3)*a3(2))/volume ! b1(1)
     391           3 :       bs(1,2) = (a2(3)*a3(1) - a2(1)*a3(3))/volume ! b1(2)
     392           3 :       bs(1,3) = (a2(1)*a3(2) - a2(2)*a3(1))/volume ! b1(3)
     393           3 :       bs(2,1) = (a3(2)*a1(3) - a3(3)*a1(2))/volume ! b2(1)
     394           3 :       bs(2,2) = (a3(3)*a1(1) - a3(1)*a1(3))/volume ! b2(2)
     395           3 :       bs(2,3) = (a3(1)*a1(2) - a3(2)*a1(1))/volume ! b2(3)
     396           3 :       bs(3,1) = (a1(2)*a2(3) - a1(3)*a2(2))/volume ! b3(1)
     397           3 :       bs(3,2) = (a1(3)*a2(1) - a1(1)*a2(3))/volume ! b3(2)
     398           3 :       bs(3,3) = (a1(1)*a2(2) - a1(2)*a2(1))/volume ! b3(3)
     399             : 
     400           3 :       RETURN
     401             :       END SUBROUTINE setab_scaled
     402             : 
     403             :       END SUBROUTINE crystal
     404             :       END MODULE m_crystal

Generated by: LCOV version 1.13