1 subroutine etotal(energia)
2 implicit real*8 (a-h,o-z)
4 include 'DIMENSIONS.ZSCOPT'
8 cMS$ATTRIBUTES C :: proc_proc
11 include 'COMMON.IOUNITS'
12 double precision energia(0:max_ene),energia1(0:max_ene+1)
18 include 'COMMON.FFIELD'
19 include 'COMMON.DERIV'
20 include 'COMMON.INTERACT'
21 include 'COMMON.SBRIDGE'
22 include 'COMMON.CHAIN'
23 cd write(iout, '(a,i2)')'Calling etotal ipot=',ipot
24 cd print *,'nnt=',nnt,' nct=',nct
26 C Compute the side-chain and electrostatic interaction energy
28 goto (101,102,103,104,105) ipot
29 C Lennard-Jones potential.
31 cd print '(a)','Exit ELJ'
33 C Lennard-Jones-Kihara potential (shifted).
36 C Berne-Pechukas potential (dilated LJ, angular dependence).
39 C Gay-Berne potential (shifted LJ, angular dependence).
42 C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
45 C Calculate electrostatic (H-bonding) energy of the main chain.
47 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
49 C Calculate excluded-volume interaction energy between peptide groups
52 call escp(evdw2,evdw2_14)
54 C Calculate the disulfide-bridge and other energy and the contributions
55 C from other distance constraints.
56 cd print *,'Calling EHPB'
58 cd print *,'EHPB exitted succesfully.'
60 C Calculate the virtual-bond-angle energy.
63 cd print *,'Bend energy finished.'
65 C Calculate the SC local energy.
68 cd print *,'SCLOC energy finished.'
70 C Calculate the virtual-bond torsional energy.
72 cd print *,'nterm=',nterm
73 call etor(etors,edihcnstr)
75 C 6/23/01 Calculate double-torsional energy
79 C 12/1/95 Multi-body terms
83 if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0
84 & .or. wturn6.gt.0.0d0) then
85 c print *,"calling multibody_eello"
86 call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1)
87 c write (*,*) 'n_corr=',n_corr,' n_corr1=',n_corr1
88 c print *,ecorr,ecorr5,ecorr6,eturn6
90 if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then
91 call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
93 C call multibody(ecorr)
97 C scale large componenets
101 eello_turn3_scal=100.0
102 eello_turn4_scal=100.0
114 ecorr5=ecorr5/ecorr5_scal
115 eel_loc=eel_loc/eel_loc_scal
116 eello_turn3=eello_turn3/eello_turn3_scal
117 eello_turn4=eello_turn4/eello_turn4_scal
118 eturn6=eturn6/eturn6_scal
119 ecorr6=ecorr6/ecorr6_scal
121 if (fgprocs.gt.1) then
122 cd call enerprint(evdw,evdw1,evdw2,ees,ebe,escloc,etors,ehpb,
123 cd & edihcnstr,ecorr,eel_loc,eello_turn4,etot)
133 energia(10)=edihcnstr
137 energia(14)=eello_turn3
138 energia(15)=eello_turn4
143 energia1(i)=energia(i)
145 cd write (iout,*) 'BossID=',BossID,' MyGroup=',MyGroup
146 cd write (*,*) 'BossID=',BossID,' MyGroup=',MyGroup
147 cd write (*,*) 'Processor',MyID,' calls MP_REDUCE in ENERGY',
148 cd & ' BossID=',BossID,' MyGroup=',MyGroup
149 call mp_reduce(energia1(1),energia(1),msglen,BossID,d_vadd,
151 cd write (iout,*) 'Processor',MyID,' Reduce finished'
161 edihcnstr=energia(10)
165 eello_turn3=energia(14)
166 eello_turn4=energia(15)
170 c if (MyID.eq.BossID) then
172 etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
173 & +wang*ebe+wtor*etors+wscloc*escloc
174 & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5
175 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
176 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
185 energia(8)=eello_turn3
186 energia(9)=eello_turn4
193 energia(16)=edihcnstr
198 idumm=proc_proc(etot,i)
200 call proc_proc(etot,i)
202 if(i.eq.1)energia(0)=1.0d+99
208 C Sum up the components of the Cartesian gradient.
212 gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+
213 & welec*gelc(j,i)+wstrain*ghpbc(j,i)+
214 & wcorr*gradcorr(j,i)+
215 & wel_loc*gel_loc(j,i)/eel_loc_scal+
216 & wturn3*gcorr3_turn(j,i)/eello_turn3_scal+
217 & wturn4*gcorr4_turn(j,i)/eello_turn4_scal+
218 & wcorr5*gradcorr5(j,i)/ecorr5_scal+
219 & wcorr6*gradcorr6(j,i)/ecorr6_scal+
220 & wturn6*gcorr6_turn(j,i)/eturn6_scal
221 gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
222 & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)
224 cd print '(i3,9(1pe12.4))',i,(gvdwc(k,i),k=1,3),(gelc(k,i),k=1,3),
225 cd & (gradc(k,i),k=1,3)
230 cd write (iout,*) i,g_corr5_loc(i)
231 gloc(i,icg)=gloc(i,icg)+wcorr*gcorr_loc(i)
232 & +wcorr5*g_corr5_loc(i)/ecorr5_scal
233 & +wcorr6*g_corr6_loc(i)/ecorr6_scal
234 & +wturn4*gel_loc_turn4(i)/eello_turn4_scal
235 & +wturn3*gel_loc_turn3(i)/eello_turn3_scal
236 & +wturn6*gel_loc_turn6(i)/eturn6_scal
237 & +wel_loc*gel_loc_loc(i)/eel_loc_scal
240 cd print*,evdw,wsc,evdw2,wscp,ees+evdw1,welec,ebe,wang,
241 cd & escloc,wscloc,etors,wtor,ehpb,wstrain,nss,ebr,etot
242 cd call enerprint(energia(0))
247 C------------------------------------------------------------------------
248 subroutine enerprint(energia)
249 implicit real*8 (a-h,o-z)
251 include 'DIMENSIONS.ZSCOPT'
252 include 'COMMON.IOUNITS'
253 include 'COMMON.FFIELD'
254 include 'COMMON.SBRIDGE'
255 double precision energia(0:max_ene)
264 eello_turn3=energia(8)
265 eello_turn4=energia(9)
266 eello_turn6=energia(10)
272 edihcnstr=energia(16)
273 write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,ebe,wang,
274 & escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain,
276 & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
277 & eello_turn4,wturn4,eello_turn6,wturn6,edihcnstr,ebr*nss,etot
278 10 format (/'Virtual-chain energies:'//
279 & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
280 & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
281 & 'EES= ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p)'/
282 & 'EBE= ',1pE16.6,' WEIGHT=',1pD16.6,' (bending)'/
283 & 'ESC= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/
284 & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/
285 & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/
286 & 'EHBP= ',1pE16.6,' WEIGHT=',1pD16.6,
287 & ' (SS bridges & dist. cnstr.)'/
288 & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
289 & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
290 & 'ECORR6=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
291 & 'EELLO= ',1pE16.6,' WEIGHT=',1pD16.6,' (electrostatic-local)'/
292 & 'ETURN3=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 3rd order)'/
293 & 'ETURN4=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 4th order)'/
294 & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
295 & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
296 & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
297 & 'ETOT= ',1pE16.6,' (total)')
300 C-----------------------------------------------------------------------
303 C This subroutine calculates the interaction energy of nonbonded side chains
304 C assuming the LJ potential of interaction.
306 implicit real*8 (a-h,o-z)
308 include 'DIMENSIONS.ZSCOPT'
309 parameter (accur=1.0d-10)
312 include 'COMMON.LOCAL'
313 include 'COMMON.CHAIN'
314 include 'COMMON.DERIV'
315 include 'COMMON.INTERACT'
316 include 'COMMON.TORSION'
317 include 'COMMON.ENEPS'
318 include 'COMMON.SBRIDGE'
319 include 'COMMON.NAMES'
320 include 'COMMON.IOUNITS'
321 include 'COMMON.CONTACTS'
325 cd print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
328 eneps_temp(j,i)=0.0d0
341 C Calculate SC interaction energy.
344 cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
345 cd & 'iend=',iend(i,iint)
346 do j=istart(i,iint),iend(i,iint)
351 C Change 12/1/95 to calculate four-body interactions
352 rij=xj*xj+yj*yj+zj*zj
354 c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
355 eps0ij=eps(itypi,itypj)
357 e1=fac*fac*aa(itypi,itypj)
358 e2=fac*bb(itypi,itypj)
360 ij=icant(itypi,itypj)
361 eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij)
362 eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij
363 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
364 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
365 cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)')
366 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
367 cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
368 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
372 C Calculate the components of the gradient in DC and X
374 fac=-rrij*(e1+evdwij)
379 gvdwx(k,i)=gvdwx(k,i)-gg(k)
380 gvdwx(k,j)=gvdwx(k,j)+gg(k)
384 gvdwc(l,k)=gvdwc(l,k)+gg(l)
389 C 12/1/95, revised on 5/20/97
391 C Calculate the contact function. The ith column of the array JCONT will
392 C contain the numbers of atoms that make contacts with the atom I (of numbers
393 C greater than I). The arrays FACONT and GACONT will contain the values of
394 C the contact function and its derivative.
396 C Uncomment next line, if the correlation interactions include EVDW explicitly.
397 c if (j.gt.i+1 .and. evdwij.le.0.0D0) then
398 C Uncomment next line, if the correlation interactions are contact function only
399 if (j.gt.i+1.and. eps0ij.gt.0.0D0) then
401 sigij=sigma(itypi,itypj)
402 r0ij=rs0(itypi,itypj)
404 C Check whether the SC's are not too far to make a contact.
407 call gcont(rij,rcut,1.0d0,0.2d0*rcut,fcont,fprimcont)
408 C Add a new contact, if the SC's are close enough, but not too close (r<sigma).
410 if (fcont.gt.0.0D0) then
411 C If the SC-SC distance if close to sigma, apply spline.
412 cAdam call gcont(-rij,-1.03d0*sigij,2.0d0*sigij,1.0d0,
413 cAdam & fcont1,fprimcont1)
414 cAdam fcont1=1.0d0-fcont1
415 cAdam if (fcont1.gt.0.0d0) then
416 cAdam fprimcont=fprimcont*fcont1+fcont*fprimcont1
417 cAdam fcont=fcont*fcont1
419 C Uncomment following 4 lines to have the geometric average of the epsilon0's
420 cga eps0ij=1.0d0/dsqrt(eps0ij)
422 cga gg(k)=gg(k)*eps0ij
424 cga eps0ij=-evdwij*eps0ij
425 C Uncomment for AL's type of SC correlation interactions.
427 num_conti=num_conti+1
429 facont(num_conti,i)=fcont*eps0ij
430 fprimcont=eps0ij*fprimcont/rij
432 cAdam gacont(1,num_conti,i)=-fprimcont*xj+fcont*gg(1)
433 cAdam gacont(2,num_conti,i)=-fprimcont*yj+fcont*gg(2)
434 cAdam gacont(3,num_conti,i)=-fprimcont*zj+fcont*gg(3)
435 C Uncomment following 3 lines for Skolnick's type of SC correlation.
436 gacont(1,num_conti,i)=-fprimcont*xj
437 gacont(2,num_conti,i)=-fprimcont*yj
438 gacont(3,num_conti,i)=-fprimcont*zj
439 cd write (iout,'(2i5,2f10.5)') i,j,rij,facont(num_conti,i)
440 cd write (iout,'(2i3,3f10.5)')
441 cd & i,j,(gacont(kk,num_conti,i),kk=1,3)
447 num_cont(i)=num_conti
452 gvdwc(j,i)=expon*gvdwc(j,i)
453 gvdwx(j,i)=expon*gvdwx(j,i)
457 C******************************************************************************
461 C To save time, the factor of EXPON has been extracted from ALL components
462 C of GVDWC and GRADX. Remember to multiply them by this factor before further
465 C******************************************************************************
468 C-----------------------------------------------------------------------------
469 subroutine eljk(evdw)
471 C This subroutine calculates the interaction energy of nonbonded side chains
472 C assuming the LJK potential of interaction.
474 implicit real*8 (a-h,o-z)
476 include 'DIMENSIONS.ZSCOPT'
479 include 'COMMON.LOCAL'
480 include 'COMMON.CHAIN'
481 include 'COMMON.DERIV'
482 include 'COMMON.INTERACT'
483 include 'COMMON.ENEPS'
484 include 'COMMON.IOUNITS'
485 include 'COMMON.NAMES'
490 c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
493 eneps_temp(j,i)=0.0d0
504 C Calculate SC interaction energy.
507 do j=istart(i,iint),iend(i,iint)
512 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
514 e_augm=augm(itypi,itypj)*fac_augm
517 r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
518 fac=r_shift_inv**expon
519 e1=fac*fac*aa(itypi,itypj)
520 e2=fac*bb(itypi,itypj)
522 ij=icant(itypi,itypj)
523 eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm)
524 & /dabs(eps(itypi,itypj))
525 eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps(itypi,itypj)
526 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
527 cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
528 cd write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
529 cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
530 cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
531 cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
532 cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
536 C Calculate the components of the gradient in DC and X
538 fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
543 gvdwx(k,i)=gvdwx(k,i)-gg(k)
544 gvdwx(k,j)=gvdwx(k,j)+gg(k)
548 gvdwc(l,k)=gvdwc(l,k)+gg(l)
558 gvdwc(j,i)=expon*gvdwc(j,i)
559 gvdwx(j,i)=expon*gvdwx(j,i)
565 C-----------------------------------------------------------------------------
568 C This subroutine calculates the interaction energy of nonbonded side chains
569 C assuming the Berne-Pechukas potential of interaction.
571 implicit real*8 (a-h,o-z)
573 include 'DIMENSIONS.ZSCOPT'
576 include 'COMMON.LOCAL'
577 include 'COMMON.CHAIN'
578 include 'COMMON.DERIV'
579 include 'COMMON.NAMES'
580 include 'COMMON.INTERACT'
581 include 'COMMON.ENEPS'
582 include 'COMMON.IOUNITS'
583 include 'COMMON.CALC'
585 c double precision rrsave(maxdim)
591 eneps_temp(j,i)=0.0d0
595 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
597 c if (icall.eq.0) then
609 dxi=dc_norm(1,nres+i)
610 dyi=dc_norm(2,nres+i)
611 dzi=dc_norm(3,nres+i)
612 dsci_inv=dsc_inv(itypi)
614 C Calculate SC interaction energy.
617 do j=istart(i,iint),iend(i,iint)
620 dscj_inv=dsc_inv(itypj)
621 chi1=chi(itypi,itypj)
622 chi2=chi(itypj,itypi)
629 alf12=0.5D0*(alf1+alf2)
630 C For diagnostics only!!!
643 dxj=dc_norm(1,nres+j)
644 dyj=dc_norm(2,nres+j)
645 dzj=dc_norm(3,nres+j)
646 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
647 cd if (icall.eq.0) then
653 C Calculate the angle-dependent terms of energy & contributions to derivatives.
655 C Calculate whole angle-dependent part of epsilon and contributions
657 fac=(rrij*sigsq)**expon2
658 e1=fac*fac*aa(itypi,itypj)
659 e2=fac*bb(itypi,itypj)
660 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
661 eps2der=evdwij*eps3rt
662 eps3der=evdwij*eps2rt
663 evdwij=evdwij*eps2rt*eps3rt
664 ij=icant(itypi,itypj)
665 aux=eps1*eps2rt**2*eps3rt**2
666 eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux
667 & /dabs(eps(itypi,itypj))
668 eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj)
672 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
673 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
674 cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
675 cd & restyp(itypi),i,restyp(itypj),j,
676 cd & epsi,sigm,chi1,chi2,chip1,chip2,
677 cd & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
678 cd & om1,om2,om12,1.0D0/dsqrt(rrij),
681 C Calculate gradient components.
682 e1=e1*eps1*eps2rt**2*eps3rt**2
683 fac=-expon*(e1+evdwij)
686 C Calculate radial part of the gradient
690 C Calculate the angular part of the gradient and sum add the contributions
691 C to the appropriate components of the Cartesian gradient.
700 C-----------------------------------------------------------------------------
703 C This subroutine calculates the interaction energy of nonbonded side chains
704 C assuming the Gay-Berne potential of interaction.
706 implicit real*8 (a-h,o-z)
708 include 'DIMENSIONS.ZSCOPT'
711 include 'COMMON.LOCAL'
712 include 'COMMON.CHAIN'
713 include 'COMMON.DERIV'
714 include 'COMMON.NAMES'
715 include 'COMMON.INTERACT'
716 include 'COMMON.ENEPS'
717 include 'COMMON.IOUNITS'
718 include 'COMMON.CALC'
725 eneps_temp(j,i)=0.0d0
729 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
732 c if (icall.gt.0) lprn=.true.
740 dxi=dc_norm(1,nres+i)
741 dyi=dc_norm(2,nres+i)
742 dzi=dc_norm(3,nres+i)
743 dsci_inv=dsc_inv(itypi)
745 C Calculate SC interaction energy.
748 do j=istart(i,iint),iend(i,iint)
751 dscj_inv=dsc_inv(itypj)
752 sig0ij=sigma(itypi,itypj)
753 chi1=chi(itypi,itypj)
754 chi2=chi(itypj,itypi)
761 alf12=0.5D0*(alf1+alf2)
762 C For diagnostics only!!!
775 dxj=dc_norm(1,nres+j)
776 dyj=dc_norm(2,nres+j)
777 dzj=dc_norm(3,nres+j)
778 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
780 C Calculate angle-dependent terms of energy and contributions to their
784 sig=sig0ij*dsqrt(sigsq)
785 rij_shift=1.0D0/rij-sig+sig0ij
786 C I hate to put IF's in the loops, but here don't have another choice!!!!
787 if (rij_shift.le.0.0D0) then
792 c---------------------------------------------------------------
793 rij_shift=1.0D0/rij_shift
795 e1=fac*fac*aa(itypi,itypj)
796 e2=fac*bb(itypi,itypj)
797 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
798 eps2der=evdwij*eps3rt
799 eps3der=evdwij*eps2rt
800 evdwij=evdwij*eps2rt*eps3rt
802 ij=icant(itypi,itypj)
803 aux=eps1*eps2rt**2*eps3rt**2
804 eneps_temp(1,ij)=eneps_temp(1,ij)+aux*e1
805 & /dabs(eps(itypi,itypj))
806 eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj)
807 c write (iout,*) "i",i," j",j," itypi",itypi," itypj",itypj,
808 c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)),
809 c & aux*e2/eps(itypi,itypj)
811 c sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
812 c epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
813 c write (iout,'(2(a3,i3,2x),17(0pf7.3))')
814 c & restyp(itypi),i,restyp(itypj),j,
815 c & epsi,sigm,chi1,chi2,chip1,chip2,
816 c & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
817 c & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
821 C Calculate gradient components.
822 e1=e1*eps1*eps2rt**2*eps3rt**2
823 fac=-expon*(e1+evdwij)*rij_shift
826 C Calculate the radial part of the gradient
830 C Calculate angular part of the gradient.
837 C-----------------------------------------------------------------------------
838 subroutine egbv(evdw)
840 C This subroutine calculates the interaction energy of nonbonded side chains
841 C assuming the Gay-Berne-Vorobjev potential of interaction.
843 implicit real*8 (a-h,o-z)
845 include 'DIMENSIONS.ZSCOPT'
848 include 'COMMON.LOCAL'
849 include 'COMMON.CHAIN'
850 include 'COMMON.DERIV'
851 include 'COMMON.NAMES'
852 include 'COMMON.INTERACT'
853 include 'COMMON.ENEPS'
854 include 'COMMON.IOUNITS'
855 include 'COMMON.CALC'
862 eneps_temp(j,i)=0.0d0
866 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
869 c if (icall.gt.0) lprn=.true.
877 dxi=dc_norm(1,nres+i)
878 dyi=dc_norm(2,nres+i)
879 dzi=dc_norm(3,nres+i)
880 dsci_inv=dsc_inv(itypi)
882 C Calculate SC interaction energy.
885 do j=istart(i,iint),iend(i,iint)
888 dscj_inv=dsc_inv(itypj)
889 sig0ij=sigma(itypi,itypj)
891 chi1=chi(itypi,itypj)
892 chi2=chi(itypj,itypi)
899 alf12=0.5D0*(alf1+alf2)
900 C For diagnostics only!!!
913 dxj=dc_norm(1,nres+j)
914 dyj=dc_norm(2,nres+j)
915 dzj=dc_norm(3,nres+j)
916 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
918 C Calculate angle-dependent terms of energy and contributions to their
922 sig=sig0ij*dsqrt(sigsq)
923 rij_shift=1.0D0/rij-sig+r0ij
924 C I hate to put IF's in the loops, but here don't have another choice!!!!
925 if (rij_shift.le.0.0D0) then
930 c---------------------------------------------------------------
931 rij_shift=1.0D0/rij_shift
933 e1=fac*fac*aa(itypi,itypj)
934 e2=fac*bb(itypi,itypj)
935 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
936 eps2der=evdwij*eps3rt
937 eps3der=evdwij*eps2rt
939 e_augm=augm(itypi,itypj)*fac_augm
940 evdwij=evdwij*eps2rt*eps3rt
941 evdw=evdw+evdwij+e_augm
942 ij=icant(itypi,itypj)
943 aux=eps1*eps2rt**2*eps3rt**2
944 eneps_temp(1,ij)=eneps_temp(1,ij)+aux*(e1+e_augm)
945 & /dabs(eps(itypi,itypj))
946 eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj)
947 c eneps_temp(ij)=eneps_temp(ij)
948 c & +(evdwij+e_augm)/eps(itypi,itypj)
950 c sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
951 c epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
952 c write (iout,'(2(a3,i3,2x),17(0pf7.3))')
953 c & restyp(itypi),i,restyp(itypj),j,
954 c & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
955 c & chi1,chi2,chip1,chip2,
956 c & eps1,eps2rt**2,eps3rt**2,
957 c & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
961 C Calculate gradient components.
962 e1=e1*eps1*eps2rt**2*eps3rt**2
963 fac=-expon*(e1+evdwij)*rij_shift
965 fac=rij*fac-2*expon*rrij*e_augm
966 C Calculate the radial part of the gradient
970 C Calculate angular part of the gradient.
977 C-----------------------------------------------------------------------------
978 subroutine sc_angular
979 C Calculate eps1,eps2,eps3,sigma, and parts of their derivatives in om1,om2,
980 C om12. Called by ebp, egb, and egbv.
982 include 'COMMON.CALC'
986 om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
987 om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
988 om12=dxi*dxj+dyi*dyj+dzi*dzj
990 C Calculate eps1(om12) and its derivative in om12
991 faceps1=1.0D0-om12*chiom12
992 faceps1_inv=1.0D0/faceps1
993 eps1=dsqrt(faceps1_inv)
994 C Following variable is eps1*deps1/dom12
995 eps1_om12=faceps1_inv*chiom12
996 C Calculate sigma(om1,om2,om12) and the derivatives of sigma**2 in om1,om2,
1001 facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
1002 sigsq=1.0D0-facsig*faceps1_inv
1003 sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv
1004 sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv
1005 sigsq_om12=-chi12*(om1om2*faceps1-om12*facsig)*faceps1_inv**2
1006 C Calculate eps2 and its derivatives in om1, om2, and om12.
1009 chipom12=chip12*om12
1010 facp=1.0D0-om12*chipom12
1012 facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
1013 C Following variable is the square root of eps2
1014 eps2rt=1.0D0-facp1*facp_inv
1015 C Following three variables are the derivatives of the square root of eps
1016 C in om1, om2, and om12.
1017 eps2rt_om1=-4.0D0*(chipom1-chipom12*om2)*facp_inv
1018 eps2rt_om2=-4.0D0*(chipom2-chipom12*om1)*facp_inv
1019 eps2rt_om12=4.0D0*chip12*(om1om2*facp-om12*facp1)*facp_inv**2
1020 C Evaluate the "asymmetric" factor in the VDW constant, eps3
1021 eps3rt=1.0D0-alf1*om1+alf2*om2-alf12*om12
1022 C Calculate whole angle-dependent part of epsilon and contributions
1023 C to its derivatives
1026 C----------------------------------------------------------------------------
1028 implicit real*8 (a-h,o-z)
1029 include 'DIMENSIONS'
1030 include 'COMMON.CHAIN'
1031 include 'COMMON.DERIV'
1032 include 'COMMON.CALC'
1033 double precision dcosom1(3),dcosom2(3)
1034 eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
1035 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
1036 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
1037 & -2.0D0*alf12*eps3der+sigder*sigsq_om12
1039 dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
1040 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
1043 gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
1046 gvdwx(k,i)=gvdwx(k,i)-gg(k)
1047 & +(eom12*dc_norm(k,nres+j)+eom1*erij(k))*dsci_inv
1048 gvdwx(k,j)=gvdwx(k,j)+gg(k)
1049 & +(eom12*dc_norm(k,nres+i)+eom2*erij(k))*dscj_inv
1052 C Calculate the components of the gradient in DC and X
1056 gvdwc(l,k)=gvdwc(l,k)+gg(l)
1061 c------------------------------------------------------------------------------
1062 subroutine vec_and_deriv
1063 implicit real*8 (a-h,o-z)
1064 include 'DIMENSIONS'
1065 include 'COMMON.IOUNITS'
1066 include 'COMMON.GEO'
1067 include 'COMMON.VAR'
1068 include 'COMMON.LOCAL'
1069 include 'COMMON.CHAIN'
1070 include 'COMMON.VECTORS'
1071 include 'COMMON.DERIV'
1072 include 'COMMON.INTERACT'
1073 dimension uyder(3,3,2),uzder(3,3,2)
1074 C Compute the local reference systems. For reference system (i), the
1075 C X-axis points from CA(i) to CA(i+1), the Y axis is in the
1076 C CA(i)-CA(i+1)-CA(i+2) plane, and the Z axis is perpendicular to this plane.
1078 if (i.eq.nres-1 .or. itel(i+1).eq.0) then
1079 C Case of the last full residue
1080 C Compute the Z-axis
1081 call vecpr(dc_norm(1,i),dc_norm(1,i-1),uz(1,i))
1082 costh=dcos(pi-theta(nres))
1083 fac=1.0d0/dsqrt(1.0d0-costh*costh)
1088 C Compute the derivatives of uz
1090 uzder(2,1,1)=-dc_norm(3,i-1)
1091 uzder(3,1,1)= dc_norm(2,i-1)
1092 uzder(1,2,1)= dc_norm(3,i-1)
1094 uzder(3,2,1)=-dc_norm(1,i-1)
1095 uzder(1,3,1)=-dc_norm(2,i-1)
1096 uzder(2,3,1)= dc_norm(1,i-1)
1099 uzder(2,1,2)= dc_norm(3,i)
1100 uzder(3,1,2)=-dc_norm(2,i)
1101 uzder(1,2,2)=-dc_norm(3,i)
1103 uzder(3,2,2)= dc_norm(1,i)
1104 uzder(1,3,2)= dc_norm(2,i)
1105 uzder(2,3,2)=-dc_norm(1,i)
1108 C Compute the Y-axis
1111 uy(k,i)=fac*(dc_norm(k,i-1)-costh*dc_norm(k,i))
1114 C Compute the derivatives of uy
1117 uyder(k,j,1)=2*dc_norm(k,i-1)*dc_norm(j,i)
1118 & -dc_norm(k,i)*dc_norm(j,i-1)
1119 uyder(k,j,2)=-dc_norm(j,i)*dc_norm(k,i)
1121 uyder(j,j,1)=uyder(j,j,1)-costh
1122 uyder(j,j,2)=1.0d0+uyder(j,j,2)
1127 uygrad(l,k,j,i)=uyder(l,k,j)
1128 uzgrad(l,k,j,i)=uzder(l,k,j)
1132 call unormderiv(uy(1,i),uyder(1,1,1),facy,uygrad(1,1,1,i))
1133 call unormderiv(uy(1,i),uyder(1,1,2),facy,uygrad(1,1,2,i))
1134 call unormderiv(uz(1,i),uzder(1,1,1),fac,uzgrad(1,1,1,i))
1135 call unormderiv(uz(1,i),uzder(1,1,2),fac,uzgrad(1,1,2,i))
1139 C Compute the Z-axis
1140 call vecpr(dc_norm(1,i),dc_norm(1,i+1),uz(1,i))
1141 costh=dcos(pi-theta(i+2))
1142 fac=1.0d0/dsqrt(1.0d0-costh*costh)
1147 C Compute the derivatives of uz
1149 uzder(2,1,1)=-dc_norm(3,i+1)
1150 uzder(3,1,1)= dc_norm(2,i+1)
1151 uzder(1,2,1)= dc_norm(3,i+1)
1153 uzder(3,2,1)=-dc_norm(1,i+1)
1154 uzder(1,3,1)=-dc_norm(2,i+1)
1155 uzder(2,3,1)= dc_norm(1,i+1)
1158 uzder(2,1,2)= dc_norm(3,i)
1159 uzder(3,1,2)=-dc_norm(2,i)
1160 uzder(1,2,2)=-dc_norm(3,i)
1162 uzder(3,2,2)= dc_norm(1,i)
1163 uzder(1,3,2)= dc_norm(2,i)
1164 uzder(2,3,2)=-dc_norm(1,i)
1167 C Compute the Y-axis
1170 uy(k,i)=facy*(dc_norm(k,i+1)-costh*dc_norm(k,i))
1173 C Compute the derivatives of uy
1176 uyder(k,j,1)=2*dc_norm(k,i+1)*dc_norm(j,i)
1177 & -dc_norm(k,i)*dc_norm(j,i+1)
1178 uyder(k,j,2)=-dc_norm(j,i)*dc_norm(k,i)
1180 uyder(j,j,1)=uyder(j,j,1)-costh
1181 uyder(j,j,2)=1.0d0+uyder(j,j,2)
1186 uygrad(l,k,j,i)=uyder(l,k,j)
1187 uzgrad(l,k,j,i)=uzder(l,k,j)
1191 call unormderiv(uy(1,i),uyder(1,1,1),facy,uygrad(1,1,1,i))
1192 call unormderiv(uy(1,i),uyder(1,1,2),facy,uygrad(1,1,2,i))
1193 call unormderiv(uz(1,i),uzder(1,1,1),fac,uzgrad(1,1,1,i))
1194 call unormderiv(uz(1,i),uzder(1,1,2),fac,uzgrad(1,1,2,i))
1203 uygrad(l,k,j,i)=vblinv*uygrad(l,k,j,i)
1204 uzgrad(l,k,j,i)=vblinv*uzgrad(l,k,j,i)
1212 C-----------------------------------------------------------------------------
1213 subroutine vec_and_deriv_test
1214 implicit real*8 (a-h,o-z)
1215 include 'DIMENSIONS'
1216 include 'COMMON.IOUNITS'
1217 include 'COMMON.GEO'
1218 include 'COMMON.VAR'
1219 include 'COMMON.LOCAL'
1220 include 'COMMON.CHAIN'
1221 include 'COMMON.VECTORS'
1222 dimension uyder(3,3,2),uzder(3,3,2)
1223 C Compute the local reference systems. For reference system (i), the
1224 C X-axis points from CA(i) to CA(i+1), the Y axis is in the
1225 C CA(i)-CA(i+1)-CA(i+2) plane, and the Z axis is perpendicular to this plane.
1227 if (i.eq.nres-1) then
1228 C Case of the last full residue
1229 C Compute the Z-axis
1230 call vecpr(dc_norm(1,i),dc_norm(1,i-1),uz(1,i))
1231 costh=dcos(pi-theta(nres))
1232 fac=1.0d0/dsqrt(1.0d0-costh*costh)
1233 c write (iout,*) 'fac',fac,
1234 c & 1.0d0/dsqrt(scalar(uz(1,i),uz(1,i)))
1235 fac=1.0d0/dsqrt(scalar(uz(1,i),uz(1,i)))
1239 C Compute the derivatives of uz
1241 uzder(2,1,1)=-dc_norm(3,i-1)
1242 uzder(3,1,1)= dc_norm(2,i-1)
1243 uzder(1,2,1)= dc_norm(3,i-1)
1245 uzder(3,2,1)=-dc_norm(1,i-1)
1246 uzder(1,3,1)=-dc_norm(2,i-1)
1247 uzder(2,3,1)= dc_norm(1,i-1)
1250 uzder(2,1,2)= dc_norm(3,i)
1251 uzder(3,1,2)=-dc_norm(2,i)
1252 uzder(1,2,2)=-dc_norm(3,i)
1254 uzder(3,2,2)= dc_norm(1,i)
1255 uzder(1,3,2)= dc_norm(2,i)
1256 uzder(2,3,2)=-dc_norm(1,i)
1258 C Compute the Y-axis
1260 uy(k,i)=fac*(dc_norm(k,i-1)-costh*dc_norm(k,i))
1263 facy=1.0d0/dsqrt(scalar(dc_norm(1,i),dc_norm(1,i))*
1264 & (scalar(dc_norm(1,i-1),dc_norm(1,i-1))**2-
1265 & scalar(dc_norm(1,i),dc_norm(1,i-1))**2))
1267 c uy(k,i)=facy*(dc_norm(k,i+1)-costh*dc_norm(k,i))
1270 & dc_norm(k,i-1)*scalar(dc_norm(1,i),dc_norm(1,i))
1271 & -scalar(dc_norm(1,i),dc_norm(1,i-1))*dc_norm(k,i)
1274 c write (iout,*) 'facy',facy,
1275 c & 1.0d0/dsqrt(scalar(uy(1,i),uy(1,i)))
1276 facy=1.0d0/dsqrt(scalar(uy(1,i),uy(1,i)))
1278 uy(k,i)=facy*uy(k,i)
1280 C Compute the derivatives of uy
1283 uyder(k,j,1)=2*dc_norm(k,i-1)*dc_norm(j,i)
1284 & -dc_norm(k,i)*dc_norm(j,i-1)
1285 uyder(k,j,2)=-dc_norm(j,i)*dc_norm(k,i)
1287 c uyder(j,j,1)=uyder(j,j,1)-costh
1288 c uyder(j,j,2)=1.0d0+uyder(j,j,2)
1289 uyder(j,j,1)=uyder(j,j,1)
1290 & -scalar(dc_norm(1,i),dc_norm(1,i-1))
1291 uyder(j,j,2)=scalar(dc_norm(1,i),dc_norm(1,i))
1297 uygrad(l,k,j,i)=uyder(l,k,j)
1298 uzgrad(l,k,j,i)=uzder(l,k,j)
1302 call unormderiv(uy(1,i),uyder(1,1,1),facy,uygrad(1,1,1,i))
1303 call unormderiv(uy(1,i),uyder(1,1,2),facy,uygrad(1,1,2,i))
1304 call unormderiv(uz(1,i),uzder(1,1,1),fac,uzgrad(1,1,1,i))
1305 call unormderiv(uz(1,i),uzder(1,1,2),fac,uzgrad(1,1,2,i))
1308 C Compute the Z-axis
1309 call vecpr(dc_norm(1,i),dc_norm(1,i+1),uz(1,i))
1310 costh=dcos(pi-theta(i+2))
1311 fac=1.0d0/dsqrt(1.0d0-costh*costh)
1312 fac=1.0d0/dsqrt(scalar(uz(1,i),uz(1,i)))
1316 C Compute the derivatives of uz
1318 uzder(2,1,1)=-dc_norm(3,i+1)
1319 uzder(3,1,1)= dc_norm(2,i+1)
1320 uzder(1,2,1)= dc_norm(3,i+1)
1322 uzder(3,2,1)=-dc_norm(1,i+1)
1323 uzder(1,3,1)=-dc_norm(2,i+1)
1324 uzder(2,3,1)= dc_norm(1,i+1)
1327 uzder(2,1,2)= dc_norm(3,i)
1328 uzder(3,1,2)=-dc_norm(2,i)
1329 uzder(1,2,2)=-dc_norm(3,i)
1331 uzder(3,2,2)= dc_norm(1,i)
1332 uzder(1,3,2)= dc_norm(2,i)
1333 uzder(2,3,2)=-dc_norm(1,i)
1335 C Compute the Y-axis
1337 facy=1.0d0/dsqrt(scalar(dc_norm(1,i),dc_norm(1,i))*
1338 & (scalar(dc_norm(1,i+1),dc_norm(1,i+1))**2-
1339 & scalar(dc_norm(1,i),dc_norm(1,i+1))**2))
1341 c uy(k,i)=facy*(dc_norm(k,i+1)-costh*dc_norm(k,i))
1344 & dc_norm(k,i+1)*scalar(dc_norm(1,i),dc_norm(1,i))
1345 & -scalar(dc_norm(1,i),dc_norm(1,i+1))*dc_norm(k,i)
1348 c write (iout,*) 'facy',facy,
1349 c & 1.0d0/dsqrt(scalar(uy(1,i),uy(1,i)))
1350 facy=1.0d0/dsqrt(scalar(uy(1,i),uy(1,i)))
1352 uy(k,i)=facy*uy(k,i)
1354 C Compute the derivatives of uy
1357 uyder(k,j,1)=2*dc_norm(k,i+1)*dc_norm(j,i)
1358 & -dc_norm(k,i)*dc_norm(j,i+1)
1359 uyder(k,j,2)=-dc_norm(j,i)*dc_norm(k,i)
1361 c uyder(j,j,1)=uyder(j,j,1)-costh
1362 c uyder(j,j,2)=1.0d0+uyder(j,j,2)
1363 uyder(j,j,1)=uyder(j,j,1)
1364 & -scalar(dc_norm(1,i),dc_norm(1,i+1))
1365 uyder(j,j,2)=scalar(dc_norm(1,i),dc_norm(1,i))
1371 uygrad(l,k,j,i)=uyder(l,k,j)
1372 uzgrad(l,k,j,i)=uzder(l,k,j)
1376 call unormderiv(uy(1,i),uyder(1,1,1),facy,uygrad(1,1,1,i))
1377 call unormderiv(uy(1,i),uyder(1,1,2),facy,uygrad(1,1,2,i))
1378 call unormderiv(uz(1,i),uzder(1,1,1),fac,uzgrad(1,1,1,i))
1379 call unormderiv(uz(1,i),uzder(1,1,2),fac,uzgrad(1,1,2,i))
1386 uygrad(l,k,j,i)=vblinv*uygrad(l,k,j,i)
1387 uzgrad(l,k,j,i)=vblinv*uzgrad(l,k,j,i)
1394 C-----------------------------------------------------------------------------
1395 subroutine check_vecgrad
1396 implicit real*8 (a-h,o-z)
1397 include 'DIMENSIONS'
1398 include 'COMMON.IOUNITS'
1399 include 'COMMON.GEO'
1400 include 'COMMON.VAR'
1401 include 'COMMON.LOCAL'
1402 include 'COMMON.CHAIN'
1403 include 'COMMON.VECTORS'
1404 dimension uygradt(3,3,2,maxres),uzgradt(3,3,2,maxres)
1405 dimension uyt(3,maxres),uzt(3,maxres)
1406 dimension uygradn(3,3,2),uzgradn(3,3,2),erij(3)
1407 double precision delta /1.0d-7/
1410 crc write(iout,'(2i5,2(3f10.5,5x))') i,1,dc_norm(:,i)
1411 crc write(iout,'(2i5,2(3f10.5,5x))') i,2,uy(:,i)
1412 crc write(iout,'(2i5,2(3f10.5,5x)/)')i,3,uz(:,i)
1413 cd write(iout,'(2i5,2(3f10.5,5x))') i,1,
1414 cd & (dc_norm(if90,i),if90=1,3)
1415 cd write(iout,'(2i5,2(3f10.5,5x))') i,2,(uy(if90,i),if90=1,3)
1416 cd write(iout,'(2i5,2(3f10.5,5x)/)')i,3,(uz(if90,i),if90=1,3)
1417 cd write(iout,'(a)')
1423 uygradt(l,k,j,i)=uygrad(l,k,j,i)
1424 uzgradt(l,k,j,i)=uzgrad(l,k,j,i)
1429 call vec_and_deriv_test
1437 cd write (iout,*) 'i=',i
1439 erij(k)=dc_norm(k,i)
1443 dc_norm(k,i)=erij(k)
1445 dc_norm(j,i)=dc_norm(j,i)+delta
1446 c fac=dsqrt(scalar(dc_norm(1,i),dc_norm(1,i)))
1448 c dc_norm(k,i)=dc_norm(k,i)/fac
1450 c write (iout,*) (dc_norm(k,i),k=1,3)
1451 c write (iout,*) (erij(k),k=1,3)
1452 call vec_and_deriv_test
1454 uygradn(k,j,1)=(uy(k,i)-uyt(k,i))/delta
1455 uygradn(k,j,2)=(uy(k,i-1)-uyt(k,i-1))/delta
1456 uzgradn(k,j,1)=(uz(k,i)-uzt(k,i))/delta
1457 uzgradn(k,j,2)=(uz(k,i-1)-uzt(k,i-1))/delta
1459 c write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)')
1460 c & j,(uzgradt(k,j,1,i),k=1,3),(uzgradn(k,j,1),k=1,3),
1461 c & (uzgradt(k,j,2,i-1),k=1,3),(uzgradn(k,j,2),k=1,3)
1464 dc_norm(k,i)=erij(k)
1467 cd write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)')
1468 cd & k,(uygradt(k,l,1,i),l=1,3),(uygradn(k,l,1),l=1,3),
1469 cd & (uygradt(k,l,2,i-1),l=1,3),(uygradn(k,l,2),l=1,3)
1470 cd write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)')
1471 cd & k,(uzgradt(k,l,1,i),l=1,3),(uzgradn(k,l,1),l=1,3),
1472 cd & (uzgradt(k,l,2,i-1),l=1,3),(uzgradn(k,l,2),l=1,3)
1473 cd write (iout,'(a)')
1478 C--------------------------------------------------------------------------
1479 subroutine set_matrices
1480 implicit real*8 (a-h,o-z)
1481 include 'DIMENSIONS'
1482 include 'DIMENSIONS.ZSCOPT'
1483 include 'COMMON.IOUNITS'
1484 include 'COMMON.GEO'
1485 include 'COMMON.VAR'
1486 include 'COMMON.LOCAL'
1487 include 'COMMON.CHAIN'
1488 include 'COMMON.DERIV'
1489 include 'COMMON.INTERACT'
1490 include 'COMMON.CONTACTS'
1491 include 'COMMON.TORSION'
1492 include 'COMMON.VECTORS'
1493 include 'COMMON.FFIELD'
1494 double precision auxvec(2),auxmat(2,2)
1496 C Compute the virtual-bond-torsional-angle dependent quantities needed
1497 C to calculate the el-loc multibody terms of various order.
1500 if (i .lt. nres+1) then
1537 if (i .gt. 3 .and. i .lt. nres+1) then
1538 obrot_der(1,i-2)=-sin1
1539 obrot_der(2,i-2)= cos1
1540 Ugder(1,1,i-2)= sin1
1541 Ugder(1,2,i-2)=-cos1
1542 Ugder(2,1,i-2)=-cos1
1543 Ugder(2,2,i-2)=-sin1
1546 obrot2_der(1,i-2)=-dwasin2
1547 obrot2_der(2,i-2)= dwacos2
1548 Ug2der(1,1,i-2)= dwasin2
1549 Ug2der(1,2,i-2)=-dwacos2
1550 Ug2der(2,1,i-2)=-dwacos2
1551 Ug2der(2,2,i-2)=-dwasin2
1553 obrot_der(1,i-2)=0.0d0
1554 obrot_der(2,i-2)=0.0d0
1555 Ugder(1,1,i-2)=0.0d0
1556 Ugder(1,2,i-2)=0.0d0
1557 Ugder(2,1,i-2)=0.0d0
1558 Ugder(2,2,i-2)=0.0d0
1559 obrot2_der(1,i-2)=0.0d0
1560 obrot2_der(2,i-2)=0.0d0
1561 Ug2der(1,1,i-2)=0.0d0
1562 Ug2der(1,2,i-2)=0.0d0
1563 Ug2der(2,1,i-2)=0.0d0
1564 Ug2der(2,2,i-2)=0.0d0
1566 if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
1567 iti = itortyp(itype(i-2))
1571 if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
1572 iti1 = itortyp(itype(i-1))
1576 cd write (iout,*) '*******i',i,' iti1',iti
1577 cd write (iout,*) 'b1',b1(:,iti)
1578 cd write (iout,*) 'b2',b2(:,iti)
1579 cd write (iout,*) 'Ug',Ug(:,:,i-2)
1580 if (i .gt. iatel_s+2) then
1581 call matvec2(Ug(1,1,i-2),b2(1,iti),Ub2(1,i-2))
1582 call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2))
1583 call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2))
1584 call matmat2(DD(1,1,iti),Ug(1,1,i-2),DUg(1,1,i-2))
1585 call matmat2(Dtilde(1,1,iti),Ug2(1,1,i-2),DtUg2(1,1,i-2))
1586 call matvec2(Ctilde(1,1,iti1),obrot(1,i-2),Ctobr(1,i-2))
1587 call matvec2(Dtilde(1,1,iti),obrot2(1,i-2),Dtobr2(1,i-2))
1597 DtUg2(l,k,i-2)=0.0d0
1601 call matvec2(Ugder(1,1,i-2),b2(1,iti),Ub2der(1,i-2))
1602 call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2))
1603 call matmat2(CC(1,1,iti1),Ugder(1,1,i-2),CUgder(1,1,i-2))
1604 call matmat2(DD(1,1,iti),Ugder(1,1,i-2),DUgder(1,1,i-2))
1605 call matmat2(Dtilde(1,1,iti),Ug2der(1,1,i-2),DtUg2der(1,1,i-2))
1606 call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2))
1607 call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2))
1609 muder(k,i-2)=Ub2der(k,i-2)
1611 if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
1612 iti1 = itortyp(itype(i-1))
1617 mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
1619 C Vectors and matrices dependent on a single virtual-bond dihedral.
1620 call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1))
1621 call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2))
1622 call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2))
1623 call matvec2(CC(1,1,iti1),Ub2(1,i-2),CUgb2(1,i-2))
1624 call matvec2(CC(1,1,iti1),Ub2der(1,i-2),CUgb2der(1,i-2))
1625 call matmat2(EUg(1,1,i-2),CC(1,1,iti1),EUgC(1,1,i-2))
1626 call matmat2(EUgder(1,1,i-2),CC(1,1,iti1),EUgCder(1,1,i-2))
1627 call matmat2(EUg(1,1,i-2),DD(1,1,iti1),EUgD(1,1,i-2))
1628 call matmat2(EUgder(1,1,i-2),DD(1,1,iti1),EUgDder(1,1,i-2))
1629 cd write (iout,*) 'i',i,' mu ',(mu(k,i-2),k=1,2),
1630 cd & ' mu1',(b1(k,i-2),k=1,2),' mu2',(Ub2(k,i-2),k=1,2)
1632 C Matrices dependent on two consecutive virtual-bond dihedrals.
1633 C The order of matrices is from left to right.
1635 call matmat2(DtUg2(1,1,i-1),EUg(1,1,i),DtUg2EUg(1,1,i))
1636 call matmat2(DtUg2der(1,1,i-1),EUg(1,1,i),DtUg2EUgder(1,1,1,i))
1637 call matmat2(DtUg2(1,1,i-1),EUgder(1,1,i),DtUg2EUgder(1,1,2,i))
1638 call transpose2(DtUg2(1,1,i-1),auxmat(1,1))
1639 call matmat2(auxmat(1,1),EUg(1,1,i),Ug2DtEUg(1,1,i))
1640 call matmat2(auxmat(1,1),EUgder(1,1,i),Ug2DtEUgder(1,1,2,i))
1641 call transpose2(DtUg2der(1,1,i-1),auxmat(1,1))
1642 call matmat2(auxmat(1,1),EUg(1,1,i),Ug2DtEUgder(1,1,1,i))
1645 cd iti = itortyp(itype(i))
1648 cd write (iout,'(2f10.5,5x,2f10.5,5x,2f10.5)')
1649 cd & (EE(j,k,iti),k=1,2),(Ug(j,k,i),k=1,2),(EUg(j,k,i),k=1,2)
1654 C--------------------------------------------------------------------------
1655 subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
1657 C This subroutine calculates the average interaction energy and its gradient
1658 C in the virtual-bond vectors between non-adjacent peptide groups, based on
1659 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
1660 C The potential depends both on the distance of peptide-group centers and on
1661 C the orientation of the CA-CA virtual bonds.
1663 implicit real*8 (a-h,o-z)
1664 include 'DIMENSIONS'
1665 include 'DIMENSIONS.ZSCOPT'
1666 include 'COMMON.CONTROL'
1667 include 'COMMON.IOUNITS'
1668 include 'COMMON.GEO'
1669 include 'COMMON.VAR'
1670 include 'COMMON.LOCAL'
1671 include 'COMMON.CHAIN'
1672 include 'COMMON.DERIV'
1673 include 'COMMON.INTERACT'
1674 include 'COMMON.CONTACTS'
1675 include 'COMMON.TORSION'
1676 include 'COMMON.VECTORS'
1677 include 'COMMON.FFIELD'
1678 dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
1679 & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
1680 double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
1681 & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
1682 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1
1683 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
1684 double precision scal_el /0.5d0/
1686 C 13-go grudnia roku pamietnego...
1687 double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
1688 & 0.0d0,1.0d0,0.0d0,
1689 & 0.0d0,0.0d0,1.0d0/
1690 cd write(iout,*) 'In EELEC'
1692 cd write(iout,*) 'Type',i
1693 cd write(iout,*) 'B1',B1(:,i)
1694 cd write(iout,*) 'B2',B2(:,i)
1695 cd write(iout,*) 'CC',CC(:,:,i)
1696 cd write(iout,*) 'DD',DD(:,:,i)
1697 cd write(iout,*) 'EE',EE(:,:,i)
1699 cd call check_vecgrad
1701 if (icheckgrad.eq.1) then
1703 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
1705 dc_norm(k,i)=dc(k,i)*fac
1707 c write (iout,*) 'i',i,' fac',fac
1710 if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1711 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
1712 & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
1713 cd if (wel_loc.gt.0.0d0) then
1714 if (icheckgrad.eq.1) then
1715 call vec_and_deriv_test
1722 cd write (iout,*) 'i=',i
1724 cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
1727 cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
1728 cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
1741 cd print '(a)','Enter EELEC'
1742 cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
1744 gel_loc_loc(i)=0.0d0
1747 do i=iatel_s,iatel_e
1748 if (itel(i).eq.0) goto 1215
1752 dx_normi=dc_norm(1,i)
1753 dy_normi=dc_norm(2,i)
1754 dz_normi=dc_norm(3,i)
1755 xmedi=c(1,i)+0.5d0*dxi
1756 ymedi=c(2,i)+0.5d0*dyi
1757 zmedi=c(3,i)+0.5d0*dzi
1759 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
1760 do j=ielstart(i),ielend(i)
1761 if (itel(j).eq.0) goto 1216
1765 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1766 aaa=app(iteli,itelj)
1767 bbb=bpp(iteli,itelj)
1768 C Diagnostics only!!!
1774 ael6i=ael6(iteli,itelj)
1775 ael3i=ael3(iteli,itelj)
1779 dx_normj=dc_norm(1,j)
1780 dy_normj=dc_norm(2,j)
1781 dz_normj=dc_norm(3,j)
1782 xj=c(1,j)+0.5D0*dxj-xmedi
1783 yj=c(2,j)+0.5D0*dyj-ymedi
1784 zj=c(3,j)+0.5D0*dzj-zmedi
1785 rij=xj*xj+yj*yj+zj*zj
1791 cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
1792 cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
1793 cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
1794 fac=cosa-3.0D0*cosb*cosg
1796 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
1797 if (j.eq.i+2) ev1=scal_el*ev1
1802 el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
1805 c write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij
1806 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
1807 ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
1810 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
1811 cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
1812 cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
1813 cd & xmedi,ymedi,zmedi,xj,yj,zj
1815 C Calculate contributions to the Cartesian gradient.
1820 fac=-3*rrmij*(facvdw+facvdw+facel)
1826 * Radial derivatives. First process both termini of the fragment (i,j)
1833 gelc(k,i)=gelc(k,i)+ghalf
1834 gelc(k,j)=gelc(k,j)+ghalf
1837 * Loop over residues i+1 thru j-1.
1841 gelc(l,k)=gelc(l,k)+ggg(l)
1847 ecosa=2.0D0*fac3*fac1+fac4
1850 ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
1851 ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
1853 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
1854 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
1856 cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
1857 cd & (dcosg(k),k=1,3)
1859 ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1863 gelc(k,i)=gelc(k,i)+ghalf
1864 & +(ecosa*dc_norm(k,j)+ecosb*erij(k))*vblinv
1865 gelc(k,j)=gelc(k,j)+ghalf
1866 & +(ecosa*dc_norm(k,i)+ecosg*erij(k))*vblinv
1870 gelc(l,k)=gelc(l,k)+ggg(l)
1875 IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
1876 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
1877 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
1879 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
1880 C energy of a peptide unit is assumed in the form of a second-order
1881 C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
1882 C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
1883 C are computed for EVERY pair of non-contiguous peptide groups.
1885 if (j.lt.nres-1) then
1896 muij(kkk)=mu(k,i)*mu(l,j)
1899 cd write (iout,*) 'EELEC: i',i,' j',j
1900 cd write (iout,*) 'j',j,' j1',j1,' j2',j2
1901 cd write(iout,*) 'muij',muij
1902 ury=scalar(uy(1,i),erij)
1903 urz=scalar(uz(1,i),erij)
1904 vry=scalar(uy(1,j),erij)
1905 vrz=scalar(uz(1,j),erij)
1906 a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
1907 a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
1908 a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
1909 a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
1910 C For diagnostics only
1915 fac=dsqrt(-ael6i)*r3ij
1916 cd write (2,*) 'fac=',fac
1917 C For diagnostics only
1923 cd write (iout,'(4i5,4f10.5)')
1924 cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1925 cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
1926 cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') (uy(k,i),k=1,3),
1927 cd & (uz(k,i),k=1,3),(uy(k,j),k=1,3),(uz(k,j),k=1,3)
1928 cd write (iout,'(4f10.5)')
1929 cd & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
1930 cd & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
1931 cd write (iout,'(4f10.5)') ury,urz,vry,vrz
1932 cd write (iout,'(2i3,9f10.5/)') i,j,
1933 cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
1935 C Derivatives of the elements of A in virtual-bond vectors
1936 call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
1943 uryg(k,1)=scalar(erder(1,k),uy(1,i))
1944 uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
1945 uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
1946 urzg(k,1)=scalar(erder(1,k),uz(1,i))
1947 urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
1948 urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
1949 vryg(k,1)=scalar(erder(1,k),uy(1,j))
1950 vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
1951 vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
1952 vrzg(k,1)=scalar(erder(1,k),uz(1,j))
1953 vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
1954 vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
1964 C Compute radial contributions to the gradient
1986 C Add the contributions coming from er
1989 agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
1990 agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
1991 agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
1992 agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
1995 C Derivatives in DC(i)
1996 ghalf1=0.5d0*agg(k,1)
1997 ghalf2=0.5d0*agg(k,2)
1998 ghalf3=0.5d0*agg(k,3)
1999 ghalf4=0.5d0*agg(k,4)
2000 aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
2001 & -3.0d0*uryg(k,2)*vry)+ghalf1
2002 aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
2003 & -3.0d0*uryg(k,2)*vrz)+ghalf2
2004 aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
2005 & -3.0d0*urzg(k,2)*vry)+ghalf3
2006 aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
2007 & -3.0d0*urzg(k,2)*vrz)+ghalf4
2008 C Derivatives in DC(i+1)
2009 aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
2010 & -3.0d0*uryg(k,3)*vry)+agg(k,1)
2011 aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
2012 & -3.0d0*uryg(k,3)*vrz)+agg(k,2)
2013 aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
2014 & -3.0d0*urzg(k,3)*vry)+agg(k,3)
2015 aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
2016 & -3.0d0*urzg(k,3)*vrz)+agg(k,4)
2017 C Derivatives in DC(j)
2018 aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
2019 & -3.0d0*vryg(k,2)*ury)+ghalf1
2020 aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
2021 & -3.0d0*vrzg(k,2)*ury)+ghalf2
2022 aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
2023 & -3.0d0*vryg(k,2)*urz)+ghalf3
2024 aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
2025 & -3.0d0*vrzg(k,2)*urz)+ghalf4
2026 C Derivatives in DC(j+1) or DC(nres-1)
2027 aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
2028 & -3.0d0*vryg(k,3)*ury)
2029 aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
2030 & -3.0d0*vrzg(k,3)*ury)
2031 aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
2032 & -3.0d0*vryg(k,3)*urz)
2033 aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
2034 & -3.0d0*vrzg(k,3)*urz)
2039 C Derivatives in DC(i+1)
2040 cd aggi1(k,1)=agg(k,1)
2041 cd aggi1(k,2)=agg(k,2)
2042 cd aggi1(k,3)=agg(k,3)
2043 cd aggi1(k,4)=agg(k,4)
2044 C Derivatives in DC(j)
2049 C Derivatives in DC(j+1)
2054 if (j.eq.nres-1 .and. i.lt.j-2) then
2056 aggj1(k,l)=aggj1(k,l)+agg(k,l)
2057 cd aggj1(k,l)=agg(k,l)
2063 C Check the loc-el terms by numerical integration
2073 aggi(k,l)=-aggi(k,l)
2074 aggi1(k,l)=-aggi1(k,l)
2075 aggj(k,l)=-aggj(k,l)
2076 aggj1(k,l)=-aggj1(k,l)
2079 if (j.lt.nres-1) then
2085 aggi(k,l)=-aggi(k,l)
2086 aggi1(k,l)=-aggi1(k,l)
2087 aggj(k,l)=-aggj(k,l)
2088 aggj1(k,l)=-aggj1(k,l)
2099 aggi(k,l)=-aggi(k,l)
2100 aggi1(k,l)=-aggi1(k,l)
2101 aggj(k,l)=-aggj(k,l)
2102 aggj1(k,l)=-aggj1(k,l)
2108 IF (wel_loc.gt.0.0d0) THEN
2109 C Contribution to the local-electrostatic energy coming from the i-j pair
2110 eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
2112 cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
2113 cd write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
2114 eel_loc=eel_loc+eel_loc_ij
2115 C Partial derivatives in virtual-bond dihedral angles gamma
2118 & gel_loc_loc(i-1)=gel_loc_loc(i-1)+
2119 & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
2120 & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
2121 gel_loc_loc(j-1)=gel_loc_loc(j-1)+
2122 & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
2123 & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
2124 cd call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij)
2125 cd write(iout,*) 'agg ',agg
2126 cd write(iout,*) 'aggi ',aggi
2127 cd write(iout,*) 'aggi1',aggi1
2128 cd write(iout,*) 'aggj ',aggj
2129 cd write(iout,*) 'aggj1',aggj1
2131 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
2133 ggg(l)=agg(l,1)*muij(1)+
2134 & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
2138 gel_loc(l,k)=gel_loc(l,k)+ggg(l)
2141 C Remaining derivatives of eello
2143 gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
2144 & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
2145 gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
2146 & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
2147 gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
2148 & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
2149 gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
2150 & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
2154 if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
2155 C Contributions from turns
2160 call eturn34(i,j,eello_turn3,eello_turn4)
2162 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
2163 if (j.gt.i+1 .and. num_conti.le.maxconts) then
2165 C Calculate the contact function. The ith column of the array JCONT will
2166 C contain the numbers of atoms that make contacts with the atom I (of numbers
2167 C greater than I). The arrays FACONT and GACONT will contain the values of
2168 C the contact function and its derivative.
2169 c r0ij=1.02D0*rpp(iteli,itelj)
2170 c r0ij=1.11D0*rpp(iteli,itelj)
2171 r0ij=2.20D0*rpp(iteli,itelj)
2172 c r0ij=1.55D0*rpp(iteli,itelj)
2173 call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
2174 if (fcont.gt.0.0D0) then
2175 num_conti=num_conti+1
2176 if (num_conti.gt.maxconts) then
2177 write (iout,*) 'WARNING - max. # of contacts exceeded;',
2178 & ' will skip next contacts for this conf.'
2180 jcont_hb(num_conti,i)=j
2181 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
2182 & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
2183 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
2185 d_cont(num_conti,i)=rij
2186 cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
2187 C --- Electrostatic-interaction matrix ---
2188 a_chuj(1,1,num_conti,i)=a22
2189 a_chuj(1,2,num_conti,i)=a23
2190 a_chuj(2,1,num_conti,i)=a32
2191 a_chuj(2,2,num_conti,i)=a33
2192 C --- Gradient of rij
2194 grij_hb_cont(kkk,num_conti,i)=erij(kkk)
2197 c a_chuj(1,1,num_conti,i)=-0.61d0
2198 c a_chuj(1,2,num_conti,i)= 0.4d0
2199 c a_chuj(2,1,num_conti,i)= 0.65d0
2200 c a_chuj(2,2,num_conti,i)= 0.50d0
2201 c else if (i.eq.2) then
2202 c a_chuj(1,1,num_conti,i)= 0.0d0
2203 c a_chuj(1,2,num_conti,i)= 0.0d0
2204 c a_chuj(2,1,num_conti,i)= 0.0d0
2205 c a_chuj(2,2,num_conti,i)= 0.0d0
2207 C --- and its gradients
2208 cd write (iout,*) 'i',i,' j',j
2210 cd write (iout,*) 'iii 1 kkk',kkk
2211 cd write (iout,*) agg(kkk,:)
2214 cd write (iout,*) 'iii 2 kkk',kkk
2215 cd write (iout,*) aggi(kkk,:)
2218 cd write (iout,*) 'iii 3 kkk',kkk
2219 cd write (iout,*) aggi1(kkk,:)
2222 cd write (iout,*) 'iii 4 kkk',kkk
2223 cd write (iout,*) aggj(kkk,:)
2226 cd write (iout,*) 'iii 5 kkk',kkk
2227 cd write (iout,*) aggj1(kkk,:)
2234 a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
2235 a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
2236 a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
2237 a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
2238 a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
2240 c a_chuj_der(k,l,m,mm,num_conti,i)=0.0d0
2246 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
2247 C Calculate contact energies
2249 wij=cosa-3.0D0*cosb*cosg
2252 c fac3=dsqrt(-ael6i)/r0ij**3
2253 fac3=dsqrt(-ael6i)*r3ij
2254 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
2255 ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
2257 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
2258 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
2259 C Diagnostics. Comment out or remove after debugging!
2260 c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
2261 c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
2262 c ees0m(num_conti,i)=0.0D0
2264 c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
2265 c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
2266 facont_hb(num_conti,i)=fcont
2268 C Angular derivatives of the contact function
2269 ees0pij1=fac3/ees0pij
2270 ees0mij1=fac3/ees0mij
2271 fac3p=-3.0D0*fac3*rrmij
2272 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
2273 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
2275 ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
2276 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
2277 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
2278 ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
2279 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
2280 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
2281 ecosap=ecosa1+ecosa2
2282 ecosbp=ecosb1+ecosb2
2283 ecosgp=ecosg1+ecosg2
2284 ecosam=ecosa1-ecosa2
2285 ecosbm=ecosb1-ecosb2
2286 ecosgm=ecosg1-ecosg2
2295 fprimcont=fprimcont/rij
2296 cd facont_hb(num_conti,i)=1.0D0
2297 C Following line is for diagnostics.
2300 dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
2301 dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
2304 gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
2305 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
2307 gggp(1)=gggp(1)+ees0pijp*xj
2308 gggp(2)=gggp(2)+ees0pijp*yj
2309 gggp(3)=gggp(3)+ees0pijp*zj
2310 gggm(1)=gggm(1)+ees0mijp*xj
2311 gggm(2)=gggm(2)+ees0mijp*yj
2312 gggm(3)=gggm(3)+ees0mijp*zj
2313 C Derivatives due to the contact function
2314 gacont_hbr(1,num_conti,i)=fprimcont*xj
2315 gacont_hbr(2,num_conti,i)=fprimcont*yj
2316 gacont_hbr(3,num_conti,i)=fprimcont*zj
2318 ghalfp=0.5D0*gggp(k)
2319 ghalfm=0.5D0*gggm(k)
2320 gacontp_hb1(k,num_conti,i)=ghalfp
2321 & +(ecosap*dc_norm(k,j)+ecosbp*erij(k))*vblinv
2322 gacontp_hb2(k,num_conti,i)=ghalfp
2323 & +(ecosap*dc_norm(k,i)+ecosgp*erij(k))*vblinv
2324 gacontp_hb3(k,num_conti,i)=gggp(k)
2325 gacontm_hb1(k,num_conti,i)=ghalfm
2326 & +(ecosam*dc_norm(k,j)+ecosbm*erij(k))*vblinv
2327 gacontm_hb2(k,num_conti,i)=ghalfm
2328 & +(ecosam*dc_norm(k,i)+ecosgm*erij(k))*vblinv
2329 gacontm_hb3(k,num_conti,i)=gggm(k)
2332 C Diagnostics. Comment out or remove after debugging!
2334 cdiag gacontp_hb1(k,num_conti,i)=0.0D0
2335 cdiag gacontp_hb2(k,num_conti,i)=0.0D0
2336 cdiag gacontp_hb3(k,num_conti,i)=0.0D0
2337 cdiag gacontm_hb1(k,num_conti,i)=0.0D0
2338 cdiag gacontm_hb2(k,num_conti,i)=0.0D0
2339 cdiag gacontm_hb3(k,num_conti,i)=0.0D0
2342 endif ! num_conti.le.maxconts
2347 num_cont_hb(i)=num_conti
2351 cd write (iout,'(i3,3f10.5,5x,3f10.5)')
2352 cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
2354 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
2355 ccc eel_loc=eel_loc+eello_turn3
2358 C-----------------------------------------------------------------------------
2359 subroutine eturn34(i,j,eello_turn3,eello_turn4)
2360 C Third- and fourth-order contributions from turns
2361 implicit real*8 (a-h,o-z)
2362 include 'DIMENSIONS'
2363 include 'DIMENSIONS.ZSCOPT'
2364 include 'COMMON.IOUNITS'
2365 include 'COMMON.GEO'
2366 include 'COMMON.VAR'
2367 include 'COMMON.LOCAL'
2368 include 'COMMON.CHAIN'
2369 include 'COMMON.DERIV'
2370 include 'COMMON.INTERACT'
2371 include 'COMMON.CONTACTS'
2372 include 'COMMON.TORSION'
2373 include 'COMMON.VECTORS'
2374 include 'COMMON.FFIELD'
2376 double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
2377 & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
2378 & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2)
2379 double precision agg(3,4),aggi(3,4),aggi1(3,4),
2380 & aggj(3,4),aggj1(3,4),a_temp(2,2)
2381 common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
2383 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2385 C Third-order contributions
2392 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2393 cd call checkint_turn3(i,a_temp,eello_turn3_num)
2394 call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1))
2395 call transpose2(auxmat(1,1),auxmat1(1,1))
2396 call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
2397 eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
2398 cd write (2,*) 'i,',i,' j',j,'eello_turn3',
2399 cd & 0.5d0*(pizda(1,1)+pizda(2,2)),
2400 cd & ' eello_turn3_num',4*eello_turn3_num
2402 C Derivatives in gamma(i)
2403 call matmat2(EUgder(1,1,i+1),EUg(1,1,i+2),auxmat2(1,1))
2404 call transpose2(auxmat2(1,1),pizda(1,1))
2405 call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1))
2406 gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2))
2407 C Derivatives in gamma(i+1)
2408 call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1))
2409 call transpose2(auxmat2(1,1),pizda(1,1))
2410 call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1))
2411 gel_loc_turn3(i+1)=gel_loc_turn3(i+1)
2412 & +0.5d0*(pizda(1,1)+pizda(2,2))
2413 C Cartesian derivatives
2415 a_temp(1,1)=aggi(l,1)
2416 a_temp(1,2)=aggi(l,2)
2417 a_temp(2,1)=aggi(l,3)
2418 a_temp(2,2)=aggi(l,4)
2419 call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
2420 gcorr3_turn(l,i)=gcorr3_turn(l,i)
2421 & +0.5d0*(pizda(1,1)+pizda(2,2))
2422 a_temp(1,1)=aggi1(l,1)
2423 a_temp(1,2)=aggi1(l,2)
2424 a_temp(2,1)=aggi1(l,3)
2425 a_temp(2,2)=aggi1(l,4)
2426 call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
2427 gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1)
2428 & +0.5d0*(pizda(1,1)+pizda(2,2))
2429 a_temp(1,1)=aggj(l,1)
2430 a_temp(1,2)=aggj(l,2)
2431 a_temp(2,1)=aggj(l,3)
2432 a_temp(2,2)=aggj(l,4)
2433 call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
2434 gcorr3_turn(l,j)=gcorr3_turn(l,j)
2435 & +0.5d0*(pizda(1,1)+pizda(2,2))
2436 a_temp(1,1)=aggj1(l,1)
2437 a_temp(1,2)=aggj1(l,2)
2438 a_temp(2,1)=aggj1(l,3)
2439 a_temp(2,2)=aggj1(l,4)
2440 call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
2441 gcorr3_turn(l,j1)=gcorr3_turn(l,j1)
2442 & +0.5d0*(pizda(1,1)+pizda(2,2))
2445 else if (j.eq.i+3) then
2446 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2448 C Fourth-order contributions
2456 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2457 cd call checkint_turn4(i,a_temp,eello_turn4_num)
2458 iti1=itortyp(itype(i+1))
2459 iti2=itortyp(itype(i+2))
2460 iti3=itortyp(itype(i+3))
2461 call transpose2(EUg(1,1,i+1),e1t(1,1))
2462 call transpose2(Eug(1,1,i+2),e2t(1,1))
2463 call transpose2(Eug(1,1,i+3),e3t(1,1))
2464 call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
2465 call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
2466 s1=scalar2(b1(1,iti2),auxvec(1))
2467 call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
2468 call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
2469 s2=scalar2(b1(1,iti1),auxvec(1))
2470 call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
2471 call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
2472 s3=0.5d0*(pizda(1,1)+pizda(2,2))
2473 eello_turn4=eello_turn4-(s1+s2+s3)
2474 cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
2475 cd & ' eello_turn4_num',8*eello_turn4_num
2476 C Derivatives in gamma(i)
2478 call transpose2(EUgder(1,1,i+1),e1tder(1,1))
2479 call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
2480 call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1))
2481 s1=scalar2(b1(1,iti2),auxvec(1))
2482 call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1))
2483 s3=0.5d0*(pizda(1,1)+pizda(2,2))
2484 gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3)
2485 C Derivatives in gamma(i+1)
2486 call transpose2(EUgder(1,1,i+2),e2tder(1,1))
2487 call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1))
2488 s2=scalar2(b1(1,iti1),auxvec(1))
2489 call matmat2(ae3(1,1),e2tder(1,1),auxmat(1,1))
2490 call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1))
2491 s3=0.5d0*(pizda(1,1)+pizda(2,2))
2492 gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3)
2493 C Derivatives in gamma(i+2)
2494 call transpose2(EUgder(1,1,i+3),e3tder(1,1))
2495 call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1))
2496 s1=scalar2(b1(1,iti2),auxvec(1))
2497 call matmat2(a_temp(1,1),e3tder(1,1),auxmat(1,1))
2498 call matvec2(auxmat(1,1),Ub2(1,i+2),auxvec(1))
2499 s2=scalar2(b1(1,iti1),auxvec(1))
2500 call matmat2(auxmat(1,1),e2t(1,1),auxmat(1,1))
2501 call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1))
2502 s3=0.5d0*(pizda(1,1)+pizda(2,2))
2503 gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3)
2504 C Cartesian derivatives
2505 C Derivatives of this turn contributions in DC(i+2)
2506 if (j.lt.nres-1) then
2508 a_temp(1,1)=agg(l,1)
2509 a_temp(1,2)=agg(l,2)
2510 a_temp(2,1)=agg(l,3)
2511 a_temp(2,2)=agg(l,4)
2512 call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
2513 call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
2514 s1=scalar2(b1(1,iti2),auxvec(1))
2515 call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
2516 call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
2517 s2=scalar2(b1(1,iti1),auxvec(1))
2518 call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
2519 call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
2520 s3=0.5d0*(pizda(1,1)+pizda(2,2))
2522 gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3)
2525 C Remaining derivatives of this turn contribution
2527 a_temp(1,1)=aggi(l,1)
2528 a_temp(1,2)=aggi(l,2)
2529 a_temp(2,1)=aggi(l,3)
2530 a_temp(2,2)=aggi(l,4)
2531 call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
2532 call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
2533 s1=scalar2(b1(1,iti2),auxvec(1))
2534 call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
2535 call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
2536 s2=scalar2(b1(1,iti1),auxvec(1))
2537 call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
2538 call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
2539 s3=0.5d0*(pizda(1,1)+pizda(2,2))
2540 gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3)
2541 a_temp(1,1)=aggi1(l,1)
2542 a_temp(1,2)=aggi1(l,2)
2543 a_temp(2,1)=aggi1(l,3)
2544 a_temp(2,2)=aggi1(l,4)
2545 call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
2546 call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
2547 s1=scalar2(b1(1,iti2),auxvec(1))
2548 call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
2549 call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
2550 s2=scalar2(b1(1,iti1),auxvec(1))
2551 call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
2552 call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
2553 s3=0.5d0*(pizda(1,1)+pizda(2,2))
2554 gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3)
2555 a_temp(1,1)=aggj(l,1)
2556 a_temp(1,2)=aggj(l,2)
2557 a_temp(2,1)=aggj(l,3)
2558 a_temp(2,2)=aggj(l,4)
2559 call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
2560 call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
2561 s1=scalar2(b1(1,iti2),auxvec(1))
2562 call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
2563 call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
2564 s2=scalar2(b1(1,iti1),auxvec(1))
2565 call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
2566 call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
2567 s3=0.5d0*(pizda(1,1)+pizda(2,2))
2568 gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3)
2569 a_temp(1,1)=aggj1(l,1)
2570 a_temp(1,2)=aggj1(l,2)
2571 a_temp(2,1)=aggj1(l,3)
2572 a_temp(2,2)=aggj1(l,4)
2573 call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
2574 call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
2575 s1=scalar2(b1(1,iti2),auxvec(1))
2576 call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
2577 call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
2578 s2=scalar2(b1(1,iti1),auxvec(1))
2579 call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
2580 call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
2581 s3=0.5d0*(pizda(1,1)+pizda(2,2))
2582 gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
2588 C-----------------------------------------------------------------------------
2589 subroutine vecpr(u,v,w)
2590 implicit real*8(a-h,o-z)
2591 dimension u(3),v(3),w(3)
2592 w(1)=u(2)*v(3)-u(3)*v(2)
2593 w(2)=-u(1)*v(3)+u(3)*v(1)
2594 w(3)=u(1)*v(2)-u(2)*v(1)
2597 C-----------------------------------------------------------------------------
2598 subroutine unormderiv(u,ugrad,unorm,ungrad)
2599 C This subroutine computes the derivatives of a normalized vector u, given
2600 C the derivatives computed without normalization conditions, ugrad. Returns
2603 double precision u(3),ugrad(3,3),unorm,ungrad(3,3)
2604 double precision vec(3)
2605 double precision scalar
2607 c write (2,*) 'ugrad',ugrad
2610 vec(i)=scalar(ugrad(1,i),u(1))
2612 c write (2,*) 'vec',vec
2615 ungrad(j,i)=(ugrad(j,i)-u(j)*vec(i))*unorm
2618 c write (2,*) 'ungrad',ungrad
2621 C-----------------------------------------------------------------------------
2622 subroutine escp(evdw2,evdw2_14)
2624 C This subroutine calculates the excluded-volume interaction energy between
2625 C peptide-group centers and side chains and its gradient in virtual-bond and
2626 C side-chain vectors.
2628 implicit real*8 (a-h,o-z)
2629 include 'DIMENSIONS'
2630 include 'DIMENSIONS.ZSCOPT'
2631 include 'COMMON.GEO'
2632 include 'COMMON.VAR'
2633 include 'COMMON.LOCAL'
2634 include 'COMMON.CHAIN'
2635 include 'COMMON.DERIV'
2636 include 'COMMON.INTERACT'
2637 include 'COMMON.FFIELD'
2638 include 'COMMON.IOUNITS'
2642 cd print '(a)','Enter ESCP'
2643 c write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e,
2644 c & ' scal14',scal14
2645 do i=iatscp_s,iatscp_e
2647 c write (iout,*) "i",i," iteli",iteli," nscp_gr",nscp_gr(i),
2648 c & " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
2649 if (iteli.eq.0) goto 1225
2650 xi=0.5D0*(c(1,i)+c(1,i+1))
2651 yi=0.5D0*(c(2,i)+c(2,i+1))
2652 zi=0.5D0*(c(3,i)+c(3,i+1))
2654 do iint=1,nscp_gr(i)
2656 do j=iscpstart(i,iint),iscpend(i,iint)
2658 C Uncomment following three lines for SC-p interactions
2662 C Uncomment following three lines for Ca-p interactions
2666 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
2668 e1=fac*fac*aad(itypj,iteli)
2669 e2=fac*bad(itypj,iteli)
2670 if (iabs(j-i) .le. 2) then
2673 evdw2_14=evdw2_14+e1+e2
2676 c write (iout,*) i,j,evdwij
2680 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
2682 fac=-(evdwij+e1)*rrij
2687 cd write (iout,*) 'j<i'
2688 C Uncomment following three lines for SC-p interactions
2690 c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
2693 cd write (iout,*) 'j>i'
2696 C Uncomment following line for SC-p interactions
2697 c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
2701 gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
2705 cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
2706 cd write (iout,*) ggg(1),ggg(2),ggg(3)
2709 gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
2719 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
2720 gradx_scp(j,i)=expon*gradx_scp(j,i)
2723 C******************************************************************************
2727 C To save time the factor EXPON has been extracted from ALL components
2728 C of GVDWC and GRADX. Remember to multiply them by this factor before further
2731 C******************************************************************************
2734 C--------------------------------------------------------------------------
2735 subroutine edis(ehpb)
2737 C Evaluate bridge-strain energy and its gradient in virtual-bond and SC vectors.
2739 implicit real*8 (a-h,o-z)
2740 include 'DIMENSIONS'
2741 include 'COMMON.SBRIDGE'
2742 include 'COMMON.CHAIN'
2743 include 'COMMON.DERIV'
2744 include 'COMMON.VAR'
2747 cd print *,'edis: nhpb=',nhpb,' fbr=',fbr
2748 cd print *,'link_start=',link_start,' link_end=',link_end
2749 if (link_end.eq.0) return
2750 do i=link_start,link_end
2751 C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
2752 C CA-CA distance used in regularization of structure.
2755 C iii and jjj point to the residues for which the distance is assigned.
2756 if (ii.gt.nres) then
2763 C Calculate the distance between the two points and its difference from the
2767 C Get the force constant corresponding to this distance.
2769 C Calculate the contribution to energy.
2770 ehpb=ehpb+waga*rdis*rdis
2772 C Evaluate gradient.
2775 cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
2776 cd & ' waga=',waga,' fac=',fac
2778 ggg(j)=fac*(c(j,jj)-c(j,ii))
2780 cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
2781 C If this is a SC-SC distace, we need to calculate the contributions to the
2782 C Cartesian gradient in the SC vectors (ghpbx).
2785 ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
2786 ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
2791 ghpbc(k,j)=ghpbc(k,j)+ggg(k)
2798 C--------------------------------------------------------------------------
2799 subroutine ebend(etheta)
2801 C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
2802 C angles gamma and its derivatives in consecutive thetas and gammas.
2804 implicit real*8 (a-h,o-z)
2805 include 'DIMENSIONS'
2806 include 'DIMENSIONS.ZSCOPT'
2807 include 'COMMON.LOCAL'
2808 include 'COMMON.GEO'
2809 include 'COMMON.INTERACT'
2810 include 'COMMON.DERIV'
2811 include 'COMMON.VAR'
2812 include 'COMMON.CHAIN'
2813 include 'COMMON.IOUNITS'
2814 include 'COMMON.NAMES'
2815 include 'COMMON.FFIELD'
2816 common /calcthet/ term1,term2,termm,diffak,ratak,
2817 & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,
2818 & delthe0,sig0inv,sigtc,sigsqtc,delthec,it
2819 double precision y(2),z(2)
2821 time11=dexp(-2*time)
2824 c write (iout,*) "nres",nres
2825 c write (*,'(a,i2)') 'EBEND ICG=',icg
2826 c write (iout,*) ithet_start,ithet_end
2827 do i=ithet_start,ithet_end
2828 C Zero the energy function and its derivative at 0 or pi.
2829 call splinthet(theta(i),0.5d0*delta,ss,ssd)
2831 if (i.gt.ithet_start .and.
2832 & (itel(i-1).eq.0 .or. itel(i-2).eq.0)) goto 1215
2833 if (i.gt.3 .and. (i.le.4 .or. itel(i-3).ne.0)) then
2841 if (i.lt.nres .and. itel(i).ne.0) then
2849 C Calculate the "mean" value of theta from the part of the distribution
2850 C dependent on the adjacent virtual-bond-valence angles (gamma1 & gamma2).
2851 C In following comments this theta will be referred to as t_c.
2852 thet_pred_mean=0.0d0
2856 thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
2858 c write (iout,*) "thet_pred_mean",thet_pred_mean
2859 dthett=thet_pred_mean*ssd
2860 thet_pred_mean=thet_pred_mean*ss+a0thet(it)
2861 c write (iout,*) "thet_pred_mean",thet_pred_mean
2862 C Derivatives of the "mean" values in gamma1 and gamma2.
2863 dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss
2864 dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss
2865 if (theta(i).gt.pi-delta) then
2866 call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0,
2868 call mixder(pi-delta,thet_pred_mean,theta0(it),fprim_tc0)
2869 call theteng(pi,thet_pred_mean,theta0(it),f1,fprim1,E_tc1)
2870 call spline1(theta(i),pi-delta,delta,f0,f1,fprim0,ethetai,
2872 call spline2(theta(i),pi-delta,delta,E_tc0,E_tc1,fprim_tc0,
2874 else if (theta(i).lt.delta) then
2875 call theteng(delta,thet_pred_mean,theta0(it),f0,fprim0,E_tc0)
2876 call theteng(0.0d0,thet_pred_mean,theta0(it),f1,fprim1,E_tc1)
2877 call spline1(theta(i),delta,-delta,f0,f1,fprim0,ethetai,
2879 call mixder(delta,thet_pred_mean,theta0(it),fprim_tc0)
2880 call spline2(theta(i),delta,-delta,E_tc0,E_tc1,fprim_tc0,
2883 call theteng(theta(i),thet_pred_mean,theta0(it),ethetai,
2886 etheta=etheta+ethetai
2887 c write (iout,'(2i3,3f8.3,f10.5)') i,it,rad2deg*theta(i),
2888 c & rad2deg*phii,rad2deg*phii1,ethetai
2889 if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
2890 if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
2891 gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
2894 C Ufff.... We've done all this!!!
2897 C---------------------------------------------------------------------------
2898 subroutine theteng(thetai,thet_pred_mean,theta0i,ethetai,E_theta,
2900 implicit real*8 (a-h,o-z)
2901 include 'DIMENSIONS'
2902 include 'COMMON.LOCAL'
2903 include 'COMMON.IOUNITS'
2904 common /calcthet/ term1,term2,termm,diffak,ratak,
2905 & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,
2906 & delthe0,sig0inv,sigtc,sigsqtc,delthec,it
2907 C Calculate the contributions to both Gaussian lobes.
2908 C 6/6/97 - Deform the Gaussians using the factor of 1/(1+time)
2909 C The "polynomial part" of the "standard deviation" of this part of
2913 sig=sig*thet_pred_mean+polthet(j,it)
2915 C Derivative of the "interior part" of the "standard deviation of the"
2916 C gamma-dependent Gaussian lobe in t_c.
2917 sigtc=3*polthet(3,it)
2919 sigtc=sigtc*thet_pred_mean+j*polthet(j,it)
2922 C Set the parameters of both Gaussian lobes of the distribution.
2923 C "Standard deviation" of the gamma-dependent Gaussian lobe (sigtc)
2924 fac=sig*sig+sigc0(it)
2927 C Following variable (sigsqtc) is -(1/2)d[sigma(t_c)**(-2))]/dt_c
2928 sigsqtc=-4.0D0*sigcsq*sigtc
2929 c print *,i,sig,sigtc,sigsqtc
2930 C Following variable (sigtc) is d[sigma(t_c)]/dt_c
2931 sigtc=-sigtc/(fac*fac)
2932 C Following variable is sigma(t_c)**(-2)
2933 sigcsq=sigcsq*sigcsq
2935 sig0inv=1.0D0/sig0i**2
2936 delthec=thetai-thet_pred_mean
2937 delthe0=thetai-theta0i
2938 term1=-0.5D0*sigcsq*delthec*delthec
2939 term2=-0.5D0*sig0inv*delthe0*delthe0
2940 C Following fuzzy logic is to avoid underflows in dexp and subsequent INFs and
2941 C NaNs in taking the logarithm. We extract the largest exponent which is added
2942 C to the energy (this being the log of the distribution) at the end of energy
2943 C term evaluation for this virtual-bond angle.
2944 if (term1.gt.term2) then
2946 term2=dexp(term2-termm)
2950 term1=dexp(term1-termm)
2953 C The ratio between the gamma-independent and gamma-dependent lobes of
2954 C the distribution is a Gaussian function of thet_pred_mean too.
2955 diffak=gthet(2,it)-thet_pred_mean
2956 ratak=diffak/gthet(3,it)**2
2957 ak=dexp(gthet(1,it)-0.5D0*diffak*ratak)
2958 C Let's differentiate it in thet_pred_mean NOW.
2960 C Now put together the distribution terms to make complete distribution.
2961 termexp=term1+ak*term2
2962 termpre=sigc+ak*sig0i
2963 C Contribution of the bending energy from this theta is just the -log of
2964 C the sum of the contributions from the two lobes and the pre-exponential
2965 C factor. Simple enough, isn't it?
2966 ethetai=(-dlog(termexp)-termm+dlog(termpre))
2967 C NOW the derivatives!!!
2968 C 6/6/97 Take into account the deformation.
2969 E_theta=(delthec*sigcsq*term1
2970 & +ak*delthe0*sig0inv*term2)/termexp
2971 E_tc=((sigtc+aktc*sig0i)/termpre
2972 & -((delthec*sigcsq+delthec*delthec*sigsqtc)*term1+
2973 & aktc*term2)/termexp)
2976 c-----------------------------------------------------------------------------
2977 subroutine mixder(thetai,thet_pred_mean,theta0i,E_tc_t)
2978 implicit real*8 (a-h,o-z)
2979 include 'DIMENSIONS'
2980 include 'COMMON.LOCAL'
2981 include 'COMMON.IOUNITS'
2982 common /calcthet/ term1,term2,termm,diffak,ratak,
2983 & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,
2984 & delthe0,sig0inv,sigtc,sigsqtc,delthec,it
2985 delthec=thetai-thet_pred_mean
2986 delthe0=thetai-theta0i
2987 C "Thank you" to MAPLE (probably spared one day of hand-differentiation).
2988 t3 = thetai-thet_pred_mean
2992 t14 = t12+t6*sigsqtc
2994 t21 = thetai-theta0i
3000 E_tc_t = -((sigcsq+2.D0*t3*sigsqtc)*t9-t14*sigcsq*t3*t16*t9
3001 & -aktc*sig0inv*t27)/t32+(t14*t9+aktc*t26)/t40
3002 & *(-t12*t9-ak*sig0inv*t27)
3005 c-----------------------------------------------------------------------------
3006 subroutine esc(escloc)
3007 C Calculate the local energy of a side chain and its derivatives in the
3008 C corresponding virtual-bond valence angles THETA and the spherical angles
3010 implicit real*8 (a-h,o-z)
3011 include 'DIMENSIONS'
3012 include 'DIMENSIONS.ZSCOPT'
3013 include 'COMMON.GEO'
3014 include 'COMMON.LOCAL'
3015 include 'COMMON.VAR'
3016 include 'COMMON.INTERACT'
3017 include 'COMMON.DERIV'
3018 include 'COMMON.CHAIN'
3019 include 'COMMON.IOUNITS'
3020 include 'COMMON.NAMES'
3021 include 'COMMON.FFIELD'
3022 double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
3023 & ddersc0(3),ddummy(3),xtemp(3),temp(3)
3024 common /sccalc/ time11,time12,time112,theti,it,nlobit
3027 c write (iout,'(a)') 'ESC'
3028 do i=loc_start,loc_end
3030 if (it.eq.10) goto 1
3032 c print *,'i=',i,' it=',it,' nlobit=',nlobit
3033 c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
3034 theti=theta(i+1)-pipol
3039 if (x(2).gt.pi-delta) then
3043 call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
3045 call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
3046 call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
3048 call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
3049 & ddersc0(1),dersc(1))
3050 call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
3051 & ddersc0(3),dersc(3))
3053 call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
3055 call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
3056 call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
3057 & dersc0(2),esclocbi,dersc02)
3058 call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
3060 call splinthet(x(2),0.5d0*delta,ss,ssd)
3065 dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
3067 dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
3068 c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
3070 escloci=ss*escloci+(1.0d0-ss)*esclocbi
3072 c write (iout,*) escloci
3073 else if (x(2).lt.delta) then
3077 call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
3079 call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
3080 call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
3082 call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
3083 & ddersc0(1),dersc(1))
3084 call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
3085 & ddersc0(3),dersc(3))
3087 call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
3089 call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
3090 call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
3091 & dersc0(2),esclocbi,dersc02)
3092 call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
3097 call splinthet(x(2),0.5d0*delta,ss,ssd)
3099 dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
3101 dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
3102 c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
3104 escloci=ss*escloci+(1.0d0-ss)*esclocbi
3105 c write (iout,*) escloci
3107 call enesc(x,escloci,dersc,ddummy,.false.)
3110 escloc=escloc+escloci
3111 c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
3113 gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
3115 gloc(ialph(i,1),icg)=wscloc*dersc(2)
3116 gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
3121 C---------------------------------------------------------------------------
3122 subroutine enesc(x,escloci,dersc,ddersc,mixed)
3123 implicit real*8 (a-h,o-z)
3124 include 'DIMENSIONS'
3125 include 'COMMON.GEO'
3126 include 'COMMON.LOCAL'
3127 include 'COMMON.IOUNITS'
3128 common /sccalc/ time11,time12,time112,theti,it,nlobit
3129 double precision x(3),z(3),Ax(3,maxlob,-1:1),dersc(3),ddersc(3)
3130 double precision contr(maxlob,-1:1)
3132 c write (iout,*) 'it=',it,' nlobit=',nlobit
3136 if (mixed) ddersc(j)=0.0d0
3140 C Because of periodicity of the dependence of the SC energy in omega we have
3141 C to add up the contributions from x(3)-2*pi, x(3), and x(3+2*pi).
3142 C To avoid underflows, first compute & store the exponents.
3150 z(k)=x(k)-censc(k,j,it)
3155 Axk=Axk+gaussc(l,k,j,it)*z(l)
3161 expfac=expfac+Ax(k,j,iii)*z(k)
3169 C As in the case of ebend, we want to avoid underflows in exponentiation and
3170 C subsequent NaNs and INFs in energy calculation.
3171 C Find the largest exponent
3175 if (emin.gt.contr(j,iii)) emin=contr(j,iii)
3179 cd print *,'it=',it,' emin=',emin
3181 C Compute the contribution to SC energy and derivatives
3185 expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin)
3186 cd print *,'j=',j,' expfac=',expfac
3187 escloc_i=escloc_i+expfac
3189 dersc(k)=dersc(k)+Ax(k,j,iii)*expfac
3193 ddersc(k)=ddersc(k)+(-Ax(2,j,iii)*Ax(k,j,iii)
3194 & +gaussc(k,2,j,it))*expfac
3201 dersc(1)=dersc(1)/cos(theti)**2
3202 ddersc(1)=ddersc(1)/cos(theti)**2
3205 escloci=-(dlog(escloc_i)-emin)
3207 dersc(j)=dersc(j)/escloc_i
3211 ddersc(j)=(ddersc(j)/escloc_i+dersc(2)*dersc(j))
3216 C------------------------------------------------------------------------------
3217 subroutine enesc_bound(x,escloci,dersc,dersc12,mixed)
3218 implicit real*8 (a-h,o-z)
3219 include 'DIMENSIONS'
3220 include 'COMMON.GEO'
3221 include 'COMMON.LOCAL'
3222 include 'COMMON.IOUNITS'
3223 common /sccalc/ time11,time12,time112,theti,it,nlobit
3224 double precision x(3),z(3),Ax(3,maxlob),dersc(3)
3225 double precision contr(maxlob)
3236 z(k)=x(k)-censc(k,j,it)
3242 Axk=Axk+gaussc(l,k,j,it)*z(l)
3248 expfac=expfac+Ax(k,j)*z(k)
3253 C As in the case of ebend, we want to avoid underflows in exponentiation and
3254 C subsequent NaNs and INFs in energy calculation.
3255 C Find the largest exponent
3258 if (emin.gt.contr(j)) emin=contr(j)
3262 C Compute the contribution to SC energy and derivatives
3266 expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin)
3267 escloc_i=escloc_i+expfac
3269 dersc(k)=dersc(k)+Ax(k,j)*expfac
3271 if (mixed) dersc12=dersc12+(-Ax(2,j)*Ax(1,j)
3272 & +gaussc(1,2,j,it))*expfac
3276 dersc(1)=dersc(1)/cos(theti)**2
3277 dersc12=dersc12/cos(theti)**2
3278 escloci=-(dlog(escloc_i)-emin)
3280 dersc(j)=dersc(j)/escloc_i
3282 if (mixed) dersc12=(dersc12/escloc_i+dersc(2)*dersc(1))
3285 c------------------------------------------------------------------------------
3286 subroutine gcont(rij,r0ij,eps0ij,delta,fcont,fprimcont)
3288 C This procedure calculates two-body contact function g(rij) and its derivative:
3291 C g(rij) = esp0ij*(-0.9375*x+0.625*x**3-0.1875*x**5) ! -1 =< x =< 1
3294 C where x=(rij-r0ij)/delta
3296 C rij - interbody distance, r0ij - contact distance, eps0ij - contact energy
3299 double precision rij,r0ij,eps0ij,fcont,fprimcont
3300 double precision x,x2,x4,delta
3304 if (x.lt.-1.0D0) then
3307 else if (x.le.1.0D0) then
3310 fcont=eps0ij*(x*(-0.9375D0+0.6250D0*x2-0.1875D0*x4)+0.5D0)
3311 fprimcont=eps0ij * (-0.9375D0+1.8750D0*x2-0.9375D0*x4)/delta
3318 c------------------------------------------------------------------------------
3319 subroutine splinthet(theti,delta,ss,ssder)
3320 implicit real*8 (a-h,o-z)
3321 include 'DIMENSIONS'
3322 include 'COMMON.VAR'
3323 include 'COMMON.GEO'
3326 if (theti.gt.pipol) then
3327 call gcont(theti,thetup,1.0d0,delta,ss,ssder)
3329 call gcont(-theti,-thetlow,1.0d0,delta,ss,ssder)
3334 c------------------------------------------------------------------------------
3335 subroutine spline1(x,x0,delta,f0,f1,fprim0,f,fprim)
3337 double precision x,x0,delta,f0,f1,fprim0,f,fprim
3338 double precision ksi,ksi2,ksi3,a1,a2,a3
3339 a1=fprim0*delta/(f1-f0)
3345 f=f0+(f1-f0)*ksi*(a1+ksi*(a2+a3*ksi))
3346 fprim=(f1-f0)/delta*(a1+ksi*(2*a2+3*ksi*a3))
3349 c------------------------------------------------------------------------------
3350 subroutine spline2(x,x0,delta,f0x,f1x,fprim0x,fx)
3352 double precision x,x0,delta,f0x,f1x,fprim0x,fx
3353 double precision ksi,ksi2,ksi3,a1,a2,a3
3358 a2=3*(f1x-f0x)-2*fprim0x*delta
3359 a3=fprim0x*delta-2*(f1x-f0x)
3360 fx=f0x+a1*ksi+a2*ksi2+a3*ksi3
3363 C-----------------------------------------------------------------------------
3365 C-----------------------------------------------------------------------------
3366 subroutine etor(etors,edihcnstr)
3367 implicit real*8 (a-h,o-z)
3368 include 'DIMENSIONS'
3369 include 'DIMENSIONS.ZSCOPT'
3370 include 'COMMON.VAR'
3371 include 'COMMON.GEO'
3372 include 'COMMON.LOCAL'
3373 include 'COMMON.TORSION'
3374 include 'COMMON.INTERACT'
3375 include 'COMMON.DERIV'
3376 include 'COMMON.CHAIN'
3377 include 'COMMON.NAMES'
3378 include 'COMMON.IOUNITS'
3379 include 'COMMON.FFIELD'
3380 include 'COMMON.TORCNSTR'
3382 C Set lprn=.true. for debugging
3386 do i=iphi_start,iphi_end
3387 itori=itortyp(itype(i-2))
3388 itori1=itortyp(itype(i-1))
3391 C Proline-Proline pair is a special case...
3392 if (itori.eq.3 .and. itori1.eq.3) then
3393 if (phii.gt.-dwapi3) then
3395 fac=1.0D0/(1.0D0-cosphi)
3396 etorsi=v1(1,3,3)*fac
3397 etorsi=etorsi+etorsi
3398 etors=etors+etorsi-v1(1,3,3)
3399 gloci=gloci-3*fac*etorsi*dsin(3*phii)
3402 v1ij=v1(j+1,itori,itori1)
3403 v2ij=v2(j+1,itori,itori1)
3406 etors=etors+v1ij*cosphi+v2ij*sinphi+dabs(v1ij)+dabs(v2ij)
3407 gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
3411 v1ij=v1(j,itori,itori1)
3412 v2ij=v2(j,itori,itori1)
3415 etors=etors+v1ij*cosphi+v2ij*sinphi+dabs(v1ij)+dabs(v2ij)
3416 gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
3420 & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
3421 & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
3422 & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
3423 gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
3424 c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
3426 ! 6/20/98 - dihedral angle constraints
3429 itori=idih_constr(i)
3432 if (difi.gt.drange(i)) then
3434 edihcnstr=edihcnstr+0.25d0*ftors*difi**4
3435 gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
3436 else if (difi.lt.-drange(i)) then
3438 edihcnstr=edihcnstr+0.25d0*ftors*difi**4
3439 gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
3441 ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
3442 ! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
3444 ! write (iout,*) 'edihcnstr',edihcnstr
3447 c------------------------------------------------------------------------------
3449 subroutine etor(etors,edihcnstr)
3450 implicit real*8 (a-h,o-z)
3451 include 'DIMENSIONS'
3452 include 'DIMENSIONS.ZSCOPT'
3453 include 'COMMON.VAR'
3454 include 'COMMON.GEO'
3455 include 'COMMON.LOCAL'
3456 include 'COMMON.TORSION'
3457 include 'COMMON.INTERACT'
3458 include 'COMMON.DERIV'
3459 include 'COMMON.CHAIN'
3460 include 'COMMON.NAMES'
3461 include 'COMMON.IOUNITS'
3462 include 'COMMON.FFIELD'
3463 include 'COMMON.TORCNSTR'
3465 C Set lprn=.true. for debugging
3469 do i=iphi_start,iphi_end
3470 if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215
3471 itori=itortyp(itype(i-2))
3472 itori1=itortyp(itype(i-1))
3475 C Regular cosine and sine terms
3476 do j=1,nterm(itori,itori1)
3477 v1ij=v1(j,itori,itori1)
3478 v2ij=v2(j,itori,itori1)
3481 etors=etors+v1ij*cosphi+v2ij*sinphi
3482 gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
3486 C E = SUM ----------------------------------- - v1
3487 C [v2 cos(phi/2)+v3 sin(phi/2)]^2 + 1
3489 cosphi=dcos(0.5d0*phii)
3490 sinphi=dsin(0.5d0*phii)
3491 do j=1,nlor(itori,itori1)
3492 vl1ij=vlor1(j,itori,itori1)
3493 vl2ij=vlor2(j,itori,itori1)
3494 vl3ij=vlor3(j,itori,itori1)
3495 pom=vl2ij*cosphi+vl3ij*sinphi
3496 pom1=1.0d0/(pom*pom+1.0d0)
3497 etors=etors+vl1ij*pom1
3499 gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
3501 C Subtract the constant term
3502 etors=etors-v0(itori,itori1)
3504 & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
3505 & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
3506 & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
3507 gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
3508 c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
3511 ! 6/20/98 - dihedral angle constraints
3515 itori=idih_constr(i)
3518 if (difi.gt.drange(i)) then
3520 edihcnstr=edihcnstr+0.25d0*ftors*difi**4
3521 gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
3522 else if (difi.lt.-drange(i)) then
3524 edihcnstr=edihcnstr+0.25d0*ftors*difi**4
3525 gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
3527 ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
3528 ! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
3530 ! write (iout,*) 'edihcnstr',edihcnstr
3533 c----------------------------------------------------------------------------
3534 subroutine etor_d(etors_d)
3535 C 6/23/01 Compute double torsional energy
3536 implicit real*8 (a-h,o-z)
3537 include 'DIMENSIONS'
3538 include 'DIMENSIONS.ZSCOPT'
3539 include 'COMMON.VAR'
3540 include 'COMMON.GEO'
3541 include 'COMMON.LOCAL'
3542 include 'COMMON.TORSION'
3543 include 'COMMON.INTERACT'
3544 include 'COMMON.DERIV'
3545 include 'COMMON.CHAIN'
3546 include 'COMMON.NAMES'
3547 include 'COMMON.IOUNITS'
3548 include 'COMMON.FFIELD'
3549 include 'COMMON.TORCNSTR'
3551 C Set lprn=.true. for debugging
3555 do i=iphi_start,iphi_end-1
3556 if (itel(i-2).eq.0 .or. itel(i-1).eq.0 .or. itel(i).eq.0)
3558 itori=itortyp(itype(i-2))
3559 itori1=itortyp(itype(i-1))
3560 itori2=itortyp(itype(i))
3565 C Regular cosine and sine terms
3566 do j=1,ntermd_1(itori,itori1,itori2)
3567 v1cij=v1c(1,j,itori,itori1,itori2)
3568 v1sij=v1s(1,j,itori,itori1,itori2)
3569 v2cij=v1c(2,j,itori,itori1,itori2)
3570 v2sij=v1s(2,j,itori,itori1,itori2)
3571 cosphi1=dcos(j*phii)
3572 sinphi1=dsin(j*phii)
3573 cosphi2=dcos(j*phii1)
3574 sinphi2=dsin(j*phii1)
3575 etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+
3576 & v2cij*cosphi2+v2sij*sinphi2
3577 gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
3578 gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
3580 do k=2,ntermd_2(itori,itori1,itori2)
3582 v1cdij = v2c(k,l,itori,itori1,itori2)
3583 v2cdij = v2c(l,k,itori,itori1,itori2)
3584 v1sdij = v2s(k,l,itori,itori1,itori2)
3585 v2sdij = v2s(l,k,itori,itori1,itori2)
3586 cosphi1p2=dcos(l*phii+(k-l)*phii1)
3587 cosphi1m2=dcos(l*phii-(k-l)*phii1)
3588 sinphi1p2=dsin(l*phii+(k-l)*phii1)
3589 sinphi1m2=dsin(l*phii-(k-l)*phii1)
3590 etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+
3591 & v1sdij*sinphi1p2+v2sdij*sinphi1m2
3592 gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2
3593 & -v1cdij*sinphi1p2-v2cdij*sinphi1m2)
3594 gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2
3595 & -v1cdij*sinphi1p2+v2cdij*sinphi1m2)
3598 gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
3599 gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
3605 c------------------------------------------------------------------------------
3606 subroutine multibody(ecorr)
3607 C This subroutine calculates multi-body contributions to energy following
3608 C the idea of Skolnick et al. If side chains I and J make a contact and
3609 C at the same time side chains I+1 and J+1 make a contact, an extra
3610 C contribution equal to sqrt(eps(i,j)*eps(i+1,j+1)) is added.
3611 implicit real*8 (a-h,o-z)
3612 include 'DIMENSIONS'
3613 include 'COMMON.IOUNITS'
3614 include 'COMMON.DERIV'
3615 include 'COMMON.INTERACT'
3616 include 'COMMON.CONTACTS'
3617 double precision gx(3),gx1(3)
3620 C Set lprn=.true. for debugging
3624 write (iout,'(a)') 'Contact function values:'
3626 write (iout,'(i2,20(1x,i2,f10.5))')
3627 & i,(jcont(j,i),facont(j,i),j=1,num_cont(i))
3642 num_conti=num_cont(i)
3643 num_conti1=num_cont(i1)
3648 if (j1.eq.j+ishift .or. j1.eq.j-ishift) then
3649 cd write(iout,*)'i=',i,' j=',j,' i1=',i1,' j1=',j1,
3650 cd & ' ishift=',ishift
3651 C Contacts I--J and I+ISHIFT--J+-ISHIFT1 occur simultaneously.
3652 C The system gains extra energy.
3653 ecorr=ecorr+esccorr(i,j,i1,j1,jj,kk)
3654 endif ! j1==j+-ishift
3663 c------------------------------------------------------------------------------
3664 double precision function esccorr(i,j,k,l,jj,kk)
3665 implicit real*8 (a-h,o-z)
3666 include 'DIMENSIONS'
3667 include 'COMMON.IOUNITS'
3668 include 'COMMON.DERIV'
3669 include 'COMMON.INTERACT'
3670 include 'COMMON.CONTACTS'
3671 double precision gx(3),gx1(3)
3676 cd write (iout,'(4i5,3f10.5)') i,j,k,l,eij,ekl,-eij*ekl
3677 C Calculate the multi-body contribution to energy.
3678 C Calculate multi-body contributions to the gradient.
3679 cd write (iout,'(2(2i3,3f10.5))')i,j,(gacont(m,jj,i),m=1,3),
3680 cd & k,l,(gacont(m,kk,k),m=1,3)
3682 gx(m) =ekl*gacont(m,jj,i)
3683 gx1(m)=eij*gacont(m,kk,k)
3684 gradxorr(m,i)=gradxorr(m,i)-gx(m)
3685 gradxorr(m,j)=gradxorr(m,j)+gx(m)
3686 gradxorr(m,k)=gradxorr(m,k)-gx1(m)
3687 gradxorr(m,l)=gradxorr(m,l)+gx1(m)
3691 gradcorr(ll,m)=gradcorr(ll,m)+gx(ll)
3696 gradcorr(ll,m)=gradcorr(ll,m)+gx1(ll)
3702 c------------------------------------------------------------------------------
3704 subroutine pack_buffer(dimen1,dimen2,atom,indx,buffer)
3705 implicit real*8 (a-h,o-z)
3706 include 'DIMENSIONS'
3707 integer dimen1,dimen2,atom,indx
3708 double precision buffer(dimen1,dimen2)
3709 double precision zapas
3710 common /contacts_hb/ zapas(3,20,maxres,7),
3711 & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres),
3712 & num_cont_hb(maxres),jcont_hb(20,maxres)
3713 num_kont=num_cont_hb(atom)
3717 buffer(i,indx+(k-1)*3+j)=zapas(j,i,atom,k)
3720 buffer(i,indx+22)=facont_hb(i,atom)
3721 buffer(i,indx+23)=ees0p(i,atom)
3722 buffer(i,indx+24)=ees0m(i,atom)
3723 buffer(i,indx+25)=dfloat(jcont_hb(i,atom))
3725 buffer(1,indx+26)=dfloat(num_kont)
3728 c------------------------------------------------------------------------------
3729 subroutine unpack_buffer(dimen1,dimen2,atom,indx,buffer)
3730 implicit real*8 (a-h,o-z)
3731 include 'DIMENSIONS'
3732 integer dimen1,dimen2,atom,indx
3733 double precision buffer(dimen1,dimen2)
3734 double precision zapas
3735 common /contacts_hb/ zapas(3,20,maxres,7),
3736 & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres),
3737 & num_cont_hb(maxres),jcont_hb(20,maxres)
3738 num_kont=buffer(1,indx+26)
3739 num_kont_old=num_cont_hb(atom)
3740 num_cont_hb(atom)=num_kont+num_kont_old
3745 zapas(j,ii,atom,k)=buffer(i,indx+(k-1)*3+j)
3748 facont_hb(ii,atom)=buffer(i,indx+22)
3749 ees0p(ii,atom)=buffer(i,indx+23)
3750 ees0m(ii,atom)=buffer(i,indx+24)
3751 jcont_hb(ii,atom)=buffer(i,indx+25)
3755 c------------------------------------------------------------------------------
3757 subroutine multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
3758 C This subroutine calculates multi-body contributions to hydrogen-bonding
3759 implicit real*8 (a-h,o-z)
3760 include 'DIMENSIONS'
3761 include 'DIMENSIONS.ZSCOPT'
3762 include 'COMMON.IOUNITS'
3764 include 'COMMON.INFO'
3766 include 'COMMON.FFIELD'
3767 include 'COMMON.DERIV'
3768 include 'COMMON.INTERACT'
3769 include 'COMMON.CONTACTS'
3771 parameter (max_cont=maxconts)
3772 parameter (max_dim=2*(8*3+2))
3773 parameter (msglen1=max_cont*max_dim*4)
3774 parameter (msglen2=2*msglen1)
3775 integer source,CorrelType,CorrelID,Error
3776 double precision buffer(max_cont,max_dim)
3778 double precision gx(3),gx1(3)
3781 C Set lprn=.true. for debugging
3786 if (fgProcs.le.1) goto 30
3788 write (iout,'(a)') 'Contact function values:'
3790 write (iout,'(2i3,50(1x,i2,f5.2))')
3791 & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
3792 & j=1,num_cont_hb(i))
3795 C Caution! Following code assumes that electrostatic interactions concerning
3796 C a given atom are split among at most two processors!
3806 cd write (iout,*) 'MyRank',MyRank,' mm',mm
3809 cd write (iout,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone
3810 if (MyRank.gt.0) then
3811 C Send correlation contributions to the preceding processor
3813 nn=num_cont_hb(iatel_s)
3814 call pack_buffer(max_cont,max_dim,iatel_s,0,buffer)
3815 cd write (iout,*) 'The BUFFER array:'
3817 cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,26)
3819 if (ielstart(iatel_s).gt.iatel_s+ispp) then
3821 call pack_buffer(max_cont,max_dim,iatel_s+1,26,buffer)
3822 C Clear the contacts of the atom passed to the neighboring processor
3823 nn=num_cont_hb(iatel_s+1)
3825 cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j+26),j=1,26)
3827 num_cont_hb(iatel_s)=0
3829 cd write (iout,*) 'Processor ',MyID,MyRank,
3830 cd & ' is sending correlation contribution to processor',MyID-1,
3831 cd & ' msglen=',msglen
3832 cd write (*,*) 'Processor ',MyID,MyRank,
3833 cd & ' is sending correlation contribution to processor',MyID-1,
3834 cd & ' msglen=',msglen,' CorrelType=',CorrelType
3835 call mp_bsend(buffer,msglen,MyID-1,CorrelType,CorrelID)
3836 cd write (iout,*) 'Processor ',MyID,
3837 cd & ' has sent correlation contribution to processor',MyID-1,
3838 cd & ' msglen=',msglen,' CorrelID=',CorrelID
3839 cd write (*,*) 'Processor ',MyID,
3840 cd & ' has sent correlation contribution to processor',MyID-1,
3841 cd & ' msglen=',msglen,' CorrelID=',CorrelID
3843 endif ! (MyRank.gt.0)
3847 cd write (iout,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone
3848 if (MyRank.lt.fgProcs-1) then
3849 C Receive correlation contributions from the next processor
3851 if (ielend(iatel_e).lt.nct-1) msglen=msglen2
3852 cd write (iout,*) 'Processor',MyID,
3853 cd & ' is receiving correlation contribution from processor',MyID+1,
3854 cd & ' msglen=',msglen,' CorrelType=',CorrelType
3855 cd write (*,*) 'Processor',MyID,
3856 cd & ' is receiving correlation contribution from processor',MyID+1,
3857 cd & ' msglen=',msglen,' CorrelType=',CorrelType
3859 do while (nbytes.le.0)
3860 call mp_probe(MyID+1,CorrelType,nbytes)
3862 cd print *,'Processor',MyID,' msglen',msglen,' nbytes',nbytes
3863 call mp_brecv(buffer,msglen,MyID+1,CorrelType,nbytes)
3864 cd write (iout,*) 'Processor',MyID,
3865 cd & ' has received correlation contribution from processor',MyID+1,
3866 cd & ' msglen=',msglen,' nbytes=',nbytes
3867 cd write (iout,*) 'The received BUFFER array:'
3869 cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,52)
3871 if (msglen.eq.msglen1) then
3872 call unpack_buffer(max_cont,max_dim,iatel_e+1,0,buffer)
3873 else if (msglen.eq.msglen2) then
3874 call unpack_buffer(max_cont,max_dim,iatel_e,0,buffer)
3875 call unpack_buffer(max_cont,max_dim,iatel_e+1,26,buffer)
3878 & 'ERROR!!!! message length changed while processing correlations.'
3880 & 'ERROR!!!! message length changed while processing correlations.'
3881 call mp_stopall(Error)
3882 endif ! msglen.eq.msglen1
3883 endif ! MyRank.lt.fgProcs-1
3890 write (iout,'(a)') 'Contact function values:'
3892 write (iout,'(2i3,50(1x,i2,f5.2))')
3893 & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
3894 & j=1,num_cont_hb(i))
3898 C Remove the loop below after debugging !!!
3905 C Calculate the local-electrostatic correlation terms
3906 do i=iatel_s,iatel_e+1
3908 num_conti=num_cont_hb(i)
3909 num_conti1=num_cont_hb(i+1)
3914 c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
3915 c & ' jj=',jj,' kk=',kk
3916 if (j1.eq.j+1 .or. j1.eq.j-1) then
3917 C Contacts I-J and (I+1)-(J+1) or (I+1)-(J-1) occur simultaneously.
3918 C The system gains extra energy.
3919 ecorr=ecorr+ehbcorr(i,j,i+1,j1,jj,kk,0.72D0,0.32D0)
3921 else if (j1.eq.j) then
3922 C Contacts I-J and I-(J+1) occur simultaneously.
3923 C The system loses extra energy.
3924 c ecorr=ecorr+ehbcorr(i,j,i+1,j,jj,kk,0.60D0,-0.40D0)
3929 c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
3930 c & ' jj=',jj,' kk=',kk
3932 C Contacts I-J and (I+1)-J occur simultaneously.
3933 C The system loses extra energy.
3934 c ecorr=ecorr+ehbcorr(i,j,i,j+1,jj,kk,0.60D0,-0.40D0)
3941 c------------------------------------------------------------------------------
3942 subroutine multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,
3944 C This subroutine calculates multi-body contributions to hydrogen-bonding
3945 implicit real*8 (a-h,o-z)
3946 include 'DIMENSIONS'
3947 include 'DIMENSIONS.ZSCOPT'
3948 include 'COMMON.IOUNITS'
3950 include 'COMMON.INFO'
3952 include 'COMMON.FFIELD'
3953 include 'COMMON.DERIV'
3954 include 'COMMON.INTERACT'
3955 include 'COMMON.CONTACTS'
3957 parameter (max_cont=maxconts)
3958 parameter (max_dim=2*(8*3+2))
3959 parameter (msglen1=max_cont*max_dim*4)
3960 parameter (msglen2=2*msglen1)
3961 integer source,CorrelType,CorrelID,Error
3962 double precision buffer(max_cont,max_dim)
3964 double precision gx(3),gx1(3)
3967 C Set lprn=.true. for debugging
3973 if (fgProcs.le.1) goto 30
3975 write (iout,'(a)') 'Contact function values:'
3977 write (iout,'(2i3,50(1x,i2,f5.2))')
3978 & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
3979 & j=1,num_cont_hb(i))
3982 C Caution! Following code assumes that electrostatic interactions concerning
3983 C a given atom are split among at most two processors!
3993 cd write (iout,*) 'MyRank',MyRank,' mm',mm
3996 cd write (iout,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone
3997 if (MyRank.gt.0) then
3998 C Send correlation contributions to the preceding processor
4000 nn=num_cont_hb(iatel_s)
4001 call pack_buffer(max_cont,max_dim,iatel_s,0,buffer)
4002 cd write (iout,*) 'The BUFFER array:'
4004 cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,26)
4006 if (ielstart(iatel_s).gt.iatel_s+ispp) then
4008 call pack_buffer(max_cont,max_dim,iatel_s+1,26,buffer)
4009 C Clear the contacts of the atom passed to the neighboring processor
4010 nn=num_cont_hb(iatel_s+1)
4012 cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j+26),j=1,26)
4014 num_cont_hb(iatel_s)=0
4016 cd write (iout,*) 'Processor ',MyID,MyRank,
4017 cd & ' is sending correlation contribution to processor',MyID-1,
4018 cd & ' msglen=',msglen
4019 cd write (*,*) 'Processor ',MyID,MyRank,
4020 cd & ' is sending correlation contribution to processor',MyID-1,
4021 cd & ' msglen=',msglen,' CorrelType=',CorrelType
4022 call mp_bsend(buffer,msglen,MyID-1,CorrelType,CorrelID)
4023 cd write (iout,*) 'Processor ',MyID,
4024 cd & ' has sent correlation contribution to processor',MyID-1,
4025 cd & ' msglen=',msglen,' CorrelID=',CorrelID
4026 cd write (*,*) 'Processor ',MyID,
4027 cd & ' has sent correlation contribution to processor',MyID-1,
4028 cd & ' msglen=',msglen,' CorrelID=',CorrelID
4030 endif ! (MyRank.gt.0)
4034 cd write (iout,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone
4035 if (MyRank.lt.fgProcs-1) then
4036 C Receive correlation contributions from the next processor
4038 if (ielend(iatel_e).lt.nct-1) msglen=msglen2
4039 cd write (iout,*) 'Processor',MyID,
4040 cd & ' is receiving correlation contribution from processor',MyID+1,
4041 cd & ' msglen=',msglen,' CorrelType=',CorrelType
4042 cd write (*,*) 'Processor',MyID,
4043 cd & ' is receiving correlation contribution from processor',MyID+1,
4044 cd & ' msglen=',msglen,' CorrelType=',CorrelType
4046 do while (nbytes.le.0)
4047 call mp_probe(MyID+1,CorrelType,nbytes)
4049 cd print *,'Processor',MyID,' msglen',msglen,' nbytes',nbytes
4050 call mp_brecv(buffer,msglen,MyID+1,CorrelType,nbytes)
4051 cd write (iout,*) 'Processor',MyID,
4052 cd & ' has received correlation contribution from processor',MyID+1,
4053 cd & ' msglen=',msglen,' nbytes=',nbytes
4054 cd write (iout,*) 'The received BUFFER array:'
4056 cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,52)
4058 if (msglen.eq.msglen1) then
4059 call unpack_buffer(max_cont,max_dim,iatel_e+1,0,buffer)
4060 else if (msglen.eq.msglen2) then
4061 call unpack_buffer(max_cont,max_dim,iatel_e,0,buffer)
4062 call unpack_buffer(max_cont,max_dim,iatel_e+1,26,buffer)
4065 & 'ERROR!!!! message length changed while processing correlations.'
4067 & 'ERROR!!!! message length changed while processing correlations.'
4068 call mp_stopall(Error)
4069 endif ! msglen.eq.msglen1
4070 endif ! MyRank.lt.fgProcs-1
4077 write (iout,'(a)') 'Contact function values:'
4079 write (iout,'(2i3,50(1x,i2,f5.2))')
4080 & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
4081 & j=1,num_cont_hb(i))
4087 C Remove the loop below after debugging !!!
4094 C Calculate the dipole-dipole interaction energies
4095 if (wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) then
4096 do i=iatel_s,iatel_e+1
4097 num_conti=num_cont_hb(i)
4104 C Calculate the local-electrostatic correlation terms
4105 do i=iatel_s,iatel_e+1
4107 num_conti=num_cont_hb(i)
4108 num_conti1=num_cont_hb(i+1)
4113 c write (*,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
4114 c & ' jj=',jj,' kk=',kk
4115 if (j1.eq.j+1 .or. j1.eq.j-1) then
4116 C Contacts I-J and (I+1)-(J+1) or (I+1)-(J-1) occur simultaneously.
4117 C The system gains extra energy.
4119 sqd1=dsqrt(d_cont(jj,i))
4120 sqd2=dsqrt(d_cont(kk,i1))
4121 sred_geom = sqd1*sqd2
4122 IF (sred_geom.lt.cutoff_corr) THEN
4123 call gcont(sred_geom,r0_corr,1.0D0,delt_corr,
4125 c write (*,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
4126 c & ' jj=',jj,' kk=',kk
4127 fac_prim1=0.5d0*sqd2/sqd1*fprimcont
4128 fac_prim2=0.5d0*sqd1/sqd2*fprimcont
4130 g_contij(l,1)=fac_prim1*grij_hb_cont(l,jj,i)
4131 g_contij(l,2)=fac_prim2*grij_hb_cont(l,kk,i1)
4134 cd write (iout,*) 'sred_geom=',sred_geom,
4135 cd & ' ekont=',ekont,' fprim=',fprimcont
4136 call calc_eello(i,j,i+1,j1,jj,kk)
4137 if (wcorr4.gt.0.0d0)
4138 & ecorr=ecorr+eello4(i,j,i+1,j1,jj,kk)
4139 if (wcorr5.gt.0.0d0)
4140 & ecorr5=ecorr5+eello5(i,j,i+1,j1,jj,kk)
4141 c print *,"wcorr5",ecorr5
4142 cd write(2,*)'wcorr6',wcorr6,' wturn6',wturn6
4143 cd write(2,*)'ijkl',i,j,i+1,j1
4144 if (wcorr6.gt.0.0d0 .and. (j.ne.i+4 .or. j1.ne.i+3
4145 & .or. wturn6.eq.0.0d0))then
4146 cd write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
4147 ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk)
4148 cd write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
4149 cd & 'ecorr6=',ecorr6
4150 cd write (iout,'(4e15.5)') sred_geom,
4151 cd & dabs(eello4(i,j,i+1,j1,jj,kk)),
4152 cd & dabs(eello5(i,j,i+1,j1,jj,kk)),
4153 cd & dabs(eello6(i,j,i+1,j1,jj,kk))
4154 else if (wturn6.gt.0.0d0
4155 & .and. (j.eq.i+4 .and. j1.eq.i+3)) then
4156 cd write (iout,*) '******eturn6: i,j,i+1,j1',i,j,i+1,j1
4157 eturn6=eturn6+eello_turn6(i,jj,kk)
4158 cd write (2,*) 'multibody_eello:eturn6',eturn6
4162 else if (j1.eq.j) then
4163 C Contacts I-J and I-(J+1) occur simultaneously.
4164 C The system loses extra energy.
4165 c ecorr=ecorr+ehbcorr(i,j,i+1,j,jj,kk,0.60D0,-0.40D0)
4170 c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
4171 c & ' jj=',jj,' kk=',kk
4173 C Contacts I-J and (I+1)-J occur simultaneously.
4174 C The system loses extra energy.
4175 c ecorr=ecorr+ehbcorr(i,j,i,j+1,jj,kk,0.60D0,-0.40D0)
4182 c------------------------------------------------------------------------------
4183 double precision function ehbcorr(i,j,k,l,jj,kk,coeffp,coeffm)
4184 implicit real*8 (a-h,o-z)
4185 include 'DIMENSIONS'
4186 include 'COMMON.IOUNITS'
4187 include 'COMMON.DERIV'
4188 include 'COMMON.INTERACT'
4189 include 'COMMON.CONTACTS'
4190 double precision gx(3),gx1(3)
4200 ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
4201 cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
4202 C Following 4 lines for diagnostics.
4207 c write (iout,*)'Contacts have occurred for peptide groups',i,j,
4209 c write (iout,*)'Contacts have occurred for peptide groups',
4210 c & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
4211 c & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
4212 C Calculate the multi-body contribution to energy.
4213 ecorr=ecorr+ekont*ees
4215 C Calculate multi-body contributions to the gradient.
4217 ghalf=0.5D0*ees*ekl*gacont_hbr(ll,jj,i)
4218 gradcorr(ll,i)=gradcorr(ll,i)+ghalf
4219 & -ekont*(coeffp*ees0pkl*gacontp_hb1(ll,jj,i)+
4220 & coeffm*ees0mkl*gacontm_hb1(ll,jj,i))
4221 gradcorr(ll,j)=gradcorr(ll,j)+ghalf
4222 & -ekont*(coeffp*ees0pkl*gacontp_hb2(ll,jj,i)+
4223 & coeffm*ees0mkl*gacontm_hb2(ll,jj,i))
4224 ghalf=0.5D0*ees*eij*gacont_hbr(ll,kk,k)
4225 gradcorr(ll,k)=gradcorr(ll,k)+ghalf
4226 & -ekont*(coeffp*ees0pij*gacontp_hb1(ll,kk,k)+
4227 & coeffm*ees0mij*gacontm_hb1(ll,kk,k))
4228 gradcorr(ll,l)=gradcorr(ll,l)+ghalf
4229 & -ekont*(coeffp*ees0pij*gacontp_hb2(ll,kk,k)+
4230 & coeffm*ees0mij*gacontm_hb2(ll,kk,k))
4234 gradcorr(ll,m)=gradcorr(ll,m)+
4235 & ees*ekl*gacont_hbr(ll,jj,i)-
4236 & ekont*(coeffp*ees0pkl*gacontp_hb3(ll,jj,i)+
4237 & coeffm*ees0mkl*gacontm_hb3(ll,jj,i))
4242 gradcorr(ll,m)=gradcorr(ll,m)+
4243 & ees*eij*gacont_hbr(ll,kk,k)-
4244 & ekont*(coeffp*ees0pij*gacontp_hb3(ll,kk,k)+
4245 & coeffm*ees0mij*gacontm_hb3(ll,kk,k))
4252 C---------------------------------------------------------------------------
4253 subroutine dipole(i,j,jj)
4254 implicit real*8 (a-h,o-z)
4255 include 'DIMENSIONS'
4256 include 'DIMENSIONS.ZSCOPT'
4257 include 'COMMON.IOUNITS'
4258 include 'COMMON.CHAIN'
4259 include 'COMMON.FFIELD'
4260 include 'COMMON.DERIV'
4261 include 'COMMON.INTERACT'
4262 include 'COMMON.CONTACTS'
4263 include 'COMMON.TORSION'
4264 include 'COMMON.VAR'
4265 include 'COMMON.GEO'
4266 dimension dipi(2,2),dipj(2,2),dipderi(2),dipderj(2),auxvec(2),
4268 iti1 = itortyp(itype(i+1))
4269 if (j.lt.nres-1) then
4270 itj1 = itortyp(itype(j+1))
4275 dipi(iii,1)=Ub2(iii,i)
4276 dipderi(iii)=Ub2der(iii,i)
4277 dipi(iii,2)=b1(iii,iti1)
4278 dipj(iii,1)=Ub2(iii,j)
4279 dipderj(iii)=Ub2der(iii,j)
4280 dipj(iii,2)=b1(iii,itj1)
4284 call matvec2(a_chuj(1,1,jj,i),dipj(1,iii),auxvec(1))
4287 dip(kkk,jj,i)=scalar2(dipi(1,jjj),auxvec(1))
4290 if (.not.calc_grad) return
4295 call matvec2(a_chuj_der(1,1,lll,kkk,jj,i),dipj(1,iii),
4299 dipderx(lll,kkk,mmm,jj,i)=scalar2(dipi(1,jjj),auxvec(1))
4304 call transpose2(a_chuj(1,1,jj,i),auxmat(1,1))
4305 call matvec2(auxmat(1,1),dipderi(1),auxvec(1))
4307 dipderg(iii,jj,i)=scalar2(auxvec(1),dipj(1,iii))
4309 call matvec2(a_chuj(1,1,jj,i),dipderj(1),auxvec(1))
4311 dipderg(iii+2,jj,i)=scalar2(auxvec(1),dipi(1,iii))
4315 C---------------------------------------------------------------------------
4316 subroutine calc_eello(i,j,k,l,jj,kk)
4318 C This subroutine computes matrices and vectors needed to calculate
4319 C the fourth-, fifth-, and sixth-order local-electrostatic terms.
4321 implicit real*8 (a-h,o-z)
4322 include 'DIMENSIONS'
4323 include 'DIMENSIONS.ZSCOPT'
4324 include 'COMMON.IOUNITS'
4325 include 'COMMON.CHAIN'
4326 include 'COMMON.DERIV'
4327 include 'COMMON.INTERACT'
4328 include 'COMMON.CONTACTS'
4329 include 'COMMON.TORSION'
4330 include 'COMMON.VAR'
4331 include 'COMMON.GEO'
4332 include 'COMMON.FFIELD'
4333 double precision aa1(2,2),aa2(2,2),aa1t(2,2),aa2t(2,2),
4334 & aa1tder(2,2,3,5),aa2tder(2,2,3,5),auxmat(2,2)
4337 cd write (iout,*) 'calc_eello: i=',i,' j=',j,' k=',k,' l=',l,
4338 cd & ' jj=',jj,' kk=',kk
4339 cd if (i.ne.2 .or. j.ne.4 .or. k.ne.3 .or. l.ne.5) return
4342 aa1(iii,jjj)=a_chuj(iii,jjj,jj,i)
4343 aa2(iii,jjj)=a_chuj(iii,jjj,kk,k)
4346 call transpose2(aa1(1,1),aa1t(1,1))
4347 call transpose2(aa2(1,1),aa2t(1,1))
4350 call transpose2(a_chuj_der(1,1,lll,kkk,jj,i),
4351 & aa1tder(1,1,lll,kkk))
4352 call transpose2(a_chuj_der(1,1,lll,kkk,kk,k),
4353 & aa2tder(1,1,lll,kkk))
4357 C parallel orientation of the two CA-CA-CA frames.
4359 iti=itortyp(itype(i))
4363 itk1=itortyp(itype(k+1))
4364 itj=itortyp(itype(j))
4365 if (l.lt.nres-1) then
4366 itl1=itortyp(itype(l+1))
4370 C A1 kernel(j+1) A2T
4372 cd write (iout,'(3f10.5,5x,3f10.5)')
4373 cd & (EUg(iii,jjj,k),jjj=1,2),(EUg(iii,jjj,l),jjj=1,2)
4375 call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i),
4376 & aa2tder(1,1,1,1),1,.false.,EUg(1,1,l),EUgder(1,1,l),
4377 & AEA(1,1,1),AEAderg(1,1,1),AEAderx(1,1,1,1,1,1))
4378 C Following matrices are needed only for 6-th order cumulants
4379 IF (wcorr6.gt.0.0d0) THEN
4380 call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i),
4381 & aa2tder(1,1,1,1),1,.false.,EUgC(1,1,l),EUgCder(1,1,l),
4382 & AECA(1,1,1),AECAderg(1,1,1),AECAderx(1,1,1,1,1,1))
4383 call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i),
4384 & aa2tder(1,1,1,1),2,.false.,Ug2DtEUg(1,1,l),
4385 & Ug2DtEUgder(1,1,1,l),ADtEA(1,1,1),ADtEAderg(1,1,1,1),
4386 & ADtEAderx(1,1,1,1,1,1))
4388 call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i),
4389 & aa2tder(1,1,1,1),2,.false.,DtUg2EUg(1,1,l),
4390 & DtUg2EUgder(1,1,1,l),ADtEA1(1,1,1),ADtEA1derg(1,1,1,1),
4391 & ADtEA1derx(1,1,1,1,1,1))
4393 C End 6-th order cumulants
4396 cd write (2,*) 'In calc_eello6'
4398 cd write (2,*) 'iii=',iii
4400 cd write (2,*) 'kkk=',kkk
4402 cd write (2,'(3(2f10.5),5x)')
4403 cd & ((ADtEA1derx(jjj,mmm,lll,kkk,iii,1),mmm=1,2),lll=1,3)
4408 call transpose2(EUgder(1,1,k),auxmat(1,1))
4409 call matmat2(auxmat(1,1),AEA(1,1,1),EAEAderg(1,1,1,1))
4410 call transpose2(EUg(1,1,k),auxmat(1,1))
4411 call matmat2(auxmat(1,1),AEA(1,1,1),EAEA(1,1,1))
4412 call matmat2(auxmat(1,1),AEAderg(1,1,1),EAEAderg(1,1,2,1))
4416 call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,1),
4417 & EAEAderx(1,1,lll,kkk,iii,1))
4421 C A1T kernel(i+1) A2
4422 call kernel(aa1t(1,1),aa2(1,1),aa1tder(1,1,1,1),
4423 & a_chuj_der(1,1,1,1,kk,k),1,.false.,EUg(1,1,k),EUgder(1,1,k),
4424 & AEA(1,1,2),AEAderg(1,1,2),AEAderx(1,1,1,1,1,2))
4425 C Following matrices are needed only for 6-th order cumulants
4426 IF (wcorr6.gt.0.0d0) THEN
4427 call kernel(aa1t(1,1),aa2(1,1),aa1tder(1,1,1,1),
4428 & a_chuj_der(1,1,1,1,kk,k),1,.false.,EUgC(1,1,k),EUgCder(1,1,k),
4429 & AECA(1,1,2),AECAderg(1,1,2),AECAderx(1,1,1,1,1,2))
4430 call kernel(aa1t(1,1),aa2(1,1),aa1tder(1,1,1,1),
4431 & a_chuj_der(1,1,1,1,kk,k),2,.false.,Ug2DtEUg(1,1,k),
4432 & Ug2DtEUgder(1,1,1,k),ADtEA(1,1,2),ADtEAderg(1,1,1,2),
4433 & ADtEAderx(1,1,1,1,1,2))
4434 call kernel(aa1t(1,1),aa2(1,1),aa1tder(1,1,1,1),
4435 & a_chuj_der(1,1,1,1,kk,k),2,.false.,DtUg2EUg(1,1,k),
4436 & DtUg2EUgder(1,1,1,k),ADtEA1(1,1,2),ADtEA1derg(1,1,1,2),
4437 & ADtEA1derx(1,1,1,1,1,2))
4439 C End 6-th order cumulants
4440 call transpose2(EUgder(1,1,l),auxmat(1,1))
4441 call matmat2(auxmat(1,1),AEA(1,1,2),EAEAderg(1,1,1,2))
4442 call transpose2(EUg(1,1,l),auxmat(1,1))
4443 call matmat2(auxmat(1,1),AEA(1,1,2),EAEA(1,1,2))
4444 call matmat2(auxmat(1,1),AEAderg(1,1,2),EAEAderg(1,1,2,2))
4448 call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,2),
4449 & EAEAderx(1,1,lll,kkk,iii,2))
4454 C Calculate the vectors and their derivatives in virtual-bond dihedral angles.
4455 C They are needed only when the fifth- or the sixth-order cumulants are
4457 IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) THEN
4458 call transpose2(AEA(1,1,1),auxmat(1,1))
4459 call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1))
4460 call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1))
4461 call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1))
4462 call transpose2(AEAderg(1,1,1),auxmat(1,1))
4463 call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1))
4464 call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1))
4465 call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1))
4466 call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1))
4467 call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1))
4468 call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1))
4469 call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1))
4470 call transpose2(AEA(1,1,2),auxmat(1,1))
4471 call matvec2(auxmat(1,1),b1(1,itj),AEAb1(1,1,2))
4472 call matvec2(auxmat(1,1),Ub2(1,j),AEAb2(1,1,2))
4473 call matvec2(auxmat(1,1),Ub2der(1,j),AEAb2derg(1,2,1,2))
4474 call transpose2(AEAderg(1,1,2),auxmat(1,1))
4475 call matvec2(auxmat(1,1),b1(1,itj),AEAb1derg(1,1,2))
4476 call matvec2(auxmat(1,1),Ub2(1,j),AEAb2derg(1,1,1,2))
4477 call matvec2(AEA(1,1,2),b1(1,itl1),AEAb1(1,2,2))
4478 call matvec2(AEAderg(1,1,2),b1(1,itl1),AEAb1derg(1,2,2))
4479 call matvec2(AEA(1,1,2),Ub2(1,l+1),AEAb2(1,2,2))
4480 call matvec2(AEAderg(1,1,2),Ub2(1,l+1),AEAb2derg(1,1,2,2))
4481 call matvec2(AEA(1,1,2),Ub2der(1,l+1),AEAb2derg(1,2,2,2))
4482 C Calculate the Cartesian derivatives of the vectors.
4486 call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1))
4487 call matvec2(auxmat(1,1),b1(1,iti),
4488 & AEAb1derx(1,lll,kkk,iii,1,1))
4489 call matvec2(auxmat(1,1),Ub2(1,i),
4490 & AEAb2derx(1,lll,kkk,iii,1,1))
4491 call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
4492 & AEAb1derx(1,lll,kkk,iii,2,1))
4493 call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1),
4494 & AEAb2derx(1,lll,kkk,iii,2,1))
4495 call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1))
4496 call matvec2(auxmat(1,1),b1(1,itj),
4497 & AEAb1derx(1,lll,kkk,iii,1,2))
4498 call matvec2(auxmat(1,1),Ub2(1,j),
4499 & AEAb2derx(1,lll,kkk,iii,1,2))
4500 call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itl1),
4501 & AEAb1derx(1,lll,kkk,iii,2,2))
4502 call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,l+1),
4503 & AEAb2derx(1,lll,kkk,iii,2,2))
4510 C Antiparallel orientation of the two CA-CA-CA frames.
4512 iti=itortyp(itype(i))
4516 itk1=itortyp(itype(k+1))
4517 itl=itortyp(itype(l))
4518 itj=itortyp(itype(j))
4519 if (j.lt.nres-1) then
4520 itj1=itortyp(itype(j+1))
4524 C A2 kernel(j-1)T A1T
4525 call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i),
4526 & aa2tder(1,1,1,1),1,.true.,EUg(1,1,j),EUgder(1,1,j),
4527 & AEA(1,1,1),AEAderg(1,1,1),AEAderx(1,1,1,1,1,1))
4528 C Following matrices are needed only for 6-th order cumulants
4529 IF (wcorr6.gt.0.0d0 .or. (wturn6.gt.0.0d0 .and.
4530 & j.eq.i+4 .and. l.eq.i+3)) THEN
4531 call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i),
4532 & aa2tder(1,1,1,1),1,.true.,EUgC(1,1,j),EUgCder(1,1,j),
4533 & AECA(1,1,1),AECAderg(1,1,1),AECAderx(1,1,1,1,1,1))
4534 call kernel(aa2(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i),
4535 & aa2tder(1,1,1,1),2,.true.,Ug2DtEUg(1,1,j),
4536 & Ug2DtEUgder(1,1,1,j),ADtEA(1,1,1),ADtEAderg(1,1,1,1),
4537 & ADtEAderx(1,1,1,1,1,1))
4538 call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i),
4539 & aa2tder(1,1,1,1),2,.true.,DtUg2EUg(1,1,j),
4540 & DtUg2EUgder(1,1,1,j),ADtEA1(1,1,1),ADtEA1derg(1,1,1,1),
4541 & ADtEA1derx(1,1,1,1,1,1))
4543 C End 6-th order cumulants
4544 call transpose2(EUgder(1,1,k),auxmat(1,1))
4545 call matmat2(auxmat(1,1),AEA(1,1,1),EAEAderg(1,1,1,1))
4546 call transpose2(EUg(1,1,k),auxmat(1,1))
4547 call matmat2(auxmat(1,1),AEA(1,1,1),EAEA(1,1,1))
4548 call matmat2(auxmat(1,1),AEAderg(1,1,1),EAEAderg(1,1,2,1))
4552 call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,1),
4553 & EAEAderx(1,1,lll,kkk,iii,1))
4557 C A2T kernel(i+1)T A1
4558 call kernel(aa2t(1,1),aa1(1,1),aa2tder(1,1,1,1),
4559 & a_chuj_der(1,1,1,1,jj,i),1,.true.,EUg(1,1,k),EUgder(1,1,k),
4560 & AEA(1,1,2),AEAderg(1,1,2),AEAderx(1,1,1,1,1,2))
4561 C Following matrices are needed only for 6-th order cumulants
4562 IF (wcorr6.gt.0.0d0 .or. (wturn6.gt.0.0d0 .and.
4563 & j.eq.i+4 .and. l.eq.i+3)) THEN
4564 call kernel(aa2t(1,1),aa1(1,1),aa2tder(1,1,1,1),
4565 & a_chuj_der(1,1,1,1,jj,i),1,.true.,EUgC(1,1,k),EUgCder(1,1,k),
4566 & AECA(1,1,2),AECAderg(1,1,2),AECAderx(1,1,1,1,1,2))
4567 call kernel(aa2t(1,1),aa1(1,1),aa2tder(1,1,1,1),
4568 & a_chuj_der(1,1,1,1,jj,i),2,.true.,Ug2DtEUg(1,1,k),
4569 & Ug2DtEUgder(1,1,1,k),ADtEA(1,1,2),ADtEAderg(1,1,1,2),
4570 & ADtEAderx(1,1,1,1,1,2))
4571 call kernel(aa2t(1,1),aa1(1,1),aa2tder(1,1,1,1),
4572 & a_chuj_der(1,1,1,1,jj,i),2,.true.,DtUg2EUg(1,1,k),
4573 & DtUg2EUgder(1,1,1,k),ADtEA1(1,1,2),ADtEA1derg(1,1,1,2),
4574 & ADtEA1derx(1,1,1,1,1,2))
4576 C End 6-th order cumulants
4577 call transpose2(EUgder(1,1,j),auxmat(1,1))
4578 call matmat2(auxmat(1,1),AEA(1,1,1),EAEAderg(1,1,2,2))
4579 call transpose2(EUg(1,1,j),auxmat(1,1))
4580 call matmat2(auxmat(1,1),AEA(1,1,2),EAEA(1,1,2))
4581 call matmat2(auxmat(1,1),AEAderg(1,1,2),EAEAderg(1,1,2,2))
4585 call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,2),
4586 & EAEAderx(1,1,lll,kkk,iii,2))
4591 C Calculate the vectors and their derivatives in virtual-bond dihedral angles.
4592 C They are needed only when the fifth- or the sixth-order cumulants are
4594 IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 .or.
4595 & (wturn6.gt.0.0d0 .and. j.eq.i+4 .and. l.eq.i+3)) THEN
4596 call transpose2(AEA(1,1,1),auxmat(1,1))
4597 call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1))
4598 call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1))
4599 call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1))
4600 call transpose2(AEAderg(1,1,1),auxmat(1,1))
4601 call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1))
4602 call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1))
4603 call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1))
4604 call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1))
4605 call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1))
4606 call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1))
4607 call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1))
4608 call transpose2(AEA(1,1,2),auxmat(1,1))
4609 call matvec2(auxmat(1,1),b1(1,itj1),AEAb1(1,1,2))
4610 call matvec2(auxmat(1,1),Ub2(1,l),AEAb2(1,1,2))
4611 call matvec2(auxmat(1,1),Ub2der(1,l),AEAb2derg(1,2,1,2))
4612 call transpose2(AEAderg(1,1,2),auxmat(1,1))
4613 call matvec2(auxmat(1,1),b1(1,itl),AEAb1(1,1,2))
4614 call matvec2(auxmat(1,1),Ub2(1,l),AEAb2derg(1,1,1,2))
4615 call matvec2(AEA(1,1,2),b1(1,itj1),AEAb1(1,2,2))
4616 call matvec2(AEAderg(1,1,2),b1(1,itj1),AEAb1derg(1,2,2))
4617 call matvec2(AEA(1,1,2),Ub2(1,j),AEAb2(1,2,2))
4618 call matvec2(AEAderg(1,1,2),Ub2(1,j),AEAb2derg(1,1,2,2))
4619 call matvec2(AEA(1,1,2),Ub2der(1,j),AEAb2derg(1,2,2,2))
4620 C Calculate the Cartesian derivatives of the vectors.
4624 call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1))
4625 call matvec2(auxmat(1,1),b1(1,iti),
4626 & AEAb1derx(1,lll,kkk,iii,1,1))
4627 call matvec2(auxmat(1,1),Ub2(1,i),
4628 & AEAb2derx(1,lll,kkk,iii,1,1))
4629 call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
4630 & AEAb1derx(1,lll,kkk,iii,2,1))
4631 call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1),
4632 & AEAb2derx(1,lll,kkk,iii,2,1))
4633 call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1))
4634 call matvec2(auxmat(1,1),b1(1,itl),
4635 & AEAb1derx(1,lll,kkk,iii,1,2))
4636 call matvec2(auxmat(1,1),Ub2(1,l),
4637 & AEAb2derx(1,lll,kkk,iii,1,2))
4638 call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itj1),
4639 & AEAb1derx(1,lll,kkk,iii,2,2))
4640 call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,j),
4641 & AEAb2derx(1,lll,kkk,iii,2,2))
4650 C---------------------------------------------------------------------------
4651 subroutine kernel(aa1,aa2t,aa1derx,aa2tderx,nderg,transp,
4652 & KK,KKderg,AKA,AKAderg,AKAderx)
4656 double precision aa1(2,2),aa2t(2,2),aa1derx(2,2,3,5),
4657 & aa2tderx(2,2,3,5),KK(2,2),KKderg(2,2,nderg),AKA(2,2),
4658 & AKAderg(2,2,nderg),AKAderx(2,2,3,5,2)
4663 call prodmat3(aa1(1,1),aa2t(1,1),KK(1,1),transp,AKA(1,1))
4665 call prodmat3(aa1(1,1),aa2t(1,1),KKderg(1,1,iii),transp,
4668 cd if (lprn) write (2,*) 'In kernel'
4670 cd if (lprn) write (2,*) 'kkk=',kkk
4672 call prodmat3(aa1derx(1,1,lll,kkk),aa2t(1,1),
4673 & KK(1,1),transp,AKAderx(1,1,lll,kkk,1))
4675 cd write (2,*) 'lll=',lll
4676 cd write (2,*) 'iii=1'
4678 cd write (2,'(3(2f10.5),5x)')
4679 cd & (AKAderx(jjj,mmm,lll,kkk,1),mmm=1,2)
4682 call prodmat3(aa1(1,1),aa2tderx(1,1,lll,kkk),
4683 & KK(1,1),transp,AKAderx(1,1,lll,kkk,2))
4685 cd write (2,*) 'lll=',lll
4686 cd write (2,*) 'iii=2'
4688 cd write (2,'(3(2f10.5),5x)')
4689 cd & (AKAderx(jjj,mmm,lll,kkk,2),mmm=1,2)
4696 C---------------------------------------------------------------------------
4697 double precision function eello4(i,j,k,l,jj,kk)
4698 implicit real*8 (a-h,o-z)
4699 include 'DIMENSIONS'
4700 include 'COMMON.IOUNITS'
4701 include 'COMMON.CHAIN'
4702 include 'COMMON.DERIV'
4703 include 'COMMON.INTERACT'
4704 include 'COMMON.CONTACTS'
4705 include 'COMMON.TORSION'
4706 include 'COMMON.VAR'
4707 include 'COMMON.GEO'
4708 double precision pizda(2,2),ggg1(3),ggg2(3)
4709 cd if (i.ne.1 .or. j.ne.5 .or. k.ne.2 .or.l.ne.4) then
4713 cd print *,'eello4:',i,j,k,l,jj,kk
4714 cd write (2,*) 'i',i,' j',j,' k',k,' l',l
4715 cd call checkint4(i,j,k,l,jj,kk,eel4_num)
4716 cold eij=facont_hb(jj,i)
4717 cold ekl=facont_hb(kk,k)
4719 eel4=-EAEA(1,1,1)-EAEA(2,2,1)
4721 cd eel41=-EAEA(1,1,2)-EAEA(2,2,2)
4722 gcorr_loc(k-1)=gcorr_loc(k-1)
4723 & -ekont*(EAEAderg(1,1,1,1)+EAEAderg(2,2,1,1))
4725 gcorr_loc(l-1)=gcorr_loc(l-1)
4726 & -ekont*(EAEAderg(1,1,2,1)+EAEAderg(2,2,2,1))
4728 gcorr_loc(j-1)=gcorr_loc(j-1)
4729 & -ekont*(EAEAderg(1,1,2,1)+EAEAderg(2,2,2,1))
4734 derx(lll,kkk,iii)=-EAEAderx(1,1,lll,kkk,iii,1)
4735 & -EAEAderx(2,2,lll,kkk,iii,1)
4736 cd derx(lll,kkk,iii)=0.0d0
4740 cd gcorr_loc(l-1)=0.0d0
4741 cd gcorr_loc(j-1)=0.0d0
4742 cd gcorr_loc(k-1)=0.0d0
4744 cd write (iout,*)'Contacts have occurred for peptide groups',
4745 cd & i,j,' fcont:',eij,' eij',' and ',k,l,
4746 cd & ' fcont ',ekl,' eel4=',eel4,' eel4_num',16*eel4_num
4747 if (j.lt.nres-1) then
4754 if (l.lt.nres-1) then
4762 cold ghalf=0.5d0*eel4*ekl*gacont_hbr(ll,jj,i)
4763 ggg1(ll)=eel4*g_contij(ll,1)
4764 ggg2(ll)=eel4*g_contij(ll,2)
4765 ghalf=0.5d0*ggg1(ll)
4767 gradcorr(ll,i)=gradcorr(ll,i)+ghalf+ekont*derx(ll,2,1)
4768 gradcorr(ll,i+1)=gradcorr(ll,i+1)+ekont*derx(ll,3,1)
4769 gradcorr(ll,j)=gradcorr(ll,j)+ghalf+ekont*derx(ll,4,1)
4770 gradcorr(ll,j1)=gradcorr(ll,j1)+ekont*derx(ll,5,1)
4771 cold ghalf=0.5d0*eel4*eij*gacont_hbr(ll,kk,k)
4772 ghalf=0.5d0*ggg2(ll)
4774 gradcorr(ll,k)=gradcorr(ll,k)+ghalf+ekont*derx(ll,2,2)
4775 gradcorr(ll,k+1)=gradcorr(ll,k+1)+ekont*derx(ll,3,2)
4776 gradcorr(ll,l)=gradcorr(ll,l)+ghalf+ekont*derx(ll,4,2)
4777 gradcorr(ll,l1)=gradcorr(ll,l1)+ekont*derx(ll,5,2)
4782 cold gradcorr(ll,m)=gradcorr(ll,m)+eel4*ekl*gacont_hbr(ll,jj,i)
4783 gradcorr(ll,m)=gradcorr(ll,m)+ggg1(ll)
4788 cold gradcorr(ll,m)=gradcorr(ll,m)+eel4*eij*gacont_hbr(ll,kk,k)
4789 gradcorr(ll,m)=gradcorr(ll,m)+ggg2(ll)
4795 gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,1)
4800 gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,2)
4804 cd write (2,*) iii,gcorr_loc(iii)
4808 cd write (2,*) 'ekont',ekont
4809 cd write (iout,*) 'eello4',ekont*eel4
4812 C---------------------------------------------------------------------------
4813 double precision function eello5(i,j,k,l,jj,kk)
4814 implicit real*8 (a-h,o-z)
4815 include 'DIMENSIONS'
4816 include 'COMMON.IOUNITS'
4817 include 'COMMON.CHAIN'
4818 include 'COMMON.DERIV'
4819 include 'COMMON.INTERACT'
4820 include 'COMMON.CONTACTS'
4821 include 'COMMON.TORSION'
4822 include 'COMMON.VAR'
4823 include 'COMMON.GEO'
4824 double precision pizda(2,2),auxmat(2,2),auxmat1(2,2),vv(2)
4825 double precision ggg1(3),ggg2(3)
4826 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4831 C /l\ / \ \ / \ / \ / C
4832 C / \ / \ \ / \ / \ / C
4833 C j| o |l1 | o | o| o | | o |o C
4834 C \ |/k\| |/ \| / |/ \| |/ \| C
4835 C \i/ \ / \ / / \ / \ C
4837 C (I) (II) (III) (IV) C
4839 C eello5_1 eello5_2 eello5_3 eello5_4 C
4841 C Antiparallel chains C
4844 C /j\ / \ \ / \ / \ / C
4845 C / \ / \ \ / \ / \ / C
4846 C j1| o |l | o | o| o | | o |o C
4847 C \ |/k\| |/ \| / |/ \| |/ \| C
4848 C \i/ \ / \ / / \ / \ C
4850 C (I) (II) (III) (IV) C
4852 C eello5_1 eello5_2 eello5_3 eello5_4 C
4854 C o denotes a local interaction, vertical lines an electrostatic interaction. C
4856 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4857 cd if (i.ne.2 .or. j.ne.6 .or. k.ne.3 .or. l.ne.5) then
4862 cd & 'EELLO5: Contacts have occurred for peptide groups',i,j,
4864 itk=itortyp(itype(k))
4865 itl=itortyp(itype(l))
4866 itj=itortyp(itype(j))
4871 cd call checkint5(i,j,k,l,jj,kk,eel5_1_num,eel5_2_num,
4872 cd & eel5_3_num,eel5_4_num)
4876 derx(lll,kkk,iii)=0.0d0
4880 cd eij=facont_hb(jj,i)
4881 cd ekl=facont_hb(kk,k)
4883 cd write (iout,*)'Contacts have occurred for peptide groups',
4884 cd & i,j,' fcont:',eij,' eij',' and ',k,l
4886 C Contribution from the graph I.
4887 cd write (2,*) 'AEA ',AEA(1,1,1),AEA(2,1,1),AEA(1,2,1),AEA(2,2,1)
4888 cd write (2,*) 'AEAb2',AEAb2(1,1,1),AEAb2(2,1,1)
4889 call transpose2(EUg(1,1,k),auxmat(1,1))
4890 call matmat2(AEA(1,1,1),auxmat(1,1),pizda(1,1))
4891 vv(1)=pizda(1,1)-pizda(2,2)
4892 vv(2)=pizda(1,2)+pizda(2,1)
4893 eello5_1=scalar2(AEAb2(1,1,1),Ub2(1,k))
4894 & +0.5d0*scalar2(vv(1),Dtobr2(1,i))
4896 C Explicit gradient in virtual-dihedral angles.
4897 if (i.gt.1) g_corr5_loc(i-1)=g_corr5_loc(i-1)
4898 & +ekont*(scalar2(AEAb2derg(1,2,1,1),Ub2(1,k))
4899 & +0.5d0*scalar2(vv(1),Dtobr2der(1,i)))
4900 call transpose2(EUgder(1,1,k),auxmat1(1,1))
4901 call matmat2(AEA(1,1,1),auxmat1(1,1),pizda(1,1))
4902 vv(1)=pizda(1,1)-pizda(2,2)
4903 vv(2)=pizda(1,2)+pizda(2,1)
4904 g_corr5_loc(k-1)=g_corr5_loc(k-1)
4905 & +ekont*(scalar2(AEAb2(1,1,1),Ub2der(1,k))
4906 & +0.5d0*scalar2(vv(1),Dtobr2(1,i)))
4907 call matmat2(AEAderg(1,1,1),auxmat(1,1),pizda(1,1))
4908 vv(1)=pizda(1,1)-pizda(2,2)
4909 vv(2)=pizda(1,2)+pizda(2,1)
4911 if (l.lt.nres-1) g_corr5_loc(l-1)=g_corr5_loc(l-1)
4912 & +ekont*(scalar2(AEAb2derg(1,1,1,1),Ub2(1,k))
4913 & +0.5d0*scalar2(vv(1),Dtobr2(1,i)))
4915 if (j.lt.nres-1) g_corr5_loc(j-1)=g_corr5_loc(j-1)
4916 & +ekont*(scalar2(AEAb2derg(1,1,1,1),Ub2(1,k))
4917 & +0.5d0*scalar2(vv(1),Dtobr2(1,i)))
4919 C Cartesian gradient
4923 call matmat2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1),
4925 vv(1)=pizda(1,1)-pizda(2,2)
4926 vv(2)=pizda(1,2)+pizda(2,1)
4927 derx(lll,kkk,iii)=derx(lll,kkk,iii)
4928 & +scalar2(AEAb2derx(1,lll,kkk,iii,1,1),Ub2(1,k))
4929 & +0.5d0*scalar2(vv(1),Dtobr2(1,i))
4936 C Contribution from graph II
4937 call transpose2(EE(1,1,itk),auxmat(1,1))
4938 call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1))
4939 vv(1)=pizda(1,1)+pizda(2,2)
4940 vv(2)=pizda(2,1)-pizda(1,2)
4941 eello5_2=scalar2(AEAb1(1,2,1),b1(1,itk))
4942 & -0.5d0*scalar2(vv(1),Ctobr(1,k))
4944 C Explicit gradient in virtual-dihedral angles.
4945 g_corr5_loc(k-1)=g_corr5_loc(k-1)
4946 & -0.5d0*ekont*scalar2(vv(1),Ctobrder(1,k))
4947 call matmat2(auxmat(1,1),AEAderg(1,1,1),pizda(1,1))
4948 vv(1)=pizda(1,1)+pizda(2,2)
4949 vv(2)=pizda(2,1)-pizda(1,2)
4951 g_corr5_loc(l-1)=g_corr5_loc(l-1)
4952 & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk))
4953 & -0.5d0*scalar2(vv(1),Ctobr(1,k)))
4955 g_corr5_loc(j-1)=g_corr5_loc(j-1)
4956 & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk))
4957 & -0.5d0*scalar2(vv(1),Ctobr(1,k)))
4959 C Cartesian gradient
4963 call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,1),
4965 vv(1)=pizda(1,1)+pizda(2,2)
4966 vv(2)=pizda(2,1)-pizda(1,2)
4967 derx(lll,kkk,iii)=derx(lll,kkk,iii)
4968 & +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,itk))
4969 & -0.5d0*scalar2(vv(1),Ctobr(1,k))
4978 C Parallel orientation
4979 C Contribution from graph III
4980 call transpose2(EUg(1,1,l),auxmat(1,1))
4981 call matmat2(AEA(1,1,2),auxmat(1,1),pizda(1,1))
4982 vv(1)=pizda(1,1)-pizda(2,2)
4983 vv(2)=pizda(1,2)+pizda(2,1)
4984 eello5_3=scalar2(AEAb2(1,1,2),Ub2(1,l))
4985 & +0.5d0*scalar2(vv(1),Dtobr2(1,j))
4987 C Explicit gradient in virtual-dihedral angles.
4988 g_corr5_loc(j-1)=g_corr5_loc(j-1)
4989 & +ekont*(scalar2(AEAb2derg(1,2,1,2),Ub2(1,l))
4990 & +0.5d0*scalar2(vv(1),Dtobr2der(1,j)))
4991 call matmat2(AEAderg(1,1,2),auxmat(1,1),pizda(1,1))
4992 vv(1)=pizda(1,1)-pizda(2,2)
4993 vv(2)=pizda(1,2)+pizda(2,1)
4994 g_corr5_loc(k-1)=g_corr5_loc(k-1)
4995 & +ekont*(scalar2(AEAb2derg(1,1,1,2),Ub2(1,l))
4996 & +0.5d0*scalar2(vv(1),Dtobr2(1,j)))
4997 call transpose2(EUgder(1,1,l),auxmat1(1,1))
4998 call matmat2(AEA(1,1,2),auxmat1(1,1),pizda(1,1))
4999 vv(1)=pizda(1,1)-pizda(2,2)
5000 vv(2)=pizda(1,2)+pizda(2,1)
5001 g_corr5_loc(l-1)=g_corr5_loc(l-1)
5002 & +ekont*(scalar2(AEAb2(1,1,2),Ub2der(1,l))
5003 & +0.5d0*scalar2(vv(1),Dtobr2(1,j)))
5004 C Cartesian gradient
5008 call matmat2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1),
5010 vv(1)=pizda(1,1)-pizda(2,2)
5011 vv(2)=pizda(1,2)+pizda(2,1)
5012 derx(lll,kkk,iii)=derx(lll,kkk,iii)
5013 & +scalar2(AEAb2derx(1,lll,kkk,iii,1,2),Ub2(1,l))
5014 & +0.5d0*scalar2(vv(1),Dtobr2(1,j))
5020 C Contribution from graph IV
5022 call transpose2(EE(1,1,itl),auxmat(1,1))
5023 call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
5024 vv(1)=pizda(1,1)+pizda(2,2)
5025 vv(2)=pizda(2,1)-pizda(1,2)
5026 eello5_4=scalar2(AEAb1(1,2,2),b1(1,itl))
5027 & -0.5d0*scalar2(vv(1),Ctobr(1,l))
5029 C Explicit gradient in virtual-dihedral angles.
5030 g_corr5_loc(l-1)=g_corr5_loc(l-1)
5031 & -0.5d0*ekont*scalar2(vv(1),Ctobrder(1,l))
5032 call matmat2(auxmat(1,1),AEAderg(1,1,2),pizda(1,1))
5033 vv(1)=pizda(1,1)+pizda(2,2)
5034 vv(2)=pizda(2,1)-pizda(1,2)
5035 g_corr5_loc(k-1)=g_corr5_loc(k-1)
5036 & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itl))
5037 & -0.5d0*scalar2(vv(1),Ctobr(1,l)))
5038 C Cartesian gradient
5042 call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,2),
5044 vv(1)=pizda(1,1)+pizda(2,2)
5045 vv(2)=pizda(2,1)-pizda(1,2)
5046 derx(lll,kkk,iii)=derx(lll,kkk,iii)
5047 & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itl))
5048 & -0.5d0*scalar2(vv(1),Ctobr(1,l))
5054 C Antiparallel orientation
5055 C Contribution from graph III
5057 call transpose2(EUg(1,1,j),auxmat(1,1))
5058 call matmat2(AEA(1,1,2),auxmat(1,1),pizda(1,1))
5059 vv(1)=pizda(1,1)-pizda(2,2)
5060 vv(2)=pizda(1,2)+pizda(2,1)
5061 eello5_3=scalar2(AEAb2(1,1,2),Ub2(1,j))
5062 & +0.5d0*scalar2(vv(1),Dtobr2(1,l))
5064 C Explicit gradient in virtual-dihedral angles.
5065 g_corr5_loc(l-1)=g_corr5_loc(l-1)
5066 & +ekont*(scalar2(AEAb2derg(1,2,1,2),Ub2(1,j))
5067 & +0.5d0*scalar2(vv(1),Dtobr2der(1,l)))
5068 call matmat2(AEAderg(1,1,2),auxmat(1,1),pizda(1,1))
5069 vv(1)=pizda(1,1)-pizda(2,2)
5070 vv(2)=pizda(1,2)+pizda(2,1)
5071 g_corr5_loc(k-1)=g_corr5_loc(k-1)
5072 & +ekont*(scalar2(AEAb2derg(1,1,1,2),Ub2(1,j))
5073 & +0.5d0*scalar2(vv(1),Dtobr2(1,l)))
5074 call transpose2(EUgder(1,1,j),auxmat1(1,1))
5075 call matmat2(AEA(1,1,2),auxmat1(1,1),pizda(1,1))
5076 vv(1)=pizda(1,1)-pizda(2,2)
5077 vv(2)=pizda(1,2)+pizda(2,1)
5078 g_corr5_loc(j-1)=g_corr5_loc(j-1)
5079 & +ekont*(scalar2(AEAb2(1,1,2),Ub2der(1,j))
5080 & +0.5d0*scalar2(vv(1),Dtobr2(1,l)))
5081 C Cartesian gradient
5085 call matmat2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1),
5087 vv(1)=pizda(1,1)-pizda(2,2)
5088 vv(2)=pizda(1,2)+pizda(2,1)
5089 derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)
5090 & +scalar2(AEAb2derx(1,lll,kkk,iii,1,2),Ub2(1,j))
5091 & +0.5d0*scalar2(vv(1),Dtobr2(1,l))
5097 C Contribution from graph IV
5099 call transpose2(EE(1,1,itj),auxmat(1,1))
5100 call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
5101 vv(1)=pizda(1,1)+pizda(2,2)
5102 vv(2)=pizda(2,1)-pizda(1,2)
5103 eello5_4=scalar2(AEAb1(1,2,2),b1(1,itj))
5104 & -0.5d0*scalar2(vv(1),Ctobr(1,j))
5106 C Explicit gradient in virtual-dihedral angles.
5107 g_corr5_loc(j-1)=g_corr5_loc(j-1)
5108 & -0.5d0*ekont*scalar2(vv(1),Ctobrder(1,j))
5109 call matmat2(auxmat(1,1),AEAderg(1,1,2),pizda(1,1))
5110 vv(1)=pizda(1,1)+pizda(2,2)
5111 vv(2)=pizda(2,1)-pizda(1,2)
5112 g_corr5_loc(k-1)=g_corr5_loc(k-1)
5113 & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itj))
5114 & -0.5d0*scalar2(vv(1),Ctobr(1,j)))
5115 C Cartesian gradient
5119 call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,2),
5121 vv(1)=pizda(1,1)+pizda(2,2)
5122 vv(2)=pizda(2,1)-pizda(1,2)
5123 derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)
5124 & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itj))
5125 & -0.5d0*scalar2(vv(1),Ctobr(1,j))
5132 eel5=eello5_1+eello5_2+eello5_3+eello5_4
5133 cd if (i.eq.2 .and. j.eq.8 .and. k.eq.3 .and. l.eq.7) then
5134 cd write (2,*) 'ijkl',i,j,k,l
5135 cd write (2,*) 'eello5_1',eello5_1,' eello5_2',eello5_2,
5136 cd & ' eello5_3',eello5_3,' eello5_4',eello5_4
5138 cd write(iout,*) 'eello5_1',eello5_1,' eel5_1_num',16*eel5_1_num
5139 cd write(iout,*) 'eello5_2',eello5_2,' eel5_2_num',16*eel5_2_num
5140 cd write(iout,*) 'eello5_3',eello5_3,' eel5_3_num',16*eel5_3_num
5141 cd write(iout,*) 'eello5_4',eello5_4,' eel5_4_num',16*eel5_4_num
5143 if (j.lt.nres-1) then
5150 if (l.lt.nres-1) then
5160 cd write (2,*) 'eij',eij,' ekl',ekl,' ekont',ekont
5162 ggg1(ll)=eel5*g_contij(ll,1)
5163 ggg2(ll)=eel5*g_contij(ll,2)
5164 cold ghalf=0.5d0*eel5*ekl*gacont_hbr(ll,jj,i)
5165 ghalf=0.5d0*ggg1(ll)
5167 gradcorr5(ll,i)=gradcorr5(ll,i)+ghalf+ekont*derx(ll,2,1)
5168 gradcorr5(ll,i+1)=gradcorr5(ll,i+1)+ekont*derx(ll,3,1)
5169 gradcorr5(ll,j)=gradcorr5(ll,j)+ghalf+ekont*derx(ll,4,1)
5170 gradcorr5(ll,j1)=gradcorr5(ll,j1)+ekont*derx(ll,5,1)
5171 cold ghalf=0.5d0*eel5*eij*gacont_hbr(ll,kk,k)
5172 ghalf=0.5d0*ggg2(ll)
5174 gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2)
5175 gradcorr5(ll,k+1)=gradcorr5(ll,k+1)+ekont*derx(ll,3,2)
5176 gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2)
5177 gradcorr5(ll,l1)=gradcorr5(ll,l1)+ekont*derx(ll,5,2)
5182 cold gradcorr5(ll,m)=gradcorr5(ll,m)+eel5*ekl*gacont_hbr(ll,jj,i)
5183 gradcorr5(ll,m)=gradcorr5(ll,m)+ggg1(ll)
5188 cold gradcorr5(ll,m)=gradcorr5(ll,m)+eel5*eij*gacont_hbr(ll,kk,k)
5189 gradcorr5(ll,m)=gradcorr5(ll,m)+ggg2(ll)
5195 gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,1)
5200 gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,2)
5204 cd write (2,*) iii,g_corr5_loc(iii)
5208 cd write (2,*) 'ekont',ekont
5209 cd write (iout,*) 'eello5',ekont*eel5
5212 c--------------------------------------------------------------------------
5213 double precision function eello6(i,j,k,l,jj,kk)
5214 implicit real*8 (a-h,o-z)
5215 include 'DIMENSIONS'
5216 include 'DIMENSIONS.ZSCOPT'
5217 include 'COMMON.IOUNITS'
5218 include 'COMMON.CHAIN'
5219 include 'COMMON.DERIV'
5220 include 'COMMON.INTERACT'
5221 include 'COMMON.CONTACTS'
5222 include 'COMMON.TORSION'
5223 include 'COMMON.VAR'
5224 include 'COMMON.GEO'
5225 include 'COMMON.FFIELD'
5226 double precision ggg1(3),ggg2(3)
5227 cd if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then
5232 cd & 'EELLO6: Contacts have occurred for peptide groups',i,j,
5240 cd call checkint6(i,j,k,l,jj,kk,eel6_1_num,eel6_2_num,
5241 cd & eel6_3_num,eel6_4_num,eel6_5_num,eel6_6_num)
5245 derx(lll,kkk,iii)=0.0d0
5249 cd eij=facont_hb(jj,i)
5250 cd ekl=facont_hb(kk,k)
5256 eello6_1=eello6_graph1(i,j,k,l,1,.false.)
5257 eello6_2=eello6_graph1(j,i,l,k,2,.false.)
5258 eello6_3=eello6_graph2(i,j,k,l,jj,kk,.false.)
5259 eello6_4=eello6_graph4(i,j,k,l,jj,kk,1,.false.)
5260 eello6_5=eello6_graph4(j,i,l,k,jj,kk,2,.false.)
5261 eello6_6=eello6_graph3(i,j,k,l,jj,kk,.false.)
5263 eello6_1=eello6_graph1(i,j,k,l,1,.false.)
5264 eello6_2=eello6_graph1(l,k,j,i,2,.true.)
5265 eello6_3=eello6_graph2(i,l,k,j,jj,kk,.true.)
5266 eello6_4=eello6_graph4(i,j,k,l,jj,kk,1,.false.)
5267 if (wturn6.eq.0.0d0 .or. j.ne.i+4) then
5268 eello6_5=eello6_graph4(l,k,j,i,kk,jj,2,.true.)
5272 eello6_6=eello6_graph3(i,l,k,j,jj,kk,.true.)
5274 C If turn contributions are considered, they will be handled separately.
5275 eel6=eello6_1+eello6_2+eello6_3+eello6_4+eello6_5+eello6_6
5276 cd write(iout,*) 'eello6_1',eello6_1,' eel6_1_num',16*eel6_1_num
5277 cd write(iout,*) 'eello6_2',eello6_2,' eel6_2_num',16*eel6_2_num
5278 cd write(iout,*) 'eello6_3',eello6_3,' eel6_3_num',16*eel6_3_num
5279 cd write(iout,*) 'eello6_4',eello6_4,' eel6_4_num',16*eel6_4_num
5280 cd write(iout,*) 'eello6_5',eello6_5,' eel6_5_num',16*eel6_5_num
5281 cd write(iout,*) 'eello6_6',eello6_6,' eel6_6_num',16*eel6_6_num
5284 if (j.lt.nres-1) then
5291 if (l.lt.nres-1) then
5299 ggg1(ll)=eel6*g_contij(ll,1)
5300 ggg2(ll)=eel6*g_contij(ll,2)
5301 cold ghalf=0.5d0*eel6*ekl*gacont_hbr(ll,jj,i)
5302 ghalf=0.5d0*ggg1(ll)
5304 gradcorr6(ll,i)=gradcorr6(ll,i)+ghalf+ekont*derx(ll,2,1)
5305 gradcorr6(ll,i+1)=gradcorr6(ll,i+1)+ekont*derx(ll,3,1)
5306 gradcorr6(ll,j)=gradcorr6(ll,j)+ghalf+ekont*derx(ll,4,1)
5307 gradcorr6(ll,j1)=gradcorr6(ll,j1)+ekont*derx(ll,5,1)
5308 ghalf=0.5d0*ggg2(ll)
5309 cold ghalf=0.5d0*eel6*eij*gacont_hbr(ll,kk,k)
5311 gradcorr6(ll,k)=gradcorr6(ll,k)+ghalf+ekont*derx(ll,2,2)
5312 gradcorr6(ll,k+1)=gradcorr6(ll,k+1)+ekont*derx(ll,3,2)
5313 gradcorr6(ll,l)=gradcorr6(ll,l)+ghalf+ekont*derx(ll,4,2)
5314 gradcorr6(ll,l1)=gradcorr6(ll,l1)+ekont*derx(ll,5,2)
5319 cold gradcorr6(ll,m)=gradcorr6(ll,m)+eel6*ekl*gacont_hbr(ll,jj,i)
5320 gradcorr6(ll,m)=gradcorr6(ll,m)+ggg1(ll)
5325 cold gradcorr6(ll,m)=gradcorr6(ll,m)+eel6*eij*gacont_hbr(ll,kk,k)
5326 gradcorr6(ll,m)=gradcorr6(ll,m)+ggg2(ll)
5332 gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,1)
5337 gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,2)
5341 cd write (2,*) iii,g_corr6_loc(iii)
5345 cd write (2,*) 'ekont',ekont
5346 cd write (iout,*) 'eello6',ekont*eel6
5349 c--------------------------------------------------------------------------
5350 double precision function eello6_graph1(i,j,k,l,imat,swap)
5351 implicit real*8 (a-h,o-z)
5352 include 'DIMENSIONS'
5353 include 'COMMON.IOUNITS'
5354 include 'COMMON.CHAIN'
5355 include 'COMMON.DERIV'
5356 include 'COMMON.INTERACT'
5357 include 'COMMON.CONTACTS'
5358 include 'COMMON.TORSION'
5359 include 'COMMON.VAR'
5360 include 'COMMON.GEO'
5361 double precision vv(2),vv1(2),pizda(2,2),auxmat(2,2),pizda1(2,2)
5365 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5367 C Parallel Antiparallel
5373 C \ j|/k\| / \ |/k\|l /
5378 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5379 itk=itortyp(itype(k))
5380 s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i))
5381 s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k))
5382 s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k))
5383 call transpose2(EUgC(1,1,k),auxmat(1,1))
5384 call matmat2(AEA(1,1,imat),auxmat(1,1),pizda1(1,1))
5385 vv1(1)=pizda1(1,1)-pizda1(2,2)
5386 vv1(2)=pizda1(1,2)+pizda1(2,1)
5387 s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i))
5388 vv(1)=AEAb1(1,2,imat)*b1(1,itk)-AEAb1(2,2,imat)*b1(2,itk)
5389 vv(2)=AEAb1(1,2,imat)*b1(2,itk)+AEAb1(2,2,imat)*b1(1,itk)
5390 s5=scalar2(vv(1),Dtobr2(1,i))
5391 cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5
5392 eello6_graph1=-0.5d0*(s1+s2+s3+s4+s5)
5393 if (.not. calc_grad) return
5394 if (i.gt.1) g_corr6_loc(i-1)=g_corr6_loc(i-1)
5395 & -0.5d0*ekont*(scalar2(AEAb1(1,2,imat),CUgb2der(1,i))
5396 & -scalar2(AEAb2derg(1,2,1,imat),Ug2Db1t(1,k))
5397 & +scalar2(AEAb2derg(1,2,1,imat),CUgb2(1,k))
5398 & +0.5d0*scalar2(vv1(1),Dtobr2der(1,i))
5399 & +scalar2(vv(1),Dtobr2der(1,i)))
5400 call matmat2(AEAderg(1,1,imat),auxmat(1,1),pizda1(1,1))
5401 vv1(1)=pizda1(1,1)-pizda1(2,2)
5402 vv1(2)=pizda1(1,2)+pizda1(2,1)
5403 vv(1)=AEAb1derg(1,2,imat)*b1(1,itk)-AEAb1derg(2,2,imat)*b1(2,itk)
5404 vv(2)=AEAb1derg(1,2,imat)*b1(2,itk)+AEAb1derg(2,2,imat)*b1(1,itk)
5406 g_corr6_loc(l-1)=g_corr6_loc(l-1)
5407 & +ekont*(-0.5d0*(scalar2(AEAb1derg(1,2,imat),CUgb2(1,i))
5408 & -scalar2(AEAb2derg(1,1,1,imat),Ug2Db1t(1,k))
5409 & +scalar2(AEAb2derg(1,1,1,imat),CUgb2(1,k))
5410 & +0.5d0*scalar2(vv1(1),Dtobr2(1,i))+scalar2(vv(1),Dtobr2(1,i))))
5412 g_corr6_loc(j-1)=g_corr6_loc(j-1)
5413 & +ekont*(-0.5d0*(scalar2(AEAb1derg(1,2,imat),CUgb2(1,i))
5414 & -scalar2(AEAb2derg(1,1,1,imat),Ug2Db1t(1,k))
5415 & +scalar2(AEAb2derg(1,1,1,imat),CUgb2(1,k))
5416 & +0.5d0*scalar2(vv1(1),Dtobr2(1,i))+scalar2(vv(1),Dtobr2(1,i))))
5418 call transpose2(EUgCder(1,1,k),auxmat(1,1))
5419 call matmat2(AEA(1,1,imat),auxmat(1,1),pizda1(1,1))
5420 vv1(1)=pizda1(1,1)-pizda1(2,2)
5421 vv1(2)=pizda1(1,2)+pizda1(2,1)
5422 if (k.gt.1) g_corr6_loc(k-1)=g_corr6_loc(k-1)
5423 & +ekont*(-0.5d0*(-scalar2(AEAb2(1,1,imat),Ug2Db1tder(1,k))
5424 & +scalar2(AEAb2(1,1,imat),CUgb2der(1,k))
5425 & +0.5d0*scalar2(vv1(1),Dtobr2(1,i))))
5434 s1= scalar2(AEAb1derx(1,lll,kkk,iii,2,imat),CUgb2(1,i))
5435 s2=-scalar2(AEAb2derx(1,lll,kkk,iii,1,imat),Ug2Db1t(1,k))
5436 s3= scalar2(AEAb2derx(1,lll,kkk,iii,1,imat),CUgb2(1,k))
5437 call transpose2(EUgC(1,1,k),auxmat(1,1))
5438 call matmat2(AEAderx(1,1,lll,kkk,iii,imat),auxmat(1,1),
5440 vv1(1)=pizda1(1,1)-pizda1(2,2)
5441 vv1(2)=pizda1(1,2)+pizda1(2,1)
5442 s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i))
5443 vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,itk)
5444 & -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,itk)
5445 vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,itk)
5446 & +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,itk)
5447 s5=scalar2(vv(1),Dtobr2(1,i))
5448 derx(lll,kkk,ind)=derx(lll,kkk,ind)-0.5d0*(s1+s2+s3+s4+s5)
5454 c----------------------------------------------------------------------------
5455 double precision function eello6_graph2(i,j,k,l,jj,kk,swap)
5456 implicit real*8 (a-h,o-z)
5457 include 'DIMENSIONS'
5458 include 'COMMON.IOUNITS'
5459 include 'COMMON.CHAIN'
5460 include 'COMMON.DERIV'
5461 include 'COMMON.INTERACT'
5462 include 'COMMON.CONTACTS'
5463 include 'COMMON.TORSION'
5464 include 'COMMON.VAR'
5465 include 'COMMON.GEO'
5467 double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2),
5468 & auxvec1(2),auxvec2(1),auxmat1(2,2)
5471 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5473 C Parallel Antiparallel
5484 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5485 cd write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l
5486 C AL 7/4/01 s1 would occur in the sixth-order moment,
5487 C but not in a cluster cumulant
5489 s1=dip(1,jj,i)*dip(1,kk,k)
5491 call matvec2(ADtEA1(1,1,1),Ub2(1,k),auxvec(1))
5492 s2=-0.5d0*scalar2(Ub2(1,i),auxvec(1))
5493 call matvec2(ADtEA(1,1,2),Ub2(1,l),auxvec1(1))
5494 s3=-0.5d0*scalar2(Ub2(1,j),auxvec1(1))
5495 call transpose2(EUg(1,1,k),auxmat(1,1))
5496 call matmat2(ADtEA1(1,1,1),auxmat(1,1),pizda(1,1))
5497 vv(1)=pizda(1,1)-pizda(2,2)
5498 vv(2)=pizda(1,2)+pizda(2,1)
5499 s4=-0.25d0*scalar2(vv(1),Dtobr2(1,i))
5500 cd write (2,*) 'eello6_graph2:','s1',s1,' s2',s2,' s3',s3,' s4',s4
5502 eello6_graph2=-(s1+s2+s3+s4)
5504 eello6_graph2=-(s2+s3+s4)
5507 if (.not. calc_grad) return
5508 C Derivatives in gamma(i-1)
5511 s1=dipderg(1,jj,i)*dip(1,kk,k)
5513 s2=-0.5d0*scalar2(Ub2der(1,i),auxvec(1))
5514 call matvec2(ADtEAderg(1,1,1,2),Ub2(1,l),auxvec2(1))
5515 s3=-0.5d0*scalar2(Ub2(1,j),auxvec2(1))
5516 s4=-0.25d0*scalar2(vv(1),Dtobr2der(1,i))
5518 g_corr6_loc(i-1)=g_corr6_loc(i-1)-ekont*(s1+s2+s3+s4)
5520 g_corr6_loc(i-1)=g_corr6_loc(i-1)-ekont*(s2+s3+s4)
5522 c g_corr6_loc(i-1)=g_corr6_loc(i-1)-s3
5524 C Derivatives in gamma(k-1)
5526 s1=dip(1,jj,i)*dipderg(1,kk,k)
5528 call matvec2(ADtEA1(1,1,1),Ub2der(1,k),auxvec2(1))
5529 s2=-0.5d0*scalar2(Ub2(1,i),auxvec2(1))
5530 call matvec2(ADtEAderg(1,1,2,2),Ub2(1,l),auxvec2(1))
5531 s3=-0.5d0*scalar2(Ub2(1,j),auxvec2(1))
5532 call transpose2(EUgder(1,1,k),auxmat1(1,1))
5533 call matmat2(ADtEA1(1,1,1),auxmat1(1,1),pizda(1,1))
5534 vv(1)=pizda(1,1)-pizda(2,2)
5535 vv(2)=pizda(1,2)+pizda(2,1)
5536 s4=-0.25d0*scalar2(vv(1),Dtobr2(1,i))
5538 g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s1+s2+s3+s4)
5540 g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s2+s3+s4)
5542 c g_corr6_loc(k-1)=g_corr6_loc(k-1)-s3
5543 C Derivatives in gamma(j-1) or gamma(l-1)
5546 s1=dipderg(3,jj,i)*dip(1,kk,k)
5548 call matvec2(ADtEA1derg(1,1,1,1),Ub2(1,k),auxvec2(1))
5549 s2=-0.5d0*scalar2(Ub2(1,i),auxvec2(1))
5550 s3=-0.5d0*scalar2(Ub2der(1,j),auxvec1(1))
5551 call matmat2(ADtEA1derg(1,1,1,1),auxmat(1,1),pizda(1,1))
5552 vv(1)=pizda(1,1)-pizda(2,2)
5553 vv(2)=pizda(1,2)+pizda(2,1)
5554 s4=-0.25d0*scalar2(vv(1),Dtobr2(1,i))
5557 g_corr6_loc(l-1)=g_corr6_loc(l-1)-ekont*s1
5559 g_corr6_loc(j-1)=g_corr6_loc(j-1)-ekont*s1
5562 g_corr6_loc(j-1)=g_corr6_loc(j-1)-ekont*(s2+s3+s4)
5563 c g_corr6_loc(j-1)=g_corr6_loc(j-1)-s3
5565 C Derivatives in gamma(l-1) or gamma(j-1)
5568 s1=dip(1,jj,i)*dipderg(3,kk,k)
5570 call matvec2(ADtEA1derg(1,1,2,1),Ub2(1,k),auxvec2(1))
5571 s2=-0.5d0*scalar2(Ub2(1,i),auxvec2(1))
5572 call matvec2(ADtEA(1,1,2),Ub2der(1,l),auxvec2(1))
5573 s3=-0.5d0*scalar2(Ub2(1,j),auxvec2(1))
5574 call matmat2(ADtEA1derg(1,1,2,1),auxmat(1,1),pizda(1,1))
5575 vv(1)=pizda(1,1)-pizda(2,2)
5576 vv(2)=pizda(1,2)+pizda(2,1)
5577 s4=-0.25d0*scalar2(vv(1),Dtobr2(1,i))
5580 g_corr6_loc(j-1)=g_corr6_loc(j-1)-ekont*s1
5582 g_corr6_loc(l-1)=g_corr6_loc(l-1)-ekont*s1
5585 g_corr6_loc(l-1)=g_corr6_loc(l-1)-ekont*(s2+s3+s4)
5586 c g_corr6_loc(l-1)=g_corr6_loc(l-1)-s3
5588 C Cartesian derivatives.
5590 write (2,*) 'In eello6_graph2'
5592 write (2,*) 'iii=',iii
5594 write (2,*) 'kkk=',kkk
5596 write (2,'(3(2f10.5),5x)')
5597 & ((ADtEA1derx(jjj,mmm,lll,kkk,iii,1),mmm=1,2),lll=1,3)
5607 s1=dipderx(lll,kkk,1,jj,i)*dip(1,kk,k)
5609 s1=dip(1,jj,i)*dipderx(lll,kkk,1,kk,k)
5612 call matvec2(ADtEA1derx(1,1,lll,kkk,iii,1),Ub2(1,k),
5614 s2=-0.5d0*scalar2(Ub2(1,i),auxvec(1))
5615 call matvec2(ADtEAderx(1,1,lll,kkk,iii,2),Ub2(1,l),
5617 s3=-0.5d0*scalar2(Ub2(1,j),auxvec(1))
5618 call transpose2(EUg(1,1,k),auxmat(1,1))
5619 call matmat2(ADtEA1derx(1,1,lll,kkk,iii,1),auxmat(1,1),
5621 vv(1)=pizda(1,1)-pizda(2,2)
5622 vv(2)=pizda(1,2)+pizda(2,1)
5623 s4=-0.25d0*scalar2(vv(1),Dtobr2(1,i))
5624 cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4',s4
5626 derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s1+s2+s4)
5628 derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s2+s4)
5631 derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)-s3
5633 derx(lll,kkk,iii)=derx(lll,kkk,iii)-s3
5640 c----------------------------------------------------------------------------
5641 double precision function eello6_graph3(i,j,k,l,jj,kk,swap)
5642 implicit real*8 (a-h,o-z)
5643 include 'DIMENSIONS'
5644 include 'COMMON.IOUNITS'
5645 include 'COMMON.CHAIN'
5646 include 'COMMON.DERIV'
5647 include 'COMMON.INTERACT'
5648 include 'COMMON.CONTACTS'
5649 include 'COMMON.TORSION'
5650 include 'COMMON.VAR'
5651 include 'COMMON.GEO'
5652 double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
5654 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5656 C Parallel Antiparallel
5667 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5669 C 4/7/01 AL Component s1 was removed, because it pertains to the respective
5670 C energy moment and not to the cluster cumulant.
5671 iti=itortyp(itype(i))
5672 if (j.lt.nres-1) then
5673 itj1=itortyp(itype(j+1))
5677 itk=itortyp(itype(k))
5678 itk1=itortyp(itype(k+1))
5679 if (l.lt.nres-1) then
5680 itl1=itortyp(itype(l+1))
5685 s1=dip(4,jj,i)*dip(4,kk,k)
5687 call matvec2(AECA(1,1,1),b1(1,itk1),auxvec(1))
5688 s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
5689 call matvec2(AECA(1,1,2),b1(1,itl1),auxvec(1))
5690 s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
5691 call transpose2(EE(1,1,itk),auxmat(1,1))
5692 call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1))
5693 vv(1)=pizda(1,1)+pizda(2,2)
5694 vv(2)=pizda(2,1)-pizda(1,2)
5695 s4=-0.25d0*scalar2(vv(1),Ctobr(1,k))
5696 cd write (2,*) 'eello6_graph3:','s1',s1,' s2',s2,' s3',s3,' s4',s4
5698 eello6_graph3=-(s1+s2+s3+s4)
5700 eello6_graph3=-(s2+s3+s4)
5703 if (.not. calc_grad) return
5704 C Derivatives in gamma(k-1)
5705 call matvec2(AECAderg(1,1,2),b1(1,itl1),auxvec(1))
5706 s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
5707 s4=-0.25d0*scalar2(vv(1),Ctobrder(1,k))
5708 g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s3+s4)
5709 C Derivatives in gamma(l-1)
5710 call matvec2(AECAderg(1,1,1),b1(1,itk1),auxvec(1))
5711 s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
5712 call matmat2(auxmat(1,1),AECAderg(1,1,1),pizda(1,1))
5713 vv(1)=pizda(1,1)+pizda(2,2)
5714 vv(2)=pizda(2,1)-pizda(1,2)
5715 s4=-0.25d0*scalar2(vv(1),Ctobr(1,k))
5716 g_corr6_loc(l-1)=g_corr6_loc(l-1)-ekont*(s2+s4)
5717 C Cartesian derivatives.
5723 s1=dipderx(lll,kkk,4,jj,i)*dip(4,kk,k)
5725 s1=dip(4,jj,i)*dipderx(lll,kkk,4,kk,k)
5728 call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
5730 s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
5731 call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,itl1),
5733 s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
5734 call matmat2(auxmat(1,1),AECAderx(1,1,lll,kkk,iii,1),
5736 vv(1)=pizda(1,1)+pizda(2,2)
5737 vv(2)=pizda(2,1)-pizda(1,2)
5738 s4=-0.25d0*scalar2(vv(1),Ctobr(1,k))
5740 derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s1+s2+s4)
5742 derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s2+s4)
5745 derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)-s3
5747 derx(lll,kkk,iii)=derx(lll,kkk,iii)-s3
5749 c derx(lll,kkk,iii)=derx(lll,kkk,iii)-s4
5755 c----------------------------------------------------------------------------
5756 double precision function eello6_graph4(i,j,k,l,jj,kk,imat,swap)
5757 implicit real*8 (a-h,o-z)
5758 include 'DIMENSIONS'
5759 include 'DIMENSIONS.ZSCOPT'
5760 include 'COMMON.IOUNITS'
5761 include 'COMMON.CHAIN'
5762 include 'COMMON.DERIV'
5763 include 'COMMON.INTERACT'
5764 include 'COMMON.CONTACTS'
5765 include 'COMMON.TORSION'
5766 include 'COMMON.VAR'
5767 include 'COMMON.GEO'
5768 include 'COMMON.FFIELD'
5769 double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2),
5770 & auxvec1(2),auxmat1(2,2)
5772 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5774 C Parallel Antiparallel
5785 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5787 C 4/7/01 AL Component s1 was removed, because it pertains to the respective
5788 C energy moment and not to the cluster cumulant.
5789 cd write (2,*) 'eello_graph4: wturn6',wturn6
5790 iti=itortyp(itype(i))
5791 itj=itortyp(itype(j))
5792 if (j.lt.nres-1) then
5793 itj1=itortyp(itype(j+1))
5797 itk=itortyp(itype(k))
5798 if (k.lt.nres-1) then
5799 itk1=itortyp(itype(k+1))
5803 itl=itortyp(itype(l))
5804 if (l.lt.nres-1) then
5805 itl1=itortyp(itype(l+1))
5809 cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l
5810 cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk,
5811 cd & ' itl',itl,' itl1',itl1
5814 s1=dip(3,jj,i)*dip(3,kk,k)
5816 s1=dip(2,jj,j)*dip(2,kk,l)
5819 call matvec2(AECA(1,1,imat),Ub2(1,k),auxvec(1))
5820 s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
5822 call matvec2(ADtEA1(1,1,3-imat),b1(1,itj1),auxvec1(1))
5823 s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
5825 call matvec2(ADtEA1(1,1,3-imat),b1(1,itl1),auxvec1(1))
5826 s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
5828 call transpose2(EUg(1,1,k),auxmat(1,1))
5829 call matmat2(AECA(1,1,imat),auxmat(1,1),pizda(1,1))
5830 vv(1)=pizda(1,1)-pizda(2,2)
5831 vv(2)=pizda(2,1)+pizda(1,2)
5832 s4=0.25d0*scalar2(vv(1),Dtobr2(1,i))
5833 cd write (2,*) 'eello6_graph4:','s1',s1,' s2',s2,' s3',s3,' s4',s4
5835 eello6_graph4=-(s1+s2+s3+s4)
5837 eello6_graph4=-(s2+s3+s4)
5839 if (.not. calc_grad) return
5840 C Derivatives in gamma(i-1)
5844 s1=dipderg(2,jj,i)*dip(3,kk,k)
5846 s1=dipderg(4,jj,j)*dip(2,kk,l)
5849 s2=0.5d0*scalar2(Ub2der(1,i),auxvec(1))
5851 call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itj1),auxvec1(1))
5852 s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
5854 call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itl1),auxvec1(1))
5855 s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
5857 s4=0.25d0*scalar2(vv(1),Dtobr2der(1,i))
5858 if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then
5859 cd write (2,*) 'turn6 derivatives'
5861 gel_loc_turn6(i-1)=gel_loc_turn6(i-1)-ekont*(s1+s2+s3+s4)
5863 gel_loc_turn6(i-1)=gel_loc_turn6(i-1)-ekont*(s2+s3+s4)
5867 g_corr6_loc(i-1)=g_corr6_loc(i-1)-ekont*(s1+s2+s3+s4)
5869 g_corr6_loc(i-1)=g_corr6_loc(i-1)-ekont*(s2+s3+s4)
5873 C Derivatives in gamma(k-1)
5876 s1=dip(3,jj,i)*dipderg(2,kk,k)
5878 s1=dip(2,jj,j)*dipderg(4,kk,l)
5881 call matvec2(AECA(1,1,imat),Ub2der(1,k),auxvec1(1))
5882 s2=0.5d0*scalar2(Ub2(1,i),auxvec1(1))
5884 call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itj1),auxvec1(1))
5885 s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
5887 call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itl1),auxvec1(1))
5888 s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
5890 call transpose2(EUgder(1,1,k),auxmat1(1,1))
5891 call matmat2(AECA(1,1,imat),auxmat1(1,1),pizda(1,1))
5892 vv(1)=pizda(1,1)-pizda(2,2)
5893 vv(2)=pizda(2,1)+pizda(1,2)
5894 s4=0.25d0*scalar2(vv(1),Dtobr2(1,i))
5895 if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then
5897 gel_loc_turn6(k-1)=gel_loc_turn6(k-1)-ekont*(s1+s2+s3+s4)
5899 gel_loc_turn6(k-1)=gel_loc_turn6(k-1)-ekont*(s2+s3+s4)
5903 g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s1+s2+s3+s4)
5905 g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s2+s3+s4)
5908 C Derivatives in gamma(j-1) or gamma(l-1)
5909 if (l.eq.j+1 .and. l.gt.1) then
5910 call matvec2(AECAderg(1,1,imat),Ub2(1,k),auxvec(1))
5911 s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
5912 call matmat2(AECAderg(1,1,imat),auxmat(1,1),pizda(1,1))
5913 vv(1)=pizda(1,1)-pizda(2,2)
5914 vv(2)=pizda(2,1)+pizda(1,2)
5915 s4=0.25d0*scalar2(vv(1),Dtobr2(1,i))
5916 g_corr6_loc(l-1)=g_corr6_loc(l-1)-ekont*(s2+s4)
5917 else if (j.gt.1) then
5918 call matvec2(AECAderg(1,1,imat),Ub2(1,k),auxvec(1))
5919 s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
5920 call matmat2(AECAderg(1,1,imat),auxmat(1,1),pizda(1,1))
5921 vv(1)=pizda(1,1)-pizda(2,2)
5922 vv(2)=pizda(2,1)+pizda(1,2)
5923 s4=0.25d0*scalar2(vv(1),Dtobr2(1,i))
5924 if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then
5925 gel_loc_turn6(j-1)=gel_loc_turn6(j-1)-ekont*(s2+s4)
5927 g_corr6_loc(j-1)=g_corr6_loc(j-1)-ekont*(s2+s4)
5930 C Cartesian derivatives.
5937 s1=dipderx(lll,kkk,3,jj,i)*dip(3,kk,k)
5939 s1=dipderx(lll,kkk,2,jj,j)*dip(2,kk,l)
5943 s1=dip(3,jj,i)*dipderx(lll,kkk,3,kk,k)
5945 s1=dip(2,jj,j)*dipderx(lll,kkk,2,kk,l)
5949 call matvec2(AECAderx(1,1,lll,kkk,iii,imat),Ub2(1,k),
5951 s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
5953 call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat),
5954 & b1(1,itj1),auxvec(1))
5955 s3=-0.5d0*scalar2(b1(1,itj),auxvec(1))
5957 call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat),
5958 & b1(1,itl1),auxvec(1))
5959 s3=-0.5d0*scalar2(b1(1,itl),auxvec(1))
5961 call matmat2(AECAderx(1,1,lll,kkk,iii,imat),auxmat(1,1),
5963 vv(1)=pizda(1,1)-pizda(2,2)
5964 vv(2)=pizda(2,1)+pizda(1,2)
5965 s4=0.25d0*scalar2(vv(1),Dtobr2(1,i))
5967 if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then
5969 derx_turn(lll,kkk,3-iii)=derx_turn(lll,kkk,3-iii)
5972 derx_turn(lll,kkk,3-iii)=derx_turn(lll,kkk,3-iii)
5975 derx_turn(lll,kkk,iii)=derx_turn(lll,kkk,iii)-s3
5978 derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)-(s1+s2+s4)
5980 derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)-(s2+s4)
5982 derx(lll,kkk,iii)=derx(lll,kkk,iii)-s3
5986 derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s1+s2+s4)
5988 derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s2+s4)
5991 derx(lll,kkk,iii)=derx(lll,kkk,iii)-s3
5993 derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)-s3
6001 c----------------------------------------------------------------------------
6002 double precision function eello_turn6(i,jj,kk)
6003 implicit real*8 (a-h,o-z)
6004 include 'DIMENSIONS'
6005 include 'COMMON.IOUNITS'
6006 include 'COMMON.CHAIN'
6007 include 'COMMON.DERIV'
6008 include 'COMMON.INTERACT'
6009 include 'COMMON.CONTACTS'
6010 include 'COMMON.TORSION'
6011 include 'COMMON.VAR'
6012 include 'COMMON.GEO'
6013 double precision vtemp1(2),vtemp2(2),vtemp3(2),vtemp4(2),
6014 & atemp(2,2),auxmat(2,2),achuj_temp(2,2),gtemp(2,2),gvec(2),
6016 double precision vtemp1d(2),vtemp2d(2),vtemp3d(2),vtemp4d(2),
6017 & atempd(2,2),auxmatd(2,2),achuj_tempd(2,2),gtempd(2,2),gvecd(2)
6018 C 4/7/01 AL Components s1, s8, and s13 were removed, because they pertain to
6019 C the respective energy moment and not to the cluster cumulant.
6024 iti=itortyp(itype(i))
6025 itk=itortyp(itype(k))
6026 itk1=itortyp(itype(k+1))
6027 itl=itortyp(itype(l))
6028 itj=itortyp(itype(j))
6029 cd write (2,*) 'itk',itk,' itk1',itk1,' itl',itl,' itj',itj
6030 cd write (2,*) 'i',i,' k',k,' j',j,' l',l
6031 cd if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then
6036 cd & 'EELLO6: Contacts have occurred for peptide groups',i,j,
6038 cd call checkint_turn6(i,jj,kk,eel_turn6_num)
6042 derx_turn(lll,kkk,iii)=0.0d0
6049 eello6_5=eello6_graph4(l,k,j,i,kk,jj,2,.true.)
6051 cd write (2,*) 'eello6_5',eello6_5
6053 call transpose2(AEA(1,1,1),auxmat(1,1))
6054 call matmat2(EUg(1,1,i+1),auxmat(1,1),auxmat(1,1))
6055 ss1=scalar2(Ub2(1,i+2),b1(1,itl))
6056 s1 = (auxmat(1,1)+auxmat(2,2))*ss1
6058 call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1))
6059 call matvec2(AEA(1,1,1),vtemp1(1),vtemp1(1))
6060 s2 = scalar2(b1(1,itk),vtemp1(1))
6062 call transpose2(AEA(1,1,2),atemp(1,1))
6063 call matmat2(atemp(1,1),EUg(1,1,i+4),atemp(1,1))
6064 call matvec2(Ug2(1,1,i+2),dd(1,1,itk1),vtemp2(1))
6065 s8 = -(atemp(1,1)+atemp(2,2))*scalar2(cc(1,1,itl),vtemp2(1))
6067 call matmat2(EUg(1,1,i+3),AEA(1,1,2),auxmat(1,1))
6068 call matvec2(auxmat(1,1),Ub2(1,i+4),vtemp3(1))
6069 s12 = scalar2(Ub2(1,i+2),vtemp3(1))
6071 call transpose2(a_chuj(1,1,kk,i+1),achuj_temp(1,1))
6072 call matmat2(achuj_temp(1,1),EUg(1,1,i+2),gtemp(1,1))
6073 call matmat2(gtemp(1,1),EUg(1,1,i+3),gtemp(1,1))
6074 call matvec2(a_chuj(1,1,jj,i),Ub2(1,i+4),vtemp4(1))
6075 ss13 = scalar2(b1(1,itk),vtemp4(1))
6076 s13 = (gtemp(1,1)+gtemp(2,2))*ss13
6078 c write (2,*) 's1,s2,s8,s12,s13',s1,s2,s8,s12,s13
6084 eel_turn6 = eello6_5 - 0.5d0*(s1+s2+s12+s8+s13)
6086 C Derivatives in gamma(i+2)
6088 call transpose2(AEA(1,1,1),auxmatd(1,1))
6089 call matmat2(EUgder(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
6090 s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1
6091 call transpose2(AEAderg(1,1,2),atempd(1,1))
6092 call matmat2(atempd(1,1),EUg(1,1,i+4),atempd(1,1))
6093 s8d = -(atempd(1,1)+atempd(2,2))*scalar2(cc(1,1,itl),vtemp2(1))
6095 call matmat2(EUg(1,1,i+3),AEAderg(1,1,2),auxmatd(1,1))
6096 call matvec2(auxmatd(1,1),Ub2(1,i+4),vtemp3d(1))
6097 s12d = scalar2(Ub2(1,i+2),vtemp3d(1))
6103 gel_loc_turn6(i)=gel_loc_turn6(i)-0.5d0*ekont*(s1d+s8d+s12d)
6104 C Derivatives in gamma(i+3)
6106 call transpose2(AEA(1,1,1),auxmatd(1,1))
6107 call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
6108 ss1d=scalar2(Ub2der(1,i+2),b1(1,itl))
6109 s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1d
6111 call matvec2(EUgder(1,1,i+2),b1(1,itl),vtemp1d(1))
6112 call matvec2(AEA(1,1,1),vtemp1d(1),vtemp1d(1))
6113 s2d = scalar2(b1(1,itk),vtemp1d(1))
6115 call matvec2(Ug2der(1,1,i+2),dd(1,1,itk1),vtemp2d(1))
6116 s8d = -(atemp(1,1)+atemp(2,2))*scalar2(cc(1,1,itl),vtemp2d(1))
6118 s12d = scalar2(Ub2der(1,i+2),vtemp3(1))
6120 call matmat2(achuj_temp(1,1),EUgder(1,1,i+2),gtempd(1,1))
6121 call matmat2(gtempd(1,1),EUg(1,1,i+3),gtempd(1,1))
6122 s13d = (gtempd(1,1)+gtempd(2,2))*ss13
6130 gel_loc_turn6(i+1)=gel_loc_turn6(i+1)
6131 & -0.5d0*ekont*(s1d+s2d+s8d+s12d+s13d)
6133 gel_loc_turn6(i+1)=gel_loc_turn6(i+1)
6134 & -0.5d0*ekont*(s2d+s12d)
6136 C Derivatives in gamma(i+4)
6137 call matmat2(EUgder(1,1,i+3),AEA(1,1,2),auxmatd(1,1))
6138 call matvec2(auxmatd(1,1),Ub2(1,i+4),vtemp3d(1))
6139 s12d = scalar2(Ub2(1,i+2),vtemp3d(1))
6141 call matmat2(achuj_temp(1,1),EUg(1,1,i+2),gtempd(1,1))
6142 call matmat2(gtempd(1,1),EUgder(1,1,i+3),gtempd(1,1))
6143 s13d = (gtempd(1,1)+gtempd(2,2))*ss13
6151 gel_loc_turn6(i+2)=gel_loc_turn6(i+2)-0.5d0*ekont*(s12d+s13d)
6153 gel_loc_turn6(i+2)=gel_loc_turn6(i+2)-0.5d0*ekont*(s12d)
6155 C Derivatives in gamma(i+5)
6157 call transpose2(AEAderg(1,1,1),auxmatd(1,1))
6158 call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
6159 s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1
6161 call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1d(1))
6162 call matvec2(AEAderg(1,1,1),vtemp1d(1),vtemp1d(1))
6163 s2d = scalar2(b1(1,itk),vtemp1d(1))
6165 call transpose2(AEA(1,1,2),atempd(1,1))
6166 call matmat2(atempd(1,1),EUgder(1,1,i+4),atempd(1,1))
6167 s8d = -(atempd(1,1)+atempd(2,2))*scalar2(cc(1,1,itl),vtemp2(1))
6169 call matvec2(auxmat(1,1),Ub2der(1,i+4),vtemp3d(1))
6170 s12d = scalar2(Ub2(1,i+2),vtemp3d(1))
6172 call matvec2(a_chuj(1,1,jj,i),Ub2der(1,i+4),vtemp4d(1))
6173 ss13d = scalar2(b1(1,itk),vtemp4d(1))
6174 s13d = (gtemp(1,1)+gtemp(2,2))*ss13d
6182 gel_loc_turn6(i+3)=gel_loc_turn6(i+3)
6183 & -0.5d0*ekont*(s1d+s2d+s8d+s12d+s13d)
6185 gel_loc_turn6(i+3)=gel_loc_turn6(i+3)
6186 & -0.5d0*ekont*(s2d+s12d)
6188 C Cartesian derivatives
6193 call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmatd(1,1))
6194 call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
6195 s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1
6197 call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1))
6198 call matvec2(AEAderx(1,1,lll,kkk,iii,1),vtemp1(1),
6200 s2d = scalar2(b1(1,itk),vtemp1d(1))
6202 call transpose2(AEAderx(1,1,lll,kkk,iii,2),atempd(1,1))
6203 call matmat2(atempd(1,1),EUg(1,1,i+4),atempd(1,1))
6204 s8d = -(atempd(1,1)+atempd(2,2))*
6205 & scalar2(cc(1,1,itl),vtemp2(1))
6207 call matmat2(EUg(1,1,i+3),AEAderx(1,1,lll,kkk,iii,2),
6209 call matvec2(auxmatd(1,1),Ub2(1,i+4),vtemp3d(1))
6210 s12d = scalar2(Ub2(1,i+2),vtemp3d(1))
6217 derx_turn(lll,kkk,iii) = derx_turn(lll,kkk,iii)
6220 derx_turn(lll,kkk,iii) = derx_turn(lll,kkk,iii)
6224 derx_turn(lll,kkk,3-iii) = derx_turn(lll,kkk,3-iii)
6225 & - 0.5d0*(s8d+s12d)
6227 derx_turn(lll,kkk,3-iii) = derx_turn(lll,kkk,3-iii)
6236 call transpose2(a_chuj_der(1,1,lll,kkk,kk,i+1),
6238 call matmat2(achuj_tempd(1,1),EUg(1,1,i+2),gtempd(1,1))
6239 call matmat2(gtempd(1,1),EUg(1,1,i+3),gtempd(1,1))
6240 s13d=(gtempd(1,1)+gtempd(2,2))*ss13
6241 derx_turn(lll,kkk,2) = derx_turn(lll,kkk,2)-0.5d0*s13d
6242 call matvec2(a_chuj_der(1,1,lll,kkk,jj,i),Ub2(1,i+4),
6244 ss13d = scalar2(b1(1,itk),vtemp4d(1))
6245 s13d = (gtemp(1,1)+gtemp(2,2))*ss13d
6246 derx_turn(lll,kkk,1) = derx_turn(lll,kkk,1)-0.5d0*s13d
6250 cd write(iout,*) 'eel6_turn6',eel_turn6,' eel_turn6_num',
6251 cd & 16*eel_turn6_num
6253 if (j.lt.nres-1) then
6260 if (l.lt.nres-1) then
6268 ggg1(ll)=eel_turn6*g_contij(ll,1)
6269 ggg2(ll)=eel_turn6*g_contij(ll,2)
6270 ghalf=0.5d0*ggg1(ll)
6272 gcorr6_turn(ll,i)=gcorr6_turn(ll,i)+ghalf
6273 & +ekont*derx_turn(ll,2,1)
6274 gcorr6_turn(ll,i+1)=gcorr6_turn(ll,i+1)+ekont*derx_turn(ll,3,1)
6275 gcorr6_turn(ll,j)=gcorr6_turn(ll,j)+ghalf
6276 & +ekont*derx_turn(ll,4,1)
6277 gcorr6_turn(ll,j1)=gcorr6_turn(ll,j1)+ekont*derx_turn(ll,5,1)
6278 ghalf=0.5d0*ggg2(ll)
6280 gcorr6_turn(ll,k)=gcorr6_turn(ll,k)+ghalf
6281 & +ekont*derx_turn(ll,2,2)
6282 gcorr6_turn(ll,k+1)=gcorr6_turn(ll,k+1)+ekont*derx_turn(ll,3,2)
6283 gcorr6_turn(ll,l)=gcorr6_turn(ll,l)+ghalf
6284 & +ekont*derx_turn(ll,4,2)
6285 gcorr6_turn(ll,l1)=gcorr6_turn(ll,l1)+ekont*derx_turn(ll,5,2)
6290 gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg1(ll)
6295 gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg2(ll)
6301 gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,1)
6306 gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,2)
6310 cd write (2,*) iii,g_corr6_loc(iii)
6313 eello_turn6=ekont*eel_turn6
6314 cd write (2,*) 'ekont',ekont
6315 cd write (2,*) 'eel_turn6',ekont*eel_turn6
6318 crc-------------------------------------------------
6319 SUBROUTINE MATVEC2(A1,V1,V2)
6320 implicit real*8 (a-h,o-z)
6321 include 'DIMENSIONS'
6322 DIMENSION A1(2,2),V1(2),V2(2)
6326 c 3 VI=VI+A1(I,K)*V1(K)
6330 vaux1=a1(1,1)*v1(1)+a1(1,2)*v1(2)
6331 vaux2=a1(2,1)*v1(1)+a1(2,2)*v1(2)
6336 C---------------------------------------
6337 SUBROUTINE MATMAT2(A1,A2,A3)
6338 implicit real*8 (a-h,o-z)
6339 include 'DIMENSIONS'
6340 DIMENSION A1(2,2),A2(2,2),A3(2,2)
6341 c DIMENSION AI3(2,2)
6345 c A3IJ=A3IJ+A1(I,K)*A2(K,J)
6351 ai3_11=a1(1,1)*a2(1,1)+a1(1,2)*a2(2,1)
6352 ai3_12=a1(1,1)*a2(1,2)+a1(1,2)*a2(2,2)
6353 ai3_21=a1(2,1)*a2(1,1)+a1(2,2)*a2(2,1)
6354 ai3_22=a1(2,1)*a2(1,2)+a1(2,2)*a2(2,2)
6362 c-------------------------------------------------------------------------
6363 double precision function scalar2(u,v)
6365 double precision u(2),v(2)
6368 scalar2=u(1)*v(1)+u(2)*v(2)
6372 C-----------------------------------------------------------------------------
6374 subroutine transpose2(a,at)
6376 double precision a(2,2),at(2,2)
6383 c--------------------------------------------------------------------------
6384 subroutine transpose(n,a,at)
6387 double precision a(n,n),at(n,n)
6395 C---------------------------------------------------------------------------
6396 subroutine prodmat3(a1,a2,kk,transp,prod)
6399 double precision a1(2,2),a2(2,2),a2t(2,2),kk(2,2),prod(2,2)
6401 crc double precision auxmat(2,2),prod_(2,2)
6404 crc call transpose2(kk(1,1),auxmat(1,1))
6405 crc call matmat2(a1(1,1),auxmat(1,1),auxmat(1,1))
6406 crc call matmat2(auxmat(1,1),a2(1,1),prod_(1,1))
6408 prod(1,1)=(a1(1,1)*kk(1,1)+a1(1,2)*kk(1,2))*a2(1,1)
6409 & +(a1(1,1)*kk(2,1)+a1(1,2)*kk(2,2))*a2(2,1)
6410 prod(1,2)=(a1(1,1)*kk(1,1)+a1(1,2)*kk(1,2))*a2(1,2)
6411 & +(a1(1,1)*kk(2,1)+a1(1,2)*kk(2,2))*a2(2,2)
6412 prod(2,1)=(a1(2,1)*kk(1,1)+a1(2,2)*kk(1,2))*a2(1,1)
6413 & +(a1(2,1)*kk(2,1)+a1(2,2)*kk(2,2))*a2(2,1)
6414 prod(2,2)=(a1(2,1)*kk(1,1)+a1(2,2)*kk(1,2))*a2(1,2)
6415 & +(a1(2,1)*kk(2,1)+a1(2,2)*kk(2,2))*a2(2,2)
6418 crc call matmat2(a1(1,1),kk(1,1),auxmat(1,1))
6419 crc call matmat2(auxmat(1,1),a2(1,1),prod_(1,1))
6421 prod(1,1)=(a1(1,1)*kk(1,1)+a1(1,2)*kk(2,1))*a2(1,1)
6422 & +(a1(1,1)*kk(1,2)+a1(1,2)*kk(2,2))*a2(2,1)
6423 prod(1,2)=(a1(1,1)*kk(1,1)+a1(1,2)*kk(2,1))*a2(1,2)
6424 & +(a1(1,1)*kk(1,2)+a1(1,2)*kk(2,2))*a2(2,2)
6425 prod(2,1)=(a1(2,1)*kk(1,1)+a1(2,2)*kk(2,1))*a2(1,1)
6426 & +(a1(2,1)*kk(1,2)+a1(2,2)*kk(2,2))*a2(2,1)
6427 prod(2,2)=(a1(2,1)*kk(1,1)+a1(2,2)*kk(2,1))*a2(1,2)
6428 & +(a1(2,1)*kk(1,2)+a1(2,2)*kk(2,2))*a2(2,2)
6431 c call transpose2(a2(1,1),a2t(1,1))
6434 crc print *,((prod_(i,j),i=1,2),j=1,2)
6435 crc print *,((prod(i,j),i=1,2),j=1,2)
6439 C-----------------------------------------------------------------------------
6440 double precision function scalar(u,v)
6442 double precision u(3),v(3)