restraints from contact prediction for CSA version
[unres.git] / source / unres / src_CSA / energy_p_new_barrier.F
index 1a9044f..2121254 100644 (file)
@@ -125,13 +125,13 @@ C
   107 continue
       
 C     JUYONG for dfa test!
-      call edfad(edfadis)
+      if (wdfa_dist.gt.0) call edfad(edfadis)
 c      print*, 'edfad is finished!', edfadis
-      call edfat(edfator)
+      if (wdfa_tor.gt.0) call edfat(edfator)
 c      print*, 'edfat is finished!', edfator
-      call edfan(edfanei)
+      if (wdfa_nei.gt.0) call edfan(edfanei)
 c      print*, 'edfan is finished!', edfanei
-      call edfab(edfabet)
+      if (wdfa_beta.gt.0) call edfab(edfabet)
 c      print*, 'edfab is finished!', edfabet
 C      stop
 C     JUYONG
@@ -1082,6 +1082,9 @@ C
       include 'COMMON.NAMES'
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       dimension gg(3)
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
@@ -1991,6 +1994,9 @@ C
       include 'COMMON.NAMES'
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       dimension gg(3)
 cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
       evdw=0.0D0
@@ -2063,6 +2069,9 @@ C
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -2421,6 +2430,9 @@ C--------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -2875,6 +2887,9 @@ C
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -3037,6 +3052,9 @@ C-------------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -3680,6 +3698,9 @@ C Third- and fourth-order contributions from turns
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -3780,6 +3801,9 @@ C Third- and fourth-order contributions from turns
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
@@ -4215,26 +4239,72 @@ C iii and jjj point to the residues for which the distance is assigned.
           iii=ii
           jjj=jj
         endif
-cd        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
+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 (i.le.nss) then
+         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
+          dd=dist(ii,jj)
+          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  
+          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
 C Calculate the distance between the two points and its difference from the
 C target distance.
         dd=dist(ii,jj)
+          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)
 C Calculate the contribution to energy.
         ehpb=ehpb+waga*rdis*rdis
+c            write (iout,*) "alpha reg",dd,waga*rdis*rdis
 C
 C Evaluate gradient.
 C
         fac=waga*rdis/dd
+          endif
 cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
 cd   &   ' waga=',waga,' fac=',fac
         do j=1,3
@@ -5842,28 +5912,52 @@ c        amino-acid residues.
 C Set lprn=.true. for debugging
       lprn=.false.
 c      lprn=.true.
-c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
+c      write (iout,*) "EBACK_SC_COR",itau_start,itau_end
       esccor=0.0D0
-      do i=iphi_start,iphi_end
+      do i=itau_start,itau_end
         esccor_ii=0.0D0
-        itori=itype(i-2)
-        itori1=itype(i-1)
+        isccori=isccortyp(itype(i-2))
+        isccori1=isccortyp(itype(i-1))
+c      write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
         phii=phi(i)
+        do intertyp=1,3 !intertyp
+cc Added 09 May 2012 (Adasko)
+cc  Intertyp means interaction type of backbone mainchain correlation: 
+c   1 = SC...Ca...Ca...Ca
+c   2 = Ca...Ca...Ca...SC
+c   3 = SC...Ca...Ca...SCi
         gloci=0.0D0
-        do j=1,nterm_sccor
-          v1ij=v1sccor(j,itori,itori1)
-          v2ij=v2sccor(j,itori,itori1)
-          cosphi=dcos(j*phii)
-          sinphi=dsin(j*phii)
+        if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
+     &      (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or.
+     &      (itype(i-1).eq.ntyp1)))
+     &    .or. ((intertyp.eq.1).and.((itype(i-2).eq.10)
+     &     .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)
+     &     .or.(itype(i).eq.ntyp1)))
+     &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
+     &      (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &      (itype(i-3).eq.ntyp1)))) cycle
+        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle
+        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1))
+     & cycle
+       do j=1,nterm_sccor(isccori,isccori1)
+          v1ij=v1sccor(j,intertyp,isccori,isccori1)
+          v2ij=v2sccor(j,intertyp,isccori,isccori1)
+          cosphi=dcos(j*tauangle(intertyp,i))
+          sinphi=dsin(j*tauangle(intertyp,i))
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
+c      write (iout,*) "EBACK_SC_COR",i,esccor,intertyp
+        gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
         if (lprn)
      &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
-     &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
-     &  (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6)
-        gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
+     &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,isccori,isccori1,
+     &  (v1sccor(j,intertyp,isccori,isccori1),j=1,6)
+     & ,(v2sccor(j,intertyp,isccori,isccori1),j=1,6)
+C        gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
+       enddo !intertyp
       enddo
+
       return
       end
 c----------------------------------------------------------------------------
@@ -5878,6 +5972,9 @@ C contribution equal to sqrt(eps(i,j)*eps(i+1,j+1)) is added.
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       double precision gx(3),gx1(3)
       logical lprn
 
@@ -5932,6 +6029,9 @@ c------------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       double precision gx(3),gx1(3)
       logical lprn
       lprn=.false.
@@ -5985,6 +6085,9 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.CONTROL'
       include 'COMMON.LOCAL'
       double precision gx(3),gx1(3),time00
@@ -6277,6 +6380,9 @@ c------------------------------------------------------------------------------
       parameter (max_cont=maxconts)
       parameter (max_dim=26)
       include "COMMON.CONTACTS"
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       double precision zapas(max_dim,maxconts,max_fg_procs),
      &  zapas_recv(max_dim,maxconts,max_fg_procs)
       common /przechowalnia/ zapas
@@ -6348,6 +6454,9 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding
       include 'COMMON.LOCAL'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.CHAIN'
       include 'COMMON.CONTROL'
       double precision gx(3),gx1(3)
@@ -6718,6 +6827,9 @@ c------------------------------------------------------------------------------
       parameter (max_cont=maxconts)
       parameter (max_dim=70)
       include "COMMON.CONTACTS"
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       double precision zapas(max_dim,maxconts,max_fg_procs),
      &  zapas_recv(max_dim,maxconts,max_fg_procs)
       common /przechowalnia/ zapas
@@ -6771,6 +6883,9 @@ c------------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       double precision gx(3),gx1(3)
       logical lprn
       lprn=.false.
@@ -6859,6 +6974,9 @@ C---------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -6924,6 +7042,9 @@ C
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -7302,6 +7423,9 @@ C---------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -7414,6 +7538,9 @@ C---------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -7818,6 +7945,9 @@ c--------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -7958,6 +8088,9 @@ c--------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -8062,12 +8195,15 @@ c----------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
       logical swap
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2),
-     & auxvec1(2),auxvec2(1),auxmat1(2,2)
+     & auxvec1(2),auxvec2(2),auxmat1(2,2)
       logical lprn
       common /kutas/ lprn
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
@@ -8247,6 +8383,9 @@ c----------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -8362,6 +8501,9 @@ c----------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
@@ -8606,6 +8748,9 @@ c----------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+#ifdef MOMENT
+      include 'COMMON.CONTACTS.MOMENT'
+#endif 
       include 'COMMON.TORSION'
       include 'COMMON.VAR'
       include 'COMMON.GEO'