module density_profile_modulated_def_2d implicit none real*8, parameter :: R = 0.5d0 ! particle radius real*8, parameter :: Amp = 0.45d0 ! amplitude of the force real*8, parameter :: h_period = 0.45d0 ! period integer, parameter :: n_period = 3 ! not one because we need more to have a good resolution of q at low values real*8, parameter :: h_box = n_period*h_period ! box size real*8, parameter :: eta_in = 0.55d0 ! initial packing fraction (in the bulk) real*8, parameter :: eta_fin = 0.55d0 ! final packing fraction (in the bulk) integer, parameter :: niterations = 0 ! number of interations. From one iteration to the other phi increases untill the it reaches phi_fin integer, parameter :: ndft = 1 ! different type of density profile theory (1 == Eq.(9) of Roth et al. J. Chem. Phys. 136, 081101 (2012) ) real*8, parameter :: pi = 3.14159265358979323846264338327950288d0 integer, parameter :: n = 4096 ! number of points in the direction x (parallel to the modulation) must be a power of 2! real*8, parameter :: dx = h_box/n ! discretization in the x direction real*8, parameter :: dqx = 1.d0/h_box ! la rho avra` valori solo per: 2*pi / h_period integer, parameter :: nR = INT(R/dx + 0.0000000001d0) ! if -nR <= i <= nR --> -R <= x <= R real*8, parameter :: beta = 1.D0 ! inverse temperature real*8, parameter :: eps_old = 1.d0 ! ? !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!! GLOBAL QUANTITES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real*8 x(0:n-1) real*8 qx(0:n-1) real*8 wx(4,0:n-1) real*8 wqx(4,0:n-1) real*8 rhox(0:n-1) real*8 rhoqx(0:n-1) real*8 wd(4,0:n-1) real*8 :: eta, rho_g, rho_g_in real*8 :: pd(0:5,0:n-1), pdw(0:5,0:n-1) real*8 :: p(0:n-1), f(0:n-1), p1(0:n-1), pd1(0:n-1), p2(0:n-1), pd2(0:n-1), pdwsum(0:n-1) real*8 :: V(0:n-1) complex*16 :: h1(0:n-1), h2(0:n-1),h3(0:n-1) real*8 :: rho_new(0:n-1) save contains Subroutine print_rho () implicit none integer i open (21, file='rhoV_n.dat', status='unknown') !~ call contascorri(i,21) !~ write (unit=21, fmt="(a,f12.8)") "# rho_bulk :", rho_g !~ write (unit=21, fmt="(a,f12.8,a,f4.1,a,f5.3)") "# packing fraction :", eta, " period :",h_period, " Amp :",Amp do i = 0,n-1 WRITE (21,'(3es14.6)') x(i), rhox(i), V(i) enddo !~ write(21,*) ' ' close (21) end Subroutine print_rho Subroutine print_x (n,ff) implicit none integer i,n real*8 ff(0:n-1) do i = 0,n-1,8 write (*,'(I4,A1,8(1X,F10.4))') i,')', & ff(i),ff(i+1),ff(i+2),ff(i+3),ff(i+4),ff(i+5),ff(i+6),ff(i+7) enddo write(*,*) '-------------------------------------------------------' end Subroutine print_x Subroutine define_grid () implicit none integer i do i = 0,n-1 x(i) = -n/2*dx + i*dx qx(i) = i*dqx enddo end Subroutine define_grid Subroutine define_weight_planar_2d () implicit none integer i complex*16 y(0:n-1),ffty(0:n-1) ! Define weight functions ! wx(1,:) == w2 == w[2,2]*R ! wx(2,:) == w(0) == w0*2*pi*R ! wx(3,:) == w(1) ! wx(4,:) == w[1,1] wx(:,:) = 0.d0 wx(1,n/2) = 2*R wx(2,n/2) = 2.d0 wx(3,n/2) = 0.d0 wx(4,n/2) = 0.d0 do i = n/2+1,n/2+nR wx (1,i) = 2*dsqrt(R**2 - x(i)**2) ! real and even wx (2,i) = 2*R/ dsqrt(R**2 - x(i)**2) ! real and even wx (3,i) = 2*x(i)/ dsqrt(R**2 - x(i)**2) ! real and odd wx (4,i) = 2*x(i)*x(i)/R /dsqrt(R**2 - x(i)**2) ! real and even wx (1,n-i) = wx(1,i) wx (2,n-i) = wx(2,i) wx (3,n-i) = -wx(3,i) wx (4,n-i) = wx(4,i) ! cure divergences if (i.eq.n/2+nR .and. x(i).ge.R) then wx (2,i) = 0.d0 wx (2,n-i) = 0.d0 wx (3,i) = 0.d0 wx (3,n-i) = 0.d0 wx (4,i) = 0.d0 wx (4,n-i) = 0.d0 endif enddo ! and theirs FT (FT(f) = int f(x) exp (2 pi i q x) dx wqx (1,0) = pi*R/2 wqx (2,0) = 2*pi*R wqx (3,0) = 0.d0 wqx (4,0) = pi*R do i = 1,n/2-1 wqx (1,i) = R*BESSJ1(2*pi*qx(i)*R) / qx(i) ! real and even wqx (2,i) = 2*pi*R*BESSJ0(2*pi*qx(i)*R) ! real and even wqx (3,i) = 2*pi*R*BESSJ1(2*pi*qx(i)*R) ! imaginary and odd wqx (4,i) = 2*pi*R*BESSJ0(2*pi*qx(i)*R)-BESSJ1(2*pi*qx(i)*R)/qx(i) ! real and even ! we have to symmetrize to take into account symmetries wqx (1,n-i) = wqx (1,i) wqx (2,n-i) = wqx (2,i) wqx (3,n-i) = -wqx (3,i) wqx (4,n-i) = wqx (4,i) enddo wqx (1,n/2) = R*BESSJ1(2*pi*qx(n/2)*R) / qx(n/2) ! real and even wqx (2,n/2) = 2*pi*R*BESSJ0(2*pi*qx(n/2)*R) ! real and even wqx (3,n/2) = 0. ! imaginary and odd wqx (4,n/2) = 2*pi*R*BESSJ0(2*pi*qx(n/2)*R)-BESSJ1(2*pi*qx(n/2)*R)/qx(n/2) ! real and even ! better to do the FFT because with this function there are problems... !~ do i=0,n-1 !~ y(i) = dcmplx(wx(1,i),0.d0) ! real and even !~ enddo !~ call fft_mik1 (n,x,y,qx,ffty,1,1,1) !~ wqx(1,:) = dreal(ffty(:)) ! real and even !~ do i=0,n-1 !~ y(i) = dcmplx(wx(2,i),0.d0) ! real and even !~ enddo !~ call fft_mik1 (n,x,y,qx,ffty,1,1,1) !~ wqx(2,:) = dreal(ffty(:)) ! real and even !~ do i=0,n-1 !~ y(i) = dcmplx(wx(3,i),0.d0) ! real and odd !~ enddo !~ call fft_mik1 (n,x,y,qx,ffty,1,0,1) !~ wqx(3,:) = dimag(ffty(:)) ! imaginary and odd !~ do i=0,n-1 !~ y(i) = dcmplx(wx(4,i),0.d0) ! real and even !~ enddo !~ call fft_mik1 (n,x,y,qx,ffty,1,1,1) !~ wqx(4,:) = dreal(ffty(:)) ! real and even end Subroutine define_weight_planar_2d Subroutine define_potential () implicit none integer i do i = 0,n-1 V(i) = Amp*dcos(2*pi*x(i)/h_period) enddo end Subroutine define_potential FUNCTION BESSJ0 (X) IMPLICIT NONE REAL*8 X,BESSJ0,AX,FR,FS,Z,FP,FQ,XX ! This subroutine calculates the First Kind Bessel Function of ! order 0, for any real number X. The polynomial approximation by ! series of Chebyshev polynomials is used for 0 i) then c1 = g(j) g(j) = g(i) g(i) = c1 endif ! Determination of j: Use previous j and add bits from the "left" side ! In case of overflow subtract the bits m = n/2 do while (m > 1 .AND. j >= m) j = j - m m = m/2 enddo j = j + m enddo ! Danielson-Lanczos nstep = 1 do while (nstep < n) do i = 0, nstep-1 ex1 = EXP( nrev*(i)*(0,1)*pi/nstep ) ex2 = EXP( nrev*(i+nstep)*(0,1)*pi/nstep ) do j = i, n-1, 2*nstep h(j) = g(j) + g(j+nstep) * ex1 h(j+nstep) = g(j) + g(j+nstep) * ex2 enddo enddo nstep = 2*nstep g = h enddo end subroutine fft_c subroutine contascorri (N,a) implicit none ! conta le linee N del file aperto alla unità a integer N,a,status character (len=100) line N=0 do read (a,*,iostat=status) line if (status.lt.0) exit N = N+1 enddo end subroutine contascorri end module density_profile_modulated_def_2d