| 查看: 279 | 回复: 0 | |||
[交流]
【讨论】麻烦各位高手给看一下我这程序有问题吗？多谢
|
|
麻烦高手帮忙给看一下我这段程序有问题吗？为什么开始几十步的时候，能量基本保持不变，但是后来，动能急剧增加呢，势能也在增加，有劳各位了，多谢了 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 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 ] |
回复此楼» 猜你喜欢
售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急
已经有5人回复
-
找导师
已经有4人回复
-
我就是申请一个面上项目而已，这评审意见是按照杰青的条件评的吧？
已经有11人回复
-
有多少人是今天查系统知道结果的？
已经有18人回复
-
售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急
已经有4人回复
-
售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急
已经有4人回复
-
售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急
已经有4人回复
-
售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急
已经有4人回复
-
售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急
已经有5人回复
-
售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急
已经有4人回复
投票:
