How to use moltemplate to construct crystals, such as SiC

[S1mplecool]

Hello, Andrew.
I'm a student who just started using moltemplate. I found that the examples in official website of moltemplate do not include crystal construction, such as hexagonal crystal of SiC. This confuses me, and I can't build it after looking for many methods. Can you give me a simple example for me to learn? Please reply me, thank you very much.

[jewettaij]

Hi S1mplecool

Here are a couple of examples showing how to build simple crystals with moltemplate:

The first example (FCC aluminum) is probably the most simple and useful:

1. Aluminum (FCC)
   a) File that defines the unit cell: al_cell.lt
   b) File that defines a small crystal (10x10x10): system.lt
   c) Pictures and description: README.md
   d) Command to run moltemplate.sh:

moltemplate.sh system.lt    # (add "-vmd" if VMD is installed)

In this example, the forces between atoms are determined by a many-body force field ("pair_style eam/alloy"). No explicit bonds between atoms are used.

Here is a more fun version of this example

2. Water Ice (hexagonal):
   a) File that defines a water molecule: spce.lt
   b) File that defines the unit cell of 8 water molecules: spce_ice_rect8.lt
   c) File that defines a small crystal (3x2x2): system.lt
   d) Pictures and description: README.md
   In this example, the forces between atoms in the same molecule are determined by bonds (and RATTLE constraints), and the forces between atoms in different molecules are determined by Lennard-Jones and electrostatic forces.

3. Graphene (2-D hexagonal crystal):
   a) File that defines the 2-D unit cell: graphene.lt
   b) Instructions how to build a 2-D graphene crystal are here
   c) Pictures of graphene and carbon nanotubes built this way: README.md

In this example (and the next example: graphite), there are no bonds between atoms. It is possible to use many-body force fields to simulate the motion of these atoms. However at low temperatures (eg. 300K), I would not expect the atoms to move enough to deform or break the graphene. For this reason, I decided to make the graphene rigid. (To do that, I used "fix rigid".) The method used to do this is different depending on whether you are running your simulations at constant volume or constant pressure. (See the run.in.nvt and run.in.npt files for this example.)

4. Graphite (3-D hexagonal crystal):
   a) File that defines the 3-D unit cell: graphite.lt
   b) Instructions how to build a 3-D graphite crystal are here

Note 1: If you want to put impurity substitutions at random places in the crystal, this is possible, but it is more difficult. Probably the easiest way to do that in this case would be to write a short for-loop (in python or even BASH, for example) that generates a list of LAMMPS "set" commands to change the type of random atoms in the crystal. This might look something like this:

:
set atom 1486 type 1
set atom 1487 type 2
set atom 1488 type 2
set atom 1489 type 1
set atom 1490 type 1
set atom 1491 type 1
set atom 1492 type 1
set atom 1493 type 2
set atom 1494 type 1
set atom 1495 type 1
set atom 1496 type 1
set atom 1497 type 1
:

Save this text in a file (eg "random_types.in"), and read them into your main LAMMPS input script (eg "run.in") using the "include" command in LAMMPS (eg. "include random_types.in"). I'm not sure if the positions of the impurities should be truly random, however. (They might cluster together.)

Note 2: If you want to have explicit bonds between atoms in different unit cells, this is difficult to do using moltemplate (or VMD/topotools, or just about any other software). You would have to add these bonds manually. (There is a way to do this automatically, but I did not want to spend the time to implement it. It gets a little bit tricky when you have periodic boundary conditions. I can share this with you if you want.)

I hope at least one of these examples is useful to you.