X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_base.f90;h=0e1a986274d921551755bd4d63a8d29c37bcd804;hb=8d2ab9ba185dbbc31bfe3d1c66d7e1c9d632463b;hp=f5c7bbf741cef65daf16854e68267b95645532f0;hpb=35f220f409bd5d21be33a402d79da2c23d3e0c3a;p=unres4.git diff --git a/source/unres/io_base.f90 b/source/unres/io_base.f90 index f5c7bbf..0e1a986 100644 --- a/source/unres/io_base.f90 +++ b/source/unres/io_base.f90 @@ -5,10 +5,18 @@ implicit none !----------------------------------------------------------------------------- ! Max. number of AA residues - integer,parameter :: maxres=4000!1200 + integer,parameter :: maxres=6000!1200 ! Appr. max. number of interaction sites integer,parameter :: maxres2=2*maxres +! parameter (maxres6=6*maxres) +! parameter (mmaxres2=(maxres2*(maxres2+1)/2)) !----------------------------------------------------------------------------- +! Max. number of S-S bridges +! integer,parameter :: maxss=20 +!----------------------------------------------------------------------------- +! Max. number of derivatives of virtual-bond and side-chain vectors in theta +! or phi. +! integer,parameter :: maxdim=(maxres-1)*(maxres-2)/2 !----------------------------------------------------------------------------- ! ! @@ -60,6 +68,7 @@ write(iout,*) ' iss:',(iss(i),i=1,ns) ! Check whether the specified bridging residues are cystines. do i=1,ns + write(iout,*) i,iss(i) if (itype(iss(i)).ne.1) then if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') & 'Do you REALLY think that the residue ',& @@ -75,9 +84,13 @@ #endif endif enddo + if (dyn_ss) then + if(.not.allocated(ihpb)) allocate(ihpb(ns)) + if(.not.allocated(jhpb)) allocate(jhpb(ns)) + endif ! Read preformed bridges. if (ns.gt.0) then - read (inp,*) nss + read (inp,*) nss if (nss.gt.0) then if(.not.allocated(ihpb)) allocate(ihpb(nss)) if(.not.allocated(jhpb)) allocate(jhpb(nss)) @@ -133,7 +146,7 @@ enddo endif endif -!write(iout,*) "end read_bridge" +! write(iout,*) "end read_bridge" return end subroutine read_bridge !----------------------------------------------------------------------------- @@ -151,12 +164,11 @@ ! include 'COMMON.CONTROL' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' +! Read coordinates from input ! !el local variables integer :: l,k,j,i,kanal -! Read coordinates from input -! read(kanal,'(8f10.5)',end=10,err=10) & ((c(l,k),l=1,3),k=1,nres),& ((c(l,k+nres),l=1,3),k=nnt,nct) @@ -241,7 +253,8 @@ end subroutine read_threadbase !----------------------------------------------------------------------------- #ifdef WHAM_RUN - subroutine read_angles(kanal,iscor,energ,iprot,*) +!el subroutine read_angles(kanal,iscor,energ,iprot,*) + subroutine read_angles(kanal,*) use geometry_data use energy_data @@ -261,6 +274,8 @@ subroutine read_angles(kanal,*) use geometry_data + ! use energy + ! use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.GEO' @@ -409,9 +424,34 @@ card=card(:ilen(card)+1)//karta return end subroutine card_concat +!---------------------------------------------------------------------------- + subroutine read_afminp + use geometry_data + use energy_data + use control_data, only:out1file + use MPI_data + character*320 afmcard + integer i + print *, "wchodze" + call card_concat(afmcard,.true.) + call readi(afmcard,"BEG",afmbeg,0) + call readi(afmcard,"END",afmend,0) + call reada(afmcard,"FORCE",forceAFMconst,0.0d0) + call reada(afmcard,"VEL",velAFMconst,0.0d0) + print *,'FORCE=' ,forceAFMconst +!------ NOW PROPERTIES FOR AFM + distafminit=0.0d0 + do i=1,3 + distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit + enddo + distafminit=dsqrt(distafminit) + print *,'initdist',distafminit + return + end subroutine read_afminp !----------------------------------------------------------------------------- subroutine read_dist_constr use MPI_data + ! use control use geometry, only: dist use geometry_data use control_data @@ -438,6 +478,16 @@ ! write (iout,*) "Calling read_dist_constr" ! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup ! call flush(iout) + if(.not.allocated(ihpb)) allocate(ihpb(maxdim)) + if(.not.allocated(jhpb)) allocate(jhpb(maxdim)) + if(.not.allocated(dhpb)) allocate(dhpb(maxdim)) + if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim)) + if(.not.allocated(forcon)) allocate(forcon(maxdim)) + if(.not.allocated(fordepth)) allocate(fordepth(maxdim)) + if ((genconstr.gt.0).and.(constr_dist.eq.11)) then + call gen_dist_constr2 + go to 1712 + endif call card_concat(controlcard,.true.) call readi(controlcard,"NFRAG",nfrag_,0) call readi(controlcard,"NPAIR",npair_,0) @@ -456,11 +506,11 @@ ! do i=1,npair_ ! write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) ! enddo - if(.not.allocated(ihpb)) allocate(ihpb(maxdim)) - if(.not.allocated(jhpb)) allocate(jhpb(maxdim)) - if(.not.allocated(dhpb)) allocate(dhpb(maxdim)) - if(.not.allocated(forcon)) allocate(forcon(maxdim)) - +! if(.not.allocated(ihpb)) allocate(ihpb(maxdim)) +! if(.not.allocated(jhpb)) allocate(jhpb(maxdim)) +! if(.not.allocated(dhpb)) allocate(dhpb(maxdim)) +! if(.not.allocated(forcon)) allocate(forcon(maxdim)) + call flush(iout) do i=1,nfrag_ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup @@ -535,10 +585,29 @@ endif enddo do i=1,ndist_ - read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) +! read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) +! if (forcon(nhpb+1).gt.0.0d0) then +! nhpb=nhpb+1 +! dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + if (constr_dist.eq.11) then + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & + ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) + fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) + else +!C print *,"in else" + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & + ibecarb(i),forcon(nhpb+1) + endif if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 - dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + if (ibecarb(i).gt.0) then + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + endif + if (dhpb(nhpb).eq.0.0d0) & + dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + endif + #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",& @@ -547,12 +616,60 @@ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif - endif enddo + 1712 continue call flush(iout) return end subroutine read_dist_constr !----------------------------------------------------------------------------- + subroutine gen_dist_constr2 + use MPI_data + ! use control + use geometry, only: dist + use geometry_data + use control_data + use energy_data + integer :: i,j + real(kind=8) :: distance + if (constr_dist.eq.11) then + do i=nstart_sup,nstart_sup+nsup-1 + do j=i+2,nstart_sup+nsup-1 + distance=dist(i,j) + if (distance.le.15.0) then + nhpb=nhpb+1 + ihpb(nhpb)=i+nstart_seq-nstart_sup + jhpb(nhpb)=j+nstart_seq-nstart_sup + forcon(nhpb)=sqrt(0.04*distance) + fordepth(nhpb)=sqrt(40.0/distance) + dhpb(nhpb)=distance-0.1d0 + dhpb1(nhpb)=distance+0.1d0 + +#ifdef MPI + if (.not.out1file .or. me.eq.king) & + write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", & + nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) +#else + write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", & + nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) +#endif + endif + enddo + enddo + else + do i=nstart_sup,nstart_sup+nsup-1 + do j=i+2,nstart_sup+nsup-1 + nhpb=nhpb+1 + ihpb(nhpb)=i+nstart_seq-nstart_sup + jhpb(nhpb)=j+nstart_seq-nstart_sup + forcon(nhpb)=weidis + dhpb(nhpb)=dist(i,j) + enddo + enddo + endif + return + end subroutine gen_dist_constr2 + +!----------------------------------------------------------------------------- #ifdef WINIFL subroutine flush(iu) return @@ -833,6 +950,7 @@ use geometry_data, only: c,nres use energy_data + ! use control use compare_data use MD_data ! implicit real*8 (a-h,o-z) @@ -849,7 +967,7 @@ character*(*) :: tytul character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/) integer,dimension(nres) :: ica !(maxres) - + integer iti1 !el local variables integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit real(kind=8) :: etot @@ -941,9 +1059,11 @@ ires=0 do i=nnt,nct iti=itype(i) + iti1=itype(i+1) + if ((iti.eq.ntyp1).and.(iti1.eq.ntyp1)) cycle if (iti.eq.ntyp1) then ichain=ichain+1 - ires=0 +! ires=0 write (iunit,'(a)') 'TER' else ires=ires+1 @@ -992,6 +1112,7 @@ ! format. use geometry_data, only: c use energy_data + ! use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' @@ -1081,10 +1202,14 @@ return end subroutine intout !----------------------------------------------------------------------------- +#ifdef CLUSTER + subroutine briefout(it,ener,free)!,plik) +#else subroutine briefout(it,ener) - +#endif use geometry_data use energy_data + ! use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' @@ -1098,10 +1223,9 @@ ! print '(a,i5)',intname,igeom !el local variables integer :: i,it - real(kind=8) :: ener -#ifdef WHAM_RUN - integer :: iii -#endif + real(kind=8) :: ener,free +! character(len=80) :: plik +! integer :: iii #if defined(AIX) || defined(PGI) open (igeom,file=intname,position='append') @@ -1109,13 +1233,20 @@ open (igeom,file=intname,access='append') #endif #ifdef WHAM_RUN - iii=igeom +! iii=igeom igeom=iout #endif + print *,nss IF (NSS.LE.9) THEN +#ifdef CLUSTER + WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS) + ELSE + WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8) +#else WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS) ELSE WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9) +#endif WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS) ENDIF ! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES) @@ -1226,5 +1357,79 @@ end subroutine reads #endif !----------------------------------------------------------------------------- +! permut.F +!----------------------------------------------------------------------------- + subroutine permut(isym) + + use geometry_data, only: tabperm +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.LOCAL' +! include 'COMMON.VAR' +! include 'COMMON.CHAIN' +! include 'COMMON.INTERACT' +! include 'COMMON.IOUNITS' +! include 'COMMON.GEO' +! include 'COMMON.NAMES' +! include 'COMMON.CONTROL' + + integer :: n,isym +! logical nextp +!el external nextp + integer,dimension(isym) :: a +! parameter(n=symetr) +!el local variables + integer :: kkk,i + + n=isym + if (n.eq.1) then + tabperm(1,1)=1 + return + endif + kkk=0 + do i=1,n + a(i)=i + enddo + 10 print *,(a(i),i=1,n) + kkk=kkk+1 + do i=1,n + tabperm(kkk,i)=a(i) +! write (iout,*) "tututu", kkk + enddo + if(nextp(n,a)) go to 10 + return + end subroutine permut +!----------------------------------------------------------------------------- + logical function nextp(n,a) + + integer :: n,i,j,k,t +! logical :: nextp + integer,dimension(n) :: a + i=n-1 + 10 if(a(i).lt.a(i+1)) go to 20 + i=i-1 + if(i.eq.0) go to 20 + go to 10 + 20 j=i+1 + k=n + 30 t=a(j) + a(j)=a(k) + a(k)=t + j=j+1 + k=k-1 + if(j.lt.k) go to 30 + j=i + if(j.ne.0) go to 40 + nextp=.false. + return + 40 j=j+1 + if(a(j).lt.a(i)) go to 40 + t=a(i) + a(i)=a(j) + a(j)=t + nextp=.true. + return + end function nextp +!----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module io_base