1 implicit real*8 (a-h,o-z)
4 include 'COMMON.INTERACT'
5 include 'COMMON.SBRIDGE'
7 real*4 prec,time,potE,uconst,t_bath,qfrag(100)
9 character*80 arg,seqfile,pdbfile
10 character*3 sequenc(maxres)
12 character*8 onethree,cfreq,cntraj,citraj
16 logical oneletter,iblnk
19 integer lshift(maxchain),lmin(maxchain)
20 double precision shift_coord(3,maxchain),cm(3,maxchain),ccm(3),
24 if (iargc().lt.8) then
26 & "Usage: xdrf2pdb-m one/three seqfile boxx boxy boxz cxfile ",
27 & "ntraj itraj [pdbfile] [freq]"
30 call getarg(1,onethree)
31 onethree = ucase(onethree)
32 if (onethree.eq.'ONE') then
34 else if (onethree.eq.'THREE') then
37 print *,"ONE or THREE must be specified"
39 call getarg(2,seqfile)
40 open (1,file=seqfile,status='old')
42 read(1,'(80a1)',end=10,err=10) (sequenc(i)(1:1),i=1,maxres)
46 do while (.not.iblnk(sequenc(i+1)(1:1)))
51 itype(i)=rescode(i,sequenc(i),1)
54 read(1,'(20(a3,1x))',end=11,err=11) (sequenc(i),i=1,maxres)
58 do while (.not.iblnk(sequenc(i+1)(1:1)))
63 itype(i)=rescode(i,sequenc(i),0)
66 print '(20(a3,1x))',(sequenc(i),i=1,nres)
68 call seq2chains(nres,itype,nchain,chain_length,chain_border,
70 print *,"nchain",nchain
71 print *,"chain_border, chain_length"
74 nrestot=nrestot+chain_length(i)
75 print *,chain_border(1,i),chain_border(2,i),chain_length(i)
77 print *,"nrestot",nrestot
79 read(arg,*) boxsize(1)
81 read(arg,*) boxsize(2)
83 read(arg,*) boxsize(3)
85 iext = index(arg,'.cx') - 1
87 print *,"Error - not a cx file"
94 if (iargc().gt.5) then
95 call getarg(6,pdbfile)
97 write(licz,'(bz,i3.3)') itraj
98 pdbfile=arg(:iext)//'_'//licz//'.pdb'
100 if (iargc().gt.6) then
104 c print *,"ifreq",ifreq," ntraj",ntraj," itraj",itraj
107 if (itype(1).eq.ntyp1) nnt = 2
109 if (itype(nres).eq.ntyp1) nct = nres-1
110 print *,"nnt",nnt," nct",nct
111 call xdrfopen(ixdrf,arg, "r", iret)
114 call xdrffloat(ixdrf, time, iret)
117 call xdrffloat(ixdrf, potE, iret)
118 call xdrffloat(ixdrf, uconst, iret)
120 c call xdrffloat(ixdrf, uconst_back, iret)
122 call xdrffloat(ixdrf, t_bath, iret)
123 call xdrfint(ixdrf, nss, iret)
125 call xdrfint(ixdrf, ihpb(j), iret)
126 call xdrfint(ixdrf, jhpb(j), iret)
128 call xdrfint(ixdrf, nfrag, iret)
129 call xdrfint(ixdrf, iset, iret)
131 call xdrffloat(ixdrf, qfrag(i), iret)
134 c print *,'nss=',nss,'nfrag=',nfrag,'iset=',iset
136 call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)
139 c write (*,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
140 c write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
141 c write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
142 c write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize)
143 if (mod(kk/ntraj,ifreq).eq.0 .and. mod(kk,ntraj).eq.itraj) then
144 if (isize .ne. nres+nct-nnt+1) then
145 print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1
156 c(j,i+nres)=coord(j,ii+nres)
159 c Bring all chains to same box and to a common origin
161 c print *,"ichain",ichain
165 do i=chain_border(1,ichain),chain_border(2,ichain)
168 cm(j,ichain)=cm(j,ichain)+c(j,i)
172 cm(j,ichain)=cm(j,ichain)/chain_length(ichain)
177 c write (*,'(i5,3f10.5)') i,(cm(j,i),j=1,3)
185 lshift(ichain)=mod(ii,3)-1
188 c print '(10i3)',(lshift(ichain),ichain=1,nchain)
192 sumodl=sumodl+dabs(cm(j,ichain)+lshift(ichain)*boxsize(j)
193 & -cm(j,jchain)-lshift(jchain)*boxsize(j))
196 c print *,"sumodl",sumodl," sumodl_min",sumodl_min
197 if (sumodl.lt.sumodl_min) then
203 cm(j,ichain)=cm(j,ichain)+lmin(ichain)*boxsize(j)
204 shift_coord(j,ichain)=lmin(ichain)*boxsize(j)
209 ccm(j)=ccm(j)+cm(j,ichain)*chain_length(ichain)
213 ccm(j)=ccm(j)/nrestot
215 c write (*,'(a,3f10.5)') "ccm ",(ccm(j),j=1,3)
216 c write (*,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3)
217 c write (*,'(a,3f10.5)') "diff ",(ccm(j)-ccm_prev(j),j=1,3)
218 c write (*,*) "shift_coord"
220 c write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
227 if (ccm(j)-ccm_prev(j).gt.boxsize(j)/2) then
229 shift_coord(j,ichain)=shift_coord(j,ichain)-boxsize(j)
231 ccm(j)=ccm(j)-boxsize(j)
232 else if (ccm(j)-ccm_prev(j).lt.-boxsize(j)/2) then
234 shift_coord(j,ichain)=shift_coord(j,ichain)+boxsize(j)
236 ccm(j)=ccm(j)+boxsize(j)
254 shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j)
258 c write (*,*) "shift_coord"
260 c write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
263 do i=chain_border(1,ichain),chain_border(2,ichain)
265 c(j,i)=c(j,i)+shift_coord(j,ichain)
266 c(j,i+nres)=c(j,i+nres)+shift_coord(j,ichain)
271 write (tytul,'(a1,i6,a8,f6.1,a6,f10.1)') "#",kk," t_bath ",
272 & t_bath, " time ",time
273 call pdbout(etot,tytul,9)