source/wham/src-NEWSC/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 <<<<<<< HEAD
   7 =======
   8 #ifdef MPI
   9       include "mpif.h"
  10       include "COMMON.MPI"
  11 #endif
  12 >>>>>>> e183793... Added src_MD-M-newcorr (Adasko's source) and src-NEWSC of WHAM (with Momo's SCSC potentials)
  13       integer MaxTraj
  14       parameter (MaxTraj=2050)
  15       include 'COMMON.CHAIN'
  16       include 'COMMON.INTERACT'
  17       include 'COMMON.NAMES'
  18       include 'COMMON.IOUNITS'
  19       include 'COMMON.HEADER'
  20       include 'COMMON.SBRIDGE'
  21       include 'COMMON.PROTFILES'
  22       include 'COMMON.OBCINKA'
  23       include 'COMMON.FREE'
  24       include 'COMMON.VAR'
  25       include 'COMMON.GEO'
  26       include 'COMMON.PROT'
  27       character*64 nazwa,bprotfile_temp
  28       real*4 rtime,rpotE,ruconst,rt_bath,rprop(maxQ)
  29       double precision time
  30       integer iret,itmp,itraj,ntraj
  31       real xoord(3,maxres2+2),prec
  32       integer nstep(0:MaxTraj-1)
  33       integer ilen
  34       external ilen
  35       integer ii,jj(maxslice),kk(maxslice),ll(maxslice),mm(maxslice)
  36       integer is(MaxSlice),ie(MaxSlice),nrec_slice
  37       double precision ts(MaxSlice),te(MaxSlice),time_slice
  38       integer slice
  39       logical conf_check
  40 <<<<<<< HEAD
  41 =======
  42       character*4 lt_bath
  43       character*256 pdbfilename
  44       character*50 tytul
  45 >>>>>>> e183793... Added src_MD-M-newcorr (Adasko's source) and src-NEWSC of WHAM (with Momo's SCSC potentials)
  46       call set_slices(is,ie,ts,te,iR,ib,iparm)
  47
  48       do i=1,nQ
  49         rprop(i)=0.0d0
  50       enddo
  51       do i=0,MaxTraj-1
  52         nstep(i)=0
  53       enddo
  54       ntraj=0
  55       it=0
  56       iret=1
  57 #if (defined(AIX) && !defined(JUBL))
  58       call xdrfopen_(ixdrf,nazwa, "r", iret)
  59 #else
  60       call xdrfopen(ixdrf,nazwa, "r", iret)
  61 #endif
  62       if (iret.eq.0) return1
  63
  64       islice1=1
  65       call opentmp(islice1,ientout,bprotfile_temp)
  66 c      print *,"bumbum"
  67       do while (iret.gt.0)
  68
  69 #if (defined(AIX) && !defined(JUBL))
  70 #ifdef DEBUG
  71       write (iout,*) "ii",ii," itraj",itraj," it",it
  72 #endif
  73       call xdrffloat_(ixdrf, rtime, iret)
  74       call xdrffloat_(ixdrf, rpotE, iret)
  75 #ifdef DEBUG
  76       write (iout,*) "rtime",rtime," rpotE",rpotE," iret",iret
  77 #endif
  78       call flush(iout)
  79       call xdrffloat_(ixdrf, ruconst, iret)
  80       call xdrffloat_(ixdrf, rt_bath, iret)
  81       call xdrfint_(ixdrf, nss, iret)
  82 #ifdef DEBUG
  83       write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
  84 #endif
  85       do j=1,nss
  86         call xdrfint_(ixdrf, ihpb(j), iret)
  87         call xdrfint_(ixdrf, jhpb(j), iret)
  88       enddo
  89       call xdrfint_(ixdrf, nprop, iret)
  90       if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep)
  91      &  call xdrfint(ixdrf, iset, iret)
  92       do i=1,nprop
  93         call xdrffloat_(ixdrf, rprop(i), iret)
  94       enddo
  95 #else
  96 #ifdef DEBUG
  97       write (iout,*) "ii",ii," itraj",itraj," it",it
  98 #endif
  99       call xdrffloat(ixdrf, rtime, iret)
 100       call xdrffloat(ixdrf, rpotE, iret)
 101 #ifdef DEBUG
 102       write (iout,*) "rtime",rtime," rpotE",rpotE," iret",iret
 103 #endif
 104       call flush(iout)
 105       call xdrffloat(ixdrf, ruconst, iret)
 106       call xdrffloat(ixdrf, rt_bath, iret)
 107       call xdrfint(ixdrf, nss, iret)
 108 #ifdef DEBUG
 109       write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
 110 #endif
 111       do j=1,nss
 112         call xdrfint(ixdrf, ihpb(j), iret)
 113         call xdrfint(ixdrf, jhpb(j), iret)
 114       enddo
 115       call xdrfint(ixdrf, nprop, iret)
 116 c      write (iout,*) "nprop",nprop
 117       if (it.gt.0 .and. nprop.ne.nprop_prev) then
 118         write (iout,*) "Warning previous nprop",nprop_prev,
 119      &   " current",nprop
 120         nprop=nprop_prev
 121       else
 122         nprop_prev=nprop
 123       endif
 124       call flush(iout)
 125       if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep)
 126      &  call xdrfint(ixdrf, iset, iret)
 127       do i=1,nprop
 128         call xdrffloat(ixdrf, rprop(i), iret)
 129       enddo
 130 #endif
 131       if (iret.eq.0) exit
 132       itraj=mod(it,totraj(iR,iparm))
 133       if (iset.eq.0) iset = 1
 134       call flush(iout)
 135       it=it+1
 136       if (itraj.gt.ntraj) ntraj=itraj
 137       nstep(itraj)=nstep(itraj)+1
 138 c      rprop(2)=dsqrt(rprop(2))
 139 c      rprop(3)=dsqrt(rprop(3))
 140 #ifdef DEBUG
 141        write (iout,*) "umbrella ",umbrella
 142        write (iout,*) rtime,rpotE,rt_bath,nss,
 143      &     (ihpb(j),jhpb(j),j=1,nss),(rprop(j),j=1,nprop)
 144        write (iout,*) "nprop",nprop," iset",iset," myparm",myparm
 145        call flush(iout)
 146 #endif
 147       prec=10000.0
 148
 149       itmp=0
 150 #if (defined(AIX) && !defined(JUBL))
 151       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 152 #else
 153       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 154 #endif
 155 #ifdef DEBUG
 156       write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,itmp)
 157 #endif
 158       if (iret.eq.0) exit
 159       if (itmp .ne. nres + nct - nnt + 1) then
 160         write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
 161         call flush(iout)
 162         exit
 163       endif
 164
 165       time=rtime
 166 c      write (iout,*) "calling slice"
 167 c      call flush(iout)
 168       islice=slice(nstep(itraj),time,is,ie,ts,te)
 169 c      write (iout,*) "islice",islice
 170 c      call flush(iout)
 171
 172       do i=1,nres
 173         do j=1,3
 174           c(j,i)=xoord(j,i)
 175         enddo
 176       enddo
 177       do i=1,nct-nnt+1
 178         do j=1,3
 179           c(j,i+nres+nnt-1)=xoord(j,i+nres)
 180         enddo
 181       enddo
 182
 183       if (islice.gt.0 .and. islice.le.nslice .and. (.not.separate_parset
 184      &    .or. iset.eq.myparm)) then
 185         ii=ii+1
 186         kk(islice)=kk(islice)+1
 187         mm(islice)=mm(islice)+1
 188         if (mod(nstep(itraj),isampl(iparm)).eq.0 .and.
 189      &     conf_check(ll(islice)+1,1)) then
 190           if (replica(iparm)) then
 191              rt_bath=1.0d0/(rt_bath*1.987D-3)
 192              do i=1,nT_h(iparm)
 193                if (abs(real(beta_h(i,iparm))-rt_bath).lt.1.0e-4) then
 194                  iib = i
 195                  goto 22
 196                endif
 197              enddo
 198   22         continue
 199              if (i.gt.nT_h(iparm)) then
 200                write (iout,*) "Error - temperature of conformation",
 201      &         ii,1.0d0/(rt_bath*1.987D-3),
 202      &         " does not match any of the list"
 203                write (iout,*)
 204      &          1.0d0/(rt_bath*1.987D-3),
 205      &          (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 206                call flush(iout)
 207 c               exit
 208 c               call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 209                ii=ii-1
 210                kk(islice)=kk(islice)-1
 211                mm(islice)=mm(islice)-1
 212                goto 112
 213              endif
 214           else
 215             iib = ib
 216           endif
 217
 218           efree=0.0d0
 219           jj(islice)=jj(islice)+1
 220           if (umbrella(iparm)) then
 221             snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
 222           else if (hamil_rep) then
 223             snk(1,iib,iparm,islice)=snk(1,iib,iparm,islice)+1
 224           else
 225             snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
 226           endif
 227           ll(islice)=ll(islice)+1
 228 #ifdef DEBUG
 229           write (iout,*) "Writing conformation, record",ll(islice)
 230           write (iout,*) "ib",ib," iib",iib
 231           write (iout,*) "ntraj",ntraj," itraj",itraj,
 232      &      " nstep",nstep(itraj)
 233           write (iout,*) "pote",rpotE," time",rtime
 234 c          if (replica(iparm)) then
 235 c            write (iout,*) "TEMP",1.0d0/(rt_bath*1.987D-3)
 236 c            write (iout,*) "TEMP list"
 237 c            write (iout,*)
 238 c     &       (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 239 c          endif
 240           write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 241 c          write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
 242 c          write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
 243           call flush(iout)
 244 #endif
 245           if (islice.ne.islice1) then
 246 c            write (iout,*) "islice",islice," islice1",islice1
 247             close(ientout)
 248 c            write (iout,*) "Closing file ",
 249 c     &          bprotfile_temp(:ilen(bprotfile_temp))
 250             call opentmp(islice,ientout,bprotfile_temp)
 251 c            write (iout,*) "Opening file ",
 252 c     &          bprotfile_temp(:ilen(bprotfile_temp))
 253             islice1=islice
 254           endif
 255           if (umbrella(iparm)) then
 256             write(ientout,rec=ll(islice))
 257      &        ((xoord(l,k),l=1,3),k=1,nres),
 258      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 259      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 260      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 261      &        iset,iib,iparm
 262           else if (hamil_rep) then
 263             write(ientout,rec=ll(islice))
 264      &        ((xoord(l,k),l=1,3),k=1,nres),
 265      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 266      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 267      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 268      &        iR,iib,iset
 269           else
 270             write(ientout,rec=ll(islice))
 271      &        ((xoord(l,k),l=1,3),k=1,nres),
 272      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 273      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 274      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 275      &        iR,iib,iparm
 276           endif
 277 <<<<<<< HEAD
 278 =======
 279 #ifdef PDBOUT
 280 #ifdef MPI
 281           if (me.eq.Master) then
 282 #endif
 283           write (iout,*) "PDBOUT"
 284           write (iout,*) "temperature",1.0d0/(rt_bath*1.987D-3)
 285           call flush(iout)
 286           write (lt_bath,'(f4.0)') 1.0d0/(rt_bath*1.987D-3)
 287           write (iout,*) "lt_bath ",lt_bath
 288           pdbfilename=prefix(:ilen(prefix))//"_"//lt_bath//"pdb"
 289           write (iout,*) "pdb ",pdbfilename
 290           call flush(iout)
 291           open(ipdb,file=pdbfilename,position="append")
 292 c          write (tytul,'("Conformation",i10," T=",f5.1)')
 293 c     &      kk(islice),rt_bath
 294           call pdbout(kk(islice),1.0d0/(rt_bath*1.987D-3),
 295      &      efree+0.0d0,rpotE+0.0d0,efree+0.0d0,rmsdev+0.0d0)
 296           close(ipdb)
 297 #ifdef MPI
 298           endif
 299 #endif
 300 #endif
 301 >>>>>>> e183793... Added src_MD-M-newcorr (Adasko's source) and src-NEWSC of WHAM (with Momo's SCSC potentials)
 302 #ifdef DEBUG
 303           call int_from_cart1(.false.)
 304           write (iout,*) "Writing conformation, record",ll(islice)
 305           write (iout,*) "Cartesian coordinates"
 306           write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 307           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 308           write (iout,*) "Internal coordinates"
 309           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 310           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 311           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 312           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 313           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 314           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 315           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 316 c          write (iout,'(8f10.5)') (rprop(j),j=1,nQ)
 317           write (iout,'(16i5)') iscor
 318           call flush(iout)
 319 #endif
 320         endif
 321       endif
 322
 323   112 continue
 324
 325       enddo
 326       close(ientout)
 327 #if (defined(AIX) && !defined(JUBL))
 328       call xdrfclose_(ixdrf, iret)
 329 #else
 330       call xdrfclose(ixdrf, iret)
 331 #endif
 332       write (iout,'(i10," trajectories found in file.")') ntraj+1
 333       write (iout,'(a)') "Numbers of steps in trajectories:"
 334       write (iout,'(8i10)') (nstep(i),i=0,ntraj)
 335       write (iout,*) ii," conformations read from file",
 336      &   nazwa(:ilen(nazwa))
 337       do islice=1,nslice
 338         write (iout,*) mm(islice)," conformations read so far, slice",
 339      &    islice
 340         write (iout,*) ll(islice),
 341      &  " conformations stored so far, slice",islice
 342       enddo
 343       call flush(iout)
 344       return
 345       end