page@swan.ulowell.edu (Bob Page) (12/29/88)
Submitted-by: creutz@bnlux0.bnl.gov (michael creutz)
Posting-number: Volume 2, Issue 107
Archive-name: applications/assyising.1
# This is a shell archive.
# Remove everything above and including the cut line.
# Then run the rest of the file through sh.
#----cut here-----cut here-----cut here-----cut here----#
#!/bin/sh
# shar: Shell Archiver
# Run the following text with /bin/sh to create:
# assyising.bas
# This archive created: Wed Dec 28 16:16:55 1988
cat << \SHAR_EOF > assyising.bas
' program assyising.bas
'This Amiga Basic program calls an assembly language subroutine
'to simulate the two dimensional Ising model. Run it as
'an ordinary Amiga Basic program.
'The model represents a two dimensional magnet.
'The two colors represent domains of spin up and spin down.
'Increasing the parameter f will make the system more
'random, while decreasing it will cause the system to
'magnetize with one color dominating.
' The updating is done by 32 "demons" of two bits each.
' The demons run around the lattice flipping those spins
' for which they have enough energy. The algorithm is
' discussed in M. Creutz, Phys. Rev. Lett. 50, 1411 (1983)
' and G. Bhanot, M. Creutz and H. Neuberger, Nuclear
' Physics 235[FS11], 417 (1984).
' Michael Creutz
' creutz@bnlux0.bnl.gov
CLEAR ,30000&
WINDOW 1,,(20,100)-(200,180)
depth=1 ' if >1 then recent history on lower planes
SCREEN 1,320,200,depth,1
WINDOW 2,,,,1
PALETTE 0,0,0,0
PALETTE 1,.2,0,.8
GOSUB initcode
xmin=6 : ymin=6 ' upper left corner of display
sx1=160 ' x size, should be multiple of 32
sy1=175 ' y size
f!=.4 ' initial energy proportional to f**2
sx=sx1-1 : sy=sy1-1
LOCATE 1,23
PRINT "2-d Ising model"
LOCATE 3,23
PRINT "Lattice size"
LOCATE 4,24
PRINT sx1;"by";sy1
LOCATE 6,33
PRINT "sweeps"
LOCATE 9,24
PRINT "magnetization:"
LOCATE 12,24
PRINT "demon energy:"
LOCATE 15,24
PRINT "beta="
dimlat=2+depth*(sy+1)*INT((sx+16)/16)
latvol=(sy+1)*(sx+1)
DIM lat%(dimlat)
' draw initial pattern
DIM pat%(1)
pat%(0)=&H5555
pat%(1)=&HAAAA
PATTERN ,pat%
LINE (xmin,ymin)-STEP(sx*f!,sy*f!),,bf
GET (xmin,ymin)-(xmin+sx,ymin+sy),lat%
pat%(0)=&HFFFF
pat%(1)=&HFFFF
PATTERN ,pat%
LINE (xmin-4,ymin-4)-STEP(sx+8,sy+8),1,bf ' put frame around picture
' IMPORTANT: Predefine all variables so arrays don't move unexpectedly.
magnetization!=0!
energy!=0!
n=0
GOSUB findbeta
program&=0
alat&=0
program&=VARPTR(code%(0))
alat&=VARPTR(lat%(0))
loop:
code%(30)=0 ' turn off measuring
FOR n=1 TO 49
CALL program&(alat&)
PUT (xmin,ymin),lat%,PSET
NEXT n
code%(30)=-1 ' turn on measuring
CALL program&(alat&)
PUT (xmin,ymin),lat%,PSET
LOCATE 6,23
PRINT USING "#,###,###"; CVL(MKI$(code%(12))+MKI$(code%(13)))
magnetization!= (1!*CVL(MKI$(code%(16))+MKI$(code%(17))))/latvol
LOCATE 10,27
PRINT USING "#.####"; magnetization!
energy!= (1!*CVL(MKI$(code%(14))+MKI$(code%(15))))/latvol
LOCATE 13,27
PRINT USING "#.####"; energy!
GOSUB findbeta
LOCATE 15,29
PRINT USING "#.####"; beta!
GOTO loop
END
findbeta:
' Solve for beta given that the probability for a given
' energy is proportional TO EXP(-4*beta*energy)
' From exact solution the critical beta=log(1+sqr(2))/2.
bmin!=0! : bmax!=2! : beta!=.44
FOR nbeta=1 TO 25
ex!=EXP(-4*beta!)
ex2!=ex!*ex!
ex3!=ex!*ex2!
etrial!=(ex!+2*ex2!+3*ex3!)/(1!+ex!+ex2!+ex3!)
IF etrial!>energy! THEN
bmin!=beta!
beta!=.5*(beta!+bmax!)
ELSE
bmax!=beta!
beta!=.5*(bmin!+beta!)
END IF
NEXT nbeta
RETURN
initcode:
DIM code%( 232)
' ;ISING SIMULATION, LATTICE OBTAINED FROM BASIC GET
' ; CALL PROGRAM&(ALAT&)
code%( 0)=&H48E7 ' ISING MOVEM.L A0-A6/D0-D7,-(SP) ;SAVE REGISTERS
code%( 1)=&HFFFE '
code%( 2)=&H6000 ' BRA START
code%( 3)=&H38 '
code%( 4)=&H0 ' WORKSPACE DC.L 0 ; DEMON1 CODE%(4-5)
code%( 5)=&H0 '
code%( 6)=&H0 ' DC.L 0 ; DEMON2 CODE%(6-7)
code%( 7)=&H0 '
code%( 8)=&H5555 ' DC.L $55555555 ; PARITY MASK, CODE%(8-9)
code%( 9)=&H5555 '
code%( 10)=&H0 ' DC.L $44 ; BYTEJUMP, SHOULD BE MULTIPLE OF 4, CODE%(10-11)
code%( 11)=&H44 '
code%( 12)=&H0 ' DC.L 0 ; SWEEP COUNTER, CODE%(12-13)
code%( 13)=&H0 '
code%( 14)=&H0 ' DC.L 0 ; SUM OF DEMON ENERGIES OVER SWEEP, CODE%(14-15)
code%( 15)=&H0 '
code%( 16)=&H0 ' DC.L 0 ; TOTAL SPINS SET, CODE%(16-17)
code%( 17)=&H0 '
code%( 18)=&H0 ' DC.L 0 ; LATVOL 18
code%( 19)=&H0 '
code%( 20)=&H0 ' DC.L 0 ; ROWLENGTH 20
code%( 21)=&H0 '
code%( 22)=&H0 ' DC.L 0 ; LATSTART 22
code%( 23)=&H0 '
code%( 24)=&H0 ' DC.L 0 ; LATEND 24
code%( 25)=&H0 '
code%( 26)=&H0 ' DC.L 0 ; ROWSTART 26
code%( 27)=&H0 '
code%( 28)=&H0 ' DC.L 0 ; ROWEND 28
code%( 29)=&H0 '
code%( 30)=&HFFFF ' DC.W -1 ; MEASURE FLAG 30
' LAT EQU 64 ; (1+7+8)*4
' DEMON1 EQU 0
' DEMON2 EQU 4 ;SHIFTS FROM DEMON1
' MASK EQU 8
' JUMP EQU 12
' SWEEPS EQU 16
' ENERGY EQU 20
' MAG EQU 24
' LATVOL EQU 28
' ROWLENGTH EQU 32
' LATSTART EQU 36
' LATEND EQU 40
' ROWSTART EQU 44
' ROWEND EQU 48
' MEASURE EQU 52
code%( 31)=&H4DFA ' START LEA WORKSPACE(PC),A6
code%( 32)=&HFFC8 '
code%( 33)=&H42AE ' CLR.L ENERGY(A6)
code%( 34)=&H14 '
code%( 35)=&H42AE ' CLR.L MAG(A6)
code%( 36)=&H18 '
code%( 37)=&H206F ' MOVEA.L LAT(SP),A0 ;LAT ADDRESS IN A0
code%( 38)=&H40 '
code%( 39)=&H4280 ' CLR.L D0
code%( 40)=&H3010 ' MOVE.W (A0),D0 ; X SIZE IN D0
code%( 41)=&H240 ' ANDI.W #$FFE0,D0 ; ROUND DOWN TO LONGWORD
code%( 42)=&HFFE0 '
code%( 43)=&H30C0 ' MOVE.W D0,(A0)+ ; SAVE ROUNDED SIZE
code%( 44)=&H3218 ' MOVE.W (A0)+,D1 ; Y SIZE IN D1
code%( 45)=&H3418 ' MOVE.W (A0)+,D2 ; DEPTH IN D2
code%( 46)=&HE640 ' ASR.W #3,D0 ; BYTES PER ROW
code%( 47)=&HC2C0 ' MULU D0,D1 ; TOTAL NO OF BYTES PER PLANE IN D1
code%( 48)=&H2D41 ' MOVE.L D1,LATVOL(A6)
code%( 49)=&H1C '
code%( 50)=&H2D40 ' MOVE.L D0,ROWLENGTH(A6)
code%( 51)=&H20 '
code%( 52)=&H2D48 ' MOVE.L A0,LATSTART(A6)
code%( 53)=&H24 '
code%( 54)=&H2D48 ' MOVE.L A0,ROWSTART(A6)
code%( 55)=&H2C '
code%( 56)=&H2248 ' MOVEA.L A0,A1
code%( 57)=&HD3C1 ' ADDA.L D1,A1 ;END OF LATTICE
code%( 58)=&H2D49 ' MOVE.L A1,LATEND(A6)
code%( 59)=&H28 '
code%( 60)=&H2248 ' MOVEA.L A0,A1
code%( 61)=&HD3C0 ' ADDA.L D0,A1
code%( 62)=&H2D49 ' MOVE.L A1,ROWEND(A6)
code%( 63)=&H30 '
' ; AT THIS POINT:
' ; D0=BYTES PER ROW
' ; D1=BYTES PER PLANE
' ; D2=DEPTH
' ; A0=START OF LAT
' ; A6=ADDRESS OF DEMON1 AND WORKSPACE
' ; SHIFT TO LOWER PLANES IF DEPTH > 1
code%( 64)=&HC42 ' CMP.W #1,D2 ;IS DEPTH 1
code%( 65)=&H1 '
code%( 66)=&H6F16 ' BLE.S ONEPLANE
code%( 67)=&HC4C1 ' MULU D1,D2 ;TOTAL BYTES IN ALL PLANES
code%( 68)=&H2248 ' MOVEA.L A0,A1
code%( 69)=&H2448 ' MOVEA.L A0,A2
code%( 70)=&HD3C2 ' ADDA.L D2,A1 ;END OF LAST PLANE
code%( 71)=&H9481 ' SUB.L D1,D2 ;BYTES IN BOTTOM PLANES
code%( 72)=&HD5C2 ' ADDA.L D2,A2 ;END OF NEXT TO LAST PLANE
code%( 73)=&HE482 ' ASR.L #2,D2 ;LONG WORDS IN LOOP
code%( 74)=&H6002 ' BRA.S SHIFTLOOP
code%( 75)=&H2322 ' SHIFTLAT MOVE.L -(A2),-(A1) ; SHIFT WORDS
code%( 76)=&H51CA ' SHIFTLOOP DBF D2,SHIFTLAT
code%( 77)=&HFFFC '
' ; SET UP FOR MAIN LOOP
code%( 78)=&H2401 ' ONEPLANE MOVE.L D1,D2 ;USE D2 FOR LOOPCOUNTER
code%( 79)=&H7E1F ' MOVEQ.L #31,D7 ; USED FOR END CHECKS
code%( 80)=&H7C00 ' MOVEQ.L #0,D6 ; "
code%( 81)=&HE482 ' ASR.L #2,D2 ; D2=LOOP COUNTER, WORK WITH LONG WORDS
code%( 82)=&H2A6E ' MOVEA.L LATEND(A6),A5 ;KEEP LATEND IN A5
code%( 83)=&H28 '
code%( 84)=&H6000 ' BRA STARTLOOP
code%( 85)=&HD8 '
' ; SUBTRACT NEIGHBORS XOR OLD FROM NOT DEMONS INTO D3,D4,D5
' ; USE D0 FOR OLD, D1 AND D6 FOR SCRATCH
code%( 86)=&H2010 ' LOOP MOVE.L (A0),D0 ;OLD SPINS
code%( 87)=&H2616 ' MOVE.L (A6),D3 ;START WITH OLD DEMONS
code%( 88)=&H282E ' MOVE.L DEMON2(A6),D4
code%( 89)=&H4 '
code%( 90)=&H4683 ' NOT.L D3
code%( 91)=&H4684 ' NOT.L D4
code%( 92)=&H2A00 ' MOVE.L D0,D5 ; COPY OLD TO D5
code%( 93)=&HE38D ' LSL.L #1,D5
code%( 94)=&H4A11 ' TST.B (A1) ; BYTE ON RIGHT
code%( 95)=&H6C02 ' BGE.S NORIGHT
code%( 96)=&HDC5 ' BSET.B D6,D5 ; EAST NEIGHBOR IN D5
code%( 97)=&HB185 ' NORIGHT EOR.L D0,D5 ;XOR WITH OLD
code%( 98)=&HBB83 ' EOR.L D5,D3
code%( 99)=&HCA83 ' AND.L D3,D5 ;BORROW
code%( 100)=&HBB84 ' EOR.L D5,D4
code%( 101)=&HCA84 ' AND.L D4,D5 ;BORROW ; DONE WITH EAST #1
code%( 102)=&H2200 ' MOVE.L D0,D1
code%( 103)=&HE289 ' LSR.L #1,D1
code%( 104)=&HD12 ' BTST.B D6,(A2) ; TEST WEST BYTE FOR LOW ORDER BIT
code%( 105)=&H6702 ' BEQ.S NOLEFT
code%( 106)=&HFC1 ' BSET.L D7,D1 ; WEST IN D1
code%( 107)=&HB181 ' NOLEFT EOR.L D0,D1 ;XOR WITH OLD
code%( 108)=&HB383 ' EOR.L D1,D3
code%( 109)=&HC283 ' AND.L D3,D1 ;BORROW
code%( 110)=&HB384 ' EOR.L D1,D4
code%( 111)=&HC284 ' AND.L D4,D1 ;BORROW
code%( 112)=&HB385 ' EOR.L D1,D5 ;WEST DONE #2
code%( 113)=&H2213 ' MOVE.L (A3),D1 ;SOUTH IN D1
code%( 114)=&HB181 ' EOR.L D0,D1 ;XOR WITH OLD
code%( 115)=&HB383 ' EOR.L D1,D3
code%( 116)=&HC283 ' AND.L D3,D1
code%( 117)=&HB384 ' EOR.L D1,D4
code%( 118)=&HC284 ' AND.L D4,D1
code%( 119)=&HB385 ' EOR.L D1,D5 ;SOUTH DONE #3
code%( 120)=&H2214 ' MOVE.L (A4),D1 ;NORTH IN D1
code%( 121)=&HB181 ' EOR.L D0,D1 ;XOR WITH OLD
code%( 122)=&HB383 ' EOR.L D1,D3
code%( 123)=&HC283 ' AND.L D3,D1
code%( 124)=&HB384 ' EOR.L D1,D4
code%( 125)=&HC284 ' AND.L D4,D1
code%( 126)=&HB385 ' EOR.L D1,D5 ;NORTH DONE #4
code%( 127)=&H4683 ' NOT.L D3
code%( 128)=&HB985 ' EOR.L D4,D5 ; D5=REJECT
code%( 129)=&H8AAE ' OR.L MASK(A6),D5 ; OR WITH MASK
code%( 130)=&H8 '
' ; D3 AND D4 ARE NEW DEMONS
' ; D5=REJECT
' ; D0=OLD SPINS
code%( 131)=&H2205 ' MOVE.L D5,D1
code%( 132)=&H4681 ' NOT.L D1 ;ACCEPT
code%( 133)=&HB390 ' EOR.L D1,(A0) ; NEW SPINS
code%( 134)=&HCB96 ' AND.L D5,(A6)
code%( 135)=&HC681 ' AND.L D1,D3
code%( 136)=&H8796 ' OR.L D3,(A6) ;NEW DEMON1
code%( 137)=&HCBAE ' AND.L D5,DEMON2(A6)
code%( 138)=&H4 '
code%( 139)=&HC881 ' AND.L D1,D4
code%( 140)=&H89AE ' OR.L D4,DEMON2(A6) ;NEW DEMON2
code%( 141)=&H4 '
code%( 142)=&H46AE ' NOT.L MASK(A6) ;INVERT MASK
code%( 143)=&H8 '
code%( 144)=&H6D00 ' BLT LOOP ;REPEAT UPDATING FOR OTHER PARITY BITS
code%( 145)=&HFF8A '
' ; MAKE MEASUREMENTS:
code%( 146)=&H4A6E ' TST.W MEASURE(A6)
code%( 147)=&H34 '
code%( 148)=&H672C ' BEQ.S HOP ; IF MEASURE=0 THEN DON't measure things
code%( 149)=&HB380 ' EOR.L D1,D0 ; NEW SPINS
code%( 150)=&H2616 ' MOVE.L (A6),D3 ;DEMON1
code%( 151)=&H282E ' MOVE.L DEMON2(A6),D4 ;DEMON2
code%( 152)=&H4 '
code%( 153)=&H7A00 ' MOVEQ.L #0,D5
code%( 154)=&HE288 ' COUNTLOOP LSR.L #1,D0
code%( 155)=&HDB86 ' ADDX.L D6,D5 ; COUNT SET BITS IN LATTICE
code%( 156)=&H4A80 ' TST.L D0
code%( 157)=&H66F8 ' BNE.S COUNTLOOP
code%( 158)=&HDBAE ' ADD.L D5,MAG(A6) ; ACCUMULATE MAGNETIZATION
code%( 159)=&H18 '
code%( 160)=&HE28C ' D2COUNT LSR.L #1,D4
code%( 161)=&HD186 ' ADDX.L D6,D0 ; COUNT BITS IN DEMON2
code%( 162)=&H4A84 ' TST.L D4
code%( 163)=&H66F8 ' BNE.S D2COUNT
code%( 164)=&HE388 ' LSL.L #1,D0 ;DEMON2 BITS TIMES TWO
code%( 165)=&HE28B ' D1COUNT LSR.L #1,D3
code%( 166)=&HD186 ' ADDX.L D6,D0 ; COUNT BITS IN DEMON1
code%( 167)=&H4A83 ' TST.L D3
code%( 168)=&H66F8 ' BNE.S D1COUNT
code%( 169)=&HD1AE ' ADD.L D0,ENERGY(A6) ; ACCUMULATE ENERGY
code%( 170)=&H14 '
code%( 171)=&HD1EE ' HOP ADDA.L JUMP(A6),A0 ;JUMP TO NEW SITES
code%( 172)=&HC '
code%( 173)=&HB1EE ' ROWCHECK CMPA.L ROWEND(A6),A0 ;ADJUST ROW POINTERS
code%( 174)=&H30 '
code%( 175)=&H6D0E ' BLT.S ENDCHECK
code%( 176)=&H202E ' MOVE.L ROWLENGTH(A6),D0
code%( 177)=&H20 '
code%( 178)=&HD1AE ' ADD.L D0,ROWSTART(A6)
code%( 179)=&H2C '
code%( 180)=&HD1AE ' ADD.L D0,ROWEND(A6)
code%( 181)=&H30 '
code%( 182)=&H60EC ' BRA.S ROWCHECK
code%( 183)=&HB1CD ' ENDCHECK CMPA.L A5,A0 ;CHECK IF OFF END; A5=LATEND
code%( 184)=&H6D10 ' BLT.S STARTLOOP
code%( 185)=&H202E ' MOVE.L LATVOL(A6),D0
code%( 186)=&H1C '
code%( 187)=&H91C0 ' SUBA.L D0,A0
code%( 188)=&H91AE ' SUB.L D0,ROWSTART(A6)
code%( 189)=&H2C '
code%( 190)=&H91AE ' SUB.L D0,ROWEND(A6)
code%( 191)=&H30 '
code%( 192)=&H60EC ' BRA.S ENDCHECK
' ;FIND NEIGHBOR ADDRESSES WITH PERIODIC BOUNDARIES
code%( 193)=&H2248 ' STARTLOOP MOVEA.L A0,A1
code%( 194)=&H2448 ' MOVEA.L A0,A2
code%( 195)=&H2648 ' MOVEA.L A0,A3
code%( 196)=&H2848 ' MOVEA.L A0,A4
code%( 197)=&HD3FC ' ADDA.L #4,A1 ;EAST
code%( 198)=&H0 '
code%( 199)=&H4 '
code%( 200)=&H95FC ' SUBA.L #1,A2 ;WORD WEST
code%( 201)=&H0 '
code%( 202)=&H1 '
code%( 203)=&HD7EE ' ADDA.L ROWLENGTH(A6),A3 ;SOUTH
code%( 204)=&H20 '
code%( 205)=&H99EE ' SUBA.L ROWLENGTH(A6),A4 ;NORTH
code%( 206)=&H20 '
code%( 207)=&HB3EE ' CMPA.L ROWEND(A6),A1 ; IS EAST IN CORRECT ROW?
code%( 208)=&H30 '
code%( 209)=&H6D04 ' BLT.S EASTOK
code%( 210)=&H93EE ' SUBA.L ROWLENGTH(A6),A1
code%( 211)=&H20 '
code%( 212)=&HB5EE ' EASTOK CMPA.L ROWSTART(A6),A2 ; IS WEST?
code%( 213)=&H2C '
code%( 214)=&H6C04 ' BGE.S WESTOK
code%( 215)=&HD5EE ' ADDA.L ROWLENGTH(A6),A2
code%( 216)=&H20 '
code%( 217)=&HB7CD ' WESTOK CMPA.L A5,A3 ; IS SOUTH OK?
code%( 218)=&H6D04 ' BLT.S SOUTHOK
code%( 219)=&H97EE ' SUBA.L LATVOL(A6),A3
code%( 220)=&H1C '
code%( 221)=&HB9EE ' SOUTHOK CMPA.L LATSTART(A6),A4 ; IS NORTH?
code%( 222)=&H24 '
code%( 223)=&H6C04 ' BGE.S NORTHOK
code%( 224)=&HD9EE ' ADDA.L LATVOL(A6),A4
code%( 225)=&H1C '
code%( 226)=&H51CA ' NORTHOK DBF D2,LOOP
code%( 227)=&HFEE6 '
code%( 228)=&H52AE ' ADDQ.L #1,SWEEPS(A6)
code%( 229)=&H10 '
code%( 230)=&H4CDF ' MOVEM.L (SP)+,A0-A6/D0-D7 ;RESTORE REGISTERS
code%( 231)=&H7FFF '
code%( 232)=&H4E75 ' RTS
' END
RETURN
SHAR_EOF
# End of shell archive
exit 0
--
Bob Page, U of Lowell CS Dept. page@swan.ulowell.edu ulowell!page
Have five nice days.