write_vasp Subroutine

subroutine write_vasp(self, unit, comment_line)
   class(structure_type), intent(in) :: self
   integer, intent(in) :: unit
   character(len=*), intent(in), optional :: comment_line
   integer :: i, j, izp
   integer, allocatable :: kinds(:), species(:)
   real(wp), allocatable :: inv_lat(:, :)
   real(wp), allocatable :: abc(:, :)

   allocate(species(self%nat))
   allocate(kinds(self%nat), source=1)

   j = 0
   izp = 0
   do i = 1, self%nat
      if (izp.eq.self%id(i)) then
         kinds(j) = kinds(j)+1
      else
         j = j+1
         izp = self%id(i)
         species(j) = self%id(i)
      end if
   end do

   ! use vasp 5.x format
   if (present(comment_line)) then
      write(unit, '(a)') comment_line
   else
      if (allocated(self%comment)) then
         write(unit, '(a)') self%comment
      else
         write(unit, '(a)')
      end if
   end if

   ! scaling factor for lattice parameters is always one
   write(unit, '(f20.14)') self%info%scale
   ! write the lattice parameters
   if (any(self%periodic)) then
      if (size(self%lattice, 2) == 3) then
         write(unit, '(3f20.14)') self%lattice
      else
         write(unit, '(3f20.14)') spread(0.0_wp, 1, 9)
      end if
   else
      write(unit, '(3f20.14)') spread(0.0_wp, 1, 9)
   end if

   do i = 1, j
      write(unit, '(1x, a)', advance='no') self%sym(species(i))
   end do
   write(unit, '(a)')

   ! write the count of the consecutive atom types
   do i = 1, j
      write(unit, '(1x, i0)', advance='no') kinds(i)
   end do
   write(unit, '(a)')
   deallocate(kinds, species)

   if (self%info%selective) write(unit, '("Selective")')

   ! we write cartesian coordinates
   if (any(shape(self%lattice) /= [3, 3]) .or. self%info%cartesian) then
      write(unit, '("Cartesian")')

      ! now write the cartesian coordinates
      do i = 1, self%nat
         write(unit, '(3f20.14)') self%xyz(:, i)*autoaa/self%info%scale
      end do
   else
      write(unit, '("Direct")')
      inv_lat = matinv_3x3(self%lattice)
      abc = matmul(inv_lat, self%xyz)

      ! now write the fractional coordinates
      do i = 1, self%nat
         write(unit, '(3f20.14)') abc(:, i)
      end do
   end if

end subroutine write_vasp