vendredi 25 mars 2016

How to keep the distance between $n$ particles within a certain range?

So, I am working on a problem in Molecular Dynamics and I need to randomly generate a position array for np particles within a box of size [-L,L] x [-L,L]. In fact, I need to generate the x-array for the x-coordinates with x(1) = 0 and the y-array for the y-coordinates with y(1)=y(2) =0. I need the particles to be such that the distances between neighboring particles are within some range (e.g: 0.9 <= r <= 1.1) like in the following picture: enter image description here

However in my code I get something like this: enter image description here

See how the red lines are larger than what I want.

My code is

!Generating random coordinates for the particles
x(1) = 0.0d0
y(1) = 0.0d0
y(2) = 0.0d0
!-------------------------------------------------------------------------  
! translation and rotaion of the whole system were froze (saving 4 degrees of
! freedome)
! x(1) = 0.0d0; y(1) = 0.0d0    fix one particle in the origin
! y(2) = 0.0d0              fix the second particle on the x-axis
!-------------------------------------------------------------------------

rmatrix = 100.0
minv = 0.0
maxv = 10
iter0 = 0
101 DO WHILE(maxv >= 1.1 .OR. minv <= 0.9) 
iter0 = iter0 + 1
PRINT *, iter0
    CALL init_random_seed()
    DO i = 2, np
        CALL RANDOM_NUMBER(w1)
        x(i) = 10 * w1 - 5
    END DO
    DO i = 3, np
        CALL RANDOM_NUMBER(w2)
        y(i) = 10 * w2 - 5
    END DO
! rmatrix contains the distances between all particles 
    DO i = 1, np
        DO j = 1, np
            IF(j .NE. i) THEN
                xij = x(i) - x(j)
                yij = y(i) - y(j)
                rij = SQRT(xij * xij + yij * yij)
                rmatrix(i,j) = rij
            END IF
        END DO
    END DO
    minv = MINVAL(rmatrix)
    DO i = 1, np
        DO j = 1, np
            IF(j .NE. i) THEN
                maxv = MIN(maxv, rmatrix(i,j))
            END IF
        END DO
        IF(maxv >= 1.1) THEN
            GOTO 101
        END IF
    END DO
END DO

CONTANIS
      SUBROUTINE init_random_seed()
        INTEGER :: i, n, clock
        INTEGER, DIMENSION(:), ALLOCATABLE :: seed

        CALL RANDOM_SEED(size = n)
        ALLOCATE(seed(n))

        CALL SYSTEM_CLOCK(COUNT=clock)

        seed = clock + 37 * (/ (i - 1, i = 1, n) /)

        CALL RANDOM_SEED(PUT = seed)

      END SUBROUTINE init_random_seed




Aucun commentaire:

Enregistrer un commentaire