* This file is mc8.inp.
* For Monte Carlo movement of domains of AraC.
* Input is coordinates of dimerization domains and
*   two strands of DNA and DBD positioned so as to
*   be in a possible position wrt dimerization domain.
* Outputs pdb file with coordinates of both domains
*  of the two subunits and the four DNA strands.
* Header of pdb contains final energy and number
*  of translations, rotations, and Metropolis jumps.
* Useage charmm seed=xy link=value1 HH=value2
*    HT=value3 <mc8.inp >mc8.out
* where xy is seed for random number generator,
*  link is linker distance between domains, 10-20 Angstroms
*  HH is head-to-head distance, typically 50 to 60
*  and HT is head-to-tail, and tail-to-tail distance.
*

!***********************************************************************!
!Initialize, read in the files and molecules, create second DBD + DNA

! Open and read modified rtf.
open read card unit 20 name "top_all27_prot_dna.rtf" 
read rtf card unit 20
close unit 20

! Open and read parameter file.
open read card unit 20 name "par_all27_prot_na.prm" 
read param card unit 20
close unit 20

! Read sequence from pdb file and generate segment
open read unit 12 card name "dna1a.pdb" 
read sequence pdb unit 12
generate dnaa setup first 5ter last 3ter 
rewind unit 12

open read unit 13 card name "dna2a.pdb"
read sequence pdb unit 13
generate dnab setup first 5ter last 3ter
rewind unit 13

open read unit 14 card name "prota.pdb"
read sequence pdb unit 14
generate dbd1 setup
rewind unit 14

open read unit 21 card name chaina1.pdb 
read sequence pdb unit 21
generate dim1 setup 
rewind unit 21

open read unit 22 card name chainb1.pdb 
read sequence pdb unit 22
generate dim2 setup
rewind unit 22

! Read coordinates and close the units
read coordinate pdb offset -400 unit 12
close unit 12

read coordinate pdb append offset -424 unit 13
close unit 13

read coordinate pdb append offset -8 unit 14
close unit 14

read coordinate pdb append offset -6 unit 21
close unit 21

read coordinate pdb append offset -6 unit 22
close unit 22

!Generate missing atoms
ic fill preserve
ic parameter
ic build
hbuild


!Create second DBD + DNA
generate dnac duplicate dnaa setup
generate dnad duplicate dnab setup
generate dbd2 duplicate dbd1 setup

!Put copy of first DBD domain + DNA in comparison coordinate set as initial reference
coordinate copy comp select segid dnaa .or. segid dnab .or. segid dbd1 end

system -
"echo Finished initial setup > /dev/tty"

!*********************************************************************!

!Set up the main outer loop
set j 0
label outloop

!Set up random number generation
calculate ijseed @seed * (@j + 3141)
rand uniform scale 200 offset -100 iseed @ijseed

!*********************************************************************!

!Begining of coarse-grained Monte Carlo loop

!Initialize operations counters
set rotops 0
set transops 0
set mcops 0

set i 0
label ctoploop

!Beginning translation step
!Obtain coordinates of first DNA binding domain + DNA
coordinate copy select segid dnaa .or. segid dnab .or. segid dbd1 end

!Perform random translation on first DNA binding domain + DNA
coordinate translate XDIR ?rand YDIR ?rand ZDIR ?rand -
  select segid dnaa .or. segid dnab .or. segid dbd1 end

!Setup to goto and return from from routine for generating symm subunit and noe constraints
set returnpoint ctranslation

!Go to the dupnoe routine (duplicate symmetrical subunit and set noe constraints)
goto dupnoe

!Come back from routine
label ctranslation

!Check for having lower energy
update
gete
set newener ?ENER

!First cycle special treatment
if @I LE 0 set refener ?ENER

calculate diff @refener - @newener

!The Metropolis section

!If new coordinates yield higher energy, jump to higher section, otherwise do lower stuff
if diff LT 0 goto chighera   
coordinate copy comp select segid dnaa .or. segid dnab .or. segid dbd1 end
set refener ?ENER
increment transops
goto cproceeda

label chighera

!Forget it if energy difference is large
if @diff LT -5 goto cproceeda

calculate boltz exp ( (-1 * @DIFF) / 0.6 )
calculate randone ( ( ?RAND + 1 ) * 0.5 ) 
if @randone GT @boltz goto cproceeda
set refener ?ENER
increment mcops
coordinate copy comp select segid dnaa .or. segid dnab .or. segid dbd1 end

label cproceeda

system -
"echo Loop @I coarse tr: current reference energy @REFENER energy ?ENER > /dev/tty"

!Finished translation step and keeping or not keeping of new coordinates

!***********************************************************************!

!Begin rotation step of coarse-grained Monte Carlo

!Obtain coordinates of first domain
coordinate copy select segid dnaa .or. segid dnab .or. segid dbd1 end

!Perform random rotation
set xrot ?rand 
set yrot ?rand 
set zrot ?rand 
set angle ?rand
coordinate statistics select segid dnaa .or. segid dnab .or. segid dbd1 end 

coordinate rotate -
  XDIR @xrot YDIR @yrot ZDIR @zrot -
  XCEN ?XAVE  YCEN ?YAVE ZCEN  ?ZAVE -
  PHI @angle select segid dnaa .or. segid dnab .or. segid dbd1 end

!Setup to goto and return from from routine for generating symm subunit and noe constraints
set returnpoint crotation

!Go to dupnoe routine
goto dupnoe

!Return point from noe routine
label crotation

!Same old energy comparison stuff again
update
gete
set newener ?ENER

calculate diff @refener - @newener

!The Metropolis section

!If new coordinates yield higher energy, jump to higher section, otherwise do lower stuff
if diff LT 0 goto chigherb   
coordinate copy comp select segid dnaa .or. segid dnab .or. segid dbd1 end
set refener ?ENER
increment rotops
goto cproceedb

label chigherb

!Forget it if the energy difference is large, no chance of a Metropolis step
if @diff LT -5 goto cproceedb

calculate boltz exp ( (-1 * @DIFF) / 0.6 )
calculate randone ( ( ?RAND + 1 ) * 0.5 ) 
if @randone GT @boltz goto cproceedb
set refener ?ENER
increment mcops
coordinate copy comp select segid dnaa .or. segid dnab .or. segid dbd1 end

label cproceedb

system -
"echo Loop @I coarse rot: current reference energy @REFENER energy ?ENER > /dev/tty"

!End of MC loop
increment i
if @i lt 2000 goto ctoploop

!***********************************************************************!
!Begining of fine-grained Monte Carlo section

!Set up random number generation for the fine-grained loop
calculate ijseed @seed * (@j + 1)
rand uniform scale 1 offset -0.5 iseed @ijseed

!Initialize operations counters
set rotops 0
set transops 0
set mcops 0 

set i 0
label ftoploop

!Obtain coordinates of first DNA binding domain + DNA
coordinate copy select segid dnaa .or. segid dnab .or. segid dbd1 end

!Perform random translation on first DNA binding domain + DNA
coordinate translate XDIR ?rand YDIR ?rand ZDIR ?rand -
  select segid dnaa .or. segid dnab .or. segid dbd1 end

!Setup to goto and return from from routine for generating symm subunit and noe constraints
set returnpoint ftranslation

!Go to the dupnoe routine (duplicate symmetrical subunit and set noe constraints)
goto dupnoe

!Come back from routine
label ftranslation

!Check for having lower energy
update
gete
set newener ?ENER

!First cycle special treatment
if @I LE 0 set refener ?ENER

calculate diff @refener - @newener

!The Metropolis section

!If new coordinates yield higher energy, jump to higher section, otherwise do lower stuff
if diff LT 0 goto fhighera   
coordinate copy comp select segid dnaa .or. segid dnab .or. segid dbd1 end
set refener ?ENER
increment transops
goto fproceeda

label fhighera

!Forget it if energy difference is large
if @diff LT -5 goto fproceeda

calculate boltz exp ( (-1 * @DIFF) / 0.6 )
calculate randone ( ( ?RAND +.5 ) * 0.5 ) 
if @randone GT @boltz goto fproceeda
set refener ?ENER
increment mcops
coordinate copy comp select segid dnaa .or. segid dnab .or. segid dbd1 end

label fproceeda

system -
"echo Loop @I fine tr: current reference energy @REFENER energy ?ENER > /dev/tty"

!Finished translation step and keeping or not keeping of new coordinates

!******************************************************************!

!Begin rotation step of fine-grained Monte Carlo

!Obtain coordinates of first domain
coordinate copy select segid dnaa .or. segid dnab .or. segid dbd1 end

!Perform random rotation
set xrot ?rand 
set yrot ?rand 
set zrot ?rand 
set angle ?rand
calculate angle ( @angle ) 
coordinate statistics select segid dnaa .or. segid dnab .or. segid dbd1 end 
coordinate rotate -
  XDIR @xrot YDIR @yrot ZDIR @zrot -
  XCEN ?XAVE  YCEN ?YAVE ZCEN  ?ZAVE -
  PHI @angle select segid dnaa .or. segid dnab .or. segid dbd1 end

!Setup to goto and return from from routine for generating symm subunit and noe constraints
set returnpoint frotation

!Go to dupnoe routine
goto dupnoe

!Return point from noe routine
label frotation

!Same old energy comparison stuff again
update
gete
set newener ?ENER

calculate diff @refener - @newener

!The Metropolis section

!If new coordinates yield higher energy, jump to higher section, otherwise do lower stuff
if diff LT 0 goto fhigherb   
coordinate copy comp select segid dnaa .or. segid dnab .or. segid dbd1 end
set refener ?ENER
increment rotops
goto fproceedb

label fhigherb

!Forget it if the energy difference is large, no chance of a Metropolis step
if @diff LT -5 goto fproceedb

calculate boltz exp ( (-1 * @DIFF) / 0.6 )
calculate randone ( ( ?RAND + 1 ) * 0.5 ) 
if @randone GT @boltz goto fproceedb
set refener ?ENER
increment mcops
coordinate copy comp select segid dnaa .or. segid dnab .or. segid dbd1 end

label fproceedb

system -
"echo Loop @I fine rot: current reference energy @REFENER energy ?ENER > /dev/tty"

!End of MC loop
increment i
if @i lt 500 goto ftoploop

!*************************************************************************!
!Output coordinates of fine-grained MC search

!Get the reference coordinates
coordinate copy select segid dnaa .or. segid dnab .or. segid dbd1 end

!Generate symmetrical subunit
coordinate duplicate select segid dnaa end select segid dnac end
coordinate duplicate select segid dnab end select segid dnad end
coordinate duplicate select segid dbd1 end select segid dbd2 end

!rotate 180 about z
coordinate rotate matrix select segid dnac .or. segid dnad .or. segid dbd2 end
-1 0 0
0 -1 0
0 0 1

!Update energy so its value can be put in the output pdb file
update
gete

!Write one pdb with all coordinates and statistics of search
open write unit 15 card name fine@j.pdb
write coordinates pdb select all end unit 15
* Shifted by MC fine, Energy ?ENER , Trans @TRANSOPS, Rot @ROTOPS, Jumps @MCOPS
*  iseed @ijseed linker @link head-to-head @hh head-to-tail & tail-to-tail @ht
*
close unit 15

increment j
if @j LT 5 goto outloop

system -
"echo Normal end, system energy @REFENER  > /dev/tty"

stop

!****************************************************************************!
!****************************************************************************!
!Routine for generating symm. subunit and setting noe constraints

label dupnoe

coordinate duplicate select segid dnaa end select segid dnac end
coordinate duplicate select segid dnab end select segid dnad end
coordinate duplicate select segid dbd1 end select segid dbd2 end

!rotate 180 about z
coordinate rotate matrix select segid dnac .or. segid dnad .or. segid dbd2 end
-1 0 0
0 -1 0
0 0 1

!Set noe constraints
noe
    reset  

!Holding N-termini of the dbd's near the C-termini of dd's.
    assign sele segid dim1 .and. resid 168 .and. type ca end sele segid dbd1 .and. resid 9 .and. type ca end -
    kmin 50 rmin 5 kmax 50 rmax @link

    assign sele segid dim2 .and. resid 168 .and. type ca end sele segid dbd2 .and. resid 9 .and. type ca end -
    kmin 50 rmin 5 kmax 50 rmax @link   

!Head-to-head distance constraint
    assign sele segid dnac .and. resid 401 .and. type O4' end sele segid dnaa .and. resid 401 .and. type O4' end -
    kmin 50 rmin @HH kmax 50 rmax @HH

!Tail-to-tail distance constraint
    assign sele segid dnaa .and. resid 423 .and. type O4' end sele segid dnac .and. resid 423 .and. type O4' end -
    kmin 50 rmin @HT kmax 50 rmax @HT

!Head-to-tail distance constraints
    assign sele segid dnaa .and. resid 401 .and. type O4' end sele segid dnac .and. resid 423 .and. type O4' end -
    kmin 50 rmin @HT kmax 50 rmax @HT

    assign sele segid dnac .and. resid 401 .and. type O4' end sele segid dnaa .and. resid 423 .and. type O4' end -
    kmin 50 rmin @HT kmax 50 rmax @HT
    end

!Return to point where routine was called
goto @returnpoint

!**************************************************************!