View source for ICME 2012 POT-GEN-6

[[ICME_2012_HW1| < Homework 1]]

Author(s): [https://sites.google.com/site/lsiliyanage/home Laalitha Liyanage]

==License==
By using the codes provided here you accept the the [[CodeRepository:License Agreement|Mississippi State University's license agreement]]. Please read the agreement carefully before usage.

==Overview==

A method to introduce uncertainties to goal values are presented here. Uncertainty is introduce through the generation of random numbers that are in a normal distribution between the uncertainty specified.
For each set of goal values a best potential is selected and the response recorded. Therefore at the end a 1000 (or any number one wishes) potentials and responses can be evaluated.

==Modifying the script==
The sensitive parameters and their ranges should be the same as [[ICME 2012 POT-GEN-3|script-3]]. The database from sampling the parameter space should be used (output from [[ICME 2012 POT-GEN-4|script-4]]). The experimental and DFT values and the default potential used should be the same as [[ICME 2012 POT-GEN-3|script-3]].

== Running the script == 
*Copy the script to a linux file and name it ''any_name.m''. 
*Copy and save '''[[ICME_2012_POT-GEN-2#Auxiliary_scripts|auxiliary scripts]]''' with provided names.

To run the script you can execute the following command in the command prompt or insert it in a pbs command script
<pre>matlab -nodisplay -nosplash < <script_name>.m</pre>

The script will call [[ICME 2012 POT-GEN-1|script-1]]. Therefore the number of processors allocated in the pbs command script should correspond with the number of processors specified in the mpirun command in [[ICME 2012 POT-GEN-1|script-1]].

==Output of script==
The output is a database consisting of the potentials and their responses. Also it will print the mean and standard deviation of the predicted goal value.

==The script==

<pre>
%% Effect of uncertainty in Goal values
%Load data from potential space sampling
clear all
load('Al_LHS_v1.mat')
%% 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]);

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

%% Further analysis

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

</pre>