* Calculate free energy
* for each residu.
*
bomblev -2 

$SET$

! Topology, parameters and PSF
$LIB$

!Initial psf 
$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

! Results files
open unit 11 write formatted name @path/gbmv/gaz-decompo.dat ! for gas energies
write title unit 11
*# Interaction calculation
*# Resnum  Segid  Resname  DGgas   gas/2
*
close unit 11

open unit 11 write formatted name @path/gbmv/solv-decompo.dat ! for solv energies
write title unit 11
*# Solv calculation
*# Part GBEN
*
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

! Allow to know the number of residue of the first part
set res = 0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!inte between each resid of the first part and all the resid of the !
!second part                                                        !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!loop begins at the first atom of the fp
set j = 0
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

NBONDS  cutnb 999 ctonnb 979 ctofnb 989 -
vswitch shift cdie eps 1 !-

SKIPe  ALL EXCLude ELECtrostatic VDW
inte sele (resid @residue .and. segid @segid) end sele @sp end
! this is also the total int. energy, for consistency check
set a = ?ENER
set b = ?VDW
set c = ?ELEC
calc d = @c/2

! let's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/gbmv/gaz-decompo.dat
write title unit 12
*  ?SELRESI  ?SELSEGI ?SELRESN  ?ELEC @d
*
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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!inte between each resid of the second part and all the resid of the !
!first part                                                          !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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

NBONDS  cutnb 999 ctonnb 979 ctofnb 989 -
vswitch shift cdie eps 1 !-

SKIPe  ALL EXCLude ELECtrostatic VDW
inte sele  @fp end sele (segid @segid .and. resid @residue) end

! this is also the total int. energy, for consistency check
set a = ?ENER
set b = ?VDW
set c = ?ELEC
calc d = @c/2

! let's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/gbmv/gaz-decompo.dat
write title unit 12
*  ?SELRESI  ?SELSEGI ?SELRESN  ?ELEC @d
*
close unit 12
FORMAT
!Now go on with next residue
incr k by @resiatom
if k le  @lastatomsp goto  loop2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!This part allows to calculate the DG gas of each aminoacid         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
NBONDS  cutnb 999 ctonnb 979 ctofnb 989 -
vswitch shift cdie eps 1 !-

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's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/gbmv/gaz-decompo.dat
write title unit 12
*  Total = ?ELEC
*
close unit 12
FORMAT


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Calculation of the born radii of each aminoacid of the whole complex !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
scalar wmain = radii

GBMV beta -20 weight fixa epsilon 80 dn 1.0 watr 1.4 geom  -
     tol 1e-8 bufr 0.5 mem 10 cuta 20 hsx1 -0.125 hsx2 0.25 -
     alfrq 1 emp 1.5 P4 0.0 P6 8.0 P3 0.70 ONX 1.9 OFFX 2.1 -
     WTYP 2 NPHI 38 SHIFT -0.102 SLOPE 0.9085 CORR 1

Energy ctonnb 979 ctofnb 989 cutnb 999

scalar wmain = GBalpha

OPEN WRITE formatted unit 27 name @path/gbmv/@name_cx-alpha.crd
WRITE coor card unit 27
* All Hydrogen Coordinates Set @cx with alpha
*
CLOSE unit 27

! let's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/gbmv/solv-decompo.dat
write title unit 12
*  CX ?GBEN
*
close unit 12
FORMAT


Gbmv clear


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Calculation of the born radii of each aminoacid for the part A       !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
stream @path/inp/gbmv_fp.str

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Calculation of the born radii of each aminoacid for the part B       !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
stream @path/inp/gbmv_sp.str



stop