cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc program DPD_POLY_2D cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c This code is for a single homopolymer in a solvent in 2D c by dissipative particle dynamics (DPD). c c We do the following (main points): c c - generate the initial configuration c - equilibrate the system c - do the dynamics in equilibrium c -- calculate the temperature of the system c -- for the chain, calculate the radius of gyration squared c as well as the end-to-end distance squared c -- determine the velocity autocorrelation function which c is related to the diffusion coefficient of the chain c - close the simulation c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c NOTE: The input parameters (that can be changed) are c given on lines 146 - 161 below. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit none ! MAXIMUM SIZE OF SOME VARIABLES: integer ! ------------------------------------------ + Npmax, ! Maximum number of particles in a system + maxnab, ! Maximum size of the Verlet table + maxt_vacf, ! Maximum number of data points for the vacf + Nclen, ! Maximum number of monomers in a chain + Ifile ! Numbering of snapshots starts from Ifile: ! If Ifort = 50, then the first snapshot is ! fort.50, the second one is fort.51 etc. parameter( Npmax = 10000 ) parameter( maxnab = 100*Npmax ) parameter( maxt_vacf = 100 ) parameter( Nclen = 20 ) parameter( Ifile = 50 ) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc integer ! Defines the two particle types + n_flag(Npmax) ! ('1' water, '2' monomer) double precision + rX(Npmax), ! Positions of particles (x coordinate) + rY(Npmax), ! Positions of particles (y coordinate) + rX0(Npmax), ! Saved positions + rY0(Npmax), ! Saved positions + vX(Npmax), ! Velocities of particles (x coordinate) + vY(Npmax), ! Velocities of particles (y coordinate) + FcX(Npmax), ! Conservative forces + FcY(Npmax), ! + FdX(Npmax), ! Dissipative forces + FdY(Npmax), ! + FrX(Npmax), ! Random forces + FrY(Npmax) ! integer + it_vacf, ! For the velocity autocorrelation function + ivacf_counter, + it_bound(maxt_vacf + 1), + ivacf_loop(-maxt_vacf:maxt_vacf) double precision + velX(maxt_vacf), ! For the velocity autocorrelation function + velY(maxt_vacf), + dvacf(0:(maxt_vacf-1)), + dcs_vacf(0:(maxt_vacf-1)) integer + npoint(Npmax), ! For the Verlet neighbor tables + list(maxnab), + nVsize, nVintv double precision + alpha(2,2), ! Conservative force strengths + gamma, ! Dissipative force strength + sigma ! Random force strength logical + update ! Logical variable for the Verlet tables integer + Np, ! Total particle number + maxstp, ! Number of time steps after equilibration + incstp, ! Number of steps between samples + N_water, ! Number of water particles + len_chain, ! Number of monomers in a polymer chain + isteps, ! Number of time steps for equilibration + issstp ! Number of steps between snapshots double precision + dLx,dLy, ! Size of the simulation box in 2D + dLxINV, dLyINV, ! Inverse system size + rho, rho0, ! Density + dNp, ! Total particle number + dNpINV, ! Inverse particle number + rfac, ! Scaling factor of the random force + Tfac, ! Scaling factor of temperature + dr_gyr, ! Radius of gyration + dr_ee, ! End-to-end distance of the chain + dt, ! Time step for the molecular dynamics + tmax, ! Time of the simulation after equilibration + vCM, ! Center-of-mass velocity + temp, ! Temperature + pm, ! Particle mass + pmINV, ! Inverse particle mass + skin, ! Skin radius for the Verlet tables + skinSQ, ! Skin radius squared + dispmx, ! Displacement length + d_bond ! Harmonic force constant for the chain integer + iseed ! For the random number generator real ! For the random number generator + R2S, R2INIS, dummy integer ! Various integer variables + i, ip, it, + nconf double precision ! Various double precision variables + vol, dtH, dtRTH c --------------------------------------------------------------- c c Here we define the INPUT PARAMETERS for the model study c c --------------------------------------------------------------- c c NOTE! There are some restrictions: c c - "len_chain" should not be larger than "Nclen" c - The total number of particles should not exceed "Npmax". c In practice this means that c " rho0 * dLx * dLy " should be smaller than or equal to "Npmax". c c --------------------------------------------------------------- dLx = 16.d0 ! system size in x-direction dLy = 16.d0 ! system size in y-direction rho0 = 3.d0 ! total density of all particles len_chain = 20 ! size of the chain d_bond = 10.d0 ! harmonic bonding btw. nn beads in a chain pm = 1.d0 ! particle mass (identical for all particles) alpha(1,1) = 10.d0 ! water - water interaction alpha(1,2) = 3.d0 ! water - monomer alpha(2,2) = 10.d0 ! monomer - monomer sigma = 3.d0 ! dpd parameter iseed = 619099 ! seed for a random number generator dt = 1.0d-2 ! time step for the integrator isteps = 500 ! number of steps to equilibrate the system maxstp = 2500 ! number of time steps after equilibration incstp = 20 ! how often do we calculate quantities issstp = 100 ! how often do we take snapshots c --------------------------------------------------------------- OPEN(3,STATUS='NEW', ! Open file for the general data + FILE='data.dat', + FORM='FORMATTED', + ACCESS='SEQUENTIAL') OPEN(21,STATUS='NEW', ! Open file for the radius of gyration + FILE='rg2.dat', + FORM='FORMATTED', + ACCESS='SEQUENTIAL') if((maxstp/incstp).ge.maxt_vacf) then OPEN(24,STATUS='NEW', ! Open file for the velocity + FILE='vcf.dat', ! autocorrelation function + FORM='FORMATTED', + ACCESS='SEQUENTIAL') end if !-------------------- ! Fix some constants: !-------------------- dLxINV = 1.d0/dLx ! Inverse system size dLyINV = 1.d0/dLy tmax = dt*dble(maxstp) ! Time of the production run vol = dLx*dLy ! Volume (area) of the system Np = nint( rho0*vol ) ! Total number of particles if (Np.gt.Npmax) + stop 'Fatal Error: Np > Npmax' if (len_chain.gt.Nclen) + stop 'Fatal Error: len_chain > Nclen' dNp = dble(Np) dNpINV = 1.d0/dNp ! Inverse particle number rho = dNp/vol ! Total density of the system N_water = Np - len_chain ! Number of water particles pmINV = 1.d0/pm ! Inverse particle mass alpha(2,1) = alpha(1,2) ! Symmetric conservative interaction gamma = ! The strength of the dissipative force + 0.5d0*sigma*sigma ! (using the fluctuation-dissipation theorem) dtH = 0.5d0*dt ! Time steps for the integrator dtRTH = 0.5d0*dsqrt(dt) skin = 2.0d0 ! Skin radius for the Verlet neighbor table skinSQ = skin*skin rfac = dsqrt(3.d0) ! Scaling of random forces Tfac = ! Factor for temperature control + 1.d0/(2.d0*dNp - 2.d0) do it=0,(maxt_vacf - 1) ! Here we initialize a few variables dvacf(it) = 0.d0 ! to examine the decay of the end do ! velocity correlation function of ! the polymer. This is related to do i=1,maxt_vacf ! the tracer diffusion coefficient ivacf_loop(i) = i ! of the polymer motion. end do do i=-maxt_vacf,0 ivacf_loop(i) = i + maxt_vacf end do it_vacf = 0 ivacf_counter = 0 do it=1,maxt_vacf it_bound(it) = it end do it_bound(maxt_vacf + 1) = 1 dummy = R2INIS( iseed ) ! Initialize the random number generator do ip=1,N_water ! Fix the particle types. n_flag(ip) = 1 ! Type '1' corresponds to solvent (water), end do ! while the type '2' represents monomers ! in a polymer chain. do ip = (N_water + 1), Np n_flag(ip) = 2 end do it = 0 ! Set time t = 0 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Time to start the actual DPD simulation c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc call GENCON(Npmax, ! Generate the initial configuration + rX,rY,vX,vY, + N_water,Np, + len_chain, + dLx,dLy, + dLxINV,dLyINV, + dNpINV,pm,pmINV, + Tfac) call R_GYR(Npmax, ! Then calculate the radius of gyration + Nclen,rX,rY, ! and the end-to-end distance at + Np,N_water, ! time t = 0 + len_chain,dLx,dLy, + dLxINV,dLyINV, + dr_gyr,dr_ee) write(21,755) + it, dr_gyr, dr_ee ! Write the results for the radius of call FLUSH(21) ! gyration and end-to-end distance update = .true. ! Make sure that the Verlet table ! will be constructed call FCDR(Npmax, ! Calculate all forces for + maxnab,rX,rY, ! the initial configuration. + rX0,rY0,vX,vY, + FcX,FcY,FdX,FdY, + FrX,FrY,npoint, + list,n_flag,alpha, + Np,N_water,len_chain, + dLx,dLy,dLxINV,dLyINV, + rfac,update,nVintv, + nVsize,gamma,sigma, + d_bond,skinSQ) nVintv = 0 c --------------------------------------------------------------- do it = 1, isteps ! EQUILIBRATE THE SYSTEM! call INTV(Npmax, ! Integrate velocities over dt/2 + Np,vX,vY, + FcX,FcY,FdX,FdY, + FrX,FrY,pmINV, + dtH,dtRTH) call INTR(Npmax, ! Integrate particle positions over dt + Np,rX,rY, + vX,vY,dt) call CHECK(Npmax, ! Check if we should update the + Np,rX,rY,rX0,rY0, ! Verlet neighbor table + skin,update) if (update) then call FOLDB(Npmax, ! Particles have wandered outside + Np,rX,rY, ! the simulation box, so let us + dLx,dLy, ! shift them back inside + dLxINV,dLyINV) endif call FCDR(Npmax, ! Calculate all forces + maxnab,rX,rY, ! (conservative, dissipative, and + rX0,rY0,vX,vY, ! random ones) + FcX,FcY,FdX,FdY, + FrX,FrY,npoint, + list,n_flag,alpha, + Np,N_water,len_chain, + dLx,dLy,dLxINV,dLyINV, + rfac,update,nVintv, + nVsize,gamma, + sigma,d_bond, + skinSQ) call INTV(Npmax,Np, ! Integrate velocities over dt/2 + vX,vY,FcX,FcY, ! (Note: here we use the so-called + FdX,FdY,FrX,FrY, ! DPD-VV integration scheme) + pmINV,dtH,dtRTH) call FDR(Npmax,maxnab, ! Calculate dissipative forces + rX,rY,vX,vY, + FdX,FdY,npoint, + list,Np,dLx,dLy, + dLxINV,dLyINV, + gamma) if (mod(it,incstp).eq.0) then call VCMTP(Npmax, + Np,vX,vY, ! Time to take samples of the + dNpINV,pm, ! center-of-mass velocity and + pmINV,Tfac, ! the temperature of the system + vCM,temp) write(3,1000) ! Write the results + it, + sngl(dble(it)*dt), + vCM, sngl(temp) call FLUSH(3) call R_GYR(Npmax, ! Calculate the radius of gyration + Nclen,rX,rY, ! squared and the end-to-end distance + Np,N_water, + len_chain, + dLx,dLy, + dLxINV,dLyINV, + dr_gyr,dr_ee) write(21,755) ! Write the results + it, dr_gyr, + dr_ee call FLUSH(21) endif enddo ! EQUILIBRATION HAS BEEN DONE c --------------------------------------------------------------- nconf = 0 ! Initialize counter for configurations call snaps(Npmax,Np, ! Take a snapshot of the polymer chain + N_water,Nclen, + len_chain,nconf, + dLx,dLy, + dLxINV,dLyINV, + rX,rY,Ifile) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Then do the dynamics in equilibrium: c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 755 format(1x,1i7,1x,1f10.6,1x,1f10.6) 1000 format(i10, 1x,f10.4, 1x,d12.6, 1x,f10.6) do it = (isteps + 1), (isteps + maxstp) call INTV(Npmax, ! Integrate velocities over dt/2 + Np,vX,vY, + FcX,FcY,FdX,FdY, + FrX,FrY,pmINV, + dtH,dtRTH) call INTR(Npmax, ! Integrate particle positions over dt + Np,rX,rY, + vX,vY,dt) call CHECK(Npmax, ! Check if we should update the + Np,rX,rY,rX0,rY0, ! Verlet neighbor table + skin,update) if (update) then call FOLDB(Npmax, ! Particles have wandered outside + Np,rX,rY, ! the simulation box, so let us + dLx,dLy, ! shift them back inside + dLxINV,dLyINV) endif call FCDR(Npmax, ! Calculate all forces + maxnab,rX,rY, ! (conservative, dissipative, and + rX0,rY0,vX,vY, ! random forces) + FcX,FcY,FdX,FdY, + FrX,FrY,npoint, + list,n_flag,alpha, + Np,N_water,len_chain, + dLx,dLy,dLxINV,dLyINV, + rfac,update,nVintv, + nVsize,gamma,sigma, + d_bond,skinSQ) call INTV(Npmax,Np, ! Integrate velocities over dt/2 + vX,vY,FcX,FcY, + FdX,FdY,FrX,FrY, + pmINV,dtH,dtRTH) call FDR(Npmax,maxnab, ! Calculate dissipative forces + rX,rY,vX,vY, + FdX,FdY,npoint, + list,Np,dLx,dLy, + dLxINV,dLyINV, + gamma) if (mod(it,incstp).eq.0) then call VCMTP(Npmax, + Np,vX,vY, ! Take samples of the center-of-mass + dNpINV,pm, ! velocity and the temperature of + pmINV,Tfac, ! the system + vCM,temp) write(3,1000) ! Write the results + it, + sngl( dble(it)*dt ), + vCM, sngl(temp) call FLUSH(3) call R_GYR(Npmax, ! Calculate the radius of gyration squared + Nclen,rX,rY, ! and the end-to-end distance squared + Np,N_water, + len_chain, + dLx,dLy, + dLxINV,dLyINV, + dr_gyr,dr_ee) write(21,755) it, ! Write the results + dr_gyr, dr_ee call FLUSH(21) call VACF(Npmax, ! Update the velocity autocorrelation + Np,N_water, ! functions that are being calculated + len_chain, + maxt_vacf,it_vacf, + ivacf_counter,it_bound,ivacf_loop, + vX,vY,velX,velY,dvacf) endif if (mod(it,issstp).eq.0) then nconf = nconf + 1 call snaps(Npmax, ! Take a snapshot of the polymer chain + Np,N_water, + Nclen,len_chain, + nconf,dLx,dLy, + dLxINV,dLyINV,rX,rY,Ifile) endif enddo cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c The dynamics in equilibrium has been completed. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !-------------------------------------- if(ivacf_counter. ! The purpose here is to calculate the + ge.maxt_vacf) ! velocity autocorrelation function and + then ! then write the final averaged result ! to a file !-------------------------------------- do it=0,(maxt_vacf-1) dvacf(it) = dvacf(it)/dfloat(ivacf_counter + - maxt_vacf + 1) end do dcs_vacf(0) = dvacf(0) * (1.d0/6.d0) * + dt * dfloat(incstp) do it=1,(maxt_vacf-1) dcs_vacf(it) = dcs_vacf(it-1) + + dvacf(it) * (1.d0/3.d0) * dt * dfloat(incstp) end do do it=0,(maxt_vacf-1) ! Write the results write(24,*) + sngl(dt*(it*incstp)), + sngl(dvacf(it)), + sngl(dcs_vacf(it)) end do end if close(3) close(21) close(24) stop 'Groovy!' ! The DPD simulation is over. end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c Generate the initial configuration c (positions and velocities). c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine GENCON(Npmax,rX,rY,vX,vY,N_water,Np,len_chain, + dLx,dLy,dLxINV,dLyINV,dNpINV,pm,pmINV,Tfac) implicit none integer + Npmax, N_water, Np, len_chain double precision + rX(Npmax), rY(Npmax), + vX(Npmax), vY(Npmax), + dLx, dLy, dLxINV, dLyINV, + dNpINV, pm, pmINV, Tfac, + drX_tmp, drY_tmp, r_dist, d_scale, + rXi, rYi, vCMX, vCMY, vv, tscal, temper integer + n_flag(Npmax), + ip, ip2 real + R2S ! Generate random initial positions do ip = 1, N_water ! for the water particles rX(ip) = + (R2S()-0.5)*dLx rY(ip) = + (R2S()-0.5)*dLy enddo ip = N_water + 1 ! Generate random initial positions ! for the monomers of the chain rX(ip) = + (R2S()-0.5)*dLx rY(ip) = + (R2S()-0.5)*dLy do ip2 = 1,(len_chain - 1) drX_tmp = ( R2S() - 0.5 ) drY_tmp = ( R2S() - 0.5 ) r_dist = dsqrt( drX_tmp**2 + drY_tmp**2 ) d_scale = 0.5d0/r_dist rX(ip+ip2) = rX(ip+ip2-1) + (drX_tmp * d_scale) rY(ip+ip2) = rY(ip+ip2-1) + (drY_tmp * d_scale) end do do ip = 1, Np ! Make sure that all particles are ! inside the simulation box rXi = rX(ip) rYi = rY(ip) rX(ip) = rXi - dLx*dnint( rXi*dLxINV ) rY(ip) = rYi - dLy*dnint( rYi*dLyINV ) enddo vCMX = 0.d0 ! Then generate random initial vCMY = 0.d0 ! velocities for the particles do ip = 1, Np vv = dble( R2S() - 0.5 ) vX(ip) = vv vCMX = vCMX + pm*vv vv = dble( R2S() - 0.5 ) vY(ip) = vv vCMY = vCMY + pm*vv enddo vCMX = vCMX*dNpINV*pmINV ! Satisfy the condition that the vCMY = vCMY*dNpINV*pmINV ! center-of-mass velocity vanishes do ip = 1, Np vX(ip) = vX(ip) - vCMX vY(ip) = vY(ip) - vCMY enddo temper = 0.d0 ! Finally rescale velocities such that do ip = 1, Np ! the kinetic temperature is equal to one. temper = temper + + pm*( vX(ip)*vX(ip) + + vY(ip)*vY(ip) ) enddo temper = Tfac*temper tscal = 1.d0/dsqrt( temper ) do ip = 1, Np vX(ip) = vX(ip)*tscal vY(ip) = vY(ip)*tscal enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c Shift (fold) particles back into primary simulation box. c c NOTE: Calls of FOLDB are allowed only when the Verlet list c has to be reconstructed. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine FOLDB(Npmax,Np,rX,rY,dLx,dLy,dLxINV,dLyINV) implicit none integer + Npmax, Np, ip double precision + rX(Npmax), rY(Npmax), + dLx, dLy, dLxINV, dLyINV, + rXi, rYi do ip = 1, Np rXi = rX(ip) rYi = rY(ip) rX(ip) = rXi - dLx*dnint( rXi*dLxINV ) rY(ip) = rYi - dLy*dnint( rYi*dLyINV ) enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c Integrate velocities using a time step of "hcd" for c conservative and dissipative forces, and a time step c of "hr" for random forces. The routine integrates c velocities over a time increment of dt/2, and therefore c integrates velocities in two parts. This approach is based c on the DPD-VV integration scheme (see the reference: c [G. Besold, I. Vattulainen, M. Karttunen, and J.M. Polson, c Phys. Rev. E Rapid Comm. vol. 62, R7611 (2000)]). c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine INTV (Npmax,Np,vX,vY, + FcX,FcY,FdX,FdY,FrX,FrY,pmINV,hcd,hr) implicit none integer + Npmax, Np, ip double precision + vX(Npmax), vY(Npmax), + FcX(Npmax), FcY(Npmax), + FdX(Npmax), FdY(Npmax), + FrX(Npmax), FrY(Npmax), + pmINV, hcd, hr do ip = 1, Np vX(ip) = vX(ip) + pmINV*( (FcX(ip)+FdX(ip))*hcd + FrX(ip)*hr ) vY(ip) = vY(ip) + pmINV*( (FcY(ip)+FdY(ip))*hcd + FrY(ip)*hr ) enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c Integrate the positions of the particles. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine INTR(Npmax,Np,rX,rY,vX,vY,h) implicit none integer + Npmax, Np, ip double precision + rX(Npmax), rY(Npmax), + vX(Npmax), vY(Npmax), h do ip = 1, Np rX(ip) = rX(ip) + vX(ip)*h rY(ip) = rY(ip) + vY(ip)*h enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c Calculate all forces: conservative, dissipative, and random c forces. c c If needed (the logical variable update = .true.), c then update the Verlet neighbor table. Otherwise use c the present table. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine FCDR(Npmax,maxnab,rX,rY,rX0,rY0,vX,vY, + FcX,FcY,FdX,FdY,FrX,FrY,npoint,list,n_flag,alpha, + Np,N_water,len_chain,dLx,dLy,dLxINV,dLyINV,rfac, + update,nVintv,nVsize,gamma,sigma,d_bond,skinSQ) implicit none integer + Npmax, maxnab, + npoint(Npmax), list(maxnab), n_flag(Npmax), + Np, N_water, len_chain, nVintv, nVsize double precision + rX(Npmax), rY(Npmax), + rX0(Npmax), rY0(Npmax), + vX(Npmax), vY(Npmax), + FcX(Npmax), FcY(Npmax), + FdX(Npmax), FdY(Npmax), + FrX(Npmax), FrY(Npmax), + alpha(2,2), rfac, d_bond, + dLx, dLy, dLxINV, dLyINV, + gamma, sigma, skinSQ logical + update integer + ip, jp, nlist, nfi, nfj, + jpbeg, jpend, jpnab double precision + rXi, rYi, vXi, vYi, + FcXi, FcYi, FdXi, FdYi, FrXi, FrYi, + rXij, rYij, rijSQ, rij, rijINV, omega, + eXij, eYij, vXij, vYij, Fcfac, Fdfac, Frfac, + FcXij, FcYij, FdXij, FdYij, FrXij, FrYij real + R2S do ip = 1, Np ! Initialize all forces FcX(ip) = 0.d0 FcY(ip) = 0.d0 FdX(ip) = 0.d0 FdY(ip) = 0.d0 FrX(ip) = 0.d0 FrY(ip) = 0.d0 enddo if ( update ) then ! Check if the Verlet list should ! be updated nVintv = 1 do ip = 1, Np ! Ok, let us do it. rX0(ip) = rX(ip) ! Save the particle positions. rY0(ip) = rY(ip) enddo nlist = 0 do ip = 1, Np-1 ! Update the Verlet list and calculate ! all forces npoint(ip) = nlist + 1 rXi = rX(ip) ! Position of particle "ip" (x coordinate) rYi = rY(ip) vXi = vX(ip) ! Velocity of particle "ip" vYi = vY(ip) FcXi = FcX(ip) ! Conservative forces acting on particle "ip" FcYi = FcY(ip) FdXi = FdX(ip) ! Dissipative forces FdYi = FdY(ip) FrXi = FrX(ip) ! Random forces FrYi = FrY(ip) nfi = n_flag(ip) ! Type of particle "ip" do jp = (ip+1), Np ! Go over all pairs within the Verlet list rXij = rXi - rX(jp) ! Vector from "jp" to "ip" rYij = rYi - rY(jp) rXij = rXij - dLx*dnint( rXij*dLxINV ) rYij = rYij - dLy*dnint( rYij*dLyINV ) rijSQ = ! Distance between two particles squared + rXij*rXij + + rYij*rYij if ( rijSQ .lt. skinSQ ) then ! The particles are neighbors nlist = nlist + 1 list(nlist) = jp if ( nlist .eq. maxnab ) stop 'list too small' if ( rijSQ .lt. 1.d0 ) then ! The two particles interact rij = dsqrt( rijSQ ) ! Distance between the particles rijINV = 1.d0/rij omega = 1.d0 - rij eXij = rXij*rijINV eYij = rYij*rijINV vXij = vXi - vX(jp) ! Relative velocity vYij = vYi - vY(jp) nfj = n_flag(jp) ! Type of particle "jp" ('1' or '2') Fcfac = omega ! Calculate the pairwise concervative force FcXij = + Fcfac * eXij * alpha(nfi,nfj) FcYij = + Fcfac * eYij * alpha(nfi,nfj) Fdfac = ! Calculate the pairwise dissipative force + omega*omega*( eXij*vXij + + eYij*vYij ) FdXij = Fdfac*eXij FdYij = Fdfac*eYij Frfac = ! Calculate the pairwise random force + omega * rfac * ( 2.*R2S() - 1. ) FrXij = Frfac*eXij FrYij = Frfac*eYij FcXi = FcXi + FcXij ! Update forces FcYi = FcYi + FcYij FdXi = FdXi + FdXij FdYi = FdYi + FdYij FrXi = FrXi + FrXij FrYi = FrYi + FrYij FcX(jp) = FcX(jp) - FcXij FcY(jp) = FcY(jp) - FcYij FdX(jp) = FdX(jp) - FdXij FdY(jp) = FdY(jp) - FdYij FrX(jp) = FrX(jp) - FrXij FrY(jp) = FrY(jp) - FrYij endif endif enddo FcX(ip) = FcXi FcY(ip) = FcYi FdX(ip) = FdXi FdY(ip) = FdYi FrX(ip) = FrXi FrY(ip) = FrYi enddo npoint(Np) = nlist + 1 nVsize = nlist + 1 else ! There is no need to update the Verlet table, ! so let us use the current list and calculate nVintv = nVintv + 1 ! all forces do ip = 1, Np-1 jpbeg = npoint(ip) jpend = npoint(ip+1) - 1 if (jpbeg.le.jpend) then rXi = rX(ip) rYi = rY(ip) vXi = vX(ip) vYi = vY(ip) FcXi = FcX(ip) FcYi = FcY(ip) FdXi = FdX(ip) FdYi = FdY(ip) FrXi = FrX(ip) FrYi = FrY(ip) nfi = n_flag(ip) do jpnab = jpbeg, jpend jp = list(jpnab) rXij = rXi - rX(jp) rYij = rYi - rY(jp) rXij = rXij - dLx*dnint( rXij*dLxINV ) rYij = rYij - dLy*dnint( rYij*dLyINV ) rijSQ = rXij*rXij + rYij*rYij if (rijSQ.lt.1.d0) then rij = dsqrt( rijSQ ) rijINV = 1.d0/rij omega = 1.d0 - rij eXij = rXij*rijINV eYij = rYij*rijINV vXij = vXi - vX(jp) vYij = vYi - vY(jp) nfj = n_flag(jp) Fcfac = omega FcXij = Fcfac * eXij * alpha(nfi,nfj) FcYij = Fcfac * eYij * alpha(nfi,nfj) Fdfac = omega*omega*( eXij*vXij + + eYij*vYij ) FdXij = Fdfac*eXij FdYij = Fdfac*eYij Frfac = omega * rfac*( 2.*R2S() - 1. ) FrXij = Frfac*eXij FrYij = Frfac*eYij FcXi = FcXi + FcXij FcYi = FcYi + FcYij FdXi = FdXi + FdXij FdYi = FdYi + FdYij FrXi = FrXi + FrXij FrYi = FrYi + FrYij FcX(jp) = FcX(jp) - FcXij FcY(jp) = FcY(jp) - FcYij FdX(jp) = FdX(jp) - FdXij FdY(jp) = FdY(jp) - FdYij FrX(jp) = FrX(jp) - FrXij FrY(jp) = FrY(jp) - FrYij endif enddo FcX(ip) = FcXi FcY(ip) = FcYi FdX(ip) = FdXi FdY(ip) = FdYi FrX(ip) = FrXi FrY(ip) = FrYi endif enddo endif do ip = 1, Np ! Multiply remaining forces with prefactors FdX(ip) = -gamma*FdX(ip) FdY(ip) = -gamma*FdY(ip) FrX(ip) = sigma*FrX(ip) FrY(ip) = sigma*FrY(ip) enddo ip = N_water + 1 ! Calculate conservative forces due to ! harmonic springs between adjacent monomers do jp = ip,(ip + len_chain - 2) rXij = rX(jp) - rX(jp + 1) rYij = rY(jp) - rY(jp + 1) rXij = rXij - dLx*dnint( rXij*dLxINV ) rYij = rYij - dLy*dnint( rYij*dLyINV ) FcX(jp) = FcX(jp) - d_bond * rXij FcY(jp) = FcY(jp) - d_bond * rYij FcX(jp+1) = FcX(jp+1) + d_bond * rXij FcY(jp+1) = FcY(jp+1) + d_bond * rYij end do return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c Calculate dissipative forces only. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine FDR(Npmax,maxnab,rX,rY,vX,vY,FdX,FdY, + npoint,list,Np,dLx,dLy,dLxINV,dLyINV,gamma) implicit none integer + Npmax, maxnab, + npoint(Npmax), list(maxnab), Np double precision + rX(Npmax), rY(Npmax), + vX(Npmax), vY(Npmax), + FdX(Npmax), FdY(Npmax), + dLx, dLy, dLxINV, dLyINV, gamma integer + ip, jp, + jpbeg, jpend, jpnab double precision + rXi, rYi, vXi, vYi, + FdXi, FdYi, FdXij, FdYij, + rXij, rYij, rijSQ, rij, rijINV, omega, + eXij, eYij, vXij, vYij, Fdfac do ip = 1, Np FdX(ip) = 0.d0 FdY(ip) = 0.d0 enddo do ip = 1, Np-1 jpbeg = npoint(ip) jpend = npoint(ip+1) - 1 if ( jpbeg .le. jpend ) then rXi = rX(ip) rYi = rY(ip) vXi = vX(ip) vYi = vY(ip) FdXi = FdX(ip) FdYi = FdY(ip) do jpnab = jpbeg, jpend jp = list(jpnab) rXij = rXi - rX(jp) rYij = rYi - rY(jp) rXij = rXij - dLx*dnint( rXij*dLxINV ) rYij = rYij - dLy*dnint( rYij*dLyINV ) rijSQ = rXij*rXij + rYij*rYij if ( rijSQ .lt. 1.d0 ) then rij = dsqrt( rijSQ ) rijINV = 1.d0/rij omega = 1.d0 - rij eXij = rXij*rijINV eYij = rYij*rijINV vXij = vXi - vX(jp) vYij = vYi - vY(jp) Fdfac = omega*omega*( eXij*vXij + + eYij*vYij ) FdXij = Fdfac*eXij FdYij = Fdfac*eYij FdXi = FdXi + FdXij FdYi = FdYi + FdYij FdX(jp) = FdX(jp) - FdXij FdY(jp) = FdY(jp) - FdYij endif enddo FdX(ip) = FdXi FdY(ip) = FdYi endif enddo do ip = 1, Np ! Multiply forces with prefactors FdX(ip) = -gamma*FdX(ip) FdY(ip) = -gamma*FdY(ip) enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c Check if the Verlet neighbor table should be updated. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine CHECK(Npmax,Np,rX,rY,rX0,rY0,skin,update) implicit none integer + Npmax, Np, ip double precision + rX(Npmax), rY(Npmax), + rX0(Npmax), rY0(Npmax), + skin, dispmx logical + update dispmx = 0.d0 do ip=1, Np dispmx = dmax1( dabs( rX(ip) - rX0(ip) ), dispmx ) dispmx = dmax1( dabs( rY(ip) - rY0(ip) ), dispmx ) enddo dispmx = 2.d0*dsqrt( 3.d0*dispmx*dispmx ) update = (dispmx.gt.(skin-1.d0)) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c Calculate the center-of-mass velocity of the system c (to guarantee that momentum conservation is satisfied), c and the instantaneous kinetic temperature of the system c based on the velocoties of the particles. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine VCMTP(Npmax,Np,vX,vY, + dNpINV,pm,pmINV,Tfac,vCM,temp) implicit none integer + Npmax, Np, ip double precision + vX(Npmax), vY(Npmax), + vCM, dNpINV, pm, pmINV, Tfac, temp, + vXi, vYi, vCMX, vCMY vCMX = 0.d0 vCMY = 0.d0 temp = 0.d0 do ip = 1, Np vXi = vX(ip) vYi = vY(ip) vCMX = vCMX + pm*vXi vCMY = vCMY + pm*vYi temp = temp + pm*( vXi*vXi + vYi*vYi ) enddo vCM = dNpINV*pmINV*dsqrt(vCMX*vCMX + vCMY*vCMY ) temp = Tfac*temp return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c Calculate c - the radius of gyration squared c - the end-to-end distance squared c of the polymer chain. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine R_GYR(Npmax,Nclen,rX,rY, + Np,N_water,len_chain,dLx,dLy,dLxINV,dLyINV, + dr_gyr,dr_ee) implicit none integer + Npmax, Nclen, Np, N_water, len_chain double precision + rX(Npmax), rY(Npmax), + dLx, dLy, dLxINV, dLyINV, dr_gyr, dr_ee, + drX_dpl(Nclen), drY_dpl(Nclen), + drX_cm, drY_cm, drX_g, drY_g, + disp_X, disp_Y, drX_tm, drY_tm, + drX_ee, drY_ee integer + j, lenc drX_cm = 0.d0 drY_cm = 0.d0 drX_g = 0.d0 drY_g = 0.d0 drX_dpl(1) = 0.d0 drY_dpl(1) = 0.d0 ! Determine the positions of the ! monomers in a chain with respect do j = 2,len_chain ! to the first monomer disp_X = rX(N_water+j) - rX(N_water+j-1) disp_Y = rY(N_water+j) - rY(N_water+j-1) disp_X = disp_X - dLx*dnint( disp_X * dLxINV ) disp_Y = disp_Y - dLy*dnint( disp_Y * dLyINV ) drX_dpl(j) = drX_dpl(j-1) + disp_X drY_dpl(j) = drY_dpl(j-1) + disp_Y end do ! Then calculate the radius of gyration do j = 1,len_chain ! over the whole chain drX_cm = drX_cm + drX_dpl(j) drY_cm = drY_cm + drY_dpl(j) end do drX_cm = drX_cm / len_chain drY_cm = drY_cm / len_chain do j = 1,len_chain drX_tm = drX_dpl(j) - drX_cm drY_tm = drY_dpl(j) - drY_cm drX_tm = drX_tm - dLx*dnint( drX_tm * dLxINV ) drY_tm = drY_tm - dLy*dnint( drY_tm * dLyINV ) drX_g = drX_g + drX_tm * drX_tm drY_g = drY_g + drY_tm * drY_tm end do drX_g = drX_g / len_chain drY_g = drY_g / len_chain dr_gyr = drX_g + drY_g ! Radius of gyration squared (whole chain) lenc = len_chain ! Then calculate the squared drX_ee = drX_dpl(lenc) ! end-to-end distance drY_ee = drY_dpl(lenc) drX_ee = drX_ee * drX_ee drY_ee = drY_ee * drY_ee dr_ee = drX_ee + drY_ee ! Squared end-to-end distance return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c Take a snapshot of the current configuration of c the polymer chain. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine snaps(Npmax,Np,N_water,Nclen,len_chain,nconf, + dLx,dLy,dLxINV,dLyINV,rX,rY,Ifile) implicit none integer + Npmax, Np, N_water, Nclen, len_chain, nconf, + ip, j, Ifile double precision + rX(Npmax), rY(Npmax), + drX_dpl(Nclen), drY_dpl(Nclen), + disp_X, disp_Y, dLx, dLy, dLxINV, dLyINV ip = N_water + 1 drX_dpl(1) = rX(ip) ! The starting position is the position drY_dpl(1) = rY(ip) ! of one of the end monomers do j = 1,(len_chain - 1) disp_X = rX(ip + j) - rX(ip + j - 1) disp_Y = rY(ip + j) - rY(ip + j - 1) disp_X = disp_X - dLx*dnint( disp_X * dLxINV ) disp_Y = disp_Y - dLy*dnint( disp_Y * dLyINV ) drX_dpl(j+1) = ! The positions of other monomers are + drX_dpl(j) + ! written with respect to the starting + disp_X ! position. Thus it is possible that drY_dpl(j+1) = ! some of the monomers reside outside + drY_dpl(j) + ! the simulation box. + disp_Y end do 925 format(1x,1f12.6,1x,1f12.6) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Here we write the (x,y)-coordinates of the monomers c of the polymer chain. The files will be called "fort.abc" c where "abc" is a number ranging from Ifile to (Ifile + 49). c The parameter Ifile is defined to have a value of 50, thus c the maximum number of snapshots is also 50. c c The if-loop below is not needed on most computer architectures. c However, on the Linux-workstations used during the summer c school, the Fortran77 compilers (f77 and g77) limit the c maximum unit number to be 99 (fort.99 is ok, while fort.100 c leads to an error message and core dumped). if(nconf.lt.50) then do j=1,len_chain write(nconf+Ifile,925) drX_dpl(j), drY_dpl(j) end do close(nconf+Ifile) endif return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c Update the velocity autocorrelation function for the c study of the tracer diffusion coefficient to characterize c the motion of the polymer chain along the 2D plane. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine VACF(Npmax,Np,N_water,len_chain,maxt_vacf, + it_vacf,ivacf_counter,it_bound,ivacf_loop, + vX,vY,velX,velY,dvacf) implicit none integer + Npmax, Np, N_water, len_chain, + maxt_vacf, it_vacf, + it_bound(maxt_vacf + 1), + ivacf_loop(-maxt_vacf:maxt_vacf), + ivacf_counter, ihist, it_diff, ip double precision + vX(Npmax), vY(Npmax), + velX(maxt_vacf), velY(maxt_vacf), + dvacf(0:(maxt_vacf-1)) it_vacf = ! Rotate the label of the register + it_bound(it_vacf + 1) velX(it_vacf) = 0.d0 ! Calculate the velocity of the velY(it_vacf) = 0.d0 ! center-of-mass of the polymer chain do ip=(N_water+1),Np velX(it_vacf) = velX(it_vacf) + vX(ip) velY(it_vacf) = velY(it_vacf) + vY(ip) end do velX(it_vacf) = velX(it_vacf) / len_chain velY(it_vacf) = velY(it_vacf) / len_chain ivacf_counter = ivacf_counter + 1 if(ivacf_counter.ge.maxt_vacf) then do it_diff=0,(maxt_vacf - 1) ihist = ivacf_loop( it_vacf - it_diff ) dvacf(it_diff) = dvacf(it_diff) + + velX(it_vacf) * + velX(ihist) + ! Update the velocity autocorrelation function + velY(it_vacf) * + velY(ihist) end do end if return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE: c c The pseudorandom number generator. c This subroutine has been taken as it is. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c*********************************************************************** real function R2S() c*************************************************************(GB 12/95) c c Portable long-period (ca. 2.3 * 10^18) random number generator of c L'ECUYER [1] with BAYS-DURHAM shuffle [2] and added safeguards as c published by PRESS et al. [3] as "ran2" in 1992. In this version c (called "R2S" for distinction) no argument is needed, and the c initialization can be done by an entry point "R2INIS(iseedS)" with c any positive integer "iseedS" as seed. The internal state corres- c ponding to iseedS = 1 is used if no initialization with "R2INIS" c is done before the first call of "R2S". c c "R2S" returns a uniform random number in the interval ]0., 1.[ c (NOTE: The endpoint values are excluded!) c c c * INITIALIZATION of "R2S": rdummy = R2INIS(iseedS) c c * GENERATION of a random number: r = R2S() c c c NOTE: c c * "rdummy" is a dummy variable of type REAL. c c * No variable declaractions are necessary in the calling module. c c * Parameter RNMX=1.-EPS should approximate the largest floating c point value that is less than 1. (EPS for a specific machine c can be determined with subroutine "MACHAR" from chapt. 20.1 of c ref [3]). (For "IBM RISC System/6000" workstations with "AIX XL c Fortran/6000", subroutine MACHAR gives EPS=1.192092896E-07 ). c c c REFERENCES: c c [1] P. L'ECUYER, Communications of the ACM, vol. 31 (1988) 742-774. c [2] in D.E. KNUTH, "Seminumerical Algorithms" (2nd. ed.), vol. 2 of c "The Art of Computer Programming", Addison-Wesley, Reading, MA c (1981) paragraphs 3.2-3.3 . c [3] W.H. PRESS, S.A. TEUKOLSKY, W.T. VETTERLING, and B.P. FLANNERY, c "Numerical Recipes in FORTRAN" (2nd ed.), Cambridge University c Press, Cambridge (1992), chapt. 7.1 . c c c TEST OUTPUT (first 35 numbers for iseed=1, in row-wise sequence): c c 0.285381 0.253358 0.093469 0.608497 0.903420 0.195873 0.462954 c 0.939021 0.127216 0.415931 0.533739 0.107446 0.399671 0.650371 c 0.027072 0.076975 0.068986 0.851946 0.637346 0.573097 0.902278 c 0.887676 0.372177 0.347516 0.207896 0.569131 0.364677 0.392418 c 0.366707 0.432149 0.153942 0.626891 0.749454 0.341041 0.830393 c c*********************************************************************** parameter ( IM1 = 2147483563, * IM2 = 2147483399, * AM = 1./IM1, * IMM1 = IM1-1, * IA1 = 40014, * IA2 = 40692, * IQ1 = 53668, * IQ2 = 52774, * IR1 = 12211, * IR2 = 3791, * NTAB = 32, * NDIV = 1+IMM1/NTAB, * EPS = 1.2e-7, * RNMX = 1.-EPS ) integer ivS(NTAB) save ivS, iyS, idum1S, idum2S data idum1S/1720212868/, idum2S/1/, iyS/1720212868/ data ivS/ 1720212868, 1392842846, 1031324961, 718590712, * 82237802, 1816996195, 1529538438, 1789446856, * 156648835, 52437849, 1441478319, 36906150, * 1269685686, 1644535938, 394503142, 310212663, * 1596049480, 7553450, 322224693, 445508654, * 28884682, 643161691, 407948861, 479214492, * 2124954851, 612891482, 112933431, 1814689225, * 53445315, 1904850491, 1695805043, 1860990862 / c*********************************************************************** c c---- Compute "idum1S" by idum1S = mod( IA1*idum1S, IM1 ) without c overflow by SCHRAGE's method (see ref. [3]): kS = idum1S / IQ1 idum1S = IA1*( idum1S - kS*IQ1 ) - kS*IR1 if ( idum1S .lt. 0 ) idum1S = idum1S + IM1 c---- Compute "idum2S" likewise: kS = idum2S / IQ2 idum2S = IA2*( idum2S - kS*IQ2 ) - kS*IR2 if ( idum2S .lt. 0 ) idum2S = idum2S + IM2 c---- "jS" will be in the range [1 (1) NTAB] : jS = 1 + iyS/NDIV c---- Here "idum1S" is shuffled, and "idum1S" and "idum2S" are combined c to produce output: iyS = ivS(jS) - idum2S ivS(jS) = idum1S if ( iyS .lt. 1 ) iyS = iyS + IMM1 c---- Because users don't expect endpoint values: R2S = min( AM*iyS, RNMX ) return c>>>> Initialization: >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> entry R2INIS (iseedS) c---- Be sure to prevent a negative "iseedS" or iseedS = 0 : idum1S = max( iabs(iseedS), 1 ) idum2S = idum1S c---- Load the shuffle table (after 8 warm-ups): do 10 jS= NTAB+8, 1, -1 kS = idum1S / IQ1 idum1S = IA1*( idum1S - kS*IQ1 ) - kS*IR1 if ( idum1S .lt. 0 ) idum1S = idum1S + IM1 if ( jS .le. NTAB ) ivS(jS) = idum1S 10 continue iyS = ivS(1) R2INIS = iyS return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc