minor tweaks

data/atefe/mor.py

@@ -2,20 +2,25 @@
 from openmm import *
 from openmm.unit import *
 from sys import stdout
-psf = CharmmPsfFile('/data/numerik/ag_cmd/arostami/charmm/charmm-gui(6)/charmm-gui-2501706787/openmm/step5_input.psf')
+
+a,b,c = (10,10,16.5)
+psf = CharmmPsfFile('/data/numerik/ag_cmd/arostami/charmm/charmm-gui(6)/charmm-gui-2501706787/openmm/step5_input.psf',
+                    unitCellDimensions=(a,b,c))
 pdb = PDBFile('/data/numerik/ag_cmd/arostami/charmm/charmm-gui(6)/charmm-gui-2501706787/openmm/step5_input.pdb')
 
 params = CharmmParameterSet('/data/numerik/ag_cmd/arostami/charmm/charmm-gui(6)/charmm-gui-2501706787/toppar/par_all36m_prot.prm', '/data/numerik/ag_cmd/arostami/charmm/charmm-gui(6)/charmm-gui-2501706787/toppar/top_all36_prot.rtf',
                        '/data/numerik/ag_cmd/arostami/charmm/charmm-gui(6)/charmm-gui-2501706787/toppar/par_all36_lipid.prm', '/data/numerik/ag_cmd/arostami/charmm/charmm-gui(6)/charmm-gui-2501706787/toppar/top_all36_lipid.rtf',
                        '/data/numerik/ag_cmd/arostami/charmm/charmm-gui(6)/charmm-gui-2501706787/toppar/toppar_water_ions.str', '/data/numerik/ag_cmd/arostami/charmm/charmm-gui(6)/charmm-gui-2501706787/7v7/7v7.prm', '/data/numerik/ag_cmd/arostami/charmm/charmm-gui(6)/charmm-gui-2501706787/7v7/7v7.rtf')
-                       
-system = psf.createSystem(params, nonbondedMethod=CutoffNonPeriodic,
+
+system = psf.createSystem(params, nonbondedMethod=CutoffPeriodic,
          nonbondedCutoff=1*nanometer, constraints=HBonds)
 
-integrator = LangevinIntegrator(310*kelvin,   
-                                1/picosecond, 
+integrator = LangevinIntegrator(310*kelvin,
+                                1/picosecond,
                                 0.002*picoseconds)
 
-simulation = Simulation(psf.topology, system, integrator)
+platform = Platform.getPlatformByName('CUDA')
+
+simulation = Simulation(psf.topology, system, integrator, platform)
 simulation.context.setPositions(pdb.positions)
 simulation

scripts/atefe.jl

@@ -6,24 +6,40 @@ end
 
 @time "loading" sim = OpenMMSimulation(py="data/atefe/mor.py", steps=10_000)
 
+@time "minimizing" OpenMM.minimize!(sim)
+
+
 NFEATURES = 3000
 NX = 1_000
 NK = 1
 NTRAIN = 30_000
-TESTRUN = false
 
-if TESTRUN
-    NX = 2
-    NK = 1
-end
+#begin
+#    @time xs = laggedtrajectory(sim, 1_000)
+#    @time JLD2.save("data/atafe/traj.jld2", "xs", xs)
+#end
+
+featurizer = OpenMM.FeaturesPairs(sample(OpenMM.calpha_pairs(sim.pysim), NFEATURES, replace=false))
+
+
+
 
-begin
-    @time xs = laggedtrajectory(sim, 1_000)
-    @time JLD2.save("data/atafe/traj.jld2", "xs", xs)
-end
 
-featurizer = FeaturesPairs(sample(OpenMM.calpha_pairs(sim.pysim), NFEATURES, replace=false))
 data = @time "sampling" trajectorydata_bursts(sim, NX, NK; featurizer)
-iso = Iso(data, opt=NesterovRegularized(), gpu=CUDA.has_cuda())
+iso = Iso(data, opt=NesterovRegularized())
+
+ISOKANN.save("atefe.iso", iso)
+
+@time "training" run!(iso, NTRAIN)
+
+
+ISOKANN.save("atefe.iso", iso)
+
+
+## scrapboard
 
-@time "training" run!(iso, 100)
+function lastsoluteindex(sim)
+    map(sim.pysim.topology.atoms()) do a
+        !(a.residue.name in ["POPC", "HOH", "SOD", "CLA"])
+    end |> findlast
+end

src/iso.jl

@@ -282,7 +282,7 @@ end
 Save the complete Iso object to a JLD2 file """
 function save(path::String, iso::Iso)
     iso = cpu(iso)
-    JLD2.save(path, "iso", iso)
+    JLD2.jldsave(path; iso)
 end
 
 """
@@ -294,9 +294,6 @@ An OpenMMSimulation will be reconstructed anew from the saved pdb file.
 """
 function load(path::String)
     iso = JLD2.load(path, "iso")
-    if iso.data.sim isa OpenMMSimulation
-        @warn "currently only the OpenMM `defaultsystem` is supported, any custom OpenMM simulations will not be reconstructed"
-    end
     return iso
 end
 

src/molutils.jl

@@ -139,7 +139,7 @@ end
 """
     save_trajectory(filename, coords::AbstractMatrix; top::String)
 
-save the trajectory given in `coords` to `filename` with the topology provided by the file `top`
+save the trajectory given in `coords` to `filename` with the topology provided by the file `top` using mdtraj.
 """
 function save_trajectory(filename, coords::AbstractMatrix; top::String)
     mdtraj = pyimport_conda("mdtraj", "mdtraj", "conda-forge")
@@ -172,7 +172,12 @@ function readchemfile(traj::Chemfiles.Trajectory, frames)
     xs = Array{Float32}(undef, length(Chemfiles.positions(frame)), length(frames))
     for (i, s) in enumerate(frames)
         Chemfiles.read_step!(traj, s - 1, frame)
-        xs[:, i] .= Chemfiles.positions(frame).data |> vec
+        try
+            xs[:, i] .= Chemfiles.positions(frame).data |> vec
+        catch e
+            i == size(xs, 2) || rethrow(e) ## hadle connect record which is not read properly by chemfiles
+            xs = xs[:, 1:end-1]
+        end
     end
     xs ./= 10 # convert from Angstrom to nm
     return xs
@@ -184,11 +189,15 @@ readchemfile(traj::Chemfiles.Trajectory, frames::Colon=:) =
 readchemfile(traj::Chemfiles.Trajectory, frame::Int) =
     readchemfile(traj, frame:frame) |> vec
 
+"""
+    writechemfile(filename, data::Array{<:Any,2}; source)
+
+Save the coordinates in `data` to `filename` with `source` as template using the Chemfiles library"""
 function writechemfile(filename, data::Array{<:Any,2}; source)
     trajectory = Chemfiles.Trajectory(source, 'r')
     try
         frame = Chemfiles.read(trajectory)
-        trajectory = Chemfiles.Trajectory(filename, 'w', uppercase(split(filename, ".")[end]))
+        trajectory = Chemfiles.Trajectory(filename, 'w')#, uppercase(split(filename, ".")[end]))
         for i in 1:size(data, 2)
             Chemfiles.positions(frame) .= reshape(data[:, i], 3, :) .* 10 # convert from nm to Angstrom
             write(trajectory, frame)

src/simulation.jl

@@ -29,6 +29,10 @@ end
 trajectory(sim::IsoSimulation, steps) = error("not implemented")
 laggedtrajectory(sim::IsoSimulation, nx) = error("not implemented")
 
+#laggedtrajectory(sim::OpenMMSimulation, n_lags, steps_per_lag=steps(sim); x0=getcoords(sim))
+#    trajectory(sim, n_lags * steps_per_lag; saveevery=steps_per_lag, x0)
+
+
 #TODO:
 
 
@@ -239,7 +243,7 @@ Simulate a single long trajectory of `steps` times the lagtime and start `nk` bu
 
 
 x0---x----x---
-    / |  / | 
+    / |  / |
     y y  y y
 """
 function trajectorydata_bursts(sim::IsoSimulation, steps, nk; x0=getcoords(sim), kwargs...)

test/runtests.jl

@@ -11,6 +11,10 @@ else
     @info "No functional GPU found. Skipping GPU tests"
 end
 
+@time @testset "workflow run" verbose = true begin
+    include("workflow.jl")
+end
+
 function with_possible_broken_domain(f)
     try
         r = f()