1 subroutine elecont(lprint,ncont,icont,ist,ien)
4 include 'DIMENSIONS.ZSCOPT'
5 include 'DIMENSIONS.COMPAR'
6 include 'COMMON.IOUNITS'
8 include 'COMMON.INTERACT'
9 include 'COMMON.FFIELD'
10 include 'COMMON.NAMES'
11 include 'COMMON.LOCAL'
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, rij,zj_temp,xj_temp,yj_temp,
18 & sscale,sscagrad,dist_temp,xj_safe,yj_safe,zj_safe,dist_init
19 double precision elpp6c(2,2),elpp3c(2,2),ael6c(2,2),ael3c(2,2),
21 double precision elcutoff,elecutoff_14
22 integer ncont,icont(2,maxcont),xshift,yshift,zshift,isubchap
23 double precision econt(maxcont)
25 * Load the constants of peptide bond - peptide bond interactions.
26 * Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
27 * proline) - determined by averaging ECEPP energy.
31 c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
32 c data rpp / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
33 data elpp6c /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
34 data elpp3c / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
35 data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
38 if (lprint) write (iout,'(a)')
39 & "Constants of electrostatic interaction energy expression."
43 appc(i,j)=epp(i,j)*rri*rri
44 bppc(i,j)=-2.0*epp(i,j)*rri
45 ael6c(i,j)=elpp6c(i,j)*4.2**6
46 ael3c(i,j)=elpp3c(i,j)*4.2**3
48 & write (iout,'(2i2,4e15.4)') i,j,appc(i,j),bppc(i,j),ael6c(i,j),
63 xmedi=mod(xmedi,boxxsize)
64 if (xmedi.lt.0) xmedi=xmedi+boxxsize
65 ymedi=mod(ymedi,boxysize)
66 if (ymedi.lt.0) ymedi=ymedi+boxysize
67 zmedi=mod(zmedi,boxzsize)
68 if (zmedi.lt.0) zmedi=zmedi+boxzsize
73 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
74 if (iteli.eq.2 .and. itelj.eq.2
75 & .or.iteli.eq.0 .or.itelj.eq.0) goto 4
78 ael6i=ael6c(iteli,itelj)
79 ael3i=ael3c(iteli,itelj)
87 if (xj.lt.0) xj=xj+boxxsize
89 if (yj.lt.0) yj=yj+boxysize
91 if (zj.lt.0) zj=zj+boxzsize
92 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
100 xj=xj_safe+xshift*boxxsize
101 yj=yj_safe+yshift*boxysize
102 zj=zj_safe+zshift*boxzsize
103 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
104 if(dist_temp.lt.dist_init) then
114 if (isubchap.eq.1) then
123 rij=xj*xj+yj*yj+zj*zj
124 sss=sscale(sqrt(rij))
125 sssgrad=sscagrad(sqrt(rij))
126 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
131 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
132 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
133 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
134 fac=cosa-3.0*cosb*cosg
140 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
143 if (j.gt.i+2 .and. eesij.le.elcutoff .or.
144 & j.eq.i+2 .and. eesij.le.elecutoff_14) then
155 write (iout,*) 'Total average electrostatic energy: ',ees
156 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
158 write (iout,*) 'Electrostatic contacts before pruning: '
164 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
165 & i,restyp(it1),i1,restyp(it2),i2,econt(i)
168 c For given residues keep only the contacts with the greatest energy.
170 do while (i.lt.ncont)
176 do while (j.lt.ncont)
178 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
179 & ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
180 c write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
181 c & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
182 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
183 if (ic1.eq.icont(1,j)) then
185 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
186 & .and. iabs(icont(1,k)-ic1).le.2 .and.
187 & econt(k).lt.econt(j) ) goto 21
189 else if (ic2.eq.icont(2,j) ) then
191 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
192 & .and. iabs(icont(2,k)-ic2).le.2 .and.
193 & econt(k).lt.econt(j) ) goto 21
198 icont(1,k-1)=icont(1,k)
199 icont(2,k-1)=icont(2,k)
204 c write (iout,*) "ncont",ncont
206 c write (iout,*) icont(1,k),icont(2,k)
209 else if (econt(j).gt.ene .and. ic2.ne.ic1+2)
211 if (ic1.eq.icont(1,j)) then
213 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
214 & .and. iabs(icont(1,k)-icont(1,j)).le.2 .and.
215 & econt(k).lt.econt(i) ) goto 21
217 else if (ic2.eq.icont(2,j) ) then
219 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
220 & .and. iabs(icont(2,k)-icont(2,j)).le.2 .and.
221 & econt(k).lt.econt(i) ) goto 21
226 icont(1,k-1)=icont(1,k)
227 icont(2,k-1)=icont(2,k)
231 c write (iout,*) "ncont",ncont
233 c write (iout,*) icont(1,k),icont(2,k)
244 write (iout,*) 'Electrostatic contacts after pruning: '
250 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
251 & i,restyp(it1),i1,restyp(it2),i2,econt(i)