[Fortran] 纯文本查看 复制代码
!
! PURPOSE: Entry point for the console application.
!
!****************************************************************************
module ini
implicit none
real*8,parameter::T=10.,wx=6.,wy=6.,wz=11., pi = 3.141592653589793239
integer,parameter::N=3
!T是温度，单位为微K，N是原子数，6，6，11刻画了原子的构型
contains
function ran() !returns random number between 0 - 1
implicit none
integer , save :: flag = 0
double precision :: ran
if(flag==0) then
call random_seed()
flag = 1
endif
call random_number(ran) ! built in fortran 90 random number function
end function ran
function normal(mean,sigma)
implicit none
integer :: flag
double precision, parameter :: pi = 3.141592653589793239
double precision :: u1, u2, y1, y2, normal, mean, sigma
save flag
data flag /0/
u1 = ran(); u2 = ran()
if (flag.eq.0) then
y1 = sqrt(-2.0d0*log(u1))*cos(2.0d0*pi*u2)
normal = mean + sigma*y1
flag = 1
else
y2 = sqrt(-2.0d0*log(u1))*sin(2.0d0*pi*u2)
normal = mean + sigma*y2
flag = 0
endif
end function normal
!The above codes are made in Fortran 90 language, if you have any question, you may write to sealin2008@hotmail.com
subroutine position0(R0)
implicit none
real*8::R0(N,3),x,y
integer ::i
do i=1,N
do
x=ran()
y=ran()
if((2*x-1)**2+(2*y-1)**2<=1) exit
end do
R0(i,1)=wx*(2.0*x-1.0)
R0(i,2)=wy*(2.0*y-1.0)
end do
do i=1,N
R0(i,3)=normal(0.D0,wz/2.0)
end do
end subroutine position0
subroutine Velocity(V)
implicit none
real*8::V(N,3),r1,r2,r3
integer::i
do i=1,N
r1=ran()*2.*pi
r2=ran()*2.-1.0
r3=normal(0.d0,dsqrt(1.3806504E4*T/(87*1.65053878)))
V(i,1)=r3*sqrt(1-r2*r2)*cos(r1)
V(i,2)=r3*sqrt(1-r2*r2)*sin(r1)
V(i,3)=r3*r2
end do
end subroutine Velocity
function position(i,t,R0,V)
implicit none
real*8::R0(N,3),V(N,3),t,position(3)
integer::i
position(1)=R0(i,1)+t*V(i,1)
position(2)=R0(i,2)+t*V(i,2)
position(3)=R0(i,3)+t*V(i,3)
end function position
end module ini
program main
use ini
implicit none
real*8::R0(N,3),V(N,3)
integer ::i
call position0(R0)
call Velocity(V)
Open( 12 , File = 'out.txt' )
Do i = 1 , N
write( 12 , * ) R0(i,1),R0(i,2),R0(i,3)
End Do
Do i = 1 , N
write( 12 , * ) V(i,1),V(i,2),V(i,3)
End Do
Close( 12 )
end program main
[Fortran] 纯文本查看 复制代码
R0 = position0()
V = Velocity()
Function position0() result(R0)
implicit none
real*8::R0(N,3),x,y
integer ::i
do i=1,N
do
x=ran()
y=ran()
if((2*x-1)**2+(2*y-1)**2<=1) exit
end do
R0(i,1)=wx*(2.0*x-1.0)
R0(i,2)=wy*(2.0*y-1.0)
end do
do i=1,N
R0(i,3)=normal(0.D0,wz/2.0)
end do
end Function position0
Function Velocity() result(V)
implicit none
real*8::V(N,3),r1,r2,r3
integer::i
do i=1,N
r1=ran()*2.*pi
r2=ran()*2.-1.0
r3=normal(0.d0,dsqrt(1.3806504E4*T/(87*1.65053878)))
V(i,1)=r3*sqrt(1-r2*r2)*cos(r1)
V(i,2)=r3*sqrt(1-r2*r2)*sin(r1)
V(i,3)=r3*r2
end do
end Function Velocity