C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%   Subroutine "recell2d" is used to rewrite the structure %%%%
C%%%%   file with standardized unit cell using the spglib,     %%%%
C%%%%   which is a library for finding and handling crystal    %%%%
C%%%%   symmetries for 2D materials                            %%%%
C%%%%                        REFERENCE                         %%%%
C%%%%   1\ https://atztogo.github.io/spglib/                   %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C                         Standardized Unit Cell
C 
C
C        >>>>>>>>>>>>>>>>>>  Trianglar system  <<<<<<<<<<<<<<<<<<
C      
C        b axis is taken as the unique axis
C      
C        alpha = 90 degree and bata = 90 degree
C        gamma > 90 degree
C 
C        a > b
C       
C        a is set along +x direction of Cartesian coordinates
C        b is set in x-y plane of Cartesian coordinates      
C        c is set along +z direction of Cartesian coordinates
C        
C
C        >>>>>>>>>>>>>>>>>  Rectanglar system  <<<<<<<<<<<<<<<<<      
C
C      
C        alpha = beta = gamma = 90 degree
C
C        a > b
C      
C        a is set along +x direction of Cartesian coordinates
C        b is set along +y direction of Cartesian coordinates
C        c is set along +z direction of Cartesian coordinates
C
C
C        >>>>>>>>>>>>>>>>>>  Square system  <<<<<<<<<<<<<<<<<<      
C      
C
C        alpha = beta = gamma = 90 degree
C      
C        a = b
C      
C        a is set along +x direction of Cartesian coordinates
C        b is set along +y direction of Cartesian coordinates
C        c is set along +z direction of Cartesian coordinates
C
C      
C        >>>>>>>>>>>>>>>>>>   Hexagon system   <<<<<<<<<<<<<<<<<<      
C
C      
C        alpha = beta = 90 degree
C        gamma = 120 degree
C      
C        a = b
C      
C        a is set along +x direction of Cartesian coordinates
C        b is set in x-y plane of Cartesian coordinates
C        c is set along +z direction of Cartesian coordinates
C
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


C        module defs_basis
C          implicit none
C          integer, parameter :: dp=kind(1.0d0)
C        end module defs_basis

        subroutine recell2d

        use defs_basis
        use spglib_f08
        use constant
        implicit none        

        double precision alat,pripos,privect
        dimension privect(3,3),pripos(1000,3)
        character*10 title,fix
        character*5 sname
        character*10 inpos   
        dimension sname(10),fix(1000,3)
        integer i,j,k,m,ntype,natomi,nn
        integer Ncoord,Ndyn,atom_types,anum
        dimension natomi(10),atom_types(1000),anum(11) 

        character*11 symbol
        character*15 outpos
        character*10 pricoord
        integer posmode
        integer space_group, refine_cell
        real*8 lattice, positions
        dimension lattice(3,3), positions(3,1000)

        double precision temp, temppos, tempvect, recipvect
        double precision lega, legb, legc
        double precision legar, legbr, legcr
        double precision cos_alpha, cos_bata, cos_gamma
        double precision cos_alphar, cos_batar, cos_gammar
        double precision cos_angle
        dimension temppos(1000,3), tempvect(5,5), recipvect(3,3)

        inpos='INPOS'

        call readpos(inpos,
     &               title, 
     &               alat, 
     &               privect, 
     &               sname,
     &               natomi,               
     &               ntype,
     &               nn,  
     &               Ndyn,
     &               Ncoord,
     &               pripos,
     &               fix,
     &               atom_types)

C       output the initial structure

        if (Ncoord .EQ. 0) THEN
                pricoord='Direct'
        else if (Ncoord .EQ. 1) THEN
                pricoord='Cartesian'
        else
                open(unit=19,position='Append',FILE='RELAS') 
                write(19,*) 'ERROR!! NOT DIRECT OR CARTESIAN!!!'
                close(19)
                stop
        endif

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

        write(19,'(a22)') 'The initial structure:'
        write(19,*)

          write(19,'(a4,f8.5)') '    ',alat
          do i=1,3
             write(19,"(3f20.10)")(privect(i,j),j=1,3)
          enddo
          write(19,'(a4,10A5)') '    ',(sname(i),i=1,ntype)
          write(19,'(a4,10I4)') '    ',(natomi(i), i=1,ntype)
          write(19,'(a4,A6)') '    ',pricoord
          do i=1, nn
             write(19,"(3f20.10)") (pripos(i,j),j=1,3)
          enddo
        
        close(19)        

C       end of outputing the initial structure          

        if (Ncoord .EQ. 1 .OR. Ncoord .EQ. 0) THEN
                call directpos(alat,
     &                        privect,
     &                        nn,
     &                        Ncoord,
     &                        pripos)
                Ncoord=0
        else
                open(unit=19,position='Append',FILE='RELAS')
                write(19,*) 'ERROR!! NOT DIRECT OR CARTESIAN!!!'
                close(19)
                stop
        endif


C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%    To get the space group of the structure "INPOS"      %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        do i=1,3
           do j=1,3
              lattice(i,j)=privect(i,j)
           enddo
        enddo

        do i=1,nn
           do j=1,3
              positions(j,i)=pripos(i,j)
           enddo
        enddo

        space_group = spg_get_international( symbol, lattice, 
     &                positions, atom_types, nn, symprec )

        if ( space_group .EQ. 0) then
                open(unit=19,position='Append',FILE='RELAS')
                write(19,*) 'ERROR!! Cannot get space group!!!'
                close(19)
                stop
        endif

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%%%%%   To refine the new cell with standardized unit cell %%% 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        refine_cell = spg_refine_cell( lattice, positions, atom_types,
     &                nn, symprec)

        if ( refine_cell .EQ. 0) then
                open(unit=19,position='Append',FILE='RELAS') 
                write(19,*) 'ERROR!! Cannot recell the structure!!!'
                close(19)
                stop
        endif

        do i=1,3
           do j=1,3
              privect(i,j)=lattice(i,j)
           enddo
        enddo

C       refine_cell is the atom number in the recelled structure

        nn=refine_cell

        do i=1, ntype
           natomi(i)=0
           do j=1, nn
              if ( atom_types(j) .EQ. i) then
                   natomi(i)=natomi(i)+1
              endif
           enddo
        enddo

        m=1

        do i=1,ntype
           do j=1,nn
              if ( atom_types(j) .EQ. i) then
                   do k=1,3
                      pripos(m,k)=positions(k,j)
                   enddo
                   m=m+1 
              endif              
           enddo
        enddo

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%%%%%          Othorhombic lattice: b < a < c              %%% 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        if (space_group .GE. 16 .AND. space_group .LE. 74) then

                do i=1,2                
                   if ( privect(3,3) .LT. privect(i,i) ) then
                           temp=privect(3,3)
                           privect(3,3)=privect(i,i)
                           privect(i,i)=temp

                           do j=1,nn
                              temp=pripos(j,3)
                              pripos(j,3)=pripos(j,i)
                              pripos(j,i)=temp
                           enddo
                   endif
                enddo

                if ( privect(1,1) .LT. privect(2,2) ) then
                        temp=privect(1,1)
                        privect(1,1)=privect(2,2)
                        privect(2,2)=temp

                        do j=1,nn
                           temp=pripos(j,1)
                           pripos(j,1)=pripos(j,2)
                           pripos(j,2)=temp
                        enddo
                endif               

        endif

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%%%%%              Monoclinic lattice: c <-> b             %%% 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       if (space_group .GE. 3 .AND. space_group .LE. 15) then
                
                temp=privect(2,2)
                privect(2,2)=privect(3,3)
                privect(3,3)=temp

                privect(2,1)=privect(3,1)
                privect(3,1)=0.0000
                     
                do j=1,nn
                   temp=pripos(j,3)
                   pripos(j,3)=pripos(j,2)
                   pripos(j,2)=temp
                enddo   

       endif 
        
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%%%%%              Monoclinic lattice: b < a               %%% 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        if (space_group .GE. 3 .AND. space_group .LE. 15) then        

                lega=sqrt(privect(1,1)**2+privect(1,2)**2
     &                                   +privect(1,3)**2)
                legb=sqrt(privect(2,1)**2+privect(2,2)**2
     &                                   +privect(2,3)**2)

                if ( lega .LT. legb ) then
                       privect(1,1)=legb

                       privect(2,2)=lega*privect(2,2)/legb
                       privect(2,1)=lega*privect(2,1)/legb
                       
                       do j=1,nn
                           temp=pripos(j,1)
                           pripos(j,1)=pripos(j,2)
                           pripos(j,2)=temp
                       enddo
               endif
                
        endif

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%%%%% Create the POSCAR for energy and stress calculation %%%%
C%%%%%%%%% The default output file is "RECELL"                 %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C       output the recelled structure

        if (Ncoord .EQ. 0) THEN
                pricoord='Direct'
        else if (Ncoord .EQ. 1) THEN
                pricoord='Cartesian'
        else
                open(unit=19,position='Append',FILE='RELAS') 
                write(19,*) 'ERROR!! NOT DIRECT OR CARTESIAN!!!'
                close(19)
                stop
        endif

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

        write(19,'(a23)') 'The recelled structure:'
        write(19,*) 

          write(19,'(a4,f8.5)') '    ',alat
          do i=1,3
             write(19,"(3f20.10)")(privect(i,j),j=1,3)
          enddo
          write(19,'(a4,10A5)') '    ',(sname(i),i=1,ntype)
          write(19,'(a4,10I4)') '    ',(natomi(i), i=1,ntype)
          write(19,'(a4,A6)') '    ',pricoord
          do i=1, nn
             write(19,"(3f20.10)") (pripos(i,j),j=1,3)
          enddo

        close(19)

C       end of outputing the recelled structure   


        posmode=00
        outpos='RECELL'

        call writepos(posmode,
     &                outpos,                           
     &                title, 
     &                alat, 
     &                privect, 
     &                sname,
     &                natomi,                
     &                ntype,
     &                nn,
     &                Ndyn,
     &                Ncoord,        
     &                pripos,
     &                fix,
     &                atom_types)   

        end 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                    End of recell2d.f                    %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
