Merge pull request #72 from m3g/testing_with_julia

Testing with julia

.github/workflows/actions.yml

@@ -40,8 +40,5 @@ jobs:
     - name: Run
       run: |
           cd testing
-          ./test.sh "water_box.inp" "packmol.log" "Success!"
-          ./test.sh "ieee_signaling.inp" "packmol.log" "Success!"
-          ./test.sh "water_box_failed.inp" "packmol.log" "FORCED"
-          cd ..
+          ./test.sh
       shell: bash

.gitignore

@@ -11,3 +11,6 @@ testing/ieee_signaling_box.pdb
 testing/water_box.pdb
 testing/water_box.pdb_FORCED
 testing/packmol.log
+testing/Manifest.toml
+testing/output.pdb
+testing/output.pdb_FORCED

testing/Project.toml

@@ -0,0 +1,3 @@
+[deps]
+CellListMap = "69e1c6dd-3888-40e6-b3c8-31ac5f578864"
+PDBTools = "e29189f1-7114-4dbd-93d0-c5673a921a58"

testing/input_files/bilayer.inp

@@ -0,0 +1,40 @@
+#
+# Lipid double layer with water above and below
+#
+
+tolerance 2.0 
+filetype pdb
+output output.pdb
+
+structure ./structure_files/water.pdb 
+  number 50
+  inside box 0. 0. -10. 40. 40. 0.
+end structure
+
+structure ./structure_files/water.pdb
+  number 50
+  inside box 0. 0. 28. 40. 40. 38.
+end structure
+
+structure ./structure_files/palmitoil.pdb 
+  number 10
+  inside box 0. 0. 0. 40. 40. 14.
+  atoms 31 32
+    below plane 0. 0. 1. 2.
+  end atoms
+  atoms 1 2
+    above plane 0. 0. 1. 12.
+  end atoms
+end structure 
+
+structure ./structure_files/palmitoil.pdb 
+  number 10
+  inside box 0. 0. 14. 40. 40. 28.
+  atoms 1 2
+    below plane 0. 0. 1. 16.
+  end atoms
+  atoms 31 32
+    above plane 0. 0. 1. 26 
+  end atoms
+end structure 
+

testing/ieee_signaling.inp → testing/input_files/ieee_signaling.inp

@@ -1,22 +1,22 @@
 tolerance 6.0
-output ieee_signaling_box.pdb
+output output.pdb
 filetype pdb
-structure ieee_signaling.pdb
+structure ./structure_files/ieee_signaling.pdb
    number 6916
    inside sphere 0. 0. 0. 385.
    outside sphere 0. 0. 0. 370.
 end structure
-structure ieee_signaling.pdb
+structure ./structure_files/ieee_signaling.pdb
    number 7377
    inside sphere 0. 0. 0. 315.
    outside sphere 0. 0. 0. 300.
 end structure
-structure ieee_signaling.pdb
+structure ./structure_files/ieee_signaling.pdb
    number 4783
    inside sphere 0. 0. 0. 255.
    outside sphere 0. 0. 0. 240.
 end structure
-structure ieee_signaling.pdb
+structure ./structure_files/ieee_signaling.pdb
    number 2766
    inside sphere 0. 0. 0. 205.
    outside sphere 0. 0. 0. 195.

testing/input_files/mixture.inp

@@ -0,0 +1,18 @@
+#
+# A mixture of water and urea
+#
+
+tolerance 2.0
+filetype pdb
+output output.pdb
+
+structure ./structure_files/water.pdb 
+  number 1000 
+  inside box 0. 0. 0. 40. 40. 40. 
+end structure
+
+structure ./structure_files/urea.pdb
+  number 400
+  inside box 0. 0. 0. 40. 40. 40. 
+end structure
+

testing/input_files/solvprotein.inp

@@ -0,0 +1,32 @@
+#
+# Protein solvated by water and ions
+#
+
+tolerance 2.0
+
+structure ./structure_files/protein.pdb
+  resnumbers 0
+  number 1 
+  fixed 0. 0. 0. 0. 0. 0.
+  centerofmass
+end structure
+
+structure ./structure_files/water.pdb
+  resnumbers 2
+  number 1000 
+  inside sphere 0. 0. 0. 50.
+end structure
+
+structure ./structure_files/CLA.pdb
+  resnumbers 2
+  number 20
+  inside sphere 0. 0. 0. 50.
+end structure
+
+structure ./structure_files/SOD.pdb
+  resnumbers 2
+  number 30
+  inside sphere 0. 0. 0. 50.
+end structure
+
+output output.pdb

testing/input_files/spherical.inp

@@ -0,0 +1,42 @@
+#
+# Double layered vesicle with water inside and outside
+#
+
+tolerance 2.0 
+
+structure ./structure_files/palmitoil.pdb 
+  resnumbers 2
+  number 20
+  atoms 37
+    inside sphere 0. 0. 0. 14. 
+  end atoms
+  atoms 5
+    outside sphere 0. 0. 0. 26.
+  end atoms
+end structure 
+
+structure ./structure_files/palmitoil.pdb 
+  resnumbers 2
+  number 30
+  atoms 5
+    inside sphere 0. 0. 0. 29. 
+  end atoms
+  atoms 37
+    outside sphere 0. 0. 0. 41.
+  end atoms
+end structure 
+
+structure ./structure_files/water.pdb
+   resnumbers 2
+   number 50
+   inside sphere 0. 0. 0. 13.
+end structure
+
+structure ./structure_files/water.pdb
+   resnumbers 2
+   number 100
+   inside box -47.5 -47.5 -47.5 47.5 47.5 47.5
+   outside sphere 0. 0. 0. 43.
+end structure
+
+output output.pdb

testing/water_box.inp → testing/input_files/water_box.inp

@@ -4,9 +4,9 @@
 
 tolerance 2.0
 filetype pdb
-output water_box.pdb
+output output.pdb
 
-structure water.pdb
+structure ./structure_files/water.pdb
   number 1000
   inside box 0. 0. 0. 40. 40. 40.
 end structure

testing/water_box_failed.inp → testing/input_files/water_box_failed.inp

@@ -4,12 +4,12 @@
 
 tolerance 2.0
 filetype pdb
-output water_box.pdb
+output output.pdb
 
 maxit 2
 nloop 1
 
-structure water.pdb
+structure ./structure_files/water.pdb
   number 2
   inside box 0. 0. 0. 3.0 3.0 1.0
 end structure

testing/runtests.jl

@@ -0,0 +1,71 @@
+import Pkg; Pkg.activate(@__DIR__); Pkg.instantiate()
+using PDBTools
+using CellListMap
+
+# 
+# Note: this tests depend on each molecule being identified as a single residue
+# in the structure. Thus, the PDB file should have a single residue for each molecule,
+# for instance, a protein must have a single residue name for all its atoms.
+# The residue numbers are overwritten by the `resnumbers` option of the input
+# packmol files .
+#
+
+struct MinimumDistance d::Float64 end
+function update_mind(i, j, d2, pdb, md::MinimumDistance)
+    residue(pdb[i]) == residue(pdb[j]) && return md
+    MinimumDistance(min(sqrt(d2), md.d))
+end
+CellListMap.reducer(md1::MinimumDistance, md2::MinimumDistance) = MinimumDistance(min(md1.d,md2.d))
+CellListMap.copy_output(md::MinimumDistance) = md
+CellListMap.reset_output(::MinimumDistance) = MinimumDistance(+Inf)
+
+function check_mind(input_file::String)
+    tolerance = nothing
+    output_name = nothing
+    unitcell = nothing
+    precision = 0.01
+    for line in eachline(input_file)
+        line = strip(line)
+        isempty(line) && continue
+        keyword, values... = split(line)
+        keyword == "tolerance" && (tolerance = parse(Float64, values[1]))
+        keyword == "output" && (output_name = values[1])
+        keyword == "pbc" && (unitcell = (values[1], values[2], values[3]))
+        keyword == "precision" && (precision = values[1])
+    end
+    if isnothing(tolerance) || isnothing(output_name)
+        error("tolerance or output not found")
+    end
+    pdb = readPDB(string(output_name))
+    sys = ParticleSystem(
+        positions = coor.(pdb),
+        unitcell = unitcell,
+        cutoff = tolerance * 1.5,
+        output = MinimumDistance(+Inf),
+    )
+    mind = map_pairwise((x,y,i,j,d2,md) -> update_mind(i, j, d2, pdb, md), sys)
+    if (mind.d < (1 - precision) * tolerance)
+        error("""\n
+
+            Packing reported success, but minimum distance is not correct for $input_file
+            Obtained minimum-distance = $(mind.d) for tolerance $tolerance and precision $precision.
+            
+        """)
+    end
+    return nothing
+end 
+
+if !isinteractive()
+    packmol = joinpath(@__DIR__,"..","packmol")
+    for input_test in ARGS
+        log = IOBuffer()
+        run(pipeline(`$packmol`; stdin=input_test, stdout=log))
+        if occursin("Success!", String(take!(log)))
+            check_mind(input_test)
+        else
+            error("Failed packing for $input_test")
+        end
+    end
+end
+
+

testing/structure_files/CLA.pdb

@@ -0,0 +1,2 @@
+HEADER
+HETATM    1 CLA  CLA    2        0.000   0.000   0.000

testing/structure_files/SOD.pdb

@@ -0,0 +1,2 @@
+HEADER
+HETATM    1 SOD  SOD    2        0.000   0.000   0.000