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