1 C-----------------------------------------------------------------------
2 double precision function sscalelip(r)
3 double precision r,gamm
4 include "COMMON.SPLITELE"
5 C if(r.lt.r_cut-rlamb) then
7 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
8 C gamm=(r-(r_cut-rlamb))/rlamb
9 sscalelip=1.0d0+r*r*(2.0d0*r-3.0d0)
15 C-----------------------------------------------------------------------
16 double precision function sscagradlip(r)
17 double precision r,gamm
18 include "COMMON.SPLITELE"
19 C if(r.lt.r_cut-rlamb) then
21 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
22 C gamm=(r-(r_cut-rlamb))/rlamb
23 sscagradlip=r*(6.0d0*r-6.0d0)
30 C-----------------------------------------------------------------------
31 double precision function sscale(r)
32 double precision r,gamm
33 include "COMMON.SPLITELE"
34 if(r.lt.r_cut-rlamb) then
36 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
37 gamm=(r-(r_cut-rlamb))/rlamb
38 sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
44 C-----------------------------------------------------------------------
45 C-----------------------------------------------------------------------
46 double precision function sscagrad(r)
47 double precision r,gamm
48 include "COMMON.SPLITELE"
49 if(r.lt.r_cut-rlamb) then
51 else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
52 gamm=(r-(r_cut-rlamb))/rlamb
53 sscagrad=gamm*(6*gamm-6.0d0)/rlamb
59 C-----------------------------------------------------------------------
60 subroutine elj_long(evdw)
62 C This subroutine calculates the interaction energy of nonbonded side chains
63 C assuming the LJ potential of interaction.
65 implicit real*8 (a-h,o-z)
67 parameter (accur=1.0d-10)
70 include 'COMMON.LOCAL'
71 include 'COMMON.CHAIN'
72 include 'COMMON.DERIV'
73 include 'COMMON.INTERACT'
74 include 'COMMON.TORSION'
75 include 'COMMON.SBRIDGE'
76 include 'COMMON.NAMES'
77 include 'COMMON.IOUNITS'
78 include 'COMMON.CONTACTS'
80 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
84 if (itypi.eq.ntyp1) cycle
90 C Calculate SC interaction energy.
93 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
94 cd & 'iend=',iend(i,iint)
95 do j=istart(i,iint),iend(i,iint)
97 if (itypj.eq.ntyp1) cycle
101 rij=xj*xj+yj*yj+zj*zj
102 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
103 if (sss.lt.1.0d0) then
105 eps0ij=eps(itypi,itypj)
110 evdw=evdw+(1.0d0-sss)*evdwij
112 C Calculate the components of the gradient in DC and X
114 fac=-rrij*(e1+evdwij)*(1.0d0-sss)
119 gvdwx(k,i)=gvdwx(k,i)-gg(k)
120 gvdwx(k,j)=gvdwx(k,j)+gg(k)
121 gvdwc(k,i)=gvdwc(k,i)-gg(k)
122 gvdwc(k,j)=gvdwc(k,j)+gg(k)
130 gvdwc(j,i)=expon*gvdwc(j,i)
131 gvdwx(j,i)=expon*gvdwx(j,i)
134 C******************************************************************************
138 C To save time, the factor of EXPON has been extracted from ALL components
139 C of GVDWC and GRADX. Remember to multiply them by this factor before further
142 C******************************************************************************
145 C-----------------------------------------------------------------------
146 subroutine elj_short(evdw)
148 C This subroutine calculates the interaction energy of nonbonded side chains
149 C assuming the LJ potential of interaction.
151 implicit real*8 (a-h,o-z)
153 parameter (accur=1.0d-10)
156 include 'COMMON.LOCAL'
157 include 'COMMON.CHAIN'
158 include 'COMMON.DERIV'
159 include 'COMMON.INTERACT'
160 include 'COMMON.TORSION'
161 include 'COMMON.SBRIDGE'
162 include 'COMMON.NAMES'
163 include 'COMMON.IOUNITS'
164 include 'COMMON.CONTACTS'
166 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
170 if (itypi.eq.ntyp1) cycle
178 C Calculate SC interaction energy.
181 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
182 cd & 'iend=',iend(i,iint)
183 do j=istart(i,iint),iend(i,iint)
185 if (itypj.eq.ntyp1) cycle
189 C Change 12/1/95 to calculate four-body interactions
190 rij=xj*xj+yj*yj+zj*zj
191 sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
192 if (sss.gt.0.0d0) then
194 eps0ij=eps(itypi,itypj)
201 C Calculate the components of the gradient in DC and X
203 fac=-rrij*(e1+evdwij)*sss
208 gvdwx(k,i)=gvdwx(k,i)-gg(k)
209 gvdwx(k,j)=gvdwx(k,j)+gg(k)
210 gvdwc(k,i)=gvdwc(k,i)-gg(k)
211 gvdwc(k,j)=gvdwc(k,j)+gg(k)
219 gvdwc(j,i)=expon*gvdwc(j,i)
220 gvdwx(j,i)=expon*gvdwx(j,i)
223 C******************************************************************************
227 C To save time, the factor of EXPON has been extracted from ALL components
228 C of GVDWC and GRADX. Remember to multiply them by this factor before further
231 C******************************************************************************
234 C-----------------------------------------------------------------------------
235 subroutine eljk_long(evdw)
237 C This subroutine calculates the interaction energy of nonbonded side chains
238 C assuming the LJK potential of interaction.
240 implicit real*8 (a-h,o-z)
244 include 'COMMON.LOCAL'
245 include 'COMMON.CHAIN'
246 include 'COMMON.DERIV'
247 include 'COMMON.INTERACT'
248 include 'COMMON.IOUNITS'
249 include 'COMMON.NAMES'
252 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
256 if (itypi.eq.ntyp1) cycle
262 C Calculate SC interaction energy.
265 do j=istart(i,iint),iend(i,iint)
267 if (itypj.eq.ntyp1) cycle
271 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
273 e_augm=augm(itypi,itypj)*fac_augm
276 sss=sscale(rij/sigma(itypi,itypj))
277 if (sss.lt.1.0d0) then
278 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
279 fac=r_shift_inv**expon
283 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
284 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
285 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
286 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
287 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
288 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
289 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
290 evdw=evdw+(1.0d0-sss)*evdwij
292 C Calculate the components of the gradient in DC and X
294 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
300 gvdwx(k,i)=gvdwx(k,i)-gg(k)
301 gvdwx(k,j)=gvdwx(k,j)+gg(k)
302 gvdwc(k,i)=gvdwc(k,i)-gg(k)
303 gvdwc(k,j)=gvdwc(k,j)+gg(k)
311 gvdwc(j,i)=expon*gvdwc(j,i)
312 gvdwx(j,i)=expon*gvdwx(j,i)
317 C-----------------------------------------------------------------------------
318 subroutine eljk_short(evdw)
320 C This subroutine calculates the interaction energy of nonbonded side chains
321 C assuming the LJK potential of interaction.
323 implicit real*8 (a-h,o-z)
327 include 'COMMON.LOCAL'
328 include 'COMMON.CHAIN'
329 include 'COMMON.DERIV'
330 include 'COMMON.INTERACT'
331 include 'COMMON.IOUNITS'
332 include 'COMMON.NAMES'
335 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
339 if (itypi.eq.ntyp1) cycle
345 C Calculate SC interaction energy.
348 do j=istart(i,iint),iend(i,iint)
350 if (itypj.eq.ntyp1) cycle
354 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
356 e_augm=augm(itypi,itypj)*fac_augm
359 sss=sscale(rij/sigma(itypi,itypj))
360 if (sss.gt.0.0d0) then
361 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
362 fac=r_shift_inv**expon
366 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
367 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
368 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
369 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
370 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
371 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
372 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
375 C Calculate the components of the gradient in DC and X
377 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
383 gvdwx(k,i)=gvdwx(k,i)-gg(k)
384 gvdwx(k,j)=gvdwx(k,j)+gg(k)
385 gvdwc(k,i)=gvdwc(k,i)-gg(k)
386 gvdwc(k,j)=gvdwc(k,j)+gg(k)
394 gvdwc(j,i)=expon*gvdwc(j,i)
395 gvdwx(j,i)=expon*gvdwx(j,i)
400 C-----------------------------------------------------------------------------
401 subroutine ebp_long(evdw)
403 C This subroutine calculates the interaction energy of nonbonded side chains
404 C assuming the Berne-Pechukas potential of interaction.
406 implicit real*8 (a-h,o-z)
410 include 'COMMON.LOCAL'
411 include 'COMMON.CHAIN'
412 include 'COMMON.DERIV'
413 include 'COMMON.NAMES'
414 include 'COMMON.INTERACT'
415 include 'COMMON.IOUNITS'
416 include 'COMMON.CALC'
418 c double precision rrsave(maxdim)
421 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
423 c if (icall.eq.0) then
431 if (itypi.eq.ntyp1) cycle
436 dxi=dc_norm(1,nres+i)
437 dyi=dc_norm(2,nres+i)
438 dzi=dc_norm(3,nres+i)
439 c dsci_inv=dsc_inv(itypi)
440 dsci_inv=vbld_inv(i+nres)
442 C Calculate SC interaction energy.
445 do j=istart(i,iint),iend(i,iint)
448 if (itypj.eq.ntyp1) cycle
449 c dscj_inv=dsc_inv(itypj)
450 dscj_inv=vbld_inv(j+nres)
451 chi1=chi(itypi,itypj)
452 chi2=chi(itypj,itypi)
459 alf12=0.5D0*(alf1+alf2)
463 dxj=dc_norm(1,nres+j)
464 dyj=dc_norm(2,nres+j)
465 dzj=dc_norm(3,nres+j)
466 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
468 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
470 if (sss.lt.1.0d0) then
472 C Calculate the angle-dependent terms of energy & contributions to derivatives.
474 C Calculate whole angle-dependent part of epsilon and contributions
476 fac=(rrij*sigsq)**expon2
479 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
480 eps2der=evdwij*eps3rt
481 eps3der=evdwij*eps2rt
482 evdwij=evdwij*eps2rt*eps3rt
483 evdw=evdw+evdwij*(1.0d0-sss)
485 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
487 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
488 cd & restyp(itypi),i,restyp(itypj),j,
489 cd & epsi,sigm,chi1,chi2,chip1,chip2,
490 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
491 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
494 C Calculate gradient components.
495 e1=e1*eps1*eps2rt**2*eps3rt**2
496 fac=-expon*(e1+evdwij)
499 C Calculate radial part of the gradient
503 C Calculate the angular part of the gradient and sum add the contributions
504 C to the appropriate components of the Cartesian gradient.
505 call sc_grad_scale(1.0d0-sss)
513 C-----------------------------------------------------------------------------
514 subroutine ebp_short(evdw)
516 C This subroutine calculates the interaction energy of nonbonded side chains
517 C assuming the Berne-Pechukas potential of interaction.
519 implicit real*8 (a-h,o-z)
523 include 'COMMON.LOCAL'
524 include 'COMMON.CHAIN'
525 include 'COMMON.DERIV'
526 include 'COMMON.NAMES'
527 include 'COMMON.INTERACT'
528 include 'COMMON.IOUNITS'
529 include 'COMMON.CALC'
531 c double precision rrsave(maxdim)
534 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
536 c if (icall.eq.0) then
544 if (itypi.eq.ntyp1) cycle
549 dxi=dc_norm(1,nres+i)
550 dyi=dc_norm(2,nres+i)
551 dzi=dc_norm(3,nres+i)
552 c dsci_inv=dsc_inv(itypi)
553 dsci_inv=vbld_inv(i+nres)
555 C Calculate SC interaction energy.
558 do j=istart(i,iint),iend(i,iint)
561 if (itypj.eq.ntyp1) cycle
562 c dscj_inv=dsc_inv(itypj)
563 dscj_inv=vbld_inv(j+nres)
564 chi1=chi(itypi,itypj)
565 chi2=chi(itypj,itypi)
572 alf12=0.5D0*(alf1+alf2)
576 dxj=dc_norm(1,nres+j)
577 dyj=dc_norm(2,nres+j)
578 dzj=dc_norm(3,nres+j)
579 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
581 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
583 if (sss.gt.0.0d0) then
585 C Calculate the angle-dependent terms of energy & contributions to derivatives.
587 C Calculate whole angle-dependent part of epsilon and contributions
589 fac=(rrij*sigsq)**expon2
592 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
593 eps2der=evdwij*eps3rt
594 eps3der=evdwij*eps2rt
595 evdwij=evdwij*eps2rt*eps3rt
598 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
600 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
601 cd & restyp(itypi),i,restyp(itypj),j,
602 cd & epsi,sigm,chi1,chi2,chip1,chip2,
603 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
604 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
607 C Calculate gradient components.
608 e1=e1*eps1*eps2rt**2*eps3rt**2
609 fac=-expon*(e1+evdwij)
612 C Calculate radial part of the gradient
616 C Calculate the angular part of the gradient and sum add the contributions
617 C to the appropriate components of the Cartesian gradient.
618 call sc_grad_scale(sss)
626 C-----------------------------------------------------------------------------
627 subroutine egb_long(evdw)
629 C This subroutine calculates the interaction energy of nonbonded side chains
630 C assuming the Gay-Berne potential of interaction.
632 implicit real*8 (a-h,o-z)
636 include 'COMMON.LOCAL'
637 include 'COMMON.CHAIN'
638 include 'COMMON.DERIV'
639 include 'COMMON.NAMES'
640 include 'COMMON.INTERACT'
641 include 'COMMON.IOUNITS'
642 include 'COMMON.CALC'
643 include 'COMMON.CONTROL'
646 ccccc energy_dec=.false.
647 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
650 c if (icall.eq.0) lprn=.false.
654 if (itypi.eq.ntyp1) cycle
660 if (xi.lt.0) xi=xi+boxxsize
662 if (yi.lt.0) yi=yi+boxysize
664 if (zi.lt.0) zi=zi+boxzsize
665 if ((zi.gt.bordlipbot)
666 &.and.(zi.lt.bordliptop)) then
667 C the energy transfer exist
668 if (zi.lt.buflipbot) then
669 C what fraction I am in
671 & ((zi-bordlipbot)/lipbufthick)
672 C lipbufthick is thickenes of lipid buffore
673 sslipi=sscalelip(fracinbuf)
674 ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
675 elseif (zi.gt.bufliptop) then
676 fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
677 sslipi=sscalelip(fracinbuf)
678 ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
688 dxi=dc_norm(1,nres+i)
689 dyi=dc_norm(2,nres+i)
690 dzi=dc_norm(3,nres+i)
691 c dsci_inv=dsc_inv(itypi)
692 dsci_inv=vbld_inv(i+nres)
693 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
694 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
696 C Calculate SC interaction energy.
699 do j=istart(i,iint),iend(i,iint)
702 if (itypj.eq.ntyp1) cycle
703 c dscj_inv=dsc_inv(itypj)
704 dscj_inv=vbld_inv(j+nres)
705 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
706 c & 1.0d0/vbld(j+nres)
707 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
708 sig0ij=sigma(itypi,itypj)
709 chi1=chi(itypi,itypj)
710 chi2=chi(itypj,itypi)
717 alf12=0.5D0*(alf1+alf2)
722 if (xj.lt.0) xj=xj+boxxsize
724 if (yj.lt.0) yj=yj+boxysize
726 if (zj.lt.0) zj=zj+boxzsize
727 if ((zj.gt.bordlipbot)
728 &.and.(zj.lt.bordliptop)) then
729 C the energy transfer exist
730 if (zj.lt.buflipbot) then
731 C what fraction I am in
733 & ((zj-bordlipbot)/lipbufthick)
734 C lipbufthick is thickenes of lipid buffore
735 sslipj=sscalelip(fracinbuf)
736 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
737 elseif (zj.gt.bufliptop) then
738 fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
739 sslipj=sscalelip(fracinbuf)
740 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
750 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
751 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
752 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
753 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
755 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
763 xj=xj_safe+xshift*boxxsize
764 yj=yj_safe+yshift*boxysize
765 zj=zj_safe+zshift*boxzsize
766 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
767 if(dist_temp.lt.dist_init) then
777 if (subchap.eq.1) then
786 dxj=dc_norm(1,nres+j)
787 dyj=dc_norm(2,nres+j)
788 dzj=dc_norm(3,nres+j)
789 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
791 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
792 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
793 if (sss.lt.1.0d0) then
795 C Calculate angle-dependent terms of energy and contributions to their
799 sig=sig0ij*dsqrt(sigsq)
800 rij_shift=1.0D0/rij-sig+sig0ij
801 c for diagnostics; uncomment
802 c rij_shift=1.2*sig0ij
803 C I hate to put IF's in the loops, but here don't have another choice!!!!
804 if (rij_shift.le.0.0D0) then
806 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
807 cd & restyp(itypi),i,restyp(itypj),j,
808 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
812 c---------------------------------------------------------------
813 rij_shift=1.0D0/rij_shift
818 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
819 eps2der=evdwij*eps3rt
820 eps3der=evdwij*eps2rt
821 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
822 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
823 evdwij=evdwij*eps2rt*eps3rt
824 evdw=evdw+evdwij*(1.0d0-sss)
826 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
828 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
829 & restyp(itypi),i,restyp(itypj),j,
830 & epsi,sigm,chi1,chi2,chip1,chip2,
831 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
832 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
836 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
839 C Calculate gradient components.
840 e1=e1*eps1*eps2rt**2*eps3rt**2
841 fac=-expon*(e1+evdwij)*rij_shift
844 fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij
846 C Calculate the radial part of the gradient
850 gg_lipi(3)=eps1*(eps2rt*eps2rt)
851 &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
852 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
853 &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
854 gg_lipj(3)=ssgradlipj*gg_lipi(3)
855 gg_lipi(3)=gg_lipi(3)*ssgradlipi
856 C Calculate angular part of the gradient.
857 call sc_grad_scale(1.0d0-sss)
862 c write (iout,*) "Number of loop steps in EGB:",ind
863 cccc energy_dec=.false.
866 C-----------------------------------------------------------------------------
867 subroutine egb_short(evdw)
869 C This subroutine calculates the interaction energy of nonbonded side chains
870 C assuming the Gay-Berne potential of interaction.
872 implicit real*8 (a-h,o-z)
876 include 'COMMON.LOCAL'
877 include 'COMMON.CHAIN'
878 include 'COMMON.DERIV'
879 include 'COMMON.NAMES'
880 include 'COMMON.INTERACT'
881 include 'COMMON.IOUNITS'
882 include 'COMMON.CALC'
883 include 'COMMON.CONTROL'
885 integer xshift,yshift,zshift
887 ccccc energy_dec=.false.
888 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
891 c if (icall.eq.0) lprn=.false.
895 if (itypi.eq.ntyp1) cycle
901 if (xi.lt.0) xi=xi+boxxsize
903 if (yi.lt.0) yi=yi+boxysize
905 if (zi.lt.0) zi=zi+boxzsize
906 if ((zi.gt.bordlipbot)
907 &.and.(zi.lt.bordliptop)) then
908 C the energy transfer exist
909 if (zi.lt.buflipbot) then
910 C what fraction I am in
912 & ((zi-bordlipbot)/lipbufthick)
913 C lipbufthick is thickenes of lipid buffore
914 sslipi=sscalelip(fracinbuf)
915 ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
916 elseif (zi.gt.bufliptop) then
917 fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
918 sslipi=sscalelip(fracinbuf)
919 ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
929 dxi=dc_norm(1,nres+i)
930 dyi=dc_norm(2,nres+i)
931 dzi=dc_norm(3,nres+i)
932 c dsci_inv=dsc_inv(itypi)
933 dsci_inv=vbld_inv(i+nres)
934 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
935 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
937 C Calculate SC interaction energy.
940 do j=istart(i,iint),iend(i,iint)
943 if (itypj.eq.ntyp1) cycle
944 c dscj_inv=dsc_inv(itypj)
945 dscj_inv=vbld_inv(j+nres)
946 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
947 c & 1.0d0/vbld(j+nres)
948 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
949 sig0ij=sigma(itypi,itypj)
950 chi1=chi(itypi,itypj)
951 chi2=chi(itypj,itypi)
958 alf12=0.5D0*(alf1+alf2)
963 if (xj.lt.0) xj=xj+boxxsize
965 if (yj.lt.0) yj=yj+boxysize
967 if (zj.lt.0) zj=zj+boxzsize
968 if ((zj.gt.bordlipbot)
969 &.and.(zj.lt.bordliptop)) then
970 C the energy transfer exist
971 if (zj.lt.buflipbot) then
972 C what fraction I am in
974 & ((zj-bordlipbot)/lipbufthick)
975 C lipbufthick is thickenes of lipid buffore
976 sslipj=sscalelip(fracinbuf)
977 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
978 elseif (zj.gt.bufliptop) then
979 fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
980 sslipj=sscalelip(fracinbuf)
981 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
990 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
991 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
992 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
993 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
994 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1002 xj=xj_safe+xshift*boxxsize
1003 yj=yj_safe+yshift*boxysize
1004 zj=zj_safe+zshift*boxzsize
1005 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1006 if(dist_temp.lt.dist_init) then
1016 if (subchap.eq.1) then
1025 dxj=dc_norm(1,nres+j)
1026 dyj=dc_norm(2,nres+j)
1027 dzj=dc_norm(3,nres+j)
1028 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1030 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1031 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
1032 if (sss.gt.0.0d0) then
1034 C Calculate angle-dependent terms of energy and contributions to their
1038 sig=sig0ij*dsqrt(sigsq)
1039 rij_shift=1.0D0/rij-sig+sig0ij
1040 c for diagnostics; uncomment
1041 c rij_shift=1.2*sig0ij
1042 C I hate to put IF's in the loops, but here don't have another choice!!!!
1043 if (rij_shift.le.0.0D0) then
1045 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1046 cd & restyp(itypi),i,restyp(itypj),j,
1047 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1051 c---------------------------------------------------------------
1052 rij_shift=1.0D0/rij_shift
1053 fac=rij_shift**expon
1057 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1058 eps2der=evdwij*eps3rt
1059 eps3der=evdwij*eps2rt
1060 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1061 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1062 evdwij=evdwij*eps2rt*eps3rt
1063 evdw=evdw+evdwij*sss
1065 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1067 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1068 & restyp(itypi),i,restyp(itypj),j,
1069 & epsi,sigm,chi1,chi2,chip1,chip2,
1070 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1071 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1075 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1078 C Calculate gradient components.
1079 e1=e1*eps1*eps2rt**2*eps3rt**2
1080 fac=-expon*(e1+evdwij)*rij_shift
1083 fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij
1085 C Calculate the radial part of the gradient
1089 gg_lipi(3)=eps1*(eps2rt*eps2rt)
1090 &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
1091 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
1092 &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
1093 gg_lipj(3)=ssgradlipj*gg_lipi(3)
1094 gg_lipi(3)=gg_lipi(3)*ssgradlipi
1095 C Calculate angular part of the gradient.
1096 call sc_grad_scale(sss)
1101 C write (iout,*) "Number of loop steps in EGB:",ind
1102 cccc energy_dec=.false.
1105 C-----------------------------------------------------------------------------
1106 subroutine egbv_long(evdw)
1108 C This subroutine calculates the interaction energy of nonbonded side chains
1109 C assuming the Gay-Berne-Vorobjev potential of interaction.
1111 implicit real*8 (a-h,o-z)
1112 include 'DIMENSIONS'
1113 include 'COMMON.GEO'
1114 include 'COMMON.VAR'
1115 include 'COMMON.LOCAL'
1116 include 'COMMON.CHAIN'
1117 include 'COMMON.DERIV'
1118 include 'COMMON.NAMES'
1119 include 'COMMON.INTERACT'
1120 include 'COMMON.IOUNITS'
1121 include 'COMMON.CALC'
1122 common /srutu/ icall
1125 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1128 c if (icall.eq.0) lprn=.true.
1130 do i=iatsc_s,iatsc_e
1132 if (itypi.eq.ntyp1) cycle
1137 dxi=dc_norm(1,nres+i)
1138 dyi=dc_norm(2,nres+i)
1139 dzi=dc_norm(3,nres+i)
1140 c dsci_inv=dsc_inv(itypi)
1141 dsci_inv=vbld_inv(i+nres)
1143 C Calculate SC interaction energy.
1145 do iint=1,nint_gr(i)
1146 do j=istart(i,iint),iend(i,iint)
1149 if (itypj.eq.ntyp1) cycle
1150 c dscj_inv=dsc_inv(itypj)
1151 dscj_inv=vbld_inv(j+nres)
1152 sig0ij=sigma(itypi,itypj)
1153 r0ij=r0(itypi,itypj)
1154 chi1=chi(itypi,itypj)
1155 chi2=chi(itypj,itypi)
1162 alf12=0.5D0*(alf1+alf2)
1166 dxj=dc_norm(1,nres+j)
1167 dyj=dc_norm(2,nres+j)
1168 dzj=dc_norm(3,nres+j)
1169 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1172 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1174 if (sss.lt.1.0d0) then
1176 C Calculate angle-dependent terms of energy and contributions to their
1180 sig=sig0ij*dsqrt(sigsq)
1181 rij_shift=1.0D0/rij-sig+r0ij
1182 C I hate to put IF's in the loops, but here don't have another choice!!!!
1183 if (rij_shift.le.0.0D0) then
1188 c---------------------------------------------------------------
1189 rij_shift=1.0D0/rij_shift
1190 fac=rij_shift**expon
1193 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1194 eps2der=evdwij*eps3rt
1195 eps3der=evdwij*eps2rt
1196 fac_augm=rrij**expon
1197 e_augm=augm(itypi,itypj)*fac_augm
1198 evdwij=evdwij*eps2rt*eps3rt
1199 evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
1201 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1203 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1204 & restyp(itypi),i,restyp(itypj),j,
1205 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1206 & chi1,chi2,chip1,chip2,
1207 & eps1,eps2rt**2,eps3rt**2,
1208 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1211 C Calculate gradient components.
1212 e1=e1*eps1*eps2rt**2*eps3rt**2
1213 fac=-expon*(e1+evdwij)*rij_shift
1215 fac=rij*fac-2*expon*rrij*e_augm
1216 C Calculate the radial part of the gradient
1220 C Calculate angular part of the gradient.
1221 call sc_grad_scale(1.0d0-sss)
1227 C-----------------------------------------------------------------------------
1228 subroutine egbv_short(evdw)
1230 C This subroutine calculates the interaction energy of nonbonded side chains
1231 C assuming the Gay-Berne-Vorobjev potential of interaction.
1233 implicit real*8 (a-h,o-z)
1234 include 'DIMENSIONS'
1235 include 'COMMON.GEO'
1236 include 'COMMON.VAR'
1237 include 'COMMON.LOCAL'
1238 include 'COMMON.CHAIN'
1239 include 'COMMON.DERIV'
1240 include 'COMMON.NAMES'
1241 include 'COMMON.INTERACT'
1242 include 'COMMON.IOUNITS'
1243 include 'COMMON.CALC'
1244 common /srutu/ icall
1247 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1250 c if (icall.eq.0) lprn=.true.
1252 do i=iatsc_s,iatsc_e
1254 if (itypi.eq.ntyp1) cycle
1259 dxi=dc_norm(1,nres+i)
1260 dyi=dc_norm(2,nres+i)
1261 dzi=dc_norm(3,nres+i)
1262 c dsci_inv=dsc_inv(itypi)
1263 dsci_inv=vbld_inv(i+nres)
1265 C Calculate SC interaction energy.
1267 do iint=1,nint_gr(i)
1268 do j=istart(i,iint),iend(i,iint)
1271 if (itypj.eq.ntyp1) cycle
1272 c dscj_inv=dsc_inv(itypj)
1273 dscj_inv=vbld_inv(j+nres)
1274 sig0ij=sigma(itypi,itypj)
1275 r0ij=r0(itypi,itypj)
1276 chi1=chi(itypi,itypj)
1277 chi2=chi(itypj,itypi)
1284 alf12=0.5D0*(alf1+alf2)
1288 dxj=dc_norm(1,nres+j)
1289 dyj=dc_norm(2,nres+j)
1290 dzj=dc_norm(3,nres+j)
1291 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1294 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
1296 if (sss.gt.0.0d0) then
1298 C Calculate angle-dependent terms of energy and contributions to their
1302 sig=sig0ij*dsqrt(sigsq)
1303 rij_shift=1.0D0/rij-sig+r0ij
1304 C I hate to put IF's in the loops, but here don't have another choice!!!!
1305 if (rij_shift.le.0.0D0) then
1310 c---------------------------------------------------------------
1311 rij_shift=1.0D0/rij_shift
1312 fac=rij_shift**expon
1315 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1316 eps2der=evdwij*eps3rt
1317 eps3der=evdwij*eps2rt
1318 fac_augm=rrij**expon
1319 e_augm=augm(itypi,itypj)*fac_augm
1320 evdwij=evdwij*eps2rt*eps3rt
1321 evdw=evdw+(evdwij+e_augm)*sss
1323 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1325 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1326 & restyp(itypi),i,restyp(itypj),j,
1327 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1328 & chi1,chi2,chip1,chip2,
1329 & eps1,eps2rt**2,eps3rt**2,
1330 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1333 C Calculate gradient components.
1334 e1=e1*eps1*eps2rt**2*eps3rt**2
1335 fac=-expon*(e1+evdwij)*rij_shift
1337 fac=rij*fac-2*expon*rrij*e_augm
1338 C Calculate the radial part of the gradient
1342 C Calculate angular part of the gradient.
1343 call sc_grad_scale(sss)
1349 C----------------------------------------------------------------------------
1350 subroutine sc_grad_scale(scalfac)
1351 implicit real*8 (a-h,o-z)
1352 include 'DIMENSIONS'
1353 include 'COMMON.CHAIN'
1354 include 'COMMON.DERIV'
1355 include 'COMMON.CALC'
1356 include 'COMMON.IOUNITS'
1357 double precision dcosom1(3),dcosom2(3)
1358 double precision scalfac
1359 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1360 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1361 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1362 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1366 c eom12=evdwij*eps1_om12
1368 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1369 c & " sigder",sigder
1370 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1371 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1373 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1374 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1377 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1378 gg_lipi(k)=gg_lipi(k)*scalfac
1379 gg_lipj(k)=gg_lipj(k)*scalfac
1381 c write (iout,*) "gg",(gg(k),k=1,3)
1383 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
1384 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1385 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1386 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
1387 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1388 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1389 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1390 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1391 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1392 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1395 C Calculate the components of the gradient in DC and X
1398 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
1399 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
1403 C--------------------------------------------------------------------------
1404 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1406 C This subroutine calculates the average interaction energy and its gradient
1407 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1408 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1409 C The potential depends both on the distance of peptide-group centers and on
1410 C the orientation of the CA-CA virtual bonds.
1412 implicit real*8 (a-h,o-z)
1416 include 'DIMENSIONS'
1417 include 'COMMON.CONTROL'
1418 include 'COMMON.SETUP'
1419 include 'COMMON.IOUNITS'
1420 include 'COMMON.GEO'
1421 include 'COMMON.VAR'
1422 include 'COMMON.LOCAL'
1423 include 'COMMON.CHAIN'
1424 include 'COMMON.DERIV'
1425 include 'COMMON.INTERACT'
1426 include 'COMMON.CONTACTS'
1427 include 'COMMON.TORSION'
1428 include 'COMMON.VECTORS'
1429 include 'COMMON.FFIELD'
1430 include 'COMMON.TIME1'
1431 include 'COMMON.SHIELD'
1432 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1433 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1434 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1435 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1436 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1437 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1439 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1441 double precision scal_el /1.0d0/
1443 double precision scal_el /0.5d0/
1446 C 13-go grudnia roku pamietnego...
1447 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1448 & 0.0d0,1.0d0,0.0d0,
1449 & 0.0d0,0.0d0,1.0d0/
1450 cd write(iout,*) 'In EELEC'
1452 cd write(iout,*) 'Type',i
1453 cd write(iout,*) 'B1',B1(:,i)
1454 cd write(iout,*) 'B2',B2(:,i)
1455 cd write(iout,*) 'CC',CC(:,:,i)
1456 cd write(iout,*) 'DD',DD(:,:,i)
1457 cd write(iout,*) 'EE',EE(:,:,i)
1459 cd call check_vecgrad
1461 C print *,"WCHODZE3"
1462 if (icheckgrad.eq.1) then
1464 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1466 dc_norm(k,i)=dc(k,i)*fac
1468 c write (iout,*) 'i',i,' fac',fac
1471 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1472 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1473 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1474 c call vec_and_deriv
1480 time_mat=time_mat+MPI_Wtime()-time01
1484 cd write (iout,*) 'i=',i
1486 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1489 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1490 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1503 cd print '(a)','Enter EELEC'
1504 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1506 gel_loc_loc(i)=0.0d0
1511 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1513 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1515 do i=iturn3_start,iturn3_end
1516 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1517 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1518 C & .or. itype(i-1).eq.ntyp1
1519 C & .or. itype(i+4).eq.ntyp1
1524 dx_normi=dc_norm(1,i)
1525 dy_normi=dc_norm(2,i)
1526 dz_normi=dc_norm(3,i)
1527 xmedi=c(1,i)+0.5d0*dxi
1528 ymedi=c(2,i)+0.5d0*dyi
1529 zmedi=c(3,i)+0.5d0*dzi
1530 xmedi=mod(xmedi,boxxsize)
1531 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1532 ymedi=mod(ymedi,boxysize)
1533 if (ymedi.lt.0) ymedi=ymedi+boxysize
1534 zmedi=mod(zmedi,boxzsize)
1535 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1537 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1538 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1539 num_cont_hb(i)=num_conti
1541 do i=iturn4_start,iturn4_end
1542 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1543 & .or. itype(i+3).eq.ntyp1
1544 & .or. itype(i+4).eq.ntyp1
1545 C & .or. itype(i+5).eq.ntyp1
1546 C & .or. itype(i-1).eq.ntyp1
1551 dx_normi=dc_norm(1,i)
1552 dy_normi=dc_norm(2,i)
1553 dz_normi=dc_norm(3,i)
1554 xmedi=c(1,i)+0.5d0*dxi
1555 ymedi=c(2,i)+0.5d0*dyi
1556 zmedi=c(3,i)+0.5d0*dzi
1557 xmedi=mod(xmedi,boxxsize)
1558 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1559 ymedi=mod(ymedi,boxysize)
1560 if (ymedi.lt.0) ymedi=ymedi+boxysize
1561 zmedi=mod(zmedi,boxzsize)
1562 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1563 num_conti=num_cont_hb(i)
1564 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1565 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1566 & call eturn4(i,eello_turn4)
1567 num_cont_hb(i)=num_conti
1570 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1572 do i=iatel_s,iatel_e
1573 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1574 C & .or. itype(i+2).eq.ntyp1
1575 C & .or. itype(i-1).eq.ntyp1
1580 dx_normi=dc_norm(1,i)
1581 dy_normi=dc_norm(2,i)
1582 dz_normi=dc_norm(3,i)
1583 xmedi=c(1,i)+0.5d0*dxi
1584 ymedi=c(2,i)+0.5d0*dyi
1585 zmedi=c(3,i)+0.5d0*dzi
1586 xmedi=mod(xmedi,boxxsize)
1587 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1588 ymedi=mod(ymedi,boxysize)
1589 if (ymedi.lt.0) ymedi=ymedi+boxysize
1590 zmedi=mod(zmedi,boxzsize)
1591 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1592 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1593 num_conti=num_cont_hb(i)
1594 do j=ielstart(i),ielend(i)
1595 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1596 C & .or.itype(j+2).eq.ntyp1
1597 C & .or.itype(j-1).eq.ntyp1
1599 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1601 num_cont_hb(i)=num_conti
1603 c write (iout,*) "Number of loop steps in EELEC:",ind
1605 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1606 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1608 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1609 ccc eel_loc=eel_loc+eello_turn3
1610 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1613 C-------------------------------------------------------------------------------
1614 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1615 implicit real*8 (a-h,o-z)
1616 include 'DIMENSIONS'
1620 include 'COMMON.CONTROL'
1621 include 'COMMON.IOUNITS'
1622 include 'COMMON.GEO'
1623 include 'COMMON.VAR'
1624 include 'COMMON.LOCAL'
1625 include 'COMMON.CHAIN'
1626 include 'COMMON.DERIV'
1627 include 'COMMON.INTERACT'
1628 include 'COMMON.CONTACTS'
1629 include 'COMMON.TORSION'
1630 include 'COMMON.VECTORS'
1631 include 'COMMON.FFIELD'
1632 include 'COMMON.TIME1'
1633 include 'COMMON.SHIELD'
1634 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1635 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1636 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1637 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1638 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1639 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1641 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1643 double precision scal_el /1.0d0/
1645 double precision scal_el /0.5d0/
1648 C 13-go grudnia roku pamietnego...
1649 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1650 & 0.0d0,1.0d0,0.0d0,
1651 & 0.0d0,0.0d0,1.0d0/
1652 integer xhift,yshift,zshift
1653 c time00=MPI_Wtime()
1654 cd write (iout,*) "eelecij",i,j
1655 C print *,"WCHODZE2"
1659 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1660 aaa=app(iteli,itelj)
1661 bbb=bpp(iteli,itelj)
1662 ael6i=ael6(iteli,itelj)
1663 ael3i=ael3(iteli,itelj)
1667 dx_normj=dc_norm(1,j)
1668 dy_normj=dc_norm(2,j)
1669 dz_normj=dc_norm(3,j)
1674 if (xj.lt.0) xj=xj+boxxsize
1676 if (yj.lt.0) yj=yj+boxysize
1678 if (zj.lt.0) zj=zj+boxzsize
1679 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1687 xj=xj_safe+xshift*boxxsize
1688 yj=yj_safe+yshift*boxysize
1689 zj=zj_safe+zshift*boxzsize
1690 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1691 if(dist_temp.lt.dist_init) then
1701 if (isubchap.eq.1) then
1711 rij=xj*xj+yj*yj+zj*zj
1715 c For extracting the short-range part of Evdwpp
1717 sssgrad=sscagrad(rij)
1720 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1721 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1722 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1723 fac=cosa-3.0D0*cosb*cosg
1725 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1726 if (j.eq.i+2) ev1=scal_el*ev1
1731 if (shield_mode.eq.0) then
1735 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1737 el1=el1*fac_shield(i)**2*fac_shield(j)**2
1738 el2=el2*fac_shield(i)**2*fac_shield(j)**2
1740 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1741 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1743 evdw1=evdw1+evdwij*(1.0d0-sss)
1744 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1745 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1746 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1747 cd & xmedi,ymedi,zmedi,xj,yj,zj
1749 if (energy_dec) then
1750 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
1751 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1755 C Calculate contributions to the Cartesian gradient.
1758 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
1759 facel=-3*rrmij*(el1+eesij)
1765 * Radial derivatives. First process both termini of the fragment (i,j)
1770 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1771 & (shield_mode.gt.0)) then
1773 do ilist=1,ishield_list(i)
1774 iresshield=shield_list(ilist,i)
1776 rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
1778 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1780 & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
1781 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1782 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1783 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1784 C if (iresshield.gt.i) then
1785 C do ishi=i+1,iresshield-1
1786 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1787 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1791 C do ishi=iresshield,i
1792 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1793 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1799 do ilist=1,ishield_list(j)
1800 iresshield=shield_list(ilist,j)
1802 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1804 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1806 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
1807 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1812 gshieldc(k,i)=gshieldc(k,i)+
1813 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
1814 gshieldc(k,j)=gshieldc(k,j)+
1815 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
1816 gshieldc(k,i-1)=gshieldc(k,i-1)+
1817 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
1818 gshieldc(k,j-1)=gshieldc(k,j-1)+
1819 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
1825 c ghalf=0.5D0*ggg(k)
1826 c gelc(k,i)=gelc(k,i)+ghalf
1827 c gelc(k,j)=gelc(k,j)+ghalf
1829 c 9/28/08 AL Gradient compotents will be summed only at the end
1831 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1832 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1835 * Loop over residues i+1 thru j-1.
1839 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1842 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj
1843 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj
1844 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj
1846 c ghalf=0.5D0*ggg(k)
1847 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1848 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1850 c 9/28/08 AL Gradient compotents will be summed only at the end
1852 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1853 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1856 * Loop over residues i+1 thru j-1.
1860 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1864 facvdw=ev1+evdwij*(1.0d0-sss)
1867 fac=-3*rrmij*(facvdw+facvdw+facel)
1872 * Radial derivatives. First process both termini of the fragment (i,j)
1878 c ghalf=0.5D0*ggg(k)
1879 c gelc(k,i)=gelc(k,i)+ghalf
1880 c gelc(k,j)=gelc(k,j)+ghalf
1882 c 9/28/08 AL Gradient compotents will be summed only at the end
1884 gelc_long(k,j)=gelc(k,j)+ggg(k)
1885 gelc_long(k,i)=gelc(k,i)-ggg(k)
1888 * Loop over residues i+1 thru j-1.
1892 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1895 c 9/28/08 AL Gradient compotents will be summed only at the end
1899 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj
1900 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj
1901 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj
1903 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1904 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1910 ecosa=2.0D0*fac3*fac1+fac4
1913 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1914 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1916 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1917 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1919 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1920 cd & (dcosg(k),k=1,3)
1922 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
1923 & fac_shield(i)**2*fac_shield(j)**2
1927 c ghalf=0.5D0*ggg(k)
1928 c gelc(k,i)=gelc(k,i)+ghalf
1929 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1930 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1931 c gelc(k,j)=gelc(k,j)+ghalf
1932 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1933 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1937 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1942 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1943 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
1944 & *fac_shield(i)**2*fac_shield(j)**2
1947 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1948 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
1949 & *fac_shield(i)**2*fac_shield(j)**2
1950 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1951 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1953 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1954 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1955 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1957 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1958 C energy of a peptide unit is assumed in the form of a second-order
1959 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1960 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1961 C are computed for EVERY pair of non-contiguous peptide groups.
1963 if (j.lt.nres-1) then
1974 muij(kkk)=mu(k,i)*mu(l,j)
1977 cd write (iout,*) 'EELEC: i',i,' j',j
1978 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1979 cd write(iout,*) 'muij',muij
1980 ury=scalar(uy(1,i),erij)
1981 urz=scalar(uz(1,i),erij)
1982 vry=scalar(uy(1,j),erij)
1983 vrz=scalar(uz(1,j),erij)
1984 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1985 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1986 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1987 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1988 fac=dsqrt(-ael6i)*r3ij
1993 cd write (iout,'(4i5,4f10.5)')
1994 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1995 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1996 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
1997 cd & uy(:,j),uz(:,j)
1998 cd write (iout,'(4f10.5)')
1999 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
2000 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
2001 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
2002 cd write (iout,'(9f10.5/)')
2003 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
2004 C Derivatives of the elements of A in virtual-bond vectors
2005 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
2007 uryg(k,1)=scalar(erder(1,k),uy(1,i))
2008 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
2009 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
2010 urzg(k,1)=scalar(erder(1,k),uz(1,i))
2011 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
2012 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
2013 vryg(k,1)=scalar(erder(1,k),uy(1,j))
2014 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
2015 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
2016 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
2017 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
2018 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
2020 C Compute radial contributions to the gradient
2038 C Add the contributions coming from er
2041 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
2042 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
2043 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
2044 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
2047 C Derivatives in DC(i)
2048 cgrad ghalf1=0.5d0*agg(k,1)
2049 cgrad ghalf2=0.5d0*agg(k,2)
2050 cgrad ghalf3=0.5d0*agg(k,3)
2051 cgrad ghalf4=0.5d0*agg(k,4)
2052 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2053 & -3.0d0*uryg(k,2)*vry)!+ghalf1
2054 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2055 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
2056 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2057 & -3.0d0*urzg(k,2)*vry)!+ghalf3
2058 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2059 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
2060 C Derivatives in DC(i+1)
2061 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2062 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2063 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2064 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2065 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2066 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2067 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2068 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2069 C Derivatives in DC(j)
2070 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2071 & -3.0d0*vryg(k,2)*ury)!+ghalf1
2072 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2073 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
2074 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2075 & -3.0d0*vryg(k,2)*urz)!+ghalf3
2076 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2077 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
2078 C Derivatives in DC(j+1) or DC(nres-1)
2079 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2080 & -3.0d0*vryg(k,3)*ury)
2081 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2082 & -3.0d0*vrzg(k,3)*ury)
2083 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2084 & -3.0d0*vryg(k,3)*urz)
2085 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2086 & -3.0d0*vrzg(k,3)*urz)
2087 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
2089 cgrad aggj1(k,l)=aggj1(k,l)+agg(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)
2108 if (j.lt.nres-1) then
2114 aggi(k,l)=-aggi(k,l)
2115 aggi1(k,l)=-aggi1(k,l)
2116 aggj(k,l)=-aggj(k,l)
2117 aggj1(k,l)=-aggj1(k,l)
2128 aggi(k,l)=-aggi(k,l)
2129 aggi1(k,l)=-aggi1(k,l)
2130 aggj(k,l)=-aggj(k,l)
2131 aggj1(k,l)=-aggj1(k,l)
2136 IF (wel_loc.gt.0.0d0) THEN
2137 C Contribution to the local-electrostatic energy coming from the i-j pair
2138 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2140 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2142 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2143 & 'eelloc',i,j,eel_loc_ij
2146 if (shield_mode.eq.0) then
2153 eel_loc_ij=eel_loc_ij
2154 & *fac_shield(i)*fac_shield(j)
2155 eel_loc=eel_loc+eel_loc_ij
2157 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2158 & (shield_mode.gt.0)) then
2161 do ilist=1,ishield_list(i)
2162 iresshield=shield_list(ilist,i)
2164 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2167 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2169 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2170 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2174 do ilist=1,ishield_list(j)
2175 iresshield=shield_list(ilist,j)
2177 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2180 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2182 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2183 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2190 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2191 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2192 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2193 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2194 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2195 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2196 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2197 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2201 C Partial derivatives in virtual-bond dihedral angles gamma
2203 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2204 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2205 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2206 & *fac_shield(i)*fac_shield(j)
2208 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2209 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2210 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2211 & *fac_shield(i)*fac_shield(j)
2213 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2215 ggg(l)=(agg(l,1)*muij(1)+
2216 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2217 & *fac_shield(i)*fac_shield(j)
2219 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2220 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2221 cgrad ghalf=0.5d0*ggg(l)
2222 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
2223 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
2227 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2230 C Remaining derivatives of eello
2232 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2233 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2234 & *fac_shield(i)*fac_shield(j)
2236 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2237 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2238 & *fac_shield(i)*fac_shield(j)
2240 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2241 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2242 & *fac_shield(i)*fac_shield(j)
2244 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2245 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2246 & *fac_shield(i)*fac_shield(j)
2250 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2251 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
2252 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2253 & .and. num_conti.le.maxconts) then
2254 c write (iout,*) i,j," entered corr"
2256 C Calculate the contact function. The ith column of the array JCONT will
2257 C contain the numbers of atoms that make contacts with the atom I (of numbers
2258 C greater than I). The arrays FACONT and GACONT will contain the values of
2259 C the contact function and its derivative.
2260 c r0ij=1.02D0*rpp(iteli,itelj)
2261 c r0ij=1.11D0*rpp(iteli,itelj)
2262 r0ij=2.20D0*rpp(iteli,itelj)
2263 c r0ij=1.55D0*rpp(iteli,itelj)
2264 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2265 if (fcont.gt.0.0D0) then
2266 num_conti=num_conti+1
2267 if (num_conti.gt.maxconts) then
2268 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2269 & ' will skip next contacts for this conf.'
2271 jcont_hb(num_conti,i)=j
2272 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
2273 cd & " jcont_hb",jcont_hb(num_conti,i)
2274 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2275 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2276 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2278 d_cont(num_conti,i)=rij
2279 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2280 C --- Electrostatic-interaction matrix ---
2281 a_chuj(1,1,num_conti,i)=a22
2282 a_chuj(1,2,num_conti,i)=a23
2283 a_chuj(2,1,num_conti,i)=a32
2284 a_chuj(2,2,num_conti,i)=a33
2285 C --- Gradient of rij
2287 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2294 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2295 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2296 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2297 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2298 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2303 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2304 C Calculate contact energies
2306 wij=cosa-3.0D0*cosb*cosg
2309 c fac3=dsqrt(-ael6i)/r0ij**3
2310 fac3=dsqrt(-ael6i)*r3ij
2311 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2312 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2313 if (ees0tmp.gt.0) then
2314 ees0pij=dsqrt(ees0tmp)
2318 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2319 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2320 if (ees0tmp.gt.0) then
2321 ees0mij=dsqrt(ees0tmp)
2326 if (shield_mode.eq.0) then
2330 ees0plist(num_conti,i)=j
2331 C fac_shield(i)=0.4d0
2332 C fac_shield(j)=0.6d0
2334 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2335 & *fac_shield(i)*fac_shield(j)
2336 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2337 & *fac_shield(i)*fac_shield(j)
2339 C Diagnostics. Comment out or remove after debugging!
2340 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2341 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2342 c ees0m(num_conti,i)=0.0D0
2344 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2345 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2346 C Angular derivatives of the contact function
2347 ees0pij1=fac3/ees0pij
2348 ees0mij1=fac3/ees0mij
2349 fac3p=-3.0D0*fac3*rrmij
2350 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2351 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2353 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2354 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2355 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2356 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2357 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2358 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2359 ecosap=ecosa1+ecosa2
2360 ecosbp=ecosb1+ecosb2
2361 ecosgp=ecosg1+ecosg2
2362 ecosam=ecosa1-ecosa2
2363 ecosbm=ecosb1-ecosb2
2364 ecosgm=ecosg1-ecosg2
2373 facont_hb(num_conti,i)=fcont
2374 fprimcont=fprimcont/rij
2375 cd facont_hb(num_conti,i)=1.0D0
2376 C Following line is for diagnostics.
2379 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2380 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2383 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2384 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2386 gggp(1)=gggp(1)+ees0pijp*xj
2387 gggp(2)=gggp(2)+ees0pijp*yj
2388 gggp(3)=gggp(3)+ees0pijp*zj
2389 gggm(1)=gggm(1)+ees0mijp*xj
2390 gggm(2)=gggm(2)+ees0mijp*yj
2391 gggm(3)=gggm(3)+ees0mijp*zj
2392 C Derivatives due to the contact function
2393 gacont_hbr(1,num_conti,i)=fprimcont*xj
2394 gacont_hbr(2,num_conti,i)=fprimcont*yj
2395 gacont_hbr(3,num_conti,i)=fprimcont*zj
2398 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2399 c following the change of gradient-summation algorithm.
2401 cgrad ghalfp=0.5D0*gggp(k)
2402 cgrad ghalfm=0.5D0*gggm(k)
2403 gacontp_hb1(k,num_conti,i)=!ghalfp
2404 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2405 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2406 & *fac_shield(i)*fac_shield(j)
2408 gacontp_hb2(k,num_conti,i)=!ghalfp
2409 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2410 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2411 & *fac_shield(i)*fac_shield(j)
2413 gacontp_hb3(k,num_conti,i)=gggp(k)
2414 & *fac_shield(i)*fac_shield(j)
2416 gacontm_hb1(k,num_conti,i)=!ghalfm
2417 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2418 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2419 & *fac_shield(i)*fac_shield(j)
2421 gacontm_hb2(k,num_conti,i)=!ghalfm
2422 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2423 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2424 & *fac_shield(i)*fac_shield(j)
2426 gacontm_hb3(k,num_conti,i)=gggm(k)
2427 & *fac_shield(i)*fac_shield(j)
2431 endif ! num_conti.le.maxconts
2434 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2437 ghalf=0.5d0*agg(l,k)
2438 aggi(l,k)=aggi(l,k)+ghalf
2439 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2440 aggj(l,k)=aggj(l,k)+ghalf
2443 if (j.eq.nres-1 .and. i.lt.j-2) then
2446 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2451 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2454 C-----------------------------------------------------------------------
2455 subroutine evdwpp_short(evdw1)
2459 implicit real*8 (a-h,o-z)
2460 include 'DIMENSIONS'
2461 include 'COMMON.CONTROL'
2462 include 'COMMON.IOUNITS'
2463 include 'COMMON.GEO'
2464 include 'COMMON.VAR'
2465 include 'COMMON.LOCAL'
2466 include 'COMMON.CHAIN'
2467 include 'COMMON.DERIV'
2468 include 'COMMON.INTERACT'
2469 include 'COMMON.CONTACTS'
2470 include 'COMMON.TORSION'
2471 include 'COMMON.VECTORS'
2472 include 'COMMON.FFIELD'
2474 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2476 double precision scal_el /1.0d0/
2478 double precision scal_el /0.5d0/
2482 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2483 c & " iatel_e_vdw",iatel_e_vdw
2485 do i=iatel_s_vdw,iatel_e_vdw
2486 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2490 dx_normi=dc_norm(1,i)
2491 dy_normi=dc_norm(2,i)
2492 dz_normi=dc_norm(3,i)
2493 xmedi=c(1,i)+0.5d0*dxi
2494 ymedi=c(2,i)+0.5d0*dyi
2495 zmedi=c(3,i)+0.5d0*dzi
2496 xmedi=mod(xmedi,boxxsize)
2497 if (xmedi.lt.0) xmedi=xmedi+boxxsize
2498 ymedi=mod(ymedi,boxysize)
2499 if (ymedi.lt.0) ymedi=ymedi+boxysize
2500 zmedi=mod(zmedi,boxzsize)
2501 if (zmedi.lt.0) zmedi=zmedi+boxzsize
2503 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2504 c & ' ielend',ielend_vdw(i)
2506 do j=ielstart_vdw(i),ielend_vdw(i)
2507 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2511 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2512 aaa=app(iteli,itelj)
2513 bbb=bpp(iteli,itelj)
2517 dx_normj=dc_norm(1,j)
2518 dy_normj=dc_norm(2,j)
2519 dz_normj=dc_norm(3,j)
2524 if (xj.lt.0) xj=xj+boxxsize
2526 if (yj.lt.0) yj=yj+boxysize
2528 if (zj.lt.0) zj=zj+boxzsize
2529 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2537 xj=xj_safe+xshift*boxxsize
2538 yj=yj_safe+yshift*boxysize
2539 zj=zj_safe+zshift*boxzsize
2540 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2541 if(dist_temp.lt.dist_init) then
2551 if (isubchap.eq.1) then
2560 rij=xj*xj+yj*yj+zj*zj
2564 sssgrad=sscagrad(rij)
2565 if (sss.gt.0.0d0) then
2570 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2571 if (j.eq.i+2) ev1=scal_el*ev1
2574 if (energy_dec) then
2575 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2577 evdw1=evdw1+evdwij*sss
2579 C Calculate contributions to the Cartesian gradient.
2581 facvdw=-6*rrmij*(ev1+evdwij)*sss
2582 ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
2583 ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
2584 ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
2589 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2590 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2597 C-----------------------------------------------------------------------------
2598 subroutine escp_long(evdw2,evdw2_14)
2600 C This subroutine calculates the excluded-volume interaction energy between
2601 C peptide-group centers and side chains and its gradient in virtual-bond and
2602 C side-chain vectors.
2604 implicit real*8 (a-h,o-z)
2605 include 'DIMENSIONS'
2606 include 'COMMON.GEO'
2607 include 'COMMON.VAR'
2608 include 'COMMON.LOCAL'
2609 include 'COMMON.CHAIN'
2610 include 'COMMON.DERIV'
2611 include 'COMMON.INTERACT'
2612 include 'COMMON.FFIELD'
2613 include 'COMMON.IOUNITS'
2614 include 'COMMON.CONTROL'
2616 integer xshift,yshift,zshift
2619 CD print '(a)','Enter ESCP KURWA'
2620 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2621 do i=iatscp_s,iatscp_e
2622 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2624 xi=0.5D0*(c(1,i)+c(1,i+1))
2625 yi=0.5D0*(c(2,i)+c(2,i+1))
2626 zi=0.5D0*(c(3,i)+c(3,i+1))
2628 if (xi.lt.0) xi=xi+boxxsize
2630 if (yi.lt.0) yi=yi+boxysize
2632 if (zi.lt.0) zi=zi+boxzsize
2633 do iint=1,nscp_gr(i)
2635 do j=iscpstart(i,iint),iscpend(i,iint)
2637 if (itypj.eq.ntyp1) cycle
2638 C Uncomment following three lines for SC-p interactions
2642 C Uncomment following three lines for Ca-p interactions
2647 if (xj.lt.0) xj=xj+boxxsize
2649 if (yj.lt.0) yj=yj+boxysize
2651 if (zj.lt.0) zj=zj+boxzsize
2653 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2661 xj=xj_safe+xshift*boxxsize
2662 yj=yj_safe+yshift*boxysize
2663 zj=zj_safe+zshift*boxzsize
2664 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2665 if(dist_temp.lt.dist_init) then
2675 if (subchap.eq.1) then
2685 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2686 sss=sscale(1.0d0/(dsqrt(rrij)))
2687 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
2688 if (sss.lt.1.0d0) then
2690 e1=fac*fac*aad(itypj,iteli)
2691 e2=fac*bad(itypj,iteli)
2692 if (iabs(j-i) .le. 2) then
2695 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
2698 evdw2=evdw2+evdwij*(1.0d0-sss)
2699 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2700 & 'evdw2',i,j,sss,evdwij
2702 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2705 fac=-(evdwij+e1)*rrij*(1.0d0-sss)
2706 fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/expon
2710 C Uncomment following three lines for SC-p interactions
2712 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2714 C Uncomment following line for SC-p interactions
2715 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2717 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2718 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2727 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2728 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2729 gradx_scp(j,i)=expon*gradx_scp(j,i)
2732 C******************************************************************************
2736 C To save time the factor EXPON has been extracted from ALL components
2737 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2740 C******************************************************************************
2743 C-----------------------------------------------------------------------------
2744 subroutine escp_short(evdw2,evdw2_14)
2746 C This subroutine calculates the excluded-volume interaction energy between
2747 C peptide-group centers and side chains and its gradient in virtual-bond and
2748 C side-chain vectors.
2750 implicit real*8 (a-h,o-z)
2751 include 'DIMENSIONS'
2752 include 'COMMON.GEO'
2753 include 'COMMON.VAR'
2754 include 'COMMON.LOCAL'
2755 include 'COMMON.CHAIN'
2756 include 'COMMON.DERIV'
2757 include 'COMMON.INTERACT'
2758 include 'COMMON.FFIELD'
2759 include 'COMMON.IOUNITS'
2760 include 'COMMON.CONTROL'
2762 integer xshift,yshift,zshift
2765 cd print '(a)','Enter ESCP'
2766 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2767 do i=iatscp_s,iatscp_e
2768 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2770 xi=0.5D0*(c(1,i)+c(1,i+1))
2771 yi=0.5D0*(c(2,i)+c(2,i+1))
2772 zi=0.5D0*(c(3,i)+c(3,i+1))
2774 if (xi.lt.0) xi=xi+boxxsize
2776 if (yi.lt.0) yi=yi+boxysize
2778 if (zi.lt.0) zi=zi+boxzsize
2780 do iint=1,nscp_gr(i)
2782 do j=iscpstart(i,iint),iscpend(i,iint)
2784 if (itypj.eq.ntyp1) cycle
2785 C Uncomment following three lines for SC-p interactions
2789 C Uncomment following three lines for Ca-p interactions
2794 if (xj.lt.0) xj=xj+boxxsize
2796 if (yj.lt.0) yj=yj+boxysize
2798 if (zj.lt.0) zj=zj+boxzsize
2799 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2807 xj=xj_safe+xshift*boxxsize
2808 yj=yj_safe+yshift*boxysize
2809 zj=zj_safe+zshift*boxzsize
2810 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2811 if(dist_temp.lt.dist_init) then
2821 if (subchap.eq.1) then
2830 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2831 sss=sscale(1.0d0/(dsqrt(rrij)))
2832 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
2833 if (sss.gt.0.0d0) then
2836 e1=fac*fac*aad(itypj,iteli)
2837 e2=fac*bad(itypj,iteli)
2838 if (iabs(j-i) .le. 2) then
2841 evdw2_14=evdw2_14+(e1+e2)*sss
2844 evdw2=evdw2+evdwij*sss
2845 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2846 & 'evdw2',i,j,sss,evdwij
2849 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2851 fac=-(evdwij+e1)*rrij*sss
2852 fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon
2856 C Uncomment following three lines for SC-p interactions
2858 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2860 C Uncomment following line for SC-p interactions
2861 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2863 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2864 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2873 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2874 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2875 gradx_scp(j,i)=expon*gradx_scp(j,i)
2878 C******************************************************************************
2882 C To save time the factor EXPON has been extracted from ALL components
2883 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2886 C******************************************************************************