From 66a287108621ff549a0441345a88575022033439 Mon Sep 17 00:00:00 2001 From: NehZio Date: Thu, 3 Mar 2022 16:37:44 +0100 Subject: [PATCH] Added plane orientation --- crystal_mec.py | 42 ++- src/__pycache__/out.cpython-38.pyc | Bin 0 -> 5308 bytes src/__pycache__/out.cpython-39.pyc | Bin 0 -> 5314 bytes src/__pycache__/read_input.cpython-38.pyc | Bin 0 -> 5107 bytes src/__pycache__/read_input.cpython-39.pyc | Bin 0 -> 5068 bytes src/__pycache__/utils.cpython-38.pyc | Bin 0 -> 8091 bytes src/__pycache__/utils.cpython-39.pyc | Bin 0 -> 8202 bytes src/read_input.py | 26 +- src/read_input.py~ | 196 +++++++++++ src/utils.py | 10 +- src/utils.py~ | 390 ++++++++++++++++++++++ 11 files changed, 658 insertions(+), 6 deletions(-) create mode 100644 src/__pycache__/out.cpython-38.pyc create mode 100644 src/__pycache__/out.cpython-39.pyc create mode 100644 src/__pycache__/read_input.cpython-38.pyc create mode 100644 src/__pycache__/read_input.cpython-39.pyc create mode 100644 src/__pycache__/utils.cpython-38.pyc create mode 100644 src/__pycache__/utils.cpython-39.pyc create mode 100644 src/read_input.py~ create mode 100644 src/utils.py~ diff --git a/crystal_mec.py b/crystal_mec.py index 42c7fc8..e4fcadb 100644 --- a/crystal_mec.py +++ b/crystal_mec.py @@ -31,7 +31,7 @@ if __name__=='__main__': inputFile = sys.argv[1] # Reads all the parameters from the input file - rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translation = read_input(inputFile) + rB , rPP, center, X, Y, Z, xOy, yOz, xOz, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translation = read_input(inputFile) if verbose > 0: out_input_param(rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translation) @@ -79,25 +79,59 @@ if __name__=='__main__': # Orienting the big cell + if xOy != []: + a = find_center(xOy[0], coordinates) + b = a + w = [a] + while np.absolute ( np.absolute( np.dot( a / np.linalg.norm(a) , b / np.linalg.norm(a) ) ) - 1 ) <= 1e-6: + b = find_center(xOy[1], coordinates, without=w) + w.append(b) + c = np.cross(a, b) + k = [0,0,1] + M = rotation_matrix(c, k) + coordinates = rotate(M, coordinates) + if xOz != []: + a = find_center(xOz[0], coordinates) + b = a + w = [a] + while np.absolute ( np.absolute( np.dot( a / np.linalg.norm(a) , b / np.linalg.norm(a) ) ) - 1 ) <= 1e-6: + b = find_center(xOz[1], coordinates, without=w) + w.append(b) + c = np.cross(a, b) + k = [0,1,0] + M = rotation_matrix(c, k) + coordinates = rotate(M, coordinates) + if yOz != []: + a = find_center(yOz[0], coordinates) + b = a + w = [a] + while np.absolute ( np.absolute( np.dot( a / np.linalg.norm(a) , b / np.linalg.norm(a) ) ) - 1 ) <= 1e-6: + b = find_center(yOz[1], coordinates, without=w) + w.append(b) + c = np.cross(a, b) + k = [1,0,0] + M = rotation_matrix(c, k) + coordinates = rotate(M, coordinates) if X != []: k = [1,0,0] xVec = find_center(X,coordinates) - M = rotation_matrix(k,xVec) + M = rotation_matrix(xVec, k) coordinates = rotate(M, coordinates) + if Y != []: k = [0,1,0] yVec = find_center(Y,coordinates) - M = rotation_matrix(k,yVec) + M = rotation_matrix(yVec, k) coordinates = rotate(M, coordinates) if Z != []: k = [0,0,1] zVec = find_center(Z,coordinates) - M = rotation_matrix(k,zVec) + M = rotation_matrix(zVec, k) coordinates = rotate(M, coordinates) diff --git a/src/__pycache__/out.cpython-38.pyc b/src/__pycache__/out.cpython-38.pyc new file mode 100644 index 0000000000000000000000000000000000000000..0caae17558dfc52a3b08c27cf76c9b49a3b602a7 GIT binary patch literal 5308 zcma)AO>Er873T0~f2beZiY)8LN!Zj%wTiuxoJet1$BiWaG^s7SjZ!$<#s;-Rtr&61 z)sS0R1Q#f(0L>+U3-nav;N+H5b16{t)KgLP+!N7TiUK~ir=X{PZ^)%3xo*M*hi~4z z_vX#aH*enTyF){!g6HqQ`;)CqD9V>qIs9m-yoonlMZp!$Jf*?pS8b^Hs-D)?8#+_C z#`O=ChQWt;nO8n!jS{bN6Pz+1;Dh*9ID1E_4L^n*WxdAI3U|ZE_M1*~l34k37T*fq zHN4@Bt|$@uOv+=*0b@j~iBiy7!U|f?+l^$f;FXf0f-Wb+1zkx_6m&H?Sy zqY+Cg2aFHzsG|x$UDWhZ<&2Wo#`u||{C>}on(+@9`t@bM<4HA{oKW}#Ka*3#JJQyc znvYH8L!?Dg==a2d2HkpZ3f@d++++$Pl()&aF!h40kC}+)t2eZ<9>+5C>J^vHltFSw?qd!!|+wkLC(uxWzW75h? z$>bq_)!p-va-aY`EYS+RUqbKOQjf-;>76^K#|u56EAQ`*(3_Te`ZM)DDfEm&@88gS zS?U?j)caeemrU^Q(66QKXb&}EL30K?Lf?c$-vmE|K8&iA_o9irkgdsHN<}ZTM@FUX zF6X^G2iseae?`iVq86cczNpa*U7#9zI9>GcS8(h3ni4lhSB9lJ-Mw z7f8?Ne=N77i*v9upW7ksQ9EO>b57b>cp4FX{J`>8nPsjZvRqt`U%cIoI^D?Hay`d# zL+jA9Y)yw|JiTg1+m^7o+da|>K?qn2OZA2MEr{I-oh}b9bb`q7BiHtf{!A~-OBWA+1w*KPEDIG>eAZlv zg0|bVIzi}W+SyQ0F2^&Hp1zowk>koZ$H~WZa>S$N#pDn!e~-)0f~&s=8rA zvxSqEJD5R7_$dqH8YHc)!1IDV*Kc7?WeaBB^DMj5acq%USPWA=a@$Tkfbqa2iqK3I zoG6;(^LXTVkBebLw=tsljhptSe3o{#&_1U8-ZxM_+Fj%8HCk9<`NV5S z*DeyaA9@mc+4^ldw1l4PaJkPG>#t(=Xcp<46Y=>C;Xo%J?h_>ZZ#>JP0#%k3|CD;~ zh@gdCnDr+dXSeHm?k4c>$k}SjB@Zn-l2y-byHQB%PE#3n2&AFa4XO8QKmPy=ge|lk za5H}EDdQb1WM*Yc*sV5})@lZU;I5C(Ld(q;mts&~p#ES#Hsibh$J8bTre~N5Eovnj zmJK{o<$g2pfK`WQmT!>dCo?v^NI1d1@q-A6gIzy#Iflsl-gRVvdlKKmZ3p7ZaXCj{ zyz~rwX{EpPsTs(~ejSr4%D`M&@%y2}pv{ArKv1k>%KfsN-!3FMEi>k%3me__G=Ypu?l`KvsHGIonzx{9u#qA*s%PL zveMU3C)o>!Budcgzr*Y-YA2-yqF!R+JbXzt+w)S@51IzKmOn%IR`F8I4^5mUk&-aZ z4%D9UD`ksum1`eX66Go~f?g@nxxS-{S#I>ozRpX%a>PDRlM>w@C)!R){P`R9LUQF#72Js^ zFFCQEvb`pqulaPAy1}2H9nwqNhc`*aQ&sl1C zbc)Fa@d{X}dD#QBL|V9h)wm&UeL6?@Ct~&j;KRk&|*F& z%?+HRL&Ax=!iQ!TpcNOuj*V&8GR=54Kl!nx7v{5VE(7?}6;tLFf)F9daIKQ+NJxBO zr<0btzPsCXXwx>MREynCYDjMdAt^PQUJyF9p>*Jo^Nk9nb?Zp9#6<`-jJ<8QiLFeL zP^4@tWnM$W< z5jG5J1RG&i)lm-1Z8OMbR6$qOCJnLw(2aTkFExaJK!M$WTP${iPCJ3?-1tQ689;t1 zq1(W3dgVla#CSO|2){Z0h*OamC$#cjCBdx?H@hm(SjKKx#f~w7+)~aTm{20@q62mY zXqNT}-YA6bHM#DwM(Oq~t4O(klI(fVCk}u24_9(JE_hkG zW{TI~N(Lx`YrHVMlxLw?-cIuU9R3)cCE%t6V#>^#CYH$7_fVu-7>SxG#)$C-G0J(a zkXBuvI}dJRr-<)U1LI;vej|w{@g^!6oG9>Ta2m|ulpR+-*x@(|AW0nrihv+sio$Az zjj|;sZldi-Y@jGw{Qpr$)MzO^7SSCHF^{uxlt9B@dMYrba;$aMazMA=FM3)6?u63Q z@zwA(fCxP?fV>h{fwB_Uc68C;IwCJm$kILF3U~w-DjX;Qu6~&qd7V)`{H$n~b#bp~ z&*bk#+9N$FK_A#1N1WxA9Yx%uK6dXVWigwuM@DWVvs`Lv5@WNwqlo1!DnEiQ^r|1L z7%ToYlDEz(&;ElUFVA7TRs08|l8-cQS)#w1zC}`f)ryR$0N}fbNVW< zL4RQZtNp+zTaH%f5&(n^!vh)L1RVh4brd*FH^a0<*BhUwr7#eY!_!K3aYIxy5y0Cm zhi*b*5c0LEglJ>1@1q1KoeZKz1^tsIrl@!uMWbXln@$H=(av$0PGreOG3fh!gW&-b zv*^ckb(D<(sTF3*m{$URe?rlaomD9=(ofKtQ&Z(3A&-(;H8uQhyR$C`C@7ewRiruf vZsdAlHgrlPWZ0DWmnIH5_+&FxF1G{T^_(kodWJLs=n@;L@`#3ffl zZe$T$pojw0mjEu%Q;~y{TP|%b1&W?}D)I;PMD&)Tg^z84qUfo=H{{ZiTsPr@GjAT> zd-G=IH*emoH85Z(c>esGKUm7RqI^M(!;gl>ckl))D7eCztJImS)w+sWb+wjW*O|gK zu79A^3w(f=c=-cXFY*dEz$x)Q-jBM>*&9lA@CocF>s6MNxf6tz*RUHC#LAv2)MdP@ zc!Tr0qJ-!(E{!M$j1jHIN=|Dr%V|CPUWog1UNIiX=~6tH)8%+5rz`P^oHpW*D^i1nBXuZIP zLl&107$4kGhZTN0Z|TFzsFJlt_~=o7ujf!r`3DUBda~cKxDro{D}0=fX4K%0w6&#X zV^jG6X^|BA9XimUTkB53n~9VgPhy1fdpyYpM*$1c<&%8sLzO9y*e=IA1Nj%F_r?+L z(^k)WYMqhZAL8fUEBqxnYX&@+mEN0QHf!klukc=l-GLqbkt)6iKfWfdsIW33t(=V~ z4*9F@UI>)~1?XY1mg~I_y{}6>8h@&H?wB6W^?!G6`q0Fikgdu-ig_QiN5-V> zUe5Y>0k$zd{<4%GMoXY|A#V|;FH#FVoX&gr9Nb#AqQuRxl_D!Xjb2a56_}%5N&Atu z`|Er*Mlw6PI14-TnH}ODwKD=c=cJv5XYtTO4=n#HwagX7l}qc4NjXv$a) z{g%@(+kW7r+UZbGUWsNTJ$)%PBgd6;j+0O5+=zyajcr@l=JYFpIep17r!T{oM0J9? zW(hkfwK0RX@Ddh8RY;m!zU%sXj@QJT$}gCC*EOwn+qOh%VKGSb&}rFGAI1ZdC;}r< zaGq$k$D>ondt40ax`h!%uidgXZI`@V-L^#2CThLY+O&n~Z<(I$G`BbXj$Fgr)&r;2 zX_=e8$4%Gvn&Eczs_coBu%3&xOIwo63+YPUN{MN6srF6ask#an*@A|Jb=9Zp z%rv8)Q15NwH?bqr{)BDsb{yB)#9pxF%r)fl2c{Lurt7qvFdz)jv<7VgY+!Z*>iznA z?_#mAowj{$M6W+%yp2Uot!xRa*}^iL4c`~s@z7adI@v-~Jn9M59_&X(^w$3|wMjwh zIc5TrzLE{o0ye2}zu~*Ut<6)*H_7tT8Jk`t%wP|BehBozt`|5QL*za0+A`QZjdj7c z4e={cDMMnk^c;(@WeS5g2- zns#+6QN>E4iMw}`Vnc43guRupw-Z)R%Dtse3|~h18<8q*9#*B+DM2GlvpOr zO4w$?8c9LM(?q8jniOQ%gt={hk4}T6VAC3>j2m!1>BCID?xjaSQjrbG?}rCBZI5EA zFOo_#C;JXzJYi14-u+Y&FT=C{EX{5EEql(j{iYS*l;Cq~expM$3FcOPZqKcXeWbwd z#oO0c=K|4~!_aE&{bUIBeaFLDc3++?MHmYnqEJ}B{FSr1*D5oUmGL{zR@e!4j*YQ- zP{f&GgR&lG#V@7K9FizPtM?AFDYQ;V3q+k|;v#%WG|P1p)$9R z^s$255#=aD>j~Ry&{^Ai=*Yz#M|lBnFpQ#!alF&tWig|^exdH_Z=_(7h^nfITALSD z2!6_vf=efw{2(rXl^9oDAiv?a+SfMF{wH1~n%)ymXd_DGajACI^&6HOTm#XQ7juvm zw2)88pADR+1HukF!h>c9fE6_|R+x57!-!_H3;^4CVLsjBGI&4RFl2rq2n~V^SIdcx z#KZ$`+DWnFIlCR3Hf$qIw8&{E1?jCAglJM|xPD+)2a>)+&ezM7-mN3y60b<1y=|v~ zZA?*6By1~TZe7E~h%2cHa=u=cGqvV4Vmc@&Cjuv03w zgQ}~i*dVA=*aWkxj&e|LnSM5-iW`t@&=7kM-K+cX%A4%(Kwu}}E{mO@(=OmTFMOnR z3qXA_raQqeyQNrv%y=m-5Oy={5vMHiO-SY4a*Vqh?spX+v4owlg56>OwZ)9zH?D-( zK?iIUNS5{p+o?_wIw@u36zvR^EN&IDov|Xd2IVd}li^m(!~c5u#t)Zo-oC!RQ6IYd z#`Pbq-MM>b1FhGWZ`Jjc<&Ezmtg}0Hb=6JSs$BO-qqO_BnWtVrN%}nK5r^yj!;OrN za$cIY8Da^pq<|v0MhnwRSvH#H?Slo;P8 zMk&h^l8WPT`@t>j6hZ%iB#etSSwn%OX?;^TQP5A})StpBJ*qsggE169lG+ax0YShN zh14<|W=l-m!M7uEfg))BFGvj$qp5URNOv*BJI=;X0u7&cRbWcxSnKr50o{Q=>uNE$ z<4RXYt)VUe5qewz@`_vq%8Fdu(M6r>h`Kx`OZR{);1O6TbD#vcdOt2?ZAR_jll;4E zi+lO^RQ_J5J=WtQ^nu+m#8_V5QN%s!WA|QM60G%h3j127s_3cp&4Opx`fFMS+uaGf0Yb&GC3r415vVJSnFaI7BrA z0ld|;=|&{ze63a_MC<)MA0;>?7-z}mBo$Ls+(uC^T8)O?Mn<%A9H!%GqLB~!e$QZV z0L3)=FO>c%4>zUk*^vagC=HZhT00Bp98eyp08|9oA)sNPOH6N< zfJT7I%xI4SjRB2=+yu}h&=k-#0}YJU^!5xhUieM@ve2Flh#)L|q_yW*radp_kLd|v z*&V%;tVq|A4s51_{zTdARpFCHM_8Y$%kN02#Q>9_LPg|CqSi`y{TD4J*u_ zB`~!F=1+T=mj=WvDa=2J%5k3MPC!2!Vdbab$z8NqPv{?2vY*Ce$6~TZLUtTvHC~Fz zPQ+wW3E4>{yBd?7ipi!EveQa-J0?35lg%V#XO-;JFcR%-;@%wG)A?MOZ6xi5$BOwNVI$MMc2sh~l$u zIpT+0nTJmR_73>XBSlZ%YbU|mlIVv6st;uS9(@nJHIMB)cs9gFT1W#TM5+R1Zv;+ zi%}oMP+w7~&l0GY5~zO}#8&rt4D*JH>eU2hDTbLE%zyW{G0X*pxu3usNnrAU?EGE~ z^HqiU=LBXsf%&IFFr(_}=)1p2RPL7mKg~PQ{qw6U>>{=ZeU|2nF#nc5Tf42BqHtUY zZST)d>?~Hn&QiQ(T-6NhDDdKHwtnCSI)5Tz+gFwB z55eZ~1pjM_9{_(W!T-AAcjwOWCH4l>PP7-!y(x%1%mMa|6E(rYsP&)dWOhGFAm0My zEO(>mo}YKf>hue;iu>bvS&br@LTxO4)0yO@O?8)T^J%I(jp_WA}8{qAc|s042w&mBt}G8jEXTaE+)jJm=e=sW@`rC z%hn#$x^H7dl5qLqdvD!eC`e<|3btjs>DznElWA`+Xzm3v)3k!X=B~`T(TZAx9@T8 zar4P`-LcK4#r?WlGlOm0RO03qvTvE!e6#MGt{0e=xojZ7S%`wO!*8`GoK?us_kH zw)~o;b#U9Y;=044s)9go>7cleT0_s_@o<96ywA@Z=X?-zT?!~ppY*gU?_n5sM*TYkA#nW(*I&{k^7T4*& z;lneqI-?sCuIdhdgQ>bJwkiCMfp^ic_EeG4lJ3fS&fFp3vlC94T5(E@&_l1?RtB)?NGOG*6RDV+cN>j984n9sk=64 zKgEu3tqlZy)*l51(aHJS-jfYAb8`W=Tpt{@2aFQNC3>drx#qrIg`pb_$6f-4z%dvY z8{t?xQ$QCEu2_ek*YLlFF!V~N>(3UqYBFW-|JZhU5$->yuRN!VdA`P!;!&X+ce0$4tn9Sm= z#I(wkAMkn;PEDsC$kdkOS%J*IhkJ@jv`jTQ^pVv4Lth&90SM_$uPM`36CjyPK?-ZO z)E#hq4!o5aZZ{pP3dmK*^X-KzGN+8Z4||!^_-(24hY$G(iqiKbeJtrZpTetrxbvc% z8V0KSRZUBS!N^zbg#Vs|9l+;t1%J{Iei6A7BN4(bHSOOi0jS@(4 zJ!QBKbtAXGj?0|49#%R+)X)!&-c;zQI}L>TulMwa)*~MYQA>0Jjv+$$Y!NVrfscM* uR2JbJ@-RxMrKs`y5c36%FT+W#V@bmI5Uw(nU2S-5&#~7qLq2XZS^YmfTDq_R literal 0 HcmV?d00001 diff --git a/src/__pycache__/read_input.cpython-39.pyc b/src/__pycache__/read_input.cpython-39.pyc new file mode 100644 index 0000000000000000000000000000000000000000..ac029ebfe32c42453c462dd5203890d252c7a953 GIT binary patch literal 5068 zcmcIoO>7)V74GWp`Ty~F?2Lb7cVL~3!H!K>76@MNDsjTD5UWUu0b8>^)vOR+ao}V)tPnzCytlosU{8BMdq5mHq7vc&65uQ36d}QT)$W<@#16BFozbhV ze)al$uijtvj9DmT75MD_+SsOZit-n{iGEbLc@4-Lhk+`T*vh6tRjR$KY!XWDDvSCF ztSE0*30_V=m*6@LbQ8$CpMfVy7o@vtUFm8-I#2>A36!F0Hx1VeP!=c$ln2=YP!VW^ zYTXjhC{UT|-7%nXpbH>Z0h$1s1e&TVH1R&^PE-BNZ}Jzp?u@T+g{BVG?kr7rFY=2$ za;(tIu6C4Wf9JD3^;qd?$BJCh;X23X>dLOlKG}PZ&(Z9Zg7SRkKw*ETIX*Mstto>! z&F`v5guT=$`3c|_rj$e2{dW{+Nmu;j5oxay1;2iZ7K8N)NWgu%mpKOh1I1BSm8eEK zA`yQSLsTP(KkOrxhD02dh#&H~UY2H$K|3wev3Fq?8)&T-(>^X`{}hqE5RuhmvK5e3 zSvDd&5s^*AWGAKU%MsbBh-@+@J1u3mBeF9Q*;Gt+R?5B;c%qw*tzCpQja7nJ%Q!>4 zWfsiIEch$W^dJM#Baj_32mHKW;B$xMkZdS#t#{}7Je{v7G%L0IgJ-A#d6DxWmV6## zjn6vekRNc7cckw9F4PP2pT+FDB<*Sq_SY%X&qa2ag=k-q`J=V%ezZT&FaFouc>$0U z?2}M4?CR+qpY9X<70Gb$(r=%W;|qYQv*}ny7Wx^fx5xW7UOqn~kBIa0fmZsW#B^eq z#R%rkkeIh4m{%m`Z(^7u5zLQ=!E_!+5ErFC%Q3`KA2HD`4u$xB1n~A{V`_um!#c4rMd1{ zjDHRIDF(eS?1i1nQZH1)A&FRp-CSB#*SgbSNl=f;_P_d;d^xn_IlcTPscmVXjeaFk z=aracFH6}Cu=qlZ|5eHV9`MIw{I5y=M}Y@=BlPQ3Jyy@`a+Bb7sC#<-SoW)@pEyeG zy%R(J1|X-|z0kMk#~dO$bw)({zB?nLfgjWDjlLh(&h}%d9jyC(%(_?LjD9a%joSwN z)Be_3;uFxP(URvx_B6{voBq(Zsh?%>*oIn<(t9@8G}|j3!x_$J;H z?>&YJSm#E(W5aWgIpHeTxXu$i$x}SdGd#<4JkJZf$VYgIkMc4f&08B7TiJui?}Uf5$SUxKTs)bz{LZ8lK^}zF`^%rrn?hnrtwW zHugQ_>{~2>aCiSlxM1%s*S8$MVcIcdk*osOVzJj)bfN}y7sxBYP>0vYx6)n0$xqZh z{jTzX!p4C8fhyFitAaG-ogALLva;hgtrgpH>!t^LrYjq6b-xKpyp@|SwN`Gjw&$Dn z^0#l@T=7_Sg;^#I{CBn07TGns;rUg!*;-!%4gUmEI&%@3HQTM4wzqzDAh?-)v}HbO zc%N<$1DhGX<40~M@OBc0S8*~0!)CJm*=V>i(9rp8G~5a`h;*3H|93yV9_nC8oTXq8 z8dSy!Ceq2j;lt7}dqyu-ip((#fXuOcm%vXIc;yW8Oa^6<2(97x=SJxPJS?Ns`PNVc zf(7X8H@7WjxHTwC*$DbI7*ybbx0$kGrBSsQ{ESnEgx)r3P_RY~>K+xHf?)NV&`m(+ zg2(WuC(FSAy~b|sn9!sw-xP_u*=(9D4+|&gn$EqU_Jf_na@_s;PN3T}>Wu@->6?J< z3_6ZsHyjJJpI~Ko)`xE+=pATw`P21Iv#A)8q!PXm++Z}>lg=C6`E_< z1p5Ifz{a34Yz6J>ndEtHgNju1dBy!pa6?es!Cq?2tc!$o@V4c!0<1kDFFqkltcczk zu8}|Mt^RRggnrANZVI(^uH1lp3C@D+Mzi7d&xc*B{m^C9SYF1{VR?CRROoKYg5#=c zdREP~sYqi_#E^=F=d(r&E-ky^i$u+KOQ|?Frp_1VUQNZHc7W z0!S(n5Wc!4G#fS!e}yO_#jKWXRsp$cyPmapNo1svx8b}JD!U;xcIOT&qbT{VAP)rD zWcaU#$Ofk?%E=}mt6t`_&?&TaSxSTKXY=q7GjPB^PM|?SPThdBk1h0XZ(~%@X>T}p zf+IX5ZzMkt9qpUp;VdJ?UVtYzv(+7|y88y6$%52`j&9miZI4wDm-vzFz+Nx1gYyOJ z&THZ6t?J~{e$(o#npJhRNb)2NUzL>LTb9@28S|`MJ%20%Mww(`ERiCaR4Xvzy$U1h zKx~viis~uDJ?ciC(Roznym?XT2vCE4sN|aDq3$FQ>c8CAA6SokBt$Kd3XB0l@N6EW y(H7hXdSyvIjXd-cYRRh%50OQct-?jlBSFG<4JIR%UTeB^-?r8B|gaUj{OZI3sg1k2sJ8Ifpna=j8(8oGi*k#Cf?Smk}3a zNuEJGCs*Vu;(56y&mvxs>+&4pqC793LA)r%SFFl~Cs>iSRi_n>{yZe;+rSrXA<$Uy zSxXBF*Q}w9yv=!p$UDe8oJWYf8(T+0V>io|b|_r6v}@vH%=wYKr36NO!U0WOu@k$~ zOWdfh;)<7uS|T1KVlNT7SYA=!XRCRn>pqpFe^i^r z$9qqZ(Aam$(ss$BT}gq^*3Pkz_A>U$Id>-KjdIc+wdKCHM!PhHHYG@f!l6$ya7!-L zrWdzuKzdx~=;=DlKALyZZ_VPdx#c;FdZI6nJt;3)RMT*V8J)q|;~X69P`MP)VN_g2 zB2pxf#Hn|q#EF`nN}w`SN}@?bJ@IO--cBuXAB1r&@fx*uyH@pZGz}cvm-wtqFh}hd z8ex3D9=2Nd+qGCV_lhWt=sYd^fqP!OD%OAOdTIsvdY($tzl(?|&&O~J3v5`!e2inz zj&`+&jb#q~VL@AOTKB#;oYQlpn!|a-P|4u}VyNPo!wTNC9*cYaa3yy1iWIAssXMf>gMS^a#@>;y^Lphh z!H4s^So5(MuBD}`Go@$snm&tFucJ1j*J)l-R-~|NO2Hi0Nlqq>an5OK8|d|8tYo;Z z&;4I}qV;*S&FTXB+W(IkpV7}qZ`V>k(yRK+ZU*rgNNSa(^P~Q!aZaDXnHUKnMbf#q zHe7)`^ZE=)L9eVq0>XN1_if1Xg3h!2J{dl%FJPo+Y1as`V}#f%LTH4Mf9&AI0+#Hj z0JBMZ8xtN?6c~*XH+)oUC4Q~f3p+CLRISqp6Q>n+7$O0JO4araQhKU8KpP(j*B+ox z)KzgP6Q|ehRW2skM%W2eE$*snA=+=hJ<3f$RV`u-YKekn3Q7p7cIORzzk+Wg@$22L zl2J9ct!nkJb~kD&P8K}PQrNbpxIy;KaY8$76MDS&XOGn{vj^V z&p9vfv)3J0tcVpKAL=iOd2t!Fm+S!T=Iy{KiWS>+a^fOteik3Kcp2ZdG*9tm$3^QE z^n1Zo*D&7Nzm;X6H0Q;IBbJc#Gx{V=(fWeC`vU~2NMj{7eV#Z z=g_B`F$IyT&@gObFH!0W4d|Q7lWLydZgwQ}oJKl8yNK3d*)GT;>G~S$`4Vz#;)1=>+h%(&{IOX9;hzHS;;_z@FAmNAu94ctI1O+9Us0nw$ph}aHKzqb= zP#R2?W~g)(l<{e;F-GRpsF|7mz8RJKpBfjlJ*7294Ni{Alr^ifsck}0n58}gIVXM# zcG+qqZl|l-iC0(MC<1D$Imvx>~>OC`%ScoE+DXMV4hnLMPOr26Y-~sdi7bzGg@fedLqckbv2@dpt^$O>26EsY<6hJ&sj98+o)n4CA_@Nt)e#x41uMA-^rnl}lC9oG$;{5_ z(lJKxo5-93!ArC^*aKEP0>VSc10V$`BwM1I5j4WL+vZw1vk`_6GGdhS9(ehF)Y}PF zXf|}>Xmj;JeBadnCHkdcP6)6@7sg*jS)JBC`XgsOzNzfNcVHj500Yd^Q+#aSy9rPh z8URa{$#DjEnGBsVYRki5w7+3CW>670d46P z;W569jDaWEw~rf5_7HR8J`>>XW$Wx{f3!;M+}DeA=fbMi zaJMQRL2K!!=L|o8Ec!3<&F*i91@y*RMvjj5_gJWhV05CN8_9$2cB}CDyEf`JV^0@j z{|MLwmzRgjE7BM&&AHS2f}GNwV_(b26=eatg=Xzy4!FT(UEIx3YAneMLnF!MROXN* z>c3C2#Hzlp5t@FoaD&#l$!jfQwN*@e+Awdm z``7!I;rd;)SsE7PJhYzHH;0k!{=1O2dkFusZA1STw$v{`87GdHee4r*NndUJ+s+Sv z|1bXY)7Ng9+x{<6!ktqRB|fWjX&SHrD~0;-{d$Eog|aiSA$W@iEqX;SiPEUcJ7$b)Yw&5{aU>q_CR>U zy?WKHVO(_y^Cw~#3|BR4t&>PJ$4(qnxj~OZe}lFW5fZp@dx-&g9>82gE?{s9go{%8 z2(EpxB+39RBCti9h$@j#AKwz{%b>TWrDf3E<@0o|yj77kjx(!U5LCdfl6)D9dEKxnCpa@-PNeh~!L zf+BcZY7;Z5FHo?CplWZlKxSvpw&I%ubkv`oL8WY@@^Src)klj1cTsmDz% zFr0`G&|+p}vp{XIcZk+jUq<;1<`dM8*-;67VfFTh#^)>bDS~_8~vGbriP|@WZt!};6ie5uv)WeW<;_~L1aGI%4s90tj zH=sD|B;NbvT@!D+)y13N#03{v;c#`1yX`2M8_&s{o%M|k+io?hb_-Y5CgO*P>rH$Q znYd5(s?77^iPaZrZ)Gghg5}QJbpOJ&6n4J^TQ1rK#4kWSmxTH)w3$)P916yk`dwt8 zWS~m0HaImrDUm`tq>ylRc({X`8QJUaOxx?wLu*ndkCh8<6bcAp&Xy~ndKuJ#j;D$h z!~y41#TjOa6vK&sJI+cEPo2nZJ1dw;SaI$se+-M$&Jo-h(RUr3GX&@b?XcqTeakoi zI$%Dx0fFCJN_AEF)pFeJmfM5+PI)Z+sJx9DI+Akd3I?_(zC{E!BkB3um z^d3lHN?yk>K_CS42ML7RCK#eSHlX4W?&NfzPbn?Os7LG_1`HE?S<~J*JX8$r{-u72 zfdcRM?8gZm5|KMie7}YY_;i+Mof
=`kVrQtEC`Ot{(9=y`S@@}1$h4eS5%N-klGtM za|Ai!?E_a=iT_m`!|p(6*HG`G{=IQM|G#Xe9u}wcR6SbE)C2C-cgFQR!%Te!^?yB8 zf3kmo`hOhP)BLmJXHma4RX?+Se(X;X1H||iT#N@&N^T&l{aqiV@p%wMT-A?ILVIS^ zK^mF>`kg-c=V1MWODJU*;vy%;IS;aPwF3eU?c_1&UkoFh`UgC#=FMEu2WAYJV>gz^ z97ydtSjC@W6(6K?Q~6}>g3jQqJ{*#>%;ASC{h=-MOvyh%Sz)S-DDfvIkBbTIDbMSr z)z1;dMo3qTIikWdYoYzl<7&$fx1hpP`(xVw0oL~kt&ye0{|?$F^jp}fIvdyajISbA zn;^SSWd~Jf?>166cUyk{%nTU_&-*>fCT^m literal 0 HcmV?d00001 diff --git a/src/__pycache__/utils.cpython-39.pyc b/src/__pycache__/utils.cpython-39.pyc new file mode 100644 index 0000000000000000000000000000000000000000..9300972f0696a318a254f573a9a6db339cb353ae GIT binary patch literal 8202 zcmb7JO>7*=b?&PEnVvs#I2;ZsN~@h!>_8NDrCoWm-YAN!HDwuAcDaJ=7{(y7Fq&!( zN1Wg49VAx2$ z_qt~|Ln(2>LD%1_-&e1`_f@e`EaoiyeQwiz_~SP$>t~di{gcSt#25ZK0%a-DvFd_j zS(hB!bz4bge`wVm?vPm5oc6R*_4xk~*heLcFBJSFMfn&oCovt4TAQ{JBWbw}~&@LSQiC zGnNq)u32M=yyQGW#!2F1jQP~r(gM9c?_m0SJhs2lX6@XBUi z1y?rpei$`6uf6xB%`2g9UKvI0PI!5+AN#5eL9ZFmqCy{m<(33~QZFFg^r$5LBWV^N zuRTV>VBHnV$cja)QUakg_OVcM8EfU7GnaEGIVC4`xvey)m!?ps0+oDUSCMb#sWkn&fSB^U4Yx7DmNm{t*al-8 z$GBKn=D-^ljCISp|IxT;ie#GO1;jAP@gicF;<#i4JfvhqkLq3&NOdiI9)~13PVf+@h zRanDo!*#4@C9)&OrmYH#?{C+RV7xfM0Oytf}f99BlpNNc~d<@@L~T6 z=6o#1Yia4~TyT#5I*m)liWF8&DICXj(o>Q~KW7cq_09S*W-?wkXaBD) zQTrU~W=#QY<^N;Em&{Ab-L>>D&8k`16Np!!sa4j_FNdE-IkSR2Ni>8MY3KghxC(vd z%?fG3RM(&ZVI9h$gf7pUJnQd=)bw>`_A0llePZ3dC!k|y3O-rJ~_09w`b=3>fD@X_4|C=f53b{LucKb$Bw^?-F$8* z0bjM$=4!CHSNM#*!u$OS@AoSw`>nhdfD=|?5C6mD0eP7j$a}`%N4VGe!^-!*{7+9l z{6Z`?HwuhKu@meyIm5u*;ddCbR+8PxDf7lzdXsEfT~`?9CVq2WeO??;1J)z_eFf0vDfVPwF>LG zZQW>ot-smOQEa1IO`b{<_0R;CbsPQ1Vib$LSbR`-qyDX+6N`2%9>r3>f$wVF)^A*m zrE1`Vi&gb(FZf`Sdyi!!iiNH_`tkmg(=hHoY_}N3&taUfgTNAwU6#J!AL0W2>~jJ? zxo$h6Dykkn$S;WnaT(G}(nq}o>Dwhym5!Yg7a;jfd{pBtd{@&v#g}adwX0}%(b2z! z@#ZKBZ5??j_78WEV5Gl^(=3&6Eu43zET6Y*>pHxKl4%)E=d{dD%W&S3vY}%@Pc$(S zUjK%zFCiW0S^e!^BMQP~ugHPTIJs;D^9VbT4ug}!Y(8d0E=^XYVh2AE4d*91Wf z^t)hy`Yp7nXOemFjHcWhTs^)s&Z1| z8>5+w@FnN6m zx{kdLj%%kCJH5W{#%@#h!w_Jun@nm|Khle+riqPi^dOP{J=&!RPv^f%x5)QU)})zF{whdJuP+@q6s!WbVjz6F zhlle1n?P$}fW34sInLmMlc61kv@-US)s3?e`2k?P0figqm>kiy=R3Im6d3J2B0sxQ z-<#Gfpe9|4T;Ta1;llBAQ#3iy?R(0@CCNtX_3^Z^?JD8lgdBvMC)oc-_Fy6ceg`D25mKlJ~=L0c{J3CJx|>cAZ0^zX@t5hA2<;K$pB# zIR9M<+0Do`rN}!18o~b+;r~k118cGP{EDEVbj?`NGIAxgfOW!}b}@#&iA!8*Pf%*A z!}F;Qms4#+hme0lTEwhAF$j}(YT-J~^Cr)=gtb<%cX^O~qijpq{L+?w8wSXRugTx? zHUAkop|7?6ZRcly`gi~N>5bcoza^{Hbj}i={X3NUT?Efa5R3P5`9C@TbZ}x9zA$*4 zxZ)ZL^i2xhL9kKb;}rXY207wh2zEJ1#@Sv!`bzJ%ZnWxNo8wsY;!G-!{YlLezl@i$ zW z)G@20E3pkpr$ne?S1qv#|VvkWfeoAt@Eg_$JEd@ScEn%8V*#tMIb+VdC%7 z#p5ABvNr=CI&zMg5hYw%*~!Dn{g`+Uy*1+5Bw^9SmOo~t;RF+IK|zTEDea*@M&Dg< zCgQ*2Oj_biOZ>+(_`W|D zE$c7k&;cQu`Y(D~z5M!%o!a#dSaGx89o$$#@BfZ|tS7(0C#4s(TRRW?BOPA8-sv|R zo$v+{liJCvAp;A|YL1~ZMPLe+*~Wz^3VN~o0r}O~-R|`91bE_n3yHD0xEgoExHuh? zxjNe#A8OfY*JTGc*f!!vh?{MEkC?EZtWzhiji+W`!u-M-Eexl2LDF>$w^W?@GR|>H z77$;AX_kfFMV&drB($8K((fY!0|Q0ENt2Pl3lkZmO~wdkhj%ZsaIV=t&RYF%8|x}$Ebzm57EMONwwDc6K` zGlGkAmKVwB;qT3c62PA_uE*U<+gkOiBS-PP-h7-1C&T3jzE&>vKSrY#j^|HPT7L+fe&a zqP|DSp4Q<~Z`AAr4gCO(Mx7A1oOei%3Di4$w%C_KgXfXcd+1}%@6pMh>hDfMiftms z`0y+-_VD3q3gG?(S5qIa=u%%ZnIq_tpF%)^kMRR6uF+)>(r?U2ojEB^o8Fs|dUMhY zq<=jl_2;BnNPjRR&CN;kkp6T=T9}g-Azho1E}TgDwLirOkl`=Fp|~KUSPxFGw@1{m zGW><{j*zESGRh#0od?ZJU;J}8?Sw-pWq08qhl@d)bD=B8*r42{$iwK0e9fTZ@)v8nuI?1IT)Z$2E8W6a@)t2>DJ5)<){P*#{Jqt$$r zm%_Yutzy5j!@Z4=h5 zY}M_}tB=8G$oiSVazLf~C?!gJBBIhKlhUcEZa;pX>T(a$Zm5*%P763yxA(bp+BlID zB+d29@+k#xy|#JoU;p{7^)JkDZ=z)znS@SHo#zZn(xGLixv?;Fa(%q?c7sUo*Q;r& z{UlKQkej?y*2}Zy31Ze4N4@rN6nupUJbBrACei_3Emhq`?e_=z4h=?&sN0R`4JOr; zx*ChYX*ehmaTrCxgUL+p5gOb>Tf7X@54h^s^g2v$z%Trmd*_72lL(kK`C%aOXV8)u z_C>VHiPupI+PIHiH_EZw8+8Z!iATAKBHg2a8AX*+^i0gpt^5MV&no=Fz!!8T;!Hpj zzeM=XAjQ4UERTH$(+U527Y{+5;0BhK@QC!aJMS0#-}kfrihs_}BYwv(`Ro3wzx3bu CnwI7O literal 0 HcmV?d00001 diff --git a/src/read_input.py b/src/read_input.py index 0e0bb37..9dd2f91 100644 --- a/src/read_input.py +++ b/src/read_input.py @@ -9,6 +9,9 @@ def read_input(inputFile): X = [] Y = [] Z = [] + xOy = [] + xOz = [] + yOz = [] symmetry = [] outputFile = "" pattern = [] @@ -67,6 +70,27 @@ def read_input(inputFile): elif ls[0].casefold() == 'z_axis': ls.pop(0) Z = [i for i in ls] + elif ls[0].casefold() == "xoy": + line = f.readline() + ls = line.split() + xOy.append(ls) + line = f.readline() + ls = line.split() + xOy.append(ls) + elif ls[0].casefold() == "xoz": + line = f.readline() + ls = line.split() + xOz.append(ls) + line = f.readline() + ls = line.split() + xOz.append(ls) + elif ls[0].casefold() == "yoz": + line = f.readline() + ls = line.split() + yOz.append(ls) + line = f.readline() + ls = line.split() + yOz.append(ls) elif ls[0].casefold() == 'symmetry': ls.pop(0) symmetry = [i for i in ls] @@ -192,4 +216,4 @@ def read_input(inputFile): print("Bad input : missing the keyword -- %s --"%t) sys.exit() - return rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translate + return rB , rPP, center, X, Y, Z, xOy, yOz, xOz, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translate diff --git a/src/read_input.py~ b/src/read_input.py~ new file mode 100644 index 0000000..1682874 --- /dev/null +++ b/src/read_input.py~ @@ -0,0 +1,196 @@ +import sys + +# Parses the input file +def read_input(inputFile): + + rB = 0.0 + rPP = 0.0 + center = [] + X = [] + Y = [] + Z = [] + symmetry = [] + outputFile = "" + pattern = [] + npattern = [] + atoms = [] + dist = [] + a = 0.0 + b = 0.0 + c = 0.0 + alpha = 90.0 + beta = 90.0 + gamma = 90.0 + showBath = False + evjen = False + showFrag = False + notInPseudo = [] + notInFrag = [] + symGenerator = [] + generator = [] + translate = [0.0,0.0,0.0] + + checkInput = {'bath':False,'pseudo':False,'output':False,'pattern':False,'npattern':False,'a':False,'b':False,'c':False,'atoms':False,'symmetry_generator':False,'generator':False} + + f = open(inputFile,'r') + + line = 'x' + + while line.casefold() != 'end_of_input': + line = f.readline().strip() + ls = line.split() + print(ls) + if ls == []: + continue + if ls[0].casefold() in checkInput: + checkInput[ls[0].casefold()] = True + if ls[0].casefold() == 'bath': + try: + rB = float(ls[1]) + except ValueError: + print("Error while parsing the input file : %s is not a valid bath radius value"%(ls[1])) + sys.exit() + elif ls[0].casefold() == 'pseudo': + try: + rPP = float(ls[1]) + except ValueError: + print("Error while parsing the input file : %s is not a valid pseudopotential radius value"%(ls[1])) + sys.exit() + elif ls[0].casefold() == 'center': + ls.pop(0) + center = [i for i in ls] + elif ls[0].casefold() == 'x_axis': + ls.pop(0) + X = [i for i in ls] + elif ls[0].casefold() == 'y_axis': + ls.pop(0) + Y = [i for i in ls] + elif ls[0].casefold() == 'z_axis': + ls.pop(0) + Z = [i for i in ls] + elif ls[0].casefold() == 'symmetry': + ls.pop(0) + symmetry = [i for i in ls] + elif ls[0].casefold() == 'output': + outputFile = ls[1] + elif ls[0].casefold() == 'pattern': + line = f.readline() + while line.strip().casefold() != 'end': + pattern.append([]) + ls = line.split() + for i in range(len(ls)): + if i%2 == 0: + pattern[-1].append(int(ls[i])) + else: + pattern[-1].append(ls[i]) + line = f.readline() + elif ls[0].casefold() == 'npattern': + ls.pop(0) + try: + npattern = [int(i) for i in ls] + except ValueError: + print("Error while parsing the input file : the number of patterns is not valid %s"%(line)) + sys.exit() + elif ls[0].casefold() == 'lattice': + line = f.readline() + while line.strip().casefold() != 'end': + ls = line.split() + if ls[0].casefold() in checkInput: + checkInput[ls[0].casefold()] = True + if ls[0].casefold() == 'a': + try: + a = float(ls[1]) + except ValueError: + print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1]) + sys.exit() + elif ls[0].casefold() == 'b': + try: + b = float(ls[1]) + except ValueError: + print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1]) + sys.exit() + elif ls[0].casefold() == 'c': + try: + c = float(ls[1]) + except ValueError: + print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1]) + sys.exit() + elif ls[0].casefold() == 'alpha': + try: + alpha = float(ls[1]) + except ValueError: + print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1]) + sys.exit() + elif ls[0].casefold() == 'beta': + try: + beta = float(ls[1]) + except ValueError: + print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1]) + sys.exit() + elif ls[0].casefold() == 'gamma': + try: + gamma = float(ls[1]) + except ValueError: + print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1]) + sys.exit() + line = f.readline() + elif ls[0].casefold() == 'atoms': + line = f.readline() + while line.strip().casefold() != 'end': + ls = line.split() + if(len(ls)) != 4: + print("Error while parsing the input file : not enough values given for the atom in line %s"%line) + sys.exit() + try: + atoms.append([ls[0], float(ls[1]), int(ls[2]), float(ls[3])]) + except ValueError: + print("Error while parsing the input file : bad value for the atom %s"%line) + line = f.readline() + elif ls[0].casefold() == 'show_bath': + showBath = True + elif ls[0].casefold() == 'translate': + ls.pop(0) + try: + translate = [float(ls[0]),float(ls[1]),float(ls[2])] + except ValueError: + print("Error while parsing the input file : the translation vector is not valid %s"%line) + sys.exit() + elif ls[0].casefold() == 'not_in_pseudo': + ls.pop(0) + notInPseudo = [i for i in ls] + elif ls[0].casefold() == 'show_frag': + showFrag = True + elif ls[0].casefold() == 'evjen': + evjen = True + elif ls[0].casefold() == 'symmetry_generator': + line = f.readline().replace("'","") + while line.strip().casefold() != 'end': + symGenerator.append(line.split(',')) + line = f.readline().replace("'","") + elif ls[0].casefold() == 'generator': + line = f.readline() + while line.strip().casefold() != 'end': + ls = line.split() + try: + generator.append([ls[0],float(ls[1]),float(ls[2]),float(ls[3])]) + except ValueError: + print("Error while parsing the input file : bad value for the generator atom %s"%line) + sys.exit() + line = f.readline() + elif ls[0].casefold() == 'not_in_frag': + line = f.readline() + while line.strip().casefold() != 'end': + ls = line.split() + try: + notInFrag.append([float(ls[0]),float(ls[1]),float(ls[2])]) + except ValueError: + print("Error while parsing the input file : bad value for the atom %s"%line) + sys.exit() + line = f.readline() + f.close() + for t in checkInput: + if checkInput[t] == False: + print("Bad input : missing the keyword -- %s --"%t) + sys.exit() + + return rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translate diff --git a/src/utils.py b/src/utils.py index 512f712..0a89a33 100644 --- a/src/utils.py +++ b/src/utils.py @@ -119,7 +119,7 @@ def translate(v,coordinates): # Finds the point at the center of the given atoms that are the # closest to the origin -def find_center(centerList, coordinates): +def find_center(centerList, coordinates, without=[]): centers = [] for i in range(len(centerList)): @@ -129,6 +129,14 @@ def find_center(centerList, coordinates): c.append(distance(c,[0,0,0])) # Computing the distance to the origin for at in coordinates: + w = True + for i in without: + d = distance(at, i) + if d < 1e-6: + w = False + break + if not w: + continue if at[3] in centerList: centers = sorted(centers, key=operator.itemgetter(3)) # Sorting the list with respect to the distance to the origin d = distance(at,[0,0,0]) diff --git a/src/utils.py~ b/src/utils.py~ new file mode 100644 index 0000000..348e604 --- /dev/null +++ b/src/utils.py~ @@ -0,0 +1,390 @@ +import numpy as np +import operator +import sys + +def distance(a,b): + # Returns the 3D distance between a and b where + # a and b are array where x, y and z are at the + # position 0, 1 and 2 + + x = a[0]-b[0] + y = a[1]-b[1] + z = a[2]-b[2] + + return np.sqrt(x**2+y**2+z**2) + +def get_cell_matrix(a,b,c,alpha,beta,gamma): + # Computing the volume of the primitive cell + omega = a*b*c*np.sqrt(1-np.cos(alpha)**2-np.cos(beta)**2-np.cos(gamma)**2+2*np.cos(alpha)*np.cos(beta)*np.cos(gamma)) + + # Computing the matrix + M = [ + [a,b*np.cos(gamma),c*np.cos(beta)], + [0,b*np.sin(gamma),c*(np.cos(alpha)-np.cos(beta)*np.cos(gamma))/(np.sin(gamma))], + [0,0,omega/(a*b*np.sin(gamma))] + ] + return M + +def big_cell(generator,symGenerator,a,b,c,alpha,beta,gamma,nA,nB,nC): + coords = [] + + # Computing the matrix converting fractional to cartesian + fracToCart = get_cell_matrix(a,b,c,alpha,beta,gamma) + + for gen in generator: + x = gen[1] + y = gen[2] + z = gen[3] + + for sym in symGenerator: + u = eval(sym[0]) + v = eval(sym[1]) + w = eval(sym[2]) + + # Making sure the value is within the range [0,1] + u = u + 1*(u<0) - 1*(u>1) + v = v + 1*(v<0) - 1*(v>1) + w = w + 1*(w<0) - 1*(w>1) + + coords.append([u,v,w,gen[0]]) + + # Deleting the redundant atoms + toDel = [] + for i in range(len(coords)-1): + for j in range(i+1,len(coords)): + # Computing the distance using the minimum image convention + # as described in Appendix B equation 9 of + # "Statistical Mechanics : Theory and Molecular Simulations + # Mark E. Tuckerman" + + r1 = np.array(coords[i][:3]) + r2 = np.array(coords[j][:3]) + r12 = r1-r2 + da = np.sqrt(r12[0]**2+r12[1]**2+r12[2]**2) + r12 = r12 - np.round(r12) + db = da - np.sqrt(r12[0]**2+r12[1]**2+r12[2]**2) + r12 = np.matmul(fracToCart,r12) + d = np.sqrt(r12[0]**2+r12[1]**2+r12[2]**2) + + if(d<1e-2): + # We check if we don't already want to delete this atom + if j not in toDel: + toDel.append(j) + + toDel = sorted(toDel) + + # We delete the atoms in the list + for i in range(len(toDel)): + coords.pop(toDel[i]-i) + + newCoords = [] + + # We replicate the cell nA, nB, nC times + for at in coords: + newCoords.append([at[0],at[1],at[2],at[3]]) + for a in range(1,nA): + newCoords.append([at[0]+a,at[1],at[2],at[3]]) + for b in range(1,nB): + newCoords.append([at[0]+a,at[1]+b,at[2],at[3]]) + for c in range(1,nC): + newCoords.append([at[0]+a,at[1]+b,at[2]+c,at[3]]) + for c in range(1,nC): + newCoords.append([at[0]+a,at[1],at[2]+c,at[3]]) + for b in range(1,nB): + newCoords.append([at[0],at[1]+b,at[2],at[3]]) + for c in range(1,nC): + newCoords.append([at[0],at[1]+b,at[2]+c,at[3]]) + for c in range(1,nC): + newCoords.append([at[0],at[1],at[2]+c,at[3]]) + + # Now we convert the fractionnal coordinates to cartesian coordinates + coords = [] + + for at in newCoords: + r = [at[0],at[1],at[2]] + rxyz = np.matmul(fracToCart,r) + coords.append([rxyz[0],rxyz[1],rxyz[2],at[3],'C']) + + + # Returns the list of the atoms [x,y,z,label,second_label] + return coords + +# Translates all the coordinates with the vector v +def translate(v,coordinates): + for c in coordinates: + c[0] += v[0] + c[1] += v[1] + c[2] += v[2] + return coordinates + +# Finds the point at the center of the given atoms that are the +# closest to the origin +def find_center(centerList, coordinates): + + centers = [] + for i in range(len(centerList)): + centers.append([100,100,100]) # Setting a large value for each center + + for c in centers: + c.append(distance(c,[0,0,0])) # Computing the distance to the origin + + for at in coordinates: + if at[3] in centerList: + centers = sorted(centers, key=operator.itemgetter(3)) # Sorting the list with respect to the distance to the origin + d = distance(at,[0,0,0]) + if d <= centers[-1][-1] and d > 0.0: + centers[-1] = [at[0],at[1],at[2],d] + + center = np.mean(centers,axis=0)[:3] # Computing the barycenter + + return center + +# Defines a rotation matrix that will put r1 at the position r2 +def rotation_matrix(r1,r2): + + r1 = np.array(r1)/np.linalg.norm(r1) + r2 = np.array(r2)/np.linalg.norm(r2) + + # Computing the cross product which is the vector around which + # the rotation is done + crossProduct = np.cross(r1,r2) + crossProduct = crossProduct/np.linalg.norm(crossProduct) + + # Computing the angle of the rotation + a = np.arccos(np.dot(r1,r2)) + + c = np.cos(a) + s = np.sin(a) + x = crossProduct[0] + y = crossProduct[1] + z = crossProduct[2] + + M = [ + [x**2*(1-c)+c,x*y*(1-c)-z*s,x*z*(1-c)+y*s], + [x*y*(1-c)+z*s,y**2*(1-c)+c,y*z*(1-c)-x*s], + [x*z*(1-c)-y*s,y*z*(1-c)+x*s,z**2*(1-c)+c] + ] + + return M + +# Rotates all the coordinates using the rotation matric M +def rotate(M,coordinates): + for i in range(len(coordinates)): + r = [coordinates[i][0],coordinates[i][1],coordinates[i][2]] + rV = np.matmul(M,r) + coordinates[i][0] = rV[0] + coordinates[i][1] = rV[1] + coordinates[i][2] = rV[2] + + return coordinates + +# Cuts a sphere centered on the origin in the coordinates +def cut_sphere(coordinates,r): + sphere = [] + for i in range(len(coordinates)): + if distance(coordinates[i],[0,0,0]) <= r: + sphere.append(coordinates[i]) + + return sphere + +# Finds the fragment in the coordinates +def find_fragment(coordinates, patterns, npatterns,notInFrag): + + inFrag = [] + + for n in range(len(patterns)): + pattern = patterns[n] + npattern = npatterns[n] + for i in range(npattern): + c = [100,100,100] + dc = distance([0,0,0],c) + + inPattern = [] + # Finding the closest atom of the first type in the pattern + for at in coordinates: + if at[3] == pattern[1]: + d = distance([0,0,0],at) + if d > 10: + break + if d < dc : + accept = True + for exc in notInFrag: + d = distance(exc,at) + if d < 1e-5: + accept = False + if accept and coordinates.index(at) not in inFrag: + c = [at[0],at[1],at[2],0.0, coordinates.index(at)] + dc = distance([0,0,0],c) + # Finding the rest of the pattern around the atom previously found + atIn = [] + for j in range(0,len(pattern),2): + d = distance(c,[100,100,100]) + # Initializing the atoms + for k in range(pattern[j]): + atIn.append([100,100,100,d]) + + for at in coordinates: + if distance(at,[0,0,0]) > 10: + break + if at[3] == pattern[j+1]: + atIn = sorted(atIn,key=operator.itemgetter(3)) + d = distance(at,c) + trial = [at[0],at[1],at[2],d,coordinates.index(at)] + if d < atIn[-1][3] and trial not in atIn: + accept = True + for exc in notInFrag: + d = distance(exc,trial) + if d < 1e-5: + accept = False + if accept: + atIn[-1] = trial + for at in atIn: + inPattern.append(at[4]) + + for at in inPattern: + if at not in inFrag: + inFrag.append(at) + + for at in inFrag: + coordinates[at][4] = 'O' + + return len(inFrag), coordinates + + +# Finds the pseudopotential layer around +# the fragment +def find_pseudo(coordinates, rPP, notInPseudo): + + for at in coordinates: + if at[4] != 'O': + continue + for i in range(len(coordinates)): + if coordinates[i][4] != 'C': + continue + d = distance(at,coordinates[i]) + if d < rPP: + coordinates[i][4] = 'Cl' + + return coordinates + +# Creates lists containing the neighbours of each +# atom +def find_neighbours(coordinates, atoms): + neighbourList = [[] for i in range(len(coordinates))] + + atoms = np.array(atoms).flatten() + + for i in range(len(coordinates)-1): + for j in range(i+1,len(coordinates)): + li = coordinates[i][3] # Label of the atom i + lj = coordinates[j][3] # Label of the atom j + + ii = np.where(atoms==li)[0] + jj = np.where(atoms==lj)[0] + + ci = float(atoms[ii+1]) # Charge of the atom i + cj = float(atoms[jj+1]) # Charge of the atom j + + if ci*cj < 0: # Checking if the charges have opposite signs + d = distance(coordinates[i],coordinates[j]) + + if d < float(atoms[ii+3]) and d < float(atoms[jj+3]): + neighbourList[i].append(j) + neighbourList[j].append(i) + return neighbourList + +# For each atom, finds if it has the correct number of neighbours, +# if not, modify its charge +def evjen_charges(coordinates,atoms): + neighbourList = find_neighbours(coordinates,atoms) + + atoms = np.array(atoms).flatten() + + charges = [] + + for i in range(len(coordinates)): + li = coordinates[i][3] + ii = np.where(atoms==li)[0] + + nr = len(neighbourList[i]) + nt = int(atoms[ii+2]) + ci = float(atoms[ii+1]) + + if nr > nt: + print("Error : too much neighbours for atom n°%i, count %i neighbours where it should have a maximum of %i"%(i,nr,nt)) + sys.exit() + charges.append(ci*nr/nt) + + return charges + +# Computes the nuclear repulsion +def nuclear_repulsion(coordinates,charges): + + rep = 0.0 + + for i in range(len(coordinates)-1): + for j in range(i+1,len(coordinates)): + rij = distance(coordinates[i],coordinates[j]) + ci = charges[i] + cj = charges[j] + + if(rij < 1): + print(i,j,"\n",coordinates[i],"\n",coordinates[j],"\n",rij) + + rep += (ci*cj)/rij + return rep + +# Computes the symmetry in the whole system +def compute_symmetry(coordinates,charges,symmetry): + symmetrizedCoordinates = [] + symmetrizedCharges = [] + uniqueIndexList = [] # The list containing the indexes of the unique atoms + + treated = [] # Will store the index of the atoms already treated + + symOp = [] + + # Storing all the symmetry operations + for s in symmetry: + if s == 'C2x': + symOp.append(np.array([1,-1,-1])) + elif s == 'C2y': + symOp.append(np.array([-1,1,-1])) + elif s == 'C2z': + symOp.append(np.array([-1,-1,1])) + elif s == 'xOy': + symOp.append(np.array([1,1,-1])) + elif s == 'xOz': + symOp.append(np.array([1,-1,1])) + elif s == 'yOz': + symOp.append(np.array([-1,1,1])) + elif s == 'i': + symOp.append(np.array([-1,-1,-1])) + + for i in range(len(coordinates)): + print(i) + if i in treated: + continue + + treated.append(i) + at1 = np.array(coordinates[i][:3]) + symmetrizedCoordinates.append(coordinates[i]) + symmetrizedCharges.append(charges[i]) + uniqueIndexList.append(len(symmetrizedCoordinates)-1) + + for j in range(len(coordinates)): + if j in treated or coordinates[i][3] != coordinates[j][3]: + continue + + at2 = np.array(coordinates[j][:3]) + + for s in symOp: + if distance(at2, at1*s) < 5: + if distance(at1,at1*s) > 1e-4 and distance(at2,at1*s) < 1e-4: # Checking if op.at1 != at1 and that op.at2 = at1 + p = at1*s + treated.append(j) + symmetrizedCoordinates.append([p[0],p[1],p[2],coordinates[i][3],coordinates[i][4]]) + symmetrizedCharges.append(charges[i]) + break + + return symmetrizedCoordinates,symmetrizedCharges,uniqueIndexList