1

最急降下アルゴリズムを使用して、2 次元の 400 原子系の総エネルギーを最小化するプログラムを作成しようとしています。私のプログラムの一般的な考え方は次のとおりです。

  1. 原子座標 (x, y) を取得する
  2. ランダムに原子を選ぶ
  3. その原子にかかる力の x 成分と y 成分を計算します
  4. x と y の位置の変化、dx と dy を計算します。
  5. 新しい座標 (x+dx、y+dy) を生成し、配列を更新します
  6. 各原子にかかる力が ~0 になるまで、手順 2 ～ 5 を繰り返します。

原子 210 にかかる初期の力の大きさは大きいため、系がどれだけ収束に近づいているかを示す良い指標となります。力が許容範囲内にあるときに反復が停止するように、コードをまだ修正していません。そうは言っても、私のコードは原子 210 に対する力の x 成分を出力して、力が 0 に向かっているかどうかを確認できるようにします。

コードを実行すると、座標配列が更新されていないようです (上記の手順 5)。質問をこの Web サイトに投稿するか、物理学の Web サイトに投稿するかは正確にはわかりませんでしたが、Fortran 77 で配列を更新するための技術的な問題が含まれていると思います。これがこの Web サイトの範囲外である場合は申し訳ありません。どこを向いたらいいのかわからなかった。皆様のご協力に感謝いたします。私の注釈が明確でない場合、または誰かがより多くの情報を必要とする場合はお知らせください。これがコードです。

  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