HW1 Script 3

Script 3

Potential_Space_Sampling.m

%% Sampling the potential space
clear all
fID_Progress = fopen('Progress_Script', 'a');
fprintf(fID_Progress, 'Potential_Space_Sampling_Begin\n');
tic
%Specify sensitive parameters and their bounds
%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];

% Specify no. of samples
ndiv = 100;
% Specify no. of sensitive parameters
npar = 12;
% Design LHS
X = lhsdesign(ndiv,npar,'criterion','correlation');

%Initialize database for responses
Y = zeros([size(X,1),15]);

%% Initialize potential parameters to default in 01_Al_Coarse_pot_eval.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);

%% Run LAMMPS and get responses for sample potentials
filename1 = 'library.meam';
filename2 = 'Al.meam';

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

    p = parameter;

    % Make changes
    
    p.alpha     = low(1)+(high(1)-low(1))*X(i,1);
    p.asub      = low(1)+(high(2)-low(2))*X(i,2);
    p.b0        = low(1)+(high(3)-low(3))*X(i,3);
    p.b1        = low(1)+(high(4)-low(4))*X(i,4);
    p.b2        = low(1)+(high(5)-low(5))*X(i,5);
    p.b3        = low(1)+(high(6)-low(6))*X(i,6);
    p.t0        = low(1)+(high(7)-low(7))*X(i,7);
    p.t1        = low(1)+(high(8)-low(8))*X(i,8);
    p.t2        = low(1)+(high(9)-low(9))*X(i,9);
    p.t3        = low(1)+(high(10)-low(10))*X(i,10);
    p.cmax      = low(1)+(high(11)-low(11))*X(i,11);
    p.cmin      = low(1)+(high(12)-low(12))*X(i,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(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

save('Al_LHS_v1.mat','X','Y','default','low','high')

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

fclose('all');

%{
%% Further analysis...

clear all
load('Al_LHS_v1.mat')

% Select potential based on criteria and objective function

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

%Example weight
w = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];

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

X(nmin(1),:);
Y(nmin(1),:);

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

i=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);
%}