! #include "copyright.h" #include "../include/dprec.fh" #define REQUIRE(e) if(.not.(e)) call croak(__FILE__,__LINE__) #define ASSERT(e) if(.not.(e)) call croak(__FILE__,__LINE__) module decomp !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ! --- Module for energy decomposition ! ! Holger Gohlke ! 5.1.2004 ! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ integer, public, dimension(:), allocatable :: jgroup, index, irespw integer, private, parameter :: ndectype = 28 ! number of "enum" variables following integer, private, parameter :: & ! DECOMP contributions for sidechains sideintsel = 1, & ! DECOMP -ene int sel sideintind = 2, & ! DECOMP ene int ind sidevdwsel = 3, & ! DECOMP -ene vdw sel sidevdwind = 4, & ! DECOMP ene vdw ind sidevdwdir = 5, & ! DECOMP ene vdw dir sideeelsel = 6, & ! DECOMP -ene eel sel sideeelind = 7, & ! DECOMP ene eel ind sideeeldir = 8, & ! DECOMP ene eel dir sidepolsel = 9, & ! DECOMP -ene pol sel sidepolind = 10, & ! DECOMP ene pol ind sidepoldir = 11, & ! DECOMP ene pol dir sidesassel = 12, & ! DECOMP -ene sas sel sidesasind = 13, & ! DECOMP ene sas ind sidesasdir = 14, & ! DECOMP ene sas dir ! DECOMP contributions for backbone backintsel = 15, & ! DECOMP -ene int sel backintind = 16, & ! DECOMP ene int ind backvdwsel = 17, & ! DECOMP -ene vdw sel backvdwind = 18, & ! DECOMP ene vdw ind backvdwdir = 19, & ! DECOMP ene vdw dir backeelsel = 20, & ! DECOMP -ene eel sel backeelind = 21, & ! DECOMP ene eel ind backeeldir = 22, & ! DECOMP ene eel dir backpolsel = 23, & ! DECOMP -ene pol sel backpolind = 24, & ! DECOMP ene pol ind backpoldir = 25, & ! DECOMP ene pol dir backsassel = 26, & ! DECOMP -ene sas sel backsasind = 27, & ! DECOMP ene sas ind backsasdir = 28 ! DECOMP ene sas dir integer, private, dimension(ndectype) :: ndecind ! index into the dec array integer, private :: ndecno ! length of the dec array minus 1 _REAL_, private, dimension(:), allocatable :: dec ! stores decomposition values contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ allocates scratch space for integer arrays of module decomp subroutine allocate_int_decomp( natom, nres ) implicit none integer, intent(in) :: natom integer, intent(in) :: nres integer ier allocate( jgroup(natom), stat = ier ) REQUIRE( ier == 0 ) jgroup(:) = 0 allocate( index(nres), stat = ier ) REQUIRE( ier == 0 ) index(:) = 0 allocate( irespw(nres), stat = ier ) REQUIRE( ier == 0 ) irespw(:) = 0 return end subroutine allocate_int_decomp !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ allocates scratch space for real array of module decomp subroutine allocate_real_decomp( nin ) implicit none integer, intent(in) :: nin integer i, j, ncnt, n, ier ndecno = nin n = nin + 1 ! +1 for "rest" storage allocate( dec(ndectype * n), stat = ier ) REQUIRE( ier == 0 ) ! Init ndecind and dec arrays ncnt = 1 do i=1,ndectype ndecind(i) = ncnt do j=1,n dec(ncnt) = 0.0d0 ncnt = ncnt + 1 end do end do return end subroutine allocate_real_decomp !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ deallocates scratch space for integer arrays of module decomp subroutine deallocate_int_decomp( ) implicit none integer ier if ( allocated( jgroup ) ) then deallocate( jgroup, stat = ier ) REQUIRE( ier == 0 ) else ASSERT( .false. ) ! cannot deallocate un-allocated array end if if ( allocated( index ) ) then deallocate( index, stat = ier ) REQUIRE( ier == 0 ) else ASSERT( .false. ) ! cannot deallocate un-allocated array end if if ( allocated( irespw ) ) then deallocate( irespw, stat = ier ) REQUIRE( ier == 0 ) else ASSERT( .false. ) ! cannot deallocate un-allocated array end if return end subroutine deallocate_int_decomp !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ deallocates scratch space for real array of module decomp subroutine deallocate_real_decomp( ) implicit none integer ier if ( allocated( dec ) ) then deallocate( dec, stat = ier ) REQUIRE( ier == 0 ) else ASSERT( .false. ) ! cannot deallocate un-allocated array end if return end subroutine deallocate_real_decomp !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine decpair here] subroutine decpair(nty,nat1,nat2,fval) ! Accumulates contributions fval into ENE_XXX_SEL/IND/DIR, ! where XXX stands for EGB, VDW, EEL or ! INT in the case of bond energies ! Holger Gohlke ! 12.11.2001 use memory_module, only : nres, npdec implicit none _REAL_ fval,hfval integer nty, ntype integer nat1, nat2 integer nres1, nres2 integer nind1, nind2 integer nssel, nsind, nsdir integer nbsel, nbind, nbdir integer ntmp logical isside1, isside2 logical isprot1, isprot2 logical pairwise ! mjhsieh: warnings eliminator nssel = -1; nsdir = -1; nsind = -1 nbsel = -1; nbind = -1; nbdir = -1 hfval = 0.5d0 * fval if(nty < 0) then pairwise = .true. ntype = -nty else pairwise = .false. ntype = nty end if if (ntype == 1) then ! --- GB decomposition nssel = ndecind(sidepolsel) nsind = ndecind(sidepolind) nsdir = ndecind(sidepoldir) nbsel = ndecind(backpolsel) nbind = ndecind(backpolind) nbdir = ndecind(backpoldir) else if (ntype == 2) then ! --- EEL decomposition ! (+ 1-4 EEL if idecomp == 2 or 4; ! this follows decomp. for mm_pbsa) nssel = ndecind(sideeelsel) nsind = ndecind(sideeelind) nsdir = ndecind(sideeeldir) nbsel = ndecind(backeelsel) nbind = ndecind(backeelind) nbdir = ndecind(backeeldir) else if (ntype == 3) then ! --- VDW decomposition ! (+ 1-4 VDW if idecomp == 2 or 4; ! this follows decomp for mm_pbsa) nssel = ndecind(sidevdwsel) nsind = ndecind(sidevdwind) nsdir = ndecind(sidevdwdir) nbsel = ndecind(backvdwsel) nbind = ndecind(backvdwind) nbdir = ndecind(backvdwdir) else if (ntype == 4) then ! --- BOND, UREY-BRADLEY decomposition ! (+ 1-4 EEL + 1-4 VDW if idecomp == 1 or 3; ! this follows decomp as done within sander) nssel = ndecind(sideintsel) nsind = ndecind(sideintind) nsdir = nsind ! "Hack" that assures that cases such nbsel = ndecind(backintsel) ! as nres1 = 0, nres2 != 0 will work nbind = ndecind(backintind) nbdir = nbind ! dito else write(6,'(a,i6)') 'Wrong input for ntype: ',ntype call mexit(6,1) end if ! (ntype == 1) ! --- Determine if side or back and if prot or lig isside1 = .true. isprot1 = .true. nres1 = jgroup(nat1) if(nres1 < 0) then isprot1 = .false. nres1 = -nres1 end if if(nres1 > nres) then isside1 = .false. nres1 = nres1 - nres end if isside2 = .true. isprot2 = .true. nres2 = jgroup(nat2) if(nres2 < 0) then isprot2 = .false. nres2 = -nres2 end if if(nres2 > nres) then isside2 = .false. nres2 = nres2 - nres end if ! --- Echo result ! write(6,*) nat1,nres1,isside1,isprot1, & ! nat2,nres2,isside2,isprot2 ! --- Decompose if (nres1 == nres2) then ! --- self-energy of residue if (isside1) then nind1 = nssel else nind1 = nbsel end if if (isside2) then nind2 = nssel else nind2 = nbsel end if else if (( isprot1 .and. isprot2) .or. & (.not.isprot1 .and. .not.isprot2)) then ! --- indirect energy from within own molecule if (isside1) then nind1 = nsind else nind1 = nbind end if if (isside2) then nind2 = nsind else nind2 = nbind end if else ! --- direct energy between molecules if (isside1) then nind1 = nsdir else nind1 = nbdir end if if (isside2) then nind2 = nsdir else nind2 = nbdir end if end if end if ! (nres1 == nres2) ! --- Accumulate results if(nind1 == 0 .or. nind2 == 0) then write(6,'(a)') 'NIND1 or NIND2 wrong in DECPAIR.' write(6,'(a,i6)') 'NTYPE = ',ntype write(6,'(8i6)') nat1,nres1,isside1,isprot1, & nat2,nres2,isside2,isprot2 call mexit(6,1) end if if(nres1 > 0 .and. nres2 > 0) then ! --- Interactions "wanted" to be considered if(pairwise) then ! --- Pairwise decomp ntmp = (index(nres1) - 1) * npdec + index(nres2) dec(nind1 + ntmp - 1) = dec(nind1 + ntmp - 1) + hfval ntmp = (index(nres2) - 1) * npdec + index(nres1) dec(nind2 + ntmp - 1) = dec(nind2 + ntmp - 1) + hfval else ! --- Residue-based decomp dec(nind1 + nres1 - 1) = dec(nind1 + nres1 - 1) + hfval dec(nind2 + nres2 - 1) = dec(nind2 + nres2 - 1) + hfval end if else ! --- Collect as "rest" dec(nind1 + ndecno) = dec(nind1 + ndecno) + hfval dec(nind2 + ndecno) = dec(nind2 + ndecno) + hfval end if return end subroutine decpair !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine decangle here] subroutine decangle(nat1,nat2,nat3,fval) ! Accumulates contributions fval into ENE_INT_SEL/IND, ! for ANGLE energies ! Holger Gohlke ! 12.11.2001 use memory_module, only : nres implicit none _REAL_ fval,hfval integer nat1, nat2, nat3 integer nres1, nres2, nres3 integer nind1, nind2, nind3 logical isside1, isside2, isside3 integer nssel, nsind, nsdir integer nbsel, nbind, nbdir hfval = 0.3333333333333333334d0 * fval nssel = ndecind(sideintsel) nsind = ndecind(sideintind) nsdir = 0 nbsel = ndecind(backintsel) nbind = ndecind(backintind) nbdir = 0 ! --- Determine if side or back and if prot or lig isside1 = .true. nres1 = iabs(jgroup(nat1)) if(nres1 > nres) then isside1 = .false. nres1 = nres1 - nres end if isside2 = .true. nres2 = iabs(jgroup(nat2)) if(nres2 > nres) then isside2 = .false. nres2 = nres2 - nres end if isside3 = .true. nres3 = iabs(jgroup(nat3)) if(nres3 > nres) then isside3 = .false. nres3 = nres3 - nres end if ! --- Echo result ! write(6,*) nat1,nres1,isside1, ! + nat2,nres2,isside2, ! + nat3,nres3,isside3 ! --- Decompose if ((nres1 == nres2) .and. (nres1 == nres3)) then ! --- self-energy of residue if (isside1) then nind1 = nssel else nind1 = nbsel end if if (isside2) then nind2 = nssel else nind2 = nbsel end if if (isside3) then nind3 = nssel else nind3 = nbsel end if else ! --- indirect energy from within own molecule if (isside1) then nind1 = nsind else nind1 = nbind end if if (isside2) then nind2 = nsind else nind2 = nbind end if if (isside3) then nind3 = nsind else nind3 = nbind end if end if ! ((nres1 == nres2) .and. (nres1 == nres3)) ! --- Accumulate results if(nind1 == 0 .or. nind2 == 0 .or. nind3 == 0) then write(6,'(a)') 'NIND1 or NIND2 or NIND3 wrong in DECPAIR.' write(6,'(9i6)') nat1,nres1,isside1, & nat2,nres2,isside2, & nat3,nres3,isside3 call mexit(6,1) end if if(nres1 > 0 .and. nres2 > 0 .and. nres3 > 0) then ! --- Interactions "wanted" to be considered dec(nind1 + nres1 - 1) = dec(nind1 + nres1 - 1) + hfval dec(nind2 + nres2 - 1) = dec(nind2 + nres2 - 1) + hfval dec(nind3 + nres3 - 1) = dec(nind3 + nres3 - 1) + hfval else ! --- Collect as "rest" dec(nind1 + ndecno) = dec(nind1 + ndecno) + hfval dec(nind2 + ndecno) = dec(nind2 + ndecno) + hfval dec(nind3 + ndecno) = dec(nind3 + ndecno) + hfval end if return end subroutine decangle !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine decphi here] subroutine decphi(nat1,nat2,nat3,nat4,fval) ! Accumulates contributions fval into ENE_INT_SEL/IND, ! for PHI energies ! Holger Gohlke ! 12.11.2001 use memory_module, only : nres implicit none _REAL_ fval,hfval integer nat1, nat2, nat3, nat4 integer nres1, nres2, nres3, nres4 integer nind1, nind2, nind3, nind4 logical isside1, isside2, isside3, isside4 integer nssel, nsind, nsdir integer nbsel, nbind, nbdir hfval = 0.25d0 * fval nssel = ndecind(sideintsel) nsind = ndecind(sideintind) nsdir = 0 nbsel = ndecind(backintsel) nbind = ndecind(backintind) nbdir = 0 ! --- Determine if side or back and if prot or lig isside1 = .true. nres1 = iabs(jgroup(nat1)) if(nres1 > nres) then isside1 = .false. nres1 = nres1 - nres end if isside2 = .true. nres2 = iabs(jgroup(nat2)) if(nres2 > nres) then isside2 = .false. nres2 = nres2 - nres end if isside3 = .true. nres3 = iabs(jgroup(nat3)) if(nres3 > nres) then isside3 = .false. nres3 = nres3 - nres end if isside4 = .true. nres4 = iabs(jgroup(nat4)) if(nres4 > nres) then isside4 = .false. nres4 = nres4 - nres end if ! --- Echo result ! write(6,*) nat1,nres1,isside1, ! + nat2,nres2,isside2, ! + nat3,nres3,isside3, ! + nat4,nres4,isside4 ! --- Decompose if ((nres1 == nres2) .and. & (nres1 == nres3) .and. & (nres1 == nres4)) then ! --- self-energy of residue if (isside1) then nind1 = nssel else nind1 = nbsel end if if (isside2) then nind2 = nssel else nind2 = nbsel end if if (isside3) then nind3 = nssel else nind3 = nbsel end if if (isside4) then nind4 = nssel else nind4 = nbsel end if else ! --- indirect energy from within own molecule if (isside1) then nind1 = nsind else nind1 = nbind end if if (isside2) then nind2 = nsind else nind2 = nbind end if if (isside3) then nind3 = nsind else nind3 = nbind end if if (isside4) then nind4 = nsind else nind4 = nbind end if end if ! ((nres1 == nres2) .and. ! --- Accumulate results if(nind1 == 0 .or. nind2 == 0 .or. & nind3 == 0 .or. nind4 == 0) then write(6,'(a)') & 'NIND1 or NIND2 or NIND3 or NIND4 wrong in DECPAIR.' write(6,'(12i6)') nat1,nres1,isside1, & nat2,nres2,isside2, & nat3,nres3,isside3, & nat4,nres4,isside4 call mexit(6,1) end if if(nres1 > 0 .and. nres2 > 0 .and. & nres3 > 0 .and. nres4 > 0) then ! --- Interactions "wanted" to be considered dec(nind1 + nres1 - 1) = dec(nind1 + nres1 - 1) + hfval dec(nind2 + nres2 - 1) = dec(nind2 + nres2 - 1) + hfval dec(nind3 + nres3 - 1) = dec(nind3 + nres3 - 1) + hfval dec(nind4 + nres4 - 1) = dec(nind4 + nres4 - 1) + hfval else ! --- Collect as "rest" dec(nind1 + ndecno) = dec(nind1 + ndecno) + hfval dec(nind2 + ndecno) = dec(nind2 + ndecno) + hfval dec(nind3 + ndecno) = dec(nind3 + ndecno) + hfval dec(nind4 + ndecno) = dec(nind4 + ndecno) + hfval end if return end subroutine decphi !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine decsasa here] subroutine decsasa(nty,nat1,nat2,nat3,fval) ! Accumulates contributions fval into ENE_SAS_SEL/IND/DIR ! Holger Gohlke ! 12.11.2001 use memory_module, only : npdec, nres implicit none _REAL_ fval,hfval,addfval integer nty, ntype integer nat1, nat2, nat3 integer nres1, nres2, nres3 integer nind1, nind2, nind3 integer ntmp logical isside1, isside2, isside3 logical isprot1, isprot2, isprot3 logical pairwise integer nssel, nsind, nsdir integer nbsel, nbind, nbdir ! mjhsieh: warning eliminator ntype = -1; nres3 = -1; isprot3 = .false. addfval = -1; nind3 = -1; nind1 = -1 nres2 = -1; nind2 = -1; pairwise = .false. hfval = 0.5d0 * fval nssel = ndecind(sidesassel) nsind = ndecind(sidesasind) nsdir = ndecind(sidesasdir) nbsel = ndecind(backsassel) nbind = ndecind(backsasind) nbdir = ndecind(backsasdir) ! --- Check NTYPE if(nty == 0 .or. abs(nty) >= 4) then write(6,'(a,i6)') 'Wrong input for NTYPE: ',nty call mexit(6,1) else if(nty < 0) then pairwise = .true. ntype = -nty else pairwise = .false. ntype = nty end if ! --- Determine if side or back and if prot or lig isside1 = .true. isprot1 = .true. nres1 = jgroup(nat1) if(nres1 < 0) then isprot1 = .false. nres1 = -nres1 end if if(nres1 > nres) then isside1 = .false. nres1 = nres1 - nres end if if(ntype >= 2) then isside2 = .true. isprot2 = .true. nres2 = jgroup(nat2) if(nres2 < 0) then isprot2 = .false. nres2 = -nres2 end if if(nres2 > nres) then isside2 = .false. nres2 = nres2 - nres end if if(ntype == 3) then isside3 = .true. isprot3 = .true. nres3 = jgroup(nat3) if(nres3 < 0) then isprot3 = .false. nres3 = -nres3 end if if(nres3 > nres) then isside3 = .false. nres3 = nres3 - nres end if end if end if ! --- Echo result ! if(ntype.eq.1) then ! write(6,*) nat1,nres1,isside1, isprot1, fval ! else if(ntype.eq.2) then ! write(6,*) nat1,nres1,isside1, isprot1, & ! nat2,nres2,isside2, isprot2, & ! fval ! else if(ntype.eq.3) then ! write(6,*) nat1,nres1,isside1, isprot1, & ! nat2,nres2,isside2, isprot2, & ! nat3,nres3,isside3, isprot3, & ! fval ! end if ! --- Decompose ! here: nat1 solely determines if sidechain or backbone if(ntype == 1) then ! --- self-energy addfval = fval if (isside1) then nind1 = nssel else nind1 = nbsel end if else if(ntype >= 2) then if(ntype == 2) then addfval = fval else if(ntype == 3) then addfval = hfval end if if(nres2 == nres1) then ! --- self-energy if (isside1) then nind2 = nssel else nind2 = nbsel end if else if (( isprot1 .and. isprot2) .or. & (.not.isprot1 .and. .not.isprot2)) then ! --- indirect energy if (isside1) then nind2 = nsind else nind2 = nbind end if else ! --- direct energy if (isside1) then nind2 = nsdir else nind2 = nbdir end if end if if(ntype == 3) then if(nres3 == nres1) then ! --- self-energy if (isside1) then nind3 = nssel else nind3 = nbsel end if else if (( isprot1 .and. isprot3) .or. & (.not.isprot1 .and. .not.isprot3)) then ! --- indirect energy if (isside1) then nind3 = nsind else nind3 = nbind end if else ! --- direct energy if (isside1) then nind3 = nsdir else nind3 = nbdir end if end if end if end if ! --- Accumulate results if((ntype == 1 .and. nres1 > 0) .or. & (ntype == 2 .and. nres1 > 0 .and. nres2 > 0) .or. & (ntype == 3 .and. & nres1 > 0 .and. nres2 > 0 .and. nres3 > 0)) then ! --- Interactions "wanted" to be considered ! here: all contributions are collected for nres1(!) if(ntype == 1) then if(pairwise) then ntmp = (index(nres1) - 1) * npdec + index(nres1) dec(nind1 + ntmp - 1) = dec(nind1 + ntmp - 1) + addfval else dec(nind1 + nres1 - 1) = dec(nind1 + nres1 - 1) + addfval end if else if(ntype >= 2) then if(pairwise) then ntmp = (index(nres1) - 1) * npdec + index(nres2) dec(nind2 + ntmp - 1) = dec(nind2 + ntmp - 1) + addfval else dec(nind2 + nres1 - 1) = dec(nind2 + nres1 - 1) + addfval end if if(ntype == 3) then if(pairwise) then ntmp = (index(nres1) - 1) * npdec + index(nres3) dec(nind3 + ntmp - 1) = dec(nind3 + ntmp - 1) + addfval else dec(nind3 + nres1 - 1) = dec(nind3 + nres1 - 1) + addfval end if end if end if else ! --- Collect as "rest" if(ntype == 1) then dec(nind1 + ndecno) = dec(nind1 + ndecno) + addfval else if(ntype >= 2) then dec(nind2 + ndecno) = dec(nind2 + ndecno) + addfval if(ntype == 3) then dec(nind3 + ndecno) = dec(nind3 + ndecno) + addfval end if end if end if ! ((ntype == 1 .and. nres1 > 0) .or. return end subroutine decsasa !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine checkdec here] subroutine checkdec(idecomp) ! Sums contents in ENE_XXX_SEL/IND/DIR ! Holger Gohlke ! 13.11.2001 use memory_module, only : npdec, nres implicit none integer, intent(in) :: idecomp _REAL_ & esintsel, esintind, esvdwsel, esvdwind, esvdwdir, & eseelsel, eseelind, eseeldir, espolsel, espolind, & espoldir, essassel, essasind, essasdir _REAL_ & ebintsel, ebintind, ebvdwsel, ebvdwind, ebvdwdir, & ebeelsel, ebeelind, ebeeldir, ebpolsel, ebpolind, & ebpoldir, ebsassel, ebsasind, ebsasdir integer i, ndec integer nsintsel, nsintind, nsvdwsel, nsvdwind, nsvdwdir, & nseelsel, nseelind, nseeldir, nspolsel, nspolind, & nspoldir, nssassel, nssasind, nssasdir integer nbintsel, nbintind, nbvdwsel, nbvdwind, nbvdwdir, & nbeelsel, nbeelind, nbeeldir, nbpolsel, nbpolind, & nbpoldir, nbsassel, nbsasind, nbsasdir nsintsel = ndecind(sideintsel) nsintind = ndecind(sideintind) nsvdwsel = ndecind(sidevdwsel) nsvdwind = ndecind(sidevdwind) nsvdwdir = ndecind(sidevdwdir) nseelsel = ndecind(sideeelsel) nseelind = ndecind(sideeelind) nseeldir = ndecind(sideeeldir) nspolsel = ndecind(sidepolsel) nspolind = ndecind(sidepolind) nspoldir = ndecind(sidepoldir) nssassel = ndecind(sidesassel) nssasind = ndecind(sidesasind) nssasdir = ndecind(sidesasdir) nbintsel = ndecind(backintsel) nbintind = ndecind(backintind) nbvdwsel = ndecind(backvdwsel) nbvdwind = ndecind(backvdwind) nbvdwdir = ndecind(backvdwdir) nbeelsel = ndecind(backeelsel) nbeelind = ndecind(backeelind) nbeeldir = ndecind(backeeldir) nbpolsel = ndecind(backpolsel) nbpolind = ndecind(backpolind) nbpoldir = ndecind(backpoldir) nbsassel = ndecind(backsassel) nbsasind = ndecind(backsasind) nbsasdir = ndecind(backsasdir) esintsel = 0.0d0 esintind = 0.0d0 esvdwsel = 0.0d0 esvdwind = 0.0d0 esvdwdir = 0.0d0 eseelsel = 0.0d0 eseelind = 0.0d0 eseeldir = 0.0d0 espolsel = 0.0d0 espolind = 0.0d0 espoldir = 0.0d0 essassel = 0.0d0 essasind = 0.0d0 essasdir = 0.0d0 ebintsel = 0.0d0 ebintind = 0.0d0 ebvdwsel = 0.0d0 ebvdwind = 0.0d0 ebvdwdir = 0.0d0 ebeelsel = 0.0d0 ebeelind = 0.0d0 ebeeldir = 0.0d0 ebpolsel = 0.0d0 ebpolind = 0.0d0 ebpoldir = 0.0d0 ebsassel = 0.0d0 ebsasind = 0.0d0 ebsasdir = 0.0d0 if(idecomp < 3) then ndec = nres else if(idecomp > 2) then ndec = npdec*npdec endif do i=0,ndec-1 ! here: not till nres -> keeps "rest" apart esintsel = esintsel + dec(nsintsel + i) esintind = esintind + dec(nsintind + i) esvdwsel = esvdwsel + dec(nsvdwsel + i) esvdwind = esvdwind + dec(nsvdwind + i) esvdwdir = esvdwdir + dec(nsvdwdir + i) eseelsel = eseelsel + dec(nseelsel + i) eseelind = eseelind + dec(nseelind + i) eseeldir = eseeldir + dec(nseeldir + i) espolsel = espolsel + dec(nspolsel + i) espolind = espolind + dec(nspolind + i) espoldir = espoldir + dec(nspoldir + i) essassel = essassel + dec(nssassel + i) essasind = essasind + dec(nssasind + i) essasdir = essasdir + dec(nssasdir + i) ebintsel = ebintsel + dec(nbintsel + i) ebintind = ebintind + dec(nbintind + i) ebvdwsel = ebvdwsel + dec(nbvdwsel + i) ebvdwind = ebvdwind + dec(nbvdwind + i) ebvdwdir = ebvdwdir + dec(nbvdwdir + i) ebeelsel = ebeelsel + dec(nbeelsel + i) ebeelind = ebeelind + dec(nbeelind + i) ebeeldir = ebeeldir + dec(nbeeldir + i) ebpolsel = ebpolsel + dec(nbpolsel + i) ebpolind = ebpolind + dec(nbpolind + i) ebpoldir = ebpoldir + dec(nbpoldir + i) ebsassel = ebsassel + dec(nbsassel + i) ebsasind = ebsasind + dec(nbsasind + i) ebsasdir = ebsasdir + dec(nbsasdir + i) end do ! --- Output total energies (w/ rest) write(6,300) write(6,350) esintsel + esintind + ebintsel + ebintind + & dec(nsintsel+ndecno) + dec(nsintind+ndecno) + & dec(nbintsel+ndecno) + dec(nbintind+ndecno) write(6,360) esvdwsel + esvdwind + esvdwdir + & ebvdwsel + ebvdwind + ebvdwdir + & dec(nsvdwsel+ndecno) + dec(nsvdwind+ndecno) + & dec(nsvdwdir+ndecno) + & dec(nbvdwsel+ndecno) + dec(nbvdwind+ndecno) + & dec(nbvdwdir+ndecno), & eseelsel + eseelind + eseeldir + & ebeelsel + ebeelind + ebeeldir + & dec(nseelsel+ndecno) + dec(nseelind+ndecno) + & dec(nseeldir+ndecno) + & dec(nbeelsel+ndecno) + dec(nbeelind+ndecno) + & dec(nbeeldir+ndecno) write(6,370) espolsel + espolind + espoldir + & ebpolsel + ebpolind + ebpoldir + & dec(nspolsel+ndecno) + dec(nspolind+ndecno) + & dec(nspoldir+ndecno) + & dec(nbpolsel+ndecno) + dec(nbpolind+ndecno) + & dec(nbpoldir+ndecno), & essassel + essasind + essasdir + & ebsassel + ebsasind + ebsasdir + & dec(nssassel+ndecno) + dec(nssasind+ndecno) + & dec(nssasdir+ndecno) + & dec(nbsassel+ndecno) + dec(nbsasind+ndecno) + & dec(nbsasdir+ndecno) ! --- Output self energies (w/o rest) write(6,302) write(6,350) esintsel + ebintsel write(6,360) esvdwsel + ebvdwsel, & eseelsel + ebeelsel write(6,370) espolsel + ebpolsel, & essassel + ebsassel ! --- Output indirect energies (w/o rest) write(6,304) write(6,350) esintind + ebintind write(6,360) esvdwind + ebvdwind, & eseelind + ebeelind write(6,370) espolind + ebpolind, & essasind + ebsasind ! --- Output direct energies (w/o rest) write(6,306) write(6,350) 0.0d0 ! No direct interactions for internal energies write(6,360) esvdwdir + ebvdwdir, & eseeldir + ebeeldir write(6,370) espoldir + ebpoldir, & essasdir + ebsasdir ! --- Output rest energies write(6,308) write(6,350) dec(nsintsel+ndecno) + dec(nsintind+ndecno) + & dec(nbintsel+ndecno) + dec(nbintind+ndecno) write(6,360) dec(nsvdwsel+ndecno) + dec(nsvdwind+ndecno) + & dec(nsvdwdir+ndecno) + & dec(nbvdwsel+ndecno) + dec(nbvdwind+ndecno) + & dec(nbvdwdir+ndecno), & dec(nseelsel+ndecno) + dec(nseelind+ndecno) + & dec(nseeldir+ndecno) + & dec(nbeelsel+ndecno) + dec(nbeelind+ndecno) + & dec(nbeeldir+ndecno) write(6,370) dec(nspolsel+ndecno) + dec(nspolind+ndecno) + & dec(nspoldir+ndecno) + & dec(nbpolsel+ndecno) + dec(nbpolind+ndecno) + & dec(nbpoldir+ndecno), & dec(nssassel+ndecno) + dec(nssasind+ndecno) + & dec(nssasdir+ndecno) + & dec(nbsassel+ndecno) + dec(nbsasind+ndecno) + & dec(nbsasdir+ndecno) write(6,'(/)') 300 format(/ /20x,'CHECK DECOMP - TOTAL ENERGIES (w/ REST)',/) 302 format(/ /20x,'CHECK DECOMP - SELF ENERGIES (w/o REST)',/) 304 format(/ /20x,'CHECK DECOMP - INDIRECT ENERGIES (w/o REST)',/) 306 format(/ /20x,'CHECK DECOMP - DIRECT ENERGIES (w/o REST)',/) 308 format(/ /20x,'CHECK DECOMP - REST ENERGIES',/) 350 format(1x,'INTERNAL= ',f13.4) 360 format(1x,'VDWAALS = ',f13.4,2x,'EEL = ',f13.4) 370 format(1x,'EGB = ',f13.4,2x,'ESURF = ',f13.4) return end subroutine checkdec !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine printdec here] subroutine printdec(ix) ! Outputs contents in ENE_XXX_SEL/IND/DIR with respect to residues ! Holger Gohlke ! 13.11.2001 use memory_module, only : i02, ibellygp, nres implicit none _REAL_ energy dimension energy(5) integer i, j, ipresst, igst, istart integer ix(*) integer nssel, nsind, nsdir integer nbsel, nbind, nbdir dimension nssel(5), nsind(5), nsdir(5) dimension nbsel(5), nbind(5), nbdir(5) ipresst = i02 igst = ibellygp nssel(1) = ndecind(sideintsel) nsind(1) = ndecind(sideintind) nsdir(1) = 0 nssel(2) = ndecind(sidevdwsel) nsind(2) = ndecind(sidevdwind) nsdir(2) = ndecind(sidevdwdir) nssel(3) = ndecind(sideeelsel) nsind(3) = ndecind(sideeelind) nsdir(3) = ndecind(sideeeldir) nssel(4) = ndecind(sidepolsel) nsind(4) = ndecind(sidepolind) nsdir(4) = ndecind(sidepoldir) nssel(5) = ndecind(sidesassel) nsind(5) = ndecind(sidesasind) nsdir(5) = ndecind(sidesasdir) nbsel(1) = ndecind(backintsel) nbind(1) = ndecind(backintind) nbdir(1) = 0 nbsel(2) = ndecind(backvdwsel) nbind(2) = ndecind(backvdwind) nbdir(2) = ndecind(backvdwdir) nbsel(3) = ndecind(backeelsel) nbind(3) = ndecind(backeelind) nbdir(3) = ndecind(backeeldir) nbsel(4) = ndecind(backpolsel) nbind(4) = ndecind(backpolind) nbdir(4) = ndecind(backpoldir) nbsel(5) = ndecind(backsassel) nbind(5) = ndecind(backsasind) nbdir(5) = ndecind(backsasdir) ! --- Output total energies write(6,300) write(6,350) write(6,355) do i=0,nres-1 istart = ix(ipresst + i) if(ix(igst + istart - 1) > 0) then do j=1,5 energy(j) = 0.0d0 if(nssel(j) > 0 .and. nbsel(j) > 0) & energy(j) = energy(j) + dec(nssel(j)+i) + dec(nbsel(j)+i) if(nsind(j) > 0 .and. nbind(j) > 0) & energy(j) = energy(j) + dec(nsind(j)+i) + dec(nbind(j)+i) if(nsdir(j) > 0 .and. nbdir(j) > 0) & energy(j) = energy(j) + dec(nsdir(j)+i) + dec(nbdir(j)+i) end do write(6,360) i+1,(energy(j),j=1,5) end if end do ! --- Output sidechain energies write(6,302) write(6,350) write(6,355) do i=0,nres-1 istart = ix(ipresst + i) if(ix(igst + istart - 1) > 0) then do j=1,5 energy(j) = 0.0d0 if(nssel(j) > 0) & energy(j) = energy(j) + dec(nssel(j)+i) if(nsind(j) > 0) & energy(j) = energy(j) + dec(nsind(j)+i) if(nsdir(j) > 0) & energy(j) = energy(j) + dec(nsdir(j)+i) end do write(6,362) i+1,(energy(j),j=1,5) end if end do ! --- Output backbone energies write(6,304) write(6,350) write(6,355) do i=0,nres-1 istart = ix(ipresst + i) if(ix(igst + istart - 1) > 0) then do j=1,5 energy(j) = 0.0d0 if(nbsel(j) > 0) & energy(j) = energy(j) + dec(nbsel(j)+i) if(nbind(j) > 0) & energy(j) = energy(j) + dec(nbind(j)+i) if(nbdir(j) > 0) & energy(j) = energy(j) + dec(nbdir(j)+i) end do write(6,364) i+1,(energy(j),j=1,5) end if end do write(6,'(/)') 300 format(/ /20x,'PRINT DECOMP - TOTAL ENERGIES',/) 302 format(/ /20x,'PRINT DECOMP - SIDECHAIN ENERGIES',/) 304 format(/ /20x,'PRINT DECOMP - BACKBONE ENERGIES',/) 350 format(3x,1x,'resid |internal |vdw |eel |pol |sas') 355 format(6('==========')) 360 format('TDC',1x,i6,5(1x,f9.3)) 362 format('SDC',1x,i6,5(1x,f9.3)) 364 format('BDC',1x,i6,5(1x,f9.3)) return end subroutine printdec !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine printpdec here] subroutine printpdec(ix) ! Outputs contents in ENE_XXX_SEL/IND/DIR with ! respect to pairwise decomp ! Holger Gohlke ! 29.11.2001 use memory_module, only : i02, ibellygp, npdec, nres implicit none _REAL_ energy dimension energy(5) integer i, j, k, ipresst, igst, istart integer ix(*) integer nssel, nsind, nsdir integer nbsel, nbind, nbdir integer ntmp dimension nssel(5), nsind(5), nsdir(5) dimension nbsel(5), nbind(5), nbdir(5) ipresst = i02 igst = ibellygp nssel(1) = ndecind(sideintsel) nsind(1) = ndecind(sideintind) nsdir(1) = 0 nssel(2) = ndecind(sidevdwsel) nsind(2) = ndecind(sidevdwind) nsdir(2) = ndecind(sidevdwdir) nssel(3) = ndecind(sideeelsel) nsind(3) = ndecind(sideeelind) nsdir(3) = ndecind(sideeeldir) nssel(4) = ndecind(sidepolsel) nsind(4) = ndecind(sidepolind) nsdir(4) = ndecind(sidepoldir) nssel(5) = ndecind(sidesassel) nsind(5) = ndecind(sidesasind) nsdir(5) = ndecind(sidesasdir) nbsel(1) = ndecind(backintsel) nbind(1) = ndecind(backintind) nbdir(1) = 0 nbsel(2) = ndecind(backvdwsel) nbind(2) = ndecind(backvdwind) nbdir(2) = ndecind(backvdwdir) nbsel(3) = ndecind(backeelsel) nbind(3) = ndecind(backeelind) nbdir(3) = ndecind(backeeldir) nbsel(4) = ndecind(backpolsel) nbind(4) = ndecind(backpolind) nbdir(4) = ndecind(backpoldir) nbsel(5) = ndecind(backsassel) nbind(5) = ndecind(backsasind) nbdir(5) = ndecind(backsasdir) ! --- Output total energies print *, dec(186) write(6,300) write(6,350) write(6,355) do i=0,nres-1 istart = ix(ipresst + i) if(ix(igst + istart - 1) > 0 .and. & index(i + 1) > 0) then do k=0,nres-1 istart = ix(ipresst + k) if(ix(igst + istart - 1) > 0 .and. index(k + 1) > 0) then ntmp = (index(i + 1) - 1) * npdec + index(k + 1) - 1 do j=1,5 energy(j) = 0.0d0 if(nssel(j) > 0 .and. nbsel(j) > 0) & energy(j) = energy(j) + & dec(nssel(j)+ntmp) + dec(nbsel(j)+ntmp) if(nsind(j) > 0 .and. nbind(j) > 0) & energy(j) = energy(j) + & dec(nsind(j)+ntmp) + dec(nbind(j)+ntmp) if(nsdir(j) > 0 .and. nbdir(j) > 0) & energy(j) = energy(j) + & dec(nsdir(j)+ntmp) + dec(nbdir(j)+ntmp) end do write(6,360) i+1,k+1,(energy(j),j=1,5) end if end do end if end do ! --- Output sidechain energies write(6,302) write(6,350) write(6,355) do i=0,nres-1 istart = ix(ipresst + i) if(ix(igst + istart - 1) > 0 .and. index(i + 1) > 0) then do k=0,nres-1 istart = ix(ipresst + k) if(ix(igst + istart - 1) > 0 .and. index(k + 1) > 0) then ntmp = (index(i + 1) - 1) * npdec + index(k + 1) - 1 do j=1,5 energy(j) = 0.0d0 if(nssel(j) > 0) & energy(j) = energy(j) + dec(nssel(j)+ntmp) if(nsind(j) > 0) & energy(j) = energy(j) + dec(nsind(j)+ntmp) if(nsdir(j) > 0) & energy(j) = energy(j) + dec(nsdir(j)+ntmp) end do write(6,362) i+1,k+1,(energy(j),j=1,5) end if end do end if end do ! --- Output backbone energies write(6,304) write(6,350) write(6,355) do i=0,nres-1 istart = ix(ipresst + i) if(ix(igst + istart - 1) > 0 .and. index(i + 1) > 0) then do k=0,nres-1 istart = ix(ipresst + k) if(ix(igst + istart - 1) > 0 .and. index(k + 1) > 0) then ntmp = (index(i + 1) - 1) * npdec + index(k + 1) - 1 do j=1,5 energy(j) = 0.0d0 if(nbsel(j) > 0) & energy(j) = energy(j) + dec(nbsel(j)+ntmp) if(nbind(j) > 0) & energy(j) = energy(j) + dec(nbind(j)+ntmp) if(nbdir(j) > 0) & energy(j) = energy(j) + dec(nbdir(j)+ntmp) end do write(6,364) i+1,k+1,(energy(j),j=1,5) end if end do end if end do write(6,'(/)') 300 format(/ /20x,'PRINT PAIR DECOMP - TOTAL ENERGIES',/) 302 format(/ /20x,'PRINT PAIR DECOMP - SIDECHAIN ENERGIES',/) 304 format(/ /20x,'PRINT PAIR DECOMP - BACKBONE ENERGIES',/) 350 format(3x,1x,'resid1 ->resid2 ', & '|internal |vdw |eel |pol |sas') 355 format(7('==========')) 360 format('TDC',1x,i7,'->',i7,5(1x,f9.3)) 362 format('SDC',1x,i7,'->',i7,5(1x,f9.3)) 364 format('BDC',1x,i7,'->',i7,5(1x,f9.3)) return end subroutine printpdec end module decomp