Molecular Dynamics code and Latex

report/MolDynBib.bib

@@ -0,0 +1,50 @@
+@article{deGroot,
+  title = {Computational Biomolecular Dynamics Group},
+  author = {B. de Groot},
+  publisher = {Max Planck Institute for Biophysical Chemistery},
+  url = {http://www3.mpibpc.mpg.de/groups/de_groot/}
+}
+@book{thijssen,
+  author    = {J.M. Thijssen}, 
+  title     = {Computational Physics},
+  publisher = {Cambridge University Press},
+  year      = 2007,
+  volume    = 4,
+  series    = 10,
+  address   = {Cambridge},
+  edition   = 3,
+  month     = 7,
+  isbn      = {3257227892}
+}
+@article{verlet,
+  title = {Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules},
+  author = {Verlet, Loup},
+  journal = {Phys. Rev.},
+  volume = {159},
+  issue = {1},
+  pages = {98--103},
+  numpages = {0},
+  year = {1967},
+  month = {Jul},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRev.159.98},
+  url = {http://link.aps.org/doi/10.1103/PhysRev.159.98}
+}
+@article{crystal,
+  title= {Crystal Structures},
+  author = {Academic Resource Center},
+  publisher = {Illinois Institute for Technology},
+  url = {https://iit.edu/arc/workshops/pdfs/Crystal_Structures.pdf}
+}
+@article{Birkhoff,
+author = {Birkhoff, George D.}, 
+title = {Proof of the Ergodic Theorem},
+volume = {17}, 
+number = {12}, 
+pages = {656-660}, 
+year = {1931}, 
+doi = {10.1073/pnas.17.2.656}, 
+URL = {http://www.pnas.org/content/17/12/656.short}, 
+eprint = {http://www.pnas.org/content/17/12/656.full.pdf+html}, 
+journal = {Proceedings of the National Academy of Sciences} 
+}

src/correlationmodule.py

@@ -0,0 +1,47 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import pylab
+
+def cor(npdim,N,distances,nbins,finalbins,plotflag,Ttarg,density):
+#find the correlation function of the system
+  dmax = np.sqrt(5.)*npdim
+  bin_vec = np.zeros(nbins, dtype = float)
+
+  for i in range(N):
+    for j in range(i+1,N):
+      bin_num = int(distances[i][j]*nbins/dmax)
+      if bin_num < nbins/2.:
+        bin_vec[bin_num] = bin_vec[bin_num] + 1
+
+
+
+
+#normalize based on the volume of the radial shell
+#keeping order R^2*dR -> V=4*PI*R^2*dR
+  dR = dmax/nbins
+  Rin = np.zeros((nbins), dtype=float)
+  for bin_num in range(nbins):
+    Rout = (bin_num + 2)*dR
+    Rin[bin_num] = (bin_num + 1)*dR
+    bin_vec[bin_num] = bin_vec[bin_num]/(4*np.pi*Rin[bin_num]**2*dR)
+
+
+  finalbins[0] = 0.0001 # to make sure that the plot starts at r=0  
+  if plotflag == 1:
+    plt.plot(Rin,finalbins)#,width=dR)
+    plt.xlim([0,max(Rin)*(5/8.)])
+    plt.ylabel('g(r)')
+    plt.xlabel('r')
+    plt.title('Correlationfunction')
+    plt.text(max(Rin)*(1/4.),max(finalbins)-0.1,r'$T$=%s, $\rho$=%s'%(Ttarg,density)) # x and y values for position of text are choosen for 864 particles
+    plt.show()
+    print max(Rin)
+    plt.savefig('correlationfunctionT%sRho%s.jpg'%(Ttarg,density))
+    np.savetxt("cordataT%srho%sN%sRin.txt"%(Ttarg,density,N),Rin)
+    np.savetxt("cordataT%srho%sN%sfinbin.txt"%(Ttarg,density,N),finalbins)
+
+
+
+  return bin_vec
+
+

src/dataplotter.py

@@ -0,0 +1,95 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Importing the data
+
+Rin = np.loadtxt("cordataT0.5rho1.2N864Rin.txt")
+finalbins = np.loadtxt("cordataT0.5rho1.2N864finbin.txt")
+prestime = np.loadtxt("presdataT3rho0.3N864n100lpnum10000prestime.txt")
+averagedp = np.loadtxt("presdataT3rho0.3N864n100lpnum10000averagedp.txt")
+kenarray = np.loadtxt("kenT3rho0.5N864n100lpnum10000.txt")
+potarray = np.loadtxt("potT1rho0.88N864n100lpnum10000.txt")
+totenarray = np.loadtxt("totenT1rho0.88N864n100lpnum10000.txt")
+
+
+# Values for the different parameters, need to be equal to those mentioned in the filenames.
+
+Ttarg = 1
+density = 0.88
+n = 100 #width of pressure block
+N = 864 #number of particles
+lpnum = 10000
+cutoff = int(lpnum/5.)
+
+
+# Plotmodules
+
+def corplotter(Rin,finalbins,Ttarg,density):
+  plt.plot(Rin,finalbins)#,width=dR)
+  plt.xlim([0,6])
+  plt.ylabel('g(r)')
+  plt.xlabel('r')
+  plt.title('Correlationfunction')
+  plt.text(4.2,max(finalbins)-0.1,r'$T$=%s, $\rho$=%s'%(Ttarg,density), fontsize=18) # x and y values for position of text are choosen for 864 particles
+  plt.show()
+  return
+
+
+def presplotter(prestime,averagedp,n):
+  fig = plt.figure()
+  ax1 = fig.add_subplot(2,1,1)
+  ax2 = fig.add_subplot(1,1,1)
+  fig.subplots_adjust(hspace=.35)
+  ax1.plot(range(len(prestime)),prestime)
+  ax2.plot(range(len(averagedp)),averagedp)
+  ax1.set_ylabel('pressure')
+  ax1.set_xlabel('time')
+  ax2.set_ylabel('pressure')
+  ax2.set_xlabel('time')
+  ax1.set_title('pressure evolution')
+  ax2.set_title('pressure evolution with pressure blocks')
+  ax2.text(1,2.715,'blocksize = %s'%(n),fontsize=12)
+  plt.show()
+  return
+
+def energyplotter(kenarray,potarray,totenarray):
+  plt.plot(range(len(kenarray)),kenarray)
+  plt.plot(range(len(potarray)),potarray)
+  plt.plot(range(len(totenarray)),totenarray)
+  plt.title('Energy evaluation')  
+  plt.xlabel('time')
+  plt.ylabel('energy')
+  plt.text(6800,-6500,r'N=%s, T=%s, $\rho$=%s'%(N,Ttarg,density), fontsize = 12)
+  plt.show()
+  return
+
+def renormplotter(kenarray,totenarray):
+  plt.plot(range(len(kenarray[0:2500])),kenarray[0:2500])
+  plt.plot(range(len(totenarray[0:2500])),totenarray[0:2500])
+  plt.title('Energy renormalization')  
+  plt.xlabel('time')
+  plt.ylabel('energy')
+  plt.text(1600,-5000,r'N=%s, T=%s, $\rho$=%s'%(N,Ttarg,density), fontsize = 12)
+  plt.show()
+  return
+
+
+# Calling the plot you want to make  
+
+#renormplotter(kenarray,totenarray)
+#energyplotter(kenarray,potarray,totenarray)
+#presplotter(prestime,averagedp,n)
+
+
+
+# Calculating the errors of the pressure and temperature
+
+mnp = np.mean(prestime)
+sdp = np.std(prestime)
+sdpb = np.std(averagedp)
+sdompb =  np.std(averagedp)/np.sqrt(len(averagedp))
+mnT = np.mean(kenarray[cutoff:9999])*(2/(3.*N))
+sdT = np.std(kenarray[cutoff:9999])*(2/(3.*N))
+sdomT = sdT/np.sqrt(len(kenarray[cutoff:9999]))
+
+print 'mnT=',mnT,'sdT=',sdT, 'sdomT=',sdomT, 'mnp=',mnp,'sdp=',sdp,'sdpb=',sdpb,'sdompb=',sdompb 

src/forces4.f90

@@ -0,0 +1,47 @@
+
+      subroutine calc_forces(N,pos,forces,pot,l,distances,presvirial)
+
+      !implicit none
+      integer, intent(in)  :: N
+      real(8), intent(in)  :: pos(N,3),l
+      real(8), intent(inout)  :: forces(N,3),pot, distances(N,N)
+      real(8), intent(inout) :: presvirial(1)        
+!f2py intent (in,out) :: forces, pot, distances, presvirial
+      real(8), parameter :: rc =3.2
+      real(8) :: dr_vec(3),partialforce(3),dr2,F,partialpot,cnt
+      integer :: i,j
+
+
+      presvirial = 0.0
+      forces(:,:) = 0.0 
+      pot = 0.0
+      distances(:,:) = 0.0
+      cnt = 0.0
+
+      !calculating the radial distance between each pair of particles
+      !Afterwards, do forces
+      do j = 1, N
+        do i = j+1, N
+          dr_vec = pos(i,:) - pos(j,:)
+          dr_vec = dr_vec - nint(dr_vec/l)*l
+          dr2 = dot_product(dr_vec,dr_vec)
+          distances(i,j) = sqrt(dr2)
+          distances(j,i) = sqrt(dr2)
+          if (dr2 < rc*rc) then
+            dr2 = 1/dr2
+            f = 2*dr2**7 - dr2**4
+            partialforce(:) = 24*dr_vec*f
+            forces(j,:) = forces(j,:) - partialforce(:)
+            forces(i,:) = forces(i,:) + partialforce(:)
+            presvirial = presvirial + dot_product(dr_vec,partialforce)
+            partialpot = 4*(dr2**6 - dr2**3)
+            pot = pot + partialpot
+            cnt = cnt + 1
+          end if
+        enddo
+      enddo
+      presvirial = presvirial/cnt
+
+      end subroutine calc_forces
+
+

src/init_sys.py

@@ -0,0 +1,58 @@
+#Initial momenta module
+import numpy as np
+import random as rn
+
+def init_mom(N, temp):
+  momenta = np.random.normal(0.0, np.sqrt(temp), (N,3))
+  avmom = np.zeros((1,3),dtype=float)
+  avmom = np.sum(momenta, axis=0)/N
+  momenta = momenta - avmom
+  return momenta
+
+def init_pos(N,density):
+  npdim = int(round((N/4)**(1./3.)))
+#  print npdim, 'npdim'
+  l = (N/density)**(1.0/3)
+  a = 0.5*l/npdim
+#  print npdim, l, a
+  pos = np.zeros((N,3), dtype = float)
+  cnt = 0
+  for k in range(npdim):
+    for j in range(npdim):
+      for i in range(npdim):
+        pos[cnt,:]   = [0+2*a*i,0+2*a*j,0+2*a*k]
+        pos[cnt+1,:] = [a+2*a*i,a+2*a*j,0+2*a*k]
+        pos[cnt+2,:] = [a+2*a*i,0+2*a*j,a+2*a*k]
+        pos[cnt+3,:] = [0+2*a*i,a+2*a*j,a+2*a*k]
+        cnt = cnt + 4
+
+  return pos, l, npdim
+
+def init_forc(N):
+  forces = np.zeros((N,3), dtype = float)
+  return forces
+
+def init_dist(N):
+  distances = np.zeros((N,N), dtype = float)
+  return distances
+
+def init_bins(nbins,lpnum):
+  bin_vec_tot = np.zeros(nbins, dtype = float)
+  finalbins = np.zeros(nbins, dtype = float)
+  return bin_vec_tot,finalbins
+
+def init_toten(lpnum):
+  toten = np.zeros((lpnum,1), dtype = float)
+  return toten
+
+def init_presvirialtime(lpnum):
+  presvirial = np.zeros((lpnum,1), dtype = float)
+  return presvirial
+
+def init_kenarray(lpnum):
+  kenarray = np.zeros((lpnum,1), dtype = float)
+  return kenarray
+
+def init_potarray(lpnum):
+  potarray = np.zeros((lpnum,1), dtype = float)
+  return potarray