1

我正在尝试编写一个程序，使用最速下降算法最小化二维 400 原子系统的总能量。我的程序的总体思路如下：

  1. 获取原子坐标 (x, y)
  2. 随机选择一个原子
  3. 计算作用在该原子上的力的 x 和 y 分量
  4. 计算 x 和 y 位置的变化，dx 和 dy
  5. 生成新坐标（x+dx, y+dy）并更新数组
  6. 重复步骤 2-5，直到每个原子上的力变为 ~0

由于原子 210 上的初始力的大小很大，因此它是系统接近收敛的良好标志。我还没有修复代码，以便在力在某个容差范围内时停止迭代。话虽如此，我的代码在原子 210 上打印了力的 x 分量，以便我可以查看力是否趋于 0。

当我运行我的代码时，坐标数组似乎没有更新（上面的第 5 步）。我不确定是在这个网站还是物理网站上发布我的问题，但是，我相信我的问题涉及在 Fortran 77 中更新数组的一些技术性。如果这超出了本网站的范围，我很抱歉。我只是不知道该转向哪里。感谢大家提供的任何帮助。如果我的注释不清楚或者是否有人需要更多信息，请告诉我。这是代码。

  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
4

0 回答 0