I_Eckart

function I_Eckart

dihedral=input('\n Please input the dihedral angle specified by 4 atoms: ');

symmetry=input('\n Please input the symmetry number of the torsional potential: ');

angle_interval=input('\n Please input the stepsize in degrees between each geometry optimization: ');

angles_degrees=input('\n Please input a vector of angles in degrees to calculate I_Eckart: ');

mass_matrix=loadascii('mass.txt');
cartesian_matrix=loadascii('cart.txt');

N=length(mass_matrix);

intervals=2*angle_interval;
angles=pi/180*angles_degrees;

center_of_mass_matrix=center_of_mass_frame(mass_matrix,cartesian_matrix);

[gaussian_matrix,gaussian_rot_matrix]=gaussian_frame(center_of_mass_matrix,dihedral);

[gaussian_matrix,permutations]=permutation(gaussian_matrix,symmetry);

mass_weight_matrix=mass_weight_frame(mass_matrix,gaussian_matrix);

[a,eckart_rot_matrix]=eckart(mass_weight_matrix,intervals);

diff_a_of_phi=a_prime_of_phi(a);

mu_ss=inverse_inertia(a,diff_a_of_phi,permutations);

coef_inertia=four_coef_phi(1./mu_ss');

max_exponent=(length(coef_inertia)-1)/2;

exponents=[0 linspace(-1,-max_exponent,max_exponent) linspace(max_exponent,1,max_exponent)]';

for j=1:length(angles)

    inertia_vector(j)=real(sum(coef_inertia.*exp(i*exponents*angles(j))));

end

disp(' ')
disp(' For the torsional angles you entered (column 1), ')
disp(' the Eckart inertias in amu*Angstrom^2 (column 2) are:')
disp(' ')
disp(' Torsional Angle (degrees)    Eckart Inertia (amu*Angstrom^2)')
disp(' ')
fprintf(1, '  %14.2f %34.10G \n',[angles_degrees;inertia_vector])