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