MATLAB Answers

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

29 views (last 30 days)
Nismo NismoR
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','LAMMPS readable file ');
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
% --------------------------------------------------------

  0 Comments

Sign in to comment.

Accepted Answer

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

  31 Comments

Nismo NismoR
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
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
darova on 15 Feb 2020
I don't want to do all your homework for you. You probably should start HERE

Sign in to comment.

More Answers (1)

Nismo NismoR
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

  2 Comments

darova
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
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

Sign in to comment.

Sign in to answer this question.