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

bomblev -2

! Code for this job

$SET$

! Topology, parameters and PSF
$LIB$

!Initial psf 
open unit 1 card read name @path/scr/@name_complex.psf
read psf card unit 1
close unit 1


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


!Result file
open unit 11 write formatted name @path/scr/@name.vdw.inte
write title unit 11
*# JOBNAME: @name
*# 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 
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 ) end 
incr tot by ?ENER

! let's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/scr/@name.vdw.inte
!if ?ENER NE 0.0001 then 
write title unit 12
*  ?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/@name.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 ) end 

incr tot by ?ENER

! let's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/scr/@name.vdw.inte
if ?ENER NE 0.0001 then 
write title unit 12
*  ?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/@name.vdw.inte
write title unit 12
* sum of individual = @tot   TOTAL energy = @a
* Van der Waals   =  @b      Electrostatic = @c
*
close unit 12


stop  !pfew