Adam's changes
[unres.git] / source / wham / src-restraints-PM / elecont.f
1       subroutine elecont(lprint,ncont,icont,ist,ien)
2       implicit none
3       include 'DIMENSIONS'
4       include 'DIMENSIONS.ZSCOPT'
5       include 'DIMENSIONS.COMPAR'
6       include 'COMMON.IOUNITS'
7       include 'COMMON.CHAIN'
8       include 'COMMON.INTERACT'
9       include 'COMMON.FFIELD'
10       include 'COMMON.NAMES'
11       include 'COMMON.LOCAL'
12       logical lprint
13       integer i,j,k,ist,ien,iteli,itelj,ind,i1,i2,it1,it2,ic1,ic2
14       double precision rri,xi,yi,zi,dxi,dyi,dzi,xmedi,ymedi,zmedi,
15      &  xj,yj,zj,dxj,dyj,dzj,aaa,bbb,ael6i,ael3i,rrmij,rmij,r3ij,r6ij,
16      &  vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,evdwij,el1,el2,
17      &  eesij,ees,evdw,ene
18       double precision elpp6c(2,2),elpp3c(2,2),ael6c(2,2),ael3c(2,2),
19      &  appc(2,2),bppc(2,2)
20       double precision elcutoff,elecutoff_14
21       integer ncont,icont(2,maxcont)
22       double precision econt(maxcont)
23 *
24 * Load the constants of peptide bond - peptide bond interactions.
25 * Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
26 * proline) - determined by averaging ECEPP energy.      
27 *
28 * as of 7/06/91.
29 *
30 c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
31 c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
32       data elpp6c  /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
33       data elpp3c  / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
34       data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
35       ees=0.0d0
36       evdw=0.0d0
37       if (lprint) write (iout,'(a)') 
38      &  "Constants of electrostatic interaction energy expression."
39       do i=1,2
40         do j=1,2
41         rri=rpp(i,j)**6
42         appc(i,j)=epp(i,j)*rri*rri 
43         bppc(i,j)=-2.0*epp(i,j)*rri
44         ael6c(i,j)=elpp6c(i,j)*4.2**6
45         ael3c(i,j)=elpp3c(i,j)*4.2**3
46         if (lprint)
47      &  write (iout,'(2i2,4e15.4)') i,j,appc(i,j),bppc(i,j),ael6c(i,j),
48      &                               ael3c(i,j)
49         enddo
50       enddo
51       ncont=0
52       do 1 i=ist,ien-2
53         xi=c(1,i)
54         yi=c(2,i)
55         zi=c(3,i)
56         dxi=c(1,i+1)-c(1,i)
57         dyi=c(2,i+1)-c(2,i)
58         dzi=c(3,i+1)-c(3,i)
59         xmedi=xi+0.5*dxi
60         ymedi=yi+0.5*dyi
61         zmedi=zi+0.5*dzi
62         do 4 j=i+2,ien-1
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=appc(iteli,itelj)
69           bbb=bppc(iteli,itelj)
70           ael6i=ael6c(iteli,itelj)
71           ael3i=ael3c(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-xmedi
76           yj=c(2,j)+0.5*dyj-ymedi
77           zj=c(3,j)+0.5*dzj-zmedi
78           rrmij=1.0/(xj*xj+yj*yj+zj*zj)
79           rmij=sqrt(rrmij)
80           r3ij=rrmij*rmij
81           r6ij=r3ij*r3ij  
82           vrmij=vblinv*rmij
83           cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2      
84           cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
85           cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
86           fac=cosa-3.0*cosb*cosg
87           ev1=aaa*r6ij*r6ij
88           ev2=bbb*r6ij
89           fac3=ael6i*r6ij
90           fac4=ael3i*r3ij
91           evdwij=ev1+ev2
92           el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
93           el2=fac4*fac       
94           eesij=el1+el2
95           if (j.gt.i+2 .and. eesij.le.elcutoff .or.
96      &        j.eq.i+2 .and. eesij.le.elecutoff_14) then
97              ncont=ncont+1
98              icont(1,ncont)=i
99              icont(2,ncont)=j
100              econt(ncont)=eesij
101           endif
102           ees=ees+eesij
103           evdw=evdw+evdwij
104     4   continue
105     1 continue
106       if (lprint) then
107         write (iout,*) 'Total average electrostatic energy: ',ees
108         write (iout,*) 'VDW energy between peptide-group centers: ',evdw
109         write (iout,*)
110         write (iout,*) 'Electrostatic contacts before pruning: '
111         do i=1,ncont
112           i1=icont(1,i)
113           i2=icont(2,i)
114           it1=itype(i1)
115           it2=itype(i2)
116           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
117      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
118         enddo
119       endif
120 c For given residues keep only the contacts with the greatest energy.
121       i=0
122       do while (i.lt.ncont)
123         i=i+1
124         ene=econt(i)
125         ic1=icont(1,i)
126         ic2=icont(2,i)
127         j=i
128         do while (j.lt.ncont)
129           j=j+1
130           if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
131      &        ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
132 c            write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
133 c     &       " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
134             if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
135               if (ic1.eq.icont(1,j)) then
136                 do k=1,ncont
137                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
138      &               .and. iabs(icont(1,k)-ic1).le.2 .and. 
139      &               econt(k).lt.econt(j) ) goto 21 
140                 enddo
141               else if (ic2.eq.icont(2,j) ) then
142                 do k=1,ncont
143                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
144      &               .and. iabs(icont(2,k)-ic2).le.2 .and. 
145      &               econt(k).lt.econt(j) ) goto 21 
146                 enddo
147               endif
148 c Remove ith contact
149               do k=i+1,ncont
150                 icont(1,k-1)=icont(1,k)
151                 icont(2,k-1)=icont(2,k)
152                 econt(k-1)=econt(k) 
153               enddo
154               i=i-1
155               ncont=ncont-1
156 c              write (iout,*) "ncont",ncont
157 c              do k=1,ncont
158 c                write (iout,*) icont(1,k),icont(2,k)
159 c              enddo
160               goto 20
161             else if (econt(j).gt.ene .and. ic2.ne.ic1+2) 
162      &      then
163               if (ic1.eq.icont(1,j)) then
164                 do k=1,ncont
165                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
166      &               .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. 
167      &               econt(k).lt.econt(i) ) goto 21 
168                 enddo
169               else if (ic2.eq.icont(2,j) ) then
170                 do k=1,ncont
171                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
172      &               .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. 
173      &               econt(k).lt.econt(i) ) goto 21 
174                 enddo
175               endif
176 c Remove jth contact
177               do k=j+1,ncont
178                 icont(1,k-1)=icont(1,k)
179                 icont(2,k-1)=icont(2,k)
180                 econt(k-1)=econt(k) 
181               enddo
182               ncont=ncont-1
183 c              write (iout,*) "ncont",ncont
184 c              do k=1,ncont
185 c                write (iout,*) icont(1,k),icont(2,k)
186 c              enddo
187               j=j-1
188             endif   
189           endif
190    21     continue
191         enddo
192    20   continue
193       enddo
194       if (lprint) then
195         write (iout,*)
196         write (iout,*) 'Electrostatic contacts after pruning: '
197         do i=1,ncont
198           i1=icont(1,i)
199           i2=icont(2,i)
200           it1=itype(i1)
201           it2=itype(i2)
202           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
203      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
204         enddo
205       endif
206       return
207       end