Eastsheng's Wiki

LAMMPS计算TIP4P/ICE水温度分布时，低于系统平均温度

2024-03-28 18:37:28

[toc]

计算温度分布脚本

1
2
3
4
5
6
7
8
9
10
compute   ke  all ke/atom

variable    kB equal 1.3806504e-23    # [J/K] Boltzmann
variable    kcalpermol2J equal 4186.0/6.02214e23

variable temp atom c_ke*${kcalpermol2J}/(1.5*${kB})

compute tt all chunk/atom bin/1d z lower 0.1 units reduced
fix out_temp all ave/chunk 10 100 1000 tt v_temp file temp.dat

TIP4P/ICE力场时温度分布低于thermo的平均温度

  • z方向温度分布~173 K
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# Chunk-averaged data for fix T_out_relax and group file
# Timestep Number-of-chunks Total-count
# Chunk Coord1 Ncount v_temp
200000 10 1536
  1 0.05 154.066 172.915
  2 0.15 155.458 172.365
  3 0.25 145.723 172.41
  4 0.35 155.457 172.784
  5 0.45 154.917 173.494
  6 0.55 154.212 173.457
  7 0.65 156.178 173.831
  8 0.75 153.307 174.22
  9 0.85 155.066 174.519
  10 0.95 151.616 173.927

  • thermo温度~260 K
1
2
3
4
5
Step Density Temp Press Pxx Pyy Pzz PotEng KinEng TotEng 
       0   0.56728221          260   -1857.7524    1720.1132   -2612.0574    -4681.313   -6692.1318    792.83586   -5899.2959 
    1000    0.7154115    246.54399   -734.00106   -1037.3156   -1329.2594    164.57188   -6160.1241    751.80353   -5408.3206 
    2000    0.8533722    266.16031   -79.064947   -267.96115    108.42078   -77.654467   -6469.8195     811.6209   -5658.1986 
    3000   0.96464464    260.20011    1225.1164    296.83731    1470.1031    1908.4087   -6652.5907    793.44607   -5859.1446

计算其他体系时温度分布正常

  • 如300K下甲烷气体
1
2
3
4
5
6
7
8
9
10
11
12
13
14
# Chunk-averaged data for fix T_out_relax and group file
# Timestep Number-of-chunks Total-count
# Chunk Coord1 Ncount v_temp
200000 10 5000
  1 0.05 501.048 301.611
  2 0.15 507.387 298.271
  3 0.25 500.772 299.352
  4 0.35 495.566 295.747
  5 0.45 499.28 298.474
  6 0.55 497.176 303.487
  7 0.65 495.479 303.706
  8 0.75 497.552 300.502
  9 0.85 500.8 299.323
  10 0.95 504.94 300.134

可能原因排除

  • 是TIP4P/ICE力场的原因？

  • 是shake的原因？√

    • shake导致计算的温度分布低于平均温度

    • 取消掉fix shake后：

      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      # Chunk-averaged data for fix T_out_relax and group file
      # Timestep Number-of-chunks Total-count
      # Chunk Coord1 Ncount v_temp
      200000 10 1536
        1 0.05 157.316 258.756
        2 0.15 155.085 259.871
        3 0.25 157.547 261.236
        4 0.35 141.184 260.59
        5 0.45 162.571 260.871
        6 0.55 147.449 259.799
        7 0.65 150.527 258.551
        8 0.75 160.536 259.461
        9 0.85 146.207 259.575
        10 0.95 157.578 259.584

    • 但是，SPC、SPCE、TIP4P、TIP4P/ICE等都是刚性模型

解决方法

采用temp/chunk

  • If spatially binned chunks contain some number of water molecules and fix shake is used to make each molecule rigid, then you could calculate a temperature with six degrees of freedom (DOF) (three translational, three rotational) per molecule by setting adof to 2.0.

  • 修改后脚本

    1
    2
    3
    compute tt all chunk/atom bin/1d z lower 0.05 units reduced
    compute myChunk all temp/chunk tt temp adof 2.0
    fix out_temp all ave/chunk 100 1 100 c_myChunk[1] file temp.dat mode vector

  • 300K的水温度分布结果：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    # TimeStep Number-of-rows
    # Row c_tts[1]
    10000 20
    1 299.841
    2 299.94
    3 297.445
    4 297.264
    5 299.993
    6 295.925
    7 294.42
    8 300.127
    9 294.811
    10 296.218
    11 289.603
    12 297.217
    13 295.78
    14 294.731
    15 305.435
    16 302.281
    17 300.403
    18 303.721
    19 305.458
    20 306.914

直接在计算温度时修改自由度的影响

1
2
3
4
5
6
7
8
9
10
compute   ke  all ke/atom

variable    kB equal 1.3806504e-23    # [J/K] Boltzmann
variable    kcalpermol2J equal 4186.0/6.02214e23

variable temp atom c_ke*${kcalpermol2J}/(${kB})

compute tt all chunk/atom bin/1d z lower 0.1 units reduced
fix out_temp all ave/chunk 10 100 1000 tt v_temp file temp.dat

延申问题

  • 对于水/固体的温度分布仍然有问题（保持shake时）：
    • 如果采用compute ke/atom，并正常计算，则会低估水温度、固体正常；
    • 采用compute ke/atom，分别计算固体和水部分温度，并在计算水温度时修改自由度影响，。
    • 如果采用compute temp/chunk则会高估固体温度、水正常；
  • 可能的解决方法：选用柔性水分子力场；
  • 输出两个温度分布，一个采用compute ke/atom固体温度，一个采用compute temp/chunk输出水温度

参考：

Tags: LAMMPS TIP4P 温度分布
← 论文复现-机器学习预测甲烷水合物平衡条件 LBM/PALABOS：PALABOS介绍与安装 →