2019-04-30 12:16:32 +02:00
###########################################################
# #
# __ __ #
# ___( o)> Author : Léo GASPARD <(o )___ #
# \ <_. ) Twitter : @leo_gaspard ( ._> / #
# `---' Mail : leo.gaspard@outlook.fr `---' #
# #
# #
###########################################################
2019-04-17 15:06:53 +02:00
import scipy . optimize as optimize
2019-03-15 11:04:19 +01:00
import numpy as np
import os
import operator
2019-04-17 15:06:53 +02:00
import sys
2019-04-27 16:50:13 +02:00
import datetime
2019-03-15 11:04:19 +01:00
2019-04-30 11:57:26 +02:00
###### GLOBAL VARIABLES #######
2019-03-15 11:04:19 +01:00
dr = 0.001
da = 0.01
DA = 1.5
2019-04-17 15:06:53 +02:00
rBath = ' x '
rPP = ' x '
center = ' x '
xAxis = ' x '
yAxis = ' x '
zAxis = ' x '
sym = ' x '
output_file = ' x '
2019-05-09 08:19:00 +02:00
pattern = [ ]
npattern = [ ]
2019-04-17 15:06:53 +02:00
atoms = ' x '
dist = ' x '
lattice = ' x '
a = ' x '
b = ' x '
c = ' x '
alpha = ' x '
beta = ' x '
gamma = ' x '
visu = 0
seefrag = 0
opti = 0
trsl = ' x '
notIn = ' x '
2019-04-19 09:24:55 +02:00
evj = 0
2019-04-27 16:50:13 +02:00
norep = 0
2019-04-30 11:57:26 +02:00
symop = [ ]
generator = [ ]
2019-05-27 14:31:09 +02:00
noinfrag = [ ]
2019-04-30 11:57:26 +02:00
##############################
def big_cell ( na , nb , nc ) :
coords = [ ]
newCoords = [ ]
newNewCoords = [ ]
for i in generator :
x = i [ 1 ]
y = i [ 2 ]
z = i [ 3 ]
for j in symop :
coords . append ( [ eval ( j [ 0 ] ) , eval ( j [ 1 ] ) , eval ( j [ 2 ] ) , i [ 0 ] ] )
for i in coords :
i [ 0 ] * = a
i [ 1 ] * = b
i [ 2 ] * = c
2019-05-08 14:49:36 +02:00
if i [ 0 ] > a :
i [ 0 ] - = a
if i [ 1 ] > b :
i [ 1 ] - = b
if i [ 2 ] > c :
i [ 2 ] - = c
2019-04-30 11:57:26 +02:00
if i [ 0 ] < 0 :
i [ 0 ] + = a
if i [ 1 ] < 0 :
i [ 1 ] + = b
if i [ 2 ] < 0 :
i [ 2 ] + = c
for i in coords :
if i not in newCoords :
2019-06-26 18:40:47 +02:00
if i [ 0 ] < = a and i [ 1 ] < = b and i [ 2 ] < = c :
2019-04-30 11:57:26 +02:00
newCoords . append ( i )
for i in newCoords :
newNewCoords . append ( i )
for j in range ( 1 , na ) :
newNewCoords . append ( [ i [ 0 ] + a * j , i [ 1 ] , i [ 2 ] , i [ 3 ] ] )
for k in range ( 1 , nb ) :
newNewCoords . append ( [ i [ 0 ] + a * j , i [ 1 ] + b * k , i [ 2 ] , i [ 3 ] ] )
for l in range ( 1 , nc ) :
newNewCoords . append ( [ i [ 0 ] + a * j , i [ 1 ] + b * k , i [ 2 ] + c * l , i [ 3 ] ] )
for k in range ( 1 , nc ) :
newNewCoords . append ( [ i [ 0 ] + a * j , i [ 1 ] , i [ 2 ] + c * k , i [ 3 ] ] )
for j in range ( 1 , nb ) :
newNewCoords . append ( [ i [ 0 ] , i [ 1 ] + b * j , i [ 2 ] , i [ 3 ] ] )
for k in range ( 1 , nc ) :
newNewCoords . append ( [ i [ 0 ] , i [ 1 ] + b * j , i [ 2 ] + c * k , i [ 3 ] ] )
for j in range ( 1 , nc ) :
newNewCoords . append ( [ i [ 0 ] , i [ 1 ] , i [ 2 ] + c * j , i [ 3 ] ] )
return newNewCoords
def printProgressBar ( start , now , iteration , total , prefix = ' ' , suffix = ' ' , decimals = 3 , length = 100 , fill = ' \u2588 ' ) :
2019-04-27 16:50:13 +02:00
"""
Call in a loop to create terminal progress bar
@params :
iteration - Required : current iteration ( Int )
total - Required : total iterations ( Int )
prefix - Optional : prefix string ( Str )
suffix - Optional : suffix string ( Str )
decimals - Optional : positive number of decimals in percent complete ( Int )
length - Optional : character length of bar ( Int )
fill - Optional : bar fill character ( Str )
"""
percent = ( " { 0:. " + str ( decimals ) + " f} " ) . format ( 100 * ( iteration / float ( total ) ) )
filledLength = int ( length * iteration / / total )
bar = fill * filledLength + ' - ' * ( length - filledLength )
2019-04-28 12:57:44 +02:00
dif = now - start
2019-04-27 16:50:13 +02:00
2019-04-30 11:57:26 +02:00
print ( ' \r %30s | %s | %7s %% %s Elapsed time : %s ' % ( prefix , bar , percent , suffix , str ( dif ) ) , end = ' \r ' )
2019-04-27 16:50:13 +02:00
# Print New Line on Complete
if iteration == total :
print ( )
2019-04-17 15:06:53 +02:00
def read_input ( inputFile ) :
global rBath
global rPP
global center
global xAxis
global yAxis
global zAxis
global sym
global output_file
global pattern
global npattern
global atoms
global dist
global lattice
global a
global b
global c
global alpha
global beta
global gamma
global visu
2019-04-19 09:24:55 +02:00
global evj
2019-04-17 15:06:53 +02:00
global seefrag
global opti
global trsl
global notIn
2019-04-27 16:50:13 +02:00
global norep
2019-04-30 11:57:26 +02:00
global symop
global generator
2019-04-17 15:06:53 +02:00
f = open ( inputFile , ' r ' )
line = ' x '
while line != [ ' END_OF_INPUT ' ] :
line = f . readline ( )
line = line . split ( )
if line == [ ] :
continue
elif line [ 0 ] == ' BATH ' :
rBath = float ( line [ 1 ] )
print ( " Bath radius : %f \n " % ( rBath ) )
elif line [ 0 ] == ' PSEUDO ' :
rPP = float ( line [ 1 ] )
print ( " First Shell radius : %f \n " % ( rPP ) )
elif line [ 0 ] == ' CENTER ' :
center = [ ]
for i in range ( 1 , len ( line ) ) :
center . append ( line [ i ] )
elif line [ 0 ] == ' X_AXIS ' :
xAxis = [ ]
for i in range ( 1 , len ( line ) ) :
xAxis . append ( line [ i ] )
elif line [ 0 ] == ' Y_AXIS ' :
yAxis = [ ]
for i in range ( 1 , len ( line ) ) :
yAxis . append ( line [ i ] )
elif line [ 0 ] == ' Z_AXIS ' :
zAxis = [ ]
for i in range ( 1 , len ( line ) ) :
zAxis . append ( line [ i ] )
2019-05-06 13:02:09 +02:00
elif line [ 0 ] == ' SYMMETRY ' :
2019-04-17 15:06:53 +02:00
print ( " Will treat symmetry " )
sym = [ ]
for i in range ( 1 , len ( line ) ) :
sym . append ( line [ i ] )
elif line [ 0 ] == ' OUTPUT ' :
output_file = line [ 1 ]
elif line [ 0 ] == ' PATTERN ' :
pattern = [ ]
2019-05-09 08:19:00 +02:00
line = f . readline ( )
while line . strip ( ) != ' NRETTAP ' :
pattern . append ( [ ] )
line = line . split ( )
for i in range ( len ( line ) ) :
if i % 2 == 0 :
pattern [ - 1 ] . append ( int ( line [ i ] ) )
else :
pattern [ - 1 ] . append ( line [ i ] )
line = f . readline ( )
2019-04-17 15:06:53 +02:00
elif line [ 0 ] == ' NPATTERN ' :
2019-05-09 08:19:00 +02:00
for i in range ( 1 , len ( line ) ) :
npattern . append ( int ( line [ i ] ) )
2019-04-17 15:06:53 +02:00
elif line [ 0 ] == ' LATTICE ' :
a = float ( f . readline ( ) . split ( ) [ 1 ] )
b = float ( f . readline ( ) . split ( ) [ 1 ] )
c = float ( f . readline ( ) . split ( ) [ 1 ] )
alpha = float ( f . readline ( ) . split ( ) [ 1 ] )
beta = float ( f . readline ( ) . split ( ) [ 1 ] )
gamma = float ( f . readline ( ) . split ( ) [ 1 ] )
2019-04-19 09:24:55 +02:00
print ( " Lattice parameter : \n a = %f \n b = %f \n c = %f \n alpha = %f \n beta = %f \n gamma = %f \n " % ( a , b , c , alpha , beta , gamma ) )
2019-04-17 15:06:53 +02:00
elif line [ 0 ] == ' ATOMS ' :
atoms = [ ]
for i in range ( 1 , len ( line ) ) :
if i % 4 == 1 :
atoms . append ( line [ i ] )
elif i % 4 == 2 :
atoms . append ( float ( line [ i ] ) )
elif i % 4 == 3 :
atoms . append ( int ( line [ i ] ) )
elif i % 4 == 0 :
atoms . append ( float ( line [ i ] ) )
elif line [ 0 ] == ' DIST ' :
dist = [ ]
line = f . readline ( )
while line . strip ( ) != ' TSID ' :
line = line . split ( )
dist . append ( [ line [ 0 ] , line [ 1 ] , float ( line [ 2 ] ) ] )
line = f . readline ( )
elif line [ 0 ] == ' COLOR ' :
visu = 1
2019-04-30 11:57:26 +02:00
elif line [ 0 ] == ' NOCOLOR ' :
visu = 2
2019-04-17 15:06:53 +02:00
elif line [ 0 ] == ' TRANSLATE ' :
trsl = [ float ( line [ 1 ] ) , float ( line [ 2 ] ) , float ( line [ 3 ] ) ]
elif line [ 0 ] == ' NOTINPP ' :
notIn = [ ]
for i in range ( 1 , len ( line ) ) :
notIn . append ( line [ i ] )
elif line [ 0 ] == ' OPTIMIZATION ' :
opti = 1
elif line [ 0 ] == ' SEEFRAG ' :
seefrag = 1
2019-04-19 09:24:55 +02:00
elif line [ 0 ] == ' EVJEN ' :
evj = 1
2019-04-27 16:50:13 +02:00
elif line [ 0 ] == ' NOREP ' :
norep = 1
2019-04-30 11:57:26 +02:00
elif line [ 0 ] == ' SYMOP ' :
line = f . readline ( )
while line . strip ( ) != ' POMYS ' :
symop . append ( line . split ( ' , ' ) )
line = f . readline ( )
elif line [ 0 ] == ' GENERATOR ' :
line = f . readline ( )
while line . strip ( ) != ' ROTARENEG ' :
gen = line . split ( )
generator . append ( [ gen [ 0 ] , float ( gen [ 1 ] ) , float ( gen [ 2 ] ) , float ( gen [ 3 ] ) ] )
line = f . readline ( )
2019-05-27 14:31:09 +02:00
elif line [ 0 ] == ' NOINFRAG ' :
line = f . readline ( )
while line . strip ( ) != ' GARFNION ' :
noin = line . split ( )
noinfrag . append ( [ noin [ 0 ] , float ( noin [ 1 ] ) , float ( noin [ 2 ] ) , float ( noin [ 3 ] ) ] )
line = f . readline ( )
2019-03-15 11:04:19 +01:00
def parse ( fileName ) :
f = open ( fileName , ' r ' )
line = 0
dataTable = [ ]
f . readline ( )
f . readline ( )
while True :
line = f . readline ( ) . strip ( ) . split ( )
if line == [ ] :
break
line [ 1 ] = float ( line [ 1 ] )
line [ 2 ] = float ( line [ 2 ] )
line [ 3 ] = float ( line [ 3 ] )
dataTable . append ( line )
return dataTable
2019-04-17 15:06:53 +02:00
def distance ( u , v ) : #Return the distance between two given points
2019-03-15 11:04:19 +01:00
return np . sqrt ( ( u [ 0 ] - v [ 0 ] ) * * 2 + ( u [ 1 ] - v [ 1 ] ) * * 2 + ( u [ 2 ] - v [ 2 ] ) * * 2 )
2019-04-17 15:06:53 +02:00
def vect_product ( u , v ) : #Return the vectorial product between two vectors
2019-03-15 11:04:19 +01:00
return [ u [ 1 ] * v [ 2 ] - u [ 2 ] * v [ 1 ] , u [ 2 ] * v [ 0 ] - u [ 0 ] * v [ 2 ] , u [ 0 ] * v [ 1 ] - u [ 1 ] * v [ 0 ] ]
2019-04-17 15:06:53 +02:00
def dot_product ( u , v ) : #Return the dot product between two vectors
2019-03-15 11:04:19 +01:00
return ( u [ 0 ] * v [ 0 ] + u [ 1 ] * v [ 1 ] + u [ 2 ] * v [ 2 ] )
2019-04-17 15:06:53 +02:00
def normalize ( u ) : #Normalize a vector
2019-03-15 11:04:19 +01:00
norm = np . sqrt ( u [ 0 ] * * 2 + u [ 1 ] * * 2 + u [ 2 ] * * 2 )
2019-04-17 15:06:53 +02:00
if norm != 0 :
return [ u [ 0 ] / norm , u [ 1 ] / norm , u [ 2 ] / norm ]
else :
return u
2019-03-15 11:04:19 +01:00
2019-04-17 15:06:53 +02:00
def rot_matrix ( oldAxis , newAxis ) : #Build a rotation matrix to change director vector
2019-03-15 11:04:19 +01:00
newAxis = normalize ( newAxis )
vp = normalize ( vect_product ( oldAxis , newAxis ) )
angle = np . arccos ( dot_product ( oldAxis , newAxis ) )
rMat = [ [ np . cos ( angle ) + vp [ 0 ] * vp [ 0 ] * ( 1 - np . cos ( angle ) ) , vp [ 0 ] * vp [ 1 ] * ( 1 - np . cos ( angle ) ) - vp [ 2 ] * np . sin ( angle ) , vp [ 0 ] * vp [ 2 ] * ( 1 - np . cos ( angle ) ) + vp [ 1 ] * np . sin ( angle ) ] , [ vp [ 0 ] * vp [ 1 ] * ( 1 - np . cos ( angle ) ) + vp [ 2 ] * np . sin ( angle ) , vp [ 1 ] * vp [ 1 ] * ( 1 - np . cos ( angle ) ) + np . cos ( angle ) , vp [ 1 ] * vp [ 2 ] * ( 1 - np . cos ( angle ) ) - vp [ 0 ] * np . sin ( angle ) ] , [ vp [ 0 ] * vp [ 2 ] * ( 1 - np . cos ( angle ) ) - vp [ 1 ] * np . sin ( angle ) , vp [ 1 ] * vp [ 2 ] * ( 1 - np . cos ( angle ) ) + vp [ 0 ] * np . sin ( angle ) , np . cos ( angle ) + vp [ 2 ] * vp [ 2 ] * ( 1 - np . cos ( angle ) ) ] ]
##This is the formula for the rotation matrix to change axis
return rMat
2019-04-17 15:06:53 +02:00
def rotation ( coord , rMat ) : #Apply a rotation to coordinates
2019-03-15 11:04:19 +01:00
newCoord = [ ]
for i in range ( len ( coord ) ) :
newX = coord [ i ] [ 0 ] * rMat [ 0 ] [ 0 ] + coord [ i ] [ 1 ] * rMat [ 1 ] [ 0 ] + coord [ i ] [ 2 ] * rMat [ 2 ] [ 0 ]
newY = coord [ i ] [ 0 ] * rMat [ 0 ] [ 1 ] + coord [ i ] [ 1 ] * rMat [ 1 ] [ 1 ] + coord [ i ] [ 2 ] * rMat [ 2 ] [ 1 ]
newZ = coord [ i ] [ 0 ] * rMat [ 0 ] [ 2 ] + coord [ i ] [ 1 ] * rMat [ 1 ] [ 2 ] + coord [ i ] [ 2 ] * rMat [ 2 ] [ 2 ]
newCoord . append ( [ newX , newY , newZ ] )
return newCoord
2019-04-17 15:06:53 +02:00
def cut_bath ( rBath , coords ) : #Select which atoms are in the bath with distance constraint (sphere)
2019-03-15 11:04:19 +01:00
bath = [ ]
2019-04-27 16:50:13 +02:00
start = datetime . datetime . now ( )
2019-03-15 11:04:19 +01:00
for i in range ( len ( coords ) ) :
2019-04-27 16:50:13 +02:00
now = datetime . datetime . now ( )
printProgressBar ( start , now , i + 1 , len ( coords ) , prefix = ' Cutting the bath ' , length = 50 )
2019-03-15 11:04:19 +01:00
if distance ( [ 0 , 0 , 0 ] , [ coords [ i ] [ 0 ] , coords [ i ] [ 1 ] , coords [ i ] [ 2 ] ] ) < = rBath :
bath . append ( [ coords [ i ] [ 0 ] , coords [ i ] [ 1 ] , coords [ i ] [ 2 ] , coords [ i ] [ 3 ] , ' C ' ] )
2019-05-06 13:02:09 +02:00
else :
printProgressBar ( start , now , 1 , 1 , prefix = ' Cutting the bath ' , length = 50 )
return bath
2019-03-15 11:04:19 +01:00
2019-04-17 15:06:53 +02:00
def set_pp ( rPP , coords , notIn ) : #Select which atoms are in the first shell of pseudopotential
2019-03-15 11:04:19 +01:00
pp = [ ]
2019-04-27 16:50:13 +02:00
start = datetime . datetime . now ( )
2019-03-15 11:04:19 +01:00
for i in range ( len ( coords ) ) :
2019-04-27 16:50:13 +02:00
now = datetime . datetime . now ( )
printProgressBar ( start , now , i + 1 , len ( coords ) , prefix = ' Finding the first shell ' , length = 50 )
2019-03-15 11:04:19 +01:00
for j in range ( len ( coords ) ) :
2019-05-09 08:19:00 +02:00
if coords [ i ] [ 4 ] != ' O ' or dist_zero ( coords [ j ] ) > 2 * ( max ( [ atoms [ k ] for k in range ( 3 , len ( atoms ) , 4 ) ] ) * sum ( npattern ) + rPP ) :
2019-05-06 13:02:09 +02:00
break
2019-03-15 11:04:19 +01:00
if coords [ i ] [ 4 ] == ' O ' :
2019-03-25 15:12:20 +01:00
if coords [ j ] [ 4 ] == ' C ' and coords [ j ] [ 3 ] not in notIn :
2019-03-15 11:04:19 +01:00
if distance ( [ coords [ i ] [ 0 ] , coords [ i ] [ 1 ] , coords [ i ] [ 2 ] ] , [ coords [ j ] [ 0 ] , coords [ j ] [ 1 ] , coords [ j ] [ 2 ] ] ) < = rPP :
coords [ j ] [ 4 ] == ' Cl '
pp . append ( j )
for i in pp :
coords [ i ] [ 4 ] = ' Cl '
return coords
2019-03-25 15:12:20 +01:00
def find_frag ( pattern , n , coords ) : #We mark the atoms in the bath corresponding to
#the fragment according to user input
inFrag = [ ]
2019-04-27 16:50:13 +02:00
start = datetime . datetime . now ( )
2019-05-09 08:19:00 +02:00
inPattern = [ ]
2019-03-25 15:12:20 +01:00
for k in range ( n ) :
closest = [ 100 , 100 , 100 ]
2019-04-27 16:50:13 +02:00
now = datetime . datetime . now ( )
printProgressBar ( start , now , k + 1 , n , prefix = ' Finding the fragment ' , length = 50 )
2019-03-19 14:52:21 +01:00
for j in coords :
2019-03-25 15:12:20 +01:00
if j [ 3 ] == pattern [ 1 ] :
if distance ( j , [ 0 , 0 , 0 ] ) < distance ( [ 0 , 0 , 0 ] , closest ) and [ j [ 0 ] , j [ 1 ] , j [ 2 ] , distance ( j , j ) , coords . index ( j ) ] not in inFrag :
2019-05-27 14:31:09 +02:00
notin = 0
for noin in noinfrag :
if distance ( j , [ noin [ 1 ] , noin [ 2 ] , noin [ 3 ] ] ) < da :
notin = 1
if notin == 0 :
closest = [ j [ 0 ] , j [ 1 ] , j [ 2 ] , distance ( j , j ) , coords . index ( j ) ]
2019-05-09 08:19:00 +02:00
for i in range ( len ( pattern ) / / 2 ) :
2019-03-25 15:12:20 +01:00
inPattern = [ closest ]
for j in range ( int ( pattern [ 2 * i ] ) ) :
inPattern . append ( [ 100 , 100 , 100 , distance ( [ 100 , 100 , 100 ] , closest ) ] )
for j in coords :
inPattern = sorted ( inPattern , key = operator . itemgetter ( 3 ) )
if j [ 3 ] == pattern [ 2 * i + 1 ] :
if distance ( j , closest ) < = inPattern [ - 1 ] [ 3 ] :
inPattern [ - 1 ] = [ j [ 0 ] , j [ 1 ] , j [ 2 ] , distance ( j , closest ) , coords . index ( j ) ]
for j in inPattern :
inFrag . append ( j )
for j in inFrag :
coords [ j [ 4 ] ] [ 4 ] = ' O '
2019-03-15 11:04:19 +01:00
return coords
2019-04-17 15:06:53 +02:00
def symmetry ( coord , atoms , charges , operations ) : #Find symmetry elements in the coordinates
2019-03-15 11:04:19 +01:00
newCoord = [ ]
2019-04-27 16:50:13 +02:00
total = len ( coord )
start = datetime . datetime . now ( )
2019-04-30 11:57:26 +02:00
progress = 0
2019-05-06 13:02:09 +02:00
coord = sorted ( coord , key = dist_zero )
2019-03-15 11:04:19 +01:00
while coord != [ ] :
toDel = [ ]
newCoord . append ( coord [ 0 ] ) #Add the atom to a new list
name = atoms [ 0 ]
a = newCoord [ - 1 ] [ : ] #label a = E
b = [ - a [ 0 ] , a [ 1 ] , a [ 2 ] ] #label b = yOz mirror plan
c = [ - a [ 0 ] , - a [ 1 ] , a [ 2 ] ] #label c = C2 rotation around z
d = [ a [ 0 ] , - a [ 1 ] , a [ 2 ] ] #label d = xOz mirror plan
e = [ a [ 0 ] , a [ 1 ] , - a [ 2 ] ] #label e = xOy mirror plan
f = [ - a [ 0 ] , a [ 1 ] , - a [ 2 ] ] #label f = C2 rotation around y
g = [ a [ 0 ] , - a [ 1 ] , - a [ 2 ] ] #label g = C2 rotation around x
h = [ - a [ 0 ] , - a [ 1 ] , - a [ 2 ] ] #label h = i
newCoord [ - 1 ] . append ( name + " a " )
newCoord [ - 1 ] . append ( charges [ 0 ] )
2019-04-30 11:57:26 +02:00
progress + = 1
2019-03-15 11:04:19 +01:00
del atoms [ 0 ] #Delete from old list
del coord [ 0 ]
del charges [ 0 ]
2019-04-30 11:57:26 +02:00
now = datetime . datetime . now ( )
printProgressBar ( start , now , progress , total , prefix = ' Treating Symmetry ' , length = 50 , decimals = 3 )
2019-03-15 11:04:19 +01:00
for t in coord :
index = coord . index ( t )
2019-05-06 13:02:09 +02:00
if np . absolute ( dist_zero ( a ) - dist_zero ( t ) ) > da :
break
2019-03-15 11:04:19 +01:00
if name == atoms [ index ] : #Check if it is the same atom
if distance ( t , a ) == 0 :
print ( " ERROR : Twice the same atom " )
if ' xOz ' in operations :
if distance ( t , d ) < = da :
newCoord . append ( d )
newCoord [ - 1 ] . append ( name + ' d ' )
newCoord [ - 1 ] . append ( charges [ index ] )
2019-06-12 16:44:46 +02:00
if index not in toDel :
progress + = 1
toDel . append ( index )
2019-03-15 11:04:19 +01:00
elif da < distance ( t , d ) and distance ( t , d ) < DA :
2019-04-27 16:50:13 +02:00
print ( " Error : This atom should not be there " , distance ( t , d ) , t , d , a , charges [ index ] )
print ( " Are you sure about the xOz symmetry operation ? " )
2019-03-15 11:04:19 +01:00
break
if ' yOz ' in operations :
if distance ( t , b ) < = da :
newCoord . append ( b )
newCoord [ - 1 ] . append ( name + ' b ' )
newCoord [ - 1 ] . append ( charges [ index ] )
2019-06-12 16:44:46 +02:00
if index not in toDel :
progress + = 1
toDel . append ( index )
2019-03-15 11:04:19 +01:00
elif da < distance ( t , b ) and distance ( t , b ) < DA :
2019-04-27 16:50:13 +02:00
print ( " Error : This atom should not be there " , distance ( t , d ) , t , d , a , charges [ index ] )
2019-05-09 14:23:29 +02:00
print ( " Are you sure about the yOz symmetry operation ? " )
2019-03-15 11:04:19 +01:00
break
if ' C2z ' in operations :
if distance ( t , c ) < = da :
newCoord . append ( c )
newCoord [ - 1 ] . append ( name + ' c ' )
newCoord [ - 1 ] . append ( charges [ index ] )
2019-06-12 16:44:46 +02:00
if index not in toDel :
progress + = 1
toDel . append ( index )
2019-03-15 11:04:19 +01:00
elif da < distance ( t , c ) and distance ( t , c ) < DA :
2019-04-27 16:50:13 +02:00
print ( " Error : This atom should not be there " , distance ( t , d ) , t , d , a , charges [ index ] )
2019-05-09 14:23:29 +02:00
print ( " Are you sure about the C2z symmetry operation ? " )
2019-03-15 11:04:19 +01:00
break
if ' xOy ' in operations :
if distance ( t , e ) < = da :
newCoord . append ( e )
newCoord [ - 1 ] . append ( name + ' e ' )
newCoord [ - 1 ] . append ( charges [ index ] )
2019-06-12 16:44:46 +02:00
if index not in toDel :
progress + = 1
toDel . append ( index )
2019-03-15 11:04:19 +01:00
elif da < distance ( t , e ) and distance ( t , e ) < DA :
2019-04-27 16:50:13 +02:00
print ( " Error : This atom should not be there " , distance ( t , d ) , t , d , a , charges [ index ] )
2019-05-09 14:23:29 +02:00
print ( " Are you sure about the xOy symmetry operation ? " )
2019-03-15 11:04:19 +01:00
break
if ' C2y ' in operations :
if distance ( t , f ) < = da :
newCoord . append ( f )
newCoord [ - 1 ] . append ( name + ' f ' )
newCoord [ - 1 ] . append ( charges [ index ] )
2019-06-12 16:44:46 +02:00
if index not in toDel :
progress + = 1
toDel . append ( index )
2019-03-15 11:04:19 +01:00
elif da < distance ( t , f ) and distance ( t , f ) < DA :
2019-04-27 16:50:13 +02:00
print ( " Error : This atom should not be there " , distance ( t , d ) , t , d , a , charges [ index ] )
2019-05-09 14:23:29 +02:00
print ( " Are you sure about the C2y symmetry operation ? " )
2019-03-15 11:04:19 +01:00
break
if ' C2x ' in operations :
if distance ( t , g ) < = da :
newCoord . append ( g )
newCoord [ - 1 ] . append ( name + ' g ' )
newCoord [ - 1 ] . append ( charges [ index ] )
2019-06-12 16:44:46 +02:00
if index not in toDel :
progress + = 1
toDel . append ( index )
2019-03-15 11:04:19 +01:00
elif da < distance ( t , g ) and distance ( t , g ) < DA :
2019-04-27 16:50:13 +02:00
print ( " Error : This atom should not be there " , distance ( t , d ) , t , d , a , charges [ index ] )
2019-05-09 14:23:29 +02:00
print ( " Are you sure about the C2x symmetry operation ? " )
2019-03-15 11:04:19 +01:00
break
if ' i ' in operations :
if distance ( t , h ) < = da :
newCoord . append ( h )
newCoord [ - 1 ] . append ( name + ' h ' )
newCoord [ - 1 ] . append ( charges [ index ] )
2019-06-12 16:44:46 +02:00
if index not in toDel :
progress + = 1
toDel . append ( index )
2019-03-15 11:04:19 +01:00
elif da < distance ( t , h ) and distance ( t , h ) < DA :
2019-04-27 16:50:13 +02:00
print ( " Error : This atom should not be there " , distance ( t , d ) , t , d , a , charges [ index ] )
2019-05-09 14:23:29 +02:00
print ( " Are you sure about the i symmetry operation ? " )
2019-03-15 11:04:19 +01:00
break
2019-05-06 13:02:09 +02:00
now = datetime . datetime . now ( )
printProgressBar ( start , now , progress , total , prefix = ' Treating Symmetry ' , length = 50 , decimals = 3 )
2019-03-15 11:04:19 +01:00
for m in range ( len ( toDel ) ) : #We delete the atoms seen in the simmetry
del coord [ toDel [ m ] - m ]
del atoms [ toDel [ m ] - m ]
del charges [ toDel [ m ] - m ]
return newCoord
def write_input ( fragCoord , ppCoord , bathCoord , fileName , sym ) :
g = open ( fileName , ' w ' )
2019-04-17 15:06:53 +02:00
if sym != ' x ' :
g . write ( ' FRAGMENT \n ' )
g . write ( ' LABEL X Y Z CHARGE \n ' )
2019-03-15 11:04:19 +01:00
for i in range ( len ( fragCoord ) ) :
if fragCoord [ i ] [ 3 ] [ - 1 ] == ' a ' :
2019-05-06 13:02:09 +02:00
g . write ( ' %8s % 7.3f % 7.3f % 7.3f % 8.5f \n ' % ( fragCoord [ i ] [ 3 ] [ : - 1 ] + str ( i + 1 ) , fragCoord [ i ] [ 0 ] , fragCoord [ i ] [ 1 ] , fragCoord [ i ] [ 2 ] , fragCoord [ i ] [ 4 ] ) )
2019-03-15 11:04:19 +01:00
g . write ( " \n \n " )
2019-04-17 15:06:53 +02:00
g . write ( ' PSEUDO \n ' )
g . write ( ' LABEL X Y Z CHARGE \n ' )
2019-03-15 11:04:19 +01:00
for i in range ( len ( ppCoord ) ) :
if ppCoord [ i ] [ 3 ] [ - 1 ] == ' a ' :
2019-05-06 13:02:09 +02:00
g . write ( ' %8s % 7.3f % 7.3f % 7.3f % 8.5f \n ' % ( ppCoord [ i ] [ 3 ] [ : - 1 ] + str ( i + 1 ) , ppCoord [ i ] [ 0 ] , ppCoord [ i ] [ 1 ] , ppCoord [ i ] [ 2 ] , ppCoord [ i ] [ 4 ] ) )
2019-03-15 11:04:19 +01:00
g . write ( " \n \n " )
2019-04-17 15:06:53 +02:00
g . write ( ' CHARGES \n ' )
g . write ( ' LABEL X Y Z CHARGE \n ' )
2019-03-15 11:04:19 +01:00
for i in range ( len ( bathCoord ) ) :
2019-04-27 16:50:13 +02:00
if bathCoord [ i ] [ 3 ] [ - 1 ] == ' a ' :
2019-05-06 13:02:09 +02:00
g . write ( ' %8s % 7.3f % 7.3f % 7.3f % 8.5f \n ' % ( bathCoord [ i ] [ 3 ] [ : - 1 ] + str ( i + 1 ) , bathCoord [ i ] [ 0 ] , bathCoord [ i ] [ 1 ] , bathCoord [ i ] [ 2 ] , bathCoord [ i ] [ 4 ] ) )
2019-04-17 15:06:53 +02:00
if sym == ' x ' :
g . write ( ' FRAGMENT \n ' )
g . write ( ' LABEL X Y Z CHARGE \n ' )
2019-03-15 11:04:19 +01:00
for i in range ( len ( fragCoord ) ) :
2019-04-19 09:24:55 +02:00
g . write ( ' %8s % 7.3f % 7.3f % 7.3f % 8.5f \n ' % ( fragCoord [ i ] [ 3 ] + str ( i + 1 ) , fragCoord [ i ] [ 0 ] , fragCoord [ i ] [ 1 ] , fragCoord [ i ] [ 2 ] , fragCoord [ i ] [ 4 ] ) )
2019-03-15 11:04:19 +01:00
g . write ( ' \n \n ' )
2019-04-17 15:06:53 +02:00
g . write ( ' PSEUDO \n ' )
g . write ( ' LABEL X Y Z CHARGE \n ' )
2019-03-15 11:04:19 +01:00
for i in range ( len ( ppCoord ) ) :
2019-04-19 09:24:55 +02:00
g . write ( ' %8s % 7.3f % 7.3f % 7.3f % 8.5f \n ' % ( ppCoord [ i ] [ 3 ] + str ( i + 1 ) , ppCoord [ i ] [ 0 ] , ppCoord [ i ] [ 1 ] , ppCoord [ i ] [ 2 ] , ppCoord [ i ] [ 4 ] ) )
2019-03-15 11:04:19 +01:00
g . write ( ' \n \n ' )
2019-04-17 15:06:53 +02:00
g . write ( ' CHARGES \n ' )
g . write ( ' LABEL X Y Z CHARGE \n ' )
2019-03-15 11:04:19 +01:00
for i in range ( len ( bathCoord ) ) :
2019-04-19 09:24:55 +02:00
g . write ( ' %8s % 7.3f % 7.3f % 7.3f % 8.5f \n ' % ( bathCoord [ i ] [ 3 ] + str ( i + 1 ) , bathCoord [ i ] [ 0 ] , bathCoord [ i ] [ 1 ] , bathCoord [ i ] [ 2 ] , bathCoord [ i ] [ 4 ] ) )
2019-03-15 11:04:19 +01:00
g . close ( )
2019-04-17 15:06:53 +02:00
def translation ( vec , coords ) :
2019-03-15 11:04:19 +01:00
for i in range ( len ( coords ) ) :
coords [ i ] [ 0 ] - = vec [ 0 ]
coords [ i ] [ 1 ] - = vec [ 1 ]
coords [ i ] [ 2 ] - = vec [ 2 ]
return coords
2019-04-17 15:06:53 +02:00
def get_charge ( charge , numbers , const ) : #Return the square of the charge of the atoms
result = 0
for i in const :
result + = i [ 0 ] * i [ 1 ]
for i in range ( len ( numbers ) ) :
result + = numbers [ i ] * charge [ i ]
return result * * 2
def optimization ( coords ) :
numbers = [ ]
const = [ ]
charge = [ ]
nneighbour = [ ]
for atom in coords :
ini = 0
if atom [ 5 ] == ' full ' :
for i in const :
if i [ 0 ] == atoms [ atoms . index ( atom [ 3 ] ) + 1 ] :
i [ 1 ] + = 1
ini = 1
if ini == 0 :
const . append ( [ atoms [ atoms . index ( atom [ 3 ] ) + 1 ] , 1 ] )
else :
ini = 0
for i in range ( len ( charge ) ) :
if nneighbour [ i ] == atom [ 5 ] and charge [ i ] == atoms [ atoms . index ( atom [ 3 ] ) + 1 ] :
numbers [ i ] + = 1
ini = 1
if ini == 0 :
nneighbour . append ( atom [ 5 ] )
numbers . append ( 1 )
charge . append ( atoms [ atoms . index ( atom [ 3 ] ) + 1 ] )
results = optimize . minimize ( get_charge , charge , args = ( numbers , const ) ) #Scipy built in method thad uses gradient descent to find the local minima of a given function, here it works with the square of the total charge (so that the minima will be at 0)
print ( ' CHARGE OPTIMIZED CHANGE ( % ) ' )
newCharges = results . x
for i in range ( len ( charge ) ) :
print ( ' % 7.5f % 7.5f % 3.2f \n ' % ( charge [ i ] , newCharges [ i ] , 100 - newCharges [ i ] / charge [ i ] * 100 ) )
for atom in coords :
if atom [ 5 ] == ' full ' :
atom [ 5 ] = atoms [ atoms . index ( atom [ 3 ] ) + 1 ]
else :
for i in range ( len ( charge ) ) :
if atom [ 5 ] == nneighbour [ i ] and charge [ i ] == atoms [ atoms . index ( atom [ 3 ] ) + 1 ] :
atom [ 5 ] = newCharges [ i ]
return coords
def count_neighbours ( coords ) :
2019-04-27 16:50:13 +02:00
start = datetime . datetime . now ( )
2019-04-28 12:57:44 +02:00
2019-05-06 13:02:09 +02:00
coords = sorted ( coords , key = get_xyz )
count = 0
2019-04-17 15:06:53 +02:00
for i in coords :
2019-04-28 12:57:44 +02:00
i . append ( 0 )
for i in range ( len ( coords ) - 1 ) :
for j in range ( i + 1 , len ( coords ) ) :
2019-05-06 13:02:09 +02:00
if np . absolute ( coords [ j ] [ 0 ] - coords [ i ] [ 0 ] ) > max ( [ atoms [ i ] for i in range ( 3 , len ( atoms ) , 4 ) ] ) :
count + = len ( coords ) - j
printProgressBar ( start , now , count , ( len ( coords ) * ( len ( coords ) - 1 ) ) / 2 , prefix = ' Counting neighbours ' , length = 50 )
break
elif np . absolute ( coords [ j ] [ 1 ] - coords [ j ] [ 1 ] ) > max ( [ atoms [ i ] for i in range ( 3 , len ( atoms ) , 4 ) ] ) :
count + = len ( coords ) - j
printProgressBar ( start , now , count , ( len ( coords ) * ( len ( coords ) - 1 ) ) / 2 , prefix = ' Counting neighbours ' , length = 50 )
break
elif np . absolute ( coords [ j ] [ 2 ] - coords [ j ] [ 2 ] ) > max ( [ atoms [ i ] for i in range ( 3 , len ( atoms ) , 4 ) ] ) :
count + = len ( coords ) - j
printProgressBar ( start , now , count , ( len ( coords ) * ( len ( coords ) - 1 ) ) / 2 , prefix = ' Counting neighbours ' , length = 50 )
break
elif coords [ i ] [ 5 ] == atoms [ atoms . index ( coords [ i ] [ 3 ] ) + 2 ] :
count + = len ( coords ) - j
printProgressBar ( start , now , count , ( len ( coords ) * ( len ( coords ) - 1 ) ) / 2 , prefix = ' Counting neighbours ' , length = 50 )
break
2019-04-28 12:57:44 +02:00
count + = 1
2019-04-27 16:50:13 +02:00
now = datetime . datetime . now ( )
2019-04-28 12:57:44 +02:00
printProgressBar ( start , now , count , ( len ( coords ) * ( len ( coords ) - 1 ) ) / 2 , prefix = ' Counting neighbours ' , length = 50 )
if distance ( coords [ i ] , coords [ j ] ) < atoms [ atoms . index ( coords [ i ] [ 3 ] ) + 3 ] and coords [ i ] [ 3 ] != coords [ j ] [ 3 ] :
coords [ i ] [ 5 ] + = 1
coords [ j ] [ 5 ] + = 1
2019-05-06 13:02:09 +02:00
2019-04-28 12:57:44 +02:00
for i in coords :
if i [ 5 ] == atoms [ atoms . index ( i [ 3 ] ) + 2 ] :
i [ 5 ] == ' full '
2019-04-17 15:06:53 +02:00
return coords
2019-05-06 13:02:09 +02:00
def get_xyz ( a ) :
return ( a [ 0 ] , a [ 1 ] , a [ 2 ] )
2019-04-19 09:24:55 +02:00
def evjen ( coords ) :
for i in range ( len ( coords ) ) :
if coords [ i ] [ 5 ] == ' full ' :
coords [ i ] [ 5 ] = atoms [ atoms . index ( coords [ i ] [ 3 ] ) + 1 ]
else :
coords [ i ] [ 5 ] = atoms [ atoms . index ( coords [ i ] [ 3 ] ) + 1 ] * coords [ i ] [ 5 ] / atoms [ atoms . index ( coords [ i ] [ 3 ] ) + 2 ]
return coords
2019-04-17 15:06:53 +02:00
2019-05-06 13:02:09 +02:00
def dist_zero ( atom ) :
return np . sqrt ( atom [ 0 ] * * 2 + atom [ 1 ] * * 2 + atom [ 2 ] * * 2 )
2019-03-15 11:04:19 +01:00
def main ( ) :
2019-04-30 11:57:26 +02:00
print ( " Input file is : %s " % ( sys . argv [ 1 ] ) )
2019-04-17 15:06:53 +02:00
read_input ( sys . argv [ 1 ] )
2019-03-15 11:04:19 +01:00
2019-04-30 11:57:26 +02:00
2019-05-29 10:08:44 +02:00
nA = int ( np . floor ( 2 * rBath / a ) + 2 ) #We chose the number of time we need to replicate
nB = int ( np . floor ( 2 * rBath / ( b * np . sin ( np . radians ( gamma ) ) ) ) + 2 ) #to be able to cut the bath
nC = int ( np . floor ( 2 * rBath / ( c * np . sin ( np . radians ( beta ) ) ) ) + 2 )
2019-03-15 11:04:19 +01:00
2019-04-30 11:57:26 +02:00
coords = big_cell ( nA , nB , nC )
2019-03-15 11:04:19 +01:00
2019-06-26 18:40:47 +02:00
2019-03-15 11:04:19 +01:00
coords = translation ( [ nA * a / 2 , nB * b / 2 , nC * c / 2 ] , coords ) #Putting the origin at the center of the cell
2019-06-26 18:40:47 +02:00
2019-03-15 11:04:19 +01:00
2019-03-18 14:36:09 +01:00
2019-03-15 11:04:19 +01:00
labels = [ i [ 3 ] for i in coords ]
2019-04-30 11:57:26 +02:00
2019-04-17 15:06:53 +02:00
if trsl != ' x ' :
2019-04-05 17:55:18 +02:00
translation ( trsl , coords )
2019-03-18 14:36:09 +01:00
centers = [ ] #input, and translating the coordinates
for i in range ( len ( center ) ) :
centers . append ( [ 100 , 100 , 100 ] )
2019-03-15 11:04:19 +01:00
2019-03-18 14:36:09 +01:00
for i in centers :
i . append ( distance ( [ 0 , 0 , 0 ] , i ) )
2019-03-15 11:04:19 +01:00
2019-03-18 14:36:09 +01:00
for i in coords :
centers = sorted ( centers , key = operator . itemgetter ( 3 ) )
if i [ 3 ] in center :
if distance ( i , [ 0 , 0 , 0 ] ) < = centers [ - 1 ] [ 3 ] :
centers [ - 1 ] = [ i [ 0 ] , i [ 1 ] , i [ 2 ] , distance ( i , [ 0 , 0 , 0 ] ) ]
newOgn = np . mean ( np . array ( centers ) , axis = 0 )
newOgn = [ newOgn [ 0 ] , newOgn [ 1 ] , newOgn [ 2 ] ]
coords = translation ( newOgn , coords )
axis = [ ' x ' , ' y ' , ' z ' ]
for k in axis :
2019-04-17 15:06:53 +02:00
if k == ' x ' : #Loop to find the 3 axis
nAxis = xAxis
if k == ' y ' :
nAxis = yAxis
if k == ' z ' :
nAxis = zAxis
2019-03-18 14:36:09 +01:00
nAxiss = [ ] #according to user input, and
#rotating the coordinates
2019-03-19 14:52:21 +01:00
if nAxis [ 0 ] == ' x ' :
continue
2019-03-18 14:36:09 +01:00
for i in range ( len ( nAxis ) ) :
nAxiss . append ( [ 100 , 100 , 100 ] )
for i in nAxiss :
i . append ( distance ( [ 0 , 0 , 0 ] , i ) )
for i in coords :
nAxiss = sorted ( nAxiss , key = operator . itemgetter ( 3 ) )
2019-05-08 14:49:36 +02:00
if i [ 3 ] in nAxis and dist_zero ( i ) > da :
2019-03-18 14:36:09 +01:00
if distance ( i , [ 0 , 0 , 0 ] ) < = nAxiss [ - 1 ] [ 3 ] :
nAxiss [ - 1 ] = [ i [ 0 ] , i [ 1 ] , i [ 2 ] , distance ( i , [ 0 , 0 , 0 ] ) ]
newN = np . mean ( np . array ( nAxiss ) , axis = 0 )
newN = [ newN [ 0 ] , newN [ 1 ] , newN [ 2 ] ]
if k == ' x ' :
oldN = [ 1 , 0 , 0 ]
elif k == ' y ' :
oldN = [ 0 , 1 , 0 ]
elif k == ' z ' :
oldN = [ 0 , 0 , 1 ]
rMat = rot_matrix ( oldN , newN )
coords = rotation ( coords , rMat )
coords = [ [ coords [ i ] [ 0 ] , coords [ i ] [ 1 ] , coords [ i ] [ 2 ] , labels [ i ] ] for i in range ( len ( coords ) ) ]
2019-06-26 18:40:47 +02:00
2019-03-18 14:36:09 +01:00
#We now have one big cell oriented and centered as we want
2019-03-15 11:04:19 +01:00
2019-05-06 13:02:09 +02:00
coords = sorted ( coords , key = dist_zero )
2019-03-15 11:04:19 +01:00
coords = cut_bath ( rBath , coords )
2019-05-09 08:19:00 +02:00
for i in range ( len ( pattern ) ) :
coords = find_frag ( pattern [ i ] , npattern [ i ] , coords )
2019-05-06 13:02:09 +02:00
coords = sorted ( coords , key = operator . itemgetter ( 4 ) )
2019-03-25 15:12:20 +01:00
coords = set_pp ( rPP , coords , notIn )
2019-05-06 13:02:09 +02:00
coords = sorted ( coords , key = dist_zero , reverse = True )
2019-03-25 15:12:20 +01:00
2019-04-19 09:24:55 +02:00
if evj == 1 :
coords = count_neighbours ( coords )
coords = evjen ( coords )
elif opti == 1 :
2019-04-17 15:06:53 +02:00
coords = count_neighbours ( coords )
coords = optimization ( coords )
else :
for i in coords :
i . append ( atoms [ atoms . index ( i [ 3 ] ) + 1 ] )
2019-05-06 13:02:09 +02:00
# coords = sorted(coords,key=operator.itemgetter(5))
# frag = sorted([i for i in coords if i[4] == 'O'],key=operator.itemgetter(3))
# pp = sorted([i for i in coords if i[4] == 'Cl'],key=operator.itemgetter(3))
# bath = sorted([i for i in coords if i[4] == 'C'],key=operator.itemgetter(3))
frag = sorted ( [ i for i in coords if i [ 4 ] == ' O ' ] , key = dist_zero )
pp = sorted ( [ i for i in coords if i [ 4 ] == ' Cl ' ] , key = dist_zero )
bath = sorted ( [ i for i in coords if i [ 4 ] == ' C ' ] , key = dist_zero )
2019-03-15 11:04:19 +01:00
2019-04-17 15:06:53 +02:00
if seefrag == 1 :
2019-03-25 15:12:20 +01:00
g = open ( ' tmp.xyz ' , ' w ' )
g . write ( ' %i \n \n ' % len ( frag ) )
for j in frag :
2019-05-10 10:18:58 +02:00
g . write ( ' %s % 6.5f % 6.5f % 6.5f \n ' % ( j [ 3 ] , j [ 0 ] , j [ 1 ] , j [ 2 ] ) )
2019-03-25 15:12:20 +01:00
g . close ( )
os . system ( ' avogadro tmp.xyz ' )
2019-04-05 17:55:18 +02:00
os . system ( ' rm tmp.xyz ' )
2019-05-02 11:50:56 +02:00
ch = 0
for i in range ( len ( coords ) ) :
ch + = coords [ i ] [ 5 ]
print ( " Total charge : % 8.5f " % ch )
2019-03-25 15:12:20 +01:00
2019-04-17 15:06:53 +02:00
if sym != ' x ' :
2019-04-27 16:50:13 +02:00
if norep == 0 :
rep = 0
2019-04-30 11:57:26 +02:00
prog = 0
2019-04-27 16:50:13 +02:00
start = datetime . datetime . now ( )
for i in range ( len ( coords ) - 1 ) :
for j in range ( i + 1 , len ( coords ) ) :
2019-04-30 11:57:26 +02:00
prog + = 1
2019-04-27 16:50:13 +02:00
now = datetime . datetime . now ( )
2019-04-30 11:57:26 +02:00
printProgressBar ( start , now , prog , ( ( len ( coords ) - 1 ) * len ( coords ) ) / 2 , prefix = ' Calculating nuclear repulsion ' , length = 50 )
2019-04-27 16:50:13 +02:00
rep + = ( coords [ i ] [ 5 ] * coords [ j ] [ 5 ] ) / distance ( coords [ i ] , coords [ j ] )
print ( " Nuclear repulsion before symmetry : %f " % rep )
2019-03-15 11:04:19 +01:00
2019-04-17 15:06:53 +02:00
frag = symmetry ( [ [ i [ 0 ] , i [ 1 ] , i [ 2 ] ] for i in frag ] , [ i [ 3 ] for i in frag ] , [ i [ 5 ] for i in frag ] , sym )
pp = symmetry ( [ [ i [ 0 ] , i [ 1 ] , i [ 2 ] ] for i in pp ] , [ i [ 3 ] for i in pp ] , [ i [ 5 ] for i in pp ] , sym )
bath = symmetry ( [ [ i [ 0 ] , i [ 1 ] , i [ 2 ] ] for i in bath ] , [ i [ 3 ] for i in bath ] , [ i [ 5 ] for i in bath ] , sym )
2019-03-15 11:04:19 +01:00
coords = frag + pp + bath
2019-04-27 16:50:13 +02:00
if norep == 0 :
start = datetime . datetime . now ( )
rep = 0
2019-04-30 11:57:26 +02:00
prog = 0
2019-04-27 16:50:13 +02:00
for i in range ( len ( coords ) - 1 ) :
for j in range ( i + 1 , len ( coords ) ) :
2019-04-30 11:57:26 +02:00
prog + = 1
2019-04-27 16:50:13 +02:00
now = datetime . datetime . now ( )
2019-04-30 11:57:26 +02:00
printProgressBar ( start , now , prog , ( ( len ( coords ) - 1 ) * len ( coords ) ) / 2 , prefix = ' Calculating nuclear repulsion ' , length = 50 )
2019-04-27 16:50:13 +02:00
rep + = ( coords [ i ] [ 4 ] * coords [ j ] [ 4 ] ) / distance ( coords [ i ] , coords [ j ] )
print ( " Nuclear repulsion after symmetry : %f " % rep )
2019-04-17 15:06:53 +02:00
else :
2019-03-15 11:04:19 +01:00
frag = [ [ i [ 0 ] , i [ 1 ] , i [ 2 ] , i [ 3 ] , i [ 5 ] ] for i in frag ]
pp = [ [ i [ 0 ] , i [ 1 ] , i [ 2 ] , i [ 3 ] , i [ 5 ] ] for i in pp ]
bath = [ [ i [ 0 ] , i [ 1 ] , i [ 2 ] , i [ 3 ] , i [ 5 ] ] for i in bath ]
2019-05-06 13:02:09 +02:00
frag = sorted ( frag , key = operator . itemgetter ( 3 ) )
pp = sorted ( pp , key = operator . itemgetter ( 3 ) )
bath = sorted ( bath , key = operator . itemgetter ( 3 ) )
2019-04-17 15:06:53 +02:00
write_input ( frag , pp , bath , output_file , sym )
2019-03-15 11:04:19 +01:00
2019-04-17 15:06:53 +02:00
if visu != 0 :
2019-03-15 11:04:19 +01:00
g = open ( ' tmp.xyz ' , ' w ' )
2019-04-17 15:06:53 +02:00
if visu == 1 :
2019-03-19 14:52:21 +01:00
g . write ( ' %i \n \n ' % ( len ( frag ) + len ( pp ) + len ( bath ) ) )
for i in frag :
2019-04-17 15:06:53 +02:00
g . write ( ' O % 7.3f % 7.3f % 7.3f \n ' % ( i [ 0 ] , i [ 1 ] , i [ 2 ] ) )
2019-03-19 14:52:21 +01:00
for i in pp :
2019-04-17 15:06:53 +02:00
g . write ( ' Cl % 7.3f % 7.3f % 7.3f \n ' % ( i [ 0 ] , i [ 1 ] , i [ 2 ] ) )
2019-03-19 14:52:21 +01:00
for i in bath :
2019-04-17 15:06:53 +02:00
g . write ( ' C % 7.3f % 7.3f % 7.3f \n ' % ( i [ 0 ] , i [ 1 ] , i [ 2 ] ) )
elif visu == 2 :
2019-03-19 14:52:21 +01:00
coords = frag + pp + bath
2019-04-17 15:06:53 +02:00
if sym != ' x ' :
2019-03-19 14:52:21 +01:00
for i in coords :
l = [ ' a ' , ' b ' , ' c ' , ' d ' , ' e ' , ' f ' , ' g ' , ' h ' ]
if i [ 3 ] [ - 1 ] in l :
i [ 3 ] = i [ 3 ] [ : - 1 ]
g . write ( ' %i \n \n ' % len ( coords ) )
for i in coords :
2019-04-17 15:06:53 +02:00
g . write ( ' %s % 7.3f % 7.3f % 7.3f \n ' % ( i [ 3 ] , i [ 0 ] , i [ 1 ] , i [ 2 ] ) )
2019-03-15 11:04:19 +01:00
g . close ( )
os . system ( ' avogadro tmp.xyz ' )
2019-03-25 15:12:20 +01:00
os . system ( ' rm tmp.xyz ' )
2019-03-15 11:04:19 +01:00
2019-05-06 13:02:09 +02:00
2019-03-15 11:04:19 +01:00
main ( )