X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Fenergy_p_new.F.org;fp=source%2Fwham%2Fsrc-NEWSC%2Fenergy_p_new.F.org;h=8f99a1612a2bfceb4ac937d443fd4f7ff3f2894e;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/wham/src-NEWSC/energy_p_new.F.org b/source/wham/src-NEWSC/energy_p_new.F.org new file mode 100755 index 0000000..8f99a16 --- /dev/null +++ b/source/wham/src-NEWSC/energy_p_new.F.org @@ -0,0 +1,6452 @@ + subroutine etotal(energia) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + + external proc_proc +#ifdef WINPGI +cMS$ATTRIBUTES C :: proc_proc +#endif + + include 'COMMON.IOUNITS' + double precision energia(0:max_ene),energia1(0:max_ene+1) +#ifdef MPL + include 'COMMON.INFO' + external d_vadd + integer ready +#endif + include 'COMMON.FFIELD' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' +cd write(iout, '(a,i2)')'Calling etotal ipot=',ipot +cd print *,'nnt=',nnt,' nct=',nct +C +C Compute the side-chain and electrostatic interaction energy +C + goto (101,102,103,104,105) ipot +C Lennard-Jones potential. + 101 call elj(evdw) +cd print '(a)','Exit ELJ' + goto 106 +C Lennard-Jones-Kihara potential (shifted). + 102 call eljk(evdw) + goto 106 +C Berne-Pechukas potential (dilated LJ, angular dependence). + 103 call ebp(evdw) + goto 106 +C Gay-Berne potential (shifted LJ, angular dependence). + 104 call egb(evdw) + goto 106 +C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence). + 105 call egbv(evdw) +C +C Calculate electrostatic (H-bonding) energy of the main chain. +C + 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) +C +C Calculate excluded-volume interaction energy between peptide groups +C and side chains. +C + call escp(evdw2,evdw2_14) +C +C Calculate the disulfide-bridge and other energy and the contributions +C from other distance constraints. +cd print *,'Calling EHPB' + call edis(ehpb) +cd print *,'EHPB exitted succesfully.' +C +C Calculate the virtual-bond-angle energy. +C + call ebend(ebe) +cd print *,'Bend energy finished.' +C +C Calculate the SC local energy. +C + call esc(escloc) +cd print *,'SCLOC energy finished.' +C +C Calculate the virtual-bond torsional energy. +C +cd print *,'nterm=',nterm + call etor(etors,edihcnstr) +C +C 6/23/01 Calculate double-torsional energy +C + call etor_d(etors_d) +C +C 12/1/95 Multi-body terms +C + n_corr=0 + n_corr1=0 + if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 + & .or. wturn6.gt.0.0d0) then +c print *,"calling multibody_eello" + call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1) +c write (*,*) 'n_corr=',n_corr,' n_corr1=',n_corr1 +c print *,ecorr,ecorr5,ecorr6,eturn6 + endif + if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then + call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1) + endif +C call multibody(ecorr) +C +C Sum the energies +C +C scale large componenets +#ifdef SCALE + ecorr5_scal=1000.0 + eel_loc_scal=100.0 + eello_turn3_scal=100.0 + eello_turn4_scal=100.0 + eturn6_scal=1000.0 + ecorr6_scal=1000.0 +#else + ecorr5_scal=1.0 + eel_loc_scal=1.0 + eello_turn3_scal=1.0 + eello_turn4_scal=1.0 + eturn6_scal=1.0 + ecorr6_scal=1.0 +#endif + + ecorr5=ecorr5/ecorr5_scal + eel_loc=eel_loc/eel_loc_scal + eello_turn3=eello_turn3/eello_turn3_scal + eello_turn4=eello_turn4/eello_turn4_scal + eturn6=eturn6/eturn6_scal + ecorr6=ecorr6/ecorr6_scal +#ifdef MPL + if (fgprocs.gt.1) then +cd call enerprint(evdw,evdw1,evdw2,ees,ebe,escloc,etors,ehpb, +cd & edihcnstr,ecorr,eel_loc,eello_turn4,etot) + energia(1)=evdw + energia(2)=evdw2 + energia(3)=ees + energia(4)=evdw1 + energia(5)=ecorr + energia(6)=etors + energia(7)=ebe + energia(8)=escloc + energia(9)=ehpb + energia(10)=edihcnstr + energia(11)=eel_loc + energia(12)=ecorr5 + energia(13)=ecorr6 + energia(14)=eello_turn3 + energia(15)=eello_turn4 + energia(16)=eturn6 + energia(17)=etors_d + msglen=80 + do i=1,15 + energia1(i)=energia(i) + enddo +cd write (iout,*) 'BossID=',BossID,' MyGroup=',MyGroup +cd write (*,*) 'BossID=',BossID,' MyGroup=',MyGroup +cd write (*,*) 'Processor',MyID,' calls MP_REDUCE in ENERGY', +cd & ' BossID=',BossID,' MyGroup=',MyGroup + call mp_reduce(energia1(1),energia(1),msglen,BossID,d_vadd, + & fgGroupID) +cd write (iout,*) 'Processor',MyID,' Reduce finished' + evdw=energia(1) + evdw2=energia(2) + ees=energia(3) + evdw1=energia(4) + ecorr=energia(5) + etors=energia(6) + ebe=energia(7) + escloc=energia(8) + ehpb=energia(9) + edihcnstr=energia(10) + eel_loc=energia(11) + ecorr5=energia(12) + ecorr6=energia(13) + eello_turn3=energia(14) + eello_turn4=energia(15) + eturn6=energia(16) + etors_d=energia(17) + endif +c if (MyID.eq.BossID) then +#endif + etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) + & +wang*ebe+wtor*etors+wscloc*escloc + & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 + & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 + & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d + energia(0)=etot + energia(1)=evdw + energia(2)=evdw2 + energia(3)=ees+evdw1 + energia(4)=ecorr + energia(5)=ecorr5 + energia(6)=ecorr6 + energia(7)=eel_loc + energia(8)=eello_turn3 + energia(9)=eello_turn4 + energia(10)=eturn6 + energia(11)=ebe + energia(12)=escloc + energia(13)=etors + energia(14)=etors_d + energia(15)=ehpb + energia(16)=edihcnstr + energia(17)=evdw2_14 +c detecting NaNQ + i=0 +#ifdef WINPGI + idumm=proc_proc(etot,i) +#else + call proc_proc(etot,i) +#endif + if(i.eq.1)energia(0)=1.0d+99 +#ifdef MPL +c endif +#endif + if (calc_grad) then +C +C Sum up the components of the Cartesian gradient. +C + do i=1,nct + do j=1,3 + gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+ + & welec*gelc(j,i)+wstrain*ghpbc(j,i)+ + & wcorr*gradcorr(j,i)+ + & wel_loc*gel_loc(j,i)/eel_loc_scal+ + & wturn3*gcorr3_turn(j,i)/eello_turn3_scal+ + & wturn4*gcorr4_turn(j,i)/eello_turn4_scal+ + & wcorr5*gradcorr5(j,i)/ecorr5_scal+ + & wcorr6*gradcorr6(j,i)/ecorr6_scal+ + & wturn6*gcorr6_turn(j,i)/eturn6_scal + gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ + & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i) + enddo +cd print '(i3,9(1pe12.4))',i,(gvdwc(k,i),k=1,3),(gelc(k,i),k=1,3), +cd & (gradc(k,i),k=1,3) + enddo + + + do i=1,nres-3 +cd write (iout,*) i,g_corr5_loc(i) + gloc(i,icg)=gloc(i,icg)+wcorr*gcorr_loc(i) + & +wcorr5*g_corr5_loc(i)/ecorr5_scal + & +wcorr6*g_corr6_loc(i)/ecorr6_scal + & +wturn4*gel_loc_turn4(i)/eello_turn4_scal + & +wturn3*gel_loc_turn3(i)/eello_turn3_scal + & +wturn6*gel_loc_turn6(i)/eturn6_scal + & +wel_loc*gel_loc_loc(i)/eel_loc_scal + enddo + endif +cd print*,evdw,wsc,evdw2,wscp,ees+evdw1,welec,ebe,wang, +cd & escloc,wscloc,etors,wtor,ehpb,wstrain,nss,ebr,etot +cd call enerprint(energia(0)) +cd call intout +cd stop + return + end +C------------------------------------------------------------------------ + subroutine enerprint(energia) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.IOUNITS' + include 'COMMON.FFIELD' + include 'COMMON.SBRIDGE' + double precision energia(0:max_ene) + etot=energia(0) + evdw=energia(1) + evdw2=energia(2) + ees=energia(3) + ecorr=energia(4) + ecorr5=energia(5) + ecorr6=energia(6) + eel_loc=energia(7) + eello_turn3=energia(8) + eello_turn4=energia(9) + eello_turn6=energia(10) + ebe=energia(11) + escloc=energia(12) + etors=energia(13) + etors_d=energia(14) + ehpb=energia(15) + edihcnstr=energia(16) + write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,ebe,wang, + & escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain, + & ecorr,wcorr, + & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, + & eello_turn4,wturn4,eello_turn6,wturn6,edihcnstr,ebr*nss,etot + 10 format (/'Virtual-chain energies:'// + & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ + & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ + & 'EES= ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p)'/ + & 'EBE= ',1pE16.6,' WEIGHT=',1pD16.6,' (bending)'/ + & 'ESC= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/ + & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/ + & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/ + & 'EHBP= ',1pE16.6,' WEIGHT=',1pD16.6, + & ' (SS bridges & dist. cnstr.)'/ + & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ + & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ + & 'ECORR6=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ + & 'EELLO= ',1pE16.6,' WEIGHT=',1pD16.6,' (electrostatic-local)'/ + & 'ETURN3=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 3rd order)'/ + & 'ETURN4=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 4th order)'/ + & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ + & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ + & 'ETOT= ',1pE16.6,' (total)') + return + end +C----------------------------------------------------------------------- + subroutine elj(evdw) +C +C This subroutine calculates the interaction energy of nonbonded side chains +C assuming the LJ potential of interaction. +C + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + parameter (accur=1.0d-10) + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.TORSION' + include 'COMMON.ENEPS' + include 'COMMON.SBRIDGE' + include 'COMMON.NAMES' + include 'COMMON.IOUNITS' + include 'COMMON.CONTACTS' + dimension gg(3) + integer icant + external icant +cd print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon + do i=1,210 + do j=1,2 + eneps_temp(j,i)=0.0d0 + enddo + enddo + evdw=0.0D0 + do i=iatsc_s,iatsc_e + itypi=itype(i) + itypi1=itype(i+1) + xi=c(1,nres+i) + yi=c(2,nres+i) + zi=c(3,nres+i) +C Change 12/1/95 + num_conti=0 +C +C Calculate SC interaction energy. +C + do iint=1,nint_gr(i) +cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), +cd & 'iend=',iend(i,iint) + do j=istart(i,iint),iend(i,iint) + itypj=itype(j) + xj=c(1,nres+j)-xi + yj=c(2,nres+j)-yi + zj=c(3,nres+j)-zi +C Change 12/1/95 to calculate four-body interactions + rij=xj*xj+yj*yj+zj*zj + rrij=1.0D0/rij +c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj + eps0ij=eps(itypi,itypj) + fac=rrij**expon2 + e1=fac*fac*aa(itypi,itypj) + e2=fac*bb(itypi,itypj) + evdwij=e1+e2 + ij=icant(itypi,itypj) + eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij) + eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij +cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) +cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) +cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') +cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), +cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm, +cd & (c(k,i),k=1,3),(c(k,j),k=1,3) + evdw=evdw+evdwij + if (calc_grad) then +C +C Calculate the components of the gradient in DC and X +C + fac=-rrij*(e1+evdwij) + gg(1)=xj*fac + gg(2)=yj*fac + gg(3)=zj*fac + do k=1,3 + gvdwx(k,i)=gvdwx(k,i)-gg(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k) + enddo + do k=i,j-1 + do l=1,3 + gvdwc(l,k)=gvdwc(l,k)+gg(l) + enddo + enddo + endif +C +C 12/1/95, revised on 5/20/97 +C +C Calculate the contact function. The ith column of the array JCONT will +C contain the numbers of atoms that make contacts with the atom I (of numbers +C greater than I). The arrays FACONT and GACONT will contain the values of +C the contact function and its derivative. +C +C Uncomment next line, if the correlation interactions include EVDW explicitly. +c if (j.gt.i+1 .and. evdwij.le.0.0D0) then +C Uncomment next line, if the correlation interactions are contact function only + if (j.gt.i+1.and. eps0ij.gt.0.0D0) then + rij=dsqrt(rij) + sigij=sigma(itypi,itypj) + r0ij=rs0(itypi,itypj) +C +C Check whether the SC's are not too far to make a contact. +C + rcut=1.5d0*r0ij + call gcont(rij,rcut,1.0d0,0.2d0*rcut,fcont,fprimcont) +C Add a new contact, if the SC's are close enough, but not too close (ri' + do k=1,3 + ggg(k)=-ggg(k) +C Uncomment following line for SC-p interactions +c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k) + enddo + endif + do k=1,3 + gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k) + enddo + kstart=min0(i+1,j) + kend=max0(i-1,j-1) +cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend +cd write (iout,*) ggg(1),ggg(2),ggg(3) + do k=kstart,kend + do l=1,3 + gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l) + enddo + enddo + endif + enddo + enddo ! iint + 1225 continue + enddo ! i + do i=1,nct + do j=1,3 + gvdwc_scp(j,i)=expon*gvdwc_scp(j,i) + gradx_scp(j,i)=expon*gradx_scp(j,i) + enddo + enddo +C****************************************************************************** +C +C N O T E !!! +C +C To save time the factor EXPON has been extracted from ALL components +C of GVDWC and GRADX. Remember to multiply them by this factor before further +C use! +C +C****************************************************************************** + return + end +C-------------------------------------------------------------------------- + subroutine edis(ehpb) +C +C Evaluate bridge-strain energy and its gradient in virtual-bond and SC vectors. +C + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.VAR' + dimension ggg(3) + ehpb=0.0D0 +cd print *,'edis: nhpb=',nhpb,' fbr=',fbr +cd print *,'link_start=',link_start,' link_end=',link_end + if (link_end.eq.0) return + do i=link_start,link_end +C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a +C CA-CA distance used in regularization of structure. + ii=ihpb(i) + jj=jhpb(i) +C iii and jjj point to the residues for which the distance is assigned. + if (ii.gt.nres) then + iii=ii-nres + jjj=jj-nres + else + iii=ii + jjj=jj + endif +C Calculate the distance between the two points and its difference from the +C target distance. + dd=dist(ii,jj) + rdis=dd-dhpb(i) +C Get the force constant corresponding to this distance. + waga=forcon(i) +C Calculate the contribution to energy. + ehpb=ehpb+waga*rdis*rdis +C +C Evaluate gradient. +C + fac=waga*rdis/dd +cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, +cd & ' waga=',waga,' fac=',fac + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo +cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) +C If this is a SC-SC distace, we need to calculate the contributions to the +C Cartesian gradient in the SC vectors (ghpbx). + if (iii.lt.ii) then + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + endif + do j=iii,jjj-1 + do k=1,3 + ghpbc(k,j)=ghpbc(k,j)+ggg(k) + enddo + enddo + enddo + ehpb=0.5D0*ehpb + return + end +C-------------------------------------------------------------------------- + subroutine ebend(etheta) +C +C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral +C angles gamma and its derivatives in consecutive thetas and gammas. +C + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.LOCAL' + include 'COMMON.GEO' + include 'COMMON.INTERACT' + include 'COMMON.DERIV' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + include 'COMMON.FFIELD' + common /calcthet/ term1,term2,termm,diffak,ratak, + & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, + & delthe0,sig0inv,sigtc,sigsqtc,delthec,it + double precision y(2),z(2) + delta=0.02d0*pi + time11=dexp(-2*time) + time12=1.0d0 + etheta=0.0D0 +c write (iout,*) "nres",nres +c write (*,'(a,i2)') 'EBEND ICG=',icg +c write (iout,*) ithet_start,ithet_end + do i=ithet_start,ithet_end +C Zero the energy function and its derivative at 0 or pi. + call splinthet(theta(i),0.5d0*delta,ss,ssd) + it=itype(i-1) + if (i.gt.ithet_start .and. + & (itel(i-1).eq.0 .or. itel(i-2).eq.0)) goto 1215 + if (i.gt.3 .and. (i.le.4 .or. itel(i-3).ne.0)) then + phii=phi(i) + y(1)=dcos(phii) + y(2)=dsin(phii) + else + y(1)=0.0D0 + y(2)=0.0D0 + endif + if (i.lt.nres .and. itel(i).ne.0) then + phii1=phi(i+1) + z(1)=dcos(phii1) + z(2)=dsin(phii1) + else + z(1)=0.0D0 + z(2)=0.0D0 + endif +C Calculate the "mean" value of theta from the part of the distribution +C dependent on the adjacent virtual-bond-valence angles (gamma1 & gamma2). +C In following comments this theta will be referred to as t_c. + thet_pred_mean=0.0d0 + do k=1,2 + athetk=athet(k,it) + bthetk=bthet(k,it) + thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k) + enddo +c write (iout,*) "thet_pred_mean",thet_pred_mean + dthett=thet_pred_mean*ssd + thet_pred_mean=thet_pred_mean*ss+a0thet(it) +c write (iout,*) "thet_pred_mean",thet_pred_mean +C Derivatives of the "mean" values in gamma1 and gamma2. + dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss + dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss + if (theta(i).gt.pi-delta) then + call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0, + & E_tc0) + call mixder(pi-delta,thet_pred_mean,theta0(it),fprim_tc0) + call theteng(pi,thet_pred_mean,theta0(it),f1,fprim1,E_tc1) + call spline1(theta(i),pi-delta,delta,f0,f1,fprim0,ethetai, + & E_theta) + call spline2(theta(i),pi-delta,delta,E_tc0,E_tc1,fprim_tc0, + & E_tc) + else if (theta(i).lt.delta) then + call theteng(delta,thet_pred_mean,theta0(it),f0,fprim0,E_tc0) + call theteng(0.0d0,thet_pred_mean,theta0(it),f1,fprim1,E_tc1) + call spline1(theta(i),delta,-delta,f0,f1,fprim0,ethetai, + & E_theta) + call mixder(delta,thet_pred_mean,theta0(it),fprim_tc0) + call spline2(theta(i),delta,-delta,E_tc0,E_tc1,fprim_tc0, + & E_tc) + else + call theteng(theta(i),thet_pred_mean,theta0(it),ethetai, + & E_theta,E_tc) + endif + etheta=etheta+ethetai +c write (iout,'(2i3,3f8.3,f10.5)') i,it,rad2deg*theta(i), +c & rad2deg*phii,rad2deg*phii1,ethetai + if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 + if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 + gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) + 1215 continue + enddo +C Ufff.... We've done all this!!! + return + end +C--------------------------------------------------------------------------- + subroutine theteng(thetai,thet_pred_mean,theta0i,ethetai,E_theta, + & E_tc) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.LOCAL' + include 'COMMON.IOUNITS' + common /calcthet/ term1,term2,termm,diffak,ratak, + & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, + & delthe0,sig0inv,sigtc,sigsqtc,delthec,it +C Calculate the contributions to both Gaussian lobes. +C 6/6/97 - Deform the Gaussians using the factor of 1/(1+time) +C The "polynomial part" of the "standard deviation" of this part of +C the distribution. + sig=polthet(3,it) + do j=2,0,-1 + sig=sig*thet_pred_mean+polthet(j,it) + enddo +C Derivative of the "interior part" of the "standard deviation of the" +C gamma-dependent Gaussian lobe in t_c. + sigtc=3*polthet(3,it) + do j=2,1,-1 + sigtc=sigtc*thet_pred_mean+j*polthet(j,it) + enddo + sigtc=sig*sigtc +C Set the parameters of both Gaussian lobes of the distribution. +C "Standard deviation" of the gamma-dependent Gaussian lobe (sigtc) + fac=sig*sig+sigc0(it) + sigcsq=fac+fac + sigc=1.0D0/sigcsq +C Following variable (sigsqtc) is -(1/2)d[sigma(t_c)**(-2))]/dt_c + sigsqtc=-4.0D0*sigcsq*sigtc +c print *,i,sig,sigtc,sigsqtc +C Following variable (sigtc) is d[sigma(t_c)]/dt_c + sigtc=-sigtc/(fac*fac) +C Following variable is sigma(t_c)**(-2) + sigcsq=sigcsq*sigcsq + sig0i=sig0(it) + sig0inv=1.0D0/sig0i**2 + delthec=thetai-thet_pred_mean + delthe0=thetai-theta0i + term1=-0.5D0*sigcsq*delthec*delthec + term2=-0.5D0*sig0inv*delthe0*delthe0 +C Following fuzzy logic is to avoid underflows in dexp and subsequent INFs and +C NaNs in taking the logarithm. We extract the largest exponent which is added +C to the energy (this being the log of the distribution) at the end of energy +C term evaluation for this virtual-bond angle. + if (term1.gt.term2) then + termm=term1 + term2=dexp(term2-termm) + term1=1.0d0 + else + termm=term2 + term1=dexp(term1-termm) + term2=1.0d0 + endif +C The ratio between the gamma-independent and gamma-dependent lobes of +C the distribution is a Gaussian function of thet_pred_mean too. + diffak=gthet(2,it)-thet_pred_mean + ratak=diffak/gthet(3,it)**2 + ak=dexp(gthet(1,it)-0.5D0*diffak*ratak) +C Let's differentiate it in thet_pred_mean NOW. + aktc=ak*ratak +C Now put together the distribution terms to make complete distribution. + termexp=term1+ak*term2 + termpre=sigc+ak*sig0i +C Contribution of the bending energy from this theta is just the -log of +C the sum of the contributions from the two lobes and the pre-exponential +C factor. Simple enough, isn't it? + ethetai=(-dlog(termexp)-termm+dlog(termpre)) +C NOW the derivatives!!! +C 6/6/97 Take into account the deformation. + E_theta=(delthec*sigcsq*term1 + & +ak*delthe0*sig0inv*term2)/termexp + E_tc=((sigtc+aktc*sig0i)/termpre + & -((delthec*sigcsq+delthec*delthec*sigsqtc)*term1+ + & aktc*term2)/termexp) + return + end +c----------------------------------------------------------------------------- + subroutine mixder(thetai,thet_pred_mean,theta0i,E_tc_t) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.LOCAL' + include 'COMMON.IOUNITS' + common /calcthet/ term1,term2,termm,diffak,ratak, + & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, + & delthe0,sig0inv,sigtc,sigsqtc,delthec,it + delthec=thetai-thet_pred_mean + delthe0=thetai-theta0i +C "Thank you" to MAPLE (probably spared one day of hand-differentiation). + t3 = thetai-thet_pred_mean + t6 = t3**2 + t9 = term1 + t12 = t3*sigcsq + t14 = t12+t6*sigsqtc + t16 = 1.0d0 + t21 = thetai-theta0i + t23 = t21**2 + t26 = term2 + t27 = t21*t26 + t32 = termexp + t40 = t32**2 + E_tc_t = -((sigcsq+2.D0*t3*sigsqtc)*t9-t14*sigcsq*t3*t16*t9 + & -aktc*sig0inv*t27)/t32+(t14*t9+aktc*t26)/t40 + & *(-t12*t9-ak*sig0inv*t27) + return + end +c----------------------------------------------------------------------------- + subroutine esc(escloc) +C Calculate the local energy of a side chain and its derivatives in the +C corresponding virtual-bond valence angles THETA and the spherical angles +C ALPHA and OMEGA. + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.INTERACT' + include 'COMMON.DERIV' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + include 'COMMON.FFIELD' + double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3), + & ddersc0(3),ddummy(3),xtemp(3),temp(3) + common /sccalc/ time11,time12,time112,theti,it,nlobit + delta=0.02d0*pi + escloc=0.0D0 +c write (iout,'(a)') 'ESC' + do i=loc_start,loc_end + it=itype(i) + if (it.eq.10) goto 1 + nlobit=nlob(it) +c print *,'i=',i,' it=',it,' nlobit=',nlobit +c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad + theti=theta(i+1)-pipol + x(1)=dtan(theti) + x(2)=alph(i) + x(3)=omeg(i) + + if (x(2).gt.pi-delta) then + xtemp(1)=x(1) + xtemp(2)=pi-delta + xtemp(3)=x(3) + call enesc(xtemp,escloci0,dersc0,ddersc0,.true.) + xtemp(2)=pi + call enesc(xtemp,escloci1,dersc1,ddummy,.false.) + call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2), + & escloci,dersc(2)) + call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1), + & ddersc0(1),dersc(1)) + call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3), + & ddersc0(3),dersc(3)) + xtemp(2)=pi-delta + call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.) + xtemp(2)=pi + call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.) + call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1, + & dersc0(2),esclocbi,dersc02) + call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1), + & dersc12,dersc01) + call splinthet(x(2),0.5d0*delta,ss,ssd) + dersc0(1)=dersc01 + dersc0(2)=dersc02 + dersc0(3)=0.0d0 + do k=1,3 + dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k) + enddo + dersc(2)=dersc(2)+ssd*(escloci-esclocbi) +c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, +c & esclocbi,ss,ssd + escloci=ss*escloci+(1.0d0-ss)*esclocbi +c escloci=esclocbi +c write (iout,*) escloci + else if (x(2).lt.delta) then + xtemp(1)=x(1) + xtemp(2)=delta + xtemp(3)=x(3) + call enesc(xtemp,escloci0,dersc0,ddersc0,.true.) + xtemp(2)=0.0d0 + call enesc(xtemp,escloci1,dersc1,ddummy,.false.) + call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2), + & escloci,dersc(2)) + call spline2(x(2),delta,-delta,dersc0(1),dersc1(1), + & ddersc0(1),dersc(1)) + call spline2(x(2),delta,-delta,dersc0(3),dersc1(3), + & ddersc0(3),dersc(3)) + xtemp(2)=delta + call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.) + xtemp(2)=0.0d0 + call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.) + call spline1(x(2),delta,-delta,esclocbi0,esclocbi1, + & dersc0(2),esclocbi,dersc02) + call spline2(x(2),delta,-delta,dersc0(1),dersc1(1), + & dersc12,dersc01) + dersc0(1)=dersc01 + dersc0(2)=dersc02 + dersc0(3)=0.0d0 + call splinthet(x(2),0.5d0*delta,ss,ssd) + do k=1,3 + dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k) + enddo + dersc(2)=dersc(2)+ssd*(escloci-esclocbi) +c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, +c & esclocbi,ss,ssd + escloci=ss*escloci+(1.0d0-ss)*esclocbi +c write (iout,*) escloci + else + call enesc(x,escloci,dersc,ddummy,.false.) + endif + + escloc=escloc+escloci +c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc + + gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ + & wscloc*dersc(1) + gloc(ialph(i,1),icg)=wscloc*dersc(2) + gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3) + 1 continue + enddo + return + end +C--------------------------------------------------------------------------- + subroutine enesc(x,escloci,dersc,ddersc,mixed) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.IOUNITS' + common /sccalc/ time11,time12,time112,theti,it,nlobit + double precision x(3),z(3),Ax(3,maxlob,-1:1),dersc(3),ddersc(3) + double precision contr(maxlob,-1:1) + logical mixed +c write (iout,*) 'it=',it,' nlobit=',nlobit + escloc_i=0.0D0 + do j=1,3 + dersc(j)=0.0D0 + if (mixed) ddersc(j)=0.0d0 + enddo + x3=x(3) + +C Because of periodicity of the dependence of the SC energy in omega we have +C to add up the contributions from x(3)-2*pi, x(3), and x(3+2*pi). +C To avoid underflows, first compute & store the exponents. + + do iii=-1,1 + + x(3)=x3+iii*dwapi + + do j=1,nlobit + do k=1,3 + z(k)=x(k)-censc(k,j,it) + enddo + do k=1,3 + Axk=0.0D0 + do l=1,3 + Axk=Axk+gaussc(l,k,j,it)*z(l) + enddo + Ax(k,j,iii)=Axk + enddo + expfac=0.0D0 + do k=1,3 + expfac=expfac+Ax(k,j,iii)*z(k) + enddo + contr(j,iii)=expfac + enddo ! j + + enddo ! iii + + x(3)=x3 +C As in the case of ebend, we want to avoid underflows in exponentiation and +C subsequent NaNs and INFs in energy calculation. +C Find the largest exponent + emin=contr(1,-1) + do iii=-1,1 + do j=1,nlobit + if (emin.gt.contr(j,iii)) emin=contr(j,iii) + enddo + enddo + emin=0.5D0*emin +cd print *,'it=',it,' emin=',emin + +C Compute the contribution to SC energy and derivatives + do iii=-1,1 + + do j=1,nlobit + expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin) +cd print *,'j=',j,' expfac=',expfac + escloc_i=escloc_i+expfac + do k=1,3 + dersc(k)=dersc(k)+Ax(k,j,iii)*expfac + enddo + if (mixed) then + do k=1,3,2 + ddersc(k)=ddersc(k)+(-Ax(2,j,iii)*Ax(k,j,iii) + & +gaussc(k,2,j,it))*expfac + enddo + endif + enddo + + enddo ! iii + + dersc(1)=dersc(1)/cos(theti)**2 + ddersc(1)=ddersc(1)/cos(theti)**2 + ddersc(3)=ddersc(3) + + escloci=-(dlog(escloc_i)-emin) + do j=1,3 + dersc(j)=dersc(j)/escloc_i + enddo + if (mixed) then + do j=1,3,2 + ddersc(j)=(ddersc(j)/escloc_i+dersc(2)*dersc(j)) + enddo + endif + return + end +C------------------------------------------------------------------------------ + subroutine enesc_bound(x,escloci,dersc,dersc12,mixed) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.IOUNITS' + common /sccalc/ time11,time12,time112,theti,it,nlobit + double precision x(3),z(3),Ax(3,maxlob),dersc(3) + double precision contr(maxlob) + logical mixed + + escloc_i=0.0D0 + + do j=1,3 + dersc(j)=0.0D0 + enddo + + do j=1,nlobit + do k=1,2 + z(k)=x(k)-censc(k,j,it) + enddo + z(3)=dwapi + do k=1,3 + Axk=0.0D0 + do l=1,3 + Axk=Axk+gaussc(l,k,j,it)*z(l) + enddo + Ax(k,j)=Axk + enddo + expfac=0.0D0 + do k=1,3 + expfac=expfac+Ax(k,j)*z(k) + enddo + contr(j)=expfac + enddo ! j + +C As in the case of ebend, we want to avoid underflows in exponentiation and +C subsequent NaNs and INFs in energy calculation. +C Find the largest exponent + emin=contr(1) + do j=1,nlobit + if (emin.gt.contr(j)) emin=contr(j) + enddo + emin=0.5D0*emin + +C Compute the contribution to SC energy and derivatives + + dersc12=0.0d0 + do j=1,nlobit + expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin) + escloc_i=escloc_i+expfac + do k=1,2 + dersc(k)=dersc(k)+Ax(k,j)*expfac + enddo + if (mixed) dersc12=dersc12+(-Ax(2,j)*Ax(1,j) + & +gaussc(1,2,j,it))*expfac + dersc(3)=0.0d0 + enddo + + dersc(1)=dersc(1)/cos(theti)**2 + dersc12=dersc12/cos(theti)**2 + escloci=-(dlog(escloc_i)-emin) + do j=1,2 + dersc(j)=dersc(j)/escloc_i + enddo + if (mixed) dersc12=(dersc12/escloc_i+dersc(2)*dersc(1)) + return + end +c------------------------------------------------------------------------------ + subroutine gcont(rij,r0ij,eps0ij,delta,fcont,fprimcont) +C +C This procedure calculates two-body contact function g(rij) and its derivative: +C +C eps0ij ! x < -1 +C g(rij) = esp0ij*(-0.9375*x+0.625*x**3-0.1875*x**5) ! -1 =< x =< 1 +C 0 ! x > 1 +C +C where x=(rij-r0ij)/delta +C +C rij - interbody distance, r0ij - contact distance, eps0ij - contact energy +C + implicit none + double precision rij,r0ij,eps0ij,fcont,fprimcont + double precision x,x2,x4,delta +c delta=0.02D0*r0ij +c delta=0.2D0*r0ij + x=(rij-r0ij)/delta + if (x.lt.-1.0D0) then + fcont=eps0ij + fprimcont=0.0D0 + else if (x.le.1.0D0) then + x2=x*x + x4=x2*x2 + fcont=eps0ij*(x*(-0.9375D0+0.6250D0*x2-0.1875D0*x4)+0.5D0) + fprimcont=eps0ij * (-0.9375D0+1.8750D0*x2-0.9375D0*x4)/delta + else + fcont=0.0D0 + fprimcont=0.0D0 + endif + return + end +c------------------------------------------------------------------------------ + subroutine splinthet(theti,delta,ss,ssder) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.VAR' + include 'COMMON.GEO' + thetup=pi-delta + thetlow=delta + if (theti.gt.pipol) then + call gcont(theti,thetup,1.0d0,delta,ss,ssder) + else + call gcont(-theti,-thetlow,1.0d0,delta,ss,ssder) + ssder=-ssder + endif + return + end +c------------------------------------------------------------------------------ + subroutine spline1(x,x0,delta,f0,f1,fprim0,f,fprim) + implicit none + double precision x,x0,delta,f0,f1,fprim0,f,fprim + double precision ksi,ksi2,ksi3,a1,a2,a3 + a1=fprim0*delta/(f1-f0) + a2=3.0d0-2.0d0*a1 + a3=a1-2.0d0 + ksi=(x-x0)/delta + ksi2=ksi*ksi + ksi3=ksi2*ksi + f=f0+(f1-f0)*ksi*(a1+ksi*(a2+a3*ksi)) + fprim=(f1-f0)/delta*(a1+ksi*(2*a2+3*ksi*a3)) + return + end +c------------------------------------------------------------------------------ + subroutine spline2(x,x0,delta,f0x,f1x,fprim0x,fx) + implicit none + double precision x,x0,delta,f0x,f1x,fprim0x,fx + double precision ksi,ksi2,ksi3,a1,a2,a3 + ksi=(x-x0)/delta + ksi2=ksi*ksi + ksi3=ksi2*ksi + a1=fprim0x*delta + a2=3*(f1x-f0x)-2*fprim0x*delta + a3=fprim0x*delta-2*(f1x-f0x) + fx=f0x+a1*ksi+a2*ksi2+a3*ksi3 + return + end +C----------------------------------------------------------------------------- +#ifdef CRYST_TOR +C----------------------------------------------------------------------------- + subroutine etor(etors,edihcnstr) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.VAR' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.TORSION' + include 'COMMON.INTERACT' + include 'COMMON.DERIV' + include 'COMMON.CHAIN' + include 'COMMON.NAMES' + include 'COMMON.IOUNITS' + include 'COMMON.FFIELD' + include 'COMMON.TORCNSTR' + logical lprn +C Set lprn=.true. for debugging + lprn=.false. +c lprn=.true. + etors=0.0D0 + do i=iphi_start,iphi_end + itori=itortyp(itype(i-2)) + itori1=itortyp(itype(i-1)) + phii=phi(i) + gloci=0.0D0 +C Proline-Proline pair is a special case... + if (itori.eq.3 .and. itori1.eq.3) then + if (phii.gt.-dwapi3) then + cosphi=dcos(3*phii) + fac=1.0D0/(1.0D0-cosphi) + etorsi=v1(1,3,3)*fac + etorsi=etorsi+etorsi + etors=etors+etorsi-v1(1,3,3) + gloci=gloci-3*fac*etorsi*dsin(3*phii) + endif + do j=1,3 + v1ij=v1(j+1,itori,itori1) + v2ij=v2(j+1,itori,itori1) + cosphi=dcos(j*phii) + sinphi=dsin(j*phii) + etors=etors+v1ij*cosphi+v2ij*sinphi+dabs(v1ij)+dabs(v2ij) + gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) + enddo + else + do j=1,nterm_old + v1ij=v1(j,itori,itori1) + v2ij=v2(j,itori,itori1) + cosphi=dcos(j*phii) + sinphi=dsin(j*phii) + etors=etors+v1ij*cosphi+v2ij*sinphi+dabs(v1ij)+dabs(v2ij) + gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) + enddo + endif + if (lprn) + & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') + & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, + & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6) + gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci +c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) + enddo +! 6/20/98 - dihedral angle constraints + edihcnstr=0.0d0 + do i=1,ndih_constr + itori=idih_constr(i) + phii=phi(itori) + difi=phii-phi0(i) + if (difi.gt.drange(i)) then + difi=difi-drange(i) + edihcnstr=edihcnstr+0.25d0*ftors*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + edihcnstr=edihcnstr+0.25d0*ftors*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + endif +! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, +! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) + enddo +! write (iout,*) 'edihcnstr',edihcnstr + return + end +c------------------------------------------------------------------------------ +#else + subroutine etor(etors,edihcnstr) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.VAR' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.TORSION' + include 'COMMON.INTERACT' + include 'COMMON.DERIV' + include 'COMMON.CHAIN' + include 'COMMON.NAMES' + include 'COMMON.IOUNITS' + include 'COMMON.FFIELD' + include 'COMMON.TORCNSTR' + logical lprn +C Set lprn=.true. for debugging + lprn=.false. +c lprn=.true. + etors=0.0D0 + do i=iphi_start,iphi_end + if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215 + itori=itortyp(itype(i-2)) + itori1=itortyp(itype(i-1)) + phii=phi(i) + gloci=0.0D0 +C Regular cosine and sine terms + do j=1,nterm(itori,itori1) + v1ij=v1(j,itori,itori1) + v2ij=v2(j,itori,itori1) + cosphi=dcos(j*phii) + sinphi=dsin(j*phii) + etors=etors+v1ij*cosphi+v2ij*sinphi + gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) + enddo +C Lorentz terms +C v1 +C E = SUM ----------------------------------- - v1 +C [v2 cos(phi/2)+v3 sin(phi/2)]^2 + 1 +C + cosphi=dcos(0.5d0*phii) + sinphi=dsin(0.5d0*phii) + do j=1,nlor(itori,itori1) + vl1ij=vlor1(j,itori,itori1) + vl2ij=vlor2(j,itori,itori1) + vl3ij=vlor3(j,itori,itori1) + pom=vl2ij*cosphi+vl3ij*sinphi + pom1=1.0d0/(pom*pom+1.0d0) + etors=etors+vl1ij*pom1 + pom=-pom*pom1*pom1 + gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom + enddo +C Subtract the constant term + etors=etors-v0(itori,itori1) + if (lprn) + & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') + & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, + & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6) + gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci +c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) + 1215 continue + enddo +! 6/20/98 - dihedral angle constraints + edihcnstr=0.0d0 + do i=1,ndih_constr + print *,"i",i + itori=idih_constr(i) + phii=phi(itori) + difi=phii-phi0(i) + if (difi.gt.drange(i)) then + difi=difi-drange(i) + edihcnstr=edihcnstr+0.25d0*ftors*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + edihcnstr=edihcnstr+0.25d0*ftors*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + endif +! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, +! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) + enddo +! write (iout,*) 'edihcnstr',edihcnstr + return + end +c---------------------------------------------------------------------------- + subroutine etor_d(etors_d) +C 6/23/01 Compute double torsional energy + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.VAR' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.TORSION' + include 'COMMON.INTERACT' + include 'COMMON.DERIV' + include 'COMMON.CHAIN' + include 'COMMON.NAMES' + include 'COMMON.IOUNITS' + include 'COMMON.FFIELD' + include 'COMMON.TORCNSTR' + logical lprn +C Set lprn=.true. for debugging + lprn=.false. +c lprn=.true. + etors_d=0.0D0 + do i=iphi_start,iphi_end-1 + if (itel(i-2).eq.0 .or. itel(i-1).eq.0 .or. itel(i).eq.0) + & goto 1215 + itori=itortyp(itype(i-2)) + itori1=itortyp(itype(i-1)) + itori2=itortyp(itype(i)) + phii=phi(i) + phii1=phi(i+1) + gloci1=0.0D0 + gloci2=0.0D0 +C Regular cosine and sine terms + do j=1,ntermd_1(itori,itori1,itori2) + v1cij=v1c(1,j,itori,itori1,itori2) + v1sij=v1s(1,j,itori,itori1,itori2) + v2cij=v1c(2,j,itori,itori1,itori2) + v2sij=v1s(2,j,itori,itori1,itori2) + cosphi1=dcos(j*phii) + sinphi1=dsin(j*phii) + cosphi2=dcos(j*phii1) + sinphi2=dsin(j*phii1) + etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+ + & v2cij*cosphi2+v2sij*sinphi2 + gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1) + gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2) + enddo + do k=2,ntermd_2(itori,itori1,itori2) + do l=1,k-1 + v1cdij = v2c(k,l,itori,itori1,itori2) + v2cdij = v2c(l,k,itori,itori1,itori2) + v1sdij = v2s(k,l,itori,itori1,itori2) + v2sdij = v2s(l,k,itori,itori1,itori2) + cosphi1p2=dcos(l*phii+(k-l)*phii1) + cosphi1m2=dcos(l*phii-(k-l)*phii1) + sinphi1p2=dsin(l*phii+(k-l)*phii1) + sinphi1m2=dsin(l*phii-(k-l)*phii1) + etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+ + & v1sdij*sinphi1p2+v2sdij*sinphi1m2 + gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2 + & -v1cdij*sinphi1p2-v2cdij*sinphi1m2) + gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2 + & -v1cdij*sinphi1p2+v2cdij*sinphi1m2) + enddo + enddo + gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1 + gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2 + 1215 continue + enddo + return + end +#endif +c------------------------------------------------------------------------------ + subroutine multibody(ecorr) +C This subroutine calculates multi-body contributions to energy following +C the idea of Skolnick et al. If side chains I and J make a contact and +C at the same time side chains I+1 and J+1 make a contact, an extra +C contribution equal to sqrt(eps(i,j)*eps(i+1,j+1)) is added. + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + double precision gx(3),gx1(3) + logical lprn + +C Set lprn=.true. for debugging + lprn=.false. + + if (lprn) then + write (iout,'(a)') 'Contact function values:' + do i=nnt,nct-2 + write (iout,'(i2,20(1x,i2,f10.5))') + & i,(jcont(j,i),facont(j,i),j=1,num_cont(i)) + enddo + endif + ecorr=0.0D0 + do i=nnt,nct + do j=1,3 + gradcorr(j,i)=0.0D0 + gradxorr(j,i)=0.0D0 + enddo + enddo + do i=nnt,nct-2 + + DO ISHIFT = 3,4 + + i1=i+ishift + num_conti=num_cont(i) + num_conti1=num_cont(i1) + do jj=1,num_conti + j=jcont(jj,i) + do kk=1,num_conti1 + j1=jcont(kk,i1) + if (j1.eq.j+ishift .or. j1.eq.j-ishift) then +cd write(iout,*)'i=',i,' j=',j,' i1=',i1,' j1=',j1, +cd & ' ishift=',ishift +C Contacts I--J and I+ISHIFT--J+-ISHIFT1 occur simultaneously. +C The system gains extra energy. + ecorr=ecorr+esccorr(i,j,i1,j1,jj,kk) + endif ! j1==j+-ishift + enddo ! kk + enddo ! jj + + ENDDO ! ISHIFT + + enddo ! i + return + end +c------------------------------------------------------------------------------ + double precision function esccorr(i,j,k,l,jj,kk) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + double precision gx(3),gx1(3) + logical lprn + lprn=.false. + eij=facont(jj,i) + ekl=facont(kk,k) +cd write (iout,'(4i5,3f10.5)') i,j,k,l,eij,ekl,-eij*ekl +C Calculate the multi-body contribution to energy. +C Calculate multi-body contributions to the gradient. +cd write (iout,'(2(2i3,3f10.5))')i,j,(gacont(m,jj,i),m=1,3), +cd & k,l,(gacont(m,kk,k),m=1,3) + do m=1,3 + gx(m) =ekl*gacont(m,jj,i) + gx1(m)=eij*gacont(m,kk,k) + gradxorr(m,i)=gradxorr(m,i)-gx(m) + gradxorr(m,j)=gradxorr(m,j)+gx(m) + gradxorr(m,k)=gradxorr(m,k)-gx1(m) + gradxorr(m,l)=gradxorr(m,l)+gx1(m) + enddo + do m=i,j-1 + do ll=1,3 + gradcorr(ll,m)=gradcorr(ll,m)+gx(ll) + enddo + enddo + do m=k,l-1 + do ll=1,3 + gradcorr(ll,m)=gradcorr(ll,m)+gx1(ll) + enddo + enddo + esccorr=-eij*ekl + return + end +c------------------------------------------------------------------------------ +#ifdef MPL + subroutine pack_buffer(dimen1,dimen2,atom,indx,buffer) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + integer dimen1,dimen2,atom,indx + double precision buffer(dimen1,dimen2) + double precision zapas + common /contacts_hb/ zapas(3,20,maxres,7), + & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres), + & num_cont_hb(maxres),jcont_hb(20,maxres) + num_kont=num_cont_hb(atom) + do i=1,num_kont + do k=1,7 + do j=1,3 + buffer(i,indx+(k-1)*3+j)=zapas(j,i,atom,k) + enddo ! j + enddo ! k + buffer(i,indx+22)=facont_hb(i,atom) + buffer(i,indx+23)=ees0p(i,atom) + buffer(i,indx+24)=ees0m(i,atom) + buffer(i,indx+25)=dfloat(jcont_hb(i,atom)) + enddo ! i + buffer(1,indx+26)=dfloat(num_kont) + return + end +c------------------------------------------------------------------------------ + subroutine unpack_buffer(dimen1,dimen2,atom,indx,buffer) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + integer dimen1,dimen2,atom,indx + double precision buffer(dimen1,dimen2) + double precision zapas + common /contacts_hb/ zapas(3,20,maxres,7), + & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres), + & num_cont_hb(maxres),jcont_hb(20,maxres) + num_kont=buffer(1,indx+26) + num_kont_old=num_cont_hb(atom) + num_cont_hb(atom)=num_kont+num_kont_old + do i=1,num_kont + ii=i+num_kont_old + do k=1,7 + do j=1,3 + zapas(j,ii,atom,k)=buffer(i,indx+(k-1)*3+j) + enddo ! j + enddo ! k + facont_hb(ii,atom)=buffer(i,indx+22) + ees0p(ii,atom)=buffer(i,indx+23) + ees0m(ii,atom)=buffer(i,indx+24) + jcont_hb(ii,atom)=buffer(i,indx+25) + enddo ! i + return + end +c------------------------------------------------------------------------------ +#endif + subroutine multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1) +C This subroutine calculates multi-body contributions to hydrogen-bonding + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.IOUNITS' +#ifdef MPL + include 'COMMON.INFO' +#endif + include 'COMMON.FFIELD' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' +#ifdef MPL + parameter (max_cont=maxconts) + parameter (max_dim=2*(8*3+2)) + parameter (msglen1=max_cont*max_dim*4) + parameter (msglen2=2*msglen1) + integer source,CorrelType,CorrelID,Error + double precision buffer(max_cont,max_dim) +#endif + double precision gx(3),gx1(3) + logical lprn,ldone + +C Set lprn=.true. for debugging + lprn=.false. +#ifdef MPL + n_corr=0 + n_corr1=0 + if (fgProcs.le.1) goto 30 + if (lprn) then + write (iout,'(a)') 'Contact function values:' + do i=nnt,nct-2 + write (iout,'(2i3,50(1x,i2,f5.2))') + & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), + & j=1,num_cont_hb(i)) + enddo + endif +C Caution! Following code assumes that electrostatic interactions concerning +C a given atom are split among at most two processors! + CorrelType=477 + CorrelID=MyID+1 + ldone=.false. + do i=1,max_cont + do j=1,max_dim + buffer(i,j)=0.0D0 + enddo + enddo + mm=mod(MyRank,2) +cd write (iout,*) 'MyRank',MyRank,' mm',mm + if (mm) 20,20,10 + 10 continue +cd write (iout,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone + if (MyRank.gt.0) then +C Send correlation contributions to the preceding processor + msglen=msglen1 + nn=num_cont_hb(iatel_s) + call pack_buffer(max_cont,max_dim,iatel_s,0,buffer) +cd write (iout,*) 'The BUFFER array:' +cd do i=1,nn +cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,26) +cd enddo + if (ielstart(iatel_s).gt.iatel_s+ispp) then + msglen=msglen2 + call pack_buffer(max_cont,max_dim,iatel_s+1,26,buffer) +C Clear the contacts of the atom passed to the neighboring processor + nn=num_cont_hb(iatel_s+1) +cd do i=1,nn +cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j+26),j=1,26) +cd enddo + num_cont_hb(iatel_s)=0 + endif +cd write (iout,*) 'Processor ',MyID,MyRank, +cd & ' is sending correlation contribution to processor',MyID-1, +cd & ' msglen=',msglen +cd write (*,*) 'Processor ',MyID,MyRank, +cd & ' is sending correlation contribution to processor',MyID-1, +cd & ' msglen=',msglen,' CorrelType=',CorrelType + call mp_bsend(buffer,msglen,MyID-1,CorrelType,CorrelID) +cd write (iout,*) 'Processor ',MyID, +cd & ' has sent correlation contribution to processor',MyID-1, +cd & ' msglen=',msglen,' CorrelID=',CorrelID +cd write (*,*) 'Processor ',MyID, +cd & ' has sent correlation contribution to processor',MyID-1, +cd & ' msglen=',msglen,' CorrelID=',CorrelID + msglen=msglen1 + endif ! (MyRank.gt.0) + if (ldone) goto 30 + ldone=.true. + 20 continue +cd write (iout,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone + if (MyRank.lt.fgProcs-1) then +C Receive correlation contributions from the next processor + msglen=msglen1 + if (ielend(iatel_e).lt.nct-1) msglen=msglen2 +cd write (iout,*) 'Processor',MyID, +cd & ' is receiving correlation contribution from processor',MyID+1, +cd & ' msglen=',msglen,' CorrelType=',CorrelType +cd write (*,*) 'Processor',MyID, +cd & ' is receiving correlation contribution from processor',MyID+1, +cd & ' msglen=',msglen,' CorrelType=',CorrelType + nbytes=-1 + do while (nbytes.le.0) + call mp_probe(MyID+1,CorrelType,nbytes) + enddo +cd print *,'Processor',MyID,' msglen',msglen,' nbytes',nbytes + call mp_brecv(buffer,msglen,MyID+1,CorrelType,nbytes) +cd write (iout,*) 'Processor',MyID, +cd & ' has received correlation contribution from processor',MyID+1, +cd & ' msglen=',msglen,' nbytes=',nbytes +cd write (iout,*) 'The received BUFFER array:' +cd do i=1,max_cont +cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,52) +cd enddo + if (msglen.eq.msglen1) then + call unpack_buffer(max_cont,max_dim,iatel_e+1,0,buffer) + else if (msglen.eq.msglen2) then + call unpack_buffer(max_cont,max_dim,iatel_e,0,buffer) + call unpack_buffer(max_cont,max_dim,iatel_e+1,26,buffer) + else + write (iout,*) + & 'ERROR!!!! message length changed while processing correlations.' + write (*,*) + & 'ERROR!!!! message length changed while processing correlations.' + call mp_stopall(Error) + endif ! msglen.eq.msglen1 + endif ! MyRank.lt.fgProcs-1 + if (ldone) goto 30 + ldone=.true. + goto 10 + 30 continue +#endif + if (lprn) then + write (iout,'(a)') 'Contact function values:' + do i=nnt,nct-2 + write (iout,'(2i3,50(1x,i2,f5.2))') + & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), + & j=1,num_cont_hb(i)) + enddo + endif + ecorr=0.0D0 +C Remove the loop below after debugging !!! + do i=nnt,nct + do j=1,3 + gradcorr(j,i)=0.0D0 + gradxorr(j,i)=0.0D0 + enddo + enddo +C Calculate the local-electrostatic correlation terms + do i=iatel_s,iatel_e+1 + i1=i+1 + num_conti=num_cont_hb(i) + num_conti1=num_cont_hb(i+1) + do jj=1,num_conti + j=jcont_hb(jj,i) + do kk=1,num_conti1 + j1=jcont_hb(kk,i1) +c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, +c & ' jj=',jj,' kk=',kk + if (j1.eq.j+1 .or. j1.eq.j-1) then +C Contacts I-J and (I+1)-(J+1) or (I+1)-(J-1) occur simultaneously. +C The system gains extra energy. + ecorr=ecorr+ehbcorr(i,j,i+1,j1,jj,kk,0.72D0,0.32D0) + n_corr=n_corr+1 + else if (j1.eq.j) then +C Contacts I-J and I-(J+1) occur simultaneously. +C The system loses extra energy. +c ecorr=ecorr+ehbcorr(i,j,i+1,j,jj,kk,0.60D0,-0.40D0) + endif + enddo ! kk + do kk=1,num_conti + j1=jcont_hb(kk,i) +c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, +c & ' jj=',jj,' kk=',kk + if (j1.eq.j+1) then +C Contacts I-J and (I+1)-J occur simultaneously. +C The system loses extra energy. +c ecorr=ecorr+ehbcorr(i,j,i,j+1,jj,kk,0.60D0,-0.40D0) + endif ! j1==j+1 + enddo ! kk + enddo ! jj + enddo ! i + return + end +c------------------------------------------------------------------------------ + subroutine multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr, + & n_corr1) +C This subroutine calculates multi-body contributions to hydrogen-bonding + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.IOUNITS' +#ifdef MPL + include 'COMMON.INFO' +#endif + include 'COMMON.FFIELD' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' +#ifdef MPL + parameter (max_cont=maxconts) + parameter (max_dim=2*(8*3+2)) + parameter (msglen1=max_cont*max_dim*4) + parameter (msglen2=2*msglen1) + integer source,CorrelType,CorrelID,Error + double precision buffer(max_cont,max_dim) +#endif + double precision gx(3),gx1(3) + logical lprn,ldone + +C Set lprn=.true. for debugging + lprn=.false. + eturn6=0.0d0 +#ifdef MPL + n_corr=0 + n_corr1=0 + if (fgProcs.le.1) goto 30 + if (lprn) then + write (iout,'(a)') 'Contact function values:' + do i=nnt,nct-2 + write (iout,'(2i3,50(1x,i2,f5.2))') + & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), + & j=1,num_cont_hb(i)) + enddo + endif +C Caution! Following code assumes that electrostatic interactions concerning +C a given atom are split among at most two processors! + CorrelType=477 + CorrelID=MyID+1 + ldone=.false. + do i=1,max_cont + do j=1,max_dim + buffer(i,j)=0.0D0 + enddo + enddo + mm=mod(MyRank,2) +cd write (iout,*) 'MyRank',MyRank,' mm',mm + if (mm) 20,20,10 + 10 continue +cd write (iout,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone + if (MyRank.gt.0) then +C Send correlation contributions to the preceding processor + msglen=msglen1 + nn=num_cont_hb(iatel_s) + call pack_buffer(max_cont,max_dim,iatel_s,0,buffer) +cd write (iout,*) 'The BUFFER array:' +cd do i=1,nn +cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,26) +cd enddo + if (ielstart(iatel_s).gt.iatel_s+ispp) then + msglen=msglen2 + call pack_buffer(max_cont,max_dim,iatel_s+1,26,buffer) +C Clear the contacts of the atom passed to the neighboring processor + nn=num_cont_hb(iatel_s+1) +cd do i=1,nn +cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j+26),j=1,26) +cd enddo + num_cont_hb(iatel_s)=0 + endif +cd write (iout,*) 'Processor ',MyID,MyRank, +cd & ' is sending correlation contribution to processor',MyID-1, +cd & ' msglen=',msglen +cd write (*,*) 'Processor ',MyID,MyRank, +cd & ' is sending correlation contribution to processor',MyID-1, +cd & ' msglen=',msglen,' CorrelType=',CorrelType + call mp_bsend(buffer,msglen,MyID-1,CorrelType,CorrelID) +cd write (iout,*) 'Processor ',MyID, +cd & ' has sent correlation contribution to processor',MyID-1, +cd & ' msglen=',msglen,' CorrelID=',CorrelID +cd write (*,*) 'Processor ',MyID, +cd & ' has sent correlation contribution to processor',MyID-1, +cd & ' msglen=',msglen,' CorrelID=',CorrelID + msglen=msglen1 + endif ! (MyRank.gt.0) + if (ldone) goto 30 + ldone=.true. + 20 continue +cd write (iout,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone + if (MyRank.lt.fgProcs-1) then +C Receive correlation contributions from the next processor + msglen=msglen1 + if (ielend(iatel_e).lt.nct-1) msglen=msglen2 +cd write (iout,*) 'Processor',MyID, +cd & ' is receiving correlation contribution from processor',MyID+1, +cd & ' msglen=',msglen,' CorrelType=',CorrelType +cd write (*,*) 'Processor',MyID, +cd & ' is receiving correlation contribution from processor',MyID+1, +cd & ' msglen=',msglen,' CorrelType=',CorrelType + nbytes=-1 + do while (nbytes.le.0) + call mp_probe(MyID+1,CorrelType,nbytes) + enddo +cd print *,'Processor',MyID,' msglen',msglen,' nbytes',nbytes + call mp_brecv(buffer,msglen,MyID+1,CorrelType,nbytes) +cd write (iout,*) 'Processor',MyID, +cd & ' has received correlation contribution from processor',MyID+1, +cd & ' msglen=',msglen,' nbytes=',nbytes +cd write (iout,*) 'The received BUFFER array:' +cd do i=1,max_cont +cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,52) +cd enddo + if (msglen.eq.msglen1) then + call unpack_buffer(max_cont,max_dim,iatel_e+1,0,buffer) + else if (msglen.eq.msglen2) then + call unpack_buffer(max_cont,max_dim,iatel_e,0,buffer) + call unpack_buffer(max_cont,max_dim,iatel_e+1,26,buffer) + else + write (iout,*) + & 'ERROR!!!! message length changed while processing correlations.' + write (*,*) + & 'ERROR!!!! message length changed while processing correlations.' + call mp_stopall(Error) + endif ! msglen.eq.msglen1 + endif ! MyRank.lt.fgProcs-1 + if (ldone) goto 30 + ldone=.true. + goto 10 + 30 continue +#endif + if (lprn) then + write (iout,'(a)') 'Contact function values:' + do i=nnt,nct-2 + write (iout,'(2i3,50(1x,i2,f5.2))') + & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), + & j=1,num_cont_hb(i)) + enddo + endif + ecorr=0.0D0 + ecorr5=0.0d0 + ecorr6=0.0d0 +C Remove the loop below after debugging !!! + do i=nnt,nct + do j=1,3 + gradcorr(j,i)=0.0D0 + gradxorr(j,i)=0.0D0 + enddo + enddo +C Calculate the dipole-dipole interaction energies + if (wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) then + do i=iatel_s,iatel_e+1 + num_conti=num_cont_hb(i) + do jj=1,num_conti + j=jcont_hb(jj,i) + call dipole(i,j,jj) + enddo + enddo + endif +C Calculate the local-electrostatic correlation terms + do i=iatel_s,iatel_e+1 + i1=i+1 + num_conti=num_cont_hb(i) + num_conti1=num_cont_hb(i+1) + do jj=1,num_conti + j=jcont_hb(jj,i) + do kk=1,num_conti1 + j1=jcont_hb(kk,i1) +c write (*,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, +c & ' jj=',jj,' kk=',kk + if (j1.eq.j+1 .or. j1.eq.j-1) then +C Contacts I-J and (I+1)-(J+1) or (I+1)-(J-1) occur simultaneously. +C The system gains extra energy. + n_corr=n_corr+1 + sqd1=dsqrt(d_cont(jj,i)) + sqd2=dsqrt(d_cont(kk,i1)) + sred_geom = sqd1*sqd2 + IF (sred_geom.lt.cutoff_corr) THEN + call gcont(sred_geom,r0_corr,1.0D0,delt_corr, + & ekont,fprimcont) +c write (*,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, +c & ' jj=',jj,' kk=',kk + fac_prim1=0.5d0*sqd2/sqd1*fprimcont + fac_prim2=0.5d0*sqd1/sqd2*fprimcont + do l=1,3 + g_contij(l,1)=fac_prim1*grij_hb_cont(l,jj,i) + g_contij(l,2)=fac_prim2*grij_hb_cont(l,kk,i1) + enddo + n_corr1=n_corr1+1 +cd write (iout,*) 'sred_geom=',sred_geom, +cd & ' ekont=',ekont,' fprim=',fprimcont + call calc_eello(i,j,i+1,j1,jj,kk) + if (wcorr4.gt.0.0d0) + & ecorr=ecorr+eello4(i,j,i+1,j1,jj,kk) + if (wcorr5.gt.0.0d0) + & ecorr5=ecorr5+eello5(i,j,i+1,j1,jj,kk) +c print *,"wcorr5",ecorr5 +cd write(2,*)'wcorr6',wcorr6,' wturn6',wturn6 +cd write(2,*)'ijkl',i,j,i+1,j1 + if (wcorr6.gt.0.0d0 .and. (j.ne.i+4 .or. j1.ne.i+3 + & .or. wturn6.eq.0.0d0))then +cd write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1 + ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk) +cd write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5, +cd & 'ecorr6=',ecorr6 +cd write (iout,'(4e15.5)') sred_geom, +cd & dabs(eello4(i,j,i+1,j1,jj,kk)), +cd & dabs(eello5(i,j,i+1,j1,jj,kk)), +cd & dabs(eello6(i,j,i+1,j1,jj,kk)) + else if (wturn6.gt.0.0d0 + & .and. (j.eq.i+4 .and. j1.eq.i+3)) then +cd write (iout,*) '******eturn6: i,j,i+1,j1',i,j,i+1,j1 + eturn6=eturn6+eello_turn6(i,jj,kk) +cd write (2,*) 'multibody_eello:eturn6',eturn6 + endif + ENDIF +1111 continue + else if (j1.eq.j) then +C Contacts I-J and I-(J+1) occur simultaneously. +C The system loses extra energy. +c ecorr=ecorr+ehbcorr(i,j,i+1,j,jj,kk,0.60D0,-0.40D0) + endif + enddo ! kk + do kk=1,num_conti + j1=jcont_hb(kk,i) +c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, +c & ' jj=',jj,' kk=',kk + if (j1.eq.j+1) then +C Contacts I-J and (I+1)-J occur simultaneously. +C The system loses extra energy. +c ecorr=ecorr+ehbcorr(i,j,i,j+1,jj,kk,0.60D0,-0.40D0) + endif ! j1==j+1 + enddo ! kk + enddo ! jj + enddo ! i + return + end +c------------------------------------------------------------------------------ + double precision function ehbcorr(i,j,k,l,jj,kk,coeffp,coeffm) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + double precision gx(3),gx1(3) + logical lprn + lprn=.false. + eij=facont_hb(jj,i) + ekl=facont_hb(kk,k) + ees0pij=ees0p(jj,i) + ees0pkl=ees0p(kk,k) + ees0mij=ees0m(jj,i) + ees0mkl=ees0m(kk,k) + ekont=eij*ekl + ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl) +cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl) +C Following 4 lines for diagnostics. +cd ees0pkl=0.0D0 +cd ees0pij=1.0D0 +cd ees0mkl=0.0D0 +cd ees0mij=1.0D0 +c write (iout,*)'Contacts have occurred for peptide groups',i,j, +c & ' and',k,l +c write (iout,*)'Contacts have occurred for peptide groups', +c & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l +c & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees +C Calculate the multi-body contribution to energy. + ecorr=ecorr+ekont*ees + if (calc_grad) then +C Calculate multi-body contributions to the gradient. + do ll=1,3 + ghalf=0.5D0*ees*ekl*gacont_hbr(ll,jj,i) + gradcorr(ll,i)=gradcorr(ll,i)+ghalf + & -ekont*(coeffp*ees0pkl*gacontp_hb1(ll,jj,i)+ + & coeffm*ees0mkl*gacontm_hb1(ll,jj,i)) + gradcorr(ll,j)=gradcorr(ll,j)+ghalf + & -ekont*(coeffp*ees0pkl*gacontp_hb2(ll,jj,i)+ + & coeffm*ees0mkl*gacontm_hb2(ll,jj,i)) + ghalf=0.5D0*ees*eij*gacont_hbr(ll,kk,k) + gradcorr(ll,k)=gradcorr(ll,k)+ghalf + & -ekont*(coeffp*ees0pij*gacontp_hb1(ll,kk,k)+ + & coeffm*ees0mij*gacontm_hb1(ll,kk,k)) + gradcorr(ll,l)=gradcorr(ll,l)+ghalf + & -ekont*(coeffp*ees0pij*gacontp_hb2(ll,kk,k)+ + & coeffm*ees0mij*gacontm_hb2(ll,kk,k)) + enddo + do m=i+1,j-1 + do ll=1,3 + gradcorr(ll,m)=gradcorr(ll,m)+ + & ees*ekl*gacont_hbr(ll,jj,i)- + & ekont*(coeffp*ees0pkl*gacontp_hb3(ll,jj,i)+ + & coeffm*ees0mkl*gacontm_hb3(ll,jj,i)) + enddo + enddo + do m=k+1,l-1 + do ll=1,3 + gradcorr(ll,m)=gradcorr(ll,m)+ + & ees*eij*gacont_hbr(ll,kk,k)- + & ekont*(coeffp*ees0pij*gacontp_hb3(ll,kk,k)+ + & coeffm*ees0mij*gacontm_hb3(ll,kk,k)) + enddo + enddo + endif + ehbcorr=ekont*ees + return + end +C--------------------------------------------------------------------------- + subroutine dipole(i,j,jj) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.FFIELD' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + include 'COMMON.TORSION' + include 'COMMON.VAR' + include 'COMMON.GEO' + dimension dipi(2,2),dipj(2,2),dipderi(2),dipderj(2),auxvec(2), + & auxmat(2,2) + iti1 = itortyp(itype(i+1)) + if (j.lt.nres-1) then + itj1 = itortyp(itype(j+1)) + else + itj1=ntortyp+1 + endif + do iii=1,2 + dipi(iii,1)=Ub2(iii,i) + dipderi(iii)=Ub2der(iii,i) + dipi(iii,2)=b1(iii,iti1) + dipj(iii,1)=Ub2(iii,j) + dipderj(iii)=Ub2der(iii,j) + dipj(iii,2)=b1(iii,itj1) + enddo + kkk=0 + do iii=1,2 + call matvec2(a_chuj(1,1,jj,i),dipj(1,iii),auxvec(1)) + do jjj=1,2 + kkk=kkk+1 + dip(kkk,jj,i)=scalar2(dipi(1,jjj),auxvec(1)) + enddo + enddo + if (.not.calc_grad) return + do kkk=1,5 + do lll=1,3 + mmm=0 + do iii=1,2 + call matvec2(a_chuj_der(1,1,lll,kkk,jj,i),dipj(1,iii), + & auxvec(1)) + do jjj=1,2 + mmm=mmm+1 + dipderx(lll,kkk,mmm,jj,i)=scalar2(dipi(1,jjj),auxvec(1)) + enddo + enddo + enddo + enddo + call transpose2(a_chuj(1,1,jj,i),auxmat(1,1)) + call matvec2(auxmat(1,1),dipderi(1),auxvec(1)) + do iii=1,2 + dipderg(iii,jj,i)=scalar2(auxvec(1),dipj(1,iii)) + enddo + call matvec2(a_chuj(1,1,jj,i),dipderj(1),auxvec(1)) + do iii=1,2 + dipderg(iii+2,jj,i)=scalar2(auxvec(1),dipi(1,iii)) + enddo + return + end +C--------------------------------------------------------------------------- + subroutine calc_eello(i,j,k,l,jj,kk) +C +C This subroutine computes matrices and vectors needed to calculate +C the fourth-, fifth-, and sixth-order local-electrostatic terms. +C + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + include 'COMMON.TORSION' + include 'COMMON.VAR' + include 'COMMON.GEO' + include 'COMMON.FFIELD' + double precision aa1(2,2),aa2(2,2),aa1t(2,2),aa2t(2,2), + & aa1tder(2,2,3,5),aa2tder(2,2,3,5),auxmat(2,2) + logical lprn + common /kutas/ lprn +cd write (iout,*) 'calc_eello: i=',i,' j=',j,' k=',k,' l=',l, +cd & ' jj=',jj,' kk=',kk +cd if (i.ne.2 .or. j.ne.4 .or. k.ne.3 .or. l.ne.5) return + do iii=1,2 + do jjj=1,2 + aa1(iii,jjj)=a_chuj(iii,jjj,jj,i) + aa2(iii,jjj)=a_chuj(iii,jjj,kk,k) + enddo + enddo + call transpose2(aa1(1,1),aa1t(1,1)) + call transpose2(aa2(1,1),aa2t(1,1)) + do kkk=1,5 + do lll=1,3 + call transpose2(a_chuj_der(1,1,lll,kkk,jj,i), + & aa1tder(1,1,lll,kkk)) + call transpose2(a_chuj_der(1,1,lll,kkk,kk,k), + & aa2tder(1,1,lll,kkk)) + enddo + enddo + if (l.eq.j+1) then +C parallel orientation of the two CA-CA-CA frames. + if (i.gt.1) then + iti=itortyp(itype(i)) + else + iti=ntortyp+1 + endif + itk1=itortyp(itype(k+1)) + itj=itortyp(itype(j)) + if (l.lt.nres-1) then + itl1=itortyp(itype(l+1)) + else + itl1=ntortyp+1 + endif +C A1 kernel(j+1) A2T +cd do iii=1,2 +cd write (iout,'(3f10.5,5x,3f10.5)') +cd & (EUg(iii,jjj,k),jjj=1,2),(EUg(iii,jjj,l),jjj=1,2) +cd enddo + call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), + & aa2tder(1,1,1,1),1,.false.,EUg(1,1,l),EUgder(1,1,l), + & AEA(1,1,1),AEAderg(1,1,1),AEAderx(1,1,1,1,1,1)) +C Following matrices are needed only for 6-th order cumulants + IF (wcorr6.gt.0.0d0) THEN + call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), + & aa2tder(1,1,1,1),1,.false.,EUgC(1,1,l),EUgCder(1,1,l), + & AECA(1,1,1),AECAderg(1,1,1),AECAderx(1,1,1,1,1,1)) + call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), + & aa2tder(1,1,1,1),2,.false.,Ug2DtEUg(1,1,l), + & Ug2DtEUgder(1,1,1,l),ADtEA(1,1,1),ADtEAderg(1,1,1,1), + & ADtEAderx(1,1,1,1,1,1)) + lprn=.false. + call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), + & aa2tder(1,1,1,1),2,.false.,DtUg2EUg(1,1,l), + & DtUg2EUgder(1,1,1,l),ADtEA1(1,1,1),ADtEA1derg(1,1,1,1), + & ADtEA1derx(1,1,1,1,1,1)) + ENDIF +C End 6-th order cumulants +cd lprn=.false. +cd if (lprn) then +cd write (2,*) 'In calc_eello6' +cd do iii=1,2 +cd write (2,*) 'iii=',iii +cd do kkk=1,5 +cd write (2,*) 'kkk=',kkk +cd do jjj=1,2 +cd write (2,'(3(2f10.5),5x)') +cd & ((ADtEA1derx(jjj,mmm,lll,kkk,iii,1),mmm=1,2),lll=1,3) +cd enddo +cd enddo +cd enddo +cd endif + call transpose2(EUgder(1,1,k),auxmat(1,1)) + call matmat2(auxmat(1,1),AEA(1,1,1),EAEAderg(1,1,1,1)) + call transpose2(EUg(1,1,k),auxmat(1,1)) + call matmat2(auxmat(1,1),AEA(1,1,1),EAEA(1,1,1)) + call matmat2(auxmat(1,1),AEAderg(1,1,1),EAEAderg(1,1,2,1)) + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,1), + & EAEAderx(1,1,lll,kkk,iii,1)) + enddo + enddo + enddo +C A1T kernel(i+1) A2 + call kernel(aa1t(1,1),aa2(1,1),aa1tder(1,1,1,1), + & a_chuj_der(1,1,1,1,kk,k),1,.false.,EUg(1,1,k),EUgder(1,1,k), + & AEA(1,1,2),AEAderg(1,1,2),AEAderx(1,1,1,1,1,2)) +C Following matrices are needed only for 6-th order cumulants + IF (wcorr6.gt.0.0d0) THEN + call kernel(aa1t(1,1),aa2(1,1),aa1tder(1,1,1,1), + & a_chuj_der(1,1,1,1,kk,k),1,.false.,EUgC(1,1,k),EUgCder(1,1,k), + & AECA(1,1,2),AECAderg(1,1,2),AECAderx(1,1,1,1,1,2)) + call kernel(aa1t(1,1),aa2(1,1),aa1tder(1,1,1,1), + & a_chuj_der(1,1,1,1,kk,k),2,.false.,Ug2DtEUg(1,1,k), + & Ug2DtEUgder(1,1,1,k),ADtEA(1,1,2),ADtEAderg(1,1,1,2), + & ADtEAderx(1,1,1,1,1,2)) + call kernel(aa1t(1,1),aa2(1,1),aa1tder(1,1,1,1), + & a_chuj_der(1,1,1,1,kk,k),2,.false.,DtUg2EUg(1,1,k), + & DtUg2EUgder(1,1,1,k),ADtEA1(1,1,2),ADtEA1derg(1,1,1,2), + & ADtEA1derx(1,1,1,1,1,2)) + ENDIF +C End 6-th order cumulants + call transpose2(EUgder(1,1,l),auxmat(1,1)) + call matmat2(auxmat(1,1),AEA(1,1,2),EAEAderg(1,1,1,2)) + call transpose2(EUg(1,1,l),auxmat(1,1)) + call matmat2(auxmat(1,1),AEA(1,1,2),EAEA(1,1,2)) + call matmat2(auxmat(1,1),AEAderg(1,1,2),EAEAderg(1,1,2,2)) + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,2), + & EAEAderx(1,1,lll,kkk,iii,2)) + enddo + enddo + enddo +C AEAb1 and AEAb2 +C Calculate the vectors and their derivatives in virtual-bond dihedral angles. +C They are needed only when the fifth- or the sixth-order cumulants are +C indluded. + IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) THEN + call transpose2(AEA(1,1,1),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1)) + call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1)) + call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1)) + call transpose2(AEAderg(1,1,1),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1)) + call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1)) + call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1)) + call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1)) + call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1)) + call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1)) + call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1)) + call transpose2(AEA(1,1,2),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,itj),AEAb1(1,1,2)) + call matvec2(auxmat(1,1),Ub2(1,j),AEAb2(1,1,2)) + call matvec2(auxmat(1,1),Ub2der(1,j),AEAb2derg(1,2,1,2)) + call transpose2(AEAderg(1,1,2),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,itj),AEAb1derg(1,1,2)) + call matvec2(auxmat(1,1),Ub2(1,j),AEAb2derg(1,1,1,2)) + call matvec2(AEA(1,1,2),b1(1,itl1),AEAb1(1,2,2)) + call matvec2(AEAderg(1,1,2),b1(1,itl1),AEAb1derg(1,2,2)) + call matvec2(AEA(1,1,2),Ub2(1,l+1),AEAb2(1,2,2)) + call matvec2(AEAderg(1,1,2),Ub2(1,l+1),AEAb2derg(1,1,2,2)) + call matvec2(AEA(1,1,2),Ub2der(1,l+1),AEAb2derg(1,2,2,2)) +C Calculate the Cartesian derivatives of the vectors. + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,iti), + & AEAb1derx(1,lll,kkk,iii,1,1)) + call matvec2(auxmat(1,1),Ub2(1,i), + & AEAb2derx(1,lll,kkk,iii,1,1)) + call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1), + & AEAb1derx(1,lll,kkk,iii,2,1)) + call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1), + & AEAb2derx(1,lll,kkk,iii,2,1)) + call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,itj), + & AEAb1derx(1,lll,kkk,iii,1,2)) + call matvec2(auxmat(1,1),Ub2(1,j), + & AEAb2derx(1,lll,kkk,iii,1,2)) + call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itl1), + & AEAb1derx(1,lll,kkk,iii,2,2)) + call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,l+1), + & AEAb2derx(1,lll,kkk,iii,2,2)) + enddo + enddo + enddo + ENDIF +C End vectors + else +C Antiparallel orientation of the two CA-CA-CA frames. + if (i.gt.1) then + iti=itortyp(itype(i)) + else + iti=ntortyp+1 + endif + itk1=itortyp(itype(k+1)) + itl=itortyp(itype(l)) + itj=itortyp(itype(j)) + if (j.lt.nres-1) then + itj1=itortyp(itype(j+1)) + else + itj1=ntortyp+1 + endif +C A2 kernel(j-1)T A1T + call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), + & aa2tder(1,1,1,1),1,.true.,EUg(1,1,j),EUgder(1,1,j), + & AEA(1,1,1),AEAderg(1,1,1),AEAderx(1,1,1,1,1,1)) +C Following matrices are needed only for 6-th order cumulants + IF (wcorr6.gt.0.0d0 .or. (wturn6.gt.0.0d0 .and. + & j.eq.i+4 .and. l.eq.i+3)) THEN + call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), + & aa2tder(1,1,1,1),1,.true.,EUgC(1,1,j),EUgCder(1,1,j), + & AECA(1,1,1),AECAderg(1,1,1),AECAderx(1,1,1,1,1,1)) + call kernel(aa2(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), + & aa2tder(1,1,1,1),2,.true.,Ug2DtEUg(1,1,j), + & Ug2DtEUgder(1,1,1,j),ADtEA(1,1,1),ADtEAderg(1,1,1,1), + & ADtEAderx(1,1,1,1,1,1)) + call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), + & aa2tder(1,1,1,1),2,.true.,DtUg2EUg(1,1,j), + & DtUg2EUgder(1,1,1,j),ADtEA1(1,1,1),ADtEA1derg(1,1,1,1), + & ADtEA1derx(1,1,1,1,1,1)) + ENDIF +C End 6-th order cumulants + call transpose2(EUgder(1,1,k),auxmat(1,1)) + call matmat2(auxmat(1,1),AEA(1,1,1),EAEAderg(1,1,1,1)) + call transpose2(EUg(1,1,k),auxmat(1,1)) + call matmat2(auxmat(1,1),AEA(1,1,1),EAEA(1,1,1)) + call matmat2(auxmat(1,1),AEAderg(1,1,1),EAEAderg(1,1,2,1)) + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,1), + & EAEAderx(1,1,lll,kkk,iii,1)) + enddo + enddo + enddo +C A2T kernel(i+1)T A1 + call kernel(aa2t(1,1),aa1(1,1),aa2tder(1,1,1,1), + & a_chuj_der(1,1,1,1,jj,i),1,.true.,EUg(1,1,k),EUgder(1,1,k), + & AEA(1,1,2),AEAderg(1,1,2),AEAderx(1,1,1,1,1,2)) +C Following matrices are needed only for 6-th order cumulants + IF (wcorr6.gt.0.0d0 .or. (wturn6.gt.0.0d0 .and. + & j.eq.i+4 .and. l.eq.i+3)) THEN + call kernel(aa2t(1,1),aa1(1,1),aa2tder(1,1,1,1), + & a_chuj_der(1,1,1,1,jj,i),1,.true.,EUgC(1,1,k),EUgCder(1,1,k), + & AECA(1,1,2),AECAderg(1,1,2),AECAderx(1,1,1,1,1,2)) + call kernel(aa2t(1,1),aa1(1,1),aa2tder(1,1,1,1), + & a_chuj_der(1,1,1,1,jj,i),2,.true.,Ug2DtEUg(1,1,k), + & Ug2DtEUgder(1,1,1,k),ADtEA(1,1,2),ADtEAderg(1,1,1,2), + & ADtEAderx(1,1,1,1,1,2)) + call kernel(aa2t(1,1),aa1(1,1),aa2tder(1,1,1,1), + & a_chuj_der(1,1,1,1,jj,i),2,.true.,DtUg2EUg(1,1,k), + & DtUg2EUgder(1,1,1,k),ADtEA1(1,1,2),ADtEA1derg(1,1,1,2), + & ADtEA1derx(1,1,1,1,1,2)) + ENDIF +C End 6-th order cumulants + call transpose2(EUgder(1,1,j),auxmat(1,1)) + call matmat2(auxmat(1,1),AEA(1,1,1),EAEAderg(1,1,2,2)) + call transpose2(EUg(1,1,j),auxmat(1,1)) + call matmat2(auxmat(1,1),AEA(1,1,2),EAEA(1,1,2)) + call matmat2(auxmat(1,1),AEAderg(1,1,2),EAEAderg(1,1,2,2)) + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,2), + & EAEAderx(1,1,lll,kkk,iii,2)) + enddo + enddo + enddo +C AEAb1 and AEAb2 +C Calculate the vectors and their derivatives in virtual-bond dihedral angles. +C They are needed only when the fifth- or the sixth-order cumulants are +C indluded. + IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 .or. + & (wturn6.gt.0.0d0 .and. j.eq.i+4 .and. l.eq.i+3)) THEN + call transpose2(AEA(1,1,1),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1)) + call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1)) + call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1)) + call transpose2(AEAderg(1,1,1),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1)) + call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1)) + call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1)) + call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1)) + call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1)) + call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1)) + call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1)) + call transpose2(AEA(1,1,2),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,itj1),AEAb1(1,1,2)) + call matvec2(auxmat(1,1),Ub2(1,l),AEAb2(1,1,2)) + call matvec2(auxmat(1,1),Ub2der(1,l),AEAb2derg(1,2,1,2)) + call transpose2(AEAderg(1,1,2),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,itl),AEAb1(1,1,2)) + call matvec2(auxmat(1,1),Ub2(1,l),AEAb2derg(1,1,1,2)) + call matvec2(AEA(1,1,2),b1(1,itj1),AEAb1(1,2,2)) + call matvec2(AEAderg(1,1,2),b1(1,itj1),AEAb1derg(1,2,2)) + call matvec2(AEA(1,1,2),Ub2(1,j),AEAb2(1,2,2)) + call matvec2(AEAderg(1,1,2),Ub2(1,j),AEAb2derg(1,1,2,2)) + call matvec2(AEA(1,1,2),Ub2der(1,j),AEAb2derg(1,2,2,2)) +C Calculate the Cartesian derivatives of the vectors. + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,iti), + & AEAb1derx(1,lll,kkk,iii,1,1)) + call matvec2(auxmat(1,1),Ub2(1,i), + & AEAb2derx(1,lll,kkk,iii,1,1)) + call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1), + & AEAb1derx(1,lll,kkk,iii,2,1)) + call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1), + & AEAb2derx(1,lll,kkk,iii,2,1)) + call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1)) + call matvec2(auxmat(1,1),b1(1,itl), + & AEAb1derx(1,lll,kkk,iii,1,2)) + call matvec2(auxmat(1,1),Ub2(1,l), + & AEAb2derx(1,lll,kkk,iii,1,2)) + call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itj1), + & AEAb1derx(1,lll,kkk,iii,2,2)) + call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,j), + & AEAb2derx(1,lll,kkk,iii,2,2)) + enddo + enddo + enddo + ENDIF +C End vectors + endif + return + end +C--------------------------------------------------------------------------- + subroutine kernel(aa1,aa2t,aa1derx,aa2tderx,nderg,transp, + & KK,KKderg,AKA,AKAderg,AKAderx) + implicit none + integer nderg + logical transp + double precision aa1(2,2),aa2t(2,2),aa1derx(2,2,3,5), + & aa2tderx(2,2,3,5),KK(2,2),KKderg(2,2,nderg),AKA(2,2), + & AKAderg(2,2,nderg),AKAderx(2,2,3,5,2) + integer iii,kkk,lll + integer jjj,mmm + logical lprn + common /kutas/ lprn + call prodmat3(aa1(1,1),aa2t(1,1),KK(1,1),transp,AKA(1,1)) + do iii=1,nderg + call prodmat3(aa1(1,1),aa2t(1,1),KKderg(1,1,iii),transp, + & AKAderg(1,1,iii)) + enddo +cd if (lprn) write (2,*) 'In kernel' + do kkk=1,5 +cd if (lprn) write (2,*) 'kkk=',kkk + do lll=1,3 + call prodmat3(aa1derx(1,1,lll,kkk),aa2t(1,1), + & KK(1,1),transp,AKAderx(1,1,lll,kkk,1)) +cd if (lprn) then +cd write (2,*) 'lll=',lll +cd write (2,*) 'iii=1' +cd do jjj=1,2 +cd write (2,'(3(2f10.5),5x)') +cd & (AKAderx(jjj,mmm,lll,kkk,1),mmm=1,2) +cd enddo +cd endif + call prodmat3(aa1(1,1),aa2tderx(1,1,lll,kkk), + & KK(1,1),transp,AKAderx(1,1,lll,kkk,2)) +cd if (lprn) then +cd write (2,*) 'lll=',lll +cd write (2,*) 'iii=2' +cd do jjj=1,2 +cd write (2,'(3(2f10.5),5x)') +cd & (AKAderx(jjj,mmm,lll,kkk,2),mmm=1,2) +cd enddo +cd endif + enddo + enddo + return + end +C--------------------------------------------------------------------------- + double precision function eello4(i,j,k,l,jj,kk) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + include 'COMMON.TORSION' + include 'COMMON.VAR' + include 'COMMON.GEO' + double precision pizda(2,2),ggg1(3),ggg2(3) +cd if (i.ne.1 .or. j.ne.5 .or. k.ne.2 .or.l.ne.4) then +cd eello4=0.0d0 +cd return +cd endif +cd print *,'eello4:',i,j,k,l,jj,kk +cd write (2,*) 'i',i,' j',j,' k',k,' l',l +cd call checkint4(i,j,k,l,jj,kk,eel4_num) +cold eij=facont_hb(jj,i) +cold ekl=facont_hb(kk,k) +cold ekont=eij*ekl + eel4=-EAEA(1,1,1)-EAEA(2,2,1) + if (calc_grad) then +cd eel41=-EAEA(1,1,2)-EAEA(2,2,2) + gcorr_loc(k-1)=gcorr_loc(k-1) + & -ekont*(EAEAderg(1,1,1,1)+EAEAderg(2,2,1,1)) + if (l.eq.j+1) then + gcorr_loc(l-1)=gcorr_loc(l-1) + & -ekont*(EAEAderg(1,1,2,1)+EAEAderg(2,2,2,1)) + else + gcorr_loc(j-1)=gcorr_loc(j-1) + & -ekont*(EAEAderg(1,1,2,1)+EAEAderg(2,2,2,1)) + endif + do iii=1,2 + do kkk=1,5 + do lll=1,3 + derx(lll,kkk,iii)=-EAEAderx(1,1,lll,kkk,iii,1) + & -EAEAderx(2,2,lll,kkk,iii,1) +cd derx(lll,kkk,iii)=0.0d0 + enddo + enddo + enddo +cd gcorr_loc(l-1)=0.0d0 +cd gcorr_loc(j-1)=0.0d0 +cd gcorr_loc(k-1)=0.0d0 +cd eel4=1.0d0 +cd write (iout,*)'Contacts have occurred for peptide groups', +cd & i,j,' fcont:',eij,' eij',' and ',k,l, +cd & ' fcont ',ekl,' eel4=',eel4,' eel4_num',16*eel4_num + if (j.lt.nres-1) then + j1=j+1 + j2=j-1 + else + j1=j-1 + j2=j-2 + endif + if (l.lt.nres-1) then + l1=l+1 + l2=l-1 + else + l1=l-1 + l2=l-2 + endif + do ll=1,3 +cold ghalf=0.5d0*eel4*ekl*gacont_hbr(ll,jj,i) + ggg1(ll)=eel4*g_contij(ll,1) + ggg2(ll)=eel4*g_contij(ll,2) + ghalf=0.5d0*ggg1(ll) +cd ghalf=0.0d0 + gradcorr(ll,i)=gradcorr(ll,i)+ghalf+ekont*derx(ll,2,1) + gradcorr(ll,i+1)=gradcorr(ll,i+1)+ekont*derx(ll,3,1) + gradcorr(ll,j)=gradcorr(ll,j)+ghalf+ekont*derx(ll,4,1) + gradcorr(ll,j1)=gradcorr(ll,j1)+ekont*derx(ll,5,1) +cold ghalf=0.5d0*eel4*eij*gacont_hbr(ll,kk,k) + ghalf=0.5d0*ggg2(ll) +cd ghalf=0.0d0 + gradcorr(ll,k)=gradcorr(ll,k)+ghalf+ekont*derx(ll,2,2) + gradcorr(ll,k+1)=gradcorr(ll,k+1)+ekont*derx(ll,3,2) + gradcorr(ll,l)=gradcorr(ll,l)+ghalf+ekont*derx(ll,4,2) + gradcorr(ll,l1)=gradcorr(ll,l1)+ekont*derx(ll,5,2) + enddo +cd goto 1112 + do m=i+1,j-1 + do ll=1,3 +cold gradcorr(ll,m)=gradcorr(ll,m)+eel4*ekl*gacont_hbr(ll,jj,i) + gradcorr(ll,m)=gradcorr(ll,m)+ggg1(ll) + enddo + enddo + do m=k+1,l-1 + do ll=1,3 +cold gradcorr(ll,m)=gradcorr(ll,m)+eel4*eij*gacont_hbr(ll,kk,k) + gradcorr(ll,m)=gradcorr(ll,m)+ggg2(ll) + enddo + enddo +1112 continue + do m=i+2,j2 + do ll=1,3 + gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,1) + enddo + enddo + do m=k+2,l2 + do ll=1,3 + gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,2) + enddo + enddo +cd do iii=1,nres-3 +cd write (2,*) iii,gcorr_loc(iii) +cd enddo + endif + eello4=ekont*eel4 +cd write (2,*) 'ekont',ekont +cd write (iout,*) 'eello4',ekont*eel4 + return + end +C--------------------------------------------------------------------------- + double precision function eello5(i,j,k,l,jj,kk) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + include 'COMMON.TORSION' + include 'COMMON.VAR' + include 'COMMON.GEO' + double precision pizda(2,2),auxmat(2,2),auxmat1(2,2),vv(2) + double precision ggg1(3),ggg2(3) +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C C +C Parallel chains C +C C +C o o o o C +C /l\ / \ \ / \ / \ / C +C / \ / \ \ / \ / \ / C +C j| o |l1 | o | o| o | | o |o C +C \ |/k\| |/ \| / |/ \| |/ \| C +C \i/ \ / \ / / \ / \ C +C o k1 o C +C (I) (II) (III) (IV) C +C C +C eello5_1 eello5_2 eello5_3 eello5_4 C +C C +C Antiparallel chains C +C C +C o o o o C +C /j\ / \ \ / \ / \ / C +C / \ / \ \ / \ / \ / C +C j1| o |l | o | o| o | | o |o C +C \ |/k\| |/ \| / |/ \| |/ \| C +C \i/ \ / \ / / \ / \ C +C o k1 o C +C (I) (II) (III) (IV) C +C C +C eello5_1 eello5_2 eello5_3 eello5_4 C +C C +C o denotes a local interaction, vertical lines an electrostatic interaction. C +C C +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +cd if (i.ne.2 .or. j.ne.6 .or. k.ne.3 .or. l.ne.5) then +cd eello5=0.0d0 +cd return +cd endif +cd write (iout,*) +cd & 'EELLO5: Contacts have occurred for peptide groups',i,j, +cd & ' and',k,l + itk=itortyp(itype(k)) + itl=itortyp(itype(l)) + itj=itortyp(itype(j)) + eello5_1=0.0d0 + eello5_2=0.0d0 + eello5_3=0.0d0 + eello5_4=0.0d0 +cd call checkint5(i,j,k,l,jj,kk,eel5_1_num,eel5_2_num, +cd & eel5_3_num,eel5_4_num) + do iii=1,2 + do kkk=1,5 + do lll=1,3 + derx(lll,kkk,iii)=0.0d0 + enddo + enddo + enddo +cd eij=facont_hb(jj,i) +cd ekl=facont_hb(kk,k) +cd ekont=eij*ekl +cd write (iout,*)'Contacts have occurred for peptide groups', +cd & i,j,' fcont:',eij,' eij',' and ',k,l +cd goto 1111 +C Contribution from the graph I. +cd write (2,*) 'AEA ',AEA(1,1,1),AEA(2,1,1),AEA(1,2,1),AEA(2,2,1) +cd write (2,*) 'AEAb2',AEAb2(1,1,1),AEAb2(2,1,1) + call transpose2(EUg(1,1,k),auxmat(1,1)) + call matmat2(AEA(1,1,1),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + eello5_1=scalar2(AEAb2(1,1,1),Ub2(1,k)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,i)) + if (calc_grad) then +C Explicit gradient in virtual-dihedral angles. + if (i.gt.1) g_corr5_loc(i-1)=g_corr5_loc(i-1) + & +ekont*(scalar2(AEAb2derg(1,2,1,1),Ub2(1,k)) + & +0.5d0*scalar2(vv(1),Dtobr2der(1,i))) + call transpose2(EUgder(1,1,k),auxmat1(1,1)) + call matmat2(AEA(1,1,1),auxmat1(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + g_corr5_loc(k-1)=g_corr5_loc(k-1) + & +ekont*(scalar2(AEAb2(1,1,1),Ub2der(1,k)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,i))) + call matmat2(AEAderg(1,1,1),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + if (l.eq.j+1) then + if (l.lt.nres-1) g_corr5_loc(l-1)=g_corr5_loc(l-1) + & +ekont*(scalar2(AEAb2derg(1,1,1,1),Ub2(1,k)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,i))) + else + if (j.lt.nres-1) g_corr5_loc(j-1)=g_corr5_loc(j-1) + & +ekont*(scalar2(AEAb2derg(1,1,1,1),Ub2(1,k)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,i))) + endif +C Cartesian gradient + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call matmat2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1), + & pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + derx(lll,kkk,iii)=derx(lll,kkk,iii) + & +scalar2(AEAb2derx(1,lll,kkk,iii,1,1),Ub2(1,k)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,i)) + enddo + enddo + enddo +c goto 1112 + endif +c1111 continue +C Contribution from graph II + call transpose2(EE(1,1,itk),auxmat(1,1)) + call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + eello5_2=scalar2(AEAb1(1,2,1),b1(1,itk)) + & -0.5d0*scalar2(vv(1),Ctobr(1,k)) + if (calc_grad) then +C Explicit gradient in virtual-dihedral angles. + g_corr5_loc(k-1)=g_corr5_loc(k-1) + & -0.5d0*ekont*scalar2(vv(1),Ctobrder(1,k)) + call matmat2(auxmat(1,1),AEAderg(1,1,1),pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + if (l.eq.j+1) then + g_corr5_loc(l-1)=g_corr5_loc(l-1) + & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk)) + & -0.5d0*scalar2(vv(1),Ctobr(1,k))) + else + g_corr5_loc(j-1)=g_corr5_loc(j-1) + & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk)) + & -0.5d0*scalar2(vv(1),Ctobr(1,k))) + endif +C Cartesian gradient + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,1), + & pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + derx(lll,kkk,iii)=derx(lll,kkk,iii) + & +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,itk)) + & -0.5d0*scalar2(vv(1),Ctobr(1,k)) + enddo + enddo + enddo +cd goto 1112 + endif +cd1111 continue + if (l.eq.j+1) then +cd goto 1110 +C Parallel orientation +C Contribution from graph III + call transpose2(EUg(1,1,l),auxmat(1,1)) + call matmat2(AEA(1,1,2),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + eello5_3=scalar2(AEAb2(1,1,2),Ub2(1,l)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,j)) + if (calc_grad) then +C Explicit gradient in virtual-dihedral angles. + g_corr5_loc(j-1)=g_corr5_loc(j-1) + & +ekont*(scalar2(AEAb2derg(1,2,1,2),Ub2(1,l)) + & +0.5d0*scalar2(vv(1),Dtobr2der(1,j))) + call matmat2(AEAderg(1,1,2),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + g_corr5_loc(k-1)=g_corr5_loc(k-1) + & +ekont*(scalar2(AEAb2derg(1,1,1,2),Ub2(1,l)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,j))) + call transpose2(EUgder(1,1,l),auxmat1(1,1)) + call matmat2(AEA(1,1,2),auxmat1(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + g_corr5_loc(l-1)=g_corr5_loc(l-1) + & +ekont*(scalar2(AEAb2(1,1,2),Ub2der(1,l)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,j))) +C Cartesian gradient + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call matmat2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1), + & pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + derx(lll,kkk,iii)=derx(lll,kkk,iii) + & +scalar2(AEAb2derx(1,lll,kkk,iii,1,2),Ub2(1,l)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,j)) + enddo + enddo + enddo +cd goto 1112 + endif +C Contribution from graph IV +cd1110 continue + call transpose2(EE(1,1,itl),auxmat(1,1)) + call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + eello5_4=scalar2(AEAb1(1,2,2),b1(1,itl)) + & -0.5d0*scalar2(vv(1),Ctobr(1,l)) + if (calc_grad) then +C Explicit gradient in virtual-dihedral angles. + g_corr5_loc(l-1)=g_corr5_loc(l-1) + & -0.5d0*ekont*scalar2(vv(1),Ctobrder(1,l)) + call matmat2(auxmat(1,1),AEAderg(1,1,2),pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + g_corr5_loc(k-1)=g_corr5_loc(k-1) + & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itl)) + & -0.5d0*scalar2(vv(1),Ctobr(1,l))) +C Cartesian gradient + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,2), + & pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + derx(lll,kkk,iii)=derx(lll,kkk,iii) + & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itl)) + & -0.5d0*scalar2(vv(1),Ctobr(1,l)) + enddo + enddo + enddo + endif + else +C Antiparallel orientation +C Contribution from graph III +c goto 1110 + call transpose2(EUg(1,1,j),auxmat(1,1)) + call matmat2(AEA(1,1,2),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + eello5_3=scalar2(AEAb2(1,1,2),Ub2(1,j)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,l)) + if (calc_grad) then +C Explicit gradient in virtual-dihedral angles. + g_corr5_loc(l-1)=g_corr5_loc(l-1) + & +ekont*(scalar2(AEAb2derg(1,2,1,2),Ub2(1,j)) + & +0.5d0*scalar2(vv(1),Dtobr2der(1,l))) + call matmat2(AEAderg(1,1,2),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + g_corr5_loc(k-1)=g_corr5_loc(k-1) + & +ekont*(scalar2(AEAb2derg(1,1,1,2),Ub2(1,j)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,l))) + call transpose2(EUgder(1,1,j),auxmat1(1,1)) + call matmat2(AEA(1,1,2),auxmat1(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + g_corr5_loc(j-1)=g_corr5_loc(j-1) + & +ekont*(scalar2(AEAb2(1,1,2),Ub2der(1,j)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,l))) +C Cartesian gradient + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call matmat2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1), + & pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii) + & +scalar2(AEAb2derx(1,lll,kkk,iii,1,2),Ub2(1,j)) + & +0.5d0*scalar2(vv(1),Dtobr2(1,l)) + enddo + enddo + enddo +cd goto 1112 + endif +C Contribution from graph IV +1110 continue + call transpose2(EE(1,1,itj),auxmat(1,1)) + call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + eello5_4=scalar2(AEAb1(1,2,2),b1(1,itj)) + & -0.5d0*scalar2(vv(1),Ctobr(1,j)) + if (calc_grad) then +C Explicit gradient in virtual-dihedral angles. + g_corr5_loc(j-1)=g_corr5_loc(j-1) + & -0.5d0*ekont*scalar2(vv(1),Ctobrder(1,j)) + call matmat2(auxmat(1,1),AEAderg(1,1,2),pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + g_corr5_loc(k-1)=g_corr5_loc(k-1) + & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itj)) + & -0.5d0*scalar2(vv(1),Ctobr(1,j))) +C Cartesian gradient + do iii=1,2 + do kkk=1,5 + do lll=1,3 + call matmat2(auxmat(1,1),AEAderx(1,1,lll,kkk,iii,2), + & pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii) + & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itj)) + & -0.5d0*scalar2(vv(1),Ctobr(1,j)) + enddo + enddo + enddo + endif + endif +1112 continue + eel5=eello5_1+eello5_2+eello5_3+eello5_4 +cd if (i.eq.2 .and. j.eq.8 .and. k.eq.3 .and. l.eq.7) then +cd write (2,*) 'ijkl',i,j,k,l +cd write (2,*) 'eello5_1',eello5_1,' eello5_2',eello5_2, +cd & ' eello5_3',eello5_3,' eello5_4',eello5_4 +cd endif +cd write(iout,*) 'eello5_1',eello5_1,' eel5_1_num',16*eel5_1_num +cd write(iout,*) 'eello5_2',eello5_2,' eel5_2_num',16*eel5_2_num +cd write(iout,*) 'eello5_3',eello5_3,' eel5_3_num',16*eel5_3_num +cd write(iout,*) 'eello5_4',eello5_4,' eel5_4_num',16*eel5_4_num + if (calc_grad) then + if (j.lt.nres-1) then + j1=j+1 + j2=j-1 + else + j1=j-1 + j2=j-2 + endif + if (l.lt.nres-1) then + l1=l+1 + l2=l-1 + else + l1=l-1 + l2=l-2 + endif +cd eij=1.0d0 +cd ekl=1.0d0 +cd ekont=1.0d0 +cd write (2,*) 'eij',eij,' ekl',ekl,' ekont',ekont + do ll=1,3 + ggg1(ll)=eel5*g_contij(ll,1) + ggg2(ll)=eel5*g_contij(ll,2) +cold ghalf=0.5d0*eel5*ekl*gacont_hbr(ll,jj,i) + ghalf=0.5d0*ggg1(ll) +cd ghalf=0.0d0 + gradcorr5(ll,i)=gradcorr5(ll,i)+ghalf+ekont*derx(ll,2,1) + gradcorr5(ll,i+1)=gradcorr5(ll,i+1)+ekont*derx(ll,3,1) + gradcorr5(ll,j)=gradcorr5(ll,j)+ghalf+ekont*derx(ll,4,1) + gradcorr5(ll,j1)=gradcorr5(ll,j1)+ekont*derx(ll,5,1) +cold ghalf=0.5d0*eel5*eij*gacont_hbr(ll,kk,k) + ghalf=0.5d0*ggg2(ll) +cd ghalf=0.0d0 + gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2) + gradcorr5(ll,k+1)=gradcorr5(ll,k+1)+ekont*derx(ll,3,2) + gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2) + gradcorr5(ll,l1)=gradcorr5(ll,l1)+ekont*derx(ll,5,2) + enddo +cd goto 1112 + do m=i+1,j-1 + do ll=1,3 +cold gradcorr5(ll,m)=gradcorr5(ll,m)+eel5*ekl*gacont_hbr(ll,jj,i) + gradcorr5(ll,m)=gradcorr5(ll,m)+ggg1(ll) + enddo + enddo + do m=k+1,l-1 + do ll=1,3 +cold gradcorr5(ll,m)=gradcorr5(ll,m)+eel5*eij*gacont_hbr(ll,kk,k) + gradcorr5(ll,m)=gradcorr5(ll,m)+ggg2(ll) + enddo + enddo +c1112 continue + do m=i+2,j2 + do ll=1,3 + gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,1) + enddo + enddo + do m=k+2,l2 + do ll=1,3 + gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,2) + enddo + enddo +cd do iii=1,nres-3 +cd write (2,*) iii,g_corr5_loc(iii) +cd enddo + endif + eello5=ekont*eel5 +cd write (2,*) 'ekont',ekont +cd write (iout,*) 'eello5',ekont*eel5 + return + end +c-------------------------------------------------------------------------- + double precision function eello6(i,j,k,l,jj,kk) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + include 'COMMON.TORSION' + include 'COMMON.VAR' + include 'COMMON.GEO' + include 'COMMON.FFIELD' + double precision ggg1(3),ggg2(3) +cd if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then +cd eello6=0.0d0 +cd return +cd endif +cd write (iout,*) +cd & 'EELLO6: Contacts have occurred for peptide groups',i,j, +cd & ' and',k,l + eello6_1=0.0d0 + eello6_2=0.0d0 + eello6_3=0.0d0 + eello6_4=0.0d0 + eello6_5=0.0d0 + eello6_6=0.0d0 +cd call checkint6(i,j,k,l,jj,kk,eel6_1_num,eel6_2_num, +cd & eel6_3_num,eel6_4_num,eel6_5_num,eel6_6_num) + do iii=1,2 + do kkk=1,5 + do lll=1,3 + derx(lll,kkk,iii)=0.0d0 + enddo + enddo + enddo +cd eij=facont_hb(jj,i) +cd ekl=facont_hb(kk,k) +cd ekont=eij*ekl +cd eij=1.0d0 +cd ekl=1.0d0 +cd ekont=1.0d0 + if (l.eq.j+1) then + eello6_1=eello6_graph1(i,j,k,l,1,.false.) + eello6_2=eello6_graph1(j,i,l,k,2,.false.) + eello6_3=eello6_graph2(i,j,k,l,jj,kk,.false.) + eello6_4=eello6_graph4(i,j,k,l,jj,kk,1,.false.) + eello6_5=eello6_graph4(j,i,l,k,jj,kk,2,.false.) + eello6_6=eello6_graph3(i,j,k,l,jj,kk,.false.) + else + eello6_1=eello6_graph1(i,j,k,l,1,.false.) + eello6_2=eello6_graph1(l,k,j,i,2,.true.) + eello6_3=eello6_graph2(i,l,k,j,jj,kk,.true.) + eello6_4=eello6_graph4(i,j,k,l,jj,kk,1,.false.) + if (wturn6.eq.0.0d0 .or. j.ne.i+4) then + eello6_5=eello6_graph4(l,k,j,i,kk,jj,2,.true.) + else + eello6_5=0.0d0 + endif + eello6_6=eello6_graph3(i,l,k,j,jj,kk,.true.) + endif +C If turn contributions are considered, they will be handled separately. + eel6=eello6_1+eello6_2+eello6_3+eello6_4+eello6_5+eello6_6 +cd write(iout,*) 'eello6_1',eello6_1,' eel6_1_num',16*eel6_1_num +cd write(iout,*) 'eello6_2',eello6_2,' eel6_2_num',16*eel6_2_num +cd write(iout,*) 'eello6_3',eello6_3,' eel6_3_num',16*eel6_3_num +cd write(iout,*) 'eello6_4',eello6_4,' eel6_4_num',16*eel6_4_num +cd write(iout,*) 'eello6_5',eello6_5,' eel6_5_num',16*eel6_5_num +cd write(iout,*) 'eello6_6',eello6_6,' eel6_6_num',16*eel6_6_num +cd goto 1112 + if (calc_grad) then + if (j.lt.nres-1) then + j1=j+1 + j2=j-1 + else + j1=j-1 + j2=j-2 + endif + if (l.lt.nres-1) then + l1=l+1 + l2=l-1 + else + l1=l-1 + l2=l-2 + endif + do ll=1,3 + ggg1(ll)=eel6*g_contij(ll,1) + ggg2(ll)=eel6*g_contij(ll,2) +cold ghalf=0.5d0*eel6*ekl*gacont_hbr(ll,jj,i) + ghalf=0.5d0*ggg1(ll) +cd ghalf=0.0d0 + gradcorr6(ll,i)=gradcorr6(ll,i)+ghalf+ekont*derx(ll,2,1) + gradcorr6(ll,i+1)=gradcorr6(ll,i+1)+ekont*derx(ll,3,1) + gradcorr6(ll,j)=gradcorr6(ll,j)+ghalf+ekont*derx(ll,4,1) + gradcorr6(ll,j1)=gradcorr6(ll,j1)+ekont*derx(ll,5,1) + ghalf=0.5d0*ggg2(ll) +cold ghalf=0.5d0*eel6*eij*gacont_hbr(ll,kk,k) +cd ghalf=0.0d0 + gradcorr6(ll,k)=gradcorr6(ll,k)+ghalf+ekont*derx(ll,2,2) + gradcorr6(ll,k+1)=gradcorr6(ll,k+1)+ekont*derx(ll,3,2) + gradcorr6(ll,l)=gradcorr6(ll,l)+ghalf+ekont*derx(ll,4,2) + gradcorr6(ll,l1)=gradcorr6(ll,l1)+ekont*derx(ll,5,2) + enddo +cd goto 1112 + do m=i+1,j-1 + do ll=1,3 +cold gradcorr6(ll,m)=gradcorr6(ll,m)+eel6*ekl*gacont_hbr(ll,jj,i) + gradcorr6(ll,m)=gradcorr6(ll,m)+ggg1(ll) + enddo + enddo + do m=k+1,l-1 + do ll=1,3 +cold gradcorr6(ll,m)=gradcorr6(ll,m)+eel6*eij*gacont_hbr(ll,kk,k) + gradcorr6(ll,m)=gradcorr6(ll,m)+ggg2(ll) + enddo + enddo +1112 continue + do m=i+2,j2 + do ll=1,3 + gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,1) + enddo + enddo + do m=k+2,l2 + do ll=1,3 + gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,2) + enddo + enddo +cd do iii=1,nres-3 +cd write (2,*) iii,g_corr6_loc(iii) +cd enddo + endif + eello6=ekont*eel6 +cd write (2,*) 'ekont',ekont +cd write (iout,*) 'eello6',ekont*eel6 + return + end +c-------------------------------------------------------------------------- + double precision function eello6_graph1(i,j,k,l,imat,swap) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + include 'COMMON.TORSION' + include 'COMMON.VAR' + include 'COMMON.GEO' + double precision vv(2),vv1(2),pizda(2,2),auxmat(2,2),pizda1(2,2) + logical swap + logical lprn + common /kutas/ lprn +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C +C Parallel Antiparallel +C +C o o +C /l\ /j\ +C / \ / \ +C /| o | | o |\ +C \ j|/k\| / \ |/k\|l / +C \ / \ / \ / \ / +C o o o o +C i i +C +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + itk=itortyp(itype(k)) + s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i)) + s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k)) + s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k)) + call transpose2(EUgC(1,1,k),auxmat(1,1)) + call matmat2(AEA(1,1,imat),auxmat(1,1),pizda1(1,1)) + vv1(1)=pizda1(1,1)-pizda1(2,2) + vv1(2)=pizda1(1,2)+pizda1(2,1) + s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i)) + vv(1)=AEAb1(1,2,imat)*b1(1,itk)-AEAb1(2,2,imat)*b1(2,itk) + vv(2)=AEAb1(1,2,imat)*b1(2,itk)+AEAb1(2,2,imat)*b1(1,itk) + s5=scalar2(vv(1),Dtobr2(1,i)) +cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5 + eello6_graph1=-0.5d0*(s1+s2+s3+s4+s5) + if (.not. calc_grad) return + if (i.gt.1) g_corr6_loc(i-1)=g_corr6_loc(i-1) + & -0.5d0*ekont*(scalar2(AEAb1(1,2,imat),CUgb2der(1,i)) + & -scalar2(AEAb2derg(1,2,1,imat),Ug2Db1t(1,k)) + & +scalar2(AEAb2derg(1,2,1,imat),CUgb2(1,k)) + & +0.5d0*scalar2(vv1(1),Dtobr2der(1,i)) + & +scalar2(vv(1),Dtobr2der(1,i))) + call matmat2(AEAderg(1,1,imat),auxmat(1,1),pizda1(1,1)) + vv1(1)=pizda1(1,1)-pizda1(2,2) + vv1(2)=pizda1(1,2)+pizda1(2,1) + vv(1)=AEAb1derg(1,2,imat)*b1(1,itk)-AEAb1derg(2,2,imat)*b1(2,itk) + vv(2)=AEAb1derg(1,2,imat)*b1(2,itk)+AEAb1derg(2,2,imat)*b1(1,itk) + if (l.eq.j+1) then + g_corr6_loc(l-1)=g_corr6_loc(l-1) + & +ekont*(-0.5d0*(scalar2(AEAb1derg(1,2,imat),CUgb2(1,i)) + & -scalar2(AEAb2derg(1,1,1,imat),Ug2Db1t(1,k)) + & +scalar2(AEAb2derg(1,1,1,imat),CUgb2(1,k)) + & +0.5d0*scalar2(vv1(1),Dtobr2(1,i))+scalar2(vv(1),Dtobr2(1,i)))) + else + g_corr6_loc(j-1)=g_corr6_loc(j-1) + & +ekont*(-0.5d0*(scalar2(AEAb1derg(1,2,imat),CUgb2(1,i)) + & -scalar2(AEAb2derg(1,1,1,imat),Ug2Db1t(1,k)) + & +scalar2(AEAb2derg(1,1,1,imat),CUgb2(1,k)) + & +0.5d0*scalar2(vv1(1),Dtobr2(1,i))+scalar2(vv(1),Dtobr2(1,i)))) + endif + call transpose2(EUgCder(1,1,k),auxmat(1,1)) + call matmat2(AEA(1,1,imat),auxmat(1,1),pizda1(1,1)) + vv1(1)=pizda1(1,1)-pizda1(2,2) + vv1(2)=pizda1(1,2)+pizda1(2,1) + if (k.gt.1) g_corr6_loc(k-1)=g_corr6_loc(k-1) + & +ekont*(-0.5d0*(-scalar2(AEAb2(1,1,imat),Ug2Db1tder(1,k)) + & +scalar2(AEAb2(1,1,imat),CUgb2der(1,k)) + & +0.5d0*scalar2(vv1(1),Dtobr2(1,i)))) + do iii=1,2 + if (swap) then + ind=3-iii + else + ind=iii + endif + do kkk=1,5 + do lll=1,3 + s1= scalar2(AEAb1derx(1,lll,kkk,iii,2,imat),CUgb2(1,i)) + s2=-scalar2(AEAb2derx(1,lll,kkk,iii,1,imat),Ug2Db1t(1,k)) + s3= scalar2(AEAb2derx(1,lll,kkk,iii,1,imat),CUgb2(1,k)) + call transpose2(EUgC(1,1,k),auxmat(1,1)) + call matmat2(AEAderx(1,1,lll,kkk,iii,imat),auxmat(1,1), + & pizda1(1,1)) + vv1(1)=pizda1(1,1)-pizda1(2,2) + vv1(2)=pizda1(1,2)+pizda1(2,1) + s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i)) + vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,itk) + & -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,itk) + vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,itk) + & +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,itk) + s5=scalar2(vv(1),Dtobr2(1,i)) + derx(lll,kkk,ind)=derx(lll,kkk,ind)-0.5d0*(s1+s2+s3+s4+s5) + enddo + enddo + enddo + return + end +c---------------------------------------------------------------------------- + double precision function eello6_graph2(i,j,k,l,jj,kk,swap) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + include 'COMMON.TORSION' + include 'COMMON.VAR' + include 'COMMON.GEO' + logical swap + double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2), + & auxvec1(2),auxvec2(1),auxmat1(2,2) + logical lprn + common /kutas/ lprn +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C +C Parallel Antiparallel +C +C o o +C \ /l\ /j\ / +C \ / \ / \ / +C o| o | | o |o +C \ j|/k\| \ |/k\|l +C \ / \ \ / \ +C o o +C i i +C +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +cd write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l +C AL 7/4/01 s1 would occur in the sixth-order moment, +C but not in a cluster cumulant +#ifdef MOMENT + s1=dip(1,jj,i)*dip(1,kk,k) +#endif + call matvec2(ADtEA1(1,1,1),Ub2(1,k),auxvec(1)) + s2=-0.5d0*scalar2(Ub2(1,i),auxvec(1)) + call matvec2(ADtEA(1,1,2),Ub2(1,l),auxvec1(1)) + s3=-0.5d0*scalar2(Ub2(1,j),auxvec1(1)) + call transpose2(EUg(1,1,k),auxmat(1,1)) + call matmat2(ADtEA1(1,1,1),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + s4=-0.25d0*scalar2(vv(1),Dtobr2(1,i)) +cd write (2,*) 'eello6_graph2:','s1',s1,' s2',s2,' s3',s3,' s4',s4 +#ifdef MOMENT + eello6_graph2=-(s1+s2+s3+s4) +#else + eello6_graph2=-(s2+s3+s4) +#endif +c eello6_graph2=-s3 + if (.not. calc_grad) return +C Derivatives in gamma(i-1) + if (i.gt.1) then +#ifdef MOMENT + s1=dipderg(1,jj,i)*dip(1,kk,k) +#endif + s2=-0.5d0*scalar2(Ub2der(1,i),auxvec(1)) + call matvec2(ADtEAderg(1,1,1,2),Ub2(1,l),auxvec2(1)) + s3=-0.5d0*scalar2(Ub2(1,j),auxvec2(1)) + s4=-0.25d0*scalar2(vv(1),Dtobr2der(1,i)) +#ifdef MOMENT + g_corr6_loc(i-1)=g_corr6_loc(i-1)-ekont*(s1+s2+s3+s4) +#else + g_corr6_loc(i-1)=g_corr6_loc(i-1)-ekont*(s2+s3+s4) +#endif +c g_corr6_loc(i-1)=g_corr6_loc(i-1)-s3 + endif +C Derivatives in gamma(k-1) +#ifdef MOMENT + s1=dip(1,jj,i)*dipderg(1,kk,k) +#endif + call matvec2(ADtEA1(1,1,1),Ub2der(1,k),auxvec2(1)) + s2=-0.5d0*scalar2(Ub2(1,i),auxvec2(1)) + call matvec2(ADtEAderg(1,1,2,2),Ub2(1,l),auxvec2(1)) + s3=-0.5d0*scalar2(Ub2(1,j),auxvec2(1)) + call transpose2(EUgder(1,1,k),auxmat1(1,1)) + call matmat2(ADtEA1(1,1,1),auxmat1(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + s4=-0.25d0*scalar2(vv(1),Dtobr2(1,i)) +#ifdef MOMENT + g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s1+s2+s3+s4) +#else + g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s2+s3+s4) +#endif +c g_corr6_loc(k-1)=g_corr6_loc(k-1)-s3 +C Derivatives in gamma(j-1) or gamma(l-1) + if (j.gt.1) then +#ifdef MOMENT + s1=dipderg(3,jj,i)*dip(1,kk,k) +#endif + call matvec2(ADtEA1derg(1,1,1,1),Ub2(1,k),auxvec2(1)) + s2=-0.5d0*scalar2(Ub2(1,i),auxvec2(1)) + s3=-0.5d0*scalar2(Ub2der(1,j),auxvec1(1)) + call matmat2(ADtEA1derg(1,1,1,1),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + s4=-0.25d0*scalar2(vv(1),Dtobr2(1,i)) +#ifdef MOMENT + if (swap) then + g_corr6_loc(l-1)=g_corr6_loc(l-1)-ekont*s1 + else + g_corr6_loc(j-1)=g_corr6_loc(j-1)-ekont*s1 + endif +#endif + g_corr6_loc(j-1)=g_corr6_loc(j-1)-ekont*(s2+s3+s4) +c g_corr6_loc(j-1)=g_corr6_loc(j-1)-s3 + endif +C Derivatives in gamma(l-1) or gamma(j-1) + if (l.gt.1) then +#ifdef MOMENT + s1=dip(1,jj,i)*dipderg(3,kk,k) +#endif + call matvec2(ADtEA1derg(1,1,2,1),Ub2(1,k),auxvec2(1)) + s2=-0.5d0*scalar2(Ub2(1,i),auxvec2(1)) + call matvec2(ADtEA(1,1,2),Ub2der(1,l),auxvec2(1)) + s3=-0.5d0*scalar2(Ub2(1,j),auxvec2(1)) + call matmat2(ADtEA1derg(1,1,2,1),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + s4=-0.25d0*scalar2(vv(1),Dtobr2(1,i)) +#ifdef MOMENT + if (swap) then + g_corr6_loc(j-1)=g_corr6_loc(j-1)-ekont*s1 + else + g_corr6_loc(l-1)=g_corr6_loc(l-1)-ekont*s1 + endif +#endif + g_corr6_loc(l-1)=g_corr6_loc(l-1)-ekont*(s2+s3+s4) +c g_corr6_loc(l-1)=g_corr6_loc(l-1)-s3 + endif +C Cartesian derivatives. + if (lprn) then + write (2,*) 'In eello6_graph2' + do iii=1,2 + write (2,*) 'iii=',iii + do kkk=1,5 + write (2,*) 'kkk=',kkk + do jjj=1,2 + write (2,'(3(2f10.5),5x)') + & ((ADtEA1derx(jjj,mmm,lll,kkk,iii,1),mmm=1,2),lll=1,3) + enddo + enddo + enddo + endif + do iii=1,2 + do kkk=1,5 + do lll=1,3 +#ifdef MOMENT + if (iii.eq.1) then + s1=dipderx(lll,kkk,1,jj,i)*dip(1,kk,k) + else + s1=dip(1,jj,i)*dipderx(lll,kkk,1,kk,k) + endif +#endif + call matvec2(ADtEA1derx(1,1,lll,kkk,iii,1),Ub2(1,k), + & auxvec(1)) + s2=-0.5d0*scalar2(Ub2(1,i),auxvec(1)) + call matvec2(ADtEAderx(1,1,lll,kkk,iii,2),Ub2(1,l), + & auxvec(1)) + s3=-0.5d0*scalar2(Ub2(1,j),auxvec(1)) + call transpose2(EUg(1,1,k),auxmat(1,1)) + call matmat2(ADtEA1derx(1,1,lll,kkk,iii,1),auxmat(1,1), + & pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(1,2)+pizda(2,1) + s4=-0.25d0*scalar2(vv(1),Dtobr2(1,i)) +cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4',s4 +#ifdef MOMENT + derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s1+s2+s4) +#else + derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s2+s4) +#endif + if (swap) then + derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)-s3 + else + derx(lll,kkk,iii)=derx(lll,kkk,iii)-s3 + endif + enddo + enddo + enddo + return + end +c---------------------------------------------------------------------------- + double precision function eello6_graph3(i,j,k,l,jj,kk,swap) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + include 'COMMON.TORSION' + include 'COMMON.VAR' + include 'COMMON.GEO' + double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2) + logical swap +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C +C Parallel Antiparallel +C +C o o +C /l\ / \ /j\ +C / \ / \ / \ +C /| o |o o| o |\ +C j|/k\| / |/k\|l / +C / \ / / \ / +C / o / o +C i i +C +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C +C 4/7/01 AL Component s1 was removed, because it pertains to the respective +C energy moment and not to the cluster cumulant. + iti=itortyp(itype(i)) + if (j.lt.nres-1) then + itj1=itortyp(itype(j+1)) + else + itj1=ntortyp+1 + endif + itk=itortyp(itype(k)) + itk1=itortyp(itype(k+1)) + if (l.lt.nres-1) then + itl1=itortyp(itype(l+1)) + else + itl1=ntortyp+1 + endif +#ifdef MOMENT + s1=dip(4,jj,i)*dip(4,kk,k) +#endif + call matvec2(AECA(1,1,1),b1(1,itk1),auxvec(1)) + s2=0.5d0*scalar2(b1(1,itk),auxvec(1)) + call matvec2(AECA(1,1,2),b1(1,itl1),auxvec(1)) + s3=0.5d0*scalar2(b1(1,itj1),auxvec(1)) + call transpose2(EE(1,1,itk),auxmat(1,1)) + call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + s4=-0.25d0*scalar2(vv(1),Ctobr(1,k)) +cd write (2,*) 'eello6_graph3:','s1',s1,' s2',s2,' s3',s3,' s4',s4 +#ifdef MOMENT + eello6_graph3=-(s1+s2+s3+s4) +#else + eello6_graph3=-(s2+s3+s4) +#endif +c eello6_graph3=-s4 + if (.not. calc_grad) return +C Derivatives in gamma(k-1) + call matvec2(AECAderg(1,1,2),b1(1,itl1),auxvec(1)) + s3=0.5d0*scalar2(b1(1,itj1),auxvec(1)) + s4=-0.25d0*scalar2(vv(1),Ctobrder(1,k)) + g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s3+s4) +C Derivatives in gamma(l-1) + call matvec2(AECAderg(1,1,1),b1(1,itk1),auxvec(1)) + s2=0.5d0*scalar2(b1(1,itk),auxvec(1)) + call matmat2(auxmat(1,1),AECAderg(1,1,1),pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + s4=-0.25d0*scalar2(vv(1),Ctobr(1,k)) + g_corr6_loc(l-1)=g_corr6_loc(l-1)-ekont*(s2+s4) +C Cartesian derivatives. + do iii=1,2 + do kkk=1,5 + do lll=1,3 +#ifdef MOMENT + if (iii.eq.1) then + s1=dipderx(lll,kkk,4,jj,i)*dip(4,kk,k) + else + s1=dip(4,jj,i)*dipderx(lll,kkk,4,kk,k) + endif +#endif + call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,itk1), + & auxvec(1)) + s2=0.5d0*scalar2(b1(1,itk),auxvec(1)) + call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,itl1), + & auxvec(1)) + s3=0.5d0*scalar2(b1(1,itj1),auxvec(1)) + call matmat2(auxmat(1,1),AECAderx(1,1,lll,kkk,iii,1), + & pizda(1,1)) + vv(1)=pizda(1,1)+pizda(2,2) + vv(2)=pizda(2,1)-pizda(1,2) + s4=-0.25d0*scalar2(vv(1),Ctobr(1,k)) +#ifdef MOMENT + derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s1+s2+s4) +#else + derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s2+s4) +#endif + if (swap) then + derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)-s3 + else + derx(lll,kkk,iii)=derx(lll,kkk,iii)-s3 + endif +c derx(lll,kkk,iii)=derx(lll,kkk,iii)-s4 + enddo + enddo + enddo + return + end +c---------------------------------------------------------------------------- + double precision function eello6_graph4(i,j,k,l,jj,kk,imat,swap) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + include 'COMMON.TORSION' + include 'COMMON.VAR' + include 'COMMON.GEO' + include 'COMMON.FFIELD' + double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2), + & auxvec1(2),auxmat1(2,2) + logical swap +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C +C Parallel Antiparallel +C +C o o +C /l\ / \ /j\ +C / \ / \ / \ +C /| o |o o| o |\ +C \ j|/k\| \ |/k\|l +C \ / \ \ / \ +C o \ o \ +C i i +C +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C +C 4/7/01 AL Component s1 was removed, because it pertains to the respective +C energy moment and not to the cluster cumulant. +cd write (2,*) 'eello_graph4: wturn6',wturn6 + iti=itortyp(itype(i)) + itj=itortyp(itype(j)) + if (j.lt.nres-1) then + itj1=itortyp(itype(j+1)) + else + itj1=ntortyp+1 + endif + itk=itortyp(itype(k)) + if (k.lt.nres-1) then + itk1=itortyp(itype(k+1)) + else + itk1=ntortyp+1 + endif + itl=itortyp(itype(l)) + if (l.lt.nres-1) then + itl1=itortyp(itype(l+1)) + else + itl1=ntortyp+1 + endif +cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l +cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk, +cd & ' itl',itl,' itl1',itl1 +#ifdef MOMENT + if (imat.eq.1) then + s1=dip(3,jj,i)*dip(3,kk,k) + else + s1=dip(2,jj,j)*dip(2,kk,l) + endif +#endif + call matvec2(AECA(1,1,imat),Ub2(1,k),auxvec(1)) + s2=0.5d0*scalar2(Ub2(1,i),auxvec(1)) + if (j.eq.l+1) then + call matvec2(ADtEA1(1,1,3-imat),b1(1,itj1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1)) + else + call matvec2(ADtEA1(1,1,3-imat),b1(1,itl1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1)) + endif + call transpose2(EUg(1,1,k),auxmat(1,1)) + call matmat2(AECA(1,1,imat),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(2,1)+pizda(1,2) + s4=0.25d0*scalar2(vv(1),Dtobr2(1,i)) +cd write (2,*) 'eello6_graph4:','s1',s1,' s2',s2,' s3',s3,' s4',s4 +#ifdef MOMENT + eello6_graph4=-(s1+s2+s3+s4) +#else + eello6_graph4=-(s2+s3+s4) +#endif + if (.not. calc_grad) return +C Derivatives in gamma(i-1) + if (i.gt.1) then +#ifdef MOMENT + if (imat.eq.1) then + s1=dipderg(2,jj,i)*dip(3,kk,k) + else + s1=dipderg(4,jj,j)*dip(2,kk,l) + endif +#endif + s2=0.5d0*scalar2(Ub2der(1,i),auxvec(1)) + if (j.eq.l+1) then + call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itj1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1)) + else + call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itl1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1)) + endif + s4=0.25d0*scalar2(vv(1),Dtobr2der(1,i)) + if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then +cd write (2,*) 'turn6 derivatives' +#ifdef MOMENT + gel_loc_turn6(i-1)=gel_loc_turn6(i-1)-ekont*(s1+s2+s3+s4) +#else + gel_loc_turn6(i-1)=gel_loc_turn6(i-1)-ekont*(s2+s3+s4) +#endif + else +#ifdef MOMENT + g_corr6_loc(i-1)=g_corr6_loc(i-1)-ekont*(s1+s2+s3+s4) +#else + g_corr6_loc(i-1)=g_corr6_loc(i-1)-ekont*(s2+s3+s4) +#endif + endif + endif +C Derivatives in gamma(k-1) +#ifdef MOMENT + if (imat.eq.1) then + s1=dip(3,jj,i)*dipderg(2,kk,k) + else + s1=dip(2,jj,j)*dipderg(4,kk,l) + endif +#endif + call matvec2(AECA(1,1,imat),Ub2der(1,k),auxvec1(1)) + s2=0.5d0*scalar2(Ub2(1,i),auxvec1(1)) + if (j.eq.l+1) then + call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itj1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1)) + else + call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itl1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1)) + endif + call transpose2(EUgder(1,1,k),auxmat1(1,1)) + call matmat2(AECA(1,1,imat),auxmat1(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(2,1)+pizda(1,2) + s4=0.25d0*scalar2(vv(1),Dtobr2(1,i)) + if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then +#ifdef MOMENT + gel_loc_turn6(k-1)=gel_loc_turn6(k-1)-ekont*(s1+s2+s3+s4) +#else + gel_loc_turn6(k-1)=gel_loc_turn6(k-1)-ekont*(s2+s3+s4) +#endif + else +#ifdef MOMENT + g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s1+s2+s3+s4) +#else + g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s2+s3+s4) +#endif + endif +C Derivatives in gamma(j-1) or gamma(l-1) + if (l.eq.j+1 .and. l.gt.1) then + call matvec2(AECAderg(1,1,imat),Ub2(1,k),auxvec(1)) + s2=0.5d0*scalar2(Ub2(1,i),auxvec(1)) + call matmat2(AECAderg(1,1,imat),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(2,1)+pizda(1,2) + s4=0.25d0*scalar2(vv(1),Dtobr2(1,i)) + g_corr6_loc(l-1)=g_corr6_loc(l-1)-ekont*(s2+s4) + else if (j.gt.1) then + call matvec2(AECAderg(1,1,imat),Ub2(1,k),auxvec(1)) + s2=0.5d0*scalar2(Ub2(1,i),auxvec(1)) + call matmat2(AECAderg(1,1,imat),auxmat(1,1),pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(2,1)+pizda(1,2) + s4=0.25d0*scalar2(vv(1),Dtobr2(1,i)) + if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then + gel_loc_turn6(j-1)=gel_loc_turn6(j-1)-ekont*(s2+s4) + else + g_corr6_loc(j-1)=g_corr6_loc(j-1)-ekont*(s2+s4) + endif + endif +C Cartesian derivatives. + do iii=1,2 + do kkk=1,5 + do lll=1,3 +#ifdef MOMENT + if (iii.eq.1) then + if (imat.eq.1) then + s1=dipderx(lll,kkk,3,jj,i)*dip(3,kk,k) + else + s1=dipderx(lll,kkk,2,jj,j)*dip(2,kk,l) + endif + else + if (imat.eq.1) then + s1=dip(3,jj,i)*dipderx(lll,kkk,3,kk,k) + else + s1=dip(2,jj,j)*dipderx(lll,kkk,2,kk,l) + endif + endif +#endif + call matvec2(AECAderx(1,1,lll,kkk,iii,imat),Ub2(1,k), + & auxvec(1)) + s2=0.5d0*scalar2(Ub2(1,i),auxvec(1)) + if (j.eq.l+1) then + call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat), + & b1(1,itj1),auxvec(1)) + s3=-0.5d0*scalar2(b1(1,itj),auxvec(1)) + else + call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat), + & b1(1,itl1),auxvec(1)) + s3=-0.5d0*scalar2(b1(1,itl),auxvec(1)) + endif + call matmat2(AECAderx(1,1,lll,kkk,iii,imat),auxmat(1,1), + & pizda(1,1)) + vv(1)=pizda(1,1)-pizda(2,2) + vv(2)=pizda(2,1)+pizda(1,2) + s4=0.25d0*scalar2(vv(1),Dtobr2(1,i)) + if (swap) then + if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then +#ifdef MOMENT + derx_turn(lll,kkk,3-iii)=derx_turn(lll,kkk,3-iii) + & -(s1+s2+s4) +#else + derx_turn(lll,kkk,3-iii)=derx_turn(lll,kkk,3-iii) + & -(s2+s4) +#endif + derx_turn(lll,kkk,iii)=derx_turn(lll,kkk,iii)-s3 + else +#ifdef MOMENT + derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)-(s1+s2+s4) +#else + derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)-(s2+s4) +#endif + derx(lll,kkk,iii)=derx(lll,kkk,iii)-s3 + endif + else +#ifdef MOMENT + derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s1+s2+s4) +#else + derx(lll,kkk,iii)=derx(lll,kkk,iii)-(s2+s4) +#endif + if (l.eq.j+1) then + derx(lll,kkk,iii)=derx(lll,kkk,iii)-s3 + else + derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)-s3 + endif + endif + enddo + enddo + enddo + return + end +c---------------------------------------------------------------------------- + double precision function eello_turn6(i,jj,kk) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.INTERACT' + include 'COMMON.CONTACTS' + include 'COMMON.TORSION' + include 'COMMON.VAR' + include 'COMMON.GEO' + double precision vtemp1(2),vtemp2(2),vtemp3(2),vtemp4(2), + & atemp(2,2),auxmat(2,2),achuj_temp(2,2),gtemp(2,2),gvec(2), + & ggg1(3),ggg2(3) + double precision vtemp1d(2),vtemp2d(2),vtemp3d(2),vtemp4d(2), + & atempd(2,2),auxmatd(2,2),achuj_tempd(2,2),gtempd(2,2),gvecd(2) +C 4/7/01 AL Components s1, s8, and s13 were removed, because they pertain to +C the respective energy moment and not to the cluster cumulant. + eello_turn6=0.0d0 + j=i+4 + k=i+1 + l=i+3 + iti=itortyp(itype(i)) + itk=itortyp(itype(k)) + itk1=itortyp(itype(k+1)) + itl=itortyp(itype(l)) + itj=itortyp(itype(j)) +cd write (2,*) 'itk',itk,' itk1',itk1,' itl',itl,' itj',itj +cd write (2,*) 'i',i,' k',k,' j',j,' l',l +cd if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then +cd eello6=0.0d0 +cd return +cd endif +cd write (iout,*) +cd & 'EELLO6: Contacts have occurred for peptide groups',i,j, +cd & ' and',k,l +cd call checkint_turn6(i,jj,kk,eel_turn6_num) + do iii=1,2 + do kkk=1,5 + do lll=1,3 + derx_turn(lll,kkk,iii)=0.0d0 + enddo + enddo + enddo +cd eij=1.0d0 +cd ekl=1.0d0 +cd ekont=1.0d0 + eello6_5=eello6_graph4(l,k,j,i,kk,jj,2,.true.) +cd eello6_5=0.0d0 +cd write (2,*) 'eello6_5',eello6_5 +#ifdef MOMENT + call transpose2(AEA(1,1,1),auxmat(1,1)) + call matmat2(EUg(1,1,i+1),auxmat(1,1),auxmat(1,1)) + ss1=scalar2(Ub2(1,i+2),b1(1,itl)) + s1 = (auxmat(1,1)+auxmat(2,2))*ss1 +#endif + call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1)) + call matvec2(AEA(1,1,1),vtemp1(1),vtemp1(1)) + s2 = scalar2(b1(1,itk),vtemp1(1)) +#ifdef MOMENT + call transpose2(AEA(1,1,2),atemp(1,1)) + call matmat2(atemp(1,1),EUg(1,1,i+4),atemp(1,1)) + call matvec2(Ug2(1,1,i+2),dd(1,1,itk1),vtemp2(1)) + s8 = -(atemp(1,1)+atemp(2,2))*scalar2(cc(1,1,itl),vtemp2(1)) +#endif + call matmat2(EUg(1,1,i+3),AEA(1,1,2),auxmat(1,1)) + call matvec2(auxmat(1,1),Ub2(1,i+4),vtemp3(1)) + s12 = scalar2(Ub2(1,i+2),vtemp3(1)) +#ifdef MOMENT + call transpose2(a_chuj(1,1,kk,i+1),achuj_temp(1,1)) + call matmat2(achuj_temp(1,1),EUg(1,1,i+2),gtemp(1,1)) + call matmat2(gtemp(1,1),EUg(1,1,i+3),gtemp(1,1)) + call matvec2(a_chuj(1,1,jj,i),Ub2(1,i+4),vtemp4(1)) + ss13 = scalar2(b1(1,itk),vtemp4(1)) + s13 = (gtemp(1,1)+gtemp(2,2))*ss13 +#endif +c write (2,*) 's1,s2,s8,s12,s13',s1,s2,s8,s12,s13 +c s1=0.0d0 +c s2=0.0d0 +c s8=0.0d0 +c s12=0.0d0 +c s13=0.0d0 + eel_turn6 = eello6_5 - 0.5d0*(s1+s2+s12+s8+s13) + if (calc_grad) then +C Derivatives in gamma(i+2) +#ifdef MOMENT + call transpose2(AEA(1,1,1),auxmatd(1,1)) + call matmat2(EUgder(1,1,i+1),auxmatd(1,1),auxmatd(1,1)) + s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1 + call transpose2(AEAderg(1,1,2),atempd(1,1)) + call matmat2(atempd(1,1),EUg(1,1,i+4),atempd(1,1)) + s8d = -(atempd(1,1)+atempd(2,2))*scalar2(cc(1,1,itl),vtemp2(1)) +#endif + call matmat2(EUg(1,1,i+3),AEAderg(1,1,2),auxmatd(1,1)) + call matvec2(auxmatd(1,1),Ub2(1,i+4),vtemp3d(1)) + s12d = scalar2(Ub2(1,i+2),vtemp3d(1)) +c s1d=0.0d0 +c s2d=0.0d0 +c s8d=0.0d0 +c s12d=0.0d0 +c s13d=0.0d0 + gel_loc_turn6(i)=gel_loc_turn6(i)-0.5d0*ekont*(s1d+s8d+s12d) +C Derivatives in gamma(i+3) +#ifdef MOMENT + call transpose2(AEA(1,1,1),auxmatd(1,1)) + call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1)) + ss1d=scalar2(Ub2der(1,i+2),b1(1,itl)) + s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1d +#endif + call matvec2(EUgder(1,1,i+2),b1(1,itl),vtemp1d(1)) + call matvec2(AEA(1,1,1),vtemp1d(1),vtemp1d(1)) + s2d = scalar2(b1(1,itk),vtemp1d(1)) +#ifdef MOMENT + call matvec2(Ug2der(1,1,i+2),dd(1,1,itk1),vtemp2d(1)) + s8d = -(atemp(1,1)+atemp(2,2))*scalar2(cc(1,1,itl),vtemp2d(1)) +#endif + s12d = scalar2(Ub2der(1,i+2),vtemp3(1)) +#ifdef MOMENT + call matmat2(achuj_temp(1,1),EUgder(1,1,i+2),gtempd(1,1)) + call matmat2(gtempd(1,1),EUg(1,1,i+3),gtempd(1,1)) + s13d = (gtempd(1,1)+gtempd(2,2))*ss13 +#endif +c s1d=0.0d0 +c s2d=0.0d0 +c s8d=0.0d0 +c s12d=0.0d0 +c s13d=0.0d0 +#ifdef MOMENT + gel_loc_turn6(i+1)=gel_loc_turn6(i+1) + & -0.5d0*ekont*(s1d+s2d+s8d+s12d+s13d) +#else + gel_loc_turn6(i+1)=gel_loc_turn6(i+1) + & -0.5d0*ekont*(s2d+s12d) +#endif +C Derivatives in gamma(i+4) + call matmat2(EUgder(1,1,i+3),AEA(1,1,2),auxmatd(1,1)) + call matvec2(auxmatd(1,1),Ub2(1,i+4),vtemp3d(1)) + s12d = scalar2(Ub2(1,i+2),vtemp3d(1)) +#ifdef MOMENT + call matmat2(achuj_temp(1,1),EUg(1,1,i+2),gtempd(1,1)) + call matmat2(gtempd(1,1),EUgder(1,1,i+3),gtempd(1,1)) + s13d = (gtempd(1,1)+gtempd(2,2))*ss13 +#endif +c s1d=0.0d0 +c s2d=0.0d0 +c s8d=0.0d0 +C s12d=0.0d0 +c s13d=0.0d0 +#ifdef MOMENT + gel_loc_turn6(i+2)=gel_loc_turn6(i+2)-0.5d0*ekont*(s12d+s13d) +#else + gel_loc_turn6(i+2)=gel_loc_turn6(i+2)-0.5d0*ekont*(s12d) +#endif +C Derivatives in gamma(i+5) +#ifdef MOMENT + call transpose2(AEAderg(1,1,1),auxmatd(1,1)) + call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1)) + s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1 +#endif + call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1d(1)) + call matvec2(AEAderg(1,1,1),vtemp1d(1),vtemp1d(1)) + s2d = scalar2(b1(1,itk),vtemp1d(1)) +#ifdef MOMENT + call transpose2(AEA(1,1,2),atempd(1,1)) + call matmat2(atempd(1,1),EUgder(1,1,i+4),atempd(1,1)) + s8d = -(atempd(1,1)+atempd(2,2))*scalar2(cc(1,1,itl),vtemp2(1)) +#endif + call matvec2(auxmat(1,1),Ub2der(1,i+4),vtemp3d(1)) + s12d = scalar2(Ub2(1,i+2),vtemp3d(1)) +#ifdef MOMENT + call matvec2(a_chuj(1,1,jj,i),Ub2der(1,i+4),vtemp4d(1)) + ss13d = scalar2(b1(1,itk),vtemp4d(1)) + s13d = (gtemp(1,1)+gtemp(2,2))*ss13d +#endif +c s1d=0.0d0 +c s2d=0.0d0 +c s8d=0.0d0 +c s12d=0.0d0 +c s13d=0.0d0 +#ifdef MOMENT + gel_loc_turn6(i+3)=gel_loc_turn6(i+3) + & -0.5d0*ekont*(s1d+s2d+s8d+s12d+s13d) +#else + gel_loc_turn6(i+3)=gel_loc_turn6(i+3) + & -0.5d0*ekont*(s2d+s12d) +#endif +C Cartesian derivatives + do iii=1,2 + do kkk=1,5 + do lll=1,3 +#ifdef MOMENT + call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmatd(1,1)) + call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1)) + s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1 +#endif + call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1)) + call matvec2(AEAderx(1,1,lll,kkk,iii,1),vtemp1(1), + & vtemp1d(1)) + s2d = scalar2(b1(1,itk),vtemp1d(1)) +#ifdef MOMENT + call transpose2(AEAderx(1,1,lll,kkk,iii,2),atempd(1,1)) + call matmat2(atempd(1,1),EUg(1,1,i+4),atempd(1,1)) + s8d = -(atempd(1,1)+atempd(2,2))* + & scalar2(cc(1,1,itl),vtemp2(1)) +#endif + call matmat2(EUg(1,1,i+3),AEAderx(1,1,lll,kkk,iii,2), + & auxmatd(1,1)) + call matvec2(auxmatd(1,1),Ub2(1,i+4),vtemp3d(1)) + s12d = scalar2(Ub2(1,i+2),vtemp3d(1)) +c s1d=0.0d0 +c s2d=0.0d0 +c s8d=0.0d0 +c s12d=0.0d0 +c s13d=0.0d0 +#ifdef MOMENT + derx_turn(lll,kkk,iii) = derx_turn(lll,kkk,iii) + & - 0.5d0*(s1d+s2d) +#else + derx_turn(lll,kkk,iii) = derx_turn(lll,kkk,iii) + & - 0.5d0*s2d +#endif +#ifdef MOMENT + derx_turn(lll,kkk,3-iii) = derx_turn(lll,kkk,3-iii) + & - 0.5d0*(s8d+s12d) +#else + derx_turn(lll,kkk,3-iii) = derx_turn(lll,kkk,3-iii) + & - 0.5d0*s12d +#endif + enddo + enddo + enddo +#ifdef MOMENT + do kkk=1,5 + do lll=1,3 + call transpose2(a_chuj_der(1,1,lll,kkk,kk,i+1), + & achuj_tempd(1,1)) + call matmat2(achuj_tempd(1,1),EUg(1,1,i+2),gtempd(1,1)) + call matmat2(gtempd(1,1),EUg(1,1,i+3),gtempd(1,1)) + s13d=(gtempd(1,1)+gtempd(2,2))*ss13 + derx_turn(lll,kkk,2) = derx_turn(lll,kkk,2)-0.5d0*s13d + call matvec2(a_chuj_der(1,1,lll,kkk,jj,i),Ub2(1,i+4), + & vtemp4d(1)) + ss13d = scalar2(b1(1,itk),vtemp4d(1)) + s13d = (gtemp(1,1)+gtemp(2,2))*ss13d + derx_turn(lll,kkk,1) = derx_turn(lll,kkk,1)-0.5d0*s13d + enddo + enddo +#endif +cd write(iout,*) 'eel6_turn6',eel_turn6,' eel_turn6_num', +cd & 16*eel_turn6_num +cd goto 1112 + if (j.lt.nres-1) then + j1=j+1 + j2=j-1 + else + j1=j-1 + j2=j-2 + endif + if (l.lt.nres-1) then + l1=l+1 + l2=l-1 + else + l1=l-1 + l2=l-2 + endif + do ll=1,3 + ggg1(ll)=eel_turn6*g_contij(ll,1) + ggg2(ll)=eel_turn6*g_contij(ll,2) + ghalf=0.5d0*ggg1(ll) +cd ghalf=0.0d0 + gcorr6_turn(ll,i)=gcorr6_turn(ll,i)+ghalf + & +ekont*derx_turn(ll,2,1) + gcorr6_turn(ll,i+1)=gcorr6_turn(ll,i+1)+ekont*derx_turn(ll,3,1) + gcorr6_turn(ll,j)=gcorr6_turn(ll,j)+ghalf + & +ekont*derx_turn(ll,4,1) + gcorr6_turn(ll,j1)=gcorr6_turn(ll,j1)+ekont*derx_turn(ll,5,1) + ghalf=0.5d0*ggg2(ll) +cd ghalf=0.0d0 + gcorr6_turn(ll,k)=gcorr6_turn(ll,k)+ghalf + & +ekont*derx_turn(ll,2,2) + gcorr6_turn(ll,k+1)=gcorr6_turn(ll,k+1)+ekont*derx_turn(ll,3,2) + gcorr6_turn(ll,l)=gcorr6_turn(ll,l)+ghalf + & +ekont*derx_turn(ll,4,2) + gcorr6_turn(ll,l1)=gcorr6_turn(ll,l1)+ekont*derx_turn(ll,5,2) + enddo +cd goto 1112 + do m=i+1,j-1 + do ll=1,3 + gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg1(ll) + enddo + enddo + do m=k+1,l-1 + do ll=1,3 + gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg2(ll) + enddo + enddo +1112 continue + do m=i+2,j2 + do ll=1,3 + gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,1) + enddo + enddo + do m=k+2,l2 + do ll=1,3 + gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,2) + enddo + enddo +cd do iii=1,nres-3 +cd write (2,*) iii,g_corr6_loc(iii) +cd enddo + endif + eello_turn6=ekont*eel_turn6 +cd write (2,*) 'ekont',ekont +cd write (2,*) 'eel_turn6',ekont*eel_turn6 + return + end +crc------------------------------------------------- + SUBROUTINE MATVEC2(A1,V1,V2) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + DIMENSION A1(2,2),V1(2),V2(2) +c DO 1 I=1,2 +c VI=0.0 +c DO 3 K=1,2 +c 3 VI=VI+A1(I,K)*V1(K) +c Vaux(I)=VI +c 1 CONTINUE + + vaux1=a1(1,1)*v1(1)+a1(1,2)*v1(2) + vaux2=a1(2,1)*v1(1)+a1(2,2)*v1(2) + + v2(1)=vaux1 + v2(2)=vaux2 + END +C--------------------------------------- + SUBROUTINE MATMAT2(A1,A2,A3) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + DIMENSION A1(2,2),A2(2,2),A3(2,2) +c DIMENSION AI3(2,2) +c DO J=1,2 +c A3IJ=0.0 +c DO K=1,2 +c A3IJ=A3IJ+A1(I,K)*A2(K,J) +c enddo +c A3(I,J)=A3IJ +c enddo +c enddo + + ai3_11=a1(1,1)*a2(1,1)+a1(1,2)*a2(2,1) + ai3_12=a1(1,1)*a2(1,2)+a1(1,2)*a2(2,2) + ai3_21=a1(2,1)*a2(1,1)+a1(2,2)*a2(2,1) + ai3_22=a1(2,1)*a2(1,2)+a1(2,2)*a2(2,2) + + A3(1,1)=AI3_11 + A3(2,1)=AI3_21 + A3(1,2)=AI3_12 + A3(2,2)=AI3_22 + END + +c------------------------------------------------------------------------- + double precision function scalar2(u,v) + implicit none + double precision u(2),v(2) + double precision sc + integer i + scalar2=u(1)*v(1)+u(2)*v(2) + return + end + +C----------------------------------------------------------------------------- + + subroutine transpose2(a,at) + implicit none + double precision a(2,2),at(2,2) + at(1,1)=a(1,1) + at(1,2)=a(2,1) + at(2,1)=a(1,2) + at(2,2)=a(2,2) + return + end +c-------------------------------------------------------------------------- + subroutine transpose(n,a,at) + implicit none + integer n,i,j + double precision a(n,n),at(n,n) + do i=1,n + do j=1,n + at(j,i)=a(i,j) + enddo + enddo + return + end +C--------------------------------------------------------------------------- + subroutine prodmat3(a1,a2,kk,transp,prod) + implicit none + integer i,j + double precision a1(2,2),a2(2,2),a2t(2,2),kk(2,2),prod(2,2) + logical transp +crc double precision auxmat(2,2),prod_(2,2) + + if (transp) then +crc call transpose2(kk(1,1),auxmat(1,1)) +crc call matmat2(a1(1,1),auxmat(1,1),auxmat(1,1)) +crc call matmat2(auxmat(1,1),a2(1,1),prod_(1,1)) + + prod(1,1)=(a1(1,1)*kk(1,1)+a1(1,2)*kk(1,2))*a2(1,1) + & +(a1(1,1)*kk(2,1)+a1(1,2)*kk(2,2))*a2(2,1) + prod(1,2)=(a1(1,1)*kk(1,1)+a1(1,2)*kk(1,2))*a2(1,2) + & +(a1(1,1)*kk(2,1)+a1(1,2)*kk(2,2))*a2(2,2) + prod(2,1)=(a1(2,1)*kk(1,1)+a1(2,2)*kk(1,2))*a2(1,1) + & +(a1(2,1)*kk(2,1)+a1(2,2)*kk(2,2))*a2(2,1) + prod(2,2)=(a1(2,1)*kk(1,1)+a1(2,2)*kk(1,2))*a2(1,2) + & +(a1(2,1)*kk(2,1)+a1(2,2)*kk(2,2))*a2(2,2) + + else +crc call matmat2(a1(1,1),kk(1,1),auxmat(1,1)) +crc call matmat2(auxmat(1,1),a2(1,1),prod_(1,1)) + + prod(1,1)=(a1(1,1)*kk(1,1)+a1(1,2)*kk(2,1))*a2(1,1) + & +(a1(1,1)*kk(1,2)+a1(1,2)*kk(2,2))*a2(2,1) + prod(1,2)=(a1(1,1)*kk(1,1)+a1(1,2)*kk(2,1))*a2(1,2) + & +(a1(1,1)*kk(1,2)+a1(1,2)*kk(2,2))*a2(2,2) + prod(2,1)=(a1(2,1)*kk(1,1)+a1(2,2)*kk(2,1))*a2(1,1) + & +(a1(2,1)*kk(1,2)+a1(2,2)*kk(2,2))*a2(2,1) + prod(2,2)=(a1(2,1)*kk(1,1)+a1(2,2)*kk(2,1))*a2(1,2) + & +(a1(2,1)*kk(1,2)+a1(2,2)*kk(2,2))*a2(2,2) + + endif +c call transpose2(a2(1,1),a2t(1,1)) + +crc print *,transp +crc print *,((prod_(i,j),i=1,2),j=1,2) +crc print *,((prod(i,j),i=1,2),j=1,2) + + return + end +C----------------------------------------------------------------------------- + double precision function scalar(u,v) + implicit none + double precision u(3),v(3) + double precision sc + integer i + sc=0.0d0 + do i=1,3 + sc=sc+u(i)*v(i) + enddo + scalar=sc + return + end +