Culo

/*
   Metropolis Monte Carlo simulation of a 2D Ising system

   Cameron F. Abrams

   Written for the course CHE 800-002, Molecular Simulation
   Spring 0304

   compile using "gcc -o ising ising.c -lm -lgsl"
   (assumes the GNU Scientific Library is installed)

   runs as "./ising -L <sidelength(20)> -T <temperature(1.0)> \
                    -nc <numcycles(1e6)> -fs <samplefreq(100)> \
            -s <seed(?)>"

   For example, to run a 20x20 system at T = 0.5 for 1e7 cycles
   sampling every 100 cycles, the command looks like

          ./ising -L 20 -T 0.5 -nc 1e7 -fs 100

   The default values are shown in parentheses above.

   You must have the GNU Scientific Library installed; see
   the coursenotes to learn how to do this.

   The Hamiltonian is

   H = -J sum_<ij> s_i * s_j,

   where "sum_<ij>" is the sum over all unique
   nearest neighbor pairs, s_i = +/- 1, and J
   is an arbitrary "coupling" parameter having
   units of energy.  We work in a unit system
   where energy is measured in units of J and
   temperature is nondimensionalized as (k_B*T)/J.

   Drexel University, Department of Chemical Engineering
   Philadelphia
   (c) 2004
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include <gsl/gsl_rng.h>

//Computa el cambio de energía debido al giro del spin i,j,
//que corresponde a la energía asociada a los primeros vecinos
//No computa diagonales
int E( int * F, int L, int i, int j)
{
    //(i - 1) < 0 ? L - 1 : i - 1 es una expresión que da L - 1 si (i - 1) < 0
    //es decir se fue de la grilla (condiciones periodicas). Si no da (i - 1).
    //La energía la
    return -2 * F[i*L+j] * (F[((i - 1) < 0 ? L*(L - 1) : L*(i - 1))+ j + F[((i+1) % L)*L+j] +
              F[i*L + ((j-1) < 0 ? L - 1 : j - 1)] + F[i*L+((j + 1) % L)]);
}

//Función que mide la magnetización y la energía por sitio
//No devuelve nada, se pasa los valores de magnetización m y energía e
void sampleo(int *F, int L, double *m, double *e)
{
    int i,j;
    *m=0.0;
    *e=0.0;
    for (i = 0; i < L; i++)
    {
        for (j = 0; j < L; j++)
        {
            *m += (double)F[i*L+j];
            //Como sumo sobre todos los índices, con calcular la energía con los vecinos "posteriores"
            //es suficiente
            *e -= (double) (F[i*L+j])*(F[i*L+(j+1)% L]+F[((i+1)%L)*L+j]);
        }
    }
    *m /= (L*L);
    *e /= (L*L);
}

/* Randomly assigns all spins */
void init(int *F, int L)
{
  int i,j;
  /* Visit each position (i,j) in the lattice */
  for (i=0; i<L; i++)
  {
    for (j=0; j<L; j++)
    {
        //Número entero al azar entre -1 y 1
        F[i*L+j]= (int) rand() % 2 - 1;
        if (F[i*L+j] == 0)
        {
            F[i*L+j] = 1;
        }
    }
  }
}

int main(int argc, char * argv[])
{
    //Parámetros de sistema
    int L = 20; //Tamaño de la red
    int N = L*L; //Cantidad de spines
    double beta = 1.0; //Temperatura adimensional  beta = J/(k*T)

    //Parametro de corrida
    int nIter = 2000000; //Pasos de MC, cada uno tiene N cambios de spin
    int nTerm = 1000000; //Pasos hasta termalizar
    int fSamp = 1000;  //Cantidad de pasos entre mediciones

    //Variables internas de cálculo
    int nSamp = 0; //Número de muestras, para obtener el promedio de mag y energia
    int de;   //Delta de energía
    double p; //Factor de aceptación
    double r; //Número aleatorio, para comprar con la aceptación
    int i,j,a,c;  //Contadores

    //Observables
    double m = 0.0, msum=0.0;    //Magentización
    double e=0.0, esum=0.0;    //Energía

  //Obtencion de parámetros
    for (i=1;i<argc;i++)
    {
        if (!strcmp(argv[i],"-L"))
        {
            L = atoi(argv[++i]);
        }
        else if (!strcmp(argv[i],"-T"))
        {
            beta = atof(argv[++i]);
        }
        else if (!strcmp(argv[i],"-N"))
        {
            nIter = atoi(argv[++i]);
        }
        else if (!strcmp(argv[i],"-fs"))
        {
            fSamp = atoi(argv[++i]);
        }
        //else if (!strcmp(argv[i],"-s")) Seed = (unsigned long)atoi(argv[++i]);
    }
    //Toma del tiempo desde el epoch una semilla para el PRNG
    srand(time(NULL));
    int F[L][L];  //Red
    // Corrida de MC
    for (c = 0; c < nIter; c++)
    {
        //Por cada paso de MC hace N giros de spin
        for (a = 0;a < N;a++)
        {
            i = (int)rand() % L;
            j = (int)rand() % L;
            de = E(F,L,i,j); //Delta de Enegía al posiblemente cambiar el estado
            p = exp(de*beta); //Probabilidad de aceptación
            r = rand() % 1; //Número aleatorio entre 0 y 1
            if (p > r)
            {   //Acepta cambio -> lo invierte
                F[i][j] *= -1;
            }
        }
        /* Sample and accumulate averages */
        if ((c % fSamp) == 0 && c > nTerm)
        {
            sampleo(F, L, &m, &e);
            //fprintf(stdout,"%i %.5le %.5le\n", c, s, e);
            //fflush(stdout);
            msum += m;
            esum+=e;
            nSamp++;
        }
    }
    //sprintf(str, "\n%.5lf %.5lf %.5lf\n", beta, ssum/nSamp, esum/nSamp);
    printf("%.5lf %.5lf %.5lf \n", beta, msum/nSamp, esum/nSamp);
    //fputs(str,fp);
    //fclose(fp);
    //fprintf(stdout,"# The average magnetization is %.5lf\n",ssum/nSamp);
    //fprintf(stdout,"# The average energy per spin is %.5lf\n",esum/nSamp);
    //fprintf(stdout,"# Program ends.\n");
}