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