feram: MD simulator for bulk and thin-film ferroelectrics

feram logo

Table of Contents:

Homepage, download, compilation and publication

feram's homepage is http://loto.sourceforge.net/feram/ .

You can freely download a tarball package of feram (feram-X.YY.ZZ.tar.xz) from http://sourceforge.net/projects/loto/files/feram/ . feram is GPLed free software.

Step-by-step procedure of compilation is described in INSTALL.

Publications related to feram is listed in doc/publication.html.

What is feram?

feram is a fast molecular dynamics (MD) simulator for bulk and thin-film ferroelectrics and relaxors.

Features and keywords

In Physics




feram reads parameter files specified as command-line arguments.

$ OMP_NUM_THREADS=6 SOMEWHERE/feram foo.feram bar.feram baz.feram

Optionally, you can make MPI version of feram, feram_mpi. feram_mpi processes input files in parallel.

$ OMP_NUM_THREADS=16 mpiexec -np 3 SOMEWHERE/feram_mpi foo.feram bar.feram baz.feram

The number of threads for each process should be given with environment variable OMP_NUM_THREADS. The most efficient value for OMP_NUM_THREADS is depends on the system size (N = L_x L_y L_z), the method of calculation and your computer. Latest Intel Xeon machines with GNU/Linux OS give high performance to feram. feram does not require large memory, but requires broad memory band width for high performance computing. Up to system size of 64x64x64, execution on 1 slot of CPU may gives good performance. For lager system size calculations, multi-slot NUMA execution may gives better performance. Use taskset(1) or numactl(8) to retrieve or set a process's CPU affinity.

Input files

parameter input file (.feram file)

The parameter file for feram is a text file consisting of comment lines and 'tag = value(s)' lines. Filenames of parameter files may be foo123 or foo123.feram. In these cases, filenames of output files will be foo123.*,e.g. foo123.log, foo123.0000050000.cood, etc.

How to determine the parameters is described in parameters/parameters.html.


Lines beginning with '#' are ignored. Blank lines are also ignored.

# This is a comment line.

# Here are two more
# comment lines.
shebang: "#!/usr/local/bin/feram" or "#!./feram"

You can write a shebang at the first line of a .feram input file and execute the file as a script.


You must put ' = ', space-equal-space, between tag and value(s) as:

tag = 1.0

tag = -2.0 -3.0 -4.0
tag =  5.0  6.0  7.0

Followings are currently available tags.


verbose tag determines how much messages will be written to the .log file. (Details of .log file is described below.)

verbose = 0

(NOT IMPLEMENTED YET) .log file will not be created.

verbose = 1

1 is the default value for verbose tag. Energies of the first and last iterations will be reported in the .log file.

verbose = 2

When 2 is given, energies will be reported in the .log file for EVERY iterations.

verbose = 3

When 3 is given, detailed messages will be written into the .log file.

verbose = 4

When 4 is given, more detailed messages will be written into the .log file.

method = 'md'

'md' is for a molecular-dynamics simulation in the canonical ensemble using the Nosé-Poincaré thermostat.

method = 'vs'

'vs' is for a molecular-dynamics simulation in the canonical ensemble using the velocity scaling thermostat.

method = 'lf'

'lf' is for a molecular-dynamics simulation in the microcanonical ensemble using the leap-frog method.

method = 'mc'

'mc' is for a Monte Carlo simulation, but it is NOT IMPLEMENTED YET.

method = 'hl'

This 'hl' feature will be obsolete. For simulations of hysteresis loops, use n_E_wave_period and E_wave_type instead.


Pressure in GPa unit.

GPa = -5.0

You can apply uniaxial normal stress by giving three values to GPa.

GPa = 0.0 0.0 -3.0

Temperature in Kelvin.

kelvin = 100

Effective mass of u(R), the length of dipole on the unit cell R. In atomic mass unit.


If acoustic_mass_amu is not set, acoustic displacements, i.e. inhomogeneous strains, are optimized according to the {u(R)}. If it is set to the effective mass of acoustic displacements

acoustic_mass_amu = 46.44

molecular dynamics is performed on acoustic displacements. Its unit is atomic mass unit. Generally, it is average of mass of 5 atoms. You cannot use method = 'md' together with acoustic_mass_amu, currently.


Q_Nose for Nose-Poincare thermostat. For large system and high temperature, use large value. For Small system and low temperature, use small value.

Q_Nose = 14.4

Specifies the structure of the system. 'bulk' for infinitely periodic bulk 'film' for free standing thin film. 'epit' for epitaxially strained thin film. See epi_strain tag.

bulk_or_film = 'epit'

System size, L_x, L_y and L_z. They must be equal or larger than 4.

L = 32 32 4

Padding for speeding up the calculations. You do not have to give this value, if you cannot understand it. Best padding_y is system-dependent, but for 64x64x64 and 96x96x96

padding_y = 5

for 128x128x128

padding_y = 1

may be good. Its default value is 3.


0: There is no dead layer. 1: Single dead layer is on the bottom electrode. 2: There are dead layers both on the bottom and top electrodes.

gap_id = 1

Thickness of dead layer is always 1 unit cell.


Polarization of dead layer(s) in Angstrom unit.

gap_dipole_u = 0.0 0.0 0.16

Lattice constant for the perovskite ferroelectrics in Angstrom unit.

a0 = 3.99

Epitaxial strain.

epi_strain = -0.01

Note that epi_strain is effective only when bulk_or_film = 'epit'. See src/optimize-homo-strain.F for more details.


Time step in pico second.

dt = 0.002
n_thermalize, n_average, n_coord_freq, n_hl_freq

The number of time steps of thermalizing, averaging, frequency of taking snapshots, and frequency of writing to .hl file.

n_thermalize = 40000
n_average    = 10000
n_coord_freq = 50000
n_hl_freq    =  5000
n_E_wave_period, E_wave_type

n_E_wave_period is the period of alternating external electric field which have external_E_field amplitude. Set 'triangular_sin' or 'triangular_cos' to E_wave_type. See example files in src/27example-BaTiO3-new-param-E_wave/ . n_E_wave_period must be a positive integer divisible by 4 for 'triangular_sin' or 'triangular_cos'. You can also use 'ramping_off' and 'ramping_on'. n_E_wave_period must no be divisible by 4 for 'ramping_on' and 'ramping_off' cases.

Figure 1: Four options of E_wave_type.

You can restart and continue a hysteresis-loop calculation from the n_hysteresis_loop_continue-th iteration.

n_hysteresis_loop_continue = 20000

If this tag is set as

coord_directory = 'coord_files'

.coord file will be written into coord_files/foo123.0000020000.coord . For example, use this tag to store huge .coord files into a local HDD.

If you do not want to get any .coord files, set it as

coord_directory = 'never'

If this tag is set as

distribution_directory = 'distribution_files'

.distribution and .distribution3d files will be written into distribution_files/foo123.distribution and distribution_files/foo123.distribution3d . For example, use this tag to store the huge .distribution3d file into a local HDD.

If you do not want to get any .distribution and .distribution3d files, set it as

distribution_directory = 'never'

This can also reduce computational time.


Vector of the external electric field in the unit of [V/Angstrom].

external_E_field = 0.00 0.00 -0.01

External electric field will be constant in this value when method is 'md', 'vs' or 'lf'. Or, if E_wave_type is specified, E will be altering in period n_E_wave_period.

P_kappa2, P_alpha, P_gamma

Coefficients for the 4-th order polynomial.

P_kappa2 =    5.502  [eV/Angstrom^2] # P_4(u) = kappa2*u^2 + alpha*u^4
P_alpha  =  110.4    [eV/Angstrom^4] #  + gamma*(u_y*u_z+u_z*u_x+u_x*u_y),
P_gamma  = -163.1    [eV/Angstrom^4] # where u^2 = u_x^2 + u_y^2 + u_z^2
P_k1, P_k2, P_k3

Coefficients for the 6-th order part


Coefficient for u^8


Short-range inter-site interaction coefficients

j = -2.648 3.894 0.898 -0.789 0.562 0.358 0.179    j(i) [eV/Angstrom^2]
B11, B12, B44

Elastic Constants

B11 = 126.
B12 =  44.9
B44 =  50.3  [eV]
B1xx, B1yy, B4yz

Elastic Coupling

B1xx = -211.   [eV/Angstrom^2]
B1yy =  -19.3  [eV/Angstrom^2]
B4yz =   -7.75 [eV/Angstrom^2]

Two integers for Marsaglia-Tsang 64-bit universal RNG (random number generator).

seed = 123384 8984847

Its default values are 123456789 987654321. It may be a good idea to generate the seed with jot(1) command.

$ jot -r 2 1 2147483648

You can also use $((RANDOM)) of bash(1).

$ echo $((RANDOM*65536+RANDOM)) $((RANDOM*65536+RANDOM))
2146723680 1863789493
init_dipo_avg, init_dipo_dev

Values for initial random setting of dipoles. The normal-distributed random number generator N(mu,sigma^2) is used.

For example, if you only have four dipoles and

init_dipo_avg = 0.10 0.00 0.00   [Angstrom]  # Average (mu)      of initial dipole displacements
init_dipo_dev = 0.02 0.02 0.02   [Angstrom]  # Deviation (sigma) of initial dipole displacements,

the dipoles may be

0.09  0.01  0.02
0.11 -0.02  0.01
0.08  0.02 -0.03
0.12 -0.01  0.00 .

If there is a .restart file in the same directory, there values are not used and the .restart file will be read.


Effective charge per site.

Z_star     = 9.956
epsilon_inf   = 5.24

This is for debug. You do not have to set this tag. Arbitrary kappa for Ewald summation in src/dipole-dipole-long-range.F.

kappa = 0.15


plot_dispersion = .true.

data files for plotting dispersion relation will be written. Default is .false..

For plot_dispersion = .true., L_x, L_y and L_z must be a same even number.

L = 32 32 32

may be a good idea.

See an example of src/08example-BaTiO3-dispersion/ or try

$ tar xf SOMEWHERE/feram-X.YY.ZZ.tar.xz
$ cd feram-X.YY.ZZ
$ mkdir builddirectory
$ cd builddirectory
$ ../configure && make -j4
$ make check TESTS='dispersion_check.sh dispersion_32x32x32.inhomo-K.gp'
$ less dispersion_32x32x32.feram   # the input file
$ less dispersion_32x32x32.gp      # the GNUPLOT script
$ gv src/dispersion_32x32x32.long+short.energy.eps &
$ gv src/dispersion_32x32x32.long+short.eps &
$ gv src/dispersion_32x32x32.long+short.interaction.eps &
$ gv src/dispersion_32x32x32.inhomo-K.eps &
                  # Previews of plots. You may also use evince(1) instead of gv(1).

Compare the plots with Fig.3 in [Takeshi Nishimatsu el al.: Phys. Rev. B 82, 134106 (2010), http://dx.doi.org/10.1103/PhysRevB.82.134106 ] and Fig.5 in parameters/parameters.html.


If feram is executed as

$ feram foo.feram bar.feram

and, in bar.feram,

continue = .true.

is declared, simulation of bar.feram uses results of foo.feram as initial configurations. The .restart file will not be used. Default value of continue is .false..

.restart input file

If you executed feram with a parameter file "foo123.feram" as

$ feram foo123.feram

feram searches for the "foo123.restart" file as an initial configurations. If there is no such file, feram starts from a random configurations.

.quadratic input file

Quadratic inter-atomic force constant (IFC) matrices can be read from .quadratic file. Set all j1-j7 to ZERO. Use P_kappa2 for *elevation*. (Lx/2+1)*Ly*Lz 3x3 matrices are required. See src/25example-BaTiO3-read-quadratic/ for more details.

.defects input file

If you executed feram with a parameter file "foo123.feram" as

$ feram foo123.feram

feram searches for the "foo123.defects" file as information of defects. If there is no such file, feram do not introduce defects.

The .defects input file gives positions and fixed dipole moments u for simulations of defects in ferroelectrics.

0  0  0   0.1 -0.1  0.7
1  0  0   0.1  0.1  0.7
2  0  0   0.1  0.2  0.7
position_x position_y position_z    fixed_u_x fixed_u_y fixed_u_z [Angstrom]

Use src/feram_defects_maker.rb to make this .defects file.

.localfield input file

If you executed feram with a parameter file "foo123.feram" as

$ feram foo123.feram

feram searches for the "foo123.localfield" file as information of local electric field. Format of .localfield file is:

 4  6 13   0.00 0.00 0.10
ix iy iz   Ex   Ey   Ez

The unit of local field is V/Angstrom. If feram cannot find both ionic.configuration and .localfield, it will set the local field to ZERO for all sites.

Use attached feram_random_field_generator for generating random fields.

$ ./feram_random_field_generator Lx Ly Lz  mu_x mu_y mu_z  sigma_x sigma_y sigma_z  seed1 seed2
$ ./feram_random_field_generator 32 32 243  0.0 0.0 0.0  0.02 0.02 0.02  123456789 987654321

.modulation file

.modulation file gives acoustic modulation, a(R), of each site. Energy from acoustic modulation is

E_acoustic_modulation = modulation_constant Sigma_R_alpha eta_alpha(R) a(R) .

FFTW wisdom file

By preparing an FFTW wisdom file in your current directory, you can reduce the total computational time. Sometimes, however, this trick fails and the computational times increases. You should optimize the computational condition with trial-and-error.

You can generate an FFTW wisdom file with `feram_fftw_wisdom` command, which will be built in the src/ directory besides `feram` command. In the src/18example-benchmark/ directory, for example, you can generate an FFTW wisdom file of system size Lx*Ly*Lz=32*32*243 for the forward.feram input file,

$ w   # Before executing feram_fftw_wisdom, confirm that load average is zero and there is no background job.
$ OMP_NUM_THREADS=6 ../feram_fftw_wisdom 1000 32 32 243 3
  It takes a few minutes or a few hours depending on the system size.
  '1000' in the command line arguments is the number of iterations of
  FFT benchmarks. So, do not set a large number for it.
$ mv wisdom_new wisdom
$ OMP_NUM_THREADS=6 ../feram forward.feram
$ grep FFTW_WISDOM forward.log
  feram_common.F: 47: FFTW_WISDOM: Successfully imported FFTW wisdom in current directory.
$ grep TIMING_REPORT forward.log

Command line arguments for feram_fftw_wisdom are

feram_fftw_wisdom <N_TIMES> <Lx> <Ly> <Lz> <padding_y>

Check http://www.fftw.org/doc/Other-Important-Topics.html for more details on FFTW wisdom. If you are using the MKL for FFTW3 wrappers, you do not have to prepare the wisdom file.

Output files

If you execute the feram like `feram foo123.feram`, filenames of output files are starting from 'foo123', e.g. foo123.avg, foo123.0000050000.coord, foo123.param.gp, foo123.dipole-dipole-long.dat, etc.


To foo123.log, feram reports the energies of each iteration. The "verbose" tag determines how much messages will be written into the .log file.

For debugging, it may be a good idea to use "verbose = 2" and plot energies in the .log file. For example:

$ ./feram zzznp.feram
$ gnuplot zzznp.gp
$ ./feram zzzlf.feram
$ gnuplot zzzlf.gp

In the end of the .log file, you will have a TIMING_REPORT.

   T =  177.1 [K] END ============================================================================
  molecular-dynamics.F:158: END: molecular_dynamics().
feram_common.F:104: DATE_AND_TIME: 2016-06-24T23:16:26.498+0900
 t_initialization, t_simulation,  t_total, n_threads
       0.47          39.77          40.24    2     #TIMING_REPORT
feram_common.F:117: FINISH: love && peace && free_software


After n_thermalize iterations of thermalization, feram averages propertqies for n_average iterations. For example, by combining .avg files of some calculations altering temperature (kelvin tag), you can see temperature dependence of properties. Please see the end of src/average_module.F.


Distribution function of u_alpha(R) (alpha=x,y,z)

D_alpha(u) = (1/N) Sum_R delta(u_alpha(R)-u)

are averaged in n_average iterations and reported in foo123.distribution. From 1st to 4th column, u D_x(u) D_y(u) D_z(u) are written. Please see src/average_module.F.


Distribution function of vector u(R)

D(u) = (1/N) Sum_R delta(u(R)-u)

are averaged in n_average iterations and reported in foo123.distribution3d.


In the case of method = 'hl', a simulation of hysteresis loop, one .hl file will be created instead of a .avg file. In .hl file, properties will be recorded in every n_hl_freq iterations. Please see src/hysteresis_loop.F.


Snapshot of dipoles at nnnnnnnnnn-th iteration. Snapshots are taken at evry n_coord_freq iterations. This file can be visualized with feram_slicer.rb (Ruby script) and feram_cross_section_q.sh, feram_cross_section_p.sh, and feram_cross_section_dVddi.sh (Bourne Shell scripts). feram_cross_section_p.sh and feram_cross_section_dVddi.sh are symbolic links or copy of feram_cross_section_q.sh. Output of this file can be suppressed with

coord_directory = 'never'

in the .feram input file.


Visualized slice of a film at z=Lz/4, at the nnnnnnnnnn-th iteration and at evry n_coord_freq iterations. Output of this .slice.eps file can be suppressed with

slice_directory = 'never'

in the .feram input file.


Parameter file for GNUPLOT. rename it to param.gp, then use it with cross_section_q.sh etc.


foo123.dipole-dipole-long+short.dat, foo123.dipole-dipole-long.dat and foo123.dipole-dipole-short.dat will be made if plot_dispersion = .true. and contain optical phonon dispersion. foo123.inhomo-K.dat contains acoustic phonon dispersion.

A line in foo123.dipole-dipole-long+short.dat may be

   0.59375  16   3   0 -0.2022348055959993E+001  0.1270295721445465E+001  0.3455139983790274E+002

that is

   length  index_of_k-point  3_eigenvalues_of_IFC_matrix.

Eigenvectors are also written in the lines starting with "#".

For more details, see descriptions of plot_dispersion above or source code of src/dipole-dipole.F and src/print-eigenvalues.F .

local.field output file

The calculated local field from the given ionic configuration (ionic.configuration file) is stored in this file. This file will be reloaded in consequent temperature calculations and the other calculations.


Tools for analyzing and plotting are prepared in the src/ directory. They are start with `feram_`. They will be also installed into /usr/local/bin/ by default, if you invoke `make install`. Most of them are Ruby scripts and shell scripts. Some of them are written in Fortran language. Basically, usages are described in the head of source code.

feram_transition_detector.rb and feram_transition_sorter.rb

`feram_transition_detector.rb` reads cooling.avg or heating.avg and reports transition temperatures by marking `Tc`, then `feram_transition_sorter.rb` can sort the results. See src/34example-BST/README.md for more details.


feram_slicer.rb reads a .coord file and visualizes the supercell as a slice with +z-polarized red empty squares and -z-polarized blue filled squares. feram_slicer.rb is a Ruby script file which writes an Encapsulated PostScript file (EPSF, .eps file) directly. You can specify the z-height as the second command-line argument.

feram_slicer.rb zzz16.0000000060.coord
feram_slicer.rb zzz16.0000000060.coord 8   # ---> zzz16.0000000060-slice-z008.eps


feram_cross_section_q.sh reads a .coord file and draws a cross section with arrows representing dipoles. Prepare param.gp as a symlink. If you specify zzz16.0000000060.coord, you will get an EPS file, zzz16.0000000060-q-x.eps for example.

   ln -s zzz16.param.gp param.gp
   feram_cross_section_q.sh zzz16.0000000060.coord   # You will get zzz16.0000000060-q-x.eps.
   feram_cross_section_q.sh 150K.0000050000.coord 4.0 10 y 0.9
Usage: ./feram_cross_section_q.sh coord-file [FACTOR] [CONST_Alpha] [Alpha] [ratio] [max_z]
 [FACTOR]      u*[FACTOR] will be the length of each arrow.
 [CONST_Alpha] Visulaize cross sections of alpha=[CONST_Alpha]. Default value: 8.
 [Alpha]       Alpha=x,y,z. Default value: x.
 [ratio]       Optional argument to keep the shape of unitcell square exactly.
               You may want to use this argument when Lx=Ly!=Lz. Default value: 0.7231.
 [max_z]       It is useful for vertical cross section of a thin-film.


You can find input files and their simulated results in following directories in the software package.


Energies vs MD step in cooling simulation

Figure 2: Energies vs MD steps in cooling simulation.

Temperature dependence of strain and susceptibility of bulk BaTiO3

Figure 3: Temperature dependence of three components of strain and relative dielectric constant of bulk BaTiO3. Three phase transitions are clearly seen. Susceptibilities are calculated from fluctuations of dipoles.

Temperature dependence of hysteresis loops for bulk BaTiO3

Figure 4: Temperature dependence of hysteresis loops for bulk BaTiO3. The phase transition from paraelectric phase to ferroelectric phase is clearly seen.

Hysteresis loops for epitaxially constrained and "free" BaTiO3 film capacitors

Figure 5: Calculated hysteresis loops for capacitors with (a) epitaxially constrained films, and (b) "free" films of various thickness l and with dead layer d at 100 K.

Mailing list

Development of feram

GNU Autotools

GNU Autotools (autoconf and automake) are required for development of feram. Most Linux distributions have packages of GNU Autotools. In the case of Debian or Ubuntu, you can install GNU Autotools with

$ sudo apt-get install autoconf automake libtool autoconf-doc libtool-doc

Getting current source tree, compilation and execution

Get current source tree with svn(1). Anonymous checkout can be done without --username=YourUsername option.

$ svn checkout --username=YourUsername https://svn.code.sf.net/p/loto/code/feram/trunk feram-trunk
$ cd feram-trunk
$ ls -l
$ autoreconf -v
$ automake --add-missing
$ autoreconf -v
$ ls -l
$ ./configure --help
$ ./configure
$ emacs src/feram.F   # Edit any source code here.
$ make
$ make check

Useful svn commands:

$ svn --help
$ svn update
$ svn stat
$ svn diff
$ svn log

How to execute:

$ ./feram foo1.feram foo2.feram foo3.feram &
$ mpiexec -np 3 ./feram_mpi foo1.feram foo2.feram foo3.feram &


index.en.html is automatically generated from README.en with ULMUL http://t-nissie.users.sourceforge.net/ULMUL/ . After editing README.ja (Japanese) or README.en (English), transform them with

$ make index.en.html index.ja.html



Copyright © 2007-2016 by Takeshi Nishimatsu

feram is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY. You can copy, modify and redistribute feram, but only under the conditions described in the GNU General Public License (the "GPL"). For more detail, see COPYING.

You can make a DONATION to Takeshi Nishimatsu through Tohoku University http://www.rpip.tohoku.ac.jp/main/kihukin.html (in Japanese).

The author is grateful if you would kindly refer the name of this program "feram" (all lowercase letters) and its URL http://loto.sourceforge.net/feram/ in your papers. Also, please cite our articles, [Takeshi Nishimatsu, Umesh V. Waghmare, Yoshiyuki Kawazoe, and David Vanderbilt: Phys. Rev. B 78 (2008) 104104] and [Jaita Paul, Takeshi Nishimatsu, Yoshiyuki Kawazoe, and Umesh V. Waghmare: Phys. Rev. Lett. 99 (2007) 077601].

Cited figures

Authors and/or publishers of the papers have their copyrights for the following cited figures.


Takeshi Nishimatsu (t_nissie{at}yahoo.co.jp)


This project was partially supported by:

http://loto.sourceforge.net/feram/ is hosted by SourceForge.net Logo .

Copyright © 2017 Takeshi Nishimatsu