function [cpi] = read_fhi_psp(name)
% 
% [cpi] = read_fhi_psp(name)
%
% Read in the pseudopotential in FHI98MD format.
%
%  cpi.[]
%	type = fhi,
%	Zv, lmax
%	ngrid(:), clog(:),
%	RS(:), rPsiS(:), VionS(:),
%	RP(:), rPsiP(:), VionP(:),
%	RD(:), rPsiD(:), VionD(:),
%	core(:,1:4);

cpi.type = 'fhi';

fid = fopen(name,'r');

line = fgetl(fid);
data = sscanf(line,'%g',2);
cpi.Zv = data(1);
cpi.lmax = data(2);

if (cpi.lmax > 3)
  error('Do not know how to handle more than s, p, and d potentials');
end

%% skip ten lines of numbers for parameter-fitted psp.
for i=1:10,
  fgetl(fid);
end

%% for each angular momentum, read in data
for i=1:cpi.lmax,
  data = fscanf(fid,'%g',2);
  cpi.ngrid(i) = data(1);
  cpi.clog(i) = data(2);
  data = fscanf(fid,'%g',4*cpi.ngrid(i));
  data2 = zeros(4,cpi.ngrid(i));
  data2(:) = data;
  data2 = data2';
  if i==1
    cpi.RS = data2(:,2);
    cpi.rPsiS = data2(:,3);
    cpi.VionS = data2(:,4);
  elseif i==2
    cpi.RP = data2(:,2);
    cpi.rPsiP = data2(:,3);
    cpi.VionP = data2(:,4);
  else
    cpi.RD = data2(:,2);
    cpi.rPsiD = data2(:,3);
    cpi.VionD = data2(:,4);
  end    
end

%% read in core electron density
[data, count] = fscanf(fid,'%g',4*cpi.ngrid(cpi.lmax));
if count==(4*cpi.ngrid(cpi.lmax))
  cpi.core = zeros(4,cpi.ngrid(cpi.lmax));
  cpi.core(:) = data;
  cpi.core = cpi.core';
else
  cpi.core = '';
end

cpi.i = 1:1:cpi.ngrid(cpi.lmax);

