Unreasonable water temperatures

Dear all,

I am doing simulation about water phase change on solid surface. The system setup is shown below. The solid layer at the bottom is fixed to avoid penetration. The second to five layers are set as phantom atoms which act as heat source in the temperature ramp period. Spc/e model is employed for the water molecules.

image.png

Here’s the problem. In the initial equilibration period, I applied a Berendsen thermostat to both the solid wall and water. It works well as their temperature was stable around 300K. (The temperature history is shown below at the bottom. Among these, c_wtemp is for water temperature, c_stemp for solid temperature, c_ptemp for phatom atoms temperature.)

But then when I changed the water molecules to NVE ensemble with NVT for solid wall unchanged, the water temperature went up to near 400K!!! My intuition is that the water temperature would fluatuate around 300K.

Afterwards, Nose-Hoover thermostat with 300K was applied to both the water and solid copper atoms. Water dropped to around 300K while the solid temperature dropped to about 280K, which is also unreasonable.

Next, in the nonequilibrium period, the temperature of phantom atoms were set as 700K with langevin thermostat. Due to energy transfer, the water was also heated. However, its temperature went up to more than 700K which is unreasonale one more time.

I tried to equiblirate the system for a longer time but the results are still the same. I also compared my spc/e water model parameters with examples and literatures and nothing is wrong. My input script and forcefield settings files are attached for your reference. Data file was too big that I cannot upload it but can be generated using the system.in and molecule.spce files. It would be really appreaciated if anyone could point out my mistakes.

Thank you so much in advance.

Best regards,
Huaqiang

Time-averaged data for fix ave_equil

TimeStep c_wtemp c_wtempcom c_stemp c_ptemp v_wateng c_peng c_NV c_NL c_NV2 c_flux[1] c_flux[2] c_flux[3]

50000 299.77 299.584 299.736 295.325 -69334.9 -82737.4 4.725 22495.3 15.6844 0.00543093 -0.0031244 -0.00237987

100000 300.049 300.029 300 297.667 -69542.5 -82957.5 3.8156 22496.2 12.307 0.0165766 0.00826255 -0.000999008

150000 300.068 300.051 300.002 300.263 -69521.9 -82937.7 6.1218 22493.9 14.2596 0.00494952 0.0165203 0.000129047

200000 300.038 300.017 300.006 299.196 -69511.6 -82926.1 7.8786 22492.1 16.371 0.00941122 -0.00155875 -0.00386503

Berendsen 250000 300.05 300.03 300.007 297.886 -69512.2 -82927.2 9.3256 22490.7 17.2088 0.0165929 -0.00903798 0.00217158
thermostat 300000 300.055 300.037 300.009 299.689 -69506.6 -82921.9 9.7596 22490.2 18.5956 -0.00784755 0.0145744 -0.000636559

350000 300.066 300.039 300.001 296.428 -69497.8 -82913.5 9.0586 22490.9 17.0374 0.00597621 -0.00849835 -0.00269738

400000 300.049 300.035 300.007 296.388 -69490 -82905 9.6596 22490.3 18.1366 0.00616239 -0.00164509 -3.99688e-05

450000 300.04 300.023 300.007 300.527 -69486.3 -82900.8 9.702 22490.3 18.1164 -0.00902277 0.0089902 -0.0050541

500000 300.068 300.05 300 304.778 -69514 -82929.8 7.856 22492.1 16.2416 -0.000545675 0.00769304 0.00256577

550000 306.208 306.195 300.013 296.489 -68567 -82257.3 6.8224 22493.2 15.725 -0.00457548 0.000869905 -0.00236908

600000 319.577 319.563 300.032 297.126 -66511.9 -80799.9 8.1848 22491.8 19.058 0.00971271 0.00319303 -0.00510984

650000 331.64 331.623 300.05 293.201 -64637.6 -79465 11.017 22489 23.664 0.00194841 0.010399 -0.00753428

700000 342.133 342.119 300.069 294.596 -63040.4 -78336.9 16.6684 22483.3 29.6428 -0.00706581 0.00484214 -0.0161666

NVE 750000 351.522 351.504 300.074 293.518 -61591.7 -77307.9 28.13 22471.9 43.6686 0.000925071 0.00900138 -0.00943204
for water 800000 359.96 359.932 300.096 289.144 -60273.8 -76367.3 36.0242 22464 52.8628 -0.000535201 -0.00299176 -0.0126955

850000 368.77 368.743 300.099 287.976 -58922.5 -75409.9 47.03 22453 65.4018 0.0106733 0.00260604 -0.0208987

900000 375.974 375.962 300.099 288.718 -57847.9 -74657.4 59.08 22440.9 78.9062 -6.07862e-05 0.00291083 -0.017027

950000 382.925 382.907 300.117 284.156 -56768.7 -73889 65.7782 22434.2 87.5514 0.0159226 0.00179776 -0.0191682

1000000 389.239 389.218 300.127 283.897 -55785.5 -73188 79.118 22420.9 102.275 -0.00148225 -0.0015981 -0.0219305

1050000 322.78 322.681 242.271 231.73 -65737.5 -80168.7 69.707 22430.3 81.7444 0.0019977 0.00466713 -0.0181628

1100000 313.949 313.926 263.757 256.648 -67128.6 -81165 64.6534 22435.3 75.5494 0.00853491 0.00328957 -0.01307

1150000 310.47 310.447 272.588 265.257 -67704.6 -81585.5 56.1712 22443.8 66.1684 -0.0153683 -0.00205204 -0.00454909

1200000 307.541 307.517 280.645 274.327 -68178.8 -81928.7 52.609 22447.4 62.5066 0.000250151 0.0167834 -0.0137123

Nose-Hoover 1250000 306.536 306.515 282.983 275.428 -68382.9 -82087.9 37.656 22462.3 46.7558 0.000397961 0.0043129 -0.00136455
thermostat 1300000 306.844 306.829 282.445 277.133 -68370.9 -82089.6 26.246 22473.8 34.8958 0.00292211 -0.00711526 -0.00855466

1350000 306.513 306.489 283.005 281.201 -68454.1 -82158 19.4082 22480.6 29.1596 -0.00938407 -0.00159307 -0.00660982

1400000 306.63 306.598 283.089 279.495 -68449.8 -82159 19.957 22480 29.9068 0.00500791 0.0132423 -0.00333734

1450000 307.572 307.548 280.672 276.157 -68284.2 -82035.5 20.3876 22479.6 29.0646 0.00318251 0.00889593 -0.00646374

1500000 307.288 307.269 280.824 273.688 -68312 -82050.6 21.9588 22478 31.5498 -0.00396037 -0.00922864 -0.0122016

1550000 324.024 324.007 498.224 696.98 -65799.8 -80286.6 15.2142 22484.8 25.5126 -0.00780476 0.00521678 0.27386

1600000 367.557 367.536 562.839 699.453 -59238 -75671.2 22.1876 22477.8 36.4606 0.0215501 -0.0233665 0.171368

1650000 413.207 413.179 589.147 700.434 -52291.5 -70765.6 52.4516 22447.5 75.3834 0.00326484 0.0150883 0.14894

1700000 454.726 454.698 607.887 701.322 -45763.5 -66093.9 160.426 22339.6 191.479 0.0244065 -0.0150139 0.116178

1750000 491.284 491.236 619.181 699.939 -39415 -61379.9 367.288 22132.7 413.903 0.00735649 0.0126466 0.0949383

1800000 523.119 523.042 635.778 699.959 -33379.4 -56767.6 700.981 21799 768.455 0.0256375 0.00555611 0.088359

1850000 547.693 547.568 651.663 701.436 -27506.5 -51993.4 1255.4 21244.6 1415.05 0.0133475 0.00614807 0.0700585

1900000 565.969 565.526 658.215 700.102 -22028.5 -47332.5 2247.34 20252.7 2711.51 0.0201717 0.00447774 0.0586792

1950000 580.146 577.587 664.521 700.394 -16904.1 -42841.9 4262.92 18237.1 5073.95 0.000113979 0.0226446 0.0347155

2000000 589.097 585.088 664.922 700.974 -12503.5 -38841.5 6642.5 15857.5 7727.88 0.00878033 -0.00166608 0.0498048

2050000 598.097 596.367 668.039 700.27 -9125.22 -35865.6 8204.39 14295.6 9692.06 -0.0206725 -0.00254996 0.0305262

2100000 607.982 607.808 679.811 700.677 -6293.35 -33475.7 9201.41 13298.6 11184.4 -5.38116e-05 0.0154385 0.0229525

2150000 616.582 616.283 691.272 699.399 -3804.14 -31371 10893.5 11606.5 13077.8 0.000591747 -0.0168724 0.00488508

2200000 628.934 628.668 701.126 700.238 -1348.84 -29468 12095.4 10404.6 14567 0.0243336 0.00272982 0.00373849

2250000 646.225 645.831 701.537 700.209 750.215 -28142 12312.1 10187.9 14965.6 -0.00192003 0.0083245 -0.0104636

Nonequilbrium 2300000 656.71 656.489 706.098 700.092 1904.86 -27456.1 13159.3 9340.73 15437.2 0.0112473 -0.0075944 -0.00942512
period 2350000 672.569 672.398 716.428 700.39 3378.82 -26691.2 13515.7 8984.34 15494.5 -0.00502918 -0.0067747 -0.0143113

forcefield.system (1.05 KB)

run.in (5.04 KB)

system.in (2.29 KB)

molecule.spce (598 Bytes)

your inputs are needlessly complex and convoluted and thus difficult to read and figure out. for problems with a large data file, please always try to simplify the system and reproduce the issue with (much) fewer atoms (a compressed data file of a smaller system may fit into the mailing list size limits). if this is a principal problem, it should show with a smaller system just as well, and debugging and testing is much easier, faster and less wasteful. in fact, almost all MD simulation studies are best first prototyped with a smaller system before scaling up the the real problem set size. LAMMPS can make this very easy through the replicate command.

furthermore, the habit of defining or computing lots of constants as variables at the beginning and then using the variables throughout the input - while convenient in production simulations - make it difficult to read and understand and thus much more effort to debug.

also, your system.in file does not run, because it is missing a variable definition.

the first thing i would suggest to try is to reduce the timestep to 1fs or less. 2.0fs is rather aggressive for a system that is not yet well equilibrated and inhomogeneous. also, the fact that you are picking up kinetic energy gradually is an indication, that your energy conservation is insufficient. also, you have to watch out for the ordering of fixes when using fix shake.

axel.