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,ikont
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
90 c do i=iatsc_s,iatsc_e
91 do ikont=g_listscsc_start,g_listscsc_end
95 if (itypi.eq.ntyp1) cycle
96 itypi1=iabs(itype(i+1))
101 C Calculate SC interaction energy.
103 c do iint=1,nint_gr(i)
104 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
105 cd & 'iend=',iend(i,iint)
106 c do j=istart(i,iint),iend(i,iint)
108 if (itypj.eq.ntyp1) cycle
112 rij=xj*xj+yj*yj+zj*zj
114 eps0ij=eps(itypi,itypj)
115 sss1=sscale(sqrij,r_cut_int)
116 if (sss1.eq.0.0d0) cycle
117 sssgrad1=sscagrad(sqrij,r_cut_int)
119 & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
120 sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
121 if (sss.lt.1.0d0) then
127 evdw=evdw+(1.0d0-sss)*sss1*evdwij/sqrij/expon
129 C Calculate the components of the gradient in DC and X
131 fac=-rrij*(e1+evdwij)*(1.0d0-sss)*sss1
132 & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
133 & +(1.0d0-sss)*sssgrad1)/sqrij
138 gvdwx(k,i)=gvdwx(k,i)-gg(k)
139 gvdwx(k,j)=gvdwx(k,j)+gg(k)
140 gvdwc(k,i)=gvdwc(k,i)-gg(k)
141 gvdwc(k,j)=gvdwc(k,j)+gg(k)
149 gvdwc(j,i)=expon*gvdwc(j,i)
150 gvdwx(j,i)=expon*gvdwx(j,i)
153 C******************************************************************************
157 C To save time, the factor of EXPON has been extracted from ALL components
158 C of GVDWC and GRADX. Remember to multiply them by this factor before further
161 C******************************************************************************
164 C-----------------------------------------------------------------------
165 subroutine elj_short(evdw)
167 C This subroutine calculates the interaction energy of nonbonded side chains
168 C assuming the LJ potential of interaction.
174 include 'COMMON.LOCAL'
175 include 'COMMON.CHAIN'
176 include 'COMMON.DERIV'
177 include 'COMMON.INTERACT'
178 include 'COMMON.TORSION'
179 include 'COMMON.SBRIDGE'
180 include 'COMMON.NAMES'
181 include 'COMMON.IOUNITS'
182 include "COMMON.SPLITELE"
183 c include 'COMMON.CONTACTS'
184 double precision gg(3)
185 double precision evdw,evdwij
186 integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont
187 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
188 & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
189 double precision sscale,sscagrad
190 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
192 c do i=iatsc_s,iatsc_e
193 do ikont=g_listscsc_start,g_listscsc_end
194 i=newcontlisti(ikont)
195 j=newcontlistj(ikont)
197 if (itypi.eq.ntyp1) cycle
198 itypi1=iabs(itype(i+1))
205 C Calculate SC interaction energy.
207 c do iint=1,nint_gr(i)
208 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
209 cd & 'iend=',iend(i,iint)
210 c do j=istart(i,iint),iend(i,iint)
212 if (itypj.eq.ntyp1) cycle
216 C Change 12/1/95 to calculate four-body interactions
217 rij=xj*xj+yj*yj+zj*zj
219 sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
220 if (sss.gt.0.0d0) then
222 & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
224 eps0ij=eps(itypi,itypj)
231 C Calculate the components of the gradient in DC and X
233 fac=-rrij*(e1+evdwij)*sss+evdwij*sssgrad/sqrij/expon
238 gvdwx(k,i)=gvdwx(k,i)-gg(k)
239 gvdwx(k,j)=gvdwx(k,j)+gg(k)
240 gvdwc(k,i)=gvdwc(k,i)-gg(k)
241 gvdwc(k,j)=gvdwc(k,j)+gg(k)
249 gvdwc(j,i)=expon*gvdwc(j,i)
250 gvdwx(j,i)=expon*gvdwx(j,i)
253 C******************************************************************************
257 C To save time, the factor of EXPON has been extracted from ALL components
258 C of GVDWC and GRADX. Remember to multiply them by this factor before further
261 C******************************************************************************
264 C-----------------------------------------------------------------------------
265 subroutine eljk_long(evdw)
267 C This subroutine calculates the interaction energy of nonbonded side chains
268 C assuming the LJK potential of interaction.
274 include 'COMMON.LOCAL'
275 include 'COMMON.CHAIN'
276 include 'COMMON.DERIV'
277 include 'COMMON.INTERACT'
278 include 'COMMON.IOUNITS'
279 include 'COMMON.NAMES'
280 include "COMMON.SPLITELE"
281 double precision gg(3)
282 double precision evdw,evdwij
283 integer i,j,k,itypi,itypj,itypi1,iint,ikont
284 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
285 & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
287 double precision sscale,sscagrad
288 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
290 c do i=iatsc_s,iatsc_e
291 do ikont=g_listscsc_start,g_listscsc_end
292 i=newcontlisti(ikont)
293 j=newcontlistj(ikont)
295 if (itypi.eq.ntyp1) cycle
296 itypi1=iabs(itype(i+1))
301 C Calculate SC interaction energy.
303 c do iint=1,nint_gr(i)
304 c do j=istart(i,iint),iend(i,iint)
306 if (itypj.eq.ntyp1) cycle
310 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
312 e_augm=augm(itypi,itypj)*fac_augm
315 sss1=sscale(rij,r_cut_int)
316 if (sss1.eq.0.0d0) cycle
317 sssgrad1=sscagrad(rij,r_cut_int)
318 sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
319 if (sss.lt.1.0d0) then
321 & sscagrad(rij/sigma(itypi,itypj),r_cut_respa)
322 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
323 fac=r_shift_inv**expon
327 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
328 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
329 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
330 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
331 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
332 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
333 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
334 evdw=evdw+(1.0d0-sss)*sss1*evdwij
336 C Calculate the components of the gradient in DC and X
338 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
339 fac=fac*(1.0d0-sss)*sss1
340 & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
341 & +(1.0d0-sss)*sssgrad1)*r_inv_ij/expon
346 gvdwx(k,i)=gvdwx(k,i)-gg(k)
347 gvdwx(k,j)=gvdwx(k,j)+gg(k)
348 gvdwc(k,i)=gvdwc(k,i)-gg(k)
349 gvdwc(k,j)=gvdwc(k,j)+gg(k)
357 gvdwc(j,i)=expon*gvdwc(j,i)
358 gvdwx(j,i)=expon*gvdwx(j,i)
363 C-----------------------------------------------------------------------------
364 subroutine eljk_short(evdw)
366 C This subroutine calculates the interaction energy of nonbonded side chains
367 C assuming the LJK potential of interaction.
369 implicit real*8 (a-h,o-z)
373 include 'COMMON.LOCAL'
374 include 'COMMON.CHAIN'
375 include 'COMMON.DERIV'
376 include 'COMMON.INTERACT'
377 include 'COMMON.IOUNITS'
378 include 'COMMON.NAMES'
379 include "COMMON.SPLITELE"
380 double precision gg(3)
381 double precision evdw,evdwij
382 integer i,j,k,itypi,itypj,itypi1,iint,ikont
383 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
384 & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
386 double precision sscale,sscagrad
387 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
389 c do i=iatsc_s,iatsc_e
390 do ikont=g_listscsc_start,g_listscsc_end
391 i=newcontlisti(ikont)
392 j=newcontlistj(ikont)
394 if (itypi.eq.ntyp1) cycle
395 itypi1=iabs(itype(i+1))
400 C Calculate SC interaction energy.
402 c do iint=1,nint_gr(i)
403 c do j=istart(i,iint),iend(i,iint)
405 if (itypj.eq.ntyp1) cycle
409 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
411 e_augm=augm(itypi,itypj)*fac_augm
414 sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
415 if (sss.gt.0.0d0) then
416 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
417 fac=r_shift_inv**expon
421 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
422 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
423 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
424 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
425 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
426 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
427 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
430 C Calculate the components of the gradient in DC and X
432 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
433 & +evdwij*sssgrad/sigma(itypi,itypj)*r_inv_ij/expon
439 gvdwx(k,i)=gvdwx(k,i)-gg(k)
440 gvdwx(k,j)=gvdwx(k,j)+gg(k)
441 gvdwc(k,i)=gvdwc(k,i)-gg(k)
442 gvdwc(k,j)=gvdwc(k,j)+gg(k)
450 gvdwc(j,i)=expon*gvdwc(j,i)
451 gvdwx(j,i)=expon*gvdwx(j,i)
456 C-----------------------------------------------------------------------------
457 subroutine ebp_long(evdw)
459 C This subroutine calculates the interaction energy of nonbonded side chains
460 C assuming the Berne-Pechukas potential of interaction.
466 include 'COMMON.LOCAL'
467 include 'COMMON.CHAIN'
468 include 'COMMON.DERIV'
469 include 'COMMON.NAMES'
470 include 'COMMON.INTERACT'
471 include 'COMMON.IOUNITS'
472 include 'COMMON.CALC'
473 include "COMMON.SPLITELE"
476 double precision evdw
477 integer itypi,itypj,itypi1,iint,ind,ikont
478 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
479 double precision sss1,sssgrad1
480 double precision sscale,sscagrad
481 c double precision rrsave(maxdim)
484 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
486 c if (icall.eq.0) then
492 c do i=iatsc_s,iatsc_e
493 do ikont=g_listscsc_start,g_listscsc_end
494 i=newcontlisti(ikont)
495 j=newcontlistj(ikont)
497 if (itypi.eq.ntyp1) cycle
498 itypi1=iabs(itype(i+1))
502 dxi=dc_norm(1,nres+i)
503 dyi=dc_norm(2,nres+i)
504 dzi=dc_norm(3,nres+i)
505 c dsci_inv=dsc_inv(itypi)
506 dsci_inv=vbld_inv(i+nres)
508 C Calculate SC interaction energy.
510 c do iint=1,nint_gr(i)
511 c do j=istart(i,iint),iend(i,iint)
514 if (itypj.eq.ntyp1) cycle
515 c dscj_inv=dsc_inv(itypj)
516 dscj_inv=vbld_inv(j+nres)
517 chi1=chi(itypi,itypj)
518 chi2=chi(itypj,itypi)
525 alf12=0.5D0*(alf1+alf2)
529 dxj=dc_norm(1,nres+j)
530 dyj=dc_norm(2,nres+j)
531 dzj=dc_norm(3,nres+j)
532 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
534 sss1=sscale(1.0d0/rij,r_cut_int)
535 if (sss1.eq.0.0d0) cycle
536 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
537 if (sss.lt.1.0d0) then
539 & sscagrad(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
540 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
541 C Calculate the angle-dependent terms of energy & contributions to derivatives.
543 C Calculate whole angle-dependent part of epsilon and contributions
545 fac=(rrij*sigsq)**expon2
548 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
549 eps2der=evdwij*eps3rt
550 eps3der=evdwij*eps2rt
551 evdwij=evdwij*eps2rt*eps3rt
552 evdw=evdw+evdwij*(1.0d0-sss)*sss1
554 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
556 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
557 cd & restyp(itypi),i,restyp(itypj),j,
558 cd & epsi,sigm,chi1,chi2,chip1,chip2,
559 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
560 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
563 C Calculate gradient components.
564 e1=e1*eps1*eps2rt**2*eps3rt**2
565 fac=-expon*(e1+evdwij)
567 fac=(fac+evdwij*(sss1/(1.0d0-sss)*sssgrad/
568 & sigmaii(itypi,itypj)+(1.0d0-sss)/sss1*sssgrad1))*rij
569 C Calculate radial part of the gradient
573 C Calculate the angular part of the gradient and sum add the contributions
574 C to the appropriate components of the Cartesian gradient.
575 call sc_grad_scale((1.0d0-sss)*sss1)
583 C-----------------------------------------------------------------------------
584 subroutine ebp_short(evdw)
586 C This subroutine calculates the interaction energy of nonbonded side chains
587 C assuming the Berne-Pechukas potential of interaction.
589 implicit real*8 (a-h,o-z)
593 include 'COMMON.LOCAL'
594 include 'COMMON.CHAIN'
595 include 'COMMON.DERIV'
596 include 'COMMON.NAMES'
597 include 'COMMON.INTERACT'
598 include 'COMMON.IOUNITS'
599 include 'COMMON.CALC'
600 include "COMMON.SPLITELE"
603 double precision evdw
604 integer itypi,itypj,itypi1,iint,ind,ikont
605 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
606 double precision sscale,sscagrad
607 c double precision rrsave(maxdim)
610 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
612 c if (icall.eq.0) then
618 c do i=iatsc_s,iatsc_e
619 do ikont=g_listscsc_start,g_listscsc_end
620 i=newcontlisti(ikont)
621 j=newcontlistj(ikont)
623 if (itypi.eq.ntyp1) cycle
624 itypi1=iabs(itype(i+1))
628 dxi=dc_norm(1,nres+i)
629 dyi=dc_norm(2,nres+i)
630 dzi=dc_norm(3,nres+i)
631 c dsci_inv=dsc_inv(itypi)
632 dsci_inv=vbld_inv(i+nres)
634 C Calculate SC interaction energy.
636 c do iint=1,nint_gr(i)
637 c do j=istart(i,iint),iend(i,iint)
640 if (itypj.eq.ntyp1) cycle
641 c dscj_inv=dsc_inv(itypj)
642 dscj_inv=vbld_inv(j+nres)
643 chi1=chi(itypi,itypj)
644 chi2=chi(itypj,itypi)
651 alf12=0.5D0*(alf1+alf2)
655 dxj=dc_norm(1,nres+j)
656 dyj=dc_norm(2,nres+j)
657 dzj=dc_norm(3,nres+j)
658 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
660 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
662 if (sss.gt.0.0d0) then
664 C Calculate the angle-dependent terms of energy & contributions to derivatives.
666 C Calculate whole angle-dependent part of epsilon and contributions
668 fac=(rrij*sigsq)**expon2
671 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
672 eps2der=evdwij*eps3rt
673 eps3der=evdwij*eps2rt
674 evdwij=evdwij*eps2rt*eps3rt
677 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
679 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
680 cd & restyp(itypi),i,restyp(itypj),j,
681 cd & epsi,sigm,chi1,chi2,chip1,chip2,
682 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
683 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
686 C Calculate gradient components.
687 e1=e1*eps1*eps2rt**2*eps3rt**2
688 fac=-expon*(e1+evdwij)
690 fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rrij
691 C Calculate radial part of the gradient
695 C Calculate the angular part of the gradient and sum add the contributions
696 C to the appropriate components of the Cartesian gradient.
697 call sc_grad_scale(sss)
705 C-----------------------------------------------------------------------------
706 subroutine egb_long(evdw)
708 C This subroutine calculates the interaction energy of nonbonded side chains
709 C assuming the Gay-Berne potential of interaction.
715 include 'COMMON.LOCAL'
716 include 'COMMON.CHAIN'
717 include 'COMMON.DERIV'
718 include 'COMMON.NAMES'
719 include 'COMMON.INTERACT'
720 include 'COMMON.IOUNITS'
721 include 'COMMON.CALC'
722 include 'COMMON.CONTROL'
723 include "COMMON.SPLITELE"
725 integer xshift,yshift,zshift
726 double precision evdw
727 integer itypi,itypj,itypi1,iint,ind,ikont
728 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
729 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
730 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
731 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
732 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
733 double precision subchap,sss1,sssgrad1
735 ccccc energy_dec=.false.
736 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
739 c if (icall.eq.0) lprn=.false.
741 c do i=iatsc_s,iatsc_e
742 do ikont=g_listscsc_start,g_listscsc_end
743 i=newcontlisti(ikont)
744 j=newcontlistj(ikont)
746 if (itypi.eq.ntyp1) cycle
747 itypi1=iabs(itype(i+1))
752 if (xi.lt.0) xi=xi+boxxsize
754 if (yi.lt.0) yi=yi+boxysize
756 if (zi.lt.0) zi=zi+boxzsize
757 dxi=dc_norm(1,nres+i)
758 dyi=dc_norm(2,nres+i)
759 dzi=dc_norm(3,nres+i)
760 c dsci_inv=dsc_inv(itypi)
761 dsci_inv=vbld_inv(i+nres)
762 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
763 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
765 C Calculate SC interaction energy.
767 c do iint=1,nint_gr(i)
768 c do j=istart(i,iint),iend(i,iint)
771 if (itypj.eq.ntyp1) cycle
772 c dscj_inv=dsc_inv(itypj)
773 dscj_inv=vbld_inv(j+nres)
774 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
775 c & 1.0d0/vbld(j+nres)
776 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
777 sig0ij=sigma(itypi,itypj)
778 chi1=chi(itypi,itypj)
779 chi2=chi(itypj,itypi)
786 alf12=0.5D0*(alf1+alf2)
791 if (xj.lt.0) xj=xj+boxxsize
793 if (yj.lt.0) yj=yj+boxysize
795 if (zj.lt.0) zj=zj+boxzsize
796 if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
797 C the energy transfer exist
798 if (zj.lt.buflipbot) then
799 C what fraction I am in
800 fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
801 C lipbufthick is thickenes of lipid buffore
802 sslipj=sscalelip(fracinbuf)
803 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
804 else if (zi.gt.bufliptop) then
805 fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
806 sslipj=sscalelip(fracinbuf)
807 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
816 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
817 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
818 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
819 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
820 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
828 xj=xj_safe+xshift*boxxsize
829 yj=yj_safe+yshift*boxysize
830 zj=zj_safe+zshift*boxzsize
831 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
832 if (dist_temp.lt.dist_init) then
842 if (subchap.eq.1) then
851 dxj=dc_norm(1,nres+j)
852 dyj=dc_norm(2,nres+j)
853 dzj=dc_norm(3,nres+j)
854 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
856 sss1=sscale(1.0d0/rij,r_cut_int)
857 if (sss1.eq.0.0d0) cycle
858 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
859 if (sss.lt.1.0d0) then
860 C Calculate angle-dependent terms of energy and contributions to their
863 & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
864 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
867 sig=sig0ij*dsqrt(sigsq)
868 rij_shift=1.0D0/rij-sig+sig0ij
869 c for diagnostics; uncomment
870 c rij_shift=1.2*sig0ij
871 C I hate to put IF's in the loops, but here don't have another choice!!!!
872 if (rij_shift.le.0.0D0) then
874 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
875 cd & restyp(itypi),i,restyp(itypj),j,
876 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
880 c---------------------------------------------------------------
881 rij_shift=1.0D0/rij_shift
885 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
886 eps2der=evdwij*eps3rt
887 eps3der=evdwij*eps2rt
888 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
889 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
890 evdwij=evdwij*eps2rt*eps3rt
891 evdw=evdw+evdwij*(1.0d0-sss)*sss1
893 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
895 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
896 & restyp(itypi),i,restyp(itypj),j,
897 & epsi,sigm,chi1,chi2,chip1,chip2,
898 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
899 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
903 if (energy_dec) write (iout,'(a6,2i5,4f10.5)')
904 & 'evdw',i,j,rij,sss,sss1,evdwij
906 C Calculate gradient components.
907 e1=e1*eps1*eps2rt**2*eps3rt**2
908 fac=-expon*(e1+evdwij)*rij_shift
910 fac=(fac+evdwij*(-sss1*sssgrad/(1.0d0-sss)
911 & /sigmaii(itypi,itypj)+(1.0d0-sss)*sssgrad1/sss1))*rij
913 C Calculate the radial part of the gradient
917 gg_lipi(3)=ssgradlipi*evdwij
918 gg_lipj(3)=ssgradlipj*evdwij
919 C Calculate angular part of the gradient.
920 call sc_grad_scale((1.0d0-sss)*sss1)
925 c write (iout,*) "Number of loop steps in EGB:",ind
926 cccc energy_dec=.false.
929 C-----------------------------------------------------------------------------
930 subroutine egb_short(evdw)
932 C This subroutine calculates the interaction energy of nonbonded side chains
933 C assuming the Gay-Berne potential of interaction.
939 include 'COMMON.LOCAL'
940 include 'COMMON.CHAIN'
941 include 'COMMON.DERIV'
942 include 'COMMON.NAMES'
943 include 'COMMON.INTERACT'
944 include 'COMMON.IOUNITS'
945 include 'COMMON.CALC'
946 include 'COMMON.CONTROL'
947 include "COMMON.SPLITELE"
949 integer xshift,yshift,zshift
950 double precision evdw
951 integer itypi,itypj,itypi1,iint,ind,ikont
952 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
953 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
954 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
955 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
956 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
957 double precision subchap
959 ccccc energy_dec=.false.
960 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
963 c if (icall.eq.0) lprn=.false.
965 c do i=iatsc_s,iatsc_e
966 do ikont=g_listscsc_start,g_listscsc_end
967 i=newcontlisti(ikont)
968 j=newcontlistj(ikont)
970 if (itypi.eq.ntyp1) cycle
971 itypi1=iabs(itype(i+1))
976 if (xi.lt.0) xi=xi+boxxsize
978 if (yi.lt.0) yi=yi+boxysize
980 if (zi.lt.0) zi=zi+boxzsize
981 dxi=dc_norm(1,nres+i)
982 dyi=dc_norm(2,nres+i)
983 dzi=dc_norm(3,nres+i)
984 c dsci_inv=dsc_inv(itypi)
985 dsci_inv=vbld_inv(i+nres)
986 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
987 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
989 C Calculate SC interaction energy.
991 c do iint=1,nint_gr(i)
992 c do j=istart(i,iint),iend(i,iint)
995 if (itypj.eq.ntyp1) cycle
996 c dscj_inv=dsc_inv(itypj)
997 dscj_inv=vbld_inv(j+nres)
998 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
999 c & 1.0d0/vbld(j+nres)
1000 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1001 sig0ij=sigma(itypi,itypj)
1002 chi1=chi(itypi,itypj)
1003 chi2=chi(itypj,itypi)
1010 alf12=0.5D0*(alf1+alf2)
1015 if (xj.lt.0) xj=xj+boxxsize
1017 if (yj.lt.0) yj=yj+boxysize
1019 if (zj.lt.0) zj=zj+boxzsize
1020 if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
1021 C the energy transfer exist
1022 if (zj.lt.buflipbot) then
1023 C what fraction I am in
1024 fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
1025 C lipbufthick is thickenes of lipid buffore
1026 sslipj=sscalelip(fracinbuf)
1027 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
1028 elseif (zi.gt.bufliptop) then
1029 fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
1030 sslipj=sscalelip(fracinbuf)
1031 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
1040 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1041 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
1042 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1043 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
1044 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1052 xj=xj_safe+xshift*boxxsize
1053 yj=yj_safe+yshift*boxysize
1054 zj=zj_safe+zshift*boxzsize
1055 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1056 if(dist_temp.lt.dist_init) then
1066 if (subchap.eq.1) then
1075 dxj=dc_norm(1,nres+j)
1076 dyj=dc_norm(2,nres+j)
1077 dzj=dc_norm(3,nres+j)
1078 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1080 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1081 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1082 if (sss.gt.0.0d0) then
1084 C Calculate angle-dependent terms of energy and contributions to their
1088 sig=sig0ij*dsqrt(sigsq)
1089 rij_shift=1.0D0/rij-sig+sig0ij
1090 c for diagnostics; uncomment
1091 c rij_shift=1.2*sig0ij
1092 C I hate to put IF's in the loops, but here don't have another choice!!!!
1093 if (rij_shift.le.0.0D0) then
1095 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1096 cd & restyp(itypi),i,restyp(itypj),j,
1097 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1101 c---------------------------------------------------------------
1102 rij_shift=1.0D0/rij_shift
1103 fac=rij_shift**expon
1106 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1107 eps2der=evdwij*eps3rt
1108 eps3der=evdwij*eps2rt
1109 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1110 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1111 evdwij=evdwij*eps2rt*eps3rt
1112 evdw=evdw+evdwij*sss
1114 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1116 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1117 & restyp(itypi),i,restyp(itypj),j,
1118 & epsi,sigm,chi1,chi2,chip1,chip2,
1119 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1120 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1124 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1127 C Calculate gradient components.
1128 e1=e1*eps1*eps2rt**2*eps3rt**2
1129 fac=-expon*(e1+evdwij)*rij_shift
1131 fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rij
1133 C Calculate the radial part of the gradient
1137 gg_lipi(3)=ssgradlipi*evdwij
1138 gg_lipj(3)=ssgradlipj*evdwij
1139 C Calculate angular part of the gradient.
1140 call sc_grad_scale(sss)
1145 c write (iout,*) "Number of loop steps in EGB:",ind
1146 cccc energy_dec=.false.
1149 C-----------------------------------------------------------------------------
1150 subroutine egbv_long(evdw)
1152 C This subroutine calculates the interaction energy of nonbonded side chains
1153 C assuming the Gay-Berne-Vorobjev potential of interaction.
1155 implicit real*8 (a-h,o-z)
1156 include 'DIMENSIONS'
1157 include 'COMMON.GEO'
1158 include 'COMMON.VAR'
1159 include 'COMMON.LOCAL'
1160 include 'COMMON.CHAIN'
1161 include 'COMMON.DERIV'
1162 include 'COMMON.NAMES'
1163 include 'COMMON.INTERACT'
1164 include 'COMMON.IOUNITS'
1165 include 'COMMON.CALC'
1166 include "COMMON.SPLITELE"
1168 common /srutu/ icall
1170 integer itypi,itypj,itypi1,iint,ind,ikont
1171 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1172 & xi,yi,zi,fac_augm,e_augm
1173 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1174 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1175 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1176 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1177 double precision sss1,sssgrad1
1179 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1182 c if (icall.eq.0) lprn=.true.
1184 c do i=iatsc_s,iatsc_e
1185 do ikont=g_listscsc_start,g_listscsc_end
1186 i=newcontlisti(ikont)
1187 j=newcontlistj(ikont)
1188 itypi=iabs(itype(i))
1189 if (itypi.eq.ntyp1) cycle
1190 itypi1=iabs(itype(i+1))
1194 dxi=dc_norm(1,nres+i)
1195 dyi=dc_norm(2,nres+i)
1196 dzi=dc_norm(3,nres+i)
1197 c dsci_inv=dsc_inv(itypi)
1198 dsci_inv=vbld_inv(i+nres)
1200 C Calculate SC interaction energy.
1202 c do iint=1,nint_gr(i)
1203 c do j=istart(i,iint),iend(i,iint)
1205 itypj=iabs(itype(j))
1206 if (itypj.eq.ntyp1) cycle
1207 c dscj_inv=dsc_inv(itypj)
1208 dscj_inv=vbld_inv(j+nres)
1209 sig0ij=sigma(itypi,itypj)
1210 r0ij=r0(itypi,itypj)
1211 chi1=chi(itypi,itypj)
1212 chi2=chi(itypj,itypi)
1219 alf12=0.5D0*(alf1+alf2)
1223 dxj=dc_norm(1,nres+j)
1224 dyj=dc_norm(2,nres+j)
1225 dzj=dc_norm(3,nres+j)
1226 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1229 sss1=sscale(1.0d0/rij,r_cut_int)
1230 if (sss1.eq.0.0d0) cycle
1232 if (sss.lt.1.0d0) then
1234 & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1235 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
1237 C Calculate angle-dependent terms of energy and contributions to their
1241 sig=sig0ij*dsqrt(sigsq)
1242 rij_shift=1.0D0/rij-sig+r0ij
1243 C I hate to put IF's in the loops, but here don't have another choice!!!!
1244 if (rij_shift.le.0.0D0) then
1249 c---------------------------------------------------------------
1250 rij_shift=1.0D0/rij_shift
1251 fac=rij_shift**expon
1254 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1255 eps2der=evdwij*eps3rt
1256 eps3der=evdwij*eps2rt
1257 fac_augm=rrij**expon
1258 e_augm=augm(itypi,itypj)*fac_augm
1259 evdwij=evdwij*eps2rt*eps3rt
1260 evdw=evdw+(evdwij+e_augm)*sss1*(1.0d0-sss)
1262 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1264 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1265 & restyp(itypi),i,restyp(itypj),j,
1266 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1267 & chi1,chi2,chip1,chip2,
1268 & eps1,eps2rt**2,eps3rt**2,
1269 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1272 C Calculate gradient components.
1273 e1=e1*eps1*eps2rt**2*eps3rt**2
1274 fac=-expon*(e1+evdwij)*rij_shift
1276 fac=rij*fac-2*expon*rrij*e_augm
1277 fac=fac+(evdwij+e_augm)*
1278 & (-sss1*sssgrad/(1.0d0-sss)/sigmaii(itypi,itypj)
1279 & +(1.0d0-sss)*sssgrad1/sss1)*rij
1280 C Calculate the radial part of the gradient
1284 C Calculate angular part of the gradient.
1285 call sc_grad_scale((1.0d0-sss)*sss1)
1291 C-----------------------------------------------------------------------------
1292 subroutine egbv_short(evdw)
1294 C This subroutine calculates the interaction energy of nonbonded side chains
1295 C assuming the Gay-Berne-Vorobjev potential of interaction.
1298 include 'DIMENSIONS'
1299 include 'COMMON.GEO'
1300 include 'COMMON.VAR'
1301 include 'COMMON.LOCAL'
1302 include 'COMMON.CHAIN'
1303 include 'COMMON.DERIV'
1304 include 'COMMON.NAMES'
1305 include 'COMMON.INTERACT'
1306 include 'COMMON.IOUNITS'
1307 include 'COMMON.CALC'
1308 include "COMMON.SPLITELE"
1310 common /srutu/ icall
1312 integer itypi,itypj,itypi1,iint,ind,ikont
1313 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1314 & xi,yi,zi,fac_augm,e_augm
1315 double precision evdw
1316 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1317 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1318 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1319 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1321 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1324 c if (icall.eq.0) lprn=.true.
1326 c do i=iatsc_s,iatsc_e
1327 do ikont=g_listscsc_start,g_listscsc_end
1328 i=newcontlisti(ikont)
1329 j=newcontlistj(ikont)
1330 itypi=iabs(itype(i))
1331 if (itypi.eq.ntyp1) cycle
1332 itypi1=iabs(itype(i+1))
1336 dxi=dc_norm(1,nres+i)
1337 dyi=dc_norm(2,nres+i)
1338 dzi=dc_norm(3,nres+i)
1339 c dsci_inv=dsc_inv(itypi)
1340 dsci_inv=vbld_inv(i+nres)
1342 C Calculate SC interaction energy.
1344 c do iint=1,nint_gr(i)
1345 c do j=istart(i,iint),iend(i,iint)
1347 itypj=iabs(itype(j))
1348 if (itypj.eq.ntyp1) cycle
1349 c dscj_inv=dsc_inv(itypj)
1350 dscj_inv=vbld_inv(j+nres)
1351 sig0ij=sigma(itypi,itypj)
1352 r0ij=r0(itypi,itypj)
1353 chi1=chi(itypi,itypj)
1354 chi2=chi(itypj,itypi)
1361 alf12=0.5D0*(alf1+alf2)
1365 dxj=dc_norm(1,nres+j)
1366 dyj=dc_norm(2,nres+j)
1367 dzj=dc_norm(3,nres+j)
1368 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1371 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1373 if (sss.gt.0.0d0) then
1375 C Calculate angle-dependent terms of energy and contributions to their
1379 sig=sig0ij*dsqrt(sigsq)
1380 rij_shift=1.0D0/rij-sig+r0ij
1381 C I hate to put IF's in the loops, but here don't have another choice!!!!
1382 if (rij_shift.le.0.0D0) then
1387 c---------------------------------------------------------------
1388 rij_shift=1.0D0/rij_shift
1389 fac=rij_shift**expon
1392 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1393 eps2der=evdwij*eps3rt
1394 eps3der=evdwij*eps2rt
1395 fac_augm=rrij**expon
1396 e_augm=augm(itypi,itypj)*fac_augm
1397 evdwij=evdwij*eps2rt*eps3rt
1398 evdw=evdw+(evdwij+e_augm)*sss
1400 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1402 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1403 & restyp(itypi),i,restyp(itypj),j,
1404 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1405 & chi1,chi2,chip1,chip2,
1406 & eps1,eps2rt**2,eps3rt**2,
1407 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1410 C Calculate gradient components.
1411 e1=e1*eps1*eps2rt**2*eps3rt**2
1412 fac=-expon*(e1+evdwij)*rij_shift
1414 fac=rij*fac-2*expon*rrij*e_augm+
1415 & (evdwij+e_augm)*sssgrad/sigmaii(itypi,itypj)/sss*rij
1416 C Calculate the radial part of the gradient
1420 C Calculate angular part of the gradient.
1421 call sc_grad_scale(sss)
1427 C----------------------------------------------------------------------------
1428 subroutine sc_grad_scale(scalfac)
1429 implicit real*8 (a-h,o-z)
1430 include 'DIMENSIONS'
1431 include 'COMMON.CHAIN'
1432 include 'COMMON.DERIV'
1433 include 'COMMON.CALC'
1434 include 'COMMON.IOUNITS'
1435 include "COMMON.SPLITELE"
1436 double precision dcosom1(3),dcosom2(3)
1437 double precision scalfac
1438 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1439 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1440 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1441 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1445 c eom12=evdwij*eps1_om12
1447 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1448 c & " sigder",sigder
1449 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1450 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1452 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1453 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1456 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1458 c write (iout,*) "gg",(gg(k),k=1,3)
1460 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
1461 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1462 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1463 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
1464 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1465 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1466 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1467 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1468 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1469 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1472 C Calculate the components of the gradient in DC and X
1475 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
1476 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
1480 C--------------------------------------------------------------------------
1481 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1483 C This subroutine calculates the average interaction energy and its gradient
1484 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1485 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1486 C The potential depends both on the distance of peptide-group centers and on
1487 C the orientation of the CA-CA virtual bonds.
1489 implicit real*8 (a-h,o-z)
1493 include 'DIMENSIONS'
1494 include 'COMMON.CONTROL'
1495 include 'COMMON.SETUP'
1496 include 'COMMON.IOUNITS'
1497 include 'COMMON.GEO'
1498 include 'COMMON.VAR'
1499 include 'COMMON.LOCAL'
1500 include 'COMMON.CHAIN'
1501 include 'COMMON.DERIV'
1502 include 'COMMON.INTERACT'
1504 include 'COMMON.CONTACTS'
1505 include 'COMMON.CONTMAT'
1507 include 'COMMON.CORRMAT'
1508 include 'COMMON.TORSION'
1509 include 'COMMON.VECTORS'
1510 include 'COMMON.FFIELD'
1511 include 'COMMON.TIME1'
1512 include 'COMMON.SHIELD'
1513 include "COMMON.SPLITELE"
1515 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1516 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1517 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1518 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1519 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1520 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1522 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1524 double precision scal_el /1.0d0/
1526 double precision scal_el /0.5d0/
1529 C 13-go grudnia roku pamietnego...
1530 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1531 & 0.0d0,1.0d0,0.0d0,
1532 & 0.0d0,0.0d0,1.0d0/
1533 cd write(iout,*) 'In EELEC'
1535 cd write(iout,*) 'Type',i
1536 cd write(iout,*) 'B1',B1(:,i)
1537 cd write(iout,*) 'B2',B2(:,i)
1538 cd write(iout,*) 'CC',CC(:,:,i)
1539 cd write(iout,*) 'DD',DD(:,:,i)
1540 cd write(iout,*) 'EE',EE(:,:,i)
1542 cd call check_vecgrad
1544 C print *,"WCHODZE3"
1545 if (icheckgrad.eq.1) then
1547 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1549 dc_norm(k,i)=dc(k,i)*fac
1551 c write (iout,*) 'i',i,' fac',fac
1554 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1555 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1556 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1557 c call vec_and_deriv
1563 time_mat=time_mat+MPI_Wtime()-time01
1567 cd write (iout,*) 'i=',i
1569 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1572 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1573 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1588 cd print '(a)','Enter EELEC'
1589 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1591 gel_loc_loc(i)=0.0d0
1596 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1598 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1600 do i=iturn3_start,iturn3_end
1601 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1602 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1603 C & .or. itype(i-1).eq.ntyp1
1604 C & .or. itype(i+4).eq.ntyp1
1609 dx_normi=dc_norm(1,i)
1610 dy_normi=dc_norm(2,i)
1611 dz_normi=dc_norm(3,i)
1612 xmedi=c(1,i)+0.5d0*dxi
1613 ymedi=c(2,i)+0.5d0*dyi
1614 zmedi=c(3,i)+0.5d0*dzi
1615 xmedi=mod(xmedi,boxxsize)
1616 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1617 ymedi=mod(ymedi,boxysize)
1618 if (ymedi.lt.0) ymedi=ymedi+boxysize
1619 zmedi=mod(zmedi,boxzsize)
1620 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1622 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1623 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1625 num_cont_hb(i)=num_conti
1628 do i=iturn4_start,iturn4_end
1629 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1630 & .or. itype(i+3).eq.ntyp1
1631 & .or. itype(i+4).eq.ntyp1
1632 C & .or. itype(i+5).eq.ntyp1
1633 C & .or. itype(i-1).eq.ntyp1
1638 dx_normi=dc_norm(1,i)
1639 dy_normi=dc_norm(2,i)
1640 dz_normi=dc_norm(3,i)
1641 xmedi=c(1,i)+0.5d0*dxi
1642 ymedi=c(2,i)+0.5d0*dyi
1643 zmedi=c(3,i)+0.5d0*dzi
1644 xmedi=mod(xmedi,boxxsize)
1645 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1646 ymedi=mod(ymedi,boxysize)
1647 if (ymedi.lt.0) ymedi=ymedi+boxysize
1648 zmedi=mod(zmedi,boxzsize)
1649 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1651 num_conti=num_cont_hb(i)
1653 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1654 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1655 & call eturn4(i,eello_turn4)
1657 num_cont_hb(i)=num_conti
1661 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1663 c do i=iatel_s,iatel_e
1664 do ikont=g_listpp_start,g_listpp_end
1665 i=newcontlistppi(ikont)
1666 j=newcontlistppj(ikont)
1667 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1668 C & .or. itype(i+2).eq.ntyp1
1669 C & .or. itype(i-1).eq.ntyp1
1674 dx_normi=dc_norm(1,i)
1675 dy_normi=dc_norm(2,i)
1676 dz_normi=dc_norm(3,i)
1677 xmedi=c(1,i)+0.5d0*dxi
1678 ymedi=c(2,i)+0.5d0*dyi
1679 zmedi=c(3,i)+0.5d0*dzi
1680 xmedi=mod(xmedi,boxxsize)
1681 if (xmedi.lt.0) xmedi=xmedi+boxxsize
1682 ymedi=mod(ymedi,boxysize)
1683 if (ymedi.lt.0) ymedi=ymedi+boxysize
1684 zmedi=mod(zmedi,boxzsize)
1685 if (zmedi.lt.0) zmedi=zmedi+boxzsize
1686 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
1688 num_conti=num_cont_hb(i)
1690 c do j=ielstart(i),ielend(i)
1691 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1692 C & .or.itype(j+2).eq.ntyp1
1693 C & .or.itype(j-1).eq.ntyp1
1695 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1698 num_cont_hb(i)=num_conti
1701 c write (iout,*) "Number of loop steps in EELEC:",ind
1703 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1704 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1706 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1707 ccc eel_loc=eel_loc+eello_turn3
1708 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1711 C-------------------------------------------------------------------------------
1712 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1714 include 'DIMENSIONS'
1718 include 'COMMON.CONTROL'
1719 include 'COMMON.IOUNITS'
1720 include 'COMMON.GEO'
1721 include 'COMMON.VAR'
1722 include 'COMMON.LOCAL'
1723 include 'COMMON.CHAIN'
1724 include 'COMMON.DERIV'
1725 include 'COMMON.INTERACT'
1727 include 'COMMON.CONTACTS'
1728 include 'COMMON.CONTMAT'
1730 include 'COMMON.CORRMAT'
1731 include 'COMMON.TORSION'
1732 include 'COMMON.VECTORS'
1733 include 'COMMON.FFIELD'
1734 include 'COMMON.TIME1'
1735 include 'COMMON.SHIELD'
1736 include "COMMON.SPLITELE"
1737 integer xshift,yshift,zshift
1738 double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1739 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1740 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1741 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
1742 & gmuij2(4),gmuji2(4)
1743 integer j1,j2,num_conti
1744 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1745 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1747 integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ind,itypi,itypj
1748 integer ilist,iresshield
1749 double precision rlocshield
1750 double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
1751 double precision ees,evdw1,eel_loc,aaa,bbb,ael3i
1752 double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,
1753 & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,
1754 & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa,
1755 & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,
1756 & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,
1757 & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,
1758 & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
1759 & ecosgp,ecosam,ecosbm,ecosgm,ghalf,geel_loc_ij,geel_loc_ji,
1760 & dxi,dyi,dzi,a22,a23,a32,a33
1761 double precision dist_init,xmedi,ymedi,zmedi,xj_safe,yj_safe,
1762 & zj_safe,xj_temp,yj_temp,zj_temp,dist_temp,dx_normi,dy_normi,
1764 double precision sss1,sssgrad1
1765 double precision sscale,sscagrad
1766 double precision scalar
1768 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1770 double precision scal_el /1.0d0/
1772 double precision scal_el /0.5d0/
1775 C 13-go grudnia roku pamietnego...
1776 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1777 & 0.0d0,1.0d0,0.0d0,
1778 & 0.0d0,0.0d0,1.0d0/
1779 c time00=MPI_Wtime()
1780 cd write (iout,*) "eelecij",i,j
1781 C print *,"WCHODZE2"
1785 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1786 aaa=app(iteli,itelj)
1787 bbb=bpp(iteli,itelj)
1788 ael6i=ael6(iteli,itelj)
1789 ael3i=ael3(iteli,itelj)
1793 dx_normj=dc_norm(1,j)
1794 dy_normj=dc_norm(2,j)
1795 dz_normj=dc_norm(3,j)
1800 if (xj.lt.0) xj=xj+boxxsize
1802 if (yj.lt.0) yj=yj+boxysize
1804 if (zj.lt.0) zj=zj+boxzsize
1805 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1813 xj=xj_safe+xshift*boxxsize
1814 yj=yj_safe+yshift*boxysize
1815 zj=zj_safe+zshift*boxzsize
1816 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1817 if(dist_temp.lt.dist_init) then
1827 if (isubchap.eq.1) then
1837 rij=xj*xj+yj*yj+zj*zj
1841 c For extracting the short-range part of Evdwpp
1842 sss1=sscale(rij,r_cut_int)
1843 if (sss1.eq.0.0d0) return
1844 sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
1845 sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
1846 sssgrad1=sscagrad(rij,r_cut_int)
1849 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1850 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1851 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1852 fac=cosa-3.0D0*cosb*cosg
1854 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1855 if (j.eq.i+2) ev1=scal_el*ev1
1860 if (shield_mode.eq.0) then
1864 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1866 el1=el1*fac_shield(i)**2*fac_shield(j)**2
1867 el2=el2*fac_shield(i)**2*fac_shield(j)**2
1869 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1870 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1872 evdw1=evdw1+evdwij*(1.0d0-sss)*sss1
1873 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1874 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1875 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1876 cd & xmedi,ymedi,zmedi,xj,yj,zj
1878 if (energy_dec) then
1879 write (iout,'(a6,2i5,0pf7.3,2f7.3)')
1880 & 'evdw1',i,j,evdwij,sss,sss1
1881 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1885 C Calculate contributions to the Cartesian gradient.
1888 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
1889 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1890 facel=-3*rrmij*(el1+eesij)*sss1
1891 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1897 * Radial derivatives. First process both termini of the fragment (i,j)
1899 aux=facel+sssgrad1*(1.0d0-sss)*eesij*rmij
1900 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1907 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1908 & (shield_mode.gt.0)) then
1910 do ilist=1,ishield_list(i)
1911 iresshield=shield_list(ilist,i)
1913 rlocshield=grad_shield_side(k,ilist,i)*eesij*sss1
1914 & /fac_shield(i)*2.0*sss1
1915 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1917 & +grad_shield_loc(k,ilist,i)*eesij*sss1/fac_shield(i)*2.0
1919 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1920 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1921 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1922 C if (iresshield.gt.i) then
1923 C do ishi=i+1,iresshield-1
1924 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1925 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1929 C do ishi=iresshield,i
1930 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1931 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1937 do ilist=1,ishield_list(j)
1938 iresshield=shield_list(ilist,j)
1940 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1942 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1944 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss1
1945 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1950 gshieldc(k,i)=gshieldc(k,i)+
1951 & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
1952 gshieldc(k,j)=gshieldc(k,j)+
1953 & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
1954 gshieldc(k,i-1)=gshieldc(k,i-1)+
1955 & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
1956 gshieldc(k,j-1)=gshieldc(k,j-1)+
1957 & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
1963 c ghalf=0.5D0*ggg(k)
1964 c gelc(k,i)=gelc(k,i)+ghalf
1965 c gelc(k,j)=gelc(k,j)+ghalf
1967 c 9/28/08 AL Gradient compotents will be summed only at the end
1969 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
1970 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
1972 c gelc_long(3,i)=gelc_long(3,i)+
1973 c ssgradlipi*eesij/2.0d0*lipscale**2*sss1
1975 * Loop over residues i+1 thru j-1.
1979 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
1983 & (-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij
1984 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1989 c ghalf=0.5D0*ggg(k)
1990 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
1991 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
1993 c 9/28/08 AL Gradient compotents will be summed only at the end
1995 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
1996 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
1999 * Loop over residues i+1 thru j-1.
2003 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2007 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
2008 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2009 facel=-3*rrmij*(el1+eesij)*sss1
2010 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2012 c facvdw=ev1+evdwij*(1.0d0-sss)*sss1
2015 fac=-3*rrmij*(facvdw+facvdw+facel)
2020 * Radial derivatives. First process both termini of the fragment (i,j)
2022 aux=fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj))
2024 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2032 c ghalf=0.5D0*ggg(k)
2033 c gelc(k,i)=gelc(k,i)+ghalf
2034 c gelc(k,j)=gelc(k,j)+ghalf
2036 c 9/28/08 AL Gradient compotents will be summed only at the end
2038 gelc_long(k,j)=gelc(k,j)+ggg(k)
2039 gelc_long(k,i)=gelc(k,i)-ggg(k)
2042 * Loop over residues i+1 thru j-1.
2046 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
2049 c 9/28/08 AL Gradient compotents will be summed only at the end
2054 & (-sssgrad*sss1/rpp(iteli,itelj)+sssgrad1*(1.0d0-sss))*rmij*evdwij
2059 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2060 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2066 ecosa=2.0D0*fac3*fac1+fac4
2069 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
2070 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
2072 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2073 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2075 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
2076 cd & (dcosg(k),k=1,3)
2078 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1
2079 & *fac_shield(i)**2*fac_shield(j)**2
2080 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2084 c ghalf=0.5D0*ggg(k)
2085 c gelc(k,i)=gelc(k,i)+ghalf
2086 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
2087 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2088 c gelc(k,j)=gelc(k,j)+ghalf
2089 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2090 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2094 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
2099 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
2100 & +ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss1
2101 & *fac_shield(i)**2*fac_shield(j)**2
2102 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2105 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2106 & +ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss1
2107 & *fac_shield(i)**2*fac_shield(j)**2
2108 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2109 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
2110 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
2112 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
2113 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
2114 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2116 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
2117 C energy of a peptide unit is assumed in the form of a second-order
2118 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
2119 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
2120 C are computed for EVERY pair of non-contiguous peptide groups.
2122 if (j.lt.nres-1) then
2133 muij(kkk)=mu(k,i)*mu(l,j)
2135 gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
2136 c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
2137 gmuij2(kkk)=gUb2(k,i)*mu(l,j)
2138 gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
2139 c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
2140 gmuji2(kkk)=mu(k,i)*gUb2(l,j)
2144 cd write (iout,*) 'EELEC: i',i,' j',j
2145 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
2146 cd write(iout,*) 'muij',muij
2147 ury=scalar(uy(1,i),erij)
2148 urz=scalar(uz(1,i),erij)
2149 vry=scalar(uy(1,j),erij)
2150 vrz=scalar(uz(1,j),erij)
2151 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
2152 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
2153 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
2154 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
2155 fac=dsqrt(-ael6i)*r3ij
2160 cd write (iout,'(4i5,4f10.5)')
2161 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
2162 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
2163 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
2164 cd & uy(:,j),uz(:,j)
2165 cd write (iout,'(4f10.5)')
2166 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
2167 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
2168 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
2169 cd write (iout,'(9f10.5/)')
2170 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
2171 C Derivatives of the elements of A in virtual-bond vectors
2172 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
2174 uryg(k,1)=scalar(erder(1,k),uy(1,i))
2175 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
2176 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
2177 urzg(k,1)=scalar(erder(1,k),uz(1,i))
2178 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
2179 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
2180 vryg(k,1)=scalar(erder(1,k),uy(1,j))
2181 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
2182 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
2183 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
2184 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
2185 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
2187 C Compute radial contributions to the gradient
2205 C Add the contributions coming from er
2208 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
2209 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
2210 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
2211 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
2214 C Derivatives in DC(i)
2215 cgrad ghalf1=0.5d0*agg(k,1)
2216 cgrad ghalf2=0.5d0*agg(k,2)
2217 cgrad ghalf3=0.5d0*agg(k,3)
2218 cgrad ghalf4=0.5d0*agg(k,4)
2219 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2220 & -3.0d0*uryg(k,2)*vry)!+ghalf1
2221 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2222 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
2223 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2224 & -3.0d0*urzg(k,2)*vry)!+ghalf3
2225 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2226 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
2227 C Derivatives in DC(i+1)
2228 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2229 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2230 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2231 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2232 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2233 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2234 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2235 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2236 C Derivatives in DC(j)
2237 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2238 & -3.0d0*vryg(k,2)*ury)!+ghalf1
2239 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2240 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
2241 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2242 & -3.0d0*vryg(k,2)*urz)!+ghalf3
2243 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2244 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
2245 C Derivatives in DC(j+1) or DC(nres-1)
2246 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2247 & -3.0d0*vryg(k,3)*ury)
2248 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2249 & -3.0d0*vrzg(k,3)*ury)
2250 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2251 & -3.0d0*vryg(k,3)*urz)
2252 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2253 & -3.0d0*vrzg(k,3)*urz)
2254 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
2256 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
2269 aggi(k,l)=-aggi(k,l)
2270 aggi1(k,l)=-aggi1(k,l)
2271 aggj(k,l)=-aggj(k,l)
2272 aggj1(k,l)=-aggj1(k,l)
2275 if (j.lt.nres-1) then
2281 aggi(k,l)=-aggi(k,l)
2282 aggi1(k,l)=-aggi1(k,l)
2283 aggj(k,l)=-aggj(k,l)
2284 aggj1(k,l)=-aggj1(k,l)
2295 aggi(k,l)=-aggi(k,l)
2296 aggi1(k,l)=-aggi1(k,l)
2297 aggj(k,l)=-aggj(k,l)
2298 aggj1(k,l)=-aggj1(k,l)
2303 IF (wel_loc.gt.0.0d0) THEN
2304 C Contribution to the local-electrostatic energy coming from the i-j pair
2305 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2307 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2309 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2310 & 'eelloc',i,j,eel_loc_ij
2313 if (shield_mode.eq.0) then
2320 eel_loc_ij=eel_loc_ij
2321 & *fac_shield(i)*fac_shield(j)*sss1
2322 eel_loc=eel_loc+eel_loc_ij
2324 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2325 & (shield_mode.gt.0)) then
2328 do ilist=1,ishield_list(i)
2329 iresshield=shield_list(ilist,i)
2331 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2334 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2336 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2337 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2341 do ilist=1,ishield_list(j)
2342 iresshield=shield_list(ilist,j)
2344 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2347 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2349 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2350 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2357 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2358 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2359 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2360 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2361 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2362 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2363 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2364 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2369 geel_loc_ij=(a22*gmuij1(1)
2373 & *fac_shield(i)*fac_shield(j)*sss1
2374 c write(iout,*) "derivative over thatai"
2375 c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
2377 gloc(nphi+i,icg)=gloc(nphi+i,icg)+
2378 & geel_loc_ij*wel_loc
2379 c write(iout,*) "derivative over thatai-1"
2380 c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
2387 gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
2388 & geel_loc_ij*wel_loc
2389 & *fac_shield(i)*fac_shield(j)*sss1
2391 c Derivative over j residue
2392 geel_loc_ji=a22*gmuji1(1)
2396 c write(iout,*) "derivative over thataj"
2397 c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
2400 gloc(nphi+j,icg)=gloc(nphi+j,icg)+
2401 & geel_loc_ji*wel_loc
2402 & *fac_shield(i)*fac_shield(j)*sss1
2409 c write(iout,*) "derivative over thataj-1"
2410 c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
2412 gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
2413 & geel_loc_ji*wel_loc
2414 & *fac_shield(i)*fac_shield(j)*sss1
2416 cC Paral derivatives in virtual-bond dihedral angles gamma
2418 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2419 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2420 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2421 & *fac_shield(i)*fac_shield(j)*sss1
2422 c & *fac_shield(i)*fac_shield(j)
2423 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2426 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2427 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2428 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2429 & *fac_shield(i)*fac_shield(j)*sss1
2430 c & *fac_shield(i)*fac_shield(j)
2431 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2433 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2434 aux=eel_loc_ij/sss1*sssgrad1*rmij
2439 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+
2440 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2441 & *fac_shield(i)*fac_shield(j)*sss1
2442 c & *fac_shield(i)*fac_shield(j)
2443 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2445 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2446 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2447 cgrad ghalf=0.5d0*ggg(l)
2448 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
2449 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
2453 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2456 C Remaining derivatives of eello
2457 c gel_loc_long(3,j)=gel_loc_long(3,j)+ &
2458 c ssgradlipj*eel_loc_ij/2.0d0*lipscale/ &
2459 c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2461 c gel_loc_long(3,i)=gel_loc_long(3,i)+ &
2462 c ssgradlipi*eel_loc_ij/2.0d0*lipscale/ &
2463 c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2466 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2467 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2468 & *fac_shield(i)*fac_shield(j)*sss1
2469 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2471 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2472 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2473 & *fac_shield(i)*fac_shield(j)*sss1
2474 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2476 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2477 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2478 & *fac_shield(i)*fac_shield(j)*sss1
2479 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2481 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2482 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2483 & *fac_shield(i)*fac_shield(j)*sss1
2484 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2489 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2490 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
2491 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2492 & .and. num_conti.le.maxconts) then
2493 c write (iout,*) i,j," entered corr"
2495 C Calculate the contact function. The ith column of the array JCONT will
2496 C contain the numbers of atoms that make contacts with the atom I (of numbers
2497 C greater than I). The arrays FACONT and GACONT will contain the values of
2498 C the contact function and its derivative.
2499 c r0ij=1.02D0*rpp(iteli,itelj)
2500 c r0ij=1.11D0*rpp(iteli,itelj)
2501 r0ij=2.20D0*rpp(iteli,itelj)
2502 c r0ij=1.55D0*rpp(iteli,itelj)
2503 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2504 if (fcont.gt.0.0D0) then
2505 num_conti=num_conti+1
2506 if (num_conti.gt.maxconts) then
2507 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2508 & ' will skip next contacts for this conf.'
2510 jcont_hb(num_conti,i)=j
2511 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
2512 cd " jcont_hb",jcont_hb(num_conti,i)
2513 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2514 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2515 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2517 d_cont(num_conti,i)=rij
2518 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2519 C --- Electrostatic-interaction matrix ---
2520 a_chuj(1,1,num_conti,i)=a22
2521 a_chuj(1,2,num_conti,i)=a23
2522 a_chuj(2,1,num_conti,i)=a32
2523 a_chuj(2,2,num_conti,i)=a33
2524 C --- Gradient of rij
2526 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2533 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2534 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2535 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2536 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2537 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2542 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2543 C Calculate contact energies
2545 wij=cosa-3.0D0*cosb*cosg
2548 c fac3=dsqrt(-ael6i)/r0ij**3
2549 fac3=dsqrt(-ael6i)*r3ij
2550 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2551 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2552 if (ees0tmp.gt.0) then
2553 ees0pij=dsqrt(ees0tmp)
2557 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2558 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2559 if (ees0tmp.gt.0) then
2560 ees0mij=dsqrt(ees0tmp)
2565 if (shield_mode.eq.0) then
2569 ees0plist(num_conti,i)=j
2570 C fac_shield(i)=0.4d0
2571 C fac_shield(j)=0.6d0
2573 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2574 & *fac_shield(i)*fac_shield(j)*sss1
2575 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2576 & *fac_shield(i)*fac_shield(j)*sss1
2578 C Diagnostics. Comment out or remove after debugging!
2579 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2580 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2581 c ees0m(num_conti,i)=0.0D0
2583 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2584 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2585 C Angular derivatives of the contact function
2586 ees0pij1=fac3/ees0pij
2587 ees0mij1=fac3/ees0mij
2588 fac3p=-3.0D0*fac3*rrmij
2589 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2590 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2592 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2593 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2594 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2595 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2596 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2597 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2598 ecosap=ecosa1+ecosa2
2599 ecosbp=ecosb1+ecosb2
2600 ecosgp=ecosg1+ecosg2
2601 ecosam=ecosa1-ecosa2
2602 ecosbm=ecosb1-ecosb2
2603 ecosgm=ecosg1-ecosg2
2612 facont_hb(num_conti,i)=fcont
2613 fprimcont=fprimcont/rij
2614 cd facont_hb(num_conti,i)=1.0D0
2615 C Following line is for diagnostics.
2618 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2619 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2622 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2623 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2625 gggp(1)=gggp(1)+ees0pijp*xj
2626 & +ees0p(num_conti,i)/sss1*rmij*xj*sssgrad1
2627 gggp(2)=gggp(2)+ees0pijp*yj
2628 & +ees0p(num_conti,i)/sss1*rmij*yj*sssgrad1
2629 gggp(3)=gggp(3)+ees0pijp*zj
2630 & +ees0p(num_conti,i)/sss1*rmij*zj*sssgrad1
2631 gggm(1)=gggm(1)+ees0mijp*xj
2632 & +ees0m(num_conti,i)/sss1*rmij*xj*sssgrad1
2633 gggm(2)=gggm(2)+ees0mijp*yj
2634 & +ees0m(num_conti,i)/sss1*rmij*yj*sssgrad1
2635 gggm(3)=gggm(3)+ees0mijp*zj
2636 & +ees0m(num_conti,i)/sss1*rmij*zj*sssgrad1
2637 C Derivatives due to the contact function
2638 gacont_hbr(1,num_conti,i)=fprimcont*xj
2639 gacont_hbr(2,num_conti,i)=fprimcont*yj
2640 gacont_hbr(3,num_conti,i)=fprimcont*zj
2643 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2644 c following the change of gradient-summation algorithm.
2646 cgrad ghalfp=0.5D0*gggp(k)
2647 cgrad ghalfm=0.5D0*gggm(k)
2648 gacontp_hb1(k,num_conti,i)=!ghalfp
2649 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2650 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2651 & *sss1*fac_shield(i)*fac_shield(j)
2653 gacontp_hb2(k,num_conti,i)=!ghalfp
2654 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2655 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2656 & *sss1*fac_shield(i)*fac_shield(j)
2658 gacontp_hb3(k,num_conti,i)=gggp(k)
2659 & *sss1*fac_shield(i)*fac_shield(j)
2661 gacontm_hb1(k,num_conti,i)=!ghalfm
2662 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2663 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2664 & *sss1*fac_shield(i)*fac_shield(j)
2666 gacontm_hb2(k,num_conti,i)=!ghalfm
2667 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2668 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2669 & *sss1*fac_shield(i)*fac_shield(j)
2671 gacontm_hb3(k,num_conti,i)=gggm(k)
2672 & *sss1*fac_shield(i)*fac_shield(j)
2676 endif ! num_conti.le.maxconts
2680 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2683 ghalf=0.5d0*agg(l,k)
2684 aggi(l,k)=aggi(l,k)+ghalf
2685 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2686 aggj(l,k)=aggj(l,k)+ghalf
2689 if (j.eq.nres-1 .and. i.lt.j-2) then
2692 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2697 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2700 C-----------------------------------------------------------------------
2701 subroutine evdwpp_short(evdw1)
2706 include 'DIMENSIONS'
2707 include 'COMMON.CONTROL'
2708 include 'COMMON.IOUNITS'
2709 include 'COMMON.GEO'
2710 include 'COMMON.VAR'
2711 include 'COMMON.LOCAL'
2712 include 'COMMON.CHAIN'
2713 include 'COMMON.DERIV'
2714 include 'COMMON.INTERACT'
2715 c include 'COMMON.CONTACTS'
2716 include 'COMMON.TORSION'
2717 include 'COMMON.VECTORS'
2718 include 'COMMON.FFIELD'
2719 include "COMMON.SPLITELE"
2720 double precision ggg(3)
2721 integer xshift,yshift,zshift
2722 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2724 double precision scal_el /1.0d0/
2726 double precision scal_el /0.5d0/
2728 c write (iout,*) "evdwpp_short"
2729 integer i,j,k,iteli,itelj,num_conti,ind,isubchap
2730 double precision dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb
2731 double precision xj,yj,zj,rij,rrmij,r3ij,r6ij,evdw1,
2732 & dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
2733 & dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw
2734 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
2735 & dist_temp, dist_init,sss_grad
2736 double precision sscale,sscagrad
2740 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2741 c & " iatel_e_vdw",iatel_e_vdw
2743 c do i=iatel_s_vdw,iatel_e_vdw
2744 do ikont=g_listpp_vdw_start,g_listpp_vdw_end
2745 i=newcontlistpp_vdwi(ikont)
2746 j=newcontlistpp_vdwj(ikont)
2747 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2751 dx_normi=dc_norm(1,i)
2752 dy_normi=dc_norm(2,i)
2753 dz_normi=dc_norm(3,i)
2754 xmedi=c(1,i)+0.5d0*dxi
2755 ymedi=c(2,i)+0.5d0*dyi
2756 zmedi=c(3,i)+0.5d0*dzi
2757 xmedi=mod(xmedi,boxxsize)
2758 if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
2759 ymedi=mod(ymedi,boxysize)
2760 if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
2761 zmedi=mod(zmedi,boxzsize)
2762 if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
2764 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2765 c & ' ielend',ielend_vdw(i)
2767 c do j=ielstart_vdw(i),ielend_vdw(i)
2768 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2772 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2773 aaa=app(iteli,itelj)
2774 bbb=bpp(iteli,itelj)
2778 dx_normj=dc_norm(1,j)
2779 dy_normj=dc_norm(2,j)
2780 dz_normj=dc_norm(3,j)
2785 if (xj.lt.0) xj=xj+boxxsize
2787 if (yj.lt.0) yj=yj+boxysize
2789 if (zj.lt.0) zj=zj+boxzsize
2790 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2798 xj=xj_safe+xshift*boxxsize
2799 yj=yj_safe+yshift*boxysize
2800 zj=zj_safe+zshift*boxzsize
2801 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
2802 if(dist_temp.lt.dist_init) then
2812 if (isubchap.eq.1) then
2821 rij=xj*xj+yj*yj+zj*zj
2824 c sss=sscale(rij/rpp(iteli,itelj))
2825 c sssgrad=sscagrad(rij/rpp(iteli,itelj))
2826 sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
2827 sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
2828 if (sss.gt.0.0d0) then
2833 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2834 if (j.eq.i+2) ev1=scal_el*ev1
2837 if (energy_dec) then
2838 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2840 evdw1=evdw1+evdwij*sss
2841 if (energy_dec) write (iout,'(a10,2i5,0pf7.3)')
2842 & 'evdw1_sum',i,j,evdw1
2844 C Calculate contributions to the Cartesian gradient.
2846 facvdw=-6*rrmij*(ev1+evdwij)*sss
2847 ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
2848 ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
2849 ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
2854 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2855 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2862 C-----------------------------------------------------------------------------
2863 subroutine escp_long(evdw2,evdw2_14)
2865 C This subroutine calculates the excluded-volume interaction energy between
2866 C peptide-group centers and side chains and its gradient in virtual-bond and
2867 C side-chain vectors.
2869 implicit real*8 (a-h,o-z)
2870 include 'DIMENSIONS'
2871 include 'COMMON.GEO'
2872 include 'COMMON.VAR'
2873 include 'COMMON.LOCAL'
2874 include 'COMMON.CHAIN'
2875 include 'COMMON.DERIV'
2876 include 'COMMON.INTERACT'
2877 include 'COMMON.FFIELD'
2878 include 'COMMON.IOUNITS'
2879 include 'COMMON.CONTROL'
2880 include "COMMON.SPLITELE"
2881 logical lprint_short
2882 common /shortcheck/ lprint_short
2883 double precision ggg(3)
2884 integer xshift,yshift,zshift
2885 integer i,iint,j,k,iteli,itypj,subchap
2886 double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
2888 double precision evdw2,evdw2_14,evdwij
2889 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
2890 & dist_temp, dist_init
2891 double precision sscale,sscagrad
2893 if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
2896 CD print '(a)','Enter ESCP KURWA'
2897 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2899 c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
2900 c & ' iatscp_e=',iatscp_e
2901 c do i=iatscp_s,iatscp_e
2902 do ikont=g_listscp_start,g_listscp_end
2903 i=newcontlistscpi(ikont)
2904 j=newcontlistscpj(ikont)
2905 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2907 xi=0.5D0*(c(1,i)+c(1,i+1))
2908 yi=0.5D0*(c(2,i)+c(2,i+1))
2909 zi=0.5D0*(c(3,i)+c(3,i+1))
2911 if (xi.lt.0) xi=xi+boxxsize
2913 if (yi.lt.0) yi=yi+boxysize
2915 if (zi.lt.0) zi=zi+boxzsize
2917 c do iint=1,nscp_gr(i)
2919 c do j=iscpstart(i,iint),iscpend(i,iint)
2920 itypj=iabs(itype(j))
2921 if (itypj.eq.ntyp1) cycle
2922 C Uncomment following three lines for SC-p interactions
2926 C Uncomment following three lines for Ca-p interactions
2932 if (xj.lt.0) xj=xj+boxxsize
2934 if (yj.lt.0) yj=yj+boxysize
2936 if (zj.lt.0) zj=zj+boxzsize
2938 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2946 xj=xj_safe+xshift*boxxsize
2947 yj=yj_safe+yshift*boxysize
2948 zj=zj_safe+zshift*boxzsize
2949 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
2950 if(dist_temp.lt.dist_init) then
2960 if (subchap.eq.1) then
2970 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2972 sss1=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
2973 if (sss1.eq.0) cycle
2974 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
2976 & sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
2977 sssgrad1=sscagrad(1.0d0/dsqrt(rrij),r_cut_int)
2978 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2979 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2980 if (sss.lt.1.0d0) then
2982 e1=fac*fac*aad(itypj,iteli)
2983 e2=fac*bad(itypj,iteli)
2984 if (iabs(j-i) .le. 2) then
2987 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)*sss1
2990 evdw2=evdw2+evdwij*(1.0d0-sss)*sss1
2991 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2992 & 'evdw2',i,j,sss,evdwij
2994 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2997 fac=-(evdwij+e1)*rrij*(1.0d0-sss)*sss1
2998 fac=fac+evdwij*dsqrt(rrij)*(-sssgrad/rscp(itypj,iteli)
3003 C Uncomment following three lines for SC-p interactions
3005 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3007 C Uncomment following line for SC-p interactions
3008 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3010 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
3011 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
3020 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
3021 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
3022 gradx_scp(j,i)=expon*gradx_scp(j,i)
3025 C******************************************************************************
3029 C To save time the factor EXPON has been extracted from ALL components
3030 C of GVDWC and GRADX. Remember to multiply them by this factor before further
3033 C******************************************************************************
3036 C-----------------------------------------------------------------------------
3037 subroutine escp_short(evdw2,evdw2_14)
3039 C This subroutine calculates the excluded-volume interaction energy between
3040 C peptide-group centers and side chains and its gradient in virtual-bond and
3041 C side-chain vectors.
3044 include 'DIMENSIONS'
3045 include 'COMMON.GEO'
3046 include 'COMMON.VAR'
3047 include 'COMMON.LOCAL'
3048 include 'COMMON.CHAIN'
3049 include 'COMMON.DERIV'
3050 include 'COMMON.INTERACT'
3051 include 'COMMON.FFIELD'
3052 include 'COMMON.IOUNITS'
3053 include 'COMMON.CONTROL'
3054 include "COMMON.SPLITELE"
3055 integer xshift,yshift,zshift
3056 logical lprint_short
3057 common /shortcheck/ lprint_short
3058 integer i,iint,j,k,iteli,itypj,subchap
3059 double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
3061 double precision evdw2,evdw2_14,evdwij
3062 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
3063 & dist_temp, dist_init
3064 double precision ggg(3)
3065 double precision sscale,sscagrad
3068 cd print '(a)','Enter ESCP'
3070 c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
3071 c & ' iatscp_e=',iatscp_e
3072 if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
3073 do i=iatscp_s,iatscp_e
3074 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
3076 xi=0.5D0*(c(1,i)+c(1,i+1))
3077 yi=0.5D0*(c(2,i)+c(2,i+1))
3078 zi=0.5D0*(c(3,i)+c(3,i+1))
3080 if (xi.lt.0) xi=xi+boxxsize
3082 if (yi.lt.0) yi=yi+boxysize
3084 if (zi.lt.0) zi=zi+boxzsize
3087 c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
3088 c & " nscp_gr",nscp_gr(i)
3089 do iint=1,nscp_gr(i)
3091 do j=iscpstart(i,iint),iscpend(i,iint)
3092 itypj=iabs(itype(j))
3094 c & write (iout,*) "j",j," itypj",itypj
3095 if (itypj.eq.ntyp1) cycle
3096 C Uncomment following three lines for SC-p interactions
3100 C Uncomment following three lines for Ca-p interactions
3106 if (xj.lt.0) xj=xj+boxxsize
3108 if (yj.lt.0) yj=yj+boxysize
3110 if (zj.lt.0) zj=zj+boxzsize
3112 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
3113 c if (lprint_short) then
3114 c write (iout,*) i,j,xi,yi,zi,xj,yj,zj
3115 c write (iout,*) "dist_init",dsqrt(dist_init)
3124 xj=xj_safe+xshift*boxxsize
3125 yj=yj_safe+yshift*boxysize
3126 zj=zj_safe+zshift*boxzsize
3127 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
3128 if(dist_temp.lt.dist_init) then
3138 c if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp)
3139 if (subchap.eq.1) then
3148 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
3149 c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
3150 c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
3151 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
3152 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),
3154 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
3155 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
3156 c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),
3157 c & " subchap",subchap," sss",sss
3158 if (sss.gt.0.0d0) then
3161 e1=fac*fac*aad(itypj,iteli)
3162 e2=fac*bad(itypj,iteli)
3163 if (iabs(j-i) .le. 2) then
3166 evdw2_14=evdw2_14+(e1+e2)*sss
3169 evdw2=evdw2+evdwij*sss
3170 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
3171 & 'evdw2',i,j,sss,evdwij
3173 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
3175 fac=-(evdwij+e1)*rrij*sss
3176 fac=fac+evdwij*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)/expon
3180 C Uncomment following three lines for SC-p interactions
3182 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3184 C Uncomment following line for SC-p interactions
3185 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3187 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
3188 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
3197 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
3198 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
3199 gradx_scp(j,i)=expon*gradx_scp(j,i)
3202 C******************************************************************************
3206 C To save time the factor EXPON has been extracted from ALL components
3207 C of GVDWC and GRADX. Remember to multiply them by this factor before further
3210 C******************************************************************************