HW1 Script 5

Script 5

Weighting_and_multi_objective_optimization.m

%% Determining weights for target properties
%Load data from potential space sampling
clear all
load('Al_LHS_v1.mat')

fID_Progress = fopen('Progress_Script', 'a');
fprintf(fID_Progress, 'Weight_and_Multi-Objective_Optimization_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

%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);
%load('parameter.mat')
%modify parameter

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

%% Multi-objective optimization of response surface

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

% Weighting
%Default weight
default_w = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];
%No. of points
ndiv = 10;
%Weight range for each target
lo  = zeros([1 15]);
hi = [10,10,10,10,10,10,10,10,10,10,10,10,10,10,10];
incr = (hi-lo)/ndiv;

clear Z Z_norm

count = 0;
for i = 1:length(lo)
    for j = lo(i):incr(i):hi(i)
        count = count + 1;
        Z(count,:) = default_w;
        Z(count,i) = j;
        
    end
end

Z_norm=Z;

for i = 1:size(Z_norm,1)
    for j=1:size(Z_norm,2)
        Z_norm(i,j)=Z_norm(i,j)/goal(j);
        Z_norm(i,j)=Z_norm(i,j)/sum(Z_norm(i,:));
    end
end

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

load('parameter.mat');

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

for j=1:size(Z,1)
    
    
    clear w x0 lb ub
    w= Z_norm(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];    



    goal = [0.29 152.34 63.10 175.7 140.7 0.359 10.270 0.010 9.710 4.91 6.63 4.38 2.140 1.250];
    fprintf('Goal: %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n',goal);
    fprintf('Xsol: %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n',Y(j,:));
    
    save('Al_weights.mat','X','Y','Z','Z_norm','xsol','lo','hi','incr');
    
    fprintf(fID_Progress, '%d\t%f\n', j, toc);
    
end


fprintf(fID_Progress, '\t%f\n', toc);
fprintf(fID_Progress, 'Weight_and_Multi-Objective_Optimization_End\n');

fclose('all');

%{
%% Further analysis...

clear all
load('Al_weights.mat')

% Select potential based on criteria and objective function

%Modify goal values according to your groups DFT results
%Order 
%Cohesive energy, Lattice parameter, volume, c11, c12, c44, Bulk
%modulus surface energies of 100, 110, 111, intrinsic generalized stacking
%fault energy, extrinsic generalized stacking fault energy, vacancy
%formation energy, octahedral interstitial energy, tetrahedral interstitial
%energy.
%Example goal
goal = [3.43, 4.05, 16.61,107,61,28,82,890,960,780,133,133,0.5,2.8,3,3];
%Example weight
w = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];

clear fobj
for i = 1:length(n)
    fobj(i,:) = w.*(Y(n(i),:) - goal).^2;
end
fobj = sum(fobj,2);
nmin = find(fobj == min(fobj));

Y(n(nmin(1)),:)
X(n(nmin(1)),:)
Z(n(nmin),:)

%If Y(n(nmin(1)),:) is not good enough make
%default_w = Z(n(nmin),:) and run again

% Make New potential files from X(i,:)

i=n(nmin(1));
load('parameter.mat')
p = parameter;

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

write_lmpmeam('library.meam.new','Al.meam.new',p);


%}