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:         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
  
negative.code:                 dum = M_ZERO
                end where
                !write(unit,'(4(g20.12))') (dum(ir)*new_m%r(ir)*M_TWO, ir = 1, new_m%np)
                !write(unit,'(4(g20.12))') (dum(ir)*new_m%r(ir), ir = 1, new_m%np)
                !why not just use our own data
                write(unit,'(4(g20.12))') (info%psp_v(ir,i)*info%m%r(ir), ir = 1, info%m%np)
                !write(unit,'(4(g20.12))') (dum(ir), ir = 1, new_m%np)!!edit here
              end if
            end do
          end do
   
          !Up pseudopotentials
          if (n_up /= 0) then
            do l = 0, 3
  