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'
645 integer xshift,yshift,zshift
647 ccccc energy_dec=.false.
648 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
651 c if (icall.eq.0) lprn=.false.
655 if (itypi.eq.ntyp1) cycle
661 if (xi.lt.0) xi=xi+boxxsize
663 if (yi.lt.0) yi=yi+boxysize
665 if (zi.lt.0) zi=zi+boxzsize
666 dxi=dc_norm(1,nres+i)
667 dyi=dc_norm(2,nres+i)
668 dzi=dc_norm(3,nres+i)
669 c dsci_inv=dsc_inv(itypi)
670 dsci_inv=vbld_inv(i+nres)
671 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
672 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
674 C Calculate SC interaction energy.
677 do j=istart(i,iint),iend(i,iint)
680 if (itypj.eq.ntyp1) cycle
681 c dscj_inv=dsc_inv(itypj)
682 dscj_inv=vbld_inv(j+nres)
683 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
684 c & 1.0d0/vbld(j+nres)
685 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
686 sig0ij=sigma(itypi,itypj)
687 chi1=chi(itypi,itypj)
688 chi2=chi(itypj,itypi)
695 alf12=0.5D0*(alf1+alf2)
700 if (xj.lt.0) xj=xj+boxxsize
702 if (yj.lt.0) yj=yj+boxysize
704 if (zj.lt.0) zj=zj+boxzsize
705 if ((zj.gt.bordlipbot)
706 &.and.(zj.lt.bordliptop)) then
707 C the energy transfer exist
708 if (zj.lt.buflipbot) then
709 C what fraction I am in
711 & ((positi-bordlipbot)/lipbufthick)
712 C lipbufthick is thickenes of lipid buffore
713 sslipj=sscalelip(fracinbuf)
714 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
715 elseif (zi.gt.bufliptop) then
716 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
717 sslipj=sscalelip(fracinbuf)
718 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
727 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
728 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
729 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
730 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
732 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
740 xj=xj_safe+xshift*boxxsize
741 yj=yj_safe+yshift*boxysize
742 zj=zj_safe+zshift*boxzsize
743 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
744 if(dist_temp.lt.dist_init) then
754 if (subchap.eq.1) then
763 dxj=dc_norm(1,nres+j)
764 dyj=dc_norm(2,nres+j)
765 dzj=dc_norm(3,nres+j)
766 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
768 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
769 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
770 if (sss.lt.1.0d0) then
772 C Calculate angle-dependent terms of energy and contributions to their
776 sig=sig0ij*dsqrt(sigsq)
777 rij_shift=1.0D0/rij-sig+sig0ij
778 c for diagnostics; uncomment
779 c rij_shift=1.2*sig0ij
780 C I hate to put IF's in the loops, but here don't have another choice!!!!
781 if (rij_shift.le.0.0D0) then
783 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
784 cd & restyp(itypi),i,restyp(itypj),j,
785 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
789 c---------------------------------------------------------------
790 rij_shift=1.0D0/rij_shift
794 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
795 eps2der=evdwij*eps3rt
796 eps3der=evdwij*eps2rt
797 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
798 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
799 evdwij=evdwij*eps2rt*eps3rt
800 evdw=evdw+evdwij*(1.0d0-sss)
802 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
804 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
805 & restyp(itypi),i,restyp(itypj),j,
806 & epsi,sigm,chi1,chi2,chip1,chip2,
807 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
808 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
812 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
815 C Calculate gradient components.
816 e1=e1*eps1*eps2rt**2*eps3rt**2
817 fac=-expon*(e1+evdwij)*rij_shift
820 fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij
822 C Calculate the radial part of the gradient
826 gg_lipi(3)=ssgradlipi*evdwij
827 gg_lipj(3)=ssgradlipj*evdwij
828 C Calculate angular part of the gradient.
829 call sc_grad_scale(1.0d0-sss)
834 c write (iout,*) "Number of loop steps in EGB:",ind
835 cccc energy_dec=.false.
838 C-----------------------------------------------------------------------------
839 subroutine egb_short(evdw)
841 C This subroutine calculates the interaction energy of nonbonded side chains
842 C assuming the Gay-Berne potential of interaction.
844 implicit real*8 (a-h,o-z)
848 include 'COMMON.LOCAL'
849 include 'COMMON.CHAIN'
850 include 'COMMON.DERIV'
851 include 'COMMON.NAMES'
852 include 'COMMON.INTERACT'
853 include 'COMMON.IOUNITS'
854 include 'COMMON.CALC'
855 include 'COMMON.CONTROL'
857 integer xshift,yshift,zshift
859 ccccc energy_dec=.false.
860 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
863 c if (icall.eq.0) lprn=.false.
867 if (itypi.eq.ntyp1) cycle
873 if (xi.lt.0) xi=xi+boxxsize
875 if (yi.lt.0) yi=yi+boxysize
877 if (zi.lt.0) zi=zi+boxzsize
878 dxi=dc_norm(1,nres+i)
879 dyi=dc_norm(2,nres+i)
880 dzi=dc_norm(3,nres+i)
881 c dsci_inv=dsc_inv(itypi)
882 dsci_inv=vbld_inv(i+nres)
883 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
884 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
886 C Calculate SC interaction energy.
889 do j=istart(i,iint),iend(i,iint)
892 if (itypj.eq.ntyp1) cycle
893 c dscj_inv=dsc_inv(itypj)
894 dscj_inv=vbld_inv(j+nres)
895 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
896 c & 1.0d0/vbld(j+nres)
897 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
898 sig0ij=sigma(itypi,itypj)
899 chi1=chi(itypi,itypj)
900 chi2=chi(itypj,itypi)
907 alf12=0.5D0*(alf1+alf2)
912 if (xj.lt.0) xj=xj+boxxsize
914 if (yj.lt.0) yj=yj+boxysize
916 if (zj.lt.0) zj=zj+boxzsize
917 if ((zj.gt.bordlipbot)
918 &.and.(zj.lt.bordliptop)) then
919 C the energy transfer exist
920 if (zj.lt.buflipbot) then
921 C what fraction I am in
923 & ((positi-bordlipbot)/lipbufthick)
924 C lipbufthick is thickenes of lipid buffore
925 sslipj=sscalelip(fracinbuf)
926 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
927 elseif (zi.gt.bufliptop) then
928 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
929 sslipj=sscalelip(fracinbuf)
930 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
939 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
940 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
941 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
942 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
943 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
951 xj=xj_safe+xshift*boxxsize
952 yj=yj_safe+yshift*boxysize
953 zj=zj_safe+zshift*boxzsize
954 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
955 if(dist_temp.lt.dist_init) then
965 if (subchap.eq.1) then
974 dxj=dc_norm(1,nres+j)
975 dyj=dc_norm(2,nres+j)
976 dzj=dc_norm(3,nres+j)
977 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
979 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
980 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
981 if (sss.gt.0.0d0) then
983 C Calculate angle-dependent terms of energy and contributions to their
987 sig=sig0ij*dsqrt(sigsq)
988 rij_shift=1.0D0/rij-sig+sig0ij
989 c for diagnostics; uncomment
990 c rij_shift=1.2*sig0ij
991 C I hate to put IF's in the loops, but here don't have another choice!!!!
992 if (rij_shift.le.0.0D0) then
994 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
995 cd & restyp(itypi),i,restyp(itypj),j,
996 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1000 c---------------------------------------------------------------
1001 rij_shift=1.0D0/rij_shift
1002 fac=rij_shift**expon
1005 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1006 eps2der=evdwij*eps3rt
1007 eps3der=evdwij*eps2rt
1008 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1009 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1010 evdwij=evdwij*eps2rt*eps3rt
1011 evdw=evdw+evdwij*sss
1013 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1015 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1016 & restyp(itypi),i,restyp(itypj),j,
1017 & epsi,sigm,chi1,chi2,chip1,chip2,
1018 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1019 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1023 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1026 C Calculate gradient components.
1027 e1=e1*eps1*eps2rt**2*eps3rt**2
1028 fac=-expon*(e1+evdwij)*rij_shift
1031 fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij
1033 C Calculate the radial part of the gradient
1037 gg_lipi(3)=ssgradlipi*evdwij
1038 gg_lipj(3)=ssgradlipj*evdwij
1039 C Calculate angular part of the gradient.
1040 call sc_grad_scale(sss)
1045 c write (iout,*) "Number of loop steps in EGB:",ind
1046 cccc energy_dec=.false.
1049 C-----------------------------------------------------------------------------
1050 subroutine egbv_long(evdw)
1052 C This subroutine calculates the interaction energy of nonbonded side chains
1053 C assuming the Gay-Berne-Vorobjev potential of interaction.
1055 implicit real*8 (a-h,o-z)
1056 include 'DIMENSIONS'
1057 include 'COMMON.GEO'
1058 include 'COMMON.VAR'
1059 include 'COMMON.LOCAL'
1060 include 'COMMON.CHAIN'
1061 include 'COMMON.DERIV'
1062 include 'COMMON.NAMES'
1063 include 'COMMON.INTERACT'
1064 include 'COMMON.IOUNITS'
1065 include 'COMMON.CALC'
1066 common /srutu/ icall
1069 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1072 c if (icall.eq.0) lprn=.true.
1074 do i=iatsc_s,iatsc_e
1076 if (itypi.eq.ntyp1) cycle
1081 dxi=dc_norm(1,nres+i)
1082 dyi=dc_norm(2,nres+i)
1083 dzi=dc_norm(3,nres+i)
1084 c dsci_inv=dsc_inv(itypi)
1085 dsci_inv=vbld_inv(i+nres)
1087 C Calculate SC interaction energy.
1089 do iint=1,nint_gr(i)
1090 do j=istart(i,iint),iend(i,iint)
1093 if (itypj.eq.ntyp1) cycle
1094 c dscj_inv=dsc_inv(itypj)
1095 dscj_inv=vbld_inv(j+nres)
1096 sig0ij=sigma(itypi,itypj)
1097 r0ij=r0(itypi,itypj)
1098 chi1=chi(itypi,itypj)
1099 chi2=chi(itypj,itypi)
1106 alf12=0.5D0*(alf1+alf2)
1110 dxj=dc_norm(1,nres+j)
1111 dyj=dc_norm(2,nres+j)
1112 dzj=dc_norm(3,nres+j)
1113 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1116 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1118 if (sss.lt.1.0d0) then
1120 C Calculate angle-dependent terms of energy and contributions to their
1124 sig=sig0ij*dsqrt(sigsq)
1125 rij_shift=1.0D0/rij-sig+r0ij
1126 C I hate to put IF's in the loops, but here don't have another choice!!!!
1127 if (rij_shift.le.0.0D0) then
1132 c---------------------------------------------------------------
1133 rij_shift=1.0D0/rij_shift
1134 fac=rij_shift**expon
1137 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1138 eps2der=evdwij*eps3rt
1139 eps3der=evdwij*eps2rt
1140 fac_augm=rrij**expon
1141 e_augm=augm(itypi,itypj)*fac_augm
1142 evdwij=evdwij*eps2rt*eps3rt
1143 evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
1145 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1147 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1148 & restyp(itypi),i,restyp(itypj),j,
1149 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1150 & chi1,chi2,chip1,chip2,
1151 & eps1,eps2rt**2,eps3rt**2,
1152 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1155 C Calculate gradient components.
1156 e1=e1*eps1*eps2rt**2*eps3rt**2
1157 fac=-expon*(e1+evdwij)*rij_shift
1159 fac=rij*fac-2*expon*rrij*e_augm
1160 C Calculate the radial part of the gradient
1164 C Calculate angular part of the gradient.
1165 call sc_grad_scale(1.0d0-sss)
1171 C-----------------------------------------------------------------------------
1172 subroutine egbv_short(evdw)
1174 C This subroutine calculates the interaction energy of nonbonded side chains
1175 C assuming the Gay-Berne-Vorobjev potential of interaction.
1177 implicit real*8 (a-h,o-z)
1178 include 'DIMENSIONS'
1179 include 'COMMON.GEO'
1180 include 'COMMON.VAR'
1181 include 'COMMON.LOCAL'
1182 include 'COMMON.CHAIN'
1183 include 'COMMON.DERIV'
1184 include 'COMMON.NAMES'
1185 include 'COMMON.INTERACT'
1186 include 'COMMON.IOUNITS'
1187 include 'COMMON.CALC'
1188 common /srutu/ icall
1191 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1194 c if (icall.eq.0) lprn=.true.
1196 do i=iatsc_s,iatsc_e
1198 if (itypi.eq.ntyp1) cycle
1203 dxi=dc_norm(1,nres+i)
1204 dyi=dc_norm(2,nres+i)
1205 dzi=dc_norm(3,nres+i)
1206 c dsci_inv=dsc_inv(itypi)
1207 dsci_inv=vbld_inv(i+nres)
1209 C Calculate SC interaction energy.
1211 do iint=1,nint_gr(i)
1212 do j=istart(i,iint),iend(i,iint)
1215 if (itypj.eq.ntyp1) cycle
1216 c dscj_inv=dsc_inv(itypj)
1217 dscj_inv=vbld_inv(j+nres)
1218 sig0ij=sigma(itypi,itypj)
1219 r0ij=r0(itypi,itypj)
1220 chi1=chi(itypi,itypj)
1221 chi2=chi(itypj,itypi)
1228 alf12=0.5D0*(alf1+alf2)
1232 dxj=dc_norm(1,nres+j)
1233 dyj=dc_norm(2,nres+j)
1234 dzj=dc_norm(3,nres+j)
1235 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1238 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1240 if (sss.gt.0.0d0) then
1242 C Calculate angle-dependent terms of energy and contributions to their
1246 sig=sig0ij*dsqrt(sigsq)
1247 rij_shift=1.0D0/rij-sig+r0ij
1248 C I hate to put IF's in the loops, but here don't have another choice!!!!
1249 if (rij_shift.le.0.0D0) then
1254 c---------------------------------------------------------------
1255 rij_shift=1.0D0/rij_shift
1256 fac=rij_shift**expon
1259 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1260 eps2der=evdwij*eps3rt
1261 eps3der=evdwij*eps2rt
1262 fac_augm=rrij**expon
1263 e_augm=augm(itypi,itypj)*fac_augm
1264 evdwij=evdwij*eps2rt*eps3rt
1265 evdw=evdw+(evdwij+e_augm)*sss
1267 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1269 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1270 & restyp(itypi),i,restyp(itypj),j,
1271 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1272 & chi1,chi2,chip1,chip2,
1273 & eps1,eps2rt**2,eps3rt**2,
1274 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1277 C Calculate gradient components.
1278 e1=e1*eps1*eps2rt**2*eps3rt**2
1279 fac=-expon*(e1+evdwij)*rij_shift
1281 fac=rij*fac-2*expon*rrij*e_augm
1282 C Calculate the radial part of the gradient
1286 C Calculate angular part of the gradient.
1287 call sc_grad_scale(sss)
1293 C----------------------------------------------------------------------------
1294 subroutine sc_grad_scale(scalfac)
1295 implicit real*8 (a-h,o-z)
1296 include 'DIMENSIONS'
1297 include 'COMMON.CHAIN'
1298 include 'COMMON.DERIV'
1299 include 'COMMON.CALC'
1300 include 'COMMON.IOUNITS'
1301 double precision dcosom1(3),dcosom2(3)
1302 double precision scalfac
1303 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1304 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1305 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1306 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1310 c eom12=evdwij*eps1_om12
1312 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1313 c & " sigder",sigder
1314 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1315 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1317 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1318 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1321 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1323 c write (iout,*) "gg",(gg(k),k=1,3)
1325 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
1326 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1327 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1328 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
1329 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1330 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1331 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1332 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1333 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1334 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1337 C Calculate the components of the gradient in DC and X
1340 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
1341 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
1345 C--------------------------------------------------------------------------
1346 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1348 C This subroutine calculates the average interaction energy and its gradient
1349 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1350 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1351 C The potential depends both on the distance of peptide-group centers and on
1352 C the orientation of the CA-CA virtual bonds.
1354 implicit real*8 (a-h,o-z)
1358 include 'DIMENSIONS'
1359 include 'COMMON.CONTROL'
1360 include 'COMMON.SETUP'
1361 include 'COMMON.IOUNITS'
1362 include 'COMMON.GEO'
1363 include 'COMMON.VAR'
1364 include 'COMMON.LOCAL'
1365 include 'COMMON.CHAIN'
1366 include 'COMMON.DERIV'
1367 include 'COMMON.INTERACT'
1368 include 'COMMON.CONTACTS'
1369 include 'COMMON.TORSION'
1370 include 'COMMON.VECTORS'
1371 include 'COMMON.FFIELD'
1372 include 'COMMON.TIME1'
1373 include 'COMMON.SHIELD'
1374 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1375 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1376 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1377 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1378 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1379 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1381 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1383 double precision scal_el /1.0d0/
1385 double precision scal_el /0.5d0/
1388 C 13-go grudnia roku pamietnego...
1389 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1390 & 0.0d0,1.0d0,0.0d0,
1391 & 0.0d0,0.0d0,1.0d0/
1392 cd write(iout,*) 'In EELEC'
1394 cd write(iout,*) 'Type',i
1395 cd write(iout,*) 'B1',B1(:,i)
1396 cd write(iout,*) 'B2',B2(:,i)
1397 cd write(iout,*) 'CC',CC(:,:,i)
1398 cd write(iout,*) 'DD',DD(:,:,i)
1399 cd write(iout,*) 'EE',EE(:,:,i)
1401 cd call check_vecgrad
1403 C print *,"WCHODZE3"
1404 if (icheckgrad.eq.1) then
1406 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1408 dc_norm(k,i)=dc(k,i)*fac
1410 c write (iout,*) 'i',i,' fac',fac
1413 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1414 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1415 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1416 c call vec_and_deriv
1422 time_mat=time_mat+MPI_Wtime()-time01
1426 cd write (iout,*) 'i=',i
1428 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1431 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1432 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1445 cd print '(a)','Enter EELEC'
1446 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1448 gel_loc_loc(i)=0.0d0
1453 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1455 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1457 do i=iturn3_start,iturn3_end
1458 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1459 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1460 C & .or. itype(i-1).eq.ntyp1
1461 C & .or. itype(i+4).eq.ntyp1
1466 dx_normi=dc_norm(1,i)
1467 dy_normi=dc_norm(2,i)
1468 dz_normi=dc_norm(3,i)
1469 xmedi=c(1,i)+0.5d0*dxi
1470 ymedi=c(2,i)+0.5d0*dyi
1471 zmedi=c(3,i)+0.5d0*dzi
1472 xmedi=mod(xmedi,boxxsize)
1473 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1474 ymedi=mod(ymedi,boxysize)
1475 if (ymedi.lt.0) ymedi=ymedi+boxysize
1476 zmedi=mod(zmedi,boxzsize)
1477 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1479 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1480 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1481 num_cont_hb(i)=num_conti
1483 do i=iturn4_start,iturn4_end
1484 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1485 & .or. itype(i+3).eq.ntyp1
1486 & .or. itype(i+4).eq.ntyp1
1487 C & .or. itype(i+5).eq.ntyp1
1488 C & .or. itype(i-1).eq.ntyp1
1493 dx_normi=dc_norm(1,i)
1494 dy_normi=dc_norm(2,i)
1495 dz_normi=dc_norm(3,i)
1496 xmedi=c(1,i)+0.5d0*dxi
1497 ymedi=c(2,i)+0.5d0*dyi
1498 zmedi=c(3,i)+0.5d0*dzi
1499 xmedi=mod(xmedi,boxxsize)
1500 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1501 ymedi=mod(ymedi,boxysize)
1502 if (ymedi.lt.0) ymedi=ymedi+boxysize
1503 zmedi=mod(zmedi,boxzsize)
1504 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1505 num_conti=num_cont_hb(i)
1506 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1507 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1508 & call eturn4(i,eello_turn4)
1509 num_cont_hb(i)=num_conti
1512 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1514 do i=iatel_s,iatel_e
1515 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1516 C & .or. itype(i+2).eq.ntyp1
1517 C & .or. itype(i-1).eq.ntyp1
1522 dx_normi=dc_norm(1,i)
1523 dy_normi=dc_norm(2,i)
1524 dz_normi=dc_norm(3,i)
1525 xmedi=c(1,i)+0.5d0*dxi
1526 ymedi=c(2,i)+0.5d0*dyi
1527 zmedi=c(3,i)+0.5d0*dzi
1528 xmedi=mod(xmedi,boxxsize)
1529 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1530 ymedi=mod(ymedi,boxysize)
1531 if (ymedi.lt.0) ymedi=ymedi+boxysize
1532 zmedi=mod(zmedi,boxzsize)
1533 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1534 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1535 num_conti=num_cont_hb(i)
1536 do j=ielstart(i),ielend(i)
1537 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1538 C & .or.itype(j+2).eq.ntyp1
1539 C & .or.itype(j-1).eq.ntyp1
1541 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1543 num_cont_hb(i)=num_conti
1545 c write (iout,*) "Number of loop steps in EELEC:",ind
1547 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1548 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1550 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1551 ccc eel_loc=eel_loc+eello_turn3
1552 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1555 C-------------------------------------------------------------------------------
1556 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1557 implicit real*8 (a-h,o-z)
1558 include 'DIMENSIONS'
1562 include 'COMMON.CONTROL'
1563 include 'COMMON.IOUNITS'
1564 include 'COMMON.GEO'
1565 include 'COMMON.VAR'
1566 include 'COMMON.LOCAL'
1567 include 'COMMON.CHAIN'
1568 include 'COMMON.DERIV'
1569 include 'COMMON.INTERACT'
1570 include 'COMMON.CONTACTS'
1571 include 'COMMON.TORSION'
1572 include 'COMMON.VECTORS'
1573 include 'COMMON.FFIELD'
1574 include 'COMMON.TIME1'
1575 include 'COMMON.SHIELD'
1576 integer xshift,yshift,zshift
1577 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1578 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1579 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1580 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
1581 & gmuij2(4),gmuji2(4)
1582 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1583 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1585 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1587 double precision scal_el /1.0d0/
1589 double precision scal_el /0.5d0/
1592 C 13-go grudnia roku pamietnego...
1593 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1594 & 0.0d0,1.0d0,0.0d0,
1595 & 0.0d0,0.0d0,1.0d0/
1596 c time00=MPI_Wtime()
1597 cd write (iout,*) "eelecij",i,j
1598 C print *,"WCHODZE2"
1602 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1603 aaa=app(iteli,itelj)
1604 bbb=bpp(iteli,itelj)
1605 ael6i=ael6(iteli,itelj)
1606 ael3i=ael3(iteli,itelj)
1610 dx_normj=dc_norm(1,j)
1611 dy_normj=dc_norm(2,j)
1612 dz_normj=dc_norm(3,j)
1617 if (xj.lt.0) xj=xj+boxxsize
1619 if (yj.lt.0) yj=yj+boxysize
1621 if (zj.lt.0) zj=zj+boxzsize
1622 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1630 xj=xj_safe+xshift*boxxsize
1631 yj=yj_safe+yshift*boxysize
1632 zj=zj_safe+zshift*boxzsize
1633 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1634 if(dist_temp.lt.dist_init) then
1644 if (isubchap.eq.1) then
1654 rij=xj*xj+yj*yj+zj*zj
1658 c For extracting the short-range part of Evdwpp
1659 sss=sscale(rij/rpp(iteli,itelj))
1660 sssgrad=sscagrad(rij/rpp(iteli,itelj))
1663 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1664 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1665 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1666 fac=cosa-3.0D0*cosb*cosg
1668 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1669 if (j.eq.i+2) ev1=scal_el*ev1
1674 if (shield_mode.eq.0) then
1678 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1680 el1=el1*fac_shield(i)**2*fac_shield(j)**2
1681 el2=el2*fac_shield(i)**2*fac_shield(j)**2
1683 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1684 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1686 evdw1=evdw1+evdwij*(1.0d0-sss)
1687 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1688 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1689 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1690 cd & xmedi,ymedi,zmedi,xj,yj,zj
1692 if (energy_dec) then
1693 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1694 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1698 C Calculate contributions to the Cartesian gradient.
1701 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1702 facel=-3*rrmij*(el1+eesij)
1708 * Radial derivatives. First process both termini of the fragment (i,j)
1713 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1714 & (shield_mode.gt.0)) then
1716 do ilist=1,ishield_list(i)
1717 iresshield=shield_list(ilist,i)
1719 rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
1721 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1723 & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
1724 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1725 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1726 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1727 C if (iresshield.gt.i) then
1728 C do ishi=i+1,iresshield-1
1729 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1730 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1734 C do ishi=iresshield,i
1735 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1736 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1742 do ilist=1,ishield_list(j)
1743 iresshield=shield_list(ilist,j)
1745 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1747 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1749 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
1750 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1755 gshieldc(k,i)=gshieldc(k,i)+
1756 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
1757 gshieldc(k,j)=gshieldc(k,j)+
1758 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
1759 gshieldc(k,i-1)=gshieldc(k,i-1)+
1760 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
1761 gshieldc(k,j-1)=gshieldc(k,j-1)+
1762 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
1768 c ghalf=0.5D0*ggg(k)
1769 c gelc(k,i)=gelc(k,i)+ghalf
1770 c gelc(k,j)=gelc(k,j)+ghalf
1772 c 9/28/08 AL Gradient compotents will be summed only at the end
1774 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1775 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1778 * Loop over residues i+1 thru j-1.
1782 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1785 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1786 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1787 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1789 c ghalf=0.5D0*ggg(k)
1790 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1791 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1793 c 9/28/08 AL Gradient compotents will be summed only at the end
1795 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1796 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1799 * Loop over residues i+1 thru j-1.
1803 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1807 facvdw=ev1+evdwij*(1.0d0-sss)
1810 fac=-3*rrmij*(facvdw+facvdw+facel)
1815 * Radial derivatives. First process both termini of the fragment (i,j)
1821 c ghalf=0.5D0*ggg(k)
1822 c gelc(k,i)=gelc(k,i)+ghalf
1823 c gelc(k,j)=gelc(k,j)+ghalf
1825 c 9/28/08 AL Gradient compotents will be summed only at the end
1827 gelc_long(k,j)=gelc(k,j)+ggg(k)
1828 gelc_long(k,i)=gelc(k,i)-ggg(k)
1831 * Loop over residues i+1 thru j-1.
1835 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1838 c 9/28/08 AL Gradient compotents will be summed only at the end
1842 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1843 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1844 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1846 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1847 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1853 ecosa=2.0D0*fac3*fac1+fac4
1856 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1857 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1859 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1860 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1862 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1863 cd & (dcosg(k),k=1,3)
1865 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
1866 & fac_shield(i)**2*fac_shield(j)**2
1870 c ghalf=0.5D0*ggg(k)
1871 c gelc(k,i)=gelc(k,i)+ghalf
1872 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1873 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1874 c gelc(k,j)=gelc(k,j)+ghalf
1875 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1876 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1880 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1885 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1886 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
1887 & *fac_shield(i)**2*fac_shield(j)**2
1890 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1891 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
1892 & *fac_shield(i)**2*fac_shield(j)**2
1893 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1894 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1896 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1897 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1898 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1900 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1901 C energy of a peptide unit is assumed in the form of a second-order
1902 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1903 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1904 C are computed for EVERY pair of non-contiguous peptide groups.
1906 if (j.lt.nres-1) then
1917 muij(kkk)=mu(k,i)*mu(l,j)
1919 gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
1920 c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
1921 gmuij2(kkk)=gUb2(k,i)*mu(l,j)
1922 gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
1923 c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
1924 gmuji2(kkk)=mu(k,i)*gUb2(l,j)
1928 cd write (iout,*) 'EELEC: i',i,' j',j
1929 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1930 cd write(iout,*) 'muij',muij
1931 ury=scalar(uy(1,i),erij)
1932 urz=scalar(uz(1,i),erij)
1933 vry=scalar(uy(1,j),erij)
1934 vrz=scalar(uz(1,j),erij)
1935 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1936 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1937 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1938 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1939 fac=dsqrt(-ael6i)*r3ij
1944 cd write (iout,'(4i5,4f10.5)')
1945 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1946 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1947 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1948 cd & uy(:,j),uz(:,j)
1949 cd write (iout,'(4f10.5)')
1950 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1951 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1952 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
1953 cd write (iout,'(9f10.5/)')
1954 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1955 C Derivatives of the elements of A in virtual-bond vectors
1956 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1958 uryg(k,1)=scalar(erder(1,k),uy(1,i))
1959 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1960 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1961 urzg(k,1)=scalar(erder(1,k),uz(1,i))
1962 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1963 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1964 vryg(k,1)=scalar(erder(1,k),uy(1,j))
1965 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1966 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1967 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1968 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1969 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1971 C Compute radial contributions to the gradient
1989 C Add the contributions coming from er
1992 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1993 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1994 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1995 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1998 C Derivatives in DC(i)
1999 cgrad ghalf1=0.5d0*agg(k,1)
2000 cgrad ghalf2=0.5d0*agg(k,2)
2001 cgrad ghalf3=0.5d0*agg(k,3)
2002 cgrad ghalf4=0.5d0*agg(k,4)
2003 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2004 & -3.0d0*uryg(k,2)*vry)!+ghalf1
2005 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2006 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
2007 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2008 & -3.0d0*urzg(k,2)*vry)!+ghalf3
2009 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2010 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
2011 C Derivatives in DC(i+1)
2012 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2013 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2014 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2015 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2016 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2017 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2018 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2019 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2020 C Derivatives in DC(j)
2021 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2022 & -3.0d0*vryg(k,2)*ury)!+ghalf1
2023 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2024 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
2025 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2026 & -3.0d0*vryg(k,2)*urz)!+ghalf3
2027 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2028 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
2029 C Derivatives in DC(j+1) or DC(nres-1)
2030 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2031 & -3.0d0*vryg(k,3)*ury)
2032 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2033 & -3.0d0*vrzg(k,3)*ury)
2034 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2035 & -3.0d0*vryg(k,3)*urz)
2036 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2037 & -3.0d0*vrzg(k,3)*urz)
2038 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
2040 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
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)
2059 if (j.lt.nres-1) then
2065 aggi(k,l)=-aggi(k,l)
2066 aggi1(k,l)=-aggi1(k,l)
2067 aggj(k,l)=-aggj(k,l)
2068 aggj1(k,l)=-aggj1(k,l)
2079 aggi(k,l)=-aggi(k,l)
2080 aggi1(k,l)=-aggi1(k,l)
2081 aggj(k,l)=-aggj(k,l)
2082 aggj1(k,l)=-aggj1(k,l)
2087 IF (wel_loc.gt.0.0d0) THEN
2088 C Contribution to the local-electrostatic energy coming from the i-j pair
2089 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2091 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2093 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2094 & 'eelloc',i,j,eel_loc_ij
2097 if (shield_mode.eq.0) then
2104 eel_loc_ij=eel_loc_ij
2105 & *fac_shield(i)*fac_shield(j)
2106 eel_loc=eel_loc+eel_loc_ij
2108 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2109 & (shield_mode.gt.0)) then
2112 do ilist=1,ishield_list(i)
2113 iresshield=shield_list(ilist,i)
2115 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2118 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2120 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2121 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2125 do ilist=1,ishield_list(j)
2126 iresshield=shield_list(ilist,j)
2128 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2131 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2133 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2134 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2141 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2142 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2143 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2144 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2145 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2146 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2147 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2148 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2153 geel_loc_ij=(a22*gmuij1(1)
2157 & *fac_shield(i)*fac_shield(j)
2158 c write(iout,*) "derivative over thatai"
2159 c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
2161 gloc(nphi+i,icg)=gloc(nphi+i,icg)+
2162 & geel_loc_ij*wel_loc
2163 c write(iout,*) "derivative over thatai-1"
2164 c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
2171 gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
2172 & geel_loc_ij*wel_loc
2173 & *fac_shield(i)*fac_shield(j)
2175 c Derivative over j residue
2176 geel_loc_ji=a22*gmuji1(1)
2180 c write(iout,*) "derivative over thataj"
2181 c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
2184 gloc(nphi+j,icg)=gloc(nphi+j,icg)+
2185 & geel_loc_ji*wel_loc
2186 & *fac_shield(i)*fac_shield(j)
2193 c write(iout,*) "derivative over thataj-1"
2194 c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
2196 gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
2197 & geel_loc_ji*wel_loc
2198 & *fac_shield(i)*fac_shield(j)
2200 cC Partial derivatives in virtual-bond dihedral angles gamma
2202 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2203 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2204 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2205 & *fac_shield(i)*fac_shield(j)
2207 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2208 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2209 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2210 & *fac_shield(i)*fac_shield(j)
2212 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2214 ggg(l)=(agg(l,1)*muij(1)+
2215 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2216 & *fac_shield(i)*fac_shield(j)
2218 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2219 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2220 cgrad ghalf=0.5d0*ggg(l)
2221 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
2222 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
2226 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2229 C Remaining derivatives of eello
2231 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2232 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2233 & *fac_shield(i)*fac_shield(j)
2235 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2236 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2237 & *fac_shield(i)*fac_shield(j)
2239 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2240 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2241 & *fac_shield(i)*fac_shield(j)
2243 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2244 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2245 & *fac_shield(i)*fac_shield(j)
2249 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2250 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
2251 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2252 & .and. num_conti.le.maxconts) then
2253 c write (iout,*) i,j," entered corr"
2255 C Calculate the contact function. The ith column of the array JCONT will
2256 C contain the numbers of atoms that make contacts with the atom I (of numbers
2257 C greater than I). The arrays FACONT and GACONT will contain the values of
2258 C the contact function and its derivative.
2259 c r0ij=1.02D0*rpp(iteli,itelj)
2260 c r0ij=1.11D0*rpp(iteli,itelj)
2261 r0ij=2.20D0*rpp(iteli,itelj)
2262 c r0ij=1.55D0*rpp(iteli,itelj)
2263 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2264 if (fcont.gt.0.0D0) then
2265 num_conti=num_conti+1
2266 if (num_conti.gt.maxconts) then
2267 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2268 & ' will skip next contacts for this conf.'
2270 jcont_hb(num_conti,i)=j
2271 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
2272 cd & " jcont_hb",jcont_hb(num_conti,i)
2273 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2274 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2275 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2277 d_cont(num_conti,i)=rij
2278 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2279 C --- Electrostatic-interaction matrix ---
2280 a_chuj(1,1,num_conti,i)=a22
2281 a_chuj(1,2,num_conti,i)=a23
2282 a_chuj(2,1,num_conti,i)=a32
2283 a_chuj(2,2,num_conti,i)=a33
2284 C --- Gradient of rij
2286 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2293 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2294 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2295 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2296 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2297 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2302 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2303 C Calculate contact energies
2305 wij=cosa-3.0D0*cosb*cosg
2308 c fac3=dsqrt(-ael6i)/r0ij**3
2309 fac3=dsqrt(-ael6i)*r3ij
2310 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2311 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2312 if (ees0tmp.gt.0) then
2313 ees0pij=dsqrt(ees0tmp)
2317 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2318 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2319 if (ees0tmp.gt.0) then
2320 ees0mij=dsqrt(ees0tmp)
2325 if (shield_mode.eq.0) then
2329 ees0plist(num_conti,i)=j
2330 C fac_shield(i)=0.4d0
2331 C fac_shield(j)=0.6d0
2333 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2334 & *fac_shield(i)*fac_shield(j)
2335 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2336 & *fac_shield(i)*fac_shield(j)
2338 C Diagnostics. Comment out or remove after debugging!
2339 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2340 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2341 c ees0m(num_conti,i)=0.0D0
2343 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2344 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2345 C Angular derivatives of the contact function
2346 ees0pij1=fac3/ees0pij
2347 ees0mij1=fac3/ees0mij
2348 fac3p=-3.0D0*fac3*rrmij
2349 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2350 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2352 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2353 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2354 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2355 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2356 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2357 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2358 ecosap=ecosa1+ecosa2
2359 ecosbp=ecosb1+ecosb2
2360 ecosgp=ecosg1+ecosg2
2361 ecosam=ecosa1-ecosa2
2362 ecosbm=ecosb1-ecosb2
2363 ecosgm=ecosg1-ecosg2
2372 facont_hb(num_conti,i)=fcont
2373 fprimcont=fprimcont/rij
2374 cd facont_hb(num_conti,i)=1.0D0
2375 C Following line is for diagnostics.
2378 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2379 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2382 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2383 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2385 gggp(1)=gggp(1)+ees0pijp*xj
2386 gggp(2)=gggp(2)+ees0pijp*yj
2387 gggp(3)=gggp(3)+ees0pijp*zj
2388 gggm(1)=gggm(1)+ees0mijp*xj
2389 gggm(2)=gggm(2)+ees0mijp*yj
2390 gggm(3)=gggm(3)+ees0mijp*zj
2391 C Derivatives due to the contact function
2392 gacont_hbr(1,num_conti,i)=fprimcont*xj
2393 gacont_hbr(2,num_conti,i)=fprimcont*yj
2394 gacont_hbr(3,num_conti,i)=fprimcont*zj
2397 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2398 c following the change of gradient-summation algorithm.
2400 cgrad ghalfp=0.5D0*gggp(k)
2401 cgrad ghalfm=0.5D0*gggm(k)
2402 gacontp_hb1(k,num_conti,i)=!ghalfp
2403 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2404 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2405 & *fac_shield(i)*fac_shield(j)
2407 gacontp_hb2(k,num_conti,i)=!ghalfp
2408 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2409 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2410 & *fac_shield(i)*fac_shield(j)
2412 gacontp_hb3(k,num_conti,i)=gggp(k)
2413 & *fac_shield(i)*fac_shield(j)
2415 gacontm_hb1(k,num_conti,i)=!ghalfm
2416 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2417 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2418 & *fac_shield(i)*fac_shield(j)
2420 gacontm_hb2(k,num_conti,i)=!ghalfm
2421 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2422 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2423 & *fac_shield(i)*fac_shield(j)
2425 gacontm_hb3(k,num_conti,i)=gggm(k)
2426 & *fac_shield(i)*fac_shield(j)
2430 endif ! num_conti.le.maxconts
2433 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2436 ghalf=0.5d0*agg(l,k)
2437 aggi(l,k)=aggi(l,k)+ghalf
2438 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2439 aggj(l,k)=aggj(l,k)+ghalf
2442 if (j.eq.nres-1 .and. i.lt.j-2) then
2445 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2450 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2453 C-----------------------------------------------------------------------
2454 subroutine evdwpp_short(evdw1)
2458 implicit real*8 (a-h,o-z)
2459 include 'DIMENSIONS'
2460 include 'COMMON.CONTROL'
2461 include 'COMMON.IOUNITS'
2462 include 'COMMON.GEO'
2463 include 'COMMON.VAR'
2464 include 'COMMON.LOCAL'
2465 include 'COMMON.CHAIN'
2466 include 'COMMON.DERIV'
2467 include 'COMMON.INTERACT'
2468 include 'COMMON.CONTACTS'
2469 include 'COMMON.TORSION'
2470 include 'COMMON.VECTORS'
2471 include 'COMMON.FFIELD'
2473 integer xshift,yshift,zshift
2474 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2476 double precision scal_el /1.0d0/
2478 double precision scal_el /0.5d0/
2480 c write (iout,*) "evdwpp_short"
2483 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2484 c & " iatel_e_vdw",iatel_e_vdw
2486 do i=iatel_s_vdw,iatel_e_vdw
2487 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2491 dx_normi=dc_norm(1,i)
2492 dy_normi=dc_norm(2,i)
2493 dz_normi=dc_norm(3,i)
2494 xmedi=c(1,i)+0.5d0*dxi
2495 ymedi=c(2,i)+0.5d0*dyi
2496 zmedi=c(3,i)+0.5d0*dzi
2497 xmedi=mod(xmedi,boxxsize)
2498 if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
2499 ymedi=mod(ymedi,boxysize)
2500 if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
2501 zmedi=mod(zmedi,boxzsize)
2502 if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
2504 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2505 c & ' ielend',ielend_vdw(i)
2507 do j=ielstart_vdw(i),ielend_vdw(i)
2508 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2512 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2513 aaa=app(iteli,itelj)
2514 bbb=bpp(iteli,itelj)
2518 dx_normj=dc_norm(1,j)
2519 dy_normj=dc_norm(2,j)
2520 dz_normj=dc_norm(3,j)
2525 if (xj.lt.0) xj=xj+boxxsize
2527 if (yj.lt.0) yj=yj+boxysize
2529 if (zj.lt.0) zj=zj+boxzsize
2530 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2538 xj=xj_safe+xshift*boxxsize
2539 yj=yj_safe+yshift*boxysize
2540 zj=zj_safe+zshift*boxzsize
2541 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2542 if(dist_temp.lt.dist_init) then
2552 if (isubchap.eq.1) then
2561 rij=xj*xj+yj*yj+zj*zj
2564 c sss=sscale(rij/rpp(iteli,itelj))
2565 c sssgrad=sscagrad(rij/rpp(iteli,itelj))
2567 sssgrad=sscagrad(rij)
2568 if (sss.gt.0.0d0) then
2573 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2574 if (j.eq.i+2) ev1=scal_el*ev1
2577 if (energy_dec) then
2578 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2580 evdw1=evdw1+evdwij*sss
2581 if (energy_dec) write (iout,'(a10,2i5,0pf7.3)')
2582 & 'evdw1_sum',i,j,evdw1
2584 C Calculate contributions to the Cartesian gradient.
2586 facvdw=-6*rrmij*(ev1+evdwij)*sss
2587 ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2588 ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2589 ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2594 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2595 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2602 C-----------------------------------------------------------------------------
2603 subroutine escp_long(evdw2,evdw2_14)
2605 C This subroutine calculates the excluded-volume interaction energy between
2606 C peptide-group centers and side chains and its gradient in virtual-bond and
2607 C side-chain vectors.
2609 implicit real*8 (a-h,o-z)
2610 include 'DIMENSIONS'
2611 include 'COMMON.GEO'
2612 include 'COMMON.VAR'
2613 include 'COMMON.LOCAL'
2614 include 'COMMON.CHAIN'
2615 include 'COMMON.DERIV'
2616 include 'COMMON.INTERACT'
2617 include 'COMMON.FFIELD'
2618 include 'COMMON.IOUNITS'
2619 include 'COMMON.CONTROL'
2620 logical lprint_short
2621 common /shortcheck/ lprint_short
2623 integer xshift,yshift,zshift
2624 if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
2627 CD print '(a)','Enter ESCP KURWA'
2628 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2630 c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
2631 c & ' iatscp_e=',iatscp_e
2632 do i=iatscp_s,iatscp_e
2633 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2635 xi=0.5D0*(c(1,i)+c(1,i+1))
2636 yi=0.5D0*(c(2,i)+c(2,i+1))
2637 zi=0.5D0*(c(3,i)+c(3,i+1))
2639 if (xi.lt.0) xi=xi+boxxsize
2641 if (yi.lt.0) yi=yi+boxysize
2643 if (zi.lt.0) zi=zi+boxzsize
2644 do iint=1,nscp_gr(i)
2646 do j=iscpstart(i,iint),iscpend(i,iint)
2648 if (itypj.eq.ntyp1) cycle
2649 C Uncomment following three lines for SC-p interactions
2653 C Uncomment following three lines for Ca-p interactions
2659 if (xj.lt.0) xj=xj+boxxsize
2661 if (yj.lt.0) yj=yj+boxysize
2663 if (zj.lt.0) zj=zj+boxzsize
2665 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2673 xj=xj_safe+xshift*boxxsize
2674 yj=yj_safe+yshift*boxysize
2675 zj=zj_safe+zshift*boxzsize
2676 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2677 if(dist_temp.lt.dist_init) then
2687 if (subchap.eq.1) then
2697 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2699 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2700 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2701 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2702 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2703 if (sss.lt.1.0d0) then
2705 e1=fac*fac*aad(itypj,iteli)
2706 e2=fac*bad(itypj,iteli)
2707 if (iabs(j-i) .le. 2) then
2710 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2713 evdw2=evdw2+evdwij*(1.0d0-sss)
2714 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2715 & 'evdw2',i,j,sss,evdwij
2717 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2720 fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2721 fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2725 C Uncomment following three lines for SC-p interactions
2727 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2729 C Uncomment following line for SC-p interactions
2730 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2732 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2733 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2742 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2743 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2744 gradx_scp(j,i)=expon*gradx_scp(j,i)
2747 C******************************************************************************
2751 C To save time the factor EXPON has been extracted from ALL components
2752 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2755 C******************************************************************************
2758 C-----------------------------------------------------------------------------
2759 subroutine escp_short(evdw2,evdw2_14)
2761 C This subroutine calculates the excluded-volume interaction energy between
2762 C peptide-group centers and side chains and its gradient in virtual-bond and
2763 C side-chain vectors.
2765 implicit real*8 (a-h,o-z)
2766 include 'DIMENSIONS'
2767 include 'COMMON.GEO'
2768 include 'COMMON.VAR'
2769 include 'COMMON.LOCAL'
2770 include 'COMMON.CHAIN'
2771 include 'COMMON.DERIV'
2772 include 'COMMON.INTERACT'
2773 include 'COMMON.FFIELD'
2774 include 'COMMON.IOUNITS'
2775 include 'COMMON.CONTROL'
2776 integer xshift,yshift,zshift
2777 logical lprint_short
2778 common /shortcheck/ lprint_short
2782 cd print '(a)','Enter ESCP'
2784 c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
2785 c & ' iatscp_e=',iatscp_e
2786 if (energy_dec) write (iout,*) "escp_short:",r_cut,rlamb
2787 do i=iatscp_s,iatscp_e
2788 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2790 xi=0.5D0*(c(1,i)+c(1,i+1))
2791 yi=0.5D0*(c(2,i)+c(2,i+1))
2792 zi=0.5D0*(c(3,i)+c(3,i+1))
2794 if (xi.lt.0) xi=xi+boxxsize
2796 if (yi.lt.0) yi=yi+boxysize
2798 if (zi.lt.0) zi=zi+boxzsize
2801 c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
2802 c & " nscp_gr",nscp_gr(i)
2803 do iint=1,nscp_gr(i)
2805 do j=iscpstart(i,iint),iscpend(i,iint)
2808 c & write (iout,*) "j",j," itypj",itypj
2809 if (itypj.eq.ntyp1) cycle
2810 C Uncomment following three lines for SC-p interactions
2814 C Uncomment following three lines for Ca-p interactions
2820 if (xj.lt.0) xj=xj+boxxsize
2822 if (yj.lt.0) yj=yj+boxysize
2824 if (zj.lt.0) zj=zj+boxzsize
2826 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2827 c if (lprint_short) then
2828 c write (iout,*) i,j,xi,yi,zi,xj,yj,zj
2829 c write (iout,*) "dist_init",dsqrt(dist_init)
2838 xj=xj_safe+xshift*boxxsize
2839 yj=yj_safe+yshift*boxysize
2840 zj=zj_safe+zshift*boxzsize
2841 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2842 if(dist_temp.lt.dist_init) then
2852 c if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp)
2853 if (subchap.eq.1) then
2862 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2863 c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2864 c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2865 sss=sscale(1.0d0/(dsqrt(rrij)))
2866 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
2867 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2868 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2869 c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),
2870 c & " subchap",subchap," sss",sss
2871 if (sss.gt.0.0d0) then
2874 e1=fac*fac*aad(itypj,iteli)
2875 e2=fac*bad(itypj,iteli)
2876 if (iabs(j-i) .le. 2) then
2879 evdw2_14=evdw2_14+(e1+e2)*sss
2882 evdw2=evdw2+evdwij*sss
2883 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2884 & 'evdw2',i,j,sss,evdwij
2886 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2888 fac=-(evdwij+e1)*rrij*sss
2889 fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2893 C Uncomment following three lines for SC-p interactions
2895 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2897 C Uncomment following line for SC-p interactions
2898 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2900 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2901 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2910 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2911 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2912 gradx_scp(j,i)=expon*gradx_scp(j,i)
2915 C******************************************************************************
2919 C To save time the factor EXPON has been extracted from ALL components
2920 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2923 C******************************************************************************