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