C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C-->   subroutine "kpta" is used to create the KPOINTS file    %%%%
C-->   according the "K-Point Per Reciprocal Atom"             %%%%
C-->                        REFERENCE                          %%%%
C-->   1\ Walle A V D, Ceder G. Automating first-principles    %%%%
C-->      phase diagram calculations[J]. Journal of Phase      %%%%
C-->      Equilibria, 2002, 23(4):348-359.                     %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        subroutine kpta(kppra, kscheme)

        integer kppra
        character*1 kscheme
                
        double precision alat,pripos,privect
        dimension privect(3,3),pripos(1000,3)
        character*10 title,fix
        character*1 pricoord,seldyn
        character*5 sname
        character*10 inpos   
        dimension sname(10),fix(1000,3)
        integer i,j,k,ntype,natomi,nn
        integer Ncoord,Ndyn,atom_types,anum
        dimension natomi(10),atom_types(1000),anum(11)
        
        integer imesh
        real, parameter:: PI=3.1415926       
        double precision bestrf, zero_tolerance
        double precision proj, fmesh, dotnorm, crossnorm
        double precision recipvect, tempvect, crossvect
        double precision normalizer, nb_kpts

        dimension proj(3), imesh(3), fmesh(3), IS(3), JS(3)
        dimension dotnorm(3), crossnorm(3)
        dimension crossvect(3,3), recipvect(5,5), tempvect(3,3)

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%  1) KPPRA stands for "K-Point Per Reciprocal Atom".      %%%%
C%%%%     This is a way to automatically set the k-point       %%%%
C%%%%     mesh for a number of similar systems. Here is        %%%%
C%%%%     how it works. If you know that you need 500          %%%%
C%%%%     k-points for a 2-atom structure, then you type       %%%% 
C%%%%     in 1000. If you now try a similar structure with     %%%%
C%%%%     4 atoms, the code automatically use 250 k-points.    %%%%
C%%%%     The mesh along the three axes is automatically       %%%%
C%%%%     chosen to make the mesh as uniform as possible.      %%%%
C%%%%  2) kscheme specifies the type mesh to use: Gamma-       %%%%
C%%%%     centered or Monkorst-Pack. ("m" "M" or "g" "G")      %%%% 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        inpos='RECELL'

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

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

        write(19,*) 'After reading, the three lattice vectors:'
        do i=1,3
           write(19,'(3f20.15)') (privect(i,j),j=1,3)
        enddo
        write(19,'(10A5)') (sname(i),i=1,ntype)
        write(19,'(10I4)') (natomi(i), i=1,ntype)

        write(19,*) '  -----------------------------------------'

        write(19,'(a32,i5)') ' K-Point Per Reciprocal Atom is ', kppra
        write(19,*) '  -----------------------------------------'

        close(19)

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%         define the reciprocal matrix of privect          %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        do i=1,3
           do j=1,3
              tempvect(i,j)=privect(i,j)
           enddo
        enddo

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

        do i=1,3
           do j=1,3
              recipvect(i,j)=2*PI*tempvect(j,i)
           enddo
        enddo

        do i=4,5
           do j=1,3
              recipvect(i,j)=recipvect(i-3,j)
           enddo
        enddo

        do i=1,5
           do j=4,5
              recipvect(i,j)=recipvect(i,j-3)
           enddo
        enddo

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%                     determin the proj(i)                 %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        do i=1,3
           do j=1,3
              crossvect(i,j)=recipvect(i+1,j+1)*recipvect(i+2,j+2)
     &                        -recipvect(i+1,j+2)*recipvect(i+2,j+1)
           enddo
        enddo

        do i=1,3
           dotnorm(i)=0.000
           do j=1,3
           dotnorm(i)=dotnorm(i)+recipvect(i,j)*crossvect(i,j)
           enddo
           dotnorm(i)=ABS(dotnorm(i))
        enddo

        do i=1,3
           crossnorm(i)=0.0000
           do j=1,3
              crossnorm(i)=crossnorm(i)+crossvect(i,j)*crossvect(i,j)
           enddo
           crossnorm(i)=SQRT(crossnorm(i))
        enddo

        do i=1,3
           proj(i)=dotnorm(i)/crossnorm(i)
        enddo

        nb_kpts=kppra/nn
        normalizer=(proj(1)*proj(2)*proj(3)/nb_kpts)**(1.00/3.00)
        do i=1,3
           imesh(i)=NINT(proj(i)/normalizer)
           fmesh(i)=proj(i)/normalizer-REAL(imesh(i))
        enddo

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%                     Floor and Ceiling                    %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if(kscheme.NE.'m' .AND. kscheme.NE.'g' 
     &     .AND. kscheme.NE.'M' .AND. kscheme.NE.'G') then
                open(unit=19,position='Append',FILE='RELAS')
                write(19,*) 'Error!!! NOT "M" OR "G" FOR KSCHEME'
                close(19)
                stop
        endif

        if(kscheme.EQ.'m' .OR. kscheme.EQ.'M') then 
                do i=1,3
                   if( MOD(imesh(i),2) .EQ. 1 ) then
                           imesh(i)=imesh(i)+1
                           fmesh(i)=0.0000
                   endif
                enddo
        endif

        zero_tolerance=1E-4

        DO 200, WHILE ((imesh(1)*imesh(2)*imesh(3)) .LT. nb_kpts)
        
        bestrf=0.0000

        do i=1,3
           if (fmesh(i) .GT. bestrf) then
                   bestrf=fmesh(i)
           endif
        enddo

        do i=1,3
           if (ABS(fmesh(i)-bestrf) .LT. zero_tolerance) then
                   imesh(i)=imesh(i)+1
                   if(kscheme.EQ.'m' .OR. kscheme.EQ.'M') then
                        imesh(i)=imesh(i)+1
                   endif
                   fmesh(i)=0.0000
           endif
        enddo
200     CONTINUE

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%                 writing the KPOINTS file                 %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        open(unit=19,position='Append',FILE='RELAS') 

        write(19,'(A14)')  'Automatic mesh'
        write(19,'(A3)') '  0'
        if(kscheme.EQ.'m' .OR. kscheme.EQ.'M') then
                write(19,'(A14)')  'Monkhorst-pack'
        else
                write(19,'(A14)')   'Gamma-centered'
        endif
        write(19,'(3I3)') (imesh(i), i=1,3)
        write(19,'(A9)') '  0  0  0'

        close(19)


        open(unit = 30, file = 'NEWKPT', status = 'unknown')

        write(30,'(A14)')  'Automatic mesh'
        write(30,'(A3)') '  0'
        if(kscheme.EQ.'m' .OR. kscheme.EQ.'M') then
                write(30,'(A14)')  'Monkhorst-pack'
        else
                write(30,'(A14)')   'Gamma-centered'
        endif
        write(30,'(3I3)') (imesh(i), i=1,3)
        write(30,'(A9)') '  0  0  0'

        close(30)

        end 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%                       End of kpta.f                      %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
