edis ctest multichain
[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F
index 8932003..37f648f 100644 (file)
@@ -5192,7 +5192,7 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
-      dimension ggg(3)
+      dimension ggg(3),ggg_peak(3,20)
       ehpb=0.0D0
       do i=1,3
        ggg(i)=0.0d0
@@ -5200,8 +5200,69 @@ C
 C      write (iout,*) ,"link_end",link_end,constr_dist
 cd      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
 c      write(iout,*)'link_start=',link_start,' link_end=',link_end,
-c     &  " constr_dist",constr_dist
-      if (link_end.eq.0) return
+c     &  " constr_dist",constr_dist," link_start_peak",link_start_peak,
+c     &  " link_end_peak",link_end_peak
+      if (link_end.eq.0.and.link_end_peak.eq.0) return
+      if (link_end_peak.ne.0) then
+      do i=link_start_peak,link_end_peak
+        ehpb_peak=0.0d0
+c        print *,"i",i," link_end_peak",link_end_peak," ipeak",
+c     &   ipeak(1,i),ipeak(2,i)
+        do ip=ipeak(1,i),ipeak(2,i)
+          ii=ihpb_peak(ip)
+          jj=jhpb_peak(ip)
+          dd=dist(ii,jj)
+          iip=ip-ipeak(1,i)+1
+C iii and jjj point to the residues for which the distance is assigned.
+          if (ii.gt.nres) then
+            iii=ii-nres
+            jjj=jj-nres 
+          else
+            iii=ii
+            jjj=jj
+          endif
+          aux=rlornmr1(dd,dhpb_peak(ip),dhpb1_peak(ip),forcon_peak(ip))
+          aux=dexp(-scal_peak*aux)
+          ehpb_peak=ehpb_peak+aux
+          fac=rlornmr1prim(dd,dhpb_peak(ip),dhpb1_peak(ip),
+     &      forcon_peak(ip))*aux/dd
+          do j=1,3
+            ggg_peak(j,iip)=fac*(c(j,jj)-c(j,ii))
+          enddo
+          if (energy_dec) write (iout,'(a6,3i5,6f10.3,i5)')
+     &      "edisL",i,ii,jj,dd,dhpb_peak(ip),dhpb1_peak(ip),
+     &      forcon_peak(ip),fordepth_peak(ip),ehpb_peak
+        enddo
+c        write (iout,*) "ehpb_peak",ehpb_peak," scal_peak",scal_peak
+        ehpb=ehpb-fordepth_peak(ipeak(1,i))*dlog(ehpb_peak)/scal_peak
+        do ip=ipeak(1,i),ipeak(2,i)
+          iip=ip-ipeak(1,i)+1
+          do j=1,3
+            ggg(j)=ggg_peak(j,iip)/ehpb_peak
+          enddo
+          ii=ihpb_peak(ip)
+          jj=jhpb_peak(ip)
+C iii and jjj point to the residues for which the distance is assigned.
+          if (ii.gt.nres) then
+            iii=ii-nres
+            jjj=jj-nres 
+          else
+            iii=ii
+            jjj=jj
+          endif
+          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 k=1,3
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+          enddo
+        enddo
+      enddo
+      endif
       do i=link_start,link_end
 C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
 C CA-CA distance used in regularization of structure.
@@ -5228,7 +5289,7 @@ C 15/02/13 CC dynamic SSbond - additional check
           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
+           ehpb=ehpb+eij
          endif
 cd          write (iout,*) "eij",eij
 cd   &   ' waga=',waga,' fac=',fac
@@ -11328,46 +11389,46 @@ c CA CA
       enddo ! i
 #ifdef MPI
       if (nfgtasks.gt.1) then 
-        call MPI_Reduce(Pcalc(1),Pcalc_(1),nsaxs,MPI_DOUBLE_PRECISION,
-     &    MPI_SUM,king,FG_COMM,IERR)
-        if (fg_rank.eq.king) then
+       call MPI_AllReduce(Pcalc(1),Pcalc_(1),nsaxs,MPI_DOUBLE_PRECISION,
+     &    MPI_SUM,FG_COMM,IERR)
+c        if (fg_rank.eq.king) then
           do k=1,nsaxs
             Pcalc(k) = Pcalc_(k)
           enddo
-        endif
-        call MPI_Reduce(PgradC(k,1,1),PgradC_(k,1,1),3*maxsaxs*nres,
-     &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
-        if (fg_rank.eq.king) then
-          do i=1,nres
-            do l=1,3
-              do k=1,nsaxs
-                PgradC(k,l,i) = PgradC_(k,l,i)
-              enddo
-            enddo
-          enddo
-        endif
+c        endif
+c        call MPI_AllReduce(PgradC(k,1,1),PgradC_(k,1,1),3*maxsaxs*nres,
+c     &    MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR)
+c        if (fg_rank.eq.king) then
+c          do i=1,nres
+c            do l=1,3
+c              do k=1,nsaxs
+c                PgradC(k,l,i) = PgradC_(k,l,i)
+c              enddo
+c            enddo
+c          enddo
+c        endif
 #ifdef ALLSAXS
-        call MPI_Reduce(PgradX(k,1,1),PgradX_(k,1,1),3*maxsaxs*nres,
-     &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
-        if (fg_rank.eq.king) then
-          do i=1,nres
-            do l=1,3
-              do k=1,nsaxs
-                PgradX(k,l,i) = PgradX_(k,l,i)
-              enddo
-            enddo
-          enddo
-        endif
+c        call MPI_AllReduce(PgradX(k,1,1),PgradX_(k,1,1),3*maxsaxs*nres,
+c     &    MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR)
+c        if (fg_rank.eq.king) then
+c          do i=1,nres
+c            do l=1,3
+c              do k=1,nsaxs
+c                PgradX(k,l,i) = PgradX_(k,l,i)
+c              enddo
+c            enddo
+c          enddo
+c        endif
 #endif
       endif
 #endif
-#ifdef MPI
-      if (fg_rank.eq.king) then
-#endif
       Cnorm = 0.0d0
       do k=1,nsaxs
         Cnorm = Cnorm + Pcalc(k)
       enddo
+#ifdef MPI
+      if (fg_rank.eq.king) then
+#endif
       Esaxs_constr = dlog(Cnorm)-wsaxs0
       do k=1,nsaxs
         if (Pcalc(k).gt.0.0d0) 
@@ -11379,6 +11440,11 @@ c CA CA
 #ifdef DEBUG
       write (iout,*) "Cnorm",Cnorm," Esaxs_constr",Esaxs_constr
 #endif
+#ifdef MPI
+      endif
+#endif
+      gsaxsC=0.0d0
+      gsaxsX=0.0d0
       do i=nnt,nct
         do l=1,3
           auxC=0.0d0
@@ -11403,7 +11469,7 @@ c     *     " gradX",wsaxs*(auxX - auxX1/Cnorm)
         enddo
       enddo
 #ifdef MPI
-      endif
+c      endif
 #endif
       return
       end