X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=4814c1de0c0eecae4f9e1da00c6c76cf5d558af3;hb=b913c8425e6d7f8ed27a8d1eda35daae04f7ae70;hp=030de74baa63e988b875036a7d3af0aab6e0c3b7;hpb=f53d2f56ffdf788a1429ecd19bfb5f3de055204d;p=unres.git diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F index 030de74..4814c1d 100644 --- a/source/cluster/wham/src-M/energy_p_new.F +++ b/source/cluster/wham/src-M/energy_p_new.F @@ -124,7 +124,6 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t #endif energia(0)=etot energia(1)=evdw -c call enerprint(energia(0),frac) #ifdef SCP14 energia(2)=evdw2-evdw2_14 energia(17)=evdw2_14 @@ -229,7 +228,8 @@ 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) +c & +wsccor*fact(1)*gsccor_loc(i) +c ROZNICA Z WHAMem enddo endif return @@ -360,11 +360,19 @@ C integer icant external icant cd print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon +c ROZNICA DODANE Z WHAM +c do i=1,210 +c do j=1,2 +c eneps_temp(j,i)=0.0d0 +c enddo +c enddo +cROZNICA + evdw=0.0D0 evdw_t=0.0d0 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) @@ -379,7 +387,7 @@ 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)) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -393,6 +401,11 @@ c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj e2=fac*bb(itypi,itypj) evdwij=e1+e2 ij=icant(itypi,itypj) +c ROZNICA z WHAM +c eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij) +c eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij +c + 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)/)') @@ -529,7 +542,7 @@ c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw_t=0.0d0 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) @@ -540,7 +553,7 @@ C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@ -633,7 +646,7 @@ c endif ind=0 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) @@ -649,7 +662,7 @@ C do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) @@ -762,7 +775,7 @@ c if (icall.gt.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) @@ -778,7 +791,7 @@ C do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) @@ -901,7 +914,7 @@ c if (icall.gt.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) - if (itypi.eq.21) cycle + if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) @@ -917,7 +930,7 @@ C do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) - if (itypj.eq.21) cycle + if (itypj.eq.ntyp1) cycle dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) r0ij=r0(itypi,itypj) @@ -1806,7 +1819,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) @@ -1820,7 +1833,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) @@ -2556,7 +2569,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 @@ -2757,7 +2770,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)) @@ -2770,7 +2783,7 @@ c & " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i)) do j=iscpstart(i,iint),iscpend(i,iint) itypj=iabs(itype(j)) - if (itypj.eq.21) cycle + 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 @@ -3029,7 +3042,7 @@ c estr=0.0d0 estr1=0.0d0 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) @@ -3053,7 +3066,7 @@ c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included c do i=nnt,nct iti=iabs(itype(i)) - if (iti.ne.10 .and. iti.ne.21) then + if (iti.ne.10 .and. iti.ne.ntyp1) then nbi=nbondterm(iti) if (nbi.eq.1) then diff=vbld(i+nres)-vbldsc0(1,iti) @@ -3122,14 +3135,14 @@ C & delthe0,sig0inv,sigtc,sigsqtc,delthec,it double precision y(2),z(2) delta=0.02d0*pi - time11=dexp(-2*time) - time12=1.0d0 +c time11=dexp(-2*time) +c 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 - 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) @@ -3145,11 +3158,11 @@ C Zero the energy function and its derivative at 0 or pi. ichir21=isign(1,itype(i)) ichir22=isign(1,itype(i)) endif - 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) - icrc=0 - call proc_proc(phii,icrc) +c icrc=0 +c call proc_proc(phii,icrc) if (icrc.eq.1) phii=150.0 #else phii=phi(i) @@ -3160,11 +3173,11 @@ 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 - call proc_proc(phii1,icrc) +c icrc=0 +c call proc_proc(phii1,icrc) if (icrc.eq.1) phii1=150.0 phii1=pinorm(phii1) z(1)=cos(phii1) @@ -3232,7 +3245,7 @@ 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 +c 1215 continue enddo C Ufff.... We've done all this!!! return @@ -3375,38 +3388,43 @@ 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 +c if (itype(i-1).eq.ntyp1) cycle + if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + &(itype(i).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) CC Ta zmina jest niewlasciwa - 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 .and. itype(i-2).ne.21) 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 +c ityp1=nthetyp+1 do k=1,nsingle + ityp1=ithetyp((itype(i-2))) cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.21) 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 @@ -3414,14 +3432,15 @@ CC Ta zmina jest niewlasciwa #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 +c ityp3=nthetyp+1 + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -3430,7 +3449,7 @@ CC Ta zmina jest niewlasciwa 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) @@ -3452,11 +3471,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 @@ -3475,24 +3495,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) @@ -3500,28 +3520,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) @@ -3536,7 +3557,8 @@ c call flush(iout) etheta=etheta+ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 - gloc(nphi+i-2,icg)=wang*dethetai +c gloc(nphi+i-2,icg)=wang*dethetai + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai enddo return end @@ -3567,7 +3589,7 @@ 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(iabs(it)) c print *,'i=',i,' it=',it,' nlobit=',nlobit @@ -3860,7 +3882,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))) @@ -3887,7 +3909,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) @@ -3909,7 +3931,7 @@ c do j = 1,3 xx = xx + x_prime(j)*dc_norm(j,i+nres) yy = yy + y_prime(j)*dc_norm(j,i+nres) - zz = zz + dsign(1.0,itype(i))*z_prime(j)*dc_norm(j,i+nres) + zz = zz + z_prime(j)*dc_norm(j,i+nres) enddo xxtab(i)=xx @@ -3927,7 +3949,8 @@ 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)) +c 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 @@ -4098,8 +4121,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)) @@ -4251,8 +4276,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) @@ -4336,8 +4361,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 if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215 if (iabs(itype(i)).eq.20) then iblock=2 @@ -4435,8 +4460,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)) @@ -4517,26 +4542,48 @@ 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 +c gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) + enddo +c write (iout,*) "EBACK_SC_COR",i,esccor,intertyp +c gloc_sc(intertyp,i-3)=gloc_sc(intertyp,i-3)+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) + & (v1sccor(j,1,itori,itori1),j=1,6), + & (v2sccor(j,1,itori,itori1),j=1,6) gsccor_loc(i-3)=gloci + enddo !intertyp enddo return end @@ -5205,7 +5252,11 @@ C--------------------------------------------------------------------------- & auxmat(2,2) iti1 = itortyp(itype(i+1)) if (j.lt.nres-1) then - itj1 = itortyp(itype(j+1)) + if (itype(j).le.ntyp) then + itj1 = itortyp(itype(j+1)) + else + itj1=ntortyp+1 + endif else itj1=ntortyp+1 endif @@ -5293,14 +5344,16 @@ cd if (i.ne.2 .or. j.ne.4 .or. k.ne.3 .or. l.ne.5) return enddo if (l.eq.j+1) then C parallel orientation of the two CA-CA-CA frames. - if (i.gt.1) then +c if (i.gt.1) then + if (i.gt.1 .and. itype(i).le.ntyp) 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 +c if (l.lt.nres-1) then + if (l.lt.nres-1 .and. itype(l+1).le.ntyp) then itl1=itortyp(itype(l+1)) else itl1=ntortyp+1 @@ -5446,7 +5499,8 @@ C Calculate the Cartesian derivatives of the vectors. C End vectors else C Antiparallel orientation of the two CA-CA-CA frames. - if (i.gt.1) then +c if (i.gt.1) then + if (i.gt.1 .and. itype(i).le.ntyp) then iti=itortyp(itype(i)) else iti=ntortyp+1 @@ -5454,7 +5508,8 @@ C Antiparallel orientation of the two CA-CA-CA frames. itk1=itortyp(itype(k+1)) itl=itortyp(itype(l)) itj=itortyp(itype(j)) - if (j.lt.nres-1) then +c if (j.lt.nres-1) then + if (j.lt.nres-1 .and. itype(j+1).le.ntyp) then itj1=itortyp(itype(j+1)) else itj1=ntortyp+1 @@ -6407,7 +6462,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 @@ -6612,14 +6667,16 @@ 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 +c if (j.lt.nres-1) then + if (j.lt.nres-1 .and. itype(j+1).le.ntyp) 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 +c if (l.lt.nres-1) then + if (l.lt.nres-1 .and. itype(l+1).le.ntyp) then itl1=itortyp(itype(l+1)) else itl1=ntortyp+1 @@ -6732,13 +6789,15 @@ 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 +c if (j.lt.nres-1) then + if (j.lt.nres-1 .and. itype(j+1).le.ntyp) then itj1=itortyp(itype(j+1)) else itj1=ntortyp+1 endif itk=itortyp(itype(k)) - if (k.lt.nres-1) then +c if (k.lt.nres-1) then + if (k.lt.nres-1 .and. itype(k+1).le.ntyp) then itk1=itortyp(itype(k+1)) else itk1=ntortyp+1