求助用MDanalysis计算温度分布计算出的结果很小

Flash001

大家好！我想用MDanalysis计算温度分布  基本思想就是切片   然后计算每片的总动能，然后根据公式计算温度分布  下面是我的代码：
import MDAnalysis as mda
import numpy as np

u=mda.Universe("system-NVT.tpr","system-NVT.trr")

# 获取系统边界（x轴最小和最大值）
x_min = u.atoms.positions[:, 0].min()
x_max = u.atoms.positions[:, 0].max()

# 计算每片的宽度
slice_width = (x_max - x_min) / 20

# 准备数据结构存储结果
slices_mass_energy = [{'mass': 0, 'total_kinetic_energy': 0.0,'velocity': np.zeros(3), 'atom_count': 0,'center_x': 0.0, 'water_molecule_count': 0} for _ in range(20)]

# 切换到最后一帧
u.trajectory[-1]
print('Time', u.trajectory.time)
for i in range(20):
        x_start = x_min + i * slice_width
        x_end = x_start + slice_width
        # 计算每个分片的中心点 x 坐标
center_x = (x_start + x_end) / 2

        # 输出当前分片范围及其中心点
print(f"Processing slice {i + 1}: x range = [{x_start:.3f}, {x_end:.3f}], center_x = {center_x:.3f}")
        # 选择x轴在指定范围内的原子
slice_atoms = u.select_atoms(f"prop x >= {x_start} and prop x < {x_end}")
        # 输出选中原子的数量
print(f"Number of atoms in slice {i + 1}: {len(slice_atoms)}")

        # 选择分片中的水分子数量
water_molecules = u.select_atoms(f"resname SOL and prop x >= {x_start} and prop x < {x_end}")
        water_molecule_count = len(water_molecules.residues)  # 获取水分子的数量（按残基计数）
print(f"Water molecule count in slice {i + 1}: {water_molecule_count}")

        if len(slice_atoms) > 0:
        # 获取原子的质量和速度
masses = slice_atoms.masses * 1.66053906660e-27  # 将道尔顿转换为千克
velocities = slice_atoms.velocities * 100   # 将速度从Å/ps转换为m/s
        # 输出原子的质量和速度
print(f"Original masses (amu): {slice_atoms.masses}")
         print(f"Converted masses (kg): {masses}")
         print(f"Velocities (m/s): {velocities}")

        # 计算每个原子的动能
kinetic_energies = 0.5 * masses * np.sum(velocities ** 2, axis=1)
         print(f"Kinetic_energies (J):{kinetic_energies}")



        # # 计算每片的总质量和总动能
total_mass = masses.sum()
         total_kinetic_energy = kinetic_energies.sum()

        # 将结果存储在数据结构中
slices_mass_energy['mass'] = total_mass
         slices_mass_energy['total_kinetic_energy'] = total_kinetic_energy
         slices_mass_energy['atom_count'] = len(slice_atoms)  # 获取每片的原子数目
slices_mass_energy['center_x'] = center_x  # 中心点x坐标
slices_mass_energy['water_molecule_count'] = water_molecule_count  # 获取每片的水分子数
else:
            print(f"No atoms found in slice {i + 1}")

        # 计算每个分片的温度
for i, slice_data in enumerate(slices_mass_energy):
        total_kinetic_energy = slice_data['total_kinetic_energy']  # 使用总动能
atom_count = slice_data['atom_count']  # 获取原子数
center_x = slice_data['center_x']
        water_molecule_count = slice_data['water_molecule_count']  # 获取水分子数
# 动态计算自由度
constrained_degrees_of_freedom = 3 * water_molecule_count  # 每个水分子有3个约束自由度
degrees_of_freedom = 3 * atom_count - constrained_degrees_of_freedom - 3  # 总自由度减去约束自由度和系统质心自由度
# 计算温度，使用总动能
Temperature = (2 * total_kinetic_energy) / (degrees_of_freedom * 1.38e-23) if atom_count > 0 else 0.0

            # 输出每片的总质量（kg）、总动能（J）、原子数目、温度（K）及中心点 x 坐标
print(f"Slice {i + 1}:")
        print(f"  Center X (Å): {center_x:.3f}")
        print(f"  Total Mass (kg): {slice_data['mass']:.3e}")
        print(f"  Total Kinetic Energy (J): {total_kinetic_energy:.3e}")
        print(f"  Atom Count: {atom_count}")
        print(f"  Water Molecule Count: {water_molecule_count}")
        print(f"  Temperature (K): {Temperature:.3f}")我问了几个人别人说代码乍一看没啥问题  但是计算出的结果很小  有没有大佬帮忙看一下！！！感谢！！

zdworld

为什么都换算成标准单位。。。用原来的单位直接KJ/mol不是很方便吗，你这一个1e-23下去误差直接炸了，小数除小数和小数减小数都会导致结果出问题

Flash001

sobereva 发表于 2024-9-3 18:44
如置顶的新社员必读贴、论坛首页的公告栏、版头的红色大字非常明确所示，求助帖必须在帖子标题明确体现出此 ...

感谢社长 下次注意  麻烦社长能帮忙找找问题吗  这个问题卡好几天了

sobereva

如置顶的新社员必读贴、论坛首页的公告栏、版头的红色大字非常明确所示，求助帖必须在帖子标题明确体现出此帖内容是求助或提问，并清楚、准确反映出帖子具体内容，避免有任何歧义和含糊性，仔细看http://bbs.keinsci.com/thread-9348-1-1.html。我已把你的不恰当标题 “Gromacs计算温度分布” 改了，以后务必注意

发帖时恰当排版，别人都很难注意到你要问什么