Issue with Temperature Increase During Vapor Deposition in NVT Simulation

LAMMPS LAMMPS Beginners
,
1

Dear LAMMPS Community,

I am currently attempting to deposit vapor molecules onto a substrate using the “deposit” command within an NVT simulation. However, I have observed that the temperature of the vapor increases over time during the simulation.

I would like to know whether the deposited atoms are being added to the thermostat. If they are, could you please help me understand why the temperature is still increasing?

I have attached my input file for reference.

units metal
dimension 3
atom_style full
boundary p p f

read_data data.system

#########Smooth Surface#########
pair_style hybrid lj/cut 15 sw
pair_coeff * * sw Si.sw Si Si
pair_coeff 1 2 lj/cut 0.0052 3.19 15
pair_coeff 1 1 lj/cut 0.0018 3.19 15

bond_style harmonic
bond_coeff 1 20.33777008 1.400000

angle_style harmonic
angle_coeff 1 2.731939264 120.137800 # CA CA CA

dihedral_style charmm
dihedral_coeff 1 0.628779672 2 180 0 # X CA CA X

neighbor 0.3 bin
neigh_modify delay 0 check yes

variable lx equal lx
variable ly equal ly
variable lz equal lz

region pressbox block EDGE EDGE EDGE EDGE 50 150
region slab block EDGE EDGE EDGE EDGE 150 EDGE
region sbox block EDGE EDGE EDGE EDGE 0 EDGE
group liquid dynamic all region sbox
group pressvapor dynamic all region pressbox

group surface type 1
group drop id 1:120
group addatoms region slab

velocity all create 300 432423

compute atomPE all pe/atom
variable pot_eng atom c_atomPE
compute atomKE all ke/atom
variable atomEng atom c_atomPE+c_atomKE
compute peratom pressvapor stress/atom NULL
compute p all reduce sum c_peratom[1] c_peratom[2] c_peratom[3]
variable press equal -(c_p[1]+c_p[2]+c_p[3])/(lxly100)

minimize 0.0 1.0e-8 1000 100000

timestep 0.001
fix freeze surface setforce 0.0 0.0 0.0
fix center all momentum 1 linear 1 1 1 angular
thermo 1000
thermo_style custom step temp pe ke etotal enthalpy press v_press pxx pyy pzz pxy pxz pyz vol

fix 4 addatoms deposit 1000 2 10000 12345 region slab near 5.0 vx -1 1 vy -1 1 vz -1 1
fix zwalls all wall/reflect zlo EDGE zhi EDGE

compute mdtemp liquid temp
compute_modify mdtemp dynamic/dof yes
fix 1 liquid nvt temp 300.0 300.0 0.1
fix_modify 1 temp mdtemp

dump 4 all custom 10000 system.atom.eng id v_atomEng
dump_modify 4 sort id
dump 1 all custom 10000 system.lammpstrj id type x y z vx vy vz
dump_modify 1 sort id
run 1000000
write_restart system.restart
write_data system.data

Best regards,
Moid

2

Hi @Moid_Khan,

Please have a look at this post to help you formatting your question and post to make it easier to read.

Your input is complex and it is hard to understand what is going on without more details on what you are trying to do. There are a lot of potential source for your issue in your script. As I understand it, the full box is defined in the region sbox so all atoms are incorporated in the liquid group. You also do not say how you “observe” that the temperature of the vapor increases. The use of the setforce 0.0 0.0 0.0 fix combined with NVT integration is also something that looks fishy to me.

But again, these are only guesses based on the few information available to me. Also where does your model comes from? Is this hybrid model from a paper about Si and Ca? With intramolecular interactions?

3

Dear Germain,

Thank you very much for your response.

I am attempting to simulate a water droplet on a graphene like substrate at z=−0.001 in a vapor environment. To achieve this, I have defined a group “liquid” that includes all atoms except the substrate. I am running this simulation in LAMMPS (version 29 Sep 2021). However, I have noticed in the output file that the temperature of the system increases linearly as atoms are inserted through the deposit command. Once all atoms are inserted, the temperature stabilizes at the increased value.

I am using the setforce 0.0 0.0 0.0 fix in combination with NVT integration to observe the dynamics of the “liquid” while keeping the substrate (surface atoms) fixed, as you mentioned. I am not including the substrate atoms in the NVT time integration.

Below are a few lines from the output file where the second column shows the temperature:

Step Temp PotEng KinEng TotEng Enthalpy Press v_press  Volume 
       0      300   -639.32969     1845.641    1206.3113   -24328.546   -700.99036           -0     58362213 
   10000    300.68728   -625.62513    1849.8692    1224.2441   -24304.486   -700.82213 0.00050319577    58362213 
   20000    301.06554   -620.12946    1852.1963    1232.0669   -24304.871   -701.04747  0.017610689    58362213 
   30000     301.6025   -614.56895    1855.4998    1240.9308   -24294.972   -701.01907  0.014078395   58362213 
   40000    302.64903   -610.06667    1861.9382    1251.8715   -24280.007   -700.90858   0.53048035  58362213 
   50000    303.02418   -607.29429    1864.2461    1256.9518    -24261.75   -700.54685     1.057573   58362213 
   60000    303.95923   -605.69781    1869.9987    1264.3009   -24248.993   -700.39838    1.4380654   58362213 
   70000    304.39614   -603.07646    1872.6866    1269.6102   -24246.471   -700.47491    2.2036335   58362213 
   80000    305.08863   -603.26657    1876.9469    1273.6803   -24236.016   -700.29963    2.4627704  58362213

I hope this clarifies my question.

Best regards,
Moid

4

The model in the provided input has nothing to do with that. Why bring in SW Si model? How can you do such a simulation with (apparently) only 2 different atom types? This makes no sense. Also It was impossible to guess that there is a substrate surface below z=0 from the provided script.

This is a very old version. Consider upgrading.

Then why monitor the temperature of the whole system and not only the mdtemp value? Also, why are you giving an initial velocity to your substrate if it is not supposed to move?

Overall I think your input has several weird things going on and the best way to solve your issue is to take time to simplify and think on your model, your methodology and its implementation in your script.

5

Quite apart from the issues @Germain has raised, I don’t see how this behaviour is unusual.

Water’s enthalpy of condensation is just under 10 kcal/mol, and you have been apparently running a simulation which continuously condenses high pressure water vapour into a liquid droplet. All of that heat has to go somewhere. The fact that your total energy has only risen by 70 kcal/mol after all that shows that the thermostat is already working very hard to remove the excess heat. You can of course run a “tighter” thermostat if you’d like.

1 Like
6

Thank you for your feedback. Let me clarify a few points:

I am using a coarse-grained water model where one atom represents one water molecule. There are only two types of atoms in the simulation: type 1 is the substrate, and type 2 is the coarse-grained water model. The SW potential has been used for interactions between water molecules. “Si” just denotes the atom name to call the parameter from the Si.sw file.

I am not monitoring the temperature of the whole system, as you can see I am using “fix_modify 1 temp mdtemp” to compute the temperature of the liquid group.
The initial velocity assigned to the substrate is a standard practice to ensure a proper starting condition, even though the substrate itself is not meant to move.

Best,
Moid

7

Dear @srtee

Thank you very much for your response. I appreciate your insightful approach to solving the issue. As you mentioned, the excess heat should be dissipated into the thermal bath, but I’m still not reaching the desired temperature even with a 5*timestep coupling.

To investigate the inconsistency, I set up an empty simulation box and added water molecules using the “deposit” command. The temperature starts at 0 because there were no atoms at t=0 and remains at 0 over time, not reaching the temperature set in the NVT run.

The output file looks like this:

Step Temp PotEng KinEng TotEng Enthalpy Pxx Pyy Pzz Pxy Pxz Pyz Volume
       0            0            0           -0            0            0            0            0            0            0            0            0       490000
   10000            0  -0.41239565           -0  -0.41239565   -0.4558501    10.443606    9.6577471    2.9854012    2.0508643    2.1996175    3.9281301       490000
   20000            0   -2.9838889           -0   -2.9838889   -2.4387559     10.85628     18.70707    24.311416    3.9177672   -2.2437789    10.304582       490000
   30000            0   -9.5407791           -0   -9.5407791   -8.9372795    31.048025    25.580314    28.180007    6.2523544   -4.2157236    4.1245605       490000
   40000            0    -25.07975           -0    -25.07975   -33.839955    18.202811    9.7931186   -16.436661    -12.82674    16.854723   -9.7451211       490000
   50000            0   -54.259374           -0   -54.259374   -64.606824    18.318886    5.0451892     -9.09772     7.075048    -7.786005    2.7420807       490000
   60000            0   -93.011191           -0   -93.011191   -104.73747    16.368273    9.1744867    18.547345    18.866539    -2.101244   -18.038211       490000
   70000            0   -132.48906           -0   -132.48906   -141.97639    39.624061    33.397455     6.706831    19.536277   -13.460371    23.334785       490000
   80000            0   -163.74316           -0   -163.74316   -183.72051   -13.845382   -5.9511298    27.036284   -19.038243    2.7391047   -6.8187243       490000
   90000            0   -200.93216           -0   -200.93216   -221.52723    4.3676411   -1.3280351    24.369078    4.5684369    16.076015   -16.803354       490000
  100000            0   -248.81819           -0   -248.81819   -264.04308    34.687728    47.588331    33.346309   -20.721421    14.419523   -27.815409       490000

These are the two warnings I am getting in the output file:

WARNING: Temperature for fix modify is not for group all (src/fix_nh.cpp:1433)

WARNING: One or more dynamic groups may not be updated at correct point in timestep (src/fix_group.cpp:171)

In my understanding, “compute_modify mdtemp dynamic/dof yes” might not be working in my case.

Here, I am attaching the input files for this simulation
lammps.input (1.0 KB)
Si.sw (783 Bytes)

8

How do you want anyone to guess that you are using a modified version of SW potential, if you modify a file of parameters derived from Si, by keeping the filename and the atom names the same and saying nothing about it? The reason to keep such obvious error and confusion sources in your protocol is beyond my mind but, as the saying goes, you do you. Even the comments of the file you are using (which you posted in a separated answer) are wrong since they were only partly modified. :roll_eyes:

You modified the temperature of the fix but the temperature output in your thermo info is still the temperature of the full system, including the velocities of your substrate atoms. My remark was about adding the c_mdtemp value to the thermo_style command.

Except if you are using thermostat that depends on velocity (such as DPD), this is actually useless if your atoms do not move. Even so, your velocity distribution will actually be wrong if it is not allowed to equilibrate.

1 Like
9

I’m still not sure of the scientific value of this simulation. I still don’t know why a real-life water droplet would stay at room temperature while high-pressure vapour rapidly continuously condenses into it – and physics is the only reason why an MD simulation is useful.

Furthermore, all your simulations so far show that the thermostat is removing heat extremely efficiently. You may simply need to use a setpoint below 300K to keep the steady state temperature at 300K.

The main point of curiosity in your new results is that, if you look carefully, the thermo readout is showing zero kinetic energy despite your fix deposit command which should assign an initial random velocity. This might be undefined behaviour from fix nvt. I say “undefined behaviour” rather than a “bug” because the Nose-Hoover thermostat needs initial values for its equations of motion, and the temperature of a set of zero atoms is … undefinable.

Some things you can try (really, for your own curiosity):

As @Germain has suggested, track the output of your mdtemp compute, not just the default thermo temperature.

Track the thermostat’s removed heat using the econserve thermo keyword.

Use the fix langevin thermostat which does not depend on a “thermometer”. (Set the time constant to something like a hundred timesteps – without an actual physical analogue, like I said at the start, there really isn’t a “correct” value for that parameter.)

2 Likes
10

Dear @Germain

How do you want anyone to guess that you are using a modified version of SW potential, if you modify a file of parameters derived from Si, by keeping the filename and the atom names the same and saying nothing about it? The reason to keep such obvious error and confusion sources in your protocol is beyond my mind but, as the saying goes, you do you. Even the comments of the file you are using (which you posted in a separated answer) are wrong since they were only partly modified. :roll_eyes:

Thank you for your patience and for tolerating my poor description of the input file.
I used the SW potential for the mW water model, which is simply a set of modified parameters derived from Si. I apologize for not changing the file name, which caused some confusion and inconvenience for you.

You modified the temperature of the fix but the temperature output in your thermo info is still the temperature of the full system, including the velocities of your substrate atoms. My remark was about adding the c_mdtemp value to the thermo_style command.

Yes, that’s the point I was missing. Now it’s giving the correct temperature for the ‘liquid’ group.
I’m attaching the output file here:

Step Temp PotEng KinEng TotEng c_mdtemp Enthalpy Pxx Pyy Pzz Pxy Pxz Pyz Volume 
       0            0            0           -0            0            0            0            0            0            0            0            0            0       490000 
   10000            0  -0.41239565           -0  -0.41239565     280.9727   -0.4558501    10.443606    9.6577471    2.9854012    2.0508643    2.1996175    3.9281301       490000 
   20000            0   -2.9838889           -0   -2.9838889    288.48633   -2.4387559     10.85628     18.70707    24.311416    3.9177672   -2.2437789    10.304582       490000 
   30000            0   -9.5407791           -0   -9.5407791    312.12851   -8.9372795    31.048025    25.580314    28.180007    6.2523544   -4.2157236    4.1245605       490000 
   40000            0    -25.07975           -0    -25.07975    289.05429   -33.839955    18.202811    9.7931186   -16.436661    -12.82674    16.854723   -9.7451211       490000 
   50000            0   -54.259374           -0   -54.259374    274.45768   -64.606824    18.318886    5.0451892     -9.09772     7.075048    -7.786005    2.7420807       490000 
   60000            0   -93.011191           -0   -93.011191    314.25213   -104.73747    16.368273    9.1744867    18.547345    18.866539    -2.101244   -18.038211       490000 
   70000            0   -132.48906           -0   -132.48906    292.44031   -141.97639    39.624061    33.397455     6.706831    19.536277   -13.460371    23.334785       490000 
   80000            0   -163.74316           -0   -163.74316    300.86623   -183.72051   -13.845382   -5.9511298    27.036284   -19.038243    2.7391047   -6.8187243       490000 
   90000            0   -200.93216           -0   -200.93216    301.91386   -221.52723    4.3676411   -1.3280351    24.369078    4.5684369    16.076015   -16.803354       490000 
  100000            0   -248.81819           -0   -248.81819    313.77434   -264.04308    34.687728    47.588331    33.346309   -20.721421    14.419523   -27.815409       490000 
  110000            0   -275.57301           -0   -275.57301     302.1324    -293.3341    14.634389    29.207906    37.205507   -2.8102005   -1.3592484    4.5718559       490000 
  120000            0   -300.39822           -0   -300.39822    296.16712   -319.97723   -1.2218595    41.268167    17.997005   -2.4057149    4.0791402   -3.1887602       490000 
  130000            0   -321.31314           -0   -321.31314    300.01897    -350.8024   -26.280158   -26.577781    16.941902   -11.762189    9.7983646   -23.273869       490000 
  140000            0   -333.32864           -0   -333.32864    306.89969   -362.03263    9.5134787   -10.010571   -21.905712   -17.732873    -4.057041    21.482749       490000 
  150000            0   -342.28446           -0   -342.28446     300.8258   -358.10418     71.59044   -4.9765082     33.33534   -14.203522   -6.6767313    34.892186       490000

Once again, thank you very much for the suggestion.
However, I’m still puzzled as to why the system (only the deposited molecules in this case) shows a temperature and kinetic energy of 0, while the correct values are being printed for the group defined as the full system. I can observe the dynamics of the atoms, so this discrepancy is confusing to me.

11

This behaviour (rarely encountered) matches the thermo_style documentation:

By default [temperature-related quantities in the thermo output are calculated] by using a temperature compute which is created when LAMMPS starts up, as if this command had been issued:

compute thermo_temp all temp

See the compute temp command for details. Note that the ID of this compute is thermo_temp and the group is all. You can change the attributes of this temperature (e.g. its degrees-of-freedom) via the compute_modify command. Alternatively, you can directly assign a new compute (that calculates temperature) which you have defined, to be used for calculating any thermodynamic quantity that requires a temperature. This is done via the thermo_modify command.

and

The kinetic energy of the system ke is inferred from the temperature of the system with 1/2 kB T of energy for each degree of freedom. Thus, using different compute commands for calculating temperature, via the thermo_modify temp command, may yield different kinetic energies, since different computes that calculate temperature can subtract out different non-thermal components of velocity and/or include different degrees of freedom (translational, rotational, etc).

So that explains the output:

By default compute temp calculates the number of degrees of freedom of N atoms in D dimensions as D(N-1). Minus one for the centre of mass which is stationary for momentum-conserving dynamics (although … maybe that shouldn’t be assumed within compute temp).

When a 3D simulation starts with zero atoms in group all, the thermo_temp compute thus assumes that there are negative three degrees of freedom – I’m pretty sure this triggers internal code that outputs a temperature and kinetic energy of zero, rather than outright giving up with a divide-by-zero error. Since thermo_temp has not been told to dynamically update its DoFs, this state of affairs continues throughout the run, leading to what you observed.

Note that the temperature compute already internally sums up all kinetic energies to calculate the temperature. The thermo ke output reuses this information to save calculation time – for very large systems, each new kinetic energy parallel sum will invoke one round of communications across tens or hundreds of nodes.

1 Like
12

Dear @srtee

Thank you for bringing up the concern about the simulation setup.

This wasn’t meant to reflect real conditions; I simply wanted to magnify the effect I was observing by applying high pressure or a high deposition rate, allowing me to observe the behavior within a shorter simulation time.

As you can see below, the thermostat is working properly and is effectively removing or adding energy to the system.

Step Temp PotEng KinEng TotEng c_mdtemp Ecouple Econserve Enthalpy Pxx Pyy Pzz Pxy Pxz Pyz Volume
       0            0            0           -0            0            0            0            0            0            0            0            0            0            0            0       490000
   10000            0  -0.41239565           -0  -0.41239565     280.9727   -26.826386   -27.238782   -0.4558501    10.443606    9.6577471    2.9854012    2.0508643    2.1996175    3.9281301       490000
   20000            0   -2.9838889           -0   -2.9838889    288.48633    -56.87814   -59.862029   -2.4387559     10.85628     18.70707    24.311416    3.9177672   -2.2437789    10.304582       490000
   30000            0   -9.5407791           -0   -9.5407791    312.12851   -85.385212   -94.925991   -8.9372795    31.048025    25.580314    28.180007    6.2523544   -4.2157236    4.1245605       490000
   40000            0    -25.07975           -0    -25.07975    289.05429   -101.88318   -126.96293   -33.839955    18.202811    9.7931186   -16.436661    -12.82674    16.854723   -9.7451211       490000
   50000            0   -54.259374           -0   -54.259374    274.45768   -100.79428   -155.05365   -64.606824    18.318886    5.0451892     -9.09772     7.075048    -7.786005    2.7420807       490000
   60000            0   -93.011191           -0   -93.011191    314.25213   -89.479332   -182.49052   -104.73747    16.368273    9.1744867    18.547345    18.866539    -2.101244   -18.038211       490000
   70000            0   -132.48906           -0   -132.48906    292.44031   -66.191484   -198.68054   -141.97639    39.624061    33.397455     6.706831    19.536277   -13.460371    23.334785       490000
   80000            0   -163.74316           -0   -163.74316    300.86623   -50.967958   -214.71112   -183.72051   -13.845382   -5.9511298    27.036284   -19.038243    2.7391047   -6.8187243       490000
   90000            0   -200.93216           -0   -200.93216    301.91386   -24.028432   -224.96059   -221.52723    4.3676411   -1.3280351    24.369078    4.5684369    16.076015   -16.803354       490000
  100000            0   -248.81819           -0   -248.81819    313.77434    14.940947   -233.87725   -264.04308    34.687728    47.588331    33.346309   -20.721421    14.419523   -27.815409       490000
  110000            0   -275.57301           -0   -275.57301     302.1324    43.150717    -232.4223    -293.3341    14.634389    29.207906    37.205507   -2.8102005   -1.3592484    4.5718559       490000
  120000            0   -300.39822           -0   -300.39822    296.16712    68.783659   -231.61457   -319.97723   -1.2218595    41.268167    17.997005   -2.4057149    4.0791402   -3.1887602       490000
  130000            0   -321.31314           -0   -321.31314    300.01897    89.197601   -232.11554    -350.8024   -26.280158   -26.577781    16.941902   -11.762189    9.7983646   -23.273869       490000
  140000            0   -333.32864           -0   -333.32864    306.89969    100.31944    -233.0092   -362.03263    9.5134787   -10.010571   -21.905712   -17.732873    -4.057041    21.482749       490000

See the compute temp command for details. Note that the ID of this compute is thermo_temp and the group is all. You can change the attributes of this temperature (e.g. its degrees-of-freedom) via the compute_modify command.

compute_modify             thermo_temp   dynamic/dof   yes
fix                        1 liquid nvt temp 300.0 300.0 0.1
fix_modify                 1  temp thermo_temp

This is working. Output file is following

> Step Temp PotEng KinEng TotEng Ecouple Econserve Enthalpy Pxx Pyy Pzz Pxy Pxz Pyz Volume
>        0            0            0           -0            0            0            0            0            0            0            0            0            0            0       490000
>    10000    317.32356  -0.17545501    4.0607117    3.8852567   -27.224204   -23.338947    6.2771089     7.410838    9.3403339    6.7110905     3.508029   -5.1854798   -3.9186431       490000
>    20000    286.62435   -1.7542168     7.372772    5.6185551    -57.46734   -51.848785    9.4883188    9.2657267    16.124712    12.569016    5.8320695   -2.1967416   -9.9729194       490000
>    30000     276.3485   -8.0093941    10.680533    2.6711391   -84.929506   -82.258367    7.3522393    23.534989    13.725381    8.6576842   0.25817776    3.1333478     -10.9667       490000
>    40000    311.38201   -25.421418    16.059464   -9.3619538   -101.77955   -111.14151   -5.9587349    5.3261578    16.046703    12.010143  -0.32428887 -0.080267319    9.9565463       490000
>    50000    297.22031   -48.665886    19.170952   -29.494934   -107.37839   -136.87332    -25.39502     28.71399    3.6811303    7.8219409    13.653969     11.09787   -5.7700444       490000
>    60000    302.73068   -82.879441    23.439476   -59.439965   -98.443434    -157.8834   -64.079044   -4.1187759   -2.6018189   -38.785264    3.5539362   -20.008458    10.674036       490000
>    70000    296.66381   -125.78877    26.804418   -98.984355   -74.116475   -173.10083   -101.81881   0.12295793   -26.925753   -1.0011098     9.579856   -4.0386459    8.5765963       490000
>    80000    296.87049   -167.78344    30.660444     -137.123   -46.021366   -183.14436   -132.89252    11.779554    9.1195093    20.598664   -4.5334383   -12.589437   -5.1380494       490000
>    90000    293.04411   -199.87521    34.053153   -165.82206   -24.104255   -189.92632   -167.37768     11.53074   -12.502411   -14.287784    -7.244236   -1.7439608    13.367925       490000
>   100000    299.65575   -247.38264    38.694812   -208.68783    15.842943   -192.84489   -220.92945   -26.521709   -60.787186   -32.772122   -7.3035457    5.9448595    9.2022907       490000
>   110000     289.2472   -277.64849    37.350747   -240.29775    47.462956   -192.83479   -241.35494   -8.5985435    1.5338507   -3.3055313     10.30012   -18.826557    9.1295486       490000
>   120000    301.91749    -299.9041    38.986873   -260.91722    68.078296   -192.83893   -261.67418    6.9346657   -7.0649779   -7.2948499   -25.122295   -29.822728      2.71265       490000
>   130000    303.32964   -314.48831    39.169225   -275.31909    82.480226   -192.83886   -282.77087   -25.464176   -27.287258   -20.344926   -23.321589   0.93990769    29.818334       490000

Similarly, I must update other thermodynamic quantities (press, pe …).

Thank you very much for clarifying the issue in detail. I appreciate the efforts both you and @Germain have put into solving this problem.

1 Like
13

I am working to maintain constant vapor pressure in my simulation by substituting the following line in the provided LAMMPS input file (from my first message in this thread):

fix 4 addatoms deposit 1000 2 10000 12345 region slab near 5.0 vx -1 1 vy -1 1 vz -1 1

with this line:

fix 4 addatoms gcmc 1 10 10 2 456543 300.0 -2.46 3.0 pressure 0.005 region slab

However, the substitution doesn’t seem to work as expected. The simulation isn’t terminated, but there’s no output in the thermo or any other files, even when I attempt to dump output files at each step. Despite this, the simulation continues to run.

Below is the output log file I received after running the GCMC-substituted original input file:

dynamic group liquid defined
dynamic group pressvapor defined
47476 atoms in group surface
120 atoms in group drop
238 atoms in group addatoms
WARNING: Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions (src/pair.cpp:244)
WARNING: One or more dynamic groups may not be updated at correct point in timestep (src/fix_group.cpp:171)
WARNING: One or more dynamic groups may not be updated at correct point in timestep (src/fix_group.cpp:171)
WARNING: Fix gcmc using full_energy option (src/MC/fix_gcmc.cpp:483)

To find out the source of the error, I ran a simulation using GCMC for half of the simulation box for an LJ particle. That simulation worked properly.

Input file for LJ fluid with half box GCMC and full NVT

# GCMC for LJ simple fluid, with dynamics
# T = 2.0
# rho ~ 0.5
# p ~ 1.5
# mu_ex ~ 0.0
# comparable to Frenkel and Smit GCMC Case Study, Figure 5.8

# variables modifiable using -var command line switch

variable        mu index -1.25
variable        temp index 2.0
variable        disp index 1.0
variable        lbox index 50

# global model settings

units           lj
atom_style      atomic
pair_style      lj/cut 3.0
pair_modify     tail no # turn of to avoid triggering full_energy
boundary        p p f

# box

region          box block 0 ${lbox} 0 ${lbox} 0 ${lbox}
create_box      1 box

region          slab block 0 ${lbox} 0 ${lbox} 25 ${lbox}
# lj parameters

pair_coeff      * * 1.0 1.0
mass            * 1.0

# we recommend setting up a dedicated group for gcmc

group           gcmcgroup type 1

# gcmc

fix             mygcmc gcmcgroup gcmc 1 100 1 1 29494 ${temp} ${mu} ${disp} pressure 0.005 region slab
fix             1 all nvt temp 2.0 2.0 0.1
fix             zwalls all wall/reflect zlo EDGE zhi EDGE
# atom count

variable        type1 atom "type==1"
group           type1 dynamic gcmcgroup var type1
variable        n1 equal count(type1)

# averaging

variable        rho equal density
variable        p equal press
variable        nugget equal 1.0e-8
variable        lambda equal 1.0
variable        muex equal ${mu}-${temp}*ln(density*${lambda}+${nugget})
fix             ave all ave/time 10 100 1000 v_rho v_p v_muex v_n1 ave one file rho_vs_p.dat
variable        rhoav equal f_ave[1]
variable        pav equal f_ave[2]
variable        muexav equal f_ave[3]
variable        n1av equal f_ave[4]

# output
compute_modify  thermo_temp dynamic/dof yes
#compute         thermo_press all pressure thermo_temp
variable        tacc equal f_mygcmc[2]/(f_mygcmc[1]+${nugget})
variable        iacc equal f_mygcmc[4]/(f_mygcmc[3]+${nugget})
variable        dacc equal f_mygcmc[6]/(f_mygcmc[5]+${nugget})
compute_modify  thermo_temp dynamic/dof yes
thermo_style    custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav v_n1av
thermo          1000

dump            1 all custom 10 test.gcmc.lj.lammpstrj id type x y z vx vy vz
dump_modify     1 sort id

# run

run             1000000

The output of the above script is as follows:

Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav v_n1av
       0            0            0            0           -0            0        0            0            0            0            0            0            0            0
    1000    2.0739961 0.0030743144 -0.022606211    3.0944463     0.001504      188   0.95755403   0.96179802   0.99106256   0.00131592 0.0026341975    12.021816       123.34
    2000    2.3342505 0.0033075063 -0.006491215     3.481815     0.001432      179   0.95757449   0.96145439   0.98620337   0.00143392 0.0030825456    11.849581       138.86
    3000    2.2802809  0.003478746 -0.0090035573    3.4026989     0.001544      193   0.95767396   0.96297486   0.98659966   0.00159952 0.0034801434    11.630191       159.28
    4000    2.1393647  0.003633831  -0.01273557    3.1941212      0.00172      215   0.95859382   0.96179193   0.98561508   0.00169712 0.0036851333    11.510703       171.48
    5000    2.2091034 0.0042540784 -0.0061327643    3.3000186     0.001944      243   0.95965029   0.96074081   0.98543786   0.00183944  0.004097954    11.349656       187.51
    6000    2.2170991 0.0043717764 -0.014395972    3.3122927     0.001992      249   0.95832141   0.96123798    0.9856127    0.0020192 0.0045293704    11.162913       211.88
    7000     2.160243 0.0040297167 -0.012031101    3.2266342     0.001888      236   0.95817141   0.96167988   0.98551754   0.00207568 0.0045988313    11.106981       217.82
    8000    2.2457954 0.0047532472 -0.017596743     3.355981      0.00212      265   0.95749699     0.962152   0.98602446   0.00210512 0.0046866174    11.079109       222.15
    9000    2.3657757 0.0047493316 -0.0051657931    3.5346372     0.002024      253   0.95775656   0.96181846   0.98669033   0.00208464  0.004637359    11.098776       220.31
   10000    2.2723047 0.0049489149 -0.010956748    3.3960626       0.0022      275   0.95828042   0.96099142   0.98707625   0.00213952 0.0048525838    11.046637       225.93
   11000    2.2086118 0.0046713666 -0.012487406     3.300556     0.002144      268   0.95869052   0.96079813   0.98623812   0.00218368 0.0048409925    11.005383       232.68
   12000    2.2143607 0.0051171688 -0.012820724    3.3097625     0.002256      282   0.95878175   0.96043096   0.98598789   0.00222464 0.0047906312    10.968668       237.35
   13000    2.1015738 0.0051070322 -0.004834693    3.1418528       0.0024      300   0.95858722   0.96033104   0.98598775   0.00226688 0.0048525867    10.931029        241.8
   14000    2.1314644 0.0047510957  -0.01314817    3.1855704       0.0022      275   0.95892711   0.96004203   0.

I would appreciate any advice or insights on how to resolve this issue.

14

As per forum rules – this is a new issue; please start a new thread. :slight_smile:

15

Sure, I am raising a new issue with the subject Issues with GCMC and NVT in Vapor-Substrate Simulation in LAMMPS