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 fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
89 double precision sscale,sscagrad,sscagradlip,sscalelip
90 double precision boxshift
91 double precision gg_lipi(3),gg_lipj(3)
92 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
96 c do i=iatsc_s,iatsc_e
97 do ikont=g_listscsc_start,g_listscsc_end
101 if (itypi.eq.ntyp1) cycle
102 itypi1=iabs(itype(i+1))
106 call to_box(xi,yi,zi)
107 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
109 C Calculate SC interaction energy.
111 c do iint=1,nint_gr(i)
112 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
113 cd & 'iend=',iend(i,iint)
114 c do j=istart(i,iint),iend(i,iint)
116 if (itypj.eq.ntyp1) cycle
120 call to_box(xj,yj,zj)
121 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
122 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
123 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
124 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
125 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
126 xj=boxshift(xj-xi,boxxsize)
127 yj=boxshift(yj-yi,boxysize)
128 zj=boxshift(zj-zi,boxzsize)
129 rij=xj*xj+yj*yj+zj*zj
131 eps0ij=eps(itypi,itypj)
132 sss1=sscale(sqrij,r_cut_int)
133 if (sss1.eq.0.0d0) cycle
134 sssgrad1=sscagrad(sqrij,r_cut_int)
136 & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
137 sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
138 if (sss.lt.1.0d0) then
145 evdw=evdw+(1.0d0-sss)*sss1*evdwij/sqrij/expon
147 C Calculate the components of the gradient in DC and X
149 fac=-rrij*(e1+evdwij)*(1.0d0-sss)*sss1
150 & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
151 & +(1.0d0-sss)*sssgrad1)/sqrij
155 gg_lipi(3)=(sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
156 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
157 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon
158 gg_lipj(3)=ssgradlipj*gg_lipi(3)
159 gg_lipi(3)=gg_lipi(3)*ssgradlipi
161 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
162 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
163 gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k)
164 gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k)
172 gvdwc(j,i)=expon*gvdwc(j,i)
173 gvdwx(j,i)=expon*gvdwx(j,i)
176 C******************************************************************************
180 C To save time, the factor of EXPON has been extracted from ALL components
181 C of GVDWC and GRADX. Remember to multiply them by this factor before further
184 C******************************************************************************
187 C-----------------------------------------------------------------------
188 subroutine elj_short(evdw)
190 C This subroutine calculates the interaction energy of nonbonded side chains
191 C assuming the LJ potential of interaction.
197 include 'COMMON.LOCAL'
198 include 'COMMON.CHAIN'
199 include 'COMMON.DERIV'
200 include 'COMMON.INTERACT'
201 include 'COMMON.TORSION'
202 include 'COMMON.SBRIDGE'
203 include 'COMMON.NAMES'
204 include 'COMMON.IOUNITS'
205 include "COMMON.SPLITELE"
206 c include 'COMMON.CONTACTS'
207 double precision gg(3)
208 double precision evdw,evdwij
209 integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont
210 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
211 & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
212 double precision sscale,sscagrad,sscagradlip,sscalelip
213 double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
215 double precision boxshift
216 double precision gg_lipi(3),gg_lipj(3)
217 c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
221 c do i=iatsc_s,iatsc_e
222 do ikont=g_listscsc_start,g_listscsc_end
223 i=newcontlisti(ikont)
224 j=newcontlistj(ikont)
226 if (itypi.eq.ntyp1) cycle
227 itypi1=iabs(itype(i+1))
231 call to_box(xi,yi,zi)
232 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
236 C Calculate SC interaction energy.
238 c do iint=1,nint_gr(i)
239 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
240 cd & 'iend=',iend(i,iint)
241 c do j=istart(i,iint),iend(i,iint)
243 if (itypj.eq.ntyp1) cycle
247 call to_box(xj,yj,zj)
248 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
249 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
250 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
251 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
252 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
253 xj=boxshift(xj-xi,boxxsize)
254 yj=boxshift(yj-yi,boxysize)
255 zj=boxshift(zj-zi,boxzsize)
256 C Change 12/1/95 to calculate four-body interactions
257 rij=xj*xj+yj*yj+zj*zj
259 sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
260 if (sss.gt.0.0d0) then
262 & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
264 eps0ij=eps(itypi,itypj)
272 C Calculate the components of the gradient in DC and X
274 fac=-rrij*(e1+evdwij)*sss+evdwij*sssgrad/sqrij/expon
278 gg_lipi(3)=(sss/2.0d0*(faclip*faclip*
279 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
280 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon
281 gg_lipj(3)=ssgradlipj*gg_lipi(3)
282 gg_lipi(3)=gg_lipi(3)*ssgradlipi
284 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
285 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
286 gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k)
287 gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k)
295 gvdwc(j,i)=expon*gvdwc(j,i)
296 gvdwx(j,i)=expon*gvdwx(j,i)
299 C******************************************************************************
303 C To save time, the factor of EXPON has been extracted from ALL components
304 C of GVDWC and GRADX. Remember to multiply them by this factor before further
307 C******************************************************************************
310 C-----------------------------------------------------------------------------
311 subroutine eljk_long(evdw)
313 C This subroutine calculates the interaction energy of nonbonded side chains
314 C assuming the LJK potential of interaction.
320 include 'COMMON.LOCAL'
321 include 'COMMON.CHAIN'
322 include 'COMMON.DERIV'
323 include 'COMMON.INTERACT'
324 include 'COMMON.IOUNITS'
325 include 'COMMON.NAMES'
326 include "COMMON.SPLITELE"
327 double precision gg(3)
328 double precision evdw,evdwij
329 integer i,j,k,itypi,itypj,itypi1,iint,ikont
330 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
331 & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
333 double precision sscale,sscagrad,sscagradlip,sscalelip
334 double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
336 double precision boxshift
337 double precision gg_lipi(3),gg_lipj(3)
338 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
342 c do i=iatsc_s,iatsc_e
343 do ikont=g_listscsc_start,g_listscsc_end
344 i=newcontlisti(ikont)
345 j=newcontlistj(ikont)
347 if (itypi.eq.ntyp1) cycle
348 itypi1=iabs(itype(i+1))
352 call to_box(xi,yi,zi)
353 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
355 C Calculate SC interaction energy.
357 c do iint=1,nint_gr(i)
358 c do j=istart(i,iint),iend(i,iint)
360 if (itypj.eq.ntyp1) cycle
364 call to_box(xj,yj,zj)
365 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
366 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
367 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
368 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
369 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
370 xj=boxshift(xj-xi,boxxsize)
371 yj=boxshift(yj-yi,boxysize)
372 zj=boxshift(zj-zi,boxzsize)
373 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
375 e_augm=augm(itypi,itypj)*fac_augm
378 sss1=sscale(rij,r_cut_int)
379 if (sss1.eq.0.0d0) cycle
380 sssgrad1=sscagrad(rij,r_cut_int)
381 sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
382 if (sss.lt.1.0d0) then
384 & sscagrad(rij/sigma(itypi,itypj),r_cut_respa)
385 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
386 fac=r_shift_inv**expon
391 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
392 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
393 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
394 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
395 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
396 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
397 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
398 evdw=evdw+(1.0d0-sss)*sss1*evdwij
400 C Calculate the components of the gradient in DC and X
402 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
403 fac=fac*(1.0d0-sss)*sss1
404 & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
405 & +(1.0d0-sss)*sssgrad1)*r_inv_ij/expon
409 gg_lipi(3)=(sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
410 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
411 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon
412 gg_lipj(3)=ssgradlipj*gg_lipi(3)
413 gg_lipi(3)=gg_lipi(3)*ssgradlipi
415 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
416 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
417 gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k)
418 gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k)
426 gvdwc(j,i)=expon*gvdwc(j,i)
427 gvdwx(j,i)=expon*gvdwx(j,i)
432 C-----------------------------------------------------------------------------
433 subroutine eljk_short(evdw)
435 C This subroutine calculates the interaction energy of nonbonded side chains
436 C assuming the LJK potential of interaction.
438 implicit real*8 (a-h,o-z)
442 include 'COMMON.LOCAL'
443 include 'COMMON.CHAIN'
444 include 'COMMON.DERIV'
445 include 'COMMON.INTERACT'
446 include 'COMMON.IOUNITS'
447 include 'COMMON.NAMES'
448 include "COMMON.SPLITELE"
449 double precision gg(3)
450 double precision evdw,evdwij
451 integer i,j,k,itypi,itypj,itypi1,iint,ikont
452 double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
453 & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
455 double precision sscale,sscagrad,sscagradlip,sscalelip
456 double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
458 double precision boxshift
459 double precision gg_lipi(3),gg_lipj(3)
460 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
464 c do i=iatsc_s,iatsc_e
465 do ikont=g_listscsc_start,g_listscsc_end
466 i=newcontlisti(ikont)
467 j=newcontlistj(ikont)
469 if (itypi.eq.ntyp1) cycle
470 itypi1=iabs(itype(i+1))
474 call to_box(xi,yi,zi)
475 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
477 C Calculate SC interaction energy.
479 c do iint=1,nint_gr(i)
480 c do j=istart(i,iint),iend(i,iint)
482 if (itypj.eq.ntyp1) cycle
486 call to_box(xj,yj,zj)
487 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
488 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
489 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
490 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
491 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
492 xj=boxshift(xj-xi,boxxsize)
493 yj=boxshift(yj-yi,boxysize)
494 zj=boxshift(zj-zi,boxzsize)
495 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
497 e_augm=augm(itypi,itypj)*fac_augm
500 sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
501 if (sss.gt.0.0d0) then
502 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
503 fac=r_shift_inv**expon
508 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
509 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
510 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
511 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
512 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
513 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
514 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
517 C Calculate the components of the gradient in DC and X
519 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
520 & +evdwij*sssgrad/sigma(itypi,itypj)*r_inv_ij/expon
525 gg_lipi(3)=(sss/2.0d0*(faclip*faclip*
526 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
527 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon
528 gg_lipj(3)=ssgradlipj*gg_lipi(3)
529 gg_lipi(3)=gg_lipi(3)*ssgradlipi
531 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
532 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
533 gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k)
534 gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k)
542 gvdwc(j,i)=expon*gvdwc(j,i)
543 gvdwx(j,i)=expon*gvdwx(j,i)
548 C-----------------------------------------------------------------------------
549 subroutine ebp_long(evdw)
551 C This subroutine calculates the interaction energy of nonbonded side chains
552 C assuming the Berne-Pechukas potential of interaction.
558 include 'COMMON.LOCAL'
559 include 'COMMON.CHAIN'
560 include 'COMMON.DERIV'
561 include 'COMMON.NAMES'
562 include 'COMMON.INTERACT'
563 include 'COMMON.IOUNITS'
564 include 'COMMON.CALC'
565 include "COMMON.SPLITELE"
568 double precision evdw
569 integer itypi,itypj,itypi1,iint,ind,ikont
570 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
571 double precision sss1,sssgrad1
572 double precision sscale,sscagrad,sscagradlip,sscalelip
573 double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
575 double precision boxshift
576 c double precision rrsave(maxdim)
579 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
582 c if (icall.eq.0) then
588 c do i=iatsc_s,iatsc_e
589 do ikont=g_listscsc_start,g_listscsc_end
590 i=newcontlisti(ikont)
591 j=newcontlistj(ikont)
593 if (itypi.eq.ntyp1) cycle
594 itypi1=iabs(itype(i+1))
598 call to_box(xi,yi,zi)
599 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
600 dxi=dc_norm(1,nres+i)
601 dyi=dc_norm(2,nres+i)
602 dzi=dc_norm(3,nres+i)
603 c dsci_inv=dsc_inv(itypi)
604 dsci_inv=vbld_inv(i+nres)
606 C Calculate SC interaction energy.
608 c do iint=1,nint_gr(i)
609 c do j=istart(i,iint),iend(i,iint)
612 if (itypj.eq.ntyp1) cycle
613 c dscj_inv=dsc_inv(itypj)
614 dscj_inv=vbld_inv(j+nres)
615 chi1=chi(itypi,itypj)
616 chi2=chi(itypj,itypi)
623 alf12=0.5D0*(alf1+alf2)
627 call to_box(xj,yj,zj)
628 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
629 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
630 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
631 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
632 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
633 xj=boxshift(xj-xi,boxxsize)
634 yj=boxshift(yj-yi,boxysize)
635 zj=boxshift(zj-zi,boxzsize)
636 dxj=dc_norm(1,nres+j)
637 dyj=dc_norm(2,nres+j)
638 dzj=dc_norm(3,nres+j)
639 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
641 sss1=sscale(1.0d0/rij,r_cut_int)
642 if (sss1.eq.0.0d0) cycle
643 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
644 if (sss.lt.1.0d0) then
646 & sscagrad(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
647 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
648 C Calculate the angle-dependent terms of energy & contributions to derivatives.
650 C Calculate whole angle-dependent part of epsilon and contributions
652 fac=(rrij*sigsq)**expon2
656 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
657 eps2der=evdwij*eps3rt
658 eps3der=evdwij*eps2rt
659 evdwij=evdwij*eps2rt*eps3rt
660 evdw=evdw+evdwij*(1.0d0-sss)*sss1
662 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
664 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
665 cd & restyp(itypi),i,restyp(itypj),j,
666 cd & epsi,sigm,chi1,chi2,chip1,chip2,
667 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
668 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
671 C Calculate gradient components.
672 e1=e1*eps1*eps2rt**2*eps3rt**2
673 fac=-expon*(e1+evdwij)
675 fac=(fac+evdwij*(sss1/(1.0d0-sss)*sssgrad/
676 & sigmaii(itypi,itypj)+(1.0d0-sss)/sss1*sssgrad1))*rij
677 C Calculate radial part of the gradient
681 gg_lipi(3)=eps1*(eps2rt*eps2rt)
682 & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
683 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
684 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
685 gg_lipj(3)=ssgradlipj*gg_lipi(3)
686 gg_lipi(3)=gg_lipi(3)*ssgradlipi
687 C Calculate the angular part of the gradient and sum add the contributions
688 C to the appropriate components of the Cartesian gradient.
689 call sc_grad_scale((1.0d0-sss)*sss1)
697 C-----------------------------------------------------------------------------
698 subroutine ebp_short(evdw)
700 C This subroutine calculates the interaction energy of nonbonded side chains
701 C assuming the Berne-Pechukas potential of interaction.
703 implicit real*8 (a-h,o-z)
707 include 'COMMON.LOCAL'
708 include 'COMMON.CHAIN'
709 include 'COMMON.DERIV'
710 include 'COMMON.NAMES'
711 include 'COMMON.INTERACT'
712 include 'COMMON.IOUNITS'
713 include 'COMMON.CALC'
714 include "COMMON.SPLITELE"
717 double precision evdw
718 integer itypi,itypj,itypi1,iint,ind,ikont
719 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
720 double precision sscale,sscagrad,sscagradlip,sscalelip
721 double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
723 double precision boxshift
724 c double precision rrsave(maxdim)
729 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
730 c if (icall.eq.0) then
736 c do i=iatsc_s,iatsc_e
737 do ikont=g_listscsc_start,g_listscsc_end
738 i=newcontlisti(ikont)
739 j=newcontlistj(ikont)
741 if (itypi.eq.ntyp1) cycle
742 itypi1=iabs(itype(i+1))
746 call to_box(xi,yi,zi)
747 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
748 dxi=dc_norm(1,nres+i)
749 dyi=dc_norm(2,nres+i)
750 dzi=dc_norm(3,nres+i)
751 c dsci_inv=dsc_inv(itypi)
752 dsci_inv=vbld_inv(i+nres)
754 C Calculate SC interaction energy.
756 c do iint=1,nint_gr(i)
757 c do j=istart(i,iint),iend(i,iint)
760 if (itypj.eq.ntyp1) cycle
761 c dscj_inv=dsc_inv(itypj)
762 dscj_inv=vbld_inv(j+nres)
763 chi1=chi(itypi,itypj)
764 chi2=chi(itypj,itypi)
771 alf12=0.5D0*(alf1+alf2)
775 call to_box(xj,yj,zj)
776 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
777 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
778 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
779 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
780 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
781 xj=boxshift(xj-xi,boxxsize)
782 yj=boxshift(yj-yi,boxysize)
783 zj=boxshift(zj-zi,boxzsize)
784 dxj=dc_norm(1,nres+j)
785 dyj=dc_norm(2,nres+j)
786 dzj=dc_norm(3,nres+j)
787 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
789 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
791 if (sss.gt.0.0d0) then
793 C Calculate the angle-dependent terms of energy & contributions to derivatives.
795 C Calculate whole angle-dependent part of epsilon and contributions
797 fac=(rrij*sigsq)**expon2
801 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
802 eps2der=evdwij*eps3rt
803 eps3der=evdwij*eps2rt
804 evdwij=evdwij*eps2rt*eps3rt
807 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
809 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
810 cd & restyp(itypi),i,restyp(itypj),j,
811 cd & epsi,sigm,chi1,chi2,chip1,chip2,
812 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
813 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
816 C Calculate gradient components.
817 e1=e1*eps1*eps2rt**2*eps3rt**2
818 fac=-expon*(e1+evdwij)
820 fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rrij
821 C Calculate radial part of the gradient
825 gg_lipi(3)=(sss/2.0d0*(faclip*faclip*
826 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
827 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon
828 gg_lipj(3)=ssgradlipj*gg_lipi(3)
829 gg_lipi(3)=gg_lipi(3)*ssgradlipi
831 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
832 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
833 gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k)
834 gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k)
836 C Calculate the angular part of the gradient and sum add the contributions
837 C to the appropriate components of the Cartesian gradient.
838 call sc_grad_scale(sss)
846 C-----------------------------------------------------------------------------
847 subroutine egb_long(evdw)
849 C This subroutine calculates the interaction energy of nonbonded side chains
850 C assuming the Gay-Berne potential of interaction.
856 include 'COMMON.LOCAL'
857 include 'COMMON.CHAIN'
858 include 'COMMON.DERIV'
859 include 'COMMON.NAMES'
860 include 'COMMON.INTERACT'
861 include 'COMMON.IOUNITS'
862 include 'COMMON.CALC'
863 include 'COMMON.CONTROL'
864 include "COMMON.SPLITELE"
866 integer xshift,yshift,zshift
867 double precision evdw
868 integer itypi,itypj,itypi1,iint,ind,ikont
869 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
870 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
871 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
872 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
873 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
874 double precision subchap,sss1,sssgrad1
875 double precision boxshift
877 ccccc energy_dec=.false.
878 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
882 c if (icall.eq.0) lprn=.false.
884 c do i=iatsc_s,iatsc_e
885 do ikont=g_listscsc_start,g_listscsc_end
886 i=newcontlisti(ikont)
887 j=newcontlistj(ikont)
889 if (itypi.eq.ntyp1) cycle
890 itypi1=iabs(itype(i+1))
894 call to_box(xi,yi,zi)
895 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
896 dxi=dc_norm(1,nres+i)
897 dyi=dc_norm(2,nres+i)
898 dzi=dc_norm(3,nres+i)
899 c dsci_inv=dsc_inv(itypi)
900 dsci_inv=vbld_inv(i+nres)
901 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
902 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
904 C Calculate SC interaction energy.
906 c do iint=1,nint_gr(i)
907 c do j=istart(i,iint),iend(i,iint)
910 if (itypj.eq.ntyp1) cycle
911 c dscj_inv=dsc_inv(itypj)
912 dscj_inv=vbld_inv(j+nres)
913 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
914 c & 1.0d0/vbld(j+nres)
915 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
916 sig0ij=sigma(itypi,itypj)
917 chi1=chi(itypi,itypj)
918 chi2=chi(itypj,itypi)
925 alf12=0.5D0*(alf1+alf2)
929 call to_box(xj,yj,zj)
930 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
931 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
932 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
933 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
934 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
935 xj=boxshift(xj-xi,boxxsize)
936 yj=boxshift(yj-yi,boxysize)
937 zj=boxshift(zj-zi,boxzsize)
938 dxj=dc_norm(1,nres+j)
939 dyj=dc_norm(2,nres+j)
940 dzj=dc_norm(3,nres+j)
941 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
943 sss1=sscale(1.0d0/rij,r_cut_int)
944 if (sss1.eq.0.0d0) cycle
945 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
946 if (sss.lt.1.0d0) then
947 C Calculate angle-dependent terms of energy and contributions to their
950 & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
951 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
954 sig=sig0ij*dsqrt(sigsq)
955 rij_shift=1.0D0/rij-sig+sig0ij
956 c for diagnostics; uncomment
957 c rij_shift=1.2*sig0ij
958 C I hate to put IF's in the loops, but here don't have another choice!!!!
959 if (rij_shift.le.0.0D0) then
961 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
962 cd & restyp(itypi),i,restyp(itypj),j,
963 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
967 c---------------------------------------------------------------
968 rij_shift=1.0D0/rij_shift
973 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
974 eps2der=evdwij*eps3rt
975 eps3der=evdwij*eps2rt
976 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
977 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
978 evdwij=evdwij*eps2rt*eps3rt
979 evdw=evdw+evdwij*(1.0d0-sss)*sss1
981 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
983 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
984 & restyp(itypi),i,restyp(itypj),j,
985 & epsi,sigm,chi1,chi2,chip1,chip2,
986 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
987 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
991 if (energy_dec) write (iout,'(a,2i5,5f10.5,e15.5)')
992 & 'r sss evdw',i,j,1.0d0/rij,sss1,sss,sslipi,sslipj,evdwij
994 C Calculate gradient components.
995 e1=e1*eps1*eps2rt**2*eps3rt**2
996 fac=-expon*(e1+evdwij)*rij_shift
998 fac=(fac+evdwij*(-sss1*sssgrad/(1.0d0-sss)
999 & /sigmaii(itypi,itypj)+(1.0d0-sss)*sssgrad1/sss1))*rij
1001 C Calculate the radial part of the gradient
1005 gg_lipi(3)=eps1*(eps2rt*eps2rt)
1006 & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
1007 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
1008 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
1009 gg_lipj(3)=ssgradlipj*gg_lipi(3)
1010 gg_lipi(3)=gg_lipi(3)*ssgradlipi
1012 C Calculate angular part of the gradient.
1013 call sc_grad_scale((1.0d0-sss)*sss1)
1018 c write (iout,*) "Number of loop steps in EGB:",ind
1019 cccc energy_dec=.false.
1022 C-----------------------------------------------------------------------------
1023 subroutine egb_short(evdw)
1025 C This subroutine calculates the interaction energy of nonbonded side chains
1026 C assuming the Gay-Berne potential of interaction.
1029 include 'DIMENSIONS'
1030 include 'COMMON.GEO'
1031 include 'COMMON.VAR'
1032 include 'COMMON.LOCAL'
1033 include 'COMMON.CHAIN'
1034 include 'COMMON.DERIV'
1035 include 'COMMON.NAMES'
1036 include 'COMMON.INTERACT'
1037 include 'COMMON.IOUNITS'
1038 include 'COMMON.CALC'
1039 include 'COMMON.CONTROL'
1040 include "COMMON.SPLITELE"
1042 double precision evdw
1043 integer itypi,itypj,itypi1,iint,ind,ikont
1044 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
1045 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1046 & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
1047 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1048 double precision boxshift
1050 ccccc energy_dec=.false.
1051 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1055 c if (icall.eq.0) lprn=.false.
1057 c do i=iatsc_s,iatsc_e
1058 do ikont=g_listscsc_start,g_listscsc_end
1059 i=newcontlisti(ikont)
1060 j=newcontlistj(ikont)
1061 itypi=iabs(itype(i))
1062 if (itypi.eq.ntyp1) cycle
1063 itypi1=iabs(itype(i+1))
1067 call to_box(xi,yi,zi)
1068 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
1069 dxi=dc_norm(1,nres+i)
1070 dyi=dc_norm(2,nres+i)
1071 dzi=dc_norm(3,nres+i)
1072 c dsci_inv=dsc_inv(itypi)
1073 dsci_inv=vbld_inv(i+nres)
1074 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1075 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1077 C Calculate SC interaction energy.
1079 c do iint=1,nint_gr(i)
1080 c do j=istart(i,iint),iend(i,iint)
1082 itypj=iabs(itype(j))
1083 if (itypj.eq.ntyp1) cycle
1084 c dscj_inv=dsc_inv(itypj)
1085 dscj_inv=vbld_inv(j+nres)
1086 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1087 c & 1.0d0/vbld(j+nres)
1088 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1089 sig0ij=sigma(itypi,itypj)
1090 chi1=chi(itypi,itypj)
1091 chi2=chi(itypj,itypi)
1098 alf12=0.5D0*(alf1+alf2)
1102 call to_box(xj,yj,zj)
1103 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
1104 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1105 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1106 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1107 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1108 c write (iout,*) "aa bb",aa_lip(itypi,itypj),
1109 c & bb_lip(itypi,itypj),aa_aq(itypi,itypj),
1110 c & bb_aq(itypi,itypj),aa,bb
1111 c write (iout,*) (sslipi+sslipj)/2.0d0,
1112 c & (2.0d0-sslipi-sslipj)/2.0d0
1113 xj=boxshift(xj-xi,boxxsize)
1114 yj=boxshift(yj-yi,boxysize)
1115 zj=boxshift(zj-zi,boxzsize)
1116 dxj=dc_norm(1,nres+j)
1117 dyj=dc_norm(2,nres+j)
1118 dzj=dc_norm(3,nres+j)
1119 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1121 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1122 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1123 if (sss.gt.0.0d0) then
1125 C Calculate angle-dependent terms of energy and contributions to their
1129 sig=sig0ij*dsqrt(sigsq)
1130 rij_shift=1.0D0/rij-sig+sig0ij
1131 c for diagnostics; uncomment
1132 c rij_shift=1.2*sig0ij
1133 C I hate to put IF's in the loops, but here don't have another choice!!!!
1134 if (rij_shift.le.0.0D0) then
1136 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1137 cd & restyp(itypi),i,restyp(itypj),j,
1138 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1142 c---------------------------------------------------------------
1143 rij_shift=1.0D0/rij_shift
1144 fac=rij_shift**expon
1148 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1149 eps2der=evdwij*eps3rt
1150 eps3der=evdwij*eps2rt
1151 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1152 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1153 evdwij=evdwij*eps2rt*eps3rt
1154 evdw=evdw+evdwij*sss
1156 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1158 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1159 & restyp(itypi),i,restyp(itypj),j,
1160 & epsi,sigm,chi1,chi2,chip1,chip2,
1161 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1162 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1166 if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)')
1167 & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij
1169 C Calculate gradient components.
1170 e1=e1*eps1*eps2rt**2*eps3rt**2
1171 fac=-expon*(e1+evdwij)*rij_shift
1173 fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rij
1175 C Calculate the radial part of the gradient
1179 gg_lipi(3)=eps1*(eps2rt*eps2rt)
1180 & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
1181 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
1182 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
1183 gg_lipj(3)=ssgradlipj*gg_lipi(3)
1184 gg_lipi(3)=gg_lipi(3)*ssgradlipi
1185 c write (iout,*) "gglip",i,j,gg_lipi,gg_lipj
1186 C Calculate angular part of the gradient.
1187 call sc_grad_scale(sss)
1192 c write (iout,*) "Number of loop steps in EGB:",ind
1193 cccc energy_dec=.false.
1196 C-----------------------------------------------------------------------------
1197 subroutine egbv_long(evdw)
1199 C This subroutine calculates the interaction energy of nonbonded side chains
1200 C assuming the Gay-Berne-Vorobjev potential of interaction.
1202 implicit real*8 (a-h,o-z)
1203 include 'DIMENSIONS'
1204 include 'COMMON.GEO'
1205 include 'COMMON.VAR'
1206 include 'COMMON.LOCAL'
1207 include 'COMMON.CHAIN'
1208 include 'COMMON.DERIV'
1209 include 'COMMON.NAMES'
1210 include 'COMMON.INTERACT'
1211 include 'COMMON.IOUNITS'
1212 include 'COMMON.CALC'
1213 include "COMMON.SPLITELE"
1215 common /srutu/ icall
1217 integer itypi,itypj,itypi1,iint,ind,ikont
1218 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1219 & xi,yi,zi,fac_augm,e_augm
1220 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1221 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1222 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1223 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1224 double precision sss1,sssgrad1
1228 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1230 c if (icall.eq.0) lprn=.true.
1232 c do i=iatsc_s,iatsc_e
1233 do ikont=g_listscsc_start,g_listscsc_end
1234 i=newcontlisti(ikont)
1235 j=newcontlistj(ikont)
1236 itypi=iabs(itype(i))
1237 if (itypi.eq.ntyp1) cycle
1238 itypi1=iabs(itype(i+1))
1242 call to_box(xi,yi,zi)
1243 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
1244 dxi=dc_norm(1,nres+i)
1245 dyi=dc_norm(2,nres+i)
1246 dzi=dc_norm(3,nres+i)
1247 c dsci_inv=dsc_inv(itypi)
1248 dsci_inv=vbld_inv(i+nres)
1250 C Calculate SC interaction energy.
1252 c do iint=1,nint_gr(i)
1253 c do j=istart(i,iint),iend(i,iint)
1255 itypj=iabs(itype(j))
1256 if (itypj.eq.ntyp1) cycle
1257 c dscj_inv=dsc_inv(itypj)
1258 dscj_inv=vbld_inv(j+nres)
1259 sig0ij=sigma(itypi,itypj)
1260 r0ij=r0(itypi,itypj)
1261 chi1=chi(itypi,itypj)
1262 chi2=chi(itypj,itypi)
1269 alf12=0.5D0*(alf1+alf2)
1273 call to_box(xj,yj,zj)
1274 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
1275 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1276 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1277 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1278 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1279 xj=boxshift(xj-xi,boxxsize)
1280 yj=boxshift(yj-yi,boxysize)
1281 zj=boxshift(zj-zi,boxzsize)
1282 dxj=dc_norm(1,nres+j)
1283 dyj=dc_norm(2,nres+j)
1284 dzj=dc_norm(3,nres+j)
1285 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1288 sss1=sscale(1.0d0/rij,r_cut_int)
1289 if (sss1.eq.0.0d0) cycle
1291 if (sss.lt.1.0d0) then
1293 & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1294 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
1296 C Calculate angle-dependent terms of energy and contributions to their
1300 sig=sig0ij*dsqrt(sigsq)
1301 rij_shift=1.0D0/rij-sig+r0ij
1302 C I hate to put IF's in the loops, but here don't have another choice!!!!
1303 if (rij_shift.le.0.0D0) then
1308 c---------------------------------------------------------------
1309 rij_shift=1.0D0/rij_shift
1310 fac=rij_shift**expon
1314 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1315 eps2der=evdwij*eps3rt
1316 eps3der=evdwij*eps2rt
1317 fac_augm=rrij**expon
1318 e_augm=augm(itypi,itypj)*fac_augm
1319 evdwij=evdwij*eps2rt*eps3rt
1320 evdw=evdw+(evdwij+e_augm)*sss1*(1.0d0-sss)
1322 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1324 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1325 & restyp(itypi),i,restyp(itypj),j,
1326 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1327 & chi1,chi2,chip1,chip2,
1328 & eps1,eps2rt**2,eps3rt**2,
1329 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1332 C Calculate gradient components.
1333 e1=e1*eps1*eps2rt**2*eps3rt**2
1334 fac=-expon*(e1+evdwij)*rij_shift
1336 fac=rij*fac-2*expon*rrij*e_augm
1337 fac=fac+(evdwij+e_augm)*
1338 & (-sss1*sssgrad/(1.0d0-sss)/sigmaii(itypi,itypj)
1339 & +(1.0d0-sss)*sssgrad1/sss1)*rij
1340 C Calculate the radial part of the gradient
1344 gg_lipi(3)=eps1*(eps2rt*eps2rt)
1345 & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
1346 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
1347 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
1348 gg_lipj(3)=ssgradlipj*gg_lipi(3)
1349 gg_lipi(3)=gg_lipi(3)*ssgradlipi
1350 C Calculate angular part of the gradient.
1351 call sc_grad_scale((1.0d0-sss)*sss1)
1357 C-----------------------------------------------------------------------------
1358 subroutine egbv_short(evdw)
1360 C This subroutine calculates the interaction energy of nonbonded side chains
1361 C assuming the Gay-Berne-Vorobjev potential of interaction.
1364 include 'DIMENSIONS'
1365 include 'COMMON.GEO'
1366 include 'COMMON.VAR'
1367 include 'COMMON.LOCAL'
1368 include 'COMMON.CHAIN'
1369 include 'COMMON.DERIV'
1370 include 'COMMON.NAMES'
1371 include 'COMMON.INTERACT'
1372 include 'COMMON.IOUNITS'
1373 include 'COMMON.CALC'
1374 include "COMMON.SPLITELE"
1376 common /srutu/ icall
1378 integer itypi,itypj,itypi1,iint,ind,ikont
1379 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1380 & xi,yi,zi,fac_augm,e_augm
1381 double precision evdw
1382 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1383 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1384 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1385 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1386 double precision boxshift
1390 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1392 c if (icall.eq.0) lprn=.true.
1394 c do i=iatsc_s,iatsc_e
1395 do ikont=g_listscsc_start,g_listscsc_end
1396 i=newcontlisti(ikont)
1397 j=newcontlistj(ikont)
1398 itypi=iabs(itype(i))
1399 if (itypi.eq.ntyp1) cycle
1400 itypi1=iabs(itype(i+1))
1404 dxi=dc_norm(1,nres+i)
1405 dyi=dc_norm(2,nres+i)
1406 dzi=dc_norm(3,nres+i)
1407 call to_box(xi,yi,zi)
1408 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
1409 c dsci_inv=dsc_inv(itypi)
1410 dsci_inv=vbld_inv(i+nres)
1412 C Calculate SC interaction energy.
1414 c do iint=1,nint_gr(i)
1415 c do j=istart(i,iint),iend(i,iint)
1417 itypj=iabs(itype(j))
1418 if (itypj.eq.ntyp1) cycle
1419 c dscj_inv=dsc_inv(itypj)
1420 dscj_inv=vbld_inv(j+nres)
1421 sig0ij=sigma(itypi,itypj)
1422 r0ij=r0(itypi,itypj)
1423 chi1=chi(itypi,itypj)
1424 chi2=chi(itypj,itypi)
1431 alf12=0.5D0*(alf1+alf2)
1435 call to_box(xj,yj,zj)
1436 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
1437 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1438 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1439 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1440 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1441 xj=boxshift(xj-xi,boxxsize)
1442 yj=boxshift(yj-yi,boxysize)
1443 zj=boxshift(zj-zi,boxzsize)
1444 dxj=dc_norm(1,nres+j)
1445 dyj=dc_norm(2,nres+j)
1446 dzj=dc_norm(3,nres+j)
1447 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1450 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1452 if (sss.gt.0.0d0) then
1454 C Calculate angle-dependent terms of energy and contributions to their
1458 sig=sig0ij*dsqrt(sigsq)
1459 rij_shift=1.0D0/rij-sig+r0ij
1460 C I hate to put IF's in the loops, but here don't have another choice!!!!
1461 if (rij_shift.le.0.0D0) then
1466 c---------------------------------------------------------------
1467 rij_shift=1.0D0/rij_shift
1468 fac=rij_shift**expon
1472 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1473 eps2der=evdwij*eps3rt
1474 eps3der=evdwij*eps2rt
1475 fac_augm=rrij**expon
1476 e_augm=augm(itypi,itypj)*fac_augm
1477 evdwij=evdwij*eps2rt*eps3rt
1478 evdw=evdw+(evdwij+e_augm)*sss
1480 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1482 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1483 & restyp(itypi),i,restyp(itypj),j,
1484 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1485 & chi1,chi2,chip1,chip2,
1486 & eps1,eps2rt**2,eps3rt**2,
1487 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1490 C Calculate gradient components.
1491 e1=e1*eps1*eps2rt**2*eps3rt**2
1492 fac=-expon*(e1+evdwij)*rij_shift
1494 fac=rij*fac-2*expon*rrij*e_augm+
1495 & (evdwij+e_augm)*sssgrad/sigmaii(itypi,itypj)/sss*rij
1496 C Calculate the radial part of the gradient
1500 gg_lipi(3)=eps1*(eps2rt*eps2rt)
1501 & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
1502 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
1503 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
1504 gg_lipj(3)=ssgradlipj*gg_lipi(3)
1505 gg_lipi(3)=gg_lipi(3)*ssgradlipi
1506 C Calculate angular part of the gradient.
1507 call sc_grad_scale(sss)
1513 C----------------------------------------------------------------------------
1514 subroutine sc_grad_scale(scalfac)
1515 implicit real*8 (a-h,o-z)
1516 include 'DIMENSIONS'
1517 include 'COMMON.CHAIN'
1518 include 'COMMON.DERIV'
1519 include 'COMMON.CALC'
1520 include 'COMMON.IOUNITS'
1521 include "COMMON.SPLITELE"
1522 double precision dcosom1(3),dcosom2(3)
1523 double precision scalfac
1524 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1525 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1526 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1527 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1531 c eom12=evdwij*eps1_om12
1533 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1534 c & " sigder",sigder
1535 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1536 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1538 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1539 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1542 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1544 c write (iout,*) "gg",(gg(k),k=1,3)
1546 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
1547 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1548 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1549 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
1550 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1551 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1552 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1553 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1554 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1555 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1558 C Calculate the components of the gradient in DC and X
1560 c write (iout,*) "scgrad gglip",i,j,gg_lipi,gg_lipj
1562 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
1563 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
1567 C--------------------------------------------------------------------------
1568 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1570 C This subroutine calculates the average interaction energy and its gradient
1571 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1572 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1573 C The potential depends both on the distance of peptide-group centers and on
1574 C the orientation of the CA-CA virtual bonds.
1576 implicit real*8 (a-h,o-z)
1580 include 'DIMENSIONS'
1581 include 'COMMON.CONTROL'
1582 include 'COMMON.SETUP'
1583 include 'COMMON.IOUNITS'
1584 include 'COMMON.GEO'
1585 include 'COMMON.VAR'
1586 include 'COMMON.LOCAL'
1587 include 'COMMON.CHAIN'
1588 include 'COMMON.DERIV'
1589 include 'COMMON.INTERACT'
1591 include 'COMMON.CONTACTS'
1592 include 'COMMON.CONTMAT'
1594 include 'COMMON.CORRMAT'
1595 include 'COMMON.TORSION'
1596 include 'COMMON.VECTORS'
1597 include 'COMMON.FFIELD'
1598 include 'COMMON.TIME1'
1599 include 'COMMON.SHIELD'
1600 include "COMMON.SPLITELE"
1602 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1603 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1604 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1605 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1606 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1607 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1609 double precision sslipi,sslipj,ssgradlipi,ssgradlipj
1610 common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj
1611 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1613 double precision scal_el /1.0d0/
1615 double precision scal_el /0.5d0/
1618 C 13-go grudnia roku pamietnego...
1619 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1620 & 0.0d0,1.0d0,0.0d0,
1621 & 0.0d0,0.0d0,1.0d0/
1622 cd write(iout,*) 'In EELEC'
1624 cd write(iout,*) 'Type',i
1625 cd write(iout,*) 'B1',B1(:,i)
1626 cd write(iout,*) 'B2',B2(:,i)
1627 cd write(iout,*) 'CC',CC(:,:,i)
1628 cd write(iout,*) 'DD',DD(:,:,i)
1629 cd write(iout,*) 'EE',EE(:,:,i)
1631 cd call check_vecgrad
1633 C print *,"WCHODZE3"
1634 if (icheckgrad.eq.1) then
1636 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1638 dc_norm(k,i)=dc(k,i)*fac
1640 c write (iout,*) 'i',i,' fac',fac
1643 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1644 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1645 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1646 c call vec_and_deriv
1652 time_mat=time_mat+MPI_Wtime()-time01
1656 cd write (iout,*) 'i=',i
1658 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1661 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1662 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1677 cd print '(a)','Enter EELEC'
1678 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1680 gel_loc_loc(i)=0.0d0
1685 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1687 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1689 do i=iturn3_start,iturn3_end
1690 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1691 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1692 C & .or. itype(i-1).eq.ntyp1
1693 C & .or. itype(i+4).eq.ntyp1
1698 dx_normi=dc_norm(1,i)
1699 dy_normi=dc_norm(2,i)
1700 dz_normi=dc_norm(3,i)
1701 xmedi=c(1,i)+0.5d0*dxi
1702 ymedi=c(2,i)+0.5d0*dyi
1703 zmedi=c(3,i)+0.5d0*dzi
1704 call to_box(xmedi,ymedi,zmedi)
1705 call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
1707 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1708 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1710 num_cont_hb(i)=num_conti
1713 do i=iturn4_start,iturn4_end
1714 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1715 & .or. itype(i+3).eq.ntyp1
1716 & .or. itype(i+4).eq.ntyp1
1717 C & .or. itype(i+5).eq.ntyp1
1718 C & .or. itype(i-1).eq.ntyp1
1723 dx_normi=dc_norm(1,i)
1724 dy_normi=dc_norm(2,i)
1725 dz_normi=dc_norm(3,i)
1726 xmedi=c(1,i)+0.5d0*dxi
1727 ymedi=c(2,i)+0.5d0*dyi
1728 zmedi=c(3,i)+0.5d0*dzi
1729 call to_box(xmedi,ymedi,zmedi)
1730 call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
1732 num_conti=num_cont_hb(i)
1734 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1735 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1736 & call eturn4(i,eello_turn4)
1738 num_cont_hb(i)=num_conti
1742 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1744 c do i=iatel_s,iatel_e
1745 do ikont=g_listpp_start,g_listpp_end
1746 i=newcontlistppi(ikont)
1747 j=newcontlistppj(ikont)
1748 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1749 C & .or. itype(i+2).eq.ntyp1
1750 C & .or. itype(i-1).eq.ntyp1
1755 dx_normi=dc_norm(1,i)
1756 dy_normi=dc_norm(2,i)
1757 dz_normi=dc_norm(3,i)
1758 xmedi=c(1,i)+0.5d0*dxi
1759 ymedi=c(2,i)+0.5d0*dyi
1760 zmedi=c(3,i)+0.5d0*dzi
1761 call to_box(xmedi,ymedi,zmedi)
1762 call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
1763 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
1765 num_conti=num_cont_hb(i)
1767 c do j=ielstart(i),ielend(i)
1768 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1769 C & .or.itype(j+2).eq.ntyp1
1770 C & .or.itype(j-1).eq.ntyp1
1772 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1775 num_cont_hb(i)=num_conti
1778 c write (iout,*) "Number of loop steps in EELEC:",ind
1780 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1781 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1783 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1784 ccc eel_loc=eel_loc+eello_turn3
1785 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1788 C-------------------------------------------------------------------------------
1789 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1791 include 'DIMENSIONS'
1795 include 'COMMON.CONTROL'
1796 include 'COMMON.IOUNITS'
1797 include 'COMMON.GEO'
1798 include 'COMMON.VAR'
1799 include 'COMMON.LOCAL'
1800 include 'COMMON.CHAIN'
1801 include 'COMMON.DERIV'
1802 include 'COMMON.INTERACT'
1804 include 'COMMON.CONTACTS'
1805 include 'COMMON.CONTMAT'
1807 include 'COMMON.CORRMAT'
1808 include 'COMMON.TORSION'
1809 include 'COMMON.VECTORS'
1810 include 'COMMON.FFIELD'
1811 include 'COMMON.TIME1'
1812 include 'COMMON.SHIELD'
1813 include "COMMON.SPLITELE"
1814 integer xshift,yshift,zshift
1815 double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1816 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1817 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1818 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
1819 & gmuij2(4),gmuji2(4)
1820 integer j1,j2,num_conti
1821 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1822 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1824 integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ind,itypi,itypj
1825 integer ilist,iresshield
1826 double precision rlocshield
1827 double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
1828 double precision ees,evdw1,eel_loc,aaa,bbb,ael3i
1829 double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,
1830 & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,
1831 & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa,
1832 & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,
1833 & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,
1834 & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,
1835 & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
1836 & ecosgp,ecosam,ecosbm,ecosgm,ghalf,geel_loc_ij,geel_loc_ji,
1837 & dxi,dyi,dzi,a22,a23,a32,a33
1838 double precision dist_init,xmedi,ymedi,zmedi,xj_safe,yj_safe,
1839 & zj_safe,xj_temp,yj_temp,zj_temp,dist_temp,dx_normi,dy_normi,
1841 double precision sss1,sssgrad1
1842 double precision sscale,sscagrad
1843 double precision scalar
1844 double precision boxshift
1845 double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij,
1847 common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij
1848 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1850 double precision scal_el /1.0d0/
1852 double precision scal_el /0.5d0/
1855 C 13-go grudnia roku pamietnego...
1856 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1857 & 0.0d0,1.0d0,0.0d0,
1858 & 0.0d0,0.0d0,1.0d0/
1859 c time00=MPI_Wtime()
1860 cd write (iout,*) "eelecij",i,j
1861 C print *,"WCHODZE2"
1865 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1866 aaa=app(iteli,itelj)
1867 bbb=bpp(iteli,itelj)
1868 ael6i=ael6(iteli,itelj)
1869 ael3i=ael3(iteli,itelj)
1873 dx_normj=dc_norm(1,j)
1874 dy_normj=dc_norm(2,j)
1875 dz_normj=dc_norm(3,j)
1879 call to_box(xj,yj,zj)
1880 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
1881 faclipij=(sslipi+sslipj)/2.0d0*lipscale+1.0d0
1882 faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
1883 xj=boxshift(xj-xmedi,boxxsize)
1884 yj=boxshift(yj-ymedi,boxysize)
1885 zj=boxshift(zj-zmedi,boxzsize)
1886 rij=xj*xj+yj*yj+zj*zj
1890 c For extracting the short-range part of Evdwpp
1891 sss1=sscale(rij,r_cut_int)
1892 if (sss1.eq.0.0d0) return
1893 sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
1894 sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
1895 sssgrad1=sscagrad(rij,r_cut_int)
1898 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1899 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1900 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1901 fac=cosa-3.0D0*cosb*cosg
1903 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1904 if (j.eq.i+2) ev1=scal_el*ev1
1909 if (shield_mode.eq.0) then
1913 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1915 el1=el1*fac_shield(i)**2*fac_shield(j)**2
1916 el2=el2*fac_shield(i)**2*fac_shield(j)**2
1918 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1919 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1920 ees=ees+eesij*sss1*faclipij2
1921 evdw1=evdw1+evdwij*(1.0d0-sss)*sss1*faclipij2
1922 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1923 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1924 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1925 cd & xmedi,ymedi,zmedi,xj,yj,zj
1927 if (energy_dec) then
1928 write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,5f10.5)')
1929 & 'evdw1',i,j,evdwij,iteli,itelj,aaa,sss,sss1,sssgrad,sssgrad1,rij
1930 write (iout,'(a6,2i5,0pf7.3,6f8.5)') 'ees',i,j,eesij,
1931 & fac_shield(i),fac_shield(j),sslipi,sslipj,faclipij,faclipij2
1935 C Calculate contributions to the Cartesian gradient.
1938 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
1939 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1940 facel=-3*rrmij*(el1+eesij)*sss1
1941 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1947 * Radial derivatives. First process both termini of the fragment (i,j)
1949 c old aux=(facel+sssgrad1*(1.0d0-sss)*eesij*rmij)*faclipij2
1950 aux=(facel+sssgrad1*eesij*rmij)*faclipij2
1951 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1958 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1959 & (shield_mode.gt.0)) then
1961 do ilist=1,ishield_list(i)
1962 iresshield=shield_list(ilist,i)
1964 rlocshield=grad_shield_side(k,ilist,i)*eesij*sss1
1965 & /fac_shield(i)*2.0*sss1
1966 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1968 & +grad_shield_loc(k,ilist,i)*eesij*sss1/fac_shield(i)*2.0
1970 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1971 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1972 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1973 C if (iresshield.gt.i) then
1974 C do ishi=i+1,iresshield-1
1975 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1976 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1980 C do ishi=iresshield,i
1981 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1982 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1988 do ilist=1,ishield_list(j)
1989 iresshield=shield_list(ilist,j)
1991 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1993 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1995 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss1
1996 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
2001 gshieldc(k,i)=gshieldc(k,i)+
2002 & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
2003 gshieldc(k,j)=gshieldc(k,j)+
2004 & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
2005 gshieldc(k,i-1)=gshieldc(k,i-1)+
2006 & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
2007 gshieldc(k,j-1)=gshieldc(k,j-1)+
2008 & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
2014 c ghalf=0.5D0*ggg(k)
2015 c gelc(k,i)=gelc(k,i)+ghalf
2016 c gelc(k,j)=gelc(k,j)+ghalf
2018 c 9/28/08 AL Gradient compotents will be summed only at the end
2020 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
2021 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
2023 gelc_long(3,j)=gelc_long(3,j)+
2024 & ssgradlipj*eesij/2.0d0*lipscale**2*sss1
2025 gelc_long(3,i)=gelc_long(3,i)+
2026 & ssgradlipi*eesij/2.0d0*lipscale**2*sss1
2027 c gelc_long(3,i)=gelc_long(3,i)+
2028 c ssgradlipi*eesij/2.0d0*lipscale**2*sss1
2030 * Loop over residues i+1 thru j-1.
2034 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
2038 &(-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij)
2044 c ghalf=0.5D0*ggg(k)
2045 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2046 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2048 c 9/28/08 AL Gradient compotents will be summed only at the end
2050 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2051 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2053 !C Lipidic part for scaling weight
2054 gvdwpp(3,j)=gvdwpp(3,j)+
2055 & sss1*(1.0d0-sss)*ssgradlipj*evdwij/2.0d0*lipscale**2
2056 gvdwpp(3,i)=gvdwpp(3,i)+
2057 & sss1*(1.0d0-sss)*ssgradlipi*evdwij/2.0d0*lipscale**2
2059 * Loop over residues i+1 thru j-1.
2063 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2067 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
2068 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2069 facel=-3*rrmij*(el1+eesij)*sss1
2070 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2072 c facvdw=ev1+evdwij*(1.0d0-sss)*sss1
2075 fac=-3*rrmij*(facvdw+facvdw+facel)
2080 * Radial derivatives. First process both termini of the fragment (i,j)
2082 aux=(fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj))
2083 & *eesij*rmij)*faclipij2
2084 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2092 c ghalf=0.5D0*ggg(k)
2093 c gelc(k,i)=gelc(k,i)+ghalf
2094 c gelc(k,j)=gelc(k,j)+ghalf
2096 c 9/28/08 AL Gradient compotents will be summed only at the end
2098 gelc_long(k,j)=gelc(k,j)+ggg(k)
2099 gelc_long(k,i)=gelc(k,i)-ggg(k)
2102 * Loop over residues i+1 thru j-1.
2106 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
2109 c 9/28/08 AL Gradient compotents will be summed only at the end
2114 & (-sssgrad*sss1/rpp(iteli,itelj)+sssgrad1*(1.0d0-sss))*rmij*evdwij
2119 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2120 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2122 gvdwpp(3,j)=gvdwpp(3,j)+
2123 & sss1*ssgradlipj*evdwij/2.0d0*lipscale**2
2124 gvdwpp(3,i)=gvdwpp(3,i)+
2125 & sss1*ssgradlipi*evdwij/2.0d0*lipscale**2
2130 ecosa=2.0D0*fac3*fac1+fac4
2133 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
2134 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
2136 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2137 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2139 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
2140 cd & (dcosg(k),k=1,3)
2142 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1
2143 & *fac_shield(i)**2*fac_shield(j)**2*faclipij2
2144 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2148 c ghalf=0.5D0*ggg(k)
2149 c gelc(k,i)=gelc(k,i)+ghalf
2150 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
2151 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2152 c gelc(k,j)=gelc(k,j)+ghalf
2153 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2154 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2158 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
2163 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
2164 & +ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss1
2165 & *fac_shield(i)**2*fac_shield(j)**2*faclipij2
2166 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2169 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2170 & +ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss1
2171 & *fac_shield(i)**2*fac_shield(j)**2*faclipij2
2172 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2173 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
2174 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
2176 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
2177 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
2178 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2180 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
2181 C energy of a peptide unit is assumed in the form of a second-order
2182 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
2183 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
2184 C are computed for EVERY pair of non-contiguous peptide groups.
2186 if (j.lt.nres-1) then
2197 muij(kkk)=mu(k,i)*mu(l,j)
2199 gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
2200 c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
2201 gmuij2(kkk)=gUb2(k,i)*mu(l,j)
2202 gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
2203 c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
2204 gmuji2(kkk)=mu(k,i)*gUb2(l,j)
2208 cd write (iout,*) 'EELEC: i',i,' j',j
2209 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
2210 cd write(iout,*) 'muij',muij
2211 ury=scalar(uy(1,i),erij)
2212 urz=scalar(uz(1,i),erij)
2213 vry=scalar(uy(1,j),erij)
2214 vrz=scalar(uz(1,j),erij)
2215 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
2216 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
2217 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
2218 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
2219 fac=dsqrt(-ael6i)*r3ij
2224 cd write (iout,'(4i5,4f10.5)')
2225 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
2226 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
2227 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
2228 cd & uy(:,j),uz(:,j)
2229 cd write (iout,'(4f10.5)')
2230 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
2231 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
2232 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
2233 cd write (iout,'(9f10.5/)')
2234 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
2235 C Derivatives of the elements of A in virtual-bond vectors
2236 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
2238 uryg(k,1)=scalar(erder(1,k),uy(1,i))
2239 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
2240 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
2241 urzg(k,1)=scalar(erder(1,k),uz(1,i))
2242 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
2243 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
2244 vryg(k,1)=scalar(erder(1,k),uy(1,j))
2245 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
2246 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
2247 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
2248 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
2249 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
2251 C Compute radial contributions to the gradient
2269 C Add the contributions coming from er
2272 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
2273 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
2274 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
2275 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
2278 C Derivatives in DC(i)
2279 cgrad ghalf1=0.5d0*agg(k,1)
2280 cgrad ghalf2=0.5d0*agg(k,2)
2281 cgrad ghalf3=0.5d0*agg(k,3)
2282 cgrad ghalf4=0.5d0*agg(k,4)
2283 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2284 & -3.0d0*uryg(k,2)*vry)!+ghalf1
2285 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2286 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
2287 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2288 & -3.0d0*urzg(k,2)*vry)!+ghalf3
2289 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2290 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
2291 C Derivatives in DC(i+1)
2292 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2293 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2294 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2295 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2296 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2297 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2298 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2299 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2300 C Derivatives in DC(j)
2301 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2302 & -3.0d0*vryg(k,2)*ury)!+ghalf1
2303 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2304 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
2305 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2306 & -3.0d0*vryg(k,2)*urz)!+ghalf3
2307 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2308 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
2309 C Derivatives in DC(j+1) or DC(nres-1)
2310 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2311 & -3.0d0*vryg(k,3)*ury)
2312 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2313 & -3.0d0*vrzg(k,3)*ury)
2314 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2315 & -3.0d0*vryg(k,3)*urz)
2316 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2317 & -3.0d0*vrzg(k,3)*urz)
2318 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
2320 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
2333 aggi(k,l)=-aggi(k,l)
2334 aggi1(k,l)=-aggi1(k,l)
2335 aggj(k,l)=-aggj(k,l)
2336 aggj1(k,l)=-aggj1(k,l)
2339 if (j.lt.nres-1) then
2345 aggi(k,l)=-aggi(k,l)
2346 aggi1(k,l)=-aggi1(k,l)
2347 aggj(k,l)=-aggj(k,l)
2348 aggj1(k,l)=-aggj1(k,l)
2359 aggi(k,l)=-aggi(k,l)
2360 aggi1(k,l)=-aggi1(k,l)
2361 aggj(k,l)=-aggj(k,l)
2362 aggj1(k,l)=-aggj1(k,l)
2367 IF (wel_loc.gt.0.0d0) THEN
2368 C Contribution to the local-electrostatic energy coming from the i-j pair
2369 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2371 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2373 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2374 & 'eelloc',i,j,eel_loc_ij
2377 if (shield_mode.eq.0) then
2384 eel_loc_ij=eel_loc_ij
2385 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2386 eel_loc=eel_loc+eel_loc_ij
2388 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2389 & (shield_mode.gt.0)) then
2392 do ilist=1,ishield_list(i)
2393 iresshield=shield_list(ilist,i)
2395 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2398 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2400 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2401 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2405 do ilist=1,ishield_list(j)
2406 iresshield=shield_list(ilist,j)
2408 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2411 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2413 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2414 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2421 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2422 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2423 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2424 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2425 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2426 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2427 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2428 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2433 geel_loc_ij=(a22*gmuij1(1)
2437 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2438 c write(iout,*) "derivative over thatai"
2439 c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
2441 gloc(nphi+i,icg)=gloc(nphi+i,icg)+
2442 & geel_loc_ij*wel_loc
2443 c write(iout,*) "derivative over thatai-1"
2444 c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
2451 gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
2452 & geel_loc_ij*wel_loc
2453 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2455 c Derivative over j residue
2456 geel_loc_ji=a22*gmuji1(1)
2460 c write(iout,*) "derivative over thataj"
2461 c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
2464 gloc(nphi+j,icg)=gloc(nphi+j,icg)+
2465 & geel_loc_ji*wel_loc
2466 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2473 c write(iout,*) "derivative over thataj-1"
2474 c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
2476 gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
2477 & geel_loc_ji*wel_loc
2478 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2480 cC Paral derivatives in virtual-bond dihedral angles gamma
2482 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2483 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2484 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2485 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2486 c & *fac_shield(i)*fac_shield(j)
2487 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2490 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2491 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2492 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2493 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2494 c & *fac_shield(i)*fac_shield(j)
2495 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2497 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2498 aux=eel_loc_ij/sss1*sssgrad1*rmij
2503 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+
2504 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2505 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2506 c & *fac_shield(i)*fac_shield(j)
2507 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2509 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2510 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2511 cgrad ghalf=0.5d0*ggg(l)
2512 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
2513 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
2515 gel_loc_long(3,j)=gel_loc_long(3,j)+
2516 & ssgradlipj*eel_loc_ij/2.0d0*lipscale/faclipij
2517 gel_loc_long(3,i)=gel_loc_long(3,i)+
2518 & ssgradlipi*eel_loc_ij/2.0d0*lipscale/faclipij
2521 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2524 C Remaining derivatives of eello
2525 c gel_loc_long(3,j)=gel_loc_long(3,j)+ &
2526 c ssgradlipj*eel_loc_ij/2.0d0*lipscale/ &
2527 c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2529 c gel_loc_long(3,i)=gel_loc_long(3,i)+ &
2530 c ssgradlipi*eel_loc_ij/2.0d0*lipscale/ &
2531 c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2534 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2535 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2536 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2537 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2539 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2540 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2541 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2542 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2544 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2545 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2546 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2547 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2549 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2550 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2551 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2552 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2557 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2558 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
2559 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2560 & .and. num_conti.le.maxconts) then
2561 c write (iout,*) i,j," entered corr"
2563 C Calculate the contact function. The ith column of the array JCONT will
2564 C contain the numbers of atoms that make contacts with the atom I (of numbers
2565 C greater than I). The arrays FACONT and GACONT will contain the values of
2566 C the contact function and its derivative.
2567 c r0ij=1.02D0*rpp(iteli,itelj)
2568 c r0ij=1.11D0*rpp(iteli,itelj)
2569 r0ij=2.20D0*rpp(iteli,itelj)
2570 c r0ij=1.55D0*rpp(iteli,itelj)
2571 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2572 if (fcont.gt.0.0D0) then
2573 num_conti=num_conti+1
2574 if (num_conti.gt.maxconts) then
2575 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2576 & ' will skip next contacts for this conf.'
2578 jcont_hb(num_conti,i)=j
2579 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
2580 cd " jcont_hb",jcont_hb(num_conti,i)
2581 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2582 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2583 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2585 d_cont(num_conti,i)=rij
2586 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2587 C --- Electrostatic-interaction matrix ---
2588 a_chuj(1,1,num_conti,i)=a22
2589 a_chuj(1,2,num_conti,i)=a23
2590 a_chuj(2,1,num_conti,i)=a32
2591 a_chuj(2,2,num_conti,i)=a33
2592 C --- Gradient of rij
2594 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2601 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2602 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2603 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2604 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2605 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2610 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2611 C Calculate contact energies
2613 wij=cosa-3.0D0*cosb*cosg
2616 c fac3=dsqrt(-ael6i)/r0ij**3
2617 fac3=dsqrt(-ael6i)*r3ij
2618 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2619 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2620 if (ees0tmp.gt.0) then
2621 ees0pij=dsqrt(ees0tmp)
2625 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2626 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2627 if (ees0tmp.gt.0) then
2628 ees0mij=dsqrt(ees0tmp)
2633 if (shield_mode.eq.0) then
2637 ees0plist(num_conti,i)=j
2638 C fac_shield(i)=0.4d0
2639 C fac_shield(j)=0.6d0
2641 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2642 & *fac_shield(i)*fac_shield(j)*sss1
2643 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2644 & *fac_shield(i)*fac_shield(j)*sss1
2646 C Diagnostics. Comment out or remove after debugging!
2647 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2648 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2649 c ees0m(num_conti,i)=0.0D0
2651 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2652 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2653 C Angular derivatives of the contact function
2654 ees0pij1=fac3/ees0pij
2655 ees0mij1=fac3/ees0mij
2656 fac3p=-3.0D0*fac3*rrmij
2657 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2658 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2660 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2661 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2662 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2663 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2664 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2665 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2666 ecosap=ecosa1+ecosa2
2667 ecosbp=ecosb1+ecosb2
2668 ecosgp=ecosg1+ecosg2
2669 ecosam=ecosa1-ecosa2
2670 ecosbm=ecosb1-ecosb2
2671 ecosgm=ecosg1-ecosg2
2680 facont_hb(num_conti,i)=fcont
2681 fprimcont=fprimcont/rij
2682 cd facont_hb(num_conti,i)=1.0D0
2683 C Following line is for diagnostics.
2686 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2687 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2690 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2691 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2693 gggp(1)=gggp(1)+ees0pijp*xj
2694 & +ees0p(num_conti,i)/sss1*rmij*xj*sssgrad1
2695 gggp(2)=gggp(2)+ees0pijp*yj
2696 & +ees0p(num_conti,i)/sss1*rmij*yj*sssgrad1
2697 gggp(3)=gggp(3)+ees0pijp*zj
2698 & +ees0p(num_conti,i)/sss1*rmij*zj*sssgrad1
2699 gggm(1)=gggm(1)+ees0mijp*xj
2700 & +ees0m(num_conti,i)/sss1*rmij*xj*sssgrad1
2701 gggm(2)=gggm(2)+ees0mijp*yj
2702 & +ees0m(num_conti,i)/sss1*rmij*yj*sssgrad1
2703 gggm(3)=gggm(3)+ees0mijp*zj
2704 & +ees0m(num_conti,i)/sss1*rmij*zj*sssgrad1
2705 C Derivatives due to the contact function
2706 gacont_hbr(1,num_conti,i)=fprimcont*xj
2707 gacont_hbr(2,num_conti,i)=fprimcont*yj
2708 gacont_hbr(3,num_conti,i)=fprimcont*zj
2711 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2712 c following the change of gradient-summation algorithm.
2714 cgrad ghalfp=0.5D0*gggp(k)
2715 cgrad ghalfm=0.5D0*gggm(k)
2716 gacontp_hb1(k,num_conti,i)=!ghalfp
2717 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2718 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2719 & *sss1*fac_shield(i)*fac_shield(j)
2721 gacontp_hb2(k,num_conti,i)=!ghalfp
2722 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2723 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2724 & *sss1*fac_shield(i)*fac_shield(j)
2726 gacontp_hb3(k,num_conti,i)=gggp(k)
2727 & *sss1*fac_shield(i)*fac_shield(j)
2729 gacontm_hb1(k,num_conti,i)=!ghalfm
2730 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2731 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2732 & *sss1*fac_shield(i)*fac_shield(j)
2734 gacontm_hb2(k,num_conti,i)=!ghalfm
2735 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2736 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2737 & *sss1*fac_shield(i)*fac_shield(j)
2739 gacontm_hb3(k,num_conti,i)=gggm(k)
2740 & *sss1*fac_shield(i)*fac_shield(j)
2744 endif ! num_conti.le.maxconts
2748 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2751 ghalf=0.5d0*agg(l,k)
2752 aggi(l,k)=aggi(l,k)+ghalf
2753 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2754 aggj(l,k)=aggj(l,k)+ghalf
2757 if (j.eq.nres-1 .and. i.lt.j-2) then
2760 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2765 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2768 C-----------------------------------------------------------------------
2769 subroutine evdwpp_short(evdw1)
2774 include 'DIMENSIONS'
2775 include 'COMMON.CONTROL'
2776 include 'COMMON.IOUNITS'
2777 include 'COMMON.GEO'
2778 include 'COMMON.VAR'
2779 include 'COMMON.LOCAL'
2780 include 'COMMON.CHAIN'
2781 include 'COMMON.DERIV'
2782 include 'COMMON.INTERACT'
2783 c include 'COMMON.CONTACTS'
2784 include 'COMMON.TORSION'
2785 include 'COMMON.VECTORS'
2786 include 'COMMON.FFIELD'
2787 include "COMMON.SPLITELE"
2788 double precision ggg(3)
2789 integer xshift,yshift,zshift
2790 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2792 double precision scal_el /1.0d0/
2794 double precision scal_el /0.5d0/
2796 c write (iout,*) "evdwpp_short"
2797 integer i,j,k,iteli,itelj,num_conti,ind,isubchap
2798 double precision dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb
2799 double precision xj,yj,zj,rij,rrmij,r3ij,r6ij,evdw1,
2800 & dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
2801 & dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw
2802 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
2803 & dist_temp, dist_init,sss_grad
2804 double precision sscale,sscagrad
2805 double precision sslipi,ssgradlipi,sslipj,ssgradlipj
2806 double precision boxshift
2808 double precision faclipij2
2811 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2812 c & " iatel_e_vdw",iatel_e_vdw
2814 c do i=iatel_s_vdw,iatel_e_vdw
2815 do ikont=g_listpp_vdw_start,g_listpp_vdw_end
2816 i=newcontlistpp_vdwi(ikont)
2817 j=newcontlistpp_vdwj(ikont)
2818 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2822 dx_normi=dc_norm(1,i)
2823 dy_normi=dc_norm(2,i)
2824 dz_normi=dc_norm(3,i)
2825 xmedi=c(1,i)+0.5d0*dxi
2826 ymedi=c(2,i)+0.5d0*dyi
2827 zmedi=c(3,i)+0.5d0*dzi
2828 call to_box(xmedi,ymedi,zmedi)
2829 call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
2831 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2832 c & ' ielend',ielend_vdw(i)
2834 c do j=ielstart_vdw(i),ielend_vdw(i)
2835 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2839 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2840 aaa=app(iteli,itelj)
2841 bbb=bpp(iteli,itelj)
2845 dx_normj=dc_norm(1,j)
2846 dy_normj=dc_norm(2,j)
2847 dz_normj=dc_norm(3,j)
2851 call to_box(xj,yj,zj)
2852 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
2853 faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
2854 xj=boxshift(xj-xmedi,boxxsize)
2855 yj=boxshift(yj-ymedi,boxysize)
2856 zj=boxshift(zj-zmedi,boxzsize)
2857 rij=xj*xj+yj*yj+zj*zj
2860 c sss=sscale(rij/rpp(iteli,itelj))
2861 c sssgrad=sscagrad(rij/rpp(iteli,itelj))
2862 sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
2863 sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
2864 if (sss.gt.0.0d0) then
2869 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2870 if (j.eq.i+2) ev1=scal_el*ev1
2873 if (energy_dec) then
2874 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2876 evdw1=evdw1+evdwij*sss*faclipij2
2877 if (energy_dec) write (iout,'(a10,2i5,0pf7.3)')
2878 & 'evdw1_sum',i,j,evdw1
2880 C Calculate contributions to the Cartesian gradient.
2882 facvdw=(-6*rrmij*(ev1+evdwij)*sss+sssgrad*rmij*evdwij/
2883 & rpp(iteli,itelj))*faclipij2
2891 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2892 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2894 !C Lipidic part for scaling weight
2895 gvdwpp(3,j)=gvdwpp(3,j)+
2896 & sss*ssgradlipj*evdwij/2.0d0*lipscale**2
2897 gvdwpp(3,i)=gvdwpp(3,i)+
2898 & sss*ssgradlipi*evdwij/2.0d0*lipscale**2
2904 C-----------------------------------------------------------------------------
2905 subroutine escp_long(evdw2,evdw2_14)
2907 C This subroutine calculates the excluded-volume interaction energy between
2908 C peptide-group centers and side chains and its gradient in virtual-bond and
2909 C side-chain vectors.
2911 implicit real*8 (a-h,o-z)
2912 include 'DIMENSIONS'
2913 include 'COMMON.GEO'
2914 include 'COMMON.VAR'
2915 include 'COMMON.LOCAL'
2916 include 'COMMON.CHAIN'
2917 include 'COMMON.DERIV'
2918 include 'COMMON.INTERACT'
2919 include 'COMMON.FFIELD'
2920 include 'COMMON.IOUNITS'
2921 include 'COMMON.CONTROL'
2922 include "COMMON.SPLITELE"
2923 logical lprint_short
2924 common /shortcheck/ lprint_short
2925 double precision ggg(3)
2926 integer i,iint,j,k,iteli,itypj,subchap
2927 double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
2929 double precision evdw2,evdw2_14,evdwij
2930 double precision sscale,sscagrad
2931 double precision boxshift
2933 if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
2936 CD print '(a)','Enter ESCP KURWA'
2937 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2939 c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
2940 c & ' iatscp_e=',iatscp_e
2941 c do i=iatscp_s,iatscp_e
2942 do ikont=g_listscp_start,g_listscp_end
2943 i=newcontlistscpi(ikont)
2944 j=newcontlistscpj(ikont)
2945 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2947 xi=0.5D0*(c(1,i)+c(1,i+1))
2948 yi=0.5D0*(c(2,i)+c(2,i+1))
2949 zi=0.5D0*(c(3,i)+c(3,i+1))
2950 call to_box(xi,yi,zi)
2952 c do iint=1,nscp_gr(i)
2954 c do j=iscpstart(i,iint),iscpend(i,iint)
2955 itypj=iabs(itype(j))
2956 if (itypj.eq.ntyp1) cycle
2957 C Uncomment following three lines for SC-p interactions
2961 C Uncomment following three lines for Ca-p interactions
2966 call to_box(xj,yj,zj)
2967 xj=boxshift(xj-xi,boxxsize)
2968 yj=boxshift(yj-yi,boxysize)
2969 zj=boxshift(zj-zi,boxzsize)
2971 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2973 sss1=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
2974 if (sss1.eq.0) cycle
2975 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
2977 & sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
2978 sssgrad1=sscagrad(1.0d0/dsqrt(rrij),r_cut_int)
2979 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
2980 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
2981 if (sss.lt.1.0d0) then
2983 e1=fac*fac*aad(itypj,iteli)
2984 e2=fac*bad(itypj,iteli)
2985 if (iabs(j-i) .le. 2) then
2988 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)*sss1
2991 evdw2=evdw2+evdwij*(1.0d0-sss)*sss1
2992 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
2993 & 'evdw2',i,j,sss,evdwij
2995 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2998 fac=-(evdwij+e1)*rrij*(1.0d0-sss)*sss1
2999 fac=fac+evdwij*dsqrt(rrij)*(-sssgrad/rscp(itypj,iteli)
3004 C Uncomment following three lines for SC-p interactions
3006 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3008 C Uncomment following line for SC-p interactions
3009 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3011 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
3012 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
3021 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
3022 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
3023 gradx_scp(j,i)=expon*gradx_scp(j,i)
3026 C******************************************************************************
3030 C To save time the factor EXPON has been extracted from ALL components
3031 C of GVDWC and GRADX. Remember to multiply them by this factor before further
3034 C******************************************************************************
3037 C-----------------------------------------------------------------------------
3038 subroutine escp_short(evdw2,evdw2_14)
3040 C This subroutine calculates the excluded-volume interaction energy between
3041 C peptide-group centers and side chains and its gradient in virtual-bond and
3042 C side-chain vectors.
3045 include 'DIMENSIONS'
3046 include 'COMMON.GEO'
3047 include 'COMMON.VAR'
3048 include 'COMMON.LOCAL'
3049 include 'COMMON.CHAIN'
3050 include 'COMMON.DERIV'
3051 include 'COMMON.INTERACT'
3052 include 'COMMON.FFIELD'
3053 include 'COMMON.IOUNITS'
3054 include 'COMMON.CONTROL'
3055 include "COMMON.SPLITELE"
3056 integer xshift,yshift,zshift
3057 logical lprint_short
3058 common /shortcheck/ lprint_short
3059 integer i,iint,j,k,iteli,itypj,subchap
3060 double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
3062 double precision evdw2,evdw2_14,evdwij
3063 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
3064 & dist_temp, dist_init
3065 double precision ggg(3)
3066 double precision sscale,sscagrad
3067 double precision boxshift
3070 cd print '(a)','Enter ESCP'
3072 c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
3073 c & ' iatscp_e=',iatscp_e
3074 if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
3075 do i=iatscp_s,iatscp_e
3076 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
3078 xi=0.5D0*(c(1,i)+c(1,i+1))
3079 yi=0.5D0*(c(2,i)+c(2,i+1))
3080 zi=0.5D0*(c(3,i)+c(3,i+1))
3081 call to_box(xi,yi,zi)
3084 c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
3085 c & " nscp_gr",nscp_gr(i)
3086 do iint=1,nscp_gr(i)
3088 do j=iscpstart(i,iint),iscpend(i,iint)
3089 itypj=iabs(itype(j))
3091 c & write (iout,*) "j",j," itypj",itypj
3092 if (itypj.eq.ntyp1) cycle
3093 C Uncomment following three lines for SC-p interactions
3097 C Uncomment following three lines for Ca-p interactions
3105 call to_box(xj,yj,zj)
3106 xj=boxshift(xj-xi,boxxsize)
3107 yj=boxshift(yj-yi,boxysize)
3108 zj=boxshift(zj-zi,boxzsize)
3109 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
3110 c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
3111 c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
3112 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
3113 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),
3115 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
3116 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
3117 c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),
3118 c & " subchap",subchap," sss",sss
3119 if (sss.gt.0.0d0) then
3122 e1=fac*fac*aad(itypj,iteli)
3123 e2=fac*bad(itypj,iteli)
3124 if (iabs(j-i) .le. 2) then
3127 evdw2_14=evdw2_14+(e1+e2)*sss
3130 evdw2=evdw2+evdwij*sss
3131 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
3132 & 'evdw2',i,j,sss,evdwij
3134 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
3136 fac=-(evdwij+e1)*rrij*sss
3137 fac=fac+evdwij*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)/expon
3141 C Uncomment following three lines for SC-p interactions
3143 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3145 C Uncomment following line for SC-p interactions
3146 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3148 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
3149 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
3158 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
3159 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
3160 gradx_scp(j,i)=expon*gradx_scp(j,i)
3163 C******************************************************************************
3167 C To save time the factor EXPON has been extracted from ALL components
3168 C of GVDWC and GRADX. Remember to multiply them by this factor before further
3171 C******************************************************************************