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_long,g_listscsc_end_long
98 i=newcontlisti_long(ikont)
99 j=newcontlistj_long(ikont)
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_short,g_listscsc_end_short
223 i=newcontlisti_short(ikont)
224 j=newcontlistj_short(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_long,g_listscsc_end_long
344 i=newcontlisti_long(ikont)
345 j=newcontlistj_long(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_short,g_listscsc_end_short
466 i=newcontlisti_short(ikont)
467 j=newcontlistj_short(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_long,g_listscsc_end_long
590 i=newcontlisti_long(ikont)
591 j=newcontlistj_long(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_short,g_listscsc_end_short
738 i=newcontlisti_short(ikont)
739 j=newcontlistj_short(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,rij1
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
886 & write(2,*) "g_listscsc_start_long,g_listscsc_end_long",
887 & g_listscsc_start_long,g_listscsc_end_long
888 do ikont=g_listscsc_start_long,g_listscsc_end_long
889 i=newcontlisti_long(ikont)
890 j=newcontlistj_long(ikont)
892 if (itypi.eq.ntyp1) cycle
893 itypi1=iabs(itype(i+1))
897 call to_box(xi,yi,zi)
898 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
899 dxi=dc_norm(1,nres+i)
900 dyi=dc_norm(2,nres+i)
901 dzi=dc_norm(3,nres+i)
902 c dsci_inv=dsc_inv(itypi)
903 dsci_inv=vbld_inv(i+nres)
904 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
905 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
907 C Calculate SC interaction energy.
909 c do iint=1,nint_gr(i)
910 c do j=istart(i,iint),iend(i,iint)
913 if (itypj.eq.ntyp1) cycle
914 c dscj_inv=dsc_inv(itypj)
915 dscj_inv=vbld_inv(j+nres)
916 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
917 c & 1.0d0/vbld(j+nres)
918 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
919 sig0ij=sigma(itypi,itypj)
920 chi1=chi(itypi,itypj)
921 chi2=chi(itypj,itypi)
928 alf12=0.5D0*(alf1+alf2)
932 call to_box(xj,yj,zj)
933 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
934 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
935 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
936 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
937 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
938 xj=boxshift(xj-xi,boxxsize)
939 yj=boxshift(yj-yi,boxysize)
940 zj=boxshift(zj-zi,boxzsize)
941 dxj=dc_norm(1,nres+j)
942 dyj=dc_norm(2,nres+j)
943 dzj=dc_norm(3,nres+j)
944 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
947 c sss1=sscale(1.0d0/rij,r_cut_int)
948 sss1=sscale(rij1,r_cut_int)
949 if (sss1.eq.0.0d0) cycle
950 rij1=rij1/sigmaii(itypi,itypj)
951 sss=sscale(rij1,r_cut_respa)
952 c sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
953 if (sss.lt.1.0d0) then
954 C Calculate angle-dependent terms of energy and contributions to their
957 & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
958 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
961 sig=sig0ij*dsqrt(sigsq)
962 rij_shift=1.0D0/rij-sig+sig0ij
963 c for diagnostics; uncomment
964 c rij_shift=1.2*sig0ij
965 C I hate to put IF's in the loops, but here don't have another choice!!!!
966 if (rij_shift.le.0.0D0) then
968 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
969 cd & restyp(itypi),i,restyp(itypj),j,
970 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
974 c---------------------------------------------------------------
975 rij_shift=1.0D0/rij_shift
980 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
981 eps2der=evdwij*eps3rt
982 eps3der=evdwij*eps2rt
983 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
984 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
985 evdwij=evdwij*eps2rt*eps3rt
986 evdw=evdw+evdwij*(1.0d0-sss)*sss1
988 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
990 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
991 & restyp(itypi),i,restyp(itypj),j,
992 & epsi,sigm,chi1,chi2,chip1,chip2,
993 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
994 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
998 if (energy_dec) write (iout,'(a,2i5,5f10.5,e15.5)')
999 & 'r sss evdw',i,j,1.0d0/rij,sss1,sss,sslipi,sslipj,evdwij
1001 C Calculate gradient components.
1002 e1=e1*eps1*eps2rt**2*eps3rt**2
1003 fac=-expon*(e1+evdwij)*rij_shift
1005 fac=(fac+evdwij*(-sss1*sssgrad/(1.0d0-sss)
1006 & /sigmaii(itypi,itypj)+(1.0d0-sss)*sssgrad1/sss1))*rij
1008 C Calculate the radial part of the gradient
1012 gg_lipi(3)=eps1*(eps2rt*eps2rt)
1013 & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
1014 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
1015 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
1016 gg_lipj(3)=ssgradlipj*gg_lipi(3)
1017 gg_lipi(3)=gg_lipi(3)*ssgradlipi
1019 C Calculate angular part of the gradient.
1020 call sc_grad_scale((1.0d0-sss)*sss1)
1025 c write (iout,*) "Number of loop steps in EGB:",ind
1026 cccc energy_dec=.false.
1029 C-----------------------------------------------------------------------------
1030 subroutine egb_short(evdw)
1032 C This subroutine calculates the interaction energy of nonbonded side chains
1033 C assuming the Gay-Berne potential of interaction.
1037 include 'DIMENSIONS'
1038 include 'COMMON.GEO'
1039 include 'COMMON.VAR'
1040 include 'COMMON.LOCAL'
1041 include 'COMMON.CHAIN'
1042 include 'COMMON.DERIV'
1043 include 'COMMON.NAMES'
1044 include 'COMMON.INTERACT'
1045 include 'COMMON.IOUNITS'
1046 include 'COMMON.CALC'
1047 include 'COMMON.CONTROL'
1048 include "COMMON.SPLITELE"
1049 include 'COMMON.TIME1'
1051 double precision evdw
1052 integer itypi,itypj,itypi1,iint,ind,ikont
1053 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
1054 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1055 & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
1056 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1057 double precision boxshift
1058 double precision time01
1059 c time01=MPI_Wtime()
1061 ccccc energy_dec=.false.
1062 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1066 c if (icall.eq.0) lprn=.false.
1068 c do i=iatsc_s,iatsc_e
1070 & write(2,*) "g_listscsc_start_short,g_listscsc_end_short",
1071 & g_listscsc_start_short,g_listscsc_end_short
1072 do ikont=g_listscsc_start_short,g_listscsc_end_short
1073 i=newcontlisti_short(ikont)
1074 j=newcontlistj_short(ikont)
1075 itypi=iabs(itype(i))
1076 if (itypi.eq.ntyp1) cycle
1077 itypi1=iabs(itype(i+1))
1081 call to_box(xi,yi,zi)
1082 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
1083 dxi=dc_norm(1,nres+i)
1084 dyi=dc_norm(2,nres+i)
1085 dzi=dc_norm(3,nres+i)
1086 c dsci_inv=dsc_inv(itypi)
1087 dsci_inv=vbld_inv(i+nres)
1088 c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
1089 c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
1091 C Calculate SC interaction energy.
1093 c do iint=1,nint_gr(i)
1094 c do j=istart(i,iint),iend(i,iint)
1096 itypj=iabs(itype(j))
1097 if (itypj.eq.ntyp1) cycle
1098 c dscj_inv=dsc_inv(itypj)
1099 dscj_inv=vbld_inv(j+nres)
1100 c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
1101 c & 1.0d0/vbld(j+nres)
1102 c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
1103 sig0ij=sigma(itypi,itypj)
1104 chi1=chi(itypi,itypj)
1105 chi2=chi(itypj,itypi)
1112 alf12=0.5D0*(alf1+alf2)
1116 call to_box(xj,yj,zj)
1117 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
1118 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1119 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1120 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1121 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1122 c write (iout,*) "aa bb",aa_lip(itypi,itypj),
1123 c & bb_lip(itypi,itypj),aa_aq(itypi,itypj),
1124 c & bb_aq(itypi,itypj),aa,bb
1125 c write (iout,*) (sslipi+sslipj)/2.0d0,
1126 c & (2.0d0-sslipi-sslipj)/2.0d0
1127 xj=boxshift(xj-xi,boxxsize)
1128 yj=boxshift(yj-yi,boxysize)
1129 zj=boxshift(zj-zi,boxzsize)
1130 dxj=dc_norm(1,nres+j)
1131 dyj=dc_norm(2,nres+j)
1132 dzj=dc_norm(3,nres+j)
1133 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1135 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1136 if (sss.gt.0.0d0) then
1137 sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1139 C Calculate angle-dependent terms of energy and contributions to their
1143 sig=sig0ij*dsqrt(sigsq)
1144 rij_shift=1.0D0/rij-sig+sig0ij
1145 c for diagnostics; uncomment
1146 c rij_shift=1.2*sig0ij
1147 C I hate to put IF's in the loops, but here don't have another choice!!!!
1148 if (rij_shift.le.0.0D0) then
1150 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1151 cd & restyp(itypi),i,restyp(itypj),j,
1152 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
1156 c---------------------------------------------------------------
1157 rij_shift=1.0D0/rij_shift
1158 fac=rij_shift**expon
1162 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1163 eps2der=evdwij*eps3rt
1164 eps3der=evdwij*eps2rt
1165 c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
1166 c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
1167 evdwij=evdwij*eps2rt*eps3rt
1168 evdw=evdw+evdwij*sss
1170 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1172 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1173 & restyp(itypi),i,restyp(itypj),j,
1174 & epsi,sigm,chi1,chi2,chip1,chip2,
1175 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
1176 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1180 if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)')
1181 & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij
1183 C Calculate gradient components.
1184 e1=e1*eps1*eps2rt**2*eps3rt**2
1185 fac=-expon*(e1+evdwij)*rij_shift
1187 fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rij
1189 C Calculate the radial part of the gradient
1193 gg_lipi(3)=eps1*(eps2rt*eps2rt)
1194 & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
1195 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
1196 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
1197 gg_lipj(3)=ssgradlipj*gg_lipi(3)
1198 gg_lipi(3)=gg_lipi(3)*ssgradlipi
1199 c write (iout,*) "gglip",i,j,gg_lipi,gg_lipj
1200 C Calculate angular part of the gradient.
1201 call sc_grad_scale(sss)
1206 c time_evdw_short=time_evdw_short+MPI_Wtime()-time01
1207 c write (iout,*) "Number of loop steps in EGB:",ind
1208 cccc energy_dec=.false.
1211 C-----------------------------------------------------------------------------
1212 subroutine egbv_long(evdw)
1214 C This subroutine calculates the interaction energy of nonbonded side chains
1215 C assuming the Gay-Berne-Vorobjev potential of interaction.
1217 implicit real*8 (a-h,o-z)
1218 include 'DIMENSIONS'
1219 include 'COMMON.GEO'
1220 include 'COMMON.VAR'
1221 include 'COMMON.LOCAL'
1222 include 'COMMON.CHAIN'
1223 include 'COMMON.DERIV'
1224 include 'COMMON.NAMES'
1225 include 'COMMON.INTERACT'
1226 include 'COMMON.IOUNITS'
1227 include 'COMMON.CALC'
1228 include "COMMON.SPLITELE"
1230 common /srutu/ icall
1232 integer itypi,itypj,itypi1,iint,ind,ikont
1233 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1234 & xi,yi,zi,fac_augm,e_augm
1235 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1236 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1237 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1238 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1239 double precision sss1,sssgrad1
1243 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1245 c if (icall.eq.0) lprn=.true.
1247 c do i=iatsc_s,iatsc_e
1248 do ikont=g_listscsc_start_long,g_listscsc_end_long
1249 i=newcontlisti_long(ikont)
1250 j=newcontlistj_long(ikont)
1251 itypi=iabs(itype(i))
1252 if (itypi.eq.ntyp1) cycle
1253 itypi1=iabs(itype(i+1))
1257 call to_box(xi,yi,zi)
1258 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
1259 dxi=dc_norm(1,nres+i)
1260 dyi=dc_norm(2,nres+i)
1261 dzi=dc_norm(3,nres+i)
1262 c dsci_inv=dsc_inv(itypi)
1263 dsci_inv=vbld_inv(i+nres)
1265 C Calculate SC interaction energy.
1267 c do iint=1,nint_gr(i)
1268 c do j=istart(i,iint),iend(i,iint)
1270 itypj=iabs(itype(j))
1271 if (itypj.eq.ntyp1) cycle
1272 c dscj_inv=dsc_inv(itypj)
1273 dscj_inv=vbld_inv(j+nres)
1274 sig0ij=sigma(itypi,itypj)
1275 r0ij=r0(itypi,itypj)
1276 chi1=chi(itypi,itypj)
1277 chi2=chi(itypj,itypi)
1284 alf12=0.5D0*(alf1+alf2)
1288 call to_box(xj,yj,zj)
1289 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
1290 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1291 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1292 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1293 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1294 xj=boxshift(xj-xi,boxxsize)
1295 yj=boxshift(yj-yi,boxysize)
1296 zj=boxshift(zj-zi,boxzsize)
1297 dxj=dc_norm(1,nres+j)
1298 dyj=dc_norm(2,nres+j)
1299 dzj=dc_norm(3,nres+j)
1300 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1303 sss1=sscale(1.0d0/rij,r_cut_int)
1304 if (sss1.eq.0.0d0) cycle
1306 if (sss.lt.1.0d0) then
1308 & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
1309 sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
1311 C Calculate angle-dependent terms of energy and contributions to their
1315 sig=sig0ij*dsqrt(sigsq)
1316 rij_shift=1.0D0/rij-sig+r0ij
1317 C I hate to put IF's in the loops, but here don't have another choice!!!!
1318 if (rij_shift.le.0.0D0) then
1323 c---------------------------------------------------------------
1324 rij_shift=1.0D0/rij_shift
1325 fac=rij_shift**expon
1329 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1330 eps2der=evdwij*eps3rt
1331 eps3der=evdwij*eps2rt
1332 fac_augm=rrij**expon
1333 e_augm=augm(itypi,itypj)*fac_augm
1334 evdwij=evdwij*eps2rt*eps3rt
1335 evdw=evdw+(evdwij+e_augm)*sss1*(1.0d0-sss)
1337 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1339 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1340 & restyp(itypi),i,restyp(itypj),j,
1341 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1342 & chi1,chi2,chip1,chip2,
1343 & eps1,eps2rt**2,eps3rt**2,
1344 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1347 C Calculate gradient components.
1348 e1=e1*eps1*eps2rt**2*eps3rt**2
1349 fac=-expon*(e1+evdwij)*rij_shift
1351 fac=rij*fac-2*expon*rrij*e_augm
1352 fac=fac+(evdwij+e_augm)*
1353 & (-sss1*sssgrad/(1.0d0-sss)/sigmaii(itypi,itypj)
1354 & +(1.0d0-sss)*sssgrad1/sss1)*rij
1355 C Calculate the radial part of the gradient
1359 gg_lipi(3)=eps1*(eps2rt*eps2rt)
1360 & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
1361 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
1362 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
1363 gg_lipj(3)=ssgradlipj*gg_lipi(3)
1364 gg_lipi(3)=gg_lipi(3)*ssgradlipi
1365 C Calculate angular part of the gradient.
1366 call sc_grad_scale((1.0d0-sss)*sss1)
1372 C-----------------------------------------------------------------------------
1373 subroutine egbv_short(evdw)
1375 C This subroutine calculates the interaction energy of nonbonded side chains
1376 C assuming the Gay-Berne-Vorobjev potential of interaction.
1379 include 'DIMENSIONS'
1380 include 'COMMON.GEO'
1381 include 'COMMON.VAR'
1382 include 'COMMON.LOCAL'
1383 include 'COMMON.CHAIN'
1384 include 'COMMON.DERIV'
1385 include 'COMMON.NAMES'
1386 include 'COMMON.INTERACT'
1387 include 'COMMON.IOUNITS'
1388 include 'COMMON.CALC'
1389 include "COMMON.SPLITELE"
1391 common /srutu/ icall
1393 integer itypi,itypj,itypi1,iint,ind,ikont
1394 double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
1395 & xi,yi,zi,fac_augm,e_augm
1396 double precision evdw
1397 double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
1398 & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
1399 & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
1400 double precision dist,sscale,sscagrad,sscagradlip,sscalelip
1401 double precision boxshift
1405 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
1407 c if (icall.eq.0) lprn=.true.
1409 c do i=iatsc_s,iatsc_e
1410 do ikont=g_listscsc_start_short,g_listscsc_end_short
1411 i=newcontlisti_short(ikont)
1412 j=newcontlistj_short(ikont)
1413 itypi=iabs(itype(i))
1414 if (itypi.eq.ntyp1) cycle
1415 itypi1=iabs(itype(i+1))
1419 dxi=dc_norm(1,nres+i)
1420 dyi=dc_norm(2,nres+i)
1421 dzi=dc_norm(3,nres+i)
1422 call to_box(xi,yi,zi)
1423 call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
1424 c dsci_inv=dsc_inv(itypi)
1425 dsci_inv=vbld_inv(i+nres)
1427 C Calculate SC interaction energy.
1429 c do iint=1,nint_gr(i)
1430 c do j=istart(i,iint),iend(i,iint)
1432 itypj=iabs(itype(j))
1433 if (itypj.eq.ntyp1) cycle
1434 c dscj_inv=dsc_inv(itypj)
1435 dscj_inv=vbld_inv(j+nres)
1436 sig0ij=sigma(itypi,itypj)
1437 r0ij=r0(itypi,itypj)
1438 chi1=chi(itypi,itypj)
1439 chi2=chi(itypj,itypi)
1446 alf12=0.5D0*(alf1+alf2)
1450 call to_box(xj,yj,zj)
1451 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
1452 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1453 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1454 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
1455 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
1456 xj=boxshift(xj-xi,boxxsize)
1457 yj=boxshift(yj-yi,boxysize)
1458 zj=boxshift(zj-zi,boxzsize)
1459 dxj=dc_norm(1,nres+j)
1460 dyj=dc_norm(2,nres+j)
1461 dzj=dc_norm(3,nres+j)
1462 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1465 sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
1467 if (sss.gt.0.0d0) then
1469 C Calculate angle-dependent terms of energy and contributions to their
1473 sig=sig0ij*dsqrt(sigsq)
1474 rij_shift=1.0D0/rij-sig+r0ij
1475 C I hate to put IF's in the loops, but here don't have another choice!!!!
1476 if (rij_shift.le.0.0D0) then
1481 c---------------------------------------------------------------
1482 rij_shift=1.0D0/rij_shift
1483 fac=rij_shift**expon
1487 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
1488 eps2der=evdwij*eps3rt
1489 eps3der=evdwij*eps2rt
1490 fac_augm=rrij**expon
1491 e_augm=augm(itypi,itypj)*fac_augm
1492 evdwij=evdwij*eps2rt*eps3rt
1493 evdw=evdw+(evdwij+e_augm)*sss
1495 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
1497 write (iout,'(2(a3,i3,2x),17(0pf7.3))')
1498 & restyp(itypi),i,restyp(itypj),j,
1499 & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
1500 & chi1,chi2,chip1,chip2,
1501 & eps1,eps2rt**2,eps3rt**2,
1502 & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
1505 C Calculate gradient components.
1506 e1=e1*eps1*eps2rt**2*eps3rt**2
1507 fac=-expon*(e1+evdwij)*rij_shift
1509 fac=rij*fac-2*expon*rrij*e_augm+
1510 & (evdwij+e_augm)*sssgrad/sigmaii(itypi,itypj)/sss*rij
1511 C Calculate the radial part of the gradient
1515 gg_lipi(3)=eps1*(eps2rt*eps2rt)
1516 & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
1517 & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
1518 & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
1519 gg_lipj(3)=ssgradlipj*gg_lipi(3)
1520 gg_lipi(3)=gg_lipi(3)*ssgradlipi
1521 C Calculate angular part of the gradient.
1522 call sc_grad_scale(sss)
1528 C----------------------------------------------------------------------------
1529 subroutine sc_grad_scale(scalfac)
1530 implicit real*8 (a-h,o-z)
1531 include 'DIMENSIONS'
1532 include 'COMMON.CHAIN'
1533 include 'COMMON.DERIV'
1534 include 'COMMON.CALC'
1535 include 'COMMON.IOUNITS'
1536 include "COMMON.SPLITELE"
1537 double precision dcosom1(3),dcosom2(3)
1538 double precision scalfac
1539 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1540 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1541 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1542 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1546 c eom12=evdwij*eps1_om12
1548 c write (iout,*) "eps2der",eps2der," eps3der",eps3der,
1549 c & " sigder",sigder
1550 c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
1551 c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
1553 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1554 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1557 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
1559 c write (iout,*) "gg",(gg(k),k=1,3)
1561 gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
1562 & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1563 & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
1564 gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
1565 & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1566 & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
1567 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
1568 c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
1569 c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
1570 c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
1573 C Calculate the components of the gradient in DC and X
1575 c write (iout,*) "scgrad gglip",i,j,gg_lipi,gg_lipj
1577 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
1578 gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
1582 C--------------------------------------------------------------------------
1583 subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1585 C This subroutine calculates the average interaction energy and its gradient
1586 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1587 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1588 C The potential depends both on the distance of peptide-group centers and on
1589 C the orientation of the CA-CA virtual bonds.
1591 implicit real*8 (a-h,o-z)
1595 include 'DIMENSIONS'
1596 include 'COMMON.CONTROL'
1597 include 'COMMON.SETUP'
1598 include 'COMMON.IOUNITS'
1599 include 'COMMON.GEO'
1600 include 'COMMON.VAR'
1601 include 'COMMON.LOCAL'
1602 include 'COMMON.CHAIN'
1603 include 'COMMON.DERIV'
1604 include 'COMMON.INTERACT'
1606 include 'COMMON.CONTACTS'
1607 include 'COMMON.CONTMAT'
1609 include 'COMMON.CORRMAT'
1610 include 'COMMON.TORSION'
1611 include 'COMMON.VECTORS'
1612 include 'COMMON.FFIELD'
1613 include 'COMMON.TIME1'
1614 include 'COMMON.SHIELD'
1615 include "COMMON.SPLITELE"
1617 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1618 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1619 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1620 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1621 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1622 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1624 double precision sslipi,sslipj,ssgradlipi,ssgradlipj
1625 common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj
1626 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1628 double precision scal_el /1.0d0/
1630 double precision scal_el /0.5d0/
1633 C 13-go grudnia roku pamietnego...
1634 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1635 & 0.0d0,1.0d0,0.0d0,
1636 & 0.0d0,0.0d0,1.0d0/
1637 cd write(iout,*) 'In EELEC'
1639 cd write(iout,*) 'Type',i
1640 cd write(iout,*) 'B1',B1(:,i)
1641 cd write(iout,*) 'B2',B2(:,i)
1642 cd write(iout,*) 'CC',CC(:,:,i)
1643 cd write(iout,*) 'DD',DD(:,:,i)
1644 cd write(iout,*) 'EE',EE(:,:,i)
1646 cd call check_vecgrad
1648 C print *,"WCHODZE3"
1649 if (icheckgrad.eq.1) then
1651 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1653 dc_norm(k,i)=dc(k,i)*fac
1655 c write (iout,*) 'i',i,' fac',fac
1658 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1659 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1660 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1661 c call vec_and_deriv
1667 time_mat=time_mat+MPI_Wtime()-time01
1671 cd write (iout,*) 'i=',i
1673 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1676 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1677 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1692 cd print '(a)','Enter EELEC'
1693 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1695 gel_loc_loc(i)=0.0d0
1700 c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
1702 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
1704 do i=iturn3_start,iturn3_end
1705 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
1706 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
1707 C & .or. itype(i-1).eq.ntyp1
1708 C & .or. itype(i+4).eq.ntyp1
1713 dx_normi=dc_norm(1,i)
1714 dy_normi=dc_norm(2,i)
1715 dz_normi=dc_norm(3,i)
1716 xmedi=c(1,i)+0.5d0*dxi
1717 ymedi=c(2,i)+0.5d0*dyi
1718 zmedi=c(3,i)+0.5d0*dzi
1719 call to_box(xmedi,ymedi,zmedi)
1720 call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
1722 call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
1723 if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
1725 num_cont_hb(i)=num_conti
1728 do i=iturn4_start,iturn4_end
1729 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1730 & .or. itype(i+3).eq.ntyp1
1731 & .or. itype(i+4).eq.ntyp1
1732 C & .or. itype(i+5).eq.ntyp1
1733 C & .or. itype(i-1).eq.ntyp1
1738 dx_normi=dc_norm(1,i)
1739 dy_normi=dc_norm(2,i)
1740 dz_normi=dc_norm(3,i)
1741 xmedi=c(1,i)+0.5d0*dxi
1742 ymedi=c(2,i)+0.5d0*dyi
1743 zmedi=c(3,i)+0.5d0*dzi
1744 call to_box(xmedi,ymedi,zmedi)
1745 call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
1747 num_conti=num_cont_hb(i)
1749 call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
1750 if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
1751 & call eturn4(i,eello_turn4)
1753 num_cont_hb(i)=num_conti
1757 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
1759 c do i=iatel_s,iatel_e
1761 & write(iout,*) "g_listpp_start,g_listpp_end",
1762 & g_listpp_start,g_listpp_end
1763 do ikont=g_listpp_start,g_listpp_end
1764 i=newcontlistppi(ikont)
1765 j=newcontlistppj(ikont)
1766 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
1767 C & .or. itype(i+2).eq.ntyp1
1768 C & .or. itype(i-1).eq.ntyp1
1773 dx_normi=dc_norm(1,i)
1774 dy_normi=dc_norm(2,i)
1775 dz_normi=dc_norm(3,i)
1776 xmedi=c(1,i)+0.5d0*dxi
1777 ymedi=c(2,i)+0.5d0*dyi
1778 zmedi=c(3,i)+0.5d0*dzi
1779 call to_box(xmedi,ymedi,zmedi)
1780 call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
1781 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
1783 num_conti=num_cont_hb(i)
1785 c do j=ielstart(i),ielend(i)
1786 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
1787 C & .or.itype(j+2).eq.ntyp1
1788 C & .or.itype(j-1).eq.ntyp1
1790 call eelecij_scale(i,j,ees,evdw1,eel_loc)
1793 num_cont_hb(i)=num_conti
1796 c write (iout,*) "Number of loop steps in EELEC:",ind
1798 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
1799 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
1801 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
1802 ccc eel_loc=eel_loc+eello_turn3
1803 cd print *,"Processor",fg_rank," t_eelecij",t_eelecij
1806 C-------------------------------------------------------------------------------
1807 subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
1809 include 'DIMENSIONS'
1813 include 'COMMON.CONTROL'
1814 include 'COMMON.IOUNITS'
1815 include 'COMMON.GEO'
1816 include 'COMMON.VAR'
1817 include 'COMMON.LOCAL'
1818 include 'COMMON.CHAIN'
1819 include 'COMMON.DERIV'
1820 include 'COMMON.INTERACT'
1822 include 'COMMON.CONTACTS'
1823 include 'COMMON.CONTMAT'
1825 include 'COMMON.CORRMAT'
1826 include 'COMMON.TORSION'
1827 include 'COMMON.VECTORS'
1828 include 'COMMON.FFIELD'
1829 include 'COMMON.TIME1'
1830 include 'COMMON.SHIELD'
1831 include "COMMON.SPLITELE"
1832 integer xshift,yshift,zshift
1833 double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1834 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1835 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1836 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
1837 & gmuij2(4),gmuji2(4)
1838 integer j1,j2,num_conti
1839 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
1840 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
1842 integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ind,itypi,itypj
1843 integer ilist,iresshield
1844 double precision rlocshield
1845 double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
1846 double precision ees,evdw1,eel_loc,aaa,bbb,ael3i
1847 double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,
1848 & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,
1849 & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa,
1850 & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,
1851 & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,
1852 & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,
1853 & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
1854 & ecosgp,ecosam,ecosbm,ecosgm,ghalf,geel_loc_ij,geel_loc_ji,
1855 & dxi,dyi,dzi,a22,a23,a32,a33
1856 double precision dist_init,xmedi,ymedi,zmedi,xj_safe,yj_safe,
1857 & zj_safe,xj_temp,yj_temp,zj_temp,dist_temp,dx_normi,dy_normi,
1859 double precision sss1,sssgrad1
1860 double precision sscale,sscagrad
1861 double precision scalar
1862 double precision boxshift
1863 double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij,
1865 common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij
1866 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1868 double precision scal_el /1.0d0/
1870 double precision scal_el /0.5d0/
1873 C 13-go grudnia roku pamietnego...
1874 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1875 & 0.0d0,1.0d0,0.0d0,
1876 & 0.0d0,0.0d0,1.0d0/
1877 c time00=MPI_Wtime()
1878 cd write (iout,*) "eelecij",i,j
1879 C print *,"WCHODZE2"
1883 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1884 aaa=app(iteli,itelj)
1885 bbb=bpp(iteli,itelj)
1886 ael6i=ael6(iteli,itelj)
1887 ael3i=ael3(iteli,itelj)
1891 dx_normj=dc_norm(1,j)
1892 dy_normj=dc_norm(2,j)
1893 dz_normj=dc_norm(3,j)
1897 call to_box(xj,yj,zj)
1898 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
1899 faclipij=(sslipi+sslipj)/2.0d0*lipscale+1.0d0
1900 faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
1901 xj=boxshift(xj-xmedi,boxxsize)
1902 yj=boxshift(yj-ymedi,boxysize)
1903 zj=boxshift(zj-zmedi,boxzsize)
1904 rij=xj*xj+yj*yj+zj*zj
1908 c For extracting the short-range part of Evdwpp
1909 sss1=sscale(rij,r_cut_int)
1910 if (sss1.eq.0.0d0) return
1911 sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
1912 sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
1913 sssgrad1=sscagrad(rij,r_cut_int)
1916 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1917 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1918 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1919 fac=cosa-3.0D0*cosb*cosg
1921 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1922 if (j.eq.i+2) ev1=scal_el*ev1
1927 if (shield_mode.eq.0) then
1931 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1933 el1=el1*fac_shield(i)**2*fac_shield(j)**2
1934 el2=el2*fac_shield(i)**2*fac_shield(j)**2
1936 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1937 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1938 ees=ees+eesij*sss1*faclipij2
1939 evdw1=evdw1+evdwij*(1.0d0-sss)*sss1*faclipij2
1940 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1941 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1942 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1943 cd & xmedi,ymedi,zmedi,xj,yj,zj
1945 if (energy_dec) then
1946 write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,5f10.5)')
1947 & 'evdw1',i,j,evdwij,iteli,itelj,aaa,sss,sss1,sssgrad,sssgrad1,rij
1948 write (iout,'(a6,2i5,0pf7.3,6f8.5)') 'ees',i,j,eesij,
1949 & fac_shield(i),fac_shield(j),sslipi,sslipj,faclipij,faclipij2
1953 C Calculate contributions to the Cartesian gradient.
1956 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
1957 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1958 facel=-3*rrmij*(el1+eesij)*sss1
1959 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1965 * Radial derivatives. First process both termini of the fragment (i,j)
1967 c old aux=(facel+sssgrad1*(1.0d0-sss)*eesij*rmij)*faclipij2
1968 aux=(facel+sssgrad1*eesij*rmij)*faclipij2
1969 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
1976 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1977 & (shield_mode.gt.0)) then
1979 do ilist=1,ishield_list(i)
1980 iresshield=shield_list(ilist,i)
1982 rlocshield=grad_shield_side(k,ilist,i)*eesij*sss1
1983 & /fac_shield(i)*2.0*sss1
1984 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1986 & +grad_shield_loc(k,ilist,i)*eesij*sss1/fac_shield(i)*2.0
1988 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1989 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1990 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1991 C if (iresshield.gt.i) then
1992 C do ishi=i+1,iresshield-1
1993 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1994 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1998 C do ishi=iresshield,i
1999 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
2000 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
2006 do ilist=1,ishield_list(j)
2007 iresshield=shield_list(ilist,j)
2009 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
2011 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
2013 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss1
2014 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
2019 gshieldc(k,i)=gshieldc(k,i)+
2020 & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
2021 gshieldc(k,j)=gshieldc(k,j)+
2022 & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
2023 gshieldc(k,i-1)=gshieldc(k,i-1)+
2024 & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
2025 gshieldc(k,j-1)=gshieldc(k,j-1)+
2026 & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
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_long(k,j)+ggg(k)
2039 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
2041 gelc_long(3,j)=gelc_long(3,j)+
2042 & ssgradlipj*eesij/2.0d0*lipscale**2*sss1
2043 gelc_long(3,i)=gelc_long(3,i)+
2044 & ssgradlipi*eesij/2.0d0*lipscale**2*sss1
2045 c gelc_long(3,i)=gelc_long(3,i)+
2046 c ssgradlipi*eesij/2.0d0*lipscale**2*sss1
2048 * Loop over residues i+1 thru j-1.
2052 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
2056 &(-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij)
2062 c ghalf=0.5D0*ggg(k)
2063 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
2064 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
2066 c 9/28/08 AL Gradient compotents will be summed only at the end
2068 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2069 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2071 !C Lipidic part for scaling weight
2072 gvdwpp(3,j)=gvdwpp(3,j)+
2073 & sss1*(1.0d0-sss)*ssgradlipj*evdwij/2.0d0*lipscale**2
2074 gvdwpp(3,i)=gvdwpp(3,i)+
2075 & sss1*(1.0d0-sss)*ssgradlipi*evdwij/2.0d0*lipscale**2
2077 * Loop over residues i+1 thru j-1.
2081 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
2085 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
2086 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2087 facel=-3*rrmij*(el1+eesij)*sss1
2088 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2090 c facvdw=ev1+evdwij*(1.0d0-sss)*sss1
2093 fac=-3*rrmij*(facvdw+facvdw+facel)
2098 * Radial derivatives. First process both termini of the fragment (i,j)
2100 aux=(fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj))
2101 & *eesij*rmij)*faclipij2
2102 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2110 c ghalf=0.5D0*ggg(k)
2111 c gelc(k,i)=gelc(k,i)+ghalf
2112 c gelc(k,j)=gelc(k,j)+ghalf
2114 c 9/28/08 AL Gradient compotents will be summed only at the end
2116 gelc_long(k,j)=gelc(k,j)+ggg(k)
2117 gelc_long(k,i)=gelc(k,i)-ggg(k)
2120 * Loop over residues i+1 thru j-1.
2124 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
2127 c 9/28/08 AL Gradient compotents will be summed only at the end
2132 & (-sssgrad*sss1/rpp(iteli,itelj)+sssgrad1*(1.0d0-sss))*rmij*evdwij
2137 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2138 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2140 gvdwpp(3,j)=gvdwpp(3,j)+
2141 & sss1*ssgradlipj*evdwij/2.0d0*lipscale**2
2142 gvdwpp(3,i)=gvdwpp(3,i)+
2143 & sss1*ssgradlipi*evdwij/2.0d0*lipscale**2
2148 ecosa=2.0D0*fac3*fac1+fac4
2151 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
2152 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
2154 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2155 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2157 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
2158 cd & (dcosg(k),k=1,3)
2160 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1
2161 & *fac_shield(i)**2*fac_shield(j)**2*faclipij2
2162 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2166 c ghalf=0.5D0*ggg(k)
2167 c gelc(k,i)=gelc(k,i)+ghalf
2168 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
2169 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2170 c gelc(k,j)=gelc(k,j)+ghalf
2171 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2172 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2176 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
2181 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
2182 & +ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss1
2183 & *fac_shield(i)**2*fac_shield(j)**2*faclipij2
2184 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2187 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
2188 & +ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss1
2189 & *fac_shield(i)**2*fac_shield(j)**2*faclipij2
2190 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
2191 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
2192 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
2194 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
2195 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
2196 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2198 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
2199 C energy of a peptide unit is assumed in the form of a second-order
2200 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
2201 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
2202 C are computed for EVERY pair of non-contiguous peptide groups.
2204 if (j.lt.nres-1) then
2215 muij(kkk)=mu(k,i)*mu(l,j)
2217 gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
2218 c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
2219 gmuij2(kkk)=gUb2(k,i)*mu(l,j)
2220 gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
2221 c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
2222 gmuji2(kkk)=mu(k,i)*gUb2(l,j)
2226 cd write (iout,*) 'EELEC: i',i,' j',j
2227 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
2228 cd write(iout,*) 'muij',muij
2229 ury=scalar(uy(1,i),erij)
2230 urz=scalar(uz(1,i),erij)
2231 vry=scalar(uy(1,j),erij)
2232 vrz=scalar(uz(1,j),erij)
2233 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
2234 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
2235 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
2236 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
2237 fac=dsqrt(-ael6i)*r3ij
2242 cd write (iout,'(4i5,4f10.5)')
2243 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
2244 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
2245 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
2246 cd & uy(:,j),uz(:,j)
2247 cd write (iout,'(4f10.5)')
2248 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
2249 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
2250 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
2251 cd write (iout,'(9f10.5/)')
2252 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
2253 C Derivatives of the elements of A in virtual-bond vectors
2254 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
2256 uryg(k,1)=scalar(erder(1,k),uy(1,i))
2257 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
2258 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
2259 urzg(k,1)=scalar(erder(1,k),uz(1,i))
2260 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
2261 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
2262 vryg(k,1)=scalar(erder(1,k),uy(1,j))
2263 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
2264 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
2265 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
2266 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
2267 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
2269 C Compute radial contributions to the gradient
2287 C Add the contributions coming from er
2290 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
2291 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
2292 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
2293 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
2296 C Derivatives in DC(i)
2297 cgrad ghalf1=0.5d0*agg(k,1)
2298 cgrad ghalf2=0.5d0*agg(k,2)
2299 cgrad ghalf3=0.5d0*agg(k,3)
2300 cgrad ghalf4=0.5d0*agg(k,4)
2301 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2302 & -3.0d0*uryg(k,2)*vry)!+ghalf1
2303 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2304 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
2305 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2306 & -3.0d0*urzg(k,2)*vry)!+ghalf3
2307 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2308 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
2309 C Derivatives in DC(i+1)
2310 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2311 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
2312 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2313 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
2314 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2315 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
2316 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2317 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
2318 C Derivatives in DC(j)
2319 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2320 & -3.0d0*vryg(k,2)*ury)!+ghalf1
2321 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2322 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
2323 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2324 & -3.0d0*vryg(k,2)*urz)!+ghalf3
2325 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2326 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
2327 C Derivatives in DC(j+1) or DC(nres-1)
2328 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2329 & -3.0d0*vryg(k,3)*ury)
2330 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2331 & -3.0d0*vrzg(k,3)*ury)
2332 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2333 & -3.0d0*vryg(k,3)*urz)
2334 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2335 & -3.0d0*vrzg(k,3)*urz)
2336 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
2338 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
2351 aggi(k,l)=-aggi(k,l)
2352 aggi1(k,l)=-aggi1(k,l)
2353 aggj(k,l)=-aggj(k,l)
2354 aggj1(k,l)=-aggj1(k,l)
2357 if (j.lt.nres-1) then
2363 aggi(k,l)=-aggi(k,l)
2364 aggi1(k,l)=-aggi1(k,l)
2365 aggj(k,l)=-aggj(k,l)
2366 aggj1(k,l)=-aggj1(k,l)
2377 aggi(k,l)=-aggi(k,l)
2378 aggi1(k,l)=-aggi1(k,l)
2379 aggj(k,l)=-aggj(k,l)
2380 aggj1(k,l)=-aggj1(k,l)
2385 IF (wel_loc.gt.0.0d0) THEN
2386 C Contribution to the local-electrostatic energy coming from the i-j pair
2387 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2389 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2391 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
2392 & 'eelloc',i,j,eel_loc_ij
2395 if (shield_mode.eq.0) then
2402 eel_loc_ij=eel_loc_ij
2403 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2404 eel_loc=eel_loc+eel_loc_ij
2406 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2407 & (shield_mode.gt.0)) then
2410 do ilist=1,ishield_list(i)
2411 iresshield=shield_list(ilist,i)
2413 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
2416 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2418 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
2419 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2423 do ilist=1,ishield_list(j)
2424 iresshield=shield_list(ilist,j)
2426 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
2429 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
2431 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
2432 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
2439 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
2440 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2441 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
2442 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2443 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
2444 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
2445 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
2446 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
2451 geel_loc_ij=(a22*gmuij1(1)
2455 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2456 c write(iout,*) "derivative over thatai"
2457 c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
2459 gloc(nphi+i,icg)=gloc(nphi+i,icg)+
2460 & geel_loc_ij*wel_loc
2461 c write(iout,*) "derivative over thatai-1"
2462 c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
2469 gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
2470 & geel_loc_ij*wel_loc
2471 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2473 c Derivative over j residue
2474 geel_loc_ji=a22*gmuji1(1)
2478 c write(iout,*) "derivative over thataj"
2479 c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
2482 gloc(nphi+j,icg)=gloc(nphi+j,icg)+
2483 & geel_loc_ji*wel_loc
2484 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2491 c write(iout,*) "derivative over thataj-1"
2492 c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
2494 gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
2495 & geel_loc_ji*wel_loc
2496 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2498 cC Paral derivatives in virtual-bond dihedral angles gamma
2500 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2501 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2502 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
2503 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2504 c & *fac_shield(i)*fac_shield(j)
2505 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2508 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2509 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2510 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
2511 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2512 c & *fac_shield(i)*fac_shield(j)
2513 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2515 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2516 aux=eel_loc_ij/sss1*sssgrad1*rmij
2521 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+
2522 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
2523 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2524 c & *fac_shield(i)*fac_shield(j)
2525 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2527 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
2528 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
2529 cgrad ghalf=0.5d0*ggg(l)
2530 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
2531 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
2533 gel_loc_long(3,j)=gel_loc_long(3,j)+
2534 & ssgradlipj*eel_loc_ij/2.0d0*lipscale/faclipij
2535 gel_loc_long(3,i)=gel_loc_long(3,i)+
2536 & ssgradlipi*eel_loc_ij/2.0d0*lipscale/faclipij
2539 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2542 C Remaining derivatives of eello
2543 c gel_loc_long(3,j)=gel_loc_long(3,j)+ &
2544 c ssgradlipj*eel_loc_ij/2.0d0*lipscale/ &
2545 c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2547 c gel_loc_long(3,i)=gel_loc_long(3,i)+ &
2548 c ssgradlipi*eel_loc_ij/2.0d0*lipscale/ &
2549 c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
2552 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
2553 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
2554 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2555 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2557 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
2558 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
2559 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2560 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2562 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
2563 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
2564 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2565 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2567 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
2568 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
2569 & *fac_shield(i)*fac_shield(j)*sss1*faclipij
2570 c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
2575 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2576 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
2577 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
2578 & .and. num_conti.le.maxconts) then
2579 c write (iout,*) i,j," entered corr"
2581 C Calculate the contact function. The ith column of the array JCONT will
2582 C contain the numbers of atoms that make contacts with the atom I (of numbers
2583 C greater than I). The arrays FACONT and GACONT will contain the values of
2584 C the contact function and its derivative.
2585 c r0ij=1.02D0*rpp(iteli,itelj)
2586 c r0ij=1.11D0*rpp(iteli,itelj)
2587 r0ij=2.20D0*rpp(iteli,itelj)
2588 c r0ij=1.55D0*rpp(iteli,itelj)
2589 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2590 if (fcont.gt.0.0D0) then
2591 num_conti=num_conti+1
2592 if (num_conti.gt.maxconts) then
2593 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2594 & ' will skip next contacts for this conf.'
2596 jcont_hb(num_conti,i)=j
2597 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
2598 cd " jcont_hb",jcont_hb(num_conti,i)
2599 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2600 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2601 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2603 d_cont(num_conti,i)=rij
2604 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2605 C --- Electrostatic-interaction matrix ---
2606 a_chuj(1,1,num_conti,i)=a22
2607 a_chuj(1,2,num_conti,i)=a23
2608 a_chuj(2,1,num_conti,i)=a32
2609 a_chuj(2,2,num_conti,i)=a33
2610 C --- Gradient of rij
2612 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2619 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2620 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2621 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2622 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2623 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2628 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2629 C Calculate contact energies
2631 wij=cosa-3.0D0*cosb*cosg
2634 c fac3=dsqrt(-ael6i)/r0ij**3
2635 fac3=dsqrt(-ael6i)*r3ij
2636 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2637 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
2638 if (ees0tmp.gt.0) then
2639 ees0pij=dsqrt(ees0tmp)
2643 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2644 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
2645 if (ees0tmp.gt.0) then
2646 ees0mij=dsqrt(ees0tmp)
2651 if (shield_mode.eq.0) then
2655 ees0plist(num_conti,i)=j
2656 C fac_shield(i)=0.4d0
2657 C fac_shield(j)=0.6d0
2659 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2660 & *fac_shield(i)*fac_shield(j)*sss1
2661 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2662 & *fac_shield(i)*fac_shield(j)*sss1
2664 C Diagnostics. Comment out or remove after debugging!
2665 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2666 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2667 c ees0m(num_conti,i)=0.0D0
2669 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2670 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2671 C Angular derivatives of the contact function
2672 ees0pij1=fac3/ees0pij
2673 ees0mij1=fac3/ees0mij
2674 fac3p=-3.0D0*fac3*rrmij
2675 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2676 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2678 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2679 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2680 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2681 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2682 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2683 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2684 ecosap=ecosa1+ecosa2
2685 ecosbp=ecosb1+ecosb2
2686 ecosgp=ecosg1+ecosg2
2687 ecosam=ecosa1-ecosa2
2688 ecosbm=ecosb1-ecosb2
2689 ecosgm=ecosg1-ecosg2
2698 facont_hb(num_conti,i)=fcont
2699 fprimcont=fprimcont/rij
2700 cd facont_hb(num_conti,i)=1.0D0
2701 C Following line is for diagnostics.
2704 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2705 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2708 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2709 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2711 gggp(1)=gggp(1)+ees0pijp*xj
2712 & +ees0p(num_conti,i)/sss1*rmij*xj*sssgrad1
2713 gggp(2)=gggp(2)+ees0pijp*yj
2714 & +ees0p(num_conti,i)/sss1*rmij*yj*sssgrad1
2715 gggp(3)=gggp(3)+ees0pijp*zj
2716 & +ees0p(num_conti,i)/sss1*rmij*zj*sssgrad1
2717 gggm(1)=gggm(1)+ees0mijp*xj
2718 & +ees0m(num_conti,i)/sss1*rmij*xj*sssgrad1
2719 gggm(2)=gggm(2)+ees0mijp*yj
2720 & +ees0m(num_conti,i)/sss1*rmij*yj*sssgrad1
2721 gggm(3)=gggm(3)+ees0mijp*zj
2722 & +ees0m(num_conti,i)/sss1*rmij*zj*sssgrad1
2723 C Derivatives due to the contact function
2724 gacont_hbr(1,num_conti,i)=fprimcont*xj
2725 gacont_hbr(2,num_conti,i)=fprimcont*yj
2726 gacont_hbr(3,num_conti,i)=fprimcont*zj
2729 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
2730 c following the change of gradient-summation algorithm.
2732 cgrad ghalfp=0.5D0*gggp(k)
2733 cgrad ghalfm=0.5D0*gggm(k)
2734 gacontp_hb1(k,num_conti,i)=!ghalfp
2735 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
2736 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2737 & *sss1*fac_shield(i)*fac_shield(j)
2739 gacontp_hb2(k,num_conti,i)=!ghalfp
2740 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
2741 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2742 & *sss1*fac_shield(i)*fac_shield(j)
2744 gacontp_hb3(k,num_conti,i)=gggp(k)
2745 & *sss1*fac_shield(i)*fac_shield(j)
2747 gacontm_hb1(k,num_conti,i)=!ghalfm
2748 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
2749 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
2750 & *sss1*fac_shield(i)*fac_shield(j)
2752 gacontm_hb2(k,num_conti,i)=!ghalfm
2753 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
2754 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
2755 & *sss1*fac_shield(i)*fac_shield(j)
2757 gacontm_hb3(k,num_conti,i)=gggm(k)
2758 & *sss1*fac_shield(i)*fac_shield(j)
2762 endif ! num_conti.le.maxconts
2766 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2769 ghalf=0.5d0*agg(l,k)
2770 aggi(l,k)=aggi(l,k)+ghalf
2771 aggi1(l,k)=aggi1(l,k)+agg(l,k)
2772 aggj(l,k)=aggj(l,k)+ghalf
2775 if (j.eq.nres-1 .and. i.lt.j-2) then
2778 aggj1(l,k)=aggj1(l,k)+agg(l,k)
2783 c t_eelecij=t_eelecij+MPI_Wtime()-time00
2786 C-----------------------------------------------------------------------
2787 subroutine evdwpp_short(evdw1)
2792 include 'DIMENSIONS'
2793 include 'COMMON.CONTROL'
2794 include 'COMMON.IOUNITS'
2795 include 'COMMON.GEO'
2796 include 'COMMON.VAR'
2797 include 'COMMON.LOCAL'
2798 include 'COMMON.CHAIN'
2799 include 'COMMON.DERIV'
2800 include 'COMMON.INTERACT'
2801 c include 'COMMON.CONTACTS'
2802 include 'COMMON.TORSION'
2803 include 'COMMON.VECTORS'
2804 include 'COMMON.FFIELD'
2805 include "COMMON.SPLITELE"
2806 double precision ggg(3)
2807 integer xshift,yshift,zshift
2808 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
2810 double precision scal_el /1.0d0/
2812 double precision scal_el /0.5d0/
2814 c write (iout,*) "evdwpp_short"
2815 integer i,j,k,iteli,itelj,num_conti,ind,isubchap
2816 double precision dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb
2817 double precision xj,yj,zj,rij,rrmij,r3ij,r6ij,evdw1,
2818 & dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
2819 & dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw
2820 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
2821 & dist_temp, dist_init,sss_grad
2822 double precision sscale,sscagrad
2823 double precision sslipi,ssgradlipi,sslipj,ssgradlipj
2824 double precision boxshift
2826 double precision faclipij2
2829 c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
2830 c & " iatel_e_vdw",iatel_e_vdw
2832 c do i=iatel_s_vdw,iatel_e_vdw
2834 & write(iout,*) "g_listpp_vdw_start_short,g_listpp_vdw_end_short",
2835 & g_listpp_vdw_start_short,g_listpp_vdw_end_short
2836 do ikont=g_listpp_vdw_start_short,g_listpp_vdw_end_short
2837 i=newcontlistpp_vdwi_short(ikont)
2838 j=newcontlistpp_vdwj_short(ikont)
2839 if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
2843 dx_normi=dc_norm(1,i)
2844 dy_normi=dc_norm(2,i)
2845 dz_normi=dc_norm(3,i)
2846 xmedi=c(1,i)+0.5d0*dxi
2847 ymedi=c(2,i)+0.5d0*dyi
2848 zmedi=c(3,i)+0.5d0*dzi
2849 call to_box(xmedi,ymedi,zmedi)
2850 call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
2852 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
2853 c & ' ielend',ielend_vdw(i)
2855 c do j=ielstart_vdw(i),ielend_vdw(i)
2856 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
2860 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
2861 aaa=app(iteli,itelj)
2862 bbb=bpp(iteli,itelj)
2866 dx_normj=dc_norm(1,j)
2867 dy_normj=dc_norm(2,j)
2868 dz_normj=dc_norm(3,j)
2872 call to_box(xj,yj,zj)
2873 call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
2874 faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
2875 xj=boxshift(xj-xmedi,boxxsize)
2876 yj=boxshift(yj-ymedi,boxysize)
2877 zj=boxshift(zj-zmedi,boxzsize)
2878 rij=xj*xj+yj*yj+zj*zj
2881 c sss=sscale(rij/rpp(iteli,itelj))
2882 c sssgrad=sscagrad(rij/rpp(iteli,itelj))
2883 sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
2884 sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
2885 if (sss.gt.0.0d0) then
2890 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
2891 if (j.eq.i+2) ev1=scal_el*ev1
2894 if (energy_dec) then
2895 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
2897 evdw1=evdw1+evdwij*sss*faclipij2
2898 if (energy_dec) write (iout,'(a10,2i5,0pf7.3)')
2899 & 'evdw1_sum',i,j,evdw1
2901 C Calculate contributions to the Cartesian gradient.
2903 facvdw=(-6*rrmij*(ev1+evdwij)*sss+sssgrad*rmij*evdwij/
2904 & rpp(iteli,itelj))*faclipij2
2912 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
2913 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
2915 !C Lipidic part for scaling weight
2916 gvdwpp(3,j)=gvdwpp(3,j)+
2917 & sss*ssgradlipj*evdwij/2.0d0*lipscale**2
2918 gvdwpp(3,i)=gvdwpp(3,i)+
2919 & sss*ssgradlipi*evdwij/2.0d0*lipscale**2
2925 C-----------------------------------------------------------------------------
2926 subroutine escp_long(evdw2,evdw2_14)
2928 C This subroutine calculates the excluded-volume interaction energy between
2929 C peptide-group centers and side chains and its gradient in virtual-bond and
2930 C side-chain vectors.
2932 implicit real*8 (a-h,o-z)
2933 include 'DIMENSIONS'
2934 include 'COMMON.GEO'
2935 include 'COMMON.VAR'
2936 include 'COMMON.LOCAL'
2937 include 'COMMON.CHAIN'
2938 include 'COMMON.DERIV'
2939 include 'COMMON.INTERACT'
2940 include 'COMMON.FFIELD'
2941 include 'COMMON.IOUNITS'
2942 include 'COMMON.CONTROL'
2943 include "COMMON.SPLITELE"
2944 logical lprint_short
2945 common /shortcheck/ lprint_short
2946 double precision ggg(3)
2947 integer i,iint,j,k,iteli,itypj,subchap
2948 double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
2950 double precision evdw2,evdw2_14,evdwij
2951 double precision sscale,sscagrad
2952 double precision boxshift
2954 if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
2957 CD print '(a)','Enter ESCP KURWA'
2958 cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
2960 c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
2961 c & ' iatscp_e=',iatscp_e
2962 c do i=iatscp_s,iatscp_e
2964 & write(iout,*)"g_listscp_start_long,g_listscp_end_long",
2965 & g_listscp_start_long,g_listscp_end_long
2966 do ikont=g_listscp_start_long,g_listscp_end_long
2967 i=newcontlistscpi_long(ikont)
2968 j=newcontlistscpj_long(ikont)
2969 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2971 xi=0.5D0*(c(1,i)+c(1,i+1))
2972 yi=0.5D0*(c(2,i)+c(2,i+1))
2973 zi=0.5D0*(c(3,i)+c(3,i+1))
2974 call to_box(xi,yi,zi)
2976 c do iint=1,nscp_gr(i)
2978 c do j=iscpstart(i,iint),iscpend(i,iint)
2979 itypj=iabs(itype(j))
2980 if (itypj.eq.ntyp1) cycle
2981 C Uncomment following three lines for SC-p interactions
2985 C Uncomment following three lines for Ca-p interactions
2990 call to_box(xj,yj,zj)
2991 xj=boxshift(xj-xi,boxxsize)
2992 yj=boxshift(yj-yi,boxysize)
2993 zj=boxshift(zj-zi,boxzsize)
2995 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2997 sss1=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
2998 if (sss1.eq.0) cycle
2999 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
3001 & sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
3002 sssgrad1=sscagrad(1.0d0/dsqrt(rrij),r_cut_int)
3003 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
3004 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
3005 if (sss.lt.1.0d0) then
3007 e1=fac*fac*aad(itypj,iteli)
3008 e2=fac*bad(itypj,iteli)
3009 if (iabs(j-i) .le. 2) then
3012 evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)*sss1
3015 evdw2=evdw2+evdwij*(1.0d0-sss)*sss1
3016 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
3017 & 'evdw2',i,j,sss,evdwij
3019 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
3022 fac=-(evdwij+e1)*rrij*(1.0d0-sss)*sss1
3023 fac=fac+evdwij*dsqrt(rrij)*(-sssgrad/rscp(itypj,iteli)
3028 C Uncomment following three lines for SC-p interactions
3030 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3032 C Uncomment following line for SC-p interactions
3033 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3035 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
3036 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
3045 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
3046 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
3047 gradx_scp(j,i)=expon*gradx_scp(j,i)
3050 C******************************************************************************
3054 C To save time the factor EXPON has been extracted from ALL components
3055 C of GVDWC and GRADX. Remember to multiply them by this factor before further
3058 C******************************************************************************
3061 C-----------------------------------------------------------------------------
3062 subroutine escp_short(evdw2,evdw2_14)
3064 C This subroutine calculates the excluded-volume interaction energy between
3065 C peptide-group centers and side chains and its gradient in virtual-bond and
3066 C side-chain vectors.
3069 include 'DIMENSIONS'
3070 include 'COMMON.GEO'
3071 include 'COMMON.VAR'
3072 include 'COMMON.LOCAL'
3073 include 'COMMON.CHAIN'
3074 include 'COMMON.DERIV'
3075 include 'COMMON.INTERACT'
3076 include 'COMMON.FFIELD'
3077 include 'COMMON.IOUNITS'
3078 include 'COMMON.CONTROL'
3079 include "COMMON.SPLITELE"
3080 integer xshift,yshift,zshift
3081 logical lprint_short
3082 common /shortcheck/ lprint_short
3083 integer i,iint,j,k,iteli,itypj,subchap
3084 double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
3086 double precision evdw2,evdw2_14,evdwij
3087 double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
3088 & dist_temp, dist_init
3089 double precision ggg(3)
3090 double precision sscale,sscagrad
3091 double precision boxshift
3095 cd print '(a)','Enter ESCP'
3097 c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
3098 c & ' iatscp_e=',iatscp_e
3099 c if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
3101 & write(iout,*) "g_listscp_start_short,g_listscp_end_short",
3102 & g_listscp_start_short,g_listscp_end_short
3103 do ikont=g_listscp_start_short,g_listscp_end_short
3104 i=newcontlistscpi_short(ikont)
3105 j=newcontlistscpj_short(ikont)
3106 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
3108 xi=0.5D0*(c(1,i)+c(1,i+1))
3109 yi=0.5D0*(c(2,i)+c(2,i+1))
3110 zi=0.5D0*(c(3,i)+c(3,i+1))
3111 call to_box(xi,yi,zi)
3114 c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
3115 c & " nscp_gr",nscp_gr(i)
3116 c do iint=1,nscp_gr(i)
3118 c do j=iscpstart(i,iint),iscpend(i,iint)
3119 itypj=iabs(itype(j))
3121 c & write (iout,*) "j",j," itypj",itypj
3122 if (itypj.eq.ntyp1) cycle
3123 C Uncomment following three lines for SC-p interactions
3127 C Uncomment following three lines for Ca-p interactions
3135 call to_box(xj,yj,zj)
3136 xj=boxshift(xj-xi,boxxsize)
3137 yj=boxshift(yj-yi,boxysize)
3138 zj=boxshift(zj-zi,boxzsize)
3139 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
3140 c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
3141 c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
3142 sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
3143 sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),
3145 if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
3146 & " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
3147 c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),
3148 c & " subchap",subchap," sss",sss
3149 if (sss.gt.0.0d0) then
3152 e1=fac*fac*aad(itypj,iteli)
3153 e2=fac*bad(itypj,iteli)
3154 if (iabs(j-i) .le. 2) then
3157 evdw2_14=evdw2_14+(e1+e2)*sss
3160 evdw2=evdw2+evdwij*sss
3161 if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
3162 & 'evdw2',i,j,sss,evdwij
3164 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
3166 fac=-(evdwij+e1)*rrij*sss
3167 fac=fac+evdwij*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)/expon
3171 C Uncomment following three lines for SC-p interactions
3173 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3175 C Uncomment following line for SC-p interactions
3176 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
3178 gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
3179 gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
3188 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
3189 gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i)
3190 gradx_scp(j,i)=expon*gradx_scp(j,i)
3193 C******************************************************************************
3197 C To save time the factor EXPON has been extracted from ALL components
3198 C of GVDWC and GRADX. Remember to multiply them by this factor before further
3201 C******************************************************************************