#!/usr/bin/env python
# -*- coding: iso-8859-1 -*-
from __future__ import print_function, division
import  sys
if( sys.version_info[0] == 2 ):
    range = range

import  qm3.mol

m = qm3.mol.molecule( "pdb" )


HSD = [ ("A",100) ]
HSE = [ ("A",200) ]
HSP = [ ("A",300) ]


for i,j in HSD:
    for k in m.indx[i][j].values():
        m.resn[k] = "HID"

for i,j in HSE:
    for k in m.indx[i][j].values():
        m.resn[k] = "HIE"

for i,j in HSP:
    for k in m.indx[i][j].values():
        m.resn[k] = "HIP"

for i in range( len( m.res_lim ) - 1 ):
    if( m.resn[m.res_lim[i]] == "ILE" ):
        try:
            m.labl[m.indx[m.segn[m.res_lim[i]]][m.resi[m.res_lim[i]]]["CD"]] = "CD1"
        except:
            pass

f = file( "tleap.inp", "wt" )
f.write( """source leaprc.protein.ff14SB
source leaprc.gaff
source leaprc.water.tip3p
#loadamberprep custom.prepi
#loadamberparams custom.frcmod
""" )

t = "mol = combine {"

for i in range( len( m.seg_lim ) - 1 ):
    w = m.segn[m.res_lim[m.seg_lim[i]]]
    f.write( "s_%s = {"%( w ) )
    for j in range( m.seg_lim[i], m.seg_lim[i+1] ):
        f.write( " %s"%( m.resn[m.res_lim[j]] ) )
        if( ( j + 1 ) % 50 == 0 ):
            f.write( "\n" )
    f.write( " }\n" )
    m.pdb_write( "tleap_pdb.%s"%( w ), sele = range( m.res_lim[m.seg_lim[i]], m.res_lim[m.seg_lim[i+1]] ) )
    f.write( "p_%s = loadpdbusingseq tleap_pdb.%s s_%s\n"%( w, w, w ) )
    t += " p_%s"%( w )

f.write( t + " }\n" )

f.write( """
#mol = loadpdb pdb

solvatebox mol TIP3PBOX 14
#addions mol Cl- 0
#addions mol Na+ 0
saveamberparm mol prmtop inpcrd
savepdb mol pdb
quit """ )
f.close()
print( "\n\ntleap -s -f tleap.inp\n\n" )
