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)
15 integer xshift,yshift,zshift
17 * Load the constants of peptide bond - peptide bond interactions.
18 * Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
19 * proline) - determined by averaging ECEPP energy.
23 c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
24 data rpp_ / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
25 data elpp_6 /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
26 data elpp_3 / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
27 data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
28 if (lprint) write (iout,'(a)')
29 & "Constants of electrostatic interaction energy expression."
33 app_(i,j)=epp(i,j)*rri*rri
34 bpp_(i,j)=-2.0*epp(i,j)*rri
35 ael6_(i,j)=elpp_6(i,j)*4.2**6
36 ael3_(i,j)=elpp_3(i,j)*4.2**3
38 & write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),
46 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) goto 1
56 c write (iout,*) "i",xmedi,ymedi,zmedi
57 xmedi=mod(xmedi,boxxsize)
58 if (xmedi.lt.0) xmedi=xmedi+boxxsize
59 ymedi=mod(ymedi,boxysize)
60 if (ymedi.lt.0) ymedi=ymedi+boxysize
61 zmedi=mod(zmedi,boxzsize)
62 if (zmedi.lt.0) zmedi=zmedi+boxzsize
63 c write (iout,*) "i",xmedi,ymedi,zmedi
65 c write (iout,*) "i",i," j",j
66 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
70 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
71 if (iteli.eq.2 .and. itelj.eq.2) goto 4
74 ael6_i=ael6_(iteli,itelj)
75 ael3_i=ael3_(iteli,itelj)
82 c write (iout,*) "j",xj,yj,zj
84 if (xj.lt.0) xj=xj+boxxsize
86 if (yj.lt.0) yj=yj+boxysize
88 if (zj.lt.0) zj=zj+boxzsize
89 c write (iout,*) "j",xj,yj,zj
90 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
91 c write (iout,*) "dist",dsqrt(dist_init)
99 xj=xj_safe+xshift*boxxsize
100 yj=yj_safe+yshift*boxysize
101 zj=zj_safe+zshift*boxzsize
102 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
103 c write (iout,*) "shift",xshift,yshift,zshift," dist_temp",
104 c & dist_temp," dist_init",dist_init
105 if(dist_temp.lt.dist_init) then
115 if (isubchap.eq.1) then
124 rij=xj*xj+yj*yj+zj*zj
125 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
130 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
131 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
132 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
133 fac=cosa-3.0*cosb*cosg
139 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
142 if (j.gt.i+2 .and. eesij.le.elcutoff .or.
143 & j.eq.i+2 .and. eesij.le.elecutoff_14) then
150 c write (iout,*) "i"," j",j," rij",dsqrt(rij)," eesij",eesij
154 write (iout,*) 'Total average electrostatic energy: ',ees
155 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
157 write (iout,*) 'Electrostatic contacts before pruning: '
163 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
164 & i,restyp(it1),i1,restyp(it2),i2,econt(i)
167 c For given residues keep only the contacts with the greatest energy.
169 do while (i.lt.ncont)
175 do while (j.lt.ncont)
177 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
178 & ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
179 c write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
180 c & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
181 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
182 if (ic1.eq.icont(1,j)) then
184 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
185 & .and. iabs(icont(1,k)-ic1).le.2 .and.
186 & econt(k).lt.econt(j) ) goto 21
188 else if (ic2.eq.icont(2,j) ) then
190 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
191 & .and. iabs(icont(2,k)-ic2).le.2 .and.
192 & econt(k).lt.econt(j) ) goto 21
197 icont(1,k-1)=icont(1,k)
198 icont(2,k-1)=icont(2,k)
203 c write (iout,*) "ncont",ncont
205 c write (iout,*) icont(1,k),icont(2,k)
208 else if (econt(j).gt.ene .and. ic2.ne.ic1+2)
210 if (ic1.eq.icont(1,j)) then
212 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
213 & .and. iabs(icont(1,k)-icont(1,j)).le.2 .and.
214 & econt(k).lt.econt(i) ) goto 21
216 else if (ic2.eq.icont(2,j) ) then
218 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
219 & .and. iabs(icont(2,k)-icont(2,j)).le.2 .and.
220 & econt(k).lt.econt(i) ) goto 21
225 icont(1,k-1)=icont(1,k)
226 icont(2,k-1)=icont(2,k)
230 c write (iout,*) "ncont",ncont
232 c write (iout,*) icont(1,k),icont(2,k)
243 write (iout,*) 'Electrostatic contacts after pruning: '
249 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
250 & i,restyp(it1),i1,restyp(it2),i2,econt(i)
255 c--------------------------------------------
256 subroutine secondary2(lprint)
257 implicit real*8 (a-h,o-z)
259 include 'COMMON.CHAIN'
260 include 'COMMON.IOUNITS'
261 include 'COMMON.FRAG'
264 include 'COMMON.CONTROL'
265 integer ncont,icont(2,maxcont),isec(maxres,4),nsec(maxres)
266 logical lprint,not_done,freeres
267 double precision p1,p2
270 cc???? if(.not.dccart) call chainbuild
271 cd call write_pdb(99,'sec structure',0d0)
281 call elecont(lprint,ncont,icont)
283 c finding parallel beta
284 cd write (iout,*) '------- looking for parallel beta -----------'
290 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
293 cd write (iout,*) i1,j1
299 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and.
300 & freeres(i1,j1,nsec,isec)) goto 5
304 cd write (iout,*) i1,j1,not_done
308 if (i1-ii1.gt.1) then
312 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
313 & nbeta,ii1,i1,jj1,j1
316 bfrag(1,nbfrag)=ii1+1
318 bfrag(3,nbfrag)=jj1+1
319 bfrag(4,nbfrag)=min0(j1+1,nres)
323 isec(ij,nsec(ij))=nbeta
327 isec(ij,nsec(ij))=nbeta
333 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
334 & "DefPropRes 'strand",nstrand,
335 & "' 'num = ",ii1-1,"..",i1-1,"'"
337 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
338 & "DefPropRes 'strand",nstrand,
339 & "' 'num = ",ii1-1,"..",i1-1,"'"
343 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
344 & "DefPropRes 'strand",nstrand,
345 & "' 'num = ",jj1-1,"..",j1-1,"'"
347 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
348 & "DefPropRes 'strand",nstrand,
349 & "' 'num = ",jj1-1,"..",j1-1,"'"
352 & "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
358 c finding alpha or 310 helix
366 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
370 & ((p1.ge.10.and.p1.le.80).or.i1.le.2).and.
371 & ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
372 cd if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
373 co if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
376 if (nsec(ii1).eq.0) then
385 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
391 if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80)
393 cd write (iout,*) i1,j1,not_done,p1,p2
396 if (j1-ii1.gt.5) then
398 cd write (iout,*)'helix',nhelix,ii1,j1
408 write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
409 if (nhelix.le.9) then
410 write(12,'(a17,i1,a9,i3,a2,i3,a1)')
411 & "DefPropRes 'helix",nhelix,
412 & "' 'num = ",ii1-1,"..",j1-2,"'"
414 write(12,'(a17,i2,a9,i3,a2,i3,a1)')
415 & "DefPropRes 'helix",nhelix,
416 & "' 'num = ",ii1-1,"..",j1-2,"'"
423 if (nhelix.gt.0.and.lprint) then
424 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
426 if (nhelix.le.9) then
427 write(12,'(a8,i1,$)') " | helix",i
429 write(12,'(a8,i2,$)') " | helix",i
436 c finding antiparallel beta
437 cd write (iout,*) '--------- looking for antiparallel beta ---------'
442 if (freeres(i1,j1,nsec,isec)) then
445 cd write (iout,*) i1,j1
452 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
453 & freeres(i1,j1,nsec,isec)) goto 6
457 cd write (iout,*) i1,j1,not_done
461 if (i1-ii1.gt.1) then
465 bfrag(2,nbfrag)=min0(i1+1,nres)
466 bfrag(3,nbfrag)=min0(jj1+1,nres)
473 if (nsec(ij).le.2) then
474 isec(ij,nsec(ij))=nbeta
480 if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
481 isec(ij,nsec(ij))=nbeta
487 write (iout,'(a,i3,4i4)')'antiparallel beta',
488 & nbeta,ii1-1,i1,jj1,j1-1
490 if (nstrand.le.9) then
491 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
492 & "DefPropRes 'strand",nstrand,
493 & "' 'num = ",ii1-2,"..",i1-1,"'"
495 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
496 & "DefPropRes 'strand",nstrand,
497 & "' 'num = ",ii1-2,"..",i1-1,"'"
500 if (nstrand.le.9) then
501 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
502 & "DefPropRes 'strand",nstrand,
503 & "' 'num = ",j1-2,"..",jj1-1,"'"
505 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
506 & "DefPropRes 'strand",nstrand,
507 & "' 'num = ",j1-2,"..",jj1-1,"'"
510 & "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
516 if (nstrand.gt.0.and.lprint) then
517 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
520 write(12,'(a9,i1,$)') " | strand",i
522 write(12,'(a9,i2,$)') " | strand",i
531 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
532 write(12,'(a20)') "XMacStand ribbon.mac"
535 write(iout,*) 'UNRES seq:'
537 write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
541 write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
547 c-------------------------------------------------
548 logical function freeres(i,j,nsec,isec)
549 implicit real*8 (a-h,o-z)
551 integer isec(maxres,4),nsec(maxres)
554 if (nsec(i).lt.0.or.nsec(j).lt.0) return
555 if (nsec(i).gt.1.or.nsec(j).gt.1) return
558 if (isec(i,k).eq.isec(j,l)) return