Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- /*
- 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");
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement