feram: molecular dynamics simulator for bulk and thin-film ferroelectrics and relaxors

feram logo mark

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 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.coord, 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.


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. Minimum messages will be written to the .log file.

verbose = 2

When 2 is given for verbose tag, messages will be written to .log file in every iterations.

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'

'hl' is for a simulation of hysteresis loop. In this case, tags of n_average and external_E_field have special meanings (see Fig. 1). .hl file,e.g. foo123.hl, will be created. See the directory of src/11example-BaTiO3-new-hysteresis-loop/.

Figure 1: Special meanings of n_average and external_E_field, in the case of method = 'hl', a simulation of hysteresis loop. After thermalization, the external electric field is gradually altered from external_E_field to -external_E_field.

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' togeter 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

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/ .


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 uiniversal 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
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..

.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

Qudratic 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/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 random_field_generator for generating ramdom fields.

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

ionic.configuration file

(Not implemented yet.) This file is required for the calculation of the local field.

0.0 0.0 0.0   +3.0
0.5 0.0 0.0   -2.0
0.0 0.5 0.0   -2.0
0.0 0.0 0.5   -2.0
position_x position_y position_z    effective_charge_of_the_ion

Position must be given in the unit of a0.

FFTW wisdom file

By preparing an FFTW wisdom file in your current directory, you can reduce initial computational time, i.e. computational time before the "Ready for simulation!!!" message. You can generate an FFTW wisdom file with `fft_check` 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 fft_check, confirm that load average is zero and there is no background job.
$ OMP_NUM_THREADS=6 ../fft_check 1000 32 32 243
  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.
$ w
$ OMP_NUM_THREADS=6 ../fft_check 100 300 300 300
$ sudo cp wisdom_new /etc/fftw/wisdom
$ rm wisdom wisdom_new
$ OMP_NUM_THREADS=6 ../feram epit300x300x300.feram &
$ grep FFTW_WISDOM epit300x300x300.log
 feram_common.F: 52: FFTW_WISDOM: FFTW wisdom is not in current directory.
 feram_common.F: 55: FFTW_WISDOM: Successfully imported FFTW system wisdom, /etc/fftw/wisdom.

Here, /etc/fftw/wisdom is the system wisdom. In this case, you can use this /etc/fftw/wisdom for both Lx*Ly*Lz=32*32*243 and Lx*Ly*Lz=300*300*300 in the OMP_NUM_THREADS=6 environment.

See src/fft_check.html for more details on fft_check.

Check http://www.fftw.org/doc/Other-Important-Topics.html for more details on FFTW wisdom.

If you are not using the FFTW library but other FFT libraries such as MKL, MATRIX/MPP, SSL II, ACML, etc., 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.320K0050000.coord, foo123.dipole-dipole-long.dat, etc.


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

"verbose = 0" (NOT IMPLEMENTED YET). The .log file will not be created.

"verbose = 1" is the default value for verbose tag. Minimum messages will be written to the .log file.

"verbose = 2". feram reports the energies of EVERY iterations.

For debugging, using "verbose = 2", src/plot.gp can plot the energies.

$ feram foo123.feram
$ gnuplot plot.gp


After n_thermalize iterations of thermalization, feram averages properties 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. This file can be visualized with slice.rb (Ruby script) and cross-section-q.sh, cross-section-p.sh, and cross-section-dVddi.sh (Bourne Shell scripts). cross-section-p.sh and cross-section-dVddi.sh are symbolic links or copy of cross-section-q.sh.


Parameter file for GNUPLOT. rename it to param.gp, then use it with cross-section-q.sh etc.





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.


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://ulmul.rubyforge.org/ . After editing README.ja (Japanese) or README.en (English), transform them with

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



Copyright © 2007-2014 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 and 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], in your papers.

Cited figures

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


Takeshi Nishimatsu (t-nissie{at}imr.tohoku.ac.jp)


This project is/was partially supported by:

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

Copyright © 2014 Takeshi Nishimatsu