juliamolsim / julip.jl Goto Github PK

The current TightBinding.jl crashes with the efficient "no-copy" implementation of positions. No clue whatsoever what is going on. I reverted to the old implementation which makes a double (triple?) copy. This works for now, but this needs a fix.

The tag name "v.0.1.0" is not of the appropriate SemVer form (vX.Y.Z).
cc: @cortner

Either

• Implement minimalistic visualisation using ThreeJS and/or GLVisualise

or

• get chemview to work under Julia

The following seems to work for the Travis-CI tests, but may not be optimal

ENV["PYTHON"]=""
Pkg.add("PyCall")
Pkg.add("Conda")
using Conda
Conda.add("NumPy")
Conda.add("SciPy")
Conda.add("pip")
pip = joinpath(Conda.BINDIR, "pip")
run(`$pip install ase`)
run(`$pip install -e git+https://github.com/libAtoms/matscipy#egg=matscipy`)
Pkg.clone(pwd())
Pkg.build("JuLIP")
Pkg.test("JuLIP"; coverage=true)

@jameskermode wrote in another issue:

@cortner the low-level neighbour_list C call looks like this:

_matscipy.neighbour_list(quantities, a.cell,
                                    np.linalg.inv(a.cell.T), a.pbc,
                                    a.positions, cutoff, *args)

so it can already be called without an ASE Atoms object, but doing it without Python at all is much harder, as the Python API is very convenient for building up variable length return arrays etc.

introduce temporary arrays for Atoms, and implement more clever calculators that precompute and store information to save compute time, e.g., any PairPotentialCalculator should compute energy and forces as the same time and store in an array that gets wiped after each positions update

Hi Christoph,

When I try testsolve.jl (with include("testsolve.jl"), I got the following error message,

same test but large and with Exp preconditioner

ERROR: LoadError: MethodError: no method matching PyAMG.AMGSolver{PyAMG.RugeStuben}(::SparseMatrixCSC{Float64,Int64}; tol=1.0e-7)

I tried to use "using PyAMG" first also with include("runtests.jl"), but those gave me the same error.

Best regards,
Lei

Why is this happening????? Maybe cell means something different from what I thought?

using JuLIP
using JuLIP.Potentials

println("-------------------------------------------------")
println("   Variable Cell Test")
println("-------------------------------------------------")
calc = lennardjones(r0=JuLIP.ASE.rnn("Al"))
at = Atoms("Al", pbc=(true,true,true)) * 2   # cubic=true,
set_calculator!(at, calc)
set_constraint!(at, VariableCell(at))

# compute energy before rotation
E0 = energy(at)

# generate a rotation matrix
U = [0 1.0 0; -1.0 0 0; 0 0 0]
Q = expm(0.1 * U)
@assert abs(det(Q) - 1.0) < 1e-10
@assert vecnorm(Q'*Q - eye(3)) < 1e-10

# deform cell
F0 = cell(at)
X0 = positions(at)
F = Q * F0
X = [ JMat(Q) * x for x in X0 ]
set_cell!(at, F)
set_positions!(at, X)

# compute energy after rotation
Eh = energy(at)
@show E0, Eh
# output
#    (E0,Eh) = (-126.23143471478076,-112.77456617599405)

I am thinking of slowly start transitioning to Julia v0.6 in preparation for v0.7 = v1.0. Apparently jumping a version is not advised. Any concerns?

title says it. Link to NearestNeighours.jl

James brought up this question when I started work on SaddleSearch.jl: what should be in JuLIP and what shouldn't be? More generally, should JuLIP become a complete environment, or do we want to follow the Julia convention of many smaller packages.

using @GRAD as opposed to grad has a huge performance penalty (type instability? check this) even though they should be equivalent.

Collect some thoughts here before going ahead . . .

CC @jameskermode

I thought, we'd collect some ideas and action points to discuss how to and the consequences of a pure Julia Atoms object.

TODO

• pure Julia neighbourlist (EDIT)

ReverseDiffSource will never be able to resolve the complex layers of function calls, but in particular the neighbourlist. The most likely candidates for more general usage of AD are ReverseDiffPrototype.jl and Autodiff.jl. But those (as well as ForwardDiff) require generic types instead of Float64 being hard-coded. So to take advantage of AD packages some fairly serious rewriting is needed unfortunately.

ReverseDiffSource.jl now supports variable length inputs.

• fix JuLIP to be able to ReverseDiffSource differentiate all potentials
• test performance of automatic versus manual
• probably rewrite all potentials to use ReverseDiffSource.

@jameskermode

FYI, this is now implemented, about to be merged, let me know if there is a bug.

julia> using JuLIP
WARNING: 'gpaw' Python module cannot be imported: JuLIP.DFT disabled

julia> at = bulk("Si", cubic=true)

julia> x, y, z = xyz(at)
([0.0,1.3575,0.0,1.3575,2.715,4.0725,2.715,4.0725],[0.0,1.3575,2.715,4.0725,0.0,1.3575,2.715,4.0725],[0.0,1.3575,2.715,4.0725,2.715,4.0725,0.0,1.3575])

julia> x
8-element SubArray{Float64,1,Array{Float64,2},Tuple{Int64,Colon},true}:
 0.0
 1.3575
 0.0
 1.3575
 2.715
 4.0725
 2.715
 4.0725

and then

julia> x[1] = 0.1
0.1

julia> set_positions!(at, x, y, z)

julia> positions(at)[1]
3-element StaticArrays.SVector{3,Float64}:
 0.1
 0.0
 0.0

to make it easy to implement much more general constraints, Including bound constraints, box constraints, all sorts of useful stuff.

WARNING: adsite.jl could not be included; most likely some AD package is missing;
at the moment it needs ForwardDiff, ReverseDiffPrototype
In [3]:

Probably as PairPotential(ex::Expr).

Cf. Vectorize.jl

It would be a shame though to go back to the old days of vectorizing everything.

It occurred to me that #20 would also lead to a natural way to incorporate dynamics into the framework. The idea would be to define (e.g.) a Dynamics <: DofManager (ignoring cell shape dynamics for now) which returns

  dofs(at) = [positions[:]; momenta[:]] 

as opposed to what happens at the moment

 dofs(at) = positions[:]

is already on 15 December. I am inclined to immediately move to v0.7 = v1.0, then I won't develop much (at all) for v0.6. I will just test whether the main dependencies are immediately being updated to work on v0.7

create a python module that implements a JuliaCalculator or JulipCalculator from ASE

StaticArrays has the advantage that arrays are subtypes of AbstractArray, which has a number of advantages. Consider switching. But this shows also how crucial it will be to keep the usage of fixedsize arrays in a way that allows switching of packages, since it is not clear what the future of each of these packages is.

Right now, I've written energy as

   energy <- energy - pressure * det( defm(at) ) 

which is surely wrong, since there should be some form of normalisation against size of the of the cell; deformed / undeformed maybe? But what does this mean here?

I looked in ASE, but there doesn't seem to be any pressure option in the UnitCellFilter?

CC @jameskermode

. . . in JuLIP.Visualise

Could be as simple as this:

import JuLIP.ASE

function read(filename)
    atoms = JuLIP.ASE.ase_io.read(filename)
    return JuLIP.ASE.ASEAtoms(atoms)
end

But perhaps name shouldn't clash with builtin read()?

There is an interesting new "feature" in v0.6 : "In version 0.6, Julia now keeps track of all the callers of a given function (called backedges) so it can go back and invalidate the previous code with any updated methods."

See e.g.

• https://discourse.julialang.org/t/running-in-world-age-x-while-current-world-is-y-errors/5871
• https://discourse.julialang.org/t/applicable-method-may-be-too-new-error/2733
• https://discourse.julialang.org/t/help-with-solving-a-world-age-error-in-package-on-julia-master/1126

The is explained in detail in Redefining Methods

Related from #68 :
preconditioner matrix assembly in v0.6 crashes; cf. preconditioners.jl / testsolve.jl; I am now guessing that the problem has something to do with either the precomputation of the refactor and then redifinition; or with the @d, @dd macros. Will need to revisit them or get rid of them. Or quite possible there is a problem with the way the analytic potentials are defined. . .

I've run into a case where I construct the neighbourlist first with a cutoff r1 and then with a different cutoff r2, but on the same set of positions. In that case we don't recompute. The canonical solution is probably to store the neighbour lists with keys of the form (:nlist, rcut).

This is just a reminder to go back and try JuLIP with the augmented Lagrangian implementation of constrained optimisation from https://github.com/cortner/ConstrainedOptim.jl

Related to issue #33

The distinction between calculators and Potentials looks more and more arbitrary. E.g., we have a PairCalculator which has as a single argument the pair potential.

Probably want to refactor code a little so that a Potential is a calculator and then use dispatch to choose the correct versions of energy, force, etc. E.g.

abstract PairPotential <: AbstractCalculator
type AnalyticPairPotential <: PairPotential  ... end 
function energy(pot::PairPotential, at::AbstractAtoms)
   # requires only that `pot` has expected pair potential behaviour
end

CC @jameskermode

using JuLIP, JuLIP.Potentials, PyPlot
V = PairPotential(:( (r-1.0)^2 - 1.0 )) * StepFunction(2.0)
r = linspace(0, 3, 101)
plot(r, V.(r), "b-")
axis([0, 3, -1.1, 1.1])

S Johnson pointed out that EAM can simply be written as f2(sum(f1, r)). What else can be made this short?

The aim is to define a distance between images along the string at the net of “jumps" across the PBC as it evolves.

So, the issue is that I have two configurations X1, X2, and I want to compute the difference X2 - X1 but need to do it in a continuous way.

One atom of X1 may cross the unit cell boundary during dynamics/relaxation without the same atom crossing in X2.

computing min configuration distance:

def undo_pbc_jumps(at):
    if not at.has('prev_pos'):
        at.set_array('prev_pos', at.get_positions())
    prev_pos = at.get_array('prev_pos')
    p = at.get_positions().T
    g = np.linalg.inv(at.get_cell().T)
    p -= np.dot(at.cell.T, np.floor(np.dot(g, (p - prev_pos.T))+0.5))        
    at.positions[:] = p.T
    at.set_array('prev_pos', p.T)
 
def pbc_diff(a, b):
    aa = a.copy()
    undo_pbc_jumps(aa)
    bb = aa.copy()
    bb.set_positions(b.positions)
    undo_pbc_jumps(bb)
    return aa.get_positions() - bb.get_positions()

It looks as if, at the moment, we make two copies of a positions (or similar) array when converting between ASE and JuLIP. First copying from Julia to Python and then transposing; but the memory layout of the original and final are the same, so it should be possible to convert between the JuLIP and ASE arrays and never copy them.

I would like to be able to create multiple subtypes of ASE Calculator, e.g. to allow dispatching differently for DFT-based calculators.

How about renaming the current ASECalculator to AbstractASECalculator and adding the following:

ASECalculator <: AbstractASECalculator <: AbstractCalculator
GPAWCalculator <: AbstractASECalculator <: AbstractCalculator

The tag name "v0.1" is not of the appropriate SemVer form (vX.Y.Z).
cc: @cortner

should be generalised to compute a reasonable rnn estimate for an arbitrary set of positions

cf Logging.jl

I would like to be able to store additional per-atom scalar and vector and tensor properties inside the Atoms type. Examples include local energy, masses, atomic numbers, momenta and/or velocities, local stresses.

Would an at.arrays::Dict{String,Any} that maps names to values be a reasonable solution? Should 1xN, 3xN and 3x3xN properties live in different places?

The following script demonstrates the error. We tracked it down as far as set_position_dofs!().

atoms_broken_relaxed.xyz.txt

using PyCall
using JuLIP
@pyimport atomistica
@pyimport ase

atoms = ASEAtoms(ASE.ase_io.read("atoms_broken_relaxed.xyz"))
calc = ASECalculator(atomistica.TersoffScr())
set_calculator!(atoms, calc)

set_constraint!(atoms, FixedCell(atoms))
println(energy(atoms))
println(energy(atoms, dofs(atoms)))

Results are:

julia> include("run.jl")
-2040.7844479506173
NaN

i.e. second call to energy(atoms, dofs(atoms)) gives a NaN result.

Am I doing something obviously stupid here? I'm trying to plot the dimer energy curve for LennardJones but get ~zero for all bond lengths. Same code works OK with StillingerWeber (after adjusting linspace range).

using JuLIP
using JuLIP.Potentials

si = ASEAtoms("Si2")
set_cell!(si, 100.0*eye(3))
set_pbc!(si, true)

p = positions(si) |> mat
sw = LennardJones(1.0, 1.0)
#sw = StillingerWeber()
set_calculator!(si, sw)
forces(si)

r = linspace(0.5, 1.5, 10)
e = []
for d in r
    p[1,1] = d
    set_positions!(si, p)
    push!(e, energy(si))
end
(r, e)

Output:

(linspace(0.5,1.5,10),Any[-3.02097e-69,-3.02097e-69,-3.02097e-69,-3.02097e-69,-3.02097e-69,-3.02097e-69,-3.02097e-69,-3.02097e-69,-3.02097e-69,-3.02097e-69])

This issue is to discuss anything related to the following variable cell implementation:

We transform F = exp(U), then the following happens:

• F is a rotation iff. U = - U’ (e.g. one can get rid of rotational motion by just imposing the U = U’)
• F is diagonal iff U is diagonal (for isotropic constraint)
• det(F) = exp( trace(U) ) (so volume constraint becomes linear, trace(U) = log(volume); so trivial to impose in optimisation)
• maybe the best one: F is non-degenerate (det(F) \neq 0) for ALL U (so no increment can degenerate the cell)

Unfortunately, in initial tests this turned out to be not robust.

Since making a Potential a calculator, JuLIP got confused about what grad means. It can mean either - forces or it can mean evaluate_d * R/r. I suggest to rename -forces = gradient and keep grad = evaluate_d * R/r (or similar usage).

You are overriding show(x) to call display(x), which isn't how these functions are intended to be customized. See the manual on custom pretty printing.

To customize REPL output, you should override Base.show(io::IO, x::MyType) and/or Base.show(io::IO, ::MIME"text/plain", x::MyType), as explained in the manual.

main barrier to doing it is to check that none of the implemented potentials assume that they only get data within the cut-off

for harmonic potentials and similar

. . . in particular without needing to reference to Optim.jl