EAM implementation and dimension error about molly.jl

[kangpyo]

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

[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!