How to cross-link amorphous polymer chains

LAMMPS LAMMPS Beginners
,
1

Dears,

I want to obtain a cross-linked polyethylene (XLPE) system. First, I constructed an amorphous system with 3 polyethylene chains of 10 carbon atoms, under PCFF forcefield. And then converted it into the data file using msi2lmp tool. Then I attempted to cross-link the chains in LAMMPS by fix bond/create command. However, the chains seem to not cross-link with each other as the cumulative number of the bonds created remains 0 all the time (see the variable v_N1 in the log file). What are the problems that prevent the carbon atoms to cross-link and how to solve them? Thank you very much!

Besides, I also want to ask, when I convert the 3-chain polyethylene model exported from the Materials Studio into the data file, why there are 3 atom types (2 carbon atom types and 1 hydrogen atom type) rather than 2 atom types (1 carbon atom type and 1 hydrogen atom type)? And the two carbon atom types seem to have the same attributes.

The following are the input file and part of the data file.

input file:

cross-link of polyethylene

units real
newton on off
boundary p p p
atom_style full
neighbor 2.0 bin
neigh_modify delay 0 every 5 check yes

pair_style lj/class2/coul/long 10
bond_style class2
angle_style class2
dihedral_style class2
improper_style class2
kspace_style pppm 1e-4

read_data 3chainsPE.data

dump 1 all atom 1000 crosslink.lammpstrj
thermo 1000

velocity all create 300 428459 dist gaussian

timestep 1.0

fix 1 all npt temp 300.0 450.0 100.0 iso 1.0 150.0 1000.0
fix 2 all bond/create 1 1 3 1 2 prob 1 1423 # C3-C2
fix 3 all bond/create 1 1 1 1 4 prob 1 1423 # C3-C3
fix 4 all bond/create 1 3 3 1 4 prob 1 1423 # C2-C2

variable N1 equal f_2[2]+f_3[2]+f_4[2]
thermo_style custom step temp press etotal ebond epair v_N1

run 50000
write_data XLPE.data

part of the data file:

LAMMPS data file. msi2lmp v3.9.9 / 05 Nov 2018 / CGCMM for 3chainsPE

there are three pe chains in the amorphous system, with each chain having 10 carbon # atoms. 1-C3, 3-C2, 2-H

 96 atoms
 93 bonds
180 angles
243 dihedrals
120 impropers

3 atom types
4 bond types
7 angle types
7 dihedral types
6 improper types

 0.696054580     9.856854580 xlo xhi
 1.817578638    10.978378638 ylo yhi
-0.485494923     8.675305077 zlo zhi

Masses

1 12.011150 # c3
2 1.007970 # hc
3 12.011150 # c2

Pair Coeffs # lj/class2/coul/long

1 0.0540000000 4.0100000000 # c3
2 0.0200000000 2.9950000000 # hc
3 0.0540000000 4.0100000000 # c2

Bond Coeffs # class2

1 1.1010 345.0000 -691.8900 844.6000 # c3-hc
2 1.5300 299.6700 -501.7700 679.8100 # c3-c2
3 1.1010 345.0000 -691.8900 844.6000 # hc-c2
4 1.5300 299.6700 -501.7700 679.8100 # c2-c2

Angle Coeffs # class2

1 107.6600 39.6410 -12.9210 -2.4318 # hc-c3-hc
2 110.7700 41.4530 -10.6040 5.1290 # hc-c3-c2
3 107.6600 39.6410 -12.9210 -2.4318 # hc-c2-hc
4 110.7700 41.4530 -10.6040 5.1290 # c3-c2-hc
5 110.7700 41.4530 -10.6040 5.1290 # hc-c2-c2
6 112.6700 39.5160 -7.4430 -9.5583 # c3-c2-c2
7 112.6700 39.5160 -7.4430 -9.5583 # c2-c2-c2

Dihedral Coeffs # class2

1 -0.1432 0.0000 0.0617 0.0000 -0.1083 0.0000# hc-c3-c2-hc
2 0.0000 0.0000 0.0316 0.0000 -0.1681 0.0000# hc-c3-c2-c2
3 -0.1432 0.0000 0.0617 0.0000 -0.1083 0.0000# hc-c2-c2-hc
4 0.0000 0.0000 0.0316 0.0000 -0.1681 0.0000# hc-c2-c2-c2
5 0.0000 0.0000 0.0316 0.0000 -0.1681 0.0000# c3-c2-c2-hc
6 0.0000 0.0000 0.0514 0.0000 -0.1430 0.0000# c3-c2-c2-c2
7 0.0000 0.0000 0.0514 0.0000 -0.1430 0.0000# c2-c2-c2-c2

Improper Coeffs # class2

1 0.0000 0.0000
2 0.0000 0.0000
3 0.0000 0.0000
4 0.0000 0.0000
5 0.0000 0.0000
6 0.0000 0.0000

BondBond Coeffs

1 5.3316 1.1010 1.1010
2 3.3872 1.1010 1.5300
3 5.3316 1.1010 1.1010
4 3.3872 1.5300 1.1010
5 3.3872 1.1010 1.5300
6 0.0000 1.5300 1.5300
7 0.0000 1.5300 1.5300

BondAngle Coeffs

1 18.1030 18.1030 1.1010 1.1010
2 11.4210 20.7540 1.1010 1.5300
3 18.1030 18.1030 1.1010 1.1010
4 20.7540 11.4210 1.5300 1.1010
5 11.4210 20.7540 1.1010 1.5300
6 8.0160 8.0160 1.5300 1.5300
7 8.0160 8.0160 1.5300 1.5300

AngleAngle Coeffs

1 -0.3157 -0.3157 -0.3157 107.6600 107.6600 107.6600
2 0.2738 0.2738 -0.4825 107.6600 110.7700 110.7700
3 0.2738 -0.4825 0.2738 110.7700 107.6600 110.7700
4 0.2738 0.2738 -0.4825 107.6600 110.7700 110.7700
5 0.1184 -1.3199 -1.3199 110.7700 110.7700 112.6700
6 -1.3199 0.1184 -1.3199 110.7700 112.6700 110.7700

AngleAngleTorsion Coeffs

1 -12.5640 110.7700 110.7700
2 -16.1640 110.7700 112.6700
3 -12.5640 110.7700 110.7700
4 -16.1640 110.7700 112.6700
5 -16.1640 112.6700 110.7700
6 -22.0450 112.6700 112.6700
7 -22.0450 112.6700 112.6700

EndBondTorsion Coeffs

1 0.2130 0.3120 0.0777 0.2130 0.3120 0.0777 1.1010 1.1010
2 0.0814 0.0591 0.2219 0.2486 0.2422 -0.0925 1.1010 1.5300
3 0.2130 0.3120 0.0777 0.2130 0.3120 0.0777 1.1010 1.1010
4 0.0814 0.0591 0.2219 0.2486 0.2422 -0.0925 1.1010 1.5300
5 0.2486 0.2422 -0.0925 0.0814 0.0591 0.2219 1.5300 1.1010
6 -0.0732 0.0000 0.0000 -0.0732 0.0000 0.0000 1.5300 1.5300
7 -0.0732 0.0000 0.0000 -0.0732 0.0000 0.0000 1.5300 1.5300

MiddleBondTorsion Coeffs

1 -14.2610 -0.5322 -0.4864 1.5300
2 -14.8790 -3.6581 -0.3138 1.5300
3 -14.2610 -0.5322 -0.4864 1.5300
4 -14.8790 -3.6581 -0.3138 1.5300
5 -14.8790 -3.6581 -0.3138 1.5300
6 -17.7870 -7.1877 0.0000 1.5300
7 -17.7870 -7.1877 0.0000 1.5300

BondBond13 Coeffs

1 0.0000 1.1010 1.1010
2 0.0000 1.1010 1.5300
3 0.0000 1.1010 1.1010
4 0.0000 1.1010 1.5300
5 0.0000 1.5300 1.1010
6 0.0000 1.5300 1.5300
7 0.0000 1.5300 1.5300

AngleTorsion Coeffs

1 -0.8085 0.5569 -0.2466 -0.8085 0.5569 -0.2466 110.7700 110.7700
2 0.3113 0.4516 -0.1988 -0.2454 0.0000 -0.1136 110.7700 112.6700
3 -0.8085 0.5569 -0.2466 -0.8085 0.5569 -0.2466 110.7700 110.7700
4 0.3113 0.4516 -0.1988 -0.2454 0.0000 -0.1136 110.7700 112.6700
5 -0.2454 0.0000 -0.1136 0.3113 0.4516 -0.1988 112.6700 110.7700
6 0.3886 -0.3139 0.1389 0.3886 -0.3139 0.1389 112.6700 112.6700
7 0.3886 -0.3139 0.1389 0.3886 -0.3139 0.1389 112.6700 112.6700

Atoms # full

  1      1   1 -0.159000     6.604461002     3.586085686     5.364481687   0   1   0 # c3
  2      1   2  0.053000     6.894786610     4.022718417     6.330580459   0   1   0 # hc
  3      1   2  0.053000     7.314779609     3.950771599     4.612033539   0   1   0 # hc
  4      1   2  0.053000     5.606022917     3.969091269     5.104165507   0   1   0 # hc
  5      2   3 -0.106000     6.608589871     2.061913304     5.418187323   0   1   0 # c2
  6      2   2  0.053000     6.456460023    10.831237834     4.398340784   0   0   0 # hc
  7      2   2  0.053000     5.745224134    10.879364391     6.009828321   0   0   0 # hc
  8      3   3 -0.106000     7.903265530    10.647792740     6.005024093   0   0   0 # c2
  9      3   2  0.053000     8.048862342     1.863693803     7.033479610   0   1   0 # hc
 10      3   2  0.053000     8.760883737     1.871024730     5.428396743   0   1   0 # hc

…
…

2

That is a matter of the details of the force field. Atom are typed because of their environment and with molecular systems, this also impacts the derivation of bond, angle, dihedral, or improper terms.

There may be cases, where in one of those the difference between c2 and c3 types matter, even though it doesn’t for your specific molecule. Similarly, you need to keep in mind that when you are cross-linking your monomers, that atom types may change and thus also parameters (or not).

3

It looks like your carbon atoms can’t get close enough to each other to crosslink, given your fix bond/create length of 1. The atoms already bonded to a given carbon atom surround it at an equilibrium distance of 1 and change.

Besides that, how realistic is what you’re simulating? You’re trying to add a 5th bond to fully bonded carbon atoms.

4

When crosslinking polymers using class2 force fields, you may want to consider using ‘fix bond/react’ instead of ‘fix bond/create’ for more control over bonding changes and to correctly update angles, dihedrals, etc.

5

Thank you for your reply!

Actually, I have considered this reason. And I tried to increase the R_{min} of fix bond/create command to solve the problem. However, it did not work even if I set the R_{min} to 2. When I attempted to increase the R_min furtherly, the code went into errors, which seemed to be caused due to $R_{min} \ge$ the cutoff distance of the pairwise interaction. Should I increase the cutoff distance as well as the R_{min} to make the distance between carbon atoms closer than R_{min}? Or can I take some other actions to make the carbon atoms get closer during simulation, like changing the temperature or pressure by fix npt command? Could you please give me some advice to solve the problem?

My ultimate goal is to simulate the chemical reaction of XLPE polymer by fix reaxff command. The cross-linking of the PE is just an intermediate step. After the cross-linking, I will delete all the original settings about the forcefield in the data file, and then add the reaxff forcefield by potential file. Are they the right and appropriate procedures to realize my ultimate objective?

Thank you very much!

6

Thank you for your explanation!

How about the cross-link problem. Will the c2 and c3 types influence the cross-linking of my monomers or should I adjust the settings for fix bond/create command in the case of two carbon atom types?

This is just a test code so I only used 3 PE chains of 10 carbon atoms. Later I will construct a much larger PE system to cross-link. Thus, I used the Materials Studio rather than the LAMMPS itself for geometric modeling and chose the PCFF forcefield. Should I change the forcefield, like CVFF? What can I do to obtain a cross-linked polyethylene system? I want to simulate its chemical reaction in LAMMPS with reaxff forcefield.

Could you please give me some suggestions? Thank you very much!

7

Sorry, but this is not my area of research (way back when I did do some research). I have some ideas of how I would approach this, but without practical experience and checking them, I probably would do more harm than good by sharing those.

I suggest you study the workflow descriptions of similar publications and discuss with your adviser(s)/tutor(s)/mentor(s). Please keep in mind that for system preparation it doesn’t matter how realistic the system is during preparation for as long as you have a realistic and physically and chemically meaningful and correct system at the end. If you want to continue with ReaxFF, then it doesn’t matter much how accurate and detailed your force field is for as long as it gets you into the right ballpark.

Please note the comment by @jrgissing. Fix bond/create is more suitable for bead-spring models and united atom force fields (where the addition/removal of hydrogens is just a change in atom type). For more complex reactions where you need to delete hydrogens, you need a better tool like fix bond/react.

And note the commend from @Michael_Jacobs reminding you of creating a system with correct chemistry. If you don’t achieve that, your ReaxFF simulated system will “fall apart” or do other crazy things.

8

Thank you for your suggestion!

I read the introduction about the fix bond/react command in the manual. It seems to be much harder to use bond/react than to bond/create.

Do you mean that it is the complicated class 2 forcefield that increases the complexity of the simulation of polymer cross-linking, which leads to the failure of cross-link simulation? So can I still use the bond/create command to simulate the cross-linking if I choose a less complicated forcefield, like CVFF or something else?

By the way, during the cross-linking of the polyethylene, shall I have to break the bonds between carbon atoms and the hydrogen atoms first (e.g., with fix bond/break command), then the carbon atoms in the polyethylene chains can be cross-linked by fix bond/create command?

Look forward to hearing from you. Thank you very much!

9

Thank you for your precious suggestion!

Do you mean that I can directly connect some carbon atom pairs with the sketching tool in Materials Studio to cross-link polyethylene chains manually? Then I can also obtain a chemically meaningful cross-linked polyethylene structure (as long as the geometry optimization has been done) to implement ReaxFF simulation?

10

Questions that start like this always irritate me, since they mean that either somebody has not paid attention or is speculating beyond what I have stated. In the first case, you just need to read my suggestions again and follow them instead of discussing, in the second case you will have to construct a test case and see if it works out. That is just the nature of doing research. Adding my opinion or restating what I have already stated is not providing much help. I already pointed out, that I have little to none practical experience with this kind of research so any more specific comments from me would have a high risk of being inaccurate (or worse).

11

Sorry for my inappropriate question. And thank you very much for all your suggestions and help!