This example requires minimal computational effort and should be reproducible on just about any modern desktop system.
As this is a simulation of an isolated molecule, we will use a single k-point ( 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.
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 molecule is at its equilibrium bond length and so
will only begin to oscillate if we assign initial velocities. We will 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 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.
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.
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.
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 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.
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.
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 bond. A plot of the bond length as a function of time is shown in figure 4.3 along with its 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.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.