1 C-----------------------------------------------------------------------
2 double precision function sscalelip(r)
4 double precision r,gamm
5 include "COMMON.SPLITELE"
6 C if(r.lt.r_cut-rlamb) then
8 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
9 C gamm=(r-(r_cut-rlamb))/rlamb
10 sscalelip=1.0d0+r*r*(2*r-3.0d0)
16 C-----------------------------------------------------------------------
17 double precision function sscagradlip(r)
19 double precision r,gamm
20 include "COMMON.SPLITELE"
21 C if(r.lt.r_cut-rlamb) then
23 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
24 C gamm=(r-(r_cut-rlamb))/rlamb
25 sscagradlip=r*(6*r-6.0d0)
32 C-----------------------------------------------------------------------
33 double precision function sscale(r,r_cut)
35 double precision r,r_cut,gamm
36 include "COMMON.SPLITELE"
37 if(r.lt.r_cut-rlamb) then
39 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
40 gamm=(r-(r_cut-rlamb))/rlamb
41 sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
47 C-----------------------------------------------------------------------
48 double precision function sscagrad(r,r_cut)
50 double precision r,r_cut,gamm
51 include "COMMON.SPLITELE"
52 if(r.lt.r_cut-rlamb) then
54 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
55 gamm=(r-(r_cut-rlamb))/rlamb
56 sscagrad=gamm*(6*gamm-6.0d0)/rlamb
62 C-----------------------------------------------------------------------
63 subroutine elj_long(evdw)
65 C This subroutine calculates the interaction energy of nonbonded side chains
66 C assuming the LJ potential of interaction.
72 include 'COMMON.LOCAL'
73 include 'COMMON.CHAIN'
74 include 'COMMON.DERIV'
75 include 'COMMON.INTERACT'
76 include 'COMMON.TORSION'
77 include 'COMMON.SBRIDGE'
78 include 'COMMON.NAMES'
79 include 'COMMON.IOUNITS'
80 include "COMMON.SPLITELE"
81 c include 'COMMON.CONTACTS'
82 double precision gg(3)
83 double precision evdw,evdwij
84 integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont
85 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
86 & sigij,r0ij,rcut,sss1,sssgrad1,sqrij
87 double precision sscale,sscagrad
88 double precision boxshift
89 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
91 c do i=iatsc_s,iatsc_e
92 do ikont=g_listscsc_start,g_listscsc_end
96 if (itypi.eq.ntyp1) cycle
97 itypi1=iabs(itype(i+1))
101 call to_box(xi,yi,zi)
103 C Calculate SC interaction energy.
105 c do iint=1,nint_gr(i)
106 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
107 cd & 'iend=',iend(i,iint)
108 c do j=istart(i,iint),iend(i,iint)
110 if (itypj.eq.ntyp1) cycle
114 call to_box(xj,yj,zj)
115 xj=boxshift(xj-xi,boxxsize)
116 yj=boxshift(yj-yi,boxysize)
117 zj=boxshift(zj-zi,boxzsize)
118 rij=xj*xj+yj*yj+zj*zj
120 eps0ij=eps(itypi,itypj)
121 sss1=sscale(sqrij,r_cut_int)
122 if (sss1.eq.0.0d0) cycle
123 sssgrad1=sscagrad(sqrij,r_cut_int)
125 & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
126 sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
127 if (sss.lt.1.0d0) then
133 evdw=evdw+(1.0d0-sss)*sss1*evdwij/sqrij/expon
135 C Calculate the components of the gradient in DC and X
137 fac=-rrij*(e1+evdwij)*(1.0d0-sss)*sss1
138 & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
139 & +(1.0d0-sss)*sssgrad1)/sqrij
144 gvdwx(k,i)=gvdwx(k,i)-gg(k)
145 gvdwx(k,j)=gvdwx(k,j)+gg(k)
146 gvdwc(k,i)=gvdwc(k,i)-gg(k)
147 gvdwc(k,j)=gvdwc(k,j)+gg(k)
155 gvdwc(j,i)=expon*gvdwc(j,i)
156 gvdwx(j,i)=expon*gvdwx(j,i)
159 C******************************************************************************
163 C To save time, the factor of EXPON has been extracted from ALL components
164 C of GVDWC and GRADX. Remember to multiply them by this factor before further
167 C******************************************************************************
170 C-----------------------------------------------------------------------
171 subroutine elj_short(evdw)
173 C This subroutine calculates the interaction energy of nonbonded side chains
174 C assuming the LJ potential of interaction.
180 include 'COMMON.LOCAL'
181 include 'COMMON.CHAIN'
182 include 'COMMON.DERIV'
183 include 'COMMON.INTERACT'
184 include 'COMMON.TORSION'
185 include 'COMMON.SBRIDGE'
186 include 'COMMON.NAMES'
187 include 'COMMON.IOUNITS'
188 include "COMMON.SPLITELE"
189 c include 'COMMON.CONTACTS'
190 double precision gg(3)
191 double precision evdw,evdwij
192 integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont
193 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
194 & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
195 double precision sscale,sscagrad
196 double precision boxshift
197 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
199 c do i=iatsc_s,iatsc_e
200 do ikont=g_listscsc_start,g_listscsc_end
201 i=newcontlisti(ikont)
202 j=newcontlistj(ikont)
204 if (itypi.eq.ntyp1) cycle
205 itypi1=iabs(itype(i+1))
209 call to_box(xi,yi,zi)
213 C Calculate SC interaction energy.
215 c do iint=1,nint_gr(i)
216 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
217 cd & 'iend=',iend(i,iint)
218 c do j=istart(i,iint),iend(i,iint)
220 if (itypj.eq.ntyp1) cycle
224 call to_box(xj,yj,zj)
225 xj=boxshift(xj-xi,boxxsize)
226 yj=boxshift(yj-yi,boxysize)
227 zj=boxshift(zj-zi,boxzsize)
228 C Change 12/1/95 to calculate four-body interactions
229 rij=xj*xj+yj*yj+zj*zj
231 sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
232 if (sss.gt.0.0d0) then
234 & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
236 eps0ij=eps(itypi,itypj)
243 C Calculate the components of the gradient in DC and X
245 fac=-rrij*(e1+evdwij)*sss+evdwij*sssgrad/sqrij/expon
250 gvdwx(k,i)=gvdwx(k,i)-gg(k)
251 gvdwx(k,j)=gvdwx(k,j)+gg(k)
252 gvdwc(k,i)=gvdwc(k,i)-gg(k)
253 gvdwc(k,j)=gvdwc(k,j)+gg(k)
261 gvdwc(j,i)=expon*gvdwc(j,i)
262 gvdwx(j,i)=expon*gvdwx(j,i)
265 C******************************************************************************
269 C To save time, the factor of EXPON has been extracted from ALL components
270 C of GVDWC and GRADX. Remember to multiply them by this factor before further
273 C******************************************************************************
276 C-----------------------------------------------------------------------------
277 subroutine eljk_long(evdw)
279 C This subroutine calculates the interaction energy of nonbonded side chains
280 C assuming the LJK potential of interaction.
286 include 'COMMON.LOCAL'
287 include 'COMMON.CHAIN'
288 include 'COMMON.DERIV'
289 include 'COMMON.INTERACT'
290 include 'COMMON.IOUNITS'
291 include 'COMMON.NAMES'
292 include "COMMON.SPLITELE"
293 double precision gg(3)
294 double precision evdw,evdwij
295 integer i,j,k,itypi,itypj,itypi1,iint,ikont
296 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
297 & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
299 double precision sscale,sscagrad
300 double precision boxshift
301 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
303 c do i=iatsc_s,iatsc_e
304 do ikont=g_listscsc_start,g_listscsc_end
305 i=newcontlisti(ikont)
306 j=newcontlistj(ikont)
308 if (itypi.eq.ntyp1) cycle
309 itypi1=iabs(itype(i+1))
313 call to_box(xi,yi,zi)
315 C Calculate SC interaction energy.
317 c do iint=1,nint_gr(i)
318 c do j=istart(i,iint),iend(i,iint)
320 if (itypj.eq.ntyp1) cycle
324 call to_box(xj,yj,zj)
325 xj=boxshift(xj-xi,boxxsize)
326 yj=boxshift(yj-yi,boxysize)
327 zj=boxshift(zj-zi,boxzsize)
328 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
330 e_augm=augm(itypi,itypj)*fac_augm
333 sss1=sscale(rij,r_cut_int)
334 if (sss1.eq.0.0d0) cycle
335 sssgrad1=sscagrad(rij,r_cut_int)
336 sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
337 if (sss.lt.1.0d0) then
339 & sscagrad(rij/sigma(itypi,itypj),r_cut_respa)
340 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
341 fac=r_shift_inv**expon
345 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
346 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
347 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
348 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
349 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
350 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
351 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
352 evdw=evdw+(1.0d0-sss)*sss1*evdwij
354 C Calculate the components of the gradient in DC and X
356 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
357 fac=fac*(1.0d0-sss)*sss1
358 & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
359 & +(1.0d0-sss)*sssgrad1)*r_inv_ij/expon
364 gvdwx(k,i)=gvdwx(k,i)-gg(k)
365 gvdwx(k,j)=gvdwx(k,j)+gg(k)
366 gvdwc(k,i)=gvdwc(k,i)-gg(k)
367 gvdwc(k,j)=gvdwc(k,j)+gg(k)
375 gvdwc(j,i)=expon*gvdwc(j,i)
376 gvdwx(j,i)=expon*gvdwx(j,i)
381 C-----------------------------------------------------------------------------
382 subroutine eljk_short(evdw)
384 C This subroutine calculates the interaction energy of nonbonded side chains
385 C assuming the LJK potential of interaction.
387 implicit real*8 (a-h,o-z)
391 include 'COMMON.LOCAL'
392 include 'COMMON.CHAIN'
393 include 'COMMON.DERIV'
394 include 'COMMON.INTERACT'
395 include 'COMMON.IOUNITS'
396 include 'COMMON.NAMES'
397 include "COMMON.SPLITELE"
398 double precision gg(3)
399 double precision evdw,evdwij
400 integer i,j,k,itypi,itypj,itypi1,iint,ikont
401 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
402 & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
404 double precision sscale,sscagrad
405 double precision boxshift
406 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
408 c do i=iatsc_s,iatsc_e
409 do ikont=g_listscsc_start,g_listscsc_end
410 i=newcontlisti(ikont)
411 j=newcontlistj(ikont)
413 if (itypi.eq.ntyp1) cycle
414 itypi1=iabs(itype(i+1))
418 call to_box(xi,yi,zi)
420 C Calculate SC interaction energy.
422 c do iint=1,nint_gr(i)
423 c do j=istart(i,iint),iend(i,iint)
425 if (itypj.eq.ntyp1) cycle
429 call to_box(xj,yj,zj)
430 xj=boxshift(xj-xi,boxxsize)
431 yj=boxshift(yj-yi,boxysize)
432 zj=boxshift(zj-zi,boxzsize)
433 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
435 e_augm=augm(itypi,itypj)*fac_augm
438 sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
439 if (sss.gt.0.0d0) then
440 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
441 fac=r_shift_inv**expon
445 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
446 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
447 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
448 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
449 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
450 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
451 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
454 C Calculate the components of the gradient in DC and X
456 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
457 & +evdwij*sssgrad/sigma(itypi,itypj)*r_inv_ij/expon
463 gvdwx(k,i)=gvdwx(k,i)-gg(k)
464 gvdwx(k,j)=gvdwx(k,j)+gg(k)
465 gvdwc(k,i)=gvdwc(k,i)-gg(k)
466 gvdwc(k,j)=gvdwc(k,j)+gg(k)
474 gvdwc(j,i)=expon*gvdwc(j,i)
475 gvdwx(j,i)=expon*gvdwx(j,i)
480 C-----------------------------------------------------------------------------
481 subroutine ebp_long(evdw)
483 C This subroutine calculates the interaction energy of nonbonded side chains
484 C assuming the Berne-Pechukas potential of interaction.
490 include 'COMMON.LOCAL'
491 include 'COMMON.CHAIN'
492 include 'COMMON.DERIV'
493 include 'COMMON.NAMES'
494 include 'COMMON.INTERACT'
495 include 'COMMON.IOUNITS'
496 include 'COMMON.CALC'
497 include "COMMON.SPLITELE"
500 double precision evdw
501 integer itypi,itypj,itypi1,iint,ind,ikont
502 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
503 double precision sss1,sssgrad1
504 double precision sscale,sscagrad
505 double precision boxshift
506 c double precision rrsave(maxdim)
509 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
511 c if (icall.eq.0) then
517 c do i=iatsc_s,iatsc_e
518 do ikont=g_listscsc_start,g_listscsc_end
519 i=newcontlisti(ikont)
520 j=newcontlistj(ikont)
522 if (itypi.eq.ntyp1) cycle
523 itypi1=iabs(itype(i+1))
527 call to_box(xi,yi,zi)
528 dxi=dc_norm(1,nres+i)
529 dyi=dc_norm(2,nres+i)
530 dzi=dc_norm(3,nres+i)
531 c dsci_inv=dsc_inv(itypi)
532 dsci_inv=vbld_inv(i+nres)
534 C Calculate SC interaction energy.
536 c do iint=1,nint_gr(i)
537 c do j=istart(i,iint),iend(i,iint)
540 if (itypj.eq.ntyp1) cycle
541 c dscj_inv=dsc_inv(itypj)
542 dscj_inv=vbld_inv(j+nres)
543 chi1=chi(itypi,itypj)
544 chi2=chi(itypj,itypi)
551 alf12=0.5D0*(alf1+alf2)
555 call to_box(xj,yj,zj)
556 xj=boxshift(xj-xi,boxxsize)
557 yj=boxshift(yj-yi,boxysize)
558 zj=boxshift(zj-zi,boxzsize)
559 dxj=dc_norm(1,nres+j)
560 dyj=dc_norm(2,nres+j)
561 dzj=dc_norm(3,nres+j)
562 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
564 sss1=sscale(1.0d0/rij,r_cut_int)
565 if (sss1.eq.0.0d0) cycle
566 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
567 if (sss.lt.1.0d0) then
569 & sscagrad(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
570 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
571 C Calculate the angle-dependent terms of energy & contributions to derivatives.
573 C Calculate whole angle-dependent part of epsilon and contributions
575 fac=(rrij*sigsq)**expon2
578 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
579 eps2der=evdwij*eps3rt
580 eps3der=evdwij*eps2rt
581 evdwij=evdwij*eps2rt*eps3rt
582 evdw=evdw+evdwij*(1.0d0-sss)*sss1
584 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
586 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
587 cd & restyp(itypi),i,restyp(itypj),j,
588 cd & epsi,sigm,chi1,chi2,chip1,chip2,
589 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
590 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
593 C Calculate gradient components.
594 e1=e1*eps1*eps2rt**2*eps3rt**2
595 fac=-expon*(e1+evdwij)
597 fac=(fac+evdwij*(sss1/(1.0d0-sss)*sssgrad/
598 & sigmaii(itypi,itypj)+(1.0d0-sss)/sss1*sssgrad1))*rij
599 C Calculate radial part of the gradient
603 C Calculate the angular part of the gradient and sum add the contributions
604 C to the appropriate components of the Cartesian gradient.
605 call sc_grad_scale((1.0d0-sss)*sss1)
613 C-----------------------------------------------------------------------------
614 subroutine ebp_short(evdw)
616 C This subroutine calculates the interaction energy of nonbonded side chains
617 C assuming the Berne-Pechukas potential of interaction.
619 implicit real*8 (a-h,o-z)
623 include 'COMMON.LOCAL'
624 include 'COMMON.CHAIN'
625 include 'COMMON.DERIV'
626 include 'COMMON.NAMES'
627 include 'COMMON.INTERACT'
628 include 'COMMON.IOUNITS'
629 include 'COMMON.CALC'
630 include "COMMON.SPLITELE"
633 double precision evdw
634 integer itypi,itypj,itypi1,iint,ind,ikont
635 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
636 double precision sscale,sscagrad
637 double precision boxshift
638 c double precision rrsave(maxdim)
641 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
642 c if (icall.eq.0) then
648 c do i=iatsc_s,iatsc_e
649 do ikont=g_listscsc_start,g_listscsc_end
650 i=newcontlisti(ikont)
651 j=newcontlistj(ikont)
653 if (itypi.eq.ntyp1) cycle
654 itypi1=iabs(itype(i+1))
658 call to_box(xi,yi,zi)
659 dxi=dc_norm(1,nres+i)
660 dyi=dc_norm(2,nres+i)
661 dzi=dc_norm(3,nres+i)
662 c dsci_inv=dsc_inv(itypi)
663 dsci_inv=vbld_inv(i+nres)
665 C Calculate SC interaction energy.
667 c do iint=1,nint_gr(i)
668 c do j=istart(i,iint),iend(i,iint)
671 if (itypj.eq.ntyp1) cycle
672 c dscj_inv=dsc_inv(itypj)
673 dscj_inv=vbld_inv(j+nres)
674 chi1=chi(itypi,itypj)
675 chi2=chi(itypj,itypi)
682 alf12=0.5D0*(alf1+alf2)
686 call to_box(xj,yj,zj)
687 xj=boxshift(xj-xi,boxxsize)
688 yj=boxshift(yj-yi,boxysize)
689 zj=boxshift(zj-zi,boxzsize)
690 dxj=dc_norm(1,nres+j)
691 dyj=dc_norm(2,nres+j)
692 dzj=dc_norm(3,nres+j)
693 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
695 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
697 if (sss.gt.0.0d0) then
699 C Calculate the angle-dependent terms of energy & contributions to derivatives.
701 C Calculate whole angle-dependent part of epsilon and contributions
703 fac=(rrij*sigsq)**expon2
706 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
707 eps2der=evdwij*eps3rt
708 eps3der=evdwij*eps2rt
709 evdwij=evdwij*eps2rt*eps3rt
712 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
714 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
715 cd & restyp(itypi),i,restyp(itypj),j,
716 cd & epsi,sigm,chi1,chi2,chip1,chip2,
717 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
718 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
721 C Calculate gradient components.
722 e1=e1*eps1*eps2rt**2*eps3rt**2
723 fac=-expon*(e1+evdwij)
725 fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rrij
726 C Calculate radial part of the gradient
730 C Calculate the angular part of the gradient and sum add the contributions
731 C to the appropriate components of the Cartesian gradient.
732 call sc_grad_scale(sss)
740 C-----------------------------------------------------------------------------
741 subroutine egb_long(evdw)
743 C This subroutine calculates the interaction energy of nonbonded side chains
744 C assuming the Gay-Berne potential of interaction.
750 include 'COMMON.LOCAL'
751 include 'COMMON.CHAIN'
752 include 'COMMON.DERIV'
753 include 'COMMON.NAMES'
754 include 'COMMON.INTERACT'
755 include 'COMMON.IOUNITS'
756 include 'COMMON.CALC'
757 include 'COMMON.CONTROL'
758 include "COMMON.SPLITELE"
760 integer xshift,yshift,zshift
761 double precision evdw
762 integer itypi,itypj,itypi1,iint,ind,ikont
763 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
764 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
765 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
766 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
767 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
768 double precision subchap,sss1,sssgrad1
769 double precision boxshift
771 ccccc energy_dec=.false.
772 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
775 c if (icall.eq.0) lprn=.false.
777 c do i=iatsc_s,iatsc_e
778 do ikont=g_listscsc_start,g_listscsc_end
779 i=newcontlisti(ikont)
780 j=newcontlistj(ikont)
782 if (itypi.eq.ntyp1) cycle
783 itypi1=iabs(itype(i+1))
787 call to_box(xi,yi,zi)
788 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
789 dxi=dc_norm(1,nres+i)
790 dyi=dc_norm(2,nres+i)
791 dzi=dc_norm(3,nres+i)
792 c dsci_inv=dsc_inv(itypi)
793 dsci_inv=vbld_inv(i+nres)
794 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
795 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
797 C Calculate SC interaction energy.
799 c do iint=1,nint_gr(i)
800 c do j=istart(i,iint),iend(i,iint)
803 if (itypj.eq.ntyp1) cycle
804 c dscj_inv=dsc_inv(itypj)
805 dscj_inv=vbld_inv(j+nres)
806 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
807 c & 1.0d0/vbld(j+nres)
808 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
809 sig0ij=sigma(itypi,itypj)
810 chi1=chi(itypi,itypj)
811 chi2=chi(itypj,itypi)
818 alf12=0.5D0*(alf1+alf2)
822 call to_box(xj,yj,zj)
823 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
824 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
825 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
826 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
827 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
828 xj=boxshift(xj-xi,boxxsize)
829 yj=boxshift(yj-yi,boxysize)
830 zj=boxshift(zj-zi,boxzsize)
831 dxj=dc_norm(1,nres+j)
832 dyj=dc_norm(2,nres+j)
833 dzj=dc_norm(3,nres+j)
834 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
836 sss1=sscale(1.0d0/rij,r_cut_int)
837 if (sss1.eq.0.0d0) cycle
838 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
839 if (sss.lt.1.0d0) then
840 C Calculate angle-dependent terms of energy and contributions to their
843 & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
844 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
847 sig=sig0ij*dsqrt(sigsq)
848 rij_shift=1.0D0/rij-sig+sig0ij
849 c for diagnostics; uncomment
850 c rij_shift=1.2*sig0ij
851 C I hate to put IF's in the loops, but here don't have another choice!!!!
852 if (rij_shift.le.0.0D0) then
854 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
855 cd & restyp(itypi),i,restyp(itypj),j,
856 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
860 c---------------------------------------------------------------
861 rij_shift=1.0D0/rij_shift
865 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
866 eps2der=evdwij*eps3rt
867 eps3der=evdwij*eps2rt
868 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
869 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
870 evdwij=evdwij*eps2rt*eps3rt
871 evdw=evdw+evdwij*(1.0d0-sss)*sss1
873 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
875 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
876 & restyp(itypi),i,restyp(itypj),j,
877 & epsi,sigm,chi1,chi2,chip1,chip2,
878 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
879 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
883 if (energy_dec) write (iout,'(a6,2i5,4f10.5)')
884 & 'evdw',i,j,rij,sss,sss1,evdwij
886 C Calculate gradient components.
887 e1=e1*eps1*eps2rt**2*eps3rt**2
888 fac=-expon*(e1+evdwij)*rij_shift
890 fac=(fac+evdwij*(-sss1*sssgrad/(1.0d0-sss)
891 & /sigmaii(itypi,itypj)+(1.0d0-sss)*sssgrad1/sss1))*rij
893 C Calculate the radial part of the gradient
897 gg_lipi(3)=ssgradlipi*evdwij
898 gg_lipj(3)=ssgradlipj*evdwij
899 C Calculate angular part of the gradient.
900 call sc_grad_scale((1.0d0-sss)*sss1)
905 c write (iout,*) "Number of loop steps in EGB:",ind
906 cccc energy_dec=.false.
909 C-----------------------------------------------------------------------------
910 subroutine egb_short(evdw)
912 C This subroutine calculates the interaction energy of nonbonded side chains
913 C assuming the Gay-Berne potential of interaction.
919 include 'COMMON.LOCAL'
920 include 'COMMON.CHAIN'
921 include 'COMMON.DERIV'
922 include 'COMMON.NAMES'
923 include 'COMMON.INTERACT'
924 include 'COMMON.IOUNITS'
925 include 'COMMON.CALC'
926 include 'COMMON.CONTROL'
927 include "COMMON.SPLITELE"
929 double precision evdw
930 integer itypi,itypj,itypi1,iint,ind,ikont
931 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
932 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
933 & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
934 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
935 double precision boxshift
937 ccccc energy_dec=.false.
938 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
941 c if (icall.eq.0) lprn=.false.
943 c do i=iatsc_s,iatsc_e
944 do ikont=g_listscsc_start,g_listscsc_end
945 i=newcontlisti(ikont)
946 j=newcontlistj(ikont)
948 if (itypi.eq.ntyp1) cycle
949 itypi1=iabs(itype(i+1))
953 call to_box(xi,yi,zi)
954 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
955 dxi=dc_norm(1,nres+i)
956 dyi=dc_norm(2,nres+i)
957 dzi=dc_norm(3,nres+i)
958 c dsci_inv=dsc_inv(itypi)
959 dsci_inv=vbld_inv(i+nres)
960 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
961 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
963 C Calculate SC interaction energy.
965 c do iint=1,nint_gr(i)
966 c do j=istart(i,iint),iend(i,iint)
969 if (itypj.eq.ntyp1) cycle
970 c dscj_inv=dsc_inv(itypj)
971 dscj_inv=vbld_inv(j+nres)
972 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
973 c & 1.0d0/vbld(j+nres)
974 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
975 sig0ij=sigma(itypi,itypj)
976 chi1=chi(itypi,itypj)
977 chi2=chi(itypj,itypi)
984 alf12=0.5D0*(alf1+alf2)
988 call to_box(xj,yj,zj)
989 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
990 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
991 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
992 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
993 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
994 xj=boxshift(xj-xi,boxxsize)
995 yj=boxshift(yj-yi,boxysize)
996 zj=boxshift(zj-zi,boxzsize)
997 dxj=dc_norm(1,nres+j)
998 dyj=dc_norm(2,nres+j)
999 dzj=dc_norm(3,nres+j)
1000 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1002 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1003 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1004 if (sss.gt.0.0d0) then
1006 C Calculate angle-dependent terms of energy and contributions to their
1010 sig=sig0ij*dsqrt(sigsq)
1011 rij_shift=1.0D0/rij-sig+sig0ij
1012 c for diagnostics; uncomment
1013 c rij_shift=1.2*sig0ij
1014 C I hate to put IF's in the loops, but here don't have another choice!!!!
1015 if (rij_shift.le.0.0D0) then
1017 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1018 cd & restyp(itypi),i,restyp(itypj),j,
1019 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1023 c---------------------------------------------------------------
1024 rij_shift=1.0D0/rij_shift
1025 fac=rij_shift**expon
1028 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1029 eps2der=evdwij*eps3rt
1030 eps3der=evdwij*eps2rt
1031 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1032 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1033 evdwij=evdwij*eps2rt*eps3rt
1034 evdw=evdw+evdwij*sss
1036 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1038 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1039 & restyp(itypi),i,restyp(itypj),j,
1040 & epsi,sigm,chi1,chi2,chip1,chip2,
1041 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1042 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1046 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1049 C Calculate gradient components.
1050 e1=e1*eps1*eps2rt**2*eps3rt**2
1051 fac=-expon*(e1+evdwij)*rij_shift
1053 fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rij
1055 C Calculate the radial part of the gradient
1059 gg_lipi(3)=ssgradlipi*evdwij
1060 gg_lipj(3)=ssgradlipj*evdwij
1061 C Calculate angular part of the gradient.
1062 call sc_grad_scale(sss)
1067 c write (iout,*) "Number of loop steps in EGB:",ind
1068 cccc energy_dec=.false.
1071 C-----------------------------------------------------------------------------
1072 subroutine egbv_long(evdw)
1074 C This subroutine calculates the interaction energy of nonbonded side chains
1075 C assuming the Gay-Berne-Vorobjev potential of interaction.
1077 implicit real*8 (a-h,o-z)
1078 include 'DIMENSIONS'
1079 include 'COMMON.GEO'
1080 include 'COMMON.VAR'
1081 include 'COMMON.LOCAL'
1082 include 'COMMON.CHAIN'
1083 include 'COMMON.DERIV'
1084 include 'COMMON.NAMES'
1085 include 'COMMON.INTERACT'
1086 include 'COMMON.IOUNITS'
1087 include 'COMMON.CALC'
1088 include "COMMON.SPLITELE"
1090 common /srutu/ icall
1092 integer itypi,itypj,itypi1,iint,ind,ikont
1093 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1094 & xi,yi,zi,fac_augm,e_augm
1095 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1096 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1097 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1098 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1099 double precision sss1,sssgrad1
1101 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1103 c if (icall.eq.0) lprn=.true.
1105 c do i=iatsc_s,iatsc_e
1106 do ikont=g_listscsc_start,g_listscsc_end
1107 i=newcontlisti(ikont)
1108 j=newcontlistj(ikont)
1109 itypi=iabs(itype(i))
1110 if (itypi.eq.ntyp1) cycle
1111 itypi1=iabs(itype(i+1))
1115 call to_box(xi,yi,zi)
1116 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
1117 dxi=dc_norm(1,nres+i)
1118 dyi=dc_norm(2,nres+i)
1119 dzi=dc_norm(3,nres+i)
1120 c dsci_inv=dsc_inv(itypi)
1121 dsci_inv=vbld_inv(i+nres)
1123 C Calculate SC interaction energy.
1125 c do iint=1,nint_gr(i)
1126 c do j=istart(i,iint),iend(i,iint)
1128 itypj=iabs(itype(j))
1129 if (itypj.eq.ntyp1) cycle
1130 c dscj_inv=dsc_inv(itypj)
1131 dscj_inv=vbld_inv(j+nres)
1132 sig0ij=sigma(itypi,itypj)
1133 r0ij=r0(itypi,itypj)
1134 chi1=chi(itypi,itypj)
1135 chi2=chi(itypj,itypi)
1142 alf12=0.5D0*(alf1+alf2)
1146 call to_box(xj,yj,zj)
1147 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
1148 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1149 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1150 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1151 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1152 xj=boxshift(xj-xi,boxxsize)
1153 yj=boxshift(yj-yi,boxysize)
1154 zj=boxshift(zj-zi,boxzsize)
1155 dxj=dc_norm(1,nres+j)
1156 dyj=dc_norm(2,nres+j)
1157 dzj=dc_norm(3,nres+j)
1158 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1161 sss1=sscale(1.0d0/rij,r_cut_int)
1162 if (sss1.eq.0.0d0) cycle
1164 if (sss.lt.1.0d0) then
1166 & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1167 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
1169 C Calculate angle-dependent terms of energy and contributions to their
1173 sig=sig0ij*dsqrt(sigsq)
1174 rij_shift=1.0D0/rij-sig+r0ij
1175 C I hate to put IF's in the loops, but here don't have another choice!!!!
1176 if (rij_shift.le.0.0D0) then
1181 c---------------------------------------------------------------
1182 rij_shift=1.0D0/rij_shift
1183 fac=rij_shift**expon
1186 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1187 eps2der=evdwij*eps3rt
1188 eps3der=evdwij*eps2rt
1189 fac_augm=rrij**expon
1190 e_augm=augm(itypi,itypj)*fac_augm
1191 evdwij=evdwij*eps2rt*eps3rt
1192 evdw=evdw+(evdwij+e_augm)*sss1*(1.0d0-sss)
1194 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1196 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1197 & restyp(itypi),i,restyp(itypj),j,
1198 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1199 & chi1,chi2,chip1,chip2,
1200 & eps1,eps2rt**2,eps3rt**2,
1201 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1204 C Calculate gradient components.
1205 e1=e1*eps1*eps2rt**2*eps3rt**2
1206 fac=-expon*(e1+evdwij)*rij_shift
1208 fac=rij*fac-2*expon*rrij*e_augm
1209 fac=fac+(evdwij+e_augm)*
1210 & (-sss1*sssgrad/(1.0d0-sss)/sigmaii(itypi,itypj)
1211 & +(1.0d0-sss)*sssgrad1/sss1)*rij
1212 C Calculate the radial part of the gradient
1216 C Calculate angular part of the gradient.
1217 call sc_grad_scale((1.0d0-sss)*sss1)
1223 C-----------------------------------------------------------------------------
1224 subroutine egbv_short(evdw)
1226 C This subroutine calculates the interaction energy of nonbonded side chains
1227 C assuming the Gay-Berne-Vorobjev potential of interaction.
1230 include 'DIMENSIONS'
1231 include 'COMMON.GEO'
1232 include 'COMMON.VAR'
1233 include 'COMMON.LOCAL'
1234 include 'COMMON.CHAIN'
1235 include 'COMMON.DERIV'
1236 include 'COMMON.NAMES'
1237 include 'COMMON.INTERACT'
1238 include 'COMMON.IOUNITS'
1239 include 'COMMON.CALC'
1240 include "COMMON.SPLITELE"
1242 common /srutu/ icall
1244 integer itypi,itypj,itypi1,iint,ind,ikont
1245 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1246 & xi,yi,zi,fac_augm,e_augm
1247 double precision evdw
1248 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1249 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1250 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1251 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1252 double precision boxshift
1254 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1256 c if (icall.eq.0) lprn=.true.
1258 c do i=iatsc_s,iatsc_e
1259 do ikont=g_listscsc_start,g_listscsc_end
1260 i=newcontlisti(ikont)
1261 j=newcontlistj(ikont)
1262 itypi=iabs(itype(i))
1263 if (itypi.eq.ntyp1) cycle
1264 itypi1=iabs(itype(i+1))
1268 dxi=dc_norm(1,nres+i)
1269 dyi=dc_norm(2,nres+i)
1270 dzi=dc_norm(3,nres+i)
1271 call to_box(xi,yi,zi)
1272 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
1273 c dsci_inv=dsc_inv(itypi)
1274 dsci_inv=vbld_inv(i+nres)
1276 C Calculate SC interaction energy.
1278 c do iint=1,nint_gr(i)
1279 c do j=istart(i,iint),iend(i,iint)
1281 itypj=iabs(itype(j))
1282 if (itypj.eq.ntyp1) cycle
1283 c dscj_inv=dsc_inv(itypj)
1284 dscj_inv=vbld_inv(j+nres)
1285 sig0ij=sigma(itypi,itypj)
1286 r0ij=r0(itypi,itypj)
1287 chi1=chi(itypi,itypj)
1288 chi2=chi(itypj,itypi)
1295 alf12=0.5D0*(alf1+alf2)
1299 call to_box(xj,yj,zj)
1300 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
1301 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1302 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1303 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1304 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1305 xj=boxshift(xj-xi,boxxsize)
1306 yj=boxshift(yj-yi,boxysize)
1307 zj=boxshift(zj-zi,boxzsize)
1308 dxj=dc_norm(1,nres+j)
1309 dyj=dc_norm(2,nres+j)
1310 dzj=dc_norm(3,nres+j)
1311 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1314 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1316 if (sss.gt.0.0d0) then
1318 C Calculate angle-dependent terms of energy and contributions to their
1322 sig=sig0ij*dsqrt(sigsq)
1323 rij_shift=1.0D0/rij-sig+r0ij
1324 C I hate to put IF's in the loops, but here don't have another choice!!!!
1325 if (rij_shift.le.0.0D0) then
1330 c---------------------------------------------------------------
1331 rij_shift=1.0D0/rij_shift
1332 fac=rij_shift**expon
1335 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1336 eps2der=evdwij*eps3rt
1337 eps3der=evdwij*eps2rt
1338 fac_augm=rrij**expon
1339 e_augm=augm(itypi,itypj)*fac_augm
1340 evdwij=evdwij*eps2rt*eps3rt
1341 evdw=evdw+(evdwij+e_augm)*sss
1343 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1345 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1346 & restyp(itypi),i,restyp(itypj),j,
1347 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1348 & chi1,chi2,chip1,chip2,
1349 & eps1,eps2rt**2,eps3rt**2,
1350 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1353 C Calculate gradient components.
1354 e1=e1*eps1*eps2rt**2*eps3rt**2
1355 fac=-expon*(e1+evdwij)*rij_shift
1357 fac=rij*fac-2*expon*rrij*e_augm+
1358 & (evdwij+e_augm)*sssgrad/sigmaii(itypi,itypj)/sss*rij
1359 C Calculate the radial part of the gradient
1363 C Calculate angular part of the gradient.
1364 call sc_grad_scale(sss)
1370 C----------------------------------------------------------------------------
1371 subroutine sc_grad_scale(scalfac)
1372 implicit real*8 (a-h,o-z)
1373 include 'DIMENSIONS'
1374 include 'COMMON.CHAIN'
1375 include 'COMMON.DERIV'
1376 include 'COMMON.CALC'
1377 include 'COMMON.IOUNITS'
1378 include "COMMON.SPLITELE"
1379 double precision dcosom1(3),dcosom2(3)
1380 double precision scalfac
1381 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1382 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1383 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1384 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1388 c eom12=evdwij*eps1_om12
1390 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1391 c & " sigder",sigder
1392 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1393 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1395 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1396 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1399 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1401 c write (iout,*) "gg",(gg(k),k=1,3)
1403 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
1404 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1405 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1406 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
1407 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1408 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1409 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1410 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1411 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1412 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1415 C Calculate the components of the gradient in DC and X
1418 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
1419 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
1423 C--------------------------------------------------------------------------
1424 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1426 C This subroutine calculates the average interaction energy and its gradient
1427 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1428 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1429 C The potential depends both on the distance of peptide-group centers and on
1430 C the orientation of the CA-CA virtual bonds.
1432 implicit real*8 (a-h,o-z)
1436 include 'DIMENSIONS'
1437 include 'COMMON.CONTROL'
1438 include 'COMMON.SETUP'
1439 include 'COMMON.IOUNITS'
1440 include 'COMMON.GEO'
1441 include 'COMMON.VAR'
1442 include 'COMMON.LOCAL'
1443 include 'COMMON.CHAIN'
1444 include 'COMMON.DERIV'
1445 include 'COMMON.INTERACT'
1447 include 'COMMON.CONTACTS'
1448 include 'COMMON.CONTMAT'
1450 include 'COMMON.CORRMAT'
1451 include 'COMMON.TORSION'
1452 include 'COMMON.VECTORS'
1453 include 'COMMON.FFIELD'
1454 include 'COMMON.TIME1'
1455 include 'COMMON.SHIELD'
1456 include "COMMON.SPLITELE"
1458 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1459 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1460 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1461 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1462 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1463 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1465 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1467 double precision scal_el /1.0d0/
1469 double precision scal_el /0.5d0/
1472 C 13-go grudnia roku pamietnego...
1473 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1474 & 0.0d0,1.0d0,0.0d0,
1475 & 0.0d0,0.0d0,1.0d0/
1476 cd write(iout,*) 'In EELEC'
1478 cd write(iout,*) 'Type',i
1479 cd write(iout,*) 'B1',B1(:,i)
1480 cd write(iout,*) 'B2',B2(:,i)
1481 cd write(iout,*) 'CC',CC(:,:,i)
1482 cd write(iout,*) 'DD',DD(:,:,i)
1483 cd write(iout,*) 'EE',EE(:,:,i)
1485 cd call check_vecgrad
1487 C print *,"WCHODZE3"
1488 if (icheckgrad.eq.1) then
1490 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1492 dc_norm(k,i)=dc(k,i)*fac
1494 c write (iout,*) 'i',i,' fac',fac
1497 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1498 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1499 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1500 c call vec_and_deriv
1506 time_mat=time_mat+MPI_Wtime()-time01
1510 cd write (iout,*) 'i=',i
1512 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1515 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1516 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1531 cd print '(a)','Enter EELEC'
1532 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1534 gel_loc_loc(i)=0.0d0
1539 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1541 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1543 do i=iturn3_start,iturn3_end
1544 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1545 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1546 C & .or. itype(i-1).eq.ntyp1
1547 C & .or. itype(i+4).eq.ntyp1
1552 dx_normi=dc_norm(1,i)
1553 dy_normi=dc_norm(2,i)
1554 dz_normi=dc_norm(3,i)
1555 xmedi=c(1,i)+0.5d0*dxi
1556 ymedi=c(2,i)+0.5d0*dyi
1557 zmedi=c(3,i)+0.5d0*dzi
1558 call to_box(xmedi,ymedi,zmedi)
1560 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1561 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1563 num_cont_hb(i)=num_conti
1566 do i=iturn4_start,iturn4_end
1567 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1568 & .or. itype(i+3).eq.ntyp1
1569 & .or. itype(i+4).eq.ntyp1
1570 C & .or. itype(i+5).eq.ntyp1
1571 C & .or. itype(i-1).eq.ntyp1
1576 dx_normi=dc_norm(1,i)
1577 dy_normi=dc_norm(2,i)
1578 dz_normi=dc_norm(3,i)
1579 xmedi=c(1,i)+0.5d0*dxi
1580 ymedi=c(2,i)+0.5d0*dyi
1581 zmedi=c(3,i)+0.5d0*dzi
1582 call to_box(xmedi,ymedi,zmedi)
1584 num_conti=num_cont_hb(i)
1586 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1587 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1588 & call eturn4(i,eello_turn4)
1590 num_cont_hb(i)=num_conti
1594 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1596 c do i=iatel_s,iatel_e
1597 do ikont=g_listpp_start,g_listpp_end
1598 i=newcontlistppi(ikont)
1599 j=newcontlistppj(ikont)
1600 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1601 C & .or. itype(i+2).eq.ntyp1
1602 C & .or. itype(i-1).eq.ntyp1
1607 dx_normi=dc_norm(1,i)
1608 dy_normi=dc_norm(2,i)
1609 dz_normi=dc_norm(3,i)
1610 xmedi=c(1,i)+0.5d0*dxi
1611 ymedi=c(2,i)+0.5d0*dyi
1612 zmedi=c(3,i)+0.5d0*dzi
1613 call to_box(xmedi,ymedi,zmedi)
1614 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
1616 num_conti=num_cont_hb(i)
1618 c do j=ielstart(i),ielend(i)
1619 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1620 C & .or.itype(j+2).eq.ntyp1
1621 C & .or.itype(j-1).eq.ntyp1
1623 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1626 num_cont_hb(i)=num_conti
1629 c write (iout,*) "Number of loop steps in EELEC:",ind
1631 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1632 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1634 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1635 ccc eel_loc=eel_loc+eello_turn3
1636 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1639 C-------------------------------------------------------------------------------
1640 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1642 include 'DIMENSIONS'
1646 include 'COMMON.CONTROL'
1647 include 'COMMON.IOUNITS'
1648 include 'COMMON.GEO'
1649 include 'COMMON.VAR'
1650 include 'COMMON.LOCAL'
1651 include 'COMMON.CHAIN'
1652 include 'COMMON.DERIV'
1653 include 'COMMON.INTERACT'
1655 include 'COMMON.CONTACTS'
1656 include 'COMMON.CONTMAT'
1658 include 'COMMON.CORRMAT'
1659 include 'COMMON.TORSION'
1660 include 'COMMON.VECTORS'
1661 include 'COMMON.FFIELD'
1662 include 'COMMON.TIME1'
1663 include 'COMMON.SHIELD'
1664 include "COMMON.SPLITELE"
1665 integer xshift,yshift,zshift
1666 double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1667 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1668 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1669 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
1670 & gmuij2(4),gmuji2(4)
1671 integer j1,j2,num_conti
1672 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1673 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1675 integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ind,itypi,itypj
1676 integer ilist,iresshield
1677 double precision rlocshield
1678 double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
1679 double precision ees,evdw1,eel_loc,aaa,bbb,ael3i
1680 double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,
1681 & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,
1682 & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa,
1683 & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,
1684 & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,
1685 & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,
1686 & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
1687 & ecosgp,ecosam,ecosbm,ecosgm,ghalf,geel_loc_ij,geel_loc_ji,
1688 & dxi,dyi,dzi,a22,a23,a32,a33
1689 double precision dist_init,xmedi,ymedi,zmedi,xj_safe,yj_safe,
1690 & zj_safe,xj_temp,yj_temp,zj_temp,dist_temp,dx_normi,dy_normi,
1692 double precision sss1,sssgrad1
1693 double precision sscale,sscagrad
1694 double precision scalar
1695 double precision boxshift
1697 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1699 double precision scal_el /1.0d0/
1701 double precision scal_el /0.5d0/
1704 C 13-go grudnia roku pamietnego...
1705 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1706 & 0.0d0,1.0d0,0.0d0,
1707 & 0.0d0,0.0d0,1.0d0/
1708 c time00=MPI_Wtime()
1709 cd write (iout,*) "eelecij",i,j
1710 C print *,"WCHODZE2"
1714 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1715 aaa=app(iteli,itelj)
1716 bbb=bpp(iteli,itelj)
1717 ael6i=ael6(iteli,itelj)
1718 ael3i=ael3(iteli,itelj)
1722 dx_normj=dc_norm(1,j)
1723 dy_normj=dc_norm(2,j)
1724 dz_normj=dc_norm(3,j)
1728 call to_box(xj,yj,zj)
1729 xj=boxshift(xj-xmedi,boxxsize)
1730 yj=boxshift(yj-ymedi,boxysize)
1731 zj=boxshift(zj-zmedi,boxzsize)
1732 rij=xj*xj+yj*yj+zj*zj
1736 c For extracting the short-range part of Evdwpp
1737 sss1=sscale(rij,r_cut_int)
1738 if (sss1.eq.0.0d0) return
1739 sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
1740 sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
1741 sssgrad1=sscagrad(rij,r_cut_int)
1744 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1745 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1746 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1747 fac=cosa-3.0D0*cosb*cosg
1749 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1750 if (j.eq.i+2) ev1=scal_el*ev1
1755 if (shield_mode.eq.0) then
1759 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1761 el1=el1*fac_shield(i)**2*fac_shield(j)**2
1762 el2=el2*fac_shield(i)**2*fac_shield(j)**2
1764 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1765 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1767 evdw1=evdw1+evdwij*(1.0d0-sss)*sss1
1768 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1769 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1770 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1771 cd & xmedi,ymedi,zmedi,xj,yj,zj
1773 if (energy_dec) then
1774 write (iout,'(a6,2i5,0pf7.3,2f7.3)')
1775 & 'evdw1',i,j,evdwij,sss,sss1
1776 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1780 C Calculate contributions to the Cartesian gradient.
1783 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
1784 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1785 facel=-3*rrmij*(el1+eesij)*sss1
1786 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1792 * Radial derivatives. First process both termini of the fragment (i,j)
1794 aux=facel+sssgrad1*(1.0d0-sss)*eesij*rmij
1795 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1802 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1803 & (shield_mode.gt.0)) then
1805 do ilist=1,ishield_list(i)
1806 iresshield=shield_list(ilist,i)
1808 rlocshield=grad_shield_side(k,ilist,i)*eesij*sss1
1809 & /fac_shield(i)*2.0*sss1
1810 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1812 & +grad_shield_loc(k,ilist,i)*eesij*sss1/fac_shield(i)*2.0
1814 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1815 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1816 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1817 C if (iresshield.gt.i) then
1818 C do ishi=i+1,iresshield-1
1819 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1820 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1824 C do ishi=iresshield,i
1825 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1826 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1832 do ilist=1,ishield_list(j)
1833 iresshield=shield_list(ilist,j)
1835 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1837 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1839 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss1
1840 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1845 gshieldc(k,i)=gshieldc(k,i)+
1846 & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
1847 gshieldc(k,j)=gshieldc(k,j)+
1848 & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
1849 gshieldc(k,i-1)=gshieldc(k,i-1)+
1850 & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
1851 gshieldc(k,j-1)=gshieldc(k,j-1)+
1852 & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
1858 c ghalf=0.5D0*ggg(k)
1859 c gelc(k,i)=gelc(k,i)+ghalf
1860 c gelc(k,j)=gelc(k,j)+ghalf
1862 c 9/28/08 AL Gradient compotents will be summed only at the end
1864 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1865 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1867 c gelc_long(3,i)=gelc_long(3,i)+
1868 c ssgradlipi*eesij/2.0d0*lipscale**2*sss1
1870 * Loop over residues i+1 thru j-1.
1874 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1878 & (-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij
1879 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1884 c ghalf=0.5D0*ggg(k)
1885 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1886 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1888 c 9/28/08 AL Gradient compotents will be summed only at the end
1890 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1891 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1894 * Loop over residues i+1 thru j-1.
1898 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1902 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
1903 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1904 facel=-3*rrmij*(el1+eesij)*sss1
1905 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1907 c facvdw=ev1+evdwij*(1.0d0-sss)*sss1
1910 fac=-3*rrmij*(facvdw+facvdw+facel)
1915 * Radial derivatives. First process both termini of the fragment (i,j)
1917 aux=fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj))
1919 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1927 c ghalf=0.5D0*ggg(k)
1928 c gelc(k,i)=gelc(k,i)+ghalf
1929 c gelc(k,j)=gelc(k,j)+ghalf
1931 c 9/28/08 AL Gradient compotents will be summed only at the end
1933 gelc_long(k,j)=gelc(k,j)+ggg(k)
1934 gelc_long(k,i)=gelc(k,i)-ggg(k)
1937 * Loop over residues i+1 thru j-1.
1941 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1944 c 9/28/08 AL Gradient compotents will be summed only at the end
1949 & (-sssgrad*sss1/rpp(iteli,itelj)+sssgrad1*(1.0d0-sss))*rmij*evdwij
1954 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1955 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1961 ecosa=2.0D0*fac3*fac1+fac4
1964 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1965 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1967 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1968 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1970 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1971 cd & (dcosg(k),k=1,3)
1973 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1
1974 & *fac_shield(i)**2*fac_shield(j)**2
1975 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1979 c ghalf=0.5D0*ggg(k)
1980 c gelc(k,i)=gelc(k,i)+ghalf
1981 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1982 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1983 c gelc(k,j)=gelc(k,j)+ghalf
1984 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1985 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1989 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1994 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1995 & +ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss1
1996 & *fac_shield(i)**2*fac_shield(j)**2
1997 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2000 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2001 & +ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss1
2002 & *fac_shield(i)**2*fac_shield(j)**2
2003 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2004 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
2005 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
2007 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
2008 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
2009 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2011 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
2012 C energy of a peptide unit is assumed in the form of a second-order
2013 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
2014 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
2015 C are computed for EVERY pair of non-contiguous peptide groups.
2017 if (j.lt.nres-1) then
2028 muij(kkk)=mu(k,i)*mu(l,j)
2030 gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
2031 c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
2032 gmuij2(kkk)=gUb2(k,i)*mu(l,j)
2033 gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
2034 c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
2035 gmuji2(kkk)=mu(k,i)*gUb2(l,j)
2039 cd write (iout,*) 'EELEC: i',i,' j',j
2040 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
2041 cd write(iout,*) 'muij',muij
2042 ury=scalar(uy(1,i),erij)
2043 urz=scalar(uz(1,i),erij)
2044 vry=scalar(uy(1,j),erij)
2045 vrz=scalar(uz(1,j),erij)
2046 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
2047 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
2048 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
2049 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
2050 fac=dsqrt(-ael6i)*r3ij
2055 cd write (iout,'(4i5,4f10.5)')
2056 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
2057 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
2058 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
2059 cd & uy(:,j),uz(:,j)
2060 cd write (iout,'(4f10.5)')
2061 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
2062 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
2063 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
2064 cd write (iout,'(9f10.5/)')
2065 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
2066 C Derivatives of the elements of A in virtual-bond vectors
2067 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
2069 uryg(k,1)=scalar(erder(1,k),uy(1,i))
2070 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
2071 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
2072 urzg(k,1)=scalar(erder(1,k),uz(1,i))
2073 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
2074 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
2075 vryg(k,1)=scalar(erder(1,k),uy(1,j))
2076 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
2077 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
2078 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
2079 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
2080 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
2082 C Compute radial contributions to the gradient
2100 C Add the contributions coming from er
2103 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
2104 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
2105 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
2106 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
2109 C Derivatives in DC(i)
2110 cgrad ghalf1=0.5d0*agg(k,1)
2111 cgrad ghalf2=0.5d0*agg(k,2)
2112 cgrad ghalf3=0.5d0*agg(k,3)
2113 cgrad ghalf4=0.5d0*agg(k,4)
2114 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2115 & -3.0d0*uryg(k,2)*vry)!+ghalf1
2116 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2117 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
2118 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2119 & -3.0d0*urzg(k,2)*vry)!+ghalf3
2120 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2121 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
2122 C Derivatives in DC(i+1)
2123 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2124 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2125 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2126 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2127 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2128 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2129 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2130 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2131 C Derivatives in DC(j)
2132 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2133 & -3.0d0*vryg(k,2)*ury)!+ghalf1
2134 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2135 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
2136 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2137 & -3.0d0*vryg(k,2)*urz)!+ghalf3
2138 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2139 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
2140 C Derivatives in DC(j+1) or DC(nres-1)
2141 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2142 & -3.0d0*vryg(k,3)*ury)
2143 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2144 & -3.0d0*vrzg(k,3)*ury)
2145 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2146 & -3.0d0*vryg(k,3)*urz)
2147 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2148 & -3.0d0*vrzg(k,3)*urz)
2149 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
2151 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
2164 aggi(k,l)=-aggi(k,l)
2165 aggi1(k,l)=-aggi1(k,l)
2166 aggj(k,l)=-aggj(k,l)
2167 aggj1(k,l)=-aggj1(k,l)
2170 if (j.lt.nres-1) then
2176 aggi(k,l)=-aggi(k,l)
2177 aggi1(k,l)=-aggi1(k,l)
2178 aggj(k,l)=-aggj(k,l)
2179 aggj1(k,l)=-aggj1(k,l)
2190 aggi(k,l)=-aggi(k,l)
2191 aggi1(k,l)=-aggi1(k,l)
2192 aggj(k,l)=-aggj(k,l)
2193 aggj1(k,l)=-aggj1(k,l)
2198 IF (wel_loc.gt.0.0d0) THEN
2199 C Contribution to the local-electrostatic energy coming from the i-j pair
2200 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2202 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2204 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2205 & 'eelloc',i,j,eel_loc_ij
2208 if (shield_mode.eq.0) then
2215 eel_loc_ij=eel_loc_ij
2216 & *fac_shield(i)*fac_shield(j)*sss1
2217 eel_loc=eel_loc+eel_loc_ij
2219 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2220 & (shield_mode.gt.0)) then
2223 do ilist=1,ishield_list(i)
2224 iresshield=shield_list(ilist,i)
2226 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2229 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2231 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2232 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2236 do ilist=1,ishield_list(j)
2237 iresshield=shield_list(ilist,j)
2239 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2242 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2244 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2245 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2252 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2253 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2254 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2255 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2256 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2257 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2258 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2259 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2264 geel_loc_ij=(a22*gmuij1(1)
2268 & *fac_shield(i)*fac_shield(j)*sss1
2269 c write(iout,*) "derivative over thatai"
2270 c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
2272 gloc(nphi+i,icg)=gloc(nphi+i,icg)+
2273 & geel_loc_ij*wel_loc
2274 c write(iout,*) "derivative over thatai-1"
2275 c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
2282 gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
2283 & geel_loc_ij*wel_loc
2284 & *fac_shield(i)*fac_shield(j)*sss1
2286 c Derivative over j residue
2287 geel_loc_ji=a22*gmuji1(1)
2291 c write(iout,*) "derivative over thataj"
2292 c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
2295 gloc(nphi+j,icg)=gloc(nphi+j,icg)+
2296 & geel_loc_ji*wel_loc
2297 & *fac_shield(i)*fac_shield(j)*sss1
2304 c write(iout,*) "derivative over thataj-1"
2305 c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
2307 gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
2308 & geel_loc_ji*wel_loc
2309 & *fac_shield(i)*fac_shield(j)*sss1
2311 cC Paral derivatives in virtual-bond dihedral angles gamma
2313 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2314 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2315 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2316 & *fac_shield(i)*fac_shield(j)*sss1
2317 c & *fac_shield(i)*fac_shield(j)
2318 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2321 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2322 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2323 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2324 & *fac_shield(i)*fac_shield(j)*sss1
2325 c & *fac_shield(i)*fac_shield(j)
2326 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2328 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2329 aux=eel_loc_ij/sss1*sssgrad1*rmij
2334 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+
2335 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2336 & *fac_shield(i)*fac_shield(j)*sss1
2337 c & *fac_shield(i)*fac_shield(j)
2338 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2340 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2341 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2342 cgrad ghalf=0.5d0*ggg(l)
2343 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
2344 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
2348 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2351 C Remaining derivatives of eello
2352 c gel_loc_long(3,j)=gel_loc_long(3,j)+ &
2353 c ssgradlipj*eel_loc_ij/2.0d0*lipscale/ &
2354 c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2356 c gel_loc_long(3,i)=gel_loc_long(3,i)+ &
2357 c ssgradlipi*eel_loc_ij/2.0d0*lipscale/ &
2358 c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2361 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2362 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2363 & *fac_shield(i)*fac_shield(j)*sss1
2364 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2366 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2367 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2368 & *fac_shield(i)*fac_shield(j)*sss1
2369 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2371 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2372 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2373 & *fac_shield(i)*fac_shield(j)*sss1
2374 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2376 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2377 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2378 & *fac_shield(i)*fac_shield(j)*sss1
2379 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2384 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2385 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
2386 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2387 & .and. num_conti.le.maxconts) then
2388 c write (iout,*) i,j," entered corr"
2390 C Calculate the contact function. The ith column of the array JCONT will
2391 C contain the numbers of atoms that make contacts with the atom I (of numbers
2392 C greater than I). The arrays FACONT and GACONT will contain the values of
2393 C the contact function and its derivative.
2394 c r0ij=1.02D0*rpp(iteli,itelj)
2395 c r0ij=1.11D0*rpp(iteli,itelj)
2396 r0ij=2.20D0*rpp(iteli,itelj)
2397 c r0ij=1.55D0*rpp(iteli,itelj)
2398 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2399 if (fcont.gt.0.0D0) then
2400 num_conti=num_conti+1
2401 if (num_conti.gt.maxconts) then
2402 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2403 & ' will skip next contacts for this conf.'
2405 jcont_hb(num_conti,i)=j
2406 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
2407 cd " jcont_hb",jcont_hb(num_conti,i)
2408 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2409 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2410 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2412 d_cont(num_conti,i)=rij
2413 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2414 C --- Electrostatic-interaction matrix ---
2415 a_chuj(1,1,num_conti,i)=a22
2416 a_chuj(1,2,num_conti,i)=a23
2417 a_chuj(2,1,num_conti,i)=a32
2418 a_chuj(2,2,num_conti,i)=a33
2419 C --- Gradient of rij
2421 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2428 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2429 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2430 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2431 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2432 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2437 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2438 C Calculate contact energies
2440 wij=cosa-3.0D0*cosb*cosg
2443 c fac3=dsqrt(-ael6i)/r0ij**3
2444 fac3=dsqrt(-ael6i)*r3ij
2445 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2446 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2447 if (ees0tmp.gt.0) then
2448 ees0pij=dsqrt(ees0tmp)
2452 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2453 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2454 if (ees0tmp.gt.0) then
2455 ees0mij=dsqrt(ees0tmp)
2460 if (shield_mode.eq.0) then
2464 ees0plist(num_conti,i)=j
2465 C fac_shield(i)=0.4d0
2466 C fac_shield(j)=0.6d0
2468 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2469 & *fac_shield(i)*fac_shield(j)*sss1
2470 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2471 & *fac_shield(i)*fac_shield(j)*sss1
2473 C Diagnostics. Comment out or remove after debugging!
2474 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2475 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2476 c ees0m(num_conti,i)=0.0D0
2478 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2479 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2480 C Angular derivatives of the contact function
2481 ees0pij1=fac3/ees0pij
2482 ees0mij1=fac3/ees0mij
2483 fac3p=-3.0D0*fac3*rrmij
2484 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2485 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2487 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2488 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2489 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2490 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2491 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2492 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2493 ecosap=ecosa1+ecosa2
2494 ecosbp=ecosb1+ecosb2
2495 ecosgp=ecosg1+ecosg2
2496 ecosam=ecosa1-ecosa2
2497 ecosbm=ecosb1-ecosb2
2498 ecosgm=ecosg1-ecosg2
2507 facont_hb(num_conti,i)=fcont
2508 fprimcont=fprimcont/rij
2509 cd facont_hb(num_conti,i)=1.0D0
2510 C Following line is for diagnostics.
2513 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2514 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2517 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2518 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2520 gggp(1)=gggp(1)+ees0pijp*xj
2521 & +ees0p(num_conti,i)/sss1*rmij*xj*sssgrad1
2522 gggp(2)=gggp(2)+ees0pijp*yj
2523 & +ees0p(num_conti,i)/sss1*rmij*yj*sssgrad1
2524 gggp(3)=gggp(3)+ees0pijp*zj
2525 & +ees0p(num_conti,i)/sss1*rmij*zj*sssgrad1
2526 gggm(1)=gggm(1)+ees0mijp*xj
2527 & +ees0m(num_conti,i)/sss1*rmij*xj*sssgrad1
2528 gggm(2)=gggm(2)+ees0mijp*yj
2529 & +ees0m(num_conti,i)/sss1*rmij*yj*sssgrad1
2530 gggm(3)=gggm(3)+ees0mijp*zj
2531 & +ees0m(num_conti,i)/sss1*rmij*zj*sssgrad1
2532 C Derivatives due to the contact function
2533 gacont_hbr(1,num_conti,i)=fprimcont*xj
2534 gacont_hbr(2,num_conti,i)=fprimcont*yj
2535 gacont_hbr(3,num_conti,i)=fprimcont*zj
2538 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2539 c following the change of gradient-summation algorithm.
2541 cgrad ghalfp=0.5D0*gggp(k)
2542 cgrad ghalfm=0.5D0*gggm(k)
2543 gacontp_hb1(k,num_conti,i)=!ghalfp
2544 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2545 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2546 & *sss1*fac_shield(i)*fac_shield(j)
2548 gacontp_hb2(k,num_conti,i)=!ghalfp
2549 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2550 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2551 & *sss1*fac_shield(i)*fac_shield(j)
2553 gacontp_hb3(k,num_conti,i)=gggp(k)
2554 & *sss1*fac_shield(i)*fac_shield(j)
2556 gacontm_hb1(k,num_conti,i)=!ghalfm
2557 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2558 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2559 & *sss1*fac_shield(i)*fac_shield(j)
2561 gacontm_hb2(k,num_conti,i)=!ghalfm
2562 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2563 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2564 & *sss1*fac_shield(i)*fac_shield(j)
2566 gacontm_hb3(k,num_conti,i)=gggm(k)
2567 & *sss1*fac_shield(i)*fac_shield(j)
2571 endif ! num_conti.le.maxconts
2575 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2578 ghalf=0.5d0*agg(l,k)
2579 aggi(l,k)=aggi(l,k)+ghalf
2580 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2581 aggj(l,k)=aggj(l,k)+ghalf
2584 if (j.eq.nres-1 .and. i.lt.j-2) then
2587 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2592 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2595 C-----------------------------------------------------------------------
2596 subroutine evdwpp_short(evdw1)
2601 include 'DIMENSIONS'
2602 include 'COMMON.CONTROL'
2603 include 'COMMON.IOUNITS'
2604 include 'COMMON.GEO'
2605 include 'COMMON.VAR'
2606 include 'COMMON.LOCAL'
2607 include 'COMMON.CHAIN'
2608 include 'COMMON.DERIV'
2609 include 'COMMON.INTERACT'
2610 c include 'COMMON.CONTACTS'
2611 include 'COMMON.TORSION'
2612 include 'COMMON.VECTORS'
2613 include 'COMMON.FFIELD'
2614 include "COMMON.SPLITELE"
2615 double precision ggg(3)
2616 integer xshift,yshift,zshift
2617 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2619 double precision scal_el /1.0d0/
2621 double precision scal_el /0.5d0/
2623 c write (iout,*) "evdwpp_short"
2624 integer i,j,k,iteli,itelj,num_conti,ind,isubchap
2625 double precision dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb
2626 double precision xj,yj,zj,rij,rrmij,r3ij,r6ij,evdw1,
2627 & dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
2628 & dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw
2629 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
2630 & dist_temp, dist_init,sss_grad
2631 double precision sscale,sscagrad
2632 double precision boxshift
2636 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2637 c & " iatel_e_vdw",iatel_e_vdw
2639 c do i=iatel_s_vdw,iatel_e_vdw
2640 do ikont=g_listpp_vdw_start,g_listpp_vdw_end
2641 i=newcontlistpp_vdwi(ikont)
2642 j=newcontlistpp_vdwj(ikont)
2643 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2647 dx_normi=dc_norm(1,i)
2648 dy_normi=dc_norm(2,i)
2649 dz_normi=dc_norm(3,i)
2650 xmedi=c(1,i)+0.5d0*dxi
2651 ymedi=c(2,i)+0.5d0*dyi
2652 zmedi=c(3,i)+0.5d0*dzi
2653 call to_box(xmedi,ymedi,zmedi)
2655 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2656 c & ' ielend',ielend_vdw(i)
2658 c do j=ielstart_vdw(i),ielend_vdw(i)
2659 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2663 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2664 aaa=app(iteli,itelj)
2665 bbb=bpp(iteli,itelj)
2669 dx_normj=dc_norm(1,j)
2670 dy_normj=dc_norm(2,j)
2671 dz_normj=dc_norm(3,j)
2675 call to_box(xj,yj,zj)
2676 xj=boxshift(xj-xmedi,boxxsize)
2677 yj=boxshift(yj-ymedi,boxysize)
2678 zj=boxshift(zj-zmedi,boxzsize)
2679 rij=xj*xj+yj*yj+zj*zj
2682 c sss=sscale(rij/rpp(iteli,itelj))
2683 c sssgrad=sscagrad(rij/rpp(iteli,itelj))
2684 sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
2685 sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
2686 if (sss.gt.0.0d0) then
2691 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2692 if (j.eq.i+2) ev1=scal_el*ev1
2695 if (energy_dec) then
2696 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2698 evdw1=evdw1+evdwij*sss
2699 if (energy_dec) write (iout,'(a10,2i5,0pf7.3)')
2700 & 'evdw1_sum',i,j,evdw1
2702 C Calculate contributions to the Cartesian gradient.
2704 facvdw=-6*rrmij*(ev1+evdwij)*sss
2705 ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2706 ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2707 ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2712 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2713 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2720 C-----------------------------------------------------------------------------
2721 subroutine escp_long(evdw2,evdw2_14)
2723 C This subroutine calculates the excluded-volume interaction energy between
2724 C peptide-group centers and side chains and its gradient in virtual-bond and
2725 C side-chain vectors.
2727 implicit real*8 (a-h,o-z)
2728 include 'DIMENSIONS'
2729 include 'COMMON.GEO'
2730 include 'COMMON.VAR'
2731 include 'COMMON.LOCAL'
2732 include 'COMMON.CHAIN'
2733 include 'COMMON.DERIV'
2734 include 'COMMON.INTERACT'
2735 include 'COMMON.FFIELD'
2736 include 'COMMON.IOUNITS'
2737 include 'COMMON.CONTROL'
2738 include "COMMON.SPLITELE"
2739 logical lprint_short
2740 common /shortcheck/ lprint_short
2741 double precision ggg(3)
2742 integer i,iint,j,k,iteli,itypj,subchap
2743 double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
2745 double precision evdw2,evdw2_14,evdwij
2746 double precision sscale,sscagrad
2747 double precision boxshift
2749 if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
2752 CD print '(a)','Enter ESCP KURWA'
2753 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2755 c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
2756 c & ' iatscp_e=',iatscp_e
2757 c do i=iatscp_s,iatscp_e
2758 do ikont=g_listscp_start,g_listscp_end
2759 i=newcontlistscpi(ikont)
2760 j=newcontlistscpj(ikont)
2761 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2763 xi=0.5D0*(c(1,i)+c(1,i+1))
2764 yi=0.5D0*(c(2,i)+c(2,i+1))
2765 zi=0.5D0*(c(3,i)+c(3,i+1))
2766 call to_box(xi,yi,zi)
2768 c do iint=1,nscp_gr(i)
2770 c do j=iscpstart(i,iint),iscpend(i,iint)
2771 itypj=iabs(itype(j))
2772 if (itypj.eq.ntyp1) cycle
2773 C Uncomment following three lines for SC-p interactions
2777 C Uncomment following three lines for Ca-p interactions
2782 call to_box(xj,yj,zj)
2783 xj=boxshift(xj-xi,boxxsize)
2784 yj=boxshift(yj-yi,boxysize)
2785 zj=boxshift(zj-zi,boxzsize)
2787 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2789 sss1=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
2790 if (sss1.eq.0) cycle
2791 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
2793 & sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
2794 sssgrad1=sscagrad(1.0d0/dsqrt(rrij),r_cut_int)
2795 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2796 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2797 if (sss.lt.1.0d0) then
2799 e1=fac*fac*aad(itypj,iteli)
2800 e2=fac*bad(itypj,iteli)
2801 if (iabs(j-i) .le. 2) then
2804 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)*sss1
2807 evdw2=evdw2+evdwij*(1.0d0-sss)*sss1
2808 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2809 & 'evdw2',i,j,sss,evdwij
2811 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2814 fac=-(evdwij+e1)*rrij*(1.0d0-sss)*sss1
2815 fac=fac+evdwij*dsqrt(rrij)*(-sssgrad/rscp(itypj,iteli)
2820 C Uncomment following three lines for SC-p interactions
2822 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2824 C Uncomment following line for SC-p interactions
2825 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2827 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2828 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2837 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2838 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2839 gradx_scp(j,i)=expon*gradx_scp(j,i)
2842 C******************************************************************************
2846 C To save time the factor EXPON has been extracted from ALL components
2847 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2850 C******************************************************************************
2853 C-----------------------------------------------------------------------------
2854 subroutine escp_short(evdw2,evdw2_14)
2856 C This subroutine calculates the excluded-volume interaction energy between
2857 C peptide-group centers and side chains and its gradient in virtual-bond and
2858 C side-chain vectors.
2861 include 'DIMENSIONS'
2862 include 'COMMON.GEO'
2863 include 'COMMON.VAR'
2864 include 'COMMON.LOCAL'
2865 include 'COMMON.CHAIN'
2866 include 'COMMON.DERIV'
2867 include 'COMMON.INTERACT'
2868 include 'COMMON.FFIELD'
2869 include 'COMMON.IOUNITS'
2870 include 'COMMON.CONTROL'
2871 include "COMMON.SPLITELE"
2872 integer xshift,yshift,zshift
2873 logical lprint_short
2874 common /shortcheck/ lprint_short
2875 integer i,iint,j,k,iteli,itypj,subchap
2876 double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
2878 double precision evdw2,evdw2_14,evdwij
2879 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
2880 & dist_temp, dist_init
2881 double precision ggg(3)
2882 double precision sscale,sscagrad
2883 double precision boxshift
2886 cd print '(a)','Enter ESCP'
2888 c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
2889 c & ' iatscp_e=',iatscp_e
2890 if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
2891 do i=iatscp_s,iatscp_e
2892 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2894 xi=0.5D0*(c(1,i)+c(1,i+1))
2895 yi=0.5D0*(c(2,i)+c(2,i+1))
2896 zi=0.5D0*(c(3,i)+c(3,i+1))
2897 call to_box(xi,yi,zi)
2900 c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
2901 c & " nscp_gr",nscp_gr(i)
2902 do iint=1,nscp_gr(i)
2904 do j=iscpstart(i,iint),iscpend(i,iint)
2905 itypj=iabs(itype(j))
2907 c & write (iout,*) "j",j," itypj",itypj
2908 if (itypj.eq.ntyp1) cycle
2909 C Uncomment following three lines for SC-p interactions
2913 C Uncomment following three lines for Ca-p interactions
2921 call to_box(xj,yj,zj)
2922 xj=boxshift(xj-xi,boxxsize)
2923 yj=boxshift(yj-yi,boxysize)
2924 zj=boxshift(zj-zi,boxzsize)
2925 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2926 c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2927 c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2928 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
2929 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),
2931 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2932 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2933 c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),
2934 c & " subchap",subchap," sss",sss
2935 if (sss.gt.0.0d0) then
2938 e1=fac*fac*aad(itypj,iteli)
2939 e2=fac*bad(itypj,iteli)
2940 if (iabs(j-i) .le. 2) then
2943 evdw2_14=evdw2_14+(e1+e2)*sss
2946 evdw2=evdw2+evdwij*sss
2947 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2948 & 'evdw2',i,j,sss,evdwij
2950 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2952 fac=-(evdwij+e1)*rrij*sss
2953 fac=fac+evdwij*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)/expon
2957 C Uncomment following three lines for SC-p interactions
2959 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2961 C Uncomment following line for SC-p interactions
2962 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2964 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2965 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2974 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2975 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2976 gradx_scp(j,i)=expon*gradx_scp(j,i)
2979 C******************************************************************************
2983 C To save time the factor EXPON has been extracted from ALL components
2984 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2987 C******************************************************************************