File Formats and Interpretation of Data

FAQ on interpretation of data:

• Why is my total energy so much different from a Gaussian calculation?
• What is the meaning of the energies reported in a molecular dynamics simulation?
• What do GNMAX, GNORM, and CNSTR mean in a geometry optimization?
• All IR intensities are zero in VIB.LOG when I try to calculate the IR of NH$_4^+$ ion by CPMD
• Why can some atoms move far away from the central box in a long MD simulation of bulk liquid using periodic boundary conditions?
• Why did the electron energy continuously increase during the MD simulation of my metal?
• What is the meaning of the files APT, POLARIZATION and POLARIZABILITY resulting from my linear response calculation of Raman spectra?
• What is the meaning of columns 2 to 7 in the DIPOLE file?
• How do I use the binary RESTART file across different platforms with different encoding?

Q: Why is my total energy so much different from a Gaussian calculation?

A: With CPMD you are using pseudopotentials to describe the atoms. Since the total energy describes only the interactions between the pseudocores and the valence electrons (and some core electrons in the case of so-called semi-core pseudopotentials), you are missing the contribution of the core electrons and the full core charges of a regular all-electron calculation. Energy differences between two configurations, on the other hand, should be comparable, provided you use the same number of atoms, the same plane wave cutoff, the same pseudopotentials, and the same supercell geometry in the CPMD calculation.

Q: In a molecular dynamics simulation, CPMD prints out a list of energies for each integration step. Does anyone know the meaning of the individual values.

A: Some explanations to the energy terms:

EKINC

fictitious kinetic energy of the electrons in a.u. this quantity should oscillate but not increase during a simulation.

TEMPP

Temperature of the ions, calculated from the kinetic energy of the ions (EKIONS).

EKS

Kohn-Sham energy (the equivalent of the potential energy in classical MD).

ECLASSIC

= EKS + EKIONS

EHAM

= ECLASSIC + EKINC. Hamiltonian energy, this is the conserved quantity, depending on the time step and the electron mass, this might oscillate but should not drift.

DIS

mean square displacement of the ions with respect to the initial positions. Gives some information on the diffusion.

You can modify the list of individual energies to be displayed with the PRINT ENERGY keyword.

Q: What do GNMAX, GNORM and CNSTR in a geometry optimization mean?

A:

GNMAX

max${}_{I,a}$ ($\vert F_{Ia}\vert$ ) = largest absolute component ($a=x,y,z$ ) of the force on any atom $I$ .

GNORM

$\left< F_{I}^2\right>_{I}$ = average force on the atoms $I$

CNSTR

max ${}_{I,a} { F^{constr}_{Ia} }$ = largest absolute component ($a=x,y,z$ ) of force due to constraints on any atom $I$ .

Q: I found all the IR intensities in VIB.log file are zero when I try to calculate the IR of NH${}_4^+$ ion by CPMD.

 Harmonic frequencies (cm**-1), IR intensities (KM/Mole),
 Raman scattering activities (A**4/AMU), Raman depolarization ratios,
 reduced masses (AMU), force constants (mDyne/A) and normal coordinates:
                     1                      2                      3
                    ?A                     ?A                     ?A
 Frequencies --   142.9800               188.9340               237.2614
 Red. masses --     0.0000                 0.0000                 0.0000
 Frc consts  --     0.0000                 0.0000                 0.0000
 IR Inten    --     0.0000                 0.0000                 0.0000
 Raman Activ --     0.0000                 0.0000                 0.0000
 Depolar     --     0.0000                 0.0000                 0.0000
 Atom AN      X      Y      Z        X      Y      Z        X      Y      Z
   1   7     0.00   0.00   0.00     0.00   0.00   0.00     0.00   0.00   0.00
   2   1     0.00  -0.35  -0.50    -0.35   0.00   0.00    -0.50   0.00   0.00
   3   1     0.00  -0.35   0.50    -0.35   0.00   0.00     0.50   0.00   0.00
   4   1     0.00   0.35   0.00     0.35   0.00   0.50     0.00  -0.50   0.00
   5   1     0.00   0.35   0.00     0.35   0.00  -0.50     0.00   0.50   0.00
                     4                      5                      6

A: That's not a problem of your calculation. The keyword VIBRATIONAL ANALYSIS does not calculate intensities. The calculation of intensities is currently not possible in CPMD. The intensities in the 'VIBx.log' files are arbitrarily set to zero. The entries have to be there so that visualisation programs, that are able to read output of the Gaussian program, can be also used to visualize the CPMD results.

Q: I am trying to simulate a bulk liquid in CPMD and supposing that periodic boundary conditions are built into the program. But after several thousand MD steps, I found some particles are far away from the central simulation box.

Why it is so if periodic boundary conditions (PBC) on particle coordinates are imposed in all three directions?

A: If you are not using the

SYMMETRY
  0

options your calculations are actually using periodic boundary conditions (PBC). PBC are imposed within CPMD for all calculations. However, the particle positions are not folded back to the original computational box. The reason for this is that most people prefer to have ``smooth'' trajectories without jumps of particles. This allows for easier tracking of special particles and nicer graphics. In addition it is easy (with a little script) to apply PBC afterwards yourself, if needed.

Q: I am trying to simulate a bulk sodium and I found electron energy is increasing continuously and it is in the range of 0.07 a.u. at the end of 20000 steps.

A: Sodium is a metal, and therefore missing an important feature that allows for stable CP dynamics: the band gap. Using Nosé thermostats (on electrons and ions) it might still be possible to perform meaningful CP simulations [205].

The choice of parameters for the thermostats, however, will be nontrivial, highly system dependent and require extensive testing. Without thermostats you will have strong coupling between electronic degrees of freedom and ionic degrees of freedom. Adiabaticity is not maintained and a steady increase of the fictitious kinetic energy will occur.

Q: I have computed RAMAN by LINEAR RESPONSE, and get three files: APT, POLARIZATION and POLARIZABILITY with lots of data in these files. I want to know the meaning of the data, please give me some answer in detail.

A: The POLARIZABILITY file simply contains the polarizability tensor of the whole system in atomic units. The POLARIZATION file contains the total dipole moment (electronic + ionic) of the whole system in atomic units. As for the file APT, it contains the atomic polar tensors for each atom in the system. The atomic polar tensor is the derivative of the forces on the atoms with respect to an applied external electric field. Equivalently it is, from a Maxwell relation, the derivative of the total dipole of the system with respect to the nuclei positions. It is thus an important ingredient of the calculation of infrared spectra intensities, for example used in an harmonic approximation. The trace of this tensor is the so-called Born charge of the considered atom. The data is arranged in the following order (still in a.u.): the APT tensor is $\frac{\mathrm{d}F_{I,i}}{\mathrm{d}E_{j}}$ where $F_{I,i}$ is the force on atom I along $i=x,y,z$ and $E_j$ is the electric field along $j=x,y,z$ . $(I,i)$ are the indices of the $3N$ atoms lines in the APT file, one atom after the other, and $j$ is the column index in the APT file.

Q: I was wondering what columns 2 to 7 in the DIPOLE file correspond to? When I run CPMD v.3.5.1, columns 2 to 4 come out identical to columns 5 to 7 respectively. When I run with CPMD v.3.4.1, the columns come out different. Is there an explanation for this?

A: Columns 2 to 4 in the DIPOLE file are the electronic contribution to the dipole moment, columns 5 to 7 are the total (electronic + ionic) dipole moment. All dipole moments are divided by the volume of the box.

In CPMD version 3.5.1 we have changed the reference point of the calculation. Now the reference point is chosen such that the ionic contribution is zero and the electronic contribution minimal (=total dipole). This avoids a problem that occasionally was seen in older versions. The electronic dipole is calculated modulo(2$\pi$ /L). Now if the electronic dipole became too large, because the ionic contribution was large (bad choice of reference point) the total dipole made jumps of 2$\pi$ .

Q: As you know, the cpmd RESTART file is saved as binary. But I want to change it to ASCII and vice versa, because I use several machines of different architecture, for example COMPAQ, IBM, and LINUX machine. Please help me with any comments.

A: The code to read and write the RESTART file is in the files rv30.F and wv30.F. Feel free to implement an ASCII version of the restart, but be aware that the file will be huge.

But you may not need to do that. Let's say you decide to use big-endian binary encoding (this is what e.g. IBM, Sun and SGI machines do natively).
With Compaq machines there is a compiler flag, -convert, which you could set to big_endian (we only have here linuxalpha, but the compaq compiler should be essentially the same).

On a Linux PC you can use the use the -Mbyteswapio or the -byteswapio flag, if you have the PGI compiler.

For the Intel compiler (ifc/ifort/efc) you simply set the environment variable F_UFMTENDIAN to big (i.e.
'export F_UFMTENDIAN=big'
if you are in a bourne/korn shell and
'setenv F_UFMTENDIAN big'
if you are in a (t)csh).

Now even your cpmd executables will read and write big-endian restart files.

Check your compiler documentation for more details (search for endian).