/* (Quantenmechanische) MC-Simulation eines 'Quantendots' Numerische Methoden der Vielteilchenphysik WS01/02 Roland Winkler und Guenter Zwicknagel */ #include #include #include #define PI 3.14159265359 int np, anzmon; double omega2, coul, temp; double u_eff (double **partx, double **party); double upot (double *partx, double *party); void force (double *x, double *y, double **f); double dforce (); double ham_op (double **partx, double **party, double **f); float ran1 (int *idum); float gasdev (int *idum); void mswap (double ***x, double ***y); double **dmatrix(int nrl, int nrh, int ncl, int nch); /* #define ZUFALL (ran1 (&idum) - 0.5) */ #define ZUFALL gasdev (&idum) /*****************/ /* Hauptprogramm */ /*****************/ int main (int argc, char *argv[]) { extern int np, anzmon; extern double omega2, coul, temp; int anzges, idum, i, j, k, schritt, akzep; double zeta, rskala, upot_alt, upot_neu, du, ham, ham_sum; double **partx, **party, **xn, **yn, **frc; int startsamp = 20000, anzmess = 1000; FILE *daten_file; if (argc == 7) { sscanf (argv [1], "%lf", &omega2); /* omega2 */ sscanf (argv [2], "%lf", &coul); /* coul */ sscanf (argv [3], "%lf", &temp); /* Temperatur */ sscanf (argv [4], "%i", &np); /* Teilchenzahl */ sscanf (argv [5], "%i", &anzmon); /* Zahl der Monomere */ sscanf (argv [6], "%i", &anzges); /* Zahl der Schritte */ } else { printf ("%s: Wrong number of arguments\n", argv [0]); exit (1); } /* Felder alloziieren */ partx = dmatrix (0, anzmon, 1, np); party = dmatrix (0, anzmon, 1, np); xn = dmatrix (0, anzmon, 1, np); yn = dmatrix (0, anzmon, 1, np); frc = dmatrix (1, np, 1, 2); /* Hilsfeld fuer Erwartungswert */ /* Die charakteristische Laengenskala */ rskala = 4.0 / sqrt (omega2) + 5.0 * pow (coul/omega2, 1.0/3.0); /* Random-Generator initialisieren */ idum = - fabs (rskala*np*temp + 123.) - 0.5; printf ("Initialisierung: %f\n", ran1 (&idum)); /* Teilchenorte initialisieren */ zeta = 0.3 * rskala; for (j = 1; j <= np; j++) { partx [0][j] = zeta * ZUFALL; party [0][j] = zeta * ZUFALL; } /* Idee: Die einzelnen Monomere liegen im Phasenraum eng zusammen */ zeta = 0.01 * zeta; for (i = 1; i <= anzmon; i++) { for (j = 1; j <= np; j++) { partx [i][j] = partx [0][j] + zeta * ZUFALL; party [i][j] = party [0][j] + zeta * ZUFALL; } } /* Berechnung des Potentials fuer die Anfangskonfiguration */ upot_alt = u_eff (partx, party); ham = ham_op (partx, party, frc); printf ("Startpotential %e\n", ham); /* Ausgabe initialisieren */ daten_file = fopen ("datqdclass.dat", "w"); fprintf (daten_file, "# omega2: %f\n", omega2); fprintf (daten_file, "# coul: %f\n", coul); fprintf (daten_file, "# temp: %f\n", temp); fprintf (daten_file, "# Teilchenzahl: %i\n", np); fprintf (daten_file, "# Zahl der Monomere: %i\n", anzmon); fprintf (daten_file, "# Zahl der Schritte: %i\n", anzges); akzep = 0; zeta = 0.02 * rskala; ham_sum = 0.0; k = 0; /* Beginn des Schleifenprogramms */ for (schritt = 1; schritt <= anzges; schritt++) { /* Erzeugen einer neuen Konfiguration */ /* Idee: Variiere abwechselnd immer nur ein Monomer */ for (i = 0; i <= anzmon; i++) { for (j = 1; j <= np; j++) { xn [i][j] = partx [i][j]; yn [i][j] = party [i][j]; } } k = (k+1) % (anzmon+1); for (j = 1; j <= np; j++) { xn [k][j] = partx [k][j] + zeta * ZUFALL; yn [k][j] = party [k][j] + zeta * ZUFALL; } /* Berechnung der Aenderung der pot. Energie */ upot_neu = u_eff (xn, yn); du = upot_neu - upot_alt; /* Befrage das Orakel */ if (exp (-du / temp) > ran1 (&idum)) { akzep ++; upot_alt = upot_neu; mswap (&partx, &xn); mswap (&party, &yn); /* Nur wenn sich die Konfiguration geaendert hat, wird H ausgewertet */ ham = ham_op (partx, party, frc); } /* Observable bestimmen */ if (schritt > startsamp) { ham_sum += ham; } if ((schritt % anzmess) == 0) { /* Datenausgabe */ if (schritt > startsamp) { fprintf (daten_file, "%i %13.6e %13.6e %13.6f\n", schritt, ham_sum / (np * (schritt - startsamp)), (double) akzep / anzmess, zeta); } /* Anpassung des zeta-Parameters */ if ((double) akzep / anzmess > 0.5) zeta *= 1.05; else zeta *= 0.95; akzep = 0; } } printf ("Letztes Potential %e\n", ham); return 0; } /************************/ /* Effektives Potential */ /************************/ double u_eff (double **partx, double **party) { extern int np, anzmon; extern double temp; double u, upoly, rx, ry, efak, pfak; int i, j; /* Einmal haben wir immer das echte Potential */ pfak = 1.0 / (anzmon+1); u = pfak * upot (partx [0], party[0]); if (anzmon > 0) { efak = 0.5 * (anzmon+1) * temp * temp; for (i = 1; i <= anzmon; i++) { /* Harmonische Polymer-Wechselwirkung */ upoly = 0.0; for (j = 1; j <= np; j++) { rx = partx [i][j] - partx [i-1][j]; ry = party [i][j] - party [i-1][j]; upoly += rx * rx + ry * ry; } u += efak * upoly + pfak * upot (partx [i], party[i]); } /* Ring-Polymer */ upoly = 0.0; for (j = 1; j <= np; j++) { rx = partx [0][j] - partx [anzmon][j]; ry = party [0][j] - party [anzmon][j]; upoly += rx * rx + ry * ry; } u += efak * upoly; } return u; } /************************/ /* Echte Wechselwirkung */ /************************/ double upot (double *partx, double *party) { extern int np; extern double omega2, coul; double u = 0.; int i, j; /* harmonischer Anteil */ for (i = 1; i <= np; i++) u += partx [i] * partx [i] + party [i] * party [i]; u *= 0.5 * omega2; /* Coulomb-Anteil */ if (coul != 0.0) { for (i = 1; i <= np-1; i++) { for (j = i+1; j <= np; j++) { double rx = partx [i] - partx [j]; double ry = party [i] - party [j]; u += coul / sqrt (rx*rx + ry*ry); } } } return u; } /*********************/ /* Hamilton-Operator */ /*********************/ double ham_op (double **partx, double **party, double **frc) { extern int np, anzmon; extern double temp; double fak = (anzmon+1) * temp; double ham = 0.0; int i; if (anzmon > 0) { /* Kraefte aus Polymer-Wechselwirkung */ for (i = 1; i <= np; i++) { frc [i][1] = fak * (partx [0][i] - partx [1][i]); frc [i][2] = fak * (party [0][i] - party [1][i]); } /* Kraefte der echten Wechselwirkung draufaddieren */ force (partx [0], party [0], frc); /* Quadrieren und aufsummieren */ for (i = 1; i <= np; i++) ham += frc [i][1] * frc [i][1] + frc [i][2] * frc [i][2]; } ham = pow ((anzmon+1) * temp / (4 * PI * PI), (double) np * (anzmon+1)) * (- 0.5 * (ham + np * ( - 2.0 * fak - dforce () / fak)) + upot (partx [0], party [0])); return ham; } /***********/ /* Kraefte */ /***********/ void force (double *x, double *y, double **frc) { extern int np, anzmon; extern double omega2, coul; int i, j; double fx, fy, rx, ry, sq, den; double ofak = omega2 * (anzmon+1) * temp; double cfak = coul / ((anzmon+1) * temp); /* Kraefte beim harmonischen Oszillator */ for (i = 1; i <= np; i++) { frc [i][1] += ofak * x[i]; frc [i][2] += ofak * y[i]; } /* Kraefte aus Coulomb-Wechselwirkung */ if (coul != 0.0) { for (i = 1; i <= np; i++) { for (j = 1; j < i; j++) { rx = x [i] - x [j]; ry = y [i] - y [j]; sq = sqrt (rx * rx + ry * ry); den = cfak / (sq * sq * sq); fx = rx * den; fy = ry * den; /* Actio = Reactio */ frc [i][1] -= fx; frc [i][2] -= fy; frc [j][1] += fx; frc [j][2] += fy; } } } return; } /***********************************/ /* Zweite Ableitung des Potentials */ /***********************************/ double dforce () { extern double omega2; return 2.0 * omega2; } /*****************/ /* Swap matrices */ /*****************/ void mswap (double ***x, double ***y) { double **temp; temp = *x; *x = *y; *y = temp; } /***********************************************************************/ /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ /***********************************************************************/ #define NR_END 1 double **dmatrix(int nrl, int nrh, int ncl, int nch) { int i, nrow=nrh-nrl+1,ncol=nch-ncl+1; double **m; /* allocate pointers to rows */ m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*))); m += NR_END; m -= nrl; /* allocate rows and set pointers to them */ m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double))); m[nrl] += NR_END; m[nrl] -= ncl; for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; /* return pointer to array of pointers to rows */ return m; } #undef NR_END /**************************************************/ /* Zufallszahlengenerator aus "Numercial Recipes" */ /**************************************************/ #define M1 259200 #define IA1 7141 #define IC1 54773 #define RM1 (1.0/M1) #define M2 134456 #define IA2 8121 #define IC2 28411 #define RM2 (1.0/M2) #define M3 243000 #define IA3 4561 #define IC3 51349 float ran1(int *idum) { static long ix1,ix2,ix3; static float r[98]; float temp; static int iff=0; int j; if (*idum < 0 || iff == 0) { iff=1; ix1=(IC1-(*idum)) % M1; ix1=(IA1*ix1+IC1) % M1; ix2=ix1 % M2; ix1=(IA1*ix1+IC1) % M1; ix3=ix1 % M3; for (j=1;j<=97;j++) { ix1=(IA1*ix1+IC1) % M1; ix2=(IA2*ix2+IC2) % M2; r[j]=(ix1+ix2*RM2)*RM1; } *idum=1; } ix1=(IA1*ix1+IC1) % M1; ix2=(IA2*ix2+IC2) % M2; ix3=(IA3*ix3+IC3) % M3; j=1 + ((97*ix3)/M3); temp=r[j]; r[j]=(ix1+ix2*RM2)*RM1; return temp; } #undef M1 #undef IA1 #undef IC1 #undef RM1 #undef M2 #undef IA2 #undef IC2 #undef RM2 #undef M3 #undef IA3 #undef IC3 /*********************************/ /* Gauss-verteilte Zufallszahlen */ /*********************************/ float gasdev (int *idum) { static int iset=0; static float gset; float fac,r,v1,v2; if (iset == 0) { do { v1=2.0*ran1(idum)-1.0; v2=2.0*ran1(idum)-1.0; r=v1*v1+v2*v2; } while (r >= 1.0); fac=sqrt(-2.0*log(r)/r); gset=v1*fac; iset=1; return v2*fac; } else { iset=0; return gset; } }