Acellera / htmd

HTMD: Programming Environment for Molecular Discovery
https://software.acellera.com/docs/latest/htmd/index.html
Other
254 stars 58 forks source link

Is it possible to insert a protein in a membrane with LPSA? #993

Closed phisanti closed 2 years ago

phisanti commented 3 years ago

I tried to make a membrane with Charmm-gui based on LPS on one leaflet. Then, I would like to insert a protein and generate the files for a simulation with openMM. Is htmd capable of building pdb/psf or prmtop/inpcrd files? If so, what would be the topologies/param arguments that I have to pass to the build function?

stefdoerr commented 3 years ago

Yes it can produce PSF. The PRMTOP maybe no because the lipids might not exist in AMBER. If you send me the membrane file I can take a look. Sometimes with weird lipids you need to add additional FF files when building

phisanti commented 3 years ago

Here it is a sample file: lps_membrane.zip

I tried to add some *.str files [toppar_all36_lipid_cardiolipin.str, toppar_all36_lipid_lps.str, toppar_all36_lipid_bacterial.str] to the topos argument list of charmm.build(). I got the files from the charmm-gui website where I prepared the membrane itself. However, I get all kinds of errors regarding the LPSA molecule and the PVCL2 molecule. Mainly, that LPS seems to be named as LPSPA in the str file (renaming the residues seems to help), but with PVCL2, I cannot move ahead.

phisanti commented 3 years ago

One last thing, in case you find something wrong with the file itself, the membrane was build following this tutorial: https://www.youtube.com/watch?v=KNEhbHbCYGM

phisanti commented 3 years ago

Any update on the issue or should I opt for a different alternative? My only goal is to remove some water from the bottom layer so that the simulation would run faster, and maybe change the charge state of the bottom layer lipids.

stefdoerr commented 3 years ago

Sorry I am a bit busy so I didn't have time to look at it yet. From my experience with 4-character lipid names you need to alias them to some other name.

I tried to load your PDB but the file seems to have weird Op and Cc atoms/elements. Could you give me the PSF as well which CHARMM-GUI produces with the right elements?

ATOM  42684 CC11 PVCL   31      89.501  45.189  81.597  1.00  0.00      MEMBCc 
ATOM  28887 OPA4 ECLI    1       1.489  36.010 102.883  1.00  0.00      LPSAOp  
phisanti commented 3 years ago

Don't worry, I am very thankful that you are indeed helping me. Here are the files:

pdb_psf.zip

The charmm-gui procedures, I do not know exactly what do you mean. If it is the charmm terminal commands, then I have not clue as I a a newcomer in the world of MD. If it is just the list of click options used in the charmm-gui website, I just followed this protocol: https://www.youtube.com/watch?v=KNEhbHbCYGM

Basically, it shows how to build a bacterial membrane with a protein inserted there. Here is the job retriever where I believe it's possible to re-check all the options: ID: 1366784860

Regarding the custom residues, I indeed tried to rename them, but I still got an error.

stefdoerr commented 3 years ago

Does this membrane contain every lipid in existence? :laughing:

Ok so you can nearly build it with the following commands:

mol = Molecule("step5_charmm2omm.pdb")
streams = ["str/carb/toppar_all36_carb_imlab.str", "str/lipid/toppar_all36_lipid_lps.str", "str/lipid/toppar_all36_lipid_bacterial.str"]
prms = charmm.defaultParam() + ["par/par_all36_carb.prm"]
topo = charmm.defaultTopo() + ["top/top_all36_carb.rtf"]
bmol = charmm.build(mol, topo=topo, param=prms, aliasresidues={"ECLIPA": "ECLI", "PVCL2": "PVCL"}, stream=streams)

It manages to build the system but then it crashes because it cannot load the PSF and PDB file at the same time. The reason for this is that the PSF contains 5-letter atom names which are not allowed in the PDB format so my check fails on it.

If you don't care about ionization you can use the above command and ignore the error. It will produce the correct PSF/PDB files.

For the error in reading PSF+PDB I will need to make a patch to moleculekit but this might take a few days.

phisanti commented 3 years ago

Thank you very much, I really appreciate it. I know the molecule is quite complicated, so this was a big challenging for me. Could I ask one last question out of curiosity, what is the difference between topology and stream?

stefdoerr commented 3 years ago

Stream files contain both topology and parameters in one file. You can even manually split them into prm and rtf.

phisanti commented 3 years ago

Sorry, it took me a bit long to test the solution. Now, it seems the problem with the naming is gone. I had to change the order of the elements in the dictionary tho. However, the build command is still throwing an error related with the naming. This time seems to be a discrepancy between psf vs pbd. Here it is the error: Same number of atoms but different atom names read from topology file [['/notebooks/test/structure.psf', 0]]

TopologyInconsistencyError: "Same number of atoms but different atom names read from topology file [['/notebooks/test/structure.psf', 0]]"

And here it is the code:


from htmd.ui import *
import os
import itertools
from moleculekit.util import opm
config(viewer='ngl')
config(njobs=6) # Use 4 cores
wd = os.getcwd()

mol = Molecule("step5_charmm2omm.pdb", validateElements=False)

streams = ["./storage/toppar/toppar_all36_carb_imlab.str", "./storage/toppar/toppar_all36_lipid_lps.str", "./storage/toppar/toppar_all36_lipid_bacterial.str"]
prms = charmm.defaultParam() + ["storage/toppar/par_all36_carb.prm"]
topo = charmm.defaultTopo() + ["storage/toppar/top_all36_carb.rtf", './storage/toppar/toppar_water_ions.str']
mol_s = autoSegment(mol)
bmol = charmm.build(mol_s, topo=topo, param=prms, aliasresidues={"ECLI": "ECLIPA", "PVCL": "PVCL2"}, stream=streams, outdir='test/')

I am missing something?

stefdoerr commented 3 years ago

No, it's the error I mentioned above. I can try to make a fix for it but for the moment you can find the built system (without ions though) in the build folder. I'll try to fix the issue with this error as soon as possible

stefdoerr commented 3 years ago

Ok so this version of molecule.py fixes the issue: https://github.com/Acellera/moleculekit/blob/master/moleculekit/molecule.py Currently I'm a bit busy trying to make the releasing working again so I cannot make an updated version of moleculekit but feel free to replace the molecule.py file in your miniconda with the one I linked. That should allow you to circumvent the issue.

phisanti commented 3 years ago

I have just tried it. Seems to work fine as far as I remove the ionization. Thank you very much, you have indeed gone the extra mile to help me!

stefdoerr commented 3 years ago

You are welcome :) With the updated molecule.py you can even ionize fine by the way. But it will take me maybe some days to release it officially