X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2FssMD.F;h=6df4b405702580beaecebd0ab39443332285a061;hb=54453b9e3cc18886558071b1bce406b21e6566cf;hp=bf1a9b33ea477dada8140e3a08d4bcf7127052c3;hpb=29858260fdd1ede4a6a430ee7883f998abaac045;p=unres.git

diff --git a/source/unres/src_MD/ssMD.F b/source/unres/src_MD/ssMD.F
index bf1a9b3..6df4b40 100644
--- a/source/unres/src_MD/ssMD.F
+++ b/source/unres/src_MD/ssMD.F
@@ -582,13 +582,14 @@ cmc      write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
         endif
       enddo
 
+#ifdef MPI
       if (nfgtasks.gt.1)then
 
         call MPI_Reduce(newnss,g_newnss,1,
      &    MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
         call MPI_Gather(newnss,1,MPI_INTEGER,
      &                  i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
-        displ(0)=0
+C        displ(0)=0
         do i=1,nfgtasks-1,1
           displ(i)=i_newnss(i-1)+displ(i-1)
         enddo
@@ -609,6 +610,7 @@ c         print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
          enddo
         endif
       endif
+#endif
 
       diff=newnss-nss
 
@@ -622,7 +624,8 @@ cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
         enddo
 #ifndef CLUST
 #ifndef WHAM
-        if (.not.found) write(iout,'(a15,f12.2,f8.1,2i5)')
+        if (.not.found.and.fg_rank.eq.0) 
+     &      write(iout,'(a15,f12.2,f8.1,2i5)')
      &       "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
 #endif
 #endif
@@ -636,7 +639,8 @@ cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
         enddo
 #ifndef CLUST
 #ifndef WHAM
-        if (.not.found) write(iout,'(a15,f12.2,f8.1,2i5)')
+        if (.not.found.and.fg_rank.eq.0) 
+     &      write(iout,'(a15,f12.2,f8.1,2i5)')
      &       "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
 #endif
 #endif
@@ -651,41 +655,6 @@ cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
       return
       end
 
-c----------------------------------------------------------------------------
-
-#ifdef WHAM
-      subroutine read_ssHist
-      implicit none
-
-c     Includes
-      include 'DIMENSIONS'
-      include "DIMENSIONS.FREE"
-      include 'COMMON.FREE'
-
-c     Local variables
-      integer i,j
-      character*80 controlcard
-
-      do i=1,dyn_nssHist
-        call card_concat(controlcard,.true.)
-        read(controlcard,*)
-     &       dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
-      enddo
-
-      return
-      end
-#endif
-
-c----------------------------------------------------------------------------
-
-
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
-C-----------------------------------------------------------------------------
 
 c$$$c-----------------------------------------------------------------------------
 c$$$
@@ -1945,3 +1914,150 @@ c$$$      return
 c$$$      end                                                               
 c$$$
 c$$$C-----------------------------------------------------------------------------
+         subroutine triple_ssbond_ene(resi,resj,resk,eij)
+      include 'DIMENSIONS'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.LOCAL'
+      include 'COMMON.INTERACT'
+      include 'COMMON.VAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CALC'
+#ifndef CLUST
+#ifndef WHAM
+      include 'COMMON.MD'
+#endif
+#endif
+
+c     External functions
+      double precision h_base
+      external h_base
+
+c     Input arguments
+      integer resi,resj,resk
+
+c     Output arguments
+      double precision eij,eij1,eij2,eij3
+
+c     Local variables
+      logical havebond
+c      integer itypi,itypj,k,l
+      double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
+      double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
+      double precision xik,yik,zik,xjk,yjk,zjk
+      double precision sig0ij,ljd,sig,fac,e1,e2
+      double precision dcosom1(3),dcosom2(3),ed
+      double precision pom1,pom2
+      double precision ljA,ljB,ljXs
+      double precision d_ljB(1:3)
+      double precision ssA,ssB,ssC,ssXs
+      double precision ssxm,ljxm,ssm,ljm
+      double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
+
+      i=resi
+      j=resj
+      k=resk
+C      write(iout,*) resi,resj,resk
+      itypi=itype(i)
+      dxi=dc_norm(1,nres+i)
+      dyi=dc_norm(2,nres+i)
+      dzi=dc_norm(3,nres+i)
+      dsci_inv=vbld_inv(i+nres)
+      xi=c(1,nres+i)
+      yi=c(2,nres+i)
+      zi=c(3,nres+i)
+
+      itypj=itype(j)
+      xj=c(1,nres+j)
+      yj=c(2,nres+j)
+      zj=c(3,nres+j)
+      
+      dxj=dc_norm(1,nres+j)
+      dyj=dc_norm(2,nres+j)
+      dzj=dc_norm(3,nres+j)
+      dscj_inv=vbld_inv(j+nres)
+      itypk=itype(k)
+      xk=c(1,nres+k)
+      yk=c(2,nres+k)
+      zk=c(3,nres+k)
+      
+      dxk=dc_norm(1,nres+k)
+      dyk=dc_norm(2,nres+k)
+      dzk=dc_norm(3,nres+k)
+      dscj_inv=vbld_inv(k+nres)
+      xij=xj-xi
+      xik=xk-xi
+      xjk=xk-xj
+      yij=yj-yi
+      yik=yk-yi
+      yjk=yk-yj
+      zij=zj-zi
+      zik=zk-zi
+      zjk=zk-zj
+      rrij=(xij*xij+yij*yij+zij*zij)
+      rij=dsqrt(rrij)  ! sc_angular needs rij to really be the inverse
+      rrik=(xik*xik+yik*yik+zik*zik)
+      rik=dsqrt(rrik)
+      rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
+      rjk=dsqrt(rrjk)
+C there are three combination of distances for each trisulfide bonds
+C The first case the ith atom is the center
+C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
+C distance y is second distance the a,b,c,d are parameters derived for
+C this problem d parameter was set as a penalty currenlty set to 1.
+      eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
+C second case jth atom is center
+      eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
+C the third case kth atom is the center
+      eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
+C      eij2=0.0
+C      eij3=0.0
+C      eij1=0.0
+      eij=eij1+eij2+eij3
+C      write(iout,*)i,j,k,eij
+C The energy penalty calculated now time for the gradient part 
+C derivative over rij
+      fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
+     &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))  
+            gg(1)=xij*fac/rij
+            gg(2)=yij*fac/rij
+            gg(3)=zij*fac/rij
+      do m=1,3
+        gvdwx(m,i)=gvdwx(m,i)-gg(m)
+        gvdwx(m,j)=gvdwx(m,j)+gg(m)
+      enddo
+      do l=1,3
+        gvdwc(l,i)=gvdwc(l,i)-gg(l)
+        gvdwc(l,j)=gvdwc(l,j)+gg(l)
+      enddo
+C now derivative over rik
+      fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
+     &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
+            gg(1)=xik*fac/rik
+            gg(2)=yik*fac/rik
+            gg(3)=zik*fac/rik
+      do m=1,3
+        gvdwx(m,i)=gvdwx(m,i)-gg(m)
+        gvdwx(m,k)=gvdwx(m,k)+gg(m)
+      enddo
+      do l=1,3
+        gvdwc(l,i)=gvdwc(l,i)-gg(l)
+        gvdwc(l,k)=gvdwc(l,k)+gg(l)
+      enddo
+C now derivative over rjk
+      fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
+     &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
+            gg(1)=xjk*fac/rjk
+            gg(2)=yjk*fac/rjk
+            gg(3)=zjk*fac/rjk
+      do m=1,3
+        gvdwx(m,j)=gvdwx(m,j)-gg(m)
+        gvdwx(m,k)=gvdwx(m,k)+gg(m)
+      enddo
+      do l=1,3
+        gvdwc(l,j)=gvdwc(l,j)-gg(l)
+        gvdwc(l,k)=gvdwc(l,k)+gg(l)
+      enddo
+      return
+      end