| program main |
| use, intrinsic :: iso_c_binding, only: sp=>c_float, dp=>c_double |
| implicit none |
| integer*4, PARAMETER :: L = 32 |
| integer*2 :: S(1:L, 1:L) |
| real*8 :: MAGNE, ENERG, TEMP, E, DIFE, DELTA, M, SUMA, W(-8:8), genrand_real2 |
| integer*4 :: PBC(0:L+1), I, J, IMC, MCTOT, MCINI, MCD, IPAS, N, SEED, NSEED, SEED0 |
| character :: TEMPSTRING*100, str*20 |
| character(:), allocatable :: NOM |
| |
| integer :: SUMI |
| real(dp) :: SUME, SUME2, SUMM, SUMAM, SUMM2, VARE, VARM |
| |
| if (command_argument_count() < 1) then |
| print *, "Error: please supply a command-line argument with the" |
| print *, "temperature, via './mc1 {TEMPERATURE}', where {TEMPERATURE}" |
| print *, "is the temperature." |
| stop |
| endif |
| |
| ! Inicialitzem algunes variables sobre el problema |
| call get_command_argument(1, TEMPSTRING) ! Ex: 2.0d0 |
| TEMPSTRING = trim(adjustl(TEMPSTRING)) |
| read (TEMPSTRING, *) TEMP |
| NSEED = 10 |
| SEED0 = 117654 |
| MCTOT = 10000 |
| MCINI = 1000 |
| MCD = 10 |
| N = L*L |
| NOM = "out_data/" // "SIM-L-" // trim(str(L)) // "-TEMP-" // & |
| trim(str(int(1000*TEMP))) // "-MCTOT-" // trim(str(MCTOT)) |
| |
| ! Obrim el fitxer on escriurem els resultats a cada iteració |
| open(unit = 12, file = NOM // ".out") |
| |
| ! Inicialitzem variable que ens fa latent la periodicitat |
| PBC(0) = L |
| PBC(L + 1) = 1 |
| do I = 1, L |
| PBC(I) = I |
| enddo |
| |
| ! Cache dels valors de l'exponencial |
| do I = -8, 8 |
| W(I) = exp(-float(I)/TEMP) |
| enddo |
| |
| 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, PBC) |
| M = MAGNE(S, L) |
| print *, "Energia inicial:", E, "Magnet. inicial:", M |
| 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 = INT(genrand_real2()*L) + 1 |
| J = INT(genrand_real2()*L) + 1 |
| SUMA = S(PBC(I - 1), J) + S(PBC(I + 1), J) + S(I, PBC(J - 1)) + S(I, PBC(J + 1)) |
| DIFE = 2*SUMA*S(I, J) |
| |
| if (DIFE > 0) then |
| DELTA = genrand_real2() |
| 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, 1000) == 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*SUME |
| VARM = SUMM2 - SUMM*SUMM |
| |
| ! Guardem a un fiter els promitjos |
| open(unit = 13, file = NOM // ".res") |
| write(13, *) L, TEMP, SUMI, SUME, SUME2, VARE, SUMM, SUMAM, SUMM2, VARM |
| |
| ! Tanquem els fitxers de sortida |
| close(12) |
| close(13) |
| endprogram main |