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 write(iout,*) "long",i,itypi,j,itypj," rij",1.0d0/rij,
703 & " 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 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 write(iout,*) "short",i,itypi,j,itypj," rij",1.0d0/rij,
856 & " 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 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
1332 do i=iatel_s,iatel_e
1336 dx_normi=dc_norm(1,i)
1337 dy_normi=dc_norm(2,i)
1338 dz_normi=dc_norm(3,i)
1339 xmedi=c(1,i)+0.5d0*dxi
1340 ymedi=c(2,i)+0.5d0*dyi
1341 zmedi=c(3,i)+0.5d0*dzi
1343 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1344 do j=ielstart(i),ielend(i)
1348 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1349 aaa=app(iteli,itelj)
1350 bbb=bpp(iteli,itelj)
1351 ael6i=ael6(iteli,itelj)
1352 ael3i=ael3(iteli,itelj)
1353 C Diagnostics only!!!
1362 dx_normj=dc_norm(1,j)
1363 dy_normj=dc_norm(2,j)
1364 dz_normj=dc_norm(3,j)
1365 xj=c(1,j)+0.5D0*dxj-xmedi
1366 yj=c(2,j)+0.5D0*dyj-ymedi
1367 zj=c(3,j)+0.5D0*dzj-zmedi
1368 rij=xj*xj+yj*yj+zj*zj
1372 c For extracting the short-range part of Evdwpp
1373 sss=sscale(rij/rpp(iteli,itelj))
1377 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1378 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1379 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1380 fac=cosa-3.0D0*cosb*cosg
1382 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1383 if (j.eq.i+2) ev1=scal_el*ev1
1388 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1391 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1392 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1394 evdw1=evdw1+evdwij*(1.0d0-sss)
1395 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1396 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1397 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1398 cd & xmedi,ymedi,zmedi,xj,yj,zj
1400 if (energy_dec) then
1401 write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
1402 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1406 C Calculate contributions to the Cartesian gradient.
1409 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1410 facel=-3*rrmij*(el1+eesij)
1416 * Radial derivatives. First process both termini of the fragment (i,j)
1423 gelc(k,i)=gelc(k,i)+ghalf
1424 gelc(k,j)=gelc(k,j)+ghalf
1427 * Loop over residues i+1 thru j-1.
1431 gelc(l,k)=gelc(l,k)+ggg(l)
1439 gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1440 gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1443 * Loop over residues i+1 thru j-1.
1447 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1451 facvdw=(ev1+evdwij)*(1.0d0-sss)
1454 fac=-3*rrmij*(facvdw+facvdw+facel)
1459 * Radial derivatives. First process both termini of the fragment (i,j)
1466 gelc(k,i)=gelc(k,i)+ghalf
1467 gelc(k,j)=gelc(k,j)+ghalf
1470 * Loop over residues i+1 thru j-1.
1474 gelc(l,k)=gelc(l,k)+ggg(l)
1481 ecosa=2.0D0*fac3*fac1+fac4
1484 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1485 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1487 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1488 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1490 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1491 cd & (dcosg(k),k=1,3)
1493 ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1497 gelc(k,i)=gelc(k,i)+ghalf
1498 & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1499 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1500 gelc(k,j)=gelc(k,j)+ghalf
1501 & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1502 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1506 gelc(l,k)=gelc(l,k)+ggg(l)
1510 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1511 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1512 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1514 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1515 C energy of a peptide unit is assumed in the form of a second-order
1516 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1517 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1518 C are computed for EVERY pair of non-contiguous peptide groups.
1520 if (j.lt.nres-1) then
1531 muij(kkk)=mu(k,i)*mu(l,j)
1534 cd write (iout,*) 'EELEC: i',i,' j',j
1535 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1536 cd write(iout,*) 'muij',muij
1537 ury=scalar(uy(1,i),erij)
1538 urz=scalar(uz(1,i),erij)
1539 vry=scalar(uy(1,j),erij)
1540 vrz=scalar(uz(1,j),erij)
1541 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1542 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1543 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1544 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1545 C For diagnostics only
1550 fac=dsqrt(-ael6i)*r3ij
1551 cd write (2,*) 'fac=',fac
1552 C For diagnostics only
1558 cd write (iout,'(4i5,4f10.5)')
1559 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1560 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1561 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1562 cd & uy(:,j),uz(:,j)
1563 cd write (iout,'(4f10.5)')
1564 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1565 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1566 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
1567 cd write (iout,'(9f10.5/)')
1568 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1569 C Derivatives of the elements of A in virtual-bond vectors
1570 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1577 uryg(k,1)=scalar(erder(1,k),uy(1,i))
1578 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1579 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1580 urzg(k,1)=scalar(erder(1,k),uz(1,i))
1581 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1582 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1583 vryg(k,1)=scalar(erder(1,k),uy(1,j))
1584 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1585 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1586 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1587 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1588 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1598 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 ghalf1=0.5d0*agg(k,1)
1631 ghalf2=0.5d0*agg(k,2)
1632 ghalf3=0.5d0*agg(k,3)
1633 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)
1673 C Derivatives in DC(i+1)
1674 cd aggi1(k,1)=agg(k,1)
1675 cd aggi1(k,2)=agg(k,2)
1676 cd aggi1(k,3)=agg(k,3)
1677 cd aggi1(k,4)=agg(k,4)
1678 C Derivatives in DC(j)
1683 C Derivatives in DC(j+1)
1688 if (j.eq.nres-1 .and. i.lt.j-2) then
1690 aggj1(k,l)=aggj1(k,l)+agg(k,l)
1691 cd aggj1(k,l)=agg(k,l)
1696 C Check the loc-el terms by numerical integration
1706 aggi(k,l)=-aggi(k,l)
1707 aggi1(k,l)=-aggi1(k,l)
1708 aggj(k,l)=-aggj(k,l)
1709 aggj1(k,l)=-aggj1(k,l)
1712 if (j.lt.nres-1) then
1718 aggi(k,l)=-aggi(k,l)
1719 aggi1(k,l)=-aggi1(k,l)
1720 aggj(k,l)=-aggj(k,l)
1721 aggj1(k,l)=-aggj1(k,l)
1732 aggi(k,l)=-aggi(k,l)
1733 aggi1(k,l)=-aggi1(k,l)
1734 aggj(k,l)=-aggj(k,l)
1735 aggj1(k,l)=-aggj1(k,l)
1741 IF (wel_loc.gt.0.0d0) THEN
1742 C Contribution to the local-electrostatic energy coming from the i-j pair
1743 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
1745 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
1747 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1748 & 'eelloc',i,j,eel_loc_ij
1750 eel_loc=eel_loc+eel_loc_ij
1751 C Partial derivatives in virtual-bond dihedral angles gamma
1753 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
1754 & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1755 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1756 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
1757 & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1758 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1759 cd call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij)
1760 cd write(iout,*) 'agg ',agg
1761 cd write(iout,*) 'aggi ',aggi
1762 cd write(iout,*) 'aggi1',aggi1
1763 cd write(iout,*) 'aggj ',aggj
1764 cd write(iout,*) 'aggj1',aggj1
1766 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
1768 ggg(l)=agg(l,1)*muij(1)+
1769 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1773 gel_loc(l,k)=gel_loc(l,k)+ggg(l)
1776 C Remaining derivatives of eello
1778 gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1779 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1780 gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1781 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1782 gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1783 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1784 gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1785 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1788 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
1789 C Contributions from turns
1794 call eturn34(i,j,eello_turn3,eello_turn4)
1796 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
1797 if (j.gt.i+1 .and. num_conti.le.maxconts) then
1799 C Calculate the contact function. The ith column of the array JCONT will
1800 C contain the numbers of atoms that make contacts with the atom I (of numbers
1801 C greater than I). The arrays FACONT and GACONT will contain the values of
1802 C the contact function and its derivative.
1803 c r0ij=1.02D0*rpp(iteli,itelj)
1804 c r0ij=1.11D0*rpp(iteli,itelj)
1805 r0ij=2.20D0*rpp(iteli,itelj)
1806 c r0ij=1.55D0*rpp(iteli,itelj)
1807 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
1808 if (fcont.gt.0.0D0) then
1809 num_conti=num_conti+1
1810 if (num_conti.gt.maxconts) then
1811 write (iout,*) 'WARNING - max. # of contacts exceeded;',
1812 & ' will skip next contacts for this conf.'
1814 jcont_hb(num_conti,i)=j
1815 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
1816 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1817 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
1819 d_cont(num_conti,i)=rij
1820 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
1821 C --- Electrostatic-interaction matrix ---
1822 a_chuj(1,1,num_conti,i)=a22
1823 a_chuj(1,2,num_conti,i)=a23
1824 a_chuj(2,1,num_conti,i)=a32
1825 a_chuj(2,2,num_conti,i)=a33
1826 C --- Gradient of rij
1828 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
1831 c a_chuj(1,1,num_conti,i)=-0.61d0
1832 c a_chuj(1,2,num_conti,i)= 0.4d0
1833 c a_chuj(2,1,num_conti,i)= 0.65d0
1834 c a_chuj(2,2,num_conti,i)= 0.50d0
1835 c else if (i.eq.2) then
1836 c a_chuj(1,1,num_conti,i)= 0.0d0
1837 c a_chuj(1,2,num_conti,i)= 0.0d0
1838 c a_chuj(2,1,num_conti,i)= 0.0d0
1839 c a_chuj(2,2,num_conti,i)= 0.0d0
1841 C --- and its gradients
1842 cd write (iout,*) 'i',i,' j',j
1844 cd write (iout,*) 'iii 1 kkk',kkk
1845 cd write (iout,*) agg(kkk,:)
1848 cd write (iout,*) 'iii 2 kkk',kkk
1849 cd write (iout,*) aggi(kkk,:)
1852 cd write (iout,*) 'iii 3 kkk',kkk
1853 cd write (iout,*) aggi1(kkk,:)
1856 cd write (iout,*) 'iii 4 kkk',kkk
1857 cd write (iout,*) aggj(kkk,:)
1860 cd write (iout,*) 'iii 5 kkk',kkk
1861 cd write (iout,*) aggj1(kkk,:)
1868 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
1869 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
1870 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
1871 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
1872 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
1874 c a_chuj_der(k,l,m,mm,num_conti,i)=0.0d0
1880 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
1881 C Calculate contact energies
1883 wij=cosa-3.0D0*cosb*cosg
1886 c fac3=dsqrt(-ael6i)/r0ij**3
1887 fac3=dsqrt(-ael6i)*r3ij
1888 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
1889 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
1890 if (ees0tmp.gt.0) then
1891 ees0pij=dsqrt(ees0tmp)
1895 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
1896 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
1897 if (ees0tmp.gt.0) then
1898 ees0mij=dsqrt(ees0tmp)
1903 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
1904 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
1905 C Diagnostics. Comment out or remove after debugging!
1906 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
1907 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
1908 c ees0m(num_conti,i)=0.0D0
1910 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
1911 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
1912 C Angular derivatives of the contact function
1913 ees0pij1=fac3/ees0pij
1914 ees0mij1=fac3/ees0mij
1915 fac3p=-3.0D0*fac3*rrmij
1916 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
1917 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
1919 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
1920 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
1921 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
1922 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
1923 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
1924 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
1925 ecosap=ecosa1+ecosa2
1926 ecosbp=ecosb1+ecosb2
1927 ecosgp=ecosg1+ecosg2
1928 ecosam=ecosa1-ecosa2
1929 ecosbm=ecosb1-ecosb2
1930 ecosgm=ecosg1-ecosg2
1939 facont_hb(num_conti,i)=fcont
1940 fprimcont=fprimcont/rij
1941 cd facont_hb(num_conti,i)=1.0D0
1942 C Following line is for diagnostics.
1945 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1946 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1949 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
1950 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
1952 gggp(1)=gggp(1)+ees0pijp*xj
1953 gggp(2)=gggp(2)+ees0pijp*yj
1954 gggp(3)=gggp(3)+ees0pijp*zj
1955 gggm(1)=gggm(1)+ees0mijp*xj
1956 gggm(2)=gggm(2)+ees0mijp*yj
1957 gggm(3)=gggm(3)+ees0mijp*zj
1958 C Derivatives due to the contact function
1959 gacont_hbr(1,num_conti,i)=fprimcont*xj
1960 gacont_hbr(2,num_conti,i)=fprimcont*yj
1961 gacont_hbr(3,num_conti,i)=fprimcont*zj
1963 ghalfp=0.5D0*gggp(k)
1964 ghalfm=0.5D0*gggm(k)
1965 gacontp_hb1(k,num_conti,i)=ghalfp
1966 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
1967 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1968 gacontp_hb2(k,num_conti,i)=ghalfp
1969 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
1970 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1971 gacontp_hb3(k,num_conti,i)=gggp(k)
1972 gacontm_hb1(k,num_conti,i)=ghalfm
1973 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
1974 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1975 gacontm_hb2(k,num_conti,i)=ghalfm
1976 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
1977 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1978 gacontm_hb3(k,num_conti,i)=gggm(k)
1980 C Diagnostics. Comment out or remove after debugging!
1982 cdiag gacontp_hb1(k,num_conti,i)=0.0D0
1983 cdiag gacontp_hb2(k,num_conti,i)=0.0D0
1984 cdiag gacontp_hb3(k,num_conti,i)=0.0D0
1985 cdiag gacontm_hb1(k,num_conti,i)=0.0D0
1986 cdiag gacontm_hb2(k,num_conti,i)=0.0D0
1987 cdiag gacontm_hb3(k,num_conti,i)=0.0D0
1990 endif ! num_conti.le.maxconts
1994 num_cont_hb(i)=num_conti
1997 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1998 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
2000 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
2001 ccc eel_loc=eel_loc+eello_turn3
2004 C-----------------------------------------------------------------------
2005 subroutine evdwpp_long(evdw1)
2009 implicit real*8 (a-h,o-z)
2010 include 'DIMENSIONS'
2011 include 'COMMON.CONTROL'
2012 include 'COMMON.IOUNITS'
2013 include 'COMMON.GEO'
2014 include 'COMMON.VAR'
2015 include 'COMMON.LOCAL'
2016 include 'COMMON.CHAIN'
2017 include 'COMMON.DERIV'
2018 include 'COMMON.INTERACT'
2019 include 'COMMON.CONTACTS'
2020 include 'COMMON.TORSION'
2021 include 'COMMON.VECTORS'
2022 include 'COMMON.FFIELD'
2024 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2026 double precision scal_el /1.0d0/
2028 double precision scal_el /0.5d0/
2031 do i=iatel_s,iatel_e
2035 dx_normi=dc_norm(1,i)
2036 dy_normi=dc_norm(2,i)
2037 dz_normi=dc_norm(3,i)
2038 xmedi=c(1,i)+0.5d0*dxi
2039 ymedi=c(2,i)+0.5d0*dyi
2040 zmedi=c(3,i)+0.5d0*dzi
2042 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
2043 do j=ielstart(i),ielend(i)
2047 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2048 aaa=app(iteli,itelj)
2049 bbb=bpp(iteli,itelj)
2053 dx_normj=dc_norm(1,j)
2054 dy_normj=dc_norm(2,j)
2055 dz_normj=dc_norm(3,j)
2056 xj=c(1,j)+0.5D0*dxj-xmedi
2057 yj=c(2,j)+0.5D0*dyj-ymedi
2058 zj=c(3,j)+0.5D0*dzj-zmedi
2059 rij=xj*xj+yj*yj+zj*zj
2062 sss=sscale(rij/rpp(iteli,itelj))
2063 if (sss.lt.1.0d0) then
2068 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2069 if (j.eq.i+2) ev1=scal_el*ev1
2072 if (energy_dec) then
2073 write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
2074 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
2076 evdw1=evdw1+evdwij*(1.0d0-sss)
2078 C Calculate contributions to the Cartesian gradient.
2080 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
2087 gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2088 gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2091 * Loop over residues i+1 thru j-1.
2095 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2103 C-----------------------------------------------------------------------
2104 subroutine evdwpp_short(evdw1)
2108 implicit real*8 (a-h,o-z)
2109 include 'DIMENSIONS'
2110 include 'COMMON.CONTROL'
2111 include 'COMMON.IOUNITS'
2112 include 'COMMON.GEO'
2113 include 'COMMON.VAR'
2114 include 'COMMON.LOCAL'
2115 include 'COMMON.CHAIN'
2116 include 'COMMON.DERIV'
2117 include 'COMMON.INTERACT'
2118 include 'COMMON.CONTACTS'
2119 include 'COMMON.TORSION'
2120 include 'COMMON.VECTORS'
2121 include 'COMMON.FFIELD'
2123 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2125 double precision scal_el /1.0d0/
2127 double precision scal_el /0.5d0/
2130 do i=iatel_s,iatel_e
2134 dx_normi=dc_norm(1,i)
2135 dy_normi=dc_norm(2,i)
2136 dz_normi=dc_norm(3,i)
2137 xmedi=c(1,i)+0.5d0*dxi
2138 ymedi=c(2,i)+0.5d0*dyi
2139 zmedi=c(3,i)+0.5d0*dzi
2141 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
2142 do j=ielstart(i),ielend(i)
2146 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2147 aaa=app(iteli,itelj)
2148 bbb=bpp(iteli,itelj)
2152 dx_normj=dc_norm(1,j)
2153 dy_normj=dc_norm(2,j)
2154 dz_normj=dc_norm(3,j)
2155 xj=c(1,j)+0.5D0*dxj-xmedi
2156 yj=c(2,j)+0.5D0*dyj-ymedi
2157 zj=c(3,j)+0.5D0*dzj-zmedi
2158 rij=xj*xj+yj*yj+zj*zj
2161 sss=sscale(rij/rpp(iteli,itelj))
2162 if (sss.gt.0.0d0) then
2167 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2168 if (j.eq.i+2) ev1=scal_el*ev1
2171 if (energy_dec) then
2172 write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
2173 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
2175 evdw1=evdw1+evdwij*sss
2177 C Calculate contributions to the Cartesian gradient.
2179 facvdw=-6*rrmij*(ev1+evdwij)*sss
2186 gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2187 gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2190 * Loop over residues i+1 thru j-1.
2194 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2202 C-----------------------------------------------------------------------------
2203 subroutine escp_long(evdw2,evdw2_14)
2205 C This subroutine calculates the excluded-volume interaction energy between
2206 C peptide-group centers and side chains and its gradient in virtual-bond and
2207 C side-chain vectors.
2209 implicit real*8 (a-h,o-z)
2210 include 'DIMENSIONS'
2211 include 'COMMON.GEO'
2212 include 'COMMON.VAR'
2213 include 'COMMON.LOCAL'
2214 include 'COMMON.CHAIN'
2215 include 'COMMON.DERIV'
2216 include 'COMMON.INTERACT'
2217 include 'COMMON.FFIELD'
2218 include 'COMMON.IOUNITS'
2219 include 'COMMON.CONTROL'
2223 cd print '(a)','Enter ESCP'
2224 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2225 do i=iatscp_s,iatscp_e
2227 xi=0.5D0*(c(1,i)+c(1,i+1))
2228 yi=0.5D0*(c(2,i)+c(2,i+1))
2229 zi=0.5D0*(c(3,i)+c(3,i+1))
2231 do iint=1,nscp_gr(i)
2233 do j=iscpstart(i,iint),iscpend(i,iint)
2235 C Uncomment following three lines for SC-p interactions
2239 C Uncomment following three lines for Ca-p interactions
2243 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2245 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2247 if (sss.lt.1.0d0) then
2250 e1=fac*fac*aad(itypj,iteli)
2251 e2=fac*bad(itypj,iteli)
2252 if (iabs(j-i) .le. 2) then
2255 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2258 evdw2=evdw2+evdwij*(1.0d0-sss)
2259 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2260 & 'evdw2',i,j,evdwij
2262 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2264 fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2269 cd write (iout,*) 'j<i'
2270 C Uncomment following three lines for SC-p interactions
2272 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2275 cd write (iout,*) 'j>i'
2278 C Uncomment following line for SC-p interactions
2279 c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
2283 gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
2287 cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
2288 cd write (iout,*) ggg(1),ggg(2),ggg(3)
2291 gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
2303 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2304 gradx_scp(j,i)=expon*gradx_scp(j,i)
2307 C******************************************************************************
2311 C To save time the factor EXPON has been extracted from ALL components
2312 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2315 C******************************************************************************
2318 C-----------------------------------------------------------------------------
2319 subroutine escp_short(evdw2,evdw2_14)
2321 C This subroutine calculates the excluded-volume interaction energy between
2322 C peptide-group centers and side chains and its gradient in virtual-bond and
2323 C side-chain vectors.
2325 implicit real*8 (a-h,o-z)
2326 include 'DIMENSIONS'
2327 include 'COMMON.GEO'
2328 include 'COMMON.VAR'
2329 include 'COMMON.LOCAL'
2330 include 'COMMON.CHAIN'
2331 include 'COMMON.DERIV'
2332 include 'COMMON.INTERACT'
2333 include 'COMMON.FFIELD'
2334 include 'COMMON.IOUNITS'
2335 include 'COMMON.CONTROL'
2339 cd print '(a)','Enter ESCP'
2340 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2341 do i=iatscp_s,iatscp_e
2343 xi=0.5D0*(c(1,i)+c(1,i+1))
2344 yi=0.5D0*(c(2,i)+c(2,i+1))
2345 zi=0.5D0*(c(3,i)+c(3,i+1))
2347 do iint=1,nscp_gr(i)
2349 do j=iscpstart(i,iint),iscpend(i,iint)
2351 C Uncomment following three lines for SC-p interactions
2355 C Uncomment following three lines for Ca-p interactions
2359 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2361 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2363 if (sss.gt.0.0d0) then
2366 e1=fac*fac*aad(itypj,iteli)
2367 e2=fac*bad(itypj,iteli)
2368 if (iabs(j-i) .le. 2) then
2371 evdw2_14=evdw2_14+(e1+e2)*sss
2374 evdw2=evdw2+evdwij*sss
2375 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2376 & 'evdw2',i,j,evdwij
2378 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2380 fac=-(evdwij+e1)*rrij*sss
2385 cd write (iout,*) 'j<i'
2386 C Uncomment following three lines for SC-p interactions
2388 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2391 cd write (iout,*) 'j>i'
2394 C Uncomment following line for SC-p interactions
2395 c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
2399 gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
2403 cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
2404 cd write (iout,*) ggg(1),ggg(2),ggg(3)
2407 gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
2419 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2420 gradx_scp(j,i)=expon*gradx_scp(j,i)
2423 C******************************************************************************
2427 C To save time the factor EXPON has been extracted from ALL components
2428 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2431 C******************************************************************************