function tunable_potentials_3D
clear;clc;
format long
%---- This script computes and print 3D tunable potentials deveopled by Rajan et al.
%---- 1. computes A1, B1, C1, ... depending on r1, r2, r3, phi(r2), which satisfy the continuity
%---- of thbetween Mores function, spline 1 and spline 2.


alpha = 13/(2^(5/6)); % --- the values defined here is taken from Rajan and Curtin ----

energy = -1*[0.318 0.266 0.207 0.143 0.077 0.009]; % these are desired values of phi(r2) ---- This values are again taken from Rajan et. Curtin
radius_3 = [1.391 1.391 1.391 1.360 1.304 1.247]; % values taken from Rajan and Curtin
[d1 d2]=size(energy);
parameters = [];
for j=1:d2
	% ----- Connection points between the curves of interest --------
	x1 = 1.05;
	x2 = sqrt(3/2);
	x3 = radius_3(1,j);

	% ------ y values of the interest ------
	y1 = (1-exp(-alpha*(x1-1)))^2-1; % ----eq(1) Morse potential and its vale for x1
	y1_slope = 2*alpha*(1-exp(-alpha*(x1-1)))*exp(-alpha*(x1-1)); % ---- First derivative of the eq(1)
	%%% ----- change this!!!
	y2 = energy(1,j); % ----- comes from the input
	y3 = 0;
	y3_slope = 0;
	%-------------------------------------------------------------------------------------------

	Y = transpose([y1_slope, y1, y2, 0, 0, y2, y3, y3_slope]);
	M = [3*x1^2 2*x1 1 0 0 0 0 0;...
		 x1^3 x1^2 x1 1 0 0 0 0; ...
		 x2^3 x2^2 x2 1 0 0 0 0;...
		 3*x2^2 2*x2 1 0 -3*x2^2 -2*x2 -1 0; ...
		6*x2 2 0 0 -6*x2 -2 0 0; ...
		0 0 0 0 x2^3 x2^2 x2 1; ...
		0 0 0 0 x3^3 x3^2 x3 1; ...
		0 0 0 0 3*x3^2 2*x3 1 0];

	COEF = inv(M)*Y;

	%------- Coefficients for the desired potential ------


	A1 = COEF(1,1);
	B1 = COEF(2,1);
	C1 = COEF(3,1);
	D1 = COEF(4,1);

	A2 = COEF(5,1);
	B2 = COEF(6,1);
	C2 = COEF(7,1);
	D2 = COEF(8,1);
	parameters1 = [A1 B1 C1 D1 A2 B2 C2 D2];
	parameters = [parameters;parameters1];
	% ---------------------------------------------------

	% ------- Computation of the potential and the force ------
	R1 = 0.00000001;
	r1 = x1;
	r2 = x2;
	r3 = radius_3(1,j);
	rnn = sqrt(2);

	delta=1.9/3000;
	r = 0.1:delta:rnn;
	r = transpose([R1,r]);

	[d1 d2]=size(r);

	% ---- Energy part -----
		E = [];
		for i = 1:d1
			if r(i,1) <=r1
				E1 = (1-exp(-alpha*(r(i,1)-1)))^2-1;;
				E = [E;E1];
			elseif r(i,1) <=r2
				E1 = A1*(r(i,1))^3 + B1*(r(i,1))^2 + C1*(r(i,1)) + D1;
				E = [E;E1];
			elseif r(i,1) <=r3
				E1 = A2*(r(i,1))^3 + B2*(r(i,1))^2 + C2*(r(i,1)) + D2;
				E = [E;E1];
			elseif r(i,1) > r3
				E1 = 0;
				E = [E;E1];
			end
		end

	% ---- Force part -----
		F = [];
		for i = 1:d1
			if r(i,1) <=r1
				F1 = 2*alpha*(1-exp(-alpha*(r(i,1)-1)))*exp(-alpha*(r(i,1)-1));
				F = [F;F1];
			elseif r(i,1) <=r2
				F1 = 3*A1*(r(i,1))^2 + 2*B1*(r(i,1)) + C1;
				F = [F;F1];
			elseif r(i,1) <=r3
				F1 = 3*A2*(r(i,1))^2 + 2*B2*(r(i,1)) + C2;
				F = [F;F1];
			elseif r(i,1) > r3
				F1 = 0;
				F = [F;F1];
			end
		end
		F1=F;
		F=-1*F; % ---- negative/ positive force definition .....


	% -------- Writing the potential into file ------------------
	txtfile=['Pair_potential_FCC_3D_' num2str(j) '.txt'];
	fileID=fopen(txtfile,'w');
	fprintf(fileID,'\n');
	fprintf(fileID,'\n');
	fprintf(fileID,'# Tabulated potential, r1 = 1.05, r2 = sqrt(3/2), r_nn = sqrt(2) , alpha = 13/2^(5/6)\n');
	fprintf(fileID,'\n');
	fprintf(fileID,'dummy\n');
	fprintf(fileID,'N 3002\n');
	fprintf(fileID,'\n');
	for i=1:d1
		fprintf(fileID,'%d %10.8f %10.8f %10.8f\n', i, r(i,1), E(i,1), F(i,1));
	end
	fclose(fileID);

	figure(1)
	plot(r,E,'-','LineWidth',1.5)
	axis([0.8 1.5 -1.01 1.0])
	set(gca,'FontSize',14)
	set(gca,'FontName','Times New Roman')
	set(gca,'XTick',[0.8:0.1:1.5],'LineWidth',1.5)
	set(gca,'YTick',[-1:0.5:1.0],'LineWidth',1.5)
	xlabel('Interatomic separation, r','FontSize',18)
	ylabel('Energy','FontSize',18)
	hold on

	figure(2)
	plot(r,F1,'-','LineWidth',1.5)
	axis([0.8 1.5 -1.01 7.0])
	set(gca,'FontSize',14)
	set(gca,'FontName','Times New Roman')
	set(gca,'XTick',[0.8:0.1:1.5],'LineWidth',1.5)
	set(gca,'YTick',[-1:1:7.0],'LineWidth',1.5)
	xlabel('Interatomic separation, r','FontSize',18)
	ylabel('(Negative)Force','FontSize',18)
	hold on

end
hold off

% --- Printing parameters in txt file for further use -----

[s1 s2] = size(parameters);
fileID=fopen('3D_parameters_tunable_potentials.txt','w');
fprintf(fileID,'# -- Parameters for the 3D potentials with tunable ductility\n');
fprintf(fileID,'# -- Ordered from the most ductile to the most brittle\n#\n#\n');
fprintf(fileID,'alpha = %20.15f\n',alpha)
fprintf(fileID,'r1    = %20.15f\n',x1)
fprintf(fileID,'r2    = %20.15f\n',x2)
fprintf(fileID,'A1 B1 C1 D1 A2 B2 C2 D2\n');
for jj = 1:s1
		for ii =1:s2
		fprintf(fileID,'%20.15f ',parameters(jj,ii));
		end
		fprintf(fileID,' \n');
end
fclose(fileID);