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)
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)
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
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
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
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/bb)**(1.0D0/6.0D0)
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
592 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
593 eps2der=evdwij*eps3rt
594 eps3der=evdwij*eps2rt
595 evdwij=evdwij*eps2rt*eps3rt
598 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
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 if ((zj.gt.bordlipbot)
705 &.and.(zj.lt.bordliptop)) then
706 C the energy transfer exist
707 if (zj.lt.buflipbot) then
708 C what fraction I am in
710 & ((positi-bordlipbot)/lipbufthick)
711 C lipbufthick is thickenes of lipid buffore
712 sslipj=sscalelip(fracinbuf)
713 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
714 elseif (zi.gt.bufliptop) then
715 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
716 sslipj=sscalelip(fracinbuf)
717 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
726 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
727 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
728 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
729 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
731 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
739 xj=xj_safe+xshift*boxxsize
740 yj=yj_safe+yshift*boxysize
741 zj=zj_safe+zshift*boxzsize
742 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
743 if(dist_temp.lt.dist_init) then
753 if (subchap.eq.1) then
762 dxj=dc_norm(1,nres+j)
763 dyj=dc_norm(2,nres+j)
764 dzj=dc_norm(3,nres+j)
765 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
767 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
768 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
769 if (sss.lt.1.0d0) then
771 C Calculate angle-dependent terms of energy and contributions to their
775 sig=sig0ij*dsqrt(sigsq)
776 rij_shift=1.0D0/rij-sig+sig0ij
777 c for diagnostics; uncomment
778 c rij_shift=1.2*sig0ij
779 C I hate to put IF's in the loops, but here don't have another choice!!!!
780 if (rij_shift.le.0.0D0) then
782 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
783 cd & restyp(itypi),i,restyp(itypj),j,
784 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
788 c---------------------------------------------------------------
789 rij_shift=1.0D0/rij_shift
793 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
794 eps2der=evdwij*eps3rt
795 eps3der=evdwij*eps2rt
796 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
797 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
798 evdwij=evdwij*eps2rt*eps3rt
799 evdw=evdw+evdwij*(1.0d0-sss)
801 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
803 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
804 & restyp(itypi),i,restyp(itypj),j,
805 & epsi,sigm,chi1,chi2,chip1,chip2,
806 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
807 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
811 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
814 C Calculate gradient components.
815 e1=e1*eps1*eps2rt**2*eps3rt**2
816 fac=-expon*(e1+evdwij)*rij_shift
819 fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij
821 C Calculate the radial part of the gradient
825 gg_lipi(3)=ssgradlipi*evdwij
826 gg_lipj(3)=ssgradlipj*evdwij
827 C Calculate angular part of the gradient.
828 call sc_grad_scale(1.0d0-sss)
833 c write (iout,*) "Number of loop steps in EGB:",ind
834 cccc energy_dec=.false.
837 C-----------------------------------------------------------------------------
838 subroutine egb_short(evdw)
840 C This subroutine calculates the interaction energy of nonbonded side chains
841 C assuming the Gay-Berne potential of interaction.
843 implicit real*8 (a-h,o-z)
847 include 'COMMON.LOCAL'
848 include 'COMMON.CHAIN'
849 include 'COMMON.DERIV'
850 include 'COMMON.NAMES'
851 include 'COMMON.INTERACT'
852 include 'COMMON.IOUNITS'
853 include 'COMMON.CALC'
854 include 'COMMON.CONTROL'
857 ccccc energy_dec=.false.
858 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
861 c if (icall.eq.0) lprn=.false.
865 if (itypi.eq.ntyp1) cycle
871 if (xi.lt.0) xi=xi+boxxsize
873 if (yi.lt.0) yi=yi+boxysize
875 if (zi.lt.0) zi=zi+boxzsize
876 dxi=dc_norm(1,nres+i)
877 dyi=dc_norm(2,nres+i)
878 dzi=dc_norm(3,nres+i)
879 c dsci_inv=dsc_inv(itypi)
880 dsci_inv=vbld_inv(i+nres)
881 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
882 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
884 C Calculate SC interaction energy.
887 do j=istart(i,iint),iend(i,iint)
890 if (itypj.eq.ntyp1) cycle
891 c dscj_inv=dsc_inv(itypj)
892 dscj_inv=vbld_inv(j+nres)
893 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
894 c & 1.0d0/vbld(j+nres)
895 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
896 sig0ij=sigma(itypi,itypj)
897 chi1=chi(itypi,itypj)
898 chi2=chi(itypj,itypi)
905 alf12=0.5D0*(alf1+alf2)
910 if (xj.lt.0) xj=xj+boxxsize
912 if (yj.lt.0) yj=yj+boxysize
914 if (zj.lt.0) zj=zj+boxzsize
915 if ((zj.gt.bordlipbot)
916 &.and.(zj.lt.bordliptop)) then
917 C the energy transfer exist
918 if (zj.lt.buflipbot) then
919 C what fraction I am in
921 & ((positi-bordlipbot)/lipbufthick)
922 C lipbufthick is thickenes of lipid buffore
923 sslipj=sscalelip(fracinbuf)
924 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
925 elseif (zi.gt.bufliptop) then
926 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
927 sslipj=sscalelip(fracinbuf)
928 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
937 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
938 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
939 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
940 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
941 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
949 xj=xj_safe+xshift*boxxsize
950 yj=yj_safe+yshift*boxysize
951 zj=zj_safe+zshift*boxzsize
952 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
953 if(dist_temp.lt.dist_init) then
963 if (subchap.eq.1) then
972 dxj=dc_norm(1,nres+j)
973 dyj=dc_norm(2,nres+j)
974 dzj=dc_norm(3,nres+j)
975 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
977 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
978 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
979 if (sss.gt.0.0d0) then
981 C Calculate angle-dependent terms of energy and contributions to their
985 sig=sig0ij*dsqrt(sigsq)
986 rij_shift=1.0D0/rij-sig+sig0ij
987 c for diagnostics; uncomment
988 c rij_shift=1.2*sig0ij
989 C I hate to put IF's in the loops, but here don't have another choice!!!!
990 if (rij_shift.le.0.0D0) then
992 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
993 cd & restyp(itypi),i,restyp(itypj),j,
994 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
998 c---------------------------------------------------------------
999 rij_shift=1.0D0/rij_shift
1000 fac=rij_shift**expon
1003 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1004 eps2der=evdwij*eps3rt
1005 eps3der=evdwij*eps2rt
1006 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1007 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1008 evdwij=evdwij*eps2rt*eps3rt
1009 evdw=evdw+evdwij*sss
1011 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1013 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1014 & restyp(itypi),i,restyp(itypj),j,
1015 & epsi,sigm,chi1,chi2,chip1,chip2,
1016 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1017 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1021 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1024 C Calculate gradient components.
1025 e1=e1*eps1*eps2rt**2*eps3rt**2
1026 fac=-expon*(e1+evdwij)*rij_shift
1029 fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij
1031 C Calculate the radial part of the gradient
1035 gg_lipi(3)=ssgradlipi*evdwij
1036 gg_lipj(3)=ssgradlipj*evdwij
1037 C Calculate angular part of the gradient.
1038 call sc_grad_scale(sss)
1043 c write (iout,*) "Number of loop steps in EGB:",ind
1044 cccc energy_dec=.false.
1047 C-----------------------------------------------------------------------------
1048 subroutine egbv_long(evdw)
1050 C This subroutine calculates the interaction energy of nonbonded side chains
1051 C assuming the Gay-Berne-Vorobjev potential of interaction.
1053 implicit real*8 (a-h,o-z)
1054 include 'DIMENSIONS'
1055 include 'COMMON.GEO'
1056 include 'COMMON.VAR'
1057 include 'COMMON.LOCAL'
1058 include 'COMMON.CHAIN'
1059 include 'COMMON.DERIV'
1060 include 'COMMON.NAMES'
1061 include 'COMMON.INTERACT'
1062 include 'COMMON.IOUNITS'
1063 include 'COMMON.CALC'
1064 common /srutu/ icall
1067 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1070 c if (icall.eq.0) lprn=.true.
1072 do i=iatsc_s,iatsc_e
1074 if (itypi.eq.ntyp1) cycle
1079 dxi=dc_norm(1,nres+i)
1080 dyi=dc_norm(2,nres+i)
1081 dzi=dc_norm(3,nres+i)
1082 c dsci_inv=dsc_inv(itypi)
1083 dsci_inv=vbld_inv(i+nres)
1085 C Calculate SC interaction energy.
1087 do iint=1,nint_gr(i)
1088 do j=istart(i,iint),iend(i,iint)
1091 if (itypj.eq.ntyp1) cycle
1092 c dscj_inv=dsc_inv(itypj)
1093 dscj_inv=vbld_inv(j+nres)
1094 sig0ij=sigma(itypi,itypj)
1095 r0ij=r0(itypi,itypj)
1096 chi1=chi(itypi,itypj)
1097 chi2=chi(itypj,itypi)
1104 alf12=0.5D0*(alf1+alf2)
1108 dxj=dc_norm(1,nres+j)
1109 dyj=dc_norm(2,nres+j)
1110 dzj=dc_norm(3,nres+j)
1111 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1114 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1116 if (sss.lt.1.0d0) then
1118 C Calculate angle-dependent terms of energy and contributions to their
1122 sig=sig0ij*dsqrt(sigsq)
1123 rij_shift=1.0D0/rij-sig+r0ij
1124 C I hate to put IF's in the loops, but here don't have another choice!!!!
1125 if (rij_shift.le.0.0D0) then
1130 c---------------------------------------------------------------
1131 rij_shift=1.0D0/rij_shift
1132 fac=rij_shift**expon
1135 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1136 eps2der=evdwij*eps3rt
1137 eps3der=evdwij*eps2rt
1138 fac_augm=rrij**expon
1139 e_augm=augm(itypi,itypj)*fac_augm
1140 evdwij=evdwij*eps2rt*eps3rt
1141 evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
1143 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1145 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1146 & restyp(itypi),i,restyp(itypj),j,
1147 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1148 & chi1,chi2,chip1,chip2,
1149 & eps1,eps2rt**2,eps3rt**2,
1150 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1153 C Calculate gradient components.
1154 e1=e1*eps1*eps2rt**2*eps3rt**2
1155 fac=-expon*(e1+evdwij)*rij_shift
1157 fac=rij*fac-2*expon*rrij*e_augm
1158 C Calculate the radial part of the gradient
1162 C Calculate angular part of the gradient.
1163 call sc_grad_scale(1.0d0-sss)
1169 C-----------------------------------------------------------------------------
1170 subroutine egbv_short(evdw)
1172 C This subroutine calculates the interaction energy of nonbonded side chains
1173 C assuming the Gay-Berne-Vorobjev potential of interaction.
1175 implicit real*8 (a-h,o-z)
1176 include 'DIMENSIONS'
1177 include 'COMMON.GEO'
1178 include 'COMMON.VAR'
1179 include 'COMMON.LOCAL'
1180 include 'COMMON.CHAIN'
1181 include 'COMMON.DERIV'
1182 include 'COMMON.NAMES'
1183 include 'COMMON.INTERACT'
1184 include 'COMMON.IOUNITS'
1185 include 'COMMON.CALC'
1186 common /srutu/ icall
1189 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1192 c if (icall.eq.0) lprn=.true.
1194 do i=iatsc_s,iatsc_e
1196 if (itypi.eq.ntyp1) cycle
1201 dxi=dc_norm(1,nres+i)
1202 dyi=dc_norm(2,nres+i)
1203 dzi=dc_norm(3,nres+i)
1204 c dsci_inv=dsc_inv(itypi)
1205 dsci_inv=vbld_inv(i+nres)
1207 C Calculate SC interaction energy.
1209 do iint=1,nint_gr(i)
1210 do j=istart(i,iint),iend(i,iint)
1213 if (itypj.eq.ntyp1) cycle
1214 c dscj_inv=dsc_inv(itypj)
1215 dscj_inv=vbld_inv(j+nres)
1216 sig0ij=sigma(itypi,itypj)
1217 r0ij=r0(itypi,itypj)
1218 chi1=chi(itypi,itypj)
1219 chi2=chi(itypj,itypi)
1226 alf12=0.5D0*(alf1+alf2)
1230 dxj=dc_norm(1,nres+j)
1231 dyj=dc_norm(2,nres+j)
1232 dzj=dc_norm(3,nres+j)
1233 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1236 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1238 if (sss.gt.0.0d0) then
1240 C Calculate angle-dependent terms of energy and contributions to their
1244 sig=sig0ij*dsqrt(sigsq)
1245 rij_shift=1.0D0/rij-sig+r0ij
1246 C I hate to put IF's in the loops, but here don't have another choice!!!!
1247 if (rij_shift.le.0.0D0) then
1252 c---------------------------------------------------------------
1253 rij_shift=1.0D0/rij_shift
1254 fac=rij_shift**expon
1257 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1258 eps2der=evdwij*eps3rt
1259 eps3der=evdwij*eps2rt
1260 fac_augm=rrij**expon
1261 e_augm=augm(itypi,itypj)*fac_augm
1262 evdwij=evdwij*eps2rt*eps3rt
1263 evdw=evdw+(evdwij+e_augm)*sss
1265 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1267 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1268 & restyp(itypi),i,restyp(itypj),j,
1269 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1270 & chi1,chi2,chip1,chip2,
1271 & eps1,eps2rt**2,eps3rt**2,
1272 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1275 C Calculate gradient components.
1276 e1=e1*eps1*eps2rt**2*eps3rt**2
1277 fac=-expon*(e1+evdwij)*rij_shift
1279 fac=rij*fac-2*expon*rrij*e_augm
1280 C Calculate the radial part of the gradient
1284 C Calculate angular part of the gradient.
1285 call sc_grad_scale(sss)
1291 C----------------------------------------------------------------------------
1292 subroutine sc_grad_scale(scalfac)
1293 implicit real*8 (a-h,o-z)
1294 include 'DIMENSIONS'
1295 include 'COMMON.CHAIN'
1296 include 'COMMON.DERIV'
1297 include 'COMMON.CALC'
1298 include 'COMMON.IOUNITS'
1299 double precision dcosom1(3),dcosom2(3)
1300 double precision scalfac
1301 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1302 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1303 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1304 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1308 c eom12=evdwij*eps1_om12
1310 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1311 c & " sigder",sigder
1312 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1313 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1315 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1316 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1319 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1321 c write (iout,*) "gg",(gg(k),k=1,3)
1323 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
1324 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1325 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1326 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
1327 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1328 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1329 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1330 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1331 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1332 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1335 C Calculate the components of the gradient in DC and X
1338 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
1339 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
1343 C--------------------------------------------------------------------------
1344 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1346 C This subroutine calculates the average interaction energy and its gradient
1347 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1348 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1349 C The potential depends both on the distance of peptide-group centers and on
1350 C the orientation of the CA-CA virtual bonds.
1352 implicit real*8 (a-h,o-z)
1356 include 'DIMENSIONS'
1357 include 'COMMON.CONTROL'
1358 include 'COMMON.SETUP'
1359 include 'COMMON.IOUNITS'
1360 include 'COMMON.GEO'
1361 include 'COMMON.VAR'
1362 include 'COMMON.LOCAL'
1363 include 'COMMON.CHAIN'
1364 include 'COMMON.DERIV'
1365 include 'COMMON.INTERACT'
1366 include 'COMMON.CONTACTS'
1367 include 'COMMON.TORSION'
1368 include 'COMMON.VECTORS'
1369 include 'COMMON.FFIELD'
1370 include 'COMMON.TIME1'
1371 include 'COMMON.SHIELD'
1372 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1373 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1374 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1375 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1376 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1377 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1379 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1381 double precision scal_el /1.0d0/
1383 double precision scal_el /0.5d0/
1386 C 13-go grudnia roku pamietnego...
1387 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1388 & 0.0d0,1.0d0,0.0d0,
1389 & 0.0d0,0.0d0,1.0d0/
1390 cd write(iout,*) 'In EELEC'
1392 cd write(iout,*) 'Type',i
1393 cd write(iout,*) 'B1',B1(:,i)
1394 cd write(iout,*) 'B2',B2(:,i)
1395 cd write(iout,*) 'CC',CC(:,:,i)
1396 cd write(iout,*) 'DD',DD(:,:,i)
1397 cd write(iout,*) 'EE',EE(:,:,i)
1399 cd call check_vecgrad
1401 C print *,"WCHODZE3"
1402 if (icheckgrad.eq.1) then
1404 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1406 dc_norm(k,i)=dc(k,i)*fac
1408 c write (iout,*) 'i',i,' fac',fac
1411 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1412 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1413 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1414 c call vec_and_deriv
1420 time_mat=time_mat+MPI_Wtime()-time01
1424 cd write (iout,*) 'i=',i
1426 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1429 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1430 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1443 cd print '(a)','Enter EELEC'
1444 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1446 gel_loc_loc(i)=0.0d0
1451 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1453 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1455 do i=iturn3_start,iturn3_end
1456 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1457 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1458 & .or. itype(i-1).eq.ntyp1
1459 & .or. itype(i+4).eq.ntyp1
1464 dx_normi=dc_norm(1,i)
1465 dy_normi=dc_norm(2,i)
1466 dz_normi=dc_norm(3,i)
1467 xmedi=c(1,i)+0.5d0*dxi
1468 ymedi=c(2,i)+0.5d0*dyi
1469 zmedi=c(3,i)+0.5d0*dzi
1470 xmedi=mod(xmedi,boxxsize)
1471 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1472 ymedi=mod(ymedi,boxysize)
1473 if (ymedi.lt.0) ymedi=ymedi+boxysize
1474 zmedi=mod(zmedi,boxzsize)
1475 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1477 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1478 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1479 num_cont_hb(i)=num_conti
1481 do i=iturn4_start,iturn4_end
1482 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1483 & .or. itype(i+3).eq.ntyp1
1484 & .or. itype(i+4).eq.ntyp1
1485 & .or. itype(i+5).eq.ntyp1
1486 & .or. itype(i-1).eq.ntyp1
1491 dx_normi=dc_norm(1,i)
1492 dy_normi=dc_norm(2,i)
1493 dz_normi=dc_norm(3,i)
1494 xmedi=c(1,i)+0.5d0*dxi
1495 ymedi=c(2,i)+0.5d0*dyi
1496 zmedi=c(3,i)+0.5d0*dzi
1497 xmedi=mod(xmedi,boxxsize)
1498 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1499 ymedi=mod(ymedi,boxysize)
1500 if (ymedi.lt.0) ymedi=ymedi+boxysize
1501 zmedi=mod(zmedi,boxzsize)
1502 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1503 num_conti=num_cont_hb(i)
1504 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1505 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1506 & call eturn4(i,eello_turn4)
1507 num_cont_hb(i)=num_conti
1510 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1512 do i=iatel_s,iatel_e
1513 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1514 & .or. itype(i+2).eq.ntyp1
1515 & .or. itype(i-1).eq.ntyp1
1520 dx_normi=dc_norm(1,i)
1521 dy_normi=dc_norm(2,i)
1522 dz_normi=dc_norm(3,i)
1523 xmedi=c(1,i)+0.5d0*dxi
1524 ymedi=c(2,i)+0.5d0*dyi
1525 zmedi=c(3,i)+0.5d0*dzi
1526 xmedi=mod(xmedi,boxxsize)
1527 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1528 ymedi=mod(ymedi,boxysize)
1529 if (ymedi.lt.0) ymedi=ymedi+boxysize
1530 zmedi=mod(zmedi,boxzsize)
1531 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1532 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1533 num_conti=num_cont_hb(i)
1534 do j=ielstart(i),ielend(i)
1535 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1536 & .or.itype(j+2).eq.ntyp1
1537 & .or.itype(j-1).eq.ntyp1
1539 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1541 num_cont_hb(i)=num_conti
1543 c write (iout,*) "Number of loop steps in EELEC:",ind
1545 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1546 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1548 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1549 ccc eel_loc=eel_loc+eello_turn3
1550 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1553 C-------------------------------------------------------------------------------
1554 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1555 implicit real*8 (a-h,o-z)
1556 include 'DIMENSIONS'
1560 include 'COMMON.CONTROL'
1561 include 'COMMON.IOUNITS'
1562 include 'COMMON.GEO'
1563 include 'COMMON.VAR'
1564 include 'COMMON.LOCAL'
1565 include 'COMMON.CHAIN'
1566 include 'COMMON.DERIV'
1567 include 'COMMON.INTERACT'
1568 include 'COMMON.CONTACTS'
1569 include 'COMMON.TORSION'
1570 include 'COMMON.VECTORS'
1571 include 'COMMON.FFIELD'
1572 include 'COMMON.TIME1'
1573 include 'COMMON.SHIELD'
1574 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1575 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1576 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1577 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1578 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1579 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1581 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1583 double precision scal_el /1.0d0/
1585 double precision scal_el /0.5d0/
1588 C 13-go grudnia roku pamietnego...
1589 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1590 & 0.0d0,1.0d0,0.0d0,
1591 & 0.0d0,0.0d0,1.0d0/
1592 c time00=MPI_Wtime()
1593 cd write (iout,*) "eelecij",i,j
1594 C print *,"WCHODZE2"
1598 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1599 aaa=app(iteli,itelj)
1600 bbb=bpp(iteli,itelj)
1601 ael6i=ael6(iteli,itelj)
1602 ael3i=ael3(iteli,itelj)
1606 dx_normj=dc_norm(1,j)
1607 dy_normj=dc_norm(2,j)
1608 dz_normj=dc_norm(3,j)
1613 if (xj.lt.0) xj=xj+boxxsize
1615 if (yj.lt.0) yj=yj+boxysize
1617 if (zj.lt.0) zj=zj+boxzsize
1618 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1626 xj=xj_safe+xshift*boxxsize
1627 yj=yj_safe+yshift*boxysize
1628 zj=zj_safe+zshift*boxzsize
1629 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1630 if(dist_temp.lt.dist_init) then
1640 if (isubchap.eq.1) then
1650 rij=xj*xj+yj*yj+zj*zj
1654 c For extracting the short-range part of Evdwpp
1655 sss=sscale(rij/rpp(iteli,itelj))
1656 sssgrad=sscagrad(rij/rpp(iteli,itelj))
1659 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1660 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1661 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1662 fac=cosa-3.0D0*cosb*cosg
1664 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1665 if (j.eq.i+2) ev1=scal_el*ev1
1670 if (shield_mode.eq.0) then
1674 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1676 el1=el1*fac_shield(i)**2*fac_shield(j)**2
1677 el2=el2*fac_shield(i)**2*fac_shield(j)**2
1679 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1680 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1682 evdw1=evdw1+evdwij*(1.0d0-sss)
1683 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1684 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1685 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1686 cd & xmedi,ymedi,zmedi,xj,yj,zj
1688 if (energy_dec) then
1689 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1690 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1694 C Calculate contributions to the Cartesian gradient.
1697 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1698 facel=-3*rrmij*(el1+eesij)
1704 * Radial derivatives. First process both termini of the fragment (i,j)
1709 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1710 & (shield_mode.gt.0)) then
1712 do ilist=1,ishield_list(i)
1713 iresshield=shield_list(ilist,i)
1715 rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
1717 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1719 & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
1720 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1721 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1722 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1723 C if (iresshield.gt.i) then
1724 C do ishi=i+1,iresshield-1
1725 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1726 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1730 C do ishi=iresshield,i
1731 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1732 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1738 do ilist=1,ishield_list(j)
1739 iresshield=shield_list(ilist,j)
1741 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1743 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1745 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
1746 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1751 gshieldc(k,i)=gshieldc(k,i)+
1752 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
1753 gshieldc(k,j)=gshieldc(k,j)+
1754 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
1755 gshieldc(k,i-1)=gshieldc(k,i-1)+
1756 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
1757 gshieldc(k,j-1)=gshieldc(k,j-1)+
1758 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
1764 c ghalf=0.5D0*ggg(k)
1765 c gelc(k,i)=gelc(k,i)+ghalf
1766 c gelc(k,j)=gelc(k,j)+ghalf
1768 c 9/28/08 AL Gradient compotents will be summed only at the end
1770 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1771 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1774 * Loop over residues i+1 thru j-1.
1778 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1781 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1782 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1783 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1785 c ghalf=0.5D0*ggg(k)
1786 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1787 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1789 c 9/28/08 AL Gradient compotents will be summed only at the end
1791 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1792 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1795 * Loop over residues i+1 thru j-1.
1799 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1803 facvdw=ev1+evdwij*(1.0d0-sss)
1806 fac=-3*rrmij*(facvdw+facvdw+facel)
1811 * Radial derivatives. First process both termini of the fragment (i,j)
1817 c ghalf=0.5D0*ggg(k)
1818 c gelc(k,i)=gelc(k,i)+ghalf
1819 c gelc(k,j)=gelc(k,j)+ghalf
1821 c 9/28/08 AL Gradient compotents will be summed only at the end
1823 gelc_long(k,j)=gelc(k,j)+ggg(k)
1824 gelc_long(k,i)=gelc(k,i)-ggg(k)
1827 * Loop over residues i+1 thru j-1.
1831 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1834 c 9/28/08 AL Gradient compotents will be summed only at the end
1838 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1839 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1840 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1842 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1843 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1849 ecosa=2.0D0*fac3*fac1+fac4
1852 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1853 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1855 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1856 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1858 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1859 cd & (dcosg(k),k=1,3)
1861 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
1862 & fac_shield(i)**2*fac_shield(j)**2
1866 c ghalf=0.5D0*ggg(k)
1867 c gelc(k,i)=gelc(k,i)+ghalf
1868 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1869 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1870 c gelc(k,j)=gelc(k,j)+ghalf
1871 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1872 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1876 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1881 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1882 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
1883 & *fac_shield(i)**2*fac_shield(j)**2
1886 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1887 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
1888 & *fac_shield(i)**2*fac_shield(j)**2
1889 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1890 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1892 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1893 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1894 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1896 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1897 C energy of a peptide unit is assumed in the form of a second-order
1898 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1899 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1900 C are computed for EVERY pair of non-contiguous peptide groups.
1902 if (j.lt.nres-1) then
1913 muij(kkk)=mu(k,i)*mu(l,j)
1916 cd write (iout,*) 'EELEC: i',i,' j',j
1917 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1918 cd write(iout,*) 'muij',muij
1919 ury=scalar(uy(1,i),erij)
1920 urz=scalar(uz(1,i),erij)
1921 vry=scalar(uy(1,j),erij)
1922 vrz=scalar(uz(1,j),erij)
1923 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1924 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1925 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1926 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1927 fac=dsqrt(-ael6i)*r3ij
1932 cd write (iout,'(4i5,4f10.5)')
1933 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1934 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1935 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1936 cd & uy(:,j),uz(:,j)
1937 cd write (iout,'(4f10.5)')
1938 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1939 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1940 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
1941 cd write (iout,'(9f10.5/)')
1942 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1943 C Derivatives of the elements of A in virtual-bond vectors
1944 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1946 uryg(k,1)=scalar(erder(1,k),uy(1,i))
1947 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1948 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1949 urzg(k,1)=scalar(erder(1,k),uz(1,i))
1950 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1951 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1952 vryg(k,1)=scalar(erder(1,k),uy(1,j))
1953 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1954 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1955 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1956 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1957 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1959 C Compute radial contributions to the gradient
1977 C Add the contributions coming from er
1980 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1981 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1982 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1983 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1986 C Derivatives in DC(i)
1987 cgrad ghalf1=0.5d0*agg(k,1)
1988 cgrad ghalf2=0.5d0*agg(k,2)
1989 cgrad ghalf3=0.5d0*agg(k,3)
1990 cgrad ghalf4=0.5d0*agg(k,4)
1991 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
1992 & -3.0d0*uryg(k,2)*vry)!+ghalf1
1993 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
1994 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
1995 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
1996 & -3.0d0*urzg(k,2)*vry)!+ghalf3
1997 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
1998 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
1999 C Derivatives in DC(i+1)
2000 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2001 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2002 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2003 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2004 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2005 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2006 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2007 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2008 C Derivatives in DC(j)
2009 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2010 & -3.0d0*vryg(k,2)*ury)!+ghalf1
2011 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2012 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
2013 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2014 & -3.0d0*vryg(k,2)*urz)!+ghalf3
2015 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2016 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
2017 C Derivatives in DC(j+1) or DC(nres-1)
2018 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2019 & -3.0d0*vryg(k,3)*ury)
2020 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2021 & -3.0d0*vrzg(k,3)*ury)
2022 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2023 & -3.0d0*vryg(k,3)*urz)
2024 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2025 & -3.0d0*vrzg(k,3)*urz)
2026 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
2028 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
2041 aggi(k,l)=-aggi(k,l)
2042 aggi1(k,l)=-aggi1(k,l)
2043 aggj(k,l)=-aggj(k,l)
2044 aggj1(k,l)=-aggj1(k,l)
2047 if (j.lt.nres-1) then
2053 aggi(k,l)=-aggi(k,l)
2054 aggi1(k,l)=-aggi1(k,l)
2055 aggj(k,l)=-aggj(k,l)
2056 aggj1(k,l)=-aggj1(k,l)
2067 aggi(k,l)=-aggi(k,l)
2068 aggi1(k,l)=-aggi1(k,l)
2069 aggj(k,l)=-aggj(k,l)
2070 aggj1(k,l)=-aggj1(k,l)
2075 IF (wel_loc.gt.0.0d0) THEN
2076 C Contribution to the local-electrostatic energy coming from the i-j pair
2077 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2079 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2081 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2082 & 'eelloc',i,j,eel_loc_ij
2085 if (shield_mode.eq.0) then
2092 eel_loc_ij=eel_loc_ij
2093 & *fac_shield(i)*fac_shield(j)
2094 eel_loc=eel_loc+eel_loc_ij
2096 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2097 & (shield_mode.gt.0)) then
2100 do ilist=1,ishield_list(i)
2101 iresshield=shield_list(ilist,i)
2103 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2106 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2108 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2109 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2113 do ilist=1,ishield_list(j)
2114 iresshield=shield_list(ilist,j)
2116 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2119 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2121 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2122 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2129 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2130 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2131 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2132 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2133 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2134 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2135 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2136 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2140 C Partial derivatives in virtual-bond dihedral angles gamma
2142 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2143 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2144 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2145 & *fac_shield(i)*fac_shield(j)
2147 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2148 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2149 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2150 & *fac_shield(i)*fac_shield(j)
2152 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2154 ggg(l)=(agg(l,1)*muij(1)+
2155 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2156 & *fac_shield(i)*fac_shield(j)
2158 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2159 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2160 cgrad ghalf=0.5d0*ggg(l)
2161 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
2162 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
2166 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2169 C Remaining derivatives of eello
2171 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2172 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2173 & *fac_shield(i)*fac_shield(j)
2175 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2176 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2177 & *fac_shield(i)*fac_shield(j)
2179 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2180 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2181 & *fac_shield(i)*fac_shield(j)
2183 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2184 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2185 & *fac_shield(i)*fac_shield(j)
2189 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2190 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
2191 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2192 & .and. num_conti.le.maxconts) then
2193 c write (iout,*) i,j," entered corr"
2195 C Calculate the contact function. The ith column of the array JCONT will
2196 C contain the numbers of atoms that make contacts with the atom I (of numbers
2197 C greater than I). The arrays FACONT and GACONT will contain the values of
2198 C the contact function and its derivative.
2199 c r0ij=1.02D0*rpp(iteli,itelj)
2200 c r0ij=1.11D0*rpp(iteli,itelj)
2201 r0ij=2.20D0*rpp(iteli,itelj)
2202 c r0ij=1.55D0*rpp(iteli,itelj)
2203 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2204 if (fcont.gt.0.0D0) then
2205 num_conti=num_conti+1
2206 if (num_conti.gt.maxconts) then
2207 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2208 & ' will skip next contacts for this conf.'
2210 jcont_hb(num_conti,i)=j
2211 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
2212 cd & " jcont_hb",jcont_hb(num_conti,i)
2213 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2214 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2215 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2217 d_cont(num_conti,i)=rij
2218 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2219 C --- Electrostatic-interaction matrix ---
2220 a_chuj(1,1,num_conti,i)=a22
2221 a_chuj(1,2,num_conti,i)=a23
2222 a_chuj(2,1,num_conti,i)=a32
2223 a_chuj(2,2,num_conti,i)=a33
2224 C --- Gradient of rij
2226 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2233 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2234 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2235 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2236 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2237 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2242 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2243 C Calculate contact energies
2245 wij=cosa-3.0D0*cosb*cosg
2248 c fac3=dsqrt(-ael6i)/r0ij**3
2249 fac3=dsqrt(-ael6i)*r3ij
2250 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2251 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2252 if (ees0tmp.gt.0) then
2253 ees0pij=dsqrt(ees0tmp)
2257 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2258 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2259 if (ees0tmp.gt.0) then
2260 ees0mij=dsqrt(ees0tmp)
2265 if (shield_mode.eq.0) then
2269 ees0plist(num_conti,i)=j
2270 C fac_shield(i)=0.4d0
2271 C fac_shield(j)=0.6d0
2273 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2274 & *fac_shield(i)*fac_shield(j)
2275 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2276 & *fac_shield(i)*fac_shield(j)
2278 C Diagnostics. Comment out or remove after debugging!
2279 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2280 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2281 c ees0m(num_conti,i)=0.0D0
2283 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2284 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2285 C Angular derivatives of the contact function
2286 ees0pij1=fac3/ees0pij
2287 ees0mij1=fac3/ees0mij
2288 fac3p=-3.0D0*fac3*rrmij
2289 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2290 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2292 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2293 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2294 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2295 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2296 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2297 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2298 ecosap=ecosa1+ecosa2
2299 ecosbp=ecosb1+ecosb2
2300 ecosgp=ecosg1+ecosg2
2301 ecosam=ecosa1-ecosa2
2302 ecosbm=ecosb1-ecosb2
2303 ecosgm=ecosg1-ecosg2
2312 facont_hb(num_conti,i)=fcont
2313 fprimcont=fprimcont/rij
2314 cd facont_hb(num_conti,i)=1.0D0
2315 C Following line is for diagnostics.
2318 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2319 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2322 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2323 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2325 gggp(1)=gggp(1)+ees0pijp*xj
2326 gggp(2)=gggp(2)+ees0pijp*yj
2327 gggp(3)=gggp(3)+ees0pijp*zj
2328 gggm(1)=gggm(1)+ees0mijp*xj
2329 gggm(2)=gggm(2)+ees0mijp*yj
2330 gggm(3)=gggm(3)+ees0mijp*zj
2331 C Derivatives due to the contact function
2332 gacont_hbr(1,num_conti,i)=fprimcont*xj
2333 gacont_hbr(2,num_conti,i)=fprimcont*yj
2334 gacont_hbr(3,num_conti,i)=fprimcont*zj
2337 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2338 c following the change of gradient-summation algorithm.
2340 cgrad ghalfp=0.5D0*gggp(k)
2341 cgrad ghalfm=0.5D0*gggm(k)
2342 gacontp_hb1(k,num_conti,i)=!ghalfp
2343 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2344 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2345 & *fac_shield(i)*fac_shield(j)
2347 gacontp_hb2(k,num_conti,i)=!ghalfp
2348 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2349 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2350 & *fac_shield(i)*fac_shield(j)
2352 gacontp_hb3(k,num_conti,i)=gggp(k)
2353 & *fac_shield(i)*fac_shield(j)
2355 gacontm_hb1(k,num_conti,i)=!ghalfm
2356 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2357 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2358 & *fac_shield(i)*fac_shield(j)
2360 gacontm_hb2(k,num_conti,i)=!ghalfm
2361 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2362 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2363 & *fac_shield(i)*fac_shield(j)
2365 gacontm_hb3(k,num_conti,i)=gggm(k)
2366 & *fac_shield(i)*fac_shield(j)
2370 endif ! num_conti.le.maxconts
2373 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2376 ghalf=0.5d0*agg(l,k)
2377 aggi(l,k)=aggi(l,k)+ghalf
2378 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2379 aggj(l,k)=aggj(l,k)+ghalf
2382 if (j.eq.nres-1 .and. i.lt.j-2) then
2385 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2390 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2393 C-----------------------------------------------------------------------
2394 subroutine evdwpp_short(evdw1)
2398 implicit real*8 (a-h,o-z)
2399 include 'DIMENSIONS'
2400 include 'COMMON.CONTROL'
2401 include 'COMMON.IOUNITS'
2402 include 'COMMON.GEO'
2403 include 'COMMON.VAR'
2404 include 'COMMON.LOCAL'
2405 include 'COMMON.CHAIN'
2406 include 'COMMON.DERIV'
2407 include 'COMMON.INTERACT'
2408 include 'COMMON.CONTACTS'
2409 include 'COMMON.TORSION'
2410 include 'COMMON.VECTORS'
2411 include 'COMMON.FFIELD'
2413 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2415 double precision scal_el /1.0d0/
2417 double precision scal_el /0.5d0/
2421 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2422 c & " iatel_e_vdw",iatel_e_vdw
2424 do i=iatel_s_vdw,iatel_e_vdw
2425 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2429 dx_normi=dc_norm(1,i)
2430 dy_normi=dc_norm(2,i)
2431 dz_normi=dc_norm(3,i)
2432 xmedi=c(1,i)+0.5d0*dxi
2433 ymedi=c(2,i)+0.5d0*dyi
2434 zmedi=c(3,i)+0.5d0*dzi
2435 xmedi=mod(xmedi,boxxsize)
2436 if (xmedi.lt.0) xmedi=xmedi+boxxsize
2437 ymedi=mod(ymedi,boxysize)
2438 if (ymedi.lt.0) ymedi=ymedi+boxysize
2439 zmedi=mod(zmedi,boxzsize)
2440 if (zmedi.lt.0) zmedi=zmedi+boxzsize
2442 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2443 c & ' ielend',ielend_vdw(i)
2445 do j=ielstart_vdw(i),ielend_vdw(i)
2446 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2450 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2451 aaa=app(iteli,itelj)
2452 bbb=bpp(iteli,itelj)
2456 dx_normj=dc_norm(1,j)
2457 dy_normj=dc_norm(2,j)
2458 dz_normj=dc_norm(3,j)
2463 if (xj.lt.0) xj=xj+boxxsize
2465 if (yj.lt.0) yj=yj+boxysize
2467 if (zj.lt.0) zj=zj+boxzsize
2468 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2476 xj=xj_safe+xshift*boxxsize
2477 yj=yj_safe+yshift*boxysize
2478 zj=zj_safe+zshift*boxzsize
2479 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2480 if(dist_temp.lt.dist_init) then
2490 if (isubchap.eq.1) then
2499 rij=xj*xj+yj*yj+zj*zj
2502 sss=sscale(rij/rpp(iteli,itelj))
2503 sssgrad=sscagrad(rij/rpp(iteli,itelj))
2504 if (sss.gt.0.0d0) then
2509 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2510 if (j.eq.i+2) ev1=scal_el*ev1
2513 if (energy_dec) then
2514 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2516 evdw1=evdw1+evdwij*sss
2518 C Calculate contributions to the Cartesian gradient.
2520 facvdw=-6*rrmij*(ev1+evdwij)*sss
2521 ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2522 ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2523 ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2528 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2529 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2536 C-----------------------------------------------------------------------------
2537 subroutine escp_long(evdw2,evdw2_14)
2539 C This subroutine calculates the excluded-volume interaction energy between
2540 C peptide-group centers and side chains and its gradient in virtual-bond and
2541 C side-chain vectors.
2543 implicit real*8 (a-h,o-z)
2544 include 'DIMENSIONS'
2545 include 'COMMON.GEO'
2546 include 'COMMON.VAR'
2547 include 'COMMON.LOCAL'
2548 include 'COMMON.CHAIN'
2549 include 'COMMON.DERIV'
2550 include 'COMMON.INTERACT'
2551 include 'COMMON.FFIELD'
2552 include 'COMMON.IOUNITS'
2553 include 'COMMON.CONTROL'
2557 CD print '(a)','Enter ESCP KURWA'
2558 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2559 do i=iatscp_s,iatscp_e
2560 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2562 xi=0.5D0*(c(1,i)+c(1,i+1))
2563 yi=0.5D0*(c(2,i)+c(2,i+1))
2564 zi=0.5D0*(c(3,i)+c(3,i+1))
2566 if (xi.lt.0) xi=xi+boxxsize
2568 if (yi.lt.0) yi=yi+boxysize
2570 if (zi.lt.0) zi=zi+boxzsize
2571 do iint=1,nscp_gr(i)
2573 do j=iscpstart(i,iint),iscpend(i,iint)
2575 if (itypj.eq.ntyp1) cycle
2576 C Uncomment following three lines for SC-p interactions
2580 C Uncomment following three lines for Ca-p interactions
2584 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2592 xj=xj_safe+xshift*boxxsize
2593 yj=yj_safe+yshift*boxysize
2594 zj=zj_safe+zshift*boxzsize
2595 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2596 if(dist_temp.lt.dist_init) then
2606 if (subchap.eq.1) then
2616 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2618 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2619 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2620 if (sss.lt.1.0d0) then
2622 e1=fac*fac*aad(itypj,iteli)
2623 e2=fac*bad(itypj,iteli)
2624 if (iabs(j-i) .le. 2) then
2627 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2630 evdw2=evdw2+evdwij*(1.0d0-sss)
2631 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2632 & 'evdw2',i,j,sss,evdwij
2634 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2637 fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2638 fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2642 C Uncomment following three lines for SC-p interactions
2644 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2646 C Uncomment following line for SC-p interactions
2647 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2649 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2650 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2659 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2660 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2661 gradx_scp(j,i)=expon*gradx_scp(j,i)
2664 C******************************************************************************
2668 C To save time the factor EXPON has been extracted from ALL components
2669 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2672 C******************************************************************************
2675 C-----------------------------------------------------------------------------
2676 subroutine escp_short(evdw2,evdw2_14)
2678 C This subroutine calculates the excluded-volume interaction energy between
2679 C peptide-group centers and side chains and its gradient in virtual-bond and
2680 C side-chain vectors.
2682 implicit real*8 (a-h,o-z)
2683 include 'DIMENSIONS'
2684 include 'COMMON.GEO'
2685 include 'COMMON.VAR'
2686 include 'COMMON.LOCAL'
2687 include 'COMMON.CHAIN'
2688 include 'COMMON.DERIV'
2689 include 'COMMON.INTERACT'
2690 include 'COMMON.FFIELD'
2691 include 'COMMON.IOUNITS'
2692 include 'COMMON.CONTROL'
2696 cd print '(a)','Enter ESCP'
2697 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2698 do i=iatscp_s,iatscp_e
2699 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2701 xi=0.5D0*(c(1,i)+c(1,i+1))
2702 yi=0.5D0*(c(2,i)+c(2,i+1))
2703 zi=0.5D0*(c(3,i)+c(3,i+1))
2705 if (xi.lt.0) xi=xi+boxxsize
2707 if (yi.lt.0) yi=yi+boxysize
2709 if (zi.lt.0) zi=zi+boxzsize
2711 do iint=1,nscp_gr(i)
2713 do j=iscpstart(i,iint),iscpend(i,iint)
2715 if (itypj.eq.ntyp1) cycle
2716 C Uncomment following three lines for SC-p interactions
2720 C Uncomment following three lines for Ca-p interactions
2724 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2732 xj=xj_safe+xshift*boxxsize
2733 yj=yj_safe+yshift*boxysize
2734 zj=zj_safe+zshift*boxzsize
2735 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2736 if(dist_temp.lt.dist_init) then
2746 if (subchap.eq.1) then
2755 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2756 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2757 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2758 if (sss.gt.0.0d0) then
2761 e1=fac*fac*aad(itypj,iteli)
2762 e2=fac*bad(itypj,iteli)
2763 if (iabs(j-i) .le. 2) then
2766 evdw2_14=evdw2_14+(e1+e2)*sss
2769 evdw2=evdw2+evdwij*sss
2770 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2771 & 'evdw2',i,j,sss,evdwij
2773 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2775 fac=-(evdwij+e1)*rrij*sss
2776 fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2780 C Uncomment following three lines for SC-p interactions
2782 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2784 C Uncomment following line for SC-p interactions
2785 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2787 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2788 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2797 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2798 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2799 gradx_scp(j,i)=expon*gradx_scp(j,i)
2802 C******************************************************************************
2806 C To save time the factor EXPON has been extracted from ALL components
2807 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2810 C******************************************************************************