* Calculate interaction energy
* lysozyme and each residu of vhh
*

bomblev -2

! Code for this job
$SET$

! Topology, parameters and PSF
open unit 1 card read name lib/top_all27_prot_na.rtf
read RTF card unit 1
close unit 1

open unit 1 card read name lib/par_all27_prot_na.prm
read PARA card unit 1
close unit 1

stream lib/aco.rtf

!Initial psf 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!puts @path/scr/@inp.dyna_desolv.psf when it is the beginning psf!
!puts @path/scr/@out/@out.dyna_desolv.psf when it is a new psf   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
open unit 1 card read name @path/scr/@inp.dyna_desolv.psf 
read psf card unit 1
close unit 1

! Initial coordinates
open unit 1 card read name @path/scr/@out/@out-cha.crd.his
read coor card unit 1
close unit 1

define kma sele ( type ca2 .or. type ha21 .or. type ha22 ) end
define ct2 sele ( type nt .or. type ht1 .or. type ht2 ) end
define ct3 sele ( type nt .or. type hnt .or. type cat .or. type ht1 .or. type ht2 .or. type ht3 ) end
define ace sele ( type cay .or. type hy1 .or. type hy2 .or. type hy3 .or. type cy .or. type oy ) end

!Result file
open unit 11 write formatted name @path/scr/@out/@out.vdw.inte
write title unit 11
*# JOBNAME: @out
*# Interaction Energy Calculation
*# Residue   VdW.           Elec.          TOT.
*
close unit 11

!Allows to know the first and the last atom number of the fp
define first sele @fp end
set firstatomfp  ?SELATOM
set lastatomfp  ?SELATOM
incr lastatomfp by ?NSEL
decr lastatomfp by 1

!Allows to know the first and the last atom number of the sp
define second sele @sp end
set firstatomsp ?SELATOM
set lastatomsp ?SELATOM
incr lastatomsp by ?NSEL
decr lastatomsp by 1


NBONDS  cutnb 13.5 ctonnb 8.5 ctofnb 12.5 -
vswitch shift cdie eps 1 NBXM 0


SKIPe  ALL EXCLude ELECtrostatic VDW 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!puts the good keyword to select the 2 parts of the system!
!the order of the selection isn't important               !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
inte sele @fp end sele @sp end

! this is also the total int. energy, for consistency check
set a = ?ENER
set b = ?VDW
set c = ?ELEC
! let first put the total energy to zero
set tot 0

! Now let's loop over residues; this will calculate
! the interaction energy between the minimum and
! each protein residue 

set j 0

!loop begins at the first atom of the fp
incr j by @firstatomfp
label loop1

!allows to know the segid and the residue of the atom j and the number of atom in the residue
define atom sele bynum @j end
set residue ?SELRESI
set segid ?SELSEGI
define resi sele resid @residue .and. segid @segid end
set resiatom ?NSEL

inte sele @sp end sele ( resid @residue .and. segid @segid .and. (type N .or. type HN* .or. type CA .or. type HA .or. type C .or. type O .or. kma .or. ct2 .or. ace )) end
incr tot by ?ENER
! let's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/scr/@out/@out.vdw.inte
!if ?ENER NE 0.0001 then
write title unit 12
*BB  ?SELRESI  ?SELRESN    ?VDW    ?ELEC    ?ENER   @tot
*

inte sele @sp end sele ( resid @residue .and. segid @segid .and. .not. (type N .or. type HN* .or. type CA .or. type HA .or. type C .or. type O  .or. kma .or. ct2 .or. ace )) end
incr tot by ?ENER
! let's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/scr/@out/@out.vdw.inte
!if ?ENER NE 0.0001 then
write title unit 12
*SC  ?SELRESI  ?SELRESN    ?VDW    ?ELEC    ?ENER   @tot
*


close unit 12
FORMAT
! Now go on with next residue (first atom of the next residue)
incr j by @resiatom
if j le @lastatomfp goto loop1


! Compare the 2 values of total energy, just to check
open unit 12 append formatted name @path/scr/@out/@out.vdw.inte
write title unit 12
* sum of individual = @tot   TOTAL energy = @a
* Van der Waals   =  @b      Electrostatic = @c
*
close unit 12


set tot 0

set k = 0
incr k by @firstatomsp
label loop2

define atom sele bynum @k end
set residue ?SELRESI
set segid ?SELSEGI
define resi sele resid @residue .and. segid @segid end
set resiatom ?NSEL

inte sele (@fp) end sele ( resid @residue .and. segid @segid .and. (type N .or. type HN* .or. type CA .or. type HA .or. type C .or. type O .or. kma .or. ct2 .or. ace )) end
incr tot by ?ENER

! let's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/scr/@out/@out.vdw.inte
if ?ENER NE 0.0001 then
write title unit 12
*BB  ?SELRESI  ?SELRESN    ?VDW    ?ELEC    ?ENER   @tot
*

inte sele (@fp) end sele ( resid @residue .and. segid @segid .and. .not. (type N .or. type HN* .or. type CA .or. type HA .or. type C .or. type O .or. kma .or. ct2 .or. ace )) end
incr tot by ?ENER

! let's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/scr/@out/@out.vdw.inte
if ?ENER NE 0.0001 then
write title unit 12
*SC  ?SELRESI  ?SELRESN    ?VDW    ?ELEC    ?ENER   @tot
*



close unit 12
FORMAT
!Now go on with next residue

incr k by @resiatom
if k le  @lastatomsp goto  loop2

! Compare the 2 values of total energy, just to check
open unit 12 append formatted name @path/scr/@out/@out.vdw.inte
write title unit 12
* sum of individual = @tot   TOTAL energy = @a
* Van der Waals   =  @b      Electrostatic = @c
*
close unit 12

stop  !pfew