X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fcompare.F90;h=1dfe01e013fa88a24d0644f0e6325a1ff240445c;hb=4db3a5f0b00eb4e0dc343054d4560d5e76f548f5;hp=b65e57cb7a73c5ab657bfdb6c02de28a79f4d6db;hpb=df0d15d1ab81e01e177d3d39354e72364b294e1c;p=unres4.git diff --git a/source/unres/compare.F90 b/source/unres/compare.F90 index b65e57c..1dfe01e 100644 --- a/source/unres/compare.F90 +++ b/source/unres/compare.F90 @@ -37,18 +37,29 @@ ! include 'COMMON.NAMES' real(kind=8) :: facont=1.569D0 ! facont = (2/(1-sqrt(1-1/4)))**(1/6) integer :: ncont - integer,dimension(2,12*nres) :: icont!(2,12*nres) !(2,maxcont) (maxcont=12*maxres) + integer,dimension(2,100*nres) :: icont!(2,100*nres) !(2,maxcont) (maxcont=12*maxres) logical :: lprint !el local variables real(kind=8) :: co,rcomp - integer :: kkk,i,j,i1,i2,it1,it2,iti,itj + integer :: kkk,i,j,i1,i2,it1,it2,iti,itj,inum,jnum ncont=0 kkk=3 - do i=nnt+kkk,nct - iti=iabs(itype(i)) + do i=nnt_molec(1)+kkk,nct_molec(1) + iti=iabs(itype(i,molnum(i))) + if (molnum(i).lt.3) then + inum=i+nres + else + inum=i + endif + do j=nnt,i-kkk - itj=iabs(itype(j)) + itj=iabs(itype(j,molnum(i))) + if (molnum(j).lt.3) then + jnum=j+nres + else + jnum=j + endif if (ipot.ne.4) then ! rcomp=sigmaii(iti,itj)+1.0D0 rcomp=facont*sigmaii(iti,itj) @@ -58,8 +69,9 @@ endif ! rcomp=6.5D0 ! print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j) - if (dist(nres+i,nres+j).lt.rcomp) then + if (dist(inum,jnum).lt.rcomp) then ncont=ncont+1 + if (ncont.gt.nres*100) ncont=nres*100 icont(1,ncont)=i icont(2,ncont)=j endif @@ -70,10 +82,10 @@ do i=1,ncont i1=icont(1,i) i2=icont(2,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,1) + it2=itype(i2,1) write (iout,'(i3,2x,a,i4,2x,a,i4)') & - i,restyp(it1),i1,restyp(it2),i2 + i,restyp(it1,1),i1,restyp(it2,1),i2 enddo endif co = 0.0d0 @@ -90,7 +102,7 @@ ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' integer :: ncont,ncont_ref - integer,dimension(2,12*nres) :: icont,icont_ref !(2,12*nres) (2,maxcont) (maxcont=12*maxres) + integer,dimension(2,100*nres) :: icont,icont_ref !(2,100*nres) (2,maxcont) (maxcont=12*maxres) !el local variables integer :: i,j,nmatch nmatch=0 @@ -117,7 +129,7 @@ ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' integer :: ncont,ncont_ref - integer,dimension(2,12*nres) :: icont,icont_ref !(2,12*nres) (2,maxcont) (maxcont=12*maxres) + integer,dimension(2,100*nres) :: icont,icont_ref !(2,100*nres) (2,maxcont) (maxcont=12*maxres) !el local variables integer :: i,j,nmatch nmatch=0 @@ -149,19 +161,20 @@ ! include 'COMMON.FFIELD' ! include 'COMMON.NAMES' integer :: ncont - integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres) + integer,dimension(:,:),allocatable :: icont !(2,maxcont) (maxcont=12*maxres) integer :: nharp - integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3) + integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3) logical :: lprint,not_done real(kind=8) :: rcomp=6.0d0 !el local variables integer :: i,j,kkk,k,i1,i2,it1,it2,j1,ii1,jj1 -! allocate(icont(2,12*nres)) - + if (.not.allocated(icont)) then + allocate(icont(2,100*nres_molec(1)+1)) + endif ncont=0 kkk=0 ! print *,'nnt=',nnt,' nct=',nct - do i=nnt,nct-3 + do i=nnt_molec(1),nct_molec(1)-3 do k=1,3 c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1)) enddo @@ -181,10 +194,10 @@ do i=1,ncont i1=icont(1,i) i2=icont(2,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,1) + it2=itype(i2,1) write (iout,'(i3,2x,a,i4,2x,a,i4)') & - i,restyp(it1),i1,restyp(it2),i2 + i,restyp(it1,1),i1,restyp(it2,1),i2 enddo endif ! finding hairpins @@ -230,10 +243,10 @@ ii1=iharp(3,i) jj1=iharp(4,i) write (iout,*) - write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=i1,ii1) - write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=j1,jj1,-1) + write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=i1,ii1) + write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=j1,jj1,-1) ! do k=jj1,j1,-1 -! write (iout,'(a,i3,$)') restyp(itype(k)),k +! write (iout,'(a,i3,$)') restyp(itype(k,1)),k ! enddo enddo endif @@ -257,8 +270,8 @@ real(kind=8) :: ael6_i,ael3_i real(kind=8),dimension(2,2) :: app_,bpp_,rpp_ integer :: ncont - integer,dimension(2,12*nres) :: icont !(2,12*nres)(2,maxcont) (maxcont=12*maxres) - real(kind=8),dimension(12*nres) :: econt !(maxcont) + integer,dimension(:,:),allocatable :: icont !(2,100*nres)(2,maxcont) (maxcont=12*maxres) + real(kind=8),dimension(:),allocatable :: econt !(maxcont) !el local variables integer :: i,j,k,iteli,itelj,i1,i2,it1,it2,ic1,ic2 real(kind=8) :: elcutoff,elecutoff_14,rri,ees,evdw @@ -279,8 +292,13 @@ data elpp_6 /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/ data elpp_3 / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/ -!el allocate(econt(12*nres)) !(maxcont) - +!el allocate(econt(100*nres)) !(maxcont) + if (.not.allocated(icont)) then + allocate(icont(2,100*nres_molec(1)+1)) + endif + if (.not.allocated(econt)) then + allocate(econt(100*nres_molec(1)+1)) + endif elcutoff = -0.3d0 elecutoff_14 = -0.5d0 if (lprint) write (iout,'(a)') & @@ -300,8 +318,9 @@ ncont=0 ees=0.0 evdw=0.0 - do 1 i=nnt,nct-2 - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) goto 1 +! print *, "nntt,nct",nnt,nct-2 + do 1 i=nnt_molec(1),nct_molec(1)-2 + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1 xi=c(1,i) yi=c(2,i) zi=c(3,i) @@ -312,7 +331,7 @@ ymedi=yi+0.5*dyi zmedi=zi+0.5*dzi do 4 j=i+2,nct-1 - if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4 + if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) goto 4 iteli=itel(i) itelj=itel(j) if (j.eq.i+2 .and. itelj.eq.2) iteli=2 @@ -363,10 +382,10 @@ do i=1,ncont i1=icont(1,i) i2=icont(2,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,1) + it2=itype(i2,1) write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') & - i,restyp(it1),i1,restyp(it2),i2,econt(i) + i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i) enddo endif ! For given residues keep only the contacts with the greatest energy. @@ -449,10 +468,10 @@ do i=1,ncont i1=icont(1,i) i2=icont(2,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,1) + it2=itype(i2,1) write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') & - i,restyp(it1),i1,restyp(it2),i2,econt(i) + i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i) enddo endif return @@ -470,16 +489,18 @@ ! include 'COMMON.CONTROL' integer :: ncont,i,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,nhelix,& iii1,jjj1 - integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres) - integer,dimension(nres,4) :: isec !(maxres,4) + integer,dimension(:,:),allocatable :: icont !(2,maxcont) (maxcont=12*maxres) + integer,dimension(nres,0:4) :: isec !(maxres,4) integer,dimension(nres) :: nsec !(maxres) logical :: lprint,not_done !,freeres real(kind=8) :: p1,p2 !el external freeres -!el allocate(icont(2,12*nres),isec(nres,4),nsec(nres)) - - if(.not.dccart) call chainbuild +!el allocate(icont(2,100*nres),isec(nres,4),nsec(nres)) + if (.not.allocated(icont)) then + allocate(icont(2,100*nres+1)) + endif + if(.not.dccart) call chainbuild_cart if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3) !d call write_pdb(99,'sec structure',0d0) ncont=0 @@ -492,6 +513,8 @@ enddo call elecont(lprint,ncont,icont) +! print *,"after elecont" + if (nres_molec(1).eq.0) return ! finding parallel beta !d write (iout,*) '------- looking for parallel beta -----------' @@ -742,7 +765,7 @@ write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'" write(12,'(a20)') "XMacStand ribbon.mac" - + if (nres_molec(1).eq.0) return write(iout,*) 'UNRES seq:' do j=1,nbfrag write(iout,*) 'beta ',(bfrag(i,j),i=1,4) @@ -815,6 +838,7 @@ ! & obr,non_conv) ! rms=dsqrt(rms) call rmsd(rms) +! print *,"before contact" !elte(iout,*) "rms_nacc before contact" call contact(.false.,ncont,icont,co) frac=contact_fract(ncont,ncont_ref,icont,icont_ref) @@ -846,7 +870,7 @@ !el local variables real(kind=8) :: drms,rminroz,roznica - integer :: i,j,iatom,kkk,iti,k + integer :: i,j,iatom,kkk,iti,k,mnum !el allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres @@ -861,7 +885,7 @@ do kkk=1,nperm ! do i=nz_start,nz_end ! iatom=iatom+1 -! iti=itype(i) +! iti=itype(i,1) ! do k=1,3 ! ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup) ! crefcopy(k,iatom,kkk)=cref(k,i,kkk) @@ -877,14 +901,18 @@ ! else ! do kkk=1,nperm iatom=0 +! print *,nz_start,nz_end,nstart_seq-nstart_sup do i=nz_start,nz_end +! iatom=iatom+1 + iti=itype(i+nstart_seq-nstart_sup,molnum(i)) + if (iti.eq.ntyp1_molec(molnum(i))) cycle iatom=iatom+1 - iti=itype(i) + mnum=molnum(i+nstart_seq-nstart_sup) do k=1,3 ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup) crefcopy(k,iatom)=cref(k,i,kkk) enddo - if (iz_sc.eq.1.and.iti.ne.10) then + if (iz_sc.eq.1.and.iti.ne.10.and.mnum.lt.4) then iatom=iatom+1 do k=1,3 ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup) @@ -959,7 +987,7 @@ iatom=0 do i=nz_start,nz_end iatom=iatom+1 - iti=itype(i) + iti=itype(i,1) do k=1,3 ccopy(k,iatom)=c(k,i) crefcopy(k,iatom)=crefjlee(k,i) @@ -2353,9 +2381,9 @@ ! include 'COMMON.IOUNITS' ! include 'COMMON.DISTFIT' - integer :: ncont,icont(2,nres*nres/2),isec(nres,3) + integer :: ncont,isec(nres,3) logical :: lprint,not_done - real(kind=4) :: dcont(nres*nres/2),d + real(kind=4) :: d real(kind=4) :: rcomp = 7.0 real(kind=4) :: rbeta = 5.2 real(kind=4) :: ralfa = 5.2 @@ -2363,7 +2391,15 @@ real(kind=8),dimension(3) :: xpi,xpj integer :: i,k,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,iii1,jjj1,& nhelix - call chainbuild + integer, dimension(:,:),allocatable :: icont + real(kind=4),dimension(:),allocatable :: dcont + if (.not.allocated(icont)) then + allocate(icont(2,100*nres_molec(1)+1)) + endif + if (.not.allocated(dcont)) then + allocate(dcont(100*nres_molec(1)+1)) + endif + call chainbuild_cart !d call write_pdb(99,'sec structure',0d0) ncont=0 nbfrag=0 @@ -2390,6 +2426,7 @@ (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) + & (xpi(3)-xpj(3))*(xpi(3)-xpj(3)) if ( d.lt.rcomp*rcomp) then + if (ncont.gt.(100*nres_molec(1)+1)) ncont=100*nres_molec(1)+1 ncont=ncont+1 icont(1,ncont)=i icont(2,ncont)=j @@ -3935,20 +3972,20 @@ c(j,i)=cart_base(j,i+ist,ii) ! cref(j,i)=c(j,i) enddo -!d write (iout,'(a,i4,3f10.5)') restyp(itype(i)),i,(c(j,i),j=1,3) +!d write (iout,'(a,i4,3f10.5)') restyp(itype(i,1)),i,(c(j,i),j=1,3) enddo !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr, !d non_conv) !d write (iout,'(a,f10.5)') !d & 'Initial RMS deviation from reference structure:',rms - if (itype(nres).eq.ntyp1) then + if (itype(nres,1).eq.ntyp1) then do j=1,3 dcj=c(j,nres-2)-c(j,nres-3) c(j,nres)=c(j,nres-1)+dcj c(j,2*nres)=c(j,nres) enddo endif - if (itype(1).eq.ntyp1) then + if (itype(1,1).eq.ntyp1) then do j=1,3 dcj=c(j,4)-c(j,3) c(j,1)=c(j,2)-dcj @@ -4178,9 +4215,9 @@ link_end0=link_end link_end=min0(link_end,nss) do i=nnt,nct - if (itype(i).ne.10) then -!d print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1) - call gen_side(itype(i),theta(i+1),alph(i),omeg(i),fail) + if (itype(i,1).ne.10) then +!d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1) + call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail,1) endif enddo call chainbuild @@ -4191,12 +4228,12 @@ glycine=.true. do while(glycine) ind_sc=iran_num(nnt,nct) - glycine=(itype(ind_sc).eq.10) + glycine=(itype(ind_sc,1).eq.10) enddo alph0=alph(ind_sc) omeg0=omeg(ind_sc) - call gen_side(itype(ind_sc),theta(ind_sc+1),alph(ind_sc),& - omeg(ind_sc),fail) + call gen_side(itype(ind_sc,1),theta(ind_sc+1),alph(ind_sc),& + omeg(ind_sc),fail,1) call chainbuild call etotal(energia) !d write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))') @@ -4530,19 +4567,21 @@ if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3)) !(4,maxres/3) if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3) ! COMMON /WAGI/ - allocate(w(maxres22),d0(maxres22)) !(maxres22) +#ifndef FIVEDIAG + if(.not.allocated(w)) allocate(w(maxres22),d0(maxres22)) !(maxres22) ! COMMON /POCHODNE/ !el allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES) - allocate(DDD(maxres22)) !(maxres22) - allocate(H(nres,nres)) !(MAXRES,MAXRES) - allocate(XX(nres)) !(MAXRES) + if (.not.allocated(DDD)) allocate(DDD(maxres22)) !(maxres22) + if (.not.allocated(H)) allocate(H(nres,nres)) !(MAXRES,MAXRES) +#endif + if (.not.allocated(XX)) allocate(XX(nres)) !(MAXRES) ! COMMON /frozen/ - allocate(mask(nres)) !(maxres) + if (.not.allocated(mask)) allocate(mask(nres)) !(maxres) ! common.thread ! common /thread/ - allocate(iexam(2,maxthread),ipatt(2,maxthread)) !(2,maxthread) + if (.not.allocated(iexam))allocate(iexam(2,maxthread),ipatt(2,maxthread)) !(2,maxthread) ! common /thread1/ - allocate(ener0(n_ene+2,maxthread),ener(n_ene+2,maxthread)) !(n_ene+2,maxthread) + if (.not.allocated(ener0)) allocate(ener0(n_ene+2,maxthread),ener(n_ene+2,maxthread)) !(n_ene+2,maxthread) return end subroutine alloc_compare_arrays