ПРИЛОЖЕНИЕ Ж. РАСПЕЧАТКИ ПРОГРАММ НА ЯЗЫКЕ ФОРТРАН: ЧАСТЬ 2
Ниже приводятся распечатки перевода восьми программ, приведенных в ч. 2 текста, с языка True BASIC на Фортран-77. Программы пропускались иа машине VAX 11/750. Поскольку язык Фортран разрабатывался в те годы, когда машинная графика еще не получила широкого распространения, то стандартный Фортран не содержит никаких графических инструкций и подпрограмм. В качестве примера использования графического пакета в Фортране в наших программах проиллюстрировано применение PLOT 10 (фирма Tektronix)—распространенного пакета, имеющегося на машине VAX и других компьютерах. На большинстве компьютеров PLOT 10 может также вызываться из Паскаля. Краткое описание графических подпрограмм нз пакета PLOT 10, задействованных в наших программах, приведено в табл. П.
ГЛАВА 10
PROGRAM integ
* вычисляется интеграл от f(x) в пределах от х = а до х = b CALL initial(a,b, h, n)
CALL rectang!e(a, b, h, n, area)
CALL output(area)
STOP
END
SUBROUTINE initia!( a, b,h,n)
• a - нижний предел интегрирования, b - верхний предел а = 0.0
Ь = 0.5*( 3.14159)
WRITE(6,*) 'число интервалов = '
READ( 5, *) п
h = (Ь - а)/п
RETURN
END
SUBROUTINE rectangle( a, b, h. n, area) x = a sum = 0
DO 100 i = 0,n - 1
sum = sum + f(x)
x = x + h 100 CONTINUE
area - sum*h
RETURN
END
SUBROUTINE output(area) WRITE(6,13) area 13 FORMAT(2x, 'площадь = \F12.7) RETURN END
FUNCTION f(x) f = cos(x) RETURN END
ГЛАВА 11
PROGRAM rwalk
• моделирование Моите-Карло одномерного случайного блуждания DIMENSION prob(-64:64)
CALL start(p,N,ntrial, iseed) DO 100 itrial = l.ntria! CALL walk( ix, p, N, iseed)
* сбор данных после N шагов CALL data(ix,xcum, x2cum, prob)
100 CONTINUE
CALL aver( N, ntrial, xcum, x2cum, prob)
STOP
END
SUBROUTINE start( р, N, ntrial, iseed)
WRITE(6,*) 'начальное случайное число (положительное целое) = ' READ( 5,*) iseed
число испытаний ntrial = 100
вероятность шага вправо р = 0.5
WRITE(6,*) 'максимальное число шагов = '
READ(5,*) N
RETURN
END
SUBROUTINE wa!k(ix,p,N, iseed)
N-шаговое блуждание
начальное положение для каждого испытания ix = 0
DO 100 istep = I.N
IF (ran( iseed). le.p) then
ix = ix + 1 ELSE
ix = ix - 1 ENDIF 100 CONTINUE RETURN END
SUBROUTINE data( ix, xcum, x2cum, prob)
DIMENSION prob(-64:64)
xcum = xcum + ix
x2cum = x2cum + ix*ix
prob(ix) = prob(ix) + 1
RETURN
END
SUBROUTINE aver( N, ntrial, xcum, x2cum, prob)* средние значения для N-шагового блуждания
DIMENSION prob(-64:64)
znorm = 1.0/ntrial
xbar = xcum*znorm
x2bar = x2cum*znorm
DO 10 i = -N,N
prob(i) = prob( i)*znorm WRITE(6,13) i,prob(i) 13 FORMAT(2x,i6,f10.5) 10 CONTINUE
var = x2bar - xbar*xbar
sigma = sqrt(var)
WRITE(6,*) 'среднее смещение = ', xbar
WRITE(6,*) 'sigma = ', sigma
RETURN
END
ГЛАВА 12
PROGRAM percl
DIMENSION lat( 100,100), r( 100,100), icl( 100,100), np( 5000)DIMENSION mass(5000)CALL start(L, p, ntrial, iseed)DO 100 itrial = 1, ntrial35 CALL assign(L,r)
CALL occupy( L, p, lat, r)
CALL cluster( L, lat, icl, np)
проверка вертикального протекания CALL span(L, icl, np, l.ivert)
проверка горизонтального протекания CALL span(L, icl, np, 2, ihori)
IF ((ivert + ihori). eq.0) then
нет протекания ни в одном направлении ispan = 0
ELSEIF ((ivert*ihori).eq.0) then
протекание только в одном направлении ispan = ivert + ihori
IF (ivert.ne.O) spvert = spvert + 1 IF( ihori. ne. 0) sphori = sphori + 1 ELSE
* соединение в обоих направлениях IF (ivert.eq. ihori) then
iipan = ivert spvert = spvert + 1 sphori = sphori + 1 ELSE
* два соединяющих кластера —> это испытание пропускаем GO ТО 35
ENDIF ENDIF
вычисление массы каждого кластера CALL iize( L, id, np, mass, labmax)
вычисление среднего размера иесоединяющих кластеров CALL mean( labmax, mass, ispan, clsize, sum)
dtot = dtot + clsize
вычисление вероятности принадлежности соединяющему кластеру CALL Pinf(mass, ispan, sum, zPinf)
zPtot = zPtot + zPinf
вычисление радиуса циркуляции наибольшего несоед. кластера CALL radius( L, id, np, mass, labmax, ispan, radgyr)
radtot = radtot + radgyr 100 CONTINUE
CALL output( ntrial, spvert, sphori,dtot, zPtot, radtot)
STOP
END
SUBROUTINE start( L, p. ntrial, iseed)
WRITE(6,«) 'введите линейный размер решетки'
READ( 5, *) L
WRITE(6,*) 'введите вероятность занятия ячейки' READ{5,*) р
WRITE(6,*) 'enter seed, positive integer' READ(5,*) iseed
WRITE(6,«) 'введите число испытаний'
READ( 5, *) ntrial
RETURN
END
SUBROUTINE assign(L,r)
присваивание каждой ячейке решетки случайного числа DIMENSION R( 100,100)
DO 10 i = 1.L DO 10 j = 1.L
r{i,j) = ran(iseed) 10 CONTINUE RETURN END
SUBROUTINE occupy(L,p,lat,r)
занятие ячеек решетки DIMENSION Ы( 100,100), r{ 100,100) DO 10 i = 1,L
DO 10 j = I.L
IF (r(i,j).le.p) then
lat(i.j) = -1 ELSE
lat(i,j) = 0 ENDIF 10 CONTINUE RETURN END
SUBROUTINE cluster( L, lat, id, np)
маркировка ячеек
метки хранятся в массиве id, массив правильных меток пр DIMENSION lat( 100,100), icl( 100,100), np( 5000)
DO 10 k = 1,1000 np(k) = 0 10 CONTINUE
DO 15 i = I.L DO 15 j = 1,L icKi.i) = 0 15 CONTINUE nclust = 0
IF (lat(1,1).lt.0) CALL newcl(nclust,id, np, 1,1) DO 20 i = 2,L
IF (lat(i,1).lt.0) then ileft = i - 1
IF (lat(ileft,1).lt.0) then icl(i,1) = icl(ileft,1) ELSE
CALL newcl( nclust, icl, np,i,1) ENDIF ENDIF
20 CONTINUE
DO 30 j = 2,L
IF (lat(1,j).lt.0) then jdown = j - 1 IF (lat(1, jdown). It. 0) then icl(l.j) = icl(1, jdown) ELSE
CALL newcl( nclust, icl, np,1,j) ENDIF ENDIF
DO 40 i = 2,L
IF (lat(i.j).lt.O) then jdown = j - 1 ileft = i - 1
IF (icl(i,jdown) + icl(ileft,j).eq.O) then
CALL newcl( nclust, icl, np,i,j) ELSE
CALL neigh( lat, icl, np, i, j) ENDIF ENDIF
40 CONTINUE 30 CONTINUE
RETURN
END
SUBROUTINE newcl( nclust, icl, np, i, j)* инициируется новый кластер, указанный в подпрограмме cluster
DIMENSION icl( 100,100), пр{ 5000) nclust = nclust + 1 icl(i,j) = nclust np{ nclust) = 0 RETURN END
SUBROUTINE neigh( lat, icl, np, i, j)
присваивание кпастеру меток соседей DIMENSION lat(100,100), icl( 100,100), np( 5000) jdown = j - 1
Heft = i - 1
IF (lat(i,jdown)*lat(ileft,j).gt.0) then
CALL choice{icl, np, i,j, ileft,jdown)
RETURN ENDIF
IF (icl(i, jdown). gt.0) then
icl(i,j) = icl(i,jdown)
RETURN ENDIF
icl(i.j) = icl( ileft, j)
RETURN
END
SUBROUTINE choice(icl, np, i, j, ileft, jdown)
присваивание кпастеру меток при наличии выбора соседей DIMENSION icl( 100,100), пр( 5000)
IF (icl( ileft,j).eq.icl(i,jdown)) then
icl(i,j) = icl(ileft,j) ELSE
* соседи имеют разные метки nleft = icl( ileft, j)
ndown = icl(i, jdown) CALL proper(np, nleft) CALL proper( np, ndown) kmax = max0( nleft, ndown) kmin = min0( nleft, ndown) icl(i.j) = kmin
IF (kmax.ne.kmin)np(kmax) = kmin ENDIF RETURN END
SUBROUTINE properfnp, label)
определяется массив правильных меток пр DIMENSION пр{5000)
10 IF (пр{ label). eq.O) RETURN label = np( label) GO TO 10 END
SUBROUTINE span( L, id, np, idir, ispan)
выясняется существование соединяющего кластера проверкой
наличия ячеек с одинаковой правильной меткой на верхней и
нижней (idir = 1 - вертикальное протекание) ипи левой и
правой (idir * 1 - горизонтальное протекание) границах
решетки. На выходе: ispan = метке соединяющего кластера DIMENSION icl(100,100),np(5000)
DO 100 И = 1.L
IF (idir.eq. 1) then
iarg = i1
jarg = 1 ELSE
iarg = 1
jarg = И ENDIF
IF (icl( iarg, jarg). gt.O) then nl = icl( iarg, jarg) CALL proper(np, nl) DO 20 i2 = 1,L
IF (idir.eq.1) then iarg = i2 jarg = L ELSE
iarg = L jarg = i2 ENDIF
IF (icl(iarg,jarg).gt.O) then n2 = icl(iarg, jarg) CALL proper( np, n2) IF (n1.eq.n2) then ispan = nl RETURN ENDIF ENDIF
20 CONTINUE ENDIF
100 CONTINUE ispan = 0 RETURN END
SUBROUTINE size{ L, id, np, mass, labmax) DIMENSION icl( 100,100), np( 5000), mass (5000) DO 10 к = 1,1000 mass(k) = 0 10 CONTINUE labmax = 0 DO 100 i = 1,L DO 100 j = 1,L
IF (icl(i,j).gt.0) then label = icl(i,j) CALL proper( np, label) IF (label, gt. labmax) labmax = label mass( label) = mass( label) + 1 ENDIF 100 CONTINUE RETURN END
SUBROUTINE mean( labmax, mass, ispan, clsize, sum)
вычисление среднего размера несоединяющих кластеров DIMENSION mass(SOOO)
sum = О sum2 = О
DO 10 n = 1, labmax IF (n.ne. ispan) then
sum = sum + mass(n)
sum2 = sum2 + mass(n)**2 ENDIF
10 CONTINUE
clsize = sum2/sum
RETURN
END
SUBROUTINE Pinf( mass, ispan, sum, zPinf)
вычисление вероятности принадлежности соединяющему кластеру DIMENSION mass(SOOO)
IF (ispan. eq.O) then
* соединяющего кластера нет zPinf =0.0
RETURN ENDIF
zPinf = mass(ispan)/(sum + mass( ispan))
RETURN
END
SUBROUTINE radius( L, id, np, mass, labmax, ispan, radgyr)
расчет радиуса циркуляции наибольшего несоединяющего кластера
сначала находится наибольших несоединяющий кластер DIMENSION icl( 100,100), пр{ 5000), mass( 5000)
maxcl = 0 xcm = 0.0 уст =0.0 r2 = 0.0
DO 50 к = 1,labmax
IF ((mass(k).ne. 0).and.(k.ne. ispan)) then IF (mass(k).gt.maxcl) then maxcl = mass(k) maxk = к ENDIF ENDIF
50 CONTINUE
вычисление r**2 и центра масс DO 100 i = I.L
DO 100 j = I.L
IF (icl(i.j).gt.O) then label = icl(i.j) CALL proper(np, label) IF (label, eq. maxk) then xcm = xcm + i ycm = ycm + j r2 = r2 + i*i + j*j ENDIF ENDIF 100 CONTINUE
вычисляем радиус циркуляции (вычитая центр масс) r2cm = (хст*хст + ycm*ycm)/(maxcl*maxcl) radgyr = sqrt( r2/maxcl - г2ст)
RETURN END
SUBROUTINE output( ntrial, spvert, sphori, cltot, zPtot, radtot)
spvert = spvert/ntrial
sphori = sphori/ntrial
cltot = cltot/ntrial
zPtot = zPtot/ntrial
radtot = radtot/ntrial
WRITE(6,*) 'доля вертикального протекания = '.spvert WRITE(6,*) 'доля горизонтального протекания = '.sphori WRITE(6,*) 'средний размер кластера (иесоедиияющего) = '.cltot WRITE(6,*) 'вероятность принадлежности соед. кластеру = '.zPtot WRITE(6,*) 'радиус циркуляции наибольшего кластера = '.radtot RETURN END
PROGRAM invasion
расчет оккупирующего перколяционного кластера, рисование
кластера, печать относительной занятости и Р(г) DIMENSION г{ 0:100,0:20), регх( 500), регу( 500) CALL setup{ Lx, Ly, iseed)
CALL assign( Lx, Ly, r, perx, pery, iseed)
CALL invade{ perx, pery, r, Lx, Ly)
CALL aver( Lx, Ly, r)
STOP
END
SUBROUTINE sehip(Lx,Ly, iseed)
WRITE(6,*) 'Размер решетки в направлении у = ' READ(5,*) Ly
WRITE(6,*) 'начальное случайное число = ' READ(5,*) iseed Lx = 2*Ly
задание оконных координат CALL initt(1200)
CALL dwindo( 0.0, float( Lx+1), 0.0, float( Ly+2))
RETURN
END
SUBROUTINE assign( Lx, Ly, r, perx, pery, iseed)
присваивание случайных чисел каждой ячейке
занятие первого столбца
ячейки второго столбца образуют начальный лериметр DIMENSION г( 0:100,0:20), регх( 500), регу( 500)
DO 200 j = 1,Ly
занимается первый столбец
г{1.|) = 1 CALL box(l.j)
присваивание случайных чисел остатку j-й строки DO 100 i = 2,Lx
r(i,j) = ran( iseed) 100 CONTINUE
для ячеек периметра r(i,j) устанавливается больше 2 r(2,j) = 2* r(2,j)
sort perimeter sites
CALL sort(perx, pery, r,2, j, j) 200 CONTINUE RETURN END
SUBROUTINE sort(perx,pery,r, iO, jO, nper)
двоичная сортировка: список делится пополам; определяется, к
какой половине относится новое число; этот список тоже делится
пополам и определяется, в какой половине новое число. Процесс
продолжается до определения точного места нового числа DIMENSION г( 0:100,0:20), регх( 500), регу( 500)
если только одна ячейка периметра IF (nper.eq. 1) then
perx(nper) = iO pery(nper) = jO RETURN ENDIF
если новая ячейка меньше всех предыдущих ячеек периметра i = perx(nper-l)
j = pery(nper-1)
IF (r(i0,j0).lt.r(i,j)) then
perx(nper) = iO
pery(nper) = jO
RETURN ENDIF
начало сортировки Шелла, k2 - середина списка, k1,k3 - концы И = 1
k3 = nper - 1 к2 = (kl + кЗ)/2 100 CONTINUE
цикл до отыскания точного положения i = perx(k2)
j = pery(k2)
определение в какой половине списка находится новая ячейка
IF (r(iO,jO).gt.r{i,j)) then
кЗ = к2 ELSE
И = к2 ENDIF
новая середина к2 = (к1 + кЗ)/2
IF ((k1.eq.k2).or.(k2.eq.k3)) then
точное место найдено и равно кЗ
сдвигаем все ячейки выше кЗ на одну позицию вверх DO 10 кк = прег,кЗ+1,-1
perx(kk) = perx(kk-l) pery{kk) = регу(кк-1) 10 CONTINUE
регх(кЗ) = 10 регу(кЗ) = J0 RETURN ENDIF
GO TO 100
RETURN
END
SUBROUTINE sort2( perx, pery, r, iO, j0, nper)
стандартная сортировка методом вставок, элементы perx и pery с
наибольшим аргументом соответствуют наименьшему значению DIMENSION г( 0:100,0:20), регх( 500), регу( 500)
DO 200 к = 1,прег-1 i = регх(к) j = регу(к)
IF (r{i0,j0) .gt. r(i,j)) then
* вставляем новую ячейку DO 100 kk = nper,k+1,-1
perx(kk) = perx(kk-l) pery(kk) = pery(kk-l) 100 CONTINUE perx(k) = iO pery(k) = JO RETURN ENDIF
200 CONTINUE
новая ячейка меньше всех предыдущих ячеек периметра perx(k) = iO
pery(k) = JO
RETURN
END
SUBROUTINE invade( perx, pery, r, Lx, Ly)
массивы nnx and nny указывают положение ближайших соседей
относительно каждой ячейки
DIMENSION г( 0:100,0:20), регх( 500), регу( 500) DIMENSION nnx(4),nny(4) DATA ппх/1,-1,0,0/ DATA ппу/0,0,1,-1/ nper = Ly
DO 1000 idurnrny = 1,5000 i = perx(nper) j = pery(nper) nper = nper - 1
• занятые ячейки имеют значения между 1 и 2 r{i.j) = r<i,j) - 1
CALL box(i,j) DO 100 nn = 1,4
находим новые ячейки периметра nx = i + nnx(nn)
ny = j + nny(nn)
периодические краевые условия no у IF (ny .gt. Ly) then
ny = 1 ELSEIF (ny.lt.1) then
ny = Ly ENDIF
новая ячейка периметра
IF (r{nx,ny).lt.1) then
r{nx, ny) = r{nx,ny) + 2 nper = nper + 1
CALL sort( perx, pery, r,nx,ny, nper) ENDIF
100 CONTINUE
• выход по достижении кластером правой границы IF (i.gt.Lx) RETURN
1000 CONTINUE RETURN END
SUBROUTINE aver{Lx,Ly,r)
нахождение доли занятых ячеек н распределения вероятности
занятия для средней половины решетки DIMENSION г[ 0:100,0:20) DIMENSION р(0:20),п${0:20)
Lmin = Lx/3 Lmax = 2*Lmin
число ячеек в средней половине п = (Lmax - Lmin + 1)*Ly
DO 200 i = Lmin, Lmax DO 100 j = 1,Ly
iarg = 20*(amod(r(i, j),1.0)) ns(iarg) = ns(iarg) + 1 IF ((r(i,j) .gt. 1).and.(r(i,j) .It. 2) ) then occupied = occupied + 1
• плотность вероятности P(r) p(iarg) = p(iarg) + 1
ENDIF
100 CONTINUE 200 CONTINUE
для вывода результатов в файл по каналу 9 ждем нажатия клавиши fraction = occupied/n
CALL finitt(1,1)
WRITE(9,») 'доля оккупированных ячеек = ',fraction DO 300 i = 0,20 rr = i/20.0
IF (ns(i) .gt. 0) WRITE(9,*) rr, p(i)/ns(i) 300 CONTINUE
RETURN END
SUBROUTINE box(i,j)
вычерчивание прямоугольника, состоящего нз 11 полос хО = i - 0.5
у0 = j - 0.5 DO 10 k = 0,10
у = у0 + k*0.1
CALL movea{x0,Y)
• проводим прямую относитепьно х0,у CALL drawr(1.0,0.0)
10 CONTINUE CALL tsend RETURN END
ГЛАВА 14
PROGRAM entropy
вычисляется энтропия методом подсчета совпадений Ma DIMENSION rnleft( 10), mright( 10), rnicro( 0:2000)
CALL start( nl, nr, mleft, mright, micro, nexch, iseed)
обмен частнц
CALL exch( nl, nr, nexch, mleft, mright, micro, iseed)
расчет частоты совпадений н энтропии CALL output{ nexch, micro)
STOP END
SUBROUTINE start( nl, nr, mleft, mright, micro, nexch, iseed)
ввод параметров н выбор начальной конфигурации частиц DIMENSION mleft( 10), mright( 10), micro( 0:2000) WRITE(6,*) 'полное число частиц = '
READ( 5, *) N
WRITE(6,*) 'число частиц слева = ' READ( 5, *) nl
WRITE(6,*) 'начальное случайное число (целое положительное) = ' READ( 5, •) iseed
число частиц с правой стороны nr = N - nl
micro(O) = 0 DO 100 il = l.nl
список номеров частиц с левой стороны mleft(il) = il
начальное микросостояние micro(O) = micro(0)*2 + 2
100 CONTINUE
DO 200 ir = 1,nr
список номеров частиц с правой стороны mright(ir) = ir + nl
200 CONTINUE
WRITE(6,*) 'число обменов = '
READ( 5, •) nexch
RETURN
END
SUBROUTINE exch(nl, nr,nexch, mleft, mright, micro, iseed)
обмен частицы слева с номером i left
с частицей справа с номером fright DIMENSION mleft( 10), mright( 10), micro( 0:2000) DO 100 iexch = 1, nexch
индексы массивов выбираем случайным образом i left = int(ran(iseed)*nl + 1)
iright = int(ran(iseed)*nr + 1) jleft = mleft(ileft) jright= mright( iright)
новый номер частицы в левом массиве mlef»(ileft) = jright
новый номер частицы в правом массиве mright( iright) = jleH
определяем новое мнкросостояние micro(iexch) = micro(iexch - 1) + 2**jrighi micro(iexch) = micro(iexch) - 2**jleft
100 CONTINUE RETURN END
SUBROUTINE output(nexch, micro)
вычисляются частота совпадений и энтропия
полное число сравнений DIMENSION micro( 0:2000) ncomp = nexch*(nexch - 1)/2
сравниваем микросостояния DO 200 iexch = l.nexch - 1
DO 100 jexch = iexch + 1,nexch
IF (micro( iexch). eq.micro( jexch)) ncoin = ncoin + 1 100 CONTINUE 200 CONTINUE
частота совпадений
rate = float(ncoin)/float( ncomp) IF (rate.gt.0) S = alog(1.0/rate) WRITE(6,*) 'оценка энтропии = ', S RETURN END
ГЛАВА 15
PROGRAM conduct
алгоритм демона для одномерной модели Изинга
подвод и отвод тепла предусмотрен в концевых узлах DIMENSION s( 1000), edemon( 1000), edcum( 1000), zmgcum( 1000) CALL start( N, nmcs, s, zJ, iseed, nheat)
DO 100 imcs = l.nmcs DO 50 i = 1.N-2
istep = (imcs - 1)*(N-2) + i
IF ( mod( istep, nheat). eq.0) CALL heat( N, zJ, s, edemon) CALL change( N, edemon, s, zJ, accept, iseed) 50 CONTINUE
CALL data( N, s, edemon, zmgcum, edcum) 100 CONTINUE
CALL averg( N, nmcs, zmgcum, edcum, accept)
STOP
END
SUBROUTINE start( N, nmcs, s, zJ, iseed, nheat)
ввод параметров с клавиатуры, запись на устройство 9 DIMENSION s(1000)
WRITE(6,») 'число спинов = ' READ(5,*) N
WRITE(9,*) 'число спинов = ',N
WRITE(6,*) 'число шагов Монте-Карло на спин = ' READ(5,*) nmcs
WRITE(9,*) 'число шагов Монте-Карло на спин = '.nmcs WRITE(6,») 'константа обменного взаимодействия = ' READ( 5, •) zJ
WRITE(9,*) 'константа обменного взаимодействия = ',zJ WRITE(6,*) 'начальное случайное число = ' READ(5,*) iseed
WRITE(9,*) 'начальное случайное число = ',iseed WRITE(6,*) 'время между изменением концевых спинов = ' READ(5,») nheat
WRITE(9,*) 'время между изменением концевых спинов = ',nheat
DO 100 i = 1,N
• задание случайной начальной конфигурации IF (ran( iseed). gt. 0.5) then
s(i) = 1 ELSE
s(i) = -1 END IF
100 CONTINUE RETURN END
SUBROUTINE change( N, edemon, s, zJ, accept, iseed) DIMENSION s( 1000), edemon( 1000)
вычисляется случайный спин от 2-го до (N-l)-ro ispin = int(ran(iseed)*(N-2) + 2)
дииамнка опрокидывания спина, пробное изменение энергии = de de = 2*zJ*s(ispin)*(s(ispin-1) + s(ispin+1))
IF (de.le.edemon(ispin)) then
s( ispin) = -s( ispin)
accept = accept + 1
edemon( ispin) = edemon( ispin) - de ENDIF RETURN END
SUBROUTINE data( N, s, edemon, zmgcum, edcum)
накопление данных
DIMENSION s( 1000), edemon( 1000), edcum( 1000), zmgcumj 1000) DO 10 i = 1.N
edcum(i) = edcum(i) + edemon(i)
zmgcum(i) = zmgcum(i) + s(i) 10 CONTINUE
edemon(l) = 0 edemon(n) = 0 RETURN END
SUBROUTINE averg( N, nmcs, zmgcum, edcum, accept)
• нормирование накопленных данных н печать результатов DIMENSION zmgcum( 1000), edcum( 1000)
znorm = 1.0/nmcs
accept = accept*znorm/(N-2)
WRITE(6,*) 'коэффициент принятия = '.accept
WRITE(9,*) 'коэффициент принятия = '.accept
* величина edcum для спина 1 (ипи N) равна скорости поступления
• (отвода) тепла
WRITE(6,*) 'скорость подвода тепла = ',edcum(1)*znorm WRITE(6,*) 'скорость отвода тепла = ',edcum(N)*znorm DO 10 i = 2.N-1
edave = edcum( i)*znorm
IF (edave. ne. 0.0) temp = 4.0/alog(1.0 + 4.0/edave) zmag = zmgcum(i)*znorm WRITE(6,13) i, edave, temp, zmag WRITE(9,13) i, edave, temp, zmag 13 FORMAT(1x,i4,3(2x,f13.5)) 10 CONTINUE RETURN END
SUBROUTINE heat(N,zJ,s, edemon)
* производится попытка добавить (отвести) тепло у спина 1 (N) DIMENSION s( 1000), edemon( 1000)
IF ((s(1)*s(2).gt.0).and.(s(N)*s(N-1).lt.0)) then edemon(l) = edemon(1) + 2*zJ s(1) = -s(1)
edemon(N) = edemon(N) - 2*zJ
s(N) = -s(N) ENDIF RETURN END
ГЛАВА 16
PROGRAM Ising2
алгоритм Метрополиса для двумерной модели Изинга DIMENSION spin(32,32),w(-4:4)
CALL start( N, L, T, nmcs, spin, E, M, w, iseed) DO 100 imcs = 1,nmcs
CALL Melrop{ N, L, spin, E, M, w, ratio, iseed) CALL data( E, M, ecum, e2cum, mcum, m2cum) 100 CONTINUE
CALL output{N, nmcs,ecum, e2cum, mcum, m2cum, ratio)
STOP
END
SUBROUTINE start(N, L,T, nmcs, spin, E,M,w, iseed) DIMENSION spin(32,32),w(-4 :4) WRITE(6,*) 'линейный размер решетки = ' READ'S,*) L
число спинов N = L»L
WRITE(6,*) 'число шагов Монте-Карло на спин = ' READ(5,*) nmcs
WRITE(6,*) 'приведенная температура = ' READ( 5, *) Т
WRITE(6,*) 'начальное случайное число = ' READ{ 5, *) iseed
случайная начальная конфигурация DO 200 j = 1,L
DO 100 i = 1,L
IF (ran(iseed) .It. 0.5) then
spin(i,j) = 1 ELSE
spin(i,j) = -1 ENDIF
• суммарная намагниченность M = M + spin(i.j)
100 CONTINUE 200 CONTINUE
* вычисление начальной энергии Е DO 300 j = 1,L
IF (j .eq. L) then
iup = 1 ELSE
iup = j + 1 ENDIF
DO 400 i = 1,L
IF (i .eq. L) then
iright = 1 ELSE
iright = i + I ENDIF
sum = spin(i,iup) + spin( iright, j)
• полная энергия
E = E - spin( i, j)*sum 400 CONTINUE 300 CONTINUE
• вычисление вероятностей перехода
* индекс массива w равен сумме спинов четырех соседей е4 = ехр(-4.0/Т)
е8 = е4*е4
• аргумент функции ехр = изменению энергии w(4) = е8
w(-4) = е8 w(2) = е4 w(-2) = е4 RETURN END
SUBROUTINE Metrop( N, L, spin, E,M,w, ratio, iseed) DIMENSION spin(32,32),w(-4 :4) DO 100 ispin = 1,N
• случайный разыгрыш координат спина i = int(L*ran( iseed) + 1)
j = int(L*ran( iseed) + 1)
находим значения соседних спинов
периодические краевые условия CALL periodic(i, j,L, spin, isum)
IF ($pin(i,j)*i$um . le. 0 ) then
CALL accept( i, j, M, E, isum, spin, ratio)
ELSEIF (ran(iseed) .It. w(isum) ) then CALL accept{ i, j, M, E, isum, spin, ratio)
ENDIF
100 CONTINUE RETURN END
SUBROUTINE periodic(i,j,L,spin,isum)
находится сумма значений соседних спинов;
используются периодические краевые условия DIMENSION spin(32,32)
IF (i .eq. 1) then left = spin(L, j) ELSE
left = spin(i - 1,j) ENDIF
IF (i .eq. L) then
iright = spin(1,j) ELSE
iright = spin(i + 1,j) ENDIF
IF (j .eq. 1) then
idown = spin( i, L) ELSE
idown = spin(i,j - 1) ENDIF
IF (j .eq. L) then iup = spin(i,1) ELSE
iup = spin(i,j + 1) ENDIF
isum = left + iright + iup + idown
RETURN
END
SUBROUTINE ассер*{ i, j, М, Е, isum, spin, ratio)
DIMENSION spin(32,32)
spin(i.j) = -spin(i,j)
ratio = ratio + 1
M = M + 2*spin(i,j)
E = E - 2*spin{i,j)*isum
RETURN
END
SUBROUTINE data( E, M, ecum, e2cum, mcum, m2cum)
накапливание данных после каждого шага Монте-Карло на спин ecum = ecum + Е
e2cum = e2cum + Е*Е mcum = mcum + M m2cum = m2cum + M*M RETURN END
SUBROUTINE output( N, nmcs, ecum, e2cum, mcum, m2cum, ratio) znorm = 1/float(nmcs*N)
средние на спин ratio = ratio*znorm eave = ecum*znorm e2ave = e2cum*znorm save = mcum*znorm s2ave = m2cum*znorm
WRITE{6,*) 'коэффициент принятия = '.ratio WRITE{6,*) "средняя энергия на спин - '.eave WRITE{6,*) 'средний квадрат энергии на спин = ',e2ave WRITE{6,*) 'средняя намагниченность - ', save WRITE{6,*) 'средний квадрат намагниченности = ',s2ave RETURN END
ГЛАВА 17
PROGRAM wp
график движения квантовомеханического волнового пакета DIMENSION аг{ -100:100), ai( -100:100), prob( -100:100) DIMENSION prl(-100:100,-100:100),pim(-100:100,-100:100) CALL param( zkO, width, xO, xmin, h, dk, dt, tmax, m, N)
CALL packet( width, xO, m, dk, ar, ai)
CALL free{m,N,zkO,dk,h,ar,ai,prl,pim) CALL barrie( m, N, zkO, dk, h. ar, ai, prl, pim) ntime = tmax/dt
DO 10 itime = 0, ntime t = itime*dt
CALL move(t, m, N, zkO, dk, h, prl, pim, prob) CALL plot{ prob, xmi n, xmax, N, h) 10 CONTINUE
конец рисования, выставление курсора в позицию 1,1 CALL finitt(1,1)
STOP END
SUBROUTINE param( zkO, width, xO, xmin, h, dk, dt, tmax, m, N) WRITE{6,*) 'максимальное значение x = ' READ{ 5, *) xmax
WRITE{6,») 'минимальное значение x = '
READ{ 5, *) xmin
WRITE{6,*) 'размер шага = '
READ( 5, *) h
N = количество вычисляемых кординат N = (xmax - xmin)/h
WRITE(6,*) 'центральный волновой вектор = ' READ( 5, *) zkO
WRITE{6,*) 'число волновых векторов = ' READ{5,») m
WRITE(6,*) 'ширина волнового пакета в х-пространстве = ' READ{5,») width
WRITE(6,*) 'начальная координата пакета = ' READ{ 5, *) xO
dk = 2.0/{width*m) m = m/2 N = N/2
WRITE{6,*) 'шаг no времени = ' READ{ 5, *) dt
WRITE{6,*) 'полное время = ' READ{ 5, *) tmax
WRITE{6,*) 'максимальное значение квадрата psi = READ( 5, *) ymax CALL initt{1200)
CALL dwindo{xmin,xmax,-0.1,Ymax)
RETURN
END
SUBROUTINE packet( width, xO, m, dk, ar, ai)* вычисляются коэффициенты Фурье
DIMENSION ar{-100:100), ai{-100:100) sum2 = 0.0 DO 100 ik = -m,m dzk = ik*dk
a = exp{-dzk*dzk*width*width/4.0) ar{ik) = a*co${-dzk*x0) ai(ik) = a*sin(-dzk*x0)
$um2 = $um2 + ar{ik)*ar(ik) + ai(ik)*ai(ik) 100 CONTINUE
znorm = dk*$qrt(1.0/($um2)) DO 200 ik = -m,m
ar(ik) = ar(ik)*znorm ai(ik) = ai{ik)*znorm 200 CONTINUE RETURN END
SUBROUTINE free( m, N, zkO, dk, h, ar, ai, prl, pirn) DIMENSION prl{ -100:100, -100:100), pim( -100:100, -100:100) DIMENSION ar{-100:100),ai(-100:100) DO 100 ik = -m,m zk = zkO + ik*dk DO 50 ix = -N,N arg = zk*ix*dx fr = cos(arg) fi = sin(arg)
prl(ik, ix) = fr*ar(ik) - (i*ai(ik) pim(ik, ix) = fr*ai(ik) + fi*ar{ik) 50 CONTINUE 100 CONTINUE RETURN END
SUBROUTINE move{ », m, N, zkO, dk, h, prl, pirn, prob)DIMENSION prl{-100:100,-100:100), pim{-100:100,-100:100)DIMENSION prob( -100:100)DO 100 ix = -N, Npsir =0.0psii = 0.0DO 10 ik = -m,mzk = zkO + ik*dkarg = -0.5*zk*zk*tfr = cos(arg)fi = sin(arg)* вещественная и мнимая части psi
psir = piir + fr*prl{ ik, ix) - fi*pim(ik, ix) psii = psii + fi*prl( ik, ix) + fr*pim(ik, ix) 10 CONTINUE
prob(ix) = psir*psir + psii* psii 100 CONTINUE RETURN END
SUBROUTINE barrie( m, N, zkO, dk, h, ar, ai, prl, pirn) DIMENSION ar{-100:100),ai{-100:100)
DIMENSION prl{-100:100,-100:100), pim{-100:100,-100:100)
вычисление a(k)(exp(ik1x) + b exp(-iklx)) для x > 0
и a(k)*c*exp(ik2x) для x > 0 WRITE(6,*) 'высота барьера = ' READ{5,*) VO
DO 100 ik = -m, m zk1 = zkO + ik*dk E = 0.5*zk1*zk1 IF (E.gt.VO) then
zk2 = $qrt{2.0*(E - VO))
br = (zk1 - zk2)/{zk1 + zk2)
bi = 0.0
cr = 2.0*zk1/{zk1 + zk2) ci = 0.0 ELSE
zk2 = $qrt{2.0*{V0 - E)) denom = zk1*zk1 + zk2*zk2 br = (zk1*zk1 - zk2*zk2)/denom bi = 2.0*zk1*zk2/denom cr = 2.0*zk1*zk1/denom ci = -2.0*zk1*zk2/denom ENDIF
DO 40 ix = -Nx,0 arg = zk1*ix*h cs = cos(arg) si = sin(arg)
* далее использовано sin(-kx) =-sin(kx) и cos(-kx) = cos(kx) fr = cs + br*cs + bi*si
fi = si + bi*cs - br*si prl(ik, ix) = fr*ar(ik) - fi*ai(ik) pim(ik, ix) = fr*ai(ik) + fi*ar(ik) 40 CONTINUE
DO 60 ix = 1,N arg = zk2*ix*h IF (E.gt.VO) then
cs = cos(arg)
si = sin(arg) ELSE
cs = exp(-arg)
si = 0.0 ENDIF
(r - cr*cs - ci*si fi = ci *cs + cr*si pri(ik, ix) = fr*ar(ik) - fi*ai(ik) pim(ik, ix) = fr*ai(ik) + (i*ar(ik) 60 CONTINUE 100 CONTINUE RETURN END
SUBROUTINE plot(prob,xmin,N,h)
стирание экрана CALL newpag
нанесение делений DO 10 i = -N,N
x = i*h
CALL movea(x,0.0) CALL drawa(x.0.01*h) 10 CONTINUE
проведение горизонтальной прямой, после цикла х = xmax CALL movea(xmin, 0.0)
CALL drawa(x,0.0)
построение графика квадрата волновой функции CALL movea{-N*h,prob(-N))
DO 20 i = -N, N x = i*h
CALL drawa(x,prob(i)) 20 CONTINUE RETURN END