Merge branch 'devel' into feature-ga
[unres.git] / source / wham / src / energy_p_new.F
index 052d8c3..dfef8c3 100644 (file)
@@ -2869,16 +2869,16 @@ C Evaluate bridge-strain energy and its gradient in virtual-bond and SC vectors.
 C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
-      include 'DIMENSIONS.ZSCOPT'
       include 'COMMON.SBRIDGE'
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.VAR'
       include 'COMMON.INTERACT'
+      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
+cd      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
+cd      write(iout,*)'link_start=',link_start,' 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
@@ -2893,43 +2893,85 @@ C iii and jjj point to the residues for which the distance is assigned.
           iii=ii
           jjj=jj
         endif
+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
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
+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)
-        rdis=dd-dhpb(i)
+          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)
+            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
+            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
-          ggg(j)=fac*(c(j,jj)-c(j,ii))
-        enddo
+            do j=1,3
+              ggg(j)=fac*(c(j,jj)-c(j,ii))
+            enddo
 cd      print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
 C If this is a SC-SC distance, we need to calculate the contributions to the
 C Cartesian gradient in the SC vectors (ghpbx).
-        if (iii.lt.ii) then
+          if (iii.lt.ii) then
           do j=1,3
             ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
             ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
           enddo
-        endif
-        do j=iii,jjj-1
+          endif
           do k=1,3
-            ghpbc(k,j)=ghpbc(k,j)+ggg(k)
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
           enddo
-        enddo
         endif
       enddo
       ehpb=0.5D0*ehpb
@@ -6280,18 +6322,18 @@ c--------------------------------------------------------------------------
       logical lprn
       common /kutas/ lprn
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                              
-C      Parallel       Antiparallel
-C                                             
-C          o             o         
-C         /l\           /j\       
-C        /   \         /   \      
-C       /| o |         | o |\     
-C     \ j|/k\|  /   \  |/k\|l /   
-C      \ /   \ /     \ /   \ /    
-C       o     o       o     o                
-C       i             i                     
-C
+C                                                                              C
+C      Parallel       Antiparallel                                             C
+C                                                                              C
+C          o             o                                                     C
+C         /l\           /j\                                                    C 
+C        /   \         /   \                                                   C
+C       /| o |         | o |\                                                  C
+C     \ j|/k\|  /   \  |/k\|l /                                                C
+C      \ /   \ /     \ /   \ /                                                 C
+C       o     o       o     o                                                  C
+C       i             i                                                        C
+C                                                                              C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       itk=itortyp(itype(k))
       s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i))
@@ -6387,18 +6429,18 @@ c----------------------------------------------------------------------------
       logical lprn
       common /kutas/ lprn
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                              
-C      Parallel       Antiparallel
-C                                             
-C          o             o         
-C     \   /l\           /j\   /   
-C      \ /   \         /   \ /    
-C       o| o |         | o |o     
-C     \ j|/k\|      \  |/k\|l     
-C      \ /   \       \ /   \      
-C       o             o                      
-C       i             i                     
-C
+C                                                                              C 
+C      Parallel       Antiparallel                                             C
+C                                                                              C
+C          o             o                                                     C
+C     \   /l\           /j\   /                                                C
+C      \ /   \         /   \ /                                                 C
+C       o| o |         | o |o                                                  C
+C     \ j|/k\|      \  |/k\|l                                                  C
+C      \ /   \       \ /   \                                                   C
+C       o             o                                                        C
+C       i             i                                                        C
+C                                                                              C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 cd      write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l
 C AL 7/4/01 s1 would occur in the sixth-order moment, 
@@ -6571,18 +6613,18 @@ c----------------------------------------------------------------------------
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                              
-C      Parallel       Antiparallel
-C                                             
-C          o             o         
-C         /l\   /   \   /j\       
-C        /   \ /     \ /   \      
-C       /| o |o       o| o |\     
-C       j|/k\|  /      |/k\|l /   
-C        /   \ /       /   \ /    
-C       /     o       /     o                
-C       i             i                     
-C
+C                                                                              C
+C      Parallel       Antiparallel                                             C
+C                                                                              C
+C          o             o                                                     C
+C         /l\   /   \   /j\                                                    C
+C        /   \ /     \ /   \                                                   C
+C       /| o |o       o| o |\                                                  C
+C       j|/k\|  /      |/k\|l /                                                C
+C        /   \ /       /   \ /                                                 C
+C       /     o       /     o                                                  C
+C       i             i                                                        C
+C                                                                              C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
@@ -6689,18 +6731,18 @@ c----------------------------------------------------------------------------
      & auxvec1(2),auxmat1(2,2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                              
-C      Parallel       Antiparallel
-C                                             
-C          o             o         
-C         /l\   /   \   /j\       
-C        /   \ /     \ /   \      
-C       /| o |o       o| o |\     
-C     \ j|/k\|      \  |/k\|l     
-C      \ /   \       \ /   \      
-C       o     \       o     \                
-C       i             i                     
-C
+C                                                                              C
+C      Parallel       Antiparallel                                             C
+C                                                                              C
+C          o             o                                                     C 
+C         /l\   /   \   /j\                                                    C
+C        /   \ /     \ /   \                                                   C
+C       /| o |o       o| o |\                                                  C
+C     \ j|/k\|      \  |/k\|l                                                  C
+C      \ /   \       \ /   \                                                   C
+C       o     \       o     \                                                  C
+C       i             i                                                        C
+C                                                                              C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective