clear
clc
unitMass = (1/1.6021766208e-19)*1e-20*1e30*1.660539e-27; % eV*A^-2*fs^2
hbar = 4.1356676969/2/pi; % eV*fs
% Cu bulk
v_sound = 4760*1e10*1e-15; % A*fs^-1, annealed
density = 8.96*1e3*unitMass/1.660539e-27*1e-30; % eV*A^-5*fs^2
%% Create/Read Potential Energy
relativeMass = 10.81;
atom_num = 1;
load('p1_ho_B.mat');
load('atom_B.mat');
distance1o = atom_B(:,1);
energyo = atom_B(:,2);
dipoleo = atom_B(:,3);
range_ii = find(p1_ho.range_cutoff==8);
disp(relativeMass(atom_num))
mass = unitMass*relativeMass(atom_num);
temperature = 300; % K
steplength = 1e-3;
max_vibrationalNum = 15;
hfunc = p1_ho.hfunc;
parameter_energy = p1_ho.parameters1_min{range_ii}(end,:);
range_cutoff_energy = p1_ho.range_cutoff(range_ii);
range_fit_energy = p1_ho.range_fit;
x = (distance1o(1):steplength:distance1o(end)).';
x_cutoff_energy = distance1o(range_fit_energy(range_cutoff_energy));
x_diff_energy = abs(x-x_cutoff_energy);
loca = find(x_diff_energy==min(x_diff_energy));
yenergy = hfunc(parameter_energy,x);
yenergy = [x, yenergy];
%% Create/Read Dipole Moment
ydipole = interp1(distance1o,dipoleo,x,'spline');
ydipole = [x, ydipole];
%% Obtain Vibrational Wavefunction And Eigen Energy
[position,potentialEnergy,wavefunc,eigenEnergy] = ...
obtainVibrationalWavefunction( ...
yenergy,hbar,mass,max_vibrationalNum,steplength);
eigenEnergy_10 = eigenEnergy(2) -eigenEnergy(1);
%%
% figure; plot(position,potentialEnergy);
% hold on; plot(distance1o,energyo,'o');
%
% plot(position,wavefunc(:,end));
%% figure; plot((0:length(eigenEnergy)-1).',eigenEnergy);
% ylim([min(min(energyo),min(wavefunc(:,end)))-0.1,max(max(energyo),max(wavefunc(:,end)))+0.1])
%% Obtain Dipole Moment On Vibrational States
% unit Debye
dipoleMoment = obtainDipoleMoment(ydipole,wavefunc);
dipoleMoment_gs = dipoleMoment(1);
%%
% figure; plot(position,ydipole(:,2))
% hold on; plot(distance1o,dipoleo,'o')
% figure; bar((0:max_vibrationalNum-1).',dipoleMoment);
%% Obtain Transition Rate Using Fermi Gold Rule
% unit fs^-1
transitionRate = obtainTransitionRate(position,potentialEnergy, ...
wavefunc,eigenEnergy,temperature,hbar,v_sound,density);
%%
% figure
% rho_0 = exp(-(eigenEnergy-eigenEnergy(1))/(temperature*0.026/300));
% bar(rho_0/sum(rho_0))
%%
clearvars -except eigenEnergy dipoleMoment transitionRate temperature