ICME 2012 POT-GEN-5
From EVOCD
Author(s): Laalitha Liyanage
Contents |
License
By using the codes provided here you accept the the Mississippi State University's license agreement. Please read the agreement carefully before usage.
Overview
Here the goal values are weighted to find the best potential parameters that fit most of the goal values. A multi-objective function is used to evaluate the response of the potentials.
Modification
The sensitive parameters and their ranges found in script-2 should be used here. Also the name of the data file generated in script-3 should be the same file name in the first load command. The experimental and DFT values and the default potential used should be the same as script-3.
Running the script
- Copy the script to a linux file and name it any_name.m.
- Copy and save auxiliary scripts with provided names.
- Copy and save mat file with provided name.
- Copy and save objfun1.m with provided name.
To run the script you can execute the following command in the command prompt or insert it in a pbs command script
matlab -nodisplay -nosplash < <script_name>.m
The script will call 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 script-1.
Output
The output is a database that has responses to combination different weights. Using the further analysis section you could generate the best potential file from the database. The new files will be named library.meam.new and Al.meam.new.
Weighting and multi-objective optimization
%% Determining weights for target properties
%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
%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');
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');
end
%{
%% 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);
%}
