corrections
[unres.git] / source / unres / src-HCD-5D / ssMD.F
index aa938b5..26807a0 100644 (file)
@@ -84,10 +84,13 @@ ct           rij=ran_number(rmin,rmax)
 C-----------------------------------------------------------------------------
 
       subroutine dyn_ssbond_ene(resi,resj,eij)
-c      implicit none
-
-c     Includes
+      implicit none
       include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+      include 'COMMON.SETUP'
+      integer ierr
+#endif 
       include 'COMMON.SBRIDGE'
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
@@ -96,6 +99,7 @@ c     Includes
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
+      include 'COMMON.NAMES'
 #ifndef CLUST
 #ifndef WHAM
       include 'COMMON.MD'
@@ -128,6 +132,9 @@ c      integer itypi,itypj,k,l
       double precision omega,delta_inv,deltasq_inv,fac1,fac2
 c-------FIRST METHOD
       double precision xm,d_xm(1:3)
+      double precision sslipi,sslipj,ssgradlipi,ssgradlipj
+      integer ici,icj,itypi,itypj
+      double precision boxshift,sscale,sscagrad
 c-------END FIRST METHOD
 c-------SECOND METHOD
 c$$$      double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
@@ -144,125 +151,52 @@ c-------END TESTING CODE
 
       i=resi
       j=resj
-
+      ici=icys(i)
+      icj=icys(j)
+      if (ici.eq.0 .or. icj.eq.0) then
+#ifdef MPI
+        write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)') 
+     &  "Processor",me," attempt to create",
+     &  " a disulfide link between non-cysteine residues ",restyp(i),i,
+     &  restyp(j),j
+        call MPI_Abort(MPI_COMM_WORLD,ierr)
+#else
+        write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)') 
+     &  "Processor",me," attempt to create",
+     &  " a disulfide link between non-cysteine residues ",restyp(i),i,
+     &  restyp(j),j
+        stop
+#endif
+      endif
       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)
-          xi=dmod(xi,boxxsize)
-          if (xi.lt.0) xi=xi+boxxsize
-          yi=dmod(yi,boxysize)
-          if (yi.lt.0) yi=yi+boxysize
-          zi=dmod(zi,boxzsize)
-          if (zi.lt.0) zi=zi+boxzsize
+      xi=c(1,nres+i)
+      yi=c(2,nres+i)
+      zi=c(3,nres+i)
+      call to_box(xi,yi,zi)
 C define scaling factor for lipids
 
 C        if (positi.le.0) positi=positi+boxzsize
 C        print *,i
 C first for peptide groups
 c for each residue check if it is in lipid or lipid water border area
-       if ((zi.gt.bordlipbot)
-     &.and.(zi.lt.bordliptop)) then
-C the energy transfer exist
-        if (zi.lt.buflipbot) then
-C what fraction I am in
-         fracinbuf=1.0d0-
-     &        ((positi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
-         sslipi=sscalelip(fracinbuf)
-         ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
-        elseif (zi.gt.bufliptop) then
-         fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
-         sslipi=sscalelip(fracinbuf)
-         ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
-        else
-         sslipi=1.0d0
-         ssgradlipi=0.0
-        endif
-       else
-         sslipi=0.0d0
-         ssgradlipi=0.0
-       endif
+      call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
       itypj=itype(j)
-            xj=c(1,nres+j)
-            yj=c(2,nres+j)
-            zj=c(3,nres+j)
-          xj=dmod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=dmod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=dmod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-       if ((zj.gt.bordlipbot)
-     &.and.(zj.lt.bordliptop)) then
-C the energy transfer exist
-        if (zj.lt.buflipbot) then
-C what fraction I am in
-         fracinbuf=1.0d0-
-     &        ((positi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
-         sslipj=sscalelip(fracinbuf)
-         ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
-        elseif (zi.gt.bufliptop) then
-         fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
-         sslipj=sscalelip(fracinbuf)
-         ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
-        else
-         sslipj=1.0d0
-         ssgradlipj=0.0
-        endif
-       else
-         sslipj=0.0d0
-         ssgradlipj=0.0
-       endif
+      xj=c(1,nres+j)
+      yj=c(2,nres+j)
+      zj=c(3,nres+j)
+      call to_box(xj,yj,zj)
+      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
-
-      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      subchap=0
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      subchap=0
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            subchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (subchap.eq.1) then
-          xj=xj_temp-xi
-          yj=yj_temp-yi
-          zj=zj_temp-zi
-       else
-          xj=xj_safe-xi
-          yj=yj_safe-yi
-          zj=zj_safe-zi
-       endif
-
-C     xj=c(1,nres+j)-c(1,nres+i)
-C      yj=c(2,nres+j)-c(2,nres+i)
-C      zj=c(3,nres+j)-c(3,nres+i)
+      xj=boxshift(xj-xi,boxxsize)
+      yj=boxshift(yj-yi,boxysize)
+      zj=boxshift(zj-zi,boxzsize)
       dxj=dc_norm(1,nres+j)
       dyj=dc_norm(2,nres+j)
       dzj=dc_norm(3,nres+j)
@@ -522,9 +456,10 @@ c     &         "SSBOND_E_FORM",totT,t_bath,i,j
 c        endif
 #endif
 #endif
-        dyn_ssbond_ij(i,j)=eij
-      else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
-        dyn_ssbond_ij(i,j)=1.0d300
+        dyn_ssbond_ij(ici,icj)=eij
+      else if (.not.havebond .and. dyn_ssbond_ij(ici,icj).lt.1.0d300) 
+     &then
+        dyn_ssbond_ij(ici,icj)=1.0d300
 #ifndef CLUST
 #ifndef WHAM
 c        write(iout,'(a15,f12.2,f8.1,2i5)')
@@ -646,16 +581,16 @@ c     Includes
 c     Local variables
       double precision emin
       integer i,j,imin
-      integer diff,allflag(maxdim),allnss,
-     &     allihpb(maxdim),alljhpb(maxdim),
-     &     newnss,newihpb(maxdim),newjhpb(maxdim)
+      integer diff,allflag(maxss),allnss,
+     &     allihpb(maxss),alljhpb(maxss),
+     &     newnss,newihpb(maxss),newjhpb(maxss)
       logical found
       integer i_newnss(max_fg_procs),displ(0:max_fg_procs)
-      integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
+      integer g_newihpb(maxss),g_newjhpb(maxss),g_newnss
 
       allnss=0
-      do i=1,nres-1
-        do j=i+1,nres
+      do i=1,ns-1
+        do j=i+1,ns
           if (dyn_ssbond_ij(i,j).lt.1.0d300) then
             allnss=allnss+1
             allflag(allnss)=0
@@ -2037,7 +1972,7 @@ c$$$      end
 c$$$
 c$$$C-----------------------------------------------------------------------------
 c$$$C-----------------------------------------------------------------------------
-         subroutine triple_ssbond_ene(resi,resj,resk,eij)
+      subroutine triple_ssbond_ene(resi,resj,resk,eij)
       include 'DIMENSIONS'
       include 'COMMON.SBRIDGE'
       include 'COMMON.CHAIN'