Merge pull request #25 from m3g/fixing_coordnum

Fixing coordnum

src/mddf.jl

@@ -176,13 +176,13 @@ function mddf(trajectory::Trajectory, options::Options = Options(); coordination
             # Read frame coordinates
             lock(read_lock) do
                 # skip frames if stride > 1
-                while (iframe + 1) % options.stride != 0
+                while (iframe - options.firstframe + 1) % options.stride != 0
                     nextframe!(trajectory)
                     iframe += 1
                 end
                 # Read frame for computing 
-                iframe += 1
                 nextframe!(trajectory)
+                iframe += 1
                 # The solute coordinates must be read in intermediate arrays, because the 
                 # solute molecules will be considered one at a time in the computation of the
                 # minimum distances
@@ -494,6 +494,11 @@ end
     R2 = mddf(traj2, Options(lastframe=2, frame_weights=[0.5, 0.25]))
     @test R2.md_count ≈ R1.md_count
     @test R2.volume.total ≈ R1.volume.total
+    # Varying weights with stride
+    R1 = mddf(traj1, Options(firstframe=2, lastframe=3))
+    R2 = mddf(traj1, Options(lastframe=3, stride=2, frame_weights=[1.0, 100.0, 1.0]))
+    @test R2.md_count ≈ R1.md_count
+    @test R2.volume.total ≈ R1.volume.total
 
 end
 

src/results.jl

@@ -402,6 +402,7 @@ function finalresults!(R::Result, options::Options, trajectory::Trajectory)
     #
     # Computing the distribution functions and KB integrals, from the MDDF and from the RDF
     #
+    warn = false
     for ibin = 1:R.nbins
         # For the MDDF
         if R.md_count_random[ibin] > 0.0
@@ -412,6 +413,11 @@ function finalresults!(R::Result, options::Options, trajectory::Trajectory)
             for j = 1:trajectory.solvent.natomspermol
                 R.solvent_atom[ibin, j] = R.solvent_atom[ibin, j] / R.md_count_random[ibin]
             end
+        else
+            if !warn && !options.silent
+                @warn "Ideal-gas histogram bins with zero samples. Increase n_random_samples or trajectory length."
+                warn = true
+            end
         end
         if ibin == 1
             R.coordination_number[ibin] = R.md_count[ibin]

src/tools/contributions.jl

@@ -83,7 +83,7 @@ function contributions(
     warning = true,
 )
     (warning && s.nmols > 1) && warning_nmols_types()
-    indexes = [atom.index for atom in atoms]
+    indexes = PDBTools.index.(atoms)
     # Check which types of atoms belong to this selection
     selected_types = which_types(s, indexes, warning = warning)
     return contributions(s, atom_contributions, selected_types)

src/tools/coordination_number.jl

@@ -73,6 +73,10 @@ function coordination_number(
     return coordination_number(R, group_contributions)
 end
 function coordination_number(R::Result, group_contributions::Vector{Float64})
+    # obs: we accumulate group coordination numbers, but at the end we 
+    # divide by random count of, so we the group contributions are the contributions 
+    # to the MDDF. To recover the coordination number we need to multiply back
+    # the group contributions by the random count.
     cn = cumsum(
         group_contributions[i] * R.md_count_random[i] for
         i in eachindex(group_contributions)
@@ -82,16 +86,30 @@ end
 
 @testitem "coordination_number" begin
     using ComplexMixtures, PDBTools
-    import ComplexMixtures.Testing: test_dir
-
-    pdb = readPDB("$test_dir/data/NAMD/structure.pdb")
-    R = load("$test_dir/data/NAMD/protein_tmao.json"; legacy_warning = false)
-    solute = Selection(PDBTools.select(pdb, "protein"), nmols = 1)
-    cn = coordination_number(solute, R.solute_atom, R, PDBTools.select(pdb, "residue 50"))
-    @test cn[findlast(<(5), R.d)] ≈ 0.24999999999999997 atol = 1e-10
-
-    pdb = readPDB("$test_dir/data/NAMD/Protein_in_Glycerol/system.pdb")
-    R = load("$test_dir/data/NAMD/Protein_in_Glycerol/protein_water.json"; legacy_warning=false)
+    using ComplexMixtures.Testing
+    dir = "$(Testing.data_dir)/NAMD"
+    atoms = readPDB("$dir/structure.pdb")
+    protein = Selection(select(atoms, "protein"), nmols = 1)
+    tmao = Selection(select(atoms, "resname TMAO"), natomspermol = 14)
+    options = Options(lastframe = 1, seed = 321, StableRNG = true, nthreads = 1, silent = true, n_random_samples=200)
+    traj = Trajectory("$dir/trajectory.dcd", protein, tmao)
+    R = mddf(traj, options)
+
+    # Test self-consistency
+    @test sum(R.solute_atom) ≈ sum(R.solvent_atom)
+    @test coordination_number(R, contributions(protein, R.solute_atom, select(atoms, "protein"))) == R.coordination_number
+
+    # Checked with vmd: same residue as (resname TMAO and within 3.0 of protein)
+    @test R.coordination_number[findfirst(>(3), R.d)] == 7.0
+    @test R.coordination_number[findfirst(>(5), R.d)] == 14.0
+
+    # THR93 forms a hydrogen bond with TMAO in this frame
+    thr93 = select(atoms, "protein and residue 93")
+    @test last(coordination_number(R, contributions(protein, R.solute_atom, thr93))) == 1
+
+    # Consistency of a read-out result
+    pdb = readPDB("$(Testing.data_dir)/NAMD/Protein_in_Glycerol/system.pdb")
+    R = load("$(Testing.data_dir)/NAMD/Protein_in_Glycerol/protein_water.json"; legacy_warning=false)
     group = PDBTools.select(pdb, "protein")
     solute = Selection(group, nmols = 1)
     cn = coordination_number(solute, R.solute_atom, R, group)

src/updatecounters.jl

@@ -43,7 +43,7 @@ function updatecounters!(R::Result, system::AbstractPeriodicSystem, frame_weight
         R.md_count_random[ibin] += frame_weight
         if md.ref_atom_within_cutoff
             ibin = setbin(md.d_ref_atom, R.options.binstep)
-            R.rdf_count_random[ibin] += 1
+            R.rdf_count_random[ibin] += frame_weight
         end
     end
     return R