source/maxlik/src_FPy/readrtns_MP.F

   1       subroutine read_general_data(*)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5 #ifdef MPI
   6       include "mpif.h"
   7       include "COMMON.MPI"
   8       integer ierror,kolor,klucz
   9 #endif
  10       include "COMMON.WEIGHTS"
  11       include "COMMON.NAMES"
  12       include "COMMON.VMCPAR"
  13       include "COMMON.TORSION"
  14       include "COMMON.INTERACT"
  15       include "COMMON.ENERGIES"
  16       include "COMMON.MINPAR"
  17       include "COMMON.IOUNITS"
  18       include "COMMON.TIME1"
  19       include "COMMON.PROTFILES"
  20       include "COMMON.CHAIN"
  21       include "COMMON.CLASSES"
  22       include "COMMON.THERMAL"
  23       include "COMMON.FFIELD"
  24       include "COMMON.OPTIM"
  25       include "COMMON.CONTROL"
  26       include "COMMON.SCCOR"
  27       include "COMMON.SPLITELE"
  28       character*800 controlcard
  29       integer i,j,k,l,ii,n_ene_found,ist1,iet1,ist2,iet2,ls,le
  30       integer ind,itype1,itype2,itypf,itypsc,itypp,itypt1,itypt2
  31       integer ilen,lenpot,lenpre
  32       external ilen
  33       character*3 liczba,liczba1
  34 #ifndef ISNAN
  35       external proc_proc
  36 #endif
  37       character*16 ucase
  38       character*16 key
  39       external ucase
  40       double precision xchuj,weitemp,weitemp_low,weitemp_up
  41       character*80 item(7)
  42       integer nitem
  43       integer rescode
  44
  45 c      write (iout,*) "Enter read_general_data"
  46 c      call flush(iout)
  47
  48       lenpot=ilen(pot)
  49       lenpre=ilen(prefix)
  50       do i=1,max_ene
  51         ename(i)=" "
  52       enddo
  53 C Read first record: seed and number of energy components
  54       call card_concat(controlcard,.true.)
  55 c      write (iout,*) "card1",controlcard
  56 c      call flush(iout)
  57       call readi(controlcard,"ISEED",iseed,0)
  58       if (iseed.eq.0) stop "Seed is zero!"
  59 c      print *,me," iseed",iseed
  60       call readi(controlcard,"NPARMSET",nparmset,1)
  61 c      print *,me," nparmset",nparmset
  62 #ifdef MPI
  63 c Split processor pool if multiple parameter sets are treated
  64       if (nparmset.eq.1) then
  65         WHAM_COMM = ALL_COMM
  66 c        print *,me," opening ",outname," opened"
  67         if (me.eq.Master)open(iout,file=outname,status='unknown')
  68         myparm=1
  69 c        print *,me," outname ",outname," opened"
  70       else
  71         if (nprocs.lt.nparmset) then
  72           if (me.eq.Master) write (*,*)
  73      & "*** Cannot split parameter sets for fewer processors than sets",
  74      &  nprocs,nparmset
  75           call MPI_Finalize(ierror)
  76           stop
  77         endif
  78 c        write (iout,*) "nparmset",nparmset
  79         nprocs = nprocs/nparmset
  80         kolor = me/nprocs
  81         klucz = mod(me,nprocs)
  82 c        write (*,*) "My old rank",me," kolor",kolor," klucz",klucz
  83         call MPI_Comm_split(ALL_COMM,kolor,klucz,WHAM_COMM,ierror)
  84         call MPI_Comm_size(WHAM_COMM,nprocs,ierror)
  85         call MPI_Comm_rank(WHAM_COMM,me,ierror)
  86 c        write (*,*) "My new rank",me," comm size",nprocs
  87 c        write (*,*) "ALL_COMM",ALL_COMM
  88 c     &   " WHAM_COMM",WHAM_COMM
  89         myparm=kolor+1
  90 c        write (*,*) "My parameter set is",myparm
  91         write(liczba,'(bz,i3.3)') me
  92         write(liczba1,'(bz,i3.3)') myparm
  93         outname=prefix(:lenpre)//'.out_par'//liczba1//'_'//
  94      &    pot(:lenpot)//liczba
  95         if (me.eq.Master)open(iout,file=outname,status='unknown')
  96       endif
  97 #endif
  98 c      print *,me," checkpoint 1"
  99       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 100       call readi(controlcard,'TORMODE',tor_mode,0)
 101 c      print *,me," tor_mode",tor_mode
 102       call readi(controlcard,"N_ENE",n_ene,max_ene)
 103 c      print *,"iseed",iseed," n_ene",n_ene
 104       call readi(controlcard,"NPARMSET",nparmset,1)
 105       geom_and_wham_weights =
 106      &   index(controlcard,"GEOM_AND_WHAM_WEIGHTS").gt.0
 107 c      write (iout,*) "GEOM_AND_WHAM_WEIGHTS",geom_and_wham_weights
 108       if (iseed.gt.0) iseed=-iseed
 109       call vrndst(iseed)
 110       out_newe=index(controlcard,"OUT_NEWE").gt.0
 111       unres_pdb = index(controlcard,"UNRES_PDB").gt.0
 112 c      write (iout,*) "UNRES_PDB ",unres_pdb
 113 c Energy calculation settings section
 114       call readi(controlcard,'TORMODE',tor_mode,0)
 115 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 116        call reada(controlcard,'BOXX',boxxsize,100.0d0)
 117        call reada(controlcard,'BOXY',boxysize,100.0d0)
 118        call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 119 c      print *,me," checkpoint 2"
 120 c Cutoff range for interactions
 121       call reada(controlcard,"R_CUT",r_cut,15.0d0)
 122       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 123       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 124       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 125       if (lipthick.gt.0.0d0) then
 126        bordliptop=(boxzsize+lipthick)/2.0
 127        bordlipbot=bordliptop-lipthick
 128 C      endif
 129       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)
 130      & .and.me.eq.Master)
 131      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 132       buflipbot=bordlipbot+lipbufthick
 133       bufliptop=bordliptop-lipbufthick
 134       if ((lipbufthick*2.0d0).gt.lipthick.and.me.eq.Master)
 135      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 136       endif
 137 c      write(iout,*) "bordliptop=",bordliptop
 138 c      write(iout,*) "bordlipbot=",bordlipbot
 139 c      write(iout,*) "bufliptop=",bufliptop
 140 c      write(iout,*) "buflipbot=",buflipbot
 141       call readi(controlcard,'SYM',symetr,1)
 142 c      print *,me," checkpoint 3"
 143 c      write (iout,*) "n_ene",n_ene
 144       call flush(iout)
 145 c Energy-term-weights section
 146       wname(4)="WCORRH"
 147 C Read third record: weights
 148       do i=1,max_ene
 149         ww(i)=ww0(i)
 150       enddo
 151 c      print *,me," checkpoint 4"
 152 C Read fourth record: masks
 153       call card_concat(controlcard,.true.)
 154 c      write (iout,*) "card2",controlcard
 155       do i=1,n_ene
 156         key = "MASK_"//wname(i)(2:ilen(wname(i)))
 157         call readi(controlcard,key(:ilen(key)),imask(i),0)
 158       enddo
 159 c      print *,me," checkpoint 5"
 160 C Read fifth record: lower limits of weights
 161       call card_concat(controlcard,.true.)
 162 c      write (iout,*) "card3",controlcard
 163       do i=1,n_ene
 164         key = "WLOW_"//wname(i)(2:ilen(wname(i)))
 165         call reada(controlcard,key(:ilen(key)),ww_low(i),ww(i))
 166       enddo
 167 C Read sixth record: upper limits of weights
 168       call card_concat(controlcard,.true.)
 169 c      write (iout,*) "card4",controlcard
 170       do i=1,n_ene
 171         key = "WUP_"//wname(i)(2:ilen(wname(i)))
 172         call reada(controlcard,key(:ilen(key)),ww_up(i),ww(i))
 173       enddo
 174 C Read seventh record: VMC parameters
 175       call card_concat(controlcard,.true.)
 176 c      write (iout,*) "card5",controlcard
 177       call readi(controlcard,"MODE",mode,3)
 178       call readi(controlcard,"MAXSTEP_SCAN",maxstep_scan,50)
 179       call reada(controlcard,"RSTEP_SCAN",step_scan,1.0d-1)
 180       call readi(controlcard,"READ_STAT",read_stat,3)
 181       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 182 c      init_ene = index(controlcard,"INIT_ENE").gt.0 .and. read_stat.gt.1
 183       init_ene = .true.
 184       call readi(controlcard,"NMCM",nmcm,0)
 185       call readi(controlcard,"MAXSCAN",maxscan,0)
 186       call readi(controlcard,"MAXMIN",maxmin,100)
 187       call readi(controlcard,"MAXFUN",maxfun,100)
 188       call reada(controlcard,"TOLF",tolf,1.0d-6)
 189       call reada(controlcard,"RTOLF",rtolf,1.0d-6)
 190       out_minim=0
 191       if (index(controlcard,"OUT_MINIM").gt.0) out_minim=iout
 192       print_ini=0
 193       if (index(controlcard,"PRINT_INI").gt.0) print_ini=1
 194       print_fin=0
 195       if (index(controlcard,"PRINT_FIN").gt.0) print_fin=1
 196       print_stat=0
 197       if (index(controlcard,"PRINT_STAT").gt.0) print_stat=1
 198       call reada(controlcard,"RSTIME",rstime,9.0d2)
 199       call reads(controlcard,"MINIMIZER",minimizer,"SUMSL")
 200       call readi(controlcard,"OPT_MODE",opt_mode,0)
 201       mod_other_params=index(controlcard,"OPTIMIZE_OTHER").gt.0
 202       if (read_stat.eq.0 .and. mod_other_params) then
 203         if (me.eq.Master)
 204      &    write (iout,*) "Error: only optimization of energy-term",
 205      &    " weights can be performed with READ_STAT=",read_stat
 206         call flush(iout)
 207         return1
 208       endif
 209       if (index(controlcard,"RESTART").gt.0) then
 210         restart_flag=.true.
 211       else
 212         restart_flag=.false.
 213       endif
 214 c      print *,me," checkpoint 6"
 215       return
 216       end
 217 c-----------------------------------------------------------------------
 218       subroutine read_optim_parm(*)
 219       implicit none
 220       include "DIMENSIONS"
 221       include "DIMENSIONS.ZSCOPT"
 222 #ifdef MPI
 223       include "mpif.h"
 224       include "COMMON.MPI"
 225       integer ierror,kolor,klucz
 226 #endif
 227       include "COMMON.WEIGHTS"
 228       include "COMMON.NAMES"
 229       include "COMMON.VMCPAR"
 230       include "COMMON.TORSION"
 231       include "COMMON.INTERACT"
 232       include "COMMON.ENERGIES"
 233       include "COMMON.MINPAR"
 234       include "COMMON.IOUNITS"
 235       include "COMMON.TIME1"
 236       include "COMMON.PROTFILES"
 237       include "COMMON.CHAIN"
 238       include "COMMON.CLASSES"
 239       include "COMMON.THERMAL"
 240       include "COMMON.FFIELD"
 241       include "COMMON.OPTIM"
 242       include "COMMON.CONTROL"
 243       include "COMMON.SCCOR"
 244       character*800 controlcard
 245       character*4 reskind
 246       integer i,j,k,l,ii,n_ene_found,ist1,iet1,ist2,iet2,ls,le
 247       integer ind,ind1,ind2,itype1,itype2,itypf,itypsc,itypp,
 248      &   itypt1,itypt2,masktemp,iblock,iaux
 249       integer ilen,lenpot,lenpre
 250       external ilen
 251       double precision aux
 252       character*3 liczba,liczba1
 253 #ifndef ISNAN
 254       external proc_proc
 255 #endif
 256       character*16 ucase
 257       character*16 key
 258       external ucase
 259       double precision xchuj,weitemp,weitemp_low,weitemp_up
 260       character*80 item(7)
 261       character*3 typf
 262       integer nitem
 263       integer rescode
 264
 265       lenpot=ilen(pot)
 266       lenpre=ilen(prefix)
 267       if (me.eq.Master)
 268      &  write (iout,*) "MOD_OTHER_PARAMS ",mod_other_params
 269       do iblock=1,2
 270       do itypt1=-ntyp,ntyp
 271         do itypt2=-ntyp,ntyp
 272           mask_tor(0,itypt1,itypt2,iblock)=0
 273           weitor(0,itypt1,itypt2,iblock)=1.0d0
 274           weitor_low(0,itypt1,itypt2,iblock)=1.0d0
 275           weitor_up(0,itypt1,itypt2,iblock)=1.0d0
 276         enddo
 277       enddo
 278       do itypt1=-ntyp,ntyp
 279         do itypt2=-ntyp,ntyp
 280           do l=0,3
 281             mask_tor(l,itypt1,itypt2,iblock)=0
 282             weitor(l,itypt1,itypt2,iblock)=1.0
 283             weitor_low(l,itypt1,itypt2,iblock)=1.0
 284             weitor_up(l,itypt1,itypt2,iblock)=1.0
 285           enddo
 286         enddo
 287       enddo
 288       enddo
 289 #ifdef DEBUG
 290           write (iout,*) "ntyp",ntyp
 291           do itypt1=-ntyp,ntyp
 292             do itypt2=-ntyp,ntyp
 293               write (iout,*) "itypt1",itypt1," itypt2",itypt2,
 294      &      "weitor",weitor(0,itypt1,itypt2,1),eitor(0,itypt1,itypt2,2)
 295             enddo
 296           enddo
 297 #endif
 298       if (mod_other_params) then
 299 c        write (*,*)
 300 c     &"Internal parameters of UNRES energy components will be optimized"
 301         call card_concat(controlcard,.true.)
 302 c        write (iout,*) "mod_side ",controlcard
 303         call flush(iout)
 304         mod_side = (index(controlcard,"MOD_SIDE").gt.0)
 305         if (mod_side) then
 306           nsingle_sc=0
 307           npair_sc=0
 308           call card_concat(controlcard,.true.)
 309           do while ( index(controlcard,"END").eq.0 )
 310             call split_string(controlcard,item(1),4,nitem)
 311             if (nitem.eq.1 .or. item(2)(:1).eq."*") then
 312               nsingle_sc=nsingle_sc+1
 313               ityp_ssc(nsingle_sc)=rescode(1,item(1),0)
 314               if (nitem.lt.3) then
 315                 epss_low(nsingle_sc)=0.02d0
 316               else
 317                 read (item(3),*) epss_low(nsingle_sc)
 318               endif
 319               if (nitem.lt.4) then
 320                 epss_up(nsingle_sc)=0.0d0
 321               else
 322                 read (item(4),*) epss_up(nsingle_sc)
 323               endif
 324             else
 325               npair_sc=npair_sc+1
 326               ityp_psc(1,npair_sc)=rescode(1,item(1),0)
 327               ityp_psc(2,npair_sc)=rescode(1,item(2),0)
 328               if (nitem.lt.3) then
 329                 epsp_low(npair_sc)=0.02d0
 330               else
 331                 read (item(3),*) epsp_low(npair_sc)
 332               endif
 333               if (nitem.lt.4) then
 334                 epsp_up(npair_sc)=0.0d0
 335               else
 336                 read (item(4),*) epsp_up(npair_sc)
 337               endif
 338             endif
 339             call card_concat(controlcard,.true.)
 340           enddo
 341           if (nsingle_sc+npair_sc.eq.0) mod_side=.false.
 342           call card_concat(controlcard,.true.)
 343         endif
 344         mod_side_other=
 345      &    index(controlcard,"MOD_SIDE_OTHER").gt.0
 346 c        write (iout,*) "mod_side_other ",controlcard,mod_side_other
 347         if (mod_side_other) then
 348           mod_side_other=.false.
 349           do i=1,ntyp
 350             do j=1,5
 351               mask_sigma(j,i)=0
 352             enddo
 353           enddo
 354           call card_concat(controlcard,.true.)
 355 c          write (iout,*) "mod_side_oter ",controlcard
 356           do while ( index(controlcard,"END").eq.0 )
 357            call reads(controlcard,"RESKIND",reskind,"   ")
 358            itypsc=rescode(1,reskind,0)
 359            if (itypsc.lt.1 .or. itypsc.gt.20) then
 360               if (me.eq.Master) then
 361               write (iout,*) "Error in SC optimization data;",
 362      &          " SC type must not be dummy, type is" ,restyp," ",itypsc
 363               write (*,*) "Error in SC optimization data;",
 364      &          " SC type must not be dummy, type is" ,restyp," ",itypsc
 365               endif
 366               return1
 367            endif
 368            call readi(controlcard,"MASK_SIGMA0",mask_sigma(1,itypsc),0)
 369            call readi(controlcard,"MASK_SIGMAII",mask_sigma(2,itypsc),0)
 370            call readi(controlcard,"MASK_CHIP",mask_sigma(3,itypsc),0)
 371            call readi(controlcard,"MASK_ALP",mask_sigma(4,itypsc),0)
 372            call reada(controlcard,"SIGMA0_LOW",sigma_low(1,itypsc),0.d0)
 373            call reada(controlcard,"SIGMA0_UP",sigma_up(1,itypsc),0.0d0)
 374            call reada(controlcard,"SIGMAII_LOW",sigma_low(2,itypsc),
 375      &        0.0d0)
 376            call reada(controlcard,"SIGMAII_UP",sigma_up(2,itypsc),0.0d0)
 377            call reada(controlcard,"CHIP_LOW",sigma_low(3,itypsc),0.d0)
 378            call reada(controlcard,"CHIP_UP",sigma_up(3,itypsc),0.0d0)
 379            call reada(controlcard,"ALP_LOW",sigma_low(4,itypsc),0.d0)
 380            call reada(controlcard,"ALP_UP",sigma_up(4,itypsc),0.0d0)
 381            do k=1,4
 382              if (sigma_low(k,itypsc).eq.0.0d0 .and.
 383      &           sigma_up(k,itypsc).eq.0.0d0) mask_sigma(k,itypsc)=0
 384            enddo
 385            do k=1,4
 386              mod_side_other=mod_side_other.or.mask_sigma(k,itypsc).gt.0
 387            enddo
 388            if (me.eq.Master)
 389      &      write (iout,'(a4,i3,4(i2,2f8.3))') reskind,itypsc,
 390      &      (mask_sigma(k,itypsc),
 391      &      sigma_low(k,itypsc),sigma_up(k,itypsc),k=1,4)
 392            call card_concat(controlcard,.true.)
 393 c           write (iout,*) "mod_side_oter ",controlcard
 394           enddo
 395           if (me.eq.Master) then
 396           write (iout,*) "Optimization flags of one-body SC parameters"
 397           do i=1,ntyp
 398            write (iout,'(a4,i3,4(i2,2f8.3))') restyp(i),i,
 399      &      (mask_sigma(k,i),sigma_low(k,i),sigma_up(k,i),k=1,4)
 400           enddo
 401           endif
 402           call card_concat(controlcard,.true.)
 403         endif
 404 c        write (iout,*) "mod_side_other ",mod_side_other
 405 c        write (iout,*) "mod_tor ",controlcard
 406 c        call flush(iout)
 407         mod_tor = index(controlcard,"MOD_TOR").gt.0
 408         if (mod_tor) then
 409           do iblock=1,2
 410           do i=-ntortyp,ntortyp
 411             do j=-ntortyp,ntortyp
 412               mask_tor(0,i,j,iblock)=0
 413               weitor(0,i,j,iblock)=1.0d0
 414               weitor_low(0,i,j,iblock)=0.0d0
 415               weitor_up(0,i,j,iblock)=2.0d0
 416             enddo
 417           enddo
 418           enddo
 419           call card_concat(controlcard,.true.)
 420 c          write (iout,*) controlcard
 421           do while ( index(controlcard,"END").eq.0 )
 422             call split_string(controlcard,item(1),7,nitem)
 423             if (nitem.lt.2) then
 424               if (me.eq.Master)
 425      &        write (*,*) "Error in torsional optimization data;",
 426      &            " must specify both residues and type"
 427               return1
 428             endif
 429             weitemp=1.0d0
 430             weitemp_low=0.0d0
 431             weitemp_up=2.0d0
 432             masktemp=1
 433             iblock=1
 434 c            write (iout,*) "item3 ",item(3)," item4 ",item(4),
 435 c     &        " item5",item(5)
 436 c            call flush(iout)
 437             if (nitem.gt.2) read(item(3),*) masktemp
 438             if (nitem.gt.3) read(item(4),*) weitemp
 439             if (nitem.gt.4) read(item(5),*) weitemp_low
 440             if (nitem.gt.5) read(item(6),*) weitemp_up
 441             if (nitem.gt.6) read(item(7),*) iblock
 442 c            write (iout,*) controlcard
 443 c            write (iout,*) item(1)," ",item(2),weitemp,
 444 c     &         weitemp_low,weitemp_up
 445             if (index(item(1),"*").gt.0) then
 446               ist1=1
 447               iet1=ntyp
 448             else
 449               ist1=itortyp(rescode(1,item(1),0))
 450               iet1=ist1
 451             endif
 452             if (index(item(2),"*").gt.0) then
 453               ist2=1
 454               iet2=ntyp
 455             else
 456               ist2=itortyp(rescode(1,item(2),0))
 457               iet2=ist2
 458             endif
 459 c            write (iout,*) "ist1",ist1," iet1",iet1," ist2",ist2,
 460 c     &         " iet2",iet2
 461 c            call flush(iout)
 462             do itypt1=ist1,iet1
 463               do itypt2=ist2,iet2
 464 c                write (iout,*) "itypt1",itypt1," itypt2",itypt2,
 465 c     &            weitemp,weitemp_low,weitemp_up
 466                 mask_tor(0,itypt1,itypt2,iblock)=masktemp
 467                 weitor(0,itypt1,itypt2,iblock)=weitemp
 468                 weitor_low(0,itypt1,itypt2,iblock)=weitemp_low
 469                 weitor_up(0,itypt1,itypt2,iblock)=weitemp_up
 470 c                write (iout,*) "itypt1",itypt1," itypt2",itypt2,
 471 c     &            mask_tor(0,itypt1,itypt2,iblock),
 472 c     &            weitor(0,itypt1,itypt2,iblock),
 473 c     &            weitor_low(0,itypt1,itypt2,iblock),
 474 c     &            weitor_up(0,itypt1,itypt2,iblock)
 475               enddo
 476             enddo
 477             call card_concat(controlcard,.true.)
 478 c            write (iout,*) controlcard
 479           enddo
 480 #ifdef NEWCORR
 481           if (tor_mode.gt.1) then
 482             if (me.eq.Master) write (iout,*) "TORMODE is",tor_mode,
 483      &    " torsional constants are computed from energy map expansion."
 484           endif
 485 #endif
 486 #ifdef DEBUG
 487           if (me.eq.Master)
 488      &    write (iout,*) "Initialized torsional parameters:"
 489           do iblock=1,2
 490           do itypt1=-ntortyp,ntortyp
 491             do itypt2=-ntortyp,ntortyp
 492               if (me.eq.Master)write(iout,'(3i5,3f10.5)') itypt1,itypt2,
 493      &         mask_tor(0,itypt1,itypt2,iblock),
 494      &         weitor(0,itypt1,itypt2,iblock),
 495      &         weitor_low(0,itypt1,itypt2,iblock),
 496      &         weitor_up(0,itypt1,itypt2,iblock)
 497             enddo
 498           enddo
 499           enddo
 500 #endif
 501           if (tor_mode.eq.1) then
 502 c Exchange indices because the residue order in new torsionals is reversed
 503           do iblock=1,2
 504           do itypt1=-ntortyp,ntortyp
 505             do itypt2=itypt1+1,ntortyp
 506               iaux=mask_tor(0,itypt1,itypt2,iblock)
 507               mask_tor(0,itypt1,itypt2,iblock)=
 508      &          mask_tor(0,itypt2,itypt1,iblock)
 509               mask_tor(0,itypt2,itypt1,iblock)=iaux
 510               aux=weitor(0,itypt1,itypt2,iblock)
 511               weitor(0,itypt1,itypt2,iblock)=
 512      &          weitor(0,itypt2,itypt1,iblock)
 513               weitor(0,itypt2,itypt1,iblock)=aux
 514               aux=weitor_low(0,itypt1,itypt2,iblock)
 515               weitor_low(0,itypt1,itypt2,iblock)=
 516      &          weitor_low(0,itypt2,itypt1,iblock)
 517               weitor_low(0,itypt2,itypt1,iblock)=aux
 518               aux=weitor_up(0,itypt1,itypt2,iblock)
 519               weitor_up(0,itypt1,itypt2,iblock)=
 520      &          weitor_up(0,itypt2,itypt1,iblock)
 521               weitor_up(0,itypt2,itypt1,iblock)=aux
 522             enddo
 523           enddo
 524           enddo
 525           endif
 526           call card_concat(controlcard,.true.)
 527         endif
 528 c        write (iout,*) "mod_sccor ",controlcard
 529         call flush(iout)
 530         mod_sccor = index(controlcard,"MOD_SCCOR").gt.0
 531         if (mod_sccor) then
 532           call card_concat(controlcard,.true.)
 533           do iblock=1,2
 534           do l=1,3
 535           do i=-nsccortyp,nsccortyp
 536             do j=-nsccortyp,nsccortyp
 537               mask_tor(l,i,j,iblock)=0
 538               weitor(l,i,j,iblock)=1.0d0
 539               weitor_low(l,i,j,iblock)=0.0d0
 540               weitor_up(l,i,j,iblock)=2.0d0
 541             enddo
 542           enddo
 543           enddo
 544           enddo
 545           do while ( index(controlcard,"END").eq.0 )
 546             call split_string(controlcard,item(1),7,nitem)
 547             if (nitem.lt.3) then
 548               if (me.eq.Master)
 549      &        write (*,*) "Error in torsional optimization data;",
 550      &            " must specify both residues and type"
 551               return1
 552             endif
 553             weitemp=1.0d0
 554             weitemp_low=0.0d0
 555             weitemp_up=0.0d0
 556             if (nitem.gt.3) read(item(4),*) masktemp
 557             if (nitem.gt.4) read(item(5),*) weitemp
 558             if (nitem.gt.5) read(item(6),*) weitemp_low
 559             if (nitem.gt.6) read(item(7),*) weitemp_up
 560             if (index(item(1),"*").gt.0) then
 561               ist1=1
 562               iet1=ntyp
 563             else
 564               ist1=isccortyp(rescode(1,item(1),0))
 565               iet1=ist1
 566             endif
 567             if (index(item(2),"*").gt.0) then
 568               ist2=1
 569               iet2=ntyp
 570             else
 571               ist2=isccortyp(rescode(1,item(2),0))
 572               iet2=ist2
 573             endif
 574             if (index(item(3),"*").gt.0) then
 575               ls=1
 576               le=3
 577             else
 578               read(item(3),*) ls
 579               le=ls
 580             endif
 581             do itypt1=ist1,iet1
 582               do itypt2=ist2,iet2
 583                 do l=ls,le
 584                   mask_tor(l,itypt1,itypt2,1)=masktemp
 585                   weitor(l,itypt1,itypt2,1)=weitemp
 586                   weitor_low(l,itypt1,itypt2,1)=weitemp_low
 587                   weitor_up(l,itypt1,itypt2,1)=weitemp_up
 588                 enddo
 589               enddo
 590             enddo
 591             call card_concat(controlcard,.true.)
 592           enddo
 593           call card_concat(controlcard,.true.)
 594         endif
 595 #define DEBUG
 596 #ifdef DEBUG
 597         if (me.eq.Master)
 598      &   write (iout,*) "Optimizable sidechain-torsional parameters:"
 599         do itypt1=-nsccortyp,nsccortyp
 600           do itypt2=-nsccortyp,nsccortyp
 601             do l=1,3
 602               if (mask_tor(l,itypt1,itypt2,1).gt.0 .and. me.eq.Master)
 603      &         write (iout,'(4i5,3f10.5)') itypt1,itypt2,l,
 604      &         mask_tor(l,itypt1,itypt2,1),weitor(l,itypt1,itypt2,1),
 605      &        weitor_low(l,itypt1,itypt2,1),weitor_up(l,itypt1,itypt2,1)
 606             enddo
 607           enddo
 608         enddo
 609 #endif
 610 #undef DEBUG
 611 c        write (iout,*) "mod_fourier ",controlcard
 612 c        call flush(iout)
 613         mod_fourier(nloctyp)=index(controlcard,"MOD_FOURIER")
 614 #ifdef NEWCORR
 615         if (mod_fourier(nloctyp).gt.0) then
 616           mod_fourier(nloctyp)=0
 617           do i=0,nloctyp-1
 618             mod_fourier(i)=0
 619             do ii=1,3
 620               do j=1,2
 621                 mask_bnew1(ii,j,i)=0
 622                 mask_bnew2(ii,j,i)=0
 623                 mask_ccnew(ii,j,i)=0
 624                 mask_ddnew(ii,j,i)=0
 625               enddo
 626             enddo
 627             do ii=1,2
 628               do j=1,2
 629                 do k=1,2
 630                   mask_eenew(ii,j,k,i)=0
 631                 enddo
 632               enddo
 633             enddo
 634             do ii=1,3
 635               mask_e0new(ii,i)=0
 636             enddo
 637           enddo
 638           call card_concat(controlcard,.true.)
 639           do while ( index(controlcard,"END").eq.0 )
 640 c            write(iout,*) controlcard
 641             call reads(controlcard,"TYPF",typf,"   ")
 642             itypf=rescode(1,typf,0)
 643 c            write (iout,*) "TYPF ",typf," itypf",itypf
 644             if (itypf.eq.0 .or. itypf.gt.ntyp) then
 645               if (me.eq.Master)
 646      &        write (*,*) "Error specifying fourier parms"
 647               stop
 648             endif
 649             itypf=itype2loc(itypf)
 650 c            write (iout,*) "local type",itypf
 651
 652             if (index(controlcard,"B1_LOW").gt.0) then
 653
 654             call readi(controlcard,"IND",ind,0)
 655             call readi(controlcard,"COEFF",ii,0)
 656             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 657               if (me.eq.Master) write (iout,*)
 658      &          "Error specifying B1, components undefined",ind,ii
 659               stop
 660             endif
 661             mask_bnew1(ii,ind,itypf)=1
 662             call reada(controlcard,"B1_LOW",bnew1_low(ii,ind,itypf),
 663      &        0.1d0)
 664             call reada(controlcard,"B1_UP",bnew1_up(ii,ind,itypf),
 665      &        0.0d0)
 666             mod_fourier(itypf)=mod_fourier(itypf)
 667      &        +mask_bnew1(ii,ind,itypf)
 668
 669             else if (index(controlcard,"B2_LOW").gt.0) then
 670
 671             mask_bnew2(ii,ind,itypf)=1
 672             call readi(controlcard,"IND",ind,0)
 673             call readi(controlcard,"COEFF",ii,0)
 674             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 675               if (me.eq.Master) write (iout,*)
 676      &          "Error specifying B2, components undefined",ind,ii
 677               stop
 678             endif
 679             call reada(controlcard,"B2_LOW",bnew2_low(ii,ind,itypf),
 680      &        0.1d0)
 681             call reada(controlcard,"B2_UP",bnew2_up(ii,ind,itypf),
 682      &        0.0d0)
 683             mod_fourier(itypf)=mod_fourier(itypf)
 684      &        +mask_bnew2(ii,ind,itypf)
 685
 686             else if (index(controlcard,"C_LOW").gt.0) then
 687
 688             call readi(controlcard,"IND",ind,0)
 689             call readi(controlcard,"COEFF",ii,0)
 690             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 691               if (me.eq.Master) write (iout,*)
 692      &          "Error specifying C, components undefined",ind,ii
 693               stop
 694             endif
 695             mask_ccnew(ii,ind,itypf)=1
 696             call reada(controlcard,"C_LOW",ccnew_low(ii,ind,itypf),.1d0)
 697             call reada(controlcard,"C_UP",ccnew_up(ii,ind,itypf),0.0d0)
 698             mod_fourier(itypf)=mod_fourier(itypf)
 699      &        +mask_ccnew(ii,ind,itypf)
 700
 701             else if (index(controlcard,"D_LOW").gt.0) then
 702
 703             call readi(controlcard,"IND",ind,0)
 704             call readi(controlcard,"COEFF",ii,0)
 705             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 706               if (me.eq.Master) write (iout,*)
 707      &          "Error specifying D, components undefined",ind,ii
 708               stop
 709             endif
 710             mask_ddnew(ii,ind,itypf)=1
 711             call reada(controlcard,"D_LOW",ddnew_low(ii,ind,itypf),
 712      &        0.1d0)
 713             call reada(controlcard,"D_UP",ddnew_up(ii,ind,itypf),
 714      &        0.0d0)
 715             mod_fourier(itypf)=mod_fourier(itypf)
 716      &        +mask_ddnew(ii,ind,itypf)
 717
 718             else if (index(controlcard,"E0_LOW").gt.0) then
 719
 720             call readi(controlcard,"COEFF",ii,0)
 721             if (ii.eq.0 .or. ii.gt.3) then
 722               if (me.eq.Master) write (iout,*)
 723      &          "Error specifying E0, components undefined",ii
 724               stop
 725             endif
 726             mask_e0new(ii,itypf)=1
 727             call reada(controlcard,"E0_LOW",e0new_low(ii,itypf),
 728      &        0.1d0)
 729             call reada(controlcard,"E0_UP",e0new_up(ii,itypf),
 730      &        0.0d0)
 731             mod_fourier(itypf)=mod_fourier(itypf)
 732      &        +mask_e0new(ii,itypf)
 733
 734             else if (index(controlcard,"E_LOW").gt.0) then
 735
 736             call readi(controlcard,"IND1",ind1,0)
 737             call readi(controlcard,"IND2",ind2,0)
 738             call readi(controlcard,"COEFF",ii,0)
 739             if (ind1.eq.0 .or. ind2.eq.0 .or. ii.eq.0 .or. ii.gt.2) then
 740               if (me.eq.Master) write (iout,*)
 741      &          "Error specifying E, components undefined",ind1,ind2,ii
 742               stop
 743             endif
 744             mask_eenew(ii,ind1,ind2,itypf)=1
 745             call reada(controlcard,"E_LOW",
 746      &        eenew_low(ii,ind1,ind2,itypf),0.1d0)
 747             call reada(controlcard,"E_UP",
 748      &        eenew_up(ii,ind1,ind2,itypf),0.0d0)
 749             mod_fourier(itypf)=mod_fourier(itypf)
 750      &        +mask_eenew(ii,ind1,ind2,itypf)
 751             endif
 752             call card_concat(controlcard,.true.)
 753           enddo
 754           call card_concat(controlcard,.true.)
 755         endif
 756         do itypf=0,nloctyp-1
 757 c          write (iout,*) "itypf",itypf," mod_fourier",mod_fourier(itypf)
 758           mod_fourier(nloctyp)=mod_fourier(nloctyp)+mod_fourier(itypf)
 759         enddo
 760         if (me.eq.Master)
 761      &   write (iout,*) "mod_fourier all",mod_fourier(nloctyp)
 762 #else
 763           call card_concat(controlcard,.true.)
 764           do while ( index(controlcard,"END").eq.0 )
 765             call readi(controlcard,"TYPF",typf,"   ")
 766             itypf=refcode(typf)
 767             if (itypf.eq.0 .or. itypf.gt.ntyp) then
 768               if (me.eq.Master)
 769      &         write (*,*) "Error specifying fourier parms"
 770               stop
 771             endif
 772             itypf=itype2loc(itypf)
 773             call readi(controlcard,"IND",ind,0)
 774             call reada(controlcard,"B_LOW",b_low(ind,itypf),0.1d0)
 775             call reada(controlcard,"B_UP",b_up(ind,itypf),0.0d0)
 776             mask_fourier(ind,itypf)=1
 777             mod_fourier(itypf)=mod_fourier(itypf)
 778      &         +mask_fourier(ind,itypf)
 779             mod_fourier(nloctyp)=mod_fourier(nloctyp)
 780      &         +mask_fourier(ind,itypf)
 781             call card_concat(controlcard,.true.)
 782           enddo
 783           call card_concat(controlcard,.true.)
 784         endif
 785 #endif
 786         do i=1,2
 787           do j=1,2
 788             do k=1,2
 789               mask_elec(i,j,k)=0
 790             enddo
 791           enddo
 792         enddo
 793         mod_elec=index(controlcard,"MOD_ELEC").gt.0
 794         if (mod_elec) then
 795           mod_elec=.false.
 796           call card_concat(controlcard,.true.)
 797           do while ( index(controlcard,"END").eq.0 )
 798             call readi(controlcard,"TYPE1",itype1,0)
 799             call readi(controlcard,"TYPE2",itype2,0)
 800             if (itype1.eq.0 .or. itype1.gt.2 .or. itype2.eq.0
 801      &           .or. itype2.gt.2) then
 802               if (me.eq.Master)
 803      &         write (iout,*) "Error specifying elec parms"
 804               stop
 805             endif
 806             if (index(controlcard,"EPP").gt.0) then
 807               mod_elec=.true.
 808               mask_elec(itype1,itype2,1)=1
 809               mask_elec(itype2,itype1,1)=1
 810               call reada(controlcard,"EPP_LOW",epp_low(itype1,itype2),
 811      &          0.1d0)
 812               epp_low(itype2,itype1)=epp_low(itype1,itype2)
 813               call reada(controlcard,"EPP_UP",epp_up(itype1,itype2),
 814      &          0.0d0)
 815               epp_up(itype2,itype1)=epp_up(itype1,itype2)
 816             endif
 817             if (index(controlcard,"RPP").gt.0) then
 818               mod_elec=.true.
 819               mask_elec(itype1,itype2,2)=1
 820               mask_elec(itype2,itype1,2)=1
 821               call reada(controlcard,"RPP_LOW",rpp_low(itype1,itype2),
 822      &          0.1d0)
 823               rpp_low(itype2,itype1)=rpp_low(itype1,itype2)
 824               call reada(controlcard,"RPP_UP",rpp_up(itype1,itype2),
 825      &          0.0d0)
 826               rpp_up(itype2,itype1)=rpp_up(itype1,itype2)
 827             endif
 828             if (index(controlcard,"ELPP6").gt.0) then
 829               mod_elec=.true.
 830               mask_elec(itype1,itype2,3)=1
 831               mask_elec(itype2,itype1,3)=1
 832               call reada(controlcard,"ELPP6_LOW",
 833      &          elpp6_low(itype1,itype2),0.1d0)
 834               elpp6_low(itype2,itype1)=elpp6_low(itype1,itype2)
 835               call reada(controlcard,"ELPP6_UP",
 836      &          elpp6_up(itype1,itype2),0.0d0)
 837               elpp6_up(itype2,itype1)=elpp6_up(itype1,itype2)
 838             endif
 839             if (index(controlcard,"ELPP3").gt.0) then
 840               mod_elec=.true.
 841               mask_elec(itype1,itype2,4)=1
 842               mask_elec(itype2,itype1,4)=1
 843               call reada(controlcard,"ELPP3_LOW",
 844      &          elpp3_low(itype1,itype2),0.1d0)
 845               elpp3_low(itype2,itype1)=elpp3_low(itype1,itype2)
 846               call reada(controlcard,"ELPP3_UP",
 847      &          elpp3_up(itype1,itype2),0.0d0)
 848               elpp3_up(itype2,itype1)=elpp3_up(itype1,itype2)
 849             endif
 850             call card_concat(controlcard,.true.)
 851           enddo
 852           call card_concat(controlcard,.true.)
 853         endif
 854         do i=1,20
 855           do j=1,2
 856             do k=1,2
 857               mask_scp(i,j,k)=0
 858             enddo
 859           enddo
 860         enddo
 861         mod_scp=index(controlcard,"MOD_SCP").gt.0
 862         if (mod_scp) then
 863           mod_scp=.false.
 864           call card_concat(controlcard,.true.)
 865           do while (index(controlcard,"END").eq.0)
 866             if (index(controlcard,"EPSCP").gt.0) then
 867               mod_scp=.true.
 868               call readi(controlcard,"ITYPSC",itypsc,0)
 869               call readi(controlcard,"ITYPP",itypp,0)
 870               if (itypsc.gt.20 .or. itypp.eq.0 .or. itypp.gt.2) then
 871                 if (me.eq.Master)
 872      &           write (iout,*) "Error specifying scp parms"
 873                 stop
 874               endif
 875               mask_scp(itypsc,itypp,1)=1
 876               call reada(controlcard,"EPSCP_LOW",
 877      &          epscp_low(itypsc,itypp),0.1d0)
 878               call reada(controlcard,"EPSCP_UP",
 879      &          epscp_up(itypsc,itypp),0.0d0)
 880             endif
 881             if (index(controlcard,"RSCP").gt.0) then
 882               mod_scp=.true.
 883               call readi(controlcard,"ITYPSC",itypsc,0)
 884               call readi(controlcard,"ITYPP",itypp,0)
 885               mask_scp(itypsc,itypp,2)=1
 886               call readi(controlcard,"ITYPSC",itypsc,0)
 887               call readi(controlcard,"ITYPP",itypp,0)
 888               if (itypsc.gt.20 .or. itypp.eq.0 .or. itypp.gt.2) then
 889                 if (me.eq.Master)
 890      &           write (iout,*) "Error specifying scp parms"
 891                 stop
 892               endif
 893               call reada(controlcard,"RSCP_LOW",
 894      &          rscp_low(itypsc,itypp),0.1d0)
 895 c              write(iout,*)itypsc,itypp,rscp_low(itypsc,itypp)
 896               call reada(controlcard,"RSCP_UP",
 897      &          rscp_up(itypsc,itypp),0.0d0)
 898 c              write(iout,*)itypsc,itypp,rscp_up(itypsc,itypp)
 899             endif
 900             call card_concat(controlcard,.true.)
 901           enddo
 902         endif
 903 c        write (iout,*) "END ",controlcard
 904         call flush(iout)
 905         mod_other_params=
 906      &   mod_fourier(nloctyp).gt.0 .or. mod_elec .or. mod_scp
 907      &     .or. mod_sccor
 908         if (read_stat.lt.2. .and. mod_other_params) then
 909           if (me.eq.Master)
 910      &    write (iout,*)"Error - only weights and sidechain parameters",
 911      &     " can be optimized with READ_STAT=",read_stat
 912           call flush(iout)
 913           return1
 914         endif
 915         init_ene = init_ene .or. read_stat.eq.2
 916       endif
 917 c      write (*,*) "MOD_OTHER_PARAMS ",mod_other_params
 918 c      write (*,*) "MOD_SIDE ",mod_side," MOD_FOURIER",
 919 c     & mod_fourier(nloctyp),
 920 c     & " MOD_ELEC ",mod_elec," MOD_SCP ",mod_scp
 921       init_ene=init_ene .or. mod_other_params
 922 c      write (iout,*) "init_ene",init_ene
 923 c      write (iout,*)
 924       call flush(iout)
 925       return
 926       end
 927 c-----------------------------------------------------------------------------
 928       subroutine print_general_data(*)
 929       implicit none
 930       include "DIMENSIONS"
 931       include "DIMENSIONS.ZSCOPT"
 932 #ifdef MPI
 933       include "mpif.h"
 934       include "COMMON.MPI"
 935       integer ierror,kolor,klucz
 936 #endif
 937       include "COMMON.WEIGHTS"
 938       include "COMMON.NAMES"
 939       include "COMMON.VMCPAR"
 940       include "COMMON.TORSION"
 941       include "COMMON.INTERACT"
 942       include "COMMON.ENERGIES"
 943       include "COMMON.MINPAR"
 944       include "COMMON.IOUNITS"
 945       include "COMMON.TIME1"
 946       include "COMMON.PROTFILES"
 947       include "COMMON.CHAIN"
 948       include "COMMON.CLASSES"
 949       include "COMMON.THERMAL"
 950       include "COMMON.FFIELD"
 951       include "COMMON.OPTIM"
 952       include "COMMON.CONTROL"
 953       include "COMMON.SCCOR"
 954       character*800 controlcard
 955       integer i,j,k,l,ii,n_ene_found,itypt,jtypt
 956       integer ind,itype1,itype2,itypf,itypsc,itypp
 957       integer ilen,lenpot,lenpre
 958       external ilen
 959       character*3 liczba,liczba1
 960       if (mode.eq.1) then
 961         if (me.eq.Master)
 962      &    write (iout,*) "Single-point target function evaluation"
 963       else if (mode.eq.2) then
 964         if (me.eq.Master)
 965      &    write (iout,*) "Grid search of the parameter space"
 966       else if (mode.eq.3) then
 967         if (me.eq.Master)
 968      &    write (iout,*) "Minimization of the target function"
 969       else
 970         if (me.eq.Master)
 971      &    write (iout,*) "Wrong MODE",mode,", should be 1, 2, or 3"
 972         call flush(iout)
 973         return1
 974       endif
 975       if (me.eq.Master) write (iout,*) "RESCALE_MODE",rescale_mode
 976       mod_other_params=index(controlcard,"OPTIMIZE_OTHER").gt.0
 977       if (read_stat.eq.0 .and. mod_other_params) then
 978         if (me.eq.Master) then
 979         write (iout,*) "Error: only optimization of energy-term",
 980      &    " weights can be performed with READ_STAT=",read_stat
 981         call flush(iout)
 982         endif
 983         return1
 984       endif
 985       if (mode.eq.3 .and. me.eq.Master) then
 986         write (iout,*) "Parameters of the SUMSL procedure:"
 987         write (iout,'(a,t15,i5)') "MAXMIN",maxmin
 988         write (iout,'(a,t15,i5)') "MAXFUN",maxfun
 989         write (iout,'(a,t15,e15.7)') "TOLF",tolf
 990         write (iout,'(a,t15,e15.7)') "RTOLF",rtolf
 991       endif
 992       if (mod_other_params) then
 993         if (me.eq.Master) write (iout,*)
 994      &"Internal parameters of UNRES energy components will be optimized"
 995         call card_concat(controlcard,.true.)
 996         mod_side = (index(controlcard,"MOD_SIDE").gt.0)
 997         if (mod_side .and. me.eq.Master) then
 998           write (iout,*) "SC epsilons to be optimized:"
 999           write (iout,*) "Single: eps(i,j)=eps(i,j)+(x(i)+x(j))/2"
1000           do i=1,nsingle_sc
1001             write (iout,*) restyp(ityp_ssc(i)),epss_low(i),epss_up(i)
1002           enddo
1003           write (iout,*) "Pairs: eps(i,j)=eps(i,j)+x(i,j):"
1004           do i=1,npair_sc
1005             write (iout,*) restyp(ityp_psc(1,i)),
1006      &       restyp(ityp_psc(2,i)),epsp_low(i),epsp_up(i)
1007           enddo
1008         endif
1009         if (mod_sccor .and. me.eq.Master) then
1010           write (iout,*)"Torsional weights/coefficients to be optimized"
1011           write(iout,'(a)') "TYP1 TYP2    WEIGHT",
1012      &     "   LOWER BOUND UPPER BOUND"
1013           do itypt=-nsccortyp,nsccortyp
1014             do jtypt=-nsccortyp,nsccortyp
1015               do l=1,3
1016               if (mask_tor(l,itypt,jtypt,1).gt.0) then
1017                 write(iout,'(3a4,3f10.5)')l,restyp(itypt),restyp(jtypt),
1018      &           weitor(l,itypt,jtypt,1),weitor_low(l,itypt,jtypt,1),
1019      &           weitor_up(l,itypt,jtypt,1)
1020               endif
1021               enddo
1022             enddo
1023           enddo
1024         endif
1025         if (mod_fourier(nloctyp).gt.0 .and. me.eq.Master) then
1026           write (iout,*) "Fourier coefficients to be optimized"
1027           do itypf=0,nloctyp-1
1028             if (mod_fourier(itypf).gt.0) then
1029               write (iout,'(3a,i2)') "Type ",restyp(iloctyp(itypf)),
1030      &           " number of coeffs",mod_fourier(itypf)
1031               write(iout,'(a)') "NAME  COEFF    LOWER BOUND UPPER BOUND"
1032 #ifdef NEWCORR
1033              do j=1,2
1034                do k=1,3
1035                  if (mask_bnew1(k,j,itypf).gt.0)
1036      &            write (iout,'(2hB1,2i2,f10.5)') k,j,bnew1(k,j,itypf),
1037      &             bnew1_low(k,j,itypf),bnew1_up(k,j,itypf)
1038                enddo
1039              enddo
1040              do j=1,2
1041                do k=1,3
1042                  if (mask_bnew2(k,j,itypf).gt.0)
1043      &            write (iout,'(2hB2,2i2,3f11.5)') k,j,bnew2(k,j,itypf),
1044      &             bnew2_low(k,j,itypf),bnew2_up(k,j,itypf)
1045                enddo
1046              enddo
1047              do j=1,2
1048                do k=1,3
1049                  if (mask_ccnew(k,j,itypf).gt.0)
1050      &            write (iout,'(2hCC,2i2,3f11.5)') k,j,ccnew(k,j,itypf),
1051      &             ccnew_low(k,j,itypf),ccnew_up(k,j,itypf)
1052                enddo
1053              enddo
1054              do j=1,2
1055                do k=1,3
1056                  if (mask_ddnew(k,j,itypf).gt.0)
1057      &            write (iout,'(2hDD,2i2,3f11.5)') k,j,ddnew(k,j,itypf),
1058      &             ddnew_low(k,j,itypf),ddnew_up(k,j,itypf)
1059                enddo
1060              enddo
1061              do j=1,2
1062                if (mask_e0new(j,itypf).gt.0)
1063      &            write (iout,'(2hE0,i2,3f11.5)') j,e0new(j,itypf),
1064      &             e0new_low(j,itypf),e0new_up(j,itypf)
1065              enddo
1066              do j=1,2
1067                do k=1,2
1068                  do l=1,2
1069                    if (mask_eenew(l,k,j,itypf).gt.0)
1070      &              write (iout,'(2hEE,3i2,3f11.5)')
1071      &               l,k,j,eenew(l,k,j,itypf),eenew_low(l,k,j,itypf),
1072      &               eenew_up(l,k,j,itypf)
1073                  enddo
1074                enddo
1075              enddo
1076 #else
1077               do i=2,13
1078                 if (mask_fourier(i,itypf).gt.0) then
1079                   write (iout,'(1hB,i2,3f11.5)')
1080      &            i,b(i,itypf),b_low(i,itypf),b_up(i,itypf)
1081                 endif
1082               enddo
1083 #endif
1084             endif
1085           enddo
1086         endif
1087         if (mod_elec .and. Me.eq.Master) then
1088           write (iout,*)
1089           write (iout,*) "Electrostatic parameters to be optimized"
1090           do itype1=1,2
1091             do itype2=1,itype1
1092               if (mask_elec(itype1,itype2,1).gt.0)
1093      &        write (iout,'(2i3," EPP",f8.3," EPP_LOW",f8.3,
1094      &          " EPP_UP",f8.3)')
1095      &          itype1,itype2,epp(itype1,itype2),epp_low(itype1,itype2),
1096      &          epp_up(itype1,itype2)
1097               if (mask_elec(itype1,itype2,2).gt.0)
1098      &        write (iout,'(2i3," RPP",f8.3," RPP_LOW",f8.3,
1099      &          " RPP_UP",f8.3)')
1100      &          itype1,itype2,rpp(itype1,itype2),rpp_low(itype1,itype2),
1101      &          rpp_up(itype1,itype2)
1102               if (mask_elec(itype1,itype2,3).gt.0)
1103      &        write (iout,'(2i3," ELPP6",f8.3," ELPP6_LOW",f8.3,
1104      &          " ELPP6_UP",f8.3)')
1105      &          itype1,itype2,elpp6(itype1,itype2),
1106      &          elpp6_low(itype1,itype2),elpp6_up(itype1,itype2)
1107               if (mask_elec(itype1,itype2,4).gt.0)
1108      &        write (iout,'(2i3," ELPP3",f8.3," ELPP3_LOW",f8.3,
1109      &          " ELPP3_UP",f8.3)')
1110      &          itype1,itype2,elpp3(itype1,itype2),
1111      &          elpp3_low(itype1,itype2),elpp3_up(itype1,itype2)
1112             enddo
1113           enddo
1114         endif
1115         if (mod_scp .and. me.eq.Master) then
1116           mod_scp=.false.
1117           write (iout,*)
1118           write (iout,*) "SCP parameters to be optimized:"
1119           if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
1120      &        mask_scp(0,2,2) .gt. 0) then
1121             write (iout,*)
1122      & "SCP parameters are assumed to depend on peptide group type only"
1123             do j=1,2
1124               if (mask_scp(0,j,1).gt.0)
1125      &          write (iout,'(i3," EPSCP",f8.3," EPSCP_LOW",f8.3,
1126      &           " EPSCP_UP",f8.3)') j,eps_scp(1,j),epscp_low(0,j),
1127      &           epscp_up(0,j)
1128               if (mask_scp(0,j,2).gt.0)
1129      &          write (iout,'(i3," RSCP",f8.3," RSCP_LOW",f8.3,
1130      &           " RSCP_UP",f8.3)') j,rscp(1,j),rscp_low(0,j),
1131      &           rscp_up(0,j)
1132             enddo
1133           else
1134             do i=1,20
1135               do j=1,2
1136                 if (mask_scp(i,j,1).gt.0)
1137      &            write (iout,'(2i3," EPSCP",f8.3," EPSCP_LOW",f8.3,
1138      &             " EPSCP_UP",f8.3)') i,j,eps_scp(i,j),epscp_low(i,j),
1139      &             epscp_up(i,j)
1140                 if (mask_scp(i,j,2).gt.0)
1141      &            write (iout,'(2i3," RSCP",f8.3," RSCP_LOW",f8.3,
1142      &             " RSCP_UP",f8.3)') i,j,rscp(i,j),rscp_low(i,j),
1143      &             rscp_up(i,j)
1144               enddo
1145             enddo
1146           endif
1147         endif
1148       endif
1149       if (me.eq.Master)
1150      &  write (iout,*) "MOD_OTHER_PARAMS ",mod_other_params
1151       if (me.eq.Master)
1152      &  write (iout,*) "MOD_SIDE ",mod_side," MOD_FOURIER",
1153      & mod_fourier(nloctyp),
1154      & " MOD_ELEC ",mod_elec," MOD_SCP ",mod_scp
1155       init_ene=init_ene .or. mod_other_params
1156       if (me.eq.Master) write (iout,*) "init_ene",init_ene
1157       if (me.eq.Master) write (iout,*)
1158       call flush(iout)
1159       return
1160       end
1161 c-----------------------------------------------------------------------------
1162       subroutine read_protein_data(*)
1163       implicit none
1164       include "DIMENSIONS"
1165       include "DIMENSIONS.ZSCOPT"
1166 #ifdef MPI
1167       include "mpif.h"
1168       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1169       include "COMMON.MPI"
1170 #endif
1171       include "COMMON.CONTROL"
1172       include "COMMON.CHAIN"
1173       include "COMMON.CLASSES"
1174       include "COMMON.COMPAR"
1175       include "COMMON.ENERGIES"
1176       include "COMMON.IOUNITS"
1177       include "COMMON.PROTFILES"
1178       include "COMMON.PROTNAME"
1179       include "COMMON.VMCPAR"
1180       include "COMMON.OPTIM"
1181       include "COMMON.WEIGHTS"
1182       include "COMMON.XBOUND"
1183       include "COMMON.NAMES"
1184       include "COMMON.ALLPROT"
1185       include "COMMON.THERMAL"
1186       include "COMMON.TIME1"
1187       include "COMMON.WEIGHTDER"
1188       character*64 nazwa,key
1189       character*16000 controlcard
1190       character*16000 all_protfiles
1191       integer i,j,k,kk,l,iprot,ii,if,ib,ibatch,nn,nn1,iww,maskcheck,
1192      &  il,inat,ilevel,weightread,jfrag,jclass,nfragm,iparm
1193       integer nrec,nlines,iscor
1194       double precision energ,wangnorm(maxT),wq(maxT),wrms(maxT),
1195      &  wrgy(maxT),wsign(maxT),wangnat(maxT),wqnat(maxT),wrmsnat(maxT),
1196      &  wrgynat(maxT),wcorangnorm(maxT),wcorq(maxT),
1197      &  wcorrms(maxT),wcorrgy(maxT),wcorsign(maxT),wcorangnat(maxT),
1198      &  wcorqnat(maxT),wcorrmsnat(maxT),wcorrgynat(maxT),
1199      &  angnormlow(maxT),qlow(maxT),rmslow(maxT),
1200      &  rgylow(maxT),signlow(maxT),angnatlow(maxT),
1201      &  qnatlow(maxT),rmsnatlow(maxT),rgynatlow(maxT),
1202      &  angcorlow(maxT),qcorlow(maxT),rmscorlow(maxT),rgycorlow(maxT),
1203      &  signcorlow(maxT),angcornatlow(maxT),qcornatlow(maxT),
1204      &  rmscornatlow(maxT),rgycornatlow(maxT),signcornatlow(maxT),
1205      &  angnormup(maxT),qup(maxT),rmsup(maxT),rgyup(maxT),signup(maxT),
1206      &  angnatup(maxT),qnatup(maxT),rmsnatup(maxT),rgynatup(maxT),
1207      &  angcorup(maxT),qcorup(maxT),rmscorup(maxT),rgycorup(maxT),
1208      &  signcorup(maxT),
1209      &  angcornatup(maxT),qcornatup(maxT),rmscornatup(maxT),
1210      &  rgycornatup(maxT),signcornatup(MaxT)
1211       integer ilen,iroof
1212       external ilen,iroof
1213       double precision ebia(maxprot),rjunk
1214       character*10 liczba
1215       character*64 zeros /
1216      &'0000000000000000000000000000000000000000000000000000000000000000'
1217      & /
1218       logical lerr
1219 c      print *,"Processor",me," calls read_protein_data"
1220
1221 C Read seventh record: general data of proteins used for calibration
1222       call card_concat(controlcard,.true.)
1223 c      write(2, *)controlcard
1224       call readi(controlcard,"NPROT",nprot,0)
1225       pdbref=index(controlcard,"PDBREF").gt.0
1226       print_refstr=me.eq.Master.and.
1227      & index(controlcard,"PRINT_REFSTR").gt.0
1228       if (nprot.eq.0) then
1229         if (me.eq.Master) then
1230         write(iout,*) "Error: at least one protein must be specified."
1231         call flush(iout)
1232         endif
1233         return1
1234       endif
1235       protname(0)="all"
1236       do iprot=1,nprot
1237         read (inp,'(a)') protname(iprot)
1238         if (me.eq.Master) then
1239         write (iout,*)
1240         write (iout,*) "Reading data of protein",iprot," named ",
1241      &   protname(iprot)
1242         endif
1243         call card_concat(controlcard,.true.)
1244         call reada(controlcard,"ENECUT_MIN",enecut_min(iprot),15.0d0)
1245         call reada(controlcard,"ENECUT_MAX",enecut_max(iprot),100.0d0)
1246         call reada(controlcard,"ENECUT",enecut(iprot),enecut_min(iprot))
1247         if (enecut(iprot).lt.enecut_min(iprot))
1248      &    enecut(iprot)=enecut_min(iprot)
1249         if (enecut_max(iprot).le.enecut_min(iprot))
1250      &    enecut_max(iprot)=2*enecut_min(iprot)
1251         if (me.eq.Master)
1252      &  write (iout,'(3(a,f9.1))') "ENECUT",enecut(iprot)," ENECUT_MIN",
1253      &   enecut_min(iprot)," ENECUT_MAX",enecut_max(iprot)
1254         call readi(controlcard,"HEFAC",hefac(iprot),50)
1255         call readi(controlcard,"HTFAC",htfac(iprot),50)
1256         call readi(controlcard,"HELOW",hemax_low(iprot),20)
1257         call readi(controlcard,"HTLOW",htmax_low(iprot),20)
1258         if (me.eq.Master) write (iout,*) "iprot",iprot,
1259      &   " hefac",hefac(iprot)," helow",hemax_low(iprot),
1260      &   " htfac",htfac(iprot)," htlow",htmax_low(iprot)
1261 c 7/27/2013 AL Read maximum likelihood data
1262         call card_concat(controlcard,.true.)
1263         call readi(controlcard,"NBETA_L",nbeta(1,iprot),0)
1264         if (me.eq.Master) write (iout,*) "NBETA_L",nbeta(1,iprot)
1265         caonly(iprot)=index(controlcard,"CAONLY").gt.0
1266         sconly(iprot)=index(controlcard,"SCONLY").gt.0
1267         rmscomp(iprot)=index(controlcard,"RMSCOMP").gt.0
1268         call reada(controlcard,"SIGMA",sigma2(iprot),4.0d0)
1269         if (me.eq.Master)
1270      &   write (iout,*) "RMSCOMP",rmscomp," SIGMA",sigma2(iprot),
1271      &     " CAONLY ",caonly(iprot)," SCONLY",sconly(iprot)
1272         do ib=1,nbeta(1,iprot)
1273           read(inp,*)betaT(ib,1,iprot),weilik(ib,iprot),
1274      &      pdbfile(ib,iprot)
1275           if (me.eq.Master)
1276      &    write (iout,*) i,betaT(ib,1,iprot),weilik(ib,iprot),
1277      &      pdbfile(ib,iprot)
1278         enddo
1279 c 7/27/2013 AL Read heat-capacity data
1280         call card_concat(controlcard,.true.)
1281         call readi(controlcard,"NBETA_CV",nbeta(2,iprot),0)
1282         call reada(controlcard,"WCV",wcv(iprot),1.0d0)
1283         call reada(controlcard,"BASE",heatbase(iprot),0.0d0)
1284         if (me.eq.Master)
1285      &   write (iout,*) "NBETA_CV",nbeta(2,iprot)," WCV",wcv(iprot)
1286         do ib=1,nbeta(2,iprot)
1287           read(inp,*) betaT(ib,2,iprot),target_cv(ib,iprot),
1288      &      weiCv(ib,iprot)
1289           weiCv(ib,iprot)=weiCv(ib,iprot)*wcv(iprot)
1290           if (me.eq.Master)
1291      &    write (iout,*) betaT(ib,2,iprot),target_cv(ib,iprot),
1292      &      weiCv(ib,iprot)
1293         enddo
1294         if (me.eq.Master) then
1295         write (iout,*) "Experimental heat capacity curve"
1296         do ib=1,nbeta(2,iprot)
1297           write (iout,'(f6.2,2f10.5)') betaT(ib,2,iprot),
1298      &     target_cv(ib,iprot),weiCv(ib,iprot)
1299         enddo
1300         write (iout,'(a,f10.5)') "Baseline",heatbase(iprot)
1301         call flush(iout)
1302         endif
1303 c Conformation-dependent averages
1304         call card_concat(controlcard,.true.)
1305         call readi(controlcard,"NATLIKE",natlike(iprot),0)
1306         do i=1,natlike(iprot)
1307           call card_concat(controlcard,.true.)
1308           call reada(controlcard,"WNAT",wnat(i,iprot),1.0d0)
1309           call readi(controlcard,"NUMNAT",numnat(i,iprot),1)
1310           call readi(controlcard,"NATDIM",natdim(i,iprot),1)
1311           do ib=1,nbeta(i+2,iprot)
1312             read(inp,*)betaT(ib,i+2,iprot),(weinat(k,ib,i,iprot),
1313      &             nuexp(k,ib,i,iprot),k=1,natdim(i,iprot))
1314           enddo
1315         enddo
1316         do i=1,natlike(iprot)+2
1317           do j=1,nbeta(i,iprot)
1318             betaT(j,i,iprot)=1.0d0/(Rgas*betaT(j,i,iprot))
1319             if (me.eq.Master)
1320      &      write (2,*) "R i",i," j",j," beta",betaT(j,i,iprot)
1321           enddo
1322         enddo
1323         do iparm=1,nparmset
1324
1325 C Read names of files with the data for protein IPROT
1326         call card_concat(controlcard,.false.)
1327
1328         if (iparm.eq.myparm) then
1329         call split_string(controlcard,protfiles(1,iprot),
1330      &     maxfile_prot,nfile_prot(iprot))
1331 #ifdef DEBUG
1332         write(iout,*)"iprot",iprot," nfile",nfile_prot(iprot)
1333         write(iout,*)
1334      &     (protfiles(i,iprot),i=1,nfile_prot(iprot))
1335 #endif
1336         endif
1337
1338         enddo ! iparm
1339
1340 c Read molecular information of the protein
1341         call molread_zs(iprot)
1342 c 3/31/04 AL Read the reference structures of protein iprot
1343 c        print *,"Calling read_ref_structure"
1344         call read_ref_structure(iprot,*11)
1345 #ifdef DEBUG
1346         write (iout,*) "Protein",iprot," left READ_REF_STRUCTURE"
1347 #endif
1348       enddo ! iprot
1349       return
1350    11 return1
1351       end
1352 c-------------------------------------------------------------------------------
1353       subroutine read_database(*)
1354       implicit none
1355       include "DIMENSIONS"
1356       include "DIMENSIONS.ZSCOPT"
1357 #ifdef MPI
1358       include "mpif.h"
1359       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1360       include "COMMON.MPI"
1361 #endif
1362       include "COMMON.CHAIN"
1363       include "COMMON.INTERACT"
1364       include "COMMON.CLASSES"
1365       include "COMMON.ENERGIES"
1366       include "COMMON.IOUNITS"
1367       include "COMMON.PROTFILES"
1368       include "COMMON.PROTNAME"
1369       include "COMMON.VMCPAR"
1370       include "COMMON.WEIGHTS"
1371       include "COMMON.NAMES"
1372       include "COMMON.ALLPROT"
1373       include "COMMON.WEIGHTDER"
1374       include "COMMON.VAR"
1375       include "COMMON.SBRIDGE"
1376       include "COMMON.GEO"
1377       include "COMMON.COMPAR"
1378       include "COMMON.OPTIM"
1379       character*128 nazwa
1380       character*16000 controlcard
1381       character*16000 all_protfiles
1382       character*3 liczba,liczba1
1383       integer i,j,jj,jjj,jj_old,icount,k,kk,l,iprot,ii,if,ib,ibatch,
1384      &  nn,nn1,inan
1385       integer ixdrf,iret,itmp
1386       integer nrec,nlines,iscor
1387       double precision energ,t_acq,tcpu
1388       integer ilen,iroof
1389       external ilen,iroof
1390       double precision rjunk
1391       integer ntot_all(0:maxprot,0:maxprocs-1)
1392       logical lerr
1393       double precision energia(0:max_ene),etot
1394       real*4 csingle(3,maxres2),reini,refree,rmsdev,prec
1395       integer Previous,Next
1396 c      print *,"Processor",me," calls read_protein_data"
1397 #ifdef MPI
1398       if (me.eq.master) then
1399         Previous=MPI_PROC_NULL
1400       else
1401         Previous=me-1
1402       endif
1403       if (me.eq.nprocs-1) then
1404         Next=MPI_PROC_NULL
1405       else
1406         Next=me+1
1407       endif
1408 #endif
1409 c Set the scratchfile names
1410       write (liczba,'(bz,i3.3)') me
1411       write (liczba1,'(bz,i3.3)') myparm
1412       do iprot=1,nprot
1413 c 1/27/05 AL Change stored coordinates to single precision and don't store
1414 c         energy components in the binary databases.
1415         lenrec(iprot)=12*(nres_zs(iprot)+nct_zs(iprot)-nnt_zs(iprot)+1)
1416      &   +4*(2*nss_zs(1,iprot)+1)+8*natlike(iprot)*maxdimnat+16
1417 c 4/13/04 AL Add space for similarity measure
1418         lenrec_ene(iprot) = (2*nntyp+3*n_ene+2)*8
1419      &   +8*natlike(iprot)*maxdimnat
1420 #ifdef DEBUG
1421         call flush(iout)
1422         write (iout,*) "Protein i"," lenrec",lenrec(iprot)
1423         write (iout,*) "lenrec_ene",lenrec_ene(iprot)
1424         call flush(iout)
1425 #endif
1426         bprotfiles(iprot) = scratchdir(:ilen(scratchdir))//
1427      &       "/"//protname(iprot)(:ilen(protname(iprot)))//
1428      &       "_par"//liczba1//"_"//liczba//".xbin"
1429         benefiles(iprot) = scratchdir(:ilen(scratchdir))//
1430      &       "/"//protname(iprot)(:ilen(protname(iprot)))//
1431      &       "_par"//liczba1//"_"//liczba//".enebin"
1432 c        write (iout,*) "scratchfile ",
1433 c     &    bprotfiles(iprot)(:ilen(bprotfiles(iprot)))
1434       enddo
1435       call flush(iout)
1436
1437         do i=0,nprot
1438           ntot(i)=0
1439         enddo
1440         do iprot=1,nprot
1441           jj_old=1
1442           call restore_molinfo(iprot)
1443           open (ientout,file=bprotfiles(iprot),status="unknown",
1444      &      form="unformatted",access="direct",recl=lenrec(iprot))
1445           open (istat,file=benefiles(iprot),status="unknown",
1446      &      form="unformatted",access="direct",recl=lenrec_ene(iprot))
1447 c Read conformations from binary DA files (one per batch) and write them to
1448 c a binary DA scratchfile.
1449           jj=0
1450           jjj=0
1451           write (liczba,'(bz,i3.3)') me
1452 #ifdef MPI
1453           IF (ME.EQ.MASTER) THEN
1454 c Only the master reads the database; it'll send it to the other procs
1455 c through a ring.
1456 #endif
1457             t_acq = tcpu()
1458             icount=0
1459             itmp=0
1460             do if=1,nfile_prot(iprot)
1461               nazwa=protfiles(if,iprot)
1462      &          (:ilen(protfiles(if,iprot)))//".cx"
1463 #ifdef DEBUG
1464               write (iout,*) "Opening file ",nazwa(:ilen(nazwa))
1465 #endif
1466 #if (defined(AIX) && !defined(JUBL))
1467               call xdrfopen_(ixdrf,nazwa, "r", iret)
1468 #else
1469               call xdrfopen(ixdrf,nazwa, "r", iret)
1470 #endif
1471               if (iret.eq.0) goto 1111
1472               do i=1,maxstr
1473
1474                 prec=10000.0
1475 #if (defined(AIX) && !defined(JUBL))
1476                 call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
1477                 if (iret.eq.0) goto 102
1478                 call xdrfint_(ixdrf, nss, iret)
1479                 if (iret.eq.0) goto 102
1480                 do j=1,nss
1481                   call xdrfint_(ixdrf, ihpb(j), iret)
1482                   if (iret.eq.0) goto 102
1483                   call xdrfint_(ixdrf, jhpb(j), iret)
1484                   if (iret.eq.0) goto 102
1485                 enddo
1486                 call xdrffloat_(ixdrf,reini,iret)
1487                 if (iret.eq.0) goto 102
1488                 call xdrffloat_(ixdrf,refree,iret)
1489                 if (iret.eq.0) goto 102
1490                 call xdrfint_(ixdrf,natlik,iret)
1491                 if (iret.eq.0) goto 102
1492                 do j=1,natlik
1493                   call xdrfint(ixdrf,natliktemp(j),iret)
1494                   if (iret.eq.0) goto 102
1495                   do k=1,natliktemp(j)
1496                     call xdrffloat(ixdrf,nutemp(k,j),iret)
1497                     if (iret.eq.0) goto 102
1498                   enddo
1499                 enddo
1500 #else
1501 c                write (0,*) "me",me," iprot",iprot," i",i
1502                 call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret)
1503                 if (iret.eq.0) goto 102
1504                 call xdrfint(ixdrf, nss, iret)
1505                 if (iret.eq.0) goto 102
1506                 do k=1,nss
1507                   call xdrfint(ixdrf, ihpb(k), iret)
1508                   if (iret.eq.0) goto 102
1509                   call xdrfint(ixdrf, jhpb(k), iret)
1510                   if (iret.eq.0) goto 102
1511                 enddo
1512                 call xdrffloat(ixdrf,reini,iret)
1513                 if (iret.eq.0) goto 102
1514                 call xdrffloat(ixdrf,refree,iret)
1515                 if (iret.eq.0) goto 102
1516 #ifdef FINAL
1517                 call xdrfint(ixdrf,natlik,iret)
1518                 if (iret.eq.0) goto 102
1519                 do j=1,natlik
1520                   call xdrfint(ixdrf,natliktemp(j),iret)
1521                   if (iret.eq.0) goto 102
1522                   do k=1,natliktemp(j)
1523                     call xdrffloat(ixdrf,nutemp(k,j),iret)
1524                     if (iret.eq.0) goto 102
1525                   enddo
1526                 enddo
1527 #else
1528                 natlik=0
1529                 call xdrffloat(ixdrf,rmsdev,iret)
1530                 if (iret.eq.0) goto 102
1531 c                write (2,*) "rmsdev",rmsdev
1532                 call xdrfint(ixdrf,iscor,iret)
1533                 if (iret.eq.0) goto 102
1534 c                write (2,*) "iscor",iscor
1535 #endif
1536 #endif
1537                 eini(jj+1,iprot)=reini
1538                 entfac(jj+1,iprot)=refree
1539                 do k=1,nres
1540                   do l=1,3
1541                     c(l,k)=csingle(l,k)
1542                   enddo
1543                 enddo
1544                 do k=nnt,nct
1545                   do l=1,3
1546                     c(l,nres+k)=csingle(l,nres+k-nnt+1)
1547                   enddo
1548                 enddo
1549 #ifdef DEBUG
1550                 write (iout,'(5hREAD ,i5,2f15.4)')
1551      &           jj+1,eini(jj+1,iprot),entfac(jj+1,iprot)
1552                 write (iout,*) "natlik",natlik
1553                 do l=1,natlik
1554                   write (iout,*) "natliketemp",natliktemp(l)
1555                   write(iout,*) (nutemp(k,l),k=1,natliktemp(l))
1556                 enddo
1557                 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1558                 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1559                 call flush(iout)
1560 #endif
1561                 call add_new_cconf(jjj,jj,jj_old,icount,iprot,
1562      &              Next)
1563               enddo
1564   102         continue
1565               if (me.eq.Master)
1566      &         write (iout,*) "Protein ",protname(iprot),
1567      &          i-1," conformations read from DA file ",
1568      &          nazwa(:ilen(nazwa))
1569               if (me.eq.Master)
1570      &        write (iout,*) jj," conformations read so far"
1571 #if (defined(AIX) && !defined(JUBL))
1572               call xdrfclose_(ixdrf,nazwa,iret)
1573 #else
1574               call xdrfclose(ixdrf,nazwa,iret)
1575 #endif
1576 c              print *,"file ",nazwa(:ilen(nazwa))," closed"
1577             enddo
1578 #ifdef MPI
1579 #ifdef DEBUG
1580             write (iout,*) "jj_old",jj_old," jj",jj
1581 #endif
1582             call write_and_send_cconf(icount,jj_old,jj,iprot,Next)
1583             if (icount.gt.0) call MPI_Send(0,1,MPI_INTEGER,Next,570,
1584      &             WHAM_COMM,IERROR)
1585             jj_old=jj+1
1586 #else
1587             call write_and_send_cconf(icount,jj_old,jj,iprot,Next)
1588 #endif
1589             t_acq = tcpu() - t_acq
1590             if (me.eq.Master)
1591      &       write (iout,*) "Processor",me," protein",iprot,
1592      &          " batch",ibatch," time for conformation read/send",t_acq
1593 #ifdef MPI
1594           ELSE
1595 c A worker gets the confs from the master and sends them to its neighbor
1596             t_acq = tcpu()
1597             call receive_and_pass_cconf(icount,jj_old,jj,iprot,
1598      &         Previous,Next)
1599             t_acq = tcpu() - t_acq
1600 c            write (iout,*) "Processor",me," protein",iprot,
1601 c     &          " batch",ibatch,
1602 c     &          " time for conformation receive/send",t_acq
1603           ENDIF
1604 #endif
1605           ntot(iprot)=jj
1606           if (me.eq.Master) write (iout,*) "Protein",
1607      &     protname(iprot)(:ilen(protname(iprot))),", ",ntot(iprot),
1608      &   " conformatons read ",ntot(iprot)," conformations written to ",
1609      &    bprotfiles(iprot)(:ilen(bprotfiles(iprot)))
1610           ntot(0) = ntot(0)+ntot(iprot)
1611           close(ientout)
1612           close(istat)
1613
1614         enddo ! iprot
1615         if (me.eq.Master)
1616      &   write(iout,*)"A total of",ntot(0)," conformations read."
1617 #ifdef MPI
1618 c Check if everyone has the same number of conformations
1619       call MPI_Allgather(ntot(0),maxprot+1,MPI_INTEGER,
1620      &  ntot_all(0,0),maxprot+1,MPI_INTEGER,ALL_Comm,IERROR)
1621       lerr=.false.
1622       do i=0,nprocs-1
1623         if (i.ne.me) then
1624           do j=1,nprot
1625             if (ntot(j).ne.ntot_all(j,i)) then
1626               if (me.eq.Master)
1627      &         write (iout,*) "Number of conformations at processor",i,
1628      &         " for protein",j," differs from that at processor",me,
1629      &         ntot(j),ntot_all(j,i)
1630               lerr = .true.
1631             endif
1632           enddo
1633         endif
1634       enddo
1635 #ifdef WEICLUST
1636 c----------- Temporary; reading probs from external file
1637       open (88,file='1LE1_wham_last_2.rms',status='old')
1638       do i=1,ntot(1)
1639         read (88,*) ii,weirms(i,1)
1640       enddo
1641       do i=1,ntot(1)
1642         if (me.eq.Master) write (iout,*) "i",i," weirms",weirms(i,1)
1643       enddo
1644       close(88)
1645       call MPI_Bcast(weirms(1,1), ntot(1), MPI_Double_Precision,
1646      &        Master, ALL_COMM, IERROR)
1647 c----------- end temportary stuff
1648 #endif
1649       if (lerr) then
1650         if (me.eq.Master) then
1651         write (iout,*)
1652         write (iout,*) "Number of conformations read by processors"
1653         write (iout,'(10x,7a10)') (protname(i),i=0,nprot)
1654         write (iout,*)
1655         do i=0,nprocs-1
1656           write (iout,'(8i10)') i,(ntot_all(j,i),j=0,nprot)
1657         enddo
1658         write (iout,*) "Calculation terminated."
1659         call flush(iout)
1660         endif
1661         return1
1662       endif
1663       return
1664 #endif
1665  1111 if (me.eq.Master)
1666      &write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
1667       call flush(iout)
1668 #ifdef MPI
1669       call MPI_Abort(ALL_COMM,IERROR,ERRCODE)
1670 #endif
1671       return1
1672       end
1673 c------------------------------------------------------------------------------
1674       subroutine add_new_cconf(jjj,jj,jj_old,icount,iprot,Next)
1675       implicit none
1676       include "DIMENSIONS"
1677       include "DIMENSIONS.ZSCOPT"
1678       include "COMMON.CHAIN"
1679       include "COMMON.INTERACT"
1680       include "COMMON.LOCAL"
1681       include "COMMON.CLASSES"
1682       include "COMMON.ENERGIES"
1683       include "COMMON.IOUNITS"
1684       include "COMMON.PROTFILES"
1685       include "COMMON.PROTNAME"
1686       include "COMMON.VMCPAR"
1687       include "COMMON.WEIGHTS"
1688       include "COMMON.NAMES"
1689       include "COMMON.ALLPROT"
1690       include "COMMON.WEIGHTDER"
1691       include "COMMON.VAR"
1692       include "COMMON.SBRIDGE"
1693       include "COMMON.GEO"
1694       include "COMMON.COMPAR"
1695 #ifdef ISNAN
1696 c      integer isnan
1697 #endif
1698       integer i,j,jj,jjj,jj_old,icount,k,kk,l,iprot,ii,ib,ibatch,
1699      &  nn,nn1,inan,Next,itj
1700       double precision etot,energia(0:max_ene)
1701       jjj=jjj+1
1702 c      if (entfac(jj+1,iprot).gt.-10.0d0
1703 c     &     .or. entfac(jj+1,iprot).lt.-150.0d0) then
1704 c        write (iout,*) "Entropy factor out of range for conformation",
1705 c     &    jjj,entfac(jj+1,iprot),", conformation skipped."
1706 c        return
1707 c      endif
1708       call int_from_cart1(.false.)
1709       do j=nnt+1,nct
1710         if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
1711           write (iout,*) "nnt",nnt," nct",nct
1712           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j)
1713           write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
1714           write (iout,*) "The Cartesian geometry is:"
1715           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1716           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1717           write (iout,*) "The internal geometry is:"
1718           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1719           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1720           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1721           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1722           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1723           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1724           write (iout,*)
1725      &      "This conformation WILL NOT be added to the database."
1726           return
1727         endif
1728       enddo
1729       do j=nnt,nct
1730         itj=itype(j)
1731         if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
1732           write (iout,*) "nnt",nnt," nct",nct
1733           write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
1734           write (iout,*) "The Cartesian geometry is:"
1735           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1736           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1737           write (iout,*) "The internal geometry is:"
1738           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1739           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1740           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1741           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1742           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1743           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1744           write (iout,*)
1745      &      "This conformation WILL NOT be added to the database."
1746           return
1747         endif
1748       enddo
1749       do j=3,nres
1750         if (theta(j).le.0.0d0) then
1751           write (iout,*)
1752      &      "Zero theta angle(s) in conformation",ii
1753           write (iout,*) "The Cartesian geometry is:"
1754           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1755           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1756           write (iout,*) "The internal geometry is:"
1757           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1758           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1759           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1760           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1761           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1762           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1763           write (iout,*)
1764      &      "This conformation WILL NOT be added to the database."
1765           return
1766         endif
1767         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
1768       enddo
1769       if (.not. init_ene) then
1770       etot=0.0d0
1771       do j=1,n_ene
1772         etot=etot+ww(j)*enetb(icount+1,j,iprot)
1773       enddo
1774       inan=0
1775 c detecting NaNQ
1776 #ifdef ISNAN
1777 #ifdef AIX
1778       if (isnan(etot).ne.0) inan=1
1779 #else
1780       if (isnan(real(etot))) inan=1
1781 #endif
1782 #else
1783 #ifdef WINPGI
1784       idumm=proc_proc(etot,inan)
1785 #else
1786       call proc_proc(etot,inan)
1787 #endif
1788 #endif
1789       else
1790       inan=0
1791       endif
1792       if (inan.eq.1) then
1793         write (iout,*) "NaNs detected in some of the energy",
1794      &    " components for protein",iprot," batch",ibatch,
1795      &    " conformation",jjj
1796         write (iout,*) "etot",etot
1797         write (iout,*) "The Cartesian geometry is:"
1798         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
1799         write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
1800         write (iout,*) "The internal geometry is:"
1801         write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1802         write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1803         write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1804         write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1805         write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1806         write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1807         write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1808         write (iout,*) "The components of the energy are:"
1809         energia(0)=etot
1810         do k=1,n_ene
1811           energia(k)=enetb(jj+1,k,iprot)
1812         enddo
1813         call enerprint(energia(0))
1814         write (iout,*)
1815      &    "This conformation WILL NOT be added to the database."
1816         call flush(iout)
1817       else
1818         jj=jj+1
1819 #ifdef DEBUG
1820         write (iout,'(e15.5,16i5)') entfac(icount+1,iprot),
1821      &        iscore(icount+1,0,iprot),ibatch
1822         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
1823         write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
1824         write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1825         write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1826         write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1827         write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1828         write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1829         write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1830         write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1831         write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
1832         write (iout,'(5e15.5)') (enetb(icount+1,j,iprot),j=1,n_ene)
1833 c        write (iout,'(2e15.5)') (eneps(1,j,icount+1,iprot),
1834 c     &        eneps(2,j,icount+1,iprot),j=1,nntyp)
1835 #endif
1836 c        write (iout,*) "First nu in readrtms"
1837         call flush(iout)
1838         icount=icount+1
1839         list_conf(jj,iprot)=jj
1840         call store_cconf_from_file(jj,icount,iprot)
1841         do j=1,natlike(iprot)
1842           do k=1,natlik
1843             if (k.eq.numnat(j,iprot)) then
1844               do l=1,natdim(j,iprot)
1845                 nu(l,k,j,iprot)=nutemp(l,k)
1846               enddo
1847               exit
1848             endif
1849           enddo
1850         enddo
1851 c        write (iout,*) "End First nu in readrtms"
1852         call flush(iout)
1853         if (icount.eq.maxstr_proc) then
1854 #ifdef DEBUG
1855           write (iout,* ) "jj_old",jj_old," jj",jj
1856           write (iout,*) "Calling write_and_send_cconf"
1857           call flush(iout)
1858 #endif
1859           call write_and_send_cconf(icount,jj_old,jj,iprot,Next)
1860           jj_old=jj+1
1861 #ifdef DEBUG
1862           write (iout,*) "Exited write_and_send_cconf"
1863           call flush(iout)
1864 #endif
1865           icount=0
1866         endif
1867       endif
1868       return
1869       end
1870 c------------------------------------------------------------------------------
1871       subroutine store_cconf_from_file(jj,icount,iprot)
1872       implicit none
1873       include "DIMENSIONS"
1874       include "DIMENSIONS.ZSCOPT"
1875       include "COMMON.CHAIN"
1876       include "COMMON.SBRIDGE"
1877       include "COMMON.INTERACT"
1878       include "COMMON.IOUNITS"
1879       include "COMMON.CLASSES"
1880       include "COMMON.ALLPROT"
1881       include "COMMON.VAR"
1882       integer i,j,jj,icount,ibatch,iprot
1883 c Store the conformation that has been read in
1884       do i=1,2*nres
1885         do j=1,3
1886           c_zs(j,i,icount,iprot)=c(j,i)
1887         enddo
1888       enddo
1889       nss_zs(icount,iprot)=nss
1890       do i=1,nss
1891         ihpb_zs(i,icount,iprot)=ihpb(i)
1892         jhpb_zs(i,icount,iprot)=jhpb(i)
1893       enddo
1894       return
1895       end
1896 c------------------------------------------------------------------------------
1897       subroutine write_and_send_cconf(icount,jj_old,jj,iprot,Next)
1898       implicit none
1899       include "DIMENSIONS"
1900       include "DIMENSIONS.ZSCOPT"
1901 #ifdef MPI
1902       include "mpif.h"
1903       integer IERROR
1904       include "COMMON.MPI"
1905 #endif
1906       include "COMMON.CHAIN"
1907       include "COMMON.SBRIDGE"
1908       include "COMMON.INTERACT"
1909       include "COMMON.IOUNITS"
1910       include "COMMON.CLASSES"
1911       include "COMMON.VAR"
1912       include "COMMON.ALLPROT"
1913       include "COMMON.ENERGIES"
1914       include "COMMON.WEIGHTDER"
1915       include "COMMON.OPTIM"
1916       include "COMMON.COMPAR"
1917       integer icount,jj_old,jj,Next,iprot
1918 c Write the structures to a scratch file
1919 #ifdef MPI
1920 c Master sends the portion of conformations that have been read in to the neighbor
1921 #ifdef DEBUG
1922       write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
1923       call flush(iout)
1924 #endif
1925       call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM,IERROR)
1926 #ifdef DEBUG
1927       write (iout,*) "Processor",me," Next",next," sent icount=",icount
1928       write (iout,*) "Processor",me," jj_old",jj_old," jj",jj
1929       call flush(iout)
1930 #endif
1931       if (icount.gt.0) then
1932       call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER,
1933      &    Next,571,WHAM_COMM,IERROR)
1934       call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
1935      &    Next,572,WHAM_COMM,IERROR)
1936       if (.not. init_ene) call MPI_Send(enetb(jj_old,1,iprot),
1937      &    maxstr*n_ene,
1938      &    MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,IERROR)
1939       call MPI_Send(nu(1,1,jj_old,iprot),
1940      &    maxdimnat*natlike(iprot)*icount,
1941      &    MPI_DOUBLE_PRECISION,
1942      &    Next,577,WHAM_COMM,IERROR)
1943       call MPI_Send(eini(jj_old,iprot),icount,
1944      &    MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR)
1945       call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
1946      &    Next,579,WHAM_COMM,IERROR)
1947       call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2,
1948      &    MPI_REAL,Next,580,WHAM_COMM,IERROR)
1949       if (.not. init_ene) call MPI_Send(eneps(1,1,1,iprot),
1950      &    2*icount*nntyp,
1951      &    MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR)
1952       endif
1953 #endif
1954       write (2,*) "Calling dawrite_ccoords"
1955       call flush(iout)
1956       call dawrite_ccoords(iprot,jj_old,jj,ientout)
1957       call dawrite_ene(iprot,jj_old,jj,istat)
1958       return
1959       end
1960 c------------------------------------------------------------------------------
1961 #ifdef MPI
1962       subroutine receive_and_pass_cconf(icount,jj_old,jj,iprot,Previous,
1963      &  Next)
1964       implicit none
1965       include "DIMENSIONS"
1966       include "DIMENSIONS.ZSCOPT"
1967       include "mpif.h"
1968       integer IERROR,STATUS(MPI_STATUS_SIZE)
1969       include "COMMON.MPI"
1970       include "COMMON.CHAIN"
1971       include "COMMON.SBRIDGE"
1972       include "COMMON.INTERACT"
1973       include "COMMON.IOUNITS"
1974       include "COMMON.CLASSES"
1975       include "COMMON.ALLPROT"
1976       include "COMMON.ENERGIES"
1977       include "COMMON.VAR"
1978       include "COMMON.GEO"
1979       include "COMMON.WEIGHTDER"
1980       include "COMMON.COMPAR"
1981       include "COMMON.OPTIM"
1982       integer i,j,k,icount,jj_old,jj,iprot,Previous,Next
1983       icount=1
1984 #ifdef DEBUG
1985       write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
1986       call flush(iout)
1987 #endif
1988       do while (icount.gt.0)
1989 c      call MPI_Probe(Previous,570,WHAM_COMM,STATUS,IERROR)
1990       call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,WHAM_COMM,
1991      &     STATUS,IERROR)
1992 #ifdef DEBUG
1993       write (iout,*)"Processor",me," previous",previous," icount",icount
1994       call flush(iout)
1995 #endif
1996       call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM,
1997      &     IERROR)
1998 #ifdef DEBUG
1999       write (iout,*) "Processor",me," icount",icount
2000       call flush(iout)
2001 #endif
2002       if (icount.eq.0) return
2003       call MPI_Recv(nss_zs(1,iprot),icount,MPI_INTEGER,
2004      &    Previous,571,WHAM_COMM,STATUS,IERROR)
2005       call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER,
2006      &  Next,571,WHAM_COMM,IERROR)
2007       call MPI_Recv(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2008      &    Previous,572,WHAM_COMM,STATUS,IERROR)
2009       call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2010      &  Next,572,WHAM_COMM,IERROR)
2011       if (.not. init_ene) then
2012         call MPI_Recv(enetb(jj_old,1,iprot),maxstr*n_ene,
2013      &   MPI_DOUBLE_PRECISION,Previous,576,WHAM_COMM,STATUS,IERROR)
2014         call MPI_Send(enetb(jj_old,1,iprot),maxstr*n_ene,
2015      &   MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,STATUS,IERROR)
2016       endif
2017       call MPI_Recv(nu(1,1,jj_old,iprot),
2018      &  maxdimnat*natlike(iprot)*icount,
2019      &  MPI_DOUBLE_PRECISION,
2020      &  Previous,577,WHAM_COMM,STATUS,IERROR)
2021       call MPI_Send(nu(1,1,jj_old,iprot),
2022      &  maxdimnat*natlike(iprot)*icount,
2023      &  MPI_DOUBLE_PRECISION,
2024      &  Next,577,WHAM_COMM,IERROR)
2025       call MPI_Recv(eini(jj_old,iprot),icount,
2026      &  MPI_DOUBLE_PRECISION,Previous,578,WHAM_COMM,STATUS,IERROR)
2027       call MPI_Send(eini(jj_old,iprot),icount,
2028      &  MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR)
2029       call MPI_Recv(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2030      &  Previous,579,WHAM_COMM,STATUS,IERROR)
2031       call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2032      &  Next,579,WHAM_COMM,IERROR)
2033       call MPI_Recv(c_zs(1,1,1,iprot),3*icount*maxres2,
2034      &  MPI_REAL,Previous,580,WHAM_COMM,STATUS,IERROR)
2035       call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2,
2036      &  MPI_REAL,Next,580,WHAM_COMM,IERROR)
2037       if (.not. init_ene) then
2038         call MPI_Recv(eneps(1,1,1,iprot),2*icount*nntyp,
2039      &   MPI_DOUBLE_PRECISION,Previous,582,WHAM_COMM,STATUS,IERROR)
2040         call MPI_Send(eneps(1,1,1,iprot),2*icount*nntyp,
2041      &   MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR)
2042       endif
2043       jj=jj_old+icount-1
2044       do i=jj_old,jj
2045         list_conf(i,iprot)=i
2046       enddo
2047       call dawrite_ccoords(iprot,jj_old,jj,ientout)
2048       call dawrite_ene(iprot,jj_old,jj,istat)
2049       jj_old=jj+1
2050 #ifdef DEBUG
2051       write (iout,*) "Protein",iprot
2052       write (iout,*) "Processor",me," received",icount," conformations"
2053       do i=1,icount
2054         write (iout,'(8f10.4)') ((c_zs(l,k,i,iprot),l=1,3,k=1,nres)
2055         write (iout,'(8f10.4)')((c_zs(l,k,i+nres,iprot),l=1,3,k=nnt,nct)
2056         write (iout,'(16i5)') nss_zs(i,iprot),(ihpb_zs(k,i,iprot),
2057      &    jhpb_zs(k,i,iprot),k=1,nss_zs(i,iprot))
2058         write (iout,'(5e15.5)') (enetb(i,j,iprot),j=1,n_ene)
2059         write (iout,'(2e15.5)') (eneps(1,j,i,iprot),
2060      &        eneps(2,j,i,iprot),j=1,nntyp)
2061         write (iout,'(e15.5,16i5)') entfac(i),iscore(i,0,iprot),
2062      &        kbatch(i,iprot)
2063       enddo
2064 #endif
2065       enddo
2066       return
2067       end
2068 #endif
2069 c------------------------------------------------------------------------------
2070       subroutine read_thermal
2071       implicit none
2072       include "DIMENSIONS"
2073       include "DIMENSIONS.ZSCOPT"
2074 #ifdef MPI
2075       include "mpif.h"
2076       include "COMMON.MPI"
2077 #endif
2078       include "COMMON.CHAIN"
2079       include "COMMON.SBRIDGE"
2080       include "COMMON.INTERACT"
2081       include "COMMON.IOUNITS"
2082       include "COMMON.CLASSES"
2083       include "COMMON.VAR"
2084       include "COMMON.THERMAL"
2085       character*800 controlcard
2086       call card_concat(controlcard,.true.)
2087       call readi(controlcard,"NGRIDT",NGridT,200)
2088       call reada(controlcard,"DELTAT",deltaT,5.0d0)
2089       call reada(controlcard,"T0",GridT0,2.0d2)
2090       if (me.eq.Master) then
2091       write (iout,*) "Parameters of thermal curves"
2092       write (iout,*) "NGRIDT",NGridT," DELTAT",deltaT," T0",GridT0
2093       endif
2094       return
2095       end
2096 c------------------------------------------------------------------------------
2097       subroutine daread_ene(iprot,istart_conf,iend_conf)
2098       implicit none
2099       include "DIMENSIONS"
2100       include "DIMENSIONS.ZSCOPT"
2101 #ifdef MPI
2102       include "mpif.h"
2103       include "COMMON.MPI"
2104 #endif
2105       include "COMMON.CHAIN"
2106       include "COMMON.CLASSES"
2107       include "COMMON.ENERGIES"
2108       include "COMMON.IOUNITS"
2109       include "COMMON.PROTFILES"
2110       include "COMMON.ALLPROT"
2111       include "COMMON.WEIGHTDER"
2112       include "COMMON.COMPAR"
2113       include "COMMON.VMCPAR"
2114       integer iprot,istart_conf,iend_conf
2115       integer i,ii,iii,j,k
2116 #ifdef DEBUG
2117       write (iout,*) "Calling DAREAD_ENE"
2118 #endif
2119 c      write (iout,*) "Reading: n_ene",n_ene
2120 c      call flush(iout)
2121 c      write (iout,*) "DAREAD_ENE",istart_conf,iend_conf
2122 c
2123 c Read conformations off a DA scratchfile.
2124 c
2125       do ii=istart_conf,iend_conf
2126         iii=list_conf(ii,iprot)
2127         i = ii - istart_conf + 1
2128         read(ientin,rec=iii) (enetb(i,j,iprot),j=1,n_ene),
2129      &     (enetb_orig(i,j,iprot),j=1,n_ene),
2130      &     (enetb_oorig(i,j,iprot),j=1,n_ene),
2131      &     eini(ii,iprot),entfac(ii,iprot),
2132      &     (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp),
2133      &     ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2134 #ifdef DEBUG
2135         write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot),
2136      &  entfac(ii,iprot),
2137         write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene)
2138         write (iout,'(18(1pe12.4))')
2139      &  ((nu(k,j,i,iprot)k=1,maxdimnat),j=1,natlike(iprot))
2140         call flush(iout)
2141 #endif
2142       enddo
2143
2144       return
2145       end
2146 c------------------------------------------------------------------------------
2147       subroutine dawrite_ene(iprot,istart_conf,iend_conf,unit_out)
2148       implicit none
2149       include "DIMENSIONS"
2150       include "DIMENSIONS.ZSCOPT"
2151 #ifdef MPI
2152       include "mpif.h"
2153       include "COMMON.MPI"
2154 #endif
2155       include "COMMON.CHAIN"
2156       include "COMMON.CLASSES"
2157       include "COMMON.ENERGIES"
2158       include "COMMON.IOUNITS"
2159       include "COMMON.PROTFILES"
2160       include "COMMON.ALLPROT"
2161       include "COMMON.WEIGHTDER"
2162       include "COMMON.VMCPAR"
2163       include "COMMON.COMPAR"
2164       integer iprot,istart_conf,iend_conf,unit_out
2165       integer i,ii,iii,j,k
2166 c      write (iout,*) "Writing: n_ene",n_ene
2167 c      call flush(iout)
2168 c      write (iout,*) "DAWRITE_ENE",istart_conf,iend_conf
2169 c
2170 c Write conformations to a DA scratchfile.
2171 c
2172       do ii=istart_conf,iend_conf
2173         iii=list_conf(ii,iprot)
2174         i = ii - istart_conf + 1
2175         write(unit_out,rec=iii) (enetb(i,j,iprot),j=1,n_ene),
2176      &     (enetb_orig(i,j,iprot),j=1,n_ene),
2177      &     (enetb_oorig(i,j,iprot),j=1,n_ene),
2178      &     eini(ii,iprot),entfac(ii,iprot),
2179      &     (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp),
2180      &     ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2181 #ifdef DEBUG
2182         write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot),
2183      &  entfac(ii,iprot)
2184         write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene)
2185         write (iout,'(18(1pe12.4))')
2186      &  ((nu(kj,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2187         call flush(iout)
2188 #endif
2189       enddo
2190
2191       return
2192       end
2193 c------------------------------------------------------------------------------
2194       subroutine daread_ccoords(iprot,istart_conf,iend_conf)
2195       implicit none
2196       include "DIMENSIONS"
2197       include "DIMENSIONS.ZSCOPT"
2198 #ifdef MPI
2199       include "mpif.h"
2200       include "COMMON.MPI"
2201 #endif
2202       include "COMMON.CHAIN"
2203       include "COMMON.CLASSES"
2204       include "COMMON.ENERGIES"
2205       include "COMMON.IOUNITS"
2206       include "COMMON.PROTFILES"
2207       include "COMMON.ALLPROT"
2208       include "COMMON.INTERACT"
2209       include "COMMON.VAR"
2210       include "COMMON.SBRIDGE"
2211       include "COMMON.GEO"
2212       include "COMMON.COMPAR"
2213       include "COMMON.VMCPAR"
2214       include "COMMON.WEIGHTDER"
2215       integer iprot,istart_conf,iend_conf
2216       integer i,j,k,ij,ii,iii
2217       integer len
2218       real*4 csingle(3,maxres2)
2219       character*16 form,acc
2220       character*32 nam
2221 c
2222 c Read conformations off a DA scratchfile.
2223 c
2224 #ifdef DEBUG
2225       write (iout,*) "DAREAD_COORDS"
2226       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2227       write (iout,*) "lenrec",lenrec(iprot)
2228       inquire(unit=ientin,name=nam,recl=len,form=form,access=acc)
2229       write (iout,*) "len=",len," form=",form," acc=",acc
2230       write (iout,*) "nam=",nam
2231       call flush(iout)
2232 #endif
2233       do ii=istart_conf,iend_conf
2234         iii=list_conf(ii,iprot)
2235         ij = ii - istart_conf + 1
2236 #ifdef DEBUG
2237         write (iout,*) "Reading binary file, record",iii," ii",ii
2238         call flush(iout)
2239 #endif
2240         read(icbase,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2241      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2242      &    nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot),
2243      &    entfac(ii,iprot),
2244      &    ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot))
2245         do i=1,2*nres
2246           do j=1,3
2247             c(j,i)=csingle(j,i)
2248           enddo
2249         enddo
2250 #ifdef DEBUG
2251         write (iout,*) "iprot",iprot," ii",ii
2252         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2253         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2254         write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot)
2255         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2256         call flush(iout)
2257 #endif
2258         call store_ccoords(iprot,ii-istart_conf+1)
2259       enddo
2260       return
2261       end
2262 c------------------------------------------------------------------------------
2263       subroutine dawrite_ccoords(iprot,istart_conf,iend_conf,unit_out)
2264       implicit none
2265       include "DIMENSIONS"
2266       include "DIMENSIONS.ZSCOPT"
2267 #ifdef MPI
2268       include "mpif.h"
2269       include "COMMON.MPI"
2270 #endif
2271       include "COMMON.CHAIN"
2272       include "COMMON.INTERACT"
2273       include "COMMON.CLASSES"
2274       include "COMMON.ENERGIES"
2275       include "COMMON.IOUNITS"
2276       include "COMMON.PROTFILES"
2277       include "COMMON.ALLPROT"
2278       include "COMMON.VAR"
2279       include "COMMON.SBRIDGE"
2280       include "COMMON.GEO"
2281       include "COMMON.COMPAR"
2282       include "COMMON.WEIGHTDER"
2283       include "COMMON.VMCPAR"
2284       real*4 csingle(3,maxres2)
2285       integer iprot,istart_conf,iend_conf
2286       integer i,j,k,ii,ij,iii,unit_out
2287       integer len
2288       character*16 form,acc
2289       character*32 nam
2290 c
2291 c Write conformations to a DA scratchfile.
2292 c
2293 #ifdef DEBUG
2294       write (iout,*) "DAWRITE_COORDS"
2295       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2296       write (iout,*) "lenrec",lenrec(iprot)
2297       inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
2298       write (iout,*) "len=",len," form=",form," acc=",acc
2299       write (iout,*) "nam=",nam
2300       call flush(iout)
2301 #endif
2302       do ii=istart_conf,iend_conf
2303         iii=list_conf(ii,iprot)
2304         ij = ii - istart_conf + 1
2305         call restore_ccoords(iprot,ii-istart_conf+1)
2306 #ifdef DEBUG
2307         write (iout,*) "Writing binary file, record",iii," ii",ii
2308         call flush(iout)
2309 #endif
2310         do i=1,2*nres
2311           do j=1,3
2312             csingle(j,i)=c(j,i)
2313           enddo
2314         enddo
2315         write(unit_out,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2316      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2317      &    nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot),
2318      &    entfac(ii,iprot),
2319      &    ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot))
2320 #ifdef DEBUG
2321         write (iout,*) "kbatch",kbatch(ii,iprot)," iscore",
2322      &    iscore(ii,0,iprot)
2323         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2324         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2325         write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot)
2326         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2327         call flush(iout)
2328 #endif
2329       enddo
2330       return
2331       end
2332 c------------------------------------------------------------------------------
2333       subroutine store_ccoords(iprot,ii)
2334       implicit none
2335       include "DIMENSIONS"
2336       include "DIMENSIONS.ZSCOPT"
2337       include "COMMON.VAR"
2338       include "COMMON.CHAIN"
2339       include "COMMON.ALLPROT"
2340       include "COMMON.IOUNITS"
2341       include "COMMON.GEO"
2342       include "COMMON.SBRIDGE"
2343       double precision thetnorm
2344       integer iprot,i,j,ii
2345       do i=1,nres_zs(iprot)
2346         do j=1,3
2347           c_zs(j,i,ii,iprot)=c(j,i)
2348         enddo
2349       enddo
2350       do i=nnt_zs(iprot),nct_zs(iprot)
2351         do j=1,3
2352           c_zs(j,i+nres,ii,iprot)=c(j,i+nres)
2353         enddo
2354       enddo
2355 c 5/7/02 AL: store sbridge info
2356       nss_zs(ii,iprot)=nss
2357       do i=1,nss
2358         ihpb_zs(i,ii,iprot)=ihpb(i)
2359         jhpb_zs(i,ii,iprot)=jhpb(i)
2360       enddo
2361       return
2362       end
2363 c------------------------------------------------------------------------------
2364       subroutine restore_ccoords(iprot,ii)
2365       implicit none
2366       include "DIMENSIONS"
2367       include "DIMENSIONS.ZSCOPT"
2368       include "COMMON.INTERACT"
2369       include "COMMON.VAR"
2370       include "COMMON.ALLPROT"
2371       include "COMMON.SBRIDGE"
2372       include "COMMON.CHAIN"
2373       include "COMMON.IOUNITS"
2374       integer iprot,i,j,ii
2375       do i=1,nres_zs(iprot)
2376         do j=1,3
2377           c(j,i)=c_zs(j,i,ii,iprot)
2378         enddo
2379       enddo
2380       do i=nnt_zs(iprot),nct_zs(iprot)
2381         do j=1,3
2382           c(j,i+nres)=c_zs(j,i+nres,ii,iprot)
2383         enddo
2384       enddo
2385 c 5/7/02 AL: restore sbridge info
2386       nss=nss_zs(ii,iprot)
2387       do i=1,nss
2388         ihpb(i)=ihpb_zs(i,ii,iprot)+nres
2389         jhpb(i)=jhpb_zs(i,ii,iprot)+nres
2390 c        forcon(i)=fbr
2391 c        dhpb(i)=dbr
2392       enddo
2393 #ifdef DEBUG
2394         write (iout,*) "restore_ccoords",ii
2395         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2396         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2397         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2398         call flush(iout)
2399 #endif
2400       return
2401       end
2402 c------------------------------------------------------------------------------
2403       subroutine card_concat(card,to_upper)
2404       implicit none
2405       include 'DIMENSIONS.ZSCOPT'
2406       include "COMMON.IOUNITS"
2407       character*(*) card
2408       character*80 karta,ucase
2409       logical to_upper
2410       integer ilen
2411       external ilen
2412       read (inp,'(a)') karta
2413       if (to_upper) karta=ucase(karta)
2414       card=' '
2415       do while (karta(80:80).eq.'&')
2416         card=card(:ilen(card)+1)//karta(:79)
2417         read (inp,'(a)') karta
2418         if (to_upper) karta=ucase(karta)
2419       enddo
2420       card=card(:ilen(card)+1)//karta
2421       return
2422       end
2423 c------------------------------------------------------------------------------
2424       subroutine readi(rekord,lancuch,wartosc,default)
2425       implicit none
2426       character*(*) rekord,lancuch
2427       integer wartosc,default
2428       integer ilen,iread
2429       external ilen
2430       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2431       if (iread.eq.0) then
2432         wartosc=default
2433         return
2434       endif
2435       iread=iread+ilen(lancuch)+1
2436       read (rekord(iread:),*) wartosc
2437       return
2438       end
2439 c----------------------------------------------------------------------------
2440       subroutine reada(rekord,lancuch,wartosc,default)
2441       implicit none
2442       character*(*) rekord,lancuch
2443       character*80 aux
2444       double precision wartosc,default
2445       integer ilen,iread
2446       external ilen
2447       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2448       if (iread.eq.0) then
2449         wartosc=default
2450         return
2451       endif
2452       iread=iread+ilen(lancuch)+1
2453       read (rekord(iread:),*) wartosc
2454       return
2455       end
2456 c----------------------------------------------------------------------------
2457       subroutine multreadi(rekord,lancuch,tablica,dim,default)
2458       implicit none
2459       integer dim,i
2460       integer tablica(dim),default
2461       character*(*) rekord,lancuch
2462       character*80 aux
2463       integer ilen,iread
2464       external ilen
2465       do i=1,dim
2466         tablica(i)=default
2467       enddo
2468       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2469       if (iread.eq.0) return
2470       iread=iread+ilen(lancuch)+1
2471       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2472    10 return
2473       end
2474 c----------------------------------------------------------------------------
2475       subroutine multreada(rekord,lancuch,tablica,dim,default)
2476       implicit none
2477       integer dim,i
2478       double precision tablica(dim),default
2479       character*(*) rekord,lancuch
2480       character*80 aux
2481       integer ilen,iread
2482       external ilen
2483       do i=1,dim
2484         tablica(i)=default
2485       enddo
2486       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2487       if (iread.eq.0) return
2488       iread=iread+ilen(lancuch)+1
2489       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2490    10 return
2491       end
2492 c----------------------------------------------------------------------------
2493       subroutine reads(rekord,lancuch,wartosc,default)
2494       implicit none
2495       character*(*) rekord,lancuch,wartosc,default
2496       character*80 aux
2497       integer ilen,lenlan,lenrec,iread,ireade
2498       external ilen
2499       logical iblnk
2500       external iblnk
2501       lenlan=ilen(lancuch)
2502       lenrec=ilen(rekord)
2503       iread=index(rekord,lancuch(:lenlan)//"=")
2504 c      print *,"rekord",rekord," lancuch",lancuch
2505 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2506       if (iread.eq.0) then
2507         wartosc=default
2508         return
2509       endif
2510       iread=iread+lenlan+1
2511 c      print *,"iread",iread
2512 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2513       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2514         iread=iread+1
2515 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2516       enddo
2517 c      print *,"iread",iread
2518       if (iread.gt.lenrec) then
2519          wartosc=default
2520         return
2521       endif
2522       ireade=iread+1
2523 c      print *,"ireade",ireade
2524       do while (ireade.lt.lenrec .and.
2525      &   .not.iblnk(rekord(ireade:ireade)))
2526         ireade=ireade+1
2527       enddo
2528       wartosc=rekord(iread:ireade)
2529       return
2530       end
2531 c----------------------------------------------------------------------------
2532       subroutine multreads(rekord,lancuch,tablica,dim,default)
2533       implicit none
2534       integer dim,i
2535       character*(*) rekord,lancuch,tablica(dim),default
2536       character*80 aux
2537       integer ilen,lenlan,lenrec,iread,ireade
2538       external ilen
2539       logical iblnk
2540       external iblnk
2541       do i=1,dim
2542         tablica(i)=default
2543       enddo
2544       lenlan=ilen(lancuch)
2545       lenrec=ilen(rekord)
2546       iread=index(rekord,lancuch(:lenlan)//"=")
2547 c      print *,"rekord",rekord," lancuch",lancuch
2548 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2549       if (iread.eq.0) return
2550       iread=iread+lenlan+1
2551       do i=1,dim
2552 c      print *,"iread",iread
2553 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2554       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2555         iread=iread+1
2556 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2557       enddo
2558 c      print *,"iread",iread
2559       if (iread.gt.lenrec) return
2560       ireade=iread+1
2561 c      print *,"ireade",ireade
2562       do while (ireade.lt.lenrec .and.
2563      &   .not.iblnk(rekord(ireade:ireade)))
2564         ireade=ireade+1
2565       enddo
2566       tablica(i)=rekord(iread:ireade)
2567       iread=ireade+1
2568       enddo
2569       end
2570 c----------------------------------------------------------------------------
2571       subroutine split_string(rekord,tablica,dim,nsub)
2572       implicit none
2573       integer dim,nsub,i,ii,ll,kk
2574       character*(*) tablica(dim)
2575       character*(*) rekord
2576       integer ilen
2577       external ilen
2578       do i=1,dim
2579         tablica(i)=" "
2580       enddo
2581       ii=1
2582       ll = ilen(rekord)
2583       nsub=0
2584       do i=1,dim
2585 C Find the start of term name
2586         kk = 0
2587         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
2588           ii = ii+1
2589         enddo
2590 C Parse the name into TABLICA(i) until blank found
2591         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
2592           kk = kk+1
2593           tablica(i)(kk:kk)=rekord(ii:ii)
2594           ii = ii+1
2595         enddo
2596         if (kk.gt.0) nsub=nsub+1
2597         if (ii.gt.ll) return
2598       enddo
2599       return
2600       end
2601 c-------------------------------------------------------------------------------
2602       integer function iroof(n,m)
2603       ii = n/m
2604       if (ii*m .lt. n) ii=ii+1
2605       iroof = ii
2606       return
2607       end