为什么 If 语句在蒙特卡洛模拟中不起作用？答案

【问题标题】：Why is If statement not working in Monte Carlo simulation?为什么 If 语句在蒙特卡洛模拟中不起作用？
【发布时间】：2022-06-10 18:11:47
【问题描述】：

我对 python 还很陌生，我正在尝试使用蒙特卡罗方法运行分子动力学模拟，其中我构建了一个对称系统并稍微扰乱一个随机粒子并计算系统的能量变化。我正在尝试实施一个 If 语句来拒绝概率上不可能的能量变化。但结果并不是拒绝不可能的系统。我附上了计算能量变化的代码，我做错了什么？结果，energy_new 变成了energy_initial，即使概率小于随机生成的数字，我不希望发生这种情况。

#Initialising positions

def initialise():
  
  global arr

  for i in range(nc):
    for j in range(nc):
      for k in range(nc):
        arr = np.append(arr,[i*a,j*a,k*a])
  
  arr.shape = (len(arr)//3,3)


#Calculating Energy
def Energy():
  global arr,L,rc
  ulj = 0
  rc2 = rc*rc
  for i in range(len(arr)-1):
    for j in range(i+1,len(arr)):
      dx,dy,dz = arr[i] - arr[j]
      dx,dy,dz = dx-L*round(dx/L),dy-L*round(dy/L),dz-L*round(dz/L) #Minimum Image convention
      r2 = dx*dx + dy*dy + dz*dz
      if r2 < rc2 and r2!=0:
        r6 = r2*r2*r2
        ulj += (1/r6)*(1/r6 - 1.0)
  ulj = 4*ulj
  return ulj


def loop():
      global Total_energy, count, arr
      energy_initial=Energy()   #initial energy
      print( "energy_initial=")
      print( energy_initial)
    
    #selecting and displacing a random particle
      random_particle = random.randint(0,len(arr)-1)
      
    #(random displacement) will be between -1 to 1
      _x = arr[random_particle,0] + random.random()-2
      _y = arr[random_particle,1] + random.random()-2
      _z = arr[random_particle,2] + random.random()-2
      displaced_particle = np.array([_x,_y,_z])
      arr_new=arr
      arr_new[random_particle]=displaced_particle
      
      energy_new=Energy()   #new energy
      print( "energy_new=")
      print(energy_new)
      
      dE= energy_new-energy_initial
      print( "dE=")
      print(dE)
    
      prob=math.exp(-beta*dE)   #calculating probability of move happening
      print("probability=")
      print(prob)
      random_no=random.uniform(0,1)
      print("random number")
      print(random_no)
    
      if (random_no > min(1.0,prob)):
        arr = arr
      else:
        Total_energy += energy_new
        count += 1
        arr = arr_new

【问题讨论】：

你需要调试你的代码。 This article 提供了一些很好的提示来帮助您入门。
我有一些注意事项。 arr_new = arr 不会创建数组的副本。它只是创建了对 SAME 数组的另一个引用。您对arr_new 所做的任何事情也将在arr 中可见。此外，您不需要Energy 中的全局语句。只有在为变量分配新值时才需要全局。
全局变量通常是个坏主意。函数应该接受其输入作为参数，并使用return 语句返回其结果。 initialize 应该创建 arr 并返回它，让调用者决定把它放在哪里。与loop 相同：arr 应该是输入，其他三个应该返回。
arr = arr 的意义何在？
也许先用print()看看你在变量中有什么以及执行了哪部分代码——它被称为"print debuging"

标签： python simulation physics montecarlo chemistry

【解决方案1】：

import numpy as np

a = np.ones((2,2))
b = a

b*=2
print("both matrices will be the same, since when you change b, you also change a")
print(a)
print(b)

a = np.ones((2,2))
b = np.copy(a)

b*=2
print("both matrices will be different, b is not identical to a")
print(a)
print(b)

正如 Tim Roberts 所指出的，您的行 arr_new = arr 将两个数组设置为相等。试试arr_new = np.copy(arr)

【讨论】：