How to modify the code to create a Non-periodic structure
4 views (last 30 days)
Show older comments
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
Accepted Answer
darova
on 14 Feb 2020
I used polyxpoly to find intersection point dc
Then just scaled distance. The result:
See attached script
31 Comments
More Answers (1)
Nismo NismoR
on 15 Feb 2020
2 Comments
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
See Also
Categories
Find more on Surfaces, Volumes, and Polygons in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!