xplor.parseArguments()

seq= " ACE ASP ALA ALA DAP ALA ALA "

import psfGen
psfGen.residueTypes['protein'].append("DAP")

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


xplor.command("""

topology
  residue DAP
  
  group
   atom N type=NH1 charge=-0.360 end
   atom HN type=H charge= 0.260 end
  group
   atom CA type=CT charge= 0.000 end
   atom HA type=HA charge= 0.100 end
  group
   atom CB type=CT charge=-0.200 end
   atom HB1 type=HA charge= 0.100 end
   atom HB2 type=HA charge= 0.100 end
  group
   atom NZ type=NH3 charge=-0.810 end
   atom HZ type=HC charge= 0.435 end
  group
   atom C type=C charge= 0.480 end
   atom O type=O charge=-0.480 end
  
   bond N HN
   bond N CA bond CA HA
   bond CA CB bond CB HB1 bond CB HB2
   bond CB NZ bond NZ HZ
   bond CA C
   bond C O
  
   improper HA N C CB !stereo CA
   improper HB1 HB2 CA NZ !stereo CB
  
   dihedral NZ CB CA N
  
  end
  
  presidue CLOS
  
   modify    atom -NZ   type=NH1 end
   modify    atom -HZ   type=H end

   delete atom +OD2 end
  
   add bond -NZ +CG
  
   add angle +CB +CG -NZ
!   add angle +OD2 +CG -NZ
   add angle +CG -NZ -CB
   add angle +CG -NZ -HZ
  
  end
end
""")

psfGen.seqToPSF(seq, amidate_cterm=True)


xplor.command("""

 patch CLOS reference=-=(resid 5) reference=+=(resid 2) end
""")

# in my last email, I forgot I need change ASP from L to D. 
psfGen.dAmino(resid=2)


xplor.command("write psf output=isodd.psf end")

import protocol
protocol.initRandomSeed()
protocol.addUnknownAtoms(dyn_stepsize=0.012,dyn_numStepMul=2)

for i in range(4):
 try:
     protocol.fixupCovalentGeom(useVDW=1,maxIters=400)
     break
 except Exception, e:
     if e.args[0].startswith("Covalent geometry still violated"):
         continue
     raise
 pass

from pdbTool import PDBTool
PDBTool("isodd.pdb").write()