<jupyter_start><jupyter_text>Building crystal and nanoparticle atomistic models  Table of Content 1. [General setups](setups)2. [Periodic structure builders](2.)   2.1. [Reading xyz inputs](2.1.)    2.2. [Creating supercells](2.2.)    2.3. [Atom-by-atom construction](2.3.)    2.4. [All atoms together - using crystal coordinates](2.4.)    2.5. [All atoms together - using Cartesian coordinates](2.5.)  3. [Determining connectivities](3.)   3.1. [Non-periodic case](3.1.)    3.2. [Periodic case (Warning: a bit buggy for now)](3.2.)4. [Quantum Dot building](4.)   4.1. [With capping](4.1.)    4.2. [Without capping](4.2.)    4.3. [Streamlined building of QDs](4.3.)    A. Learning objectives- to build extended periodic systems and quantum dots starting from xyz files- to build extended periodic systems and quantum dots starting from scratch- to create QDs of needed size- to replace surface atoms of of QE with other atoms- to detect connectivities in complex structures and print that information B. Use cases- building beriodic structures and supercells- building quantum dots- molecular structure visualization C. Functions- `libra_py`    - `autoconnect`    - [`autoconnect`](autoconnect-1)     - [`find_undercoordinated_atoms`](find_undercoordinated_atoms-1)  - `build`    - [`read_xyz`](read_xyz-1) | [also here](read_xyz-2)    - [`make_xyz`](make_xyz-1) | [also here](make_xyz-2)    - [`read_xyz_crystal`](read_xyz_crystal-1)     - [`generate_replicas_xyz`](generate_replicas_xyz-1)     - [`generate_replicas_xyz2`](generate_replicas_xyz2-1) | [also here](generate_replicas_xyz2-2)    - [`crop_sphere_xyz`](crop_sphere_xyz-1)     - [`crop_sphere_xyz2`](crop_sphere_xyz2-1) | [also here](crop_sphere_xyz2-2)    - [`crop_sphere_xyz3`](crop_sphere_xyz3-1) | [also here](crop_sphere_xyz3-2)    - [`add_atom_to_system`](add_atom_to_system-1)    - [`add_atoms_to_system`](add_atoms_to_system-1) | [also here](add_atoms_to_system-2)  | [also here](add_atoms_to_system-3)     D. Classes and class members- `liblibra::libchemobjects::libchemsys`  - `System`    - [`System`](System-1)    - [`get_xyz`](get_xyz-1) | [also here](get_xyz-2)    - [`init_box`](init_box-1)    - [`print_ent`](print_ent-1) | [also here](print_ent-2)      1. General setups[Back to TOC](TOC)Let's import all the needed libraries, including the `py3Dmol` library for visualizing the structures we generate.<jupyter_code>import math
import liblibra_core
from liblibra_core import *
from libra_py import build
from libra_py import units
from libra_py import autoconnect

import py3Dmol
import matplotlib.pyplot as plt
#from matplotlib.pyplot import figure
%matplotlib inline

plt.rc('axes', titlesize=24)      # fontsize of the axes title
plt.rc('axes', labelsize=20)      # fontsize of the x and y labels
plt.rc('legend', fontsize=20)     # legend fontsize
plt.rc('xtick', labelsize=16)    # fontsize of the tick labels
plt.rc('ytick', labelsize=16)    # fontsize of the tick labels

plt.rc('figure.subplot', left=0.2)
plt.rc('figure.subplot', right=0.95)
plt.rc('figure.subplot', bottom=0.13)
plt.rc('figure.subplot', top=0.88)

colors = {}
colors.update({"11": "#8b1a0e"})  # red       
colors.update({"12": "#FF4500"})  # orangered 
colors.update({"13": "#B22222"})  # firebrick 
colors.update({"14": "#DC143C"})  # crimson   

colors.update({"21": "#5e9c36"})  # green
colors.update({"22": "#006400"})  # darkgreen  
colors.update({"23": "#228B22"})  # forestgreen
colors.update({"24": "#808000"})  # olive      

colors.update({"31": "#8A2BE2"})  # blueviolet
colors.update({"32": "#00008B"})  # darkblue  

colors.update({"41": "#2F4F4F"})  # darkslategray<jupyter_output>/home/alexey/Conda/Miniconda3/envs/libra/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: to-Python converter for std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > already registered; second conversion method ignored.
  return f(*args, **kwds)
/home/alexey/Conda/Miniconda3/envs/libra/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: to-Python converter for boost::python::detail::container_element<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, unsigned long, boost::python::detail::final_vector_derived_policies<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, false> > already registered; second conversion method ignored.
  return f(*args, **kwds)
/home/alexey/Conda/Miniconda3/envs/libra/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: to-Python converter for std::vector<std::v[...]<jupyter_text>2. Periodic structure builders[Back to TOC](TOC) 2.1. Reading xyz inputs[Back to TOC](TOC)Define a unit cell parameters - periodic translations of the minimal unit cell (provided as an input file)<jupyter_code>a = VECTOR(5.833, 0.0, 0.0) * units.Angst
b = VECTOR(0.0, 5.833, 0.0) * units.Angst
c = VECTOR(0.0, 0.0, 5.833) * units.Angst<jupyter_output><empty_output><jupyter_text>Read the minimal unit cell information from the xyz file using:<jupyter_code>help(build.read_xyz)
L0, R0 = build.read_xyz("CdSe.xyz")
print(L0, R0)<jupyter_output>['Cd', 'Cd', 'Cd', 'Cd', 'Se', 'Se', 'Se', 'Se'] [<liblibra_core.VECTOR object at 0x7f5a57e958b8>, <liblibra_core.VECTOR object at 0x7f5a57e95648>, <liblibra_core.VECTOR object at 0x7f5a57e95e68>, <liblibra_core.VECTOR object at 0x7f5a57e95920>, <liblibra_core.VECTOR object at 0x7f5a57e95988>, <liblibra_core.VECTOR object at 0x7f5a57e95850>, <liblibra_core.VECTOR object at 0x7f5a57e95a58>, <liblibra_core.VECTOR object at 0x7f5a57e95ac0>]<jupyter_text>Now, using atomic labels `L0` and coordinates `R0`, we can create the string file that represents the xyz file format, using:<jupyter_code>help(build.make_xyz)<jupyter_output>Help on function make_xyz in module libra_py.build:

make_xyz(L, R, inp_units=0, out_units=1)
    Convert atomic labels and coordinates to a string formatted according to xyz
    
    Args:
        L ( list of N strings ): the labels of N atoms in the system
        R ( list of N VECTORs ): the coordinates of N atoms in the system
        inp_units ( int ): defines the units of variables stored in R:
    
            * 0 - Bohr ( default )
            * 1 - Angstrom 
    
        out_units ( int ): defines the units of the coordinates written in xyz file:
    
            * 0 - Bohr 
            * 1 - Angstrom ( default )
    
    Returns: 
        string: a string representing the xyz file of the system<jupyter_text>In principle, for the purpose of visualization, we could have simply read the xyz file as a string and used the string in the py3Dmol viewer. However, the current way also gives us the digital representation of atomic information (atom element labels and coordinates), which we will be using in the following steps.<jupyter_code>xyz = build.make_xyz(L0, R0)
view = py3Dmol.view(width=800,height=400)  # linked=False,viewergrid=(3,2)
view.setBackgroundColor('0xeeeeee')
view.zoomTo()
view.addModel(xyz,'xyz',{'vibrate': {'frames':10,'amplitude':1.0}})
view.setStyle({'sphere':{'colorscheme':'Jmol'}})
view.show()

print(xyz)<jupyter_output><empty_output><jupyter_text>The `make_xyz` function could be used as a first step before storing the molecular structure information as an xyz file - just produce the string and then save it to the file.  2.2. Creating supercells[Back to TOC](TOC)Now, lets create a super-cell by replicating the minimal cell in 3 directionsWe need large enough super-cell so that the QDs we are going to create fits inside itso, replicate 6 times in each directionThis will generate the labels and coordinates of the new super-system<jupyter_code>help(build.generate_replicas_xyz2)
Nx, Ny, Nz = 6, 6, 6 
L1, R1 = build.generate_replicas_xyz2(L0, R0, a, b, c, Nx, Ny, Nz)

#print( L1, len(R1))

xyz = build.make_xyz(L1, R1)
view = py3Dmol.view(width=800,height=400)  # linked=False,viewergrid=(3,2)
view.setBackgroundColor('0xeeeeee')
view.zoomTo()
view.addModel(xyz,'xyz',{'vibrate': {'frames':10,'amplitude':1.0}})
view.setStyle({'sphere':{'colorscheme':'Jmol'}})
view.show()<jupyter_output><empty_output><jupyter_text>The instructions in the commented line would have resulted in a lot of text output, since the newly generated lists L1 and R1 are quite large.   2.3. Atom-by-atom construction[Back to TOC](TOC)We can also construct systems by adding atoms to an object of the class `System` defined in the core of Libra package. The `System` class represent a general atomistic system and has mutlitude of methods and data members.Here, we are using the `add_atom_to_system` function that does such an addition:<jupyter_code>help(build.add_atom_to_system)<jupyter_output>Help on function add_atom_to_system in module libra_py.build:

add_atom_to_system(syst, coords, MaxCoords, Nx, Ny, Nz, a, b, c, shift, elt, mass, scl, max_coord, rnd, inp_units=0)
    This function adds an atom and its periordic replicas to a (chemical) System object. 
    The momenta of the atoms are also initialized and are drawn from the normal distribution along
    each of the directions (Cartesian x, y, and z), no center of mass momentum subtraction is made however.
    This is good for constructing lattices.
    
    Args:
        syst ( System object ): the chemical system to which we are going to add the atoms
        coords ( list of VECTORs ): coordinates of all atoms in a system
        MaxCoords (list of ints): Maximal coordination number of each atom
        Nx ( int ): how many times to replicate the atom along a direction
        Ny ( int ): how many times to replicate the atom along b direction
        Nz ( int ): how many times to replicate the atom along c direction
[...]<jupyter_text>We will add only 2 basis atoms (in the crystallochemical terms),and replicate each of them well in all 3 directions. Below, we identify directions of crystal periodicity vectors, which determine the shape of the resulting structure.<jupyter_code>syst1 = System()
rnd = Random()

# Unit cell parameters
a = VECTOR(5.0, 5.0, 0.0) * units.Angst
b = VECTOR(0.0, 5.0, 0.0) * units.Angst
c = VECTOR(0.0, 0.0, 5.0) * units.Angst

# Create a Chemical system
Nx, Ny, Nz = 5, 5, 5  # 

R, MaxCoord = [], []
build.add_atom_to_system(syst1, R, MaxCoord, Nx,Ny,Nz, a, b, c, VECTOR(0.0, 0.0, 0.0) * units.Angst, "Na", 2000.0, VECTOR(1.0, 1.0, 1.0), 6, rnd)
build.add_atom_to_system(syst1, R, MaxCoord, Nx,Ny,Nz, a, b, c, VECTOR(2.5, 2.5, 2.5) * units.Angst, "Cl", 2000.0, VECTOR(1.0, 1.0, 1.0), 6, rnd)<jupyter_output><empty_output><jupyter_text>As a result, we obtained a cupercell of the NaCl system. Note: this is not true NaCl crystal structure, but we had to pick somethihg for the atom elements.To return the xyz string representing the system, we can use the `get_xyz` function, whic prints the state of the system into file in XYZ format.     std::string System::get_xyz(int fold,std::string pbc_type,int frame)    Here:    * **fold** Controlls the folding of the coordinates into the unit-cell   if fold==1 - will output coordinates folded into simulation box* **pbc_type** The parameter controlling the periodicity (when and if folding) of the unit cell  Can take values: "a", "b", "c", "ab", "ac", "bc", and "abc"* **frame** The integer index (in MD - the MD step, for instance) to be printed in the label of each record<jupyter_code>xyz = syst1.get_xyz(0, "abc", 0)
view = py3Dmol.view(width=800,height=400)  # linked=False,viewergrid=(3,2)
view.setBackgroundColor('0xeeeeee')
view.zoomTo()
view.addModel(xyz,'xyz',{'vibrate': {'frames':10,'amplitude':1.0}})
view.setStyle({'sphere':{'colorscheme':'Jmol'}})
view.show()<jupyter_output><empty_output><jupyter_text>2.4. All atoms together - using crystal coordinates[Back to TOC](TOC)Alternatively, to create a supercell of a crystal, we can read a file in the xyz format, create replicas, and add them to the `System` object.Here, we are going to read a file containing unit cell of the MAPI system in crystal coordinates system. To do this, we use:<jupyter_code>help(build.read_xyz_crystal)
a = VECTOR(16.9843, 0.0, 0.0) # a.u.
b = VECTOR(0.0, 16.9843, 0.0) # a.u.
c = VECTOR(0.0, 0.0, 16.9843) # a.u.

L0, R0raw = build.read_xyz_crystal("perov_madjet.xyz", a, b, c, 1, 0)  # in Bohr

print(L0)
for r in R0raw:
    print(r.x, r.y, r.z)<jupyter_output>['C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'Pb', 'Pb', 'Pb', 'Pb', 'Pb', 'Pb', 'Pb', 'Pb', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I', 'I']
6.487873926104878 7.943980795283024 7.445056766278085
6.844447225710289 7.943046811195377 23.706608782335106
6.491808855628773 23.83864396239548 7.445056766278085
6.8483821552341855 23.83771318787515 23.706608782335106
22.535710483591227 7.943980795283024 7.445056766278085
22.89228378319664 7.943046811195377 23.706608782335106
22.539645413115124 23.83864396239548 7.445056766278085
22.896218712720533 23.83771318787515 23.706608782335106
9.971608896484868 7.946917549373045 8.90758918921385
10.32818219609028 [...]<jupyter_text>Note that since we were reading the atomic coordiantes are fractional coordinates in teerms of unit cell lattice parameters, we had to defined the latter. The output Cartesian coordinates of the resulting atoms to in Bohrs (a.u.), wheres the lattice translational vectors are give in Angstroms. Thus we set the parameters to the above function as `inp_units=1` and `out_units=0`.<jupyter_code>L0, R0raw = build.read_xyz_crystal("perov_madjet.xyz", a, b, c)  # in Bohr

R0 = []
for r in R0raw:
    R0.append( r )   # coordinates in a.u.
    print( r.x, r.y, r.z)<jupyter_output>3.4332352753100004 4.20377390242 3.93975465735
3.62192575302 4.20327965929 12.544998015760001
3.43531755049 12.6148680291 3.93975465735
3.6240080282000005 12.614375484400002 12.544998015760001
11.925385275310001 4.20377390242 3.93975465735
12.11407575302 4.20327965929 12.544998015760001
11.92746755049 12.6148680291 3.93975465735
12.1161580282 12.614375484400002 12.544998015760001
5.276748562770001 4.2053279658700005 4.7136935413200005
5.465439040480001 4.20483372274 13.318936899730001
5.27883083795 12.616422092550001 4.7136935413200005
5.46752131566 12.61592784942 13.318936899730001
13.768898562770001 4.2053279658700005 4.7136935413200005
13.95758904048 4.20483372274 13.318936899730001
13.77098083795 12.616422092550001 4.7136935413200005
13.959671315660001 12.61592784942 13.318936899730001
3.5380199141600004 4.20239477726 2.45918906717
3.72671039187 4.20190223256 11.064449409880002
3.54010218934 12.613473618070001 2.45918906717
3.72879436548 12.61297937494 11.064449409880002
12.0301699[...]<jupyter_text>As before, we first create the needed number of replicas of the unit cell, then we define the required "masses" variales for all atoms, and finally we add them all to the `System` object. Note, we use the `add_atoms_to_system`, instead of the `add_atom_to_system` function considered above:<jupyter_code>help(build.add_atoms_to_system)<jupyter_output>Help on function add_atoms_to_system in module libra_py.build:

add_atoms_to_system(syst, L, R, scl, mass, univ_file)
    Args:    
        syst ( System object ): the chemical system to which we are going to add the atoms
        L ( list of strings ): element names    
        R ( list of VECTOR objects ): coordinates of the particles [ Bohrs ]
        scl ( VECTOR ): momentum rescaling factors along each direction [a.u. of momentum]
        mass ( list of doubles ): atomic masses [a.u. of mass]
        univ_file ( string ): name of the file that contains the Universe properties
    
    Returns:
        None: but adds new atoms to the System object, modifies the syst variable<jupyter_text>To setup the supercell properties of the created `System` object, we use the `init_box` member function.    // This function automatically determines the size of the box    // based on the system size. The box will be cubical or parallelepiped    void System::init_box()        void System::init_box(VECTOR tv1,VECTOR tv2,VECTOR tv3){    /**      tv1 The vector defining the "a" unit cell direction      tv2 The vector defining the "b" unit cell direction      tv3 The vector defining the "c" unit cell direction     This function creats the box based on argument vectors    */           void System::init_box(double maxx,double maxy,double maxz){    /**      maxx The length of the "a" (Cartesian X) unit cell dimension      maxy The length of the "b" (Cartesian Y) unit cell dimension      maxz The length of the "c" (Cartesian Z) unit cell dimension     This function creates the parralelepiped box with dimensions     given by the arguments         */<jupyter_code>Nx, Ny, Nz = 1, 1, 1
L1, R1 = build.generate_replicas_xyz2(L0, R0, a, b, c, Nx, Ny, Nz)

masses = []
for i in range(len(L1)):
    masses.append(1.0)

# Create the system
syst2 = System()
build.add_atoms_to_system(syst2, L1, R1, VECTOR(1.0, 1.0, 1.0), masses, "elements.dat")
syst2.init_box(a*Nx, b*Ny, c*Nz)<jupyter_output><empty_output><jupyter_text>Finally, we can print out the coordinates of atoms and the supercell translation vectors into the .ent format file recognized by packages such as VMD using the `print_ent` function, which prints the state of the system into a file in Brookhaven PDB format (ENT). The `print_ent` comese with the following signatures:    void System::print_ent(std::string filename)        void System::print_ent(std::string filename,int fold,std::string pbc_type)        void System::print_ent(std::string filename,boost::python::list atoms_list)    void System::print_ent(std::string filename,boost::python::list atoms_list,int fold,std::string pbc_type)    Here:* **filename** is the name of the file where the info will be printed out* **fold** controlls the folding of the coordinates into the unit-cell    if fold==1 is will output coordinates folded into simulation box, if 0 - periodicity is not accounted for* **pbc_type** the parameter controlling the periodicity (if `fold==1`) of the unit cell  Can take values: "a", "b", "c", "ab", "ac", "bc", and "abc"* **atoms_list** is the list of the atom IDs for the atoms to be printed out (atom IDs start from 0)The versions that don't check for `fold` and `pbc_type` options, it is assumes that `fold==1` and `pbc_type=="abc"`. The versions that don't check for `atoms_list` assume the operation is done for the entire system.We will test all 4 version to potentially see some difference:<jupyter_code># And print it out in the .ent format. Note, this will only print out the coordinates, and not the connectivities
# The entire system with the 3D folding
syst2.print_ent("perov_madjet-1.ent") 

# The entire system folded only in a and b directions
syst2.print_ent("perov_madjet-2.ent", 1, "ab") 

# Print only methylammonium molecules, fold in 3 dimentions
syst2.print_ent("perov_madjet-3.ent", list(range(0, 64)) ) 

# Print only methylammonium molecules, but fold only in c direction
syst2.print_ent("perov_madjet-4.ent", list(range(0, 64)), 1, "c")<jupyter_output><empty_output><jupyter_text>Having the periodicity information in this files is important for visualizing the periodic images of the supercell in the VMD programThe visualized structures are as follows (this was done outside of this tutorial file):|     |  Case 1 |  Case 3 | |-----|---------|---------|| Original | ![](pm-1-orig.png)   |   ![](pm-3-orig.png) || VMD replicas | ![](pm-1-xyz.png)  | ![](pm-3-xyz.png)  | Alternatively, we can obtain the geometries in the xyz format using the `get_xyz` funtion, as before<jupyter_code>xyz = syst2.get_xyz(0, "abc", 0)
view = py3Dmol.view(width=800,height=400)  # linked=False,viewergrid=(3,2)
view.setBackgroundColor('0xeeeeee')
view.zoomTo()
view.addModel(xyz,'xyz',{'vibrate': {'frames':10,'amplitude':1.0}})
view.setStyle({'sphere':{'colorscheme':'Jmol'}})
view.show()<jupyter_output><empty_output><jupyter_text>2.5. All atoms together - using Cartesian coordinates[Back to TOC](TOC)Here, we use the `read_xyz` function we already encountered. We use it to read in the atomic coordinates of the unit cell (in Angstrom)We then set the unit cell vectors (also in Angstrom)Then, we generate a 2x1x1 supercell, add atoms to `System` object, setup the supercell box, and print the coordinates in the .ent file format.<jupyter_code>L0, R0raw = build.read_xyz("perov_weili.xyz")  # in Bohr

a = VECTOR(8.8009, 0.0,      0.0)*units.Angst     # conversion to a.u.
b = VECTOR(0.0,    8.009,    0.0)*units.Angst     # conversion to a.u.
c = VECTOR(0.0,    0.0,  12.6857)*units.Angst     # conversion to a.u.

R0 = []
for r in R0raw:
    R0.append( r )   # coordinates in a.u.
    #print(r.x, r.y, r.z)
    

Nx, Ny, Nz = 1, 1, 1
L1, R1 = build.generate_replicas_xyz2(L0, R0, a, b, c, Nx, Ny, Nz, inp_units=0, out_units=0)


masses = []
for i in range(len(L1)):
    masses.append(1.0)

# Create the system
syst3 = System()
build.add_atoms_to_system(syst3, L1, R1, VECTOR(1.0, 1.0, 1.0), masses, "elements.dat")
syst3.init_box(a*Nx,b*Ny,c*Nz)

# And print it out in the .ent format. Note, this will only print out the coordinates, and not the connectivities
syst3.print_ent("perov_weili-0.ent", 0, "abc")     
syst3.print_ent("perov_weili-1.ent", 1, "abc")<jupyter_output><empty_output><jupyter_text>Note that in this case we demonstrated the difference between doing and not doing the folding. It makes the original pictures different (because of the way VMD draws the bonds is affected by the atomic coordination), but the physical meaning of files is the same. The visual differences become negligible when periodic images are shown.After VMD visualization, we obtain:|              |   No folding        |   With folding      ||--------------|---------------------|---------------------|| Original     | ![](pwl-0-orig.png) | ![](pwl-1-orig.png) || VMD replicas | ![](pwl-0-xyz.png)  | ![](pwl-1-xyz.png)  |<jupyter_code>xyz = syst3.get_xyz(0, "abc", 0)
view = py3Dmol.view(width=800,height=400)  # linked=False,viewergrid=(3,2)
view.setBackgroundColor('0xeeeeee')
view.zoomTo()
view.addModel(xyz,'xyz',{'vibrate': {'frames':10,'amplitude':1.0}})
view.setStyle({'sphere':{'colorscheme':'Jmol'}})
view.show()<jupyter_output><empty_output><jupyter_text>3. Determining connectivities[Back to TOC](TOC)In this section, we demonstrate the utilities for finding the connectivities between atoms, which is important when it comes to quantum dot building and atom replacements discussed below. The connectivities are also important for setting up the classical force field calculations, because the declaration that two atoms are connected together makes the interaction to be of the "bonded" type as opposed to the "non-bonded"We consider than simple case of Si. Let's load the xyz file<jupyter_code>L, R = build.read_xyz("Si.xyz")
print( L, R)
xyz = build.make_xyz(L, R)
view = py3Dmol.view(width=800,height=400)  # linked=False,viewergrid=(3,2)
view.setBackgroundColor('0xeeeeee')
view.zoomTo()
view.addModel(xyz,'xyz',{'vibrate': {'frames':10,'amplitude':1.0}})
view.setStyle({'sphere':{'colorscheme':'Jmol'}})
view.show()

print( xyz)<jupyter_output><empty_output><jupyter_text>Next, we define the coordination number for all atom types (elements).The `PT_coord` dictionary defines the maximal coordinations for each element and the following loop uses this info to setup the maximal coordination numbers for all atoms<jupyter_code>PT_coord = { "Si": 4 }
MaxCoord = []
for i in L:
    MaxCoord.append(PT_coord[i])    # maximal coordination number of each atom in super-cell<jupyter_output><empty_output><jupyter_text>Next, we use the `autoconnect` function to determine the connectivities<jupyter_code>help(autoconnect.autoconnect)<jupyter_output>Help on function autoconnect in module libra_py.autoconnect:

autoconnect(R, MaxCoord, params)
    Args:
        R ( list of VECTOR objects ): The atomic coordinates of the system [ units: arbitrary ]
        MaxCoord ( list of ints ): Maximal coorination numbers of each atom 
        params ( dictionary ): Parameters controlling the execution of this function
    
            * **params["Rcut"]** ( double ): The maximal radius of connectivity [ units: same as R, default: 0.0 ]
            * **params["tv1"]** ( VECTOR ): unit cell translation vector a [ units: same as R ]
            * **params["tv2"]** ( VECTOR ): unit cell translation vector b [ units: same as R ]
            * **params["tv3"]** ( VECTOR ): unit cell translation vector c [ units: same as R ]
            * **params["pbc_opt"]** ( string ): what type of periodicity to assume. 
                Options: "a", "b", "c", "ab", "ac", "bc", "abc", "none" [ default: None ]
            * **params["opt"]** ( int ): Option to obe[...]<jupyter_text>3.1. Non-periodic case[Back to TOC](TOC)<jupyter_code>params = {"Rcut":5.001 * units.Angst, "pbc_opt":"none", "opt":0, "verbosity":0  }

res, line, pairs = autoconnect.autoconnect(R, MaxCoord, params)<jupyter_output><empty_output><jupyter_text>Here we looked for all pairs of atoms that are closer to each other that 5.001 Angstrom.We treat the loaded structure as non-periodic, so the atoms do not "see" the periodic images which they may be connected to. Also, the atoms are not connected if the connections would disobey the maximal coordination number specified.The first returned result shows the connectivity of each atom<jupyter_code>print(res)<jupyter_output>[[0, [1, 2, 3]], [1, [0, 2, 3]], [2, [0, 1, 3]], [3, [0, 1, 2]]]<jupyter_text>The second one shows the connectivities in a special format that can be added to the .ent file and understandable by Libra's force field setups. This format follows to the one produced by the HyperChem software<jupyter_code>print(line)<jupyter_output>CONECT     1      2      3      4 
CONECT     2      1      3      4 
CONECT     3      1      2      4 
CONECT     4      1      2      3<jupyter_text>Finally, the 2 last returned variables shows the periodic translation of first and second atoms in the pair needed for the connectivity to take place:<jupyter_code>for pair in pairs:
    t0, t1 = pair[2], pair[3]
    print(pair[0], pair[1], t0.x, t0.y, t0.z,  t1.x, t1.y, t1.z)<jupyter_output>0 1 0.0 0.0 0.0 0.0 0.0 0.0
0 2 0.0 0.0 0.0 0.0 0.0 0.0
0 3 0.0 0.0 0.0 0.0 0.0 0.0
1 2 0.0 0.0 0.0 0.0 0.0 0.0
1 3 0.0 0.0 0.0 0.0 0.0 0.0
2 3 0.0 0.0 0.0 0.0 0.0 0.0<jupyter_text>Together with the information on the maximal coordination allowed for all atoms, the `pairs` result can be used to determine the undercoordinated atoms, using the `find_undercoordinated_atoms` function.<jupyter_code>help(autoconnect.find_undercoordinated_atoms)
out = autoconnect.find_undercoordinated_atoms(res, MaxCoord)

print(out)<jupyter_output>[[0, 1], [1, 1], [2, 1], [3, 1]]<jupyter_text>Here, we see that there are in total 4 undercoodinated atoms (indiced 0 through 3). Each of them is within 5.001 Angstrom to 3 other atoms (each other here). As a result, each atom has 1 unsaturated valency (dangling bond).  3.2. Periodic case (Warning: a bit buggy for now)[Back to TOC](TOC)Let's now consider what happens when we allow connections with the periodic images.In this case, we also need to provide information about the unit cell (or supercell) translation vectors.<jupyter_code>a = VECTOR( 5.4307100000000000,    0.0000000000000000,    0.0000000000000000 ) * units.Angst
b = VECTOR( 0.0000000000000000,    5.4307100000000000,    0.0000000000000000 ) * units.Angst
c = VECTOR( 0.0000000000000000,    0.0000000000000000,    5.4307100000000000 ) * units.Angst

params = {"Rcut":5.001 * units.Angst,
          "pbc_opt":"abc", "tv1":a, "tv2":b, "tv3":c,
          "opt":0, "verbosity":0  }

res, line, pairs = autoconnect.autoconnect(R, MaxCoord, params)

print(res)
for pair in pairs:
    t0, t1 = pair[2], pair[3]
    print(pair[0], pair[1], t0.x, t0.y, t0.z,  t1.x, t1.y, t1.z)
out = autoconnect.find_undercoordinated_atoms(res, MaxCoord)

print(out)<jupyter_output>[]<jupyter_text>In this case, we find not undercoordinated atoms. However, the connectivity is yet to be fixed!  4. Quantum Dot building[Back to TOC](TOC)In this section, we will construct quantum dots. We will start by making the supercells and then cropping the created slabs to spheres of needed radii, then optionally replace some of the atoms on the QD surface, while cropping. 4.1. With capping[Back to TOC](TOC)We start by creating a supercell, as we did in the [section 2](2.).<jupyter_code>a = VECTOR(5.833, 0.0, 0.0) * units.Angst
b = VECTOR(0.0, 5.833, 0.0) * units.Angst
c = VECTOR(0.0, 0.0, 5.833) * units.Angst

L0, R0 = build.read_xyz("CdSe.xyz")
L1, R1 = build.generate_replicas_xyz2(L0, R0, a, b, c, 6, 6, 6)<jupyter_output><empty_output><jupyter_text>Next, we'll be connecting atoms in the super-cell to each other, based on their separation distancesand the maximal coordination numbers acceptedAt this point, we are going to define the maximal coordination number of each atom type present in the system (here, Cd and Se), as well as for each additional elements - e.g. those which we are going to capthe QD with (here, H)<jupyter_code>PT_coord = { "Cd": 4, "Se": 4, "H":1 }
MaxCoord1 = []
masses = []
for i in L1:
    MaxCoord1.append(PT_coord[i])    # maximal coordination number of each atom in super-cell
    masses.append(1.0)               # masses of each atom in the super-cell, do not matter here<jupyter_output><empty_output><jupyter_text>The parameter that defines the atomic connectivities: if the atoms are separatedthe the distance shorter than this, and if their actual coordination number has notreached the maximal one yet, we will assume they are covelently bonded (in FF there willbe a harmonic potential term for instance)<jupyter_code>Rbond = 3.001 * units.Angst<jupyter_output><empty_output><jupyter_text>Now connect all atoms in the supercell. The "none" indicates that there is no periodic boundary conditions to check when doing such connections<jupyter_code>params = {"Rcut":Rbond, 
          "tv1":Nx*a, "tv2":Ny*b, "tv3":Nz*c,
          "pbc_opt":"none", "opt":0, "verbosity":0
         }
res, line, pairs = autoconnect.autoconnect(R1, MaxCoord1, params)<jupyter_output><empty_output><jupyter_text>The radius of the QD<jupyter_code>Rdot = 7.0 * units.Angst<jupyter_output><empty_output><jupyter_text>Finally, we can carve a QD out of the super-cell, using `crop_sphere_xyz3` function<jupyter_code>help(build.crop_sphere_xyz3)<jupyter_output>Help on function crop_sphere_xyz3 in module libra_py.build:

crop_sphere_xyz3(L, R, Rcut, pairs, new_L)
    This function removes all atoms that are outside of a sphere of 
    Rcut radius. The sphere is centered on GEOMETRIC center of the 
    system
    
    Args: 
        L ( list of strings ): element names, this list will be trimmed accordingly
        R ( VECTORList or list of VECTOR objects): coordinates of the particles [ Bohr ]
        Rcut ( double ): the radius of the sphere from the center of the QD [ Bohrs ]
        pairs ( list of [int, int, VECTOR, VECTOR] lists ): integers describe the indices for 
            the connected atoms, the VECTORs describe the translation vector by which the atoms
            should be translated before considering connected (e.g. for periodic boundary conditions)
            these VECTOR objects are not used in the present function, but the format
            is kept such that the output of other functions could be used in this function.
  [...]<jupyter_text>All the atoms that are outside the shepre of Rdot radiusare removed, except for the atoms that are bound (as encoded in the "pairs" argument) to the atoms insidethe sphere. The bound outside atoms' identities will be switched according to which atoms there are connected.This is encoded in the last (dictionary) argument. Here, it means that all the atoms outside the sphere thatare connected to the Cd inside will become H, all the outside atoms that are connected to the inside Se willbe converted to Cd.In other words, all Se atoms will be capped by Cd and all Cd atoms will be capped with H.The L and R variables will contain the labels and coordinates of the final (QD) atoms.<jupyter_code>L, R = build.crop_sphere_xyz3(L1, R1, Rdot, pairs, {"Cd":"H", "Se":"Cd"})

xyz = build.make_xyz(L, R)
view = py3Dmol.view(width=800,height=400)  # linked=False,viewergrid=(3,2)
view.setBackgroundColor('0xeeeeee')
view.zoomTo()
view.addModel(xyz,'xyz',{'vibrate': {'frames':10,'amplitude':1.0}})
view.setStyle({'sphere':{'colorscheme':'Jmol'}})
view.show()<jupyter_output><empty_output><jupyter_text>Now, we are going to auto-connect all the atoms in the final QDDefine the maximal coordination for the atoms<jupyter_code>MaxCoord2 = []
for i in L:
    MaxCoord2.append(PT_coord[i])   # maximal coordination number of each atom in the resulting QD<jupyter_output><empty_output><jupyter_text>Determine the connectivities of the new set of atoms<jupyter_code>params = {"Rcut":Rbond, 
          "tv1":Nx*a, "tv2":Ny*b, "tv3":Nz*c,
          "pbc_opt":"none", "opt":0, "verbosity":0
         }
res, line, pairs = autoconnect.autoconnect(R, MaxCoord2, params)

print("Autodetermined connections: ")
print(line)<jupyter_output>Autodetermined connections: 
CONECT     1     32     37     45     53 
CONECT     2     34     40     46     54 
CONECT     3     35     41     48     55 
CONECT     4     36     42     49     56 
CONECT     5     43     50     68     69 
CONECT     6     44     58     71     72 
CONECT     7     52     60     75     76 
CONECT     8     29     55     80     81 
CONECT     9     30     37     46     56 
CONECT    10     31     48     57     83 
CONECT    11     32     39     49     58 
CONECT    12     33     41     59     86 
CONECT    13     34     42     51     60 
CONECT    14     36     44     52     89 
CONECT    15     29     46     92     93 
CONECT    16     30     47     54     94 
CONECT    17     31     37     49     55 
CONECT    18     32     38     50     56 
CONECT    19     33     40     51     98 
CONECT    20     35     42     52     59 
CONECT    21     36     43     60    101 
CONECT    22     29     37    103    104 
CONECT    23     30     38     53    105 
CONEC[...]<jupyter_text>Create the chemical system object and save it in the .ent form suitable for force field calculations<jupyter_code># Create an empty chemical system object
syst = System()

# Add all the atoms of the final QD to the chemical system object
build.add_atoms_to_system(syst, L, R, VECTOR(1.0, 1.0, 1.0), masses, "elements.dat")

# And print it out in the .ent format. Note, this will only print out the coordinates, and not the connectivities
syst.print_ent("CdSe-qd-capped.ent") 

# We print out the connectivities afterwards manually
f = open("CdSe-qd-capped.ent", "a")
f.write("END\n")
f.write(line)
f.write("END\n")
f.close()<jupyter_output><empty_output><jupyter_text>Here, we visualize the VMD picture of the created QD![](CdSe-qd-capped.png)  4.2. Without capping[Back to TOC](TOC)This example shows how to build a CdSe quantum dot (QD) of a given size.Here, we don't cap the surface atoms, but still generate the connectivity map for the final QD to save the file in the .ent format for futher use in MD simulationsFirst, let's repeat many steps of the workflow above:<jupyter_code>a = VECTOR(5.833, 0.0, 0.0) * units.Angst
b = VECTOR(0.0, 5.833, 0.0) * units.Angst
c = VECTOR(0.0, 0.0, 5.833) * units.Angst

L0, R0 = build.read_xyz("CdSe.xyz")
L1, R1 = build.generate_replicas_xyz2(L0, R0, a, b, c, 6, 6, 6)

PT_coord = { "Cd": 4, "Se": 4, "H":1 }
MaxCoord1 = []
masses = []
for i in L1:
    MaxCoord1.append(PT_coord[i])    # maximal coordination number of each atom in super-cell
    masses.append(1.0)      
    
Rbond = 3.001 * units.Angst
Rdot = 7.0 * units.Angst

# Finally, we can carve a QD out of the super-cell. All the atoms that are outside the shepre of Rdot radius
# are removed. The L and R variables will contain the labels and coordinates of the final (QD) atoms.
L, R = build.crop_sphere_xyz2(L1, R1, Rdot)<jupyter_output><empty_output><jupyter_text>Note that we simply skipped the step of determining the connectivities and atom substitutions before we do the cropping.For this reason, we use the version `2` of the cropping, instead of the version `3`. It does not need any information about the connectivities and the atom replacement rules.<jupyter_code>help(build.crop_sphere_xyz2)<jupyter_output>Help on function crop_sphere_xyz2 in module libra_py.build:

crop_sphere_xyz2(L, R, Rcut)
    This function removes all atoms that are outside of a sphere of 
    Rcut radius. The sphere is centered on GEOMETRIC center of the 
    system
    
    Args: 
        L ( list of strings ): element names, this list will be trimmed accordingly
        R ( VECTORList or list of VECTOR objects): coordinates of the particles [ Bohr ]
        Rcut ( double ): the radius of the sphere from the center of the QD [ Bohrs ]
    
    Returns:
        tuple: (lab, coords), where:
    
            * lab ( list of N strings ): the labels of all remaining atoms
            * coords ( list of N VECTORs ): the coordinates of all remaining atoms<jupyter_text>The rest of the workflow is as before<jupyter_code># Now, we are going to auto-connect all the atoms in the final QD
# Define the maximal coordination for the atoms 
MaxCoord2 = []
for i in L:
    MaxCoord2.append(PT_coord[i])   # maximal coordination number of each atom in the resulting QD

    
params = {"Rcut":Rbond, 
          "tv1":Nx*a, "tv2":Ny*b, "tv3":Nz*c,
          "pbc_opt":"none", "opt":0, "verbosity":0
         }

# Determine the connectivities of the new set of atoms
res, line, pairs = autoconnect.autoconnect(R, MaxCoord2, params)
print("Autodetermined connections: ")
print(line)


syst = System()
build.add_atoms_to_system(syst, L, R, VECTOR(1.0, 1.0, 1.0), masses, "elements.dat")
syst.print_ent("CdSe-qd-uncapped.ent") 
f = open("CdSe-qd-uncapped.ent", "a")
f.write("END\n")
f.write(line)
f.write("END\n")
f.close()<jupyter_output>Autodetermined connections: 
CONECT     1     32     37     45     53 
CONECT     2     34     40     46     54 
CONECT     3     35     41     48     55 
CONECT     4     36     42     49     56 
CONECT     5     43     50 
CONECT     6     44     58 
CONECT     7     52     60 
CONECT     8     29     55 
CONECT     9     30     37     46     56 
CONECT    10     31     48     57 
CONECT    11     32     39     49     58 
CONECT    12     33     41     59 
CONECT    13     34     42     51     60 
CONECT    14     36     44     52 
CONECT    15     29     46 
CONECT    16     30     47     54 
CONECT    17     31     37     49     55 
CONECT    18     32     38     50     56 
CONECT    19     33     40     51 
CONECT    20     35     42     52     59 
CONECT    21     36     43     60 
CONECT    22     29     37 
CONECT    23     30     38     53 
CONECT    24     31     39     45 
CONECT    25     33     42     46     55 
CONECT    26     34     43     47     56 
CONECT    27     35[...]<jupyter_text>Visualize the produced QD![](CdSe-qd-uncapped.png)  4.3. Streamlined building of QDs[Back to TOC](TOC)This example demonstrates how to build a quantum dot (QD) of a given radius.Here, we first generate the xyz with the supercell, using `generate_replicas_xyz`<jupyter_code>help(build.generate_replicas_xyz)<jupyter_output>Help on function generate_replicas_xyz in module libra_py.build:

generate_replicas_xyz(tv1, tv2, tv3, rep1, rep2, rep3, filename, outfile, inp_units=0, out_units=0)
    This function generates a lattice of coordinates and an array of 
    corresponding atomic labels by replicating a given set of labeled
    atoms periodically in up to 3 dimensions (with a variable number of 
    replicas in each direction). This function first reads in the 
    atomic coordinates (supplied in an .xyz format) and then writes 
    the resulting data into another file (.xyz format)
    
    Args:
        tv1 ( VECTOR ): translation vector in direction 1 [defined by inp_units] 
        tv2 ( VECTOR ): translation vector in direction 2 [defined by inp_units] 
        tv3 ( VECTOR ): translation vector in direction 3 [defined by inp_units] 
        rep1 ( int ): the number of replications along the vector tv1, not counting the original cell
        rep2 ( int ): the number of replications along the vector t[...]<jupyter_text>And then use even simpler than before cropping function `crop_sphere_xyz`<jupyter_code>help(build.crop_sphere_xyz)
## Si QD
a = [     5.4307100000000000,    0.0000000000000000,    0.0000000000000000 ]
b = [     0.0000000000000000,    5.4307100000000000,    0.0000000000000000 ]
c = [     0.0000000000000000,    0.0000000000000000,    5.4307100000000000 ]
build.generate_replicas_xyz(a, b, c, 5, 5, 5 , "Si.xyz", "si-5.xyz")
build.crop_sphere_xyz("si-5.xyz", "si_5.xyz", 5.0)<jupyter_output>Molecule<jupyter_text>Visualizing in VMD:|  The cubic supercell   |   The carved out sphere ||------------------------|-------------------------|| ![](si-5.png)          | ![](si-5-qd.png)        |<jupyter_code>## CdSe QD
a = [5.833, 0.0, 0.0]
b = [0.0, 5.833, 0.0]
c = [0.0, 0.0, 5.833]
build.generate_replicas_xyz(a, b, c, 20, 20, 20 , "CdSe.xyz", "cdse-20.xyz")
label, R = build.crop_sphere_xyz("cdse-20.xyz", "cdse_15.xyz", 15.0)

label, R = build.crop_sphere_xyz("cdse-20.xyz", "cdse_25.xyz", 25.0)

label, R = build.crop_sphere_xyz("cdse-20.xyz", "cdse_35.xyz", 35.0)<jupyter_output>Molecule 
 Molecule 
 Molecule
{'fraction_non_alphanumeric': 0.11423583709974412, 'fraction_numerical': 0.0760217136571251, 'mean_word_length': 3.8366874855424475, 'pattern_counts': {'":': 52, '<': 178, '<?xml version=': 0, '>': 177, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 18, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '42', 'dir': 'jupyter-structured-clean-dedup', 'id': '536209', 'license': 'CC0-1.0', 'n_tokens_mistral': 16394, 'n_tokens_neox': 14181, 'n_words': 4746, 'path': '12_molecular_builders/1_crystal_and_qd_builder/tutorial.ipynb', 'repo_name': 'langerest/Tutorials_Libra', 'size': '41817'}