Dear fellow stackoverflow users,
I'm a computational chemist and I have a geometry problem. I have a bunch of coordinates that define a molecular surface and I would like to derive the outward normal vectors of this surface. It seems that the surface approximates the properties of a manifold when I look at it, though the coordinate points were not derived using this framework explicitly. I have to make it clear also that in the general case, molecular surfaces are not always convex hull and can have some degree of concavity. What they do not have is discontinuities, the surfaces are smooth by construction. But since I don't know what to do of these mathematical specifications, I tried to devise an algorithm for the general problem.
As a preliminary remark, for each point of the surface, I can determine the position of the nearest atom. So, for each point, I have these xyz coordinates available as well. The algorithm takes the following form:
1. Computing the distance matrix between each of the available points (which scales, unavoidably, to the square of the number of points but it remains reasonable for my cases using numpy)
2. Extracting the two nearest neighbors of each point
3. Use this triplet of points to generate two vectors centered on each point
4. Get the normal vector based on the cross product of these two vectors followed by its normalization
5. Calculate the vector between the point and its underlying atom
6. If this vector makes an angle below 90° with the normal vector, this vector is inwards pointing and thus it is replaced by its opposite
The result of this full procedure is somewhat okay, but there are still various vectors that, when I visually check the result using matplotlib, are somewhat parallel to the surface. Here is the matplotlib result for the water molecule:
Here is the molecular surface of water for comparison (where you can find the underlying atoms). Ignore the color coding of the surface, it is color coded by surface charge, which is irrelevant for the present discussion.
This surface is obtained by a third party software from which I have no access to the code. I can only visualize it, and I cannot have access to the smoothing procedures used within it for final rendering.
As the image suggests it, the surface is very smooth and thus I expect the normal vectors to account for this "smoothness", which they do imperfectly. I need the normal vectors to actually reflect the smoothness of the surface because the roughness of the surface depicted by the present normal vectors has a noticeable impact on the quality of my subsequent computations based on these normal vectors. Does anyone has any idea about what I could do, under a reasonable computational time, to fix this issue?
Here is a working code that will reproduce my first figure:
import math
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from sklearn.metrics.pairwise import euclidean_distances
coord = np.array([[ -1.2873729481345813 , 0.03256614731449952 , 1.5416924157851974 ],
[ 0.01652948394293976 , 1.163819319621563 , 1.6678642152662158 ],
[ 0.2794741299695759 , -0.5617278297487918 , 2.0029980474538105 ],
[ -0.6610946883884103 , -1.520269012229838 , 0.8599487353786675 ],
[ -1.0054760643749452 , 1.3940480132966795 , 0.33773540172063504 ],
[ 1.4883666878869983 , 0.43621600821755363 , 1.1450836953714032 ],
[ 1.093267571959524 , -1.3195317617899807 , 0.5499808804367919 ],
[ -0.044527698166979345 , -1.2654977063529413 , -0.7625313980846089 ],
[ 0.5831126343857715 , 1.5608391229347571 , -0.025329599172539626 ],
[ -1.0148260960520252 , 0.5209590347086124 , 1.6888058951547953 ],
[ -0.5081521680198725 , 0.9028613364971467 , 1.774471036317154 ],
[ -0.8515308971351676 , -0.19088994302497722 , 1.8836904989342724 ],
[ -0.28882376105739577 , -0.3548423909103548 , 2.05954224229103 ],
[ -1.2492850863979417 , -0.6364567978568106 , 1.3978142132864995 ],
[ -1.041946944653105 , -1.1796586713507826 , 1.095177216946184 ],
[ -1.5436272899575374 , -0.22919366151705667 , 1.1247655324014234 ],
[ -1.4689928122303186 , 0.5000569131417527 , 1.143407587573163 ],
[ -1.290329873679581 , 1.0646014214476844 , 0.8016157211074683 ],
[ 0.14732873180147785 , 0.6686847769257502 , 1.979342311827131 ],
[ 0.2131468214849169 , 0.056951984120299164 , 2.1073021163341092 ],
[ -0.0337308931325195 , -0.9942670995597854 , 1.8046165534409335 ],
[ -0.39947905888129415 , -1.37104542496434 , 1.3601926538530802 ],
[ -1.0634807007597644 , -1.3502061551644402 , 0.46773962277667314 ],
[ -0.7567553825092089 , 1.504835877060738 , 0.749651591341389 ],
[ -0.4472405474525535 , 1.4079282750475794 , 1.2824860152857014 ],
[ 1.5587883382335772 , -0.17395630632391745 , 1.1074452857884236 ],
[ 1.4235308457149791 , -0.8271279096055679 , 0.8993584919303867 ],
[ 0.7692195046037488 , -1.5167605398337778 , 0.1442407045573779 ],
[ 0.32445030395311525 , -1.4651179890494386 , -0.4390957189920336 ],
[ 0.00020957327211999695 , -0.9406835964416262 , -1.0384819480983047 ],
[ -0.025419507383719626 , 1.1326596798290633 , -0.892662493220727 ],
[ 0.273179306851356 , 1.4187473807855993 , -0.5317469722738123 ],
[ 1.0466981946916247 , 1.36444283710828 , 0.43534667343027367 ],
[ 1.3436060983763205 , 0.9793761802108056 , 0.8419277824771877 ],
[ 0.6151824124478109 , 0.9956899926113054 , 1.6618898508556756 ],
[ 1.1365292964079234 , 0.8011459463250883 , 1.4138732980146593 ],
[ 1.2547802693224217 , 0.08060290685487881 , 1.575156048146257 ],
[ 0.8061825016438882 , -0.24826596542943638 , 1.9004575307208522 ],
[ 0.7233645669036894 , -0.8826763855961071 , 1.6883835169082952 ],
[ 0.9497767017034661 , -1.2104420562255824 , 1.1703947277344628 ],
[ 0.5485586943697519 , -1.5967100787262565 , 0.7302017476622693 ],
[ -0.057205092501839166 , -1.6791781005999753 , 0.7696455114375087 ],
[ -0.48640523757119286 , -1.644761057333236 , 0.272618697792796 ],
[ -0.25748236328111623 , -1.501405645515218 , -0.39720000180351417 ],
[ -0.548591180200772 , 1.5966973503597166 , 0.07283122808381894 ],
[ 0.03137820171609954 , 1.6808890342631952 , 0.03811707406017944 ],
[ 0.5162256751875325 , 1.631621931751996 , 0.5739653241581716 ],
[ 0.29520943803145566 , 1.496670085663698 , 1.1960154986767626 ],
[ -0.4357510193390936 , 0.3045200764909755 , 2.0372933983710304 ],
[ -0.7252024085145093 , -0.8680072662231473 , 1.697293877794835 ],
[ -1.4026930025843594 , -0.8680775688445073 , 0.8886609086751069 ],
[ -1.0256665134561849 , 1.0589029487344046 , 1.2875897568848411 ],
[ 0.7343131121864293 , 0.4025541238612941 , 1.9038880538445522 ],
[ 0.2960177800611157 , -1.4085266358073394 , 1.3432377515245404 ],
[ -0.612740309720231 , -1.5270987389589377 , -0.09945502108597855 ],
[ -0.1766515084369774 , 1.6876887908829954 , 0.68242908228023 ],
[ 1.2576856987740817 , -0.6008401432185712 , 1.4092999027282396 ],
[ 0.20077831851211708 , -1.6825464196730353 , 0.10625820785053844 ],
[ -0.27284992140625597 , 1.4578679023663585 , -0.46947226287431315 ],
[ 0.9025274704849868 , 1.2861804462963813 , 1.101223178352364 ],
[ 1.7839350486036138 , -0.019389958495239716 , -1.0076276425199053 ],
[ 1.7377417775538346 , -0.8515596284340875 , -0.06551032812067904 ],
[ 1.9308978238593717 , 0.4526510335304934 , 0.15333098082293778 ],
[ 1.2640412923943216 , 1.1135989051613437 , -0.6487291165436305 ],
[ 0.7110361800083296 , 0.3011211057247756 , -1.4642623430903987 ],
[ 0.9695876990906258 , -0.9216814194628465 , -1.0943999970506841 ],
[ 2.04200929018949 , -0.17973742001943738 , -0.3635312878501747 ],
[ 1.8622593656466528 , 0.5635148250993117 , -0.6107377370213111 ],
[ 1.3955962627400398 , 0.5519618371910718 , -1.1944706185224625 ],
[ 1.2607629333920816 , -0.21551781723753682 , -1.3829616822422597 ],
[ 1.6583970787010358 , -0.6949731322511498 , -0.8399375243260477 ],
[ 1.9444143764920114 , -0.2454332569517364 , 0.2874392721860158 ],
[ 1.636157003293636 , 0.961650115510226 , -0.12326332301091819 ],
[ 0.8313542025004879 , 0.894165587734687 , -1.1420580452198033 ],
[ 0.6430205539552906 , -0.46024492898875324 , -1.4104511496936394 ],
[ 1.2947193164727608 , -1.1400972481068032 , -0.5315305572334722 ],
[ -1.7839430005914738 , 0.019376416779039715 , -1.0076136023161453 ],
[ -1.7377420093346745 , 0.8515508599214876 , -0.06550088225767904 ],
[ -1.9308962617200118 , -0.45265869341099335 , 0.15334861627561777 ],
[ -1.2640463026705615 , -1.1136106280858835 , -0.6487135972817705 ],
[ -0.7110478738279696 , -0.3011369604867556 , -1.4642554706297384 ],
[ -0.9695963617672257 , 0.9216674385272464 , -1.0943972008635638 ],
[ -2.04201196413603 , 0.17972714175629736 , -0.36351594480525473 ],
[ -1.8622640647650528 , -0.5635263554023318 , -0.6107201020978111 ],
[ -1.3956057456456397 , -0.5519763250811119 , -1.1944568661926225 ],
[ -1.2607739609741015 , 0.21550237470677686 , -1.3829529227257198 ],
[ -1.6584036564084357 , 0.6949604403980297 , -0.8399279355844478 ],
[ -1.9444117157749716 , 0.24542627653835639 , 0.2874534822565558 ],
[ -1.636157708161396 , -0.961659176659366 , -0.12324552404161819 ],
[ -0.8313632557119278 , -0.8941798105055468 , -1.1420471832711232 ],
[ -0.6430318069679907 , 0.46022934675447325 , -1.4104486921817194 ],
[ -1.294723366816481 , 1.1400861183930433 , -0.5315262031404322 ],
[ -4.5154929399999335e-06 , -0.6700687207417302 , -1.1844736267236027 ],
[ 0.16109703335223763 , -0.6271001397877708 , -1.2186543868869022 ],
[ -0.16110646810245766 , -0.6271000117262108 , -1.2186531586601221 ],
[ -4.0048342399999414e-06 , 0.6700542836529702 , -1.1844770214133027 ],
[ 0.16109751120177762 , 0.6270854068873908 , -1.218657563554442 ],
[ -0.16110599025291766 , 0.6270855243653509 , -1.2186563353276623 ],
[ 0.36994498278851456 , 0.7707305304675487 , -1.169504657558063 ],
[ 0.7025740208231097 , 1.1378228594549835 , -0.7470950439135691 ],
[ 1.1902461772283626 , 1.1877239069058625 , -0.12779283816255813 ],
[ 1.5961784560313765 , 0.748033999209189 , 0.38770812025235435 ],
[ 1.7530425415210142 , -3.3999814999999503e-06 , 0.5869144868197315 ],
[ 1.5961778861045166 , -0.748041688194589 , 0.3877119097103343 ],
[ 1.1902452718013825 , -1.1877338978242624 , -0.12778682138595815 ],
[ 0.7025731540262697 , -1.1378356163972434 , -0.7470892795558292 ],
[ 0.3699443953987146 , -0.7707451739365088 , -1.1695007532680228 ],
[ 0.27881228877701597 , 1.0004381797559854 , -0.9277686907765464 ],
[ 0.7663313983885688 , 1.3243776625508408 , -0.3086616683059155 ],
[ 1.2752305145620413 , 1.102396662519724 , 0.337597800103595 ],
[ 1.5957376215740167 , 0.42599740483075377 , 0.7446167357487491 ],
[ 1.5957372971866766 , -0.4260032850789138 , 0.7446188937447891 ],
[ 1.2752296742242015 , -1.102404360501184 , 0.3376033845401351 ],
[ 0.7663303887131288 , -1.3243882472092006 , -0.3086549593618755 ],
[ 0.27881152675781595 , -1.0004515288506655 , -0.9277636228196864 ],
[ -0.36995352216617455 , -0.7707449252219087 , -1.169497906808803 ],
[ -0.7025791448730497 , -1.1378352248040433 , -0.747083649080629 ],
[ -1.1902463164027026 , -1.1877331892522427 , -0.12777707971133812 ],
[ -1.5961744570181167 , -0.748040609725749 , 0.3877250868215143 ],
[ -1.7530369449133343 , -2.0638019999999696e-06 , 0.5869289937602514 ],
[ -1.5961738870912567 , 0.748035353909989 , 0.38772129736353433 ],
[ -1.1902454109757226 , 1.1877250123628826 , -0.12778309701711812 ],
[ -0.7025782780762097 , 1.1378235389221032 , -0.747089413438369 ],
[ -0.3699529347763746 , 0.7707308453296488 , -1.1695018110988429 ],
[ -0.27881865163733593 , -1.0004514007891054 , -0.9277613558125664 ],
[ -0.7663327657896889 , -1.3243878630245205 , -0.3086485986182755 ],
[ -1.2752266880614613 , -1.102403584723304 , 0.33761404540041506 ],
[ -1.5957305300328366 , -0.42600214945863374 , 0.7446322698276491 ],
[ -1.5957302056454967 , 0.42599870132175377 , 0.7446301118316091 ],
[ -1.2752258482528014 , 1.102397829890804 , 0.33760846043469506 ],
[ -0.7663317561142488 , 1.3243784467956008 , -0.3086553075623155 ],
[ -0.2788178890889559 , 1.0004384766259653 , -0.9277664242986065 ],
[ -0.23018790924333662 , 0.4635384934491132 , -1.3322307085678806 ],
[ -0.23018827914015663 , -0.20734291124417695 , -1.3622983801385802 ],
[ 0.23017815010577664 , 0.20732757296187695 , -1.3623011720922602 ],
[ 0.23017800828553664 , -0.46355367932757324 , -1.3322301015984204 ]])
atoms = np.asarray([[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -1.21073950301e-06 , 0.0 , 0.40150605663 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ -0.764815694412 , 0.0 , -0.200752550787 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ],
[ 0.764816906209 , 0.0 , -0.200753505843 ]])
#transpose coordinate array
xyz = np.transpose(coord)
#establish distance matrix
n = len(coord)
dist_vectors = list()
dist_vectors.append(xyz[0])
dist_vectors.append(xyz[1])
dist_vectors.append(xyz[2])
dist_vectors = map(list, zip(*dist_vectors))
d = euclidean_distances(dist_vectors, dist_vectors)
#find two nearest neighbors for each segment
x1 = np.zeros((n))
x2 = np.zeros((n))
y1 = np.zeros((n))
y2 = np.zeros((n))
z1 = np.zeros((n))
z2 = np.zeros((n))
for i in range(n):
x_copy = xyz[0]
y_copy = xyz[1]
z_copy = xyz[2]
d1 = np.delete(d[i], i) #removes distance between segment and itself
x_copy = np.delete(x_copy, i)
y_copy = np.delete(y_copy, i)
z_copy = np.delete(z_copy, i)
j1 = np.argmin(d1) #get indice of minimum distance
x1[i] = x_copy[j1]
y1[i] = y_copy[j1]
z1[i] = z_copy[j1]
d2 = np.delete(d1, j1) #removes minimum distance
x_copy = np.delete(x_copy, j1)
y_copy = np.delete(y_copy, j1)
z_copy = np.delete(z_copy, j1)
j2 = np.argmin(d2) #get indice of second minimum distance
x2[i] = x_copy[j2]
y2[i] = y_copy[j2]
z2[i] = z_copy[j2]
#compute normal vector for each segment based on cross product
normal = list()
forGraphs = list()
for i in range(n):
#make vectors for cross product
v1 = np.zeros((3))
v1[0] = x1[i] - coord[i][0]
v1[1] = y1[i] - coord[i][1]
v1[2] = z1[i] - coord[i][2]
v2 = np.zeros((3))
v2[0] = x2[i] - coord[i][0]
v2[1] = y2[i] - coord[i][1]
v2[2] = z2[i] - coord[i][2]
#make cross product and normalize (normal vector should have a unit norm)
nv = np.cross(v1, v2)
nv = nv / np.linalg.norm(nv)
normal.append(nv)
#check of outwards pointing
atv = np.zeros((3))
atv[0] = atoms[i][0] - coord[i][0]
atv[1] = atoms[i][1] - coord[i][1]
atv[2] = atoms[i][2] - coord[i][2]
th_check = math.acos(np.dot(nv, atv) / (np.linalg.norm(nv) * np.linalg.norm(atv)))
if th_check < (math.pi / 2): #if inwards pointing (i. e. pointing towards underlying atom), normal vector is replaced by its opposite
nv[0] = -nv[0]
nv[1] = -nv[1]
nv[2] = -nv[2]
forGraphs.append(np.array([coord[i][0],coord[i][1],coord[i][2],nv[0],nv[1], nv[2]]))
#plot normal vectors (for checkup)
forGraphs = np.asarray(forGraphs)
X, Y, Z, U, V, W = zip(*forGraphs)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.quiver(X, Y, Z, U, V, W)
ax.set_xlim([min(xyz[0])- 1, max(xyz[0]) + 1])
ax.set_ylim([min(xyz[1])- 1, max(xyz[1]) + 1])
ax.set_zlim([min(xyz[2])- 1, max(xyz[2]) + 1])
plt.show()
The first 280 lines of this code are mostly dedicated to the coordinate tables necessary to reproduce the result. The most important part of this code is from line 282 to line 355, where the algorithm I just outlined is implemented.
Thanks in advance for your kind help!
