Skip to content

Commit

Permalink
Merge pull request #64 from m3g/calls_to_coordination_number
Browse files Browse the repository at this point in the history
Calls to coordination number
  • Loading branch information
lmiq authored Dec 8, 2024
2 parents e431dea + 0384786 commit 71e9e52
Show file tree
Hide file tree
Showing 5 changed files with 69 additions and 15 deletions.
32 changes: 24 additions & 8 deletions src/Trajectory.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,7 @@
Trajectory constructor data type.
Defaults to reading with the Chemfiles infrastructure, except for DCD and PDB trajectory
files, if the "PDBTraj" option is provided.
Defaults to reading with the Chemfiles infrastructure, except for DCD and PDB trajectory files, if the "PDBTraj" option is provided.
See memory issue (https://github.com/chemfiles/Chemfiles.jl/issues/44)
Expand All @@ -21,6 +20,9 @@ end
stream(traj::Trajectory) = traj.stream.st
set_stream!(traj::Trajectory, st) = traj.stream.st = st

# Trajectory formats
const trajectory_formats = String[]

# Include specific predefined trajectory types
include("./trajectory_formats/ChemFiles.jl")
include("./trajectory_formats/NamdDCD.jl")
Expand All @@ -33,20 +35,28 @@ function Trajectory(
format::String="",
chemfiles=false,
)
if !isempty(format) && !(format in trajectory_formats)
throw(ArgumentError("""\n
Trajectory format not properly set: $format
Available trajectory formats: $(join(trajectory_formats, ", "))
"""))
end
filename = expanduser(filename) # expand tilde on Unix systems, to username
if !chemfiles && (format == "dcd" || split(filename, '.')[end] == "dcd")
trajectory = NamdDCD(filename, solute, solvent)
elseif !chemfiles && format == "PDBTraj"
trajectory = PDBTraj(filename, solute, solvent)
file_extension = filename[findlast(==('.'), filename)+1:end]
trajectory = if !chemfiles && (format == "dcd" || file_extension == "dcd")
NamdDCD(filename, solute, solvent)
elseif !chemfiles && (format == "PDBTraj" || file_extension =="pdb")
PDBTraj(filename, solute, solvent)
else
trajectory = ChemFile(filename, solute, solvent, format=format)
ChemFile(filename, solute, solvent, format=format)
end
return trajectory
end

# If only one selection is provided, assume that the solute and the solvent are the same
Trajectory(filename::String, solvent::AtomSelection; format::String="", chemfiles=false) =
Trajectory(filename, solvent, solvent, format=format, chemfiles=chemfiles)
Trajectory(filename, solvent, solvent; format, chemfiles)

#=
convert_unitcell(unitcell::Union{SVector{3}, SMatrix{3,3}})
Expand Down Expand Up @@ -108,6 +118,9 @@ end
protein = AtomSelection(select(atoms, "protein"), nmols=1)
tmao = AtomSelection(select(atoms, "resname TMAO"), natomspermol=14)

# Throw if trajectory type is not available
@test_throws ArgumentError Trajectory("$(Testing.data_dir)/PDB/trajectory.pdb", protein, tmao; format="XXX")

# NAMD DCD file
traj = Trajectory("$(Testing.data_dir)/NAMD/trajectory.dcd", protein, tmao)
@test traj.nframes == 20
Expand All @@ -124,6 +137,9 @@ end
@test ComplexMixtures.natoms(traj.solute) == 1463
@test ComplexMixtures.natoms(traj.solvent) == 2534
@test ComplexMixtures.convert_unitcell(ComplexMixtures.getunitcell(traj)) SVector(84.47962951660156, 84.47962951660156, 84.47962951660156)
# Test determining format from file extension
traj = Trajectory("$(Testing.data_dir)/PDB/trajectory.pdb", protein, tmao)
@test ComplexMixtures.natoms(traj.solute) == 1463

# Chemfiles with NAMD
traj = Trajectory(
Expand Down
43 changes: 36 additions & 7 deletions src/mddf.jl
Original file line number Diff line number Diff line change
Expand Up @@ -183,7 +183,7 @@ function mddf(
trajectory_file::String,
solute::AtomSelection,
solvent::AtomSelection,
options::Options;
options::Options=Options();
trajectory_format::String="",
chemfiles::Bool=false,
kargs...)
Expand All @@ -194,7 +194,7 @@ end
function mddf(
trajectory_file::String,
solute_and_solvent::AtomSelection,
options::Options;
options::Options=Options();
trajectory_format::String="",
chemfiles::Bool=false,
kargs...
Expand Down Expand Up @@ -528,7 +528,7 @@ function coordination_number(
trajectory_file::String,
solute::AtomSelection,
solvent::AtomSelection,
options::Options;
options::Options=Options();
trajectory_format::String="",
chemfiles::Bool=false,
kargs...)
Expand All @@ -540,7 +540,7 @@ end
function coordination_number(
trajectory_file::String,
solute_and_solvent::AtomSelection,
options::Options;
options::Options=Options();
trajectory_format::String="",
chemfiles::Bool=false,
kargs...
Expand Down Expand Up @@ -608,9 +608,6 @@ end
# Test wrong frame_weights input
@test_throws ArgumentError mddf(trajectory_file, protein, water, Options(); frame_weights=[1.0], trajectory_format)

# The coordination_number(::String, args...; kargs...) is a placeholder for the docs only
@test_throws ArgumentError coordination_number("string.txt", 1.0)

# Self correlation
atoms = readPDB("$data_dir/toy/self_monoatomic.pdb")
atom = AtomSelection(select(atoms, "resname WAT and model 1"), natomspermol=1)
Expand Down Expand Up @@ -709,6 +706,38 @@ end
@test R.density.solvent 2 / R.volume.total
@test R.density.solvent_bulk 0.5 / R.volume.bulk
end
end

@testitem "mddf - available methods" begin
using ComplexMixtures
using PDBTools: readPDB, select
using ComplexMixtures.Testing: data_dir
atoms = readPDB("$data_dir/toy/cross.pdb")
protein = AtomSelection(select(atoms, "protein and model 1"), nmols=1)
water = AtomSelection(select(atoms, "resname WAT and model 1"), natomspermol=3)
trajectory_file = "$data_dir/toy/cross.pdb"

atom = AtomSelection(select(atoms, "resname WAT and model 1"), natomspermol=1)
t = Trajectory(trajectory_file, atom; format="PDBTraj")

options = Options()
trajectory_file = "$data_dir/toy/self_monoatomic.pdb"
atom = AtomSelection(select(atoms, "resname WAT and model 1"), natomspermol=1)
r1 = mddf(trajectory_file, atom, atom)
r2 = mddf(trajectory_file, atom, atom, Options())
@test coordination_number(r1) coordination_number(r2)
r1 = mddf(trajectory_file, atom)
r2 = mddf(trajectory_file, atom, Options())
@test coordination_number(r1) coordination_number(r2)
r1 = coordination_number(trajectory_file, atom, atom)
r2 = coordination_number(trajectory_file, atom, atom, Options())
@test coordination_number(r1) coordination_number(r2)
r1 = coordination_number(trajectory_file, atom)
r2 = coordination_number(trajectory_file, atom, Options())
@test coordination_number(r1) coordination_number(r2)

# The coordination_number(::String, args...; kargs...) is a placeholder for the docs only
@test_throws ArgumentError coordination_number("string.txt", 1.0)

end

Expand Down
3 changes: 3 additions & 0 deletions src/trajectory_formats/ChemFiles.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,9 @@
#
import Chemfiles

# Add this format to list of available formats
push!(trajectory_formats, "Chemfiles")

#"""
#
#$(TYPEDEF)
Expand Down
3 changes: 3 additions & 0 deletions src/trajectory_formats/NamdDCD.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,6 @@
# Add this format to list of available formats
push!(trajectory_formats, "dcd")

#
# Structure to contain DCD trajectories produces with Namd.
#
Expand Down
3 changes: 3 additions & 0 deletions src/trajectory_formats/PDBTraj.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,6 @@
# Add this format to list of available files
push!(trajectory_formats, "PDBTraj")

#"""
#
#$(TYPEDEF)
Expand Down

0 comments on commit 71e9e52

Please sign in to comment.