In [ ]:

    

import Optim
using LujiaLt
using LujiaLt.Potentials, LujiaLt.Plotting, LujiaLt.Solve


    

In [ ]:

    

# vacancy defect: create geometry and visualise
at = Atm(V=ToyEAMPotential(), Ra=8.1, defect=:vacancy)
plot(at, axis = [-5.1, 5.1, -3.1, 3.1])
# optimise the geometry; this uses Precon.jl (Hager Zhang CG) with a preconditioner to speed up convergence
Ysol = solve(at, display_result = true)
# plot the equilibrated positions
plot(at, X=Ysol, axis = [-5.1, 5.1, -3.1, 3.1]);


    

In [ ]:

    

# same for interstitial defect: 
# ToyEAM seems to be unstable here (?!?), so we use LJ instead
import Optim
at = Atm(V=LennardJonesPotential(), Ra=8.1, defect=:interstitial)
plot(at, axis = [-5.1, 5.1, -3.1, 3.1], rnn = 1.1)
Ysol = solve(at, display_result = true, Optimiser=Optim.ConjugateGradient, randomise=0.01)
plot(at, X=Ysol, axis = [-5.1, 5.1, -3.1, 3.1], rnn = 1.3);


    

In [ ]:

    

# for fun - try a different initial condition
at = Atm(V=LennardJonesPotential(), Ra=8.1, defect=:interstitial)
Y = at.geom.tri.X
I0 = size(Y, 2)  # interstitial site
I1 = find( (abs(Y[2,:]) .< 0.1) .* (abs(Y[1,:]) .< 0.1) )[1]  #  left neighbour
I2 = find( (abs(Y[2,:]) .< 0.1) .* (abs(Y[1,:]-1.0) .< 0.1) )[1]  # right neighbour
I3 = find( (abs(Y[1,:]-0.5) .< 0.1) .* (abs(Y[2,:]+sqrt(3)/2) .< 0.1) ) # bottom neighbour
Y[2,I0] -= sqrt(3)/4
Y[1,I1] -= 0.1
Y[2,I2] += 0.1
Y[2,I3] -= 0.2
at.geom.tri.X = Y
plot(at, axis = [-5.1, 5.1, -3.1, 3.1], rnn=1.1)
Ysol = solve(at, display_result = true, Optimiser=Optim.ConjugateGradient, show_trace=false)
plot(at, X=Ysol, axis = [-5.1, 5.1, -3.1, 3.1], rnn=1.3)


    

In [ ]:

    

}