1 / 38

Tipico schema in un codice

Tipico schema in un codice. fornite dall'esterno (file contenente una configurazione). posizioni iniziali. create dal codice stesso in base alle richieste dello user ("crea" cristallo fcc, bcc, diamante, etc.). Tipico schema in un codice.

balthasar-lali
Download Presentation

Tipico schema in un codice

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Tipico schema in un codice fornite dall'esterno (file contenente una configurazione) posizioni iniziali create dal codice stesso in base alle richieste dello user ("crea" cristallo fcc, bcc, diamante, etc.)

  2. Tipico schema in un codice x(1),y(1),z(1),x(2),y(2),z(2), ..,x(natoms),y(natoms),z(natoms) posizioni iniziali Gli atomi sono sempre numerati in qualche modo. La particolare numerazione non è di per sè significativa 2 3 1 4 6 7 5 8 9 10 12 11

  3. Tipico schema in un codice lista dei vicini: voglio sapere per ogni atomo i quanti - nvicini(i) - atomi rientrano nel suo raggio di cutoff configurazione nota Errore gravissimo e comune: cerco i vicini di 6 e trovo (se metto un cutoff a primi vicini) 2,5,7,10. Totale: nvicini(6)=4. Se per esempio devo fare un ciclo in cui calcolo la distanza tra 6 e i suoi vicini, non devo fare do j=1,nvicini(6) dx=x(6)-x(j) .....!!!!calcolerei le distanze tra 6 e 1, 6 e 2, etc. Invece devo fare: dx=x(6)-x(ivicini(6,j)) dove ivicini(6,1)=2; ivicini(6,2)=5, etc. 2 3 1 4 E' chiave la matrice intera ivicini(i,j) che mi dice che il vicino j-simo dell'atomo i è l'atomo ivicini(i,j) 6 7 5 8 9 10 12 11

  4. Implementazione calcolo energia

  5. Metodo diretto Notare che potremmo fare il ciclo in metà del tempo ... do i=1,natoms do j=1,natoms rij=(x(i)-x(j))*(x(i)-x(j))+(y(i) .... rij=sqrt(rij) if(i.ne.j)then ene=ene+..... endif end do end do

  6. Metodo basato sulla lista dei vicini call vicini(natoms,nvicmax,nvicini, ivicini) do i=1,natoms do k=1,nvicini(i) j=ivicini(i,k) rij=(x(i)-x(j))*(x(i)-x(j))+(y(i) .... rij=sqrt(rij) ene=ene+..... end do end do

  7. Se introduco il raccordo polinomiale ... call vicini(natoms,nvicmax,nvicini, ivicini) do i=1,natoms do k=1,nvicini(i) j=ivicini(i,k) rij=(x(i)-x(j))*(x(i)-x(j))+(y(i) .... rij=sqrt(rij) ecoppia=0.d0 if(rij.le.rprimo)then ecoppia=4.d0*eps ….. elseif(rij.lt.rc)then ecoppia=polinomio endif epot=epot+ecoppia end do end do epot=0.5d0*epot

  8. Notare che detto così non è molto furbo call vicini(natoms,nvicmax,nvicini, ivicini) do i=1,natoms do k=1,nvicini(i) j=ivicini(i,k) rij=(x(i)-x(j))*(x(i)-x(j))+(y(i) .... rij=sqrt(rij) ene=ene+..... end do end do doppio ciclo (natoms x natoms) doppio ciclo (natoms x nvicini) ... a menoche non sipossaevitarediricalcolare la listadeivicini

  9. Calcolo delle forze: calcoliamo la forza lungo x che agisce sull'atomo k forza su x che agisce sull'atomo k per effetto dell'atomo j

  10. Calcolo delle forze: calcoliamo la forza lungo x che agisce sull'atomo k, dovuta all'atomo j

  11. Calcolo delle forze: calcoliamo la forza che agisce sull'atomo k, dovuta all'atomo j (la forzatotale su k è subito ottenuta sommando su j) E' così facile perchè è a coppie (e anche perchè la forma funzionale è quella semplicissima di Lennard-Jones)

  12. Calcolo delle forze: come si traduce tutto ciò in un codice? forza su x che agisce sull'atomo k per effetto del suo vicino j-simo Alla fine mi interessa quella totale che agisce su k: basta sommare su j..

  13. Non ci dimentichiamo di due complicazioni Se ho PBC ......... Nel calcolo sia dell'energia sia delle forze devo utilizzare il raccordo (e la sua derivata) invece dell'espressione standard del potenziale!!!

  14. Calcolo delle forze: calcoliamo la forza lungo x che agisce sull'atomo k, dovuta all'atomo j nel raccordo

  15. Calcolo delle forze: calcoliamo la forza lungo x che agisce sull'atomo k, dovuta all'atomo j nel raccordo

  16. r’<rc Dopo un po’ di conti …………………..

  17. Calcolo delle forze: calcoliamo numericamente la forza lungo x che agisce sull'atomo k

  18. Calcolo delle forze: calcoliamo numericamente la forza lungo x che agisce sull'atomo k • Nel codice: • salvo le posizione ufficiali in vettori xsave,ysave,zsave • 2) x(k)=x(k)+h; calcolo l’energia e la salvo in eplus • 3) x(k)=xsave(k)-h; calcolo l’energia e la salvo in eminus • 4)fx(k)=-(eplus-eminus)/2h • Facciamolo nel caso dell’atomo che si allontana dalla superficie, al variare di z, confrontando forze analitiche e numeriche.

  19. Il calcolo delle forze permette immediatamente di effettuare simulazioni di STATICA molecolare Statica molecolare: configurazioni (locali) di minima energia U

  20. Calcolo delle forze: calcoliamo numericamente la forza lungo x che agisce sull'atomo k Statica molecolare: configurazioni (locali) di minima energia U Mi muovo piano piano lungo la direzione della forza

  21. Calcolo delle forze: calcoliamo numericamente la forza lungo x che agisce sull'atomo k Statica molecolare: configurazioni (locali) di minima energia U Mi muovo piano piano lungo la direzione della forza

  22. Calcolo delle forze: calcoliamo numericamente la forza lungo x che agisce sull'atomo k Statica molecolare: configurazioni (locali) di minima energia U Mi muovo piano piano lungo la direzione della forza

  23. Calcolo delle forze: calcoliamo numericamente la forza lungo x che agisce sull'atomo k Statica molecolare: configurazioni (locali) di minima energia U Mi muovo piano piano lungo la direzione della forza. Sono nel minimo!

  24. In due dimensioni … mi muovo in x lungo fx, in y lungo fy V(x,y) Metodo della discesa più rapida (“steepest descent”)

  25. Schema implementazione steepest descent in tante dimensioni Versione semplificata fmax = max(|fx1|,......|fzN|) calcolata in force Problema: Se ho forze grandi devo muovermi di poco se no la direzione scelta è poco rappresentativa! C troppo grande: non convergo, e/o oscillo, e/o il sistema “esplode” C troppo piccolo: ci metto troppo Devo variare C ...

  26. Salvo posizioni in xsave

  27. Calcolo delle forze: calcoliamo numericamente la forza lungo x che agisce sull'atomo k Criteri di convergenza L’energia varia di poco tra un passo e l’altro (livello 0); La norma della forza è al di sotto di una certa soglia (livello 1); La forza massima su un atomo è sotto ad una certa soglia (livello 2); Devono valere tutti e tre assieme (livello 3). Forze piccole? < 0.001 eV/Å

  28. Calcolo delle forze: calcoliamo numericamente la forza lungo x che agisce sull'atomo k Occhio! Lo steepest descent è un metodo che permette di trovare solo i minimi locali!!!!

  29. Example: adatom positions on a fcc(001) Surface. Top view.

  30. Let us “deposit” few adatoms

  31. Let us find the closest minimum configuration

  32. Let us find the closest minimum configuration

  33. Let us find the closest minimum configuration

  34. Let us find the closest minimum configuration

  35. Let us find the closest minimum configuration

  36. Let us find the closest minimum configuration (Also: surface relaxations!!)

More Related