X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=47583a4e0791d0fcb764701c868fb9a2b3893257;hb=050ed1574772b7d2bdee4225feda68c312beca75;hp=e943ce5c3ebaee31375293623675d31ca3ba4939;hpb=2613142be7b5226fb002eaceb389e3f26766aa65;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 e943ce5..47583a4 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 @@ -774,7 +780,6 @@ c enddo do i=1,4*nres glocbuf(i)=gloc(i,icg) enddo -#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc before reduce" do i=1,nres @@ -783,7 +788,6 @@ c enddo enddo enddo #endif -#undef DEBUG do i=1,nres do j=1,3 gloc_scbuf(j,i)=gloc_sc(j,i,icg) @@ -802,7 +806,6 @@ c enddo call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 -#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc after reduce" do i=1,nres @@ -811,7 +814,6 @@ c enddo enddo enddo #endif -#undef DEBUG #ifdef DEBUG write (iout,*) "gloc after reduce" do i=1,4*nres @@ -1056,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)'/ @@ -1137,8 +1139,8 @@ C c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1151,7 +1153,7 @@ C 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=iabs(itype(j)) + itypj=itype(j) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -1314,8 +1316,8 @@ C c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1324,7 +1326,7 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=iabs(itype(j)) + itypj=itype(j) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -1431,8 +1433,8 @@ c else c endif ind=0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1556,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. @@ -1567,8 +1570,8 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.false. ind=0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1584,8 +1587,14 @@ 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' + ELSE ind=ind+1 - itypj=iabs(itype(j)) + itypj=itype(j) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, @@ -1693,6 +1702,7 @@ C Calculate angular part of the gradient. #else call sc_grad #endif + ENDIF ! dyn_ss enddo ! j enddo ! iint enddo ! i @@ -1726,8 +1736,8 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1742,7 +1752,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=iabs(itype(j)) + itypj=itype(j) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) @@ -2046,8 +2056,8 @@ C cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -2058,7 +2068,7 @@ C 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=iabs(itype(j)) + itypj=itype(j) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -3013,6 +3023,9 @@ C C Loop over i,i+2 and i,i+3 pairs of the peptide groups C do i=iturn3_start,iturn3_end +C if (itype(i).eq.21 .or. itype(i+1).eq.21 +C & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21.or.itype(i+4).eq.21) +C & cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3028,6 +3041,10 @@ C num_cont_hb(i)=num_conti enddo do i=iturn4_start,iturn4_end +C if (itype(i).eq.21 .or. itype(i+1).eq.21 +C & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21.or.itype(i+4).eq.21 +C & .or. itype(i+5).eq.21) +C & cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3046,6 +3063,8 @@ c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c do i=iatel_s,iatel_e +C if (itype(i).eq.21 .or. itype(i+1).eq.21 +C &.or.itype(i+2)) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3058,6 +3077,8 @@ c c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) +C if (itype(j).eq.21 .or. itype(j+1).eq.21 +C &.or.itype(j+2)) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j num_cont_hb(i)=num_conti @@ -4059,7 +4080,7 @@ cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=iabs(itype(j)) + itypj=itype(j) C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi @@ -4153,7 +4174,7 @@ cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=iabs(itype(j)) + itypj=itype(j) C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi @@ -4270,10 +4291,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. iabs(itype(iii)).eq.1 .and. iabs(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 @@ -4375,7 +4401,7 @@ C include 'COMMON.VAR' include 'COMMON.IOUNITS' double precision erij(3),dcosom1(3),dcosom2(3),gg(3) - itypi=iabs(itype(i)) + itypi=itype(i) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -4384,7 +4410,7 @@ C dzi=dc_norm(3,nres+i) c dsci_inv=dsc_inv(itypi) dsci_inv=vbld_inv(nres+i) - itypj=iabs(itype(j)) + itypj=itype(j) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(nres+j) xj=c(1,nres+j)-xi @@ -4412,7 +4438,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, @@ -4478,7 +4504,7 @@ c c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included c do i=ibond_start,ibond_end - iti=iabs(itype(i)) + iti=itype(i) if (iti.ne.10) then nbi=nbondterm(iti) if (nbi.eq.1) then @@ -4553,19 +4579,7 @@ c write (*,'(a,i2)') 'EBEND ICG=',icg 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)) - 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 + it=itype(i-1) if (i.gt.3) then #ifdef OSF phii=phi(i) @@ -4599,27 +4613,15 @@ 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,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 + athetk=athet(k,it) + bthetk=bthet(k,it) 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,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 + 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) @@ -4787,36 +4789,39 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end + if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + &(itype(i).eq.ntyp1)) cycle dethetai=0.0d0 dephii=0.0d0 dephii1=0.0d0 theti2=0.5d0*theta(i) - ityp2=ithetyp(iabs(itype(i-1))) + ityp2=ithetyp(itype(i-1)) do k=1,nntheterm coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3) then +C if (i.gt.3) then + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 #else phii=phi(i) #endif - ityp1=ithetyp(iabs(itype(i-2))) + ityp1=ithetyp(itype(i-2)) do k=1,nsingle cosph1(k)=dcos(k*phii) sinph1(k)=dsin(k*phii) enddo else phii=0.0d0 - ityp1=nthetyp+1 + ityp1=ithetyp(itype(i-2)) do k=1,nsingle cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif - if (i.lt.nres) then + if ((i.lt.nres).and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4824,14 +4829,14 @@ C #else phii1=phi(i+1) #endif - ityp3=ithetyp(iabs(itype(i))) + ityp3=ithetyp(itype(i)) do k=1,nsingle cosph2(k)=dcos(k*phii1) sinph2(k)=dsin(k*phii1) enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp(itype(i)) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -4937,8 +4942,8 @@ C enddo enddo 10 continue - if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') - & i,theta(i)*rad2deg,phii*rad2deg, + if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)') + & 'ebe', i,theta(i)*rad2deg,phii*rad2deg, & phii1*rad2deg,ethetai etheta=etheta+ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii @@ -4975,7 +4980,7 @@ c write (iout,'(a)') 'ESC' do i=loc_start,loc_end it=itype(i) if (it.eq.10) goto 1 - nlobit=nlob(iabs(it)) + 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 @@ -5132,11 +5137,11 @@ C Compute the contribution to SC energy and derivatives do j=1,nlobit #ifdef OSF - adexp=bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin + adexp=bsc(j,it)-0.5D0*contr(j,iii)+emin if(adexp.ne.adexp) adexp=1.0 expfac=dexp(adexp) #else - expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin) + expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin) #endif cd print *,'j=',j,' expfac=',expfac escloc_i=escloc_i+expfac @@ -5218,7 +5223,7 @@ C Compute the contribution to SC energy and derivatives dersc12=0.0d0 do j=1,nlobit - expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin) + 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 @@ -5797,17 +5802,12 @@ c lprn=.true. etors_ii=0.0D0 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) - if (iabs(itype(i)).eq.20) then - iblock=2 - else - iblock=1 - endif phii=phi(i) gloci=0.0D0 C Regular cosine and sine terms - do j=1,nterm(itori,itori1,iblock) - v1ij=v1(j,itori,itori1,iblock) - v2ij=v2(j,itori,itori1,iblock) + 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 @@ -5822,7 +5822,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,iblock) + do j=1,nlor(itori,itori1) vl1ij=vlor1(j,itori,itori1) vl2ij=vlor2(j,itori,itori1) vl3ij=vlor3(j,itori,itori1) @@ -5835,14 +5835,13 @@ C gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom enddo C Subtract the constant term - etors=etors-v0(itori,itori1,iblock) + etors=etors-v0(itori,itori1) if (energy_dec) write (iout,'(a6,i5,0pf7.3)') - & 'etor',i,etors_ii-v0(itori,itori1,iblock) + & 'etor',i,etors_ii-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,iblock),j=1,6), - & (v2(j,itori,itori1,iblock),j=1,6) + & (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 @@ -5897,17 +5896,15 @@ c lprn=.true. itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) - iblock=1 - if (iabs(itype(i+1)).eq.20) iblock=2 phii=phi(i) phii1=phi(i+1) gloci1=0.0D0 gloci2=0.0D0 - 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) + 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) @@ -5917,12 +5914,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,iblock) + do k=2,ntermd_2(itori,itori1,itori2) do l=1,k-1 - 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) + 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) @@ -5971,18 +5968,22 @@ c lprn=.true. c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor esccor=0.0D0 do i=itau_start,itau_end +C do i=42,42 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) + cccc Added 9 May 2012 cc Tauangle is torsional engle depending on the value of first digit c(see comment below) cc Omicron is flat angle depending on the value of first digit c(see comment below) - +C print *,i,tauangle(1,i) - do intertyp=1,3 !intertyp +c do intertyp=1,3 !intertyp + do intertyp=2,2 !intertyp cc Added 09 May 2012 (Adasko) cc Intertyp means interaction type of backbone mainchain correlation: c 1 = SC...Ca...Ca...Ca @@ -6007,6 +6008,7 @@ c 3 = SC...Ca...Ca...SCi esccor=esccor+v1ij*cosphi+v2ij*sinphi gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo +C print *,i,tauangle(1,i),gloci gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci c write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi c &gloc_sc(intertyp,i-3,icg) @@ -6019,7 +6021,8 @@ c &gloc_sc(intertyp,i-3,icg) enddo !intertyp enddo c do i=1,nres -c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc(i,icg) +c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc_sc(2,i,icg), +c & gloc_sc(3,i,icg) c enddo return end @@ -8224,7 +8227,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