Isolated Nitrogen molecule (NVE)

This example requires minimal computational effort and should be reproducible on just about any modern desktop system.

The DFT calculation

As this is a simulation of an isolated molecule, we will use a single k-point ($\Gamma$ point) in the electronic structure calculation. The cell file that we will use is given below.

!Example isolated and strained N2 molecule

%BLOCK LATTICE_ABC ! In Angstroms
  6.2  6.2  6.2
  90.0 90.0 90.0
%ENDBLOCK LATTICE_ABC


%block POSITIONS_FRAC
N   0.406932   0.500000   0.500000
N   0.593068   0.500000   0.500000
%endblock POSITIONS_FRAC

%BLOCK SPECIES_POT
N N_00.usp
%ENDBLOCK SPECIES_POT

%BLOCK KPOINTS_LIST
0.0 0.0 0.0  1
%ENDBLOCK KPOINTS_LIST

The energy is converged w.r.t. basis set size to within 0.004 eV at a cut-off energy of 320 eV. Convergence w.r.t. the size of the super-cell indicates similar accuracy with a cell dimension of 6.2 Å. The ionic positions have then been contracted to the resulting equilibrium positions using the BFGS geometry optimiser.

Figure 4.1: Charge density of the N$_{2}$ molecule converged w.r.t. basis set size and cell volume.

Finding a Suitable Time-step

We now identify a suitable time-step for the md calculation. The ultimate aim is to simulate the longest possible time with the smallest amount of computational effort. The choice of time-step is therefore a compromise between a large value which will require less MD steps for a given simulated time, and a small value for which better re-use of the previous wave-function will be made allowing each time-step to be calculated faster. We must also make sure that the time-step is sufficiently small that the dynamics correctly conserve the total (Hamiltonian) energy of the system.

As a first guess, we will examine the performance with a time-step of 1 fs using the following parameters file. Note that the N$_{2}$ molecule is at its equilibrium bond length and so will only begin to oscillate if we assign initial velocities. We will set set md_temperature to 300 K to give the ions some initial velocity.

comment          = example NVE dynamics input file
task             = molecular dynamics

cut_off_energy   = 320 eV
fix_occupancy    = true
xc_functional    = LDA

md_ensemble      = NVE
md_delta_t       = 1.0 fs
md_num_iter      = 1000
md_temperature   = 300 K

rand_seed = 1

We have used a fixed random number seed to ensure results will match those given here. Also add

fix_com = .true.

to the cell file. This will keep the N$_{2}$ molecule in the centre of our cell. After running with a time-step of 1 fs, repeat the calculation using time-steps of 0.5 fs and 2 fs. It is useful to compare the conservation of the Hamiltonian for each of the three runs graphically as in figure 4.2.

Figure 4.2: Hamiltonian and potential energy for the N$_{2}$ molecule at 300 K using three different time-steps.

Examination of this figure indicates a small but noticeable drift from the initial value of the Hamiltonian in the 2 fs case. No advantage in the cpu time per time-step is gained by decreasing to 0.5 fs, and therefore we conclude that the optimal time-step is around our initial guess of 1 fs.

Eliminating Rotations

In the above MD simulations, velocities were assigned randomly to match an initial temperature, including components not along the direction defined by the N-N bond. This leads to rotation of the molecule as is easily confirmed by visualisation. Two methods for eliminating these rotations are given below.

Specifying Initial Conditions

We can easily eliminate rotations by either specifying the initial velocities manually to only include components in a single direction, or alternatively by specifying zero initial temperature (and hence zero velocity) and applying a stretch to the N$_{2}$ bond, i.e.

md_temperature = 0 K

in the parameters file, and in the cell file make the following alteration.

%block POSITIONS_FRAC
N   0.403932   0.500000   0.500000
N   0.596068   0.500000   0.500000
%endblock POSITIONS_FRAC

This stretches the bond by 0.372 Åin the x-direction only. If we now run this simulation, we should see no rotations.

Using Constraints

A second way of ensuring that the atoms only move in the x-direction, is to specify constraints for the y and z co-ordinates of each atom. This requires specification of four linear constraints in the cell file. We will also keep the centre of mass constrained as before. Note that we cannot allow the fix_com keyword to do this as it will generate constraints for all three components of the centre of mass, including those we have already constrained by specifying fixed y and z. We therefore specify the centre of mass constraint as manual constraint number 5 as below.

%BLOCK IONIC_CONSTRAINTS
1 N 1  0 1 0
2 N 2  0 1 0
3 N 1  0 0 1
4 N 2  0 0 1
5 N 1  1 0 0 
5 N 2  1 0 0
%ENDBLOCK IONIC_CONSTRAINTS

fix_com = .false.

This method is preferred over the initial conditions method, as it avoids having to determine the specific initial conditions that would cause the molecule to oscillate at the desired temperature.

Results

After running the simulation for 1000 time-steps, we are in a position to extract useful data from the output. The mechanics of this will be discussed later and can involve either converting the .md file into a format suitable for an analysis program, or using UNIX utilities like grep and awk to extract specific data.

As an example we will examine some properties of the N$_{2}$ bond. A plot of the bond length as a function of time is shown in figure 4.3 along with its Fourier transform.

Figure 4.3: Length of the N$_{2}$ bond as a function of time, and the corresponding Fourier Transform.

From these plots it is trivial to extract that the average bond length is 1.294 Å and that the bond frequency is 0.0725 fs$^{-1}$.If we were interested in taking this example further, we could examine the temperature dependence of the bond length and frequency, and determine the bond dissociation temperature.