help on input file

[Nick-Mul]

Hi Joe,
I was wondering if I could get some pointers on input file? Tutorial files work fine but could get 1AKI from openmm tutorials to work.

I am using openmm/pdbfixer for the ligand prep. (version 8.2.0 and 1.11)

from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout
from openmm.app import PDBFile
from pdbfixer import PDBFixer

fixer = PDBFixer(filename='1AKI_clean.pdb')
fixer.findMissingResidues()
fixer.findNonstandardResidues()
fixer.replaceNonstandardResidues()
fixer.removeHeterogens(True)
fixer.findMissingAtoms()
fixer.addMissingAtoms()
fixer.addMissingHydrogens(7.0)
fixer.addSolvent(Vec3(5, 5, 5)*nanometer, positiveIon='Na+', ionicStrength=0.1*molar)
PDBFile.writeFile(fixer.topology, fixer.positions, open('output.pdb', 'w'))
pdb = PDBFile("output.pdb")

However I get this error

ERROR: LoadError: KeyError: key "H2" not found
Stacktrace:
 [1] getindex
   @ ./dict.jl:477 [inlined]
 [2] System(coord_file::String, force_field::MolecularForceField{Float64, Quantity{Float64, 𝐌 𝐍⁻¹, Unitful.FreeUnits{(g, mol⁻¹), 𝐌 𝐍⁻¹, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, Quantity{Float64, 𝐋² 𝐌 𝐍⁻¹ 𝐓⁻², Unitful.FreeUnits{(kJ, mol⁻¹), 𝐋² 𝐌 𝐍⁻¹ 𝐓⁻², nothing}}, Quantity{Float64, 𝐌 𝐍⁻¹ 𝐓⁻², Unitful.FreeUnits{(kJ, nm⁻², mol⁻¹), 𝐌 𝐍⁻¹ 𝐓⁻², nothing}}}; boundary::Nothing, velocities::Nothing, loggers::Tuple{}, units::Bool, array_type::Type{Array}, dist_cutoff::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, dist_neighbors::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, center_coords::Bool, neighbor_finder_type::Nothing, data::Nothing, implicit_solvent::Nothing, kappa::Quantity{Float64, 𝐋⁻¹, Unitful.FreeUnits{(nm⁻¹,), 𝐋⁻¹, nothing}}, rename_terminal_res::Bool)
   @ Molly ~/.julia/packages/Molly/86RTP/src/setup.jl:574
 [3] top-level scope

when trying to use the "output.csv"

using Molly

data_dir = joinpath(dirname(pathof(Molly)), "..", "data")

ff = MolecularForceField(
    joinpath(data_dir, "force_fields", "ff99SBildn.xml"),
    joinpath(data_dir, "force_fields", "tip3p_standard.xml"),
    joinpath(data_dir, "force_fields", "his.xml"),
)

sys = System(
     "output.pdb",
    ff; rename_terminal_res=false
)

[jgreener64]

You probably want rename_terminal_res=true for a standard protein.

If this doesn't work, could you attach the 1AKI_clean.pdb and output.pdb files?

[Nick-Mul]

Thanks for the reply, I've save the files here https://github.com/Nick-Mul/Julia_md.git

[jgreener64]

Unfortunately, we currently use atom names to match residue templates so atom names in the PDB file need to match the force field files exactly. This will hopefully change in future to use an approach like OpenMM, where residue templates are matched by bonding topology rather than atom name.

In this case you can change atom 2 to be called H1 and set rename_terminal_res=true but then you run into an issue with the HIS residue, which you might need to rename to one of the histidine variants in the ff99SBildn.xml file.

I just pushed a better error for this: 3beadb8.

[Nick-Mul]

Thanks I'll look into this, as I asking I've just began to play around with meta's UMA using ASE, can I use xyz format with molly?

[jgreener64]

You can convert an AtomsBase.jl system to a Molly System with something like System(ab_sys; force_units=u"kJ/Å", energy_units=u"kJ"), and an AtomsBase.jl system can be created from XYZ files using AtomsIO.jl. Let me know if you run into any issues.

[jgreener64]

System setup is much more flexible in Molly v0.23.1, hopefully this is fixed.