【How to make a block copolymer with nonstandard PGLeu monomers】

[kurokawaikki]

Dear developers,

I am new to the CG model simulation. I was using all atom MD simulation with Desmond. Currently, I would like to study the process of endosome escape, but all atom simulation is too slow and fail to serve the purpose.
My model basically contains 2 parts. The first one is the endosome lipid bilayer (I obtained from this reference: A molecular view on the escape of lipoplexed DNA from the endosome). The second part is my polymers. The polymer contains 2 parts: one is PEG with leucine in side chain and the other one is PEG tail (It is shown in the following figure).
https://i.ibb.co/nkQRSNL/Untitled.png

I would like to know how to parametrize it and generate a topology file that can be used in the martini force field.
Because polyply is powerful in making the CG model, I would like to build my polymer model with it. However, my block copolymer contains nonstandard monomer. How could I obtain and parameterize the monomer of PEG with Leu in side chain?
Finally, is there
a way to import the nonstandard monomer model into polyply and use it and PEG monomer in the polyply to make a PEG-PGLeu copolymer with polyply?
Could you give me some guides and hints on how to walk through the process?

Thank you very much for your great help and your time.

[fgrunewald]

@kurokawaikki we are always happy to help with it. I try my best walking you through this. First you need to choose Martini2 or Martini3 as the force-field. Please note that the current PEG martini3 does not work well in combination with membranes, so using Martini2 might be better. However, the LEU is pretty bad in Martini2. We are currently working on improving PEG for Martini3. I'd expect and updated version to be released at most in 4-5 weeks from now. As soon as there is some working option for Martini3 PEG I can share it with you, if you like. The rest of the answer is written assuming you want to use Martini3.

Please also note that you should verify the properties of PGLeu independently from the PEG either against atomistic simulation data or experimental data. But don't worry about that now first we obtain some parameters.

I understand you already have an all-atom simulation of your polymer or at least a fragment of it? If that is the case, what you need to do is derive the bonded parameters for the missing monomer. In particular you need the bond and an angle between the main backbone and the side chain as well as bonds and angles between PEG and the modified PEG bead. The values you can obtain by mapping the trajectory essentially following this tutorial.

It might also be useful for you to read through his issue, which deals with a similar question.

To generate a preliminary itp file you need a .ff input file for the new monomer. Let's give it the residue name 'PGLeu'. To clarify I would treat the PG-LEU unit as one monomer i.e. 3 beads following the Martini3 mapping. The syntax of the .ff file is documented to some extent here. For you case I've written a template file as shown below. Note the places with '<some_text>' you need to replace with the parameters obtained from the mapping.

[ moleculetype ]
; molname       nrexcl
PGLeu                1

[ atoms ]
;id type resnr residu atom cgnr   charge
 1   SN2r 1     PGLeu    PG      1       0    ; perhaps this bead type needs to be changed
 2   P2   1     PGLeu    BB       2      0    ; perhaps this bead type needs to be changed
 3   SC2  1     PGLeu    SC1     3      0

[ bonds ]
 2    3    1       0.363     7500
[ angles ]
1    2     3    1   <angle> <force_constant>

;;;;; linking the PGLeu together
[ link ]
resname "PGLeu"
[ bonds ]
PG  >PG   1  <length> <force-const> ; this introduces a bond between to bonded PELeu monomers connecting the PG beads

[ link ]
resname "PGLeu"
[ angles ]
PG  >PG  >>PG  1  <angle> <force_constant> {"group": "BB angles"}

[ link ]
resname "PGLeu"
[ angles ]
PG  >PG  >SC1  1  <angle> <force_constant> {"group": "BB-SC angles"} ; angle between side-chain and backbone

;;;;; link between PEG and PGLEU

[ link ]
resname "PGLeu|PEO"
[ bonds ]
EC  >PG   1  <length> <force-const> {"comment": "PEG-PGLeu bond"}
[ angles ]
EC   EC  >PG  1 <angle> <force_constant> {"group": "PEG-PGLeu angles"}
EC >PG >>PG 1 <angle> <force_constant> {"group": "PEG-PGLeu angles"}
EC >PG >SC1 1 <angle> <force_constant> {"group": "PEG-PGLeu angles"}

To generate the itp file save this file as PGLeu.ff and the use it with polyply as follows:

polyply gen_params  -lib martini3 -f PGLeu.ff -seq PEO:10 PGLeu:10 -o init.itp -name PEO-b-PGLeu 

Now all you need to do is fill in the missing parameters based on the atomistic simulations. I would also recommend adjusting the bead-types of the PG fragment as well as the LEU fragment. However, once you have obtained the bonded parameters we can take that step. Simply reply to this issue.

Please try of you can generate the itp and let me know if you have more questions.

[kurokawaikki]

Thank you very much for your guidance.
I am preparing to run the all atom MD of my PEG-PGLeu polymer.
My understanding so far is that I use the MD data of the fragment of polymer (one PEG connected to one PGLeu) and map to the CG model. Calculate the distribution for bond length, angle and dihedrals for each atoms. However, could you give me some hint on how to derive force constant? Or suggest some materials for me to read.
Finally, I can take the derived parameter and fill them into the .ff file and generate the initial itp for my polymer.

Thank you again for your great help.

[kurokawaikki]

I have now progressed to the step of parameterization of PEGLeu monomer.
I now got the atomistic results that are ready for mapping to CG model.

My question is how can I map my molecule to CG model?
In the template that you gave me, I can found that 3 type of beads were used.
However, I am wondering how do you choose the beads for a specific molecule.
Additionally, I cannot find the type SN2r in the martini 3 building block table.
I have checked the tutorial that you gave (writing .ff input files). However, I still not fully understand how to write each definition correctly, if I would like to adapt to a new monomer (for example PEGLeu to PEGGly).

In my PEGLeu case, if I choose 4-1 mapping, I will need three beads. However, I have no idea on how to choose the bead type from the building block. Most important of all, what is the relation between the CG model that I fit and the beads that you chose here?
Can you give me more explanation on chose of bead type and atom?

[ atoms ]
How can I know what to input in the atom? Is there any tutorial on writing this?
;id type resnr residu atom cgnr charge
1 SN2r 1 PGLeu PG 1 0 ; perhaps this bead type needs to be changed
2 P2 1 PGLeu BB 2 0 ; perhaps this bead type needs to be changed
3 SC2 1 PGLeu SC1 3 0

I am completely new to CG model and gromacs and I am sorry for lots of questions.
Thank you very much for your great help.

[fgrunewald]

Hi,

Let me try to answer your questions:

My question is how can I map my molecule to CG model?

The best way I would say is by similarity to what already exists. For PGLeu I've checked how many beads the Leu part of the amino-acid has and then simply added one bead for the backbone as that is most similar to PEG. We also have some more defined mapping rules in the SI of the Martini3 paper (see page 36)

In my PEGLeu case, if I choose 4-1 mapping, I will need three beads. However, 
I have no idea on how to choose the bead type from the building block. 
Most important of all, what is the relation between the CG model that I fit and the beads that you chose here?
Can you give me more explanation on chose of bead type and atom?

I would follow the following 3 guiding principles:

• look for similar molecules already defined in Martini and take the bead-types of those fragments. For example, LEU is already defined for Martini proteins so the side-chain should have the same bead-types
• otherwise look for the fragment in the Martini bead-type fragment table and choose the closest fragment.
• the third way is to choose the bead-type that matches best some properties from experiment or AA simulations.

In principle, you should always validate your bead assignment by comparing to experimental and/or simulation data of either the polymer or the fragment. For example, if you have an atomistic simulation of your polymer, you can compute the radius of gyration from that. Then try different beads types for the backbone to see which matches the RG data best.

 How can I know what to input in the atom? 

I don't really understand what you mean by this questions.

[kurokawaikki]

Thank you very much!! I now have some sense on how to find the appropriate bead type.
Do we have some tables listing all the building blocks as well as proteins in martini 3? I found the paper have a table of small molecules and building block, but I can not find the table with amino acid.

For example, I cannot find the SN2r in the building block table.
If there is no such info, maybe I can help to make a summarized table stating all beads type.

I now understand that SN2r -> PEG, [P2 and SC2] -> Leu

Because I am now in the mapping stage in the parametrization, I try to fully understand the meaning that you write.

Thank you very much for your great help! I can try to proceed to next step.
If have further question, I will ask you again here.
Thank you again for continuously helping to figure out how to proceed step by step.

After I finish this, if you needed, maybe I can help to write a small tutorial for beginners like me.

[fgrunewald]

@kurokawaikki any help with making tutorials for beginners is greatly appreciated. Hope the parametrization goes well!

[kurokawaikki]

@fgrunewald Thank you very much!

I encounter some problem on how to get the right parameter of my polymer.
I am wondering if you could provide me with some advice.

I followed the tutorial that you gave me. I choose to use Bartender xtb way to find the initial parameters, because I do not know how to guess the initial value on my own.

The following figures are the bond and angle distributions. Red is the CG model that Bartender gave me and the blue line is the atomistic simulation data parameter with ATB.
The CG beads model is shown as follow.

The bond distribution:

The angle distribution:

I have tried many times to change the value generated by Bartender to match the average bond length and angle that were calculated from atomistic simulation data.
However, I still cannot get a good match...
Could you give me some advice on this? Additionally, how can I obtain the force constant for bond length, angle and dihedrals? I can run quantum mechanic MD calculation with Amsterdam Modeling Suite or Maestro. But I do not know how to extract the needed information.
Thank you very much!

[kurokawaikki]

To be clearer, here is my CG itp file generated with Bartender.

`; Topology by Bartender - www.github.com/rmera/bartender
; Please cite the Bartender reference: XXXXXXXXXXX
[ moleculetype ]
; Name nrexcl
PGLE 1
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 SN3r 0 PGLE C1 1 0
2 SN3r 0 PGLE C2 2 0
3 SN3r 0 PGLE C3 3 0
4 SN3r 0 PGLE C4 4 0
5 SN3r 0 PGLE C5 5 0
6 SN2r 0 PGLE C6 6 0
7 P2 0 PGLE L61 7 0
8 TN6d 0 PGLE L62 8 1
9 C1 0 PGLE L63 9 0
10 SP2 0 PGLE C7 10 0
11 P2 0 PGLE L71 11 0
12 TN6d 0 PGLE L72 12 1
13 C1 0 PGLE L73 13 0

[bonds]
; i j funct length force.c.
1 2 1 0.35612 2072.45 ; rmsd: 1.30
2 3 1 0.33917 1965.26 ; rmsd: 1.58
3 4 1 0.34044 2010.50 ; rmsd: 1.40
4 5 1 0.33276 1672.18 ; rmsd: 1.50
5 6 1 0.336 1370.19 ; rmsd: 1.60
6 7 1 0.309 2139.18 ; rmsd: 1.99
7 8 1 0.346 15362.05 ; rmsd: 0.83
8 9 1 0.350 6271.65 ; rmsd: 2.46
6 10 1 0.263 3604.95 ; rmsd: 1.64
10 11 1 0.271 4276.60 ; rmsd: 1.76
11 12 1 0.35780 14041.80 ; rmsd: 1.01
12 13 1 0.357 9126.62 ; rmsd: 2.16
[constraints]
; i j funct length
[angles]
; i j k funct angle force_constant
1 2 3 1 130.30 16.94 ; rmsd: 1.03
2 3 4 1 128.54 16.93 ; rmsd: 1.05
3 4 5 1 121.132 18.55 ; rmsd: 0.94
4 5 6 1 123.17 16.89 ; rmsd: 0.65
5 6 7 1 70.6237 80.25 ; rmsd: 0.89
6 7 8 1 127.386 21.75 ; rmsd: 1.26
7 8 9 1 88.1504 50.27 ; rmsd: 1.86
5 6 10 1 126.06 16.37 ; rmsd: 1.68
7 6 10 1 129.135 12.52 ; rmsd: 1.64
6 10 11 1 84.1168 70.93 ; rmsd: 1.23
10 11 12 1 126.83 35.84 ; rmsd: 1.06
11 12 13 1 80.9223 66.02 ; rmsd: 2.11
;; 1 2 3 2 120.17 27.70 ; Harmonic-Cos (Gromos96) potential rmsd: 1.17
;; 2 3 4 2 118.67 27.25 ; Harmonic-Cos (Gromos96) potential rmsd: 1.42
;; 3 4 5 2 117.36 31.61 ; Harmonic-Cos (Gromos96) potential rmsd: 1.47
;; 4 5 6 2 113.46 26.83 ; Harmonic-Cos (Gromos96) potential rmsd: 1.19
;; 5 6 7 2 81.27 99.62 ; Harmonic-Cos (Gromos96) potential rmsd: 1.07
;; 6 7 8 2 108.24 39.68 ; Harmonic-Cos (Gromos96) potential rmsd: 1.61
;; 7 8 9 2 86.61 70.44 ; Harmonic-Cos (Gromos96) potential rmsd: 1.98
;; 5 6 10 2 116.08 27.71 ; Harmonic-Cos (Gromos96) potential rmsd: 1.83
;; 7 6 10 2 113.12 22.20 ; Harmonic-Cos (Gromos96) potential rmsd: 1.83
;; 6 10 11 2 94.19 97.71 ; Harmonic-Cos (Gromos96) potential rmsd: 1.34
;; 10 11 12 2 118.17 60.73 ; Harmonic-Cos (Gromos96) potential rmsd: 2.06
;; 11 12 13 2 87.19 82.99 ; Harmonic-Cos (Gromos96) potential rmsd: 2.16
[angles]
; i j k funct angle force_constant
; ReB
[dihedrals]
; i j k l funct phase kd pn
;Improper
; i j k l funct angle kd `

I modified the Leu part. Originally, Leu was represented with 2 beads. I changed to 3 beads and take on of the beads contains amine to be positively charged.
Because the positively charged amine and the negatively charged PS lipid are very important in my simulation, I want to represent them closely with CG model.

Could you give me some advice on how to derive correct parameters especially on the force constant?
How to calculate the force constant?

Thank you very much.

[fgrunewald]

@kurokawaikki I have too look at the data and need to ask some questions about bartender to some experts in our group. I'm not very familiar with the software. Then I get back to you. Overall it, however, doesn't look bad. The discrepancy in bond length can be related to the fact that bartender doesn't sample the atomistic dihedral conformations.

[fgrunewald]

@kurokawaikki let's start with a small checklist:

• Did you include the hydrogen atoms in your mapping (i.e. center-of-geometry)?
• Did you map the PEG fragment as CH2-O-CH2?
• Did bartender give any log-output, warning etc.?

Secondly regarding the LEU, you have to be very careful remapping these residues. Typically it is better to change the bead-type when there is a charge present. Like the C/N terminals of proteins have the same number of beads, but just a different bead type.

[kurokawaikki]

Did you include the hydrogen atoms in your mapping (i.e. center-of-geometry)?

-> Yes I have selected the H atom when I building the model in the CG builder.

Did you map the PEG fragment as CH2-O-CH2?

-> Yes, I mapped it as CH2-O-CH2. But the last one at the tail is CH2-OH

Did bartender give any log-output, warning etc.?

-> There is no log output and no warning. But I noticed that some of the xtb data and fitted function (generated by Bartender) did not matched in some of the output figures. For example, the first one is good fitted function and the second one is not.

2. Most of the poor fitted functions are in the cos angle output figure.
   
   

Because I did not add any information on dihedrals (my polymer contains single bond in back bone so I think it is not necessary), will this be the problem?

Secondly regarding the LEU, you have to be very careful remapping these residues. Typically it is better to change the bead-type when there is a charge present. Like the C/N terminals of proteins have the same number of beads, but just a different bead type.

I change the Leu to three beads... So I have to change back to two beads and change the bead type of beads containing C/N (positively charged). Which bead type do you suggest?

Additionally, could you tell me in general way how do you derive the force constant without Bartender? In the tutorial, it only stated that you should try and error to fit the parameter to atomistic simulation results. I understand that I can fit the bond angel and bond length and even dihedrals, but I cannot find any tutorial on how to derive force constant. Could you also give me some advice on that?

Thank you very much!

[fgrunewald]

@kurokawaikki I understand. My opinion is that for a flexible fragment like polymer bartender's sampling of the degrees of freedom is not good enough. That's why you see a difference in bond length and angles. On top of that come the strange fits that you show. Probably using a classical approach of running an all-atom simulation, mapping it and than matching the distributions is better. At least for the PEG part you can take the parameters from the library.

Having said that I think your biggest problem is to guess force-constants. The way it is usually done is the following:

• run an all-atom simulation of the fragment using any force-field you like
• map the trajectory using gmx trjconv or see below for some other tools
• compute the angle, and bond distributions using gromacs tools like 'gmx angle' and 'gmx distance' from the mapped trajectory
• set the bond length and average angle to the average value of the distributions
• guess a force-constant for bonds and angles
• iteratively run a CG simulation, analyze the angles and bonds again and adjust the force-constants

As a first guess: Bond force constants typically are between 2000-7000 anything about 9000 should probably be a constraint. Angles typically range from 10-100 in terms of force-constant.

There are other tools like pycgmap or fast_forward that can do the mapping and extract distributions. In principle SwarmCG can automatically derive bond and angle force-constants using a machine learning approach. However, SwarmCG can be resource intensive and is perhaps not the most easy to use. But it could work in your case.

[fgrunewald]

@ricalessandri this issue I mentioned

[kurokawaikki]

Dear Dr. @fgrunewald,
I followed your guide using SwarmCG. I successfully get the optimized parameter from the SwarmCG. It needs some modification in the code of SwarmCG for letting me running on my purpose. The following is the figure showing the fitting.

I have some question in the following:

1. I am wondering that normally how we can define good enough fitting of CG to AA model.

2. Additionally, because you mentioned that it is better to keep the same number of beads. In your experience, how will we choose the beads for positively charged C/N groups? In my feeling that I will try to change SC2 directly to TN6d (but maybe not tiny beads will work better?). I will build another CG model on that and compare the goodness of fitting.

3. Most important of all, how can I have a beautiful figures or movies with my simulated trajectory or CG molecules? What program are you usually using for making publishable figures on papers? In my experience with Desmond, I can use Maestro, PyMol or VMD for drawing. I have tried to import my CG model.gro into VMD, but it only showed up as small dots without linkage between beads.

I think after I got the right parameter, I can use it to build my polymer with polyply.
Then, build the membrane system with insane.py. If I encounter some problem along the road, can I ask you question here also?

Finally, as for the tutorial for the beginners on Martini CG, I can write my experience on this topic and explain along the steps.

Thank you very much again.

[fgrunewald]

Hi @kurokawaikki,

I think the SwarmCG distributions look good. Very good job on that! To answer your questions:

1. You try to achieve the best fit possible by visually inspecting. Sometimes that is impossible, especially when the distributions are bimodal. For example, I would perhaps want to optimize angle 7 and 12 a bit more. When the distributions are bimodal it is often better to make the angle or bond weaker such that the CG distribution covers both beaks to at least FWHM or completely. Simply set the reference angle to the average and reduce the force-constant. Other than that there is no clear cut definition. Just make sure your model runs stable with 20fs time-step.

2. For changing from neutral to charged beads you typically keep the size, which is dictated by the number of heavy atoms. The bead type follows from the chemical nature of the fragment. You can have a look at the protein terminals for inspiration, which are Q5 in both cases. Otherwise the SI of the paper contains a more detailed list. If you have all parameters collected I can also have a final look to see, if there is anything that could be improved. Note that tiny beads are no per se more precise just because the mapping resolution is higher. Choosing the bead that best matches the volume of the underlying fragment is the name of the game.

3. We use VMD but you need special input to visualize the bonds. Either you use gmx editconf to write connect records or you use the cg_bonds script.

When you're happy with the parameters you can use polyply to make the itp files and coordinates for any length of polymer. In case you have any questions on how to write the .ff input files let me know in this issue. If you want we can add the parameters to the library down the road as well. I also would be super happy about a tutorial. Perhaps you can open a discussion thread about the tutorial here. Then we keep it separated from this particular case.

You can use insane or TS2CG to make the membrane. I would use TS2CG in case you want to make a vesicle. If you have questions about making membranes/vesicles I'm happy to help. However, I would ask you to raise a new issue here and tag me in it. Then either I or someone else will help you.

Let me know if it all works out in the end either way. I'm curious now.

Cheers,
Fabian

[kurokawaikki]

Dear Dr. @fgrunewald ,

I have done a parameterization on both model. I increase the model a little bit to include PG > PG >> PG parameters.
In the first one, I use the same number of beads as in the Martini 3 library for Leu. Only change the beads type of protonated amine part to Q5.
In the second one, I use the same model that I used in the previous one.

The parameter of first one is:

[ moleculetype ]
; molname        nrexcl
PGLE             1
[ atoms ]
; id type resnr residue   atom  cgnr    charge     mass

1    SN3r     0    PGLE     C1  1      0.00000     54.00
2    SN3r     0    PGLE     C2  2      0.00000     54.00
3    SN3r     0    PGLE     C3  3      0.00000     54.00
4    SN3r     0    PGLE     C4  4      0.00000     54.00
5    SN3r     0    PGLE     C5  5      0.00000     54.00
6    SN2r     0    PGLE     P1  6      0.00000     54.00
7      P2     0    PGLE    L11  7      0.00000     72.00
8      Q5     0    PGLE    L12  8      1.00000     72.00
9    SN2r     0    PGLE     P2  9      0.00000     54.00
10     P2     0    PGLE    L21  10     0.00000     72.00
11     Q5     0    PGLE    L22  11     1.00000     72.00
12   SN2r     0    PGLE     P3  12     0.00000     54.00
13     P2     0    PGLE    L31  13     0.00000     72.00
14     Q5     0    PGLE    L32  14     1.00000     72.00
15    SP2     0    PGLE     P4  15     0.00000     54.00
16     P2     0    PGLE    L41  16     0.00000     72.00
17     Q5     0    PGLE    L42  17     1.00000     72.00


[ bonds ]
;   i     j   funct   length   force.c.

; bond type 1
    1     2       1    0.336  5555.89           ; 1

; bond type 2
    2     3       1    0.337  6948.63           ; 2

; bond type 3
    3     4       1    0.336  3789.17           ; 3

; bond type 4
    4     5       1    0.334  3229.54           ; 4

; bond type 5
    5     6       1    0.312  2561.31           ; 5

; bond type 6
    6     7       1    0.349  2872.24           ; 6

; bond type 7
    7     8       1    0.354  7457.72           ; 7

; bond type 8
    6     9       1    0.332  2007.60           ; 8

; bond type 9
    9    10       1    0.350  2859.17           ; 9

; bond type 10
   10    11       1    0.375  2183.87           ; 10

; bond type 11
    9    12       1    0.341  2472.32           ; 11

; bond type 12
   12    13       1    0.366  1991.99           ; 12

; bond type 13
   13    14       1    0.353  2918.89           ; 13

; bond type 14
   12    15       1    0.299  10436.35           ; 14

; bond type 15
   15    16       1    0.328  4510.66           ; 15

; bond type 16
   16    17       1    0.387  1555.59           ; 16


[ angles ]
;   i     j     k   funct     angle   force.c.

; angle type 1
    1     2     3       1    150.07      7.75           ; 1

; angle type 2
    2     3     4       1    137.68     15.26           ; 2

; angle type 3
    3     4     5       1    132.61      7.49           ; 3

; angle type 4
    4     5     6       1    142.88      7.53           ; 4

; angle type 5
    5     6     7       1     67.89    108.03           ; 5

; angle type 6
    6     7     8       1    135.33     19.30           ; 6

; angle type 7
    7     6     9       1    147.98     12.66           ; 7

; angle type 8
    5     6     9       1    105.53      4.98           ; 8

; angle type 9
    6     9    10       1     73.67     75.94           ; 9

; angle type 10
    9    10    11       1    113.67     17.14           ; 10

; angle type 11
    6     9    12       1    132.94     11.50           ; 11

; angle type 12
   10     9    12       1    137.23     17.61           ; 12

; angle type 13
    9    12    15       1    166.40      5.16           ; 13

; angle type 14
    9    12    13       1     64.14     48.89           ; 14

; angle type 15
   12    13    14       1    115.71     12.70           ; 15

; angle type 16
   13    12    15       1    135.23     15.65           ; 16

; angle type 17
   12    15    16       1     78.29    118.48           ; 17

; angle type 18
   15    16    17       1    129.52     31.56           ; 18

The goodness of fitting is shown in the figure:

I have tried different combination of force constant, length, angle, but I can not get further improvement.

The parameter for second model is:

[ moleculetype ]
; molname        nrexcl
PGLE             1

[ atoms ]
; id type resnr residue   atom  cgnr    charge     mass

1    SN3r     0    PGLE     C1  1      0.00000     54.00
2    SN3r     0    PGLE     C2  2      0.00000     54.00
3    SN3r     0    PGLE     C3  3      0.00000     54.00
4    SN3r     0    PGLE     C4  4      0.00000     54.00
5    SN3r     0    PGLE     C5  5      0.00000     54.00
6    SN2r     0    PGLE     P1  6      0.00000     54.00
7      P2     0    PGLE    L11  7      0.00000     72.00
8    TN6d     0    PGLE    L12  8      1.00000     36.00
9      C1     0    PGLE    L13  9      0.00000     72.00
10   SN2r     0    PGLE     P2  10     0.00000     54.00
11     P2     0    PGLE    L21  11     0.00000     72.00
12   TN6d     0    PGLE    L22  12     1.00000     36.00
13     C1     0    PGLE    L23  13     0.00000     72.00
14   SN2r     0    PGLE     P3  14     0.00000     54.00
15     P2     0    PGLE    L31  15     0.00000     72.00
16   TN6d     0    PGLE    L32  16     1.00000     36.00
17     C1     0    PGLE    L33  17     0.00000     72.00
18    SP2     0    PGLE     P4  18     0.00000     54.00
19     P2     0    PGLE    L41  19     0.00000     72.00
20   TN6d     0    PGLE    L42  20     1.00000     36.00
21     C1     0    PGLE    L43  21     0.00000     72.00


[ bonds ]
;   i     j   funct   length   force.c.

; bond type 1
    1     2       1    0.340  5561.50           ; 1

; bond type 2
    2     3       1    0.336  8672.95           ; 2

; bond type 3
    3     4       1    0.333  3275.97           ; 3

; bond type 4
    4     5       1    0.341  2968.34           ; 4

; bond type 5
    5     6       1    0.328  3466.86           ; 5

; bond type 6
    6     7       1    0.360  2015.48           ; 6

; bond type 7
    7     8       1    0.292  11437.23           ; 7

; bond type 8
    8     9       1    0.311  5334.60           ; 8

; bond type 9
    6    10       1    0.317  2742.85           ; 9

; bond type 10
   10    11       1    0.354  2315.34           ; 10

; bond type 11
   11    12       1    0.274  9070.97           ; 11

; bond type 12
   12    13       1    0.330  9758.50           ; 12

; bond type 13
   10    14       1    0.322  2290.21           ; 13

; bond type 14
   14    15       1    0.356  2021.23           ; 14

; bond type 15
   15    16       1    0.272  10525.03           ; 15

; bond type 16
   16    17       1    0.325  15393.71           ; 16

; bond type 17
   14    18       1    0.290  9106.59           ; 17

; bond type 18
   18    19       1    0.309  2533.98           ; 18

; bond type 19
   19    20       1    0.282  10184.45           ; 19

; bond type 20
   20    21       1    0.321  10114.37           ; 20


[ angles ]
;   i     j     k   funct     angle   force.c.

; angle type 1
    1     2     3       1    140.43     11.25           ; 1

; angle type 2
    2     3     4       1    158.80     14.05           ; 2

; angle type 3
    3     4     5       1    151.24      8.22           ; 3

; angle type 4
    4     5     6       1    145.36      8.79           ; 4

; angle type 5
    5     6     7       1     69.81     61.77           ; 5

; angle type 6
    6     7     8       1    141.51     16.10           ; 6

; angle type 7
    7     8     9       1     85.63    110.40           ; 7

; angle type 8
    5     6    10       1    123.14      4.01           ; 8

; angle type 9
    7     6    10       1    146.27      9.10           ; 9

; angle type 10
    6    10    11       1     76.99     47.56           ; 10

; angle type 11
   10    11    12       1    136.39     20.75           ; 11

; angle type 12
   11    12    13       1     76.11    533.06           ; 12

; angle type 13
    6    10    14       1    120.36      5.06           ; 13

; angle type 14
   11    10    14       1    135.30     19.47           ; 14

; angle type 15
   10    14    15       1     73.56    178.97           ; 15

; angle type 16
   14    15    16       1    134.62     38.92           ; 16

; angle type 17
   15    16    17       1     68.11    276.67           ; 17

; angle type 18
   10    11    12       1    141.18     26.66           ; 18

; angle type 19
   10    14    18       1    129.47     10.72           ; 19

; angle type 20
   15    14    18       1    134.98     16.27           ; 20

; angle type 21
   14    18    19       1     80.52     66.12           ; 21

; angle type 22
   18    19    20       1    142.33     18.65           ; 22

; angle type 23
   19    20    21       1     65.54    271.53           ; 23

The goodness of fitting is shown here:

I think the first one is a little bit better.

Additionally, do we have the list for all bead types corresponding to chemical structures? Am I missing something in the papers, because I did not find the whole list in martini 3 paper and SI?

Finally, here is the ff file that I wrote. It seems working...
I have change the (3) bead type to Q5.
And I change the angle for PG >PG >SC1 to PG >PG >BB, because I only fit this angle for now...
Will this affect the results?

[ moleculetype ]
; molname       nrexcl

PGLeu                1

[ atoms ]

;id type resnr residu atom cgnr   charge

 1   SN2r 1     PGLeu    PG      1       0    ; perhaps this bead type needs to be changed

 2   P2   1     PGLeu    BB       2      0    ; perhaps this bead type needs to be changed

 3   Q5  1     PGLeu    SC1     3      0


[ bonds ]

 2    3    1       0.354     7457.72

[ angles ]

1    2     3    1   135.33 19.30


;;;;; linking the PGLeu together

[ link ]

resname "PGLeu"

[ bonds ]

PG  >PG   1 0.341 2472.32 ; this introduces a bond between to bonded PELeu monomers connecting the PG beads



[ link ]

resname "PGLeu"

[ angles ]

PG  >PG  >>PG  1  132.94 11.50 {"group": "BB angles"}



[ link ]

resname "PGLeu"

[ angles ]

PG  >PG  >BB  1  64.14 48.89 {"group": "BB-SC angles"} ; angle between side-chain and backbone



;;;;; link between PEG and PGLEU



[ link ]

resname "PGLeu|PEO"

[ bonds ]

EC  >PG   1  0.312 2561.31 {"comment": "PEG-PGLeu bond"}

[ angles ]

EC   EC  >PG  1 142.88 7.53 {"group": "PEG-PGLeu angles"}

EC >PG >>PG 1 105.53 4.98 {"group": "PEG-PGLeu angles"}

EC >PG >BB 1 67.89 108.03 {"group": "PEG-PGLeu angles"}

I read the documentation for FF. However, the syntax is different from the sample that you provided. Is there some place I can find more documentation?

Thank you very much for your great help.
I will start up a discussion there to write a tutorial on how to get parameterization in a simple way and use it with Polyply.

The final initial.itp is:

; /home/ikki/MD/polyply-env/bin/polyply gen_params -lib martini3 -f pegleu.ff -seq PEO:10 PGLeu:10 -o init.itp -name PEO-b-PGLeu

; Please cite the following papers:
; Souza, P C T; Alessandri, R; Barnoud, J; Thallmair, S; Faustino, I; Grünewald, F; Patmanidis, I; Abdizadeh, H; Bruininks, B M H; Wassenaar, T A; Kroon, P C; Melcr, J; Nieto, V; Corradi, V; Khan, H M; Domański, J; Javanainen, M; Martinez-Seara, H; Reuter, N; Best, R B; Vattulainen, I; Monticelli, L; Periole, X; Tieleman, D P; de Vries, A H; Marrink, S J;  Nature Methods 2021; 10.1038/s41592-021-01098-3
; Grunewald, F; Alessandri, R; Kroon, P C; Monticelli, L; Souza, P C; Marrink, S J;  Nature Communications 2022; 10.1038/s41467-021-27627-4

[ moleculetype ]
PGLE 1

[ atoms ]
 1 SN3r  1 PEO   EC   1 0.0 45.0
 2 SN3r  2 PEO   EC   2 0.0 45.0
 3 SN3r  3 PEO   EC   3 0.0 45.0
 4 SN3r  4 PEO   EC   4 0.0 45.0
 5 SN3r  5 PEO   EC   5 0.0 45.0
 6 SN3r  6 PEO   EC   6 0.0 45.0
 7 SN3r  7 PEO   EC   7 0.0 45.0
 8 SN3r  8 PEO   EC   8 0.0 45.0
 9 SN3r  9 PEO   EC   9 0.0 45.0
10 SN3r 10 PEO   EC  10 0.0 45.0
11 SN2r 11 PGLeu PG  11 0.0     
12 P2   11 PGLeu BB  12 0.0     
13 Q5   11 PGLeu SC1 13 1.0     
14 SN2r 12 PGLeu PG  14 0.0     
15 P2   12 PGLeu BB  15 0.0     
16 Q5   12 PGLeu SC1 16 1.0     
17 SN2r 13 PGLeu PG  17 0.0     
18 P2   13 PGLeu BB  18 0.0     
19 Q5   13 PGLeu SC1 19 1.0     
20 SN2r 14 PGLeu PG  20 0.0     
21 P2   14 PGLeu BB  21 0.0     
22 Q5   14 PGLeu SC1 22 1.0     
23 SN2r 15 PGLeu PG  23 0.0     
24 P2   15 PGLeu BB  24 0.0     
25 Q5   15 PGLeu SC1 25 1.0     
26 SN2r 16 PGLeu PG  26 0.0     
27 P2   16 PGLeu BB  27 0.0     
28 Q5   16 PGLeu SC1 28 1.0     
29 SN2r 17 PGLeu PG  29 0.0     
30 P2   17 PGLeu BB  30 0.0     
31 Q5   17 PGLeu SC1 31 1.0     
32 SN2r 18 PGLeu PG  32 0.0     
33 P2   18 PGLeu BB  33 0.0     
34 Q5   18 PGLeu SC1 34 1.0     
35 SN2r 19 PGLeu PG  35 0.0     
36 P2   19 PGLeu BB  36 0.0     
37 Q5   19 PGLeu SC1 37 1.0     
38 SN2r 20 PGLeu PG  38 0.0     
39 P2   20 PGLeu BB  39 0.0     
40 Q5   20 PGLeu SC1 40 1.0     

[ bonds ]
11 12 1 0.350 2859.17
12 13 1 0.354 7457.72
14 15 1 0.350 2859.17
15 16 1 0.354 7457.72
17 18 1 0.350 2859.17
18 19 1 0.354 7457.72
20 21 1 0.350 2859.17
21 22 1 0.354 7457.72
23 24 1 0.350 2859.17
24 25 1 0.354 7457.72
26 27 1 0.350 2859.17
27 28 1 0.354 7457.72
29 30 1 0.350 2859.17
30 31 1 0.354 7457.72
32 33 1 0.350 2859.17
33 34 1 0.354 7457.72
35 36 1 0.350 2859.17
36 37 1 0.354 7457.72
38 39 1 0.350 2859.17
39 40 1 0.354 7457.72
 1  2 1 0.36 7000 ; central-bond
 2  3 1 0.36 7000 ; central-bond
 3  4 1 0.36 7000 ; central-bond
 4  5 1 0.36 7000 ; central-bond
 5  6 1 0.36 7000 ; central-bond
 6  7 1 0.36 7000 ; central-bond
 7  8 1 0.36 7000 ; central-bond
 8  9 1 0.36 7000 ; central-bond
 9 10 1 0.36 7000 ; central-bond
11 14 1 0.341 2472.32
14 17 1 0.341 2472.32
17 20 1 0.341 2472.32
20 23 1 0.341 2472.32
23 26 1 0.341 2472.32
26 29 1 0.341 2472.32
29 32 1 0.341 2472.32
32 35 1 0.341 2472.32
35 38 1 0.341 2472.32
10 11 1 0.312 2561.31 ; PEG-PGLeu bond

[ angles ]
11 12 13 1 135.33 19.30
14 15 16 1 135.33 19.30
17 18 19 1 135.33 19.30
20 21 22 1 135.33 19.30
23 24 25 1 135.33 19.30
26 27 28 1 135.33 19.30
29 30 31 1 135.33 19.30
32 33 34 1 135.33 19.30
35 36 37 1 135.33 19.30
38 39 40 1 135.33 19.30
 1  2  3 2 123.00 80 ; central-angle
 2  3  4 2 123.00 80 ; central-angle
 3  4  5 2 123.00 80 ; central-angle
 4  5  6 2 123.00 80 ; central-angle
 5  6  7 2 123.00 80 ; central-angle
 6  7  8 2 123.00 80 ; central-angle
 7  8  9 2 123.00 80 ; central-angle
 8  9 10 2 123.00 80 ; central-angle

; BB angles
11 14 17 1 132.94 11.50
14 17 20 1 132.94 11.50
17 20 23 1 132.94 11.50
20 23 26 1 132.94 11.50
23 26 29 1 132.94 11.50
26 29 32 1 132.94 11.50
29 32 35 1 132.94 11.50
32 35 38 1 132.94 11.50

; BB-SC angles
11 14 15 1 64.14 48.89
14 17 18 1 64.14 48.89
17 20 21 1 64.14 48.89
20 23 24 1 64.14 48.89
23 26 27 1 64.14 48.89
26 29 30 1 64.14 48.89
29 32 33 1 64.14 48.89
32 35 36 1 64.14 48.89
35 38 39 1 64.14 48.89

; PEG-PGLeu angles
 9 10 11 1 142.88 7.53
10 11 14 1 105.53 4.98
10 11 12 1 67.89 108.03

[ dihedrals ]
#ifndef TI
 1  2  3  4 11 8 0.6570 -1.3278 -0.43661278 1.0808704 0.680055
 2  3  4  5 11 8 0.6570 -1.3278 -0.43661278 1.0808704 0.680055
 3  4  5  6 11 8 0.6570 -1.3278 -0.43661278 1.0808704 0.680055
 4  5  6  7 11 8 0.6570 -1.3278 -0.43661278 1.0808704 0.680055
 5  6  7  8 11 8 0.6570 -1.3278 -0.43661278 1.0808704 0.680055
 6  7  8  9 11 8 0.6570 -1.3278 -0.43661278 1.0808704 0.680055
 7  8  9 10 11 8 0.6570 -1.3278 -0.43661278 1.0808704 0.680055
#endif

[fgrunewald]

Hi @kurokawaikki,

Your SwarmCG results look good to me for either mapping and I agree the first one is a little better.

Additionally, do we have the list for all bead types corresponding to chemical structures?
Am I missing something in the papers, because I did not find the whole list in martini 3 paper and SI?

We have some more detailed tables internally related to currently unpublished molecules mostly. It is anyways best to verify bead-type assignments against some form of experimental data.

And I change the angle for PG >PG >SC1 to PG >PG >BB, because I only fit this angle for now...
Will this affect the results?

I think only specifying the PG PG BB angle is fine.

I read the documentation for FF. However, the syntax is different from the sample that you provided. 
Is there some place I can find more documentation?

The FF syntax has a very large range with many specialized short-cuts. The documentation provided at the moment gives a conservative description that will always work. However, with the upcoming documentation of martinize2/vermouth more information on the FF format will become available. If you have any specific questions feel free to ask here. But I take your model is working for now.

A last remark perhaps. Be aware that the Membrane Affinity of the current PEG model is overestimated. We have developed a preliminary update to the Martini interaction table that is currently validated. For more details on that you can write me an email. Unfortunately I cannot disclose this publicly yet.

Cheers,
Fabian

[kurokawaikki]

Dear Dr. @fgrunewald,

Thank you very much. I thought my model is working. So, I will go to the next step on making a membrane system and the polymeric micelles to test the endosomal escape properties. If I have problem on that, can I ask you on the different github issues?

In addition, I will start to writing the tutorial.

Thank you very much again! Hopefully, it will work at the end~

[fgrunewald]

@kurokawaikki Happy to help with this! I hope it works for you. If you have any problems feel free to raise another issue. I'm closing this one for now. Also thanks for putting in the tutorial draft. If you want to expand it, I think the outline is quite complete already.