WHAM for Lorentzian like constrains fixed - DEBUG OFF
[unres.git] / source / wham / src-M / energy_p_new.F
index beae087..8cf29b2 100644 (file)
@@ -819,8 +819,8 @@ C
             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
               call dyn_ssbond_ene(i,j,evdwij)
               evdw=evdw+evdwij
-            write (iout,'(a6,2i5,0pf7.3,a3,2f10.3)')
-     &                        'evdw',i,j,evdwij,' ss',evdw,evdw_t
+C            write (iout,'(a6,2i5,0pf7.3,a3,2f10.3)')
+C     &                        'evdw',i,j,evdwij,' ss',evdw,evdw_t
 C triple bond artifac removal
              do k=j+1,iend(i,iint)
 C search over all next residues
@@ -831,8 +831,8 @@ C              write(iout,*) 'k=',k
 C call the energy function that removes the artifical triple disulfide
 C bond the soubroutine is located in ssMD.F
               evdw=evdw+evdwij
-             write (iout,'(a6,2i5,0pf7.3,a3,2f10.3)')
-     &                        'evdw',i,j,evdwij,'tss',evdw,evdw_t
+C             write (iout,'(a6,2i5,0pf7.3,a3,2f10.3)')
+C     &                        'evdw',i,j,evdwij,'tss',evdw,evdw_t
               endif!dyn_ss_mask(k)
              enddo! k
             ELSE
@@ -929,7 +929,7 @@ C Calculate the radial part of the gradient
 C Calculate angular part of the gradient.
             call sc_grad
             endif
-            write(iout,*)  "partial sum", evdw, evdw_t
+C            write(iout,*)  "partial sum", evdw, evdw_t
             ENDIF    ! dyn_ss            
           enddo      ! j
         enddo        ! iint
       include 'COMMON.DERIV'
       include 'COMMON.VAR'
       include 'COMMON.INTERACT'
+      include 'COMMON.CONTROL'
+      include 'COMMON.IOUNITS'
       dimension ggg(3)
       ehpb=0.0D0
 cd    print *,'edis: nhpb=',nhpb,' fbr=',fbr
 cd    print *,'link_start=',link_start,' link_end=',link_end
+C      write(iout,*) link_end, "link_end"
       if (link_end.eq.0) return
       do i=link_start,link_end
 C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
@@ -2969,28 +2972,96 @@ C 24/11/03 AL: SS bridges handled separately because of introducing a specific
 C    distance and angle dependent SS bond potential.
 C        if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. 
 C     & iabs(itype(jjj)).eq.1) then
-
+C       write(iout,*) constr_dist,"const"
        if (.not.dyn_ss .and. i.le.nss) then
          if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
      & iabs(itype(jjj)).eq.1) then
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
-           endif
-        else
-C Calculate the distance between the two points and its difference from the
-C target distance.
-        dd=dist(ii,jj)
-        rdis=dd-dhpb(i)
+           endif !ii.gt.neres
+        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
+C            ehpb=ehpb+fordepth(i)**4.0d0
+C     &          *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+            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
+C          write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
+C     &    ehpb,fordepth(i),dd
+C            write(iout,*) ehpb,"atu?"
+C            ehpb,"tu?"
+C            fac=fordepth(i)**4.0d0
+C     &          *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(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
+c            write (iout,*) "beta nmr",
+c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+          else
+            dd=dist(ii,jj)
+            rdis=dd-dhpb(i)
+C Get the force constant corresponding to this distance.
+            waga=forcon(i)
+C Calculate the contribution to energy.
+            ehpb=ehpb+waga*rdis*rdis
+c            write (iout,*) "beta reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+            fac=waga*rdis/dd
+          endif !end dhpb1(i).gt.0
+          endif !end const_dist=11
+          do j=1,3
+            ggg(j)=fac*(c(j,jj)-c(j,ii))
+          enddo
+          do j=1,3
+            ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+            ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+          enddo
+          do k=1,3
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+          enddo
+        else !ii.gt.nres
+C          write(iout,*) "before"
+          dd=dist(ii,jj)
+C          write(iout,*) "after",dd
+          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
+C            ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i))
+C            fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd
+C            print *,ehpb,"tu?"
+C            write(iout,*) ehpb,"btu?",
+C     & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i)
+C          write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
+C     &    ehpb,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
+c            write (iout,*) "alph nmr",
+c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+          else
+            rdis=dd-dhpb(i)
 C Get the force constant corresponding to this distance.
-        waga=forcon(i)
+            waga=forcon(i)
 C Calculate the contribution to energy.
-        ehpb=ehpb+waga*rdis*rdis
+            ehpb=ehpb+waga*rdis*rdis
+c            write (iout,*) "alpha reg",dd,waga*rdis*rdis
 C
 C Evaluate gradient.
 C
-        fac=waga*rdis/dd
-cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
-cd   &   ' waga=',waga,' fac=',fac
+            fac=waga*rdis/dd
+          endif
+          endif
+
         do j=1,3
           ggg(j)=fac*(c(j,jj)-c(j,ii))
         enddo
@@ -3010,7 +3081,7 @@ C Cartesian gradient in the SC vectors (ghpbx).
         enddo
         endif
       enddo
-      ehpb=0.5D0*ehpb
+      if (constr_dist.ne.11) ehpb=0.5D0*ehpb
       return
       end
 C--------------------------------------------------------------------------
@@ -4659,7 +4730,7 @@ c   3 = SC...Ca...Ca...SCi
            esccor=esccor+v1ij*cosphi+v2ij*sinphi
            gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
          enddo
-      write (iout,*)"EBACK_SC_COR",esccor,i
+C      write (iout,*)"EBACK_SC_COR",esccor,i
 c      write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp,
 c     & nterm_sccor(isccori,isccori1),isccori,isccori1
 c        gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci