X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Felecont.f;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Felecont.f;h=a962630aedda87c6143d5cbca2082ec0c9def691;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git

diff --git a/source/unres/src_MD-M-newcorr/elecont.f b/source/unres/src_MD-M-newcorr/elecont.f
new file mode 100644
index 0000000..a962630
--- /dev/null
+++ b/source/unres/src_MD-M-newcorr/elecont.f
@@ -0,0 +1,511 @@
+      subroutine elecont(lprint,ncont,icont)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.LOCAL'
+      include 'COMMON.FFIELD'
+      include 'COMMON.NAMES'
+      logical lprint
+      double precision elpp_6(2,2),elpp_3(2,2),ael6_(2,2),ael3_(2,2)
+      double precision app_(2,2),bpp_(2,2),rpp_(2,2)
+      integer ncont,icont(2,maxcont)
+      double precision econt(maxcont)
+*
+* Load the constants of peptide bond - peptide bond interactions.
+* Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
+* proline) - determined by averaging ECEPP energy.      
+*
+* as of 7/06/91.
+*
+c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
+      data rpp_    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
+      data elpp_6  /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
+      data elpp_3  / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
+      data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
+      if (lprint) write (iout,'(a)') 
+     &  "Constants of electrostatic interaction energy expression."
+      do i=1,2
+        do j=1,2
+        rri=rpp_(i,j)**6
+        app_(i,j)=epp(i,j)*rri*rri 
+        bpp_(i,j)=-2.0*epp(i,j)*rri
+        ael6_(i,j)=elpp_6(i,j)*4.2**6
+        ael3_(i,j)=elpp_3(i,j)*4.2**3
+        if (lprint)
+     &  write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),
+     &                               ael3_(i,j)
+        enddo
+      enddo
+      ncont=0
+      ees=0.0
+      evdw=0.0
+      do 1 i=nnt,nct-2
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) goto 1
+        xi=c(1,i)
+        yi=c(2,i)
+        zi=c(3,i)
+        dxi=c(1,i+1)-c(1,i)
+        dyi=c(2,i+1)-c(2,i)
+        dzi=c(3,i+1)-c(3,i)
+        xmedi=xi+0.5*dxi
+        ymedi=yi+0.5*dyi
+        zmedi=zi+0.5*dzi
+        do 4 j=i+2,nct-1
+          if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
+          ind=ind+1
+          iteli=itel(i)
+          itelj=itel(j)
+          if (j.eq.i+2 .and. itelj.eq.2) iteli=2
+          if (iteli.eq.2 .and. itelj.eq.2) goto 4
+          aaa=app_(iteli,itelj)
+          bbb=bpp_(iteli,itelj)
+          ael6_i=ael6_(iteli,itelj)
+          ael3_i=ael3_(iteli,itelj) 
+          dxj=c(1,j+1)-c(1,j)
+          dyj=c(2,j+1)-c(2,j)
+          dzj=c(3,j+1)-c(3,j)
+          xj=c(1,j)+0.5*dxj-xmedi
+          yj=c(2,j)+0.5*dyj-ymedi
+          zj=c(3,j)+0.5*dzj-zmedi
+          rrmij=1.0/(xj*xj+yj*yj+zj*zj)
+          rmij=sqrt(rrmij)
+          r3ij=rrmij*rmij
+          r6ij=r3ij*r3ij  
+          vrmij=vblinv*rmij
+          cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2      
+          cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
+          cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
+          fac=cosa-3.0*cosb*cosg
+          ev1=aaa*r6ij*r6ij
+          ev2=bbb*r6ij
+          fac3=ael6_i*r6ij
+          fac4=ael3_i*r3ij
+          evdwij=ev1+ev2
+          el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
+          el2=fac4*fac       
+          eesij=el1+el2
+          if (j.gt.i+2 .and. eesij.le.elcutoff .or.
+     &        j.eq.i+2 .and. eesij.le.elecutoff_14) then
+             ncont=ncont+1
+             icont(1,ncont)=i
+             icont(2,ncont)=j
+	     econt(ncont)=eesij
+          endif
+          ees=ees+eesij
+          evdw=evdw+evdwij
+    4   continue
+    1 continue
+      if (lprint) then
+        write (iout,*) 'Total average electrostatic energy: ',ees
+        write (iout,*) 'VDW energy between peptide-group centers: ',evdw
+        write (iout,*)
+        write (iout,*) 'Electrostatic contacts before pruning: '
+        do i=1,ncont
+          i1=icont(1,i)
+          i2=icont(2,i)
+          it1=itype(i1)
+          it2=itype(i2)
+          write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
+     &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
+        enddo
+      endif
+c For given residues keep only the contacts with the greatest energy.
+      i=0
+      do while (i.lt.ncont)
+        i=i+1
+        ene=econt(i)
+        ic1=icont(1,i)
+        ic2=icont(2,i)
+        j=i
+        do while (j.lt.ncont)
+          j=j+1
+          if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
+     &        ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
+c            write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
+c     &       " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
+            if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
+              if (ic1.eq.icont(1,j)) then
+                do k=1,ncont
+                  if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
+     &               .and. iabs(icont(1,k)-ic1).le.2 .and. 
+     &               econt(k).lt.econt(j) ) goto 21 
+                enddo
+              else if (ic2.eq.icont(2,j) ) then
+                do k=1,ncont
+                  if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
+     &               .and. iabs(icont(2,k)-ic2).le.2 .and. 
+     &               econt(k).lt.econt(j) ) goto 21 
+                enddo
+              endif
+c Remove ith contact
+              do k=i+1,ncont
+                icont(1,k-1)=icont(1,k)
+                icont(2,k-1)=icont(2,k)
+                econt(k-1)=econt(k) 
+              enddo
+              i=i-1
+              ncont=ncont-1
+c              write (iout,*) "ncont",ncont
+c              do k=1,ncont
+c                write (iout,*) icont(1,k),icont(2,k)
+c              enddo
+              goto 20
+            else if (econt(j).gt.ene .and. ic2.ne.ic1+2) 
+     &      then
+              if (ic1.eq.icont(1,j)) then
+                do k=1,ncont
+                  if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
+     &               .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. 
+     &               econt(k).lt.econt(i) ) goto 21 
+                enddo
+              else if (ic2.eq.icont(2,j) ) then
+                do k=1,ncont
+                  if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
+     &               .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. 
+     &               econt(k).lt.econt(i) ) goto 21 
+                enddo
+              endif
+c Remove jth contact
+              do k=j+1,ncont
+                icont(1,k-1)=icont(1,k)
+                icont(2,k-1)=icont(2,k)
+                econt(k-1)=econt(k) 
+              enddo
+              ncont=ncont-1
+c              write (iout,*) "ncont",ncont
+c              do k=1,ncont
+c                write (iout,*) icont(1,k),icont(2,k)
+c              enddo
+              j=j-1
+            endif   
+          endif
+   21     continue
+        enddo
+   20   continue
+      enddo
+      if (lprint) then
+        write (iout,*)
+        write (iout,*) 'Electrostatic contacts after pruning: '
+        do i=1,ncont
+          i1=icont(1,i)
+          i2=icont(2,i)
+          it1=itype(i1)
+          it2=itype(i2)
+          write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
+     &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
+        enddo
+      endif
+      return
+      end
+c--------------------------------------------
+      subroutine secondary2(lprint)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.DISTFIT'
+      include 'COMMON.VAR'
+      include 'COMMON.GEO'
+      include 'COMMON.CONTROL'
+      integer ncont,icont(2,maxcont),isec(maxres,4),nsec(maxres)
+      logical lprint,not_done,freeres
+      double precision p1,p2
+      external freeres
+
+      if(.not.dccart) call chainbuild
+cd      call write_pdb(99,'sec structure',0d0)
+      ncont=0
+      nbfrag=0
+      nhfrag=0
+      do i=1,nres
+        isec(i,1)=0
+        isec(i,2)=0
+        nsec(i)=0
+      enddo
+
+      call elecont(lprint,ncont,icont)
+
+c finding parallel beta
+cd      write (iout,*) '------- looking for parallel beta -----------'
+      nbeta=0
+      nstrand=0
+      do i=1,ncont
+        i1=icont(1,i)
+        j1=icont(2,i)
+        if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
+          ii1=i1
+          jj1=j1
+cd          write (iout,*) i1,j1
+          not_done=.true.
+          do while (not_done)
+           i1=i1+1
+           j1=j1+1
+            do j=1,ncont
+              if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and.
+     &             freeres(i1,j1,nsec,isec)) goto 5
+            enddo
+            not_done=.false.
+  5         continue
+cd            write (iout,*) i1,j1,not_done
+          enddo
+          j1=j1-1
+          i1=i1-1
+          if (i1-ii1.gt.1) then
+            ii1=max0(ii1-1,1)
+            jj1=max0(jj1-1,1)
+            nbeta=nbeta+1
+            if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
+     &               nbeta,ii1,i1,jj1,j1
+
+            nbfrag=nbfrag+1
+            bfrag(1,nbfrag)=ii1+1
+            bfrag(2,nbfrag)=i1+1
+            bfrag(3,nbfrag)=jj1+1
+            bfrag(4,nbfrag)=min0(j1+1,nres) 
+
+            do ij=ii1,i1
+             nsec(ij)=nsec(ij)+1
+             isec(ij,nsec(ij))=nbeta
+            enddo
+            do ij=jj1,j1
+             nsec(ij)=nsec(ij)+1
+             isec(ij,nsec(ij))=nbeta
+            enddo
+
+           if(lprint) then 
+            nstrand=nstrand+1
+            if (nbeta.le.9) then
+              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",ii1-1,"..",i1-1,"'"
+            else
+              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",ii1-1,"..",i1-1,"'"
+            endif
+            nstrand=nstrand+1
+            if (nbeta.le.9) then
+              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",jj1-1,"..",j1-1,"'"
+            else
+              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",jj1-1,"..",j1-1,"'"
+            endif
+              write(12,'(a8,4i4)')
+     &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
+           endif
+          endif
+        endif
+      enddo
+
+c finding alpha or 310 helix
+
+      nhelix=0
+      do i=1,ncont
+        i1=icont(1,i)
+        j1=icont(2,i)
+        p1=phi(i1+2)*rad2deg
+        p2=0.0
+        if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
+
+
+        if (j1.eq.i1+3 .and. 
+     &       ((p1.ge.10.and.p1.le.80).or.i1.le.2).and.
+     &       ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
+cd          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
+co          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
+          ii1=i1
+          jj1=j1
+          if (nsec(ii1).eq.0) then 
+            not_done=.true.
+          else
+            not_done=.false.
+          endif
+          do while (not_done)
+            i1=i1+1
+            j1=j1+1
+            do j=1,ncont
+              if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
+            enddo
+            not_done=.false.
+  10        continue
+            p1=phi(i1+2)*rad2deg
+            p2=phi(j1+2)*rad2deg
+            if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) 
+     &                              not_done=.false.
+cd           write (iout,*) i1,j1,not_done,p1,p2
+          enddo
+          j1=j1+1
+          if (j1-ii1.gt.5) then
+            nhelix=nhelix+1
+cd            write (iout,*)'helix',nhelix,ii1,j1
+
+            nhfrag=nhfrag+1
+            hfrag(1,nhfrag)=ii1
+            hfrag(2,nhfrag)=j1
+
+            do ij=ii1,j1
+             nsec(ij)=-1
+            enddo
+           if (lprint) then
+            write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
+            if (nhelix.le.9) then
+              write(12,'(a17,i1,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'helix",nhelix,
+     &          "' 'num = ",ii1-1,"..",j1-2,"'"
+            else
+              write(12,'(a17,i2,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'helix",nhelix,
+     &          "' 'num = ",ii1-1,"..",j1-2,"'"
+            endif
+           endif
+          endif
+        endif
+      enddo
+       
+      if (nhelix.gt.0.and.lprint) then
+        write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
+        do i=2,nhelix
+         if (nhelix.le.9) then
+          write(12,'(a8,i1,$)') " | helix",i
+         else
+          write(12,'(a8,i2,$)') " | helix",i
+         endif
+        enddo
+        write(12,'(a1)') "'"
+      endif
+
+
+c finding antiparallel beta
+cd      write (iout,*) '--------- looking for antiparallel beta ---------'
+
+      do i=1,ncont
+        i1=icont(1,i)
+        j1=icont(2,i)
+        if (freeres(i1,j1,nsec,isec)) then
+          ii1=i1
+          jj1=j1
+cd          write (iout,*) i1,j1
+
+          not_done=.true.
+          do while (not_done)
+           i1=i1+1
+           j1=j1-1
+            do j=1,ncont
+              if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
+     &             freeres(i1,j1,nsec,isec)) goto 6
+            enddo
+            not_done=.false.
+  6         continue
+cd            write (iout,*) i1,j1,not_done
+          enddo
+          i1=i1-1
+          j1=j1+1
+          if (i1-ii1.gt.1) then
+
+            nbfrag=nbfrag+1
+            bfrag(1,nbfrag)=ii1
+            bfrag(2,nbfrag)=min0(i1+1,nres)
+            bfrag(3,nbfrag)=min0(jj1+1,nres)
+            bfrag(4,nbfrag)=j1
+
+            nbeta=nbeta+1
+            iii1=max0(ii1-1,1)
+            do ij=iii1,i1
+             nsec(ij)=nsec(ij)+1
+             if (nsec(ij).le.2) then
+              isec(ij,nsec(ij))=nbeta
+             endif
+            enddo
+            jjj1=max0(j1-1,1)  
+            do ij=jjj1,jj1
+             nsec(ij)=nsec(ij)+1
+             if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
+              isec(ij,nsec(ij))=nbeta
+             endif
+            enddo
+
+
+           if (lprint) then
+            write (iout,'(a,i3,4i4)')'antiparallel beta',
+     &                   nbeta,ii1-1,i1,jj1,j1-1
+            nstrand=nstrand+1
+            if (nstrand.le.9) then
+              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",ii1-2,"..",i1-1,"'"
+            else
+              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",ii1-2,"..",i1-1,"'"
+            endif
+            nstrand=nstrand+1
+            if (nstrand.le.9) then
+              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",j1-2,"..",jj1-1,"'"
+            else
+              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",j1-2,"..",jj1-1,"'"
+            endif
+              write(12,'(a8,4i4)')
+     &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
+           endif
+          endif
+        endif
+      enddo
+
+      if (nstrand.gt.0.and.lprint) then
+        write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
+        do i=2,nstrand
+         if (i.le.9) then
+          write(12,'(a9,i1,$)') " | strand",i
+         else
+          write(12,'(a9,i2,$)') " | strand",i
+         endif
+        enddo
+        write(12,'(a1)') "'"
+      endif
+
+       
+
+      if (lprint) then
+       write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
+       write(12,'(a20)') "XMacStand ribbon.mac"
+         
+        
+       write(iout,*) 'UNRES seq:'
+       do j=1,nbfrag
+        write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
+       enddo
+
+       do j=1,nhfrag
+        write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
+       enddo
+      endif       
+
+      return
+      end
+c-------------------------------------------------
+      logical function freeres(i,j,nsec,isec)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      integer isec(maxres,4),nsec(maxres)
+      freeres=.false.
+
+      if (nsec(i).lt.0.or.nsec(j).lt.0) return
+      if (nsec(i).gt.1.or.nsec(j).gt.1) return
+      do k=1,nsec(i)
+        do l=1,nsec(j)
+          if (isec(i,k).eq.isec(j,l)) return
+        enddo
+      enddo
+      freeres=.true.
+      return
+      end
+