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 xmedi=mod(xmedi,boxxsize) if (xmedi.lt.0) xmedi=xmedi+boxxsize ymedi=mod(ymedi,boxysize) if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize 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 yj=c(2,j)+0.5*dyj zj=c(3,j)+0.5*dzj xj=mod(xj,boxxsize) if (xj.lt.0) xj=xj+boxxsize yj=mod(yj,boxysize) if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj zj_safe=zj isubchap=0 do xshift=-1,1 do yshift=-1,1 do zshift=-1,1 xj=xj_safe+xshift*boxxsize yj=yj_safe+yshift*boxysize zj=zj_safe+zshift*boxzsize dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 if(dist_temp.lt.dist_init) then dist_init=dist_temp xj_temp=xj yj_temp=yj zj_temp=zj isubchap=1 endif enddo enddo enddo if (isubchap.eq.1) then xj=xj_temp-xmedi yj=yj_temp-ymedi zj=zj_temp-zmedi else xj=xj_safe-xmedi yj=yj_safe-ymedi zj=zj_safe-zmedi endif rij=xj*xj+yj*yj+zj*zj sss=sscale(sqrt(rij)) sssgrad=sscagrad(sqrt(rij)) 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*sss 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 cc???? 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