Fortran Coder

查看: 183|回复: 1

[通用算法] 随机数生成无变化

[复制链接]

68

帖子

25

主题

0

精华

专家

F 币
321 元
贡献
197 点
发表于 2024-6-21 11:39:27 | 显示全部楼层 |阅读模式
大家好,

这是我写的一段小程序,用以在一个四方形的盒子里生成一定数量的结构,这些结构彼此相同,但在盒子里生成的位置是随机的,这些结构的朝向也是随机的(也就是在三维空间中,随机的旋转了一定的角度)。

程序如下。
[Fortran] 纯文本查看 复制代码
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

......
......
......
也就是说,程序并没有生成一个可以落在盒子里的结构。而且每次运行源随机数后,生成的结构,都是相同的。

请问是哪里出了问题呢?

能麻烦大家给些修改建议吗?

谢谢,盼复。

1995

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1481 元
贡献
636 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2024-6-21 13:05:29 | 显示全部楼层
random_seed 在主程序开始调用一次,中间不要再调用
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-7-23 19:13

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表