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