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);