[lammps-users] H2 GCMC with Furan Based Polymers Conundrum


I am attempting to perform gas solubility calculations of united-atom molecular hydrogen in a polymer system of PEF (polyethylene furanoate) using the GCMC fix in LAMMPS.

I am using the 3Mar20 version of LAMMPS and created the polymer structure and force field parameters using the OPLS forcefield and Moltemplate.
The polymer structure was equilibrated via an annealing procedure using the NPT ensemble.
I then appended FF parameters for UA H2 to the polymer parameter list and made a molecule file for insertions.

I’ve reached a roadblock with the mailing list archives in searching for answers to these weird results.

No matter what I try, whether it be pressure, temperature, tilmestep, insertion frequency, etc, there is almost never more than only a single molecule of H2 inside the polymer pores at one time which doesn’t feel right considering how small of a molecule it is.
The only thing I can think of is it may have something to do with the “intra_energy” option that I haven’t been using out of lack of understanding how to apply it in this situation based on the LAMMPS manual.

Below is an example of an input file I am using. I can provide the force field and data files if needed.

variables available on command line

variable mu index -6.9
variable disp index 0.5
variable temp index 298.15
variable press index 10.0

global model settings

units real
atom_style full
boundary p p p

include polymer-structure/system.in.init
read_data polymer-structure/system_after_min.data extra/atom/types 1
include polymer-structure/system.in.settings
include polymer-structure/system.in.charges
molecule h2mol h2-temp.txt


mass 907 2.01565

UA H2 model

pair_coeff 907 907 0.07293 2.958

MD settings

group h2 type 907
group ptf type 13 408 407 411 412 511 508 507 509 96 97 211 212 210
neighbor 2.0 bin
neigh_modify every 1 delay 10 check yes

fix fxnpt all npt temp {temp} {temp} 100 iso {press} {press} 1000

dynamically update fix nvt temperature ndof

fix_modify fxnpt dynamic/dof yes


variable tfac equal 1.5

fix mygcmc h2 gcmc 1000 100 0 0 68541 {temp} {mu} {disp} mol & h2mol tfac_insert {tfac} pressure ${press}

atom counts

variable hydrogen atom “type==907”
group hydrogen dynamic h2 var hydrogen
variable nH equal count(hydrogen)


variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)

dynamically update default temperature ndof

compute_modify thermo_temp dynamic/dof yes

dump 1 all custom 100 traj-h2ptfgcmc.lammpstrj id mol type x y z ix iy iz
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_nH
thermo 1000
timestep 1.0


run 10000
write_data system_after_min.data

Any and all help/advice would be greatly appreciated, thanks!
Micah Welsch
University of Kansas

I’m CCing Aidan (GCMC expert) to see if he has any ideas.