1 C-----------------------------------------------------------------------
2 double precision function sscalelip(r)
3 double precision r,gamm
4 include "COMMON.SPLITELE"
5 C if(r.lt.r_cut-rlamb) then
7 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
8 C gamm=(r-(r_cut-rlamb))/rlamb
9 sscalelip=1.0d0+r*r*(2*r-3.0d0)
15 C-----------------------------------------------------------------------
16 double precision function sscagradlip(r)
17 double precision r,gamm
18 include "COMMON.SPLITELE"
19 C if(r.lt.r_cut-rlamb) then
21 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
22 C gamm=(r-(r_cut-rlamb))/rlamb
23 sscagradlip=r*(6*r-6.0d0)
30 C-----------------------------------------------------------------------
31 double precision function sscale(r)
32 double precision r,gamm
33 include "COMMON.SPLITELE"
34 if(r.lt.r_cut-rlamb) then
36 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
37 gamm=(r-(r_cut-rlamb))/rlamb
38 sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
44 C-----------------------------------------------------------------------
45 C-----------------------------------------------------------------------
46 double precision function sscagrad(r)
47 double precision r,gamm
48 include "COMMON.SPLITELE"
49 if(r.lt.r_cut-rlamb) then
51 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
52 gamm=(r-(r_cut-rlamb))/rlamb
53 sscagrad=gamm*(6*gamm-6.0d0)/rlamb
59 C-----------------------------------------------------------------------
60 subroutine elj_long(evdw)
62 C This subroutine calculates the interaction energy of nonbonded side chains
63 C assuming the LJ potential of interaction.
65 implicit real*8 (a-h,o-z)
67 parameter (accur=1.0d-10)
70 include 'COMMON.LOCAL'
71 include 'COMMON.CHAIN'
72 include 'COMMON.DERIV'
73 include 'COMMON.INTERACT'
74 include 'COMMON.TORSION'
75 include 'COMMON.SBRIDGE'
76 include 'COMMON.NAMES'
77 include 'COMMON.IOUNITS'
78 include 'COMMON.CONTACTS'
80 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
84 if (itypi.eq.ntyp1) cycle
90 C Calculate SC interaction energy.
93 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
94 cd & 'iend=',iend(i,iint)
95 do j=istart(i,iint),iend(i,iint)
97 if (itypj.eq.ntyp1) cycle
101 rij=xj*xj+yj*yj+zj*zj
102 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
103 if (sss.lt.1.0d0) then
105 eps0ij=eps(itypi,itypj)
107 e1=fac*fac*aa(itypi,itypj)
108 e2=fac*bb(itypi,itypj)
110 evdw=evdw+(1.0d0-sss)*evdwij
112 C Calculate the components of the gradient in DC and X
114 fac=-rrij*(e1+evdwij)*(1.0d0-sss)
119 gvdwx(k,i)=gvdwx(k,i)-gg(k)
120 gvdwx(k,j)=gvdwx(k,j)+gg(k)
121 gvdwc(k,i)=gvdwc(k,i)-gg(k)
122 gvdwc(k,j)=gvdwc(k,j)+gg(k)
130 gvdwc(j,i)=expon*gvdwc(j,i)
131 gvdwx(j,i)=expon*gvdwx(j,i)
134 C******************************************************************************
138 C To save time, the factor of EXPON has been extracted from ALL components
139 C of GVDWC and GRADX. Remember to multiply them by this factor before further
142 C******************************************************************************
145 C-----------------------------------------------------------------------
146 subroutine elj_short(evdw)
148 C This subroutine calculates the interaction energy of nonbonded side chains
149 C assuming the LJ potential of interaction.
151 implicit real*8 (a-h,o-z)
153 parameter (accur=1.0d-10)
156 include 'COMMON.LOCAL'
157 include 'COMMON.CHAIN'
158 include 'COMMON.DERIV'
159 include 'COMMON.INTERACT'
160 include 'COMMON.TORSION'
161 include 'COMMON.SBRIDGE'
162 include 'COMMON.NAMES'
163 include 'COMMON.IOUNITS'
164 include 'COMMON.CONTACTS'
166 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
170 if (itypi.eq.ntyp1) cycle
178 C Calculate SC interaction energy.
181 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
182 cd & 'iend=',iend(i,iint)
183 do j=istart(i,iint),iend(i,iint)
185 if (itypj.eq.ntyp1) cycle
189 C Change 12/1/95 to calculate four-body interactions
190 rij=xj*xj+yj*yj+zj*zj
191 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
192 if (sss.gt.0.0d0) then
194 eps0ij=eps(itypi,itypj)
196 e1=fac*fac*aa(itypi,itypj)
197 e2=fac*bb(itypi,itypj)
201 C Calculate the components of the gradient in DC and X
203 fac=-rrij*(e1+evdwij)*sss
208 gvdwx(k,i)=gvdwx(k,i)-gg(k)
209 gvdwx(k,j)=gvdwx(k,j)+gg(k)
210 gvdwc(k,i)=gvdwc(k,i)-gg(k)
211 gvdwc(k,j)=gvdwc(k,j)+gg(k)
219 gvdwc(j,i)=expon*gvdwc(j,i)
220 gvdwx(j,i)=expon*gvdwx(j,i)
223 C******************************************************************************
227 C To save time, the factor of EXPON has been extracted from ALL components
228 C of GVDWC and GRADX. Remember to multiply them by this factor before further
231 C******************************************************************************
234 C-----------------------------------------------------------------------------
235 subroutine eljk_long(evdw)
237 C This subroutine calculates the interaction energy of nonbonded side chains
238 C assuming the LJK potential of interaction.
240 implicit real*8 (a-h,o-z)
244 include 'COMMON.LOCAL'
245 include 'COMMON.CHAIN'
246 include 'COMMON.DERIV'
247 include 'COMMON.INTERACT'
248 include 'COMMON.IOUNITS'
249 include 'COMMON.NAMES'
252 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
256 if (itypi.eq.ntyp1) cycle
262 C Calculate SC interaction energy.
265 do j=istart(i,iint),iend(i,iint)
267 if (itypj.eq.ntyp1) cycle
271 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
273 e_augm=augm(itypi,itypj)*fac_augm
276 sss=sscale(rij/sigma(itypi,itypj))
277 if (sss.lt.1.0d0) then
278 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
279 fac=r_shift_inv**expon
280 e1=fac*fac*aa(itypi,itypj)
281 e2=fac*bb(itypi,itypj)
283 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
284 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
285 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
286 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
287 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
288 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
289 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
290 evdw=evdw+(1.0d0-sss)*evdwij
292 C Calculate the components of the gradient in DC and X
294 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
300 gvdwx(k,i)=gvdwx(k,i)-gg(k)
301 gvdwx(k,j)=gvdwx(k,j)+gg(k)
302 gvdwc(k,i)=gvdwc(k,i)-gg(k)
303 gvdwc(k,j)=gvdwc(k,j)+gg(k)
311 gvdwc(j,i)=expon*gvdwc(j,i)
312 gvdwx(j,i)=expon*gvdwx(j,i)
317 C-----------------------------------------------------------------------------
318 subroutine eljk_short(evdw)
320 C This subroutine calculates the interaction energy of nonbonded side chains
321 C assuming the LJK potential of interaction.
323 implicit real*8 (a-h,o-z)
327 include 'COMMON.LOCAL'
328 include 'COMMON.CHAIN'
329 include 'COMMON.DERIV'
330 include 'COMMON.INTERACT'
331 include 'COMMON.IOUNITS'
332 include 'COMMON.NAMES'
335 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
339 if (itypi.eq.ntyp1) cycle
345 C Calculate SC interaction energy.
348 do j=istart(i,iint),iend(i,iint)
350 if (itypj.eq.ntyp1) cycle
354 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
356 e_augm=augm(itypi,itypj)*fac_augm
359 sss=sscale(rij/sigma(itypi,itypj))
360 if (sss.gt.0.0d0) then
361 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
362 fac=r_shift_inv**expon
363 e1=fac*fac*aa(itypi,itypj)
364 e2=fac*bb(itypi,itypj)
366 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
367 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
368 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
369 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
370 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
371 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
372 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
375 C Calculate the components of the gradient in DC and X
377 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
383 gvdwx(k,i)=gvdwx(k,i)-gg(k)
384 gvdwx(k,j)=gvdwx(k,j)+gg(k)
385 gvdwc(k,i)=gvdwc(k,i)-gg(k)
386 gvdwc(k,j)=gvdwc(k,j)+gg(k)
394 gvdwc(j,i)=expon*gvdwc(j,i)
395 gvdwx(j,i)=expon*gvdwx(j,i)
400 C-----------------------------------------------------------------------------
401 subroutine ebp_long(evdw)
403 C This subroutine calculates the interaction energy of nonbonded side chains
404 C assuming the Berne-Pechukas potential of interaction.
406 implicit real*8 (a-h,o-z)
410 include 'COMMON.LOCAL'
411 include 'COMMON.CHAIN'
412 include 'COMMON.DERIV'
413 include 'COMMON.NAMES'
414 include 'COMMON.INTERACT'
415 include 'COMMON.IOUNITS'
416 include 'COMMON.CALC'
418 c double precision rrsave(maxdim)
421 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
423 c if (icall.eq.0) then
431 if (itypi.eq.ntyp1) cycle
436 dxi=dc_norm(1,nres+i)
437 dyi=dc_norm(2,nres+i)
438 dzi=dc_norm(3,nres+i)
439 c dsci_inv=dsc_inv(itypi)
440 dsci_inv=vbld_inv(i+nres)
442 C Calculate SC interaction energy.
445 do j=istart(i,iint),iend(i,iint)
448 if (itypj.eq.ntyp1) cycle
449 c dscj_inv=dsc_inv(itypj)
450 dscj_inv=vbld_inv(j+nres)
451 chi1=chi(itypi,itypj)
452 chi2=chi(itypj,itypi)
459 alf12=0.5D0*(alf1+alf2)
463 dxj=dc_norm(1,nres+j)
464 dyj=dc_norm(2,nres+j)
465 dzj=dc_norm(3,nres+j)
466 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
468 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
470 if (sss.lt.1.0d0) then
472 C Calculate the angle-dependent terms of energy & contributions to derivatives.
474 C Calculate whole angle-dependent part of epsilon and contributions
476 fac=(rrij*sigsq)**expon2
477 e1=fac*fac*aa(itypi,itypj)
478 e2=fac*bb(itypi,itypj)
479 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
480 eps2der=evdwij*eps3rt
481 eps3der=evdwij*eps2rt
482 evdwij=evdwij*eps2rt*eps3rt
483 evdw=evdw+evdwij*(1.0d0-sss)
485 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
486 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
487 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
488 cd & restyp(itypi),i,restyp(itypj),j,
489 cd & epsi,sigm,chi1,chi2,chip1,chip2,
490 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
491 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
494 C Calculate gradient components.
495 e1=e1*eps1*eps2rt**2*eps3rt**2
496 fac=-expon*(e1+evdwij)
499 C Calculate radial part of the gradient
503 C Calculate the angular part of the gradient and sum add the contributions
504 C to the appropriate components of the Cartesian gradient.
505 call sc_grad_scale(1.0d0-sss)
513 C-----------------------------------------------------------------------------
514 subroutine ebp_short(evdw)
516 C This subroutine calculates the interaction energy of nonbonded side chains
517 C assuming the Berne-Pechukas potential of interaction.
519 implicit real*8 (a-h,o-z)
523 include 'COMMON.LOCAL'
524 include 'COMMON.CHAIN'
525 include 'COMMON.DERIV'
526 include 'COMMON.NAMES'
527 include 'COMMON.INTERACT'
528 include 'COMMON.IOUNITS'
529 include 'COMMON.CALC'
531 c double precision rrsave(maxdim)
534 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
536 c if (icall.eq.0) then
544 if (itypi.eq.ntyp1) cycle
549 dxi=dc_norm(1,nres+i)
550 dyi=dc_norm(2,nres+i)
551 dzi=dc_norm(3,nres+i)
552 c dsci_inv=dsc_inv(itypi)
553 dsci_inv=vbld_inv(i+nres)
555 C Calculate SC interaction energy.
558 do j=istart(i,iint),iend(i,iint)
561 if (itypj.eq.ntyp1) cycle
562 c dscj_inv=dsc_inv(itypj)
563 dscj_inv=vbld_inv(j+nres)
564 chi1=chi(itypi,itypj)
565 chi2=chi(itypj,itypi)
572 alf12=0.5D0*(alf1+alf2)
576 dxj=dc_norm(1,nres+j)
577 dyj=dc_norm(2,nres+j)
578 dzj=dc_norm(3,nres+j)
579 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
581 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
583 if (sss.gt.0.0d0) then
585 C Calculate the angle-dependent terms of energy & contributions to derivatives.
587 C Calculate whole angle-dependent part of epsilon and contributions
589 fac=(rrij*sigsq)**expon2
590 e1=fac*fac*aa(itypi,itypj)
591 e2=fac*bb(itypi,itypj)
592 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
593 eps2der=evdwij*eps3rt
594 eps3der=evdwij*eps2rt
595 evdwij=evdwij*eps2rt*eps3rt
598 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
599 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
600 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
601 cd & restyp(itypi),i,restyp(itypj),j,
602 cd & epsi,sigm,chi1,chi2,chip1,chip2,
603 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
604 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
607 C Calculate gradient components.
608 e1=e1*eps1*eps2rt**2*eps3rt**2
609 fac=-expon*(e1+evdwij)
612 C Calculate radial part of the gradient
616 C Calculate the angular part of the gradient and sum add the contributions
617 C to the appropriate components of the Cartesian gradient.
618 call sc_grad_scale(sss)
626 C-----------------------------------------------------------------------------
627 subroutine egb_long(evdw)
629 C This subroutine calculates the interaction energy of nonbonded side chains
630 C assuming the Gay-Berne potential of interaction.
632 implicit real*8 (a-h,o-z)
636 include 'COMMON.LOCAL'
637 include 'COMMON.CHAIN'
638 include 'COMMON.DERIV'
639 include 'COMMON.NAMES'
640 include 'COMMON.INTERACT'
641 include 'COMMON.IOUNITS'
642 include 'COMMON.CALC'
643 include 'COMMON.CONTROL'
646 ccccc energy_dec=.false.
647 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
650 c if (icall.eq.0) lprn=.false.
654 if (itypi.eq.ntyp1) cycle
660 if (xi.lt.0) xi=xi+boxxsize
662 if (yi.lt.0) yi=yi+boxysize
664 if (zi.lt.0) zi=zi+boxzsize
665 dxi=dc_norm(1,nres+i)
666 dyi=dc_norm(2,nres+i)
667 dzi=dc_norm(3,nres+i)
668 c dsci_inv=dsc_inv(itypi)
669 dsci_inv=vbld_inv(i+nres)
670 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
671 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
673 C Calculate SC interaction energy.
676 do j=istart(i,iint),iend(i,iint)
679 if (itypj.eq.ntyp1) cycle
680 c dscj_inv=dsc_inv(itypj)
681 dscj_inv=vbld_inv(j+nres)
682 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
683 c & 1.0d0/vbld(j+nres)
684 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
685 sig0ij=sigma(itypi,itypj)
686 chi1=chi(itypi,itypj)
687 chi2=chi(itypj,itypi)
694 alf12=0.5D0*(alf1+alf2)
699 if (xj.lt.0) xj=xj+boxxsize
701 if (yj.lt.0) yj=yj+boxysize
703 if (zj.lt.0) zj=zj+boxzsize
704 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
712 xj=xj_safe+xshift*boxxsize
713 yj=yj_safe+yshift*boxysize
714 zj=zj_safe+zshift*boxzsize
715 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
716 if(dist_temp.lt.dist_init) then
726 if (subchap.eq.1) then
735 dxj=dc_norm(1,nres+j)
736 dyj=dc_norm(2,nres+j)
737 dzj=dc_norm(3,nres+j)
738 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
740 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
741 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
742 if (sss.lt.1.0d0) then
744 C Calculate angle-dependent terms of energy and contributions to their
748 sig=sig0ij*dsqrt(sigsq)
749 rij_shift=1.0D0/rij-sig+sig0ij
750 c for diagnostics; uncomment
751 c rij_shift=1.2*sig0ij
752 C I hate to put IF's in the loops, but here don't have another choice!!!!
753 if (rij_shift.le.0.0D0) then
755 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
756 cd & restyp(itypi),i,restyp(itypj),j,
757 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
761 c---------------------------------------------------------------
762 rij_shift=1.0D0/rij_shift
764 e1=fac*fac*aa(itypi,itypj)
765 e2=fac*bb(itypi,itypj)
766 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
767 eps2der=evdwij*eps3rt
768 eps3der=evdwij*eps2rt
769 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
770 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
771 evdwij=evdwij*eps2rt*eps3rt
772 evdw=evdw+evdwij*(1.0d0-sss)
774 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
775 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
776 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
777 & restyp(itypi),i,restyp(itypj),j,
778 & epsi,sigm,chi1,chi2,chip1,chip2,
779 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
780 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
784 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
787 C Calculate gradient components.
788 e1=e1*eps1*eps2rt**2*eps3rt**2
789 fac=-expon*(e1+evdwij)*rij_shift
792 fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij
794 C Calculate the radial part of the gradient
798 C Calculate angular part of the gradient.
799 call sc_grad_scale(1.0d0-sss)
804 c write (iout,*) "Number of loop steps in EGB:",ind
805 cccc energy_dec=.false.
808 C-----------------------------------------------------------------------------
809 subroutine egb_short(evdw)
811 C This subroutine calculates the interaction energy of nonbonded side chains
812 C assuming the Gay-Berne potential of interaction.
814 implicit real*8 (a-h,o-z)
818 include 'COMMON.LOCAL'
819 include 'COMMON.CHAIN'
820 include 'COMMON.DERIV'
821 include 'COMMON.NAMES'
822 include 'COMMON.INTERACT'
823 include 'COMMON.IOUNITS'
824 include 'COMMON.CALC'
825 include 'COMMON.CONTROL'
828 ccccc energy_dec=.false.
829 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
832 c if (icall.eq.0) lprn=.false.
836 if (itypi.eq.ntyp1) cycle
842 if (xi.lt.0) xi=xi+boxxsize
844 if (yi.lt.0) yi=yi+boxysize
846 if (zi.lt.0) zi=zi+boxzsize
847 dxi=dc_norm(1,nres+i)
848 dyi=dc_norm(2,nres+i)
849 dzi=dc_norm(3,nres+i)
850 c dsci_inv=dsc_inv(itypi)
851 dsci_inv=vbld_inv(i+nres)
852 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
853 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
855 C Calculate SC interaction energy.
858 do j=istart(i,iint),iend(i,iint)
861 if (itypj.eq.ntyp1) cycle
862 c dscj_inv=dsc_inv(itypj)
863 dscj_inv=vbld_inv(j+nres)
864 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
865 c & 1.0d0/vbld(j+nres)
866 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
867 sig0ij=sigma(itypi,itypj)
868 chi1=chi(itypi,itypj)
869 chi2=chi(itypj,itypi)
876 alf12=0.5D0*(alf1+alf2)
881 if (xj.lt.0) xj=xj+boxxsize
883 if (yj.lt.0) yj=yj+boxysize
885 if (zj.lt.0) zj=zj+boxzsize
886 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
894 xj=xj_safe+xshift*boxxsize
895 yj=yj_safe+yshift*boxysize
896 zj=zj_safe+zshift*boxzsize
897 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
898 if(dist_temp.lt.dist_init) then
908 if (subchap.eq.1) then
917 dxj=dc_norm(1,nres+j)
918 dyj=dc_norm(2,nres+j)
919 dzj=dc_norm(3,nres+j)
920 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
922 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
923 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
924 if (sss.gt.0.0d0) then
926 C Calculate angle-dependent terms of energy and contributions to their
930 sig=sig0ij*dsqrt(sigsq)
931 rij_shift=1.0D0/rij-sig+sig0ij
932 c for diagnostics; uncomment
933 c rij_shift=1.2*sig0ij
934 C I hate to put IF's in the loops, but here don't have another choice!!!!
935 if (rij_shift.le.0.0D0) then
937 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
938 cd & restyp(itypi),i,restyp(itypj),j,
939 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
943 c---------------------------------------------------------------
944 rij_shift=1.0D0/rij_shift
946 e1=fac*fac*aa(itypi,itypj)
947 e2=fac*bb(itypi,itypj)
948 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
949 eps2der=evdwij*eps3rt
950 eps3der=evdwij*eps2rt
951 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
952 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
953 evdwij=evdwij*eps2rt*eps3rt
956 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
957 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
958 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
959 & restyp(itypi),i,restyp(itypj),j,
960 & epsi,sigm,chi1,chi2,chip1,chip2,
961 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
962 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
966 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
969 C Calculate gradient components.
970 e1=e1*eps1*eps2rt**2*eps3rt**2
971 fac=-expon*(e1+evdwij)*rij_shift
974 fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij
976 C Calculate the radial part of the gradient
980 C Calculate angular part of the gradient.
981 call sc_grad_scale(sss)
986 c write (iout,*) "Number of loop steps in EGB:",ind
987 cccc energy_dec=.false.
990 C-----------------------------------------------------------------------------
991 subroutine egbv_long(evdw)
993 C This subroutine calculates the interaction energy of nonbonded side chains
994 C assuming the Gay-Berne-Vorobjev potential of interaction.
996 implicit real*8 (a-h,o-z)
1000 include 'COMMON.LOCAL'
1001 include 'COMMON.CHAIN'
1002 include 'COMMON.DERIV'
1003 include 'COMMON.NAMES'
1004 include 'COMMON.INTERACT'
1005 include 'COMMON.IOUNITS'
1006 include 'COMMON.CALC'
1007 common /srutu/ icall
1010 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1013 c if (icall.eq.0) lprn=.true.
1015 do i=iatsc_s,iatsc_e
1017 if (itypi.eq.ntyp1) cycle
1022 dxi=dc_norm(1,nres+i)
1023 dyi=dc_norm(2,nres+i)
1024 dzi=dc_norm(3,nres+i)
1025 c dsci_inv=dsc_inv(itypi)
1026 dsci_inv=vbld_inv(i+nres)
1028 C Calculate SC interaction energy.
1030 do iint=1,nint_gr(i)
1031 do j=istart(i,iint),iend(i,iint)
1034 if (itypj.eq.ntyp1) cycle
1035 c dscj_inv=dsc_inv(itypj)
1036 dscj_inv=vbld_inv(j+nres)
1037 sig0ij=sigma(itypi,itypj)
1038 r0ij=r0(itypi,itypj)
1039 chi1=chi(itypi,itypj)
1040 chi2=chi(itypj,itypi)
1047 alf12=0.5D0*(alf1+alf2)
1051 dxj=dc_norm(1,nres+j)
1052 dyj=dc_norm(2,nres+j)
1053 dzj=dc_norm(3,nres+j)
1054 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1057 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1059 if (sss.lt.1.0d0) then
1061 C Calculate angle-dependent terms of energy and contributions to their
1065 sig=sig0ij*dsqrt(sigsq)
1066 rij_shift=1.0D0/rij-sig+r0ij
1067 C I hate to put IF's in the loops, but here don't have another choice!!!!
1068 if (rij_shift.le.0.0D0) then
1073 c---------------------------------------------------------------
1074 rij_shift=1.0D0/rij_shift
1075 fac=rij_shift**expon
1076 e1=fac*fac*aa(itypi,itypj)
1077 e2=fac*bb(itypi,itypj)
1078 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1079 eps2der=evdwij*eps3rt
1080 eps3der=evdwij*eps2rt
1081 fac_augm=rrij**expon
1082 e_augm=augm(itypi,itypj)*fac_augm
1083 evdwij=evdwij*eps2rt*eps3rt
1084 evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
1086 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1087 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1088 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1089 & restyp(itypi),i,restyp(itypj),j,
1090 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1091 & chi1,chi2,chip1,chip2,
1092 & eps1,eps2rt**2,eps3rt**2,
1093 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1096 C Calculate gradient components.
1097 e1=e1*eps1*eps2rt**2*eps3rt**2
1098 fac=-expon*(e1+evdwij)*rij_shift
1100 fac=rij*fac-2*expon*rrij*e_augm
1101 C Calculate the radial part of the gradient
1105 C Calculate angular part of the gradient.
1106 call sc_grad_scale(1.0d0-sss)
1112 C-----------------------------------------------------------------------------
1113 subroutine egbv_short(evdw)
1115 C This subroutine calculates the interaction energy of nonbonded side chains
1116 C assuming the Gay-Berne-Vorobjev potential of interaction.
1118 implicit real*8 (a-h,o-z)
1119 include 'DIMENSIONS'
1120 include 'COMMON.GEO'
1121 include 'COMMON.VAR'
1122 include 'COMMON.LOCAL'
1123 include 'COMMON.CHAIN'
1124 include 'COMMON.DERIV'
1125 include 'COMMON.NAMES'
1126 include 'COMMON.INTERACT'
1127 include 'COMMON.IOUNITS'
1128 include 'COMMON.CALC'
1129 common /srutu/ icall
1132 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1135 c if (icall.eq.0) lprn=.true.
1137 do i=iatsc_s,iatsc_e
1139 if (itypi.eq.ntyp1) cycle
1144 dxi=dc_norm(1,nres+i)
1145 dyi=dc_norm(2,nres+i)
1146 dzi=dc_norm(3,nres+i)
1147 c dsci_inv=dsc_inv(itypi)
1148 dsci_inv=vbld_inv(i+nres)
1150 C Calculate SC interaction energy.
1152 do iint=1,nint_gr(i)
1153 do j=istart(i,iint),iend(i,iint)
1156 if (itypj.eq.ntyp1) cycle
1157 c dscj_inv=dsc_inv(itypj)
1158 dscj_inv=vbld_inv(j+nres)
1159 sig0ij=sigma(itypi,itypj)
1160 r0ij=r0(itypi,itypj)
1161 chi1=chi(itypi,itypj)
1162 chi2=chi(itypj,itypi)
1169 alf12=0.5D0*(alf1+alf2)
1173 dxj=dc_norm(1,nres+j)
1174 dyj=dc_norm(2,nres+j)
1175 dzj=dc_norm(3,nres+j)
1176 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1179 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1181 if (sss.gt.0.0d0) then
1183 C Calculate angle-dependent terms of energy and contributions to their
1187 sig=sig0ij*dsqrt(sigsq)
1188 rij_shift=1.0D0/rij-sig+r0ij
1189 C I hate to put IF's in the loops, but here don't have another choice!!!!
1190 if (rij_shift.le.0.0D0) then
1195 c---------------------------------------------------------------
1196 rij_shift=1.0D0/rij_shift
1197 fac=rij_shift**expon
1198 e1=fac*fac*aa(itypi,itypj)
1199 e2=fac*bb(itypi,itypj)
1200 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1201 eps2der=evdwij*eps3rt
1202 eps3der=evdwij*eps2rt
1203 fac_augm=rrij**expon
1204 e_augm=augm(itypi,itypj)*fac_augm
1205 evdwij=evdwij*eps2rt*eps3rt
1206 evdw=evdw+(evdwij+e_augm)*sss
1208 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
1209 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
1210 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1211 & restyp(itypi),i,restyp(itypj),j,
1212 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1213 & chi1,chi2,chip1,chip2,
1214 & eps1,eps2rt**2,eps3rt**2,
1215 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1218 C Calculate gradient components.
1219 e1=e1*eps1*eps2rt**2*eps3rt**2
1220 fac=-expon*(e1+evdwij)*rij_shift
1222 fac=rij*fac-2*expon*rrij*e_augm
1223 C Calculate the radial part of the gradient
1227 C Calculate angular part of the gradient.
1228 call sc_grad_scale(sss)
1234 C----------------------------------------------------------------------------
1235 subroutine sc_grad_scale(scalfac)
1236 implicit real*8 (a-h,o-z)
1237 include 'DIMENSIONS'
1238 include 'COMMON.CHAIN'
1239 include 'COMMON.DERIV'
1240 include 'COMMON.CALC'
1241 include 'COMMON.IOUNITS'
1242 double precision dcosom1(3),dcosom2(3)
1243 double precision scalfac
1244 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1245 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1246 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1247 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1251 c eom12=evdwij*eps1_om12
1253 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1254 c & " sigder",sigder
1255 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1256 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1258 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1259 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1262 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1264 c write (iout,*) "gg",(gg(k),k=1,3)
1266 gvdwx(k,i)=gvdwx(k,i)-gg(k)
1267 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1268 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1269 gvdwx(k,j)=gvdwx(k,j)+gg(k)
1270 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1271 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1272 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1273 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1274 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1275 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1278 C Calculate the components of the gradient in DC and X
1281 gvdwc(l,i)=gvdwc(l,i)-gg(l)
1282 gvdwc(l,j)=gvdwc(l,j)+gg(l)
1286 C--------------------------------------------------------------------------
1287 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1289 C This subroutine calculates the average interaction energy and its gradient
1290 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1291 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1292 C The potential depends both on the distance of peptide-group centers and on
1293 C the orientation of the CA-CA virtual bonds.
1295 implicit real*8 (a-h,o-z)
1299 include 'DIMENSIONS'
1300 include 'COMMON.CONTROL'
1301 include 'COMMON.SETUP'
1302 include 'COMMON.IOUNITS'
1303 include 'COMMON.GEO'
1304 include 'COMMON.VAR'
1305 include 'COMMON.LOCAL'
1306 include 'COMMON.CHAIN'
1307 include 'COMMON.DERIV'
1308 include 'COMMON.INTERACT'
1309 include 'COMMON.CONTACTS'
1310 include 'COMMON.TORSION'
1311 include 'COMMON.VECTORS'
1312 include 'COMMON.FFIELD'
1313 include 'COMMON.TIME1'
1314 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1315 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1316 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1317 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1318 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1319 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1321 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1323 double precision scal_el /1.0d0/
1325 double precision scal_el /0.5d0/
1328 C 13-go grudnia roku pamietnego...
1329 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1330 & 0.0d0,1.0d0,0.0d0,
1331 & 0.0d0,0.0d0,1.0d0/
1332 cd write(iout,*) 'In EELEC'
1334 cd write(iout,*) 'Type',i
1335 cd write(iout,*) 'B1',B1(:,i)
1336 cd write(iout,*) 'B2',B2(:,i)
1337 cd write(iout,*) 'CC',CC(:,:,i)
1338 cd write(iout,*) 'DD',DD(:,:,i)
1339 cd write(iout,*) 'EE',EE(:,:,i)
1341 cd call check_vecgrad
1343 C print *,"WCHODZE3"
1344 if (icheckgrad.eq.1) then
1346 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1348 dc_norm(k,i)=dc(k,i)*fac
1350 c write (iout,*) 'i',i,' fac',fac
1353 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1354 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1355 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1356 c call vec_and_deriv
1362 time_mat=time_mat+MPI_Wtime()-time01
1366 cd write (iout,*) 'i=',i
1368 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1371 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1372 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1385 cd print '(a)','Enter EELEC'
1386 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1388 gel_loc_loc(i)=0.0d0
1393 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1395 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1397 do i=iturn3_start,iturn3_end
1398 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1399 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1400 & .or. itype(i-1).eq.ntyp1
1401 & .or. itype(i+4).eq.ntyp1
1406 dx_normi=dc_norm(1,i)
1407 dy_normi=dc_norm(2,i)
1408 dz_normi=dc_norm(3,i)
1409 xmedi=c(1,i)+0.5d0*dxi
1410 ymedi=c(2,i)+0.5d0*dyi
1411 zmedi=c(3,i)+0.5d0*dzi
1412 xmedi=mod(xmedi,boxxsize)
1413 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1414 ymedi=mod(ymedi,boxysize)
1415 if (ymedi.lt.0) ymedi=ymedi+boxysize
1416 zmedi=mod(zmedi,boxzsize)
1417 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1419 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1420 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1421 num_cont_hb(i)=num_conti
1423 do i=iturn4_start,iturn4_end
1424 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1425 & .or. itype(i+3).eq.ntyp1
1426 & .or. itype(i+4).eq.ntyp1
1427 & .or. itype(i+5).eq.ntyp1
1428 & .or. itype(i-1).eq.ntyp1
1433 dx_normi=dc_norm(1,i)
1434 dy_normi=dc_norm(2,i)
1435 dz_normi=dc_norm(3,i)
1436 xmedi=c(1,i)+0.5d0*dxi
1437 ymedi=c(2,i)+0.5d0*dyi
1438 zmedi=c(3,i)+0.5d0*dzi
1439 xmedi=mod(xmedi,boxxsize)
1440 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1441 ymedi=mod(ymedi,boxysize)
1442 if (ymedi.lt.0) ymedi=ymedi+boxysize
1443 zmedi=mod(zmedi,boxzsize)
1444 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1445 num_conti=num_cont_hb(i)
1446 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1447 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1448 & call eturn4(i,eello_turn4)
1449 num_cont_hb(i)=num_conti
1452 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1454 do i=iatel_s,iatel_e
1455 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1456 & .or. itype(i+2).eq.ntyp1
1457 & .or. itype(i-1).eq.ntyp1
1462 dx_normi=dc_norm(1,i)
1463 dy_normi=dc_norm(2,i)
1464 dz_normi=dc_norm(3,i)
1465 xmedi=c(1,i)+0.5d0*dxi
1466 ymedi=c(2,i)+0.5d0*dyi
1467 zmedi=c(3,i)+0.5d0*dzi
1468 xmedi=mod(xmedi,boxxsize)
1469 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1470 ymedi=mod(ymedi,boxysize)
1471 if (ymedi.lt.0) ymedi=ymedi+boxysize
1472 zmedi=mod(zmedi,boxzsize)
1473 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1474 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1475 num_conti=num_cont_hb(i)
1476 do j=ielstart(i),ielend(i)
1477 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1478 & .or.itype(j+2).eq.ntyp1
1479 & .or.itype(j-1).eq.ntyp1
1481 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1483 num_cont_hb(i)=num_conti
1485 c write (iout,*) "Number of loop steps in EELEC:",ind
1487 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1488 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1490 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1491 ccc eel_loc=eel_loc+eello_turn3
1492 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1495 C-------------------------------------------------------------------------------
1496 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1497 implicit real*8 (a-h,o-z)
1498 include 'DIMENSIONS'
1502 include 'COMMON.CONTROL'
1503 include 'COMMON.IOUNITS'
1504 include 'COMMON.GEO'
1505 include 'COMMON.VAR'
1506 include 'COMMON.LOCAL'
1507 include 'COMMON.CHAIN'
1508 include 'COMMON.DERIV'
1509 include 'COMMON.INTERACT'
1510 include 'COMMON.CONTACTS'
1511 include 'COMMON.TORSION'
1512 include 'COMMON.VECTORS'
1513 include 'COMMON.FFIELD'
1514 include 'COMMON.TIME1'
1515 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1516 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1517 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1518 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1519 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1520 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1522 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1524 double precision scal_el /1.0d0/
1526 double precision scal_el /0.5d0/
1529 C 13-go grudnia roku pamietnego...
1530 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1531 & 0.0d0,1.0d0,0.0d0,
1532 & 0.0d0,0.0d0,1.0d0/
1533 c time00=MPI_Wtime()
1534 cd write (iout,*) "eelecij",i,j
1535 C print *,"WCHODZE2"
1539 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1540 aaa=app(iteli,itelj)
1541 bbb=bpp(iteli,itelj)
1542 ael6i=ael6(iteli,itelj)
1543 ael3i=ael3(iteli,itelj)
1547 dx_normj=dc_norm(1,j)
1548 dy_normj=dc_norm(2,j)
1549 dz_normj=dc_norm(3,j)
1554 if (xj.lt.0) xj=xj+boxxsize
1556 if (yj.lt.0) yj=yj+boxysize
1558 if (zj.lt.0) zj=zj+boxzsize
1559 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1567 xj=xj_safe+xshift*boxxsize
1568 yj=yj_safe+yshift*boxysize
1569 zj=zj_safe+zshift*boxzsize
1570 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1571 if(dist_temp.lt.dist_init) then
1581 if (isubchap.eq.1) then
1591 rij=xj*xj+yj*yj+zj*zj
1595 c For extracting the short-range part of Evdwpp
1596 sss=sscale(rij/rpp(iteli,itelj))
1597 sssgrad=sscagrad(rij/rpp(iteli,itelj))
1600 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1601 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1602 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1603 fac=cosa-3.0D0*cosb*cosg
1605 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1606 if (j.eq.i+2) ev1=scal_el*ev1
1611 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1614 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1615 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1617 evdw1=evdw1+evdwij*(1.0d0-sss)
1618 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1619 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1620 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1621 cd & xmedi,ymedi,zmedi,xj,yj,zj
1623 if (energy_dec) then
1624 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1625 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1629 C Calculate contributions to the Cartesian gradient.
1632 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1633 facel=-3*rrmij*(el1+eesij)
1639 * Radial derivatives. First process both termini of the fragment (i,j)
1645 c ghalf=0.5D0*ggg(k)
1646 c gelc(k,i)=gelc(k,i)+ghalf
1647 c gelc(k,j)=gelc(k,j)+ghalf
1649 c 9/28/08 AL Gradient compotents will be summed only at the end
1651 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1652 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1655 * Loop over residues i+1 thru j-1.
1659 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1662 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1663 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1664 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1666 c ghalf=0.5D0*ggg(k)
1667 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1668 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1670 c 9/28/08 AL Gradient compotents will be summed only at the end
1672 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1673 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1676 * Loop over residues i+1 thru j-1.
1680 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1684 facvdw=ev1+evdwij*(1.0d0-sss)
1687 fac=-3*rrmij*(facvdw+facvdw+facel)
1692 * Radial derivatives. First process both termini of the fragment (i,j)
1698 c ghalf=0.5D0*ggg(k)
1699 c gelc(k,i)=gelc(k,i)+ghalf
1700 c gelc(k,j)=gelc(k,j)+ghalf
1702 c 9/28/08 AL Gradient compotents will be summed only at the end
1704 gelc_long(k,j)=gelc(k,j)+ggg(k)
1705 gelc_long(k,i)=gelc(k,i)-ggg(k)
1708 * Loop over residues i+1 thru j-1.
1712 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1715 c 9/28/08 AL Gradient compotents will be summed only at the end
1719 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1720 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1721 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1723 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1724 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1730 ecosa=2.0D0*fac3*fac1+fac4
1733 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1734 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1736 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1737 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1739 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1740 cd & (dcosg(k),k=1,3)
1742 ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1745 c ghalf=0.5D0*ggg(k)
1746 c gelc(k,i)=gelc(k,i)+ghalf
1747 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1748 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1749 c gelc(k,j)=gelc(k,j)+ghalf
1750 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1751 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1755 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1760 & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1761 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1763 & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1764 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1765 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1766 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1768 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1769 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1770 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1772 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1773 C energy of a peptide unit is assumed in the form of a second-order
1774 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1775 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1776 C are computed for EVERY pair of non-contiguous peptide groups.
1778 if (j.lt.nres-1) then
1789 muij(kkk)=mu(k,i)*mu(l,j)
1792 cd write (iout,*) 'EELEC: i',i,' j',j
1793 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1794 cd write(iout,*) 'muij',muij
1795 ury=scalar(uy(1,i),erij)
1796 urz=scalar(uz(1,i),erij)
1797 vry=scalar(uy(1,j),erij)
1798 vrz=scalar(uz(1,j),erij)
1799 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1800 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1801 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1802 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1803 fac=dsqrt(-ael6i)*r3ij
1808 cd write (iout,'(4i5,4f10.5)')
1809 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1810 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1811 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1812 cd & uy(:,j),uz(:,j)
1813 cd write (iout,'(4f10.5)')
1814 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1815 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1816 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
1817 cd write (iout,'(9f10.5/)')
1818 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1819 C Derivatives of the elements of A in virtual-bond vectors
1820 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1822 uryg(k,1)=scalar(erder(1,k),uy(1,i))
1823 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1824 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1825 urzg(k,1)=scalar(erder(1,k),uz(1,i))
1826 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1827 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1828 vryg(k,1)=scalar(erder(1,k),uy(1,j))
1829 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1830 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1831 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1832 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1833 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1835 C Compute radial contributions to the gradient
1853 C Add the contributions coming from er
1856 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1857 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1858 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1859 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1862 C Derivatives in DC(i)
1863 cgrad ghalf1=0.5d0*agg(k,1)
1864 cgrad ghalf2=0.5d0*agg(k,2)
1865 cgrad ghalf3=0.5d0*agg(k,3)
1866 cgrad ghalf4=0.5d0*agg(k,4)
1867 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1868 & -3.0d0*uryg(k,2)*vry)!+ghalf1
1869 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1870 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
1871 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1872 & -3.0d0*urzg(k,2)*vry)!+ghalf3
1873 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1874 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
1875 C Derivatives in DC(i+1)
1876 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
1877 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
1878 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
1879 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
1880 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
1881 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
1882 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
1883 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
1884 C Derivatives in DC(j)
1885 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
1886 & -3.0d0*vryg(k,2)*ury)!+ghalf1
1887 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
1888 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
1889 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
1890 & -3.0d0*vryg(k,2)*urz)!+ghalf3
1891 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
1892 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
1893 C Derivatives in DC(j+1) or DC(nres-1)
1894 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
1895 & -3.0d0*vryg(k,3)*ury)
1896 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
1897 & -3.0d0*vrzg(k,3)*ury)
1898 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
1899 & -3.0d0*vryg(k,3)*urz)
1900 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
1901 & -3.0d0*vrzg(k,3)*urz)
1902 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
1904 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
1917 aggi(k,l)=-aggi(k,l)
1918 aggi1(k,l)=-aggi1(k,l)
1919 aggj(k,l)=-aggj(k,l)
1920 aggj1(k,l)=-aggj1(k,l)
1923 if (j.lt.nres-1) then
1929 aggi(k,l)=-aggi(k,l)
1930 aggi1(k,l)=-aggi1(k,l)
1931 aggj(k,l)=-aggj(k,l)
1932 aggj1(k,l)=-aggj1(k,l)
1943 aggi(k,l)=-aggi(k,l)
1944 aggi1(k,l)=-aggi1(k,l)
1945 aggj(k,l)=-aggj(k,l)
1946 aggj1(k,l)=-aggj1(k,l)
1951 IF (wel_loc.gt.0.0d0) THEN
1952 C Contribution to the local-electrostatic energy coming from the i-j pair
1953 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
1955 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
1957 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1958 & 'eelloc',i,j,eel_loc_ij
1960 eel_loc=eel_loc+eel_loc_ij
1961 C Partial derivatives in virtual-bond dihedral angles gamma
1963 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
1964 & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1965 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1966 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
1967 & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1968 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1969 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
1971 ggg(l)=agg(l,1)*muij(1)+
1972 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1973 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
1974 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
1975 cgrad ghalf=0.5d0*ggg(l)
1976 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
1977 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
1981 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
1984 C Remaining derivatives of eello
1986 gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1987 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1988 gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1989 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1990 gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1991 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1992 gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1993 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1996 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
1997 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
1998 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
1999 & .and. num_conti.le.maxconts) then
2000 c write (iout,*) i,j," entered corr"
2002 C Calculate the contact function. The ith column of the array JCONT will
2003 C contain the numbers of atoms that make contacts with the atom I (of numbers
2004 C greater than I). The arrays FACONT and GACONT will contain the values of
2005 C the contact function and its derivative.
2006 c r0ij=1.02D0*rpp(iteli,itelj)
2007 c r0ij=1.11D0*rpp(iteli,itelj)
2008 r0ij=2.20D0*rpp(iteli,itelj)
2009 c r0ij=1.55D0*rpp(iteli,itelj)
2010 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2011 if (fcont.gt.0.0D0) then
2012 num_conti=num_conti+1
2013 if (num_conti.gt.maxconts) then
2014 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2015 & ' will skip next contacts for this conf.'
2017 jcont_hb(num_conti,i)=j
2018 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
2019 cd & " jcont_hb",jcont_hb(num_conti,i)
2020 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2021 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2022 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2024 d_cont(num_conti,i)=rij
2025 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2026 C --- Electrostatic-interaction matrix ---
2027 a_chuj(1,1,num_conti,i)=a22
2028 a_chuj(1,2,num_conti,i)=a23
2029 a_chuj(2,1,num_conti,i)=a32
2030 a_chuj(2,2,num_conti,i)=a33
2031 C --- Gradient of rij
2033 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2040 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2041 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2042 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2043 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2044 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2049 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2050 C Calculate contact energies
2052 wij=cosa-3.0D0*cosb*cosg
2055 c fac3=dsqrt(-ael6i)/r0ij**3
2056 fac3=dsqrt(-ael6i)*r3ij
2057 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2058 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2059 if (ees0tmp.gt.0) then
2060 ees0pij=dsqrt(ees0tmp)
2064 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2065 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2066 if (ees0tmp.gt.0) then
2067 ees0mij=dsqrt(ees0tmp)
2072 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2073 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2074 C Diagnostics. Comment out or remove after debugging!
2075 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2076 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2077 c ees0m(num_conti,i)=0.0D0
2079 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2080 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2081 C Angular derivatives of the contact function
2082 ees0pij1=fac3/ees0pij
2083 ees0mij1=fac3/ees0mij
2084 fac3p=-3.0D0*fac3*rrmij
2085 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2086 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2088 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2089 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2090 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2091 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2092 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2093 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2094 ecosap=ecosa1+ecosa2
2095 ecosbp=ecosb1+ecosb2
2096 ecosgp=ecosg1+ecosg2
2097 ecosam=ecosa1-ecosa2
2098 ecosbm=ecosb1-ecosb2
2099 ecosgm=ecosg1-ecosg2
2108 facont_hb(num_conti,i)=fcont
2109 fprimcont=fprimcont/rij
2110 cd facont_hb(num_conti,i)=1.0D0
2111 C Following line is for diagnostics.
2114 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2115 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2118 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2119 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2121 gggp(1)=gggp(1)+ees0pijp*xj
2122 gggp(2)=gggp(2)+ees0pijp*yj
2123 gggp(3)=gggp(3)+ees0pijp*zj
2124 gggm(1)=gggm(1)+ees0mijp*xj
2125 gggm(2)=gggm(2)+ees0mijp*yj
2126 gggm(3)=gggm(3)+ees0mijp*zj
2127 C Derivatives due to the contact function
2128 gacont_hbr(1,num_conti,i)=fprimcont*xj
2129 gacont_hbr(2,num_conti,i)=fprimcont*yj
2130 gacont_hbr(3,num_conti,i)=fprimcont*zj
2133 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2134 c following the change of gradient-summation algorithm.
2136 cgrad ghalfp=0.5D0*gggp(k)
2137 cgrad ghalfm=0.5D0*gggm(k)
2138 gacontp_hb1(k,num_conti,i)=!ghalfp
2139 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2140 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2141 gacontp_hb2(k,num_conti,i)=!ghalfp
2142 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2143 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2144 gacontp_hb3(k,num_conti,i)=gggp(k)
2145 gacontm_hb1(k,num_conti,i)=!ghalfm
2146 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2147 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2148 gacontm_hb2(k,num_conti,i)=!ghalfm
2149 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2150 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2151 gacontm_hb3(k,num_conti,i)=gggm(k)
2154 endif ! num_conti.le.maxconts
2157 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2160 ghalf=0.5d0*agg(l,k)
2161 aggi(l,k)=aggi(l,k)+ghalf
2162 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2163 aggj(l,k)=aggj(l,k)+ghalf
2166 if (j.eq.nres-1 .and. i.lt.j-2) then
2169 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2174 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2177 C-----------------------------------------------------------------------
2178 subroutine evdwpp_short(evdw1)
2182 implicit real*8 (a-h,o-z)
2183 include 'DIMENSIONS'
2184 include 'COMMON.CONTROL'
2185 include 'COMMON.IOUNITS'
2186 include 'COMMON.GEO'
2187 include 'COMMON.VAR'
2188 include 'COMMON.LOCAL'
2189 include 'COMMON.CHAIN'
2190 include 'COMMON.DERIV'
2191 include 'COMMON.INTERACT'
2192 include 'COMMON.CONTACTS'
2193 include 'COMMON.TORSION'
2194 include 'COMMON.VECTORS'
2195 include 'COMMON.FFIELD'
2197 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2199 double precision scal_el /1.0d0/
2201 double precision scal_el /0.5d0/
2205 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2206 c & " iatel_e_vdw",iatel_e_vdw
2208 do i=iatel_s_vdw,iatel_e_vdw
2209 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2213 dx_normi=dc_norm(1,i)
2214 dy_normi=dc_norm(2,i)
2215 dz_normi=dc_norm(3,i)
2216 xmedi=c(1,i)+0.5d0*dxi
2217 ymedi=c(2,i)+0.5d0*dyi
2218 zmedi=c(3,i)+0.5d0*dzi
2219 xmedi=mod(xmedi,boxxsize)
2220 if (xmedi.lt.0) xmedi=xmedi+boxxsize
2221 ymedi=mod(ymedi,boxysize)
2222 if (ymedi.lt.0) ymedi=ymedi+boxysize
2223 zmedi=mod(zmedi,boxzsize)
2224 if (zmedi.lt.0) zmedi=zmedi+boxzsize
2226 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2227 c & ' ielend',ielend_vdw(i)
2229 do j=ielstart_vdw(i),ielend_vdw(i)
2230 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2234 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2235 aaa=app(iteli,itelj)
2236 bbb=bpp(iteli,itelj)
2240 dx_normj=dc_norm(1,j)
2241 dy_normj=dc_norm(2,j)
2242 dz_normj=dc_norm(3,j)
2247 if (xj.lt.0) xj=xj+boxxsize
2249 if (yj.lt.0) yj=yj+boxysize
2251 if (zj.lt.0) zj=zj+boxzsize
2252 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2260 xj=xj_safe+xshift*boxxsize
2261 yj=yj_safe+yshift*boxysize
2262 zj=zj_safe+zshift*boxzsize
2263 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2264 if(dist_temp.lt.dist_init) then
2274 if (isubchap.eq.1) then
2283 rij=xj*xj+yj*yj+zj*zj
2286 sss=sscale(rij/rpp(iteli,itelj))
2287 sssgrad=sscagrad(rij/rpp(iteli,itelj))
2288 if (sss.gt.0.0d0) then
2293 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2294 if (j.eq.i+2) ev1=scal_el*ev1
2297 if (energy_dec) then
2298 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2300 evdw1=evdw1+evdwij*sss
2302 C Calculate contributions to the Cartesian gradient.
2304 facvdw=-6*rrmij*(ev1+evdwij)*sss
2305 ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2306 ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2307 ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2312 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2313 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2320 C-----------------------------------------------------------------------------
2321 subroutine escp_long(evdw2,evdw2_14)
2323 C This subroutine calculates the excluded-volume interaction energy between
2324 C peptide-group centers and side chains and its gradient in virtual-bond and
2325 C side-chain vectors.
2327 implicit real*8 (a-h,o-z)
2328 include 'DIMENSIONS'
2329 include 'COMMON.GEO'
2330 include 'COMMON.VAR'
2331 include 'COMMON.LOCAL'
2332 include 'COMMON.CHAIN'
2333 include 'COMMON.DERIV'
2334 include 'COMMON.INTERACT'
2335 include 'COMMON.FFIELD'
2336 include 'COMMON.IOUNITS'
2337 include 'COMMON.CONTROL'
2341 CD print '(a)','Enter ESCP KURWA'
2342 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2343 do i=iatscp_s,iatscp_e
2344 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2346 xi=0.5D0*(c(1,i)+c(1,i+1))
2347 yi=0.5D0*(c(2,i)+c(2,i+1))
2348 zi=0.5D0*(c(3,i)+c(3,i+1))
2350 if (xi.lt.0) xi=xi+boxxsize
2352 if (yi.lt.0) yi=yi+boxysize
2354 if (zi.lt.0) zi=zi+boxzsize
2355 do iint=1,nscp_gr(i)
2357 do j=iscpstart(i,iint),iscpend(i,iint)
2359 if (itypj.eq.ntyp1) cycle
2360 C Uncomment following three lines for SC-p interactions
2364 C Uncomment following three lines for Ca-p interactions
2368 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2376 xj=xj_safe+xshift*boxxsize
2377 yj=yj_safe+yshift*boxysize
2378 zj=zj_safe+zshift*boxzsize
2379 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2380 if(dist_temp.lt.dist_init) then
2390 if (subchap.eq.1) then
2400 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2402 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2403 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2404 if (sss.lt.1.0d0) then
2406 e1=fac*fac*aad(itypj,iteli)
2407 e2=fac*bad(itypj,iteli)
2408 if (iabs(j-i) .le. 2) then
2411 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2414 evdw2=evdw2+evdwij*(1.0d0-sss)
2415 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2416 & 'evdw2',i,j,sss,evdwij
2418 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2421 fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2422 fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2426 C Uncomment following three lines for SC-p interactions
2428 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2430 C Uncomment following line for SC-p interactions
2431 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2433 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2434 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2443 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2444 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2445 gradx_scp(j,i)=expon*gradx_scp(j,i)
2448 C******************************************************************************
2452 C To save time the factor EXPON has been extracted from ALL components
2453 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2456 C******************************************************************************
2459 C-----------------------------------------------------------------------------
2460 subroutine escp_short(evdw2,evdw2_14)
2462 C This subroutine calculates the excluded-volume interaction energy between
2463 C peptide-group centers and side chains and its gradient in virtual-bond and
2464 C side-chain vectors.
2466 implicit real*8 (a-h,o-z)
2467 include 'DIMENSIONS'
2468 include 'COMMON.GEO'
2469 include 'COMMON.VAR'
2470 include 'COMMON.LOCAL'
2471 include 'COMMON.CHAIN'
2472 include 'COMMON.DERIV'
2473 include 'COMMON.INTERACT'
2474 include 'COMMON.FFIELD'
2475 include 'COMMON.IOUNITS'
2476 include 'COMMON.CONTROL'
2480 cd print '(a)','Enter ESCP'
2481 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2482 do i=iatscp_s,iatscp_e
2483 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2485 xi=0.5D0*(c(1,i)+c(1,i+1))
2486 yi=0.5D0*(c(2,i)+c(2,i+1))
2487 zi=0.5D0*(c(3,i)+c(3,i+1))
2489 if (xi.lt.0) xi=xi+boxxsize
2491 if (yi.lt.0) yi=yi+boxysize
2493 if (zi.lt.0) zi=zi+boxzsize
2495 do iint=1,nscp_gr(i)
2497 do j=iscpstart(i,iint),iscpend(i,iint)
2499 if (itypj.eq.ntyp1) cycle
2500 C Uncomment following three lines for SC-p interactions
2504 C Uncomment following three lines for Ca-p interactions
2508 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2516 xj=xj_safe+xshift*boxxsize
2517 yj=yj_safe+yshift*boxysize
2518 zj=zj_safe+zshift*boxzsize
2519 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2520 if(dist_temp.lt.dist_init) then
2530 if (subchap.eq.1) then
2539 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2540 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2541 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2542 if (sss.gt.0.0d0) then
2545 e1=fac*fac*aad(itypj,iteli)
2546 e2=fac*bad(itypj,iteli)
2547 if (iabs(j-i) .le. 2) then
2550 evdw2_14=evdw2_14+(e1+e2)*sss
2553 evdw2=evdw2+evdwij*sss
2554 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2555 & 'evdw2',i,j,sss,evdwij
2557 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2559 fac=-(evdwij+e1)*rrij*sss
2560 fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2564 C Uncomment following three lines for SC-p interactions
2566 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2568 C Uncomment following line for SC-p interactions
2569 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2571 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2572 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2581 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2582 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2583 gradx_scp(j,i)=expon*gradx_scp(j,i)
2586 C******************************************************************************
2590 C To save time the factor EXPON has been extracted from ALL components
2591 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2594 C******************************************************************************