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
60 e1=fac*fac*aa(itypi,itypj)
61 e2=fac*bb(itypi,itypj)
63 evdw=evdw+(1.0d0-sss)*evdwij
65 C Calculate the components of the gradient in DC and X
67 fac=-rrij*(e1+evdwij)*(1.0d0-sss)
72 gvdwx(k,i)=gvdwx(k,i)-gg(k)
73 gvdwx(k,j)=gvdwx(k,j)+gg(k)
77 gvdwc(l,k)=gvdwc(l,k)+gg(l)
86 gvdwc(j,i)=expon*gvdwc(j,i)
87 gvdwx(j,i)=expon*gvdwx(j,i)
90 C******************************************************************************
94 C To save time, the factor of EXPON has been extracted from ALL components
95 C of GVDWC and GRADX. Remember to multiply them by this factor before further
98 C******************************************************************************
101 C-----------------------------------------------------------------------
102 subroutine elj_short(evdw)
104 C This subroutine calculates the interaction energy of nonbonded side chains
105 C assuming the LJ potential of interaction.
107 implicit real*8 (a-h,o-z)
109 parameter (accur=1.0d-10)
112 include 'COMMON.LOCAL'
113 include 'COMMON.CHAIN'
114 include 'COMMON.DERIV'
115 include 'COMMON.INTERACT'
116 include 'COMMON.TORSION'
117 include 'COMMON.SBRIDGE'
118 include 'COMMON.NAMES'
119 include 'COMMON.IOUNITS'
120 include 'COMMON.CONTACTS'
122 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 rij=xj*xj+yj*yj+zj*zj
142 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
143 if (sss.gt.0.0d0) then
146 e1=fac*fac*aa(itypi,itypj)
147 e2=fac*bb(itypi,itypj)
151 C Calculate the components of the gradient in DC and X
153 fac=-rrij*(e1+evdwij)*sss
158 gvdwx(k,i)=gvdwx(k,i)-gg(k)
159 gvdwx(k,j)=gvdwx(k,j)+gg(k)
163 gvdwc(l,k)=gvdwc(l,k)+gg(l)
172 gvdwc(j,i)=expon*gvdwc(j,i)
173 gvdwx(j,i)=expon*gvdwx(j,i)
176 C******************************************************************************
180 C To save time, the factor of EXPON has been extracted from ALL components
181 C of GVDWC and GRADX. Remember to multiply them by this factor before further
184 C******************************************************************************
187 C-----------------------------------------------------------------------------
188 subroutine eljk_long(evdw)
190 C This subroutine calculates the interaction energy of nonbonded side chains
191 C assuming the LJK potential of interaction.
193 implicit real*8 (a-h,o-z)
197 include 'COMMON.LOCAL'
198 include 'COMMON.CHAIN'
199 include 'COMMON.DERIV'
200 include 'COMMON.INTERACT'
201 include 'COMMON.IOUNITS'
202 include 'COMMON.NAMES'
205 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
214 C Calculate SC interaction energy.
217 do j=istart(i,iint),iend(i,iint)
222 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
224 e_augm=augm(itypi,itypj)*fac_augm
227 sss=sscale(rij/sigma(itypi,itypj))
229 if (sss.lt.1.0d0) then
231 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
232 fac=r_shift_inv**expon
233 e1=fac*fac*aa(itypi,itypj)
234 e2=fac*bb(itypi,itypj)
236 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
237 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
238 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
239 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
240 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
241 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
242 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
243 evdw=evdw+evdwij*(1.0d0-sss)
245 C Calculate the components of the gradient in DC and X
247 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
253 gvdwx(k,i)=gvdwx(k,i)-gg(k)
254 gvdwx(k,j)=gvdwx(k,j)+gg(k)
258 gvdwc(l,k)=gvdwc(l,k)+gg(l)
269 gvdwc(j,i)=expon*gvdwc(j,i)
270 gvdwx(j,i)=expon*gvdwx(j,i)
275 C-----------------------------------------------------------------------------
276 subroutine eljk_short(evdw)
278 C This subroutine calculates the interaction energy of nonbonded side chains
279 C assuming the LJK potential of interaction.
281 implicit real*8 (a-h,o-z)
285 include 'COMMON.LOCAL'
286 include 'COMMON.CHAIN'
287 include 'COMMON.DERIV'
288 include 'COMMON.INTERACT'
289 include 'COMMON.IOUNITS'
290 include 'COMMON.NAMES'
293 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
302 C Calculate SC interaction energy.
305 do j=istart(i,iint),iend(i,iint)
310 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
312 e_augm=augm(itypi,itypj)*fac_augm
315 sss=sscale(rij/sigma(itypi,itypj))
317 if (sss.gt.0.0d0) then
319 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
320 fac=r_shift_inv**expon
321 e1=fac*fac*aa(itypi,itypj)
322 e2=fac*bb(itypi,itypj)
324 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
325 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
326 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
327 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
328 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
329 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
330 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
333 C Calculate the components of the gradient in DC and X
335 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
341 gvdwx(k,i)=gvdwx(k,i)-gg(k)
342 gvdwx(k,j)=gvdwx(k,j)+gg(k)
346 gvdwc(l,k)=gvdwc(l,k)+gg(l)
357 gvdwc(j,i)=expon*gvdwc(j,i)
358 gvdwx(j,i)=expon*gvdwx(j,i)
363 C-----------------------------------------------------------------------------
364 subroutine ebp_long(evdw)
366 C This subroutine calculates the interaction energy of nonbonded side chains
367 C assuming the Berne-Pechukas potential of interaction.
369 implicit real*8 (a-h,o-z)
373 include 'COMMON.LOCAL'
374 include 'COMMON.CHAIN'
375 include 'COMMON.DERIV'
376 include 'COMMON.NAMES'
377 include 'COMMON.INTERACT'
378 include 'COMMON.IOUNITS'
379 include 'COMMON.CALC'
381 c double precision rrsave(maxdim)
384 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
386 c if (icall.eq.0) then
398 dxi=dc_norm(1,nres+i)
399 dyi=dc_norm(2,nres+i)
400 dzi=dc_norm(3,nres+i)
401 c dsci_inv=dsc_inv(itypi)
402 dsci_inv=vbld_inv(i+nres)
404 C Calculate SC interaction energy.
407 do j=istart(i,iint),iend(i,iint)
410 c dscj_inv=dsc_inv(itypj)
411 dscj_inv=vbld_inv(j+nres)
412 chi1=chi(itypi,itypj)
413 chi2=chi(itypj,itypi)
420 alf12=0.5D0*(alf1+alf2)
421 C For diagnostics only!!!
434 dxj=dc_norm(1,nres+j)
435 dyj=dc_norm(2,nres+j)
436 dzj=dc_norm(3,nres+j)
437 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
438 cd if (icall.eq.0) then
444 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
446 if (sss.lt.1.0d0) then
448 C Calculate the angle-dependent terms of energy & contributions to derivatives.
450 C Calculate whole angle-dependent part of epsilon and contributions
452 fac=(rrij*sigsq)**expon2
453 e1=fac*fac*aa(itypi,itypj)
454 e2=fac*bb(itypi,itypj)
455 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
456 eps2der=evdwij*eps3rt
457 eps3der=evdwij*eps2rt
458 evdwij=evdwij*eps2rt*eps3rt
459 evdw=evdw+evdwij*(1.0d0-sss)
461 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
462 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
463 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
464 cd & restyp(itypi),i,restyp(itypj),j,
465 cd & epsi,sigm,chi1,chi2,chip1,chip2,
466 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
467 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
470 C Calculate gradient components.
471 e1=e1*eps1*eps2rt**2*eps3rt**2
472 fac=-expon*(e1+evdwij)
475 C Calculate radial part of the gradient
479 C Calculate the angular part of the gradient and sum add the contributions
480 C to the appropriate components of the Cartesian gradient.
481 call sc_grad_scale(1.0d0-sss)
491 C-----------------------------------------------------------------------------
492 subroutine ebp_short(evdw)
494 C This subroutine calculates the interaction energy of nonbonded side chains
495 C assuming the Berne-Pechukas potential of interaction.
497 implicit real*8 (a-h,o-z)
501 include 'COMMON.LOCAL'
502 include 'COMMON.CHAIN'
503 include 'COMMON.DERIV'
504 include 'COMMON.NAMES'
505 include 'COMMON.INTERACT'
506 include 'COMMON.IOUNITS'
507 include 'COMMON.CALC'
509 c double precision rrsave(maxdim)
512 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
514 c if (icall.eq.0) then
526 dxi=dc_norm(1,nres+i)
527 dyi=dc_norm(2,nres+i)
528 dzi=dc_norm(3,nres+i)
529 c dsci_inv=dsc_inv(itypi)
530 dsci_inv=vbld_inv(i+nres)
532 C Calculate SC interaction energy.
535 do j=istart(i,iint),iend(i,iint)
538 c dscj_inv=dsc_inv(itypj)
539 dscj_inv=vbld_inv(j+nres)
540 chi1=chi(itypi,itypj)
541 chi2=chi(itypj,itypi)
548 alf12=0.5D0*(alf1+alf2)
549 C For diagnostics only!!!
562 dxj=dc_norm(1,nres+j)
563 dyj=dc_norm(2,nres+j)
564 dzj=dc_norm(3,nres+j)
565 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
566 cd if (icall.eq.0) then
572 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
574 if (sss.gt.0.0d0) then
576 C Calculate the angle-dependent terms of energy & contributions to derivatives.
578 C Calculate whole angle-dependent part of epsilon and contributions
580 fac=(rrij*sigsq)**expon2
581 e1=fac*fac*aa(itypi,itypj)
582 e2=fac*bb(itypi,itypj)
583 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
584 eps2der=evdwij*eps3rt
585 eps3der=evdwij*eps2rt
586 evdwij=evdwij*eps2rt*eps3rt
589 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
590 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
591 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
592 cd & restyp(itypi),i,restyp(itypj),j,
593 cd & epsi,sigm,chi1,chi2,chip1,chip2,
594 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
595 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
598 C Calculate gradient components.
599 e1=e1*eps1*eps2rt**2*eps3rt**2
600 fac=-expon*(e1+evdwij)
603 C Calculate radial part of the gradient
607 C Calculate the angular part of the gradient and sum add the contributions
608 C to the appropriate components of the Cartesian gradient.
609 call sc_grad_scale(sss)
619 C-----------------------------------------------------------------------------
620 subroutine egb_long(evdw)
622 C This subroutine calculates the interaction energy of nonbonded side chains
623 C assuming the Gay-Berne potential of interaction.
625 implicit real*8 (a-h,o-z)
629 include 'COMMON.LOCAL'
630 include 'COMMON.CHAIN'
631 include 'COMMON.DERIV'
632 include 'COMMON.NAMES'
633 include 'COMMON.INTERACT'
634 include 'COMMON.IOUNITS'
635 include 'COMMON.CALC'
636 include 'COMMON.CONTROL'
639 ccccc energy_dec=.false.
640 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
643 c if (icall.eq.0) lprn=.false.
651 dxi=dc_norm(1,nres+i)
652 dyi=dc_norm(2,nres+i)
653 dzi=dc_norm(3,nres+i)
654 c dsci_inv=dsc_inv(itypi)
655 dsci_inv=vbld_inv(i+nres)
656 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
657 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
659 C Calculate SC interaction energy.
662 do j=istart(i,iint),iend(i,iint)
665 c dscj_inv=dsc_inv(itypj)
666 dscj_inv=vbld_inv(j+nres)
667 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
668 c & 1.0d0/vbld(j+nres)
669 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
670 sig0ij=sigma(itypi,itypj)
671 chi1=chi(itypi,itypj)
672 chi2=chi(itypj,itypi)
679 alf12=0.5D0*(alf1+alf2)
680 C For diagnostics only!!!
693 dxj=dc_norm(1,nres+j)
694 dyj=dc_norm(2,nres+j)
695 dzj=dc_norm(3,nres+j)
696 c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
697 c write (iout,*) "j",j," dc_norm",
698 c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
699 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
701 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
702 c write(iout,*) "long",i,itypi,j,itypj," rij",1.0d0/rij,
703 c & " sigmaii",sigmaii(itypi,itypj)," sss",sss
705 if (sss.lt.1.0d0) then
707 C Calculate angle-dependent terms of energy and contributions to their
711 sig=sig0ij*dsqrt(sigsq)
712 rij_shift=1.0D0/rij-sig+sig0ij
713 c for diagnostics; uncomment
714 c rij_shift=1.2*sig0ij
715 C I hate to put IF's in the loops, but here don't have another choice!!!!
716 if (rij_shift.le.0.0D0) then
718 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
719 cd & restyp(itypi),i,restyp(itypj),j,
720 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
724 c---------------------------------------------------------------
725 rij_shift=1.0D0/rij_shift
727 e1=fac*fac*aa(itypi,itypj)
728 e2=fac*bb(itypi,itypj)
729 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
730 eps2der=evdwij*eps3rt
731 eps3der=evdwij*eps2rt
732 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
733 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
734 evdwij=evdwij*eps2rt*eps3rt
735 evdw=evdw+evdwij*(1.0d0-sss)
736 c write (iout,*) "evdwij",evdwij," evdw",evdw
738 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
739 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
740 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
741 & restyp(itypi),i,restyp(itypj),j,
742 & epsi,sigm,chi1,chi2,chip1,chip2,
743 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
744 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
748 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
751 C Calculate gradient components.
752 e1=e1*eps1*eps2rt**2*eps3rt**2
753 fac=-expon*(e1+evdwij)*rij_shift
757 C Calculate the radial part of the gradient
761 C Calculate angular part of the gradient.
762 call sc_grad_scale(1.0d0-sss)
769 cccc energy_dec=.false.
772 C-----------------------------------------------------------------------------
773 subroutine egb_short(evdw)
775 C This subroutine calculates the interaction energy of nonbonded side chains
776 C assuming the Gay-Berne potential of interaction.
778 implicit real*8 (a-h,o-z)
782 include 'COMMON.LOCAL'
783 include 'COMMON.CHAIN'
784 include 'COMMON.DERIV'
785 include 'COMMON.NAMES'
786 include 'COMMON.INTERACT'
787 include 'COMMON.IOUNITS'
788 include 'COMMON.CALC'
789 include 'COMMON.CONTROL'
792 ccccc energy_dec=.false.
793 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
796 c if (icall.eq.0) lprn=.false.
804 dxi=dc_norm(1,nres+i)
805 dyi=dc_norm(2,nres+i)
806 dzi=dc_norm(3,nres+i)
807 c dsci_inv=dsc_inv(itypi)
808 dsci_inv=vbld_inv(i+nres)
809 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
810 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
812 C Calculate SC interaction energy.
815 do j=istart(i,iint),iend(i,iint)
818 c dscj_inv=dsc_inv(itypj)
819 dscj_inv=vbld_inv(j+nres)
820 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
821 c & 1.0d0/vbld(j+nres)
822 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
823 sig0ij=sigma(itypi,itypj)
824 chi1=chi(itypi,itypj)
825 chi2=chi(itypj,itypi)
832 alf12=0.5D0*(alf1+alf2)
833 C For diagnostics only!!!
846 dxj=dc_norm(1,nres+j)
847 dyj=dc_norm(2,nres+j)
848 dzj=dc_norm(3,nres+j)
849 c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
850 c write (iout,*) "j",j," dc_norm",
851 c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
852 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
854 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
855 c write(iout,*) "short",i,itypi,j,itypj," rij",1.0d0/rij,
856 c & " sigmaii",sigmaii(itypi,itypj)," sss",sss
857 if (sss.gt.0.0d0) then
859 C Calculate angle-dependent terms of energy and contributions to their
863 sig=sig0ij*dsqrt(sigsq)
864 rij_shift=1.0D0/rij-sig+sig0ij
865 c for diagnostics; uncomment
866 c rij_shift=1.2*sig0ij
867 C I hate to put IF's in the loops, but here don't have another choice!!!!
868 if (rij_shift.le.0.0D0) then
870 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
871 cd & restyp(itypi),i,restyp(itypj),j,
872 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
876 c---------------------------------------------------------------
877 rij_shift=1.0D0/rij_shift
879 e1=fac*fac*aa(itypi,itypj)
880 e2=fac*bb(itypi,itypj)
881 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
882 eps2der=evdwij*eps3rt
883 eps3der=evdwij*eps2rt
884 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
885 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
886 evdwij=evdwij*eps2rt*eps3rt
888 c write (iout,*) "evdwij",evdwij," evdw",evdw
890 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
891 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
892 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
893 & restyp(itypi),i,restyp(itypj),j,
894 & epsi,sigm,chi1,chi2,chip1,chip2,
895 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
896 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
900 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
903 C Calculate gradient components.
904 e1=e1*eps1*eps2rt**2*eps3rt**2
905 fac=-expon*(e1+evdwij)*rij_shift
909 C Calculate the radial part of the gradient
913 C Calculate angular part of the gradient.
914 call sc_grad_scale(sss)
921 cccc energy_dec=.false.
924 C-----------------------------------------------------------------------------
925 subroutine egbv_long(evdw)
927 C This subroutine calculates the interaction energy of nonbonded side chains
928 C assuming the Gay-Berne-Vorobjev potential of interaction.
930 implicit real*8 (a-h,o-z)
934 include 'COMMON.LOCAL'
935 include 'COMMON.CHAIN'
936 include 'COMMON.DERIV'
937 include 'COMMON.NAMES'
938 include 'COMMON.INTERACT'
939 include 'COMMON.IOUNITS'
940 include 'COMMON.CALC'
944 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
947 c if (icall.eq.0) lprn=.true.
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)
961 C Calculate SC interaction energy.
964 do j=istart(i,iint),iend(i,iint)
967 c dscj_inv=dsc_inv(itypj)
968 dscj_inv=vbld_inv(j+nres)
969 sig0ij=sigma(itypi,itypj)
971 chi1=chi(itypi,itypj)
972 chi2=chi(itypj,itypi)
979 alf12=0.5D0*(alf1+alf2)
980 C For diagnostics only!!!
993 dxj=dc_norm(1,nres+j)
994 dyj=dc_norm(2,nres+j)
995 dzj=dc_norm(3,nres+j)
996 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
999 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1001 if (sss.lt.1.0d0) then
1003 C Calculate angle-dependent terms of energy and contributions to their
1007 sig=sig0ij*dsqrt(sigsq)
1008 rij_shift=1.0D0/rij-sig+r0ij
1009 C I hate to put IF's in the loops, but here don't have another choice!!!!
1010 if (rij_shift.le.0.0D0) then
1015 c---------------------------------------------------------------
1016 rij_shift=1.0D0/rij_shift
1017 fac=rij_shift**expon
1018 e1=fac*fac*aa(itypi,itypj)
1019 e2=fac*bb(itypi,itypj)
1020 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1021 eps2der=evdwij*eps3rt
1022 eps3der=evdwij*eps2rt
1023 fac_augm=rrij**expon
1024 e_augm=augm(itypi,itypj)*fac_augm
1025 evdwij=evdwij*eps2rt*eps3rt
1026 evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
1028 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1029 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1030 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1031 & restyp(itypi),i,restyp(itypj),j,
1032 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1033 & chi1,chi2,chip1,chip2,
1034 & eps1,eps2rt**2,eps3rt**2,
1035 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1038 C Calculate gradient components.
1039 e1=e1*eps1*eps2rt**2*eps3rt**2
1040 fac=-expon*(e1+evdwij)*rij_shift
1042 fac=rij*fac-2*expon*rrij*e_augm
1043 C Calculate the radial part of the gradient
1047 C Calculate angular part of the gradient.
1048 call sc_grad_scale(1.0d0-sss)
1056 C-----------------------------------------------------------------------------
1057 subroutine egbv_short(evdw)
1059 C This subroutine calculates the interaction energy of nonbonded side chains
1060 C assuming the Gay-Berne-Vorobjev potential of interaction.
1062 implicit real*8 (a-h,o-z)
1063 include 'DIMENSIONS'
1064 include 'COMMON.GEO'
1065 include 'COMMON.VAR'
1066 include 'COMMON.LOCAL'
1067 include 'COMMON.CHAIN'
1068 include 'COMMON.DERIV'
1069 include 'COMMON.NAMES'
1070 include 'COMMON.INTERACT'
1071 include 'COMMON.IOUNITS'
1072 include 'COMMON.CALC'
1073 common /srutu/ icall
1076 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1079 c if (icall.eq.0) lprn=.true.
1081 do i=iatsc_s,iatsc_e
1087 dxi=dc_norm(1,nres+i)
1088 dyi=dc_norm(2,nres+i)
1089 dzi=dc_norm(3,nres+i)
1090 c dsci_inv=dsc_inv(itypi)
1091 dsci_inv=vbld_inv(i+nres)
1093 C Calculate SC interaction energy.
1095 do iint=1,nint_gr(i)
1096 do j=istart(i,iint),iend(i,iint)
1099 c dscj_inv=dsc_inv(itypj)
1100 dscj_inv=vbld_inv(j+nres)
1101 sig0ij=sigma(itypi,itypj)
1102 r0ij=r0(itypi,itypj)
1103 chi1=chi(itypi,itypj)
1104 chi2=chi(itypj,itypi)
1111 alf12=0.5D0*(alf1+alf2)
1112 C For diagnostics only!!!
1125 dxj=dc_norm(1,nres+j)
1126 dyj=dc_norm(2,nres+j)
1127 dzj=dc_norm(3,nres+j)
1128 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1131 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1133 if (sss.gt.0.0d0) then
1135 C Calculate angle-dependent terms of energy and contributions to their
1139 sig=sig0ij*dsqrt(sigsq)
1140 rij_shift=1.0D0/rij-sig+r0ij
1141 C I hate to put IF's in the loops, but here don't have another choice!!!!
1142 if (rij_shift.le.0.0D0) then
1147 c---------------------------------------------------------------
1148 rij_shift=1.0D0/rij_shift
1149 fac=rij_shift**expon
1150 e1=fac*fac*aa(itypi,itypj)
1151 e2=fac*bb(itypi,itypj)
1152 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1153 eps2der=evdwij*eps3rt
1154 eps3der=evdwij*eps2rt
1155 fac_augm=rrij**expon
1156 e_augm=augm(itypi,itypj)*fac_augm
1157 evdwij=evdwij*eps2rt*eps3rt
1158 evdw=evdw+(evdwij+e_augm)*sss
1160 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1161 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1162 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1163 & restyp(itypi),i,restyp(itypj),j,
1164 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1165 & chi1,chi2,chip1,chip2,
1166 & eps1,eps2rt**2,eps3rt**2,
1167 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1170 C Calculate gradient components.
1171 e1=e1*eps1*eps2rt**2*eps3rt**2
1172 fac=-expon*(e1+evdwij)*rij_shift
1174 fac=rij*fac-2*expon*rrij*e_augm
1175 C Calculate the radial part of the gradient
1179 C Calculate angular part of the gradient.
1180 call sc_grad_scale(sss)
1188 C----------------------------------------------------------------------------
1189 subroutine sc_grad_scale(scalfac)
1190 implicit real*8 (a-h,o-z)
1191 include 'DIMENSIONS'
1192 include 'COMMON.CHAIN'
1193 include 'COMMON.DERIV'
1194 include 'COMMON.CALC'
1195 include 'COMMON.IOUNITS'
1196 double precision dcosom1(3),dcosom2(3)
1197 double precision scalfac
1198 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1199 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1200 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1201 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1205 c eom12=evdwij*eps1_om12
1207 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1208 c & " sigder",sigder
1209 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1210 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1212 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1213 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1216 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1218 c write (iout,*) "gg",(gg(k),k=1,3)
1220 gvdwx(k,i)=gvdwx(k,i)-gg(k)
1221 & +((eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1222 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv)*scalfac
1223 gvdwx(k,j)=gvdwx(k,j)+gg(k)
1224 & +((eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1225 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv)*scalfac
1226 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1227 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1228 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1229 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1232 C Calculate the components of the gradient in DC and X
1236 gvdwc(l,k)=gvdwc(l,k)+gg(l)
1241 C--------------------------------------------------------------------------
1242 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1244 C This subroutine calculates the average interaction energy and its gradient
1245 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1246 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1247 C The potential depends both on the distance of peptide-group centers and on
1248 C the orientation of the CA-CA virtual bonds.
1250 implicit real*8 (a-h,o-z)
1251 include 'DIMENSIONS'
1252 include 'COMMON.CONTROL'
1253 include 'COMMON.IOUNITS'
1254 include 'COMMON.GEO'
1255 include 'COMMON.VAR'
1256 include 'COMMON.LOCAL'
1257 include 'COMMON.CHAIN'
1258 include 'COMMON.DERIV'
1259 include 'COMMON.INTERACT'
1260 include 'COMMON.CONTACTS'
1261 include 'COMMON.TORSION'
1262 include 'COMMON.VECTORS'
1263 include 'COMMON.FFIELD'
1264 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1265 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1266 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1267 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1268 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
1269 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1271 double precision scal_el /1.0d0/
1273 double precision scal_el /0.5d0/
1276 C 13-go grudnia roku pamietnego...
1277 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1278 & 0.0d0,1.0d0,0.0d0,
1279 & 0.0d0,0.0d0,1.0d0/
1280 cd write(iout,*) 'In EELEC'
1282 cd write(iout,*) 'Type',i
1283 cd write(iout,*) 'B1',B1(:,i)
1284 cd write(iout,*) 'B2',B2(:,i)
1285 cd write(iout,*) 'CC',CC(:,:,i)
1286 cd write(iout,*) 'DD',DD(:,:,i)
1287 cd write(iout,*) 'EE',EE(:,:,i)
1289 cd call check_vecgrad
1291 if (icheckgrad.eq.1) then
1293 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1295 dc_norm(k,i)=dc(k,i)*fac
1297 c write (iout,*) 'i',i,' fac',fac
1300 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1301 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1302 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1303 c call vec_and_deriv
1307 cd write (iout,*) 'i=',i
1309 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1312 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1313 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1326 cd print '(a)','Enter EELEC'
1327 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1329 gel_loc_loc(i)=0.0d0
1333 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1334 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1337 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1339 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1341 do i=iturn3_start,iturn3_end
1345 dx_normi=dc_norm(1,i)
1346 dy_normi=dc_norm(2,i)
1347 dz_normi=dc_norm(3,i)
1348 xmedi=c(1,i)+0.5d0*dxi
1349 ymedi=c(2,i)+0.5d0*dyi
1350 zmedi=c(3,i)+0.5d0*dzi
1352 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1353 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1354 num_cont_hb(i)=num_conti
1356 do i=iturn4_start,iturn4_end
1360 dx_normi=dc_norm(1,i)
1361 dy_normi=dc_norm(2,i)
1362 dz_normi=dc_norm(3,i)
1363 xmedi=c(1,i)+0.5d0*dxi
1364 ymedi=c(2,i)+0.5d0*dyi
1365 zmedi=c(3,i)+0.5d0*dzi
1367 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1368 if (wturn4.gt.0.0d0) call eturn4(i,eello_turn4)
1369 num_cont_hb(i)=num_cont_hb(i)+num_conti
1372 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1374 do i=iatel_s,iatel_e
1378 dx_normi=dc_norm(1,i)
1379 dy_normi=dc_norm(2,i)
1380 dz_normi=dc_norm(3,i)
1381 xmedi=c(1,i)+0.5d0*dxi
1382 ymedi=c(2,i)+0.5d0*dyi
1383 zmedi=c(3,i)+0.5d0*dzi
1385 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1386 do j=ielstart(i),ielend(i)
1387 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1389 num_cont_hb(i)=num_cont_hb(i)+num_conti
1393 C-------------------------------------------------------------------------------
1394 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1396 C This subroutine calculates the average interaction energy and its gradient
1397 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1398 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1399 C The potential depends both on the distance of peptide-group centers and on
1400 C the orientation of the CA-CA virtual bonds.
1402 implicit real*8 (a-h,o-z)
1403 include 'DIMENSIONS'
1404 include 'COMMON.CONTROL'
1405 include 'COMMON.IOUNITS'
1406 include 'COMMON.GEO'
1407 include 'COMMON.VAR'
1408 include 'COMMON.LOCAL'
1409 include 'COMMON.CHAIN'
1410 include 'COMMON.DERIV'
1411 include 'COMMON.INTERACT'
1412 include 'COMMON.CONTACTS'
1413 include 'COMMON.TORSION'
1414 include 'COMMON.VECTORS'
1415 include 'COMMON.FFIELD'
1416 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1417 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1418 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1419 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1420 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1421 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1423 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1425 double precision scal_el /1.0d0/
1427 double precision scal_el /0.5d0/
1430 C 13-go grudnia roku pamietnego...
1431 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1432 & 0.0d0,1.0d0,0.0d0,
1433 & 0.0d0,0.0d0,1.0d0/
1437 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1438 aaa=app(iteli,itelj)
1439 bbb=bpp(iteli,itelj)
1440 ael6i=ael6(iteli,itelj)
1441 ael3i=ael3(iteli,itelj)
1442 C Diagnostics only!!!
1451 dx_normj=dc_norm(1,j)
1452 dy_normj=dc_norm(2,j)
1453 dz_normj=dc_norm(3,j)
1454 xj=c(1,j)+0.5D0*dxj-xmedi
1455 yj=c(2,j)+0.5D0*dyj-ymedi
1456 zj=c(3,j)+0.5D0*dzj-zmedi
1457 rij=xj*xj+yj*yj+zj*zj
1461 c For extracting the short-range part of Evdwpp
1462 sss=sscale(rij/rpp(iteli,itelj))
1466 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1467 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1468 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1469 fac=cosa-3.0D0*cosb*cosg
1471 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1472 if (j.eq.i+2) ev1=scal_el*ev1
1477 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1480 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1481 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1483 evdw1=evdw1+evdwij*(1.0d0-sss)
1484 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1485 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1486 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1487 cd & xmedi,ymedi,zmedi,xj,yj,zj
1489 if (energy_dec) then
1490 write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
1491 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1495 C Calculate contributions to the Cartesian gradient.
1498 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1499 facel=-3*rrmij*(el1+eesij)
1505 * Radial derivatives. First process both termini of the fragment (i,j)
1512 gelc(k,i)=gelc(k,i)+ghalf
1513 gelc(k,j)=gelc(k,j)+ghalf
1516 * Loop over residues i+1 thru j-1.
1520 gelc(l,k)=gelc(l,k)+ggg(l)
1528 gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1529 gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1532 * Loop over residues i+1 thru j-1.
1536 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1540 facvdw=(ev1+evdwij)*(1.0d0-sss)
1543 fac=-3*rrmij*(facvdw+facvdw+facel)
1548 * Radial derivatives. First process both termini of the fragment (i,j)
1555 gelc(k,i)=gelc(k,i)+ghalf
1556 gelc(k,j)=gelc(k,j)+ghalf
1559 * Loop over residues i+1 thru j-1.
1563 gelc(l,k)=gelc(l,k)+ggg(l)
1570 ecosa=2.0D0*fac3*fac1+fac4
1573 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1574 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1576 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1577 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1579 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1580 cd & (dcosg(k),k=1,3)
1582 ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1586 gelc(k,i)=gelc(k,i)+ghalf
1587 & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1588 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1589 gelc(k,j)=gelc(k,j)+ghalf
1590 & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1591 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1595 gelc(l,k)=gelc(l,k)+ggg(l)
1599 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1600 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1601 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1603 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1604 C energy of a peptide unit is assumed in the form of a second-order
1605 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1606 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1607 C are computed for EVERY pair of non-contiguous peptide groups.
1609 if (j.lt.nres-1) then
1620 muij(kkk)=mu(k,i)*mu(l,j)
1623 cd write (iout,*) 'EELEC: i',i,' j',j
1624 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1625 cd write(iout,*) 'muij',muij
1626 ury=scalar(uy(1,i),erij)
1627 urz=scalar(uz(1,i),erij)
1628 vry=scalar(uy(1,j),erij)
1629 vrz=scalar(uz(1,j),erij)
1630 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1631 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1632 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1633 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1634 C For diagnostics only
1639 fac=dsqrt(-ael6i)*r3ij
1640 cd write (2,*) 'fac=',fac
1641 C For diagnostics only
1647 cd write (iout,'(4i5,4f10.5)')
1648 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1649 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1650 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1651 cd & uy(:,j),uz(:,j)
1652 cd write (iout,'(4f10.5)')
1653 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1654 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1655 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
1656 cd write (iout,'(9f10.5/)')
1657 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1658 C Derivatives of the elements of A in virtual-bond vectors
1659 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1666 uryg(k,1)=scalar(erder(1,k),uy(1,i))
1667 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1668 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1669 urzg(k,1)=scalar(erder(1,k),uz(1,i))
1670 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1671 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1672 vryg(k,1)=scalar(erder(1,k),uy(1,j))
1673 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1674 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1675 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1676 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1677 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1687 C Compute radial contributions to the gradient
1709 C Add the contributions coming from er
1712 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1713 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1714 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1715 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1718 C Derivatives in DC(i)
1719 ghalf1=0.5d0*agg(k,1)
1720 ghalf2=0.5d0*agg(k,2)
1721 ghalf3=0.5d0*agg(k,3)
1722 ghalf4=0.5d0*agg(k,4)
1723 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1724 & -3.0d0*uryg(k,2)*vry)+ghalf1
1725 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1726 & -3.0d0*uryg(k,2)*vrz)+ghalf2
1727 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1728 & -3.0d0*urzg(k,2)*vry)+ghalf3
1729 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1730 & -3.0d0*urzg(k,2)*vrz)+ghalf4
1731 C Derivatives in DC(i+1)
1732 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
1733 & -3.0d0*uryg(k,3)*vry)+agg(k,1)
1734 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
1735 & -3.0d0*uryg(k,3)*vrz)+agg(k,2)
1736 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
1737 & -3.0d0*urzg(k,3)*vry)+agg(k,3)
1738 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
1739 & -3.0d0*urzg(k,3)*vrz)+agg(k,4)
1740 C Derivatives in DC(j)
1741 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
1742 & -3.0d0*vryg(k,2)*ury)+ghalf1
1743 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
1744 & -3.0d0*vrzg(k,2)*ury)+ghalf2
1745 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
1746 & -3.0d0*vryg(k,2)*urz)+ghalf3
1747 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
1748 & -3.0d0*vrzg(k,2)*urz)+ghalf4
1749 C Derivatives in DC(j+1) or DC(nres-1)
1750 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
1751 & -3.0d0*vryg(k,3)*ury)
1752 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
1753 & -3.0d0*vrzg(k,3)*ury)
1754 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
1755 & -3.0d0*vryg(k,3)*urz)
1756 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
1757 & -3.0d0*vrzg(k,3)*urz)
1762 C Derivatives in DC(i+1)
1763 cd aggi1(k,1)=agg(k,1)
1764 cd aggi1(k,2)=agg(k,2)
1765 cd aggi1(k,3)=agg(k,3)
1766 cd aggi1(k,4)=agg(k,4)
1767 C Derivatives in DC(j)
1772 C Derivatives in DC(j+1)
1777 if (j.eq.nres-1 .and. i.lt.j-2) then
1779 aggj1(k,l)=aggj1(k,l)+agg(k,l)
1780 cd aggj1(k,l)=agg(k,l)
1785 C Check the loc-el terms by numerical integration
1795 aggi(k,l)=-aggi(k,l)
1796 aggi1(k,l)=-aggi1(k,l)
1797 aggj(k,l)=-aggj(k,l)
1798 aggj1(k,l)=-aggj1(k,l)
1801 if (j.lt.nres-1) then
1807 aggi(k,l)=-aggi(k,l)
1808 aggi1(k,l)=-aggi1(k,l)
1809 aggj(k,l)=-aggj(k,l)
1810 aggj1(k,l)=-aggj1(k,l)
1821 aggi(k,l)=-aggi(k,l)
1822 aggi1(k,l)=-aggi1(k,l)
1823 aggj(k,l)=-aggj(k,l)
1824 aggj1(k,l)=-aggj1(k,l)
1830 IF (wel_loc.gt.0.0d0) THEN
1831 C Contribution to the local-electrostatic energy coming from the i-j pair
1832 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
1834 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
1836 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1837 & 'eelloc',i,j,eel_loc_ij
1839 eel_loc=eel_loc+eel_loc_ij
1840 C Partial derivatives in virtual-bond dihedral angles gamma
1842 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
1843 & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1844 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1845 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
1846 & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1847 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1848 cd call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij)
1849 cd write(iout,*) 'agg ',agg
1850 cd write(iout,*) 'aggi ',aggi
1851 cd write(iout,*) 'aggi1',aggi1
1852 cd write(iout,*) 'aggj ',aggj
1853 cd write(iout,*) 'aggj1',aggj1
1855 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
1857 ggg(l)=agg(l,1)*muij(1)+
1858 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1862 gel_loc(l,k)=gel_loc(l,k)+ggg(l)
1865 C Remaining derivatives of eello
1867 gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1868 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1869 gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1870 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1871 gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1872 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1873 gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1874 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1877 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
1878 if (j.gt.i+1 .and. num_conti.le.maxconts) then
1880 C Calculate the contact function. The ith column of the array JCONT will
1881 C contain the numbers of atoms that make contacts with the atom I (of numbers
1882 C greater than I). The arrays FACONT and GACONT will contain the values of
1883 C the contact function and its derivative.
1884 c r0ij=1.02D0*rpp(iteli,itelj)
1885 c r0ij=1.11D0*rpp(iteli,itelj)
1886 r0ij=2.20D0*rpp(iteli,itelj)
1887 c r0ij=1.55D0*rpp(iteli,itelj)
1888 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
1889 if (fcont.gt.0.0D0) then
1890 num_conti=num_conti+1
1891 if (num_conti.gt.maxconts) then
1892 write (iout,*) 'WARNING - max. # of contacts exceeded;',
1893 & ' will skip next contacts for this conf.'
1895 jcont_hb(num_conti,i)=j
1896 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
1897 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1898 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
1900 d_cont(num_conti,i)=rij
1901 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
1902 C --- Electrostatic-interaction matrix ---
1903 a_chuj(1,1,num_conti,i)=a22
1904 a_chuj(1,2,num_conti,i)=a23
1905 a_chuj(2,1,num_conti,i)=a32
1906 a_chuj(2,2,num_conti,i)=a33
1907 C --- Gradient of rij
1909 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
1912 c a_chuj(1,1,num_conti,i)=-0.61d0
1913 c a_chuj(1,2,num_conti,i)= 0.4d0
1914 c a_chuj(2,1,num_conti,i)= 0.65d0
1915 c a_chuj(2,2,num_conti,i)= 0.50d0
1916 c else if (i.eq.2) then
1917 c a_chuj(1,1,num_conti,i)= 0.0d0
1918 c a_chuj(1,2,num_conti,i)= 0.0d0
1919 c a_chuj(2,1,num_conti,i)= 0.0d0
1920 c a_chuj(2,2,num_conti,i)= 0.0d0
1922 C --- and its gradients
1923 cd write (iout,*) 'i',i,' j',j
1925 cd write (iout,*) 'iii 1 kkk',kkk
1926 cd write (iout,*) agg(kkk,:)
1929 cd write (iout,*) 'iii 2 kkk',kkk
1930 cd write (iout,*) aggi(kkk,:)
1933 cd write (iout,*) 'iii 3 kkk',kkk
1934 cd write (iout,*) aggi1(kkk,:)
1937 cd write (iout,*) 'iii 4 kkk',kkk
1938 cd write (iout,*) aggj(kkk,:)
1941 cd write (iout,*) 'iii 5 kkk',kkk
1942 cd write (iout,*) aggj1(kkk,:)
1949 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
1950 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
1951 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
1952 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
1953 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
1955 c a_chuj_der(k,l,m,mm,num_conti,i)=0.0d0
1961 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
1962 C Calculate contact energies
1964 wij=cosa-3.0D0*cosb*cosg
1967 c fac3=dsqrt(-ael6i)/r0ij**3
1968 fac3=dsqrt(-ael6i)*r3ij
1969 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
1970 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
1971 if (ees0tmp.gt.0) then
1972 ees0pij=dsqrt(ees0tmp)
1976 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
1977 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
1978 if (ees0tmp.gt.0) then
1979 ees0mij=dsqrt(ees0tmp)
1984 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
1985 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
1986 C Diagnostics. Comment out or remove after debugging!
1987 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
1988 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
1989 c ees0m(num_conti,i)=0.0D0
1991 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
1992 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
1993 C Angular derivatives of the contact function
1994 ees0pij1=fac3/ees0pij
1995 ees0mij1=fac3/ees0mij
1996 fac3p=-3.0D0*fac3*rrmij
1997 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
1998 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2000 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2001 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2002 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2003 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2004 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2005 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2006 ecosap=ecosa1+ecosa2
2007 ecosbp=ecosb1+ecosb2
2008 ecosgp=ecosg1+ecosg2
2009 ecosam=ecosa1-ecosa2
2010 ecosbm=ecosb1-ecosb2
2011 ecosgm=ecosg1-ecosg2
2020 facont_hb(num_conti,i)=fcont
2021 fprimcont=fprimcont/rij
2022 cd facont_hb(num_conti,i)=1.0D0
2023 C Following line is for diagnostics.
2026 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2027 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2030 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2031 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2033 gggp(1)=gggp(1)+ees0pijp*xj
2034 gggp(2)=gggp(2)+ees0pijp*yj
2035 gggp(3)=gggp(3)+ees0pijp*zj
2036 gggm(1)=gggm(1)+ees0mijp*xj
2037 gggm(2)=gggm(2)+ees0mijp*yj
2038 gggm(3)=gggm(3)+ees0mijp*zj
2039 C Derivatives due to the contact function
2040 gacont_hbr(1,num_conti,i)=fprimcont*xj
2041 gacont_hbr(2,num_conti,i)=fprimcont*yj
2042 gacont_hbr(3,num_conti,i)=fprimcont*zj
2044 ghalfp=0.5D0*gggp(k)
2045 ghalfm=0.5D0*gggm(k)
2046 gacontp_hb1(k,num_conti,i)=ghalfp
2047 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2048 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2049 gacontp_hb2(k,num_conti,i)=ghalfp
2050 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2051 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2052 gacontp_hb3(k,num_conti,i)=gggp(k)
2053 gacontm_hb1(k,num_conti,i)=ghalfm
2054 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2055 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2056 gacontm_hb2(k,num_conti,i)=ghalfm
2057 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2058 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2059 gacontm_hb3(k,num_conti,i)=gggm(k)
2061 C Diagnostics. Comment out or remove after debugging!
2063 cdiag gacontp_hb1(k,num_conti,i)=0.0D0
2064 cdiag gacontp_hb2(k,num_conti,i)=0.0D0
2065 cdiag gacontp_hb3(k,num_conti,i)=0.0D0
2066 cdiag gacontm_hb1(k,num_conti,i)=0.0D0
2067 cdiag gacontm_hb2(k,num_conti,i)=0.0D0
2068 cdiag gacontm_hb3(k,num_conti,i)=0.0D0
2071 endif ! num_conti.le.maxconts
2076 C-----------------------------------------------------------------------
2077 subroutine evdwpp_long(evdw1)
2081 implicit real*8 (a-h,o-z)
2082 include 'DIMENSIONS'
2083 include 'COMMON.CONTROL'
2084 include 'COMMON.IOUNITS'
2085 include 'COMMON.GEO'
2086 include 'COMMON.VAR'
2087 include 'COMMON.LOCAL'
2088 include 'COMMON.CHAIN'
2089 include 'COMMON.DERIV'
2090 include 'COMMON.INTERACT'
2091 include 'COMMON.CONTACTS'
2092 include 'COMMON.TORSION'
2093 include 'COMMON.VECTORS'
2094 include 'COMMON.FFIELD'
2096 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2098 double precision scal_el /1.0d0/
2100 double precision scal_el /0.5d0/
2103 do i=iatel_s,iatel_e
2107 dx_normi=dc_norm(1,i)
2108 dy_normi=dc_norm(2,i)
2109 dz_normi=dc_norm(3,i)
2110 xmedi=c(1,i)+0.5d0*dxi
2111 ymedi=c(2,i)+0.5d0*dyi
2112 zmedi=c(3,i)+0.5d0*dzi
2114 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
2115 do j=ielstart(i),ielend(i)
2119 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2120 aaa=app(iteli,itelj)
2121 bbb=bpp(iteli,itelj)
2125 dx_normj=dc_norm(1,j)
2126 dy_normj=dc_norm(2,j)
2127 dz_normj=dc_norm(3,j)
2128 xj=c(1,j)+0.5D0*dxj-xmedi
2129 yj=c(2,j)+0.5D0*dyj-ymedi
2130 zj=c(3,j)+0.5D0*dzj-zmedi
2131 rij=xj*xj+yj*yj+zj*zj
2134 sss=sscale(rij/rpp(iteli,itelj))
2135 if (sss.lt.1.0d0) then
2140 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2141 if (j.eq.i+2) ev1=scal_el*ev1
2144 if (energy_dec) then
2145 write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
2146 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
2148 evdw1=evdw1+evdwij*(1.0d0-sss)
2150 C Calculate contributions to the Cartesian gradient.
2152 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
2159 gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2160 gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2163 * Loop over residues i+1 thru j-1.
2167 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2175 C-----------------------------------------------------------------------
2176 subroutine evdwpp_short(evdw1)
2180 implicit real*8 (a-h,o-z)
2181 include 'DIMENSIONS'
2182 include 'COMMON.CONTROL'
2183 include 'COMMON.IOUNITS'
2184 include 'COMMON.GEO'
2185 include 'COMMON.VAR'
2186 include 'COMMON.LOCAL'
2187 include 'COMMON.CHAIN'
2188 include 'COMMON.DERIV'
2189 include 'COMMON.INTERACT'
2190 include 'COMMON.CONTACTS'
2191 include 'COMMON.TORSION'
2192 include 'COMMON.VECTORS'
2193 include 'COMMON.FFIELD'
2195 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2197 double precision scal_el /1.0d0/
2199 double precision scal_el /0.5d0/
2202 do i=iatel_s,iatel_e
2206 dx_normi=dc_norm(1,i)
2207 dy_normi=dc_norm(2,i)
2208 dz_normi=dc_norm(3,i)
2209 xmedi=c(1,i)+0.5d0*dxi
2210 ymedi=c(2,i)+0.5d0*dyi
2211 zmedi=c(3,i)+0.5d0*dzi
2213 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
2214 do j=ielstart(i),ielend(i)
2218 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2219 aaa=app(iteli,itelj)
2220 bbb=bpp(iteli,itelj)
2224 dx_normj=dc_norm(1,j)
2225 dy_normj=dc_norm(2,j)
2226 dz_normj=dc_norm(3,j)
2227 xj=c(1,j)+0.5D0*dxj-xmedi
2228 yj=c(2,j)+0.5D0*dyj-ymedi
2229 zj=c(3,j)+0.5D0*dzj-zmedi
2230 rij=xj*xj+yj*yj+zj*zj
2233 sss=sscale(rij/rpp(iteli,itelj))
2234 if (sss.gt.0.0d0) then
2239 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2240 if (j.eq.i+2) ev1=scal_el*ev1
2243 if (energy_dec) then
2244 write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
2245 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
2247 evdw1=evdw1+evdwij*sss
2249 C Calculate contributions to the Cartesian gradient.
2251 facvdw=-6*rrmij*(ev1+evdwij)*sss
2258 gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2259 gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2262 * Loop over residues i+1 thru j-1.
2266 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2274 C-----------------------------------------------------------------------------
2275 subroutine escp_long(evdw2,evdw2_14)
2277 C This subroutine calculates the excluded-volume interaction energy between
2278 C peptide-group centers and side chains and its gradient in virtual-bond and
2279 C side-chain vectors.
2281 implicit real*8 (a-h,o-z)
2282 include 'DIMENSIONS'
2283 include 'COMMON.GEO'
2284 include 'COMMON.VAR'
2285 include 'COMMON.LOCAL'
2286 include 'COMMON.CHAIN'
2287 include 'COMMON.DERIV'
2288 include 'COMMON.INTERACT'
2289 include 'COMMON.FFIELD'
2290 include 'COMMON.IOUNITS'
2291 include 'COMMON.CONTROL'
2295 cd print '(a)','Enter ESCP'
2296 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2297 do i=iatscp_s,iatscp_e
2299 xi=0.5D0*(c(1,i)+c(1,i+1))
2300 yi=0.5D0*(c(2,i)+c(2,i+1))
2301 zi=0.5D0*(c(3,i)+c(3,i+1))
2303 do iint=1,nscp_gr(i)
2305 do j=iscpstart(i,iint),iscpend(i,iint)
2307 C Uncomment following three lines for SC-p interactions
2311 C Uncomment following three lines for Ca-p interactions
2315 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2317 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2319 if (sss.lt.1.0d0) then
2322 e1=fac*fac*aad(itypj,iteli)
2323 e2=fac*bad(itypj,iteli)
2324 if (iabs(j-i) .le. 2) then
2327 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2330 evdw2=evdw2+evdwij*(1.0d0-sss)
2331 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2332 & 'evdw2',i,j,evdwij
2334 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2336 fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2341 cd write (iout,*) 'j<i'
2342 C Uncomment following three lines for SC-p interactions
2344 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2347 cd write (iout,*) 'j>i'
2350 C Uncomment following line for SC-p interactions
2351 c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
2355 gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
2359 cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
2360 cd write (iout,*) ggg(1),ggg(2),ggg(3)
2363 gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
2375 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2376 gradx_scp(j,i)=expon*gradx_scp(j,i)
2379 C******************************************************************************
2383 C To save time the factor EXPON has been extracted from ALL components
2384 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2387 C******************************************************************************
2390 C-----------------------------------------------------------------------------
2391 subroutine escp_short(evdw2,evdw2_14)
2393 C This subroutine calculates the excluded-volume interaction energy between
2394 C peptide-group centers and side chains and its gradient in virtual-bond and
2395 C side-chain vectors.
2397 implicit real*8 (a-h,o-z)
2398 include 'DIMENSIONS'
2399 include 'COMMON.GEO'
2400 include 'COMMON.VAR'
2401 include 'COMMON.LOCAL'
2402 include 'COMMON.CHAIN'
2403 include 'COMMON.DERIV'
2404 include 'COMMON.INTERACT'
2405 include 'COMMON.FFIELD'
2406 include 'COMMON.IOUNITS'
2407 include 'COMMON.CONTROL'
2411 cd print '(a)','Enter ESCP'
2412 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2413 do i=iatscp_s,iatscp_e
2415 xi=0.5D0*(c(1,i)+c(1,i+1))
2416 yi=0.5D0*(c(2,i)+c(2,i+1))
2417 zi=0.5D0*(c(3,i)+c(3,i+1))
2419 do iint=1,nscp_gr(i)
2421 do j=iscpstart(i,iint),iscpend(i,iint)
2423 C Uncomment following three lines for SC-p interactions
2427 C Uncomment following three lines for Ca-p interactions
2431 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2433 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2435 if (sss.gt.0.0d0) then
2438 e1=fac*fac*aad(itypj,iteli)
2439 e2=fac*bad(itypj,iteli)
2440 if (iabs(j-i) .le. 2) then
2443 evdw2_14=evdw2_14+(e1+e2)*sss
2446 evdw2=evdw2+evdwij*sss
2447 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2448 & 'evdw2',i,j,evdwij
2450 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2452 fac=-(evdwij+e1)*rrij*sss
2457 cd write (iout,*) 'j<i'
2458 C Uncomment following three lines for SC-p interactions
2460 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2463 cd write (iout,*) 'j>i'
2466 C Uncomment following line for SC-p interactions
2467 c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
2471 gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
2475 cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
2476 cd write (iout,*) ggg(1),ggg(2),ggg(3)
2479 gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
2491 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2492 gradx_scp(j,i)=expon*gradx_scp(j,i)
2495 C******************************************************************************
2499 C To save time the factor EXPON has been extracted from ALL components
2500 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2503 C******************************************************************************