* This test case is aiming to enable the relatively automated
* evaluation  of test case output for the HFB method under TReK.
* By Ilja V. Khavrutskii 01/30/2007 (Adv: Charles L. Brooks III)
*

! Test to ensure the functionality exists in this build of CHARMM.
! Use exctly this format so that it can be parsed by test facilities.
if ?hfb .ne. 1 then
 echo HFB functionality absent: Test NOT performed
 stop
endif

! In datadir.def we define a variable pnode that should be
! used to "protect" all calls to prnlev
! additionally, the variable tol is defined to control
! differences in computed and reference variables.
set tol = 1.0e-4
stream datadir.def

! Format for prnlev shown here
prnlev 5 @pnode

label BEGINSETUP=0 ! This demarks the beginning of the setup region
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! In general use standard topology/parameter files, or include these
! directly in the input script to eliminate an abundance of files in
! the data directory. Test cases should be as self-contained as possible

! read parameter and topology files
open unit 10 read form name @0top_all22_prot.inp
read rtf unit 10 card
close unit 10

open unit 10 read form name @0par_all22_prot.inp
read param unit 10 card
close unit 10

set tolermsd 5.0e-5
set tolenerg 2.0e-2
set name ala
set nbpar atom cdie eps 1.0 cutnb 21 ctofnb 18 ctonnb 16 switch vswitch
set nbeads 32
set print 100
set sdstep 200
set abstep 1000

read sequ card
* Sequence of ala
*
1
ALA
generate ala first ACE last CT3 setup warn

IC PARA
IC SEED 1 CY 1 N 1 CA
IC BUILD

define rcs select type n .or. type c .or. type ca .or. type nt .or. type cy end

mini sd nstep @sdstep nprint @print @nbpar
mini abnr nstep @abstep nprint @print @nbpar

!This minimization gives the C7eq structure

coor orient mass
coor copy comp

open unit 10 write card name @9/ala_c7eq.chr
!open unit 10 write card name @9/ala_c7eq.pdb
write unit 10 coor card
!write unit 10 coor pdb
close unit 10

!do the rotation
coor axis sele resid 1 .and. type CA end sele resid 1 .and. type C end
coor rotate axis phi - 140.0 -
sele (type O .or. type NT .or. type HNT .or. type HT* .or. type CAT) end

coor axis sele resid 1 .and. type N end sele resid 1 .and. type CA end
coor rotate axis phi 150.0 -
sele (type CAY .or. type HY* .or. type OY .or. type CY .or. type HN) end

mini sd nstep @sdstep nprint @print @nbpar
mini abnr nstep @abstep nprint @print @nbpar

!This minimization gives the C7ax structure

coor orient mass rms sele rcs end

open unit 10 write card name @9/ala_c7ax.chr
!open unit 10 write card name @9/ala_c7ax.pdb
write unit 10 coor card
!write unit 10 coor pdb
close unit 10

! ready to start trek run #1
trek maxpoint @nbeads select rcs end mass
   traj read
     @9/ala_c7eq.chr !reactant
     @9/ala_c7ax.chr !product
   done
   fp trnc 24 init
   ! writing the path to a trajectory file
   traj write name @9/@name_init.dcd restart @9/@name.cpr 
quit

! ready to start trek run #2
trek maxpoint @nbeads select rcs end mass
   traj read
     @9/ala_c7eq.chr !reactant
     @9/ala_c7ax.chr !product
   done
   fp trnc 24 dens
   ! writing the path to a trajectory file
   traj write name @9/@name_dens.dcd restart @9/@name.cpr 
quit

!write out charmm coordinates from the generated trajectory
!with just reactant and product, dens gives the same result as init
open unit 11 read unform name @9/@name_init.dcd
traj iread 11 nunits 1 nfile @nbeads skip 1

set f 1
label tloop1
 traj read

 open unit 10 write card name @9/@name_init_r@{f}.chr
 write unit 10 coor card
 close unit 10

 incr f by 1
if @f .le. @nbeads then goto tloop1
close unit 11

!SHORTCUT FOR TEST DEVELOPMENT
!goto shortcut

!write out charmm coordinates from the generated trajectory
!with just reactant and product, dens gives the same result as init
open unit 11 read unform name @9/@name_dens.dcd
traj iread 11 nunits 1 nfile @nbeads skip 1

set meanrmsd 0.0
set f 1
label tloop2
 open unit 10 read card name @9/@name_init_r@{f}.chr
 read unit 10 coor card comp
 close unit 10

 traj read

 coor orient rms mass sele rcs end comp
 incr meanrmsd by ?rms

 open unit 10 write card name @9/@name_dens_r@{f}.chr
 write unit 10 coor card
 close unit 10

 incr f by 1
if @f .le. @nbeads then goto tloop2
close unit 11

calc meanrmsd = @meanrmsd/@nbeads

!Insert test here with some tolerance.
if @meanrmsd .gt. @tolermsd then
   echo TEST=1 FAILED
   stop
else
   echo TEST=1 PASSED
endif

!Removing a few random points from the path of 32 beads
!to demonstrate the robustness of the interpolation.
calc nfil = @nbeads - 4

open unit 11 write unform name @9/@name_trunc.dcd
traj iwrite 11 nwrite 1 nfile @nfil skip 1

set f 1
label wloop1

 if @f .eq. 5 then goto skip
 if @f .eq. 10 then goto skip
 if @f .eq. 17 then goto skip
 if @f .eq. 26 then goto skip

 open unit 10 read card name @9/@name_dens_r@{f}.chr
 read unit 10 coor card
 close unit 10

 traj write

label skip

 incr f by 1
if @f .le. @nbeads then goto wloop1
close unit 11

! ready to start trek run #3
trek maxpoint @nbeads select rcs end mass
   ! reading the truncated initial-guess path
   traj read name @9/@name_trunc.dcd begin 1 skip 1
   fp trnc 24 dens
   ! writing the path to a trajectory file
   traj write name @9/@name_full.dcd restart @9/@name.cpr 
quit

!write out charmm coordinates from the generated trajectory
!with just reactant and product, dens gives the same result as init
open unit 11 read unform name @9/@name_full.dcd
traj iread 11 nunits 1 nfile @nbeads skip 1

set meanrmsd 0.0
set f 1
label tloop3
 open unit 10 read card name @9/@name_init_r@{f}.chr
 read unit 10 coor card comp
 close unit 10

 traj read

 coor orient rms mass sele rcs end comp
 incr meanrmsd by ?rms

 incr f by 1
if @f .le. @nbeads then goto tloop3
close unit 11

calc meanrmsd = @meanrmsd/@nbeads

!Insert test here with some tolerance.
if @meanrmsd .gt. @tolermsd then
   echo TEST=2 FAILED
   stop
else
   echo TEST=2 PASSED
endif

!Now let's test the optimization with the primitive FP option
!We will use the evolution engine that is a slight modification 
!of the one written by Charles L. Brooks III.

set stp 1
set hfbstp 10

label fphfbloop

if stp .eq. 1 then
 set reference @name_init
 set evolution @name_mini
else
 set reference @name_refe
 set evolution @name_mini
endif

open unit 11 read unform name @9/@reference.dcd
open unit 12 write unform name @9/@evolution.dcd
traj iread 11 iwrite 12 nunit 1 nfile @nbeads skip 1

set f 1
label loopbeads1
    traj read
    coor copy comp
!    cons harm force 10 mass select rcs end
    cons harm force 10 mass select rcs end comp
    mini sd nstep @sdstep nprint @print @nbpar
    mini abnr nstep @abstep nprint @print @nbpar tolgrd 1.0e-5
    cons harm clear
!    cons harm force 0 select all end
    traj write
    incr f by 1
if f .le. @nbeads goto loopbeads1

set reference @name_refe

! ready to start trek run #4
trek maxpoint @nbeads select rcs end mass
   ! reading the truncated initial-guess path
   traj read name @9/@evolution.dcd begin 1 skip 1
   fp trnc 24
   ! writing the path to a trajectory file
   traj write name @9/@reference.dcd restart @9/@name.cpr
quit

 incr stp by 1
if stp .le. @hfbstp then goto fphfbloop

! Testing the progress in terms of mean rmsd
open unit 11 read unform name @9/@reference.dcd
traj iread 11 nunits 1 nfile @nbeads skip 1

set meanrmsd 0.0
set f 1
label tloop4
 open unit 10 read card name @9/@name_init_r@{f}.chr
 read unit 10 coor card comp
 close unit 10

 traj read

 coor orient rms mass sele rcs end comp
 incr meanrmsd by ?rms

 incr f by 1
if @f .le. @nbeads then goto tloop4
close unit 11

calc meanrmsd = @meanrmsd/@nbeads

!Insert test here with some tolerance.
calc diff = abs(@meanrmsd - 7.429156E-02)
if @diff .gt. @tolermsd then
   echo TEST=3 FAILED
   stop
else
   echo TEST=3 PASSED
endif

!Now let's test the optimization with the gradient enhanced WRKP option
!We will use the evolution engine that is a slight modification 
!of the one written by Charles L. Brooks III.

set stp 1
set hfbstp 10
set evolution @name_mini
set newref @name_refe
calc maxbeads = @nbeads * 2

label wrkphfbloop

if stp .eq. 1 then
 set oldref @name_init
endif

open unit 11 read unform name @9/@oldref.dcd
open unit 12 write unform name @9/@evolution.dcd
traj iread 11 iwrite 12 nunit 1 nfile @nbeads skip 1

set f 1
label loopbeads2
    traj read
    coor copy comp
!    cons harm force 10 mass select rcs end
    cons harm force 10 mass select rcs end comp
    mini sd nstep @sdstep nprint @print @nbpar
    mini abnr nstep @abstep nprint @print @nbpar tolgrd 1.0e-5
    cons harm clear
!    cons harm force 0 select all end
    traj write
    incr f by 1
if f .le. @nbeads goto loopbeads2

! ready to start trek run #5
trek maxpoint @maxbeads select rcs end mass
   ! reading the truncated initial-guess path
   traj read name @9/@evolution.dcd begin 1 skip 1
   traj read name @9/@oldref.dcd begin 1 skip 1
   wrkp trnc 24 bead @nbeads grid @maxbeads rtyp 1 rfrc 10.0 step 0.0
   ! writing the path to a trajectory file
   traj write name @9/@newref.dcd restart @9/@name.cpr
quit

set oldref @newref

 incr stp by 1
if stp .le. @hfbstp then goto wrkphfbloop

! Testing the progress in terms of mean rmsd
open unit 11 read unform name @9/@newref.dcd
traj iread 11 nunits 1 nfile @nbeads skip 1

set meanrmsd 0.0
set f 1
label tloop5
 open unit 10 read card name @9/@name_init_r@{f}.chr
 read unit 10 coor card comp
 close unit 10

 traj read

 coor orient rms mass sele rcs end comp
 incr meanrmsd by ?rms

 incr f by 1
if @f .le. @nbeads then goto tloop5
close unit 11

calc meanrmsd = @meanrmsd/@nbeads

!Insert test here with some tolerance.
calc diff = abs(@meanrmsd - 7.429156E-02)
if @diff .gt. @tolermsd then
   echo TEST=4 FAILED
   stop
else
   echo TEST=4 PASSED
endif

!label shortcut

!Now let's test the optimization with the gradient enhanced WRKP 
!and the SD option (non-zero step size parameter)
!We will use the evolution engine that is a slight modification 
!of the one written by Charles L. Brooks III.

set stp 1
set hfbstp 10
calc maxbeads = @nbeads * 2

label wrkphfbloop2

calc nxt = @stp + 1
set evolution @name_mini_s@stp
set newref @name_refe_s@nxt
if stp .eq. 1 then
 set oldref @name_init
else
 set oldref @name_refe_s@stp
endif

open unit 11 read unform name @9/@oldref.dcd
open unit 12 write unform name @9/@evolution.dcd
traj iread 11 iwrite 12 nunit 1 nfile @nbeads skip 1

set f 1
label loopbeads3
    traj read
    coor copy comp
!    cons harm force 10 mass select rcs end
    cons harm force 10 mass select rcs end comp
    mini sd nstep @sdstep nprint @print @nbpar
    mini abnr nstep @abstep nprint @print @nbpar tolgrd 1.0e-5
    cons harm clear
!    cons harm force 0 select all end
    traj write
    incr f by 1
if f .le. @nbeads goto loopbeads3

! ready to start trek run #6
trek maxpoint @maxbeads select rcs end mass
   ! reading the truncated initial-guess path
   traj read name @9/@evolution.dcd begin 1 skip 1
   traj read name @9/@oldref.dcd begin 1 skip 1
!   wrkp trnc 24 rtyp 1 rfrc 10.0 step 0.0000125
   wrkp trnc 24 rtyp 1 rfrc 10.0 step 0.004
   ! writing the path to a trajectory file
   traj write name @9/@newref.dcd restart @9/@name.cpr
quit

 incr stp by 1
if stp .le. @hfbstp then goto wrkphfbloop2

! Testing the progress in terms of mean rmsd
open unit 11 read unform name @9/@newref.dcd
traj iread 11 nunits 1 nfile @nbeads skip 1

set meanrmsd 0.0
set f 1
label tloop6
 open unit 10 read card name @9/@name_init_r@{f}.chr
 read unit 10 coor card comp
 close unit 10

 traj read

 coor orient rms mass sele rcs end comp
 incr meanrmsd by ?rms

 incr f by 1
if @f .le. @nbeads then goto tloop6
close unit 11

calc meanrmsd = @meanrmsd/@nbeads

!Insert test here with some tolerance.
!calc diff = abs(@meanrmsd - 9.360906E-02)
calc diff = abs(@meanrmsd - 9.8475E-02)
if @diff .gt. @tolermsd then
   echo TEST=5 FAILED
   stop
else
   echo TEST=5 PASSED
endif

!Finally let's test the energy profile. We'll use correl for comparison
!of the two profiles (from energy optimization with rcs fixed and 
!from the reversible work line integral)

open write card unit 14 name @9/reversible_work.dat

calc maxbeads = @nbeads * 8 !just for test, 4 otherwise
calc ngrid = @nbeads * 4
! ready to start trek run #6
trek maxpoint @maxbeads select rcs end mass
   ! reading the truncated initial-guess path
   traj read name @9/@evolution.dcd begin 1 skip 1
   traj read name @9/@oldref.dcd begin 1 skip 1
!   wrkp trnc 24 bead @nbeads grid @ngrid rtyp 1 rfrc 10.0 step 0.0000125 iprf 14 ptrj
   wrkp trnc 24 bead @nbeads grid @ngrid rtyp 1 rfrc 10.0 step 0.004 iprf 14 ptrj
   ! writing the path to a trajectory file
   traj write name @9/@name_profile.dcd restart @9/@name.cpr
quit

!write out charmm coordinates from the generated trajectory
!with just reactant and product, dens gives the same result as init
open write card unit 15 name @9/energy_profile.dat
open unit 11 read unform name @9/@name_profile.dcd
traj iread 11 nunits 1 nfile @ngrid skip 1

set f 1
label tloop7
 traj read
 cons fix sele rcs end
 mini sd nstep @sdstep nprint @print @nbpar
 mini abnr nstep @abstep nprint @print @nbpar tolgrd 1.0e-5
 cons fix sele none end
 energy @nbpar
 if f .eq. 1 then
  set eref ?ener
 endif
 calc erel = ?ener - @eref

write title unit 15
* @f @erel
*
 incr f by 1
if @f .le. @ngrid then goto tloop7
close unit 11

open unit 14 read formatted name @9/reversible_work.dat
open unit 15 read formatted name @9/energy_profile.dat
set  mxts 10
set  mxat @maxbeads
set  mxtt @maxbeads
correl maxtime @mxtt maxseries @mxts maxatoms @mxat noupdate
     enter wrk ener
     enter prf ener
     read wrk unit 14 dumb column 2
     read prf unit 15 dumb column 2
     mantime wrk mult -1
     mantime prf add wrk
     mantime prf square
set average ?aver
     show all
end

!Insert a simple test here with some tolerance.
calc ermsd = sqrt(@average)
if @ermsd .gt. @tolenerg then
   echo TEST=6 FAILED
   stop
else
   echo TEST=6 PASSED
endif

!This is it.
stop