1 subroutine elecont(lprint,ncont,icont)
2 implicit real*8 (a-h,o-z)
4 include 'COMMON.IOUNITS'
6 include 'COMMON.INTERACT'
8 include 'COMMON.FFIELD'
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)
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.
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."
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
37 & write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),
45 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) goto 1
55 xmedi=mod(xmedi,boxxsize)
56 if (xmedi.lt.0) xmedi=xmedi+boxxsize
57 ymedi=mod(ymedi,boxysize)
58 if (ymedi.lt.0) ymedi=ymedi+boxysize
59 zmedi=mod(zmedi,boxzsize)
60 if (zmedi.lt.0) zmedi=zmedi+boxzsize
62 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
66 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
67 if (iteli.eq.2 .and. itelj.eq.2) goto 4
70 ael6_i=ael6_(iteli,itelj)
71 ael3_i=ael3_(iteli,itelj)
79 if (xj.lt.0) xj=xj+boxxsize
81 if (yj.lt.0) yj=yj+boxysize
83 if (zj.lt.0) zj=zj+boxzsize
84 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
85 write (iout,*) "i",i,xi,yi,zi," j",j,xj,yj,xj,"dist",
94 xj=xj_safe+xshift*boxxsize
95 yj=yj_safe+yshift*boxysize
96 zj=zj_safe+zshift*boxzsize
97 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
98 if(dist_temp.lt.dist_init) then
108 if (isubchap.eq.1) then
117 rij=xj*xj+yj*yj+zj*zj
118 sss=sscale(sqrt(rij))
119 sssgrad=sscagrad(sqrt(rij))
120 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
125 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
126 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
127 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
128 fac=cosa-3.0*cosb*cosg
134 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
137 if (j.gt.i+2 .and. eesij.le.elcutoff .or.
138 & j.eq.i+2 .and. eesij.le.elecutoff_14) then
146 write (iout,*) "i"," j",j," rij",dsqrt(rij)," eesij",eesij
150 write (iout,*) 'Total average electrostatic energy: ',ees
151 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
153 write (iout,*) 'Electrostatic contacts before pruning: '
159 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
160 & i,restyp(it1),i1,restyp(it2),i2,econt(i)
163 c For given residues keep only the contacts with the greatest energy.
165 do while (i.lt.ncont)
171 do while (j.lt.ncont)
173 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
174 & ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
175 c write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
176 c & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
177 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
178 if (ic1.eq.icont(1,j)) then
180 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
181 & .and. iabs(icont(1,k)-ic1).le.2 .and.
182 & econt(k).lt.econt(j) ) goto 21
184 else if (ic2.eq.icont(2,j) ) then
186 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
187 & .and. iabs(icont(2,k)-ic2).le.2 .and.
188 & econt(k).lt.econt(j) ) goto 21
193 icont(1,k-1)=icont(1,k)
194 icont(2,k-1)=icont(2,k)
199 c write (iout,*) "ncont",ncont
201 c write (iout,*) icont(1,k),icont(2,k)
204 else if (econt(j).gt.ene .and. ic2.ne.ic1+2)
206 if (ic1.eq.icont(1,j)) then
208 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
209 & .and. iabs(icont(1,k)-icont(1,j)).le.2 .and.
210 & econt(k).lt.econt(i) ) goto 21
212 else if (ic2.eq.icont(2,j) ) then
214 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
215 & .and. iabs(icont(2,k)-icont(2,j)).le.2 .and.
216 & econt(k).lt.econt(i) ) goto 21
221 icont(1,k-1)=icont(1,k)
222 icont(2,k-1)=icont(2,k)
226 c write (iout,*) "ncont",ncont
228 c write (iout,*) icont(1,k),icont(2,k)
239 write (iout,*) 'Electrostatic contacts after pruning: '
245 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
246 & i,restyp(it1),i1,restyp(it2),i2,econt(i)
251 c--------------------------------------------
252 subroutine secondary2(lprint)
253 implicit real*8 (a-h,o-z)
255 include 'COMMON.CHAIN'
256 include 'COMMON.IOUNITS'
257 include 'COMMON.FRAG'
260 include 'COMMON.CONTROL'
261 integer ncont,icont(2,maxcont),isec(maxres,4),nsec(maxres)
262 logical lprint,not_done,freeres
263 double precision p1,p2
266 cc???? if(.not.dccart) call chainbuild
267 cd call write_pdb(99,'sec structure',0d0)
277 call elecont(lprint,ncont,icont)
279 c finding parallel beta
280 cd write (iout,*) '------- looking for parallel beta -----------'
286 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
289 cd write (iout,*) i1,j1
295 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and.
296 & freeres(i1,j1,nsec,isec)) goto 5
300 cd write (iout,*) i1,j1,not_done
304 if (i1-ii1.gt.1) then
308 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
309 & nbeta,ii1,i1,jj1,j1
312 bfrag(1,nbfrag)=ii1+1
314 bfrag(3,nbfrag)=jj1+1
315 bfrag(4,nbfrag)=min0(j1+1,nres)
319 isec(ij,nsec(ij))=nbeta
323 isec(ij,nsec(ij))=nbeta
329 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
330 & "DefPropRes 'strand",nstrand,
331 & "' 'num = ",ii1-1,"..",i1-1,"'"
333 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
334 & "DefPropRes 'strand",nstrand,
335 & "' 'num = ",ii1-1,"..",i1-1,"'"
339 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
340 & "DefPropRes 'strand",nstrand,
341 & "' 'num = ",jj1-1,"..",j1-1,"'"
343 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
344 & "DefPropRes 'strand",nstrand,
345 & "' 'num = ",jj1-1,"..",j1-1,"'"
348 & "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
354 c finding alpha or 310 helix
362 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
366 & ((p1.ge.10.and.p1.le.80).or.i1.le.2).and.
367 & ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
368 cd if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
369 co if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
372 if (nsec(ii1).eq.0) then
381 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
387 if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80)
389 cd write (iout,*) i1,j1,not_done,p1,p2
392 if (j1-ii1.gt.5) then
394 cd write (iout,*)'helix',nhelix,ii1,j1
404 write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
405 if (nhelix.le.9) then
406 write(12,'(a17,i1,a9,i3,a2,i3,a1)')
407 & "DefPropRes 'helix",nhelix,
408 & "' 'num = ",ii1-1,"..",j1-2,"'"
410 write(12,'(a17,i2,a9,i3,a2,i3,a1)')
411 & "DefPropRes 'helix",nhelix,
412 & "' 'num = ",ii1-1,"..",j1-2,"'"
419 if (nhelix.gt.0.and.lprint) then
420 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
422 if (nhelix.le.9) then
423 write(12,'(a8,i1,$)') " | helix",i
425 write(12,'(a8,i2,$)') " | helix",i
432 c finding antiparallel beta
433 cd write (iout,*) '--------- looking for antiparallel beta ---------'
438 if (freeres(i1,j1,nsec,isec)) then
441 cd write (iout,*) i1,j1
448 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
449 & freeres(i1,j1,nsec,isec)) goto 6
453 cd write (iout,*) i1,j1,not_done
457 if (i1-ii1.gt.1) then
461 bfrag(2,nbfrag)=min0(i1+1,nres)
462 bfrag(3,nbfrag)=min0(jj1+1,nres)
469 if (nsec(ij).le.2) then
470 isec(ij,nsec(ij))=nbeta
476 if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
477 isec(ij,nsec(ij))=nbeta
483 write (iout,'(a,i3,4i4)')'antiparallel beta',
484 & nbeta,ii1-1,i1,jj1,j1-1
486 if (nstrand.le.9) then
487 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
488 & "DefPropRes 'strand",nstrand,
489 & "' 'num = ",ii1-2,"..",i1-1,"'"
491 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
492 & "DefPropRes 'strand",nstrand,
493 & "' 'num = ",ii1-2,"..",i1-1,"'"
496 if (nstrand.le.9) then
497 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
498 & "DefPropRes 'strand",nstrand,
499 & "' 'num = ",j1-2,"..",jj1-1,"'"
501 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
502 & "DefPropRes 'strand",nstrand,
503 & "' 'num = ",j1-2,"..",jj1-1,"'"
506 & "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
512 if (nstrand.gt.0.and.lprint) then
513 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
516 write(12,'(a9,i1,$)') " | strand",i
518 write(12,'(a9,i2,$)') " | strand",i
527 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
528 write(12,'(a20)') "XMacStand ribbon.mac"
531 write(iout,*) 'UNRES seq:'
533 write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
537 write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
543 c-------------------------------------------------
544 logical function freeres(i,j,nsec,isec)
545 implicit real*8 (a-h,o-z)
547 integer isec(maxres,4),nsec(maxres)
550 if (nsec(i).lt.0.or.nsec(j).lt.0) return
551 if (nsec(i).gt.1.or.nsec(j).gt.1) return
554 if (isec(i,k).eq.isec(j,l)) return