Statement functions which make use of system dependent features such as bit manipulation.
INTEGER*4 DUMARG, E
LOGICAL ODD
E( DUMARG ) = SIGN( 1, ISHFT( DUMARG, 31))
ODD( DUMARG ) = BTEST( DUMARG, 0 )
C set up a phase space point that corresponds to a face centred cubic
C crystal with random velocities constrained by the temperature and zero
C drift. 3 dimensional version
SUBROUTINE FCC
INCLUDE (PARAM)
DOUBLE PRECISION X(NP),Y(NP),Z(NP),PX(NP), PY(NP), PZ(NP), LATSIZ
DOUBLE PRECISION SUMPX, SUMPY, SUMPZ, SUMPSQ
INTEGER SEED, I,J,K,M, IJ, NCCC
COMMON /GAMMAG/ X,Y,Z,PX,PY,PZ
C These parameters are used if an initial colour current=J0 is required
DOUBLE PRECISION C, CALJ, SCALE
INCLUDE (STFUN)
LATSIZ = CUBSIZ * (NP/4) ** (-1/3.)
NCCC = (NP/4) ** (1./3) + 0.1
C make primitive cell
C 1st particle
X(1) = 0.25 * LATSIZ
Y(1) = 0.25 * LATSIZ
Z(1) = 0.25 * LATSIZ
C 2nd particle
X(2) = 0.75 * LATSIZ
Y(2) = 0.75 * LATSIZ
Z(2) = 0.25 * LATSIZ
C 3rd particle
X(3) = 0.25 * LATSIZ
Y(3) = 0.75 * LATSIZ
Z(3) = 0.75 * LATSIZ
C 4th particle
X(4) = 0.75 * LATSIZ
Y(4) = 0.25 * LATSIZ
Z(4) = 0.75 * LATSIZ
C Replicate primitive cell
C M = 6 * the number of the cell generated. I,J,K run through the cube i
C different dimensions, and IJ runs through the four particles in each c
C The first run through the loop is redundant
M = 0
DO 1 I=1,NCCC
DO 2 J=1,NCCC
DO 3 K=1,NCCC
DO 4 IJ = 1,4
X(IJ+M) = X(IJ) + (K-1) * LATSIZ
Y(IJ+M) = Y(IJ) + (J-1) * LATSIZ
Z(IJ+M) = Z(IJ) + (I-1) * LATSIZ
4 CONTINUE
M = M + 4
3 CONTINUE
2 CONTINUE
1 CONTINUE
C Generate random momenta
SUMPX = 0
SUMPY = 0
SUMPZ = 0
SUMPSQ = 0
SEED = 10101
DO 10 I=1,SEED0
CALL RAND(SEED,PX(1))
10 CONTINUE
DO 5 I=1,NP
CALL RAND( SEED, PX(I) )
CALL RAND( SEED, PY(I) )
CALL RAND( SEED, PZ(I) )
SUMPX = SUMPX + PX(I)
SUMPY = SUMPY + PY(I)
SUMPZ = SUMPZ + PZ(I)
SUMPSQ = SUMPSQ + PX(I)**2 + PY(I)**2 + PZ(I)**2
5 CONTINUE
C calculate centre of mass kinetic energy
SUMPSQ = SUMPSQ - 1./NP * (SUMPX**2 + SUMPY**2 + SUMPZ**2)
C Renormalize so that drift is zero and the kinetic energy is equal to t
C temperature
DO 6 I=1,NP
PX(I) = (PX(I) - SUMPX / NP) *
& SQRT( 3 * NP * TEMP / SUMPSQ)
PY(I) = (PY(I) - SUMPY / NP) *
& SQRT( 3 * NP * TEMP / SUMPSQ)
PZ(I) = (PZ(I) - SUMPZ / NP) *
& SQRT( 3 * NP * TEMP / SUMPSQ)
6 CONTINUE
C If an initial current is required, then adjust the momenta, and rescal
C adjust the kinetic energy
IF (CURRNT) THEN
C Compute initial colour current
C = CALJ(PX)
SCALE = SQRT( (J0**2-3*TEMP)/(C**2-3*TEMP))
C set up initial colour current
DO 7 I=1,NP
PX(I) = (PX(I) - E(I) * C) * SCALE + E(I) * J0
PY(I) = PY(I) * SCALE
PZ(I) = PZ(I) * SCALE
7 CONTINUE
ENDIF
RETURN
END
SUBROUTINE RAND(ISEED,R)
DOUBLE PRECISION DSEED,D2P31M, D2P31, R
DATA D2P31M/2147483647.D0/
DATA D2P31/2147483648.D0/
DSEED=ISEED
DSEED=DMOD(16807.D0*DSEED,D2P31M)
R=DSEED/D2P31
ISEED=INT(DSEED)
RETURN
END
DOUBLE PRECISION FUNCTION CALJ(PX)
INCLUDE (PARAM)
DOUBLE PRECISION PX(NP), J
INTEGER I
INCLUDE (STFUN)
J = 0
DO 1 I= 1,NP
J = J + PX(I) * E(I)
1 CONTINUE
CALJ = J / NP
RETURN
END
C return the kinetic energy as a function of phase space
DOUBLE PRECISION FUNCTION K(P)
INCLUDE (PARAM)
DOUBLE PRECISION SUMPSQ, P( D*NP )
INTEGER I
SUMPSQ = 0
DO 1 I=1,D*NP
SUMPSQ = SUMPSQ + P(I)**2
1 CONTINUE
K = SUMPSQ / 2
RETURN
END