X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=26eccc3bd38f7fc3bd8bff4846f2d699c60057c5;hb=12f81ba32aeccd1c2983dae882fa969432a5c33f;hp=734acc1551f5a31a691636251dc6f31791aeff23;hpb=b4662dbad52b91578a5cda22124037728093c6ed;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 734acc1..26eccc3 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -434,7 +434,7 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef MPI include 'mpif.h' double precision gradbufc(3,maxres),gradbufx(3,maxres), - & glocbuf(4*maxres),gradbufc_sum(3,maxres) + & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres) #endif include 'COMMON.SETUP' include 'COMMON.IOUNITS' @@ -447,6 +447,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.CONTROL' include 'COMMON.TIME1' include 'COMMON.MAXGRAD' + include 'COMMON.SCCOR' #ifdef TIMING time01=MPI_Wtime() #endif @@ -689,7 +690,6 @@ c enddo & +wturn3*gel_loc_turn3(i) & +wturn6*gel_loc_turn6(i) & +wel_loc*gel_loc_loc(i) - & +wsccor*gsccor_loc(i) enddo #ifdef DEBUG write (iout,*) "gloc after adding corr" @@ -708,6 +708,21 @@ 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 + do j=1,1 + write (iout,*) i,j,gloc_sc(j,i,icg) + enddo + enddo +#endif +#undef DEBUG + do i=1,nres + do j=1,3 + gloc_scbuf(j,i)=gloc_sc(j,i,icg) + enddo + enddo time00=MPI_Wtime() call MPI_Barrier(FG_COMM,IERR) time_barrier_g=time_barrier_g+MPI_Wtime()-time00 @@ -719,6 +734,19 @@ c enddo call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 + 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 + do j=1,1 + write (iout,*) i,j,gloc_sc(j,i,icg) + enddo + enddo +#endif +#undef DEBUG #ifdef DEBUG write (iout,*) "gloc after reduce" do i=1,4*nres @@ -1026,7 +1054,7 @@ c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=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) @@ -1041,7 +1069,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 @@ -1179,7 +1207,7 @@ c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=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) @@ -1190,7 +1218,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 @@ -1272,7 +1300,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) @@ -1289,7 +1317,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 c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) @@ -1392,7 +1420,7 @@ c if (icall.eq.0) lprn=.false. 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) @@ -1411,7 +1439,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 c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, @@ -1537,7 +1565,7 @@ c if (icall.eq.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) @@ -1554,7 +1582,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 c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) @@ -1785,7 +1813,7 @@ cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct evdw=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) @@ -1798,7 +1826,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 @@ -1866,7 +1894,7 @@ cd write(iout,*) 'In EELEC_soft_sphere' eello_turn4=0.0d0 ind=0 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 dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -1876,7 +1904,7 @@ cd write(iout,*) 'In EELEC_soft_sphere' 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 ind=ind+1 iteli=itel(i) itelj=itel(j) @@ -2755,8 +2783,8 @@ C C Loop over i,i+2 and i,i+3 pairs of the peptide groups C do i=iturn3_start,iturn3_end - if (itype(i).eq.21 .or. itype(i+1).eq.21 - & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21) cycle + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 + & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2772,9 +2800,9 @@ C num_cont_hb(i)=num_conti enddo do i=iturn4_start,iturn4_end - if (itype(i).eq.21 .or. itype(i+1).eq.21 - & .or. itype(i+3).eq.21 - & .or. itype(i+4).eq.21) cycle + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 + & .or. itype(i+3).eq.ntyp1 + & .or. itype(i+4).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2786,7 +2814,7 @@ C zmedi=c(3,i)+0.5d0*dzi num_conti=num_cont_hb(i) call eelecij(i,i+3,ees,evdw1,eel_loc) - if (wturn4.gt.0.0d0 .and. itype(i+2).ne.21) + if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & call eturn4(i,eello_turn4) num_cont_hb(i)=num_conti enddo ! i @@ -2794,7 +2822,7 @@ 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 - if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2808,7 +2836,7 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) c write (iout,*) i,j,itype(i),itype(j) - if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle + if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j num_cont_hb(i)=num_conti @@ -3802,7 +3830,7 @@ C cd print '(a)','Enter ESCP' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e 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) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) @@ -3811,7 +3839,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) - if (itype(j).eq.21) cycle + if (itype(j).eq.ntyp1) cycle itypj=iabs(itype(j)) C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi @@ -3898,7 +3926,7 @@ C cd print '(a)','Enter ESCP' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e 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) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) @@ -3908,7 +3936,7 @@ cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e 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 @@ -4180,7 +4208,7 @@ c estr=0.0d0 estr1=0.0d0 do i=ibondp_start,ibondp_end - 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) @@ -4205,7 +4233,7 @@ 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)) - 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) @@ -4278,7 +4306,7 @@ c time12=1.0d0 etheta=0.0D0 c write (*,'(a,i2)') 'EBEND ICG=',icg 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) @@ -4295,7 +4323,7 @@ C Zero the energy function and its derivative at 0 or pi. 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) if (phii.ne.phii) phii=150.0 @@ -4308,7 +4336,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) if (phii1.ne.phii1) phii1=150.0 @@ -4516,24 +4544,27 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 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(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-2).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))) +C propagation of chirality for glycine type do k=1,nsingle cosph1(k)=dcos(k*phii) sinph1(k)=dsin(k*phii) @@ -4546,7 +4577,7 @@ C 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 @@ -4554,7 +4585,7 @@ 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) @@ -4567,7 +4598,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) @@ -4580,20 +4611,21 @@ C sinph1ph2(k,l)=scl-csl enddo enddo - if (lprn) then +c if (lprn) then write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2, & " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1 - write (iout,*) "coskt and sinkt" - do k=1,nntheterm - write (iout,*) k,coskt(k),sinkt(k) - enddo - endif +c write (iout,*) "coskt and sinkt" +c do k=1,nntheterm +c write (iout,*) k,coskt(k),sinkt(k) +c enddo +c 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 @@ -4612,24 +4644,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) @@ -4637,28 +4669,29 @@ 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) @@ -4667,9 +4700,12 @@ C enddo enddo 10 continue - if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') +c lprn1=.true. +c if (lprn1) + write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') & i,theta(i)*rad2deg,phii*rad2deg, & phii1*rad2deg,ethetai +c lprn1=.false. 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 @@ -4704,7 +4740,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 @@ -5003,7 +5039,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))) @@ -5052,7 +5088,8 @@ 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 + z_prime(j)*dc_norm(j,i+nres) + zz = zz + dsign(1.0,dfloat(itype(i))) + & *z_prime(j)*dc_norm(j,i+nres) enddo xxtab(i)=xx @@ -5062,7 +5099,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 @@ -5070,7 +5107,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,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 @@ -5431,8 +5468,8 @@ c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end etors_ii=0.0D0 - 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) @@ -5528,8 +5565,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 etors_ii=0.0D0 if (iabs(itype(i)).eq.20) then iblock=2 @@ -5628,9 +5665,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.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 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -5706,29 +5744,53 @@ 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=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=iabs(itype(i-2)) - itori1=iabs(itype(i-1)) + isccori=isccortyp(itype(i-2)) + isccori1=isccortyp(itype(i-1)) +c write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1) 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) + 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 + 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) + & 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 + return end c----------------------------------------------------------------------------