+++ /dev/null
- 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
- 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
- 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
-
projects / unres.git / blobdiff