#*********************************************************************************
#* Copyright (C) 2017-2019 Brendan Smith, Ekadashi Pradhan, 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:: autoconnect
:platform: Unix, Windows
:synopsis: This module implements functions to determine the connectivity
matrix in molecular systems.
.. moduleauthor:: Brendan Smith, Ekadashi Pradhan, Alexey V. Akimov
"""
import os
import sys
import math
import copy
import unittest
if sys.platform=="cygwin":
from cyglibra_core import *
elif sys.platform=="linux" or sys.platform=="linux2":
from liblibra_core import *
#import common_utils as comn
import util.libutil as comn
[docs]def 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 obey maximal coordination number:
- 0: the connected atoms in a pair have to obey both the coordination numbers
Some of the atoms may stay undercoordinated. This is the consistent scheme in
the sense that if A is connected to B, then B is necessarily connected to A [default]
- 1: treat the MaxCoord as the minimal coordination number each atom should attain.
Some atoms may be overcoordinated.
* **params["verbosity"]** ( int ): Flag to control the amount of the printout [ default: 0 ]
Returns:
tuple: ( res, line, unsorted_pairs_out ), where:
* res ( list of lists [ res[i][0], res[0][1] ] ), where
* res[i][0] ( int ): index of the atom
* res[i][1] ( list of ints ): indices of the atoms that are connected to res[i][0]
* line ( string ): the connectivity information in the .ent format
* unsorted_pairs_out ( same type as ```res```): same as ```res```, but not sorted
"""
critical_params = [ ]
default_params = { "Rcut":0.0, "pbc_opt":"none", "opt":0,
"tv1":VECTOR(1.0e+10, 0.0, 0.0),
"tv2":VECTOR(0.0, 1.0e+10, 0.0),
"tv3":VECTOR(0.0, 0.0, 1.0e+10),
"verbosity":0
}
comn.check_input(params, default_params, critical_params)
Rcut = params["Rcut"]
pbc_opt = params["pbc_opt"]
verbosity = params["verbosity"]
opt = params["opt"]
tv1 = params["tv1"]
tv2 = params["tv2"]
tv3 = params["tv3"]
# Preprocessing - sanity check
N = len(R)
if len(MaxCoord) != N:
print("Error: The length of the MaxCoord list should be the same as the lenght of the R list!")
print("Exiting now...")
sys.exit(0)
unsorted_pairs = []
mapping = []
periodicity = []
if pbc_opt not in ["a", "b", "c", "ab", "ac", "bc", "abc","none"]:
print("Error: pbc_opt ", pbc_opt, " is not recognized")
sys.exit(0)
transl_a = [0.0]
transl_b = [0.0]
transl_c = [0.0]
if pbc_opt in ["a", "ab", "ac", "abc"]:
transl_a = [-1.0, 0.0, 1.0]
if pbc_opt in ["b", "ab", "bc", "abc"]:
transl_b = [-1.0, 0.0, 1.0]
if pbc_opt in ["c", "ac", "bc", "abc"]:
transl_c = [-1.0, 0.0, 1.0]
# Distances between all the pairs
count = 0
for i in range(0,N):
for j in range(i+1,N):
for n1 in transl_a:
for n2 in transl_b:
for n3 in transl_c:
T = n1 * tv1 + n2 * tv2 + n3 * tv3
r = (R[i]-R[j]-T).length()
if r<Rcut:
unsorted_pairs.append([count, r])
mapping.append([i,j])
periodicity.append([n1,n2,n3])
count += 1
# Sort all the pairs according to the interparticle distance
sorted_pairs = merge_sort(unsorted_pairs)
unsorted_pairs_out = []
for it in unsorted_pairs:
unsorted_pairs_out.append( [ mapping[it[0]][0], mapping[it[0]][1],
VECTOR(0.0, 0.0, 0.0),
VECTOR( periodicity[it[0]][0] * tv1 + periodicity[it[0]][1] * tv2 + periodicity[it[0]][2] * tv3 ) ]
)
if verbosity==1:
print("Pairs of atoms separated no more than by Rcut")
print(unsorted_pairs)
print("Formatted printout:")
for it in unsorted_pairs:
print(F"Atoms {mapping[it[0]][0]} and {mapping[it[0]][1]} are separated by\
{it[1]:8.5f}. Periodic translation applied = {periodicity[it[0]][0]:8.5f},\
{periodicity[it[0]][1]:8.5f}, {periodicity[it[0]][2]:8.5f} )" )
print("Pairs of atoms separated no more than by Rcut, sorted by the distance")
print(sorted_pairs)
print("Formatted printout:")
for it in sorted_pairs:
print(F"Atoms {mapping[it[0]][0]} and {mapping[it[0]][1]} are separated by\
{it[1]:8.5f}. Periodic translation applied = {periodicity[it[0]][0]:8.5f},\
{periodicity[it[0]][1]:8.5f}, {periodicity[it[0]][2]:8.5f} )" )
# Initialize the results variables
act_coord = [0]*N # actual coordination numbers
num_sat = 0 # the number of valence-saturated atoms
res = []
for i in range(0,N):
res.append([i,[]])
# Now lets pick the pairs with the minimal distance, depending on the option...
for it in sorted_pairs:
if num_sat < N:
i,j = mapping[it[0]][0], mapping[it[0]][1]
if opt==0: #... no overcoordinated atoms
if act_coord[i]<MaxCoord[i] and act_coord[j]<MaxCoord[j]:
res[i][1].append(j)
res[j][1].append(i)
act_coord[i] += 1
act_coord[j] += 1
if act_coord[i]==MaxCoord[i]:
num_sat += 1
if act_coord[j]==MaxCoord[j]:
num_sat += 1
elif opt==1: # ... no undercoordinated atoms
if act_coord[i]<MaxCoord[i]:
res[i][1].append(j)
act_coord[i] += 1
if act_coord[j]<MaxCoord[j]:
res[j][1].append(i)
act_coord[j] += 1
# Order the indices of the atoms to which each atom is connected (for testing and convenience)
for i in range(0,N):
res[i][1] = sorted(res[i][1])
line = ""
for it in res:
line = line + "CONECT %5i " % (it[0]+1)
for it2 in it[1]:
line = line + " %5i " % (it2+1)
line = line + "\n"
if verbosity==1:
print("res = ", res)
print("Formatted output:")
print(line)
return res, line, unsorted_pairs_out
[docs]def find_undercoordinated_atoms(res, MaxCoord):
"""
Args:
res ( list of lists [ res[i][0], res[0][1] ] ), where
* res[i][0] ( int ): index of the atom
* res[i][1] ( list of ints ): indices of the atoms that are connected to res[i][0]
This would typically be the first output of the ```autoconnect``` function
MaxCoord ( list of ints ): Maximal coorination numbers of each atom
Returns:
(list): out, such that:
out[i][0] - is the index of the i-th undercoordinated atom
out[i][1] - is the number of dangling bonds on the i-th atom
"""
nat = len(MaxCoord) # the total number of atoms
out = []
for i in range(0,nat):
diff = MaxCoord[i] - len(res[i][1])
if diff > 0:
out.append([i, diff])
return out
def example_1():
rnd = Random()
N = 10
R = []
MaxCoord = []
for i in range(0,N):
if rnd.uniform(0.0, 1.0)<0.5:
MaxCoord.append(3)
else:
MaxCoord.append(3)
R.append(VECTOR(rnd.normal(), rnd.normal(), rnd.normal()))
Rcut = 3.0
print("Runing autoconnect(R, MaxCoord, Rcut, 0, 1) with : ")
print("R = ", R)
print("MaxCoord = ", MaxCoord)
print("Rcut = ", Rcut)
autoconnect(R, MaxCoord, {"Rcut":Rcut, "opt":0, "verbosity":1})
[docs]class TestAutoconnect(unittest.TestCase):
""" Docs
"""
[docs] def test_1(self):
""" Diatomic"""
R = []
MaxCoord = []
R.append(VECTOR(0.0, 0.0, 0.0)); MaxCoord.append(1) # Atom 1
R.append(VECTOR(1.0, 0.0, 0.0)); MaxCoord.append(1) # Atom 2
res, lines = autoconnect(R, MaxCoord, {"Rcut":1.5})
self.assertEqual(res, [[0, [1]], [1, [0]]] )
res, lines = autoconnect(R, MaxCoord, {"Rcut":0.5})
self.assertEqual(res, [[0,[]], [1, []]] )
[docs] def test_2(self):
""" Triatomic """
R = []
MaxCoord = []
R.append(VECTOR(0.0, 0.0, 0.0)); MaxCoord.append(2) # Atom 1
R.append(VECTOR(1.0, 0.0, 0.0)); MaxCoord.append(2) # Atom 2
R.append(VECTOR(1.2, 0.0, 0.0)); MaxCoord.append(2) # Atom 2
res, lines = autoconnect(R, MaxCoord, {"Rcut":1.5})
self.assertEqual(res, [[0, [1,2]], [1, [0,2]], [2, [0,1]] ] )
res, lines = autoconnect(R, MaxCoord, {"Rcut":0.5})
print(res)
self.assertEqual(res, [ [0, []], [1, [2]], [2, [1]] ] )
if __name__=='__main__':
unittest.main()
#example_1()