added new directory for lecture notes

gp.jl

@@ -4,6 +4,7 @@ export gpr, CovG, CovL, CovGsin, CovLsin, CovDot, CovDotWeight
 
 CovG = ( (x1, x2, sig::Float64) -> exp(- sumabs2(x1-x2)/(2.0*sig^2)))
 CovGsin = ( (x1, x2, sig::Float64) -> exp(- 2.0*sumabs2(sin((x1-x2)/2.0))/sig^2))
+CovGsin1 = ( (x1, x2, sig::Float64) -> exp(- 2.0*sumabs2(sin((x1[1]-x2[1])/2.0))/(4.0*sig)^2)*exp(- sumabs2(x1[2:end]-x2[2:end])/(2.0*sig^2)))
 CovL = ( (x1, x2, sig::Float64) -> exp(- sumabs(x1-x2)/sig))
 CovLsin = ( (x1, x2, sig::Float64) -> exp(- sumabs(sin((x1-x2)/2.0))/sig))
 CovDot = ( (x1, x2, sig::Float64) -> (x1 ⋅ x2)^sig)

lj.jl

@@ -3,7 +3,7 @@ module lj
 import pbc;
 
 global conf = 0.0;    # strength of confining potential 
-global cutoff = 3.0;  # cutoff of the LJ potential for the periodic case, 0 means no cutoff
+global cutoff = 0.0;  # cutoff of the LJ potential for the periodic case, 0 means no cutoff
 global cutoff2 = cutoff^2;
 
 global Vljshift = 0.0;
@@ -55,6 +55,28 @@ function energy(x)
     return e+(rand()-0.5)*1e-12;
 end
 
+#
+#
+# force on particles
+# x is (N,dim) array of positions
+#
+function force(x)
+    N = size(x,2)
+
+    f = zeros(x)
+    for i = 1:N
+        for j=i+1:N
+            dr = x[:,j]-x[:,i]
+            r = norm(dr)
+            df = Vljder(r)*(dr/r)
+            f[:,i] += df
+            f[:,j] -= df
+        end
+    end
+
+    return f
+end
+
 #
 # energy difference of LJ system upon moving particle i by vector dx
 #
@@ -90,6 +112,19 @@ function Vlj(r)
     return (1/r)^12-(1/r)^6-Vljshift;
 end
 
+#
+# derivative of LJ potential
+#
+
+function Vljder(r)
+    if cutoff >0.0 && r > cutoff
+        return 0.0
+    end
+    return -12.0*(1/r)^13+6.0*(1/r)^7
+end
+
+
+
 #
 # LJ pair potential operating on r^2
 #

md.jl

@@ -0,0 +1,28 @@
+module md
+global f_old
+
+export md_step
+
+
+
+function md_langevin_step_baoab(q, p, f, m, dt, kT, gamma)
+
+    global f_old
+    #baoab discretisation of Langevin dynamics
+
+    if !isdefined(f_old)
+        f_old = f
+    end
+
+    p = p + (dt/2)*f_old
+    q = q + (dt/2)*p./m
+    p = p * exp(-dt * gamma) + sqrt(kT*(1-exp(-2*dt*gamma))*m).*randn(size(p))
+    q = q + (dt/2)*p./m
+    p = p + (dt/2)*f
+    f_old = f
+end
+
+#alias
+md_step = md_langevin_step_baoab
+
+end