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