Fixed the "***" issue in PDB output format by shifting the center of each molecule...
[unres.git] / source / wham / src / cxread.F.org
1       subroutine cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,*)
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'DIMENSIONS.ZSCOPT'
5       include 'DIMENSIONS.FREE'
6       integer MaxTraj
7       parameter (MaxTraj=2050)
8       include 'COMMON.CONTROL'
9       include 'COMMON.CHAIN'
10       include 'COMMON.INTERACT'
11       include 'COMMON.NAMES'
12       include 'COMMON.IOUNITS'
13       include 'COMMON.HEADER'
14       include 'COMMON.SBRIDGE'
15       include 'COMMON.PROTFILES'
16       include 'COMMON.OBCINKA'
17       include 'COMMON.FREE'
18       include 'COMMON.VAR'
19       include 'COMMON.GEO'
20       include 'COMMON.PROT'
21       character*64 nazwa,bprotfile_temp
22       real*4 rtime,rpotE,ruconst,rt_bath,rprop(maxQ)
23       double precision time
24       integer iret,itmp,itraj,ntraj
25       real xoord(3,maxres2+2),prec
26       integer nstep(0:MaxTraj-1)
27       integer ilen
28       external ilen
29       integer ii,jj(maxslice),kk(maxslice),ll(maxslice),mm(maxslice)
30       integer is(MaxSlice),ie(MaxSlice),nrec_slice
31       double precision ts(MaxSlice),te(MaxSlice),time_slice
32       integer slice
33       logical conf_check
34       call set_slices(is,ie,ts,te,iR,ib,iparm)
35
36       do i=1,nQ
37         rprop(i)=0.0d0
38       enddo
39       do i=0,MaxTraj-1
40         nstep(i)=0
41       enddo
42       ntraj=0
43       it=0
44       iret=1
45 #if (defined(AIX) && !defined(JUBL))
46       call xdrfopen_(ixdrf,nazwa, "r", iret)
47 #else
48       call xdrfopen(ixdrf,nazwa, "r", iret)
49 #endif
50       if (iret.eq.0) return1
51
52       islice1=1
53       call opentmp(islice1,ientout,bprotfile_temp)
54 c      print *,"bumbum"
55       do while (iret.gt.0) 
56
57 #if (defined(AIX) && !defined(JUBL))
58 #ifdef DEBUG
59       write (iout,*) "ii",ii," itraj",itraj," it",it
60 #endif
61       call xdrffloat_(ixdrf, rtime, iret)
62       call xdrffloat_(ixdrf, rpotE, iret)
63 #ifdef DEBUG
64       write (iout,*) "rtime",rtime," rpotE",rpotE," iret",iret
65 #endif
66       call flush(iout)
67       call xdrffloat_(ixdrf, ruconst, iret)
68       call xdrffloat_(ixdrf, rt_bath, iret)
69       call xdrfint_(ixdrf, nss, iret)
70 #ifdef DEBUG
71       write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
72 #endif
73       do j=1,nss
74        if (dyn_ss) then
75         call xdrfint_(ixdrf, idssb(j), iret)
76         call xdrfint_(ixdrf, jdssb(j), iret)
77        idssb(j)=idssb(j)-nres
78        jdssb(j)=jdssb(j)-nres
79        else
80         call xdrfint_(ixdrf, ihpb(j), iret)
81         call xdrfint_(ixdrf, jhpb(j), iret)
82        endif
83       enddo
84       call xdrfint_(ixdrf, nprop, iret)
85       if (umbrella(iparm) .or. homol_nset.gt.1 .or. read_iset(iparm) 
86      &  .or. hamil_rep) 
87      &  call xdrfint(ixdrf, iset, iret)
88       do i=1,nprop
89         call xdrffloat_(ixdrf, rprop(i), iret)
90       enddo
91 #else
92 #ifdef DEBUG
93       write (iout,*) "ii",ii," itraj",itraj," it",it
94 #endif
95       call xdrffloat(ixdrf, rtime, iret)
96       call xdrffloat(ixdrf, rpotE, iret)
97 #ifdef DEBUG
98       write (iout,*) "rtime",rtime," rpotE",rpotE," iret",iret
99 #endif
100       call flush(iout)
101       call xdrffloat(ixdrf, ruconst, iret)
102       call xdrffloat(ixdrf, rt_bath, iret)
103       call xdrfint(ixdrf, nss, iret)
104 #ifdef DEBUG
105       write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
106 #endif
107       do j=1,nss
108        if (dyn_ss) then
109         call xdrfint(ixdrf, idssb(j), iret)
110         call xdrfint(ixdrf, jdssb(j), iret)
111 cc        idssb(j)=idssb(j)-nres
112 cc        jdssb(j)=jdssb(j)-nres
113 cc        write(iout,*) idssb(j),jdssb(j)
114        else
115         call xdrfint(ixdrf, ihpb(j), iret)
116         call xdrfint(ixdrf, jhpb(j), iret)
117        endif
118       enddo
119       call xdrfint(ixdrf, nprop, iret)
120 c      write (iout,*) "nprop",nprop
121       if (it.gt.0 .and. nprop.ne.nprop_prev) then
122         write (iout,*) "Warning previous nprop",nprop_prev,
123      &   " current",nprop
124         nprop=nprop_prev
125       else
126         nprop_prev=nprop
127       endif
128       call flush(iout)
129       if (umbrella(iparm) .or. homol_nset.gt.1 .or. read_iset(iparm) 
130      &  .or. hamil_rep) 
131      &  call xdrfint(ixdrf, iset, iret)
132       do i=1,nprop
133         call xdrffloat(ixdrf, rprop(i), iret)
134       enddo
135 #endif
136       if (iret.eq.0) exit
137       itraj=mod(it,totraj(iR,iparm))
138       if (iset.eq.0) iset = 1
139       call flush(iout)
140       it=it+1
141       if (itraj.gt.ntraj) ntraj=itraj
142       nstep(itraj)=nstep(itraj)+1
143 c      rprop(2)=dsqrt(rprop(2))
144 c      rprop(3)=dsqrt(rprop(3))
145 #ifdef DEBUG
146        write (iout,*) "umbrella ",umbrella
147        write (iout,*) rtime,rpotE,rt_bath,nss,
148      &     (ihpb(j),jhpb(j),j=1,nss),(rprop(j),j=1,nprop)
149        write (iout,*) "nprop",nprop," iset",iset," myparm",myparm
150        call flush(iout)
151 #endif
152       prec=10000.0
153
154       itmp=0
155 #if (defined(AIX) && !defined(JUBL))
156       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
157 #else
158       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
159 #endif
160 #ifdef DEBUG
161       write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,itmp)
162 #endif
163       if (iret.eq.0) exit
164       if (itmp .ne. nres + nct - nnt + 1) then
165         write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
166         call flush(iout)
167         exit
168       endif
169
170       time=rtime
171 c      write (iout,*) "calling slice"
172 c      call flush(iout)
173       islice=slice(nstep(itraj),time,is,ie,ts,te)
174 c      write (iout,*) "islice",islice
175 c      call flush(iout)
176
177       do i=1,nres
178         do j=1,3
179           c(j,i)=xoord(j,i)
180         enddo
181       enddo
182       do i=1,nct-nnt+1
183         do j=1,3
184           c(j,i+nres+nnt-1)=xoord(j,i+nres)
185         enddo
186       enddo
187
188       if (islice.gt.0 .and. islice.le.nslice .and. (.not.separate_parset
189      &    .or. iset.eq.myparm)) then
190         ii=ii+1
191         kk(islice)=kk(islice)+1
192         mm(islice)=mm(islice)+1
193         if (mod(nstep(itraj),isampl(iparm)).eq.0 .and. 
194      &     conf_check(ll(islice)+1,1)) then
195           if (replica(iparm)) then
196              rt_bath=1.0d0/(rt_bath*1.987D-3)
197              do i=1,nT_h(iparm)
198                if (abs(real(beta_h(i,iparm))-rt_bath).lt.1.0e-4) then
199                  iib = i
200                  goto 22
201                endif
202              enddo
203   22         continue
204              if (i.gt.nT_h(iparm)) then
205                write (iout,*) "Error - temperature of conformation",
206      &         ii,1.0d0/(rt_bath*1.987D-3),
207      &         " does not match any of the list"
208                write (iout,*)
209      &          1.0d0/(rt_bath*1.987D-3),
210      &          (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
211                call flush(iout)
212 c               exit
213 c               call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
214                ii=ii-1
215                kk(islice)=kk(islice)-1
216                mm(islice)=mm(islice)-1
217                goto 112
218              endif
219           else
220             iib = ib
221           endif
222
223           efree=0.0d0
224           jj(islice)=jj(islice)+1
225           if (umbrella(iparm)) then
226             snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
227           else if (hamil_rep) then
228             snk(1,iib,iparm,islice)=snk(1,iib,iparm,islice)+1
229           else
230             snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
231           endif
232           ll(islice)=ll(islice)+1
233 #ifdef DEBUG
234           write (iout,*) "Writing conformation, record",ll(islice)
235           write (iout,*) "ib",ib," iib",iib
236           write (iout,*) "ntraj",ntraj," itraj",itraj,
237      &      " nstep",nstep(itraj)
238           write (iout,*) "pote",rpotE," time",rtime
239 c          if (replica(iparm)) then
240 c            write (iout,*) "TEMP",1.0d0/(rt_bath*1.987D-3)
241 c            write (iout,*) "TEMP list"
242 c            write (iout,*)
243 c     &       (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
244 c          endif
245           write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
246 c          write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
247 c          write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
248           call flush(iout)
249 #endif
250           if (islice.ne.islice1) then
251 c            write (iout,*) "islice",islice," islice1",islice1
252             close(ientout) 
253 c            write (iout,*) "Closing file ",
254 c     &          bprotfile_temp(:ilen(bprotfile_temp))
255             call opentmp(islice,ientout,bprotfile_temp)
256 c            write (iout,*) "Opening file ",
257 c     &          bprotfile_temp(:ilen(bprotfile_temp))
258             islice1=islice
259           endif
260           if (umbrella(iparm) .or. homol_nset.gt.1) then
261             write(ientout,rec=ll(islice))
262      &        ((xoord(l,k),l=1,3),k=1,nres),
263      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
264      &        nss,(ihpb(k),jhpb(k),k=1,nss),
265      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
266      &        iset,iib,iparm
267           else if (hamil_rep) then
268             write(ientout,rec=ll(islice))
269      &        ((xoord(l,k),l=1,3),k=1,nres),
270      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
271      &        nss,(ihpb(k),jhpb(k),k=1,nss),
272      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
273      &        iR,iib,iset
274           else
275             write(ientout,rec=ll(islice))
276      &        ((xoord(l,k),l=1,3),k=1,nres),
277      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
278      &        nss,(ihpb(k),jhpb(k),k=1,nss),
279      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
280      &        iR,iib,iparm
281           endif
282 #ifdef DEBUG
283           write (iout,*) " constr_homology",constr_homology,
284      &      " ll",ll(islice)," iset",iset
285           call int_from_cart1(.false.)
286           write (iout,*) "Writing conformation, record",ll(islice)
287           write (iout,*) "Cartesian coordinates"
288           write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
289           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
290           write (iout,*) "Internal coordinates"
291           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
292           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
293           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
294           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
295           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
296           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
297           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
298 c          write (iout,'(8f10.5)') (rprop(j),j=1,nQ)
299           write (iout,'(16i5)') iscor
300           call flush(iout)
301 #endif
302         endif 
303       endif
304
305   112 continue
306
307       enddo
308       close(ientout)
309 #if (defined(AIX) && !defined(JUBL))
310       call xdrfclose_(ixdrf, iret)
311 #else
312       call xdrfclose(ixdrf, iret)
313 #endif
314       write (iout,'(i10," trajectories found in file.")') ntraj+1
315       write (iout,'(a)') "Numbers of steps in trajectories:"
316       write (iout,'(8i10)') (nstep(i),i=0,ntraj)
317       write (iout,*) ii," conformations read from file",
318      &   nazwa(:ilen(nazwa))
319       do islice=1,nslice
320         write (iout,*) mm(islice)," conformations read so far, slice",
321      &    islice
322         write (iout,*) ll(islice),
323      &  " conformations stored so far, slice",islice
324       enddo
325       call flush(iout)
326       return
327       end