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