LCOV - code coverage report
Current view: top level - inpgen - lattice2.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 17 273 6.2 %
Date: 2019-09-08 04:53:50 Functions: 1 3 33.3 %

          Line data    Source code
       1             :       MODULE m_lattice
       2             :       use m_juDFT
       3             : !---------------------------------------------------------------------!
       4             : ! eventually easy input of all 14 Bravais lattices in agreement with  !
       5             : ! 'International Tables of Crystallography'                           !
       6             : ! table 2.1.1, ( page13 in 3rd Ed. )                                  !
       7             : !---------------------------------------------------------------------!
       8             :       CONTAINS
       9           0 :       SUBROUTINE lattice2( 
      10           0 :      >                    buffer,xl_buffer,errfh,bfh,nline,
      11             :      <                    a1,a2,a3,aa,scale,mat,i_c,ios )
      12             : 
      13             :       USE m_constants
      14             :       IMPLICIT NONE
      15             : 
      16             : !==> Arguments
      17             :       INTEGER, INTENT (IN) :: errfh,bfh,nline
      18             :       INTEGER, INTENT (IN)                 :: xl_buffer
      19             :       CHARACTER(len=xl_buffer), INTENT(IN) :: buffer
      20             :       REAL,    INTENT (OUT) :: a1(3),a2(3),a3(3)
      21             :       REAL,    INTENT (OUT) :: aa                ! overall scaling constant
      22             :       REAL,    INTENT (OUT) :: scale(3),mat(3,3) ! for trigonal lattices
      23             :       INTEGER, INTENT (OUT) :: i_c,ios
      24             : 
      25             : !==> Local Variables
      26             :       CHARACTER(len=40) :: latsys
      27             :       REAL    :: a0,a_rho
      28             :       REAL    :: a,b,c
      29             :       REAL    :: alpha,beta,gamma
      30             :       REAL    :: c1(3),c2(3),c3(3)
      31             :       REAL    :: ar,br,cr,b1,b2,am(3,3)
      32             :       REAL    :: ca,cb,at
      33             :       INTEGER :: i,j,err,i1,i2
      34             :       LOGICAL :: noangles
      35             : 
      36             :       REAL, PARAMETER :: eps = 1.0e-7,
      37             :      &              sqrt2  =  1.4142135623730950,
      38             :      &              sqrt3  =  1.7320508075688773,
      39             :      &              sqrt32 =  0.86602540378444,
      40             :      &             msqrt32 = -0.86602540378444
      41             : 
      42             :       REAL :: lmat(3,3,8)
      43             :       DATA  lmat /  1.0,  0.0,  0.0,      ! 1: primitive     : P
      44             :      &              0.0,  1.0,  0.0,
      45             :      &              0.0,  0.0,  1.0,  
      46             :      +              0.0,  0.5,  0.5,      ! 2: face-centered : F
      47             :      &              0.5,  0.0,  0.5,
      48             :      &              0.5,  0.5,  0.0,
      49             :      +             -0.5,  0.5,  0.5,      ! 3: body-centered : I
      50             :      &              0.5, -0.5,  0.5,
      51             :      &              0.5,  0.5, -0.5,
      52             :      +            0.5, msqrt32, 0.0,      ! 4: hexagonal-P   : hP, hcp
      53             :      &            0.5,  sqrt32, 0.0,
      54             :      &            0.0,     0.0, 1.0,   
      55             :      +              0.0, -1.0,  1.0,      ! 5: hexagonal-R   : hR, trigonal
      56             :      &           sqrt32,  0.5,  1.0,
      57             :      &          msqrt32,  0.5,  1.0, 
      58             :      +              0.5, -0.5,  0.0,      ! 6: base-centered: S (C)
      59             :      &              0.5,  0.5,  0.0,
      60             :      &              0.0,  0.0,  1.0,
      61             :      +              0.5,  0.0, -0.5,      ! 7: base-centered: B
      62             :      &              0.0,  1.0,  0.0,
      63             :      &              0.5,  0.0,  0.5,
      64             :      +              1.0,  0.0,  0.0,      ! 8: base-centered: A
      65             :      &              0.0,  0.5,  0.5,
      66             :      &              0.0, -0.5,  0.5/
      67             : 
      68             : !===> namelists
      69             :       NAMELIST /lattice/ latsys,a0,a,b,c,alpha,beta,gamma
      70             : 
      71           0 :       noangles = .false.
      72           0 :       latsys = ' ' ; a0 = 0.0   
      73           0 :       a = 0.0      ; b = 0.0    ; c = 0.0 
      74           0 :       alpha = 0.0  ; beta = 0.0 ; gamma = 0.0   
      75           0 :       scale = 0.0  ; mat = 0.0
      76             :  
      77           0 :       READ (bfh,lattice,err=911,end=911,iostat=ios)
      78             :       
      79           0 :       IF ( abs(a0) < eps ) a0 = 1.0
      80           0 :       IF ( abs(a)  < eps ) a  = 1.0
      81           0 :       IF ( abs(b)  < eps ) b  = a
      82           0 :       IF ( abs(c)  < eps ) c  = a
      83           0 :       IF ( abs(alpha)  < eps ) alpha  = 90.0
      84           0 :       IF ( abs(beta)   < eps ) beta   = 90.0
      85           0 :       IF ( abs(gamma)  < eps ) gamma  = 90.0
      86             : 
      87           0 :       IF ( alpha > pi_const ) THEN     ! deg
      88           0 :         ar = alpha * pi_const / 180.0 
      89           0 :         br = beta  * pi_const / 180.0 
      90           0 :         cr = gamma * pi_const / 180.0 
      91             :       ELSE                       ! radians
      92             :         ar = alpha
      93             :         br = beta
      94             :         cr = gamma
      95             :       ENDIF
      96             : 
      97           0 :       scale(1) = a
      98           0 :       scale(2) = b
      99           0 :       scale(3) = c
     100             : 
     101           0 :       latsys = ADJUSTL(latsys)
     102             : 
     103             : !===>  1: cubic-P          (cP) sc
     104             : 
     105           0 :       IF ( latsys =='cubic-P'.OR.latsys =='cP'.OR.latsys =='sc'.OR.
     106             :      &     latsys =='simple-cubic' ) THEN
     107             : 
     108           0 :         noangles=.true.
     109           0 :         i_c = 1
     110           0 :         a1 = lmat(:,1,i_c) * scale(:)
     111           0 :         a2 = lmat(:,2,i_c) * scale(:)
     112           0 :         a3 = lmat(:,3,i_c) * scale(:)
     113           0 :         scale = 1.0
     114             : 
     115             :         IF ( a.NE.b .OR. a.NE.c ) err = 11
     116             :         IF ( ar.NE.br .OR. ar.NE.cr .OR. ar.NE.(pi_const/2.0) ) err = 12
     117             : 
     118             : !===>  2: cubic-F          (cF) fcc
     119             : 
     120           0 :       ELSEIF ( latsys =='cubic-F'.OR.latsys =='cF'.OR.latsys =='fcc'.OR.
     121             :      &         latsys =='face-centered-cubic' ) THEN
     122             : 
     123           0 :         noangles=.true.
     124           0 :         i_c = 2
     125           0 :         a1 = lmat(:,1,i_c) * scale(:)
     126           0 :         a2 = lmat(:,2,i_c) * scale(:)
     127           0 :         a3 = lmat(:,3,i_c) * scale(:)
     128           0 :         scale = 1.0
     129             : 
     130             :         IF ( a.NE.b .OR. a.NE.c ) err = 21
     131             : 
     132             : !===>  3: cubic-I          (cI) bcc
     133             : 
     134           0 :       ELSEIF ( latsys =='cubic-I'.OR.latsys =='cI'.OR.latsys =='bcc'.OR.
     135             :      &         latsys =='body-centered-cubic' ) THEN
     136             : 
     137           0 :         noangles=.true.
     138           0 :         i_c = 3
     139           0 :         a1 = lmat(:,1,i_c) * scale(:)
     140           0 :         a2 = lmat(:,2,i_c) * scale(:)
     141           0 :         a3 = lmat(:,3,i_c) * scale(:)
     142           0 :         scale = 1.0
     143             : 
     144             :         IF ( a.NE.b .OR. a.NE.c ) err = 31
     145             : 
     146             : !===>  4: hexagonal-P      (hP) hcp
     147             : 
     148             :       ELSEIF ( latsys =='hexagonal-P'.OR.latsys =='hP'.OR.latsys =='hcp'
     149           0 :      &                               .OR.latsys =='hexagonal' ) THEN
     150             : 
     151           0 :         noangles=.true.
     152           0 :         i_c = 4
     153           0 :         a1 = lmat(:,1,i_c) * scale(:)
     154           0 :         a2 = lmat(:,2,i_c) * scale(:)
     155           0 :         a3 = lmat(:,3,i_c) * scale(:)
     156           0 :         scale = 1.0
     157             : 
     158             :         IF ( a.NE.b ) err = 41
     159             : 
     160             : !===>  4.1 : hexagonal-P   60 degrees variant
     161             : 
     162           0 :       ELSEIF ( latsys =='hdp' ) THEN
     163             : 
     164           0 :         noangles=.true.
     165           0 :         i_c = 4
     166           0 :         a1 =  lmat((/2,1,3/),1,i_c) * scale(:)
     167           0 :         a2 = -lmat((/2,1,3/),2,i_c) * scale(:)
     168           0 :         a3 = lmat(:,3,i_c) * scale(:)
     169           0 :         scale = 1.0
     170           0 :         i_c = 9
     171             : 
     172             :         IF ( a.NE.b ) err = 41
     173             : 
     174             : !===>  5.1: hexagonal-R      (hR) trigonal,( rhombohedral )
     175             : 
     176             :       ELSEIF ( latsys =='hexagonal-R'.OR.latsys =='hR'.OR.latsys =='r'.
     177             :      &      OR.latsys =='R'.OR.latsys =='rhombohedral'.OR.
     178           0 :      &         latsys =='rho'.OR.latsys =='trigonal' ) THEN
     179             : 
     180           0 :         noangles=.false.
     181           0 :         i_c = 5
     182           0 :         a1 = lmat(:,1,i_c)
     183           0 :         a2 = lmat(:,2,i_c)
     184           0 :         a3 = lmat(:,3,i_c)
     185             :         
     186           0 :         IF ( a.NE.b ) err = 51
     187           0 :         IF ( b.NE.c ) err = 51
     188             :         IF ( alpha.EQ.0.0  .OR. 
     189             :      &       alpha.NE.beta .OR. alpha.NE.gamma ) err = 52
     190             : 
     191           0 :         at = sqrt( 2.0 / 3.0 * ( 1 - cos(ar) ) )
     192           0 :         am(:,1) = a1 ; am(:,2) = a2 ; am(:,3) = a3
     193           0 :         am(1,:) = at * am(1,:)
     194           0 :         am(2,:) = at * am(2,:)
     195           0 :         am(3,:) = cos(asin(at)) * am(3,:)
     196           0 :         a1 = am(:,1)*a ; a2 = am(:,2)*b ; a3 = am(:,3)*c
     197           0 :         scale = 1.0
     198             : 
     199           0 :         CALL angles( am )
     200             :  
     201             : !===>  5.2: hexagonal-R      (hR) trigonal,( rhombohedral )
     202             : 
     203             :       ELSEIF (latsys =='hexagonal-R2'.OR.latsys =='hR2'.OR.latsys =='r2'
     204           0 :      &    .OR.latsys =='R2'.OR.latsys =='rhombohedral2'.OR.
     205             :      &        latsys =='trigonal2' ) THEN
     206             : 
     207           0 :         noangles=.false.
     208           0 :         i_c = 5
     209           0 :         a1 = lmat(:,1,i_c)
     210           0 :         a2 = lmat(:,2,i_c)
     211           0 :         a3 = lmat(:,3,i_c)
     212             : 
     213           0 :         IF ( a.NE.b ) err = 53
     214             :         IF ( alpha.NE.90 ) err = 54
     215             :         IF ( alpha.NE.beta .OR. gamma.NE.120 ) err = 54
     216             : 
     217           0 :         mat(:,1) = a*lmat(:,1,4) ! to transfer atom
     218           0 :         mat(:,2) = a*lmat(:,2,4) ! positions hex -> trig
     219           0 :         mat(:,3) = c*lmat(:,3,4) ! in struct_inp.f
     220             : 
     221             : ! transform hex -> rho
     222           0 :         a_rho = sqrt( 3*a*a + c*c)/3.0
     223           0 :         ar = acos( 1. - 9./(6.+2*(c/a)**2) )
     224           0 :         at = sqrt( 2.0 / 3.0 * ( 1 - cos(ar) ) )
     225             : 
     226           0 :         a = a_rho * at
     227           0 :         b = a_rho * at
     228           0 :         c = a_rho * cos( asin(at) )
     229             : 
     230           0 :         am(:,1) = a1 ; am(:,2) = a2 ; am(:,3) = a3
     231           0 :         am(1,:) = a*am(1,:)
     232           0 :         am(2,:) = b*am(2,:)
     233           0 :         am(3,:) = c*am(3,:)
     234           0 :         a1 = am(:,1) ; a2 = am(:,2) ; a3 = am(:,3)
     235           0 :         scale = 1.0
     236             : 
     237           0 :         CALL angles( am )
     238             : 
     239             : !===>  6: tetragonal-P     (tP) st
     240             : 
     241             :       ELSEIF ( latsys =='tetragonal-P'.OR.latsys =='st'.OR.latsys =='tP'
     242           0 :      &     .OR.latsys =='simple-tetragonal' ) THEN
     243             : 
     244           0 :         noangles=.true.
     245           0 :         i_c = 1
     246           0 :         a1 = lmat(:,1,i_c) * scale(:)
     247           0 :         a2 = lmat(:,2,i_c) * scale(:)
     248           0 :         a3 = lmat(:,3,i_c) * scale(:)
     249           0 :         scale = 1.0
     250             :           
     251             :         IF ( a.NE.b ) err = 61
     252             :         IF ( ar.NE.br .OR. ar.NE.cr .OR. ar.NE.(pi_const/2.0)  ) err= 62
     253             : 
     254             : !===>  7: tetragonal-I     (tI) bct
     255             : 
     256             :       ELSEIF (latsys =='tetragonal-I'.OR.latsys =='tI'.OR.latsys =='bct'
     257           0 :      &    .OR.latsys =='body-centered-tetragonal' ) THEN
     258             : 
     259           0 :         noangles=.true.
     260           0 :         i_c = 3
     261           0 :         a1 = lmat(:,1,i_c) * scale(:)
     262           0 :         a2 = lmat(:,2,i_c) * scale(:)
     263           0 :         a3 = lmat(:,3,i_c) * scale(:)
     264           0 :         scale = 1.0
     265             :           
     266             :         IF ( a.NE.b ) err = 61
     267             : 
     268             : !===>  8: orthorhombic-P   (oP) 
     269             : 
     270           0 :       ELSEIF ( latsys =='orthorhombic-P'.OR.latsys =='oP'.OR.
     271             :      &         latsys =='simple-orthorhombic' ) THEN
     272             : 
     273           0 :         noangles=.true.
     274           0 :         i_c = 1
     275           0 :         a1 = lmat(:,1,i_c) * scale(:)
     276           0 :         a2 = lmat(:,2,i_c) * scale(:)
     277           0 :         a3 = lmat(:,3,i_c) * scale(:)
     278           0 :         scale = 1.0
     279             :           
     280             : !===>  9: orthorhombic-F   (oF) 
     281             : 
     282             :       ELSEIF ( latsys =='orthorhombic-F'.OR.latsys =='oF'.OR.
     283           0 :      &         latsys =='orF'.OR.
     284             :      &         latsys =='face-centered-orthorhombic' ) THEN
     285             : 
     286           0 :         noangles=.true.
     287           0 :         i_c = 2
     288           0 :         a1 = lmat(:,1,i_c) * scale(:)
     289           0 :         a2 = lmat(:,2,i_c) * scale(:)
     290           0 :         a3 = lmat(:,3,i_c) * scale(:)
     291           0 :         scale = 1.0
     292             :           
     293             : !===> 10: orthorhombic-I   (oI) 
     294             : 
     295             :       ELSEIF ( latsys =='orthorhombic-I'.OR.latsys =='oI'.OR.
     296           0 :      &         latsys =='orI'.OR.
     297             :      &         latsys =='body-centered-orthorhombic' ) THEN
     298             : 
     299           0 :         noangles=.true.
     300           0 :         i_c = 3
     301           0 :         a1 = lmat(:,1,i_c) * scale(:)
     302           0 :         a2 = lmat(:,2,i_c) * scale(:)
     303           0 :         a3 = lmat(:,3,i_c) * scale(:)
     304           0 :         scale = 1.0
     305             :           
     306             : !===> 11: orthorhombic-S   (oS) (oC) 
     307             : 
     308             :       ELSEIF ( latsys =='orthorhombic-S'.OR.latsys =='orthorhombic-C'.or
     309           0 :      &        .latsys =='oS'.OR.latsys =='oC'.OR.latsys =='orC'.OR.
     310             :      &         latsys =='base-centered-orthorhombic' ) THEN
     311             : 
     312           0 :         noangles=.true.
     313           0 :         i_c = 6
     314           0 :         a1 = lmat(:,1,i_c) * scale(:)
     315           0 :         a2 = lmat(:,2,i_c) * scale(:)
     316           0 :         a3 = lmat(:,3,i_c) * scale(:)
     317           0 :         scale = 1.0
     318             :           
     319             : !===> 11a: orthorhombic-A   (oA)
     320             : 
     321             :       ELSEIF ( latsys =='orthorhombic-A'.OR.latsys =='oA'.OR
     322           0 :      &        .latsys =='orA'.OR.
     323             :      &         latsys =='base-centered-orthorhombic2' ) THEN
     324             : 
     325           0 :         noangles=.true.
     326           0 :         i_c = 8
     327           0 :         a1 = lmat(:,1,i_c) * scale(:)
     328           0 :         a2 = lmat(:,2,i_c) * scale(:)
     329           0 :         a3 = lmat(:,3,i_c) * scale(:)
     330           0 :         scale = 1.0
     331             : 
     332             : !===> 11b: orthorhombic-B   (oB)
     333             : 
     334             :       ELSEIF ( latsys =='orthorhombic-B'.OR.latsys =='oB'.OR
     335           0 :      &        .latsys =='orB'.OR.
     336             :      &         latsys =='base-centered-orthorhombic3' ) THEN
     337             : 
     338           0 :         noangles=.true.
     339           0 :         i_c = 7
     340           0 :         a1 = lmat(:,1,i_c) * scale(:)
     341           0 :         a2 = lmat(:,2,i_c) * scale(:)
     342           0 :         a3 = lmat(:,3,i_c) * scale(:)
     343           0 :         scale = 1.0
     344             : 
     345             : !===> 12: monoclinic-P     (mP) 
     346             :       ELSEIF ( latsys =='monoclinic-P'.OR.latsys =='mP'.OR
     347           0 :      &        .latsys =='moP'.OR.
     348             :      &         latsys =='simple-monoclinic' ) THEN
     349             : 
     350           0 :         noangles=.false.
     351           0 :         i_c = 1
     352             : 
     353           0 :         IF ( (abs(alpha-90.0)<eps).AND.(abs(beta-90.0)<eps) ) THEN
     354           0 :            IF ( ABS(gamma - 90.0) <eps )  CALL juDFT_error
     355           0 :      +          ("no monoclinic angle!",calledby ="lattice2")
     356             :         ELSE 
     357             :            CALL juDFT_error
     358             :      +          ("lattice2: Please take gamma as monoclinic angle!"
     359           0 :      +          ,calledby ="lattice2")
     360             :         ENDIF  
     361           0 :         CALL brvmat ( alpha, beta, gamma, am )
     362           0 :         a1 = a*am(:,1)
     363           0 :         a2 = b*am(:,2)
     364           0 :         a3 = c*am(:,3) 
     365           0 :         scale = 1.0
     366             : 
     367           0 :         CALL angles( am )
     368             : 
     369             : !===> 13: monoclinic-C (mC) 
     370             :       ELSEIF ( latsys =='monoclinic-C'.OR.latsys =='mC'.OR
     371           0 :      &        .latsys =='moC'.OR.
     372             :      &         latsys =='centered-monoclinic' ) THEN
     373             : 
     374             :          CALL juDFT_error
     375             :      +    ("Monoclinic setup",hint
     376             :      +        ="Please use gamma as monoclinic angle"//
     377           0 :      +        " and center on A or B !",calledby ="lattice2")
     378             : 
     379             : !===> 13a monoclinic-A (mA)
     380             :       ELSEIF ( latsys =='monoclinic-A'.OR.latsys =='mA'.OR
     381           0 :      &        .latsys =='moA'.OR.
     382             :      &         latsys =='centered-monoclinic2' ) THEN
     383             : 
     384           0 :         noangles=.false.
     385           0 :         IF ( (abs(alpha-90.0)<eps).AND.(abs(beta-90.0)<eps) ) THEN
     386           0 :            IF ( abs(gamma - 90.0) <eps )  CALL juDFT_error
     387           0 :      +          ("no monoclinic angle!",calledby ="lattice2")
     388             :         ELSE
     389             :            CALL juDFT_error("Please take gamma as monoclinic angle!"
     390           0 :      +          ,calledby ="lattice2")
     391             :         ENDIF
     392           0 :         CALL brvmat ( alpha, beta, gamma, am )
     393           0 :         am(:,1) = a * am(:,1)
     394           0 :         am(:,2) = b * am(:,2)
     395           0 :         am(:,3) = c * am(:,3)
     396           0 :         i_c = 8
     397           0 :         am = matmul ( am, lmat(:,:,i_c) )
     398           0 :         a1 = am(:,1)
     399           0 :         a2 = am(:,2)
     400           0 :         a3 = am(:,3)
     401           0 :         CALL angles( am )
     402           0 :         scale = 1.0 
     403             : 
     404             : !===> 13b monoclinic-B (mB)
     405             :       ELSEIF ( latsys =='monoclinic-B'.OR.latsys =='mB'.OR
     406           0 :      &        .latsys =='moB'.OR.
     407             :      &         latsys =='centered-monoclinic3' ) THEN
     408             : 
     409           0 :         noangles=.false.
     410           0 :         IF ( (abs(alpha-90.0)<eps).AND.(abs(beta-90.0)<eps) ) THEN
     411           0 :            IF ( ABS(gamma - 90.0) <eps )  CALL juDFT_error
     412           0 :      +          ("no monoclinic angle!",calledby ="lattice2")
     413             :         ELSE
     414             :            CALL juDFT_error("Please take gamma as monoclinic angle!"
     415           0 :      +          ,calledby ="lattice2")
     416             :         ENDIF
     417           0 :         CALL brvmat ( alpha, beta, gamma, am )
     418           0 :         am(:,1) = a * am(:,1)
     419           0 :         am(:,2) = b * am(:,2)
     420           0 :         am(:,3) = c * am(:,3)
     421           0 :         i_c = 7
     422           0 :         am = matmul ( am, lmat(:,:,i_c) )
     423           0 :         a1 = am(:,1)
     424           0 :         a2 = am(:,2)
     425           0 :         a3 = am(:,3)
     426           0 :         CALL angles( am )
     427           0 :         scale = 1.0 
     428             : 
     429             : !===> 13c monoclinic-I (mI)
     430             :       ELSEIF ( latsys =='monoclinic-I'.OR.latsys =='mI'.OR
     431           0 :      &        .latsys =='moI' ) THEN
     432             : 
     433           0 :         noangles=.false.
     434           0 :         IF ( (abs(alpha-90.0)<eps).AND.(abs(beta-90.0)<eps) ) THEN
     435           0 :            IF ( ABS(gamma - 90.0) <eps )  CALL juDFT_error
     436           0 :      +          ("no monoclinic angle!",calledby ="lattice2")
     437             :         ELSE
     438             :            CALL juDFT_error("Please take gamma as monoclinic angle!"
     439           0 :      +          ,calledby ="lattice2")
     440             :         ENDIF
     441           0 :         CALL brvmat ( alpha, beta, gamma, am )
     442           0 :         am(:,1) = a * am(:,1)
     443           0 :         am(:,2) = b * am(:,2)
     444           0 :         am(:,3) = c * am(:,3)
     445           0 :         i_c = 3
     446           0 :         am = matmul ( am, lmat(:,:,i_c) )
     447           0 :         a1 = am(:,1)
     448           0 :         a2 = am(:,2)
     449           0 :         a3 = am(:,3)
     450           0 :         CALL angles( am )
     451           0 :         scale = 1.0 
     452             : 
     453             : !===> 14: triclinic        (aP) 
     454             : 
     455           0 :       ELSEIF ( latsys =='aP' .OR. latsys =='triclinic' .OR.
     456             :      &         latsys =='tcl' )  THEN
     457             : 
     458           0 :         noangles=.false.
     459           0 :         i_c = 1
     460           0 :         CALL brvmat ( alpha, beta, gamma, am )
     461           0 :         a1 = a*am(:,1)
     462           0 :         a2 = b*am(:,2)
     463           0 :         a3 = c*am(:,3)
     464           0 :         CALL angles( am )
     465           0 :         scale = 1.0
     466             :           
     467             :       ELSE
     468           0 :           WRITE (errfh,*)
     469           0 :           WRITE (errfh,*) '*** unknown lattice system latsys=',latsys,
     470           0 :      &                    'in line',nline,'.'
     471             :           WRITE (errfh,*) 
     472           0 :      &    '***  No reason for panic, it is probably because'
     473             :           WRITE (errfh,*) 
     474           0 :      & '***  subroutine lattice2 in struct_input.f is unfinished. ***'
     475           0 :         ios = 1
     476           0 :         RETURN
     477             :       ENDIF
     478             : 
     479           0 :       IF ( noangles .AND.
     480             :      &    ( abs(alpha - 90.0 ) > eps .OR.
     481             :      &      abs(beta  - 90.0 ) > eps .OR.
     482             :      &      abs(gamma - 90.0 ) > eps )    ) THEN
     483           0 :         WRITE (errfh,*)
     484           0 :         WRITE (errfh,*) 'ERROR in &lattice ... /. ',
     485           0 :      & 'For the given lattice system all angles should be 90deg.'
     486           0 :         WRITE (errfh,*)
     487             :       ENDIF
     488             : 
     489           0 :       IF (abs(scale(1)) < eps) THEN
     490             :         CALL juDFT_error("Illegal program section reached!"
     491           0 :      +          ,calledby ="lattice2")
     492           0 :         scale(1) = a
     493           0 :         scale(2) = b
     494           0 :         scale(3) = c
     495             :       ENDIF
     496           0 :       aa = a0
     497             : 
     498             :  911  CONTINUE
     499             : 
     500             :       RETURN
     501             :       END SUBROUTINE lattice2
     502             : !
     503             : !**********************************************************************!
     504             : !
     505           6 :       SUBROUTINE angles( am )
     506             : !----------------------------------------------------------------------!
     507             : !     given an bravais matrix (am), calculate length of basis vectors  !
     508             : !     and angles beween them ; write results to standard output        !
     509             : !----------------------------------------------------------------------!
     510             : 
     511             :       USE m_constants
     512             :       IMPLICIT NONE
     513             : 
     514             :       REAL, INTENT(IN) :: am(3,3)
     515             : 
     516             :       REAL     ca,al(3)
     517             :       INTEGER  i,j,i1,i2
     518             : 
     519           6 :       al = 0.0
     520          24 :       DO j = 1, 3
     521         126 :         DO i = 1, 3
     522          72 :           al(j) = al(j) + am(i,j)*am(i,j)
     523             :         ENDDO
     524          24 :         al(j) = sqrt(al(j))
     525             :       ENDDO
     526             : 
     527          24 :       DO j = 1, 3
     528             :         WRITE (6,'("vector ",i1," : ",3f9.5,5x," length : ",f9.5)') 
     529          24 :      &                                              j,am(:,j),al(j)
     530             :       ENDDO
     531             : 
     532          18 :       DO i1 = 1, 2 
     533          36 :         DO i2 = i1+1, 3
     534          18 :           ca = 0.0
     535          72 :           DO i = 1, 3
     536          72 :             ca = ca + am(i,i1)*am(i,i2)
     537             :           ENDDO
     538          18 :           ca  = ca/(al(i1)*al(i2))
     539          18 :           ca = acos(ca)*180/pi_const
     540             :     
     541             :           WRITE (6,'("angle between vectors (",i1,",",i1,") =",f9.5)') 
     542          30 :      &                                                        i1,i2,ca
     543             :         ENDDO
     544             :       ENDDO
     545             : 
     546           6 :       END SUBROUTINE angles
     547             : !
     548             : !**********************************************************************!
     549             : !
     550           0 :       SUBROUTINE brvmat ( alpha, beta, gamma, am )
     551             : !----------------------------------------------------------------------!
     552             : !     given the angles alpha, beta and gamma, set up an matrix 'am'    !
     553             : !     with 3 unit vectors and these angles between them. The first     !
     554             : !     unit vector poins in (1,0,0) direction                gb`05      !
     555             : !----------------------------------------------------------------------!
     556             : 
     557             :       USE m_constants
     558             :       IMPLICIT NONE
     559             : 
     560             :       REAL, INTENT (IN)  :: alpha, beta, gamma
     561             :       REAL, INTENT (OUT) :: am(3,3)
     562             : 
     563             :       REAL ca,cb,cg,sg,c1,c2
     564             : 
     565           0 :       ca = cos(alpha*pi_const/180); cb = cos(beta*pi_const/180)
     566           0 :       cg = cos(gamma*pi_const/180); sg = sin(gamma*pi_const/180)
     567           0 :       c1 = (ca - cg*cb ) / sg
     568           0 :       c2 = sqrt( 1 - cb**2 - c1**2 ) 
     569             : 
     570           0 :       am(1,1) = 1.0 ; am(2,1) = 0.0 ; am(3,1) = 0.0
     571           0 :       am(1,2) = cg  ; am(2,2) = sg  ; am(3,2) = 0.0
     572           0 :       am(1,3) = cb  ; am(2,3) = c1  ; am(3,3) = c2
     573             : 
     574           0 :       END SUBROUTINE brvmat
     575             : 
     576             :       END MODULE m_lattice

Generated by: LCOV version 1.13