| avm99963 | 802583e | 2022-05-12 10:56:56 +0200 | [diff] [blame] | 1 | program main |
| 2 | use, intrinsic :: iso_c_binding, only: sp=>c_float, dp=>c_double |
| 3 | implicit none |
| 4 | integer*4, PARAMETER :: L = 32 |
| 5 | integer*2 :: S(1:L, 1:L) |
| 6 | real*8 :: MAGNE, ENERG, TEMP, E, DIFE, DELTA, M, SUMA, W(-8:8), genrand_real2 |
| 7 | integer*4 :: PBC(0:L+1), I, J, IMC, MCTOT, MCINI, MCD, IPAS, N, SEED, NSEED, SEED0 |
| 8 | character :: TEMPSTRING*100, str*20 |
| 9 | character(:), allocatable :: NOM |
| 10 | |
| 11 | integer :: SUMI |
| 12 | real(dp) :: SUME, SUME2, SUMM, SUMAM, SUMM2, VARE, VARM |
| 13 | |
| 14 | if (command_argument_count() < 1) then |
| 15 | print *, "Error: please supply a command-line argument with the" |
| 16 | print *, "temperature, via './mc1 {TEMPERATURE}', where {TEMPERATURE}" |
| 17 | print *, "is the temperature." |
| 18 | stop |
| 19 | endif |
| 20 | |
| 21 | ! Inicialitzem algunes variables sobre el problema |
| 22 | call get_command_argument(1, TEMPSTRING) ! Ex: 2.0d0 |
| 23 | TEMPSTRING = trim(adjustl(TEMPSTRING)) |
| 24 | read (TEMPSTRING, *) TEMP |
| 25 | NSEED = 10 |
| 26 | SEED0 = 117654 |
| 27 | MCTOT = 10000 |
| 28 | MCINI = 1000 |
| 29 | MCD = 10 |
| 30 | N = L*L |
| 31 | NOM = "out_data/" // "SIM-L-" // trim(str(L)) // "-TEMP-" // & |
| 32 | trim(str(int(1000*TEMP))) // "-MCTOT-" // trim(str(MCTOT)) |
| 33 | |
| 34 | ! Obrim el fitxer on escriurem els resultats a cada iteració |
| 35 | open(unit = 12, file = NOM // ".out") |
| 36 | |
| 37 | ! Inicialitzem variable que ens fa latent la periodicitat |
| 38 | PBC(0) = L |
| 39 | PBC(L + 1) = 1 |
| 40 | do I = 1, L |
| 41 | PBC(I) = I |
| 42 | enddo |
| 43 | |
| 44 | ! Cache dels valors de l'exponencial |
| 45 | do I = -8, 8 |
| 46 | W(I) = exp(-float(I)/TEMP) |
| 47 | enddo |
| 48 | |
| 49 | do SEED = SEED0, SEED0 + NSEED - 1 |
| 50 | print *, "Seed: ", SEED |
| 51 | |
| 52 | ! Inicialitzem la matriu d'spins aleatòriament |
| 53 | call WRITECONFIGSEED(S, L, SEED) |
| 54 | |
| 55 | E = ENERG(S, L, PBC) |
| 56 | M = MAGNE(S, L) |
| 57 | print *, "Energia inicial:", E, "Magnet. inicial:", M |
| 58 | write (12, *) "0", E, M |
| 59 | |
| 60 | ! Iterem amb el mètode de Montecarlo |
| 61 | do IMC = 1, MCTOT |
| 62 | do IPAS = 1, N |
| 63 | ! LOS LOOPS HACEN COSAS |
| 64 | I = INT(genrand_real2()*L) + 1 |
| 65 | J = INT(genrand_real2()*L) + 1 |
| 66 | SUMA = S(PBC(I - 1), J) + S(PBC(I + 1), J) + S(I, PBC(J - 1)) + S(I, PBC(J + 1)) |
| 67 | DIFE = 2*SUMA*S(I, J) |
| 68 | |
| 69 | if (DIFE > 0) then |
| 70 | DELTA = genrand_real2() |
| 71 | if (DELTA >= W(int(DIFE))) then |
| 72 | ! NO ho acceptem |
| 73 | cycle |
| 74 | endif |
| 75 | endif |
| 76 | |
| 77 | ! Sí ho acceptem: |
| 78 | S(I, J) = -S(I, J) |
| 79 | E = E + DIFE |
| 80 | M = M + 2*S(I, J) |
| 81 | enddo |
| 82 | |
| 83 | if (mod(IMC, 1000) == 0) then |
| 84 | print *, "Iter:", IMC, "Energia:", E, "Magn:", INT(M) |
| 85 | endif |
| 86 | write (12, *) IMC, E, M |
| 87 | |
| 88 | if (IMC > MCINI .and. mod(IMC, MCD) == 0) then |
| 89 | SUMI = SUMI + 1 |
| 90 | |
| 91 | SUME = SUME + E |
| 92 | SUME2 = SUME2 + E*E |
| 93 | |
| 94 | SUMM = SUMM + M |
| 95 | SUMAM = SUMAM + abs(M) |
| 96 | SUMM2 = SUMM2 + M*M |
| 97 | endif |
| 98 | enddo |
| 99 | enddo |
| 100 | |
| 101 | ! Normalitzem els promitjos |
| 102 | SUME = SUME/SUMI |
| 103 | SUME2 = SUME2/SUMI |
| 104 | SUMM = SUMM/SUMI |
| 105 | SUMAM = SUMAM/SUMI |
| 106 | SUMM2 = SUMM2/SUMI |
| 107 | VARE = SUME2 - SUME*SUME |
| 108 | VARM = SUMM2 - SUMM*SUMM |
| 109 | |
| 110 | ! Guardem a un fiter els promitjos |
| 111 | open(unit = 13, file = NOM // ".res") |
| 112 | write(13, *) L, TEMP, SUMI, SUME, SUME2, VARE, SUMM, SUMAM, SUMM2, VARM |
| 113 | |
| 114 | ! Tanquem els fitxers de sortida |
| 115 | close(12) |
| 116 | close(13) |
| 117 | endprogram main |