X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=8f932cca8c6b411e2e0b20f82a69e8e2ebff5c33;hb=7675ab8ec26554f2fd3f2e7e427b177254872a45;hp=bd7cbaed113500eacd9e7d636635b3704dbb2cbf;hpb=cea715862444dcb47c99e54db632e259844221b8;p=unres.git diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index bd7cbae..8f932cc 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -228,7 +228,6 @@ C & +wturn3*fact(2)*gel_loc_turn3(i) & +wturn6*fact(5)*gel_loc_turn6(i) & +wel_loc*fact(2)*gel_loc_loc(i) - & +wsccor*fact(1)*gsccor_loc(i) enddo endif return @@ -368,9 +367,9 @@ cd print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 evdw_t=0.0d0 do i=iatsc_s,iatsc_e - itypi=itype(i) - if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi=iabs(itype(i)) + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -383,8 +382,8 @@ 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=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -541,9 +540,9 @@ c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 evdw_t=0.0d0 do i=iatsc_s,iatsc_e - itypi=itype(i) - if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi=iabs(itype(i)) + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -552,8 +551,8 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -654,9 +653,9 @@ c else c endif ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) - if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi=iabs(itype(i)) + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -670,8 +669,8 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) @@ -792,9 +791,9 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.gt.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) - if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi=iabs(itype(i)) + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -808,8 +807,8 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) @@ -942,9 +941,9 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.gt.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) - if (itypi.eq.21) cycle - itypi1=itype(i+1) + itypi=iabs(itype(i)) + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -958,8 +957,8 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) r0ij=r0(itypi,itypj) @@ -1853,7 +1852,7 @@ cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e gcorr_loc(i)=0.0d0 enddo do i=iatel_s,iatel_e - if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle if (itel(i).eq.0) goto 1215 dxi=dc(1,i) dyi=dc(2,i) @@ -1867,7 +1866,7 @@ cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) - if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle + if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle if (itel(j).eq.0) goto 1216 ind=ind+1 iteli=itel(i) @@ -1917,10 +1916,15 @@ C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) ees=ees+eesij evdw1=evdw1+evdwij -cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') -cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, -cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, -cd & xmedi,ymedi,zmedi,xj,yj,zj +c write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') +c &'evdw1',i,j,evdwij +c &,iteli,itelj,aaa,evdw1 + +c write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij +c write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') +c & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, +c & 1.0D0/dsqrt(rrmij),evdwij,eesij, +c & xmedi,ymedi,zmedi,xj,yj,zj C C Calculate contributions to the Cartesian gradient. C @@ -2266,8 +2270,10 @@ C Check the loc-el terms by numerical integration C Contribution to the local-electrostatic energy coming from the i-j pair eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) -cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij -cd write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) +c write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij +c write (iout,'(a6,2i5,0pf7.3)') +c & 'eelloc',i,j,eel_loc_ij +c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) eel_loc=eel_loc+eel_loc_ij C Partial derivatives in virtual-bond dihedral angles gamma if (calc_grad) then @@ -2603,7 +2609,7 @@ C Cartesian derivatives & +0.5d0*(pizda(1,1)+pizda(2,2)) enddo endif - else if (j.eq.i+3 .and. itype(i+2).ne.21) then + else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Fourth-order contributions @@ -2804,7 +2810,7 @@ cd print '(a)','Enter ESCP' c write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e, c & ' scal14',scal14 do i=iatscp_s,iatscp_e - if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) c write (iout,*) "i",i," iteli",iteli," nscp_gr",nscp_gr(i), c & " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i)) @@ -2816,8 +2822,8 @@ c & " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i)) do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=itype(j) - if (itypj.eq.21) cycle + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi @@ -2836,7 +2842,9 @@ C Uncomment following three lines for Ca-p interactions evdw2_14=evdw2_14+e1+e2 endif evdwij=e1+e2 -c write (iout,*) i,j,evdwij +c write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') +c & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli), +c & bad(itypj,iteli) evdw2=evdw2+evdwij if (calc_grad) then C @@ -2927,7 +2935,8 @@ C iii and jjj point to the residues for which the distance is assigned. endif 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 + if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. + & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij else @@ -2987,7 +2996,7 @@ C include 'COMMON.VAR' include 'COMMON.IOUNITS' double precision erij(3),dcosom1(3),dcosom2(3),gg(3) - itypi=itype(i) + itypi=iabs(itype(i)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -2995,7 +3004,7 @@ C dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) dsci_inv=dsc_inv(itypi) - itypj=itype(j) + itypj=iabs(itype(j)) dscj_inv=dsc_inv(itypj) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -3073,9 +3082,10 @@ c logical energy_dec /.false./ double precision u(3),ud(3) estr=0.0d0 - write (iout,*) "distchainmax",distchainmax + estr1=0.0d0 +c write (iout,*) "distchainmax",distchainmax do i=nnt+1,nct - if (itype(i-1).eq.21 .or. itype(i).eq.21) then + if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) do j=1,3 gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) @@ -3094,13 +3104,13 @@ c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff endif enddo - estr=0.5d0*AKP*estr + estr=0.5d0*AKP*estr+estr1 c c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included c do i=nnt,nct - iti=itype(i) - if (iti.ne.10 .and. iti.ne.21) then + iti=iabs(itype(i)) + if (iti.ne.10 .and. iti.ne.ntyp1) then nbi=nbondterm(iti) if (nbi.eq.1) then diff=vbld(i+nres)-vbldsc0(1,iti) @@ -3176,11 +3186,24 @@ 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 - if (itype(i-1).eq.21) cycle + if (itype(i-1).eq.ntyp1) cycle 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.3 .and. itype(i-2).ne.21) then + 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 .and. itype(i-2).ne.ntyp1) then #ifdef OSF phii=phi(i) icrc=0 @@ -3195,7 +3218,7 @@ C Zero the energy function and its derivative at 0 or pi. y(1)=0.0D0 y(2)=0.0D0 endif - if (i.lt.nres .and. itype(i).ne.21) then + if (i.lt.nres .and. itype(i).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) icrc=0 @@ -3217,8 +3240,12 @@ 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 c write (iout,*) "thet_pred_mean",thet_pred_mean @@ -3226,8 +3253,16 @@ c write (iout,*) "thet_pred_mean",thet_pred_mean 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 + 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) @@ -3398,24 +3433,26 @@ C etheta=0.0D0 c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) do i=ithet_start,ithet_end - if (itype(i-1).eq.21) cycle + if (itype(i-1).eq.ntyp1) cycle + if (iabs(itype(i+1)).eq.20) iblock=2 + if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 dephii=0.0d0 dephii1=0.0d0 theti2=0.5d0*theta(i) - ityp2=ithetyp(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 .and. itype(i-2).ne.21) then + if (i.gt.3 .and. itype(i-2).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 #else phii=phi(i) #endif - ityp1=ithetyp(itype(i-2)) + ityp1=ithetyp((itype(i-2))) do k=1,nsingle cosph1(k)=dcos(k*phii) sinph1(k)=dsin(k*phii) @@ -3428,7 +3465,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.21) then + if (i.lt.nres .and. itype(i).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -3436,7 +3473,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) #else phii1=phi(i+1) #endif - ityp3=ithetyp(itype(i)) + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=dcos(k*phii1) sinph2(k)=dsin(k*phii1) @@ -3452,7 +3489,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) c write (iout,*) "i",i," ityp1",itype(i-2),ityp1, c & " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3 c call flush(iout) - 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) @@ -3474,11 +3511,12 @@ c call flush(iout) 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 @@ -3497,24 +3535,24 @@ c call flush(iout) 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) @@ -3522,28 +3560,29 @@ c call flush(iout) 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) @@ -3589,9 +3628,9 @@ C ALPHA and OMEGA. c write (iout,'(a)') 'ESC' do i=loc_start,loc_end it=itype(i) - if (it.eq.21) cycle + if (it.eq.ntyp1) cycle if (it.eq.10) goto 1 - nlobit=nlob(it) + nlobit=nlob(iabs(it)) c print *,'i=',i,' it=',it,' nlobit=',nlobit c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad theti=theta(i+1)-pipol @@ -3746,7 +3785,7 @@ 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) + expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin) cd print *,'j=',j,' expfac=',expfac escloc_i=escloc_i+expfac do k=1,3 @@ -3827,7 +3866,7 @@ 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) + expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin) escloc_i=escloc_i+expfac do k=1,2 dersc(k)=dersc(k)+Ax(k,j)*expfac @@ -3882,7 +3921,7 @@ C delta=0.02d0*pi escloc=0.0D0 do i=loc_start,loc_end - if (itype(i).eq.21) cycle + if (itype(i).eq.ntyp1) cycle costtab(i+1) =dcos(theta(i+1)) sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1)) cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1))) @@ -3891,7 +3930,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 @@ -3909,7 +3948,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) @@ -3941,7 +3980,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 @@ -3949,7 +3988,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.0d0,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 @@ -3992,6 +4031,8 @@ c & dscp1,dscp2,sumene c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) escloc = escloc + sumene c write (2,*) "escloc",escloc +c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i), +c & zz,xx,yy if (.not. calc_grad) goto 1 #ifdef DEBUG C @@ -4120,8 +4161,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)) @@ -4273,8 +4316,8 @@ C Set lprn=.true. for debugging c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end - if (itype(i-2).eq.21 .or. itype(i-1).eq.21 - & .or. itype(i).eq.21) cycle + if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 + & .or. itype(i).eq.ntyp1) cycle itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) phii=phi(i) @@ -4358,17 +4401,22 @@ C Set lprn=.true. for debugging c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end - if (itype(i-2).eq.21 .or. itype(i-1).eq.21 - & .or. itype(i).eq.21) cycle + if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 + & .or. itype(i).eq.ntyp1) cycle if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215 + 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 @@ -4381,22 +4429,24 @@ 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) pom=vl2ij*cosphi+vl3ij*sinphi pom1=1.0d0/(pom*pom+1.0d0) etors=etors+vl1ij*pom1 +c if (energy_dec) etors_ii=etors_ii+ +c & 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) + etors=etors-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,1),j=1,6),(v2(j,itori,itori1,1),j=1,6) gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) 1215 continue @@ -4452,8 +4502,8 @@ C Set lprn=.true. for debugging c lprn=.true. etors_d=0.0D0 do i=iphi_start,iphi_end-1 - if (itype(i-2).eq.21 .or. itype(i-1).eq.21 - & .or. itype(i).eq.21 .or. itype(i+1).eq.21) cycle + if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 + & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle if (itel(i-2).eq.0 .or. itel(i-1).eq.0 .or. itel(i).eq.0) & goto 1215 itori=itortyp(itype(i-2)) @@ -4463,12 +4513,14 @@ c lprn=.true. phii1=phi(i+1) gloci1=0.0D0 gloci2=0.0D0 + 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) - 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) + 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) @@ -4478,12 +4530,12 @@ C Regular cosine and sine terms 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) @@ -4493,7 +4545,7 @@ C Regular cosine and sine terms gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2 & -v1cdij*sinphi1p2-v2cdij*sinphi1m2) gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2 - & -v1cdij*sinphi1p2+v2cdij*sinphi1m2) + & -v1cdij*sinphi1p2+v2cdij*sinphi1m2) enddo enddo gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*fact2*gloci1 @@ -4532,26 +4584,49 @@ C Set lprn=.true. for debugging c lprn=.true. c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor esccor=0.0D0 - do i=iphi_start,iphi_end - if (itype(i-2).eq.21 .or. itype(i-1).eq.21) cycle + do i=itau_start,itau_end + if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle esccor_ii=0.0D0 - itori=itype(i-2) - itori1=itype(i-1) + isccori=isccortyp(itype(i-2)) + isccori1=isccortyp(itype(i-1)) phii=phi(i) + do intertyp=1,3 !intertyp +cc Added 09 May 2012 (Adasko) +cc Intertyp means interaction type of backbone mainchain correlation: +c 1 = SC...Ca...Ca...Ca +c 2 = Ca...Ca...Ca...SC +c 3 = SC...Ca...Ca...SCi gloci=0.0D0 - do j=1,nterm_sccor - v1ij=v1sccor(j,itori,itori1) - v2ij=v2sccor(j,itori,itori1) - cosphi=dcos(j*phii) - sinphi=dsin(j*phii) - esccor=esccor+v1ij*cosphi+v2ij*sinphi - gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) - enddo + if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. + & (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.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.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) + v2ij=v2sccor(j,intertyp,isccori,isccori1) + cosphi=dcos(j*tauangle(intertyp,i)) + sinphi=dsin(j*tauangle(intertyp,i)) + esccor=esccor+v1ij*cosphi+v2ij*sinphi + gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) + enddo +c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp, +c & nterm_sccor(isccori,isccori1),isccori,isccori1 +c gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci 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,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6) - gsccor_loc(i-3)=gloci + & (v1sccor(j,1,itori,itori1),j=1,6) + & ,(v2sccor(j,1,itori,itori1),j=1,6) +c gsccor_loc(i-3)=gloci + enddo !intertyp enddo return end @@ -4660,9 +4735,9 @@ c------------------------------------------------------------------------------ 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) + common /contacts_hb/ zapas(3,ntyp,maxres,7), + & facont_hb(ntyp,maxres),ees0p(ntyp,maxres),ees0m(ntyp,maxres), + & num_cont_hb(maxres),jcont_hb(ntyp,maxres) num_kont=num_cont_hb(atom) do i=1,num_kont do k=1,7 @@ -4685,9 +4760,10 @@ c------------------------------------------------------------------------------ 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) + common /contacts_hb/ zapas(3,ntyp,maxres,7), + & facont_hb(ntyp,maxres),ees0p(ntyp,maxres), + & ees0m(ntyp,maxres), + & num_cont_hb(maxres),jcont_hb(ntyp,maxres) num_kont=buffer(1,indx+26) num_kont_old=num_cont_hb(atom) num_cont_hb(atom)=num_kont+num_kont_old