258 Molecular dynamics simulations
simply by shifting the bits ofNXandNYone position to the right (least significant
direction) and it can therefore easily be checked in the program whether the parents
of the squares under consideration are neighbours or not.
The calculation of the multipole moments in each box(Eq. (8.140))is best done
in a loop over the particles, recording its contribution to all the multipole coefficients
of the to which square it belongs. Also, the calculation of the interactions
(Eq. (8.139))can be done in a loop over the particles, by executing for each particle
a loop over the interaction list of the square to which it belongs.
Proceeding this way, it is not necessary to keep for each square a list of the
particles belonging to it. However, at the finest level, the interactions between
particles within the same square and between particles in neighbouring boxes must
be calculated directly so only for the last step do we need such a list for each square.
If you want to economise on memory, you might create a linked list for each square
containing the indices of the particles in it, but for a test you may use static arrays.
Compare the results for the tree code with those of a direct calculation, varying
the number of terms in the multipole expansion.
8.12 [C] In this problem we consider a simulation of a methane molecule using the
Langevin approach. Methane consists of a carbon atom sitting at the centre of a
tetrahedron whose vertices are occupied by four H atoms. The C–H bond has a
preferred interatomic distance of 2.104a 0. The stretch-potential associated with the
bond length varies as
Vstretch=^12 κ(l−l 0 )^2 ; l 0 =2.104a 0.
The force constantκhas the valueκ=0.30 (in atomic units). This force acts on
both the carbon and the hydrogen atoms and is directed along the C–H bond.
The preferred H–C–H angle is 109◦and the potential for this bending angle is
Vbend=−λcos(φ−φ 0 )^2 ; φ 0 = 109 ◦,
with a force constantλ=0.74. This force lies in the H–C–H plane, and acts on the
two H atoms and on the C atom. The forces on the H-atoms are perpendicular to the
C–H bonds, and the bending force on the C atom is directed along the bisecting line
of the H–C–H angle.
These two ‘force fields’, bending and stretching, specify the force on each of the
atoms. To find the forces, given the positionrCof the carbon atom and the four
positionsrHof the hydrogen atoms, you calculate first the forces on the hydrogen
atoms only. The stretch forces can easily be found by calculating the vector
rCH=rH−rC. The bending force is slightly more difficult. Denoting the two
hydrogen atoms of a H–C–H chain as H1 and H2, calulaterCH1andrCH1. Then
calculate the dot product between these two vectors. From this, the cosine of the
bending angle can be found. Moreover, the direction of the force can be found from
the cross-product ofrCH1andrCH1: the bending force on H1 is then perpendicular to
this cross productandto the vectorrCH1, and similarly for H2. Knowing the forces
on the hydrogen atoms, you can calculate their sum. The force on the carbon atom is