HW1 Script 6

Script 6

Evaluation_of_uncertainties_in_goal_values.m

%% Effect of uncertainty in Goal values
%Load data from potential space sampling
clear all
load('Al_LHS_v1.mat')

fID_Progress = fopen('Progress_Script', 'a');
fprintf(fID_Progress, 'Evaluation_of_Uncertainties_in_Goal_Values_Begin\n');
tic
%% Specify sensitive parameters and their bounds
%------------THIS SHOULD BE THE SAME AS IN 02_Al_LHS_run.m---------------------------------
%Example specifies all parameters and maximum ranges
%X(1:12)    = [alpha,    asub,  b0 ,  b1 ,  b2 , b3  , t0 ,   t1 ,  t2  ,t3  , cmax, cmin]
low         = [3.78 , 0.8,  0.0,  0.0,  0.0, 0.0 , -10,  -10 ,   -10,  -10,  2.0,  0.0]; 
high        = [5.67, 1.2,   10,   10,   10,  10 ,  10,   10 ,   110,   10,  2.8,  2.0];

%% Initialize potential parameters to default in 02_Al_LHS_run.m
%------------THIS SHOULD BE THE SAME AS IN 02_Al_LHS_run.m---------------------------------
%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);

%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];

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

filename1 = 'library.meam.init';
filename2 = 'Al.meam.init';
write_lmpmeam(filename1,filename2,parameter);

%% Apply uncertainty
%Set goal values at the beginning of Potential Space Evaluation Script
load('goal.mat'); %Load values for goal
default_goal = goal;

pname = {'Ecoh', 'a0', 'V0', 'C11', 'C12', 'C44', 'E100', 'E110', 'E111', 'GSFE_I', 'GSFE_E', 'Evac', 'Eoct', 'Etetra'};

% Weighting
% Use weighting determined from 04_Al_weights.m
default_w = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];

%Introduce uncertainty to goal value ------------ ONE AT A TIME-------------
%Through a truncated normal distribution a 1000 random values are generated
k = 1; %specifies goal to which uncertainty is applied
g_u = normrnd(goal(k),0.1/3.0,1,1000);%first value is mean (experimental),
                                      %next is standard deviation (use 1/3 of the uncertainty percentage 
                                      %to generate 99% of no. within the range)

%Generate 1000 goal value combinations to get 1000 potentials
clear G
count = 0;
for i = 1:length(g_u)
        count = count + 1;
        G(count,:) = default_goal;
        G(count,k) = g_u(i);
end


Y = zeros([size(G,1),15]);
xsol = zeros([size(G,1),size(X,2)]);
filename1 = 'library.meam';
filename2 = 'Al.meam';

load('parameter.mat');

goal = zeros([1 15]);

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

for j=1:size(G,1)
    
    
    clear w x0 lb ub
    w= default_w;
    goal = G(j,:);
    save('objfunw.mat','w','goal')
    
    %optimization
    
    load('Optim.mat')
    x0(1:size(X,2)) = 0.5;
    lb = zeros(size(x0));
    ub = ones(size(x0));
    optimproblem.x0 = x0;
    optimproblem.lb = lb;
    optimproblem.ub = ub;
    xsol(j,:) = fmincon(optimproblem);
    
    % Make changes
    p = parameter;

    p.alpha     = low(1)+(high(1)-low(1))*xsol(j,1);
    p.asub      = low(1)+(high(2)-low(2))*xsol(j,2);
    p.b0        = low(1)+(high(3)-low(3))*xsol(j,3);
    p.b1        = low(1)+(high(4)-low(4))*xsol(j,4);
    p.b2        = low(1)+(high(5)-low(5))*xsol(j,5);
    p.b3        = low(1)+(high(6)-low(6))*xsol(j,6);
    p.t0        = low(1)+(high(7)-low(7))*xsol(j,7);
    p.t1        = low(1)+(high(8)-low(8))*xsol(j,8);
    p.t2        = low(1)+(high(9)-low(9))*xsol(j,9);
    p.t3        = low(1)+(high(10)-low(10))*xsol(j,10);
    p.cmax      = low(1)+(high(11)-low(11))*xsol(j,11);
    p.cmin      = low(1)+(high(12)-low(12))*xsol(j,12);
    
    write_lmpmeam(filename1,filename2,p);
        
    system('./bb_script.py');
    
    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');
    
    Y(j,:) = [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('Goal: %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n',G(j,:));
    fprintf('Xsol: %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n',Y(j,:));
    
    fname = strcat('Al_Uncertainty_',cell2mat(pname(k)),'.mat');
    save(fname,'X','Y','G');
    
    fprintf(fID_Progress, '%d\t%f\n', j, toc);

end

%% Further analysis

m = mean(Y(:,k));
sd = std(Y(:,k));

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

fclose('all');