我正在编写代码来执行混合蒙特卡罗分子动力学。为此，我需要一个代码来根据麦克斯韦分布初始化所有粒子的速度。代码如下

   subroutine thermalize
        implicit none
        double precision :: MassInv
        double precision :: RATIO
        integer :: i

        !--- temp is temperature set point value
        STDtemp = sqrt(temp)
        do i = 1 ,np
           MassInv = 1.0d0/(mass(position(i)%type))
           MassInv = sqrt(MassInv*kboltz*kcal2amu)

  !---Gaussian() is a functions that returns gaussian random variable n(0,1)            
           v(i)%x = Gaussian()*STDtemp*MassInv*ratio
           v(i)%y = Gaussian()*STDtemp*MassInv*ratio
           v(i)%z = Gaussian()*STDtemp*MassInv*ratio


        enddo
      end subroutine thermalize

测量温度

T = 2.0*KE/(kboltz*(3*np - NCONS)

其中 KE 是动能，np 是原子数，kbotlz 是波尔兹曼常数。我试图在 temp=380 的情况下运行一些测试，并发现我的热化功能给了我初始动能，这些动能转化为 T=550 K 左右的温度。当我关闭摇动，因此 ncons = 0 时，产生的温度由thermoize 是 379。所以我想我需要修改 thermoize 命令以考虑受约束的自由度

为了试图解释抖动，在固定温度下存在抖动和不存在抖动的 KE 比率为

1/2kT*(3np-ncons)   = 1/2kT*(3np)

划分和重新排列

KE_shake/KE_noshake = (3np-ncons)/np

因此，我正在考虑乘以该比率的 sqrt。

v(i)%x = Gaussian()*STDtemp*MassInv*ratio*sqrt((3np-ncons)/np)
v(i)%y = Gaussian()*STDtemp*MassInv*ratio*sqrt((3np-ncons)/np)
v(i)%z = Gaussian()*STDtemp*MassInv*ratio*sqrt((3np-ncons)/np)

这给了我正确的温度（~380）。然而，当我热化我的粒子，并运行之前平衡流体的 NVE 轨迹（温度应该在 380 左右波动）时，温度下降到 300 左右。我不知道这在统计上是否正确。

有什么想法吗？