Hello everyone,
I am a new user of LAMMPS. In the context of my PhD thesis, in a first time I want to reproduce the structure simulated by Vashishta et al. (Interaction potentials for alumina and molecular dynamics simulations of amorphous and liquid alumina, Journal of Applied Physics, 2008) of an amorphous alumina.
As I can’t attach files to my request, I will put some extracts from my codes.
Firstly, I have the potential file : AlO.vashishta
I am sure the file is well written, the units are in metals.
DATE: 26-03-2024 Benjamin
These entries are in LAMMPS “metal” units:
H = eV*Angstroms^eta; Zi, Zj = |e| (e = electronic charge);
lambda1, lambda4, rc, r0, gamma = Angstroms;
D = eVAngstroms^4; W = eVAngstroms^6; B = eV;
other quantities are unitless
element1 element2 element3
H eta Zi Zj lambda1 D lambda4
W rc B gamma r0 C cos(theta)
Al Al Al 12.7506 7.0 1.5237 1.5237 5.0 0.0 3.75
0.0 6.0 0.0 0.0 0.0 0.0 0.0
O O O 564.7334 7.0 -1.0158 -1.0158 5.0 44.5797 3.75
79.2884 6.0 0.0 0.0 0.0 0.0 0.0
O Al Al 249.3108 9.0 -1.0158 1.5237 5.0 50.1522 3.75
0.0 6.0 8.1149 1.0 2.90 10.0 -0.33331
Al O O 249.3108 9.0 1.5237 -1.0158 5.0 50.1522 3.75
0.0 6.0 12.4844 1.0 2.90 10.0 0.0
Al O Al 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Al Al O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O Al O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O O Al 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Then I run a Matlab code : INPUT_GENERATOR_Al2O3.m (code below)
%% INPUT GENERATOR
clc, clear all, hold off, close all, tic
% Size of the analysed box
nx=6;
ny=8;
nz=6;
% Lattice parameters
a=4.76;
b=a;
c=12.99;
Lx=nxa; % Angstrom
Ly=nybsqrt(3)/2.;
Lz=nzc;
Charges
qAl=1.5237; charge of Al [e]
qO=-1.0158; % charge of O [e]
atomtypes=2;
massAl = 26.981539; % atomok molaris tomege [g/mol]
massO = 15.999400;
% Unit coordinates
coorAl=[1 qAl 0.00000 0.00000 4.57572
1 qAl 0.00000 0.00000 1.92093
1 qAl 0.00000 0.00000 8.41758
1 qAl 0.00000 0.00000 11.07237
1 qAl 2.38010 1.37415 8.90682
1 qAl 2.38010 1.37415 6.25203
1 qAl 2.38010 1.37415 12.74868
1 qAl 2.38010 1.37415 2.41017
1 qAl 0.00000 2.74830 0.24462
1 qAl 0.00000 2.74830 10.58313
1 qAl 0.00000 2.74830 4.08648
1 qAl 0.00000 2.74830 6.74127];
coorO=[ 2 qO 1.45776 0.00000 3.24832
2 qO -0.72888 1.26246 3.24832
2 qO 1.65122 2.85999 3.24832
2 qO 3.30244 0.00000 9.74498
2 qO -1.65122 2.85999 9.74498
2 qO 0.72888 1.26246 9.74498
2 qO 3.83786 1.37415 7.57942
2 qO 1.65122 2.63661 7.57942
2 qO 1.65122 0.11169 7.57942
2 qO 0.92234 1.37415 1.08277
2 qO 3.10898 0.11169 1.08277
2 qO 3.10898 2.63661 1.08277
2 qO 1.45776 2.74830 11.91052
2 qO -0.72888 4.01076 11.91052
2 qO -0.72888 1.48584 11.91052
2 qO -1.45776 2.74830 5.41387
2 qO 0.72888 1.48584 5.41387
2 qO 0.72888 4.01076 5.41387];
figure()
set(gcf,‘color’,‘w’)
scatter3(coorAl(:,3), coorAl(:,4), coorAl(:,5),‘or’,‘markerfacecolor’,‘r’)
hold on
scatter3(coorO(:,3), coorO(:,4), coorO(:,5),‘ob’,‘markerfacecolor’,‘b’)
axis equal
axis([0 Lx 0 Ly 0 Lz])
xlabel(‘x’)
ylabel(‘y’)
zlabel(‘z’)
%COORD=zeros(nxnynz*30,3);
ATOMS=[];
for i=1:nx+1
for j=1:ny+1
for k=1:nz+1
ATOMS=[ATOMS; [coorAl;coorO]+repmat([0 0 (i-1)*a-(j-1)*a/2 (j-1)bsqrt(3)/2 (k-1)*c],[30 1])];
end
end
end
idAl = find(ATOMS(:,1)==1);
idO = find(ATOMS(:,1)==2);
figure()
set(gcf,‘color’,‘w’)
scatter3(ATOMS(idAl,3), ATOMS(idAl,4), ATOMS(idAl,5),‘or’,‘markerfacecolor’,‘r’)
hold on
scatter3(ATOMS(idO,3), ATOMS(idO,4), ATOMS(idO,5),‘ob’,‘markerfacecolor’,‘b’)
axis equal
%axis([0 Lx 0 Ly 0 Lz])
xlabel(‘x’)
ylabel(‘y’)
zlabel(‘z’)
for i=1:length(ATOMS(:,1))
dist=sqrt(sum((ATOMS([1:i-1,i+1:end],3:5)-repmat(ATOMS(i,3:5),[length(ATOMS(:,1))-1,1])).^2,2));
mdist(i)=min(dist);
end
min(mdist)
return
% Tidying the coordinates
ATOMS(ATOMS(:,3)<-1e-4,3)=ATOMS(ATOMS(:,3)<-1e-4,3)+Lx;
ATOMS(ATOMS(:,3)>Lx-1e-6,:)=[];
ATOMS(ATOMS(:,4)>Lx-1e-6,:)=[];
nAtom=nxnynz*30;
%scatter3(ATOMS(:,3), ATOMS(:,4), ATOMS(:,5))
idAl = find(ATOMS(:,1)==1);
idO = find(ATOMS(:,1)==2);
figure()
set(gcf,‘color’,‘w’)
scatter3(ATOMS(idAl,3), ATOMS(idAl,4), ATOMS(idAl,5),‘or’,‘markerfacecolor’,‘r’)
hold on
scatter3(ATOMS(idO,3), ATOMS(idO,4), ATOMS(idO,5),‘ob’,‘markerfacecolor’,‘b’)
axis equal
axis([0 Lx 0 Ly 0 Lz])
xlabel(‘x’)
ylabel(‘y’)
zlabel(‘z’)
VELO=zeros(nAtom,3);
% INPUT File Export
file_1 = fopen([‘Al2O3_input_n’ num2str(nx) num2str(ny) num2str(nz) ‘.dat’],‘wt’);
fprintf(file_1,‘LAMMPS data file from restart file: timestep = 0, procs = 1\n\n’);
fprintf(file_1,‘%1.0d’,nAtom);
fprintf(file_1,’ atoms’);
fprintf(file_1,‘\n’);
fprintf(file_1,‘\n’);
fprintf(file_1,‘%1.0g’,atomtypes);
fprintf(file_1,’ atom types’);
fprintf(file_1,‘\n’);
fprintf(file_1,‘\n’);
fprintf(file_1,‘%1.4g’,0);
fprintf(file_1,’ ‘);
fprintf(file_1,’%1.4g’,Lx);
fprintf(file_1,’ ');
fprintf(file_1,‘xlo xhi\n’);
fprintf(file_1,‘%1.4g’,0);
fprintf(file_1,’ ‘);
fprintf(file_1,’%1.4g’,Ly);
fprintf(file_1,’ ');
fprintf(file_1,‘ylo yhi\n’);
fprintf(file_1,‘%1.4g’,0);
fprintf(file_1,’ ‘);
fprintf(file_1,’%1.4g’,Lz);
fprintf(file_1,’ ');
fprintf(file_1,‘zlo zhi\n\n’);
fprintf(file_1,‘Atoms\n\n’);
for int=1:nAtom
fprintf(file_1,‘%1.0d’,int);
fprintf(file_1,’ ‘);
fprintf(file_1,’1.0d',ATOMS(int,1));
fprintf(file_1,' ');
fprintf(file_1,‘3.4g',ATOMS(int,2));
fprintf(file_1,’ ‘);
fprintf(file_1,’%3.6g’,ATOMS(int,3));
fprintf(file_1,’ ‘);
fprintf(file_1,’%3.6g’,ATOMS(int,4));
fprintf(file_1,’ ‘);
fprintf(file_1,’%3.6g’,ATOMS(int,5));
fprintf(file_1,‘\n’);
end
fprintf(file_1,‘\n\n’);
fprintf(file_1,‘Masses\n\n’);
fprintf(file_1,‘%1.0d’,1);
fprintf(file_1,’ ‘);
fprintf(file_1,’%3.6g’,massAl);
fprintf(file_1,‘\n’);
fprintf(file_1,‘%1.0d’,2);
fprintf(file_1,’ ‘);
fprintf(file_1,’%3.6g’,massO);
fprintf(file_1,‘\n\n’);
fprintf(file_1,‘Velocities\n\n’);
for int=1:nAtom
fprintf(file_1,‘%1.0d’,int);
fprintf(file_1,’ ‘);
fprintf(file_1,’%3.6g’,VELO(int,1));
fprintf(file_1,’ ‘);
fprintf(file_1,’%3.6g’,VELO(int,2));
fprintf(file_1,’ ‘);
fprintf(file_1,’%3.6g’,VELO(int,3));
fprintf(file_1,’ ‘);
fprintf(file_1,’\n’);
end
fclose(file_1);
This file is used to write the file Al2O3_input_n686.dat in order to describe the geometry of the system with 8640 atoms, as Vashishta et al. did. The beginning of this file is at follows :
LAMMPS data file from restart file: timestep = 0, procs = 1
8640 atoms
2 atom types
0 28.56 xlo xhi
0 32.98 ylo yhi
0 77.94 zlo zhi
Atoms
1 1 0 0 4.57572
2 1 0 0 1.92093
3 1 0 0 8.41758
4 1 0 0 11.0724
5 1 2.3801 1.37415 8.90682
6 1 2.3801 1.37415 6.25203
7 1 2.3801 1.37415 12.7487
8 1 2.3801 1.37415 2.41017
etc. (until atom 8640 + the velocities that are set to 0 in all dimensions)
Then I have the lammps file : in_NPT_Vashishta_Al2O3_Amorphous.test
This file was originally used for another type of system that was not described by Vashishta potential.
Sample generation of an amorphous alumina with a Vashishta potential parametrization
INPUT
Simulation size
variable nx equal 6 # Simulation cube size [A]
variable ny equal 8 # Simulation cube size [A]
variable nz equal 6 # Simulation cube size [A]
Temperature
variable Trun equal 1e-5 # Temperature de reference
variable Thigh equal 2330 # Temperature de chauffe
variable RateCool equal 10 # [K/ps] # vitesse de refroidissement
Cutoff
variable coff equal 6
Time control
variable dt equal 1e-3 # pas de temps
variable dumpf equal 50 #
variable runtime equal 100 # [ps]
variable runstep equal round({runtime}/{dt}) # nombre de pas de calcul
-------------------------------------------------------------------------
LAMMPS commands
-------------------------------------------------------------------------
Phisical environment
units metal #
atom_style atomic # atomic system with charges
boundary p p p # periodic along x, y, z
read_data Al2O3_input_n${nx}{ny}{nz}.dat # Specify atom coordinates
-------------------------- Vashishta potentials --------------------#
pair_style vashishta
pair_coeff * * AlO.vashishta Al O
Neighbor list building
neighbor 2 bin # skin (atom to store: cutoff+skin) style (bin or nsq), bin faster than nsq
neigh_modify every 1 delay 0 check yes one 1500 # building pairwise neighbor lists
update of neighboor list, every step, max 1500 neighboor atoms
----------------------- Variable defintion --------------------------------
Basic variables
variable s equal step # Timestep
variable tm equal time # Mechanical time
Environmental variables
variable T equal temp
variable at equal atoms
variable V equal vol
variable p equal press
variable dens equal density
Energy variables
variable Up equal pe
variable Uk equal ke
variable Utot equal etotal
variable Ental equal enthalpy
------------------------ Dump atomic positions -------------------------------------
thermo 500 # output every 500 timesteps
thermo_style custom step time temp atoms vol density press pe # variable to output
Affichage commande en ligne
fix Journal all print 500 "{s} {tm} {T} {at} {V} {dens} {p} {Uk} {Up} {Utot} {Ental}" &
file Journal_NPT_cf{coff}.dat screen no title &
“Timestep Time Temperature Number_of_atoms Volume Density Pressure Kinetic_energy Potential_energy Total_energy Enthalpy”
------------------------ Enviromental settings -------------------------------------
Environment (NVT, NPT, NVE)
Environmental settings
timestep ${dt}
reset_timestep 0
minimize 0.0 1.0e-8 10000 1000000 # energy minimization of the system
------------------------ Heating the sample at Thigh -------------------------------------
fix 1 all nvt temp {Thigh} {Thigh} 0.001 # NVT time integration Nose/Hoover
fix 2 all press/berendsen iso 1 1 2.0 modulus 35000 # pressure control by Berendsen barostat
run 2000 #number of steps
fix 1 all nvt temp {Thigh} {Thigh} 2 # NVT
run {runstep}
dump Coord_max all custom {dumpf} dump_max_NPT_cf${coff}.dat id type x y z # snapshot of quantities
run 1000
undump Coord_max
------------------------ Quenching ------------------------
variable run_1 equal round(({Thigh}-{Trun})/{RateCool}/{dt})
fix 1 all nvt temp {Thigh} {Trun} 2 # NVT
run ${run_1}
----------------------------- End of quenching -------------------------
fix 1 all nvt temp {Trun} {Trun} 2 # NVT
run {runstep}
dump Coord_0K all custom {dumpf} dump_0K_NPT_cf${coff}.dat id type x y z
run 1000
undump Coord_0K
minimize 0.0 1.0e-8 10000 1000000 # energy minimization of the system
write_restart Geom_0K_min_NPT_cf${coff}.res
unfix 1
unfix 2
----------------------------- Verification -------------------------
fix 1 all nve # NVE
run {runstep}
dump Coord_min all custom {dumpf} dump_minK_NPT_cf${coff}.dat id type x y z
run 1000
undump Coord_min
write_restart Geom_0K_ready_NPT_cf${coff}.res
As I try to run tha lammps file on the Command Prompt, it prints me this :
LAMMPS (2 Aug 2023 - Update 2)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
Loaded 1 plugins from C:\Users\guita\AppData\Local\LAMMPS 64-bit 2Aug2023\plugins
Reading data file …
orthogonal box = (0 0 0) to (28.56 32.98 77.94)
1 by 1 by 1 MPI processor grid
reading atoms …
8640 atoms
reading velocities …
8640 velocities
read_data CPU = 0.030 seconds
Reading vashishta potential file AlO.vashishta with DATE: 26-03-2024
Neighbor list info …
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 1500, page size: 100000
master list distance cutoff = 8
ghost atom cutoff = 8
binsize = 4, bins = 8 9 20
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair vashishta, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Setting up cg style minimization …
Unit style : metal
Current step : 0
Per MPI rank memory allocation (min/avg/max) = 17.39 | 17.39 | 17.39 Mbytes
** Step Time Temp Atoms Volume Density Press PotEng**
** 0 0 0 8640 73412.372 3.9942194 nan nan**
ERROR on proc 0: Non-numeric atom coords - simulation unstable (src/OPENMP/domain_omp.cpp:58)
Last command: minimize 0.0 1.0e-8 10000 1000000 # energy minimization of the system
With the error at the end and the ‘nan’ for pressure and energy of the system.
I have already made several researches to try to solve the problem.
I tried the delete delete_atoms overlap 0.1 all all command and it works as it delete 1440 atoms from the geometry but I don’t want to delete atoms to reproduce Vashishta’s conditions !
In the Matlab code I also used the following loop :
for i=1:length(ATOMS(:,1))
dist=sqrt(sum((ATOMS([1:i-1,i+1:end],3:5)-repmat(ATOMS(i,3:5),[length(ATOMS(:,1))-1,1])).^2,2));
mdist(i)=min(dist);
end
min(mdist)
The minimal distance found using this is 1.854 A, which corresponds exactly to the distance of the Al-O bonding.
I know it makes a long message but I can’t attach files to my request so I tried to give you all the elements necessary to see what could be wrong with the files that I use.
Thank you a lot for reading my request and trying to help me.
Benjamin