*  Analysis of GCN4 association free energy from post-processing of
*  298 K trajectory of the protein dimer GCN4 leucine zipper from 2zta.pdb
*  Analysis uses the CHARMM param19 generalized Born module to estimate
*  the free energy for 2H => H:H.
*  This file performs the following functions:
*  1) reads parameter and topology files (toph19/param19)
*  2) reads trajectory a frame at a time and computes the approximate
**    binding energy
*  3) write a trajectories of the "binding" free energy vs time to a file
*

if @?N2A eq 0 set N2A = false

!  read in the topology and parameter files

open unit 1 read form name "$toppar/toph19.inp"
read rtf card unit 1
close unit 1

open unit 1 read form name "$toppar/param19.inp"
read param card unit 1
close unit 1

!  Now read in and generate the sequence.
read sequence card
*  Sequence for GCN4 Leucine Zipper coiled-coil dimer
*
  35
ACE ARG MET LYS GLN LEU GLU ASP LYS VAL GLU GLU LEU
LEU SER LYS ASN TYR HIS LEU GLU ASN GLU VAL ALA ARG
LEU LYS LYS LEU VAL GLY GLU ARG CBX

generate zipa setup  !use setup to build internal coordinate tables

!  Now read in and generate the sequence.
read sequence card
*  Sequence for GCN4 Leucine Zipper coiled-coil dimer
*
  35
ACE ARG MET LYS GLN LEU GLU ASP LYS VAL GLU GLU LEU
LEU SER LYS ASN TYR HIS LEU GLU ASN GLU VAL ALA ARG
LEU LYS LYS LEU VAL GLY GLU ARG CBX

generate zipb setup  !use setup to build internal coordinate tables

if N2A ne true open unit 30 write form name 2zta_bind.dat
if N2A eq true open unit 30 write form name 2zta_bind_N2A.dat
write title unit 30
*# Time    Ecomplex   Ehelices  DEbind   DEcum
*

open unit 1 write unform name "crd/traj.dcd"
open unit 2 read form name "crd/2zta_d1-60p-sub.hex"
dynamics unformat input 2 output 1
close unit 2

open unit 1 read unform name "crd/traj.dcd"


traj iread 1 nread 1 skip 1200
set count = 1
set DEsum = 0
label read

   traj read

   if count ne 1 goto skipfirst
      !***Optimized GB parameters for version 22 protein parameters 
      !***w/ polar hydrogen radii set to 1.0 A
      define polarH select property radii .le. 1.0 end !  Select all hydrogen atoms w/
                                                 ! vdW radii .le. 1.0 A
      scalar wmain = radii
      scalar wmain set 1.0 select polarH end

      GBorn P1 0.4380 P2 0.1664 P3 0.0138 P4 8.84109 P5 1.0 Lambda 0.7160 -
            Epsilon 80.0 Weight
      if N2A eq true define ala select ( type c .or. type o .or. -
                            type n .or. type h .or. type ca .or. type cb ) end
      if N2A eq true define asns select ( (segid zipa .and. resid 17) -
                                   .and. .not. ala ) -
                                   .or. ( (segid zipb .and. resid 17) -
                                   .and. .not. ala ) end

      Update cutnb 999 switch vswitch

   label skipfirst

   if N2A ne true Energy inbfrq 0
   if N2A eq true Intere inbfrq 0 select .not. asns end
   Set E2 = ?ener

!  "Trick" to get energy of dissociated complex:  translate to great distance.
   coor trans xdir 1000 select segid zipa end

   if N2A ne true Energy inbfrq 0
   if N2A eq true Intere inbfrq 0 select .not. asns end
   Set E1 = ?ener

   Calc DE = @E2 - @E1
   Calc DEsum = @DEsum + @DE
   Calc DEcum = @DEsum / @count
   
   write title unit 30
* @count  @E2   @E1   @DE   @DEcum

   incr count by 1
if count le 500 goto read


system "rm crd/traj.crd"
stop