Acellera / htmd

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

Patches and cofactors in charmm #3

Closed noeliaferruz closed 8 years ago

noeliaferruz commented 8 years ago

Hello guys,

I am not finding an easy way to add patches or parameterize cofactors whose parameters are included in charmm already (like phosphates (TP2), hemoglobine (HEME), INO1.. etc)

The problem comes from the fact that these RESIs are (only) included in the charmm_27prot_na ff. This has two consequences:

  1. (probably) many of the most standard patches and cofactors haven't been reparameterized since the release of C27.
  2. Adding patches using the lastest ff requieres sourcing the ff of the patch (RESI), which also depends on the charmm27na ff atomtypes. You can of course separate the RESI you want to apply along with the atomtypes it depends on (select those manually, and are different depending on the residue) in a separate file.

I was wondering if this workaround in 2), which is error-prone and requires a lot of manual intervention could be easily handled with htmd. Phosphates, hemoglobines and other molecules are present in everyday simulations. The other question is if these RESIs wouldnt be better described with parameterise.

Let me know what you think.

mj-harvey commented 8 years ago

Adreasing your last point, I'm not confident that parametrise in its current state can correctly parameterize a modified residue. The method's in principle capable of it, though, as evidenced by the gammp paper. Let me take a look over the weekend and see. M On 24 Mar 2016 15:47, "Noelia Ferruz" notifications@github.com wrote:

Hello guys,

I am not finding an easy way to add patches or parameterize cofactors whose parameters are included in charmm already (like phosphates (TP2), hemoglobine (HEME), INO1.. etc)

The problem comes from the fact that these RESIs are (only) included in the charmm_27prot_na ff. This has two consequences:

  1. (probably) many of the most standard patches and cofactors haven't been reparameterized since the release of C27.
  2. Adding patches using the lastest ff requieres sourcing the ff of the patch (RESI), which also depends on the charmm27na ff atomtypes. You can of course separate the RESI you want to apply along with the atomtypes it depends on (select those manually, and are different depending on the residue) in a separate file.

I was wondering if this workaround in 2), which is error-prone and requires a lot of manual intervention could be easily handled with htmd. Phosphates, hemoglobines and other molecules are present in everyday simulations. The other question is if these RESIs wouldnt be better described with parameterise.

Let me know what you think.

— You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3

giadefa commented 8 years ago

We have done it already On 24 Mar 2016 22:35, "M J Harvey" notifications@github.com wrote:

Adreasing your last point, I'm not confident that parametrise in its current state can correctly parameterize a modified residue. The method's in principle capable of it, though, as evidenced by the gammp paper. Let me take a look over the weekend and see. M On 24 Mar 2016 15:47, "Noelia Ferruz" notifications@github.com wrote:

Hello guys,

I am not finding an easy way to add patches or parameterize cofactors whose parameters are included in charmm already (like phosphates (TP2), hemoglobine (HEME), INO1.. etc)

The problem comes from the fact that these RESIs are (only) included in the charmm_27prot_na ff. This has two consequences:

  1. (probably) many of the most standard patches and cofactors haven't been reparameterized since the release of C27.
  2. Adding patches using the lastest ff requieres sourcing the ff of the patch (RESI), which also depends on the charmm27na ff atomtypes. You can of course separate the RESI you want to apply along with the atomtypes it depends on (select those manually, and are different depending on the residue) in a separate file.

I was wondering if this workaround in 2), which is error-prone and requires a lot of manual intervention could be easily handled with htmd. Phosphates, hemoglobines and other molecules are present in everyday simulations. The other question is if these RESIs wouldnt be better described with parameterise.

Let me know what you think.

— You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3

— You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3#issuecomment-201034895

mj-harvey commented 8 years ago

Which 'it'? On 25 Mar 2016 09:23, "giadefa" notifications@github.com wrote:

We have done it already On 24 Mar 2016 22:35, "M J Harvey" notifications@github.com wrote:

Adreasing your last point, I'm not confident that parametrise in its current state can correctly parameterize a modified residue. The method's in principle capable of it, though, as evidenced by the gammp paper. Let me take a look over the weekend and see. M On 24 Mar 2016 15:47, "Noelia Ferruz" notifications@github.com wrote:

Hello guys,

I am not finding an easy way to add patches or parameterize cofactors whose parameters are included in charmm already (like phosphates (TP2), hemoglobine (HEME), INO1.. etc)

The problem comes from the fact that these RESIs are (only) included in the charmm_27prot_na ff. This has two consequences:

  1. (probably) many of the most standard patches and cofactors haven't been reparameterized since the release of C27.
  2. Adding patches using the lastest ff requieres sourcing the ff of the patch (RESI), which also depends on the charmm27na ff atomtypes. You can of course separate the RESI you want to apply along with the atomtypes it depends on (select those manually, and are different depending on the residue) in a separate file.

I was wondering if this workaround in 2), which is error-prone and requires a lot of manual intervention could be easily handled with htmd. Phosphates, hemoglobines and other molecules are present in everyday simulations. The other question is if these RESIs wouldnt be better described with parameterise.

Let me know what you think.

— You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3

— You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3#issuecomment-201034895

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3#issuecomment-201196570

stefdoerr commented 8 years ago

I added CHARMM stream files in the new multiscalelab/htmd. Do charmm.listFiles() and check them out since they have heme parametrizations and other interesting stuff.

noeliaferruz commented 8 years ago

Fantastic!! I just saw it. How does htmd read str files? (As they contain both prm and rtf in one single file, do we have to split them before building?)

stefdoerr commented 8 years ago

Yes, it splits them into some temp files and then copies them into your build directory. I have tested only the splitting though so no guarantees the rest of the pipeline works. Do tell me though ;) Should be clear enough since the files should appear in the build dir.

On Tue, Apr 5, 2016 at 6:17 PM, Noelia Ferruz notifications@github.com wrote:

Fantastic!! I just saw it. How does htmd read str files? (As they contain both prm and rtf in one single file, do we have to split them before building?)

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3#issuecomment-205878905

giadefa commented 8 years ago

pratically they are usable directly by charmmBuilt

On 5 April 2016 at 19:18, Stefan notifications@github.com wrote:

Yes, it splits them into some temp files and then copies them into your build directory. I have tested only the splitting though so no guarantees the rest of the pipeline works. Do tell me though ;) Should be clear enough since the files should appear in the build dir.

On Tue, Apr 5, 2016 at 6:17 PM, Noelia Ferruz notifications@github.com wrote:

Fantastic!! I just saw it. How does htmd read str files? (As they contain both prm and rtf in one single file, do we have to split them before building?)

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3#issuecomment-205878905

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3#issuecomment-205905535

http://www.acellera.com

   <https://twitter.com/acellera>

https://www.youtube.com/user/acelleracom https://www.linkedin.com/company/2133167?trk=tyah&trkInfo=clickedVertical%3Acompany%2CclickedEntityId%3A2133167%2Cidx%3A2-1-2%2CtarId%3A1448018583204%2Ctas%3Aacellera https://www.acellera.com/md-simulation-blog-news/ http://is.gd/1eXkbS

giadefa commented 8 years ago

Can this be closed @noeliaferruz ?

noeliaferruz commented 8 years ago

Hola!

Sorry, I hadn't quite tested it. I have just tried to build the case with the phosphate patch. It's nothing important, I'm just trying to see if I can actually build a system with phosphates. These are the lines I'm using:

./1.build.py 5hx8.pdb 8 0 [Where 1.build.py looks like:]

system = Molecule(sys.argv[1])
system.filter("chain A and (protein or water)")
systemCaps = autoSegment(system, "protein", "P")
systemCaps.mutateResidue('resname PTR', 'TYR')
s = float(sys.argv[2])
d = maxDistance( systemCaps, "all" )
d=d+s
solvated=solvate(systemCaps, minmax=[[-d, -d, -d], [d,d,d]])
topos= [ "top/top_all36_prot.rtf",  "top/top_water_ions.rtf"]
params=[ "par/par_all36_prot_mod.prm", "par/par_water_ions.prm", "str/na/toppar_all36_na_model.str"]
built= charmm.build( solvated, topo=topos, param=params, saltconc=float(sys.argv[3]), outdir='build', patches='patch TP2 A:1035')

---------------------------------------------------------------------------------------------------------------
Videos from the HTMD2015 workshops are available on the Acellera youtube channel: https://www.youtube.com/user/acelleralive

You are on the latest HTMD version (1.0.16).
Solvating: 100% (27/27) [#######################################################################################################################################################################################################] eta 00:01 |
Traceback (most recent call last):
  File "./build.py", line 25, in <module>
    built= charmm.build( solvated, topo=topos, param=params, saltconc=float(sys.argv[3]), outdir='build', patches='patch TP2 A:1035')
  File "/nfs/grid/software/hpcc/apps/Linux-x86_64-RHEL6/acellera/current/python/lib/python3.5/site-packages/htmd/builder/charmm.py", line 127, in build
    caps = _defaultCaps(mol)
  File "/nfs/grid/software/hpcc/apps/Linux-x86_64-RHEL6/acellera/current/python/lib/python3.5/site-packages/htmd/builder/charmm.py", line 412, in _defaultCaps
    raise NameError('Segments {} contain both protein and non-protein atoms. Please assign separate segments to them.'.format(intersection))
NameError: Segments ['P2'] contain both protein and non-protein atoms. Please assign separate segments to them.

Visual inspection reveals that the phosphate is still there. I know @stefdoerr fixed this exact issue in Molecule.py some weeks ago, so I don't know. ptr

So, no don't close it yet, I'd like to find a way to add phosphates and other cofactors!

Thanks, Noelia

stefdoerr commented 8 years ago

No, I think the issue was that VMD atomselect was recognizing phosphates as "protein". Seems to be the case here at least. I don't know a fix for that, other than making your autoSegment atomselection fancier to exclude the phosphates.

noeliaferruz commented 8 years ago

But its mutateResidue who's not removing the entire sidechain of PTR, no?

stefdoerr commented 8 years ago

Ah if that's the case (which should be fixed), just swap the order. First do stuff like mutating etc and then autoSegment. Even if the phosphates are not removed (which they should), they will not be bonded any more to the protein so then the autoSegment should work correct.

stefdoerr commented 8 years ago

You are sure the phosphates belong to the old PTR residues?

noeliaferruz commented 8 years ago

Ok, I'll do that. Yes, mutateResidue keeps P, OP1 and OP2 from the original PTR residue but renames them to TYR.

stefdoerr commented 8 years ago

Wow. That means that the VMD atomselection recognizes those atoms as backbone. Can you verify that? Because I remove all atoms with 'resname PTR and not backbone'. I somehow assumed this would cover all weird atoms attached to a residue.

noeliaferruz commented 8 years ago

hahah yes! Exactly, resname PTR and not backbone removes the sidechain, and one oxygen of the phosphate. It keeps the backbone and the phosphate, and two of its three oxygen. xD. amazing. screenshot-short

j3mdamas commented 8 years ago

Maybe this could be sent to the VMD mailing list?

On Thu, Apr 21, 2016 at 6:08 PM, Noelia Ferruz notifications@github.com wrote:

hahah yes! Exactly, resname PTR and not backbone removes the sidechain, and one oxygen of the phosphate. It keeps the backbone and the phosphate, and two of its three oxygen. xD. amazing. [image: screenshot] https://cloud.githubusercontent.com/assets/12035024/14716171/c4ce9e34-07b9-11e6-9523-134c263ab0e4.png

— You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3#issuecomment-212989284

stefdoerr commented 8 years ago

Kill me... :P Well that is somehow beyond my capabilities, haha. If anyone can find a generalized atomselect (that works on all residues) and removes them I will implement it.

stefdoerr commented 8 years ago

Yeah we could try our luck there. @noeliaferruz can you send the residue (and maybe one or two before and after for "context") by mail so that I can ask them?

noeliaferruz commented 8 years ago

The backbone atoms are always N CA C O, no?

j3mdamas commented 8 years ago

I follow their mailing list, it's quite active.

On Thu, Apr 21, 2016 at 6:13 PM, Stefan notifications@github.com wrote:

Yeah we could try our luck there. @noeliaferruz https://github.com/noeliaferruz can you send the residue (and maybe one or two before and after for "context") by mail so that I can ask them?

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3#issuecomment-212990710

stefdoerr commented 8 years ago

Hm Noelia, the protein structure was it by any chance written by HTMD at some point? We had a small very weird issue with the PDB writer. Can you send me like the original original pdb structure of this residue?

stefdoerr commented 8 years ago

Apparently in the PDB format there is a significance on where you put atom names. So carbon alphas should be written in columns 13-16 as "..CA" but calcium should be ".CA." where dots indicate spaces. So if VMD assigns significance to that it could have problems.

noeliaferruz commented 8 years ago

The original pdb: 5hx8.pdb 5hx8.txt The one just wrote with htmd: protein.txt

\ But I'm using the first for the selections in the picture **

stefdoerr commented 8 years ago

Ok thanks. It's not our error. Happens also if you get the PDB from the PDB website directly. Will ask the VMD mailing list.

noeliaferruz commented 8 years ago

Would in the meanwhileresname PTR and not name N CA C O work?

stefdoerr commented 8 years ago

yes, do it manually. Do mol.remove('resname PTR and not name N CA C O') and then mol.set('resname', 'TYR', sel='resname PTR')

noeliaferruz commented 8 years ago

yep. Continuing...

stefdoerr commented 8 years ago

Ok apparently backbone in VMD works for both protein and nucleic acids. In this case it thought that your P, PO etc were nucleic acid backbone. Seems like I can substitute the backbone selection with C CA N O without loss of generality so I will do that.

j3mdamas commented 8 years ago

Maybe create a new selection protein-backbone instead of substitution?

stefdoerr commented 8 years ago

There is no such atomselection in VMD. I don't think I can invent atomselection shortcuts either.

stefdoerr commented 8 years ago

Not without modifying the VMD source code at least which I don't see as a real need. It would end up doing the same thing anyways.

j3mdamas commented 8 years ago

forget about it. you were talking about your specific application not HTMD :)

On Fri, Apr 22, 2016 at 11:43 AM, Stefan notifications@github.com wrote:

Not without modifying the VMD source code at least which I don't see as a real need. It would end up doing the same thing anyways.

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/Acellera/htmd/issues/3#issuecomment-213359079

giadefa commented 8 years ago

let's talk about it when Noelia is in BCN

j3mdamas commented 8 years ago

@noeliaferruz, let's see this on monday.

j3mdamas commented 8 years ago

we fixed this when Noelia was here, but now there's only a general problem with patching that's getting solved (#97)

stefdoerr commented 8 years ago

So we can close this then. All issues are fixed for this as far as I can see (atom naming, backbone selections)

j3mdamas commented 8 years ago

yes, only issue 97 left.