版块导航: 正在加载中...

登录注册

应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜

北京石油化工学院2026年研究生招生接收调剂公告

返回列表

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 ]

回复此楼

» 猜你喜欢

一志愿211，化学学硕，310分，本科重点双非，求调剂已经有13人回复
285求调剂已经有8人回复
22408 266求调剂已经有6人回复
264求调剂已经有8人回复
材料工程专业日语生求调剂已经有9人回复
297求调剂已经有15人回复
277、学硕，求调剂数一104，已经有12人回复
机械专硕273请求调剂已经有3人回复
080100力学316求调剂已经有3人回复
288环境专硕,求调材料方向已经有39人回复

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

查看全部散金贴

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

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

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