# How to modify the code to create a Non-periodic structure

29 views (last 30 days)
Nismo NismoR on 14 Feb 2020
Commented: Nismo NismoR on 15 Feb 2020
The current code exhibits the periodic structure. The goal is to create a non-periodic structure.
close all
clc
C_mass = 12.011;
latt_c = 1.40; % length between nearest C-C
% Define length of unit cell of graphene :
ux = latt_c*cos(pi/3)+latt_c+latt_c*cos(pi/3)+latt_c;
uy = latt_c*sin(pi/3)+latt_c*sin(pi/3);
% Define coordinates of the unit cell:
u = latt_c*[0.0 0.0
cos(pi/3) sin(pi/3)
1+cos(pi/3) sin(pi/3)
1+2*cos(pi/3) 0.0];
no_atom_unit = length(u(:,1));
% =========================================================================
% ---------------------------------------------------------------------
% Define number of repeatation in X and Y direction to generate Graphene
% Structure :
% ----------------------------------------------------------------------
% Repeatation in X direction for Left Side
num_xdir = input('Enter the number of repeatation in X : ') ;
% % Repeatation in Y direction
num_ydir = input('Enter the number of repeatation in Y : ') ;
% ----------------------------------------------------------------------
% Total number of atoms in graphene
tot_atoms = num_xdir*num_ydir*no_atom_unit;
% Allocate a matrix with all zero entries to store carbon atoms for
% graphene
store_Catoms = zeros(tot_atoms,3);
% ================== Length of the Graphene ===============
tot_Lx = ux*num_xdir; % Total Length in X Direction
tot_Ly = uy*num_ydir; % Total Length in Y Direction
% -------------------------------------------
% Define Graphene Coordinates
% -------------------------------------------
C_ID = 0.0; % Initialize the Carbon Atom
cord_z = 1 ; % All carbon atoms have same z coordinate. Define the z coord.
for j = 1:num_ydir
for i = 1:num_xdir
for c = 1:no_atom_unit
C_ID = C_ID+1; % Update the coordinate of carbon atom
store_Catoms(C_ID,1) = u(c,1)+(i-1)*ux;
store_Catoms(C_ID,2) = u(c,2)+(j-1)*uy;
store_Catoms(C_ID,3) = cord_z;
end
end
end
% ===============================================================
% Write in the format for Molecular Dynamics (MD) simulation
% ===============================================================
filename = ['data.graphene_MD'];
fid = fopen(filename,'wt');
fprintf(fid,'%s\n','');
fprintf(fid,'\t\t%4i %s\n',tot_atoms,'atoms');
fprintf(fid,'%s\n','');
fprintf(fid,'\t\t%s\n','1 atom types');
fprintf(fid,'%s\n','');
% -----------------------------------
x_lo = min(store_Catoms(:,1));
x_hi = x_lo + tot_Lx;
y_lo = min(store_Catoms(:,2));
y_hi = y_lo + tot_Ly;
z_lo = cord_z - 4;
z_hi = cord_z + 4;
% -----------------------------------
fprintf(fid,'\t%12.8f %12.8f \t%s\n',x_lo,x_hi,'xlo xhi');
fprintf(fid,'\t%12.8f %12.8f \t%s\n',y_lo,y_hi,'ylo yhi');
fprintf(fid,'\t%12.8f %12.8f \t%s\n',z_lo,z_hi,'zlo zhi');
fprintf(fid,'%s\n','');
fprintf(fid,'%s\n','Masses');
fprintf(fid,'%s\n','');
fprintf(fid,'\t\t%s %12.8f\n','1',C_mass);
fprintf(fid,'%s\n','');
fprintf(fid,'%s\n','Atoms');
fprintf(fid,'%s\n','');
% --------------------------------------------------------
% Define the position of Carbon :
% -------------------------------
num = 1;
for i = 1:tot_atoms
fprintf(fid,'\t %4i %4i %12.8f %12.8f %12.8f\n',num,1,store_Catoms(i,1),store_Catoms(i,2),store_Catoms(i,3));
num = num + 1;
end
% --------------------------------------------------------

darova on 14 Feb 2020
I used polyxpoly to find intersection point dc
Then just scaled distance. The result:
See attached script

Nismo NismoR on 15 Feb 2020
Like I said, I am doing this for the first assignemnt. And, clearly you know alot better. It would have been easier if this would work rather than just exchanging many comments. because not only I did not learn, but also being more confused and feel I wasted both of our times.
I clearly have no idea and little hints does not work. Thank you
Nismo NismoR on 15 Feb 2020
I am following every hint and when placing them in the Matlab, It does not produce the output.
darova on 15 Feb 2020
I don't want to do all your homework for you. You probably should start HERE

Nismo NismoR on 15 Feb 2020
Thanks for the link. I have done most of the work including other parts that are not mentioned here. It's just that portion that drives me crazy

darova on 15 Feb 2020
try
ix = (x-x0).^2 + (y-y0).^2 < R^2; % indices of points inside
x(ix) = []; % remove points inside a circle
y(ix) = []; % remove points inside a circle
Nismo NismoR on 15 Feb 2020
I see now. I have been defining (or trying to remove the points in y-direction) it wrong. I tried it and it worked. Thank you.
Now, am not sure why this does not work in Ovito
One other way that i tried is:
C = [0.50*(x_lo+x_hi),0.50*(y_lo+y_hi)];
C_x = C(1);
C_y = C(2);
L = x_hi - x_lo;
R = 0.15*L;
New_DefCatoms = [];
ID_new = 0;
for i = 1:tot_atoms
i_x = store_atoms(i,1);
i_y = store_atoms(i,2);
i_z = store_atoms(i,3);
distance = sqrt((i_x-C_x)^2+(i_y-C_y)^2);
end