C--------------------------------------------------------------%%%%
C-->     Subroutine "derive" is used to calculate the elastic  %%%%
C-->     constants by quadratic polynomial fit according the   %%%%
C-->     file "OSZICAR" obtained by first principle            %%%%
C-->     calculation                                           %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine derive

        use constant
        implicit none             
        
        integer i, j, k, spgmode, nelastic, countnum
        integer spgnum, engerror
        dimension engerror(21)
        double precision aver_x, aver_y
        double precision minpoint, engxy, slopey
        double precision datax, dataxx, datay, dataxy, fitxy
        dimension engxy(21,ndef)
        dimension datax(ndef), dataxx(ndef), datay(ndef), slopey(ndef)
        dimension dataxy(21,ndef),fitxy(21,ndef)
        double precision confxy, tempxy
        dimension confxy(21), tempxy(3)
        double precision E0, volume 
        double precision DT1, DT2, DT3

        character*10 inpos

        inpos='RECELL'

        call symmetry(spgnum, inpos)
        call Vstruc(volume)

C       Triclinic system
C       Point group classes: 1, -1
C       (21 independent elastic constants)        

        if (spgnum .GE. 1 .AND. spgnum .LE. 2) then
                nelastic=21
                spgmode=1
        endif

C       Monoclinic system
C       Point group classes: m, 2, 2/m
C       (13 independent elastic constants)  

        if (spgnum .GE. 3 .AND. spgnum .LE. 15) then
                nelastic=13
                spgmode=2
        endif

C       Orthorhombic system
C       Point group classes: 222, mm2, 2/m2/m2/m
C       (9 independent elastic constants)         

        if (spgnum .GE. 16 .AND. spgnum .LE. 74) then
                nelastic=9
                spgmode=3
        endif

C       Tetragonal II system
C       Point group classes: 4, -4, 4/m
C       (7 independent elastic constants)         

        if (spgnum .GE. 75 .AND. spgnum .LE. 88) then
                nelastic=7
                spgmode=4
        endif

C       Tetragonal I system
C       Point group classes: 422, 4mm, -42m, 4/m2/m/2/m
C       (6 independent elastic constants)         

        if (spgnum .GE. 89 .AND. spgnum .LE. 142) then
                nelastic=6
                spgmode=5
        endif

C       Trigonal II system
C       Point group classes: 3, -3
C       (7 independent elastic constants)         

        if (spgnum .GE. 143 .AND. spgnum .LE. 148) then
                nelastic=7
                spgmode=6
        endif

C       Trigonal I system
C       Point group classes: 32, 3m, -32/m
C       (6 independent elastic constants)         

        if (spgnum .GE. 149 .AND. spgnum .LE. 167) then
                nelastic=6
                spgmode=7
        endif

C       Hexagonal system
C       Point group classes: 6, -6, 6/m, 622, 6mm, -62m, 6/m2/m2/m
C       (5 independent elastic constants)         

        if (spgnum .GE. 168 .AND. spgnum .LE. 194) then
                nelastic=5
                spgmode=8
        endif

C       Cubic system
C       Point group classes: 23, 2/m-3, 432, -43m, 4/m-32/m
C       (3 independent elastic constants)         

        if (spgnum .GE. 195 .AND. spgnum .LE. 230) then
                nelastic=3
                spgmode=9
        endif

        do i=1,ndef
           dataxx(i)=defdata(i)
        enddo

C        open(10,file='elastic.data')
C        do i=1,nelastic
C           read(10,*) (dataxy(i,j),j=1,ndef)
C           E0=minval(dataxy)
C           do j=1,ndef
C              ! Unit: GPa
C              dataxy(i,j)=(dataxy(i,j)-E0)/volume*160.2000
C           enddo
C        enddo
C        close(10)

        call readoszi(engxy, nelastic)

        do i=1,nelastic
           E0=minval(engxy(i,1:ndef))
           do j=1,ndef
              ! Unit: GPa
              dataxy(i,j)=(engxy(i,j)-E0)/volume*160.2000
           enddo
        enddo


        do i=1,nelastic

           ! To delete the bad points
           k=1
           do j=1,ndef
              if (dataxy(i,j) .GE. 0.0 .AND. dataxy(i,j) .LE. 1.0) then
                   datay(k)=dataxy(i,j)
                   datax(k)=dataxx(j)
                   k=k+1
              endif
           enddo

           ! To avoid not enough data ( <= 3 )

           ! if ( k-1 .LE. 3 ) THEN
           !         write(*,*) 'Not enough data!!!'
           ! endif

           call HPIR(datax,datay,tempxy,k-1,3,DT1,DT2,DT3)
           confxy(i)=tempxy(3)
           
C          To solve the fitness data of dataxy

           aver_x=sum(datax(1:(k-1)))/(k-1)

           do j=1,ndef
              fitxy(i,j)=tempxy(3)*(dataxx(j)-aver_x)**2
     &                  +tempxy(2)*(dataxx(j)-aver_x)
     &                  +tempxy(1)
           enddo

C          To avoid the unstabled distortion

           engerror(i)=0
           
C           ! WARNING!! NOT positive quadratic curve
C           if ( confxy(i) .LT. 0.0 ) THEN
C                   engerror(i)=3
C           endif

           ! WARNING!! The lowest energy is < -0.004 or > 0.004 
           minpoint=-1/2*tempxy(2)/tempxy(3)
           if ( minpoint .LE. -0.004 .OR. minpoint .GE. 0.004 ) THEN
                   engerror(i)=1
           endif

           ! WARNING!! NOT quadratic curve or BAD points
           do j=2,(k-1)
              slopey(j-1)=(datay(j)-datay(j-1))/(datax(j)-datax(j-1))
           enddo

           countnum=0
           do j=2,(k-1)
              if ( slopey(j-1)*slopey(j) .LT. 0 ) THEN
                      countnum=countnum+1
              endif
           enddo

           if ( countnum .GE. 3 ) THEN
                   engerror(i)=2
           endif

           ! ERROR!! Explained Sum of Squares > 0.1!!!
           if (DT1 .GE. 0.1) then
                engerror(i)=-1
           endif           

        enddo

        call writeeng( dataxx,
     &                 engxy,
     &                 dataxy,
     &                 fitxy,
     &                 volume,
     &                 engerror,
     &                 3,
     &                 spgnum,
     &                 spgmode,
     &                 nelastic)

        do i=1,nelastic
           if ( engerror(i) .NE. 0 ) then
                stop
           endif
        enddo

        select case (spgmode)
               case(1)
                 call Tric(confxy)
               case(2)
                 call Monoc(confxy)
               case(3)
                 call Orthor(confxy)
               case(4)
                 call Tetrag7(confxy)
               case(5)
                 call Tetrag6(confxy)
               case(6)
                 call Trig7(confxy)
               case(7)
                 call Trig6(confxy)
               case(8)
                 call Hexag(confxy)
               case(9)
                 call Cubic(confxy)                
        end select

        end
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                 End of mian program                     %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Triclinic system
C%%%%%               (21 independent elastic constants)
C%%%%%                 c11  c12  c13  c14  c15  c16
C%%%%%                 c12  c22  c23  c24  c25  c26
C%%%%%                 c13  c23  c33  c34  c35  c36
C%%%%%                 c14  c24  c34  c44  c45  c46
C%%%%%                 c15  c25  c35  c45  c55  c56
C%%%%%                 c16  c26  c36  c46  c56  c66
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       
        subroutine Tric(confxy)
        integer i, j
        dimension confxy(21)
        dimension Celas(6,6), Selas(6,6)
        double precision confxy, Celas, Selas
        
        do i=1,6
           do j=1,6
              Celas(i,j)=0.0000
           enddo
        enddo 
        
        Celas(1,1)=2*confxy(1)
        Celas(2,2)=2*confxy(2)
        Celas(3,3)=2*confxy(3)
        Celas(4,4)=2*confxy(4)
        Celas(5,5)=2*confxy(5)
        Celas(6,6)=2*confxy(6)

        Celas(1,2)=confxy(7)-confxy(1)-confxy(2)
        Celas(1,3)=confxy(8)-confxy(1)-confxy(3)
        Celas(1,4)=confxy(9)-confxy(1)-confxy(4)
        Celas(1,5)=confxy(10)-confxy(1)-confxy(5)        
        Celas(1,6)=confxy(11)-confxy(1)-confxy(6)       
        Celas(2,3)=confxy(12)-confxy(2)-confxy(3)
        Celas(2,4)=confxy(13)-confxy(2)-confxy(4)
        Celas(2,5)=confxy(14)-confxy(2)-confxy(5)        
        Celas(2,6)=confxy(15)-confxy(2)-confxy(6)
        Celas(3,4)=confxy(16)-confxy(3)-confxy(4)
        Celas(3,5)=confxy(17)-confxy(3)-confxy(5)
        Celas(3,6)=confxy(18)-confxy(3)-confxy(6)        
        Celas(4,5)=confxy(19)-confxy(4)-confxy(5)
        Celas(4,6)=confxy(20)-confxy(4)-confxy(6)
        Celas(5,6)=confxy(21)-confxy(5)-confxy(6)

        Celas(2,1)=Celas(1,2)
        Celas(3,1)=Celas(1,3)
        Celas(4,1)=Celas(1,4)
        Celas(5,1)=Celas(1,5)
        Celas(6,1)=Celas(1,6)
        Celas(3,2)=Celas(2,3)
        Celas(4,2)=Celas(2,4)
        Celas(5,2)=Celas(2,5)
        Celas(6,2)=Celas(2,6)
        Celas(4,3)=Celas(3,4)
        Celas(5,3)=Celas(3,5)
        Celas(6,3)=Celas(3,6)
        Celas(5,4)=Celas(4,5)
        Celas(6,4)=Celas(4,6)
        Celas(6,5)=Celas(5,6)        

C        open(unit=19,position='Append',FILE='RELAS')          

C        write(19,*) 'The elastic constants is that: '
C        write(19,*)
C        write(19,'(a83)') 'c11 c22 c33 c44 c55 c66 c12 c13 c14 c15
C     & c16 c23 c24 c25 c26 c34 c35 c36 c45 c46 c56'
C        write(19,'(21f8.2)') Celas(1,1), Celas(2,2), Celas(3,3),
C     &             Celas(4,4), Celas(5,5), Celas(6,6), Celas(1,2),
C     &             Celas(1,3), Celas(1,4), Celas(1,5), Celas(1,6),
C     &             Celas(2,3), Celas(2,4), Celas(2,5), Celas(2,6),
C     &             Celas(3,4), Celas(3,5), Celas(3,6), Celas(4,5),
C     &             Celas(4,6), Celas(5,6)
C
C        close(19)

        call VRH(Celas)

        end        

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Monoclinic system
C%%%%%               (13 independent elastic constants)
C%%%%%                 c11  c12  c13    0  c15    0
C%%%%%                 c12  c22  c23    0  c25    0
C%%%%%                 c13  c23  c33    0  c35    0
C%%%%%                   0    0    0  c44    0  c46
C%%%%%                 c15  c25  c35    0  c55    0
C%%%%%                   0    0    0  c46    0  c66
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       
        subroutine Monoc(confxy)
        integer i, j
        dimension confxy(20)
        dimension Celas(6,6), Selas(6,6)
        double precision confxy, Celas, Selas
        
        do i=1,6
           do j=1,6
              Celas(i,j)=0.0000
           enddo
        enddo 
        
        Celas(1,1)=2*confxy(1)
        Celas(2,2)=2*confxy(2)
        Celas(3,3)=2*confxy(3)
        Celas(4,4)=2*confxy(4)
        Celas(5,5)=2*confxy(5)
        Celas(6,6)=2*confxy(6)

        Celas(1,2)=confxy(7)-confxy(1)-confxy(2)
        Celas(2,3)=confxy(8)-confxy(2)-confxy(3)
        Celas(1,3)=confxy(9)-confxy(1)-confxy(3)
        Celas(1,5)=confxy(10)-confxy(1)-confxy(5)
        Celas(2,5)=confxy(11)-confxy(2)-confxy(5)
        Celas(3,5)=confxy(12)-confxy(3)-confxy(5)
        Celas(4,6)=confxy(13)-confxy(4)-confxy(6)

        Celas(2,1)=Celas(1,2)
        Celas(3,2)=Celas(2,3)
        Celas(3,1)=Celas(1,3)
        Celas(5,1)=Celas(1,5)
        Celas(5,2)=Celas(2,5)
        Celas(5,3)=Celas(3,5)
        Celas(6,4)=Celas(4,6)

C        open(unit=19,position='Append',FILE='RELAS')  

C        write(19,*) 'The elastic constants is that: '
C        write(19,*)
C        write(19,*) 'c11 c22 c33 c44 c55 c66 c12 c13 c15 c23 ', 
C     &             'c25 c35 c46'
C        write(19,'(13f8.2)') Celas(1,1), Celas(2,2), Celas(3,3),
C     &             Celas(4,4), Celas(5,5), Celas(6,6), Celas(1,2),
C     &             Celas(1,3), Celas(1,5), Celas(2,3), Celas(2,5),
C     &             Celas(3,5), Celas(4,6)
C
C        close(19)

        call VRH(Celas)

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Orthorhombic system
C%%%%%               (9 independent elastic constants)
C%%%%%                 c11  c12  c13    0    0    0
C%%%%%                 c12  c22  c23    0    0    0
C%%%%%                 c13  c23  c33    0    0    0
C%%%%%                   0    0    0  c44    0    0
C%%%%%                   0    0    0    0  c55    0
C%%%%%                   0    0    0    0    0  c66
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        subroutine Orthor(confxy)
        integer i, j
        dimension confxy(20)
        dimension Celas(6,6), Selas(6,6)
        double precision confxy, Celas, Selas
        
        do i=1,6
           do j=1,6
              Celas(i,j)=0.0000
           enddo
        enddo 
        
        Celas(1,1)=2*confxy(1)
        Celas(2,2)=2*confxy(2)
        Celas(3,3)=2*confxy(3)
        Celas(4,4)=2*confxy(4)
        Celas(5,5)=2*confxy(5)
        Celas(6,6)=2*confxy(6)

        Celas(1,2)=confxy(7)-confxy(1)-confxy(2)
        Celas(2,3)=confxy(8)-confxy(2)-confxy(3)
        Celas(1,3)=confxy(9)-confxy(1)-confxy(3)

        Celas(2,1)=Celas(1,2)
        Celas(3,2)=Celas(2,3)
        Celas(3,1)=Celas(1,3)

C        open(unit=19,position='Append',FILE='RELAS')  
C
C        write(19,*) 'The elastic constants is that: '
C        write(19,*)
C        write(19,*) 'c11 c22 c33 c44 c55 c66 c12 c13 c23'
C        write(19,'(9f8.2)') Celas(1,1), Celas(2,2), Celas(3,3),
C     &             Celas(4,4), Celas(5,5), Celas(6,6), Celas(1,2),
C     &             Celas(1,3), Celas(2,3)
C
C        close(19)

        call VRH(Celas)

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Tetragonal II system
C%%%%%               (7 independent elastic constants)
C%%%%%                 c11  c12  c13    0    0  c16
C%%%%%                 c12  c11  c13    0    0 -c16
C%%%%%                 c13  c13  c33    0    0    0
C%%%%%                   0    0    0  c44    0    0
C%%%%%                   0    0    0    0  c44    0
C%%%%%                 c16 -c16    0    0    0  c66
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        subroutine Tetrag7(confxy)
        integer i, j
        dimension confxy(20)
        dimension Celas(6,6), Selas(6,6)
        double precision confxy, Celas, Selas
        
        do i=1,6
           do j=1,6
              Celas(i,j)=0.0000
           enddo
        enddo 
        
        Celas(1,1)=confxy(1)-(confxy(5)-2.00*confxy(6)+confxy(3))
        Celas(2,2)=Celas(1,1)
        Celas(3,3)=confxy(3)*2.00
        Celas(4,4)=confxy(4)
        Celas(5,5)=Celas(4,4)
        Celas(6,6)=confxy(2)*2.00

        Celas(1,2)=confxy(5)-2.00*confxy(6)+confxy(3)
        Celas(1,3)=(confxy(5)-confxy(3)-confxy(1))/2.00
        Celas(1,6)=confxy(7)-Celas(1,1)/2.00-Celas(6,6)/2.00
        Celas(2,3)=Celas(1,3)
        Celas(2,6)=-Celas(1,6)

        Celas(2,1)=Celas(1,2)
        Celas(3,2)=Celas(2,3)
        Celas(3,1)=Celas(1,3)
        Celas(6,1)=Celas(1,6)
        Celas(6,2)=Celas(2,6)

C        open(unit=19,position='Append',FILE='RELAS')  
C
C        write(19,*) 'The elastic constants is that: '
C        write(19,*)
C        write(19,*) 'c11 c33 c44 c66 c12 c13 c16'
C        write(19,'(7f8.2)') Celas(1,1), Celas(3,3), Celas(4,4),
C     &             Celas(6,6), Celas(1,2), Celas(1,3), Celas(1,6)
C
C        close(19)

        call VRH(Celas)

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Tetragonal I system
C%%%%%               (6 independent elastic constants)
C%%%%%                 c11  c12  c13    0    0    0
C%%%%%                 c12  c11  c13    0    0    0
C%%%%%                 c13  c13  c33    0    0    0
C%%%%%                   0    0    0  c44    0    0
C%%%%%                   0    0    0    0  c44    0
C%%%%%                   0    0    0    0    0  c66
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
        subroutine Tetrag6(confxy)
        integer i, j
        dimension confxy(20)
        dimension Celas(6,6), Selas(6,6)
        double precision confxy, Celas, Selas
        
        do i=1,6
           do j=1,6
              Celas(i,j)=0.0000
           enddo
        enddo 
        
        Celas(1,1)=confxy(1)-(confxy(5)-2.00*confxy(6)+confxy(3))
        Celas(2,2)=Celas(1,1)
        Celas(3,3)=confxy(3)*2.00
        Celas(4,4)=confxy(4)
        Celas(5,5)=Celas(4,4)
        Celas(6,6)=confxy(2)*2.00

        Celas(1,2)=confxy(5)-2.00*confxy(6)+confxy(3)
        Celas(1,3)=(confxy(5)-confxy(3)-confxy(1))/2.00
        Celas(2,3)=Celas(1,3)       
        
        Celas(2,1)=Celas(1,2)
        Celas(3,2)=Celas(2,3)
        Celas(3,1)=Celas(1,3)

C        open(unit=19,position='Append',FILE='RELAS')  
C
C        write(19,*) 'The elastic constants is that: '
C        write(19,*)
C        write(19,*) 'c11 c33 c44 c66 c12 c13'
C        write(19,'(6f8.2)') Celas(1,1), Celas(3,3), Celas(4,4),
C     &             Celas(6,6), Celas(1,2), Celas(1,3)
C
C        close(19)

        call VRH(Celas)

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Rhombohedral II system
C%%%%%               (7 independent elastic constants)
C%%%%%                 c11  c12  c13  c14  c15    0
C%%%%%                 c12  c11  c13 -c14 -c15    0
C%%%%%                 c13  c13  c33    0    0    0
C%%%%%                 c14 -c14    0  c44    0 -c15
C%%%%%                 c15 -c15    0    0  c44  c14
C%%%%%                   0    0    0 -c15  c14  c66=(c11-c12)/2
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
        subroutine Trig7(confxy)
        integer i, j
        dimension confxy(20)
        dimension Celas(6,6), Selas(6,6)
        double precision confxy, Celas, Selas
        
        do i=1,6
           do j=1,6
              Celas(i,j)=0.0000
           enddo
        enddo 
        
        Celas(1,1)=(confxy(1)+4*confxy(2))/2.00
        Celas(2,2)=Celas(1,1)
        Celas(3,3)=confxy(3)*2
        Celas(4,4)=confxy(4)
        Celas(5,5)=Celas(4,4)
        Celas(6,6)=2.00*confxy(2)

        Celas(1,2)=(confxy(1)-4*confxy(2))/2.00
        Celas(1,3)=(confxy(5)-confxy(1)-confxy(3))/2.00
        Celas(2,3)=Celas(1,3)
        Celas(1,4)=confxy(6)-confxy(4)/2.00-confxy(2)
        Celas(1,5)=-(confxy(7)-Celas(1,1)/2.00-Celas(4,4)/2.00)
        Celas(4,6)=-Celas(1,5)
        Celas(2,4)=-Celas(1,4)
        Celas(2,5)=-Celas(1,5)
        Celas(5,6)=Celas(1,4)

        Celas(2,1)=Celas(1,2)
        Celas(3,2)=Celas(2,3)
        Celas(3,1)=Celas(1,3)
        Celas(4,1)=Celas(1,4)
        Celas(5,1)=Celas(1,5)
        Celas(4,2)=Celas(2,4)
        Celas(5,2)=Celas(2,5)
        Celas(6,5)=Celas(5,6)
        Celas(6,4)=Celas(4,6)

C        open(unit=19,position='Append',FILE='RELAS')  
C
C        write(19,*) 'The elastic constants is that: '
C        write(19,*)
C        write(19,*) 'c11 c33 c44 c12 c13 c14 c15'
C        write(19,'(7f8.2)') Celas(1,1), Celas(3,3), Celas(4,4),
C     &             Celas(1,2), Celas(1,3), Celas(1,4), Celas(1,5)
C
C        close(19)

        call VRH(Celas)

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Rhombohedral I system
C%%%%%               (6 independent elastic constants)
C%%%%%                 c11  c12  c13  c14    0    0
C%%%%%                 c12  c11  c13 -c14    0    0
C%%%%%                 c13  c13  c33    0    0    0
C%%%%%                 c14 -c14    0  c44    0    0
C%%%%%                   0    0    0    0  c44  c14
C%%%%%                   0    0    0    0  c14  c66=(c11-c12)/2
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
        subroutine Trig6(confxy)
        integer i, j
        dimension confxy(20)
        dimension Celas(6,6), Selas(6,6)
        double precision confxy, Celas, Selas
        
        do i=1,6
           do j=1,6
              Celas(i,j)=0.0000
           enddo
        enddo 
        
        Celas(1,1)=(confxy(1)+4.00*confxy(2))/2.00
        Celas(2,2)=Celas(1,1)
        Celas(3,3)=confxy(3)*2.00
        Celas(4,4)=confxy(4)
        Celas(5,5)=Celas(4,4)
        Celas(6,6)=2.00*confxy(2)

        Celas(1,2)=(confxy(1)-4.00*confxy(2))/2.00
        Celas(1,3)=(confxy(5)-confxy(1)-confxy(3))/2.00
        Celas(1,4)=confxy(6)-confxy(4)/2.00-confxy(2)
        Celas(2,3)=Celas(1,3)
        Celas(2,4)=-Celas(1,4)
        Celas(5,6)=Celas(1,4)

        Celas(2,1)=Celas(1,2)
        Celas(3,2)=Celas(2,3)
        Celas(3,1)=Celas(1,3)
        Celas(4,1)=Celas(1,4)
        Celas(4,2)=Celas(2,4)
        Celas(6,5)=Celas(5,6)

C        open(unit=19,position='Append',FILE='RELAS')
C
C        write(19,*) 'The elastic constants is that: '
C        write(19,*)
C        write(19,*) 'c11 c33 c44 c12 c13 c14'
C        write(19,'(6f8.2)') Celas(1,1), Celas(3,3), Celas(4,4),
C     &                      Celas(1,2), Celas(1,3), Celas(1,4)
C
C        close(19)

        call VRH(Celas)

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                       Hexagonal system
C%%%%%               (5 independent elastic constants)
C%%%%%                 c11  c12  c13    0    0    0
C%%%%%                 c12  c11  c13    0    0    0
C%%%%%                 c13  c13  c33    0    0    0
C%%%%%                   0    0    0  c44    0    0
C%%%%%                   0    0    0    0  c44    0
C%%%%%                   0    0    0    0    0  c66=(c11-c12)/2
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine Hexag(confxy)
        integer i, j
        dimension confxy(20)
        dimension Celas(6,6), Selas(6,6)
        double precision confxy, Celas, Selas
        
        do i=1,6
           do j=1,6
              Celas(i,j)=0.0000
           enddo
        enddo 
        
        Celas(1,1)=(confxy(1)+4.00*confxy(2))/2.00
        Celas(2,2)=Celas(1,1)
        Celas(3,3)=confxy(3)*2.00
        Celas(4,4)=confxy(4)
        Celas(5,5)=Celas(4,4)
        Celas(6,6)=2.00*confxy(2)

        Celas(1,2)=(confxy(1)-4.00*confxy(2))/2.00
        Celas(1,3)=(confxy(5)-confxy(1)-confxy(3))/2.00
        Celas(2,3)=Celas(1,3)

        Celas(2,1)=Celas(1,2)
        Celas(3,2)=Celas(2,3)
        Celas(3,1)=Celas(1,3)

C        open(unit=19,position='Append',FILE='RELAS')  
C
C        write(19,*) 'The elastic constants is that: '
C        write(19,*)
C        write(19,*) 'c11 c33 c44 c12 c13'
C        write(19,'(5f8.2)') Celas(1,1), Celas(3,3), Celas(4,4),
C     &                      Celas(1,2), Celas(1,3)

C        close(19)

        call VRH(Celas)

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                          Cubic system
C%%%%%               (3 independent elastic constants)
C%%%%%                 c11  c12  c12    0    0    0
C%%%%%                 c12  c11  c12    0    0    0
C%%%%%                 c12  c12  c11    0    0    0
C%%%%%                   0    0    0  c44    0    0
C%%%%%                   0    0    0    0  c44    0
C%%%%%                   0    0    0    0    0  c44
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine Cubic(confxy)
        integer i, j
        dimension confxy(20)
        dimension Celas(6,6), Selas(6,6)
        double precision confxy, Celas, Selas
        
        do i=1,6
           do j=1,6
              Celas(i,j)=0.0000
           enddo
        enddo 
        
        Celas(1,1)=2.00*confxy(2)-2.00/3.00*confxy(3)
        Celas(2,2)=Celas(1,1)
        Celas(3,3)=Celas(1,1)
        Celas(4,4)=2.00/3.00*confxy(1)
        Celas(5,5)=Celas(4,4)
        Celas(6,6)=Celas(4,4)

        Celas(1,2)=2.00/3.00*confxy(3)-confxy(2)
        Celas(2,3)=Celas(1,2)
        Celas(1,3)=Celas(1,2)

        Celas(2,1)=Celas(1,2)
        Celas(3,2)=Celas(2,3)
        Celas(3,1)=Celas(1,3)

C        open(unit=19,position='Append',FILE='RELAS')  
C
C        write(19,*) 'The elastic constants is that: '
C        write(19,*)
C        write(19,*) 'c11 c44 c12'
C        write(19,'(3f8.2)') Celas(1,1), Celas(4,4), Celas(1,2)
C
C        close(19)

        call VRH(Celas)

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%     subroutine "VHR" is used to calculation the         %%%%
C%%%%%     elasticity modulus according elastic constants      %%%%
C%%%%%     using the Voigt-Reuss-Hill approximations           %%%%
C%%%%%                        REFERENCE                        %%%%
C%%%%%     1\ Hill R. The elastic behavior of crystalline      %%%%
C%%%%%        aggregate[J]. 1952, 65.                          %%%%
C%%%%%     2\ Ranganathan S I, Ostoja-Starzewski M. Universal  %%%%
C%%%%%        elastic anisotropy index[J]. Physical Review     %%%%
C%%%%%        Letters, 2008, 101(5): 055504.                   %%%%
C%%%%%     3\ Chung D H, Buessem W R. The elastic anisotropy   %%%%
C%%%%%        of crystals[J]. Journal of Applied Physics,      %%%%
C%%%%%        1967, 38(5): 2010-2012.                          %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine VRH(Celas)
        integer i, j, L
        integer Stable
        dimension Celas(6,6), Selas(6,6)
        double precision Celas, Selas
        dimension IS(6),JS(6)
        double precision IS, JS
        double precision Pv, Qv, Rv, Pr, Qr, Rr
        double precision Ev, Gv, Kv, Muv
        double precision Er, Gr, Kr, Mur
        double precision Eh, Gh, Kh, Muh
        double precision Ac, Au      

        do i=1,6
           do j=1,6
              Selas(i,j)=Celas(i,j)
           enddo
        enddo

        call brinv(Selas,6,L,IS,JS)

C       Young's(E), shear(G) and bulk(K) moduli and Poisson ratio(Mu) 

C       E=9KG/(3K+G)
C       Mu=(3K-2G)/(2(3K+G))

C%%%%%  Voigt

C       9Kv=(c11+c22+c33)+2(c12+c23+c13)
C       15Gv=(c11+c22+c33)-(c12+c23+c13)+3(c44+c55+c66)      

        Pv=Celas(1,1)+Celas(2,2)+Celas(3,3);
        Qv=Celas(1,2)+Celas(1,3)+Celas(2,3);
        Rv=Celas(4,4)+Celas(5,5)+Celas(6,6);

        Ev=((Pv+2.00*Qv)*(Pv-Qv+3.00*Rv))/(3.00*(2.00*Pv+3.00*Qv+Rv));
        Gv=(Pv-Qv+3.00*Rv)/15.00;
        Kv=Ev*Gv/(3.00*(3.00*Gv-Ev));
        Muv=(Ev/(2.00*Gv))-1.00;    

C%%%%%  Reuss

C       1/Kr=(s11+s22+s33)+2(s12+s23+s13)
C       15/Gr=4(s11+s22+s33)-4(s12+s23+s13)+3(s44+s55+s66)

        Pr=Selas(1,1)+Selas(2,2)+Selas(3,3);
        Qr=Selas(1,2)+Selas(1,3)+Selas(2,3);
        Rr=Selas(4,4)+Selas(5,5)+Selas(6,6);

        Er=15.00/(3.00*Pr+2.00*Qr+Rr);
        Gr=15.00/(4.00*(Pr-Qr)+3.00*Rr);
        Kr=Er*Gr/(3.00*(3.00*Gr-Er));
        Mur=(Er/(2.00*Gr))-1.00;

C%%%%%  Hill
        Eh=(Ev+Er)/2.00
        Gh=(Gv+Gr)/2.00
        Kh=(Kv+Kr)/2.00
        Muh=(Muv+Mur)/2.00

C%%%%%  the Pugh ratio (G/B)

!        write(19,*) '>>> the Pugh ratio (G/B) <<<'
!        write(19,"(f8.4)") Gh/Kh

C%%%%%  the Cauchy pressure (Pc=C12-C44)

!        write(19,*) '>>> the Cauchy pressure (Pc=C12-C44) <<<'
!        write(19,"(f8.2)") Celas(1,2)-Celas(4,4)


C%%%%%  Elastic Anisotropy Index

C       Chung-Buessem Anisotropy Index: Ac=(Gv-Gr)/(Gv+Gr)
C       Universal Elastic Anisotropy Index: Au=10*Ac/(1-Ac)

        Ac=(Gv-Gr)/(Gv+Gr)
        Au=5*Gv/Gr+Kv/Kr-6

C%%%%%  Elastic Stability Conditions

        call stability( Stable, Celas )

C       write the elastic data

        open(unit = 49, file = 'ELADAT', status = 'unknown')

        write(49,*) "Elastic tensor: "
        write(49,*) 

        do i=1,6
           write(49,"(6f8.2)") (Celas(i,j), j=1,6)
        enddo       

        write(49,*)
        write(49,*) "Compliance tensor: "
        write(49,*)

        do i=1,6
            write(49,"(6f10.6)") (Selas(i,j), j=1,6)
        enddo 

        write(49,*) 
        write(49,*) "Young's, shear and bulk moduli and Poisson ratio"
        write(49,*)
        write(49,"(a22, 3f8.2, f8.3)") '   Voigt approximate: ', 
     &                                                 Ev, Gv, Kv, Muv

        write(49,"(a22, 3f8.2, f8.3)") '   Reuss approximate: ', 
     &                                                 Er, Gr, Kr, Mur

        write(49,"(a22, 3f8.2, f8.3)") '   Hill approximate : ', 
     &                                                 Eh, Gh, Kh, Muh

        write(49,*)
        write(49,"(a19, f8.4)") ' Pugh ratio (G/K): ', Gh/Kh
        write(49,"(a31, f8.2)") ' Cauchy pressure (Pc=C12-C44): ', 
     &                                           Celas(1,2)-Celas(4,4)

        write(49,*)
        write(49,*) 'Elastic Anisotropy Index'
        write(49,*)
        write(49,'(a34,a22,f8.4)') '   Chung-Buessem Anisotropy Index ',
     &                         '(Ac=(Gv-Gr)/(Gv+Gr)): ', Ac
        write(49,'(a32,a6,a22,f8.4)') '   Universal Elastic Anisotropy '
     &                      , 'Index ', '(Au=5*Gv/Gr+Kv/Kr-6): ', Au


C       Stable=0: Stable
C       Stable=1: UNstable 

        write(49,*)
        if ( Stable .EQ. 0 ) THEN
                write(49,*) 'Elastic Stability Conditions:  Stable'
        elseif ( Stable .EQ. 1 ) THEN
                write(49,*) 'Elastic Stability Conditions:  UNstable'
        endif

        close(49)

        open(unit=19,position='Append',FILE='RELAS') 
        write(19,*) 'SUCCESSFUL'
        close(19)

C       End of writing

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                Elastic Stability Conditions             %%%% 
C%%%%%    The generic necessary and sufficient criterion that  %%%%
C%%%%%    all eigenvalues of C be position                     %%%%
C%%%%%                        REFERENCE                        %%%%
C%%%%%    1\ Mouhat F, Coudert F X. Necessary and sufficient   %%%%
C%%%%%       elastic stability conditions in various crystal   %%%%
C%%%%%       systems[J]. Physical Review B, 2014, 90(22).      %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        subroutine stability( Stable, Celas )
        integer i, j, Stable
        dimension Celas(6,6)
        dimension Eig(6,6), Vector(6,6)
        double precision EPS
        double precision Celas, Eig, Vector
        
        do i=1,6
           do j=1,6
              Eig(i,j)=Celas(i,j)
           enddo
        enddo

        EPS=0.01 
        call cjcbj(Eig,6,EPS,Vector)

C       Stable=0: Stable
C       Stable=1: UNstable        

        Stable=0
        do i=1,6
           if (Eig(i,i) .LE. 0) then
C                open(unit=49,position='Append',FILE='EDATA')
C                   
C                Write(49,*)    
C                Write(49,*) "WARNING!!! The elastic UNstability"
C                Write(49,*) 
C                Write(49,*) "   All eigenvalues of C is:"
C                Write(49,*)
C
C                do j=1,6
C                   Write(49,'(6f10.2)') Eig(j,j)
C                enddo
C       
C                close(49)

                 Stable=1

           endif
        enddo

        end        

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                     End of derive.f                     %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
