HW1 Script 2

Script 2

Potential_Space_Evaluatin.m

%% Coarse potential space evaluation

clear all
fID_Progress = fopen('Progress_Script', 'a');
fprintf(fID_Progress, 'Potential_Space_Evaluation_Begin\n');
tic
diary on
%Expermental data
Ec = 3.43;
%DFT data - Data from your calculations of E-V curves
a0 = 4.05;
B = 82/160.21765;
V0 = a0^3/4;

%calculate alpha
alpha = sqrt(9*B*V0/Ec);

%Modify goal values according to your groups DFT results
goal(1) = 3.43;    % Cohesive energy
goal(2) = 4.05;    % Lattice parameter
goal(3) = 16.61;  % Volume
goal(4) = 107;    % c11
goal(5) = 61;    % c12
goal(6) = 28;    % c44
goal(7) = 82;       % Bulk modulus
goal(8) = 890;      % Surface energy of 100
goal(9) = 960;      % Surface energy of 110
goal(10) = 780;     % Surface energy of 111
goal(11) = 133;     % Intrinsic generalized stacking fault energy
goal(12) = 133;     % Extrinsic generalized stacking fault energy
goal(13) = 0.5;     % Vacancy formation energy
goal(14) = 2.8;     % Octahedral interstitial energy
goal(15) = 3.68;       % Tetrahedral interstitial energy

save('goal.mat', 'goal'); %Saves the goal values so they can be used later

%Specify initial potential and bounds for parameters
%X(1:12)    = [alpha,    asub,  b0 ,  b1 ,  b2 , b3  , t0 ,   t1 ,  t2  , t3  , cmax, cmin]
default     = [alpha,    1.07, 2.21, 2.20,  6.0, 2.2 , 1.0, -1.78, -2.21, 8.01,  2.8,  2.0];
low         = [alpha*.8 , 0.8,  0.0,  0.0,  0.0, 0.0 , -10,  -10 ,   -10,  -10,  2.0,  0.0];
high        = [alpha*1.2, 1.2,   10,   10,   10,  10 ,  10,   10 ,   110,   10,  2.8,  2.0];

%Define no. of points
N = 10;
%Calculate interval
incr = (high-low)/N;

%Construct the potential database
clear X
count = 0;
for i = 1:length(low)
    for j = low(i):incr(i):high(i)
        count = count + 1;
        X(count,:) = default;
        X(count,i) = j;
    end
end
    

% Specify initial starting parameters

clear parameter

parameter.elt       = 'Al';
parameter.lat       = 'fcc';
parameter.z         = 12.0;
parameter.ielement  = 13.0;
parameter.atwt      = 26.9815;
parameter.rozero    = 1.0;
parameter.ibar      = 0;
parameter.rcut      = 5.0;
parameter.alat      = a0;
parameter.esub      = Ec;

parameter.alpha     = default(1);
parameter.asub      = default(2);
parameter.b0        = default(3);
parameter.b1        = default(4);
parameter.b2        = default(5);
parameter.b3        = default(6);
parameter.t0        = default(7);
parameter.t1        = default(8);
parameter.t2        = default(9);
parameter.t3        = default(10);
parameter.cmax      = default(11);
parameter.cmin      = default(12);

save('parameter.mat','parameter');

%Write out initial potential files
init_filename_lib = 'library.meam.init';
init_filename_alloy = 'Al.meam.init';

write_lmpmeam(init_filename_lib,init_filename_alloy, parameter);
%% Run Lammps for each potential and obtain responses

%Intialize database to store responses
N_res = 15; % No. of responses per calculation (No. of target properties)
Y = zeros([size(X,1),N_res]);

%These two file names must be same with the pair_coeff command in the bb_script.py
filename_lib = 'library.meam';
filename_alloy = 'Al.meam';

fprintf(fID_Progress, '\t%f\n', toc);
fprintf(fID_Progress, 'Begin_Loop\n');
'1' 
toc
for i = 1:2

    p = parameter;

    % Make changes
    
    p.alpha     = X(i,1);
    p.asub      = X(i,2);
    p.b0        = X(i,3);
    p.b1        = X(i,4);
    p.b2        = X(i,5);
    p.b3        = X(i,6);
    p.t0        = X(i,7);
    p.t1        = X(i,8);
    p.t2        = X(i,9);
    p.t3        = X(i,10);
    p.cmax      = X(i,11);
    p.cmin      = X(i,12);

    write_lmpmeam(filename_lib,filename_alloy,parameter);
'2'
toc        
    system('./bb_script.py');
'3'
toc
    p_Ec        = find_energy('responses.txt', 'E_coh');
    p_a0        = find_energy('responses.txt', 'a_0');
    p_V0        = find_energy('responses.txt', 'V_0');
    p_c11       = find_energy('responses.txt', 'C11');
    p_c12       = find_energy('responses.txt', 'C12');    
    p_c44       = find_energy('responses.txt', 'C44');
    p_B         = find_energy('responses.txt', 'B');
    p_E_100     = find_energy('responses.txt', 'E_100');
    p_E_110     = find_energy('responses.txt', 'E_110');
    p_E_111     = find_energy('responses.txt', 'E_111');
    p_GSFE_I    = find_energy('responses.txt', 'GSFE_I');
    p_GSFE_E    = find_energy('responses.txt', 'GSFE_E');
    p_Evac      = find_energy('responses.txt', 'E_vac');
    p_Eoct      = find_energy('responses.txt', 'E_Oct');
    p_Etet      = find_energy('responses.txt', 'E_Tetra');
'4'
toc
    Y(i,:) = [p_Ec, p_a0, p_V0, p_c11, p_c12, p_c44, p_B, p_E_100, p_E_110,... 
              p_E_111, p_GSFE_I, p_GSFE_E, p_Evac, p_Eoct, p_Etet];
    
    fprintf(fID_Progress, '%d\t%f\n', i, toc);
    
end
'5'
toc
save('Al_coarse_pot_eval_data.mat','X','Y','default','low','high','incr')

%% Generate files for plots

clear all
load ('Al_coarse_pot_eval_data.mat')
N = 10; %No. of points defined at the beginning
fID_Progress = fopen('Progress_Script', 'a');
%Goal values must be assigned at the beginning of the file
load('goal.mat'); %Load the saved goal values

pname = {'alpha', 'a_sub', 'b0', 'b1', 'b2', 'b3', 't0', 't1', 't2', ...
        't3', 'cmax', 'cmin'};
    

index = 1;
for i = 1:length(low)
    fname = strcat(pname(i), '_para.txt');
    fid = fopen(cell2mat(fname),'wt');
    count = 1;
    while count <= (N+1)
        fprintf(fid,'%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n',X(index,i), ...
            Y(index,1)-goal(1),Y(index,2)-goal(2),Y(index,3)-goal(3),Y(index,4)-goal(4),Y(index,5)-goal(5),...
            Y(index,6)-goal(6),Y(index,7)-goal(8),Y(index,9)-goal(9),Y(index,10)-goal(10),Y(index,11)-goal(11),...
            Y(index,12)-goal(12),Y(index,13)-goal(13),Y(index,14)-goal(14),Y(index,15)-goal(15));
        
        count = count+1;
        index = index + 1;  
    end
    fclose(fid);
    
end

fprintf(fID_Progress, '\t%f\n', toc);
fprintf(fID_Progress, 'Potential_Space_Evaluation_End\n');

fclose('all');