* This tutorial example illustrates how to set-up and build
* a peptide in CHARMM.
* Note first we will "regularize" the peptide and then solvate it
* using the MMTSB Tool Set tool convpdb.pl
* In an accompanying example we will do all this with three mmtsb commands.
*

! Read in the parameter and topology files.
! We'll use the $CHARMMDATADIR as our source for top/param files

open unit 1 read form name "$CHARMMDATA/top_all22_prot_cmap.inp"
read rtf card unit 1
close unit 1

open unit 1 read form name "$CHARMMDATA/par_all22_prot_cmap.inp"
read param card unit 1
close unit 1

! We can read the sequence from the pdb file we created
open unit 1 read form name 1rnu_cpep.pdb
read sequ pdb unit 1

! Or we can add it explicitly - here we use the former
!read sequ card unit 5
!* Sequence for 13 residue N-terminal fragment of RNAse A
!* C-peptide
!*
!13
!LYS GLU THR ALA ALA ALA LYS PHE GLU ARG GLN HSD MET

generate pro0 first none last ct3 setup

rewind unit 1
read coor pdb unit 1 resi  ! Read coordinates from pdb file using resid field
close unit 1

! Since the truncated peptide will not have the appropriate blocking groups, 
! we add them here. Use ic build and then hbuild.
ic param
ic build
coor init select hydrogen end
! For consistency with defaults in the mmtsb tool minCHARMM.pl, we use
! this hbuild set-up
hbuild atom cdie eps 80.0 cutnb 10.0 ctofnb 7.5 ctonnb 6.5 shift vshift bygr

! Now we will do a rough minimization to "regularize" the structure and to
! remove any bad contacts from crystal structure prior to solvating the peptide
! We will harmonically restrain all peptide heavy atoms to ensure that
! the conformation does not change significantly from its initial state.

! Copy current coordinates to comparison set to see the motion.
coor copy compare

constrain harmonic force 5 mass select .not. hydrogen end

mini sd nstep 100 rdie switch vswitch cutnb 12 ctonnb 8 ctofnb 11 inbfrq -1

coor orie rms select .not. hydrogen end

open unit 1 write form name 1rnu_cpep_min.pdb
write coor pdb unit 1


stop