scientificprotocols authored about 3 years ago

Authors: Chen Weijin, Zheng Yue & Wang Biao


In this protocol we present the typical procedures of conducting computational experiments on the domain structure of low dimensional ferroelectrics using a personal computer (PC). Based on an example of the phase-field model (PFM) for ferroelectric nanoplatelet, a computational experiment consisting of physical modeling, programming, simulation, result analysis and post processing, will be illustrated. Our physical model appropriately takes into account the effects of boundary conditions, inhomogeneous elastic and electrostatic fields, ambient temperature, size and shape of the low dimensional ferroelectrics, which makes it convenient to conduct a series of virtual experiments, e.g., applying external mechanical or electric field to affect the domain structure of nanostructures.


Low dimensional ferroelectrics, i.e., ferroelectrics with extreme geometrical confinement along some dimensions, are under active investigations for their important roles in device miniaturization and functional device development. (1,2) People are particularly concerned about the regular responses of domain structures in low dimensional ferroelectrics to the application of external fields. Nevertheless, to synthesize ferroelectric nanostructures of regular size and shape is still a difficult issue at present, not to say the experimental probe and control of nanoscale ferroelectric domain structures. A natural way to come out the predicament is to conduct computational experiments.

With appropriate modeling and simulating the physical processes, computational experiments can provide us convincing results to test the correctness of theory, interpret the results of hands-on experiments, and more importantly, to predict new regularity and phenomenon. Especially for situations where hands-on experiments are not feasible or quite expensive, conducting computational experiments is an attractive alternative. This is why recently computational physics has achieved great advance and becomes an important and individual subject of physics other than just an assistant tool of theoretical physics.

A complete computational experiment consists of several parts, including physical modeling, programming, simulation, result analysis and post processing. Physical processes are represented through numerically solving the governing equations, which depend on the theoretical description of the processes. Accordingly, different simulation methods could be adopted. For ferroelectrics, a quantity of simulation methods has been developed, including thermodynamic ones like phase-field model (PFM) (3,4), atomistic level simulations, e.g., molecular dynamic simulation (MD) (5) and effective Hamiltonian simulation (EH) (6), and first-principle calculations (FPC) (7), etc.

Basing on an example of phase-field model for ferroelectric nanoplatelet, in this protocol we present the procedures of conducting computational experiments on low dimensional ferroelectrics using a personal computer (PC). The detail of the phase-field model and derivation of finite element method to solve the electromechanical fields can be found in the associated publication.


  1. A personal computer (PC) with a typical operating system (Windows/Linux).
  2. C language compiler (the program is assumed written in C language), e.g., Visual C++ ( for Windows systems, or GNU Compiler Collection ( for Linux systems.
  3. MATLAB ( (version 7.1 or above).


I. Physical modeling of low dimensional ferroelectrics

  1. Determine your model system. You should determine the dimensionality of the system (i.e., 0-D nanoplatelet 1-D nanowire and 2-D nanofilm), material, and electric and mechanical boundary conditions, which are recommended to be written as optional tags in your program. In this demo, we consider a PbTiO3 ferroelectric nanoplatelet under surface traction and open-circuit condition as shown in Figure 1.
  2. Construct the system’s free energy. In general, the free energy of a ferroelectric is a functional of polarization, and should incorporate the effects of polarization inhomogeneity, elastic and electric fields. It can be written as a sum of the bulk free energy (including a Landau-type potential, gradient energy, elastic energy and electric energy) and surface free energy.
    • Tip: The construction of free energy is very important in phase-field model of ferroelectrics. The governing equation, i.e., the Time Dependent Ginzburg–Landau (TDGL) equation, describes that the evolution of domain structure toward its equilibrium is driven by the decrease of the system’s free energy.
    • Caution: The free energy should be generally constructed under the thermodynamic framework.
  3. Design your experiment. In our case, the experiment is to apply mechanical loads to the nanoplatelet to affect its domain structure.

II. Programming

The flow diagram of phase-field simulation at a given condition (e.g., fixed temperature, open-circuited and mechanical load) is as shown in Figure 2. Followings are the main functional parts of the corresponding program.

  1. Initialization. In this part of program, the parameters and variables needed for simulation would be defined and initialized. It is recommended to write this part in a way that it can either use the default values of the parameters and variables or read them from input files.
    • Tip: For numerical accuracy and convenience, it is recommended to make the parameters and variables dimensionless so that their values fall into a suitable range.
  2. Grid division of the system. This process would generate some tables, from which you can easily find out the information of the nodes and elements, such as its coordinates, adjacent nodes and elements. In our case, we make the following tables,
    • Tnode——Each row stores a node’s number, coordinates and adjacent elements,
    • Tnodes——Each row stores a node’s number and its adjacent nodes,
    • Telement——Each row stores an element’s number, coordinates and its nodes.
      • Tip: These tables can largely simplify the calculation of element stiffness matrices, element node displacement/potential vectors, and the assemblage of global stiffness matrices and displacement/potential vectors.
  3. Calculation of element stiffness matrices and the assemblage of global stiffness matrices. In most cases, the stiffness matrices (including elastic and electric) can be considered unchanged during the process, therefore the calculation only need to be done by once.
    • Tip: According to the boundary conditions, the computation size and accuracy requirement, the electric and elastic fields can be solved by different methods, such as conjugate gradient method (CGM), finite-element method (FEM), and fast Fourier transformation (FFT). For our case of a ferroelectric nanoplatelet with a moderate size, finite-element method is suitable to solve the electric and elastic fields.
  4. Calculation of the elastic field. Calculate the element node displacement vectors, and assemble them into the global node displacement vector. Solve the equation [Ku]{U}={Fu}. This should be done at each step of polarization evolution.
  5. Calculation of the electric field. Calculate the element node potential vectors, and assemble them into the global node potential vector. Solve the equation [Kφ]{Ф}={Fφ}. This should be done at each step of polarization evolution.
    • Tip: Iterative methods, e.g. the Gauss-Seidel iteration method, would be suitable to equations [Ku]{U}={Fu} and [Kφ]{Ф}={Fφ}.
  6. Calculation of the evolution force and polarization field at next step. Simple explicit difference methods or semi-implicit Fourier-spectral algorithms8 can be used to solve the governing equation, i.e., the TDGL equation.
  7. Error analysis. Calculate the error between the polarization field at this step and next step. If the error is small enough, end program; otherwise, repeat steps 3-6.
  8. Input and output. In our case, the input includes a parameter file storing the value of parameters and a configuration file storing the initial polarization field. The output includes a log file recording the monitoring information during the simulation and files storing the polarization field at selected steps.

III. Simulation

We would like to simulate how the mechanical load affects the formation of domain structure in the nanoplatelet.

  1. Prepare a set of parameter files. In our case, the parameter files differ only in the value of surface traction.
  2. For each parameter files, compile and run the program.
    • Tip: You can also make some modifications on the original program, so that the whole simulation can be fulfilled by running the program by once.

IV. Result analysis & Post processing

At this stage, we assume that there is an output file storing the polarization field, e.g., p_final.txt. The file stores the node’s number, its coordinates and the polarization components. Now we will show how to use MATLAB ( to visualize the domain structures.

  1. Start MATLAB.
  2. Import file p_final.txt into the Workspace as shown in Figure 3.
  3. Run the following commands to obtain a vector plot the domain structure as shown in Figure 4a,
    • quiver3(pfinal(:,3),pfinal(:,2),pfinal(:,4),pfinal(:,6),pfinal(:,5),pfinal(:,7)); axis equal; axis off;
  4. Run the following commands to add a color plot of the polarization magnitude to Figure 4a as shown in Figure 4b,
    • meshing=[10 3 30];
    • nx=meshing(1)+1;
    • ny=meshing(2)+1;
    • nz=meshing(3)+1;
    • m=0;
    • for i=1:nx
    • for j=1:ny
    • for k=1:nz
    • m=m+1;
    • P(i,j,k)= (pfinal(m,5)*pfinal(m,5)+pfinal(m,6)*pfinal(m,6)+pfinal(m,7)*pfinal(m,7))^0.5;
    • end
    • end
    • end
    • Xslice=[0,ny-1];Yslice=[0,nx-1];Zslice=[0,nz-1];
    • hold on;h =slice(x,y,z,P,Xslice,Yslice,Zslice);
    • set(h,’LineStyle’,’none’);
  5. Change the colormap to your preference. For example, run the following commands to as shown in Figure 4c,
    • colorround=[1 0 1];colorfloor=[0 1 0];
    • n_color=64;
    • for i=1:n_color
    • map(i,1)=colorround(1)+(i-1)*(colorfloor(1)-colorround(1))/(ncolor-1);
    • map(i,2)=colorround(2)+(i-1)*(colorfloor(2)-colorround(2))/(ncolor-1);
    • map(i,3)=colorround(3)+(i-1)*(colorfloor(3)-colorround(3))/(ncolor-1);
    • end
    • colormap(map);


In our case of a 10×3×30 meshing nanoplatelet at a given surface traction and open-circuited condition, the simulation time of polarization evolution from random toward equilibrium ranges from tens of minutes to several hours on a PC (3.00 GHz, 1.75GB internal memory). In general, the time of a simulation would strongly depend on the size of system, the algorithms, the required calculation accuracy, the computer performance and the simulated process. When phase transition happens, the simulation time would be much longer.


I. How can I make sure whether the results are right?

  1. Repeat some solved problems with your own program.
  2. Compare the results of FEM with those of other methods, such as finite difference method.
  3. Analyze whether the results are physical.

II. If the polarization field does converge to a suitable range:

  1. Check whether time step is too large.
  2. Check the other parameters.
  3. Check whether the elastic and electric field are calculated correctly.

Tip: Write the program based a former one, and make sure the former program is correct.

III. If the calculation of the elastic/electric field does not converge:

  1. Check the correctness of the tables that used to generate the stiffness matrices and force vectors.
  2. Check whether you have appropriately handled the boundary conditions.
  3. Try other algorithms.

Anticipated Results

Results of the experiment are presented in Figure 5. The results indicate that mechanical loads can affect the nucleation and balance of the vortices in the nanoplatelet, and lead to various equilibrium vortex domain structures.


  1. Scott, J. F. Nanoferroelectrics: statics and dynamics. J. Phys.: Condens. Matter 18, R361 (2006).
  2. Rørvik, P. M., Grande, T. & Einarsrud, M.-A. One-dimensional nanostructures of ferroelectric perovskites. Adv. Mater. 23, 4007-4034 (2011).
  3. Chen, L. Q. Phase-Field Method of Phase Transitions/Domain Structures in Ferroelectric Thin Films: A Review. J. Am. Ceram. Soc. 91, 1835–1844 (2008).
  4. Chen, W. J, Zheng, Y. & Wang, B. Phase field simulations of stress controlling the vortex domain structures in ferroelectric nanosheets. Appl. Phys. Lett. 100, 062901 (2012).
  5. Phillpot, S. R., Sinnott, S. B. & Asthagiri, Atomic-Level Simulation of Ferroelectricity in Oxides: Current Status and Opportunities. A. Annu. Rev. Mater. Res. 37, 239 (2007).
  6. Zhong, W., Vanderbilt, D & Rabe, K. M. First-principles theory of ferroelectric phase transitions for perovskites: The case of BaTiO3. Phys. Rev. B 52, 6301 (1995).
  7. Shimada, T., Umeno, Y. & Kitamura, T. Ab initio study of stress-induced domain switching in PbTiO3. Phys. Rev. B 77, 094105 (2008).
  8. Chen, L.Q., Shen J. Applications of semi-implicit Fourier-spectral method to phase field equations. Computer Physics Communications 108, 147 (1998).


Helpful discussions with Dr. D. C. Ma and Dr. S. P. Lin are gratefully acknowledged. The authors also acknowledge the financial support of the National Natural Science Foundation of China (NSFC) (Nos. 10902128, 11232015, 50802026, 10972239). Y. Z. also thanks support by the Fundamental Research Funds for the Central Universities to Micro&Nano Physics and Mechanics Research Laboratory, New Century Excellent Talents in University, Research Fund for the Doctoral Program of Higher Education, Fok Ying Tung Foundation, Natural Science Funds for Distinguished Young Scholar of Guangdong Province of China, and Educational Commission of Guangdong Province of China.


Figure 1 : Scheme of the model system.

Download Figure 1

A ferroelectric nanoplatelet under surface traction and open-circuited condition.

Figure 2: Flow diagram of phase-field simulation.

Download Figure 2

The polarization evolution is simulated at a given condition, e.g., T=300K, open-circuited and zero mechanical load. Finite element methods are used to solve the elastic and electric fields.

Figure 3: Scheme of importing file p_final.txt into MATLAB and illustration of the data.

Download Figure 3

Each row of the p_final matrix stores a node number, node coordinates and polarization components.

Figure 4: Some result of visualizing the domain structure of a nanoplatelet by using MATLAB.

Download Figure 4

(a) A vector plot the domain structure using quiver3() function. (b) Adding a color plot of the polarization magnitude using slice() function. (c) A plot after changing the colormap to your preference.

Figure 5: Anticipated results.

Download Figure 5

Snapshots of the domain morphologies during the polarization evolution for a nanoplatelet under different surface traction but with the same initial random perturbation at T=300K.

Associated Publications

Vortex Domain Structure in Ferroelectric Nanoplatelets and Control of its Transformation by Mechanical Load. W. J. Chen, Yue Zheng, and Biao Wang. Scientific Reports 2 () 12/11/2012 doi:10.1038/srep00796

Author information

Chen Weijin, Zheng Yue & Wang Biao, State Key Laboratory of Optoelectronic Materials and Technologies, Sun Yat-sen University, Guangzhou 510275, China. School of Physics and Engineering, Sun Yat-sen University, Guangzhou 510275, China.

Correspondence to: Zheng Yue ([email protected]) Wang Biao ([email protected])

Source: Protocol Exchange (2012) doi:10.1038/protex.2012.054. Originally published online 16 November 2012.

Average rating 0 ratings