devbox
黑客灵魂
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

ACPYPE中FAILED: [Errno 2] No such file or directory: 'FFF_AC.prmtop'的解决方法_acpype failed: [errno 2] no such file or directory

作者：黑客灵魂 | 2024-07-30 09:42:04

赞

踩

acpype failed: [errno 2] no such file or directory: 'tmp

大家如果在github上面下载ACPYPE （https://github.com/t-/acpype），放入到LINUX上面安装后，百分百报错，出现的问题如下：

[root@pdynamo test]# ../acpype.py -i FFF.pdb
============================================================================
| ACPYPE: AnteChamber PYthon Parser interfacE v. 0 0 Rev: 0 (c) 2020 AWSdS |
============================================================================
==> ... charge set to 0
==> ... converting pdb input file to mol2 input file
==> * Babel OK *
==> Executing Antechamber...
==> * Antechamber OK *
++++++++++start_quote+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
/bin/sh: -c: line 0: syntax error near unexpected token `('
/bin/sh: -c: line 0: `which: no parmchk in (/usr/local/bin:/home/peter/Amber/amber18/bin:/home/peter/gmx2019.4/bin:/home/peter/fftw-3.3.8/bin:/home/peter/openmpi-3.1.3/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin/intel64:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/libfabric/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin/intel64:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/libfabric/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin/intel64:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/libfabric/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/bin:/home/peter/intel/debugger_2019/gdb/intel64/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/home/peter/Amber/amber18/bin:/root/bin) -i FFF_bcc_gaff.mol2 -f mol2 -o FFF_AC.frcmod'
++++++++++end_quote+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ERROR: Parmchk failed
ERROR: Tleap failed
==> ... trying Sleap
==> Executing Sleap...
++++++++++start_quote+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

++++++++++end_quote+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++start_quote+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
/bin/sh: -c: line 0: syntax error near unexpected token `('
/bin/sh: -c: line 0: `which: no sleap in (/usr/local/bin:/home/peter/Amber/amber18/bin:/home/peter/gmx2019.4/bin:/home/peter/fftw-3.3.8/bin:/home/peter/openmpi-3.1.3/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin/intel64:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/libfabric/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin/intel64:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/libfabric/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin/intel64:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/libfabric/bin:/home/peter/intel/compilers_and_libraries_2019.5.281/linux/mpi/intel64/bin:/home/peter/intel/debugger_2019/gdb/intel64/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/home/peter/Amber/amber18/bin:/root/bin) -f sleap.in'

++++++++++end_quote+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ERROR: Sleap failed
==> Removing temporary files...
ACPYPE FAILED: [Errno 2] No such file or directory: 'FFF_AC.prmtop'
Total time of execution: 3m 8s

其原因是这个源文件有毛病，请大家把底下这个复制，然后覆盖acpepy.py这个文件，直接就好用。

这个文件感谢计算化学领域国内最强的卢天大神提供：http://bbs.keinsci.com/thread-14021-1-1.html

  1. #!/usr/bin/env python
  2.  
  3. from __future__ import print_function
  4. from datetime import datetime
  5. from shutil import copy2
  6. from shutil import rmtree
  7. import traceback
  8. import signal
  9. import time
  10. import optparse
  11. import math
  12. import os
  13. import pickle
  14. import sys
  15. import subprocess as sub
  16. import re
  17.  
  18. """
  19.     Requirements: Python 2.6 or higher or Python 3.x
  20.                   Antechamber (from AmberTools preferably)
  21.                   OpenBabel (optional, but strongly recommended)
  23.     This code is released under GNU General Public License V3.
  25.           <<<  NO WARRANTY AT ALL!!!  >>>
  27.     It was inspired by:
  29.     - amb2gmx.pl (Eric Sorin, David Mobley and John Chodera)
  30.       and depends on Antechamber and Openbabel
  32.     - YASARA Autosmiles:
  33.       http://www.yasara.org/autosmiles.htm (Elmar Krieger)
  35.     - topolbuild (Bruce Ray)
  37.     - xplo2d (G.J. Kleywegt)
  39.     For Antechamber, please cite:
  40.     1.  Wang, J., Wang, W., Kollman P. A.; Case, D. A. "Automatic atom type and
  41.         bond type perception in molecular mechanical calculations". Journal of
  42.         Molecular Graphics and Modelling , 25, 2006, 247260.
  43.     2.  Wang, J., Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A.
  44.         "Development and testing of a general AMBER force field". Journal of
  45.         Computational Chemistry, 25, 2004, 1157-1174.
  47.     If you use this code, I am glad if you cite:
  49.     SOUSA DA SILVA, A. W. & VRANKEN, W. F.
  50.     ACPYPE - AnteChamber PYthon Parser interfacE.
  51.     BMC Research Notes 2012, 5:367 doi:10.1186/1756-0500-5-367
  52.     http://www.biomedcentral.com/1756-0500/5/367
  54.     Alan Wilter Sousa da Silva, D.Sc.
  55.     Bioinformatician, UniProt, EMBL-EBI
  56.     Hinxton, Cambridge CB10 1SD, UK.
  57.     >>http://www.ebi.ac.uk/~awilter<<
  59.     alanwilter _at_ gmail _dot_ com
  60. """
  61.  
  62. svnId = '$Id$'
  63. try:
  64.     svnRev, svnDate, svnTime = svnId.split()[2:5]
  65. except:
  66.     svnRev, svnDate, svnTime = '0', '0', '0'
  67. year = datetime.today().year
  68. tag = "%s %s Rev: %s" % (svnDate, svnTime, svnRev)
  69.  
  70. lineHeader = \
  71.     '''
  72. | ACPYPE: AnteChamber PYthon Parser interfacE v. %s (c) %s AWSdS |
  73. ''' % (tag, year)
  74. frameLine = (len(lineHeader) - 2) * '='
  75. header = '%s%s%s' % (frameLine, lineHeader, frameLine)
  76.  
  77. #    TODO:
  78. #        Howto Charmm and Amber with NAMD
  79. #        Howto build topology for a modified amino acid
  80. #        CYANA topology files
  81.  
  82. # List of Topology Formats created by acpype so far:
  83. outTopols = ['gmx', 'cns', 'charmm']
  84. qDict = {'mopac': 0, 'divcon': 1, 'sqm': 2}
  85.  
  86. # Residues that are not solute, to be avoided when balancing charges in
  87. # amb2gmx mode
  88. ionOrSolResNameList = ['Cl-', 'Na+', 'K+', 'CIO', 'Cs+', 'IB', 'Li+', 'MG2',
  89.                        'Rb+', 'WAT', 'MOH', 'NMA']
  90.  
  91. # leapAmberFile = 'leaprc.ff99SB'  # 'leaprc.ff10' and 'leaprc.ff99bsc0' has extra Atom Types not in parm99.dat
  92. leapAmberFile = 'leaprc.protein.ff14SB'  # 'leaprc.ff14SB'
  93.  
  94. # "qm_theory='AM1', grms_tol=0.0002, maxcyc=999, tight_p_conv=1, scfconv=1.d-10,"
  95. # "AM1 ANALYT MMOK GEO-OK PRECISE"
  96.  
  97. cal = 4.184
  98. Pi = 3.141593
  99. qConv = 18.222281775  # 18.2223
  100. radPi = 57.295780  # 180/Pi
  101. maxDist = 3.0
  102. minDist = 0.5
  103. maxDist2 = maxDist ** 2  # squared Ang.
  104. minDist2 = minDist ** 2  # squared Ang.
  105. diffTol = 0.01
  106.  
  107. dictAmbAtomType2AmbGmxCode = \
  108.     {'BR': '1', 'C': '2', 'CA': '3', 'CB': '4', 'CC': '5', 'CK': '6', 'CM': '7', 'CN': '8', 'CQ': '9',
  109.      'CR': '10', 'CT': '11', 'CV': '12', 'CW': '13', 'C*': '14', 'Ca': '15', 'F': '16', 'H': '17',
  110.      'HC': '18', 'H1': '19', 'H2': '20', 'H3': '21', 'HA': '22', 'H4': '23', 'H5': '24', 'HO': '25',
  111.      'HS': '26', 'HW': '27', 'HP': '28', 'I': '29', 'Cl': '30', 'Na': '31', 'IB': '32', 'Mg': '33',
  112.      'N': '34', 'NA': '35', 'NB': '36', 'NC': '37', 'N2': '38', 'N3': '39', 'N*': '40', 'O': '41',
  113.      'OW': '42', 'OH': '43', 'OS': '44', 'O2': '45', 'P': '46', 'S': '47', 'SH': '48', 'CU': '49',
  114.      'FE': '50', 'K': '51', 'Rb': '52', 'Cs': '53', 'Li': '56', 'Zn': '57', 'Sr': '58', 'Ba': '59',
  115.      'MCH3A': 'MCH3A', 'MCH3B': 'MCH3B', 'MNH2': 'MNH2', 'MNH3': 'MNH3', 'MW': 'MW'}
  116.  
  117. dictOplsAtomType2OplsGmxCode = \
  118.     {'Ac3+': ['697'], 'Am3+': ['699'], 'Ar': ['Ar', '097'], 'Ba2+': ['414'],
  119.      'Br': ['722', '730'], 'Br-': ['402'],
  120.      'CT': ['064', '076', '122', '135', '136', '137', '138', '139', '148', '149', '152', '157', '158', '159', '161', '173', '174', '175', '181', '182', '183', '184', '206', '207', '208', '209', '210', '211', '212', '213', '214', '215', '216', '217', '218', '219', '220', '223', '224', '225', '229', '230', '242', '243', '244', '256', '257', '258', '259', '273', '274', '275', '276', '291', '292', '293', '294', '297', '305', '306', '307', '308', '331', '371', '373', '375', '391', '396', '421', '431', '443', '448', '453', '455', '458', '461', '468', '476', '482', '484', '486', '490', '491', '492', '498', '499', '505', '515', '516', '645', '670', '671', '672', '673', '674', '675', '676', '677', '678', '679', '680', '681', '701', '725', '747', '748', '755', '756', '757', '758', '762', '764', '765', '766', '774', '775', '776', '782', '783', '903', '904', '905', '906', '907', '908', '912', '913', '914', '915', '942', '943', '944', '945', '951', '957', '959', '960', '961', '962', '963', '964'],
  121.      'CA': ['053', '145', '147', '166', '199', '221', '228', '260', '263', '266', '302', '312', '315', '317', '336', '351', '362', '380', '457', '460', '463', '472', '488', '521', '522', '523', '528', '532', '533', '538', '539', '551', '582', '590', '591', '592', '593', '604', '605', '606', '607', '608', '609', '610', '611', '612', '625', '644', '647', '648', '649', '650', '651', '652', '714', '716', '718', '720', '724', '727', '729', '731', '735', '736', '737', '738', '739', '742', '752', '768', '916', '917', '918'],
  122.      'C3': ['007', '010', '036', '039', '063', '065', '067', '068', '069', '070', '080', '088', '090', '092', '096', '106', '107', '109', '126', '132', '415', '418', '425', '429'],
  123.      'C': ['001', '017', '026', '058', '095', '131', '231', '234', '235', '247', '252', '267', '320', '322', '334', '366', '378', '470', '471', '772', '952'],
  124.      'C2': ['005', '009', '015', '016', '019', '022', '027', '028', '031', '034', '037', '056', '057', '061', '071', '081', '089', '091', '093', '110'],
  125.      'CT_2': ['223B', '224B', '225B', '246', '283', '284', '285', '292B', '293B', '295', '298', '299', '906B', '912B'],
  126.      'CM': ['141', '142', '143', '227', '323', '324', '337', '338', '381', '382', '517', '518', '708'],
  127.      'CW': ['508', '514', '543', '552', '561', '567', '575', '583', '588', '637'],
  128.      'CB': ['050', '349', '350', '364', '365', '501', '595', '623', '624'],
  129.      'CH': ['006', '008', '014', '025', '029', '030', '060', '073'],
  130.      'CZ': ['261', '423', '754', '925', '927', '928', '929', '931'],
  131.      'CO': ['189', '191', '193', '195', '197', '198'],
  132.      'C_2': ['232', '233', '277', '280', '465'],
  133.      'CR': ['506', '509', '558', '572', '634'],
  134.      'CQ': ['347', '531', '621', '642'],
  135.      'CV': ['507', '560', '574', '636'],
  136.      'CY': ['711', '712', '713', '733'],
  137.      'CS': ['544', '568', '589'],
  138.      'CK': ['353', '627'], 'CN': ['502', '594'], 'CP': ['043', '048'], 'CU': ['550', '581'],
  139.      'CT_3': ['245', '296'], 'C=': ['150', '178'], 'CD': ['011', '075'],
  140.      'C4': ['066'], 'C7': ['077'], 'C8': ['074'], 'C9': ['072'], 'CX': ['510'],
  141.      'C!': ['145B'], 'C*': ['500'], 'C+': ['700'], 'C_3': ['271'],
  142.      'CC': ['045'], 'CF': ['044'], 'CG': ['049'], 'CT_4': ['160'],
  143.      'Ca2+': ['412'],
  144.      'Cl': ['123', '151', '226', '264'],
  145.      'Cl-': ['401', '709'],
  146.      'Cs+': ['410'], 'Cu2+': ['Cu2+'], 'Eu3+': ['705'],
  147.      'F': ['164', '719', '721', '726', '728', '786', '956', '965'],
  148.      'F-': ['400'], 'Fe2+': ['Fe2+'], 'Gd3+': ['706'],
  149.      'HA': ['146', '316', '318', '389', '524', '525', '526', '529', '534', '535', '536', '540', '541', '546', '547', '554', '555', '556', '563', '564', '565', '569', '570', '576', '577', '578', '584', '585', '586', '597', '598', '599', '600', '601', '602', '613', '614', '615', '616', '617', '618', '619', '629', '630', '631', '638', '639', '640', '643', '653', '654', '655', '656', '715', '717', '740', '741', '746'],
  150.      'HC': ['140', '144', '153', '156', '165', '176', '185', '190', '192', '194', '196', '255', '279', '282', '329', '330', '332', '344', '372', '374', '376', '392', '416', '419', '422', '426', '430', '432', '444', '449', '454', '456', '459', '462', '469', '477', '483', '485', '487', '702', '710', '759', '763', '777', '778', '779', '784', '911', '926', '930', '950', '958'],
  151.      'H': ['004', '013', '041', '047', '128', '240', '241', '250', '254', '314', '325', '327', '339', '342', '343', '357', '358', '360', '367', '369', '383', '385', '387', '388', '428', '479', '481', '504', '513', '545', '553', '562', '596', '632', '744', '745', '909', '910'],
  152.      'H3': ['021', '052', '055', '104', '105', '289', '290', '301', '304', '310', '941', '955'],
  153.      'HO': ['024', '079', '155', '163', '168', '170', '172', '188', '270', '435'],
  154.      'HS': ['033', '086', '087', '204', '205'],
  155.      'HW': ['112', '114', '117', '119', '796'],
  156.      'H4': ['345', '390'], 'H5': ['355', '359'],
  157.      'He': ['130'], 'I': ['732'], 'I-': ['403'], 'K+': ['408'], 'Kr': ['098'],
  158.      'LP': ['433', '797'], 'La3+': ['703'], 'Li+': ['404', '406'],
  159.      'MCH3A': ['MCH3A'], 'MCH3B': ['MCH3B'], 'MNH2': ['MNH2'], 'MNH3': ['MNH3'],
  160.      'MW': ['MW', '115'], 'Mg2+': ['411'],
  161.      'NA': ['040', '046', '319', '321', '333', '354', '361', '377', '379', '503', '512', '542', '548', '557', '587', '628'],
  162.      'NC': ['311', '335', '346', '348', '363', '520', '527', '530', '537', '603', '620', '622', '641', '646'],
  163.      'N': ['003', '012', '094', '237', '238', '239', '249', '251', '265', '478', '480', '787'],
  164.      'N3': ['020', '101', '102', '103', '286', '287', '288', '309', '427', '940', '953'],
  165.      'N2': ['051', '054', '300', '303', '313', '341', '356', '368', '386', '743'],
  166.      'NB': ['042', '352', '511', '549', '559', '573', '580', '626', '635'],
  167.      'N*': ['319B', '333B', '354B', '377B'],
  168.      'NT': ['127', '900', '901', '902'],
  169.      'NZ': ['262', '424', '750', '753'],
  170.      'NO': ['760', '767'], 'NY': ['749', '751'],
  171.      'Na+': ['405', '407'], 'Nd3+': ['704'], 'Ne': ['129'],
  172.      'OS': ['062', '108', '179', '180', '186', '395', '442', '447', '452', '467', '473', '566', '571', '579', '773'],
  173.      'O': ['002', '059', '236', '248', '253', '326', '328', '340', '370', '384', '771', '788'],
  174.      'OH': ['023', '078', '154', '162', '167', '169', '171', '187', '268', '420', '434'],
  175.      'O2': ['018', '125', '272', '394', '441', '446', '451', '954'],
  176.      'OW': ['111', '113', '116', '118', '795'],
  177.      'O_2': ['278', '281', '466'],
  178.      'OY': ['475', '494', '497'],
  179.      'OL': ['120'], 'ON': ['761'], 'OU': ['437'], 'O_3': ['269'],
  180.      'P': ['393', '440', '445', '450', '785'],
  181.      'P+': ['781'], 'Rb+': ['409'],
  182.      'S': ['035', '038', '084', '085', '124', '202', '203', '222', '633'],
  183.      'SH': ['032', '082', '083', '200', '201', '417', '734'],
  184.      'SI': ['SI'], 'SY': ['474'], 'SY2': ['493'], 'SZ': ['496'], 'Sr2+': ['413'],
  185.      'Th4+': ['698'], 'U': ['436'], 'Xe': ['099'], 'Yb3+': ['707'], 'Zn2+': ['Zn2+']}
  186.  
  187. # reverse dictOplsAtomType2OplsGmxCode
  188. oplsCode2AtomTypeDict = {}
  189. for k, v in list(dictOplsAtomType2OplsGmxCode.items()):
  190.     for code in v:
  191.         oplsCode2AtomTypeDict[code] = k
  192. #        if code in oplsCode2AtomTypeDict.keys():
  193. #            oplsCode2AtomTypeDict[code].append(k)
  194. #        else:
  195. #            oplsCode2AtomTypeDict[code] = [k]
  196.  
  197. # Cross dictAmbAtomType2AmbGmxCode with dictOplsAtomType2OplsGmxCode & add H1,HP,H2
  198. dictAtomTypeAmb2OplsGmxCode = {'H1': ['140', '1.00800'], 'HP': ['140', '1.00800'], 'H2': ['140', '1.00800']}
  199. dictOplsMass = {'SY2': ['32.06000'], 'Zn2+': ['65.37000'], 'CQ': ['12.01100'], 'CP': ['12.01100'], 'Nd3+': ['144.24000'], 'Br-': ['79.90400'], 'Cu2+': ['63.54600'], 'Br': ['79.90400'], 'H': ['1.00800'], 'P': ['30.97376'], 'Sr2+': ['87.62000'], 'ON': ['15.99940'], 'OL': ['0.00000'], 'OH': ['15.99940'], 'OY': ['15.99940'], 'OW': ['15.99940'], 'OU': ['15.99940'], 'OS': ['15.99940'], 'Am3+': ['243.06000'], 'HS': ['1.00800'], 'HW': ['1.00800'], 'HO': ['1.00800'], 'HC': ['1.00800'], 'HA': ['1.00800'], 'O2': ['15.99940'], 'Ca2+': ['40.08000'], 'Th4+': ['232.04000'], 'He': ['4.00260'], 'C': ['12.01100'], 'Cs+': ['132.90540'], 'O': ['15.99940'], 'Gd3+': ['157.25000'], 'S': ['32.06000'], 'P+': ['30.97376'], 'La3+': ['138.91000'], 'H3': ['1.00800'], 'H4': ['1.00800'], 'MNH2': ['0.00000'], 'MW': ['0.00000'], 'NB': ['14.00670'], 'K+': ['39.09830'], 'Ne': ['20.17970'], 'Rb+': ['85.46780'], 'C+': ['12.01100'], 'C*': ['12.01100'], 'NO': ['14.00670'], 'CT_4': ['12.01100'], 'NA': ['14.00670'], 'C!': ['12.01100'], 'NC': ['14.00670'], 'NZ': ['14.00670'], 'CT_2': ['12.01100'], 'CT_3': ['12.01100'], 'NY': ['14.00670'], 'C9': ['14.02700'], 'C8': ['13.01900'], 'C=': ['12.01100'], 'Yb3+': ['173.04000'], 'C3': ['15.03500', '12.01100'], 'C2': ['14.02700'], 'C7': ['12.01100'], 'C4': ['16.04300'], 'CK': ['12.01100'], 'Cl-': ['35.45300'], 'N*': ['14.00670'], 'CH': ['13.01900'], 'CO': ['12.01100'], 'CN': ['12.01100'], 'CM': ['12.01100'], 'F': ['18.99840'], 'CC': ['12.01100'], 'CB': ['12.01100'], 'CA': ['12.01100'], 'CG': ['12.01100'], 'CF': ['12.01100'], 'N': ['14.00670'], 'CZ': ['12.01100'], 'CY': ['12.01100'], 'CX': ['12.01100'], 'Ac3+': ['227.03000'], 'CS': ['12.01100'], 'CR': ['12.01100'], 'N2': ['14.00670'], 'N3': ['14.00670'], 'CW': ['12.01100'], 'CV': ['12.01100'], 'CU': ['12.01100'], 'CT': ['12.01100'], 'SZ': ['32.06000'], 'SY': ['32.06000'], 'Cl': ['35.45300'], 'NT': ['14.00670'], 'O_2': ['15.99940'], 'Xe': ['131.29300'], 'SI': ['28.08000'], 'SH': ['32.06000'], 'Eu3+': ['151.96000'], 'F-': ['18.99840'], 'MNH3': ['0.00000'], 'H5': ['1.00800'], 'C_3': ['12.01100'], 'C_2': ['12.01100'], 'I-': ['126.90450'], 'LP': ['0.00000'], 'I': ['126.90450'], 'Na+': ['22.98977'], 'Li+': ['6.94100'], 'U': ['0.00000'], 'MCH3A': ['0.00000'], 'MCH3B': ['0.00000'], 'CD': ['13.01900', '12.01100'], 'O_3': ['15.99940'], 'Kr': ['83.79800'], 'Fe2+': ['55.84700'], 'Ar': ['39.94800'], 'Mg2+': ['24.30500'], 'Ba2+': ['137.33000']}
  200. for ambKey in dictAmbAtomType2AmbGmxCode:
  201.     if ambKey in dictOplsAtomType2OplsGmxCode:
  202.         dictAtomTypeAmb2OplsGmxCode[ambKey] = dictOplsAtomType2OplsGmxCode[ambKey] + list(dictOplsMass[ambKey])
  203.  
  204. # learnt from 22 residues test.
  205. dictAtomTypeAmb2OplsGmxCode = {'HS': ['204', '1.008'], 'HP': ['140', '1.008'], 'HO': ['155', '168', '1.008'], 'HC': ['140', '1.008'], 'HA': ['146', '1.008'], 'O2': ['272', '15.9994'], 'C*': ['500', '12.011'], 'NA': ['503', '512', '14.0067'], 'NB': ['511', '14.0067'], 'CB': ['501', '12.011'], 'C': ['235', '271', '12.011'], 'CN': ['502', '12.011'], 'CM': ['302', '12.011'], 'CC': ['507', '508', '510', '12.011'], 'H': ['240', '241', '290', '301', '304', '310', '504', '513', '1.008'], 'CA': ['145', '166', '12.011'], 'O': ['236', '15.9994'], 'N': ['237', '238', '239', '14.0067'], 'S': ['202', '32.06'], 'CR': ['506', '509', '12.011'], 'N2': ['300', '303', '14.0067'], 'N3': ['287', '309', '14.0067'], 'CW': ['508', '510', '514', '12.011'], 'CV': ['507', '12.011'], 'CT': ['135', '136', '137', '149', '157', '158', '206', '209', '210', '223B', '224B', '245', '246', '274', '283', '284', '285', '292', '292B', '293B', '296', '307', '308', '505', '12.011'], 'OH': ['154', '167', '15.9994'], 'H1': ['140', '1.008'], 'H4': ['146', '1.008'], 'H5': ['146', '1.008'], 'SH': ['200', '32.06']}
  206.  
  207. # learnt from 22 residues test.
  208. dictAtomTypeGaff2OplsGmxCode = {'cc': ['500', '506', '507', '508', '514', '12.011'], 'ca': ['145', '166', '501', '502', '12.011'], 'h1': ['140', '1.008'], 'h4': ['146', '1.008'], 'h5': ['146', '1.008'], 'cz': ['302', '12.011'], 'c2': ['509', '510', '12.011'], 'nh': ['300', '303', '14.0067'], 'ha': ['146', '1.008'], 'na': ['503', '512', '14.0067'], 'nc': ['511', '14.0067'], 'nd': ['511', '14.0067'], 'hx': ['140', '1.008'], 'hs': ['204', '1.008'], 'hn': ['240', '241', '290', '301', '304', '310', '504', '513', '1.008'], 'ho': ['155', '168', '1.008'], 'c3': ['135', '136', '137', '149', '157', '158', '206', '209', '210', '223B', '224B', '245', '246', '274', '283', '284', '285', '292', '292B', '293B', '296', '307', '308', '505', '12.011'], 'hc': ['140', '1.008'], 'cd': ['500', '506', '507', '508', '514', '12.011'], 'c': ['235', '271', '12.011'], 'oh': ['154', '167', '15.9994'], 'ss': ['202', '32.06'], 'o': ['236', '272', '15.9994'], 'n': ['237', '238', '239', '14.0067'], 'sh': ['200', '32.06'], 'n4': ['287', '309', '14.0067']}
  209.  
  210. # draft
  211. atomTypeAmber2oplsDict = {'HS': ['HS'], 'HP': ['HC'], 'HO': ['HO'], 'HC': ['HC'],
  212.                           'HA': ['HA'], 'O2': ['O2'], 'C*': ['C*'], 'NA': ['NA'],
  213.                           'NB': ['NB'], 'CB': ['CB'], 'CN': ['CN'], 'CV': ['CV'],
  214.                           'CM': ['CA'], 'CA': ['CA'], 'CR': ['CR'], 'OH': ['OH'],
  215.                           'H1': ['HC'], 'H4': ['HA'], 'N2': ['N2'], 'N3': ['N3'],
  216.                           'H5': ['HA'], 'SH': ['SH'], 'N': ['N'], 'S': ['S'], 'O': ['O'],
  217.                           'C': ['C', 'C_3'], 'CW': ['CW', 'CX'], 'H': ['H', 'H3'],
  218.                           'CC': ['CX', 'CW', 'CV'], 'CT': ['CT', 'CT_2', 'CT_3']}
  219.  
  220. # draft
  221. a2oD = {'amber99_2': ['opls_235', 'opls_271'],
  222.         'amber99_3': ['opls_302', 'opls_145'],
  223.         'amber99_5': ['opls_507', 'opls_508', 'opls_510'],
  224.         'amber99_11': ['opls_209', 'opls_158', 'opls_283', 'opls_223B', 'opls_293B',
  225.                        'opls_284', 'opls_292B', 'opls_274', 'opls_136', 'opls_135',
  226.                        'opls_292', 'opls_157', 'opls_206', 'opls_137', 'opls_505',
  227.                        'opls_224B', 'opls_307', 'opls_308', 'opls_210', 'opls_149'],
  228.         'amber99_13': ['opls_514'],
  229.         'amber99_14': ['opls_500'],
  230.         'amber99_17': ['opls_504', 'opls_241', 'opls_240', 'opls_290', 'opls_301',
  231.                        'opls_310', 'opls_304', 'opls_513'],
  232.         'amber99_18': ['opls_140'],
  233.         'amber99_19': ['opls_140'],
  234.         'amber99_22': ['opls_146'],
  235.         'amber99_23': ['opls_146'],
  236.         'amber99_25': ['opls_155'],
  237.         'amber99_26': ['opls_204'],
  238.         'amber99_28': ['opls_140'],
  239.         'amber99_34': ['opls_238', 'opls_239', 'opls_237'],
  240.         'amber99_35': ['opls_512', 'opls_503'],
  241.         'amber99_36': ['opls_511'],
  242.         'amber99_38': ['opls_300', 'opls_303'],
  243.         'amber99_39': ['opls_309', 'opls_287'],
  244.         'amber99_41': ['opls_236'],
  245.         'amber99_43': ['opls_154'],
  246.         'amber99_45': ['opls_272'],
  247.         'amber99_47': ['opls_202'],
  248.         'amber99_48': ['opls_200'],
  249.         }
  250.  
  251. global pid
  252. pid = 0
  253.  
  254. head = "%s created by acpype (Rev: " + svnRev + ") on %s\n"
  255.  
  256. date = datetime.now().ctime()
  257.  
  258. usage = \
  259.     """
  260.     acpype -i _file_ [-c _string_] [-n _int_] [-m _int_] [-a _string_] [-f] etc. or
  261.     acpype -p _prmtop_ -x _inpcrd_ [-d]"""
  262.  
  263. epilog = \
  264.     """
  266.     output: assuming 'root' is the basename of either the top input file,
  267.             the 3-letter residue name or user defined (-b option)
  268.     root_bcc_gaff.mol2:  final mol2 file with 'bcc' charges and 'gaff' atom type
  269.     root_AC.inpcrd    :  coord file for AMBER
  270.     root_AC.prmtop    :  topology and parameter file for AMBER
  271.     root_AC.lib       :  residue library file for AMBER
  272.     root_AC.frcmod    :  modified force field parameters
  273.     root_GMX.gro      :  coord file for GROMACS
  274.     root_GMX.top      :  topology file for GROMACS
  275.     root_GMX.itp      :  molecule unit topology and parameter file for GROMACS
  276.     root_GMX_OPLS.itp :  OPLS/AA mol unit topol & par file for GROMACS (experimental!)
  277.     em.mdp, md.mdp    :  run parameters file for GROMACS
  278.     root_NEW.pdb      :  final pdb file generated by ACPYPE
  279.     root_CNS.top      :  topology file for CNS/XPLOR
  280.     root_CNS.par      :  parameter file for CNS/XPLOR
  281.     root_CNS.inp      :  run parameters file for CNS/XPLOR
  282.     root_CHARMM.rtf   :  topology file for CHARMM
  283.     root_CHARMM.prm   :  parameter file for CHARMM
  284.     root_CHARMM.inp   :  run parameters file for CHARMM"""
  285.  
  286. SLEAP_TEMPLATE = \
  287.     """
  288. source %(leapAmberFile)s
  289. source %(leapGaffFile)s
  290. set default fastbld on
  291. #set default disulfide auto
  292. %(res)s = loadpdb %(baseOrg)s.pdb
  293. #check %(res)s
  294. saveamberparm %(res)s %(acBase)s.prmtop %(acBase)s.inpcrd
  295. saveoff %(res)s %(acBase)s.lib
  296. quit
  297. """
  298.  
  299. TLEAP_TEMPLATE = \
  300.     """
  301. verbosity 1
  302. source %(leapAmberFile)s
  303. source %(leapGaffFile)s
  304. mods = loadamberparams %(acBase)s.frcmod
  305. %(res)s = loadmol2 %(acMol2FileName)s
  306. check %(res)s
  307. saveamberparm %(res)s %(acBase)s.prmtop %(acBase)s.inpcrd
  308. saveoff %(res)s %(acBase)s.lib
  309. quit
  310. """
  311.  
  312.  
  313. def dotproduct(aa, bb):
  314.     return sum((a * b) for a, b in zip(aa, bb))
  315.  
  316.  
  317. def crosproduct(a, b):
  318.     c = [a[1] * b[2] - a[2] * b[1],
  319.          a[2] * b[0] - a[0] * b[2],
  320.          a[0] * b[1] - a[1] * b[0]]
  321.     return c
  322.  
  323.  
  324. def length(v):
  325.     return math.sqrt(dotproduct(v, v))
  326.  
  327.  
  328. def vec_sub(aa, bb):
  329.     return [a - b for a, b in zip(aa, bb)]
  330.  
  331.  
  332. def imprDihAngle(a, b, c, d):
  333.     ba = vec_sub(a, b)
  334.     bc = vec_sub(c, b)
  335.     cb = vec_sub(b, c)
  336.     cd = vec_sub(d, c)
  337.     n1 = crosproduct(ba, bc)
  338.     n2 = crosproduct(cb, cd)
  339.     angle = math.acos(dotproduct(n1, n2) / (length(n1) * length(n2))) * 180 / Pi
  340.     cp = crosproduct(n1, n2)
  341.     if (dotproduct(cp, bc) < 0):
  342.         angle = -1 * angle
  343.     return angle
  344.  
  345.  
  346. def invalidArgs(text=None):
  347.     if text:
  348.         print('ERROR: ' + text)
  349.     sys.exit(1)
  350.  
  351. # verNum = string.split(sys.version)[0]
  352. verNum = re.sub('[^0-9\.]', '', sys.version.split()[0])
  353. version = verNum.split('.')  # string.split(verNum, ".")
  354. verList = list(map(int, version))
  355. if verList < [2, 6, 0]:
  356.     invalidArgs(text="Python version %s\n       Sorry, you need python 2.6 or higher" % verNum)
  357.  
  358. try:
  359.     set()
  360. except NameError:
  361.     from sets import Set as set  # @UnresolvedImport
  362.  
  363.  
  364. def elapsedTime(seconds, suffixes=['y', 'w', 'd', 'h', 'm', 's'], add_s=False, separator=' '):
  365.     """
  366.     Takes an amount of seconds and turns it into a human-readable amount of time.
  367.     """
  368.     # the formatted time string to be returned
  369.     time = []
  370.  
  371.     # the pieces of time to iterate over (days, hours, minutes, etc)
  372.     # - the first piece in each tuple is the suffix (d, h, w)
  373.     # - the second piece is the length in seconds (a day is 60s * 60m * 24h)
  374.     parts = [(suffixes[0], 60 * 60 * 24 * 7 * 52),
  375.              (suffixes[1], 60 * 60 * 24 * 7),
  376.              (suffixes[2], 60 * 60 * 24),
  377.              (suffixes[3], 60 * 60),
  378.              (suffixes[4], 60),
  379.              (suffixes[5], 1)]
  380.  
  381.     # for each time piece, grab the value and remaining seconds, and add it to
  382.     # the time string
  383.     for suffix, length in parts:
  384.         value = seconds // length
  385.         if value > 0:
  386.             seconds = seconds % length
  387.             time.append('%s%s' % (str(value),
  388.                                   (suffix, (suffix, suffix + 's')[value > 1])[add_s]))
  389.         if seconds < 1:
  390.             break
  391.  
  392.     return separator.join(time)
  393.  
  394.  
  395. def splitBlock(dat):
  396.     '''split a amber parm dat file in blocks
  397.        0 = mass, 1 = extra + bond, 2 = angle, 3 = dihedral, 4 = improp, 5 = hbond
  398.        6 = equiv nbon, 7 = nbon, 8 = END, 9 = etc.
  399.     '''
  400.     dict_ = {}
  401.     count = 0
  402.     for line in dat:
  403.         line = line.rstrip()
  404.         if count in dict_:
  405.             dict_[count].append(line)
  406.         else:
  407.             dict_[count] = [line]
  408.         if not line:
  409.             count += 1
  410.     return dict_
  411.  
  412.  
  413. def getParCode(line):
  414.     key = line.replace(' -', '-').replace('- ', '-').split()[0]
  415.     return key
  416.  
  417.  
  418. def parseFrcmod(lista):
  419.     heads = ['MASS', 'BOND', 'ANGL', 'DIHE', 'IMPR', 'HBON', 'NONB']
  420.     dict_ = {}
  421.     for line in lista[1:]:
  422.         line = line.strip()
  423.         if line[:4] in heads:
  424.             head = line[:4]
  425.             dict_[head] = []
  426.             dd = {}
  427.             continue
  428.         elif line:
  429.             key = line.replace(' -', '-').replace('- ', '-').split()[0]
  430.             if key in dd:
  431.                 if not dd[key].count(line):
  432.                     dd[key].append(line)
  433.             else:
  434.                 dd[key] = [line]
  435.             dict_[head] = dd
  436.     for k in dict_.keys():
  437.         if not dict_[k]:
  438.             dict_.pop(k)
  439.     return dict_
  440.  
  441.  
  442. def parmMerge(fdat1, fdat2, frcmod=False):
  443.     '''merge two amber parm dat/frcmod files and save in /tmp'''
  444.     name1 = os.path.basename(fdat1).split('.dat')[0]
  445.     if frcmod:
  446.         name2 = os.path.basename(fdat2).split('.')[1]
  447.     else:
  448.         name2 = os.path.basename(fdat2).split('.dat')[0]
  449.     mname = '/tmp/' + name1 + name2 + '.dat'
  450.     mdatFile = open(mname, 'w')
  451.     mdat = ['merged %s %s' % (name1, name2)]
  452.  
  453.     # if os.path.exists(mname): return mname
  454.     dat1 = splitBlock(open(fdat1).readlines())
  455.  
  456.     if frcmod:
  457.         dHeads = {'MASS': 0, 'BOND': 1, 'ANGL': 2, 'DIHE': 3, 'IMPR': 4, 'HBON': 5, 'NONB': 7}
  458.         dat2 = parseFrcmod(open(fdat2).readlines())  # dict
  459.         for k in dat2:
  460.             for parEntry in dat2[k]:
  461.                 idFirst = None
  462.                 for line in dat1[dHeads[k]][:]:
  463.                     if line:
  464.                         key = line.replace(' -', '-').replace('- ', '-').split()[0]
  465.                         if key == parEntry:
  466.                             if not idFirst:
  467.                                 idFirst = dat1[dHeads[k]].index(line)
  468.                             dat1[dHeads[k]].remove(line)
  469.                 rev = dat2[k][parEntry][:]
  470.                 rev.reverse()
  471.                 if idFirst is None:
  472.                     idFirst = 0
  473.                 for ll in rev:
  474.                     if dHeads[k] in [0, 1, 7]:  # MASS has title in index 0 and so BOND, NONB
  475.                         dat1[dHeads[k]].insert(idFirst + 1, ll)
  476.                     else:
  477.                         dat1[dHeads[k]].insert(idFirst, ll)
  478.         dat1[0][0] = mdat[0]
  479.         for k in dat1:
  480.             for line in dat1[k]:
  481.                 mdatFile.write(line + '\n')
  482.         return mname
  483.  
  484.     dat2 = splitBlock(open(fdat2).readlines())
  485.     for k in dat1.keys()[:8]:
  486.         if k == 0:
  487.             lines = dat1[k][1:-1] + dat2[k][1:-1] + ['']
  488.             for line in lines:
  489.                 mdat.append(line)
  490.         if k == 1:
  491.             for i in dat1[k]:
  492.                 if '-' in i:
  493.                     id1 = dat1[k].index(i)
  494.                     break
  495.             for j in dat2[k]:
  496.                 if '-' in j:
  497.                     id2 = dat2[k].index(j)
  498.                     break
  499.             l1 = dat1[k][:id1]
  500.             l2 = dat2[k][:id2]
  501.             line = ''
  502.             for item in l1 + l2:
  503.                 line += item.strip() + ' '
  504.             mdat.append(line)
  505.             lines = dat1[k][id1:-1] + dat2[k][id2:-1] + ['']
  506.             for line in lines:
  507.                 mdat.append(line)
  508.         if k in [2, 3, 4, 5, 6]:  # angles, p dih, imp dih
  509.             lines = dat1[k][:-1] + dat2[k][:-1] + ['']
  510.             for line in lines:
  511.                 mdat.append(line)
  512.         if k == 7:
  513.             lines = dat1[k][:-1] + dat2[k][1:-1] + ['']
  514.             for line in lines:
  515.                 mdat.append(line)
  516.     for k in dat1.keys()[8:]:
  517.         for line in dat1[k]:
  518.             mdat.append(line)
  519.     for k in dat2.keys()[9:]:
  520.         for line in dat2[k]:
  521.             mdat.append(line)
  522.     for line in mdat:
  523.         mdatFile.write(line + '\n')
  524.     mdatFile.close()
  525.  
  526.     return mname
  527.  
  528.  
  529. def _getoutput(cmd):
  530.     '''to simulate commands.getoutput in order to work with python 2.6 up to 3.x'''
  531.     out = sub.Popen(cmd, shell=True, stderr=sub.STDOUT, stdout=sub.PIPE).communicate()[0][:-1]
  532.     try:
  533.         o = str(out.decode())
  534.     except:
  535.         o = str(out)
  536.     return o
  537.  
  538.  
  539. class AbstractTopol(object):
  540.  
  541.     """
  542.         Super class to build topologies
  543.     """
  544.  
  545.     def __init__(self):
  546.         if self.__class__ is AbstractTopol:
  547.             raise TypeError("Attempt to create istance of abstract class AbstractTopol")
  548.  
  549.     def printDebug(self, text=''):
  550.         if self.debug:
  551.             print('DEBUG: %s' % text)
  552.  
  553.     def printWarn(self, text=''):
  554.         if self.verbose:
  555.             print('WARNING: %s' % text)
  556.  
  557.     def printError(self, text=''):
  558.         if self.verbose:
  559.             print('ERROR: %s' % text)
  560.  
  561.     def printMess(self, text=''):
  562.         if self.verbose:
  563.             print('==> %s' % text)
  564.  
  565.     def printQuoted(self, text=''):
  566.         if self.verbose:
  567.             print(10 * '+' + 'start_quote' + 59 * '+')
  568.             print(text)
  569.             print(10 * '+' + 'end_quote' + 61 * '+')
  570.  
  571.     def guessCharge(self):
  572.         """
  573.             Guess the charge of a system based on antechamber
  574.             Returns None in case of error
  575.         """
  576.         done = False
  577.         error = False
  578.         charge = self.chargeVal
  579.         localDir = os.path.abspath('.')
  580.         if not os.path.exists(self.tmpDir):
  581.             os.mkdir(self.tmpDir)
  582.         if not os.path.exists(os.path.join(self.tmpDir, self.inputFile)):
  583.             copy2(self.absInputFile, self.tmpDir)
  584.         os.chdir(self.tmpDir)
  585.  
  586.         if self.chargeType == 'user':
  587.             if self.ext == '.mol2':
  588.                 self.printMess("Reading user's charges from mol2 file...")
  589.                 charge = self.readMol2TotalCharge(self.inputFile)
  590.                 done = True
  591.             else:
  592.                 self.printWarn("cannot read charges from a PDB file")
  593.                 self.printWarn("using now 'bcc' method for charge")
  594.  
  595.         if self.chargeVal is None and not done:
  596.             self.printWarn("no charge value given, trying to guess one...")
  597.             mol2FileForGuessCharge = self.inputFile
  598.             if self.ext == ".pdb":
  599.                 cmd = '%s -ipdb %s -omol2 %s.mol2' % (self.babelExe, self.inputFile,
  600.                                                       self.baseName)
  601.                 self.printDebug("guessCharge: " + cmd)
  602.                 out = _getoutput(cmd)
  603.                 self.printDebug(out)
  604.                 mol2FileForGuessCharge = os.path.abspath(self.baseName + ".mol2")
  605.                 in_mol = 'mol2'
  606.             else:
  607.                 in_mol = self.ext[1:]
  608.                 if in_mol == 'mol':
  609.                     in_mol = 'mdl'
  610.  
  611.             cmd = '%s -i %s -fi %s -o tmp -fo mol2 -c gas -pf y' % \
  612.                 (self.acExe, mol2FileForGuessCharge, in_mol)
  613.  
  614.             if self.debug:
  615.                 self.printMess("Debugging...")
  616.                 cmd = cmd.replace('-pf y', '-pf n')
  617.                 print(cmd)
  618.  
  619.             log = _getoutput(cmd).strip()
  620.  
  621.             if os.path.exists('tmp'):
  622.                 charge = self.readMol2TotalCharge('tmp')
  623.             else:
  624.                 try:
  625.                     charge = float(log.strip().split('equal to the total charge (')[-1].split(') based on Gasteiger atom type, exit')[0])
  626.                 except:
  627.                     error = True
  628.  
  629.             if error:
  630.                 self.printError("guessCharge failed")
  631.                 os.chdir(localDir)
  632.                 self.printQuoted(log)
  633.                 self.printMess("Trying with net charge = 0 ...")
  634. #                self.chargeVal = 0
  635.                 return None
  636.         charge = float(charge)
  637.         charge2 = int(round(charge))
  638.         drift = abs(charge2 - charge)
  639.         self.printDebug("Net charge drift '%3.6f'" % drift)
  640.         if drift > diffTol:
  641.             self.printError("Net charge drift '%3.5f' bigger than tolerance '%3.5f'" % (drift, diffTol))
  642.             if not self.force:
  643.                 sys.exit(1)
  644.         self.chargeVal = str(charge2)
  645.         self.printMess("... charge set to %i" % charge2)
  646.         os.chdir(localDir)
  647.  
  648.     def setResNameCheckCoords(self):
  649.         """Set a 3 letter residue name
  650.            and check coords duplication
  651.         """
  652.         exit_ = False
  653.         localDir = os.path.abspath('.')
  654.         if not os.path.exists(self.tmpDir):
  655.             os.mkdir(self.tmpDir)
  656.         # if not os.path.exists(os.path.join(tmpDir, self.inputFile)):
  657.         copy2(self.absInputFile, self.tmpDir)
  658.         os.chdir(self.tmpDir)
  659.  
  660.         exten = self.ext[1:]
  661.         if self.ext == '.pdb':
  662.             tmpFile = open(self.inputFile, 'r')
  663.         else:
  664.             if exten == 'mol':
  665.                 exten = 'mdl'
  666.             cmd = '%s -i %s -fi %s -o tmp -fo ac -pf y' % \
  667.                 (self.acExe, self.inputFile, exten)
  668.             self.printDebug(cmd)
  669.             out = _getoutput(cmd)
  670.             if not out.isspace():
  671.                 self.printDebug(out)
  672.             try:
  673.                 tmpFile = open('tmp', 'r')
  674.             except:
  675.                 rmtree(self.tmpDir)
  676.                 raise
  677.  
  678.         tmpData = tmpFile.readlines()
  679.         residues = set()
  680.         coords = {}
  681.         for line in tmpData:
  682.             if 'ATOM  ' in line or 'HETATM' in line:
  683.                 residues.add(line[17:20])
  684.                 at = line[0:17]
  685.                 cs = line[30:54]
  686.                 if cs in coords:
  687.                     coords[cs].append(at)
  688.                 else:
  689.                     coords[cs] = [at]
  690.         # self.printDebug(coords)
  691.  
  692.         if len(residues) > 1:
  693.             self.printError("more than one residue detected '%s'" % str(residues))
  694.             self.printError("verify your input file '%s'. Aborting ..." % self.inputFile)
  695.             sys.exit(1)
  696.  
  697.         dups = ""
  698.         shortd = ""
  699.         longd = ""
  700.         longSet = set()
  701.         id_ = 0
  702.         items = list(coords.items())
  703.         l = len(items)
  704.         for item in items:
  705.             id_ += 1
  706.             if len(item[1]) > 1:  # if True means atoms with same coordinates
  707.                 for i in item[1]:
  708.                     dups += "%s %s\n" % (i, item[0])
  709.  
  710. #        for i in xrange(0,len(data),f):
  711. #            fdata += (data[i:i+f])+' '
  712.  
  713.             for id2 in range(id_, l):
  714.                 item2 = items[id2]
  715.                 c1 = list(map(float, [item[0][i:i + 8] for i in range(0, 24, 8)]))
  716.                 c2 = list(map(float, [item2[0][i:i + 8] for i in range(0, 24, 8)]))
  717.                 dist2 = self.distance(c1, c2)
  718.                 if dist2 < minDist2:
  719.                     dist = math.sqrt(dist2)
  720.                     shortd += "%8.5f       %s %s\n" % (dist, item[1], item2[1])
  721.                 if dist2 < maxDist2:  # and not longOK:
  722.                     longSet.add(str(item[1]))
  723.                     longSet.add(str(item2[1]))
  724.             if str(item[1]) not in longSet and l > 1:
  725.                 longd += "%s\n" % item[1]
  726.  
  727.         if dups:
  728.             self.printError("Atoms with same coordinates in '%s'!" % self.inputFile)
  729.             self.printQuoted(dups[:-1])
  730.             exit_ = True
  731.  
  732.         if shortd:
  733.             self.printError("Atoms TOO close (< %s Ang.)" % minDist)
  734.             self.printQuoted("Dist (Ang.)    Atoms\n" + shortd[:-1])
  735.             exit_ = True
  736.  
  737.         if longd:
  738.             self.printError("Atoms TOO alone (> %s Ang.)" % maxDist)
  739.             self.printQuoted(longd[:-1])
  740.             exit_ = True
  741.  
  742.         if exit_:
  743.             if self.force:
  744.                 self.printWarn("You chose to proceed anyway with '-f' option. GOOD LUCK!")
  745.             else:
  746.                 self.printError("Use '-f' option if you want to proceed anyway. Aborting ...")
  747.                 rmtree(self.tmpDir)
  748.                 sys.exit(1)
  749.  
  750.         resname = list(residues)[0].strip()
  751.         newresname = resname
  752.  
  753.         # To avoid resname likes: 001 (all numbers), 1e2 (sci number), ADD : reserved terms for leap
  754.         leapWords = ['_cmd_options_', '_types_', 'add', 'addAtomTypes', 'addIons',
  755.                      'addIons2', 'addPath', 'addPdbAtomMap', 'addPdbResMap', 'alias',
  756.                      'alignAxes', 'bond', 'bondByDistance', 'center', 'charge',
  757.                      'check', 'clearPdbAtomMap', 'clearPdbResMap', 'clearVariables',
  758.                      'combine', 'copy', 'createAtom', 'createParmset', 'createResidue',
  759.                      'createUnit', 'crossLink', 'debugOff', 'debugOn', 'debugStatus',
  760.                      'deleteBond', 'deleteOffLibEntry', 'deleteRestraint', 'desc',
  761.                      'deSelect', 'displayPdbAtomMap', 'displayPdbResMap', 'edit',
  762.                      'flip', 'groupSelectedAtoms', 'help', 'impose', 'list', 'listOff',
  763.                      'loadAmberParams', 'loadAmberPrep', 'loadMol2', 'loadOff',
  764.                      'loadPdb', 'loadPdbUsingSeq', 'logFile', 'matchVariables',
  765.                      'measureGeom', 'quit', 'relax', 'remove', 'restrainAngle',
  766.                      'restrainBond', 'restrainTorsion', 'saveAmberParm',
  767.                      'saveAmberParmPert', 'saveAmberParmPol', 'saveAmberParmPolPert',
  768.                      'saveAmberPrep', 'saveMol2', 'saveOff', 'saveOffParm', 'savePdb',
  769.                      'scaleCharges', 'select', 'sequence', 'set', 'setBox', 'solvateBox',
  770.                      'solvateCap', 'solvateDontClip', 'solvateOct', 'solvateShell',
  771.                      'source', 'transform', 'translate', 'verbosity', 'zMatrix']
  772.         isLeapWord = False
  773.         for word in leapWords:
  774.             if resname.upper().startswith(word.upper()):
  775.                 self.printDebug("Residue name is a reserved word: '%s'" % word.upper())
  776.                 isLeapWord = True
  777.         try:
  778.             float(resname)
  779.             self.printDebug("Residue name is a 'number': '%s'" % resname)
  780.             isNumber = True
  781.         except ValueError:
  782.             isNumber = False
  783.  
  784.         if resname[0].isdigit() or isNumber or isLeapWord:
  785.             newresname = 'R' + resname
  786.         if not resname.isalnum():
  787.             newresname = 'MOL'
  788.         if newresname != resname:
  789.             self.printWarn("In %s.lib, residue name will be '%s' instead of '%s' elsewhere"
  790.                            % (self.acBaseName, newresname, resname))
  791.  
  792.         self.resName = newresname
  793.  
  794.         os.chdir(localDir)
  795.         self.printDebug("setResNameCheckCoords done")
  796.  
  797.     def distance(self, c1, c2):
  798.         # print c1, c2
  799.         dist2 = (c1[0] - c2[0]) ** 2 + (c1[1] - c2[1]) ** 2 + (c1[0] - c2[0]) ** 2 + \
  800.                 (c1[2] - c2[2]) ** 2
  801.         # dist2 = math.sqrt(dist2)
  802.         return dist2
  803.  
  804.     def readMol2TotalCharge(self, mol2File):
  805.         """Reads the charges in given mol2 file and returns the total
  806.         """
  807.         charge = 0.0
  808.         ll = []
  809.         cmd = '%s -i %s -fi mol2 -o tmp -fo mol2 -c wc -cf tmp.crg -pf y' % \
  810.             (self.acExe, mol2File)
  811.         if self.debug:
  812.             self.printMess("Debugging...")
  813.             cmd = cmd.replace('-pf y', '-pf n')
  814.  
  815.         self.printDebug(cmd)
  816.  
  817.         log = _getoutput(cmd)
  818.  
  819.         if os.path.exists('tmp.crg'):
  820.             tmpFile = open('tmp.crg', 'r')
  821.             tmpData = tmpFile.readlines()
  822.             for line in tmpData:
  823.                 ll += line.split()
  824.             charge = sum(map(float, ll))
  825.         if not log.isspace() and self.debug:
  826.             self.printQuoted(log)
  827.  
  828.         self.printDebug("readMol2TotalCharge: " + str(charge))
  829.  
  830.         return charge
  831.  
  832.     def execAntechamber(self, chargeType=None, atomType=None):
  833.         """ AmaberTools 1.3
  834.             To call Antechamber and execute it
  836. Usage: antechamber -i  input file name
  837.                    -fi input file format
  838.                    -o  output file name
  839.                    -fo output file format
  840.                    -c  charge method
  841.                    -cf charge file name
  842.                    -nc net molecular charge (int)
  843.                    -a  additional file name
  844.                    -fa additional file format
  845.                    -ao additional file operation
  846.                         crd  only read in coordinate
  847.                         crg only read in charge
  848.                         name   only read in atom name
  849.                         type   only read in atom type
  850.                         bond   only read in bond type
  851.                    -m  multiplicity (2S+1), default is 1
  852.                    -rn residue name, overrides input file, default is MOL
  853.                    -rf residue toplogy file name in prep input file,
  854.                               default is molecule.res
  855.                    -ch check file name for gaussian, default is molecule
  856.                    -ek mopac or sqm keyword, inside quotes
  857.                    -gk gaussian keyword, inside quotes
  858.                    -df am1-bcc flag, 2 - use sqm(default); 0 - use mopac
  859.                    -at atom type, can be gaff (default), amber, bcc and sybyl
  860.                    -du fix duplicate atom names: yes(y)[default] or no(n)
  861.                    -j  atom type and bond type prediction index, default is 4
  862.                         0     no assignment
  863.                         1     atom type
  864.                         2     full  bond types
  865.                         3     part  bond types
  866.                         4     atom and full bond type
  867.                         5     atom and part bond type
  868.                    -s  status information: 0(brief), 1(default) or 2(verbose)
  869.                    -pf remove intermediate files: yes(y) or no(n)[default]
  870.                    -i -o -fi and -fo must appear; others are optional
  872.                       List of the File Formats
  874.          file format type  abbre. index | file format type abbre. index
  875.         ---------------------------------------------------------------
  876.         Antechamber        ac       1  | Sybyl Mol2         mol2    2
  877.         PDB                pdb      3  | Modified PDB       mpdb    4
  878.         AMBER PREP (int)   prepi    5  | AMBER PREP (car)   prepc   6
  879.         Gaussian Z-Matrix  gzmat    7  | Gaussian Cartesian gcrt    8
  880.         Mopac Internal     mopint   9  | Mopac Cartesian    mopcrt 10
  881.         Gaussian Output    gout    11  | Mopac Output       mopout 12
  882.         Alchemy            alc     13  | CSD                csd    14
  883.         MDL                mdl     15  | Hyper              hin    16
  884.         AMBER Restart      rst     17  | Jaguar Cartesian   jcrt   18
  885.         Jaguar Z-Matrix    jzmat   19  | Jaguar Output      jout   20
  886.         Divcon Input       divcrt  21  | Divcon Output      divout 22
  887.         Charmm             charmm  23  | SQM Output         sqmout 24
  888.         --------------------------------------------------------------
  890.                 AMBER restart file can only be read in as additional file.
  892.                       List of the Charge Methods
  894.         charge method     abbre.  index | charge method    abbre. index
  895.         ----------------------------------------------------------------
  896.         RESP               resp     1  |  AM1-BCC            bcc     2
  897.         CM1                cm1      3  |  CM2                cm2     4
  898.         ESP (Kollman)      esp      5  |  Mulliken           mul     6
  899.         Gasteiger          gas      7  |  Read in charge     rc      8
  900.         Write out charge   wc       9  |  Delete Charge      dc     10
  901.         ----------------------------------------------------------------
  902. a        """
  903.         global pid
  904.  
  905.         self.printMess("Executing Antechamber...")
  906.  
  907.         self.makeDir()
  908.  
  909.         ct = chargeType or self.chargeType
  910.         at = atomType or self.atomType
  911.         if 'amber2' in at:
  912.             at = 'amber'
  913.  
  914.         if ct == 'user':
  915.             ct = ''
  916.         else:
  917.             ct = '-c %s' % ct
  918.  
  919.         exten = self.ext[1:]
  920.         if exten == 'mol':
  921.             exten = 'mdl'
  922.  
  923.         cmd = '%s -i %s -fi %s -o %s -fo mol2 %s -nc %s -m %s -s 2 -df %i -at\
  924.  %s -pf y %s' % (self.acExe, self.inputFile, exten, self.acMol2FileName,
  925.                  ct, self.chargeVal, self.multiplicity, self.qFlag, at,
  926.                  self.ekFlag)
  927.  
  928.         if self.debug:
  929.             self.printMess("Debugging...")
  930.             cmd = cmd.replace('-pf y', '-pf n')
  931.  
  932.         self.printDebug(cmd)
  933.  
  934.         if os.path.exists(self.acMol2FileName) and not self.force:
  935.             self.printMess("AC output file present... doing nothing")
  936.         else:
  937.             try:
  938.                 os.remove(self.acMol2FileName)
  939.             except:
  940.                 pass
  941.             signal.signal(signal.SIGALRM, self.signal_handler)
  942.             signal.alarm(self.timeTol)
  943.             p = sub.Popen(cmd, shell=True, stderr=sub.STDOUT, stdout=sub.PIPE)
  944.             pid = p.pid
  945.  
  946.             out = str(p.communicate()[0].decode())  # p.stdout.read()
  947.             self.acLog = out
  948.  
  949.         if os.path.exists(self.acMol2FileName):
  950.             self.printMess("* Antechamber OK *")
  951.         else:
  952.             self.printQuoted(self.acLog)
  953.             return True
  954.  
  955.     def signal_handler(self, _signum, _frame):  # , pid = 0):
  956.         global pid
  957.         pids = self.job_pids_family(pid)
  958.         self.printDebug("PID: %s, PIDS: %s" % (pid, pids))
  959.         self.printMess("Timed out! Process %s killed, max exec time (%ss) exceeded"
  960.                        % (pids, self.timeTol))
  961.         # os.system('kill -15 %s' % pids)
  962.         for i in pids.split():
  963.             os.kill(int(i), 15)
  964.         raise Exception("Semi-QM taking too long to finish... aborting!")
  965.  
  966.     def job_pids_family(self, jpid):
  967.         '''INTERNAL: Return all job processes (PIDs)'''
  968.         pid = repr(jpid)
  969.         dict_pids = {}
  970.         pids = [pid]
  971.         cmd = "ps -A -o uid,pid,ppid|grep %i" % os.getuid()
  972.         out = _getoutput(cmd).split('\n')  # getoutput("ps -A -o uid,pid,ppid|grep %i" % os.getuid()).split('\n')
  973.         for item in out:
  974.             vec = item.split()
  975.             dict_pids[vec[2]] = vec[1]
  976.         while 1:
  977.             try:
  978.                 pid = dict_pids[pid]
  979.                 pids.append(pid)
  980.             except KeyError:
  981.                 break
  982.         return ' '.join(pids)
  983.  
  984.     def delOutputFiles(self):
  985.         delFiles = ['mopac.in', 'tleap.in', 'sleap.in', 'fixbo.log',
  986.                     'addhs.log', 'ac_tmp_ot.mol2', 'frcmod.ac_tmp', 'fragment.mol2',
  987.                     self.tmpDir]  # , 'divcon.pdb', 'mopac.pdb', 'mopac.out'] #'leap.log'
  988.         self.printMess("Removing temporary files...")
  989.         for file_ in delFiles:
  990.             file_ = os.path.join(self.absHomeDir, file_)
  991.             if os.path.exists(file_):
  992.                 if os.path.isdir(file_):
  993.                     rmtree(file_)
  994.                 else:
  995.                     os.remove(file_)
  996.  
  997.     def checkXyzAndTopFiles(self):
  998.         fileXyz = self.acXyzFileName
  999.         fileTop = self.acTopFileName
  1000.         if os.path.exists(fileXyz) and os.path.exists(fileTop):
  1001.             # self.acXyz = fileXyz
  1002.             # self.acTop = fileTop
  1003.             return True
  1004.         return False
  1005.  
  1006.     def execSleap(self):
  1007.         global pid
  1008.         self.makeDir()
  1009.  
  1010.         if self.ext == '.mol2':
  1011.             self.printWarn("Sleap doesn't work with mol2 files yet...")
  1012.             return True
  1013.  
  1014.         if self.chargeType != 'bcc':
  1015.             self.printWarn("Sleap works only with bcc charge method")
  1016.             return True
  1017.  
  1018.         if self.atomType != 'gaff':
  1019.             self.printWarn("Sleap works only with gaff atom type")
  1020.             return True
  1021.  
  1022.         sleapScpt = SLEAP_TEMPLATE % self.acParDict
  1023.  
  1024.         fp = open('sleap.in', 'w')
  1025.         fp.write(sleapScpt)
  1026.         fp.close()
  1027.  
  1028.         cmd = '%s -f sleap.in' % self.sleapExe
  1029.  
  1030.         if self.checkXyzAndTopFiles() and not self.force:
  1031.             self.printMess("Topologies files already present... doing nothing")
  1032.         else:
  1033.             try:
  1034.                 os.remove(self.acTopFileName)
  1035.                 os.remove(self.acXyzFileName)
  1036.             except:
  1037.                 pass
  1038.             self.printMess("Executing Sleap...")
  1039.             self.printDebug(cmd)
  1040.  
  1041.             p = sub.Popen(cmd, shell=True, stderr=sub.STDOUT, stdout=sub.PIPE)
  1042.             pid = p.pid
  1043.             signal.signal(signal.SIGALRM, self.signal_handler)
  1044.             signal.alarm(self.timeTol)
  1045.  
  1046.             out = str(p.communicate()[0].decode())  # p.stdout.read()
  1047.  
  1048.             self.sleapLog = out
  1049.             self.checkLeapLog(self.sleapLog)
  1050.  
  1051.             if self.checkXyzAndTopFiles():
  1052.                 self.printMess(" * Sleap OK *")
  1053.             else:
  1054.                 self.printQuoted(self.sleapLog)
  1055.                 return True
  1056.  
  1057.     def execTleap(self):
  1058.  
  1059.         fail = False
  1060.  
  1061.         self.makeDir()
  1062.  
  1063.         if self.ext == ".pdb":
  1064.             self.printMess('... converting pdb input file to mol2 input file')
  1065.             if self.convertPdbToMol2():
  1066.                 self.printError("convertPdbToMol2 failed")
  1067.  
  1068.         # print self.chargeVal
  1069.  
  1070.         if self.execAntechamber():
  1071.             self.printError("Antechamber failed")
  1072.             fail = True
  1073.             # sys.exit(1)
  1074.  
  1075.         if self.execParmchk():
  1076.             self.printError("Parmchk failed")
  1077.             fail = True
  1078.             # sys.exit(1)
  1079.  
  1080.         if fail:
  1081.             return True
  1082.  
  1083.         tleapScpt = TLEAP_TEMPLATE % self.acParDict
  1084.  
  1085.         fp = open('tleap.in', 'w')
  1086.         fp.write(tleapScpt)
  1087.         fp.close()
  1088.  
  1089.         cmd = '%s -f tleap.in' % self.tleapExe
  1090.  
  1091.         if self.checkXyzAndTopFiles() and not self.force:
  1092.             self.printMess("Topologies files already present... doing nothing")
  1093.         else:
  1094.             try:
  1095.                 os.remove(self.acTopFileName)
  1096.                 os.remove(self.acXyzFileName)
  1097.             except:
  1098.                 pass
  1099.             self.printMess("Executing Tleap...")
  1100.             self.printDebug(cmd)
  1101.             self.tleapLog = _getoutput(cmd)
  1102.             self.checkLeapLog(self.tleapLog)
  1103.  
  1104.         if self.checkXyzAndTopFiles():
  1105.             self.printMess("* Tleap OK *")
  1106.         else:
  1107.             self.printQuoted(self.tleapLog)
  1108.             return True
  1109.  
  1110.     def checkLeapLog(self, log):
  1111.         log = log.splitlines(True)
  1112.         check = ''
  1113.         block = False
  1114.         for line in log:
  1115.             # print "*"+line+"*"
  1116.             if "Checking '" in line:
  1117.                 # check += line
  1118.                 block = True
  1119.             if "Checking Unit." in line:
  1120.                 block = False
  1121.             if block:
  1122.                 check += line
  1123.         self.printQuoted(check[:-1])
  1124.  
  1125.     def locateDat(self, aFile):
  1126.         '''locate a file pertinent to $AMBERHOME/dat/leap/parm/'''
  1127.         amberhome = os.environ.get('AMBERHOME')
  1128.         if amberhome:
  1129.             aFileF = os.path.join(amberhome, 'dat/leap/parm', aFile)
  1130.             if os.path.exists(aFileF):
  1131.                 return aFileF
  1132.         aFileF = os.path.join(os.path.dirname(self.acExe), '../dat/leap/parm', aFile)
  1133.         if os.path.exists(aFileF):
  1134.             return aFileF
  1135.         return None
  1136.  
  1137.     def execParmchk(self):
  1138.  
  1139.         self.makeDir()
  1140.         cmd = '%s -i %s -f mol2 -o %s' % (self.parmchkExe, self.acMol2FileName,
  1141.                                           self.acFrcmodFileName)
  1142.  
  1143.         if 'amber' in self.atomType:
  1144.             gaffFile = self.locateDat(self.gaffDatfile)
  1145.             parmfile = self.locateDat('parm10.dat')
  1146.             frcmodffxxSB = self.locateDat('frcmod.ff14SB')
  1147.             # frcmodparmbsc0 = self.locateDat('frcmod.parmbsc0')
  1148.             parmGaffFile = parmMerge(parmfile, gaffFile)
  1149.             parmGaffffxxSBFile = parmMerge(parmGaffFile, frcmodffxxSB, frcmod=True)
  1150.             # parm99gaffff99SBparmbsc0File = parmMerge(parm99gaffff99SBFile, frcmodparmbsc0, frcmod = True)
  1151.             # parm10file = self.locateDat('parm10.dat') # PARM99 + frcmod.ff99SB + frcmod.parmbsc0 in AmberTools 1.4
  1152.  
  1153.             cmd += ' -p %s' % parmGaffffxxSBFile  # Ignoring BSC0
  1154.         elif 'gaff2' in self.atomType:
  1155.             cmd += ' -s 2'
  1156.  
  1157.         self.parmchkLog = _getoutput(cmd)
  1158.  
  1159.         self.printDebug(cmd)
  1160.  
  1161.         if os.path.exists(self.acFrcmodFileName):
  1162.             check = self.checkFrcmod()
  1163.             if check:
  1164.                 self.printWarn("Couldn't determine all parameters:")
  1165.                 self.printMess("From file '%s'\n" % self.acFrcmodFileName + check)
  1166.             else:
  1167.                 self.printMess("* Parmchk OK *")
  1168.         else:
  1169.             self.printQuoted(self.parmchkLog)
  1170.             return True
  1171.  
  1172.     def checkFrcmod(self):
  1173.         check = ""
  1174.         frcmodContent = open(self.acFrcmodFileName, 'r').readlines()
  1175.         for line in frcmodContent:
  1176.             if "ATTN, need revision" in line:
  1177.                 check += line
  1178.         return check
  1179.  
  1180.     def convertPdbToMol2(self):
  1181.         if self.ext == '.pdb':
  1182.             if self.execBabel():
  1183.                 self.printError("convert pdb to mol2 via babel failed")
  1184.                 return True
  1185.  
  1186.     def execBabel(self):
  1187.  
  1188.         self.makeDir()
  1189.  
  1190.         cmd = '%s -ipdb %s -omol2 %s.mol2' % (self.babelExe, self.inputFile,
  1191.                                               self.baseName)
  1192.         self.printDebug(cmd)
  1193.         self.babelLog = _getoutput(cmd)
  1194.         self.ext = '.mol2'
  1195.         self.inputFile = self.baseName + self.ext
  1196.         self.acParDict['ext'] = 'mol2'
  1197.         if os.path.exists(self.inputFile):
  1198.             self.printMess("* Babel OK *")
  1199.         else:
  1200.             self.printQuoted(self.babelLog)
  1201.             return True
  1202.  
  1203.     def makeDir(self):
  1204.  
  1205.         os.chdir(self.rootDir)
  1206.         self.absHomeDir = os.path.abspath(self.homeDir)
  1207.         if not os.path.exists(self.homeDir):
  1208.             os.mkdir(self.homeDir)
  1209.         os.chdir(self.homeDir)
  1210.         copy2(self.absInputFile, '.')
  1211.  
  1212.         return True
  1213.  
  1214.     def createACTopol(self):
  1215.         """
  1216.             If successful, Amber Top and Xyz files will be generated
  1217.         """
  1218.         # sleap = False
  1219.         if self.engine == 'sleap':
  1220.             if self.execSleap():
  1221.                 self.printError("Sleap failed")
  1222.                 self.printMess("... trying Tleap")
  1223.                 if self.execTleap():
  1224.                     self.printError("Tleap failed")
  1225.         if self.engine == 'tleap':
  1226.             if self.execTleap():
  1227.                 self.printError("Tleap failed")
  1228.                 if self.extOld == '.pdb':
  1229.                     self.printMess("... trying Sleap")
  1230.                     self.ext = self.extOld
  1231.                     self.inputFile = self.baseName + self.ext
  1232.                     if self.execSleap():
  1233.                         self.printError("Sleap failed")
  1234.         if not self.debug:
  1235.             self.delOutputFiles()
  1236.  
  1237.     def createMolTopol(self):
  1238.         """
  1239.             Create molTop obj
  1240.         """
  1241.         self.topFileData = open(self.acTopFileName, 'r').readlines()
  1242.         self.molTopol = MolTopol(self, verbose=self.verbose, debug=self.debug,
  1243.                                  gmx4=self.gmx4, disam=self.disam, direct=self.direct,
  1244.                                  is_sorted=self.sorted, chiral=self.chiral)
  1245.         if self.outTopols:
  1246.             if 'cns' in self.outTopols:
  1247.                 self.molTopol.writeCnsTopolFiles()
  1248.             if 'gmx' in self.outTopols:
  1249.                 self.molTopol.writeGromacsTopolFiles()
  1250.             if 'charmm' in self.outTopols:
  1251.                 self.writeCharmmTopolFiles()
  1252.         self.pickleSave()
  1253.  
  1254.     def pickleSave(self):
  1255.         """
  1256.             To restore:
  1257.                 from acpype import *
  1258.                 #import cPickle as pickle
  1259.                 import pickle
  1260.                 o = pickle.load(open('DDD.pkl','rb'))
  1261.                 NB: It fails to restore with ipython in Mac (Linux OK)
  1262.         """
  1263.         pklFile = self.baseName + ".pkl"
  1264.         dumpFlag = False
  1265.         if not os.path.exists(pklFile):
  1266.             mess = "Writing pickle file %s" % pklFile
  1267.             dumpFlag = True
  1268.         elif self.force:
  1269.             mess = "Overwriting pickle file %s" % pklFile
  1270.             dumpFlag = True
  1271.         else:
  1272.             mess = "Pickle file %s already present... doing nothing" % pklFile
  1273.         self.printMess(mess)
  1274.         if dumpFlag:
  1275.             with open(pklFile, "wb") as f:  # for python 2.6 or higher
  1276.                 # f = open(pklFile, "wb")
  1277.                 if verList[0] == 3:
  1278.                     pickle.dump(self, f, protocol=2, fix_imports=True)
  1279.                 else:
  1280.                     pickle.dump(self, f, protocol=2)
  1281.  
  1282.     def getFlagData(self, flag):
  1283.         """
  1284.             For a given acFileTop flag, return a list of the data related
  1285.         """
  1286.         block = False
  1287.         tFlag = '%FLAG ' + flag
  1288.         data = ''
  1289.  
  1290.         if len(self.topFileData) == 0:
  1291.             raise Exception("PRMTOP file empty?")
  1292.  
  1293.         for rawLine in self.topFileData:
  1294.             if '%COMMENT' in rawLine:
  1295.                 continue
  1296.             line = rawLine.replace('\r', '').replace('\n', '')
  1297.             if tFlag in line:
  1298.                 block = True
  1299.                 continue
  1300.             if block and '%FLAG ' in line:
  1301.                 break
  1302.             if block:
  1303.                 if '%FORMAT' in line:
  1304.                     line = line.strip().strip('%FORMAT()').split('.')[0]
  1305.                     for c in line:
  1306.                         if c.isalpha():
  1307.                             f = int(line.split(c)[1])
  1308.                             break
  1309.                     continue
  1310.                 data += line
  1311.         # data need format
  1312.         sdata = [data[i:i + f].strip() for i in range(0, len(data), f)]
  1313.         if '+' and '.' in data and flag != 'RESIDUE_LABEL':  # it's a float
  1314.             ndata = list(map(float, sdata))
  1315.         elif flag != 'RESIDUE_LABEL':
  1316.             try:  # try if it's integer
  1317.                 ndata = list(map(int, sdata))
  1318.             except:  # it's string
  1319.                 ndata = sdata
  1320.         else:
  1321.             ndata = sdata
  1322.         if flag == 'AMBER_ATOM_TYPE':
  1323.             nn = []
  1324.             ll = set()
  1325.             prefixed = False
  1326.             for ii in ndata:
  1327.                 prefixed = True
  1328.                 if ii[0].isdigit():
  1329.                     ll.add(ii)
  1330.                     ii = 'A' + ii
  1331.                 nn.append(ii)
  1332.             if prefixed and ll:
  1333.                 self.printDebug("GMX does not like atomtype starting with Digit")
  1334.                 self.printDebug("prefixing AtomType %s with 'A'." % list(ll))
  1335.             ndata = nn
  1336.         return ndata  # a list
  1337.  
  1338.     def getResidueLabel(self):
  1339.         """
  1340.             Get a 3 capital letters code from acFileTop
  1341.             Returns a list.
  1342.         """
  1343.         residueLabel = self.getFlagData('RESIDUE_LABEL')
  1344.         residueLabel = list(map(str, residueLabel))
  1345.         if residueLabel[0] != residueLabel[0].upper():
  1346.             self.printWarn("residue label '%s' in '%s' is not all UPPERCASE" %
  1347.                            (residueLabel[0], self.inputFile))
  1348.             self.printWarn("this may raise problem with some applications like CNS")
  1349.         self.residueLabel = residueLabel
  1350.  
  1351.     def getCoords(self):
  1352.         """
  1353.             For a given acFileXyz file, return a list of coords as:
  1354.             [[x1,y1,z1],[x2,y2,z2], etc.]
  1355.         """
  1356.         if len(self.xyzFileData) == 0:
  1357.             raise Exception("INPCRD file empty?")
  1358.         data = ''
  1359.         for rawLine in self.xyzFileData[2:]:
  1360.             line = rawLine.replace('\r', '').replace('\n', '')
  1361.             data += line
  1362.         l = len(data)
  1363.         ndata = list(map(float, [data[i:i + 12] for i in range(0, l, 12)]))
  1364.  
  1365.         gdata = []
  1366.         for i in range(0, len(ndata), 3):
  1367.             gdata.append([ndata[i], ndata[i + 1], ndata[i + 2]])
  1368.  
  1369.         self.printDebug("getCoords done")
  1370.  
  1371.         return gdata
  1372.  
  1373.     def getAtoms(self):
  1374.         """
  1375.             Set a list with all atoms objects build from dat in acFileTop
  1376.             Set also if molTopol atom type system is gaff or amber
  1377.             Set also list atomTypes
  1378.             Set also resid
  1379.             Set also molTopol total charge
  1380.         """
  1381.         atomNameList = self.getFlagData('ATOM_NAME')
  1382.         atomTypeNameList = self.getFlagData('AMBER_ATOM_TYPE')
  1383.         self._atomTypeNameList = atomTypeNameList
  1384.         massList = self.getFlagData('MASS')
  1385.         chargeList = self.getFlagData('CHARGE')
  1386. #        totalCharge = sum(chargeList)
  1387. #        self.printDebug('charge to be balanced: total %13.10f' % (totalCharge/qConv))
  1388.  
  1389.         resIds = self.getFlagData('RESIDUE_POINTER') + [0]
  1390.         # to guess the resId of the last residue before ion or water
  1391. #        for resTemp in self.residueLabel:
  1392. #            if resTemp in ionOrWaterResNameList:
  1393. #                lastSoluteResId = self.residueLabel.index(resTemp) - 1
  1394. #                break
  1395.         # print lastSoluteResId, self.residueLabel[lastSoluteResId]
  1396.         # uniqAtomTypeId = self.getFlagData('ATOM_TYPE_INDEX') # for LJ
  1397. #        balanceChargeList = self.balanceCharges(chargeList)
  1398.         coords = self.getCoords()
  1399.         ACOEFs, BCOEFs = self.getABCOEFs()
  1400.  
  1401.         atoms = []
  1402.         atomTypes = []
  1403.         tmpList = []  # a list with unique atom types
  1404.         totalCharge = 0.0
  1405.         countRes = 0
  1406.         id_ = 0
  1407.         FirstNonSoluteId = None
  1408.         for atomName in atomNameList:
  1409.             if atomName != atomName.upper():
  1410.                 self.printDebug("atom name '%s' HAS to be all UPPERCASE... Applying this here." %
  1411.                                 atomName)
  1412.                 atomName = atomName.upper()
  1413.             atomTypeName = atomTypeNameList[id_]
  1414.             if id_ + 1 == resIds[countRes]:
  1415.                 resid = countRes  # self.residueLabel[countRes]
  1416.                 countRes += 1
  1417.             resName = self.residueLabel[resid]
  1418.             if resName in ionOrSolResNameList and not FirstNonSoluteId:
  1419.                 FirstNonSoluteId = id_
  1420.                 # print id_, resid, resName
  1421.             mass = massList[id_]
  1422.             # charge = balanceChargeList[id_]
  1423.             charge = chargeList[id_]
  1424.             chargeConverted = charge / qConv
  1425.             totalCharge += charge
  1426.             coord = coords[id_]
  1427.             ACOEF = ACOEFs[id_]
  1428.             BCOEF = BCOEFs[id_]
  1429.             atomType = AtomType(atomTypeName, mass, ACOEF, BCOEF)
  1430.             if atomTypeName not in tmpList:
  1431.                 tmpList.append(atomTypeName)
  1432.                 atomTypes.append(atomType)
  1433.             atom = Atom(atomName, atomType, id_ + 1, resid, mass, chargeConverted, coord)
  1434.             atoms.append(atom)
  1435.             id_ += 1
  1436.  
  1437.         balanceChargeList, balanceValue, balanceIds = self.balanceCharges(chargeList, FirstNonSoluteId)
  1438.  
  1439.         for id_ in balanceIds:
  1440.             atoms[id_].charge = balanceValue / qConv
  1441.         # self.printDebug("atom ids and balanced charges: %s, %3f10" % (balanceIds, balanceValue/qConv))
  1442.  
  1443.         if atomTypeName[0].islower():
  1444.             self.atomTypeSystem = 'gaff'
  1445.         else:
  1446.             self.atomTypeSystem = 'amber'
  1447.  
  1448.         self.printDebug('Balanced TotalCharge %13.10f' % float(sum(balanceChargeList) / qConv))
  1449.         self.totalCharge = int(totalCharge)
  1450.  
  1451.         self.atoms = atoms
  1452.         self.atomTypes = atomTypes
  1453.  
  1454.         self.pbc = None
  1455.         if len(coords) == len(atoms) + 2 or len(coords) == len(atoms) * 2 + 2:
  1456.             self.pbc = [coords[-2], coords[-1]]
  1457.         self.printDebug("PBC = %s" % self.pbc)
  1458.         self.printDebug("getAtoms done")
  1459.  
  1460.     def getBonds(self):
  1461.         uniqKbList = self.getFlagData('BOND_FORCE_CONSTANT')
  1462.         uniqReqList = self.getFlagData('BOND_EQUIL_VALUE')
  1463.         bondCodeHList = self.getFlagData('BONDS_INC_HYDROGEN')
  1464.         bondCodeNonHList = self.getFlagData('BONDS_WITHOUT_HYDROGEN')
  1465.         bondCodeList = bondCodeHList + bondCodeNonHList
  1466.         bonds = []
  1467.         for i in range(0, len(bondCodeList), 3):
  1468.             idAtom1 = bondCodeList[i] // 3  # remember python starts with id 0
  1469.             idAtom2 = bondCodeList[i + 1] // 3
  1470.             bondTypeId = bondCodeList[i + 2] - 1
  1471.             atom1 = self.atoms[idAtom1]
  1472.             atom2 = self.atoms[idAtom2]
  1473.             kb = uniqKbList[bondTypeId]
  1474.             req = uniqReqList[bondTypeId]
  1475.             atoms = [atom1, atom2]
  1476.             bond = Bond(atoms, kb, req)
  1477.             bonds.append(bond)
  1478.         self.bonds = bonds
  1479.         self.printDebug("getBonds done")
  1480.  
  1481.     def getAngles(self):
  1482.         uniqKtList = self.getFlagData('ANGLE_FORCE_CONSTANT')
  1483.         uniqTeqList = self.getFlagData('ANGLE_EQUIL_VALUE')
  1484.         # for list below, true atom number = index/3 + 1
  1485.         angleCodeHList = self.getFlagData('ANGLES_INC_HYDROGEN')
  1486.         angleCodeNonHList = self.getFlagData('ANGLES_WITHOUT_HYDROGEN')
  1487.         angleCodeList = angleCodeHList + angleCodeNonHList
  1488.         angles = []
  1489.         for i in range(0, len(angleCodeList), 4):
  1490.             idAtom1 = angleCodeList[i] // 3  # remember python starts with id 0
  1491.             idAtom2 = angleCodeList[i + 1] // 3
  1492.             idAtom3 = angleCodeList[i + 2] // 3
  1493.             angleTypeId = angleCodeList[i + 3] - 1
  1494.             atom1 = self.atoms[idAtom1]
  1495.             atom2 = self.atoms[idAtom2]
  1496.             atom3 = self.atoms[idAtom3]
  1497.             kt = uniqKtList[angleTypeId]
  1498.             teq = uniqTeqList[angleTypeId]  # angle given in rad in prmtop
  1499.             atoms = [atom1, atom2, atom3]
  1500.             angle = Angle(atoms, kt, teq)
  1501.             angles.append(angle)
  1502.         self.angles = angles
  1503.         self.printDebug("getAngles done")
  1504.  
  1505.     def getDihedrals(self):
  1506.         """
  1507.             Get dihedrals (proper and imp), condensed list of prop dih and
  1508.             atomPairs
  1509.         """
  1510.         uniqKpList = self.getFlagData('DIHEDRAL_FORCE_CONSTANT')
  1511.         uniqPeriodList = self.getFlagData('DIHEDRAL_PERIODICITY')
  1512.         uniqPhaseList = self.getFlagData('DIHEDRAL_PHASE')
  1513.         # for list below, true atom number = abs(index)/3 + 1
  1514.         dihCodeHList = self.getFlagData('DIHEDRALS_INC_HYDROGEN')
  1515.         dihCodeNonHList = self.getFlagData('DIHEDRALS_WITHOUT_HYDROGEN')
  1516.         dihCodeList = dihCodeHList + dihCodeNonHList
  1517.         properDih = []
  1518.         improperDih = []
  1519.         condProperDih = []  # list of dihedrals condensed by the same quartet
  1520.         # atomPairs = []
  1521.         atomPairs = set()
  1522.         for i in range(0, len(dihCodeList), 5):
  1523.             idAtom1 = dihCodeList[i] // 3  # remember python starts with id 0
  1524.             idAtom2 = dihCodeList[i + 1] // 3
  1525.             # 3 and 4 indexes can be negative: if id3 < 0, end group interations
  1526.             # in amber are to be ignored; if id4 < 0, dihedral is improper
  1527.             idAtom3raw = dihCodeList[i + 2] // 3  # can be negative -> exclude from 1-4vdw
  1528.             idAtom4raw = dihCodeList[i + 3] // 3  # can be negative -> Improper
  1529.             idAtom3 = abs(idAtom3raw)
  1530.             idAtom4 = abs(idAtom4raw)
  1531.             dihTypeId = dihCodeList[i + 4] - 1
  1532.             atom1 = self.atoms[idAtom1]
  1533.             atom2 = self.atoms[idAtom2]
  1534.             atom3 = self.atoms[idAtom3]
  1535.             atom4 = self.atoms[idAtom4]
  1536.             kPhi = uniqKpList[dihTypeId]  # already divided by IDIVF
  1537.             period = int(uniqPeriodList[dihTypeId])  # integer
  1538.             phase = uniqPhaseList[dihTypeId]  # angle given in rad in prmtop
  1539.             if phase == kPhi == 0:
  1540.                 period = 0  # period is set to 0
  1541.             atoms = [atom1, atom2, atom3, atom4]
  1542.             dihedral = Dihedral(atoms, kPhi, period, phase)
  1543.             if idAtom4raw > 0:
  1544.                 try:
  1545.                     atomsPrev = properDih[-1].atoms
  1546.                 except:
  1547.                     atomsPrev = []
  1548.                 properDih.append(dihedral)
  1549.                 if idAtom3raw < 0 and atomsPrev == atoms:
  1550.                     condProperDih[-1].append(dihedral)
  1551.                 else:
  1552.                     condProperDih.append([dihedral])
  1553.                 pair = (atom1, atom4)
  1554.                 # if atomPairs.count(pair) == 0 and idAtom3raw > 0:
  1555.                 if idAtom3raw > 0:
  1556.                     atomPairs.add(pair)
  1557.             else:
  1558.                 improperDih.append(dihedral)
  1559.         try:
  1560.             atomPairs = sorted(atomPairs)
  1561.         except:
  1562.             pass
  1563.         self.properDihedrals = properDih
  1564.         self.improperDihedrals = improperDih
  1565.         self.condensedProperDihedrals = condProperDih  # [[],[],...]
  1566.         self.atomPairs = atomPairs  # set((atom1, atom2), ...)
  1567.         self.printDebug("getDihedrals done")
  1568.  
  1569.     def getChirals(self):
  1570.         """
  1571.             Get chiral atoms, its 4 neighbours and improper dihedral angle
  1572.         """
  1573.         self.chiralGroups = []
  1574.         if self.obchiralExe:
  1575.             # print (self.obchiralExe, os.getcwd())
  1576.             cmd = '%s %s' % (self.obchiralExe, self.inputFile)
  1577.             # print(cmd)
  1578.             out = map(int, re.findall('Atom (\d+) Is', _getoutput(cmd)))
  1579.             # print("*%s*" % out)
  1580.             chiralGroups = []
  1581.             for id_ in out:
  1582.                 atChi = self.atoms[id_ - 1]
  1583.                 quad = []
  1584.                 for bb in self.bonds:
  1585.                     bAts = bb.atoms[:]
  1586.                     if atChi in bAts:
  1587.                         bAts.remove(atChi)
  1588.                         quad.append(bAts[0])
  1589.                 if len(quad) != 4:
  1590.                     if self.chiral:
  1591.                         self.printWarn("Atom %s has less than 4 connections to 4 different atoms. It's NOT Chiral!" % atChi)
  1592.                     continue
  1593.                 v1, v2, v3, v4 = [x.coords for x in quad]
  1594.                 chiralGroups.append((atChi, quad, imprDihAngle(v1, v2, v3, v4)))
  1595.             self.chiralGroups = chiralGroups
  1596.  
  1597.     def sortAtomsForGromacs(self):
  1598.         """
  1599.             Re-sort atoms for gromacs, which expects all hydrogens to immediately
  1600.             follow the heavy atom they are bonded to and belong to the same charge
  1601.             group.
  1603.             Currently, atom mass < 1.2 is taken to denote a proton.  This behavior
  1604.             may be changed by modifying the 'is_hydrogen' function within.
  1606.             JDC 2011-02-03
  1607.         """
  1608.  
  1609.         # Build dictionary of bonded atoms.
  1610.         bonded_atoms = dict()
  1611.         for atom in self.atoms:
  1612.             bonded_atoms[atom] = list()
  1613.         for bond in self.bonds:
  1614.             [atom1, atom2] = bond.atoms
  1615.             bonded_atoms[atom1].append(atom2)
  1616.             bonded_atoms[atom2].append(atom1)
  1617.  
  1618.         # Define hydrogen and heavy atom classes.
  1619.         def is_hydrogen(atom):
  1620.             return (atom.mass < 1.2)
  1621.  
  1622.         def is_heavy(atom):
  1623.             return not is_hydrogen(atom)
  1624.  
  1625.         # Build list of sorted atoms, assigning charge groups by heavy atom.
  1626.         sorted_atoms = list()
  1627.         cgnr = 1  # charge group number: each heavy atoms is assigned its own charge group
  1628.         # First pass: add heavy atoms, followed by the hydrogens bonded to them.
  1629.         for atom in self.atoms:
  1630.             if is_heavy(atom):
  1631.                 # Append heavy atom.
  1632.                 atom.cgnr = cgnr
  1633.                 sorted_atoms.append(atom)
  1634.                 # Append all hydrogens.
  1635.                 for bonded_atom in bonded_atoms[atom]:
  1636.                     if is_hydrogen(bonded_atom) and not (bonded_atom in sorted_atoms):
  1637.                         # Append bonded hydrogen.
  1638.                         bonded_atom.cgnr = cgnr
  1639.                         sorted_atoms.append(bonded_atom)
  1640.                 cgnr += 1
  1641.  
  1642.         # Second pass: Add any remaining atoms.
  1643.         if len(sorted_atoms) < len(self.atoms):
  1644.             for atom in self.atoms:
  1645.                 if not (atom in sorted_atoms):
  1646.                     atom.cgnr = cgnr
  1647.                     sorted_atoms.append(atom)
  1648.                     cgnr += 1
  1649.  
  1650.         # Replace current list of atoms with sorted list.
  1651.         self.atoms = sorted_atoms
  1652.  
  1653.         # Renumber atoms in sorted list, starting from 1.
  1654.         for (index, atom) in enumerate(self.atoms):
  1655.             atom.id = index + 1
  1656.  
  1657.         return
  1658.  
  1659.     def setAtomPairs(self):
  1660.         """
  1661.             Set a list of pair of atoms pertinent to interaction 1-4 for vdw.
  1662.             WRONG: Deprecated
  1663.         """
  1664.         atomPairs = []
  1665.         for item in self.condensedProperDihedrals:
  1666.             dih = item[0]
  1667.             atom1 = dih.atoms[0]
  1668.             atom2 = dih.atoms[3]
  1669.             pair = [atom1, atom2]
  1670.             if atomPairs.count(pair) == 0:
  1671.                 atomPairs.append(pair)
  1672.         self.atomPairs = atomPairs  # [[atom1, atom2], ...]
  1673.         self.printDebug("atomPairs done")
  1674.  
  1675.     def getExcludedAtoms(self):
  1676.         """
  1677.             Returns a list of atoms with a list of its excluded atoms up to 3rd
  1678.             neighbour.
  1679.             It's implicitly indexed, i.e., a sequence of atoms in position n in
  1680.             the excludedAtomsList corresponds to atom n (self.atoms) and so on.
  1681.             NOT USED
  1682.         """
  1683.         excludedAtomsIdList = self.getFlagData('EXCLUDED_ATOMS_LIST')
  1684.         numberExcludedAtoms = self.getFlagData('NUMBER_EXCLUDED_ATOMS')
  1685.         atoms = self.atoms
  1686.         interval = 0
  1687.         excludedAtomsList = []
  1688.         for number in numberExcludedAtoms:
  1689.             temp = excludedAtomsIdList[interval:interval + number]
  1690.             if temp == [0]:
  1691.                 excludedAtomsList.append([])
  1692.             else:
  1693.                 excludedAtomsList.append([atoms[a - 1] for a in temp])
  1694.             interval += number
  1695.         self.excludedAtoms = excludedAtomsList
  1696.         self.printDebug("getExcludedAtoms")
  1697.  
  1698.     def balanceCharges(self, chargeList, FirstNonSoluteId=None):
  1699.         """
  1700.             Note that python is very annoying about floating points.
  1701.             Even after balance, there will always be some residue of order e-12
  1702.             to e-16, which is believed to vanished once one writes a topology
  1703.             file, say, for CNS or GMX, where floats are represented with 4 or 5
  1704.             maximum decimals.
  1705.         """
  1706.         limIds = []
  1707.         # self.printDebug(chargeList)
  1708.         total = sum(chargeList)
  1709.         totalConverted = total / qConv
  1710.         self.printDebug('charge to be balanced: total %13.10f' % (totalConverted))
  1711.         maxVal = max(chargeList[:FirstNonSoluteId])
  1712.         minVal = min(chargeList[:FirstNonSoluteId])
  1713.         if abs(maxVal) >= abs(minVal):
  1714.             lim = maxVal
  1715.         else:
  1716.             lim = minVal
  1717.         nLims = chargeList.count(lim)
  1718.         # limId = chargeList.index(lim)
  1719.         diff = totalConverted - round(totalConverted)
  1720.         fix = lim - diff * qConv / nLims
  1721.         id_ = 0
  1722.         for c in chargeList:
  1723.             if c == lim:
  1724.                 limIds.append(id_)
  1725.                 chargeList[id_] = fix
  1726.             id_ += 1
  1727.         # self.printDebug(chargeList)
  1728.         self.printDebug("balanceCharges done")
  1729.         return chargeList, fix, limIds
  1730.  
  1731.     def getABCOEFs(self):
  1732.         uniqAtomTypeIdList = self.getFlagData('ATOM_TYPE_INDEX')
  1733.         nonBonIdList = self.getFlagData('NONBONDED_PARM_INDEX')
  1734.         rawACOEFs = self.getFlagData('LENNARD_JONES_ACOEF')
  1735.         rawBCOEFs = self.getFlagData('LENNARD_JONES_BCOEF')
  1736.         # print nonBonIdList, len(nonBonIdList), rawACOEFs, len(rawACOEFs)
  1737.         ACOEFs = []
  1738.         BCOEFs = []
  1739.         ntypes = max(uniqAtomTypeIdList)
  1740.         # id_ = 0
  1741.         # for atName in self._atomTypeNameList:
  1742.         for id_ in range(len(self._atomTypeNameList)):
  1743.             # id_ = self._atomTypeNameList.index(atName)
  1744.             atomTypeId = uniqAtomTypeIdList[id_]
  1745.             index = ntypes * (atomTypeId - 1) + atomTypeId
  1746.             nonBondId = nonBonIdList[index - 1]
  1747.             # print "*****", index, ntypes, atName, id_, atomTypeId, nonBondId
  1748.             ACOEFs.append(rawACOEFs[nonBondId - 1])
  1749.             BCOEFs.append(rawBCOEFs[nonBondId - 1])
  1750.             # id_ += 1
  1751.         # print ACOEFs
  1752.         self.printDebug("getABCOEFs done")
  1753.         return ACOEFs, BCOEFs
  1754.  
  1755.     def setProperDihedralsCoef(self):
  1756.         """
  1757.             It takes self.condensedProperDihedrals and returns
  1758.             self.properDihedralsCoefRB, a reduced list of quartet atoms + RB.
  1759.             Coeficients ready for GMX (multiplied by 4.184)
  1761.             self.properDihedralsCoefRB = [ [atom1,..., atom4], C[0:5] ]
  1763.             For proper dihedrals: a quartet of atoms may appear with more than
  1764.             one set of parameters and to convert to GMX they are treated as RBs.
  1766.             The resulting coefs calculated here may look slighted different from
  1767.             the ones calculated by amb2gmx.pl because python is taken full float
  1768.             number from prmtop and not rounded numbers from rdparm.out as
  1769.             amb2gmx.pl does.
  1770.         """
  1771.         properDihedralsCoefRB = []
  1772.         properDihedralsAlphaGamma = []
  1773.         properDihedralsGmx45 = []
  1774.         for item in self.condensedProperDihedrals:
  1775.             V = 6 * [0.0]
  1776.             C = 6 * [0.0]
  1777.             for dih in item:
  1778.                 period = dih.period  # Pn
  1779.                 kPhi = dih.kPhi  # in rad
  1780.                 phaseRaw = dih.phase * radPi  # in degree
  1781.                 phase = int(phaseRaw)  # in degree
  1782.                 if period > 4 and self.gmx4:
  1783.                     self.printError("Likely trying to convert ILDN to RB, DO NOT use option '-z'")
  1784.                     sys.exit(1)
  1785.                 if phase in [0, 180]:
  1786.                     properDihedralsGmx45.append([item[0].atoms, phaseRaw, kPhi, period])
  1787.                     if self.gmx4:
  1788.                         if kPhi > 0:
  1789.                             V[period] = 2 * kPhi * cal
  1790.                         if period == 1:
  1791.                             C[0] += 0.5 * V[period]
  1792.                             if phase == 0:
  1793.                                 C[1] -= 0.5 * V[period]
  1794.                             else:
  1795.                                 C[1] += 0.5 * V[period]
  1796.                         elif period == 2:
  1797.                             if phase == 180:
  1798.                                 C[0] += V[period]
  1799.                                 C[2] -= V[period]
  1800.                             else:
  1801.                                 C[2] += V[period]
  1802.                         elif period == 3:
  1803.                             C[0] += 0.5 * V[period]
  1804.                             if phase == 0:
  1805.                                 C[1] += 1.5 * V[period]
  1806.                                 C[3] -= 2 * V[period]
  1807.                             else:
  1808.                                 C[1] -= 1.5 * V[period]
  1809.                                 C[3] += 2 * V[period]
  1810.                         elif period == 4:
  1811.                             if phase == 180:
  1812.                                 C[2] += 4 * V[period]
  1813.                                 C[4] -= 4 * V[period]
  1814.                             else:
  1815.                                 C[0] += V[period]
  1816.                                 C[2] -= 4 * V[period]
  1817.                                 C[4] += 4 * V[period]
  1818.                 else:
  1819.                     properDihedralsAlphaGamma.append([item[0].atoms, phaseRaw, kPhi, period])
  1820.                     # print phaseRaw, kPhi, period
  1821.             if phase in [0, 180]:
  1822.                 properDihedralsCoefRB.append([item[0].atoms, C])
  1823.  
  1824.         # print properDihedralsCoefRB
  1825.         # print properDihedralsAlphaGamma
  1826.  
  1827.         self.printDebug("setProperDihedralsCoef done")
  1828.  
  1829.         self.properDihedralsCoefRB = properDihedralsCoefRB
  1830.         self.properDihedralsAlphaGamma = properDihedralsAlphaGamma
  1831.         self.properDihedralsGmx45 = properDihedralsGmx45
  1832.  
  1833.     def writeCharmmTopolFiles(self):
  1834.  
  1835.         self.printMess("Writing CHARMM files\n")
  1836.  
  1837.         # self.makeDir()
  1838.  
  1839.         at = self.atomType
  1840.         self.getResidueLabel()
  1841.         res = self.resName  # self.residueLabel[0]
  1842.         # print res, self.residueLabel, type(self.residueLabel)
  1843.  
  1844.         cmd = '%s -i %s -fi mol2 -o %s -fo charmm -s 2 -at %s \
  1845.         -pf y -rn %s' % (self.acExe, self.acMol2FileName, self.charmmBase,
  1846.                          at, res)
  1847.  
  1848.         if self.debug:
  1849.             cmd = cmd.replace('-pf y', '-pf n')
  1850.             self.printDebug(cmd)
  1851.         _log = _getoutput(cmd)
  1852.  
  1853.     def writePdb(self, file_):
  1854.         """
  1855.             Write a new PDB file_ with the atom names defined by Antechamber
  1856.             Input: file_ path string
  1857.             The format generated here use is slightly different from
  1858.             http://www.wwpdb.org/documentation/format23/sect9.html respected to
  1859.             atom name
  1860.         """
  1861.         # TODO: assuming only one residue ('1')
  1862.         pdbFile = open(file_, 'w')
  1863.         fbase = os.path.basename(file_)
  1864.         pdbFile.write("REMARK " + head % (fbase, date))
  1865.         id_ = 1
  1866.         for atom in self.atoms:
  1867.             # id_ = self.atoms.index(atom) + 1
  1868.             aName = atom.atomName
  1869.             if len(aName) == 2:
  1870.                 aName = ' %s ' % aName
  1871.             elif len(aName) == 1:
  1872.                 aName = ' %s  ' % aName
  1873.             for l in aName:
  1874.                 if l.isalpha():
  1875.                     s = l
  1876.                     break
  1877.             rName = self.residueLabel[0]
  1878.             x = atom.coords[0]
  1879.             y = atom.coords[1]
  1880.             z = atom.coords[2]
  1881.             line = "%-6s%5d %4s %3s Z%4d%s%8.3f%8.3f%8.3f%6.2f%6.2f%s%2s\n" % \
  1882.                 ('ATOM', id_, aName, rName, 1, 4 * ' ', x, y, z, 1.0, 0.0, 10 * ' ', s)
  1883.             pdbFile.write(line)
  1884.             id_ += 1
  1885.         pdbFile.write('END\n')
  1886.  
  1887.     def writeGromacsTopolFiles(self, amb2gmx=False):
  1888.         """
  1889.             # from ~/Programmes/amber10/dat/leap/parm/gaff.dat
  1890.             #atom type        atomic mass        atomic polarizability        comments
  1891.             ca                12.01                 0.360                    Sp2 C in pure aromatic systems
  1892.             ha                1.008                 0.135                    H bonded to aromatic carbon
  1894.             #bonded atoms        harmonic force kcal/mol/A^2       eq. dist. Ang.  comments
  1895.             ca-ha                  344.3*                           1.087**         SOURCE3  1496    0.0024    0.0045
  1896.             * for gmx: 344.3 * 4.184 * 100 * 2 = 288110 kJ/mol/nm^2 (why factor 2?)
  1897.             ** convert Ang to nm ( div by 10) for gmx: 1.087 A = 0.1087 nm
  1898.             # CA HA         1    0.10800   307105.6 ; ged from 340. bsd on C6H6 nmodes; PHE,TRP,TYR (from ffamber99bon.itp)
  1899.             # CA-HA  367.0    1.080       changed from 340. bsd on C6H6 nmodes; PHE,TRP,TYR (from parm99.dat)
  1901.             # angle        HF kcal/mol/rad^2    eq angle degrees     comments
  1902.             ca-ca-ha        48.5*             120.01                SOURCE3 2980   0.1509   0.2511
  1903.             * to convert to gmx: 48.5 * 4.184 * 2 = 405.848 kJ/mol/rad^2 (why factor 2?)
  1904.             # CA  CA  HA           1   120.000    418.400 ; new99 (from ffamber99bon.itp)
  1905.             # CA-CA-HA    50.0      120.00 (from parm99.dat)
  1907.             # dihedral    idivf        barrier hight/2 kcal/mol  phase degrees       periodicity     comments
  1908.             X -ca-ca-X    4           14.500*                     180.000                2.000             intrpol.bsd.on C6H6
  1909.             * to convert to gmx: 14.5/4 * 4.184 * 2 (?) (yes in amb2gmx, not in topolbuild, why?) = 30.334 or 15.167 kJ/mol
  1910.             # X -CA-CA-X    4   14.50        180.0             2.         intrpol.bsd.on C6H6 (from parm99.dat)
  1911.             # X   CA  CA  X     3    30.33400     0.00000   -30.33400     0.00000     0.00000     0.00000   ; intrpol.bsd.on C6H6
  1912.             ;propers treated as RBs in GROMACS to use combine multiple AMBER torsions per quartet (from ffamber99bon.itp)
  1914.             # impr. dihedral        barrier hight/2      phase degrees       periodicity     comments
  1915.             X -X -ca-ha             1.1*                  180.                      2.                   bsd.on C6H6 nmodes
  1916.             * to convert to gmx: 1.1 * 4.184 = 4.6024 kJ/mol/rad^2
  1917.             # X -X -CA-HA         1.1          180.          2.           bsd.on C6H6 nmodes (from parm99.dat)
  1918.             # X   X   CA  HA       1      180.00     4.60240     2      ; bsd.on C6H6 nmodes
  1919.             ;impropers treated as propers in GROMACS to use correct AMBER analytical function (from ffamber99bon.itp)
  1921.             # 6-12 parms     sigma = 2 * r * 2^(-1/6)    epsilon
  1922.             # atomtype        radius Ang.                    pot. well depth kcal/mol      comments
  1923.               ha                  1.4590*                      0.0150**                         Spellmeyer
  1924.               ca                  1.9080                    0.0860                            OPLS
  1925.             * to convert to gmx:
  1926.                 sigma = 1.4590 * 2^(-1/6) * 2 = 2 * 1.29982 Ang. = 2 * 0.129982 nm  = 1.4590 * 2^(5/6)/10 =  0.259964 nm
  1927.             ** to convert to gmx: 0.0150 * 4.184 = 0.06276 kJ/mol
  1928.             # amber99_3    CA     0.0000  0.0000  A   3.39967e-01  3.59824e-01 (from ffamber99nb.itp)
  1929.             # amber99_22   HA     0.0000  0.0000  A   2.59964e-01  6.27600e-02 (from ffamber99nb.itp)
  1930.             # C*          1.9080  0.0860             Spellmeyer
  1931.             # HA          1.4590  0.0150             Spellmeyer (from parm99.dat)
  1932.             # to convert r and epsilon to ACOEF and BCOEF
  1933.             # ACOEF = sqrt(e1*e2) * (r1 + r2)^12 ; BCOEF = 2 * sqrt(e1*e2) * (r1 + r2)^6 = 2 * ACOEF/(r1+r2)^6
  1934.             # to convert ACOEF and BCOEF to r and epsilon
  1935.             # r = 0.5 * (2*ACOEF/BCOEF)^(1/6); ep = BCOEF^2/(4*ACOEF)
  1936.             # to convert ACOEF and BCOEF to sigma and epsilon (GMX)
  1937.             # sigma = (ACOEF/BCOEF)^(1/6) * 0.1 ; epsilon = 4.184 * BCOEF^2/(4*ACOEF)
  1938.             #   ca   ca       819971.66        531.10
  1939.             #   ca   ha        76245.15        104.66
  1940.             #   ha   ha         5716.30         18.52
  1942.             For proper dihedrals: a quartet of atoms may appear with more than
  1943.             one set of parameters and to convert to GMX they are treated as RBs;
  1944.             use the algorithm:
  1945.               for(my $j=$i;$j<=$lines;$j++){
  1946.                 my $period = $pn{$j};
  1947.                 if($pk{$j}>0) {
  1948.                   $V[$period] = 2*$pk{$j}*$cal;
  1949.                 }
  1950.                 # assign V values to C values as predefined #
  1951.                 if($period==1){
  1952.                   $C[0]+=0.5*$V[$period];
  1953.                   if($phase{$j}==0){
  1954.                     $C[1]-=0.5*$V[$period];
  1955.                   }else{
  1956.                     $C[1]+=0.5*$V[$period];
  1957.                   }
  1958.                 }elsif($period==2){
  1959.                   if(($phase{$j}==180)||($phase{$j}==3.14)){
  1960.                     $C[0]+=$V[$period];
  1961.                     $C[2]-=$V[$period];
  1962.                   }else{
  1963.                     $C[2]+=$V[$period];
  1964.                   }
  1965.                 }elsif($period==3){
  1966.                   $C[0]+=0.5*$V[$period];
  1967.                   if($phase{$j}==0){
  1968.                     $C[1]+=1.5*$V[$period];
  1969.                     $C[3]-=2*$V[$period];
  1970.                   }else{
  1971.                     $C[1]-=1.5*$V[$period];
  1972.                     $C[3]+=2*$V[$period];
  1973.                   }
  1974.                 }elsif($period==4){
  1975.                   if(($phase{$j}==180)||($phase{$j}==3.14)){
  1976.                     $C[2]+=4*$V[$period];
  1977.                     $C[4]-=4*$V[$period];
  1978.                   }else{
  1979.                     $C[0]+=$V[$period];
  1980.                     $C[2]-=4*$V[$period];
  1981.                     $C[4]+=4*$V[$period];
  1982.                   }
  1983.                 }
  1984.               }
  1985.         """
  1986.  
  1987.         self.printMess("Writing GROMACS files\n")
  1988.  
  1989.         self.setAtomType4Gromacs()
  1990.  
  1991.         self.writeGroFile()
  1992.  
  1993.         self.writeGromacsTop(amb2gmx=amb2gmx)
  1994.  
  1995.         self.writeMdpFiles()
  1996.  
  1997.     def setAtomType4Gromacs(self):
  1998.         """Atom types names in Gromacs TOP file are not case sensitive;
  1999.            this routine will append a '_' to lower case atom type.
  2000.            E.g.: CA and ca -> CA and ca_
  2001.         """
  2002.         if self.disam:
  2003.             self.printMess("Disambiguating lower and uppercase atomtypes in GMX top file.\n")
  2004.             self.atomTypesGromacs = self.atomTypes
  2005.             self.atomsGromacs = self.atoms
  2006.             return
  2007.  
  2008.         atNames = [at.atomTypeName for at in self.atomTypes]
  2009.         # print atNames
  2010.         delAtomTypes = []
  2011.         modAtomTypes = []
  2012.         atomTypesGromacs = []
  2013.         dictAtomTypes = {}
  2014.         for at in self.atomTypes:
  2015.             atName = at.atomTypeName
  2016.             dictAtomTypes[atName] = at
  2017.             if atName.islower() and atName.upper() in atNames:
  2018.                 # print atName, atName.upper()
  2019.                 atUpper = self.atomTypes[atNames.index(atName.upper())]
  2020.                 # print at.atomTypeName,at.mass, at.ACOEF, at.BCOEF
  2021.                 # print atUpper.atomTypeName, atUpper.mass, atUpper.ACOEF, atUpper.BCOEF
  2022.                 if at.ACOEF is atUpper.ACOEF and at.BCOEF is at.BCOEF:
  2023.                     delAtomTypes.append(atName)
  2024.                 else:
  2025.                     newAtName = atName + '_'
  2026.                     modAtomTypes.append(atName)
  2027.                     atomType = AtomType(newAtName, at.mass, at.ACOEF, at.BCOEF)
  2028.                     atomTypesGromacs.append(atomType)
  2029.                     dictAtomTypes[newAtName] = atomType
  2030.             else:
  2031.                 atomTypesGromacs.append(at)
  2032.  
  2033.         atomsGromacs = []
  2034.         for a in self.atoms:
  2035.             atName = a.atomType.atomTypeName
  2036.             if atName in delAtomTypes:
  2037.                 atom = Atom(a.atomName, dictAtomTypes[atName.upper()], a.id,
  2038.                             a.resid, a.mass, a.charge, a.coords)
  2039.                 atom.cgnr = a.cgnr
  2040.                 atomsGromacs.append(atom)
  2041.             elif atName in modAtomTypes:
  2042.                 atom = Atom(a.atomName, dictAtomTypes[atName + '_'], a.id,
  2043.                             a.resid, a.mass, a.charge, a.coords)
  2044.                 atom.cgnr = a.cgnr
  2045.                 atomsGromacs.append(atom)
  2046.             else:
  2047.                 atomsGromacs.append(a)
  2048.  
  2049.         self.atomTypesGromacs = atomTypesGromacs
  2050.         self.atomsGromacs = atomsGromacs
  2051.         # print [i.atomTypeName for i in atomTypesGromacs]
  2052.         # print modAtomTypes
  2053.         # print delAtomTypes
  2054.  
  2055.     def writeGromacsTop(self, amb2gmx=False):
  2056.         if self.atomTypeSystem == 'amber':
  2057.             d2opls = dictAtomTypeAmb2OplsGmxCode
  2058.         else:
  2059.             d2opls = dictAtomTypeGaff2OplsGmxCode
  2060.  
  2061.         topText = []
  2062.         itpText = []
  2063.         oitpText = []
  2064.         otopText = []
  2065.         top = self.baseName + '_GMX.top'
  2066.         itp = self.baseName + '_GMX.itp'
  2067.         otop = self.baseName + '_GMX_OPLS.top'
  2068.         oitp = self.baseName + '_GMX_OPLS.itp'
  2069.  
  2070.         headDefault = \
  2071.             """
  2072. [ defaults ]
  2073. ; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
  2074. 1               2               yes             0.5     0.8333
  2075. """
  2076.         headItp = \
  2077.             """
  2078. ; Include %s topology
  2079. #include "%s"
  2080. """
  2081.         headOpls = \
  2082.             """
  2083. ; Include forcefield parameters
  2084. #include "ffoplsaa.itp"
  2085. """
  2086.         headSystem = \
  2087.             """
  2088. [ system ]
  2089.  %s
  2090. """
  2091.         headMols = \
  2092.             """
  2093. [ molecules ]
  2094. ; Compound        nmols
  2095. """
  2096.         headAtomtypes = \
  2097.             """
  2098. [ atomtypes ]
  2099. ;name   bond_type     mass     charge   ptype   sigma         epsilon       Amb
  2100. """
  2101.         headAtomtypesOpls = \
  2102.             """
  2103. ; For OPLS atomtypes manual fine tuning
  2104. ; AC_at:OPLS_at:OPLScode: Possible_Aternatives (see ffoplsaa.atp and ffoplsaanb.itp)
  2105. """
  2106.         headMoleculetype = \
  2107.             """
  2108. [ moleculetype ]
  2109. ;name            nrexcl
  2110.  %-16s 3
  2111. """
  2112.         headAtoms = \
  2113.             """
  2114. [ atoms ]
  2115. ;   nr  type  resi  res  atom  cgnr     charge      mass       ; qtot   bond_type
  2116. """
  2117.         headBonds = \
  2118.             """
  2119. [ bonds ]
  2120. ;   ai     aj funct   r             k
  2121. """
  2122.         headPairs = \
  2123.             """
  2124. [ pairs ]
  2125. ;   ai     aj    funct
  2126. """
  2127.         headAngles = \
  2128.             """
  2129. [ angles ]
  2130. ;   ai     aj     ak    funct   theta         cth
  2131. """
  2132.         headProDih = \
  2133.             """
  2134. [ dihedrals ] ; propers
  2135. ; treated as RBs in GROMACS to use combine multiple AMBER torsions per quartet
  2136. ;    i      j      k      l   func    C0         C1         C2         C3         C4         C5
  2137. """
  2138.  
  2139.         headProDihAlphaGamma = """; treated as usual propers in GROMACS since Phase angle diff from 0 or 180 degrees
  2140. ;    i      j      k      l   func   phase     kd      pn
  2141. """
  2142.  
  2143.         headProDihGmx45 = \
  2144.             """
  2145. [ dihedrals ] ; propers
  2146. ; for gromacs 4.5 or higher, using funct 9
  2147. ;    i      j      k      l   func   phase     kd      pn
  2148. """
  2149.  
  2150.         headImpDih = \
  2151.             """
  2152. [ dihedrals ] ; impropers
  2153. ; treated as propers in GROMACS to use correct AMBER analytical function
  2154. ;    i      j      k      l   func   phase     kd      pn
  2155. """
  2156.         _headTopWaterTip3p = \
  2157.             """
  2158. [ bondtypes ]
  2159.   ; i    j      func       b0          kb
  2160.   OW    HW         1    0.09572   462750.4 ; TIP3P water
  2161.   HW    HW         1    0.15139   462750.4 ; TIP3P water
  2163. [ angletypes ]
  2164.   ;  i    j    k  func       th0       cth
  2165.   HW  OW  HW           1   104.520    836.800 ; TIP3P water
  2166.   HW  HW  OW           1   127.740      0.000 ; (found in crystallographic water with 3 bonds)
  2167. """
  2168.  
  2169.         _headTopWaterSpce = \
  2170.             """
  2171. [ bondtypes ]
  2172.   ; i    j      func       b0          kb
  2173.   OW    HW         1    0.1       462750.4 ; SPCE water
  2174.   HW    HW         1    0.1633    462750.4 ; SPCE water
  2176. [ angletypes ]
  2177.   ;  i    j    k  func       th0       cth
  2178.   HW  OW  HW           1   109.47      836.800 ; SPCE water
  2179.   HW  HW  OW           1   125.265     0.000 ; SPCE water
  2180. """
  2181. # NOTE: headTopWaterTip3p and headTopWaterSpce actually do NOTHING
  2182.         headNa = \
  2183.             """
  2184. [ moleculetype ]
  2185.   ; molname       nrexcl
  2186.   NA+             1
  2188. [ atoms ]
  2189.   ; id_    at type res nr  residu name     at name  cg nr  charge   mass
  2190.     1       IP      1          NA+         NA+       1      1     22.9898
  2191. """
  2192.         headCl = \
  2193.             """
  2194. [ moleculetype ]
  2195.   ; molname       nrexcl
  2196.   CL-             1
  2198. [ atoms ]
  2199.   ; id_    at type res nr  residu name     at name  cg nr  charge   mass
  2200.     1       IM      1         CL-           CL-      1     -1     35.45300
  2201. """
  2202.         headK = \
  2203.             """
  2204. [ moleculetype ]
  2205.   ; molname       nrexcl
  2206.   K+             1
  2208. [ atoms ]
  2209.   ; id_    at type res nr  residu name     at name  cg nr  charge   mass
  2210.     1       K       1          K+         K+       1      1     39.100
  2211. """
  2212.         headWaterTip3p = \
  2213.             """
  2214. [ moleculetype ]
  2215. ; molname       nrexcl ; TIP3P model
  2216.   WAT             2
  2218. [ atoms ]
  2219. ;   nr   type  resnr residue  atom   cgnr     charge       mass
  2220.      1     OW      1     WAT     O      1     -0.834   16.00000
  2221.      2     HW      1     WAT    H1      1      0.417    1.00800
  2222.      3     HW      1     WAT    H2      1      0.417    1.00800
  2224. #ifdef FLEXIBLE
  2225. [ bonds ]
  2226. ; i j   funct   length  force.c.
  2227. 1   2   1   0.09572   462750.4 0.09572   462750.4
  2228. 1   3   1   0.09572   462750.4 0.09572   462750.4
  2230. [ angles ]
  2231. ; i j   k   funct   angle   force.c.
  2232. 2   1   3   1   104.520    836.800  104.520    836.800
  2233. #else
  2234. [ settles ]
  2235. ; i j   funct   length
  2236. 1   1   0.09572 0.15139
  2238. [ exclusions ]
  2239. 1   2   3
  2240. 2   1   3
  2241. 3   1   2
  2242. #endif
  2243. """
  2244.  
  2245.         headWaterSpce = \
  2246.             """
  2247. [ moleculetype ]
  2248. ; molname       nrexcl ; SPCE model
  2249.   WAT             2
  2251. [ atoms ]
  2252. ;   nr   type  resnr residue  atom   cgnr     charge       mass
  2253.      1     OW      1     WAT     O      1     -0.8476  15.99940
  2254.      2     HW      1     WAT    H1      1      0.4238   1.00800
  2255.      3     HW      1     WAT    H2      1      0.4238   1.00800
  2257. #ifdef FLEXIBLE
  2258. [ bonds ]
  2259. ; i j   funct   length  force.c.
  2260. 1   2   1   0.1 462750.4  0.1     462750.4
  2261. 1   3   1   0.1 462750.4  0.1     462750.4
  2263. [ angles ]
  2264. ; i j   k   funct   angle   force.c.
  2265. 2   1   3   1   109.47  836.800 109.47  836.800
  2266. #else
  2267. [ settles ]
  2268. ; OW    funct   doh dhh
  2269. 1   1   0.1 0.16330
  2271. [ exclusions ]
  2272. 1   2   3
  2273. 2   1   3
  2274. 3   1   2
  2275. #endif
  2276. """
  2277.         if self.direct and amb2gmx:
  2278.             self.printMess("Converting directly from AMBER to GROMACS.\n")
  2279.  
  2280.         # Dict of ions dealt by acpype emulating amb2gmx
  2281.         ionsDict = {'Na+': headNa, 'Cl-': headCl, 'K+': headK}
  2282.         ionsSorted = []
  2283. # NOTE: headWaterTip3p and headWaterSpce actually do the real thing
  2284. #      so, skipping headTopWaterTip3p and headWaterTip3p
  2285.         # headTopWater = headTopWaterTip3p
  2286.         headWater = headWaterTip3p
  2287.  
  2288.         nWat = 0
  2289.         # topFile.write("; " + head % (top, date))
  2290.         topText.append("; " + head % (top, date))
  2291.         otopText.append("; " + head % (otop, date))
  2292.         # topFile.write(headDefault)
  2293.         topText.append(headDefault)
  2294.  
  2295.         nSolute = 0
  2296.         if not amb2gmx:
  2297.             topText.append(headItp % (itp, itp))
  2298.             otopText.append(headOpls)
  2299.             otopText.append(headItp % (itp, itp))
  2300.             itpText.append("; " + head % (itp, date))
  2301.             oitpText.append("; " + head % (oitp, date))
  2302.  
  2303.         self.printDebug("atomTypes %i" % len(self.atomTypesGromacs))
  2304.         temp = []
  2305.         otemp = []
  2306.         for aType in self.atomTypesGromacs:
  2307.             aTypeName = aType.atomTypeName
  2308.             oaCode = d2opls.get(aTypeName, ['x', '0'])[:-1]
  2309.             aTypeNameOpls = oplsCode2AtomTypeDict.get(oaCode[0], 'x')
  2310.             A = aType.ACOEF
  2311.             B = aType.BCOEF
  2312.             # one cannot infer sigma or epsilon for B = 0, assuming 0 for them
  2313.             if B == 0.0:
  2314.                 sigma, epsilon, r0, epAmber = 0, 0, 0, 0
  2315.             else:
  2316.                 r0 = 0.5 * math.pow((2 * A / B), (1.0 / 6))
  2317.                 epAmber = 0.25 * B * B / A
  2318.                 sigma = 0.1 * math.pow((A / B), (1.0 / 6))
  2319.                 epsilon = cal * epAmber
  2320.             if aTypeName == 'OW':
  2321.                 if A == 629362.166 and B == 625.267765:
  2322.                     # headTopWater = headTopWaterSpce
  2323.                     headWater = headWaterSpce
  2324.             # OW 629362.166 625.267765 spce
  2325.             # OW 581935.564 594.825035 tip3p
  2326.             #       print aTypeName, A, B
  2327.             line = " %-8s %-11s %3.5f  %3.5f   A   %13.5e %13.5e" % \
  2328.                 (aTypeName, aTypeName, 0.0, 0.0, sigma, epsilon) + \
  2329.                 " ; %4.2f  %1.4f\n" % (r0, epAmber)
  2330.             oline = "; %s:%s:opls_%s: %s\n" % (aTypeName, aTypeNameOpls, oaCode[0], repr(oaCode[1:]))
  2331.             # tmpFile.write(line)
  2332.             temp.append(line)
  2333.             otemp.append(oline)
  2334.         if amb2gmx:
  2335.             topText.append(headAtomtypes)
  2336.             topText += temp
  2337.             nWat = self.residueLabel.count('WAT')
  2338.             for ion in ionsDict:
  2339.                 nIon = self.residueLabel.count(ion)
  2340.                 if nIon > 0:
  2341.                     idIon = self.residueLabel.index(ion)
  2342.                     ionsSorted.append((idIon, nIon, ion))
  2343.             ionsSorted.sort()
  2344.         else:
  2345.             itpText.append(headAtomtypes)
  2346.             itpText += temp
  2347.             oitpText.append(headAtomtypesOpls)
  2348.             oitpText += otemp
  2349.         self.printDebug("GMX atomtypes done")
  2350.  
  2351.         if len(self.atoms) > 3 * nWat + sum([x[1] for x in ionsSorted]):
  2352.             nSolute = 1
  2353.  
  2354.         if nWat:
  2355.             # topText.append(headTopWater)
  2356.             self.printDebug("type of water '%s'" % headWater[43:48].strip())
  2357.  
  2358.         if nSolute:
  2359.             if amb2gmx:
  2360.                 topText.append(headMoleculetype % self.baseName)
  2361.             else:
  2362.                 itpText.append(headMoleculetype % self.baseName)
  2363.                 oitpText.append(headMoleculetype % self.baseName)
  2364.  
  2365.         self.printDebug("atoms %i" % len(self.atoms))
  2366.         qtot = 0.0
  2367.         count = 1
  2368.         temp = []
  2369.         otemp = []
  2370.         id2oplsATDict = {}
  2371.         for atom in self.atomsGromacs:
  2372.             resid = atom.resid
  2373.             resname = self.residueLabel[resid]
  2374.             if not self.direct:
  2375.                 if resname in list(ionsDict.keys()) + ['WAT']:
  2376.                     break
  2377.             aName = atom.atomName
  2378.             aType = atom.atomType.atomTypeName
  2379.             oItem = d2opls.get(aType, ['x', 0])
  2380.             oplsAtName = oplsCode2AtomTypeDict.get(oItem[0], 'x')
  2381.             id_ = atom.id
  2382.             id2oplsATDict[id_] = oplsAtName
  2383.             oaCode = 'opls_' + oItem[0]
  2384.             cgnr = id_
  2385.             if self.sorted:
  2386.                 cgnr = atom.cgnr  # JDC
  2387.             charge = atom.charge
  2388.             mass = atom.mass
  2389.             omass = float(oItem[-1])
  2390.             qtot += charge
  2391.             resnr = resid + 1
  2392.             line = "%6d %4s %5d %5s %5s %4d %12.6f %12.5f ; qtot %1.3f\n" % \
  2393.                 (id_, aType, resnr, resname, aName, cgnr, charge, mass, qtot)  # JDC
  2394.             oline = "%6d %4s %5d %5s %5s %4d %12.6f %12.5f ; qtot % 3.3f  %-4s\n" % \
  2395.                 (id_, oaCode, resnr, resname, aName, cgnr, charge, omass, qtot, oplsAtName)  # JDC
  2396.             count += 1
  2397.             temp.append(line)
  2398.             otemp.append(oline)
  2399.         if temp:
  2400.             if amb2gmx:
  2401.                 topText.append(headAtoms)
  2402.                 topText += temp
  2403.             else:
  2404.                 itpText.append(headAtoms)
  2405.                 itpText += temp
  2406.                 oitpText.append(headAtoms)
  2407.                 oitpText += otemp
  2408.         self.printDebug("GMX atoms done")
  2409.  
  2410.         # remove bond of water
  2411.         self.printDebug("bonds %i" % len(self.bonds))
  2412.         temp = []
  2413.         otemp = []
  2414.         for bond in self.bonds:
  2415.             res1 = self.residueLabel[bond.atoms[0].resid]
  2416.             res2 = self.residueLabel[bond.atoms[0].resid]
  2417.             if 'WAT' in [res1, res2]:
  2418.                 continue
  2419.             a1Name = bond.atoms[0].atomName
  2420.             a2Name = bond.atoms[1].atomName
  2421.             id1 = bond.atoms[0].id
  2422.             id2 = bond.atoms[1].id
  2423.             oat1 = id2oplsATDict.get(id1)
  2424.             oat2 = id2oplsATDict.get(id2)
  2425.             line = "%6i %6i %3i %13.4e %13.4e ; %6s - %-6s\n" % (id1, id2, 1,
  2426.                                                                  bond.rEq * 0.1, bond.kBond * 200 * cal, a1Name, a2Name)
  2427.             oline = "%6i %6i %3i ; %13.4e %13.4e ; %6s - %-6s %6s - %-6s\n" % \
  2428.                 (id1, id2, 1, bond.rEq * 0.1, bond.kBond * 200 * cal, a1Name,
  2429.                  a2Name, oat1, oat2)
  2430.             temp.append(line)
  2431.             otemp.append(oline)
  2432.         temp.sort()
  2433.         otemp.sort()
  2434.         if temp:
  2435.             if amb2gmx:
  2436.                 topText.append(headBonds)
  2437.                 topText += temp
  2438.             else:
  2439.                 itpText.append(headBonds)
  2440.                 itpText += temp
  2441.                 oitpText.append(headBonds)
  2442.                 oitpText += otemp
  2443.         self.printDebug("GMX bonds done")
  2444.  
  2445.         self.printDebug("atomPairs %i" % len(self.atomPairs))
  2446.         temp = []
  2447.         for pair in self.atomPairs:
  2448.             # if not printed:
  2449.             #    tmpFile.write(headPairs)
  2450.             #    printed = True
  2451.             a1Name = pair[0].atomName
  2452.             a2Name = pair[1].atomName
  2453.             id1 = pair[0].id
  2454.             id2 = pair[1].id
  2455.             # id1 = self.atoms.index(pair[0]) + 1
  2456.             # id2 = self.atoms.index(pair[1]) + 1
  2457.             line = "%6i %6i %6i ; %6s - %-6s\n" % (id1, id2, 1, a1Name,
  2458.                                                    a2Name)
  2459.             temp.append(line)
  2460.         temp.sort()
  2461.         if temp:
  2462.             if amb2gmx:
  2463.                 topText.append(headPairs)
  2464.                 topText += temp
  2465.             else:
  2466.                 itpText.append(headPairs)
  2467.                 itpText += temp
  2468.                 oitpText.append(headPairs)
  2469.                 oitpText += temp
  2470.         self.printDebug("GMX pairs done")
  2471.  
  2472.         self.printDebug("angles %i" % len(self.angles))
  2473.         temp = []
  2474.         otemp = []
  2475.         for angle in self.angles:
  2476.             a1 = angle.atoms[0].atomName
  2477.             a2 = angle.atoms[1].atomName
  2478.             a3 = angle.atoms[2].atomName
  2479.             id1 = angle.atoms[0].id
  2480.             id2 = angle.atoms[1].id
  2481.             id3 = angle.atoms[2].id
  2482.             oat1 = id2oplsATDict.get(id1)
  2483.             oat2 = id2oplsATDict.get(id2)
  2484.             oat3 = id2oplsATDict.get(id3)
  2485.             line = "%6i %6i %6i %6i %13.4e %13.4e ; %6s - %-6s - %-6s\n" % (id1, id2,
  2486.                                                                             id3, 1, angle.thetaEq * radPi, 2 * cal * angle.kTheta, a1, a2, a3)
  2487.             oline = "%6i %6i %6i %6i ; %13.4e %13.4e ; %6s - %-4s - %-6s %4s - %+4s - %-4s\n" % \
  2488.                 (id1, id2, id3, 1, angle.thetaEq * radPi, 2 * cal * angle.kTheta,
  2489.                  a1, a2, a3, oat1, oat2, oat3)
  2490.             temp.append(line)
  2491.             otemp.append(oline)
  2492.         temp.sort()
  2493.         otemp.sort()
  2494.         if temp:
  2495.             if amb2gmx:
  2496.                 topText.append(headAngles)
  2497.                 topText += temp
  2498.             else:
  2499.                 itpText.append(headAngles)
  2500.                 itpText += temp
  2501.                 oitpText.append(headAngles)
  2502.                 oitpText += otemp
  2503.         self.printDebug("GMX angles done")
  2504.  
  2505.         self.setProperDihedralsCoef()
  2506.         self.printDebug("properDihedralsCoefRB %i" % len(self.properDihedralsCoefRB))
  2507.         self.printDebug("properDihedralsAlphaGamma %i" % len(self.properDihedralsAlphaGamma))
  2508.         self.printDebug("properDihedralsGmx45 %i" % len(self.properDihedralsGmx45))
  2509.         temp = []
  2510.         otemp = []
  2511.         if self.gmx4:
  2512.             self.printMess("Writing RB dihedrals for old GMX 4.\n")
  2513.             for dih in self.properDihedralsCoefRB:
  2514.                 a1 = dih[0][0].atomName
  2515.                 a2 = dih[0][1].atomName
  2516.                 a3 = dih[0][2].atomName
  2517.                 a4 = dih[0][3].atomName
  2518.                 id1 = dih[0][0].id
  2519.                 id2 = dih[0][1].id
  2520.                 id3 = dih[0][2].id
  2521.                 id4 = dih[0][3].id
  2522.                 oat1 = id2oplsATDict.get(id1)
  2523.                 oat2 = id2oplsATDict.get(id2)
  2524.                 oat3 = id2oplsATDict.get(id3)
  2525.                 oat4 = id2oplsATDict.get(id4)
  2526.                 c0, c1, c2, c3, c4, c5 = dih[1]
  2527.                 line = \
  2528.                     "%6i %6i %6i %6i %6i %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f" % \
  2529.                     (id1, id2, id3, id4, 3, c0, c1, c2, c3, c4, c5) \
  2530.                     + " ; %6s-%6s-%6s-%6s\n" % (a1, a2, a3, a4)
  2531.                 oline = \
  2532.                     "%6i %6i %6i %6i %6i ; %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f" % \
  2533.                     (id1, id2, id3, id4, 3, c0, c1, c2, c3, c4, c5) \
  2534.                     + " ; %6s-%6s-%6s-%6s    %4s-%4s-%4s-%4s\n" % (a1, a2, a3, a4, oat1, oat2, oat3, oat4)
  2535.                 temp.append(line)
  2536.                 otemp.append(oline)
  2537.             temp.sort()
  2538.             otemp.sort()
  2539.             if temp:
  2540.                 if amb2gmx:
  2541.                     topText.append(headProDih)
  2542.                     topText += temp
  2543.                 else:
  2544.                     itpText.append(headProDih)
  2545.                     itpText += temp
  2546.                     oitpText.append(headProDih)
  2547.                     oitpText += otemp
  2548.             self.printDebug("GMX proper dihedrals done")
  2549.         else:
  2550.             self.printMess("Writing GMX dihedrals for GMX 4.5 and higher.\n")
  2551.             funct = 9  # 9
  2552.             for dih in self.properDihedralsGmx45:
  2553.                 a1 = dih[0][0].atomName
  2554.                 a2 = dih[0][1].atomName
  2555.                 a3 = dih[0][2].atomName
  2556.                 a4 = dih[0][3].atomName
  2557.                 id1 = dih[0][0].id
  2558.                 id2 = dih[0][1].id
  2559.                 id3 = dih[0][2].id
  2560.                 id4 = dih[0][3].id
  2561.                 ph = dih[1]  # phase already in degree
  2562.                 kd = dih[2] * cal  # kPhi PK
  2563.                 pn = dih[3]  # .period
  2564.                 line = "%6i %6i %6i %6i %6i %8.2f %9.5f %3i ; %6s-%6s-%6s-%6s\n" % \
  2565.                     (id1, id2, id3, id4, funct, ph, kd, pn, a1, a2, a3, a4)
  2566.                 oline = "%6i %6i %6i %6i %6i ; %8.2f %9.5f %3i ; %6s-%6s-%6s-%6s\n" % \
  2567.                     (id1, id2, id3, id4, funct, ph, kd, pn, a1, a2, a3, a4)
  2568.                 temp.append(line)
  2569.                 otemp.append(oline)
  2570.             temp.sort()
  2571.             otemp.sort()
  2572.             if temp:
  2573.                 if amb2gmx:
  2574.                     topText.append(headProDihGmx45)
  2575.                     topText += temp
  2576.                 else:
  2577.                     itpText.append(headProDihGmx45)
  2578.                     itpText += temp
  2579.                     oitpText.append(headProDihGmx45)
  2580.                     oitpText += otemp
  2581.  
  2582.         # for properDihedralsAlphaGamma
  2583.         if not self.gmx4:
  2584.             funct = 4  # 4
  2585.         else:
  2586.             funct = 1
  2587.         temp = []
  2588.         otemp = []
  2589.         for dih in self.properDihedralsAlphaGamma:
  2590.             a1 = dih[0][0].atomName
  2591.             a2 = dih[0][1].atomName
  2592.             a3 = dih[0][2].atomName
  2593.             a4 = dih[0][3].atomName
  2594.             id1 = dih[0][0].id
  2595.             id2 = dih[0][1].id
  2596.             id3 = dih[0][2].id
  2597.             id4 = dih[0][3].id
  2598.             ph = dih[1]  # phase already in degree
  2599.             kd = dih[2] * cal  # kPhi PK
  2600.             pn = dih[3]  # .period
  2601.             line = "%6i %6i %6i %6i %6i %8.2f %9.5f %3i ; %6s-%6s-%6s-%6s\n" % \
  2602.                 (id1, id2, id3, id4, funct, ph, kd, pn, a1, a2, a3, a4)
  2603.             oline = "%6i %6i %6i %6i %6i ; %8.2f %9.5f %3i ; %6s-%6s-%6s-%6s\n" % \
  2604.                 (id1, id2, id3, id4, funct, ph, kd, pn, a1, a2, a3, a4)
  2605.             temp.append(line)
  2606.             otemp.append(oline)
  2607.         temp.sort()
  2608.         otemp.sort()
  2609.         if temp:
  2610.             if amb2gmx:
  2611.                 topText.append(headProDihAlphaGamma)
  2612.                 topText += temp
  2613.             else:
  2614.                 itpText.append(headProDihAlphaGamma)
  2615.                 itpText += temp
  2616.                 oitpText.append(headProDihAlphaGamma)
  2617.                 oitpText += otemp
  2618.         self.printDebug("GMX special proper dihedrals done")
  2619.  
  2620.         self.printDebug("improperDihedrals %i" % len(self.improperDihedrals))
  2621.         temp = []
  2622.         otemp = []
  2623.         for dih in self.improperDihedrals:
  2624.             a1 = dih.atoms[0].atomName
  2625.             a2 = dih.atoms[1].atomName
  2626.             a3 = dih.atoms[2].atomName
  2627.             a4 = dih.atoms[3].atomName
  2628.             id1 = dih.atoms[0].id
  2629.             id2 = dih.atoms[1].id
  2630.             id3 = dih.atoms[2].id
  2631.             id4 = dih.atoms[3].id
  2632.             kd = dih.kPhi * cal
  2633.             pn = dih.period
  2634.             ph = dih.phase * radPi
  2635.             line = "%6i %6i %6i %6i %6i %8.2f %9.5f %3i ; %6s-%6s-%6s-%6s\n" % \
  2636.                 (id1, id2, id3, id4, funct, ph, kd, pn, a1, a2, a3, a4)
  2637.             oline = "%6i %6i %6i %6i %6i ; %8.2f %9.5f %3i ; %6s-%6s-%6s-%6s\n" % \
  2638.                 (id1, id2, id3, id4, funct, ph, kd, pn, a1, a2, a3, a4)
  2639.             temp.append(line)
  2640.             otemp.append(oline)
  2641.         temp.sort()
  2642.         otemp.sort()
  2643.         if temp:
  2644.             if amb2gmx:
  2645.                 topText.append(headImpDih)
  2646.                 topText += temp
  2647.             else:
  2648.                 itpText.append(headImpDih)
  2649.                 itpText += temp
  2650.                 oitpText.append(headImpDih)
  2651.                 oitpText += otemp
  2652.         self.printDebug("GMX improper dihedrals done")
  2653.  
  2654.         if not self.direct:
  2655.             for ion in ionsSorted:
  2656.                 topText.append(ionsDict[ion[2]])
  2657.  
  2658.             if nWat:
  2659.                 topText.append(headWater)
  2660.  
  2661.         topText.append(headSystem % (self.baseName))
  2662.         topText.append(headMols)
  2663.         otopText.append(headSystem % (self.baseName))
  2664.         otopText.append(headMols)
  2665.  
  2666.         if nSolute > 0:
  2667.             topText.append(" %-16s %-6i\n" % (self.baseName, nSolute))
  2668.             otopText.append(" %-16s %-6i\n" % (self.baseName, nSolute))
  2669.  
  2670.         if not self.direct:
  2671.             for ion in ionsSorted:
  2672.                 topText.append(" %-16s %-6i\n" % (ion[2].upper(), ion[1]))
  2673.  
  2674.             if nWat:
  2675.                 topText.append(" %-16s %-6i\n" % ('WAT', nWat))
  2676.  
  2677.         gmxDir = os.path.abspath('.')
  2678.         topFileName = os.path.join(gmxDir, top)
  2679.         topFile = open(topFileName, 'w')
  2680.         topFile.writelines(topText)
  2681.  
  2682.         if not amb2gmx:
  2683.             itpFileName = os.path.join(gmxDir, itp)
  2684.             itpFile = open(itpFileName, 'w')
  2685.             itpFile.writelines(itpText)
  2686.             oitpFileName = os.path.join(gmxDir, oitp)
  2687.             oitpFile = open(oitpFileName, 'w')
  2688.             oitpFile.writelines(oitpText)
  2689.             otopFileName = os.path.join(gmxDir, otop)
  2690.             otopFile = open(otopFileName, 'w')
  2691.             otopFile.writelines(otopText)
  2692.  
  2693.     def writeGroFile(self):
  2694.         # print "Writing GROMACS GRO file\n"
  2695.         self.printDebug("writing GRO file")
  2696.         gro = self.baseName + '_GMX.gro'
  2697.         gmxDir = os.path.abspath('.')
  2698.         groFileName = os.path.join(gmxDir, gro)
  2699.         groFile = open(groFileName, 'w')
  2700.         groFile.write(head % (gro, date))
  2701.         groFile.write(" %i\n" % len(self.atoms))
  2702.         count = 1
  2703.         for atom in self.atoms:
  2704.             coords = [c * 0.1 for c in atom.coords]
  2705.             resid = atom.resid
  2706.             line = "%5d%5s%5s%5d%8.3f%8.3f%8.3f\n" % \
  2707.                    (resid + 1, self.residueLabel[resid], atom.atomName,
  2708.                     count, coords[0], coords[1], coords[2])
  2709.             count += 1
  2710.             if count == 100000:
  2711.                 count = 0
  2712.             groFile.write(line)
  2713.         if self.pbc:
  2714.             boxX = self.pbc[0][0] * 0.1
  2715.             boxY = self.pbc[0][1] * 0.1
  2716.             boxZ = self.pbc[0][2] * 0.1
  2717.             vX = self.pbc[1][0]
  2718.             # vY = self.pbc[1][1]
  2719.             # vZ = self.pbc[1][2]
  2720.             if vX == 90.0:
  2721.                 self.printDebug("PBC triclinic")
  2722.                 text = "%11.5f %11.5f %11.5f\n" % (boxX, boxY, boxZ)
  2723.             elif round(vX, 2) == 109.47:
  2724.                 self.printDebug("PBC octahedron")
  2725.                 f1 = 0.471405  # 1/3 * sqrt(2)
  2726.                 f2 = 0.333333 * boxX
  2727.                 v22 = boxY * 2 * f1
  2728.                 v33 = boxZ * f1 * 1.73205  # f1 * sqrt(3)
  2729.                 v21 = v31 = v32 = 0.0
  2730.                 v12 = f2
  2731.                 v13 = -f2
  2732.                 v23 = f1 * boxX
  2733.                 text = "%11.5f %11.5f %11.5f %11.5f %11.5f %11.5f %11.5f %11.5f %11.5f\n" % \
  2734.                     (boxX, v22, v33, v21, v31, v12, v32, v13, v23)
  2735.         else:
  2736.             self.printDebug("Box size estimated")
  2737.             X = [a.coords[0] * 0.1 for a in self.atoms]
  2738.             Y = [a.coords[1] * 0.1 for a in self.atoms]
  2739.             Z = [a.coords[2] * 0.1 for a in self.atoms]
  2740.             boxX = max(X) - min(X)  # + 2.0 # 2.0 is double of rlist
  2741.             boxY = max(Y) - min(Y)  # + 2.0
  2742.             boxZ = max(Z) - min(Z)  # + 2.0
  2743.             text = "%11.5f %11.5f %11.5f\n" % (boxX * 20.0, boxY * 20.0, boxZ * 20.0)
  2744.         groFile.write(text)
  2745.  
  2746.     def writeMdpFiles(self):
  2747.         emMdp = """; to test
  2748. ; gmx grompp -f em.mdp -c {base}_GMX.gro -p {base}_GMX.top -o em.tpr -v
  2749. ; gmx mdrun -ntmpi 1 -v -deffnm em
  2750. integrator               = steep
  2751. nsteps                   = 500
  2752. """.format(base=self.baseName)
  2753.         mdMdp = """; to test
  2754. ; gmx grompp -f md.mdp -c em.gro -p {base}_GMX.top -o md.tpr
  2755. ; gmx mdrun -ntmpi 1 -v -deffnm md
  2756. integrator               = md
  2757. nsteps                   = 10000
  2758. """.format(base=self.baseName)
  2759.         emMdpFile = open('em.mdp', 'w')
  2760.         mdMdpFile = open('md.mdp', 'w')
  2761.         emMdpFile.write(emMdp)
  2762.         mdMdpFile.write(mdMdp)
  2763.  
  2764.     def writeCnsTopolFiles(self):
  2765.         autoAngleFlag = True
  2766.         autoDihFlag = True
  2767.         cnsDir = os.path.abspath('.')
  2768.  
  2769.         pdb = self.baseName + '_NEW.pdb'
  2770.         par = self.baseName + '_CNS.par'
  2771.         top = self.baseName + '_CNS.top'
  2772.         inp = self.baseName + '_CNS.inp'
  2773.  
  2774.         pdbFileName = os.path.join(cnsDir, pdb)
  2775.         parFileName = os.path.join(cnsDir, par)
  2776.         topFileName = os.path.join(cnsDir, top)
  2777.         inpFileName = os.path.join(cnsDir, inp)
  2778.  
  2779.         self.CnsTopFileName = topFileName
  2780.         self.CnsInpFileName = inpFileName
  2781.         self.CnsParFileName = parFileName
  2782.         self.CnsPdbFileName = pdbFileName
  2783.  
  2784.         parFile = open(parFileName, 'w')
  2785.         topFile = open(topFileName, 'w')
  2786.         inpFile = open(inpFileName, 'w')
  2787.  
  2788.         self.printMess("Writing NEW PDB file\n")
  2789.         self.writePdb(pdbFileName)
  2790.  
  2791.         self.printMess("Writing CNS/XPLOR files\n")
  2792.  
  2793.         # print "Writing CNS PAR file\n"
  2794.         parFile.write("Remarks " + head % (par, date))
  2795.         parFile.write("\nset echo=false end\n")
  2796.  
  2797.         parFile.write("\n{ Bonds: atomType1 atomType2 kb r0 }\n")
  2798.         lineSet = []
  2799.         for bond in self.bonds:
  2800.             a1Type = bond.atoms[0].atomType.atomTypeName + '_'
  2801.             a2Type = bond.atoms[1].atomType.atomTypeName + '_'
  2802.             kb = 1000.0
  2803.             if not self.allhdg:
  2804.                 kb = bond.kBond
  2805.             r0 = bond.rEq
  2806.             line = "BOND %5s %5s %8.1f %8.4f\n" % (a1Type, a2Type, kb, r0)
  2807.             lineRev = "BOND %5s %5s %8.1f %8.4f\n" % (a2Type, a1Type, kb, r0)
  2808.             if line not in lineSet:
  2809.                 if lineRev not in lineSet:
  2810.                     lineSet.append(line)
  2811.         for item in lineSet:
  2812.             parFile.write(item)
  2813.  
  2814.         parFile.write("\n{ Angles: aType1 aType2 aType3 kt t0 }\n")
  2815.         lineSet = []
  2816.         for angle in self.angles:
  2817.             a1 = angle.atoms[0].atomType.atomTypeName + '_'
  2818.             a2 = angle.atoms[1].atomType.atomTypeName + '_'
  2819.             a3 = angle.atoms[2].atomType.atomTypeName + '_'
  2820.             kt = 500.0
  2821.             if not self.allhdg:
  2822.                 kt = angle.kTheta
  2823.             t0 = angle.thetaEq * radPi
  2824.             line = "ANGLe %5s %5s %5s %8.1f %8.2f\n" % (a1, a2, a3, kt, t0)
  2825.             lineRev = "ANGLe %5s %5s %5s %8.1f %8.2f\n" % (a3, a2, a1, kt, t0)
  2826.             if line not in lineSet:
  2827.                 if lineRev not in lineSet:
  2828.                     lineSet.append(line)
  2829.         for item in lineSet:
  2830.             parFile.write(item)
  2831.  
  2832.         parFile.write("\n{ Proper Dihedrals: aType1 aType2 aType3 aType4 kt per\
  2833. iod phase }\n")
  2834.         lineSet = set()
  2835.         for item in self.condensedProperDihedrals:
  2836.             seq = ''
  2837.             id_ = 0
  2838.             for dih in item:
  2839.                 # id_ = item.index(dih)
  2840.                 l = len(item)
  2841.                 a1 = dih.atoms[0].atomType.atomTypeName + '_'
  2842.                 a2 = dih.atoms[1].atomType.atomTypeName + '_'
  2843.                 a3 = dih.atoms[2].atomType.atomTypeName + '_'
  2844.                 a4 = dih.atoms[3].atomType.atomTypeName + '_'
  2845.                 kp = 750.0
  2846.                 if not self.allhdg:
  2847.                     kp = dih.kPhi
  2848.                 p = dih.period
  2849.                 ph = dih.phase * radPi
  2850.                 if l > 1:
  2851.                     if id_ == 0:
  2852.                         line = "DIHEdral %5s %5s %5s %5s  MULT %1i %7.3f %4i %8\
  2853. .2f\n" % (a1, a2, a3, a4, l, kp, p, ph)
  2854.                     else:
  2855.                         line = "%s %7.3f %4i %8.2f\n" % (40 * " ", kp, p, ph)
  2856.                 else:
  2857.                     line = "DIHEdral %5s %5s %5s %5s %15.3f %4i %8.2f\n" % (a1,
  2858.                                                                             a2, a3, a4, kp, p, ph)
  2859.                 seq += line
  2860.                 id_ += 1
  2861.             lineSet.add(seq)
  2862.         for item in lineSet:
  2863.             parFile.write(item)
  2864.  
  2865.         parFile.write("\n{ Improper Dihedrals: aType1 aType2 aType3 aType4 kt p\
  2866. eriod phase }\n")
  2867.         lineSet = set()
  2868.         for idh in self.improperDihedrals:
  2869.             a1 = idh.atoms[0].atomType.atomTypeName + '_'
  2870.             a2 = idh.atoms[1].atomType.atomTypeName + '_'
  2871.             a3 = idh.atoms[2].atomType.atomTypeName + '_'
  2872.             a4 = idh.atoms[3].atomType.atomTypeName + '_'
  2873.             kp = 750.0
  2874.             if not self.allhdg:
  2875.                 kp = idh.kPhi
  2876.             p = idh.period
  2877.             ph = idh.phase * radPi
  2878.             line = "IMPRoper %5s %5s %5s %5s %13.1f %4i %8.2f\n" % (a1, a2, a3,
  2879.                                                                     a4, kp, p, ph)
  2880.             lineSet.add(line)
  2881.  
  2882.         if self.chiral:
  2883.             for idhc in self.chiralGroups:
  2884.                 _atc, neig, angle = idhc
  2885.                 a1 = neig[0].atomType.atomTypeName + '_'
  2886.                 a2 = neig[1].atomType.atomTypeName + '_'
  2887.                 a3 = neig[2].atomType.atomTypeName + '_'
  2888.                 a4 = neig[3].atomType.atomTypeName + '_'
  2889.                 kp = 11000.0
  2890.                 p = 0
  2891.                 ph = angle
  2892.                 line = "IMPRoper %5s %5s %5s %5s %13.1f %4i %8.2f\n" % (a1, a2, a3,
  2893.                                                                         a4, kp, p, ph)
  2894.                 lineSet.add(line)
  2895.  
  2896.         for item in lineSet:
  2897.             parFile.write(item)
  2898.  
  2899.         parFile.write("\n{ Nonbonded: Type Emin sigma; (1-4): Emin/2 sigma }\n")
  2900.         for at in self.atomTypes:
  2901.             A = at.ACOEF
  2902.             B = at.BCOEF
  2903.             atName = at.atomTypeName + '_'
  2904.             if B == 0.0:
  2905.                 sigma = epAmber = ep2 = sig2 = 0.0
  2906.             else:
  2907.                 epAmber = 0.25 * B * B / A
  2908.                 ep2 = epAmber / 2.0
  2909.                 sigma = math.pow((A / B), (1.0 / 6))
  2910.                 sig2 = sigma
  2911.             line = "NONBonded %5s %11.6f %11.6f %11.6f %11.6f\n" % (atName,
  2912.                                                                     epAmber, sigma, ep2, sig2)
  2913.             parFile.write(line)
  2914.         parFile.write("\nset echo=true end\n")
  2915.  
  2916.         # print "Writing CNS TOP file\n"
  2917.         topFile.write("Remarks " + head % (top, date))
  2918.         topFile.write("\nset echo=false end\n")
  2919.         topFile.write("\nautogenerate angles=%s dihedrals=%s end\n" %
  2920.                       (autoAngleFlag, autoDihFlag))
  2921.  
  2922.         topFile.write("\n{ atomType  mass }\n")
  2923.         for at in self.atomTypes:
  2924.             atType = at.atomTypeName + '_'
  2925.             mass = at.mass
  2926.             line = "MASS %-5s %8.3f\n" % (atType, mass)
  2927.             topFile.write(line)
  2928.  
  2929.         topFile.write("\nRESIdue %s\n" % self.residueLabel[0])
  2930.         topFile.write("\nGROUP\n")
  2931.  
  2932.         topFile.write("\n{ atomName  atomType  Charge }\n")
  2933.         for at in self.atoms:
  2934.             atName = at.atomName
  2935.             atType = at.atomType.atomTypeName + '_'
  2936.             charge = at.charge
  2937.             line = "ATOM %-5s TYPE= %-5s CHARGE= %8.4f END\n" % (atName, atType,
  2938.                                                                  charge)
  2939.             topFile.write(line)
  2940.  
  2941.         topFile.write("\n{ Bonds: atomName1  atomName2 }\n")
  2942.         for bond in self.bonds:
  2943.             a1Name = bond.atoms[0].atomName
  2944.             a2Name = bond.atoms[1].atomName
  2945.             line = "BOND %-5s %-5s\n" % (a1Name, a2Name)
  2946.             topFile.write(line)
  2947.  
  2948.         if not autoAngleFlag or 1:  # generating angles anyway
  2949.             topFile.write("\n{ Angles: atomName1 atomName2 atomName3}\n")
  2950.             for angle in self.angles:
  2951.                 a1Name = angle.atoms[0].atomName
  2952.                 a2Name = angle.atoms[1].atomName
  2953.                 a3Name = angle.atoms[2].atomName
  2954.                 line = "ANGLe %-5s %-5s %-5s\n" % (a1Name, a2Name, a3Name)
  2955.                 topFile.write(line)
  2956.  
  2957.         if not autoDihFlag or 1:  # generating angles anyway
  2958.             topFile.write("\n{ Proper Dihedrals: name1 name2 name3 name4 }\n")
  2959.             for item in self.condensedProperDihedrals:
  2960.                 for dih in item:
  2961.                     l = len(item)
  2962.                     a1Name = dih.atoms[0].atomName
  2963.                     a2Name = dih.atoms[1].atomName
  2964.                     a3Name = dih.atoms[2].atomName
  2965.                     a4Name = dih.atoms[3].atomName
  2966.                     line = "DIHEdral %-5s %-5s %-5s %-5s\n" % (a1Name, a2Name,
  2967.                                                                a3Name, a4Name)
  2968.                     break
  2969.                 topFile.write(line)
  2970.  
  2971.         topFile.write("\n{ Improper Dihedrals: aName1 aName2 aName3 aName4 }\n")
  2972.         for dih in self.improperDihedrals:
  2973.             a1Name = dih.atoms[0].atomName
  2974.             a2Name = dih.atoms[1].atomName
  2975.             a3Name = dih.atoms[2].atomName
  2976.             a4Name = dih.atoms[3].atomName
  2977.             line = "IMPRoper %-5s %-5s %-5s %-5s\n" % (a1Name, a2Name, a3Name,
  2978.                                                        a4Name)
  2979.             topFile.write(line)
  2980.  
  2981.         if self.chiral:
  2982.             for idhc in self.chiralGroups:
  2983.                 _atc, neig, angle = idhc
  2984.                 a1Name = neig[0].atomName
  2985.                 a2Name = neig[1].atomName
  2986.                 a3Name = neig[2].atomName
  2987.                 a4Name = neig[3].atomName
  2988.                 line = "IMPRoper %-5s %-5s %-5s %-5s\n" % (a1Name, a2Name, a3Name,
  2989.                                                            a4Name)
  2990.                 topFile.write(line)
  2991.  
  2992.         topFile.write("\nEND {RESIdue %s}\n" % self.residueLabel[0])
  2993.  
  2994.         topFile.write("\nset echo=true end\n")
  2995.  
  2996.         # print "Writing CNS INP file\n"
  2997.         inpFile.write("Remarks " + head % (inp, date))
  2998.         inpData = \
  2999.             """
  3000. topology
  3001.   @%(CNS_top)s
  3002. end
  3004. parameters
  3005.   @%(CNS_par)s
  3006.   nbonds
  3007.       atom cdie shift eps=1.0  e14fac=0.4   tolerance=0.5
  3008.       cutnb=9.0 ctonnb=7.5 ctofnb=8.0
  3009.       nbxmod=5 vswitch wmin 1.0
  3010.   end
  3011.   remark dielectric constant eps set to 1.0
  3012. end
  3014. flags exclude elec ? end
  3016. segment name="    "
  3017.   chain
  3018.    coordinates @%(NEW_pdb)s
  3019.   end
  3020. end
  3021. coordinates @%(NEW_pdb)s
  3022. coord copy end
  3024. ! Remarks If you want to shake up the coordinates a bit ...
  3025.  vector do (x=x+6*(rand()-0.5)) (all)
  3026.  vector do (y=y+6*(rand()-0.5)) (all)
  3027.  vector do (z=z+6*(rand()-0.5)) (all)
  3028.  write coordinates output=%(CNS_ran)s end
  3030. ! Remarks RMS diff after randomisation and before minimisation
  3031. coord rms sele=(known and not hydrogen) end
  3033. print threshold=0.02 bonds
  3034. print threshold=3.0 angles
  3035. print threshold=3.0 dihedrals
  3036. print threshold=3.0 impropers
  3038. ! Remarks Do Powell energy minimisation
  3039. minimise powell
  3040.   nstep=250 drop=40.0
  3041. end
  3043. write coordinates output=%(CNS_min)s end
  3044. write structure   output=%(CNS_psf)s end
  3046. ! constraints interaction (not hydro) (not hydro) end
  3048. print threshold=0.02 bonds
  3049. print threshold=3.0 angles
  3050. print threshold=3.0 dihedrals
  3051. print threshold=3.0 impropers
  3053. flags exclude * include vdw end energy end
  3054. distance from=(not hydro) to=(not hydro) cutoff=2.6 end
  3056. ! Remarks RMS fit after minimisation
  3057. coord fit sele=(known and not hydrogen) end
  3059. stop
  3060. """
  3061.         dictInp = {}
  3062.         dictInp['CNS_top'] = top
  3063.         dictInp['CNS_par'] = par
  3064.         dictInp['NEW_pdb'] = pdb
  3065.         dictInp['CNS_min'] = self.baseName + '_NEW_min.pdb'
  3066.         dictInp['CNS_psf'] = self.baseName + '_CNS.psf'
  3067.         dictInp['CNS_ran'] = self.baseName + '_rand.pdb'
  3068.         line = inpData % dictInp
  3069.         inpFile.write(line)
  3070.         if os.path.exists(self.obchiralExe):
  3071.             self.printDebug("chiralGroups %i" % len(self.chiralGroups))
  3072.         else:
  3073.             self.printDebug("No 'obchiral' to process chiral atoms. Consider installing http://openbabel.org")
  3074.  
  3075.  
  3076. class ACTopol(AbstractTopol):
  3077.  
  3078.     """
  3079.         Class to build the AC topologies (Antechamber AmberTools)
  3080.     """
  3081.  
  3082.     def __init__(self, inputFile, chargeType='bcc', chargeVal=None,
  3083.                  multiplicity='1', atomType='gaff', force=False, basename=None,
  3084.                  debug=False, outTopol='all', engine='tleap', allhdg=False,
  3085.                  timeTol=36000, qprog='sqm', ekFlag=None, verbose=True,
  3086.                  gmx4=False, disam=False, direct=False, is_sorted=False, chiral=False):
  3087.  
  3088.         self.debug = debug
  3089.         self.verbose = verbose
  3090.         self.gmx4 = gmx4
  3091.         self.disam = disam
  3092.         self.direct = direct
  3093.         self.sorted = is_sorted
  3094.         self.chiral = chiral
  3095.         self.inputFile = os.path.basename(inputFile)
  3096.         self.rootDir = os.path.abspath('.')
  3097.         self.absInputFile = os.path.abspath(inputFile)
  3098.         if not os.path.exists(self.absInputFile):
  3099.             self.printWarn("input file doesn't exist")
  3100.         baseOriginal, ext = os.path.splitext(self.inputFile)
  3101.         base = basename or baseOriginal
  3102.         self.baseOriginal = baseOriginal
  3103.         self.baseName = base  # name of the input file without ext.
  3104.         self.timeTol = timeTol
  3105.         self.printDebug("Max execution time tolerance is %s" % elapsedTime(self.timeTol))
  3106.         self.ext = ext
  3107.         if ekFlag == '"None"' or ekFlag is None:
  3108.             self.ekFlag = ''
  3109.         else:
  3110.             self.ekFlag = '-ek %s' % ekFlag
  3111.         self.extOld = ext
  3112.         self.homeDir = self.baseName + '.acpype'
  3113.         self.chargeType = chargeType
  3114.         self.chargeVal = chargeVal
  3115.         self.multiplicity = multiplicity
  3116.         self.atomType = atomType
  3117.         self.gaffDatfile = 'gaff.dat'
  3118.         leapGaffFile = 'leaprc.gaff'
  3119.         if '2' in self.atomType:
  3120.             leapGaffFile = 'leaprc.gaff2'
  3121.             self.gaffDatfile = 'gaff2.dat'
  3122.         self.force = force
  3123.         self.engine = engine
  3124.         self.allhdg = allhdg
  3125.         self.acExe = ''
  3126.         dirAmber = os.getenv('AMBERHOME', os.getenv('ACHOME'))
  3127.         if dirAmber:
  3128.             for ac_bin in ['bin', 'exe']:
  3129.                 ac_path = os.path.join(dirAmber, ac_bin, 'antechamber')
  3130.                 if os.path.exists(ac_path):
  3131.                     self.acExe = ac_path
  3132.                     break
  3133.         if not self.acExe:
  3134.             self.acExe = _getoutput('which antechamber') or ''  # '/Users/alan/Programmes/antechamber-1.27/exe/antechamber'
  3135.         if not os.path.exists(self.acExe):
  3136.             self.printError("no 'antechamber' executable!")
  3137.             return None
  3138.         self.tleapExe = _getoutput('which tleap') or ''
  3139.         self.sleapExe = _getoutput('which sleap') or ''
  3140.         self.parmchkExe = _getoutput('which parmchk2') or ''
  3141.         self.babelExe = _getoutput('which babel') or ''
  3142.         if not os.path.exists(self.babelExe):
  3143.             if self.ext != '.mol2' and self.ext != '.mdl':  # and self.ext != '.mol':
  3144.                 self.printError("no 'babel' executable; you need it if input is PDB")
  3145.                 self.printError("otherwise use only MOL2 or MDL file as input ... aborting!")
  3146.                 sys.exit(1)
  3147.             else:
  3148.                 self.printWarn("no 'babel' executable, no PDB file as input can be used!")
  3149.         acBase = base + '_AC'
  3150.         self.acBaseName = acBase
  3151.         self.acXyzFileName = acBase + '.inpcrd'
  3152.         self.acTopFileName = acBase + '.prmtop'
  3153.         self.acFrcmodFileName = acBase + '.frcmod'
  3154.         self.tmpDir = os.path.join(self.rootDir, '.acpype_tmp_%s' % os.path.basename(base))
  3155.         self.setResNameCheckCoords()
  3156.         self.guessCharge()
  3157.         acMol2FileName = '%s_%s_%s.mol2' % (base, chargeType, atomType)
  3158.         self.acMol2FileName = acMol2FileName
  3159.         self.charmmBase = '%s_CHARMM' % base
  3160.         # check for which version of antechamber
  3161.         if 'amber10' in self.acExe:
  3162.             if qprog == 'sqm':
  3163.                 self.printWarn("SQM is not implemented in AmberTools 1.2")
  3164.                 self.printWarn("Setting mopac for antechamber")
  3165.                 qprog = 'mopac'
  3166.             elif qprog == 'divcon':
  3167.                 if not os.path.exists(os.path.join(os.path.dirname(self.acExe), qprog)):
  3168.                     self.printWarn("DIVCON is not installed")
  3169.                     self.printWarn("Setting mopac for antechamber")
  3170.                     qprog = 'mopac'
  3171.         elif 'amber1' in self.acExe:
  3172.             if qprog == 'divcon':
  3173.                 self.printWarn("DIVCON is not implemented in AmberTools > 1.3 anymore")
  3174.                 self.printWarn("Setting sqm for antechamber")
  3175.                 qprog = 'sqm'
  3176.             elif qprog == 'mopac':
  3177.                 if not os.path.exists(os.path.join(os.path.dirname(self.acExe), qprog)):
  3178.                     self.printWarn("MOPAC is not installed")
  3179.                     self.printWarn("Setting sqm for antechamber")
  3180.                     return None
  3181.                     qprog = 'sqm'
  3182.         else:
  3183.             self.printWarn("Old version of antechamber. Strongly consider upgrading to AmberTools")
  3184.             self.printWarn("Setting mopac for antechamber")
  3185.             qprog = 'mopac'
  3186.         self.qFlag = qDict[qprog]
  3187.         self.outTopols = [outTopol]
  3188.         if outTopol == 'all':
  3189.             self.outTopols = outTopols
  3190.         self.acParDict = {'base': base, 'ext': ext[1:], 'acBase': acBase,
  3191.                           'acMol2FileName': acMol2FileName, 'res': self.resName,
  3192.                           'leapAmberFile': leapAmberFile, 'baseOrg': self.baseOriginal,
  3193.                           'leapGaffFile': leapGaffFile}
  3194.  
  3195.  
  3196. class MolTopol(ACTopol):
  3197.  
  3198.     """"
  3199.         Class to write topologies and parameters files for several applications
  3201.         http://amber.scripps.edu/formats.html (not updated to amber 10 yet)
  3203.         Parser, take information in AC xyz and top files and convert to objects
  3205.         INPUTS: acFileXyz and acFileTop
  3206.         RETURN: molTopol obj or None
  3207.     """
  3208.  
  3209.     def __init__(self, acTopolObj=None, acFileXyz=None, acFileTop=None,
  3210.                  debug=False, basename=None, verbose=True, gmx4=False,
  3211.                  disam=False, direct=False, is_sorted=False, chiral=False):
  3212.  
  3213.         self.chiral = chiral
  3214.         self.obchiralExe = _getoutput('which obchiral') or ''
  3215.         self.allhdg = False
  3216.         self.debug = debug
  3217.         self.gmx4 = gmx4
  3218.         self.disam = disam
  3219.         self.direct = direct
  3220.         self.sorted = is_sorted
  3221.         self.verbose = verbose
  3222.         self.inputFile = acFileTop
  3223.         if acTopolObj:
  3224.             if not acFileXyz:
  3225.                 acFileXyz = acTopolObj.acXyzFileName
  3226.             if not acFileTop:
  3227.                 acFileTop = acTopolObj.acTopFileName
  3228.             self._parent = acTopolObj
  3229.             self.allhdg = self._parent.allhdg
  3230.             self.debug = self._parent.debug
  3231.             self.inputFile = self._parent.inputFile
  3232.         if not os.path.exists(acFileXyz) and not os.path.exists(acFileTop):
  3233.             self.printError("Files '%s' and '%s' don't exist")
  3234.             self.printError("molTopol object won't be created")
  3235.             return None
  3236.  
  3237. #         if not os.path.exists(self.obchiralExe) and self.chiral:
  3238. #             self.printError("no 'obchiral' executable, it won't work to store non-planar improper dihedrals!")
  3239. #             self.printWarn("Consider installing http://openbabel.org")
  3240.  
  3241.         self.xyzFileData = open(acFileXyz, 'r').readlines()
  3242.         self.topFileData = open(acFileTop, 'r').readlines()
  3243.         self.printDebug("prmtop and inpcrd files loaded")
  3244.  
  3245. #        self.pointers = self.getFlagData('POINTERS')
  3246.  
  3247.         self.getResidueLabel()
  3248.         if len(self.residueLabel) > 1:
  3249.             self.baseName = basename or os.path.splitext(os.path.basename(acFileTop))[0]  # 'solute'
  3250.         else:
  3251.             self.baseName = basename or self.residueLabel[0]  # 3 caps letters
  3252.         if acTopolObj:
  3253.             self.baseName = basename or acTopolObj.baseName
  3254.         self.printDebug("basename defined = '%s'" % self.baseName)
  3255.  
  3256.         self.getAtoms()
  3257.  
  3258.         self.getBonds()
  3259.  
  3260.         self.getAngles()
  3261.  
  3262.         self.getDihedrals()
  3263.  
  3264.         self.getChirals()
  3265.         if not os.path.exists(self.obchiralExe) and self.chiral:
  3266.             self.printError("no 'obchiral' executable, it won't work to store non-planar improper dihedrals!")
  3267.             self.printWarn("Consider installing http://openbabel.org")
  3268.         elif self.chiral and not self.chiralGroups:
  3269.             self.printWarn("No chiral atoms found")
  3270.  
  3271.         # self.setAtomPairs()
  3272.  
  3273.         # self.getExcludedAtoms()
  3274.  
  3275.         # a list of FLAGS from acTopFile that matter
  3276. #        self.flags = ( 'POINTERS', 'ATOM_NAME', 'CHARGE', 'MASS', 'ATOM_TYPE_INDEX',
  3277. #                  'NUMBER_EXCLUDED_ATOMS', 'NONBONDED_PARM_INDEX',
  3278. #                  'RESIDUE_LABEL', 'BOND_FORCE_CONSTANT', 'BOND_EQUIL_VALUE',
  3279. #                  'ANGLE_FORCE_CONSTANT', 'ANGLE_EQUIL_VALUE',
  3280. #                  'DIHEDRAL_FORCE_CONSTANT', 'DIHEDRAL_PERIODICITY',
  3281. #                  'DIHEDRAL_PHASE', 'AMBER_ATOM_TYPE' )
  3282.  
  3283.         # Sort atoms for gromacs output. # JDC
  3284.         if self.sorted:
  3285.             self.printMess("Sorting atoms for gromacs ordering.\n")
  3286.             self.sortAtomsForGromacs()
  3287.  
  3288.  
  3289. class Atom(object):
  3290.  
  3291.     """
  3292.         Charges in prmtop file has to be divide by 18.2223 to convert to charge
  3293.         in units of the electron charge.
  3294.         To convert ACOEF and BCOEF to r0 (Ang.) and epsilon (kcal/mol), as seen
  3295.         in gaff.dat for example; same atom type (i = j):
  3296.             r0 = 1/2 * (2 * ACOEF/BCOEF)^(1/6)
  3297.             epsilon = 1/(4 * A) * BCOEF^2
  3298.         To convert r0 and epsilon to ACOEF and BCOEF
  3299.             ACOEF = sqrt(ep_i * ep_j) * (r0_i + r0_j)^12
  3300.             BCOEF = 2 * sqrt(ep_i * ep_j) * (r0_i + r0_j)^6
  3301.                   = 2 * ACOEF/(r0_i + r0_j)^6
  3302.         where index i and j for atom types.
  3303.         Coord is given in Ang. and mass in Atomic Mass Unit.
  3304.     """
  3305.  
  3306.     def __init__(self, atomName, atomType, id_, resid, mass, charge, coord):
  3307.         self.atomName = atomName
  3308.         self.atomType = atomType
  3309.         self.id = id_
  3310.         self.cgnr = id_
  3311.         self.resid = resid
  3312.         self.mass = mass
  3313.         self.charge = charge  # / qConv
  3314.         self.coords = coord
  3315.  
  3316.     def __str__(self):
  3317.         return '<Atom id=%s, name=%s, %s>' % (self.id, self.atomName, self.atomType)
  3318.  
  3319.     def __repr__(self):
  3320.         return '<Atom id=%s, name=%s, %s>' % (self.id, self.atomName, self.atomType)
  3321.  
  3322.  
  3323. class AtomType(object):
  3324.  
  3325.     """
  3326.         AtomType per atom in gaff or amber.
  3327.     """
  3328.  
  3329.     def __init__(self, atomTypeName, mass, ACOEF, BCOEF):
  3330.         self.atomTypeName = atomTypeName
  3331.         self.mass = mass
  3332.         self.ACOEF = ACOEF
  3333.         self.BCOEF = BCOEF
  3334.  
  3335.     def __str__(self):
  3336.         return '<AtomType=%s>' % self.atomTypeName
  3337.  
  3338.     def __repr__(self):
  3339.         return '<AtomType=%s>' % self.atomTypeName
  3340.  
  3341.  
  3342. class Bond(object):
  3343.  
  3344.     """
  3345.         attributes: pair of Atoms, spring constant (kcal/mol), dist. eq. (Ang)
  3346.     """
  3347.  
  3348.     def __init__(self, atoms, kBond, rEq):
  3349.         self.atoms = atoms
  3350.         self.kBond = kBond
  3351.         self.rEq = rEq
  3352.  
  3353.     def __str__(self):
  3354.         return '<%s, r=%s>' % (self.atoms, self.rEq)
  3355.  
  3356.     def __repr__(self):
  3357.         return '<%s, r=%s>' % (self.atoms, self.rEq)
  3358.  
  3359.  
  3360. class Angle(object):
  3361.  
  3362.     """
  3363.         attributes: 3 Atoms, spring constant (kcal/mol/rad^2), angle eq. (rad)
  3364.     """
  3365.  
  3366.     def __init__(self, atoms, kTheta, thetaEq):
  3367.         self.atoms = atoms
  3368.         self.kTheta = kTheta
  3369.         self.thetaEq = thetaEq  # rad, to convert to degree: thetaEq * 180/Pi
  3370.  
  3371.     def __str__(self):
  3372.         return '<%s, ang=%.2f>' % (self.atoms, self.thetaEq * 180 / Pi)
  3373.  
  3374.     def __repr__(self):
  3375.         return '<%s, ang=%.2f>' % (self.atoms, self.thetaEq * 180 / Pi)
  3376.  
  3377.  
  3378. class Dihedral(object):
  3379.  
  3380.     """
  3381.         attributes: 4 Atoms, spring constant (kcal/mol), periodicity,
  3382.         phase (rad)
  3383.     """
  3384.  
  3385.     def __init__(self, atoms, kPhi, period, phase):
  3386.         self.atoms = atoms
  3387.         self.kPhi = kPhi
  3388.         self.period = period
  3389.         self.phase = phase  # rad, to convert to degree: kPhi * 180/Pi
  3390.  
  3391.     def __str__(self):
  3392.         return '<%s, ang=%.2f>' % (self.atoms, self.phase * 180 / Pi)
  3393.  
  3394.     def __repr__(self):
  3395.         return '<%s, ang=%.2f>' % (self.atoms, self.phase * 180 / Pi)
  3396.  
  3397. if __name__ == '__main__':
  3398.     t0 = time.time()
  3399.     print(header)
  3400.  
  3401.     parser = optparse.OptionParser(usage=usage + epilog)
  3402.     parser.add_option('-i', '--input',
  3403.                       action="store",
  3404.                       dest='input',
  3405.                       help="input file name with either extension '.pdb', '.mdl' or '.mol2' (mandatory if -p and -x not set)",)
  3406.     parser.add_option('-b', '--basename',
  3407.                       action="store",
  3408.                       dest='basename',
  3409.                       help='a basename for the project (folder and output files)',)
  3410.     parser.add_option('-x', '--inpcrd',
  3411.                       action="store",
  3412.                       dest='inpcrd',
  3413.                       help="amber inpcrd file name (always used with -p)",)
  3414.     parser.add_option('-p', '--prmtop',
  3415.                       action="store",
  3416.                       dest='prmtop',
  3417.                       help="amber prmtop file name (always used with -x)",)
  3418.     parser.add_option('-c', '--charge_method',
  3419.                       type='choice',
  3420.                       choices=['gas', 'bcc', 'user'],
  3421.                       action="store",
  3422.                       default='bcc',
  3423.                       dest='charge_method',
  3424.                       help="charge method: gas, bcc (default), user (user's charges in mol2 file)",)
  3425.     parser.add_option('-n', '--net_charge',
  3426.                       action="store",
  3427.                       type='int',
  3428.                       default=0,
  3429.                       dest='net_charge',
  3430.                       help="net molecular charge (int), for gas default is 0",)
  3431.     parser.add_option('-m', '--multiplicity',
  3432.                       action="store",
  3433.                       type='int',
  3434.                       default=1,
  3435.                       dest='multiplicity',
  3436.                       help="multiplicity (2S+1), default is 1",)
  3437.     parser.add_option('-a', '--atom_type',
  3438.                       type='choice',
  3439.                       choices=['gaff', 'amber', 'gaff2', 'amber2'],  # , 'bcc', 'sybyl']
  3440.                       action="store",
  3441.                       default='gaff',
  3442.                       dest='atom_type',
  3443.                       help="atom type, can be gaff, gaff2, amber (AMBER14SB) or amber2 (AMBER14SB + GAFF2), default is gaff",)
  3444.     parser.add_option('-q', '--qprog',
  3445.                       type='choice',
  3446.                       choices=['mopac', 'sqm', 'divcon'],
  3447.                       action="store",
  3448.                       default='sqm',
  3449.                       dest='qprog',
  3450.                       help="am1-bcc flag, sqm (default), divcon, mopac",)
  3451.     parser.add_option('-k', '--keyword',
  3452.                       action="store",
  3453.                       dest='keyword',
  3454.                       help="mopac or sqm keyword, inside quotes",)
  3455.     parser.add_option('-f', '--force',
  3456.                       action="store_true",
  3457.                       dest='force',
  3458.                       help='force topologies recalculation anew',)
  3459.     parser.add_option('-d', '--debug',
  3460.                       action="store_true",
  3461.                       dest='debug',
  3462.                       help='for debugging purposes, keep any temporary file created',)
  3463.     parser.add_option('-o', '--outtop',
  3464.                       type='choice',
  3465.                       choices=['all'] + outTopols,
  3466.                       action="store",
  3467.                       default='all',
  3468.                       dest='outtop',
  3469.                       help="output topologies: all (default), gmx, cns or charmm",)
  3470.     parser.add_option('-z', '--gmx4',
  3471.                       action="store_true",
  3472.                       dest='gmx4',
  3473.                       help='write RB dihedrals old GMX 4.0',)
  3474.     parser.add_option('-t', '--cnstop',
  3475.                       action="store_true",
  3476.                       dest='cnstop',
  3477.                       help='write CNS topology with allhdg-like parameters (experimental)',)
  3478.     parser.add_option('-e', '--engine',
  3479.                       type='choice',
  3480.                       choices=['tleap', 'sleap'],
  3481.                       action="store",
  3482.                       default='tleap',
  3483.                       dest='engine',
  3484.                       help="engine: tleap (default) or sleap (not fully matured)",)
  3485.     parser.add_option('-s', '--max_time',
  3486.                       action="store",
  3487.                       type='int',
  3488.                       default=36000,
  3489.                       dest='max_time',
  3490.                       help="max time (in sec) tolerance for sqm/mopac, default is 10 hours",)
  3491.     parser.add_option('-y', '--ipython',
  3492.                       action="store_true",
  3493.                       dest='ipython',
  3494.                       help='start iPython interpreter',)
  3495.     parser.add_option('-w', '--verboseless',
  3496.                       action="store_false",
  3497.                       default=True,
  3498.                       dest='verboseless',
  3499.                       help='print nothing',)
  3500.     parser.add_option('-g', '--disambiguate',
  3501.                       action="store_true",
  3502.                       dest='disambiguate',
  3503.                       help='disambiguate lower and uppercase atomtypes in GMX top file',)
  3504.     parser.add_option('-u', '--direct',
  3505.                       action="store_true",
  3506.                       dest='direct',
  3507.                       help="for 'amb2gmx' mode, does a direct conversion, for any solvent",)
  3508.     parser.add_option('-l', '--sorted',
  3509.                       action="store_true",
  3510.                       dest='sorted',
  3511.                       help="sort atoms for GMX ordering",)
  3512.     parser.add_option('-j', '--chiral',
  3513.                       action="store_true",
  3514.                       dest='chiral',
  3515.                       help="create improper dihedral parameters for chiral atoms in CNS",)
  3516.  
  3517.     options, remainder = parser.parse_args()
  3518.  
  3519.     amb2gmx = False
  3520.  
  3521. #     if options.chiral:
  3522. #         options.cnstop = True
  3523.  
  3524.     if not options.input:
  3525.         amb2gmx = True
  3526.         if not options.inpcrd or not options.prmtop:
  3527.             parser.error("missing input files")
  3528.     elif options.inpcrd or options.prmtop:
  3529.         parser.error("either '-i' or ('-p', '-x'), but not both")
  3530.  
  3531.     if options.debug:
  3532.         text = "Python Version %s" % verNum
  3533.         print('DEBUG: %s' % text)
  3534.  
  3535.     if options.direct and not amb2gmx:
  3536.         parser.error("option -u is only meaningful in 'amb2gmx' mode")
  3537.  
  3538.     try:
  3539.         if amb2gmx:
  3540.             print("Converting Amber input files to Gromacs ...")
  3541.             system = MolTopol(acFileXyz=options.inpcrd, acFileTop=options.prmtop,
  3542.                               debug=options.debug, basename=options.basename,
  3543.                               verbose=options.verboseless, gmx4=options.gmx4,
  3544.                               disam=options.disambiguate, direct=options.direct,
  3545.                               is_sorted=options.sorted, chiral=options.chiral)
  3546.             system.printDebug("prmtop and inpcrd files parsed")
  3547.             system.writeGromacsTopolFiles(amb2gmx=True)
  3548.         else:
  3549.             molecule = ACTopol(options.input, chargeType=options.charge_method,
  3550.                                chargeVal=options.net_charge, debug=options.debug,
  3551.                                multiplicity=options.multiplicity, atomType=options.atom_type,
  3552.                                force=options.force, outTopol=options.outtop,
  3553.                                engine=options.engine, allhdg=options.cnstop,
  3554.                                basename=options.basename, timeTol=options.max_time,
  3555.                                qprog=options.qprog, ekFlag='''"%s"''' % options.keyword,
  3556.                                verbose=options.verboseless, gmx4=options.gmx4,
  3557.                                disam=options.disambiguate, direct=options.direct,
  3558.                                is_sorted=options.sorted, chiral=options.chiral)
  3559.  
  3560.             if not molecule.acExe:
  3561.                 molecule.printError("no 'antechamber' executable... aborting ! ")
  3562.                 hint1 = "HINT1: is 'AMBERHOME' or 'ACHOME' environment variable set?"
  3563.                 hint2 = "HINT2: is 'antechamber' in your $PATH?\n    What 'which antechamber' in your terminal says?\n    'alias' doesn't work for ACPYPE."
  3564.                 molecule.printMess(hint1)
  3565.                 molecule.printMess(hint2)
  3566.                 sys.exit(1)
  3567.  
  3568.             molecule.createACTopol()
  3569.             molecule.createMolTopol()
  3570.         acpypeFailed = False
  3571.     except:
  3572.         exceptionType, exceptionValue, exceptionTraceback = sys.exc_info()
  3573.         print("ACPYPE FAILED: %s" % exceptionValue)
  3574.         if options.debug:
  3575.             traceback.print_tb(exceptionTraceback, file=sys.stdout)
  3576.         acpypeFailed = True
  3577.  
  3578.     execTime = int(round(time.time() - t0))
  3579.     if execTime == 0:
  3580.         msg = "less than a second"
  3581.     else:
  3582.         msg = elapsedTime(execTime)
  3583.     print("Total time of execution: %s" % msg)
  3584.  
  3585.     if options.ipython:
  3586.         try:
  3587.             from IPython.Shell import IPShellEmbed  # @UnresolvedImport @UnusedImport
  3588.         except:
  3589.             from IPython.frontend.terminal.embed import InteractiveShellEmbed as IPShellEmbed  # @UnresolvedImport @Reimport
  3590.         ipshell = IPShellEmbed()
  3591.         ipshell()
  3592.  
  3593.     try:
  3594.         rmtree(molecule.tmpDir)
  3595.     except:
  3596.         pass
  3597.     if acpypeFailed:
  3598.         sys.exit(1)
  3599.     try:
  3600.         os.chdir(molecule.rootDir)
  3601.     except:
  3602.         pass

 

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/黑客灵魂/article/detail/903151
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号