X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_base.F90;h=f6058d7f3b5df2295afe767858ffb1150027fa7c;hb=08b96e38d3e1534b82b14c0a8c7948e490c6cebb;hp=198fb1f0bd37b4b0f3e238047bf02a0cc4158262;hpb=0f5304c0d44e3ec4d6538de7d274c4b1a1930049;p=unres4.git diff --git a/source/unres/io_base.F90 b/source/unres/io_base.F90 index 198fb1f..f6058d7 100644 --- a/source/unres/io_base.F90 +++ b/source/unres/io_base.F90 @@ -5,7 +5,7 @@ implicit none !----------------------------------------------------------------------------- ! Max. number of AA residues - integer,parameter :: maxres=6500!1200 + integer,parameter :: maxres=101000!1200 ! Appr. max. number of interaction sites integer,parameter :: maxres2=2*maxres ! parameter (maxres6=6*maxres) @@ -155,9 +155,9 @@ !----------------------------------------------------------------------------- subroutine read_x(kanal,*) - use geometry_data - use energy_data - use geometry, only:int_from_cart1 + use geometry_data + use energy_data + use geometry, only:int_from_cart1 ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.GEO' @@ -435,7 +435,8 @@ use control_data, only:out1file use MPI_data character*320 afmcard - integer i + real, dimension(3) ::cbeg,cend + integer i,j print *, "wchodze" call card_concat(afmcard,.true.) call readi(afmcard,"BEG",afmbeg,0) @@ -443,15 +444,62 @@ call reada(afmcard,"FORCE",forceAFMconst,0.0d0) call reada(afmcard,"VEL",velAFMconst,0.0d0) print *,'FORCE=' ,forceAFMconst + cbeg=0.0d0 + cend=0.0d0 + if (afmbeg.eq.-1) then + read(inp,*) nbegafmmat,(afmbegcentr(i),i=1,nbegafmmat) + do i=1,nbegafmmat + do j=1,3 + cbeg(j)=cbeg(j)+c(j,afmbegcentr(i))/nbegafmmat + enddo + enddo + else + do j=1,3 + cbeg(j)=c(j,afmend) + enddo + endif + if (afmend.eq.-1) then + read(inp,*) nendafmmat,(afmendcentr(i),i=1,nendafmmat) + do i=1,nendafmmat + do j=1,3 + cend(j)=cend(j)+c(j,afmendcentr(i))/nendafmmat + enddo + enddo + else + cend(j)=c(j,afmend) + endif !------ NOW PROPERTIES FOR AFM distafminit=0.0d0 do i=1,3 - distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit + distafminit=(cend(i)-cbeg(i))**2+distafminit enddo distafminit=dsqrt(distafminit) print *,'initdist',distafminit return end subroutine read_afminp +!C------------------------------------------------------------ + subroutine read_afmnano + use geometry_data + use energy_data + use control_data, only:out1file + use MPI_data + integer :: i + real*8 :: cm + read(inp,*,err=11) inanomove,(inanotab(i),i=1,inanomove),forcenanoconst + cm=0.0d0 + do i=1,inanomove + cm=cm+c(3,inanotab(i)) + enddo + cm=cm/inanomove + distnanoinit=cm-tubecenter(3) + return +11 write(iout,*)& + "error in afmnano",& + ", number of center, their index and forceconstance" + stop + return + end subroutine read_afmnano + !----------------------------------------------------------------------------- subroutine read_dist_constr use MPI_data @@ -475,9 +523,9 @@ character(len=640) :: controlcard !el local variables - integer :: i,k,j,ddjk,ii,jj,itemp + integer :: i,k,j,ii,jj,itemp,mnumkk,mnumjj,k1,j1 integer :: nfrag_,npair_,ndist_ - real(kind=8) :: dist_cut + real(kind=8) :: dist_cut,ddjk ! write (iout,*) "Calling read_dist_constr" ! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup @@ -502,11 +550,11 @@ call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0) call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0) call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0) -! write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ + write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ ! write (iout,*) "IFRAG" -! do i=1,nfrag_ -! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) -! enddo + do i=1,nfrag_ + write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) + enddo ! write (iout,*) "IPAIR" ! do i=1,npair_ ! write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) @@ -518,15 +566,25 @@ call flush(iout) do i=1,nfrag_ - if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup - if (ifrag_(2,i).gt.nstart_sup+nsup-1) & - ifrag_(2,i)=nstart_sup+nsup-1 +! if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup +! if (ifrag_(2,i).gt.nstart_sup+nsup-1) & +! ifrag_(2,i)=nstart_sup+nsup-1 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) call flush(iout) if (wfrag_(i).gt.0.0d0) then do j=ifrag_(1,i),ifrag_(2,i)-1 do k=j+1,ifrag_(2,i) -! write (iout,*) "j",j," k",k + write (iout,*) "j",j," k",k,nres + j1=j + k1=k + if (j.gt.nres) j1=j-nres + if (k.gt.nres) k1=k-nres + mnumkk=molnum(k1) + mnumjj=molnum(j1) + + if ((itype(k1,mnumkk).gt.ntyp_molec(mnumkk)).or.& + (itype(j1,mnumjj).gt.ntyp_molec(mnumjj))) cycle + write (iout,*) "j",j," k",k,itype(k1,mnumkk),itype(j1,mnumjj) ddjk=dist(j,k) if (constr_dist.eq.1) then nhpb=nhpb+1 @@ -960,7 +1018,6 @@ use geometry_data, only: c,nres,boxxsize,boxysize,boxzsize use energy_data -! use energy, only: to_box ! use control use compare_data use MD_data @@ -1081,8 +1138,10 @@ iatom=0 ichain=1 - ires=0 - do i=nnt,nct + ires=0 + write(iout,*) "TUTUT" + do i=nnt,nct +! write(iout,*), "coord",c(1,i),c(2,i),c(3,i) iti=itype(i,molnum(i)) print *,i,molnum(i) if (molnum(i+1).eq.0) then @@ -1115,6 +1174,18 @@ xi=c(1,i) yi=c(2,i) zi=c(3,i) +! xi=dmod(xi,boxxsize) +! if (xi.lt.0.0d0) xi=xi+boxxsize +! yi=dmod(yi,boxysize) +! if (yi.lt.0.0d0) yi=yi+boxysize +! zi=dmod(zi,boxzsize) +! if (zi.lt.0.0d0) zi=zi+boxzsize + write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),& + ires,xi,yi,zi,vtot(i) + elseif (molnum(i).eq.5) then + xi=c(1,i) + yi=c(2,i) + zi=c(3,i) xi=dmod(xi,boxxsize) if (xi.lt.0.0d0) xi=xi+boxxsize yi=dmod(yi,boxysize) @@ -1123,7 +1194,8 @@ if (zi.lt.0.0d0) zi=zi+boxzsize write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),& ires,xi,yi,zi,vtot(i) - else + else + write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),& ires,(c(j,i),j=1,3),vtot(i) endif @@ -1144,6 +1216,7 @@ enddo write (iunit,'(a)') 'TER' do i=nnt,nct-1 + if (molnum(i).eq.5) cycle if (itype(i,1).eq.ntyp1) cycle if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then write (iunit,30) ica(i),ica(i+1) @@ -1266,7 +1339,7 @@ do i=1,nres iti=itype(i,1) ! print *,vbld(i),"vbld(i)",i - write (iout,'(a3,i4,6f10.3)') restyp(iti,1),i,vbld(i),& + write (iout,'(a3,i6,6f10.3)') restyp(iti,1),i,vbld(i),& rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),& rad2deg*omeg(i) enddo @@ -1430,46 +1503,39 @@ !----------------------------------------------------------------------------- ! 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 + subroutine permut(isym,nperm,tabperm) +!c integer maxperm,maxsym +!c parameter (maxperm=3628800) +!c parameter (maxsym=10) +! include "DIMENSIONS" + integer n,a,tabperm,nperm,kkk,i,isym ! logical nextp -!el external nextp - integer,dimension(isym) :: a -! parameter(n=symetr) -!el local variables - integer :: kkk,i - +! external nextp + dimension a(isym),tabperm(50,5040) n=isym + nperm=1 if (n.eq.1) then tabperm(1,1)=1 return endif + do i=2,n + nperm=nperm*i + enddo kkk=0 do i=1,n a(i)=i enddo - 10 print *,(a(i),i=1,n) + 10 continue +!c print '(i3,2x,100i3)',kkk+1,(a(i),i=1,n) kkk=kkk+1 do i=1,n - tabperm(kkk,i)=a(i) -! write (iout,*) "tututu", kkk + tabperm(i,kkk)=a(i) enddo if(nextp(n,a)) go to 10 return - end subroutine permut + end subroutine + + !----------------------------------------------------------------------------- logical function nextp(n,a) @@ -1507,7 +1573,6 @@ !----------------------------------------------------------------------------- integer function rescode(iseq,nam,itype,molecule) -! use io_base, only: ucase ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.NAMES' @@ -1552,7 +1617,7 @@ write(iout,*) "Explicit LIPID not yet implemented" stop else if (molecule.eq.5) then - do i=1,ntyp1_molec(molecule) + do i=-1,ntyp1_molec(molecule) print *,restyp(i,molecule) print *,i,restyp(i,molecule)(1:2),nam(1:2) if (ucase(nam(1:2)).eq.restyp(i,molecule)(1:2)) then