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