anchor.diff:           end if
          enddo
   
  -       do l=1,nwfx
  +       do l=1,nwfsx!define the proper cutoff values of the wfs
            do n=1,grid_in%mesh
  -           if (info%m%r(n).lt.rcut(l)) then 
  -             info%wfs_rc(l)=info%m%r(n)
  +           if ( info%m%r(n) .lt. rcut(l) ) then
  +             info%wfs_rc(l) = info%m%r(n)
              endif
            enddo
          enddo
   
  -
  -!!       l=1
  -!!       n=1
  -!!       do l =1, nwfx
  -!!         do while (info%m%r(n) .lt. rcut(l) )
  -!!           n=n+1
  -!!           info%wfs_rc(l)=info%m%r(n)
  -!!         enddo
  -!!       enddo
  -!!       print *, "here"
  -!!       print *, "here"
  -!!       print *, ik(:)
  -
  -!!       do while (abs(info%psp_v(n,l)) .gt. 0.00012)
  -!!         info%wfs_rc(l)=info%m%r(n)
  -!!         n=n+1
  -!!       enddo
  -!
  -!!       do l=0,4
  -!!       !do ir=1,ndmx
  -!!       do ir=1,ndmx
  -!!       !!it should always be vnl(:,l,1), i.e. non relativistic
  -!!       !info%psp_v(ir,l)= (vnl(ir,l,2)+vpsloc(ir)*0.5)
  -!!       !info%psp_v(ir,l)= (vpsloc(ir)*0.5)
  -!!       !info%psp_v(ir,l)= (vnl(ir,l,2))
  -!!       info%psp_v(ir,l)= m%r(ir)*(vnl(ir,l,1)+vpsloc(ir))
  -!!       enddo
  -!!       enddo
         endsubroutine ps_io_type_fill
   
         subroutine ps_io_save(info)!this seems only necessary to setup the filename and assign io units
  
positive.code:        info%psp_j(i)=M_ZERO!not needed for non spin polarised 
         enddo
  
         allocate(info%psp_v(ndmx,4))
  
         do lam=1,4
           do n=1,ndmx
  !           !info%psp_v(n,lam)= (vnl(n,lam,1)+m%r(n)*0.5*vpsloc(n))
              info%psp_v(n,lam)= (vnl(n,lam,1)+0.5*vpsloc(n))
           enddo
         enddo
   
         do n=1,nwfsx !set the number of pseudo wavefunctions
           if (ocs(n) .ne. 0) then
  
negative.code:        call ps_io_save(info)
        endsubroutine siesta_output
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        function cutoff_function(r) result(j)
          real(R8), intent(in) :: r ! input
          real(R8):: j ! output
          j = exp(-5*r)!this is the correct siesta way
          if (abs(j) .lt. 1.0d-99) then
            j = 0.0d0
          endif
  
        end function cutoff_function
   
        subroutine ps_io_type_fill(grid_in,info)
         type(radial_grid_type), intent(in) :: grid_in
         type(ps_io_type) ,intent(out) :: info
         !mesh parameters
         type(mesh_type) :: m
         integer :: no_mesh_points,i,j,l,ir,ns,lam,n,jcut
         real(R8) :: rn,r1
         !generalities
         character, external :: atom_name*2
  