Dynamic atomic force microscopy simulations over a viscoelastic material

Content under Creative Commons Attribution license CC-BY 4.0 version,

Enrique A. López-Guerra.

This notebook contains atomic force microscopy (AFM) dynamic simulations for the case of a tapping mode simulation. In this example only the 1st eigenmode is excited, but you can easily modify it to excite up to three eigenmodes. The cantilever dynamics are assumed to be contained in the first three eigenmodes.

The simulation corresponds to the case of an AFM spherical tip interacting with a viscoelastic surface. The viscoelastic model is a generalized Maxwell model (Wiechert model) containing a large number of characteristic times. The contact mechanics have been implemented with respect to the classical theory of Lee and Radok: Lee, E. Ho, and Jens Rainer Maria Radok. "The contact problem for viscoelastic bodies." Journal of Applied Mechanics 27.3 (1960): 438-444.

Let's first start by importing some useful libraries


In [1]:
import sys
sys.path.append('d:\Github\pycroscopy')
from __future__ import division, absolute_import, print_function
from pycroscopy.simulation.afm_lib import dynamic_spectroscopy
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline

Inserting simulation parameters


In [2]:
fo1, fo2, fo3 = 45.0e3, 280.0e3, 17.6*45.0e3 #eigenmodes' resonance frequencies
k_m1, k_m2, k_m3 = 5.80, 210.0, 5.80*(fo3/fo1) #1st eigenmode cantilever stiffness
A1, A2, A3 = 50.0e-9, 0.0, 0.0  #target oscillating free amplitude of the 1st three eigenmodes
Q1, Q2, Q3 = 167.0,340.0, 500.0 #quality factor of the 1st three eigenmodes
          
R = 10.0e-9 #tip radius

period1, period2 = 1.0/fo1, 1.0/fo2 #oscillating period of first two eigenmodes

dt = period1/10.0e3 #simulation timestep
startprint = 3.0*Q1*period1 #starting point when results will start to get printed
simultime = startprint + 25.0*period1 #total simulation time
printstep = dt*10.0 #how often the results will be stored

z_step = 0.05*A1 #the spatial step between tapping mode runs (i.e., spatial vertical distance the cantilever base moves between tapping mode runs, the smaller this number is the more points your amplitude phase curve will have)

OPTIONAL: Pulling sample parameters corresponding to Polyisobutylene (data obtained from: Brinson, Hal F., and L. Catherine Brinson. "Polymer engineering science and viscoelasticity." An Introduction (2008).)

If you want to use this model, you should place this and the following cell right before the cell labeled as Main portion of the simulation.


In [3]:
#Sample parameters for polyisobutylene
df_G = pd.read_csv('PIB.txt', delimiter='\t', header=None)
tau = df_G.iloc[:,0].values  #relaxation times of the generalized Maxwell model
G = df_G.iloc[:,1].values #moduli of the springs in the Maxwell arms 
Ge = 0.0 #equilibrium modulus (rubbery modulus)
H = 5.0e-19 #Hammaker constant

Defining sample parameters: a simple model with two relaxation times


In [4]:
G = np.array([1.0e9,1.0e7]) #moduli of the springs in the Maxwell arms
tau = np.array([1.0/fo1/10.0,1.0/fo1]) #relaxation times of the generalized Maxwell model
Ge = 2.0e8 #equilibrium modulus (rubbery modulus)
H = 5.0e-19 #Hammaker constant

Main portion of the simulation

The simulations described correspond to dynamic AFM spectroscopy simulations. Here, the cantilever will be brought towards the sample approaching with discrete steps and will be allowed to oscillate until achieving a quasi-steady state. From the tip trajectories and other information recoreded (e.g., tip-sample force, sample position) in the steady state it will be possible to retrieve common information recorded in a dynamic spectroscopy experiment (e.g., amplitude, phase curves).


In [5]:
print('Note that this cell takes a while to execute')
import time
start_time = time.time()
amp, phase, zeq, Ediss, virial, peakF, maxdepth, t_a, tip_a, Fts_a, xb_a = dynamic_spectroscopy(G, tau, R, dt, startprint, simultime, fo1, fo2, fo3, k_m1, k_m2,k_m3, A1, A2, A3, printstep, Ge, Q1, Q2, Q3, H, z_step)
end_time = time.time()
print("total time taken this loop: %5.2f seconds"%(end_time - start_time))


Note that this cell takes a while to execute
total time taken this loop: 100.00 seconds

Plotting dynamic spectroscopy curves (Amplitude, Phase approach curves)


In [6]:
plt.plot(amp/A1, phase,'-o', color ='b')
plt.xlabel('$A_1/A_0$', fontsize =15)
plt.ylabel('Phase, deg', fontsize=15)


Out[6]:
Text(0,0.5,'Phase, deg')

In [7]:
plt.plot(zeq*1.0e9, amp*1.0e9,'-o', color ='g')
plt.xlabel('$Z_{eq}, nm$', fontsize =15)
plt.ylabel('Amplitude, nm', fontsize=15)


Out[7]:
Text(0,0.5,'Amplitude, nm')

Plotting a typical force-distance curve


In [8]:
N = 10
setpoint = amp[N]/A1
plt.plot(tip_a[N]*1.0e9, Fts_a[N]*1.0e9, color ='g', label = '$setpoint: %.2f$'%setpoint)
plt.legend(loc='best')
plt.xlabel('$tip, nm$', fontsize =15)
plt.xlim(-5,10)
plt.ylabel('$F_{ts}, nN$', fontsize=15)


Out[8]:
Text(0,0.5,'$F_{ts}, nN$')

Plotting the tip trajectory


In [9]:
N = 10
setpoint = amp[N]/A1
plt.plot(t_a*1.0e6, tip_a[N]*1.0e9, color ='g', label = '$setpoint: %.2f$'%setpoint)
plt.legend(loc='best')
plt.xlim(0,t_a[int(len(t_a)/10.0)]*1.0e6)
plt.ylim(min(tip_a[N])*5.0e9, max(tip_a[N])*1.2e9)
plt.xlabel('$tip, nm$', fontsize =15)
plt.ylabel('$time, \mu s$', fontsize=15)


Out[9]:
Text(0,0.5,'$time, \\mu s$')

Plotting more results


In [10]:
plt.figure(1)
plt.plot(amp/A1, -maxdepth*1.0e9,'-o', color ='c')
plt.xlabel('$A_1/A_0$', fontsize =15)
plt.ylabel('Maximum Penetration, nm', fontsize=15)
#plt.savefig('Max_penetration.png', bbox_inches='tight') #optional to save the figure in current path

plt.figure(2)
plt.plot(amp/A1, peakF*1.0e9,'-o', color ='m')
plt.xlabel('$A_1/A_0$', fontsize =15)
plt.ylabel('Peak Force, nN', fontsize=15)
#plt.savefig('Max_Force.png', bbox_inches='tight')  #optional line to save the figure in current path


Out[10]:
Text(0,0.5,'Peak Force, nN')

In [ ]: