求质量加权Schlitter熵Python计算代码及现有代码排查修正

Let's walk through fixing your Schlitter formula entropy calculation step by step—there are a few key issues in your code that are causing the abnormal results.

Key Issues in Your Original Code #

• Wrong constant for electron charge: You used Euler = 2.71828 (natural constant e) instead of the actual electron charge e = 1.602176634e-19 C. The Schlitter formula requires (e^2) (electron charge squared), not the square of the natural logarithm base—this is the most critical mistake throwing off your results.
• Incorrect matrix reshaping: np.resize can distort your data by repeating or truncating elements if the total count doesn't align perfectly. Since your 1407675×3 file has exactly (2055 \times 2055 = 4223025) elements, you should use np.reshape instead to safely convert the flat data into the 2055×2055 covariance matrix.
• Unnecessary imports: Libraries like pmx, numpy.testing, and math (you didn't use any of their functions here) add clutter without contributing to the calculation.
• Variable name conflict: You overwrote det (imported from numpy.linalg) with det = logdet, which is bad practice even if you don't use the original det function later.

Corrected Python Code #

import numpy as np

# Load the ASCII data file: 1407675 rows × 3 columns = 4223025 elements total
sigma_flat = np.loadtxt('covar_fit.dat').flatten()

# Reshape into the 2055×2055 covariance matrix (safe since element count matches)
cov_matrix = sigma_flat.reshape((2055, 2055))

# Physical constants (SI units)
hbar = 1.054571628e-34  # Joule-seconds
T = 310.0               # Kelvin (simulation temperature)
Kb = 1.3806504e-23      # Boltzmann constant (J/K)
e = 1.602176634e-19     # Electron charge (Coulombs)

# Calculate the scaled covariance matrix term in the Schlitter formula
scaled_cov = (Kb * T * e**2 / hbar**2) * cov_matrix
# Add identity matrix to form (I + scaled_cov)
m = np.eye(2055) + scaled_cov

# Compute log determinant (numerically stable for large matrices)
sign, log_det = np.linalg.slogdet(m)

# Calculate entropy using the Schlitter formula
entropy = 0.5 * Kb * log_det

print(f"Calculated entropy: {entropy} J/K")

Additional Notes #

• Why slogdet instead of det?: For large matrices like 2055×2055, the direct determinant can underflow or overflow floating-point numbers. slogdet computes the logarithm of the absolute determinant plus the sign, which is far more numerically stable.
• Symmetry check: Covariance matrices should be symmetric. You can add a quick check with np.allclose(cov_matrix, cov_matrix.T) to confirm your reshaping didn't introduce errors.
• Unit consistency: All constants use SI units, so the resulting entropy will be in J/K—adjust constants if you need to work with a different unit system.

内容的提问来源于stack exchange，提问作者Sree