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. You will need to copy script file au-nw-torsion.tcl to the appropriate directory.
cd ${MDPP}
bin/eam_mc-cc scripts/work/au/au-nw-torsion.tcl 0 0
The first "0" in the command line tells MD++ which job to do. The second "0" tells MD++ which directory to save the output files. 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 = "handle-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 1
Here are the energies and stress of the nanowire at various stages of the two simulation. Notice that the stress value reported here are not the true stress in the nanowire. This is because we used the virial formula that requires division by a volume, which is taken to be the volume of the simulation box (much bigger than the nanowire volume). (Original crystal has 77064 atoms. After making the cylinder 15637 atoms remain. The box is then further enlarged by a factor of 2 in both x and y directions. So the stress need to be multiplied by a factor of 77064/15637 * 2 * 2 = 19.7132) Furthermore, because the box is further elongated by (1+4*dz) when creating the handle, we need to multiply this value to the stress reported by t-PBC simulations before the comparison. The values reported here are already multiplied by 19.7132 (for both t-PBC and PBC) and also by (1+4*dz) for t-PBC only. This is done by computing the "vacuumratio" variable in the tcl script.
| t-PBC with phi=0 (energy) | t-PBC with phi=0 (tensile stress) | conventional PBC (energy) | conventional PBC (tensile stress) | |
|---|---|---|---|---|
| as created nanowire | -5.92058e4 (eV) | 4534 (MPa) | -5.92058e4 (eV) | 4534 (MPa) |
| after 1st relax | -5.96779e4 (eV) | 1231 (MPa) | -5.96779e4 (eV) | 1218 (MPa) |
| after MD (1000K) | -5.70067e4 (eV) | 2203 (MPa) | -5.70069e4 (eV) | 2218 (MPa) |
| after 2nd relax | -5.96806e4 (eV) | 392.3 (MPa) | -5.96808e4 (eV) | 396.2 (MPa) |
When the nanowire is first created, the energy and stress in t-PBC and PBC are identical (as they should). After relaxation and annealing, the energy and stress values become slightly different (because small round off error can drive the system to different final states). It is gratifying to see that their difference remain small even after the 2nd relaxation. This indicates that t-PBC (at least when phi=0) is implemented correctly.
Relax axial stress before applying torsion
Gold nanowires are known to shrink its length from the perfect lattice values due to the surface stress. Therefore, it is desirable to relax its length before subsequent torsion simulations. At zero temperature, this can be done by conjugate gradient relaxation while allowing the box length to change. At finite temperature (to allow thermal expansion), this can be done by MD simulations under the NPT ensemble. In either case, whether the t-PBC handle is made before or after the length relaxation/equilibration should not change the result (for self-consistency). In this section, we present numerical results to show that this is the case. For simplicity, we make sure all simulations use the same "handle" atoms (as in the above test case) by saving this information into a file (savehandle) and load it when applying t-PBC (readhandle).
The following table shows the energy and length when the nanowire length is allowed to relax/equilibrate under t-PBC (phi=0) and conventional PBC. The test cases are performed by the following commands.
cd ${MDPP}
bin/eam_mc-cc scripts/work/au/au-nw-torsion.tcl 2 2
bin/eam_mc-cc scripts/work/au/au-nw-torsion.tcl 3 3
| t-PBC with phi=0 (energy) | t-PBC with phi=0 (length) | conventional PBC (energy) | conventional PBC (length) | |
|---|---|---|---|---|
| as created nanowire | -5.92058428e4 (eV) | 13.4269 (nm) | -5.92058428e4 (eV) | 13.4269 (nm) |
| after 1st relax | -5.96896563e4 (eV) | 13.2722 (nm) | -5.96896497e4 (eV) | 13.2759 (nm) |
| after MD (300K) | -5.9058278e4 (eV) | 13.3188 (nm) | -5.9059614e4 (eV) | 13.3185 (nm) |
| average during 2nd half of MD | -5.9061118e4 (eV) | 13.3196 (nm) | -5.9061491e4 (eV) | 13.3200 (nm) |
| after 2nd relax | -5.9691107e4 (eV) | 13.2760 (nm) | -5.9691581e4 (eV) | 13.2755 (nm) |
By averaging the second half of the MD simulation, we find that the equilibrium length of the nanowire at T = 300 K is 13.320 nm (corresponding to a thermal expansion coefficient of 1e-5 /K). Therefore, for subsequent torsion simulations we equilibrate the nanowire by NVT simulations with the nanowire length fixed at 13.32 nm. The initial condition for this simulation is the final configuration (and velocity) of the previous NPT simulation (except that the length is scaled to 13.32 nm). (If we use a relaxed configuration as the initial condition, artificial oscillation of the axial stress will persist for a long time.) Again this can be done under either t-PBC (phi = 0) or PBC, using the following commands.
cd ${MDPP}
bin/eam_mc-cc scripts/work/au/au-nw-torsion.tcl 4 4
bin/eam_mc-cc scripts/work/au/au-nw-torsion.tcl 5 5
Here are a comparison of the final and averaged energy and stress of the NVT MD simulation. The stress in t-PBC simulation are multiplied by (1 + 4*dz) to compare with that in PBC simulation.
| t-PBC with phi=0 (energy) | t-PBC with phi=0 (stress) | conventional PBC (energy) | conventional PBC (stress) | |
|---|---|---|---|---|
| NPT equilibrated nanowire with Lz scaled to 13.32 nm | -5.9058e4 (eV) | 8.0 (MPa) | -5.9060e4 (eV) | 0.16 (MPa) |
| after MD (300K) | -5.?????e4 (eV) | ???? (MPa) | -5.9066e4 (eV) | 60 (MPa) |
| average during MD | -5.?????e4 (eV) | ???? (MPa) | -5.9062e4 (eV) | 5.6 (MPa) |
| standard deviation during MD | -5.?????e4 (eV) | ???? (MPa) | 3.8 (eV) | 31 (MPa) |
It is good to see that the gold nanowire can be equilibrated to (close-to) zero stress at 300K using both t-PBC (phi=0) and conventional PBC. Artificial stress oscillation (with amplitude of 1GPa) no longer exists in the above NVT simulations.
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.