1 subroutine eelecij_scale(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.SHIELD'
21 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
22 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
23 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
24 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
25 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
26 & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
28 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
30 double precision scal_el /1.0d0/
32 double precision scal_el /0.5d0/
35 C 13-go grudnia roku pamietnego...
36 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
40 cd write (iout,*) "eelecij",i,j
45 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
48 ael6i=ael6(iteli,itelj)
49 ael3i=ael3(iteli,itelj)
60 if (xj.lt.0) xj=xj+boxxsize
62 if (yj.lt.0) yj=yj+boxysize
64 if (zj.lt.0) zj=zj+boxzsize
65 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
73 xj=xj_safe+xshift*boxxsize
74 yj=yj_safe+yshift*boxysize
75 zj=zj_safe+zshift*boxzsize
76 dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
77 if(dist_temp.lt.dist_init) then
87 if (isubchap.eq.1) then
101 c For extracting the short-range part of Evdwpp
102 sss=sscale(rij/rpp(iteli,itelj))
103 sssgrad=sscagrad(rij/rpp(iteli,itelj))
106 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
107 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
108 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
109 fac=cosa-3.0D0*cosb*cosg
111 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
112 if (j.eq.i+2) ev1=scal_el*ev1
117 if (shield_mode.eq.0) then
121 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
123 el1=el1*fac_shield(i)**2*fac_shield(j)**2
124 el2=el2*fac_shield(i)**2*fac_shield(j)**2
126 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
127 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
129 evdw1=evdw1+evdwij*(1.0d0-sss)
130 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
131 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
132 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
133 cd & xmedi,ymedi,zmedi,xj,yj,zj
136 write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
137 write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
141 C Calculate contributions to the Cartesian gradient.
144 facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
145 facel=-3*rrmij*(el1+eesij)
151 * Radial derivatives. First process both termini of the fragment (i,j)
156 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
157 & (shield_mode.gt.0)) then
159 do ilist=1,ishield_list(i)
160 iresshield=shield_list(ilist,i)
162 rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
164 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
166 & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
167 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
168 C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
169 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
170 C if (iresshield.gt.i) then
171 C do ishi=i+1,iresshield-1
172 C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
173 C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
177 C do ishi=iresshield,i
178 C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
179 C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
185 do ilist=1,ishield_list(j)
186 iresshield=shield_list(ilist,j)
188 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
190 gshieldx(k,iresshield)=gshieldx(k,iresshield)+
192 & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
193 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
198 gshieldc(k,i)=gshieldc(k,i)+
199 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
200 gshieldc(k,j)=gshieldc(k,j)+
201 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
202 gshieldc(k,i-1)=gshieldc(k,i-1)+
203 & grad_shield(k,i)*eesij/fac_shield(i)*2.0
204 gshieldc(k,j-1)=gshieldc(k,j-1)+
205 & grad_shield(k,j)*eesij/fac_shield(j)*2.0
212 c gelc(k,i)=gelc(k,i)+ghalf
213 c gelc(k,j)=gelc(k,j)+ghalf
215 c 9/28/08 AL Gradient compotents will be summed only at the end
217 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
218 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
221 * Loop over residues i+1 thru j-1.
225 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
228 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
229 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
230 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
233 c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
234 c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
236 c 9/28/08 AL Gradient compotents will be summed only at the end
238 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
239 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
242 * Loop over residues i+1 thru j-1.
246 cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
250 facvdw=ev1+evdwij*(1.0d0-sss)
253 fac=-3*rrmij*(facvdw+facvdw+facel)
258 * Radial derivatives. First process both termini of the fragment (i,j)
265 c gelc(k,i)=gelc(k,i)+ghalf
266 c gelc(k,j)=gelc(k,j)+ghalf
268 c 9/28/08 AL Gradient compotents will be summed only at the end
270 gelc_long(k,j)=gelc(k,j)+ggg(k)
271 gelc_long(k,i)=gelc(k,i)-ggg(k)
274 * Loop over residues i+1 thru j-1.
278 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
281 c 9/28/08 AL Gradient compotents will be summed only at the end
285 ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
286 ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
287 ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
289 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
290 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
296 ecosa=2.0D0*fac3*fac1+fac4
299 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
300 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
302 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
303 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
305 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
306 cd & (dcosg(k),k=1,3)
308 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
309 & fac_shield(i)**2*fac_shield(j)**2
314 c gelc(k,i)=gelc(k,i)+ghalf
315 c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
316 c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
317 c gelc(k,j)=gelc(k,j)+ghalf
318 c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
319 c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
323 cgrad gelc(l,k)=gelc(l,k)+ggg(l)
328 & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
329 & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
330 & *fac_shield(i)**2*fac_shield(j)**2
333 & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
334 & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
335 & *fac_shield(i)**2*fac_shield(j)**2
336 gelc_long(k,j)=gelc_long(k,j)+ggg(k)
337 gelc_long(k,i)=gelc_long(k,i)-ggg(k)
339 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
340 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
341 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
343 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
344 C energy of a peptide unit is assumed in the form of a second-order
345 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
346 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
347 C are computed for EVERY pair of non-contiguous peptide groups.
349 if (j.lt.nres-1) then
360 muij(kkk)=mu(k,i)*mu(l,j)
363 cd write (iout,*) 'EELEC: i',i,' j',j
364 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
365 cd write(iout,*) 'muij',muij
366 ury=scalar(uy(1,i),erij)
367 urz=scalar(uz(1,i),erij)
368 vry=scalar(uy(1,j),erij)
369 vrz=scalar(uz(1,j),erij)
370 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
371 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
372 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
373 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
374 fac=dsqrt(-ael6i)*r3ij
379 cd write (iout,'(4i5,4f10.5)')
380 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
381 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
382 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
384 cd write (iout,'(4f10.5)')
385 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
386 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
387 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
388 cd write (iout,'(9f10.5/)')
389 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
390 C Derivatives of the elements of A in virtual-bond vectors
391 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
393 uryg(k,1)=scalar(erder(1,k),uy(1,i))
394 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
395 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
396 urzg(k,1)=scalar(erder(1,k),uz(1,i))
397 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
398 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
399 vryg(k,1)=scalar(erder(1,k),uy(1,j))
400 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
401 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
402 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
403 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
404 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
406 C Compute radial contributions to the gradient
424 C Add the contributions coming from er
427 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
428 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
429 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
430 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
433 C Derivatives in DC(i)
434 cgrad ghalf1=0.5d0*agg(k,1)
435 cgrad ghalf2=0.5d0*agg(k,2)
436 cgrad ghalf3=0.5d0*agg(k,3)
437 cgrad ghalf4=0.5d0*agg(k,4)
438 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
439 & -3.0d0*uryg(k,2)*vry)!+ghalf1
440 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
441 & -3.0d0*uryg(k,2)*vrz)!+ghalf2
442 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
443 & -3.0d0*urzg(k,2)*vry)!+ghalf3
444 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
445 & -3.0d0*urzg(k,2)*vrz)!+ghalf4
446 C Derivatives in DC(i+1)
447 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
448 & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
449 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
450 & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
451 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
452 & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
453 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
454 & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
455 C Derivatives in DC(j)
456 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
457 & -3.0d0*vryg(k,2)*ury)!+ghalf1
458 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
459 & -3.0d0*vrzg(k,2)*ury)!+ghalf2
460 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
461 & -3.0d0*vryg(k,2)*urz)!+ghalf3
462 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
463 & -3.0d0*vrzg(k,2)*urz)!+ghalf4
464 C Derivatives in DC(j+1) or DC(nres-1)
465 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
466 & -3.0d0*vryg(k,3)*ury)
467 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
468 & -3.0d0*vrzg(k,3)*ury)
469 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
470 & -3.0d0*vryg(k,3)*urz)
471 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
472 & -3.0d0*vrzg(k,3)*urz)
473 cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
475 cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
489 aggi1(k,l)=-aggi1(k,l)
491 aggj1(k,l)=-aggj1(k,l)
494 if (j.lt.nres-1) then
501 aggi1(k,l)=-aggi1(k,l)
503 aggj1(k,l)=-aggj1(k,l)
515 aggi1(k,l)=-aggi1(k,l)
517 aggj1(k,l)=-aggj1(k,l)
522 IF (wel_loc.gt.0.0d0) THEN
523 C Contribution to the local-electrostatic energy coming from the i-j pair
524 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
526 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
528 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
529 & 'eelloc',i,j,eel_loc_ij
532 if (shield_mode.eq.0) then
539 eel_loc_ij=eel_loc_ij
540 & *fac_shield(i)*fac_shield(j)
541 eel_loc=eel_loc+eel_loc_ij
543 if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
544 & (shield_mode.gt.0)) then
547 do ilist=1,ishield_list(i)
548 iresshield=shield_list(ilist,i)
550 rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
553 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
555 & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
556 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
560 do ilist=1,ishield_list(j)
561 iresshield=shield_list(ilist,j)
563 rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
566 gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
568 & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
569 gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
576 gshieldc_ll(k,i)=gshieldc_ll(k,i)+
577 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
578 gshieldc_ll(k,j)=gshieldc_ll(k,j)+
579 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
580 gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
581 & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
582 gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
583 & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
587 C Partial derivatives in virtual-bond dihedral angles gamma
589 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
590 & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
591 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
592 & *fac_shield(i)*fac_shield(j)
594 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
595 & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
596 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
597 & *fac_shield(i)*fac_shield(j)
599 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
601 ggg(l)=(agg(l,1)*muij(1)+
602 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
603 & *fac_shield(i)*fac_shield(j)
605 gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
606 gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
607 cgrad ghalf=0.5d0*ggg(l)
608 cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
609 cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
613 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
616 C Remaining derivatives of eello
618 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
619 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
620 & *fac_shield(i)*fac_shield(j)
622 gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
623 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
624 & *fac_shield(i)*fac_shield(j)
626 gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
627 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
628 & *fac_shield(i)*fac_shield(j)
630 gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
631 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
632 & *fac_shield(i)*fac_shield(j)
636 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
637 c if (j.gt.i+1 .and. num_conti.le.maxconts) then
638 if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
639 & .and. num_conti.le.maxconts) then
640 c write (iout,*) i,j," entered corr"
642 C Calculate the contact function. The ith column of the array JCONT will
643 C contain the numbers of atoms that make contacts with the atom I (of numbers
644 C greater than I). The arrays FACONT and GACONT will contain the values of
645 C the contact function and its derivative.
646 c r0ij=1.02D0*rpp(iteli,itelj)
647 c r0ij=1.11D0*rpp(iteli,itelj)
648 r0ij=2.20D0*rpp(iteli,itelj)
649 c r0ij=1.55D0*rpp(iteli,itelj)
650 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
651 if (fcont.gt.0.0D0) then
652 num_conti=num_conti+1
653 if (num_conti.gt.maxconts) then
654 write (iout,*) 'WARNING - max. # of contacts exceeded;',
655 & ' will skip next contacts for this conf.'
657 jcont_hb(num_conti,i)=j
658 cd write (iout,*) "i",i," j",j," num_conti",num_conti,
659 cd & " jcont_hb",jcont_hb(num_conti,i)
660 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
661 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
662 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
664 d_cont(num_conti,i)=rij
665 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
666 C --- Electrostatic-interaction matrix ---
667 a_chuj(1,1,num_conti,i)=a22
668 a_chuj(1,2,num_conti,i)=a23
669 a_chuj(2,1,num_conti,i)=a32
670 a_chuj(2,2,num_conti,i)=a33
671 C --- Gradient of rij
673 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
680 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
681 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
682 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
683 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
684 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
689 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
690 C Calculate contact energies
692 wij=cosa-3.0D0*cosb*cosg
695 c fac3=dsqrt(-ael6i)/r0ij**3
696 fac3=dsqrt(-ael6i)*r3ij
697 c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
698 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
699 if (ees0tmp.gt.0) then
700 ees0pij=dsqrt(ees0tmp)
704 c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
705 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
706 if (ees0tmp.gt.0) then
707 ees0mij=dsqrt(ees0tmp)
712 if (shield_mode.eq.0) then
716 ees0plist(num_conti,i)=j
717 C fac_shield(i)=0.4d0
718 C fac_shield(j)=0.6d0
720 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
721 & *fac_shield(i)*fac_shield(j)
722 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
723 & *fac_shield(i)*fac_shield(j)
725 C Diagnostics. Comment out or remove after debugging!
726 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
727 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
728 c ees0m(num_conti,i)=0.0D0
730 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
731 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
732 C Angular derivatives of the contact function
733 ees0pij1=fac3/ees0pij
734 ees0mij1=fac3/ees0mij
735 fac3p=-3.0D0*fac3*rrmij
736 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
737 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
739 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
740 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
741 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
742 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
743 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
744 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
759 facont_hb(num_conti,i)=fcont
760 fprimcont=fprimcont/rij
761 cd facont_hb(num_conti,i)=1.0D0
762 C Following line is for diagnostics.
765 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
766 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
769 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
770 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
772 gggp(1)=gggp(1)+ees0pijp*xj
773 gggp(2)=gggp(2)+ees0pijp*yj
774 gggp(3)=gggp(3)+ees0pijp*zj
775 gggm(1)=gggm(1)+ees0mijp*xj
776 gggm(2)=gggm(2)+ees0mijp*yj
777 gggm(3)=gggm(3)+ees0mijp*zj
778 C Derivatives due to the contact function
779 gacont_hbr(1,num_conti,i)=fprimcont*xj
780 gacont_hbr(2,num_conti,i)=fprimcont*yj
781 gacont_hbr(3,num_conti,i)=fprimcont*zj
784 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
785 c following the change of gradient-summation algorithm.
787 cgrad ghalfp=0.5D0*gggp(k)
788 cgrad ghalfm=0.5D0*gggm(k)
789 gacontp_hb1(k,num_conti,i)=!ghalfp
790 & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
791 & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
792 & *fac_shield(i)*fac_shield(j)
794 gacontp_hb2(k,num_conti,i)=!ghalfp
795 & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
796 & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
797 & *fac_shield(i)*fac_shield(j)
799 gacontp_hb3(k,num_conti,i)=gggp(k)
800 & *fac_shield(i)*fac_shield(j)
802 gacontm_hb1(k,num_conti,i)=!ghalfm
803 & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
804 & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
805 & *fac_shield(i)*fac_shield(j)
807 gacontm_hb2(k,num_conti,i)=!ghalfm
808 & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
809 & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
810 & *fac_shield(i)*fac_shield(j)
812 gacontm_hb3(k,num_conti,i)=gggm(k)
813 & *fac_shield(i)*fac_shield(j)
817 endif ! num_conti.le.maxconts
820 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
824 aggi(l,k)=aggi(l,k)+ghalf
825 aggi1(l,k)=aggi1(l,k)+agg(l,k)
826 aggj(l,k)=aggj(l,k)+ghalf
829 if (j.eq.nres-1 .and. i.lt.j-2) then
832 aggj1(l,k)=aggj1(l,k)+agg(l,k)
837 c t_eelecij=t_eelecij+MPI_Wtime()-time00