Build alchembed wrapper

[tomnewport]

Using alchembed could work in a couple of ways:

• Load a structure with topology
  • Use alchembed to embed a selection in the rest of the topology
  • OR
  • Use alchembed to embed a new structure and topology

a = Alchembed(structure="input.gro", topology="input.top")
a.add(structure="input2.gro", topology="input2.top")
a.add(selection="resname POPC")

The big challenge here is understanding GROMACS moltypes. Let's look at how the Protein moltype works:

[ moleculetype ]
; Name         Exclusions
Protein           1

[ atoms ]
    1    Qd     1   ALA    BB     1  1.0000 ; C
    2    P5     2   VAL    BB     2  0.0000 ; C
    3   AC2     2   VAL   SC1     3  0.0000 ; C
    4    P4     3   ALA    BB     4  0.0000 ; C
    5    P5     4   ASP    BB     5  0.0000 ; C
    6    Qa     4   ASP   SC1     6 -1.0000 ; C
    7    Nd     5   LYS    BB     7  0.0000 ; 1
    8    C3     5   LYS   SC1     8  0.0000 ; 1
    9    Qd     5   LYS   SC2     9  1.0000 ; 1
...
 2346    N0  1148   ARG   SC1  2346  0.0000 ; C
 2347    Qd  1148   ARG   SC2  2347  1.0000 ; C
 2348    Qa  1149   VAL    BB  2348 -1.0000 ; C
 2349   AC2  1149   VAL   SC1  2349  0.0000 ; C

[ bonds ]
; Backbone bonds
    1     2      1   0.35000  1250 ; ALA(C)-VAL(C)
    2     4      1   0.35000  1250 ; VAL(C)-ALA(C)
    4     5      1   0.35000  1250 ; ALA(C)-ASP(C)
...

[ constraints ]
    5     7      1   0.33000 ; ASP(C)-LYS(1)
    7    10      1   0.31000 ; LYS(1)-ALA(1)
   10    11      1   0.31000 ; ALA(1)-ASP(1)
   11    13      1   0.31000 ; ASP(1)-ASN(1)
   13    15      1   0.31000 ; ASN(1)-ALA(H)
   15    16      1   0.31000 ; ALA(H)-PHE(H)
   16    20      1   0.31000 ; PHE(H)-MET(H)
   20    22      1   0.31000 ; MET(H)-MET(H)
   22    24      1   0.31000 ; MET(H)-ILE(H)

[ angles ]
; Backbone angles
    1     2     4      2    127    20 ; ALA(C)-VAL(C)-ALA(C)
    2     4     5      2    127    20 ; VAL(C)-ALA(C)-ASP(C)
    4     5     7      2    127    20 ; ALA(C)-ASP(C)-LYS(1)
    5     7    10      2    127    20 ; ASP(C)-LYS(1)-ALA(1)
    7    10    11      2     96   700 ; LYS(1)-ALA(1)-ASP(1)
   10    11    13      2     96   700 ; ALA(1)-ASP(1)-ASN(1)
   11    13    15      2     96   700 ; ASP(1)-ASN(1)-ALA(H)
   13    15    16      2     96   700 ; ASN(1)-ALA(H)-PHE(H)
   15    16    20      2     96   700 ; ALA(H)-PHE(H)-MET(H)
   16    20    22      2     96   700 ; PHE(H)-MET(H)-MET(H)
   20    22    24      2     96   700 ; MET(H)-MET(H)-ILE(H)
   22    24    26      2     96   700 ; MET(H)-ILE(H)-CYS(H)

[ dihedrals ]
; Backbone dihedrals
    7    10    11    13      1   -120   400     1 ; LYS(1)-ALA(1)-ASP(1)-ASN(1)
   10    11    13    15      1   -120   400     1 ; ALA(1)-ASP(1)-ASN(1)-ALA(H)
   11    13    15    16      1   -120   400     1 ; ASP(1)-ASN(1)-ALA(H)-PHE(H)
   13    15    16    20      1   -120   400     1 ; ASN(1)-ALA(H)-PHE(H)-MET(H)
   15    16    20    22      1   -120   400     1 ; ALA(H)-PHE(H)-MET(H)-MET(H)
   16    20    22    24      1   -120   400     1 ; PHE(H)-MET(H)-MET(H)-ILE(H)
   20    22    24    26      1   -120   400     1 ; MET(H)-MET(H)-ILE(H)-CYS(H)
   22    24    26    28      1   -120   400     1 ; MET(H)-ILE(H)-CYS(H)-THR(H)

• The residue name seems pretty much irrelevant
• EVERYTHING ends up in this single itp file
• I'll need a decent ITP parser...
• Protein needs to be included

[tomnewport]

The final interface should work something like;

alchembed(simulation wrapper, input.gro, input.top, alchembed.mdp, molecule type, working directory )