Diazepam/xTB 速度初始化方法

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! initialize velocities uniformly
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

subroutine mdinitu(n,iat,velo,mass,Ekin)
   use xtb_setparam
   implicit none
   integer n
   real(wp) velo(3,n)
   integer iat(n)
   real(wp) mass(n)
   real(wp) Ekin
   real x,ranf

   real(wp) eperat,v,f,t,edum,f2
   integer i,j
   integer, allocatable :: iseed(:)

   call initrand

   eperat=Ekin/(3.0*n)

   do i=1,n
      f2=1
      if(iat(i).eq.1) f2=2
      v=sqrt(2*eperat/mass(i))
      f=1.0d0
      call random_number(x)
      if(x.gt.0.5)f=-1.0d0
      velo(1,i)=v*f*f2
      f=1.0d0
      call random_number(x)
      if(x.gt.0.5)f=-1.0d0
      velo(2,i)=v*f*f2
      f=1.0d0
      call random_number(x)
      if(x.gt.0.5)f=-1.0d0
      velo(3,i)=v*f*f2
   enddo

   call ekinet(n,velo,mass,edum)

   t=edum/(0.5*3*n*0.316681534524639E-05)

end subroutine mdinitu

   molmass=0
   do i=1,mol%n
      molmass=molmass+atmass(i)
      mass(i)=atmass(i)*amutoau
   enddo
   tmass=molmass*amutoau
   do i=1,mol%n
      if(mol%at(i).eq.1.and.set%md_hmass.gt.0) then
         mass(i)=dble(set%md_hmass)*amutoau
         ! k=nbo(1,i) ! atom to which H is bonded
         ! mass(k)=mass(k)-mass(i)+ams(1)*amutoau ! reduce by the increase H mass
      endif
   enddo
   molmass=molmass*amutokg

   nfreedom = 3._wp*real(mol%n,wp)  ! 注：可以理解为 nfreedom = 原子数目 * 3
   if(set%mdrtrconstr ) nfreedom = nfreedom - 6.0d0
   if(set%shake_md    ) nfreedom = nfreedom - dble(ncons)
   if(zconstr.eq.1) nfreedom = nfreedom - dble(iatf1)  ! fragment 1 in Z plane
   write(*,'('' # deg. of freedom  :'',i6  )')idint(nfreedom)

   Tinit=Tsoll  ! 注：实践发现 Tsoll 是 xtb 配置文件中 temp 量。
   if (equi) Tinit=0.3*Tsoll   ! slow equi 注：“slow equi"可能是指"slow equilibration”（慢速平衡）的简称
   
   if(.not.restart)then
      f=1
      if(.not.thermostat) f = 2
      edum=f*Tinit*0.5*kB*nfreedom
      ! 注：Tinit 暂时理解为初始温度
      ! 注：kB 在 src/mctc/mctc_constants.f90 定义为 3.166808578545117e-06，指的是 Boltzmann constant in Eh/K
      ! 注：nfreedom 暂时理解为原子数目 * 3
      
      call mdinitu(mol%n,mol%at,velo,mass,edum)
   else
      call rdmdrestart(mol%n,mol%xyz,velo)
   endif

      !> Ordinal numbers
      integer, allocatable :: at(:)

from typing import Union
from rdkit import Chem
import numpy as np


def xTB_initial_velocity(atom_list: Union[list, np.ndarray], t_init: Union[int, float], nvt=True) -> np.ndarray:
    """
    xTB_initial_velocity
    使用 xTB 的方法初始化速度
    :param atom_list: Union[list, np.ndarray]，体系原子序数表
    :param t_init: Union[int, float]，初始温度
    :param nvt=True，NVT系综
    :return: np.ndarray

    Author: Heqi Liu
    GitHub: https://github.com/metaphorme
    Time: 4 April 2024
    Licensed under the MIT License
    """
    
    if isinstance(atom_list, list):
        atom_list = np.array(atom_list)

    # 初始化常数
    amutoau = 1.660539040e-27 / 9.10938356e-31
    kB = 3.166808578545117e-06  # Boltzmann constant in Eh/K

    # 初始化原子质量
    atomic_number_atomic_mass = []
    for i in range(119):
        atomic_number_atomic_mass.append([Chem.Atom.GetMass(Chem.Atom(i)) * amutoau])
    atomic_number_atomic_mass = np.array(atomic_number_atomic_mass)  # 质量（单位：电子质量 Me），列表index索引代表第index个元素，第0个元素为0。

    nfreedom = len(atom_list) * 3  # 自由度
    f = 1 if nvt else 2  # nvt 系综

    edum = f * t_init * 0.5 * kB * nfreedom
    eperat = edum / (3.0 * len(atom_list))

    f = np.random.choice([-1, 1], size=(len(atom_list), 3))  # 均匀分布的 array
    f2 = np.where(atom_list == 1, 2, 1).reshape(-1, 1)  # 如果原子是 H，则 f2 = 2，否则 f2 = 1；并转换为列向量
    f2 = np.repeat(f2, 3, axis=1)  # 横向扩展 3 次

    masses = atomic_number_atomic_mass[atom_list]  # 原子序数对应质量
    masses = np.repeat(masses, 3, axis=1)  # 横向扩展 3 次

    v = np.sqrt(2 * eperat / masses)
    velocity = v * f * f2
    return velocity

# 使用 xTB 的速度初始化方法初始化 C2H4O2 在 300K 时的速度，NVT 系综
velocities = xTB_initial_velocities(atom_list=[8, 8, 6, 6, 1, 1, 1, 1], t_init=300, nvt=True)
np.set_printoptions(precision=14)  # xTB 速度显示为 14 位小数
print(velocities)

xtb AceticAcid.xyz --input omd.inp --omd --gfn 2