anchor.diff:         type(ps_io_type) :: info
          print *, "call the ps fill up routine"
          call ps_io_type_fill(grid_in,info)
  -       print *, "core charge"
  -       print *, info%rho_core(1)
  +       print *, "information correctly formatted, now outputting"
          call ps_io_save(info)
         endsubroutine siesta_output
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
positive.code:        !set pseudo_potentials
         info%have_psp=.true.
         info%n_psp=nwfs
         allocate(info%psp_l(nwfs))
         info%psp_l=lls(:nwfs)
         allocate(info%psp_j(nwfs))
         do i=1,nwfs,1
         info%psp_j(i)=M_ZERO!not needed for non spin polarised 
         enddo
         allocate(info%psp_v(ndmx,4))
  
         !do ns=1,nbeta
           !lam=lls(ns)
           !do n=1,ikk(ns)
           do lam=0,4
           do n=1,ndmx
             !something
             !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
  
         l=0
         n=1
  
         do while (abs(info%psp_v(n,l)) .gt. 0.00012)
           info%wfs_rc(l)=info%m%r(n)
           n=n+1
         enddo
         print *,  info%wfs_rc(:)
  
  !       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
  
negative.code: 
        end subroutine ps_io_save
    
          subroutine siesta_save(unit, parsec, info)
          !-----------------------------------------------------------------------!
          ! Write the pseudo-potential stuff using the format of the SIESTA code. !
          !-----------------------------------------------------------------------!
          integer, intent(in) :: unit
          logical, intent(in) :: parsec
          type(ps_io_type), intent(in) :: info
      
          integer :: i, l, n_dn, n_up, values(8), n, ir
          real(R8) :: j, occ
          character(3)    :: irel
          character(4)    :: icore
          character(60)   :: header
          character(70)   :: title
          type(mesh_type) :: new_m
          real(R8), allocatable :: dum(:)
      
          !Get spin/relativistic mode and number of states
          select case (info%wave_eq)
          case (DIRAC)
            irel = "rel"
            n_up = 0
            n_dn = 0
            do i = 1, info%n_psp
              if (info%psp_l(i) == 0 .or. (info%psp_l(i) /= 0 .and. &
                   info%psp_j(i) == info%psp_l(i) - M_HALF)) then
                n_dn = n_dn + 1          
              elseif (info%psp_l(i) /= 0 .and. &
                      (info%psp_j(i) == info%psp_l(i) + M_HALF)) then
                n_up = n_up + 1
      
              end if
            end do
          case (SCHRODINGER, SCALAR_REL)
            n_dn = info%n_psp
            select case (info%nspin)
            case (1)
              irel = "nrl"
              n_up = 0
            case (2)
              irel = "isp"
              n_up = n_dn
            end select
          end select
            
          call date_and_time(values=values)
          write(header,'(A,A,3X,I4.4,"/",I2.2,"/",I2.2)') "APE-QE Version-", PACKAGE_VERSION, values(1), values(2), values(3)
          select case (info%scheme)
          case (HAM)
            header = trim(header)//"   Hamann"
          case (TM)
            header = trim(header)//"   Troullier-Martins"
          end select
          
          !Core corrections mode
          if (info%nlcc) then
            icore = "pcec"
          else
            icore = "nc  "
          end if
      
          !
          title = ""
      
          do i=1, info%n_wfs
          print *, info%wfs_n(i),info%wfs_l(i)
          enddo
      
          do l = 0, 3
            n = minval(info%wfs_n, mask=info%wfs_l == l)
            !j = minval(info%wfs_j, mask=(info%wfs_l == l .and. info%wfs_n == n))
      
            do i = 1, info%n_wfs
              !if (info%wfs_n(i) == n .and. info%wfs_l(i) == l .and. info%wfs_j(i) == j) then
              if (info%wfs_n(i) == n .and. info%wfs_l(i) == l) then
                !select case (irel)
                !case('nrl')
                  print *, "Here is what we want"
                  print *, info%wfs_rc
                  write(title,'(A)') (title)
                  write(title,'(A)') trim(title)
                  write(title,'(A,A2,F5.2,"  r=",F5.2,"/")') trim(title), &
                                  info%wfs_label(i), info%wfs_occ(i, 1), info%wfs_rc(i)
       !!          case('isp')
       !!            write(title,'(A,A2,F4.2,1X,F4.2,1X,F4.2,"/")') trim(title), &
       !!                     info%wfs_label(i), info%wfs_occ(i, 1), info%wfs_occ(i, 2), &
       !!                     info%wfs_rc(i)
       !!          case ('rel')
       !!            occ = sum(info%wfs_occ(:, 1), mask=(info%wfs_l == l .and. &
       !!                 (info%wfs_j == l + M_HALF .or. info%wfs_j == l - M_HALF)))
       !!            write(title,'(A,A2,F5.2,"  r=",F5.2,"/")') trim(title), &
       !!                                           info%wfs_label(i), occ, info%wfs_rc(i)
       !!          end select
              end if
            end do
      
          end do
      
          !Mesh
          call mesh_null(new_m)
          call mesh_generation(new_m, LOG2, info%m%r(1), info%m%r(info%m%np), info%m%np) 
          
          !General info
          !write(unit,'(1X,A2,1X,A2,1X,A3,1X,A4)') info%symbol(1:2), ixc_to_icorr(info%ixc), irel, icore! we have a problem with ixc
          write(unit,'(1X,A2,1X,A2,1X,A3,1X,A4)') info%symbol(1:2), info%ixc, irel, icore
          write(unit,'(1X,A60)') header
          write(unit,'(1X,A70)') title
          write(unit,'(1X,2I3,I5,3F20.10)') n_dn, n_up, new_m%np, new_m%b, new_m%a, info%z_val
            
          !Write radial grid
          write(unit,'(" Radial grid follows")')
          write(unit,'(4(g20.12))') (new_m%r(i),i = 1, new_m%np)
          allocate(dum(new_m%np))
      
          !Down pseudopotentials
          do l = 0, 3
            if (info%wave_eq == DIRAC) then
              j = max(l - M_HALF, M_HALF)
            else
              j = M_ZERO
            end if
      
            do i = 1, info%n_psp
              if (info%psp_l(i) == l .and. info%psp_j(i) == j) then
                write(unit,'(" Down Pseudopotential follows (l on next line)")')
                write(unit,'(1X,I2)') info%psp_l(i)
                
                call mesh_transfer(info%m, info%psp_v(:, i), new_m, dum, 3)
                where (abs(dum) < 1.0E-30)
                  dum = M_ZERO
                end where
                write(unit,'(4(g20.12))') (dum(ir)*new_m%r(ir)*M_TWO, ir = 1, new_m%np)
              end if
            end do
          end do
      
          !Up pseudopotentials
          if (n_up /= 0) then
            do l = 0, 3
              do i = 1, info%n_psp
                if (l /= 0 .and. info%wave_eq == DIRAC) then
                  j = l + M_HALF
                else
                  j = M_ZERO
                end if
                
                if (info%psp_l(i) == l .and. info%psp_j(i) == j) then
                  write(unit,'(" Up Pseudopotential follows (l on next line)")')
                  write(unit,'(1X,I2)') info%psp_l(i)
                
                  call mesh_transfer(info%m, info%psp_v(:, i), new_m, dum, 3)
                  where (abs(dum) < 1.0E-30)
                    dum = M_ZERO
                  end where
                  write(unit,'(4(g20.12))') (dum(ir)*new_m%r(ir)*M_TWO, ir = 1, new_m%np)
                end if
              end do
            end do
          end if
      
          !Write core charge
          write(unit,'(" Core charge follows")')
          if (icore == "nc  ") then
            dum = M_ZERO
          else
            call mesh_transfer(info%m, info%rho_core, new_m, dum, 3)
            where (abs(dum) < 1.0E-30)
              dum = M_ZERO
            elsewhere
              dum = new_m%r**2*M_FOUR*M_PI*dum
            end where
          end if
          write(unit,'(4(g20.12))') (dum(i), i = 1, new_m%np)
      
          !Write valence charge
          write(unit,'(" Valence charge follows")')
          call mesh_transfer(info%m, info%rho_val, new_m, dum, 3)
          where (abs(dum) < 1.0E-30)
            dum = M_ZERO
          end where
          write(unit,'(4(g20.12))') (dum(i)*new_m%r(i)**2*M_FOUR*M_PI, i = 1, new_m%np)
      
          if (parsec) then
            !Write pseudo-wave-functions
            do l = 0, 3
              n = minval(info%wfs_n, mask=info%wfs_l == l)
              j = minval(info%wfs_j, mask=(info%wfs_l == l .and. info%wfs_n == n))
      
              do i = 1, info%n_wfs
                if (info%wfs_n(i) == n .and. info%wfs_l(i) == l .and. info%wfs_j(i) == j) then
      
                  call mesh_transfer(info%m, info%wfs(:, i), new_m, dum, 3)
      
                  write(unit,'(1X,A,A2)') 'Pseudo-wave-function follows (l, zelect, rc)  ',info%wfs_label(i)
                  write(unit,'(I2,F6.2,2X,F6.2)') info%wfs_l(i), &
                       info%wfs_occ(i, 1), info%wfs_rc(i)
                  write(unit,'(1P4E19.11)') (dum(ir)*new_m%r(ir), ir = 1, new_m%np)
      
                end if
              end do
      
            end do
          end if
          
          deallocate(dum)
          call mesh_end(new_m)
          
        end subroutine siesta_save
        end module siesta_pass
  