% Function : Solving the Hessian Matrix by eigs

function [] = eigensolver ( pdbid, EIGM );

disp('The PDB ID is : ')
pdbid

disp('Now loading the Hessian Matrix ... ... ')
Hessian = load ([pdbid,'_Hessian.dat']);
disp('Successfully loading the Hessian Matrix')
S = spconvert(Hessian);


options.issym=1;
options.isreal=1;
options.tol=1.0e-6;
options.maxit=1000;
[V,D] = eigs(S, EIGM+7, 'sm', options);

Dag = diag(D);

for i = 1 : EIGM+6
   eigvals(i,1) = i;
   eigvals(i,2) = Dag(EIGM-i+7);
end 

sz = size(V);
SITE= sz(1)/3;

for i = 1 : SITE   
   for j = 1 : EIGM
       eigvX(i,j) = V(3*i-2,EIGM-j+1);
       eigvY(i,j) = V(3*i-1,EIGM-j+1);
       eigvZ(i,j) = V(3*i-0,EIGM-j+1);
   end
end

N=(1:SITE)';
eigvX=[N eigvX];
eigvY=[N eigvY];
eigvZ=[N eigvZ];

disp('EigenVector-X : ');
disp(eigvX(1:10,1:5))

%------------------------- Part 3 Output the eigenvalue and eigenvector
disp('Eigenvalues : ')
disp(eigvals)

fideigvals = fopen ([pdbid,'.eigvals'],'w');

fprintf(fideigvals, '%6.0f  %+12.5E \n', eigvals');

status = fclose(fideigvals);


fideigX = fopen ([pdbid,'.slwX'],'w');
fideigY = fopen ([pdbid,'.slwY'],'w');
fideigZ = fopen ([pdbid,'.slwZ'],'w');

for i = 1 : SITE
    fprintf(fideigX,'%8.0f', eigvX(i,1));
    fprintf(fideigY,'%8.0f', eigvY(i,1));
    fprintf(fideigZ,'%8.0f', eigvZ(i,1));
    for j = 2 : EIGM+1
       fprintf(fideigX,' %+9.6f', eigvX(i,j));
       fprintf(fideigY,' %+9.6f', eigvY(i,j));
       fprintf(fideigZ,' %+9.6f', eigvZ(i,j));
    end
    fprintf(fideigX,' \n');
    fprintf(fideigY,' \n');
    fprintf(fideigZ,' \n');
end

status = fclose(fideigX);
status = fclose(fideigY);
status = fclose(fideigZ);

clear;
return;