logical oneletter,iblnk
integer rescode
external rescode
+ integer lshift(maxchain),lmin(maxchain)
+ double precision shift_coord(3,maxchain),cm(3,maxchain),ccm(3),
+ & ccm_prev(3)
ifreq=1
- if (iargc().lt.5) then
+ if (iargc().lt.8) then
print '(2a)',
- & "Usage: xdrf2pdb-m one/three seqfile cxfile ntraj itraj",
- & " [pdbfile] [freq]"
+ & "Usage: xdrf2pdb-m one/three seqfile boxx boxy boxz cxfile ",
+ & "ntraj itraj [pdbfile] [freq]"
stop
endif
call getarg(1,onethree)
print *,nres
print '(20(a3,1x))',(sequenc(i),i=1,nres)
endif
+ call seq2chains(nres,itype,nchain,chain_length,chain_border,
+ & ireschain)
+ print *,"nchain",nchain
+ print *,"chain_border, chain_length"
+ nrestot=0
+ do i=1,nchain
+ nrestot=nrestot+chain_length(i)
+ print *,chain_border(1,i),chain_border(2,i),chain_length(i)
+ enddo
+ print *,"nrestot",nrestot
call getarg(3,arg)
+ read(arg,*) boxsize(1)
+ call getarg(4,arg)
+ read(arg,*) boxsize(2)
+ call getarg(5,arg)
+ read(arg,*) boxsize(3)
+ call getarg(6,arg)
iext = index(arg,'.cx') - 1
if (iext.lt.0) then
print *,"Error - not a cx file"
+ print *,arg
stop
endif
- call getarg(4,cntraj)
+ call getarg(7,cntraj)
read (cntraj,*) ntraj
- call getarg(5,citraj)
+ call getarg(8,citraj)
read (citraj,*) itraj
- if (iargc().gt.5) then
- call getarg(6,pdbfile)
+ if (iargc().gt.8) then
+ call getarg(9,pdbfile)
else
write(licz,'(bz,i3.3)') itraj
pdbfile=arg(:iext)//'_'//licz//'.pdb'
endif
- if (iargc().gt.6) then
- call getarg(7,cfreq)
+ if (iargc().gt.9) then
+ call getarg(10,cfreq)
read (cfreq,*) ifreq
endif
c print *,"ifreq",ifreq," ntraj",ntraj," itraj",itraj
open(9,file=pdbfile)
nnt = 1
- if (itype(1).eq.21) nnt = 2
+ if (itype(1).eq.ntyp1) nnt = 2
nct=nres
- if (itype(nres).eq.21) nct = nres-1
+ if (itype(nres).eq.ntyp1) nct = nres-1
print *,"nnt",nnt," nct",nct
call xdrfopen(ixdrf,arg, "r", iret)
kk = 0
kk = kk + 1
call xdrffloat(ixdrf, potE, iret)
call xdrffloat(ixdrf, uconst, iret)
-#ifdef NEWUNRES
- call xdrffloat(ixdrf, uconst_back, iret)
-#endif
+c#ifdef NEWUNRES
+c call xdrffloat(ixdrf, uconst_back, iret)
+c#endif
call xdrffloat(ixdrf, t_bath, iret)
call xdrfint(ixdrf, nss, iret)
do j=1,nss
call xdrfint(ixdrf, jhpb(j), iret)
enddo
call xdrfint(ixdrf, nfrag, iret)
+ call xdrfint(ixdrf, iset, iret)
do i=1,nfrag
call xdrffloat(ixdrf, qfrag(i), iret)
enddo
prec=10000.0
-
+c print *,'nss=',nss,'nfrag=',nfrag,'iset=',iset
isize=0
call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)
c(j,i+nres)=coord(j,ii+nres)
enddo
enddo
+c Bring all chains to same box and to a common origin
+ do ichain=1,nchain
+c print *,"ichain",ichain
+ do j=1,3
+ cm(j,ichain)=0.0d0
+ enddo
+ do i=chain_border(1,ichain),chain_border(2,ichain)
+c print *,i,c(:,i)
+ do j=1,3
+ cm(j,ichain)=cm(j,ichain)+c(j,i)
+ enddo
+ enddo
+ do j=1,3
+ cm(j,ichain)=cm(j,ichain)/chain_length(ichain)
+ enddo
+ enddo
+c write (*,*) "CM"
+c do i=1,nchain
+c write (*,'(i5,3f10.5)') i,(cm(j,i),j=1,3)
+c enddo
+ do j=1,3
+ lmin=0
+ sumodl_min=1.0d10
+ do i=0,3**nchain-1
+ ii=i
+ do ichain=1,nchain
+ lshift(ichain)=mod(ii,3)-1
+ ii=ii/3
+ enddo
+c print '(10i3)',(lshift(ichain),ichain=1,nchain)
+ sumodl=0.0d0
+ do ichain=1,nchain
+ do jchain=1,ichain-1
+ sumodl=sumodl+dabs(cm(j,ichain)+lshift(ichain)*boxsize(j)
+ & -cm(j,jchain)-lshift(jchain)*boxsize(j))
+ enddo
+ enddo
+c print *,"sumodl",sumodl," sumodl_min",sumodl_min
+ if (sumodl.lt.sumodl_min) then
+ sumodl_min=sumodl
+ lmin=lshift
+ endif
+ enddo
+ do ichain=1,nchain
+ cm(j,ichain)=cm(j,ichain)+lmin(ichain)*boxsize(j)
+ shift_coord(j,ichain)=lmin(ichain)*boxsize(j)
+ enddo
+ enddo
+ do ichain=1,nchain
+ do j=1,3
+ ccm(j)=ccm(j)+cm(j,ichain)*chain_length(ichain)
+ enddo
+ enddo
+ do j=1,3
+ ccm(j)=ccm(j)/nrestot
+ enddo
+c write (*,'(a,3f10.5)') "ccm ",(ccm(j),j=1,3)
+c write (*,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3)
+c write (*,'(a,3f10.5)') "diff ",(ccm(j)-ccm_prev(j),j=1,3)
+c write (*,*) "shift_coord"
+c do i=1,nchain
+c write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
+c enddo
+#ifdef PREVIOUS
+ if (previous) then
+ do j=1,3
+ change=.true.
+ do while (change)
+ if (ccm(j)-ccm_prev(j).gt.boxsize(j)/2) then
+ do ichain=1,nchain
+ shift_coord(j,ichain)=shift_coord(j,ichain)-boxsize(j)
+ enddo
+ ccm(j)=ccm(j)-boxsize(j)
+ else if (ccm(j)-ccm_prev(j).lt.-boxsize(j)/2) then
+ do ichain=1,nchain
+ shift_coord(j,ichain)=shift_coord(j,ichain)+boxsize(j)
+ enddo
+ ccm(j)=ccm(j)+boxsize(j)
+ else
+ change=.false.
+ endif
+ enddo
+ enddo
+ do j=1,3
+ ccm_prev(j)=ccm(j)
+ enddo
+ else
+ previous=.true.
+ do j=1,3
+ ccm_prev(j)=ccm(j)
+ enddo
+ endif
+#else
+ do ichain=1,nchain
+ do j=1,3
+ shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j)
+ enddo
+ enddo
+#endif
+c write (*,*) "shift_coord"
+c do i=1,nchain
+c write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
+c enddo
+ do ichain=1,nchain
+ do i=chain_border(1,ichain),chain_border(2,ichain)
+ do j=1,3
+ c(j,i)=c(j,i)+shift_coord(j,ichain)
+ c(j,i+nres)=c(j,i+nres)+shift_coord(j,ichain)
+ enddo
+ enddo
+ enddo
etot=potE
- write (tytul,'(a,i6)') "Structure",kk
+ write (tytul,'(a1,i6,a8,f6.1,a6,f10.1)') "#",kk," t_bath ",
+ & t_bath, " time ",time
call pdbout(etot,tytul,9)
endif
enddo