Mod:Hunt Research Group/w63Cl cpmd

63 water molecules + 1 Cl- CPMD simulation (Notes by Ling)

1. How to judge the simulations goes well:

a. Temperature is about right (if one starts from a random configuration, and does not use thermostats, typically T goes to thousand degress in just a few steps, so the dynamics becomes unstable).

b. Monitor EKINC (kinetic energy of CP degrees of freedom) is increasing a bit, but not too rapidly.

It is the kinetic energy of the fictitious Car-Parrinello electrons. When it is zero, it means that the system is in the ground state, when non-zero it means that the system is oscillating around the ground state. If EKINC keep increasing, rather than stabilising to a given value, the system is moving away from the ground state, and the forces are no longer computed accurately.

It is possible to control EKINC using the fictitious mass of the electron in the input.

EKINC increase from 0.001 at step 1 and stabilise around 0.029 at step 1000.

Initially it should increase a bit, because energy is being redistributed between the atoms and the fictitious set of electrons. At the certain point, when the two systems are in thermal equilibrium, EKINC should become more or less constant with time (with some oscillations).

Indeed, at the beginning of the simulation some of the temperature of the atoms is transferred to the Car-Parrinello electrons.

c. Monitor EHAM, which should be conserved on average, provided you are not using a thermostat. (Indeed EHAM seems to be oscillating within 10^-5 hartree, so it is doing well.)

d. AVERAGED QUANTITIES at the end of OUTPUT file:

                             MEAN VALUE       +/-  RMS DEVIATION
                                    <x>     [<x^2>-<x>^2]**(1/2)
ELECTRON KINETIC ENERGY    0.384138E-01             0.719529E-02
IONIC TEMPERATURE                262.59                    15.36
DENSITY FUNCTIONAL ENERGY  -1095.541142             0.180823E-01
CLASSICAL ENERGY           -1095.305409             0.734628E-02
CONSERVED ENERGY           -1095.266996             0.157838E-03
NOSE ENERGY ELECTRONS          0.000000              0.00000
NOSE ENERGY IONS               0.000000              0.00000
CONSTRAINTS ENERGY             0.000000              0.00000
ION DISPLACEMENT           0.643998                 0.320104

These can be very useful, especially for the temperature (and its standard deviation) and for the conserved energy. Probably for EKINC a plot would make much clearer whether everything is going well.

2. Charged system:

specify the overall charge in the unit cell. It should be in the &SYSTEM section, something like:

CHARGE

-1

It will find out where to put the extra charge itself, depending on the ordering and localisation of the Kohn-Sham orbitals computed during the wavefunction optimisation.

3. TRAJECTORY: The trajectory is in a file called TRAJECTORY. Everything is in atomic units. Section: &CPMD (With the additional keyword XYZ the trajectory is also written in xyz-format on the file TRAJEC.xyz. But for some unknown reason, although I copied TRAJEC.xyz to my work directory on hpc in the run script, I cannot find TRAJEC.xyz).

Solution: new little code: traj2xyz.x

./traj2xyz.x < TRAJECTORY > traj.xyz

Use VMD to visualise traj.xyz and then we can see the movie.

4. ENERGIES: Plot temperature and total energy, and see when they are oscillating around an average value, without drift. This information is given in the outoput file and in a file called ENERGIES.

5. Plot total energy vs Ecut to check convergence of Ecut.

6. EMASS: fictitious electron mass. Set to 500 a.u. atm, which is the same as that in Tricia's previous input.

In Kuo04, P12997 IV. Conclusions:

"As shown by Grossman et al.12 and in this work, values of mu=800 au (for H2O) or mu= 1100 au (for D2O) are inappropriate for present day simulations of liquid water that are accessing a length of more than just a few picoseconds. However, a value of mu=400 au is sufficiently small to allow the extension of microcanonical CPMD simulations for H2O at ambient conditions to many tens of picoseconds."

=> So we decide to go leave my EMASS at 500 a.u. for the future simulations.

7. Run 3ps, and then 3ps. In the latter 3ps, analyse 20 structures including Kohn-Sham orbitals, Wannier centres and work out how to print Wannier Matrix.

7-July-09: w63Cl- 10000 steps running on hpc, estimate ~85 hrs on 4 cpus.

timestep = 5 a.u.

1 a.u. = 0.0241888428 fs

10000 steps = 50000 a.u. = 1.21 ps