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):
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:
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:
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:
Code_Aster for Windows 