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
40 if (itypi.eq.ntyp1) cycle
46 C Calculate SC interaction energy.
49 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
50 cd & 'iend=',iend(i,iint)
51 do j=istart(i,iint),iend(i,iint)
53 if (itypj.eq.ntyp1) cycle
58 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
59 if (sss.lt.1.0d0) then
61 eps0ij=eps(itypi,itypj)
63 e1=fac*fac*aa(itypi,itypj)
64 e2=fac*bb(itypi,itypj)
66 evdw=evdw+(1.0d0-sss)*evdwij
68 C Calculate the components of the gradient in DC and X
70 fac=-rrij*(e1+evdwij)*(1.0d0-sss)
75 gvdwx(k,i)=gvdwx(k,i)-gg(k)
76 gvdwx(k,j)=gvdwx(k,j)+gg(k)
77 gvdwc(k,i)=gvdwc(k,i)-gg(k)
78 gvdwc(k,j)=gvdwc(k,j)+gg(k)
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
126 if (itypi.eq.ntyp1) cycle
134 C Calculate SC interaction energy.
137 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
138 cd & 'iend=',iend(i,iint)
139 do j=istart(i,iint),iend(i,iint)
141 if (itypj.eq.ntyp1) cycle
145 C Change 12/1/95 to calculate four-body interactions
146 rij=xj*xj+yj*yj+zj*zj
147 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
148 if (sss.gt.0.0d0) then
150 eps0ij=eps(itypi,itypj)
152 e1=fac*fac*aa(itypi,itypj)
153 e2=fac*bb(itypi,itypj)
157 C Calculate the components of the gradient in DC and X
159 fac=-rrij*(e1+evdwij)*sss
164 gvdwx(k,i)=gvdwx(k,i)-gg(k)
165 gvdwx(k,j)=gvdwx(k,j)+gg(k)
166 gvdwc(k,i)=gvdwc(k,i)-gg(k)
167 gvdwc(k,j)=gvdwc(k,j)+gg(k)
175 gvdwc(j,i)=expon*gvdwc(j,i)
176 gvdwx(j,i)=expon*gvdwx(j,i)
179 C******************************************************************************
183 C To save time, the factor of EXPON has been extracted from ALL components
184 C of GVDWC and GRADX. Remember to multiply them by this factor before further
187 C******************************************************************************
190 C-----------------------------------------------------------------------------
191 subroutine eljk_long(evdw)
193 C This subroutine calculates the interaction energy of nonbonded side chains
194 C assuming the LJK potential of interaction.
196 implicit real*8 (a-h,o-z)
200 include 'COMMON.LOCAL'
201 include 'COMMON.CHAIN'
202 include 'COMMON.DERIV'
203 include 'COMMON.INTERACT'
204 include 'COMMON.IOUNITS'
205 include 'COMMON.NAMES'
208 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
212 if (itypi.eq.ntyp1) cycle
218 C Calculate SC interaction energy.
221 do j=istart(i,iint),iend(i,iint)
223 if (itypj.eq.ntyp1) cycle
227 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
229 e_augm=augm(itypi,itypj)*fac_augm
232 sss=sscale(rij/sigma(itypi,itypj))
233 if (sss.lt.1.0d0) then
234 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
235 fac=r_shift_inv**expon
236 e1=fac*fac*aa(itypi,itypj)
237 e2=fac*bb(itypi,itypj)
239 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
240 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
241 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
242 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
243 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
244 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
245 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
246 evdw=evdw+(1.0d0-sss)*evdwij
248 C Calculate the components of the gradient in DC and X
250 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
256 gvdwx(k,i)=gvdwx(k,i)-gg(k)
257 gvdwx(k,j)=gvdwx(k,j)+gg(k)
258 gvdwc(k,i)=gvdwc(k,i)-gg(k)
259 gvdwc(k,j)=gvdwc(k,j)+gg(k)
267 gvdwc(j,i)=expon*gvdwc(j,i)
268 gvdwx(j,i)=expon*gvdwx(j,i)
273 C-----------------------------------------------------------------------------
274 subroutine eljk_short(evdw)
276 C This subroutine calculates the interaction energy of nonbonded side chains
277 C assuming the LJK potential of interaction.
279 implicit real*8 (a-h,o-z)
283 include 'COMMON.LOCAL'
284 include 'COMMON.CHAIN'
285 include 'COMMON.DERIV'
286 include 'COMMON.INTERACT'
287 include 'COMMON.IOUNITS'
288 include 'COMMON.NAMES'
291 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
295 if (itypi.eq.ntyp1) cycle
301 C Calculate SC interaction energy.
304 do j=istart(i,iint),iend(i,iint)
306 if (itypj.eq.ntyp1) cycle
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))
316 if (sss.gt.0.0d0) then
317 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
318 fac=r_shift_inv**expon
319 e1=fac*fac*aa(itypi,itypj)
320 e2=fac*bb(itypi,itypj)
322 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
323 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
324 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
325 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
326 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
327 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
328 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
331 C Calculate the components of the gradient in DC and X
333 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
339 gvdwx(k,i)=gvdwx(k,i)-gg(k)
340 gvdwx(k,j)=gvdwx(k,j)+gg(k)
341 gvdwc(k,i)=gvdwc(k,i)-gg(k)
342 gvdwc(k,j)=gvdwc(k,j)+gg(k)
350 gvdwc(j,i)=expon*gvdwc(j,i)
351 gvdwx(j,i)=expon*gvdwx(j,i)
356 C-----------------------------------------------------------------------------
357 subroutine ebp_long(evdw)
359 C This subroutine calculates the interaction energy of nonbonded side chains
360 C assuming the Berne-Pechukas potential of interaction.
362 implicit real*8 (a-h,o-z)
366 include 'COMMON.LOCAL'
367 include 'COMMON.CHAIN'
368 include 'COMMON.DERIV'
369 include 'COMMON.NAMES'
370 include 'COMMON.INTERACT'
371 include 'COMMON.IOUNITS'
372 include 'COMMON.CALC'
374 c double precision rrsave(maxdim)
377 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
379 c if (icall.eq.0) then
387 if (itypi.eq.ntyp1) cycle
392 dxi=dc_norm(1,nres+i)
393 dyi=dc_norm(2,nres+i)
394 dzi=dc_norm(3,nres+i)
395 c dsci_inv=dsc_inv(itypi)
396 dsci_inv=vbld_inv(i+nres)
398 C Calculate SC interaction energy.
401 do j=istart(i,iint),iend(i,iint)
404 if (itypj.eq.ntyp1) cycle
405 c dscj_inv=dsc_inv(itypj)
406 dscj_inv=vbld_inv(j+nres)
407 chi1=chi(itypi,itypj)
408 chi2=chi(itypj,itypi)
415 alf12=0.5D0*(alf1+alf2)
419 dxj=dc_norm(1,nres+j)
420 dyj=dc_norm(2,nres+j)
421 dzj=dc_norm(3,nres+j)
422 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
424 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
426 if (sss.lt.1.0d0) then
428 C Calculate the angle-dependent terms of energy & contributions to derivatives.
430 C Calculate whole angle-dependent part of epsilon and contributions
432 fac=(rrij*sigsq)**expon2
433 e1=fac*fac*aa(itypi,itypj)
434 e2=fac*bb(itypi,itypj)
435 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
436 eps2der=evdwij*eps3rt
437 eps3der=evdwij*eps2rt
438 evdwij=evdwij*eps2rt*eps3rt
439 evdw=evdw+evdwij*(1.0d0-sss)
441 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
442 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
443 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
444 cd & restyp(itypi),i,restyp(itypj),j,
445 cd & epsi,sigm,chi1,chi2,chip1,chip2,
446 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
447 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
450 C Calculate gradient components.
451 e1=e1*eps1*eps2rt**2*eps3rt**2
452 fac=-expon*(e1+evdwij)
455 C Calculate radial part of the gradient
459 C Calculate the angular part of the gradient and sum add the contributions
460 C to the appropriate components of the Cartesian gradient.
461 call sc_grad_scale(1.0d0-sss)
469 C-----------------------------------------------------------------------------
470 subroutine ebp_short(evdw)
472 C This subroutine calculates the interaction energy of nonbonded side chains
473 C assuming the Berne-Pechukas potential of interaction.
475 implicit real*8 (a-h,o-z)
479 include 'COMMON.LOCAL'
480 include 'COMMON.CHAIN'
481 include 'COMMON.DERIV'
482 include 'COMMON.NAMES'
483 include 'COMMON.INTERACT'
484 include 'COMMON.IOUNITS'
485 include 'COMMON.CALC'
487 c double precision rrsave(maxdim)
490 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
492 c if (icall.eq.0) then
500 if (itypi.eq.ntyp1) cycle
505 dxi=dc_norm(1,nres+i)
506 dyi=dc_norm(2,nres+i)
507 dzi=dc_norm(3,nres+i)
508 c dsci_inv=dsc_inv(itypi)
509 dsci_inv=vbld_inv(i+nres)
511 C Calculate SC interaction energy.
514 do j=istart(i,iint),iend(i,iint)
517 if (itypj.eq.ntyp1) cycle
518 c dscj_inv=dsc_inv(itypj)
519 dscj_inv=vbld_inv(j+nres)
520 chi1=chi(itypi,itypj)
521 chi2=chi(itypj,itypi)
528 alf12=0.5D0*(alf1+alf2)
532 dxj=dc_norm(1,nres+j)
533 dyj=dc_norm(2,nres+j)
534 dzj=dc_norm(3,nres+j)
535 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
537 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
539 if (sss.gt.0.0d0) then
541 C Calculate the angle-dependent terms of energy & contributions to derivatives.
543 C Calculate whole angle-dependent part of epsilon and contributions
545 fac=(rrij*sigsq)**expon2
546 e1=fac*fac*aa(itypi,itypj)
547 e2=fac*bb(itypi,itypj)
548 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
549 eps2der=evdwij*eps3rt
550 eps3der=evdwij*eps2rt
551 evdwij=evdwij*eps2rt*eps3rt
554 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
555 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
556 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
557 cd & restyp(itypi),i,restyp(itypj),j,
558 cd & epsi,sigm,chi1,chi2,chip1,chip2,
559 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
560 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
563 C Calculate gradient components.
564 e1=e1*eps1*eps2rt**2*eps3rt**2
565 fac=-expon*(e1+evdwij)
568 C Calculate radial part of the gradient
572 C Calculate the angular part of the gradient and sum add the contributions
573 C to the appropriate components of the Cartesian gradient.
574 call sc_grad_scale(sss)
582 C-----------------------------------------------------------------------------
583 subroutine egb_long(evdw)
585 C This subroutine calculates the interaction energy of nonbonded side chains
586 C assuming the Gay-Berne potential of interaction.
588 implicit real*8 (a-h,o-z)
592 include 'COMMON.LOCAL'
593 include 'COMMON.CHAIN'
594 include 'COMMON.DERIV'
595 include 'COMMON.NAMES'
596 include 'COMMON.INTERACT'
597 include 'COMMON.IOUNITS'
598 include 'COMMON.CALC'
599 include 'COMMON.CONTROL'
602 ccccc energy_dec=.false.
603 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
606 c if (icall.eq.0) lprn=.false.
610 if (itypi.eq.ntyp1) cycle
615 dxi=dc_norm(1,nres+i)
616 dyi=dc_norm(2,nres+i)
617 dzi=dc_norm(3,nres+i)
618 c dsci_inv=dsc_inv(itypi)
619 dsci_inv=vbld_inv(i+nres)
620 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
621 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
623 C Calculate SC interaction energy.
626 do j=istart(i,iint),iend(i,iint)
629 if (itypj.eq.ntyp1) cycle
630 c dscj_inv=dsc_inv(itypj)
631 dscj_inv=vbld_inv(j+nres)
632 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
633 c & 1.0d0/vbld(j+nres)
634 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
635 sig0ij=sigma(itypi,itypj)
636 chi1=chi(itypi,itypj)
637 chi2=chi(itypj,itypi)
644 alf12=0.5D0*(alf1+alf2)
648 dxj=dc_norm(1,nres+j)
649 dyj=dc_norm(2,nres+j)
650 dzj=dc_norm(3,nres+j)
651 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
653 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
655 if (sss.lt.1.0d0) then
657 C Calculate angle-dependent terms of energy and contributions to their
661 sig=sig0ij*dsqrt(sigsq)
662 rij_shift=1.0D0/rij-sig+sig0ij
663 c for diagnostics; uncomment
664 c rij_shift=1.2*sig0ij
665 C I hate to put IF's in the loops, but here don't have another choice!!!!
666 if (rij_shift.le.0.0D0) then
668 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
669 cd & restyp(itypi),i,restyp(itypj),j,
670 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
674 c---------------------------------------------------------------
675 rij_shift=1.0D0/rij_shift
677 e1=fac*fac*aa(itypi,itypj)
678 e2=fac*bb(itypi,itypj)
679 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
680 eps2der=evdwij*eps3rt
681 eps3der=evdwij*eps2rt
682 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
683 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
684 evdwij=evdwij*eps2rt*eps3rt
685 evdw=evdw+evdwij*(1.0d0-sss)
687 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
688 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
689 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
690 & restyp(itypi),i,restyp(itypj),j,
691 & epsi,sigm,chi1,chi2,chip1,chip2,
692 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
693 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
697 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
700 C Calculate gradient components.
701 e1=e1*eps1*eps2rt**2*eps3rt**2
702 fac=-expon*(e1+evdwij)*rij_shift
706 C Calculate the radial part of the gradient
710 C Calculate angular part of the gradient.
711 call sc_grad_scale(1.0d0-sss)
716 c write (iout,*) "Number of loop steps in EGB:",ind
717 cccc energy_dec=.false.
720 C-----------------------------------------------------------------------------
721 subroutine egb_short(evdw)
723 C This subroutine calculates the interaction energy of nonbonded side chains
724 C assuming the Gay-Berne potential of interaction.
726 implicit real*8 (a-h,o-z)
730 include 'COMMON.LOCAL'
731 include 'COMMON.CHAIN'
732 include 'COMMON.DERIV'
733 include 'COMMON.NAMES'
734 include 'COMMON.INTERACT'
735 include 'COMMON.IOUNITS'
736 include 'COMMON.CALC'
737 include 'COMMON.CONTROL'
740 ccccc energy_dec=.false.
741 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
744 c if (icall.eq.0) lprn=.false.
748 if (itypi.eq.ntyp1) cycle
753 dxi=dc_norm(1,nres+i)
754 dyi=dc_norm(2,nres+i)
755 dzi=dc_norm(3,nres+i)
756 c dsci_inv=dsc_inv(itypi)
757 dsci_inv=vbld_inv(i+nres)
758 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
759 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
761 C Calculate SC interaction energy.
764 do j=istart(i,iint),iend(i,iint)
767 if (itypj.eq.ntyp1) cycle
768 c dscj_inv=dsc_inv(itypj)
769 dscj_inv=vbld_inv(j+nres)
770 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
771 c & 1.0d0/vbld(j+nres)
772 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
773 sig0ij=sigma(itypi,itypj)
774 chi1=chi(itypi,itypj)
775 chi2=chi(itypj,itypi)
782 alf12=0.5D0*(alf1+alf2)
786 dxj=dc_norm(1,nres+j)
787 dyj=dc_norm(2,nres+j)
788 dzj=dc_norm(3,nres+j)
789 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
791 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
793 if (sss.gt.0.0d0) then
795 C Calculate angle-dependent terms of energy and contributions to their
799 sig=sig0ij*dsqrt(sigsq)
800 rij_shift=1.0D0/rij-sig+sig0ij
801 c for diagnostics; uncomment
802 c rij_shift=1.2*sig0ij
803 C I hate to put IF's in the loops, but here don't have another choice!!!!
804 if (rij_shift.le.0.0D0) then
806 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
807 cd & restyp(itypi),i,restyp(itypj),j,
808 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
812 c---------------------------------------------------------------
813 rij_shift=1.0D0/rij_shift
815 e1=fac*fac*aa(itypi,itypj)
816 e2=fac*bb(itypi,itypj)
817 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
818 eps2der=evdwij*eps3rt
819 eps3der=evdwij*eps2rt
820 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
821 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
822 evdwij=evdwij*eps2rt*eps3rt
825 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
826 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
827 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
828 & restyp(itypi),i,restyp(itypj),j,
829 & epsi,sigm,chi1,chi2,chip1,chip2,
830 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
831 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
835 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
838 C Calculate gradient components.
839 e1=e1*eps1*eps2rt**2*eps3rt**2
840 fac=-expon*(e1+evdwij)*rij_shift
844 C Calculate the radial part of the gradient
848 C Calculate angular part of the gradient.
849 call sc_grad_scale(sss)
854 c write (iout,*) "Number of loop steps in EGB:",ind
855 cccc energy_dec=.false.
858 C-----------------------------------------------------------------------------
859 subroutine egbv_long(evdw)
861 C This subroutine calculates the interaction energy of nonbonded side chains
862 C assuming the Gay-Berne-Vorobjev potential of interaction.
864 implicit real*8 (a-h,o-z)
868 include 'COMMON.LOCAL'
869 include 'COMMON.CHAIN'
870 include 'COMMON.DERIV'
871 include 'COMMON.NAMES'
872 include 'COMMON.INTERACT'
873 include 'COMMON.IOUNITS'
874 include 'COMMON.CALC'
878 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
881 c if (icall.eq.0) lprn=.true.
885 if (itypi.eq.ntyp1) cycle
890 dxi=dc_norm(1,nres+i)
891 dyi=dc_norm(2,nres+i)
892 dzi=dc_norm(3,nres+i)
893 c dsci_inv=dsc_inv(itypi)
894 dsci_inv=vbld_inv(i+nres)
896 C Calculate SC interaction energy.
899 do j=istart(i,iint),iend(i,iint)
902 if (itypj.eq.ntyp1) cycle
903 c dscj_inv=dsc_inv(itypj)
904 dscj_inv=vbld_inv(j+nres)
905 sig0ij=sigma(itypi,itypj)
907 chi1=chi(itypi,itypj)
908 chi2=chi(itypj,itypi)
915 alf12=0.5D0*(alf1+alf2)
919 dxj=dc_norm(1,nres+j)
920 dyj=dc_norm(2,nres+j)
921 dzj=dc_norm(3,nres+j)
922 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
925 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
927 if (sss.lt.1.0d0) then
929 C Calculate angle-dependent terms of energy and contributions to their
933 sig=sig0ij*dsqrt(sigsq)
934 rij_shift=1.0D0/rij-sig+r0ij
935 C I hate to put IF's in the loops, but here don't have another choice!!!!
936 if (rij_shift.le.0.0D0) then
941 c---------------------------------------------------------------
942 rij_shift=1.0D0/rij_shift
944 e1=fac*fac*aa(itypi,itypj)
945 e2=fac*bb(itypi,itypj)
946 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
947 eps2der=evdwij*eps3rt
948 eps3der=evdwij*eps2rt
950 e_augm=augm(itypi,itypj)*fac_augm
951 evdwij=evdwij*eps2rt*eps3rt
952 evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
954 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
955 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
956 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
957 & restyp(itypi),i,restyp(itypj),j,
958 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
959 & chi1,chi2,chip1,chip2,
960 & eps1,eps2rt**2,eps3rt**2,
961 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
964 C Calculate gradient components.
965 e1=e1*eps1*eps2rt**2*eps3rt**2
966 fac=-expon*(e1+evdwij)*rij_shift
968 fac=rij*fac-2*expon*rrij*e_augm
969 C Calculate the radial part of the gradient
973 C Calculate angular part of the gradient.
974 call sc_grad_scale(1.0d0-sss)
980 C-----------------------------------------------------------------------------
981 subroutine egbv_short(evdw)
983 C This subroutine calculates the interaction energy of nonbonded side chains
984 C assuming the Gay-Berne-Vorobjev potential of interaction.
986 implicit real*8 (a-h,o-z)
990 include 'COMMON.LOCAL'
991 include 'COMMON.CHAIN'
992 include 'COMMON.DERIV'
993 include 'COMMON.NAMES'
994 include 'COMMON.INTERACT'
995 include 'COMMON.IOUNITS'
996 include 'COMMON.CALC'
1000 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1003 c if (icall.eq.0) lprn=.true.
1005 do i=iatsc_s,iatsc_e
1007 if (itypi.eq.ntyp1) cycle
1012 dxi=dc_norm(1,nres+i)
1013 dyi=dc_norm(2,nres+i)
1014 dzi=dc_norm(3,nres+i)
1015 c dsci_inv=dsc_inv(itypi)
1016 dsci_inv=vbld_inv(i+nres)
1018 C Calculate SC interaction energy.
1020 do iint=1,nint_gr(i)
1021 do j=istart(i,iint),iend(i,iint)
1024 if (itypj.eq.ntyp1) cycle
1025 c dscj_inv=dsc_inv(itypj)
1026 dscj_inv=vbld_inv(j+nres)
1027 sig0ij=sigma(itypi,itypj)
1028 r0ij=r0(itypi,itypj)
1029 chi1=chi(itypi,itypj)
1030 chi2=chi(itypj,itypi)
1037 alf12=0.5D0*(alf1+alf2)
1041 dxj=dc_norm(1,nres+j)
1042 dyj=dc_norm(2,nres+j)
1043 dzj=dc_norm(3,nres+j)
1044 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1047 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1049 if (sss.gt.0.0d0) then
1051 C Calculate angle-dependent terms of energy and contributions to their
1055 sig=sig0ij*dsqrt(sigsq)
1056 rij_shift=1.0D0/rij-sig+r0ij
1057 C I hate to put IF's in the loops, but here don't have another choice!!!!
1058 if (rij_shift.le.0.0D0) then
1063 c---------------------------------------------------------------
1064 rij_shift=1.0D0/rij_shift
1065 fac=rij_shift**expon
1066 e1=fac*fac*aa(itypi,itypj)
1067 e2=fac*bb(itypi,itypj)
1068 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1069 eps2der=evdwij*eps3rt
1070 eps3der=evdwij*eps2rt
1071 fac_augm=rrij**expon
1072 e_augm=augm(itypi,itypj)*fac_augm
1073 evdwij=evdwij*eps2rt*eps3rt
1074 evdw=evdw+(evdwij+e_augm)*sss
1076 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1077 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1078 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1079 & restyp(itypi),i,restyp(itypj),j,
1080 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1081 & chi1,chi2,chip1,chip2,
1082 & eps1,eps2rt**2,eps3rt**2,
1083 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1086 C Calculate gradient components.
1087 e1=e1*eps1*eps2rt**2*eps3rt**2
1088 fac=-expon*(e1+evdwij)*rij_shift
1090 fac=rij*fac-2*expon*rrij*e_augm
1091 C Calculate the radial part of the gradient
1095 C Calculate angular part of the gradient.
1096 call sc_grad_scale(sss)
1102 C----------------------------------------------------------------------------
1103 subroutine sc_grad_scale(scalfac)
1104 implicit real*8 (a-h,o-z)
1105 include 'DIMENSIONS'
1106 include 'COMMON.CHAIN'
1107 include 'COMMON.DERIV'
1108 include 'COMMON.CALC'
1109 include 'COMMON.IOUNITS'
1110 double precision dcosom1(3),dcosom2(3)
1111 double precision scalfac
1112 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1113 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1114 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1115 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1119 c eom12=evdwij*eps1_om12
1121 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1122 c & " sigder",sigder
1123 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1124 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1126 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1127 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1130 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1132 c write (iout,*) "gg",(gg(k),k=1,3)
1134 gvdwx(k,i)=gvdwx(k,i)-gg(k)
1135 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1136 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1137 gvdwx(k,j)=gvdwx(k,j)+gg(k)
1138 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1139 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1140 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1141 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1142 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1143 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1146 C Calculate the components of the gradient in DC and X
1149 gvdwc(l,i)=gvdwc(l,i)-gg(l)
1150 gvdwc(l,j)=gvdwc(l,j)+gg(l)
1154 C--------------------------------------------------------------------------
1155 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1157 C This subroutine calculates the average interaction energy and its gradient
1158 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1159 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1160 C The potential depends both on the distance of peptide-group centers and on
1161 C the orientation of the CA-CA virtual bonds.
1163 implicit real*8 (a-h,o-z)
1167 include 'DIMENSIONS'
1168 include 'COMMON.CONTROL'
1169 include 'COMMON.SETUP'
1170 include 'COMMON.IOUNITS'
1171 include 'COMMON.GEO'
1172 include 'COMMON.VAR'
1173 include 'COMMON.LOCAL'
1174 include 'COMMON.CHAIN'
1175 include 'COMMON.DERIV'
1176 include 'COMMON.INTERACT'
1177 include 'COMMON.CONTACTS'
1178 include 'COMMON.TORSION'
1179 include 'COMMON.VECTORS'
1180 include 'COMMON.FFIELD'
1181 include 'COMMON.TIME1'
1182 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1183 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1184 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1185 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1186 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1187 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1189 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1191 double precision scal_el /1.0d0/
1193 double precision scal_el /0.5d0/
1196 C 13-go grudnia roku pamietnego...
1197 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1198 & 0.0d0,1.0d0,0.0d0,
1199 & 0.0d0,0.0d0,1.0d0/
1200 cd write(iout,*) 'In EELEC'
1202 cd write(iout,*) 'Type',i
1203 cd write(iout,*) 'B1',B1(:,i)
1204 cd write(iout,*) 'B2',B2(:,i)
1205 cd write(iout,*) 'CC',CC(:,:,i)
1206 cd write(iout,*) 'DD',DD(:,:,i)
1207 cd write(iout,*) 'EE',EE(:,:,i)
1209 cd call check_vecgrad
1211 if (icheckgrad.eq.1) then
1213 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1215 dc_norm(k,i)=dc(k,i)*fac
1217 c write (iout,*) 'i',i,' fac',fac
1220 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1221 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1222 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1223 c call vec_and_deriv
1229 time_mat=time_mat+MPI_Wtime()-time01
1233 cd write (iout,*) 'i=',i
1235 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1238 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1239 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1252 cd print '(a)','Enter EELEC'
1253 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1255 gel_loc_loc(i)=0.0d0
1260 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1262 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1264 do i=iturn3_start,iturn3_end
1265 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1266 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
1270 dx_normi=dc_norm(1,i)
1271 dy_normi=dc_norm(2,i)
1272 dz_normi=dc_norm(3,i)
1273 xmedi=c(1,i)+0.5d0*dxi
1274 ymedi=c(2,i)+0.5d0*dyi
1275 zmedi=c(3,i)+0.5d0*dzi
1277 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1278 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1279 num_cont_hb(i)=num_conti
1281 do i=iturn4_start,iturn4_end
1282 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1283 & .or. itype(i+3).eq.ntyp1
1284 & .or. itype(i+4).eq.ntyp1) cycle
1288 dx_normi=dc_norm(1,i)
1289 dy_normi=dc_norm(2,i)
1290 dz_normi=dc_norm(3,i)
1291 xmedi=c(1,i)+0.5d0*dxi
1292 ymedi=c(2,i)+0.5d0*dyi
1293 zmedi=c(3,i)+0.5d0*dzi
1294 num_conti=num_cont_hb(i)
1295 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1296 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1297 & call eturn4(i,eello_turn4)
1298 num_cont_hb(i)=num_conti
1301 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1303 do i=iatel_s,iatel_e
1304 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
1308 dx_normi=dc_norm(1,i)
1309 dy_normi=dc_norm(2,i)
1310 dz_normi=dc_norm(3,i)
1311 xmedi=c(1,i)+0.5d0*dxi
1312 ymedi=c(2,i)+0.5d0*dyi
1313 zmedi=c(3,i)+0.5d0*dzi
1314 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1315 num_conti=num_cont_hb(i)
1316 do j=ielstart(i),ielend(i)
1317 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
1318 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1320 num_cont_hb(i)=num_conti
1322 c write (iout,*) "Number of loop steps in EELEC:",ind
1324 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1325 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1327 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1328 ccc eel_loc=eel_loc+eello_turn3
1329 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1332 C-------------------------------------------------------------------------------
1333 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1334 implicit real*8 (a-h,o-z)
1335 include 'DIMENSIONS'
1339 include 'COMMON.CONTROL'
1340 include 'COMMON.IOUNITS'
1341 include 'COMMON.GEO'
1342 include 'COMMON.VAR'
1343 include 'COMMON.LOCAL'
1344 include 'COMMON.CHAIN'
1345 include 'COMMON.DERIV'
1346 include 'COMMON.INTERACT'
1347 include 'COMMON.CONTACTS'
1348 include 'COMMON.TORSION'
1349 include 'COMMON.VECTORS'
1350 include 'COMMON.FFIELD'
1351 include 'COMMON.TIME1'
1352 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1353 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1354 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1355 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1356 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1357 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1359 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1361 double precision scal_el /1.0d0/
1363 double precision scal_el /0.5d0/
1366 C 13-go grudnia roku pamietnego...
1367 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1368 & 0.0d0,1.0d0,0.0d0,
1369 & 0.0d0,0.0d0,1.0d0/
1370 c time00=MPI_Wtime()
1371 cd write (iout,*) "eelecij",i,j
1375 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1376 aaa=app(iteli,itelj)
1377 bbb=bpp(iteli,itelj)
1378 ael6i=ael6(iteli,itelj)
1379 ael3i=ael3(iteli,itelj)
1383 dx_normj=dc_norm(1,j)
1384 dy_normj=dc_norm(2,j)
1385 dz_normj=dc_norm(3,j)
1386 xj=c(1,j)+0.5D0*dxj-xmedi
1387 yj=c(2,j)+0.5D0*dyj-ymedi
1388 zj=c(3,j)+0.5D0*dzj-zmedi
1389 rij=xj*xj+yj*yj+zj*zj
1393 c For extracting the short-range part of Evdwpp
1394 sss=sscale(rij/rpp(iteli,itelj))
1398 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1399 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1400 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1401 fac=cosa-3.0D0*cosb*cosg
1403 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1404 if (j.eq.i+2) ev1=scal_el*ev1
1409 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1412 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1413 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1415 evdw1=evdw1+evdwij*(1.0d0-sss)
1416 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1417 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1418 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1419 cd & xmedi,ymedi,zmedi,xj,yj,zj
1421 if (energy_dec) then
1422 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1423 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1427 C Calculate contributions to the Cartesian gradient.
1430 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1431 facel=-3*rrmij*(el1+eesij)
1437 * Radial derivatives. First process both termini of the fragment (i,j)
1443 c ghalf=0.5D0*ggg(k)
1444 c gelc(k,i)=gelc(k,i)+ghalf
1445 c gelc(k,j)=gelc(k,j)+ghalf
1447 c 9/28/08 AL Gradient compotents will be summed only at the end
1449 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1450 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1453 * Loop over residues i+1 thru j-1.
1457 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1464 c ghalf=0.5D0*ggg(k)
1465 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1466 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1468 c 9/28/08 AL Gradient compotents will be summed only at the end
1470 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1471 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1474 * Loop over residues i+1 thru j-1.
1478 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1482 facvdw=ev1+evdwij*(1.0d0-sss)
1485 fac=-3*rrmij*(facvdw+facvdw+facel)
1490 * Radial derivatives. First process both termini of the fragment (i,j)
1496 c ghalf=0.5D0*ggg(k)
1497 c gelc(k,i)=gelc(k,i)+ghalf
1498 c gelc(k,j)=gelc(k,j)+ghalf
1500 c 9/28/08 AL Gradient compotents will be summed only at the end
1502 gelc_long(k,j)=gelc(k,j)+ggg(k)
1503 gelc_long(k,i)=gelc(k,i)-ggg(k)
1506 * Loop over residues i+1 thru j-1.
1510 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1513 c 9/28/08 AL Gradient compotents will be summed only at the end
1518 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1519 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1525 ecosa=2.0D0*fac3*fac1+fac4
1528 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1529 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1531 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1532 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1534 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1535 cd & (dcosg(k),k=1,3)
1537 ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1540 c ghalf=0.5D0*ggg(k)
1541 c gelc(k,i)=gelc(k,i)+ghalf
1542 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1543 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1544 c gelc(k,j)=gelc(k,j)+ghalf
1545 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1546 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1550 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1555 & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1556 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1558 & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1559 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1560 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1561 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1563 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1564 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1565 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1567 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1568 C energy of a peptide unit is assumed in the form of a second-order
1569 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1570 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1571 C are computed for EVERY pair of non-contiguous peptide groups.
1573 if (j.lt.nres-1) then
1584 muij(kkk)=mu(k,i)*mu(l,j)
1587 cd write (iout,*) 'EELEC: i',i,' j',j
1588 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1589 cd write(iout,*) 'muij',muij
1590 ury=scalar(uy(1,i),erij)
1591 urz=scalar(uz(1,i),erij)
1592 vry=scalar(uy(1,j),erij)
1593 vrz=scalar(uz(1,j),erij)
1594 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1595 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1596 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1597 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1598 fac=dsqrt(-ael6i)*r3ij
1603 cd write (iout,'(4i5,4f10.5)')
1604 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1605 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1606 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1607 cd & uy(:,j),uz(:,j)
1608 cd write (iout,'(4f10.5)')
1609 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1610 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1611 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
1612 cd write (iout,'(9f10.5/)')
1613 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1614 C Derivatives of the elements of A in virtual-bond vectors
1615 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1617 uryg(k,1)=scalar(erder(1,k),uy(1,i))
1618 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1619 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1620 urzg(k,1)=scalar(erder(1,k),uz(1,i))
1621 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1622 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1623 vryg(k,1)=scalar(erder(1,k),uy(1,j))
1624 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1625 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1626 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1627 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1628 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1630 C Compute radial contributions to the gradient
1648 C Add the contributions coming from er
1651 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1652 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1653 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1654 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1657 C Derivatives in DC(i)
1658 cgrad ghalf1=0.5d0*agg(k,1)
1659 cgrad ghalf2=0.5d0*agg(k,2)
1660 cgrad ghalf3=0.5d0*agg(k,3)
1661 cgrad ghalf4=0.5d0*agg(k,4)
1662 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1663 & -3.0d0*uryg(k,2)*vry)!+ghalf1
1664 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1665 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
1666 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1667 & -3.0d0*urzg(k,2)*vry)!+ghalf3
1668 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1669 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
1670 C Derivatives in DC(i+1)
1671 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
1672 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
1673 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
1674 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
1675 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
1676 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
1677 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
1678 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
1679 C Derivatives in DC(j)
1680 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
1681 & -3.0d0*vryg(k,2)*ury)!+ghalf1
1682 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
1683 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
1684 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
1685 & -3.0d0*vryg(k,2)*urz)!+ghalf3
1686 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
1687 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
1688 C Derivatives in DC(j+1) or DC(nres-1)
1689 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
1690 & -3.0d0*vryg(k,3)*ury)
1691 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
1692 & -3.0d0*vrzg(k,3)*ury)
1693 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
1694 & -3.0d0*vryg(k,3)*urz)
1695 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
1696 & -3.0d0*vrzg(k,3)*urz)
1697 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
1699 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
1712 aggi(k,l)=-aggi(k,l)
1713 aggi1(k,l)=-aggi1(k,l)
1714 aggj(k,l)=-aggj(k,l)
1715 aggj1(k,l)=-aggj1(k,l)
1718 if (j.lt.nres-1) then
1724 aggi(k,l)=-aggi(k,l)
1725 aggi1(k,l)=-aggi1(k,l)
1726 aggj(k,l)=-aggj(k,l)
1727 aggj1(k,l)=-aggj1(k,l)
1738 aggi(k,l)=-aggi(k,l)
1739 aggi1(k,l)=-aggi1(k,l)
1740 aggj(k,l)=-aggj(k,l)
1741 aggj1(k,l)=-aggj1(k,l)
1746 IF (wel_loc.gt.0.0d0) THEN
1747 C Contribution to the local-electrostatic energy coming from the i-j pair
1748 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
1750 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
1752 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1753 & 'eelloc',i,j,eel_loc_ij
1755 eel_loc=eel_loc+eel_loc_ij
1756 C Partial derivatives in virtual-bond dihedral angles gamma
1758 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
1759 & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1760 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1761 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
1762 & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1763 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1764 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
1766 ggg(l)=agg(l,1)*muij(1)+
1767 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1768 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
1769 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
1770 cgrad ghalf=0.5d0*ggg(l)
1771 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
1772 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
1776 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
1779 C Remaining derivatives of eello
1781 gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1782 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1783 gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1784 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1785 gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1786 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1787 gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1788 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1791 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
1792 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
1793 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
1794 & .and. num_conti.le.maxconts) then
1795 c write (iout,*) i,j," entered corr"
1797 C Calculate the contact function. The ith column of the array JCONT will
1798 C contain the numbers of atoms that make contacts with the atom I (of numbers
1799 C greater than I). The arrays FACONT and GACONT will contain the values of
1800 C the contact function and its derivative.
1801 c r0ij=1.02D0*rpp(iteli,itelj)
1802 c r0ij=1.11D0*rpp(iteli,itelj)
1803 r0ij=2.20D0*rpp(iteli,itelj)
1804 c r0ij=1.55D0*rpp(iteli,itelj)
1805 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
1806 if (fcont.gt.0.0D0) then
1807 num_conti=num_conti+1
1808 if (num_conti.gt.maxconts) then
1809 write (iout,*) 'WARNING - max. # of contacts exceeded;',
1810 & ' will skip next contacts for this conf.'
1812 jcont_hb(num_conti,i)=j
1813 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
1814 cd & " jcont_hb",jcont_hb(num_conti,i)
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)
1835 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
1836 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
1837 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
1838 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
1839 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
1844 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
1845 C Calculate contact energies
1847 wij=cosa-3.0D0*cosb*cosg
1850 c fac3=dsqrt(-ael6i)/r0ij**3
1851 fac3=dsqrt(-ael6i)*r3ij
1852 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
1853 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
1854 if (ees0tmp.gt.0) then
1855 ees0pij=dsqrt(ees0tmp)
1859 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
1860 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
1861 if (ees0tmp.gt.0) then
1862 ees0mij=dsqrt(ees0tmp)
1867 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
1868 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
1869 C Diagnostics. Comment out or remove after debugging!
1870 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
1871 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
1872 c ees0m(num_conti,i)=0.0D0
1874 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
1875 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
1876 C Angular derivatives of the contact function
1877 ees0pij1=fac3/ees0pij
1878 ees0mij1=fac3/ees0mij
1879 fac3p=-3.0D0*fac3*rrmij
1880 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
1881 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
1883 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
1884 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
1885 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
1886 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
1887 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
1888 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
1889 ecosap=ecosa1+ecosa2
1890 ecosbp=ecosb1+ecosb2
1891 ecosgp=ecosg1+ecosg2
1892 ecosam=ecosa1-ecosa2
1893 ecosbm=ecosb1-ecosb2
1894 ecosgm=ecosg1-ecosg2
1903 facont_hb(num_conti,i)=fcont
1904 fprimcont=fprimcont/rij
1905 cd facont_hb(num_conti,i)=1.0D0
1906 C Following line is for diagnostics.
1909 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1910 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1913 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
1914 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
1916 gggp(1)=gggp(1)+ees0pijp*xj
1917 gggp(2)=gggp(2)+ees0pijp*yj
1918 gggp(3)=gggp(3)+ees0pijp*zj
1919 gggm(1)=gggm(1)+ees0mijp*xj
1920 gggm(2)=gggm(2)+ees0mijp*yj
1921 gggm(3)=gggm(3)+ees0mijp*zj
1922 C Derivatives due to the contact function
1923 gacont_hbr(1,num_conti,i)=fprimcont*xj
1924 gacont_hbr(2,num_conti,i)=fprimcont*yj
1925 gacont_hbr(3,num_conti,i)=fprimcont*zj
1928 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
1929 c following the change of gradient-summation algorithm.
1931 cgrad ghalfp=0.5D0*gggp(k)
1932 cgrad ghalfm=0.5D0*gggm(k)
1933 gacontp_hb1(k,num_conti,i)=!ghalfp
1934 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
1935 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1936 gacontp_hb2(k,num_conti,i)=!ghalfp
1937 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
1938 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1939 gacontp_hb3(k,num_conti,i)=gggp(k)
1940 gacontm_hb1(k,num_conti,i)=!ghalfm
1941 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
1942 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1943 gacontm_hb2(k,num_conti,i)=!ghalfm
1944 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
1945 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1946 gacontm_hb3(k,num_conti,i)=gggm(k)
1949 endif ! num_conti.le.maxconts
1952 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
1955 ghalf=0.5d0*agg(l,k)
1956 aggi(l,k)=aggi(l,k)+ghalf
1957 aggi1(l,k)=aggi1(l,k)+agg(l,k)
1958 aggj(l,k)=aggj(l,k)+ghalf
1961 if (j.eq.nres-1 .and. i.lt.j-2) then
1964 aggj1(l,k)=aggj1(l,k)+agg(l,k)
1969 c t_eelecij=t_eelecij+MPI_Wtime()-time00
1972 C-----------------------------------------------------------------------
1973 subroutine evdwpp_short(evdw1)
1977 implicit real*8 (a-h,o-z)
1978 include 'DIMENSIONS'
1979 include 'COMMON.CONTROL'
1980 include 'COMMON.IOUNITS'
1981 include 'COMMON.GEO'
1982 include 'COMMON.VAR'
1983 include 'COMMON.LOCAL'
1984 include 'COMMON.CHAIN'
1985 include 'COMMON.DERIV'
1986 include 'COMMON.INTERACT'
1987 include 'COMMON.CONTACTS'
1988 include 'COMMON.TORSION'
1989 include 'COMMON.VECTORS'
1990 include 'COMMON.FFIELD'
1992 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1994 double precision scal_el /1.0d0/
1996 double precision scal_el /0.5d0/
1999 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2000 c & " iatel_e_vdw",iatel_e_vdw
2002 do i=iatel_s_vdw,iatel_e_vdw
2003 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2007 dx_normi=dc_norm(1,i)
2008 dy_normi=dc_norm(2,i)
2009 dz_normi=dc_norm(3,i)
2010 xmedi=c(1,i)+0.5d0*dxi
2011 ymedi=c(2,i)+0.5d0*dyi
2012 zmedi=c(3,i)+0.5d0*dzi
2014 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2015 c & ' ielend',ielend_vdw(i)
2017 do j=ielstart_vdw(i),ielend_vdw(i)
2018 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2022 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2023 aaa=app(iteli,itelj)
2024 bbb=bpp(iteli,itelj)
2028 dx_normj=dc_norm(1,j)
2029 dy_normj=dc_norm(2,j)
2030 dz_normj=dc_norm(3,j)
2031 xj=c(1,j)+0.5D0*dxj-xmedi
2032 yj=c(2,j)+0.5D0*dyj-ymedi
2033 zj=c(3,j)+0.5D0*dzj-zmedi
2034 rij=xj*xj+yj*yj+zj*zj
2037 sss=sscale(rij/rpp(iteli,itelj))
2038 if (sss.gt.0.0d0) then
2043 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2044 if (j.eq.i+2) ev1=scal_el*ev1
2047 if (energy_dec) then
2048 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2050 evdw1=evdw1+evdwij*sss
2052 C Calculate contributions to the Cartesian gradient.
2054 facvdw=-6*rrmij*(ev1+evdwij)*sss
2059 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2060 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2067 C-----------------------------------------------------------------------------
2068 subroutine escp_long(evdw2,evdw2_14)
2070 C This subroutine calculates the excluded-volume interaction energy between
2071 C peptide-group centers and side chains and its gradient in virtual-bond and
2072 C side-chain vectors.
2074 implicit real*8 (a-h,o-z)
2075 include 'DIMENSIONS'
2076 include 'COMMON.GEO'
2077 include 'COMMON.VAR'
2078 include 'COMMON.LOCAL'
2079 include 'COMMON.CHAIN'
2080 include 'COMMON.DERIV'
2081 include 'COMMON.INTERACT'
2082 include 'COMMON.FFIELD'
2083 include 'COMMON.IOUNITS'
2084 include 'COMMON.CONTROL'
2088 cd print '(a)','Enter ESCP'
2089 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2090 do i=iatscp_s,iatscp_e
2091 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2093 xi=0.5D0*(c(1,i)+c(1,i+1))
2094 yi=0.5D0*(c(2,i)+c(2,i+1))
2095 zi=0.5D0*(c(3,i)+c(3,i+1))
2097 do iint=1,nscp_gr(i)
2099 do j=iscpstart(i,iint),iscpend(i,iint)
2101 if (itypj.eq.ntyp1) cycle
2102 C Uncomment following three lines for SC-p interactions
2106 C Uncomment following three lines for Ca-p interactions
2110 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2112 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2114 if (sss.lt.1.0d0) then
2117 e1=fac*fac*aad(itypj,iteli)
2118 e2=fac*bad(itypj,iteli)
2119 if (iabs(j-i) .le. 2) then
2122 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2125 evdw2=evdw2+evdwij*(1.0d0-sss)
2126 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2127 & 'evdw2',i,j,sss,evdwij
2129 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2131 fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2135 C Uncomment following three lines for SC-p interactions
2137 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2139 C Uncomment following line for SC-p interactions
2140 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2142 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2143 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2152 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2153 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2154 gradx_scp(j,i)=expon*gradx_scp(j,i)
2157 C******************************************************************************
2161 C To save time the factor EXPON has been extracted from ALL components
2162 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2165 C******************************************************************************
2168 C-----------------------------------------------------------------------------
2169 subroutine escp_short(evdw2,evdw2_14)
2171 C This subroutine calculates the excluded-volume interaction energy between
2172 C peptide-group centers and side chains and its gradient in virtual-bond and
2173 C side-chain vectors.
2175 implicit real*8 (a-h,o-z)
2176 include 'DIMENSIONS'
2177 include 'COMMON.GEO'
2178 include 'COMMON.VAR'
2179 include 'COMMON.LOCAL'
2180 include 'COMMON.CHAIN'
2181 include 'COMMON.DERIV'
2182 include 'COMMON.INTERACT'
2183 include 'COMMON.FFIELD'
2184 include 'COMMON.IOUNITS'
2185 include 'COMMON.CONTROL'
2189 cd print '(a)','Enter ESCP'
2190 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2191 do i=iatscp_s,iatscp_e
2192 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2194 xi=0.5D0*(c(1,i)+c(1,i+1))
2195 yi=0.5D0*(c(2,i)+c(2,i+1))
2196 zi=0.5D0*(c(3,i)+c(3,i+1))
2198 do iint=1,nscp_gr(i)
2200 do j=iscpstart(i,iint),iscpend(i,iint)
2202 if (itypj.eq.ntyp1) cycle
2203 C Uncomment following three lines for SC-p interactions
2207 C Uncomment following three lines for Ca-p interactions
2211 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2213 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2215 if (sss.gt.0.0d0) then
2218 e1=fac*fac*aad(itypj,iteli)
2219 e2=fac*bad(itypj,iteli)
2220 if (iabs(j-i) .le. 2) then
2223 evdw2_14=evdw2_14+(e1+e2)*sss
2226 evdw2=evdw2+evdwij*sss
2227 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2228 & 'evdw2',i,j,sss,evdwij
2230 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2232 fac=-(evdwij+e1)*rrij*sss
2236 C Uncomment following three lines for SC-p interactions
2238 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2240 C Uncomment following line for SC-p interactions
2241 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2243 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2244 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2253 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2254 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2255 gradx_scp(j,i)=expon*gradx_scp(j,i)
2258 C******************************************************************************
2262 C To save time the factor EXPON has been extracted from ALL components
2263 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2266 C******************************************************************************