Sept 20/2012 by Adam: commit changes mostly in the MINIM version
[unres.git] / source / unres / src_MD / energy_p_new_barrier.F
index 3979f64..9edadf8 100644 (file)
@@ -471,7 +471,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'mpif.h'
 #endif
       double precision gradbufc(3,maxres),gradbufx(3,maxres),
       include 'mpif.h'
 #endif
       double precision gradbufc(3,maxres),gradbufx(3,maxres),
-     &  glocbuf(4*maxres),gradbufc_sum(3,maxres)
+     &  glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
@@ -774,6 +774,21 @@ c      enddo
         do i=1,4*nres
           glocbuf(i)=gloc(i,icg)
         enddo
         do i=1,4*nres
           glocbuf(i)=gloc(i,icg)
         enddo
+#define DEBUG
+#ifdef DEBUG
+      write (iout,*) "gloc_sc before reduce"
+      do i=1,nres
+       do j=1,3
+        write (iout,*) i,j,gloc_sc(j,i,icg)
+       enddo
+      enddo
+#endif
+#undef DEBUG
+        do i=1,nres
+         do j=1,3
+          gloc_scbuf(j,i)=gloc_sc(j,i,icg)
+         enddo
+        enddo
         time00=MPI_Wtime()
         call MPI_Barrier(FG_COMM,IERR)
         time_barrier_g=time_barrier_g+MPI_Wtime()-time00
         time00=MPI_Wtime()
         call MPI_Barrier(FG_COMM,IERR)
         time_barrier_g=time_barrier_g+MPI_Wtime()-time00
@@ -784,7 +799,19 @@ c      enddo
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
+        call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres,
+     &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         time_reduce=time_reduce+MPI_Wtime()-time00
         time_reduce=time_reduce+MPI_Wtime()-time00
+#define DEBUG
+#ifdef DEBUG
+      write (iout,*) "gloc_sc after reduce"
+      do i=1,nres
+       do j=1,3
+        write (iout,*) i,j,gloc_sc(j,i,icg)
+       enddo
+      enddo
+#endif
+#undef DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc after reduce"
       do i=1,4*nres
 #ifdef DEBUG
       write (iout,*) "gloc after reduce"
       do i=1,4*nres
@@ -1038,7 +1065,7 @@ C------------------------------------------------------------------------
      & 'ESC=   ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/
      & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/
      & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/
      & 'ESC=   ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/
      & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/
      & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/
-     & 'EHBP=  ',1pE16.6,' WEIGHT=',1pD16.6,
+     & 'EHPB=  ',1pE16.6,' WEIGHT=',1pD16.6,
      & ' (SS bridges & dist. cnstr.)'/
      & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
      & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
      & ' (SS bridges & dist. cnstr.)'/
      & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
      & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
@@ -4239,49 +4266,90 @@ C iii and jjj point to the residues for which the distance is assigned.
           iii=ii
           jjj=jj
         endif
           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
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
 cd          write (iout,*) "eij",eij
 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.
         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.
 C Get the force constant corresponding to this distance.
-        waga=forcon(i)
+            waga=forcon(i)
 C Calculate the contribution to energy.
 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
 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
 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).
 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
           do j=1,3
             ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
             ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
           enddo
-        endif
+          endif
 cgrad        do j=iii,jjj-1
 cgrad          do k=1,3
 cgrad            ghpbc(k,j)=ghpbc(k,j)+ggg(k)
 cgrad          enddo
 cgrad        enddo
 cgrad        do j=iii,jjj-1
 cgrad          do k=1,3
 cgrad            ghpbc(k,j)=ghpbc(k,j)+ggg(k)
 cgrad          enddo
 cgrad        enddo
-        do k=1,3
-          ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
-          ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
-        enddo
+          do k=1,3
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+          enddo
         endif
       enddo
       ehpb=0.5D0*ehpb
         endif
       enddo
       ehpb=0.5D0*ehpb
@@ -5869,7 +5937,7 @@ C Set lprn=.true. for debugging
 c      lprn=.true.
 c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
       esccor=0.0D0
 c      lprn=.true.
 c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
       esccor=0.0D0
-      do i=iphi_start-1,iphi_end+1
+      do i=itau_start,itau_end
         esccor_ii=0.0D0
         isccori=isccortyp(itype(i-2))
         isccori1=isccortyp(itype(i-1))
         esccor_ii=0.0D0
         isccori=isccortyp(itype(i-2))
         isccori1=isccortyp(itype(i-1))
@@ -5880,13 +5948,14 @@ c(see comment below)
 cc Omicron is flat angle depending on the value of first digit 
 c(see comment below)
 
 cc Omicron is flat angle depending on the value of first digit 
 c(see comment below)
 
-        gloci=0.0D0
-        do intertyp=1,1 !intertyp
+        
+        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
 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...SC
+c   3 = SC...Ca...Ca...SCi
+        gloci=0.0D0
         if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
      &      (itype(i-1).eq.10).or.(itype(i-2).eq.21).or.
      &      (itype(i-1).eq.21)))
         if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
      &      (itype(i-1).eq.10).or.(itype(i-2).eq.21).or.
      &      (itype(i-1).eq.21)))
@@ -5894,8 +5963,9 @@ c   3 = SC...Ca...Ca...SC
      &     .or.(itype(i-2).eq.21)))
      &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
      &      (itype(i-1).eq.21)))) cycle  
      &     .or.(itype(i-2).eq.21)))
      &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
      &      (itype(i-1).eq.21)))) cycle  
-        if ((intertyp.eq.2).and.(i.le.iphi_start-1)) cycle
-        if ((intertyp.eq.1).and.(i.ge.iphi_end+1)) cycle
+        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle
+        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21))
+     & cycle
         do j=1,nterm_sccor(isccori,isccori1)
           v1ij=v1sccor(j,intertyp,isccori,isccori1)
           v2ij=v2sccor(j,intertyp,isccori,isccori1)
         do j=1,nterm_sccor(isccori,isccori1)
           v1ij=v1sccor(j,intertyp,isccori,isccori1)
           v2ij=v2sccor(j,intertyp,isccori,isccori1)
@@ -5905,7 +5975,7 @@ c   3 = SC...Ca...Ca...SC
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
-c        write (iout,*) "WTF",intertyp,i,itype(i),
+c        write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi
 c     &gloc_sc(intertyp,i-3,icg)
         if (lprn)
      &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
 c     &gloc_sc(intertyp,i-3,icg)
         if (lprn)
      &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
@@ -8024,9 +8094,9 @@ C
 C      Parallel       Antiparallel
 C                                             
 C          o             o         
 C      Parallel       Antiparallel
 C                                             
 C          o             o         
-C         /l\           /j\       
-C        /   \         /   \      
-C       /| o |         | o |\     
+C         /l\           /j\
+C        /   \         /   \
+C       /| o |         | o |\
 C     \ j|/k\|  /   \  |/k\|l /   
 C      \ /   \ /     \ /   \ /    
 C       o     o       o     o                
 C     \ j|/k\|  /   \  |/k\|l /   
 C      \ /   \ /     \ /   \ /    
 C       o     o       o     o                
@@ -8125,18 +8195,18 @@ c----------------------------------------------------------------------------
       logical lprn
       common /kutas/ lprn
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       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, 
 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, 
@@ -8307,18 +8377,18 @@ c----------------------------------------------------------------------------
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       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 
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
@@ -8424,18 +8494,18 @@ c----------------------------------------------------------------------------
      & auxvec1(2),auxmat1(2,2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      & 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 
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective