0

I'm fairly new to python and am trying to run molecular dynamics simulation using Monte Carlo method wherein I construct a symmetrical system and slightly perturb a random particle and calculate the change in energy of the system. I am trying to implement an If statement to reject probabilistically impossible energy change. But the outcome is not rejecting the improbable systems. I am attaching my codes for calculating energy change, what am I doing wrong? In the outcome, energy_new is becoming energy_initial even if the probability is smaller than the random generated number, which I don't want to happen.

#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
lleo10
  • 1
  • 1
  • 2
    You will need to debug your code. [This article](https://ericlippert.com/2014/03/05/how-to-debug-small-programs/) gives some great tips to get you started. – Code-Apprentice Feb 08 '22 at 20:17
  • I have some cautions. `arr_new = arr` does NOT create a copy of the array. It just creates another reference to the SAME array. Anything you do to `arr_new` will also be visible in `arr`. Also, you don't need the global statement in `Energy`. You only need global if you are assigning a new value to the variable. – Tim Roberts Feb 08 '22 at 20:22
  • 1
    Globals in general are a bad idea. A function should accept its inputs as parameters, and return its results with a `return` statement. `initialize` should create `arr` and return it, and let the caller decide where to put it. Same with `loop`: `arr` should be an input, and those other three should be returned. – Tim Roberts Feb 08 '22 at 20:23
  • What's the point of `arr = arr`? – Barmar Feb 08 '22 at 21:07
  • maybe first use `print()` to see what you have in variables and which part of code is executed - it is called `"print debuging"` – furas Feb 08 '22 at 21:16

1 Answers1

0
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)

As was pointed out by Tim Roberts your lines arr_new = arr sets both arrays equal. Try arr_new = np.copy(arr)

NopeRope
  • 1
  • 3