1 subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
3 C This subroutine calculates the average interaction energy and its gradient
4 C in the virtual-bond vectors between non-adjacent peptide groups, based on
5 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
6 C The potential depends both on the distance of peptide-group centers and on
7 C the orientation of the CA-CA virtual bonds.
9 implicit real*8 (a-h,o-z)
11 include 'COMMON.CONTROL'
12 include 'COMMON.IOUNITS'
15 include 'COMMON.LOCAL'
16 include 'COMMON.CHAIN'
17 include 'COMMON.DERIV'
18 include 'COMMON.INTERACT'
19 include 'COMMON.CONTACTS'
20 include 'COMMON.TORSION'
21 include 'COMMON.VECTORS'
22 include 'COMMON.FFIELD'
23 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
24 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
25 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
26 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
27 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
28 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
30 double precision scal_el /1.0d0/
32 double precision scal_el /0.5d0/
35 C 13-go grudnia roku pamietnego...
36 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
39 cd write(iout,*) 'In EELEC'
41 cd write(iout,*) 'Type',i
42 cd write(iout,*) 'B1',B1(:,i)
43 cd write(iout,*) 'B2',B2(:,i)
44 cd write(iout,*) 'CC',CC(:,:,i)
45 cd write(iout,*) 'DD',DD(:,:,i)
46 cd write(iout,*) 'EE',EE(:,:,i)
50 if (icheckgrad.eq.1) then
52 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
54 dc_norm(k,i)=dc(k,i)*fac
56 c write (iout,*) 'i',i,' fac',fac
59 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
60 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
61 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
66 cd write (iout,*) 'i=',i
68 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
71 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
72 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
85 cd print '(a)','Enter EELEC'
86 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
98 xmedi=c(1,i)+0.5d0*dxi
99 ymedi=c(2,i)+0.5d0*dyi
100 zmedi=c(3,i)+0.5d0*dzi
102 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
103 do j=ielstart(i),ielend(i)
107 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
110 ael6i=ael6(iteli,itelj)
111 ael3i=ael3(iteli,itelj)
112 C Diagnostics only!!!
121 dx_normj=dc_norm(1,j)
122 dy_normj=dc_norm(2,j)
123 dz_normj=dc_norm(3,j)
124 xj=c(1,j)+0.5D0*dxj-xmedi
125 yj=c(2,j)+0.5D0*dyj-ymedi
126 zj=c(3,j)+0.5D0*dzj-zmedi
127 rij=xj*xj+yj*yj+zj*zj
133 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
134 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
135 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
136 fac=cosa-3.0D0*cosb*cosg
138 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
139 if (j.eq.i+2) ev1=scal_el*ev1
144 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
147 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
148 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
151 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
152 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
153 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
154 cd & xmedi,ymedi,zmedi,xj,yj,zj
157 write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
158 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
162 C Calculate contributions to the Cartesian gradient.
165 facvdw=-6*rrmij*(ev1+evdwij)
166 facel=-3*rrmij*(el1+eesij)
172 * Radial derivatives. First process both termini of the fragment (i,j)
179 gelc(k,i)=gelc(k,i)+ghalf
180 gelc(k,j)=gelc(k,j)+ghalf
183 * Loop over residues i+1 thru j-1.
187 caug8 gelc(l,k)=gelc(l,k)+ggg(l)
195 gvdwpp(k,i)=gvdwpp(k,i)+ghalf
196 gvdwpp(k,j)=gvdwpp(k,j)+ghalf
199 * Loop over residues i+1 thru j-1.
203 caug8 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
210 fac=-3*rrmij*(facvdw+facvdw+facel)
215 * Radial derivatives. First process both termini of the fragment (i,j)
222 gelc(k,i)=gelc(k,i)+ghalf
223 gelc(k,j)=gelc(k,j)+ghalf
226 * Loop over residues i+1 thru j-1.
230 caug8 gelc(l,k)=gelc(l,k)+ggg(l)
237 ecosa=2.0D0*fac3*fac1+fac4
240 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
241 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
243 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
244 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
246 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
247 cd & (dcosg(k),k=1,3)
249 ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
253 gelc(k,i)=gelc(k,i)+ghalf
254 & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
255 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
256 gelc(k,j)=gelc(k,j)+ghalf
257 & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
258 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
262 caug8 gelc(l,k)=gelc(l,k)+ggg(l)
267 num_cont_hb(i)=num_conti
269 c write (iout,*) "Number of loop steps in EELEC:",ind
271 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
272 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
274 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
275 ccc eel_loc=eel_loc+eello_turn3