Personal tools
You are here: Home Tutorial The Basics of CPMD - Geometry Optimization and MD

The Basics of CPMD - Geometry Optimization and MD

In this section, we will inspect how to perform a geometry optimization and how to setup and analyze a Car-Parrinello MD.

Geometry Optimization

 

A geometry optimization is not much else than repeated single point calculations, where the positions of the atoms are updated according to the forces acting on them. The required changes in the input file are rather small (5-h2-geoopt.inp):

 

&CPMD
 OPTIMIZE GEOMETRY XYZ
 CONVERGENCE ORBITALS
  1.0d-7
 CONVERGENCE GEOMETRY
  1.0d-4
&END

 

We have replaced WAVEFUNCTION with GEOMETRY and added the suboption XYZ to have CPMD write a 'trajectory' of the optimization in a file name GEO_OPT.xyz (so it can be visualized later). Also we specify the convergence parameter for the geometry.

Now again start the CPMD program:
   

cpmd.x 5-h2-geoopt.inp> 5-h2-geoopt.out

 
This run should take a little longer, than the previous one, since we have to do multiple wavefunction optimizations.

 

 OPTIMIZATION OF IONIC POSITIONS

[...]

 CONVERGENCE CRITERIA FOR GEOMETRY OPTIMIZATION:     1.000000E-04
 GEOMETRY OPTIMIZATION BY GDIIS/BFGS
   SIZE OF GDIIS MATRIX:                                        5
GEOMETRY OPTIMIZATION IS SAVED ON FILE GEO_OPT.xyz
 EMPIRICAL INITIAL HESSIAN (DISCO PARAMETRISATION)

 

As you can see from the first part of the output file (5-h2-geoopt.out), CPMD has recognized the job type, our convergence parameter and the request to write a GEO_OPT.xyz file.

 

 ================================================================
 =                  GEOMETRY OPTIMIZATION                       =
 ================================================================
 NFI      GEMAX       CNORM           ETOT        DETOT      TCPU
  EWALD| SUM IN REAL SPACE OVER                     1* 1* 1 CELLS
   1  3.816E-02   2.886E-03      -1.096898   -1.097E+00      1.28
   2  8.628E-03   1.041E-03      -1.130803   -3.391E-02      1.33

[...]

  10  1.871E-08   5.247E-09      -1.132460   -8.509E-13      1.43

 RESTART INFORMATION WRITTEN ON FILE                  ./RESTART.1

   ATOM          COORDINATES            GRADIENTS (-FORCES)
   1  H  8.2600  7.5589  7.5589  -1.780E-02  9.179E-17  7.909E-17
   2  H  6.8578  7.5589  7.5589   1.780E-02  1.596E-16  1.396E-16
 ****************************************************************
 *** TOTAL STEP NR.    10           GEOMETRY STEP NR.      1  ***
 *** GNMAX=  1.779864E-02                ETOT=     -1.132460  ***
 *** GNORM=  1.027605E-02               DETOT=     0.000E+00  ***
 *** CNSTR=  0.000000E+00                TCPU=         13.63  ***
 ****************************************************************
   1  5.012E-03   9.718E-04      -1.131471    9.887E-04      1.34
   2  4.287E-04   1.613E-04      -1.132846   -1.375E-03      1.35
   3  1.489E-04   3.429E-05      -1.132883   -3.659E-05      1.33

 

In the following output you can see, that an almost identical wavefunction optimization takes place. After printing the positions and forces of the atoms, however, you see a small report block and then another wavefunction optimization starts. The numbers for GNMAX, GNORM, and CNSTR stand for the largest absolute component of the force on any atom, average force on the atoms, and the largest absolute component of a constraint force on the atoms respectively.

 

   ATOM          COORDINATES            GRADIENTS (-FORCES)
   1  H  8.2854  7.5589  7.5589   9.965E-05  1.105E-16  9.709E-17
   2  H  6.8324  7.5589  7.5589  -9.965E-05  1.835E-16  1.392E-16
 ****************************************************************
 *** TOTAL STEP NR.    36           GEOMETRY STEP NR.      5  ***
 *** GNMAX=  9.965023E-05 [5.98E-05]     ETOT=     -1.132896  ***
 *** GNORM=  5.753309E-05               DETOT=    -1.423E-08  ***
 *** CNSTR=  0.000000E+00                TCPU=          6.76  ***
 ****************************************************************
 ================================================================
 =              END OF GEOMETRY OPTIMIZATION                    =
 ================================================================

At the end of the geometry optimization, you can see that the forces and the total energy have significantly decreased from their start values as it is to be expected.

 

Car-Parrinello Molecular Dynamics

 

Based on the previously calculated electronic structure, we can now also start a 'real' Car-Parrinello calculation. Note that although you can start a CP-MD run from a non-converged wavefunction (e.g. by not restarting from a pre-optimized wavefunction), you will be far away from the Born-Oppenheimer surface, and thus your result will be unphysical.

 

For the CP-MD job you need a new input file, 6b-h2-md-4au.inp, which should be copied into the same directory, where you started the wavefunction optimization run. If you compare it to the previous input files, you will find, that the only changes are again only in the &CPMD section of the input file.

&CPMD
 MOLECULAR DYNAMICS CP
 RESTART WAVEFUNCTION COORDINATES LATEST
 TRAJECTORY XYZ
 TEMPERATURE
  50.0D0
 MAXSTEP
  200
 TIMESTEP
  4.0
&END

The keyword MOLECULAR DYNAMICS CP defines the job type. Furthermore we tell the CPMD program to pick up the previously calculated wavefunction and coordinates from the latest restart file (which is named RESTART.1 by default). MAXSTEP limits the MD to 200 steps and the equations of motion will be solved for a time step of 4 atomic units (~0.1 femtoseconds). The temperature of the system will be initialized to 50K via the TEMPERATURE keyword (note that this is no thermostatting). The keyword TRAJECTORY XYZ will have CPMD write the coordinates additionally to a file TRAJEC.xyz, which can be easily visualized with many molecular viewer programs.


Now start the CPMD program once more: 

 

cpmd.x 6b-h2-md-4au.out> 6b-h2-md-4au.out

 

This run should be completed in a few minutes. The output of the CPMD program is now in the file . There also are some more new files (TRAJEC.xyz, TRAJECTORY, ENERGIES), but we'll have a closer look at output file first.

 

 CAR-PARRINELLO MOLECULAR DYNAMICS

PATH TO THE RESTART FILES: ./
RESTART WITH OLD ORBITALS
RESTART WITH OLD ION POSITIONS
RESTART WITH LATEST RESTART FILE
ITERATIVE ORTHOGONALIZATION
MAXIT: 30
EPS: 1.00E-06
MAXIMUM NUMBER OF STEPS: 200 STEPS

 

The header is unchanged up to the point where the settings from the &CPMD section are printed. As you can see, the program has recognized the RESTART and the MAXSTEP keywords. (NOTE: in the CPMD code atoms are sometimes referred to as ions, which may be sometimes confusing. This is due to the pseudopotential approach, where you integrate the core electrons into the (pseudo)atom which then could be also described as an ion.)

 

 TIME STEP FOR ELECTRONS:                                  4.0000
 TIME STEP FOR IONS:                                       4.0000
 TRAJECTORIES ARE SAVED ON FILE
 TRAJEC.xyz IS SAVED ON FILE
 ELECTRON DYNAMICS: THE TEMPERATURE IS NOT CONTROLLED
 ION DYNAMICS:      THE TEMPERATURE IS NOT CONTROLLED

 

This part of the output tells us, that the TIMESTEP 4.0 keyword was recognized (the default is 5.0 a.u., cf. the wavefunction output file), that the trajectory will be recorded and that there will be no temperature control, i.e. we will do a microcanonical (NVE-ensemble) simulation.

 

 RV30| WARNING! NO WAVEFUNCTION VELOCITIES

 RESTART INFORMATION READ ON FILE                     ./RESTART.1

 

Here we get notified, that the program has read the requested data from the restart file. The warning about the missing wavefunction velocities is to be expected, since they will only be available when the restart was written by a previous Car-Parrinello MD run.

 

       NFI    EKINC   TEMPP           EKS      ECLASSIC          EHAM         DIS    TCPU
         1  0.00000    49.4      -1.13289      -1.13266      -1.13266   0.207E-05    1.37
         2  0.00000    47.7      -1.13289      -1.13266      -1.13266   0.817E-05    1.36
         3  0.00001    45.7      -1.13288      -1.13267      -1.13266   0.181E-04    1.37
         4  0.00001    43.5      -1.13288      -1.13267      -1.13266   0.314E-04    1.37
         5  0.00002    41.5      -1.13287      -1.13268      -1.13266   0.481E-04    1.37
         6  0.00002    39.8      -1.13287      -1.13268      -1.13266   0.677E-04    1.36
         7  0.00002    38.3      -1.13286      -1.13268      -1.13266   0.901E-04    1.34
         8  0.00002    37.1      -1.13286      -1.13268      -1.13266   0.115E-03    1.38
         9  0.00002    35.9      -1.13285      -1.13268      -1.13266   0.143E-03    1.36
        10  0.00002    34.8      -1.13284      -1.13268      -1.13266   0.173E-03    1.36
        11  0.00002    33.6      -1.13284      -1.13268      -1.13266   0.206E-03    1.39
        12  0.00002    32.4      -1.13283      -1.13267      -1.13266   0.240E-03    1.37
        13  0.00001    31.2      -1.13282      -1.13267      -1.13266   0.277E-03    1.37
        14  0.00001    30.0      -1.13281      -1.13267      -1.13266   0.314E-03    1.37
        15  0.00001    28.9      -1.13281      -1.13267      -1.13266   0.354E-03    1.36
        16  0.00001    27.8      -1.13280      -1.13267      -1.13266   0.394E-03    1.35
        17  0.00001    26.9      -1.13279      -1.13267      -1.13266   0.436E-03    1.37
        18  0.00001    26.1      -1.13279      -1.13267      -1.13266   0.478E-03    1.36
        19  0.00001    25.4      -1.13279      -1.13267      -1.13266   0.521E-03    1.37
        20  0.00001    25.0      -1.13278      -1.13267      -1.13266   0.564E-03    1.38
        21  0.00001    24.8      -1.13278      -1.13267      -1.13266   0.608E-03    1.37
        22  0.00001    24.8      -1.13278      -1.13267      -1.13266   0.652E-03    1.37
        23  0.00001    25.1      -1.13279      -1.13267      -1.13266   0.696E-03    1.36
        24  0.00001    25.5      -1.13279      -1.13267      -1.13266   0.741E-03    1.40
        25  0.00001    26.2      -1.13279      -1.13267      -1.13266   0.786E-03    1.36
       ...
       180  0.00002    43.1      -1.13288      -1.13268      -1.13266   0.321E-01    1.39
       181  0.00002    42.2      -1.13287      -1.13267      -1.13266   0.325E-01    1.36
       182  0.00001    41.1      -1.13287      -1.13267      -1.13266   0.328E-01    1.35
       183  0.00001    39.8      -1.13286      -1.13267      -1.13266   0.331E-01    1.36
       184  0.00001    38.4      -1.13285      -1.13267      -1.13266   0.335E-01    1.37
       185  0.00001    36.8      -1.13284      -1.13267      -1.13266   0.339E-01    1.37
       186  0.00001    35.2      -1.13284      -1.13267      -1.13266   0.342E-01    1.36
       187  0.00001    33.5      -1.13283      -1.13267      -1.13266   0.346E-01    1.37
       188  0.00001    32.0      -1.13282      -1.13267      -1.13266   0.349E-01    1.35
       189  0.00001    30.5      -1.13281      -1.13267      -1.13266   0.353E-01    1.37
       190  0.00001    29.3      -1.13281      -1.13267      -1.13266   0.357E-01    1.37
       191  0.00001    28.2      -1.13280      -1.13267      -1.13266   0.361E-01    1.36
       192  0.00001    27.4      -1.13280      -1.13267      -1.13266   0.365E-01    1.36
       193  0.00001    26.8      -1.13280      -1.13267      -1.13266   0.369E-01    1.35
       194  0.00001    26.4      -1.13279      -1.13267      -1.13266   0.372E-01    1.37
       195  0.00001    26.2      -1.13279      -1.13267      -1.13266   0.376E-01    1.35
       196  0.00001    26.1      -1.13279      -1.13267      -1.13266   0.380E-01    1.36
       197  0.00001    26.3      -1.13279      -1.13267      -1.13266   0.385E-01    1.38
       198  0.00001    26.6      -1.13279      -1.13267      -1.13266   0.389E-01    1.35
       199  0.00001    27.2      -1.13280      -1.13267      -1.13266   0.393E-01    1.36
       200  0.00001    28.0      -1.13280      -1.13267      -1.13266   0.397E-01    1.38

 

After some more output, we already discussed for the wavefunction optimization, this is now part of the energy summary for a Car-Parrinello-MD run.

 

The individual columns have the following meaning:

  • NFI           Step number (number of finite iterations)

  • EKINC      (fictitious) kinetic energy of the electronic (sub-)system

  • TEMPP      Temperature (= kinetic energy / degrees of freedom) for atoms (ions)

  • EKS          Kohn-Sham Energy, equivalent to the potential energy in classical MD

  • ECLASSIC  Equivalent to the total energy in a classical MD (ECLASSIC = EHAM - EKINC)

  • EHAM        total energy, should be conserved

  • DIS           mean squared displacement of the atoms from the initial coordinates.

  • TCPU         (CPU) time needed for this step.



At the end of the geometry optimization the hydrogen molecule was in the minimum of its potential. Now after starting the MD, we see that the initial kinetic energy added to the system is slowly converted into potential energy (cf. EKS) as the bond is elongated. After a while the molecule has reached the maximal elongation and the potential energy is converted back into kinetic energy (i.e. the temperature rises again). So we have a regular oscillation of the hydrogen molecule. You can also see, that a little bit of energy is transferred into the fictitious dynamic of the electronic degrees of freedom. For a meaningful Car-Parrinello MD this value has to be (and stay) very small (although for larger systems with more electrons, the absolute value of EKINC will be larger).

 

 ****************************************************************
 *                      AVERAGED QUANTITIES                     *
 ****************************************************************
                              MEAN VALUE       +/-  RMS DEVIATION
                                          [-^2]**(1/2)
 ELECTRON KINETIC ENERGY    0.130119E-04             0.380704E-05
 IONIC TEMPERATURE                 34.76                     7.57
 DENSITY FUNCTIONAL ENERGY     -1.132836             0.388578E-04
 CLASSICAL ENERGY              -1.132671             0.384979E-05
 CONSERVED ENERGY              -1.132658             0.583016E-07
 NOSE ENERGY ELECTRONS          0.000000              0.00000
 NOSE ENERGY IONS               0.000000              0.00000
 CONSTRAINTS ENERGY             0.000000              0.00000
 ION DISPLACEMENT           0.135129E-01             0.119249E-01
 CPU TIME                         1.3665

 

Finally we get a summary of some averages and root mean squared deviations for some of the monitored quantities. This is quite useful to detect unwanted energy drifts or too large fluctuations in the simulation.

 

If you want to visualize the motion of the hydrogen atoms, you can load the file TRAJEC.xyz directly into a molecular visualization program like gopenmol, molden, molekel, vmd or xmol.

 

Further Types of Jobs

 

There are several further types of calculations possible with CPMD, for example, but not limited to:

 

KOHN-SHAM ENERGIES
PROPERTIES
LINEAR RESPONSE
VIBRATIONAL ANALYSIS
ORBITAL HARDNESS
ELECTRONIC SPECTRA

 

These are, however, beyond the scope of this little introduction. Please check out, the rest of this tutorial, the CPMD manual, the CPMD mailing list archives, and other CPMD input examples (e.g. the CPMD test suite) for more information on how to perform them (correctly).

Document Actions