* This CHARMM input script demonstrates the use of constant force (cf)
* and constant speed (steered-molecular dynamics, smd) pulling
* experiments on an all beta protein E2lip3 (PDB: 1qjo.pdb).
* The experiments explore the differences in mechanical stability and
* mechanical stress induced unfolding of this protein as a function of
* where the attachment is made. The calculations utilize the sequence-dependent
* Go-model of Karanicolas and Brooks and mirror the earlier study of
* Brockwell et al. (NSB, 10, 731-738 (2003)).
* This input script was developed following a template provide by
* Emanuele Paci from the University of Zurich in July 2004.
* The only input to this script is whether one is running an constant-force
* pulling calculation (Method = cf) or a constant speed pulling calculation
* (Method = smd) (the default is smd), where the cantilever is attached
* (Last = 80/41), and the number of pulling trials (Ntrials = x).
**
**
* Required files: go_[at]prot.top, go_[at]prot.param, go_[at]prot.seq, 
*                 go_[at]prot.pdb
* where these files are from the MMTSB Go Model Server and the charmm
* variable prot (in [at]prot) is defined as prot = 1qjo
**
**
* Usage: charmm Method=[cf/smd] Last=[80/41] NTrials=[x] < pull_GoModel.inp > pull_GoModel.out
**
**
* Prepared by C.L. Brooks III, The Scripps Research Institute, July, 2004.
*

if @?Method eq 0  set Method = smd  !set default method to smd
if @?Last eq 0    set Last = 80     !default is to pull from end
if @?Ntrials eq 0 set Ntrials = 1   !default number of runs is 1

SET SegName  = prot     ! Name of the domain
SET prot     = 1qjo
SET Temp     = 300      ! temperature for heatbath during pulling
SET TimeStep = 0.010    ! with Go models we can use a large timestep
SET Bias     = 100      ! constant force in pN
SET Nsteps   = 100000   ! Number of timesteps
SET First    = 1
Calc Lmax = (@Last - @First) * 3.8 * 1.4 ! approximate maximum length of peptide
if Method eq smd SET Speed = 0.1       ! pulling speed with SMD (Angstroom/ps)

eten on  ! Turn on the special interaction term for Go models

! Stream in the Go model parameters
STREAM "GO_@PROT.top"
STREAM "GO_@PROT.param"
STREAM "GO_@PROT.seq"

set cnt = 1
label dopull


! Set up a random seed
system "echo 'set iseed = ' `date +%H%M%S` > seed.stream"
stream seed.stream
system "rm seed.stream"
calc iseed = 2* @iseed + 1

! Read in the protein coordinates
OPEN UNIT 1 READ FORM NAME "GO_@PROT.pdb"
READ COOR PDB UNIT 1
CLOSE UNIT 1

ENERGY

SCALAR FBETA SET 0.1 SELECT all END  ! Use friction coefficient to couple to Langevin
                                     ! heatbath.
SHAKE BOND TOL 1.0E-6 PARAM          ! Use constraint to fix bonds

COOR DIST SELE  (ATOM @SegName @First CA) END SELE (ATOM @SegName @Last CA) END
DEFINE nc SELE ((ATOM @SegName @First CA) .OR. (ATOM @SegName @Last CA)) END

OPEN WRITE FORMATTED UNIT 32 NAME @Method@Bias.@First-@Last_@cnt.rst
OPEN WRITE UNFORMATT UNIT 34 NAME @Method@Bias.@First-@Last_@cnt.dcd
OPEN WRITE FORMATT   UNIT 44 name @Method@Bias.@First-@Last_@cnt.juj

! { Pull with constant force }
if Method eq cf AFM @Method IUNJUJ 44 ALPHA @Bias SELE nc END XIMAX @Lmax
! { Pull with steered molecular dynamics }
if Method eq smd AFM @Method IUNJUJ 44 ALPHA @Bias BETA @Speed SELE nc END XIMAX @Lmax

DYNAMICS LEAP LANGEVIN START TIMESTEP @Timestep -
   NSTEP @Nsteps NPRINT 100 IPRFRQ 50000 -
   FIRSTT @Temp FINALT @Temp TWINDL -10.0 TWINDH 10.0 -
   ICHECW 1 IEQFRQ 500 IASORS 1 IASVEL 1 -
   NSAVC 1000 NSAVV 0 IUNVEL -1 -
   IUNREA 31 IUNWRI 32 KUNIT -1 IUNCRD 34 - 
   INBFRQ 0 ILBFRQ 0 TBATH @Temp ECHECK 9999999999.0 ISEED @iseed

OPEN UNIT 1 WRITE FORM NAME @Method@Bias.@First-@Last_@cnt.pdb
WRITE COOR PDB UNIT 1
CLOSE UNIT 1

afm reset ! reset afm module for next trip through code

incr cnt by 1
if cnt le @Ntrials goto dopull

STOP