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

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