diff options
Diffstat (limited to 'class')
| -rw-r--r-- | class/Makefile | 19 | ||||
| -rw-r--r-- | class/ke_compare.gnu | 15 | ||||
| -rw-r--r-- | class/ke_plot.gnu | 14 | ||||
| -rw-r--r-- | class/langdyn.f90 | 282 | ||||
| -rw-r--r-- | class/mean_plot.gnu | 20 | ||||
| -rw-r--r-- | class/movie_plot.gnu | 20 | ||||
| -rw-r--r-- | class/neighbourList.f90 | 63 | ||||
| -rw-r--r-- | class/putParticleInBox.f90 | 33 | ||||
| -rw-r--r-- | class/sortParticles.f90 | 56 |
9 files changed, 0 insertions, 522 deletions
diff --git a/class/Makefile b/class/Makefile deleted file mode 100644 index 17979f9..0000000 --- a/class/Makefile +++ /dev/null @@ -1,19 +0,0 @@ -FC = gfortran -FFLAGS = -O4 -fopenmp -fcheck=all -NPROC = $(shell nproc) -TARGET = langdyn.x - -all: clean $(TARGET) run - -$(TARGET): langdyn.f90 - $(FC) $(FFLAGS) -o $(TARGET) langdyn.f90 - -clean: - rm -f *.mod *.out *.gif $(TARGET) - -run: $(TARGET) - OMP_NUM_THREADS=$(NPROC) ./$(TARGET) - -plot: trajectories.out kinetic_energies.out - gnuplot movie_plot.gnu - gnuplot ke_plot.gnu diff --git a/class/ke_compare.gnu b/class/ke_compare.gnu deleted file mode 100644 index e09d021..0000000 --- a/class/ke_compare.gnu +++ /dev/null @@ -1,15 +0,0 @@ -set term cairolatex pdf size 5in,3in -set output "kinetic_energy.tex" -#set term x11 persist -set linetype 1 lc rgb '#E41A1C' pt 7 ps 0.3 -set linetype 2 lc rgb '#377EB8' pt 7 ps 0.3 -set linetype 3 lc rgb '#4DAF4A' pt 7 ps 0.3 -set linetype 4 lc rgb '#984EA3' pt 7 ps 0.3 - -set xlabel "Time [s]" -set ylabel "$\\sum$(Kinetic Energy) [J]" -set offsets graph 0.1, graph 0.1, graph 0.1, graph 0.1 -set key outside top horizontal center - -plot 'kinetic_without_slowdown' using 1:2 with linespoints lt 1 title "No External Force", \ - 'kinetic_with_slowdown' using 1:($2-4.50993810539) with linespoints lt 2 title "External Force" diff --git a/class/ke_plot.gnu b/class/ke_plot.gnu deleted file mode 100644 index c331f1c..0000000 --- a/class/ke_plot.gnu +++ /dev/null @@ -1,14 +0,0 @@ -#set term cairolatex pdf size 5in,3in -#set output "kinetic_energy.tex" -set term x11 persist -set linetype 1 lc rgb '#E41A1C' pt 7 ps 0.3 -set linetype 2 lc rgb '#377EB8' pt 7 ps 0.3 -set linetype 3 lc rgb '#4DAF4A' pt 7 ps 0.3 -set linetype 4 lc rgb '#984EA3' pt 7 ps 0.3 - -set xlabel "Time [s]" -set ylabel "$\\sum$(Kinetic Energy) [J]" -set offsets graph 0.1, graph 0.1, graph 0.1, graph 0.1 -set nokey - -plot 'kinetic_energies.out' using 1:2 with linespoints lt 1 notitle diff --git a/class/langdyn.f90 b/class/langdyn.f90 deleted file mode 100644 index b1880b2..0000000 --- a/class/langdyn.f90 +++ /dev/null @@ -1,282 +0,0 @@ -module globals -! Global variables -use omp_lib ! help the compiler find the OMP libraries -implicit none -integer :: n=1000 -double precision, parameter :: L=1.0d0 -double precision, parameter :: pi=2q0*asin(1q0) ! numerical constant -end module globals - -module Langevin -! Initialization and update rule for Langevin particles -use globals -implicit none -logical, allocatable, dimension(:) :: is_tracked -double precision :: dt,kT,g,m,sigma,eps,rc ! time step size and physical parameters -double precision :: pref1,pref2 ! auxiliary parameters -double precision, allocatable, dimension(:) :: x,y,vx,vy,ax,ay,vhx,vhy,x0,y0 ! particle positions, accellerations, velocities, half-step velocities, initial positions -contains -subroutine set_parameters - -! Set time step and physical parameters -dt=0.00001d0 ! time step size -kT=0.1d0 ! energy -g=0.1d0 ! drag coefficient -m=1d0 ! mass of the particles, can be normalized to 1. -sigma=1d-3 ! Potential parameters -eps=1d0 -rc=sigma*2d0**(1d0/6d0) ! Effective particle size - -! Set auxiliary parameters -pref1=g -pref2=sqrt(24d0*kT*g/dt) - -end subroutine set_parameters -subroutine initialize_particles -integer :: i, j, k, nx, ix, iy -double precision :: ran1(n),ran2(n),gr1(n),gr2(n),dist, lx, ly -logical, allocatable, dimension(:,:) :: is_filled -! Give particles initial position and velocity - - call random_seed() - call random_number(ran1) ! uses the built-in PRNG, easy but not very accurate - call random_number(ran2) - !x=L*(ran1-0.5d0) - !x0=x - !y=L*(ran2-0.5d0) - !y0=y - - !> Basic box stacking with no RNG - !nx = ceiling(sqrt(real(n))) - !dist = (0.9*L)/nx - !k = 1 - !outer: do i = 1,nx - ! do j = 1,nx - ! x(k) = i*dist-L/2 - ! y(k) = j*dist-L/2 - ! k = k + 1 - ! if (k.eq.n) exit outer - ! end do - !end do outer - - !> Box RNG on a finer mesh - nx = ceiling(L/rc) - dist = 0.9*L/nx - - allocate(is_filled(nx,nx)) - k = 1 - - !> Assuming (nx)*(nx) grid, first scale random numbers to a max of nx - do while (k .lt. n) - !> Random numbers - call random_number(lx) - call random_number(ly) - - !> Scale them to be grid points - ix = ceiling(lx*nx) - iy = ceiling(ly*nx) - - !> Check if grid is occupied - if (is_filled(ix,iy)) cycle - - !> If not, set the particle position and call it a day - x(k) = ix*rc - L/2 - y(k) = iy*rc - L/2 - - !> And mark the spot as filled - is_filled(ix,iy) = .true. - - !> Increment!! - k = k + 1 - end do - - ax=0d0 - ay=0d0 - call random_number(ran1) - call random_number(ran2) - gr1=sqrt(kT/(m))*sqrt(-2*log(ran1))*cos(2*pi*ran2) ! Box-Mueller transform - gr2=sqrt(kT/(m))*sqrt(-2*log(ran1))*sin(2*pi*ran2) - vx=gr1 - vy=gr2 - -end subroutine initialize_particles -end module Langevin - -module domainDecomposition - use globals - use Langevin - implicit none - integer, parameter :: b=4 - integer :: nbl(0:b*b-1,9) -contains - include "neighbourList.f90" - include "putParticleInBox.f90" - include "sortParticles.f90" -end module domainDecomposition - -module BC - ! Subroutines related to the boundary conditions - use globals - use langevin - implicit none -contains - subroutine impose_BC(i) - integer :: i - !> Recall we are inside an LxL box centered at the origin - - !> Top boundary - if(y(i) .GT. L/2) then - y(i) = L - y(i) - vhy(i) = -vhy(i) - end if - - !> Left boundary - if(x(i) .LT. -L/2) then - x(i) = -L - x(i) - vhx(i) = -vhx(i) - end if - - !> Bottom boundary - if(y(i) .LT. -L/2) then - y(i) = -L - y(i) - vhy(i) = -vhy(i) - end if - - !> Right boundary - if(x(i) .GT. L/2) then - x(i) = L - x(i) - vhx(i) = -vhx(i) - end if - - if(abs(x(i)).gt.L/2 .or. abs(y(i)).gt.L/2) is_tracked(i)=.false. - - end subroutine impose_BC -end module BC - -program main - use globals - use domainDecomposition -use Langevin -use BC -implicit none -integer :: i,j,lim(0:b*b,2),s,ns,p1,p2 -double precision :: t,t_max,m1,m2,rx,ry,dij,F -double precision :: wtime,begin,end,Tint,TsinceWrite -double precision, allocatable, dimension(:) :: ran1,ran2,xscrap,yscrap,vxscrap,vyscrap - -! Open files -open(11,file='trajectories.out') -open(12,file='means.out') -open(13,file='kinetic_energies.out') - -! Allocate arrays -allocate(x(n),y(n),vx(n),vy(n),ax(n),ay(n),vhx(n),vhy(n),x0(n),y0(n),is_tracked(n),ran1(n),ran2(n),xscrap(n),yscrap(n),vxscrap(n),vyscrap(n)) - -call buildNBL() - -is_tracked = .True. -t=0d0 -t_max=2.0d0 ! integration time -Tint = 0.001d0 -TsinceWrite=0d0 -call set_parameters -call initialize_particles - -begin = omp_get_wtime() -!call cpu_time(begin) - -! Conclusion: we need to re-order the loop like this: -! a. update half-step velocities -! b. update positions -! c. compute accellerations/forces -! d. update all velocities - -!$omp parallel -do while(t.lt.t_max) - ! one thread: fetch pseudo-random numbers - ! one thread: update velocity, position, impose BC - !$omp sections - !$omp section ! Write to disk - if(tSinceWrite.gt.Tint) then - !if(.true.) then - do i=1,n - write(11,*) xscrap(i),yscrap(i) - enddo - write(11,*) NEW_LINE('A') - write(12,*) t,sqrt(sum((xscrap-x)**2+(yscrap-y)**2, mask=is_tracked)/real(count(is_tracked),8)) - write(13,*) t,sum(0.5*vx**2 + 0.5*vy**2, mask=is_tracked) - tSinceWrite=0d0 - end if - xscrap=x - yscrap=y - vxscrap=vx - vyscrap=vy - ! write(12,*) t, sum(m*(vxscrap**2+vyscrap**2)/(2*n)) - - !$omp section - call random_number(ran1) - ran1=ran1-0.5d0 - call random_number(ran2) - ran2=ran2-0.5d0 - - !$omp section - vhx=vx+0.5d0*ax*dt - vhy=vy+0.5d0*ay*dt - x=x+vhx*dt - y=y+vhy*dt - do j=1,n - call impose_BC(j) - end do - call order(x,y,vx,vy,x0,y0,vhx,vhy,is_tracked,lim) ! BUGFIX 11/03: half-step velocities must also be re-ordered. - ax=0d0 ! Add forces here if any - ay=0d0 ! Add forces here if any - !$omp end sections - !$omp flush - !$omp single - ax=ax-pref1*vhx+pref2*ran1 - ay=ay-pref1*vhy+pref2*ran2 - !$omp end single - -! Our first attempt at parallelization of the code: run the computation of the distances and interaction forces on multiple threads: -!$omp do private(s,i,ns,p1,p2,rx,ry,dij,F) - do s=0,b*b-1 - do i=1,9 - if(nbl(s,i).eq.-1) exit - ns=nbl(s,i) ! BUGFIX 11/03: the loop counter was used as sector index instead of the entries of nbl. - do p1=lim(s,1),lim(s,2) - do p2=lim(ns,1),lim(ns,2) - if(p1.eq.p2) cycle - rx=x(p2)-x(p1) - ry=y(p2)-y(p1) - dij=sqrt(rx**2 + ry**2) - if(dij.lt.rc) then - F=4d0*eps*( -12d0*sigma**12/dij**13 + 6D0* sigma**6/dij**7 ) - ax(p1)=ax(p1)+F*rx/(dij*m) - ay(p1)=ay(p1)+F*ry/(dij*m) - end if - end do - end do - end do - end do - !$omp end do - !$omp single - vx=vhx+0.5d0*ax*dt - vy=vhy+0.5d0*ay*dt - t=t+dt - tSinceWrite=tSinceWrite+dt - !$omp end single -end do -!$omp end parallel - -end = omp_get_wtime() -!call cpu_time(end) -print *,'Wtime=',end-begin -if(lim(b*b,2).gt.lim(b*b,1)) print '(a5,i7,a8,i8,a11)','Lost ',lim(b*b,2)-lim(b*b,1),' out of ',n,' particles.' -! De-allocate arrays -deallocate(x,y,vx,vy,ax,ay,x0,y0,is_tracked,ran1,ran2,xscrap,yscrap,vxscrap,vyscrap) -! Close files -close(11) -close(12) -close(13) - -end program main diff --git a/class/mean_plot.gnu b/class/mean_plot.gnu deleted file mode 100644 index e86137d..0000000 --- a/class/mean_plot.gnu +++ /dev/null @@ -1,20 +0,0 @@ -#set term cairolatex pdf size 5in,3in -#set output "kinetic_energy.tex" -set term x11 persist -set linetype 1 lc rgb '#E41A1C' pt 7 ps 0.3 -set linetype 2 lc rgb '#377EB8' pt 7 ps 0.3 -set linetype 3 lc rgb '#4DAF4A' pt 7 ps 0.3 -set linetype 4 lc rgb '#984EA3' pt 7 ps 0.3 - -set xlabel "Time [s]" -set ylabel "RMSD [m]" -set logscale xy -set offsets graph 0.1, graph 0.1, graph 0.1, graph 0.1 -set nokey - -a=0 -cumsum(x)=(a=a+x,a) - -plot 'means.out' using 1:(cumsum($2)) with linespoints lt 1 notitle, \ - 'means.out' with linespoints notitle, \ - x**2 with lines dt 2 diff --git a/class/movie_plot.gnu b/class/movie_plot.gnu deleted file mode 100644 index 683753c..0000000 --- a/class/movie_plot.gnu +++ /dev/null @@ -1,20 +0,0 @@ -set term x11 noraise -set xrange [-0.6:0.6] -set yrange [-0.6:0.6] -set size square -set key outside top center -set object rectangle from -0.5,-0.5 to 0.5,0.5 dt 3 - -stats 'trajectories.out' nooutput - -do for [i=0:STATS_blocks-1] { - plot 'trajectories.out' index i with points pt 7 ps 0.5 title sprintf("Frame %i out of %i", i, STATS_blocks-1) - pause 0.0001 -} - -set term gif size 500,500 animate delay 0.001 -set output "wiggly.gif" - -do for [i=0:STATS_blocks-1:5] { - plot 'trajectories.out' index i with points pt 7 ps 0.5 title sprintf("Frame %i out of %i", i, STATS_blocks-1) -} diff --git a/class/neighbourList.f90 b/class/neighbourList.f90 deleted file mode 100644 index 8bf9c14..0000000 --- a/class/neighbourList.f90 +++ /dev/null @@ -1,63 +0,0 @@ -subroutine ONETO2(id, x, y) - implicit none - integer, intent(in) :: id - integer, intent(out) :: x, y - - y = id / b - x = mod(id, b) -end subroutine ONETO2 - -function TWOtoONE(x, y) result(id) - implicit none - integer, intent(in) :: x, y - integer :: id - - id = y * b + x -end function TWOtoONE - -function get_neighbour_ids(p) result (neighbours) - integer, intent(in) :: p !> The sequential position - integer, allocatable, dimension(:) :: neighbours !> The list of neighbours - integer :: i,j,k - integer :: x_cell, y_cell, x_test, y_test - integer :: max_list(9) - - !> Don't save the values between calls - k = 1 - max_list = -1 - - call ONETO2(p,x_cell,y_cell) - max_list(k) = p - - !print *, max_list - - !> Start by just getting all of the neighbours - do i = -1,1 - do j = -1,1 - if(i == 0 .and. j == 0) cycle !> Skip the cell itself - - !> The "test" coordinate - x_test = x_cell + i - y_test = y_cell + j - - !> If the coordinates are real (within the expected range), add them to the list: - if((x_test .ge. 0 .and. x_test .lt. b) .and. (y_test .ge. 0 .and. y_test .lt. b)) then - !> Increment the number of correct coordinates - k = k + 1 - !> Flip a zero to the new coordinate in the max list - max_list(k) = TWOtoONE(x_test, y_test) - end if - end do - end do - - !> Return a vector of all non-zero elements -! neighbours = pack(max_list, max_list>= 0.0) - neighbours = max_list ! entries -1 must be skipped as they correspond to non-existent neighbours at the edge or corner -end function get_neighbour_ids - -subroutine buildNBL() - integer :: sector - do sector=0,b*b-1 - nbl(sector,:) = get_neighbour_ids(sector) - end do -end subroutine buildNBL diff --git a/class/putParticleInBox.f90 b/class/putParticleInBox.f90 deleted file mode 100644 index d98a116..0000000 --- a/class/putParticleInBox.f90 +++ /dev/null @@ -1,33 +0,0 @@ - ! intakes coordinates and outputs index - integer function indexfxn(px, py) result(idx) - double precision, intent(in) :: px, py - double precision :: boxSize,yCoord - integer :: ix, iy - - !> Check if particle is in box - if(abs(px).GT.L/2 .OR. abs(py).GT.L/2) then - !> Particle is outside, don't track it - idx=b*b ! This is the "recycle bin" for untracked particles. - else - ! the sector size - ! M = floor(L / rc) - yCoord = - py - ! coordinate + half length of boundary [-L/2,L/2] -> [0,L] - ! divided by effective particle size to give coordinate - boxSize=L/real(b,8) - ix = floor((px + L/2.0d0)/boxSize) - iy = floor((yCoord + L/2.0d0)/boxSize) - - ! for particles on the boundary - ix = min(max(ix, 0), b-1) - iy = min(max(iy, 0), b-1) - if (ix < 0 .or. ix >= b .or. iy < 0 .or. iy >= b) then - print *, "ERROR: particle outside sector grid" - print *, "px, py =", px, py - stop - end if - - ! call rey's function - idx = TWOtoONE(ix,iy) - end if - end function indexfxn diff --git a/class/sortParticles.f90 b/class/sortParticles.f90 deleted file mode 100644 index 1645b5d..0000000 --- a/class/sortParticles.f90 +++ /dev/null @@ -1,56 +0,0 @@ -subroutine makeNR(nr,in) - integer :: particle,sector,nr(0:b*b),in(n) - nr = 0 - do particle=1,n - sector=indexfxn(x(particle),y(particle)) - if(sector.eq.b*b) is_tracked(particle)=.False. - in(particle)=sector - nr(sector)=nr(sector)+1 - end do - return -end subroutine makeNR - -subroutine order(x,y,vx,vy,x0,y0,vhx,vhy,is_tracked,lim) -! Order the list of positions by sector and find starting and final index for each sector -! In: x and y coordinates and velocities. Out: ordered lists x, y, vx and vy and array lim with one row -! for each sector, first column is start index, second is end index so that particles lim(k,1) through lim(k,2) reside in sector k. -integer :: lim(0:b*b,2),in(n),nr(0:b*b),ct(0:b*b),k -double precision :: x(n),y(n),d1(n),d2(n),d3(n),d4(n),d5(n),d6(n),d7(n),d8(n),vx(n),vy(n),x0(n),y0(n),vhx(n),vhy(n) -logical :: d9(n), is_tracked(n) - -call makeNR(nr,in) - -! Set loop limits based on the number of particles in each sector -lim(0,1)=1 -lim(0,2)=nr(0) - -do k=1,b*b - lim(k,1)=lim(k-1,2)+1 - lim(k,2)=lim(k-1,2)+nr(k) -end do - -! Re-order particle list -d1=x -d2=y -d3=vx -d4=vy -d5=x0 -d6=y0 -d7=vhx -d8=vhy -d9=is_tracked -ct=0 -do k=1,n - x(lim(in(k),1)+ct(in(k)))=d1(k) - y(lim(in(k),1)+ct(in(k)))=d2(k) - vx(lim(in(k),1)+ct(in(k)))=d3(k) - vy(lim(in(k),1)+ct(in(k)))=d4(k) - x0(lim(in(k),1)+ct(in(k)))=d5(k) - y0(lim(in(k),1)+ct(in(k)))=d6(k) - vhx(lim(in(k),1)+ct(in(k)))=d7(k) - vhy(lim(in(k),1)+ct(in(k)))=d8(k) - is_tracked(lim(in(k),1)+ct(in(k)))=d9(k) - ct(in(k))=ct(in(k))+1 -end do - -end subroutine order |