Dear all, <div><br></div><div>I would like to follow up on this question that I asked on June 10 this year.</div><div>I realize that I might have expressed myself in a very complicated way, and I therefore try to be clearer this time.</div><div><br></div><div>The problem is that when I diagonalize the Hessian printed in the cp2k output file with the following python code (The complete python code is attached)</div><div><br></div><font face="Courier New">from numpy import genfromtxt<br>import numpy as np<br>import sys<br><br>hessianFile=sys.argv[1]<br>hessian = genfromtxt(hessianFile,delimiter=',')<br>freq = np.linalg.eigvals(hessian)</font><div><br></div><div> </div><div>, I get slightly different frequencies than what is printed in the output file (and the .mol file). The frequencies are compared in the table below.</div><div>Column 2: Frequencies from CP2K output file</div><div>Column 3: Frequencies obtained when diagonalizing the Hessian with numpy</div><div><br></div><div><p>Table 1: The eight lowest frequencies. There are 48 in total.</p> <img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAUEAAADrCAYAAAD3/iQNAAAgAElEQVR4Ae2dB7wdRfXHH9JE+GAsCaCEBCRApEsLGEroEAKhkwRIQidKwEIRjECAD6GEkkAgdOm9dxBpoRdBuvSgSFMURCw4/8939Nz/7r7dszv37n13731nPp/39u7u7JTfnP3tzJmZc7qcBUPAEDAEejECXR988IH7/e9/b39NwuC9995z/BnGxWTMsCqGk8lTOE5w3RdffNGN7ru+/OUvu6FDh7oNN9zQ/krGYJ111nFdXV3+j9+GsS5jyCF4mTzqOJkcheODTM0777y+Q5Jkwa4ll1zS/e1vf0tet/MSEPjPf/7j5p9/fveVr3wl9QtUQhYdlcTf//53N2jQIJPHjmrValQG2YLr/vCHP3QrUNd3vvOd1BvdYtqFYAT+9Kc/1Ujwo48+Cn6+tz3wxz/+0Zk89rZW75n6arJlJNjENjASDANXE9SwlCy2IRBHQJOtwiT48MMPuzPPPNNdcskl7v3334/n0MZnf/3rX93dd9/t3n777dJr0RMk+K9//ctdffXV7owzznAXXXSRI08J//73v92vf/3r4LqhPH755Zcd2PRk0ASVcvzlL39xV155pbvgggvcpZde6s466yxfb+rOtWeeeSZW3Ga2bSyjJp2U1bafffaZu/XWW2u4nXfeee7111/3Khow4/yyyy7zR+SFfMsMVWgHTbZySZAX4rHHHnN77bWXW3vttf3w7mc/+5n79NNPy8SpZWkhBN/+9rfdueee69DhlRmaTYL/+Mc/3A033OB22mknt8EGG7hFF13U7bfffo4GRwdy9tlnu4UXXthdfPHFQXX785//7I455phg8mwUO01QSZsZ0d13393XkwmUNdZYw082Ufd+/fq55ZZbzt100021uj777LNNa9tG65r3fJlt+/HHHzveWeQD3Jike/DBBz3ZXXXVVW6VVVbx15dZZhl34oknus8//zyveEH3m/mOFS2IJlu5JEivb+TIkf6FIsPzzz/fHXnkka5TdFzvvPOOO/TQQ91DDz1Ue3mKApsXr5kkyMfphRdecMcdd1ztg3Tddde5RRZZxE2dOtV9+OGHbp999nF9+vRxV1xxReG68SGAPGbOnOnoQfRk0AQ1Wo5TTz3Vrb766p4U5fqdd97pVlppJT8JRX3/+c9/+vvNalvJtxnHZrVtGm6U//nnn3crrrii/4A0oz7NfMeKlleTrVwS/N3vfudWXnlld9tttxXNz+L9D4FmkiBDXUiQ3pEEeuejR492o0aN8gT25JNPusGDB7vLL79couQe6YEw1KSnQB49GTRBlXJAblOmTPG9lxdffFEu++PTTz/tllhiCf9Cv/rqq7F77XTSjLbVcEP1MWTIkKCPZTvhSVk12cokQb5Gv/nNb9xPf/pTP9RgGCI6J14U9Givvfaa19Hw5f3kk088LnyR0dEwvHzjjTdiWNE7Qf8gOpwbb7zRD7W5hq5n9uzZPk16m6TJeTRkpU3v5b777nN8cW6//XbfW33qqaeij/rflFnKxgsj4d577+029MvKi48Cw0zSoecl9Za0oscQEkwrGzjfddddHoc77rjD62SvueaazB4aupdx48b5oQ9DmlmzZrmll17aY8mL9eijj9bKnqX7YShM7/Ldd9+NVqVHfmuCKgXQXmZk9oADDnALLbSQ+9WvfuUfSWtb9GEihxyjsoCMXnjhhR5r5BBVAkNGhpT1hiq0reBGDzrZtpyvtdZa3Uiw6DsQ1R1H349rr722JqvJdtDaoMi7HNoWmmxlkiAvDcS05ppr+iHGCius4HsZt9xyiyfGBRZYwO2www5uu+22c3vssYd/URlSHnTQQW7HHXd0w4YN8/obmXCg0txDr8hiT3RVCy64oDvqqKPcaqut5r70pS/5SZf777/fD3XQXSCABMrywAMPxNImT9JGWT5t2jRP1GPGjPEkQO+H/IVEUfSiGN5555297ozeAsLw29/+1k2fPt2/NKI3y8qLtFhxPmnSJDd8+HCfzkYbbeReeumlzPYoQoIIZ7RsLBEBj3vuuccdfPDBXgdL7476Uif0OieffHLqWjo+Ottss43vxVEosIQEIWvCbrvt5pZffnmPf5ruJzoURqfY00ETVCmLvMzosZI9QUiQj9PXvvY1r9tCLfCtb30rphOlHcGT0c1WW23lsaYdWT9GT/qII47wPUmw5lnWeO6yyy6xHreUJe9YpbalLHzcFl98cXf44Yd7nMCKP85Zn8n7jgyEvgP0JAmozn7+85/794N3fJNNNnGMRng/+TDJ+wzW0TbgXSQu79Lpp5+uvst5mGfd12QrkwQlMSrx3e9+10F+hEceecR9//vf98KBEElP6LnnnvNfk8MOO8wPvzjOOeecnmTQH+69995etygTKgDCC0l6fImXXXZZ/yUiD86jw7i0tOeaay4PLkK9/fbbu/nmm8+/3G+++aYnENKjt0NA57H55pt7vR/nkN/EiRN9rxF951e/+lX/tedeWl5zzz23O+WUU3yvDCW89DLfeustr3vzmaT8K0KClCVaNvL/wQ9+4E477TS/c4KXkDzp5RKYsBg4cKDHLZklvWfuS5vwRQVHhrcQLcKY7AVE02jlUJhyaIIq5dRIkBeY0co3vvEN95Of/MRtscUWnhB5uSWANy89H09IlJ7jYost5kcSyDp6RT4eBJ7jo0LvsJ5QpbalI3D88cd7bOjYRHedcM7EEqsMCKHvgMwP0Pvm4yHvB+8mGMATyXcs2gY/+tGPPDmz8oROFTJP+aLvMp2CRoImW7kkyItEbyKqVxJivP7662vl4gUcMGCAW2+99XwFNttsMz9ryezl448/7nteDO0k0GukV0gPD6HjZSUNAtfIEyHk685wJJk2vU3JX8rD8JpA2vSmyI+vGkNzwGVfajJQNvLW8mL2lbwQjq233trRK+YLCrlqIY8EEcxzzjkns2xSLzCUwAeCjxL5UzcJCBu6sqgKAr3eUkst5b+y9I7y9LoIM2loRCn5NeOoCarkp5Egceipf/3rX/dK/ieeeMJ/XJFdCFICcsCSmn333dd/iJEtZAa8kUmRU2SDD5D0dOT5Iseqta3gRq8rqkemLpwzHOY9Q6bS3rci7wAfFQgv+X6Aq3RqpB3S2oBleMRllj/5LqMOaiRoslUXCaYRI+uMRowY4ehNJAMszvADcpMA8dGj5Br6AqbnAZ8AGEJMNArkmJU28aU8URLl5eec8tAD2nTTTVNJUIaMkGCRvAAT/agMqdGLZoU8EkQw6WHmlQ3digRUAKuuuqo74YQTauu5GKYzm4vKgcALyNCEulFOetzgy0yp9MQlPTkinK+88orXp7ZiKEw5NEGVcsrLnDYc5nk+wtSXF5KPAPWOkiD47bnnnm7jjTf2OlI+DMgacsj20UMOOcRtueWWjplUepP0WGjH0FCltqXsGm4yMcL7h+zkvW/aOxC9Rw+TjzK6adpB3k9GNck2YOTGe0w70B4Sl3PeZUYzQqChbUF8TbZySZCuLT2Pm2++uZY3X0iusSZLAl1hiE4Kz3WGEcxgopCnZ8dXWgKTLtITFNKTZ8mLr7N8DRjWpqVNzwxgpMck5ZHekpzzZUfnEe1R8aJDGgzHaQCJm5UXPS3qQ8+UwLN8sXhZEJy0kEeCPMMXLlk2JjWY3eTjwT30lhLABAKQoT55oDtE3wd+kDmkTy+I9Z18lak/LzvkcNJJJ0lSsSN1oB1Y0xXtYcYiNflEE9Ro1tSPNYLRHiuqCfSmc8wxhx/2EV/kQOQIWQFLXjIZwqGWQTaJC1Ew6YWum5EG+q1GQlXaVuoguCX3z/Ih58OKyoSAbKEWkPeRa8i+9g6AHSoGkR3eD+QNHR/vGL9Fpca1ZBtAkrQB8ifrPck32Ya+gHX802RLJUEq8uMf/9hhaWbs2LG+x4bwTJ482V8bP368f1kpEz0MehpU5thjj/WzaxwZRgIgOsK+fft6wWIox6wza5MYhvBloDuOchQhRZiZGNl///39rBxf6GTaDNtQ+EfLM2HCBA8awkuZGe5QeXpFLPbmxYG0yB+B4OtCXHQQEpe8KGu0HiiU0ZfQY+Uez/NH/ZmtRQDSQhESZFiQLBtlRGAQRnpy9G4QHLBBiQ/+lJMPya677uqxAi/5Q8D4eNGLAQcWUEMYP/zhDx36TeqAgEdnPJnhQ20guse0+jT7miao5M1LRpvRc2b9Ix9V2gFs0CH179/f6/ioF3GZ5BA5oP6QIDO+6L+QLXZAgT0TdAceeKAnAayNMGlE7xEd1S9/+Uu/YoH8+XBCFPyxhpK2YyYZ2RA9bBSjqrQtZeVDILixzpd3mw86pAMW88wzj5clZJy68F4UfQcYqfDh5iPKomx5P3iPUTvQTrxj6LqRL64l24D2pBz8od9Pe5ejH70ozkV+a7KlkiBdUSYOtt12W69kphdB95aXiWvoxxhySeBFQsnJUIOXlZdZAvcYXkB0CCxfHno5AMfXgy8wimwa6uijj/ZECBCiv0imDYkREDQWBVMeiJqXm6Oc01skQNIowSkX+fOCMJ0vdeEZiZuVF70zSIfnqeOMGTN82ln/ipAgzyJ00bL94he/8C8Z2GD5An2mYAo20hul14cOhvtMDvFHDwYsEW6Wy4ADRxZAQ6zcpw0QUMGWMkAc9MAhj1YFTVApE+VlGIWcMGEhskR7IKcs6ZBAXD5sSTmgTfgAI2forZlIohdNepAZayxFRjiyfo5JPZYOgREfEP44R31AecAdvNOGa1VpW8osuDHcj+4YQdUETqx6YEYd+eIjG32X5X3T3gFmd6PvBx8Zep18hEmfe6gpaAM+OtE2AEfO+ROZTb7LeTp4afu0oyZbKgmmJVbWNcCkhyMzSWWlW6V0ipJgVpkZ0qJ2YJ1ZbwiaoPZE/VHzRFUP5ElvnI8MPadkgOBQM9AbhBDl45SMl3be29o2DYOevKbJVktIEIGiV/jNb37TD7e1yYWeBKrsvBohQRoNFQDDuajaoewyVik9TVCbXU50oiyvYThM75BeDMNs1rGyAF90XdFysASEXiIjpCy9cDS+/O6NbSt1b9VRk62WkKAMsxlOM4xhVihtKNEqwMrKtxESpOvP8IGhgagdOhGjKNaaoEbjNes3Oj/02KgeGGIzXI5OppWVb29s27KwqzcdTbZaQoL1VqTdnmuEBNutrmWUVxPUMtK3NHovAppsGQk2US6MBMPA1QQ1LCWLbQjEEdBky0gwjlWpZ0aCYXBqghqWksU2BOIIaLJlJBjHqtSzKpIgC1dZ1tSqXSEawJqgJp9jJpadDixnsmAI5CGgyZaRYB56DdyvEgmy/o8ZT6zQsFMnuWuggWqW9qgmqMlMWmX9OlkOO28PBDTZMhJsYhtWiQRZPMyeZ3aTsCgWoaha0AQ1WlZmyVtl/TpaDvvdPghosmUk2MR2rBIJSjUx0cWypOhuEbnX6qMmqNGytdrkV7Qs9rs9ENBky0iwiW0YQoJp1od5HuMH7I9mSxj7XukBEbiOVeSQHTeQB9Zn2p0EW2n9uoniYkk3EQEjwSaCqyVdhAQxvsCG/KTVawxLsOdy/vnn9/ub2d/KFjr2srKPk72WbHCPWtDWysK9TiDB6FC4ipM7eW1g91uDgJFga3D3G8UhMSxoiOmmZFEwTxS1LM05Rh2wz4jujmcx7ImVXSx84NaArV2cs68VM2BFre52AgnaUDgpQXZeBAEjwSIoNSFOXk+Q/aaaZWnZZC/28DBOycSGnIsF7ajVXWaBsb4h5oyiJrM6gQT5mLTS+nUTxMSS7AEEjAR7AOS0LPJIkKFwEcvSEBmBPdcYAI2ei9VdyV/MTbH3lb+oySzyY7N/u+oEGQq32vq14GzH9kLASLBF7ZVHghQrzfowui4s7YjVa7HqLVZ2k+fSMyxSTSykYBoqy8x+kTSaFUcTVPKsgvXrZtXd0m0uApps2exwE7EvQoJp1oeTVq+ZFGHoi2FMTGslzzEemmd1l2EypqIwc463Ouzm5T3TRGhSk9YElQfYHdJq69epBbeLlUdAky0jwSY2XxESJPuk9WGxei0WeTGlxW4P/N9iWgsryNFzrGLnWd1lmIwpeSxRMxzGR7NY0m4iBEFJa4JKQlWwfh1UIYtcGQQ02TISbGIzFSXBJhahrZLWBLWtKmKFrRwCmmwZCTaxuYwEw8DVBDUsJYttCMQR0GTLSDCOValnRoJhcGqCGpaSxTYE4ghosmUkGMeq1DMjwTA4NUENS8liGwJxBDTZMhKMY1XqmZFgGJyaoIalZLENgTgCmmwZCcaxKvXMSDAMTk1Qw1Ky2IZAHAFNtowE41iVemYkGAanJqhhKVlsQyCOgCZbRoJxrEo9MxIMg1MT1LCUyotdZXcE5dWy81PSZMtIsIntbyQYBq4mqGEpNR67HdwRNF7L3pOCJltGgk2UAyPBMHA1QQ1LqfHY7eCOoPFa9p4UNNkyEmyiHBgJhoGrCWpYSuXFrrI7gvJq2fkpabJlJNjE9jcSDANXE9RkSuaOIImInWsIaLJlJKgh1+A9I8EwADVBlZTMHYEgYccQBDTZMhIMQTIwrpFgGGCaoEpK5o5AkLBjCAKabBkJhiAZGNdIMAwwTVBJqWx3BNgnzHJFQH5Vd0cQhm7vjq3JlpFgE2XDSDAMXE1QSalsdwQYlcVrX5orAvKDBKvsjiAM3d4dW5MtI8EmyoaRYBi4mqBKSuaOQJCwYwgCmmwZCYYgGRjXSDAMME1QJSVzRyBI2DEEAU22jARDkAyMayQYBpgmqNGUzB1BFA37XQQBTbaMBIsgWGccI8Ew4DRBDUvJYhsCcQQ02epacsklHS4eLTQHgQUWWMDNP//8Dp+5FnQEmIjAj7LJo46T3Q1HANmC69gOmQxdiy66qJsyZYo744wz7K9kDKZOnermnXdeN88887gTTzzR8M3B9/jjj3d9+/Y1eczByd7VcK5CtuA6dMrJ0LXQQgu5tddeu7ZMQJYL2HHDhjFZb7313Jxzzun/1l133YbT6/Q2WWedddx8881n8rhh47LX6bISWj9kq1+/fu6DDz5IcqDzw+G//e1v3W7YhcYRYAgsw+Evvvii8QQ7PAWGwYMGDXImjx3e0C2oHrLFcPgPf/hDt9xtYqQbJOVdsImRMCw15XVYShbbEIgjoMmWkWAcq1LPjATD4NQENSwli20IxBHQZCuIBM3UeBzYvLOiJMie2KuvvtpPnFx00UWO56LhjTfe8PevvPJKJ3/XX399qpKX5xh6R9P78MMPo8lV9rcmqNFC59UvD89oWsnfYH322We7yy67zH388cfJ23bepghoslWIBM3UeH0tX4QEmbq/4YYb3E477eQ22GADP4O133771XQXbPLfYYcdXFdXV+yPCa0nnniiW8FoK4iUZ0ivf//+7uc//7kjnaoHTVCl7Hn1Y3+xhqekk3bkuW233dYNGTLEjR8/PnU5Rdpzdq36CGiyVYgEzdR4fY2cR4L0aF544QV33HHHuU8//dRnct1117lFFlnEnXDCCd5gwIsvvuh22203t91227ntt9/eH5lp5iX985//HCsYEzGvvPKKmzx5smNXBeGaa65xK6ywgnvsscdicat4ogkq5c2rH3iCVxae9BDTAtfPP/98Pymz11572cRMGkhtfk2TrUIkKPU3U+OCRLFjHgn++9//9iQYXcAJGY4ePdqNGjXKD8feeecdT4bRHCHKyy+/vNt1SIKeUjQwvIM8f/Ob30QvV/K3JqgUWKvfM88848ATEszCUz4M0cqTJj3AAQMGuDPPPDN6y353EAKabBUmQYZt9E5GjhwZE7IOwqn0quSRYFqGDFvHjRvnDj744NSdE6R5yCGHuOeeey7t8dg12uzaa691N910k9cTxm5W8EQT1LTiFqmf4AlmaTtR6E2PGDHCDRs2zOtkzzrrLHfrrbemZWfX2hgBTbaMBJvYsPWQID23bbbZxt1///3dttrRa6FHh47vo48+yiw5w0LiHXDAAW7VVVf1PR10ZVUPmqBGyx5SPw1P0nz44YfdEkss4Xve9L6HDh3qBg4c6C644AL32WefRbO1322MgCZbRoJNbNh6SBBLx8ccc0zqRAZExst5ySWXdBsKR6vBsJBZ5OHDh/thXp8+fRyOiSDRKgdNUKPlDqmfhidpMgvMJNPhhx9ew3Tvvff2uwseeuihaLb2u40R0GTLSLCJDRtKgvjPYB83vZe0EDIUlufRH6Lvir7kcq9qR01Qs8qq1S8PTz4KkOCaa67p3n777VoW9957r1tsscXchRdeWPkPR63Q9kNFQJOtIBI0U+Mqzt1uhpAgexpnzpzpXn/9dZ8OM5bvv/9+TZdXdCicLAQTLWPHjnXTpk2r9XSScapyrglqVhmz6peHJ+mB6W233eZWXnnl2MQRM+xrrLGGn1mveu85Cxe7HkdAk63CJEiSDNNYqiHLOeLZ2FkSgaIkiH6PiRCWwjB8Yyh76qmnOoZjsqyj6FA4WQbS3mOPPdx9991XI9RknKqca4KaVca0+ml4ok+MBmaSN998c3fYYYfVLtMDZMM9PUkLnYGAJluFSFAWqC6//PJeaTx9+nSHkxoLOgJ5JCgK/l133TW2EJqF0csss0xsFp5GZKIj+WKivGc2kz8MDzCUPuecc2pmuzDhBaF2ymJpespZ9ZPechaebJ7//PPPPVa33HKLx4tnbrzxRrfSSiu5SZMm+R4zC6YxVyUfIL2V7W47INAwCfK1ZBHpVltt5ZfIjBkzptASjXYAp5llzCNBUfCDqyyGZk3fjjvu6FiqEX0JeYHZCZJcB0ge9GL4owc0a9Yst8UWW9TMdvFis5SkHYImqFL+Bx54ILN+eXhCeHwMwOrQQw+NzbCD9yabbOI22mgjN2PGDMnOjh2CgCZbhXqCHYJDj1cjjwR7vEAVz1AT1IoX3YpXcQQ02TISbGLjGQmGgasJalhKFtsQiCOgyZaRYByrUs+MBMPg1AQ1LCWLbQjEEdBkyxwtxbEq/UwsS9tSi3xo0V2ao6V8nCxGOALIVqajJUwtsT5N7NTZ8f9t9jWKxXnnnecdLeFsid+NptfpzzM5we4Nk8fyZLDTZaZo/ZAtuI4eYTJ0saUqaavOzuO2+wwPw8NkoP1lYMEFF+xmsBhC7Fp44YW9FQ2WZthfuRiw9GWuueby3ub4bfjq+GLNBfUBR8NKx8rwCcMHmWKUkeltzqxlJDvI5ZyjBxSdYHKnQjk5dFYqLGTG25zJY2e1axVqg2yZt7kWtITNDoeBrs3ghaVksQ2BOAKabBVaIsPOhajjHl5uC/kI1EOC7PrAxl2yN8R2MRwAsWUsrUufLA3+R5hgwFpyO5jWp/yaoCbrJ+dl4UV65khMUO28oyZbuSTI1HK9jms6D8qwGoWSIJaPsYCMgyTM6ktgyxyWYFZffXW31lpruYMOOih1lov4DLsxvMDWOQyEYjGZvbBcY1tZlYMmqGnlLgMv0mUrIh+LRRdd1K299to1J1dpedq19kRAky2VBHmhNEdA0b2t7QlNc0sdSoJ33nmnJ7rNNtusZqACXcbxxx/vvcdhvQcjCSjFcSbEvWTAXDykh+UYCdgSxC0CvaYqB01Q08pdBl6ky9743Xff3Q0ePNgboqUcFjoLAU22VBKk5wAJhjiu6SzoGqtNCAmCMRZf6OVhXp/eHwHbdhtvvLG74447aoXh5efaSy+9VLsmP+hBYgQAE/wSWEuFd7r33ntPLlXyqAlqssBl4RVN1xyJRdHorN+abKkkmAZDnuOatGd667WiJIj+Dxt2DFnp9W299da1niBGPzHwGdXrPf744/4a95IBh+ETJ070bjaxIYg648gjj3RYYK56z10T1Gg9y8RL0gUncyQmaHTeUZOtYBLMc1zTefDVX6MiJMgymnvuucebeeflFuvdYq8R15roAnElKYHeIQ7CMQ2fth3vrbfe8sPf5ZZbzpuFuvTSS1M9rUl6VTlqgiplbAZepG0kKAh35lGTrWASzHNc05kQ1lerPBLkhUaHx3BVwumnn+71f+IeUkgwOvTFkxxGQLNIkLSwPQhRskqeXmY7BE1QKX8z8TISbAcJqb+MmmwFkWCe45r6i9iZT+aRICbzzz33XHfggQf62UkIEH3giiuu6KZOnepee+01P4ylJ4hzcQkPPvigW3bZZb1P4WRPkMksHAVhOJQeI9aoUfhDmGkTKZJmFY6aoFK+ZuAl9TYSFCQ686jJVmESLOK4pjPhq79WeSTIi3fKKad4HxcbbrihW3/99f2Oib59+/qJD3wPP/XUU97fBUNmCZjSZ/nLk08+KZdqRxp7yy239AYb5CLm9bGanOXFTuK1+qgJKmVrBl5SZ9IWVUR0IlDu27G9EdBkqxAJsrQizxFQe0PUnNLnkWBarsxQsgRGnFnJRAdLYiSwnGPfffd13EuGl19+2fckUfJLoEfI3knuVTlogppV7kbxiqZrjsSiaHTWb022VBJkaIX+KctxjX0xdUEJJUF6I0cddZRfLD179uxa4vT4Ro0a5Y499ljvP5ghMztCCEym0DPEcRC/Wfh7xBFH+Nlj4uMwiOUyRx99dKoFjVomFfihCWpa8crAi3TNkVgaup11TZMtlQTzHNdUfclFq5sxlATBE89n9OJ4NhogOYa0w4cPd08//XTtFvHEcZA8Q7tBgKwlZJjdDo7XqZAmqLUKR36UhRcfc3MkFgG2A39qsqWSYAdi0aNVCiXBHi1cBTPTBLWCxbUitRECmmwZCTaxIY0Ew8DVBDUsJYttCMQR0GTLSDCOValnRoJhcGqCGpaSxTYE4ghosmUkGMeq1DMjwTA4NUENS8liGwJxBDTZMm9zcaxKPxPL0slFzaVn1AEJMttr3uY6oCErWAVkK9PbHDbUpkyZ4pdSsJzC/srDgF0feJqbZ555vIUYw1bHFuMRLBQ3edRxMjkKxwfZguvSLCl14XwEQ5IspbC/csLESF4AACAASURBVDFYb731vJOlOeec06277rqGb46MrbPOOm6++eYzeczByd7T8PcU2erXr1+qVXY/HMZQp4XyEWAILMNhFp5b0BHAaASOlkwedZzsbjgCyJY5WgrHreEnbGIkDEJNeR2WksU2BOIIaLJls8NxrEo9MxIMg1MT1LCULLYhEEdAk61CJMhQLupt7sMPP4znYGepCISQIMYQ2BqHbUH+Zs2a1c0SdIg3tFDvdKkV6OGLmqAmi5KHV0j92X4XlW/azUJnIaDJVi4JyubyHXbYwW/s79+/v9+Qj5l9CzoCISR47bXXuq6urtrfhAkTavb/Qr2hhXin02vQs3c1QU2WRMMrpP4snTBvikl0O+9cky2VBFHsY4Zp8uTJ7pNPPvHIXHPNNd5/RdTnRedBVk6NipIgEwHTpk1ze+65pxszZozbeeeda1ZiKAkb/It6Qwv1TldOTctJRRPUaA4aXpBaUe98jHDMm2IU2c79rclWLgnSC4kGDHNi7w4TWxZ0BIqSIBa7d9ttN3fzzTerCRbxhobNwBDvdGqGPXxTE9RoUTS8Qupv3hSjqHb2b022VBJMwsJXlmHITTfd5J18J+/beRyBIiSIDUDs/clQeJ999kn9wIB9EW9ot99+e5B3uniJW3umCaqULA+vUO98kq4czZuiINFZR022CpGgGFfFX8Wqq67qdSj4e7CgI1CEBFHw4xKT3vUWW2zh5phjDrfKKqs4PMZFQxESRH0hjplCvNNF82nlb01QpVxZeIkRWhyBhXrnk7Q5mjfFKBqd81uTrUIkKMZVMeg5YMAA16dPH3f33XenunvsHNgar0kREkzmggVoJp/Ychc1WhtKgqHe6ZLlaMW5JqhZ5RG8Tj75ZO+D5KqrrvIkWG/9zZtiFtLtfV2TrUIkGK0+TrwhwnaxVhwte0//rocEKSMm9pkcic7AFyXB66+/3q222mqFvdP1NCZafpqgas+B19ixY72ZfPSq9dbfvClqKLf3PU22gkkQB0AIHLOZNiTWBaNeEqSXHZ2RJxdIsIg3NLzTsRe8qHc6vQY9e1cTVK0k4AURoi989tln66q/eVPUEG7/e5psBZMgnuf22GMPd99999nkSI5s1EOCzMZPnz7dPfroow41RDQU8YaGM/eJEye6ot7poum3+rcmqFlli+KFTpTz0PqbN8UsdDvnuiZbuSTIyvtzzjmnZmLrxBNPdPixjQ7VOgeqcmtShARpnIsuuqiGLzPAYBzFlxebOMsvv7wbOHCgJ8l3333XF5bej3ibE8MDad7pHn/88XIr14TUNEGV7IrglVb/pHc+MGNUg1N786Yo6HbuUZOtXBJ84IEH/KylmO+ZNGmSH5p1Llzl1awICT7//PNu9OjRNTNbafhq3tCYLRVvc/RoJGR5p5P7VTxqgirlLYIXcbPqT5uAF8uSGAKz+H+rrbZy2223nZ+hZ5Z+xx13dGeddVZsYkryt2N7IqDJVi4JtmeVq1HqIiRYjZJWoxSaoFajhFaKdkVAky0jwSa2qpFgGLiaoIalZLENgTgCmmwZCcaxKvXMSDAMTk1Qw1Ky2IZAHAFNtszRUhyr0s/EsjQzlxZ0BFgGZI6WdIzsbn0IIFvmaKkFTqTM0VKYQxxztBSGlzlcKo6XOVpqkeMac7QU5hDHHC2F4SUrNuyYj5s5WqqvB93wUwyBZThsjpby4TRHS/kYWYz6EDBHS/Xh1vBTNjESBqGmvA5LyWIbAnEENNmy2eE4VqWeGQmGwakJalhKFtsQiCOgyVYwCbIr4eGHH/ab1ePZ2FkSgaIkiA07HP2IkyWOWIN57733kkm6oviHOBrqlkmLLmiCGi1SEbzYJjdz5kx35plnujxXEOZoKYpuZ/7WZCuIBBlXH3LIId7h0jvvvNOZaJVYqyIkyL5gnFiJZWk5LrTQQjE/IxSrKP4hjoZKrG7DSWmCKolreGFBBz0sH2kM1A4dOtQNGzbMbbvttu6hhx7qZpCCNFk6YY6WBN3OPWqyFUSCd955pzdYudlmmzkStaAjkEeCTJaIvwzZu8px3XXXdePHj3dYhImGIvh3sqMlDS98tGAQgT9ID0tHErB9OXLkSN+LlmscSc8cLUUR6dzfpZAgm/ixbnLQQQe5bbbZxtHbsKAjkEeCmMp68803u9llxHAtZvKj9hqL4h/iaEgvfc/f1QSV0mh4YREaUsO6zkYbbeQNJEgNUC/wUUmqF0gPEgRbCZAoBi1GjRpV87Ao9+zYvghoslWoJ4i5pgsvvNAPKVh0uPXWW8cEp32haW7J80gwLXeeQeXw3HPP1W6H4N/pjpZqoPzvRxIvsSe4wgoreJuXDHfx4cKHJequIJmOnGPCbNy4cb4NUD9Y6AwEGiJBdCxYKb7sssv8ZEgR68adAVvjtQglQbDGlSlmnsQsVgj+xO10R0vRVknDi/s4qWL4u9xyy7kZM2a4Sy+91OtTo89m/TZHS1nItPf1ukkQIUMvxXBCwumnn+4V+WLAU67bsTsCoSTI8PeCCy5wl1xySW0oHIJ/lATrdTTUvRY9d0UT1LRSpOEl8TBCO2TIELfgggv6UYxczzuao6U8hNrzviZb6nAYITv33HPdgQce6JcaQIDoA1dccUWvH3zttdfaE5EeKnUoCSaHdqH4Q4K9ydFSEi+aFb3gvffe6w2nvvLKKw43sYMHD/YjGSaNtGCOljR02vte3SSIPuWUU05xm2++ubd8vP7667tBgwa5vn37euUzfkZ48SykIxBCguCYHArXg39vcbSUhhetgLBvueWW7rzzzqs1Cu4gNtlkE+9TuHYx8cMcLSUA6bDTukkwDYfTTjvNmyFnFs2CjkAICWpDu2gugv8nn3wSvVz73VscLWXhxew4IxV8tUigRzhixAjHvbSA/vXggw92LLNhOIz6B+JkbWGRyZS0NO1atRAojQTpmeDacIMNNnCzZ8+uVi0rWJoQEqSRGLoxJMsKUfxlsXpvc7Qk2GThxezwEUcc4dZYYw137LHHegdWTDThpJ32ELzM0ZIg2TuOpZEgX8Ubb7zRf2URKAs6AiEkyLpLlPm8xFkhDX96fr3J0ZJgo+HF+j8IcOONN/ZqHBZLS49O8DJHS4Jk7ziWRoK9A67yahlCguXl2r4paYLavrWyklcBAU221NnhKhS+nctgJBjWepqghqVksQ2BOAKabBkJxrEq9cxIMAxOTVDDUrLYhkAcAU22zNFSHKvSz8SytC0lyoeWiR9ztJSPk8UIRwDZMkdL5mjJz5RW2TmPOVoq7jioyu1YxbKZoyVztORnSavukMccLeU7DKp6G1a1fOZoKbz3XMoTDIFlOMx2Lgs6AlhtYUeS7UvXcbK74QggWwyH00wA2sRIOJ6Fn7CJkcJQ+Yia8josJYttCMQR0GTLSDCOValnRoJhcGqCGpaSxTYE4ghoslWYBD/++GN3yy231JwBzZo1q7YKP56dnQkCISSY5xgo1HFSaHwpcyuPmqAmy5WHV6i85qWXzN/O2wsBTbYKk+C1114bcwY0YcIEl2eaqL1gKr+0RUgQXSEb9dMcA4keET3G2LFjvX+XtdZay7s4oFGzQmj8rHR6+romqFIWDS+2y0koKq9F05N07dieCGiyVYgEUVRPmzbN7bnnnm7MmDFu55137uYJrT2haW6pi5Age1mzHAPRm8FaCtP7eKTDcg9tsf3227vjjjsu9SPUyY6WaC0NL7HGHSKvRdJrrpRY6j2BQMMkKB7Rbr755p4ob8fkUYQEsQaT5Rjoww8/dBiuxRDAHXfcUcMFr3Nci1qPlpud7GiJOmp4vf/++x6GEHnV0ks6ZhKM7dh+CDREgpgewuKG+MPdZ599vPHP9oOh50tchATp7U2cONElHQPhCxfLJ5AfZqGiDsQff/xxf+22227rVqlOd7SUhReOlBgOh8qrlp5YnukGsl1oOwQaIkGEBG9dDMHQW80xxxxulVVW8c5s2g6JHi5wERKkSGmOgUTfGuI4iXWJIfF7GI7c7DRBjT6s4VWPvKalZ57mooi3/29NtgrpBKMQYJyyf//+burUqTY7HAUm5XdREuTRLMdAQmrRoS9m+FdaaSXvNyO6JzlKgkXipxS5pZc0QU0WLAuvZLyi8lo0vWT6dt4eCGiyFUyCVBnr0kyO4KPVQjYCRUiQ2ck0x0CYeac3iBHb1Vdf3T3zzDO1jB588EG37LLLOmZAkyTY6Y6WsvDCJaz0nmtA/e+HJq/1pJdM386rj0DpJHj33Xe7yZMnuyw/F9WHpGdKWIQEaZw0x0Cbbrqpe/vtt92zzz7r2PeI72cJmIYfOnSoe/LJJ+VS7djpjpay8NIcKWnyWk96NbDtR9sgUCoJYv59+vTp7tFHH/WK6LZBoQUFLUKCmmMgZobpbe+///5+SYxUYffdd3f77ruvQ/+VDCz5YKKFJTQStPgSpwpHTVClfBpeaY6U8uQ1ND0phx3bCwFNtnKHwzyMvkTM4+DF68QTT7ShcAEZKEKCmmMgCI1Aj2/UqFHebwb+MvD9zA4HgjgOYjePGB5Ii8+MctWDJqhSdg0v8M6T1yhe/M5LT/K1Y3sjoMlWLgk+//zzbvTo0TVTTJMmTXIYKLSQj0AREiSVNMdALJKOBkiOId/w4cPd008/XbtFHr3N0VIaXrKcJU9eo3jxm5CWXhL/GuD2oy0RaIgE27LGFSl0URKsSHFbXgxNUFteOCtAWyOgyVZuT7Cta97iwhsJhjWAJqhhKVlsQyCOgCZbRoJxrEo9MxIMg1MT1LCULLYhEEdAky1ztBTHqvQzsSwdXc9XeiYdkqA5WuqQhqxgNczRUgucLDGbzq6aeeed180zzzx+Rl1m2O2Y7lDIHC2l42Ly0jgu5mjJHC3VZver6giHcpmjJXO01Cz5NEdLLeqeMwSW4TDbsyzoCJijJR0fu1s/AuZoqX7sGnrSJkbC4NOU12EpWWxDII6AJls2OxzHqtQzI8EwODVBDUvJYhsCcQQ02QomwTfeeMOdffbZ3oxT2t7VeNa9+yyUBB955BF3//33uzRbdnmOk9gWd+aZZ/q2IW5eCI2fl14Z9zVBTUtfw4v4efeJU0SeUWVcffXVfusoW0ix+G2hvRDQZCuIBLF2jD+MIUOGuPHjx7vf//737YVED5e2KAmyfxUCW3TRRd3aa6/dzUF0nuOkhx9+2G299dZu2LBhbs011/Rmzl555ZXM2kIOIfEzEyr5hiao0azy8Mq7L2kVkWfSgvjw8bLBBht4W5pYWjczcoJiexw12SpEguzLPP/8892gQYPcXnvtVduo3x7Vb10pi5IgHxMsvQwePNjvDabBJGiOk1j7RFw+TBi2IGAUgHOcYvECJwN5hcRPPt/Mc01Qo/lqeBEv735ReWZii49J1GzcNddc410hRN0dRMtmv6uJgCZbuSSIIPDFHDBggO+tVLOK1SxVURKU0p922mlu5MiRsR42pp7SHC1hTIF7d911l++hvPnmm5KMwxETRlcxvpoM+CVZf/31XdH4yeebea4Jalq+aXhF46XdD5Fn4iY/JAyfcTWBdW8L7YOAJlu5JIg5pxEjRvihFos2zzrrLIdRTwv5CISQIL06enNJEsxynMSwF8syF198sVt55ZXdiy++WCsQ5LfUUkv5YRwvsgR+E3/FFVcsFF+e66mjJqjJMmThJfGy7jciz6SJNe+bbrrJ2ZInQbo9jpps5ZIg+qYllljC27PDph0WjQcOHOguuOACP/RqDwhaU8pGSRDSEh8jUZJjiAYJcg8ipJc+c+bMWiVnzZrlll566VTz+/QE8RFTJH4twR76oQlqsghZJCfxsu7XI88QHj2/Aw44wK266qp+ZGSmtgTp9jhqspVLgvhuWGihhRzGPKXh9957b9evXz/30EMPtQcCLSplmSSY5jjpyiuvdPjGHTt2rFt88cXdlClTPLlhdZoPF9a/kwHfvLvsskvh+Mnnm3muCWoy3yySk3hZ9+uRZ+wNgjW2HPng9OnTx2GyP9rLlnztWE0ENNlSSZBGRmjodeDvQgKOgRZbbDF34YUXmiAIKCnHMkgwz3ES2TIRwIQVukP0fbykOHRnVjktzJ49Oyh+WhrNuKYJajK/LJKTeGn3y5Bn/BtDhNFOgeRpx+oioMlWLgkyfELnFFUEMxzDITgzZfY1zG74UBI86aSTuukEQx0nsUYQ/9C0TRG9VWj87No2fkcT1GTqkFwaXhIv7T6y2qg8f/rpp77nPW3atNrISPK0Y3UR0GRLJUGqRC9j88039ybcpYr0ANmQ/Nvf/lYu2TEFgRAS5PFjjjnGbbfddo4XTUKI46Q77rjDr2fD3H40DUkreQyNn3y+7HNNUNPySsMrGi/tfqPy/NFHH7k99tjD3XfffYU+MtHy2O/WIaDJVi4J8vXE9y3OvvEvwheQdWbMFItfh9ZVrdo5FyVBWZC7/PLL+0knvPm9++67tcqlOU4SR0s0Lh8l9IGbbbaZmzBhghMHTSQgjoWY0ec3OkQtfi3TFvzQBDVanDy8tPtZ8sxideQ5iheOq1gSc84559QcjeFk7NRTT7XF0tEGaYPfmmzlkqDUj6UxrE1D1zRjxgy5bEcFgaIkKDq9rbbayg+Hx4wZ45577rlYylmOlnAstOOOOzr8FBMnGSgDPUP+IEcmWLT4yed78lwT1Gg58vDKu09aWfIcxYteHzPtqBfExJM5Gou2RPv81mSrMAm2T3WrU9KiJFidEre2JJqgtrZklnu7I6DJlpFgE1vXSDAMXE1Qw1Ky2IZAHAFNtowE41iVemYkGAanJqhhKVlsQyCOgCZb5mgpjlXpZ2JZ2pYS5UPLsha2+6WZEst/2mIYAtkIIFtLLrlkbF++xO7CfBMzi+bMpXFnLkkMzdFSGKbmaCkMr6S82Xk2fuZoyRwt1WY3ZZazikdztGSOlpoll+ZoSfq9PXxkCCzD4SK7N3q4eJXLzhwtVa5JOqZA5mipRU1pEyNhwGvK67CULLYhEEdAky2bHY5jVeqZkWAYnJqghqVksQ2BOAKabBUiQbYO4WgGc0Lyh3UTtmBZyEagKAmyXSvqyIfnoiHvfjQuuxywmcf2r7xQxBFRXhpl3tcENZpPHh6oHqJ4ao6RiuAVkl60nPa7OghospVLguzDxMlMV1dX7A8bg7J/tTpVrVZJipAgU/e4L9hpp528mXxm6/fbb7+aGSxsOGr3ozVG73HIIYf4dN55553ordhv2lRz7BSL3IMnmqBKMfLwwgFSUcdIRfCSfcjmaElaoD2PmmypJMgXEEsxu+22m7dugm8FrJysu+663ttcdKN+e0LT3FLnkSD4vvDCC+64446rWX3BXt0iiyziTe1DgPgRybpPjyga8C2y+uqre0MKNHpWYG9tlmOnrGd64romqOSfhxd4vPbaa4UdI+XhxcSWOVrqiZZvfh6abKkkiEVdHPKIRWkpKi8qpt2T1+W+Hf+LQB4Jgi8kCClJwATW6NGj/R9+nSHBtPu4Ovjkk0/kMR8HCycHHXSQ22abbWo9yVqElB9pjohSovXYJU1QKYSGF3jQC0yaEMtyjASmeXhBgvQEoyErvWgc+109BDTZUkkwrSq82Ay5klZO0uL29mt5JJiGDy/yuHHj3MEHH5y6c0Lu0wYM5wjo/zCPhbsDFoXiUzhKnGn5MKxMc+yUFrenrmmCmlUGwSMNL+qY5hipHrzIPyu9rLLZ9eogoMlWEAnyZcTCNM6nUShb0BGohwTpadCTu//++1Otdifv0yb33HOPd4PAy61ZW46WtlNIMIkHdWTYnOUYqR68tPSimNrv6iJQGgky/MXL3CWXXGJD4QLtXQ8JXnHFFd7CND2ctJC8j16WGXsJp59+up/IwiCoFjqFBJN4UOcsx0iQGSqGULyy0oNQLbQHAqWRoA2Fwxo8lASZhGIfN72btJC8z0t47rnnugMPPNDP9kKA9CLxK4y+i0mCrNAJJJjEI62u4hjpyCOP9PpCPuL14CVpS3rmaEkQaY9jKSTIC2dD4bAGDyHBDz74wLvLxPERgZlO3GNKbyPtPstgGP7iA4Y9l3iaGzRokOvbt6+3AI4fDHk+WXJIsOjQOflss841QU3mmYYHeNHbiwYmSnBJygeC37iHqAcvSVPSM0dLgkh7HDXZKqwTtKFweGMXJUH0qyj2WYrE8I7hGn4sWMxMIJ20+yyKTgZmfFnKFJ05TsaR8zRHRHKvFUdNUKPlycIrzQ82ccUxUtoHIQQvyhBNL0m40TLa72ohoMlWYRIkkQMOOMA8zAW0bR4JisJ91113jS1EZ2H64MGD/Qzvs88+69LuL7PMMt1mgOndHXXUUbHF0kyW4GQJ/yP8JsgC4CzHTgFVLDWqJqhkpOEleLCkq6hjJA0vMDNHS6U2b0sT02SrMAniyJuV+Ml1Uy2tWcUzzyNBUbjjYIlF6PTg+MMREi/y559/7v0Hp93HUVBysTTneAZk6Qt5E5g4wcnSoYceWrvG8hmctWuOnVoBrSaolEfD6+yzz/Z40Bss6hhJwwvM6PWZo6VWSEL5eWqyVZgEyy9W56eYR4Kdj0BYDTVBDUvJYhsCcQQ02TISjGNV6pmRYBicmqCGpWSxDYE4AppsGQnGsSr1zEgwDE5NUMNSstiGQBwBTbbM0VIcq9LPxLJ02sxk6Zm1eYJMVJijpTZvxIoWH9nKdLTUv39/vz5N7ATa8f9tJjaKxXnnnefmnXde/8fvRtPr9OeZ7MFE28yZMw2riO3OTm/3nqgfsgXX0SNMhq4+ffp0W56RtB1o53FbioaH4WEy0H4ysOCCC9ZWSESJsGvhhRd2I0aMqC3PkGUadvzvcpVGcGAJylxzzeXmnHNOvxylkbR6w7PIIeoDk8fGZa83yEtIHZEpRhnsNEoGrxOURbTJm3beGALoAUUnaLsL8rFkXSTb/kwe87GyGGEIIFvoBFnvnAw2O5xEpMRzmx0OA1ObwQtLyWIbAnEENNkqTIL4E0FhjW+Kxx57LJ6DnaUiUJQE2bkQdQzEc1mBXQyaIyUs0LB74rLLLvNmo9LSoVcazU9zRJT2fLOuaYIazbMIXpjMYqugKN3Z+SG9cTCi/nKPo+Y4rKp4RTGx3zoCmmzlkiACIFuRhg4d6oYNG+a23XZbf41tTBayEShCgkzdl+VIiXRomyFDhngfMGnWpWXfcBUdB2mCKigXxQuL0tHJiwkTJvhtdewHDnEcVmW8BBM75iOgyVYuCbL3lBcLSxwSsKU2cuRIsy4tgGQc80iQD4zmaIkeTzRkOQYi3vnnn+/1aewJzjKoio6yyo6DNEEFh6J4UX9MXe25555uzJgxbuedd3ZPPvmkNyuGW4iijsOqjldUNuy3joAmW7kkiM26jTbayJvUl2wYPowfP978DgsgGcc8EqQnneVoKc+RkvTyeFHpAQ4YMMCrKjKK4i8TN2kAg6Ehs2zYimx10ASVsuXhJU6WxEPizTffHKsSz4c4Dqs6XrHK2YmKgCZbuSSIbmXixIluhRVWcBjpZDiClV4s7CZ7KmopeuHNPBJMg0QcBxV1pERPnel/1BRnnHGGY1EoZqCKBNoyzRFRkWebEUcT1Kz8onhRH5xP4QNHhsL77LOPSvC0UVHHYVXDKwsTu94dAU22ckmQ5N566y0//F1uueXcjBkz3KWXXlrzdNY9O7siCNRDgvTMoo6W6I2kOVKSqX4Mry6xxBKOniN/6G0HDhzofcFkLTVhWJnliEjK3oqjJqhZ5RG8HnjgAR+Fni4faXq3mNSaY4453CqrrOJlOJkG2Baxll5VvJL1sfNsBDTZKkSCJI0tQRTurLrGvaOFfATqIcGk46AsR0ribhP/zywCjfq82HvvvV2/fv385FVaKRkWotIYPny4H0aza+juu+/ONMWflkYzrmmCmpVfEq9kvKOPPtpvl5o6dWq3kUtRa+lVxStZVzvPRkCTrVwS5Ct47733esOcKNWxLo3VY5ZgsADRQjYCoSSYdBxETyXLkRL+QV599VXfK19rrbXc22+/XSsI7bXYYov5jxVpaKFKjoM0QU2rQxKvtDhcw9o2kyNJD34hQ2FJu0p4SZnsmI+AJlu5JMjDW265pcMAgAT8X2yyySaZXtEkXm8/hpBgmuMgzZHSxhtv7HW06P++973vxfRefKzWWGMNb5U6jwSr5DhIE9SkLKXhleZoiefo5U6ePDnmd6XoUDiZb5XwSpbNzrMR0GQrlwRffvll78IRk+0SeMlQxnPPQjYCRUmQBdChjpRkJvTdd9/13tMwBy8BdcU666xTyB9MlRwHaYIqdeOYhVeaoyV0hNOnT3ePPvqon12WdIoOhSW+HKuEl5TJjvkIaLKVS4II0RFHHOF7Fscee6yfgWT2DV0LL7mFbATySFAU7vU4Upo9e3YtY/yKrLTSSm7SpEl+fRzrOtnZw+w9kyP0FtvBcZAmqFQ2Dy8mi+gNor9mppw/Pt74YE4OhckrzXFYO+FVEwD7kYuAJlu5JEjqKIYhQIZg+LeNKuFzc+/FEfJIUBTujThSEnhZGoOKgjWdzOBLoAz0EtvBcZAmqNRHw4utgtx/6aWX3OjRo72cIqt8GFjakgxZjsPaCa9knew8GwFNtgqRYHbSdkdDII8EtWd74z1NUHsjHlbn8hDQZMtIsDycu6VkJNgNEvWCJqjqg3bTEMhBQJMtI8Ec8Bq5bSQYhp4mqGEpWWxDII6AJltGgnGsSj0zEgyDUxPUsJQstiEQR0CTLfM2F8eq9DOxLJ23Xq/0jNswQSYwzNtcGzZcGxQZ2cr0NsdLKktfZFmBHf+7vKJRHFiaIRv5WarRaHqd/vyUKVPc3HPPbfL4v+U9nd7ePVk/ZAuuY11tMnSNGzeutvSFJQX2Vx4GLFcRZzD8Nmx1bFmCxRpHWYpleOl4GT7F8UGm4DrZZBAlwq7oif02BAwBQ6C3IWAk2Nta3OprCBgCajIsnQAAADlJREFUMQSMBGNw2IkhYAj0NgSMBHtbi1t9DQFDIIaAkWAMDjsxBAyB3oaAkWBva3GrryFgCMQQ+D/Ux0GNsQPuKgAAAABJRU5ErkJggg==" alt="" data-iml="1495164.834999945"><br> <p>What am I doing wrong in the diagonalization?</p><p>Thanks a lot for your help.</p><p>Regards,</p><p>Torstein Fjermestad</p><p><br></p><p><br></p> <p><br></p><br></div><div class="gmail_quote"><div dir="auto" class="gmail_attr">onsdag 10. juni 2020 kl. 15:56:53 UTC+2 skrev Torstein Fjermestad:<br/></div><blockquote class="gmail_quote" style="margin: 0 0 0 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><div dir="ltr">Dear all, <div><br></div><div>I would like to use TAMkin (<a href="http://molmod.github.io/tamkin/" target="_blank" rel="nofollow" data-saferedirecturl="https://www.google.com/url?hl=no&q=http://molmod.github.io/tamkin/&source=gmail&ust=1608498665747000&usg=AFQjCNF-zHDW2zv8LU7aADsqeN_tyJOmuw">http://molmod.github.io/tamkin/</a>) to compute thermodynamic properties from the vibrational frequencies computed with CP2K.</div><div>As I have understood it, the way TAMkin works is to read and diagonalize the Hessian matrix printed in the CP2K output.</div><div>When I compare the frequencies obtained with TAMkin with those obtained with cp2k, the values are not the same (Compare column 2 and 3 in the table below). </div><div><br></div><div>The first task is then to understand how CP2K computes the frequencies. To investigate this, I copied the Hessian matrix from the CP2K output (printed below the string "Hessian in cartesian coordinates") and pasted it into Excel and saved it as a csv file. This file I read into python and diagonalized the Hessian with numpy. The frequencies I get are different from those printed in the cp2k output file (compare column 3 and 4 in the table below). </div><div>I have assumed that the units of the Hessian matrix in the cp2k output file are Ry*bohr-2*kamu-1 (kamu = 1000*amu). This assumption is based on trial and error; it was what got me closest to the cp2k frequencies. </div><div><br></div><div><br></div><div>Table 1: The eight lowest frequencies. There are 48 in total.</div><div><table border="1" cellspacing="0" cellpadding="0" style="border-collapse:collapse;border-width:initial;border-style:none"> <tbody><tr style="height:30.0pt"> <td width="85" nowrap valign="top" style="width:63.8pt;border-width:1pt;border-style:solid;border-color:windowtext;padding:0cm 5.4pt;height:30pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">frequencies</span><u></u><u></u></p> </td> <td width="75" nowrap valign="top" style="width:56.45pt;border-top:1pt solid windowtext;border-right:1pt solid windowtext;border-bottom:1pt solid windowtext;border-left-width:initial;border-left-style:none;padding:0cm 5.4pt;height:30pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">TAMkin / <br> cm-1<u></u><u></u></span></p> </td> <td width="76" valign="top" style="width:2cm;border-top:1pt solid windowtext;border-right:1pt solid windowtext;border-bottom:1pt solid windowtext;border-left-width:initial;border-left-style:none;padding:0cm 5.4pt;height:30pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">cp2k / <br> cm-1<u></u><u></u></span></p> </td> <td width="151" valign="top" style="width:4cm;border-top:1pt solid windowtext;border-right:1pt solid windowtext;border-bottom:1pt solid windowtext;border-left-width:initial;border-left-style:none;padding:0cm 5.4pt;height:30pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">Diag., cp2k Hessian / <br> cm-1<u></u><u></u></span></p> </td> </tr> <tr style="height:15.0pt"> <td width="85" nowrap valign="top" style="width:63.8pt;border-right:1pt solid windowtext;border-bottom:1pt solid windowtext;border-left:1pt solid windowtext;border-top-width:initial;border-top-style:none;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">1<u></u><u></u></span></p> </td> <td width="75" nowrap valign="top" style="width:56.45pt;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">31.1<u></u><u></u></span></p> </td> <td width="76" nowrap valign="top" style="width:2cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">27.3<u></u><u></u></span></p> </td> <td width="151" nowrap valign="top" style="width:4cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">32.6<u></u><u></u></span></p> </td> </tr> <tr style="height:15.0pt"> <td width="85" nowrap valign="top" style="width:63.8pt;border-right:1pt solid windowtext;border-bottom:1pt solid windowtext;border-left:1pt solid windowtext;border-top-width:initial;border-top-style:none;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">2<u></u><u></u></span></p> </td> <td width="75" nowrap valign="top" style="width:56.45pt;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">38.3<u></u><u></u></span></p> </td> <td width="76" nowrap valign="top" style="width:2cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">40.8<u></u><u></u></span></p> </td> <td width="151" nowrap valign="top" style="width:4cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">40.1<u></u><u></u></span></p> </td> </tr> <tr style="height:15.0pt"> <td width="85" nowrap valign="top" style="width:63.8pt;border-right:1pt solid windowtext;border-bottom:1pt solid windowtext;border-left:1pt solid windowtext;border-top-width:initial;border-top-style:none;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">3<u></u><u></u></span></p> </td> <td width="75" nowrap valign="top" style="width:56.45pt;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">47.7<u></u><u></u></span></p> </td> <td width="76" nowrap valign="top" style="width:2cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">51.0<u></u><u></u></span></p> </td> <td width="151" nowrap valign="top" style="width:4cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">50.2<u></u><u></u></span></p> </td> </tr> <tr style="height:15.0pt"> <td width="85" nowrap valign="top" style="width:63.8pt;border-right:1pt solid windowtext;border-bottom:1pt solid windowtext;border-left:1pt solid windowtext;border-top-width:initial;border-top-style:none;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">4<u></u><u></u></span></p> </td> <td width="75" nowrap valign="top" style="width:56.45pt;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">74.7<u></u><u></u></span></p> </td> <td width="76" nowrap valign="top" style="width:2cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">74.0<u></u><u></u></span></p> </td> <td width="151" nowrap valign="top" style="width:4cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">78.2<u></u><u></u></span></p> </td> </tr> <tr style="height:15.0pt"> <td width="85" nowrap valign="top" style="width:63.8pt;border-right:1pt solid windowtext;border-bottom:1pt solid windowtext;border-left:1pt solid windowtext;border-top-width:initial;border-top-style:none;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">5<u></u><u></u></span></p> </td> <td width="75" nowrap valign="top" style="width:56.45pt;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">82.0<u></u><u></u></span></p> </td> <td width="76" nowrap valign="top" style="width:2cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">80.0<u></u><u></u></span></p> </td> <td width="151" nowrap valign="top" style="width:4cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">85.8<u></u><u></u></span></p> </td> </tr> <tr style="height:15.0pt"> <td width="85" nowrap valign="top" style="width:63.8pt;border-right:1pt solid windowtext;border-bottom:1pt solid windowtext;border-left:1pt solid windowtext;border-top-width:initial;border-top-style:none;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">6<u></u><u></u></span></p> </td> <td width="75" nowrap valign="top" style="width:56.45pt;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">111.0<u></u><u></u></span></p> </td> <td width="76" nowrap valign="top" style="width:2cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">109.9<u></u><u></u></span></p> </td> <td width="151" nowrap valign="top" style="width:4cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">116.3<u></u><u></u></span></p> </td> </tr> <tr style="height:15.0pt"> <td width="85" nowrap valign="top" style="width:63.8pt;border-right:1pt solid windowtext;border-bottom:1pt solid windowtext;border-left:1pt solid windowtext;border-top-width:initial;border-top-style:none;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">7<u></u><u></u></span></p> </td> <td width="75" nowrap valign="top" style="width:56.45pt;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">223.7<u></u><u></u></span></p> </td> <td width="76" nowrap valign="top" style="width:2cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">224.1<u></u><u></u></span></p> </td> <td width="151" nowrap valign="top" style="width:4cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">234.3<u></u><u></u></span></p> </td> </tr> <tr style="height:15.0pt"> <td width="85" nowrap valign="top" style="width:63.8pt;border-right:1pt solid windowtext;border-bottom:1pt solid windowtext;border-left:1pt solid windowtext;border-top-width:initial;border-top-style:none;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">8<u></u><u></u></span></p> </td> <td width="75" nowrap valign="top" style="width:56.45pt;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">245.6<u></u><u></u></span></p> </td> <td width="76" nowrap valign="top" style="width:2cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">246.3<u></u><u></u></span></p> </td> <td width="151" nowrap valign="top" style="width:4cm;border-top-width:initial;border-style:none solid solid none;border-left-width:initial;border-bottom-width:1pt;border-bottom-color:windowtext;border-right-width:1pt;border-right-color:windowtext;padding:0cm 5.4pt;height:15pt"> <p class="MsoNormal" style="margin-bottom:0cm;margin-bottom:.0001pt;line-height:normal"><span lang="IT">257.3<u></u><u></u></span></p> </td> </tr> </tbody></table></div><div><br></div><div>What is wrong with my calculation?</div><div>What steps should be done to go from the Hessian to the frequencies?</div><div><br></div><div>The Hessian in csv format, the python script and the cp2k output file are attached. </div><div><br></div><div>Thanks a lot for your help.</div><div><br></div><div>Regards,</div><div>Torstein Fjermestad</div><div> </div><div><br></div><div><br></div></div></blockquote></div>