* correl  
*

if @?task eq 0 goto deadformat

if @?pdbfile eq 0 goto iamdead

stream @pdbfile_dimens.stream

!*********************************************************
!     get parameters set for this run
!*********************************************************
!
   goto getparams
   label doneparams

!*********************************************************
!     get topology and force field for this run
!*********************************************************
!

!prnlev 0 node 0
wrnlev 0
open unit 1 read form name @toppar/top_all22_prot.inp
read rtf card unit 1
close unit 1

open unit 1 read form name @toppar/par_all22_prot.inp
read param card unit 1
close unit 1

prnlev 5 node 0
wrnlev 2


!*********************************************************
!     get psf for this system
!         This psf does not have the solvent
!*********************************************************
!
open unit 1 read form name @pdbfile.psf
read psf card unit 1
close unit 1


!*********************************************************
!           Add water
!*********************************************************
read sequ tip3 @nwat
generate wat setup noangl nodihe

goto @task

label atomrms
!*********************************************************
!             Atom Average rms
!*********************************************************
if @water eq true goto atrmswat
delete atom sele resn TIP3 end
set avgnam @pdbfile_e@begin-@end_nowat.avgcrd
open unit 30 unform read name @pdbfile_e@begin-@end_nowat.crd
open unit 31 form write name @avgnam
open unit 32 form write name @pdbfile_e@begin-@end_nowat.avgdmb
coor dyna comp firstu 30 nunit 1
write coor comp card unit 31
write coor comp dumb unit 32
close unit 30
close unit 31
close unit 32

goto extractrms
label doneextrms
stop

label atrmswat
set avgnam @pdbfile_e@begin-@end.avgcrd
open unit 30 unform read name @pdbfile_e@begin-@end.crd
open unit 31 form write name @avgnam
open unit 32 form write name @pdbfile_e@begin-@end.avgdmb
coor dyna comp firstu 30 nunit 1
write coor comp card unit 31
write coor comp dumb unit 32
close unit 30
close unit 31
close unit 32

goto extractrms
label doneextrms
stop

label timerms
!*********************************************************
!             time series of average rms
!*********************************************************
if @water eq true goto trmswat
delete atom sele resn TIP3 end
open unit 30 unform read name @pdbfile_e@begin-@end_nowat.crd
traj iread 30
traj read
close unit 30
open unit 30 unform read name @pdbfile_e@begin-@end_nowat.crd
open unit 40 form  write name @pdbfile_e@begin-@end_nowat.trms_dumb
open unit 41 form  write name @pdbfile_e@begin-@end_nowat.trms_card"
coor copy comp

correl maxt 1000 maxseries 100 maxatoms 1500 bycb -
     eps 1.0 cutnb @cutoff cutim @cutoff ctofnb @ctofnb ctonnb @ctonnb vswi -
     Ewald kappa 0.320 pmEwald order 4 fftx 64 ffty 64 fftz 64 
  enter rms rms orient
  traj first 30 nunit 1
  show all
  write rms unit 40 dumb
  write rms unit 41 card
end

system "grep -v "^*" 1fsc_e1-8_nowat.trms_card | tail +12 > 1fsc_e1-8_nowat.trms"

stop
label trmswat
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      Not set up for this yet
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
stop

label p2
!*********************************************************
!             P2 correlation of NH vectors
!*********************************************************
if @water eq true goto trmswat
delete atom sele resn TIP3 end

set first 1
set nres ?nres
set nres 61
set res 1
!label p2opnfil
!  calc fil = 20 + @res
!  open unit @fil form  write name @pdbfile_e@begin-@end_p2_@fil.p2
!  incr res
!  if @res le @nres goto p2opnfil

open unit 20 unform read name @pdbfile_e@begin-@end_nowat.crd
set res 1
label p2enter
correl maxt 1000 maxseries 20 maxatoms 1500 bycb -
     eps 1.0 cutnb @cutoff cutim @cutoff ctofnb @ctofnb ctonnb @ctonnb vswi -
     Ewald kappa 0.320 pmEwald order 4 fftx 64 ffty 64 fftz 64 

   define nvec sele resi @res .and. type N end
   if ?nsel ne 1 goto skp2 
   define hvec sele resi @res .and. type HN end
   if ?nsel ne 1 goto skp2 
   enter nv atom xyz sele nvec end
   enter hv atom xyz sele hvec end
   traj first 20 nunit 1
   mantime hv mult -1.
   mantime nv add hv
   mantime nv normal
   corfun nv nv p2
   open unit 30 form  write name @pdbfile_e@begin-@end_p2_@res.p2
   write corr unit 30 dumb
   if @first eq 0 goto skp2
      open unit 31 form  write name @pdbfile_e@begin-@end_p2_NUM.p2
      write corr unit 31 card
      close unit 31
      set first 0
label skp2
end

   incr res
   if @res le @nres goto p2enter
stop

label stripwat
!*********************************************************
!               STRIP WATER
!*********************************************************
open unit 30 unform read name @pdbfile_e@begin-@end.crd
  traj iread 30
  traj read
  close unit 30
set start ?start
set skip  ?skip
set nfile ?nfile
calc stop =  @skip * @nfile  + @start - @skip

open unit 30 unform read name @pdbfile_e@begin-@end.crd
open unit 20 unform write name @pdbfile_e@begin-@end_NoWat.crd

define notwat sele .not. resn TIP3 end
merge firstu 30 nunit 1 skip @skip begin @start stop @stop nfile @nfile -
      output 20 sele notwat end
close unit 20
close unit 30

open unit 30 form write name  @pdbfile_NoWat.psf
open unit 31 form write name @pdbfile_NoWat.crd
open unit 32 form write name @pdbfile_NoWat.pdb

delete atom sele .not. notwat end
write psf card unit 30
write coor card unit 31
write coor pdb unit 32
close unit 31
close unit 32
close unit 33

stop

label format
!*********************************************************
!                FORMAT
!*********************************************************
!
calc nunit = @end - @begin +1
  open unit 30 unform read name @pdbfile_e@begin.crd
  traj iread 30
  traj read
  close unit 30
set start ?start
set skip  ?skip
set nframes ?nfile
calc nfile = ?nfile * @nunit
calc stop =  @skip * @nframes * @nunit + @start - @skip

set unit 30
set current @begin
label dynopen
  open unit @unit unform read name @pdbfile_e@current.crd
  incr unit
  incr current
if @current le @end goto dynopen 
open unit 20 write unform name @pdbfile_e@begin-@end.crd
merge firstu 30 nunit @nunit skip @skip output 20 nfile @nfile
close unit 20
close unit 30

open unit 30 read unform name @pdbfile_e@begin-@end.crd
open unit 20 write form name @pdbfile_e@begin-@end.fcrd

dyna format first 30 nunit 1 output 20 -
            skip @skip begin @start stop @stop

stop
!
label unformat
!*********************************************************
!                  UNFORMAT
!*********************************************************

set unit 30
  open unit @unit unform write name @pdbfile_e@begin-@end.crd
!  open unit @unit unform write name a12

open unit 20 read form name @pdbfile_e@begin-@end.fcrd
dyna unformat input 20  output 30
stop
label unmerge
!*********************************************************
!                unmerge
!*********************************************************
!

calc nunit = @end - @begin +1
  open unit 30 unform read name @pdbfile_e@begin-@end.crd
  traj iread 30
  traj read
  close unit 30
set start ?start
set skip  ?skip
calc nframes ?nfile 
calc stop =  @skip * @nframes  + @start - @skip
calc nfile = ?nfile / @nunit

set unit 30
set current @begin
label mergopen
  open unit @unit unform write name @pdbfile_e@current.crd
  incr unit
  incr current
if @current le @end goto mergopen 
open unit 20 unform read name @pdbfile_e@begin-@end.crd

merge firstu 20 nunit 1 output 30 nfile @nfile -
            skip @skip begin @start stop @stop

stop
!===========================================================
!*********************************************************
!*******    END OF SCRIPT  subroutines follow     ********
!*********************************************************
!===========================================================

label iamdead
system "banner Must specify coordinate/psf filename"

stop


label deadformat
system "banner 'Must '"
system "banner 'specify '"
system "banner 'task='"
system "banner '  stripwat'"
system "banner '  format'"
system "banner '  unformat'"

stop

label getparams
!*********************************************************
!     get parameters set for this run
!*********************************************************
!

if @?toppar eq 0 set toppar = "$toppar"
if @?T eq 0 set T = 298
if @?begin eq 0 set begin = 1
if @?end eq 0 set end = 10
if @?current eq 0 set current = 1
if @?shape eq 0 set shape = cubic
if @?crystal eq 0 set crystal = false
if @?cutoff eq 0 set cutoff = 9
if @?status eq 0 set status start
if @?water eq 0 set water false

! Dimension of the solvent volume
Calc 7 = @size * 1.005 ! Initially set boxsize slightly larger for CPT equilibration
  if @shape ne octahedral goto sizeok
   if @crystal eq true calc 7 = @7 * sqrt(3) / 2.
  label sizeok
set 8 @7
set 9 @8

!  Set-up cutoffs for non-bonded list and energy calculations
Calc cutoff = @cutoff + 1  ! Add additional 1A to cutoff from build
set ctofnb = @cutoff
Calc ctofnb = @ctofnb - 1.5
set ctonnb = @ctofnb
Calc ctonnb = @ctonnb - 1.0

goto doneparams


label extractrms
!*********************************************************
!     fix rms file for plotting
!*********************************************************
!

system "/bin/rm -f rms_fnam rms.csh awk.file"
system "echo 'BEGIN {header = 1 }' > awk.file"
system "echo '{if ( header == 1 ){' >> awk.file"
system "echo 'if ($1 == "\*") header = 1;' >> awk.file"
system "echo '       else     header = 0; } ' >> awk.file"
system "echo 'else  print $1, $10; }' >> awk.file"

system "echo @AVGNAM | awk '{print tolower($1)}' >rms_fnam `"
system "echo ' set fnam=`cat rms_fnam`; set fnam=$fnam:r ; /bin/rm -f $fnam.rms ' >> rms.csh"
system "echo 'awk -f awk.file $fnam.avgcrd > $fnam.rms ' >> rms.csh"
system "csh rms.csh"
system "/bin/rm -f rms_fnam rms.csh awk.file"

goto doneextrms