anchor.diff:         info%psp_j(i)=M_ZERO!not needed for non spin polarised 
          enddo
   
  +
   !set the potential
          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)+vpsloc(n))
  -           !info%psp_v(n,lam)= (vnl(n,lam,1)+0.5*vpsloc(n))
  -           !info%psp_v(n,lam)= (vnl(n,lam,1)+vpsloc(n))
  -           !info%psp_v(n,lam)= (vnl(n,lam,1))
  -           !info%psp_v(n,lam)= -2*zval+vnl(n,lam,1)
  -           !info%psp_v(n,lam)= -2*zval
  +       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,ndmx
  +         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
  -       jcut=1
  -
  -       !smooth the curve
  -         !do n=jcut,ndmx
  -         do n=1,ndmx
  -           !info%psp_v(n,lam)=-2*zval/m%r(n)+cutoff_function(m%r(n)-m%r(jcut))*(vnl(n,lam,1)/m%r(n)+vpsloc(n)+2*zval/m%r(n))
  -           !info%psp_v(n,lam)=-2*zval+cutoff_function(m%r(n)-m%r(jcut))*(info%psp_v(n,lam))
  +
  +       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
   
  
positive.code:         !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
  
negative.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(unit,'(1X,2I3,I5,3F20.10)') n_dn, n_up, info%m%np, info%m%b, info%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)
          write(unit,'(4(g20.12))') (info%m%r(i),i = 1, info%m%np)
          allocate(dum(new_m%np))
      
          !Down pseudopotentials
  