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 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
69 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
70 if (iteli.eq.2 .and. itelj.eq.2) goto 4
73 ael6_i=ael6_(iteli,itelj)
74 ael3_i=ael3_(iteli,itelj)
81 c write (iout,*) "j",xj,yj,zj
83 if (xj.lt.0) xj=xj+boxxsize
85 if (yj.lt.0) yj=yj+boxysize
87 if (zj.lt.0) zj=zj+boxzsize
88 c write (iout,*) "j",xj,yj,zj
89 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
90 c write (iout,*) "dist",dsqrt(dist_init)
98 xj=xj_safe+xshift*boxxsize
99 yj=yj_safe+yshift*boxysize
100 zj=zj_safe+zshift*boxzsize
101 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
102 c write (iout,*) "shift",xshift,yshift,zshift," dist_temp",
103 c & dist_temp," dist_init",dist_init
104 if(dist_temp.lt.dist_init) then
114 if (isubchap.eq.1) then
123 rij=xj*xj+yj*yj+zj*zj
124 sss=sscale(sqrt(rij))
125 sssgrad=sscagrad(sqrt(rij))
126 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
131 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
132 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
133 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
134 fac=cosa-3.0*cosb*cosg
140 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
143 if (j.gt.i+2 .and. eesij.le.elcutoff .or.
144 & j.eq.i+2 .and. eesij.le.elecutoff_14) then
152 c write (iout,*) "i"," j",j," rij",dsqrt(rij)," eesij",eesij
156 write (iout,*) 'Total average electrostatic energy: ',ees
157 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
159 write (iout,*) 'Electrostatic contacts before pruning: '
165 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
166 & i,restyp(it1),i1,restyp(it2),i2,econt(i)
169 c For given residues keep only the contacts with the greatest energy.
171 do while (i.lt.ncont)
177 do while (j.lt.ncont)
179 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
180 & ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
181 c write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
182 c & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
183 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
184 if (ic1.eq.icont(1,j)) then
186 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
187 & .and. iabs(icont(1,k)-ic1).le.2 .and.
188 & econt(k).lt.econt(j) ) goto 21
190 else if (ic2.eq.icont(2,j) ) then
192 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
193 & .and. iabs(icont(2,k)-ic2).le.2 .and.
194 & econt(k).lt.econt(j) ) goto 21
199 icont(1,k-1)=icont(1,k)
200 icont(2,k-1)=icont(2,k)
205 c write (iout,*) "ncont",ncont
207 c write (iout,*) icont(1,k),icont(2,k)
210 else if (econt(j).gt.ene .and. ic2.ne.ic1+2)
212 if (ic1.eq.icont(1,j)) then
214 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
215 & .and. iabs(icont(1,k)-icont(1,j)).le.2 .and.
216 & econt(k).lt.econt(i) ) goto 21
218 else if (ic2.eq.icont(2,j) ) then
220 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
221 & .and. iabs(icont(2,k)-icont(2,j)).le.2 .and.
222 & econt(k).lt.econt(i) ) goto 21
227 icont(1,k-1)=icont(1,k)
228 icont(2,k-1)=icont(2,k)
232 c write (iout,*) "ncont",ncont
234 c write (iout,*) icont(1,k),icont(2,k)
245 write (iout,*) 'Electrostatic contacts after pruning: '
251 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
252 & i,restyp(it1),i1,restyp(it2),i2,econt(i)
257 c--------------------------------------------
258 subroutine secondary2(lprint)
259 implicit real*8 (a-h,o-z)
261 include 'COMMON.CHAIN'
262 include 'COMMON.IOUNITS'
263 include 'COMMON.FRAG'
266 include 'COMMON.CONTROL'
267 integer ncont,icont(2,maxcont),isec(maxres,4),nsec(maxres)
268 logical lprint,not_done,freeres
269 double precision p1,p2
272 cc???? if(.not.dccart) call chainbuild
273 cd call write_pdb(99,'sec structure',0d0)
283 call elecont(lprint,ncont,icont)
285 c finding parallel beta
286 cd write (iout,*) '------- looking for parallel beta -----------'
292 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
295 cd write (iout,*) i1,j1
301 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and.
302 & freeres(i1,j1,nsec,isec)) goto 5
306 cd write (iout,*) i1,j1,not_done
310 if (i1-ii1.gt.1) then
314 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
315 & nbeta,ii1,i1,jj1,j1
318 bfrag(1,nbfrag)=ii1+1
320 bfrag(3,nbfrag)=jj1+1
321 bfrag(4,nbfrag)=min0(j1+1,nres)
325 isec(ij,nsec(ij))=nbeta
329 isec(ij,nsec(ij))=nbeta
335 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
336 & "DefPropRes 'strand",nstrand,
337 & "' 'num = ",ii1-1,"..",i1-1,"'"
339 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
340 & "DefPropRes 'strand",nstrand,
341 & "' 'num = ",ii1-1,"..",i1-1,"'"
345 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
346 & "DefPropRes 'strand",nstrand,
347 & "' 'num = ",jj1-1,"..",j1-1,"'"
349 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
350 & "DefPropRes 'strand",nstrand,
351 & "' 'num = ",jj1-1,"..",j1-1,"'"
354 & "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
360 c finding alpha or 310 helix
368 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
372 & ((p1.ge.10.and.p1.le.80).or.i1.le.2).and.
373 & ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
374 cd if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
375 co if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
378 if (nsec(ii1).eq.0) then
387 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
393 if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80)
395 cd write (iout,*) i1,j1,not_done,p1,p2
398 if (j1-ii1.gt.5) then
400 cd write (iout,*)'helix',nhelix,ii1,j1
410 write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
411 if (nhelix.le.9) then
412 write(12,'(a17,i1,a9,i3,a2,i3,a1)')
413 & "DefPropRes 'helix",nhelix,
414 & "' 'num = ",ii1-1,"..",j1-2,"'"
416 write(12,'(a17,i2,a9,i3,a2,i3,a1)')
417 & "DefPropRes 'helix",nhelix,
418 & "' 'num = ",ii1-1,"..",j1-2,"'"
425 if (nhelix.gt.0.and.lprint) then
426 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
428 if (nhelix.le.9) then
429 write(12,'(a8,i1,$)') " | helix",i
431 write(12,'(a8,i2,$)') " | helix",i
438 c finding antiparallel beta
439 cd write (iout,*) '--------- looking for antiparallel beta ---------'
444 if (freeres(i1,j1,nsec,isec)) then
447 cd write (iout,*) i1,j1
454 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
455 & freeres(i1,j1,nsec,isec)) goto 6
459 cd write (iout,*) i1,j1,not_done
463 if (i1-ii1.gt.1) then
467 bfrag(2,nbfrag)=min0(i1+1,nres)
468 bfrag(3,nbfrag)=min0(jj1+1,nres)
475 if (nsec(ij).le.2) then
476 isec(ij,nsec(ij))=nbeta
482 if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
483 isec(ij,nsec(ij))=nbeta
489 write (iout,'(a,i3,4i4)')'antiparallel beta',
490 & nbeta,ii1-1,i1,jj1,j1-1
492 if (nstrand.le.9) then
493 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
494 & "DefPropRes 'strand",nstrand,
495 & "' 'num = ",ii1-2,"..",i1-1,"'"
497 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
498 & "DefPropRes 'strand",nstrand,
499 & "' 'num = ",ii1-2,"..",i1-1,"'"
502 if (nstrand.le.9) then
503 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
504 & "DefPropRes 'strand",nstrand,
505 & "' 'num = ",j1-2,"..",jj1-1,"'"
507 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
508 & "DefPropRes 'strand",nstrand,
509 & "' 'num = ",j1-2,"..",jj1-1,"'"
512 & "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
518 if (nstrand.gt.0.and.lprint) then
519 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
522 write(12,'(a9,i1,$)') " | strand",i
524 write(12,'(a9,i2,$)') " | strand",i
533 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
534 write(12,'(a20)') "XMacStand ribbon.mac"
537 write(iout,*) 'UNRES seq:'
539 write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
543 write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
549 c-------------------------------------------------
550 logical function freeres(i,j,nsec,isec)
551 implicit real*8 (a-h,o-z)
553 integer isec(maxres,4),nsec(maxres)
556 if (nsec(i).lt.0.or.nsec(j).lt.0) return
557 if (nsec(i).gt.1.or.nsec(j).gt.1) return
560 if (isec(i,k).eq.isec(j,l)) return