Torsion and Bending PBC in MD++

Torsion and Bending PBC in MD++

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).

The relaxed and annealed configurations are saved into .cn files. Notice that writeall = 1 is set before writecn. This instructs MD++ to save not only the scaled coordinates of each atom, but also their values in the fixed[] and image[] arrays.

Numerical verification

The nanowire in this test case is still subjected to PBC (although the implementation is now different). Hence the result should be consistent with another simulation where conventional PBC is applied, which can be tested by running,

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

Here are the energies of the nanowire at various stages of the two simulation.

	t-PBC with phi=0	conventional PBC
as created nanowire	?? (eV)	?? (eV)
after 1st relax	?? (eV)	?? (eV)
after MD	?? (eV)	?? (eV)
after 2nd relax	?? (eV)	?? (eV)

Questions and answers

In this section, we address a set of questions concerning the (pecularities of) the existing implementation of t-PBC in MD++. After the discussion, it becomes apparent that the implementation can benefit from a re-design. This section will be modified after the re-implementation is complete.

Q: Do we need to specify dz in subsequent simulations after handle has been made? A: Yes. After handle has been made, the simulation box is enlarged to make room for the handle (plus the empty space). The actual length of the nanowire Lz is only a fraction of the box length H_zz, i.e. Lz = H_zz / (1 + 4 * dz). When we copy-translate-paste an atom, the scaled coordinate should change by sz_new = sz + 1/(1 + 4 * dz). Therefore, it makes sense to make dz a permanent member variable in MDFrame (instead of being part of torquespec).

Q: Do we need to specify NP0 and NIMAGES in subsequent simulations? A: Yes, we need to do so in the current implementation. However, we should be able to eliminate this need in a new implementation, where the image[] entry of "real" atoms will be set to -1. In this way, MD++ should be able to calculate how many atoms are "real" atoms (NP0) and NIMAGES = NP - NP0. For this to work, writeall = 1 has to be set before writecn so that the image[] value will be saved into .cn file.

Q: Does t-PBC work with NPT MD simulation? A: In principle yes. We have to test it thoroughly. It seems a good idea to create utilities to write and read handle information (i.e. images[] array) so that all simulations with the same nanowire can be performed with the same handle, regardless of whether the handle is applied before or after NPT simulations. We also need to debug why in some simulations NW failure always occur at NW-handle interfaces.