Fix widom particle insertions

I am new in Lammps and I am trying to use the “fix widom” command in order to perform an insertion of an oxygen molecule in a polyethylene (PE) sample.
For both oxygen and PE, the TraPPE-UA force field is used.
The model for oxygen uses 2 atoms (charge -0.113 e) and an additional massless atom (charge 0.226 e) located at the center of the mass in order to compensate for the charges.
Below you can find the script I am using for the implementation of the widom method:

units                   real
atom_style              full
boundary                p p p

pair_style              lj/cut 9.09
pair_modify             mix arithmetic tail yes
dielectric              1.0
special_bonds   lj 0.0 0.0 0.0 coul 0.0 0.0 0.0
bond_style      harmonic
angle_style     harmonic
dihedral_style  nharmonic


group                   mobile union all
timestep                1


molecule                molox O2trappe.molecule
fix                     3 mobile widom 1 100 0 19494 298.15  mol molox

#print WIDOM results

variable my_time equal time
variable mu_ex equal "f_3[1]"
variable udiff equal "f_3[2]"
variable vol equal "f_3[3]"
fix mu_ex_out all print 1 "${my_time} ${mu_ex} ${udiff} ${vol}" file mu_ex.dat screen no

#fix                    1 mobile nvt temp 298.15 298.15 1.0
#velocity               mobile create 298.15 1680506336 mom yes rot yes dist gaussian
dump                    1 all atom 1 PE_film_MDP_T_298K_p_1_atm.dump
dump_modify             1 scale no image yes
#dump                   2 all custom 1 PE_film_MDO_T_298K_p_1_atm.veldump vx vy vz
thermo_style    custom step etotal ke pe ebond eangle edihed eimp evdwl ecoul elong etail temp press vol density
thermo                  1
thermo_modify   flush yes
restart                 1 PE_film_MDP_T_298L_p_1_atm1.rst PE_film_MDO_T_298K_p_1_atm2.rst
run                     5
write_restart   PE_film_MDO_T_298K_p_1_atm.rst

My questions are:

  1. Do I need to use a thermostat for the polyethylene? Or “fix widom” is adequate"
  2. In the manual for the fix widom it is reffered that

    …This fix computes a global vector of length 3…The vector values are the following global cumulative quantities:
    1 = average excess chemical potential on each timestep
    2 = average difference in potential energy on each timestep
    3 = volume of the insertion region

How can I print more details? For instance, I would like to know to which positions each time an insertion is performed. Or I would like to know what is the van der waals energy of the system after the insertion of the oxygen. Can I have access to this info?
I can provide you the files i am using if it is needed.
Thank you in advance.
Kind regards,

Using fix widom is an “advanced topic” within LAMMPS. You are well advised to first get well familiar with running simulations using LAMMPS before moving to such an advanced method. This will avoid a lot of problems and you can answer most questions yourself.

A thermostat by itself does nothing. The more relevant question is: Do you want the other atoms to move between insertion attempts? If yes, then you need to add time integration (e.g. fix nve, fix nvt, fix rigid) depending on the properties of the system. A thermostat would either be included or works on top of the time integration (e.g. with fix nve).

You need to modify the source code and (re)compile LAMMPS.

The alternative would be to not use fix widom at all, but rather use the LAMMPS python module and implement fix widom as a python script.

Please note that LAMMPS will refuse running MD with massless atoms unless you disable the corresponding check in the source code and (re)compile.

Thanks a lot for your answer.
Actually, I have previous experience with GROMACS, that’s why I attempted widom insertion in Lammps, but, indeed, getting well familiar with LAMMPS will help me.
Starting with the thermostat, I don’t want the other atoms to move, so as I understand, I will use only fix widom.
Regarding the source code, I checked the file fix_widom.cpp, but to me, it seems simpler to implement fix widom using a python script, as you suggested.
Finally, I had this problem with massless atoms in the past, that’s why I chose a pretty small mass for the middle oxygen atom and now the simulation is running properly.

I am also using fix widom command recently, how to use python script to handle widom.
Thank you for your consideration and have a nice day!