search

marrink-lab / polyply_1.0

Generate input parameters and coordinates for atomistic and coarse-grained simulations of polymers, ssDNA, and carbohydrates
Apache License 2.0
119 stars 21 forks source link

all atomic block copolymer building #355

Closed WeixZeng closed 2 months ago

WeixZeng commented 6 months ago

Thank you for creating such a wonderful toolkit. Trying to build an all atomic block copolymer model for GROMACS. Thought it should be available with the '-connect' function with the gen_seq, but don't understand the syntax:

,-, Would be great if you could provide some guildlines for this, one command line, let's say, for block copolymer version of PMMA and PMA like PMMA_b_PMA should be sufficient for my case.
fgrunewald commented 6 months ago

Dear @WeixZeng,

Let's say you have 100 residues each the syntax would be as follows:

polyply gen_seq -from_string A:100:1:PMA-1.0 B:100:1:PMMA-1.0 -seq A B -connect 0:1:99-0 -o PMMA_b_PMA.json -label 0:chiral:R-0.5,S-0.5 1:chiral:R-0.5,S-0.5 -name PMMA_b_PMA

The connect is interpreted as follows: Take seq block 0 ( i.e. A) and connect it to block B (i.e. 1) by putting a bond between residues 0 and 99. All counts are zero-indexed (i.e. first residue has resid 9 and 100 residues have a resid of 0 to 99).

Now you can generate a GROMOS itp file using:

polyply gen_params -lib 2016H66 -seqf PMMA_b_PMA.json -o PMMA_b_PMA.itp -name PMMA_b_PMA

We also have set up (but not released) a workflow to generate itp files for other force fields like OPLS or CHARMM36 for which monomers are not part of the library yet. Depending on your problem it might be more appropriate to use these force fields. Let me know if you are interested and I will show you how to set up those polymers.

WeixZeng commented 6 months ago

Hi @fgrunewald ,

I tried workflow you taught. It works well and can generate .gro file for the target polymer melt, but I also found some problems with the coordinates.

  1. Incorrect structure for the monomer. I tried PS and PMMA, respectively. Found only 13 atoms in one monomer, the three H atoms in the polyvinyl backbone are omitted, while all H atoms in PMMA are omitted. I've double checked the polymer library, found the mistakes should be sourced from the [2016H66]/polyvinyl_blocks.ff forcefield. polyply/data/2016H66/polyvinyl_blocks.ff Do you think it necessary to update the structure in this forcefield, and how can I help?

  2. Distance between monomers are way too long in the generated .gro files. The monomers seperated like solvant in the space, doesn't look like a polymer chain. Guess it still has something to do with the GROMOS 2016H66 forcefield. Would be great if I can try the workflow to generate other itp files like OPLS or OPLS-AA.

Many thanks!

fgrunewald commented 6 months ago

Hi @WeixZeng,

  1. There is no problem with the monomers; GROMOS is a united atom force-field where carbons with hydrogen (except aromatic ones) are represented by a single particle (e.g. CH2 -> 1 particle).

  2. You need to run an energy minimization using -Deg_polyply in the mpd file (i.e. add define = -Deq_polyply). The distances always look larger because they are rescaled in the last step of the backmapping. To obtain the final morphology you always need to run the energy minimization.

Can you try? If you decide to work with OPLS instead I'm sure @ricalessandri can help setting up the parameters.

WeixZeng commented 6 months ago

Hahh, thank you for your patience. I will try it right away.

Yes, I would like to work with OPLS maybe after working out my current issues and graduate from beginner period.

fgrunewald commented 6 months ago

@WeixZeng ok! I'll organize the parameters

ricalessandri commented 6 months ago

Yes, I can help with OPLS. Feel free to ping me here when it's time.

WeixZeng commented 6 months ago

@fgrunewald Thank you! The chains bacame much nicer after the EM process.

@ricalessandri Could you please share the workflow for genertating itp files with OPLS?

Many thanks!

WeixZeng commented 6 months ago

I checked a bit the tutorials, and found there's '1. oplsaaLigParGen' in the library list, but it only contains a few blocks at current version. If I understood the workflow correctly, to generte an itp file with OPLS-AA forcefield, it needs parameters of the blocks. Therefore, I should write the parameters first.

I actually also need to set up relevant coarse-grained polymer model and tried to follow provided tutorials, felt it quite good to simply switch a library with the gen_params function that generating itp files in different scale with the same .json sequence description.

I tried to generate the same block copolymer melt with Martini3 force field with following protocal:

  1. polyply gen_seq -from_string A:123:1:PS-1. B:129:1:PMMA-1. -seq A B -connect 0:1:122-0 -o PS_b_PMMA.json -name PS_b_PMMA

  2. polyply gen_params -lib martini3 -o PS_b_PMMA.itp -name PS_b_PMMA -seqf PS_b_PMMA_CG.json

with top file as:

include "martini_v3.0.0.itp"

include "PS_b_PMMA.itp"

[ system ] ; name Melt of PS-b-PMMA

[ molecules ] ; name number PS_b_PMMA 50 ~

  1. polyply gen_coords -p bcp_1.top -o melt_cg.gro -name melt_cg -dens 1050

But it run into an OSError says: "Molecule consistes of two disconnected parts. Make sure all atoms/particles in a molecule are connected by bonds, constraints or virual-sites"

Could you please help to figure out the problem here?

Many thanks!

ricalessandri commented 6 months ago

I'm sorry, I have been travelling but I should be able to reply to you today or max tomorrow.

WeixZeng commented 6 months ago

No worries and no hurry at all. Do take your time and enjoy your trip!

For the AA model, if oplsaaLigParGen is the right library, it would be great to share the workflow to write .ff file for a PMMA monomer with the chiral and the link site information as an example.

Cheers, Wei.

ricalessandri commented 5 months ago

Hi @WeixZeng, sorry, it's been busier than expected (and Fabian is on vacation).

[ Martini 3 ] The current polyply library doesn't have a link connecting PMMA and PS. These are the soonish-to-be-released advanced Martini 3 models for PMMA and PS, including their link. Namely, PMMA block and links:

;;;;;;;; PMMA
[ moleculetype ]
MMAC 1
;
[ atoms ]
;  nr    type       resnr  residu    atom    cgnr        charge
   1    SC3             1  MMAC       MVNL       1     0.0
   2    SN4a            1  MMAC       SC1        2     0.0
;
[ bonds ]
;  ai  aj   funct
   1     2   1     0.37 9000 {"comment": "MVNL_SC1"}

[ link ]
resname "MMAC"
[ bonds ]
MVNL     +MVNL    1     0.315 4000 {"comment": "MVNL_MVNL"}

[ link ]
resname "MMAC"
[ angles ]
MVNL  +MVNL  ++MVNL   10  140  40  {"group": "vinyl backbone", "comment": "MVNL_MVNL_MVNL"}

[ link ]
resname "MMAC"
[ angles ]
MVNL  SC1  +MVNL  2  48  800 {"comment": "MVNL_SC1_MVNL"}

PS block and links:

;;;;;;;; POLYSTYRENE
[ moleculetype ]
STYR 1

[ atoms ]
;  nr    type       resnr  residu    atom    cgnr        charge    mass
   1    TC3            1  STYR       VNL      1     0.0
   2    TC5            1  STYR       SC1      2     0.0
   3    TC5            1  STYR       SC2      3     0.0
   4    TC5            1  STYR       SC3      4     0.0

[ bonds ]
;  ai  aj   funct
   1     2   1     0.27 8000

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

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

[ angles ]
;  ai  aj  ak  funct
   1     2     3    10  136.000000  150.000000 {"comment": "BB_SC1_SC2", "ifndef": "TI", "version": "1"}
   1     2     3    1   136.000000  150.000000 {"comment": "BB_SC1_SC2", "ifdef": "TI", "version": "2"}
   1     2     4    1   82.000000  100.000000 {"comment": "BB_SC1_SC4"}

[ link ]
resname "STYR"
[ bonds ]
VNL     +VNL    1     0.275 8000  {"group": "vinly-backbone bonds", "comment": "VNL_VNL"}

[ link ]
resname "STYR"
[ dihedrals ]
VNL  +VNL  ++VNL +++VNL  11 1.30929464 6.51130136 -4.64191250 10.24767684 -7.64602475 12.34323750  {"group": "vinly-backbone dihedrals", "comment": "VNL_VNL_VNL_VNL"}

;;;;;; 1-4 pairs

[ link ]
resname "STYR"
[ atoms ]
SC1 { }
VNL { }
+VNL { }
+SC1 { }
[ pairs ]
SC1 +SC1 1 4.600000e-01 2.5e+00
[ edges ]
SC1 VNL
VNL +VNL
+VNL +SC1

And the PS-PMMA link:

[ link ]
resname "STYR|MMAC"
[ bonds ]
VNL +MVNL 1 0.42 5000 {"comment": "MVNL_VNL"}
[ angles ]
 VNL +MVNL ++MVNL  2  120 50 {"comment": "VNL_MVNL_MVNL"}

Just note that we use this nomenclature, so you need to use "MMAC" and "STYR" instead of "PMMA" and "PS":

| **MRU name** | **MRU molecule**       | **Homopolymer**         | **Homopolymer Abbr.** |
| ------------ | ---------------------- | ----------------------- | --------------------- |
| STYR         | Ethylbenzene           | Polystyrene             | PS                    |
| MMAC         | Methyl isobytyrate     | Polymethyl methacrylate | PMMA                  |

Let me know if something isn't clear.

I will soon share the OPLS parameters.

ricalessandri commented 5 months ago

@WeixZeng You can find OPLS parameters for PMMA: https://github.com/marrink-lab/polyply_1.0/pull/360/files Do you also need PMA? The workflow to generate .ff files with new parameters for OPLS is still under development so is a little bit early even for beta testing. However, I can do the generation for more polymers and share the parameters via PRs, if needed.

WeixZeng commented 5 months ago

Thank you @ricalessandri. Really helpful!

About the Martini 3. Looks like the main differences comparing with those in current library are the resname and the link connecting. Sorry, I'm still quite new at polyply either Martini, but how could I update the library (maybe at my local server) or use the new connecting parameters to generate input file with polyply?

About the OPLS-AA. I'm working on PMMA-PS system currently and might need to build some new monomer in the future, don't need PMA at this stage. I'll try a bit to figure out how to generate OPLS-AA .ff file for PS refering to the PMMA parameters you shared. Considering your development progress, it would be great if you could also share the OPLS model for PS.

ricalessandri commented 5 months ago

Regarding Martini 3, if you put all the lines I posted above from:

;;;;;;;; PMMA
[ moleculetype ]
MMAC 1
...

to:

...
[ angles ]
 VNL +MVNL ++MVNL  2  120 50 {"comment": "VNL_MVNL_MVNL"}

in one file, say PS_PMMA.ff, then you can just do:

polyply gen_params -f PS_PMMA.ff -o PS_b_PMMA.itp -name PS_b_PMMA -seqf PS_b_PMMA_CG.json

Note that I replaced the -lib martini3 with -f PS_PMMA.ff with respect to your command above (and assumed that the file "PS_PMMA.ff" is in the same folder from which you're executing the command).

Regarding OPLS, I'll share PS, too, then. And I will keep you posted on updates on the development.

ricalessandri commented 5 months ago

@WeixZeng you can find the PS .ff OPLS file at https://github.com/marrink-lab/polyply_1.0/pull/361/files

WeixZeng commented 5 months ago

@ricalessandri Then the question is to build PS/PMMA copolymers, thought the main thing needed is the PS-PMMA link.

I tried to write with the OPLS .ff files you shared, while the vinyl carbons changed to VC1 and VC2, i.e. in PS, C08→VC1 and C01→VC2; and in PMMA, C01→VC1, C02→VC2, the link parameter as follows:

[link]
resname "PS|PMMA"
[ atoms ]
[ bonds ]
 +VC2   VC1     1    0.1529    224262.400   {"comment":"intermonomer"}
[ angles ]
  VC2   VC1  +VC2     1    112.700   488.273   {"comment":"intermonomer"}
  VC2  +VC1  +VC2     1    112.700   488.273   {"comment":"intermonomer"}
  VC2  +VC1  +H01     1    110.700   313.800   {"comment":"intermonomer"}
  VC2  +VC1  +H02     1    110.700   313.800   {"comment":"intermonomer"}
 +VC2   VC1   H15     1    110.700   313.800   {"comment":"intermonomer"}
 +VC2   VC1   H16     1    110.700   313.800   {"comment":"intermonomer"}
[ dihedrals ]
[ dihedrals ]
 +C04  +VC2  +VC1   VC2     3     -4.960      6.286      1.310     -2.636     -0.000      0.000   {"comment":"intermonomer"}
 +VC2  +VC1   VC2   VC1     3      2.301     -1.464      0.837     -1.674     -0.000      0.000   {"comment":"intermonomer"}
 +C03  +VC2  +VC1   VC2     3      2.301     -1.464      0.837     -1.674     -0.000      0.000   {"comment":"intermonomer"}
  H02   VC1   VC2  +VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
 +H02  +VC1   VC2   VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
  H01   VC1   VC2  +VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
 +H01  +VC1   VC2   VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
  H16   VC1  +VC2  +VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
 +VC2   VC1   VC2   C02     3      2.301     -1.464      0.837     -1.674     -0.000      0.000   {"comment":"intermonomer"}
  H09   VC2   VC1  +VC2     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
  H15   VC1  +VC2  +VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
 +H09  +VC2   VC1   H1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
 +H09  +VC2   VC1   H2     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}  
[ pairs ]
  VC1   +VC2       1   {"comment":"intermonomer"}
  VC2   +C03       1   {"comment":"intermonomer"}
  VC2   +C04       1   {"comment":"intermonomer"}
  VC1   +H01       1   {"comment":"intermonomer"}
  VC1   +H02       1   {"comment":"intermonomer"}  
  C02   +VC2       1   {"comment":"intermonomer"}
 +VC1    H15       1   {"comment":"intermonomer"}
 +VC1    H16       1   {"comment":"intermonomer"}
 +VC2    H09       1   {"comment":"intermonomer"}  
  H01   +H09       1   {"comment":"intermonomer"}
  H02   +H09       1   {"comment":"intermonomer"}
[link]
; for bonded terms spanning three residues
resname "PS|PMMA|PMMA"
[dihedrals]
++VC1  +VC2  +VC1   VC2     3      2.301     -1.464      0.837     -1.674     -0.000      0.000   {"comment":"intermonomer"}
[pairs]
  VC2  ++VC1       1   {"comment":"intermonomer"}
[link]
resname "PMMA|PS"
[ atoms ]
[ bonds ]
 +VC1   VC2     1    0.1529    224262.400   {"comment":"intermonomer"}
[ angles ]
  VC1   VC2  +VC1     1    112.700   488.273   {"comment":"intermonomer"} 
  VC1  +VC2  +C02     1    114.000   527.184   {"comment":"intermonomer"}
  VC1  +VC2  +VC1     1    112.700   488.273   {"comment":"intermonomer"}  
  VC1  +VC2  +H09     1    110.700   313.800   {"comment":"intermonomer"}
  C03   VC2  +VC1     1    112.700   488.273   {"comment":"intermonomer"}
  C04   VC2  +VC1     1    111.100   527.184   {"comment":"intermonomer"}  
[ dihedrals ]
[ dihedrals ]
 +C07  +C02  +VC2   VC1     3      0.000      0.000      0.000     -0.000     -0.000      0.000   {"comment":"intermonomer"}
 +C03  +C02  +VC2   VC1     3      0.000      0.000      0.000     -0.000     -0.000      0.000   {"comment":"intermonomer"}
 +C02  +VC2   VC1   VC2     3      2.301     -1.464      0.837     -1.674     -0.000      0.000   {"comment":"intermonomer"}
 +VC1  +VC2   VC1   VC2     3      2.301     -1.464      0.837     -1.674     -0.000      0.000   {"comment":"intermonomer"}
 +H09  +VC2   VC1   VC2     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
 +H16  +VC1  +VC2   VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
 +H15  +VC1  +VC2   VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
  H05   C03   VC2  +VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
 +VC1   VC2   C04   O05     3      0.000      0.000      0.000     -0.000     -0.000      0.000   {"comment":"intermonomer"}
 +VC1   VC2   C04   O06     3     -1.157     -3.471      0.000      4.628     -0.000      0.000   {"comment":"intermonomer"}
 +VC2  +VC1   VC2   C04     3     -4.960      6.286      1.310     -2.636     -0.000      0.000   {"comment":"intermonomer"}
 +VC2  +VC1   VC2   C03     3      2.301     -1.464      0.837     -1.674     -0.000      0.000   {"comment":"intermonomer"}
  H04   C03   VC2  +VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
  H03   C03   VC2  +VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
 +H15  +VC1   VC2   C03     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
 +H16  +VC1   VC2   C03     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
 +H15  +VC1   VC2   C04     3     -0.209     -0.628      0.000      0.837     -0.000      0.000   {"comment":"intermonomer"}
 +H16  +VC1   VC2   C04     3     -0.209     -0.628      0.000      0.837     -0.000      0.000   {"comment":"intermonomer"}
  H01   VC1  +VC2  +C02     3      0.967      2.900      0.000     -3.866     -0.000      0.000   {"comment":"intermonomer"}
  H02   VC1  +VC2  +C02     3      0.967      2.900      0.000     -3.866     -0.000      0.000   {"comment":"intermonomer"}  
[ pairs ]
  VC2   +C02       1   {"comment":"intermonomer"}
  VC2   +VC1       1   {"comment":"intermonomer"}
  VC1   +C03       1   {"comment":"intermonomer"}
  VC1   +C07       1   {"comment":"intermonomer"}
  VC2   +H09       1   {"comment":"intermonomer"}
  VC1   +H15       1   {"comment":"intermonomer"}
  VC1   +H16       1   {"comment":"intermonomer"}  
 +VC1    H04       1   {"comment":"intermonomer"}
  C03   +H01       1   {"comment":"intermonomer"}
 +VC1    H05       1   {"comment":"intermonomer"}
  C04   +H01       1   {"comment":"intermonomer"}
  C03   +H02       1   {"comment":"intermonomer"}
  C04   +H02       1   {"comment":"intermonomer"}
  C03   +VC2       1   {"comment":"intermonomer"}
  O05   +VC1       1   {"comment":"intermonomer"}
  C04   +VC2       1   {"comment":"intermonomer"}
  O06   +VC1       1   {"comment":"intermonomer"}
 +C02    H01       1   {"comment":"intermonomer"}
 +C02    H02       1   {"comment":"intermonomer"}  
 +VC1    H03       1   {"comment":"intermonomer"}
 +VC1    H01       1   {"comment":"intermonomer"}
 +VC1    H02       1   {"comment":"intermonomer"}  
[link]
; for bonded terms spanning three residues
resname "PS|PMMA|PMMA"
[dihedrals]
++VC2  +VC1  +VC2   VC1     3      2.301     -1.464      0.837     -1.674     -0.000      0.000   {"comment":"intermonomer"}
[pairs]
  VC1  ++VC2       1   {"comment":"intermonomer"}

then tried use similar workflows as the UA or Martini3 forcefield, but doesn't work. The OPLS-AA forcefield used here was the OPLS-AA/M obtained from Jorgensen group's officail website:

https://traken.chem.yale.edu/doc/GMX_OPLSAAM.zip

Could you please help to check the PS-PMMA link and if the OPLS-AA forcefield is the right version, considering the PS and PMMA parameters you shared.

Many thanks.

fgrunewald commented 5 months ago

@WeixZeng could you share all the files including the PS PMMA base files you edited?

WeixZeng commented 5 months ago

Yes, @fgrunewald the oplsaam.ff and the PS PMMA file as attached. oplsaam.ff.zip

fgrunewald commented 5 months ago

@WeixZeng in order to find out why links do not apply you can run 

polyply gen_params -vv ...

I believe you have mixed up some atomnames in the link probably the pair section. When running polyply gen_params -seq PS:1 PMMA:1 polyply doesn't find H02 in PS, because there is none. But you have the following entires that require H02 to be present:

[ dihedrals ]
  H02   VC1   VC2  +VC1     3      0.628      1.883      0.000     -2.510     -0.000      0.000   {"comment":"intermonomer"}
[ pairs ]
  H02   +H09       1   {"comment":"intermonomer"}

Are these typos perhaps?

When running polyply gen_params -seq PMMA:1 PS:1 polyply doesn't find +H01 in PS, because there is none. But you have the following entires that require +H01 to be present in PS:

[ pairs ]
  C03   +H01       1   {"comment":"intermonomer"}

Also this one is perhaps a typo?

Other than that I think the file looks good.

fgrunewald commented 2 months ago

@WeixZeng can we close this issue or transfer it to the discussions section?

WeixZeng commented 2 months ago

yes for sure, thank you again for your help! ---- Replied Message ---- From Dr. Fabian @.> Date 05/27/2024 22:07 To @.> Cc @.>@.> Subject Re: [marrink-lab/polyply_1.0] all atomic block copolymer building (Issue #355) @WeixZeng can we close this issue or transfer it to the discussions section? — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: @.***>