Trouble with the steepest descent algorithm

Question

I am trying to write a program that minimizes the total energy of a 2-dimensional, 400 atom system using the steepest descent algorithm. The general idea of my program is the following:

1. Get the atomic coordinates (x, y)
2. Randomly choose an atom
3. Compute the x- and y-component of the force on that atom
4. Compute the change in the x and y position, dx and dy
5. Generate new coordinates (x+dx, y+dy) and update array
6. Repeat steps 2-5 until the force on each atom becomes ~0

Since the magnitude of the initial force on atom 210 is large, it is a good marker of how close to the system is to convergence. I have not yet fixed the code so that the iterations cease when the force is within some tolerance factor. With that being said, my code prints the x-component of the force on atom 210 so that I can see if the force is tending towards 0.

When I run my code, it seems that the coordinate array is not updating (step 5 above). I wasn't exactly sure whether to post my question on this website or a physics website, however, I believe my issue involves some technicality for updating arrays in Fortran 77. I am sorry if this is outside of the scope of this website. I just wasn't sure where to turn. Thank you all for any help you are able to provide. Please let me know if my annotations are not clear or if anyone needs more information. Here is the code.

  program molstat_stpdesc
  implicit none

  !variables
  integer i, j, k, count
  real rmin, sigma, coord(3,400), rcut, fx, fy, drx
  real neighbors(400,24), nofneighbors(400), dry, lambda
  real forces(2,400,400), tforces(2,400)
  parameter (sigma=3.405)
  !i: iteration #
  !j: loop variable
  !k: random atom
  !count: lets me know if I need to update the map of neighbors
  !rmin: initial separation between the atoms (in angstroms)
  !sigma: Lennard-Jones sigma parameter
  !rcut: distance beyond which the potential energy of interaction is 0
  !lambda: factor multiplying the force to generate dr, usually small
  !forces: not used in this calculation
  !tforces: array containing x- and y-component of force on an atom

  open(9,file='minimization.out',status='replace')

  rmin=2**(1.0/6.0)*sigma
  rcut=7.5
  lambda=(1.0/5.0)*rmin

  !get the starting coordinates, coord
  !coord(1,i): number identifying atom i
  !coord(2,i): x-coordinate of atom i
  !coord(3,i): y-coordinate of atom i
  call construct(rmin,coord)

  !get the starting map of neighbors
  call map(rmin,rcut,coord1,neighbors,nofneighbors)

  count=0
  write(9,'("Forces on atom 210")')
  do i=1,1000 !# of iterations
     do j=1,400 !try to sample every atom per iteration
        k=int(400*rand(0)) !choose a random atom
        if (k .eq. 0) then !because random number generator may give me 0
           continue
        else
           !get the forces
           call force(coord,nofneighbors,neighbors,forces,tforces)
           fx=tforces(1,k)  !x-component of the force on atom j
           fy=tforces(2,k)  !y-component
           drx=lambda*fx    !step size in the x-direction
           dry=lambda*fy    !y-direction
           coord(2,k)=coord(2,k)+drx !step in the x-direction
           coord(3,k)=coord(3,k)+dry !y-direction
        endif
     enddo
     write(9,'("Iteration: ",I4)') i
     !print the force of an atom in the center of the block
     write(9,'("x-component: ",F7.4)') tforces(1,210)
     write(9,'()')
     count=count+1
     if (count .eq. 10) then
        !update map of neighbors
        call map(rmin,rcut,coord1,neighbors,nofneighbors)
        count=0
     else
        continue
     endif
  enddo

  close(9)

  open(10,file='minimization.xy',status='replace')
  write(10,'("#x-coordinate     y-coordinate")')
  write(10,'("#------------     ------------")')
  do i=1,400
     write(10,'("    ",F5.2,"            ",F5.2)')
 +        coord(2,i), coord(3,i)
  enddo

  end