diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 68a9257..3de134c 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -30,7 +30,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest, macOS-latest] #windows-latest (vina is not available for windows in conda-forge),m think about use bioconda autodock-vina - python-version: [3.8, 3.9, '3.10', 3.11] + python-version: [3.8, 3.9, '3.10', 3.11, 3.12, 3.13] steps: - uses: actions/checkout@v4 - name: Set up Python ${{ matrix.python-version }} diff --git a/LICENSE b/LICENSE index 261eeb9..ee9ed41 100644 --- a/LICENSE +++ b/LICENSE @@ -1,201 +1,504 @@ - Apache License - Version 2.0, January 2004 - http://www.apache.org/licenses/ - - TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION - - 1. Definitions. - - "License" shall mean the terms and conditions for use, reproduction, - and distribution as defined by Sections 1 through 9 of this document. - - "Licensor" shall mean the copyright owner or entity authorized by - the copyright owner that is granting the License. - - "Legal Entity" shall mean the union of the acting entity and all - other entities that control, are controlled by, or are under common - control with that entity. For the purposes of this definition, - "control" means (i) the power, direct or indirect, to cause the - direction or management of such entity, whether by contract or - otherwise, or (ii) ownership of fifty percent (50%) or more of the - outstanding shares, or (iii) beneficial ownership of such entity. - - "You" (or "Your") shall mean an individual or Legal Entity - exercising permissions granted by this License. - - "Source" form shall mean the preferred form for making modifications, - including but not limited to software source code, documentation - source, and configuration files. - - "Object" form shall mean any form resulting from mechanical - transformation or translation of a Source form, including but - not limited to compiled object code, generated documentation, - and conversions to other media types. - - "Work" shall mean the work of authorship, whether in Source or - Object form, made available under the License, as indicated by a - copyright notice that is included in or attached to the work - (an example is provided in the Appendix below). - - "Derivative Works" shall mean any work, whether in Source or Object - form, that is based on (or derived from) the Work and for which the - editorial revisions, annotations, elaborations, or other modifications - represent, as a whole, an original work of authorship. For the purposes - of this License, Derivative Works shall not include works that remain - separable from, or merely link (or bind by name) to the interfaces of, - the Work and Derivative Works thereof. - - "Contribution" shall mean any work of authorship, including - the original version of the Work and any modifications or additions - to that Work or Derivative Works thereof, that is intentionally - submitted to Licensor for inclusion in the Work by the copyright owner - or by an individual or Legal Entity authorized to submit on behalf of - the copyright owner. For the purposes of this definition, "submitted" - means any form of electronic, verbal, or written communication sent - to the Licensor or its representatives, including but not limited to - communication on electronic mailing lists, source code control systems, - and issue tracking systems that are managed by, or on behalf of, the - Licensor for the purpose of discussing and improving the Work, but - excluding communication that is conspicuously marked or otherwise - designated in writing by the copyright owner as "Not a Contribution." - - "Contributor" shall mean Licensor and any individual or Legal Entity - on behalf of whom a Contribution has been received by Licensor and - subsequently incorporated within the Work. - - 2. Grant of Copyright License. Subject to the terms and conditions of - this License, each Contributor hereby grants to You a perpetual, - worldwide, non-exclusive, no-charge, royalty-free, irrevocable - copyright license to reproduce, prepare Derivative Works of, - publicly display, publicly perform, sublicense, and distribute the - Work and such Derivative Works in Source or Object form. - - 3. Grant of Patent License. Subject to the terms and conditions of - this License, each Contributor hereby grants to You a perpetual, - worldwide, non-exclusive, no-charge, royalty-free, irrevocable - (except as stated in this section) patent license to make, have made, - use, offer to sell, sell, import, and otherwise transfer the Work, - where such license applies only to those patent claims licensable - by such Contributor that are necessarily infringed by their - Contribution(s) alone or by combination of their Contribution(s) - with the Work to which such Contribution(s) was submitted. If You - institute patent litigation against any entity (including a - cross-claim or counterclaim in a lawsuit) alleging that the Work - or a Contribution incorporated within the Work constitutes direct - or contributory patent infringement, then any patent licenses - granted to You under this License for that Work shall terminate - as of the date such litigation is filed. - - 4. Redistribution. You may reproduce and distribute copies of the - Work or Derivative Works thereof in any medium, with or without - modifications, and in Source or Object form, provided that You - meet the following conditions: - - (a) You must give any other recipients of the Work or - Derivative Works a copy of this License; and - - (b) You must cause any modified files to carry prominent notices - stating that You changed the files; and - - (c) You must retain, in the Source form of any Derivative Works - that You distribute, all copyright, patent, trademark, and - attribution notices from the Source form of the Work, - excluding those notices that do not pertain to any part of - the Derivative Works; and - - (d) If the Work includes a "NOTICE" text file as part of its - distribution, then any Derivative Works that You distribute must - include a readable copy of the attribution notices contained - within such NOTICE file, excluding those notices that do not - pertain to any part of the Derivative Works, in at least one - of the following places: within a NOTICE text file distributed - as part of the Derivative Works; within the Source form or - documentation, if provided along with the Derivative Works; or, - within a display generated by the Derivative Works, if and - wherever such third-party notices normally appear. The contents - of the NOTICE file are for informational purposes only and - do not modify the License. You may add Your own attribution - notices within Derivative Works that You distribute, alongside - or as an addendum to the NOTICE text from the Work, provided - that such additional attribution notices cannot be construed - as modifying the License. - - You may add Your own copyright statement to Your modifications and - may provide additional or different license terms and conditions - for use, reproduction, or distribution of Your modifications, or - for any such Derivative Works as a whole, provided Your use, - reproduction, and distribution of the Work otherwise complies with - the conditions stated in this License. - - 5. Submission of Contributions. Unless You explicitly state otherwise, - any Contribution intentionally submitted for inclusion in the Work - by You to the Licensor shall be under the terms and conditions of - this License, without any additional terms or conditions. - Notwithstanding the above, nothing herein shall supersede or modify - the terms of any separate license agreement you may have executed - with Licensor regarding such Contributions. - - 6. Trademarks. This License does not grant permission to use the trade - names, trademarks, service marks, or product names of the Licensor, - except as required for reasonable and customary use in describing the - origin of the Work and reproducing the content of the NOTICE file. - - 7. Disclaimer of Warranty. Unless required by applicable law or - agreed to in writing, Licensor provides the Work (and each - Contributor provides its Contributions) on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or - implied, including, without limitation, any warranties or conditions - of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A - PARTICULAR PURPOSE. You are solely responsible for determining the - appropriateness of using or redistributing the Work and assume any - risks associated with Your exercise of permissions under this License. - - 8. Limitation of Liability. In no event and under no legal theory, - whether in tort (including negligence), contract, or otherwise, - unless required by applicable law (such as deliberate and grossly - negligent acts) or agreed to in writing, shall any Contributor be - liable to You for damages, including any direct, indirect, special, - incidental, or consequential damages of any character arising as a - result of this License or out of the use or inability to use the - Work (including but not limited to damages for loss of goodwill, - work stoppage, computer failure or malfunction, or any and all - other commercial damages or losses), even if such Contributor - has been advised of the possibility of such damages. - - 9. Accepting Warranty or Additional Liability. While redistributing - the Work or Derivative Works thereof, You may choose to offer, - and charge a fee for, acceptance of support, warranty, indemnity, - or other liability obligations and/or rights consistent with this - License. However, in accepting such obligations, You may act only - on Your own behalf and on Your sole responsibility, not on behalf - of any other Contributor, and only if You agree to indemnify, - defend, and hold each Contributor harmless for any liability - incurred by, or claims asserted against, such Contributor by reason - of your accepting any such warranty or additional liability. - - END OF TERMS AND CONDITIONS - - APPENDIX: How to apply the Apache License to your work. - - To apply the Apache License to your work, attach the following - boilerplate notice, with the fields enclosed by brackets "[]" - replaced with your own identifying information. (Don't include - the brackets!) The text should be enclosed in the appropriate - comment syntax for the file format. We also recommend that a - file or class name and description of purpose be included on the - same "printed page" as the copyright notice for easier - identification within third-party archives. - - Copyright [yyyy] [name of copyright owner] - - Licensed under the Apache License, Version 2.0 (the "License"); - you may not use this file except in compliance with the License. - You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - - Unless required by applicable law or agreed to in writing, software - distributed under the License is distributed on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - See the License for the specific language governing permissions and - limitations under the License. + GNU LESSER GENERAL PUBLIC LICENSE + Version 2.1, February 1999 + + Copyright (C) 1991, 1999 Free Software Foundation, Inc. + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + +[This is the first released version of the Lesser GPL. It also counts + as the successor of the GNU Library Public License, version 2, hence + the version number 2.1.] + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +Licenses are intended to guarantee your freedom to share and change +free software--to make sure the software is free for all its users. + + This license, the Lesser General Public License, applies to some +specially designated software packages--typically libraries--of the +Free Software Foundation and other authors who decide to use it. You +can use it too, but we suggest you first think carefully about whether +this license or the ordinary General Public License is the better +strategy to use in any particular case, based on the explanations below. + + When we speak of free software, we are referring to freedom of use, +not price. Our General Public Licenses are designed to make sure that +you have the freedom to distribute copies of free software (and charge +for this service if you wish); that you receive source code or can get +it if you want it; that you can change the software and use pieces of +it in new free programs; and that you are informed that you can do +these things. + + To protect your rights, we need to make restrictions that forbid +distributors to deny you these rights or to ask you to surrender these +rights. These restrictions translate to certain responsibilities for +you if you distribute copies of the library or if you modify it. + + For example, if you distribute copies of the library, whether gratis +or for a fee, you must give the recipients all the rights that we gave +you. You must make sure that they, too, receive or can get the source +code. If you link other code with the library, you must provide +complete object files to the recipients, so that they can relink them +with the library after making changes to the library and recompiling +it. And you must show them these terms so they know their rights. + + We protect your rights with a two-step method: (1) we copyright the +library, and (2) we offer you this license, which gives you legal +permission to copy, distribute and/or modify the library. + + To protect each distributor, we want to make it very clear that +there is no warranty for the free library. Also, if the library is +modified by someone else and passed on, the recipients should know +that what they have is not the original version, so that the original +author's reputation will not be affected by problems that might be +introduced by others. + + Finally, software patents pose a constant threat to the existence of +any free program. We wish to make sure that a company cannot +effectively restrict the users of a free program by obtaining a +restrictive license from a patent holder. Therefore, we insist that +any patent license obtained for a version of the library must be +consistent with the full freedom of use specified in this license. + + Most GNU software, including some libraries, is covered by the +ordinary GNU General Public License. This license, the GNU Lesser +General Public License, applies to certain designated libraries, and +is quite different from the ordinary General Public License. We use +this license for certain libraries in order to permit linking those +libraries into non-free programs. + + When a program is linked with a library, whether statically or using +a shared library, the combination of the two is legally speaking a +combined work, a derivative of the original library. The ordinary +General Public License therefore permits such linking only if the +entire combination fits its criteria of freedom. The Lesser General +Public License permits more lax criteria for linking other code with +the library. + + We call this license the "Lesser" General Public License because it +does Less to protect the user's freedom than the ordinary General +Public License. It also provides other free software developers Less +of an advantage over competing non-free programs. These disadvantages +are the reason we use the ordinary General Public License for many +libraries. However, the Lesser license provides advantages in certain +special circumstances. + + For example, on rare occasions, there may be a special need to +encourage the widest possible use of a certain library, so that it becomes +a de-facto standard. To achieve this, non-free programs must be +allowed to use the library. A more frequent case is that a free +library does the same job as widely used non-free libraries. In this +case, there is little to gain by limiting the free library to free +software only, so we use the Lesser General Public License. + + In other cases, permission to use a particular library in non-free +programs enables a greater number of people to use a large body of +free software. For example, permission to use the GNU C Library in +non-free programs enables many more people to use the whole GNU +operating system, as well as its variant, the GNU/Linux operating +system. + + Although the Lesser General Public License is Less protective of the +users' freedom, it does ensure that the user of a program that is +linked with the Library has the freedom and the wherewithal to run +that program using a modified version of the Library. + + The precise terms and conditions for copying, distribution and +modification follow. Pay close attention to the difference between a +"work based on the library" and a "work that uses the library". The +former contains code derived from the library, whereas the latter must +be combined with the library in order to run. + + GNU LESSER GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License Agreement applies to any software library or other +program which contains a notice placed by the copyright holder or +other authorized party saying it may be distributed under the terms of +this Lesser General Public License (also called "this License"). +Each licensee is addressed as "you". + + A "library" means a collection of software functions and/or data +prepared so as to be conveniently linked with application programs +(which use some of those functions and data) to form executables. + + The "Library", below, refers to any such software library or work +which has been distributed under these terms. A "work based on the +Library" means either the Library or any derivative work under +copyright law: that is to say, a work containing the Library or a +portion of it, either verbatim or with modifications and/or translated +straightforwardly into another language. (Hereinafter, translation is +included without limitation in the term "modification".) + + "Source code" for a work means the preferred form of the work for +making modifications to it. For a library, complete source code means +all the source code for all modules it contains, plus any associated +interface definition files, plus the scripts used to control compilation +and installation of the library. + + Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running a program using the Library is not restricted, and output from +such a program is covered only if its contents constitute a work based +on the Library (independent of the use of the Library in a tool for +writing it). Whether that is true depends on what the Library does +and what the program that uses the Library does. + + 1. You may copy and distribute verbatim copies of the Library's +complete source code as you receive it, in any medium, provided that +you conspicuously and appropriately publish on each copy an +appropriate copyright notice and disclaimer of warranty; keep intact +all the notices that refer to this License and to the absence of any +warranty; and distribute a copy of this License along with the +Library. + + You may charge a fee for the physical act of transferring a copy, +and you may at your option offer warranty protection in exchange for a +fee. + + 2. You may modify your copy or copies of the Library or any portion +of it, thus forming a work based on the Library, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) The modified work must itself be a software library. + + b) You must cause the files modified to carry prominent notices + stating that you changed the files and the date of any change. + + c) You must cause the whole of the work to be licensed at no + charge to all third parties under the terms of this License. + + d) If a facility in the modified Library refers to a function or a + table of data to be supplied by an application program that uses + the facility, other than as an argument passed when the facility + is invoked, then you must make a good faith effort to ensure that, + in the event an application does not supply such function or + table, the facility still operates, and performs whatever part of + its purpose remains meaningful. + + (For example, a function in a library to compute square roots has + a purpose that is entirely well-defined independent of the + application. Therefore, Subsection 2d requires that any + application-supplied function or table used by this function must + be optional: if the application does not supply it, the square + root function must still compute square roots.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Library, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Library, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote +it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Library. + +In addition, mere aggregation of another work not based on the Library +with the Library (or with a work based on the Library) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may opt to apply the terms of the ordinary GNU General Public +License instead of this License to a given copy of the Library. To do +this, you must alter all the notices that refer to this License, so +that they refer to the ordinary GNU General Public License, version 2, +instead of to this License. (If a newer version than version 2 of the +ordinary GNU General Public License has appeared, then you can specify +that version instead if you wish.) Do not make any other change in +these notices. + + Once this change is made in a given copy, it is irreversible for +that copy, so the ordinary GNU General Public License applies to all +subsequent copies and derivative works made from that copy. + + This option is useful when you wish to copy part of the code of +the Library into a program that is not a library. + + 4. You may copy and distribute the Library (or a portion or +derivative of it, under Section 2) in object code or executable form +under the terms of Sections 1 and 2 above provided that you accompany +it with the complete corresponding machine-readable source code, which +must be distributed under the terms of Sections 1 and 2 above on a +medium customarily used for software interchange. + + If distribution of object code is made by offering access to copy +from a designated place, then offering equivalent access to copy the +source code from the same place satisfies the requirement to +distribute the source code, even though third parties are not +compelled to copy the source along with the object code. + + 5. A program that contains no derivative of any portion of the +Library, but is designed to work with the Library by being compiled or +linked with it, is called a "work that uses the Library". Such a +work, in isolation, is not a derivative work of the Library, and +therefore falls outside the scope of this License. + + However, linking a "work that uses the Library" with the Library +creates an executable that is a derivative of the Library (because it +contains portions of the Library), rather than a "work that uses the +library". The executable is therefore covered by this License. +Section 6 states terms for distribution of such executables. + + When a "work that uses the Library" uses material from a header file +that is part of the Library, the object code for the work may be a +derivative work of the Library even though the source code is not. +Whether this is true is especially significant if the work can be +linked without the Library, or if the work is itself a library. The +threshold for this to be true is not precisely defined by law. + + If such an object file uses only numerical parameters, data +structure layouts and accessors, and small macros and small inline +functions (ten lines or less in length), then the use of the object +file is unrestricted, regardless of whether it is legally a derivative +work. (Executables containing this object code plus portions of the +Library will still fall under Section 6.) + + Otherwise, if the work is a derivative of the Library, you may +distribute the object code for the work under the terms of Section 6. +Any executables containing that work also fall under Section 6, +whether or not they are linked directly with the Library itself. + + 6. As an exception to the Sections above, you may also combine or +link a "work that uses the Library" with the Library to produce a +work containing portions of the Library, and distribute that work +under terms of your choice, provided that the terms permit +modification of the work for the customer's own use and reverse +engineering for debugging such modifications. + + You must give prominent notice with each copy of the work that the +Library is used in it and that the Library and its use are covered by +this License. You must supply a copy of this License. If the work +during execution displays copyright notices, you must include the +copyright notice for the Library among them, as well as a reference +directing the user to the copy of this License. Also, you must do one +of these things: + + a) Accompany the work with the complete corresponding + machine-readable source code for the Library including whatever + changes were used in the work (which must be distributed under + Sections 1 and 2 above); and, if the work is an executable linked + with the Library, with the complete machine-readable "work that + uses the Library", as object code and/or source code, so that the + user can modify the Library and then relink to produce a modified + executable containing the modified Library. (It is understood + that the user who changes the contents of definitions files in the + Library will not necessarily be able to recompile the application + to use the modified definitions.) + + b) Use a suitable shared library mechanism for linking with the + Library. A suitable mechanism is one that (1) uses at run time a + copy of the library already present on the user's computer system, + rather than copying library functions into the executable, and (2) + will operate properly with a modified version of the library, if + the user installs one, as long as the modified version is + interface-compatible with the version that the work was made with. + + c) Accompany the work with a written offer, valid for at + least three years, to give the same user the materials + specified in Subsection 6a, above, for a charge no more + than the cost of performing this distribution. + + d) If distribution of the work is made by offering access to copy + from a designated place, offer equivalent access to copy the above + specified materials from the same place. + + e) Verify that the user has already received a copy of these + materials or that you have already sent this user a copy. + + For an executable, the required form of the "work that uses the +Library" must include any data and utility programs needed for +reproducing the executable from it. However, as a special exception, +the materials to be distributed need not include anything that is +normally distributed (in either source or binary form) with the major +components (compiler, kernel, and so on) of the operating system on +which the executable runs, unless that component itself accompanies +the executable. + + It may happen that this requirement contradicts the license +restrictions of other proprietary libraries that do not normally +accompany the operating system. Such a contradiction means you cannot +use both them and the Library together in an executable that you +distribute. + + 7. You may place library facilities that are a work based on the +Library side-by-side in a single library together with other library +facilities not covered by this License, and distribute such a combined +library, provided that the separate distribution of the work based on +the Library and of the other library facilities is otherwise +permitted, and provided that you do these two things: + + a) Accompany the combined library with a copy of the same work + based on the Library, uncombined with any other library + facilities. This must be distributed under the terms of the + Sections above. + + b) Give prominent notice with the combined library of the fact + that part of it is a work based on the Library, and explaining + where to find the accompanying uncombined form of the same work. + + 8. You may not copy, modify, sublicense, link with, or distribute +the Library except as expressly provided under this License. Any +attempt otherwise to copy, modify, sublicense, link with, or +distribute the Library is void, and will automatically terminate your +rights under this License. However, parties who have received copies, +or rights, from you under this License will not have their licenses +terminated so long as such parties remain in full compliance. + + 9. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Library or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Library (or any work based on the +Library), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Library or works based on it. + + 10. Each time you redistribute the Library (or any work based on the +Library), the recipient automatically receives a license from the +original licensor to copy, distribute, link with or modify the Library +subject to these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties with +this License. + + 11. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Library at all. For example, if a patent +license would not permit royalty-free redistribution of the Library by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Library. + +If any portion of this section is held invalid or unenforceable under any +particular circumstance, the balance of the section is intended to apply, +and the section as a whole is intended to apply in other circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 12. If the distribution and/or use of the Library is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Library under this License may add +an explicit geographical distribution limitation excluding those countries, +so that distribution is permitted only in or among countries not thus +excluded. In such case, this License incorporates the limitation as if +written in the body of this License. + + 13. The Free Software Foundation may publish revised and/or new +versions of the Lesser General Public License from time to time. +Such new versions will be similar in spirit to the present version, +but may differ in detail to address new problems or concerns. + +Each version is given a distinguishing version number. If the Library +specifies a version number of this License which applies to it and +"any later version", you have the option of following the terms and +conditions either of that version or of any later version published by +the Free Software Foundation. If the Library does not specify a +license version number, you may choose any version ever published by +the Free Software Foundation. + + 14. If you wish to incorporate parts of the Library into other free +programs whose distribution conditions are incompatible with these, +write to the author to ask for permission. For software which is +copyrighted by the Free Software Foundation, write to the Free +Software Foundation; we sometimes make exceptions for this. Our +decision will be guided by the two goals of preserving the free status +of all derivatives of our free software and of promoting the sharing +and reuse of software generally. + + NO WARRANTY + + 15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO +WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW. +EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR +OTHER PARTIES PROVIDE THE LIBRARY "AS IS" WITHOUT WARRANTY OF ANY +KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE +LIBRARY IS WITH YOU. SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME +THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN +WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY +AND/OR REDISTRIBUTE THE LIBRARY AS PERMITTED ABOVE, BE LIABLE TO YOU +FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR +CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE +LIBRARY (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING +RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A +FAILURE OF THE LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF +SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH +DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Libraries + + If you develop a new library, and you want it to be of the greatest +possible use to the public, we recommend making it free software that +everyone can redistribute and change. You can do so by permitting +redistribution under these terms (or, alternatively, under the terms of the +ordinary General Public License). + + To apply these terms, attach the following notices to the library. It is +safest to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least the +"copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 + USA + +Also add information on how to contact you by electronic and paper mail. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the library, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the + library `Frob' (a library for tweaking knobs) written by James Random + Hacker. + + , 1 April 1990 + Ty Coon, President of Vice + +That's all there is to it! \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 57e542d..bfd3e7e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -44,7 +44,7 @@ keywords = [ "genetic algorithm", ] -requires-python = ">= 3.8 , < 3.12" +requires-python = ">= 3.8" dependencies = [ "crem", "tqdm", @@ -52,10 +52,8 @@ dependencies = [ "pandas", "pyyaml", "dill", - "meeko>=0.4.0,<0.6.0", - "six", + "scipy", "rdkit>=2022.3.5", - "scipy", # a meeko dependency ] [project.license] @@ -68,9 +66,9 @@ dev = ["requests", "pytest"] default-version = "1+unknown" [tool.versioningit.format] -distance = "{base_version}" -dirty = "{base_version}" -distance-dirty = "{base_version}" +distance = "{base_version}.post{distance}" +dirty = "{base_version}.post{distance}.dev0" +distance-dirty = "{base_version}.post{distance}.dev0" [tool.versioningit.vcs] method = "git" @@ -95,6 +93,7 @@ moldrug = [ "data/*/*.smi", "data/*/*.pdb", "data/*/*.pdbqt", + "meeko_light/data/params/*.json" ] diff --git a/src/moldrug/fitness.py b/src/moldrug/fitness.py index ff5fd55..46036e9 100644 --- a/src/moldrug/fitness.py +++ b/src/moldrug/fitness.py @@ -6,8 +6,8 @@ import numpy as np # from warnings import import warn -from meeko import (MoleculePreparation, PDBQTMolecule, PDBQTWriterLegacy, - RDKitMolCreate) +from moldrug.meeko_light import (MoleculePreparation, + PDBQTMolecule, PDBQTWriterLegacy, RDKitMolCreate) from rdkit import Chem from rdkit.Chem import QED, Descriptors diff --git a/src/moldrug/meeko_light/__init__.py b/src/moldrug/meeko_light/__init__.py new file mode 100644 index 0000000..4f6435b --- /dev/null +++ b/src/moldrug/meeko_light/__init__.py @@ -0,0 +1,28 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# This module was taken from Meeko 0.6.1 +# https://github.com/forlilab/meeko +# Only unnecessary functionalities were removed. + +__meeko_original_version__ = "0.6.1" + +from .preparation import MoleculePreparation +from .molecule_pdbqt import PDBQTMolecule +from .rdkit_mol_create import RDKitMolCreate +from .writer import PDBQTWriterLegacy + +import logging +from rdkit import rdBase +rdkit_logger = logging.getLogger("rdkit") +rdkit_logger.handlers[0].setLevel("WARNING") +rdkit_logger.handlers[0].setFormatter( + logging.Formatter('[RDKit] %(levelname)s:%(message)s'), +) +rdBase.LogToPythonLogger() + +__all__ = ['MoleculePreparation', + 'PDBQTMolecule', + 'RDKitMolCreate', + 'PDBQTWriterLegacy', + ] diff --git a/src/moldrug/meeko_light/atomtyper.py b/src/moldrug/meeko_light/atomtyper.py new file mode 100644 index 0000000..18d4179 --- /dev/null +++ b/src/moldrug/meeko_light/atomtyper.py @@ -0,0 +1,330 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import warnings +from typing import List + +import numpy as np + +from .utils import pdbutils + + +class AtomTyper: + + @classmethod + def type_everything( + cls, + molsetup, + atom_params, + charge_model, + offatom_params=None, + dihedral_params=None, + ): + + cls._type_atoms(molsetup, atom_params) + + # offatoms must be typed after charges, because offsites pull charge + if offatom_params is not None: + cached_offatoms = cls._cache_offatoms(molsetup, offatom_params) + coords = {atom.index: atom.coord for atom in molsetup.atoms if not atom.is_dummy} + cls._set_offatoms(molsetup, cached_offatoms, coords) + + return + + @staticmethod + def _type_atoms(molsetup, atom_params): + # ensure every "atompar" is defined in a single "smartsgroup" + ensure = {} + # go over all "smartsgroup"s + for smartsgroup in atom_params: + if smartsgroup == "comment": + continue + for line in atom_params[ + smartsgroup + ]: # line is a dict, e.g. {"smarts": "[#1][#7,#8,#9,#15,#16]","atype": "HD"} + smarts = str(line["smarts"]) + # get indices of the atoms in the smarts to which the parameters will be assigned + idxs = [ + 0 + ] # by default, the first atom in the smarts gets parameterized + if "IDX" in line: + idxs = [i for i in line["IDX"]] + # match SMARTS + hits = molsetup.find_pattern(smarts) + for atompar in line: + if atompar in ["smarts", "comment", "IDX"]: + continue + if atompar not in molsetup.atom_params: + molsetup.atom_params[atompar] = [None] * len(molsetup.atoms) + value = line[atompar] + # keep track of every "smartsgroup" that modified "atompar" + ensure.setdefault(atompar, []) + ensure[atompar].append(smartsgroup) + # Each "hit" is a tuple of atom indices that matched the smarts + # The length of each "hit" is the number of atoms in the smarts + for hit in hits: + # Multiple atoms may be targeted by a single smarts: + # For example: both oxygens in NO2 are parameterized by a single smarts pattern. + # "idxs" are 1-indices of atoms in the smarts to which parameters are to be assigned. + for idx in idxs: + if atompar == "atype": + molsetup.set_atom_type( + hit[idx], value + ) # overrides previous calls + molsetup.atom_params[atompar][hit[idx]] = value + + # guarantee that each atompar is exclusive of a single group + for atompar in ensure: + if len(set(ensure[atompar])) > 1: + msg = "%s is modified in multiple smartsgroups: %s" % ( + atompar, + set(ensure[atompar]), + ) + warnings.warn(msg) + return + + @staticmethod + def _cache_offatoms(molsetup, offatom_params): + """precalculate off-site atoms""" + cached_offatoms = {} + n_offatoms = 0 + atoms_with_offchrg = set() + # each parent atom can only be matched once in each smartsgroup + for smartsgroup in offatom_params: + if smartsgroup == "comment": + continue + tmp = {} + for line in offatom_params[smartsgroup]: + # SMARTS + smarts = str(line["smarts"]) + hits = molsetup.find_pattern(smarts) + # atom indexes in smarts string + smarts_idxs = [0] + if "IDX" in line: + smarts_idxs = [i for i in line["IDX"]] + for smarts_idx in smarts_idxs: + for hit in hits: + parent_idx = hit[smarts_idx] + tmp.setdefault( + parent_idx, [] + ) # TODO tmp[parent_idx] = [], yeah? + for offatom in line["OFFATOMS"]: + # set defaults + tmp[parent_idx].append( + { + "offatom": { + "distance": 1.0, + "x90": False, + "phi": 0.0, + "theta": 0.0, + "z": [], + "x": [], + }, + "atom_params": {}, + } + ) + for key in offatom: + if key in ["distance", "x90"]: + tmp[parent_idx][-1]["offatom"][key] = offatom[key] + # replace SMARTS indexes by the atomic index + elif key in ["z", "x"]: + for i in offatom[key]: + idx = hit[i] + tmp[parent_idx][-1]["offatom"][key].append(idx) + # convert degrees to radians + elif key in ["theta", "phi"]: + tmp[parent_idx][-1]["offatom"][key] = np.radians( + offatom[key] + ) + # ignore comments + elif key in ["comment"]: + pass + elif key == "atype": + tmp[parent_idx][-1]["atom_params"][key] = offatom[ + key + ] + elif key == "pull_charge_fraction": + if parent_idx in atoms_with_offchrg: + raise RuntimeError( + "atom %d has charge pulled more than once" + % parent_idx + ) + atoms_with_offchrg.add(parent_idx) + tmp[parent_idx][-1]["atom_params"][key] = offatom[ + key + ] + else: + pass + for parent_idx in tmp: + for offatom_dict in tmp[parent_idx]: + atom_params = offatom_dict["atom_params"] + offatom = offatom_dict["offatom"] + atomgeom = AtomicGeometry( + parent_idx, + neigh=offatom["z"], + xneigh=offatom["x"], + x90=offatom["x90"], + ) + if "pull_charge_fraction" in atom_params: + pull_charge_fraction = atom_params["pull_charge_fraction"] + else: + pull_charge_fraction = 0.0 + args = ( + atom_params["atype"], + offatom["distance"], + offatom["theta"], + offatom["phi"], + pull_charge_fraction, + ) + # number of coordinates (before adding new offatom) + cached_offatoms[n_offatoms] = (atomgeom, args) + n_offatoms += 1 + return cached_offatoms + + @staticmethod + def _set_offatoms(molsetup, cached_offatoms, coords): + """add cached offatoms""" + for _, (atomgeom, args) in cached_offatoms.items(): + (atom_type, dist, theta, phi, pull_charge_fraction) = args + offatom_coords = atomgeom.calc_point(dist, theta, phi, coords) + tmp = molsetup.get_pdbinfo(atomgeom.parent + 1) + pdbinfo = pdbutils.PDBAtomInfo( + "G", tmp.resName, tmp.resNum, tmp.icode, tmp.chain + ) + q_parent = (1 - pull_charge_fraction) * molsetup.get_charge(atomgeom.parent) + q_offsite = pull_charge_fraction * molsetup.get_charge(atomgeom.parent) + pseudo_atom = { + "coord": offatom_coords, + "anchor_list": [atomgeom.parent], + "charge": q_offsite, + "pdbinfo": pdbinfo, + "atom_type": atom_type, + "rotatable": False, + } + molsetup.atoms[atomgeom.parent].charge = q_parent + molsetup.add_pseudoatom(**pseudo_atom) + return + + +class AtomicGeometry: + """generate reference frames and add extra sites""" + + def __init__(self, parent, neigh, xneigh: List = None, x90=False, planar_tol=0.1): + """arguments are indices of atoms""" + + self.planar_tol = planar_tol # angstroms, length of neighbor vecs for z-axis + + # real atom hosting extra sites + if not isinstance(parent, int): + raise RuntimeError("parent must be int") + self.parent = parent + + # list of bonded atoms (used to define z-axis) + self.neigh = [] + for i in neigh: + if not isinstance(i, int): + raise RuntimeError("neigh indices must be int") + self.neigh.append(i) + + # list of atoms that + self.xneigh = [] + if isinstance(xneigh, List): + for i in xneigh: + if not isinstance(i, int): + raise RuntimeError("xneigh indices must be int") + self.xneigh.append(i) + + self.calc_x = len(self.xneigh) > 0 + self.x90 = x90 # y axis becomes x axis (useful to rotate in-plane by 90 deg) + + def calc_point(self, distance, theta, phi, coords): + """return coordinates of point specified in spherical coordinates""" + + z = self._calc_z(coords) + + # return pt aligned with z-axis + if phi == 0.0: + return z * distance + np.array(coords[self.parent]) + + # need x-vec if phi != 0 + elif self.calc_x is False: + raise RuntimeError("phi must be zero if X undefined") + + else: + x = self._calc_x(coords) + if self.x90: + x = np.cross(self.z, x) + y = np.cross(z, x) + pt = z * distance + pt = self.rot3D(pt, y, phi) + pt = self.rot3D(pt, z, theta) + pt += np.array(coords[self.parent]) + return pt + + def _calc_z(self, coords): + """maximize distance from neigh""" + z = np.zeros(3) + cumsum = np.zeros(3) + for i in self.neigh: + v = np.array(coords[self.parent]) - np.array(coords[i]) + cumsum += v + z += self.normalized(v) + z = self.normalized(z) + if np.sum(cumsum**2) < self.planar_tol**2: + raise RuntimeError("Refusing to place Z axis on planar atom") + return z + + def _calc_x(self, coords): + x = np.zeros(3) + for i in self.xneigh: + v = np.array(coords[self.parent]) - np.array(coords[i]) + x += self.normalized(v) + x = self.normalized(x) + return x + + @staticmethod + def rot3D(pt, ax, rad): + """ + Rotate point: + pt = (x,y,z) coordinates to be rotated + ax = vector around wich rotation is performed + rad = rotate by "rad" radians + """ + # If axis has len=0, rotate by 0.0 rad on any axis + # Make sure ax has unitary length + len_ax = (ax[0] ** 2 + ax[1] ** 2 + ax[2] ** 2) ** 0.5 + if len_ax == 0.0: + u, v, w = (1, 0, 0) + rad = 0.0 + else: + u, v, w = [i / len_ax for i in ax] + x, y, z = pt + ux, uy, uz = u * x, u * y, u * z + vx, vy, vz = v * x, v * y, v * z + wx, wy, wz = w * x, w * y, w * z + sa = np.sin(rad) + ca = np.cos(rad) + p0 = ( + u * (ux + vy + wz) + + (x * (v * v + w * w) - u * (vy + wz)) * ca + + (-wy + vz) * sa + ) + p1 = ( + v * (ux + vy + wz) + + (y * (u * u + w * w) - v * (ux + wz)) * ca + + (wx - uz) * sa + ) + p2 = ( + w * (ux + vy + wz) + + (z * (u * u + v * v) - w * (ux + vy)) * ca + + (-vx + uy) * sa + ) + return (p0, p1, p2) + + def normalized(self, vec): + l = sum([x**2 for x in vec]) ** 0.5 + if isinstance(vec, list): + return [x / l for x in vec] + else: + # should be np.array + return vec / l diff --git a/src/moldrug/meeko_light/bondtyper.py b/src/moldrug/meeko_light/bondtyper.py new file mode 100644 index 0000000..f73d4a5 --- /dev/null +++ b/src/moldrug/meeko_light/bondtyper.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Meeko bond typer +# + + +class BondTyperLegacy: + + def __call__( + self, setup, flexible_amides, rigidify_bonds_smarts, rigidify_bonds_indices + ): + """Typing atom bonds in the legacy way + + Args: + setup: MoleculeSetup object + + rigidify_bond_smarts (list): patterns to freeze bonds, e.g. conjugated carbons + """ + + amide_bonds = [ + (x[0], x[1]) for x in setup.find_pattern("[NX3]-[CX3]=[O,N,S]") + ] # includes amidines + + # tertiary amides with non-identical substituents will be allowed to rotate + tertiary_amides = [ + x for x in setup.find_pattern("[NX3]([!#1])([!#1])-[CX3]=[O,N,S]") + ] + equivalent_atoms = setup.get_equivalent_atoms() + num_amides_removed = 0 + num_amides_originally = len(amide_bonds) + for x in tertiary_amides: + r1, r2 = x[1], x[2] + if equivalent_atoms[r1] != equivalent_atoms[r2]: + amide_bonds.remove((x[0], x[3])) + num_amides_removed += 1 + assert num_amides_originally == num_amides_removed + len(amide_bonds) + + to_rigidify = set() + n_smarts = len(rigidify_bonds_smarts) + assert n_smarts == len(rigidify_bonds_indices) + for i in range(n_smarts): + a, b = rigidify_bonds_indices[i] + smarts = rigidify_bonds_smarts[i] + indices_list = setup.find_pattern(smarts) + for indices in indices_list: + atom_a = indices[a] + atom_b = indices[b] + to_rigidify.add((atom_a, atom_b)) + to_rigidify.add((atom_b, atom_a)) + + for bond_id, bond in setup.bond_info.items(): + rotatable = True + if bond_id in to_rigidify: + rotatable = False + # check if bond is amide + if ( + bond_id in amide_bonds + or (bond_id[1], bond_id[0]) in amide_bonds + ) and not flexible_amides: + rotatable = False + bond.rotatable = rotatable diff --git a/src/moldrug/meeko_light/data/params/ad4_desolv_param.json b/src/moldrug/meeko_light/data/params/ad4_desolv_param.json new file mode 100644 index 0000000..18cd8b0 --- /dev/null +++ b/src/moldrug/meeko_light/data/params/ad4_desolv_param.json @@ -0,0 +1,11 @@ +{"ad4_desolv_param": [ + {"smarts": "[*]", "ad4_sol_par": -0.00110}, + {"smarts": "[#1]", "ad4_sol_par": 0.00000}, + {"smarts": "[#1][#8,#7,#16,#9]", "ad4_sol_par": 0.00051}, + {"smarts": "[C]", "ad4_sol_par": -0.00143}, + {"smarts": "[c]", "ad4_sol_par": -0.00052}, + {"smarts": "[#7]", "ad4_sol_par": -0.00162}, + {"smarts": "[#8]", "ad4_sol_par": -0.00251}, + {"smarts": "[#16]", "ad4_sol_par": -0.00214} + ] +} diff --git a/src/moldrug/meeko_light/data/params/ad4_desolv_volume.json b/src/moldrug/meeko_light/data/params/ad4_desolv_volume.json new file mode 100644 index 0000000..0cd0beb --- /dev/null +++ b/src/moldrug/meeko_light/data/params/ad4_desolv_volume.json @@ -0,0 +1,19 @@ +{"ad4_desolv_volume": [ + {"smarts": "[#1]", "ad4_sol_vol": 0.0}, + {"smarts": "[#1][#8,#7,#16,#9]", "ad4_sol_vol": 0.0}, + {"smarts": "[#6]", "ad4_sol_vol": 33.5103}, + {"smarts": "[#7]", "ad4_sol_vol": 22.4493}, + {"smarts": "[#8]", "ad4_sol_vol": 17.1573}, + {"smarts": "[#9]", "ad4_sol_vol": 15.448 }, + {"smarts": "[#12]", "ad4_sol_vol": 1.56 }, + {"smarts": "[#15]", "ad4_sol_vol": 38.7924}, + {"smarts": "[#16]", "ad4_sol_vol": 33.5103}, + {"smarts": "[#17]", "ad4_sol_vol": 35.8235}, + {"smarts": "[#20]", "ad4_sol_vol": 2.77 }, + {"smarts": "[#25]", "ad4_sol_vol": 2.14 }, + {"smarts": "[#26]", "ad4_sol_vol": 1.84 }, + {"smarts": "[#30]", "ad4_sol_vol": 1.7 }, + {"smarts": "[#35]", "ad4_sol_vol": 42.5661}, + {"smarts": "[#53]", "ad4_sol_vol": 55.0585} + ] +} diff --git a/src/moldrug/meeko_light/data/params/ad4_hb.json b/src/moldrug/meeko_light/data/params/ad4_hb.json new file mode 100644 index 0000000..b2320e7 --- /dev/null +++ b/src/moldrug/meeko_light/data/params/ad4_hb.json @@ -0,0 +1,7 @@ +{"ad4_hb": [ + { "ad4_hb_hydrogen": true, "smarts": "[#1][#7,#8,#9,#15,#16]"}, + {"ad4_hb_rij": 1.9, "ad4_hb_epsij": 5.0, "smarts": "[#8]"}, + {"ad4_hb_rij": 1.9, "ad4_hb_epsij": 5.0, "smarts": "[#7;!+1;!+2;!$([#7X3v3][a]);!$([#7X3v3][#6X3v4])]"}, + {"ad4_hb_rij": 2.5, "ad4_hb_epsij": 1.0, "smarts": "[SX2]"} + ] +} diff --git a/src/moldrug/meeko_light/data/params/ad4_types.json b/src/moldrug/meeko_light/data/params/ad4_types.json new file mode 100644 index 0000000..1d29d03 --- /dev/null +++ b/src/moldrug/meeko_light/data/params/ad4_types.json @@ -0,0 +1,27 @@ +{"ad4_types": [ + {"smarts": "[#1]", "atype": "H", "comment": "invisible"}, + {"smarts": "[#1][#7,#8,#9,#15,#16]","atype": "HD"}, + {"smarts": "[#5]", "atype": "B"}, + {"smarts": "[C]", "atype": "C"}, + {"smarts": "[c]", "atype": "A"}, + {"smarts": "[#7]", "atype": "NA"}, + {"smarts": "[#8]", "atype": "OA"}, + {"smarts": "[#9]", "atype": "F"}, + {"smarts": "[#12]", "atype": "Mg"}, + {"smarts": "[#14]", "atype": "Si"}, + {"smarts": "[#15]", "atype": "P"}, + {"smarts": "[#16]", "atype": "S"}, + {"smarts": "[#17]", "atype": "Cl"}, + {"smarts": "[#20]", "atype": "Ca"}, + {"smarts": "[#25]", "atype": "Mn"}, + {"smarts": "[#26]", "atype": "Fe"}, + {"smarts": "[#30]", "atype": "Zn"}, + {"smarts": "[#35]", "atype": "Br"}, + {"smarts": "[#53]", "atype": "I"}, + {"smarts": "[#7X3v3][a]", "atype": "N", "comment": "pyrrole, aniline"}, + {"smarts": "[#7X3v3][#6X3v4]", "atype": "N", "comment": "amide"}, + {"smarts": "[#7X3v3][NX2]=[*]", "atype": "N", "comment": "triazene"}, + {"smarts": "[#7+1]", "atype": "N", "comment": "ammonium, pyridinium"}, + {"smarts": "[SX2]", "atype": "SA", "comment": "sulfur acceptor"} + ] +} diff --git a/src/moldrug/meeko_light/data/params/ad4_vdw.json b/src/moldrug/meeko_light/data/params/ad4_vdw.json new file mode 100644 index 0000000..ef321ca --- /dev/null +++ b/src/moldrug/meeko_light/data/params/ad4_vdw.json @@ -0,0 +1,18 @@ +{"ad4_vdw": [ + {"smarts": "[#1]", "rmin_half": 1.0, "epsilon": 0.02}, + {"smarts": "[#6]", "rmin_half": 2.0, "epsilon": 0.15}, + {"smarts": "[#7]", "rmin_half": 1.75, "epsilon": 0.16}, + {"smarts": "[#8]", "rmin_half": 1.6, "epsilon": 0.20}, + {"smarts": "[#9]", "rmin_half": 1.545, "epsilon": 0.08}, + {"smarts": "[#12]", "rmin_half": 0.65, "epsilon": 0.875}, + {"smarts": "[#15]", "rmin_half": 2.1, "epsilon": 0.20}, + {"smarts": "[#16]", "rmin_half": 2.0, "epsilon": 0.20}, + {"smarts": "[#17]", "rmin_half": 2.045, "epsilon": 0.276}, + {"smarts": "[#20]", "rmin_half": 0.99, "epsilon": 0.55}, + {"smarts": "[#25]", "rmin_half": 0.65, "epsilon": 0.875}, + {"smarts": "[#26]", "rmin_half": 0.65, "epsilon": 0.01}, + {"smarts": "[#30]", "rmin_half": 0.74, "epsilon": 0.55}, + {"smarts": "[#35]", "rmin_half": 2.165, "epsilon": 0.389}, + {"smarts": "[#53]", "rmin_half": 2.36, "epsilon": 0.55} + ] +} diff --git a/src/moldrug/meeko_light/data/params/example_offatom_charge.json b/src/moldrug/meeko_light/data/params/example_offatom_charge.json new file mode 100644 index 0000000..a11f07e --- /dev/null +++ b/src/moldrug/meeko_light/data/params/example_offatom_charge.json @@ -0,0 +1,9 @@ +{"example_offsite_charges": [ + {"smarts": "[#7X2;v3;!+](=,:[*])[*]", "IDX": [0], + "OFFATOMS": [ + {"z": [1, 2], "phi": 0, "distance": 0.2, + "atype": "OFFCHRG", "pull_charge_fraction": 1.08} + ] + } + ] +} diff --git a/src/moldrug/meeko_light/data/params/vina_params.json b/src/moldrug/meeko_light/data/params/vina_params.json new file mode 100644 index 0000000..7db36af --- /dev/null +++ b/src/moldrug/meeko_light/data/params/vina_params.json @@ -0,0 +1,31 @@ +{"vina_typing": [ + {"smarts": "[#1]", "vina_ri": null, "vina_donor": false, "vina_acceptor": false, "vina_hydrophobic": false}, + {"smarts": "[#6]", "vina_ri": 1.9, "vina_donor": false, "vina_acceptor": false, "vina_hydrophobic": true}, + {"smarts": "[#7]", "vina_ri": 1.8, "vina_donor": false, "vina_acceptor": false, "vina_hydrophobic": false}, + {"smarts": "[#8]", "vina_ri": 1.7, "vina_donor": false, "vina_acceptor": true, "vina_hydrophobic": false}, + {"smarts": "[#15]", "vina_ri": 2.1, "vina_donor": false, "vina_acceptor": false, "vina_hydrophobic": false}, + {"smarts": "[#16]", "vina_ri": 2.0, "vina_donor": false, "vina_acceptor": false, "vina_hydrophobic": false}, + + {"smarts": "[#9]", "vina_ri": 1.5, "vina_donor": false, "vina_acceptor": false, "vina_hydrophobic": true}, + {"smarts": "[#17]", "vina_ri": 1.8, "vina_donor": false, "vina_acceptor": false, "vina_hydrophobic": true}, + {"smarts": "[#35]", "vina_ri": 2.0, "vina_donor": false, "vina_acceptor": false, "vina_hydrophobic": true}, + {"smarts": "[#53]", "vina_ri": 2.2, "vina_donor": false, "vina_acceptor": false, "vina_hydrophobic": true}, + + {"smarts": "[#12]", "vina_ri": 1.2, "vina_donor": true, "vina_acceptor": false, "vina_hydrophobic": false}, + {"smarts": "[#20]", "vina_ri": 1.2, "vina_donor": true, "vina_acceptor": false, "vina_hydrophobic": false}, + {"smarts": "[#25]", "vina_ri": 1.2, "vina_donor": true, "vina_acceptor": false, "vina_hydrophobic": false}, + {"smarts": "[#26]", "vina_ri": 1.2, "vina_donor": true, "vina_acceptor": false, "vina_hydrophobic": false}, + {"smarts": "[#30]", "vina_ri": 1.2, "vina_donor": true, "vina_acceptor": false, "vina_hydrophobic": false}, + + {"smarts": "[#8][#1]", "vina_donor": true}, + {"smarts": "[#7][#1]", "vina_donor": true}, + {"smarts": "[#16][#1]", "vina_donor": true}, + {"smarts": "[#7X2v3]", "vina_acceptor": true}, + {"smarts": "[#7X3v3]", "vina_acceptor": true}, + {"smarts": "[#7;X1]", "vina_acceptor": true}, + {"smarts": "[#7X3v3][a]", "vina_acceptor": false}, + {"smarts": "[#7X3v3][#6X3v4]", "vina_acceptor": false}, + {"smarts": "[#16X1-]", "vina_acceptor": true}, + {"smarts": "[#6][#7,#8,#9,#15,#16,#17,#35,#53]", "vina_hydrophobic": false} + ] +} diff --git a/src/moldrug/meeko_light/flexibility.py b/src/moldrug/meeko_light/flexibility.py new file mode 100644 index 0000000..d2b3a6c --- /dev/null +++ b/src/moldrug/meeko_light/flexibility.py @@ -0,0 +1,435 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Meeko flexibility typer +# + +from copy import deepcopy +from typing import Tuple, List +from .utils import pdbutils + +from .molsetup import Bond + + +def _calc_max_weighted_depth( + model: dict, + seed_node: int, + bonds_to_break: Tuple[Tuple], + visited: List[int] = None, + depth: int = 0, +) -> int: + """ + Scores a flexibility model based on the depth of the nesting. The number of atoms in a rigid group increases the + weight of the bonds. + + Parameters + ---------- + model: dict + The flexibility model being moved over. + seed_node: int + Starting node index. + bonds_to_break: tuple[tuple] + An immutable list of bonds to break. + visited: list[int] + Nodes that have been visited. + depth: int + + Returns + ------- + max_value: int + The flexibility model score. + """ + glue_atoms = [] + # Adds the indices from bonds to break to the list of glue atoms + for i, j in bonds_to_break: + glue_atoms.append(i) + glue_atoms.append(j) + # Pulls the graph and members from the given flexibility model + graph = model["rigid_body_graph"] + members = model["rigid_body_members"] + if visited is None: + visited = [] + nr_atoms = len(members[seed_node]) + compensation = -1 # First atom after rotatable bond doesn't move, so doesn't weigh + # Atoms in breakable bonds count twice, because of added complexity + compensation += int(sum([i in glue_atoms for i in members[seed_node]])) + this_value = depth**2 * (nr_atoms + compensation) + max_value = this_value + visited.append(seed_node) + # Recurse on nodes that have not been visited yet. + for node in graph[seed_node]: + if node not in visited: + visited.append(node) + new_value = _calc_max_weighted_depth( + model, node, bonds_to_break, visited, depth + 1 + ) + max_value = max(max_value, new_value + this_value) + return max_value + + +def merge_terminal_atoms(flex_model: dict, not_terminal_atoms: List[int] = ()) -> None: + """ + Rotatable bonds that link to a rigid body group that contains one atom are removed because that one atom lies on the + bond axis and rotating the bond does not result in any movement of the atom. The atom after the removed rotatable + bond is merged with the rigid body group that is upstream of the removed bond. + + Parameters + ---------- + flex_model: dict + Flexibility model + not_terminal_atoms: list + A list of non-terminal atoms + + Returns + ------- + None + """ + # Get members and graph from the flexibility model + members_dict = flex_model["rigid_body_members"] + graph = flex_model["rigid_body_graph"] + remove = {} + # Loops over the members dict and adds members that can be deleted to the removal dict + for group_index, members in members_dict.items(): + if ( + len(members) == 1 + and len(graph[group_index]) == 1 + and members[0] not in not_terminal_atoms + and group_index != flex_model["root"] + ): + remove[group_index] = members[0] + # Goes through the removal dict and removes its contents from the flexibility model + for group_index, atom_index in remove.items(): + host_index = graph[group_index][0] + flex_model["rigid_body_members"].pop(group_index) + flex_model["rigid_body_members"][host_index].append(atom_index) + flex_model["rigid_body_graph"].pop(group_index) + flex_model["rigid_body_graph"][host_index].remove(group_index) + flex_model["rigid_body_connectivity"].pop((host_index, group_index)) + flex_model["rigid_body_connectivity"].pop((group_index, host_index)) + flex_model["rigid_index_by_atom"][atom_index] = host_index + # subtracts the number of members we are removing from the count + flex_model["rigid_body_count"] -= len(remove) + return + + +def get_flexibility_model( + molsetup, + root_atom_index: int = None, + break_combo_data: dict = None, +): + """ + Given a MoleculeSetup, creates a flexibility model for that MoleculeSetup. Breaks macrocyclic rings if bond break + information is provided. + + Parameters + ---------- + molsetup: RDKitMoleculeSetup + The molecule setup to generate a flexibility model for. + root_atom_index: int + + break_combo_data: dict + Data about different bond break combinations from the macrocycle typer. + + Returns + ------- + model: dict + The best flexibility model that could be generated given this MoleculeSetup and break data. + broken_bonds: list + A list of the bonds broken corresponding to the returned flexibility model. + + Raises + ------ + RuntimeError: + If the flexibility model generated ignores atoms that are not explicitly marked as needing to be ignored in the + MoleculeSetup. + """ + # If no macrocyclic rings are to be broken, we simply build the flexibility model. + if break_combo_data is None or len(break_combo_data["bond_break_combos"]) == 0: + bonds_to_break = () + unbroken_rings_bonds = [] + # Gets all bonds from the MoleculeSetup + for ring in molsetup.rings: + for bond in molsetup.get_bonds_in_ring(ring): + unbroken_rings_bonds.append(bond) + # Instantiates the flexibility model + flex_model = walk_rigid_body_graph( + molsetup, bonds_to_break, unbroken_rings_bonds + ) + # Gets the number of atoms that aren't marked to be ignored + nr_not_ignored = sum( + [not molsetup.get_is_ignore(i) for i in range(len(molsetup.atoms))] + ) + # Checks the validity of the model by ensuring that all of the atoms that should have been visited were visited, + # otherwise raises an error + if len(flex_model["visited"]) != nr_not_ignored: + molsetup.show() + msg = f"{len(flex_model['visited'])=} differs from not-ignored atoms {nr_not_ignored}" + raise RuntimeError(msg) + # Sets the model root body index + root_body_index = get_root_body_index(flex_model, root_atom_index) + flex_model["root"] = root_body_index + broken_bonds = [] + # Returns the model and the list of broken bonds + return flex_model, broken_bonds + + # The macrocycle typer enumerated rings to break, and a number of lists + # of bonds to break/delete that result in breaking of the enumerated rings. + # to be deleted to break each ring. Deleting a bond results in the remaining + # bonds within the same ring to become rotatable. + bond_break_combos = break_combo_data["bond_break_combos"] + bond_break_scores = break_combo_data["bond_break_scores"] + unbroken_rings_list = break_combo_data["unbroken_rings"] + best_model = None + best_score = float("+inf") + best_index = None + # Loops through the potential bond break combinations, generates flexibility model information for each one, and + # compares the flexibility models to determine the best model to return. + for index, _ in enumerate(bond_break_combos): + bond_break_combo = bond_break_combos[index] + bond_break_score = bond_break_scores[index] + unbroken_rings_bonds = [] + for ring in unbroken_rings_list[index]: + for bond in molsetup.get_bonds_in_ring(ring): + unbroken_rings_bonds.append(bond) + + # Generates a model + flex_model = walk_rigid_body_graph( + molsetup, bond_break_combo, unbroken_rings_bonds + ) + + # Gets the number of atoms that aren't marked to be ignored + nr_not_ignored = sum( + [not molsetup.get_is_ignore(i) for i in range(len(molsetup.atoms))] + ) + # Checks the validity of the model by ensuring that all of the atoms that should have been visited were visited, + # otherwise raises an error + if len(flex_model["visited"]) != nr_not_ignored: + msg = f"{len(flex_model['visited'])=} differs from not-ignored atoms {nr_not_ignored}" + raise RuntimeError(msg) + # Sets the model root body index + root_body_index = get_root_body_index(flex_model, root_atom_index) + flex_model["root"] = root_body_index + depth_weighted = _calc_max_weighted_depth( + flex_model, flex_model["root"], bond_break_combo + ) + # larger bond_break_score is better, larger depth is worse + # bond break score kinda disappeared in another branch (bonds are either breakable or not) + score = depth_weighted - 0.001 * bond_break_score + if score < best_score: + best_score = score + best_model = flex_model + best_index = index + + best_model["score"] = best_score + broken_bonds = list(break_combo_data["bond_break_combos"][best_index]) + # Returns the best model and the list of broken bonds + return best_model, broken_bonds + + +def get_root_body_index(model: dict, root_atom_index: int = None) -> int: + """ + Gets the index of the rigid body group in the flexibility model that contains the given root_atom_index + + Parameters + ---------- + model: dict + Flexibility model + root_atom_index: int + + + Returns + ------- + root_body_index: int + """ + + # find and return index of rigid body group that contains root_atom_index + if root_atom_index is not None: + for body_index in model["rigid_body_members"]: + if root_atom_index in model["rigid_body_members"][body_index]: + # 1-index atoms + return body_index + + # find rigid group that minimizes weighted graph depth + graph = deepcopy(model["rigid_body_graph"]) + while len(graph) > 2: # remove leaves until 1 or 2 rigid groups remain + leaves = [] + for vertex, edges in list(graph.items()): + if len(edges) == 1: + leaves.append(vertex) + for leaf in leaves: + for vertex, edges in list(graph.items()): + if leaf in edges: + edges.remove(leaf) + graph[vertex] = edges + del graph[leaf] + if len(graph) == 0: + root_body_index = 0 + elif len(graph) == 1: + root_body_index = list(graph.keys())[0] + else: + r1, r2 = list(graph.keys()) + r1_size = len(model["rigid_body_members"][r1]) + r2_size = len(model["rigid_body_members"][r2]) + if r1_size >= r2_size: + root_body_index = r1 + else: + root_body_index = r2 + return root_body_index + + +def update_closure_atoms( + molsetup, bonds_to_break: List[Tuple], glue_pseudo_atoms: dict +) -> None: + """ + Create pseudoatoms required by breaking bonds in the flexibility model + + Parameters + ---------- + molsetup: MoleculeSetup + Molecule Setup to modify and add pseudoatoms to + bonds_to_break: tuple + List of bonds to break + glue_pseudo_atoms: dict + + Returns + ------- + None + """ + # Loops through the bonds to break and adds all the bonds to + for i, bond in enumerate(bonds_to_break): + molsetup.ring_closure_info.bonds_removed.append( + bond + ) # bond is a pair of atom indices + + # calculate position and parameters of the pseudoatoms for the closure + for idx in (0, 1): + target = bond[1 - idx] + anchor = bond[0 - idx] + if glue_pseudo_atoms is None or len(glue_pseudo_atoms) == 0: + coord = molsetup.get_coord(target) + else: + coord = glue_pseudo_atoms[anchor] + anchor_info = molsetup.get_pdbinfo(anchor) + pdbinfo = pdbutils.PDBAtomInfo( + "G", + anchor_info.resName, + anchor_info.resNum, + anchor_info.icode, + anchor_info.chain, + ) + pseudo_index = molsetup.add_pseudoatom( + coord=coord, + charge=0.0, + anchor_list=[anchor], + atom_type=f"G{i}", + rotatable=False, + pdbinfo=pdbinfo, + ) + if anchor in molsetup.ring_closure_info.pseudos_by_atom: + raise RuntimeError("did not expect more than one G per atom") + molsetup.ring_closure_info.pseudos_by_atom[anchor] = pseudo_index + molsetup.set_atom_type(bond[0], "CG%d" % i) + molsetup.set_atom_type(bond[1], "CG%d" % i) + return + + +def walk_rigid_body_graph( + molsetup, + bonds_to_break: tuple, + unbroken_rings_bonds: List[Tuple], + start: int = None, + data: dict = None, +): + """ + Recursively walks through the MoleculeSetup to build a graph of rigid bodies. Uses that graph to create and + populate a flexibility model. + + Parameters + ---------- + molsetup: RDKitMoleculeSetup + MoleculeSetup to walk through + bonds_to_break: tuple[tuple] + unbroken_rings_bonds: list[tuple(int, int] + start: int + data: dict + + Returns + ------- + data: dict + A dictionary representing a flexibility model + """ + # If start is none, uses the default start which is the first non-ignored atom + if start is None: + for atom in molsetup.atoms: + if not atom.is_ignore: + start = atom.index + break + # If the start atom is marked to be ignored, returns nothing + if molsetup.get_is_ignore(start): + return + if data is None: + data = { + "visited": [], + "rigid_body_count": 0, + "rigid_index_by_atom": {}, + "rigid_body_members": {}, + "rigid_body_connectivity": {}, + "rigid_body_graph": {}, + } + data["visited"].append(start) + data["rigid_index_by_atom"][start] = data["rigid_body_count"] + rigid_index = data["rigid_body_count"] + sprouts_buffer = [] + # The while loop goes on until all atoms in `group_members` are queried + # for their neighbors. `idx` is the index within `group_members`, and + # `current` is the index of the atom being queried. When a neighbor is + # bonded to `current` by a non-rotatable/rigid bond, it is added to + # `group_members`. Remember that `idx` is not an index of an atom. + idx = 0 + group_members = [start] + while idx < len(group_members): + current = group_members[idx] + if molsetup.get_is_ignore(current): + idx += 1 + continue + for neigh in molsetup.get_neighbors(current): + if molsetup.get_is_ignore(neigh): + continue + bond_id = Bond.get_bond_id(current, neigh) + if bond_id in bonds_to_break: + continue + bond_info = None + if bond_id in molsetup.bond_info: + bond_info = molsetup.bond_info[bond_id] + if neigh in data["visited"]: + neigh_in_other_rigid_body = ( + rigid_index != data["rigid_index_by_atom"][neigh] + ) + if not bond_info.rotatable and neigh_in_other_rigid_body: + raise RuntimeError( + "Flexible bonds within rigid group. We have a problem." + ) + continue + if bond_info.rotatable and bond_id not in unbroken_rings_bonds: + sprouts_buffer.append((current, neigh)) + else: + group_members.append(neigh) + data["rigid_index_by_atom"][neigh] = rigid_index + data["visited"].append(neigh) + idx += 1 + data["rigid_body_members"][rigid_index] = list(group_members) + data["rigid_body_graph"].setdefault(rigid_index, []) + for current, neigh in sprouts_buffer: + if neigh in data["visited"]: + continue + data["rigid_body_count"] += 1 + next_rigid_index = data["rigid_body_count"] + data["rigid_body_connectivity"][rigid_index, next_rigid_index] = current, neigh + data["rigid_body_connectivity"][next_rigid_index, rigid_index] = neigh, current + data["rigid_body_graph"].setdefault(next_rigid_index, []) + data["rigid_body_graph"][rigid_index].append(next_rigid_index) + data["rigid_body_graph"][next_rigid_index].append(rigid_index) + walk_rigid_body_graph( + molsetup, bonds_to_break, unbroken_rings_bonds, neigh, data + ) + return data diff --git a/src/moldrug/meeko_light/macrocycle.py b/src/moldrug/meeko_light/macrocycle.py new file mode 100644 index 0000000..8f74f10 --- /dev/null +++ b/src/moldrug/meeko_light/macrocycle.py @@ -0,0 +1,319 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Meeko macrocycle builder +# +from typing import List, Set, Tuple + +from .molsetup import Bond + +# region +DEFAULT_MIN_RING_SIZE = 7 +DEFAULT_MAX_RING_SIZE = 33 +DEFAULT_DOUBLE_BOND_PENALTY = 50 +DEFAULT_MAX_BREAKS = 4 +# endregion + + +class FlexMacrocycle: + """ + Attributes + ---------- + _min_ring_size: int + _max_ring_size: int + _double_bond_penalty: float + max_breaks: int + setup: + breakable_rings: + """ + + def __init__( + self, + min_ring_size: int = DEFAULT_MIN_RING_SIZE, + max_ring_size: int = DEFAULT_MAX_RING_SIZE, + double_bond_penalty: float = DEFAULT_DOUBLE_BOND_PENALTY, + max_breaks: int = DEFAULT_MAX_BREAKS, + allow_break_atype_A: bool = False, + ): + """ + Initialize macrocycle typer. + + Parameters + ---------- + min_ring_size: int + Minimum size of the ring, default is 7. + max_ring_size: int + Maximum size of the ring, default is 33. + double_bond_penalty: float + max_breaks: int + allow_break_type_A: bool + Allow breaking bonds involving atoms typed A, default is False. + """ + self._min_ring_size = min_ring_size + self._max_ring_size = max_ring_size + # accept also double bonds (if nothing better is found) + self._double_bond_penalty = double_bond_penalty + self.max_breaks = max_breaks + self.allow_break_atype_A = allow_break_atype_A + + self.setup = None + self.breakable_rings = None + + def collect_rings(self, setup): + """ + Gets non-aromatic rings of desired size and lists bonds that are part of unbreakable rings. Bonds belonging to + rigid cycles can't be deleted or made rotatable even if they are part of a breakable ring. + + Parameters + ---------- + setup: RDKitMoleculeSetup + + Returns + ------- + breakable_rings: list + A list of breakable ring ids + bonds_in_rigid_cycles: set + A set of the bonds in rigid cycles + """ + breakable_rings = [] + rigid_rings = [] + for ring_id in list( + setup.rings.keys() + ): # ring_id are the atom indices in each ring + size = len(ring_id) + if setup.rings[ring_id].is_aromatic: + rigid_rings.append(ring_id) + elif size < self._min_ring_size: + rigid_rings.append(ring_id) + # do not add rings > _max_ring_size to rigid_rings + # because bonds in rigid rings will not be breakable + # and these bonds may also belong to breakable rings + elif size <= self._max_ring_size: + breakable_rings.append(ring_id) + + bonds_in_rigid_cycles = set() + for ring_atom_indices in rigid_rings: + for bond in setup.get_bonds_in_ring(ring_atom_indices): + bonds_in_rigid_cycles.add(bond) + + return breakable_rings, bonds_in_rigid_cycles + + def _score_bond(self, bond: Tuple[int]) -> int: + """ + Calculates a score for the likeliness that a bond will be broken. + + Parameters + ---------- + bond: tuple[int, int] + Input bond to score + + Returns + ------- + score: int + A score for the bond. + """ + bond = Bond.get_bond_id(bond[0], bond[1]) + if not self.setup.bond_info[bond].rotatable: + return -1 + atom_idx1, atom_idx2 = bond + for i in (atom_idx1, atom_idx2): + atype = self.setup.get_atom_type(i) + is_allowed_A = self.allow_break_atype_A and atype == "A" + if atype != "C" and not is_allowed_A: + return -1 + # historically we returned a score <= 100, that was lower for triple + # bonds, chiral atoms, conjugated bonds, and double bonds. This score + # gets combined with the graph depth score in flexibility.py that + # is lower when more consecutive torsions exist in a single "branch" + # of the torsion tree. Any positive number can be returned here. + return 100 + + def get_breakable_bonds(self, bonds_in_rigid_rings): + """ + Find breaking points for rings following the guidelines defined in [1]. + The optimal bond has the following properties: + - does not involve a chiral atom + - is not double/triple (?) + - is between two carbons + (preferably? we can now generate pseudoAtoms on the fly!) + - is a bond present only in one ring + + [1] Forli, Botta, J. Chem. Inf. Model., 2007, 47 (4) + DOI: 10.1021/ci700036j + + Parameters + ---------- + bonds_in_rigid_rings: set[tuple] + A set of bonds in rigid rings. + + Returns + ------- + breakable: dict + A dictionary of mapping breakable bonds to bond scores + """ + breakable = {} + for ring_atom_indices in self.breakable_rings: + for bond in self.setup.get_bonds_in_ring(ring_atom_indices): + score = self._score_bond(bond) + if score > 0 and bond not in bonds_in_rigid_rings: + breakable[bond] = score + return breakable + + def search_macrocycle(self, setup, delete_these_bonds: List = None): + """ + Search for macrocycles in the molecule + + Parameters + ---------- + setup: RDKitMoleculeSetup + delete_these_bonds: list + + Returns + ------- + break_combo_data: + bonds_in_rigid_rings: + """ + + if delete_these_bonds is None: + delete_these_bonds = [] + + self.setup = setup + + self.breakable_rings, bonds_in_rigid_rings = self.collect_rings(setup) + if len(delete_these_bonds) == 0: + breakable_bonds = self.get_breakable_bonds(bonds_in_rigid_rings) + else: + breakable_bonds = {} + for bond in delete_these_bonds: + bond = Bond.get_bond_id(bond[0], bond[1]) + breakable_bonds[bond] = self._score_bond(bond) + break_combo_data = self.combinatorial_break_search(breakable_bonds) + return break_combo_data, bonds_in_rigid_rings + + def combinatorial_break_search(self, breakable_bonds): + """ + Enumerate all combinations of broken bonds. Once a bond is broken, it will break one or more rings. Subsequent + bonds will be pulled from intact (unbroken) rings. The number of broken bonds may be variable. + Returns only combinations of broken bonds that break the maximum number of broken bonds. + + Parameters + ---------- + breakable_bonds: dict + A dictionary mapping breakable bonds to the bond score. + + Returns + ------- + break_combo_data: dict + A dictionary containing information about possible bond break combinations + """ + + max_breaks = self.max_breaks + break_combos = self._recursive_break( + self.breakable_rings, max_breaks, breakable_bonds, set(), [] + ) + break_combos = list(break_combos) # convert from set + max_broken_bonds = 0 + output_break_combos = [] # found new max, discard prior data + output_bond_scores = [] + output_unbroken_rings = [] + for broken_bonds in break_combos: + n_broken_bonds = len(broken_bonds) + bond_score = sum([breakable_bonds[bond] for bond in broken_bonds]) + if n_broken_bonds > max_broken_bonds: + max_broken_bonds = n_broken_bonds + output_break_combos = [] # found new max, discard prior data + output_bond_scores = [] + output_unbroken_rings = [] + if n_broken_bonds == max_broken_bonds: + output_break_combos.append(broken_bonds) + output_bond_scores.append(bond_score) + u = self.get_unbroken_rings(self.breakable_rings, broken_bonds) + output_unbroken_rings.append(u) + + break_combo_data = { + "bond_break_combos": output_break_combos, + "bond_break_scores": output_bond_scores, + "unbroken_rings": output_unbroken_rings, + } + return break_combo_data + + def _recursive_break( + self, rings, max_breaks, breakable_bonds, output: Set = None, broken_bonds: List = None + ): + """ + + Parameters + ---------- + rings: list + List of rings to check + max_breaks: int + Maxumimum number of breaks allowed + breakable_bonds: + output: set + broken_bonds: + + Returns + ------- + output: set + + """ + if output is None: + output = set() + if broken_bonds is None: + broken_bonds = [] + if max_breaks == 0: + return output + unbroken_rings = self.get_unbroken_rings(rings, broken_bonds) + atoms_in_broken_bonds = [ + atom_idx for bond in broken_bonds for atom_idx in bond + ] + for bond in breakable_bonds: + if bond[0] in atoms_in_broken_bonds or bond[1] in atoms_in_broken_bonds: + continue # each atom can be in only one broken bond + is_bond_in_ring = False + for ring in unbroken_rings: + if bond in self.setup.get_bonds_in_ring(ring): + is_bond_in_ring = True + break + if is_bond_in_ring: + current_broken_bonds = [(a, b) for (a, b) in broken_bonds + [bond]] + num_unbroken_rings = len( + self.get_unbroken_rings(rings, current_broken_bonds) + ) + data_row = tuple(sorted([(a, b) for (a, b) in current_broken_bonds])) + output.add(data_row) + if num_unbroken_rings > 0: + output = self._recursive_break( + rings, + max_breaks - 1, + breakable_bonds, + output, + current_broken_bonds, + ) + return output + + def get_unbroken_rings(self, rings, broken_bonds): + """ + + Parameters + ---------- + rings: + broken_bonds: + + Returns + ------- + unbroken: list + List of unbroken rings. + """ + unbroken = [] + for ring in rings: + is_unbroken = True + for bond in broken_bonds: + if bond in self.setup.get_bonds_in_ring( + ring + ): # consider precalculating bonds + is_unbroken = False + break # pun intended + if is_unbroken: + unbroken.append(ring) + return unbroken diff --git a/src/moldrug/meeko_light/molecule_pdbqt.py b/src/moldrug/meeko_light/molecule_pdbqt.py new file mode 100644 index 0000000..cd938f9 --- /dev/null +++ b/src/moldrug/meeko_light/molecule_pdbqt.py @@ -0,0 +1,690 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Meeko +# + +import os +from collections import defaultdict + +import numpy as np +from scipy import spatial + +from .utils.covalent_radius_table import covalent_radius +from .utils.autodock4_atom_types_elements import autodock4_atom_types_elements + + +atom_property_definitions = {'H': 'vdw', 'C': 'vdw', 'A': 'vdw', 'N': 'vdw', 'P': 'vdw', 'S': 'vdw', + 'Br': 'vdw', 'I': 'vdw', 'F': 'vdw', 'Cl': 'vdw', + 'NA': 'hb_acc', 'OA': 'hb_acc', 'SA': 'hb_acc', 'OS': 'hb_acc', 'NS': 'hb_acc', + 'HD': 'hb_don', 'HS': 'hb_don', + 'Mg': 'metal', 'Ca': 'metal', 'Fe': 'metal', 'Zn': 'metal', 'Mn': 'metal', + 'MG': 'metal', 'CA': 'metal', 'FE': 'metal', 'ZN': 'metal', 'MN': 'metal', + 'W': 'water', + 'G0': 'glue', 'G1': 'glue', 'G2': 'glue', 'G3': 'glue', + 'CG0': 'glue', 'CG1': 'glue', 'CG2': 'glue', 'CG3': 'glue'} + + +def _read_ligand_pdbqt_file(pdbqt_string, poses_to_read=-1, energy_range=-1, is_dlg=False, skip_typing=False): + i = 0 + n_poses = 0 + previous_serial = 0 + tmp_positions = [] + tmp_atoms = [] + tmp_actives = [] + tmp_pdbqt_string = '' + water_indices = {*()} + location = 'ligand' + energy_best_pose = None + is_first_pose = True + is_model = False + mol_index = -1 # incremented for each ROOT keyword + atoms_dtype = [('idx', 'i4'), ('serial', 'i4'), ('name', 'U4'), ('resid', 'i4'), + ('resname', 'U3'), ('chain', 'U1'), ('xyz', 'f4', (3)), + ('partial_charges', 'f4'), ('atom_type', 'U3')] + + atoms = None + positions = [] + + # flexible_residue is for atoms between BEGIN_RES and END_RES keywords, ligand otherwise. + # flexres assigned "ligand" if BEGIN/END RES keywords are missing + # mol_index distinguishes different ligands and flexres because ROOT keyword increments mol_index + atom_annotations = {'ligand': [], 'flexible_residue': [], 'water': [], + 'hb_acc': [], 'hb_don': [], + 'all': [], 'vdw': [], + 'glue': [], 'reactive': [], 'metal': [], + 'mol_index': {}, + } + pose_data = { + 'n_poses': None, + 'active_atoms': [], + 'free_energies': [], + 'intermolecular_energies': [], + 'internal_energies': [], + 'index_map': {}, + 'pdbqt_string': [], + 'smiles': {}, + 'smiles_index_map': {}, + 'smiles_h_parent': {}, + 'cluster_id': [], + 'rank_in_cluster': [], + 'cluster_leads_sorted': [], + 'cluster_size': [], + "mol_index_to_flexible_residue": {}, + } + + tmp_cluster_data = {} + + buffer_index_map = {} + buffer_smiles = None + buffer_smiles_index_map = [] + buffer_smiles_h_parent = [] + buffer_flexres_id = None + + lines = pdbqt_string.split('\n') + if len(lines[-1]) == 0: + lines = lines[:-1] + lines = [line + '\n' for line in lines] + for line in lines: + if is_dlg: + if line.startswith('DOCKED'): + line = line[8:] + # parse clustering + elif line.endswith('RANKING\n'): + fields = line.split() + cluster_id = int(fields[0]) + subrank = int(fields[1]) + run_id = int(fields[2]) + tmp_cluster_data[run_id] = (cluster_id, subrank) + else: + continue + + if not line.startswith(('MODEL', 'ENDMDL')): + # This is very lazy I know... + # But would you rather spend time on rebuilding the whole torsion tree and stuff + # for writing PDBQT files or drinking margarita? Energy was already spend to build + # that, so let's re-use it! + tmp_pdbqt_string += line + + if line.startswith('MODEL'): + # Reinitialize variables + i = 0 + previous_serial = 0 + tmp_positions = [] + tmp_atoms = [] + tmp_actives = [] + tmp_pdbqt_string = '' + is_model = True + mol_index = -1 # incremented for each ROOT keyword + elif line.startswith('ATOM') or line.startswith("HETATM"): + serial = int(line[6:11].strip()) + name = line[12:16].strip() + resname = line[17:20].strip() + chainid = line[21].strip() + resid = int(line[22:26].strip()) + xyz = np.array([line[30:38].strip(), line[38:46].strip(), line[46:54].strip()], dtype=float) + try: + # PDBQT files from dry.py script are stripped from their partial charges. sigh... + partial_charges = float(line[70:76].strip()) + except: + partial_charges = 0.0 + atom_type = line[77:-1].strip() + + # We are looking for gap in the serial atom numbers. Usually if they + # are not following it means that atoms are missing. This will happen with + # water molecules after using dry.py, only non-overlapping water molecules + # are kept. Also if the current serial becomes suddenly inferior than the + # previous and equal to 1, it means that we are now in another molecule/flexible + # residue. So here we are adding dummy atoms + + if (previous_serial + 1 != serial) and not (serial < previous_serial and serial == 1): + diff = serial - previous_serial - 1 + for _ in range(diff): + xyz_nan = [999.999, 999.999, 999.999] + tmp_atoms.append((i, 9999, 'XXXX', 9999, 'XXX', 'X', xyz_nan, 999.999, 'XX')) + tmp_positions.append(xyz_nan) + i += 1 + + # Once it is done, we can return to a normal life... and add existing atoms + tmp_atoms.append((i, serial, name, resid, resname, chainid, xyz, partial_charges, atom_type)) + tmp_positions.append(xyz) + tmp_actives.append(i) + + if is_first_pose: + atom_annotations["mol_index"].setdefault(mol_index, []) + atom_annotations["mol_index"][mol_index].append(i) + # We store water idx separately from the rest since their number can be variable + if atom_type != 'W': + atom_annotations[location].append(i) + atom_annotations['all'].append(i) + if not skip_typing: + atom_annotations[atom_property_definitions[atom_type]].append(i) + + if atom_type == 'W': + water_indices.update([i]) + + previous_serial = serial + i += 1 + elif line.startswith("ROOT") and is_first_pose: + mol_index += 1 + # buffers needed because REMARKS preceeds ROOT + pose_data["index_map"][mol_index] = buffer_index_map + pose_data["smiles"][mol_index] = buffer_smiles + pose_data["smiles_index_map"][mol_index] = buffer_smiles_index_map + pose_data["smiles_h_parent"][mol_index] = buffer_smiles_h_parent + pose_data["mol_index_to_flexible_residue"][mol_index] = buffer_flexres_id + buffer_index_map = {} + buffer_smiles = None + buffer_smiles_index_map = [] + buffer_smiles_h_parent = [] + buffer_flexres_id = None + elif line.startswith('REMARK INDEX MAP') and is_first_pose: + integers = [int(integer) for integer in line.split()[3:]] + if len(integers) % 2 == 1: + raise RuntimeError("Number of indices in INDEX MAP is odd") + for j in range(int(len(integers) / 2)): + buffer_index_map[integers[j*2]] = integers[j*2 + 1] + elif line.startswith('REMARK SMILES IDX') and is_first_pose: + integers = [int(integer) for integer in line.split()[3:]] + if len(integers) % 2 == 1: + raise RuntimeError("Number of indices in SMILES IDX is odd") + buffer_smiles_index_map.extend(integers) + elif line.startswith('REMARK H PARENT') and is_first_pose: + integers = [int(integer) for integer in line.split()[3:]] + if len(integers) % 2 == 1: + raise RuntimeError("Number of indices in H PARENT is odd") + buffer_smiles_h_parent.extend(integers) + elif line.startswith('REMARK SMILES') and is_first_pose: # must check after SMILES IDX + buffer_smiles = line.split()[2] + elif line.startswith('REMARK VINA RESULT') or line.startswith('USER Estimated Free Energy of Binding ='): + # Read free energy from output PDBQT files + try: + # Vina + energy = float(line.split()[3]) + except: + # AD4 + energy = float(line[45:].split()[0]) # no guarantee of space between = and number + + if energy_best_pose is None: + energy_best_pose = energy + energy_current_pose = energy + + # NOTE this assumes poses are sorted by increasing energy + diff_energy = energy_current_pose - energy_best_pose + if (energy_range <= diff_energy and energy_range != -1): + break + + pose_data['free_energies'].append(energy) + elif not is_dlg and line.startswith('REMARK INTER:'): + pose_data['intermolecular_energies'].append(float(line.split()[2])) + elif not is_dlg and line.startswith('REMARK INTRA:'): + pose_data['internal_energies'].append(float(line.split()[2])) + elif is_dlg and line.startswith('USER (1) Final Intermolecular Energy ='): + pose_data['intermolecular_energies'].append(float(line[45:].split()[0])) + elif is_dlg and line.startswith('USER (2) Final Total Internal Energy ='): + pose_data['internal_energies'].append(float(line[45:].split()[0])) + elif line.startswith('BEGIN_RES'): + location = 'flexible_residue' + buffer_flexres_id = " ".join(line.strip().split()[1:]) + elif line.startswith('END_RES'): + # We never know if there is a molecule just after the flexible residue... + location = 'ligand' + elif line.startswith('ENDMDL'): + n_poses += 1 + # After reading the first pose no need to store atom properties + # anymore, it is the same for every pose + is_first_pose = False + + tmp_atoms = np.array(tmp_atoms, dtype=atoms_dtype) + + if atoms is None: + # We store the atoms (topology) only once, since it is supposed to be + # the same for all the molecules in the PDBQT file (except when water molecules + # are involved... classic). But we will continue to compare the topology of + # the current pose with the first one seen in the PDBQT file, to be sure only + # the atom positions are changing. + atoms = tmp_atoms.copy() + else: + # Check if the molecule topology is the same for each pose + # We ignore water molecules (W) and atom type XX + columns = ['idx', 'serial', 'name', 'resid', 'resname', 'chain', 'partial_charges', 'atom_type'] + topology1 = atoms[np.isin(atoms['atom_type'], ['W', 'XX'], invert=True)][columns] + topology2 = tmp_atoms[np.isin(atoms['atom_type'], ['W', 'XX'], invert=True)][columns] + + if not np.array_equal(topology1, topology2): + error_msg = 'molecules have different topologies' + raise RuntimeError(error_msg) + + # Update information about water molecules (W) as soon as we find new ones + tmp_water_molecules_idx = tmp_atoms[tmp_atoms['atom_type'] == 'W']['idx'] + water_molecules_idx = atoms[atoms['atom_type'] == 'XX']['idx'] + new_water_molecules_idx = list(set(tmp_water_molecules_idx).intersection(water_molecules_idx)) + atoms[new_water_molecules_idx] = tmp_atoms[new_water_molecules_idx] + + positions.append(tmp_positions) + pose_data['active_atoms'].append(tmp_actives) + pose_data['pdbqt_string'].append(tmp_pdbqt_string) + + if (n_poses >= poses_to_read and poses_to_read != -1): + break + + # if here is only one molecule + # so when we reach the end of the file, we store the atoms, + # positions and actives stuff. + if not is_model: + n_poses += 1 + atoms = np.array(tmp_atoms, dtype=atoms_dtype) + positions.append(tmp_positions) + pose_data['active_atoms'].append(tmp_actives) + pose_data['pdbqt_string'].append(tmp_pdbqt_string) + + positions = np.array(positions).reshape((n_poses, atoms.shape[0], 3)) + + pose_data['n_poses'] = n_poses + + # We add indices of all the water molecules we saw + if water_indices: + atom_annotations['water'] = list(water_indices) + + # clustering + if len(tmp_cluster_data) > 0: + if len(tmp_cluster_data) != n_poses: + raise RuntimeError("Nr of poses in cluster data (%d) differs from nr of poses (%d)" % (len(tmp_cluster_data, n_poses))) + pose_data["cluster_id"] = [None] * n_poses + pose_data["rank_in_cluster"] = [None] * n_poses + pose_data["cluster_size"] = [None] * n_poses + cluster_ids = [cluster_id for _, (cluster_id, _) in tmp_cluster_data.items()] + n_clusters = max(cluster_ids) + pose_data["cluster_leads_sorted"] = [None] * n_clusters + for pose_index, (cluster_id, rank_in_cluster) in tmp_cluster_data.items(): + pose_data["cluster_id"][pose_index - 1] = cluster_id + pose_data["rank_in_cluster"][pose_index - 1] = rank_in_cluster + pose_data["cluster_size"][pose_index - 1] = cluster_ids.count(cluster_id) + if rank_in_cluster == 1: # is cluster lead + pose_data["cluster_leads_sorted"][cluster_id - 1] = pose_index - 1 + return atoms, positions, atom_annotations, pose_data + + +def _identify_bonds(atom_idx, positions, atom_types): + bonds = defaultdict(list) + KDTree = spatial.cKDTree(positions) + bond_allowance_factor = 1.1 + # If we ask more than the number of coordinates/element + # in the BHTree, we will end up with some inf values + k = 5 if len(atom_idx) > 5 else len(atom_idx) + atom_idx = np.array(atom_idx) + + # If there is only one atom, we know there won't be a single bond.. + if len(atom_idx) == 1: + return bonds + + for atom_i, position, atom_type in zip(atom_idx, positions, atom_types): + distances, indices = KDTree.query(position, k=k) + r_cov = covalent_radius[autodock4_atom_types_elements[atom_type]] + + optimal_distances = [bond_allowance_factor * (r_cov + covalent_radius[autodock4_atom_types_elements[atom_types[i]]]) for i in indices[1:]] + bonds[atom_i] = atom_idx[indices[1:][np.where(distances[1:] < optimal_distances)]].tolist() + + return bonds + + +class PDBQTMolecule: + + def __init__(self, pdbqt_string, name=None, poses_to_read=None, energy_range=None, is_dlg=False, skip_typing=False): + """PDBQTMolecule class for reading PDBQT (or dlg) files from AutoDock4, AutoDock-GPU or AutoDock-Vina + + Contains both __getitem__ and __iter__ methods, someone might lose his mind because of this. + + Args: + pdbqt_string (str): pdbqt string + name (str): name of the molecule (default: None, use filename without pdbqt suffix) + poses_to_read (int): total number of poses to read (default: None, read all) + energy_range (float): read docked poses until the maximum energy difference + from best pose is reach, for example 2.5 kcal/mol (default: Non, read all) + is_dlg (bool): input file is in dlg (AutoDock docking log) format (default: False) + skip_typing (bool, optional): Flag indicating that atomtyping should be skipped + """ + self._current_pose = 0 + self._pdbqt_filename = None + self._atoms = None + self._positions = None + self._bonds = None + self._atom_annotations = None + self._pose_data = None + self._name = name + + # Juice all the information from that PDBQT file + poses_to_read = poses_to_read if poses_to_read is not None else -1 + energy_range = energy_range if energy_range is not None else -1 + results = _read_ligand_pdbqt_file(pdbqt_string, poses_to_read, energy_range, is_dlg, skip_typing) + self._atoms, self._positions, self._atom_annotations, self._pose_data = results + + if self._atoms.shape[0] == 0: + raise RuntimeError('read 0 atoms. Consider PDBQTMolecule.from_file(fname)') + + # Build KDTrees for each pose (search closest atoms by distance) + self._KDTrees = [spatial.cKDTree(positions) for positions in self._positions] + + # Identify bonds in the ligands + if not skip_typing: + mol_atoms = self._atoms[self._atom_annotations['ligand']] + self._bonds = _identify_bonds(self._atom_annotations['ligand'], mol_atoms['xyz'], mol_atoms['atom_type']) + + """... then in the flexible residues + Since we are extracting bonds from docked poses, we might be in the situation + where the ligand reacted with one the flexible residues and we don't want to + consider them as normally bonded...""" + if self.has_flexible_residues(): + flex_atoms = self._atoms[self._atom_annotations['flexible_residue']] + self._bonds.update(_identify_bonds(self._atom_annotations['flexible_residue'], flex_atoms['xyz'], flex_atoms['atom_type'])) + + @classmethod + def from_file(cls, pdbqt_filename, name=None, poses_to_read=None, energy_range=None, is_dlg=False, skip_typing=False): + if name is None: + name = os.path.splitext(os.path.basename(pdbqt_filename))[0] + with open(pdbqt_filename) as f: + pdbqt_string = f.read() + instance = cls(pdbqt_string, name, poses_to_read, energy_range, is_dlg, skip_typing) + instance._pdbqt_filename = pdbqt_filename + return instance + + def __getitem__(self, value): + if isinstance(value, int): + if value < 0 or value >= self._positions.shape[0]: + raise IndexError('The index (%d) is out of range.' % value) + elif isinstance(value, slice): + raise TypeError('Slicing is not implemented for PDBQTMolecule object.') + else: + raise TypeError('Invalid argument type.') + + self._current_pose = value + return self + + def __iter__(self): + self._current_pose = -1 + return self + + def __next__(self): + if self._current_pose + 1 >= self._positions.shape[0]: + raise StopIteration + + self._current_pose += 1 + + return self + + def __repr__(self): + repr_str = '' + return (repr_str % (self._name, self._pose_data['n_poses'], self._atoms.shape[0])) + + @property + def name(self): + """Return the name of the molecule.""" + return self._name + + @property + def pose_id(self): + """Return the index of the current pose.""" + return self._current_pose + + @property + def score(self): + """Return the score (kcal/mol) of the current pose.""" + return self._pose_data['free_energies'][self._current_pose] + + def available_atom_properties(self, ignore_properties=None): + """Return all the available atom properties for that molecule. + + The following properties are ignored: ligand and flexible_residue + + """ + if ignore_properties is None: + ignore_properties = [] + + if not isinstance(ignore_properties, (list, tuple)): + ignore_properties = [ignore_properties] + + ignore_properties += ['ligand', 'flexible_residue', 'water'] + + return [k for k, v in self._atom_annotations.items() + if k not in ignore_properties and len(v) > 0] + + def has_flexible_residues(self): + """Tell if the molecule contains a flexible residue or not. + + Returns: + bool: True if contains flexible residues, otherwise False + + """ + if self._atom_annotations['flexible_residue']: + return True + else: + return False + + def has_water_molecules(self): + """Tell if the molecules contains water molecules or not in the current pose. + + Returns: + bool: True if contains water molecules in the current pose, otherwise False + + """ + active_atoms_idx = self._pose_data['active_atoms'][self._current_pose] + if set(self._atom_annotations['water']).intersection(active_atoms_idx): + return True + else: + return False + + def atoms(self, atom_idx=None, only_active=True): + """Return the atom i + + Args: + atom_idx (int, list): index of one or multiple atoms (0-based) + only_active (bool): return only active atoms (default: True, return only active atoms) + + Returns: + ndarray: 2d ndarray (atom_id, atom_name, resname, resid, chainid, xyz, q, t) + + """ + if atom_idx is not None: + if not isinstance(atom_idx, (list, tuple, np.ndarray)): + atom_idx = np.array(atom_idx, dtype=np.int) + else: + atom_idx = np.arange(0, self._atoms.shape[0]) + + # Get index of only the active atoms + if only_active: + active_atoms_idx = self._pose_data['active_atoms'][self._current_pose] + atom_idx = sorted(list(set(atom_idx).intersection(active_atoms_idx))) + + atoms = self._atoms[atom_idx].copy() + atoms['xyz'] = self._positions[self._current_pose, atom_idx, :] + + return atoms + + def positions(self, atom_idx=None, only_active=True): + """Return coordinates (xyz) of all atoms or a certain atom + + Args: + atom_idx (int, list): index of one or multiple atoms (0-based) + only_active (bool): return only active atoms (default: True, return only active atoms) + + Returns: + ndarray: 2d ndarray of coordinates (xyz) + + """ + return np.atleast_2d(self.atoms(atom_idx, only_active)['xyz']) + + def atoms_by_properties(self, atom_properties, only_active=True): + """Return atom based on their properties + + Args: + atom_properties (str or list): property of the atoms to retrieve + (properties: ligand, flexible_residue, vdw, hb_don, hb_acc, metal, water, reactive, glue) + only_active (bool): return only active atoms (default: True, return only active atoms) + + """ + if not isinstance(atom_properties, (list, tuple)): + atom_properties = [atom_properties] + + if len(atom_properties) > 1: + try: + atom_idx = set(self._atom_annotations[atom_properties[0]]) + + for atom_property in atom_properties[1:]: + atom_idx.intersection_update(self._atom_annotations[atom_property]) + except: + error_msg = 'Atom property %s is not valid. Valid atom properties are: %s' + raise KeyError(error_msg % (atom_property, self._atom_annotations.keys())) + + atom_idx = list(atom_idx) + else: + try: + atom_idx = self._atom_annotations[atom_properties[0]] + except: + error_msg = 'Atom property %s is not valid. Valid atom properties are: %s' + raise KeyError(error_msg % (atom_properties[0], self._atom_annotations.keys())) + + if atom_idx: + return self.atoms(atom_idx, only_active) + else: + return np.array([]) + + def closest_atoms_from_positions(self, xyz, radius, atom_properties=None, ignore=None): + """Retrieve indices of the closest atoms around a positions/coordinates + at a certain radius. + + Args: + xyz (np.ndarray): array of 3D coordinates + raidus (float): radius + atom_properties (str): property of the atoms to retrieve + (properties: ligand, flexible_residue, vdw, hb_don, hb_acc, metal, water, reactive, glue) + ignore (int or list): ignore atom for the search using atom id (0-based) + + Returns: + ndarray: 2d ndarray (atom_id, atom_name, resname, resid, chainid, xyz, q, t) + + """ + atom_idx = self._KDTrees[self._current_pose].query_ball_point(xyz, radius, p=2, return_sorted=True) + + # When nothing was found around... + if not atom_idx: + return np.array([]) + + # Handle the case when positions for of only one atom was passed in the input + try: + atom_idx = {i for j in atom_idx for i in j} + except: + atom_idx = set(atom_idx) + + if atom_properties is not None: + if not isinstance(atom_properties, (list, tuple)): + atom_properties = [atom_properties] + + try: + for atom_property in atom_properties: + atom_idx.intersection_update(self._atom_annotations[atom_property]) + except: + error_msg = 'Atom property %s is not valid. Valid atom properties are: %s' + raise KeyError(error_msg % (atom_property, self._atom_annotations.keys())) + + if ignore is not None: + if not isinstance(ignore, (list, tuple, np.ndarray)): + ignore = [ignore] + atom_idx = atom_idx.difference([i for i in ignore]) + + # Get index of only the active atoms + active_atoms_idx = self._pose_data['active_atoms'][self._current_pose] + atom_idx = list(set(atom_idx).intersection(active_atoms_idx)) + + if atom_idx: + atoms = self._atoms[atom_idx].copy() + atoms['xyz'] = self._positions[self._current_pose, atom_idx, :] + return atoms + else: + return np.array([]) + + def closest_atoms(self, atom_idx, radius, atom_properties=None): + """Retrieve indices of the closest atoms around a positions/coordinates + at a certain radius. + + Args: + atom_idx (int, list): index of one or multiple atoms (0-based) + raidus (float): radius + atom_properties (str or list): property of the atoms to retrieve + (properties: ligand, flexible_residue, vdw, hb_don, hb_acc, metal, water, reactive, glue) + + Returns: + ndarray: ndarray (atom_id, atom_name, resname, resid, chainid, xyz, q, t) + + """ + if not isinstance(atom_idx, (list, tuple)): + atom_idx = [atom_idx] + + # Get index of only the active atoms + active_atoms_idx = self._pose_data['active_atoms'][self._current_pose] + atom_idx = list(set(atom_idx).intersection(active_atoms_idx)) + + if atom_idx: + positions = self._positions[self._current_pose, atom_idx, :] + return self.closest_atoms_from_positions(positions, radius, atom_properties, atom_idx) + else: + return np.array([]) + + def neighbor_atoms(self, atom_idx): + """Return neighbor (bonded) atoms + + Args: + atom_idx (int, list): index of one or multiple atoms (0-based) + + Returns: + list_of_list: list of lists containing the neighbor (bonded) atoms (0-based) + + """ + if not isinstance(atom_idx, (list, tuple, np.ndarray)): + atom_idx = [atom_idx] + + # Get index of only the active atoms + active_atoms_idx = self._pose_data['active_atoms'][self._current_pose] + atom_idx = list(set(atom_idx).intersection(active_atoms_idx)) + + return [self._bonds[i] for i in atom_idx] + + def write_pdbqt_string(self, as_model=True): + """Write PDBQT output string of the current pose + + Args: + as_model (bool): Qdd MODEL/ENDMDL keywords to the output PDBQT string (default: True) + + Returns: + string: Description + + """ + if as_model: + pdbqt_string = 'MODEL %5d\n' % (self._current_pose + 1) + pdbqt_string += self._pose_data['pdbqt_string'][self._current_pose] + pdbqt_string += 'ENDMDL\n' + return pdbqt_string + else: + return self._pose_data['pdbqt_string'][self._current_pose] + + def write_pdbqt_file(self, output_pdbqtfilename, overwrite=False, as_model=False): + """Write PDBQT file of the current pose + + Args: + output_pdbqtfilename (str): filename of the output PDBQT file + overwrite (bool): overwrite on existing pdbqt file (default: False) + as_model (bool): Qdd MODEL/ENDMDL keywords to the output PDBQT string (default: False) + + """ + print(overwrite and os.path.isfile(output_pdbqtfilename)) + if not overwrite and os.path.isfile(output_pdbqtfilename): + raise RuntimeError('Output PDBQT file %s already exists' % output_pdbqtfilename) + + if as_model: + pdbqt_string = 'MODEL %5d\n' % (self._current_pose + 1) + pdbqt_string += self._pose_data['pdbqt_string'][self._current_pose] + pdbqt_string += 'ENDMDL\n' + else: + pdbqt_string = self._pose_data['pdbqt_string'][self._current_pose] + + with open(output_pdbqtfilename, 'w') as w: + w.write(pdbqt_string) diff --git a/src/moldrug/meeko_light/molsetup.py b/src/moldrug/meeko_light/molsetup.py new file mode 100644 index 0000000..ca717ec --- /dev/null +++ b/src/moldrug/meeko_light/molsetup.py @@ -0,0 +1,2328 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Meeko +# + +from abc import ABC, abstractmethod +from copy import deepcopy +from collections import defaultdict +from dataclasses import dataclass, field +import json +import sys +import warnings +from typing import List, Tuple, Union + +import numpy as np +from rdkit import Chem +from rdkit.Chem import rdPartialCharges +from rdkit.Chem import rdMolInterchange + +from .utils.jsonutils import rdkit_mol_from_json, tuple_to_string, string_to_tuple +from .utils import rdkitutils +from .utils import utils +from .utils.pdbutils import PDBAtomInfo + + +# region DEFAULT VALUES +DEFAULT_PDBINFO = None +DEFAULT_CHARGE = 0.0 +DEFAULT_COORD = np.array([0.0, 0.0, 0.0], dtype="float") +DEFAULT_ATOMIC_NUM = None +DEFAULT_ATOM_TYPE = None +DEFAULT_IS_IGNORE = False +DEFAULT_GRAPH = [] + +DEFAULT_BOND_ROTATABLE = False + +DEFAULT_RING_CORNER_FLIP = False +DEFAULT_RING_GRAPH = [] +DEFAULT_RING_IS_AROMATIC = False +DEFAULT_RING_CLOSURE_BONDS_REMOVED = [] +DEFAULT_RING_CLOSURE_PSEUDOS_BY_ATOM = defaultdict +# endregion + + +# region Helper Data Organization Classes +class UniqAtomParams: + """ + A helper class used to keep parameters organized in a particular way that lets them be more usable. + + Attributes + ---------- + params: list[] + can be thought of as rows + param_names: list[] + can be thought of as columns + """ + + def __init__(self): + self.params = [] # aka rows + self.param_names = [] # aka column names + + @classmethod + def from_dict(cls, dictionary): + """ + Creates an UniqAtomParams object, populates it with information from the input dictionary, then returns + the new object. + + Parameters + ---------- + dictionary: dict() + A dictionary containing the keys "params" and "param_names", where the value for "params" is parseable as + rows and the value for "param_names" contains the corresponding column data. + + Returns + ------- + A populated UniqAtomParams object + """ + uap = UniqAtomParams() + uap.params = [row.copy() for row in dictionary["params"]] + uap.param_names = dictionary["param_names"].copy() + return uap + + def get_indices_from_atom_params(self, atom_params): + """ + Retrieves the indices of specific atom parameters in the UniqAtomParams object. + + Parameters + ---------- + atom_params: dict() + A dict with keys that correspond to the param names already in the UniqAtomParams object. The values are + lists that should all be the same size, and + + Returns + ------- + A list of indices corresponding to the order of parameters in the atom_params value lists that indicates the + index of that "row" of parameters in UniqAtomParams params. + """ + nr_items = set([len(values) for key, values in atom_params.items()]) + if len(nr_items) != 1: + raise RuntimeError( + f"all lists in atom_params must have same length, got {nr_items}" + ) + if set(atom_params) != set(self.param_names): + msg = "parameter names in atom_params differ from internal ones\n" + msg += f" - in atom_params: {set(atom_params)}" + msg += f" - internal: {set(self.param_names)}" + raise RuntimeError(msg) + nr_items = nr_items.pop() + param_idxs = [] + for i in range(nr_items): + row = [atom_params[key][i] for key in self.param_names] + param_index = None + for j, existing_row in enumerate(self.params): + if row == existing_row: + param_index = j + break + param_idxs.append(param_index) + return param_idxs + + def add_parameter(self, new_param_dict): + # remove None values to avoid a column with only Nones + new_param_dict = {k: v for k, v in new_param_dict.items() if v is not None} + incoming_keys = set(new_param_dict.keys()) + existing_keys = set(self.param_names) + new_keys = incoming_keys.difference(existing_keys) + for new_key in new_keys: + self.param_names.append(new_key) + for row in self.params: + row.append(None) # fill in empty "cell" in new "column" + + new_row = [] + for key in self.param_names: + value = new_param_dict.get(key, None) + new_row.append(value) + + if len(new_keys) == 0: # try to match with existing row + for index, row in enumerate(self.params): + if row == new_row: + return index + + # if we are here, we didn't match + new_row_index = len(self.params) + self.params.append(new_row) + return new_row_index + + def add_molsetup( + self, molsetup, atom_params=None, add_atomic_nr=False, add_atom_type=False + ): + if "charge" in molsetup.atom_params or "atom_type" in molsetup.atom_params: + msg = '"charge" and "atom_type" found in molsetup.atom_params' + msg += " but are hard-coded to store molsetup.charge and" + msg += " molsetup.atom_type in the internal data structure" + raise RuntimeError(msg) + if atom_params is None: + atom_params = molsetup.atom_params + param_idxs = [] + for atom in molsetup.atoms: + if atom.is_ignore: + param_idx = None + else: + p = {k: v[atom.index] for (k, v) in molsetup.atom_params.items()} + if add_atomic_nr: + if "atomic_nr" in p: + raise RuntimeError( + "trying to add atomic_nr but it's already in atom_params" + ) + p["atomic_nr"] = atom.atomic_num + if add_atom_type: + if "atom_type" in p: + raise RuntimeError( + "trying to add atom_type but it's already in atom_params" + ) + p["atom_type"] = atom.atom_type + param_idx = self.add_parameter(p) + param_idxs.append(param_idx) + return param_idxs + + +@dataclass +class Atom: + index: int + pdbinfo: Union[str, PDBAtomInfo] = DEFAULT_PDBINFO + charge: float = DEFAULT_CHARGE + coord: np.ndarray = field(default_factory=np.ndarray) + atomic_num: int = DEFAULT_ATOMIC_NUM + atom_type: str = DEFAULT_ATOM_TYPE + is_ignore: bool = DEFAULT_IS_IGNORE + graph: List[int] = field(default_factory=list) + interaction_vectors: List[np.array] = field(default_factory=list) + + is_dummy: bool = False + is_pseudo_atom: bool = False + + @staticmethod + def from_json(obj: dict): + """ + Takes an object and attempts to deserialize it into an Atom object. + + Parameters + ---------- + obj: Object + This can be any object, but it should be a dictionary constructed by deserializing the JSON representation + of an Atom object. + + Returns + ------- + If the input is a dictionary corresponding to an Atom, will return an Atom with data populated from the + dictionary. Otherwise, returns the input object. + """ + # if the input object is not a dict, we know that it will not be parsable and is unlikely to be usable or + # safe data, so we should ignore it. + if not isinstance(obj, dict): + return obj + + # Check that all the keys we expect are in the object dictionary as a safety measure + expected_json_keys = { + "index", + "pdbinfo", + "charge", + "coord", + "atomic_num", + "atom_type", + "is_ignore", + "graph", + "interaction_vectors", + "is_dummy", + "is_pseudo_atom", + } + if set(obj.keys()) != expected_json_keys: + return obj + + # Constructs an atom object from the provided keys. + index = obj["index"] + pdbinfo = PDBAtomInfo(*obj["pdbinfo"]) + charge = obj["charge"] + coord = np.asarray(obj["coord"]) + atomic_num = obj["atomic_num"] + atom_type = obj["atom_type"] + is_ignore = obj["is_ignore"] + graph = obj["graph"] + interaction_vectors = [np.asarray(i) for i in obj["interaction_vectors"]] + is_dummy = obj["is_dummy"] + is_pseudo_atom = obj["is_pseudo_atom"] + output_atom = Atom( + index, + pdbinfo, + charge, + coord, + atomic_num, + atom_type, + is_ignore, + graph, + interaction_vectors, + is_dummy, + is_pseudo_atom, + ) + return output_atom + + +@dataclass +class Bond: + canon_id: Tuple[int] # two ints + index1: int + index2: int + rotatable: bool = DEFAULT_BOND_ROTATABLE + + def __init__( + self, + index1: int, + index2: int, + rotatable: bool = DEFAULT_BOND_ROTATABLE, + ): + self.canon_id = self.get_bond_id(index1, index2) + self.index1 = index1 + self.index2 = index2 + self.rotatable = rotatable + return + + @staticmethod + def get_bond_id(idx1: int, idx2: int): + """ + Generates a consistent, "canonical", bond id from a pair of atom indices in the graph. + + Parameters + ---------- + idx1: int + atom index of one of the atoms in the bond + idx2: int + atom index of the other atom in the bond + + Returns + ------- + canon_id: tuple + a tuple of the two indices in their canonical order. + """ + idx_min = min(idx1, idx2) + idx_max = max(idx1, idx2) + return idx_min, idx_max + + @staticmethod + def from_json(obj: dict): + """ + Takes an object and attempts to deserialize it into a Bond object. + + Parameters + ---------- + obj: Object + This can be any object, but it should be a dictionary constructed by deserializing the JSON representation + of a Bond object. + + Returns + ------- + If the input is a dictionary corresponding to a Bond, will return a Bond with data populated from the + dictionary. Otherwise, returns the input object. + """ + # if the input object is not a dict, we know that it will not be parsable and is unlikely to be usable or + # safe data, so we should ignore it. + if not isinstance(obj, dict): + return obj + + # Check that all the keys we expect are in the object dictionary as a safety measure + expected_json_keys = {"canon_id", "index1", "index2", "rotatable"} + if set(obj.keys()) != expected_json_keys: + return obj + + # Constructs a bond object from the provided keys. + index1 = obj["index1"] + index2 = obj["index2"] + rotatable = obj["rotatable"] + output_bond = Bond(index1, index2, rotatable) + return output_bond + + +@dataclass +class Ring: + ring_id: tuple + corner_flip: bool = DEFAULT_RING_CORNER_FLIP + graph: dict = field(default_factory=list) + is_aromatic: bool = DEFAULT_RING_IS_AROMATIC + + @staticmethod + def from_json(obj: dict): + """ + Takes an object and attempts to deserialize it into a Ring object. + + Parameters + ---------- + obj: Object + This can be any object, but it should be a dictionary constructed by deserializing the JSON representation + of a Ring object. + + Returns + ------- + If the input is a dictionary corresponding to a Ring, will return a Ring with data populated from the + dictionary. Otherwise, returns the input object. + """ + # if the input object is not a dict, we know that it will not be parsable and is unlikely to be usable or + # safe data, so we should ignore it. + if type(obj) is not dict: + return obj + + # Check that all the keys we expect are in the object dictionary as a safety measure + expected_json_keys = {"ring_id", "corner_flip", "graph", "is_aromatic"} + if set(obj.keys()) != expected_json_keys: + return obj + + # Constructs a Ring object from the provided keys. + ring_id = string_to_tuple(obj["ring_id"], int) + corner_flip = obj["corner_flip"] + graph = obj["graph"] + is_aromatic = obj["is_aromatic"] + output_ring = Ring(ring_id, corner_flip, graph, is_aromatic) + return output_ring + + +@dataclass +class RingClosureInfo: + bonds_removed: list = field(default_factory=list) + pseudos_by_atom: dict = DEFAULT_RING_CLOSURE_PSEUDOS_BY_ATOM + + +@dataclass +class Restraint: + atom_index: int + target_coords: Tuple[float] # three floats + kcal_per_angstrom_square: float + delay_angstroms: float + + def copy(self): + new_target_coords = ( + self.target_coords[0], + self.target_coords[1], + self.target_coords[2], + ) + new_restraint = Restraint( + self.atom_index, + new_target_coords, + self.kcal_per_angstrom_square, + self.delay_angstroms, + ) + return new_restraint + + @staticmethod + def from_json(obj: dict): + """ + Takes an object and attempts to deserialize it into a Restraint object. + + Parameters + ---------- + obj: Object + This can be any object, but it should be a dictionary constructed by deserializing the JSON representation + of a Restraint object. + + Returns + ------- + If the input is a dictionary corresponding to a Restraint, will return a Restraint with data populated from the + dictionary. Otherwise, returns the input object. + """ + # if the input object is not a dict, we know that it will not be parsable and is unlikely to be usable or + # safe data, so we should ignore it. + if type(obj) is not dict: + return obj + + # Check that all the keys we expect are in the object dictionary as a safety measure + expected_json_keys = { + "atom_index", + "target_coords", + "kcal_per_angstrom_square", + "delay_angstroms", + } + if set(obj.keys()) != expected_json_keys: + return obj + + # Constructs a Restraint object from the provided keys. + atom_index = obj["atom_index"] + target_coords = tuple(obj["target_coords"]) + kcal_per_angstrom_square = obj["kcal_per_angstrom_square"] + delay_angstroms = obj["delay_angstroms"] + output_restraint = Restraint( + atom_index, target_coords, kcal_per_angstrom_square, delay_angstroms + ) + return output_restraint + + +# endregion + + +class MoleculeSetup: + """ + Base MoleculeSetup Class, provides a way to store information about molecules for a number of purposes. + + Attributes + ---------- + name: str + is_sidechain: bool + pseudoatom_count: int + + atoms: list[Atom] + bond_info: dict[tuple, Bond] + rings: dict + ring_closure_info: RingClosureInfo + rotamers: list[dict] + + atom_params: dict + restraints: list[Restraint] + flexibility_model: dict + """ + + # region CLASS CONSTANTS + PSEUDOATOM_ATOMIC_NUM = 0 + # endregion + + def __init__(self, name: str = None, is_sidechain: bool = False): + # Molecule Setup Identity + self.name: str = name + self.is_sidechain: bool = is_sidechain + self.pseudoatom_count: int = 0 + + # Tracking atoms and bonds + self.atoms: list[Atom] = [] + self.bond_info: dict[tuple, Bond] = {} + self.rings: dict[tuple, Ring] = {} + self.ring_closure_info = RingClosureInfo([], {}) + self.rotamers: list[dict] = [] # TODO: revisit rotamer implementation + + self.atom_params: dict = {} + self.restraints: list = ( + [] + ) # TODO: determine whether restraints are being used anymore + + # TODO: redesign flexibility model to resolve some of the circular imports and to make it more structured + self.flexibility_model = None # from flexibility_model - from flexibility.py + + # region Manually Building A MoleculeSetup + def add_atom( + self, + atom_index: int = None, + overwrite: bool = False, + pdbinfo: Union[str, PDBAtomInfo] = DEFAULT_PDBINFO, + charge: float = DEFAULT_CHARGE, + coord: np.ndarray = None, + atomic_num: int = DEFAULT_ATOMIC_NUM, + atom_type: str = DEFAULT_ATOM_TYPE, + is_ignore: bool = DEFAULT_IS_IGNORE, + graph: List[int] = None, + ): + """ + Adds an atom with all the specified attributes to the MoleculeSetup, either at the specified atom index, or by + appending it to the internal list of atoms. Default values will be used for any attributes with unspecified + values. + + Parameters + ---------- + atom_index: int + atom index in the MoleculeSetup + overwrite: bool + can we overwrite other atoms may be in the same atom index as this one + pdbinfo: str + pdb string for the atom + coord: np.ndarray + the atom's coordinates + charge: float + partial charge to be loaded for the atom + atomic_num: int + the atomic number of the atom + atom_type: str + TODO: needs info + is_ignore: bool + ignore flag for the atom + graph: List[List[int]] + + Returns + ------- + None + + Raises + ------ + RuntimeException + If the user tries to overwrite an existing atom without explicitly allowing overwrites. + """ + # If atom index is specified and it would be trying to overwrite an existing atom in the atom list, raises a + # Runtime Exception + insert_disallowed = len(self.atoms) > atom_index and not overwrite + if ( + atom_index is not None + and insert_disallowed + and not self.atoms[atom_index].is_dummy + ): + raise RuntimeError( + "ADD_ATOM Error: the atom_index [%d] is already occupied (use 'overwrite' to force)" + ) + + # If atom index is not specified, appends the new atom to the end of the current atom list + if atom_index is None: + atom_index = len(self.atoms) + + # Inserts dummy atoms if a specified atom index is greater than the current length of the atom list + while atom_index > len(self.atoms): + self.atoms.append(Atom(len(self.atoms), is_dummy=True)) + + # Creates and adds new atom to the atom list + if coord is None: + coord = deepcopy(DEFAULT_COORD) + if graph is None: + graph = deepcopy(DEFAULT_GRAPH) + new_atom = Atom( + atom_index, + pdbinfo, + charge, + coord, + atomic_num, + atom_type, + is_ignore, + graph, + ) + if atom_index < len(self.atoms): + self.atoms[atom_index] = new_atom + return + self.atoms.append(new_atom) + return + + def add_pseudoatom( + self, + pdbinfo: Union[str, PDBAtomInfo] = DEFAULT_PDBINFO, + charge: float = DEFAULT_CHARGE, + coord: np.ndarray = None, + atom_type: str = DEFAULT_ATOM_TYPE, + is_ignore: bool = DEFAULT_IS_IGNORE, + anchor_list: List[int] = None, + rotatable: bool = False, + directional_vectors: List[int] = None, + ): + """ + Adds a pseudoatom with all the specified attributes to the MoleculeSetup. Default values will be used for any + attributes with unspecified values. Multiple bonds can be specified to support the centroids of aromatic rings. + If rotatable, makes the anchor atom rotatable to allow the pseudoatom movement. + + Parameters + ---------- + pdbinfo: str + PDB string for the pseudoatom. + charge: float + partial charge for the pseudoatom + coord: np.ndarray + the pseudoatom's coordinates + atom_type: str + TODO: needs info + is_ignore: bool + ignore flag for the pseudoatom + anchor_list: list[int] + a list of ints indicating the multiple bonds that can be specified as input + rotatable: bool + flag indicating if the anchor atom should be marked as rotatable to allow the pseudoatom movement. + directional_vectors + TODO: needs info + + Returns + ------- + pseudoatom_index: int + The atom_index of the added pseudoatom + + Raises + ------ + RuntimeError: + When the incorrect number of anchors of pseudoatoms are found in rigid_groups in the flexibility model + + """ + # Places the atom at the end of the atom list. + pseudoatom_index = len(self.atoms) + # Creates the atom and marks it as a pseudoatom + if coord is None: + coord = deepcopy(DEFAULT_COORD) + new_pseudoatom = Atom( + pseudoatom_index, + pdbinfo=pdbinfo, + charge=charge, + coord=coord, + atomic_num=self.PSEUDOATOM_ATOMIC_NUM, + atom_type=atom_type, + is_ignore=is_ignore, + is_pseudo_atom=True, + ) + self.atoms.append(new_pseudoatom) + # Adds bonds for all of the provided anchor atoms + if anchor_list is not None: + for anchor in anchor_list: + self.add_bond(pseudoatom_index, anchor, rotatable=rotatable) + # Adds directional vectors [Check what this is used for/if this is used] + if directional_vectors is not None: + self._add_interaction_vectors(pseudoatom_index, directional_vectors) + # If there are no specified anchor atoms, + if not self.flexibility_model or not anchor_list: + return pseudoatom_index + # TODO: revise this logic + # If there is a flexibility model in the MoleculeSetup, adds the psuedoatom to the flexibility model's rigid + # group tracking + rigid_groups_indices = [] + for anchor in anchor_list: + for rigid_index, members in self.flexibility_model[ + "rigid_body_members" + ].items(): + if anchor in members: + rigid_groups_indices.append(rigid_index) + if len(rigid_groups_indices) != 1: + raise RuntimeError( + f"anchors of pseudo atom found in {len(rigid_groups_indices)} rigid_groups (must be 1)" + ) + rigid_index = rigid_groups_indices[0] + self.flexibility_model["rigid_body_members"][rigid_index].append( + pseudoatom_index + ) + # returns the pseudoatom index + return pseudoatom_index + + def delete_atom(self, atom_index: int): + """ + Clears the atom data at a specified atom index and replaces the atom with a dummy atom. + + Parameters + ---------- + atom_index: int + atom index to replace with a dummy atom + + Returns + ------- + None + """ + blank_atom = Atom(atom_index, is_dummy=True) + self.atoms[atom_index] = blank_atom + return + + def add_bond( + self, + atom_index_1: int, + atom_index_2: int, + rotatable: bool = DEFAULT_BOND_ROTATABLE, + ) -> None: + """ + Creates a bond and adds it to all the internal data structures where atom bonds are being tracked. + + Parameters + ---------- + atom_index_1: int + Atom index of one of the atoms in the bond + atom_index_2: int + Atom index of the other atom in the bond + rotatable: bool + Indicates whether the bond is rotatable + + Returns + ------- + None + + Raises + ------ + IndexError: + When one or more the given bond atom indices do not exist in the MoleculeSetup + """ + # Checks that both of the atom indices provided are valid indices, otherwise throws an error + if len(self.atoms) <= atom_index_1 or len(self.atoms) <= atom_index_2: + raise IndexError( + "ADD_BOND: provided atom indices outside the range of atoms currently in MoleculeSetup" + ) + # Adds each atom to the other's bond graph + if atom_index_2 not in self.atoms[atom_index_1].graph: + self.atoms[atom_index_1].graph.append(atom_index_2) + if atom_index_1 not in self.atoms[atom_index_2].graph: + self.atoms[atom_index_2].graph.append(atom_index_1) + # Creates new bond object and uses its internal canonical bond id to add it to MoleculeSetup bond tracking. + new_bond = Bond(atom_index_1, atom_index_2, rotatable) + self.bond_info[new_bond.canon_id] = new_bond + return + + def delete_bond(self, atom_index_1: int, atom_index_2: int): + """ + Deletes a bond from the molecule setup. + + Parameters + ---------- + atom_index_1: int + The atom index of one of the atoms in the bond to delete + atom_index_2: int + The atom index of the other atom in the bond to delete + + Returns + ------- + None + """ + # Gets canon bond id for the bond to delete + canon_bond_id = Bond.get_bond_id(atom_index_1, atom_index_2) + # Deletes the bond from the internal bond table + del self.bond_info[canon_bond_id] + # Removes the bond from each atom's graph + self.atoms[atom_index_1].graph.remove(atom_index_2) + self.atoms[atom_index_2].graph.remove(atom_index_1) + return + + def add_rotamers( + self, index_list: List[Tuple[int]], angle_list: np.ndarray + ): + """ + Adds rotamers to the internal record of rotamers. + + Parameters + ---------- + index_list: list[(int, int, int, int)] + angle_list: np.ndarray + + Returns + ------- + None + """ + # It's unclear how this will work without the coordinates in the moleculesetup. food for thought. + # TODO: address issues with the lack of coords and add detail in function comment + rotamers = {} + for (_, idx2, idx3, _), angle in zip(index_list, angle_list): + bond_id = Bond.get_bond_id(idx2, idx3) + if bond_id in rotamers: + raise RuntimeError(f"repeated bond {bond_id}") + if not self.bond_info[bond_id].rotatable: + raise RuntimeError( + f"trying to add rotamer for non rotatable bond {bond_id}") + dihedral = 0 # TODO: fix this + rotamers[bond_id] = angle - dihedral + self.rotamers.append(rotamers) + return + + def delete_rotamers( + self, + bond_id_list: List[Tuple] = None, + index_list: List[Tuple[int]] = None, + ): + """ + Deletes rotamers from the internal list of rotamers, either by using bond ids or by generating bond ids from a + list of input indices. + + Parameters + ---------- + bond_id_list: list[tuple] + index_list: list[(int, int, int, int)] + + Returns + ------- + None + """ + # loops through the index list, generates bond ids from the provided indices and adds them to the bond_id_list + if index_list is not None: + for _, idx2, idx3, _ in index_list: + bond_id = Bond.get_bond_id(idx2, idx3) + bond_id_list.append(bond_id) + # deletes all bond_ids in bond_id_list from self.rotamers + if bond_id_list is not None: + for bond_id in bond_id_list: + if bond_id in self.rotamers: + del self.rotamers[bond_id] + return + + def _add_interaction_vectors(self, atom_index: int, vector_list: List[np.array]): + """ + Adds input vector list to the list of directional interaction vectors for the specified atom. + + Parameters + ---------- + atom_index: int + index of the atom to add the vectors to + vector_list: list[np.array] + a list of directional interaction vectors + + Returns + ------- + None + + Raises + ------ + IndexError + if the specified atom index does not exist or is a dummy atom. + """ + if atom_index > len(self.atoms) or self.atoms[atom_index].is_dummy: + raise IndexError( + "INTERACTION_VECTORS: provided atom index is out of range or is a dummy atom." + ) + for vector in vector_list: + self.atoms[atom_index].interaction_vectors.append(vector) + return + + @property + def true_atom_count(self): + """ + Counts the number of atoms in the MoleculeSetup that are not pseudo_atoms or marked as dummy atoms + + Returns + ------- + count: int + The number of atoms currently in the MoleculeSetup that are not dummy atoms or pseudo_atoms. + """ + count = 0 + for atom in self.atoms: + if not atom.is_pseudo_atom and not atom.is_dummy: + count += 1 + return count + + # this might void the graph connections and everything in bonds, might need to add a dict of the changes we're + # making and then use that save this for a future push. + def clean_atoms(self, remove_pseudoatoms: bool = False): + """ + Cleans dummy and potentially also pseudoatoms from the MoleculeSetup so only true atoms remain. Note that this + is pretty slow and should not be done often. + + Parameters + ---------- + remove_pseudoatoms: bool + Indicates if we want to remove all the pseudoatoms from the MoleculeSetup. + + Returns + ------- + The number of atoms removed from the MoleculeSetup. + """ + new_atoms = [] + removed_atom_count = 0 + # atom_index_mapping = {} + for atom in self.atoms: + if remove_pseudoatoms and atom.is_pseudo_atom: + removed_atom_count += 1 + continue + if atom.is_dummy: + removed_atom_count += 1 + continue + atom.index = atom.index - removed_atom_count + new_atoms.append(atom) + self.atoms = new_atoms + if remove_pseudoatoms: + self.pseudoatom_count = 0 + return removed_atom_count + + # endregion + + # region Getters and Setters + + def get_pdbinfo(self, atom_index: int): + """ + Retrieves the PDB Info string for the atom with the specified atom index. + + Parameters + ---------- + atom_index: int + Atom index to retrieve data for. + + Returns + ------- + pdbinfo: str + A string containing the pdb information for the atom + + Raises + ------ + IndexError: + When the provided atom index does not exist in the MoleculeSetup or the atom index does not contain + data. + """ + if atom_index > len(self.atoms) or self.atoms[atom_index].is_dummy: + raise IndexError( + "GET_PDBINFO: provided atom index is out of range or is a dummy atom" + ) + return self.atoms[atom_index].pdbinfo + + def get_charge(self, atom_index: int): + """ + Retrieves the partial charge for the atom with the specified atom index. + + Parameters + ---------- + atom_index: int + Atom index to retrieve data for. + + Returns + ------- + charge: float + The charge associated with the atom + + Raises + ------ + IndexError: + When the provided atom index does not exist in the MoleculeSetup or the atom index does not contain + data. + """ + if atom_index > len(self.atoms) or self.atoms[atom_index].is_dummy: + raise IndexError( + "GET_CHARGE: provided atom index is out of range or is a dummy atom" + ) + return self.atoms[atom_index].charge + + def get_coord(self, atom_index: int): + """ + Retrieves the coordinates for the atom with the specified atom index. + + Parameters + ---------- + atom_index: int + Atom index to retrieve data for. + + Returns + ------- + coord: np.ndarray + The coordinates associated with the atom. + + Raises + ------ + IndexError: + When the provided atom index does not exist in the MoleculeSetup or the atom index does not contain + data. + """ + if atom_index > len(self.atoms) or self.atoms[atom_index].is_dummy: + raise IndexError( + "GET_CHARGE: provided atom index is out of range or is a dummy atom" + ) + return self.atoms[atom_index].coord + + def get_atomic_num(self, atom_index: int): + """ + Retrieves the atomic number for the atom with the specified atom index. + + Parameters + ---------- + atom_index: int + Atom index to retrieve data for. + + Returns + ------- + atomic_num: int + The atomic number associated with an atom. + + Raises + ------ + IndexError: + When the provided atom index does not exist in the MoleculeSetup or the atom index does not contain + data. + """ + if atom_index > len(self.atoms) or self.atoms[atom_index].is_dummy: + raise IndexError( + "GET_ATOMIC_NUM: provided atom index is out of range or is a dummy atom" + ) + return self.atoms[atom_index].atomic_num + + def get_atom_type(self, atom_index: int): + """ + Retrieves the atom type for the atom with the specified atom index. + + Parameters + ---------- + atom_index: int + Atom index to retrieve data for. + + Returns + ------- + charge: str + The atom index associated with the atom + + Raises + ------ + IndexError: + When the provided atom index does not exist in the MoleculeSetup or the atom index does not contain data. + """ + if atom_index > len(self.atoms) or self.atoms[atom_index].is_dummy: + raise IndexError( + "GET_ATOM_TYPE: provided atom index is out of range or is a dummy atom" + ) + return self.atoms[atom_index].atom_type + + def set_atom_type(self, atom_index: int, atom_type: str) -> None: + """ + Sets the atom type for the atom with the specified atom index. + + Parameters + ---------- + atom_index: int + Atom index to set atom_type for. + atom_type + Atom type string to set. + + Returns + ------- + None + + Raises + ------ + IndexError: + When the provided atom index does not exist in the MoleculeSetup or the atom index does not contain data. + """ + if atom_index > len(self.atoms) or self.atoms[atom_index].is_dummy: + raise IndexError( + "SET_ATOM_TYPE: provided atom index is out of range or is a dummy atom" + ) + self.atoms[atom_index].atom_type = atom_type + return + + def set_atom_type_from_uniq_atom_params( + self, uniq_atom_params: UniqAtomParams, prefix: str + ): + """ + Uses a UniqAtomParams object to set the atom_type attribute for atoms in the Molecule Setup object. Adds the specified prefix + to each of the atom_type attributes pulled from UniqAtomParams. + + Parameters + ---------- + uniq_atom_params: UniqAtomParams + A uniq atom params object to extract atom_type from + prefix: string + A prefix to be appended to all the atom_type attributes + + Returns + ------- + None + """ + # Gets a mapping from parameter indices in atom_params to those in uniq_atom_params + parameter_indices = uniq_atom_params.get_indices_from_atom_params( + self.atom_params + ) + # Checks that we have the correct number of retrieved indices. + if len(parameter_indices) != len(self.atoms): + raise RuntimeError( + "Number of parameters ({len(parameter_indices)}) not equal to number of atoms in Molecule Setup ({len(self.atom_type)})" + ) + # Loops through the indices in parameter indices and sets atom types with the input prefix + for i, j in enumerate(parameter_indices): + self.atom_type[i] = f"{prefix}{j}" + return None + + def get_is_ignore(self, atom_index: int): + """ + Retrieves the is_ignore boolean for the atom with the specified atom index. + + Parameters + ---------- + atom_index: int + Atom index to retrieve data for. + + Returns + ------- + is_ignore: bool + Indicates whether a particular atom should be ignored + + Raises + ------ + IndexError: + When the provided atom index does not exist in the MoleculeSetup or the atom index does not contain + data. + """ + if atom_index > len(self.atoms) or self.atoms[atom_index].is_dummy: + raise IndexError( + "GET_IS_IGNORE: provided atom index is out of range or is a dummy atom" + ) + return self.atoms[atom_index].is_ignore + + def get_neighbors(self, atom_index: int): + """ + Retrieves the partial charge for the atom with the specified atom index. + + Parameters + ---------- + atom_index: int + Atom index to retrieve data for. + + Returns + ------- + graph: list[int] + The graph of the atoms connections to other atoms. + + Raises + ------ + IndexError: + When the provided atom index does not exist in the MoleculeSetup or the atom index does not contain + data. + """ + if atom_index > len(self.atoms) or self.atoms[atom_index].is_dummy: + raise IndexError( + "GET_GRAPH: provided atom index is out of range or is a dummy atom" + ) + return self.atoms[atom_index].graph + + def get_interaction_vectors(self, atom_index: int): + if atom_index > len(self.atoms) or self.atoms[atom_index].is_dummy: + raise IndexError( + "GET_INTERACTION_VECTORS: provided atom index is out of range or is a dummy atom" + ) + return self.atoms[atom_index].interaction_vectors + + # endregion + + def merge_terminal_atoms(self, indices) -> None: + """ + Primarily for merging hydrogens, but will merge the data for any atom or pseudoatom that is bonded to only one + other atom. + + Parameters + ---------- + indices: list + A list of indices to merge + + Returns + ------- + None + """ + for index in indices: + if len(self.get_neighbors(index)) != 1: + msg = "Atempted to merge atom %d with %d neighbors. " + msg += "Only atoms with one neighbor can be merged." + msg = msg % (index + 1, self.get_neighbors(index)) + raise RuntimeError(msg) + neighbor_index = self.get_neighbors(index)[0] + self.atoms[neighbor_index].charge += self.get_charge(index) + self.atoms[index].charge = 0.0 + self.atoms[index].is_ignore = True + return + + # NOTE: This is a candidate for moving to utils + @staticmethod + def get_bonds_in_ring(ring: Tuple) -> List[Tuple]: + """ + Takes as input a tuple of atom indices corresponding to atoms in a ring and returns a list of all the bonds ids + in the ring. + + Parameters + ---------- + ring: tuple + A list of atom indices of the atoms in a ring. + + Returns + ------- + A list of canonical bond id tuples for the bonds in the ring. + """ + bonds = [] + num_indices = len(ring) + for i in range(num_indices): + bond = (ring[i], ring[(i + 1) % num_indices]) + bond = Bond.get_bond_id(bond[0], bond[1]) + bonds.append(bond) + return bonds + + def _recursive_graph_walk( + self, idx: int, collected: List[int] = None, exclude: List[int] = None + ): + """ + Recursively walks through a molecular graph and returns bond-connected subgroups. + + Parameters + ---------- + idx: int + atom index to start the recursive walk from + collected: list[int] + a list of connected subgroups + exclude: list[int] + a list of atom indices to exclude from the final walk. + + Returns + ------- + A list of ints indicating the subgroups that are bond-connected. + """ + if collected is None: + collected = [] + if exclude is None: + exclude = [] + for neighbor in self.get_neighbors(idx): + if neighbor in collected or neighbor in exclude: + continue + collected.append(neighbor) + self._recursive_graph_walk(neighbor, collected, exclude) + return collected + + def write_coord_string(self) -> str: + """ + Constructs and returns a string of all atom and pseudoatom elements and coordinates. + + Returns + ------- + A string of all atom and pseudoatom elements and coordinates. + """ + n = len(self.atoms) + output_string = "%d\n\n" % n + for index in range(n): + element = "Ne" + if self.atoms[index].is_dummy: + continue + if not self.atoms[index].is_pseudo_atom: + element = utils.mini_periodic_table[self.atoms[index].atomic_num] + x, y, z = self.atoms[index].coord + output_string += "%3s %12.6f %12.6f %12.6f\n" % (element, x, y, z) + return output_string + + def show(self) -> None: + """ + Legacy function to print the contents of a MoleculeSetup in a human-readable format. + + Returns + ------- + None + """ + total_charge = 0 + + print("Molecule Setup\n") + print( + "==============[ ATOMS ]===================================================" + ) + print("idx | coords | charge |ign| atype | connections") + print( + "-----+----------------------------+--------+---+----------+--------------- . . . " + ) + for atom in self.atoms: + print( + "% 4d | % 8.3f % 8.3f % 8.3f | % 1.3f | %d" + % ( + atom.index, + atom.coord[0], + atom.coord[1], + atom.coord[2], + atom.charge, + atom.is_ignore, + ), + "| % -8s |" % atom.atom_type, + atom.graph, + ) + total_charge += atom.charge + print( + "-----+----------------------------+--------+---+----------+--------------- . . . " + ) + print(" TOT CHARGE: %3.3f" % total_charge) + + print("\n==============[ BONDS ]================") + # For sanity users, we won't show those keys for now + keys_to_not_show = ["type"] + for bond_id, bond in list(self.bond_info.items()): + t = ", ".join( + "%s: %s" % (i, j) + for i, j in bond.__dict__.items() + if i not in keys_to_not_show + ) + print("% 8s - " % str(bond_id), t) + print("") + return + + def to_json(self): + """ + Converts MoleculeSetup object to a JSON string. + + Returns + ------- + A JSON string representation of the MoleculeSetup object. + """ + return json.dumps(self, cls=MoleculeSetupEncoder) + + @staticmethod + def from_json(obj): + """ + Takes an object and attempts to decode it into a MoleculeSetup object. + + Parameters + ---------- + obj: Object + This can be any object, but it should be a dictionary generated by deserializing a JSON of a MoleculeSetup + object. + + Returns + ------- + If the input is a dictionary corresponding to a MoleculeSetup, will return a MoleculeSetup with data + populated from the dictionary. Otherwise, returns the input object. + """ + # if the input object is not a dict, we know that it will not be parsable and is unlikely to be usable or + # safe data, so we should ignore it. + if type(obj) is not dict: + return obj + + # checks that all the keys we expect are in the object dictionary as a safety measure. Allows other keys + # to exist in case this is being called by a subclass of MoleculeSetup that may have more keys in its dict. + expected_molsetup_keys = { + "name", + "is_sidechain", + "pseudoatom_count", + "atoms", + "bond_info", + "rings", + "ring_closure_info", + "rotamers", + "atom_params", + "restraints", + "flexibility_model", + } + for key in expected_molsetup_keys: + if key not in obj.keys(): + return obj + + # Constructs a MoleculeSetup object and restores the expected attributes + name = obj["name"] + is_sidechain = obj["is_sidechain"] + molsetup = MoleculeSetup(name, is_sidechain) + molsetup.pseudoatom_count = obj["pseudoatom_count"] + molsetup.atoms = [Atom.from_json(x) for x in obj["atoms"]] + molsetup.bond_info = { + string_to_tuple(k, int): Bond.from_json(v) + for k, v in obj["bond_info"].items() + } + molsetup.rings = { + string_to_tuple(k, int): Ring.from_json(v) for k, v in obj["rings"].items() + } + molsetup.ring_closure_info = RingClosureInfo( + obj["ring_closure_info"]["bonds_removed"], + obj["ring_closure_info"]["pseudos_by_atom"], + ) + molsetup.rotamers = obj["rotamers"] + molsetup.atom_params = obj["atom_params"] + molsetup.restraints = [Restraint.from_json(x) for x in obj["restraints"]] + molsetup.flexibility_model = obj["flexibility_model"] + if "rigid_body_connectivity" in molsetup.flexibility_model: + tuples_rigid_body_connectivity = { + string_to_tuple(k, int): string_to_tuple(v) + for k, v in molsetup.flexibility_model[ + "rigid_body_connectivity" + ].items() + } + molsetup.flexibility_model["rigid_body_connectivity"] = ( + tuples_rigid_body_connectivity + ) + if "rigid_body_graph" in molsetup.flexibility_model: + molsetup.flexibility_model["rigid_body_graph"] = { + int(k): v + for k, v in molsetup.flexibility_model["rigid_body_graph"].items() + } + if "rigid_body_members" in molsetup.flexibility_model: + molsetup.flexibility_model["rigid_body_members"] = { + int(k): v + for k, v in molsetup.flexibility_model["rigid_body_members"].items() + } + if "rigid_index_by_atom" in molsetup.flexibility_model: + molsetup.flexibility_model["rigid_index_by_atom"] = { + int(k): v + for k, v in molsetup.flexibility_model["rigid_index_by_atom"].items() + } + return molsetup + + +# region External Toolkit Support +class MoleculeSetupExternalToolkit(ABC): + """ + Additional functions and requirements to extend the MoleculeSetup class in order to use it with external toolkits + such as RDKit and OpenBabel. + + Required Attributes + ------------------- + dihedral_interactions: list + A list of fourier series [add detail] + """ + + @staticmethod + def are_fourier_series_identical(series1: list, series2: list) -> bool: + """ + Compares two fourier series represented as lists of dictionaries. + + Parameters + ---------- + series1: list[dict] + The first fourier series to compare. + series2: list[dict] + The second fourier series to compare. + + Returns + ------- + A bool indicicating whether the fourier series are equal. + """ + # Gets the indices of both series by periodicity and checks for equality + index_by_periodicity1 = { + series1[index]["periodicity"]: index for index in range(len(series1)) + } + index_by_periodicity2 = { + series2[index]["periodicity"]: index for index in range(len(series2)) + } + if index_by_periodicity1 != index_by_periodicity2: + return False + # After establishing equality of the indices, loops through periodicity abd checks that the values stored in + # each fourier series dictionary are equal. + for periodicity in index_by_periodicity1: + index1 = index_by_periodicity1[periodicity] + index2 = index_by_periodicity2[periodicity] + for key in ["k", "phase", "periodicity"]: + if series1[index1][key] != series2[index2][key]: + return False + return True + + def add_dihedral_interaction(self, fourier_series): + """ + Adds a safe copy of the input fourier series to the dihedral_interactions list if the fourier series is not + already in the list. + + Parameters + ---------- + fourier_series: list[dict] + + Returns + ------- + index: int + The index of the input fourier series in the dihedral interactions list. + """ + index = 0 + for existing_fs in self.dihedral_interactions: + if self.are_fourier_series_identical(existing_fs, fourier_series): + return index + index += 1 + safe_copy = json.loads(json.dumps(fourier_series)) + self.dihedral_interactions.append(safe_copy) + return index + + @abstractmethod + def init_atom(self, assign_charges, coords): + pass + + @abstractmethod + def init_bond(self): + pass + + @abstractmethod + def get_mol_name(self): + pass + + @abstractmethod + def find_pattern(self, smarts: str): + pass + + @abstractmethod + def get_smiles_and_order(self): + pass + + +class RDKitMoleculeSetup(MoleculeSetup, MoleculeSetupExternalToolkit): + """ + Subclass of MoleculeSetup, used to represent MoleculeSetup objects working with RDKit objects + + Attributes + ---------- + mol : rdkit.Chem.rdchem.Mol + An RDKit Mol object to base the Molecule Setup on. + modified_atom_positions: list + List of dictionaries where keys are atom indices, Used to store sets of coordinates, e.g. docked poses, as + dictionaries indexed by the atom index, because not all atoms need to have new coordinates specified. + Unspecified hydrogen positions bonded to modified heavy atom positions are to be calculated "on-the-fly". + dihedral_interactions: list[] + A list of unique fourier_series, each of which are represented as a list of dictionaries. + dihedral_partaking_atoms: dict() + a mapping from tuples of atom indices to the indices in dihedral_interactions + dihedral_labels: dict() + a mapping from tuples of atom indices to dihedral labels + atom_to_ring_id: dict() + mapping of atom index to ring id of each atom belonging to the ring + ring_corners: dict() + unclear what this is a mapping of, but is used to store corner flexibility for the rings + rmsd_symmetry_indices: tuple + Tuples of the indices of the molecule's atoms that match a substructure query. needs info. + + Methods + ------- + from_mol() + constructor for the RDKitMoleculeSetup object (consider adapting to init?) + """ + + def __init__(self, name: str = None, is_sidechain: bool = False): + super().__init__(name, is_sidechain) + self.mol = None + self.modified_atom_positions = [] + self.dihedral_interactions: list[dict] = [] + self.dihedral_partaking_atoms: dict = {} + self.dihedral_labels: dict = {} + self.atom_to_ring_id = {} + self.ring_corners = {} + self.rmsd_symmetry_indices = () + + def copy(self): + """ + Returns a copy of the current RDKitMoleculeSetup. + """ + newsetup = RDKitMoleculeSetup() + for key, value in self.__dict__.items(): + if key != "mol": + newsetup.__dict__[key] = deepcopy(value) + newsetup.mol = Chem.Mol(self.mol) # not sure how deep of a copy this is + return newsetup + + @classmethod + def from_mol( + cls, + mol: Chem.Mol, + keep_chorded_rings: bool = False, + keep_equivalent_rings: bool = False, + assign_charges: bool = True, + conformer_id: int = -1, + ): + """ + + Parameters + ---------- + mol: rdkit.Chem.rdchem.Mol + RDKit Mol object to build the RDKitMoleculeSetup from. + keep_chorded_rings: bool + keep_equivalent_rings: bool + assign_charges: bool + conformer_id: int + + Returns + ------- + molsetup: RDKitMoleculeSetup + A populated RDKitMoleculeSetup object + + Raises + ------ + ValueError: + If the RDKit Mol has implicit Hydrogens or if there are no conformers for the given RDKit Mol + """ + # Checks if the input molecule is valid + if cls.has_implicit_hydrogens(mol): + raise ValueError("RDKit molecule has implicit Hs. Need explicit Hs.") + if mol.GetNumConformers() == 0: + raise ValueError( + "RDKit molecule does not have a conformer. Need 3D coordinates." + ) + + # Gets the RDKit Conformer that we are going to load into the molecule setup + rdkit_conformer = mol.GetConformer(conformer_id) + if not rdkit_conformer.Is3D(): + warnings.warn( + "RDKit molecule not labeled as 3D. This warning won't show again." + ) + RDKitMoleculeSetup.warned_not3D = True + if mol.GetNumConformers() > 1 and conformer_id == -1: + msg = "RDKit molecule has multiple conformers. Considering only the first one." + print(msg, file=sys.stderr) + if len(Chem.GetMolFrags(mol)) != 1: + raise ValueError(f"RDKit molecule has {len(Chem.GetMolFrags(mol))} fragments. Must have 1.") + if mol.HasQuery(): + raise ValueError("RDKit molecule has query. Check exotic fields (atom or bond) in SDF.") + + # Creating and populating the molecule setup with properties from RDKit as well as calculated values from our + # functions + molsetup = cls() + molsetup.mol = mol + molsetup.atom_true_count = molsetup.get_num_mol_atoms() + molsetup.name = molsetup.get_mol_name() + coords = rdkit_conformer.GetPositions() + molsetup.init_atom(assign_charges, coords) + molsetup.init_bond() + molsetup.perceive_rings(keep_chorded_rings, keep_equivalent_rings) + molsetup.rmsd_symmetry_indices = cls.get_symmetries_for_rmsd(mol) + + # to store sets of coordinates, e.g. docked poses, as dictionaries indexed by + # the atom index, because not all atoms need to have new coordinates specified + # Unspecified hydrogen positions bonded to modified heavy atom positions + # are to be calculated "on-the-fly". + molsetup.modified_atom_positions = ( + [] + ) # list of dictionaries where keys are atom indices + + return molsetup + + @staticmethod + def remove_elements(mol, to_rm=(12, 20, 25, 26, 30)): + idx_to_rm = {} + neigh_idx_to_nr_h = {} + rm_to_neigh = {} + for atom in mol.GetAtoms(): + if atom.GetAtomicNum() in to_rm: + idx_to_rm[atom.GetIdx()] = atom.GetFormalCharge() + rm_to_neigh[atom.GetIdx()] = set() + for neigh in atom.GetNeighbors(): + n = neigh.GetNumExplicitHs() + neigh_idx_to_nr_h[neigh.GetIdx()] = n + rm_to_neigh[atom.GetIdx()].add(neigh.GetIdx()) + if not idx_to_rm: + return Chem.Mol(mol), idx_to_rm, rm_to_neigh + rwmol = Chem.EditableMol(mol) + for idx in sorted(idx_to_rm, reverse=True): + rwmol.RemoveAtom(idx) + mol = rwmol.GetMol() + for idx in neigh_idx_to_nr_h: + n = neigh_idx_to_nr_h[idx] + newidx = idx - sum([i < idx for i in idx_to_rm]) + mol.GetAtomWithIdx(newidx).SetNumExplicitHs(n + 1) + mol.UpdatePropertyCache() + Chem.SanitizeMol(mol) + mol = Chem.AddHs(mol) + return mol, idx_to_rm, rm_to_neigh + + def init_atom(self, assign_charges: bool, coords: List[np.ndarray]): + """ + Generates information about the atoms in an RDKit Mol and adds them to an RDKitMoleculeSetup. + + Parameters + ---------- + assign_charges: bool + Indicates whether we should extract/generate charges. + coords: list[np.ndarray] + Atom coordinates for the RDKit Mol. + + Returns + ------- + None + """ + # extract/generate charges + if assign_charges: + things = self.remove_elements(self.mol) + copy_mol, idx_rm_to_formal_charge, rm_to_neigh = things + for atom in copy_mol.GetAtoms(): + if atom.GetAtomicNum() == 34: + atom.SetAtomicNum(16) + rdPartialCharges.ComputeGasteigerCharges(copy_mol) + charges = [a.GetDoubleProp("_GasteigerCharge") for a in copy_mol.GetAtoms()] + if idx_rm_to_formal_charge: + ok_charges = charges.copy() + for i in sorted(idx_rm_to_formal_charge, reverse=True): + ok_charges.insert(i, 0.0) + nr_rm = len(idx_rm_to_formal_charge) + nr_added_h = copy_mol.GetNumAtoms() - self.mol.GetNumAtoms() + nr_rm + ok_charges = ok_charges[:len(ok_charges)-nr_added_h] + # print(f"{nr_added_h=}") + # print(f"{nr_rm=}") + # print(f"{idx_rm_to_formal_charge=}") + # print(f"{len(charges)=}") + # print(f"{len(ok_charges)=}") + # print(f"{copy_mol.GetNumAtoms()=}") + # print(f"{self.mol.GetNumAtoms()=}") + chrg_by_heavy_atom = {} + for i in range(nr_added_h): + added_H_idx = self.mol.GetNumAtoms() + i - nr_rm + # print(f"{added_H_idx=}") + neighs = copy_mol.GetAtomWithIdx(added_H_idx).GetNeighbors() + if len(neighs) != 1: + raise RuntimeError("H should have 1 neighbor") + if neighs[0].GetIdx() in chrg_by_heavy_atom: + raise RuntimeError("expected only 1 added H per heavy atom, maybe deleted element had double bond to this heavy atom") + chrg_by_heavy_atom[neighs[0].GetIdx()] = charges[added_H_idx] + # print(f"{chrg_by_heavy_atom=}") + for i, neighs in rm_to_neigh.items(): + # print(f"{i=}, {neighs=}") + ok_charges[i] += idx_rm_to_formal_charge[i] + for idx in neighs: + newidx = idx - sum([i <= idx for i in idx_rm_to_formal_charge]) + # print(f"{idx=} {newidx=}") + ok_charges[i] += chrg_by_heavy_atom[newidx] + charges = ok_charges + else: + charges = [0.0] * self.mol.GetNumAtoms() + # register atom + for a in self.mol.GetAtoms(): + idx = a.GetIdx() + self.add_atom( + atom_index=idx, + pdbinfo=rdkitutils.getPdbInfoNoNull(a), + charge=charges[idx], + coord=coords[idx], + atomic_num=a.GetAtomicNum(), + is_ignore=False, + ) + + def init_bond(self): + """ + Uses the RDKit mol to initialize bond info for the RDKitMoleculeSetup + + Returns + ------- + None + """ + for b in self.mol.GetBonds(): + idx1 = b.GetBeginAtomIdx() + idx2 = b.GetEndAtomIdx() + rotatable = int(b.GetBondType()) == 1 + self.add_bond(idx1, idx2, rotatable=rotatable) + + def find_pattern(self, smarts: str): + """ + Given a SMARTS pattern, finds substruct matches in the molecule. + + Parameters + ---------- + smarts: + A SMARTS string to find in the RDKit Mol object + + Returns + ------- + The substruct matches in the RDKit Mol for the given SMARTS. + """ + p = Chem.MolFromSmarts(smarts) + return self.mol.GetSubstructMatches(p) + + def get_mol_name(self): + """ + Gets the RDKit Mol's name from self.mol. + + Returns + ------- + If the mol has a name, returns the name property. + """ + if self.mol.HasProp("_Name"): + return self.mol.GetProp("_Name") + else: + return None + + # TODO: Add more inline comments and clean up this function + def get_smiles_and_order(self): + """ + Returns the SMILES string and the mapping between atom indices in the SMILES and self.molof an atom after + running RDKit's RemoveHs function. + + Returns + ------- + smiles: + order: + """ + mol_no_ignore = self.mol + + # 3D SDF files written by other toolkits (OEChem, ChemAxon) + # seem to not include the chiral flag in the bonds block, only in + # the atoms block. RDKit ignores the atoms chiral flag as per the + # spec. When reading SDF (e.g. from PubChem/PDB), + # we may need to have RDKit assign stereo from coordinates, see: + # https://sourceforge.net/p/rdkit/mailman/message/34399371/ + ps = Chem.RemoveHsParameters() + # a user reported PDBbind Mol Blocks to have hcount=1 for Hs, + # which adds a query to the RDKit H atom and then Chem.RemoveHs + # does not remove Hs with queries by default + # https://github.com/forlilab/Meeko/issues/62 + # https://github.com/rdkit/rdkit/issues/6615 + ps.removeWithQuery = True + mol_noH = Chem.RemoveHs(mol_no_ignore, ps) # imines (=NH) may become chiral + # stereo imines [H]/N=C keep [H] after RemoveHs() + # H isotopes also kept after RemoveHs() + atomic_num_mol_noH = [atom.GetAtomicNum() for atom in mol_noH.GetAtoms()] + noH_to_H = [] + parents_of_hs = {} + for index, atom in enumerate(mol_no_ignore.GetAtoms()): + if atom.GetAtomicNum() == 1: + continue + for i in range(len(noH_to_H), len(atomic_num_mol_noH)): + if atomic_num_mol_noH[i] > 1: + break + h_atom = mol_noH.GetAtomWithIdx(len(noH_to_H)) + assert h_atom.GetAtomicNum() == 1 + neighbors = h_atom.GetNeighbors() + assert len(neighbors) == 1 + parents_of_hs[len(noH_to_H)] = neighbors[0].GetIdx() + noH_to_H.append("H") + noH_to_H.append(index) + extra_hydrogens = len(atomic_num_mol_noH) - len(noH_to_H) + if extra_hydrogens > 0: + assert set(atomic_num_mol_noH[len(noH_to_H):]) == {1} + for i in range(extra_hydrogens): + h_atom = mol_noH.GetAtomWithIdx(len(noH_to_H)) + assert h_atom.GetAtomicNum() == 1 + neighbors = h_atom.GetNeighbors() + assert len(neighbors) == 1 + parents_of_hs[len(noH_to_H)] = neighbors[0].GetIdx() + noH_to_H.append("H") + + # noH_to_H has the same length as the number of atoms in mol_noH + # and each value is: + # - the index of the corresponding atom in mol, if value is integer + # - an hydrogen, if value is "H" + # now, we need to replace those "H" with integers + # "H" occur with stereo imine (e.g. [H]/N=C) and heavy Hs (e.g. [2H]) + hs_by_parent = {} + for hidx, pidx in parents_of_hs.items(): + hs_by_parent.setdefault(pidx, []) + hs_by_parent[pidx].append(hidx) + for pidx, hidxs in hs_by_parent.items(): + siblings_of_h = [ + atom + for atom in mol_no_ignore.GetAtomWithIdx(noH_to_H[pidx]).GetNeighbors() + if atom.GetAtomicNum() == 1 + ] + sortidx = [ + i + for i, j in sorted( + list(enumerate(siblings_of_h)), key=lambda x: x[1].GetIdx() + ) + ] + if len(hidxs) == len(siblings_of_h): + # This is the easy case, just map H to each other in the order they appear + for i, hidx in enumerate(hidxs): + noH_to_H[hidx] = siblings_of_h[sortidx[i]].GetIdx() + elif len(hidxs) < len(siblings_of_h): + # check hydrogen isotopes + sibling_isotopes = [ + siblings_of_h[sortidx[i]].GetIsotope() + for i in range(len(siblings_of_h)) + ] + # molnoH_isotopes = [mol_noH.GetAtomWithIdx(hidx) for hidx in hidxs] + matches = [] + for i, sibling_isotope in enumerate(sibling_isotopes): + for hidx in hidxs[len(matches):]: + if mol_noH.GetAtomWithIdx(hidx).GetIsotope() == sibling_isotope: + matches.append(i) + break + if len(matches) != len(hidxs): + raise RuntimeError( + "Number of matched isotopes %d differs from query Hs: %d" + % (len(matches), len(hidxs)) + ) + for hidx, i in zip(hidxs, matches): + noH_to_H[hidx] = siblings_of_h[sortidx[i]].GetIdx() + else: + raise RuntimeError( + "nr of Hs in mol_noH bonded to an atom exceeds nr of Hs in mol_no_ignore" + ) + + smiles = Chem.MolToSmiles(mol_noH) + order_string = mol_noH.GetProp("_smilesAtomOutputOrder") + order_string = order_string.replace(",]", "]") # remove trailing comma + order = json.loads(order_string) # mol_noH to smiles + order = list(np.argsort(order)) + order = {noH_to_H[i]: order[i] + 1 for i in range(len(order))} # 1-index + return smiles, order + + # region Ring Construction + def _is_ring_aromatic(self, ring_atom_indices: List[Tuple[int]]): + """ + Determines whether a ring is aromatic. + + Parameters + ---------- + ring_atom_indices: the atom indices in the ring. + + Returns + ------- + A boolean indicating whether this ring is aromatic. + """ + for atom_idx1, atom_idx2 in self.get_bonds_in_ring(ring_atom_indices): + bond = self.mol.GetBondBetweenAtoms(atom_idx1, atom_idx2) + if not bond.GetIsAromatic(): + return False + return True + + @staticmethod + def _construct_old_graph(atom_list: List[Atom]): + """ + To support older implementations of helper functions in Meeko, takes a list of atoms and uses it to create a + list of each atom's graph value, where the index of a graph in the list corresponds to the atom's atom_index. + + Parameters + ---------- + atom_list: list[Atom] + A list of populated Atom objects. + + Returns + ------- + A dict mapping from atom index to lists of ints, where each list of ints represents the bonds from that + atom index to other atom indices. + """ + output_graph = {} + for atom in atom_list: + output_graph[atom.index] = atom.graph + return output_graph + + def perceive_rings(self, keep_chorded_rings: bool, keep_equivalent_rings: bool): + """ + Uses Hanser-Jauffret-Kaufmann exhaustive ring detection to find the rings in the molecule + + Parameters + ---------- + keep_chorded_rings: bool + Indicates whether we want to keep chorded rings + keep_equivalent_rings: bool + Indicates whether we want to keep equivalent rings + + Returns + ------- + None + """ + old_graph = self._construct_old_graph(self.atoms) + hjk_ring_detection = utils.HJKRingDetection(old_graph) + rings = hjk_ring_detection.scan(keep_chorded_rings, keep_equivalent_rings) + for ring_atom_indices in rings: + ring_to_add = Ring(ring_atom_indices) + if self._is_ring_aromatic(ring_atom_indices): + ring_to_add.is_aromatic = True + for atom_idx in ring_atom_indices: + # TODO: add it to some sort of atom to ring id tracking thing -> add to atom data structure + ring_to_add.graph = self._recursive_graph_walk( + atom_idx, collected=[], exclude=list(ring_atom_indices) + ) + self.rings[ring_atom_indices] = ring_to_add + return + + # endregion + + def get_conformer_with_modified_positions(self, new_atom_positions): + """ + Gets a conformer with the specified new atom positions. + We operate on one conformer at a time because SetTerminalAtomPositions acts on all conformers of a molecule, + and we do not want to guarantee that all conformers require the same set of terminal atoms to be updated. + + Parameters + ---------- + new_atom_positions: + The new atom positions we want to use. + + Returns + ------- + new_conformer: + A new conformer with the input new atom positions. + """ + new_mol = Chem.Mol(self.mol) + new_conformer = Chem.Conformer(self.mol.GetConformer()) + is_set_list = [False] * self.mol.GetNumAtoms() + for atom_index, new_position in new_atom_positions.items(): + new_conformer.SetAtomPosition(atom_index, new_position) + is_set_list[atom_index] = True + new_mol.RemoveAllConformers() + new_mol.AddConformer(new_conformer, assignId=True) + for atom_index, is_set in enumerate(is_set_list): + if not is_set and new_mol.GetAtomWithIdx(atom_index).GetAtomicNum() == 1: + neighbors = new_mol.GetAtomWithIdx(atom_index).GetNeighbors() + if len(neighbors) != 1: + raise RuntimeError("Expected H to have one neighbors") + Chem.SetTerminalAtomCoords(new_mol, atom_index, neighbors[0].GetIdx()) + return new_conformer + + def get_mol_with_modified_positions(self, new_atom_positions_list=None): + """ + Modifies the stored RDKit Mol to a new set of atom positions, either those provided or the ones stored in + self.modified_atom_positions, and returns the modified Mol object. + + Parameters + ---------- + new_atom_positions_list: + New atom positions to add to the RDKit Mol object. + + Returns + ------- + new_mol: rdkit.Chem.rdchem.Mol + A new RDKit Mol object with conformers that have the desired new atom positions. + """ + if new_atom_positions_list is None: + new_atom_positions_list = self.modified_atom_positions + new_mol = Chem.Mol(self.mol) + new_mol.RemoveAllConformers() + for new_atom_positions in new_atom_positions_list: + conformer = self.get_conformer_with_modified_positions(new_atom_positions) + new_mol.AddConformer(conformer, assignId=True) + return new_mol + + def get_num_mol_atoms(self): + """ + Gets the number of atoms in the RDKit Mol object. + + Returns + ------- + Number of atoms in the RDKit Mol object. + """ + return self.mol.GetNumAtoms() + + def get_equivalent_atoms(self): + """ + + Returns + ------- + + """ + return list(Chem.CanonicalRankAtoms(self.mol, breakTies=False)) + + @staticmethod + def get_symmetries_for_rmsd(mol, max_matches=17): + mol_noHs = Chem.RemoveHs(mol) + matches = mol.GetSubstructMatches( + mol_noHs, uniquify=False, maxMatches=max_matches + ) + if len(matches) == max_matches: + if mol.HasProp("_Name"): + molname = mol.GetProp("_Name") + else: + molname = "" + print( + "warning: found the maximum nr of matches (%d) in RDKitMolSetup.get_symmetries_for_rmsd" + % max_matches + ) + print( + 'Maybe this molecule is "too" symmetric? %s' % molname, + Chem.MolToSmiles(mol_noHs), + ) + return matches + + @staticmethod + def has_implicit_hydrogens(mol): + # based on needsHs from RDKit's AddHs.cpp + for atom in mol.GetAtoms(): + nr_H_neighbors = 0 + for neighbor in atom.GetNeighbors(): + nr_H_neighbors += int(neighbor.GetAtomicNum() == 1) + if atom.GetTotalNumHs(includeNeighbors=False) > nr_H_neighbors: + return True + return False + + @staticmethod + def from_json(obj): + """ + Takes an object and attempts to decode it into an RDKitMoleculeSetup object. + + Parameters + ---------- + obj: Object + This can be any object, but it should be a dictionary generated by deserializing a JSON of an + RDKitMoleculeSetup object. + + Returns + ------- + If the input is a dictionary corresponding to a RDKitMoleculeSetup, will return an RDKitMoleculeSetup with data + populated from the dictionary. Otherwise, returns the input object. + """ + # if the input object is not a dict, we know that it will not be parsable and is unlikely to be usable or + # safe data, so we should ignore it. + if type(obj) is not dict: + return obj + + # checks that all the keys we expect are in the object dictionary as a safety measure + expected_molsetup_keys = { + "mol", + "modified_atom_positions", + "dihedral_interactions", + "dihedral_partaking_atoms", + "dihedral_labels", + "atom_to_ring_id", + "ring_corners", + "rmsd_symmetry_indices", + } + for key in expected_molsetup_keys: + if key not in obj.keys(): + return obj + + # Uses the MoleculeSetup json decoder to get a base MoleculeSetup: + base_molsetup = MoleculeSetup.from_json(obj) + # Constructs an RDKitMoleculeSetup object + rdkit_molsetup = RDKitMoleculeSetup() + # Restores attributes from the MoleculeSetup. There may be a more elegant way to do this + for key, value in base_molsetup.__dict__.items(): + rdkit_molsetup.__dict__[key] = deepcopy(value) + # Restores RDKitMoleculeSetup-specific attributes from the json dict + rdkit_molsetup.mol = rdkit_mol_from_json(obj["mol"]) + rdkit_molsetup.modified_atom_positions = [ + int(x) for x in obj["modified_atom_positions"] + ] + # TODO: Dihedral decoding may need another look + rdkit_molsetup.dihedral_interactions = obj["dihedral_interactions"] + rdkit_molsetup.dihedral_partaking_atoms = obj["dihedral_partaking_atoms"] + rdkit_molsetup.dihedral_labels = obj["dihedral_labels"] + rdkit_molsetup.atom_to_ring_id = { + int(k): [string_to_tuple(t) for t in v] + for k, v in obj["atom_to_ring_id"].items() + } + rdkit_molsetup.ring_corners = obj["ring_corners"] + rdkit_molsetup.rmsd_symmetry_indices = [ + string_to_tuple(v) for v in obj["rmsd_symmetry_indices"] + ] + return rdkit_molsetup + + +# endregion + + +# region JSON Serialization and Deserialization +class AtomEncoder(json.JSONEncoder): + """ + JSON Encoder class for Atom dataclass. + """ + + def default(self, obj): + """ + Overrides the default JSON encoder for data structures for Atom objects. + + Parameters + ---------- + obj: object + Can take any object as input, but will only create the Atom JSON format for Atom objects. + For all other objects will return the default JSON encoding. + + Returns + ------- + A JSON serializable object that represents the Atom class or the default JSONEncoder output for an + object. + """ + if isinstance(obj, Atom): + return { + "index": obj.index, + "pdbinfo": obj.pdbinfo, + "charge": obj.charge, + "coord": obj.coord.tolist(), # converts coord from numpy array to lists + "atomic_num": obj.atomic_num, + "atom_type": obj.atom_type, + "is_ignore": obj.is_ignore, + "graph": obj.graph, + "interaction_vectors": [v.tolist() for v in obj.interaction_vectors], + "is_dummy": obj.is_dummy, + "is_pseudo_atom": obj.is_pseudo_atom, + } + return json.JSONEncoder.default(self, obj) + + +class BondEncoder(json.JSONEncoder): + """ + JSON Encoder class for Bond dataclass. + """ + + def default(self, obj): + """ + Overrides the default JSON encoder for data structures for Bond objects. + + Parameters + ---------- + obj: object + Can take any object as input, but will only create the Bond JSON format for Bond objects. + For all other objects will return the default JSON encoding. + + Returns + ------- + A JSON serializable object that represents the Bond class or the default JSONEncoder output for an + object. + """ + if isinstance(obj, Bond): + return { + "canon_id": tuple_to_string(obj.canon_id), + "index1": obj.index1, + "index2": obj.index2, + "rotatable": obj.rotatable, + } + return json.JSONEncoder.default(self, obj) + + +class RingEncoder(json.JSONEncoder): + """ + JSON Encoder class for Ring dataclass. + """ + + def default(self, obj): + """ + Overrides the default JSON encoder for data structures for Ring objects. + + Parameters + ---------- + obj: object + Can take any object as input, but will only create the Ring JSON format for Ring objects. + For all other objects will return the default JSON encoding. + + Returns + ------- + A JSON serializable object that represents the Ring class or the default JSONEncoder output for an + object. + """ + if isinstance(obj, Ring): + return { + "ring_id": tuple_to_string(obj.ring_id), + "corner_flip": obj.corner_flip, + "graph": obj.graph, + "is_aromatic": obj.is_aromatic, + } + return json.JSONEncoder.default(self, obj) + + +class RestraintEncoder(json.JSONEncoder): + """ + JSON Encoder class for Restraint dataclass. + """ + + def default(self, obj): + """ + Overrides the default JSON encoder for data structures for Restraint objects. + + Parameters + ---------- + obj: object + Can take any object as input, but will only create the Restraint JSON format for Restraint objects. + For all other objects will return the default JSON encoding. + + Returns + ------- + A JSON serializable object that represents the Restraint class or the default JSONEncoder output for an + object. + """ + if isinstance(obj, Restraint): + return { + "atom_index": obj.atom_index, + "target_coords": tuple_to_string(obj.target_coords), + "kcal_per_angstrom_square": obj.kcal_per_angstrom_square, + "delay_angstroms": obj.delay_angstroms, + } + return json.JSONEncoder.default(self, obj) + + +class MoleculeSetupEncoder(json.JSONEncoder): + """ + JSON Encoder class for MoleculeSetup objects. + """ + + atom_encoder = AtomEncoder() + bond_encoder = BondEncoder() + ring_encoder = RingEncoder() + restraint_encoder = RestraintEncoder() + + def default(self, obj): + """ + Overrides the default JSON encoder for data structures for MoleculeSetup objects. + + Parameters + ---------- + obj: object + Can take any object as input, but will only create the MoleculeSetup JSON format for MoleculeSetup objects. + For all other objects will return the default JSON encoding. + + Returns + ------- + A JSON serializable object that represents the MoleculeSetup class or the default JSONEncoder output for an + object. + """ + output_dict = {} + if isinstance(obj, MoleculeSetup): + output_dict = { + "name": obj.name, + "is_sidechain": obj.is_sidechain, + "pseudoatom_count": obj.pseudoatom_count, + "atoms": [self.atom_encoder.default(x) for x in obj.atoms], + "bond_info": { + tuple_to_string(k): self.bond_encoder.default(v) + for k, v in obj.bond_info.items() + }, + "rings": { + tuple_to_string(k): self.ring_encoder.default(v) + for k, v in obj.rings.items() + }, + "ring_closure_info": obj.ring_closure_info.__dict__, + "rotamers": obj.rotamers, + "atom_params": obj.atom_params, + "restraints": [ + self.restraint_encoder.default(x) for x in obj.restraints + ], + "flexibility_model": obj.flexibility_model, + } + # Addressing some flexibility model-specific structures. + if "rigid_body_connectivity" in obj.flexibility_model: + new_rigid_body_conn_dict = { + tuple_to_string(k): v + for k, v in obj.flexibility_model["rigid_body_connectivity"].items() + } + output_dict["flexibility_model"] = { + k: ( + v + if k != "rigid_body_connectivity" + else new_rigid_body_conn_dict + ) + for k, v in obj.flexibility_model.items() + } + # If it's an RDKitMoleculeSetup, adds the RDKitMoleculeSetup attributes + if isinstance(obj, RDKitMoleculeSetup): + output_dict["mol"] = rdMolInterchange.MolToJSON(obj.mol) + output_dict["modified_atom_positions"] = obj.modified_atom_positions + output_dict["dihedral_interactions"] = obj.dihedral_interactions + output_dict["dihedral_partaking_atoms"] = obj.dihedral_partaking_atoms + output_dict["dihedral_labels"] = obj.dihedral_labels + output_dict["atom_to_ring_id"] = obj.atom_to_ring_id + output_dict["ring_corners"] = obj.ring_corners + output_dict["rmsd_symmetry_indices"] = obj.rmsd_symmetry_indices + # If nothing is in the dict, then none of the possible object types for this encoder matched and we should + # return the default JSON encoder. + if len(output_dict) == 0: + return json.JSONEncoder.default(self, obj) + else: + return output_dict + + +# endregion diff --git a/src/moldrug/meeko_light/preparation.py b/src/moldrug/meeko_light/preparation.py new file mode 100644 index 0000000..02fd0cb --- /dev/null +++ b/src/moldrug/meeko_light/preparation.py @@ -0,0 +1,584 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Meeko preparation +# + +from inspect import signature +import json +from os import linesep as eol +import pathlib +import warnings + +from rdkit import Chem + +from .molsetup import Bond +from .molsetup import RDKitMoleculeSetup +from .atomtyper import AtomTyper +from .bondtyper import BondTyperLegacy +from .macrocycle import FlexMacrocycle +from .macrocycle import DEFAULT_MIN_RING_SIZE as M_DEFAULT_MIN_RING_SIZE +from .macrocycle import DEFAULT_MAX_RING_SIZE as M_DEFAULT_MAX_RING_SIZE +from .macrocycle import DEFAULT_DOUBLE_BOND_PENALTY as M_DEFAULT_DOUBLE_BOND_PENALTY +from .flexibility import get_flexibility_model +from .flexibility import update_closure_atoms +from .flexibility import merge_terminal_atoms +from .writer import PDBQTWriterLegacy + +pkg_dir = pathlib.Path(__file__).parents[0] +params_dir = pkg_dir / "data" / "params" +# the above is controversial, see +# https://stackoverflow.com/questions/6028000/how-to-read-a-static-file-from-inside-a-python-package + +# DeprecationWarning is not displayed by default +warnings.filterwarnings("default", category=DeprecationWarning) + + +class MoleculePreparation: + """ + Attributes + ---------- + deprecated_setup_access: + merge_these_atom_types: tuple + hydrate: bool + flexible_amides: bool + rigid_macrocycles: bool + min_ring_size: int + max_ring_size: int + keep_chorded_rings: bool + keep_equivalent_rings: bool + double_bond_penalty: float + macrocycle_allow_A: bool + rigidify_bonds_smarts: list + rigidify_bonds_indices: list + + input_atom_params: + load_atom_params: + add_atom_types: + + atom_params: + """ + + packaged_params = {} + for path in params_dir.glob("*.json"): # e.g. data/params/ad4_types.json + name = path.with_suffix("").name # e.g. "ad4_types" + packaged_params[name] = path + + def __init__( + self, + merge_these_atom_types=("H",), + hydrate=False, + flexible_amides=False, + rigid_macrocycles=False, + min_ring_size=M_DEFAULT_MIN_RING_SIZE, + max_ring_size=M_DEFAULT_MAX_RING_SIZE, + keep_chorded_rings=False, + keep_equivalent_rings=False, + double_bond_penalty=M_DEFAULT_DOUBLE_BOND_PENALTY, + macrocycle_allow_A=False, + rigidify_bonds_smarts=None, + rigidify_bonds_indices=None, + input_atom_params=None, + load_atom_params="ad4_types", + add_atom_types=(), + input_offatom_params=None, + load_offatom_params=None, + charge_model="gasteiger", + add_index_map=False, + remove_smiles=False, + ): + """ + + Parameters + ---------- + merge_these_atom_types + hydrate + flexible_amides + rigid_macrocycles + min_ring_size + max_ring_size + keep_chorded_rings + keep_equivalent_rings + double_bond_penalty + macrocycle_allow_A + rigidify_bonds_smarts + rigidify_bonds_indices + input_atom_params + load_atom_params + add_atom_types + input_offatom_params + load_offatom_params + charge_model + add_index_map + remove_smiles + """ + + self.deprecated_setup_access = None + self.merge_these_atom_types = merge_these_atom_types + self.hydrate = hydrate + self.flexible_amides = flexible_amides + self.rigid_macrocycles = rigid_macrocycles + self.min_ring_size = min_ring_size + self.max_ring_size = max_ring_size + self.keep_chorded_rings = keep_chorded_rings + self.keep_equivalent_rings = keep_equivalent_rings + self.double_bond_penalty = double_bond_penalty + self.macrocycle_allow_A = macrocycle_allow_A + if rigidify_bonds_smarts is None: + rigidify_bonds_smarts = [] + self.rigidify_bonds_smarts = rigidify_bonds_smarts + if rigidify_bonds_indices is None: + rigidify_bonds_indices = [] + self.rigidify_bonds_indices = rigidify_bonds_indices + + self.input_atom_params = input_atom_params + self.load_atom_params = load_atom_params + self.add_atom_types = add_atom_types + + self.atom_params = self.get_atom_params( + input_atom_params, load_atom_params, add_atom_types, self.packaged_params + ) + + if load_offatom_params is not None: + raise NotImplementedError("load_offatom_params not implemented") + self.load_offatom_params = load_offatom_params + + allowed_charge_models = ["gasteiger", "zero"] + if charge_model not in allowed_charge_models: + raise ValueError( + "unrecognized charge_model: %s, allowed options are: %s" + % (charge_model, allowed_charge_models) + ) + + self.charge_model = charge_model + + self.dihedral_model = None + self.dihedral_params = [] + + self.add_index_map = add_index_map + self.remove_smiles = remove_smiles + + self._bond_typer = BondTyperLegacy() + self._macrocycle_typer = FlexMacrocycle( + self.min_ring_size, + self.max_ring_size, + self.double_bond_penalty, + allow_break_atype_A=self.macrocycle_allow_A, + ) + self._classes_setup = {Chem.rdchem.Mol: RDKitMoleculeSetup} + + if input_offatom_params is None: + self.offatom_params = {} + else: + self.offatom_params = input_offatom_params + + if keep_chorded_rings and not keep_equivalent_rings: + warnings.warn( + "keep_equivalent_rings=False ignored because keep_chorded_rings=True", + RuntimeWarning, + ) + + @classmethod + def from_config(cls, config): + """ + + Parameters + ---------- + config + + Returns + ------- + + """ + expected_keys = cls.get_defaults_dict().keys() + bad_keys = [k for k in config if k not in expected_keys] + if len(bad_keys) > 0: + err_msg = ( + "unexpected keys in MoleculePreparation.from_config():" + eol + ) + for key in bad_keys: + err_msg += " - %s" % key + eol + raise ValueError(err_msg) + p = cls(**config) + return p + + def calc_flex( + self, + setup, + root_atom_index=None, + not_terminal_atoms=None, + delete_ring_bonds=None, + glue_pseudo_atoms=None, + ): + """ + + Parameters + ---------- + setup + root_atom_index + not_terminal_atoms + delete_ring_bonds + glue_pseudo_atoms + + Returns + ------- + + """ + if not_terminal_atoms is None: + not_terminal_atoms = [] + if delete_ring_bonds is None: + delete_ring_bonds = [] + # 5. break macrocycles into open/linear form + if self.rigid_macrocycles: + break_combo_data = None + bonds_in_rigid_rings = set() + # every ring is rigid without macrocycle option + for ring in setup.rings: + for bond in setup.get_bonds_in_ring(ring): + bonds_in_rigid_rings.add(bond) + else: + break_combo_data, bonds_in_rigid_rings = ( + self._macrocycle_typer.search_macrocycle(setup, delete_ring_bonds) + ) + + # This must be done before calling get_flexibility_model + for bond in bonds_in_rigid_rings: + setup.bond_info[bond].rotatable = False + + flex_model, bonds_to_break = get_flexibility_model( + setup, root_atom_index, break_combo_data + ) + + # disasble rotatable bonds that rotate nothing (e.g. -CH3 without H) + # but glue atoms (i.e. CG) are manually marked as non terminal (by + # passing them in the `glue_atoms` list) to guarantee that the bond + # to a CG atom is rotatable and the G pseudo rotates + glue_atoms = [] + for pair in bonds_to_break: + for index in pair: + glue_atoms.append(index) + merge_terminal_atoms(flex_model, not_terminal_atoms + glue_atoms) + + # bond to a terminal atom, or in a ring that isn't flexible + actual_rotatable = [v for k, v in flex_model["rigid_body_connectivity"].items()] + actual_rotatable.extend(bonds_to_break) + for bond in setup.bond_info: + if bond not in actual_rotatable: + setup.bond_info[bond].rotatable = False + + # calculate torsions that would be rotatable without macrocycle breaking + if break_combo_data is not None and len(break_combo_data["bond_break_combos"]): + ring_bonds = [] + for ring in setup.rings: + for bond in setup.get_bonds_in_ring(ring): + ring_bonds.append(bond) + flex_model["torsions_org"] = 0 + for bond in setup.bond_info: + if setup.bond_info[bond].rotatable and bond not in ring_bonds: + flex_model["torsions_org"] += 1 + + setup.flexibility_model = flex_model + + # add G pseudo atoms and set CG types + update_closure_atoms(setup, bonds_to_break, glue_pseudo_atoms) + + return + + @staticmethod + def get_atom_params( + input_atom_params, load_atom_params, add_atom_types, packaged_params + ): + """ + + Parameters + ---------- + input_atom_params + load_atom_params + add_atom_types + packaged_params + + Returns + ------- + + """ + atom_params = {} + if isinstance(load_atom_params, str): + load_atom_params = [load_atom_params] + elif load_atom_params is None: + load_atom_params = () + for name in load_atom_params: + filename = None + if name in packaged_params: + filename = packaged_params[name] + elif name.endswith(".json"): + filename = name + else: + msg = ( + "names passed to 'load_atom_params' need to suffixed with .json" + + eol + ) + msg += ( + "or be the unsuffixed basename of a JSON file in %s." + % str(params_dir) + + eol + ) + msg += "name was %s" % name + raise ValueError(msg) + if filename is not None: + with open(filename) as f: + d = json.load(f) + overlapping_groups = set(atom_params).intersection(set(d)) + if len(overlapping_groups): + msg = "overlapping parameter groups: %s" % str(overlapping_groups) + raise ValueError(msg) + atom_params.update(d) + + if input_atom_params is not None: + d = json.loads(json.dumps(input_atom_params)) + overlapping_groups = set(atom_params).intersection(set(d)) + if len(overlapping_groups): # todo: remove duplicated code + msg = "overlapping parameter groups: %s" % str(overlapping_groups) + raise ValueError(msg) + params_set_here = set() + for rows in input_atom_params.values(): + for row in rows: + for param_name in row: + if param_name not in ["smarts", "comment", "IDX"]: + params_set_here.add(param_name) + params_set_before = set() + for rows in atom_params.values(): + for row in rows: + for param_name in row: + if param_name not in ["smarts", "comment", "IDX"]: + params_set_before.add(param_name) + overlap = params_set_before.intersection(params_set_here) + if len(overlap): + msg = f"input_atom_params {overlap} also set by one or more of {load_atom_params}\n" + msg += "consider setting load_atom_params=None" + raise ValueError( + f"input_atom_params {overlap} also set by one or more of {load_atom_params}" + ) + + atom_params.update(d) + + if len(add_atom_types) > 0: + group_keys = list(atom_params.keys()) + if len(group_keys) != 1: + msg = "add_atom_types is usable only when there is one group of parameters" + msg += ", but there are %d groups: %s" % ( + len(group_keys), + str(group_keys), + ) + raise RuntimeError(msg) + key = group_keys[0] + atom_params[key].extend(add_atom_types) + + return atom_params + + @property + def setup(self): + """ + + Returns + ------- + + """ + msg = "MoleculePreparation.setup is deprecated in Meeko v0.5." + msg += ( + " MoleculePreparation.prepare() returns a list of MoleculeSetup instances." + ) + warnings.warn(msg, DeprecationWarning) + if len(self.deprecated_setup_access) > 1: + raise RuntimeError( + "got multiple setups, use new api: molsetup_list = mk_prep(mol)" + ) + + return self.deprecated_setup_access[0] + + @classmethod + def get_defaults_dict(cls): + """ + + Returns + ------- + + """ + defaults = {} + sig = signature(cls) + for key in sig.parameters: + defaults[key] = sig.parameters[key].default + return defaults + + def __call__(self, *args): + return self.prepare(*args) + + def prepare( + self, + mol, + root_atom_index=None, + not_terminal_atoms=None, + delete_ring_bonds=None, + glue_pseudo_atoms=None, + conformer_id=-1, + ): + """ + Create an RDKitMoleculeSetup from an RDKit Mol object. + + Parameters + ---------- + mol: rdkit.Chem.rdchem.Mol + An RDKit Mol with explicit hydrogens and 3D coordinates. + root_atom_index: int + Used to set ROOT of torsion tree instead of searching. + not_terminal_atoms: list + Makes bonds with terminal atoms rotatable (e.g. C-Alpha carbon in flexres). + delete_ring_bonds: list[tuple[int, int]] + Bonds deleted for macrocycle flexibility. Each bond is a tuple of two ints (atom 0-indices). + glue_pseudo_atoms: dict + Mapping from parent atom indices to coordinates. + conformer_id: int + + Returns + ------- + setups: list[RDKitMoleculeSetup] + Returns a list of generated RDKitMoleculeSetups + """ + + if not_terminal_atoms is None: + not_terminal_atoms = [] + if delete_ring_bonds is None: + delete_ring_bonds = [] + if glue_pseudo_atoms is None: + glue_pseudo_atoms = {} + mol_type = type(mol) + if mol_type not in self._classes_setup: + raise TypeError( + "Molecule is not an instance of supported types: %s" % type(mol) + ) + setup_class = self._classes_setup[mol_type] + setup = setup_class.from_mol( + mol, + keep_chorded_rings=self.keep_chorded_rings, + keep_equivalent_rings=self.keep_equivalent_rings, + assign_charges=self.charge_model == "gasteiger", + conformer_id=conformer_id, + ) + + self.check_external_ring_break(setup, delete_ring_bonds, glue_pseudo_atoms) + + # 1. assign atom params + AtomTyper.type_everything( + setup, + self.atom_params, + self.charge_model, + self.offatom_params, + self.dihedral_params, + ) + + # merge hydrogens (or any terminal atoms) + indices = set() + for atype_to_merge in self.merge_these_atom_types: + for atom in setup.atoms: + if atom.atom_type == atype_to_merge: + indices.add(atom.index) + setup.merge_terminal_atoms(indices) + + # 3. assign bond types + # - all single bonds rotatable except some amides and SMARTS rigidification + # - macrocycle code breaks rings only at rotatable bonds + # - bonds in rigid rings are set as non-rotatable after the flex_model is built + self._bond_typer( + setup, + self.flexible_amides, + self.rigidify_bonds_smarts, + self.rigidify_bonds_indices, + ) + + self.calc_flex( + setup, + root_atom_index, + not_terminal_atoms, + delete_ring_bonds, + glue_pseudo_atoms, + ) + + # for a gentle introduction of the new API + self.deprecated_setup_access = [setup] + + return [setup] + + @staticmethod + def check_external_ring_break(molsetup, break_ring_bonds, glue_pseudo_atoms): + """ + + Parameters + ---------- + molsetup: RDKitMoleculeSetup + break_ring_bonds: + glue_pseudo_atoms: dict + + Returns + ------- + None + + Raises + ------ + ValueError: + If bonds are missing from the MoleculeSetup, if glue_pseudo_atoms is missing certain atom indices, and if + there is an incorrect number of coordinates in glue_pseudo_atoms. + """ + for index1, index2 in break_ring_bonds: + has_bond = Bond.get_bond_id(index1, index2) in molsetup.bond_info + if not has_bond: + raise ValueError("bond (%d, %d) not in molsetup" % (index1, index2)) + for index in (index1, index2): + if index not in glue_pseudo_atoms: + raise ValueError("missing glue pseudo for atom %d" % index) + xyz = glue_pseudo_atoms[index] + if len(xyz) != 3: + raise ValueError( + "expected 3 coordinates (got %d) for glue pseudo of atom %d" + % (len(xyz), index) + ) + return + + def write_pdbqt_string(self): + """ + Writes a PDBQT string. Deprecated in Meeko v0.5. + + Returns + ------- + + """ + msg = "MoleculePreparation.write_pdbqt_string() is deprecated in Meeko v0.5." + msg += " Pass the MoleculeSetup instance to PDBQTWriterLegacy.write_string()." + msg += ( + " MoleculePreparation.prepare() returns a list of MoleculeSetup instances." + ) + warnings.warn(msg, DeprecationWarning) + pdbqt_string, is_ok, err_msg = PDBQTWriterLegacy.write_string(self.setup) + if not is_ok: + msg = "Cannot generate PDBQT, error from PDBQTWriterLegacy:" + eol + msg += err_msg + raise RuntimeError(msg) + return pdbqt_string + + def write_pdbqt_file(self, pdbqt_filename): + """ + Writes out a pdbqt file. Deprecated in Meeko v0.5 + + Parameters + ---------- + pdbqt_filename: str + PDBQT filename to write to + + Returns + ------- + None + """ + warnings.warn( + "MoleculePreparation.write_pdbqt_file() is deprecated since Meeko v0.5", + DeprecationWarning, + ) + with open(pdbqt_filename, "w") as w: + w.write(self.write_pdbqt_string()) diff --git a/src/moldrug/meeko_light/rdkit_mol_create.py b/src/moldrug/meeko_light/rdkit_mol_create.py new file mode 100644 index 0000000..8a560a4 --- /dev/null +++ b/src/moldrug/meeko_light/rdkit_mol_create.py @@ -0,0 +1,495 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Meeko +# + + +from rdkit import Chem +from rdkit.Geometry import Point3D +from rdkit.Chem import AllChem +from io import StringIO +import json + + +def clean_extend(existing_dict, new_row): + nr_rows = [] + for key in existing_dict: + nr_rows.append(len(existing_dict[key])) + if key not in new_row: + existing_dict[key].append(None) + if len(nr_rows) == 0: # existing_dict is empty + nr_rows = 0 + elif len(set(nr_rows)) != 1: + msg = "existing_dict has different nr of items for different attributes" + raise ValueError(msg) + else: + nr_rows = set(nr_rows).pop() + for key, value in new_row.items(): + if key not in existing_dict: + existing_dict[key] = [None] * nr_rows + existing_dict[key].append(value) + return + + +class RDKitMolCreate: + + ambiguous_flexres_choices = { + "HIS": ["HIE", "HID", "HIP"], + "ASP": ["ASP", "ASH"], + "GLU": ["GLU", "GLH"], + "CYS": ["CYS", "CYM"], + "LYS": ["LYS", "LYN"], + "ARG": ["ARG", "ARG_mgltools"], + "ASN": ["ASN", "ASN_mgltools"], + "GLN": ["GLN", "GLN_mgltools"], + } + + flexres = { + "CYS": { + "smiles": "CCS", + "atom_names_in_smiles_order": ["CA", "CB", "SG"], + "h_to_parent_index": {"HG": 2}, + }, + "CYM": { + "smiles": "CC[S-]", + "atom_names_in_smiles_order": ["CA", "CB", "SG"], + "h_to_parent_index": {}, + }, + "ASP": { + "smiles": "CCC(=O)[O-]", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "OD1", "OD2"], + "h_to_parent_index": {}, + }, + "ASH": { + "smiles": "CCC(=O)O", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "OD1", "OD2"], + "h_to_parent_index": {"HD2": 4}, + }, + "GLU": { + "smiles": "CCCC(=O)[O-]", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD", "OE1", "OE2"], + "h_to_parent_index": {}, + }, + "GLH": { + "smiles": "CCCC(=O)O", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD", "OE1", "OE2"], + "h_to_parent_index": {"HE2": 5}, + }, + "PHE": { + "smiles": "CCc1ccccc1", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD1", "CE1", "CZ", "CE2", "CD2"], + "h_to_parent_index": {}, + }, + "HIE": { + "smiles": "CCc1c[nH]cn1", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD2", "NE2", "CE1", "ND1"], + "h_to_parent_index": {"HE2": 4}, + }, + "HID": { + "smiles": "CCc1cnc[nH]1", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD2", "NE2", "CE1", "ND1"], + "h_to_parent_index": {"HD1": 6}, + }, + "HIP": { + "smiles": "CCc1c[nH+]c[nH]1", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD2", "NE2", "CE1", "ND1"], + "h_to_parent_index": {"HE2": 4, "HD1": 6}, + }, + "ILE": { + "smiles": "CC(C)CC", + "atom_names_in_smiles_order": ["CA", "CB", "CG2", "CG1", "CD1"], + "h_to_parent_index": {}, + }, + "LYS": { + "smiles": "CCCCC[NH3+]", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD", "CE", "NZ"], + "h_to_parent_index": {"HZ1": 5, "HZ2": 5, "HZ3": 5}, + }, + "LYN": { + "smiles": "CCCCCN", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD", "CE", "NZ"], + "h_to_parent_index": {"HZ2": 5, "HZ3": 5}, + }, + "LEU": { + "smiles": "CCC(C)C", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD1", "CD2"], + "h_to_parent_index": {}, + }, + "MET": { + "smiles": "CCCSC", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "SD", "CE"], + "h_to_parent_index": {}, + }, + "ASN": { + "smiles": "CCC(=O)N", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "OD1", "ND2"], + "h_to_parent_index": {"HD21": 4, "HD22": 4}, + }, + "ASN_mgltools": { + "smiles": "CCC(=O)N", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "OD1", "ND2"], + "h_to_parent_index": {"1HD2": 4, "2HD2": 4}, + }, + "GLN": { + "smiles": "CCCC(=O)N", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD", "OE1", "NE2"], + "h_to_parent_index": {"HE21": 5, "HE22": 5}, + }, + "GLN_mgltools": { + "smiles": "CCCC(=O)N", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD", "OE1", "NE2"], + "h_to_parent_index": {"1HE2": 5, "2HE2": 5}, + }, + "ARG": { + "smiles": "CCCCNC(N)=[NH2+]", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD", "NE", "CZ", "NH1", "NH2"], + "h_to_parent_index": {"HE": 4, "HH11": 6, "HH12": 6, "HH21": 7, "HH22": 7}, + }, + "ARG_mgltools": { + "smiles": "CCCCNC(N)=[NH2+]", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD", "NE", "CZ", "NH1", "NH2"], + "h_to_parent_index": {"HE": 4, "1HH1": 6, "2HH1": 6, "1HH2": 7, "2HH2": 7}, + }, + "SER": { + "smiles": "CCO", + "atom_names_in_smiles_order": ["CA", "CB", "OG"], + "h_to_parent_index": {"HG": 2}, + }, + "THR": { + "smiles": "CC(C)O", + "atom_names_in_smiles_order": ["CA", "CB", "CG2", "OG1"], + "h_to_parent_index": {"HG1": 3}, + }, + "VAL": { + "smiles": "CC(C)C", + "atom_names_in_smiles_order": ["CA", "CB", "CG1", "CG2"], + "h_to_parent_index": {}, + }, + "TRP": { + "smiles": "CCc1c[nH]c2c1cccc2", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD1", "NE1", "CE2", "CD2", "CE3", "CZ3", "CH2", "CZ2"], + "h_to_parent_index": {"HE1": 4}, + }, + "TYR": { + "smiles": "CCc1ccc(cc1)O", + "atom_names_in_smiles_order": ["CA", "CB", "CG", "CD1", "CE1", "CZ", "CE2", "CD2", "OH"], + "h_to_parent_index": {"HH": 8}, + }, + } + + @classmethod + def from_pdbqt_mol( + cls, + pdbqt_mol, + only_cluster_leads=False, + keep_flexres=False, + ): + # todo: add pseudo-water (W atoms, variable nr each pose) + if only_cluster_leads and len(pdbqt_mol._pose_data["cluster_leads_sorted"]) == 0: + raise RuntimeError("no cluster_leads in pdbqt_mol but only_cluster_leads=True") + mol_list = [] + for mol_index in pdbqt_mol._atom_annotations["mol_index"]: + flexres_id = pdbqt_mol._pose_data["mol_index_to_flexible_residue"][mol_index] + if flexres_id is not None and not keep_flexres: + continue + smiles = pdbqt_mol._pose_data['smiles'][mol_index] + index_map = pdbqt_mol._pose_data['smiles_index_map'][mol_index] + h_parent = pdbqt_mol._pose_data['smiles_h_parent'][mol_index] + atom_idx = pdbqt_mol._atom_annotations["mol_index"][mol_index] + atom_is_flex = [i in pdbqt_mol._atom_annotations["flexible_residue"] for i in atom_idx] + if any(atom_is_flex) and all(atom_is_flex): + is_sidechain = True + elif any(atom_is_flex): + raise ValueError("some (but not all!) atoms of a ligand were parsed as sidechain") + else: + is_sidechain = False + + if smiles is None: # probably a flexible sidechain, but can be another ligand + residue_names = set() + atom_names = [] + for atom in pdbqt_mol.atoms(atom_idx): + residue_names.add(atom[4]) + atom_names.append(atom[2]) + if len(residue_names) == 1: + resname = residue_names.pop() + smiles, index_map, h_parent = cls.guess_flexres_smiles(resname, atom_names) + if smiles is None: # failed guessing smiles for possible flexres + mol_list.append(None) + continue + + if only_cluster_leads: + pose_ids = pdbqt_mol._pose_data["cluster_leads_sorted"] + else: + pose_ids = range(pdbqt_mol._pose_data["n_poses"]) + + mol = Chem.MolFromSmiles(smiles) + mol.SetProp("meeko", json.dumps({"is_sidechain": is_sidechain})) + coordinates_all_poses = [] + for i in pose_ids: + pdbqt_mol._current_pose = i + coordinates = pdbqt_mol.positions(atom_idx) + mol = cls.add_pose_to_mol(mol, coordinates, index_map) + coordinates_all_poses.append(coordinates) + + # add Hs only after all poses are added as conformers + # because Chem.AddHs() will affect all conformers at once + mol = cls.add_hydrogens(mol, coordinates_all_poses, h_parent) + + mol_list.append(mol) + return mol_list + + @classmethod + def guess_flexres_smiles(cls, resname, atom_names): + """ Determine a SMILES string for flexres based on atom names, + as well as the equivalent of smile_index_map and smiles_h_parent + which are written to PDBQT remarks for regular small molecules. + + Args: + resname (str): + + Returns: + smiles: SMILES string starting at C-alpha (excludes most of the backbone) + index_map: list of pairs of integers, first in pair is index in the smiles, + second is index of corresponding atom in atom_names + h_parent: list of pairs of integers, first in pair is index of a heavy atom + in the smiles, second is index of a hydrogen in atom_names. + The hydrogen is bonded to the heavy atom. + """ + + if len(set(atom_names)) != len(atom_names): + return None, None, None + candidate_resnames = cls.ambiguous_flexres_choices.get(resname, [resname]) + for resname in candidate_resnames: + is_match = False + if resname not in cls.flexres: + continue + atom_names_in_smiles_order = cls.flexres[resname]["atom_names_in_smiles_order"] + h_to_parent_index = cls.flexres[resname]["h_to_parent_index"] + expected_names = atom_names_in_smiles_order + list(h_to_parent_index.keys()) + if len(atom_names) != len(expected_names): + continue + nr_matched_atom_names = sum([int(n in atom_names) for n in expected_names]) + if nr_matched_atom_names == len(expected_names): + is_match = True + break + if not is_match: + return None, None, None + else: + smiles = cls.flexres[resname]["smiles"] + index_map = [] + for smiles_index, name in enumerate(atom_names_in_smiles_order): + index_map.append(smiles_index + 1) + index_map.append(atom_names.index(name) + 1) + h_parent = [] + for name, smiles_index in h_to_parent_index.items(): + h_parent.append(smiles_index + 1) + h_parent.append(atom_names.index(name) + 1) + return smiles, index_map, h_parent + + @classmethod + def add_pose_to_mol(cls, mol, ligand_coordinates, index_map): + """add given coordinates to given molecule as new conformer. + Index_map maps order of coordinates to order in smile string + used to generate rdkit mol + + Args: + ligand_coordinates: 2D array of shape (nr_atom, 3). + index_map: list of nr_atom pairs of integers, 1-indexed. + In each pair, the first int is the index in mol, and + the second int is the index in ligand_coordinates + + Raises: + RuntimeError: Will raise error if number of coordinates provided does not + match the number of atoms there should be coordinates for. + """ + + n_atoms = mol.GetNumAtoms() + n_mappings = int(len(index_map) / 2) + conf = Chem.Conformer(n_atoms) + if n_atoms < n_mappings: + raise RuntimeError( + "Number of atom is rdmol {n_atoms} mismatches" + "number of pairs in index map {n_mappings}!".format( + n_atoms=n_atoms, n_mappings=n_mappings)) + coord_is_set = [False] * n_atoms + for i in range(n_mappings): + pdbqt_index = int(index_map[i * 2 + 1]) - 1 + mol_index = int(index_map[i * 2]) - 1 + x, y, z = [float(coord) for coord in ligand_coordinates[pdbqt_index]] + conf.SetAtomPosition(mol_index, Point3D(x, y, z)) + coord_is_set[mol_index] = True + mol.AddConformer(conf, assignId=True) + # some hydrogens (isotopes) may have no coordinate set yet + for i, is_set in enumerate(coord_is_set): + if not is_set: + atom = mol.GetAtomWithIdx(i) + if atom.GetAtomicNum() != 1: + raise RuntimeError("Only H allowed to be in SMILES but not in PDBQT") + neigh = atom.GetNeighbors() + if len(neigh) != 1: + raise RuntimeError("Expected H to have one neighbor") + AllChem.SetTerminalAtomCoords(mol, i, neigh[0].GetIdx()) + return mol + + @staticmethod + def add_hydrogens(mol, coordinates_list, h_parent): + """Add hydrogen atoms to ligand RDKit mol, adjust the positions of + polar hydrogens to match pdbqt + """ + mol = Chem.AddHs(mol, addCoords=True) + conformers = list(mol.GetConformers()) + num_hydrogens = int(len(h_parent) / 2) + for conformer_idx, atom_coordinates in enumerate(coordinates_list): + conf = conformers[conformer_idx] + used_h = [] + for i in range(num_hydrogens): + parent_rdkit_index = h_parent[2 * i] - 1 + h_pdbqt_index = h_parent[2 * i + 1] - 1 + x, y, z = [ + float(coord) for coord in atom_coordinates[h_pdbqt_index] + ] + parent_atom = mol.GetAtomWithIdx(parent_rdkit_index) + candidate_hydrogens = [ + atom.GetIdx() for atom in parent_atom.GetNeighbors() + if atom.GetAtomicNum() == 1 + ] + for h_rdkit_index in candidate_hydrogens: + if h_rdkit_index not in used_h: + break + used_h.append(h_rdkit_index) + conf.SetAtomPosition(h_rdkit_index, Point3D(x, y, z)) + return mol + + @staticmethod + def combine_rdkit_mols(mol_list): + """Combines list of rdkit molecules into a single one + None's are ignored + returns None if input is empty list or all molecules are None + """ + combined_mol = None + props = {} + for mol in mol_list: + if mol is None: + continue + if mol.HasProp("meeko"): + data = json.loads(mol.GetProp("meeko")) + clean_extend(props, data) + if combined_mol is None: # first iteration + combined_mol = mol + else: + combined_mol = Chem.CombineMols(combined_mol, mol) + if len(props) > 0: + combined_mol.SetProp("meeko", json.dumps(props)) + return combined_mol + + @classmethod + def _verify_flexres(cls): + for resname in cls.flexres: + atom_names_in_smiles_order = cls.flexres[resname]["atom_names_in_smiles_order"] + h_to_parent_index = cls.flexres[resname]["h_to_parent_index"] + expected_names = atom_names_in_smiles_order + list(h_to_parent_index.keys()) + if len(expected_names) != len(set(expected_names)): + raise RuntimeError("repeated atom names in cls.flexres[%s]" % resname) + + @classmethod + def add_sandbox_coordinates(cls, dlgstring, rdmol, index_map, h_parent, groupname=None): + # this function does not deal with implicit H, at least not yet + index_map = [i + 1 for i in index_map] # 1-indexing like in PDBQT + h_parent = [i + 1 for i in h_parent] # 1-indexing like in PDBQT + coordinates = [] + energy = {"inter": [], "intra": [], "dlg_pose_idx": []} + is_atom_block = False + for line in dlgstring.split('\n'): + if line.startswith("Pose:") or line.startswith("Extra Pose:"): + if line.startswith("Pose:"): + pose_idx = int(line.split()[1]) + elif line.startswith("Extra Pose:"): + pose_idx = int(line.split()[2]) + if len(coordinates) > 0: + if len(coordinates[-1]) == 0: + # if pose info was missing, just delete data + energy["dlg_pose_idx"].pop(-1) + coordinates.pop(-1) + energy["dlg_pose_idx"].append(pose_idx) + coordinates.append([]) + elif line.startswith("DOCKED: USER (1) Final Intermolecular Energy ="): + energy["inter"].append(float(line.split()[7])) + elif line.startswith("DOCKED: USER (2) Final Total Internal Energy ="): + energy["intra"].append(float(line.split()[8])) + elif line.startswith("DOCKED: @ATOM"): + is_atom_block = True + elif line.startswith("DOCKED: @BOND"): + is_atom_block = False + elif is_atom_block: + fields = line.split() + name = fields[8] + if groupname is None or name == groupname: + x, y, z = float(fields[3]), float(fields[4]), float(fields[5]) + coordinates[-1].append([x, y, z]) + + if not (len(coordinates) == len(energy["inter"]) == len(energy["intra"])): + msg = "parsed energies differs from number of coordinates\n" + msg += "len(coordinates) = %d\n" % len(coordinates) + msg += "len(intra) = %d\n" % len(energy["intra"]) + msg += "len(inter) = %d\n" % len(energy["inter"]) + raise RuntimeError(msg) + + scores = [energy["inter"][i] + energy["intra"][i] for i in range(len(coordinates))] + idxsort = [pair[0] for pair in sorted(enumerate(scores), key=lambda pair: pair[1])] + sorted_coordinates = [] + for index in idxsort: + cls.add_pose_to_mol(rdmol, coordinates[index], index_map) + sorted_coordinates.append(coordinates[index]) + + rdmol = cls.add_hydrogens(rdmol, sorted_coordinates, h_parent) + + for key in energy: + energy[key] = [energy[key][i] for i in idxsort] + return rdmol, energy + + @staticmethod + def write_sd_string(pdbqt_mol, only_cluster_leads=False, keep_flexres=False): + sio = StringIO() + f = Chem.SDWriter(sio) + mol_list = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol, only_cluster_leads, keep_flexres) + failures = [i for i, mol in enumerate(mol_list) if mol is None] + combined_mol = RDKitMolCreate.combine_rdkit_mols(mol_list) + if combined_mol is None: + return "", failures + keys_map_mol_to_pdbqt = { + "free_energy": "free_energies", + "intermolecular_energy": "intermolecular_energies", + "internal_energy": "internal_energies", + "cluster_size": "cluster_size", + "cluster_id": "cluster_id", + "rank_in_cluster": "rank_in_cluster", + } + nr_poses = pdbqt_mol._pose_data["n_poses"] + if only_cluster_leads: + pose_idxs = pdbqt_mol._pose_data["cluster_leads_sorted"] + else: + pose_idxs = list(range(nr_poses)) + + available_properties = {} + for key_in_mol, key_in_pdbqt in keys_map_mol_to_pdbqt.items(): + if len(pdbqt_mol._pose_data[key_in_pdbqt]) == nr_poses: + available_properties[key_in_mol] = key_in_pdbqt + mol_level_data = json.loads(combined_mol.GetProp("meeko")) + if pdbqt_mol.name is not None: + combined_mol.SetProp("_Name", pdbqt_mol.name) + for conformer in combined_mol.GetConformers(): + i = conformer.GetId() + j = pose_idxs[i] + conformer_data = json.loads(json.dumps(mol_level_data)) + for (key_in_mol, key_in_pdbqt) in available_properties.items(): + if key_in_mol in conformer_data: + msg = "key %s conflict between combined_mol and write_sd_string" % key_in_mol + raise NotImplementedError(msg) + conformer_data[key_in_mol] = pdbqt_mol._pose_data[key_in_pdbqt][j] + if len(conformer_data): + combined_mol.SetProp("meeko", json.dumps(conformer_data)) + f.write(combined_mol, i) + f.close() + output_string = sio.getvalue() + return output_string, failures diff --git a/src/moldrug/meeko_light/utils/__init__.py b/src/moldrug/meeko_light/utils/__init__.py new file mode 100644 index 0000000..aea5294 --- /dev/null +++ b/src/moldrug/meeko_light/utils/__init__.py @@ -0,0 +1,5 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Raccoon utils +# diff --git a/src/moldrug/meeko_light/utils/autodock4_atom_types_elements.py b/src/moldrug/meeko_light/utils/autodock4_atom_types_elements.py new file mode 100644 index 0000000..664c1c4 --- /dev/null +++ b/src/moldrug/meeko_light/utils/autodock4_atom_types_elements.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Conversion AutoDock4 to elment +# + +autodock4_atom_types_elements = { + 'H': 'H', + 'HD': 'H', + 'HS': 'H', + 'C': 'C', + 'A': 'C', + 'N': 'N', + 'NA': 'N', + 'NS': 'N', + 'OA': 'O', + 'OS': 'O', + 'F': 'F', + 'Mg': 'Mg', + 'MG': 'Mg', + 'P': 'P', + 'SA': 'S', + 'S': 'S', + 'Cl': 'Cl', + 'CL': 'Cl', + 'Ca': 'Ca', + 'CA': 'Ca', + 'Mn': 'Mn', + 'MN': 'Mn', + 'Fe': 'Fe', + 'FE': 'Fe', + 'Zn': 'Zn', + 'ZN': 'Zn', + 'Br': 'Br', + 'BR': 'Br', + 'I': 'I', + 'Si': 'Si', + 'At': 'At', + 'G0': 'C', + 'G1': 'C', + 'G2': 'C', + 'G3': 'C', + 'CG0': 'C', + 'CG1': 'C', + 'CG2': 'C', + 'CG3': 'C', + 'W': 'O' +} diff --git a/src/moldrug/meeko_light/utils/covalent_radius_table.py b/src/moldrug/meeko_light/utils/covalent_radius_table.py new file mode 100644 index 0000000..f96a395 --- /dev/null +++ b/src/moldrug/meeko_light/utils/covalent_radius_table.py @@ -0,0 +1,132 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Covalent radius table +# +# Source: https://github.com/openbabel/openbabel/blob/master/src/elementtable.h +# - covalent radii (in Angstrom) 1.6 if unknown +# from https://doi.org/10.1039/b801115j +# + + +covalent_radius = { + "H": 0.31, + "D": 0.00, + "He": 0.28, + "Li": 1.28, + "Be": 0.96, + "B": 0.84, + "C": 0.76, + "N": 0.71, + "O": 0.66, + "F": 0.57, + "Ne": 0.58, + "Na": 1.66, + "Mg": 1.41, + "Al": 1.21, + "Si": 1.11, + "P": 1.07, + "S": 1.05, + "Cl": 1.02, + "Ar": 1.06, + "K ": 2.03, + "Ca": 1.76, + "Sc": 1.70, + "Ti": 1.60, + "V": 1.53, + "Cr": 1.39, + "Mn": 1.39, + "Fe": 1.32, + "Co": 1.26, + "Ni": 1.24, + "Cu": 1.32, + "Zn": 1.22, + "Ga": 1.22, + "Ge": 1.20, + "As": 1.19, + "Se": 1.20, + "Br": 1.20, + "Kr": 1.16, + "Rb": 2.20, + "Sr": 1.95, + "Y": 1.90, + "Zr": 1.75, + "Nb": 1.64, + "Mo": 1.54, + "Tc": 1.47, + "Ru": 1.46, + "Rh": 1.42, + "Pd": 1.39, + "Ag": 1.45, + "Cd": 1.44, + "In": 1.42, + "Sn": 1.39, + "Sb": 1.39, + "Te": 1.38, + "I": 1.39, + "Xe": 1.40, + "Cs": 2.44, + "Ba": 2.15, + "La": 2.07, + "Ce": 2.04, + "Pr": 2.03, + "Nd": 2.01, + "Pm": 1.99, + "Sm": 1.98, + "Eu": 1.98, + "Gd": 1.96, + "Tb": 1.94, + "Dy": 1.92, + "Ho": 1.92, + "Er": 1.89, + "Tm": 1.90, + "Yb": 1.87, + "Lu": 1.87, + "Hf": 1.75, + "Ta": 1.70, + "W": 1.62, + "Re": 1.51, + "Os": 1.44, + "Ir": 1.41, + "Pt": 1.36, + "Au": 1.36, + "Hg": 1.32, + "Tl": 1.45, + "Pb": 1.46, + "Bi": 1.48, + "Po": 1.40, + "At": 1.50, + "Rn": 1.50, + "Fr": 2.60, + "Ra": 2.21, + "Ac": 2.15, + "Th": 2.06, + "Pa": 2.00, + "U": 1.96, + "Np": 1.90, + "Pu": 1.87, + "Am": 1.80, + "Cm": 1.69, + "Bk": 1.60, + "Cf": 1.60, + "Es": 1.60, + "Fm": 1.60, + "Md": 1.60, + "No": 1.60, + "Lr": 1.60, + "Rf": 1.60, + "Db": 1.60, + "Sg": 1.60, + "Bh": 1.60, + "Hs": 1.60, + "Mt": 1.60, + "Ds": 1.60, + "Rg": 1.60, + "Cn": 1.60, + "Nh": 1.60, + "Fl": 1.60, + "Mc": 1.60, + "Lv": 1.60, + "Ts": 1.60, + "Og": 1.60 +} diff --git a/src/moldrug/meeko_light/utils/jsonutils.py b/src/moldrug/meeko_light/utils/jsonutils.py new file mode 100644 index 0000000..0f2a344 --- /dev/null +++ b/src/moldrug/meeko_light/utils/jsonutils.py @@ -0,0 +1,76 @@ +from rdkit import Chem +from rdkit.Chem import rdMolInterchange + + +SERIALIZATION_SEPARATOR_CHAR = "," + + +def rdkit_mol_from_json(json_str: str): + """ + Takes in a JSON string and attempts to use RDKit's JSON to Mols utility to extract just one RDKitMol from the + json string. If none or more than one Mols are returned, raises an error. + + Parameters + ---------- + json_str: str + A JSON string representing an RDKit Mol. + + Returns + ------- + rdkit_mol: rdkit.Chem.rdchem.Mol + An RDKit Mol object corresponding to the input JSON string + + Raises + ------ + ValueError + If no RDKitMol objects are returned, or if more than one is returned, throws a ValueError. + """ + if json_str is None: + return None + rdkit_mols = rdMolInterchange.JSONToMols(json_str) + if len(rdkit_mols) != 1: + raise ValueError( + f"Expected 1 rdkit mol from json string but got {len(rdkit_mols)}" + ) + Chem.SanitizeMol(rdkit_mols[0]) # needed to compute gasteiger charges + return rdkit_mols[0] + + +def tuple_to_string(input_tuple: tuple): + """ + Converts a tuple to a JSON serializable string. + + Parameters + ---------- + input_tuple: tuple + A tuple to convert to a JSON serializable string. + + Returns + ------- + A string representation of the tuple using the specified serialization separator character. + """ + return SERIALIZATION_SEPARATOR_CHAR.join([str(i) for i in input_tuple]) + + +def string_to_tuple(input_string: str, element_type: type = str): + """ + Takes a JSON string and converts it back to a tuple. If element type is specified, converts all elements of the + tuple to that type. + + Parameters + ---------- + input_string: str + String deserialized from JSON. + element_type: type + Data type for all of the elements of the tuple. + + Returns + ------- + A deserialized tuple with the specified element type. + """ + if element_type is not str: + return tuple( + [element_type(i) for i in input_string.split(SERIALIZATION_SEPARATOR_CHAR)] + ) + else: + return tuple(input_string) diff --git a/src/moldrug/meeko_light/utils/pdbutils.py b/src/moldrug/meeko_light/utils/pdbutils.py new file mode 100644 index 0000000..413ca99 --- /dev/null +++ b/src/moldrug/meeko_light/utils/pdbutils.py @@ -0,0 +1,5 @@ +from collections import namedtuple + +# named tuple to contain information about an atom +PDBAtomInfo = namedtuple("PDBAtomInfo", "name resName resNum icode chain") +PDBResInfo = namedtuple("PDBResInfo", "resName resNum chain") # used in obutils, maybe diff --git a/src/moldrug/meeko_light/utils/rdkitutils.py b/src/moldrug/meeko_light/utils/rdkitutils.py new file mode 100644 index 0000000..b16c43b --- /dev/null +++ b/src/moldrug/meeko_light/utils/rdkitutils.py @@ -0,0 +1,278 @@ +from rdkit import Chem +from rdkit.Chem import rdChemReactions +from .utils import mini_periodic_table +from .pdbutils import PDBAtomInfo +from rdkit.Geometry import Point3D +from rdkit.Chem import rdDetermineBonds +from typing import Tuple + +periodic_table = Chem.GetPeriodicTable() + + +""" +create new RDKIT residue + +mi = Chem.AtomPDBResidueInfo() +mi.SetResidueName('MOL') +mi.SetResidueNumber(1) +mi.SetOccupancy(0.0) +mi.SetTempFactor(0.0) + +source: https://sourceforge.net/p/rdkit/mailman/message/36404394/ +""" + + +def getPdbInfoNoNull(atom): + """extract information for populating an ATOM/HETATM line + in the PDB""" + minfo = atom.GetMonomerInfo() # same as GetPDBResidueInfo + if minfo is None: + atomic_number = atom.GetAtomicNum() + if atomic_number == 0: + name = "%-2s" % "*" + else: + name = "%-2s" % mini_periodic_table[atomic_number] + chain = " " + resNum = 1 + icode = "" + resName = "UNL" + else: + name = minfo.GetName() + chain = minfo.GetChainId() + resNum = minfo.GetResidueNumber() + icode = minfo.GetInsertionCode() + resName = minfo.GetResidueName() + return PDBAtomInfo( + name=name, resName=resName, resNum=resNum, icode=icode, chain=chain + ) + + +class Mol2MolSupplier: + """RDKit Mol2 molecule supplier. + Parameters + sanitize: perform RDKit sanitization of Mol2 molecule""" + + def __init__( + self, filename, sanitize=True, removeHs=False, cleanupSubstructures=True + ): + self.fp = open(filename, "r") + self._opts = { + "sanitize": sanitize, + "removeHs": removeHs, + "cleanupSubstructures": cleanupSubstructures, + } + self.buff = [] + + def __iter__(self): + return self + + def __next__(self): + """iterator step""" + while True: + line = self.fp.readline() + # empty line + if not line: + if len(self.buff): + # buffer full, returning last molecule + mol = Chem.MolFromMol2Block("".join(self.buff), **self._opts) + self.buff = [] + return mol + # buffer empty, stopping the iteration + self.fp.close() + raise StopIteration + if "@MOLECULE" in line: + # first molecule parsed + if len(self.buff) == 0: + self.buff.append(line) + else: + # found the next molecule, breaking to return the complete one + break + else: + # adding another line in the current molecule + self.buff.append(line) + # found a complete molecule, returning it + mol = Chem.MolFromMol2Block("".join(self.buff), **self._opts) + self.buff = [line] + return mol + + +class AtomField: + """Stores data parsed from PDB or mmCIF""" + + def __init__( + self, + atomname: str, + altloc: str, + resname: str, + chain: str, + resnum: int, + icode: str, + x: float, + y: float, + z: float, + element: str, + ): + self.atomname = atomname + self.altloc = altloc + self.resname = resname + self.chain = chain + self.resnum = resnum + self.icode = icode + self.x = x + self.y = y + self.z = z + if len(element) > 1: + element = f"{element[0].upper()}{element[1].lower()}" + else: + element = f"{element.upper()}" + self.atomic_nr = periodic_table.GetAtomicNumber(element) + + +def _build_rdkit_mol_for_altloc(atom_fields_list, wanted_altloc: str = None): + mol = Chem.EditableMol(Chem.Mol()) + mol.BeginBatchEdit() + positions = [] + idx_to_rdkit = {} + for index_list, atom in enumerate(atom_fields_list): + if wanted_altloc is not None: + if atom.altloc and atom.altloc != wanted_altloc: + # if atom.altloc is "" we still want to consider this atom + continue + rdkit_atom = Chem.Atom(atom.atomic_nr) + positions.append(Point3D(atom.x, atom.y, atom.z)) + res_info = Chem.AtomPDBResidueInfo() + res_info.SetName(atom.atomname) + res_info.SetResidueName(atom.resname) + res_info.SetResidueNumber(atom.resnum) + res_info.SetChainId(atom.chain) + res_info.SetInsertionCode(atom.icode) + rdkit_atom.SetPDBResidueInfo(res_info) + index_rdkit = mol.AddAtom(rdkit_atom) + idx_to_rdkit[index_list] = index_rdkit + mol.CommitBatchEdit() + mol = mol.GetMol() + conformer = Chem.Conformer(mol.GetNumAtoms()) + for index, position in enumerate(positions): + conformer.SetAtomPosition(index, position) + mol.AddConformer(conformer, assignId=True) + return mol, idx_to_rdkit + + +def build_one_rdkit_mol_per_altloc(atom_fields_list): + """ if no altlocs, the only key in the output dict is None + if altlocs exist, None is not a key: the keys are the altloc IDs + """ + altlocs = set([atom.altloc for atom in atom_fields_list if atom.altloc]) + rdkit_mol_dict = {} + if not altlocs: + altlocs = {None} + for altloc in altlocs: + mol, idx_to_rdkit = _build_rdkit_mol_for_altloc(atom_fields_list, altloc) + rdkit_mol_dict[altloc] = (mol, idx_to_rdkit) + return rdkit_mol_dict + + +def _aux_altloc_mol_build(atom_field_list, requested_altloc, default_altloc): + missed_altloc = False + needed_altloc = False + mols_dict = build_one_rdkit_mol_per_altloc(atom_field_list) + has_altloc = None not in mols_dict + if has_altloc and requested_altloc is None and default_altloc is None: + pdbmol = None + missed_altloc = False + needed_altloc = True + elif requested_altloc and requested_altloc in mols_dict: + pdbmol, idx_to_rdkit = mols_dict[requested_altloc] + elif requested_altloc and requested_altloc not in mols_dict: + pdbmol = None + missed_altloc = True + needed_altloc = False + elif default_altloc and default_altloc in mols_dict: + pdbmol, idx_to_rdkit = mols_dict[default_altloc] + elif has_altloc and default_altloc not in mols_dict: + pdbmol = None + missed_altloc = True + needed_altloc = False + elif not has_altloc and requested_altloc is None: + pdbmol, idx_to_rdkit = mols_dict[None] + else: + raise RuntimeError("programming bug, please post full error on github") + if pdbmol is None: + idx_to_rdkit = None + return pdbmol, idx_to_rdkit, missed_altloc, needed_altloc + else: + rdDetermineBonds.DetermineConnectivity(pdbmol) + for atom in pdbmol.GetAtoms(): + if atom.GetAtomicNum() == 7 and len(atom.GetNeighbors()) == 4: + atom.SetFormalCharge(1) + _ = Chem.SanitizeMol(pdbmol) + + return pdbmol, idx_to_rdkit, missed_altloc, needed_altloc + + +def react_and_map(reactants: Tuple[Chem.Mol], rxn: rdChemReactions.ChemicalReaction): + """ + Run a reaction and keep track of atom indices from reactants to products. + + Parameters + ---------- + reactants : tuple[Chem.Mol] + A tuple of RDKit molecule objects representing the reactants. + rxn : rdChemReactions.ChemicalReaction + The RDKit reaction object. + + Returns + ------- + list[tuple[Chem.Mol, dict[str, list[Optional[int]]]]] + A list of tuples where each tuple contains a product molecule and a dictionary. + The dictionary has keys 'atom_idx' and 'new_atom_label', which are ordered lists for product atoms: + - 'atom_idx' holds the corresponding atom indices in reactant. None for newly added atoms. + - 'new_atom_label' holds the reaction mapping number, only for newly added atoms. + """ + + # Prepare for multiple possible outcomes resulted from multiple matched reactive sites in reactant + outcomes = [] + for products in rxn.RunReactants(reactants): + # Assumes single product + product = products[0] + # For each atom, get react_atom_idx if they were in reactant + atom_idxmap = [ + atom.GetIntProp("react_atom_idx") if atom.HasProp("react_atom_idx") + else None + for atom in product.GetAtoms() + ] + # For each atom, get the rxn mapping number if the were added in the rxn + new_atom_label = [ + atom.GetIntProp("old_mapno") if atom.HasProp("old_mapno") and not atom.HasProp("react_atom_idx") + else None + for atom in product.GetAtoms() + ] + # Collect product and index_map + index_map = {"atom_idx": atom_idxmap, "new_atom_label": new_atom_label} + outcomes.append((product, index_map)) + + return outcomes + + +covalent_radius = { # from wikipedia + 1: 0.31, + 5: 0.84, + 6: 0.76, + 7: 0.71, + 8: 0.66, + 9: 0.57, + 12: 0.00, # hack to avoid bonds with metals + 14: 1.11, + 15: 1.07, + 16: 1.05, + 17: 1.02, + # 19: 2.03, + 20: 0.00, + # 24: 1.39, + 25: 0.00, # hack to avoid bonds with metals + 26: 0.00, + 30: 0.00, # hack to avoid bonds with metals + # 34: 1.20, + 35: 1.20, + 53: 1.39, +} diff --git a/src/moldrug/meeko_light/utils/utils.py b/src/moldrug/meeko_light/utils/utils.py new file mode 100644 index 0000000..1388021 --- /dev/null +++ b/src/moldrug/meeko_light/utils/utils.py @@ -0,0 +1,425 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Raccoon +# collection of useful snippets of code that's used frequently +# + +import os +from operator import itemgetter + + +mini_periodic_table = { + 1: "H", + 2: "He", + 3: "Li", + 4: "Be", + 5: "B", + 6: "C", + 7: "N", + 8: "O", + 9: "F", + 10: "Ne", + 11: "Na", + 12: "Mg", + 13: "Al", + 14: "Si", + 15: "P", + 16: "S", + 17: "Cl", + 18: "Ar", + 19: "K", + 20: "Ca", + 21: "Sc", + 22: "Ti", + 23: "V", + 24: "Cr", + 25: "Mn", + 26: "Fe", + 27: "Co", + 28: "Ni", + 29: "Cu", + 30: "Zn", + 31: "Ga", + 32: "Ge", + 33: "As", + 34: "Se", + 35: "Br", + 36: "Kr", + 37: "Rb", + 38: "Sr", + 39: "Y", + 40: "Zr", + 41: "Nb", + 42: "Mo", + 43: "Tc", + 44: "Ru", + 45: "Rh", + 46: "Pd", + 47: "Ag", + 48: "Cd", + 49: "In", + 50: "Sn", + 51: "Sb", + 52: "Te", + 53: "I", + 54: "Xe", + 55: "Cs", + 56: "Ba", + 57: "La", + 58: "Ce", + 59: "Pr", + 60: "Nd", + 61: "Pm", + 62: "Sm", + 63: "Eu", + 64: "Gd", + 65: "Tb", + 66: "Dy", + 67: "Ho", + 68: "Er", + 69: "Tm", + 70: "Yb", + 71: "Lu", + 72: "Hf", + 73: "Ta", + 74: "W", + 75: "Re", + 76: "Os", + 77: "Ir", + 78: "Pt", + 79: "Au", + 80: "Hg", + 81: "Tl", + 82: "Pb", + 83: "Bi", + 84: "Po", + 85: "At", + 86: "Rn", + 87: "Fr", + 88: "Ra", +} + + +def getNameExt(fname): + """extract name and extension from the input file, removing the dot + filename.ext -> [filename, ext] + """ + name, ext = os.path.splitext(fname) + return name, ext[1:] # .lower() + + +class HJKRingDetection: + """Implementation of the Hanser-Jauffret-Kaufmann exhaustive ring detection + algorithm: + ref: + Th. Hanser, Ph. Jauffret, and G. Kaufmann + J. Chem. Inf. Comput. Sci. 1996, 36, 1146-1152 + """ + + def __init__(self, mgraph, max_iterations=8000000): + self.mgraph = {key: [x for x in values] for (key, values) in mgraph.items()} + self.rings = [] + self._iterations = 0 + self._max_iterations = max_iterations + self._is_failed = False + + def scan(self, keep_chorded_rings=False, keep_equivalent_rings=False): + """run the full protocol for exhaustive ring detection + by default, only chordless rings are kept, and equivalent rings removed. + (equivalent rings are rings that have the same size and share the same + neighbors) + """ + self.prune() + self.build_pgraph() + self.vertices = self._get_sorted_vertices() + while self.vertices: + self._remove_vertex(self.vertices[0]) + if not keep_chorded_rings: + self.find_chordless_rings(keep_equivalent_rings) + output_rings = [] + for ring in self.rings: + output_rings.append(tuple(ring[:-1])) + return output_rings + + def _get_sorted_vertices(self): + """function to return the vertices to be removed, sorted by increasing + connectivity order (see paper)""" + vertices = ((k, len(v)) for k, v in self.mgraph.items()) + return [x[0] for x in sorted(vertices, key=itemgetter(1))] + + def prune(self): + """iteratively prune graph until there are no leafs left (nodes with only + one connection)""" + while True: + prune = [] + for node, neighbors in self.mgraph.items(): + if len(neighbors) == 1: + prune.append((node, neighbors)) + if len(prune) == 0: + break + for node, neighbors in prune: + self.mgraph.pop(node) + for n in neighbors: + self.mgraph[n].remove(node) + + def build_pgraph(self, prune=True): + """convert the M-graph (molecular graph) into the P-graph (path/bond graph)""" + self.pgraph = [] + for node, neigh in self.mgraph.items(): + for n in neigh: + # use sets for unique id + edge = set((node, n)) + if edge not in self.pgraph: + self.pgraph.append(edge) + # re-convert the edges to lists because order matters in cycle detection + self.pgraph = [list(x) for x in self.pgraph] + + def _remove_vertex(self, vertex): + """remove a vertex and join all edges connected by that vertex (this is + the REMOVE function from the paper) + """ + visited = {} + pool = [] + for path in self.pgraph: + if self._has_vertex(vertex, path): + pool.append(path) + for i, path1 in enumerate(pool): + for j, path2 in enumerate(pool): + if i == j: + continue + self._iterations += 1 + if self._iterations > self._max_iterations: + self._is_failed = True + break + pair_id = tuple(set((i, j))) + if pair_id in visited: + continue + visited[pair_id] = None + common = list(set(path1) & set(path2)) + common_count = len(common) + # check if two paths have only this vertex in common or (or + # two, if they're a cycle) + if not 1 <= common_count <= 2: + continue + # generate the joint path + joint_path = self._concatenate_path(path1, path2, vertex) + is_ring = joint_path[0] == joint_path[-1] + # if paths share more than two vertices but they're not a ring, then skip + if (common_count == 2) and not is_ring: + continue + # store the ring... + if is_ring: + self._add_ring(joint_path) + # ...or the common path + elif joint_path not in self.pgraph: + self.pgraph.append(joint_path) + # remove used paths + for p in pool: + self.pgraph.remove(p) + # remove the used vertex + self.vertices.remove(vertex) + + def _add_ring(self, ring): + """add newly found rings to the list (if not already there)""" + r = set(ring) + for candidate in self.rings: + if r == set(candidate): + return + self.rings.append(ring) + + def _has_vertex(self, vertex, edge): + """check if the vertex is part of this edge, and if true, return the + sorted edge so that the vertex is the first in the list""" + if edge[0] == vertex: + return edge + if edge[-1] == vertex: + return edge[::-1] + return None + + def _concatenate_path(self, path1, path2, v): + """concatenate two paths sharing a common vertex + a-b, c-b => a-b-c : idx1=1, idx2=1 + b-a, c-b => a-b-c : idx1=0, idx2=1 + a-b, b-c => a-b-c : idx1=1, idx2=0 + b-a, b-c => a-b-c : idx1=0, idx2=0 + """ + if not path1[-1] == v: + path1.reverse() + if not path2[0] == v: + path2.reverse() + return path1 + path2[1:] + + def _edge_in_pgraph(self, edge): + """check if edge is already in pgraph""" + e = set(edge) + for p in self.pgraph: + if e == set(p) and len(p) == len(edge): + return True + return False + + def find_chordless_rings(self, keep_equivalent_rings): + """find chordless rings: cycles in which two vertices are not connected + by an edge that does not itself belong to the cycle (Source: + https://en.wikipedia.org/wiki/Cycle_%28graph_theory%29#Chordless_cycle) + + - iterate through rings starting from the smallest ones: A,B,C,D... + - for each ring (A), find a candidate (e.g.: B) that is smaller and shares at least an edge + - for this pair, calculate the two differences (A-B and B-A) in the list of edges of each + - if ( (A-B) + (B-A) ) a smaller ring (e.g.: C), then the current ring has a chord + """ + # sort rings by the smallest to largest + self.rings.sort(key=len, reverse=False) + chordless_rings = [] + ring_edges = [] + rings_set = [set(x) for x in self.rings] + for r in self.rings: + edges = [] + for i in range(len(r) - 1): + edges.append( + tuple( + set((r[i], r[(i + 1) % len(r)])), + ) + ) + edges = sorted(edges, key=itemgetter(0)) + ring_edges.append(edges) + ring_contacts = {} + for i, r1 in enumerate(self.rings): + chordless = True + r1_edges = ring_edges[i] + ring_contacts[i] = [] + for j, r2 in enumerate(self.rings): + if i == j: + continue + if len(r2) >= len(r1): + # the candidate ring is larger than or the same size of the candidate + continue + # avoid rings that don't share at least an edge + # shared = set(r1) & set(r2) + r2_edges = ring_edges[j] + shared = set(r1_edges) & set(r2_edges) + if len(shared) < 1: + continue + ring_contacts[i].append(j) + # get edges difference (r2_edges - r1_edges) + core_edges = [x for x in r2_edges if x not in r1_edges] + chord = [x for x in r1_edges if x not in r2_edges] + # combined = chord + core_edges + ring_new = [] + for edge in chord + core_edges: + ring_new.append(edge[0]) + ring_new.append(edge[1]) + ring_new = set(ring_new) + if (ring_new in rings_set) and (len(ring_new) < len(r1) - 1): + chordless = False + break + if chordless: + chordless_rings.append(i) + ring_contacts[i] = set(ring_contacts[i]) + if not keep_equivalent_rings: + chordless_rings = self._remove_equivalent_rings( + chordless_rings, ring_contacts + ) + self.rings = [self.rings[x] for x in chordless_rings] + return + + def _remove_equivalent_rings(self, chordless_rings, ring_contacts): + """remove equivalent rings by clustering by size, then by ring neighbors. + Two rings A and B are equivalent if satisfy the following conditions: + - same size + - same neighbor ring(s) [C,D, ...] + - (A - C) == (B -C) + """ + size_clusters = {} + # cluster rings by their size + for ring_id in chordless_rings: + if len(ring_contacts[ring_id]) == 0: + continue + size = len(self.rings[ring_id]) - 1 + if size not in size_clusters: + size_clusters[size] = [] + size_clusters[size].append(ring_id) + remove = [] + # process rings of the same size + for size, ring_pool in size_clusters.items(): + for ri in ring_pool: + if ri in remove: + continue + for rj in ring_pool: + if ri == rj: + continue + common_neigh = ring_contacts[ri] & ring_contacts[rj] + for c in common_neigh: + d1 = set(self.rings[ri]) - set(self.rings[c]) + d2 = set(self.rings[rj]) - set(self.rings[c]) + if d1 == d2: + remove.append(rj) + chordless_rings = [i for i in chordless_rings if i not in set(remove)] + # for r in set(remove): + # chordless_rings.remove(r) + return chordless_rings + + +def parse_begin_res(string): + """ + "THR A 315" -> "A:315" + "THR 5" -> ":5" + " A 5" -> "A:5" + "THR A 315B" -> "A:5B" + + meeko v0.5-v0.6.0a3 formatting: + - "%s %s %s%s" % (resname, chain, resnum, icode) + - input strings have been stripped, thus, we got + exactly two whitespace in the input string. + + older (probably mgltools) format include more + whitespace: "SER A 56". + + This tries to be flexible and parse even: + "SER A1234C" + " A 123" + "SER 123C + """ + + if len(string.split()) > 3 or len(string) > 12 or len(string) < 3: + raise ValueError(f"{string} doesn't look like a BEGIN_RES flexres") + + resnum = "" + icode = "" + resname_chain = "" + got_space = False + for char in reversed(string): + + if char.isspace(): + got_space = True + + if char.isdigit() and not resname_chain and not got_space: + resnum = char + resnum + + elif not char.isdigit() and not resnum and not icode and not got_space: + icode = char + + # prevent icode from having more than 1 char + elif not char.isdigit() and not resnum and icode: + raise ValueError(f"{string=} misses resnum or has len(icode) > 1") + + # if we are here, resnum already exists, so it's resname_chain land + elif not char.isdigit() and not resname_chain: + resname_chain = char + + else: + resname_chain = char + resname_chain + + fields = resname_chain.split() + if len(fields) == 1 and len(fields[0]) <= 2: + chain = fields[0] + elif len(fields) == 1 and len(fields[0]) > 2: + chain = "" # assume fields[0] is resname and no chain + elif len(fields) == 2 and len(fields[1]) <= 2: + chain = fields[1] + elif len(fields) == 0 and resnum: + chain = "" + else: + raise ValueError(f"can't parse {string=}") + + return f"{chain}:{resnum}{icode}" diff --git a/src/moldrug/meeko_light/utils/van_der_waals_radius_table.py b/src/moldrug/meeko_light/utils/van_der_waals_radius_table.py new file mode 100644 index 0000000..f185c71 --- /dev/null +++ b/src/moldrug/meeko_light/utils/van_der_waals_radius_table.py @@ -0,0 +1,132 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# van der Waals table +# +# Source: https://github.com/openbabel/openbabel/blob/master/src/elementtable.h +# - van der Waals radii (in Angstrom) 2.0 if unknown +# from https://doi.org/10.1021/jp8111556 +# + + +van_der_waals_radius = { + "H": 1.10, + "D": 0.00, + "He": 1.40, + "Li": 1.81, + "Be": 1.53, + "B": 1.92, + "C": 1.70, + "N": 1.55, + "O": 1.52, + "F": 1.47, + "Ne": 1.54, + "Na": 2.27, + "Mg": 1.73, + "Al": 1.84, + "Si": 2.10, + "P": 1.80, + "S": 1.80, + "Cl": 1.75, + "Ar": 1.88, + "K": 2.75, + "Ca": 2.31, + "Sc": 2.30, + "Ti": 2.15, + "V": 2.05, + "Cr": 2.05, + "Mn": 2.05, + "Fe": 2.05, + "Co": 2.00, + "Ni": 2.00, + "Cu": 2.00, + "Zn": 2.10, + "Ga": 1.87, + "Ge": 2.11, + "As": 1.85, + "Se": 1.90, + "Br": 1.83, + "Kr": 2.02, + "Rb": 3.03, + "Sr": 2.49, + "Y": 2.40, + "Zr": 2.30, + "Nb": 2.15, + "Mo": 2.10, + "Tc": 2.05, + "Ru": 2.05, + "Rh": 2.00, + "Pd": 2.05, + "Ag": 2.10, + "Cd": 2.20, + "In": 2.20, + "Sn": 1.93, + "Sb": 2.17, + "Te": 2.06, + "I": 1.98, + "Xe": 2.16, + "Cs": 3.43, + "Ba": 2.68, + "La": 2.50, + "Ce": 2.48, + "Pr": 2.47, + "Nd": 2.45, + "Pm": 2.43, + "Sm": 2.42, + "Eu": 2.40, + "Gd": 2.38, + "Tb": 2.37, + "Dy": 2.35, + "Ho": 2.33, + "Er": 2.32, + "Tm": 2.30, + "Yb": 2.28, + "Lu": 2.27, + "Hf": 2.25, + "Ta": 2.20, + "W": 2.10, + "Re": 2.05, + "Os": 2.00, + "Ir": 2.00, + "Pt": 2.05, + "Au": 2.10, + "Hg": 2.05, + "Tl": 1.96, + "Pb": 2.02, + "Bi": 2.07, + "Po": 1.97, + "At": 2.02, + "Rn": 2.20, + "Fr": 3.48, + "Ra": 2.83, + "Ac": 2.00, + "Th": 2.40, + "Pa": 2.00, + "U": 2.30, + "Np": 2.00, + "Pu": 2.00, + "Am": 2.00, + "Cm": 2.00, + "Bk": 2.00, + "Cf": 2.00, + "Es": 2.00, + "Fm": 2.00, + "Md": 2.00, + "No": 2.00, + "Lr": 2.00, + "Rf": 2.00, + "Db": 2.00, + "Sg": 2.00, + "Bh": 2.00, + "Hs": 2.00, + "Mt": 2.00, + "Ds": 2.00, + "Rg": 2.00, + "Cn": 2.00, + "Nh": 2.00, + "Fl": 2.00, + "Mc": 2.00, + "Lv": 2.00, + "Ts": 2.00, + "Og": 2.00 + } diff --git a/src/moldrug/meeko_light/writer.py b/src/moldrug/meeko_light/writer.py new file mode 100644 index 0000000..3107ff9 --- /dev/null +++ b/src/moldrug/meeko_light/writer.py @@ -0,0 +1,766 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Meeko PDBQT writer +# + +import json +import math +from os import linesep as eol + +import numpy as np + +from .molsetup import Bond +from .utils import pdbutils + + +def oids_json_from_setup(molsetup, name="LigandFromMeeko"): + if len(molsetup.restraints): + raise NotImplementedError( + "molsetup has restraints but these aren't written to oids block yet" + ) + offchrg_type = "OFFCHRG" + offchrg_by_parent = {} + for atom in molsetup.atoms: + if atom.is_pseudo_atom and atom.atom_type == offchrg_type: + neigh = molsetup.get_neighbors(atom.index) + if len(neigh) != 1: + raise RuntimeError( + "offsite charge %s is bonded to: %s which has len() != 1" + % (atom.index, json.dumps(neigh)) + ) + if neigh[0] in offchrg_by_parent: + raise RuntimeError( + "atom %d has more than one offsite charge" % neigh[0] + ) + offchrg_by_parent[neigh[0]] = atom.index + output_indices_start_at_one = True + index_start = int(output_indices_start_at_one) + positions_block = "" + charges = [] + offchrg_by_oid_parent = {} + elements = [] + indexmap = {} # molsetup: oid + count_oids = 0 + for atom in molsetup.atoms: + index = atom.index + if atom.is_dummy or atom.is_pseudo_atom or atom.is_ignore: + continue + if atom.atom_type == offchrg_type: + continue # handled by offchrg_by_parent + oid_id = count_oids + index_start + indexmap[index] = count_oids + x, y, z = atom.coord + positions_block += "position.%d = (%f,%f,%f)\n" % (oid_id, x, y, z) + charges.append(atom.charge) + if index in offchrg_by_parent: + index_pseudo = offchrg_by_parent[index] + xq_abs, yq_abs, zq_abs = molsetup.atoms[index_pseudo].coord + xq_rel = xq_abs - x + yq_rel = yq_abs - y + zq_rel = zq_abs - z + offchrg_by_oid_parent[count_oids] = { + "q": molsetup.atoms[index_pseudo].charge, + "xyz": (xq_rel, yq_rel, zq_rel), + } + count_oids += 1 + element = "%s %s %d" % (name, atom.atom_type, oid_id) + elements.append(element) + + tmp = [] + for index in enumerate(charges): + if index in offchrg_by_oid_parent: + tmplist = [ + "%f" % charges[index], + "0.0", + "0.0", + "0.0", + ] # xyz relative to current elemtn + tmplist.append("%f" % offchrg_by_oid_parent[index]["q"]) + tmplist.append("%f,%f,%f" % offchrg_by_oid_parent[index]["xyz"]) + tmp.append(",".join(tmplist)) + else: + tmp.append("%f" % charges[index]) + charges_line = "import_charges = {%s}\n" % ("|".join(tmp)) + elements_line = "elements = %s\n" % (",".join(elements)) + + bonds = [[] for _ in range(count_oids)] + # bond_orders = [[] for _ in range(count_oids)] + static_links = [] + for i, j in molsetup.bond_info.keys(): + if molsetup.get_is_ignore(i) or molsetup.get_is_ignore(j): + continue + if ( + molsetup.get_atom_type(i) == offchrg_type + or molsetup.get_atom_type(j) == offchrg_type + ): + continue + oid_i = indexmap[i] + oid_j = indexmap[j] + bonds[oid_i].append("%d" % (oid_j + index_start)) + # bond_orders[oid_i].append("%d" % molsetup.bond_info[(i, j)].order) + if not molsetup.bond_info[(i, j)].rotatable: + static_links.append("%d,%d" % (oid_i + index_start, oid_j + index_start)) + bonds = [",".join(j_list) for j_list in bonds] + bonds_line = "connectivity = {%s}\n" % ("|".join(bonds)) + # bond_orders = [",".join(orders) for orders in bond_orders] + # bondorder_line = "bond_order = {%s}\n" % ("|".join(bond_orders)) + staticlinks_line = "static_links = {%s}\n" % ("|".join(static_links)) + + output = "" + output += "[Group: %s]\n" % name + output += positions_block + output += charges_line + output += elements_line + output += bonds_line + # output += bondorder_line + output += staticlinks_line + output += "number = 1\t\t// can only be 1 for the sandbox currently (but any number for classical MC)\n" + output += "group_dipole = 1\t// not relevant for sandbox but classical MC\n" + output += "rand_independent=0\t// not relevant for sandbox but classical MC\n" + output += "bond_range = 4\t\t// bond range AD default\n" + output += "\n" + output += get_dihedrals_block(molsetup, indexmap, name) + + return output, indexmap + + +def oids_block_from_setup(molsetup, name="LigandFromMeeko"): + if len(molsetup.restraints): + raise NotImplementedError( + "molsetup has restraints but these aren't written to oids block yet" + ) + offchrg_type = "OFFCHRG" + offchrg_by_parent = {} + for atom in molsetup.atoms: + if atom.is_pseudo_atom and atom.atom_type == offchrg_type: + neigh = molsetup.get_neighbors(atom.index) + if len(neigh) != 1: + raise RuntimeError( + "offsite charge %s is bonded to: %s which has len() != 1" + % (atom.index, json.dumps(neigh)) + ) + if neigh[0] in offchrg_by_parent: + raise RuntimeError( + "atom %d has more than one offsite charge" % neigh[0] + ) + offchrg_by_parent[neigh[0]] = atom.index + output_indices_start_at_one = True + index_start = int(output_indices_start_at_one) + positions_block = "" + charges = [] + offchrg_by_oid_parent = {} + elements = [] + indexmap = {} # molsetup: oid + count_oids = 0 + for atom in molsetup.atoms: + index = atom.index + if atom.is_dummy or atom.is_pseudo_atom or atom.is_ignore: + continue + if atom.atom_type == offchrg_type: + continue # handled by offchrg_by_parent + oid_id = count_oids + index_start + indexmap[index] = count_oids + x, y, z = atom.coord + positions_block += "position.%d = (%f,%f,%f)\n" % (oid_id, x, y, z) + charges.append(atom.charge) + if index in offchrg_by_parent: + index_pseudo = offchrg_by_parent[index] + xq_abs, yq_abs, zq_abs = molsetup.atoms[index_pseudo].coord + xq_rel = xq_abs - x + yq_rel = yq_abs - y + zq_rel = zq_abs - z + offchrg_by_oid_parent[count_oids] = { + "q": molsetup.get_charge(index_pseudo), + "xyz": (xq_rel, yq_rel, zq_rel), + } + count_oids += 1 + element = "%s %s %d" % (name, molsetup.get_atom_type(index), oid_id) + elements.append(element) + + tmp = [] + for index in range(len(charges)): + if index in offchrg_by_oid_parent: + tmplist = [ + "%f" % charges[index], + "0.0", + "0.0", + "0.0", + ] # xyz relative to current elemtn + tmplist.append("%f" % offchrg_by_oid_parent[index]["q"]) + tmplist.append("%f,%f,%f" % offchrg_by_oid_parent[index]["xyz"]) + tmp.append(",".join(tmplist)) + else: + tmp.append("%f" % charges[index]) + charges_line = "import_charges = {%s}\n" % ("|".join(tmp)) + elements_line = "elements = %s\n" % (",".join(elements)) + + bonds = [[] for _ in range(count_oids)] + # bond_orders = [[] for _ in range(count_oids)] + static_links = [] + for i, j in molsetup.bond_info.keys(): + if molsetup.get_is_ignore(i) or molsetup.get_is_ignore(j): + continue + if ( + molsetup.get_atom_type(i) == offchrg_type + or molsetup.get_atom_type(j) == offchrg_type + ): + continue + oid_i = indexmap[i] + oid_j = indexmap[j] + bonds[oid_i].append("%d" % (oid_j + index_start)) + # bond_orders[oid_i].append("%d" % molsetup.bond_info[(i, j)].order) + if not molsetup.bond_info[(i, j)].rotatable: + static_links.append("%d,%d" % (oid_i + index_start, oid_j + index_start)) + bonds = [",".join(j_list) for j_list in bonds] + bonds_line = "connectivity = {%s}\n" % ("|".join(bonds)) + # bond_orders = [",".join(orders) for orders in bond_orders] + # bondorder_line = "bond_order = {%s}\n" % ("|".join(bond_orders)) + staticlinks_line = "static_links = {%s}\n" % ("|".join(static_links)) + + output = "" + output += "[Group: %s]\n" % name + output += positions_block + output += charges_line + output += elements_line + output += bonds_line + # output += bondorder_line + output += staticlinks_line + output += "number = 1\t\t// can only be 1 for the sandbox currently (but any number for classical MC)\n" + output += "group_dipole = 1\t// not relevant for sandbox but classical MC\n" + output += "rand_independent=0\t// not relevant for sandbox but classical MC\n" + output += "bond_range = 4\t\t// bond range AD default\n" + output += "\n" + output += get_dihedrals_block(molsetup, indexmap, name) + + return output, indexmap + + +def get_dihedrals_block(molsetup, indexmap, name): + # molsetup.dihedral_interactions is a list of unique fourier_series + # molsetup.dihedral_partaking_atoms has tuples of atom indices as keys, and the values + # are the indices in molsetup.dihedral_interactions + # molsetup.dihedral_labels also has tuples of atom indices as keys, but the + # values are not guaranteed to be unique + + # Let's carefully use dihedral_labels to name the interactions + label_by_index = {} + atomidx_by_index = {} + for atomidx in molsetup.dihedral_partaking_atoms: + a, b, c, d = atomidx + if ( + molsetup.get_is_ignore(a) + or molsetup.get_is_ignore(b) + or molsetup.get_is_ignore(c) + or molsetup.get_is_ignore(d) + ): + continue + bond_id = Bond.get_bond_id(b, c) + if not molsetup.bond_info[bond_id].rotatable: + continue + index = molsetup.dihedral_partaking_atoms[atomidx] + atomidx_by_index.setdefault(index, set()) + atomidx_by_index[index].add(atomidx) + label = ( + molsetup.dihedral_labels[atomidx] + if atomidx in molsetup.dihedral_labels + else None + ) + if label is None: + label = "from_meeko_%d" % index + label_by_index.setdefault(index, set()) + label_by_index[index].add(label) + spent_labels = set() + for index in label_by_index: + label = "_".join(label_by_index[index]) + number = 0 + while label in spent_labels: + number += 1 + label = "_".join(label_by_index[index]) + "_v%d" % number + label_by_index[index] = label + spent_labels.add(label) + + text = "" + for index in label_by_index: + text += "[Interaction: %s, %s]\n" % (name, label_by_index[index]) + text += "type = dihedral\n" + atomidx_strings = [] + for atomidx in atomidx_by_index[index]: + string = ",".join(["%d" % (indexmap[i] + 1) for i in atomidx]) + atomidx_strings.append(string) + text += "elements = {%s}\n" % ("|".join(atomidx_strings)) + text += "parameters = %s\n" % _aux_fourier_conversion( + molsetup.dihedral_interactions[index] + ) + text += "\n" + return text + + +def _aux_fourier_conversion(fourier_series): + # convert from: + # k*(1+cos(n*theta-phase)) + # to: + # (k/2)*(1+cos(n*(theta+phase))) + # where n = periodicity + max_periodicity = max([fs["periodicity"] for fs in fourier_series]) + tmp = [(0, 0)] * max_periodicity + for fs in fourier_series: + i = fs["periodicity"] - 1 + k = 2.0 * fs["k"] + phase = -1 * fs["phase"] + tmp[i] = (k, phase) + strings = [] + periodicity = 0 + for k, phase in tmp: + periodicity += 1 + k_str = "0" + if phase == 0: + phase_str = "0" + else: + phase_str = ("%f" % (phase / np.pi)).rstrip("0").rstrip(".") + "*pi" + if phase_str == "1*pi": + phase_str = "pi" + if phase_str == "-1*pi": + phase_str = "-pi" + if periodicity != 1: + phase_str += "/%d" % periodicity + if k != 0: + k_str = "%f*4.184/60.221" % (k) + strings.append("%s,%s" % (k_str, phase_str)) + return "(" + ";".join(strings) + ")" + + +class PDBQTWriterLegacy: + + @staticmethod + def _get_pdbinfo_fitting_pdb_chars(pdbinfo): + """return strings and integers that are guaranteed + to fit within the designated chars of the PDB format""" + + atom_name = pdbinfo.name + res_name = pdbinfo.resName + res_num = pdbinfo.resNum + chain = pdbinfo.chain + if len(atom_name) > 4: + atom_name = atom_name[0:4] + if len(res_name) > 3: + res_name = res_name[0:3] + if res_num > 9999: + res_num = res_num % 10000 + if len(chain) > 1: + chain = chain[0:1] + return atom_name, res_name, res_num, chain + + @classmethod + def _make_pdbqt_line_from_molsetup(cls, setup, atom_idx, count): + """ """ + pdbinfo = setup.get_pdbinfo(atom_idx) + if pdbinfo is None: + pdbinfo = pdbutils.PDBAtomInfo("", "", 0, "") + atom_name, res_name, res_num, chain = cls._get_pdbinfo_fitting_pdb_chars( + pdbinfo + ) # TODO icode + coord = setup.get_coord(atom_idx) + atom_type = setup.get_atom_type(atom_idx) + charge = setup.get_charge(atom_idx) + pdbqt_line = cls._make_pdbqt_line( + count, atom_name, res_name, chain, res_num, coord, charge, atom_type + ) + return pdbqt_line + + @staticmethod + def _make_pdbqt_line( + count, atom_name, res_name, chain, res_num, coord, charge, atom_type, icode="" + ): + record_type = "ATOM" + alt_id = " " + occupancy = 1.0 + temp_factor = 0.0 + atom = "{:6s}{:5d} {:^4s}{:1s}{:3s} {:1s}{:4d}{:1s} {:8.3f}{:8.3f}{:8.3f}{:6.2f}{:6.2f} {:6.3f} {:<2s}" + pdbqt_line = atom.format( + record_type, + count, + atom_name, + alt_id, + res_name, + chain, + res_num, + icode, + float(coord[0]), + float(coord[1]), + float(coord[2]), + occupancy, + temp_factor, + charge, + atom_type, + ) + return pdbqt_line + + @classmethod + def _walk_graph_recursive(cls, setup, node, data, edge_start=0, first=False): + """recursive walk of rigid bodies""" + + if first: + data["pdbqt_buffer"].append("ROOT") + member_pool = sorted(setup.flexibility_model["rigid_body_members"][node]) + else: + member_pool = setup.flexibility_model["rigid_body_members"][node][:] + member_pool.remove(edge_start) + member_pool = [edge_start] + member_pool + + for member in member_pool: + if setup.get_is_ignore(member) == 1: + continue + pdbqt_line = cls._make_pdbqt_line_from_molsetup( + setup, member, data["count"] + ) + data["pdbqt_buffer"].append(pdbqt_line) + data["numbering"][member] = data["count"] # count starts at 1 + data["count"] += 1 + + if first: + data["pdbqt_buffer"].append("ENDROOT") + + data["visited"].append(node) + + for neigh in setup.flexibility_model["rigid_body_graph"][node]: + if neigh in data["visited"]: + continue + + # Write the branch + begin, next_index = setup.flexibility_model["rigid_body_connectivity"][ + node, neigh + ] + + # do not write branch (or anything downstream) if any of the two atoms + # defining the rotatable bond are ignored + if setup.get_is_ignore(begin) or setup.get_is_ignore(next_index): + continue + + begin = data["numbering"][begin] + end = data["count"] + + data["pdbqt_buffer"].append("BRANCH %3d %3d" % (begin, end)) + data = cls._walk_graph_recursive(setup, neigh, data, edge_start=next_index) + data["pdbqt_buffer"].append("ENDBRANCH %3d %3d" % (begin, end)) + + return data + + @staticmethod + def _is_molsetup_ok(setup, bad_charge_ok): + + success = True + error_msg = "" + + if len(setup.restraints): + error_msg = "molsetup has restraints but these can't be written to PDBQT" + success = False + + for atom in setup.atoms: + if atom.is_ignore: + continue + if atom.atom_type is None: + error_msg += "atom number %d has None type, mol name: %s\n" % ( + atom.index, + setup.get_mol_name(), + ) + success = False + for atom in setup.atoms: + if atom.is_ignore: + continue + if atom.atom_type is None: + error_msg += "atom number %d has None type, mol name: %s\n" % ( + atom.index, + setup.get_mol_name(), + ) + success = False + c = atom.charge + if not bad_charge_ok and ( + not isinstance(c, (float, int)) or math.isnan(c) or math.isinf(c) + ): + error_msg += ( + "atom number %d has non finite charge, mol name: %s, charge: %s\n" + % (atom.index, setup.get_mol_name(), str(c)) + ) + success = False + + return success, error_msg + + @classmethod + def write_string_from_polymer(cls, polymer): + rigid_pdbqt_string, flex_pdbqt_dict = cls.write_from_polymer( + polymer + ) + flex_pdbqt_string = "" + for pdbqt_string in flex_pdbqt_dict.values(): + flex_pdbqt_string += pdbqt_string + return rigid_pdbqt_string, flex_pdbqt_string + + @classmethod + def write_from_polymer(cls, polymer): + rigid_pdbqt_string = "" + flex_pdbqt_dict = {} + atom_count = 0 + flex_atom_count = 0 + for res_id, monomer in polymer.get_valid_monomers().items(): + chain, resnum = res_id.split(":") + if resnum[-1].isalpha(): + icode = resnum[-1] + resnum = int(resnum[:-1]) + else: + icode = "" + resnum = int(resnum) + molsetup = monomer.molsetup + resname = monomer.input_resname + if monomer.is_movable: + original_ignore = {atom.index: atom.is_ignore for atom in molsetup.atoms} + graph = molsetup.flexibility_model["rigid_body_graph"] + root = molsetup.flexibility_model["root"] + if len(graph[root]) != 1: + raise RuntimeError( + f"flexible residue {res_id} has {len(graph[root])}" + " rotatable bonds from root, but PDBQT is limited to 1" + ) + # set ignore to True for static atoms of flexible sidechains + # to exclude them from the PDBQT string + for atom_idx, is_flex in enumerate(monomer.is_flexres_atom): + molsetup.atoms[atom_idx].is_ignore = not is_flex + this_flex_pdbqt, ok, err = PDBQTWriterLegacy.write_string( + molsetup, remove_smiles=True, add_index_map=True + ) + for atom in molsetup.atoms: + atom.is_ignore = original_ignore[atom.index] + if not ok: + raise RuntimeError(err) + this_flex_pdbqt, flex_atom_count = ( + cls.adapt_pdbqt_for_autodock4_flexres( + this_flex_pdbqt, + resname, + chain, + int(resnum), + skip_rename_ca_cb=True, + atom_count=flex_atom_count, + ) + ) + flex_pdbqt_dict[res_id] = this_flex_pdbqt + + for atom_idx, atom in enumerate(molsetup.atoms): + if atom.is_ignore or monomer.is_flexres_atom[atom_idx]: + continue + atom_type = atom.atom_type + coord = atom.coord + atom_name = atom.pdbinfo.name + charge = atom.charge + atom_count += 1 + rigid_pdbqt_string += ( + cls._make_pdbqt_line( + atom_count, + atom_name, + resname, + chain, + resnum, + coord, + charge, + atom_type, + icode, + ) + + eol + ) + return rigid_pdbqt_string, flex_pdbqt_dict + + @classmethod + def write_string( + cls, setup, add_index_map=False, remove_smiles=False, bad_charge_ok=False + ): + """Output a PDBQT file as a string. + + Args: + setup: RDKitMoleculeSetup + + Returns: + str: PDBQT string of the molecule + bool: success + str: error message + """ + + success, error_msg = cls._is_molsetup_ok(setup, bad_charge_ok) + if not success: + pdbqt_string = "" + return pdbqt_string, success, error_msg + + data = { + "visited": [], + "numbering": {}, + "pdbqt_buffer": [], + "count": 1, + } + + torsdof = len(setup.flexibility_model["rigid_body_graph"]) - 1 + + if "torsions_org" in setup.flexibility_model: + torsdof_org = setup.flexibility_model["torsions_org"] + data["pdbqt_buffer"].append( + "REMARK Flexibility Score: %8.3f" % setup.flexibility_model["score"] + ) + active_tors = torsdof_org + else: + active_tors = torsdof + + data = cls._walk_graph_recursive( + setup, setup.flexibility_model["root"], data, first=True + ) + + if add_index_map: + for i, remark_line in enumerate( + cls.remark_index_map(setup, data["numbering"]) + ): + # Need to use 'insert' because data["numbering"] + # is populated in self._walk_graph_recursive. + data["pdbqt_buffer"].insert(i, remark_line) + + if not remove_smiles: + smiles, order = setup.get_smiles_and_order() + missing_h = [] # hydrogens which are not in the smiles + strings_h_parent = [] + for key in data["numbering"]: + if setup.atoms[key].is_pseudo_atom: + continue + if key not in order: + if setup.get_atomic_num(key) != 1: + error_msg += ( + "non-Hydrogen atom unexpectedely missing from smiles!?" + ) + error_msg += " (mol name: %s)\n" % setup.get_mol_name() + pdbqt_string = "" + success = False + return pdbqt_string, success, error_msg + missing_h.append(key) + parents = setup.get_neighbors(key) + parents = [ + i for i in parents if i < setup.true_atom_count + ] # exclude pseudos + if len(parents) != 1: + error_msg += ( + f"expected hydrogen {key} to be bonded to exactly one atom" + f" but it's bonded to {parents}" + ) + error_msg += " (mol name: %s)\n" % setup.get_mol_name() + pdbqt_string = "" + success = False + return pdbqt_string, success, error_msg + parent_idx = order[parents[0]] # already 1-indexed + string = " %d %d" % ( + parent_idx, + data["numbering"][key], + ) # key 0-indexed; _numbering[key] 1-indexed + strings_h_parent.append(string) + remarks_h_parent = cls.break_long_remark_lines( + strings_h_parent, "REMARK H PARENT" + ) + remark_prefix = "REMARK SMILES IDX" + remark_idxmap = cls.remark_index_map( + setup, data["numbering"], order, remark_prefix, missing_h + ) + remarks = [] + remarks.append("REMARK SMILES %s" % smiles) # break line at 79 chars? + remarks.extend(remark_idxmap) + remarks.extend(remarks_h_parent) + + for i, remark_line in enumerate(remarks): + # Need to use 'insert' because data["numbering"] + # is populated in self._walk_graph_recursive. + data["pdbqt_buffer"].insert(i, remark_line) + + # torsdof is always going to be the one of the rigid, non-macrocyclic one + data["pdbqt_buffer"].append("TORSDOF %d" % active_tors) + + pdbqt_string = eol.join(data["pdbqt_buffer"]) + eol + return pdbqt_string, success, error_msg + + @classmethod + def remark_index_map( + cls, setup, numbering, order=None, prefix="REMARK INDEX MAP", missing_h=() + ): + """write mapping of atom indices from input molecule to output PDBQT + order[ob_index(i.e. 'key')] = smiles_index + """ + + if order is None: + order = {key: key + 1 for key in numbering} # key+1 breaks OB + # max_line_length = 79 + # remark_lines = [] + # line = prefix + strings = [] + for key in numbering: + if setup.atoms[key].is_pseudo_atom: + continue + if key in missing_h: + continue + string = " %d %d" % (order[key], numbering[key]) + strings.append(string) + return cls.break_long_remark_lines(strings, prefix) + # candidate_text = " %d %d" % (order[key], self._numbering[key]) + # if (len(line) + len(candidate_text)) < max_line_length: + # line += candidate_text + # else: + # remark_lines.append(line) + # line = 'REMARK INDEX MAP' + candidate_text + # remark_lines.append(line) + # return remark_lines + + @staticmethod + def break_long_remark_lines(strings, prefix, max_line_length=79): + remarks = [prefix] + for string in strings: + if (len(remarks[-1]) + len(string)) < max_line_length: + remarks[-1] += string + else: + remarks.append(prefix + string) + return remarks + + @staticmethod + def adapt_pdbqt_for_autodock4_flexres( + pdbqt_string, res, chain, num, skip_rename_ca_cb=False, atom_count=None + ): + """adapt pdbqt_string to be compatible with AutoDock4 requirements: + - first and second atoms named CA and CB + - write BEGIN_RES / END_RES + - remove TORSDOF + this is for covalent docking (tethered) + """ + new_string = "BEGIN_RES %s %s %s" % (res, chain, num) + eol + atom_number = 0 + offset = atom_count + for line in pdbqt_string.split(eol): + if line == "": + continue + if line.startswith("TORSDOF"): + continue + if line.startswith("ATOM"): + if not skip_rename_ca_cb: + atom_number += 1 + if atom_number == 1: + line = line[:13] + "CA" + line[15:] + elif atom_number == 2: + line = line[:13] + "CB" + line[15:] + if atom_count is not None: + atom_count += 1 + n = "%5d" % atom_count + n = n[:5] + line = line[:6] + n + line[11:] + new_string += line + eol + continue + elif offset is not None and ( + line.startswith("BRANCH") or line.startswith("ENDBRANCH") + ): + keyword, i, j = line.split() + new_string += ( + f"{keyword} {int(i)+offset:3d} {int(j)+offset:3d}" + eol + ) + continue + new_string += line + eol + new_string += "END_RES %s %s %s" % (res, chain, num) + eol + if atom_count is None: + return new_string # just keeping backwards compatibility + else: + return new_string, atom_count diff --git a/src/moldrug/utils.py b/src/moldrug/utils.py index da4d34c..4c85339 100644 --- a/src/moldrug/utils.py +++ b/src/moldrug/utils.py @@ -20,8 +20,8 @@ import pandas as pd import tqdm from crem.crem import grow_mol, mutate_mol -from meeko import (MoleculePreparation, PDBQTMolecule, PDBQTWriterLegacy, - RDKitMolCreate) +from moldrug.meeko_light import (MoleculePreparation, + PDBQTMolecule, PDBQTWriterLegacy, RDKitMolCreate) from rdkit import Chem, RDLogger from rdkit.Chem import AllChem, DataStructs, Descriptors, Lipinski, rdFMCS @@ -1561,8 +1561,7 @@ def __call__(self, njobs: int = 1): probs = softmax((-self.beta * np.array(self.pop)).astype('float64')) if any(np.isnan(probs)): probs = np.nan_to_num(probs) - - + # TODO: This cycle should run in this way only if no user generetor was provided # with and if, else statment I could correct, and then the genereator functions is completlly up to the user, # then I do not need to worry in how the selection is made,