Source code for libra_py.build

#*********************************************************************************                     
#* Copyright (C) 2017-2019 Alexey V. Akimov                                                   
#*                                                                                                     
#* This file is distributed under the terms of the GNU General Public License                          
#* as published by the Free Software Foundation, either version 2 of                                   
#* the License, or (at your option) any later version.                                                 
#* See the file LICENSE in the root directory of this distribution   
#* or <http://www.gnu.org/licenses/>.          
#***********************************************************************************
"""
.. module:: build
   :platform: Unix, Windows
   :synopsis: This module implements functions for building molecular structures

.. moduleauthor:: Alexey V. Akimov

"""

import os
import sys
import math
import copy

if sys.platform=="cygwin":
    from cyglibra_core import *
elif sys.platform=="linux" or sys.platform=="linux2":
    from liblibra_core import *

from . import LoadPT
from . import units


[docs]def read_xyz(filename, inp_units=1, out_units=0): """Read an xyz file (single structure) Args: filename ( string ): the name of the xyz file to read inp_units ( int ): defines the coordinates' units in the input file: * 0 - Bohr * 1 - Angstrom ( default ) out_units ( int ): defines the output coordinates' units: * 0 - Bohr ( default ) * 1 - Angstrom Returns: tuple: (L, coords), where: * L ( list of N strings ): the labels of N atoms read * coords ( list of N VECTORs ): the coordinates of N atoms read [defined by out_units] """ f = open(filename,"r") A = f.readlines() f.close() tmp = A[0].split() Natoms = int(float(tmp[0])) L = [] coords = [] for a in A[2:]: tmp = a.split() if len(tmp)==4: x = float(tmp[1]) y = float(tmp[2]) z = float(tmp[3]) R = VECTOR(x,y,z) * units.length_converter(inp_units, out_units) L.append(tmp[0]) coords.append(R) return L, coords
[docs]def 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 """ nat = len(L) res = "%i\n" % (nat) res = res + "\n" for i in range(0,nat): res = res + "%s %12.8f %12.8f %12.8f \n" % (L[i], R[i].x/units.Angst, R[i].y/units.Angst, R[i].z/units.Angst) return res
[docs]def make_xyz_mat(E, R): """ This function returns a string in the xyz format with X, Y, Z where X,Y,Z are the coordinates Args: E ( list of ndof/3 string ): atom names (elements) of all atoms R ( MATRIX(ndof x nsteps-1) ): coordinates of all DOFs for all mid-timesteps Returns: string: A string representing an xyz file """ natoms = len(E) nsteps = R.num_of_cols res = "" for step in range(nsteps): res = res + "%3i\nStep = %6i\n" % (natoms, step) for i in range(0,natoms): x,y,z = R.get(3*i, step), R.get(3*i+1, step), R.get(3*i+2, step) res = res + "%s %5.3f %5.3f %5.3f\n" % (E[i], x,y,z) return res
[docs]def read_xyz_crystal(filename, a,b,c, inp_units=0, out_units=0): """Read an xyz file with atomic coordinates given in fractional lattice vector units Args: filename ( string ): the name of the xyz file to read a ( VECTOR ): the unit cell vector in the direction a [units defined by inp_units] b ( VECTOR ): the unit cell vector in the direction b [units defined by inp_units] c ( VECTOR ): the unit cell vector in the direction c [units defined by inp_units] inp_units ( int ): defines the units of a,b,c: * 0 - Bohr ( default ) * 1 - Angstrom out_units ( int ): defines the units of the coordinates in the ourput variable: * 0 - Bohr ( default ) * 1 - Angstrom Returns: tuple: (L, coords), where: * L ( list of N strings ): the labels of N atoms read * coords ( list of N VECTORs ): the coordinates of N atoms read [defined by out_units] """ M = MATRIX(3,3) M.init(a,b,c) M = M.T() f = open(filename,"r") A = f.readlines() f.close() tmp = A[0].split() Natoms = int(float(tmp[0])) L = [] coords = [] for a in A[2:]: tmp = a.split() if len(tmp)==4: x = float(tmp[1]) y = float(tmp[2]) z = float(tmp[3]) R = M * VECTOR(x,y,z) * units.length_converter(inp_units, out_units) L.append(tmp[0]) coords.append(R) return L, coords
[docs]def 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 tv2, not counting the original cell rep3 ( int ): the number of replications along the vector tv3, not counting the original cell filename ( string ): the name of the file containing input coordinates (original unit cell) outfile ( string ): the name of the file where the resulting lattice atoms will be written inp_units ( int ): defines the units of variables in filename, tv1, tv2, and tv3: * 0 - Bohr ( default ) * 1 - Angstrom out_units ( int ): defines the units of the coordinates written to outfile * 0 - Bohr ( default ) * 1 - Angstrom Returns: None: but a new .xyz file will be printed out Example: The example below read in a unit cell of Si (8 atoms), replicates in 5 times in x, y, and z directions (so we get a 6 x 6 x 6 supercell). The resulting structure is written to the "cube10.xyz" file >>> a = [ 5.4307100000000000, 0.0000000000000000, 0.0000000000000000 ] >>> b = [ 0.0000000000000000, 5.4307100000000000, 0.0000000000000000 ] >>> c = [ 0.0000000000000000, 0.0000000000000000, 5.4307100000000000 ] >>> generate_replicas_xyz(a, b, c, 5, 5, 5 , "Si.xyz", "cube10.xyz") """ transl = [] for x in range(0,rep1+1): for y in range(0,rep2+1): for z in range(0,rep3+1): transl.append([x,y,z]) f = open(filename,"r") A = f.readlines() f.close() tmp = A[0].split() Natoms = int(float(tmp[0])) f = open(outfile, "w+") f.write(F"{int(len(transl) * Natoms):i} \n") f.write(F"{A[1][:-1]:s} \n") L = [] coords = [] for T in transl: for a in A[2:]: tmp = a.split() if len(tmp)==4: x = float(tmp[1]) y = float(tmp[2]) z = float(tmp[3]) # Now apply PBC dx = T[0]*tv1[0] + T[1]*tv2[0] + T[2]*tv3[0] dy = T[0]*tv1[1] + T[1]*tv2[1] + T[2]*tv3[1] dz = T[0]*tv1[2] + T[1]*tv2[2] + T[2]*tv3[2] f.write(F"{tmp[0]:s} {x+dx:8.5f} {y+dy:8.5f} {z+dz:8.5f} \n") L.append(tmp[0]) coords.append(VECTOR(x+dx, y+dy, z+dz)) f.close() return L, coords
[docs]def generate_replicas_xyz2(L, R, tv1, tv2, tv3, Nx, Ny, Nz, 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) Args: L ( list of strings ): atomis labels R ( list of VECTOR objects ): initial atomic coords [defined by inp_units] 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] Nx ( int ): the number of replications along the vector tv1, not counting the original cell Ny ( int ): the number of replications along the vector tv2, not counting the original cell Nz ( int ): the number of replications along the vector tv3, not counting the original cell inp_units ( int ): defines the units of variables R, tv1, tv2, and tv3: * 0 - Bohr ( default ) * 1 - Angstrom out_units ( int ): defines the units of the coordinates returned: * 0 - Bohr ( default ) * 1 - Angstrom Returns: tuple: (lab, newR), where: * lab ( list of N strings ): the labels of all atoms and their corresponding replicas * newR ( list of N VECTORs ): the coordinates of all replicated atoms [defined by out_units] """ nat = len(R) lab, newR = [], [] scl = units.length_converter(inp_units, out_units) for i in range(0,nat): for nx in range(0,Nx): for ny in range(0,Ny): for nz in range(0,Nz): # Now apply PBC r = VECTOR(R[i] + nx * tv1 + ny * tv2 + nz * tv3) * scl lab.append(L[i]) newR.append(r) return lab, newR
[docs]def crop_sphere_xyz(infile, outfile, Rcut): """ This function reads an .xyz file with the geometry, cuts the atoms that are outside of a sphere of given radius and then prints out the remaining atoms to a new file in .xyz format Args: infile ( string ): the name of the .xyz file that contains the original coordinates outfile ( string ): the name of the new .xyz file to create (with the cropped geometry) Rcut ( double ): the radius of the sphere from the center of the QD [ same units as used in infile ] 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 Example: The example below read in a unit cell of Si (8 atoms), replicates in 5 times in x, y, and z directions (so we get a 6 x 6 x 6 supercell). The resulting structure is written to the "cube10.xyz" file. Then the file cube10.xyz is read and cropped to make a sphere of 5 Angstrom. The resulting system is then printed out to the "qd_5.xyz" file. >>> a = [ 5.4307100000000000, 0.0000000000000000, 0.0000000000000000 ] >>> b = [ 0.0000000000000000, 5.4307100000000000, 0.0000000000000000 ] >>> c = [ 0.0000000000000000, 0.0000000000000000, 5.4307100000000000 ] >>> generate_replicas_xyz(a, b, c, 5, 5, 5 , "Si.xyz", "cube10.xyz") >>> crop_sphere_xyz("cube10.xyz", "qd_5.xyz", 5.0) """ f = open(infile,"r") A = f.readlines() f.close() tmp = A[0].split() Natoms = int(float(tmp[0])) print(A[1][:-1]) # Read coordinates and compute COM L = [] R = [] Rcom = [0.0,0.0,0.0] for a in A[2:]: tmp = a.split() if len(tmp)==4: x = float(tmp[1]) y = float(tmp[2]) z = float(tmp[3]) Rcom[0] = Rcom[0] + x Rcom[1] = Rcom[1] + y Rcom[2] = Rcom[2] + z L.append(tmp[0]) R.append([x,y,z]) # Geometric center Rcom[0] = Rcom[0] / Natoms Rcom[1] = Rcom[1] / Natoms Rcom[2] = Rcom[2] / Natoms # Compute remaining number of atoms Nat = 0 indx = [] for i in range(0,Natoms): r = math.sqrt((R[i][0] - Rcom[0])**2 + (R[i][1] - Rcom[1])**2 + (R[i][2] - Rcom[2])**2) if(r<=Rcut): Nat = Nat + 1 indx.append(i) # Print resulting coordinates f1 = open(outfile,"w") f1.write(F"{Nat:5i}\n") f1.write(A[1]) coords = [] lab = [] for i in indx: coords.append(VECTOR(R[i][0], R[i][1], R[i][2])) lab.append(L[i]) f1.write(F"{L[i]} {R[i][0]:12.6f} {R[i][1]:12.6f} {R[i][2]:12.6f} \n") f1.close() return lab, coords
[docs]def 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 """ # Compute COM Rcom = VECTOR(0.0, 0.0, 0.0) for r in R: Rcom = Rcom + r # Geometric center Natoms = len(R) Rcom = Rcom / Natoms # Compute remaining number of atoms Nat = 0 indx = [] for i in range(0,Natoms): r = (R[i] - Rcom).length() if(r<=Rcut): Nat = Nat + 1 indx.append(i) # Prepare the coordinates and labels of the remaining atoms coords = [] lab = [] for i in indx: coords.append(R[i]) lab.append(L[i]) return lab, coords
[docs]def 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. new_L ( dictionary(string:string) ): a map containing the key-value pairs, where the key string corresponds to the label of the atom that is inside the cropped sphere. The value string corresponds to the new label of the atom outside the sphere. The function will use the original connectivity information to decide how to cap the dangling bonds formed at cropping. So, if the atom A is connected to atom B, atom A is inside of the sphere and B is outside the sphere. The vaule element will become a new label for atom B, if the key element corresponds to atom type A. 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 """ # Compute COM Rcom = VECTOR(0.0, 0.0, 0.0) for r in R: Rcom = Rcom + r # Geometric center Natoms = len(R) Rcom = Rcom / Natoms # Compute remaining number of atoms Nat = 0 indx = [] for i in range(0,Natoms): r = (R[i] - Rcom).length() if(r<=Rcut): Nat = Nat + 1 indx.append(i) # Prepare the coordinates and labels of the remaining atoms coords = [] lab = [] for i in indx: coords.append(R[i]) lab.append(L[i]) # Now add a number of atoms that are ourside the sphere added = [] # indices of the added atoms - we need them to avoid a multiple placing for it in pairs: if (it[0] in indx and it[1] in indx): # both atoms are within the sphere pass elif (it[0] in indx and it[1] not in indx): # it[0] is inside, it[1] is outside if it[1] not in added: new_lab = new_L[ L[it[0]] ] # new label if new_lab == "none": pass else: coords.append(R[it[1]]) lab.append(new_lab) added.append(it[1]) elif (it[1] in indx and it[0] not in indx): # it[1] is inside, it[0] is outside if it[0] not in added: new_lab = new_L[ L[it[1]] ] # new label if new_lab == "none": pass else: coords.append(R[it[0]]) lab.append(new_lab) added.append(it[0]) return lab, coords
[docs]def 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 """ rnd = Random() nat = syst.Number_of_atoms # Associate the coordinates with a molecular system U = Universe() LoadPT.Load_PT(U, univ_file, 0) sz = len(R) for i in range(0,sz): syst.CREATE_ATOM( Atom(U, {"Atom_element":L[i],"Atom_cm_x":R[i].x,"Atom_cm_y":R[i].y,"Atom_cm_z":R[i].z}) ) # Masses syst.Atoms[nat+i].Atom_RB.set_mass(mass[i]) # Momenta p = VECTOR(scl.x*rnd.normal(), scl.y*rnd.normal(), scl.z*rnd.normal()) syst.Atoms[nat+i].Atom_RB.set_momentum(p);
[docs]def 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 a ( VECTOR ): the periodicity translation vector along the a direction [Bohr] b ( VECTOR ): the periodicity translation vector along the b direction [Bohr] c ( VECTOR ): the periodicity translation vector along the c direction [Bohr] shift ( VECTOR ): essentially a coordinate of the atom in the unit cell [Bohr] elt ( string ): Element name mass ( double ): mass of the atom [a.u. of mass, 1 = mass of an electron] scl ( VECTOR ): momentum rescaling factors along each direction [a.u. of momentum] max_coord ( int ): maximal coordination number of the added atom allowed rnd ( Random ): an instance of the random number generator class inp_units ( int ): defines the units of variables stored in shift, a,b, and c * 0 - Bohr ( default ) * 1 - Angstrom Returns: None: But the following objects will be updated: * syst - the new atoms will be added * coords - the new coordinates will be added * MaxCoords - the maximal coordination numbers for new atoms will be added """ nat = syst.Number_of_atoms # Associate the coordinates with a molecular system U = Universe() LoadPT.Load_PT(U, "elements.dat", 0) i = 0 out_units = 0 # internal units in System objects are atomic units for nx in range(0,Nx): for ny in range(0,Ny): for nz in range(0,Nz): R = VECTOR(nx * a + ny * b + nz * c + shift) * units.length_converter(inp_units, out_units) coords.append(R) MaxCoords.append(max_coord) syst.CREATE_ATOM( Atom(U, {"Atom_element":elt,"Atom_cm_x":R.x,"Atom_cm_y":R.y,"Atom_cm_z":R.z}) ) # Masses syst.Atoms[nat+i].Atom_RB.set_mass(mass) # Momenta p = VECTOR(scl.x*rnd.normal(), scl.y*rnd.normal(), scl.z*rnd.normal()) syst.Atoms[i].Atom_RB.set_momentum(p); i = i + 1
[docs]def make_system(R, E, step, U): """ This function creates a chemical system (System) oject from the corrdinates and atom labels provided. All atoms of the system are made into a single group. Args: R ( MATRIX(ndof, nsteps) ): molecular coordinates at different times The convention is R.get(dof, step) is the coordinate of the ```dof```-th degree of freedom and ```step```-th timestep. Note: ndof = 3*natoms. [ a.u. = Bohrs ] E ( list of ```natoms``` strings): atomic labels step ( int ): selector of the timestep for which we want to construct the system U ( Universe object ): Universe Returns: System : the chemical system with the coordinates of all atoms at a given time. All atoms are grouped together into a single rigid body """ ndof = R.num_of_rows nat = int(ndof/3) syst = System() for at in range(0,nat): x = R.get(3*at+0, step) y = R.get(3*at+1, step) z = R.get(3*at+2, step) syst.CREATE_ATOM( Atom(U, {"Atom_element":E[at], "Atom_cm_x":x, "Atom_cm_y":y, "Atom_cm_z":z } ) ) syst.GROUP_ATOMS(list(range(1, nat+1)), 1) syst.init_fragments() return syst