C--------------------------------------------------------------%%%%
C-->     Subroutine "derive2d" is used to calculate the elastic%%%%
C-->     constants by quadratic polynomial fit according the   %%%%
C-->     file "OSZICAR" obtained by first principle            %%%%
C-->     calculation for 2D material                           %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine extract2d

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

        character*10 inpos

        inpos='RECELL'

        call symmetry(spgnum, inpos)
        call Astruc2d(area)

C       Oblique system
C       2D space group: p1, p2        
C       (6 independent elastic constants)      

        if (spgnum .GE. 3 .AND. spgnum .LE. 15) then
                nelastic=6
                spgmode=1
        endif

C       Rectanglar system
C       2D space group: pm, pg, p2mm, p2gg, p2mg, cm, c2mm        
C       (4 independent elastic constants)  

        if (spgnum .GE. 16 .AND. spgnum .LE. 74) then
                nelastic=4
                spgmode=2
        endif

C       Square system
C       2D space group: p4, p4mm, p4gm        
C       (3 independent elastic constants)

        if (spgnum .GE. 75 .AND. spgnum .LE. 142) then
                nelastic=3
                spgmode=3
        endif

C       Hexagon system
C       2D space group: p3, p6, p3m1, p31m, p6mm        
C       (2 independent elastic constants)  

        if (spgnum .GE. 143 .AND. spgnum .LE. 194) then
                nelastic=2
                spgmode=4
        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: N/m
C              dataxy(i,j)=(dataxy(i,j)-E0)/area*16.02000
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
C             Unit: N/m
              dataxy(i,j)=(engxy(i,j)-E0)/area*16.02000
           enddo
        enddo

        call writeeng4extract( dataxx,
     &                         engxy,
     &                         dataxy,
     &                         area,
     &                         2,
     &                         spgnum,
     &                         spgmode,
     &                         nelastic)

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