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,r_cut)
35 double precision r,r_cut,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,r_cut)
50 double precision r,r_cut,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.
72 include 'COMMON.LOCAL'
73 include 'COMMON.CHAIN'
74 include 'COMMON.DERIV'
75 include 'COMMON.INTERACT'
76 include 'COMMON.TORSION'
77 include 'COMMON.SBRIDGE'
78 include 'COMMON.NAMES'
79 include 'COMMON.IOUNITS'
80 include "COMMON.SPLITELE"
81 c include 'COMMON.CONTACTS'
82 double precision gg(3)
83 double precision evdw,evdwij
84 integer i,j,k,itypi,itypj,itypi1,num_conti,iint
85 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
86 & sigij,r0ij,rcut,sss1,sssgrad1,sqrij
87 double precision sscale,sscagrad
88 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
92 if (itypi.eq.ntyp1) cycle
93 itypi1=iabs(itype(i+1))
98 C Calculate SC interaction energy.
101 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
102 cd & 'iend=',iend(i,iint)
103 do j=istart(i,iint),iend(i,iint)
105 if (itypj.eq.ntyp1) cycle
109 rij=xj*xj+yj*yj+zj*zj
111 eps0ij=eps(itypi,itypj)
112 sss1=sscale(sqrij,r_cut_int)
113 if (sss1.eq.0.0d0) cycle
114 sssgrad1=sscagrad(sqrij,r_cut_int)
116 & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
117 sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
118 if (sss.lt.1.0d0) then
124 evdw=evdw+(1.0d0-sss)*sss1*evdwij/sqrij/expon
126 C Calculate the components of the gradient in DC and X
128 fac=-rrij*(e1+evdwij)*(1.0d0-sss)*sss1
129 & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
130 & +(1.0d0-sss)*sssgrad1)/sqrij
135 gvdwx(k,i)=gvdwx(k,i)-gg(k)
136 gvdwx(k,j)=gvdwx(k,j)+gg(k)
137 gvdwc(k,i)=gvdwc(k,i)-gg(k)
138 gvdwc(k,j)=gvdwc(k,j)+gg(k)
146 gvdwc(j,i)=expon*gvdwc(j,i)
147 gvdwx(j,i)=expon*gvdwx(j,i)
150 C******************************************************************************
154 C To save time, the factor of EXPON has been extracted from ALL components
155 C of GVDWC and GRADX. Remember to multiply them by this factor before further
158 C******************************************************************************
161 C-----------------------------------------------------------------------
162 subroutine elj_short(evdw)
164 C This subroutine calculates the interaction energy of nonbonded side chains
165 C assuming the LJ potential of interaction.
171 include 'COMMON.LOCAL'
172 include 'COMMON.CHAIN'
173 include 'COMMON.DERIV'
174 include 'COMMON.INTERACT'
175 include 'COMMON.TORSION'
176 include 'COMMON.SBRIDGE'
177 include 'COMMON.NAMES'
178 include 'COMMON.IOUNITS'
179 include "COMMON.SPLITELE"
180 c include 'COMMON.CONTACTS'
181 double precision gg(3)
182 double precision evdw,evdwij
183 integer i,j,k,itypi,itypj,itypi1,num_conti,iint
184 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
185 & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
186 double precision sscale,sscagrad
187 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
191 if (itypi.eq.ntyp1) cycle
192 itypi1=iabs(itype(i+1))
199 C Calculate SC interaction energy.
202 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
203 cd & 'iend=',iend(i,iint)
204 do j=istart(i,iint),iend(i,iint)
206 if (itypj.eq.ntyp1) cycle
210 C Change 12/1/95 to calculate four-body interactions
211 rij=xj*xj+yj*yj+zj*zj
213 sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
214 if (sss.gt.0.0d0) then
216 & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
218 eps0ij=eps(itypi,itypj)
225 C Calculate the components of the gradient in DC and X
227 fac=-rrij*(e1+evdwij)*sss+evdwij*sssgrad/sqrij/expon
232 gvdwx(k,i)=gvdwx(k,i)-gg(k)
233 gvdwx(k,j)=gvdwx(k,j)+gg(k)
234 gvdwc(k,i)=gvdwc(k,i)-gg(k)
235 gvdwc(k,j)=gvdwc(k,j)+gg(k)
243 gvdwc(j,i)=expon*gvdwc(j,i)
244 gvdwx(j,i)=expon*gvdwx(j,i)
247 C******************************************************************************
251 C To save time, the factor of EXPON has been extracted from ALL components
252 C of GVDWC and GRADX. Remember to multiply them by this factor before further
255 C******************************************************************************
258 C-----------------------------------------------------------------------------
259 subroutine eljk_long(evdw)
261 C This subroutine calculates the interaction energy of nonbonded side chains
262 C assuming the LJK potential of interaction.
268 include 'COMMON.LOCAL'
269 include 'COMMON.CHAIN'
270 include 'COMMON.DERIV'
271 include 'COMMON.INTERACT'
272 include 'COMMON.IOUNITS'
273 include 'COMMON.NAMES'
274 include "COMMON.SPLITELE"
275 double precision gg(3)
276 double precision evdw,evdwij
277 integer i,j,k,itypi,itypj,itypi1,iint
278 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
279 & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
281 double precision sscale,sscagrad
282 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
286 if (itypi.eq.ntyp1) cycle
287 itypi1=iabs(itype(i+1))
292 C Calculate SC interaction energy.
295 do j=istart(i,iint),iend(i,iint)
297 if (itypj.eq.ntyp1) cycle
301 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
303 e_augm=augm(itypi,itypj)*fac_augm
306 sss1=sscale(rij,r_cut_int)
307 if (sss1.eq.0.0d0) cycle
308 sssgrad1=sscagrad(rij,r_cut_int)
309 sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
310 if (sss.lt.1.0d0) then
312 & sscagrad(rij/sigma(itypi,itypj),r_cut_respa)
313 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
314 fac=r_shift_inv**expon
318 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
319 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
320 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
321 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
322 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
323 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
324 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
325 evdw=evdw+(1.0d0-sss)*sss1*evdwij
327 C Calculate the components of the gradient in DC and X
329 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
330 fac=fac*(1.0d0-sss)*sss1
331 & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
332 & +(1.0d0-sss)*sssgrad1)*r_inv_ij/expon
337 gvdwx(k,i)=gvdwx(k,i)-gg(k)
338 gvdwx(k,j)=gvdwx(k,j)+gg(k)
339 gvdwc(k,i)=gvdwc(k,i)-gg(k)
340 gvdwc(k,j)=gvdwc(k,j)+gg(k)
348 gvdwc(j,i)=expon*gvdwc(j,i)
349 gvdwx(j,i)=expon*gvdwx(j,i)
354 C-----------------------------------------------------------------------------
355 subroutine eljk_short(evdw)
357 C This subroutine calculates the interaction energy of nonbonded side chains
358 C assuming the LJK potential of interaction.
360 implicit real*8 (a-h,o-z)
364 include 'COMMON.LOCAL'
365 include 'COMMON.CHAIN'
366 include 'COMMON.DERIV'
367 include 'COMMON.INTERACT'
368 include 'COMMON.IOUNITS'
369 include 'COMMON.NAMES'
370 include "COMMON.SPLITELE"
371 double precision gg(3)
372 double precision evdw,evdwij
373 integer i,j,k,itypi,itypj,itypi1,iint
374 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
375 & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
377 double precision sscale,sscagrad
378 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
382 if (itypi.eq.ntyp1) cycle
383 itypi1=iabs(itype(i+1))
388 C Calculate SC interaction energy.
391 do j=istart(i,iint),iend(i,iint)
393 if (itypj.eq.ntyp1) cycle
397 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
399 e_augm=augm(itypi,itypj)*fac_augm
402 sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
403 if (sss.gt.0.0d0) then
404 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
405 fac=r_shift_inv**expon
409 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
410 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
411 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
412 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
413 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
414 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
415 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
418 C Calculate the components of the gradient in DC and X
420 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
421 & +evdwij*sssgrad/sigma(itypi,itypj)*r_inv_ij/expon
427 gvdwx(k,i)=gvdwx(k,i)-gg(k)
428 gvdwx(k,j)=gvdwx(k,j)+gg(k)
429 gvdwc(k,i)=gvdwc(k,i)-gg(k)
430 gvdwc(k,j)=gvdwc(k,j)+gg(k)
438 gvdwc(j,i)=expon*gvdwc(j,i)
439 gvdwx(j,i)=expon*gvdwx(j,i)
444 C-----------------------------------------------------------------------------
445 subroutine ebp_long(evdw)
447 C This subroutine calculates the interaction energy of nonbonded side chains
448 C assuming the Berne-Pechukas potential of interaction.
454 include 'COMMON.LOCAL'
455 include 'COMMON.CHAIN'
456 include 'COMMON.DERIV'
457 include 'COMMON.NAMES'
458 include 'COMMON.INTERACT'
459 include 'COMMON.IOUNITS'
460 include 'COMMON.CALC'
461 include "COMMON.SPLITELE"
464 double precision evdw
465 integer itypi,itypj,itypi1,iint,ind
466 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
467 double precision sss1,sssgrad1
468 double precision sscale,sscagrad
469 c double precision rrsave(maxdim)
472 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
474 c if (icall.eq.0) then
482 if (itypi.eq.ntyp1) cycle
483 itypi1=iabs(itype(i+1))
487 dxi=dc_norm(1,nres+i)
488 dyi=dc_norm(2,nres+i)
489 dzi=dc_norm(3,nres+i)
490 c dsci_inv=dsc_inv(itypi)
491 dsci_inv=vbld_inv(i+nres)
493 C Calculate SC interaction energy.
496 do j=istart(i,iint),iend(i,iint)
499 if (itypj.eq.ntyp1) cycle
500 c dscj_inv=dsc_inv(itypj)
501 dscj_inv=vbld_inv(j+nres)
502 chi1=chi(itypi,itypj)
503 chi2=chi(itypj,itypi)
510 alf12=0.5D0*(alf1+alf2)
514 dxj=dc_norm(1,nres+j)
515 dyj=dc_norm(2,nres+j)
516 dzj=dc_norm(3,nres+j)
517 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
519 sss1=sscale(1.0d0/rij,r_cut_int)
520 if (sss1.eq.0.0d0) cycle
521 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
522 if (sss.lt.1.0d0) then
524 & sscagrad(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
525 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
526 C Calculate the angle-dependent terms of energy & contributions to derivatives.
528 C Calculate whole angle-dependent part of epsilon and contributions
530 fac=(rrij*sigsq)**expon2
533 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
534 eps2der=evdwij*eps3rt
535 eps3der=evdwij*eps2rt
536 evdwij=evdwij*eps2rt*eps3rt
537 evdw=evdw+evdwij*(1.0d0-sss)*sss1
539 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
541 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
542 cd & restyp(itypi),i,restyp(itypj),j,
543 cd & epsi,sigm,chi1,chi2,chip1,chip2,
544 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
545 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
548 C Calculate gradient components.
549 e1=e1*eps1*eps2rt**2*eps3rt**2
550 fac=-expon*(e1+evdwij)
552 fac=(fac+evdwij*(sss1/(1.0d0-sss)*sssgrad/
553 & sigmaii(itypi,itypj)+(1.0d0-sss)/sss1*sssgrad1))*rij
554 C Calculate radial part of the gradient
558 C Calculate the angular part of the gradient and sum add the contributions
559 C to the appropriate components of the Cartesian gradient.
560 call sc_grad_scale((1.0d0-sss)*sss1)
568 C-----------------------------------------------------------------------------
569 subroutine ebp_short(evdw)
571 C This subroutine calculates the interaction energy of nonbonded side chains
572 C assuming the Berne-Pechukas potential of interaction.
574 implicit real*8 (a-h,o-z)
578 include 'COMMON.LOCAL'
579 include 'COMMON.CHAIN'
580 include 'COMMON.DERIV'
581 include 'COMMON.NAMES'
582 include 'COMMON.INTERACT'
583 include 'COMMON.IOUNITS'
584 include 'COMMON.CALC'
585 include "COMMON.SPLITELE"
588 double precision evdw
589 integer itypi,itypj,itypi1,iint,ind
590 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
591 double precision sscale,sscagrad
592 c double precision rrsave(maxdim)
595 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
597 c if (icall.eq.0) then
605 if (itypi.eq.ntyp1) cycle
606 itypi1=iabs(itype(i+1))
610 dxi=dc_norm(1,nres+i)
611 dyi=dc_norm(2,nres+i)
612 dzi=dc_norm(3,nres+i)
613 c dsci_inv=dsc_inv(itypi)
614 dsci_inv=vbld_inv(i+nres)
616 C Calculate SC interaction energy.
619 do j=istart(i,iint),iend(i,iint)
622 if (itypj.eq.ntyp1) cycle
623 c dscj_inv=dsc_inv(itypj)
624 dscj_inv=vbld_inv(j+nres)
625 chi1=chi(itypi,itypj)
626 chi2=chi(itypj,itypi)
633 alf12=0.5D0*(alf1+alf2)
637 dxj=dc_norm(1,nres+j)
638 dyj=dc_norm(2,nres+j)
639 dzj=dc_norm(3,nres+j)
640 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
642 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
644 if (sss.gt.0.0d0) then
646 C Calculate the angle-dependent terms of energy & contributions to derivatives.
648 C Calculate whole angle-dependent part of epsilon and contributions
650 fac=(rrij*sigsq)**expon2
653 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
654 eps2der=evdwij*eps3rt
655 eps3der=evdwij*eps2rt
656 evdwij=evdwij*eps2rt*eps3rt
659 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
661 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
662 cd & restyp(itypi),i,restyp(itypj),j,
663 cd & epsi,sigm,chi1,chi2,chip1,chip2,
664 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
665 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
668 C Calculate gradient components.
669 e1=e1*eps1*eps2rt**2*eps3rt**2
670 fac=-expon*(e1+evdwij)
672 fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rrij
673 C Calculate radial part of the gradient
677 C Calculate the angular part of the gradient and sum add the contributions
678 C to the appropriate components of the Cartesian gradient.
679 call sc_grad_scale(sss)
687 C-----------------------------------------------------------------------------
688 subroutine egb_long(evdw)
690 C This subroutine calculates the interaction energy of nonbonded side chains
691 C assuming the Gay-Berne potential of interaction.
697 include 'COMMON.LOCAL'
698 include 'COMMON.CHAIN'
699 include 'COMMON.DERIV'
700 include 'COMMON.NAMES'
701 include 'COMMON.INTERACT'
702 include 'COMMON.IOUNITS'
703 include 'COMMON.CALC'
704 include 'COMMON.CONTROL'
705 include "COMMON.SPLITELE"
707 integer xshift,yshift,zshift
708 double precision evdw
709 integer itypi,itypj,itypi1,iint,ind
710 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
711 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
712 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
713 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
714 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
715 double precision subchap,sss1,sssgrad1
717 ccccc energy_dec=.false.
718 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
721 c if (icall.eq.0) lprn=.false.
725 if (itypi.eq.ntyp1) cycle
726 itypi1=iabs(itype(i+1))
731 if (xi.lt.0) xi=xi+boxxsize
733 if (yi.lt.0) yi=yi+boxysize
735 if (zi.lt.0) zi=zi+boxzsize
736 dxi=dc_norm(1,nres+i)
737 dyi=dc_norm(2,nres+i)
738 dzi=dc_norm(3,nres+i)
739 c dsci_inv=dsc_inv(itypi)
740 dsci_inv=vbld_inv(i+nres)
741 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
742 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
744 C Calculate SC interaction energy.
747 do j=istart(i,iint),iend(i,iint)
750 if (itypj.eq.ntyp1) cycle
751 c dscj_inv=dsc_inv(itypj)
752 dscj_inv=vbld_inv(j+nres)
753 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
754 c & 1.0d0/vbld(j+nres)
755 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
756 sig0ij=sigma(itypi,itypj)
757 chi1=chi(itypi,itypj)
758 chi2=chi(itypj,itypi)
765 alf12=0.5D0*(alf1+alf2)
770 if (xj.lt.0) xj=xj+boxxsize
772 if (yj.lt.0) yj=yj+boxysize
774 if (zj.lt.0) zj=zj+boxzsize
775 if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
776 C the energy transfer exist
777 if (zj.lt.buflipbot) then
778 C what fraction I am in
779 fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
780 C lipbufthick is thickenes of lipid buffore
781 sslipj=sscalelip(fracinbuf)
782 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
783 else if (zi.gt.bufliptop) then
784 fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
785 sslipj=sscalelip(fracinbuf)
786 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
795 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
796 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
797 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
798 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
799 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
807 xj=xj_safe+xshift*boxxsize
808 yj=yj_safe+yshift*boxysize
809 zj=zj_safe+zshift*boxzsize
810 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
811 if (dist_temp.lt.dist_init) then
821 if (subchap.eq.1) then
830 dxj=dc_norm(1,nres+j)
831 dyj=dc_norm(2,nres+j)
832 dzj=dc_norm(3,nres+j)
833 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
835 sss1=sscale(1.0d0/rij,r_cut_int)
836 if (sss1.eq.0.0d0) cycle
837 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
838 if (sss.lt.1.0d0) then
839 C Calculate angle-dependent terms of energy and contributions to their
842 & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
843 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
846 sig=sig0ij*dsqrt(sigsq)
847 rij_shift=1.0D0/rij-sig+sig0ij
848 c for diagnostics; uncomment
849 c rij_shift=1.2*sig0ij
850 C I hate to put IF's in the loops, but here don't have another choice!!!!
851 if (rij_shift.le.0.0D0) then
853 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
854 cd & restyp(itypi),i,restyp(itypj),j,
855 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
859 c---------------------------------------------------------------
860 rij_shift=1.0D0/rij_shift
864 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
865 eps2der=evdwij*eps3rt
866 eps3der=evdwij*eps2rt
867 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
868 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
869 evdwij=evdwij*eps2rt*eps3rt
870 evdw=evdw+evdwij*(1.0d0-sss)*sss1
872 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
874 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
875 & restyp(itypi),i,restyp(itypj),j,
876 & epsi,sigm,chi1,chi2,chip1,chip2,
877 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
878 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
882 if (energy_dec) write (iout,'(a6,2i5,4f10.5)')
883 & 'evdw',i,j,rij,sss,sss1,evdwij
885 C Calculate gradient components.
886 e1=e1*eps1*eps2rt**2*eps3rt**2
887 fac=-expon*(e1+evdwij)*rij_shift
889 fac=(fac+evdwij*(-sss1*sssgrad/(1.0d0-sss)
890 & /sigmaii(itypi,itypj)+(1.0d0-sss)*sssgrad1/sss1))*rij
892 C Calculate the radial part of the gradient
896 gg_lipi(3)=ssgradlipi*evdwij
897 gg_lipj(3)=ssgradlipj*evdwij
898 C Calculate angular part of the gradient.
899 call sc_grad_scale((1.0d0-sss)*sss1)
904 c write (iout,*) "Number of loop steps in EGB:",ind
905 cccc energy_dec=.false.
908 C-----------------------------------------------------------------------------
909 subroutine egb_short(evdw)
911 C This subroutine calculates the interaction energy of nonbonded side chains
912 C assuming the Gay-Berne potential of interaction.
918 include 'COMMON.LOCAL'
919 include 'COMMON.CHAIN'
920 include 'COMMON.DERIV'
921 include 'COMMON.NAMES'
922 include 'COMMON.INTERACT'
923 include 'COMMON.IOUNITS'
924 include 'COMMON.CALC'
925 include 'COMMON.CONTROL'
926 include "COMMON.SPLITELE"
928 integer xshift,yshift,zshift
929 double precision evdw
930 integer itypi,itypj,itypi1,iint,ind
931 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
932 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
933 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
934 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
935 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
936 double precision subchap
938 ccccc energy_dec=.false.
939 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
942 c if (icall.eq.0) lprn=.false.
946 if (itypi.eq.ntyp1) cycle
947 itypi1=iabs(itype(i+1))
952 if (xi.lt.0) xi=xi+boxxsize
954 if (yi.lt.0) yi=yi+boxysize
956 if (zi.lt.0) zi=zi+boxzsize
957 dxi=dc_norm(1,nres+i)
958 dyi=dc_norm(2,nres+i)
959 dzi=dc_norm(3,nres+i)
960 c dsci_inv=dsc_inv(itypi)
961 dsci_inv=vbld_inv(i+nres)
962 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
963 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
965 C Calculate SC interaction energy.
968 do j=istart(i,iint),iend(i,iint)
971 if (itypj.eq.ntyp1) cycle
972 c dscj_inv=dsc_inv(itypj)
973 dscj_inv=vbld_inv(j+nres)
974 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
975 c & 1.0d0/vbld(j+nres)
976 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
977 sig0ij=sigma(itypi,itypj)
978 chi1=chi(itypi,itypj)
979 chi2=chi(itypj,itypi)
986 alf12=0.5D0*(alf1+alf2)
991 if (xj.lt.0) xj=xj+boxxsize
993 if (yj.lt.0) yj=yj+boxysize
995 if (zj.lt.0) zj=zj+boxzsize
996 if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
997 C the energy transfer exist
998 if (zj.lt.buflipbot) then
999 C what fraction I am in
1000 fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
1001 C lipbufthick is thickenes of lipid buffore
1002 sslipj=sscalelip(fracinbuf)
1003 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
1004 elseif (zi.gt.bufliptop) then
1005 fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
1006 sslipj=sscalelip(fracinbuf)
1007 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
1016 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1017 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
1018 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1019 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
1020 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1028 xj=xj_safe+xshift*boxxsize
1029 yj=yj_safe+yshift*boxysize
1030 zj=zj_safe+zshift*boxzsize
1031 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1032 if(dist_temp.lt.dist_init) then
1042 if (subchap.eq.1) then
1051 dxj=dc_norm(1,nres+j)
1052 dyj=dc_norm(2,nres+j)
1053 dzj=dc_norm(3,nres+j)
1054 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1056 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1057 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1058 if (sss.gt.0.0d0) then
1060 C Calculate angle-dependent terms of energy and contributions to their
1064 sig=sig0ij*dsqrt(sigsq)
1065 rij_shift=1.0D0/rij-sig+sig0ij
1066 c for diagnostics; uncomment
1067 c rij_shift=1.2*sig0ij
1068 C I hate to put IF's in the loops, but here don't have another choice!!!!
1069 if (rij_shift.le.0.0D0) then
1071 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1072 cd & restyp(itypi),i,restyp(itypj),j,
1073 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1077 c---------------------------------------------------------------
1078 rij_shift=1.0D0/rij_shift
1079 fac=rij_shift**expon
1082 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1083 eps2der=evdwij*eps3rt
1084 eps3der=evdwij*eps2rt
1085 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1086 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1087 evdwij=evdwij*eps2rt*eps3rt
1088 evdw=evdw+evdwij*sss
1090 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1092 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1093 & restyp(itypi),i,restyp(itypj),j,
1094 & epsi,sigm,chi1,chi2,chip1,chip2,
1095 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1096 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1100 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1103 C Calculate gradient components.
1104 e1=e1*eps1*eps2rt**2*eps3rt**2
1105 fac=-expon*(e1+evdwij)*rij_shift
1107 fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rij
1109 C Calculate the radial part of the gradient
1113 gg_lipi(3)=ssgradlipi*evdwij
1114 gg_lipj(3)=ssgradlipj*evdwij
1115 C Calculate angular part of the gradient.
1116 call sc_grad_scale(sss)
1121 c write (iout,*) "Number of loop steps in EGB:",ind
1122 cccc energy_dec=.false.
1125 C-----------------------------------------------------------------------------
1126 subroutine egbv_long(evdw)
1128 C This subroutine calculates the interaction energy of nonbonded side chains
1129 C assuming the Gay-Berne-Vorobjev potential of interaction.
1131 implicit real*8 (a-h,o-z)
1132 include 'DIMENSIONS'
1133 include 'COMMON.GEO'
1134 include 'COMMON.VAR'
1135 include 'COMMON.LOCAL'
1136 include 'COMMON.CHAIN'
1137 include 'COMMON.DERIV'
1138 include 'COMMON.NAMES'
1139 include 'COMMON.INTERACT'
1140 include 'COMMON.IOUNITS'
1141 include 'COMMON.CALC'
1142 include "COMMON.SPLITELE"
1144 common /srutu/ icall
1146 integer itypi,itypj,itypi1,iint,ind
1147 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1148 & xi,yi,zi,fac_augm,e_augm
1149 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1150 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1151 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1152 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1153 double precision sss1,sssgrad1
1155 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1158 c if (icall.eq.0) lprn=.true.
1160 do i=iatsc_s,iatsc_e
1161 itypi=iabs(itype(i))
1162 if (itypi.eq.ntyp1) cycle
1163 itypi1=iabs(itype(i+1))
1167 dxi=dc_norm(1,nres+i)
1168 dyi=dc_norm(2,nres+i)
1169 dzi=dc_norm(3,nres+i)
1170 c dsci_inv=dsc_inv(itypi)
1171 dsci_inv=vbld_inv(i+nres)
1173 C Calculate SC interaction energy.
1175 do iint=1,nint_gr(i)
1176 do j=istart(i,iint),iend(i,iint)
1178 itypj=iabs(itype(j))
1179 if (itypj.eq.ntyp1) cycle
1180 c dscj_inv=dsc_inv(itypj)
1181 dscj_inv=vbld_inv(j+nres)
1182 sig0ij=sigma(itypi,itypj)
1183 r0ij=r0(itypi,itypj)
1184 chi1=chi(itypi,itypj)
1185 chi2=chi(itypj,itypi)
1192 alf12=0.5D0*(alf1+alf2)
1196 dxj=dc_norm(1,nres+j)
1197 dyj=dc_norm(2,nres+j)
1198 dzj=dc_norm(3,nres+j)
1199 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1202 sss1=sscale(1.0d0/rij,r_cut_int)
1203 if (sss1.eq.0.0d0) cycle
1205 if (sss.lt.1.0d0) then
1207 & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1208 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
1210 C Calculate angle-dependent terms of energy and contributions to their
1214 sig=sig0ij*dsqrt(sigsq)
1215 rij_shift=1.0D0/rij-sig+r0ij
1216 C I hate to put IF's in the loops, but here don't have another choice!!!!
1217 if (rij_shift.le.0.0D0) then
1222 c---------------------------------------------------------------
1223 rij_shift=1.0D0/rij_shift
1224 fac=rij_shift**expon
1227 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1228 eps2der=evdwij*eps3rt
1229 eps3der=evdwij*eps2rt
1230 fac_augm=rrij**expon
1231 e_augm=augm(itypi,itypj)*fac_augm
1232 evdwij=evdwij*eps2rt*eps3rt
1233 evdw=evdw+(evdwij+e_augm)*sss1*(1.0d0-sss)
1235 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1237 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1238 & restyp(itypi),i,restyp(itypj),j,
1239 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1240 & chi1,chi2,chip1,chip2,
1241 & eps1,eps2rt**2,eps3rt**2,
1242 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1245 C Calculate gradient components.
1246 e1=e1*eps1*eps2rt**2*eps3rt**2
1247 fac=-expon*(e1+evdwij)*rij_shift
1249 fac=rij*fac-2*expon*rrij*e_augm
1250 fac=fac+(evdwij+e_augm)*
1251 & (-sss1*sssgrad/(1.0d0-sss)/sigmaii(itypi,itypj)
1252 & +(1.0d0-sss)*sssgrad1/sss1)*rij
1253 C Calculate the radial part of the gradient
1257 C Calculate angular part of the gradient.
1258 call sc_grad_scale((1.0d0-sss)*sss1)
1264 C-----------------------------------------------------------------------------
1265 subroutine egbv_short(evdw)
1267 C This subroutine calculates the interaction energy of nonbonded side chains
1268 C assuming the Gay-Berne-Vorobjev potential of interaction.
1271 include 'DIMENSIONS'
1272 include 'COMMON.GEO'
1273 include 'COMMON.VAR'
1274 include 'COMMON.LOCAL'
1275 include 'COMMON.CHAIN'
1276 include 'COMMON.DERIV'
1277 include 'COMMON.NAMES'
1278 include 'COMMON.INTERACT'
1279 include 'COMMON.IOUNITS'
1280 include 'COMMON.CALC'
1281 include "COMMON.SPLITELE"
1283 common /srutu/ icall
1285 integer itypi,itypj,itypi1,iint,ind
1286 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1287 & xi,yi,zi,fac_augm,e_augm
1288 double precision evdw
1289 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1290 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1291 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1292 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1294 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1297 c if (icall.eq.0) lprn=.true.
1299 do i=iatsc_s,iatsc_e
1300 itypi=iabs(itype(i))
1301 if (itypi.eq.ntyp1) cycle
1302 itypi1=iabs(itype(i+1))
1306 dxi=dc_norm(1,nres+i)
1307 dyi=dc_norm(2,nres+i)
1308 dzi=dc_norm(3,nres+i)
1309 c dsci_inv=dsc_inv(itypi)
1310 dsci_inv=vbld_inv(i+nres)
1312 C Calculate SC interaction energy.
1314 do iint=1,nint_gr(i)
1315 do j=istart(i,iint),iend(i,iint)
1317 itypj=iabs(itype(j))
1318 if (itypj.eq.ntyp1) cycle
1319 c dscj_inv=dsc_inv(itypj)
1320 dscj_inv=vbld_inv(j+nres)
1321 sig0ij=sigma(itypi,itypj)
1322 r0ij=r0(itypi,itypj)
1323 chi1=chi(itypi,itypj)
1324 chi2=chi(itypj,itypi)
1331 alf12=0.5D0*(alf1+alf2)
1335 dxj=dc_norm(1,nres+j)
1336 dyj=dc_norm(2,nres+j)
1337 dzj=dc_norm(3,nres+j)
1338 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1341 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1343 if (sss.gt.0.0d0) then
1345 C Calculate angle-dependent terms of energy and contributions to their
1349 sig=sig0ij*dsqrt(sigsq)
1350 rij_shift=1.0D0/rij-sig+r0ij
1351 C I hate to put IF's in the loops, but here don't have another choice!!!!
1352 if (rij_shift.le.0.0D0) then
1357 c---------------------------------------------------------------
1358 rij_shift=1.0D0/rij_shift
1359 fac=rij_shift**expon
1362 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1363 eps2der=evdwij*eps3rt
1364 eps3der=evdwij*eps2rt
1365 fac_augm=rrij**expon
1366 e_augm=augm(itypi,itypj)*fac_augm
1367 evdwij=evdwij*eps2rt*eps3rt
1368 evdw=evdw+(evdwij+e_augm)*sss
1370 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1372 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1373 & restyp(itypi),i,restyp(itypj),j,
1374 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1375 & chi1,chi2,chip1,chip2,
1376 & eps1,eps2rt**2,eps3rt**2,
1377 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1380 C Calculate gradient components.
1381 e1=e1*eps1*eps2rt**2*eps3rt**2
1382 fac=-expon*(e1+evdwij)*rij_shift
1384 fac=rij*fac-2*expon*rrij*e_augm+
1385 & (evdwij+e_augm)*sssgrad/sigmaii(itypi,itypj)/sss*rij
1386 C Calculate the radial part of the gradient
1390 C Calculate angular part of the gradient.
1391 call sc_grad_scale(sss)
1397 C----------------------------------------------------------------------------
1398 subroutine sc_grad_scale(scalfac)
1399 implicit real*8 (a-h,o-z)
1400 include 'DIMENSIONS'
1401 include 'COMMON.CHAIN'
1402 include 'COMMON.DERIV'
1403 include 'COMMON.CALC'
1404 include 'COMMON.IOUNITS'
1405 include "COMMON.SPLITELE"
1406 double precision dcosom1(3),dcosom2(3)
1407 double precision scalfac
1408 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1409 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1410 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1411 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1415 c eom12=evdwij*eps1_om12
1417 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1418 c & " sigder",sigder
1419 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1420 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1422 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1423 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1426 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1428 c write (iout,*) "gg",(gg(k),k=1,3)
1430 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
1431 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1432 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1433 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
1434 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1435 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1436 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1437 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1438 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1439 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1442 C Calculate the components of the gradient in DC and X
1445 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
1446 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
1450 C--------------------------------------------------------------------------
1451 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1453 C This subroutine calculates the average interaction energy and its gradient
1454 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1455 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1456 C The potential depends both on the distance of peptide-group centers and on
1457 C the orientation of the CA-CA virtual bonds.
1459 implicit real*8 (a-h,o-z)
1463 include 'DIMENSIONS'
1464 include 'COMMON.CONTROL'
1465 include 'COMMON.SETUP'
1466 include 'COMMON.IOUNITS'
1467 include 'COMMON.GEO'
1468 include 'COMMON.VAR'
1469 include 'COMMON.LOCAL'
1470 include 'COMMON.CHAIN'
1471 include 'COMMON.DERIV'
1472 include 'COMMON.INTERACT'
1474 include 'COMMON.CONTACTS'
1475 include 'COMMON.CONTMAT'
1477 include 'COMMON.CORRMAT'
1478 include 'COMMON.TORSION'
1479 include 'COMMON.VECTORS'
1480 include 'COMMON.FFIELD'
1481 include 'COMMON.TIME1'
1482 include 'COMMON.SHIELD'
1483 include "COMMON.SPLITELE"
1484 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1485 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1486 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1487 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1488 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1489 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1491 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1493 double precision scal_el /1.0d0/
1495 double precision scal_el /0.5d0/
1498 C 13-go grudnia roku pamietnego...
1499 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1500 & 0.0d0,1.0d0,0.0d0,
1501 & 0.0d0,0.0d0,1.0d0/
1502 cd write(iout,*) 'In EELEC'
1504 cd write(iout,*) 'Type',i
1505 cd write(iout,*) 'B1',B1(:,i)
1506 cd write(iout,*) 'B2',B2(:,i)
1507 cd write(iout,*) 'CC',CC(:,:,i)
1508 cd write(iout,*) 'DD',DD(:,:,i)
1509 cd write(iout,*) 'EE',EE(:,:,i)
1511 cd call check_vecgrad
1513 C print *,"WCHODZE3"
1514 if (icheckgrad.eq.1) then
1516 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1518 dc_norm(k,i)=dc(k,i)*fac
1520 c write (iout,*) 'i',i,' fac',fac
1523 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1524 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1525 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1526 c call vec_and_deriv
1532 time_mat=time_mat+MPI_Wtime()-time01
1536 cd write (iout,*) 'i=',i
1538 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1541 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1542 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1557 cd print '(a)','Enter EELEC'
1558 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1560 gel_loc_loc(i)=0.0d0
1565 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1567 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1569 do i=iturn3_start,iturn3_end
1570 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1571 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1572 C & .or. itype(i-1).eq.ntyp1
1573 C & .or. itype(i+4).eq.ntyp1
1578 dx_normi=dc_norm(1,i)
1579 dy_normi=dc_norm(2,i)
1580 dz_normi=dc_norm(3,i)
1581 xmedi=c(1,i)+0.5d0*dxi
1582 ymedi=c(2,i)+0.5d0*dyi
1583 zmedi=c(3,i)+0.5d0*dzi
1584 xmedi=mod(xmedi,boxxsize)
1585 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1586 ymedi=mod(ymedi,boxysize)
1587 if (ymedi.lt.0) ymedi=ymedi+boxysize
1588 zmedi=mod(zmedi,boxzsize)
1589 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1591 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1592 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1594 num_cont_hb(i)=num_conti
1597 do i=iturn4_start,iturn4_end
1598 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1599 & .or. itype(i+3).eq.ntyp1
1600 & .or. itype(i+4).eq.ntyp1
1601 C & .or. itype(i+5).eq.ntyp1
1602 C & .or. itype(i-1).eq.ntyp1
1607 dx_normi=dc_norm(1,i)
1608 dy_normi=dc_norm(2,i)
1609 dz_normi=dc_norm(3,i)
1610 xmedi=c(1,i)+0.5d0*dxi
1611 ymedi=c(2,i)+0.5d0*dyi
1612 zmedi=c(3,i)+0.5d0*dzi
1613 xmedi=mod(xmedi,boxxsize)
1614 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1615 ymedi=mod(ymedi,boxysize)
1616 if (ymedi.lt.0) ymedi=ymedi+boxysize
1617 zmedi=mod(zmedi,boxzsize)
1618 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1620 num_conti=num_cont_hb(i)
1622 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1623 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1624 & call eturn4(i,eello_turn4)
1626 num_cont_hb(i)=num_conti
1630 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1632 do i=iatel_s,iatel_e
1633 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1634 C & .or. itype(i+2).eq.ntyp1
1635 C & .or. itype(i-1).eq.ntyp1
1640 dx_normi=dc_norm(1,i)
1641 dy_normi=dc_norm(2,i)
1642 dz_normi=dc_norm(3,i)
1643 xmedi=c(1,i)+0.5d0*dxi
1644 ymedi=c(2,i)+0.5d0*dyi
1645 zmedi=c(3,i)+0.5d0*dzi
1646 xmedi=mod(xmedi,boxxsize)
1647 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1648 ymedi=mod(ymedi,boxysize)
1649 if (ymedi.lt.0) ymedi=ymedi+boxysize
1650 zmedi=mod(zmedi,boxzsize)
1651 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1652 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
1654 num_conti=num_cont_hb(i)
1656 do j=ielstart(i),ielend(i)
1657 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1658 C & .or.itype(j+2).eq.ntyp1
1659 C & .or.itype(j-1).eq.ntyp1
1661 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1664 num_cont_hb(i)=num_conti
1667 c write (iout,*) "Number of loop steps in EELEC:",ind
1669 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1670 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1672 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1673 ccc eel_loc=eel_loc+eello_turn3
1674 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1677 C-------------------------------------------------------------------------------
1678 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1680 include 'DIMENSIONS'
1684 include 'COMMON.CONTROL'
1685 include 'COMMON.IOUNITS'
1686 include 'COMMON.GEO'
1687 include 'COMMON.VAR'
1688 include 'COMMON.LOCAL'
1689 include 'COMMON.CHAIN'
1690 include 'COMMON.DERIV'
1691 include 'COMMON.INTERACT'
1693 include 'COMMON.CONTACTS'
1694 include 'COMMON.CONTMAT'
1696 include 'COMMON.CORRMAT'
1697 include 'COMMON.TORSION'
1698 include 'COMMON.VECTORS'
1699 include 'COMMON.FFIELD'
1700 include 'COMMON.TIME1'
1701 include 'COMMON.SHIELD'
1702 include "COMMON.SPLITELE"
1703 integer xshift,yshift,zshift
1704 double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1705 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1706 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1707 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
1708 & gmuij2(4),gmuji2(4)
1709 integer j1,j2,num_conti
1710 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1711 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1713 integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ind,itypi,itypj
1714 integer ilist,iresshield
1715 double precision rlocshield
1716 double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
1717 double precision ees,evdw1,eel_loc,aaa,bbb,ael3i
1718 double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,
1719 & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,
1720 & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa,
1721 & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,
1722 & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,
1723 & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,
1724 & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
1725 & ecosgp,ecosam,ecosbm,ecosgm,ghalf,geel_loc_ij,geel_loc_ji,
1726 & dxi,dyi,dzi,a22,a23,a32,a33
1727 double precision dist_init,xmedi,ymedi,zmedi,xj_safe,yj_safe,
1728 & zj_safe,xj_temp,yj_temp,zj_temp,dist_temp,dx_normi,dy_normi,
1730 double precision sss1,sssgrad1
1731 double precision sscale,sscagrad
1732 double precision scalar
1734 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1736 double precision scal_el /1.0d0/
1738 double precision scal_el /0.5d0/
1741 C 13-go grudnia roku pamietnego...
1742 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1743 & 0.0d0,1.0d0,0.0d0,
1744 & 0.0d0,0.0d0,1.0d0/
1745 c time00=MPI_Wtime()
1746 cd write (iout,*) "eelecij",i,j
1747 C print *,"WCHODZE2"
1751 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1752 aaa=app(iteli,itelj)
1753 bbb=bpp(iteli,itelj)
1754 ael6i=ael6(iteli,itelj)
1755 ael3i=ael3(iteli,itelj)
1759 dx_normj=dc_norm(1,j)
1760 dy_normj=dc_norm(2,j)
1761 dz_normj=dc_norm(3,j)
1766 if (xj.lt.0) xj=xj+boxxsize
1768 if (yj.lt.0) yj=yj+boxysize
1770 if (zj.lt.0) zj=zj+boxzsize
1771 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1779 xj=xj_safe+xshift*boxxsize
1780 yj=yj_safe+yshift*boxysize
1781 zj=zj_safe+zshift*boxzsize
1782 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1783 if(dist_temp.lt.dist_init) then
1793 if (isubchap.eq.1) then
1803 rij=xj*xj+yj*yj+zj*zj
1807 c For extracting the short-range part of Evdwpp
1808 sss1=sscale(rij,r_cut_int)
1809 if (sss1.eq.0.0d0) return
1810 sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
1811 sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
1812 sssgrad1=sscagrad(rij,r_cut_int)
1815 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1816 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1817 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1818 fac=cosa-3.0D0*cosb*cosg
1820 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1821 if (j.eq.i+2) ev1=scal_el*ev1
1826 if (shield_mode.eq.0) then
1830 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1832 el1=el1*fac_shield(i)**2*fac_shield(j)**2
1833 el2=el2*fac_shield(i)**2*fac_shield(j)**2
1835 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1836 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1838 evdw1=evdw1+evdwij*(1.0d0-sss)*sss1
1839 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1840 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1841 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1842 cd & xmedi,ymedi,zmedi,xj,yj,zj
1844 if (energy_dec) then
1845 write (iout,'(a6,2i5,0pf7.3,2f7.3)')
1846 & 'evdw1',i,j,evdwij,sss,sss1
1847 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1851 C Calculate contributions to the Cartesian gradient.
1854 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
1855 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1856 facel=-3*rrmij*(el1+eesij)*sss1
1857 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1863 * Radial derivatives. First process both termini of the fragment (i,j)
1865 aux=facel+sssgrad1*(1.0d0-sss)*eesij*rmij
1866 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1873 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1874 & (shield_mode.gt.0)) then
1876 do ilist=1,ishield_list(i)
1877 iresshield=shield_list(ilist,i)
1879 rlocshield=grad_shield_side(k,ilist,i)*eesij*sss1
1880 & /fac_shield(i)*2.0*sss1
1881 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1883 & +grad_shield_loc(k,ilist,i)*eesij*sss1/fac_shield(i)*2.0
1885 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1886 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1887 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1888 C if (iresshield.gt.i) then
1889 C do ishi=i+1,iresshield-1
1890 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1891 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1895 C do ishi=iresshield,i
1896 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1897 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1903 do ilist=1,ishield_list(j)
1904 iresshield=shield_list(ilist,j)
1906 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1908 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1910 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss1
1911 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1916 gshieldc(k,i)=gshieldc(k,i)+
1917 & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
1918 gshieldc(k,j)=gshieldc(k,j)+
1919 & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
1920 gshieldc(k,i-1)=gshieldc(k,i-1)+
1921 & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
1922 gshieldc(k,j-1)=gshieldc(k,j-1)+
1923 & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
1929 c ghalf=0.5D0*ggg(k)
1930 c gelc(k,i)=gelc(k,i)+ghalf
1931 c gelc(k,j)=gelc(k,j)+ghalf
1933 c 9/28/08 AL Gradient compotents will be summed only at the end
1935 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1936 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1938 c gelc_long(3,i)=gelc_long(3,i)+
1939 c ssgradlipi*eesij/2.0d0*lipscale**2*sss1
1941 * Loop over residues i+1 thru j-1.
1945 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1949 & (-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij
1950 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1955 c ghalf=0.5D0*ggg(k)
1956 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1957 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1959 c 9/28/08 AL Gradient compotents will be summed only at the end
1961 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1962 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1965 * Loop over residues i+1 thru j-1.
1969 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
1973 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
1974 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1975 facel=-3*rrmij*(el1+eesij)*sss1
1976 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1978 c facvdw=ev1+evdwij*(1.0d0-sss)*sss1
1981 fac=-3*rrmij*(facvdw+facvdw+facel)
1986 * Radial derivatives. First process both termini of the fragment (i,j)
1988 aux=fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj))
1990 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1998 c ghalf=0.5D0*ggg(k)
1999 c gelc(k,i)=gelc(k,i)+ghalf
2000 c gelc(k,j)=gelc(k,j)+ghalf
2002 c 9/28/08 AL Gradient compotents will be summed only at the end
2004 gelc_long(k,j)=gelc(k,j)+ggg(k)
2005 gelc_long(k,i)=gelc(k,i)-ggg(k)
2008 * Loop over residues i+1 thru j-1.
2012 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
2015 c 9/28/08 AL Gradient compotents will be summed only at the end
2020 & (-sssgrad*sss1/rpp(iteli,itelj)+sssgrad1*(1.0d0-sss))*rmij*evdwij
2025 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2026 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2032 ecosa=2.0D0*fac3*fac1+fac4
2035 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
2036 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
2038 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2039 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2041 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
2042 cd & (dcosg(k),k=1,3)
2044 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1
2045 & *fac_shield(i)**2*fac_shield(j)**2
2046 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2050 c ghalf=0.5D0*ggg(k)
2051 c gelc(k,i)=gelc(k,i)+ghalf
2052 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
2053 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2054 c gelc(k,j)=gelc(k,j)+ghalf
2055 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2056 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2060 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
2065 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
2066 & +ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss1
2067 & *fac_shield(i)**2*fac_shield(j)**2
2068 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2071 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2072 & +ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss1
2073 & *fac_shield(i)**2*fac_shield(j)**2
2074 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2075 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
2076 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
2078 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
2079 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
2080 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2082 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
2083 C energy of a peptide unit is assumed in the form of a second-order
2084 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
2085 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
2086 C are computed for EVERY pair of non-contiguous peptide groups.
2088 if (j.lt.nres-1) then
2099 muij(kkk)=mu(k,i)*mu(l,j)
2101 gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
2102 c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
2103 gmuij2(kkk)=gUb2(k,i)*mu(l,j)
2104 gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
2105 c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
2106 gmuji2(kkk)=mu(k,i)*gUb2(l,j)
2110 cd write (iout,*) 'EELEC: i',i,' j',j
2111 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
2112 cd write(iout,*) 'muij',muij
2113 ury=scalar(uy(1,i),erij)
2114 urz=scalar(uz(1,i),erij)
2115 vry=scalar(uy(1,j),erij)
2116 vrz=scalar(uz(1,j),erij)
2117 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
2118 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
2119 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
2120 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
2121 fac=dsqrt(-ael6i)*r3ij
2126 cd write (iout,'(4i5,4f10.5)')
2127 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
2128 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
2129 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
2130 cd & uy(:,j),uz(:,j)
2131 cd write (iout,'(4f10.5)')
2132 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
2133 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
2134 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
2135 cd write (iout,'(9f10.5/)')
2136 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
2137 C Derivatives of the elements of A in virtual-bond vectors
2138 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
2140 uryg(k,1)=scalar(erder(1,k),uy(1,i))
2141 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
2142 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
2143 urzg(k,1)=scalar(erder(1,k),uz(1,i))
2144 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
2145 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
2146 vryg(k,1)=scalar(erder(1,k),uy(1,j))
2147 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
2148 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
2149 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
2150 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
2151 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
2153 C Compute radial contributions to the gradient
2171 C Add the contributions coming from er
2174 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
2175 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
2176 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
2177 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
2180 C Derivatives in DC(i)
2181 cgrad ghalf1=0.5d0*agg(k,1)
2182 cgrad ghalf2=0.5d0*agg(k,2)
2183 cgrad ghalf3=0.5d0*agg(k,3)
2184 cgrad ghalf4=0.5d0*agg(k,4)
2185 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2186 & -3.0d0*uryg(k,2)*vry)!+ghalf1
2187 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2188 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
2189 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2190 & -3.0d0*urzg(k,2)*vry)!+ghalf3
2191 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2192 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
2193 C Derivatives in DC(i+1)
2194 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2195 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2196 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2197 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2198 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2199 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2200 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2201 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2202 C Derivatives in DC(j)
2203 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2204 & -3.0d0*vryg(k,2)*ury)!+ghalf1
2205 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2206 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
2207 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2208 & -3.0d0*vryg(k,2)*urz)!+ghalf3
2209 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2210 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
2211 C Derivatives in DC(j+1) or DC(nres-1)
2212 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2213 & -3.0d0*vryg(k,3)*ury)
2214 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2215 & -3.0d0*vrzg(k,3)*ury)
2216 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2217 & -3.0d0*vryg(k,3)*urz)
2218 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2219 & -3.0d0*vrzg(k,3)*urz)
2220 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
2222 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
2235 aggi(k,l)=-aggi(k,l)
2236 aggi1(k,l)=-aggi1(k,l)
2237 aggj(k,l)=-aggj(k,l)
2238 aggj1(k,l)=-aggj1(k,l)
2241 if (j.lt.nres-1) then
2247 aggi(k,l)=-aggi(k,l)
2248 aggi1(k,l)=-aggi1(k,l)
2249 aggj(k,l)=-aggj(k,l)
2250 aggj1(k,l)=-aggj1(k,l)
2261 aggi(k,l)=-aggi(k,l)
2262 aggi1(k,l)=-aggi1(k,l)
2263 aggj(k,l)=-aggj(k,l)
2264 aggj1(k,l)=-aggj1(k,l)
2269 IF (wel_loc.gt.0.0d0) THEN
2270 C Contribution to the local-electrostatic energy coming from the i-j pair
2271 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2273 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2275 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2276 & 'eelloc',i,j,eel_loc_ij
2279 if (shield_mode.eq.0) then
2286 eel_loc_ij=eel_loc_ij
2287 & *fac_shield(i)*fac_shield(j)*sss1
2288 eel_loc=eel_loc+eel_loc_ij
2290 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2291 & (shield_mode.gt.0)) then
2294 do ilist=1,ishield_list(i)
2295 iresshield=shield_list(ilist,i)
2297 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2300 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2302 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2303 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2307 do ilist=1,ishield_list(j)
2308 iresshield=shield_list(ilist,j)
2310 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2313 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2315 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2316 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2323 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2324 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2325 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2326 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2327 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2328 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2329 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2330 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2335 geel_loc_ij=(a22*gmuij1(1)
2339 & *fac_shield(i)*fac_shield(j)*sss1
2340 c write(iout,*) "derivative over thatai"
2341 c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
2343 gloc(nphi+i,icg)=gloc(nphi+i,icg)+
2344 & geel_loc_ij*wel_loc
2345 c write(iout,*) "derivative over thatai-1"
2346 c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
2353 gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
2354 & geel_loc_ij*wel_loc
2355 & *fac_shield(i)*fac_shield(j)*sss1
2357 c Derivative over j residue
2358 geel_loc_ji=a22*gmuji1(1)
2362 c write(iout,*) "derivative over thataj"
2363 c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
2366 gloc(nphi+j,icg)=gloc(nphi+j,icg)+
2367 & geel_loc_ji*wel_loc
2368 & *fac_shield(i)*fac_shield(j)*sss1
2375 c write(iout,*) "derivative over thataj-1"
2376 c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
2378 gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
2379 & geel_loc_ji*wel_loc
2380 & *fac_shield(i)*fac_shield(j)*sss1
2382 cC Paral derivatives in virtual-bond dihedral angles gamma
2384 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2385 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2386 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2387 & *fac_shield(i)*fac_shield(j)*sss1
2388 c & *fac_shield(i)*fac_shield(j)
2389 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2392 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2393 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2394 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2395 & *fac_shield(i)*fac_shield(j)*sss1
2396 c & *fac_shield(i)*fac_shield(j)
2397 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2399 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2400 aux=eel_loc_ij/sss1*sssgrad1*rmij
2405 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+
2406 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2407 & *fac_shield(i)*fac_shield(j)*sss1
2408 c & *fac_shield(i)*fac_shield(j)
2409 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2411 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2412 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2413 cgrad ghalf=0.5d0*ggg(l)
2414 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
2415 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
2419 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2422 C Remaining derivatives of eello
2423 c gel_loc_long(3,j)=gel_loc_long(3,j)+ &
2424 c ssgradlipj*eel_loc_ij/2.0d0*lipscale/ &
2425 c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2427 c gel_loc_long(3,i)=gel_loc_long(3,i)+ &
2428 c ssgradlipi*eel_loc_ij/2.0d0*lipscale/ &
2429 c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2432 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2433 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2434 & *fac_shield(i)*fac_shield(j)*sss1
2435 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2437 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2438 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2439 & *fac_shield(i)*fac_shield(j)*sss1
2440 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2442 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2443 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2444 & *fac_shield(i)*fac_shield(j)*sss1
2445 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2447 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2448 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2449 & *fac_shield(i)*fac_shield(j)*sss1
2450 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2455 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2456 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
2457 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2458 & .and. num_conti.le.maxconts) then
2459 c write (iout,*) i,j," entered corr"
2461 C Calculate the contact function. The ith column of the array JCONT will
2462 C contain the numbers of atoms that make contacts with the atom I (of numbers
2463 C greater than I). The arrays FACONT and GACONT will contain the values of
2464 C the contact function and its derivative.
2465 c r0ij=1.02D0*rpp(iteli,itelj)
2466 c r0ij=1.11D0*rpp(iteli,itelj)
2467 r0ij=2.20D0*rpp(iteli,itelj)
2468 c r0ij=1.55D0*rpp(iteli,itelj)
2469 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2470 if (fcont.gt.0.0D0) then
2471 num_conti=num_conti+1
2472 if (num_conti.gt.maxconts) then
2473 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2474 & ' will skip next contacts for this conf.'
2476 jcont_hb(num_conti,i)=j
2477 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
2478 cd " jcont_hb",jcont_hb(num_conti,i)
2479 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2480 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2481 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2483 d_cont(num_conti,i)=rij
2484 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2485 C --- Electrostatic-interaction matrix ---
2486 a_chuj(1,1,num_conti,i)=a22
2487 a_chuj(1,2,num_conti,i)=a23
2488 a_chuj(2,1,num_conti,i)=a32
2489 a_chuj(2,2,num_conti,i)=a33
2490 C --- Gradient of rij
2492 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2499 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2500 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2501 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2502 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2503 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2508 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2509 C Calculate contact energies
2511 wij=cosa-3.0D0*cosb*cosg
2514 c fac3=dsqrt(-ael6i)/r0ij**3
2515 fac3=dsqrt(-ael6i)*r3ij
2516 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2517 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2518 if (ees0tmp.gt.0) then
2519 ees0pij=dsqrt(ees0tmp)
2523 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2524 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2525 if (ees0tmp.gt.0) then
2526 ees0mij=dsqrt(ees0tmp)
2531 if (shield_mode.eq.0) then
2535 ees0plist(num_conti,i)=j
2536 C fac_shield(i)=0.4d0
2537 C fac_shield(j)=0.6d0
2539 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2540 & *fac_shield(i)*fac_shield(j)*sss1
2541 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2542 & *fac_shield(i)*fac_shield(j)*sss1
2544 C Diagnostics. Comment out or remove after debugging!
2545 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2546 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2547 c ees0m(num_conti,i)=0.0D0
2549 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2550 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2551 C Angular derivatives of the contact function
2552 ees0pij1=fac3/ees0pij
2553 ees0mij1=fac3/ees0mij
2554 fac3p=-3.0D0*fac3*rrmij
2555 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2556 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2558 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2559 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2560 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2561 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2562 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2563 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2564 ecosap=ecosa1+ecosa2
2565 ecosbp=ecosb1+ecosb2
2566 ecosgp=ecosg1+ecosg2
2567 ecosam=ecosa1-ecosa2
2568 ecosbm=ecosb1-ecosb2
2569 ecosgm=ecosg1-ecosg2
2578 facont_hb(num_conti,i)=fcont
2579 fprimcont=fprimcont/rij
2580 cd facont_hb(num_conti,i)=1.0D0
2581 C Following line is for diagnostics.
2584 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2585 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2588 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2589 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2591 gggp(1)=gggp(1)+ees0pijp*xj
2592 & +ees0p(num_conti,i)/sss1*rmij*xj*sssgrad1
2593 gggp(2)=gggp(2)+ees0pijp*yj
2594 & +ees0p(num_conti,i)/sss1*rmij*yj*sssgrad1
2595 gggp(3)=gggp(3)+ees0pijp*zj
2596 & +ees0p(num_conti,i)/sss1*rmij*zj*sssgrad1
2597 gggm(1)=gggm(1)+ees0mijp*xj
2598 & +ees0m(num_conti,i)/sss1*rmij*xj*sssgrad1
2599 gggm(2)=gggm(2)+ees0mijp*yj
2600 & +ees0m(num_conti,i)/sss1*rmij*yj*sssgrad1
2601 gggm(3)=gggm(3)+ees0mijp*zj
2602 & +ees0m(num_conti,i)/sss1*rmij*zj*sssgrad1
2603 C Derivatives due to the contact function
2604 gacont_hbr(1,num_conti,i)=fprimcont*xj
2605 gacont_hbr(2,num_conti,i)=fprimcont*yj
2606 gacont_hbr(3,num_conti,i)=fprimcont*zj
2609 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2610 c following the change of gradient-summation algorithm.
2612 cgrad ghalfp=0.5D0*gggp(k)
2613 cgrad ghalfm=0.5D0*gggm(k)
2614 gacontp_hb1(k,num_conti,i)=!ghalfp
2615 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2616 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2617 & *sss1*fac_shield(i)*fac_shield(j)
2619 gacontp_hb2(k,num_conti,i)=!ghalfp
2620 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2621 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2622 & *sss1*fac_shield(i)*fac_shield(j)
2624 gacontp_hb3(k,num_conti,i)=gggp(k)
2625 & *sss1*fac_shield(i)*fac_shield(j)
2627 gacontm_hb1(k,num_conti,i)=!ghalfm
2628 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2629 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2630 & *sss1*fac_shield(i)*fac_shield(j)
2632 gacontm_hb2(k,num_conti,i)=!ghalfm
2633 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2634 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2635 & *sss1*fac_shield(i)*fac_shield(j)
2637 gacontm_hb3(k,num_conti,i)=gggm(k)
2638 & *sss1*fac_shield(i)*fac_shield(j)
2642 endif ! num_conti.le.maxconts
2646 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2649 ghalf=0.5d0*agg(l,k)
2650 aggi(l,k)=aggi(l,k)+ghalf
2651 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2652 aggj(l,k)=aggj(l,k)+ghalf
2655 if (j.eq.nres-1 .and. i.lt.j-2) then
2658 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2663 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2666 C-----------------------------------------------------------------------
2667 subroutine evdwpp_short(evdw1)
2672 include 'DIMENSIONS'
2673 include 'COMMON.CONTROL'
2674 include 'COMMON.IOUNITS'
2675 include 'COMMON.GEO'
2676 include 'COMMON.VAR'
2677 include 'COMMON.LOCAL'
2678 include 'COMMON.CHAIN'
2679 include 'COMMON.DERIV'
2680 include 'COMMON.INTERACT'
2681 c include 'COMMON.CONTACTS'
2682 include 'COMMON.TORSION'
2683 include 'COMMON.VECTORS'
2684 include 'COMMON.FFIELD'
2685 include "COMMON.SPLITELE"
2686 double precision ggg(3)
2687 integer xshift,yshift,zshift
2688 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2690 double precision scal_el /1.0d0/
2692 double precision scal_el /0.5d0/
2694 c write (iout,*) "evdwpp_short"
2695 integer i,j,k,iteli,itelj,num_conti,ind,isubchap
2696 double precision dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb
2697 double precision xj,yj,zj,rij,rrmij,r3ij,r6ij,evdw1,
2698 & dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
2699 & dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw
2700 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
2701 & dist_temp, dist_init,sss_grad
2702 double precision sscale,sscagrad
2705 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2706 c & " iatel_e_vdw",iatel_e_vdw
2708 do i=iatel_s_vdw,iatel_e_vdw
2709 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2713 dx_normi=dc_norm(1,i)
2714 dy_normi=dc_norm(2,i)
2715 dz_normi=dc_norm(3,i)
2716 xmedi=c(1,i)+0.5d0*dxi
2717 ymedi=c(2,i)+0.5d0*dyi
2718 zmedi=c(3,i)+0.5d0*dzi
2719 xmedi=mod(xmedi,boxxsize)
2720 if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
2721 ymedi=mod(ymedi,boxysize)
2722 if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
2723 zmedi=mod(zmedi,boxzsize)
2724 if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
2726 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2727 c & ' ielend',ielend_vdw(i)
2729 do j=ielstart_vdw(i),ielend_vdw(i)
2730 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2734 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2735 aaa=app(iteli,itelj)
2736 bbb=bpp(iteli,itelj)
2740 dx_normj=dc_norm(1,j)
2741 dy_normj=dc_norm(2,j)
2742 dz_normj=dc_norm(3,j)
2747 if (xj.lt.0) xj=xj+boxxsize
2749 if (yj.lt.0) yj=yj+boxysize
2751 if (zj.lt.0) zj=zj+boxzsize
2752 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2760 xj=xj_safe+xshift*boxxsize
2761 yj=yj_safe+yshift*boxysize
2762 zj=zj_safe+zshift*boxzsize
2763 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2764 if(dist_temp.lt.dist_init) then
2774 if (isubchap.eq.1) then
2783 rij=xj*xj+yj*yj+zj*zj
2786 c sss=sscale(rij/rpp(iteli,itelj))
2787 c sssgrad=sscagrad(rij/rpp(iteli,itelj))
2788 sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
2789 sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
2790 if (sss.gt.0.0d0) then
2795 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2796 if (j.eq.i+2) ev1=scal_el*ev1
2799 if (energy_dec) then
2800 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2802 evdw1=evdw1+evdwij*sss
2803 if (energy_dec) write (iout,'(a10,2i5,0pf7.3)')
2804 & 'evdw1_sum',i,j,evdw1
2806 C Calculate contributions to the Cartesian gradient.
2808 facvdw=-6*rrmij*(ev1+evdwij)*sss
2809 ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2810 ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2811 ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2816 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2817 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2824 C-----------------------------------------------------------------------------
2825 subroutine escp_long(evdw2,evdw2_14)
2827 C This subroutine calculates the excluded-volume interaction energy between
2828 C peptide-group centers and side chains and its gradient in virtual-bond and
2829 C side-chain vectors.
2831 implicit real*8 (a-h,o-z)
2832 include 'DIMENSIONS'
2833 include 'COMMON.GEO'
2834 include 'COMMON.VAR'
2835 include 'COMMON.LOCAL'
2836 include 'COMMON.CHAIN'
2837 include 'COMMON.DERIV'
2838 include 'COMMON.INTERACT'
2839 include 'COMMON.FFIELD'
2840 include 'COMMON.IOUNITS'
2841 include 'COMMON.CONTROL'
2842 include "COMMON.SPLITELE"
2843 logical lprint_short
2844 common /shortcheck/ lprint_short
2845 double precision ggg(3)
2846 integer xshift,yshift,zshift
2847 integer i,iint,j,k,iteli,itypj,subchap
2848 double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
2850 double precision evdw2,evdw2_14,evdwij
2851 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
2852 & dist_temp, dist_init
2853 double precision sscale,sscagrad
2854 if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
2857 CD print '(a)','Enter ESCP KURWA'
2858 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2860 c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
2861 c & ' iatscp_e=',iatscp_e
2862 do i=iatscp_s,iatscp_e
2863 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2865 xi=0.5D0*(c(1,i)+c(1,i+1))
2866 yi=0.5D0*(c(2,i)+c(2,i+1))
2867 zi=0.5D0*(c(3,i)+c(3,i+1))
2869 if (xi.lt.0) xi=xi+boxxsize
2871 if (yi.lt.0) yi=yi+boxysize
2873 if (zi.lt.0) zi=zi+boxzsize
2875 do iint=1,nscp_gr(i)
2877 do j=iscpstart(i,iint),iscpend(i,iint)
2878 itypj=iabs(itype(j))
2879 if (itypj.eq.ntyp1) cycle
2880 C Uncomment following three lines for SC-p interactions
2884 C Uncomment following three lines for Ca-p interactions
2890 if (xj.lt.0) xj=xj+boxxsize
2892 if (yj.lt.0) yj=yj+boxysize
2894 if (zj.lt.0) zj=zj+boxzsize
2896 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2904 xj=xj_safe+xshift*boxxsize
2905 yj=yj_safe+yshift*boxysize
2906 zj=zj_safe+zshift*boxzsize
2907 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2908 if(dist_temp.lt.dist_init) then
2918 if (subchap.eq.1) then
2928 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2930 sss1=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
2931 if (sss1.eq.0) cycle
2932 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
2934 & sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
2935 sssgrad1=sscagrad(1.0d0/dsqrt(rrij),r_cut_int)
2936 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2937 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2938 if (sss.lt.1.0d0) then
2940 e1=fac*fac*aad(itypj,iteli)
2941 e2=fac*bad(itypj,iteli)
2942 if (iabs(j-i) .le. 2) then
2945 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)*sss1
2948 evdw2=evdw2+evdwij*(1.0d0-sss)*sss1
2949 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2950 & 'evdw2',i,j,sss,evdwij
2952 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2955 fac=-(evdwij+e1)*rrij*(1.0d0-sss)*sss1
2956 fac=fac+evdwij*dsqrt(rrij)*(-sssgrad/rscp(itypj,iteli)
2961 C Uncomment following three lines for SC-p interactions
2963 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2965 C Uncomment following line for SC-p interactions
2966 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2968 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
2969 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
2978 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2979 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
2980 gradx_scp(j,i)=expon*gradx_scp(j,i)
2983 C******************************************************************************
2987 C To save time the factor EXPON has been extracted from ALL components
2988 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2991 C******************************************************************************
2994 C-----------------------------------------------------------------------------
2995 subroutine escp_short(evdw2,evdw2_14)
2997 C This subroutine calculates the excluded-volume interaction energy between
2998 C peptide-group centers and side chains and its gradient in virtual-bond and
2999 C side-chain vectors.
3002 include 'DIMENSIONS'
3003 include 'COMMON.GEO'
3004 include 'COMMON.VAR'
3005 include 'COMMON.LOCAL'
3006 include 'COMMON.CHAIN'
3007 include 'COMMON.DERIV'
3008 include 'COMMON.INTERACT'
3009 include 'COMMON.FFIELD'
3010 include 'COMMON.IOUNITS'
3011 include 'COMMON.CONTROL'
3012 include "COMMON.SPLITELE"
3013 integer xshift,yshift,zshift
3014 logical lprint_short
3015 common /shortcheck/ lprint_short
3016 integer i,iint,j,k,iteli,itypj,subchap
3017 double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
3019 double precision evdw2,evdw2_14,evdwij
3020 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
3021 & dist_temp, dist_init
3022 double precision ggg(3)
3023 double precision sscale,sscagrad
3026 cd print '(a)','Enter ESCP'
3028 c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
3029 c & ' iatscp_e=',iatscp_e
3030 if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
3031 do i=iatscp_s,iatscp_e
3032 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
3034 xi=0.5D0*(c(1,i)+c(1,i+1))
3035 yi=0.5D0*(c(2,i)+c(2,i+1))
3036 zi=0.5D0*(c(3,i)+c(3,i+1))
3038 if (xi.lt.0) xi=xi+boxxsize
3040 if (yi.lt.0) yi=yi+boxysize
3042 if (zi.lt.0) zi=zi+boxzsize
3045 c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
3046 c & " nscp_gr",nscp_gr(i)
3047 do iint=1,nscp_gr(i)
3049 do j=iscpstart(i,iint),iscpend(i,iint)
3050 itypj=iabs(itype(j))
3052 c & write (iout,*) "j",j," itypj",itypj
3053 if (itypj.eq.ntyp1) cycle
3054 C Uncomment following three lines for SC-p interactions
3058 C Uncomment following three lines for Ca-p interactions
3064 if (xj.lt.0) xj=xj+boxxsize
3066 if (yj.lt.0) yj=yj+boxysize
3068 if (zj.lt.0) zj=zj+boxzsize
3070 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
3071 c if (lprint_short) then
3072 c write (iout,*) i,j,xi,yi,zi,xj,yj,zj
3073 c write (iout,*) "dist_init",dsqrt(dist_init)
3082 xj=xj_safe+xshift*boxxsize
3083 yj=yj_safe+yshift*boxysize
3084 zj=zj_safe+zshift*boxzsize
3085 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
3086 if(dist_temp.lt.dist_init) then
3096 c if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp)
3097 if (subchap.eq.1) then
3106 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
3107 c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
3108 c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
3109 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
3110 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),
3112 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
3113 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
3114 c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),
3115 c & " subchap",subchap," sss",sss
3116 if (sss.gt.0.0d0) then
3119 e1=fac*fac*aad(itypj,iteli)
3120 e2=fac*bad(itypj,iteli)
3121 if (iabs(j-i) .le. 2) then
3124 evdw2_14=evdw2_14+(e1+e2)*sss
3127 evdw2=evdw2+evdwij*sss
3128 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
3129 & 'evdw2',i,j,sss,evdwij
3131 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
3133 fac=-(evdwij+e1)*rrij*sss
3134 fac=fac+evdwij*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)/expon
3138 C Uncomment following three lines for SC-p interactions
3140 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3142 C Uncomment following line for SC-p interactions
3143 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3145 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
3146 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
3155 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
3156 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
3157 gradx_scp(j,i)=expon*gradx_scp(j,i)
3160 C******************************************************************************
3164 C To save time the factor EXPON has been extracted from ALL components
3165 C of GVDWC and GRADX. Remember to multiply them by this factor before further
3168 C******************************************************************************