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

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

Answer 1

[darova]

I used polyxpoly to find intersection point dc

Then just scaled distance. The result:

See attached script

[Nismo NismoR]

I doubt that this would make the results to become a non-periodic one.

The current code resuls in the periodic structure

[darova]

Can you please make a simple drawing how the result should look like?

[Nismo NismoR]

For non periodic: the final pic is attached below. the structure doesn't fill the box, whereas for period one it does

[darova]

I thought you wanted a hole at the center

[Nismo NismoR]

One practice is creating a hole and other one is creating a non-periodic. The below picture is how it shoud look with a a hole.

[Nismo NismoR]

Delete some carbon atoms somewhere in the center zone, In another word, code a circular shape in the middle to delete inside the parameters and I dont know how.

It uses the same coding and modify it to get the shape

[darova]

inpolygon can select points inside circle

[Nismo NismoR]

So, I would define a circular shape and use inpolygon to select the points inside the shape and command to delete the points inside the circle?

Thank you

[Nismo NismoR]

Also, delete points inside a circle

[darova]

exactly

[darova]

or even simpler

% [x0,y0] center of a circle
% R radius of a circle
ix = (x-x0).^2 + (y-y0).^2 < R^2;       % indices inside a circle

[Nismo NismoR]

When I place/run the "ix = (x-x0).^2 + (y-y0).^2 < R^2;", it is not recognized as any functon or variable

I have define R=2;

[Nismo NismoR]

How would I be able to integrate or include it in the code, Thank you

[darova]

I don't understand. Please exaplain more

[Nismo NismoR]

how would I use "ix = (x-x0).^2 + (y-y0).^2 < R^2; " in my code. If you look at the code all the way in the begining, where would this "ix = (x-x0).^2 + (y-y0).^2 < R^2; " go and how can i make it to work?

[darova]

An example

x = store_Catoms(:,1);
y = store_Catoms(:,2);
R = 3;                          % circle radius
x0 = mean(x);                   % centered
y0 = mean(y);

t = linspace(0,360,20);         % data for circle
xb = R*cosd(t) + x0;
yb = R*sind(t) + y0;
ix = (x-x0).^2 + (y-y0).^2 < R^2;   % indices of points inside
plot(x,y,'.b')                  % plot all poinnts
hold on
plot(x(ix),y(ix),'or')          % plot points inside a circle
plot(xb,yb,'r')                 % plot a circle
hold off

[Nismo NismoR]

That sounds amazing. Thanks for your help. But how would you delete the points in the circle.

I can use these commands "x = store_Catoms(:,1);

y = store_Catoms(:,2);

R = 3; % circle radius

x0 = mean(x); % centered

y0 = mean(y);

t = linspace(0,360,20); % data for circle

xb = R*cosd(t) + x0;

yb = R*sind(t) + y0;

ix = (x-x0).^2 + (y-y0).^2 < R^2;"

And then what would be the function the delete the points inside the circle?

[darova]

Do you understand this line?

plot(x(ix),y(ix),'or')          % plot points inside a circle

[Nismo NismoR]

Yeah, Basically all the points plotted in the circle. or all the points within the parameter

[darova]

So why are asking how to delete points inside if you know how to select them?

[Nismo NismoR]

Because by using a delete comand with the selected points, It results to an error again and does not delete the points

delete(ix) or delete (plot(x(ix),y(ix),'or'))

Also, this is the first assignment that am doing in matlab, so learning all at the same time.

[darova]

You can assign points to another variable

x1 = x(~ix);                % select points outside circle

Or you can remove data from original array

x(ix) = [];                 % remove points inside a circle

[Nismo NismoR]

Now, I get this error ""Error using plot

Vectors must be the same length. when I use this;

x = store_Catoms(:,1);

y = store_Catoms(:,2);

R = 3; % circle radius

x0 = mean(x); % centered

y0 = mean(y);

t = linspace(0,360,20); % data for circle

xb = R*cosd(t) + x0;

yb = R*sind(t) + y0;

ix = (x-x0).^2 + (y-y0).^2 < R^2; % indices of points inside

x(ix) = []; % remove points inside a circle

plot(x,y,'.b') % plot all poinnts

hold on

plot(x(ix),y(ix),'or') % plot points inside a circle

plot(xb,yb,'r') % plot a circle

hold off

[Nismo NismoR]

This is the entire code again, and it did not make any changes:

% ================================

clear all

close all

clc

% ======================

% Define Graphene System

% ----------------------

% Define Mass :

% -------------

C_mass = 12.011;

% Define lattice constants of Graphene :

% ----------------------------------------------------------------

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

x = store_Catoms(:,1);

y = store_Catoms(:,2);

R = 3; % circle radius

x0 = mean(x); % centered

y0 = mean(y);

t = linspace(0,360,20); % data for circle

xb = R*cosd(t) + x0;

yb = R*sind(t) + y0;

ix = (x-x0).^2 + (y-y0).^2 < R^2; % indices of points inside

%x(ix) = []; % remove points inside a circle

%plot(x,y,'.b') % plot all poinnts

%hold on

%plot(x(ix),y(ix),'or') % plot points inside a circle

%plot(xb,yb,'r') % plot a circle

%hold off

% ===============================================================

% 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

% --------------------------------------------------------

[Nismo NismoR]

Or I get this picture and no changes to the Ovito shape or deleted points

[darova]

It means that arrays x and y are not the same size

Now, I get this error ""Error using plot
Vectors must be the same length. when I use this;

You are removing data from x

x(ix) = [];                 % remove points inside a circle

But you didn't remove data from y

[Nismo NismoR]

When I remove it from y, it give this command problem, Is "y(iy) = []; " is incorrect?

Enter the number of repeatation in X : 5

Enter the number of repeatation in Y : 5

Unrecognized function or variable 'iy'.

Error in X (line 89)

y(iy) = []; % remove points inside a circle

[darova]

I have nothing to say. Don't know how to talk with you

[Nismo NismoR]

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]

I am following every hint and when placing them in the Matlab, It does not produce the output.

[darova]

I don't want to do all your homework for you. You probably should start HERE

Answer 2

[Nismo NismoR]

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]

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]

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