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 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 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'
1371 include 'COMMON.CONTACTS'
1372 include 'COMMON.TORSION'
1373 include 'COMMON.VECTORS'
1374 include 'COMMON.FFIELD'
1375 include 'COMMON.TIME1'
1376 include 'COMMON.SHIELD'
1377 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1378 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1379 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1380 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1381 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1382 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1384 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1386 double precision scal_el /1.0d0/
1388 double precision scal_el /0.5d0/
1391 C 13-go grudnia roku pamietnego...
1392 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1393 & 0.0d0,1.0d0,0.0d0,
1394 & 0.0d0,0.0d0,1.0d0/
1395 cd write(iout,*) 'In EELEC'
1397 cd write(iout,*) 'Type',i
1398 cd write(iout,*) 'B1',B1(:,i)
1399 cd write(iout,*) 'B2',B2(:,i)
1400 cd write(iout,*) 'CC',CC(:,:,i)
1401 cd write(iout,*) 'DD',DD(:,:,i)
1402 cd write(iout,*) 'EE',EE(:,:,i)
1404 cd call check_vecgrad
1406 C print *,"WCHODZE3"
1407 if (icheckgrad.eq.1) then
1409 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1411 dc_norm(k,i)=dc(k,i)*fac
1413 c write (iout,*) 'i',i,' fac',fac
1416 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1417 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1418 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1419 c call vec_and_deriv
1425 time_mat=time_mat+MPI_Wtime()-time01
1429 cd write (iout,*) 'i=',i
1431 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1434 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1435 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1448 cd print '(a)','Enter EELEC'
1449 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1451 gel_loc_loc(i)=0.0d0
1456 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1458 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1460 do i=iturn3_start,iturn3_end
1461 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1462 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1463 C & .or. itype(i-1).eq.ntyp1
1464 C & .or. itype(i+4).eq.ntyp1
1469 dx_normi=dc_norm(1,i)
1470 dy_normi=dc_norm(2,i)
1471 dz_normi=dc_norm(3,i)
1472 xmedi=c(1,i)+0.5d0*dxi
1473 ymedi=c(2,i)+0.5d0*dyi
1474 zmedi=c(3,i)+0.5d0*dzi
1475 xmedi=mod(xmedi,boxxsize)
1476 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1477 ymedi=mod(ymedi,boxysize)
1478 if (ymedi.lt.0) ymedi=ymedi+boxysize
1479 zmedi=mod(zmedi,boxzsize)
1480 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1482 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1483 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1484 num_cont_hb(i)=num_conti
1486 do i=iturn4_start,iturn4_end
1487 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1488 & .or. itype(i+3).eq.ntyp1
1489 & .or. itype(i+4).eq.ntyp1
1490 C & .or. itype(i+5).eq.ntyp1
1491 C & .or. itype(i-1).eq.ntyp1
1496 dx_normi=dc_norm(1,i)
1497 dy_normi=dc_norm(2,i)
1498 dz_normi=dc_norm(3,i)
1499 xmedi=c(1,i)+0.5d0*dxi
1500 ymedi=c(2,i)+0.5d0*dyi
1501 zmedi=c(3,i)+0.5d0*dzi
1502 xmedi=mod(xmedi,boxxsize)
1503 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1504 ymedi=mod(ymedi,boxysize)
1505 if (ymedi.lt.0) ymedi=ymedi+boxysize
1506 zmedi=mod(zmedi,boxzsize)
1507 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1508 num_conti=num_cont_hb(i)
1509 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1510 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1511 & call eturn4(i,eello_turn4)
1512 num_cont_hb(i)=num_conti
1515 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1517 do i=iatel_s,iatel_e
1518 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1519 C & .or. itype(i+2).eq.ntyp1
1520 C & .or. itype(i-1).eq.ntyp1
1525 dx_normi=dc_norm(1,i)
1526 dy_normi=dc_norm(2,i)
1527 dz_normi=dc_norm(3,i)
1528 xmedi=c(1,i)+0.5d0*dxi
1529 ymedi=c(2,i)+0.5d0*dyi
1530 zmedi=c(3,i)+0.5d0*dzi
1531 xmedi=mod(xmedi,boxxsize)
1532 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1533 ymedi=mod(ymedi,boxysize)
1534 if (ymedi.lt.0) ymedi=ymedi+boxysize
1535 zmedi=mod(zmedi,boxzsize)
1536 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1537 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1538 num_conti=num_cont_hb(i)
1539 do j=ielstart(i),ielend(i)
1540 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1541 C & .or.itype(j+2).eq.ntyp1
1542 C & .or.itype(j-1).eq.ntyp1
1544 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1546 num_cont_hb(i)=num_conti
1548 c write (iout,*) "Number of loop steps in EELEC:",ind
1550 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1551 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1553 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1554 ccc eel_loc=eel_loc+eello_turn3
1555 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1558 C-------------------------------------------------------------------------------
1559 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1560 implicit real*8 (a-h,o-z)
1561 include 'DIMENSIONS'
1565 include 'COMMON.CONTROL'
1566 include 'COMMON.IOUNITS'
1567 include 'COMMON.GEO'
1568 include 'COMMON.VAR'
1569 include 'COMMON.LOCAL'
1570 include 'COMMON.CHAIN'
1571 include 'COMMON.DERIV'
1572 include 'COMMON.INTERACT'
1573 include 'COMMON.CONTACTS'
1574 include 'COMMON.TORSION'
1575 include 'COMMON.VECTORS'
1576 include 'COMMON.FFIELD'
1577 include 'COMMON.TIME1'
1578 include 'COMMON.SHIELD'
1579 integer xshift,yshift,zshift
1580 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1581 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1582 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1583 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
1584 & gmuij2(4),gmuji2(4)
1585 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1586 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1588 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1590 double precision scal_el /1.0d0/
1592 double precision scal_el /0.5d0/
1595 C 13-go grudnia roku pamietnego...
1596 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1597 & 0.0d0,1.0d0,0.0d0,
1598 & 0.0d0,0.0d0,1.0d0/
1599 c time00=MPI_Wtime()
1600 cd write (iout,*) "eelecij",i,j
1601 C print *,"WCHODZE2"
1605 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1606 aaa=app(iteli,itelj)
1607 bbb=bpp(iteli,itelj)
1608 ael6i=ael6(iteli,itelj)
1609 ael3i=ael3(iteli,itelj)
1613 dx_normj=dc_norm(1,j)
1614 dy_normj=dc_norm(2,j)
1615 dz_normj=dc_norm(3,j)
1620 if (xj.lt.0) xj=xj+boxxsize
1622 if (yj.lt.0) yj=yj+boxysize
1624 if (zj.lt.0) zj=zj+boxzsize
1625 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1633 xj=xj_safe+xshift*boxxsize
1634 yj=yj_safe+yshift*boxysize
1635 zj=zj_safe+zshift*boxzsize
1636 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1637 if(dist_temp.lt.dist_init) then
1647 if (isubchap.eq.1) then
1657 rij=xj*xj+yj*yj+zj*zj
1661 c For extracting the short-range part of Evdwpp
1662 sss=sscale(rij/rpp(iteli,itelj))
1663 sssgrad=sscagrad(rij/rpp(iteli,itelj))
1666 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1667 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1668 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1669 fac=cosa-3.0D0*cosb*cosg
1671 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1672 if (j.eq.i+2) ev1=scal_el*ev1
1677 if (shield_mode.eq.0) then
1681 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1683 el1=el1*fac_shield(i)**2*fac_shield(j)**2
1684 el2=el2*fac_shield(i)**2*fac_shield(j)**2
1686 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1687 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1689 evdw1=evdw1+evdwij*(1.0d0-sss)
1690 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1691 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1692 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1693 cd & xmedi,ymedi,zmedi,xj,yj,zj
1695 if (energy_dec) then
1696 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1697 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1701 C Calculate contributions to the Cartesian gradient.
1704 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1705 facel=-3*rrmij*(el1+eesij)
1711 * Radial derivatives. First process both termini of the fragment (i,j)
1716 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1717 & (shield_mode.gt.0)) then
1719 do ilist=1,ishield_list(i)
1720 iresshield=shield_list(ilist,i)
1722 rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
1724 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1726 & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
1727 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1728 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1729 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1730 C if (iresshield.gt.i) then
1731 C do ishi=i+1,iresshield-1
1732 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1733 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1737 C do ishi=iresshield,i
1738 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1739 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1745 do ilist=1,ishield_list(j)
1746 iresshield=shield_list(ilist,j)
1748 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1750 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1752 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
1753 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1758 gshieldc(k,i)=gshieldc(k,i)+
1759 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
1760 gshieldc(k,j)=gshieldc(k,j)+
1761 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
1762 gshieldc(k,i-1)=gshieldc(k,i-1)+
1763 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
1764 gshieldc(k,j-1)=gshieldc(k,j-1)+
1765 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
1771 c ghalf=0.5D0*ggg(k)
1772 c gelc(k,i)=gelc(k,i)+ghalf
1773 c gelc(k,j)=gelc(k,j)+ghalf
1775 c 9/28/08 AL Gradient compotents will be summed only at the end
1777 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1778 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1781 * Loop over residues i+1 thru j-1.
1785 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1788 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1789 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1790 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1792 c ghalf=0.5D0*ggg(k)
1793 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1794 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1796 c 9/28/08 AL Gradient compotents will be summed only at the end
1798 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1799 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1802 * Loop over residues i+1 thru j-1.
1806 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1810 facvdw=ev1+evdwij*(1.0d0-sss)
1813 fac=-3*rrmij*(facvdw+facvdw+facel)
1818 * Radial derivatives. First process both termini of the fragment (i,j)
1824 c ghalf=0.5D0*ggg(k)
1825 c gelc(k,i)=gelc(k,i)+ghalf
1826 c gelc(k,j)=gelc(k,j)+ghalf
1828 c 9/28/08 AL Gradient compotents will be summed only at the end
1830 gelc_long(k,j)=gelc(k,j)+ggg(k)
1831 gelc_long(k,i)=gelc(k,i)-ggg(k)
1834 * Loop over residues i+1 thru j-1.
1838 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1841 c 9/28/08 AL Gradient compotents will be summed only at the end
1845 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
1846 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
1847 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
1849 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1850 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1856 ecosa=2.0D0*fac3*fac1+fac4
1859 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1860 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1862 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1863 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1865 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1866 cd & (dcosg(k),k=1,3)
1868 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
1869 & fac_shield(i)**2*fac_shield(j)**2
1873 c ghalf=0.5D0*ggg(k)
1874 c gelc(k,i)=gelc(k,i)+ghalf
1875 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1876 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1877 c gelc(k,j)=gelc(k,j)+ghalf
1878 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1879 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1883 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1888 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1889 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
1890 & *fac_shield(i)**2*fac_shield(j)**2
1893 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1894 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
1895 & *fac_shield(i)**2*fac_shield(j)**2
1896 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1897 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1899 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1900 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1901 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1903 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1904 C energy of a peptide unit is assumed in the form of a second-order
1905 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1906 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1907 C are computed for EVERY pair of non-contiguous peptide groups.
1909 if (j.lt.nres-1) then
1920 muij(kkk)=mu(k,i)*mu(l,j)
1922 gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
1923 c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
1924 gmuij2(kkk)=gUb2(k,i)*mu(l,j)
1925 gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
1926 c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
1927 gmuji2(kkk)=mu(k,i)*gUb2(l,j)
1931 cd write (iout,*) 'EELEC: i',i,' j',j
1932 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1933 cd write(iout,*) 'muij',muij
1934 ury=scalar(uy(1,i),erij)
1935 urz=scalar(uz(1,i),erij)
1936 vry=scalar(uy(1,j),erij)
1937 vrz=scalar(uz(1,j),erij)
1938 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1939 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1940 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1941 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1942 fac=dsqrt(-ael6i)*r3ij
1947 cd write (iout,'(4i5,4f10.5)')
1948 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1949 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1950 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1951 cd & uy(:,j),uz(:,j)
1952 cd write (iout,'(4f10.5)')
1953 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1954 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1955 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
1956 cd write (iout,'(9f10.5/)')
1957 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1958 C Derivatives of the elements of A in virtual-bond vectors
1959 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1961 uryg(k,1)=scalar(erder(1,k),uy(1,i))
1962 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1963 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1964 urzg(k,1)=scalar(erder(1,k),uz(1,i))
1965 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1966 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1967 vryg(k,1)=scalar(erder(1,k),uy(1,j))
1968 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1969 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1970 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1971 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1972 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1974 C Compute radial contributions to the gradient
1992 C Add the contributions coming from er
1995 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1996 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1997 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1998 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
2001 C Derivatives in DC(i)
2002 cgrad ghalf1=0.5d0*agg(k,1)
2003 cgrad ghalf2=0.5d0*agg(k,2)
2004 cgrad ghalf3=0.5d0*agg(k,3)
2005 cgrad ghalf4=0.5d0*agg(k,4)
2006 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2007 & -3.0d0*uryg(k,2)*vry)!+ghalf1
2008 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2009 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
2010 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2011 & -3.0d0*urzg(k,2)*vry)!+ghalf3
2012 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2013 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
2014 C Derivatives in DC(i+1)
2015 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2016 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2017 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2018 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2019 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2020 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2021 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2022 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2023 C Derivatives in DC(j)
2024 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2025 & -3.0d0*vryg(k,2)*ury)!+ghalf1
2026 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2027 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
2028 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2029 & -3.0d0*vryg(k,2)*urz)!+ghalf3
2030 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2031 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
2032 C Derivatives in DC(j+1) or DC(nres-1)
2033 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2034 & -3.0d0*vryg(k,3)*ury)
2035 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2036 & -3.0d0*vrzg(k,3)*ury)
2037 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2038 & -3.0d0*vryg(k,3)*urz)
2039 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2040 & -3.0d0*vrzg(k,3)*urz)
2041 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
2043 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
2056 aggi(k,l)=-aggi(k,l)
2057 aggi1(k,l)=-aggi1(k,l)
2058 aggj(k,l)=-aggj(k,l)
2059 aggj1(k,l)=-aggj1(k,l)
2062 if (j.lt.nres-1) then
2068 aggi(k,l)=-aggi(k,l)
2069 aggi1(k,l)=-aggi1(k,l)
2070 aggj(k,l)=-aggj(k,l)
2071 aggj1(k,l)=-aggj1(k,l)
2082 aggi(k,l)=-aggi(k,l)
2083 aggi1(k,l)=-aggi1(k,l)
2084 aggj(k,l)=-aggj(k,l)
2085 aggj1(k,l)=-aggj1(k,l)
2090 IF (wel_loc.gt.0.0d0) THEN
2091 C Contribution to the local-electrostatic energy coming from the i-j pair
2092 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2094 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2096 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2097 & 'eelloc',i,j,eel_loc_ij
2100 if (shield_mode.eq.0) then
2107 eel_loc_ij=eel_loc_ij
2108 & *fac_shield(i)*fac_shield(j)
2109 eel_loc=eel_loc+eel_loc_ij
2111 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2112 & (shield_mode.gt.0)) then
2115 do ilist=1,ishield_list(i)
2116 iresshield=shield_list(ilist,i)
2118 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2121 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2123 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2124 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2128 do ilist=1,ishield_list(j)
2129 iresshield=shield_list(ilist,j)
2131 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2134 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2136 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2137 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2144 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2145 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2146 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2147 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2148 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2149 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2150 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2151 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2156 geel_loc_ij=(a22*gmuij1(1)
2160 & *fac_shield(i)*fac_shield(j)
2161 c write(iout,*) "derivative over thatai"
2162 c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
2164 gloc(nphi+i,icg)=gloc(nphi+i,icg)+
2165 & geel_loc_ij*wel_loc
2166 c write(iout,*) "derivative over thatai-1"
2167 c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
2174 gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
2175 & geel_loc_ij*wel_loc
2176 & *fac_shield(i)*fac_shield(j)
2178 c Derivative over j residue
2179 geel_loc_ji=a22*gmuji1(1)
2183 c write(iout,*) "derivative over thataj"
2184 c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
2187 gloc(nphi+j,icg)=gloc(nphi+j,icg)+
2188 & geel_loc_ji*wel_loc
2189 & *fac_shield(i)*fac_shield(j)
2196 c write(iout,*) "derivative over thataj-1"
2197 c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
2199 gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
2200 & geel_loc_ji*wel_loc
2201 & *fac_shield(i)*fac_shield(j)
2203 cC Partial derivatives in virtual-bond dihedral angles gamma
2205 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2206 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2207 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2208 & *fac_shield(i)*fac_shield(j)
2210 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2211 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2212 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2213 & *fac_shield(i)*fac_shield(j)
2215 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2217 ggg(l)=(agg(l,1)*muij(1)+
2218 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2219 & *fac_shield(i)*fac_shield(j)
2221 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2222 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2223 cgrad ghalf=0.5d0*ggg(l)
2224 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
2225 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
2229 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2232 C Remaining derivatives of eello
2234 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2235 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2236 & *fac_shield(i)*fac_shield(j)
2238 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2239 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2240 & *fac_shield(i)*fac_shield(j)
2242 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2243 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2244 & *fac_shield(i)*fac_shield(j)
2246 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2247 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2248 & *fac_shield(i)*fac_shield(j)
2252 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2253 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
2254 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2255 & .and. num_conti.le.maxconts) then
2256 c write (iout,*) i,j," entered corr"
2258 C Calculate the contact function. The ith column of the array JCONT will
2259 C contain the numbers of atoms that make contacts with the atom I (of numbers
2260 C greater than I). The arrays FACONT and GACONT will contain the values of
2261 C the contact function and its derivative.
2262 c r0ij=1.02D0*rpp(iteli,itelj)
2263 c r0ij=1.11D0*rpp(iteli,itelj)
2264 r0ij=2.20D0*rpp(iteli,itelj)
2265 c r0ij=1.55D0*rpp(iteli,itelj)
2266 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2267 if (fcont.gt.0.0D0) then
2268 num_conti=num_conti+1
2269 if (num_conti.gt.maxconts) then
2270 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2271 & ' will skip next contacts for this conf.'
2273 jcont_hb(num_conti,i)=j
2274 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
2275 cd & " jcont_hb",jcont_hb(num_conti,i)
2276 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2277 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2278 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2280 d_cont(num_conti,i)=rij
2281 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2282 C --- Electrostatic-interaction matrix ---
2283 a_chuj(1,1,num_conti,i)=a22
2284 a_chuj(1,2,num_conti,i)=a23
2285 a_chuj(2,1,num_conti,i)=a32
2286 a_chuj(2,2,num_conti,i)=a33
2287 C --- Gradient of rij
2289 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2296 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2297 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2298 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2299 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2300 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2305 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2306 C Calculate contact energies
2308 wij=cosa-3.0D0*cosb*cosg
2311 c fac3=dsqrt(-ael6i)/r0ij**3
2312 fac3=dsqrt(-ael6i)*r3ij
2313 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2314 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2315 if (ees0tmp.gt.0) then
2316 ees0pij=dsqrt(ees0tmp)
2320 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2321 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2322 if (ees0tmp.gt.0) then
2323 ees0mij=dsqrt(ees0tmp)
2328 if (shield_mode.eq.0) then
2332 ees0plist(num_conti,i)=j
2333 C fac_shield(i)=0.4d0
2334 C fac_shield(j)=0.6d0
2336 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2337 & *fac_shield(i)*fac_shield(j)
2338 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2339 & *fac_shield(i)*fac_shield(j)
2341 C Diagnostics. Comment out or remove after debugging!
2342 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2343 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2344 c ees0m(num_conti,i)=0.0D0
2346 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2347 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2348 C Angular derivatives of the contact function
2349 ees0pij1=fac3/ees0pij
2350 ees0mij1=fac3/ees0mij
2351 fac3p=-3.0D0*fac3*rrmij
2352 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2353 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2355 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2356 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2357 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2358 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2359 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2360 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2361 ecosap=ecosa1+ecosa2
2362 ecosbp=ecosb1+ecosb2
2363 ecosgp=ecosg1+ecosg2
2364 ecosam=ecosa1-ecosa2
2365 ecosbm=ecosb1-ecosb2
2366 ecosgm=ecosg1-ecosg2
2375 facont_hb(num_conti,i)=fcont
2376 fprimcont=fprimcont/rij
2377 cd facont_hb(num_conti,i)=1.0D0
2378 C Following line is for diagnostics.
2381 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2382 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2385 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2386 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2388 gggp(1)=gggp(1)+ees0pijp*xj
2389 gggp(2)=gggp(2)+ees0pijp*yj
2390 gggp(3)=gggp(3)+ees0pijp*zj
2391 gggm(1)=gggm(1)+ees0mijp*xj
2392 gggm(2)=gggm(2)+ees0mijp*yj
2393 gggm(3)=gggm(3)+ees0mijp*zj
2394 C Derivatives due to the contact function
2395 gacont_hbr(1,num_conti,i)=fprimcont*xj
2396 gacont_hbr(2,num_conti,i)=fprimcont*yj
2397 gacont_hbr(3,num_conti,i)=fprimcont*zj
2400 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2401 c following the change of gradient-summation algorithm.
2403 cgrad ghalfp=0.5D0*gggp(k)
2404 cgrad ghalfm=0.5D0*gggm(k)
2405 gacontp_hb1(k,num_conti,i)=!ghalfp
2406 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2407 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2408 & *fac_shield(i)*fac_shield(j)
2410 gacontp_hb2(k,num_conti,i)=!ghalfp
2411 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2412 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2413 & *fac_shield(i)*fac_shield(j)
2415 gacontp_hb3(k,num_conti,i)=gggp(k)
2416 & *fac_shield(i)*fac_shield(j)
2418 gacontm_hb1(k,num_conti,i)=!ghalfm
2419 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2420 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2421 & *fac_shield(i)*fac_shield(j)
2423 gacontm_hb2(k,num_conti,i)=!ghalfm
2424 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2425 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2426 & *fac_shield(i)*fac_shield(j)
2428 gacontm_hb3(k,num_conti,i)=gggm(k)
2429 & *fac_shield(i)*fac_shield(j)
2433 endif ! num_conti.le.maxconts
2436 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2439 ghalf=0.5d0*agg(l,k)
2440 aggi(l,k)=aggi(l,k)+ghalf
2441 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2442 aggj(l,k)=aggj(l,k)+ghalf
2445 if (j.eq.nres-1 .and. i.lt.j-2) then
2448 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2453 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2456 C-----------------------------------------------------------------------
2457 subroutine evdwpp_short(evdw1)
2461 implicit real*8 (a-h,o-z)
2462 include 'DIMENSIONS'
2463 include 'COMMON.CONTROL'
2464 include 'COMMON.IOUNITS'
2465 include 'COMMON.GEO'
2466 include 'COMMON.VAR'
2467 include 'COMMON.LOCAL'
2468 include 'COMMON.CHAIN'
2469 include 'COMMON.DERIV'
2470 include 'COMMON.INTERACT'
2471 include 'COMMON.CONTACTS'
2472 include 'COMMON.TORSION'
2473 include 'COMMON.VECTORS'
2474 include 'COMMON.FFIELD'
2476 integer xshift,yshift,zshift
2477 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2479 double precision scal_el /1.0d0/
2481 double precision scal_el /0.5d0/
2483 c write (iout,*) "evdwpp_short"
2486 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2487 c & " iatel_e_vdw",iatel_e_vdw
2489 do i=iatel_s_vdw,iatel_e_vdw
2490 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2494 dx_normi=dc_norm(1,i)
2495 dy_normi=dc_norm(2,i)
2496 dz_normi=dc_norm(3,i)
2497 xmedi=c(1,i)+0.5d0*dxi
2498 ymedi=c(2,i)+0.5d0*dyi
2499 zmedi=c(3,i)+0.5d0*dzi
2500 xmedi=mod(xmedi,boxxsize)
2501 if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
2502 ymedi=mod(ymedi,boxysize)
2503 if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
2504 zmedi=mod(zmedi,boxzsize)
2505 if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
2507 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2508 c & ' ielend',ielend_vdw(i)
2510 do j=ielstart_vdw(i),ielend_vdw(i)
2511 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2515 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2516 aaa=app(iteli,itelj)
2517 bbb=bpp(iteli,itelj)
2521 dx_normj=dc_norm(1,j)
2522 dy_normj=dc_norm(2,j)
2523 dz_normj=dc_norm(3,j)
2528 if (xj.lt.0) xj=xj+boxxsize
2530 if (yj.lt.0) yj=yj+boxysize
2532 if (zj.lt.0) zj=zj+boxzsize
2533 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2541 xj=xj_safe+xshift*boxxsize
2542 yj=yj_safe+yshift*boxysize
2543 zj=zj_safe+zshift*boxzsize
2544 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2545 if(dist_temp.lt.dist_init) then
2555 if (isubchap.eq.1) then
2564 rij=xj*xj+yj*yj+zj*zj
2567 c sss=sscale(rij/rpp(iteli,itelj))
2568 c sssgrad=sscagrad(rij/rpp(iteli,itelj))
2570 sssgrad=sscagrad(rij)
2571 if (sss.gt.0.0d0) then
2576 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2577 if (j.eq.i+2) ev1=scal_el*ev1
2580 if (energy_dec) then
2581 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2583 evdw1=evdw1+evdwij*sss
2584 if (energy_dec) write (iout,'(a10,2i5,0pf7.3)')
2585 & 'evdw1_sum',i,j,evdw1
2587 C Calculate contributions to the Cartesian gradient.
2589 facvdw=-6*rrmij*(ev1+evdwij)*sss
2590 ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2591 ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2592 ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2597 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2598 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2605 C-----------------------------------------------------------------------------
2606 subroutine escp_long(evdw2,evdw2_14)
2608 C This subroutine calculates the excluded-volume interaction energy between
2609 C peptide-group centers and side chains and its gradient in virtual-bond and
2610 C side-chain vectors.
2612 implicit real*8 (a-h,o-z)
2613 include 'DIMENSIONS'
2614 include 'COMMON.GEO'
2615 include 'COMMON.VAR'
2616 include 'COMMON.LOCAL'
2617 include 'COMMON.CHAIN'
2618 include 'COMMON.DERIV'
2619 include 'COMMON.INTERACT'
2620 include 'COMMON.FFIELD'
2621 include 'COMMON.IOUNITS'
2622 include 'COMMON.CONTROL'
2623 logical lprint_short
2624 common /shortcheck/ lprint_short
2626 integer xshift,yshift,zshift
2627 if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
2630 CD print '(a)','Enter ESCP KURWA'
2631 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2633 c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
2634 c & ' iatscp_e=',iatscp_e
2635 do i=iatscp_s,iatscp_e
2636 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2638 xi=0.5D0*(c(1,i)+c(1,i+1))
2639 yi=0.5D0*(c(2,i)+c(2,i+1))
2640 zi=0.5D0*(c(3,i)+c(3,i+1))
2642 if (xi.lt.0) xi=xi+boxxsize
2644 if (yi.lt.0) yi=yi+boxysize
2646 if (zi.lt.0) zi=zi+boxzsize
2647 do iint=1,nscp_gr(i)
2649 do j=iscpstart(i,iint),iscpend(i,iint)
2650 itypj=iabs(itype(j))
2651 if (itypj.eq.ntyp1) cycle
2652 C Uncomment following three lines for SC-p interactions
2656 C Uncomment following three lines for Ca-p interactions
2662 if (xj.lt.0) xj=xj+boxxsize
2664 if (yj.lt.0) yj=yj+boxysize
2666 if (zj.lt.0) zj=zj+boxzsize
2668 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2676 xj=xj_safe+xshift*boxxsize
2677 yj=yj_safe+yshift*boxysize
2678 zj=zj_safe+zshift*boxzsize
2679 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2680 if(dist_temp.lt.dist_init) then
2690 if (subchap.eq.1) then
2700 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2702 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2703 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2704 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2705 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2706 if (sss.lt.1.0d0) then
2708 e1=fac*fac*aad(itypj,iteli)
2709 e2=fac*bad(itypj,iteli)
2710 if (iabs(j-i) .le. 2) then
2713 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2716 evdw2=evdw2+evdwij*(1.0d0-sss)
2717 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2718 & 'evdw2',i,j,sss,evdwij
2720 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2723 fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2724 fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2728 C Uncomment following three lines for SC-p interactions
2730 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2732 C Uncomment following line for SC-p interactions
2733 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2735 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2736 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2745 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2746 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2747 gradx_scp(j,i)=expon*gradx_scp(j,i)
2750 C******************************************************************************
2754 C To save time the factor EXPON has been extracted from ALL components
2755 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2758 C******************************************************************************
2761 C-----------------------------------------------------------------------------
2762 subroutine escp_short(evdw2,evdw2_14)
2764 C This subroutine calculates the excluded-volume interaction energy between
2765 C peptide-group centers and side chains and its gradient in virtual-bond and
2766 C side-chain vectors.
2768 implicit real*8 (a-h,o-z)
2769 include 'DIMENSIONS'
2770 include 'COMMON.GEO'
2771 include 'COMMON.VAR'
2772 include 'COMMON.LOCAL'
2773 include 'COMMON.CHAIN'
2774 include 'COMMON.DERIV'
2775 include 'COMMON.INTERACT'
2776 include 'COMMON.FFIELD'
2777 include 'COMMON.IOUNITS'
2778 include 'COMMON.CONTROL'
2779 integer xshift,yshift,zshift
2780 logical lprint_short
2781 common /shortcheck/ lprint_short
2785 cd print '(a)','Enter ESCP'
2787 c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
2788 c & ' iatscp_e=',iatscp_e
2789 if (energy_dec) write (iout,*) "escp_short:",r_cut,rlamb
2790 do i=iatscp_s,iatscp_e
2791 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2793 xi=0.5D0*(c(1,i)+c(1,i+1))
2794 yi=0.5D0*(c(2,i)+c(2,i+1))
2795 zi=0.5D0*(c(3,i)+c(3,i+1))
2797 if (xi.lt.0) xi=xi+boxxsize
2799 if (yi.lt.0) yi=yi+boxysize
2801 if (zi.lt.0) zi=zi+boxzsize
2804 c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
2805 c & " nscp_gr",nscp_gr(i)
2806 do iint=1,nscp_gr(i)
2808 do j=iscpstart(i,iint),iscpend(i,iint)
2809 itypj=iabs(itype(j))
2811 c & write (iout,*) "j",j," itypj",itypj
2812 if (itypj.eq.ntyp1) cycle
2813 C Uncomment following three lines for SC-p interactions
2817 C Uncomment following three lines for Ca-p interactions
2823 if (xj.lt.0) xj=xj+boxxsize
2825 if (yj.lt.0) yj=yj+boxysize
2827 if (zj.lt.0) zj=zj+boxzsize
2829 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2830 c if (lprint_short) then
2831 c write (iout,*) i,j,xi,yi,zi,xj,yj,zj
2832 c write (iout,*) "dist_init",dsqrt(dist_init)
2841 xj=xj_safe+xshift*boxxsize
2842 yj=yj_safe+yshift*boxysize
2843 zj=zj_safe+zshift*boxzsize
2844 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2845 if(dist_temp.lt.dist_init) then
2855 c if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp)
2856 if (subchap.eq.1) then
2865 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2866 c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2867 c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
2868 sss=sscale(1.0d0/(dsqrt(rrij)))
2869 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
2870 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2871 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2872 c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),
2873 c & " subchap",subchap," sss",sss
2874 if (sss.gt.0.0d0) then
2877 e1=fac*fac*aad(itypj,iteli)
2878 e2=fac*bad(itypj,iteli)
2879 if (iabs(j-i) .le. 2) then
2882 evdw2_14=evdw2_14+(e1+e2)*sss
2885 evdw2=evdw2+evdwij*sss
2886 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2887 & 'evdw2',i,j,sss,evdwij
2889 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2891 fac=-(evdwij+e1)*rrij*sss
2892 fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
2896 C Uncomment following three lines for SC-p interactions
2898 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2900 C Uncomment following line for SC-p interactions
2901 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2903 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2904 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2913 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2914 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2915 gradx_scp(j,i)=expon*gradx_scp(j,i)
2918 C******************************************************************************
2922 C To save time the factor EXPON has been extracted from ALL components
2923 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2926 C******************************************************************************