请问gromacs能计算材料的热导率吗？怎么计算呢？

Lisy

请问gromacs能计算材料的热导率吗？怎么计算呢？

yuzc

老老实实用LAMMPS吧，GROMACS做rNEMD麻烦死了。
脚本可以直接参考radonpy：

1. #########################################################
   
2. ## NEMD with kinetic energy exchange (RNEMD)
   
3. ##########################################################
   
4. fix             NVE all nve
   
5. fix             mp  all thermal/conductivity ${exchg} ${axis} ${slab}
   
6. 
   
7. # Generate temperature profile of layers
   
8. compute         layers all chunk/atom bin/1d ${axis} lower ${invslab} units reduced
   
9. fix             2 all ave/chunk ${Nevery} ${Nfreq} ${exchg} layers temp density/mass file ${Tprof} norm sample
   
10. 
    
11. # Output
    
12. dump            1 all custom 1000 ${dumpf} id type mol xs ys zs ix iy iz
    
13. dump            2 all xtc 1000 ${xtcf}
    
14. dump_modify     2 unwrap yes
    
15. restart         100000 ${rstf1} ${rstf2}
    
16. 
    
17. variable        heatflux   equal   (f_mp*${kcal2j}/${NA})/(2*${Jarea}*${ang2m}*${ang2m})   # J/m^2 = Ws/m^2
    
18. thermo_style    custom step time temp press enthalpy etotal ke pe ebond eangle edihed eimp evdwl ecoul elong etail vol lx ly lz density pxx pyy pzz pxy pxz pyz f_mp v_heatflux
    
19. thermo_modify   flush yes
    
20. thermo          ${exchg}
    
21. 
    
22. run             ${NStep}
    
23. 
    
24. 
    
25. """
    
26. 
    
27.         if decomp:
    
28.             in_strings += """
    
29. ##########################################################
    
30. ## Component decomposition of heat flux
    
31. ##########################################################
    
32. # heat flux preparation
    
33. compute         KE        all   ke/atom
    
34. compute         PE        all   pe/atom
    
35. 
    
36. compute         Spair     all   stress/atom NULL pair
    
37. compute         Sbond     all   stress/atom NULL bond
    
38. compute         Sangle    all   centroid/stress/atom NULL angle
    
39. compute         Sdihed    all   centroid/stress/atom NULL dihedral
    
40. compute         Simpro    all   centroid/stress/atom NULL improper
    
41. compute         Skspac    all   stress/atom NULL kspace
    
42. compute         Sfix      all   stress/atom NULL fix
    
43. """
    
44. 
    
45.             if decomp_intermol:
    
46.                 in_strings += """
    
47. compute         Spairer   all   stress/atom NULL interpair
    
48. compute         Spairra   all   stress/atom NULL intrapair
    
49. """
    
50. 
    
51.             in_strings += """
    
52. # Generate empty vector
    
53. group           empty     type      99999
    
54. compute         KENULL    empty ke/atom
    
55. compute         PENULL    empty pe/atom improper
    
56. compute         STNULL    empty stress/atom NULL improper
    
57. 
    
58. ########################   Cell half-left   ########################
    
59. ###  |//|  |  |  |  |**|  |  |  |  |  ###   |//| cold slab
    
60. ###  |//|  |  |  |  |**|  |  |  |  |  ###   |**| hot slab
    
61. ###  |//|  |  |  |  |**|  |  |  |  |  ###
    
62. ###  |//|  |  |  |  |**|  |  |  |  |  ###
    
63. ###  |//|  |  |  |  |**|  |  |  |  |  ###
    
64. ###      <--------->                  ###
    
65. ### heat flux decomposition of this reagion
    
66. #####################################################################
    
67. 
    
68. # left cell information
    
69. group           halfL     dynamic  all         region     lhalf     every  ${Nevery}  # Nevery=ave/time Nevery
    
70. 
    
71. #1st term   eivi
    
72. compute         lF1ke     halfL    heat/flux   KE         PENULL    STNULL
    
73. compute         lF1pe     halfL    heat/flux   KENULL     PE        STNULL
    
74. 
    
75. #2nd term   Sivi
    
76. compute         lFpair    halfL    heat/flux   KENULL     PENULL    Spair
    
77. compute         lFbond    halfL    heat/flux   KENULL     PENULL    Sbond
    
78. compute         lFangle   halfL    heat/flux   KENULL     PENULL    Sangle
    
79. compute         lFdihed   halfL    heat/flux   KENULL     PENULL    Sdihed
    
80. compute         lFimpro   halfL    heat/flux   KENULL     PENULL    Simpro
    
81. compute         lFkspac   halfL    heat/flux   KENULL     PENULL    Skspac
    
82. compute         lFfix     halfL    heat/flux   KENULL     PENULL    Sfix
    
83. """
    
84. 
    
85.             if decomp_intermol:
    
86.                 in_strings += """
    
87. compute         lFpairer  halfL    heat/flux   KENULL     PENULL    Spairer
    
88. compute         lFpairra  halfL    heat/flux   KENULL     PENULL    Spairra
    
89. fix             20 halfL ave/time ${Nevery} ${Nfreq} ${exchg}  c_lF1ke[${idx}] c_lF1pe[${idx}] c_lFpair[${idx}]  c_lFpairer[${idx}]  c_lFpairra[${idx}]  c_lFbond[${idx}]  c_lFangle[${idx}]  c_lFdihed[${idx}]  c_lFimpro[${idx}]  c_lFkspac[${idx}]  c_lFfix[${idx}]  file ${lJprof}
    
90. """
    
91.             else:
    
92.                 in_strings += """
    
93. fix             20 halfL ave/time ${Nevery} ${Nfreq} ${exchg}  c_lF1ke[${idx}] c_lF1pe[${idx}] c_lFpair[${idx}]  c_lFbond[${idx}]  c_lFangle[${idx}]  c_lFdihed[${idx}]  c_lFimpro[${idx}]  c_lFkspac[${idx}]  c_lFfix[${idx}]  file ${lJprof}
    
94. """
    
95. 
    
96.             in_strings += """
    
97. 
    
98. ########################   Cell half-right   #######################
    
99. ###  |//|  |  |  |  |**|  |  |  |  |  ###   |//| cold slab
    
100. ###  |//|  |  |  |  |**|  |  |  |  |  ###   |**| hot slab
     
101. ###  |//|  |  |  |  |**|  |  |  |  |  ###
     
102. ###  |//|  |  |  |  |**|  |  |  |  |  ###
     
103. ###  |//|  |  |  |  |**|  |  |  |  |  ###
     
104. ###                     <--------->   ###
     
105. ###        heat flux decomposition of this reagion
     
106. #####################################################################
     
107. 
     
108. # right cell information
     
109. group           halfR     dynamic  all         region     rhalf     every  ${Nevery}
     
110. 
     
111. #1st term   eivi
     
112. compute         rF1ke     halfR    heat/flux   KE         PENULL    STNULL
     
113. compute         rF1pe     halfR    heat/flux   KENULL     PE        STNULL
     
114. 
     
115. #2nd term   Sivi
     
116. compute         rFpair    halfR    heat/flux   KENULL     PENULL    Spair
     
117. compute         rFbond    halfR    heat/flux   KENULL     PENULL    Sbond
     
118. compute         rFangle   halfR    heat/flux   KENULL     PENULL    Sangle
     
119. compute         rFdihed   halfR    heat/flux   KENULL     PENULL    Sdihed
     
120. compute         rFimpro   halfR    heat/flux   KENULL     PENULL    Simpro
     
121. compute         rFkspac   halfR    heat/flux   KENULL     PENULL    Skspac
     
122. compute         rFfix     halfR    heat/flux   KENULL     PENULL    Sfix
     
123. """
     
124. 
     
125.             if decomp_intermol:
     
126.                 in_strings += """
     
127. compute         rFpairer  halfR    heat/flux   KENULL     PENULL    Spairer
     
128. compute         rFpairra  halfR    heat/flux   KENULL     PENULL    Spairra
     
129. fix             30 halfR ave/time ${Nevery} ${Nfreq} ${exchg}  c_rF1ke[${idx}] c_rF1pe[${idx}] c_rFpair[${idx}]  c_rFpairer[${idx}]  c_rFpairra[${idx}]  c_rFbond[${idx}]  c_rFangle[${idx}]  c_rFdihed[${idx}]  c_rFimpro[${idx}]  c_rFkspac[${idx}]  c_rFfix[${idx}]  file ${rJprof}
     
130. """
     
131.             else:
     
132.                 in_strings += """
     
133. fix             30 halfR ave/time ${Nevery} ${Nfreq} ${exchg}  c_rF1ke[${idx}] c_rF1pe[${idx}] c_rFpair[${idx}]  c_rFbond[${idx}]  c_rFangle[${idx}]  c_rFdihed[${idx}]  c_rFimpro[${idx}]  c_rFkspac[${idx}]  c_rFfix[${idx}]  file ${rJprof}
     
134. """
     
135. 
     
136.             in_strings += """
     
137. 
     
138. ##########################################################
     
139. ## RNEMD with kinetic energy exchange in decomposition
     
140. ##########################################################
     
141. thermo_style    custom step time temp press enthalpy etotal ke pe ebond eangle edihed eimp evdwl ecoul elong etail vol lx ly lz density pxx pyy pzz pxy pxz pyz f_mp v_heatflux
     
142. thermo_modify   flush yes
     
143. thermo          ${exchg}
     
144. 
     
145. run             ${NStepd}
     
146. 
     
147. 
     
148. """
     
149. 
     
150.         in_strings += """
     
151. write_dump      all custom ${ldumpf} id x y z xu yu zu vx vy vz fx fy fz modify sort id
     
152. write_data      ${ldataf}
     
153. quit
     
154. """

复制代码

Green-Kubo法：优点，准；缺点，需要很长的平衡时间；

1. ##########################################################
   
2. 
   
3. log             ${logf} append
   
4. 
   
5. units           real
   
6. atom_style      full
   
7. boundary        p p p
   
8. 
   
9. bond_style      ${bondst}  
   
10. angle_style     ${anglest}
    
11. dihedral_style  ${dihedst}
    
12. improper_style  ${improst}
    
13. 
    
14. pair_style      ${pairst} ${cutoff1} ${cutoff2}
    
15. pair_modify     mix arithmetic
    
16. special_bonds   amber
    
17. neighbor        2.0 bin
    
18. neigh_modify    delay 0 every 1 check yes
    
19. kspace_style    pppm 1e-6
    
20. 
    
21. read_data       ${dataf}
    
22. 
    
23. velocity        all create ${Ttemp} ${seed} mom yes rot yes dist gaussian
    
24. 
    
25. ##########################################################
    
26. ## Thermal conductivity calculation by Green-Kubo method
    
27. ##########################################################
    
28. timestep        ${TimeSt}
    
29. 
    
30. compute         kpKE      all ke/atom    # KE_i
    
31. compute         kpPE      all pe/atom    # PE_i
    
32. compute         kpStress  all centroid/stress/atom NULL virial  # S_i
    
33. compute         kpflux    all heat/flux kpKE kpPE kpStress
    
34. 
    
35. # x, y, z components of JE
    
36. variable        kpJx      equal c_kpflux[1]/vol
    
37. variable        kpJy      equal c_kpflux[2]/vol
    
38. variable        kpJz      equal c_kpflux[3]/vol
    
39. 
    
40. # Compute the autocorrelation function
    
41. fix             JJ all ave/correlate ${kpsample} ${kpcorrlen} ${kpdump} c_kpflux[1] c_kpflux[2] c_kpflux[3] type auto file ${autocorrf} overwrite ave running
    
42. 
    
43. variable        kpscale   equal ${conv}*(${kpsample}*dt)/${Ttemp}/${Ttemp}/vol/${kB}
    
44. variable        kappaxx   equal trap(f_JJ[3])*${kpscale}
    
45. variable        kappayy   equal trap(f_JJ[4])*${kpscale}
    
46. variable        kappazz   equal trap(f_JJ[5])*${kpscale}
    
47. variable        kappa     equal (v_kappaxx+v_kappayy+v_kappazz)/3.0   # in isotropic system, getting the average
    
48. fix             kappa     all ave/time ${kpdump} 1 ${kpdump} v_kappaxx v_kappayy v_kappazz v_kappa ave one file ${kappaf}
    
49. 
    
50. fix             NVT1 all nvt temp ${Ttemp} ${Ttemp} 100
    
51. 
    
52. # Output
    
53. dump            1 all custom 1000 ${dumpf} id type mol x y z vx vy vz
    
54. dump            2 all xtc 1000 ${xtcf}
    
55. dump_modify     2 unwrap yes
    
56. restart         100000 ${rstf1} ${rstf2}
    
57. 
    
58. thermo_style    custom step time temp press enthalpy etotal ke pe ebond eangle edihed eimp evdwl ecoul elong etail vol lx ly lz density pxx pyy pzz pxy pxz pyz v_kpJx v_kpJy v_kpJz
    
59. thermo          1000
    
60. 
    
61. run             ${NStep}
    
62. 
    
63. write_dump      all custom ${ldumpf} id x y z xu yu zu vx vy vz fx fy fz modify sort id
    
64. write_data      ${ldataf}
    
65. quit
    
66. """

复制代码