How to get proper density for Polymeric system?

Mahshid · November 5, 2023, 12:20am

Hello everyone,

I am using LAMMPS to simulate a polymeric system with 10 chains and 10 monomers. The system uses DREIDING force field, accompanied by water (F3C model) and Hydronium ions. I built the polymeric chain in Materials Studio and after energy minimization, I got the .pdb file and generated LAMMPS data file after generating .pdb file for the whole system, using “packmol”. I have prepared lammps input file, using the parameters reported in several papers. However, I am not getting the proper density for my system. According to the papers, conducting annealing process is necessary for getting the system to the proper density. However, following the same procedure, it seems impossible to get the correct density. To get to the proper density (1.8 dyne/cm3), I applied 5000 atm, using NPT ensemble. But when I remove the pressure, the density again reduces to 1.4-1.5 dyne/cm3.

Does anybody know how can I get the required density for my system?

Thank you very much

akohlmey · November 5, 2023, 12:37am

How can anybody know? You write that you did things “according to several papers”, but how can we know for certain that is the the case? All we know from your message that you claim you have applied 5000 atm pressure for some time, but we cannot know for how long and if that was done correctly. Your description suggests are rather complex workflows with multiple steps where information about the system is lost and needs to be recovered for a correct simulation.

Bottom line, if there is not enough information to (easily) reproduce your simulation and review everything that you did, there is no way to make suggestions for improvements or to identify flaws or errors and make suggestions for how to correct those.

Mahshid · November 5, 2023, 1:19am

Sorry for not mentioning the exact procedure. Here is what I have done during the simulation:

I have generated the polymer chain in Materials Studio. Then, I optimized the geometry after imposing DFT charges on the molecule. This polymer has 10 monomers.
The whole system is composed of 10 polymer chains, 100 hydronium molecules (to neutralize the system), 200 water molecules and 10 oxygen molecules. Packmol was used to build the whole system. Therefore, in total, there are 7840 atoms in the system.
Lammps data file was generated using VMD.
I prepared the lammps input file based on parameters of pairs, bonds, angles, and dihedrals reported in similar papers. For polymer chain, DREIDING force field is used and for water and hydronium, F3C water model and classical Hydronium model are used. LJ parameters are used for oxygen. lj/cut/coul/long pair style is used for polymer and lj/cut is used for other particles. Also, I used cutoff distance of 12 Angstrom.
After minimization, I run NVT at 300 K for 5 ns.
I use NPT at 300 K and 5000 atm for 5 ns, to get the proper density.
This is followed by several annealing processes, including NPT and NVT, each for 0.05 ns (as reported in papers).
Finally, I run NPT at 300 K and 1 atm for 1 ns to get the production run. But the density goes back again at this stage.

Here is the input file for lammps attached, if you need it. Also, if data file is required, I can send it to you, as I could not upload it here.

Thank you for your time,

echo both
units real
dimension 3
boundary p p p
atom_style full
read_data System_F3C.data
restart 50 restart

Pair_Coefficient …

#pair_style lj/cut/coul/long 10.0 8.0

pair_style hybrid lj/cut/coul/long 12 lj/cut 12

pair_coeff 1 1 lj/cut/coul/long 0.084400 3.883700 # 1 Cn
pair_coeff 2 2 lj/cut/coul/long 0.095100 3.898300 # 2 Cp
pair_coeff 3 3 lj/cut/coul/long 0.049600 3.395300 # 3 F
pair_coeff 4 4 lj/cut 0.010000 0.900000 # 4 HH
pair_coeff 5 5 lj/cut 0.010000 0.900000 # 5 HW
pair_coeff 6 6 lj/cut 0.184800 3.553200 # 6 OH
pair_coeff 7 7 lj/cut 0.095380 3.094000 # 7 OO
pair_coeff 8 8 lj/cut 0.184800 3.553200 # 8 OW
pair_coeff 9 9 lj/cut/coul/long 0.095700 3.404600 # 9 O_2
pair_coeff 10 10 lj/cut/coul/long 0.095700 3.404600 # 10 O_3
pair_coeff 11 11 lj/cut/coul/long 0.344000 4.030000 # 11 S

#Hydronium interactions
pair_coeff 1 4 lj/cut 0.029052 2.391850 # HH Cn
pair_coeff 2 4 lj/cut 0.030838 2.399150 # HH Cp
pair_coeff 3 4 lj/cut 0.022271 2.147650 # HH F
pair_coeff 4 7 lj/cut 0.030822 1.997000 # OO HH
pair_coeff 4 9 lj/cut 0.030935 2.152300 # HH O_2
pair_coeff 4 10 lj/cut 0.030935 2.152300 # HH O_3
pair_coeff 4 11 lj/cut 0.058652 2.465000 # HH S

pair_coeff 1 6 lj/cut 0.124888 3.718450 # Cn OH
pair_coeff 2 6 lj/cut 0.132569 3.725750 # Cn OH
pair_coeff 3 6 lj/cut 0.095740 3.474250 # F OH
pair_coeff 6 7 lj/cut 0.132499 3.323600 # OO OH
pair_coeff 6 9 lj/cut 0.132986 3.478900 # O_2 OH
pair_coeff 6 10 lj/cut 0.132986 3.478900 # O_3 OH
pair_coeff 6 11 lj/cut 0.252133 3.791600 # S OH

pair_coeff 4 6 lj/cut 0.042988 2.226600 # HH OH

#Water interactions
pair_coeff 1 5 lj/cut 0.029052 2.391850 # HW Cn
pair_coeff 2 5 lj/cut 0.030838 2.399150 # HW Cp
pair_coeff 3 5 lj/cut 0.022271 2.147650 # HW F
pair_coeff 4 5 lj/cut 0.001000 0.900000 # HW HH
pair_coeff 5 6 lj/cut 0.042988 2.226600 # HW OH
pair_coeff 5 7 lj/cut 0.030822 1.997000 # HW OO
pair_coeff 5 9 lj/cut 0.030935 2.152300 # HW O_2
pair_coeff 5 10 lj/cut 0.030935 2.152300 # HW O_3
pair_coeff 5 11 lj/cut 0.058652 2.465000 # HW S

pair_coeff 1 8 lj/cut 0.124888 3.718450 # Cn OW
pair_coeff 2 8 lj/cut 0.132569 3.725750 # Cp OW
pair_coeff 3 8 lj/cut 0.095740 3.474250 # F OW
pair_coeff 4 8 lj/cut 0.042988 2.226600 # HH OW
pair_coeff 6 8 lj/cut 0.184800 3.553200 # OH OW
pair_coeff 7 8 lj/cut 0.132499 3.323600 # OW OO
pair_coeff 8 9 lj/cut 0.132986 3.478900 # O_2 OW
pair_coeff 8 10 lj/cut 0.132986 3.478900 # O_3 OW
pair_coeff 8 11 lj/cut 0.252133 3.791600 # S OW

pair_coeff 5 8 lj/cut 0.042988 2.226600 # HW OW

#Oxygen interactions
pair_coeff 1 7 lj/cut 0.089543 3.488850 # Cn OO
pair_coeff 2 7 lj/cut 0.095050 3.496150 # Cp OO
pair_coeff 3 7 lj/cut 0.068644 3.244650 # F OO
pair_coeff 7 9 lj/cut 0.095349 3.249300 # O_2 OO
pair_coeff 7 10 lj/cut 0.095349 3.249300 # O_3 OO
pair_coeff 7 11 lj/cut 0.180776 3.562000 # S OO

#Nafion interactions
pair_coeff 1 2 lj/cut/coul/long 0.089590 3.891000 # Cn Cp
pair_coeff 1 3 lj/cut/coul/long 0.064701 3.639500 # Cn F
pair_coeff 2 3 lj/cut/coul/long 0.068680 3.646800 # Cp F
pair_coeff 1 9 lj/cut/coul/long 0.089873 3.644150 # Cn O_2
pair_coeff 2 9 lj/cut/coul/long 0.095400 3.651450 # Cp O_2
pair_coeff 1 10 lj/cut/coul/long 0.089873 3.644150 # Cn O_3
pair_coeff 2 10 lj/cut/coul/long 0.095400 3.651450 # Cp O_3
pair_coeff 1 11 lj/cut/coul/long 0.170392 3.956850 # Cn S
pair_coeff 2 11 lj/cut/coul/long 0.180871 3.964150 # Cp S
pair_coeff 3 9 lj/cut/coul/long 0.068896 3.399950 # F O_2
pair_coeff 3 10 lj/cut/coul/long 0.068896 3.399950 # F O_3
pair_coeff 3 11 lj/cut/coul/long 0.130623 3.712650 # F S
pair_coeff 9 10 lj/cut/coul/long 0.095700 3.404600 # O_2 O_3
pair_coeff 9 11 lj/cut/coul/long 0.181441 3.717300 # O_2 S
pair_coeff 10 11 lj/cut/coul/long 0.181441 3.717300 # O_3 S

#pair_modify mix arithmetic

pair_modify tail yes

#info coeffs out log

Bond_Coefficient …

bond_style harmonic

bond_coeff 1 429.32040 1.498200 # 1 Cn-Cn
bond_coeff 2 700.00000 1.530000 # 2 Cn-Cp
bond_coeff 3 605.25950 1.336000 # 3 Cn-F
bond_coeff 4 700.00000 1.530000 # 4 Cp-Cp
bond_coeff 5 605.25950 1.336000 # 5 Cp-F
bond_coeff 6 700.00000 1.420000 # 6 Cp-O_3
bond_coeff 7 700.00000 1.800000 # 7 Cp-S
bond_coeff 8 1085.9565 0.982000 # 8 HH-OH
bond_coeff 9 500.00000 1.000000 # 9 HW-OW
bond_coeff 10 119.02500 1.208000 # 10 OO-OO
bond_coeff 11 700.00000 1.480000 # 11 O_2-S

Angle_Coefficient …

angle_style harmonic

angle_coeff 1 106.2739 122.5536 # 1 Cn-Cn-Cn
angle_coeff 2 106.2739 122.5536 # 2 Cn-Cn-Cp
angle_coeff 3 100.3366 118.3191 # 3 Cn-Cn-F
angle_coeff 4 106.2739 122.5536 # 4 Cn-Cp-Cn
angle_coeff 5 100.3366 118.3191 # 5 Cn-Cp-F
angle_coeff 6 100.0000 109.4710 # 6 Cn-Cp-O_3
angle_coeff 7 100.3366 118.3191 # 7 Cp-Cn-F
angle_coeff 8 106.2739 122.5536 # 8 Cp-Cp-Cp
angle_coeff 9 100.3366 118.3191 # 9 Cp-Cp-F
angle_coeff 10 100.0000 109.4710 # 10 Cp-Cp-O_3
angle_coeff 11 100.0000 109.4710 # 11 Cp-Cp-S
angle_coeff 12 100.0000 125.2300 # 12 Cp-O_3-Cp
angle_coeff 13 350.0000 109.4710 # 13 Cp-S-O_2
angle_coeff 14 100.0000 109.4710 # 14 F-Cn-F
angle_coeff 15 100.0000 109.4710 # 15 F-Cp-F
angle_coeff 16 100.0000 109.4710 # 16 F-Cp-O_3
angle_coeff 17 100.0000 109.4710 # 17 F-Cp-S
angle_coeff 18 79.02630 113.4000 # 18 HH-OH-HH
angle_coeff 19 120.0000 109.4700 # 19 HW-OW-HW
angle_coeff 20 350.0000 115.5000 # 20 O_2-S-O_2

#Dihedral_coeficent …

dihedral_style harmonic

dihedral_coeff 1 6.43420 1 3 # 1 Cn-Cn-Cn-Cn
dihedral_coeff 2 2.00000 -1 3 # 2 Cn-Cn-Cn-Cp
dihedral_coeff 3 8.24440 1 3 # 3 Cn-Cn-Cn-F
dihedral_coeff 4 2.00000 -1 3 # 4 Cn-Cn-Cp-Cn
dihedral_coeff 5 2.00000 -1 3 # 5 Cn-Cn-Cp-F
dihedral_coeff 6 2.00000 -1 3 # 6 Cn-Cn-Cp-O_3
dihedral_coeff 7 2.00000 -1 3 # 7 Cn-Cp-O_3-Cp
dihedral_coeff 8 2.00000 -1 3 # 8 Cp-Cn-Cn-F
dihedral_coeff 9 2.00000 -1 3 # 9 Cp-Cp-Cp-F
dihedral_coeff 10 2.00000 -1 3 # 10 Cp-Cp-Cp-O_3
dihedral_coeff 11 2.00000 -1 3 # 11 Cp-Cp-O_3-Cp
dihedral_coeff 12 2.00000 -1 3 # 12 Cp-Cp-S-O_2
dihedral_coeff 13 8.08480 -1 3 # 13 F-Cn-Cn-F
dihedral_coeff 14 2.00000 -1 3 # 14 F-Cn-Cp-Cn
dihedral_coeff 15 2.00000 -1 3 # 15 F-Cn-Cp-F
dihedral_coeff 16 2.00000 -1 3 # 16 F-Cn-Cp-O_3
dihedral_coeff 17 2.00000 -1 3 # 17 F-Cp-Cp-F
dihedral_coeff 18 2.00000 -1 3 # 18 F-Cp-Cp-O_3
dihedral_coeff 19 2.00000 -1 3 # 19 F-Cp-Cp-S
dihedral_coeff 20 2.00000 -1 3 # 20 F-Cp-O_3-Cp
dihedral_coeff 21 2.00000 -1 3 # 21 F-Cp-S-O_2
dihedral_coeff 22 2.00000 -1 3 # 22 O_3-Cp-Cp-O_3
dihedral_coeff 23 2.00000 -1 3 # 23 O_3-Cp-Cp-S

#…

kspace_style pppm 1.0e-4
#kspace_modify gewald 0.1

#special_bonds lj 0.0 0.0 1.0
special_bonds dreiding

neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes

group H2O type 5 8

#group M type 1 2 4:11

fix freezewater H2O setforce 0.0 0.0 0.0

#minimize 1.0e-6 1.0e-6 100000 1000000

minimize 1.0e-8 1.0e-10 1000 10000
unfix freezewater

velocity all create 300 1287283

#fix 4 H2O shake 1.0e-4 100 0 b 6 a 12
#fix 3 all npt temp 300.0 300.0 10.0 iso 1.0 1.0 1000.0
#fix 1 all nve
fix 2 all nvt temp 300.0 300.0 100.0
#fix 2 all langevin 300.0 300.0 100.0 48279 scale 3 1.5

#Exports …

timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 10000000
undump 1
write_restart restart-1.equil
unfix 2

#…
fix 3 all npt temp 300.0 300.0 10.0 iso 5000.0 5000.0 1000.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 10000000
undump 1
write_restart restart-1.equil
unfix 3

#…

fix 4 all nvt temp 300.0 600.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 4
#…

fix 5 all nvt temp 600.0 300.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 5
#…

fix 6 all npt temp 300.0 300.0 10.0 iso 5000.0 5000.0 1000.0
timestep 0.5
thermo 100
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 6
#…

fix 4 all nvt temp 300.0 600.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 4
#…

fix 5 all nvt temp 600.0 300.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 5
#…

fix 6 all npt temp 300.0 300.0 10.0 iso 5000.0 5000.0 1000.0
timestep 0.5
thermo 100
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 6
#…

fix 4 all nvt temp 300.0 600.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 4
#…

fix 5 all nvt temp 600.0 300.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 5
#…

fix 6 all npt temp 300.0 300.0 10.0 iso 5000.0 5000.0 1000.0
timestep 0.5
thermo 100
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 6
#…

fix 7 all npt temp 300.0 300.0 10.0 iso 1.0 1.0 1000.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 2000000
undump 1
write_restart restart-1.equil
unfix 7

srtee · November 5, 2023, 4:07am

I can’t see why you would expect your system to have the same density under 1 atm of pressure as at 5000 atm of pressure, since it should compress under additional pressure (as a polymer would).

Any paper that claims to have observed that either has to explain it really well, or is probably unreliable.

Mahshid · November 5, 2023, 4:22am

Dear Srtee,

Thank you for your reply. Actually, my intention for increasing the pressure, is to get to the density of 1.8 g/cm3 for my system. But in the end, I need to run my system under atmospheric conditions but the same density (1.8 g/cm3).

When I run the simulation from the beginning at the atmospheric NPT ensemble, I hardly get to density of 0.7 g/cm3, which is way below what I should get. So, I increased pressure to build the system with the proper density. And I found that only at a pressure above 5000 atm, I can get to that density, which is not mentioned in the papers. In fact, in papers, they mostly emphasized the annealing process and it seems that they have all run the simulation at the atmospheric pressure.

In a nutshell, my target is to get to a density similar to the experimental works, and in order to get to that state under atmospheric conditions, I need to find the problem with my system.

Mahshid · November 5, 2023, 4:23am

Sorry for not mentioning the exact procedure. Here is what I have done during the simulation:

I have generated the polymer chain in Materials Studio. Then, I optimized the geometry after imposing DFT charges on the molecule. This polymer has 10 monomers.
The whole system is composed of 10 polymer chains, 100 hydronium molecules (to neutralize the system), 200 water molecules and 10 oxygen molecules. Packmol was used to build the whole system. Therefore, in total, there are 7840 atoms in the system.
Lammps data file was generated using VMD.
I prepared the lammps input file based on parameters of pairs, bonds, angles, and dihedrals reported in similar papers. For polymer chain, DREIDING force field is used and for water and hydronium, F3C water model and classical Hydronium model are used. LJ parameters are used for oxygen. lj/cut/coul/long pair style is used for polymer and lj/cut is used for other particles. Also, I used cutoff distance of 12 Angstrom.
After minimization, I run NVT at 300 K for 5 ns.
I use NPT at 300 K and 5000 atm for 5 ns, to get the proper density.
This is followed by several annealing processes, including NPT and NVT, each for 0.05 ns (as reported in papers).
Finally, I run NPT at 300 K and 1 atm for 1 ns to get the production run. But the density goes back again at this stage.

Here is the input file for lammps attached, if you need it. Also, if data file is required, I can send it to you, as I could not upload it here.

Thank you for your time,

echo both
units real
dimension 3
boundary p p p
atom_style full
read_data System_F3C.data
restart 50 restart

Pair_Coefficient …

#pair_style lj/cut/coul/long 10.0 8.0

pair_style hybrid lj/cut/coul/long 12 lj/cut 12

pair_coeff 1 1 lj/cut/coul/long 0.084400 3.883700 # 1 Cn
pair_coeff 2 2 lj/cut/coul/long 0.095100 3.898300 # 2 Cp
pair_coeff 3 3 lj/cut/coul/long 0.049600 3.395300 # 3 F
pair_coeff 4 4 lj/cut 0.010000 0.900000 # 4 HH
pair_coeff 5 5 lj/cut 0.010000 0.900000 # 5 HW
pair_coeff 6 6 lj/cut 0.184800 3.553200 # 6 OH
pair_coeff 7 7 lj/cut 0.095380 3.094000 # 7 OO
pair_coeff 8 8 lj/cut 0.184800 3.553200 # 8 OW
pair_coeff 9 9 lj/cut/coul/long 0.095700 3.404600 # 9 O_2
pair_coeff 10 10 lj/cut/coul/long 0.095700 3.404600 # 10 O_3
pair_coeff 11 11 lj/cut/coul/long 0.344000 4.030000 # 11 S

#Hydronium interactions
pair_coeff 1 4 lj/cut 0.029052 2.391850 # HH Cn
pair_coeff 2 4 lj/cut 0.030838 2.399150 # HH Cp
pair_coeff 3 4 lj/cut 0.022271 2.147650 # HH F
pair_coeff 4 7 lj/cut 0.030822 1.997000 # OO HH
pair_coeff 4 9 lj/cut 0.030935 2.152300 # HH O_2
pair_coeff 4 10 lj/cut 0.030935 2.152300 # HH O_3
pair_coeff 4 11 lj/cut 0.058652 2.465000 # HH S

pair_coeff 1 6 lj/cut 0.124888 3.718450 # Cn OH
pair_coeff 2 6 lj/cut 0.132569 3.725750 # Cn OH
pair_coeff 3 6 lj/cut 0.095740 3.474250 # F OH
pair_coeff 6 7 lj/cut 0.132499 3.323600 # OO OH
pair_coeff 6 9 lj/cut 0.132986 3.478900 # O_2 OH
pair_coeff 6 10 lj/cut 0.132986 3.478900 # O_3 OH
pair_coeff 6 11 lj/cut 0.252133 3.791600 # S OH

pair_coeff 4 6 lj/cut 0.042988 2.226600 # HH OH

#Water interactions
pair_coeff 1 5 lj/cut 0.029052 2.391850 # HW Cn
pair_coeff 2 5 lj/cut 0.030838 2.399150 # HW Cp
pair_coeff 3 5 lj/cut 0.022271 2.147650 # HW F
pair_coeff 4 5 lj/cut 0.001000 0.900000 # HW HH
pair_coeff 5 6 lj/cut 0.042988 2.226600 # HW OH
pair_coeff 5 7 lj/cut 0.030822 1.997000 # HW OO
pair_coeff 5 9 lj/cut 0.030935 2.152300 # HW O_2
pair_coeff 5 10 lj/cut 0.030935 2.152300 # HW O_3
pair_coeff 5 11 lj/cut 0.058652 2.465000 # HW S

pair_coeff 1 8 lj/cut 0.124888 3.718450 # Cn OW
pair_coeff 2 8 lj/cut 0.132569 3.725750 # Cp OW
pair_coeff 3 8 lj/cut 0.095740 3.474250 # F OW
pair_coeff 4 8 lj/cut 0.042988 2.226600 # HH OW
pair_coeff 6 8 lj/cut 0.184800 3.553200 # OH OW
pair_coeff 7 8 lj/cut 0.132499 3.323600 # OW OO
pair_coeff 8 9 lj/cut 0.132986 3.478900 # O_2 OW
pair_coeff 8 10 lj/cut 0.132986 3.478900 # O_3 OW
pair_coeff 8 11 lj/cut 0.252133 3.791600 # S OW

pair_coeff 5 8 lj/cut 0.042988 2.226600 # HW OW

#Oxygen interactions
pair_coeff 1 7 lj/cut 0.089543 3.488850 # Cn OO
pair_coeff 2 7 lj/cut 0.095050 3.496150 # Cp OO
pair_coeff 3 7 lj/cut 0.068644 3.244650 # F OO
pair_coeff 7 9 lj/cut 0.095349 3.249300 # O_2 OO
pair_coeff 7 10 lj/cut 0.095349 3.249300 # O_3 OO
pair_coeff 7 11 lj/cut 0.180776 3.562000 # S OO

#Nafion interactions
pair_coeff 1 2 lj/cut/coul/long 0.089590 3.891000 # Cn Cp
pair_coeff 1 3 lj/cut/coul/long 0.064701 3.639500 # Cn F
pair_coeff 2 3 lj/cut/coul/long 0.068680 3.646800 # Cp F
pair_coeff 1 9 lj/cut/coul/long 0.089873 3.644150 # Cn O_2
pair_coeff 2 9 lj/cut/coul/long 0.095400 3.651450 # Cp O_2
pair_coeff 1 10 lj/cut/coul/long 0.089873 3.644150 # Cn O_3
pair_coeff 2 10 lj/cut/coul/long 0.095400 3.651450 # Cp O_3
pair_coeff 1 11 lj/cut/coul/long 0.170392 3.956850 # Cn S
pair_coeff 2 11 lj/cut/coul/long 0.180871 3.964150 # Cp S
pair_coeff 3 9 lj/cut/coul/long 0.068896 3.399950 # F O_2
pair_coeff 3 10 lj/cut/coul/long 0.068896 3.399950 # F O_3
pair_coeff 3 11 lj/cut/coul/long 0.130623 3.712650 # F S
pair_coeff 9 10 lj/cut/coul/long 0.095700 3.404600 # O_2 O_3
pair_coeff 9 11 lj/cut/coul/long 0.181441 3.717300 # O_2 S
pair_coeff 10 11 lj/cut/coul/long 0.181441 3.717300 # O_3 S

#pair_modify mix arithmetic

pair_modify tail yes

#info coeffs out log

Bond_Coefficient …

bond_style harmonic

bond_coeff 1 429.32040 1.498200 # 1 Cn-Cn
bond_coeff 2 700.00000 1.530000 # 2 Cn-Cp
bond_coeff 3 605.25950 1.336000 # 3 Cn-F
bond_coeff 4 700.00000 1.530000 # 4 Cp-Cp
bond_coeff 5 605.25950 1.336000 # 5 Cp-F
bond_coeff 6 700.00000 1.420000 # 6 Cp-O_3
bond_coeff 7 700.00000 1.800000 # 7 Cp-S
bond_coeff 8 1085.9565 0.982000 # 8 HH-OH
bond_coeff 9 500.00000 1.000000 # 9 HW-OW
bond_coeff 10 119.02500 1.208000 # 10 OO-OO
bond_coeff 11 700.00000 1.480000 # 11 O_2-S

Angle_Coefficient …

angle_style harmonic

angle_coeff 1 106.2739 122.5536 # 1 Cn-Cn-Cn
angle_coeff 2 106.2739 122.5536 # 2 Cn-Cn-Cp
angle_coeff 3 100.3366 118.3191 # 3 Cn-Cn-F
angle_coeff 4 106.2739 122.5536 # 4 Cn-Cp-Cn
angle_coeff 5 100.3366 118.3191 # 5 Cn-Cp-F
angle_coeff 6 100.0000 109.4710 # 6 Cn-Cp-O_3
angle_coeff 7 100.3366 118.3191 # 7 Cp-Cn-F
angle_coeff 8 106.2739 122.5536 # 8 Cp-Cp-Cp
angle_coeff 9 100.3366 118.3191 # 9 Cp-Cp-F
angle_coeff 10 100.0000 109.4710 # 10 Cp-Cp-O_3
angle_coeff 11 100.0000 109.4710 # 11 Cp-Cp-S
angle_coeff 12 100.0000 125.2300 # 12 Cp-O_3-Cp
angle_coeff 13 350.0000 109.4710 # 13 Cp-S-O_2
angle_coeff 14 100.0000 109.4710 # 14 F-Cn-F
angle_coeff 15 100.0000 109.4710 # 15 F-Cp-F
angle_coeff 16 100.0000 109.4710 # 16 F-Cp-O_3
angle_coeff 17 100.0000 109.4710 # 17 F-Cp-S
angle_coeff 18 79.02630 113.4000 # 18 HH-OH-HH
angle_coeff 19 120.0000 109.4700 # 19 HW-OW-HW
angle_coeff 20 350.0000 115.5000 # 20 O_2-S-O_2

#Dihedral_coeficent …

dihedral_style harmonic

dihedral_coeff 1 6.43420 1 3 # 1 Cn-Cn-Cn-Cn
dihedral_coeff 2 2.00000 -1 3 # 2 Cn-Cn-Cn-Cp
dihedral_coeff 3 8.24440 1 3 # 3 Cn-Cn-Cn-F
dihedral_coeff 4 2.00000 -1 3 # 4 Cn-Cn-Cp-Cn
dihedral_coeff 5 2.00000 -1 3 # 5 Cn-Cn-Cp-F
dihedral_coeff 6 2.00000 -1 3 # 6 Cn-Cn-Cp-O_3
dihedral_coeff 7 2.00000 -1 3 # 7 Cn-Cp-O_3-Cp
dihedral_coeff 8 2.00000 -1 3 # 8 Cp-Cn-Cn-F
dihedral_coeff 9 2.00000 -1 3 # 9 Cp-Cp-Cp-F
dihedral_coeff 10 2.00000 -1 3 # 10 Cp-Cp-Cp-O_3
dihedral_coeff 11 2.00000 -1 3 # 11 Cp-Cp-O_3-Cp
dihedral_coeff 12 2.00000 -1 3 # 12 Cp-Cp-S-O_2
dihedral_coeff 13 8.08480 -1 3 # 13 F-Cn-Cn-F
dihedral_coeff 14 2.00000 -1 3 # 14 F-Cn-Cp-Cn
dihedral_coeff 15 2.00000 -1 3 # 15 F-Cn-Cp-F
dihedral_coeff 16 2.00000 -1 3 # 16 F-Cn-Cp-O_3
dihedral_coeff 17 2.00000 -1 3 # 17 F-Cp-Cp-F
dihedral_coeff 18 2.00000 -1 3 # 18 F-Cp-Cp-O_3
dihedral_coeff 19 2.00000 -1 3 # 19 F-Cp-Cp-S
dihedral_coeff 20 2.00000 -1 3 # 20 F-Cp-O_3-Cp
dihedral_coeff 21 2.00000 -1 3 # 21 F-Cp-S-O_2
dihedral_coeff 22 2.00000 -1 3 # 22 O_3-Cp-Cp-O_3
dihedral_coeff 23 2.00000 -1 3 # 23 O_3-Cp-Cp-S

#…

kspace_style pppm 1.0e-4

special_bonds dreiding

neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes

group H2O type 5 8

fix freezewater H2O setforce 0.0 0.0 0.0

minimize 1.0e-8 1.0e-10 1000 10000
unfix freezewater

velocity all create 300 1287283

fix 2 all nvt temp 300.0 300.0 100.0

#Exports …

timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 10000000
undump 1
write_restart restart-1.equil
unfix 2

#…
fix 3 all npt temp 300.0 300.0 10.0 iso 5000.0 5000.0 1000.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 10000000
undump 1
write_restart restart-1.equil
unfix 3

#…

fix 4 all nvt temp 300.0 600.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 4
#…

fix 5 all nvt temp 600.0 300.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 5
#…

fix 6 all npt temp 300.0 300.0 10.0 iso 5000.0 5000.0 1000.0
timestep 0.5
thermo 100
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 6
#…

fix 4 all nvt temp 300.0 600.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 4
#…

fix 5 all nvt temp 600.0 300.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 5
#…

fix 6 all npt temp 300.0 300.0 10.0 iso 5000.0 5000.0 1000.0
timestep 0.5
thermo 100
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 100000
undump 1
write_restart restart-1.equil
unfix 6
#…

fix 4 all nvt temp 300.0 600.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
run 100000
undump 1
write_restart restart-1.equil
unfix 4
#…

fix 5 all nvt temp 600.0 300.0 100.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
run 100000
undump 1
write_restart restart-1.equil
unfix 5
#…

fix 6 all npt temp 300.0 300.0 10.0 iso 5000.0 5000.0 1000.0
timestep 0.5
thermo 100
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
run 100000
undump 1
write_restart restart-1.equil
unfix 6
#…

fix 7 all npt temp 300.0 300.0 10.0 iso 1.0 1.0 1000.0
timestep 0.5
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 1000 System.xyz
#dump_modify 1 element CR F H O2C OR S
run 2000000
undump 1
write_restart restart-1.equil
unfix 7

srtee · November 5, 2023, 4:32am

Well if you have followed the papers’ procedures (it may help if you tell us what papers you are following) and do not achieve the same result, there are only two likely reasons:

You have made a mistake replicating what the papers said;
You did not make a mistake, but the papers’ procedures are unreliable and unrepeatable (or the authors made mistakes documenting their procedures).

Either way, it is not really a problem with LAMMPS or a problem we will be able to solve in this forum (unless it is a matter of no. 1 and there are some aspects of the papers’ procedures you don’t understand and would like help with.)

Mahshid · November 5, 2023, 4:59am

So, these are the several papers that I am following:

-I have used parameters, reported in paper #4, which is also referred by many other papers (I presume it should be reliable). Also, the following paragraph, is one of the procedures reported in papers (paper #1). I mostly have followed a similar procedure as reported in this paper (i.e. below), except for the fact that they never mentioned the pressure value.

“The initial cells are first optimized to minimize their energies by means of the steepest
descent method and the conjugate gradients minimization algorithms. After the optimization
process, the following dynamics simulation procedures are performed:
(a) MD simulation for 100 ps at 300 K in NVT ensemble.
(b) MD simulation for 100 ps at 300 K in NPT ensemble.
(c) The final structure of NPTMD simulation is heated from 300 K to 600 K,MD simulation
is performed for 50 ps for every 100 K increase in temperature in the NVT ensemble.
(d) The final structure of procedure (c) is cooled from 600 K to 300 K, MD simulation is
performed for 50 ps for every 100 K decrease in temperature in NVT ensemble.
(e) MD simulation for 100 ps at 300 K in NPT ensemble.
After repeating procedure (c) to (e) three times, the final structure was taken as the
initial configuration for equilibrating procedures. Equilibrations using this configuration
were carried out for 1 ns at 298 K, 323 K, and 353 K respectively, in the NPT ensemble. The
equilibrating procedure was followed by 1ns NVE dynamic simulation at 298 K, 323 K,
and 353 K, respectively, from which the trajectories were used for analyzing structural
properties and dynamic properties.”

-In addition, there is another procedure reported in paper #2, which I do not understand especially the section that is bolded here, as below. I do not get how they could change the LJ parameters during the simulation to change the density:

“To eliminate artificial voids and orientations of the molecules in the original models, the original models were subjected to high temperature and pressure relaxation and annealing treatments, re- spectively, in the initial stage of the simulation. For high tempera- ture and pressure relaxation processes, the l -J potential parameter εof Nafion atoms was reduced to the magnitude of 1/100, and MD simulation of 1 ns under the NPT ensemble was performed at 600 K and 100 MPa . Then the NPT ensemble under the condition of T = 300 K and P = 1 atm was used in the simula- tions and the parameters of l -J potential were slowly returned to the original values. After that, the annealing procedure was started. Under the NPT ensemble, the temperature of simulation systems was cycled five times between 300 K and 600 K and the anneal- ing rate was kept consistent at 1.5 K/ps for all simulated models to avoid introducing calculation deviations caused by structural dif- ferences at different annealing rates. It is worth noting that when the temperature of the systems reached the target tempera- ture, the NPT ensemble was switched to the NVT ensemble and the simulation was performed at this temperature for 0.2 ns to stabi- lize the configuration. Then the equilibration NPT MD simulation at T = 300 K and P = 1 atm was performed for 10 ns to reach the final equilibrium state.”

I appreciate your help and time very much,

srtee · November 5, 2023, 5:37am

The main advice I can give you now is that there is something very funny with your force field settings. A choice of

means that Coulomb, or charge, interactions are not considered for all pairs listed under lj/cut. Since those appear to be all interactions involving water molecules and hydronium ions (!!) it’s hard to see how your simulations up to now could have been valid in any way.

Once you have sorted that out, if you still do not have realistic densities for your system, you will simply have to find a force-field that is accurate for your system. The “Paper 1” on your list uses the COMPASS II force field, while the “Paper 4” on your list uses the DREIDING force field. It is impossible to know beforehand whether either of those force fields can accurately model your system.

And even if your system did finally match the experimental density, how will you know if its structure is correct? After all the density just tells you the size of the simulation box. The particles inside it could be arranged in any distribution – from water on the inside surrounded by a shell of polymers, to a polymer soup, to anything in between. So hopefully you have some concrete experimental data (like X-ray scattering) that will help you decide when you have reached a “correct” structure – especially since it looks like these systems have quite a history-dependent behaviour.

Bottom line, the simplest thing to do if you are fairly new to molecular dynamics is to follow a single paper closely (starting by trying to exactly replicate its results), and then to make small changes (making sure you know what each change does along the way) until you get to where you want. If you want to try doing anything else, that is advanced work which you will need to have a local mentor to help with.

akohlmey · November 5, 2023, 7:17am

Is this the procedure employed by DREIDING to get partial charges for the system? That seems rather questionable to me. I don’t know about DREIDING, but I know about several other force fields and there you need to either use some increment system (so no quantum mechanics at all) based on the atom typing or do some rather unsophisticated method (e.g. plain HF with a moderate to small basis set) since the force field was developed long before DFT was commonly used and they have to stick to those type of calculations to obtain parameters that remain consistent with the rest of the force field settings. Furthermore, the charges are determined with a specific method (there is not a single unique way to reduce a wavefunction into a set of point charges) and often combines the results of multiple calculations with some specific fitting procedure (e.g. RESP). This would then have to be followed by some validation calculations for the parameterized fragment only. and so on and so forth.

So what you describe as procedure represents rather a new parameterization and is not just following some paper.

This is very small and thus, especially since this concerns a polymer, you will have to contend with significant finite size effects and lack of proper phase space sampling since possible fluctuations to overcome barriers are smaller in smaller systems.
As you have already discovered, setting up a well equilibrated polymer system can be complex and difficult due to a phenomenon called “jamming”, i.e. the system gets “stuck” in a particular configuration since complex rearrangements would be needed to get out of that “hole” and those are rare events. Using high temperature and pressure temporarily is only the least sophisticated of strategies to get out of that. There are other approaches, often including some use of MC moves or through temporarily switching to a simpler, soft-core force field that allows atoms to move “through” each other to avoid the barriers between the (somewhat artificial) initial configuration and a proper equilibrated one.

This seems highly unusual.

This begs the more general question: what is your level of experience with MD simulations in general and polymers specifically? The description of your system suggests that it requires quite a lot of experience due to its complexity and non-conventional setup and polymer specific challenges in order to be performed in a meaningful way. And from your description it looks like you have put a significant effort into researching methods and models and have tried to adopt existing strategies from reputable sources, yet it still has many little bits and pieces that seem inconsistent and in conflict with general MD practices for such kind of systems. It is rarely straightforward to do such complex simulations without building a suitable level of experience, but a person with that level of experience would not ask the kind of questions that you ask. So please let us know.

That statement doesn’t make much sense to me. Either you have parameters that match your combination of atom types and that are parameterized for such type of molecule fragments or they are not. Using a mix and match approach is equivalent to a new force field creation and thus would require the detailed validation I mentioned already.

This seems to be a very small neighbor list skin for a system in real units. This looks more like it was taken from a simulation in reduced units.

Why this one?

What is a “classical Hydronium model”?

Which parameters? Taken from where? What about the O-O bond?

In general, I would have expected to see water, hydronium, and oxygen represented with rigid bonds, i.e. using either fix shake or fix rigid or similar. That is more common for an “old” force field like DREIDING (which goes back to the 1990s).

Further down the input and in general, I am sharing the concerns voiced by @srtee about the validity of your application of the force field and in particular skipping the charged interactions for (very) polar parts of your system. Please also note, that ignoring Coulomb interactions in the pairwise contributions will not cancel those in the kspace part. Furthermore, please note that a kspace convergence of 1.0e-4 is sufficient for converged forces but usually too loose for converged pressure.

Discussing the details of your equilibration protocol is rather moot at this stage with all the questions and concerns about your force field settings and procedure to obtain your starting configuration. If these initial parts are not representing your system well, then not even the most elaborate and sophisticated and complex equilbration procedure will get you a meaningful representation of the modeled system.

Mahshid · November 5, 2023, 8:51pm

Dear srtee,

Thank you very much for your reply. I have changed the pair style to lj/cut/coul/long for all the atoms. It has slightly improved the results but not significantly. And I am still struggling to get the proper density.

Also, I will try to use COMPASS force field to see if it will improve the results. Regarding getting the correct structure, I am not sure yet how can I justify it. But I was thinking of trying the same procedure for a larger system, which hopefully will prevent having unnecessary void spaces, though, I am not sure if it can get a more correct configuration.

I really appreciate your help and time,

Mahshid · November 5, 2023, 9:07pm

Dear akohlmey,

Many thanks for your response. I must confess that I am fairly new to the molecular dynamics field, yet the polymer simulations. It is only been about 8 months that I have been working in this area.

Also, regarding this question:

Which parameters? Taken from where? What about the O-O bond?

I have used the LJ parameter from this reference Molecular dynamics simulations of the surface tension of oxygen-supersaturated water | AIP Advances | AIP Publishing. Also, for bond interactions, I used 1.2080 angstroms for bond length and 11.025 Kcal/mol for bond energy.

Why this one?

I used this model because the majority of papers have used this model, the same as another model for Hydronium, which is called “classical hydronium model”, mentioned here https://doi.org/10.1021/jp036842c.
Also, I have used other water models, like tip3p, and the results have not made a significant change.

In addition, as you suggested, I will use a larger system to see if I can get a more reasonable result.

Thank you very much for your time and consideration.

srtee · November 5, 2023, 9:18pm

And have you adjusted the particle charges (in the data file) to match the force field of the water molecules and hydronium ions?

Once again, building a polymer configuration from scratch purely in the computer is a very difficult thing to accomplish. You should not be doing this without proper supervision. Of course, working on a computer, you won’t harm yourself or anyone around you with explosions or fires or anything like that. But you are at serious risk of wasting your time trying lots of unnecessary things, and in many ways time is the most precious resource of all.

If you are still learning how to be familiar with the basics of molecular dynamics, then you really should begin by taking the results of one paper and replicating it, before trying anything different. Please, especially do not try a bigger system, because you are likely to get the same inconclusive or wrong results while having spent even more time to get them.

(And finally – force fields will generally be more suitable to the systems they are meant to simulate and less suitable to other systems. You cannot take oxygen parameters from a simulation of super-saturated oxygenated water and re-use them to represent – polymer oxygens? The chemical environment of an atom controls how it responds physically – after all, you would not expect a benzene-ring carbon and a carbon inside a methane molecule to have exactly the same behaviour.)

Mahshid · November 7, 2023, 8:39pm

And have you adjusted the particle charges (in the data file) to match the force field of the water molecules and hydronium ions?

By this question, if you mean adjusting charges to get the neutral system, yes I did. In fact, each of my polymer chains has -10 charge and I have to add 10 Hydronium ions per chain to neutralize the system.

Also, I have another question. Is it possible to change the LJ parameters by two orders of magnitude and run NPT dynamics for some time? Then, again increase them to the original values for the rest of the simulation. This is one of the approaches that has been proposed to get the density and eliminate unnecessary voids in the simulation cell. But I am not sure how it can be possible without getting an error.