X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Funres%2Fsrc_CSA_DiL%2Felecont.f;fp=source%2Funres%2Fsrc_CSA_DiL%2Felecont.f;h=0000000000000000000000000000000000000000;hp=e9ed067bb8a14b5e6c74844af75a00de454c1d32;hb=de4bc5453ea46e111d936cb85e1758ed21c08fcd;hpb=b75425747e3e2b448ca5e0ef8367712e1f339124 diff --git a/source/unres/src_CSA_DiL/elecont.f b/source/unres/src_CSA_DiL/elecont.f deleted file mode 100644 index e9ed067..0000000 --- a/source/unres/src_CSA_DiL/elecont.f +++ /dev/null @@ -1,509 +0,0 @@ - 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 -