General Programming Language 가장 간단한 수치해석, essential example programs for physics [python]

컴퓨터 언어 파이썬을 활용한 수치해석 예제들을 나열해 보았다. 물리학 응용계산에서 기초적이면서도 근본적인 것들을 나열했다.

 

컴퓨터 언어로서 꼭 파이썬을 이용해야만 하는 것은 아니다. 파이썬을 이용하는 것이 많은 경우에 큰 도움이 된다.

매우 다양한 라이버러리들이 이미 개발되어 있다. 이것들을 철저하게 이용하는 것이 파이썬 언어를 활용한 프로그래밍 방법이다. 이것은 검색을 통해서 이루어 진다. 검색의 비중이 매우 높은 컴퓨터 언어가 파이썬 언어이다.

 

파이썬 언어의 철학은 다음과 같다. 문제를 풀 때, 그 문제를 푸는 가장 간단한 하나의 방법이 있다는 것이다. 

다양한 방법으로 문제를 푸는 것을 추구하지 않는다. 가장 간단한 방법으로 해당 문제를 풀고자 한다. 

따라서, 고도로 발달된 특정 문제 풀이는 결국 가장 간단한 하나의 프로그램으로 수렴하게 된다.

지속적으로 개발하면, 결국, 그렇게 된다는 것이다. 프로그램 길이는 계속해서 짧아진다.

검색을 지속적으로 하면, 계속해서 컴퓨터 프로그램을 간단하게 만들 수 있다.

놀라울 정도로 많은 라이버러리들이 이미 개발되어 있다. 반드시 검색을 해보아야만 한다.

효율적인 계산에 관심이 있을 수 있다. 이것도 검색을 통해서 해결할 수 있다.

 

예를 들어, 제법 복잡한 그림을 그린다고 하자. 대부분의 경우들에 대해서 기본틀을 인터넷에서 찾을 수 있다.

파이썬, 특히, 그림 그리기에서 아주 유용하다.  

인터넷 검색창에 아주 구체적으로 자신이 찾고자 하는 일을 넣어주면 된다.

마지막에 컴퓨터 언어 이름을 넣어주면 검색에 많은 도움이 된다.

 

풀고자 하는 문제는 최대한 구체적으로 영어로 작문하여 검색창에 넣어준다.

파이썬으로 어떻게 처리하는지를 직접 물어본다. 에러를 해결하기 위해서도 마찬가지이다.

발생한 에러 메시지를 그대로 검색창에 넣어 준다.

 

많은 경우, 컴퓨터 프로그램 개발에서 그림들을 그려야만 한다. 프로그램 개발과 동시에 그림을 그릴 수 있다.

이 때, matplotlib를 이용해서 중간 중간에 그림을 그린다. 그림으로 프로그램 중간단계의 계산 상태를 점검할 수 있다. 계산을 진행하면서 eps 형식의 그림을 저장할 수 있다. 계속해서 그림들을 만들기 때문에 문제가 발생할 수 있는데, 각각의 그림을 그릴 때 항상 close가 필요하다. 이렇게 하면, 매우 많은 그림들을 순차적으로 다 그릴 수 있다. 파일로 저장할 수 있다.

import matplotlib.pyplot as plt

plt.savefig(str1,dpi=1200)

plt.close()

 

또한, 최종 출력 파일을 pasing 해서 그림을 그릴 수도 있다. 문자열 처리에 익숙해져야만 한다.

파이썬 프로그래밍에서는 알고리듬 구현, 중간 단계 테스트, 그림 그리기, 데이터 뽑아내기, 최종 그림 그리기 등의 작업을 한 꺼번에 수행하는 것이 가능하다. 별도의 그림 그리기 작업을 수행하지 않게 된다. 매 단계마다 출력을 그림으로 저장할 수 있다.

파이썬 프로그램은 느리다. 하지만, 최종적으로 일을 정리하고 논문을 완성함에 있어서는 파이썬 언어를 활용한 작업이 느리지 않을 수도 있다. 

 

최근 딥러닝이 부각되면서 컴퓨터 언어 파이썬은 다시 한 번 상한가를 구가하게 된다. 왜냐하면 대부분의 딥러닝 라이버러리, 패키지들이 파이썬 언어로 되어 있기 때문이다.

 

이공계 분야이고 주요 컴퓨터 언어가 없다면 파이썬 언어를 추천한다.

 

https://en.wikipedia.org/wiki/Python_(programming_language)

파이썬 2.7은 3개월 정도 지나면 은퇴를 할 것이다. 파이썬 3.0 기준으로 살아 가는 것이 마땅하다.

 

몇 가지 문제들을 풀어보고자 한다.

 

선형대수학에서 만나는 문제. 행렬(A), 벡터(x), 벡터(b)로 이루어진 문제.

A x = b

cat z1.py 
import numpy as np
from scipy import linalg
b=[5., -1., 5.]
a=[[-11., 1., -4.], [9., -1., 3.], [4., 3., -1.]]
b=np.array(b)
a=np.array(a)
print( linalg.solve(a,b))

 

b=np.array([5., -1., 5.])
a=np.matrix( [[-11., 1., -4.], [9., -1., 3.], [4., 3., -1.]])
print( linalg.solve(a,b))

 

[ihlee@KRISS-TUCANA temp]$ python z1.py 
[ 2.26666667 -4.2        -8.53333333]
[ 2.26666667 -4.2        -8.53333333]

-------------------------------------------------------------------------------------------------------------------

FFT:

#!/usr/bin/env python
from numpy.fft import fft
from numpy import zeros
aa=((1., 2., 3., 2., 1. ))
aa_fft=fft(aa)/len(aa)
for i in range(len(aa)):
   print aa_fft[i]
bb=zeros(len(aa)*2,complex)
for i in range(len(aa)):
   bb[i]=aa_fft[i]
#for i in range(len(bb)):
#   print bb[i]

bb_fft=fft(bb)
for i in range(len(bb)):
   print bb_fft[i]


(1.8+0j)
(-0.42360679775-0.307768353718j)
(0.02360679775+0.0726542528005j)
(0.02360679775-0.0726542528005j)
(-0.42360679775+0.307768353718j)
(1+0j)
(1.8-5.51642065361e-16j)
(1+0j)
(1.8+0.906153718637j)
(2+0j)
(1.8+0.760845213036j)
(3+0j)
(1.8-0.760845213036j)
(2+0j)
(1.8-0.906153718637j)

----------------------------------------------
#!/usr/bin/env python
from numpy.fft import fft
from numpy import zeros
aa=((0., 1., 2., 3., 2., 1. ))
print '======'
for i in range(len(aa)):
   print aa[i]
print '======'
aa_fft=fft(aa)/len(aa)
for i in range(len(aa)):
   print aa_fft[i]
print '======'
bb=zeros(len(aa)*2,complex)
for i in range(len(aa)*2):
   bb[i]=0.+0.j
for i in range(len(aa)):
   bb[2*i]=aa_fft[i]

print '======'

bb_fft=fft(bb)
for i in range(len(bb)):
   print bb_fft[i]


======
0.0
1.0
2.0
3.0
2.0
1.0
======
(1.5+0j)
(-0.666666666667+1.48029736617e-16j)
0j
(-0.166666666667-2.90348170526e-16j)
0j
(-0.666666666667+1.4231843391e-16j)
======
======
0j
(1+6.38378239159e-16j)
(2-1.13797860024e-15j)
(3+0j)
(2-6.38378239159e-16j)
(1+1.13797860024e-15j)
0j
(1+6.38378239159e-16j)
(2-1.13797860024e-15j)
(3+0j)
(2-6.38378239159e-16j)
(1+1.13797860024e-15j)

-------------------------------------------------------------------------------------------------------------------

 

import numpy as np
from scipy import linalg
b=[5., -1., 5.]
a=[[-11., 1., -4.], [9., -1., 3.], [4., 3., -1.]]
b=np.array(b)
a=np.array(a)
print( linalg.solve(a,b))

 

b=np.array([5., -1., 5.])
a=np.matrix( [[-11., 1., -4.], [9., -1., 3.], [4., 3., -1.]])
print( linalg.solve(a,b))

 

from scipy.misc import derivative
from math import e,pi
print(derivative(np.sin,pi, dx=0.1, order=5))

 

from scipy import integrate
x2= lambda x: x**2
print(integrate.quad(x2,0.,1.))
print(1.**3/3.)

 

 python z1.py 
[ 2.26666667 -4.2        -8.53333333]
[ 2.26666667 -4.2        -8.53333333]
-0.9999966706326076
(0.33333333333333337, 3.700743415417189e-15)
0.3333333333333333

-------------------------------------------------------------------------------------------------------------------

import numpy as np
A = np.array([[3,1,-1], [1,3,-1], [-1,-1,5]])
w,v = np.linalg.eig(A)
print( w)
idx = w.argsort()[::-1]  # large to small
# idx = w.argsort()        # small to large
w = w[idx]
v = v[:,idx]
print (w)

[6. 2. 3.]
[6. 3. 2.]

 

-------------------------------------------------------------------------------------------------------------------

-------------------------------------------------------------------------------------------------------------------

파이썬에서 모듈이란?

 

프로그램을 설치하지 않고도 복사만 해도 되는 경우가 있을 수 있다.

같은 디렉토리에 파이썬 프로그램을 복사해두면 그것이 모듈이되고 불러서 사용할 수 있다.

-------------------------------------------------------------------------------------------------------------------

1차원 슈뢰딩거 방정식 풀이 (파이썬)

1D Schroedinger equation solver (Python)

 

# written by In-Ho Lee, KRISS, 26 Dec. 2001
# one-dimensional Schroedinger solver
# one-dimensional box, harmonic oscillator
# variational FD method; Phys. Rev. B 64, 193101 (2001).
#
import os,sys,time
from shutil import copyfile
from posix import getpid,rename,remove,mkdir,rmdir,system
from time import *
import numpy as np
from scipy import linalg

 

def ev_sort(eigval,eigvec):
 newval=np.zeros(len(eigval),float)                    #
 newvec=np.zeros(eigvec.shape,float)
 index=np.argsort(eigval)
 for i in index:
  newval[i]=eigval[index[i]]
  newvec[i,:]=eigvec[index[i],:]                     #
 return newval,newvec

 

def gen_T_op_vfd1(N,hh):
 print('-----vfd1-is-used-----')
 fac0=-np.pi*np.pi/2.e0                                          #
 fac1=np.pi*np.pi/4.e0
 hh2=hh*hh
 T=np.zeros((N,N),float)
 T=T+np.identity(N)*fac0/hh2*(-0.5)
 for i in range(N-1):
  T[i,i+1]=T[i+1,i]=fac1/hh2*(-0.5)
 return T

 

def gen_T_op_vfd2(N,hh):
 print('-----vfd2-is-used-----')
 fac0=-0.5-3.0*np.pi*np.pi/8.0
 fac1=np.pi*np.pi/4.0
 fac2=0.25-np.pi*np.pi/16.0
 hh2=hh*hh
 T=np.zeros((N,N),float)
 T=T+np.identity(N)*fac0/hh2*(-0.5)
 for i in range(N-1):
  T[i,i+1]=T[i+1,i]=fac1/hh2*(-0.5)
  if i < N-2:
   T[i,i+2]=T[i+2,i]=fac2/hh2*(-0.5)
 return T

 

def gen_T_op_vfd3(N,hh):
 print('-----vfd3-is-used-----')
 fac0=-5./6.-5.0*np.pi*np.pi/16.0
 fac1=1./12.+15.*np.pi*np.pi/64.0
 fac2=5.0/12.0-3.0*np.pi*np.pi/32.0
 fac3=-1.0/12.0+np.pi*np.pi/64.0
 hh2=hh*hh
 T=np.zeros((N,N),float)
 T=T+np.identity(N)*fac0/hh2*(-0.5)
 for i in range(N-1):
  T[i,i+1]=T[i+1,i]=fac1/hh2*(-0.5)
  if i < N-2:
   T[i,i+2]=T[i+2,i]=fac2/hh2*(-0.5)
   if i < N-3:
     T[i,i+3]=T[i+3,i]=fac3/hh2*(-0.5)
 return T

def gen_T_op_vfd4(N,hh):
 print('-----vfd4-is-used-----')
 fac0=-77./72.-35.0*np.pi*np.pi/128.0
 fac1=8./45.+7.*np.pi*np.pi/32.0
 fac2=23.0/45.0-7.0*np.pi*np.pi/64.0
 fac3=-8.0/45.0+np.pi*np.pi/32.0
 fac4=17.0/720.0-np.pi*np.pi/256.0
 hh2=hh*hh
 T=np.zeros((N,N),float)
 T=T+np.identity(N)*fac0/hh2*(-0.5)
 for i in range(N-1):
  T[i,i+1]=T[i+1,i]=fac1/hh2*(-0.5)
  if i < N-2:
   T[i,i+2]=T[i+2,i]=fac2/hh2*(-0.5)
   if i < N-3:
     T[i,i+3]=T[i+3,i]=fac3/hh2*(-0.5)
     if i < N-4:
       T[i,i+4]=T[i+4,i]=fac4/hh2*(-0.5)
 return T

 

def gen_T_op_vfd5(N,hh):
 print('-----vfd5-is-used-----')
 fac0=-449./360.-63.*np.pi*np.pi/256.
 fac1=4./15.+105.*np.pi*np.pi/512.
 fac2=59./105.-15.*np.pi*np.pi/128.
 fac3=-82./315.+45.*np.pi*np.pi/1024.
 fac4=311./5040.-5.*np.pi*np.pi/512.
 fac5=-2./315.+np.pi*np.pi/1024.
 hh2=hh*hh
 T=np.zeros((N,N),float)
 T=T+np.identity(N)*fac0/hh2*(-0.5)
 for i in range(N-1):
  T[i,i+1]=T[i+1,i]=fac1/hh2*(-0.5)
  if i < N-2:
   T[i,i+2]=T[i+2,i]=fac2/hh2*(-0.5)
   if i < N-3:
     T[i,i+3]=T[i+3,i]=fac3/hh2*(-0.5)
     if i < N-4:
       T[i,i+4]=T[i+4,i]=fac4/hh2*(-0.5)
       if i < N-5:
        T[i,i+5]=T[i+5,i]=fac5/hh2*(-0.5)
 return T

 

def gen_T_op_vfd6(N,hh):
 print('-----vfd6-is-used-----')
 fac0=-2497./1800.-231.*np.pi*np.pi/1024.
 fac1=26./75.+99.*np.pi*np.pi/512.
 fac2=493./840.-495.*np.pi*np.pi/4096.
 fac3=-103./315.+55.*np.pi*np.pi/1024.
 fac4=2647./25200.-33.*np.pi*np.pi/2048.
 fac5=-31./1575.+3.*np.pi*np.pi/1024.
 fac6=1./600.-np.pi*np.pi/4096.
 hh2=hh*hh
 T=np.zeros((N,N),float)
 T=T+np.identity(N)*fac0/hh2*(-0.5)
 for i in range(N-1):
  T[i,i+1]=T[i+1,i]=fac1/hh2*(-0.5)
  if i < N-2:
   T[i,i+2]=T[i+2,i]=fac2/hh2*(-0.5)
   if i < N-3:
     T[i,i+3]=T[i+3,i]=fac3/hh2*(-0.5)
     if i < N-4:
       T[i,i+4]=T[i+4,i]=fac4/hh2*(-0.5)
       if i < N-5:
        T[i,i+5]=T[i+5,i]=fac5/hh2*(-0.5)
        if i < N-6:
          T[i,i+6]=T[i+6,i]=fac6/hh2*(-0.5)
 return T

def gen_T_op(N,hh):
 hh2=hh*hh
 T=np.zeros((N,N),float)
 T=T+np.identity(N)/hh2
 for i in range(N-1):
  T[i,i+1]=T[i+1,i]=-0.5/hh2           #
 return T

 

def gen_V_op(N,xs,xf,well_width):    #
 dx=(xf-xs)/float(N-1)
 V=np.zeros((N,N),float)
 for i in range(N):
  xx=xs+(i)*dx
  if abs(xx) <= well_width/2.0e0:
   V[i,i]=0.0e0
  else:
   V[i,i]=800.0e0
 return V

 

def gen_V_op_ho(N,xs,xf,springk):   #
 dx=(xf-xs)/float(N-1)
 V=np.zeros((N,N),float)
 for i in range(N):
  xx=xs+(i)*dx
  V[i,i]=0.5e0*springk*xx*xx            #
 return V

 

def write_xmgr_data(filename,xs,xf,vec):
 dx=(xf-xs)/float(N-1)
 file=open(filename,'w')
 nvec,nel = vec.shape
# print(nvec,nel)
 for j in range(nvec):                         #
  for i in range(nel):
   xx=xs+dx*i
#  file.write('%f '  % xx)
#  file.write('%f \n' % vec[j][i])
   file.write('%f %f \n'  %( xx, vec[j][i] ))
  file.write('& \n')
  file.write('# \n')
 file.close()
 return

 

 

N=4000
#N=1500
N=1000
N=500
#N=100
#N=50
well_width=10.e0   # 10.d0 is not allowed in python
border_width=5.e0
springk=1.e0       # spring constant for harmonic oscillator potential
xs=-abs(well_width+border_width)/2.0e0
xf=-xs
hh=(xf-xs)/float(N-1)
print('----------------------------------------------')
print('1D Schroedinger solver with Python')
print('current working directory :',os.getcwd())
print('processor id :',os.getpid())
print('')
t_start=clock()
secs=time()
print('     time :',secs)
print('   gmtime :',gmtime(secs))
print('localtime :',localtime(secs))
tuple_input=localtime(secs)
tuple_output=asctime(tuple_input)
print('timetuple :',tuple_output)

print('[---------------------------------------------')
print('the number of grid points :', N)
print('xs,xf in a.u.', xs,xf)
print('dx in a.u.', hh)
print('---------------------------------------------]')

 

path_to_xmgr='xmgr'

 

#T=gen_T_op(N,hh)
#T=gen_T_op_vfd1(N,hh)
#T=gen_T_op_vfd2(N,hh)
#T=gen_T_op_vfd3(N,hh)
#T=gen_T_op_vfd4(N,hh)
#T=gen_T_op_vfd5(N,hh)
T=gen_T_op_vfd6(N,hh)

 

#V=gen_V_op(N,xs,xf,well_width)
V=gen_V_op_ho(N,xs,xf,springk)

 

H=T+V                     # Hamiltonian setup
val,vec=np.linalg.eig(H)
val,vec=ev_sort(val,vec)  # solution sorting

 

print('[---------------------------------------------')
print("6 lowest eigenvalues in a.u.")
print(val[:6])
#print val[1]/val[0]
#print val[2]/val[0]
#print val[3]/val[0]
#print val[4]/val[0]
#print val[5]/val[0]
print('---------------------------------------------]')
dat_filename='tmp1.dat'
write_xmgr_data(dat_filename,xs,xf,vec[:6]) # data file generation for Xmgr

 

t_final=clock()
print('CPU-TIME (sec):',t_final-t_start)

 

----------------------------------------------
1D Schroedinger solver with Python
current working directory : /home/ihlee/temp
processor id : 77148

 

     time : 1568443975.8922265
   gmtime : time.struct_time(tm_year=2019, tm_mon=9, tm_mday=14, tm_hour=6, tm_min=52, tm_sec=55, tm_wday=5, tm_yday=257, tm_isdst=0)
localtime : time.struct_time(tm_year=2019, tm_mon=9, tm_mday=14, tm_hour=15, tm_min=52, tm_sec=55, tm_wday=5, tm_yday=257, tm_isdst=0)
timetuple : Sat Sep 14 15:52:55 2019
[---------------------------------------------
the number of grid points : 500
xs,xf in a.u. -7.5 7.5
dx in a.u. 0.03006012024048096
---------------------------------------------]
-----vfd6-is-used-----
[---------------------------------------------
6 lowest eigenvalues in a.u.
[0.5 1.5 2.5 3.5 4.5 5.5]
---------------------------------------------]
CPU-TIME (sec): 0.28

-------------------------------------------------------------------------------------------------------------------

FFT with python

 

-------------------------------------------------------------------------------------------------------------------

winding number calculation (topological invariant) 위상 불변량 계산

 

-------------------------------------------------------------------------------------------------------------------

Tex in matplotlib [python]  그림에서 Tex 사용하기

 

-------------------------------------------------------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt
pi = np.pi

 

npoints = 360
r = 1.          # AU
dt = 1./npoints # fractions of a year
mu = 4. * pi**2  # 
x = np.zeros(npoints)
y = np.zeros(npoints)
vx = np.zeros(npoints)
vy = np.zeros(npoints)

 

# Initial Conditions
x[0] = r               # (x0 = r, y0 = 0) AU
vy[0] = np.sqrt(mu/r)  # (vx = 0, vy = sqrt(mu/r)) AU/yr

 

for i in range(0,npoints-1):
    vx[i+1]=vx[i]-4.*pi**2*x[i]/(r**3)*dt
    vy[i+1]=vy[i]-4.*pi**2*y[i]/(r**3)*dt
    x[i+1]=x[i]+vx[i+1]*dt
    y[i+1]=y[i]+vy[i+1]*dt

 

plt.plot(x, y, 'b')
plt.axis('equal')
plt.show()

 

-------------------------------------------------------------------------------------------------------------------

 

import numpy as np
import matplotlib.pyplot as plt
import os, sys
import matplotlib

matplotlib.rc('font', size=18)
matplotlib.rc('font', family='Arial')

# definition of numerical parameters
N = 51        # number of grid points
dt = 5.e-4   # time step
L = float(1)    # size of grid
nsteps = 620   # number of time steps
dx = L/(N-1)    # grid spacing
nplot = 20         # number of timesteps before plotting

r = dt/dx**2      # assuming heat diffusion coefficient == 1

# initialize matrices A, B and b array
A = np.zeros((N-2,N-2))
B = np.zeros((N-2,N-2))
b = np.zeros((N-2))
# define matrices A, B and b array
for i in range(N-2):
    if i==0:
        A[i,:] = [2+2*r if j==0 else (-r) if j==1 else 0 for j in range(N-2)]
        B[i,:] = [2-2*r if j==0 else r if j==1 else 0 for j in range(N-2)]
        b[i] = 0. #boundary condition at i=1
    elif i==N-3:
        A[i,:] = [-r if j==N-4 else 2+2*r if j==N-3 else 0 for j in range(N-2)]
        B[i,:] = [r if j==N-4 else 2-2*r if j==N-3 else 0 for j in range(N-2)]
        b[i] = 0. #boundary condition at i=N
    else:
        A[i,:] = [-r if j==i-1 or j==i+1 else 2+2*r if j==i else 0 for j in range(N-2)]
        B[i,:] = [r if j==i-1 or j==i+1 else 2-2*r if j==i else 0 for j in range(N-2)]

# initialize grid
x = np.linspace(0,1,N)
# initial condition
u = np.asarray([2*xx if xx<=0.5 else 2*(1-xx) for xx in x])
# evaluate right hand side at t=0
bb = B.dot(u[1:-1]) + b

fig = plt.figure()
plt.plot(x,u,linewidth=2)
filename = 'foo000.jpg'
#fig.set_tight_layout(True,"h_pad=1.0");
plt.tight_layout(pad=3.0)
plt.xlabel("x")
plt.ylabel("u")
plt.title("t = 0")
plt.savefig(filename,format="jpg")
plt.clf()

c = 0
for j in range(nsteps):
    print(j)
    # find solution inside domain
    u[1:-1] = np.linalg.solve(A,bb)
    # update right hand side
    bb = B.dot(u[1:-1]) + b
    if(j%nplot==0): #plot results every nplot timesteps
        plt.plot(x,u,linewidth=2)
        plt.ylim([0,1])
        filename = 'foo' + str(c+1).zfill(3) + '.jpg';
        plt.xlabel("x")
        plt.ylabel("u")
        plt.title("t = %2.2f"%(dt*(j+1)))
        plt.savefig(filename,format="jpg")
        plt.clf()
        c += 1

os.system("ffmpeg -y -i 'foo%03d.jpg' heat_equation.m4v")
os.system("rm -f *.jpg")

-------------------------------------------------------------------------------------------------------------------

x=np.array([2, 1, 4, 3, 5])

np.sort(x)

x.sort()

x= np.array([2, 1, 4, 3, 5])

i=np.argsort(x)

x[i]

-------------------------------------------------------------------------------------------------------------------

 import numpy as np
 temperature = np.random.random(1024)

 temperature
array([ 0.24249738, 0.6407523 , 0.84243584, ..., 0.90018119,
0.64844851, 0.65660748])

 

-------------------------------------------------------------------------------------------------------------------

# Import a plotting libraries and a maths library 
import matplotlib.pyplot as plt
import numpy as np

 

r = np.linspace(0.01,3.0,num=500) # Make a radius vector
epsilon = 1     # Energy minimum
sigma = 1       # Distance to zero crossing point
E_LJ = 4*epsilon*((sigma/r)**12-(sigma/r)**6) # Lennard-Jones potential

 

plt.figure(figsize=[6,6])
plt.plot(r,E_LJ,'r-',linewidth=1,label=r" $LJ\; pot$") # Red line is unshifted LJ

 

# The cutoff and shifting value
Rcutoff = 2.5
phicutoff = 4.0/(Rcutoff**12)-4.0/(Rcutoff**6) # Shifts the potential so at the cutoff the potential goes to zero

 

E_LJ_shift = E_LJ - phicutoff # Subtract the value of the potential at r=2.5

 

plt.plot(r[:415],E_LJ_shift[:415],'b-',linewidth=1,label=r"$LJ\; pot\; shifted$") # Blue line is shifted

 

# Plot formatting
plt.rc('text', usetex=True)
plt.rc('xtick', labelsize=20) 
plt.rc('ytick', labelsize=20) 
plt.title(r"$Lennard-Jones\; potential$",fontsize=20)
plt.xlim([0.0,3.0])
plt.ylim([-1.5,1.5])
plt.ylabel(r"$E_{LJ}/\epsilon$",fontsize=20)
plt.xlabel(r"$r/\sigma$",fontsize=20)
plt.legend(frameon=False,fontsize=20)
plt.axhline(0, color='grey',linestyle='--',linewidth=2)
plt.axvline(1, color='grey',linestyle='--',linewidth=2)

plt.show()

import numpy as np
def Compute_Forces(pos,acc,ene_pot,epsilon,BoxSize,DIM,N):
    # Compute forces on positions using the Lennard-Jones potential
    # Uses double nested loop which is slow O(N^2) time unsuitable for large systems
    Sij = np.zeros(DIM) # Box scaled units
    Rij = np.zeros(DIM) # Real space units

 

    # Set all variables to zero
    ene_pot = ene_pot*0.0
    acc = acc*0.0
    virial=0.0

 

    # Loop over all pairs of particles
    for i in range(N-1):
        for j in range(i+1,N): #i+1 to N ensures we do not double count
            Sij = pos[i,:]-pos[j,:] # Distance in box scaled units
            for l in range(DIM): # Periodic interactions
                if (np.abs(Sij[l])>0.5):
                    Sij[l] = Sij[l] - np.copysign(1.0,Sij[l]) # If distance is greater than 0.5  (scaled units) then subtract 0.5 to find periodic interaction distance.

 

            Rij = BoxSize*Sij # Scale the box to the real units in this case reduced LJ units
            Rsqij = np.dot(Rij,Rij) # Calculate the square of the distance

 

            if(Rsqij < Rcutoff**2):
                # Calculate LJ potential inside cutoff
                # We calculate parts of the LJ potential at a time to improve the efficieny of the computation (most important for compiled code)
                rm2 = 1.0/Rsqij # 1/r^2
                rm6 = rm2**3.0 # 1/r^6
                rm12 = rm6**2.0 # 1/r^12
                phi = epsilon*(4.0*(rm12-rm6)-phicutoff) # 4[1/r^12 - 1/r^6] - phi(Rc) - we are using the shifted LJ potential
                # The following is dphi = -(1/r)(dV/dr)
                dphi = epsilon*24.0*rm2*(2.0*rm12-rm6) # 24[2/r^14 - 1/r^8]
                ene_pot[i] = ene_pot[i]+0.5*phi # Accumulate energy
                ene_pot[j] = ene_pot[j]+0.5*phi # Accumulate energy
                virial = virial - dphi*Rsqij # Virial is needed to calculate the pressure
                acc[i,:] = acc[i,:]+dphi*Sij # Accumulate forces
                acc[j,:] = acc[j,:]-dphi*Sij # (Fji=-Fij)
    return acc, np.sum(ene_pot)/N, -virial/DIM # return the acceleration vector, potential energy and virial coefficient

 

 

def Calculate_Temperature(vel,BoxSize,DIM,N):

 

    ene_kin = 0.0

 

    for i in range(N):
        real_vel = BoxSize*vel[i,:]
        ene_kin = ene_kin + 0.5*np.dot(real_vel,real_vel)

 

    ene_kin_aver = 1.0*ene_kin/N
    temperature = 2.0*ene_kin_aver/DIM

 

    return ene_kin_aver,temperature

 

# Setting up the simulation
NSteps=10000 # Number of steps
deltat = 0.0032 # Time step in reduced time units
TRequested = 0.5# #Reduced temperature
DumpFreq = 100 # Save the position to file every DumpFreq steps
epsilon = 1.0 # LJ parameter for the energy between particles

 

# Main MD loop
def main(pos,NSteps,deltat,TRequested,DumpFreq,epsilon,BoxSize,DIM):

 

    # Vectors to store parameter values at each step
    N = np.size(pos[:,1])
    ene_kin_aver = np.ones(NSteps)
    ene_pot_aver = np.ones(NSteps)
    temperature = np.ones(NSteps)
    virial = np.ones(NSteps)
    pressure = np.ones(NSteps)
    ene_pot = np.ones(N)

 

    vel = (np.random.randn(N,DIM)-0.5)
    acc = (np.random.randn(N,DIM)-0.5)

 

    # Open file which we will save the outputs to
    f = open('traj.xyz', 'w')

 

    for k in range(NSteps):

 

        # Refold positions according to periodic boundary conditions
        for i in range(DIM):
            period = np.where(pos[:,i] > 0.5)
            pos[period,i]=pos[period,i]-1.0
            period = np.where(pos[:,i] < -0.5)
            pos[period,i]=pos[period,i]+1.0

 

        # r(t+dt) modify positions according to velocity and acceleration
        pos = pos + deltat*vel + 0.5*(deltat**2.0)*acc # Step 1

 

        # Calculate temperature
        ene_kin_aver[k],temperature[k] = Calculate_Temperature(vel,BoxSize,DIM,N)

 

        # Rescale velocities and take half step
        chi = np.sqrt(TRequested/temperature[k])
        vel = chi*vel + 0.5*deltat*acc # v(t+dt/2) Step 2

 

        # Compute forces a(t+dt),ene_pot,virial
        acc, ene_pot_aver[k], virial[k] = Compute_Forces(pos,acc,ene_pot,epsilon,BoxSize,DIM,N) # Step 3

 

        # Complete the velocity step 
        vel = vel + 0.5*deltat*acc # v(t+dt/2) Step 4

 

        # Calculate temperature
        ene_kin_aver[k],temperature[k] = Calculate_Temperature(vel,BoxSize,DIM,N)

 

        # Calculate pressure
        pressure[k]= density*temperature[k] + virial[k]/volume

 

        # Print output to file every DumpFreq number of steps
        if(k%DumpFreq==0): # The % symbol is the modulus so if the Step is a whole multiple of DumpFreq then print the values

 

            f.write("%s\n" %(N)) # Write the number of particles to file
            # Write all of the quantities at this step to the file
            f.write("Energy %s, Temperature %.5f\n" %(ene_kin_aver[k]+ene_pot_aver[k],temperature[k]))
            for n in range(N): # Write the positions to file
                f.write("X"+" ")
                for l in range(DIM):
                    f.write(str(pos[n][l]*BoxSize)+" ")
                f.write("\n")

 

            if(DIM==2):
                import matplotlib.pyplot as plt
                from IPython import display
                plt.cla()
                plt.xlim(-0.5*BoxSize,0.5*BoxSize)
                plt.ylim(-0.5*BoxSize,0.5*BoxSize)
                for i in range(N):
                    plt.plot(pos[i,0]*BoxSize,pos[i,1]*BoxSize,'o',markersize=20,)
                display.clear_output(wait=True)
                display.display(plt.gcf())
        #print(ene_kin_aver[k], ene_pot_aver[k], temperature[k], pressure[k])

 

    f.close() # Close the file

 

    return ene_kin_aver, ene_pot_aver, temperature, pressure, pos

 

--------------------------------------------------------------------------------------------------------------------

import numpy as np
from scipy.stats import pearsonr
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator, FormatStrFormatter, AutoMinorLocator
def gen_cmp1(ngrid,wrkrepl,wrkrcry,dist0,ireplica,kkdd,ispgrp0):
    plt.figure(figsize=(6,6))
    ax=plt.axes()
#   ax.set_xlabel(r'$r$'+' (\u00c5)',fontsize=22)
    ax.set_xlabel('Bond length'+' (\u00c5)',fontsize=22)
    ax.set_ylabel('Frequency (arb. unit)',fontsize=22)
    majorLocator= MultipleLocator(1)
    minorLocator= AutoMinorLocator()
    majorFormatter= FormatStrFormatter('%d')
    minorFormatter= FormatStrFormatter('%d')
    ax.xaxis.set_major_locator(majorLocator)
    ax.xaxis.set_major_formatter(majorFormatter)
    ax.xaxis.set_minor_locator(minorLocator)
    majorLocator= MultipleLocator(0.1)
    minorLocator= AutoMinorLocator()
    majorFormatter= FormatStrFormatter('%4.2f')
    minorFormatter= FormatStrFormatter('%4.2f')
    ax.yaxis.set_major_locator(majorLocator)
    ax.yaxis.set_major_formatter(majorFormatter)
    ax.yaxis.set_minor_locator(minorLocator)
    ax.tick_params(which='major', length=2, color='black')
    ax.tick_params(which='minor', length=4, color='brown')
#  

#   ax.set_facecolor("cornsilk")
#   ax.set_facecolor("ivory")
#   ax.set_facecolor("PeachPuff")
#   ax.set_facecolor("PapayaWhip")
#   ax.set_facecolor("PowderBlue")
#   ax.set_facecolor("Azure")
#   ax.set_facecolor("Snow")
#   ax.set_facecolor("Honeydew")
#   ax.set_facecolor("GhostWhite")
#   ax.set_facecolor("AliceBlue")
#   ax.set_facecolor("MintCream")
#   ax.set_facecolor("WhiteSmoke")
#   ax.set_facecolor("Seashell")
#   ax.set_facecolor("FloralWhite")
#   ax.set_facecolor("LavenderBlush")
#   ax.set_facecolor("Gainsboro")
#   ax.set_facecolor("LightGray")
#   ax.set_facecolor("LightPink")
#   ax.set_facecolor("LemonChiffon")
#   ax.set_facecolor("LightGoldenrodYellow")
#   ax.set_facecolor("LightGreen")
#   ax.set_facecolor("PaleGreen")
#   ax.set_facecolor("BlanchedAlmond")
    ax.set_facecolor("OldLace")
#
    ax.set_facecolor("ivory")
    plt.text(1.3, 0.45, "Pearson's distance = "+f'{dist0: .3f}', {'color' : 'blue', 'fontsize' : 12} )
    plt.text(1.3, 0.43, "space group number = "+str(ispgrp0), {'color' : 'red', 'fontsize' : 12} )
    plt.grid(True)
    plt.xlim(1., 7.)
    plt.ylim(0., 0.5)
    xgrd=np.linspace(0.0, 12.0, ngrid)
#   plt.plot(xgrd[:134],wrkrepl[:134], '-', lw=5, alpha=0.4)
    plt.plot(xgrd[:134],wrkrepl[:134], 'o', lw=5, alpha=0.4)
#   plt.plot(xgrd[:134],wrkrcry[:134], '--o', lw=2)
    plt.plot(xgrd[:134],wrkrcry[:134], '--b', lw=2)
    plt.tight_layout()
    str1='cmp'+str(ireplica)+'_'+str(kkdd)+'.eps'
    plt.savefig(str1,dpi=1200)
    plt.show()
#   plt.close()

 

ngrid=200
ireplica=1
kkdd=1
wrkrepl=[]
wrkrcry=[]
for i in range(ngrid):
    i1=0.+0.3*np.sin(float(i)*(12.-0.)/float(ngrid-1))
    wrkrepl.append(i1)
    i1=0.+0.3*np.cos(float(i)*(12.-0.)/float(ngrid-1))
    wrkrcry.append(i1)

 

wrkrepl=np.array(wrkrepl)
wrkrcry=np.array(wrkrcry)
corr, _ = pearsonr(wrkrepl,wrkrcry)
dist0=1.-corr
ispgrp0=230
gen_cmp1(ngrid,wrkrepl,wrkrcry,dist0,ireplica,kkdd,ispgrp0)
 

-------------------------------------------------------------------------------------------------------------------

polynomial fitting

 

>>> import numpy as np
>>> x=np.zeros(10)
>>> y=np.zeros(10)
>>> for i in range(10):
...    x[i]=i*1.
... 
>>> y[1]=1.1
>>> y[0]=1.1
>>> y[1]=3.9
>>> y[2]=11.2
>>> y[3]=21.5
>>> y[4]=34.8
>>> y[5]=51.
>>> y[6]=70.2
>>> y[7]=92.3
>>> y[8]=117.4
>>> y[9]=145.5
>>> p=np.poly1d(np.polyfit(x,y,deg=2))
>>> print(p)
       2
1.517 x + 2.483 x + 0.4927

-------------------------------------------------------------------------------------------------------------------

def angle(a,b):
    """ calculate the angle between vector a and b """
    return acos(np.dot(a,b)/np.linalg.norm(a)/np.linalg.norm(b))

 

 

print(res.x)

print(res.fun)

 

 

-------------------------------------------------------------------------------------------------------------------

#      calculate the Pearson's correlation between two variables
from numpy.random import randn
from numpy.random import seed
from scipy.stats import pearsonr
#      seed random number generator
seed(1)
#      prepare data
data1 = 20 * randn(1000) + 100
data2 = data1 + (10 * randn(1000) + 50)
#      calculate Pearson's correlation
corr, _ = pearsonr(data1, data2)
print('Pearsons correlation: %.3f' % corr)

 

-------------------------------------------------------------------------------------------------------------------

 

-------------------------------------------------------------------------------------------------------------------

 

-------------------------------------------------------------------------------------------------------------------

[출처] http://egloos.zum.com/incredible/v/7467457