From f87d3a371c6019a89dbe63753bbdb44193dc16d2 Mon Sep 17 00:00:00 2001
From: arthurfaria2 <arthur.faria2@gmail.com>
Date: Sat, 8 Jan 2022 19:11:07 +0100
Subject: [PATCH] files

---
 Readme.txt       |  11 +++
 stat_langevin.py | 216 +++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 227 insertions(+)
 create mode 100644 Readme.txt
 create mode 100644 stat_langevin.py

diff --git a/Readme.txt b/Readme.txt
new file mode 100644
index 0000000..e36aab5
--- /dev/null
+++ b/Readme.txt
@@ -0,0 +1,11 @@
+Brief files description:
+
+1. stat_langevin.py
+ 	- stochastic dynamics for a Brownian particle. Mean value and variance are calculated. Algorithm based on the publication: "https://aip.scitation.org/doi/abs/10.1063/1.4802990"
+	
+2. FT_langevin.py
+	- stochastic dynamics for a Brownian particle under a harmonic potential whose center of mass is displaced with constant velocity (see work: "https://journals.aps.org/pre/abstract/10.1103/PhysRevE.67.046102"). Fluctuation Theorem for work and heat are calculated. Algorithm based on the publication: "https://aip.scitation.org/doi/abs/10.1063/1.4802990"
+	
+3. Generalized_FP.py
+	- generalized Fokker-Planck (GenBM) equation (see .tex file for further infos). Both the generalized semiclassical distriubtion (GenBM_rho) and the distribution of a standard Brownian motion (BM_rho) are obatined considering a external harmonic potential. Algorithm based on finite diference approach to compute derivatives.
+
diff --git a/stat_langevin.py b/stat_langevin.py
new file mode 100644
index 0000000..fa3801a
--- /dev/null
+++ b/stat_langevin.py
@@ -0,0 +1,216 @@
+
+# coding: utf-8
+
+# Source
+# 
+# https://aip.scitation.org/doi/abs/10.1063/1.4802990
+# 
+# http://hockygroup.hosting.nyu.edu/exercise/langevin-dynamics.html
+Parameters to be used
+
+ks = 2
+
+#gamma**2 >> 4*k overdamped
+gamma = 3 #10 
+kBT = 1.0
+
+dt = 0.001
+tMax = 50
+
+# Sample conditions
+
+N = 10**2 #int(argv[1])
+x_init = 0.3
+v_init = 0.
+# In[1]:
+
+
+import pylab as plt
+import numpy as np
+import time as tm 
+from sys import argv
+
+
+# In[2]:
+
+
+# Harmonic force
+
+def force(x, *args):
+    x0 = 0.
+    ks = args[0]
+    f = -ks*(x)
+    return f
+
+
+# Return the time lenght
+
+def Time_len(tMax, dt):
+
+    t = 0
+    L = []
+    
+    while(t<tMax):
+        L.append(t)
+        t += dt
+
+    return len(L)
+
+
+# In[3]:
+
+
+# Stochastic evolution
+# for further explanation, read: https://aip.scitation.org/doi/abs/10.1063/1.4802990
+
+def position_update(x,v,dt):
+    x_dt = x + v*dt/2.
+    return x_dt
+
+def velocity_update(v,F,dt):
+    v_dt = v + F*dt/2.
+    return v_dt
+
+def random_velocity_update(v,gamma,kBT,dt):
+    R = np.random.normal()
+    c1 = np.exp(-gamma*dt)
+    c2 = np.sqrt(1-np.exp(-2*gamma*dt))*np.sqrt(kBT)
+    v_dt = c1*v + c2*R
+    return v_dt
+
+
+def BAOAB_method(x_init, v_init, tMax, dt, gamma, kBT, ks):
+    
+    x = x_init
+    v = v_init
+    t = 0
+    pos = []
+    time = []
+    
+    while(t<tMax):
+        
+        # part B
+        F = force(x, ks)
+        v = velocity_update(v,F,dt)
+        
+        # part a
+        x = position_update(x,v,dt)
+        
+        # part O
+        v = random_velocity_update(v,gamma,kBT,dt)
+        
+        # part A
+        x = position_update(x,v,dt)
+    
+        # part B
+        F = force(x, ks)
+        v = velocity_update(v,F,dt)
+        
+        pos.append(x)
+        time.append(t)
+        
+        t += dt
+    
+    return time, pos 
+
+
+# In[4]:
+
+
+###################### MAIN #############################
+
+# Parameters and time
+
+ks = 2
+gamma = 1 #gamma**2 >> 4*k overdamped
+kBT = 1.0
+
+dt = 0.01
+tMax = 50
+
+# Sample conditions
+
+N = 10**4 #int(argv[1])
+x_init = 0.3
+v_init = 0.
+
+####################################
+
+# Extras
+
+M_x = np.zeros(Time_len(tMax, dt))
+M_x2 = np.zeros(Time_len(tMax, dt))
+
+#########################################################
+
+start = tm.time()
+
+# Stochastic evolution for each in the sample
+
+for ii in range(N):    
+
+    time, position = BAOAB_method(x_init, v_init, tMax, dt, gamma, kBT, ks)
+        
+    M_x = M_x + np.array(position)
+    M_x2 = M_x2 + np.power(position,2)
+
+t = np.array(time)
+
+
+### Statistics
+
+# Mean
+M_x = M_x/N
+M_x2 = M_x2/N
+
+# Variance
+Var_x = M_x2 - M_x**2
+
+#########################################################
+
+end = tm.time()
+print(end-start)
+
+
+# In[5]:
+
+
+############ Printing data ################
+
+#output = np.array([M_x2, M_x, t])
+
+#data_path = "/home/fariaart/Dropbox/data_%s.txt" %N
+#data_path = "/home/des01/mbonanca/fariaart/Resultados/Doutorado/Lutz/over_data_%s.txt" %N
+#with open(data_path , "w+") as data:
+    #np.savetxt(data, output.T, fmt='%f')
+
+
+# In[6]:
+
+
+############ Plotting data ################
+
+# Variance plot
+
+plt.plot(t, Var_x, label= 'num' )
+
+plt.ylabel(r' $\left<x^2\right> - \left<x\right>^2$', fontsize = 12)
+plt.xlabel(r' $t$', fontsize = 12)
+
+plt.legend(loc='lower center')
+#plt.savefig('/home/fariaart/Dropbox/Pesquisa/Doutorado/Lutz/figuras/var_mean_over_kk/HO_brown/under_2_1.png', transparent=False)
+plt.figure()
+
+
+# Mean value plot
+
+plt.plot(t, M_x, label= 'num' )
+
+plt.ylabel(r' $\left<x\right> $', fontsize = 12)
+plt.xlabel(r' $t$', fontsize = 12)
+#plt.ylim(-0.07,0.07)
+
+plt.legend(loc='upper center')
+#plt.savefig('/home/fariaart/Dropbox/Pesquisa/Doutorado/Lutz/figuras/var_mean_over_kk/HO_brown/under_2_2.png', transparent=False)
+plt.figure()
+