J'ai calculé les valeurs propres de l'hamiltonien pour l'atome d'hydrogène 1D en unités atomiques avec la méthode de Fourier-Grid-Hamiltonian dans un joli petit programme Fortran.

Toutes les valeurs propres trouvées entre -1 et 0 (les états liés) sont enregistrées dans un fichier ligne par ligne comme ceci:

-0.50016671392950229
-0.18026105614262633
-0.11485673263086937
-4.7309305955423042E-002
-4.7077108902158216E-002

Comme le nombre de valeurs propres trouvées varie en fonction de la taille de pas utilisée par mon programme, le nombre d'entrées dans le fichier peut varier (en théorie, il y en a une infinie).

Je veux maintenant tracer les valeurs du fichier sous la forme d'une ligne parallèle à l'axe des x avec le décalage donné par les valeurs lues à partir du fichier. Je veux aussi pouvoir tracer les données uniquement jusqu'à un certain numéro de ligne, car les valeurs se rapprochent vraiment les unes des autres à mesure que vous vous éloignez de zéro et qu'elles ne peuvent plus être distinguées à l'œil nu.

(Ici, par exemple, il serait logique de tracer les quatre premières entrées, la cinquième est déjà trop proche de la précédente)

Je sais que l'on peut tracer des lignes parallèles à l'axe x avec la commande plot *offset* mais je ne sais pas comment dire à gnuplot d'utiliser les données du fichier. Jusqu'à présent, j'ai dû tracer manuellement les valeurs.

Dans un deuxième temps, je voudrais tracer les données uniquement dans une certaine plage x, plus concrète entre les points d'intersection avec le potentiel harmonique utilisé pour la solution numérique V(x) = -1/(1+abs(x))

Le résultat devrait ressembler à ceci: schéma du tracé souhaité (sosie)

Le plus proche que je suis arrivé, c'était avec

plot -1/(1+abs(x)),-0.5 title 'E0',-0.18 title 'E1', -0.11 title 'E2'

Ce qui m'a donné le résultat suivant: mon intrigue

J'espère que vous pouvez m'aider, et je suis vraiment curieux de savoir si gnuplot peut réellement faire la deuxième étape que j'ai décrite!

1
Vito 6 avril 2017 à 17:30

2 réponses

Meilleure réponse

Une autre approche peut être d'obtenir des données de "energy.dat" (comme indiqué dans la question) avec les commandes system et cat (donc en supposant un système de type Un * x ...) et de sélectionner V (x) et E à chaque x via max:

set key bottom right
set yr [-1:0.2]
set samples 1000

Edat = system( "cat energy.dat" )

max(a,b) = ( a > b ) ? a : b
V(x) = -1/(1+abs(x))

plot for [ E in Edat ] \
     max(V(x),real(E)) title sprintf("E = %8.6f", real(E)) lw 2, \
     V(x) title "V(x) = -1/(1+|x|)" lc rgb "red" lw 2

enter image description here

Si nous changeons le potentiel en V(x) = -abs(cos(x)), l'intrigue semble assez drôle (et les niveaux d'énergie ne sont bien sûr pas corrects!)

enter image description here


Plus de détails sur le script:

  • max n'est pas une fonction intégrée dans Gnuplot, mais une fonction définie par l'utilisateur ayant deux arguments formels. Ainsi, par exemple, nous pouvons le définir comme

    mymax( p, q ) = ( p > q ) ? p : q

Avec tout autre nom (et utilisez mymax dans la commande plot). Ensuite, le symbole ? est un opérateur ternaire qui donne un notation abrégée pour une construction if ... else. Dans un pseudo-code, cela fonctionne comme

function max( a, b ) {
    if ( a > b ) then
        return a
    else
        return b
    end
}

De cette façon, max(V(x),real(E)) sélectionne la valeur la plus élevée entre V(x) et real(E) pour n'importe quel x et E donné.

  • Ensuite, Edat = system( "cat energy.dat" ) dit à Gnuplot d'exécuter la commande shell "cat energy.dat" et d'affecter la sortie à une nouvelle variable Edat. Dans le cas ci-dessus, Edat devient une chaîne contenant une séquence de valeurs d'énergie lues à partir de "energy.dat". Vous pouvez vérifier le contenu de Edat par print( Edat ). Par exemple, cela peut être quelque chose comme

    Edat = "-0.11 -0.22 ... -0.5002"

  • plot for [ E in Edat ] ... boucle sur les mots contenus dans une chaîne Edat. Dans le cas ci-dessus, E prend une chaîne "-0,11", "-0,22", ..., "-0,5002" une par une. real(E) convertit cette chaîne en une valeur à virgule flottante. Il est utilisé pour passer E (une chaîne de caractères) à n'importe quelle fonction mathématique.

  • L'idée de base est de dessiner un potentiel tronqué au-dessus de E, max (V (x), E), pour chaque valeur de E. (Vous pouvez vérifier la forme d'un tel potentiel par plot max(V(x),-0.5), par exemple). Après avoir tracé de telles courbes, nous redessinons le potentiel V(x) pour le faire apparaître comme une seule courbe de potentiel avec une couleur différente.

  • set samples 1000 augmente la résolution du tracé de 1 000 points par courbe. 1000 est arbitraire, mais cela semble être suffisant pour rendre le chiffre assez lisse.

1
roygvib 7 avril 2017 à 20:35

Comme pour la première partie de votre question, vous pouvez par exemple utiliser le style de traçage xerrorbars comme:

set terminal pngcairo
set output 'fig.png'

unset key
set xr [-1:1]
set yr [-1:0]

unset bars
plot '-' u (0):($1<-0.1?$1:1/0):(1) w xerrorbars pt 0 lc rgb 'red'
-0.50016671392950229
-0.18026105614262633
-0.11485673263086937
-4.7309305955423042E-002
-4.7077108902158216E-002
e

L'idée ici est de:

  1. interpréter les énergies E comme des points de coordonnées (0, E) et attribuer à chacun d'eux une barre d'erreur x de largeur 1 (via la troisième partie de la spécification (0):($1<-0.1?$1:1/0):(1))
  2. "simuler" les lignes horizontales avec x-errorbars. À cette fin, unset bars et pt 0 s'assurent que Gnuplot n'affiche que des lignes simples.
  3. ne considérer que les énergies E<-0.1, les expressions $1<-0.1?$1:1/0 s'évalue autrement à une valeur indéfinie 1/0 ce qui a pour conséquence que rien n'est tracé pour un tel E.
  4. plot '-' avec des valeurs explicites peut bien sûr être remplacé par, par exemple, plot 'your_file.dat'

Cela produit: entrez la description de l'image ici

Pour la deuxième partie, cela dépend principalement de la complexité de votre fonction V(x). Dans le cas particulier de V(x)=-1/(1+|x|), on pourrait déduire directement qu'il est symétrique autour de x=0 et calculer les points de retournement explicitement, par exemple,

set terminal pngcairo
set output 'fig.png'

fName = 'test.dat'

unset key
set xr [-10:10]
set yr [-1:0]

unset bars

f(x) = -1 / (1+abs(x))
g(y) = (-1/y - 1)

plot \
    f(x) w l lc rgb 'black', \
    fName u (0):($1<-0.1?$1:1/0):(g($1)) w xerrorbars pt 0 lc rgb 'red', \
    fName u (0):($1<-0.1?$1:1/0):(sprintf("E%d", $0)) w labels offset 0, char 0.75

Qui donne

enter image description here

L'idée est fondamentalement la même qu'avant, seule la largeur de la barre d'erreur dépend maintenant de la coordonnée y (l'énergie). De plus, le style labels est utilisé pour produire des étiquettes explicites.

2
ewcz 6 avril 2017 à 17:02