Chemfiles contains position function, but in the example, there is a line like:
positions = positions(frame):

using Chemfiles

trajectory = Trajectory("filename.xyz")
frame = read(trajectory)

println("There are $(size(frame)) atoms in the frame")
positions = positions(frame)

# Do awesome things with the positions here !

The line conflicts the function name and the variable name:

julia> positions = positions(frame)
ERROR: cannot assign a value to variable Chemfiles.positions from module Main
Stacktrace:
 [1] top-level scope
   @ REPL[11]:1
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[7]:1

So, might be better to modify the line.

I report what I encountered.

I did the following, but I got errors for the test.

My environment

julia> versioninfo()
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin21.2.0)
  CPU: Apple M1
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, cyclone)

Error message

(@v1.7) pkg> add Chemfiles
    Updating registry at `~/.julia/registries/General.toml`
   Resolving package versions...
    Updating `~/.julia/environments/v1.7/Project.toml`
  [46823bd8] + Chemfiles v0.10.2
  No Changes to `~/.julia/environments/v1.7/Manifest.toml`

(@v1.7) pkg> test Chemfiles
     Testing Chemfiles
      Status `/private/var/folders/sb/1f2l11fn2bx_21yxlslz9kqc0000gn/T/jl_Zjx62t/Project.toml`
  [b99e7846] BinaryProvider v0.5.10
  [46823bd8] Chemfiles v0.10.2
  [34da2185] Compat v3.41.0
  [ffbed154] DocStringExtensions v0.8.6
  [78a364fa] Chemfiles_jll v0.10.2+0
  [8f399da3] Libdl `@stdlib/Libdl`
  [8dfed614] Test `@stdlib/Test`
      Status `/private/var/folders/sb/1f2l11fn2bx_21yxlslz9kqc0000gn/T/jl_Zjx62t/Manifest.toml`
  [b99e7846] BinaryProvider v0.5.10
  [46823bd8] Chemfiles v0.10.2
  [34da2185] Compat v3.41.0
  [ffbed154] DocStringExtensions v0.8.6
  [692b3bcd] JLLWrappers v1.4.0
  [21216c6a] Preferences v1.2.3
  [78a364fa] Chemfiles_jll v0.10.2+0
  [0dad84c5] ArgTools `@stdlib/ArgTools`
  [56f22d72] Artifacts `@stdlib/Artifacts`
  [2a0f44e3] Base64 `@stdlib/Base64`
  [ade2ca70] Dates `@stdlib/Dates`
  [8bb1440f] DelimitedFiles `@stdlib/DelimitedFiles`
  [8ba89e20] Distributed `@stdlib/Distributed`
  [f43a241f] Downloads `@stdlib/Downloads`
  [b77e0a4c] InteractiveUtils `@stdlib/InteractiveUtils`
  [b27032c2] LibCURL `@stdlib/LibCURL`
  [76f85450] LibGit2 `@stdlib/LibGit2`
  [8f399da3] Libdl `@stdlib/Libdl`
  [37e2e46d] LinearAlgebra `@stdlib/LinearAlgebra`
  [56ddb016] Logging `@stdlib/Logging`
  [d6f4376e] Markdown `@stdlib/Markdown`
  [a63ad114] Mmap `@stdlib/Mmap`
  [ca575930] NetworkOptions `@stdlib/NetworkOptions`
  [44cfe95a] Pkg `@stdlib/Pkg`
  [de0858da] Printf `@stdlib/Printf`
  [3fa0cd96] REPL `@stdlib/REPL`
  [9a3f8284] Random `@stdlib/Random`
  [ea8e919c] SHA `@stdlib/SHA`
  [9e88b42a] Serialization `@stdlib/Serialization`
  [1a1011a3] SharedArrays `@stdlib/SharedArrays`
  [6462fe0b] Sockets `@stdlib/Sockets`
  [2f01184e] SparseArrays `@stdlib/SparseArrays`
  [10745b16] Statistics `@stdlib/Statistics`
  [fa267f1f] TOML `@stdlib/TOML`
  [a4e569a6] Tar `@stdlib/Tar`
  [8dfed614] Test `@stdlib/Test`
  [cf7118a7] UUIDs `@stdlib/UUIDs`
  [4ec0a83e] Unicode `@stdlib/Unicode`
  [e66e0078] CompilerSupportLibraries_jll `@stdlib/CompilerSupportLibraries_jll`
  [deac9b47] LibCURL_jll `@stdlib/LibCURL_jll`
  [29816b5a] LibSSH2_jll `@stdlib/LibSSH2_jll`
  [c8ffd9c3] MbedTLS_jll `@stdlib/MbedTLS_jll`
  [14a3606d] MozillaCACerts_jll `@stdlib/MozillaCACerts_jll`
  [4536629a] OpenBLAS_jll `@stdlib/OpenBLAS_jll`
  [83775a58] Zlib_jll `@stdlib/Zlib_jll`
  [8e850b90] libblastrampoline_jll `@stdlib/libblastrampoline_jll`
  [8e850ede] nghttp2_jll `@stdlib/nghttp2_jll`
  [3f19e933] p7zip_jll `@stdlib/p7zip_jll`
Precompiling project...
  ✗ BinaryProvider
  0 dependencies successfully precompiled in 2 seconds (12 already precompiled)
  1 dependency errored. To see a full report either run `import Pkg; Pkg.precompile()` or load the package
     Testing Running tests...
ERROR: LoadError: InitError: UndefVarError: libchemfiles not defined
Stacktrace:
  [1] chfl_version
    @ ~/.julia/packages/Chemfiles/QiqBE/src/generated/cdef.jl:19 [inlined]
  [2] version
    @ ~/.julia/packages/Chemfiles/QiqBE/src/Chemfiles.jl:36 [inlined]
  [3] __init__()
    @ Chemfiles ~/.julia/packages/Chemfiles/QiqBE/src/Chemfiles.jl:145
  [4] _include_from_serialized(path::String, depmods::Vector{Any})
    @ Base ./loading.jl:768
  [5] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String)
    @ Base ./loading.jl:854
  [6] _require(pkg::Base.PkgId)
    @ Base ./loading.jl:1097
  [7] require(uuidkey::Base.PkgId)
    @ Base ./loading.jl:1013
  [8] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:997
  [9] include(fname::String)
    @ Base.MainInclude ./client.jl:451
 [10] top-level scope
    @ none:6
during initialization of module Chemfiles
in expression starting at /Users/siida/.julia/packages/Chemfiles/QiqBE/test/runtests.jl:1
ERROR: Package Chemfiles errored during testing

The upcoming version of Julia (0.7 / 1.0) produces several warnings relating to changes in the language. Unfortunately, fixing these warnings breaks Julia 0.6 support. As v0.6 is the current stable release, these warnings should not be addressed until the official release of Julia v0.7. Note that Julia v1.0 will not allow for these depreciations, so at some point support for v0.6 will need to be dropped - but now is not that time.

Using Molfiles plugin on Windows fails (cf https://ci.appveyor.com/project/Luthaf/chemfiles-jl/build/30#L105), and the log says

Chemfiles error: Could not find the 'gromacsplugin.so' shared library. Try setting the CHEMFILES_PLUGINS environement variable.

At the same time, CHEMFILES_PLUGINS is set correctly, and the gromacsplugin.so file exists ...

This is the position matrix:

julia> xyz
3×5 Matrix{Float64}
  0.756778  -0.756778  0.0         3.81242    3.23075
  0.0        0.0       0.0         0.395536   0.109453
 -0.522494  -0.538596  0.0668583  -1.76606   -2.69293

It has full precision:

julia> xyz[:,1]
3-element Vector{Float64}:
  0.7567777232273633
  0.0
 -0.5224940120543102

After the position is written to a frame, the frame still has full precision:

julia> frame
Frame(Chemfiles.CxxPointer{Chemfiles.lib.CHFL_FRAME}(Ptr{Chemfiles.lib.CHFL_FRAME} @0x00000000053ec150, false))
julia> positions(frame)[:,1]
3-element Vector{Float64}:
  0.7567777232273633
  0.0
 -0.5224940120543102

But the file the trajectory is written to loses precision:

5
Properties=species:S:1:pos:R:3
H 0.756778 0 -0.522494
H -0.756778 0 -0.538596
O 0 0 0.0668583
O 3.81242 0.395536 -1.76606
C 3.23075 0.109453 -2.69293

I'm playing with both the bindings for Julia and Python, and while working with residues on a liposome, I found that, in the Julia version, Residue and residue_for_atom are making copies of the residue instead of giving a view, like in the Python version where a raed-only access is given. The problem is that the Julia version is unusable for me because of the memory consumption that comes with doing so much copies of potentially big data structures.
Is there a workaround that I haven't been able to find? Otherwise, I think that we should be using views or giving the possibility to choose if the function does a copy or not

I saw from the docs of the C++ interface that Frame objects can be indexed to get atoms from the frame. But this is not implemented here. This will be very useful to have here also. If you think this is good, I'll send a PR for this.

Dear Guillaume, I could not reach you otherwise, so I am posting it here, although this is not an issue.

I wrote a very ugly parser of the selection of Chemfiles that converts the input
and output of a selection string from a one-based indexing of atoms to a zero-based, and
then back, such that I can select atom "1" using

"index 1"

and get back the number "1" in the list, using Chemfiles in background.

Probably what I have done is pretty ugly, but it seems to work fine. For example:

index_to_zero_based("resname ALA and index < 5 or index > 150 or index 5 7 17")

returns:

"resname ALA and index < 4 or index > 149 or index 4 6 16"

which then can be fed into the Chemfiles selector. Of course, afterwards another parser
increases by one every returning index such that the user gets the one-based indexing as a result again.

The function is just a hack, probably you would like something more professional, but maybe it is useful:

  #
  # Convert indexes to zero-based indexing to pass to Chemfiles
  #
  function indexes_to_zero_based(selection :: String)
    if occursin("index",selection)
      s = split(selection)
      ikey = 1
      while ikey <= length(s)
        if s[ikey] == "index"
          ikey += 1
          if tryparse(Int,s[ikey]) == nothing
            ikey += 1
          end
          index = tryparse(Int,s[ikey])
          while ikey <= length(s) && index != nothing
            index = index - 1
            s[ikey] = "$index"
            ikey += 1
            if ikey <= length(s)
              index = tryparse(Int,s[ikey])
            end
          end
        end
        ikey += 1
      end
      selection = ""
      for item in s
        selection = selection*"$item "
      end
      selection = rstrip(selection)
    end
    return selection
  end

This would improve the documentation greatly !

In functions like name(Atom), string(Selection) and others, we use fixed sized buffers. If the string does not fit in the buffer, then the string is truncated, which can be a surprising behavior.

Instead, we should automatically grow the buffer to make sure that we got the full string from the C API. Python is already doing this: https://github.com/chemfiles/chemfiles.py/blob/03d3cc880a88919d3f232a1b54de9a23e496977a/chemfiles/utils.py#L126-L147

These are the input and output PDBs I get from these lines of code:

in_trj = Trajectory("in.pdb")
in_frm = read(in_trj)

out_trj = Trajectory("out.pdb", 'w')
write(out_trj, in_frm)

I figured it must be something wrong with the C API, but then I realized that the downloaded chemfiles shared libraries belong to the 0.704 version, even though the Julia package is up to date.

What are the necessary steps to update bindings? Especially for functions that are already present.

Also, I noticed the Julia bindings documentation makes reference to a write!() function instead of the current write(). I guess the function should be renamed* since the Trajectory object is changed during writing. Right?

*to comply with Julia's policy on functions that change their input parameters.

edit: I also tried to do the same on the C API but I couldn't even compile the code.
Jeez... and I thought g++ errors were unhelpful; it's even worse for C.

I have attached here the DCD files required to reproduce the issue.

The problem is manifested in a a periodic box with Triclinic periodic boundary conditions.

The unit cell matrix used to generate the simulation was:

    unit_cell = [80.0 0.0 30.0
                 30.0 80.0 0.0
                 0.0 40.0 80.0]

which, for instance has a volume of det(unit_cell) = 548000.0.

When this file, t1.dcd is open in VMD, one can get the unit cell using pbc get, and the result

vmd > pbc get 
{85.440041 89.442719 85.440041 65.244987 70.806038 71.696266}

which returns lengths and angles for this unit cell.

If, however, I open the file with Chemfiles, I get:

julia> traj = Chemfiles.Trajectory("./t1.dcd");

julia> frame = Chemfiles.read_step(traj, 0);

julia> uc = Chemfiles.UnitCell(frame);

julia> matrix(uc)
3×3 Matrix{Float64}:
 85.44      0.0       0.0
  7.13576  89.1576    0.0
  7.14719   8.66822  84.6982

julia> lengths(uc)
3-element Vector{Float64}:
 85.44003745317531
 89.44271909999159
 85.44003745317531

julia> angles(uc)
3-element Vector{Float64}:
 83.81124760738278
 85.20150950650077
 85.42406614347279

julia> volume(uc)
645199.3227509324

Where clearly the matrix is not the same, the angles do not match the VMD reported angles and, clearly, the volume of the cell does not match the volume of the original cell.

If I fix the angles to the ones reported by VMD:

julia> set_angles!(uc, [65.244987, 70.806038, 71.696266]);

julia> volume(uc)
547999.9898359198

then the volume is correct.

Thus, I think there is a problem in the reading of the angles from the DCD trajectory, which clearly manifest itself when using triclinic periodic boundary conditions.

For the records, the angles are properly read if the cell is orthorhombic (file o1.dcd):

julia> traj = Chemfiles.Trajectory("./o1.dcd");

julia> frame = Chemfiles.read_step(traj, 0);

julia> uc = Chemfiles.UnitCell(frame);

julia> angles(uc)
3-element Vector{Float64}:
 90.0
 90.0
 90.0

namd.zip

With julia-1.5.3:

julia> Pkg.test("Chemfiles")

led to:

`Test Summary: | Pass Total
Generics | 1 1
Test Summary: | Pass Total
Atom | 16 16
WARNING: both Chemfiles and Base export "contains"; uses of it in module Main must be qualified
Residue: Error During Test at $HOME/.julia/packages/Chemfiles/uJPdJ/test/Residue.jl:22
Test threw exception
Expression: contains(residue, 56) == true
UndefVarError: contains not defined
Stacktrace:
[1] top-level scope at $HOME.julia/packages/Chemfiles/uJPdJ/test/Residue.jl:22
[2] top-level scope at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/Test/src/Test.jl:1115
[3] top-level scope at $HOME/.julia/packages/Chemfiles/uJPdJ/test/Residue.jl:2

Test Summary: | Pass Error Total
Residue | 12 1 13
ERROR: LoadError: LoadError: Some tests did not pass: 12 passed, 0 failed, 1 errored, 0 broken.
in expression starting at $HOME/.julia/packages/Chemfiles/uJPdJ/test/Residue.jl:1
in expression starting at $HOME/.julia/packages/Chemfiles/uJPdJ/test/runtests.jl:28
ERROR: Package Chemfiles errored during testing`

The problem was circumvented in Residue.jl (line: ) by using:
@test **Chemfiles.contains**(residue, 56) == true

Chemfiles.jl is using sphinx-julia for its documentation. But there seems to be some issues with how references to functions are handled (:jl:func:`foobar!` )

They mostly work, but sometimes an exception is raised:

Traceback (most recent call last):
  File "/usr/local/lib/python2.7/site-packages/sphinx/cmdline.py", line 304, in main
    app.build(args.force_all, filenames)
  File "/usr/local/lib/python2.7/site-packages/sphinx/application.py", line 331, in build
    self.builder.build_update()
  File "/usr/local/lib/python2.7/site-packages/sphinx/builders/__init__.py", line 333, in build_update
    'out of date' % len(to_build))
  File "/usr/local/lib/python2.7/site-packages/sphinx/builders/__init__.py", line 394, in build
    self.write(docnames, list(updated_docnames), method)
  File "/usr/local/lib/python2.7/site-packages/sphinx/builders/__init__.py", line 431, in write
    self._write_serial(sorted(docnames))
  File "/usr/local/lib/python2.7/site-packages/sphinx/builders/__init__.py", line 438, in _write_serial
    doctree = self.env.get_and_resolve_doctree(docname, self)
  File "/usr/local/lib/python2.7/site-packages/sphinx/environment/__init__.py", line 813, in get_and_resolve_doctree
    self.apply_post_transforms(doctree, docname)
  File "/usr/local/lib/python2.7/site-packages/sphinx/environment/__init__.py", line 860, in apply_post_transforms
    transformer.apply_transforms()
  File "/usr/local/lib/python2.7/site-packages/sphinx/transforms/__init__.py", line 96, in apply_transforms
    Transformer.apply_transforms(self)
  File "/usr/local/lib/python2.7/site-packages/docutils/transforms/__init__.py", line 171, in apply_transforms
    transform.apply(**kwargs)
  File "/usr/local/lib/python2.7/site-packages/sphinx/transforms/post_transforms/__init__.py", line 90, in apply
    typ, target, node, contnode)
  File "/usr/local/lib/python2.7/site-packages/sphinxjulia/juliadomain.py", line 162, in resolve_xref
    node)
  File "/usr/local/lib/python2.7/site-packages/sphinx/environment/__init__.py", line 305, in warn_node
    self._warnfunc(msg, '%s:%s' % get_source_line(node), **kwargs)
TypeError: 'NoneType' object is not callable

It looks like these exception are raised when a function is documented more than once (for example set_cell! have overloads for frame and trajectory)

At least in DCD trajectories obtained with NAMD, the dimensions of the cell are available in the trajectory file itself. I was unable to find how to retrieve this information when reading a frame. It seems that there should be a function similar to position(frame) to get these dimensions. From the documentation I could only find ways to define the cell shape, not to retrieve it from the trajectory. If there is such a function, I am sorry for opening the issue and I will appreciate an indication of the associated documentation.

Like in https://github.com/chemfiles/chemfiles.f03/blob/29ecc3f0830c82172ab27bb9b7ae94a5b5eb1083/scripts/check-docs.py.

This should only be a matter of taking the above script and adapting it to julia.

Feel free to jump in and ask for help if you want to do this!

While trying to read LAMMPS dump file, it seems that the properties are not read. I think this is related to this issue: chemfiles/chemfiles#409 in the chemfiles repo. But that seems to have been fixed there with chemfiles/chemfiles#413. A MWE using the test input given in the previous PR:

julia> using Chemfiles

julia> Trajectory("lj_dump.lammpstrj", 'r', "LAMMPS") do trj
           list_properties(read(trj))
       end
String[]

Anyway thanks for the great package!

The function is documented, but I didn't find it in the code.
Thanks!

Hello,

I am using Chemfiles.jl to work with GROMACS files and, as shown in https://chemfiles.org/chemfiles/latest/formats.html#list-of-supported-formats, I thought that .tpr files could be read using Trajectory. But, when I try to open a .tpr file, it shows that there is no current support for this extension. Is any plan to include support to this extension in the future?

Hi,

I'm using the Julia bindings to Chemfiles to read a LAMMPS dump file (file extension .lammpstrj). The topology information is stored in a LAMMPS data file (file extension .data). For some reason, when I index the positions read out from the LAMMPS dump using indices provided by evaluating a selection on the frame, I reliably but non-deterministically get segmentation faults that crash the Julia REPL. Here's a minimal working example with Julia Chemfiles 0.9.3:

using Chemfiles

top_dir = "electrode-sim-files"
traj_fname = joinpath(top_dir, "dumptest.lammpstrj")
topo_fname = joinpath(top_dir, "dumptest.data")

traj = Trajectory(traj_fname, 'r', "LAMMPS")
topo_traj = Trajectory(topo_fname, 'r', "LAMMPS Data")
topo_frame = read(topo_traj)

selstr = "not type == \"16\" or type == \"17\""
sel = Selection(selstr)

# this loop runs successfully
for ix in 1:length(traj)
    frame = read_step(traj, ix - 1)
    pos = positions(frame)
    println(sum(pos))
end

# this loop segfaults
for ix in 1:length(traj)
    frame = read_step(traj, ix - 1)
    pos = positions(frame)
    indices = evaluate(sel, frame)
    println(sum(pos[:, indices .+ 1]))
end

close(traj)

Here's the stack trace I get from the segfault, it appears to be somewhere in the getindex function for the positions array:

signal (11): Segmentation fault: 11
in expression starting at /Users/aditya/research/023_tba_electrolyte/clean.jl:22
macro expansion at ./multidimensional.jl:697 [inlined]
macro expansion at ./cartesian.jl:64 [inlined]
macro expansion at ./multidimensional.jl:694 [inlined]
_unsafe_getindex! at ./multidimensional.jl:690 [inlined]
_unsafe_getindex at ./multidimensional.jl:684
_getindex at ./multidimensional.jl:670 [inlined]
getindex at ./abstractarray.jl:981
top-level scope at /Users/aditya/research/023_tba_electrolyte/clean.jl:26

Could you please provide some advice on what might be going wrong here? The code above does not always segfault, but in the few times I've tried it, usually segfaults within the first 2-3 times running it. I'm attaching the offending LAMMPS dump and data files here as well. Thanks in advance for any help you can provide!

files.zip

Looks like appveyor is no longer building and checking the master branch.

It should either be fixed or replaced with Travis CI for windows.

@JuliaRegistrator register

I had started working on the new version, for Julia 0.6 and Chemfiles 0.7.4, but only now I have seen that your WIP bdf8eef

What's the next steps?
I can put efforts on this and help you. 👍