Skip to content

Commit a6ee262

Browse files
committedDec 5, 2023
in progress: updating coordination number testing
1 parent 2c17c13 commit a6ee262

File tree

1 file changed

+18
-5
lines changed

1 file changed

+18
-5
lines changed
 

‎src/tools/coordination_number.jl

+18-5
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,10 @@ function coordination_number(
7373
return coordination_number(R, group_contributions)
7474
end
7575
function coordination_number(R::Result, group_contributions::Vector{Float64})
76+
# obs: we accumulate group coordination numbers, but at the end we
77+
# divide by random count of, so we the group contributions are the contributions
78+
# to the MDDF. To recover the coordination number we need to multiply back
79+
# the group contributions by the random count.
7680
cn = cumsum(
7781
group_contributions[i] * R.md_count_random[i] for
7882
i in eachindex(group_contributions)
@@ -82,12 +86,21 @@ end
8286

8387
@testitem "coordination_number" begin
8488
using ComplexMixtures, PDBTools
85-
import ComplexMixtures.Testing: test_dir
89+
using ComplexMixtures.Testing
90+
dir = "$(Testing.data_dir)/NAMD"
91+
atoms = readPDB("$dir/structure.pdb")
92+
protein = Selection(select(atoms, "protein"), nmols = 1)
93+
tmao = Selection(select(atoms, "resname TMAO"), natomspermol = 14)
94+
options = Options(lastframe = 1, seed = 321, StableRNG = true, nthreads = 1, silent = true, n_random_samples=200)
95+
traj = Trajectory("$dir/trajectory.dcd", protein, tmao)
96+
R = mddf(traj, options)
97+
98+
# Test self-consistency
99+
@test sum(R.solute_atom) sum(R.solvent_atom)
100+
@test coordination_number(R, contributions(protein, R.solute_atom, select(atoms, "protein"))) == R.coordination_number
101+
102+
cn = coordination_number(protein, R.solute_atom, R, PDBTools.select(atoms, "residue 50"))
86103

87-
pdb = readPDB("$test_dir/data/NAMD/structure.pdb")
88-
R = load("$test_dir/data/NAMD/protein_tmao.json"; legacy_warning = false)
89-
solute = Selection(PDBTools.select(pdb, "protein"), nmols = 1)
90-
cn = coordination_number(solute, R.solute_atom, R, PDBTools.select(pdb, "residue 50"))
91104
@test cn[findlast(<(5), R.d)] 0.24999999999999997 atol = 1e-10
92105

93106
pdb = readPDB("$test_dir/data/NAMD/Protein_in_Glycerol/system.pdb")

0 commit comments

Comments
 (0)
Please sign in to comment.