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 extract

        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

        call writeeng4extract( dataxx,
     &                         engxy,
     &                         dataxy,
     &                         volume,
     &                         3,
     &                         spgnum,
     &                         spgmode,
     &                         nelastic)


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