HILFE! Molekulardynamik-Programm!
-
Liebe Mitglieder,
momentan sitze ich dabei ein Molekulardynamik-Programm in C zu schreiben. Leider weiss ich dabei nicht mehr weiter. Ich hoffe jemand kennt sich hier damit aus. Das Programm soll dazu dienen die Wechselwirkung von 125 Teilchen zu simulieren. Dafür verwende ich den Leap-Frog-Algorithmus. Irgendwie finde ich nicht den Fehler :(. Ich wäre über jede Hilfe dankbar.
Der Code sieht wie folgt aus:
/* Lennard-Jones */ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> int main() { int N = 125, t; //Teilchenzahl double dt = 0.001, rho = 0.02; //Zeitschritte, Dichte int h, i, j, k, l, m, Index; double rx[N], ry[N], rz[N]; //Ortsvektoren double vx[N], vy[N], vz[N]; //Geschwindigkeit double L3, L, rcut, rcut2; //Kastenlänge, Abschneideradius double rxij, rijsq, ryij, rzij; //Abstände double uges=0; //Potential double fxges[N], fyges[N], fzges[N]; //Kraft double vscale=5, b, scale; //Skalierung double o, G[3]; FILE *datei; datei = fopen("daten.txt","a+"); L3 = N/rho; L = pow(L3, 1.0/3.0); rcut = 0.5*L; rcut2 = rcut * rcut; srand((unsigned)time(NULL)); //Initialisiert den Random-Generator mit der akt. Zeit in Sek for(i=0;i<N;i++) { fxges[i]=0; fyges[i]=0; fzges[i]=0; } //Startkonfiguration der Geschwindigkeit for(i=0;i<N;i++) { vx[i]=(2.0*(double)rand()/(double)RAND_MAX-1.0)*vscale; vy[i]=(2.0*(double)rand()/(double)RAND_MAX-1.0)*vscale; vz[i]=(2.0*(double)rand()/(double)RAND_MAX-1.0)*vscale; } //Startkonfiguration der Ortskoordinaten o = pow ( N, 1.0/3.0); h = round(o); b = L/5; scale = b; //scale=1/((h-1.0)+0.5); G[0]=0; G[1]=0; G[2]=0; for(k=0;k<h;k++) { for(l=0;l<h;l++) { for(m=0;m<h;m++) { Index=k+h*l+h*h*m; rx[Index]=(k+G[0])*scale; ry[Index]=(l+G[1])*scale; rz[Index]=(m+G[2])*scale; } } } //Durchlauf für 1000 Zeitschritte for(t=1;t<1000;t++) { for(i=0;i<N-1;i++) { for(j=i+1;j<N;j++) { rxij = rx[i] - rx[j]; rxij = fabsl(rxij - L*round(rxij/L)); if(rxij<rcut) { ryij = ry[i] - ry[j]; ryij = fabsl(ryij - L*round(ryij/L)); rijsq = rxij*rxij+ryij*ryij; if(rijsq<rcut2) { rzij = rz[i] - rz[j]; rzij = fabsl(rzij - L*round(rzij/L)); rijsq = rijsq + rzij*rzij; if(rijsq<rcut2) { fxges[i] = fxges[i] + 48*rxij*(pow(1/rijsq,7) - 0.5*pow(1/rijsq,4)); fyges[i] = fyges[i] + 48*ryij*(pow(1/rijsq,7) - 0.5*pow(1/rijsq,4)); fzges[i] = fzges[i] + 48*rzij*(pow(1/rijsq,7) - 0.5*pow(1/rijsq,4)); fxges[j] = fxges[j] - 48*rxij*(pow(1/rijsq,7) - 0.5*pow(1/rijsq,4)); fyges[j] = fyges[j] - 48*ryij*(pow(1/rijsq,7) - 0.5*pow(1/rijsq,4)); fzges[j] = fzges[j] - 48*rzij*(pow(1/rijsq,7) - 0.5*pow(1/rijsq,4)); uges = uges + 4*(pow(1/rijsq,6)-pow(1/rijsq,3)); } } } } vx[i] = vx[i] + dt*fxges[i]; vy[i] = vy[i] + dt*fyges[i]; vz[i] = vz[i] + dt*fzges[i]; } for(i=0;i<N;i++) { rx[i] = rx[i] +dt*vx[i]; ry[i] = ry[i] +dt*vy[i]; rz[i] = rz[i] +dt*vz[i]; } fprintf(datei,"%d %f\n",t,uges/N); uges = 0; } return 0; }Mein Problem bei diesem Programm ist, dass es läuft, aber nicht so wie es soll. Meine Erwartungen sind, dass die Variable uges sich einem Wert ungleich von Null annähern soll. Leider passiert, dies in diesem Programm nicht. Uges geht allmählich gegen Null. Weiterhin sollte uges zu Beginn exponentiell abfallen zu dem erwarteten Wert. Laut Theorie sollte es bei uges/N bei ungefähr -0.4 liegen. Ich habe mich schon durch dutzende Molekulardynamik Bücher gearbeitet und der Aufbau der Codes entspricht allen Mustern, die ich bisher gesehen habe. Meine Vermutung ist einfach, dass ich irgendeinen dummen Fehler beim Code gemacht habe. Ich habe mich an diese Forum gewendet um sicher zu sein, dass der Code einwandfrei ist und vielleicht den ein oder anderen anzutreffen, der sich mit sowas schonmal auseinander gesetzt hat.
Viele Grüße
euer Walter White
-
p ist nicht deklariert.
Lies dir bitte unbedingt mal den ersten Link in meiner Signatur durch. Deine Frage macht es so gut wie unmöglich, dir hilfreich zu antworten.
-
Der Grund könnten deine nicht initialisierten Variable/Arrays sein.
Deine rx/ry/rz Arrays sehen mir da sehr verdächtig aus ...double rx[N] = {0}, ry[N] = {0}, rz[N] = {0}; // initialisiere Arrays mit 0 ...
-
Also p ist nicht deklariert, da es ein alter Code ist. Sorry! Ich hab es geändert. rx und co müssten alle initialisiert sein. Darauf habe ich auch mal alle Variablen getestet. Zu Beginn funktioniert auch alles, aber irgendwie spuckt mir das Programm viel zu riesige Werte aus, die gar nicht auftreten dürfen und der Abfall zur Null wie erwähnt darf auch nicht sein. Ich glaube, dass ich irgendeine Funktion falsch verwende. Eigentlich sollte an der Stelle der round()-Funktion die ANINT()-Funktion von Fortran. Ich hoffe beide Funktionen sind als gleichwertig anzusehen. Oder denke ich da falsch?
-
Jau jau, so wird das wohl sein - hilft jetzt auch nicht wirklich weiter, ne. :p
-
Vielleicht solltest du mal versuchen in worten zu erklären was der code tun soll. Einfach nur zwei seiten code und "es soll aber nicht 0 rauskommen" ist eine magere fehlerbeschreibung.
-
Also der Code simuliert 125 Teilchen in einem Box mit der Länge L. Die Wechselwirkung zwischen den Teilchen wird durch das Lennard-Jones-Potential beschrieben. Zu Beginn des Codes weise jedem Teilchen eine zufällige Geschwindigkeit und einen Ort zu. Dabei werden die Teilchen so verteilt, dass in jede Richtung den gleichen Abstand haben und zwar L/5. Die Box wird nicht ganz ausgefüllt, da auch in dem Programm virtuelle Teilchen simuliert werden, die um die Box liegen. Dafür ist auch die Zeile mit dem
ryijmin = fabs(ryij - L*lround(ryij/L));nun startet das Programm die Wechselwirkung zwischen Teilchen i zu alle anderen berechnet. Nachdem dies erfolgt ist, werden alle Teilchen nach einem bestimmte Algorithmus (hier: leapfrog) bewegt. Danach fängt alles von vorne an und zwar t-mal.
-
ein oder anderen anzutreffen, der sich mit sowas schonmal auseinander gesetzt hat
Leider haben die wenigsten sich damit auseinandergesetzt, mich eingeschlossen.
-
knivil schrieb:
ein oder anderen anzutreffen, der sich mit sowas schonmal auseinander gesetzt hat
Leider haben die wenigsten sich damit auseinandergesetzt, mich eingeschlossen.
Ich schon, jedoch ist die Fehlerbeschreibung absolut ungenügend. Es ist nicht einmal ein compilierender Code da. Aber nachdem ich das dem Threadersteller gesagt habe und sich nichts geändert habe, werde ich ihm bestimmt nicht 5 Seiten lang alles was ich wissen muss einzeln aus der Nase ziehen. Ich denke, so geht es vielen Leuten.
-
Ich meinte eher, dass es sich wahrscheinlich um einen inhaltliche Fehler handelt. D.h. Leapfrog und Lennard-Jones sind den wenigsten ein Begriff.
-
Irgendwie fehlt mir was bei der Anpassung der Geschwindigkeit. Du machst das einmal für jedes i, aber nur für das Teilchen i, nicht für j. Das heißt, da dein i nur von 0 bis N-2 läuft, dass das letzte Teilchen nie seine Geschwindigkeit ändert. Ist das so richtig?
Hast du den Algorithmus mal irgendwo in Pseudocode oder einer anderen Sprache oder so zum Vergleich?
-
knivil schrieb:
Ich meinte eher, dass es sich wahrscheinlich um einen inhaltliche Fehler handelt. D.h. Leapfrog und Lennard-Jones sind den wenigsten ein Begriff.
Wieso? Hier sind schon einige Leute anwesend, die von numerischer Simulation eine Ahnung haben.
@OP: Dir kann man so nur den Allgemeinsten aller Ratschläge geben. Schnapp Dir einen Debugger, schau Dir an, was das Programm macht, und dann versuche zu verstehen, was schief läuft. Mach Dir Kommentare in den Code und versuche vor allem mal ein einfaches Beispiel mit ein paar Freiheitsgraden auf einem Zettel durchzurechnen. Das hilft ungemein beim debuggen und auch für zukünftige Aufgaben dieser Art. Wenn Du Hilfe benötigst, dann komme mit konkreteren Fragen. Du hälst ja auch sicher eher jemandem die Tür auf, als ihm die Wocheneinkäufe die Treppe hoch zu tragen, oder?
-
Hast Du mal getestet, ob das Verhalten von der Groesse Deiner Zeitschritte abhaengt?
EDIT: Wenn sich uges immer weiter dem Wert 0 annaehert, dann heisst das, dass Deine Teilchen immer weiter auseinanderlaufen und die "Rueckprojektion" in Deine Box irgendwie nicht klappt. Ueberpruef mal, wie sich Deine Abstandsquadrate rijsq ueber die Zeit entwickeln. Die sollten ja eigentlich im Grossen und Ganzen nicht wesentlich groesser werden.
