anchor.diff:          !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
  +!        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)
  +        where (abs(info%rho_val) < 1.0E-30)
  +          info%rho_val = 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)
  +        !write(unit,'(4(g20.12))') (info%rho_val(i)*info%m%r(i)**2*M_FOUR*M_PI, i = 1, info%m%np)
  +        write(unit,'(4(g20.12))') (info%rho_val(i), i = 1, info%m%np)
       
           if (parsec) then
             !Write pseudo-wave-functions
  
positive.code:         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
  
negative.code:        info%psp_j(i)=M_ZERO!not needed for non spin polarised 
         enddo
  
  
  !set the potential
         allocate(info%psp_v(ndmx,4))
  
         do lam=0,3
           do n=1,m%np
             info%psp_v(n,lam+1)= (vnl(n,lam,1)+vpsloc(n))!correct
           enddo
         enddo
  
  
  
         !get jcut mesh index. this is the point where we should just cut off the calculated potential values 
         do lam=1,4
           do n=1,m%np
             if (abs(info%psp_v(n,lam)+2*zval/m%r(n)) .gt. 1.d-3) then
             !if (abs(info%psp_v(n,lam)) .gt. 1.d-3) then
               jcut=n
             endif
           enddo
         !put jcut in the middle of the mesh
  
         print *, "jcut values"
         print *, jcut
         print *, info%m%r(jcut)
  
         !smooth the curve!this seems all correct 4/8/11
           do n=jcut,m%np
             info%psp_v(n,lam)=-2*zval/m%r(n) + cutoff_function(m%r(n)-m%r(jcut))*(info%psp_v(n,lam)+2*zval/m%r(n))
           enddo
         enddo!end the lam loop
  
  