Vicsek Model

PROGRAM VICSEK

! T VICSEK, PRL, 75, 6, 1226-1229, (1995)

! TO RUN THE PROGRAM:

! gfortran vicsek.f95; ./a.out

IMPLICIT NONE

INTEGER, PARAMETER :: N = 300

REAL*8 PI,LENGTH,LX,LY,V0,VX,VY,ETA,RCUT

REAL*8 DX,DY,R,ARG,AVGTHETA,NOISE,SUMX,SUMY,AVGVEL

REAL*8, DIMENSION (N) :: X,Y,THETA,NUM,DEN

INTEGER I,J,TIME,TMAX,TIMESTEP

INTEGER,PARAMETER :: SEED = 99999999

CALL SRAND(SEED)

PI = 4.D0*ATAN(1.D0)

OPEN(UNIT=5,FILE='System_information.dat',STATUS='UNKNOWN')

! USER INPUTS.

LENGTH = 5.D0

LX = LENGTH/2.D0

LY = LX

V0 = 0.03D0

ETA = 0.1D0

RCUT = 1.0D0

TMAX = 100

TIMESTEP = 1

! INITIAL POSITION & VELOCITY OF THE PARTICLES.

DO I = 1,N,1

  X(I) = 2.D0*(RAND()-0.5D0)*LX

  Y(I) = 2.D0*(RAND()-0.5D0)*LY

  THETA(I) = 2.D0*RAND()*PI

  VX = V0*DCOS(THETA(I))

  VY = V0*DSIN(THETA(I))

  WRITE(10,*) X(I),Y(I),THETA(I),VX,VY

ENDDO ! I

! TIME ITERATION

AVGVEL = 0.D0

DO TIME = 1,TMAX,TIMESTEP

  SUMX = 0.D0

  SUMY = 0.D0

  DO I = 1,N,1

    ! INCLUDE Ith PARTICLE IN ARCTAN(THETA)

    NUM(I) = DSIN(THETA(I))

    DEN(I) = DCOS(THETA(I))

      DO J = 1,N,1

        ! EVALUATE ARCTAN(THETA) OVER THE CIRCLE OF RADIUS RCUT

        ! MINIMUM IMAGE CONVENTION INTRODUCED

        DX = (X(I)-X(J)) - NINT((X(I)-X(J))/(2.D0*LX))*2.D0*LX

        DY = (Y(I)-Y(J)) - NINT((Y(I)-Y(J))/(2.D0*LY))*2.D0*LY

        R = DSQRT(DX**2 + DY**2)

          IF (R <= RCUT) THEN

            NUM(I) = NUM(I) + DSIN(THETA(J))

            DEN(I) = DEN(I) + DCOS(THETA(J))

          ENDIF

      ENDDO !J

    ! EVALUATE THETA

    ARG = ABS(NUM(I)/DEN(I))

      IF (NUM(I) >= 0.D0 .AND. DEN(I) >= 0.D0) THEN

        AVGTHETA = DATAN(ARG)

        ELSEIF (NUM(I) >= 0.D0 .AND. DEN(I) <= 0.D0) THEN

        AVGTHETA = PI - DATAN(ARG)

        ELSEIF (NUM(I) <= 0.D0 .AND. DEN(I) <= 0.D0) THEN

        AVGTHETA = PI + DATAN(ARG)

        ELSEIF (NUM(I) <= 0.D0 .AND. DEN(I) >= 0.D0) THEN

        AVGTHETA = 2.D0*PI - DATAN(ARG)

      ENDIF

    ! NOISE OR TEMPERATURE

    NOISE = ETA*(RAND()-0.5D0)

    ! EVALUATE FINAL THETA

    THETA(I) = AVGTHETA + NOISE

    ! TIME UPDATE

    VX = V0*DCOS(THETA(I))

    VY = V0*DSIN(THETA(I))

    X(I) = X(I) + VX*TIMESTEP

    Y(I) = Y(I) + VY*TIMESTEP

    ! PERIODIC BOUNDARY CONDITION

    X(I) = X(I) - 2.D0*(INT(X(I)/LX))*LX

    Y(I) = Y(I) - 2.D0*(INT(Y(I)/LY))*LY

  ENDDO ! I

  DO I = 1,N,1

    ! EVALUATE AVERAGED NORMALISED VELOCITY

    SUMX = SUMX + V0*DCOS(THETA(I))

    SUMY = SUMY + V0*DSIN(THETA(I))

    WRITE(TIME+100,*) X(I),Y(I),THETA(I),VX,VY

  ENDDO

  CLOSE (TIME+100)

  AVGVEL = AVGVEL + DSQRT(SUMX*SUMX+SUMY*SUMY)/(V0*DFLOAT(N))

ENDDO ! TIME

! EVALUATE ORDER PARAMETER

AVGVEL = AVGVEL/DFLOAT(TMAX)

WRITE(5,*) "Length of box=", LENGTH, "Noise/Temperature (Eta)=", ETA, "r=", RCUT, "Density (Rho)=", LX*LY/DFLOAT(N)

WRITE(5,*) "Average Normalised Velocity (V_a)=", AVGVEL

WRITE(5,*) "Total Timestep=", TMAX

 CLOSE(5)

END PROGRAM VICSEK