Quantum Mechanics simulation using Finite Elements

This post is about an unusual kind of simulation using Code_Aster: a quantum mechanics system. In this case we will try to calculate the probability of  finding the electron somewhere in a hydrogen atom.

This means that we need to find the wave function of the electron in a potential field, by solving this problem (mathematically speaking it is an eigenfunction problem):

equation
Schrödinger equation

Since we need to solve a differential equation everywhere in the space domain and we would like to have more precision where the potential is stronger, we can try to use finite elements for that. This is not the usual way one would do this calculation in Physics or in Computational Chemistry, but it still works.

Why is this important ? Well, this problem has many, many practical applications, from  material engineering, microelectronics, protein folding, etc.

This problem is obviously a very simple one: in the hydrogen atom there is only one electron, so we can consider the potential as constant.

The expected results, from the Wikipedia page above, are this kind of shapes:

Hydrogen_Density_Plots.png
Wavefunctions of the electron in a hydrogen atom

Here is what we can get from Code_Aster, using a refined mesh and 5 minutes of computing on a laptop, you should be able to recognize some of the shapes above:

mode1clipmode1contour

mode3clip.png

mode6contourmode13clip

mode15contour.png

mode20clip.png

For a more detailed description of what is going on here, you can refer to this Code_Aster test case: Frequencies of a one-dimensional quantum harmonic oscillator

Let’s start with the mesh: we need to consider the whole space, since we do not really know in advance where the electron might be. We can use a sphere, centered around the atom nucleus. We can easily do this with Salome, this time by using a script (you will need some Python knowledge for this, but it is not that complicated):

# -*- coding: iso-8859-1 -*-

import sys
import salome
import inspect
import os

salome.salome_init()
theStudy = salome.myStudy

import salome_notebook
notebook = salome_notebook.NoteBook(theStudy)
scriptfile = inspect.getframeinfo(inspect.currentframe()).filename
scriptfolder = os.path.dirname(os.path.abspath(scriptfile))

sys.path.insert( 0, scriptfolder)

###
### GEOM component
###

import GEOM
from salome.geom import geomBuilder
import math
import SALOMEDS

geompy = geomBuilder.New(theStudy)

Sphere_1 = geompy.MakeSphereR(1)
geompy.addToStudy( Sphere_1, 'Sphere_1' )

###
### SMESH component
###

import SMESH, SALOMEDS
from salome.smesh import smeshBuilder

smesh = smeshBuilder.New(theStudy)
Mesh_1 = smesh.Mesh(Sphere_1)
NETGEN_2D3D = Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
NETGEN_3D_Parameters = NETGEN_2D3D.Parameters()
NETGEN_3D_Parameters.SetMaxSize( 0.347537 )
NETGEN_3D_Parameters.SetSecondOrder( 1 )
NETGEN_3D_Parameters.SetOptimize( 1 )
NETGEN_3D_Parameters.SetMinSize( 0.00350402 )
NETGEN_3D_Parameters.SetQuadAllowed( 0 )
NETGEN_3D_Parameters.SetFineness( 5 )
NETGEN_3D_Parameters.SetGrowthRate( 0.1 )
NETGEN_3D_Parameters.SetNbSegPerEdge( 15 )
NETGEN_3D_Parameters.SetNbSegPerRadius( 17 )
isDone = Mesh_1.Compute()

## Set names of Mesh objects
smesh.SetName(NETGEN_2D3D.GetAlgorithm(), 'NETGEN_2D3D')
smesh.SetName(NETGEN_3D_Parameters, 'NETGEN 3D Parameters')
smesh.SetName(Mesh_1.GetMesh(), 'Mesh_1')

Mesh_1.ExportMED( os.path.join(scriptfolder, 'Mesh3d10auModerate.med', 0, SMESH.MED_V2_2, 1, None ,1)

if salome.sg.hasDesktop():
    salome.sg.updateObjBrowser(1)

Ok, we can run this script from the GEOM module (File->Load Script) and we get a mesh file named “Mesh3d10auModerate.med”.

In order to save precision, since we are handling some very small numbers, we need to use the Atomic Unit system.

Now, let’s get to the Code_Aster command file, in this case it uses some Python too, you can name this file “hydrogen.comm”:


import math

# Hartree Atomic Units
me = 1.0 # [m_e]
e = 1.0 # [e]
mp = 1837.0 # [m_e]
ht = 1.0 # [hbar]
a0 = 1.0 # [bohr]
eps0 = 1.0/(4.0*math.pi)

k = +(e**2)/(4.0*math.pi*eps0) # potential multiplier
mu=me*mp/(me+mp) # reduced mass
c1=-ht**2/(2*mu) # laplacian multiplier

print "k=%s"%k
print "mu=%s"%mu
print "c1=%s"%c1

DEBUT(PAR_LOT='NON')

MAILLAGE=LIRE_MAILLAGE(FORMAT='MED',)

FPOTV = FORMULE(NOM_PARA=('X','Y','Z'),VALE='k/(sqrt(X**2+Y**2+Z**2)*c1)')

MAILINIT = MAILLAGE


# taille des mailles du maillage grossier initial 
h0 = 0.0005

# taille cible des mailles apres raffinement 
hc = 0.0001

# nombre de niveau de raffinements
# attention, en python log = ln (logarithme neperien)
#n = (log(h0)-log(hc))/log(2)
#n = 7
n = 0
nb_raff = int(n)+1

for index_raff in range(nb_raff):

    CHXY=CREA_CHAMP(OPERATION='EXTR',TYPE_CHAM='NOEU_GEOM_R',NOM_CHAM='GEOMETRIE',MAILLAGE=MAILINIT)
    CHFPOTV = CREA_CHAMP(OPERATION='AFFE',TYPE_CHAM='NOEU_NEUT_F',MAILLAGE=MAILINIT,AFFE= (_F(TOUT='OUI', NOM_CMP='X1',VALE_F=FPOTV,),))
    CHPOTV=CREA_CHAMP(OPERATION='EVAL',TYPE_CHAM='NOEU_NEUT_R',CHAM_F=CHFPOTV, CHAM_PARA=CHXY )
    DETRUIRE(INFO=1,CONCEPT=_F(NOM=(CHXY, CHFPOTV)))
    if index_raff != n:
        MAILRAF = CO('MAIL%d'%index_raff)
        MACR_ADAP_MAIL(ADAPTATION         = 'RAFFINEMENT',
                         CHAM_GD            = CHPOTV,
                         CRIT_RAFF_REL       = 0.02,
                         #CRIT_DERA_REL      = 0.02,
                         USAGE_CMP          = 'ABSOLU', 
                         MAILLAGE_N         = MAILINIT,
                         MAILLAGE_NP1       = MAILRAF,
                         NIVE_MAX=6,
                         )
        
        DETRUIRE(INFO=1,CONCEPT=_F(NOM=(MAILINIT, CHPOTV)))
        MAILINIT = MAILRAF
else:
	MAILRAF = MAILINIT

MAT=DEFI_MATERIAU(THER=_F( RHO_CP = 1.0,LAMBDA = 1.0),)

CHMAT=AFFE_MATERIAU(MAILLAGE=MAILRAF,
                    AFFE=_F(TOUT='OUI',
                            MATER=MAT,),)
                            
LAMBDA_F=DEFI_CONSTANTE( VALE=1.0 )
RHOCP_F = FORMULE(NOM_PARA=('NEUT1',),VALE='NEUT1')
MATV=DEFI_MATERIAU(THER_FO=_F( RHO_CP = RHOCP_F,LAMBDA = LAMBDA_F),)

CHMATV=AFFE_MATERIAU(MAILLAGE=MAILRAF,AFFE=_F(TOUT= 'OUI', MATER= MATV),AFFE_VARC=_F(NOM_VARC='NEUT1', CHAM_GD=CHPOTV,),)

MODE=AFFE_MODELE(MAILLAGE=MAILRAF,
                   AFFE=_F(TOUT='OUI',
                           PHENOMENE='THERMIQUE',
                           MODELISATION='3D',),)

K_ELEM=CALC_MATR_ELEM(MODELE=MODE,
                            CHAM_MATER=CHMAT,
                               OPTION='RIGI_THER',
                                )

MV_ELEM=CALC_MATR_ELEM(      MODELE=MODE,
                            CHAM_MATER=CHMATV,
                               OPTION='MASS_THER',
							   INCR_INST= 1.0,
                                )
                                                        

M_ELEM=CALC_MATR_ELEM(      MODELE=MODE,
                            CHAM_MATER=CHMAT,
                               OPTION='MASS_THER',
							   INCR_INST= 1.0,
                                )

NUM=NUME_DDL(MATR_RIGI=K_ELEM)

K_ASSE=ASSE_MATRICE( MATR_ELEM=K_ELEM,   NUME_DDL=NUM)

MV_ASSE=ASSE_MATRICE( MATR_ELEM=MV_ELEM,   NUME_DDL=NUM)

KMV_ASSE = COMB_MATR_ASSE(COMB_R = ( _F( MATR_ASSE = K_ASSE, COEF_R = 1.),
                                     _F( MATR_ASSE= MV_ASSE, COEF_R = 1.),
                                   )
                         )
                        
DETRUIRE(INFO=1,CONCEPT=_F(NOM=(K_ASSE, MV_ASSE)))                         

M_ASSE=ASSE_MATRICE( MATR_ELEM=M_ELEM,   NUME_DDL=NUM)

#nbmod01 = INFO_MODE(MATR_RIGI=KMV_ASSE,MATR_MASS=M_ASSE,FREQ=(-0.5,0.0)) 

MODES=CALC_MODES(     MATR_RIGI=KMV_ASSE,
                            MATR_MASS=M_ASSE,
							OPTION = 'BANDE',
                            CALC_FREQ = _F(FREQ = (-0.5, 0.0), ),
                            VERI_MODE=_F(STOP_ERREUR='NON'),                            
                            )
                            
LIMODE = RECU_TABLE(CO=MODES,NOM_PARA = 'FREQ')

pfreq = LIMODE.EXTR_TABLE().values()['FREQ']

print pfreq

for freq in pfreq:
    lam=(2*math.pi*freq)**2
    #print "calc_lam=%s"%lam
    print "calc_E=%s"%(c1*lam)
    
for n in range(1,5):
    radius=n**2*a0
    print "radius=%s"%radius
    wavelength=2*math.pi*n*a0
    #print "wavelength=%s"%wavelength
    expected_E=-(mu*e**4)/(32*math.pi**2*eps0**2*ht**2*n**2)
    print "expect_E=%s"%expected_E
                             
IMPR_RESU(FORMAT='MED',
          RESU=(
                _F(MAILLAGE=MAILRAF, RESULTAT=MODES,),
                #_F(MAILLAGE=MAILRAF,CHAM_GD=CHPOTV,),
                )
          )    

FIN()

There are many things going on here, but the most important ones are the material definition (used to express the potential field using a space-varying heat capacity) and a modal problem over a thermal model (this can be quite unusual but it is the corresponding problem with respect to the equations above).

You will also need a “hydrogen.export” file to keep everything together:


P actions make_etude
P actions make_etude
P consbtc oui
P corefilesize unlimited
P cpresok RESNOOK
P facmtps 1
P follow_output yes
P mem_aster 100.0
P memjob 1843200
P memory_limit 1800.0
P mode interactif
P mpi_nbcpu 1
P mpi_nbnoeud 1
P nbmaxnook 5
P ncpus 1
P soumbtc oui
P testlist verification sequential
P time_limit 60.0
P tpsjob 2
P version testing
A args 
A memjeveux 225.0
A tpmax 60.0
F comm hydrogen.comm D  1
F mmed Mesh3d10auModerate.med D  20
F resu hydrogen.resu R  8
F mess hydrogen.mess R  6
F rmed hydrogen.rmed R  80

In order to get better (more precise) results you would need a more refined mesh. Code_Aster can use an automatic refiner (by calling the MACR_ADAP_MAIL command) but it is disabled in the code above: we are still working on the Windows porting of that.

There are still few books about this techniques, here you have some references:

415D38oXGBL._SX326_BO1,204,203,200_
Finite Element and Boundary Element Applications in Quantum Mechanics by Ramdas Ram-Mohan
518At8nKaGL._SX365_BO1,204,203,200_
Molecular Quantum Mechanics by Peter Atkins and Ronald Friedman

 

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

w

Connecting to %s