EAM implementation and dimension error

[kangpyo]

Hi, I've been working on integrating the JuLIP package with AtomsBase. While I understand that a more seamless integration might be on the horizon with the planned move to EmpiricalPotentials.jl, I've proceeded with the available conversion functionality to use JuLIP with AtomsBase structures.

Here's a snippet of my code, where the "Atoms" structure seems okay to set up the environment for the EAM potential:
`position_matrix = hcat([ [atom.position.x, atom.position.y, atom.position.z] for atom in lattice ]...)
cell_matrix = Diagonal(fill(aluminum_lattice_constant, 3))
atoms = Atoms(X=position_matrix, cell=cell_matrix, pbc=true)

aluminum_eam_file = "Al99.eam.alloy"
eam = pyimport("ase.calculators.eam")
eam_julip = EAM(aluminum_eam_file)
set_calculator!(atoms, eam_julip)
full code: [https://www.dropbox.com/scl/fi/lsbjq8duccjgs8ju9ahcz/dislocation_mac_EAM_01272024.jl?rlkey=m9he0ahhmef7w8hf15inty1t8&dl=0](url) EAM potential: [https://www.dropbox.com/scl/fi/x40qqcvsukjuequk16eu8/Al99.eam.alloy?rlkey=8is9lal1ymf27g1vokul5t915&dl=0](url) However, I've encountered a stumbling block when trying to run simulations.simulate!(system, simulator, number_of_steps)` I keep receiving a DimensionError stating: "8.878172622804642e-27 g¹ᐟ² nm ps⁻¹ u⁻¹ᐟ² and 0.0 kJ g⁻¹ nm⁻¹ are not dimensionally compatible." This error emerges despite my attempts to align the dimensions and units.

I'm uncertain if this is an issue with the EAM potential implementation or something else. I've been trying to troubleshoot this without success. Any insights or suggestions would be greatly appreciated

Thank you
best
KP

[jgreener64]

Could you paste the full code here?

[kangpyo]

#this code has error of the dimension mismatch
cd(@__DIR__)
using Pkg
Pkg.activate(".")
using JuLIP, InteratomicPotentials, ASE, PyCall
using GLMakie
using Molly, LinearAlgebra

struct Point
    x::Float64
    y::Float64
    z::Float64
end

struct Atom
    position::Point
end

function create_fcc_lattice(size, lattice_constant)
    lattice = Atom[]
    for i in 1:size
        for j in 1:size
            for k in 1:size
                push!(lattice, Atom(Point(i*lattice_constant, j*lattice_constant, k*lattice_constant)))
            end
        end
    end
    return lattice
end

function plot_lattice_makie(lattice)
    x = [atom.position.x for atom in lattice]
    y = [atom.position.y for atom in lattice]
    z = [atom.position.z for atom in lattice]
    fig = Figure()
    ax = Axis3(fig[1, 1])
    scatter!(ax, x, y, z, markersize=5, color=:blue)
    ax.title = "Aluminum Lattice Structure"
    ax.xlabel = "X"
    ax.ylabel = "Y"
    ax.zlabel = "Z"
    return fig
end

# Define the lattice constant for aluminum with a different variable name
aluminum_lattice_constant = 4.05e-10  # in meters

# Now call create_fcc_lattice with the lattice size and the new lattice constant variable
lattice_size = 10
lattice = create_fcc_lattice(lattice_size, aluminum_lattice_constant)


position_matrix = hcat([ [atom.position.x, atom.position.y, atom.position.z] for atom in lattice ]...)
cell_matrix = Diagonal(fill(aluminum_lattice_constant, 3))
atoms = Atoms(X=position_matrix, cell=cell_matrix, pbc=true)

aluminum_eam_file = "Al99.eam.alloy"
eam = pyimport("ase.calculators.eam")
eam_julip = EAM(aluminum_eam_file)
set_calculator!(atoms, eam_julip)

timestep = 1e-14  # Example timestep
number_of_steps = 1000  # Example number of steps
simulator = VelocityVerlet(dt=timestep)

# You can adjust this value based on the actual element or compound you are simulating
default_mass = 26.9815u"u"  # Atomic mass unit for aluminum

# Set temperature with appropriate units
temperatures = 300u"K"  # Example temperature

# Initialize velocities with units, using a different variable name to avoid conflict
my_velocities = [random_velocity(default_mass, temperatures) for _ in 1:length(lattice)]

# Define the size of your simulation cell
# Assuming a cubic cell with the size determined by the lattice constant and lattice size
# cell_size = lattice_size * aluminum_lattice_constant
# Convert the lattice_constant to nanometers and then multiply by the lattice size
cell_size_nm = lattice_size * aluminum_lattice_constant * 1e9u"nm"

# Create a cubic boundary condition with the specified size in nanometers
boundary_condition_with_units = CubicBoundary(cell_size_nm)

using StaticArrays  # Make sure to import StaticArrays

# Use molar mass for aluminum in grams per mole
default_mass = 26.9815u"g/mol"  # Molar mass of aluminum

# Assuming default values for atom properties
default_charge = 0.0
default_sigma = 1.0u"nm"  # Adjust as necessary, in nanometers
default_epsilon = 1.0u"kJ/mol"  # Adjust as necessary, in kJ/mol

# Create an array of Molly.Atom objects
molly_atoms = [Molly.Atom(index=i, 
                          charge=default_charge, 
                          mass=default_mass, 
                          σ=default_sigma, 
                          ϵ=default_epsilon,
                          solute=false) for i in 1:length(lattice)]

# Update velocities to be consistent with units
using Unitful  # Ensure Unitful is being used

# Convert the coordinates to have units, assuming they are in nanometers
coords_with_units = [SVector(atom.position.x, atom.position.y, atom.position.z) * 1u"nm" for atom in lattice]

# Create the Molly.System object with units-consistent boundary condition
system = Molly.System(
    atoms=molly_atoms, 
    coords=coords_with_units,  # These are already in nm
    velocities=my_velocities,  # These should be in nm/ps
    boundary=boundary_condition_with_units, 
    # ... other necessary properties ...
)
simulate!(system, simulator, number_of_steps)
# Example: Check units for velocities
# Assume default_mass is in atomic mass units (u) and temperatures in Kelvin (K)
default_mass = 26.9815u"u"  # Atomic mass unit for aluminum
temperatures = 300u"K"      # Temperature in Kelvin

# Calculate velocities in nm/ps
my_velocities = [random_velocity(default_mass, temperatures) for _ in 1:length(lattice)]

# Example: Check units for boundary conditions
# If coords are in nm, boundary should also be specified in nm
cell_size_nm = lattice_size * aluminum_lattice_constant * 1e9u"nm"
boundary_condition = CubicBoundary(cell_size_nm)

# Create Molly.System with consistent units
system = Molly.System(
    atoms=molly_atoms, 
    coords=coords_with_units,  # in nm
    velocities=my_velocities,  # in nm/ps
    boundary=boundary_condition,
    # ... other necessary properties ...
)

simulate!(system, simulator, number_of_steps)

plot_lattice_makie(lattice)

[kangpyo]

This is the EAM potential. I cannot attach the file directly, so I'm providing a link.
https://www.dropbox.com/scl/fi/x40qqcvsukjuequk16eu8/Al99.eam.alloy?rlkey=8is9lal1ymf27g1vokul5t915&dl=0

[jgreener64]

timestep = 1e-14  # Example timestep

should have units, e.g.

timestep = 1e-14 * u"s"  # Example timestep

That makes simulate! work. It won't do much interesting though, since there are no interactions given to System.

It also seems like system.coords has values differing by 1E-10 nm, which is probably not what you intended.

Hope that helps!

[kangpyo]

great, it works. thank you for your help. I fixed the system.coords. it seems there is another unit error during my tweaking.