main.f90 Source File

Source Code


Source Code

#include "precompilerdefinitions"
program crystal_structure_info
!!{!src/crystal_structure_info/manual.md!}
! Prints general info about the lattice, high symmetry points and so on.
use konstanter, only: flyt,lo_tol,lo_sqtol,lo_bohr_to_A
use helpers, only: lo_frobnorm,lo_mpi_helper,lo_chop,open_file,tochar,lo_clean_fractional_coordinates,lo_determ,&
                   lo_fetch_tolerance,walltime,lo_return_unique,lo_invert3x3matrix,qsort
use options, only: lo_opts
use type_crystalstructure, only: lo_crystalstructure
use type_qpointmesh, only: lo_bandstructure
!
implicit none
! for normal phonon things
type(lo_opts) :: opts
type(lo_crystalstructure) :: uc
type(lo_bandstructure) :: bs
type(lo_mpi_helper) :: mw
!
real(flyt), dimension(3) :: v0,v1
real(flyt) :: tstart,f0 
integer :: i,j,u
character(len=1000) :: opf
!
! Get the command line arguments
!
call mw%init()
tstart=walltime()
call opts%parse()

! Read unitcell
call uc%readfromfile('infile.ucposcar',verbosity=opts%verbosity)
call uc%classify('wedge',timereversal=opts%timereversal)

write(*,*) 'Info about the crystal structure:'
write(*,*) '... I believe the basis is ',trim(uc%info%bravaislatticelongname),' (',trim(uc%info%bravaislattice),') with the basis vectors'
write(*,"(1X,'a1:',3(2X,F18.12))") uc%latticevectors(:,1)*lo_bohr_to_A  
write(*,"(1X,'a2:',3(2X,F18.12))") uc%latticevectors(:,2)*lo_bohr_to_A  
write(*,"(1X,'a3:',3(2X,F18.12))") uc%latticevectors(:,3)*lo_bohr_to_A  
write(*,*) '... and the atoms at these positions (in fractional coordinates)'
do i=1,uc%na
     write(*,'(1X,A8,3(2x,F15.13))') trim(uc%atomic_symbol(uc%species(i))),uc%r(:,i)
enddo
! Print what the conventional basis is
f0=lo_frobnorm(uc%latticevectors-uc%info%unique_conventional_basis)
if ( f0 .gt. lo_tol ) then
    write(*,*) 'I believe this is the conventional lattice:'
    write(*,"(1X,'a1:',3(2X,F18.12))") uc%info%unique_conventional_basis(:,1)*lo_bohr_to_A 
    write(*,"(1X,'a2:',3(2X,F18.12))") uc%info%unique_conventional_basis(:,2)*lo_bohr_to_A 
    write(*,"(1X,'a3:',3(2X,F18.12))") uc%info%unique_conventional_basis(:,3)*lo_bohr_to_A 
endif

!! It might not agree with my canonical form.
!f0=lo_frobnorm(uc%latticevectors-uc%info%unique_primitive_basis)
!if ( f0 .gt. lo_tol ) then
!    write(*,*) ' '
!    write(*,*) 'I like this primitive basis better:'
!    write(*,"(1X,'a1:',3(2X,F18.12))") uc%info%unique_primitive_basis(:,1)
!    write(*,"(1X,'a2:',3(2X,F18.12))") uc%info%unique_primitive_basis(:,2)
!    write(*,"(1X,'a3:',3(2X,F18.12))") uc%info%unique_primitive_basis(:,3)
!    write(*,*) '... with the atoms at these positions (in fractional coordinates)'
!    ! Convert the coordinates    
!    allocate(dum(3,uc%na))
!    do i=1,uc%na
!        ! get the current coordinate
!        v0=uc%r(:,i)
!        v0=matmul(uc%info%permutation_to_unique,v0)
!        v0=lo_chop(lo_clean_fractional_coordinates(v0),lo_sqtol)
!        dum(:,i)=v0
!        write(*,'(1X,A8,3(2x,F15.13))') trim(uc%atomic_symbol(uc%species(i))),dum(:,i)
!    enddo
!endif

! Get the high symmetry points
if ( uc%info%havewedge .eqv. .false. ) call uc%classify('wedge',timereversal=opts%timereversal)

write(*,*) ' '
write(*,*) 'Available high symmetry points:'
write(*,*) '(the ones called NP are in fact high symmetry points, but noone has bothered to name them)'
write(*,*) '      Cartesian coordinates                        Fractional coordinates'
do i=1,uc%irrw%nnodes
    v0=lo_chop(uc%irrw%r(:,i),lo_tol)
    v1=lo_chop(matmul(uc%inv_reciprocal_latticevectors,v0),lo_tol)
    write(*,"(1X,I2,1X,A4,3(1X,F13.10),3X,3(1X,F13.10))") i,uc%irrw%label(i),v0/lo_bohr_to_A,v1
enddo

call uc%write_bz_to_hdf5('outfile.brillouin_zone.hdf5')
call bs%standardpath(uc)
bs%npts=100
write(*,*) ' '
write(*,*) 'Per default, the following path will be used:'
do i=1,bs%npath
    write(*,*) trim(bs%symb_q_start(i)),' -> ',trim(bs%symb_q_end(i))
enddo
write(*,*) 'It is written in "outfile.qpoints_dispersion", modify and'
write(*,*) 'copy it to "infile.qpoints_dispersion" if you want'

! Dump a path file
u=open_file('out','outfile.qpoints_dispersion')
    opf="(A5,22X,' ! Bravais lattice type')"
    write(u,opf) uc%info%bravaislattice
    opf="(I5,22X,' ! Number of points on each path')"
    write(u,opf) 100
    opf="(I5,22X,' ! Number paths between special points')"
    write(u,opf) bs%npath
    opf="(A3,1X,A3,20X,' ! Starting and ending special point')"
    write(u,opf) bs%symb_q_start(1),bs%symb_q_end(1)
    opf="(A3,1X,A3,20X,' !')"
    do i=2,bs%npath
        write(u,opf) bs%symb_q_start(i),bs%symb_q_end(i)
    enddo
close(u)
 
if ( opts%printsymmetry ) then
    write(*,*) ' Symmetry operations:'
    do i=1,uc%sym%n
        write(*,*) 'Operation ',tochar(i),' ',uc%sym%op(i)%whatkind
        write(*,*) '        Axis:',lo_chop(uc%sym%op(i)%axis,lo_sqtol)
        if ( uc%sym%op(i)%fold .gt. 0 ) write(*,*) '       Angle: ',tochar(int(anint(360.0_flyt/uc%sym%op(i)%fold)))
        opf='        FMAP:'
        do j=1,size(uc%sym%op(i)%fmap)
            opf=trim(opf)//' '//tochar(j)//'->'//tochar(uc%sym%op(i)%fmap(j))
        enddo
        write(*,*) trim(opf)
        !write(*,*) '        FMAP:',uc%sym%op(i)%fmap
        write(*,'(1X,A,3(1X,F9.6))') ' translation:',uc%sym%op(i)%ftr
        write(*,*) ' ... in Cartesian coordinates'
        do j=1,3
            write(*,"(6(3X,F15.11))") uc%sym%op(i)%m(j,:)
            !write(*,"(3(3X,A10))") lo_fancy_real2char(uc%sym%op(i)%m(j,1)),lo_fancy_real2char(uc%sym%op(i)%m(j,2)),lo_fancy_real2char(uc%sym%op(i)%m(j,3))
        enddo
        write(*,*) ' ... in fractional coordinates'
        do j=1,3
            write(*,"(6(3X,F15.11))") uc%sym%op(i)%fm(j,:)
            !write(*,"(6(3X,A3))") lo_fancy_real2char(uc%sym%op(i)%fm(j,1)),lo_fancy_real2char(uc%sym%op(i)%fm(j,2)),lo_fancy_real2char(uc%sym%op(i)%fm(j,3))
        enddo
        write(*,*) ' '
    enddo

    ! Dump symmetry operations to Mathematica-format, in case you need to play with it
    u=open_file('out','outfile.symops_mathematica')
        write(u,*) 'ops=ConstantArray[0,{'//tochar(uc%sym%n)//',3,3}];'
        do i=1,uc%sym%n
            do j=1,3
                write(u,*) 'ops[['//tochar(i)//','//tochar(j)//&
                ']]={ '//lo_fancy_real2char(uc%sym%op(i)%m(1,j)) //','//lo_fancy_real2char(uc%sym%op(i)%m(2,j))//','//lo_fancy_real2char(uc%sym%op(i)%m(3,j))//'};'
            enddo
        enddo
    close(u)
endif

write(*,*) 'all done'
write(*,*) ' '
call mw%destroy()

contains

!!> Find the primitive cell
!subroutine find_primitive_cell(basis,r,species,flavor,primbasis,tolerance)
!    !> cell to be reduced
!    real(flyt), dimension(3,3), intent(in) :: basis
!    !> atoms in the cell
!    real(flyt), dimension(:,:), intent(in) :: r
!    !> species classification
!    integer, dimension(:), intent(in) :: species
!    !> flavor classification
!    integer, dimension(:), intent(in) :: flavor
!    !> primitive lattice vectors
!    real(flyt), dimension(3,3), intent(out) :: primbasis
!    !> tolerance
!    real(flyt), intent(in) :: tolerance
!    !
!    real(flyt), dimension(:,:), allocatable :: r0,r1
!    real(flyt), dimension(3,3) :: basis0,basis1
!    integer, dimension(:), allocatable :: species0,flavor0,species1,flavor1
!    integer :: i,na
!    !
!    na=size(r,2)
!    allocate(r0(3,na))
!    allocate(species0(na))
!    allocate(flavor0(na))
!    basis0=basis
!    r0=r
!    species0=species
!    flavor0=flavor
!
!    ! Recursively find a smaller cell
!    redloop: do
!        ! Try to find a smaller cell
!        call reduce_cell(basis0,r0,species0,flavor0,basis1,r1,species1,flavor1,tolerance)
!        if ( abs(abs(lo_determ(basis0))-abs(lo_determ(basis1))) .lt. tolerance ) then
!            ! this means no new cell, we are converged.
!            exit redloop
!        else
!            ! this means a new cell, update and try again!
!            if ( allocated(r0) )       deallocate(r0)
!            if ( allocated(species0) ) deallocate(species0)
!            if ( allocated(flavor0) )  deallocate(flavor0)
!            i=size(r1,2)
!            allocate(r0(3,i))
!            allocate(species0(i))
!            allocate(flavor0(i))
!            basis0=basis1
!            r0=r1
!            species0=species1
!            flavor0=flavor1
!            if ( allocated(r1) )       deallocate(r1)
!            if ( allocated(species1) ) deallocate(species1)
!            if ( allocated(flavor1) )  deallocate(flavor1)
!        endif
!    enddo redloop
!
!    ! Here we should be done.
!    primbasis=basis0
!
!    if ( allocated(r0) )       deallocate(r0)
!    if ( allocated(r1) )       deallocate(r1)
!    if ( allocated(species0) ) deallocate(species0)
!    if ( allocated(species1) ) deallocate(species1)
!    if ( allocated(flavor0) )  deallocate(flavor0)
!    if ( allocated(flavor1) )  deallocate(flavor1)
!end subroutine

!    !> single iteration of the cell reduction
!    subroutine reduce_cell(latticevectors,r,species,flavor,newlatticevectors,newr,newspecies,newflavor,tolerance)
!        !> original lattice vectors
!        real(flyt), dimension(3,3), intent(in) :: latticevectors
!        !> original positions
!        real(flyt), dimension(:,:), intent(in) :: r
!        !> original species
!        integer, dimension(:), intent(in) :: species
!        !> additional flavor for symmetry thing
!        integer, dimension(:), intent(in) :: flavor
!        !> tolerance
!        real(flyt), intent(in) :: tolerance
!
!        real(flyt), dimension(3,3), intent(out) :: newlatticevectors
!        real(flyt), dimension(:,:), allocatable, intent(out) :: newr
!        integer, dimension(:), allocatable, intent(out) :: newspecies
!        integer, dimension(:), allocatable,intent(out) :: newflavor
!
!        real(flyt), dimension(:,:), allocatable :: vecs,dum
!        real(flyt), dimension(3,3) :: m0,m1
!        real(flyt) :: vol,f0,tol,reltol,frtol
!        integer :: i,j,k,l,nvecs,na
!        logical :: cellvalid
!
!        ! set the tolerance
!        call lo_fetch_tolerance(tolerance,latticevectors,realspace_cart_tol=tol,realspace_fractional_tol=frtol,relative_tol=reltol)
!
!        na=size(r,2)
!        ! Get a list of possible lattice vectors: All the vectors in the cell + old latticevectors
!        l=na+3
!        allocate(dum(3,na+3))
!        dum(:,1:na)=lo_clean_fractional_coordinates(r,frtol**2)
!        do i=1,na
!            dum(:,i)=matmul(latticevectors,dum(:,i))
!        enddo
!        do i=1,3
!            dum(:,i+na)=latticevectors(:,i)
!        enddo
!        ! Reduce these to the unique. Should not be necessary, but never hurts.
!        call lo_return_unique(dum,vecs,frtol)
!        nvecs=size(vecs,2)
!
!        ! volume of original cell
!        vol=abs(lo_determ(latticevectors))
!        ! Loop over all vectors to try to find a new cell.
!        vl1: do i=1,nvecs
!        vl2: do j=i+1,nvecs
!        vl3: do k=j+1,nvecs
!            ! possible new set of lattice vectors
!            m0(:,1)=vecs(:,i)
!            m0(:,2)=vecs(:,j)
!            m0(:,3)=vecs(:,k)
!            ! some quick tests to see if there is any point in continuing:
!            f0=abs(lo_determ(m0))
!            if ( f0 .lt. tol ) cycle vl3                           ! stupid cell if volume is zero
!            if ( abs(mod(vol/f0,1.0_flyt)) .gt. reltol ) cycle vl3 ! has to be an integer division
!            if ( int(anint(vol/f0)) .eq. 1 ) cycle vl3             ! has to make the cell smaller, not the same
!            ! I can get really retarded cells from this, do the Delaunay thing on it
!            call reduce_basis_to_something_decent(m0,m1)
!            ! Now do an actual test
!            call test_cell(latticevectors,r,species,flavor,m1,cellvalid,newr,newspecies,newflavor,tolerance)
!            if ( cellvalid ) then
!                ! If I find a new cell, we will exit this routine and start over. Makes things way way faster.
!                newlatticevectors=m1
!                deallocate(vecs)
!                return
!            endif
!        enddo vl3
!        enddo vl2
!        enddo vl1
!        ! If we made it here, the cell is already the smallest. Return the same basis again.
!        newlatticevectors=latticevectors
!    end subroutine
!
!    !> test if a new cell is a smaller unit cell
!    subroutine test_cell(oldcell,oldr,oldspecies,oldflavor,cell,valid,newr,newspecies,newflavor,tolerance)
!        !> original cell
!        real(flyt), dimension(3,3), intent(in) :: oldcell
!        !> original positions
!        real(flyt), dimension(:,:), intent(in) :: oldr
!        !> original species
!        integer, dimension(:), intent(in) :: oldspecies
!        !> original flavor
!        integer, dimension(:), intent(in) :: oldflavor
!        !> new cell
!        real(flyt), dimension(3,3), intent(in) :: cell
!        !> is it a new cell?
!        logical, intent(out) :: valid
!        !> new positions in case of a valid cell
!        real(flyt), dimension(:,:), allocatable, intent(out) :: newr
!        !> new species in case of a valid cell
!        integer, dimension(:), allocatable, intent(out) :: newspecies
!        !> new flavor in case of a valid cell
!        integer, dimension(:), allocatable, intent(out) :: newflavor
!        !> tolerance
!        real(flyt), intent(in) :: tolerance
!        !
!        real(flyt), dimension(3,3) :: icell,m0
!        real(flyt), dimension(:,:), allocatable :: dumr1,dumr2
!        real(flyt) :: tol,sqtol
!        integer :: i,j,l,mult,oldna,newna
!
!        ! set the tolerance
!        call lo_fetch_tolerance(tolerance,oldcell,realspace_fractional_tol=tol)
!        sqtol=tol**2
!
!        valid=.false.
!        ! How many atoms should I get of each?
!        mult=int(anint(abs(lo_determ(oldcell)/lo_determ(cell))))
!
!        ! m0 is a matrix that converts to fractional coordinates in the new cell
!        icell=lo_invert3x3matrix( cell )
!        m0=matmul(icell,oldcell)
!
!        ! convert coordinates to fractional with respect to the new cell. Also
!        ! I decorate the positions with the species and flavor, to keep track of that.
!        oldna=size(oldr,2)
!        allocate(dumr1(5,oldna))
!        do i=1,oldna
!            dumr1(1:3,i)=lo_clean_fractional_coordinates(matmul(m0,oldr(:,i)),sqtol)
!            dumr1(4:5,i)=[oldspecies(i),oldflavor(i)]
!        enddo
!        ! Reduce this to the union of the list, the unique
!        call lo_return_unique(dumr1,dumr2,tol)
!
!        ! Did it become the correct number of atoms?
!        if ( abs(size(dumr2,2)*mult-oldna) .ne. 0 ) then
!            valid=.false.
!            return
!        endif
!
!        ! Slightly better check: Each of the atoms in the reduced cell must appear
!        ! exactly mult times in the old cell.
!        do i=1,size(dumr2,2)
!            l=0
!            do j=1,oldna
!                if ( sum(abs(dumr1(:,j)-dumr2(:,i))) .lt. tol ) l=l+1
!            enddo
!            ! kill if not true
!            if ( l .ne. mult ) then
!                valid=.false.
!                return
!            endif
!        enddo
!
!        ! ok, now I say this is valid, return the new cell
!        newna=size(dumr2,2)
!        allocate(newr(3,newna))
!        allocate(newspecies(newna))
!        allocate(newflavor(newna))
!        newr=dumr2(1:3,:)
!        newspecies=int(anint(dumr2(4,:)))
!        newflavor=int(anint(dumr2(5,:)))
!        valid=.true.
!    end subroutine

!!> reduce a set of lattice vectors to canonical form
!subroutine reduce_basis_to_something_decent(basis,canonical_basis)
!    !> original lattice vectors
!    real(flyt), dimension(3,3), intent(in) :: basis
!    !> canonical lattice vectors
!    real(flyt), dimension(3,3), intent(out) :: canonical_basis
!
!    real(flyt), dimension(3,4) :: ext_basis
!    real(flyt), dimension(3,7) :: dum
!    real(flyt), dimension(7) :: d
!    real(flyt) :: f0
!    integer, dimension(7) :: ind
!    integer :: i,j,k
!
!    ! First do the Delaunay reduction thing
!    ext_basis=0.0_flyt
!    ext_basis(:,1:3)=basis
!    do i=1,3
!        ext_basis(:,4)=ext_basis(:,4)-basis(:,i)
!    enddo
!    ! Iterate and stuff
!    iterloop: do
!        do i=1,4
!        do j=i+1,4
!            ! check scalar product
!            f0=dot_product(ext_basis(:,i),ext_basis(:,j))
!            if ( f0 .gt. lo_sqtol ) then
!                ! subtract i from all not i and j
!                do k=1,4
!                    if ( k .ne. i .and. k .ne. j ) then
!                        ext_basis(:,k)=ext_basis(:,k)+ext_basis(:,i)
!                    endif
!                enddo
!                ! flip sign of i
!                ext_basis(:,i)=-ext_basis(:,i)
!                ! start over!
!                cycle iterloop
!            endif
!        enddo
!        enddo
!        ! If we make it here, we are done!
!        exit iterloop
!    enddo iterloop
!
!    ! Look through all variants to get shortest vectors
!    dum=0.0_flyt
!    dum(:,1:4)=ext_basis
!    dum(:,5)=ext_basis(:,1)+ext_basis(:,2)
!    dum(:,6)=ext_basis(:,2)+ext_basis(:,3)
!    dum(:,7)=ext_basis(:,1)+ext_basis(:,3)
!    do i=1,7
!        d(i)=norm2(dum(:,i))
!    enddo
!    call qsort(d,ind)
!    !
!    do i=3,7
!        canonical_basis(:,1)=dum(:,ind(1))
!        canonical_basis(:,2)=dum(:,ind(2))
!        canonical_basis(:,3)=dum(:,ind(i))
!        if ( abs(lo_determ(canonical_basis)) .gt. lo_tol ) exit
!    enddo
!    if ( lo_determ(canonical_basis) .lt. 0.0_flyt ) canonical_basis=-canonical_basis
!end subroutine

! !> Check if an array has repeated values
! logical function lo_is_vector_degenerate_brute(v,thres)
!     !> array to test
!     real(flyt), dimension(:), intent(in) :: v
!     !> tolerance
!     real(flyt), intent(in) :: thres
!     !
!     integer :: i,j
!     real(flyt) :: f0,f1
!     !
!     f0=lo_huge
!     do i=1,size(v,1)
!         do j=i+1,size(v,1)
!             f1=abs(v(i)-v(j))
!             if ( f1 .lt. f0 ) f0=f1
!         enddo
!     enddo
!     !
!     if ( f0 .lt. thres ) then
!         lo_is_vector_degenerate_brute=.true.
!     else
!         lo_is_vector_degenerate_brute=.false.
!     endif
! end function
! 
! !> Check if an array has repeated values
! logical function lo_is_vector_degenerate_qsort(v,thres)
!     !> array to test
!     real(flyt), dimension(:), intent(in) :: v
!     !> tolerance
!     real(flyt), intent(in) :: thres
!     !
!     integer :: i,j
!     real(flyt) :: f0,f1
!     real(flyt), dimension(:), allocatable :: dum
!     !
!     lo_allocate(dum(size(v,1)))
!     dum=v
!     call qsort(dum)
!     !
!     lo_is_vector_degenerate_qsort=.false.
!     do i=1,size(dum,1)-1
!         if ( abs(dum(i)-dum(i+1)) .lt. thres ) then
!             lo_is_vector_degenerate_qsort=.true.
!             return
!         endif
!     enddo
! end function

!> converts a double to a character, and tries to find well defined numbers
function lo_fancy_real2char(x,tolerance) result(s)
    !> number to be converted
    real(flyt), intent(in) :: x
    !> optional tolerance
    real(flyt), intent(in), optional :: tolerance
    !> the string
    character(len=:), allocatable :: s

    real(flyt) :: tol
    real(flyt), dimension(67), parameter :: welldefined=[&
          0.0000000000000000_flyt,&
          0.1000000000000000_flyt,&
         -0.1000000000000000_flyt,&
          0.1111111111111111_flyt,&
         -0.1111111111111111_flyt,&
          0.1250000000000000_flyt,&
         -0.1250000000000000_flyt,&
          0.1428571428571428_flyt,&
         -0.1428571428571428_flyt,&
          0.1666666666666667_flyt,&
         -0.1666666666666667_flyt,&
          0.2000000000000000_flyt,&
         -0.2000000000000000_flyt,&
          0.2222222222222222_flyt,&
         -0.2222222222222222_flyt,&
          0.2500000000000000_flyt,&
         -0.2500000000000000_flyt,&
          0.2857142857142857_flyt,&
         -0.2857142857142857_flyt,&
          0.3000000000000000_flyt,&
         -0.3000000000000000_flyt,&
          0.3333333333333333_flyt,&
         -0.3333333333333333_flyt,&
          0.3750000000000000_flyt,&
         -0.3750000000000000_flyt,&
          0.4000000000000000_flyt,&
         -0.4000000000000000_flyt,&
          0.4285714285714285_flyt,&
         -0.4285714285714285_flyt,&
          0.4444444444444444_flyt,&
         -0.4444444444444444_flyt,&
          0.5000000000000000_flyt,&
         -0.5000000000000000_flyt,&
          0.5555555555555556_flyt,&
         -0.5555555555555556_flyt,&
          0.5714285714285714_flyt,&
         -0.5714285714285714_flyt,&
          0.6000000000000000_flyt,&
         -0.6000000000000000_flyt,&
          0.6250000000000000_flyt,&
         -0.6250000000000000_flyt,&
          0.6666666666666667_flyt,&
         -0.6666666666666667_flyt,&
          0.7000000000000000_flyt,&
         -0.7000000000000000_flyt,&
          0.7142857142857143_flyt,&
         -0.7142857142857143_flyt,&
          0.7500000000000000_flyt,&
         -0.7500000000000000_flyt,&
          0.7777777777777778_flyt,&
         -0.7777777777777778_flyt,&
          0.8000000000000000_flyt,&
         -0.8000000000000000_flyt,&
          0.8333333333333334_flyt,&
         -0.8333333333333334_flyt,&
          0.8571428571428571_flyt,&
         -0.8571428571428571_flyt,&
          0.8660254037844386_flyt,&
         -0.8660254037844386_flyt,&
          0.8750000000000000_flyt,&
         -0.8750000000000000_flyt,&
          0.8888888888888888_flyt,&
         -0.8888888888888888_flyt,&
          0.9000000000000000_flyt,&
         -0.9000000000000000_flyt,&
          1.0000000000000000_flyt,&
         -1.0000000000000000_flyt]
    character(len=12), dimension(67), parameter :: welldefined_string=[&
      '0          ',& !    0.0000000000000000_flyt,&
      '1/10       ',& !    0.1000000000000000_flyt,&
      '-1/10      ',& !   -0.1000000000000000_flyt,&
      'x          ',& !    0.1111111111111111_flyt,&
      'x          ',& !   -0.1111111111111111_flyt,&
      '1/8        ',& !    0.1250000000000000_flyt,&
      '/1/8       ',& !   -0.1250000000000000_flyt,&
      'x          ',& !    0.1428571428571428_flyt,&
      'x          ',& !   -0.1428571428571428_flyt,&
      '1/6        ',& !    0.1666666666666667_flyt,&
      '-1/6       ',& !   -0.1666666666666667_flyt,&
      '1/5        ',& !    0.2000000000000000_flyt,&
      '-1/5       ',& !   -0.2000000000000000_flyt,&
      'x          ',& !    0.2222222222222222_flyt,&
      'x          ',& !   -0.2222222222222222_flyt,&
      '1/4        ',& !    0.2500000000000000_flyt,&
      '-1/4       ',& !   -0.2500000000000000_flyt,&
      'x          ',& !    0.2857142857142857_flyt,&
      'x          ',& !   -0.2857142857142857_flyt,&
      'x          ',& !    0.3000000000000000_flyt,&
      'x          ',& !   -0.3000000000000000_flyt,&
      'x          ',& !    0.3333333333333333_flyt,&
      'x          ',& !   -0.3333333333333333_flyt,&
      'x          ',& !    0.3750000000000000_flyt,&
      'x          ',& !   -0.3750000000000000_flyt,&
      'x          ',& !    0.4000000000000000_flyt,&
      'x          ',& !   -0.4000000000000000_flyt,&
      'x          ',& !    0.4285714285714285_flyt,&
      'x          ',& !   -0.4285714285714285_flyt,&
      'x          ',& !    0.4444444444444444_flyt,&
      'x          ',& !   -0.4444444444444444_flyt,&
      '1/2        ',& !    0.5000000000000000_flyt,&
      '-1/2       ',& !   -0.5000000000000000_flyt,&
      'x          ',& !    0.5555555555555556_flyt,&
      'x          ',& !   -0.5555555555555556_flyt,&
      'x          ',& !    0.5714285714285714_flyt,&
      'x          ',& !   -0.5714285714285714_flyt,&
      'x          ',& !    0.6000000000000000_flyt,&
      'x          ',& !   -0.6000000000000000_flyt,&
      'x          ',& !    0.6250000000000000_flyt,&
      'x          ',& !   -0.6250000000000000_flyt,&
      'x          ',& !    0.6666666666666667_flyt,&
      'x          ',& !   -0.6666666666666667_flyt,&
      'x          ',& !    0.7000000000000000_flyt,&
      'x          ',& !   -0.7000000000000000_flyt,&
      'x          ',& !    0.7142857142857143_flyt,&
      'x          ',& !   -0.7142857142857143_flyt,&
      'x          ',& !    0.7500000000000000_flyt,&
      'x          ',& !   -0.7500000000000000_flyt,&
      'x          ',& !    0.7777777777777778_flyt,&
      'x          ',& !   -0.7777777777777778_flyt,&
      'x          ',& !    0.8000000000000000_flyt,&
      'x          ',& !   -0.8000000000000000_flyt,&
      'x          ',& !    0.8333333333333334_flyt,&
      'x          ',& !   -0.8333333333333334_flyt,&
      'x          ',& !    0.8571428571428571_flyt,&
      'x          ',& !   -0.8571428571428571_flyt,&
      'Sqrt[3]/2  ',& !    0.8660254037844386_flyt,&
      '-Sqrt[3]/2 ',& !   -0.8660254037844386_flyt,&
      'x          ',& !    0.8750000000000000_flyt,&
      'x          ',& !   -0.8750000000000000_flyt,&
      'x          ',& !    0.8888888888888888_flyt,&
      'x          ',& !   -0.8888888888888888_flyt,&
      '9/10       ',& !    0.9000000000000000_flyt,&
      '-9/10      ',& !   -0.9000000000000000_flyt,&
      '1          ',& !    1.0000000000000000_flyt,&
      '-1         '] !   -1.0000000000000000_flyt]
    integer :: i

    if ( present(tolerance) ) then
        tol=tolerance
    else
        tol=lo_sqtol
    endif
    
    do i=1,67
        if ( abs(x-welldefined(i)) .lt. tol ) then
            s=trim(welldefined_string(i))
            return
        endif
    enddo
    s='dunno'

end function

end program