Torsion and Bending PBC in MD++

Torsion and Bending PBC in MD++

Wei Cai

This tutorial describes how Torsion/Bending Periodic Boundary Conditions (t-PBC and b-PBC) are implemented in MD++. We provide several MD++ script files and describe how to use them in detail. The theoretical background is published in

Torsion and Bending Periodic Boundary Conditions for Modelling the Intrinsic Strength of Nanowires, Journal of Mechanics and Physics of Solids, 56, 3242, (2008). (PDF)

Compile executable file

First, you need to install MD++ on your computer by following the instructions on MD++ Manuals. Suppose the root directory of your MD++ installation is ~/Codes/MD++. For convenience, let us define an environment variable,

export MDPP=~/Codes/MD++

You will also need the ${MDPP}/cookies directory (not included in the standard MD++ distribution), where additional source files such as md_torsion.cpp are located.

To turn on the t-PBC and b-PBC utility, modify src/Makefile so that

 TOB = yes

In this tutorial, we will use gold nanowire (modelled by the EAM potential) as an example. We can compile the executable file (e.g. on mc-cc.stanford.edu) by

 cd ${MDPP}
 make eam build=R SYS=mc-cc

This will create executable eam_mc-cc in the ${MDPP}/bin directory.

Make handle atoms

As the example, let us run the following test case,

cd ${MDPP}
bin/eam_mc-cc scripts/work/au/au-nw-torsion.tcl 0

The important parts of this script file (activated in this run) is reproduced below.

  # Create the nanowire
  ...
  MD++ allocmultiple = 2
  makeperfcrystal $ND $NL $lattconst
  make_nanowire $R
  ...
  openwindow
  ...
  # Create handles for applying t-PBC
  set dphi 0.02
  set dz   0.1
  set phi  0
  MD++ torsionsim = 1 torquespec = \[ $dphi $dz $phi \] maketorquehandle

  #(note: we will rename torquespec to twistspec in the future)

  # Conjugate gradient relaxation
  relax_fixbox
  MD++ finalcnfile = "nw-relaxed.cn" writeall = 1 writecn

  # Molecular Dynamics simulation
  setup_md
  MD++ T_OBJ = 1000 totalsteps = 1000 initvelocity run
  MD++ finalcnfile = "handle-annealed.cn" writeall = 1 writecn

  # Conjugate gradient relaxation (again)
  relax_fixbox

  MD++ finalcnfile = "handle-annealed-relaxed.cn" writeall = 1 writecn

As MD++ executes this script, the first thing it does is to create a perfect crystal of FCC gold, and cut it into a cylinder with z-axis along the [111] direction. In this test case, the radius of the nanowire is R = 25 Angstrom. Notice that we have put allocmultiple = 2 before the makeperfcrystal command. This instructs MD++ to allocate twice as many atoms as necessary; the extra memory space will be used later to store the "handle" atoms to apply t-PBC.

The next important command is maketorquehandle. The input parameter to this command is specified in torquespec. Only the second variable (dz) is actually used. The first variable (dphi) will be used to add addition twist angle to the nanowire, after the handle has been created. The third variable (phi) specifies the total twist angle as one goes along the wire by distance Lz. Since we have not applied any twist yet, it is a good idea to set it to zero for now.

It identifies atoms whose scaled coordinate satisfies sz < -0.5 + dz, copies them and paste them at a new location sz_new = sz + 1. It identifies atoms whose scaled coordinate satisfies sz > 0.5 - dz, copies them and paste them at a new location sz_new = sz - 1. This will increase the total number of atoms, _NP. The original number of atoms is stored in _NP0. The number of new "handle" atoms is _NIMAGES = _NP - _NP0.

The new "handle" atoms have index i=_NP0.._NP-1, for which fixed[i] = 1. The corresponding entries in the images[] array stores the index of their "master" atoms. The "handle" atoms are slaves to their "masters" in that whenever the potential function is called, their positions are obtained by copying, translate (and rotate), and paste operations defined above.

To avoid "handle" atoms interacting through the PBC, the length of the simulation box, Lz (i.e. H[2][2]) is enlarged, Lz := Lz * (1+4*dz). To keep the real coordinates of all atoms, the scaled coordinates shrink to sz := sz / (1+4*dz).

There are some flexibility in the choice of dz but it must satisfy two constraints: (1) dz < 0.5 (because of allocmultiple = 2), (2) Lz * dz > Rc (the cut-off radius of the interatomic potential).

After maketorquehandle, the wire is still subjected to conventional PBC along z-axis. However, PBC is now applied by the copy-translate-paste operation which was not necessary before.

Finally, MD++ performs conjugate gradient relaxation on the nanowire (subjected to PBC), runs Molecular Dynamics (MD) simulation for a short time (to anneal it), and relaxes it again. Every step of the relaxation or MD simulation calls the call_potential function. The important parts of the latter function is reproduced below.

    if (_NIMAGES>0)
    {
       if (_TORSIONSIM) copytorqueatoms();
       if (_BENDSIM)    copybendatoms();
    }
    ...
         potential();
    ...
    if (_NIMAGES>0)
    {
        for(i=0;i<_NIMAGES;i++)
        {
             _EPOT -= _EPOT_IND[_NP0+i];
             _VIRIAL = _VIRIAL - _VIRIAL_IND[_NP0+i];
             _TORQUE -= _TORQUE_IND[_NP0+i];
             _BENDMOMENT -= _BENDMOMENT_IND[_NP0+i];
        }
    }

Before calling the "real" potential function, MD++ calls the copytorqueatoms function (better be called copyatoms_torsion). (This is because _NIMAGES is positive by now, and _TORSIONSIM has been set to 1 in the script file.) After the potential function call, MD++ subtracts the contributions "handle" atoms made to the potential energy, Virial stress, Torque and Bending Moment to avoid double-counting.

The self-consistency of the implementation can be verified by the success of the conjugate gradient relaxation, which requires forces on every atom to be consistent with the total potential energy (_EPOT).