#
# more complicated addAtoms example in which a subset of coordinates
# are fixed.
#
#
xplor.requireVersion("2.12")

successFile="addAtoms2.success" # written if script runs ok

import os
try:
    os.unlink(successFile)
except:
    pass

xplor.parseArguments() # check for typos on the command-line

seq="""
MET SER ASP ASN GLY GLU LEU GLU ASP LYS PRO PRO ALA PRO PRO VAL  
GLU HIS THR ILE HIS VAL GLY PHE ASP ALA VAL THR GLY GLU PHE THR
GLY MET PRO GLU GLN TRP ALA ARG LEU LEU GLN THR SER ASN ILE THR
LYS LEU GLU GLN LYS LYS ASN PRO GLN ALA VAL LEU ASP VAL LEU LYS
PHE TYR ASP SER ASN THR VAL LYS GLN LYS TYR LEU SER PHE THR PRO
PRO GLU LYS ASP GLY PHE PRO SER GLY THR PRO ALA LEU ASN ALA LYS
GLY THR GLU ALA PRO ALA VAL VAL THR GLU GLU GLU ASP ASP ASP GLU
GLU THR ALA PRO PRO VAL ILE ALA PRO ARG PRO ASP HIS THR LYS SER
ILE TYR THR ARG SER VAL ILE ASP PRO VAL PRO ALA PRO VAL GLY ASP
SER HIS VAL ASP GLY ALA ALA LYS SER LEU ASP LYS GLN LYS LYS LYS
THR LYS MET THR ASP GLU GLU ILE MET GLU LYS LEU ARG THR ILE VAL
SER ILE GLY ASP PRO LYS LYS LYS TYR THR ARG TYR GLU LYS ILE GLY
GLN GLY ALA SER GLY THR VAL PHE THR ALA THR ASP VAL ALA LEU GLY
LEU GLN ALA LEU GLU PHE LEU HIS ALA ASN GLN VAL ILE HIS ARG ASP
ILE LYS SER ASP ASN VAL LEU LEU GLY MET GLU GLY SER VAL LYS LEU
THR ASP PHE GLY PHE CYS ALA GLN ILE THR PRO GLU GLN SER LYS ARG
SER THR MET VAL GLY THR PRO TYR TRP MET ALA PRO GLU VAL VAL THR
ARG LYS ALA TYR GLY PRO LYS VAL ASP ILE TRP SER LEU GLY ILE MET
ALA ILE GLU MET VAL GLU GLY GLU PRO PRO TYR LEU ASN GLU ASN PRO
LEU ARG ALA LEU TYR LEU ILE ALA THR ASN GLY THR PRO GLU LEU GLN
ASN PRO GLU LYS LEU SER PRO ILE PHE ARG ASP PHE LEU ASN ARG CYS
LEU GLU MET ASP VAL GLU LYS ARG GLY SER ALA LYS GLU LEU LEU GLN
HIS PRO PHE LEU LYS LEU ALA LYS PRO LEU SER SER LEU THR PRO LEU
ILE MET ALA ALA LYS GLU ALA MET LYS SER ASN ARG
"""

simWorld.setRandomSeed(5521)

import protocol
protocol.initTopology('protein')
protocol.initParams('protein')

import psfGen
psfGen.seqToPSF(seq,seqType='prot',startResid=1)

# the input pdb has doesn't have MTSL (CYSP) tags - the sequence is
# different
#
protocol.initCoords('2A.pdb', maxUnreadEntries = None)

#erase all sidechain coords for mutated tag- to make certain they're
regenerated
#protocol.initCoords(erase=True,
#                  selection="resname CYSP and not (name C or name CA or name
#N)")

# atoms not in this selection will not be moved at all
#
#fixGeom=AtomSel('''(not (name ca or name n or name c or
#                 name o or name hn or name ha))
#                 and (resid 122 or resid 72 or resid 87 or
#                 resid 152 or resid 161 or resid 163)''')
fixGeom=AtomSel('''( name C or name CA or name N ) and (resid 1:3 or resid
134:136''')

import protocol, regularize

# using this try...except clause will allow the script to write
# pdb and psf files even if there was an error in fixupCovalentGeom
#
#
try:
    regularize.addUnknownAtoms_new()
    protocol.fixupCovalentGeom(sel=fixGeom,verbose=2,useVDW=1)
    open(successFile,"w").write("ok")
except protocol.CovalentViolation:
    pass

# write out new psf and pdb files.
#
xplor.command("write psf output=addAtoms2.psf end")
protocol.writePDB("addAtoms2.pdb")