2D replica exchange with homology constraints
[unres.git] / source / unres / src_MD / energy_p_new_barrier.F
index 21ac2a7..7cf67db 100644 (file)
@@ -27,6 +27,7 @@ cMS$ATTRIBUTES C ::  proc_proc
 #ifdef MPI      
 c      print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
 c     & " nfgtasks",nfgtasks
+      call flush(iout)
       if (nfgtasks.gt.1) then
 #ifdef MPI
         time00=MPI_Wtime()
@@ -131,6 +132,36 @@ C
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
   107 continue
+C     BARTEK for dfa test!
+      if (wdfa_dist.gt.0) then 
+        call edfad(edfadis)
+      else
+        edfadis=0
+      endif
+c      print*, 'edfad is finished!', edfadis
+      if (wdfa_tor.gt.0) then
+        call edfat(edfator)
+      else
+        edfator=0
+      endif
+c      print*, 'edfat is finished!', edfator
+      if (wdfa_nei.gt.0) then
+        call edfan(edfanei)
+      else
+        edfanei=0
+      endif    
+c      print*, 'edfan is finished!', edfanei
+      if (wdfa_beta.gt.0) then 
+        call edfab(edfabet)
+      else
+        edfabet=0
+      endif
+c      print*, 'edfab is finished!', edfabet
+cmc
+cmc Sep-06: egb takes care of dynamic ss bonds too
+cmc
+c      if (dyn_ss) call dyn_set_nss
+
 c      print *,"Processor",myrank," computed USCSC"
 #ifdef TIMING
 #ifdef MPI
@@ -223,6 +254,15 @@ cd    print *,'nterm=',nterm
        etors=0
        edihcnstr=0
       endif
+
+      if (constr_homology.ge.1) then
+        call e_modeller(ehomology_constr)
+      else
+        ehomology_constr=0.0d0
+      endif
+
+
+c      write(iout,*) ehomology_constr
 c      print *,"Processor",myrank," computed Utor"
 C
 C 6/23/01 Calculate double-torsional energy
@@ -324,8 +364,14 @@ C
       energia(21)=esccor
       energia(22)=evdw_p
       energia(23)=evdw_m
+      energia(24)=ehomology_constr
+      energia(25)=edfadis
+      energia(26)=edfator
+      energia(27)=edfanei
+      energia(28)=edfabet
 c      print *," Processor",myrank," calls SUM_ENERGY"
       call sum_energy(energia,.true.)
+      if (dyn_ss) call dyn_set_nss
 c      print *," Processor",myrank," left SUM_ENERGY"
 #ifdef TIMING
 #ifdef MPI
@@ -420,20 +466,29 @@ cMS$ATTRIBUTES C ::  proc_proc
       estr=energia(17)
       Uconst=energia(20)
       esccor=energia(21)
+      ehomology_constr=energia(24)
+      edfadis=energia(25)
+      edfator=energia(26)
+      edfanei=energia(27)
+      edfabet=energia(28)
 #ifdef SPLITELE
       etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
      & +wang*ebe+wtor*etors+wscloc*escloc
-     & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5
+     & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
-     & +wbond*estr+Uconst+wsccor*esccor
+     & +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
+     & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
+     & +wdfa_beta*edfabet    
 #else
       etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
      & +wang*ebe+wtor*etors+wscloc*escloc
-     & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5
+     & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
-     & +wbond*estr+Uconst+wsccor*esccor
+     & +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
+     & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
+     & +wdfa_beta*edfabet    
 #endif
       energia(0)=etot
 c detecting NaNQ
@@ -540,7 +595,11 @@ c      enddo
      &                wcorr5*gradcorr5_long(j,i)+
      &                wcorr6*gradcorr6_long(j,i)+
      &                wturn6*gcorr6_turn_long(j,i)+
-     &                wstrain*ghpbc(j,i)
+     &                wstrain*ghpbc(j,i)+
+     &                wdfa_dist*gdfad(j,i)+
+     &                wdfa_tor*gdfat(j,i)+
+     &                wdfa_nei*gdfan(j,i)+
+     &                wdfa_beta*gdfab(j,i)
         enddo
       enddo 
 #else
@@ -554,7 +613,11 @@ c      enddo
      &                wcorr5*gradcorr5_long(j,i)+
      &                wcorr6*gradcorr6_long(j,i)+
      &                wturn6*gcorr6_turn_long(j,i)+
-     &                wstrain*ghpbc(j,i)
+     &                wstrain*ghpbc(j,i)+
+     &                wdfa_dist*gdfad(j,i)+
+     &                wdfa_tor*gdfat(j,i)+
+     &                wdfa_nei*gdfan(j,i)+
+     &                wdfa_beta*gdfab(j,i)
         enddo
       enddo 
 #endif
@@ -570,7 +633,11 @@ c      enddo
      &                wcorr5*gradcorr5_long(j,i)+
      &                wcorr6*gradcorr6_long(j,i)+
      &                wturn6*gcorr6_turn_long(j,i)+
-     &                wstrain*ghpbc(j,i)
+     &                wstrain*ghpbc(j,i)+
+     &                wdfa_dist*gdfad(j,i)+
+     &                wdfa_tor*gdfat(j,i)+
+     &                wdfa_nei*gdfan(j,i)+
+     &                wdfa_beta*gdfab(j,i)
         enddo
       enddo 
 #endif
@@ -1042,6 +1109,13 @@ C------------------------------------------------------------------------
       estr=energia(17)
       Uconst=energia(20)
       esccor=energia(21)
+      ehomology_constr=energia(24)
+C     Bartek
+      edfadis = energia(25)
+      edfator = energia(26)
+      edfanei = energia(27)
+      edfabet = energia(28)
+
 #ifdef SPLITELE
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
      &  estr,wbond,ebe,wang,
@@ -1049,31 +1123,37 @@ C------------------------------------------------------------------------
      &  ecorr,wcorr,
      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
-     &  edihcnstr,ebr*nss,
-     &  Uconst,etot
+     &  edihcnstr,ehomology_constr, ebr*nss,
+     &  Uconst,edfadis,wdfa_dist,edfator,wdfa_tor,edfanei,wdfa_nei,
+     &  edfabet,wdfa_beta,etot
    10 format (/'Virtual-chain energies:'//
-     & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
-     & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
-     & 'EES=   ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p)'/
-     & 'EVDWPP=',1pE16.6,' WEIGHT=',1pD16.6,' (p-p VDW)'/
-     & 'ESTR=  ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching)'/
-     & 'EBE=   ',1pE16.6,' WEIGHT=',1pD16.6,' (bending)'/
-     & 'ESC=   ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/
-     & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/
-     & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/
-     & 'EHPB=  ',1pE16.6,' WEIGHT=',1pD16.6,
+     & 'EVDW=  ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-SC)'/
+     & 'EVDW2= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-p)'/
+     & 'EES=   ',1pE16.6,' WEIGHT=',1pE16.6,' (p-p)'/
+     & 'EVDWPP=',1pE16.6,' WEIGHT=',1pE16.6,' (p-p VDW)'/
+     & 'ESTR=  ',1pE16.6,' WEIGHT=',1pE16.6,' (stretching)'/
+     & 'EBE=   ',1pE16.6,' WEIGHT=',1pE16.6,' (bending)'/
+     & 'ESC=   ',1pE16.6,' WEIGHT=',1pE16.6,' (SC local)'/
+     & 'ETORS= ',1pE16.6,' WEIGHT=',1pE16.6,' (torsional)'/
+     & 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/
+     & 'EHPB=  ',1pE16.6,' WEIGHT=',1pE16.6,
      & ' (SS bridges & dist. cnstr.)'/
-     & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
-     & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
-     & 'ECORR6=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
-     & 'EELLO= ',1pE16.6,' WEIGHT=',1pD16.6,' (electrostatic-local)'/
-     & 'ETURN3=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 3rd order)'/
-     & 'ETURN4=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 4th order)'/
-     & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
-     & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
+     & 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+     & 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+     & 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+     & 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/
+     & 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/
+     & 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/
+     & 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/
+     & 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/
      & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+     & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
      & 'UCONST= ',1pE16.6,' (Constraint energy)'/ 
+     & 'EDFAD= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA distance energy)'/ 
+     & 'EDFAT= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA torsion energy)'/ 
+     & 'EDFAN= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA NCa energy)'/ 
+     & 'EDFAB= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA Beta energy)'/ 
      & 'ETOT=  ',1pE16.6,' (total)')
 #else
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
@@ -1082,7 +1162,9 @@ C------------------------------------------------------------------------
      &  ecorr,wcorr,
      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
-     &  ebr*nss,Uconst,etot
+     &  ehomology_constr,ebr*nss,Uconst,edfadis,wdfa_dist,edfator,
+     &  wdfa_tor,edfanei,wdfa_nei,edfabet,wdfa_beta,
+     &  etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
@@ -1103,8 +1185,13 @@ C------------------------------------------------------------------------
      & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
      & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
      & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+     & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
      & 'UCONST=',1pE16.6,' (Constraint energy)'/ 
+     & 'EDFAD= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA distance energy)'/ 
+     & 'EDFAT= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA torsion energy)'/ 
+     & 'EDFAN= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA NCa energy)'/ 
+     & 'EDFAB= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA Beta energy)'/ 
      & 'ETOT=  ',1pE16.6,' (total)')
 #endif
       return
@@ -1552,6 +1639,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
       include 'COMMON.CONTROL'
+      include 'COMMON.SBRIDGE'
       logical lprn
       evdw=0.0D0
 ccccc      energy_dec=.false.
@@ -1580,6 +1668,12 @@ C Calculate SC interaction energy.
 C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
+            IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+              call dyn_ssbond_ene(i,j,evdwij)
+              evdw=evdw+evdwij
+              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') 
+     &                        'evdw',i,j,evdwij,' ss'
+            ELSE
             ind=ind+1
             itypj=itype(j)
 c            dscj_inv=dsc_inv(itypj)
@@ -1689,6 +1783,7 @@ C Calculate angular part of the gradient.
 #else
             call sc_grad
 #endif
+            ENDIF    ! dyn_ss            
           enddo      ! j
         enddo        ! iint
       enddo          ! i
@@ -4266,9 +4361,15 @@ c        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
 c     &    dhpb(i),dhpb1(i),forcon(i)
 C 24/11/03 AL: SS bridges handled separately because of introducing a specific
 C    distance and angle dependent SS bond potential.
-        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+cmc        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
+        if (.not.dyn_ss .and. i.le.nss) then
+C 15/02/13 CC dynamic SSbond - additional check
+         if (ii.gt.nres 
+     &       .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then 
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
+         endif
 cd          write (iout,*) "eij",eij
         else if (ii.gt.nres .and. jj.gt.nres) then
 c Restraints from contact prediction
@@ -4407,7 +4508,7 @@ c      dscj_inv=dsc_inv(itypj)
       deltat12=om2-om1+2.0d0
       cosphi=om12-om1*om2
       eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2)
-     &  +akct*deltad*deltat12
+     &  +akct*deltad*deltat12+ebr
      &  +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
 c      write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
 c     &  " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
@@ -4617,7 +4718,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2.
      &      'ebend',i,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
-        gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
+        gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
       enddo
 C Ufff.... We've done all this!!! 
       return
           enddo
         enddo
 10      continue
-        if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') 
-     &   i,theta(i)*rad2deg,phii*rad2deg,
+c        lprn1=.true.
+        if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)') 
+     &  'ebe', i,theta(i)*rad2deg,phii*rad2deg,
      &   phii1*rad2deg,ethetai
+c        lprn1=.false.
         etheta=etheta+ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
-        gloc(nphi+i-2,icg)=wang*dethetai
+        gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
       enddo
       return
       end
@@ -5738,6 +5841,15 @@ C Proline-Proline pair is a special case...
       return
       end
 c------------------------------------------------------------------------------
+c LICZENIE WIEZOW Z ROWNANIA ENERGII MODELLERA
+      subroutine e_modeller(ehomology_constr)
+      ehomology_constr=0.0
+      write (iout,*) "!!!!!UWAGA, JESTEM W DZIWNEJ PETLI, TEST!!!!!"
+      return
+      end
+C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
+
+c------------------------------------------------------------------------------
       subroutine etor_d(etors_d)
       etors_d=0.0d0
       return
@@ -5838,6 +5950,187 @@ cd       write (iout,*) 'edihcnstr',edihcnstr
       return
       end
 c----------------------------------------------------------------------------
+c MODELLER restraint function
+      subroutine e_modeller(ehomology_constr)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+
+      integer nnn, i, j, k, ki, irec, l
+      integer katy, odleglosci, test7
+      real*8 odleg, odleg2, odleg3, kat, kat2, kat3, gdih(max_template)
+      real*8 distance(max_template),distancek(max_template),
+     &    min_odl,godl(max_template),dih_diff(max_template)
+
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      include 'COMMON.DERIV'
+      include 'COMMON.LOCAL'
+      include 'COMMON.INTERACT'
+      include 'COMMON.VAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.MD'
+      include 'COMMON.CONTROL'
+
+
+      do i=1,19
+        distancek(i)=9999999.9
+      enddo
+
+
+      odleg=0.0d0
+
+c Pseudo-energy and gradient from homology restraints (MODELLER-like
+c function)
+C AL 5/2/14 - Introduce list of restraints
+      do ii = link_start_homo,link_end_homo
+         i = ires_homo(ii)
+         j = jres_homo(ii)
+         dij=dist(i,j)
+         do k=1,constr_homology
+           distance(k)=odl(k,ii)-dij
+           distancek(k)=
+     &        0.5d0*waga_dist(iset)*distance(k)**2*sigma_odl(k,ii)
+         enddo
+         
+         min_odl=minval(distancek)
+#ifdef DEBUG
+         write (iout,*) "ij dij",i,j,dij
+         write (iout,*) "distance",(distance(k),k=1,constr_homology)
+         write (iout,*) "distancek",(distancek(k),k=1,constr_homology)
+         write (iout,* )"min_odl",min_odl
+#endif
+         odleg2=0.0d0
+         do k=1,constr_homology
+c Nie wiem po co to liczycie jeszcze raz!
+c            odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ 
+c     &              (2*(sigma_odl(i,j,k))**2))
+            godl(k)=dexp(-distancek(k)+min_odl)
+            odleg2=odleg2+godl(k)
+
+ccc       write(iout,779) i,j,k, "odleg2=",odleg2, "odleg3=", odleg3,
+ccc     & "dEXP(odleg3)=", dEXP(odleg3),"distance(i,j,k)^2=",
+ccc     & distance(i,j,k)**2, "dist(i+1,j+1)=", dist(i+1,j+1),
+ccc     & "sigma_odl(i,j,k)=", sigma_odl(i,j,k)
+
+         enddo
+#ifdef DEBUG
+         write (iout,*) "godl",(godl(k),k=1,constr_homology)
+         write (iout,*) "ii i j",ii,i,j," odleg2",odleg2
+#endif
+         odleg=odleg-dLOG(odleg2/constr_homology)+min_odl
+c Gradient
+         sum_godl=odleg2
+         sum_sgodl=0.0
+         do k=1,constr_homology
+c            godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+c     &           *waga_dist(iset))+min_odl
+           sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist(iset)
+           sum_sgodl=sum_sgodl+sgodl
+
+c            sgodl2=sgodl2+sgodl
+c      write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE1"
+c      write(iout,*) "constr_homology=",constr_homology
+c      write(iout,*) i, j, k, "TEST K"
+         enddo
+
+         grad_odl3=sum_sgodl/(sum_godl*dij)
+
+
+c      write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE2"
+c      write(iout,*) (distance(i,j,k)**2), (2*(sigma_odl(i,j,k))**2),
+c     &              (-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+
+ccc      write(iout,*) godl, sgodl, grad_odl3
+
+c          grad_odl=grad_odl+grad_odl3
+
+         do jik=1,3
+            ggodl=grad_odl3*(c(jik,i)-c(jik,j))
+ccc      write(iout,*) c(jik,i+1), c(jik,j+1), (c(jik,i+1)-c(jik,j+1))
+ccc      write(iout,746) "GRAD_ODL_1", i, j, jik, ggodl, 
+ccc     &              ghpbc(jik,i+1), ghpbc(jik,j+1)
+            ghpbc(jik,i)=ghpbc(jik,i)+ggodl
+            ghpbc(jik,j)=ghpbc(jik,j)-ggodl
+ccc      write(iout,746) "GRAD_ODL_2", i, j, jik, ggodl,
+ccc     &              ghpbc(jik,i+1), ghpbc(jik,j+1)
+
+         enddo
+ccc       write(iout,778)"TEST: odleg2=", odleg2, "DLOG(odleg2)=", 
+ccc     & dLOG(odleg2),"-odleg=", -odleg
+
+      enddo ! ii
+c Pseudo-energy and gradient from dihedral-angle restraints from
+c homology templates
+c      write (iout,*) "End of distance loop"
+c      call flush(iout)
+      kat=0.0d0
+c      write (iout,*) idihconstr_start_homo,idihconstr_end_homo
+      do i=idihconstr_start_homo,idihconstr_end_homo
+        kat2=0.0d0
+c        betai=beta(i,i+1,i+2,i+3)
+        betai = phi(i+3)
+        do k=1,constr_homology
+          dih_diff(k)=pinorm(dih(k,i)-betai)
+c          if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)=
+c     &                                   -(6.28318-dih_diff(i,k))
+c          if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
+c     &                                   6.28318+dih_diff(i,k)
+
+          kat3=-0.5d0*waga_angle(iset)*dih_diff(k)**2*sigma_dih(k,i)
+          gdih(k)=dexp(kat3)
+          kat2=kat2+gdih(k)
+c          write(iout,*) "kat2=", kat2, "exp(kat3)=", exp(kat3)
+c          write(*,*)""
+        enddo
+#ifdef DEBUG
+        write (iout,*) "i",i," betai",betai," kat2",kat2
+        write (iout,*) "gdih",(gdih(k),k=1,constr_homology)
+#endif
+        if (kat2.le.1.0d-14) cycle
+        kat=kat-dLOG(kat2/constr_homology)
+
+ccc       write(iout,778)"TEST: kat2=", kat2, "DLOG(kat2)=",
+ccc     & dLOG(kat2), "-kat=", -kat
+
+c ----------------------------------------------------------------------
+c Gradient
+c ----------------------------------------------------------------------
+
+        sum_gdih=kat2
+        sum_sgdih=0.0
+        do k=1,constr_homology
+          sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle(iset)
+          sum_sgdih=sum_sgdih+sgdih
+        enddo
+        grad_dih3=sum_sgdih/sum_gdih
+
+c      write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3
+ccc      write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3,
+ccc     & gloc(nphi+i-3,icg)
+        gloc(i,icg)=gloc(i,icg)+grad_dih3
+ccc      write(iout,747) "GRAD_KAT_2", i, nphi, icg, grad_dih3,
+ccc     & gloc(nphi+i-3,icg)
+
+      enddo
+
+
+c Total energy from homology restraints
+#ifdef DEBUG
+      write (iout,*) "odleg",odleg," kat",kat
+#endif
+      ehomology_constr=odleg+kat
+      return
+
+  748 format(a8,f12.3,a6,f12.3,a7,f12.3)
+  747 format(a12,i4,i4,i4,f8.3,f8.3)
+  746 format(a12,i4,i4,i4,f8.3,f8.3,f8.3)
+  778 format(a7,1X,f10.3,1X,a4,1X,f10.3,1X,a5,1X,f10.3)
+  779 format(i3,1X,i3,1X,i2,1X,a7,1X,f7.3,1X,a7,1X,f7.3,1X,a13,1X,
+     &       f7.3,1X,a17,1X,f9.3,1X,a10,1X,f8.3,1X,a10,1X,f8.3)
+      end
+
+c------------------------------------------------------------------------------
       subroutine etor_d(etors_d)
 C 6/23/01 Compute double torsional energy
       implicit real*8 (a-h,o-z)