X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_base.F90;h=f6058d7f3b5df2295afe767858ffb1150027fa7c;hb=refs%2Fheads%2FUCGM;hp=a8729a7c9046709e92419fad5e102077a3bc52f7;hpb=a18045ea1a2b2658ffe57821e33d4231012a77cf;p=unres4.git diff --git a/source/unres/io_base.F90 b/source/unres/io_base.F90 index a8729a7..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 @@ -958,7 +1016,7 @@ !----------------------------------------------------------------------------- subroutine pdbout(etot,tytul,iunit) - use geometry_data, only: c,nres + use geometry_data, only: c,nres,boxxsize,boxysize,boxzsize use energy_data ! use control use compare_data @@ -975,7 +1033,10 @@ ! include 'COMMON.MD' !el character(len=50) :: tytul character*(*) :: tytul - character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/) + character(len=1),dimension(58) :: chainid= (/'A','B','C','D','E','F','G','H','I','J',& + 'K','L','M','O','Q','P','R','S','T','U','W','X','Y','Z','a','b','c','d','e','f',& + 'g','h','i','j','k','l','m','n','o','p','r','s','t','u','w','x','y','z',& + '0','1','2','3','4','5','6','7','8','9'/) #ifdef XDRFPDB integer,dimension(maxres) :: ica !(maxres) #else @@ -984,13 +1045,14 @@ integer :: iti1 !el local variables integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit - real(kind=8) :: etot + real(kind=8) :: etot,xi,yi,zi integer :: nres2 nres2=2*nres #ifdef XDRFPDB if(.not.allocated(molnum)) allocate(molnum(maxres)) - molnum(:)=1 +! molnum(:)=mnumi(:) if(.not.allocated(vtot)) allocate(vtot(maxres*2)) !(maxres2) + if(.not.allocated(istype)) allocate(istype(maxres)) #else if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2) #endif @@ -1076,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 @@ -1091,6 +1155,7 @@ endif if (iti.eq.ntyp1_molec(molnum(i))) then ichain=ichain+1 + if (ichain.gt.58) ichain=1 ! ires=0 write (iunit,'(a)') 'TER' else @@ -1105,7 +1170,32 @@ if (istype(i).eq.0) istype(i)=1 write (iunit,40) iatom,sugartyp(istype(i)),restyp(iti,2), & chainid(ichain),ires,(c(j,i),j=1,3),vtot(i) + elseif (molnum(i).eq.4) 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) +! 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) + 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) else + write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),& ires,(c(j,i),j=1,3),vtot(i) endif @@ -1126,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) @@ -1248,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 @@ -1412,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) @@ -1489,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' @@ -1534,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 @@ -1562,6 +1645,29 @@ endif return end function sugarcode + integer function rescode_lip(res,ires) + character(len=3):: res + integer ires + rescode_lip=0 + if (res.eq.'NC3') rescode_lip =4 + if (res.eq.'NH3') rescode_lip =2 + if (res.eq.'PO4') rescode_lip =3 + if (res.eq.'GL1') rescode_lip =12 + if (res.eq.'GL2') rescode_lip =12 + if (res.eq.'C1A') rescode_lip =18 + if (res.eq.'C2A') rescode_lip =18 + if (res.eq.'C3A') rescode_lip =18 + if (res.eq.'C4A') rescode_lip =18 + if (res.eq.'C1B') rescode_lip =18 + if (res.eq.'C2B') rescode_lip =18 + if (res.eq.'C3B') rescode_lip =18 + if (res.eq.'C4B') rescode_lip =18 + if (res.eq.'D2A') rescode_lip =16 + if (res.eq.'D2B') rescode_lip =16 + + if (rescode_lip.eq.0) write(iout,*) "UNKNOWN RESIDUE",ires,res + return + end function !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module io_base