Hello LAMMPS users,
I am modelling KMC diffusion in an Inconel-Nickel alloy structure, and I am employing PyLAMMPS to do it. These are the main steps involved in it:
- Going through the entire atom list and extracting their ID, type and position and storing them in arrays based on their atom IDs.
e = eval(“pe”) # Evaluating the energy of the initial structure
elast = e
current_atom = L.atoms[a]
IDa = current_atom.id
TYPEID[IDa] = current_atom.type # Extracting the type of the chosen atom
XCO[IDa], YCO[IDa], ZCO[IDa] = current_atom.position # Extracting the position of the chosen atom
LID[IDa] = a # Storing the PyLAMMPS-assigned ID to array relative to the atom ID
-
A random atom (ATOM_1) is chosen as the first atom to be switched, and its position is stored in x0, y0, z0
.
x0, y0, z0 = XCO[ATOM_1], YCO[ATOM_1], ZCO[ATOM_1] -
All the neighbours of ATOM_1 within the lattice distance are gathered and written into a text file.
L.region(“vic sphere”, x0, y0, z0, “${latt} side in units box”)
L.group(“vicatoms region vic”)
L.write_dump(“vicatoms custom Atoms2.txt id type x y z”)
- Open the ‘Atoms2.txt’ file, and based on a random number, one of the atoms from this file is chosen as the atom (ATOM_2) to be switched with ATOM_1. The position of ATOM_2 is stored in x1, y1, z1.
x1 = float(str(XCO[ATOM_2]).replace(“[”,“”).replace(“]”,“”))
y1 = float(str(YCO[ATOM_2]).replace(“[”,“”).replace(“]”,“”))
z1 = float(str(ZCO[ATOM_2]).replace(“[”,“”).replace(“]”,“”))
- ATOM_1 and ATOM_2 are switched by assigning the position of ATOM_1 to ATOM_2 and vice versa.
TEMP1 = int(str(LID[ATOM_1]).replace(“[”,“”).replace(“]”,“”)) # The LID array has the PyLAMMPS IDs of the atoms, which is different from the actual IDs of the atoms
current_atom_3 = L.atoms[TEMP1]
TEMP2 = int(str(LID[ATOM_2]).replace(“[”,“”).replace(“]”,“”))
current_atom_4 = L.atoms[TEMP2]
current_atom_3.position = x1, y1, z1
current_atom_4.position = x0, y0, z0
- The XCO, YCO and ZCO position array values of ATOM_1 and ATOM_2 are updated after the switching of atoms.
XCO[ATOM_1], YCO[ATOM_1], ZCO[ATOM_1] = x1, y1, z1
XCO[ATOM_2], YCO[ATOM_2], ZCO[ATOM_2] = x0, y0, z0
- Depending on the two types of atoms that are switched, the activation energy is calculated. One run is done with the switched atoms.
if atom1_type == 1 and atom2_type == 1:
E_a = 1.03369
elif atom1_type == 1 and atom2_type == 2:
E_a = 1.16538
elif atom1_type == 1 and atom2_type == 3:
E_a = 0.98484
elif atom1_type == 2 and atom2_type == 1:
E_a = 1.16538
elif atom1_type == 2 and atom2_type == 2:
E_a = 1.33551
elif atom1_type == 2 and atom2_type == 3:
E_a = 1.10366
elif atom1_type == 3 and atom2_type == 1:
E_a = 0.98484
elif atom1_type == 3 and atom2_type == 2:
E_a = 1.10366
else:
E_a = 0.9404
L.run(1, “pre no post no”) # One run with the switched atoms
- Energy is evaluated to a precision of 4.
prec = 4
e1 = L.eval(“pe”)
e = round(e1,prec)
- Metropolis criterion to determine whether the atom switch should be accepted or not.
metro_rand = np.random.random()
if e <= elast:
elast = e
elif metro_rand <= np.exp(numatoms * (elast - e)/kBT):
elast = e
else:
X1, Y1, Z1 = current_atom_3.position
X0, Y0, Z0 = current_atom_4.position
current_atom_3.position = X0, Y0, Z0 # ATOM_1 and ATOM_2 are switched back to their original position.
current_atom_4.position = X1, Y1, Z1
I get ‘nan’ values for the evaluated energy after a certain number of iterations every time I run the script. I checked the structure after every atom switch, and it remains stable. Could anyone suggest a potential solution to this problem? Thanks in advance.
Regards,
Rajesh