随机数生成无变化
大家好,这是我写的一段小程序,用以在一个四方形的盒子里生成一定数量的结构,这些结构彼此相同,但在盒子里生成的位置是随机的,这些结构的朝向也是随机的(也就是在三维空间中,随机的旋转了一定的角度)。
程序如下。
SUBROUTINE MUL(cc,rx,ry,rz)
!分别沿x,y和z轴将结构先后随机地算转一定的角度。
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,14)
INTEGER :: i
REAL (KIND=dp) :: rx(3,3), ry(3,3), rz(3,3), cc(2,3) !rotation matrix and intermediate matrix
REAL (KIND=dp) :: pi, ag(3), rn !constant pi; angle around x, y and z axes; random number
!Constant
pi = DACOS(-1.0d0)
!
DO i = 1, 3, 1
CALL RANDOM_SEED()
CALL RANDOM_NUMBER (rn)
ag(i) = rn*2.0d0*pi
END DO
rx = 0.0d0
rx(1,1) = 1.0d0
rx(2,1) = DCOS(ag(1))
rx(2,3) = 0.0d0-DSIN(ag(1))
rx(3,2) = DSIN(ag(1))
rx(3,3) = DCOS(ag(1))
ry = 0.0d0
ry(1,1) = DCOS(ag(2))
ry(1,3) = DSIN(ag(2))
ry(2,2) = 1.0d0
ry(3,1) = 0.0d0-DSIN(ag(2))
ry(3,3) = DCOS(ag(2))
rz = 0.0d0
rz(1,1) = DCOS(ag(3))
rz(1,2) = 0.0d0-DSIN(ag(3))
rz(2,1) = DSIN(ag(3))
rz(2,2) = DCOS(ag(3))
rz(3,3) = 1.0d0
cc(2,1) = rx(1,1)*cc(1,1)+rx(1,2)*cc(1,2)+rx(1,3)*cc(1,3)
cc(2,2) = rx(2,1)*cc(1,1)+rx(2,2)*cc(1,2)+rx(2,3)*cc(1,3)
cc(2,3) = rx(3,1)*cc(1,1)+rx(3,2)*cc(1,2)+rx(3,3)*cc(1,3)
cc(1,1) = ry(1,1)*cc(2,1)+ry(1,2)*cc(2,2)+ry(1,3)*cc(2,3)
cc(1,2) = ry(2,1)*cc(2,1)+ry(2,2)*cc(2,2)+ry(2,3)*cc(2,3)
cc(1,3) = ry(3,1)*cc(2,1)+ry(3,2)*cc(2,2)+ry(3,3)*cc(2,3)
cc(2,1) = rz(1,1)*cc(1,1)+rz(1,2)*cc(1,2)+rz(1,3)*cc(1,3)
cc(2,2) = rz(2,1)*cc(1,1)+rz(2,2)*cc(1,2)+rz(2,3)*cc(1,3)
cc(2,3) = rz(3,1)*cc(1,1)+rz(3,2)*cc(1,2)+rz(3,3)*cc(1,3)
RETURN
END SUBROUTINE MUL
PROGRAM GENERATE_STRUCTURE
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,14)
INTEGER :: i, j, k, m
REAL (KIND=dp) :: lc(3,3)
REAL (KIND=dp) :: sc(8,3) !coordinate of atoms in the unit cell
REAL (KIND=dp) :: nc(8,3) !coordinate of atoms after the rotation in the unit cell
REAL (KIND=dp), ALLOCATABLE :: gc(:,:,:) !coordinate of atoms after movement in the unit cell
REAL (KIND=dp) :: pi, ag(3), rn !constant pi; angle around x, y and z axes; random number
REAL (KIND=dp) :: rx(3,3), ry(3,3), rz(3,3), cc(2,3) !rotation matrix and intermediate matrix
INTEGER :: np !number of generation points
REAL (KIND=dp) :: di(3) !distance between each atom
REAL (KIND=dp) :: gp(3) !coordinate of the generation points
!Number of generation points
np = 3
!
!Structure coordinates初始结构里的8个原子的实数坐标值
sc(1,1) = 0.0d0
sc(1,2) = 0.0d0
sc(1,3) = 0.0d0
sc(2,1) = 2.84056201420281d0
sc(2,2) = 0.0d0
sc(2,3) = 0.0d0
sc(3,1) = 7.1014051420281d-1
sc(3,2) = 1.2299995d0
sc(3,3) = 0.0d0
sc(4,1) = 2.1304215d0
sc(4,2) = 1.2299995d0
sc(4,3) = 0.0d0
sc(5,1) = 0.0d0
sc(5,2) = 2.459999d0
sc(5,3) = 0.0d0
sc(6,1) = 2.84056201420281d0
sc(6,2) = 2.459999d0
sc(6,3) = 0.0d0
sc(7,1) = 7.1014051420281d-1
sc(7,2) = 3.6899985d0
sc(7,3) = 0.0d0
sc(8,1) = 2.1304215d0
sc(8,2) = 3.6899985d0
sc(8,3) = 0.0d0
!
OPEN (UNIT=3, FILE='POSCAR', STATUS='UNKNOWN')
!Box lattice盒子的大小
lc = 0.0d0
lc(1,1) = 24.7284d0
lc(2,2) = 23.795d0
lc(3,3) = 25.0d0
!
ALLOCATE (gc(np,8,3))
gc = 0.0d0
print*,'here1'
DO i = 1, np, 1
!将结构旋转一个随机的角度
1000 DO j = 1, 8, 1
cc(1,:) = sc(j,:)
CALL MUL(cc,rx,ry,rz)
nc(j,:) = cc(2,:)
END DO
!将结构放在盒子里的随机位置上
CALL RANDOM_SEED()
CALL RANDOM_NUMBER (rn)
gp(1) = rn*lc(1,1)
CALL RANDOM_SEED()
CALL RANDOM_NUMBER (rn)
gp(2) = rn*lc(2,2)
CALL RANDOM_SEED()
CALL RANDOM_NUMBER (rn)
gp(3) = rn*lc(3,3)
!判断如果结构中有一个原子超出了盒子的范围,就跳回到1000标号的位置,重新生成结构的位置。
DO j = 1, 8, 1
gc(i,j,1) = nc(j,1)+gp(1)
gc(i,j,2) = nc(j,2)+gp(2)
gc(i,j,3) = nc(j,3)+gp(3)
print *,'here2',i,j,gc(i,j,:)
IF ((gc(i,j,1) >= lc(1,1)) .OR. (gc(i,j,2) >= lc(2,2)) .OR. (gc(i,j,3) >= lc(3,3))) GOTO 1000
END DO
print*,'here3'
! DO j = 1, i, 1
! DO k = 1, 8, 1
! DO m = 1, 8, 1
! di(1) = gc(i,k,1)-gc(i,m,1)
! di(2) = gc(i,k,2)-gc(i,m,2)
! di(3) = gc(i,k,3)-gc(i,m,3)
! di(1) = DSQRT(di(1)**2+di(3)**2+di(3)**2)
! IF (di(1) <= 1.0d-1) GOTO 1000
! END DO
! END DO
! END DO
print*,'here4'
END DO
print*,'here5'
WRITE (UNIT=3, FMT=*) 'structure'
WRITE (UNIT=3, FMT=*) 1.0
DO i = 1, 3, 1
WRITE (UNIT=3, FMT=*) lc(i,:)
END DO
WRITE (UNIT=3, FMT=*) 'C'
WRITE (UNIT=3, FMT=*) np*8
WRITE (UNIT=3, FMT=*) 'Cartesian'
DO i = 1, np, 1
DO j = 1, 8, 1
WRITE (UNIT=3, FMT=*) gc(i,j,:)
END DO
END DO
DEALLOCATE (gc)
CLOSE (UNIT=3)
STOP
END PROGRAM GENERATE_STRUCTURE
当我运行程序后,总是在检验断点出不断反复地输出下面的信息。
here2 1 1 24.668052567646274 23.736930446253826 24.938989752315429
here2 1 2 27.551495834038786 26.533275580222909 24.982543938606533
here2 1 1 24.668052567646274 23.736930446253826 24.938989752315429
here2 1 2 27.551495834038786 26.533275580222909 24.982543938606533
here2 1 1 24.668052567646274 23.736930446253826 24.938989752315429
here2 1 2 27.551495834038786 26.533275580222909 24.982543938606533
......
......
......
也就是说,程序并没有生成一个可以落在盒子里的结构。而且每次运行源随机数后,生成的结构,都是相同的。
请问是哪里出了问题呢?
能麻烦大家给些修改建议吗?
谢谢,盼复。
random_seed 在主程序开始调用一次,中间不要再调用
页:
[1]