'''define Structure object

'''

def __init__(self):

self.atom = [] # S_atom list for all atom in Str

self.energy = 0 # Energy of a Str

self.Latt = [] # lattice parameter

self.abc = self.Latt # same pointer?

self.Cell = [] # abc-in-xyz, POSCAR-abc, transfer-matrix

self.natom = 0 # number of atoms in Str

self.Ele_Name =[] # element name list for all atom in str

self.Ele_Index = [] # element index of all atom in str

# Ele_Index are also named iza in ZPLiu group？

self.Coord = [] # Coordination of each atom

self.sp ={} # atom and their number

self.sporder= {} # atom ordered by element

self.Cell = [] # lattice abc in Cartesian-3D

self.natompe= []

self.eleList =[] # specific element index in Str

self.ele_nameList = [] # specific element in Str

# relevant to force

self.Lfor = False # have allfor.arc to read or not

self.For= [] # list contains force of Str in each-atom-3D

self.stress = [] # stress list ?

# self.serial_num = 0

self.maxF = 0

self.max_stress = 0

self.frac = []

self.strlen = 0

self.centerf = {}

self._cycle=[]

self.cart = []

self.Q= []

def addatom(self, line: str, flag = 1 ):

'''add one atom from structure format file line

Args:

line(string): one line in structure file

flag(int): structure file format flag, default is 1

1 for arc and MS-xsd,

2 for ?,

3 for cat file,

4 for mol file,

5 for QM9,

6 for ? ,

'''

if flag == 1:

#for arc and ms

coord = [float(x) for x in line.split()[1:4]]

ele_symbol = line.split()[0] # get element symbol

if ele_symbol not in PT.Eletable:

# for ms.car file

for element in PT.Eletable:

if element in ele_symbol:

ele_symbol = element

break

else:

raise FooError("element in arc-structure not found!")

ele_number = PT.Eledict.get(ele_symbol,0) # get element number using symbol

single_atom = S_atom(coord, ele_number) #

# get charge

try:

single_atom.charge = float(line.split()[-2])

except :

single_atom.charge = 0.0

# append to atom list

self.atom.append(single_atom)

return

# not important for LASP-VASP training

elif flag == 2:

# for TrainStr

coord = [float(x) for x in line.split()[2:5]]

ele_number = int(line.split()[1])

self.atom.append(S_atom(coord,ele_number))

elif flag == 3:# for cat file

coord = [float(x) for x in line.split()[1:4]]

ele_symbol = str(line.split()[0])[0]

ele_number = PT.Eledict.get(ele_symbol,0)

self.atom.append(S_atom(coord,ele_number))

elif flag == 4:#for mol file

coord = [float(x) for x in line.split()[0:3]]

ele_symbol = line.split()[3]

ele_number = PT.Eledict(ele_symbol,0)

self.atom.append(S_atom(coord,ele_number))

elif flag == 5: #for QM9

coord = [float(x) for x in line.split()[1:4]]

ele_symbol = line.split()[0]

ele_number = PT.Eledict(ele_symbol,0)

self.atom.append(S_atom(coord,ele_number))

elif flag == 6:

coord = [float(x) for x in line.split()[1:4]]

ele_symbol = line.split()[-2]

ele_number = PT.Eledict(ele_symbol,0)

single_atom= S_atom(coord,ele_number)

single_atom.charge =float(line.split()[-1])

self.atom.append(single_atom)

def add_force(self, line :str, serial_num, flag=1):

'''flag 1 for arc-for file 2 for trainfor-file'''

if flag == 1:

self.atom[serial_num].force = [float(x) for x in line.split()]

elif flag == 2:

self.atom[serial_num].force= [float(x) for x in line.split()[2:5]]

def add_stress(self, line: str, flag=1):

'''flag 1 for arc-file 2 for traindata-file'''

if flag == 1: self.stress = [float(x) for x in line.split()]

elif flag == 2: self.stress = [float(x) for x in line.split()[1:7]]

def add_charge(self, base):

for atom in self.atom:

atom.charge = base[atom.ele_symbol]

def set_coord(self):

'''set self.Coord from self.atom, do before sort_by_ele'''

self.Coord = []

for atom in self.atom:

# atom: S_atom

self.Coord.append(atom.xyz)

self.Coord = np.array(self.Coord)

def set_for_list(self):

'''set self.For from self.atom, do before sort_by_ele'''

self.For = []

for atom in self.atom:

self.For.append(atom.force)

self.For = np.array(self.For)

def element_to_list(self, element_count: dict):

'''transfer element_count dict (Str.sp) to all_element list'''

element_list = []

ele_indexs = []

for element,count in element_count.items():

for i in range(count):

element_list.append(element)

ele_indexs.append(PT.Eledict[element])

all_element = element_list

all_ele_ind = ele_indexs

return all_element, all_ele_ind

def build_coord_set_from_POSCAR(self, filename="POSCAR"):

'''read POSCAR to get cell and atom coord to build a Str

'''

try:

f = open(filename, 'r')

except:

print('----No POSCAR info----')

else:

print("read POSCAR")

index = 0

atom_id = 0

coord_type_status = False

Cart = 0

for line in f:

L_list = line.strip().split()

index += 1

if index > 2 and index < 6:

# read Cell info

self.Cell.append([np.float64(x) for x in L_list])

if index == 6:

# read element name

# self.Latt = self.Cell2Latt(Cell=self.Cell)

elements = L_list

if index == 7:

# read element count and get all_element list

for i in range(len(L_list)):

self.ele_nameList.append(elements[i])

self.eleList.append(PT.Eledict[elements[i]])

self.sp[elements[i]] = int(L_list[i])

# get Ele_Name and Ele_Index

self.Ele_Name, self.Ele_Index = self.element_to_list(self.sp)

# problem: something will add line

if index > 7 and coord_type_status == False:

# read coord type

if (L_list[0][0]=="C" or L_list[0][0]=="c"):

Cart=1 # Cartesian coord

coord_type_status = True

if (L_list[0][0]=="D" or L_list[0][0]=="d"):

Cart=0# frac coord

coord_type_status = True

if coord_type_status == True:

# read atom coord

if L_list == []: break # pick the end

if len(L_list) < 3: continue

if Cart == 1:

cart_coord = [np.float64(x) for x in L_list[0:3]]

self.Coord.append(cart_coord)

self.frac = self.FracCoord()

if Cart == 0:

frac_coord = [np.float64(x) for x in L_list[0:3]]

cart_coord = np.dot([np.float64(x) for x in L_list[0:3]], self.Cell)

self.frac.append(frac_coord)

self.Coord.append(cart_coord)

atom_obj = S_atom(cart_coord, self.Ele_Index[atom_id])

self.atom.append(atom_obj)

atom_id += 1

# self.sort_atom_by_element # not need

# self.natom = atom_id

self.Latt = self.Cell2Latt()

self.set_strinfo_from_atom()

f.close()

return

def get_max_force(self):

maxf= 0

for atom in self.atom:

_tmpf = max([abs(x) for x in atom.force])

if maxf < _tmpf:

maxf = _tmpf

self.maxF = maxf

def get_max_stress(self):

self.max_stress = max(self.stress)

def set_strinfo_from_atom(self, ):

'''set eleList, natompe, ele_nameList, Ele_Index, Ele_Name from Str.atom

'''

self.natom = len(self.atom)

self.eleList,self.natompe = list(np.unique([atom.ele_num for atom in self.atom],return_counts=True))

self.ele_nameList = [PT.Eletable[ele_num-1] for ele_num in self.eleList]

self.sp = {}

self.sporder = {}

for index, ele_name in enumerate(self.ele_nameList):

order_count = 1

self.sp[ele_name] = self.natompe[index]

self.sporder[ele_name] = order_count

order_count += 1

self.nele = len(self.eleList)

self.Ele_Index = [atom.ele_num for atom in self.atom]

self.Ele_Name = [atom.ele_symbol for atom in self.atom]

def screen_upper_surf(self,):

self.upper =0

for atom in self.atom:

if atom.xyz[2] > 0.9*self.Latt[2]:

self.upper =1

break

def sort_atom_by_element(self,):

'''sort atom in stucture by atom.ele_num, for normalize sturcture'''

self.atom.sort(key =lambda X: X.ele_num)

def sort_atom_by_z(self,):

self.atom.sort(key= lambda X: X.xyz[2])

def calc_two_dim_coord(self, ):

self.xa = np.array([atom.xyz for atom in self.atom])

def calc_one_dim_coord(self, ):

self.cal_two_dim_coord()

self.xa_onedim = reduce(lambda a,b: a+b, [list(x) for x in self.xa])

def add_atom_ID(self):

for iatom,atom in enumerate(self.atom):

atom.id = iatom

def cdnt2fcnt(self, ):

'''set frac Coord used by self.fdnt'''

self.cal_two_dim_coord()

latinv = np.linalg.inv(self.Cell)

self.fdnt = [list(x) for x in np.matmul(self.xa, latinv)]

# self.fdnt

def calc_centroid(self):

# centroid = mass center

self.set_strinfo_from_atom()

xyzall =np.array([0,0,0])

for atom in self.atom:

xyzall= xyzall+np.array(atom.xyz)

self.centroid =xyzall/self.natom

def mv_str_to_boxcenter(self):

#for Ortholat

self.calc_centroid()

center =np.array([self.Latt[0]/2,self.Latt[1]/2,self.Latt[2]/2])

mv = center - self.centroid

for i,atom in enumerate(self.atom):

atom.xyz = list(self.xa[i]+mv)

def add_ortho_lat(self):

#self.calAtomnum()

if self.natom != 0:

self.cal_two_dim_coord()

max3d=np.amax(self.xa, axis=0)

min3d=np.amin(self.xa, axis=0)

delta =max(list(max3d-min3d))

a = 10

while ((a-delta) < 4):

a=a+5

self.Latt =[a,a,a,90,90,90]

def outPOSCAR(self,outfile="POSCAR"):

'''print structure to POSCAR file'''

f=open(outfile,'w')

f.write('system\n')

f.write('1.000000000000\n')

for item in self.Cell:

f.write('%12.8f %12.8f %12.8f\n'%(item[0],item[1],item[2]))

f.write("%s\n" %reduce(lambda a,b:a+b , ["%4s"%s for s in self.ele_nameList]))

f.write("%s\n" %reduce(lambda a,b:a+b , ["%4d"%num for num in self.natompe] ))

f.write('Cart\n')

for atom in self.atom:

f.write('%12.8f %12.8f %12.8f\n'%(atom.xyz[0],atom.xyz[1],atom.xyz[2]))

f.close()

def genPOTCAR(self,sourcedir,outfile="POTCAR"):

'''get needed POTCAR from sourcedir, POTCAR will be print named outfile'''

os.system('rm -rf %s'%outfile)

for ele in self.ele_nameList:

os.system('cat %s/POTCAR.%s >> %s'%(sourcedir,ele,outfile))

def genKPOINTS(self,outfile="KPOINTS", criteria=25):

'''get needed kpoints, ka = kb = kc = criteria

writed by JamesBourbon

'''

ka = int(m.ceil(criteria / float(self.Latt[0])))

kb = int(m.ceil(criteria / float(self.Latt[1])))

kc = int(m.ceil(criteria / float(self.Latt[2])))

with open(outfile, 'w') as fout:

fout.write('mesh auto\n')

fout.write('0\n')

fout.write('G\n')

fout.write('%d %d %d\n'%(ka,kb,kc))

fout.write('0 0 0')

def set_element_radius(self, dataset_name = ''):

"""Setting element radius data in structure, by JamesBourbon

Default radii dataset from

Chem. Eur. J. 2009, 15, 186–197, DOI: 10.1002/chem.200800987

Args:

dataset_name (str): Defaults using Eleradii in PeriodicTable.

Returns:

radius_dict (dict): Ele_Index for key and Ele_radii for value

"""

radius_dict = {}

if bool(dataset_name):

# read from external csv file if given

try:

radius_dataset = pd.read_csv(dataset_name)

for ele_index in set(self.Ele_Index):

radii = np.float64(radius_dataset[

radius_dataset["index"] == ele_index ]["radius(A)"])

radius_dict[ele_index] = radii

except:

for ele_index in set(self.Ele_Index):

radii = PT.Eleradii[ele_index-1] # notice!

radius_dict[ele_index] = radii

else:

for ele_index in set(self.Ele_Index):

radii = PT.Eleradii[ele_index-1]

radius_dict[ele_index] = radii

return radius_dict

def calc_dis(self, iatom1, iatom2):

'''calc periodic atom distance, considered periodic in x-np.round(x)

just consider one cell'''

self.cdnt2fcnt() # get fractional coord

vbond =np.array(self.fdnt[iatom1])- np.array(self.fdnt[iatom2])

dis = np.linalg.norm(np.matmul( np.array([x-np.round(x) for x in vbond]), self.Cell))

return dis

def calc_neighbour(self,iatom):

'''calc distance of one atom to all other atom

cannot consider trans-cell coordination.'''

dict ={}

for i in range(self.natom):

if (i!= iatom):

dis= self.calc_dis(iatom,i)

dict[i]= dis

return dict

def make_supercell(self,x_max:int,y_max:int,z_max:int,

x_min:int=0, y_min:int=0, z_min:int=0):

'''make supercell by x, y, z, running well

returns:

supercell: list of (ele_index, coordinate)

'''

supercell = []

for i in range(x_min,x_max+1):

for j in range(y_min,y_max+1):

for k in range(z_min,z_max+1):

dims = [i,j,k]

# get move vector: dims@self.Cell, moving along a,b,c

moved_str = ((self.Ele_Index[index],

coord+np.matmul(dims,self.Cell))

for index, coord in enumerate(self.Coord))

supercell += moved_str

return supercell

# need to be more : this just calc one Str, coordination pattern analysis should in allstr or other

# coding in 20220609

def coordination_pattern(self, tol=1.21,):

'''calc coordination pattern of Str, interface to coordination_pattern.py

coor_patterns example: {(Pd, ((Au, 2.4),)), (Au, ((Pd, 2.4),))}

Returns:

coor_patterns: coordination patterns format set

'''

# setting

radius_dict = self.set_element_radius()

# judge each dim is thin or not

thin_edge_list = [max(radius_dict.values())*i*2*tol+0.2 for i in range(1,5)]

shell_k_list = [1.2, 0.6, 0.4, 0.3, 0.25] # test result

shell_k = 1.2

# last k for large lattice parameter, cutoff ref: 1.0, 0.5, 0.34, 0.25, 0.20

# Judge min bond length

min_bond = min(PT.Eleradii) * tol

# 3x3x3 supercell

lat_xyz = np.diag(self.Cell) # describe a shell

supercell_333 = self.make_supercell(1,1,1,-1,-1,-1)

# distance calculated from central cell

coor_patterns = set()

for ci,central_coord in enumerate(self.Coord):

central_ele = self.Ele_Name[ci]

coor_pattern_dict = {}

for atom_info in supercell_333:

# atom_info: (ele_index, atom_xyz)

do_calc = True

square_dist_array = []

for dim, vector in enumerate(lat_xyz):

# to simplify calculation: not to far in any dim

for i, edge in enumerate(thin_edge_list):

if vector < edge:

shell_k = shell_k_list[i]

break

else:

shell_k = shell_k_list[-1]

square_dist_dim = np.power(central_coord[dim] - atom_info[1][dim], 2)

if square_dist_dim > np.power(vector*shell_k, 2):

do_calc = False

break

else:

do_calc = True

square_dist_array.append(square_dist_dim)

# calc coordination main

if do_calc:

bond_dist = np.sqrt(np.sum(square_dist_array))

central_index = self.Ele_Index[ci]

central_radii = radius_dict[central_index]

coor_radii = radius_dict[atom_info[0]]

# coor-bond main judgement

max_bond = (central_radii + coor_radii) * tol

if (bond_dist > max_bond) or (bond_dist < min_bond):

continue

else:

coor_ele = PT.Eletable[atom_info[0]-1]

coor_dist = np.round(bond_dist,1)

coor_pair = (coor_ele, coor_dist)

# get CN

coor_pattern_dict[coor_pair] = coor_pattern_dict.get(coor_pair,0) + 1

# end of coor_bond calc.

# transfer to Set(tuple) for coor-patterns format

coor_atoms = []

for coor_pair, count in coor_pattern_dict.items():

coor_atoms.append((coor_pair[0], coor_pair[1], count))

coor_atoms = tuple(coor_atoms)

coor_pattern_i = (central_ele, coor_atoms)

coor_patterns.add(coor_pattern_i)

return coor_patterns # can be directly used in CoordinationPatterns.patterns

# end of coordination pattern method embedding by JamesBourbon

def special_neighbour(self,iatom,spe_ele):

'''calc distance of one atom to all atom of one element'''

dict ={}

for i in range(self.natom):

if (i!= iatom and self.atom[i].ele_num == spe_ele):

dis= self.calc_dis(iatom,i)

dict[i]= dis

return dict

def longest_bond(self,ele1,ele2,lim):

'''calc longest bond in cell'''

self.cdnt2fcnt()

maxlist = []

for i,atom in enumerate(self.atom):

if atom.ele_num == ele1:

result = self.special_neighbour(i,ele2)

dis = [x for x in result.values() if x<lim]

if len(dis) > 0:

maxlist.append(max(dis))

longest = max(maxlist)

return longest

def shortest_bond(self,ele1,ele2,lim):

'''calc shortest bond in cell'''

self.cdnt2fcnt()

minlist = []

for i,atom in enumerate(self.atom):

if atom.ele_num == ele1:

result = self.special_neighbour(i,ele2)

dis = [x for x in result.values() ]

if len(dis) > 0:

minlist.append(min(dis))

short = min(minlist)

return short

def simple_class(self,iatom):

for i in range(self.natom):

if (i != iatom):

#dislist.append(self.cal_distance(iatom,i))

#print dis

dis2 = self.calc_dis(iatom,i)

#print dis

bondtest = bond(self.atom[iatom].ele_num,self.atom[i].ele_num,dis2)

bondorder = bondtest.judge_bondorder()

if bondorder != 0:

if self.atom[i].ele_num == 8:

self.atom[iatom].bondtype.append(101)

if self.atom[i].ele_num == 6:

self.atom[iatom].bondtype.append(11)

if self.atom[i].ele_num == 1:

self.atom[iatom].bondtype.append(1)

self.atom[iatom].bondtype.sort( )

return

def bond_search(self,iatom,ele2,flag=1):

bondlist =[]

for i in range(self.natom):

if (i != iatom) and self.atom[i].ele_num ==ele2:

if flag ==2 and i<iatom and self.atom[iatom].ele_num ==ele2: continue

#dislist.append(self.cal_distance(iatom,i))

dis = self.calc_dis(iatom,i)

bondtest = bond(self.atom[iatom].ele_num,self.atom[i].ele_num,dis)

bondorder = bondtest.judge_bondorder()

if bondorder!= 0:

bondlist.append(dis)

return bondlist

def determine_species(self,iatom,elelist):

#dislist = []

for i in range(self.natom):

if (i != iatom):

#dislist.append(self.cal_distance(iatom,i))

dis = self.calc_dis(iatom,i)

bondtest = bond(self.atom[iatom].ele_num,self.atom[i].ele_num,dis)

bondorder = bondtest.judge_bondorder()

if bondorder != 0 :#and self.atom[i].ele_num == 8:

self.atom[iatom].bondlist.append(bondorder)

#for id,iza in enumerate(elelist):

# if self.atom[i].ele_num ==iza:

# if id >0:

# self.atom[iatom].bondtype.append(bondorder+10**id)

# else:

# self.atom[iatom].bondtype.append(bondorder)

# break

if self.atom[i].ele_num == elelist[3]:

self.atom[iatom].bondtype.append(bondorder+1000)

if self.atom[i].ele_num == elelist[2]:

self.atom[iatom].bondtype.append(bondorder+100)

if self.atom[i].ele_num == elelist[1]:

self.atom[iatom].bondtype.append(bondorder+10)

if self.atom[i].ele_num == elelist[0]:

self.atom[iatom].bondtype.append(bondorder)

self.atom[iatom].species= len(self.atom[iatom].bondlist)

bond_all = sum(self.atom[iatom].bondlist)

self.atom[iatom].bondtype.sort(reverse = True )

if bond_all != 4:

self.atom[iatom].Ctype = 'radical'

else :

if self.atom[iatom].bondtype[-1] == 103 :

self.atom[iatom].Ctype = 'CO'

elif self.atom[iatom].bondtype[-1] == 102 :

if self.atom[iatom].bondtype[-2] == 101 :

self.atom[iatom].Ctype = 'acid'

else:

self.atom[iatom].Ctype = 'ketone'

elif self.atom[iatom].bondtype[-1] == 101 :

self.atom[iatom].Ctype = 'alcohol'

elif self.atom[iatom].bondtype[-1] == 13 :

self.atom[iatom].Ctype = 'alkyne'

elif self.atom[iatom].bondtype[-1] == 12 :

self.atom[iatom].Ctype = 'alkene'

elif self.atom[iatom].bondtype[-1] == 11 :

self.atom[iatom].Ctype = 'alkane'

elif self.atom[iatom].bondtype[-1] == 1 :

self.atom[iatom].Ctype = 'CH4'

#self.atom[iatom].species =self.atom[iatom].species+ bondorder

#dislist.sort()

#length = min(len(dislist),4)

#for i in range(length):

# if dislist[i]> 1.7:

# break

#self.atom[iatom].species = i + 1

self.atom[iatom].bondall = bond_all

return

'''

def calc_Ctypes(self, ):

self.calc_one_dim_coord()

cell = reduce(lambda a,b: a+b, self.Cell)

#print cell

self.c_natm = pointer(ctypes.c_int(self.natom))

self.c_xa = pointer((ctypes.c_double*len(self.xa_onedim))(*self.xa_onedim))

self.c_rv = pointer((ctypes.c_double*len(cell))(*cell))

iza = [atom.ele_num for atom in self.atom]

self.c_iza = pointer((ctypes.c_int*self.natom)(*iza))

'''

def calc_frag_charge(self):

"charge here is fakecharge, actually group id"

groupdict ={}

for i,atom in enumerate(self.atom):

#print atom.charge

try:

groupdict[atom.charge].append(atom)

except:

groupdict[atom.charge] =[]

groupdict[atom.charge].append(atom)

self.frag= groupdict

def cal_group_frag(self):

groupdict ={}

for i,atom in enumerate(self.atom):

#print atom.charge

try:

groupdict[self.group[i]].append(atom)

except:

groupdict[self.group[i]] =[]

groupdict[self.group[i]].append(atom)

self.frag= groupdict

def determine_charge(self):

self.cal_group_frag()

for item in self.frag.values():

_tmpcharge= 0

_elelist =[]

for atom in item:

_elelist.append(atom.ele_num)

_tmpcharge =_tmpcharge + atom.ele_num

#print _elelist

if _tmpcharge%2 != 0:

for atom in item:

if atom.ele_num== 7:

if _elelist ==[7,8]:

return 15

if _elelist ==[7,8,8]:

return 23

return -1

return 1

return 2

def check_min(self,flag = 1):

sqnatm = self.natom**2

self.calc_Ctypes()

program='/home7/kpl/pymodule/Lib/Lib_fillbond/checkminbond.so'

Lminstr = pointer(ctypes.c_bool(0))

bmatrix = pointer((ctypes.c_int*sqnatm)(*[0 for i in range(sqnatm)]))

bondneed = pointer((ctypes.c_int*(self.natom))(*[0 for i in range(self.natom)]))

surface = pointer((ctypes.c_int*(self.natom))(*[0 for i in range(self.natom)]))

checkmin = ctypes.cdll.LoadLibrary(program)

if flag == 0:

checkmin.judgebond_(self.c_natm,self.c_iza,self.c_xa,self.c_rv,Lminstr,bmatrix,bondneed)

elif flag ==1 :

checkmin.judgebondsurface_(self.c_natm,self.c_iza,self.c_xa,self.c_rv,Lminstr,bmatrix,bondneed,surface)

bmx = list(bmatrix.contents)

#self.bmx2D = np.array(bmx).reshape(self.natom, self.natom)

#self.bmx1D = bmx

self.bondneed = list(bondneed.contents)

self.Lminstr = bool(Lminstr.contents)

return self.Lminstr,bmx,list(bondneed.contents),list(surface.contents)

def bond_matrix(self, ):

program='/home7/kpl/pymodule/Lib/Lib_bondmatrix/bondmatrix.so'

self.calc_Ctypes()

sqnatm = self.natom**2

self.cdnt2fcnt()

onedim_fa = reduce(lambda a,b: a+b, self.fdnt)

self.c_fxa = pointer((ctypes.c_double*len(onedim_fa))(*onedim_fa))

bmatrix = pointer((ctypes.c_int*sqnatm)(*[0 for i in range(sqnatm)]))

bcal = ctypes.cdll.LoadLibrary(program)

bcal.hbondmatrix_ (self.c_natm, self.c_fxa, self.c_iza, self.c_rv, bmatrix)

return list(bmatrix.contents)

def seg_molecular (self, recal=True):

""" warning : bond matrix must first be calculated."""

program='/home7/kpl/pymodule/Lib/Lib_bondmatrix/bondmatrix.so'

if recal:

self.calc_Ctypes()

self.cdnt2fcnt()

onedim_fa = reduce(lambda a,b: a+b, self.fdnt)

self.c_fxa = pointer((ctypes.c_double*len(onedim_fa))(*onedim_fa))

self.c_bmx = pointer((ctypes.c_int*len(self.bmx1D))(*self.bmx1D))

group = pointer((ctypes.c_int*self.natom)(*[0 for i in range(self.natom)]))

scal = ctypes.cdll.LoadLibrary(program)

scal.hsegmentmol_ (self.c_natm, self.c_fxa, self.c_iza, self.c_rv, group, self.c_bmx)

return list(group.contents)

def get_pattern_atom(self,pattern):

Atom1 =pattern[0]

Atom2 =pattern[1]

mode = int(pattern[-1])

atomlist1 =[]

design = re.findall(r"\d+\.?\d*",Atom1)

if len(design)!= 0:

for item in design:

atomlist1.append(int(item)-1)

else:

for iatom,atom in enumerate(self.atom):

if atom.ele_symbol == Atom1:

atomlist1.append(iatom)

atomlist2 =[]

design = re.findall(r"\d+\.?\d*",Atom2)

if len(design)!= 0:

for item in design:

atomlist2.append(int(item)-1)

else:

atomlist2 =[]

for iatom,atom in enumerate(self.atom):

if atom.ele_symbol == Atom2:

atomlist2.append(iatom)

return atomlist1,atomlist2,mode

def get_all_bond(self):

self.allbond=[]

for i in range(self.natom):

for j in range(i+1,self.natom):

if self.bmx2D[i][j] > 0:

self.allbond.append([i,j,self.bmx2D[i][j]])

self.nbond = len(self.allbond)

'''

def get_atom_info(self):

for i in range(self.natom):

self.atom[i].expbond= 0

self.atom[i].imph= 0

for i in range(self.natom):

for j in range(i+1,self.natom):

if self.bmx2D[i][j] > 0:

if self.atom[j].ele_num !=1:

self.atom[i].expbond = self.atom[i].expbond +1

else:

self.atom[i].imph = self.atom[i].imph +1

if self.atom[i].ele_num !=1:

self.atom[j].expbond = self.atom[j].expbond +1

else:

self.atom[j].imph = self.atom[j].imph +1

'''

def calc_unsaturated_number(self):

nH =0

needH = 2

num_heavyatom = 0

for atom in self.atom:

if atom.ele_num == 1 :

nH = nH +1

else:

num_heavyatom = num_heavyatom +1

needH= needH + (8-atom.ele_num)

self.UnsNum =(needH -nH)/2

self.nheavyatom = num_heavyatom

def remove_metal_bond(self):

self.surface_atom = [0 for i in range(self.natom)]

self.bondneed = [0 for i in range(self.natom)]

self.bmxsave =np.zeros((self.natom,self.natom),dtype =np.int)

for i in range(self.natom):

for j in range(self.natom):

if self.bmx2D[i][j] > 0:

self.bmxsave[i][j] = self.bmx2D[i][j]

if self.atom[i].ele_num > 18 :#or self.atom[j].ele_num > 18:

self.bmx2D[i][j] =0

self.surface_atom[j]= 1

if self.atom[j].ele_num > 18:

self.bmx2D[i][j] =0

self.surface_atom[i]= 1

self.bmx1D = []

for line in self.bmx2D:

self.bmx1D.extend(line)

# def GetNeighbour(self,iatom):

# for i in range(self.natom):

# if self.bmx2D[iatom][i]>0:

# self.atom[iatom].neighbour.append()

#

# def GetAtomnodes(self):

#

# for iatom,atom in enumerate(self):

# self[i]

#

#"=============== the following part is designed for transfer ============="

#"must excute transfer before calling the different part function"

#

def transfer_toP_XYZ_coordStr(self):

'''Get Str: Ele_Name, Ele_Index, Coord, sp, sporder, For, maxF from atom list

Require Latt, natompe, and ele_nameList

'''

self.Cell = self.Latt2Cell()

self.Ele_Name =[]

self.Ele_Index = []

self.Coord = []

for atom in self.atom:

self.Ele_Name.append(atom.ele_symbol)

self.Ele_Index.append(atom.ele_num)

self.Coord.append(atom.xyz)

for iele,elesymbol in enumerate(self.ele_nameList):

self.sp[elesymbol] = self.natompe[iele]

self.sporder[elesymbol] = iele+1

if self.Lfor :

self.For= []

for atom in self.atom:

self.For.append(atom.force)

self.get_max_force()

def TransferToKplstr(self):

'''Get Str.atom from items

Require Latt, For, Coord, Ele_Index

'''

self.atom =[]

for i in range(self.natom):

self.atom.append(S_atom(self.Coord[i],self.Ele_Index[i]))

if self.Lfor:

self.atom[i].force= self.For[i]

self.sort_atom_by_element()

self.set_strinfo_from_atom()

if not self.Cell:

self.Latt2Cell()

#

#"=============== the following part is copyed from zpliu XYZCoord.py ============="

#"must excute TransferToXYZcoordStr before calling the following function"

#

def printCell(self):

""" formated cell print---used to generate INPUT_DEBUG"""

s = ""

for x in self.Cell: s += "%10.4f %10.4f %10.4f\n "%(x[0],x[1],x[2])

return s

def printCoord(self):

""" used to generate INPUT_DEBUG"""

s = ""

for i,xa in enumerate(self.Coord):

s += " %15.9f %15.9f %15.9f %3d\n"%(xa[0],xa[1],xa[2],self.sporder[self.Ele_Name[i]])

return s

def myfilter(self,BadStrP):

if( \

self.energy > BadStrP.HighE*float(self.natom) or \

self.energy < BadStrP.LowE*float(self.natom) or \

min(self.Latt[:3]) < BadStrP.MinLat or \

max(self.Latt[:3]) > BadStrP.MaxLat or \

max(self.Latt[3:7]) > BadStrP.MaxAngle or \

min(self.Latt[3:7]) < BadStrP.MinAngle or \

## max(self.Latt[:3]) > 5.0*min(self.Latt[:3]) or \

self.maxF > BadStrP.MaxFor):

#print '111',self.maxF,BadStrP.MaxFor

return False

else:

#print self.energy,True

# return True

if len(self.For):

# #print self.maxF,BadStrP.MaxFor

if self.maxF > BadStrP.MaxFor :

print(self.maxF, BadStrP.MaxFor)

return False

else: return True

else:

return True

def myfilter_byd(self,Badd1,Badd2,type):

L = True

x1 = [x for x in Badd1.keys()]

y1 = [y for y in Badd2.keys()]

# print 'XXX',x1, y1

for key,value in self.d.items():

# print key, value, Badd1[key], Badd2[key]

if key in x1:

if type==1 and value < Badd1[key] : L = False ; continue

if key in y1:

if type==1 and value > Badd2[key] : L = False ; continue

if key in x1 and key in y1:

if type ==2 and value > Badd1[key] and value < Badd2[key] : L = False

return L

def myfilter_byshortd(self,Badd1):

L = True

for key,value in self.d.items():

if value < Badd1 : L = False ; continue

return L

def Latt2Cell(self):

""" transform of (a,b,c,alpha,beta,gamma) to lattice vector of xyz

tips1: lattice vector is POSCAR format Coord

tips2: Matrix of lattice vector is coord_translation matrix

"""

# lower-triangle

a,b,c,alpha,beta,gamma = self.Latt[:6]

pi = 3.14159265358937932384626

alpha,beta,gamma = alpha*pi/180.0,beta*pi/180.0,gamma*pi/180.0

# bc2 = b**2 + c**2 - 2*b*c*np.cos(alpha)

h1 = a

h2 = b * np.cos(gamma)

h3 = b * np.sin(gamma)

h4 = c * np.cos(beta)

# h5 = ((h2 - h4)**2 + h3**2 + c**2 - h4**2 - bc2)/(2 * h3)

h5 = c * (np.cos(alpha)-np.cos(beta)*np.cos(gamma))/np.sin(gamma)

# modified in 20220422 by JamesBourbon

# x = c**2 - h4**2 - h5**2

h6 = np.sqrt(c**2 - h4**2 - h5**2)

cell = [[h1, 0., 0.], [h2, h3, 0.], [h4, h5, h6]]

return cell

def Cell2Latt(self):

""" transform of lattice vector to (a,b,c,alpha,beta,gamma)