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),epp_(2,2),rpp_(2,2)
13 integer ncont,icont(2,maxcont)
14 double precision econt(maxcont)
15 double precision boxshift
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 call to_box(xmedi,ymedi,zmedi)
57 c write (iout,*) "i",xmedi,ymedi,zmedi
59 c write (iout,*) "i",i," j",j
60 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
64 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
65 if (iteli.eq.2 .and. itelj.eq.2) goto 4
68 ael6_i=ael6_(iteli,itelj)
69 ael3_i=ael3_(iteli,itelj)
77 c write (iout,*) "j",xj,yj,zj
78 xj=boxshift(xj-xmedi,boxxsize)
79 yj=boxshift(yj-ymedi,boxysize)
80 zj=boxshift(zj-zmedi,boxzsize)
81 c write (iout,*) "j",xj,yj,zj
83 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
88 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
89 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
90 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
91 fac=cosa-3.0*cosb*cosg
97 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
100 if (j.gt.i+2 .and. eesij.le.elcutoff .or.
101 & j.eq.i+2 .and. eesij.le.elecutoff_14) then
108 c write (iout,*) "i",i," j",j," rij",dsqrt(rij)," eesij",eesij
112 write (iout,*) 'Total average electrostatic energy: ',ees
113 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
115 write (iout,*) 'Electrostatic contacts before pruning: '
121 write (iout,'(i7,2x,a,i6,2x,a,i6,f10.5)')
122 & i,restyp(it1),i1,restyp(it2),i2,econt(i)
125 c For given residues keep only the contacts with the greatest energy.
127 do while (i.lt.ncont)
133 do while (j.lt.ncont)
135 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
136 & ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
137 c write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
138 c & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
139 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
140 if (ic1.eq.icont(1,j)) then
142 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
143 & .and. iabs(icont(1,k)-ic1).le.2 .and.
144 & econt(k).lt.econt(j) ) goto 21
146 else if (ic2.eq.icont(2,j) ) then
148 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
149 & .and. iabs(icont(2,k)-ic2).le.2 .and.
150 & econt(k).lt.econt(j) ) goto 21
155 icont(1,k-1)=icont(1,k)
156 icont(2,k-1)=icont(2,k)
161 c write (iout,*) "ncont",ncont
163 c write (iout,*) icont(1,k),icont(2,k)
166 else if (econt(j).gt.ene .and. ic2.ne.ic1+2)
168 if (ic1.eq.icont(1,j)) then
170 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
171 & .and. iabs(icont(1,k)-icont(1,j)).le.2 .and.
172 & econt(k).lt.econt(i) ) goto 21
174 else if (ic2.eq.icont(2,j) ) then
176 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
177 & .and. iabs(icont(2,k)-icont(2,j)).le.2 .and.
178 & econt(k).lt.econt(i) ) goto 21
183 icont(1,k-1)=icont(1,k)
184 icont(2,k-1)=icont(2,k)
188 c write (iout,*) "ncont",ncont
190 c write (iout,*) icont(1,k),icont(2,k)
201 write (iout,*) 'Electrostatic contacts after pruning: '
207 write (iout,'(i7,2x,a,i6,2x,a,i6,f10.5)')
208 & i,restyp(it1),i1,restyp(it2),i2,econt(i)
213 c--------------------------------------------
214 subroutine secondary2(lprint)
215 implicit real*8 (a-h,o-z)
217 include 'COMMON.CHAIN'
218 include 'COMMON.IOUNITS'
219 include 'COMMON.FRAG'
222 include 'COMMON.CONTROL'
223 integer ncont,icont(2,maxcont),isec(maxres,4),nsec(maxres)
224 logical lprint,not_done,freeres
225 double precision p1,p2
228 cc???? if(.not.dccart) call chainbuild
229 cd call write_pdb(99,'sec structure',0d0)
239 call elecont(lprint,ncont,icont)
241 c finding parallel beta
242 cd write (iout,*) '------- looking for parallel beta -----------'
248 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
251 cd write (iout,*) i1,j1
257 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and.
258 & freeres(i1,j1,nsec,isec)) goto 5
262 cd write (iout,*) i1,j1,not_done
266 if (i1-ii1.gt.1) then
270 if(lprint)write(iout,'(a,i7,4i6)')'parallel beta',
271 & nbeta,ii1,i1,jj1,j1
274 bfrag(1,nbfrag)=ii1+1
276 bfrag(3,nbfrag)=jj1+1
277 bfrag(4,nbfrag)=min0(j1+1,nres)
281 isec(ij,nsec(ij))=nbeta
285 isec(ij,nsec(ij))=nbeta
291 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
292 & "DefPropRes 'strand",nstrand,
293 & "' 'num = ",ii1-1,"..",i1-1,"'"
295 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
296 & "DefPropRes 'strand",nstrand,
297 & "' 'num = ",ii1-1,"..",i1-1,"'"
301 write(12,'(a18,i1,a9,i3,a2,i3,a1)')
302 & "DefPropRes 'strand",nstrand,
303 & "' 'num = ",jj1-1,"..",j1-1,"'"
305 write(12,'(a18,i2,a9,i3,a2,i3,a1)')
306 & "DefPropRes 'strand",nstrand,
307 & "' 'num = ",jj1-1,"..",j1-1,"'"
310 & "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
316 c finding alpha or 310 helix
324 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
328 & ((p1.ge.10.and.p1.le.80).or.i1.le.2).and.
329 & ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
330 cd if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
331 co if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
334 if (nsec(ii1).eq.0) then
343 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
349 if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80)
351 cd write (iout,*) i1,j1,not_done,p1,p2
354 if (j1-ii1.gt.5) then
356 cd write (iout,*)'helix',nhelix,ii1,j1
366 write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
367 if (nhelix.le.9) then
368 write(12,'(a17,i1,a9,i3,a2,i3,a1)')
369 & "DefPropRes 'helix",nhelix,
370 & "' 'num = ",ii1-1,"..",j1-2,"'"
372 write(12,'(a17,i2,a9,i3,a2,i3,a1)')
373 & "DefPropRes 'helix",nhelix,
374 & "' 'num = ",ii1-1,"..",j1-2,"'"
381 if (nhelix.gt.0.and.lprint) then
382 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
384 if (nhelix.le.9) then
385 write(12,'(a8,i1,$)') " | helix",i
387 write(12,'(a8,i2,$)') " | helix",i
394 c finding antiparallel beta
395 cd write (iout,*) '--------- looking for antiparallel beta ---------'
400 if (freeres(i1,j1,nsec,isec)) then
403 cd write (iout,*) i1,j1
410 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
411 & freeres(i1,j1,nsec,isec)) goto 6
415 cd write (iout,*) i1,j1,not_done
419 if (i1-ii1.gt.1) then
423 bfrag(2,nbfrag)=min0(i1+1,nres)
424 bfrag(3,nbfrag)=min0(jj1+1,nres)
431 if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
432 isec(ij,nsec(ij))=nbeta
438 if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
439 isec(ij,nsec(ij))=nbeta
445 write (iout,'(a,i3,4i6)')'antiparallel beta',
446 & nbeta,ii1-1,i1,jj1,j1-1
448 if (nstrand.le.9) then
449 write(12,'(a18,i1,a9,i6,a2,i6,a1)')
450 & "DefPropRes 'strand",nstrand,
451 & "' 'num = ",ii1-2,"..",i1-1,"'"
453 write(12,'(a18,i2,a9,i6,a2,i6,a1)')
454 & "DefPropRes 'strand",nstrand,
455 & "' 'num = ",ii1-2,"..",i1-1,"'"
458 if (nstrand.le.9) then
459 write(12,'(a18,i1,a9,i6,a2,i6,a1)')
460 & "DefPropRes 'strand",nstrand,
461 & "' 'num = ",j1-2,"..",jj1-1,"'"
463 write(12,'(a18,i2,a9,i6,a2,i6,a1)')
464 & "DefPropRes 'strand",nstrand,
465 & "' 'num = ",j1-2,"..",jj1-1,"'"
468 & "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
474 if (nstrand.gt.0.and.lprint) then
475 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
478 write(12,'(a9,i1,$)') " | strand",i
480 write(12,'(a9,i2,$)') " | strand",i
489 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
490 write(12,'(a20)') "XMacStand ribbon.mac"
493 write(iout,*) 'UNRES seq:'
495 write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
499 write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
505 c-------------------------------------------------
506 logical function freeres(i,j,nsec,isec)
507 implicit real*8 (a-h,o-z)
509 integer isec(maxres,4),nsec(maxres)
512 if (nsec(i).lt.0.or.nsec(j).lt.0) return
513 if (nsec(i).gt.1.or.nsec(j).gt.1) return
516 if (isec(i,k).eq.isec(j,l)) return