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