1 C-----------------------------------------------------------------------
2 double precision function sscale(r)
3 double precision r,gamm
4 include "COMMON.SPLITELE"
5 if(r.lt.r_cut-rlamb) then
7 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
8 gamm=(r-(r_cut-rlamb))/rlamb
9 sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
15 C-----------------------------------------------------------------------
16 subroutine elj_long(evdw)
18 C This subroutine calculates the interaction energy of nonbonded side chains
19 C assuming the LJ potential of interaction.
21 implicit real*8 (a-h,o-z)
23 parameter (accur=1.0d-10)
26 include 'COMMON.LOCAL'
27 include 'COMMON.CHAIN'
28 include 'COMMON.DERIV'
29 include 'COMMON.INTERACT'
30 include 'COMMON.TORSION'
31 include 'COMMON.SBRIDGE'
32 include 'COMMON.NAMES'
33 include 'COMMON.IOUNITS'
34 include 'COMMON.CONTACTS'
36 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
45 C Calculate SC interaction energy.
48 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
49 cd & 'iend=',iend(i,iint)
50 do j=istart(i,iint),iend(i,iint)
56 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
57 if (sss.lt.1.0d0) then
59 eps0ij=eps(itypi,itypj)
61 e1=fac*fac*aa(itypi,itypj)
62 e2=fac*bb(itypi,itypj)
64 evdw=evdw+(1.0d0-sss)*evdwij
66 C Calculate the components of the gradient in DC and X
68 fac=-rrij*(e1+evdwij)*(1.0d0-sss)
73 gvdwx(k,i)=gvdwx(k,i)-gg(k)
74 gvdwx(k,j)=gvdwx(k,j)+gg(k)
75 gvdwc(k,i)=gvdwc(k,i)-gg(k)
76 gvdwc(k,j)=gvdwc(k,j)+gg(k)
84 gvdwc(j,i)=expon*gvdwc(j,i)
85 gvdwx(j,i)=expon*gvdwx(j,i)
88 C******************************************************************************
92 C To save time, the factor of EXPON has been extracted from ALL components
93 C of GVDWC and GRADX. Remember to multiply them by this factor before further
96 C******************************************************************************
99 C-----------------------------------------------------------------------
100 subroutine elj_short(evdw)
102 C This subroutine calculates the interaction energy of nonbonded side chains
103 C assuming the LJ potential of interaction.
105 implicit real*8 (a-h,o-z)
107 parameter (accur=1.0d-10)
110 include 'COMMON.LOCAL'
111 include 'COMMON.CHAIN'
112 include 'COMMON.DERIV'
113 include 'COMMON.INTERACT'
114 include 'COMMON.TORSION'
115 include 'COMMON.SBRIDGE'
116 include 'COMMON.NAMES'
117 include 'COMMON.IOUNITS'
118 include 'COMMON.CONTACTS'
120 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
131 C Calculate SC interaction energy.
134 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
135 cd & 'iend=',iend(i,iint)
136 do j=istart(i,iint),iend(i,iint)
141 C Change 12/1/95 to calculate four-body interactions
142 rij=xj*xj+yj*yj+zj*zj
143 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
144 if (sss.gt.0.0d0) then
146 eps0ij=eps(itypi,itypj)
148 e1=fac*fac*aa(itypi,itypj)
149 e2=fac*bb(itypi,itypj)
153 C Calculate the components of the gradient in DC and X
155 fac=-rrij*(e1+evdwij)*sss
160 gvdwx(k,i)=gvdwx(k,i)-gg(k)
161 gvdwx(k,j)=gvdwx(k,j)+gg(k)
162 gvdwc(k,i)=gvdwc(k,i)-gg(k)
163 gvdwc(k,j)=gvdwc(k,j)+gg(k)
171 gvdwc(j,i)=expon*gvdwc(j,i)
172 gvdwx(j,i)=expon*gvdwx(j,i)
175 C******************************************************************************
179 C To save time, the factor of EXPON has been extracted from ALL components
180 C of GVDWC and GRADX. Remember to multiply them by this factor before further
183 C******************************************************************************
186 C-----------------------------------------------------------------------------
187 subroutine eljk_long(evdw)
189 C This subroutine calculates the interaction energy of nonbonded side chains
190 C assuming the LJK potential of interaction.
192 implicit real*8 (a-h,o-z)
196 include 'COMMON.LOCAL'
197 include 'COMMON.CHAIN'
198 include 'COMMON.DERIV'
199 include 'COMMON.INTERACT'
200 include 'COMMON.IOUNITS'
201 include 'COMMON.NAMES'
204 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
213 C Calculate SC interaction energy.
216 do j=istart(i,iint),iend(i,iint)
221 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
223 e_augm=augm(itypi,itypj)*fac_augm
226 sss=sscale(rij/sigma(itypi,itypj))
227 if (sss.lt.1.0d0) then
228 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
229 fac=r_shift_inv**expon
230 e1=fac*fac*aa(itypi,itypj)
231 e2=fac*bb(itypi,itypj)
233 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
234 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
235 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
236 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
237 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
238 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
239 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
240 evdw=evdw+(1.0d0-sss)*evdwij
242 C Calculate the components of the gradient in DC and X
244 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
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)
267 C-----------------------------------------------------------------------------
268 subroutine eljk_short(evdw)
270 C This subroutine calculates the interaction energy of nonbonded side chains
271 C assuming the LJK potential of interaction.
273 implicit real*8 (a-h,o-z)
277 include 'COMMON.LOCAL'
278 include 'COMMON.CHAIN'
279 include 'COMMON.DERIV'
280 include 'COMMON.INTERACT'
281 include 'COMMON.IOUNITS'
282 include 'COMMON.NAMES'
285 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
294 C Calculate SC interaction energy.
297 do j=istart(i,iint),iend(i,iint)
302 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
304 e_augm=augm(itypi,itypj)*fac_augm
307 sss=sscale(rij/sigma(itypi,itypj))
308 if (sss.gt.0.0d0) then
309 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
310 fac=r_shift_inv**expon
311 e1=fac*fac*aa(itypi,itypj)
312 e2=fac*bb(itypi,itypj)
314 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
315 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
316 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
317 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
318 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
319 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
320 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
323 C Calculate the components of the gradient in DC and X
325 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
331 gvdwx(k,i)=gvdwx(k,i)-gg(k)
332 gvdwx(k,j)=gvdwx(k,j)+gg(k)
333 gvdwc(k,i)=gvdwc(k,i)-gg(k)
334 gvdwc(k,j)=gvdwc(k,j)+gg(k)
342 gvdwc(j,i)=expon*gvdwc(j,i)
343 gvdwx(j,i)=expon*gvdwx(j,i)
348 C-----------------------------------------------------------------------------
349 subroutine ebp_long(evdw)
351 C This subroutine calculates the interaction energy of nonbonded side chains
352 C assuming the Berne-Pechukas potential of interaction.
354 implicit real*8 (a-h,o-z)
358 include 'COMMON.LOCAL'
359 include 'COMMON.CHAIN'
360 include 'COMMON.DERIV'
361 include 'COMMON.NAMES'
362 include 'COMMON.INTERACT'
363 include 'COMMON.IOUNITS'
364 include 'COMMON.CALC'
366 c double precision rrsave(maxdim)
369 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
371 c if (icall.eq.0) then
383 dxi=dc_norm(1,nres+i)
384 dyi=dc_norm(2,nres+i)
385 dzi=dc_norm(3,nres+i)
386 c dsci_inv=dsc_inv(itypi)
387 dsci_inv=vbld_inv(i+nres)
389 C Calculate SC interaction energy.
392 do j=istart(i,iint),iend(i,iint)
395 c dscj_inv=dsc_inv(itypj)
396 dscj_inv=vbld_inv(j+nres)
397 chi1=chi(itypi,itypj)
398 chi2=chi(itypj,itypi)
405 alf12=0.5D0*(alf1+alf2)
409 dxj=dc_norm(1,nres+j)
410 dyj=dc_norm(2,nres+j)
411 dzj=dc_norm(3,nres+j)
412 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
414 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
416 if (sss.lt.1.0d0) then
418 C Calculate the angle-dependent terms of energy & contributions to derivatives.
420 C Calculate whole angle-dependent part of epsilon and contributions
422 fac=(rrij*sigsq)**expon2
423 e1=fac*fac*aa(itypi,itypj)
424 e2=fac*bb(itypi,itypj)
425 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
426 eps2der=evdwij*eps3rt
427 eps3der=evdwij*eps2rt
428 evdwij=evdwij*eps2rt*eps3rt
429 evdw=evdw+evdwij*(1.0d0-sss)
431 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
432 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
433 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
434 cd & restyp(itypi),i,restyp(itypj),j,
435 cd & epsi,sigm,chi1,chi2,chip1,chip2,
436 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
437 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
440 C Calculate gradient components.
441 e1=e1*eps1*eps2rt**2*eps3rt**2
442 fac=-expon*(e1+evdwij)
445 C Calculate radial part of the gradient
449 C Calculate the angular part of the gradient and sum add the contributions
450 C to the appropriate components of the Cartesian gradient.
451 call sc_grad_scale(1.0d0-sss)
459 C-----------------------------------------------------------------------------
460 subroutine ebp_short(evdw)
462 C This subroutine calculates the interaction energy of nonbonded side chains
463 C assuming the Berne-Pechukas potential of interaction.
465 implicit real*8 (a-h,o-z)
469 include 'COMMON.LOCAL'
470 include 'COMMON.CHAIN'
471 include 'COMMON.DERIV'
472 include 'COMMON.NAMES'
473 include 'COMMON.INTERACT'
474 include 'COMMON.IOUNITS'
475 include 'COMMON.CALC'
477 c double precision rrsave(maxdim)
480 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
482 c if (icall.eq.0) then
494 dxi=dc_norm(1,nres+i)
495 dyi=dc_norm(2,nres+i)
496 dzi=dc_norm(3,nres+i)
497 c dsci_inv=dsc_inv(itypi)
498 dsci_inv=vbld_inv(i+nres)
500 C Calculate SC interaction energy.
503 do j=istart(i,iint),iend(i,iint)
506 c dscj_inv=dsc_inv(itypj)
507 dscj_inv=vbld_inv(j+nres)
508 chi1=chi(itypi,itypj)
509 chi2=chi(itypj,itypi)
516 alf12=0.5D0*(alf1+alf2)
520 dxj=dc_norm(1,nres+j)
521 dyj=dc_norm(2,nres+j)
522 dzj=dc_norm(3,nres+j)
523 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
525 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
527 if (sss.gt.0.0d0) then
529 C Calculate the angle-dependent terms of energy & contributions to derivatives.
531 C Calculate whole angle-dependent part of epsilon and contributions
533 fac=(rrij*sigsq)**expon2
534 e1=fac*fac*aa(itypi,itypj)
535 e2=fac*bb(itypi,itypj)
536 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
537 eps2der=evdwij*eps3rt
538 eps3der=evdwij*eps2rt
539 evdwij=evdwij*eps2rt*eps3rt
542 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
543 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
544 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
545 cd & restyp(itypi),i,restyp(itypj),j,
546 cd & epsi,sigm,chi1,chi2,chip1,chip2,
547 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
548 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
551 C Calculate gradient components.
552 e1=e1*eps1*eps2rt**2*eps3rt**2
553 fac=-expon*(e1+evdwij)
556 C Calculate radial part of the gradient
560 C Calculate the angular part of the gradient and sum add the contributions
561 C to the appropriate components of the Cartesian gradient.
562 call sc_grad_scale(sss)
570 C-----------------------------------------------------------------------------
571 subroutine egb_long(evdw)
573 C This subroutine calculates the interaction energy of nonbonded side chains
574 C assuming the Gay-Berne potential of interaction.
576 implicit real*8 (a-h,o-z)
580 include 'COMMON.LOCAL'
581 include 'COMMON.CHAIN'
582 include 'COMMON.DERIV'
583 include 'COMMON.NAMES'
584 include 'COMMON.INTERACT'
585 include 'COMMON.IOUNITS'
586 include 'COMMON.CALC'
587 include 'COMMON.CONTROL'
590 ccccc energy_dec=.false.
591 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
594 c if (icall.eq.0) lprn=.false.
602 dxi=dc_norm(1,nres+i)
603 dyi=dc_norm(2,nres+i)
604 dzi=dc_norm(3,nres+i)
605 c dsci_inv=dsc_inv(itypi)
606 dsci_inv=vbld_inv(i+nres)
607 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
608 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
610 C Calculate SC interaction energy.
613 do j=istart(i,iint),iend(i,iint)
616 c dscj_inv=dsc_inv(itypj)
617 dscj_inv=vbld_inv(j+nres)
618 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
619 c & 1.0d0/vbld(j+nres)
620 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
621 sig0ij=sigma(itypi,itypj)
622 chi1=chi(itypi,itypj)
623 chi2=chi(itypj,itypi)
630 alf12=0.5D0*(alf1+alf2)
634 dxj=dc_norm(1,nres+j)
635 dyj=dc_norm(2,nres+j)
636 dzj=dc_norm(3,nres+j)
637 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
639 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
641 if (sss.lt.1.0d0) then
643 C Calculate angle-dependent terms of energy and contributions to their
647 sig=sig0ij*dsqrt(sigsq)
648 rij_shift=1.0D0/rij-sig+sig0ij
649 c for diagnostics; uncomment
650 c rij_shift=1.2*sig0ij
651 C I hate to put IF's in the loops, but here don't have another choice!!!!
652 if (rij_shift.le.0.0D0) then
654 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
655 cd & restyp(itypi),i,restyp(itypj),j,
656 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
660 c---------------------------------------------------------------
661 rij_shift=1.0D0/rij_shift
663 e1=fac*fac*aa(itypi,itypj)
664 e2=fac*bb(itypi,itypj)
665 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
666 eps2der=evdwij*eps3rt
667 eps3der=evdwij*eps2rt
668 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
669 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
670 evdwij=evdwij*eps2rt*eps3rt
671 evdw=evdw+evdwij*(1.0d0-sss)
673 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
674 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
675 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
676 & restyp(itypi),i,restyp(itypj),j,
677 & epsi,sigm,chi1,chi2,chip1,chip2,
678 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
679 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
683 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
686 C Calculate gradient components.
687 e1=e1*eps1*eps2rt**2*eps3rt**2
688 fac=-expon*(e1+evdwij)*rij_shift
692 C Calculate the radial part of the gradient
696 C Calculate angular part of the gradient.
697 call sc_grad_scale(1.0d0-sss)
702 c write (iout,*) "Number of loop steps in EGB:",ind
703 cccc energy_dec=.false.
706 C-----------------------------------------------------------------------------
707 subroutine egb_short(evdw)
709 C This subroutine calculates the interaction energy of nonbonded side chains
710 C assuming the Gay-Berne potential of interaction.
712 implicit real*8 (a-h,o-z)
716 include 'COMMON.LOCAL'
717 include 'COMMON.CHAIN'
718 include 'COMMON.DERIV'
719 include 'COMMON.NAMES'
720 include 'COMMON.INTERACT'
721 include 'COMMON.IOUNITS'
722 include 'COMMON.CALC'
723 include 'COMMON.CONTROL'
726 ccccc energy_dec=.false.
727 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
730 c if (icall.eq.0) lprn=.false.
738 dxi=dc_norm(1,nres+i)
739 dyi=dc_norm(2,nres+i)
740 dzi=dc_norm(3,nres+i)
741 c dsci_inv=dsc_inv(itypi)
742 dsci_inv=vbld_inv(i+nres)
743 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
744 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
746 C Calculate SC interaction energy.
749 do j=istart(i,iint),iend(i,iint)
752 c dscj_inv=dsc_inv(itypj)
753 dscj_inv=vbld_inv(j+nres)
754 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
755 c & 1.0d0/vbld(j+nres)
756 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
757 sig0ij=sigma(itypi,itypj)
758 chi1=chi(itypi,itypj)
759 chi2=chi(itypj,itypi)
766 alf12=0.5D0*(alf1+alf2)
770 dxj=dc_norm(1,nres+j)
771 dyj=dc_norm(2,nres+j)
772 dzj=dc_norm(3,nres+j)
773 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
775 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
777 if (sss.gt.0.0d0) then
779 C Calculate angle-dependent terms of energy and contributions to their
783 sig=sig0ij*dsqrt(sigsq)
784 rij_shift=1.0D0/rij-sig+sig0ij
785 c for diagnostics; uncomment
786 c rij_shift=1.2*sig0ij
787 C I hate to put IF's in the loops, but here don't have another choice!!!!
788 if (rij_shift.le.0.0D0) then
790 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
791 cd & restyp(itypi),i,restyp(itypj),j,
792 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
796 c---------------------------------------------------------------
797 rij_shift=1.0D0/rij_shift
799 e1=fac*fac*aa(itypi,itypj)
800 e2=fac*bb(itypi,itypj)
801 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
802 eps2der=evdwij*eps3rt
803 eps3der=evdwij*eps2rt
804 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
805 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
806 evdwij=evdwij*eps2rt*eps3rt
809 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
810 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
811 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
812 & restyp(itypi),i,restyp(itypj),j,
813 & epsi,sigm,chi1,chi2,chip1,chip2,
814 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
815 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
819 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
822 C Calculate gradient components.
823 e1=e1*eps1*eps2rt**2*eps3rt**2
824 fac=-expon*(e1+evdwij)*rij_shift
828 C Calculate the radial part of the gradient
832 C Calculate angular part of the gradient.
833 call sc_grad_scale(sss)
838 c write (iout,*) "Number of loop steps in EGB:",ind
839 cccc energy_dec=.false.
842 C-----------------------------------------------------------------------------
843 subroutine egbv_long(evdw)
845 C This subroutine calculates the interaction energy of nonbonded side chains
846 C assuming the Gay-Berne-Vorobjev potential of interaction.
848 implicit real*8 (a-h,o-z)
852 include 'COMMON.LOCAL'
853 include 'COMMON.CHAIN'
854 include 'COMMON.DERIV'
855 include 'COMMON.NAMES'
856 include 'COMMON.INTERACT'
857 include 'COMMON.IOUNITS'
858 include 'COMMON.CALC'
862 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
865 c if (icall.eq.0) lprn=.true.
873 dxi=dc_norm(1,nres+i)
874 dyi=dc_norm(2,nres+i)
875 dzi=dc_norm(3,nres+i)
876 c dsci_inv=dsc_inv(itypi)
877 dsci_inv=vbld_inv(i+nres)
879 C Calculate SC interaction energy.
882 do j=istart(i,iint),iend(i,iint)
885 c dscj_inv=dsc_inv(itypj)
886 dscj_inv=vbld_inv(j+nres)
887 sig0ij=sigma(itypi,itypj)
889 chi1=chi(itypi,itypj)
890 chi2=chi(itypj,itypi)
897 alf12=0.5D0*(alf1+alf2)
901 dxj=dc_norm(1,nres+j)
902 dyj=dc_norm(2,nres+j)
903 dzj=dc_norm(3,nres+j)
904 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
907 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
909 if (sss.lt.1.0d0) then
911 C Calculate angle-dependent terms of energy and contributions to their
915 sig=sig0ij*dsqrt(sigsq)
916 rij_shift=1.0D0/rij-sig+r0ij
917 C I hate to put IF's in the loops, but here don't have another choice!!!!
918 if (rij_shift.le.0.0D0) then
923 c---------------------------------------------------------------
924 rij_shift=1.0D0/rij_shift
926 e1=fac*fac*aa(itypi,itypj)
927 e2=fac*bb(itypi,itypj)
928 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
929 eps2der=evdwij*eps3rt
930 eps3der=evdwij*eps2rt
932 e_augm=augm(itypi,itypj)*fac_augm
933 evdwij=evdwij*eps2rt*eps3rt
934 evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
936 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
937 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
938 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
939 & restyp(itypi),i,restyp(itypj),j,
940 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
941 & chi1,chi2,chip1,chip2,
942 & eps1,eps2rt**2,eps3rt**2,
943 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
946 C Calculate gradient components.
947 e1=e1*eps1*eps2rt**2*eps3rt**2
948 fac=-expon*(e1+evdwij)*rij_shift
950 fac=rij*fac-2*expon*rrij*e_augm
951 C Calculate the radial part of the gradient
955 C Calculate angular part of the gradient.
956 call sc_grad_scale(1.0d0-sss)
962 C-----------------------------------------------------------------------------
963 subroutine egbv_short(evdw)
965 C This subroutine calculates the interaction energy of nonbonded side chains
966 C assuming the Gay-Berne-Vorobjev potential of interaction.
968 implicit real*8 (a-h,o-z)
972 include 'COMMON.LOCAL'
973 include 'COMMON.CHAIN'
974 include 'COMMON.DERIV'
975 include 'COMMON.NAMES'
976 include 'COMMON.INTERACT'
977 include 'COMMON.IOUNITS'
978 include 'COMMON.CALC'
982 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
985 c if (icall.eq.0) lprn=.true.
993 dxi=dc_norm(1,nres+i)
994 dyi=dc_norm(2,nres+i)
995 dzi=dc_norm(3,nres+i)
996 c dsci_inv=dsc_inv(itypi)
997 dsci_inv=vbld_inv(i+nres)
999 C Calculate SC interaction energy.
1001 do iint=1,nint_gr(i)
1002 do j=istart(i,iint),iend(i,iint)
1005 c dscj_inv=dsc_inv(itypj)
1006 dscj_inv=vbld_inv(j+nres)
1007 sig0ij=sigma(itypi,itypj)
1008 r0ij=r0(itypi,itypj)
1009 chi1=chi(itypi,itypj)
1010 chi2=chi(itypj,itypi)
1017 alf12=0.5D0*(alf1+alf2)
1021 dxj=dc_norm(1,nres+j)
1022 dyj=dc_norm(2,nres+j)
1023 dzj=dc_norm(3,nres+j)
1024 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1027 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1029 if (sss.gt.0.0d0) then
1031 C Calculate angle-dependent terms of energy and contributions to their
1035 sig=sig0ij*dsqrt(sigsq)
1036 rij_shift=1.0D0/rij-sig+r0ij
1037 C I hate to put IF's in the loops, but here don't have another choice!!!!
1038 if (rij_shift.le.0.0D0) then
1043 c---------------------------------------------------------------
1044 rij_shift=1.0D0/rij_shift
1045 fac=rij_shift**expon
1046 e1=fac*fac*aa(itypi,itypj)
1047 e2=fac*bb(itypi,itypj)
1048 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1049 eps2der=evdwij*eps3rt
1050 eps3der=evdwij*eps2rt
1051 fac_augm=rrij**expon
1052 e_augm=augm(itypi,itypj)*fac_augm
1053 evdwij=evdwij*eps2rt*eps3rt
1054 evdw=evdw+(evdwij+e_augm)*sss
1056 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1057 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1058 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1059 & restyp(itypi),i,restyp(itypj),j,
1060 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1061 & chi1,chi2,chip1,chip2,
1062 & eps1,eps2rt**2,eps3rt**2,
1063 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1066 C Calculate gradient components.
1067 e1=e1*eps1*eps2rt**2*eps3rt**2
1068 fac=-expon*(e1+evdwij)*rij_shift
1070 fac=rij*fac-2*expon*rrij*e_augm
1071 C Calculate the radial part of the gradient
1075 C Calculate angular part of the gradient.
1076 call sc_grad_scale(sss)
1082 C----------------------------------------------------------------------------
1083 subroutine sc_grad_scale(scalfac)
1084 implicit real*8 (a-h,o-z)
1085 include 'DIMENSIONS'
1086 include 'COMMON.CHAIN'
1087 include 'COMMON.DERIV'
1088 include 'COMMON.CALC'
1089 include 'COMMON.IOUNITS'
1090 double precision dcosom1(3),dcosom2(3)
1091 double precision scalfac
1092 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1093 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1094 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1095 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1099 c eom12=evdwij*eps1_om12
1101 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1102 c & " sigder",sigder
1103 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1104 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1106 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1107 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1110 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1112 c write (iout,*) "gg",(gg(k),k=1,3)
1114 gvdwx(k,i)=gvdwx(k,i)-gg(k)
1115 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1116 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1117 gvdwx(k,j)=gvdwx(k,j)+gg(k)
1118 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1119 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1120 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1121 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1122 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1123 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1126 C Calculate the components of the gradient in DC and X
1129 gvdwc(l,i)=gvdwc(l,i)-gg(l)
1130 gvdwc(l,j)=gvdwc(l,j)+gg(l)
1134 C--------------------------------------------------------------------------
1135 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1137 C This subroutine calculates the average interaction energy and its gradient
1138 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1139 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1140 C The potential depends both on the distance of peptide-group centers and on
1141 C the orientation of the CA-CA virtual bonds.
1143 implicit real*8 (a-h,o-z)
1147 include 'DIMENSIONS'
1148 include 'COMMON.CONTROL'
1149 include 'COMMON.SETUP'
1150 include 'COMMON.IOUNITS'
1151 include 'COMMON.GEO'
1152 include 'COMMON.VAR'
1153 include 'COMMON.LOCAL'
1154 include 'COMMON.CHAIN'
1155 include 'COMMON.DERIV'
1156 include 'COMMON.INTERACT'
1157 include 'COMMON.CONTACTS'
1158 include 'COMMON.TORSION'
1159 include 'COMMON.VECTORS'
1160 include 'COMMON.FFIELD'
1161 include 'COMMON.TIME1'
1162 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1163 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1164 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1165 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1166 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1167 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1169 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1171 double precision scal_el /1.0d0/
1173 double precision scal_el /0.5d0/
1176 C 13-go grudnia roku pamietnego...
1177 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1178 & 0.0d0,1.0d0,0.0d0,
1179 & 0.0d0,0.0d0,1.0d0/
1180 cd write(iout,*) 'In EELEC'
1182 cd write(iout,*) 'Type',i
1183 cd write(iout,*) 'B1',B1(:,i)
1184 cd write(iout,*) 'B2',B2(:,i)
1185 cd write(iout,*) 'CC',CC(:,:,i)
1186 cd write(iout,*) 'DD',DD(:,:,i)
1187 cd write(iout,*) 'EE',EE(:,:,i)
1189 cd call check_vecgrad
1191 if (icheckgrad.eq.1) then
1193 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1195 dc_norm(k,i)=dc(k,i)*fac
1197 c write (iout,*) 'i',i,' fac',fac
1200 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1201 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1202 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1203 c call vec_and_deriv
1209 time_mat=time_mat+MPI_Wtime()-time01
1213 cd write (iout,*) 'i=',i
1215 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1218 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1219 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1232 cd print '(a)','Enter EELEC'
1233 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1235 gel_loc_loc(i)=0.0d0
1240 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1242 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1244 do i=iturn3_start,iturn3_end
1248 dx_normi=dc_norm(1,i)
1249 dy_normi=dc_norm(2,i)
1250 dz_normi=dc_norm(3,i)
1251 xmedi=c(1,i)+0.5d0*dxi
1252 ymedi=c(2,i)+0.5d0*dyi
1253 zmedi=c(3,i)+0.5d0*dzi
1255 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1256 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1257 num_cont_hb(i)=num_conti
1259 do i=iturn4_start,iturn4_end
1263 dx_normi=dc_norm(1,i)
1264 dy_normi=dc_norm(2,i)
1265 dz_normi=dc_norm(3,i)
1266 xmedi=c(1,i)+0.5d0*dxi
1267 ymedi=c(2,i)+0.5d0*dyi
1268 zmedi=c(3,i)+0.5d0*dzi
1269 num_conti=num_cont_hb(i)
1270 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1271 if (wturn4.gt.0.0d0) call eturn4(i,eello_turn4)
1272 num_cont_hb(i)=num_conti
1275 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1277 do i=iatel_s,iatel_e
1281 dx_normi=dc_norm(1,i)
1282 dy_normi=dc_norm(2,i)
1283 dz_normi=dc_norm(3,i)
1284 xmedi=c(1,i)+0.5d0*dxi
1285 ymedi=c(2,i)+0.5d0*dyi
1286 zmedi=c(3,i)+0.5d0*dzi
1287 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1288 num_conti=num_cont_hb(i)
1289 do j=ielstart(i),ielend(i)
1290 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1292 num_cont_hb(i)=num_conti
1294 c write (iout,*) "Number of loop steps in EELEC:",ind
1296 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1297 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1299 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1300 ccc eel_loc=eel_loc+eello_turn3
1301 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1304 C-------------------------------------------------------------------------------
1305 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1306 implicit real*8 (a-h,o-z)
1307 include 'DIMENSIONS'
1311 include 'COMMON.CONTROL'
1312 include 'COMMON.IOUNITS'
1313 include 'COMMON.GEO'
1314 include 'COMMON.VAR'
1315 include 'COMMON.LOCAL'
1316 include 'COMMON.CHAIN'
1317 include 'COMMON.DERIV'
1318 include 'COMMON.INTERACT'
1319 include 'COMMON.CONTACTS'
1320 include 'COMMON.TORSION'
1321 include 'COMMON.VECTORS'
1322 include 'COMMON.FFIELD'
1323 include 'COMMON.TIME1'
1324 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1325 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1326 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1327 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1328 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1329 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1331 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1333 double precision scal_el /1.0d0/
1335 double precision scal_el /0.5d0/
1338 C 13-go grudnia roku pamietnego...
1339 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1340 & 0.0d0,1.0d0,0.0d0,
1341 & 0.0d0,0.0d0,1.0d0/
1342 c time00=MPI_Wtime()
1343 cd write (iout,*) "eelecij",i,j
1347 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1348 aaa=app(iteli,itelj)
1349 bbb=bpp(iteli,itelj)
1350 ael6i=ael6(iteli,itelj)
1351 ael3i=ael3(iteli,itelj)
1355 dx_normj=dc_norm(1,j)
1356 dy_normj=dc_norm(2,j)
1357 dz_normj=dc_norm(3,j)
1358 xj=c(1,j)+0.5D0*dxj-xmedi
1359 yj=c(2,j)+0.5D0*dyj-ymedi
1360 zj=c(3,j)+0.5D0*dzj-zmedi
1361 rij=xj*xj+yj*yj+zj*zj
1365 c For extracting the short-range part of Evdwpp
1366 sss=sscale(rij/rpp(iteli,itelj))
1370 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1371 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1372 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1373 fac=cosa-3.0D0*cosb*cosg
1375 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1376 if (j.eq.i+2) ev1=scal_el*ev1
1381 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1384 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1385 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1387 evdw1=evdw1+evdwij*(1.0d0-sss)
1388 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1389 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1390 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1391 cd & xmedi,ymedi,zmedi,xj,yj,zj
1393 if (energy_dec) then
1394 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1395 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1399 C Calculate contributions to the Cartesian gradient.
1402 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1403 facel=-3*rrmij*(el1+eesij)
1409 * Radial derivatives. First process both termini of the fragment (i,j)
1415 c ghalf=0.5D0*ggg(k)
1416 c gelc(k,i)=gelc(k,i)+ghalf
1417 c gelc(k,j)=gelc(k,j)+ghalf
1419 c 9/28/08 AL Gradient compotents will be summed only at the end
1421 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1422 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1425 * Loop over residues i+1 thru j-1.
1429 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1436 c ghalf=0.5D0*ggg(k)
1437 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1438 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1440 c 9/28/08 AL Gradient compotents will be summed only at the end
1442 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1443 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1446 * Loop over residues i+1 thru j-1.
1450 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1454 facvdw=ev1+evdwij*(1.0d0-sss)
1457 fac=-3*rrmij*(facvdw+facvdw+facel)
1462 * Radial derivatives. First process both termini of the fragment (i,j)
1468 c ghalf=0.5D0*ggg(k)
1469 c gelc(k,i)=gelc(k,i)+ghalf
1470 c gelc(k,j)=gelc(k,j)+ghalf
1472 c 9/28/08 AL Gradient compotents will be summed only at the end
1474 gelc_long(k,j)=gelc(k,j)+ggg(k)
1475 gelc_long(k,i)=gelc(k,i)-ggg(k)
1478 * Loop over residues i+1 thru j-1.
1482 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1485 c 9/28/08 AL Gradient compotents will be summed only at the end
1490 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1491 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1497 ecosa=2.0D0*fac3*fac1+fac4
1500 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1501 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1503 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1504 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1506 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1507 cd & (dcosg(k),k=1,3)
1509 ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1512 c ghalf=0.5D0*ggg(k)
1513 c gelc(k,i)=gelc(k,i)+ghalf
1514 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1515 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1516 c gelc(k,j)=gelc(k,j)+ghalf
1517 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1518 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1522 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1527 & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1528 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1530 & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1531 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1532 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1533 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1535 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1536 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1537 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1539 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1540 C energy of a peptide unit is assumed in the form of a second-order
1541 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1542 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1543 C are computed for EVERY pair of non-contiguous peptide groups.
1545 if (j.lt.nres-1) then
1556 muij(kkk)=mu(k,i)*mu(l,j)
1559 cd write (iout,*) 'EELEC: i',i,' j',j
1560 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1561 cd write(iout,*) 'muij',muij
1562 ury=scalar(uy(1,i),erij)
1563 urz=scalar(uz(1,i),erij)
1564 vry=scalar(uy(1,j),erij)
1565 vrz=scalar(uz(1,j),erij)
1566 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1567 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1568 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1569 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1570 fac=dsqrt(-ael6i)*r3ij
1575 cd write (iout,'(4i5,4f10.5)')
1576 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1577 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1578 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1579 cd & uy(:,j),uz(:,j)
1580 cd write (iout,'(4f10.5)')
1581 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1582 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1583 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
1584 cd write (iout,'(9f10.5/)')
1585 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1586 C Derivatives of the elements of A in virtual-bond vectors
1587 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1589 uryg(k,1)=scalar(erder(1,k),uy(1,i))
1590 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1591 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1592 urzg(k,1)=scalar(erder(1,k),uz(1,i))
1593 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1594 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1595 vryg(k,1)=scalar(erder(1,k),uy(1,j))
1596 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1597 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1598 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1599 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1600 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1602 C Compute radial contributions to the gradient
1620 C Add the contributions coming from er
1623 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1624 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1625 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1626 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1629 C Derivatives in DC(i)
1630 cgrad ghalf1=0.5d0*agg(k,1)
1631 cgrad ghalf2=0.5d0*agg(k,2)
1632 cgrad ghalf3=0.5d0*agg(k,3)
1633 cgrad ghalf4=0.5d0*agg(k,4)
1634 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1635 & -3.0d0*uryg(k,2)*vry)!+ghalf1
1636 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1637 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
1638 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1639 & -3.0d0*urzg(k,2)*vry)!+ghalf3
1640 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1641 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
1642 C Derivatives in DC(i+1)
1643 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
1644 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
1645 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
1646 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
1647 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
1648 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
1649 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
1650 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
1651 C Derivatives in DC(j)
1652 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
1653 & -3.0d0*vryg(k,2)*ury)!+ghalf1
1654 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
1655 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
1656 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
1657 & -3.0d0*vryg(k,2)*urz)!+ghalf3
1658 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
1659 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
1660 C Derivatives in DC(j+1) or DC(nres-1)
1661 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
1662 & -3.0d0*vryg(k,3)*ury)
1663 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
1664 & -3.0d0*vrzg(k,3)*ury)
1665 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
1666 & -3.0d0*vryg(k,3)*urz)
1667 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
1668 & -3.0d0*vrzg(k,3)*urz)
1669 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
1671 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
1684 aggi(k,l)=-aggi(k,l)
1685 aggi1(k,l)=-aggi1(k,l)
1686 aggj(k,l)=-aggj(k,l)
1687 aggj1(k,l)=-aggj1(k,l)
1690 if (j.lt.nres-1) then
1696 aggi(k,l)=-aggi(k,l)
1697 aggi1(k,l)=-aggi1(k,l)
1698 aggj(k,l)=-aggj(k,l)
1699 aggj1(k,l)=-aggj1(k,l)
1710 aggi(k,l)=-aggi(k,l)
1711 aggi1(k,l)=-aggi1(k,l)
1712 aggj(k,l)=-aggj(k,l)
1713 aggj1(k,l)=-aggj1(k,l)
1718 IF (wel_loc.gt.0.0d0) THEN
1719 C Contribution to the local-electrostatic energy coming from the i-j pair
1720 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
1722 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
1724 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1725 & 'eelloc',i,j,eel_loc_ij
1727 eel_loc=eel_loc+eel_loc_ij
1728 C Partial derivatives in virtual-bond dihedral angles gamma
1730 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
1731 & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1732 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1733 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
1734 & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1735 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1736 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
1738 ggg(l)=agg(l,1)*muij(1)+
1739 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1740 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
1741 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
1742 cgrad ghalf=0.5d0*ggg(l)
1743 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
1744 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
1748 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
1751 C Remaining derivatives of eello
1753 gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1754 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1755 gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1756 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1757 gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1758 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1759 gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1760 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1763 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
1764 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
1765 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
1766 & .and. num_conti.le.maxconts) then
1767 c write (iout,*) i,j," entered corr"
1769 C Calculate the contact function. The ith column of the array JCONT will
1770 C contain the numbers of atoms that make contacts with the atom I (of numbers
1771 C greater than I). The arrays FACONT and GACONT will contain the values of
1772 C the contact function and its derivative.
1773 c r0ij=1.02D0*rpp(iteli,itelj)
1774 c r0ij=1.11D0*rpp(iteli,itelj)
1775 r0ij=2.20D0*rpp(iteli,itelj)
1776 c r0ij=1.55D0*rpp(iteli,itelj)
1777 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
1778 if (fcont.gt.0.0D0) then
1779 num_conti=num_conti+1
1780 if (num_conti.gt.maxconts) then
1781 write (iout,*) 'WARNING - max. # of contacts exceeded;',
1782 & ' will skip next contacts for this conf.'
1784 jcont_hb(num_conti,i)=j
1785 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
1786 cd & " jcont_hb",jcont_hb(num_conti,i)
1787 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
1788 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1789 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
1791 d_cont(num_conti,i)=rij
1792 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
1793 C --- Electrostatic-interaction matrix ---
1794 a_chuj(1,1,num_conti,i)=a22
1795 a_chuj(1,2,num_conti,i)=a23
1796 a_chuj(2,1,num_conti,i)=a32
1797 a_chuj(2,2,num_conti,i)=a33
1798 C --- Gradient of rij
1800 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
1807 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
1808 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
1809 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
1810 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
1811 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
1816 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
1817 C Calculate contact energies
1819 wij=cosa-3.0D0*cosb*cosg
1822 c fac3=dsqrt(-ael6i)/r0ij**3
1823 fac3=dsqrt(-ael6i)*r3ij
1824 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
1825 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
1826 if (ees0tmp.gt.0) then
1827 ees0pij=dsqrt(ees0tmp)
1831 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
1832 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
1833 if (ees0tmp.gt.0) then
1834 ees0mij=dsqrt(ees0tmp)
1839 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
1840 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
1841 C Diagnostics. Comment out or remove after debugging!
1842 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
1843 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
1844 c ees0m(num_conti,i)=0.0D0
1846 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
1847 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
1848 C Angular derivatives of the contact function
1849 ees0pij1=fac3/ees0pij
1850 ees0mij1=fac3/ees0mij
1851 fac3p=-3.0D0*fac3*rrmij
1852 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
1853 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
1855 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
1856 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
1857 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
1858 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
1859 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
1860 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
1861 ecosap=ecosa1+ecosa2
1862 ecosbp=ecosb1+ecosb2
1863 ecosgp=ecosg1+ecosg2
1864 ecosam=ecosa1-ecosa2
1865 ecosbm=ecosb1-ecosb2
1866 ecosgm=ecosg1-ecosg2
1875 facont_hb(num_conti,i)=fcont
1876 fprimcont=fprimcont/rij
1877 cd facont_hb(num_conti,i)=1.0D0
1878 C Following line is for diagnostics.
1881 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1882 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1885 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
1886 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
1888 gggp(1)=gggp(1)+ees0pijp*xj
1889 gggp(2)=gggp(2)+ees0pijp*yj
1890 gggp(3)=gggp(3)+ees0pijp*zj
1891 gggm(1)=gggm(1)+ees0mijp*xj
1892 gggm(2)=gggm(2)+ees0mijp*yj
1893 gggm(3)=gggm(3)+ees0mijp*zj
1894 C Derivatives due to the contact function
1895 gacont_hbr(1,num_conti,i)=fprimcont*xj
1896 gacont_hbr(2,num_conti,i)=fprimcont*yj
1897 gacont_hbr(3,num_conti,i)=fprimcont*zj
1900 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
1901 c following the change of gradient-summation algorithm.
1903 cgrad ghalfp=0.5D0*gggp(k)
1904 cgrad ghalfm=0.5D0*gggm(k)
1905 gacontp_hb1(k,num_conti,i)=!ghalfp
1906 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
1907 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1908 gacontp_hb2(k,num_conti,i)=!ghalfp
1909 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
1910 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1911 gacontp_hb3(k,num_conti,i)=gggp(k)
1912 gacontm_hb1(k,num_conti,i)=!ghalfm
1913 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
1914 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1915 gacontm_hb2(k,num_conti,i)=!ghalfm
1916 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
1917 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1918 gacontm_hb3(k,num_conti,i)=gggm(k)
1921 endif ! num_conti.le.maxconts
1924 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
1927 ghalf=0.5d0*agg(l,k)
1928 aggi(l,k)=aggi(l,k)+ghalf
1929 aggi1(l,k)=aggi1(l,k)+agg(l,k)
1930 aggj(l,k)=aggj(l,k)+ghalf
1933 if (j.eq.nres-1 .and. i.lt.j-2) then
1936 aggj1(l,k)=aggj1(l,k)+agg(l,k)
1941 c t_eelecij=t_eelecij+MPI_Wtime()-time00
1944 C-----------------------------------------------------------------------
1945 subroutine evdwpp_short(evdw1)
1949 implicit real*8 (a-h,o-z)
1950 include 'DIMENSIONS'
1951 include 'COMMON.CONTROL'
1952 include 'COMMON.IOUNITS'
1953 include 'COMMON.GEO'
1954 include 'COMMON.VAR'
1955 include 'COMMON.LOCAL'
1956 include 'COMMON.CHAIN'
1957 include 'COMMON.DERIV'
1958 include 'COMMON.INTERACT'
1959 include 'COMMON.CONTACTS'
1960 include 'COMMON.TORSION'
1961 include 'COMMON.VECTORS'
1962 include 'COMMON.FFIELD'
1964 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1966 double precision scal_el /1.0d0/
1968 double precision scal_el /0.5d0/
1971 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
1972 c & " iatel_e_vdw",iatel_e_vdw
1974 do i=iatel_s_vdw,iatel_e_vdw
1978 dx_normi=dc_norm(1,i)
1979 dy_normi=dc_norm(2,i)
1980 dz_normi=dc_norm(3,i)
1981 xmedi=c(1,i)+0.5d0*dxi
1982 ymedi=c(2,i)+0.5d0*dyi
1983 zmedi=c(3,i)+0.5d0*dzi
1985 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
1986 c & ' ielend',ielend_vdw(i)
1988 do j=ielstart_vdw(i),ielend_vdw(i)
1992 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1993 aaa=app(iteli,itelj)
1994 bbb=bpp(iteli,itelj)
1998 dx_normj=dc_norm(1,j)
1999 dy_normj=dc_norm(2,j)
2000 dz_normj=dc_norm(3,j)
2001 xj=c(1,j)+0.5D0*dxj-xmedi
2002 yj=c(2,j)+0.5D0*dyj-ymedi
2003 zj=c(3,j)+0.5D0*dzj-zmedi
2004 rij=xj*xj+yj*yj+zj*zj
2007 sss=sscale(rij/rpp(iteli,itelj))
2008 if (sss.gt.0.0d0) then
2013 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2014 if (j.eq.i+2) ev1=scal_el*ev1
2017 if (energy_dec) then
2018 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2020 evdw1=evdw1+evdwij*sss
2022 C Calculate contributions to the Cartesian gradient.
2024 facvdw=-6*rrmij*(ev1+evdwij)*sss
2029 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2030 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2037 C-----------------------------------------------------------------------------
2038 subroutine escp_long(evdw2,evdw2_14)
2040 C This subroutine calculates the excluded-volume interaction energy between
2041 C peptide-group centers and side chains and its gradient in virtual-bond and
2042 C side-chain vectors.
2044 implicit real*8 (a-h,o-z)
2045 include 'DIMENSIONS'
2046 include 'COMMON.GEO'
2047 include 'COMMON.VAR'
2048 include 'COMMON.LOCAL'
2049 include 'COMMON.CHAIN'
2050 include 'COMMON.DERIV'
2051 include 'COMMON.INTERACT'
2052 include 'COMMON.FFIELD'
2053 include 'COMMON.IOUNITS'
2054 include 'COMMON.CONTROL'
2058 cd print '(a)','Enter ESCP'
2059 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2060 do i=iatscp_s,iatscp_e
2062 xi=0.5D0*(c(1,i)+c(1,i+1))
2063 yi=0.5D0*(c(2,i)+c(2,i+1))
2064 zi=0.5D0*(c(3,i)+c(3,i+1))
2066 do iint=1,nscp_gr(i)
2068 do j=iscpstart(i,iint),iscpend(i,iint)
2070 C Uncomment following three lines for SC-p interactions
2074 C Uncomment following three lines for Ca-p interactions
2078 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2080 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2082 if (sss.lt.1.0d0) then
2085 e1=fac*fac*aad(itypj,iteli)
2086 e2=fac*bad(itypj,iteli)
2087 if (iabs(j-i) .le. 2) then
2090 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2093 evdw2=evdw2+evdwij*(1.0d0-sss)
2094 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2095 & 'evdw2',i,j,evdwij
2097 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2099 fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2103 C Uncomment following three lines for SC-p interactions
2105 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2107 C Uncomment following line for SC-p interactions
2108 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2110 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2111 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2120 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2121 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2122 gradx_scp(j,i)=expon*gradx_scp(j,i)
2125 C******************************************************************************
2129 C To save time the factor EXPON has been extracted from ALL components
2130 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2133 C******************************************************************************
2136 C-----------------------------------------------------------------------------
2137 subroutine escp_short(evdw2,evdw2_14)
2139 C This subroutine calculates the excluded-volume interaction energy between
2140 C peptide-group centers and side chains and its gradient in virtual-bond and
2141 C side-chain vectors.
2143 implicit real*8 (a-h,o-z)
2144 include 'DIMENSIONS'
2145 include 'COMMON.GEO'
2146 include 'COMMON.VAR'
2147 include 'COMMON.LOCAL'
2148 include 'COMMON.CHAIN'
2149 include 'COMMON.DERIV'
2150 include 'COMMON.INTERACT'
2151 include 'COMMON.FFIELD'
2152 include 'COMMON.IOUNITS'
2153 include 'COMMON.CONTROL'
2157 cd print '(a)','Enter ESCP'
2158 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2159 do i=iatscp_s,iatscp_e
2161 xi=0.5D0*(c(1,i)+c(1,i+1))
2162 yi=0.5D0*(c(2,i)+c(2,i+1))
2163 zi=0.5D0*(c(3,i)+c(3,i+1))
2165 do iint=1,nscp_gr(i)
2167 do j=iscpstart(i,iint),iscpend(i,iint)
2169 C Uncomment following three lines for SC-p interactions
2173 C Uncomment following three lines for Ca-p interactions
2177 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2179 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2181 if (sss.gt.0.0d0) then
2184 e1=fac*fac*aad(itypj,iteli)
2185 e2=fac*bad(itypj,iteli)
2186 if (iabs(j-i) .le. 2) then
2189 evdw2_14=evdw2_14+(e1+e2)*sss
2192 evdw2=evdw2+evdwij*sss
2193 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2194 & 'evdw2',i,j,evdwij
2196 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2198 fac=-(evdwij+e1)*rrij*sss
2202 C Uncomment following three lines for SC-p interactions
2204 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2206 C Uncomment following line for SC-p interactions
2207 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2209 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2210 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2219 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2220 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2221 gradx_scp(j,i)=expon*gradx_scp(j,i)
2224 C******************************************************************************
2228 C To save time the factor EXPON has been extracted from ALL components
2229 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2232 C******************************************************************************