1 C-----------------------------------------------------------------------
2 double precision function sscalelip(r)
4 double precision r,gamm
5 include "COMMON.SPLITELE"
6 C if(r.lt.r_cut-rlamb) then
8 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
9 C gamm=(r-(r_cut-rlamb))/rlamb
10 sscalelip=1.0d0+r*r*(2*r-3.0d0)
16 C-----------------------------------------------------------------------
17 double precision function sscagradlip(r)
19 double precision r,gamm
20 include "COMMON.SPLITELE"
21 C if(r.lt.r_cut-rlamb) then
23 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
24 C gamm=(r-(r_cut-rlamb))/rlamb
25 sscagradlip=r*(6*r-6.0d0)
32 C-----------------------------------------------------------------------
33 double precision function sscale(r)
35 double precision r,gamm
36 include "COMMON.SPLITELE"
37 if(r.lt.r_cut-rlamb) then
39 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
40 gamm=(r-(r_cut-rlamb))/rlamb
41 sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
47 C-----------------------------------------------------------------------
48 double precision function sscagrad(r)
50 double precision r,gamm
51 include "COMMON.SPLITELE"
52 if(r.lt.r_cut-rlamb) then
54 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
55 gamm=(r-(r_cut-rlamb))/rlamb
56 sscagrad=gamm*(6*gamm-6.0d0)/rlamb
62 C-----------------------------------------------------------------------
63 subroutine elj_long(evdw)
65 C This subroutine calculates the interaction energy of nonbonded side chains
66 C assuming the LJ potential of interaction.
68 implicit real*8 (a-h,o-z)
70 parameter (accur=1.0d-10)
73 include 'COMMON.LOCAL'
74 include 'COMMON.CHAIN'
75 include 'COMMON.DERIV'
76 include 'COMMON.INTERACT'
77 include 'COMMON.TORSION'
78 include 'COMMON.SBRIDGE'
79 include 'COMMON.NAMES'
80 include 'COMMON.IOUNITS'
81 c include 'COMMON.CONTACTS'
83 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
87 if (itypi.eq.ntyp1) cycle
88 itypi1=iabs(itype(i+1))
93 C Calculate SC interaction energy.
96 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
97 cd & 'iend=',iend(i,iint)
98 do j=istart(i,iint),iend(i,iint)
100 if (itypj.eq.ntyp1) cycle
104 rij=xj*xj+yj*yj+zj*zj
105 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
106 if (sss.lt.1.0d0) then
108 eps0ij=eps(itypi,itypj)
113 evdw=evdw+(1.0d0-sss)*evdwij
115 C Calculate the components of the gradient in DC and X
117 fac=-rrij*(e1+evdwij)*(1.0d0-sss)
122 gvdwx(k,i)=gvdwx(k,i)-gg(k)
123 gvdwx(k,j)=gvdwx(k,j)+gg(k)
124 gvdwc(k,i)=gvdwc(k,i)-gg(k)
125 gvdwc(k,j)=gvdwc(k,j)+gg(k)
133 gvdwc(j,i)=expon*gvdwc(j,i)
134 gvdwx(j,i)=expon*gvdwx(j,i)
137 C******************************************************************************
141 C To save time, the factor of EXPON has been extracted from ALL components
142 C of GVDWC and GRADX. Remember to multiply them by this factor before further
145 C******************************************************************************
148 C-----------------------------------------------------------------------
149 subroutine elj_short(evdw)
151 C This subroutine calculates the interaction energy of nonbonded side chains
152 C assuming the LJ potential of interaction.
154 implicit real*8 (a-h,o-z)
156 parameter (accur=1.0d-10)
159 include 'COMMON.LOCAL'
160 include 'COMMON.CHAIN'
161 include 'COMMON.DERIV'
162 include 'COMMON.INTERACT'
163 include 'COMMON.TORSION'
164 include 'COMMON.SBRIDGE'
165 include 'COMMON.NAMES'
166 include 'COMMON.IOUNITS'
167 c include 'COMMON.CONTACTS'
169 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
173 if (itypi.eq.ntyp1) cycle
174 itypi1=iabs(itype(i+1))
181 C Calculate SC interaction energy.
184 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
185 cd & 'iend=',iend(i,iint)
186 do j=istart(i,iint),iend(i,iint)
188 if (itypj.eq.ntyp1) cycle
192 C Change 12/1/95 to calculate four-body interactions
193 rij=xj*xj+yj*yj+zj*zj
194 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
195 if (sss.gt.0.0d0) then
197 eps0ij=eps(itypi,itypj)
204 C Calculate the components of the gradient in DC and X
206 fac=-rrij*(e1+evdwij)*sss
211 gvdwx(k,i)=gvdwx(k,i)-gg(k)
212 gvdwx(k,j)=gvdwx(k,j)+gg(k)
213 gvdwc(k,i)=gvdwc(k,i)-gg(k)
214 gvdwc(k,j)=gvdwc(k,j)+gg(k)
222 gvdwc(j,i)=expon*gvdwc(j,i)
223 gvdwx(j,i)=expon*gvdwx(j,i)
226 C******************************************************************************
230 C To save time, the factor of EXPON has been extracted from ALL components
231 C of GVDWC and GRADX. Remember to multiply them by this factor before further
234 C******************************************************************************
237 C-----------------------------------------------------------------------------
238 subroutine eljk_long(evdw)
240 C This subroutine calculates the interaction energy of nonbonded side chains
241 C assuming the LJK potential of interaction.
243 implicit real*8 (a-h,o-z)
247 include 'COMMON.LOCAL'
248 include 'COMMON.CHAIN'
249 include 'COMMON.DERIV'
250 include 'COMMON.INTERACT'
251 include 'COMMON.IOUNITS'
252 include 'COMMON.NAMES'
255 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
259 if (itypi.eq.ntyp1) cycle
260 itypi1=iabs(itype(i+1))
265 C Calculate SC interaction energy.
268 do j=istart(i,iint),iend(i,iint)
270 if (itypj.eq.ntyp1) cycle
274 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
276 e_augm=augm(itypi,itypj)*fac_augm
279 sss=sscale(rij/sigma(itypi,itypj))
280 if (sss.lt.1.0d0) then
281 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
282 fac=r_shift_inv**expon
286 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
287 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
288 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
289 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
290 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
291 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
292 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
293 evdw=evdw+(1.0d0-sss)*evdwij
295 C Calculate the components of the gradient in DC and X
297 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
303 gvdwx(k,i)=gvdwx(k,i)-gg(k)
304 gvdwx(k,j)=gvdwx(k,j)+gg(k)
305 gvdwc(k,i)=gvdwc(k,i)-gg(k)
306 gvdwc(k,j)=gvdwc(k,j)+gg(k)
314 gvdwc(j,i)=expon*gvdwc(j,i)
315 gvdwx(j,i)=expon*gvdwx(j,i)
320 C-----------------------------------------------------------------------------
321 subroutine eljk_short(evdw)
323 C This subroutine calculates the interaction energy of nonbonded side chains
324 C assuming the LJK potential of interaction.
326 implicit real*8 (a-h,o-z)
330 include 'COMMON.LOCAL'
331 include 'COMMON.CHAIN'
332 include 'COMMON.DERIV'
333 include 'COMMON.INTERACT'
334 include 'COMMON.IOUNITS'
335 include 'COMMON.NAMES'
338 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
342 if (itypi.eq.ntyp1) cycle
343 itypi1=iabs(itype(i+1))
348 C Calculate SC interaction energy.
351 do j=istart(i,iint),iend(i,iint)
353 if (itypj.eq.ntyp1) cycle
357 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
359 e_augm=augm(itypi,itypj)*fac_augm
362 sss=sscale(rij/sigma(itypi,itypj))
363 if (sss.gt.0.0d0) then
364 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
365 fac=r_shift_inv**expon
369 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
370 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
371 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
372 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
373 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
374 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
375 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
378 C Calculate the components of the gradient in DC and X
380 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
386 gvdwx(k,i)=gvdwx(k,i)-gg(k)
387 gvdwx(k,j)=gvdwx(k,j)+gg(k)
388 gvdwc(k,i)=gvdwc(k,i)-gg(k)
389 gvdwc(k,j)=gvdwc(k,j)+gg(k)
397 gvdwc(j,i)=expon*gvdwc(j,i)
398 gvdwx(j,i)=expon*gvdwx(j,i)
403 C-----------------------------------------------------------------------------
404 subroutine ebp_long(evdw)
406 C This subroutine calculates the interaction energy of nonbonded side chains
407 C assuming the Berne-Pechukas potential of interaction.
409 implicit real*8 (a-h,o-z)
413 include 'COMMON.LOCAL'
414 include 'COMMON.CHAIN'
415 include 'COMMON.DERIV'
416 include 'COMMON.NAMES'
417 include 'COMMON.INTERACT'
418 include 'COMMON.IOUNITS'
419 include 'COMMON.CALC'
421 c double precision rrsave(maxdim)
424 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
426 c if (icall.eq.0) then
434 if (itypi.eq.ntyp1) cycle
435 itypi1=iabs(itype(i+1))
439 dxi=dc_norm(1,nres+i)
440 dyi=dc_norm(2,nres+i)
441 dzi=dc_norm(3,nres+i)
442 c dsci_inv=dsc_inv(itypi)
443 dsci_inv=vbld_inv(i+nres)
445 C Calculate SC interaction energy.
448 do j=istart(i,iint),iend(i,iint)
451 if (itypj.eq.ntyp1) cycle
452 c dscj_inv=dsc_inv(itypj)
453 dscj_inv=vbld_inv(j+nres)
454 chi1=chi(itypi,itypj)
455 chi2=chi(itypj,itypi)
462 alf12=0.5D0*(alf1+alf2)
466 dxj=dc_norm(1,nres+j)
467 dyj=dc_norm(2,nres+j)
468 dzj=dc_norm(3,nres+j)
469 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
471 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
473 if (sss.lt.1.0d0) then
475 C Calculate the angle-dependent terms of energy & contributions to derivatives.
477 C Calculate whole angle-dependent part of epsilon and contributions
479 fac=(rrij*sigsq)**expon2
482 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
483 eps2der=evdwij*eps3rt
484 eps3der=evdwij*eps2rt
485 evdwij=evdwij*eps2rt*eps3rt
486 evdw=evdw+evdwij*(1.0d0-sss)
488 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
490 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
491 cd & restyp(itypi),i,restyp(itypj),j,
492 cd & epsi,sigm,chi1,chi2,chip1,chip2,
493 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
494 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
497 C Calculate gradient components.
498 e1=e1*eps1*eps2rt**2*eps3rt**2
499 fac=-expon*(e1+evdwij)
502 C Calculate radial part of the gradient
506 C Calculate the angular part of the gradient and sum add the contributions
507 C to the appropriate components of the Cartesian gradient.
508 call sc_grad_scale(1.0d0-sss)
516 C-----------------------------------------------------------------------------
517 subroutine ebp_short(evdw)
519 C This subroutine calculates the interaction energy of nonbonded side chains
520 C assuming the Berne-Pechukas potential of interaction.
522 implicit real*8 (a-h,o-z)
526 include 'COMMON.LOCAL'
527 include 'COMMON.CHAIN'
528 include 'COMMON.DERIV'
529 include 'COMMON.NAMES'
530 include 'COMMON.INTERACT'
531 include 'COMMON.IOUNITS'
532 include 'COMMON.CALC'
534 c double precision rrsave(maxdim)
537 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
539 c if (icall.eq.0) then
547 if (itypi.eq.ntyp1) cycle
548 itypi1=iabs(itype(i+1))
552 dxi=dc_norm(1,nres+i)
553 dyi=dc_norm(2,nres+i)
554 dzi=dc_norm(3,nres+i)
555 c dsci_inv=dsc_inv(itypi)
556 dsci_inv=vbld_inv(i+nres)
558 C Calculate SC interaction energy.
561 do j=istart(i,iint),iend(i,iint)
564 if (itypj.eq.ntyp1) cycle
565 c dscj_inv=dsc_inv(itypj)
566 dscj_inv=vbld_inv(j+nres)
567 chi1=chi(itypi,itypj)
568 chi2=chi(itypj,itypi)
575 alf12=0.5D0*(alf1+alf2)
579 dxj=dc_norm(1,nres+j)
580 dyj=dc_norm(2,nres+j)
581 dzj=dc_norm(3,nres+j)
582 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
584 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
586 if (sss.gt.0.0d0) then
588 C Calculate the angle-dependent terms of energy & contributions to derivatives.
590 C Calculate whole angle-dependent part of epsilon and contributions
592 fac=(rrij*sigsq)**expon2
595 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
596 eps2der=evdwij*eps3rt
597 eps3der=evdwij*eps2rt
598 evdwij=evdwij*eps2rt*eps3rt
601 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
603 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
604 cd & restyp(itypi),i,restyp(itypj),j,
605 cd & epsi,sigm,chi1,chi2,chip1,chip2,
606 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
607 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
610 C Calculate gradient components.
611 e1=e1*eps1*eps2rt**2*eps3rt**2
612 fac=-expon*(e1+evdwij)
615 C Calculate radial part of the gradient
619 C Calculate the angular part of the gradient and sum add the contributions
620 C to the appropriate components of the Cartesian gradient.
621 call sc_grad_scale(sss)
629 C-----------------------------------------------------------------------------
630 subroutine egb_long(evdw)
632 C This subroutine calculates the interaction energy of nonbonded side chains
633 C assuming the Gay-Berne potential of interaction.
635 implicit real*8 (a-h,o-z)
639 include 'COMMON.LOCAL'
640 include 'COMMON.CHAIN'
641 include 'COMMON.DERIV'
642 include 'COMMON.NAMES'
643 include 'COMMON.INTERACT'
644 include 'COMMON.IOUNITS'
645 include 'COMMON.CALC'
646 include 'COMMON.CONTROL'
648 integer xshift,yshift,zshift
650 ccccc energy_dec=.false.
651 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
654 c if (icall.eq.0) lprn=.false.
658 if (itypi.eq.ntyp1) cycle
659 itypi1=iabs(itype(i+1))
664 if (xi.lt.0) xi=xi+boxxsize
666 if (yi.lt.0) yi=yi+boxysize
668 if (zi.lt.0) zi=zi+boxzsize
669 dxi=dc_norm(1,nres+i)
670 dyi=dc_norm(2,nres+i)
671 dzi=dc_norm(3,nres+i)
672 c dsci_inv=dsc_inv(itypi)
673 dsci_inv=vbld_inv(i+nres)
674 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
675 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
677 C Calculate SC interaction energy.
680 do j=istart(i,iint),iend(i,iint)
683 if (itypj.eq.ntyp1) cycle
684 c dscj_inv=dsc_inv(itypj)
685 dscj_inv=vbld_inv(j+nres)
686 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
687 c & 1.0d0/vbld(j+nres)
688 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
689 sig0ij=sigma(itypi,itypj)
690 chi1=chi(itypi,itypj)
691 chi2=chi(itypj,itypi)
698 alf12=0.5D0*(alf1+alf2)
703 if (xj.lt.0) xj=xj+boxxsize
705 if (yj.lt.0) yj=yj+boxysize
707 if (zj.lt.0) zj=zj+boxzsize
708 if ((zj.gt.bordlipbot)
709 &.and.(zj.lt.bordliptop)) then
710 C the energy transfer exist
711 if (zj.lt.buflipbot) then
712 C what fraction I am in
714 & ((positi-bordlipbot)/lipbufthick)
715 C lipbufthick is thickenes of lipid buffore
716 sslipj=sscalelip(fracinbuf)
717 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
718 elseif (zi.gt.bufliptop) then
719 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
720 sslipj=sscalelip(fracinbuf)
721 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
730 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
731 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
732 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
733 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
735 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
743 xj=xj_safe+xshift*boxxsize
744 yj=yj_safe+yshift*boxysize
745 zj=zj_safe+zshift*boxzsize
746 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
747 if(dist_temp.lt.dist_init) then
757 if (subchap.eq.1) then
766 dxj=dc_norm(1,nres+j)
767 dyj=dc_norm(2,nres+j)
768 dzj=dc_norm(3,nres+j)
769 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
771 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
772 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
773 if (sss.lt.1.0d0) then
775 C Calculate angle-dependent terms of energy and contributions to their
779 sig=sig0ij*dsqrt(sigsq)
780 rij_shift=1.0D0/rij-sig+sig0ij
781 c for diagnostics; uncomment
782 c rij_shift=1.2*sig0ij
783 C I hate to put IF's in the loops, but here don't have another choice!!!!
784 if (rij_shift.le.0.0D0) then
786 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
787 cd & restyp(itypi),i,restyp(itypj),j,
788 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
792 c---------------------------------------------------------------
793 rij_shift=1.0D0/rij_shift
797 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
798 eps2der=evdwij*eps3rt
799 eps3der=evdwij*eps2rt
800 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
801 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
802 evdwij=evdwij*eps2rt*eps3rt
803 evdw=evdw+evdwij*(1.0d0-sss)
805 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
807 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
808 & restyp(itypi),i,restyp(itypj),j,
809 & epsi,sigm,chi1,chi2,chip1,chip2,
810 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
811 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
815 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
818 C Calculate gradient components.
819 e1=e1*eps1*eps2rt**2*eps3rt**2
820 fac=-expon*(e1+evdwij)*rij_shift
823 fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij
825 C Calculate the radial part of the gradient
829 gg_lipi(3)=ssgradlipi*evdwij
830 gg_lipj(3)=ssgradlipj*evdwij
831 C Calculate angular part of the gradient.
832 call sc_grad_scale(1.0d0-sss)
837 c write (iout,*) "Number of loop steps in EGB:",ind
838 cccc energy_dec=.false.
841 C-----------------------------------------------------------------------------
842 subroutine egb_short(evdw)
844 C This subroutine calculates the interaction energy of nonbonded side chains
845 C assuming the Gay-Berne potential of interaction.
847 implicit real*8 (a-h,o-z)
851 include 'COMMON.LOCAL'
852 include 'COMMON.CHAIN'
853 include 'COMMON.DERIV'
854 include 'COMMON.NAMES'
855 include 'COMMON.INTERACT'
856 include 'COMMON.IOUNITS'
857 include 'COMMON.CALC'
858 include 'COMMON.CONTROL'
860 integer xshift,yshift,zshift
862 ccccc energy_dec=.false.
863 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
866 c if (icall.eq.0) lprn=.false.
870 if (itypi.eq.ntyp1) cycle
871 itypi1=iabs(itype(i+1))
876 if (xi.lt.0) xi=xi+boxxsize
878 if (yi.lt.0) yi=yi+boxysize
880 if (zi.lt.0) zi=zi+boxzsize
881 dxi=dc_norm(1,nres+i)
882 dyi=dc_norm(2,nres+i)
883 dzi=dc_norm(3,nres+i)
884 c dsci_inv=dsc_inv(itypi)
885 dsci_inv=vbld_inv(i+nres)
886 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
887 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
889 C Calculate SC interaction energy.
892 do j=istart(i,iint),iend(i,iint)
895 if (itypj.eq.ntyp1) cycle
896 c dscj_inv=dsc_inv(itypj)
897 dscj_inv=vbld_inv(j+nres)
898 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
899 c & 1.0d0/vbld(j+nres)
900 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
901 sig0ij=sigma(itypi,itypj)
902 chi1=chi(itypi,itypj)
903 chi2=chi(itypj,itypi)
910 alf12=0.5D0*(alf1+alf2)
915 if (xj.lt.0) xj=xj+boxxsize
917 if (yj.lt.0) yj=yj+boxysize
919 if (zj.lt.0) zj=zj+boxzsize
920 if ((zj.gt.bordlipbot)
921 &.and.(zj.lt.bordliptop)) then
922 C the energy transfer exist
923 if (zj.lt.buflipbot) then
924 C what fraction I am in
926 & ((positi-bordlipbot)/lipbufthick)
927 C lipbufthick is thickenes of lipid buffore
928 sslipj=sscalelip(fracinbuf)
929 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
930 elseif (zi.gt.bufliptop) then
931 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
932 sslipj=sscalelip(fracinbuf)
933 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
942 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
943 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
944 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
945 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
946 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
954 xj=xj_safe+xshift*boxxsize
955 yj=yj_safe+yshift*boxysize
956 zj=zj_safe+zshift*boxzsize
957 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
958 if(dist_temp.lt.dist_init) then
968 if (subchap.eq.1) then
977 dxj=dc_norm(1,nres+j)
978 dyj=dc_norm(2,nres+j)
979 dzj=dc_norm(3,nres+j)
980 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
982 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
983 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
984 if (sss.gt.0.0d0) then
986 C Calculate angle-dependent terms of energy and contributions to their
990 sig=sig0ij*dsqrt(sigsq)
991 rij_shift=1.0D0/rij-sig+sig0ij
992 c for diagnostics; uncomment
993 c rij_shift=1.2*sig0ij
994 C I hate to put IF's in the loops, but here don't have another choice!!!!
995 if (rij_shift.le.0.0D0) then
997 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
998 cd & restyp(itypi),i,restyp(itypj),j,
999 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1003 c---------------------------------------------------------------
1004 rij_shift=1.0D0/rij_shift
1005 fac=rij_shift**expon
1008 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1009 eps2der=evdwij*eps3rt
1010 eps3der=evdwij*eps2rt
1011 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1012 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1013 evdwij=evdwij*eps2rt*eps3rt
1014 evdw=evdw+evdwij*sss
1016 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1018 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1019 & restyp(itypi),i,restyp(itypj),j,
1020 & epsi,sigm,chi1,chi2,chip1,chip2,
1021 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1022 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1026 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1029 C Calculate gradient components.
1030 e1=e1*eps1*eps2rt**2*eps3rt**2
1031 fac=-expon*(e1+evdwij)*rij_shift
1034 fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij
1036 C Calculate the radial part of the gradient
1040 gg_lipi(3)=ssgradlipi*evdwij
1041 gg_lipj(3)=ssgradlipj*evdwij
1042 C Calculate angular part of the gradient.
1043 call sc_grad_scale(sss)
1048 c write (iout,*) "Number of loop steps in EGB:",ind
1049 cccc energy_dec=.false.
1052 C-----------------------------------------------------------------------------
1053 subroutine egbv_long(evdw)
1055 C This subroutine calculates the interaction energy of nonbonded side chains
1056 C assuming the Gay-Berne-Vorobjev potential of interaction.
1058 implicit real*8 (a-h,o-z)
1059 include 'DIMENSIONS'
1060 include 'COMMON.GEO'
1061 include 'COMMON.VAR'
1062 include 'COMMON.LOCAL'
1063 include 'COMMON.CHAIN'
1064 include 'COMMON.DERIV'
1065 include 'COMMON.NAMES'
1066 include 'COMMON.INTERACT'
1067 include 'COMMON.IOUNITS'
1068 include 'COMMON.CALC'
1069 common /srutu/ icall
1072 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1075 c if (icall.eq.0) lprn=.true.
1077 do i=iatsc_s,iatsc_e
1078 itypi=iabs(itype(i))
1079 if (itypi.eq.ntyp1) cycle
1080 itypi1=iabs(itype(i+1))
1084 dxi=dc_norm(1,nres+i)
1085 dyi=dc_norm(2,nres+i)
1086 dzi=dc_norm(3,nres+i)
1087 c dsci_inv=dsc_inv(itypi)
1088 dsci_inv=vbld_inv(i+nres)
1090 C Calculate SC interaction energy.
1092 do iint=1,nint_gr(i)
1093 do j=istart(i,iint),iend(i,iint)
1095 itypj=iabs(itype(j))
1096 if (itypj.eq.ntyp1) cycle
1097 c dscj_inv=dsc_inv(itypj)
1098 dscj_inv=vbld_inv(j+nres)
1099 sig0ij=sigma(itypi,itypj)
1100 r0ij=r0(itypi,itypj)
1101 chi1=chi(itypi,itypj)
1102 chi2=chi(itypj,itypi)
1109 alf12=0.5D0*(alf1+alf2)
1113 dxj=dc_norm(1,nres+j)
1114 dyj=dc_norm(2,nres+j)
1115 dzj=dc_norm(3,nres+j)
1116 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1119 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1121 if (sss.lt.1.0d0) then
1123 C Calculate angle-dependent terms of energy and contributions to their
1127 sig=sig0ij*dsqrt(sigsq)
1128 rij_shift=1.0D0/rij-sig+r0ij
1129 C I hate to put IF's in the loops, but here don't have another choice!!!!
1130 if (rij_shift.le.0.0D0) then
1135 c---------------------------------------------------------------
1136 rij_shift=1.0D0/rij_shift
1137 fac=rij_shift**expon
1140 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1141 eps2der=evdwij*eps3rt
1142 eps3der=evdwij*eps2rt
1143 fac_augm=rrij**expon
1144 e_augm=augm(itypi,itypj)*fac_augm
1145 evdwij=evdwij*eps2rt*eps3rt
1146 evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
1148 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1150 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1151 & restyp(itypi),i,restyp(itypj),j,
1152 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1153 & chi1,chi2,chip1,chip2,
1154 & eps1,eps2rt**2,eps3rt**2,
1155 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1158 C Calculate gradient components.
1159 e1=e1*eps1*eps2rt**2*eps3rt**2
1160 fac=-expon*(e1+evdwij)*rij_shift
1162 fac=rij*fac-2*expon*rrij*e_augm
1163 C Calculate the radial part of the gradient
1167 C Calculate angular part of the gradient.
1168 call sc_grad_scale(1.0d0-sss)
1174 C-----------------------------------------------------------------------------
1175 subroutine egbv_short(evdw)
1177 C This subroutine calculates the interaction energy of nonbonded side chains
1178 C assuming the Gay-Berne-Vorobjev potential of interaction.
1180 implicit real*8 (a-h,o-z)
1181 include 'DIMENSIONS'
1182 include 'COMMON.GEO'
1183 include 'COMMON.VAR'
1184 include 'COMMON.LOCAL'
1185 include 'COMMON.CHAIN'
1186 include 'COMMON.DERIV'
1187 include 'COMMON.NAMES'
1188 include 'COMMON.INTERACT'
1189 include 'COMMON.IOUNITS'
1190 include 'COMMON.CALC'
1191 common /srutu/ icall
1194 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1197 c if (icall.eq.0) lprn=.true.
1199 do i=iatsc_s,iatsc_e
1200 itypi=iabs(itype(i))
1201 if (itypi.eq.ntyp1) cycle
1202 itypi1=iabs(itype(i+1))
1206 dxi=dc_norm(1,nres+i)
1207 dyi=dc_norm(2,nres+i)
1208 dzi=dc_norm(3,nres+i)
1209 c dsci_inv=dsc_inv(itypi)
1210 dsci_inv=vbld_inv(i+nres)
1212 C Calculate SC interaction energy.
1214 do iint=1,nint_gr(i)
1215 do j=istart(i,iint),iend(i,iint)
1217 itypj=iabs(itype(j))
1218 if (itypj.eq.ntyp1) cycle
1219 c dscj_inv=dsc_inv(itypj)
1220 dscj_inv=vbld_inv(j+nres)
1221 sig0ij=sigma(itypi,itypj)
1222 r0ij=r0(itypi,itypj)
1223 chi1=chi(itypi,itypj)
1224 chi2=chi(itypj,itypi)
1231 alf12=0.5D0*(alf1+alf2)
1235 dxj=dc_norm(1,nres+j)
1236 dyj=dc_norm(2,nres+j)
1237 dzj=dc_norm(3,nres+j)
1238 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1241 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1243 if (sss.gt.0.0d0) then
1245 C Calculate angle-dependent terms of energy and contributions to their
1249 sig=sig0ij*dsqrt(sigsq)
1250 rij_shift=1.0D0/rij-sig+r0ij
1251 C I hate to put IF's in the loops, but here don't have another choice!!!!
1252 if (rij_shift.le.0.0D0) then
1257 c---------------------------------------------------------------
1258 rij_shift=1.0D0/rij_shift
1259 fac=rij_shift**expon
1262 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1263 eps2der=evdwij*eps3rt
1264 eps3der=evdwij*eps2rt
1265 fac_augm=rrij**expon
1266 e_augm=augm(itypi,itypj)*fac_augm
1267 evdwij=evdwij*eps2rt*eps3rt
1268 evdw=evdw+(evdwij+e_augm)*sss
1270 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1272 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1273 & restyp(itypi),i,restyp(itypj),j,
1274 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1275 & chi1,chi2,chip1,chip2,
1276 & eps1,eps2rt**2,eps3rt**2,
1277 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1280 C Calculate gradient components.
1281 e1=e1*eps1*eps2rt**2*eps3rt**2
1282 fac=-expon*(e1+evdwij)*rij_shift
1284 fac=rij*fac-2*expon*rrij*e_augm
1285 C Calculate the radial part of the gradient
1289 C Calculate angular part of the gradient.
1290 call sc_grad_scale(sss)
1296 C----------------------------------------------------------------------------
1297 subroutine sc_grad_scale(scalfac)
1298 implicit real*8 (a-h,o-z)
1299 include 'DIMENSIONS'
1300 include 'COMMON.CHAIN'
1301 include 'COMMON.DERIV'
1302 include 'COMMON.CALC'
1303 include 'COMMON.IOUNITS'
1304 double precision dcosom1(3),dcosom2(3)
1305 double precision scalfac
1306 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1307 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1308 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1309 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1313 c eom12=evdwij*eps1_om12
1315 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1316 c & " sigder",sigder
1317 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1318 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1320 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1321 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1324 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1326 c write (iout,*) "gg",(gg(k),k=1,3)
1328 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
1329 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1330 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1331 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
1332 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1333 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1334 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1335 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1336 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1337 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1340 C Calculate the components of the gradient in DC and X
1343 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
1344 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
1348 C--------------------------------------------------------------------------
1349 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1351 C This subroutine calculates the average interaction energy and its gradient
1352 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1353 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1354 C The potential depends both on the distance of peptide-group centers and on
1355 C the orientation of the CA-CA virtual bonds.
1357 implicit real*8 (a-h,o-z)
1361 include 'DIMENSIONS'
1362 include 'COMMON.CONTROL'
1363 include 'COMMON.SETUP'
1364 include 'COMMON.IOUNITS'
1365 include 'COMMON.GEO'
1366 include 'COMMON.VAR'
1367 include 'COMMON.LOCAL'
1368 include 'COMMON.CHAIN'
1369 include 'COMMON.DERIV'
1370 include 'COMMON.INTERACT'
1372 include 'COMMON.CONTACTS'
1373 include 'COMMON.CONTMAT'
1375 include 'COMMON.CORRMAT'
1376 include 'COMMON.TORSION'
1377 include 'COMMON.VECTORS'
1378 include 'COMMON.FFIELD'
1379 include 'COMMON.TIME1'
1380 include 'COMMON.SHIELD'
1381 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1382 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1383 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1384 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1385 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1386 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1388 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1390 double precision scal_el /1.0d0/
1392 double precision scal_el /0.5d0/
1395 C 13-go grudnia roku pamietnego...
1396 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1397 & 0.0d0,1.0d0,0.0d0,
1398 & 0.0d0,0.0d0,1.0d0/
1399 cd write(iout,*) 'In EELEC'
1401 cd write(iout,*) 'Type',i
1402 cd write(iout,*) 'B1',B1(:,i)
1403 cd write(iout,*) 'B2',B2(:,i)
1404 cd write(iout,*) 'CC',CC(:,:,i)
1405 cd write(iout,*) 'DD',DD(:,:,i)
1406 cd write(iout,*) 'EE',EE(:,:,i)
1408 cd call check_vecgrad
1410 C print *,"WCHODZE3"
1411 if (icheckgrad.eq.1) then
1413 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1415 dc_norm(k,i)=dc(k,i)*fac
1417 c write (iout,*) 'i',i,' fac',fac
1420 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1421 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1422 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1423 c call vec_and_deriv
1429 time_mat=time_mat+MPI_Wtime()-time01
1433 cd write (iout,*) 'i=',i
1435 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1438 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1439 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1454 cd print '(a)','Enter EELEC'
1455 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1457 gel_loc_loc(i)=0.0d0
1462 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1464 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1466 do i=iturn3_start,iturn3_end
1467 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1468 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1469 C & .or. itype(i-1).eq.ntyp1
1470 C & .or. itype(i+4).eq.ntyp1
1475 dx_normi=dc_norm(1,i)
1476 dy_normi=dc_norm(2,i)
1477 dz_normi=dc_norm(3,i)
1478 xmedi=c(1,i)+0.5d0*dxi
1479 ymedi=c(2,i)+0.5d0*dyi
1480 zmedi=c(3,i)+0.5d0*dzi
1481 xmedi=mod(xmedi,boxxsize)
1482 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1483 ymedi=mod(ymedi,boxysize)
1484 if (ymedi.lt.0) ymedi=ymedi+boxysize
1485 zmedi=mod(zmedi,boxzsize)
1486 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1488 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1489 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1491 num_cont_hb(i)=num_conti
1494 do i=iturn4_start,iturn4_end
1495 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1496 & .or. itype(i+3).eq.ntyp1
1497 & .or. itype(i+4).eq.ntyp1
1498 C & .or. itype(i+5).eq.ntyp1
1499 C & .or. itype(i-1).eq.ntyp1
1504 dx_normi=dc_norm(1,i)
1505 dy_normi=dc_norm(2,i)
1506 dz_normi=dc_norm(3,i)
1507 xmedi=c(1,i)+0.5d0*dxi
1508 ymedi=c(2,i)+0.5d0*dyi
1509 zmedi=c(3,i)+0.5d0*dzi
1510 xmedi=mod(xmedi,boxxsize)
1511 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1512 ymedi=mod(ymedi,boxysize)
1513 if (ymedi.lt.0) ymedi=ymedi+boxysize
1514 zmedi=mod(zmedi,boxzsize)
1515 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1517 num_conti=num_cont_hb(i)
1519 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1520 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1521 & call eturn4(i,eello_turn4)
1523 num_cont_hb(i)=num_conti
1527 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1529 do i=iatel_s,iatel_e
1530 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1531 C & .or. itype(i+2).eq.ntyp1
1532 C & .or. itype(i-1).eq.ntyp1
1537 dx_normi=dc_norm(1,i)
1538 dy_normi=dc_norm(2,i)
1539 dz_normi=dc_norm(3,i)
1540 xmedi=c(1,i)+0.5d0*dxi
1541 ymedi=c(2,i)+0.5d0*dyi
1542 zmedi=c(3,i)+0.5d0*dzi
1543 xmedi=mod(xmedi,boxxsize)
1544 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1545 ymedi=mod(ymedi,boxysize)
1546 if (ymedi.lt.0) ymedi=ymedi+boxysize
1547 zmedi=mod(zmedi,boxzsize)
1548 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1549 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
1551 num_conti=num_cont_hb(i)
1553 do j=ielstart(i),ielend(i)
1554 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1555 C & .or.itype(j+2).eq.ntyp1
1556 C & .or.itype(j-1).eq.ntyp1
1558 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1561 num_cont_hb(i)=num_conti
1564 c write (iout,*) "Number of loop steps in EELEC:",ind
1566 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1567 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1569 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1570 ccc eel_loc=eel_loc+eello_turn3
1571 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1574 C-------------------------------------------------------------------------------
1575 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1576 implicit real*8 (a-h,o-z)
1577 include 'DIMENSIONS'
1581 include 'COMMON.CONTROL'
1582 include 'COMMON.IOUNITS'
1583 include 'COMMON.GEO'
1584 include 'COMMON.VAR'
1585 include 'COMMON.LOCAL'
1586 include 'COMMON.CHAIN'
1587 include 'COMMON.DERIV'
1588 include 'COMMON.INTERACT'
1590 include 'COMMON.CONTACTS'
1591 include 'COMMON.CONTMAT'
1593 include 'COMMON.CORRMAT'
1594 include 'COMMON.TORSION'
1595 include 'COMMON.VECTORS'
1596 include 'COMMON.FFIELD'
1597 include 'COMMON.TIME1'
1598 include 'COMMON.SHIELD'
1599 integer xshift,yshift,zshift
1600 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1601 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1602 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1603 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
1604 & gmuij2(4),gmuji2(4)
1605 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1606 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1608 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1610 double precision scal_el /1.0d0/
1612 double precision scal_el /0.5d0/
1615 C 13-go grudnia roku pamietnego...
1616 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1617 & 0.0d0,1.0d0,0.0d0,
1618 & 0.0d0,0.0d0,1.0d0/
1619 c time00=MPI_Wtime()
1620 cd write (iout,*) "eelecij",i,j
1621 C print *,"WCHODZE2"
1625 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1626 aaa=app(iteli,itelj)
1627 bbb=bpp(iteli,itelj)
1628 ael6i=ael6(iteli,itelj)
1629 ael3i=ael3(iteli,itelj)
1633 dx_normj=dc_norm(1,j)
1634 dy_normj=dc_norm(2,j)
1635 dz_normj=dc_norm(3,j)
1640 if (xj.lt.0) xj=xj+boxxsize
1642 if (yj.lt.0) yj=yj+boxysize
1644 if (zj.lt.0) zj=zj+boxzsize
1645 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1653 xj=xj_safe+xshift*boxxsize
1654 yj=yj_safe+yshift*boxysize
1655 zj=zj_safe+zshift*boxzsize
1656 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1657 if(dist_temp.lt.dist_init) then
1667 if (isubchap.eq.1) then
1677 rij=xj*xj+yj*yj+zj*zj
1681 c For extracting the short-range part of Evdwpp
1682 sss=sscale(rij/rpp(iteli,itelj))
1683 sssgrad=sscagrad(rij/rpp(iteli,itelj))
1686 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1687 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1688 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1689 fac=cosa-3.0D0*cosb*cosg
1691 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1692 if (j.eq.i+2) ev1=scal_el*ev1
1697 if (shield_mode.eq.0) then
1701 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1703 el1=el1*fac_shield(i)**2*fac_shield(j)**2
1704 el2=el2*fac_shield(i)**2*fac_shield(j)**2
1706 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1707 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1709 evdw1=evdw1+evdwij*(1.0d0-sss)
1710 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1711 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1712 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1713 cd & xmedi,ymedi,zmedi,xj,yj,zj
1715 if (energy_dec) then
1716 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1717 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1721 C Calculate contributions to the Cartesian gradient.
1724 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1725 facel=-3*rrmij*(el1+eesij)
1731 * Radial derivatives. First process both termini of the fragment (i,j)
1736 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1737 & (shield_mode.gt.0)) then
1739 do ilist=1,ishield_list(i)
1740 iresshield=shield_list(ilist,i)
1742 rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
1744 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1746 & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
1747 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1748 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1749 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1750 C if (iresshield.gt.i) then
1751 C do ishi=i+1,iresshield-1
1752 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1753 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1757 C do ishi=iresshield,i
1758 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1759 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1765 do ilist=1,ishield_list(j)
1766 iresshield=shield_list(ilist,j)
1768 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1770 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1772 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
1773 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1778 gshieldc(k,i)=gshieldc(k,i)+
1779 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
1780 gshieldc(k,j)=gshieldc(k,j)+
1781 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
1782 gshieldc(k,i-1)=gshieldc(k,i-1)+
1783 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
1784 gshieldc(k,j-1)=gshieldc(k,j-1)+
1785 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
1791 c ghalf=0.5D0*ggg(k)
1792 c gelc(k,i)=gelc(k,i)+ghalf
1793 c gelc(k,j)=gelc(k,j)+ghalf
1795 c 9/28/08 AL Gradient compotents will be summed only at the end
1797 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1798 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1801 * Loop over residues i+1 thru j-1.
1805 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1808 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1809 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1810 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1812 c ghalf=0.5D0*ggg(k)
1813 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1814 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1816 c 9/28/08 AL Gradient compotents will be summed only at the end
1818 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1819 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1822 * Loop over residues i+1 thru j-1.
1826 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1830 facvdw=ev1+evdwij*(1.0d0-sss)
1833 fac=-3*rrmij*(facvdw+facvdw+facel)
1838 * Radial derivatives. First process both termini of the fragment (i,j)
1844 c ghalf=0.5D0*ggg(k)
1845 c gelc(k,i)=gelc(k,i)+ghalf
1846 c gelc(k,j)=gelc(k,j)+ghalf
1848 c 9/28/08 AL Gradient compotents will be summed only at the end
1850 gelc_long(k,j)=gelc(k,j)+ggg(k)
1851 gelc_long(k,i)=gelc(k,i)-ggg(k)
1854 * Loop over residues i+1 thru j-1.
1858 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1861 c 9/28/08 AL Gradient compotents will be summed only at the end
1865 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1866 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1867 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1869 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1870 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1876 ecosa=2.0D0*fac3*fac1+fac4
1879 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1880 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1882 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1883 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1885 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1886 cd & (dcosg(k),k=1,3)
1888 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
1889 & fac_shield(i)**2*fac_shield(j)**2
1893 c ghalf=0.5D0*ggg(k)
1894 c gelc(k,i)=gelc(k,i)+ghalf
1895 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1896 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1897 c gelc(k,j)=gelc(k,j)+ghalf
1898 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1899 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1903 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1908 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1909 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
1910 & *fac_shield(i)**2*fac_shield(j)**2
1913 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1914 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
1915 & *fac_shield(i)**2*fac_shield(j)**2
1916 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1917 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1919 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1920 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1921 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1923 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1924 C energy of a peptide unit is assumed in the form of a second-order
1925 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1926 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1927 C are computed for EVERY pair of non-contiguous peptide groups.
1929 if (j.lt.nres-1) then
1940 muij(kkk)=mu(k,i)*mu(l,j)
1942 gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
1943 c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
1944 gmuij2(kkk)=gUb2(k,i)*mu(l,j)
1945 gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
1946 c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
1947 gmuji2(kkk)=mu(k,i)*gUb2(l,j)
1951 cd write (iout,*) 'EELEC: i',i,' j',j
1952 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1953 cd write(iout,*) 'muij',muij
1954 ury=scalar(uy(1,i),erij)
1955 urz=scalar(uz(1,i),erij)
1956 vry=scalar(uy(1,j),erij)
1957 vrz=scalar(uz(1,j),erij)
1958 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1959 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1960 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1961 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1962 fac=dsqrt(-ael6i)*r3ij
1967 cd write (iout,'(4i5,4f10.5)')
1968 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1969 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1970 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1971 cd & uy(:,j),uz(:,j)
1972 cd write (iout,'(4f10.5)')
1973 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1974 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1975 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
1976 cd write (iout,'(9f10.5/)')
1977 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1978 C Derivatives of the elements of A in virtual-bond vectors
1979 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1981 uryg(k,1)=scalar(erder(1,k),uy(1,i))
1982 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1983 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1984 urzg(k,1)=scalar(erder(1,k),uz(1,i))
1985 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1986 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1987 vryg(k,1)=scalar(erder(1,k),uy(1,j))
1988 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1989 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1990 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1991 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1992 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1994 C Compute radial contributions to the gradient
2012 C Add the contributions coming from er
2015 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
2016 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
2017 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
2018 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
2021 C Derivatives in DC(i)
2022 cgrad ghalf1=0.5d0*agg(k,1)
2023 cgrad ghalf2=0.5d0*agg(k,2)
2024 cgrad ghalf3=0.5d0*agg(k,3)
2025 cgrad ghalf4=0.5d0*agg(k,4)
2026 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2027 & -3.0d0*uryg(k,2)*vry)!+ghalf1
2028 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2029 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
2030 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2031 & -3.0d0*urzg(k,2)*vry)!+ghalf3
2032 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2033 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
2034 C Derivatives in DC(i+1)
2035 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2036 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2037 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2038 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2039 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2040 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2041 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2042 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2043 C Derivatives in DC(j)
2044 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2045 & -3.0d0*vryg(k,2)*ury)!+ghalf1
2046 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2047 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
2048 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2049 & -3.0d0*vryg(k,2)*urz)!+ghalf3
2050 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2051 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
2052 C Derivatives in DC(j+1) or DC(nres-1)
2053 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2054 & -3.0d0*vryg(k,3)*ury)
2055 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2056 & -3.0d0*vrzg(k,3)*ury)
2057 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2058 & -3.0d0*vryg(k,3)*urz)
2059 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2060 & -3.0d0*vrzg(k,3)*urz)
2061 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
2063 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
2076 aggi(k,l)=-aggi(k,l)
2077 aggi1(k,l)=-aggi1(k,l)
2078 aggj(k,l)=-aggj(k,l)
2079 aggj1(k,l)=-aggj1(k,l)
2082 if (j.lt.nres-1) then
2088 aggi(k,l)=-aggi(k,l)
2089 aggi1(k,l)=-aggi1(k,l)
2090 aggj(k,l)=-aggj(k,l)
2091 aggj1(k,l)=-aggj1(k,l)
2102 aggi(k,l)=-aggi(k,l)
2103 aggi1(k,l)=-aggi1(k,l)
2104 aggj(k,l)=-aggj(k,l)
2105 aggj1(k,l)=-aggj1(k,l)
2110 IF (wel_loc.gt.0.0d0) THEN
2111 C Contribution to the local-electrostatic energy coming from the i-j pair
2112 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2114 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2116 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2117 & 'eelloc',i,j,eel_loc_ij
2120 if (shield_mode.eq.0) then
2127 eel_loc_ij=eel_loc_ij
2128 & *fac_shield(i)*fac_shield(j)
2129 eel_loc=eel_loc+eel_loc_ij
2131 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2132 & (shield_mode.gt.0)) then
2135 do ilist=1,ishield_list(i)
2136 iresshield=shield_list(ilist,i)
2138 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2141 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2143 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2144 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2148 do ilist=1,ishield_list(j)
2149 iresshield=shield_list(ilist,j)
2151 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2154 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2156 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2157 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2164 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2165 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2166 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2167 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2168 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2169 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2170 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2171 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2176 geel_loc_ij=(a22*gmuij1(1)
2180 & *fac_shield(i)*fac_shield(j)
2181 c write(iout,*) "derivative over thatai"
2182 c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
2184 gloc(nphi+i,icg)=gloc(nphi+i,icg)+
2185 & geel_loc_ij*wel_loc
2186 c write(iout,*) "derivative over thatai-1"
2187 c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
2194 gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
2195 & geel_loc_ij*wel_loc
2196 & *fac_shield(i)*fac_shield(j)
2198 c Derivative over j residue
2199 geel_loc_ji=a22*gmuji1(1)
2203 c write(iout,*) "derivative over thataj"
2204 c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
2207 gloc(nphi+j,icg)=gloc(nphi+j,icg)+
2208 & geel_loc_ji*wel_loc
2209 & *fac_shield(i)*fac_shield(j)
2216 c write(iout,*) "derivative over thataj-1"
2217 c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
2219 gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
2220 & geel_loc_ji*wel_loc
2221 & *fac_shield(i)*fac_shield(j)
2223 cC Partial derivatives in virtual-bond dihedral angles gamma
2225 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2226 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2227 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2228 & *fac_shield(i)*fac_shield(j)
2230 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2231 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2232 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2233 & *fac_shield(i)*fac_shield(j)
2235 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2237 ggg(l)=(agg(l,1)*muij(1)+
2238 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2239 & *fac_shield(i)*fac_shield(j)
2241 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2242 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2243 cgrad ghalf=0.5d0*ggg(l)
2244 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
2245 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
2249 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2252 C Remaining derivatives of eello
2254 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2255 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2256 & *fac_shield(i)*fac_shield(j)
2258 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2259 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2260 & *fac_shield(i)*fac_shield(j)
2262 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2263 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2264 & *fac_shield(i)*fac_shield(j)
2266 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2267 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2268 & *fac_shield(i)*fac_shield(j)
2273 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2274 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
2275 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2276 & .and. num_conti.le.maxconts) then
2277 c write (iout,*) i,j," entered corr"
2279 C Calculate the contact function. The ith column of the array JCONT will
2280 C contain the numbers of atoms that make contacts with the atom I (of numbers
2281 C greater than I). The arrays FACONT and GACONT will contain the values of
2282 C the contact function and its derivative.
2283 c r0ij=1.02D0*rpp(iteli,itelj)
2284 c r0ij=1.11D0*rpp(iteli,itelj)
2285 r0ij=2.20D0*rpp(iteli,itelj)
2286 c r0ij=1.55D0*rpp(iteli,itelj)
2287 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2288 if (fcont.gt.0.0D0) then
2289 num_conti=num_conti+1
2290 if (num_conti.gt.maxconts) then
2291 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2292 & ' will skip next contacts for this conf.'
2294 jcont_hb(num_conti,i)=j
2295 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
2296 cd & " jcont_hb",jcont_hb(num_conti,i)
2297 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2298 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2299 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2301 d_cont(num_conti,i)=rij
2302 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2303 C --- Electrostatic-interaction matrix ---
2304 a_chuj(1,1,num_conti,i)=a22
2305 a_chuj(1,2,num_conti,i)=a23
2306 a_chuj(2,1,num_conti,i)=a32
2307 a_chuj(2,2,num_conti,i)=a33
2308 C --- Gradient of rij
2310 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2317 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2318 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2319 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2320 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2321 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2326 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2327 C Calculate contact energies
2329 wij=cosa-3.0D0*cosb*cosg
2332 c fac3=dsqrt(-ael6i)/r0ij**3
2333 fac3=dsqrt(-ael6i)*r3ij
2334 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2335 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2336 if (ees0tmp.gt.0) then
2337 ees0pij=dsqrt(ees0tmp)
2341 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2342 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2343 if (ees0tmp.gt.0) then
2344 ees0mij=dsqrt(ees0tmp)
2349 if (shield_mode.eq.0) then
2353 ees0plist(num_conti,i)=j
2354 C fac_shield(i)=0.4d0
2355 C fac_shield(j)=0.6d0
2357 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2358 & *fac_shield(i)*fac_shield(j)
2359 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2360 & *fac_shield(i)*fac_shield(j)
2362 C Diagnostics. Comment out or remove after debugging!
2363 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2364 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2365 c ees0m(num_conti,i)=0.0D0
2367 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2368 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2369 C Angular derivatives of the contact function
2370 ees0pij1=fac3/ees0pij
2371 ees0mij1=fac3/ees0mij
2372 fac3p=-3.0D0*fac3*rrmij
2373 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2374 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2376 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2377 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2378 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2379 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2380 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2381 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2382 ecosap=ecosa1+ecosa2
2383 ecosbp=ecosb1+ecosb2
2384 ecosgp=ecosg1+ecosg2
2385 ecosam=ecosa1-ecosa2
2386 ecosbm=ecosb1-ecosb2
2387 ecosgm=ecosg1-ecosg2
2396 facont_hb(num_conti,i)=fcont
2397 fprimcont=fprimcont/rij
2398 cd facont_hb(num_conti,i)=1.0D0
2399 C Following line is for diagnostics.
2402 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2403 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2406 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2407 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2409 gggp(1)=gggp(1)+ees0pijp*xj
2410 gggp(2)=gggp(2)+ees0pijp*yj
2411 gggp(3)=gggp(3)+ees0pijp*zj
2412 gggm(1)=gggm(1)+ees0mijp*xj
2413 gggm(2)=gggm(2)+ees0mijp*yj
2414 gggm(3)=gggm(3)+ees0mijp*zj
2415 C Derivatives due to the contact function
2416 gacont_hbr(1,num_conti,i)=fprimcont*xj
2417 gacont_hbr(2,num_conti,i)=fprimcont*yj
2418 gacont_hbr(3,num_conti,i)=fprimcont*zj
2421 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2422 c following the change of gradient-summation algorithm.
2424 cgrad ghalfp=0.5D0*gggp(k)
2425 cgrad ghalfm=0.5D0*gggm(k)
2426 gacontp_hb1(k,num_conti,i)=!ghalfp
2427 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2428 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2429 & *fac_shield(i)*fac_shield(j)
2431 gacontp_hb2(k,num_conti,i)=!ghalfp
2432 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2433 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2434 & *fac_shield(i)*fac_shield(j)
2436 gacontp_hb3(k,num_conti,i)=gggp(k)
2437 & *fac_shield(i)*fac_shield(j)
2439 gacontm_hb1(k,num_conti,i)=!ghalfm
2440 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2441 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2442 & *fac_shield(i)*fac_shield(j)
2444 gacontm_hb2(k,num_conti,i)=!ghalfm
2445 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2446 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2447 & *fac_shield(i)*fac_shield(j)
2449 gacontm_hb3(k,num_conti,i)=gggm(k)
2450 & *fac_shield(i)*fac_shield(j)
2454 endif ! num_conti.le.maxconts
2458 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2461 ghalf=0.5d0*agg(l,k)
2462 aggi(l,k)=aggi(l,k)+ghalf
2463 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2464 aggj(l,k)=aggj(l,k)+ghalf
2467 if (j.eq.nres-1 .and. i.lt.j-2) then
2470 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2475 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2478 C-----------------------------------------------------------------------
2479 subroutine evdwpp_short(evdw1)
2483 implicit real*8 (a-h,o-z)
2484 include 'DIMENSIONS'
2485 include 'COMMON.CONTROL'
2486 include 'COMMON.IOUNITS'
2487 include 'COMMON.GEO'
2488 include 'COMMON.VAR'
2489 include 'COMMON.LOCAL'
2490 include 'COMMON.CHAIN'
2491 include 'COMMON.DERIV'
2492 include 'COMMON.INTERACT'
2493 c include 'COMMON.CONTACTS'
2494 include 'COMMON.TORSION'
2495 include 'COMMON.VECTORS'
2496 include 'COMMON.FFIELD'
2498 integer xshift,yshift,zshift
2499 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2501 double precision scal_el /1.0d0/
2503 double precision scal_el /0.5d0/
2505 c write (iout,*) "evdwpp_short"
2508 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2509 c & " iatel_e_vdw",iatel_e_vdw
2511 do i=iatel_s_vdw,iatel_e_vdw
2512 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2516 dx_normi=dc_norm(1,i)
2517 dy_normi=dc_norm(2,i)
2518 dz_normi=dc_norm(3,i)
2519 xmedi=c(1,i)+0.5d0*dxi
2520 ymedi=c(2,i)+0.5d0*dyi
2521 zmedi=c(3,i)+0.5d0*dzi
2522 xmedi=mod(xmedi,boxxsize)
2523 if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
2524 ymedi=mod(ymedi,boxysize)
2525 if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
2526 zmedi=mod(zmedi,boxzsize)
2527 if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
2529 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2530 c & ' ielend',ielend_vdw(i)
2532 do j=ielstart_vdw(i),ielend_vdw(i)
2533 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2537 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2538 aaa=app(iteli,itelj)
2539 bbb=bpp(iteli,itelj)
2543 dx_normj=dc_norm(1,j)
2544 dy_normj=dc_norm(2,j)
2545 dz_normj=dc_norm(3,j)
2550 if (xj.lt.0) xj=xj+boxxsize
2552 if (yj.lt.0) yj=yj+boxysize
2554 if (zj.lt.0) zj=zj+boxzsize
2555 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2563 xj=xj_safe+xshift*boxxsize
2564 yj=yj_safe+yshift*boxysize
2565 zj=zj_safe+zshift*boxzsize
2566 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2567 if(dist_temp.lt.dist_init) then
2577 if (isubchap.eq.1) then
2586 rij=xj*xj+yj*yj+zj*zj
2589 c sss=sscale(rij/rpp(iteli,itelj))
2590 c sssgrad=sscagrad(rij/rpp(iteli,itelj))
2592 sssgrad=sscagrad(rij)
2593 if (sss.gt.0.0d0) then
2598 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2599 if (j.eq.i+2) ev1=scal_el*ev1
2602 if (energy_dec) then
2603 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2605 evdw1=evdw1+evdwij*sss
2606 if (energy_dec) write (iout,'(a10,2i5,0pf7.3)')
2607 & 'evdw1_sum',i,j,evdw1
2609 C Calculate contributions to the Cartesian gradient.
2611 facvdw=-6*rrmij*(ev1+evdwij)*sss
2612 ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2613 ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2614 ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2619 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2620 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2627 C-----------------------------------------------------------------------------
2628 subroutine escp_long(evdw2,evdw2_14)
2630 C This subroutine calculates the excluded-volume interaction energy between
2631 C peptide-group centers and side chains and its gradient in virtual-bond and
2632 C side-chain vectors.
2634 implicit real*8 (a-h,o-z)
2635 include 'DIMENSIONS'
2636 include 'COMMON.GEO'
2637 include 'COMMON.VAR'
2638 include 'COMMON.LOCAL'
2639 include 'COMMON.CHAIN'
2640 include 'COMMON.DERIV'
2641 include 'COMMON.INTERACT'
2642 include 'COMMON.FFIELD'
2643 include 'COMMON.IOUNITS'
2644 include 'COMMON.CONTROL'
2645 logical lprint_short
2646 common /shortcheck/ lprint_short
2648 integer xshift,yshift,zshift
2649 if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
2652 CD print '(a)','Enter ESCP KURWA'
2653 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2655 c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
2656 c & ' iatscp_e=',iatscp_e
2657 do i=iatscp_s,iatscp_e
2658 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2660 xi=0.5D0*(c(1,i)+c(1,i+1))
2661 yi=0.5D0*(c(2,i)+c(2,i+1))
2662 zi=0.5D0*(c(3,i)+c(3,i+1))
2664 if (xi.lt.0) xi=xi+boxxsize
2666 if (yi.lt.0) yi=yi+boxysize
2668 if (zi.lt.0) zi=zi+boxzsize
2669 do iint=1,nscp_gr(i)
2671 do j=iscpstart(i,iint),iscpend(i,iint)
2672 itypj=iabs(itype(j))
2673 if (itypj.eq.ntyp1) cycle
2674 C Uncomment following three lines for SC-p interactions
2678 C Uncomment following three lines for Ca-p interactions
2684 if (xj.lt.0) xj=xj+boxxsize
2686 if (yj.lt.0) yj=yj+boxysize
2688 if (zj.lt.0) zj=zj+boxzsize
2690 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2698 xj=xj_safe+xshift*boxxsize
2699 yj=yj_safe+yshift*boxysize
2700 zj=zj_safe+zshift*boxzsize
2701 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2702 if(dist_temp.lt.dist_init) then
2712 if (subchap.eq.1) then
2722 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2724 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2725 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2726 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2727 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2728 if (sss.lt.1.0d0) then
2730 e1=fac*fac*aad(itypj,iteli)
2731 e2=fac*bad(itypj,iteli)
2732 if (iabs(j-i) .le. 2) then
2735 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2738 evdw2=evdw2+evdwij*(1.0d0-sss)
2739 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2740 & 'evdw2',i,j,sss,evdwij
2742 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2745 fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2746 fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2750 C Uncomment following three lines for SC-p interactions
2752 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2754 C Uncomment following line for SC-p interactions
2755 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2757 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2758 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2767 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2768 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2769 gradx_scp(j,i)=expon*gradx_scp(j,i)
2772 C******************************************************************************
2776 C To save time the factor EXPON has been extracted from ALL components
2777 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2780 C******************************************************************************
2783 C-----------------------------------------------------------------------------
2784 subroutine escp_short(evdw2,evdw2_14)
2786 C This subroutine calculates the excluded-volume interaction energy between
2787 C peptide-group centers and side chains and its gradient in virtual-bond and
2788 C side-chain vectors.
2790 implicit real*8 (a-h,o-z)
2791 include 'DIMENSIONS'
2792 include 'COMMON.GEO'
2793 include 'COMMON.VAR'
2794 include 'COMMON.LOCAL'
2795 include 'COMMON.CHAIN'
2796 include 'COMMON.DERIV'
2797 include 'COMMON.INTERACT'
2798 include 'COMMON.FFIELD'
2799 include 'COMMON.IOUNITS'
2800 include 'COMMON.CONTROL'
2801 integer xshift,yshift,zshift
2802 logical lprint_short
2803 common /shortcheck/ lprint_short
2807 cd print '(a)','Enter ESCP'
2809 c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
2810 c & ' iatscp_e=',iatscp_e
2811 if (energy_dec) write (iout,*) "escp_short:",r_cut,rlamb
2812 do i=iatscp_s,iatscp_e
2813 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2815 xi=0.5D0*(c(1,i)+c(1,i+1))
2816 yi=0.5D0*(c(2,i)+c(2,i+1))
2817 zi=0.5D0*(c(3,i)+c(3,i+1))
2819 if (xi.lt.0) xi=xi+boxxsize
2821 if (yi.lt.0) yi=yi+boxysize
2823 if (zi.lt.0) zi=zi+boxzsize
2826 c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
2827 c & " nscp_gr",nscp_gr(i)
2828 do iint=1,nscp_gr(i)
2830 do j=iscpstart(i,iint),iscpend(i,iint)
2831 itypj=iabs(itype(j))
2833 c & write (iout,*) "j",j," itypj",itypj
2834 if (itypj.eq.ntyp1) cycle
2835 C Uncomment following three lines for SC-p interactions
2839 C Uncomment following three lines for Ca-p interactions
2845 if (xj.lt.0) xj=xj+boxxsize
2847 if (yj.lt.0) yj=yj+boxysize
2849 if (zj.lt.0) zj=zj+boxzsize
2851 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2852 c if (lprint_short) then
2853 c write (iout,*) i,j,xi,yi,zi,xj,yj,zj
2854 c write (iout,*) "dist_init",dsqrt(dist_init)
2863 xj=xj_safe+xshift*boxxsize
2864 yj=yj_safe+yshift*boxysize
2865 zj=zj_safe+zshift*boxzsize
2866 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2867 if(dist_temp.lt.dist_init) then
2877 c if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp)
2878 if (subchap.eq.1) then
2887 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2888 c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2889 c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2890 sss=sscale(1.0d0/(dsqrt(rrij)))
2891 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
2892 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2893 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2894 c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),
2895 c & " subchap",subchap," sss",sss
2896 if (sss.gt.0.0d0) then
2899 e1=fac*fac*aad(itypj,iteli)
2900 e2=fac*bad(itypj,iteli)
2901 if (iabs(j-i) .le. 2) then
2904 evdw2_14=evdw2_14+(e1+e2)*sss
2907 evdw2=evdw2+evdwij*sss
2908 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2909 & 'evdw2',i,j,sss,evdwij
2911 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2913 fac=-(evdwij+e1)*rrij*sss
2914 fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2918 C Uncomment following three lines for SC-p interactions
2920 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2922 C Uncomment following line for SC-p interactions
2923 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2925 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2926 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2935 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2936 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2937 gradx_scp(j,i)=expon*gradx_scp(j,i)
2940 C******************************************************************************
2944 C To save time the factor EXPON has been extracted from ALL components
2945 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2948 C******************************************************************************