Skip to content

Instantly share code, notes, and snippets.

@astrojuanlu

astrojuanlu/LJSim.f90

Forked from anonymous/Numba_LJSim.py
Last active February 13, 2016 21:25
Show Gist options
  • Save astrojuanlu/fe4c54ffc7aebee2c877 to your computer and use it in GitHub Desktop.
Save astrojuanlu/fe4c54ffc7aebee2c877 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
genconfig.py
confile = open("configuration.dat","w")
dimension = 4.0
bounds = int(dimension/2.0)
for i in range(-bounds,bounds):
for j in range(-bounds,bounds):
for k in range(-bounds,bounds):
confile.write(" " + str(i) + " " + str(j) + " " + str(k) + "\n")
Raw
LJSim.f90
!============================================
program LJSim
implicit none
interface
subroutine detailed_ecalc(E_T,x,y,z,boxl,boxl2)
implicit none
real(kind(0.0d0)), intent(in) :: x(:), y(:), z(:)
real(kind(0.0d0)), intent(in) :: boxl, boxl2
real(kind(0.0d0)), intent(out) :: E_T
end subroutine
end interface
interface
subroutine shift_ecalc(E_Diff, npart, x, y, z, dx, dy, dz, nmove, boxl,boxl2)
implicit none
integer, intent(in) :: nmove, npart
real(kind(0.0d0)), intent(in) :: x(:), y(:), z(:)
real(kind(0.0d0)), intent(in) :: dx, dy, dz
real(kind(0.0d0)), intent(in) :: boxl,boxl2
real(kind(0.0d0)), intent(out) :: E_Diff
end subroutine
end interface
integer :: i,j, indx
integer :: npart, ncycle, nmove
real(kind(0.0d0)) :: boxl, boxl2
real(kind(0.0d0)),allocatable :: x(:), y(:), z(:)
real(kind(0.0d0)) :: E_T, E_Diff,E_Final, temp, beta
real(kind(0.0d0)) :: dx, dy, dz
real(kind(0.0d0)) :: timeStart, timeEnd
real(kind(0.0d0)) :: grnd
boxl = 22.5d0
boxl2 = boxl/2d0
temp = 1.8d0
beta = 1d0/temp
open(unit=2, file="configuration.dat")
npart = 0
do i = 1, nint(1d7)
read(2,*,end=20)
npart = npart + 1
enddo
20 continue
if(npart .eq. 0) then
stop "Configuration File is empty!"
endif
rewind(2)
allocate(x(1:npart))
allocate(y(1:npart))
allocate(z(1:npart))
do i = 1, npart
read(2,*) x(i), y(i), z(i)
enddo
ncycle = nint(npart*1d5)
E_T = 0d0
call detailed_ecalc(E_T,x,y,z,boxl,boxl2)
open(unit=35, file="Traj.xyz")
write(35,*) npart
write(35,*)
do i = 1, npart
write(35,*) "Ar", x(i), y(i), z(i)
enddo
write(*,*) "Number of Particles:", npart
write(*,*) "Initial Energy:", E_T
call cpu_time(timeStart)
do indx = 1, ncycle
nmove = floor(npart*grnd()+1d0)
dx = x(nmove) + 0.2d0*(2.0d0*grnd()-1.0d0)
dy = y(nmove) + 0.2d0*(2.0d0*grnd()-1.0d0)
dz = z(nmove) + 0.2d0*(2.0d0*grnd()-1.0d0)
if(abs(dx) .gt. boxl) then
dx = dx - sign(boxl,dx)
endif
if(abs(dy) .gt. boxl) then
dy = dy - sign(boxl,dy)
endif
if(abs(dz) .gt. boxl) then
dz = dz - sign(boxl,dz)
endif
call shift_ecalc(E_Diff, npart, x, y, z, dx, dy, dz, nmove, boxl,boxl2)
if(E_Diff .lt. 0.0d0) then
E_T = E_T + E_Diff
x(nmove) = dx
y(nmove) = dy
z(nmove) = dz
else
if(exp(-beta*E_Diff) .gt. grnd()) then
E_T = E_T + E_Diff
x(nmove) = dx
y(nmove) = dy
z(nmove) = dz
endif
endif
if(mod(indx,10000) .eq. 0) then
write(*,*) indx, E_T
endif
if(mod(indx,20000) .eq. 0) then
write(35,*) npart
write(35,*)
do i = 1, npart
write(35,*) "Ar", x(i), y(i), z(i)
enddo
! flush(35)
endif
enddo
call detailed_ecalc(E_Final,x,y,z,boxl,boxl2)
call cpu_time(timeEnd)
write(35,*) npart
write(35,*)
do i = 1, npart
write(35,*) "Ar", x(i), y(i), z(i)
enddo
close(35)
open(unit=45, file="final_configuration.dat")
do i = 1, npart
write(45,*) x(i), y(i), z(i)
enddo
write(*,*) "Culmaltive Energy:", E_T
write(*,*) "Final Energy:", E_Final
write(*,*) "Total Time:", timeEnd-timeStart
end program
!============================================
subroutine detailed_ecalc(E_T,x,y,z,boxl,boxl2)
implicit none
real(kind(0.0d0)), intent(in) :: x(:), y(:), z(:)
real(kind(0.0d0)), intent(in) :: boxl, boxl2
real(kind(0.0d0)), intent(out) :: E_T
integer :: i,j, npart
real(kind(0.0d0)) :: rx, ry, rz, r_sq, LJ
npart = size(x)
E_T = 0d0
do i = 1, npart-1
do j = i+1, npart
rx = x(i) - x(j)
ry = y(i) - y(j)
rz = z(i) - z(j)
if(abs(rx) .gt. boxl2) then
rx = rx - sign(boxl,rx)
endif
if(abs(ry) .gt. boxl2) then
ry = ry - sign(boxl,ry)
endif
if(abs(rz) .gt. boxl2) then
rz = rz - sign(boxl,rz)
endif
r_sq = rx*rx + ry*ry + rz*rz
if(r_sq .lt. 49d0) then
LJ = (1.0d0/r_sq)
LJ = LJ * LJ * LJ
LJ = 4.0d0*LJ*(LJ-1.0d0)
E_T = E_T + LJ
endif
enddo
enddo
end subroutine
!============================================
subroutine shift_ecalc(E_Diff, npart, x, y, z, dx, dy, dz, nmove, boxl,boxl2)
implicit none
integer, intent(in) :: nmove, npart
real(kind(0.0d0)), intent(in) :: x(:), y(:), z(:)
real(kind(0.0d0)), intent(in) :: dx, dy, dz
real(kind(0.0d0)), intent(in) :: boxl,boxl2
real(kind(0.0d0)), intent(out) :: E_Diff
integer :: i,j
real(kind(0.0d0)) :: rx, ry, rz, r_sq, LJ
E_Diff = 0d0
do j = 1, npart
if(j .ne. nmove) then
rx = dx - x(j)
ry = dy - y(j)
rz = dz - z(j)
if(abs(rx) .gt. boxl2) then
rx = rx - sign(boxl,rx)
endif
if(abs(ry) .gt. boxl2) then
ry = ry - sign(boxl,ry)
endif
if(abs(rz) .gt. boxl2) then
rz = rz - sign(boxl,rz)
endif
r_sq = rx*rx + ry*ry + rz*rz
if(r_sq .lt. 49d0) then
LJ = (1d0/r_sq)
LJ = LJ * LJ * LJ
LJ = 4d0*LJ*(LJ-1d0)
E_Diff = E_Diff + LJ
endif
rx = x(nmove) - x(j)
ry = y(nmove) - y(j)
rz = z(nmove) - z(j)
if(abs(rx) .gt. boxl2) then
rx = rx - sign(boxl,rx)
endif
if(abs(ry) .gt. boxl2) then
ry = ry - sign(boxl,ry)
endif
if(abs(rz) .gt. boxl2) then
rz = rz - sign(boxl,rz)
endif
r_sq = rx*rx + ry*ry + rz*rz
if(r_sq .lt. 49d0) then
LJ = (1d0/r_sq)
LJ = LJ * LJ * LJ
LJ = 4d0*LJ*(LJ-1d0)
E_Diff = E_Diff - LJ
endif
endif
enddo
end subroutine
!============================================
Raw
Numba_LJSim.py
from numba import jit
import numpy as np
import random as rng
import math
import time
@jit(cache=True)
def main():
boxl = 22.5
boxl2 = boxl/2.0
coords = np.loadtxt("configuration.dat")
print(coords.shape)
npart = coords.shape[0]
temp = 1.8
beta = 1.0/temp
print(npart)
ncycle = int(10**5 * npart)
E_T = 0.0
E_T = detailed_ecalc(npart, coords, boxl, boxl2)
traj = open("Traj.xyz","w")
print(npart,file=traj)
print("\n",file=traj)
for i in range(0,npart):
print("Ar",coords[i][0],coords[i][1],coords[i][2],file=traj)
print("Initial Energy: ", E_T)
print("Simulation Start!")
start_time = time.time()
coords = montecarlo_loop(ncycle,coords, beta, E_T, boxl, boxl2)
E_Final = detailed_ecalc(npart, coords, boxl, boxl2)
elapsed_time = time.time() - start_time
print("Total Time: ", elapsed_time)
finalconfig = open("final_configuration.dat", "w")
for i in range(0,npart):
print(coords[i][0],coords[i][1],coords[i][2],file=finalconfig)
print(npart,file=traj)
print("\n",file=traj)
for i in range(0,npart):
print("Ar",coords[i][0],coords[i][1],coords[i][2],file=traj)
@jit(cache=True)
def montecarlo_loop(ncycle,coords, beta, E_T, boxl, boxl2):
npart = coords.shape[0]
disp = np.ndarray(shape=[1,3], dtype=float)
e_diff = 0.0
nmove = 0
for indx in range(0,ncycle+1):
nmove = math.floor(npart*rng.random())
for i in range(0,3):
disp[0][i] = 0.2*(2.0*rng.random()-1.0)
if math.fabs(disp[0][i] + coords[nmove][i]) > boxl:
disp[0][i] = disp[0][i] - math.copysign(boxl,disp[0][i])
e_diff = shift_ecalc(npart,coords,nmove,disp, boxl, boxl2)
if e_diff < 0.0:
E_T += e_diff
coords[nmove][0] += disp[0][0]
coords[nmove][1] += disp[0][1]
coords[nmove][2] += disp[0][2]
else:
if math.exp(-beta*e_diff) > rng.random():
E_T += e_diff
coords[nmove][0] += disp[0][0]
coords[nmove][1] += disp[0][1]
coords[nmove][2] += disp[0][2]
# if indx%10000 == 0:
# fortprint.screenprint(indx,E_T)
# print(indx,E_T)
# if indx%100000 == 0:
# print(npart,file=traj)
# print(file=traj)
# for i in range(0,npart):
# print("Ar",coords[i][0],coords[i][1],coords[i][2],file=traj)
print("Culmaltive Energy: ", E_T)
E_Final = detailed_ecalc(npart, coords, boxl, boxl2)
print("Final Energy: ", E_Final)
return coords
@jit(nopython=True,cache=True)
def shift_ecalc(npart,coords,nmove,disp, boxl, boxl2):
energy = 0.0
for j in range(0,npart):
if j != nmove:
rx = coords[nmove][0]+disp[0][0] - coords[j][0]
ry = coords[nmove][1]+disp[0][1] - coords[j][1]
rz = coords[nmove][2]+disp[0][2] - coords[j][2]
if math.fabs(rx) > boxl2:
rx -= math.copysign(boxl,rx)
if math.fabs(ry) > boxl2:
ry -= math.copysign(boxl,ry)
if math.fabs(rz) > boxl2:
rz -= math.copysign(boxl,rz)
r_sq = rx*rx + ry*ry + rz*rz
if r_sq < 49.0:
LJ = (1.0/r_sq)
LJ = LJ*LJ*LJ
LJ = 4.0 * LJ * (LJ-1.0)
energy = energy + LJ
rx = coords[nmove][0] - coords[j][0]
ry = coords[nmove][1] - coords[j][1]
rz = coords[nmove][2] - coords[j][2]
if math.fabs(rx) > boxl2:
rx -= math.copysign(boxl,rx)
if math.fabs(ry) > boxl2:
ry -= math.copysign(boxl,ry)
if math.fabs(rz) > boxl2:
rz -= math.copysign(boxl,rz)
r_sq = rx*rx + ry*ry + rz*rz
if r_sq < 49.0:
LJ = (1.0/r_sq)
LJ = LJ*LJ*LJ
LJ = 4.0*LJ*(LJ-1.0)
energy = energy - LJ
return energy
@jit(nopython=True,cache=True)
def detailed_ecalc(npart,coords, boxl, boxl2):
energy = 0.0
for i in range(0,npart-1):
for j in range(i+1,npart):
rx = coords[i][0] - coords[j][0]
ry = coords[i][1] - coords[j][1]
rz = coords[i][2] - coords[j][2]
if math.fabs(rx) > boxl2:
rx -= math.copysign(boxl,rx)
if math.fabs(ry) > boxl2:
ry -= math.copysign(boxl,ry)
if math.fabs(rz) > boxl2:
rz -= math.copysign(boxl,rz)
r_sq = rx*rx+ ry*ry + rz*rz
if r_sq < 49.0:
LJ = (1.0/r_sq)
LJ = LJ*LJ*LJ
LJ = 4.0*LJ*(LJ-1.0)
energy = energy + LJ
return energy
main()
Raw
Random.f90
!=======================================================
double precision function grnd()
implicit none
real(kind(0.0d0)) :: r
call RANDOM_NUMBER(r)
grnd = r
end function
!=======================================================
subroutine sgrnd(seed)
implicit none
integer,intent(inout) :: seed
integer :: i,n
integer, allocatable :: tempSeed(:)
call RANDOM_SEED(size=n)
allocate(tempSeed(1:n))
tempSeed = seed + 37 * (/ (i - 1, i = 1, n) /)
call RANDOM_SEED(put=tempSeed)
deallocate(tempSeed)
end subroutine
!=======================================================
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment