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 derive2d

        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


        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

           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,
     &                 area,
     &                 engerror,
     &                 2,
     &                 spgnum,
     &                 spgmode,
     &                 nelastic)

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

        select case (spgmode)
               case(1)
                 call Tri2d(confxy)
               case(2)
                 call Rect2d(confxy)
               case(3)
                 call Squa2d(confxy)
               case(4)
                 call Hexa2d(confxy)
        end select

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

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                       Oblique system
C%%%%%               (6 independent elastic constants)
C%%%%%                        c11  c12  c16
C%%%%%                        c12  c22  c26
C%%%%%                        c16  c26  c66 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
        subroutine Tri2d(confxy)
        integer i, j
        dimension confxy(6)
        dimension Celas(3,3), Selas(3,3)
        double precision confxy, Celas, Selas
        
        do i=1,3
           do j=1,3
              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(1,2)=confxy(4)-confxy(1)-confxy(2)
        Celas(1,3)=confxy(5)-confxy(1)-confxy(3)
        Celas(2,3)=confxy(6)-confxy(2)-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        write(19,*) 'The elastic constants is that: '
C        write(19,*)
C        write(19,*) 'c11 c22 c66 c12 c16 c26 ' 
C        write(19,'(6f8.2)') Celas(1,1), Celas(2,2), Celas(3,3),
C     &             Celas(1,2), Celas(1,3), Celas(2,3)
C
C        close(19)

        call VRH2d(Celas)

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Rectanglar system
C%%%%%               (4 independent elastic constants)
C%%%%%                        c11  c12    0
C%%%%%                        c12  c22    0
C%%%%%                          0    0  c66
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        subroutine Rect2d(confxy)
        integer i, j
        dimension confxy(6)
        dimension Celas(3,3), Selas(3,3)
        double precision confxy, Celas, Selas
        
        do i=1,3
           do j=1,3
              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(1,2)=confxy(4)-confxy(1)-confxy(2)

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

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 c66 c12 ' 
C        write(19,'(4f8.2)') Celas(1,1), Celas(2,2), Celas(3,3),
C     &             Celas(1,2)
C
C        close(19)

        call VRH2d(Celas)

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                        Square system
C%%%%%               (3 independent elastic constants)
C%%%%%                        c11  c12    0
C%%%%%                        c12  c11    0
C%%%%%                          0    0  c66 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
        subroutine Squa2d(confxy)
        integer i, j
        dimension confxy(6)
        dimension Celas(3,3), Selas(3,3)
        double precision confxy, Celas, Selas
        
        do i=1,3
           do j=1,3
              Celas(i,j)=0.0000
           enddo
        enddo 
        
        Celas(1,1)=2*confxy(3)
        Celas(2,2)=Celas(1,1)
        Celas(3,3)=2*confxy(2)

        Celas(1,2)=confxy(1)-2*confxy(3)

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

        call VRH2d(Celas)

        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                       Hexagon system
C%%%%%               (2 independent elastic constants)
C%%%%%                        c11  c12    0
C%%%%%                        c12  c11    0
C%%%%%                          0    0  c66=(c11-c12)/2 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine Hexa2d(confxy)
        integer i, j
        dimension confxy(6)
        dimension Celas(3,3), Selas(3,3)
        double precision confxy, Celas, Selas
        
        do i=1,3
           do j=1,3
              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)=2.00*confxy(2)

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

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

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

        call VRH2d(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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine VRH2d(Celas)
        integer i, j, L
        integer Stable
        dimension Celas(3,3), Selas(3,3)
        double precision Celas, Selas
        dimension IS(3),JS(3)
        double precision IS, JS
        double precision Ex, Ey, Gxy, Muxy, Muyx

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

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

C       Ex=(c11*c22-c12*c21)/c22
C       Ey=(c11*c22-c12*c21)/c11
C       Gxy=c66        

        Ex=(Celas(1,1)*Celas(2,2)-Celas(2,1)*Celas(1,2))/Celas(2,2)
        Ey=(Celas(1,1)*Celas(2,2)-Celas(2,1)*Celas(1,2))/Celas(1,1)
        Gxy=Celas(3,3)

C       Muxy=c21/c22
C       Muyx=c12/c11        

        Muxy=Celas(2,1)/Celas(2,2)
        Muyx=Celas(1,2)/Celas(1,1)        

C%%%%%  Elastic Stability Conditions

        call stability2d( Stable, Celas )        

C       write the elastic data

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

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

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

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

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

        write(49,*)
        write(49,*) 'Young(Ex and Ey) and shear(Gxy) moduli'
        write(49,*)
        write(49,'(a8,f8.2)') '   Ex : ', Ex
        write(49,'(a8,f8.2)') '   Ey : ', Ey
        write(49,'(a8,f8.2)') '   Gxy: ', Gxy

        write(49,*)
        write(49,*) 'Poisson ratios(Muxy and Muyx)'
        write(49,*)
        write(49,'(a9,f8.3)')  '   Muxy: ', Muxy
        write(49,'(a9,f8.3)')  '   Muyx: ', Muyx


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 stability2d( Stable, Celas )
        integer i, j, Stable
        dimension Celas(3,3)
        dimension Eig(3,3), Vector(3,3)
        double precision EPS
        double precision Celas, Eig, Vector
        
        do i=1,3
           do j=1,3
              Eig(i,j)=Celas(i,j)
           enddo
        enddo

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

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

        Stable=0
        do i=1,3
           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,3
C                   Write(49,'(6f10.2)') Eig(j,j)
C                enddo

C                close(49)

                Stable=1   

           endif
        enddo

        end        

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                    End of derive2d.f                    %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
