N rovnakych bodovych atomov v dvojrozmernej stvorcovej krabici, interaguju
Lennard-Jonesovym potencialom. Odrazy od stien tiez cez LJ potencial, s tym,
ze ak sa to ukaze ako vhodne, nechame len odpudivu cast tohto potencialu.

Lennard-Jonesova potencialna energia medzi dvoma casticami vzdialenymi o R je
   U(R) = 4 * epsilonLJ * [ (sigmaLJ/R)^12 - (sigmaLJ/R)^6 ] .
Pre interakciu so stenami teda (ak chceme) vynechame 2. clen.

Vykresluju sa atomy, energie a aj histogram rozdelenia rychlosti. Ten sa
porovnava s teoretickym rozdelenim platnym veobecne pre sustavu velkeho
mnozstva (aj interagujucich, teda nielen pre idealny plyn) castic.

Pocita sa aj tlak a vztah medzi stavovymi velicinami sa potom sa porovnava
so vztahom (stavovou rovnicou) pre idealny plyn. Kedze plyn tu pocitany
nie je idealny, nemozme cakat rovnaky vysledok ako pre idealny. (A
samozrejme aj kvoli malemu poctu atomov.)

Ako uzitocnu literaturu odporucam napr. verejne dostupne prednasky
Prof. Tuckermana na adrese

http://www.nyu.edu/classes/tuckerman/stat.mech

Pozn.: Vypocty energii nie su sucastou Verletovho algoritmu. Robime ich
       len pre overenie presnosti integracie pohybovych rovnic. (Celkova
       energia sa musi zachovavat. Je to totiz sustava izolovana od okolia,
       comu sa v statistickej fyzike hovori mikrokanonicky subor.


V programe sa pouzivaju jednotky popisane v nasledovnom zozname. Su to
jednotky vybrate selfkonzistentne (take su aj SI), teda nie je potrebne
pouzivat vseliake prepocitavacie faktory.

tabulka zakladnych velicin (musia byt na sebe nezavisle) a ich jednotiek:

velicina ..... jednotka
------------------------
dlzka ........ sigmaLJ   (LJ parameter pre vzdialenost)
energia ...... epsilonLJ (LJ parameter pre energiu)
hmotnost ..... mass      (hmotnost jedneho atomu)
teplota ...... jednotku definujeme volbou Boltzmannovej konstanty kB.
               Zvolime kB = 1*jednotka energie / jednotka teploty.
               (To je ciselne 1).

odvodene veliciny a ich jednotky:

velicina ..... jednotka
--------------------------
sila ......... jednotka energie / jednotka dlzky = epsilonLJ/sigmaLJ
rychlost ..... sqrt(jednotka energie / jednotka hmotnosti)
atd.

Par poznamok ku grafickym proceduram: umoznuju ovladat grafiku priamo
z programu. Ak sa dobre pamatam, niekedy vsak aj tak treba do tych
grafov aj standartnym sposobom napisat nejake premenne, ktore sa budu
vykreslovat, ale moze to byt hocico (len aby tam nieco bolo). Kreslit sa
v skutocnosti bude to, co prikazeme prikazom Disp4.

SetScale(cislo grafu, xmin, xmax, ymin, ymax, spec)
        cislo grafu - od 1 do 4.
        Tato procedura vymaze okienko grafu a nakresli nove podla zadanych
        medzi. Celociselny parameter spec udava "specialitu" pouzitu
        pri oskalovani osi.
        SecScale sa vsak v tomto programe nevyskytuje, resp. je
        zakomentovana.

SetColor(cislo grafu, cislo farby)
        0 - cierna, 1 - biela, 2 - zlta, ..., 9 - tmavozelana
        Cislo farby mozme dat aj ZAPORNE. (-1, -2, ..., -9). V tom pripade
        sa pri opakovanom na kresleni znacky (na to iste miesto) znacka
        vymaze. To v programe vyuzivame, lebo atomy sa hybu, tak chceme,
        aby ich znacky z predoslych poloh zmizli.
        Cislo farby moze byt aj >= 10 (ale to sa uz opakuju.)

SetMark4(cislo grafu, cislo znacky) - ta 4 znamena, ze to bude platit pre
        kreslenie prikazom Disp4; ide o kreslenie specifikovane styrmi
        premennymi.
        centrovane znacky: 1 - kriz, 2 - obdlznik, 3 - elipsa
        necentrovane znacky: 4 - obdlznik, 5 - usecka, 6 - ciarkovana usecka
                             7 - vektor
Disp4(cislo grafu, x, y, rozmer_x, rozmer_y) - do prislusneho grafu
       vykresli znacku (urcenu predtym prikazom SetMark4) na suradnice x, y,
       pricom velkost tej znacky je dana cislami rozmer_x, rozmer_y.


poznamka k blokom v programe:
CTLR-b a nasledne posuvanie kurzoru a potvrdenie ENTER-om vytvori blok
CTRL-e expanduje blok
CTRL-s stlaci blok
CTRL-o odstrani blok
- - - - - - -  variables, constants, procedures and functions  - - - - - - -
sigmaLJ = 1.0                  ! LJ paramater (tu aj jednotka vzdialenosti)
epsilonLJ = 1.0                ! LJ parameter (tu aj jednotka energia)
mass = 1.0                     ! hmotnost kazdeho z atomov
kB = 1.0                       ! Boltzmannova konstatna

INT N = 10                     ! pocet atomov
L = 10.0                       ! velkost strany stvorcovej skatule
TempIni = 0.1  ! pociatocna teplota (Cislo s rozmerom teploty, ktore urcuje
               ! pociatocnu kin. energiu. O skutocnej pociatocnej teplote
               ! hovorit nemozme, lebo poc. stav nie je rovnovazny.)
dt = 0.0001                    ! casovy krok Verletovho algoritmu
tmax = 5000*dt
INT PeriKres = 20              ! Ako casto ma vykreslovat do grafov
                               ! (v kazdom PeriKres kroku).
LOGICAL KresliAtomy = TRUE     ! Aby sme usetrili vypoctovy cas, mozme
                               ! vypnut kreslenie atomov (FALSE).
INT PocetStlpcov = 40
INT histold[1 TO PocetStlpcov] ! "stary" histogram (kdesi je aj "novy" ako
                               ! lokalna premenna)

INT coord, atom, iatom, jatom, krok, stlpec
REAL KinEne, PotEneAtAt, PotEneAtWall, TotEne
REAL r[1 TO 2, 1 TO N], v[1 TO 2, 1 TO N], rold[1 TO 2, 1 TO N]
REAL rsave[1 TO 2, 1 TO N]
REAL forceWall[1 TO 2, 1 TO N], forceLJ[1 TO 2, 1 TO N]

PROCEDURE ForcesFromWalls(r[,], VAR forceWall[,]) ! pocita silu na kazdy
  INT coord, atom                                 ! atom od stien skatule
  REAL distance, AuxVar, forcepom

  FOR atom = 1 TO N DO FOR coord = 1 TO 2 DO
    forceWall[coord,atom] = 0.0

    distance = r[coord,atom]
    AuxVar = (sigmaLJ/distance)^6
    forcepom = 24*epsilonLJ/distance*AuxVar*(2*AuxVar-1.0)
!   Ked nechame len odpudivu interakciu, tak potom miesto predosleho riadku
    forcepom = 24*epsilonLJ/distance*AuxVar*2*AuxVar
    forceWall[coord,atom] = forceWall[coord,atom] + forcepom

    distance = L - r[coord,atom]
    AuxVar = (sigmaLJ/distance)^6
    forcepom = -24*epsilonLJ/distance*AuxVar*(2*AuxVar-1.0)
!   Ked nechame len odpudivu interakciu, tak potom miesto predosleho riadku
    forcepom = -24*epsilonLJ/distance*AuxVar*2*AuxVar
    forceWall[coord,atom] = forceWall[coord,atom] + forcepom
  END END
END

PROCEDURE ForcesLJ(r[,], VAR REAL forceLJ[,]) ! pocita silu na kazdy atom
  INT coord, iatom, jatom                     ! od vsetkych ostatnych atomov
  REAL rij2, AuxVar, forcepom

  FOR coord = 1 TO 2 DO FOR iatom = 1 TO N DO
    forceLJ[coord,iatom] = 0.0
  END END

  FOR iatom = 2 TO N DO FOR jatom = 1 TO iatom-1 DO
    rij2 = 0.0
    FOR coord = 1 TO 2 DO
      rij2 = rij2 + (r[coord,iatom]-r[coord,jatom])^2
    END
    AuxVar = (sigmaLJ*sigmaLJ/rij2)^3
    FOR coord = 1 TO 2 DO
      forcepom = (r[coord,iatom] - r[coord,jatom]) / rij2 * AuxVar
      forcepom = forcepom * (2*AuxVar-1.0)
      forceLJ[coord,iatom] = forceLJ[coord,iatom] + forcepom
      forceLJ[coord,jatom] = forceLJ[coord,jatom] - forcepom
    END
  END END

  FOR coord = 1 TO 2 DO FOR iatom = 1 TO N DO
    forceLJ[coord,iatom] = 24*epsilonLJ*forceLJ[coord,iatom]
    ! To je vysledok, teda sila (jej kartezska zlozka) na atom
    ! od vsetkych sostatnych atomov.
  END END
END

FUNCTION KinEnergy(v[,]) ! celkova kineticka energia sustavy
  INT coord, atom
  KinEnergy = 0.0
  FOR atom = 1 TO N DO FOR coord = 1 TO 2 DO
    KinEnergy = KinEnergy + 0.5*mass*v[coord,atom]*v[coord,atom]
  END END
END

FUNCTION PotEnergyAtAt(r[,]) ! potencialna energia sustavy (len cast
  INT coord, iatom, jatom    ! pochadzajuca od vzajomnej interakcie atomov);
  REAL rij2, AuxVar          ! "AtAt" je skratka od Atom-Atom
  PotEnergyAtAt = 0.0
  FOR iatom = 2 TO N DO FOR jatom = 1 TO iatom-1 DO
    rij2 = 0.0
    FOR coord = 1 TO 2 DO
      rij2 = rij2 + (r[coord,iatom]-r[coord,jatom])^2
    END
    AuxVar = (sigmaLJ*sigmaLJ/rij2)^3
    PotEnergyAtAt = PotEnergyAtAt + AuxVar*(AuxVar-1.0)
  END END
  PotEnergyAtAt = 4.0*epsilonLJ*PotEnergyAtAt
END

FUNCTION PotEnergyAtWall(r[,]) ! potencialna energia pochadzajuca
  INT atom, wall, coord        ! od interakcie atomov so stenami;
  REAL AuxVar                  ! skatula ma 4 "steny" teda, je tam cyklus
  PotEnergyAtWall = 0.0                                   ! wall = 1 TO 4
  FOR atom = 1 TO N DO FOR wall = 1 TO 4 DO
    coord = 1 + ((wall-1) DIV 2)
    IF wall MOD 2 = 1 THEN distance = r[coord,atom]
                      ELSE distance = L - r[coord,atom]
    END
    AuxVar = (sigmaLJ/distance)^6
!   PotEnergyAtWall = PotEnergyAtWall + AuxVar*(AuxVar-1.0)
!   Ked nechame len odpudivu interakciu, tak potom miesto predosleho riadku
    PotEnergyAtWall = PotEnergyAtWall + AuxVar*AuxVar
  END END
  PotEnergyAtWall = 4.0*epsilonLJ*PotEnergyAtWall
END

PROCEDURE histogram(v[,], INT histold[]) ! Tento podprogram kresli
  INT atom, stlpec                       ! histogram rozdelenia rychlosti.
  INT hist[1 TO PocetStlpcov]

  ! Najdi najvyssiu rychlost spomedzi vsetkych atomov, aby sme vedeli,
  ! aka dlha ma byt os x.
  vmax = 0.0
  FOR atom = 1 TO N DO
    vpom = v[1,atom]^2 + v[2,atom]^2
    IF vpom > vmax THEN vmax = vpom END
  END
  vmax = sqrt(vmax)

  SirkaStlpca = vmax/PocetStlpcov

  ! Vynuluj histogram.
  FOR stlpec = 1 TO PocetStlpcov DO hist[stlpec] = 0 END

  ! Zarad rychlosti jednotlivych atomov do histogramu.
  FOR atom = 1 TO N DO
    vpom = sqrt(v[1,atom]^2 + v[2,atom]^2)
    stlpec = 1 + trunc(vpom/SirkaStlpca)
    IF stlpec = PocetStlpcov + 1 THEN stlpec = stlpec-1 END
    hist[stlpec] = hist[stlpec] + 1
  END

  ! Najdi najvyssi stlpec, aby sme vedeli, aka velka ma byt os y.
  histmax = 0
  FOR stlpec = 1 TO PocetStlpcov DO
    IF hist[stlpec] > histmax THEN histmax = hist[stlpec] END
  END

  ! Prekreslenim zmaz stary histogram.
  FOR stlpec = 1 TO PocetStlpcov DO
    Disp4(4, (stlpec-1)*vmax/PocetStlpcov, 0, SirkaStlpca, histold[stlpec])
  END

  ! Vykresli histogram.
  SetScale(4, 0.0, 0.3*vmax, 0.0, 4.0*histmax, 1)
  SetMark4(4,4)
  SetColor(4,-1)
  FOR stlpec = 1 TO PocetStlpcov DO
    Disp4(4, (stlpec-1)*vmax/PocetStlpcov, 0, SirkaStlpca, hist[stlpec])
  END

  ! Pre porovnanie budeme kreslit aj histogram podla teoretickeho vztahu.

  ! Spocitaj teplotu potrebnu pre kreslenie teoretickeho histogramu.
  KinEne = 0.0
  FOR atom = 1 TO N DO
    KinEne = KinEne + 0.5*mass*(v[1,atom]^2 + v[2,atom]^2)
  END
  T = KinEne/(N*kB)

  ! Spocitaj normovaciu konstantu potrebnu pre kreslenie teoreticeho
  ! histogramu.
  suma = 0.0
  FOR stlpec = 1 TO PocetStlpcov DO
    vpom = (stlpec-0.5)/PocetStlpcov*vmax
    suma = suma + vpom * exp(-0.5*mass*vpom^2/(kB*T))
  END

  ! Vykresli teoreticky histogram.
  SetMark(4,3)
  FOR stlpec = 1 TO PocetStlpcov DO
    vpom = (stlpec-0.5)/PocetStlpcov*vmax
    hodnota = N/suma*vpom * exp(-0.5*mass*vpom^2/(kB*T))
    Disp(4,vpom,hodnota)
  END

  ! Uloz sucasny histogram do histold, aby sa pri nasledujucom vstupe
  ! do tejto procedury sucasny histogram mohol vymazat.
  FOR stlpec = 1 TO PocetStlpcov DO
    histold[stlpec] = hist[stlpec]
  END
END

- - - - - - - - - - - - - - - initial values   - - - - - - - - - - - - - - -

! Grafika pre graf 1 (zobrazenie stvorcovej krabice).
SetMark4(1,4)     ! kreslenou znackou bude pravouhly stvoruholnik
Disp4(1,0,0,1,1)  ! stvorcova krabica urcena suradnicami (0,0) a (1,1)

! Generuj pociatocne polohy tak, aby ziadne dva atomy neboli blizsie
! ako MinDistAccept ani neboli v stene ani za stenou
MinDistAccept = 0.5*sigmaLJ
FOR coord = 1 TO 2 DO
  r[coord,1] = MinDistAccept + (L-2*MinDistAccept)*rnd()
END

FOR iatom = 2 TO N DO
  REPEAT
    FOR coord = 1 TO 2 DO
      r[coord,iatom] = MinDistAccept + (L-2*MinDistAccept)*rnd()
    END
    rijsqmin = 10*L*L   ! na zaciatok nastav na velku hodnotu
    FOR jatom = 1 TO iatom-1 DO
      rijsq = (r[1,iatom]-r[1,jatom])^2 + (r[2,iatom]-r[2,jatom])^2
      IF rijsq < rijsqmin THEN rijsqmin = rijsq END
    END
  UNTIL rijsqmin > (2*MinDistAccept)^2
END
WRITELN 'POCIATOCNE POLOHY VYGENEROVANE; STLAC SPACE'; STOP

! Generuj pociatocne rychlosti normovane na teplotu TempIni.
! Ma teda platit Sum[ 0.5*mass*vi*vi, {i,1,N} ] = N * 2/2 * kB*TempIni
KinEne = 0.0
FOR coord = 1 TO 2 DO FOR atom = 1 TO N DO
  v[coord,atom] = -1.0 + 2*rnd()
  KinEne = KinEne + 0.5*mass*v[coord,atom]*v[coord,atom]
END END
FOR coord = 1 TO 2 DO FOR atom = 1 TO N DO
  v[coord,atom] = v[coord,atom]*sqrt(N*kB*TempIni/KinEne)
END END

! Spocitaj a vykresli histogram pociatocneho rozdelenia rychlosti.
! Na zaciatku potrebujeme nastavit histold.
FOR stlpec = 1 TO PocetStlpcov DO
  histold[stlpec] = 0
END
histogram(v[,], histold[])

! Uloz pociatocne polohy do rold (v prvom kroku to budeme potrebovat).
FOR coord = 1 TO 2 DO FOR atom = 1 TO N DO
  rold[coord,atom] = r[coord,atom]
END END

! Spocitaj a zakresli energie v case 0.

  ! Spocitaj energie.
KinEne = KinEnergy(v[,])
PotEneAtAt = PotEnergyAtAt(r[,])
PotEneAtWall = PotEnergyAtWall(r[,])
PotEne = PotEneAtAt + PotEneAtWall
TotEne = KinEne + PotEneAtAt + PotEneAtWall
TotEne0 = TotEne ! Potrebujeme kvoli kresleniu grafov.
WRITELN 'TotEne0 = ', TotEne0, '   STLAC SPACE'; STOP

  ! Nastav niektore parametre pre kreslenie energii v zavislosti od casu.
      ! graf 2: dve zavislosti: kineticka a potencialna energia
SetScale(2, 0.0, tmax, -0.2*TotEne0, 1.2*TotEne0, 1)
SetMark(2,1)
      ! graf 3: ku stavovej rovnici
SetScale(3, 0.0, tmax, 0.0, 1.5*TotEne0, 1)
SetMark(3,2)

  ! Zaznac jednotlive zlozky energie v case 0 do grafu 2.
SetColor(2,1)
Disp(2, 0, KinEne)
SetColor(2,2)
Disp(2, 0, PotEneAtAt + PotEneAtWall)
SetColor(2,4)
SetMark(2,2)
Disp(2, 0, TotEne)
SetMark(2,1)

! Spocitaj sily v case 0.
ForcesFromWalls(r[,],forceWall[,])
ForcesLJ(r[,],forceLJ[,])

! Sprav nulty krok (jediny, ktory sa nerobi Verletovym algoritmom).
! Dostavame polohy v case dt.
FOR coord = 1 TO 2 DO FOR atom = 1 TO N DO
  accel = (forceWall[coord,atom] + forceLJ[coord,atom])/mass
  r[coord,atom] = r[coord,atom] + v[coord,atom]*dt + 0.5*accel*dt*dt
END END

krok = 0
t = dt
- - - - - - - - - - - - - - - - -  model   - - - - - - - - - - - - - - - - -
krok = krok+1

! Vykresli atomy v case t.
IF KresliAtomy AND krok MOD PeriKres = 0 THEN
  SetMark4(1,3)      ! 3: kreselnou znackou (vo vseobecnosti) elipsa
  FOR atom = 1 TO N DO
    SetColor(1,-atom)
    Disp4(1,r[1,atom]/L,r[2,atom]/L,sigmaLJ/L,sigmaLJ/L)
  END
  Pause(0.02)
END

! Spocitaj sily v case t.
ForcesFromWalls(r[,],forceWall[,])    ! sily pre polohy v case t
ForcesLJ(r[,],forceLJ[,])

! Uloz r(t) do rsave.
!(Kedze v jadre Verletovho algoritmu si prepiseme r(t) updatovanim
! na r(t+dt), pred prepisanim si to musime odlozit.
! Do rold to vsak odlozit nemozme.)
FOR coord = 1 TO 2 DO FOR atom = 1 TO N DO
  rsave[coord,atom] = r[coord,atom]
END END

! Jadro Verletovho algoritmu:  r(t+dt) = 2*r(t) - r(t-dt) + a(t)*dt*dt
FOR coord = 1 TO 2 DO FOR atom = 1 TO N DO
  accel = (forceWall[coord,atom] + forceLJ[coord,atom])/mass
  r[coord,atom] = 2*r[coord,atom] - rold[coord,atom] + accel*dt*dt
END END

! Tieto veci rob len v tych krokoch, v ktorych vykreslujeme do grafov.
IF krok MOD PeriKres = 0 THEN

  ! Rychlosti (v case t) potrebujeme len pre vypocet kinetickej energie
  ! a pre kreslenie histogramu.
  FOR coord = 1 TO 2 DO FOR atom = 1 TO N DO
    v[coord,atom] = (r[coord,atom] - rold[coord,atom])/(2*dt)
  END END

  ! Spocitaj a vykresli histogram rozdelenia rychlosti v case t.
  histogram(v[,], histold[])

  ! Spocitaj energie v case t.
  KinEne = KinEnergy(v[,])
  PotEneAtAt = PotEnergyAtAt(rsave[,])
  PotEneAtWall = PotEnergyAtWall(rsave[,])
  PotEne = PotEneAtAt + PotEneAtWall
  TotEne = KinEne + PotEneAtAt + PotEneAtWall

  ! Zaznac jednotlive zlozky energie aj celkovu energiu v case t do grafu 2.
  SetColor(2,1)
  Disp(2, t, KinEne)
  SetColor(2,2)
  Disp(2, t, PotEneAtAt + PotEneAtWall)
  SetColor(2,4)
  SetMark(2,2)
  Disp(2, t, TotEne)
  SetMark(2,1)

  ! Spocitaj tlak v case t podla vztahu (vyplyva z teorie)
  !
  ! P = 1/(2*L*L) * Sum[ pi*pi/(2*mass) + ri * Fi , {i,1,N} ],
  !
  ! kde ri a pi su vektory polohy a hybnosti i-tej castice,
  ! Fi je vektor celkovej sily na i-tu casticu

  tlak = 0.0
  FOR atom = 1 TO N DO
    force_x = forceWall[1,atom] + forceLJ[1,atom]
    force_y = forceWall[2,atom] + forceLJ[2,atom]
    tlak = tlak + rsave[1,atom]*force_x + rsave[2,atom]*force_y
  END
  tlak = tlak/(2*L*L)
  tlak = tlak + KinEne/(L*L)

  ! Ideme porovnat, ci veliciny ako tlak, teplota a "objem" (tu plocha)
  ! vypocitane v nasej simulacii aspon priblizne splnaju stavovu rovnicu
  ! idealneho plynu: tlak*plocha = kineticka energia.
  ! (Pre trojrozmernu krabicu by bolo tlak*objem = 2/3 kinetickej energie).

  LS = tlak*L*L               ! lava strana stav. rovnice id. plynu
  PS = KinEne                 ! prava ----||----
  WRITELN 't, TotEne, LS, PS = ', t:4:-3, TotEne:5:-4, LS:5:-4, PS:5:-4

  ! graf 3: lava a prava strana stavovej rovnice idealneho plynu
  SetColor(3,1)
!  Disp(3,t,LS)
  SetColor(3,2)
!  Disp(3,t,PS)
END

t = t + dt

! Polohy r su uz pre cas t+dt. Uloz r(t) do rold.
FOR coord = 1 TO 2 DO FOR atom = 1 TO N DO
  rold[coord,atom] = rsave[coord,atom]
END END

IF KresliAtomy AND krok MOD PeriKres = 0 THEN
  ! Prekresli atomy, aby zo starych poloh zmizli.
  FOR atom = 1 TO N DO
    SetColor(1,-atom)
    Disp4(1,rold[1,atom]/L,rold[2,atom]/L,sigmaLJ/L,sigmaLJ/L)
  END
END

- - - - - - - - - - - - - - - - - - disp - - - - - - - - - - - - - - - - - -