Micro and Nano Mechanics Group - User contributions [en]

Perform Nanoindentation on Al-Mg Alloy

2022-10-12T04:47:44Z

Yfwang09: /* Retrieve the Forth-Depth Curve of Nanoindentation */

Nanoindentation on Al-Mg Alloy using LAMMPS 
<DIV>
Yifan Wang
</DIV>
 Created May 2017

This tutorial explains how to create an FCC sample of Al-Mg alloy and perform nanoindentation on the sample in LAMMPS.
<HR>

==Prerequisite==

For this simulation, we will use the LAMMPS 5 Nov 2016 version. The newer version may have changes of the commands. When compiling LAMMPS, please switch on the package MANYBODY in order to use the EAM potential. More details about compiling LAMMPS, please see the manual page (http://lammps.sandia.gov/doc/Section_start.html) for more information.

== Prepare The Sample of Al-Mg Alloy ==

In this note, we will use the following geometric parameters as an example, to show how to perform nanoindentation MD simulation in LAMMPS:

: Weight percentage: 0.6% (~0.67% atomic percentage)
: Simulation Box: 32.4x32.4x32.4nm (free surface on z-direction, pbc on x- and y-directions)
: Indenter radius: 16.2nm

=== Thermodynamics ===

According to the phase diagram of Al-Mg alloy (Fig.1), the solubility of Mg in Al is smaller than 1% at room temperature. For our case, since the atomic concentration is 0.67%, we can generate the sample by first creating pure FCC Al metal, and then randomly exchange required Al atoms to Mg.

[[Image:Al-Mg_phase_diagram.jpg | thumb | 500px | center | Fig.1 Phase Diagram of Aluminum-Magnesium (Murray, 1982)]]
[[Image:Al-Mg_size_factor.jpg | thumb | 500px | center | Fig.2 Size Factors of Al-Mg solid solution (Cai&Nix, 2016)]]

However, the size factor of Mg as the solute in Al is not negligible (Fig.2). We have to consider the volume change of the sample by adding Mg into pure Al. There are two ways to account for the size effect:

# Relax the sample under zero pressure after we generate the Al-Mg alloy.
# Estimate the volume change theoretically from the size factor.

Although the first method is more correct in the physical sense, the relaxation process will create surface steps on the free surface which will interfere the nanoindentation process. Here we estimate the volume change by correcting the lattice constant:

{|border="0" align="center"
|<math> a_{\text{alloy}} = a_{\text{Al}}(1+\chi_{\text{Mg}}l_{sf}) </math>.
|}

For now, we use the linear size factor from King's table. For a more accurate value, a constant pressure relaxation can be performed to determine the size factor.

=== Generate the Sample of Al-Mg Alloy ===

First we generate the simulation box of 32.4x32.4x32.4nm. We set up the simulation based on the lattice so that the entire simulation box is filled with atoms. First we define the lattice in LAMMPS ( indentation surface (001) ):

lattice fcc 4.053 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region BOX block -40.0 40.0 -40.0 40.0 -80.0 0.0 units lattice

Here the lattice constant 4.053 is the size after we consider the size effect, and 80 lattice distance is around 32.4nm. Note that if we change the indentation surface to other FCC planes ((110), (111), etc.), we should change the orientation in lattice definition, for example if we take (110) plane as the indentation surface, the setting should be:

lattice fcc 4.053 orient x 0.707 0 -0.707 orient y 0 1 0 orient z 0.707 0 0.707
region BOX block -40.0 40.0 -40.0 40.0 -80.0 0.0 units lattice

Then we create the simulation box and the Al atoms based on the region we defined:

create_box 2 BOX
create_atoms 1 region BOX
set region BOX type 1

And we randomly change 0.67% atoms to Mg atoms:

set region BOX type/fraction 2 0.0067 20170517 # random seed

Finally, after we add the Mg atoms, we perform relaxation to the entire sample to make the sample equilibrated. In this simulation, we use the EAM potential from Liu and Adams (1998) (https://www.ctcms.nist.gov/potentials/Download/Al-Mg-LA/almg.liu):

pair_style eam/alloy
pair_coeff * * almg.liu.eam.alloy Al Mg
minimize 1e-9 1e-9 10000 10000

We can visualize the sample as shown in Fig.3.

[[Image:Al-Mg_sample.jpg | thumb | 300px | center | Fig.3 Visualization of the Al-Mg alloy sample. Red: Al atoms; Blue: Mg atoms]]

== Perform MD Simulation of Nanoindentation on The Alloy Sample ==

Before the actual simulation, we first determine the dynamics of the indenter:

: Indenter stiffness constant K: 1.0eV/Å3
: Temperature T: 1.0K
: Indent rate: 0.81Å/ps = 81m/s
: Loading steps: 40000 (0.001ps/step)
: Indent depth: 3.24nm

Note: The units are defined by units command in LAMMPS, we use metal units (http://lammps.sandia.gov/doc/99/units.html)

=== The Force Field for the Indenter ===

In LAMMPS, a spherical indeter exerts a force of magnitude:

<math> F(r) = -K (r - R)^2 </math>

on each atom where K is the indenter stiffness constant, r is the distance from the atom to the center of the indenter, and R is the radius of the indenter. The force is repulsive and F(r)=0 for r>R.

We have tested the K value ranging from 1 to 1000. Note that with this force field function, the indenter is not an elastic sphere. The effective elastic modulus of the indenter is shown in Fig.4. We can see that the K value does not have significant influence on the effective modulus of the indenter. In our simulation, we choose K = 1.0eV/Å for all the simulation.

[[Image:Indenter_K_Modulus.png | thumb | 500px | center | Fig.4 The Effective Modulus of the Indenter in Terms of the Indent Depth.]]

=== Fix the bottom ===

First, We define 5 bottom layers as the bottom layer, these atom layers will be fixed to support the entire sample. Otherwise the sample will move downwards during the indentation.

region fixed_region block INF INF INF INF INF -75.0 units lattice
group bottom region fixed_region
group mobile subtract all bottom
fix bottom bottom setforce 0.0 0.0 0.0
velocity bottom set 0.0 0.0 0.0

Then we set up the temperature for the system:

compute new mobile temp
velocity mobile create 1.0 20110723 temp new # create temperature
fix ensemble mobile nvt temp 1.0 1.0 1.0

=== Loading and Unloading ===

Fig.5 is a diagram showing the loading process of the nanoindentation.

[[Image:example_nanoindentation.png | thumb | 150px | center | Fig.5 Loading of the nanoindentation process]]

In LAMMPS, the command used for indentation is "fix indent". This command can use the information of time steps and thus perform dynamic indentation during the MD simulation.

1. Calculate the indentation depth

variable indent_rate equal 0.81 # indent rate = 0.81Å/ps
variable depth equal step*dt*${indent_rate}

2. Calculate the z-direction position of the indenter center for loading steps:

variable indenter_R equal 162
variable z equal "v_indenter_R - v_depth"

Note that you cannot use ${depth} in this expression, otherwise only the depth value at the first time step will be substituted here and no time step information will be kept.

3. Fix the indenter. To use the standard units defined by units command, remember to add "units box" at the end of the fix command:

variable indent_K equal 1.0
fix indent all indent ${indent_K} sphere 0.0 0.0 v_z ${indenter_R} units box

Since the indenter is a "fix" command, the responsive force of the fixed constraint can be retrieved from variable "f_indent", which is a 3-component vector that represents the responsive force of this fixed constraint. In our case, we only need the force on z-direction, i.e. f_indent[3].

4. Run MD simulation to perform the loading steps:

variable loading_steps equal 40000
run ${loading_steps}

5. Similarly, when we perform the unloading steps, we recalculate the depth and z position of the indenter:

variable depth equal (2*${loading_steps}-step)*dt*${indent_rate}
variable z equal "v_indent_R - v_depth"
fix indent all indent ${indent_K} sphere 0.0 0.0 v_z ${indenter_R} units box
run ${loading_steps}

== Post processing ==

LAMMPS provides multiple ways to output the results of MD simulation. In this section, we discuss the two mostly used output: "dump" and "thermo".

=== Retrieve the dislocation structure from MD simulation ===

To obtain the dislocation structure, we need the detailed atom location information during the simulation. In LAMMPS, "dump" command is designed for the output of the atom profile.

Simply add the following lines before the "run" command, the atom profile will be dumped into dump file "dump.${timestep}.indent" for further analysis.

variable dumpfreq equal 1000
dump loading all atom ${dumpfreq} dump.*.indent
run 40000

If you want to set different dump frequency for unloading, you have to cancel the previous dump setting, and set up a different dump setting.

undump loading
dump unloading all atom 100 dump.*.unload
run 40000
undump unloading

For more settings about dump command, please read LAMMPS manual (http://lammps.sandia.gov/doc/dump.html).

After we get the atom profiles during the simulation, we can perform the dislocation analysis and get the dislocation structure, as shown in Fig.6:

[[Image:dislocation_analysis.png | thumb | 500px | center | Fig.6 dislocation nucleation during nanoindentation]]

=== Retrieve the Force-Depth Curve of Nanoindentation ===

As we discussed in section 3.3, the responsive force on the indenter is recorded in variable "f_indent[3]", the indentation depth is calculated and recorded in variable "v_depth". Therefore, we add the following lines before the "run" command in section 3.3, then we can retrieve the forth-depth curve from the simulation:

thermo 100 # output per 100 timestep
thermo_style custom step v_depth f_indent[3]

Then the information of timestep, depth and indentation force will be recorded in the log file. More details of "thermo" command can be found in LAMMPS official manual (http://lammps.sandia.gov/doc/thermo.html)

The force-depth curves for different indentation direction are shown below in Fig.7.

[[Image:FD_curve_orientation.png | thumb | 500px | center | Fig.7 Force-depth curve of nanoindation performed on different orientations]]

In addition, because we have periodic boundary condition on both x- and y- direction, there are actually infinite indenter images exerting infinite force on the sample. Therefore, the sample is not touching the indenter at very beginning, and the actual indentation process started at finite indentation depth. In the results, we manually corrected the indentation depth so that it starts from the position where the indenter is exactly touching the sample.

M11 MD++ Powered by Python

2019-10-12T03:37:19Z

Yfwang09:

MD++ with Python scripting language ([http://micro.stanford.edu/~caiwei/Forum/2004-12-12-MD++/PDF_Manuals/M11_MD++_v3.pdf Python])

commands

<pre>
mdpp.cmd(command_string)
'''Execute one command or multiple commands in MD++
Parameters
----------
command_string : string
a string that follows MD++ script command format
Returns
-------
exit_code : int
0 if success
'''
</pre>
<pre>
mdpp.get(varname)
'''Obtain the value of a variable in MD++ object
Parameters
----------
varname : string
variable name to obtain from the MD++ object
Returns
-------
value : int, long, double, or string
depends on the type of the variable
'''
</pre>
<pre>
mdpp.get_array(array_name, istart, iend, iskip=0)
'''Obtain an array from MD++ object to a python list
Parameters
----------
array_name : string
array name to obtain from the MD++ object
istart,iend,iskip : int
the starting position, end position, and number of elements to skip
Returns
-------
array : list
flattened list of requested array in MD++ object
'''
</pre>
<pre>
mdpp.set_array(input_list, array_name, istart, iskip=0)
'''Obtain an array from MD++ object to a python list
Parameters
----------
input_list : list
1D python list of flattened array values to be inserted.
array_name : string
array name to set in the MD++ object
istart,iskip : int
the starting position, and number of elements to skip
Returns
-------
exit_code : int
0 if success
'''
</pre>

example scripts

<UL>
<LI> [http://micro.stanford.edu/~caiwei/Forum/2004-12-12-MD++/Python/si.py si.py]
<LI> [http://micro.stanford.edu/~caiwei/Forum/2004-12-12-MD++/Python/si2.py si2.py]
<LI> [http://micro.stanford.edu/~caiwei/Forum/2004-12-12-MD++/Python/si-arg.py si-arg.py]
<LI> [http://micro.stanford.edu/~caiwei/Forum/2004-12-12-MD++/Python/si-temp.py si-temp.py]
</UL>

M11 MD++ Powered by Python

2019-10-12T03:36:38Z

[https://gitlab.com/micronano/ShElastic ShElastic] is the python code for the paper

Wang, Y., Zhang, X., & Cai, W. (2019). Spherical Harmonics Method for Computing the Image Stress Due to A Spherical Void. Journal of the Mechanics and Physics of Solids.

The code is tested under Ubuntu 16.04 LTS. with package python 3.6, numpy 1.15, scipy 1.0, matplotlib 3.0, jupyter, and SHTOOLS 4.3. Please see [http://micro.stanford.edu/wiki/Install_SHTOOLS Install_SHTOOLS] for more details on installing required packages

== Obtain the ShElastic package ==

Please download the package using git:

git clone https://gitlab.com/micronano/ShElastic.git

or directly from https://gitlab.com/micronano/ShElastic

== How to Install SHTOOLS on mc2 ==

On mc2 cluster, since we can only install packages for user only, we are not allowed to use sudo command to install the required packages. We need to compile the packages by our own.

For compiling shtools, the required packages are fftw3, lapack, and blas. Libraries for lapack and blas are available on mc2, we can use them directly for compiling shtools. For installing fftw3, please refer to [http://micro.stanford.edu/wiki/Install_FFTW3 Install_FFTW3] for more details.

After you followed the steps in FFTW3 installation (Serial only), you should have all the fftw3 libraries available:

~/usr/include/fftw3.h
~/usr/include/fftw3.f
~/usr/lib/libfftw3.a
~/usr/lib/libfftw3.la
~/usr/lib/libfftw3.so

And also, run the command:

ls /usr/lib64

The libraries for lapack and blas should be available here:

/usr/lib64/liblapack.a
/usr/lib64/liblapack.so
/usr/lib64/liblapack.so.3
...
/usr/lib64/libblas.a
/usr/lib64/libblas.so
/usr/lib64/libblas.so.3
...

Before we proceed, we need to make sure to have python installed.

=== Install python + scipy + SHTOOLS on local directory ===

We create a local folder .localpython for installing python:

mkdir ~/.localpython

download python source code:

mkdir ~/soft
cd ~/soft
wget http://www.python.org/ftp/python/X.X.X/Python-X.X.X.tgz
tar -xvf Python-X.X.X.tgz
cd Python-X.X.X

compile python from source and install to local directory:

make clean
./configure --enable-shared --prefix=${HOME}/.localpython
make
make install

Add the following lines to ~/.bashrc:

export PATH=${HOME}/.localpython/bin:$PATH
export PYTHONPATH=${HOME}/.localpython/bin
export LIBRARY_PATH=$LIBRARY_PATH:${HOME}/.localpython/lib
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${HOME}/.localpython/lib

if python version < 2.7.9, install pip by:

wget https://bootstrap.pypa.io/get-pip.py
python get-pip.py
export PATH=${HOME}/.local/bin:$PATH

Install SciPy Stack using pip on user's directory:

pip(pip3) install --user numpy scipy matplotlib ipython jupyter pandas sympy nose

Install pyshtools using pip:
First, we put the pre-compiled FFTW3 libraries into LIBRARY_PATH:

export LIBRARY_PATH=$LIBRARY_PATH:${HOME}/usr/lib
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${HOME}/usr/lib

Then we can install pyshtools through pip on user's directory:

pip(pip3) install --user pyshtools

=== Compile SHTOOLS ===

load the python module (we use python 2 as an example):

module load python/2.7.8 (Whichever is available)

Also make sure you have numpy available:

pip list

Then you are ready to install shtools on your mc2 account.

First, download the shtools source package from [https://github.com/SHTOOLS/SHTOOLS/releases Releases]

cd ~/soft
wget https://github.com/SHTOOLS/SHTOOLS/archive/v3.4.tar.gz
tar -zxvf v3.4.tar.gz
cd SHTOOLS-3.4

Then, indicate the libraries and compile the source code:

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:"$HOME/usr/lib"
make python2 LAPACK="-L/usr/lib64 -llapack" BLAS="-L/usr/lib64 -lblas" FFTW="-L$HOME/usr/lib -lfftw3"

The output should be something look like [https://shtools.oca.eu/shtools/www/makeout.html this]. Then your package is ready for use. If you are compiling with python 3, just change "python2" to "python3" in the command above.

To load the SHTOOLS module in python, use the following lines:

import sys
sys.path.append('$HOME/soft/SHTOOLS-3.4')
import pyshtools as shtools

where the "$HOME" has to be substituted by your own home folder path.

== Examples for using SHTOOLS in python ==

We consider a simple example where we do the spherical expansion of function <math>x^2</math> on the surface of a unit sphere.

First we import the necessary packages:

from __future__ import print_function # only necessary if using Python 2.x

import numpy as np

import sys
sys.path.append('/home/yfwang09/soft/SHTOOLS-3.4')
import pyshtools as psh

Then we generate the mesh grids for the spherical coordinates (Nx2N grids for latitude and longitude variables <math>\phi, \theta</math>; or <math>\pi/N</math> latitudinal and <math>2\pi/2N</math> longitudinal, both start from 0):

Ngrid = 100
phi = np.linspace(0, np.pi, Ngrid)
theta = np.linspace(0, 2*np.pi, 2*Ngrid)

THETA, PHI = np.meshgrid(theta, phi) # make the grid size as Nx2N

where the grid PHI and THETA both have N rows and 2N columns. Then we translate the Cartesian coordinates into the spherical meshgrids, and generate the target function into spherical coordinates:

Z = np.cos(PHI);
X = np.sin(PHI)*np.cos(THETA);
Y = np.sin(PHI)*np.sin(THETA);

X2 = X**2

Finally we expand the target function into spherical harmonic series:

cilm = np.empty([2,Ngrid/2,Ngrid/2])
filename = 'xsquare'
cilm = psh.expand.SHExpandDH(X2, sampling=2)
np.savetxt(filename+'_clm.txt', cilm[0,:,:], header='Row: l, Column: m')
np.savetxt(filename+'_slm.txt', cilm[1,:,:], header='Row: l, Column: m')

where cilm represent the coefficients of different spherical harmonic modes. cilm[0, l, m] and cilm[1, l, m] are the cosine and sine modes of spherical harmonics. The parameter sampling=2 means that we are using the grid Nx2N instead of NxN.

See the page [https://shtools.oca.eu/shtools/www/man/python/pyshexpanddh.html SHExpandDH] for the detailed description of the expanding function and the page [https://shtools.oca.eu/shtools/www/conventions.html Real Spherical Harmonics] for mathematical definition.

== Visualize spherical harmonic modes using scipy ==

We first import the necessary packages:

import numpy as np
from scipy.special import sph_harm
import matplotlib.pyplot as plt
from matplotlib import cm, colors
from mpl_toolkits.mplot3d import Axes3D

For a certain pair of parameters (m, l), we can plot the corresponding spherical harmonic mode Y(l, m) on the surface of a unit sphere using scipy and matplotlib packages. Here we use the meshgrid THETA, PHI, X, Y, and Z we generated in the previous section:

fvalues = sph_harm(m, l, THETA, PHI).real # Get the values of spherical harmonic functions
fmax, fmin = fvalues.max(), fvalues.min()
fcolors = (fvalues - fmin)/(fmax - fmin) # normalize the values into range [0, 1]

Then we plot Y(l, m) on a sphere:

fig = plt.figure(figsize=plt.figaspect(1.)) # make the axis with equal aspects
ax = fig.add_subplot(111, projection='3d') # add an subplot object with 3D plot
ax.plot_surface(x, y, z, rstride=1, cstride=1, facecolors=cm.Blues(fcolors))
plt.show()

For example, the figures for (l, m) = (2, 0), (3, 2), and (1, 1) are shown as below:
<gallery caption="Fig.1 Visualization of Spherical Harmonics" widths="200px" heights="200px" perrow="3">
Image:sph_harm_Y20.png|(a) <math> Y_{20}(\theta,\phi) </math>
Image:sph_harm_Y32.png|(b) <math> Y_{32}(\theta,\phi) </math>
Image:sph_harm_Y11.png|(c) <math> Y_{11}(\theta,\phi) </math>
</gallery>

Tutorials

2019-02-19T19:20:29Z

Yfwang09: /* Simulation Codes */

== Simulation Codes ==

{| class="wikitable" cellpadding="2"
|-
!width="200" | [[MD++ Manuals]]
!width="200" | [[ParaDiS Manuals]]
!width="300" | [[ParaDiS Workshop Notes]]
|-
!width="200" | [[DDLab Manuals]]
!width="200" | [[VASP Manuals]]
!width="300" | [[How to compile Qbox]]
|-
!width="200" | [[NetworCh Manuals]]
!width="200" | [[How to compile LAMMPS]]
!width="300" | [[How to compile pimc++]]
|-
!width="200" | [[ParaDiS tools]]
!width="200" | [[ShElastic Manuals]]
!width="300" |
|}

== Computers ==

{| class="wikitable" cellpadding="2"
|-
!width="200" | [[Unix Basics]]
!width="200" | [[Computer Setup]]
!width="200" | [[Computing Clusters | Parallel Clusters]]
|-
!width="200" | [[Install Ubuntu]]
!width="200" | [[Install FFTW3]]
!width="200" | [[Install HDF5]]
|-
!width="200" | [[Install GSL]]
!width="200" | [[Standard Tcl Library | Tcl Library]]
!width="200" | [[Install the GNU Compiler Collection (GCC)]]
|-
!width="200" | [[Running Hybrid MPI/OpenMP Simulations]]
!width="200" | [[Install SHTOOLS]]
|}

== Scientific Background ==

{| class="wikitable" cellpadding="2"
|-
!width="200" | [[Dislocations]]
|}

== Science Outreach ==

{| class="wikitable" cellpadding="2"
|-
!width="200" | [http://www.youtube.com/watch?v=SgM-Xes16Sw Outreach interview]
!width="200" | [[How to install MD++ in Ubuntu | Install MD++ in Ubuntu]]
!width="200" | [[Introduction to Molecular Dynamics Simulations of Fullerenes | MD of Fullerenes]]
|-
!width="200" | [[Images of Fullerenes]]
!width="200" | [[Android Applications and Information]]
|

|}

 

== [[Tutorial:Members_Only | Members Only]] ==

Install SHTOOLS

2019-02-19T19:11:31Z

Yfwang09: /* Visualization of spherical harmonic modes using scipy */

[https://shtools.oca.eu/shtools/ SHTOOLS] is a package for spherical harmonics transformation and reconstructions. The current version (SHTOOLS3.4) is available in python and Fortran 95. In this article, we are introducing the installation of SHTOOLS python packages on your own computer and mc2 cluster.

== How to Install SHTOOLS on Your Own Computer ==

On Mac or Linux, you will need to have fortran compiler, lapack/blas, and fftw3 available. This can be done simply by running the following commands with homebrew on Mac:

brew install fftw --with-fortran

or Linux:

sudo apt-get install libblas-dev liblapack-dev g++ gfortran libfftw3-dev tcsh

If you have installed python and the package manager pip, you can install the shtools directly by running the command:

pip install pyshtools

Then if you run command:

pip list

You will see pyshtools library is available in your python packages.

== How to Install SHTOOLS on mc2 ==

On mc2 cluster, since we can only install packages for user only, we are not allowed to use sudo command to install the required packages. We need to compile the packages by our own.

For compiling shtools, the required packages are fftw3, lapack, and blas. Libraries for lapack and blas are available on mc2, we can use them directly for compiling shtools. For installing fftw3, please refer to [http://micro.stanford.edu/wiki/Install_FFTW3 Install_FFTW3] for more details.

After you followed the steps in FFTW3 installation (Serial only), you should have all the fftw3 libraries available:

~/usr/include/fftw3.h
~/usr/include/fftw3.f
~/usr/lib/libfftw3.a
~/usr/lib/libfftw3.la
~/usr/lib/libfftw3.so

And also, run the command:

ls /usr/lib64

The libraries for lapack and blas should be available here:

/usr/lib64/liblapack.a
/usr/lib64/liblapack.so
/usr/lib64/liblapack.so.3
...
/usr/lib64/libblas.a
/usr/lib64/libblas.so
/usr/lib64/libblas.so.3
...

Before we proceed, we need to make sure to have python installed.

=== Install python + scipy + SHTOOLS on local directory ===

We create a local folder .localpython for installing python:

mkdir ~/.localpython

download python source code:

mkdir ~/soft
cd ~/soft
wget http://www.python.org/ftp/python/X.X.X/Python-X.X.X.tgz
tar -xvf Python-X.X.X.tgz
cd Python-X.X.X

compile python from source and install to local directory:

make clean
./configure --enable-shared --prefix=${HOME}/.localpython
make
make install

Add the following lines to ~/.bashrc:

export PATH=${HOME}/.localpython/bin:$PATH
export PYTHONPATH=${HOME}/.localpython/bin
export LIBRARY_PATH=$LIBRARY_PATH:${HOME}/.localpython/lib
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${HOME}/.localpython/lib

if python version < 2.7.9, install pip by:

wget https://bootstrap.pypa.io/get-pip.py
python get-pip.py
export PATH=${HOME}/.local/bin:$PATH

Install SciPy Stack using pip on user's directory:

pip(pip3) install --user numpy scipy matplotlib ipython jupyter pandas sympy nose

Install pyshtools using pip:
First, we put the pre-compiled FFTW3 libraries into LIBRARY_PATH:

export LIBRARY_PATH=$LIBRARY_PATH:${HOME}/usr/lib
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${HOME}/usr/lib

Then we can install pyshtools through pip on user's directory:

pip(pip3) install --user pyshtools

=== Compile SHTOOLS ===

load the python module (we use python 2 as an example):

module load python/2.7.8 (Whichever is available)

Also make sure you have numpy available:

pip list

Then you are ready to install shtools on your mc2 account.

First, download the shtools source package from [https://github.com/SHTOOLS/SHTOOLS/releases Releases]

cd ~/soft
wget https://github.com/SHTOOLS/SHTOOLS/archive/v3.4.tar.gz
tar -zxvf v3.4.tar.gz
cd SHTOOLS-3.4

Then, indicate the libraries and compile the source code:

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:"$HOME/usr/lib"
make python2 LAPACK="-L/usr/lib64 -llapack" BLAS="-L/usr/lib64 -lblas" FFTW="-L$HOME/usr/lib -lfftw3"

The output should be something look like [https://shtools.oca.eu/shtools/www/makeout.html this]. Then your package is ready for use. If you are compiling with python 3, just change "python2" to "python3" in the command above.

To load the SHTOOLS module in python, use the following lines:

import sys
sys.path.append('$HOME/soft/SHTOOLS-3.4')
import pyshtools as shtools

where the "$HOME" has to be substituted by your own home folder path.

== Examples for using SHTOOLS in python ==

We consider a simple example where we do the spherical expansion of function <math>x^2</math> on the surface of a unit sphere.

First we import the necessary packages:

from __future__ import print_function # only necessary if using Python 2.x

import numpy as np

import sys
sys.path.append('/home/yfwang09/soft/SHTOOLS-3.4')
import pyshtools as psh

Then we generate the mesh grids for the spherical coordinates (Nx2N grids for latitude and longitude variables <math>\phi, \theta</math>; or <math>\pi/N</math> latitudinal and <math>2\pi/2N</math> longitudinal, both start from 0):

Ngrid = 100
phi = np.linspace(0, np.pi, Ngrid)
theta = np.linspace(0, 2*np.pi, 2*Ngrid)

THETA, PHI = np.meshgrid(theta, phi) # make the grid size as Nx2N

where the grid PHI and THETA both have N rows and 2N columns. Then we translate the Cartesian coordinates into the spherical meshgrids, and generate the target function into spherical coordinates:

Z = np.cos(PHI);
X = np.sin(PHI)*np.cos(THETA);
Y = np.sin(PHI)*np.sin(THETA);

X2 = X**2

Finally we expand the target function into spherical harmonic series:

cilm = np.empty([2,Ngrid/2,Ngrid/2])
filename = 'xsquare'
cilm = psh.expand.SHExpandDH(X2, sampling=2)
np.savetxt(filename+'_clm.txt', cilm[0,:,:], header='Row: l, Column: m')
np.savetxt(filename+'_slm.txt', cilm[1,:,:], header='Row: l, Column: m')

where cilm represent the coefficients of different spherical harmonic modes. cilm[0, l, m] and cilm[1, l, m] are the cosine and sine modes of spherical harmonics. The parameter sampling=2 means that we are using the grid Nx2N instead of NxN.

See the page [https://shtools.oca.eu/shtools/www/man/python/pyshexpanddh.html SHExpandDH] for the detailed description of the expanding function and the page [https://shtools.oca.eu/shtools/www/conventions.html Real Spherical Harmonics] for mathematical definition.

== Visualize spherical harmonic modes using scipy ==

We first import the necessary packages:

import numpy as np
from scipy.special import sph_harm
import matplotlib.pyplot as plt
from matplotlib import cm, colors
from mpl_toolkits.mplot3d import Axes3D

For a certain pair of parameters (m, l), we can plot the corresponding spherical harmonic mode Y(l, m) on the surface of a unit sphere using scipy and matplotlib packages. Here we use the meshgrid THETA, PHI, X, Y, and Z we generated in the previous section:

fvalues = sph_harm(m, l, THETA, PHI).real # Get the values of spherical harmonic functions
fmax, fmin = fvalues.max(), fvalues.min()
fcolors = (fvalues - fmin)/(fmax - fmin) # normalize the values into range [0, 1]

Then we plot Y(l, m) on a sphere:

fig = plt.figure(figsize=plt.figaspect(1.)) # make the axis with equal aspects
ax = fig.add_subplot(111, projection='3d') # add an subplot object with 3D plot
ax.plot_surface(x, y, z, rstride=1, cstride=1, facecolors=cm.Blues(fcolors))
plt.show()

For example, the figures for (l, m) = (2, 0), (3, 2), and (1, 1) are shown as below:
<gallery caption="Fig.1 Visualization of Spherical Harmonics" widths="200px" heights="200px" perrow="3">
Image:sph_harm_Y20.png|(a) <math> Y_{20}(\theta,\phi) </math>
Image:sph_harm_Y32.png|(b) <math> Y_{32}(\theta,\phi) </math>
Image:sph_harm_Y11.png|(c) <math> Y_{11}(\theta,\phi) </math>
</gallery>

Install SHTOOLS

2018-01-13T01:27:25Z

Yfwang09: /* Install python + scipy + SHTOOLS on local directory */

[https://shtools.oca.eu/shtools/ SHTOOLS] is a package for spherical harmonics transformation and reconstructions. The current version (SHTOOLS3.4) is available in python and Fortran 95. In this article, we are introducing the installation of SHTOOLS python packages on your own computer and mc2 cluster.

== How to Install SHTOOLS on Your Own Computer ==

On Mac or Linux, you will need to have fortran compiler, lapack/blas, and fftw3 available. This can be done simply by running the following commands with homebrew on Mac:

brew install fftw --with-fortran

or Linux:

sudo apt-get install libblas-dev liblapack-dev g++ gfortran libfftw3-dev tcsh

If you have installed python and the package manager pip, you can install the shtools directly by running the command:

pip install pyshtools

Then if you run command:

pip list

You will see pyshtools library is available in your python packages.

== How to Install SHTOOLS on mc2 ==

On mc2 cluster, since we can only install packages for user only, we are not allowed to use sudo command to install the required packages. We need to compile the packages by our own.

For compiling shtools, the required packages are fftw3, lapack, and blas. Libraries for lapack and blas are available on mc2, we can use them directly for compiling shtools. For installing fftw3, please refer to [http://micro.stanford.edu/wiki/Install_FFTW3 Install_FFTW3] for more details.

After you followed the steps in FFTW3 installation (Serial only), you should have all the fftw3 libraries available:

~/usr/include/fftw3.h
~/usr/include/fftw3.f
~/usr/lib/libfftw3.a
~/usr/lib/libfftw3.la
~/usr/lib/libfftw3.so

And also, run the command:

ls /usr/lib64

The libraries for lapack and blas should be available here:

/usr/lib64/liblapack.a
/usr/lib64/liblapack.so
/usr/lib64/liblapack.so.3
...
/usr/lib64/libblas.a
/usr/lib64/libblas.so
/usr/lib64/libblas.so.3
...

Before we proceed, we need to make sure to have python installed.

=== Install python + scipy + SHTOOLS on local directory ===

We create a local folder .localpython for installing python:

mkdir ~/.localpython

download python source code:

mkdir ~/soft
cd ~/soft
wget http://www.python.org/ftp/python/X.X.X/Python-X.X.X.tgz
tar -xvf Python-X.X.X.tgz
cd Python-X.X.X

compile python from source and install to local directory:

make clean
./configure --enable-shared --prefix=${HOME}/.localpython
make
make install

Add the following lines to ~/.bashrc:

export PATH=${HOME}/.localpython/bin:$PATH
export PYTHONPATH=${HOME}/.localpython/bin
export LIBRARY_PATH=$LIBRARY_PATH:${HOME}/.localpython/lib
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${HOME}/.localpython/lib

if python version < 2.7.9, install pip by:

wget https://bootstrap.pypa.io/get-pip.py
python get-pip.py
export PATH=${HOME}/.local/bin:$PATH

Install SciPy Stack using pip on user's directory:

pip(pip3) install --user numpy scipy matplotlib ipython jupyter pandas sympy nose

Install pyshtools using pip:
First, we put the pre-compiled FFTW3 libraries into LIBRARY_PATH:

export LIBRARY_PATH=$LIBRARY_PATH:${HOME}/usr/lib
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${HOME}/usr/lib

Then we can install pyshtools through pip on user's directory:

pip(pip3) install --user pyshtools

=== Compile SHTOOLS ===

load the python module (we use python 2 as an example):

module load python/2.7.8 (Whichever is available)

Also make sure you have numpy available:

pip list

Then you are ready to install shtools on your mc2 account.

First, download the shtools source package from [https://github.com/SHTOOLS/SHTOOLS/releases Releases]

cd ~/soft
wget https://github.com/SHTOOLS/SHTOOLS/archive/v3.4.tar.gz
tar -zxvf v3.4.tar.gz
cd SHTOOLS-3.4

Then, indicate the libraries and compile the source code:

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:"$HOME/usr/lib"
make python2 LAPACK="-L/usr/lib64 -llapack" BLAS="-L/usr/lib64 -lblas" FFTW="-L$HOME/usr/lib -lfftw3"

The output should be something look like [https://shtools.oca.eu/shtools/www/makeout.html this]. Then your package is ready for use. If you are compiling with python 3, just change "python2" to "python3" in the command above.

To load the SHTOOLS module in python, use the following lines:

import sys
sys.path.append('$HOME/soft/SHTOOLS-3.4')
import pyshtools as shtools

where the "$HOME" has to be substituted by your own home folder path.

== Examples for using SHTOOLS in python ==

We consider a simple example where we do the spherical expansion of function <math>x^2</math> on the surface of a unit sphere.

First we import the necessary packages:

from __future__ import print_function # only necessary if using Python 2.x

import numpy as np

import sys
sys.path.append('/home/yfwang09/soft/SHTOOLS-3.4')
import pyshtools as psh

Then we generate the mesh grids for the spherical coordinates (Nx2N grids for latitude and longitude variables <math>\phi, \theta</math>; or <math>\pi/N</math> latitudinal and <math>2\pi/2N</math> longitudinal, both start from 0):

Ngrid = 100
phi = np.linspace(0, np.pi, Ngrid)
theta = np.linspace(0, 2*np.pi, 2*Ngrid)

THETA, PHI = np.meshgrid(theta, phi) # make the grid size as Nx2N

where the grid PHI and THETA both have N rows and 2N columns. Then we translate the Cartesian coordinates into the spherical meshgrids, and generate the target function into spherical coordinates:

Z = np.cos(PHI);
X = np.sin(PHI)*np.cos(THETA);
Y = np.sin(PHI)*np.sin(THETA);

X2 = X**2

Finally we expand the target function into spherical harmonic series:

cilm = np.empty([2,Ngrid/2,Ngrid/2])
filename = 'xsquare'
cilm = psh.expand.SHExpandDH(X2, sampling=2)
np.savetxt(filename+'_clm.txt', cilm[0,:,:], header='Row: l, Column: m')
np.savetxt(filename+'_slm.txt', cilm[1,:,:], header='Row: l, Column: m')

where cilm represent the coefficients of different spherical harmonic modes. cilm[0, l, m] and cilm[1, l, m] are the cosine and sine modes of spherical harmonics. The parameter sampling=2 means that we are using the grid Nx2N instead of NxN.

See the page [https://shtools.oca.eu/shtools/www/man/python/pyshexpanddh.html SHExpandDH] for the detailed description of the expanding function and the page [https://shtools.oca.eu/shtools/www/conventions.html Real Spherical Harmonics] for mathematical definition.

== Visualization of spherical harmonic modes using scipy ==

We first import the necessary packages:

import numpy as np
from scipy.special import sph_harm
import matplotlib.pyplot as plt
from matplotlib import cm, colors
from mpl_toolkits.mplot3d import Axes3D

For a certain pair of parameters (m, l), we can plot the corresponding spherical harmonic mode Y(l, m) on the surface of a unit sphere using scipy and matplotlib packages. Here we use the meshgrid THETA, PHI, X, Y, and Z we generated in the previous section:

fvalues = sph_harm(m, l, THETA, PHI).real # Get the values of spherical harmonic functions
fmax, fmin = fvalues.max(), fvalues.min()
fcolors = (fvalues - fmin)/(fmax - fmin) # normalize the values into range [0, 1]

Then we plot Y(l, m) on a sphere:

fig = plt.figure(figsize=plt.figaspect(1.)) # make the axis with equal aspects
ax = fig.add_subplot(111, projection='3d') # add an subplot object with 3D plot
ax.plot_surface(x, y, z, rstride=1, cstride=1, facecolors=cm.Blues(fcolors))
plt.show()

For example, the figures for (l, m) = (2, 0), (3, 2), and (1, 1) are shown as below:
<gallery caption="Fig.1 Visualization of Spherical Harmonics" widths="200px" heights="200px" perrow="3">
Image:sph_harm_Y20.png|(a) <math> Y_{20}(\theta,\phi) </math>
Image:sph_harm_Y32.png|(b) <math> Y_{32}(\theta,\phi) </math>
Image:sph_harm_Y11.png|(c) <math> Y_{11}(\theta,\phi) </math>
</gallery>