source/unres/src_MD-M/elecont.f

   1       subroutine elecont(lprint,ncont,icont)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.IOUNITS'
   5       include 'COMMON.CHAIN'
   6       include 'COMMON.INTERACT'
   7       include 'COMMON.LOCAL'
   8       include 'COMMON.FFIELD'
   9       include 'COMMON.NAMES'
  10       logical lprint
  11       double precision elpp_6(2,2),elpp_3(2,2),ael6_(2,2),ael3_(2,2)
  12       double precision app_(2,2),bpp_(2,2),rpp_(2,2)
  13       integer ncont,icont(2,maxcont)
  14       double precision econt(maxcont)
  15 *
  16 * Load the constants of peptide bond - peptide bond interactions.
  17 * Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
  18 * proline) - determined by averaging ECEPP energy.
  19 *
  20 * as of 7/06/91.
  21 *
  22 c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
  23       data rpp_    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
  24       data elpp_6  /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
  25       data elpp_3  / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
  26       data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
  27       if (lprint) write (iout,'(a)')
  28      &  "Constants of electrostatic interaction energy expression."
  29       do i=1,2
  30         do j=1,2
  31         rri=rpp_(i,j)**6
  32         app_(i,j)=epp(i,j)*rri*rri
  33         bpp_(i,j)=-2.0*epp(i,j)*rri
  34         ael6_(i,j)=elpp_6(i,j)*4.2**6
  35         ael3_(i,j)=elpp_3(i,j)*4.2**3
  36         if (lprint)
  37      &  write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),
  38      &                               ael3_(i,j)
  39         enddo
  40       enddo
  41       ncont=0
  42       ees=0.0
  43       evdw=0.0
  44       do 1 i=nnt,nct-2
  45         if (itype(i).eq.21 .or. itype(i+1).eq.21) goto 1
  46         xi=c(1,i)
  47         yi=c(2,i)
  48         zi=c(3,i)
  49         dxi=c(1,i+1)-c(1,i)
  50         dyi=c(2,i+1)-c(2,i)
  51         dzi=c(3,i+1)-c(3,i)
  52         xmedi=xi+0.5*dxi
  53         ymedi=yi+0.5*dyi
  54         zmedi=zi+0.5*dzi
  55         do 4 j=i+2,nct-1
  56           if (itype(j).eq.21 .or. itype(j+1).eq.21) goto 4
  57           ind=ind+1
  58           iteli=itel(i)
  59           itelj=itel(j)
  60           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
  61           if (iteli.eq.2 .and. itelj.eq.2) goto 4
  62           aaa=app_(iteli,itelj)
  63           bbb=bpp_(iteli,itelj)
  64           ael6_i=ael6_(iteli,itelj)
  65           ael3_i=ael3_(iteli,itelj)
  66           dxj=c(1,j+1)-c(1,j)
  67           dyj=c(2,j+1)-c(2,j)
  68           dzj=c(3,j+1)-c(3,j)
  69           xj=c(1,j)+0.5*dxj-xmedi
  70           yj=c(2,j)+0.5*dyj-ymedi
  71           zj=c(3,j)+0.5*dzj-zmedi
  72           rrmij=1.0/(xj*xj+yj*yj+zj*zj)
  73           rmij=sqrt(rrmij)
  74           r3ij=rrmij*rmij
  75           r6ij=r3ij*r3ij
  76           vrmij=vblinv*rmij
  77           cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
  78           cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
  79           cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
  80           fac=cosa-3.0*cosb*cosg
  81           ev1=aaa*r6ij*r6ij
  82           ev2=bbb*r6ij
  83           fac3=ael6_i*r6ij
  84           fac4=ael3_i*r3ij
  85           evdwij=ev1+ev2
  86           el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
  87           el2=fac4*fac
  88           eesij=el1+el2
  89           if (j.gt.i+2 .and. eesij.le.elcutoff .or.
  90      &        j.eq.i+2 .and. eesij.le.elecutoff_14) then
  91              ncont=ncont+1
  92              icont(1,ncont)=i
  93              icont(2,ncont)=j
  94              econt(ncont)=eesij
  95           endif
  96           ees=ees+eesij
  97           evdw=evdw+evdwij
  98     4   continue
  99     1 continue
 100       if (lprint) then
 101         write (iout,*) 'Total average electrostatic energy: ',ees
 102         write (iout,*) 'VDW energy between peptide-group centers: ',evdw
 103         write (iout,*)
 104         write (iout,*) 'Electrostatic contacts before pruning: '
 105         do i=1,ncont
 106           i1=icont(1,i)
 107           i2=icont(2,i)
 108           it1=itype(i1)
 109           it2=itype(i2)
 110           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
 111      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
 112         enddo
 113       endif
 114 c For given residues keep only the contacts with the greatest energy.
 115       i=0
 116       do while (i.lt.ncont)
 117         i=i+1
 118         ene=econt(i)
 119         ic1=icont(1,i)
 120         ic2=icont(2,i)
 121         j=i
 122         do while (j.lt.ncont)
 123           j=j+1
 124           if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
 125      &        ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
 126 c            write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
 127 c     &       " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
 128             if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
 129               if (ic1.eq.icont(1,j)) then
 130                 do k=1,ncont
 131                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
 132      &               .and. iabs(icont(1,k)-ic1).le.2 .and.
 133      &               econt(k).lt.econt(j) ) goto 21
 134                 enddo
 135               else if (ic2.eq.icont(2,j) ) then
 136                 do k=1,ncont
 137                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
 138      &               .and. iabs(icont(2,k)-ic2).le.2 .and.
 139      &               econt(k).lt.econt(j) ) goto 21
 140                 enddo
 141               endif
 142 c Remove ith contact
 143               do k=i+1,ncont
 144                 icont(1,k-1)=icont(1,k)
 145                 icont(2,k-1)=icont(2,k)
 146                 econt(k-1)=econt(k)
 147               enddo
 148               i=i-1
 149               ncont=ncont-1
 150 c              write (iout,*) "ncont",ncont
 151 c              do k=1,ncont
 152 c                write (iout,*) icont(1,k),icont(2,k)
 153 c              enddo
 154               goto 20
 155             else if (econt(j).gt.ene .and. ic2.ne.ic1+2)
 156      &      then
 157               if (ic1.eq.icont(1,j)) then
 158                 do k=1,ncont
 159                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
 160      &               .and. iabs(icont(1,k)-icont(1,j)).le.2 .and.
 161      &               econt(k).lt.econt(i) ) goto 21
 162                 enddo
 163               else if (ic2.eq.icont(2,j) ) then
 164                 do k=1,ncont
 165                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
 166      &               .and. iabs(icont(2,k)-icont(2,j)).le.2 .and.
 167      &               econt(k).lt.econt(i) ) goto 21
 168                 enddo
 169               endif
 170 c Remove jth contact
 171               do k=j+1,ncont
 172                 icont(1,k-1)=icont(1,k)
 173                 icont(2,k-1)=icont(2,k)
 174                 econt(k-1)=econt(k)
 175               enddo
 176               ncont=ncont-1
 177 c              write (iout,*) "ncont",ncont
 178 c              do k=1,ncont
 179 c                write (iout,*) icont(1,k),icont(2,k)
 180 c              enddo
 181               j=j-1
 182             endif
 183           endif
 184    21     continue
 185         enddo
 186    20   continue
 187       enddo
 188       if (lprint) then
 189         write (iout,*)
 190         write (iout,*) 'Electrostatic contacts after pruning: '
 191         do i=1,ncont
 192           i1=icont(1,i)
 193           i2=icont(2,i)
 194           it1=itype(i1)
 195           it2=itype(i2)
 196           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
 197      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
 198         enddo
 199       endif
 200       return
 201       end
 202 c--------------------------------------------
 203       subroutine secondary2(lprint)
 204       implicit real*8 (a-h,o-z)
 205       include 'DIMENSIONS'
 206       include 'COMMON.CHAIN'
 207       include 'COMMON.IOUNITS'
 208       include 'COMMON.DISTFIT'
 209       include 'COMMON.VAR'
 210       include 'COMMON.GEO'
 211       include 'COMMON.CONTROL'
 212       integer ncont,icont(2,maxcont),isec(maxres,4),nsec(maxres)
 213       logical lprint,not_done,freeres
 214       double precision p1,p2
 215       external freeres
 216
 217       if(.not.dccart) call chainbuild
 218 cd      call write_pdb(99,'sec structure',0d0)
 219       ncont=0
 220       nbfrag=0
 221       nhfrag=0
 222       do i=1,nres
 223         isec(i,1)=0
 224         isec(i,2)=0
 225         nsec(i)=0
 226       enddo
 227
 228       call elecont(lprint,ncont,icont)
 229
 230 c finding parallel beta
 231 cd      write (iout,*) '------- looking for parallel beta -----------'
 232       nbeta=0
 233       nstrand=0
 234       do i=1,ncont
 235         i1=icont(1,i)
 236         j1=icont(2,i)
 237         if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
 238           ii1=i1
 239           jj1=j1
 240 cd          write (iout,*) i1,j1
 241           not_done=.true.
 242           do while (not_done)
 243            i1=i1+1
 244            j1=j1+1
 245             do j=1,ncont
 246               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and.
 247      &             freeres(i1,j1,nsec,isec)) goto 5
 248             enddo
 249             not_done=.false.
 250   5         continue
 251 cd            write (iout,*) i1,j1,not_done
 252           enddo
 253           j1=j1-1
 254           i1=i1-1
 255           if (i1-ii1.gt.1) then
 256             ii1=max0(ii1-1,1)
 257             jj1=max0(jj1-1,1)
 258             nbeta=nbeta+1
 259             if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
 260      &               nbeta,ii1,i1,jj1,j1
 261
 262             nbfrag=nbfrag+1
 263             bfrag(1,nbfrag)=ii1+1
 264             bfrag(2,nbfrag)=i1+1
 265             bfrag(3,nbfrag)=jj1+1
 266             bfrag(4,nbfrag)=min0(j1+1,nres)
 267
 268             do ij=ii1,i1
 269              nsec(ij)=nsec(ij)+1
 270              isec(ij,nsec(ij))=nbeta
 271             enddo
 272             do ij=jj1,j1
 273              nsec(ij)=nsec(ij)+1
 274              isec(ij,nsec(ij))=nbeta
 275             enddo
 276
 277            if(lprint) then
 278             nstrand=nstrand+1
 279             if (nbeta.le.9) then
 280               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 281      &          "DefPropRes 'strand",nstrand,
 282      &          "' 'num = ",ii1-1,"..",i1-1,"'"
 283             else
 284               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 285      &          "DefPropRes 'strand",nstrand,
 286      &          "' 'num = ",ii1-1,"..",i1-1,"'"
 287             endif
 288             nstrand=nstrand+1
 289             if (nbeta.le.9) then
 290               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 291      &          "DefPropRes 'strand",nstrand,
 292      &          "' 'num = ",jj1-1,"..",j1-1,"'"
 293             else
 294               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 295      &          "DefPropRes 'strand",nstrand,
 296      &          "' 'num = ",jj1-1,"..",j1-1,"'"
 297             endif
 298               write(12,'(a8,4i4)')
 299      &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
 300            endif
 301           endif
 302         endif
 303       enddo
 304
 305 c finding alpha or 310 helix
 306
 307       nhelix=0
 308       do i=1,ncont
 309         i1=icont(1,i)
 310         j1=icont(2,i)
 311         p1=phi(i1+2)*rad2deg
 312         p2=0.0
 313         if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
 314
 315
 316         if (j1.eq.i1+3 .and.
 317      &       ((p1.ge.10.and.p1.le.80).or.i1.le.2).and.
 318      &       ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
 319 cd          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
 320 co          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
 321           ii1=i1
 322           jj1=j1
 323           if (nsec(ii1).eq.0) then
 324             not_done=.true.
 325           else
 326             not_done=.false.
 327           endif
 328           do while (not_done)
 329             i1=i1+1
 330             j1=j1+1
 331             do j=1,ncont
 332               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
 333             enddo
 334             not_done=.false.
 335   10        continue
 336             p1=phi(i1+2)*rad2deg
 337             p2=phi(j1+2)*rad2deg
 338             if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80)
 339      &                              not_done=.false.
 340 cd           write (iout,*) i1,j1,not_done,p1,p2
 341           enddo
 342           j1=j1+1
 343           if (j1-ii1.gt.5) then
 344             nhelix=nhelix+1
 345 cd            write (iout,*)'helix',nhelix,ii1,j1
 346
 347             nhfrag=nhfrag+1
 348             hfrag(1,nhfrag)=ii1
 349             hfrag(2,nhfrag)=j1
 350
 351             do ij=ii1,j1
 352              nsec(ij)=-1
 353             enddo
 354            if (lprint) then
 355             write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
 356             if (nhelix.le.9) then
 357               write(12,'(a17,i1,a9,i3,a2,i3,a1)')
 358      &          "DefPropRes 'helix",nhelix,
 359      &          "' 'num = ",ii1-1,"..",j1-2,"'"
 360             else
 361               write(12,'(a17,i2,a9,i3,a2,i3,a1)')
 362      &          "DefPropRes 'helix",nhelix,
 363      &          "' 'num = ",ii1-1,"..",j1-2,"'"
 364             endif
 365            endif
 366           endif
 367         endif
 368       enddo
 369
 370       if (nhelix.gt.0.and.lprint) then
 371         write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
 372         do i=2,nhelix
 373          if (nhelix.le.9) then
 374           write(12,'(a8,i1,$)') " | helix",i
 375          else
 376           write(12,'(a8,i2,$)') " | helix",i
 377          endif
 378         enddo
 379         write(12,'(a1)') "'"
 380       endif
 381
 382
 383 c finding antiparallel beta
 384 cd      write (iout,*) '--------- looking for antiparallel beta ---------'
 385
 386       do i=1,ncont
 387         i1=icont(1,i)
 388         j1=icont(2,i)
 389         if (freeres(i1,j1,nsec,isec)) then
 390           ii1=i1
 391           jj1=j1
 392 cd          write (iout,*) i1,j1
 393
 394           not_done=.true.
 395           do while (not_done)
 396            i1=i1+1
 397            j1=j1-1
 398             do j=1,ncont
 399               if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
 400      &             freeres(i1,j1,nsec,isec)) goto 6
 401             enddo
 402             not_done=.false.
 403   6         continue
 404 cd            write (iout,*) i1,j1,not_done
 405           enddo
 406           i1=i1-1
 407           j1=j1+1
 408           if (i1-ii1.gt.1) then
 409
 410             nbfrag=nbfrag+1
 411             bfrag(1,nbfrag)=ii1
 412             bfrag(2,nbfrag)=min0(i1+1,nres)
 413             bfrag(3,nbfrag)=min0(jj1+1,nres)
 414             bfrag(4,nbfrag)=j1
 415
 416             nbeta=nbeta+1
 417             iii1=max0(ii1-1,1)
 418             do ij=iii1,i1
 419              nsec(ij)=nsec(ij)+1
 420              if (nsec(ij).le.2) then
 421               isec(ij,nsec(ij))=nbeta
 422              endif
 423             enddo
 424             jjj1=max0(j1-1,1)
 425             do ij=jjj1,jj1
 426              nsec(ij)=nsec(ij)+1
 427              if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
 428               isec(ij,nsec(ij))=nbeta
 429              endif
 430             enddo
 431
 432
 433            if (lprint) then
 434             write (iout,'(a,i3,4i4)')'antiparallel beta',
 435      &                   nbeta,ii1-1,i1,jj1,j1-1
 436             nstrand=nstrand+1
 437             if (nstrand.le.9) then
 438               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 439      &          "DefPropRes 'strand",nstrand,
 440      &          "' 'num = ",ii1-2,"..",i1-1,"'"
 441             else
 442               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 443      &          "DefPropRes 'strand",nstrand,
 444      &          "' 'num = ",ii1-2,"..",i1-1,"'"
 445             endif
 446             nstrand=nstrand+1
 447             if (nstrand.le.9) then
 448               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 449      &          "DefPropRes 'strand",nstrand,
 450      &          "' 'num = ",j1-2,"..",jj1-1,"'"
 451             else
 452               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 453      &          "DefPropRes 'strand",nstrand,
 454      &          "' 'num = ",j1-2,"..",jj1-1,"'"
 455             endif
 456               write(12,'(a8,4i4)')
 457      &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
 458            endif
 459           endif
 460         endif
 461       enddo
 462
 463       if (nstrand.gt.0.and.lprint) then
 464         write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
 465         do i=2,nstrand
 466          if (i.le.9) then
 467           write(12,'(a9,i1,$)') " | strand",i
 468          else
 469           write(12,'(a9,i2,$)') " | strand",i
 470          endif
 471         enddo
 472         write(12,'(a1)') "'"
 473       endif
 474
 475
 476
 477       if (lprint) then
 478        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
 479        write(12,'(a20)') "XMacStand ribbon.mac"
 480
 481
 482        write(iout,*) 'UNRES seq:'
 483        do j=1,nbfrag
 484         write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
 485        enddo
 486
 487        do j=1,nhfrag
 488         write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
 489        enddo
 490       endif
 491
 492       return
 493       end
 494 c-------------------------------------------------
 495       logical function freeres(i,j,nsec,isec)
 496       implicit real*8 (a-h,o-z)
 497       include 'DIMENSIONS'
 498       integer isec(maxres,4),nsec(maxres)
 499       freeres=.false.
 500
 501       if (nsec(i).lt.0.or.nsec(j).lt.0) return
 502       if (nsec(i).gt.1.or.nsec(j).gt.1) return
 503       do k=1,nsec(i)
 504         do l=1,nsec(j)
 505           if (isec(i,k).eq.isec(j,l)) return
 506         enddo
 507       enddo
 508       freeres=.true.
 509       return
 510       end
 511