c0b4c84f9f32c213d95cfec2ec1b9f5124867d3e
[unres.git] / energy_p_new-sep_barrier.F
1 C-----------------------------------------------------------------------
2       double precision function sscalelip(r)
3       double precision r,gamm
4       include "COMMON.SPLITELE"
5 C      if(r.lt.r_cut-rlamb) then
6 C        sscale=1.0d0
7 C      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
8 C        gamm=(r-(r_cut-rlamb))/rlamb
9         sscalelip=1.0d0+r*r*(2*r-3.0d0)
10 C      else
11 C        sscale=0d0
12 C      endif
13       return
14       end
15 C-----------------------------------------------------------------------
16       double precision function sscagradlip(r)
17       double precision r,gamm
18       include "COMMON.SPLITELE"
19 C     if(r.lt.r_cut-rlamb) then
20 C        sscagrad=0.0d0
21 C      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
22 C        gamm=(r-(r_cut-rlamb))/rlamb
23         sscagradlip=r*(6*r-6.0d0)
24 C      else
25 C        sscagrad=0.0d0
26 C      endif
27       return
28       end
29
30 C-----------------------------------------------------------------------
31       double precision function sscale(r)
32       double precision r,gamm
33       include "COMMON.SPLITELE"
34       if(r.lt.r_cut-rlamb) then
35         sscale=1.0d0
36       else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
37         gamm=(r-(r_cut-rlamb))/rlamb
38         sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
39       else
40         sscale=0d0
41       endif
42       return
43       end
44 C-----------------------------------------------------------------------
45 C-----------------------------------------------------------------------
46       double precision function sscagrad(r)
47       double precision r,gamm
48       include "COMMON.SPLITELE"
49       if(r.lt.r_cut-rlamb) then
50         sscagrad=0.0d0
51       else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
52         gamm=(r-(r_cut-rlamb))/rlamb
53         sscagrad=gamm*(6*gamm-6.0d0)/rlamb
54       else
55         sscagrad=0.0d0
56       endif
57       return
58       end
59 C-----------------------------------------------------------------------
60       subroutine elj_long(evdw)
61 C
62 C This subroutine calculates the interaction energy of nonbonded side chains
63 C assuming the LJ potential of interaction.
64 C
65       implicit real*8 (a-h,o-z)
66       include 'DIMENSIONS'
67       parameter (accur=1.0d-10)
68       include 'COMMON.GEO'
69       include 'COMMON.VAR'
70       include 'COMMON.LOCAL'
71       include 'COMMON.CHAIN'
72       include 'COMMON.DERIV'
73       include 'COMMON.INTERACT'
74       include 'COMMON.TORSION'
75       include 'COMMON.SBRIDGE'
76       include 'COMMON.NAMES'
77       include 'COMMON.IOUNITS'
78       include 'COMMON.CONTACTS'
79       dimension gg(3)
80 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
81       evdw=0.0D0
82       do i=iatsc_s,iatsc_e
83         itypi=itype(i)
84         if (itypi.eq.ntyp1) cycle
85         itypi1=itype(i+1)
86         xi=c(1,nres+i)
87         yi=c(2,nres+i)
88         zi=c(3,nres+i)
89 C
90 C Calculate SC interaction energy.
91 C
92         do iint=1,nint_gr(i)
93 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
94 cd   &                  'iend=',iend(i,iint)
95           do j=istart(i,iint),iend(i,iint)
96             itypj=itype(j)
97             if (itypj.eq.ntyp1) cycle
98             xj=c(1,nres+j)-xi
99             yj=c(2,nres+j)-yi
100             zj=c(3,nres+j)-zi
101             rij=xj*xj+yj*yj+zj*zj
102             sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
103             if (sss.lt.1.0d0) then
104               rrij=1.0D0/rij
105               eps0ij=eps(itypi,itypj)
106               fac=rrij**expon2
107               e1=fac*fac*aa
108               e2=fac*bb
109               evdwij=e1+e2
110               evdw=evdw+(1.0d0-sss)*evdwij
111
112 C Calculate the components of the gradient in DC and X
113 C
114               fac=-rrij*(e1+evdwij)*(1.0d0-sss)
115               gg(1)=xj*fac
116               gg(2)=yj*fac
117               gg(3)=zj*fac
118               do k=1,3
119                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
120                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
121                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
122                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
123               enddo
124             endif
125           enddo      ! j
126         enddo        ! iint
127       enddo          ! i
128       do i=1,nct
129         do j=1,3
130           gvdwc(j,i)=expon*gvdwc(j,i)
131           gvdwx(j,i)=expon*gvdwx(j,i)
132         enddo
133       enddo
134 C******************************************************************************
135 C
136 C                              N O T E !!!
137 C
138 C To save time, the factor of EXPON has been extracted from ALL components
139 C of GVDWC and GRADX. Remember to multiply them by this factor before further 
140 C use!
141 C
142 C******************************************************************************
143       return
144       end
145 C-----------------------------------------------------------------------
146       subroutine elj_short(evdw)
147 C
148 C This subroutine calculates the interaction energy of nonbonded side chains
149 C assuming the LJ potential of interaction.
150 C
151       implicit real*8 (a-h,o-z)
152       include 'DIMENSIONS'
153       parameter (accur=1.0d-10)
154       include 'COMMON.GEO'
155       include 'COMMON.VAR'
156       include 'COMMON.LOCAL'
157       include 'COMMON.CHAIN'
158       include 'COMMON.DERIV'
159       include 'COMMON.INTERACT'
160       include 'COMMON.TORSION'
161       include 'COMMON.SBRIDGE'
162       include 'COMMON.NAMES'
163       include 'COMMON.IOUNITS'
164       include 'COMMON.CONTACTS'
165       dimension gg(3)
166 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
167       evdw=0.0D0
168       do i=iatsc_s,iatsc_e
169         itypi=itype(i)
170         if (itypi.eq.ntyp1) cycle
171         itypi1=itype(i+1)
172         xi=c(1,nres+i)
173         yi=c(2,nres+i)
174         zi=c(3,nres+i)
175 C Change 12/1/95
176         num_conti=0
177 C
178 C Calculate SC interaction energy.
179 C
180         do iint=1,nint_gr(i)
181 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
182 cd   &                  'iend=',iend(i,iint)
183           do j=istart(i,iint),iend(i,iint)
184             itypj=itype(j)
185             if (itypj.eq.ntyp1) cycle
186             xj=c(1,nres+j)-xi
187             yj=c(2,nres+j)-yi
188             zj=c(3,nres+j)-zi
189 C Change 12/1/95 to calculate four-body interactions
190             rij=xj*xj+yj*yj+zj*zj
191             sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
192             if (sss.gt.0.0d0) then
193               rrij=1.0D0/rij
194               eps0ij=eps(itypi,itypj)
195               fac=rrij**expon2
196               e1=fac*fac*aa
197               e2=fac*bb
198               evdwij=e1+e2
199               evdw=evdw+sss*evdwij
200
201 C Calculate the components of the gradient in DC and X
202 C
203               fac=-rrij*(e1+evdwij)*sss
204               gg(1)=xj*fac
205               gg(2)=yj*fac
206               gg(3)=zj*fac
207               do k=1,3
208                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
209                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
210                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
211                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
212               enddo
213             endif
214           enddo      ! j
215         enddo        ! iint
216       enddo          ! i
217       do i=1,nct
218         do j=1,3
219           gvdwc(j,i)=expon*gvdwc(j,i)
220           gvdwx(j,i)=expon*gvdwx(j,i)
221         enddo
222       enddo
223 C******************************************************************************
224 C
225 C                              N O T E !!!
226 C
227 C To save time, the factor of EXPON has been extracted from ALL components
228 C of GVDWC and GRADX. Remember to multiply them by this factor before further 
229 C use!
230 C
231 C******************************************************************************
232       return
233       end
234 C-----------------------------------------------------------------------------
235       subroutine eljk_long(evdw)
236 C
237 C This subroutine calculates the interaction energy of nonbonded side chains
238 C assuming the LJK potential of interaction.
239 C
240       implicit real*8 (a-h,o-z)
241       include 'DIMENSIONS'
242       include 'COMMON.GEO'
243       include 'COMMON.VAR'
244       include 'COMMON.LOCAL'
245       include 'COMMON.CHAIN'
246       include 'COMMON.DERIV'
247       include 'COMMON.INTERACT'
248       include 'COMMON.IOUNITS'
249       include 'COMMON.NAMES'
250       dimension gg(3)
251       logical scheck
252 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
253       evdw=0.0D0
254       do i=iatsc_s,iatsc_e
255         itypi=itype(i)
256         if (itypi.eq.ntyp1) cycle
257         itypi1=itype(i+1)
258         xi=c(1,nres+i)
259         yi=c(2,nres+i)
260         zi=c(3,nres+i)
261 C
262 C Calculate SC interaction energy.
263 C
264         do iint=1,nint_gr(i)
265           do j=istart(i,iint),iend(i,iint)
266             itypj=itype(j)
267             if (itypj.eq.ntyp1) cycle
268             xj=c(1,nres+j)-xi
269             yj=c(2,nres+j)-yi
270             zj=c(3,nres+j)-zi
271             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
272             fac_augm=rrij**expon
273             e_augm=augm(itypi,itypj)*fac_augm
274             r_inv_ij=dsqrt(rrij)
275             rij=1.0D0/r_inv_ij 
276             sss=sscale(rij/sigma(itypi,itypj))
277             if (sss.lt.1.0d0) then
278               r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
279               fac=r_shift_inv**expon
280               e1=fac*fac*aa
281               e2=fac*bb
282               evdwij=e_augm+e1+e2
283 cd            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
284 cd            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
285 cd            write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
286 cd   &          restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
287 cd   &          bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
288 cd   &          sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
289 cd   &          (c(k,i),k=1,3),(c(k,j),k=1,3)
290               evdw=evdw+(1.0d0-sss)*evdwij
291
292 C Calculate the components of the gradient in DC and X
293 C
294               fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
295               fac=fac*(1.0d0-sss)
296               gg(1)=xj*fac
297               gg(2)=yj*fac
298               gg(3)=zj*fac
299               do k=1,3
300                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
301                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
302                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
303                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
304               enddo
305             endif
306           enddo      ! j
307         enddo        ! iint
308       enddo          ! i
309       do i=1,nct
310         do j=1,3
311           gvdwc(j,i)=expon*gvdwc(j,i)
312           gvdwx(j,i)=expon*gvdwx(j,i)
313         enddo
314       enddo
315       return
316       end
317 C-----------------------------------------------------------------------------
318       subroutine eljk_short(evdw)
319 C
320 C This subroutine calculates the interaction energy of nonbonded side chains
321 C assuming the LJK potential of interaction.
322 C
323       implicit real*8 (a-h,o-z)
324       include 'DIMENSIONS'
325       include 'COMMON.GEO'
326       include 'COMMON.VAR'
327       include 'COMMON.LOCAL'
328       include 'COMMON.CHAIN'
329       include 'COMMON.DERIV'
330       include 'COMMON.INTERACT'
331       include 'COMMON.IOUNITS'
332       include 'COMMON.NAMES'
333       dimension gg(3)
334       logical scheck
335 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
336       evdw=0.0D0
337       do i=iatsc_s,iatsc_e
338         itypi=itype(i)
339         if (itypi.eq.ntyp1) cycle
340         itypi1=itype(i+1)
341         xi=c(1,nres+i)
342         yi=c(2,nres+i)
343         zi=c(3,nres+i)
344 C
345 C Calculate SC interaction energy.
346 C
347         do iint=1,nint_gr(i)
348           do j=istart(i,iint),iend(i,iint)
349             itypj=itype(j)
350             if (itypj.eq.ntyp1) cycle
351             xj=c(1,nres+j)-xi
352             yj=c(2,nres+j)-yi
353             zj=c(3,nres+j)-zi
354             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
355             fac_augm=rrij**expon
356             e_augm=augm(itypi,itypj)*fac_augm
357             r_inv_ij=dsqrt(rrij)
358             rij=1.0D0/r_inv_ij 
359             sss=sscale(rij/sigma(itypi,itypj))
360             if (sss.gt.0.0d0) then
361               r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
362               fac=r_shift_inv**expon
363               e1=fac*fac*aa
364               e2=fac*bb
365               evdwij=e_augm+e1+e2
366 cd            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
367 cd            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
368 cd            write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
369 cd   &          restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
370 cd   &          bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
371 cd   &          sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
372 cd   &          (c(k,i),k=1,3),(c(k,j),k=1,3)
373               evdw=evdw+sss*evdwij
374
375 C Calculate the components of the gradient in DC and X
376 C
377               fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
378               fac=fac*sss
379               gg(1)=xj*fac
380               gg(2)=yj*fac
381               gg(3)=zj*fac
382               do k=1,3
383                 gvdwx(k,i)=gvdwx(k,i)-gg(k)
384                 gvdwx(k,j)=gvdwx(k,j)+gg(k)
385                 gvdwc(k,i)=gvdwc(k,i)-gg(k)
386                 gvdwc(k,j)=gvdwc(k,j)+gg(k)
387               enddo
388             endif
389           enddo      ! j
390         enddo        ! iint
391       enddo          ! i
392       do i=1,nct
393         do j=1,3
394           gvdwc(j,i)=expon*gvdwc(j,i)
395           gvdwx(j,i)=expon*gvdwx(j,i)
396         enddo
397       enddo
398       return
399       end
400 C-----------------------------------------------------------------------------
401       subroutine ebp_long(evdw)
402 C
403 C This subroutine calculates the interaction energy of nonbonded side chains
404 C assuming the Berne-Pechukas potential of interaction.
405 C
406       implicit real*8 (a-h,o-z)
407       include 'DIMENSIONS'
408       include 'COMMON.GEO'
409       include 'COMMON.VAR'
410       include 'COMMON.LOCAL'
411       include 'COMMON.CHAIN'
412       include 'COMMON.DERIV'
413       include 'COMMON.NAMES'
414       include 'COMMON.INTERACT'
415       include 'COMMON.IOUNITS'
416       include 'COMMON.CALC'
417       common /srutu/ icall
418 c     double precision rrsave(maxdim)
419       logical lprn
420       evdw=0.0D0
421 c     print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
422       evdw=0.0D0
423 c     if (icall.eq.0) then
424 c       lprn=.true.
425 c     else
426         lprn=.false.
427 c     endif
428       ind=0
429       do i=iatsc_s,iatsc_e
430         itypi=itype(i)
431         if (itypi.eq.ntyp1) cycle
432         itypi1=itype(i+1)
433         xi=c(1,nres+i)
434         yi=c(2,nres+i)
435         zi=c(3,nres+i)
436         dxi=dc_norm(1,nres+i)
437         dyi=dc_norm(2,nres+i)
438         dzi=dc_norm(3,nres+i)
439 c        dsci_inv=dsc_inv(itypi)
440         dsci_inv=vbld_inv(i+nres)
441 C
442 C Calculate SC interaction energy.
443 C
444         do iint=1,nint_gr(i)
445           do j=istart(i,iint),iend(i,iint)
446             ind=ind+1
447             itypj=itype(j)
448             if (itypj.eq.ntyp1) cycle
449 c            dscj_inv=dsc_inv(itypj)
450             dscj_inv=vbld_inv(j+nres)
451             chi1=chi(itypi,itypj)
452             chi2=chi(itypj,itypi)
453             chi12=chi1*chi2
454             chip1=chip(itypi)
455             chip2=chip(itypj)
456             chip12=chip1*chip2
457             alf1=alp(itypi)
458             alf2=alp(itypj)
459             alf12=0.5D0*(alf1+alf2)
460             xj=c(1,nres+j)-xi
461             yj=c(2,nres+j)-yi
462             zj=c(3,nres+j)-zi
463             dxj=dc_norm(1,nres+j)
464             dyj=dc_norm(2,nres+j)
465             dzj=dc_norm(3,nres+j)
466             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
467             rij=dsqrt(rrij)
468             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
469
470             if (sss.lt.1.0d0) then
471
472 C Calculate the angle-dependent terms of energy & contributions to derivatives.
473               call sc_angular
474 C Calculate whole angle-dependent part of epsilon and contributions
475 C to its derivatives
476               fac=(rrij*sigsq)**expon2
477               e1=fac*fac*aa
478               e2=fac*bb
479               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
480               eps2der=evdwij*eps3rt
481               eps3der=evdwij*eps2rt
482               evdwij=evdwij*eps2rt*eps3rt
483               evdw=evdw+evdwij*(1.0d0-sss)
484               if (lprn) then
485               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
486               epsi=bb**2/aa
487 cd              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
488 cd     &          restyp(itypi),i,restyp(itypj),j,
489 cd     &          epsi,sigm,chi1,chi2,chip1,chip2,
490 cd     &          eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
491 cd     &          om1,om2,om12,1.0D0/dsqrt(rrij),
492 cd     &          evdwij
493               endif
494 C Calculate gradient components.
495               e1=e1*eps1*eps2rt**2*eps3rt**2
496               fac=-expon*(e1+evdwij)
497               sigder=fac/sigsq
498               fac=rrij*fac
499 C Calculate radial part of the gradient
500               gg(1)=xj*fac
501               gg(2)=yj*fac
502               gg(3)=zj*fac
503 C Calculate the angular part of the gradient and sum add the contributions
504 C to the appropriate components of the Cartesian gradient.
505               call sc_grad_scale(1.0d0-sss)
506             endif
507           enddo      ! j
508         enddo        ! iint
509       enddo          ! i
510 c     stop
511       return
512       end
513 C-----------------------------------------------------------------------------
514       subroutine ebp_short(evdw)
515 C
516 C This subroutine calculates the interaction energy of nonbonded side chains
517 C assuming the Berne-Pechukas potential of interaction.
518 C
519       implicit real*8 (a-h,o-z)
520       include 'DIMENSIONS'
521       include 'COMMON.GEO'
522       include 'COMMON.VAR'
523       include 'COMMON.LOCAL'
524       include 'COMMON.CHAIN'
525       include 'COMMON.DERIV'
526       include 'COMMON.NAMES'
527       include 'COMMON.INTERACT'
528       include 'COMMON.IOUNITS'
529       include 'COMMON.CALC'
530       common /srutu/ icall
531 c     double precision rrsave(maxdim)
532       logical lprn
533       evdw=0.0D0
534 c     print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
535       evdw=0.0D0
536 c     if (icall.eq.0) then
537 c       lprn=.true.
538 c     else
539         lprn=.false.
540 c     endif
541       ind=0
542       do i=iatsc_s,iatsc_e
543         itypi=itype(i)
544         if (itypi.eq.ntyp1) cycle
545         itypi1=itype(i+1)
546         xi=c(1,nres+i)
547         yi=c(2,nres+i)
548         zi=c(3,nres+i)
549         dxi=dc_norm(1,nres+i)
550         dyi=dc_norm(2,nres+i)
551         dzi=dc_norm(3,nres+i)
552 c        dsci_inv=dsc_inv(itypi)
553         dsci_inv=vbld_inv(i+nres)
554 C
555 C Calculate SC interaction energy.
556 C
557         do iint=1,nint_gr(i)
558           do j=istart(i,iint),iend(i,iint)
559             ind=ind+1
560             itypj=itype(j)
561             if (itypj.eq.ntyp1) cycle
562 c            dscj_inv=dsc_inv(itypj)
563             dscj_inv=vbld_inv(j+nres)
564             chi1=chi(itypi,itypj)
565             chi2=chi(itypj,itypi)
566             chi12=chi1*chi2
567             chip1=chip(itypi)
568             chip2=chip(itypj)
569             chip12=chip1*chip2
570             alf1=alp(itypi)
571             alf2=alp(itypj)
572             alf12=0.5D0*(alf1+alf2)
573             xj=c(1,nres+j)-xi
574             yj=c(2,nres+j)-yi
575             zj=c(3,nres+j)-zi
576             dxj=dc_norm(1,nres+j)
577             dyj=dc_norm(2,nres+j)
578             dzj=dc_norm(3,nres+j)
579             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
580             rij=dsqrt(rrij)
581             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
582
583             if (sss.gt.0.0d0) then
584
585 C Calculate the angle-dependent terms of energy & contributions to derivatives.
586               call sc_angular
587 C Calculate whole angle-dependent part of epsilon and contributions
588 C to its derivatives
589               fac=(rrij*sigsq)**expon2
590               e1=fac*fac*aa
591               e2=fac*bb
592               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
593               eps2der=evdwij*eps3rt
594               eps3der=evdwij*eps2rt
595               evdwij=evdwij*eps2rt*eps3rt
596               evdw=evdw+evdwij*sss
597               if (lprn) then
598               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
599               epsi=bb**2/aa
600 cd              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
601 cd     &          restyp(itypi),i,restyp(itypj),j,
602 cd     &          epsi,sigm,chi1,chi2,chip1,chip2,
603 cd     &          eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
604 cd     &          om1,om2,om12,1.0D0/dsqrt(rrij),
605 cd     &          evdwij
606               endif
607 C Calculate gradient components.
608               e1=e1*eps1*eps2rt**2*eps3rt**2
609               fac=-expon*(e1+evdwij)
610               sigder=fac/sigsq
611               fac=rrij*fac
612 C Calculate radial part of the gradient
613               gg(1)=xj*fac
614               gg(2)=yj*fac
615               gg(3)=zj*fac
616 C Calculate the angular part of the gradient and sum add the contributions
617 C to the appropriate components of the Cartesian gradient.
618               call sc_grad_scale(sss)
619             endif
620           enddo      ! j
621         enddo        ! iint
622       enddo          ! i
623 c     stop
624       return
625       end
626 C-----------------------------------------------------------------------------
627       subroutine egb_long(evdw)
628 C
629 C This subroutine calculates the interaction energy of nonbonded side chains
630 C assuming the Gay-Berne potential of interaction.
631 C
632       implicit real*8 (a-h,o-z)
633       include 'DIMENSIONS'
634       include 'COMMON.GEO'
635       include 'COMMON.VAR'
636       include 'COMMON.LOCAL'
637       include 'COMMON.CHAIN'
638       include 'COMMON.DERIV'
639       include 'COMMON.NAMES'
640       include 'COMMON.INTERACT'
641       include 'COMMON.IOUNITS'
642       include 'COMMON.CALC'
643       include 'COMMON.CONTROL'
644       logical lprn
645       evdw=0.0D0
646 ccccc      energy_dec=.false.
647 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
648       evdw=0.0D0
649       lprn=.false.
650 c     if (icall.eq.0) lprn=.false.
651       ind=0
652       do i=iatsc_s,iatsc_e
653         itypi=itype(i)
654         if (itypi.eq.ntyp1) cycle
655         itypi1=itype(i+1)
656         xi=c(1,nres+i)
657         yi=c(2,nres+i)
658         zi=c(3,nres+i)
659           xi=mod(xi,boxxsize)
660           if (xi.lt.0) xi=xi+boxxsize
661           yi=mod(yi,boxysize)
662           if (yi.lt.0) yi=yi+boxysize
663           zi=mod(zi,boxzsize)
664           if (zi.lt.0) zi=zi+boxzsize
665         dxi=dc_norm(1,nres+i)
666         dyi=dc_norm(2,nres+i)
667         dzi=dc_norm(3,nres+i)
668 c        dsci_inv=dsc_inv(itypi)
669         dsci_inv=vbld_inv(i+nres)
670 c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
671 c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
672 C
673 C Calculate SC interaction energy.
674 C
675         do iint=1,nint_gr(i)
676           do j=istart(i,iint),iend(i,iint)
677             ind=ind+1
678             itypj=itype(j)
679             if (itypj.eq.ntyp1) cycle
680 c            dscj_inv=dsc_inv(itypj)
681             dscj_inv=vbld_inv(j+nres)
682 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
683 c     &       1.0d0/vbld(j+nres)
684 c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
685             sig0ij=sigma(itypi,itypj)
686             chi1=chi(itypi,itypj)
687             chi2=chi(itypj,itypi)
688             chi12=chi1*chi2
689             chip1=chip(itypi)
690             chip2=chip(itypj)
691             chip12=chip1*chip2
692             alf1=alp(itypi)
693             alf2=alp(itypj)
694             alf12=0.5D0*(alf1+alf2)
695             xj=c(1,nres+j)
696             yj=c(2,nres+j)
697             zj=c(3,nres+j)
698           xj=mod(xj,boxxsize)
699           if (xj.lt.0) xj=xj+boxxsize
700           yj=mod(yj,boxysize)
701           if (yj.lt.0) yj=yj+boxysize
702           zj=mod(zj,boxzsize)
703           if (zj.lt.0) zj=zj+boxzsize
704        if ((zj.gt.bordlipbot)
705      &.and.(zj.lt.bordliptop)) then
706 C the energy transfer exist
707         if (zj.lt.buflipbot) then
708 C what fraction I am in
709          fracinbuf=1.0d0-
710      &        ((positi-bordlipbot)/lipbufthick)
711 C lipbufthick is thickenes of lipid buffore
712          sslipj=sscalelip(fracinbuf)
713          ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
714         elseif (zi.gt.bufliptop) then
715          fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
716          sslipj=sscalelip(fracinbuf)
717          ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
718         else
719          sslipj=1.0d0
720          ssgradlipj=0.0
721         endif
722        else
723          sslipj=0.0d0
724          ssgradlipj=0.0
725        endif
726       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
727      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
728       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
729      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
730
731       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
732       xj_safe=xj
733       yj_safe=yj
734       zj_safe=zj
735       subchap=0
736       do xshift=-1,1
737       do yshift=-1,1
738       do zshift=-1,1
739           xj=xj_safe+xshift*boxxsize
740           yj=yj_safe+yshift*boxysize
741           zj=zj_safe+zshift*boxzsize
742           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
743           if(dist_temp.lt.dist_init) then
744             dist_init=dist_temp
745             xj_temp=xj
746             yj_temp=yj
747             zj_temp=zj
748             subchap=1
749           endif
750        enddo
751        enddo
752        enddo
753        if (subchap.eq.1) then
754           xj=xj_temp-xi
755           yj=yj_temp-yi
756           zj=zj_temp-zi
757        else
758           xj=xj_safe-xi
759           yj=yj_safe-yi
760           zj=zj_safe-zi
761        endif
762             dxj=dc_norm(1,nres+j)
763             dyj=dc_norm(2,nres+j)
764             dzj=dc_norm(3,nres+j)
765             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
766             rij=dsqrt(rrij)
767             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
768             sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
769             if (sss.lt.1.0d0) then
770
771 C Calculate angle-dependent terms of energy and contributions to their
772 C derivatives.
773               call sc_angular
774               sigsq=1.0D0/sigsq
775               sig=sig0ij*dsqrt(sigsq)
776               rij_shift=1.0D0/rij-sig+sig0ij
777 c for diagnostics; uncomment
778 c              rij_shift=1.2*sig0ij
779 C I hate to put IF's in the loops, but here don't have another choice!!!!
780               if (rij_shift.le.0.0D0) then
781                 evdw=1.0D20
782 cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
783 cd     &          restyp(itypi),i,restyp(itypj),j,
784 cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
785                 return
786               endif
787               sigder=-sig*sigsq
788 c---------------------------------------------------------------
789               rij_shift=1.0D0/rij_shift 
790               fac=rij_shift**expon
791               e1=fac*fac*aa
792               e2=fac*bb
793               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
794               eps2der=evdwij*eps3rt
795               eps3der=evdwij*eps2rt
796 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
797 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
798               evdwij=evdwij*eps2rt*eps3rt
799               evdw=evdw+evdwij*(1.0d0-sss)
800               if (lprn) then
801               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
802               epsi=bb**2/aa
803               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
804      &          restyp(itypi),i,restyp(itypj),j,
805      &          epsi,sigm,chi1,chi2,chip1,chip2,
806      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
807      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
808      &          evdwij
809               endif
810
811               if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
812      &                        'evdw',i,j,evdwij
813
814 C Calculate gradient components.
815               e1=e1*eps1*eps2rt**2*eps3rt**2
816               fac=-expon*(e1+evdwij)*rij_shift
817               sigder=fac*sigder
818               fac=rij*fac
819             fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij
820 c              fac=0.0d0
821 C Calculate the radial part of the gradient
822               gg(1)=xj*fac
823               gg(2)=yj*fac
824               gg(3)=zj*fac
825               gg_lipi(3)=ssgradlipi*evdwij
826               gg_lipj(3)=ssgradlipj*evdwij
827 C Calculate angular part of the gradient.
828               call sc_grad_scale(1.0d0-sss)
829             endif
830           enddo      ! j
831         enddo        ! iint
832       enddo          ! i
833 c      write (iout,*) "Number of loop steps in EGB:",ind
834 cccc      energy_dec=.false.
835       return
836       end
837 C-----------------------------------------------------------------------------
838       subroutine egb_short(evdw)
839 C
840 C This subroutine calculates the interaction energy of nonbonded side chains
841 C assuming the Gay-Berne potential of interaction.
842 C
843       implicit real*8 (a-h,o-z)
844       include 'DIMENSIONS'
845       include 'COMMON.GEO'
846       include 'COMMON.VAR'
847       include 'COMMON.LOCAL'
848       include 'COMMON.CHAIN'
849       include 'COMMON.DERIV'
850       include 'COMMON.NAMES'
851       include 'COMMON.INTERACT'
852       include 'COMMON.IOUNITS'
853       include 'COMMON.CALC'
854       include 'COMMON.CONTROL'
855       logical lprn
856       evdw=0.0D0
857 ccccc      energy_dec=.false.
858 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
859       evdw=0.0D0
860       lprn=.false.
861 c     if (icall.eq.0) lprn=.false.
862       ind=0
863       do i=iatsc_s,iatsc_e
864         itypi=itype(i)
865         if (itypi.eq.ntyp1) cycle
866         itypi1=itype(i+1)
867         xi=c(1,nres+i)
868         yi=c(2,nres+i)
869         zi=c(3,nres+i)
870           xi=mod(xi,boxxsize)
871           if (xi.lt.0) xi=xi+boxxsize
872           yi=mod(yi,boxysize)
873           if (yi.lt.0) yi=yi+boxysize
874           zi=mod(zi,boxzsize)
875           if (zi.lt.0) zi=zi+boxzsize
876         dxi=dc_norm(1,nres+i)
877         dyi=dc_norm(2,nres+i)
878         dzi=dc_norm(3,nres+i)
879 c        dsci_inv=dsc_inv(itypi)
880         dsci_inv=vbld_inv(i+nres)
881 c        write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
882 c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
883 C
884 C Calculate SC interaction energy.
885 C
886         do iint=1,nint_gr(i)
887           do j=istart(i,iint),iend(i,iint)
888             ind=ind+1
889             itypj=itype(j)
890             if (itypj.eq.ntyp1) cycle
891 c            dscj_inv=dsc_inv(itypj)
892             dscj_inv=vbld_inv(j+nres)
893 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
894 c     &       1.0d0/vbld(j+nres)
895 c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
896             sig0ij=sigma(itypi,itypj)
897             chi1=chi(itypi,itypj)
898             chi2=chi(itypj,itypi)
899             chi12=chi1*chi2
900             chip1=chip(itypi)
901             chip2=chip(itypj)
902             chip12=chip1*chip2
903             alf1=alp(itypi)
904             alf2=alp(itypj)
905             alf12=0.5D0*(alf1+alf2)
906             xj=c(1,nres+j)
907             yj=c(2,nres+j)
908             zj=c(3,nres+j)
909           xj=mod(xj,boxxsize)
910           if (xj.lt.0) xj=xj+boxxsize
911           yj=mod(yj,boxysize)
912           if (yj.lt.0) yj=yj+boxysize
913           zj=mod(zj,boxzsize)
914           if (zj.lt.0) zj=zj+boxzsize
915        if ((zj.gt.bordlipbot)
916      &.and.(zj.lt.bordliptop)) then
917 C the energy transfer exist
918         if (zj.lt.buflipbot) then
919 C what fraction I am in
920          fracinbuf=1.0d0-
921      &        ((positi-bordlipbot)/lipbufthick)
922 C lipbufthick is thickenes of lipid buffore
923          sslipj=sscalelip(fracinbuf)
924          ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
925         elseif (zi.gt.bufliptop) then
926          fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
927          sslipj=sscalelip(fracinbuf)
928          ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
929         else
930          sslipj=1.0d0
931          ssgradlipj=0.0
932         endif
933        else
934          sslipj=0.0d0
935          ssgradlipj=0.0
936        endif
937       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
938      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
939       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
940      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
941       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
942       xj_safe=xj
943       yj_safe=yj
944       zj_safe=zj
945       subchap=0
946       do xshift=-1,1
947       do yshift=-1,1
948       do zshift=-1,1
949           xj=xj_safe+xshift*boxxsize
950           yj=yj_safe+yshift*boxysize
951           zj=zj_safe+zshift*boxzsize
952           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
953           if(dist_temp.lt.dist_init) then
954             dist_init=dist_temp
955             xj_temp=xj
956             yj_temp=yj
957             zj_temp=zj
958             subchap=1
959           endif
960        enddo
961        enddo
962        enddo
963        if (subchap.eq.1) then
964           xj=xj_temp-xi
965           yj=yj_temp-yi
966           zj=zj_temp-zi
967        else
968           xj=xj_safe-xi
969           yj=yj_safe-yi
970           zj=zj_safe-zi
971        endif
972             dxj=dc_norm(1,nres+j)
973             dyj=dc_norm(2,nres+j)
974             dzj=dc_norm(3,nres+j)
975             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
976             rij=dsqrt(rrij)
977             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
978             sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
979             if (sss.gt.0.0d0) then
980
981 C Calculate angle-dependent terms of energy and contributions to their
982 C derivatives.
983               call sc_angular
984               sigsq=1.0D0/sigsq
985               sig=sig0ij*dsqrt(sigsq)
986               rij_shift=1.0D0/rij-sig+sig0ij
987 c for diagnostics; uncomment
988 c              rij_shift=1.2*sig0ij
989 C I hate to put IF's in the loops, but here don't have another choice!!!!
990               if (rij_shift.le.0.0D0) then
991                 evdw=1.0D20
992 cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
993 cd     &          restyp(itypi),i,restyp(itypj),j,
994 cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
995                 return
996               endif
997               sigder=-sig*sigsq
998 c---------------------------------------------------------------
999               rij_shift=1.0D0/rij_shift 
1000               fac=rij_shift**expon
1001               e1=fac*fac*aa
1002               e2=fac*bb
1003               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1004               eps2der=evdwij*eps3rt
1005               eps3der=evdwij*eps2rt
1006 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1007 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1008               evdwij=evdwij*eps2rt*eps3rt
1009               evdw=evdw+evdwij*sss
1010               if (lprn) then
1011               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1012               epsi=bb**2/aa
1013               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1014      &          restyp(itypi),i,restyp(itypj),j,
1015      &          epsi,sigm,chi1,chi2,chip1,chip2,
1016      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1017      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1018      &          evdwij
1019               endif
1020
1021               if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
1022      &                        'evdw',i,j,evdwij
1023
1024 C Calculate gradient components.
1025               e1=e1*eps1*eps2rt**2*eps3rt**2
1026               fac=-expon*(e1+evdwij)*rij_shift
1027               sigder=fac*sigder
1028               fac=rij*fac
1029             fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij
1030 c              fac=0.0d0
1031 C Calculate the radial part of the gradient
1032               gg(1)=xj*fac
1033               gg(2)=yj*fac
1034               gg(3)=zj*fac
1035               gg_lipi(3)=ssgradlipi*evdwij
1036               gg_lipj(3)=ssgradlipj*evdwij
1037 C Calculate angular part of the gradient.
1038               call sc_grad_scale(sss)
1039             endif
1040           enddo      ! j
1041         enddo        ! iint
1042       enddo          ! i
1043 c      write (iout,*) "Number of loop steps in EGB:",ind
1044 cccc      energy_dec=.false.
1045       return
1046       end
1047 C-----------------------------------------------------------------------------
1048       subroutine egbv_long(evdw)
1049 C
1050 C This subroutine calculates the interaction energy of nonbonded side chains
1051 C assuming the Gay-Berne-Vorobjev potential of interaction.
1052 C
1053       implicit real*8 (a-h,o-z)
1054       include 'DIMENSIONS'
1055       include 'COMMON.GEO'
1056       include 'COMMON.VAR'
1057       include 'COMMON.LOCAL'
1058       include 'COMMON.CHAIN'
1059       include 'COMMON.DERIV'
1060       include 'COMMON.NAMES'
1061       include 'COMMON.INTERACT'
1062       include 'COMMON.IOUNITS'
1063       include 'COMMON.CALC'
1064       common /srutu/ icall
1065       logical lprn
1066       evdw=0.0D0
1067 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1068       evdw=0.0D0
1069       lprn=.false.
1070 c     if (icall.eq.0) lprn=.true.
1071       ind=0
1072       do i=iatsc_s,iatsc_e
1073         itypi=itype(i)
1074         if (itypi.eq.ntyp1) cycle
1075         itypi1=itype(i+1)
1076         xi=c(1,nres+i)
1077         yi=c(2,nres+i)
1078         zi=c(3,nres+i)
1079         dxi=dc_norm(1,nres+i)
1080         dyi=dc_norm(2,nres+i)
1081         dzi=dc_norm(3,nres+i)
1082 c        dsci_inv=dsc_inv(itypi)
1083         dsci_inv=vbld_inv(i+nres)
1084 C
1085 C Calculate SC interaction energy.
1086 C
1087         do iint=1,nint_gr(i)
1088           do j=istart(i,iint),iend(i,iint)
1089             ind=ind+1
1090             itypj=itype(j)
1091             if (itypj.eq.ntyp1) cycle
1092 c            dscj_inv=dsc_inv(itypj)
1093             dscj_inv=vbld_inv(j+nres)
1094             sig0ij=sigma(itypi,itypj)
1095             r0ij=r0(itypi,itypj)
1096             chi1=chi(itypi,itypj)
1097             chi2=chi(itypj,itypi)
1098             chi12=chi1*chi2
1099             chip1=chip(itypi)
1100             chip2=chip(itypj)
1101             chip12=chip1*chip2
1102             alf1=alp(itypi)
1103             alf2=alp(itypj)
1104             alf12=0.5D0*(alf1+alf2)
1105             xj=c(1,nres+j)-xi
1106             yj=c(2,nres+j)-yi
1107             zj=c(3,nres+j)-zi
1108             dxj=dc_norm(1,nres+j)
1109             dyj=dc_norm(2,nres+j)
1110             dzj=dc_norm(3,nres+j)
1111             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1112             rij=dsqrt(rrij)
1113
1114             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1115
1116             if (sss.lt.1.0d0) then
1117
1118 C Calculate angle-dependent terms of energy and contributions to their
1119 C derivatives.
1120               call sc_angular
1121               sigsq=1.0D0/sigsq
1122               sig=sig0ij*dsqrt(sigsq)
1123               rij_shift=1.0D0/rij-sig+r0ij
1124 C I hate to put IF's in the loops, but here don't have another choice!!!!
1125               if (rij_shift.le.0.0D0) then
1126                 evdw=1.0D20
1127                 return
1128               endif
1129               sigder=-sig*sigsq
1130 c---------------------------------------------------------------
1131               rij_shift=1.0D0/rij_shift 
1132               fac=rij_shift**expon
1133               e1=fac*fac*aa
1134               e2=fac*bb
1135               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1136               eps2der=evdwij*eps3rt
1137               eps3der=evdwij*eps2rt
1138               fac_augm=rrij**expon
1139               e_augm=augm(itypi,itypj)*fac_augm
1140               evdwij=evdwij*eps2rt*eps3rt
1141               evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
1142               if (lprn) then
1143               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1144               epsi=bb**2/aa
1145               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1146      &          restyp(itypi),i,restyp(itypj),j,
1147      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1148      &          chi1,chi2,chip1,chip2,
1149      &          eps1,eps2rt**2,eps3rt**2,
1150      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1151      &          evdwij+e_augm
1152               endif
1153 C Calculate gradient components.
1154               e1=e1*eps1*eps2rt**2*eps3rt**2
1155               fac=-expon*(e1+evdwij)*rij_shift
1156               sigder=fac*sigder
1157               fac=rij*fac-2*expon*rrij*e_augm
1158 C Calculate the radial part of the gradient
1159               gg(1)=xj*fac
1160               gg(2)=yj*fac
1161               gg(3)=zj*fac
1162 C Calculate angular part of the gradient.
1163               call sc_grad_scale(1.0d0-sss)
1164             endif
1165           enddo      ! j
1166         enddo        ! iint
1167       enddo          ! i
1168       end
1169 C-----------------------------------------------------------------------------
1170       subroutine egbv_short(evdw)
1171 C
1172 C This subroutine calculates the interaction energy of nonbonded side chains
1173 C assuming the Gay-Berne-Vorobjev potential of interaction.
1174 C
1175       implicit real*8 (a-h,o-z)
1176       include 'DIMENSIONS'
1177       include 'COMMON.GEO'
1178       include 'COMMON.VAR'
1179       include 'COMMON.LOCAL'
1180       include 'COMMON.CHAIN'
1181       include 'COMMON.DERIV'
1182       include 'COMMON.NAMES'
1183       include 'COMMON.INTERACT'
1184       include 'COMMON.IOUNITS'
1185       include 'COMMON.CALC'
1186       common /srutu/ icall
1187       logical lprn
1188       evdw=0.0D0
1189 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1190       evdw=0.0D0
1191       lprn=.false.
1192 c     if (icall.eq.0) lprn=.true.
1193       ind=0
1194       do i=iatsc_s,iatsc_e
1195         itypi=itype(i)
1196         if (itypi.eq.ntyp1) cycle
1197         itypi1=itype(i+1)
1198         xi=c(1,nres+i)
1199         yi=c(2,nres+i)
1200         zi=c(3,nres+i)
1201         dxi=dc_norm(1,nres+i)
1202         dyi=dc_norm(2,nres+i)
1203         dzi=dc_norm(3,nres+i)
1204 c        dsci_inv=dsc_inv(itypi)
1205         dsci_inv=vbld_inv(i+nres)
1206 C
1207 C Calculate SC interaction energy.
1208 C
1209         do iint=1,nint_gr(i)
1210           do j=istart(i,iint),iend(i,iint)
1211             ind=ind+1
1212             itypj=itype(j)
1213             if (itypj.eq.ntyp1) cycle
1214 c            dscj_inv=dsc_inv(itypj)
1215             dscj_inv=vbld_inv(j+nres)
1216             sig0ij=sigma(itypi,itypj)
1217             r0ij=r0(itypi,itypj)
1218             chi1=chi(itypi,itypj)
1219             chi2=chi(itypj,itypi)
1220             chi12=chi1*chi2
1221             chip1=chip(itypi)
1222             chip2=chip(itypj)
1223             chip12=chip1*chip2
1224             alf1=alp(itypi)
1225             alf2=alp(itypj)
1226             alf12=0.5D0*(alf1+alf2)
1227             xj=c(1,nres+j)-xi
1228             yj=c(2,nres+j)-yi
1229             zj=c(3,nres+j)-zi
1230             dxj=dc_norm(1,nres+j)
1231             dyj=dc_norm(2,nres+j)
1232             dzj=dc_norm(3,nres+j)
1233             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1234             rij=dsqrt(rrij)
1235
1236             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1237
1238             if (sss.gt.0.0d0) then
1239
1240 C Calculate angle-dependent terms of energy and contributions to their
1241 C derivatives.
1242               call sc_angular
1243               sigsq=1.0D0/sigsq
1244               sig=sig0ij*dsqrt(sigsq)
1245               rij_shift=1.0D0/rij-sig+r0ij
1246 C I hate to put IF's in the loops, but here don't have another choice!!!!
1247               if (rij_shift.le.0.0D0) then
1248                 evdw=1.0D20
1249                 return
1250               endif
1251               sigder=-sig*sigsq
1252 c---------------------------------------------------------------
1253               rij_shift=1.0D0/rij_shift 
1254               fac=rij_shift**expon
1255               e1=fac*fac*aa
1256               e2=fac*bb
1257               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1258               eps2der=evdwij*eps3rt
1259               eps3der=evdwij*eps2rt
1260               fac_augm=rrij**expon
1261               e_augm=augm(itypi,itypj)*fac_augm
1262               evdwij=evdwij*eps2rt*eps3rt
1263               evdw=evdw+(evdwij+e_augm)*sss
1264               if (lprn) then
1265               sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1266               epsi=bb**2/aa
1267               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1268      &          restyp(itypi),i,restyp(itypj),j,
1269      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1270      &          chi1,chi2,chip1,chip2,
1271      &          eps1,eps2rt**2,eps3rt**2,
1272      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1273      &          evdwij+e_augm
1274               endif
1275 C Calculate gradient components.
1276               e1=e1*eps1*eps2rt**2*eps3rt**2
1277               fac=-expon*(e1+evdwij)*rij_shift
1278               sigder=fac*sigder
1279               fac=rij*fac-2*expon*rrij*e_augm
1280 C Calculate the radial part of the gradient
1281               gg(1)=xj*fac
1282               gg(2)=yj*fac
1283               gg(3)=zj*fac
1284 C Calculate angular part of the gradient.
1285               call sc_grad_scale(sss)
1286             endif
1287           enddo      ! j
1288         enddo        ! iint
1289       enddo          ! i
1290       end
1291 C----------------------------------------------------------------------------
1292       subroutine sc_grad_scale(scalfac)
1293       implicit real*8 (a-h,o-z)
1294       include 'DIMENSIONS'
1295       include 'COMMON.CHAIN'
1296       include 'COMMON.DERIV'
1297       include 'COMMON.CALC'
1298       include 'COMMON.IOUNITS'
1299       double precision dcosom1(3),dcosom2(3)
1300       double precision scalfac
1301       eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1302       eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1303       eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1304      &     -2.0D0*alf12*eps3der+sigder*sigsq_om12
1305 c diagnostics only
1306 c      eom1=0.0d0
1307 c      eom2=0.0d0
1308 c      eom12=evdwij*eps1_om12
1309 c end diagnostics
1310 c      write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1311 c     &  " sigder",sigder
1312 c      write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1313 c      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1314       do k=1,3
1315         dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1316         dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1317       enddo
1318       do k=1,3
1319         gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1320       enddo 
1321 c      write (iout,*) "gg",(gg(k),k=1,3)
1322       do k=1,3
1323         gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
1324      &            +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1325      &          +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1326         gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
1327      &            +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1328      &          +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1329 c        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1330 c     &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1331 c        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1332 c     &            +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1333       enddo
1334
1335 C Calculate the components of the gradient in DC and X
1336 C
1337       do l=1,3
1338         gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
1339         gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
1340       enddo
1341       return
1342       end
1343 C--------------------------------------------------------------------------
1344       subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1345 C
1346 C This subroutine calculates the average interaction energy and its gradient
1347 C in the virtual-bond vectors between non-adjacent peptide groups, based on 
1348 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715. 
1349 C The potential depends both on the distance of peptide-group centers and on 
1350 C the orientation of the CA-CA virtual bonds.
1351
1352       implicit real*8 (a-h,o-z)
1353 #ifdef MPI
1354       include 'mpif.h'
1355 #endif
1356       include 'DIMENSIONS'
1357       include 'COMMON.CONTROL'
1358       include 'COMMON.SETUP'
1359       include 'COMMON.IOUNITS'
1360       include 'COMMON.GEO'
1361       include 'COMMON.VAR'
1362       include 'COMMON.LOCAL'
1363       include 'COMMON.CHAIN'
1364       include 'COMMON.DERIV'
1365       include 'COMMON.INTERACT'
1366       include 'COMMON.CONTACTS'
1367       include 'COMMON.TORSION'
1368       include 'COMMON.VECTORS'
1369       include 'COMMON.FFIELD'
1370       include 'COMMON.TIME1'
1371       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1372      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1373       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1374      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1375       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1376      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1377      &    num_conti,j1,j2
1378 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1379 #ifdef MOMENT
1380       double precision scal_el /1.0d0/
1381 #else
1382       double precision scal_el /0.5d0/
1383 #endif
1384 C 12/13/98 
1385 C 13-go grudnia roku pamietnego... 
1386       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1387      &                   0.0d0,1.0d0,0.0d0,
1388      &                   0.0d0,0.0d0,1.0d0/
1389 cd      write(iout,*) 'In EELEC'
1390 cd      do i=1,nloctyp
1391 cd        write(iout,*) 'Type',i
1392 cd        write(iout,*) 'B1',B1(:,i)
1393 cd        write(iout,*) 'B2',B2(:,i)
1394 cd        write(iout,*) 'CC',CC(:,:,i)
1395 cd        write(iout,*) 'DD',DD(:,:,i)
1396 cd        write(iout,*) 'EE',EE(:,:,i)
1397 cd      enddo
1398 cd      call check_vecgrad
1399 cd      stop
1400 C      print *,"WCHODZE3"
1401       if (icheckgrad.eq.1) then
1402         do i=1,nres-1
1403           fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1404           do k=1,3
1405             dc_norm(k,i)=dc(k,i)*fac
1406           enddo
1407 c          write (iout,*) 'i',i,' fac',fac
1408         enddo
1409       endif
1410       if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 
1411      &    .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or. 
1412      &    wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1413 c        call vec_and_deriv
1414 #ifdef TIMING
1415         time01=MPI_Wtime()
1416 #endif
1417         call set_matrices
1418 #ifdef TIMING
1419         time_mat=time_mat+MPI_Wtime()-time01
1420 #endif
1421       endif
1422 cd      do i=1,nres-1
1423 cd        write (iout,*) 'i=',i
1424 cd        do k=1,3
1425 cd        write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1426 cd        enddo
1427 cd        do k=1,3
1428 cd          write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)') 
1429 cd     &     uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1430 cd        enddo
1431 cd      enddo
1432       t_eelecij=0.0d0
1433       ees=0.0D0
1434       evdw1=0.0D0
1435       eel_loc=0.0d0 
1436       eello_turn3=0.0d0
1437       eello_turn4=0.0d0
1438       ind=0
1439       do i=1,nres
1440         num_cont_hb(i)=0
1441       enddo
1442 cd      print '(a)','Enter EELEC'
1443 cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1444       do i=1,nres
1445         gel_loc_loc(i)=0.0d0
1446         gcorr_loc(i)=0.0d0
1447       enddo
1448 c
1449 c
1450 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1451 C
1452 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1453 C
1454       do i=iturn3_start,iturn3_end
1455         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1456      &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1457      &  .or. itype(i-1).eq.ntyp1
1458      &  .or. itype(i+4).eq.ntyp1
1459      &   ) cycle
1460         dxi=dc(1,i)
1461         dyi=dc(2,i)
1462         dzi=dc(3,i)
1463         dx_normi=dc_norm(1,i)
1464         dy_normi=dc_norm(2,i)
1465         dz_normi=dc_norm(3,i)
1466         xmedi=c(1,i)+0.5d0*dxi
1467         ymedi=c(2,i)+0.5d0*dyi
1468         zmedi=c(3,i)+0.5d0*dzi
1469           xmedi=mod(xmedi,boxxsize)
1470           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1471           ymedi=mod(ymedi,boxysize)
1472           if (ymedi.lt.0) ymedi=ymedi+boxysize
1473           zmedi=mod(zmedi,boxzsize)
1474           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1475         num_conti=0
1476         call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1477         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1478         num_cont_hb(i)=num_conti
1479       enddo
1480       do i=iturn4_start,iturn4_end
1481         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1482      &    .or. itype(i+3).eq.ntyp1
1483      &    .or. itype(i+4).eq.ntyp1
1484      &    .or. itype(i+5).eq.ntyp1
1485      &    .or. itype(i-1).eq.ntyp1
1486      &    ) cycle
1487         dxi=dc(1,i)
1488         dyi=dc(2,i)
1489         dzi=dc(3,i)
1490         dx_normi=dc_norm(1,i)
1491         dy_normi=dc_norm(2,i)
1492         dz_normi=dc_norm(3,i)
1493         xmedi=c(1,i)+0.5d0*dxi
1494         ymedi=c(2,i)+0.5d0*dyi
1495         zmedi=c(3,i)+0.5d0*dzi
1496           xmedi=mod(xmedi,boxxsize)
1497           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1498           ymedi=mod(ymedi,boxysize)
1499           if (ymedi.lt.0) ymedi=ymedi+boxysize
1500           zmedi=mod(zmedi,boxzsize)
1501           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1502         num_conti=num_cont_hb(i)
1503         call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1504         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) 
1505      &    call eturn4(i,eello_turn4)
1506         num_cont_hb(i)=num_conti
1507       enddo   ! i
1508 c
1509 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1510 c
1511       do i=iatel_s,iatel_e
1512         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1513      &  .or. itype(i+2).eq.ntyp1
1514      &  .or. itype(i-1).eq.ntyp1
1515      &) cycle
1516         dxi=dc(1,i)
1517         dyi=dc(2,i)
1518         dzi=dc(3,i)
1519         dx_normi=dc_norm(1,i)
1520         dy_normi=dc_norm(2,i)
1521         dz_normi=dc_norm(3,i)
1522         xmedi=c(1,i)+0.5d0*dxi
1523         ymedi=c(2,i)+0.5d0*dyi
1524         zmedi=c(3,i)+0.5d0*dzi
1525           xmedi=mod(xmedi,boxxsize)
1526           if (xmedi.lt.0) xmedi=xmedi+boxxsize
1527           ymedi=mod(ymedi,boxysize)
1528           if (ymedi.lt.0) ymedi=ymedi+boxysize
1529           zmedi=mod(zmedi,boxzsize)
1530           if (zmedi.lt.0) zmedi=zmedi+boxzsize
1531 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1532         num_conti=num_cont_hb(i)
1533         do j=ielstart(i),ielend(i)
1534           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1535      & .or.itype(j+2).eq.ntyp1
1536      & .or.itype(j-1).eq.ntyp1
1537      &) cycle
1538           call eelecij_scale(i,j,ees,evdw1,eel_loc)
1539         enddo ! j
1540         num_cont_hb(i)=num_conti
1541       enddo   ! i
1542 c      write (iout,*) "Number of loop steps in EELEC:",ind
1543 cd      do i=1,nres
1544 cd        write (iout,'(i3,3f10.5,5x,3f10.5)') 
1545 cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1546 cd      enddo
1547 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1548 ccc      eel_loc=eel_loc+eello_turn3
1549 cd      print *,"Processor",fg_rank," t_eelecij",t_eelecij
1550       return
1551       end
1552 C-------------------------------------------------------------------------------
1553       subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1554       implicit real*8 (a-h,o-z)
1555       include 'DIMENSIONS'
1556 #ifdef MPI
1557       include "mpif.h"
1558 #endif
1559       include 'COMMON.CONTROL'
1560       include 'COMMON.IOUNITS'
1561       include 'COMMON.GEO'
1562       include 'COMMON.VAR'
1563       include 'COMMON.LOCAL'
1564       include 'COMMON.CHAIN'
1565       include 'COMMON.DERIV'
1566       include 'COMMON.INTERACT'
1567       include 'COMMON.CONTACTS'
1568       include 'COMMON.TORSION'
1569       include 'COMMON.VECTORS'
1570       include 'COMMON.FFIELD'
1571       include 'COMMON.TIME1'
1572       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1573      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1574       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1575      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1576       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1577      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1578      &    num_conti,j1,j2
1579 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1580 #ifdef MOMENT
1581       double precision scal_el /1.0d0/
1582 #else
1583       double precision scal_el /0.5d0/
1584 #endif
1585 C 12/13/98 
1586 C 13-go grudnia roku pamietnego... 
1587       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1588      &                   0.0d0,1.0d0,0.0d0,
1589      &                   0.0d0,0.0d0,1.0d0/
1590 c          time00=MPI_Wtime()
1591 cd      write (iout,*) "eelecij",i,j
1592 C      print *,"WCHODZE2"
1593           ind=ind+1
1594           iteli=itel(i)
1595           itelj=itel(j)
1596           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1597           aaa=app(iteli,itelj)
1598           bbb=bpp(iteli,itelj)
1599           ael6i=ael6(iteli,itelj)
1600           ael3i=ael3(iteli,itelj) 
1601           dxj=dc(1,j)
1602           dyj=dc(2,j)
1603           dzj=dc(3,j)
1604           dx_normj=dc_norm(1,j)
1605           dy_normj=dc_norm(2,j)
1606           dz_normj=dc_norm(3,j)
1607           xj=c(1,j)+0.5D0*dxj
1608           yj=c(2,j)+0.5D0*dyj
1609           zj=c(3,j)+0.5D0*dzj
1610          xj=mod(xj,boxxsize)
1611           if (xj.lt.0) xj=xj+boxxsize
1612           yj=mod(yj,boxysize)
1613           if (yj.lt.0) yj=yj+boxysize
1614           zj=mod(zj,boxzsize)
1615           if (zj.lt.0) zj=zj+boxzsize
1616       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1617       xj_safe=xj
1618       yj_safe=yj
1619       zj_safe=zj
1620       isubchap=0
1621       do xshift=-1,1
1622       do yshift=-1,1
1623       do zshift=-1,1
1624           xj=xj_safe+xshift*boxxsize
1625           yj=yj_safe+yshift*boxysize
1626           zj=zj_safe+zshift*boxzsize
1627           dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1628           if(dist_temp.lt.dist_init) then
1629             dist_init=dist_temp
1630             xj_temp=xj
1631             yj_temp=yj
1632             zj_temp=zj
1633             isubchap=1
1634           endif
1635        enddo
1636        enddo
1637        enddo
1638        if (isubchap.eq.1) then
1639           xj=xj_temp-xmedi
1640           yj=yj_temp-ymedi
1641           zj=zj_temp-zmedi
1642        else
1643           xj=xj_safe-xmedi
1644           yj=yj_safe-ymedi
1645           zj=zj_safe-zmedi
1646        endif
1647
1648           rij=xj*xj+yj*yj+zj*zj
1649           rrmij=1.0D0/rij
1650           rij=dsqrt(rij)
1651           rmij=1.0D0/rij
1652 c For extracting the short-range part of Evdwpp
1653           sss=sscale(rij/rpp(iteli,itelj))
1654           sssgrad=sscagrad(rij/rpp(iteli,itelj))
1655           r3ij=rrmij*rmij
1656           r6ij=r3ij*r3ij  
1657           cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1658           cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1659           cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1660           fac=cosa-3.0D0*cosb*cosg
1661           ev1=aaa*r6ij*r6ij
1662 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1663           if (j.eq.i+2) ev1=scal_el*ev1
1664           ev2=bbb*r6ij
1665           fac3=ael6i*r6ij
1666           fac4=ael3i*r3ij
1667           evdwij=ev1+ev2
1668           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1669           el2=fac4*fac       
1670           eesij=el1+el2
1671 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1672           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1673           ees=ees+eesij
1674           evdw1=evdw1+evdwij*(1.0d0-sss)
1675 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1676 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1677 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
1678 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
1679
1680           if (energy_dec) then 
1681               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1682               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1683           endif
1684
1685 C
1686 C Calculate contributions to the Cartesian gradient.
1687 C
1688 #ifdef SPLITELE
1689           facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1690           facel=-3*rrmij*(el1+eesij)
1691           fac1=fac
1692           erij(1)=xj*rmij
1693           erij(2)=yj*rmij
1694           erij(3)=zj*rmij
1695 *
1696 * Radial derivatives. First process both termini of the fragment (i,j)
1697 *
1698           ggg(1)=facel*xj
1699           ggg(2)=facel*yj
1700           ggg(3)=facel*zj
1701 c          do k=1,3
1702 c            ghalf=0.5D0*ggg(k)
1703 c            gelc(k,i)=gelc(k,i)+ghalf
1704 c            gelc(k,j)=gelc(k,j)+ghalf
1705 c          enddo
1706 c 9/28/08 AL Gradient compotents will be summed only at the end
1707           do k=1,3
1708             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1709             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1710           enddo
1711 *
1712 * Loop over residues i+1 thru j-1.
1713 *
1714 cgrad          do k=i+1,j-1
1715 cgrad            do l=1,3
1716 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1717 cgrad            enddo
1718 cgrad          enddo
1719           ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1720           ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1721           ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1722 c          do k=1,3
1723 c            ghalf=0.5D0*ggg(k)
1724 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1725 c            gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1726 c          enddo
1727 c 9/28/08 AL Gradient compotents will be summed only at the end
1728           do k=1,3
1729             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1730             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1731           enddo
1732 *
1733 * Loop over residues i+1 thru j-1.
1734 *
1735 cgrad          do k=i+1,j-1
1736 cgrad            do l=1,3
1737 cgrad              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1738 cgrad            enddo
1739 cgrad          enddo
1740 #else
1741           facvdw=ev1+evdwij*(1.0d0-sss) 
1742           facel=el1+eesij  
1743           fac1=fac
1744           fac=-3*rrmij*(facvdw+facvdw+facel)
1745           erij(1)=xj*rmij
1746           erij(2)=yj*rmij
1747           erij(3)=zj*rmij
1748 *
1749 * Radial derivatives. First process both termini of the fragment (i,j)
1750
1751           ggg(1)=fac*xj
1752           ggg(2)=fac*yj
1753           ggg(3)=fac*zj
1754 c          do k=1,3
1755 c            ghalf=0.5D0*ggg(k)
1756 c            gelc(k,i)=gelc(k,i)+ghalf
1757 c            gelc(k,j)=gelc(k,j)+ghalf
1758 c          enddo
1759 c 9/28/08 AL Gradient compotents will be summed only at the end
1760           do k=1,3
1761             gelc_long(k,j)=gelc(k,j)+ggg(k)
1762             gelc_long(k,i)=gelc(k,i)-ggg(k)
1763           enddo
1764 *
1765 * Loop over residues i+1 thru j-1.
1766 *
1767 cgrad          do k=i+1,j-1
1768 cgrad            do l=1,3
1769 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1770 cgrad            enddo
1771 cgrad          enddo
1772 c 9/28/08 AL Gradient compotents will be summed only at the end
1773 C          ggg(1)=facvdw*xj
1774 C          ggg(2)=facvdw*yj
1775 C          ggg(3)=facvdw*zj
1776           ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1777           ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1778           ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1779           do k=1,3
1780             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1781             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1782           enddo
1783 #endif
1784 *
1785 * Angular part
1786 *          
1787           ecosa=2.0D0*fac3*fac1+fac4
1788           fac4=-3.0D0*fac4
1789           fac3=-6.0D0*fac3
1790           ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1791           ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1792           do k=1,3
1793             dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1794             dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1795           enddo
1796 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1797 cd   &          (dcosg(k),k=1,3)
1798           do k=1,3
1799             ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k) 
1800           enddo
1801 c          do k=1,3
1802 c            ghalf=0.5D0*ggg(k)
1803 c            gelc(k,i)=gelc(k,i)+ghalf
1804 c     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1805 c     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1806 c            gelc(k,j)=gelc(k,j)+ghalf
1807 c     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1808 c     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1809 c          enddo
1810 cgrad          do k=i+1,j-1
1811 cgrad            do l=1,3
1812 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
1813 cgrad            enddo
1814 cgrad          enddo
1815           do k=1,3
1816             gelc(k,i)=gelc(k,i)
1817      &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1818      &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1819             gelc(k,j)=gelc(k,j)
1820      &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1821      &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1822             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1823             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1824           enddo
1825           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1826      &        .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 
1827      &        .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1828 C
1829 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction 
1830 C   energy of a peptide unit is assumed in the form of a second-order 
1831 C   Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1832 C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1833 C   are computed for EVERY pair of non-contiguous peptide groups.
1834 C
1835           if (j.lt.nres-1) then
1836             j1=j+1
1837             j2=j-1
1838           else
1839             j1=j-1
1840             j2=j-2
1841           endif
1842           kkk=0
1843           do k=1,2
1844             do l=1,2
1845               kkk=kkk+1
1846               muij(kkk)=mu(k,i)*mu(l,j)
1847             enddo
1848           enddo  
1849 cd         write (iout,*) 'EELEC: i',i,' j',j
1850 cd          write (iout,*) 'j',j,' j1',j1,' j2',j2
1851 cd          write(iout,*) 'muij',muij
1852           ury=scalar(uy(1,i),erij)
1853           urz=scalar(uz(1,i),erij)
1854           vry=scalar(uy(1,j),erij)
1855           vrz=scalar(uz(1,j),erij)
1856           a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1857           a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1858           a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1859           a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1860           fac=dsqrt(-ael6i)*r3ij
1861           a22=a22*fac
1862           a23=a23*fac
1863           a32=a32*fac
1864           a33=a33*fac
1865 cd          write (iout,'(4i5,4f10.5)')
1866 cd     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1867 cd          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1868 cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1869 cd     &      uy(:,j),uz(:,j)
1870 cd          write (iout,'(4f10.5)') 
1871 cd     &      scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1872 cd     &      scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1873 cd          write (iout,'(4f10.5)') ury,urz,vry,vrz
1874 cd           write (iout,'(9f10.5/)') 
1875 cd     &      fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1876 C Derivatives of the elements of A in virtual-bond vectors
1877           call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1878           do k=1,3
1879             uryg(k,1)=scalar(erder(1,k),uy(1,i))
1880             uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1881             uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1882             urzg(k,1)=scalar(erder(1,k),uz(1,i))
1883             urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1884             urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1885             vryg(k,1)=scalar(erder(1,k),uy(1,j))
1886             vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1887             vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1888             vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1889             vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1890             vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1891           enddo
1892 C Compute radial contributions to the gradient
1893           facr=-3.0d0*rrmij
1894           a22der=a22*facr
1895           a23der=a23*facr
1896           a32der=a32*facr
1897           a33der=a33*facr
1898           agg(1,1)=a22der*xj
1899           agg(2,1)=a22der*yj
1900           agg(3,1)=a22der*zj
1901           agg(1,2)=a23der*xj
1902           agg(2,2)=a23der*yj
1903           agg(3,2)=a23der*zj
1904           agg(1,3)=a32der*xj
1905           agg(2,3)=a32der*yj
1906           agg(3,3)=a32der*zj
1907           agg(1,4)=a33der*xj
1908           agg(2,4)=a33der*yj
1909           agg(3,4)=a33der*zj
1910 C Add the contributions coming from er
1911           fac3=-3.0d0*fac
1912           do k=1,3
1913             agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1914             agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1915             agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1916             agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1917           enddo
1918           do k=1,3
1919 C Derivatives in DC(i) 
1920 cgrad            ghalf1=0.5d0*agg(k,1)
1921 cgrad            ghalf2=0.5d0*agg(k,2)
1922 cgrad            ghalf3=0.5d0*agg(k,3)
1923 cgrad            ghalf4=0.5d0*agg(k,4)
1924             aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1925      &      -3.0d0*uryg(k,2)*vry)!+ghalf1
1926             aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1927      &      -3.0d0*uryg(k,2)*vrz)!+ghalf2
1928             aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1929      &      -3.0d0*urzg(k,2)*vry)!+ghalf3
1930             aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1931      &      -3.0d0*urzg(k,2)*vrz)!+ghalf4
1932 C Derivatives in DC(i+1)
1933             aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
1934      &      -3.0d0*uryg(k,3)*vry)!+agg(k,1)
1935             aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
1936      &      -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
1937             aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
1938      &      -3.0d0*urzg(k,3)*vry)!+agg(k,3)
1939             aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
1940      &      -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
1941 C Derivatives in DC(j)
1942             aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
1943      &      -3.0d0*vryg(k,2)*ury)!+ghalf1
1944             aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
1945      &      -3.0d0*vrzg(k,2)*ury)!+ghalf2
1946             aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
1947      &      -3.0d0*vryg(k,2)*urz)!+ghalf3
1948             aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i)) 
1949      &      -3.0d0*vrzg(k,2)*urz)!+ghalf4
1950 C Derivatives in DC(j+1) or DC(nres-1)
1951             aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
1952      &      -3.0d0*vryg(k,3)*ury)
1953             aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
1954      &      -3.0d0*vrzg(k,3)*ury)
1955             aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
1956      &      -3.0d0*vryg(k,3)*urz)
1957             aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i)) 
1958      &      -3.0d0*vrzg(k,3)*urz)
1959 cgrad            if (j.eq.nres-1 .and. i.lt.j-2) then
1960 cgrad              do l=1,4
1961 cgrad                aggj1(k,l)=aggj1(k,l)+agg(k,l)
1962 cgrad              enddo
1963 cgrad            endif
1964           enddo
1965           acipa(1,1)=a22
1966           acipa(1,2)=a23
1967           acipa(2,1)=a32
1968           acipa(2,2)=a33
1969           a22=-a22
1970           a23=-a23
1971           do l=1,2
1972             do k=1,3
1973               agg(k,l)=-agg(k,l)
1974               aggi(k,l)=-aggi(k,l)
1975               aggi1(k,l)=-aggi1(k,l)
1976               aggj(k,l)=-aggj(k,l)
1977               aggj1(k,l)=-aggj1(k,l)
1978             enddo
1979           enddo
1980           if (j.lt.nres-1) then
1981             a22=-a22
1982             a32=-a32
1983             do l=1,3,2
1984               do k=1,3
1985                 agg(k,l)=-agg(k,l)
1986                 aggi(k,l)=-aggi(k,l)
1987                 aggi1(k,l)=-aggi1(k,l)
1988                 aggj(k,l)=-aggj(k,l)
1989                 aggj1(k,l)=-aggj1(k,l)
1990               enddo
1991             enddo
1992           else
1993             a22=-a22
1994             a23=-a23
1995             a32=-a32
1996             a33=-a33
1997             do l=1,4
1998               do k=1,3
1999                 agg(k,l)=-agg(k,l)
2000                 aggi(k,l)=-aggi(k,l)
2001                 aggi1(k,l)=-aggi1(k,l)
2002                 aggj(k,l)=-aggj(k,l)
2003                 aggj1(k,l)=-aggj1(k,l)
2004               enddo
2005             enddo 
2006           endif    
2007           ENDIF ! WCORR
2008           IF (wel_loc.gt.0.0d0) THEN
2009 C Contribution to the local-electrostatic energy coming from the i-j pair
2010           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2011      &     +a33*muij(4)
2012 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2013
2014           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2015      &            'eelloc',i,j,eel_loc_ij
2016
2017           eel_loc=eel_loc+eel_loc_ij
2018 C Partial derivatives in virtual-bond dihedral angles gamma
2019           if (i.gt.1)
2020      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+ 
2021      &            a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2022      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
2023           gel_loc_loc(j-1)=gel_loc_loc(j-1)+ 
2024      &            a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2025      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
2026 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2027           do l=1,3
2028             ggg(l)=agg(l,1)*muij(1)+
2029      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
2030             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2031             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2032 cgrad            ghalf=0.5d0*ggg(l)
2033 cgrad            gel_loc(l,i)=gel_loc(l,i)+ghalf
2034 cgrad            gel_loc(l,j)=gel_loc(l,j)+ghalf
2035           enddo
2036 cgrad          do k=i+1,j2
2037 cgrad            do l=1,3
2038 cgrad              gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2039 cgrad            enddo
2040 cgrad          enddo
2041 C Remaining derivatives of eello
2042           do l=1,3
2043             gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
2044      &          aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
2045             gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
2046      &          aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
2047             gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
2048      &          aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
2049             gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
2050      &          aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
2051           enddo
2052           ENDIF
2053 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2054 c          if (j.gt.i+1 .and. num_conti.le.maxconts) then
2055           if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2056      &       .and. num_conti.le.maxconts) then
2057 c            write (iout,*) i,j," entered corr"
2058 C
2059 C Calculate the contact function. The ith column of the array JCONT will 
2060 C contain the numbers of atoms that make contacts with the atom I (of numbers
2061 C greater than I). The arrays FACONT and GACONT will contain the values of
2062 C the contact function and its derivative.
2063 c           r0ij=1.02D0*rpp(iteli,itelj)
2064 c           r0ij=1.11D0*rpp(iteli,itelj)
2065             r0ij=2.20D0*rpp(iteli,itelj)
2066 c           r0ij=1.55D0*rpp(iteli,itelj)
2067             call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2068             if (fcont.gt.0.0D0) then
2069               num_conti=num_conti+1
2070               if (num_conti.gt.maxconts) then
2071                 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2072      &                         ' will skip next contacts for this conf.'
2073               else
2074                 jcont_hb(num_conti,i)=j
2075 cd                write (iout,*) "i",i," j",j," num_conti",num_conti,
2076 cd     &           " jcont_hb",jcont_hb(num_conti,i)
2077                 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. 
2078      &          wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2079 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2080 C  terms.
2081                 d_cont(num_conti,i)=rij
2082 cd                write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2083 C     --- Electrostatic-interaction matrix --- 
2084                 a_chuj(1,1,num_conti,i)=a22
2085                 a_chuj(1,2,num_conti,i)=a23
2086                 a_chuj(2,1,num_conti,i)=a32
2087                 a_chuj(2,2,num_conti,i)=a33
2088 C     --- Gradient of rij
2089                 do kkk=1,3
2090                   grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2091                 enddo
2092                 kkll=0
2093                 do k=1,2
2094                   do l=1,2
2095                     kkll=kkll+1
2096                     do m=1,3
2097                       a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2098                       a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2099                       a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2100                       a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2101                       a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2102                     enddo
2103                   enddo
2104                 enddo
2105                 ENDIF
2106                 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2107 C Calculate contact energies
2108                 cosa4=4.0D0*cosa
2109                 wij=cosa-3.0D0*cosb*cosg
2110                 cosbg1=cosb+cosg
2111                 cosbg2=cosb-cosg
2112 c               fac3=dsqrt(-ael6i)/r0ij**3     
2113                 fac3=dsqrt(-ael6i)*r3ij
2114 c                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2115                 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2116                 if (ees0tmp.gt.0) then
2117                   ees0pij=dsqrt(ees0tmp)
2118                 else
2119                   ees0pij=0
2120                 endif
2121 c                ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2122                 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2123                 if (ees0tmp.gt.0) then
2124                   ees0mij=dsqrt(ees0tmp)
2125                 else
2126                   ees0mij=0
2127                 endif
2128 c               ees0mij=0.0D0
2129                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2130                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2131 C Diagnostics. Comment out or remove after debugging!
2132 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2133 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2134 c               ees0m(num_conti,i)=0.0D0
2135 C End diagnostics.
2136 c               write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2137 c    & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2138 C Angular derivatives of the contact function
2139                 ees0pij1=fac3/ees0pij 
2140                 ees0mij1=fac3/ees0mij
2141                 fac3p=-3.0D0*fac3*rrmij
2142                 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2143                 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2144 c               ees0mij1=0.0D0
2145                 ecosa1=       ees0pij1*( 1.0D0+0.5D0*wij)
2146                 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2147                 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2148                 ecosa2=       ees0mij1*(-1.0D0+0.5D0*wij)
2149                 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2) 
2150                 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2151                 ecosap=ecosa1+ecosa2
2152                 ecosbp=ecosb1+ecosb2
2153                 ecosgp=ecosg1+ecosg2
2154                 ecosam=ecosa1-ecosa2
2155                 ecosbm=ecosb1-ecosb2
2156                 ecosgm=ecosg1-ecosg2
2157 C Diagnostics
2158 c               ecosap=ecosa1
2159 c               ecosbp=ecosb1
2160 c               ecosgp=ecosg1
2161 c               ecosam=0.0D0
2162 c               ecosbm=0.0D0
2163 c               ecosgm=0.0D0
2164 C End diagnostics
2165                 facont_hb(num_conti,i)=fcont
2166                 fprimcont=fprimcont/rij
2167 cd              facont_hb(num_conti,i)=1.0D0
2168 C Following line is for diagnostics.
2169 cd              fprimcont=0.0D0
2170                 do k=1,3
2171                   dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2172                   dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2173                 enddo
2174                 do k=1,3
2175                   gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2176                   gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2177                 enddo
2178                 gggp(1)=gggp(1)+ees0pijp*xj
2179                 gggp(2)=gggp(2)+ees0pijp*yj
2180                 gggp(3)=gggp(3)+ees0pijp*zj
2181                 gggm(1)=gggm(1)+ees0mijp*xj
2182                 gggm(2)=gggm(2)+ees0mijp*yj
2183                 gggm(3)=gggm(3)+ees0mijp*zj
2184 C Derivatives due to the contact function
2185                 gacont_hbr(1,num_conti,i)=fprimcont*xj
2186                 gacont_hbr(2,num_conti,i)=fprimcont*yj
2187                 gacont_hbr(3,num_conti,i)=fprimcont*zj
2188                 do k=1,3
2189 c
2190 c 10/24/08 cgrad and ! comments indicate the parts of the code removed 
2191 c          following the change of gradient-summation algorithm.
2192 c
2193 cgrad                  ghalfp=0.5D0*gggp(k)
2194 cgrad                  ghalfm=0.5D0*gggm(k)
2195                   gacontp_hb1(k,num_conti,i)=!ghalfp
2196      &              +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2197      &              + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2198                   gacontp_hb2(k,num_conti,i)=!ghalfp
2199      &              +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2200      &              + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2201                   gacontp_hb3(k,num_conti,i)=gggp(k)
2202                   gacontm_hb1(k,num_conti,i)=!ghalfm
2203      &              +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2204      &              + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2205                   gacontm_hb2(k,num_conti,i)=!ghalfm
2206      &              +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2207      &              + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2208                   gacontm_hb3(k,num_conti,i)=gggm(k)
2209                 enddo
2210               ENDIF ! wcorr
2211               endif  ! num_conti.le.maxconts
2212             endif  ! fcont.gt.0
2213           endif    ! j.gt.i+1
2214           if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2215             do k=1,4
2216               do l=1,3
2217                 ghalf=0.5d0*agg(l,k)
2218                 aggi(l,k)=aggi(l,k)+ghalf
2219                 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2220                 aggj(l,k)=aggj(l,k)+ghalf
2221               enddo
2222             enddo
2223             if (j.eq.nres-1 .and. i.lt.j-2) then
2224               do k=1,4
2225                 do l=1,3
2226                   aggj1(l,k)=aggj1(l,k)+agg(l,k)
2227                 enddo
2228               enddo
2229             endif
2230           endif
2231 c          t_eelecij=t_eelecij+MPI_Wtime()-time00
2232       return
2233       end
2234 C-----------------------------------------------------------------------
2235       subroutine evdwpp_short(evdw1)
2236 C
2237 C Compute Evdwpp
2238
2239       implicit real*8 (a-h,o-z)
2240       include 'DIMENSIONS'
2241       include 'COMMON.CONTROL'
2242       include 'COMMON.IOUNITS'
2243       include 'COMMON.GEO'
2244       include 'COMMON.VAR'
2245       include 'COMMON.LOCAL'
2246       include 'COMMON.CHAIN'
2247       include 'COMMON.DERIV'
2248       include 'COMMON.INTERACT'
2249       include 'COMMON.CONTACTS'
2250       include 'COMMON.TORSION'
2251       include 'COMMON.VECTORS'
2252       include 'COMMON.FFIELD'
2253       dimension ggg(3)
2254 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2255 #ifdef MOMENT
2256       double precision scal_el /1.0d0/
2257 #else
2258       double precision scal_el /0.5d0/
2259 #endif
2260       evdw1=0.0D0
2261 C      print *,"WCHODZE"
2262 c      write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2263 c     & " iatel_e_vdw",iatel_e_vdw
2264       call flush(iout)
2265       do i=iatel_s_vdw,iatel_e_vdw
2266         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2267         dxi=dc(1,i)
2268         dyi=dc(2,i)
2269         dzi=dc(3,i)
2270         dx_normi=dc_norm(1,i)
2271         dy_normi=dc_norm(2,i)
2272         dz_normi=dc_norm(3,i)
2273         xmedi=c(1,i)+0.5d0*dxi
2274         ymedi=c(2,i)+0.5d0*dyi
2275         zmedi=c(3,i)+0.5d0*dzi
2276           xmedi=mod(xmedi,boxxsize)
2277           if (xmedi.lt.0) xmedi=xmedi+boxxsize
2278           ymedi=mod(ymedi,boxysize)
2279           if (ymedi.lt.0) ymedi=ymedi+boxysize
2280           zmedi=mod(zmedi,boxzsize)
2281           if (zmedi.lt.0) zmedi=zmedi+boxzsize
2282         num_conti=0
2283 c        write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2284 c     &   ' ielend',ielend_vdw(i)
2285         call flush(iout)
2286         do j=ielstart_vdw(i),ielend_vdw(i)
2287           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2288           ind=ind+1
2289           iteli=itel(i)
2290           itelj=itel(j)
2291           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2292           aaa=app(iteli,itelj)
2293           bbb=bpp(iteli,itelj)
2294           dxj=dc(1,j)
2295           dyj=dc(2,j)
2296           dzj=dc(3,j)
2297           dx_normj=dc_norm(1,j)
2298           dy_normj=dc_norm(2,j)
2299           dz_normj=dc_norm(3,j)
2300           xj=c(1,j)+0.5D0*dxj
2301           yj=c(2,j)+0.5D0*dyj
2302           zj=c(3,j)+0.5D0*dzj
2303           xj=mod(xj,boxxsize)
2304           if (xj.lt.0) xj=xj+boxxsize
2305           yj=mod(yj,boxysize)
2306           if (yj.lt.0) yj=yj+boxysize
2307           zj=mod(zj,boxzsize)
2308           if (zj.lt.0) zj=zj+boxzsize
2309       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2310       xj_safe=xj
2311       yj_safe=yj
2312       zj_safe=zj
2313       isubchap=0
2314       do xshift=-1,1
2315       do yshift=-1,1
2316       do zshift=-1,1
2317           xj=xj_safe+xshift*boxxsize
2318           yj=yj_safe+yshift*boxysize
2319           zj=zj_safe+zshift*boxzsize
2320           dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2321           if(dist_temp.lt.dist_init) then
2322             dist_init=dist_temp
2323             xj_temp=xj
2324             yj_temp=yj
2325             zj_temp=zj
2326             isubchap=1
2327           endif
2328        enddo
2329        enddo
2330        enddo
2331        if (isubchap.eq.1) then
2332           xj=xj_temp-xmedi
2333           yj=yj_temp-ymedi
2334           zj=zj_temp-zmedi
2335        else
2336           xj=xj_safe-xmedi
2337           yj=yj_safe-ymedi
2338           zj=zj_safe-zmedi
2339        endif
2340           rij=xj*xj+yj*yj+zj*zj
2341           rrmij=1.0D0/rij
2342           rij=dsqrt(rij)
2343           sss=sscale(rij/rpp(iteli,itelj))
2344             sssgrad=sscagrad(rij/rpp(iteli,itelj))
2345           if (sss.gt.0.0d0) then
2346             rmij=1.0D0/rij
2347             r3ij=rrmij*rmij
2348             r6ij=r3ij*r3ij  
2349             ev1=aaa*r6ij*r6ij
2350 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2351             if (j.eq.i+2) ev1=scal_el*ev1
2352             ev2=bbb*r6ij
2353             evdwij=ev1+ev2
2354             if (energy_dec) then 
2355               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2356             endif
2357             evdw1=evdw1+evdwij*sss
2358 C
2359 C Calculate contributions to the Cartesian gradient.
2360 C
2361             facvdw=-6*rrmij*(ev1+evdwij)*sss
2362           ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2363           ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2364           ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2365 C            ggg(1)=facvdw*xj
2366 C            ggg(2)=facvdw*yj
2367 C            ggg(3)=facvdw*zj
2368             do k=1,3
2369               gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2370               gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2371             enddo
2372           endif
2373         enddo ! j
2374       enddo   ! i
2375       return
2376       end
2377 C-----------------------------------------------------------------------------
2378       subroutine escp_long(evdw2,evdw2_14)
2379 C
2380 C This subroutine calculates the excluded-volume interaction energy between
2381 C peptide-group centers and side chains and its gradient in virtual-bond and
2382 C side-chain vectors.
2383 C
2384       implicit real*8 (a-h,o-z)
2385       include 'DIMENSIONS'
2386       include 'COMMON.GEO'
2387       include 'COMMON.VAR'
2388       include 'COMMON.LOCAL'
2389       include 'COMMON.CHAIN'
2390       include 'COMMON.DERIV'
2391       include 'COMMON.INTERACT'
2392       include 'COMMON.FFIELD'
2393       include 'COMMON.IOUNITS'
2394       include 'COMMON.CONTROL'
2395       dimension ggg(3)
2396       evdw2=0.0D0
2397       evdw2_14=0.0d0
2398 CD        print '(a)','Enter ESCP KURWA'
2399 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2400       do i=iatscp_s,iatscp_e
2401         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2402         iteli=itel(i)
2403         xi=0.5D0*(c(1,i)+c(1,i+1))
2404         yi=0.5D0*(c(2,i)+c(2,i+1))
2405         zi=0.5D0*(c(3,i)+c(3,i+1))
2406          xi=mod(xi,boxxsize)
2407           if (xi.lt.0) xi=xi+boxxsize
2408           yi=mod(yi,boxysize)
2409           if (yi.lt.0) yi=yi+boxysize
2410           zi=mod(zi,boxzsize)
2411           if (zi.lt.0) zi=zi+boxzsize
2412         do iint=1,nscp_gr(i)
2413
2414         do j=iscpstart(i,iint),iscpend(i,iint)
2415           itypj=itype(j)
2416           if (itypj.eq.ntyp1) cycle
2417 C Uncomment following three lines for SC-p interactions
2418 c         xj=c(1,nres+j)-xi
2419 c         yj=c(2,nres+j)-yi
2420 c         zj=c(3,nres+j)-zi
2421 C Uncomment following three lines for Ca-p interactions
2422           xj=c(1,j)
2423           yj=c(2,j)
2424           zj=c(3,j)
2425       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2426       xj_safe=xj
2427       yj_safe=yj
2428       zj_safe=zj
2429       subchap=0
2430       do xshift=-1,1
2431       do yshift=-1,1
2432       do zshift=-1,1
2433           xj=xj_safe+xshift*boxxsize
2434           yj=yj_safe+yshift*boxysize
2435           zj=zj_safe+zshift*boxzsize
2436           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2437           if(dist_temp.lt.dist_init) then
2438             dist_init=dist_temp
2439             xj_temp=xj
2440             yj_temp=yj
2441             zj_temp=zj
2442             subchap=1
2443           endif
2444        enddo
2445        enddo
2446        enddo
2447        if (subchap.eq.1) then
2448           xj=xj_temp-xi
2449           yj=yj_temp-yi
2450           zj=zj_temp-zi
2451        else
2452           xj=xj_safe-xi
2453           yj=yj_safe-yi
2454           zj=zj_safe-zi
2455        endif
2456
2457           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2458
2459           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2460           sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2461           if (sss.lt.1.0d0) then
2462             fac=rrij**expon2
2463             e1=fac*fac*aad(itypj,iteli)
2464             e2=fac*bad(itypj,iteli)
2465             if (iabs(j-i) .le. 2) then
2466               e1=scal14*e1
2467               e2=scal14*e2
2468               evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2469             endif
2470             evdwij=e1+e2
2471             evdw2=evdw2+evdwij*(1.0d0-sss)
2472             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2473      &          'evdw2',i,j,sss,evdwij
2474 C
2475 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2476 C
2477              
2478             fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2479             fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2480             ggg(1)=xj*fac
2481             ggg(2)=yj*fac
2482             ggg(3)=zj*fac
2483 C Uncomment following three lines for SC-p interactions
2484 c           do k=1,3
2485 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2486 c           enddo
2487 C Uncomment following line for SC-p interactions
2488 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2489             do k=1,3
2490               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2491               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2492             enddo
2493           endif
2494         enddo
2495
2496         enddo ! iint
2497       enddo ! i
2498       do i=1,nct
2499         do j=1,3
2500           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2501           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2502           gradx_scp(j,i)=expon*gradx_scp(j,i)
2503         enddo
2504       enddo
2505 C******************************************************************************
2506 C
2507 C                              N O T E !!!
2508 C
2509 C To save time the factor EXPON has been extracted from ALL components
2510 C of GVDWC and GRADX. Remember to multiply them by this factor before further 
2511 C use!
2512 C
2513 C******************************************************************************
2514       return
2515       end
2516 C-----------------------------------------------------------------------------
2517       subroutine escp_short(evdw2,evdw2_14)
2518 C
2519 C This subroutine calculates the excluded-volume interaction energy between
2520 C peptide-group centers and side chains and its gradient in virtual-bond and
2521 C side-chain vectors.
2522 C
2523       implicit real*8 (a-h,o-z)
2524       include 'DIMENSIONS'
2525       include 'COMMON.GEO'
2526       include 'COMMON.VAR'
2527       include 'COMMON.LOCAL'
2528       include 'COMMON.CHAIN'
2529       include 'COMMON.DERIV'
2530       include 'COMMON.INTERACT'
2531       include 'COMMON.FFIELD'
2532       include 'COMMON.IOUNITS'
2533       include 'COMMON.CONTROL'
2534       dimension ggg(3)
2535       evdw2=0.0D0
2536       evdw2_14=0.0d0
2537 cd    print '(a)','Enter ESCP'
2538 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2539       do i=iatscp_s,iatscp_e
2540         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2541         iteli=itel(i)
2542         xi=0.5D0*(c(1,i)+c(1,i+1))
2543         yi=0.5D0*(c(2,i)+c(2,i+1))
2544         zi=0.5D0*(c(3,i)+c(3,i+1))
2545          xi=mod(xi,boxxsize)
2546           if (xi.lt.0) xi=xi+boxxsize
2547           yi=mod(yi,boxysize)
2548           if (yi.lt.0) yi=yi+boxysize
2549           zi=mod(zi,boxzsize)
2550           if (zi.lt.0) zi=zi+boxzsize
2551
2552         do iint=1,nscp_gr(i)
2553
2554         do j=iscpstart(i,iint),iscpend(i,iint)
2555           itypj=itype(j)
2556           if (itypj.eq.ntyp1) cycle
2557 C Uncomment following three lines for SC-p interactions
2558 c         xj=c(1,nres+j)-xi
2559 c         yj=c(2,nres+j)-yi
2560 c         zj=c(3,nres+j)-zi
2561 C Uncomment following three lines for Ca-p interactions
2562           xj=c(1,j)
2563           yj=c(2,j)
2564           zj=c(3,j)
2565       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2566       xj_safe=xj
2567       yj_safe=yj
2568       zj_safe=zj
2569       subchap=0
2570       do xshift=-1,1
2571       do yshift=-1,1
2572       do zshift=-1,1
2573           xj=xj_safe+xshift*boxxsize
2574           yj=yj_safe+yshift*boxysize
2575           zj=zj_safe+zshift*boxzsize
2576           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2577           if(dist_temp.lt.dist_init) then
2578             dist_init=dist_temp
2579             xj_temp=xj
2580             yj_temp=yj
2581             zj_temp=zj
2582             subchap=1
2583           endif
2584        enddo
2585        enddo
2586        enddo
2587        if (subchap.eq.1) then
2588           xj=xj_temp-xi
2589           yj=yj_temp-yi
2590           zj=zj_temp-zi
2591        else
2592           xj=xj_safe-xi
2593           yj=yj_safe-yi
2594           zj=zj_safe-zi
2595        endif
2596           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2597           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2598           sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2599           if (sss.gt.0.0d0) then
2600
2601             fac=rrij**expon2
2602             e1=fac*fac*aad(itypj,iteli)
2603             e2=fac*bad(itypj,iteli)
2604             if (iabs(j-i) .le. 2) then
2605               e1=scal14*e1
2606               e2=scal14*e2
2607               evdw2_14=evdw2_14+(e1+e2)*sss
2608             endif
2609             evdwij=e1+e2
2610             evdw2=evdw2+evdwij*sss
2611             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2612      &          'evdw2',i,j,sss,evdwij
2613 C
2614 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2615 C
2616             fac=-(evdwij+e1)*rrij*sss
2617             fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2618             ggg(1)=xj*fac
2619             ggg(2)=yj*fac
2620             ggg(3)=zj*fac
2621 C Uncomment following three lines for SC-p interactions
2622 c           do k=1,3
2623 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2624 c           enddo
2625 C Uncomment following line for SC-p interactions
2626 c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2627             do k=1,3
2628               gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2629               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2630             enddo
2631           endif
2632         enddo
2633
2634         enddo ! iint
2635       enddo ! i
2636       do i=1,nct
2637         do j=1,3
2638           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2639           gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2640           gradx_scp(j,i)=expon*gradx_scp(j,i)
2641         enddo
2642       enddo
2643 C******************************************************************************
2644 C
2645 C                              N O T E !!!
2646 C
2647 C To save time the factor EXPON has been extracted from ALL components
2648 C of GVDWC and GRADX. Remember to multiply them by this factor before further 
2649 C use!
2650 C
2651 C******************************************************************************
2652       return
2653       end