X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=7cf6721213722ee81e9fea5df3df839ffbad632b;hb=c1b8e24ad8ccf66d0617b0055089ea0c1193840d;hp=03a7825c39d96155059047f86eddf9cdbbf7c40c;hpb=c2b8becc759571af9f6dd36d48a2f7b16037925b;p=unres.git diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F index 03a7825..7cf6721 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -131,6 +131,11 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +cmc +cmc Sep-06: egb takes care of dynamic ss bonds too +cmc +c if (dyn_ss) call dyn_set_nss + c print *,"Processor",myrank," computed USCSC" #ifdef TIMING #ifdef MPI @@ -326,6 +331,7 @@ C energia(23)=evdw_m c print *," Processor",myrank," calls SUM_ENERGY" call sum_energy(energia,.true.) + if (dyn_ss) call dyn_set_nss c print *," Processor",myrank," left SUM_ENERGY" #ifdef TIMING #ifdef MPI @@ -423,14 +429,14 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 + & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 + & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor @@ -1052,25 +1058,25 @@ C------------------------------------------------------------------------ & edihcnstr,ebr*nss, & Uconst,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)'/ - & 'EVDWPP=',1pE16.6,' WEIGHT=',1pD16.6,' (p-p VDW)'/ - & 'ESTR= ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching)'/ - & '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)'/ - & 'EHPB= ',1pE16.6,' WEIGHT=',1pD16.6, + & 'EVDW= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-SC)'/ + & 'EVDW2= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-p)'/ + & 'EES= ',1pE16.6,' WEIGHT=',1pE16.6,' (p-p)'/ + & 'EVDWPP=',1pE16.6,' WEIGHT=',1pE16.6,' (p-p VDW)'/ + & 'ESTR= ',1pE16.6,' WEIGHT=',1pE16.6,' (stretching)'/ + & 'EBE= ',1pE16.6,' WEIGHT=',1pE16.6,' (bending)'/ + & 'ESC= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC local)'/ + & 'ETORS= ',1pE16.6,' WEIGHT=',1pE16.6,' (torsional)'/ + & 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/ + & 'EHPB= ',1pE16.6,' WEIGHT=',1pE16.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)'/ - & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ + & 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ + & 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ + & 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/ + & 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/ + & 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/ + & 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/ + & 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/ + & 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ @@ -1443,7 +1449,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) @@ -1552,6 +1558,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.CONTROL' + include 'COMMON.SBRIDGE' logical lprn evdw=0.0D0 ccccc energy_dec=.false. @@ -1580,6 +1587,27 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) + IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + call dyn_ssbond_ene(i,j,evdwij) + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') + & 'evdw',i,j,evdwij,' ss' +C triple bond artifac removal + do k=j+1,iend(i,iint) +C search over all next residues + if (dyn_ss_mask(k)) then +C check if they are cysteins +C write(iout,*) 'k=',k + call triple_ssbond_ene(i,j,k,evdwij) +C call the energy function that removes the artifical triple disulfide +C bond the soubroutine is located in ssMD.F + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') + & 'evdw',i,j,evdwij,'tss' + endif!dyn_ss_mask(k) + enddo! k + ELSE +C cycle ind=ind+1 itypj=itype(j) c dscj_inv=dsc_inv(itypj) @@ -1689,6 +1717,7 @@ C Calculate angular part of the gradient. #else call sc_grad #endif + ENDIF ! dyn_ss enddo ! j enddo ! iint enddo ! i @@ -4266,9 +4295,15 @@ c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then +cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then +C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds + if (.not.dyn_ss .and. i.le.nss) then +C 15/02/13 CC dynamic SSbond - additional check + if (ii.gt.nres + & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij + endif cd write (iout,*) "eij",eij else if (ii.gt.nres .and. jj.gt.nres) then c Restraints from contact prediction @@ -4407,7 +4442,7 @@ c dscj_inv=dsc_inv(itypj) deltat12=om2-om1+2.0d0 cosphi=om12-om1*om2 eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2) - & +akct*deltad*deltat12 + & +akct*deltad*deltat12+ebr & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi c write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth, c & " akct",akct," deltad",deltad," deltat",deltat1,deltat2, @@ -4549,6 +4584,18 @@ c write (*,'(a,i2)') 'EBEND ICG=',icg 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) + ichir1=isign(1,itype(i-2)) + ichir2=isign(1,itype(i)) + if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1)) + if (itype(i).eq.10) ichir2=isign(1,itype(i-1)) + if (itype(i-1).eq.10) then + itype1=isign(10,itype(i-2)) + ichir11=isign(1,itype(i-2)) + ichir12=isign(1,itype(i-2)) + itype2=isign(10,itype(i)) + ichir21=isign(1,itype(i)) + ichir22=isign(1,itype(i)) + endif if (i.gt.3) then #ifdef OSF phii=phi(i) @@ -4582,15 +4629,27 @@ 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) + athetk=athet(k,it,ichir1,ichir2) + bthetk=bthet(k,it,ichir1,ichir2) + if (it.eq.10) then + athetk=athet(k,itype1,ichir11,ichir12) + bthetk=bthet(k,itype2,ichir21,ichir22) + endif thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k) enddo dthett=thet_pred_mean*ssd thet_pred_mean=thet_pred_mean*ss+a0thet(it) 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 + dthetg1=(-athet(1,it,ichir1,ichir2)*y(2) + &+athet(2,it,ichir1,ichir2)*y(1))*ss + dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2) + & +bthet(2,it,ichir1,ichir2)*z(1))*ss + if (it.eq.10) then + dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2) + &+athet(2,itype1,ichir11,ichir12)*y(1))*ss + dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2) + & +bthet(2,itype2,ichir21,ichir22)*z(1))*ss + endif if (theta(i).gt.pi-delta) then call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0, & E_tc0) @@ -4788,6 +4847,9 @@ C enddo endif if (i.lt.nres) then + + if (iabs(itype(i+1)).eq.20) iblock=2 + if (iabs(itype(i+1)).ne.20) iblock=1 #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4808,7 +4870,7 @@ C sinph2(k)=0.0d0 enddo endif - ethetai=aa0thet(ityp1,ityp2,ityp3) + ethetai=aa0thet(ityp1,ityp2,ityp3,iblock) do k=1,ndouble do l=1,k-1 ccl=cosph1(l)*cosph2(k-l) @@ -4830,11 +4892,12 @@ C enddo endif do k=1,ntheterm - ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k) - dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3) + ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k) + dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock) & *coskt(k) if (lprn) - & write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3), + & write (iout,*) "k",k, + & "aathet",aathet(k,ityp1,ityp2,ityp3,iblock), & " ethetai",ethetai enddo if (lprn) then @@ -4853,24 +4916,24 @@ C endif do m=1,ntheterm2 do k=1,nsingle - aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k) - & +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k) - & +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k) - & +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k) + aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k) + & +ccthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k) + & +ddthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k) + & +eethet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k) ethetai=ethetai+sinkt(m)*aux dethetai=dethetai+0.5d0*m*aux*coskt(m) dephii=dephii+k*sinkt(m)*( - & ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)- - & bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)) + & ccthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)- + & bbthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k)) dephii1=dephii1+k*sinkt(m)*( - & eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)- - & ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k)) + & eethet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)- + & ddthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k)) if (lprn) & write (iout,*) "m",m," k",k," bbthet", - & bbthet(k,m,ityp1,ityp2,ityp3)," ccthet", - & ccthet(k,m,ityp1,ityp2,ityp3)," ddthet", - & ddthet(k,m,ityp1,ityp2,ityp3)," eethet", - & eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai + & bbthet(k,m,ityp1,ityp2,ityp3,iblock)," ccthet", + & ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet", + & ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet", + & eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai enddo enddo if (lprn) @@ -4878,28 +4941,33 @@ C do m=1,ntheterm3 do k=2,ndouble do l=1,k-1 - aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+ - & ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+ - & ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+ - & ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l) + aux=ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+ + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)+ + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+ + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l) + ethetai=ethetai+sinkt(m)*aux dethetai=dethetai+0.5d0*m*coskt(m)*aux dephii=dephii+l*sinkt(m)*( - & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)- - & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+ - & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+ - & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) + & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)- + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+ + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+ + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)) + dephii1=dephii1+(k-l)*sinkt(m)*( - & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+ - & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+ - & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)- - & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) + &-ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+ + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+ + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)- + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)) + if (lprn) then write (iout,*) "m",m," k",k," l",l," ffthet", - & ffthet(l,k,m,ityp1,ityp2,ityp3), - & ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet", - & ggthet(l,k,m,ityp1,ityp2,ityp3), - & ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai + & ffthet(l,k,m,ityp1,ityp2,ityp3,iblock), + & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ggthet", + & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock), + & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock), + & " ethetai",ethetai + write (iout,*) cosph1ph2(l,k)*sinkt(m), & cosph1ph2(k,l)*sinkt(m), & sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m) @@ -5251,7 +5319,7 @@ C cosfac=dsqrt(cosfac2) sinfac2=0.5d0/(1.0d0-costtab(i+1)) sinfac=dsqrt(sinfac2) - it=itype(i) + it=iabs(itype(i)) if (it.eq.10) goto 1 c C Compute the axes of tghe local cartesian coordinates system; store in @@ -5269,7 +5337,7 @@ C & dc_norm(3,i+nres) y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac enddo do j = 1,3 - z_prime(j) = -uz(j,i-1) + z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i))) enddo c write (2,*) "i",i c write (2,*) "x_prime",(x_prime(j),j=1,3) @@ -5301,7 +5369,7 @@ C C Compute the energy of the ith side cbain C c write (2,*) "xx",xx," yy",yy," zz",zz - it=itype(i) + it=iabs(itype(i)) do j = 1,65 x(j) = sc_parmin(j,it) enddo @@ -5309,7 +5377,7 @@ c write (2,*) "xx",xx," yy",yy," zz",zz Cc diagnostics - remove later xx1 = dcos(alph(2)) yy1 = dsin(alph(2))*dcos(omeg(2)) - zz1 = -dsin(alph(2))*dsin(omeg(2)) + zz1 = -dsign(1.0, dfloat(itype(i)))*dsin(alph(2))*dsin(omeg(2)) write(2,'(3f8.1,3f9.3,1x,3f9.3)') & alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz, & xx1,yy1,zz1 @@ -5479,8 +5547,10 @@ c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i) dZZ_Ci1(k)=0.0d0 dZZ_Ci(k)=0.0d0 do j=1,3 - dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres) - dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres) + dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1) + & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres) + dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1) + & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres) enddo dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres)) @@ -5765,15 +5835,22 @@ C Set lprn=.true. for debugging c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end - etors_ii=0.0D0 + if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 + & .or. itype(i).eq.ntyp1) cycle + etors_ii=0.0D0 + if (iabs(itype(i)).eq.20) then + iblock=2 + else + iblock=1 + endif 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) + do j=1,nterm(itori,itori1,iblock) + v1ij=v1(j,itori,itori1,iblock) + v2ij=v2(j,itori,itori1,iblock) cosphi=dcos(j*phii) sinphi=dsin(j*phii) etors=etors+v1ij*cosphi+v2ij*sinphi @@ -5788,7 +5865,7 @@ 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) + do j=1,nlor(itori,itori1,iblock) vl1ij=vlor1(j,itori,itori1) vl2ij=vlor2(j,itori,itori1) vl3ij=vlor3(j,itori,itori1) @@ -5801,13 +5878,14 @@ C gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom enddo C Subtract the constant term - etors=etors-v0(itori,itori1) + etors=etors-v0(itori,itori1,iblock) if (energy_dec) write (iout,'(a6,i5,0pf7.3)') - & 'etor',i,etors_ii-v0(itori,itori1) + & 'etor',i,etors_ii-v0(itori,itori1,iblock) 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) + & (v1(j,itori,itori1,iblock),j=1,6), + & (v2(j,itori,itori1,iblock),j=1,6) gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) enddo @@ -5858,7 +5936,10 @@ C Set lprn=.true. for debugging lprn=.false. c lprn=.true. etors_d=0.0D0 +c write(iout,*) "a tu??" do i=iphid_start,iphid_end + if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 + & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -5866,11 +5947,15 @@ c lprn=.true. phii1=phi(i+1) gloci1=0.0D0 gloci2=0.0D0 - 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) + iblock=1 + if (iabs(itype(i+1)).eq.20) iblock=2 + +C Regular cosine and sine terms + do j=1,ntermd_1(itori,itori1,itori2,iblock) + v1cij=v1c(1,j,itori,itori1,itori2,iblock) + v1sij=v1s(1,j,itori,itori1,itori2,iblock) + v2cij=v1c(2,j,itori,itori1,itori2,iblock) + v2sij=v1s(2,j,itori,itori1,itori2,iblock) cosphi1=dcos(j*phii) sinphi1=dsin(j*phii) cosphi2=dcos(j*phii1) @@ -5880,12 +5965,12 @@ c lprn=.true. gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1) gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2) enddo - do k=2,ntermd_2(itori,itori1,itori2) + do k=2,ntermd_2(itori,itori1,itori2,iblock) 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) + v1cdij = v2c(k,l,itori,itori1,itori2,iblock) + v2cdij = v2c(l,k,itori,itori1,itori2,iblock) + v1sdij = v2s(k,l,itori,itori1,itori2,iblock) + v2sdij = v2s(l,k,itori,itori1,itori2,iblock) cosphi1p2=dcos(l*phii+(k-l)*phii1) cosphi1m2=dcos(l*phii-(k-l)*phii1) sinphi1p2=dsin(l*phii+(k-l)*phii1) @@ -5900,7 +5985,6 @@ c lprn=.true. enddo gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1 gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2 -c write (iout,*) "gloci", gloc(i-3,icg) enddo return end @@ -5931,10 +6015,11 @@ c amino-acid residues. C Set lprn=.true. for debugging lprn=.false. c lprn=.true. -c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor +c write (iout,*) "EBACK_SC_COR",itau_start,itau_end esccor=0.0D0 do i=itau_start,itau_end esccor_ii=0.0D0 + if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle isccori=isccortyp(itype(i-2)) isccori1=isccortyp(itype(i-1)) phii=phi(i) @@ -5953,14 +6038,16 @@ c 2 = Ca...Ca...Ca...SC c 3 = SC...Ca...Ca...SCi gloci=0.0D0 if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. - & (itype(i-1).eq.10).or.(itype(i-2).eq.21).or. - & (itype(i-1).eq.21))) + & (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or. + & (itype(i-1).eq.ntyp1))) & .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) - & .or.(itype(i-2).eq.21))) + & .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1) + & .or.(itype(i).eq.ntyp1))) & .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. - & (itype(i-1).eq.21)))) cycle - if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle - if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21)) + & (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + & (itype(i-3).eq.ntyp1)))) cycle + if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle + if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1)) & cycle do j=1,nterm_sccor(isccori,isccori1) v1ij=v1sccor(j,intertyp,isccori,isccori1) @@ -5975,9 +6062,9 @@ c write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi c &gloc_sc(intertyp,i-3,icg) 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, - & (v1sccor(j,intertyp,itori,itori1),j=1,6) - & ,(v2sccor(j,intertyp,itori,itori1),j=1,6) + & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,isccori,isccori1, + & (v1sccor(j,intertyp,isccori,isccori1),j=1,6) + & ,(v2sccor(j,intertyp,isccori,isccori1),j=1,6) gsccor_loc(i-3)=gsccor_loc(i-3)+gloci enddo !intertyp enddo @@ -8187,7 +8274,7 @@ c---------------------------------------------------------------------------- include 'COMMON.GEO' logical swap double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2), - & auxvec1(2),auxvec2(1),auxmat1(2,2) + & auxvec1(2),auxvec2(2),auxmat1(2,2) logical lprn common /kutas/ lprn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC