-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfind_sites.m
225 lines (218 loc) · 10.2 KB
/
find_sites.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
function [sites, finished] = find_sites(sim_data, material, nr_parts)
% Given the locations of sites determine at which site the atoms are.
% Materials known: 'argyrodite', 'latp', 'LiSnPS', 'na3ps4' and
% 'li3ps4_beta', 'MnO2_lambda'
disp('Finding site occupations for given structure')
finished = false; % Check if the function has finished correctly
%% Get the names and positions for the given material:
[names, sites_pos, supercell] = known_materials(material);
sites.names = unique(names);
sites.supercell = supercell;
%% Define some things:
dist_close = 2*sim_data.vibration_amp; % 'amplitude of vibration'
nr_sites = size(sites_pos,2);
part_size = ceil(sim_data.nr_steps/nr_parts);
lattice = sim_data.lattice;
%% The cartesian positions of the given sites:
cart_pos = zeros(3, nr_sites);
for i = 1:nr_sites
cart_pos(:,i) = frac_to_cart(sites_pos(:,i), lattice);
end
%% Names of the sites:
nr_per_name = zeros(size(sites.names, 1),1);
for i = 1:size(names, 1)
for j = 1:size(sites.names, 1)
if strcmp(sites.names{j}, names{i})
nr_per_name(j) = nr_per_name(j) + 1;
end
end
end
sites.nr_per_name = nr_per_name;
%% Add the names of transition sites:
unique_names = unique(names);
% Single names:
for i = 1:size(unique(names))
trans_name = strcat('Transition:', unique_names(i));
a = length(sites.names) + 1;
sites.names(a) = trans_name;
end
% All combinations:
for i = 1:size(unique(names))
for j = 1:size(unique(names))
trans_name = strcat('Transition:', unique_names(i), '-', unique_names(j));
a = length(sites.names) + 1;
sites.names(a) = trans_name;
end
end
%% Loop over all possible positions, to see how close they are to each other
% Should be less as the 2*dist_close to prevent overlap, otherwise it will be adjusted
for site1 = 1:nr_sites-1
pos1 = sites_pos(:,site1);
for site2 = site1+1:nr_sites
pos2 = sites_pos(:,site2);
dist = sqrt(calc_dist_sqrd_frac(pos1, pos2, lattice));
if dist < 2*dist_close % Make dist_close smaller to prevent overlapping sites:
fprintf('WARNING! Crystallographic sites are overlapping with the chosen dist_close (%d), making dist_close smaller! \n', dist_close)
dist_close = 0.5*(dist - 0.01); % -0.01 for extra certainty
fprintf('dist_close changed to %d \n', dist_close)
if dist_close < 0.25 %Angstrom
disp('ERROR! ERROR! Two crystallographic sites are within half an Angstrom of each other')
disp('ERROR! This is NOT realistic, check/change the given crystallographic site locations in known_materials.m!')
fprintf('ERROR! Sites number %d and %d are overlapping \n', site1, site2)
fprintf('ERROR! Coordinates are: (%d, %d, %d) and (%d, %d, %d) \n', ...
pos1(1), pos1(2), pos1(3), pos2(1), pos2(2), pos2(3))
return
end
end
end
end
close_sqrd = dist_close*dist_close; % to avoid taking the root everytime
% Necessary matrices:
atom_site = zeros(sim_data.nr_steps, sim_data.nr_diffusing); % A value of zero means its not at a known site
site_occup = zeros(nr_sites,nr_parts);
transitions = zeros(nr_sites);
trans_counter = 0;
all_trans = zeros(1,4);
succes = zeros(nr_sites, nr_sites, nr_parts);
%% Loop over all times and diffusing atoms:
diff_atom = 0;
part = 1;
end_step = sim_data.nr_steps;
disp('Determining the site positions of all diffusing atoms:')
for atom = sim_data.start_diff_elem:sim_data.end_diff_elem
diff_atom = diff_atom + 1;
% First find the first site this atom at:
not_found = true;
site = 0;
time_start = 1;
while not_found && time_start < end_step
if site < nr_sites
site = site + 1;
else
time_start = time_start + 1;
site = 1;
end
pos_atom = sim_data.frac_pos(:,atom,time_start);
dist = calc_dist_sqrd_frac(pos_atom, sites_pos(:, site), lattice);
if dist < close_sqrd
not_found = false;
% Atom is at this site
atom_site(time_start, diff_atom) = site;
% Increase occupancy of this site
site_occup(site, part) = site_occup(site, part) + 1;
prev_site = site; % Update previous site for next transition
end
end
% After the first site has been found:
for time = time_start+1:sim_data.nr_steps
% Divide the simulation in multiple part's for statistics:
part = ceil(time/part_size);
pos_atom = sim_data.frac_pos(:,atom,time);
% First check the site found at the last time_step, since it's the most likely position:
dist = calc_dist_sqrd_frac(pos_atom, sites_pos(:, prev_site), lattice);
if dist < close_sqrd
atom_site(time, diff_atom) = prev_site; % Atom is at this site
site_occup(prev_site, part) = site_occup(prev_site, part) + 1; % Increase occupancy of this site
trans_start = time + 1;
else %not found at the previous position:
% Loop over all sites, exit the loop once a site has been found.
not_found = true;
site = 0;
while not_found && site < nr_sites
site = site + 1;
dist = calc_dist_sqrd_frac(pos_atom, sites_pos(:, site), lattice);
if dist < close_sqrd
not_found = false;
% Atom is at this site
atom_site(time, diff_atom) = site;
% Increase occupancy of this site
site_occup(site, part) = site_occup(site, part) + 1;
end
end
% If a transition happened, remember some things:
if ~not_found %Transition happened, i.e. the atom is found at a new site
transitions(prev_site, site) = transitions(prev_site, site) + 1;
trans_counter = trans_counter + 1;
all_trans(trans_counter,1) = diff_atom; %the atom
all_trans(trans_counter,2) = prev_site; % the starting position
all_trans(trans_counter,3) = site; %the final position
all_trans(trans_counter,4) = trans_start; %the time the atom left it's previouw position
all_trans(trans_counter,5) = time; % the timestep (at the end of the transition)
succes(prev_site, site, part) = succes(prev_site, site, part) + 1;
prev_site = site;
end
end
end
fprintf('*') % One atom has been finished...
end
fprintf('\n') % go to next line in output
fprintf('Number of transitions between sites: %d \n', trans_counter)
% Also assign the (failed) transition states with a site-number and name
diff_atom = 0;
time_start = 0;
nr_trans_found = 0;
for atom = sim_data.start_diff_elem:sim_data.end_diff_elem
diff_atom = diff_atom + 1;
prev_site = 0;
% What to do with atoms starting outside a given site?
for time = time_start+1:sim_data.nr_steps
if atom_site(time, diff_atom) == 0 && prev_site == 0
% It left the site
if time ~= 1
prev_site = atom_site(time - 1, diff_atom);
trans_start = time;
end
elseif atom_site(time,diff_atom) ~= 0 && prev_site ~= 0
% It appeared at a new site
now_site = atom_site(time, diff_atom);
trans_end = time-1;
% Determine the 'trans_site_nr'
if now_site == prev_site
trans_site_nr = nr_sites + now_site;
names(trans_site_nr) = strcat('Transition:', names(now_site));
atom_site(trans_start:trans_end, diff_atom) = trans_site_nr;
prev_site = 0;
else % loop over all transitions already found
if nr_trans_found == 0 % For the first transition site:
trans_found(1,1) = prev_site;
trans_found(1,2) = now_site;
nr_trans_found = 1;
trans_nr = 1;
else
trans_nr = 0;
% Check if it has been found earlier:
for i = 1:nr_trans_found
if trans_found(i,1) == prev_site && trans_found(i,2) == now_site
trans_nr = i;
end
end
if trans_nr == 0 % Not found in previous step, so it is a NEW one
nr_trans_found = nr_trans_found + 1;
trans_found(nr_trans_found,1) = prev_site;
trans_found(nr_trans_found,2) = now_site;
trans_nr = nr_trans_found;
end
end
% Assign trans_site_nr etc.:
trans_site_nr = 2*nr_sites + trans_nr;
atom_site(trans_start:trans_end, diff_atom) = trans_site_nr;
names(trans_site_nr) = strcat('Transition:', names(prev_site), '-', names(now_site));
prev_site = 0;
end
end
end
end
% Save things to sites:
sites.site_radius = dist_close;
sites.frac_pos = sites_pos;
sites.cart_pos = cart_pos;
sites.site_names = names;
sites.occupancy = sum(site_occup, 2);
sites.occup_parts = site_occup;
sites.atoms = atom_site;
sites.transitions = transitions;
sites.all_trans = all_trans;
sites.succes = succes;
sites.nr_parts = nr_parts;
finished = true;
end