1 C-----------------------------------------------------------------------
2 double precision function sscale(r)
3 double precision r,gamm
4 include "COMMON.SPLITELE"
5 if(r.lt.r_cut-rlamb) then
7 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
8 gamm=(r-(r_cut-rlamb))/rlamb
9 sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
15 C-----------------------------------------------------------------------
16 subroutine elj_long(evdw)
18 C This subroutine calculates the interaction energy of nonbonded side chains
19 C assuming the LJ potential of interaction.
21 implicit real*8 (a-h,o-z)
23 parameter (accur=1.0d-10)
26 include 'COMMON.LOCAL'
27 include 'COMMON.CHAIN'
28 include 'COMMON.DERIV'
29 include 'COMMON.INTERACT'
30 include 'COMMON.TORSION'
31 include 'COMMON.SBRIDGE'
32 include 'COMMON.NAMES'
33 include 'COMMON.IOUNITS'
34 include 'COMMON.CONTACTS'
36 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
45 C Calculate SC interaction energy.
48 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
49 cd & 'iend=',iend(i,iint)
50 do j=istart(i,iint),iend(i,iint)
56 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
57 if (sss.lt.1.0d0) then
59 eps0ij=eps(itypi,itypj)
61 e1=fac*fac*aa(itypi,itypj)
62 e2=fac*bb(itypi,itypj)
64 evdw=evdw+(1.0d0-sss)*evdwij
66 C Calculate the components of the gradient in DC and X
68 fac=-rrij*(e1+evdwij)*(1.0d0-sss)
73 gvdwx(k,i)=gvdwx(k,i)-gg(k)
74 gvdwx(k,j)=gvdwx(k,j)+gg(k)
75 gvdwc(k,i)=gvdwc(k,i)-gg(k)
76 gvdwc(k,j)=gvdwc(k,j)+gg(k)
84 gvdwc(j,i)=expon*gvdwc(j,i)
85 gvdwx(j,i)=expon*gvdwx(j,i)
88 C******************************************************************************
92 C To save time, the factor of EXPON has been extracted from ALL components
93 C of GVDWC and GRADX. Remember to multiply them by this factor before further
96 C******************************************************************************
99 C-----------------------------------------------------------------------
100 subroutine elj_short(evdw)
102 C This subroutine calculates the interaction energy of nonbonded side chains
103 C assuming the LJ potential of interaction.
105 implicit real*8 (a-h,o-z)
107 parameter (accur=1.0d-10)
110 include 'COMMON.LOCAL'
111 include 'COMMON.CHAIN'
112 include 'COMMON.DERIV'
113 include 'COMMON.INTERACT'
114 include 'COMMON.TORSION'
115 include 'COMMON.SBRIDGE'
116 include 'COMMON.NAMES'
117 include 'COMMON.IOUNITS'
118 include 'COMMON.CONTACTS'
120 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
131 C Calculate SC interaction energy.
134 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
135 cd & 'iend=',iend(i,iint)
136 do j=istart(i,iint),iend(i,iint)
141 C Change 12/1/95 to calculate four-body interactions
142 rij=xj*xj+yj*yj+zj*zj
143 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
144 if (sss.gt.0.0d0) then
146 eps0ij=eps(itypi,itypj)
148 e1=fac*fac*aa(itypi,itypj)
149 e2=fac*bb(itypi,itypj)
153 C Calculate the components of the gradient in DC and X
155 fac=-rrij*(e1+evdwij)*sss
160 gvdwx(k,i)=gvdwx(k,i)-gg(k)
161 gvdwx(k,j)=gvdwx(k,j)+gg(k)
162 gvdwc(k,i)=gvdwc(k,i)-gg(k)
163 gvdwc(k,j)=gvdwc(k,j)+gg(k)
171 gvdwc(j,i)=expon*gvdwc(j,i)
172 gvdwx(j,i)=expon*gvdwx(j,i)
175 C******************************************************************************
179 C To save time, the factor of EXPON has been extracted from ALL components
180 C of GVDWC and GRADX. Remember to multiply them by this factor before further
183 C******************************************************************************
186 C-----------------------------------------------------------------------------
187 subroutine eljk_long(evdw)
189 C This subroutine calculates the interaction energy of nonbonded side chains
190 C assuming the LJK potential of interaction.
192 implicit real*8 (a-h,o-z)
196 include 'COMMON.LOCAL'
197 include 'COMMON.CHAIN'
198 include 'COMMON.DERIV'
199 include 'COMMON.INTERACT'
200 include 'COMMON.IOUNITS'
201 include 'COMMON.NAMES'
204 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
213 C Calculate SC interaction energy.
216 do j=istart(i,iint),iend(i,iint)
221 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
223 e_augm=augm(itypi,itypj)*fac_augm
226 sss=sscale(rij/sigma(itypi,itypj))
227 if (sss.lt.1.0d0) then
228 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
229 fac=r_shift_inv**expon
230 e1=fac*fac*aa(itypi,itypj)
231 e2=fac*bb(itypi,itypj)
233 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
234 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
235 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
236 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
237 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
238 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
239 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
240 evdw=evdw+(1.0d0-sss)*evdwij
242 C Calculate the components of the gradient in DC and X
244 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
250 gvdwx(k,i)=gvdwx(k,i)-gg(k)
251 gvdwx(k,j)=gvdwx(k,j)+gg(k)
252 gvdwc(k,i)=gvdwc(k,i)-gg(k)
253 gvdwc(k,j)=gvdwc(k,j)+gg(k)
261 gvdwc(j,i)=expon*gvdwc(j,i)
262 gvdwx(j,i)=expon*gvdwx(j,i)
267 C-----------------------------------------------------------------------------
268 subroutine eljk_short(evdw)
270 C This subroutine calculates the interaction energy of nonbonded side chains
271 C assuming the LJK potential of interaction.
273 implicit real*8 (a-h,o-z)
277 include 'COMMON.LOCAL'
278 include 'COMMON.CHAIN'
279 include 'COMMON.DERIV'
280 include 'COMMON.INTERACT'
281 include 'COMMON.IOUNITS'
282 include 'COMMON.NAMES'
285 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
294 C Calculate SC interaction energy.
297 do j=istart(i,iint),iend(i,iint)
302 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
304 e_augm=augm(itypi,itypj)*fac_augm
307 sss=sscale(rij/sigma(itypi,itypj))
308 if (sss.gt.0.0d0) then
309 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
310 fac=r_shift_inv**expon
311 e1=fac*fac*aa(itypi,itypj)
312 e2=fac*bb(itypi,itypj)
314 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
315 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
316 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
317 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
318 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
319 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
320 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
323 C Calculate the components of the gradient in DC and X
325 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
331 gvdwx(k,i)=gvdwx(k,i)-gg(k)
332 gvdwx(k,j)=gvdwx(k,j)+gg(k)
333 gvdwc(k,i)=gvdwc(k,i)-gg(k)
334 gvdwc(k,j)=gvdwc(k,j)+gg(k)
342 gvdwc(j,i)=expon*gvdwc(j,i)
343 gvdwx(j,i)=expon*gvdwx(j,i)
348 C-----------------------------------------------------------------------------
349 subroutine ebp_long(evdw)
351 C This subroutine calculates the interaction energy of nonbonded side chains
352 C assuming the Berne-Pechukas potential of interaction.
354 implicit real*8 (a-h,o-z)
358 include 'COMMON.LOCAL'
359 include 'COMMON.CHAIN'
360 include 'COMMON.DERIV'
361 include 'COMMON.NAMES'
362 include 'COMMON.INTERACT'
363 include 'COMMON.IOUNITS'
364 include 'COMMON.CALC'
366 c double precision rrsave(maxdim)
369 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
371 c if (icall.eq.0) then
383 dxi=dc_norm(1,nres+i)
384 dyi=dc_norm(2,nres+i)
385 dzi=dc_norm(3,nres+i)
386 c dsci_inv=dsc_inv(itypi)
387 dsci_inv=vbld_inv(i+nres)
389 C Calculate SC interaction energy.
392 do j=istart(i,iint),iend(i,iint)
395 c dscj_inv=dsc_inv(itypj)
396 dscj_inv=vbld_inv(j+nres)
397 chi1=chi(itypi,itypj)
398 chi2=chi(itypj,itypi)
405 alf12=0.5D0*(alf1+alf2)
409 dxj=dc_norm(1,nres+j)
410 dyj=dc_norm(2,nres+j)
411 dzj=dc_norm(3,nres+j)
412 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
414 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
416 if (sss.lt.1.0d0) then
418 C Calculate the angle-dependent terms of energy & contributions to derivatives.
420 C Calculate whole angle-dependent part of epsilon and contributions
422 fac=(rrij*sigsq)**expon2
423 e1=fac*fac*aa(itypi,itypj)
424 e2=fac*bb(itypi,itypj)
425 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
426 eps2der=evdwij*eps3rt
427 eps3der=evdwij*eps2rt
428 evdwij=evdwij*eps2rt*eps3rt
429 evdw=evdw+evdwij*(1.0d0-sss)
431 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
432 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
433 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
434 cd & restyp(itypi),i,restyp(itypj),j,
435 cd & epsi,sigm,chi1,chi2,chip1,chip2,
436 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
437 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
440 C Calculate gradient components.
441 e1=e1*eps1*eps2rt**2*eps3rt**2
442 fac=-expon*(e1+evdwij)
445 C Calculate radial part of the gradient
449 C Calculate the angular part of the gradient and sum add the contributions
450 C to the appropriate components of the Cartesian gradient.
451 call sc_grad_scale(1.0d0-sss)
459 C-----------------------------------------------------------------------------
460 subroutine ebp_short(evdw)
462 C This subroutine calculates the interaction energy of nonbonded side chains
463 C assuming the Berne-Pechukas potential of interaction.
465 implicit real*8 (a-h,o-z)
469 include 'COMMON.LOCAL'
470 include 'COMMON.CHAIN'
471 include 'COMMON.DERIV'
472 include 'COMMON.NAMES'
473 include 'COMMON.INTERACT'
474 include 'COMMON.IOUNITS'
475 include 'COMMON.CALC'
477 c double precision rrsave(maxdim)
480 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
482 c if (icall.eq.0) then
494 dxi=dc_norm(1,nres+i)
495 dyi=dc_norm(2,nres+i)
496 dzi=dc_norm(3,nres+i)
497 c dsci_inv=dsc_inv(itypi)
498 dsci_inv=vbld_inv(i+nres)
500 C Calculate SC interaction energy.
503 do j=istart(i,iint),iend(i,iint)
506 c dscj_inv=dsc_inv(itypj)
507 dscj_inv=vbld_inv(j+nres)
508 chi1=chi(itypi,itypj)
509 chi2=chi(itypj,itypi)
516 alf12=0.5D0*(alf1+alf2)
520 dxj=dc_norm(1,nres+j)
521 dyj=dc_norm(2,nres+j)
522 dzj=dc_norm(3,nres+j)
523 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
525 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
527 if (sss.gt.0.0d0) then
529 C Calculate the angle-dependent terms of energy & contributions to derivatives.
531 C Calculate whole angle-dependent part of epsilon and contributions
533 fac=(rrij*sigsq)**expon2
534 e1=fac*fac*aa(itypi,itypj)
535 e2=fac*bb(itypi,itypj)
536 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
537 eps2der=evdwij*eps3rt
538 eps3der=evdwij*eps2rt
539 evdwij=evdwij*eps2rt*eps3rt
542 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
543 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
544 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
545 cd & restyp(itypi),i,restyp(itypj),j,
546 cd & epsi,sigm,chi1,chi2,chip1,chip2,
547 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
548 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
551 C Calculate gradient components.
552 e1=e1*eps1*eps2rt**2*eps3rt**2
553 fac=-expon*(e1+evdwij)
556 C Calculate radial part of the gradient
560 C Calculate the angular part of the gradient and sum add the contributions
561 C to the appropriate components of the Cartesian gradient.
562 call sc_grad_scale(sss)
570 C-----------------------------------------------------------------------------
571 subroutine egb_long(evdw,evdw_p,evdw_m)
573 C This subroutine calculates the interaction energy of nonbonded side chains
574 C assuming the Gay-Berne potential of interaction.
576 implicit real*8 (a-h,o-z)
580 include 'COMMON.LOCAL'
581 include 'COMMON.CHAIN'
582 include 'COMMON.DERIV'
583 include 'COMMON.NAMES'
584 include 'COMMON.INTERACT'
585 include 'COMMON.IOUNITS'
586 include 'COMMON.CALC'
587 include 'COMMON.CONTROL'
589 ccccc energy_dec=.false.
590 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
595 c if (icall.eq.0) lprn=.false.
603 dxi=dc_norm(1,nres+i)
604 dyi=dc_norm(2,nres+i)
605 dzi=dc_norm(3,nres+i)
606 c dsci_inv=dsc_inv(itypi)
607 dsci_inv=vbld_inv(i+nres)
608 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
609 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
611 C Calculate SC interaction energy.
614 do j=istart(i,iint),iend(i,iint)
617 c dscj_inv=dsc_inv(itypj)
618 dscj_inv=vbld_inv(j+nres)
619 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
620 c & 1.0d0/vbld(j+nres)
621 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
622 sig0ij=sigma(itypi,itypj)
623 chi1=chi(itypi,itypj)
624 chi2=chi(itypj,itypi)
631 alf12=0.5D0*(alf1+alf2)
635 dxj=dc_norm(1,nres+j)
636 dyj=dc_norm(2,nres+j)
637 dzj=dc_norm(3,nres+j)
638 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
640 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
642 if (sss.lt.1.0d0) then
644 C Calculate angle-dependent terms of energy and contributions to their
648 sig=sig0ij*dsqrt(sigsq)
649 rij_shift=1.0D0/rij-sig+sig0ij
650 c for diagnostics; uncomment
651 c rij_shift=1.2*sig0ij
652 C I hate to put IF's in the loops, but here don't have another choice!!!!
653 if (rij_shift.le.0.0D0) then
655 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
656 cd & restyp(itypi),i,restyp(itypj),j,
657 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
661 c---------------------------------------------------------------
662 rij_shift=1.0D0/rij_shift
664 e1=fac*fac*aa(itypi,itypj)
665 e2=fac*bb(itypi,itypj)
666 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
667 eps2der=evdwij*eps3rt
668 eps3der=evdwij*eps2rt
669 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
670 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
671 evdwij=evdwij*eps2rt*eps3rt
673 if (bb(itypi,itypj).gt.0) then
674 evdw_p=evdw_p+evdwij*(1.0d0-sss)
676 evdw_m=evdw_m+evdwij*(1.0d0-sss)
679 evdw=evdw+evdwij*(1.0d0-sss)
682 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
683 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
684 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
685 & restyp(itypi),i,restyp(itypj),j,
686 & epsi,sigm,chi1,chi2,chip1,chip2,
687 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
688 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
692 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
695 C Calculate gradient components.
696 e1=e1*eps1*eps2rt**2*eps3rt**2
697 fac=-expon*(e1+evdwij)*rij_shift
701 C Calculate the radial part of the gradient
705 C Calculate angular part of the gradient.
707 if (bb(itypi,itypj).gt.0) then
708 call sc_grad_scale_T(1.0d0-sss)
710 call sc_grad_scale(1.0d0-sss)
713 call sc_grad_scale(1.0d0-sss)
719 c write (iout,*) "Number of loop steps in EGB:",ind
720 cccc energy_dec=.false.
723 C-----------------------------------------------------------------------------
724 subroutine egb_short(evdw,evdw_p,evdw_m)
726 C This subroutine calculates the interaction energy of nonbonded side chains
727 C assuming the Gay-Berne potential of interaction.
729 implicit real*8 (a-h,o-z)
733 include 'COMMON.LOCAL'
734 include 'COMMON.CHAIN'
735 include 'COMMON.DERIV'
736 include 'COMMON.NAMES'
737 include 'COMMON.INTERACT'
738 include 'COMMON.IOUNITS'
739 include 'COMMON.CALC'
740 include 'COMMON.CONTROL'
745 ccccc energy_dec=.false.
746 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
749 c if (icall.eq.0) lprn=.false.
757 dxi=dc_norm(1,nres+i)
758 dyi=dc_norm(2,nres+i)
759 dzi=dc_norm(3,nres+i)
760 c dsci_inv=dsc_inv(itypi)
761 dsci_inv=vbld_inv(i+nres)
762 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
763 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
765 C Calculate SC interaction energy.
768 do j=istart(i,iint),iend(i,iint)
771 c dscj_inv=dsc_inv(itypj)
772 dscj_inv=vbld_inv(j+nres)
773 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
774 c & 1.0d0/vbld(j+nres)
775 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
776 sig0ij=sigma(itypi,itypj)
777 chi1=chi(itypi,itypj)
778 chi2=chi(itypj,itypi)
785 alf12=0.5D0*(alf1+alf2)
789 dxj=dc_norm(1,nres+j)
790 dyj=dc_norm(2,nres+j)
791 dzj=dc_norm(3,nres+j)
792 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
794 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
796 if (sss.gt.0.0d0) then
798 C Calculate angle-dependent terms of energy and contributions to their
802 sig=sig0ij*dsqrt(sigsq)
803 rij_shift=1.0D0/rij-sig+sig0ij
804 c for diagnostics; uncomment
805 c rij_shift=1.2*sig0ij
806 C I hate to put IF's in the loops, but here don't have another choice!!!!
807 if (rij_shift.le.0.0D0) then
809 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
810 cd & restyp(itypi),i,restyp(itypj),j,
811 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
815 c---------------------------------------------------------------
816 rij_shift=1.0D0/rij_shift
818 e1=fac*fac*aa(itypi,itypj)
819 e2=fac*bb(itypi,itypj)
820 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
821 eps2der=evdwij*eps3rt
822 eps3der=evdwij*eps2rt
823 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
824 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
825 evdwij=evdwij*eps2rt*eps3rt
827 if (bb(itypi,itypj).gt.0) then
828 evdw_p=evdw_p+evdwij*sss
830 evdw_m=evdw_m+evdwij*sss
836 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
837 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
838 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
839 & restyp(itypi),i,restyp(itypj),j,
840 & epsi,sigm,chi1,chi2,chip1,chip2,
841 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
842 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
846 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
849 C Calculate gradient components.
850 e1=e1*eps1*eps2rt**2*eps3rt**2
851 fac=-expon*(e1+evdwij)*rij_shift
855 C Calculate the radial part of the gradient
859 C Calculate angular part of the gradient.
861 if (bb(itypi,itypj).gt.0) then
862 call sc_grad_scale_T(sss)
864 call sc_grad_scale(sss)
867 call sc_grad_scale(sss)
873 c write (iout,*) "Number of loop steps in EGB:",ind
874 cccc energy_dec=.false.
877 C-----------------------------------------------------------------------------
878 subroutine egbv_long(evdw)
880 C This subroutine calculates the interaction energy of nonbonded side chains
881 C assuming the Gay-Berne-Vorobjev potential of interaction.
883 implicit real*8 (a-h,o-z)
887 include 'COMMON.LOCAL'
888 include 'COMMON.CHAIN'
889 include 'COMMON.DERIV'
890 include 'COMMON.NAMES'
891 include 'COMMON.INTERACT'
892 include 'COMMON.IOUNITS'
893 include 'COMMON.CALC'
897 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
900 c if (icall.eq.0) lprn=.true.
908 dxi=dc_norm(1,nres+i)
909 dyi=dc_norm(2,nres+i)
910 dzi=dc_norm(3,nres+i)
911 c dsci_inv=dsc_inv(itypi)
912 dsci_inv=vbld_inv(i+nres)
914 C Calculate SC interaction energy.
917 do j=istart(i,iint),iend(i,iint)
920 c dscj_inv=dsc_inv(itypj)
921 dscj_inv=vbld_inv(j+nres)
922 sig0ij=sigma(itypi,itypj)
924 chi1=chi(itypi,itypj)
925 chi2=chi(itypj,itypi)
932 alf12=0.5D0*(alf1+alf2)
936 dxj=dc_norm(1,nres+j)
937 dyj=dc_norm(2,nres+j)
938 dzj=dc_norm(3,nres+j)
939 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
942 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
944 if (sss.lt.1.0d0) then
946 C Calculate angle-dependent terms of energy and contributions to their
950 sig=sig0ij*dsqrt(sigsq)
951 rij_shift=1.0D0/rij-sig+r0ij
952 C I hate to put IF's in the loops, but here don't have another choice!!!!
953 if (rij_shift.le.0.0D0) then
958 c---------------------------------------------------------------
959 rij_shift=1.0D0/rij_shift
961 e1=fac*fac*aa(itypi,itypj)
962 e2=fac*bb(itypi,itypj)
963 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
964 eps2der=evdwij*eps3rt
965 eps3der=evdwij*eps2rt
967 e_augm=augm(itypi,itypj)*fac_augm
968 evdwij=evdwij*eps2rt*eps3rt
969 evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
971 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
972 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
973 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
974 & restyp(itypi),i,restyp(itypj),j,
975 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
976 & chi1,chi2,chip1,chip2,
977 & eps1,eps2rt**2,eps3rt**2,
978 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
981 C Calculate gradient components.
982 e1=e1*eps1*eps2rt**2*eps3rt**2
983 fac=-expon*(e1+evdwij)*rij_shift
985 fac=rij*fac-2*expon*rrij*e_augm
986 C Calculate the radial part of the gradient
990 C Calculate angular part of the gradient.
991 call sc_grad_scale(1.0d0-sss)
997 C-----------------------------------------------------------------------------
998 subroutine egbv_short(evdw)
1000 C This subroutine calculates the interaction energy of nonbonded side chains
1001 C assuming the Gay-Berne-Vorobjev potential of interaction.
1003 implicit real*8 (a-h,o-z)
1004 include 'DIMENSIONS'
1005 include 'COMMON.GEO'
1006 include 'COMMON.VAR'
1007 include 'COMMON.LOCAL'
1008 include 'COMMON.CHAIN'
1009 include 'COMMON.DERIV'
1010 include 'COMMON.NAMES'
1011 include 'COMMON.INTERACT'
1012 include 'COMMON.IOUNITS'
1013 include 'COMMON.CALC'
1014 common /srutu/ icall
1017 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1020 c if (icall.eq.0) lprn=.true.
1022 do i=iatsc_s,iatsc_e
1028 dxi=dc_norm(1,nres+i)
1029 dyi=dc_norm(2,nres+i)
1030 dzi=dc_norm(3,nres+i)
1031 c dsci_inv=dsc_inv(itypi)
1032 dsci_inv=vbld_inv(i+nres)
1034 C Calculate SC interaction energy.
1036 do iint=1,nint_gr(i)
1037 do j=istart(i,iint),iend(i,iint)
1040 c dscj_inv=dsc_inv(itypj)
1041 dscj_inv=vbld_inv(j+nres)
1042 sig0ij=sigma(itypi,itypj)
1043 r0ij=r0(itypi,itypj)
1044 chi1=chi(itypi,itypj)
1045 chi2=chi(itypj,itypi)
1052 alf12=0.5D0*(alf1+alf2)
1056 dxj=dc_norm(1,nres+j)
1057 dyj=dc_norm(2,nres+j)
1058 dzj=dc_norm(3,nres+j)
1059 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1062 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1064 if (sss.gt.0.0d0) then
1066 C Calculate angle-dependent terms of energy and contributions to their
1070 sig=sig0ij*dsqrt(sigsq)
1071 rij_shift=1.0D0/rij-sig+r0ij
1072 C I hate to put IF's in the loops, but here don't have another choice!!!!
1073 if (rij_shift.le.0.0D0) then
1078 c---------------------------------------------------------------
1079 rij_shift=1.0D0/rij_shift
1080 fac=rij_shift**expon
1081 e1=fac*fac*aa(itypi,itypj)
1082 e2=fac*bb(itypi,itypj)
1083 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1084 eps2der=evdwij*eps3rt
1085 eps3der=evdwij*eps2rt
1086 fac_augm=rrij**expon
1087 e_augm=augm(itypi,itypj)*fac_augm
1088 evdwij=evdwij*eps2rt*eps3rt
1089 evdw=evdw+(evdwij+e_augm)*sss
1091 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1092 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1093 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1094 & restyp(itypi),i,restyp(itypj),j,
1095 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1096 & chi1,chi2,chip1,chip2,
1097 & eps1,eps2rt**2,eps3rt**2,
1098 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1101 C Calculate gradient components.
1102 e1=e1*eps1*eps2rt**2*eps3rt**2
1103 fac=-expon*(e1+evdwij)*rij_shift
1105 fac=rij*fac-2*expon*rrij*e_augm
1106 C Calculate the radial part of the gradient
1110 C Calculate angular part of the gradient.
1111 call sc_grad_scale(sss)
1117 C----------------------------------------------------------------------------
1118 subroutine sc_grad_scale(scalfac)
1119 implicit real*8 (a-h,o-z)
1120 include 'DIMENSIONS'
1121 include 'COMMON.CHAIN'
1122 include 'COMMON.DERIV'
1123 include 'COMMON.CALC'
1124 include 'COMMON.IOUNITS'
1125 double precision dcosom1(3),dcosom2(3)
1126 double precision scalfac
1127 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1128 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1129 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1130 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1134 c eom12=evdwij*eps1_om12
1136 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1137 c & " sigder",sigder
1138 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1139 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1141 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1142 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1145 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1147 c write (iout,*) "gg",(gg(k),k=1,3)
1149 gvdwx(k,i)=gvdwx(k,i)-gg(k)
1150 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1151 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1152 gvdwx(k,j)=gvdwx(k,j)+gg(k)
1153 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1154 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1155 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1156 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1157 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1158 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1161 C Calculate the components of the gradient in DC and X
1164 gvdwc(l,i)=gvdwc(l,i)-gg(l)
1165 gvdwc(l,j)=gvdwc(l,j)+gg(l)
1169 C----------------------------------------------------------------------------
1170 subroutine sc_grad_scale_T(scalfac)
1171 implicit real*8 (a-h,o-z)
1172 include 'DIMENSIONS'
1173 include 'COMMON.CHAIN'
1174 include 'COMMON.DERIV'
1175 include 'COMMON.CALC'
1176 include 'COMMON.IOUNITS'
1177 double precision dcosom1(3),dcosom2(3)
1178 double precision scalfac
1179 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1180 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1181 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1182 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1186 c eom12=evdwij*eps1_om12
1188 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1189 c & " sigder",sigder
1190 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1191 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1193 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1194 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1197 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1199 c write (iout,*) "gg",(gg(k),k=1,3)
1201 gvdwxT(k,i)=gvdwxT(k,i)-gg(k)
1202 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1203 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1204 gvdwxT(k,j)=gvdwxT(k,j)+gg(k)
1205 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1206 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1207 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1208 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1209 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1210 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1213 C Calculate the components of the gradient in DC and X
1216 gvdwcT(l,i)=gvdwcT(l,i)-gg(l)
1217 gvdwcT(l,j)=gvdwcT(l,j)+gg(l)
1222 C--------------------------------------------------------------------------
1223 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1225 C This subroutine calculates the average interaction energy and its gradient
1226 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1227 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1228 C The potential depends both on the distance of peptide-group centers and on
1229 C the orientation of the CA-CA virtual bonds.
1231 implicit real*8 (a-h,o-z)
1235 include 'DIMENSIONS'
1236 include 'COMMON.CONTROL'
1237 include 'COMMON.SETUP'
1238 include 'COMMON.IOUNITS'
1239 include 'COMMON.GEO'
1240 include 'COMMON.VAR'
1241 include 'COMMON.LOCAL'
1242 include 'COMMON.CHAIN'
1243 include 'COMMON.DERIV'
1244 include 'COMMON.INTERACT'
1245 include 'COMMON.CONTACTS'
1246 include 'COMMON.TORSION'
1247 include 'COMMON.VECTORS'
1248 include 'COMMON.FFIELD'
1249 include 'COMMON.TIME1'
1250 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1251 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1252 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1253 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1254 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1255 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1257 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1259 double precision scal_el /1.0d0/
1261 double precision scal_el /0.5d0/
1264 C 13-go grudnia roku pamietnego...
1265 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1266 & 0.0d0,1.0d0,0.0d0,
1267 & 0.0d0,0.0d0,1.0d0/
1268 cd write(iout,*) 'In EELEC'
1270 cd write(iout,*) 'Type',i
1271 cd write(iout,*) 'B1',B1(:,i)
1272 cd write(iout,*) 'B2',B2(:,i)
1273 cd write(iout,*) 'CC',CC(:,:,i)
1274 cd write(iout,*) 'DD',DD(:,:,i)
1275 cd write(iout,*) 'EE',EE(:,:,i)
1277 cd call check_vecgrad
1279 if (icheckgrad.eq.1) then
1281 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1283 dc_norm(k,i)=dc(k,i)*fac
1285 c write (iout,*) 'i',i,' fac',fac
1288 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1289 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1290 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1291 c call vec_and_deriv
1297 time_mat=time_mat+MPI_Wtime()-time01
1301 cd write (iout,*) 'i=',i
1303 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1306 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1307 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1320 cd print '(a)','Enter EELEC'
1321 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1323 gel_loc_loc(i)=0.0d0
1328 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1330 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1332 do i=iturn3_start,iturn3_end
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 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1344 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1345 num_cont_hb(i)=num_conti
1347 do i=iturn4_start,iturn4_end
1351 dx_normi=dc_norm(1,i)
1352 dy_normi=dc_norm(2,i)
1353 dz_normi=dc_norm(3,i)
1354 xmedi=c(1,i)+0.5d0*dxi
1355 ymedi=c(2,i)+0.5d0*dyi
1356 zmedi=c(3,i)+0.5d0*dzi
1357 num_conti=num_cont_hb(i)
1358 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1359 if (wturn4.gt.0.0d0) call eturn4(i,eello_turn4)
1360 num_cont_hb(i)=num_conti
1363 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1365 do i=iatel_s,iatel_e
1369 dx_normi=dc_norm(1,i)
1370 dy_normi=dc_norm(2,i)
1371 dz_normi=dc_norm(3,i)
1372 xmedi=c(1,i)+0.5d0*dxi
1373 ymedi=c(2,i)+0.5d0*dyi
1374 zmedi=c(3,i)+0.5d0*dzi
1375 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1376 num_conti=num_cont_hb(i)
1377 do j=ielstart(i),ielend(i)
1378 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1380 num_cont_hb(i)=num_conti
1382 c write (iout,*) "Number of loop steps in EELEC:",ind
1384 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1385 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1387 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1388 ccc eel_loc=eel_loc+eello_turn3
1389 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1392 C-------------------------------------------------------------------------------
1393 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1394 implicit real*8 (a-h,o-z)
1395 include 'DIMENSIONS'
1399 include 'COMMON.CONTROL'
1400 include 'COMMON.IOUNITS'
1401 include 'COMMON.GEO'
1402 include 'COMMON.VAR'
1403 include 'COMMON.LOCAL'
1404 include 'COMMON.CHAIN'
1405 include 'COMMON.DERIV'
1406 include 'COMMON.INTERACT'
1407 include 'COMMON.CONTACTS'
1408 include 'COMMON.TORSION'
1409 include 'COMMON.VECTORS'
1410 include 'COMMON.FFIELD'
1411 include 'COMMON.TIME1'
1412 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1413 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1414 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1415 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1416 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1417 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1419 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1421 double precision scal_el /1.0d0/
1423 double precision scal_el /0.5d0/
1426 C 13-go grudnia roku pamietnego...
1427 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1428 & 0.0d0,1.0d0,0.0d0,
1429 & 0.0d0,0.0d0,1.0d0/
1430 c time00=MPI_Wtime()
1431 cd write (iout,*) "eelecij",i,j
1435 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1436 aaa=app(iteli,itelj)
1437 bbb=bpp(iteli,itelj)
1438 ael6i=ael6(iteli,itelj)
1439 ael3i=ael3(iteli,itelj)
1443 dx_normj=dc_norm(1,j)
1444 dy_normj=dc_norm(2,j)
1445 dz_normj=dc_norm(3,j)
1446 xj=c(1,j)+0.5D0*dxj-xmedi
1447 yj=c(2,j)+0.5D0*dyj-ymedi
1448 zj=c(3,j)+0.5D0*dzj-zmedi
1449 rij=xj*xj+yj*yj+zj*zj
1453 c For extracting the short-range part of Evdwpp
1454 sss=sscale(rij/rpp(iteli,itelj))
1458 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1459 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1460 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1461 fac=cosa-3.0D0*cosb*cosg
1463 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1464 if (j.eq.i+2) ev1=scal_el*ev1
1469 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1472 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1473 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1475 evdw1=evdw1+evdwij*(1.0d0-sss)
1476 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1477 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1478 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1479 cd & xmedi,ymedi,zmedi,xj,yj,zj
1481 if (energy_dec) then
1482 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1483 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1487 C Calculate contributions to the Cartesian gradient.
1490 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1491 facel=-3*rrmij*(el1+eesij)
1497 * Radial derivatives. First process both termini of the fragment (i,j)
1503 c ghalf=0.5D0*ggg(k)
1504 c gelc(k,i)=gelc(k,i)+ghalf
1505 c gelc(k,j)=gelc(k,j)+ghalf
1507 c 9/28/08 AL Gradient compotents will be summed only at the end
1509 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1510 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1513 * Loop over residues i+1 thru j-1.
1517 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1524 c ghalf=0.5D0*ggg(k)
1525 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1526 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1528 c 9/28/08 AL Gradient compotents will be summed only at the end
1530 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1531 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1534 * Loop over residues i+1 thru j-1.
1538 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1542 facvdw=ev1+evdwij*(1.0d0-sss)
1545 fac=-3*rrmij*(facvdw+facvdw+facel)
1550 * Radial derivatives. First process both termini of the fragment (i,j)
1556 c ghalf=0.5D0*ggg(k)
1557 c gelc(k,i)=gelc(k,i)+ghalf
1558 c gelc(k,j)=gelc(k,j)+ghalf
1560 c 9/28/08 AL Gradient compotents will be summed only at the end
1562 gelc_long(k,j)=gelc(k,j)+ggg(k)
1563 gelc_long(k,i)=gelc(k,i)-ggg(k)
1566 * Loop over residues i+1 thru j-1.
1570 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1573 c 9/28/08 AL Gradient compotents will be summed only at the end
1578 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1579 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1585 ecosa=2.0D0*fac3*fac1+fac4
1588 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1589 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1591 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1592 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1594 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1595 cd & (dcosg(k),k=1,3)
1597 ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1600 c ghalf=0.5D0*ggg(k)
1601 c gelc(k,i)=gelc(k,i)+ghalf
1602 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1603 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1604 c gelc(k,j)=gelc(k,j)+ghalf
1605 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1606 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1610 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1615 & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1616 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1618 & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1619 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1620 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1621 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1623 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1624 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1625 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1627 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1628 C energy of a peptide unit is assumed in the form of a second-order
1629 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1630 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1631 C are computed for EVERY pair of non-contiguous peptide groups.
1633 if (j.lt.nres-1) then
1644 muij(kkk)=mu(k,i)*mu(l,j)
1647 cd write (iout,*) 'EELEC: i',i,' j',j
1648 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1649 cd write(iout,*) 'muij',muij
1650 ury=scalar(uy(1,i),erij)
1651 urz=scalar(uz(1,i),erij)
1652 vry=scalar(uy(1,j),erij)
1653 vrz=scalar(uz(1,j),erij)
1654 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1655 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1656 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1657 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1658 fac=dsqrt(-ael6i)*r3ij
1663 cd write (iout,'(4i5,4f10.5)')
1664 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1665 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1666 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1667 cd & uy(:,j),uz(:,j)
1668 cd write (iout,'(4f10.5)')
1669 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1670 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1671 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
1672 cd write (iout,'(9f10.5/)')
1673 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1674 C Derivatives of the elements of A in virtual-bond vectors
1675 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1677 uryg(k,1)=scalar(erder(1,k),uy(1,i))
1678 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1679 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1680 urzg(k,1)=scalar(erder(1,k),uz(1,i))
1681 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1682 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1683 vryg(k,1)=scalar(erder(1,k),uy(1,j))
1684 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1685 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1686 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1687 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1688 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1690 C Compute radial contributions to the gradient
1708 C Add the contributions coming from er
1711 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1712 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1713 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1714 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1717 C Derivatives in DC(i)
1718 cgrad ghalf1=0.5d0*agg(k,1)
1719 cgrad ghalf2=0.5d0*agg(k,2)
1720 cgrad ghalf3=0.5d0*agg(k,3)
1721 cgrad ghalf4=0.5d0*agg(k,4)
1722 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1723 & -3.0d0*uryg(k,2)*vry)!+ghalf1
1724 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1725 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
1726 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1727 & -3.0d0*urzg(k,2)*vry)!+ghalf3
1728 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1729 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
1730 C Derivatives in DC(i+1)
1731 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
1732 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
1733 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
1734 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
1735 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
1736 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
1737 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
1738 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
1739 C Derivatives in DC(j)
1740 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
1741 & -3.0d0*vryg(k,2)*ury)!+ghalf1
1742 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
1743 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
1744 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
1745 & -3.0d0*vryg(k,2)*urz)!+ghalf3
1746 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
1747 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
1748 C Derivatives in DC(j+1) or DC(nres-1)
1749 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
1750 & -3.0d0*vryg(k,3)*ury)
1751 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
1752 & -3.0d0*vrzg(k,3)*ury)
1753 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
1754 & -3.0d0*vryg(k,3)*urz)
1755 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
1756 & -3.0d0*vrzg(k,3)*urz)
1757 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
1759 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
1772 aggi(k,l)=-aggi(k,l)
1773 aggi1(k,l)=-aggi1(k,l)
1774 aggj(k,l)=-aggj(k,l)
1775 aggj1(k,l)=-aggj1(k,l)
1778 if (j.lt.nres-1) then
1784 aggi(k,l)=-aggi(k,l)
1785 aggi1(k,l)=-aggi1(k,l)
1786 aggj(k,l)=-aggj(k,l)
1787 aggj1(k,l)=-aggj1(k,l)
1798 aggi(k,l)=-aggi(k,l)
1799 aggi1(k,l)=-aggi1(k,l)
1800 aggj(k,l)=-aggj(k,l)
1801 aggj1(k,l)=-aggj1(k,l)
1806 IF (wel_loc.gt.0.0d0) THEN
1807 C Contribution to the local-electrostatic energy coming from the i-j pair
1808 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
1810 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
1812 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1813 & 'eelloc',i,j,eel_loc_ij
1815 eel_loc=eel_loc+eel_loc_ij
1816 C Partial derivatives in virtual-bond dihedral angles gamma
1818 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
1819 & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1820 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1821 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
1822 & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1823 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1824 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
1826 ggg(l)=agg(l,1)*muij(1)+
1827 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1828 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
1829 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
1830 cgrad ghalf=0.5d0*ggg(l)
1831 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
1832 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
1836 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
1839 C Remaining derivatives of eello
1841 gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1842 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1843 gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1844 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1845 gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1846 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1847 gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1848 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1851 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
1852 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
1853 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
1854 & .and. num_conti.le.maxconts) then
1855 c write (iout,*) i,j," entered corr"
1857 C Calculate the contact function. The ith column of the array JCONT will
1858 C contain the numbers of atoms that make contacts with the atom I (of numbers
1859 C greater than I). The arrays FACONT and GACONT will contain the values of
1860 C the contact function and its derivative.
1861 c r0ij=1.02D0*rpp(iteli,itelj)
1862 c r0ij=1.11D0*rpp(iteli,itelj)
1863 r0ij=2.20D0*rpp(iteli,itelj)
1864 c r0ij=1.55D0*rpp(iteli,itelj)
1865 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
1866 if (fcont.gt.0.0D0) then
1867 num_conti=num_conti+1
1868 if (num_conti.gt.maxconts) then
1869 write (iout,*) 'WARNING - max. # of contacts exceeded;',
1870 & ' will skip next contacts for this conf.'
1872 jcont_hb(num_conti,i)=j
1873 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
1874 cd & " jcont_hb",jcont_hb(num_conti,i)
1875 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
1876 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1877 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
1879 d_cont(num_conti,i)=rij
1880 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
1881 C --- Electrostatic-interaction matrix ---
1882 a_chuj(1,1,num_conti,i)=a22
1883 a_chuj(1,2,num_conti,i)=a23
1884 a_chuj(2,1,num_conti,i)=a32
1885 a_chuj(2,2,num_conti,i)=a33
1886 C --- Gradient of rij
1888 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
1895 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
1896 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
1897 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
1898 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
1899 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
1904 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
1905 C Calculate contact energies
1907 wij=cosa-3.0D0*cosb*cosg
1910 c fac3=dsqrt(-ael6i)/r0ij**3
1911 fac3=dsqrt(-ael6i)*r3ij
1912 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
1913 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
1914 if (ees0tmp.gt.0) then
1915 ees0pij=dsqrt(ees0tmp)
1919 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
1920 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
1921 if (ees0tmp.gt.0) then
1922 ees0mij=dsqrt(ees0tmp)
1927 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
1928 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
1929 C Diagnostics. Comment out or remove after debugging!
1930 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
1931 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
1932 c ees0m(num_conti,i)=0.0D0
1934 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
1935 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
1936 C Angular derivatives of the contact function
1937 ees0pij1=fac3/ees0pij
1938 ees0mij1=fac3/ees0mij
1939 fac3p=-3.0D0*fac3*rrmij
1940 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
1941 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
1943 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
1944 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
1945 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
1946 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
1947 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
1948 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
1949 ecosap=ecosa1+ecosa2
1950 ecosbp=ecosb1+ecosb2
1951 ecosgp=ecosg1+ecosg2
1952 ecosam=ecosa1-ecosa2
1953 ecosbm=ecosb1-ecosb2
1954 ecosgm=ecosg1-ecosg2
1963 facont_hb(num_conti,i)=fcont
1964 fprimcont=fprimcont/rij
1965 cd facont_hb(num_conti,i)=1.0D0
1966 C Following line is for diagnostics.
1969 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1970 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1973 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
1974 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
1976 gggp(1)=gggp(1)+ees0pijp*xj
1977 gggp(2)=gggp(2)+ees0pijp*yj
1978 gggp(3)=gggp(3)+ees0pijp*zj
1979 gggm(1)=gggm(1)+ees0mijp*xj
1980 gggm(2)=gggm(2)+ees0mijp*yj
1981 gggm(3)=gggm(3)+ees0mijp*zj
1982 C Derivatives due to the contact function
1983 gacont_hbr(1,num_conti,i)=fprimcont*xj
1984 gacont_hbr(2,num_conti,i)=fprimcont*yj
1985 gacont_hbr(3,num_conti,i)=fprimcont*zj
1988 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
1989 c following the change of gradient-summation algorithm.
1991 cgrad ghalfp=0.5D0*gggp(k)
1992 cgrad ghalfm=0.5D0*gggm(k)
1993 gacontp_hb1(k,num_conti,i)=!ghalfp
1994 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
1995 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1996 gacontp_hb2(k,num_conti,i)=!ghalfp
1997 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
1998 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1999 gacontp_hb3(k,num_conti,i)=gggp(k)
2000 gacontm_hb1(k,num_conti,i)=!ghalfm
2001 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2002 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2003 gacontm_hb2(k,num_conti,i)=!ghalfm
2004 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2005 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2006 gacontm_hb3(k,num_conti,i)=gggm(k)
2009 endif ! num_conti.le.maxconts
2012 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2015 ghalf=0.5d0*agg(l,k)
2016 aggi(l,k)=aggi(l,k)+ghalf
2017 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2018 aggj(l,k)=aggj(l,k)+ghalf
2021 if (j.eq.nres-1 .and. i.lt.j-2) then
2024 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2029 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2032 C-----------------------------------------------------------------------
2033 subroutine evdwpp_short(evdw1)
2037 implicit real*8 (a-h,o-z)
2038 include 'DIMENSIONS'
2039 include 'COMMON.CONTROL'
2040 include 'COMMON.IOUNITS'
2041 include 'COMMON.GEO'
2042 include 'COMMON.VAR'
2043 include 'COMMON.LOCAL'
2044 include 'COMMON.CHAIN'
2045 include 'COMMON.DERIV'
2046 include 'COMMON.INTERACT'
2047 include 'COMMON.CONTACTS'
2048 include 'COMMON.TORSION'
2049 include 'COMMON.VECTORS'
2050 include 'COMMON.FFIELD'
2052 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2054 double precision scal_el /1.0d0/
2056 double precision scal_el /0.5d0/
2059 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2060 c & " iatel_e_vdw",iatel_e_vdw
2062 do i=iatel_s_vdw,iatel_e_vdw
2066 dx_normi=dc_norm(1,i)
2067 dy_normi=dc_norm(2,i)
2068 dz_normi=dc_norm(3,i)
2069 xmedi=c(1,i)+0.5d0*dxi
2070 ymedi=c(2,i)+0.5d0*dyi
2071 zmedi=c(3,i)+0.5d0*dzi
2073 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2074 c & ' ielend',ielend_vdw(i)
2076 do j=ielstart_vdw(i),ielend_vdw(i)
2080 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2081 aaa=app(iteli,itelj)
2082 bbb=bpp(iteli,itelj)
2086 dx_normj=dc_norm(1,j)
2087 dy_normj=dc_norm(2,j)
2088 dz_normj=dc_norm(3,j)
2089 xj=c(1,j)+0.5D0*dxj-xmedi
2090 yj=c(2,j)+0.5D0*dyj-ymedi
2091 zj=c(3,j)+0.5D0*dzj-zmedi
2092 rij=xj*xj+yj*yj+zj*zj
2095 sss=sscale(rij/rpp(iteli,itelj))
2096 if (sss.gt.0.0d0) then
2101 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2102 if (j.eq.i+2) ev1=scal_el*ev1
2105 if (energy_dec) then
2106 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2108 evdw1=evdw1+evdwij*sss
2110 C Calculate contributions to the Cartesian gradient.
2112 facvdw=-6*rrmij*(ev1+evdwij)*sss
2117 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2118 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2125 C-----------------------------------------------------------------------------
2126 subroutine escp_long(evdw2,evdw2_14)
2128 C This subroutine calculates the excluded-volume interaction energy between
2129 C peptide-group centers and side chains and its gradient in virtual-bond and
2130 C side-chain vectors.
2132 implicit real*8 (a-h,o-z)
2133 include 'DIMENSIONS'
2134 include 'COMMON.GEO'
2135 include 'COMMON.VAR'
2136 include 'COMMON.LOCAL'
2137 include 'COMMON.CHAIN'
2138 include 'COMMON.DERIV'
2139 include 'COMMON.INTERACT'
2140 include 'COMMON.FFIELD'
2141 include 'COMMON.IOUNITS'
2142 include 'COMMON.CONTROL'
2146 cd print '(a)','Enter ESCP'
2147 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2148 do i=iatscp_s,iatscp_e
2150 xi=0.5D0*(c(1,i)+c(1,i+1))
2151 yi=0.5D0*(c(2,i)+c(2,i+1))
2152 zi=0.5D0*(c(3,i)+c(3,i+1))
2154 do iint=1,nscp_gr(i)
2156 do j=iscpstart(i,iint),iscpend(i,iint)
2158 C Uncomment following three lines for SC-p interactions
2162 C Uncomment following three lines for Ca-p interactions
2166 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2168 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2170 if (sss.lt.1.0d0) then
2173 e1=fac*fac*aad(itypj,iteli)
2174 e2=fac*bad(itypj,iteli)
2175 if (iabs(j-i) .le. 2) then
2178 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2181 evdw2=evdw2+evdwij*(1.0d0-sss)
2182 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2183 & 'evdw2',i,j,evdwij
2185 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2187 fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2191 C Uncomment following three lines for SC-p interactions
2193 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2195 C Uncomment following line for SC-p interactions
2196 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2198 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2199 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2208 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2209 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2210 gradx_scp(j,i)=expon*gradx_scp(j,i)
2213 C******************************************************************************
2217 C To save time the factor EXPON has been extracted from ALL components
2218 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2221 C******************************************************************************
2224 C-----------------------------------------------------------------------------
2225 subroutine escp_short(evdw2,evdw2_14)
2227 C This subroutine calculates the excluded-volume interaction energy between
2228 C peptide-group centers and side chains and its gradient in virtual-bond and
2229 C side-chain vectors.
2231 implicit real*8 (a-h,o-z)
2232 include 'DIMENSIONS'
2233 include 'COMMON.GEO'
2234 include 'COMMON.VAR'
2235 include 'COMMON.LOCAL'
2236 include 'COMMON.CHAIN'
2237 include 'COMMON.DERIV'
2238 include 'COMMON.INTERACT'
2239 include 'COMMON.FFIELD'
2240 include 'COMMON.IOUNITS'
2241 include 'COMMON.CONTROL'
2245 cd print '(a)','Enter ESCP'
2246 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2247 do i=iatscp_s,iatscp_e
2249 xi=0.5D0*(c(1,i)+c(1,i+1))
2250 yi=0.5D0*(c(2,i)+c(2,i+1))
2251 zi=0.5D0*(c(3,i)+c(3,i+1))
2253 do iint=1,nscp_gr(i)
2255 do j=iscpstart(i,iint),iscpend(i,iint)
2257 C Uncomment following three lines for SC-p interactions
2261 C Uncomment following three lines for Ca-p interactions
2265 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2267 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2269 if (sss.gt.0.0d0) then
2272 e1=fac*fac*aad(itypj,iteli)
2273 e2=fac*bad(itypj,iteli)
2274 if (iabs(j-i) .le. 2) then
2277 evdw2_14=evdw2_14+(e1+e2)*sss
2280 evdw2=evdw2+evdwij*sss
2281 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2282 & 'evdw2',i,j,evdwij
2284 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2286 fac=-(evdwij+e1)*rrij*sss
2290 C Uncomment following three lines for SC-p interactions
2292 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2294 C Uncomment following line for SC-p interactions
2295 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2297 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2298 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2307 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2308 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2309 gradx_scp(j,i)=expon*gradx_scp(j,i)
2312 C******************************************************************************
2316 C To save time the factor EXPON has been extracted from ALL components
2317 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2320 C******************************************************************************