How to get the Martini CG of an AA polymer?

[pranoy-ray]

I am trying to CG a molecule of polystyrene (styrene with a degree of polymerization: 20) while I have the input .gro file, the top file, and I am able to generate the AA ff file, the mapping file, and the CG ff file. I have added this polymer as a new reidue in my local martini folder. The number and type of atoms match everywhere (322) atoms for AA and 80 beads for CG. Vermouth-martinize2 runs without any error. However, it does not give me any results even after multiple hours. Is this an internal issue, or am I doing something wrong? This is where it gets frozen:

    INFO - step - Guessing the bonds.
    INFO - general - 1 molecules after guessing bonds
 WARNING - general - Mutation "cter" not found. Check target resid!
 WARNING - general - Mutation "nter" not found. Check target resid!
    INFO - step - Repairing the graph.
....

And the command I am using to get here is:
martinize2 -ff-dir force_fields -map-dir mappings -f PS.gro -x PSm.pdb -o topol.top -maxwarn 3

[pckroon]

That's indeed not really desired behaviour, for a polymer this small martinize2 shouldn't take much more than a few seconds.. There are a few steps in the pipeline whose runtime scales exponentially with the size of the residue.

• Does your input gro file separate the residues by residue number? We don't test the gro-reading very extensively, you could try converting it to a PDB, and trying with that. If that does work you found a bug I'd be interested in.
• Could you try with -v (or -vv)? That might shed some more light. In particular it will tell when residues clash, and erroneous bonds get created, which can also cause these symptoms.
• Finally, if this doesn't help, could you post your atomistic styrene definition? It may be there's missing edges.

[pranoy-ray]

Hey @pckroon, I want to clarify the fact that I am CGing my polymer (322 atoms: 20-polystyrene) as a single residue, is that a reason why it is scaling exponentially? What do you suggest I do in that case, martinize monomers and then connect them manually?

Right now, I have the ITP for 322 atoms and the GRO for the same 322 atoms, all representing 20 polystyrene in the form of a single residue. Should I relabel the styrene monomers in these GROs and ITPs to represent the building blocks of styrene that have been used? The issue with this question is that the polymer's head and tail atoms have satisfied valency by a hydrogen atom, which is removed for polymerization. The monomer unit has 16 atoms. The polymer unit is represented by 16x20 = 320 atoms in total, but to satisfy the valency of the end atoms, we add two hydrogens, bringing up the total number to 322.

Here is a plot showing how the degree of polymerization currently scales. Beyond ~112 atoms, the martinize2 time becomes infinity.

I also used the -v command as you recommended, and that had this output, and no progress from thereon:

    INFO - step - vermouth - Guessing the bonds.
    INFO - general - vermouth.processors.make_bonds - 1 molecules after guessing bonds
   DEBUG - step - vermouth - Annotating required mutations and modifications.
 WARNING - general - vermouth.processors.annotate_mut_mod - Mutation "cter" not found. Check target resid!
 WARNING - general - vermouth.processors.annotate_mut_mod - Mutation "nter" not found. Check target resid!
    INFO - step - vermouth - Repairing the graph.
   DEBUG - step - vermouth.processors.repair_graph - Making reference graph
   DEBUG - step - vermouth.processors.repair_graph - Matching residue -PS1 to its reference

[pckroon]

Yes, that is kind of expected. You should do it the same way we do for proteins, basically. Define a styrene monomer, and feed it the conformation of a polymer where each monomer has a different residue number.
As you note, the termini will be different: you'll need to take a look at modifications (c.f. proteins!). In addition you'll need to define the parameters used to link the monomers together. We call those links.
Did you take a look at the documentation on how to add residues, and the issue I referenced 2 weeks ago?

[pranoy-ray]

@pckroon, say I want 19 links for the 20 monomers. I am beading styrene into 4 CG beads. and the CG4 of every monomer needs to be connected to another monomer's CG4, and so on.

S stands for styrene and 20 stands for 20 monomers. Should the [ links ] part of my CF-ff file look like this:

[ link ]
[ bonds ]
CG4 {"resname": "S"} +CG4 {"resname": "S"} 20 0.35 4000

I am really confused and I need some clarification here. Yes, I read the documentation you shared earlier.

[pckroon]

Almost. You don't need to specify how often the link needs to be applied, that will follow automatically from the AA structure.

Bonds require 2 atoms, so: CG4 {"resname": "S"} is the first bead (beadname CG4, resname S, this residue); and +CG4 {"resname": "S"} (beadname CG4, resname S, next residue [meaning, residue number is one higher]) is the second bead. The following 3 numbers get put into the output itp file verbatim. If memory serves: the first number is the functional type, the second is the bond length, and the last is the spring constant. The exact meaning of those you'll need to find in the gromacs manual. If the link you put here gets applied it will create (for example):

[bonds]
1 5 20 0.35 4000
5 9 20 0.35 4000
9 13 20 0.35 4000
etc...

[pranoy-ray]

@pckroon My input .gro file looks like this:

PS20
320
    1S     C1     1   4.175   4.874   7.235  0.0000  0.0000  0.0000
    1S     C2     2   4.178   4.809   7.093  0.0000  0.0000  0.0000
    1S     C3     3   4.356   4.884   7.413  0.0000  0.0000  0.0000
    1S     C4     4   4.438   4.823   7.508  0.0000  0.0000  0.0000
    .....................
    9S     H6   143   5.213   4.651   5.195  0.0000  0.0000  0.0000
    9S     H7   144   5.036   4.478   5.179  0.0000  0.0000  0.0000
    9S     H8   145   4.798   4.544   5.165  0.0000  0.0000  0.0000
    10S     C1   146   4.701   4.796   4.887  0.0000  0.0000  0.0000
    10S     C2   147   4.632   4.750   4.749  0.0000  0.0000  0.0000
    10S     C3   148   4.948   4.859   4.895  0.0000  0.0000  0.0000
    10S     C4   149   5.083   4.825   4.898  0.0000  0.0000  0.0000
    .....................
    20S     H5   318   4.469   4.851   1.941  0.0000  0.0000  0.0000
    20S     H6   319   4.549   4.622   1.888  0.0000  0.0000  0.0000
    20S     H7   320   4.446   4.426   1.999  0.0000  0.0000  0.0000
    20S     H8   321   4.267   4.458   2.167  0.0000  0.0000  0.0000
9.455400000000001 9.455400000000001 9.455400000000001

The corresponding S.ff for charmm looks like this:

[ moleculetype ]
;Name        nrexcl 
S            3     

[ atoms ]
;         nr         type        resnr       residu         atom         cgnr       charge 
          1           c3            1        S              C1             1 0.0799680400 
          2           c3            1        S              C2             2 -0.088620950 
          3           ca            1        S              C3             3 -0.195149890 
          4           ca            1        S              C4             4 -0.140383920 
          5           ca            1        S              C5             5 -0.140026920 
          .....................
         15           ha            1        S              H7            15 0.1384930800 
         16           ha            1        S              H8            16 0.1470680800 


[ bonds ]
;         ai           aj       funct        b0[nm] kb[kJ mol^-1 nm^-2] 
           3            4            1 0.1387000000 400325.11999 
           4            5            1 0.1387000000 400325.11999 
           5            6            1 0.1387000000 400325.11999 
            .....................
          14            5            1 0.1087000000 288110.23999 
          15            6            1 0.1087000000 288110.23999 
          16            7            1 0.1087000000 288110.23999

The corresponding S.ff for martini3001 looks like this:

[ moleculetype ]
S  3
;
[ atoms ]
;  nr    type       resnr  residu    atom    cgnr        charge    mass 
   1    TC5            1  S        CG1       1     0.0   45
   2    TC5            1  S        CG2       2     0.0   45
   3    TC5            1  S        CG3       3     0.0   45
   4    TC3            1  S        CG4       4     0.0   45
;
[ bonds ]
;  ai  aj   funct
   1     4   1     0.27 8000

[ constraints ]
;  ai  aj  funct length
   2     3    1     0.290 {"ifndef": "FLEXIBLE"}
   3     1    1     0.290 {"ifndef": "FLEXIBLE"}
   1     2    1     0.290 {"ifndef": "FLEXIBLE"}

[ bonds ]
;  ai  aj  funct length
   2     3    1     0.290 10000 {"ifdef": "FLEXIBLE"}
   3     1    1     0.290 10000 {"ifdef": "FLEXIBLE"}
   1     2    1     0.290 10000 {"ifdef": "FLEXIBLE"}

;
[ angles ]
;  ai  aj  ak  funct
   4     1     2    1   136.000000  100.000000
   4     1     3    1   136.000000  100.000000

[ link ]
resname "S"
[ bonds ]
CG4   +CG4    1     0.27 8000

[ link ]
resname "S"
[ angles ]
CG4   +CG1     +CG2    1   136.000000  100.000000
CG4   +CG1     +CG3    1   136.000000  100.000000
CG1   CG4      +CG1    1   120.000000   25.000000
CG4   +CG1     +CG4    1    52.000000  550.000000

[pranoy-ray]

@pckroon Output log for the above:

    INFO - step - Guessing the bonds.
 WARNING - unknown-residue - Can't add bonds based on atom names for residue -0S1 because 'Residue 0S is not known to force field charmm'. Falling back to distance criteria.
 WARNING - unknown-residue - Can't add bonds based on atom names for residue -1S1 because 'Residue 1S is not known to force field charmm'. Falling back to distance criteria.
 WARNING - unknown-residue - Can't add bonds based on atom names for residue -2S1 because 'Residue 2S is not known to force field charmm'. Falling back to distance criteria.
 WARNING - unknown-residue - Can't add bonds based on atom names for residue -3S1 because 'Residue 3S is not known to force field charmm'. Falling back to distance criteria.
 WARNING - unknown-residue - Can't add bonds based on atom names for residue -4S1 because 'Residue 4S is not known to force field charmm'. Falling back to distance criteria.
 WARNING - unknown-residue - Can't add bonds based on atom names for residue -5S1 because 'Residue 5S is not known to force field charmm'. Falling back to distance criteria.
 WARNING - unknown-residue - Can't add bonds based on atom names for residue -6S1 because 'Residue 6S is not known to force field charmm'. Falling back to distance criteria.
 WARNING - unknown-residue - Can't add bonds based on atom names for residue -7S1 because 'Residue 7S is not known to force field charmm'. Falling back to distance criteria.
 WARNING - unknown-residue - Can't add bonds based on atom names for residue -8S1 because 'Residue 8S is not known to force field charmm'. Falling back to distance criteria.
 WARNING - unknown-residue - Can't add bonds based on atom names for residue -9S1 because 'Residue 9S is not known to force field charmm'. Falling back to distance criteria.
 WARNING - unknown-residue - Can't add bonds based on atom names for residue -0S2 because 'Residue 0S is not known to force field charmm'. Falling back to distance criteria.
    INFO - general - 1 molecules after guessing bonds
 WARNING - general - Mutation "cter" not found. Check target resid!
 WARNING - general - Mutation "nter" not found. Check target resid!
    INFO - step - Repairing the graph.
 WARNING - unknown-residue - Cannot recognize residue '0S' in  molecule 0. Deleting the molecule.
    INFO - step - Dealing with modifications.
    INFO - step - Read input.
    INFO - step - Creating the graph at the target resolution.
    INFO - step - Averaging the coordinates.
    INFO - step - Applying the links.
    INFO - step - Placing the charge dummies.
    INFO - step - Sorting atomids
    INFO - step - Writing output.
Traceback (most recent call last):
  File "/Users/pranoy/opt/anaconda3/envs/venv/bin/martinize2", line 1114, in <module>
    entry()
  File "/Users/pranoy/opt/anaconda3/envs/venv/bin/martinize2", line 1085, in entry
    write_gmx_topology(system,
  File "/Users/pranoy/opt/anaconda3/envs/venv/lib/python3.9/site-packages/vermouth/gmx/topology.py", line 148, in write_gmx_topology
    raise ValueError("No molecule in the system. Nothing to write.")
ValueError: No molecule in the system. Nothing to write.
Time Taken: 5.22s

I understand that the atoms 10S onwards are not identified for whatever reason; what can I do to mitigate this? Seems like a formatting issue now.

[pckroon]

How did you generate your gro file? Either that is mildly misformatted (note the inconsistent outlining of the lines below), or there is something wrong with our gro parser :)

    9S     H8   145   4.798   4.544   5.165  0.0000  0.0000  0.0000
    10S     C1   146   4.701   4.796   4.887  0.0000  0.0000  0.0000

Martinize2 does tell you what's going wrong: WARNING - unknown-residue - Can't add bonds based on atom names for residue -0S1 because 'Residue 0S is not known to force field charmm'. Falling back to distance criteria.

At-a-glance your ff file looks sane. You'll still need a mapping (.map or .mapping) file as well, but martinize2 will complain when we get there.

[pranoy-ray]

I do have a mapping file, and to note again upto 9PS there is no issue and the martinized structure for 9-polystyrene is generating correctly. About 10S and onwards I have tried moving the 10S column one space to the left as wel as tried introducing a space between 9S and the others above, however the gro parser fails both the times. Yes, I have converted into a pdb and passed it through martinize2, it didn't work out! @pckroon

[pckroon]

From the gromacs documentation on the .gro file format:

Columns contain the following information (from left to right):

• residue number (5 positions, integer)
• residue name (5 characters)
• ...

Since your gro file is misformatted I'm not surprised converting it to pdb doesn't help: the conversion will propagate the same problem.

[pranoy-ray]

@pckroon Got it, that was the error as you mentioned, it worked!

[pranoy-ray]

@pckroon How do I extend the CG version of PS20 onto a box containing 100 polymer chains in GROMACS? Keeping the intra-chain connections intact.

[pckroon]

Martinize2 does nothing (very little) with coordinates. If you provide the gro/pdb with the coordinates, Martinize2 will generate your itp files

[pckroon]

But also take a look at polyply: https://github.com/marrink-lab/polyply_1.0

[pranoy-ray]

@pckroon Are there any GROMACS commands like insert-molecule that can help?

[pckroon]

Maybe. I'm more than happy to help people with vermouth/martinize2, but I really don't have the time for anything outside those.