1 subroutine eelecij(i,j,ees,evdw1,eel_loc)
2 implicit real*8 (a-h,o-z)
7 include 'COMMON.CONTROL'
8 include 'COMMON.IOUNITS'
11 include 'COMMON.LOCAL'
12 include 'COMMON.CHAIN'
13 include 'COMMON.DERIV'
14 include 'COMMON.INTERACT'
15 include 'COMMON.CONTACTS'
16 include 'COMMON.TORSION'
17 include 'COMMON.VECTORS'
18 include 'COMMON.FFIELD'
19 include 'COMMON.TIME1'
20 include 'COMMON.SPLITELE'
21 include 'COMMON.SHIELD'
22 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
23 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
24 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
25 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
27 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
28 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
30 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
32 double precision scal_el /1.0d0/
34 double precision scal_el /0.5d0/
37 C 13-go grudnia roku pamietnego...
38 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
41 integer xshift,yshift,zshift
43 cd write (iout,*) "eelecij",i,j
47 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
50 ael6i=ael6(iteli,itelj)
51 ael3i=ael3(iteli,itelj)
58 C xj=c(1,j)+0.5D0*dxj-xmedi
59 C yj=c(2,j)+0.5D0*dyj-ymedi
60 C zj=c(3,j)+0.5D0*dzj-zmedi
65 if (xj.lt.0) xj=xj+boxxsize
67 if (yj.lt.0) yj=yj+boxysize
69 if (zj.lt.0) zj=zj+boxzsize
70 if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
71 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
79 xj=xj_safe+xshift*boxxsize
80 yj=yj_safe+yshift*boxysize
81 zj=zj_safe+zshift*boxzsize
82 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
83 if(dist_temp.lt.dist_init) then
93 if (isubchap.eq.1) then
102 C if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC
104 c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
105 c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
106 C Condition for being inside the proper box
107 c if ((xj.gt.((0.5d0)*boxxsize)).or.
108 c & (xj.lt.((-0.5d0)*boxxsize))) then
112 c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
113 c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
114 C Condition for being inside the proper box
115 c if ((yj.gt.((0.5d0)*boxysize)).or.
116 c & (yj.lt.((-0.5d0)*boxysize))) then
120 c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
121 c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
122 C Condition for being inside the proper box
123 c if ((zj.gt.((0.5d0)*boxzsize)).or.
124 c & (zj.lt.((-0.5d0)*boxzsize))) then
127 C endif !endPBC condintion
131 rij=xj*xj+yj*yj+zj*zj
133 sss=sscale(sqrt(rij))
134 sssgrad=sscagrad(sqrt(rij))
135 c if (sss.gt.0.0d0) then
141 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
142 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
143 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
144 fac=cosa-3.0D0*cosb*cosg
146 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
147 if (j.eq.i+2) ev1=scal_el*ev1
152 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
156 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
157 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
158 if (shield_mode.gt.0) then
161 el1=el1*fac_shield(i)**2*fac_shield(j)**2
162 el2=el2*fac_shield(i)**2*fac_shield(j)**2
171 evdw1=evdw1+evdwij*sss
172 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
173 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
174 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
175 cd & xmedi,ymedi,zmedi,xj,yj,zj
178 write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)')
180 &,iteli,itelj,aaa,evdw1,sss
181 write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
182 &fac_shield(i),fac_shield(j)
186 C Calculate contributions to the Cartesian gradient.
189 facvdw=-6*rrmij*(ev1+evdwij)*sss
190 facel=-3*rrmij*(el1+eesij)
197 * Radial derivatives. First process both termini of the fragment (i,j)
202 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
203 & (shield_mode.gt.0)) then
205 do ilist=1,ishield_list(i)
206 iresshield=shield_list(ilist,i)
208 rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
210 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
212 & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
213 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
214 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
215 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
216 C if (iresshield.gt.i) then
217 C do ishi=i+1,iresshield-1
218 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
219 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
223 C do ishi=iresshield,i
224 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
225 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
231 do ilist=1,ishield_list(j)
232 iresshield=shield_list(ilist,j)
234 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
236 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
238 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
239 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
241 C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
242 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
243 C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
244 C if (iresshield.gt.j) then
245 C do ishi=j+1,iresshield-1
246 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
247 C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
251 C do ishi=iresshield,j
252 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
253 C & -grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
260 gshieldc(k,i)=gshieldc(k,i)+
261 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
262 gshieldc(k,j)=gshieldc(k,j)+
263 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
264 gshieldc(k,i-1)=gshieldc(k,i-1)+
265 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
266 gshieldc(k,j-1)=gshieldc(k,j-1)+
267 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
273 c gelc(k,i)=gelc(k,i)+ghalf
274 c gelc(k,j)=gelc(k,j)+ghalf
276 c 9/28/08 AL Gradient compotents will be summed only at the end
277 C print *,"before", gelc_long(1,i), gelc_long(1,j)
279 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
280 C & +grad_shield(k,j)*eesij/fac_shield(j)
281 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
282 C & +grad_shield(k,i)*eesij/fac_shield(i)
283 C gelc_long(k,i-1)=gelc_long(k,i-1)
284 C & +grad_shield(k,i)*eesij/fac_shield(i)
285 C gelc_long(k,j-1)=gelc_long(k,j-1)
286 C & +grad_shield(k,j)*eesij/fac_shield(j)
288 C print *,"bafter", gelc_long(1,i), gelc_long(1,j)
291 * Loop over residues i+1 thru j-1.
295 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
299 ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
300 ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
301 ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
309 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
310 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
312 c 9/28/08 AL Gradient compotents will be summed only at the end
314 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
315 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
318 * Loop over residues i+1 thru j-1.
322 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
327 facvdw=(ev1+evdwij)*sss
330 fac=-3*rrmij*(facvdw+facvdw+facel)
335 * Radial derivatives. First process both termini of the fragment (i,j)
338 C+eesij*grad_shield(1,i)+eesij*grad_shield(1,j)
340 C+eesij*grad_shield(2,i)+eesij*grad_shield(2,j)
342 C+eesij*grad_shield(3,i)+eesij*grad_shield(3,j)
345 c gelc(k,i)=gelc(k,i)+ghalf
346 c gelc(k,j)=gelc(k,j)+ghalf
348 c 9/28/08 AL Gradient compotents will be summed only at the end
350 gelc_long(k,j)=gelc(k,j)+ggg(k)
351 gelc_long(k,i)=gelc(k,i)-ggg(k)
354 * Loop over residues i+1 thru j-1.
358 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
361 c 9/28/08 AL Gradient compotents will be summed only at the end
362 ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
363 ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
364 ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
366 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
367 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
373 ecosa=2.0D0*fac3*fac1+fac4
376 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
377 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
379 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
380 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
382 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
383 cd & (dcosg(k),k=1,3)
385 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
386 & fac_shield(i)**2*fac_shield(j)**2
390 c gelc(k,i)=gelc(k,i)+ghalf
391 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
392 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
393 c gelc(k,j)=gelc(k,j)+ghalf
394 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
395 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
399 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
402 C print *,"before22", gelc_long(1,i), gelc_long(1,j)
405 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
406 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
407 & *fac_shield(i)**2*fac_shield(j)**2
409 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
410 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
411 & *fac_shield(i)**2*fac_shield(j)**2
412 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
413 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
415 C print *,"before33", gelc_long(1,i), gelc_long(1,j)
419 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
420 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
421 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
423 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
424 C energy of a peptide unit is assumed in the form of a second-order
425 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
426 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
427 C are computed for EVERY pair of non-contiguous peptide groups.
430 if (j.lt.nres-1) then
442 muij(kkk)=mu(k,i)*mu(l,j)
443 c write(iout,*) 'mumu=', mu(k,i),mu(l,j),i,j,k,l
445 gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
446 c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
447 gmuij2(kkk)=gUb2(k,i)*mu(l,j)
448 gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
449 c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
450 gmuji2(kkk)=mu(k,i)*gUb2(l,j)
454 cd write (iout,*) 'EELEC: i',i,' j',j
455 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
456 cd write(iout,*) 'muij',muij
457 ury=scalar(uy(1,i),erij)
458 urz=scalar(uz(1,i),erij)
459 vry=scalar(uy(1,j),erij)
460 vrz=scalar(uz(1,j),erij)
461 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
462 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
463 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
464 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
465 fac=dsqrt(-ael6i)*r3ij
470 cd write (iout,'(4i5,4f10.5)')
471 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
472 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
473 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
475 cd write (iout,'(4f10.5)')
476 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
477 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
478 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
479 cd write (iout,'(9f10.5/)')
480 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
481 C Derivatives of the elements of A in virtual-bond vectors
482 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
484 uryg(k,1)=scalar(erder(1,k),uy(1,i))
485 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
486 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
487 urzg(k,1)=scalar(erder(1,k),uz(1,i))
488 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
489 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
490 vryg(k,1)=scalar(erder(1,k),uy(1,j))
491 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
492 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
493 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
494 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
495 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
497 C Compute radial contributions to the gradient
515 C Add the contributions coming from er
518 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
519 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
520 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
521 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
524 C Derivatives in DC(i)
525 cgrad ghalf1=0.5d0*agg(k,1)
526 cgrad ghalf2=0.5d0*agg(k,2)
527 cgrad ghalf3=0.5d0*agg(k,3)
528 cgrad ghalf4=0.5d0*agg(k,4)
529 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
530 & -3.0d0*uryg(k,2)*vry)!+ghalf1
531 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
532 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
533 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
534 & -3.0d0*urzg(k,2)*vry)!+ghalf3
535 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
536 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
537 C Derivatives in DC(i+1)
538 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
539 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
540 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
541 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
542 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
543 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
544 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
545 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
546 C Derivatives in DC(j)
547 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
548 & -3.0d0*vryg(k,2)*ury)!+ghalf1
549 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
550 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
551 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
552 & -3.0d0*vryg(k,2)*urz)!+ghalf3
553 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
554 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
555 C Derivatives in DC(j+1) or DC(nres-1)
556 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
557 & -3.0d0*vryg(k,3)*ury)
558 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
559 & -3.0d0*vrzg(k,3)*ury)
560 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
561 & -3.0d0*vryg(k,3)*urz)
562 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
563 & -3.0d0*vrzg(k,3)*urz)
564 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
566 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
580 aggi1(k,l)=-aggi1(k,l)
582 aggj1(k,l)=-aggj1(k,l)
585 if (j.lt.nres-1) then
592 aggi1(k,l)=-aggi1(k,l)
594 aggj1(k,l)=-aggj1(k,l)
606 aggi1(k,l)=-aggi1(k,l)
608 aggj1(k,l)=-aggj1(k,l)
613 IF (wel_loc.gt.0.0d0) THEN
614 C Contribution to the local-electrostatic energy coming from the i-j pair
615 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
618 write (iout,*) "muij",muij," a22",a22," a23",a23," a32",a32,
620 write (iout,*) "ij",i,j," eel_loc_ij",eel_loc_ij,
623 if (shield_mode.eq.0) then
630 eel_loc_ij=eel_loc_ij
631 & *fac_shield(i)*fac_shield(j)
632 c if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
633 c & 'eelloc',i,j,eel_loc_ij
634 C Now derivative over eel_loc
635 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
636 & (shield_mode.gt.0)) then
639 do ilist=1,ishield_list(i)
640 iresshield=shield_list(ilist,i)
642 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
645 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
647 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
648 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
652 do ilist=1,ishield_list(j)
653 iresshield=shield_list(ilist,j)
655 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
658 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
660 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
661 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
668 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
669 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
670 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
671 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
672 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
673 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
674 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
675 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
680 c write (iout,*) 'i',i,' j',j,itype(i),itype(j),
681 c & ' eel_loc_ij',eel_loc_ij
682 C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
683 C Calculate patrial derivative for theta angle
685 geel_loc_ij=(a22*gmuij1(1)
689 & *fac_shield(i)*fac_shield(j)
690 c write(iout,*) "derivative over thatai"
691 c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
693 gloc(nphi+i,icg)=gloc(nphi+i,icg)+
694 & geel_loc_ij*wel_loc
695 c write(iout,*) "derivative over thatai-1"
696 c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
703 gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
704 & geel_loc_ij*wel_loc
705 & *fac_shield(i)*fac_shield(j)
707 c Derivative over j residue
708 geel_loc_ji=a22*gmuji1(1)
712 c write(iout,*) "derivative over thataj"
713 c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
716 gloc(nphi+j,icg)=gloc(nphi+j,icg)+
717 & geel_loc_ji*wel_loc
718 & *fac_shield(i)*fac_shield(j)
725 c write(iout,*) "derivative over thataj-1"
726 c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
728 gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
729 & geel_loc_ji*wel_loc
730 & *fac_shield(i)*fac_shield(j)
732 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
734 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
735 & 'eelloc',i,j,eel_loc_ij
736 c if (eel_loc_ij.ne.0)
737 c & write (iout,'(a4,2i4,8f9.5)')'chuj',
738 c & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
740 eel_loc=eel_loc+eel_loc_ij
741 C Partial derivatives in virtual-bond dihedral angles gamma
743 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
744 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
745 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
746 & *fac_shield(i)*fac_shield(j)
748 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
749 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
750 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
751 & *fac_shield(i)*fac_shield(j)
752 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
754 ggg(l)=(agg(l,1)*muij(1)+
755 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
756 & *fac_shield(i)*fac_shield(j)
757 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
758 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
759 cgrad ghalf=0.5d0*ggg(l)
760 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
761 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
765 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
768 C Remaining derivatives of eello
770 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
771 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
772 & *fac_shield(i)*fac_shield(j)
774 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
775 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
776 & *fac_shield(i)*fac_shield(j)
778 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
779 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
780 & *fac_shield(i)*fac_shield(j)
782 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
783 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
784 & *fac_shield(i)*fac_shield(j)
788 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
789 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
790 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
791 & .and. num_conti.le.maxconts) then
792 c write (iout,*) i,j," entered corr"
794 C Calculate the contact function. The ith column of the array JCONT will
795 C contain the numbers of atoms that make contacts with the atom I (of numbers
796 C greater than I). The arrays FACONT and GACONT will contain the values of
797 C the contact function and its derivative.
798 c r0ij=1.02D0*rpp(iteli,itelj)
799 c r0ij=1.11D0*rpp(iteli,itelj)
800 r0ij=2.20D0*rpp(iteli,itelj)
801 c r0ij=1.55D0*rpp(iteli,itelj)
802 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
803 if (fcont.gt.0.0D0) then
804 num_conti=num_conti+1
805 if (num_conti.gt.maxconts) then
806 write (iout,*) 'WARNING - max. # of contacts exceeded;',
807 & ' will skip next contacts for this conf.'
809 jcont_hb(num_conti,i)=j
810 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
811 cd & " jcont_hb",jcont_hb(num_conti,i)
812 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
813 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
814 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
816 d_cont(num_conti,i)=rij
817 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
818 C --- Electrostatic-interaction matrix ---
819 a_chuj(1,1,num_conti,i)=a22
820 a_chuj(1,2,num_conti,i)=a23
821 a_chuj(2,1,num_conti,i)=a32
822 a_chuj(2,2,num_conti,i)=a33
823 C --- Gradient of rij
825 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
832 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
833 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
834 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
835 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
836 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
841 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
842 C Calculate contact energies
844 wij=cosa-3.0D0*cosb*cosg
847 c fac3=dsqrt(-ael6i)/r0ij**3
848 fac3=dsqrt(-ael6i)*r3ij
849 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
850 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
851 if (ees0tmp.gt.0) then
852 ees0pij=dsqrt(ees0tmp)
856 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
857 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
858 if (ees0tmp.gt.0) then
859 ees0mij=dsqrt(ees0tmp)
864 if (shield_mode.eq.0) then
868 ees0plist(num_conti,i)=j
869 C fac_shield(i)=0.4d0
870 C fac_shield(j)=0.6d0
872 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
873 & *fac_shield(i)*fac_shield(j)
874 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
875 & *fac_shield(i)*fac_shield(j)
876 C Diagnostics. Comment out or remove after debugging!
877 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
878 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
879 c ees0m(num_conti,i)=0.0D0
881 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
882 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
883 C Angular derivatives of the contact function
884 ees0pij1=fac3/ees0pij
885 ees0mij1=fac3/ees0mij
886 fac3p=-3.0D0*fac3*rrmij
887 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
888 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
890 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
891 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
892 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
893 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
894 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
895 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
910 facont_hb(num_conti,i)=fcont
911 fprimcont=fprimcont/rij
912 cd facont_hb(num_conti,i)=1.0D0
913 C Following line is for diagnostics.
916 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
917 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
920 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
921 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
923 gggp(1)=gggp(1)+ees0pijp*xj
924 gggp(2)=gggp(2)+ees0pijp*yj
925 gggp(3)=gggp(3)+ees0pijp*zj
926 gggm(1)=gggm(1)+ees0mijp*xj
927 gggm(2)=gggm(2)+ees0mijp*yj
928 gggm(3)=gggm(3)+ees0mijp*zj
929 C Derivatives due to the contact function
930 gacont_hbr(1,num_conti,i)=fprimcont*xj
931 gacont_hbr(2,num_conti,i)=fprimcont*yj
932 gacont_hbr(3,num_conti,i)=fprimcont*zj
935 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
936 c following the change of gradient-summation algorithm.
938 cgrad ghalfp=0.5D0*gggp(k)
939 cgrad ghalfm=0.5D0*gggm(k)
940 gacontp_hb1(k,num_conti,i)=!ghalfp
941 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
942 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
943 & *fac_shield(i)*fac_shield(j)
945 gacontp_hb2(k,num_conti,i)=!ghalfp
946 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
947 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
948 & *fac_shield(i)*fac_shield(j)
950 gacontp_hb3(k,num_conti,i)=gggp(k)
951 & *fac_shield(i)*fac_shield(j)
953 gacontm_hb1(k,num_conti,i)=!ghalfm
954 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
955 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
956 & *fac_shield(i)*fac_shield(j)
958 gacontm_hb2(k,num_conti,i)=!ghalfm
959 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
960 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
961 & *fac_shield(i)*fac_shield(j)
963 gacontm_hb3(k,num_conti,i)=gggm(k)
964 & *fac_shield(i)*fac_shield(j)
967 C Diagnostics. Comment out or remove after debugging!
969 cdiag gacontp_hb1(k,num_conti,i)=0.0D0
970 cdiag gacontp_hb2(k,num_conti,i)=0.0D0
971 cdiag gacontp_hb3(k,num_conti,i)=0.0D0
972 cdiag gacontm_hb1(k,num_conti,i)=0.0D0
973 cdiag gacontm_hb2(k,num_conti,i)=0.0D0
974 cdiag gacontm_hb3(k,num_conti,i)=0.0D0
977 endif ! num_conti.le.maxconts
980 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
984 aggi(l,k)=aggi(l,k)+ghalf
985 aggi1(l,k)=aggi1(l,k)+agg(l,k)
986 aggj(l,k)=aggj(l,k)+ghalf
989 if (j.eq.nres-1 .and. i.lt.j-2) then
992 aggj1(l,k)=aggj1(l,k)+agg(l,k)
997 c t_eelecij=t_eelecij+MPI_Wtime()-time00