homology dihederal i-3 gradient correction src_MD

[unres.git] / source / unres / src_MD / energy_p_new_barrier.F

diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F
index 8943347..efa39f2 100644 (file)
--- a/source/unres/src_MD/energy_p_new_barrier.F
+++ b/source/unres/src_MD/energy_p_new_barrier.F
@@ -263,7 +263,7 @@ cd    print *,'nterm=',nterm
        edihcnstr=0
       endif
 
-      if (constr_homology.ge.1) then
+      if (constr_homology.ge.1.and.waga_homology(iset).ne.0d0) then
         call e_modeller(ehomology_constr)
 c        print *,'iset=',iset,'me=',me,ehomology_constr,
 c     &  'Processor',fg_rank,' CG group',kolor,
@@ -822,7 +822,7 @@ c      enddo
 #endif
         enddo
       enddo 
-      if (constr_homology.gt.0) then
+      if (constr_homology.gt.0.and.waga_homology(iset).ne.0d0) then
         do i=1,nct
           do j=1,3
             gradc(j,i,icg)=gradc(j,i,icg)+duscdiff(j,i)
@@ -1781,9 +1781,10 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
      &        evdwij
             endif
 
-            if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
-     &                        'evdw',i,j,evdwij
-
+            if (energy_dec) then
+              write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij
+              call flush(iout)
+            endif
 C Calculate gradient components.
             e1=e1*eps1*eps2rt**2*eps3rt**2
             fac=-expon*(e1+evdwij)*rij_shift
@@ -4360,6 +4361,7 @@ C
       include 'COMMON.VAR'
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
+      include 'COMMON.CONTROL'
       dimension ggg(3)
       ehpb=0.0D0
 cd      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
@@ -4395,6 +4397,16 @@ cd          write (iout,*) "eij",eij
         else if (ii.gt.nres .and. jj.gt.nres) then
 c Restraints from contact prediction
           dd=dist(ii,jj)
+          if (constr_dist.eq.11) then
+            ehpb=ehpb+fordepth(i)**4.0d0
+     &          *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+            fac=fordepth(i)**4.0d0
+     &          *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+          if (energy_dec) write (iout,'(a6,2i5,f15.6,2f8.3)') 
+     &     "edisl",ii,jj,
+     &     fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)),
+     &     fordepth(i),dd
+           else
           if (dhpb1(i).gt.0.0d0) then
             ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
             fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
@@ -4412,7 +4424,8 @@ C
 C Evaluate gradient.
 C
             fac=waga*rdis/dd
-          endif  
+          endif
+          endif
           do j=1,3
             ggg(j)=fac*(c(j,jj)-c(j,ii))
           enddo
@@ -4428,6 +4441,18 @@ C
 C Calculate the distance between the two points and its difference from the
 C target distance.
           dd=dist(ii,jj)
+         if (constr_dist.eq.11) then
+            ehpb=ehpb+fordepth(i)**4.0d0
+     &           *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+            fac=fordepth(i)**4.0d0
+     &           *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+          if (energy_dec) write (iout,'(a6,2i5,f15.6,2f8.3)') 
+     7     "edisl",ii,jj,
+     &     fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)),
+     &     fordepth(i),dd
+c          if (energy_dec)
+c     &      write (iout,*) fac
+         else   
           if (dhpb1(i).gt.0.0d0) then
             ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
             fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
@@ -4445,6 +4470,7 @@ C Evaluate gradient.
 C
             fac=waga*rdis/dd
           endif
+         endif
 cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
 cd   &   ' waga=',waga,' fac=',fac
             do j=1,3
@@ -4470,7 +4496,10 @@ cgrad        enddo
           enddo
         endif
       enddo
-      ehpb=0.5D0*ehpb
+      if (constr_dist.ne.11) ehpb=0.5D0*ehpb
+c      do i=1,nres
+c        write (iout,*) "ghpbc",i,(ghpbc(j,i),j=1,3)
+c      enddo
       return
       end
 C--------------------------------------------------------------------------
@@ -4584,6 +4613,8 @@ c
       do i=ibondp_start,ibondp_end
         diff = vbld(i)-vbldp0
 c        write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
+        if (energy_dec)    write (iout,'(a7,i5,4f7.3)') 
+     &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
         estr=estr+diff*diff
         do j=1,3
           gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
@@ -4602,6 +4633,12 @@ c
             diff=vbld(i+nres)-vbldsc0(1,iti)
 c            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
 c     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
+            if (energy_dec)  then
+              write (iout,*) 
+     &         "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+     &         AKSC(1,iti),AKSC(1,iti)*diff*diff
+              call flush(iout)
+            endif
             estr=estr+0.5d0*AKSC(1,iti)*diff*diff
             do j=1,3
               gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
@@ -4894,7 +4931,7 @@ c     &     " ithet_end",ithet_end
           sinkt(k)=dsin(k*theti2)
         enddo
 C        if (i.gt.3) then
-         if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+        if (i.gt.3 .and. itype(imax0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -6020,7 +6057,7 @@ c
 c
 
 
-      do i=1,19
+      do i=1,max_template
         distancek(i)=9999999.9
       enddo
 
@@ -6061,7 +6098,18 @@ c
            endif
          enddo
          
-         min_odl=minval(distancek)
+
+c         min_odl=minval(distancek)
+         do kk=1,constr_homology
+          if(l_homo(kk,ii)) then 
+            min_odl=distancek(kk)
+            exit
+          endif
+         enddo
+         do kk=1,constr_homology
+          if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) 
+     &              min_odl=distancek(kk)
+         enddo
 c        write (iout,* )"min_odl",min_odl
 #ifdef DEBUG
          write (iout,*) "ij dij",i,j,dij
@@ -6214,7 +6262,7 @@ 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)
+        betai = phi(i)
 c       write (iout,*) "betai =",betai
         do k=1,constr_homology
           dih_diff(k)=pinorm(dih(k,i)-betai)
@@ -6261,7 +6309,7 @@ c       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
+        gloc(i-3,icg)=gloc(i-3,icg)+grad_dih3
 c        if (i.eq.25) then
 c        write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg)
 c        endif