24小时热门版块排行榜

返回列表

leigp

木虫 (正式写手)

应助: 10 (幼儿园)
金币: 1811.8
帖子: 638
在线: 225小时
虫号: 758557

[交流] 【讨论】麻烦各位高手给看一下我这程序有问题吗？多谢

麻烦高手帮忙给看一下我这段程序有问题吗？为什么开始几十步的时候，能量基本保持不变，但是后来，动能急剧增加呢，势能也在增加，有劳各位了，多谢了

SUBROUTINE L_J (n,a,ahalf,rcut,position_x,position_y,position_z,u,f_x,f_y,f_z,Ep)
IMPLICIT NONE
! the number of particles
INTEGER,INTENT(IN)::n
! the potential cutoff
REAL,INTENT(IN)::rcut
! the length and half length of box
REAL,INTENT(IN)::a,ahalf
!potential energy
REAL,INTENT(OUT)::Ep
! the position of particles
REAL(8),DIMENSION(n),INTENT(INOUT)::position_x,position_y,position_z
! the potential and force of particles
REAL,DIMENSION(n),INTENT(OUT)::u,f_x,f_y,f_z
! the contral variable of loop
INTEGER::i,j,k
! the label of particles
INTEGER::ijk
! the distance of every two particles along the three axis
REAL::distance_x,distance_y,distance_z
! the distance's square of every two particles
REAL::distance_sq
REAL::distance6,distance12,rcut_sq,rcut6,rcut12
REAL::Ec,ff
! initialize the potential and force
u=0.
Ep=0
f_x=0.
f_y=0.
f_z=0.
ijk=0
!calculate the rcut_sq
rcut_sq=rcut*rcut
rcut6=rcut_sq**3
rcut12=rcut6*rcut6
Ec=(1./rcut12-1./rcut6)
position_i: DO i=1,n-1
distance_ij:  DO j=i+1,n
distance_x=position_x(i)-position_x(j)
      distance_y=position_y(i)-position_y(j)
      distance_z=position_z(i)-position_z(j)
      !x_distance:
distance_x=distance_x-a*REAL(ANINT(distance_x/a))
      !y_distance:
distance_y=distance_y-a*REAL(ANINT(distance_y/a))
      !z_distance:
distance_z=distance_z-a*REAL(ANINT(distance_z/a))
      ! calculate the distance_sq
      distance_sq=distance_x**2+distance_y**2+distance_z**2
      ! calculate potential and force
   force:IF (distance_sq    ! calculate the distance6 and distance12
      distance6=distance_sq**3
      distance12=distance6*distance6
   ff=(distance_x/distance_sq)*((12./distance12)-(6./distance6))
   !  u(ijk)=u(ijk)+(1./distance12-1./distance6)-Ec
            f_x(i)=f_x(i)+4.*ff
            f_y(i)=f_y(i)+4.*ff
            f_z(i)=f_z(i)+4.*ff
            f_x(j)=f_x(j)-4.*ff
            f_y(j)=f_y(j)-4.*ff
            f_z(j)=f_z(j)-4.*ff
   Ep=Ep+(1./distance12-1./distance6)-Ec
   END IF force
   END DO distance_ij
  ! u(ijk)=4.*u(ijk)
END DO position_i
Ep=4.*Ep
Ep=Ep/real(n)
END SUBROUTINE L_J

[ Last edited by ghcacj on 2010-12-16 at 13:08 ]

回复此楼

» 猜你喜欢

关于核磁制样已经有3人回复
郑州大学申博已经有6人回复

» 抢金币啦！回帖就可以得到:

查看全部散金贴

1楼 2010-12-16 10:45:23

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

相关版块跳转我要订阅楼主 leigp 的主题更新