73325f2fa65e9a380e1e9406c1a826983c405528
[unres.git] / elecont.f
1       subroutine elecont(lprint,ncont,icont)
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'COMMON.IOUNITS'
5       include 'COMMON.CHAIN'
6       include 'COMMON.INTERACT'
7       include 'COMMON.LOCAL'
8       include 'COMMON.FFIELD'
9       include 'COMMON.NAMES'
10       logical lprint
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 *
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.      
19 *
20 * as of 7/06/91.
21 *
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."
29       do i=1,2
30         do j=1,2
31         rri=rpp_(i,j)**6
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
36         if (lprint)
37      &  write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),
38      &                               ael3_(i,j)
39         enddo
40       enddo
41       ncont=0
42       ees=0.0
43       evdw=0.0
44       do 1 i=nnt,nct-2
45         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) goto 1
46         xi=c(1,i)
47         yi=c(2,i)
48         zi=c(3,i)
49         dxi=c(1,i+1)-c(1,i)
50         dyi=c(2,i+1)-c(2,i)
51         dzi=c(3,i+1)-c(3,i)
52         xmedi=xi+0.5*dxi
53         ymedi=yi+0.5*dyi
54         zmedi=zi+0.5*dzi
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
61         do 4 j=i+2,nct-1
62           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
63           ind=ind+1
64           iteli=itel(i)
65           itelj=itel(j)
66           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
67           if (iteli.eq.2 .and. itelj.eq.2) goto 4
68           aaa=app_(iteli,itelj)
69           bbb=bpp_(iteli,itelj)
70           ael6_i=ael6_(iteli,itelj)
71           ael3_i=ael3_(iteli,itelj) 
72           dxj=c(1,j+1)-c(1,j)
73           dyj=c(2,j+1)-c(2,j)
74           dzj=c(3,j+1)-c(3,j)
75           xj=c(1,j)+0.5*dxj
76           yj=c(2,j)+0.5*dyj
77           zj=c(3,j)+0.5*dzj
78           xj=mod(xj,boxxsize)
79           if (xj.lt.0) xj=xj+boxxsize
80           yj=mod(yj,boxysize)
81           if (yj.lt.0) yj=yj+boxysize
82           zj=mod(zj,boxzsize)
83           if (zj.lt.0) zj=zj+boxzsize
84       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
85       xj_safe=xj
86       yj_safe=yj
87       zj_safe=zj
88       isubchap=0
89       do xshift=-1,1
90       do yshift=-1,1
91       do zshift=-1,1
92           xj=xj_safe+xshift*boxxsize
93           yj=yj_safe+yshift*boxysize
94           zj=zj_safe+zshift*boxzsize
95           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
96           if(dist_temp.lt.dist_init) then
97             dist_init=dist_temp
98             xj_temp=xj
99             yj_temp=yj
100             zj_temp=zj
101             isubchap=1
102           endif
103        enddo
104        enddo
105        enddo
106        if (isubchap.eq.1) then
107           xj=xj_temp-xmedi
108           yj=yj_temp-ymedi
109           zj=zj_temp-zmedi
110        else
111           xj=xj_safe-xmedi
112           yj=yj_safe-ymedi
113           zj=zj_safe-zmedi
114        endif
115           rij=xj*xj+yj*yj+zj*zj
116             sss=sscale(sqrt(rij))
117             sssgrad=sscagrad(sqrt(rij))
118           rrmij=1.0/(xj*xj+yj*yj+zj*zj)
119           rmij=sqrt(rrmij)
120           r3ij=rrmij*rmij
121           r6ij=r3ij*r3ij  
122           vrmij=vblinv*rmij
123           cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2      
124           cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
125           cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
126           fac=cosa-3.0*cosb*cosg
127           ev1=aaa*r6ij*r6ij
128           ev2=bbb*r6ij
129           fac3=ael6_i*r6ij
130           fac4=ael3_i*r3ij
131           evdwij=ev1+ev2
132           el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
133           el2=fac4*fac       
134           eesij=el1+el2
135           if (j.gt.i+2 .and. eesij.le.elcutoff .or.
136      &        j.eq.i+2 .and. eesij.le.elecutoff_14) then
137              ncont=ncont+1
138              icont(1,ncont)=i
139              icont(2,ncont)=j
140              econt(ncont)=eesij
141           endif
142           ees=ees+eesij
143           evdw=evdw+evdwij*sss
144     4   continue
145     1 continue
146       if (lprint) then
147         write (iout,*) 'Total average electrostatic energy: ',ees
148         write (iout,*) 'VDW energy between peptide-group centers: ',evdw
149         write (iout,*)
150         write (iout,*) 'Electrostatic contacts before pruning: '
151         do i=1,ncont
152           i1=icont(1,i)
153           i2=icont(2,i)
154           it1=itype(i1)
155           it2=itype(i2)
156           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
157      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
158         enddo
159       endif
160 c For given residues keep only the contacts with the greatest energy.
161       i=0
162       do while (i.lt.ncont)
163         i=i+1
164         ene=econt(i)
165         ic1=icont(1,i)
166         ic2=icont(2,i)
167         j=i
168         do while (j.lt.ncont)
169           j=j+1
170           if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
171      &        ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
172 c            write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
173 c     &       " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
174             if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
175               if (ic1.eq.icont(1,j)) then
176                 do k=1,ncont
177                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
178      &               .and. iabs(icont(1,k)-ic1).le.2 .and. 
179      &               econt(k).lt.econt(j) ) goto 21 
180                 enddo
181               else if (ic2.eq.icont(2,j) ) then
182                 do k=1,ncont
183                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
184      &               .and. iabs(icont(2,k)-ic2).le.2 .and. 
185      &               econt(k).lt.econt(j) ) goto 21 
186                 enddo
187               endif
188 c Remove ith contact
189               do k=i+1,ncont
190                 icont(1,k-1)=icont(1,k)
191                 icont(2,k-1)=icont(2,k)
192                 econt(k-1)=econt(k) 
193               enddo
194               i=i-1
195               ncont=ncont-1
196 c              write (iout,*) "ncont",ncont
197 c              do k=1,ncont
198 c                write (iout,*) icont(1,k),icont(2,k)
199 c              enddo
200               goto 20
201             else if (econt(j).gt.ene .and. ic2.ne.ic1+2) 
202      &      then
203               if (ic1.eq.icont(1,j)) then
204                 do k=1,ncont
205                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
206      &               .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. 
207      &               econt(k).lt.econt(i) ) goto 21 
208                 enddo
209               else if (ic2.eq.icont(2,j) ) then
210                 do k=1,ncont
211                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
212      &               .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. 
213      &               econt(k).lt.econt(i) ) goto 21 
214                 enddo
215               endif
216 c Remove jth contact
217               do k=j+1,ncont
218                 icont(1,k-1)=icont(1,k)
219                 icont(2,k-1)=icont(2,k)
220                 econt(k-1)=econt(k) 
221               enddo
222               ncont=ncont-1
223 c              write (iout,*) "ncont",ncont
224 c              do k=1,ncont
225 c                write (iout,*) icont(1,k),icont(2,k)
226 c              enddo
227               j=j-1
228             endif   
229           endif
230    21     continue
231         enddo
232    20   continue
233       enddo
234       if (lprint) then
235         write (iout,*)
236         write (iout,*) 'Electrostatic contacts after pruning: '
237         do i=1,ncont
238           i1=icont(1,i)
239           i2=icont(2,i)
240           it1=itype(i1)
241           it2=itype(i2)
242           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
243      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
244         enddo
245       endif
246       return
247       end
248 c--------------------------------------------
249       subroutine secondary2(lprint)
250       implicit real*8 (a-h,o-z)
251       include 'DIMENSIONS'
252       include 'COMMON.CHAIN'
253       include 'COMMON.IOUNITS'
254       include 'COMMON.DISTFIT'
255       include 'COMMON.VAR'
256       include 'COMMON.GEO'
257       include 'COMMON.CONTROL'
258       integer ncont,icont(2,maxcont),isec(maxres,4),nsec(maxres)
259       logical lprint,not_done,freeres
260       double precision p1,p2
261       external freeres
262
263 cc????      if(.not.dccart) call chainbuild
264 cd      call write_pdb(99,'sec structure',0d0)
265       ncont=0
266       nbfrag=0
267       nhfrag=0
268       do i=1,nres
269         isec(i,1)=0
270         isec(i,2)=0
271         nsec(i)=0
272       enddo
273
274       call elecont(lprint,ncont,icont)
275
276 c finding parallel beta
277 cd      write (iout,*) '------- looking for parallel beta -----------'
278       nbeta=0
279       nstrand=0
280       do i=1,ncont
281         i1=icont(1,i)
282         j1=icont(2,i)
283         if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
284           ii1=i1
285           jj1=j1
286 cd          write (iout,*) i1,j1
287           not_done=.true.
288           do while (not_done)
289            i1=i1+1
290            j1=j1+1
291             do j=1,ncont
292               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and.
293      &             freeres(i1,j1,nsec,isec)) goto 5
294             enddo
295             not_done=.false.
296   5         continue
297 cd            write (iout,*) i1,j1,not_done
298           enddo
299           j1=j1-1
300           i1=i1-1
301           if (i1-ii1.gt.1) then
302             ii1=max0(ii1-1,1)
303             jj1=max0(jj1-1,1)
304             nbeta=nbeta+1
305             if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
306      &               nbeta,ii1,i1,jj1,j1
307
308             nbfrag=nbfrag+1
309             bfrag(1,nbfrag)=ii1+1
310             bfrag(2,nbfrag)=i1+1
311             bfrag(3,nbfrag)=jj1+1
312             bfrag(4,nbfrag)=min0(j1+1,nres) 
313
314             do ij=ii1,i1
315              nsec(ij)=nsec(ij)+1
316              isec(ij,nsec(ij))=nbeta
317             enddo
318             do ij=jj1,j1
319              nsec(ij)=nsec(ij)+1
320              isec(ij,nsec(ij))=nbeta
321             enddo
322
323            if(lprint) then 
324             nstrand=nstrand+1
325             if (nbeta.le.9) then
326               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
327      &          "DefPropRes 'strand",nstrand,
328      &          "' 'num = ",ii1-1,"..",i1-1,"'"
329             else
330               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
331      &          "DefPropRes 'strand",nstrand,
332      &          "' 'num = ",ii1-1,"..",i1-1,"'"
333             endif
334             nstrand=nstrand+1
335             if (nbeta.le.9) then
336               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
337      &          "DefPropRes 'strand",nstrand,
338      &          "' 'num = ",jj1-1,"..",j1-1,"'"
339             else
340               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
341      &          "DefPropRes 'strand",nstrand,
342      &          "' 'num = ",jj1-1,"..",j1-1,"'"
343             endif
344               write(12,'(a8,4i4)')
345      &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
346            endif
347           endif
348         endif
349       enddo
350
351 c finding alpha or 310 helix
352
353       nhelix=0
354       do i=1,ncont
355         i1=icont(1,i)
356         j1=icont(2,i)
357         p1=phi(i1+2)*rad2deg
358         p2=0.0
359         if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
360
361
362         if (j1.eq.i1+3 .and. 
363      &       ((p1.ge.10.and.p1.le.80).or.i1.le.2).and.
364      &       ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
365 cd          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
366 co          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
367           ii1=i1
368           jj1=j1
369           if (nsec(ii1).eq.0) then 
370             not_done=.true.
371           else
372             not_done=.false.
373           endif
374           do while (not_done)
375             i1=i1+1
376             j1=j1+1
377             do j=1,ncont
378               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
379             enddo
380             not_done=.false.
381   10        continue
382             p1=phi(i1+2)*rad2deg
383             p2=phi(j1+2)*rad2deg
384             if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) 
385      &                              not_done=.false.
386 cd           write (iout,*) i1,j1,not_done,p1,p2
387           enddo
388           j1=j1+1
389           if (j1-ii1.gt.5) then
390             nhelix=nhelix+1
391 cd            write (iout,*)'helix',nhelix,ii1,j1
392
393             nhfrag=nhfrag+1
394             hfrag(1,nhfrag)=ii1
395             hfrag(2,nhfrag)=j1
396
397             do ij=ii1,j1
398              nsec(ij)=-1
399             enddo
400            if (lprint) then
401             write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
402             if (nhelix.le.9) then
403               write(12,'(a17,i1,a9,i3,a2,i3,a1)') 
404      &          "DefPropRes 'helix",nhelix,
405      &          "' 'num = ",ii1-1,"..",j1-2,"'"
406             else
407               write(12,'(a17,i2,a9,i3,a2,i3,a1)') 
408      &          "DefPropRes 'helix",nhelix,
409      &          "' 'num = ",ii1-1,"..",j1-2,"'"
410             endif
411            endif
412           endif
413         endif
414       enddo
415        
416       if (nhelix.gt.0.and.lprint) then
417         write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
418         do i=2,nhelix
419          if (nhelix.le.9) then
420           write(12,'(a8,i1,$)') " | helix",i
421          else
422           write(12,'(a8,i2,$)') " | helix",i
423          endif
424         enddo
425         write(12,'(a1)') "'"
426       endif
427
428
429 c finding antiparallel beta
430 cd      write (iout,*) '--------- looking for antiparallel beta ---------'
431
432       do i=1,ncont
433         i1=icont(1,i)
434         j1=icont(2,i)
435         if (freeres(i1,j1,nsec,isec)) then
436           ii1=i1
437           jj1=j1
438 cd          write (iout,*) i1,j1
439
440           not_done=.true.
441           do while (not_done)
442            i1=i1+1
443            j1=j1-1
444             do j=1,ncont
445               if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
446      &             freeres(i1,j1,nsec,isec)) goto 6
447             enddo
448             not_done=.false.
449   6         continue
450 cd            write (iout,*) i1,j1,not_done
451           enddo
452           i1=i1-1
453           j1=j1+1
454           if (i1-ii1.gt.1) then
455
456             nbfrag=nbfrag+1
457             bfrag(1,nbfrag)=ii1
458             bfrag(2,nbfrag)=min0(i1+1,nres)
459             bfrag(3,nbfrag)=min0(jj1+1,nres)
460             bfrag(4,nbfrag)=j1
461
462             nbeta=nbeta+1
463             iii1=max0(ii1-1,1)
464             do ij=iii1,i1
465              nsec(ij)=nsec(ij)+1
466              if (nsec(ij).le.2) then
467               isec(ij,nsec(ij))=nbeta
468              endif
469             enddo
470             jjj1=max0(j1-1,1)  
471             do ij=jjj1,jj1
472              nsec(ij)=nsec(ij)+1
473              if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
474               isec(ij,nsec(ij))=nbeta
475              endif
476             enddo
477
478
479            if (lprint) then
480             write (iout,'(a,i3,4i4)')'antiparallel beta',
481      &                   nbeta,ii1-1,i1,jj1,j1-1
482             nstrand=nstrand+1
483             if (nstrand.le.9) then
484               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
485      &          "DefPropRes 'strand",nstrand,
486      &          "' 'num = ",ii1-2,"..",i1-1,"'"
487             else
488               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
489      &          "DefPropRes 'strand",nstrand,
490      &          "' 'num = ",ii1-2,"..",i1-1,"'"
491             endif
492             nstrand=nstrand+1
493             if (nstrand.le.9) then
494               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
495      &          "DefPropRes 'strand",nstrand,
496      &          "' 'num = ",j1-2,"..",jj1-1,"'"
497             else
498               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
499      &          "DefPropRes 'strand",nstrand,
500      &          "' 'num = ",j1-2,"..",jj1-1,"'"
501             endif
502               write(12,'(a8,4i4)')
503      &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
504            endif
505           endif
506         endif
507       enddo
508
509       if (nstrand.gt.0.and.lprint) then
510         write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
511         do i=2,nstrand
512          if (i.le.9) then
513           write(12,'(a9,i1,$)') " | strand",i
514          else
515           write(12,'(a9,i2,$)') " | strand",i
516          endif
517         enddo
518         write(12,'(a1)') "'"
519       endif
520
521        
522
523       if (lprint) then
524        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
525        write(12,'(a20)') "XMacStand ribbon.mac"
526          
527         
528        write(iout,*) 'UNRES seq:'
529        do j=1,nbfrag
530         write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
531        enddo
532
533        do j=1,nhfrag
534         write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
535        enddo
536       endif       
537
538       return
539       end
540 c-------------------------------------------------
541       logical function freeres(i,j,nsec,isec)
542       implicit real*8 (a-h,o-z)
543       include 'DIMENSIONS'
544       integer isec(maxres,4),nsec(maxres)
545       freeres=.false.
546
547       if (nsec(i).lt.0.or.nsec(j).lt.0) return
548       if (nsec(i).gt.1.or.nsec(j).gt.1) return
549       do k=1,nsec(i)
550         do l=1,nsec(j)
551           if (isec(i,k).eq.isec(j,l)) return
552         enddo
553       enddo
554       freeres=.true.
555       return
556       end
557