Polymer/CO2 system not time integrating when using fix rigid/nve/small

Good afternoon,

I am trying to simulate a system of polymer mixed with some CO2 molecules in the 16Mar18 version of LAMMPS with 24 MPI threads. I’ve attached the input file and the simplified data file containing 1 polymer chain and 20 CO2 molecules. In short, when I apply the rigid/nve/small or rigid/small fixes, the LAMMPS output stops at Timestep: 2.0 without ever showing any thermodynamic data. In contrast, if I simply delete the CO2 from the system and remove the fix, thermodynamic data is produced quite easily. Are there any particular commands that I may be missing in my input file for a mixed system like mine to work correctly?

data.PEGE_CO2_min (110 KB)

in.PEGE_CO2_min (1.62 KB)

bill,

thanks for providing illustrative and complete input examples. that helps a lot in figuring out what might be a problem.

i am seeing two issues:

  1. your CO2 geometries are not “precise” enough, so i had to change the TOLERANCE setting in the rigid fix
    ​ source code and recompile ​
    for it to accept those as linear molecules.
  2. you have a rather high potential energy initial configuration and there seems to be some conflict with CO2 molecules in combination with the langevin settings, which then makes the time integrator become unstable.

​when i change the middle part of your input to do just one step without any progressions, everything looks like it is straightening out and then a stable MD seems to be following.
see below for the modified input and please give it a try.

axel.​

print .
print =====================================
print “Step2: Heating: NVE dynamics with Langevin thermostat”
print =====================================
print .

reset_timestep 0
timestep 2.0
fix fixshake poly shake 0.0001 20 0 b 3 7
fix fixnve poly nve
fix fixrigid co2 rigid/nve/small molecule
run 0
unfix fixrigid
fix fixlangevin poly langevin {itemp} {ftemp} 100.0 699483
fix fixrigid co2 rigid/nvt/small molecule temp {itemp} {ftemp} 100.0
thermo 100
thermo_style custom time temp pe etotal press density

Hello Axel,

Before I use the modified input script, I have a couple of questions. First, what do you mean “precise”? Is the CO2 angle not 180 degrees in my data file? If not, then would simply minimizing the molecule without any pair interactions turned on help in the future? Second, what did you change your TOLERANCE setting to in the fix rigid source code?

Thanks again for your help.

Hello Axel,

Before I use the modified input script, I have a couple of questions.
First, what do you mean "precise"?

​you may need to use more digits for the coordinates in the data file​.

Is the CO2 angle not 180 degrees in my data file?

​apparently not enough for the code to recognize it.​

If not, then would simply minimizing the molecule without any pair
interactions turned on help in the future?

​never tried that.​ it might work.

Second, what did you change your TOLERANCE setting to in the fix rigid
source code?

​form 1.0e-6 to 1.0e-5

axel.

Hi Axel,

So I took a look at your new version of the script once I recompiled LAMMPS with the new tolerance. It seems that your version certainly got LAMMPS to recognize that the CO2 molecules were indeed linear. However, it seemed that the temperature diverged quite badly. I did manage to combine your edits to input with some additions. (1) I added in a minimization step that removed a good chunk of any overlaps in the system. (2) I assigned velocities right off the bat for the desired final temperature. This seems to run stably for me. I’ve attached it below for anyone else who may encounter a similar issue in the future.

group co2 type 1 5
group poly subtract all co2

print .
print ==========================================
print “step1: minimization: 5000 steps CG Minimization”
print ==========================================
print .

dump dumptrj all atom 10 {molname}_min-cg.trj dump_modify dumptrj image yes scale yes thermo 100 thermo_style custom time temp pe etotal press lx ly lz density restart 10 {molname}1_sd.min {molname}2_sd.min min_style sd minimize 1.0e-10 1.0e-10 100000 1000000 min_style cg restart 10 {molname}1_cg.min {molname}2_cg.min minimize 1.0e-10 1.0e-10 100000 1000000 write_data data.{molname}_minimize
undump dumptrj

print =====================================
print “Step2: Heating: NVE dynamics with Langevin thermostat”
print =====================================
print .

reset_timestep 0
timestep 2.0
fix fixshake poly shake 0.0001 20 0 b 3 7
fix fixnve poly nve
fix fixrigid co2 rigid/nve/small molecule
run 0
unfix fixrigid
unfix fixnve

fix fixrigid co2 rigid/nvt/small molecule temp {ftemp} {ftemp} 100.0
fix fixnvt poly nvt temp {ftemp} {ftemp} 100.0
velocity all create {ftemp} 12355 run 0 velocity all scale {ftemp}
thermo 100
thermo_style custom time temp pe etotal press density


run 50000

Thanks for your help.