How to construct PEO melts at atomistic level #320

YanChen32 · 2023-05-16T14:49:14Z

YanChen32
May 16, 2023

Dear all,

I am currently working on simulating a mixed melt of PEO and PVA using the OPLSAA force field. I have already completed all the tutorials provided by polyply_1.0. However, most of the examples in polyply_1.0 are currently focused on coarse-grained models. Therefore, I have a few questions regarding how to construct all-atom models using the OPLSAA force field.

Firstly, I would like to know the process of building all-atom models based on the OPLSAA force field in polyply_1.0. Is there a comprehensive example available that demonstrates the usage of OPLSAA in polyply_1.0?

Additionally, I would like to understand how to create a melt model from scratch using a customized OPLSAA force field file, e.g., PEO or PVA. Specifically, I am interested in the steps and procedures involved in this process.

Many thanks,

Yan

fgrunewald · 2023-05-17T09:42:51Z

fgrunewald
May 17, 2023
Maintainer

Dear Yan,

Thanks for opening this discussion! We'll try our best to guide you through the process. Unfortunately, there is no comprehensive example yet. However, we've been working with OPLS quite a bit (especially @ricalessandri is the expert here). Let's start:

Creating a polymer melt with OPLS FF

1 - generating an itp file

Creating an initial simulation setup with OPLS is fairly easy and involves two steps. First you need the itp file of the polymer and secondly you need initial system coordinates. If you want to do this for PEG for example you need the following command to make the itp file, which follows the same idea s for CG systems.

polyply gen_params -lib oplsaaLigParGen -seq OHter:1 PEO:100 OHter:1 -o PEO.itp -name PEO

This command will generate an itp file for 100 monomers of PEO using the OPLS LigParGen force-field. The LigParGen only means that the charges and atom-type assignment were done using LigParGen. We'll come to that in a minute. Note that you have to use the OHter:1 residue in order to terminate the PEO chain with an alcohol group. In your particular case we can later change those to CH3 if it's more suitable.

2 - building the initial melt

To generate a melt you need a topology file of the system. In the PEG case the non-bonded parameters are provided by the LigParGen server, so the following topology file will do. Note that we aim to generate 100 polymers of PEO.

#include "oplsaa.ff/forcefield.itp"
#include "PEO.itp"

[ system ]
test
[ molecules ]
PEO 100

Next you simply generate the melt at a particular density using the following command. This should not take very long, but depends a little bit on your python version and machine. For me it took 30s.

 polyply gen_coords -p system.top -o start.gro -dens 1000 -name melt

NOTE: @ricalessandri made an edit here given that, based on this recent PR, the required .top file is simpler.

3 - preparing to run a simulation

Part of the structure generation is an energy minization. When you look at the start.gro file you will see that the PEG looks like the bonds are not correct and the structure is bad. However, polyply scales the residue coordinates and they first have to "unflod" in an energy minimization. So next you run an energy minization using the default OPLS settings. Make sure however to do the energy minimization with flexible bonds not constraints. After that you can start the melt simulation in NVT/NpT as you want.

Creating parameters for a custom polymer

In order to create itp-files for a custom polymer, you need to define the bonded parameters, charges and atom-types for the repeat unit as well as how they change when they are linked together. All this information is stored in a .ff file. Luckily for you for OPLS @ricalessandri has created a small script which allows you to generate the .ff input files you need from an itp file of a tetramer for example. So typically we generate a tetramer using LigParGen server, then apply the script itp_to_ff.py, and finally you can generate the itp file as follows:

polyply gen_params -f *.ff   ...

Watch out for tacticity

Polyply does not define tacticity, instead the tacticity has to be defined via improper dihedral angles. You do this by adding the following entry to the ff-files. Let us suppose for PVA the main backbone atoms in the monomer are called CC1 CC2and you have theHCC2andO` as the names of the hydrogen connected to the backbone carbon, where also the O group is attached. In other words the PVA repeat unit looks as follows:

          HCC2
          |
 - CC1 - CC2 -
          |
          O - H

In this case you add the following two entries to any of the .ff files:

[ link ]
resname PVA
[ impropers ]
CC2 {"chiral": "S"} CC1 {"chiral": "S"} +CC1 {"chiral": "R|S"} O {"chiral": "S"} 2 -3.5264390e+01  3.3484625e+02 {"comment":"center S", "group":"link"}
[ edges ]
CC2  CC1
CC2 O
CC2 +CC1

[ link ]
resname "PVA"
[ impropers ]
CC2 {"chiral": "R"} CC1 {"chiral": "R"} +CC1 {"chiral": "S|R"} O {"chiral": "R"} 2 3.5264390e+01  3.3484625e+02 {"comment":"center R", "group":"link"}
[ edges ]
CC2  CC1
CC2 O
CC2 +CC1

Next to generate a polymer with random tacticity you have to define the chirality of each monomer randomly. To do so call polyply gen_seq as shown below.

polyply gen_seq -from_string A:100:1:PVA-1.0 -o PVA.json -name PVA -seq A -label 0:chiral:R-0.5,S-0.5

To generate iso-tactic PVA you use this command:

-label 0:chiral:R-1.0

Note that the last part -label 0:chiral:R-1.0 defines the percentage of monomers that have R or S chirality. Any float combination is valid. If you have more special patterns like triads that can also be included, but we can discuss when that is needed. The itp file now has to be generated using this sequence file:

polyply gen_params -f <your_ff_files> -seqf PVA.json -o PVA.itp -name PVA

That's it!

If you have any question or assistance let me know. In case you want to submit your PVA/PEG to the polyply library and make it available to others, we're happy to assist. And ultimately if you like the package cite us and give us a star on GitHub.

8 replies

YanChen32 May 23, 2023
Author

For example, I extract the ITP file for PVA as
`;
; GENERATED BY LigParGen Server
; Jorgensen Lab @ Yale University
;
[ atomtypes ]
opls_136 C136 12.0110 0.000 A 3.50000E-01 2.76144E-01
opls_158 C158 12.0110 0.000 A 3.50000E-01 2.76144E-01
opls_154 O154 15.9990 0.000 A 3.12000E-01 7.11280E-01
opls_140 H140 1.0080 0.000 A 2.50000E-01 1.25520E-01
opls_140 H140 1.0080 0.000 A 2.50000E-01 1.25520E-01
opls_140 H140 1.0080 0.000 A 2.50000E-01 1.25520E-01
opls_155 H155 1.0080 0.000 A 0.00000E+00 0.00000E+00
[ moleculetype ]
; Name nrexcl
PVA 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 opls_136 1 PVA C01 1 -0.2459 12.0110
2 opls_158 1 PVA C02 1 0.1759 12.0110
3 opls_154 1 PVA O03 1 -0.6888 15.9990
4 opls_140 1 PVA H04 1 0.1114 1.0080
5 opls_140 1 PVA H05 1 0.1114 1.0080
6 opls_140 1 PVA H06 1 0.1297 1.0080
7 opls_155 1 PVA H07 1 0.4103 1.0080
[ bonds ]
2 1 1 0.1529 224262.400
3 2 1 0.1410 267776.000
4 1 1 0.1090 284512.000
5 1 1 0.1090 284512.000
6 2 1 0.1090 284512.000
7 3 1 0.0945 462750.400

[ angles ]
; ai aj ak funct c0 c1 c2 c3
1 2 3 1 109.500 418.400
1 2 6 1 110.700 313.800
2 3 7 1 108.500 460.240
2 1 4 1 110.700 313.800
2 1 5 1 110.700 313.800
4 1 5 1 107.800 276.144
3 2 6 1 109.500 292.880

[ dihedrals ]
; IMPROPER DIHEDRAL ANGLES
; ai aj ak al funct c0 c1 c2 c3 c4 c5

[ dihedrals ]
; PROPER DIHEDRAL ANGLES
; ai aj ak al funct c0 c1 c2 c3 c4 c5
6 2 1 4 3 0.628 1.883 0.000 -2.510 -0.000 0.000
6 2 1 5 3 0.628 1.883 0.000 -2.510 -0.000 0.000
4 1 2 3 3 0.979 2.937 0.000 -3.916 -0.000 0.000
5 1 2 3 3 0.979 2.937 0.000 -3.916 -0.000 0.000
7 3 2 1 3 -0.444 3.833 0.728 -4.117 -0.000 0.000
7 3 2 6 3 0.736 2.209 0.000 -2.946 -0.000 0.000

[ pairs ]
3 4 1
3 5 1
4 6 1
5 6 1
6 7 1
`
Then, how to generate the corresponding .ff input files @ricalessandri .

fgrunewald May 23, 2023
Maintainer

@YanChen32 the idea I believe is that you can use the script, which you currently cannot access, and that would extract the ff-files from the itp file you generated. The issue is that apparently we have not yet made the script publicly available. I'm sure @ricalessandri can add you to it.

Otherwise you can also manually write the ff-files. You've already extracted the main repeat unit that is the block, so you can simply use that as one part of the input. Now you only need to write how links connect the blocks. We have written a tutorial here for GROMOS but it works the same way for OPLS. Below I've written two links that should cover everything but pairs. For the terminal groups you repeat the process.

; example link you need to fill in the values for the bonds and angles
[ link ]
resname "PVA"
[ bonds ]
C02  +C01  1
[ angles ]
C02 +C01 +C02  1
C01  C02 +C01  1
[ dihedrals ]
C01 C02 +C01 +C02 3 
; I think OPLS defines these dihedrals too but make sure to check
H04  C01 C02 +C01 3
H05  C01 C02 +C01 3
H06  C02 +C01 +C02 3
H06  C02 +C01 +H04 3
H06  C02 +C01 +H05 3
O03  C02 +C01 +C02 3
O03  C02 +C01 +H05 3
O03  C02 +C01 +H05 3
H07 O03 C02 +C01 3
; don't forget the pairs

[ link ]
resname "PVA"
[ dihedrals ]
C02 +C01 +C02 ++C01 3

ricalessandri May 24, 2023
Collaborator

Hi @YanChen32 , so indeed the small script called itp_to_ff.py is not public yet. @fgrunewald and I will work to publish an open-beta version soon, at least we need to put down some decent documentation for it. In the meantime, I uploaded the correct (i.e., validated against oligomers produced by LigParGen) PVA.ff I produced with itp_to_ff via this PR. Feel free to use it, just double-check that everything is in order.

Regarding the question about the difference in charges in termini, are you referring to -CH3 termini between PEO and PVA? If you see significant differences, you should go with different termini blocks, in my opinion. However, deciding how significant is significant, might be tricky.

Final note: as Fabian mentioned, PEO in the library currently is only available OH-terminated. In case you want it -CH3-terminated, you'll need to define the CH3 block and it link to the PEO block. It should be relatively simple but let me know if you need help with it. In any case it may be good for us to have it in the library.

YanChen32 May 24, 2023
Author

Thanks a lot for the help from @ricalessandri and @fgrunewald. I have already attempted manually write the ff-files, but it's a laborious process that is susceptible to errors and omissions. Hence, I am eagerly anticipating the open-source release of the script "itp_to_ff.py" as it will greatly expedite the implementation of OPLS-AA force fields in Polyply.

Regarding the PVA system, could you assist me in generating the ff-file with -OH termini? I'm prepared to test both systems with two distinct terminal groups.

ricalessandri May 26, 2023
Collaborator

@YanChen32 I've added the OH-terminated PVA to that PR, see the files.

YanChen32 · 2023-06-05T15:45:55Z

YanChen32
Jun 5, 2023
Author

Hi @ricalessandri and @fgrunewald,

I am writing to seek your guidance on some questions regarding the use of the build_file.bld for constructing PEO and PVA in confined spaces. I have encountered specific inquiries, and I appreciate any assistance you can provide.

Firstly, I construct PEO500.itp and PVA68.itp using the following commands:

polyply gen_params -f PEO.oplsaa.LigParGen.ff -seq OHter:1 PEO:500 OHter:1 -o PEO500.itp -name PEO

and

polyply gen_params -f PVA_CH3ter.oplsaa.LigParGen.ff -seq CH3a:1 PVA:68 CH3b:1 -o PVA68.itp -name PVA

Next, I would like to confine the polymer model within a narrow slit space: [-25.0305,-5.81625, 1.75]->[25.0305,-5.81625,3.25]. Therefore, I have written the following bld file:

[ molecule ]
; molname
PEO 0 20
[ rectangle ]
; resname start stop  inside-out  x  y  z    a  b  c
 OHter      1      1     in       25.0305 5.81625 2.5  25.0305 5.81625 0.75
 PEO        2    501     in       25.0305 5.81625 2.5  25.0305 5.81625 0.75
 OHter    502    502     in       25.0305 5.81625 2.5  25.0305 5.81625 0.75

[ molecule ]
; molname
PVA 20 69
[ rectangle ]
; resname start stop  inside-out  x  y  z    a  b  c
 CH3a      1      1     in       25.0305 5.81625 2.5  25.0305 5.81625 0.75
 PVA       2     69     in       25.0305 5.81625 2.5  25.0305 5.81625 0.75
 CH3b      70    70     in       25.0305 5.81625 2.5  25.0305 5.81625 0.75

And the system file:

#include "oplsaa.ff/forcefield.itp"
#include "PEO500.itp"
#include "PVA68.itp"

[ system ]
; name
Melt of PEO and PVA

[ molecules ]
; name  number
PEO 20
PVA 49

After modeling with the following command:

polyply gen_coords -p system.top -b build_file.bld -o melt.gro -box 50.061 11.6325 5.0 -name melt

I observed that the polymer is not completely confined within the designated space. Here is an image illustrating the issue

I have spent a considerable amount of time trying to address this problem but haven't been successful. I would greatly appreciate your assistance in resolving this issue. Thank you for your time and support.

12 replies

fgrunewald Jun 13, 2023
Maintainer

It looks like it should work. Can you share part of your reference gro file?

YanChen32 Jun 13, 2023
Author

Something like this，

93916
    1OHter  OA1    1   8.230   3.182   2.5
    1OHter   C2    2   8.183   3.203   2.5
    1OHter  HA3    3   8.221   3.172   2.5
    1OHter   H4    4   8.185   3.196   2.5
    1OHter   H5    5   8.182   3.246   2.5
    2PEO    C01    6   7.846   3.299   2.5
    2PEO    O02    7   7.827   3.313   2.5
    2PEO    C03    8   7.861   3.293   2.5
... ...
   69PVA    H0541342  12.678   2.281   2.5
   69PVA    H0641343  12.732   2.205   2.5
   69PVA    H0741344  12.651   2.178   2.5
   70CH3b    C141345  12.957   2.164   2.5
   70CH3b    H141346  12.952   2.122   2.5
   70CH3b    H241347  12.998   2.162   2.5
   70CH3b    H341348  12.958   2.196   2.5
50.061 11.6325 5.0

The upper and lower bounds of the optimized strcutrue are approximately ± 11

fgrunewald Jun 15, 2023
Maintainer

@YanChen32 honestly I don't exactly know. The fact that it remains a slab indicates that the flat-bottom potential does act. Otherwise it wouldn't remain a slab. My guess is that the density of the mixture is not compatible with the dimensions. That is you try to squeeze too much matter into this narrow slit. Ultimately, the LJ potential against which you work when squeezing an incompressible material will always win over a harmonic potential even if the force-constant is really high.

You could use gmx density to measure the density along z and check if that is approximately the amorphous density averaged from pure PVA and PEO. If the value is close it simply means you have too many polymer chains in the slit.

YanChen32 Jun 25, 2023
Author

Thank you so much! You were absolutely correct in identifying the root cause of the problem: the density setting was too high. I was able to successfully resolve the issue by simply increasing the thickness of the slit.

fgrunewald Jun 26, 2023
Maintainer

Perfect! If you don't have any further questions please feel free to mark this discussion as resolved

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

How to construct PEO melts at atomistic level #320

{{title}}

Replies: 2 comments 20 replies

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

How to construct PEO melts at atomistic level #320

YanChen32 May 16, 2023

Replies: 2 comments · 20 replies

fgrunewald May 17, 2023 Maintainer

Creating a polymer melt with OPLS FF

1 - generating an itp file

2 - building the initial melt

3 - preparing to run a simulation

Creating parameters for a custom polymer

Watch out for tacticity

YanChen32 May 23, 2023 Author

fgrunewald May 23, 2023 Maintainer

ricalessandri May 24, 2023 Collaborator

YanChen32 May 24, 2023 Author

ricalessandri May 26, 2023 Collaborator

YanChen32 Jun 5, 2023 Author

fgrunewald Jun 13, 2023 Maintainer

YanChen32 Jun 13, 2023 Author

fgrunewald Jun 15, 2023 Maintainer

YanChen32 Jun 25, 2023 Author

fgrunewald Jun 26, 2023 Maintainer

YanChen32
May 16, 2023

Replies: 2 comments 20 replies

fgrunewald
May 17, 2023
Maintainer

YanChen32 May 23, 2023
Author

fgrunewald May 23, 2023
Maintainer

ricalessandri May 24, 2023
Collaborator

YanChen32 May 24, 2023
Author

ricalessandri May 26, 2023
Collaborator

YanChen32
Jun 5, 2023
Author

fgrunewald Jun 13, 2023
Maintainer

YanChen32 Jun 13, 2023
Author

fgrunewald Jun 15, 2023
Maintainer

YanChen32 Jun 25, 2023
Author

fgrunewald Jun 26, 2023
Maintainer