Afegir sessió 4 del lab de fenomens
Change-Id: I68a9c725acd6cde3fd43413fe673e5ec7794f7e6
diff --git a/quad10/fenomens/lab/p4/MC-2.f90 b/quad10/fenomens/lab/p4/MC-2.f90
new file mode 100644
index 0000000..3a2260d
--- /dev/null
+++ b/quad10/fenomens/lab/p4/MC-2.f90
@@ -0,0 +1,116 @@
+program main
+ use, intrinsic :: iso_c_binding, only: sp=>c_float, dp=>c_double
+ use iso_fortran_env
+ implicit none
+ integer*4 :: L
+ real*8 :: MAGNE, ENERG, TEMP, E, DIFE, DELTA, M, SUMA, V1, V2, V3, V4, W(-8:8)
+ integer*4 :: I, J, IMC, MCTOT, MCINI, MCD, IPAS, N, SEED, NSEED, SEED0, PBC
+ integer :: genrand_int31
+ character(100) :: NOM
+
+ integer :: SUMI
+ real(dp) :: SUME, SUME2, SUMM, SUMAM, SUMM2, VARE, VARM, CV, CHI
+
+ ! This should be replaced by L but it doesn't work dynamically, so we're
+ ! setting 64 as the upper limit which should be higher than L.
+ integer*2 :: S(1:64, 1:64)
+
+ ! Inicialitzem algunes variables sobre el problema
+ namelist /DADES/ L, NOM, TEMP, NSEED, SEED0, MCTOT, MCINI, MCD
+ open(10, file="mc2.dat")
+ read(10, DADES)
+ close(10)
+
+ N = L*L
+
+ ! Obrim el fitxer on escriurem els resultats a cada iteració
+ open(unit = 12, file = "data_out/" // trim(NOM) // ".out")
+
+ ! Cache dels valors de l'exponencial
+ do I = 0, 8
+ W(I) = (2**30 - 1 + 2**30)*exp(-float(I)/TEMP)
+ enddo
+
+ SUMI = 0
+ SUME = 0
+ SUME2 = 0
+
+ SUMM = 0
+ SUMAM = 0
+ SUMM2 = 0
+
+ do SEED = SEED0, SEED0 + NSEED - 1
+ print *, "Seed: ", SEED
+
+ ! Inicialitzem la matriu d'spins aleatòriament
+ call WRITECONFIGSEED(S, L, SEED)
+
+ E = ENERG(S, L)
+ M = MAGNE(S, L)
+ !write (12, *) "0", E, M
+
+ ! Iterem amb el mètode de Montecarlo
+ do IMC = 1, MCTOT
+ do IPAS = 1, N
+ ! LOS LOOPS HACEN COSAS
+ I = MOD(genrand_int31(), L) + 1
+ J = MOD(genrand_int31(), L) + 1
+ V1 = S(PBC(I - 1, L), J)
+ V2 = S(PBC(I + 1, L), J)
+ V3 = S(I, PBC(J - 1, L))
+ V4 = S(I, PBC(J + 1, L))
+ SUMA = V1 + V2 + V3 + V4
+ DIFE = 2*SUMA*S(I, J)
+
+ if (DIFE > 0) then
+ DELTA = genrand_int31()
+ if (DELTA >= W(int(DIFE))) then
+ ! NO ho acceptem
+ cycle
+ endif
+ endif
+
+ ! Sí ho acceptem:
+ S(I, J) = -S(I, J)
+ E = E + DIFE
+ M = M + 2*S(I, J)
+ enddo
+
+ if (mod(IMC, 10000) == 0) then
+ print *, "Iter:", IMC, "Energia:", E, "Magn:", INT(M)
+ endif
+ !write (12, *) IMC, E, M
+
+ if (IMC > MCINI .and. mod(IMC, MCD) == 0) then
+ SUMI = SUMI + 1
+
+ SUME = SUME + E
+ SUME2 = SUME2 + E*E
+
+ SUMM = SUMM + M
+ SUMAM = SUMAM + abs(M)
+ SUMM2 = SUMM2 + M*M
+ endif
+ enddo
+ enddo
+
+ ! Normalitzem els promitjos
+ SUME = SUME/SUMI
+ SUME2 = SUME2/SUMI
+ SUMM = SUMM/SUMI
+ SUMAM = SUMAM/SUMI
+ SUMM2 = SUMM2/SUMI
+ VARE = SUME2 - SUME**2
+ VARM = SUMM2 - SUMM**2
+ CV = (SUME2 - SUME**2)/real(N*TEMP**2)
+ CHI = (SUMM2 - SUMAM**2)/real(N*TEMP)
+
+ ! Guardem a un fiter els promitjos
+ open(unit = 13, file = "data_out/" // trim(NOM) // ".res")
+ !write(13, *) "L, T, <E>, <E**2>, Var(E), <M>, <M**2>, Var(M), C_V, CHI"
+ write(13, *) L, TEMP, SUME, SUME2, VARE, SUMM, SUMAM, SUMM2, VARM, CV, CHI
+
+ ! Tanquem els fitxers de sortida
+ close(12)
+ close(13)
+endprogram main