added source code
[unres.git] / source / unres / src_MD / src / old_F / eelec.F
1       subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
2 C
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.
8
9       implicit real*8 (a-h,o-z)
10       include 'DIMENSIONS'
11       include 'COMMON.CONTROL'
12       include 'COMMON.IOUNITS'
13       include 'COMMON.GEO'
14       include 'COMMON.VAR'
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
29 #ifdef MOMENT
30       double precision scal_el /1.0d0/
31 #else
32       double precision scal_el /0.5d0/
33 #endif
34 C 12/13/98 
35 C 13-go grudnia roku pamietnego... 
36       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
37      &                   0.0d0,1.0d0,0.0d0,
38      &                   0.0d0,0.0d0,1.0d0/
39 cd      write(iout,*) 'In EELEC'
40 cd      do i=1,nloctyp
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)
47 cd      enddo
48 cd      call check_vecgrad
49 cd      stop
50       if (icheckgrad.eq.1) then
51         do i=1,nres-1
52           fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
53           do k=1,3
54             dc_norm(k,i)=dc(k,i)*fac
55           enddo
56 c          write (iout,*) 'i',i,' fac',fac
57         enddo
58       endif
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
62 c        call vec_and_deriv
63         call set_matrices
64       endif
65 cd      do i=1,nres-1
66 cd        write (iout,*) 'i=',i
67 cd        do k=1,3
68 cd        write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
69 cd        enddo
70 cd        do k=1,3
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)
73 cd        enddo
74 cd      enddo
75       num_conti_hb=0
76       ees=0.0D0
77       evdw1=0.0D0
78       eel_loc=0.0d0 
79       eello_turn3=0.0d0
80       eello_turn4=0.0d0
81       ind=0
82       do i=1,nres
83         num_cont_hb(i)=0
84       enddo
85 cd      print '(a)','Enter EELEC'
86 cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
87       do i=1,nres
88         gel_loc_loc(i)=0.0d0
89         gcorr_loc(i)=0.0d0
90       enddo
91       do i=iatel_s,iatel_e
92         dxi=dc(1,i)
93         dyi=dc(2,i)
94         dzi=dc(3,i)
95         dx_normi=dc_norm(1,i)
96         dy_normi=dc_norm(2,i)
97         dz_normi=dc_norm(3,i)
98         xmedi=c(1,i)+0.5d0*dxi
99         ymedi=c(2,i)+0.5d0*dyi
100         zmedi=c(3,i)+0.5d0*dzi
101         num_conti=0
102 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
103         do j=ielstart(i),ielend(i)
104           ind=ind+1
105           iteli=itel(i)
106           itelj=itel(j)
107           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
108           aaa=app(iteli,itelj)
109           bbb=bpp(iteli,itelj)
110           ael6i=ael6(iteli,itelj)
111           ael3i=ael3(iteli,itelj) 
112 C Diagnostics only!!!
113 c         aaa=0.0D0
114 c         bbb=0.0D0
115 c         ael6i=0.0D0
116 c         ael3i=0.0D0
117 C End diagnostics
118           dxj=dc(1,j)
119           dyj=dc(2,j)
120           dzj=dc(3,j)
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
128           rrmij=1.0D0/rij
129           rij=dsqrt(rij)
130           rmij=1.0D0/rij
131           r3ij=rrmij*rmij
132           r6ij=r3ij*r3ij  
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
137           ev1=aaa*r6ij*r6ij
138 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
139           if (j.eq.i+2) ev1=scal_el*ev1
140           ev2=bbb*r6ij
141           fac3=ael6i*r6ij
142           fac4=ael3i*r3ij
143           evdwij=ev1+ev2
144           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
145           el2=fac4*fac       
146           eesij=el1+el2
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)
149           ees=ees+eesij
150           evdw1=evdw1+evdwij
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
155
156           if (energy_dec) then 
157               write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
158               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
159           endif
160
161 C
162 C Calculate contributions to the Cartesian gradient.
163 C
164 #ifdef SPLITELE
165           facvdw=-6*rrmij*(ev1+evdwij)
166           facel=-3*rrmij*(el1+eesij)
167           fac1=fac
168           erij(1)=xj*rmij
169           erij(2)=yj*rmij
170           erij(3)=zj*rmij
171 *
172 * Radial derivatives. First process both termini of the fragment (i,j)
173 *
174           ggg(1)=facel*xj
175           ggg(2)=facel*yj
176           ggg(3)=facel*zj
177           do k=1,3
178             ghalf=0.5D0*ggg(k)
179             gelc(k,i)=gelc(k,i)+ghalf
180             gelc(k,j)=gelc(k,j)+ghalf
181           enddo
182 *
183 * Loop over residues i+1 thru j-1.
184 *
185 caug8          do k=i+1,j-1
186 caug8            do l=1,3
187 caug8              gelc(l,k)=gelc(l,k)+ggg(l)
188 caug8            enddo
189 caug8          enddo
190           ggg(1)=facvdw*xj
191           ggg(2)=facvdw*yj
192           ggg(3)=facvdw*zj
193           do k=1,3
194             ghalf=0.5D0*ggg(k)
195             gvdwpp(k,i)=gvdwpp(k,i)+ghalf
196             gvdwpp(k,j)=gvdwpp(k,j)+ghalf
197           enddo
198 *
199 * Loop over residues i+1 thru j-1.
200 *
201 caug8          do k=i+1,j-1
202 caug8            do l=1,3
203 caug8              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
204 caug8            enddo
205 caug8          enddo
206 #else
207           facvdw=ev1+evdwij 
208           facel=el1+eesij  
209           fac1=fac
210           fac=-3*rrmij*(facvdw+facvdw+facel)
211           erij(1)=xj*rmij
212           erij(2)=yj*rmij
213           erij(3)=zj*rmij
214 *
215 * Radial derivatives. First process both termini of the fragment (i,j)
216
217           ggg(1)=fac*xj
218           ggg(2)=fac*yj
219           ggg(3)=fac*zj
220           do k=1,3
221             ghalf=0.5D0*ggg(k)
222             gelc(k,i)=gelc(k,i)+ghalf
223             gelc(k,j)=gelc(k,j)+ghalf
224           enddo
225 *
226 * Loop over residues i+1 thru j-1.
227 *
228 caug8          do k=i+1,j-1
229 caug8            do l=1,3
230 caug8              gelc(l,k)=gelc(l,k)+ggg(l)
231 caug8            enddo
232 caug8          enddo
233 #endif
234 *
235 * Angular part
236 *          
237           ecosa=2.0D0*fac3*fac1+fac4
238           fac4=-3.0D0*fac4
239           fac3=-6.0D0*fac3
240           ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
241           ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
242           do k=1,3
243             dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
244             dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
245           enddo
246 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
247 cd   &          (dcosg(k),k=1,3)
248           do k=1,3
249             ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k) 
250           enddo
251           do k=1,3
252             ghalf=0.5D0*ggg(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)
259           enddo
260 caug8          do k=i+1,j-1
261 caug8            do l=1,3
262 caug8              gelc(l,k)=gelc(l,k)+ggg(l)
263 caug8            enddo
264 caug8          enddo
265
266         enddo ! j
267         num_cont_hb(i)=num_conti
268       enddo   ! i
269 c      write (iout,*) "Number of loop steps in EELEC:",ind
270 cd      do i=1,nres
271 cd        write (iout,'(i3,3f10.5,5x,3f10.5)') 
272 cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
273 cd      enddo
274 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
275 ccc      eel_loc=eel_loc+eello_turn3
276       return
277       end
278