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