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