read2sigma in wham and cluster_wham
[unres.git] / source / wham / src / xread.F
1       subroutine xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
2       implicit none
3       include "DIMENSIONS"
4       include "DIMENSIONS.ZSCOPT"
5       include "DIMENSIONS.FREE"
6       integer MaxTraj
7       parameter (MaxTraj=2050)
8 #ifdef MPI
9       include "mpif.h"
10       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
11       include "COMMON.MPI"
12 #endif
13       include "COMMON.CHAIN"
14       include "COMMON.IOUNITS"
15       include "COMMON.PROTFILES"
16       include "COMMON.NAMES"
17       include "COMMON.VAR"
18       include "COMMON.GEO"
19       include "COMMON.ENEPS"
20       include "COMMON.PROT"
21       include "COMMON.INTERACT"
22       include "COMMON.FREE"
23       include "COMMON.SBRIDGE"
24       include "COMMON.OBCINKA"
25       real*4 csingle(3,maxres2)
26       character*64 nazwa
27       character*128 bprotfile_temp
28       integer i,j,k,l,ii,jj(maxslice),kk(maxslice),ll(maxslice),
29      &  mm(maxslice)
30       integer iscor,islice,islice1,slice
31       double precision energ
32       integer ilen,iroof
33       external ilen,iroof
34       double precision rmsdev,energia(0:max_ene),efree,eini,temp
35       double precision prop(maxQ)
36       integer ntot_all(0:maxprocs-1)
37       integer iparm,ib,iib,ir,nprop,nthr
38       double precision etot,time,ts(maxslice),te(maxslice)
39       integer is(maxslice),ie(maxslice),itraj,ntraj,it,iset
40       integer nstep(0:MaxTraj-1)
41       logical lerr
42
43       call set_slices(is,ie,ts,te,iR,ib,iparm)
44       do i=1,nQ
45         prop(i)=0.0d0
46       enddo
47       do i=0,MaxTraj-1
48         nstep(i)=0
49       enddo
50       ntraj=0
51       it=0
52       islice1=1
53       call opentmp(islice1,ientout,bprotfile_temp)
54       do while (.true.)
55         if (replica(iparm)) then
56           if (hamil_rep .or. umbrella(iparm)) then
57           read (ientin,*,end=1112,err=1112) time,eini,
58      &      etot,temp,nss,(ihpb(j),jhpb(j),j=1,nss),
59      &      nprop,(prop(j),j=1,nprop),iset
60           else
61           read (ientin,*,end=1112,err=1112) time,eini,
62      &      etot,temp,nss,(ihpb(j),jhpb(j),j=1,nss),
63      &      nprop,(prop(j),j=1,nprop)
64           endif
65           temp=1.0d0/(temp*1.987D-3)
66 c           write (iout,*) time,eini,etot,nss,
67 c     &     (ihpb(j),jhpb(j),j=1,nss),(prop(j),j=1,nprop)
68 c           call flush(iout)
69            do i=1,nT_h(iparm)
70              if (beta_h(i,iparm).eq.temp) then
71                iib = i
72                goto 22
73              endif
74            enddo
75   22       continue
76            if (i.gt.nT_h(iparm)) then
77              write (iout,*) "Error - temperature of conformation",
78      &       ii,1.0d0/(temp*1.987D-3),
79      &       " does not match any of the list"
80              write (iout,*)
81      &        1.0d0/(temp*1.987D-3),
82      &        (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
83              call flush(iout)
84              call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
85            endif
86         else
87            read (ientin,*,end=1112,err=1112) time,eini,
88      &       etot,nss,(ihpb(j),jhpb(j),j=1,nss),
89      &       nprop,(prop(j),j=1,nprop)
90              iib = ib
91         endif
92         itraj=mod(it,totraj(iR,iparm))
93 c        write (*,*) "ii",ii," itraj",itraj
94 c        call flush(iout)
95         it=it+1
96         if (itraj.gt.ntraj) ntraj=itraj
97         nstep(itraj)=nstep(itraj)+1
98         islice=slice(nstep(itraj),time,is,ie,ts,te)
99         read (ientin,'(8f10.5)',end=1112,err=1112)
100      &    ((csingle(l,k),l=1,3),k=1,nres),
101      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct)
102         efree=0.0d0
103         if (islice.gt.0 .and. islice.le.nslice) then
104         ii=ii+1
105         kk(islice)=kk(islice)+1
106         mm(islice)=mm(islice)+1
107         if (mod(nstep(itraj),isampl(iparm)).eq.0) then
108         jj(islice)=jj(islice)+1
109         if (hamil_rep) then
110           snk(iR,iib,iset,islice)=snk(iR,iib,iset,islice)+1
111         else if (umbrella(iparm)) then
112           snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
113         else
114           snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
115         endif
116         ll(islice)=ll(islice)+1
117 c         write (iout,*) ii,kk,jj,ll,eini,(prop(j),j=1,nprop)
118 #ifdef DEBUG
119 c        write (iout,*) "Writing conformation, record",ll(islice)
120 c        write (iout,*) "ib",ib," iib",iib
121          if (replica(iparm)) then 
122            write (iout,*) "TEMP",1.0d0/(temp*1.987D-3)
123            write (iout,*) "TEMP list"
124            write (iout,*) 
125      &      (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
126          endif
127          call flush(iout)
128 #endif
129 c         write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
130 c         write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
131 c         write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
132 c         call flush(iout)
133          if (islice.ne.islice1) then
134 c             write (iout,*) "islice",islice," islice1",islice1
135              close(ientout)
136 c             write (iout,*) "Closing file ",
137 c     &             bprotfile_temp(:ilen(bprotfile_temp))
138              call opentmp(islice,ientout,bprotfile_temp)
139 c             write (iout,*) "Opening file ",
140 c     &             bprotfile_temp(:ilen(bprotfile_temp))
141 c             call flush(iout)
142              islice1=islice
143          endif
144          write(ientout,rec=ll(islice)) 
145      &     ((csingle(l,k),l=1,3),k=1,nres),
146      &     ((csingle(l,k+nres),l=1,3),k=nnt,nct),
147      &     nss,(ihpb(k),jhpb(k),k=1,nss),
148      &     eini,efree,rmsdev,(prop(i),i=1,nQ),iR,iib,iparm
149 #ifdef DEBUG
150          do i=1,2*nres
151            do j=1,3
152              c(j,i)=csingle(j,i)
153            enddo
154          enddo
155          call int_from_cart1(.false.)
156          write (iout,*) "Writing conformation, record",ll(islice)
157          write (iout,*) "Cartesian coordinates"
158          write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
159          write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
160          write (iout,*) "Internal coordinates"
161          write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
162          write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
163          write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
164          write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
165          write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
166          write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
167          write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
168 c         write (iout,'(8f10.5)') (prop(j),j=1,nQ)
169          write (iout,'(16i5)') iscor
170          call flush(iout)
171 #endif
172          endif
173          endif
174        enddo
175  1112  continue
176        close(ientout)
177        write (iout,'(i10," trajectories found in file.")') ntraj+1
178        write (iout,'(a)') "Numbers of steps in trajectories:"
179        write (iout,'(8i10)') (nstep(i),i=0,ntraj)
180        write (iout,*) ii," conformations read from file",
181      &   nazwa(:ilen(nazwa))
182        write (iout,*) mm(islice)," conformations read so far, slice",
183      &    islice
184        write (iout,*) ll(islice)," conformations stored so far, slice",
185      &   islice
186        call flush(iout)
187        return
188        end