ICME 2012 POT-GEN-2

< Homework 1

Author(s): Laalitha Liyanage

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

The MATLAB script presented here will evaluate the potential space on the bounds set for each parameter of the Modified Embedded Atom Method potential.

Modifications to the script

Before using this script you should modify the data in the Experimental data and DFT data to the values calculated or found by your group.

Running the script

Copy the script to a linux file and name it any_name.m. Copy the auxiliary scripts with file names provided to the same directory. script-1 should be in the same directory as the current script. If you changed the name of script-1 be sure to change it in the current script as well. 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 named bb_script.py. 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.

Potential space evaluation script

%% Coarse potential space evaluation

clear all

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

%Modify goal values according to your groups DFT results
goal(1) = 3.43;    % Cohesive energy
goal(2) = 4.05;    % Lattice parameter
goal(3) = 16.61;  % Volume
goal(4) = 107;    % c11
goal(5) = 61;    % c12
goal(6) = 28;    % c44
goal(7) = 82;       % Bulk modulus
goal(8) = 890;      % Surface energy of 100
goal(9) = 960;      % Surface energy of 110
goal(10) = 780;     % Surface energy of 111
goal(11) = 133;     % Intrinsic generalized stacking fault energy
goal(12) = 133;     % Extrinsic generalized stacking fault energy
goal(13) = 0.5;     % Vacancy formation energy
goal(14) = 2.8;     % Octahedral interstitial energy
goal(15) = 3.68;       % Tetrahedral interstitial energy

save('goal.mat', 'goal'); %Saves the goal values so they can be used later

%Specify initial potential and bounds for parameters
%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];
low         = [alpha*.8 , 0.8,  0.0,  0.0,  0.0, 0.0 , -10,  -10 ,   -10,  -10,  2.0,  0.0];
high        = [alpha*1.2, 1.2,   10,   10,   10,  10 ,  10,   10 ,   110,   10,  2.8,  2.0];

%Define no. of points
N = 10;
%Calculate interval
incr = (high-low)/N;

%Construct the potential database
clear X
count = 0;
for i = 1:length(low)
    for j = low(i):incr(i):high(i)
        count = count + 1;
        X(count,:) = default;
        X(count,i) = j;
    end
end
    

% Specify initial starting parameters

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

%Write out initial potential files
init_filename_lib = 'library.meam.init';
init_filename_alloy = 'Al.meam.init';

write_lmpmeam(init_filename_lib,init_filename_alloy, parameter);

%% Run Lammps for each potential and obtain responses

%Intialize database to store responses
N_res = 15; % No. of responses per calculation (No. of target properties)
Y = zeros([size(X,1),N_res]);

%These two file names must be same with the pair_coeff command in the bb_script.py
filename_lib = 'library.meam';
filename_alloy = 'Al.meam';

for i = 1:size(X,1)

    p = parameter;

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

    write_lmpmeam(filename_lib,filename_alloy,parameter);
        
    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];

    
end

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

%% Generate files for plots

clear all
load ('Al_coarse_pot_eval_data.mat')
N = 10; %No. of points defined at the beginning

%Goal values must be assigned at the beginning of the file
load('goal.mat'); %Load the saved goal values

pname = {'alpha', 'a_sub', 'b0', 'b1', 'b2', 'b3', 't0', 't1', 't2', ...
        't3', 'cmax', 'cmin'};
    

index = 1;
for i = 1:length(low)
    fname = strcat(pname(i), '_para.txt');
    fid = fopen(cell2mat(fname),'wt');
    count = 1;
    while count <= (N+1)
        fprintf(fid,'%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n',X(index,i), ...
            Y(index,1)-goal(1),Y(index,2)-goal(2),Y(index,3)-goal(3),Y(index,4)-goal(4),Y(index,5)-goal(5),...
            Y(index,6)-goal(6),Y(index,7)-goal(8),Y(index,9)-goal(9),Y(index,10)-goal(10),Y(index,11)-goal(11),...
            Y(index,12)-goal(12),Y(index,13)-goal(13),Y(index,14)-goal(14),Y(index,15)-goal(15));
        
        count = count+1;
        index = index + 1;  
    end
    fclose(fid);
    
end

Auxiliary scripts

write_lmpmeam.m

%This function writes the file equivalent to meafile(dynamo) in LAMMPS
function write_lmpmeam(filename1,filename2,parameter)

fid = fopen(filename1,'wt');

fprintf(fid,'  #MEAM data\n'); 
fprintf(fid,'  # elt       lat     z       ielement     atwt\n');
fprintf(fid,'  # alpha     b0      b1      b2           b3       alat    esub    asub\n');
fprintf(fid,'  # t0        t1              t2           t3               rozero  ibar\n\n');

fprintf(fid,'''%s'' ''%s'' %u %u %f \n',parameter.elt, parameter.lat, ...
            parameter.z, parameter.ielement, parameter.atwt);
fprintf(fid,'%f %f %f %f %f %f %f %f \n',parameter.alpha, parameter.b0, ...
            parameter.b1, parameter.b2, parameter.b3, parameter.alat, ...
            parameter.esub, parameter.asub);
fprintf(fid,'%d %f %f %f %d %d \n', parameter.t0, parameter.t1, parameter.t2, ...
            parameter.t3, parameter.rozero, parameter.ibar);
        
fclose(fid);

fid = fopen(filename2, 'wt');

fprintf(fid,'augt1=0\nialloy=1\nnn2(1,1)=1\n');
fprintf(fid,'rc=%f \n', parameter.rcut);
fprintf(fid,'rho0(1)=%f\n',parameter.rozero);

fclose(fid); 

find_energy.m

function [var] = find_energy(outfile, str)

var = 0;
fid = fopen(outfile);

tline = fgetl(fid);
while ~feof(fid)
   matches = strfind(tline, str);
   num = length(matches);
   if num > 0
       C = textscan(tline, '%*s %f %f');
       var = C{1}; 
       %fprintf(1,'%s \n',tline);
   end
   tline = fgetl(fid);
end

fclose(fid);

Remove_files

This script will delete the output files that are created by the scripts. Do not run this script unless you plan to restart your calculations! Note that it should be run from a terminal.

rm c11_fit.log c12_fit.log c44_fit.log c11_fit.gif c12_fit.gif c44_fit.gif c11_fit.gp c12_fit.gp
rm c44_fit.gp Al.meam Al.meam.init alpha_para.txt c11_report.lammps c12_report.lammps
rm c44_report.lammps datafile dump.E0 dump.gsfe.E dump.gsfe.I dump.oct dump.surf100 dump.surf110
rm dump.surf111 dump.tet dump.vac gnuplot.report infile_relax library.meam library.meam.init log.E0
rm log.gsfe.E log.gsfe.I log.lammps log.oct log.surf100 log.surf110 log.surf111 log.tet log.vac
rm report.lammps responses.txt summary Al_coarse_pot_eval_data.mat parameter.mat Al_LHS_v1.mat
rm Al_Uncertainty_Ecoh.mat Al_weights.mat *_para.txt goal.mat objfun*.mat objfunw.mat
rm t*_para.txt log Progress_Script dump

OUTPUT of script

The script will output a Al_coarse_pot_eval_data.mat file and *_para.txt files. Each *_para.txt file will have the responses for a single parameter variation. Format of *_para.txt file will be as follows

Parameter_value_1 response_1 response_1 ... response_15
Parameter_value_2 response_1 response_1 ... response_15
.
.
.
Parameter_value_N response_1 response_1 ... response_15