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