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