source/maxlik/src_MD_T_maxlik-NEWCORR-PMF/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)
 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)
 407               ityp_psc(2,npair_sc)=rescode(1,item(2),0)
 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)
 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))
 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))
 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))
 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))
 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)
 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)
 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)
 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)
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 (vbld(j).lt.2.0d0 .or. vbld(j).gt.6.5d0) then
1988           write (iout,*) "nnt",nnt," nct",nct
1989           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j)
1990           write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
1991           write (iout,*) "The Cartesian geometry is:"
1992           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1993           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1994           write (iout,*) "The internal geometry is:"
1995           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1996           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1997           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1998           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1999           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2000           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2001           write (iout,*)
2002      &      "This conformation WILL NOT be added to the database."
2003           return
2004         endif
2005       enddo
2006       do j=nnt,nct
2007         itj=itype(j)
2008         if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
2009           write (iout,*) "nnt",nnt," nct",nct
2010           write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
2011           write (iout,*) "The Cartesian geometry is:"
2012           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
2013           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
2014           write (iout,*) "The internal geometry is:"
2015           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2016           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2017           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2018           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2019           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2020           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2021           write (iout,*)
2022      &      "This conformation WILL NOT be added to the database."
2023           return
2024         endif
2025       enddo
2026       do j=3,nres
2027         if (theta(j).le.0.0d0) then
2028           write (iout,*)
2029      &      "Zero theta angle(s) in conformation",ii
2030           write (iout,*) "The Cartesian geometry is:"
2031           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
2032           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
2033           write (iout,*) "The internal geometry is:"
2034           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2035           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2036           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2037           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2038           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2039           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2040           write (iout,*)
2041      &      "This conformation WILL NOT be added to the database."
2042           return
2043         endif
2044         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
2045       enddo
2046       if (.not. init_ene) then
2047       etot=0.0d0
2048       do j=1,n_ene
2049         etot=etot+ww(j)*enetb(icount+1,j,iprot)
2050       enddo
2051       inan=0
2052 c detecting NaNQ
2053 #ifdef ISNAN
2054 #ifdef AIX
2055       if (isnan(etot).ne.0) inan=1
2056 #else
2057       if (isnan(etot)) inan=1
2058 #endif
2059 #else
2060 #ifdef WINPGI
2061       idumm=proc_proc(etot,inan)
2062 #else
2063       call proc_proc(etot,inan)
2064 #endif
2065 #endif
2066       else
2067       inan=0
2068       endif
2069       if (inan.eq.1) then
2070         write (iout,*) "NaNs detected in some of the energy",
2071      &    " components for protein",iprot," batch",ibatch,
2072      &    " conformation",jjj
2073         write (iout,*) "etot",etot
2074         write (iout,*) "The Cartesian geometry is:"
2075         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2076         write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
2077         write (iout,*) "The internal geometry is:"
2078         write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2079         write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2080         write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2081         write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2082         write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
2083         write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2084         write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2085         write (iout,*) "The components of the energy are:"
2086         energia(0)=etot
2087         do k=1,n_ene
2088           energia(k)=enetb(jj+1,k,iprot)
2089         enddo
2090         call enerprint(energia(0))
2091         write (iout,*)
2092      &    "This conformation WILL NOT be added to the database."
2093         call flush(iout)
2094       else
2095         jj=jj+1
2096 #ifdef DEBUG
2097         write (iout,'(e15.5,16i5)') entfac(icount+1,iprot),
2098      &        iscore(icount+1,0,iprot),ibatch
2099         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2100         write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
2101         write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2102         write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2103         write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2104         write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2105         write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
2106         write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2107         write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2108         write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
2109         write (iout,'(5e15.5)') (enetb(icount+1,j,iprot),j=1,n_ene)
2110 c        write (iout,'(2e15.5)') (eneps(1,j,icount+1,iprot),
2111 c     &        eneps(2,j,icount+1,iprot),j=1,nntyp)
2112 #endif
2113 c        write (iout,*) "First nu in readrtms"
2114         call flush(iout)
2115         icount=icount+1
2116         list_conf(jj,iprot)=jj
2117         call store_cconf_from_file(jj,icount,iprot)
2118         do j=1,natlike(iprot)
2119           do k=1,natlik
2120             if (k.eq.numnat(j,iprot)) then
2121               do l=1,natdim(j,iprot)
2122                 nu(l,k,j,iprot)=nutemp(l,k)
2123               enddo
2124               exit
2125             endif
2126           enddo
2127         enddo
2128 c        write (iout,*) "End First nu in readrtms"
2129         call flush(iout)
2130         if (icount.eq.maxstr_proc) then
2131 #ifdef DEBUG
2132           write (iout,* ) "jj_old",jj_old," jj",jj
2133           write (iout,*) "Calling write_and_send_cconf"
2134           call flush(iout)
2135 #endif
2136           call write_and_send_cconf(icount,jj_old,jj,iprot,Next)
2137           jj_old=jj+1
2138 #ifdef DEBUG
2139           write (iout,*) "Exited write_and_send_cconf"
2140           call flush(iout)
2141 #endif
2142           icount=0
2143         endif
2144       endif
2145       return
2146       end
2147 c------------------------------------------------------------------------------
2148       subroutine store_cconf_from_file(jj,icount,iprot)
2149       implicit none
2150       include "DIMENSIONS"
2151       include "DIMENSIONS.ZSCOPT"
2152       include "COMMON.CHAIN"
2153       include "COMMON.SBRIDGE"
2154       include "COMMON.INTERACT"
2155       include "COMMON.IOUNITS"
2156       include "COMMON.CLASSES"
2157       include "COMMON.ALLPROT"
2158       include "COMMON.VAR"
2159       integer i,j,jj,icount,ibatch,iprot
2160 c Store the conformation that has been read in
2161       do i=1,2*nres
2162         do j=1,3
2163           c_zs(j,i,icount,iprot)=c(j,i)
2164         enddo
2165       enddo
2166       nss_zs(icount,iprot)=nss
2167       do i=1,nss
2168         ihpb_zs(i,icount,iprot)=ihpb(i)
2169         jhpb_zs(i,icount,iprot)=jhpb(i)
2170       enddo
2171       return
2172       end
2173 c------------------------------------------------------------------------------
2174       subroutine write_and_send_cconf(icount,jj_old,jj,iprot,Next)
2175       implicit none
2176       include "DIMENSIONS"
2177       include "DIMENSIONS.ZSCOPT"
2178 #ifdef MPI
2179       include "mpif.h"
2180       integer IERROR
2181       include "COMMON.MPI"
2182 #endif
2183       include "COMMON.WEIGHTS"
2184       include "COMMON.CHAIN"
2185       include "COMMON.SBRIDGE"
2186       include "COMMON.INTERACT"
2187       include "COMMON.IOUNITS"
2188       include "COMMON.CLASSES"
2189       include "COMMON.VAR"
2190       include "COMMON.ALLPROT"
2191       include "COMMON.ENERGIES"
2192       include "COMMON.WEIGHTDER"
2193       include "COMMON.OPTIM"
2194       include "COMMON.COMPAR"
2195       integer icount,jj_old,jj,Next,iprot
2196 c Write the structures to a scratch file
2197 #ifdef MPI
2198 c Master sends the portion of conformations that have been read in to the neighbor
2199 #ifdef DEBUG
2200       write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
2201       call flush(iout)
2202 #endif
2203       call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM,IERROR)
2204 #ifdef DEBUG
2205       write (iout,*) "Processor",me," Next",next," sent icount=",icount
2206       write (iout,*) "Processor",me," jj_old",jj_old," jj",jj
2207       call flush(iout)
2208 #endif
2209       if (icount.gt.0) then
2210       call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER,
2211      &    Next,571,WHAM_COMM,IERROR)
2212       call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2213      &    Next,572,WHAM_COMM,IERROR)
2214       if (.not. init_ene) call MPI_Send(enetb(jj_old,1,iprot),
2215      &    maxstr*n_ene,
2216      &    MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,IERROR)
2217       call MPI_Send(nu(1,1,jj_old,iprot),
2218      &    maxdimnat*natlike(iprot)*icount,
2219      &    MPI_DOUBLE_PRECISION,
2220      &    Next,577,WHAM_COMM,IERROR)
2221       call MPI_Send(eini(jj_old,iprot),icount,
2222      &    MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR)
2223       call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2224      &    Next,579,WHAM_COMM,IERROR)
2225       call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2,
2226      &    MPI_REAL,Next,580,WHAM_COMM,IERROR)
2227       if (.not. init_ene) call MPI_Send(eneps(1,1,1,iprot),
2228      &    2*icount*nntyp,
2229      &    MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR)
2230       endif
2231 #endif
2232       call flush(iout)
2233       call dawrite_ccoords(iprot,jj_old,jj,ientout)
2234 c Change AL 20/12/2017
2235       if (.not. mod_other_params)
2236      &call dawrite_ene(iprot,jj_old,jj,istat)
2237       return
2238       end
2239 c------------------------------------------------------------------------------
2240 #ifdef MPI
2241       subroutine receive_and_pass_cconf(icount,jj_old,jj,iprot,Previous,
2242      &  Next)
2243       implicit none
2244       include "DIMENSIONS"
2245       include "DIMENSIONS.ZSCOPT"
2246       include "mpif.h"
2247       integer IERROR,STATUS(MPI_STATUS_SIZE)
2248       include "COMMON.MPI"
2249       include "COMMON.CHAIN"
2250       include "COMMON.SBRIDGE"
2251       include "COMMON.INTERACT"
2252       include "COMMON.IOUNITS"
2253       include "COMMON.CLASSES"
2254       include "COMMON.ALLPROT"
2255       include "COMMON.ENERGIES"
2256       include "COMMON.VAR"
2257       include "COMMON.GEO"
2258       include "COMMON.WEIGHTS"
2259       include "COMMON.WEIGHTDER"
2260       include "COMMON.COMPAR"
2261       include "COMMON.OPTIM"
2262       integer i,j,k,icount,jj_old,jj,iprot,Previous,Next
2263       icount=1
2264 #ifdef DEBUG
2265       write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
2266       call flush(iout)
2267 #endif
2268       do while (icount.gt.0)
2269 c      call MPI_Probe(Previous,570,WHAM_COMM,STATUS,IERROR)
2270       call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,WHAM_COMM,
2271      &     STATUS,IERROR)
2272 #ifdef DEBUG
2273       write (iout,*)"Processor",me," previous",previous," icount",icount
2274       call flush(iout)
2275 #endif
2276       call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM,
2277      &     IERROR)
2278 #ifdef DEBUG
2279       write (iout,*) "Processor",me," icount",icount
2280       call flush(iout)
2281 #endif
2282       if (icount.eq.0) return
2283       call MPI_Recv(nss_zs(1,iprot),icount,MPI_INTEGER,
2284      &    Previous,571,WHAM_COMM,STATUS,IERROR)
2285       call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER,
2286      &  Next,571,WHAM_COMM,IERROR)
2287       call MPI_Recv(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2288      &    Previous,572,WHAM_COMM,STATUS,IERROR)
2289       call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2290      &  Next,572,WHAM_COMM,IERROR)
2291       if (.not. init_ene) then
2292         call MPI_Recv(enetb(jj_old,1,iprot),maxstr*n_ene,
2293      &   MPI_DOUBLE_PRECISION,Previous,576,WHAM_COMM,STATUS,IERROR)
2294         call MPI_Send(enetb(jj_old,1,iprot),maxstr*n_ene,
2295      &   MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,STATUS,IERROR)
2296       endif
2297       call MPI_Recv(nu(1,1,jj_old,iprot),
2298      &  maxdimnat*natlike(iprot)*icount,
2299      &  MPI_DOUBLE_PRECISION,
2300      &  Previous,577,WHAM_COMM,STATUS,IERROR)
2301       call MPI_Send(nu(1,1,jj_old,iprot),
2302      &  maxdimnat*natlike(iprot)*icount,
2303      &  MPI_DOUBLE_PRECISION,
2304      &  Next,577,WHAM_COMM,IERROR)
2305       call MPI_Recv(eini(jj_old,iprot),icount,
2306      &  MPI_DOUBLE_PRECISION,Previous,578,WHAM_COMM,STATUS,IERROR)
2307       call MPI_Send(eini(jj_old,iprot),icount,
2308      &  MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR)
2309       call MPI_Recv(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2310      &  Previous,579,WHAM_COMM,STATUS,IERROR)
2311       call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2312      &  Next,579,WHAM_COMM,IERROR)
2313       call MPI_Recv(c_zs(1,1,1,iprot),3*icount*maxres2,
2314      &  MPI_REAL,Previous,580,WHAM_COMM,STATUS,IERROR)
2315       call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2,
2316      &  MPI_REAL,Next,580,WHAM_COMM,IERROR)
2317       if (.not. init_ene) then
2318         call MPI_Recv(eneps(1,1,1,iprot),2*icount*nntyp,
2319      &   MPI_DOUBLE_PRECISION,Previous,582,WHAM_COMM,STATUS,IERROR)
2320         call MPI_Send(eneps(1,1,1,iprot),2*icount*nntyp,
2321      &   MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR)
2322       endif
2323       jj=jj_old+icount-1
2324       do i=jj_old,jj
2325         list_conf(i,iprot)=i
2326       enddo
2327       call dawrite_ccoords(iprot,jj_old,jj,ientout)
2328 c Change AL 12/20/2017
2329       if (.not. mod_other_params)
2330      &call dawrite_ene(iprot,jj_old,jj,istat)
2331       jj_old=jj+1
2332 #ifdef DEBUG
2333       write (iout,*) "Protein",iprot
2334       write (iout,*) "Processor",me," received",icount," conformations"
2335       do i=1,icount
2336         write (iout,'(8f10.4)') ((c_zs(l,k,i,iprot),l=1,3,k=1,nres)
2337         write (iout,'(8f10.4)')((c_zs(l,k,i+nres,iprot),l=1,3,k=nnt,nct)
2338         write (iout,'(16i5)') nss_zs(i,iprot),(ihpb_zs(k,i,iprot),
2339      &    jhpb_zs(k,i,iprot),k=1,nss_zs(i,iprot))
2340         write (iout,'(5e15.5)') (enetb(i,j,iprot),j=1,n_ene)
2341         write (iout,'(2e15.5)') (eneps(1,j,i,iprot),
2342      &        eneps(2,j,i,iprot),j=1,nntyp)
2343         write (iout,'(e15.5,16i5)') entfac(i),iscore(i,0,iprot),
2344      &        kbatch(i,iprot)
2345       enddo
2346 #endif
2347       enddo
2348       return
2349       end
2350 #endif
2351 c------------------------------------------------------------------------------
2352       subroutine read_thermal
2353       implicit none
2354       include "DIMENSIONS"
2355       include "DIMENSIONS.ZSCOPT"
2356       include "COMMON.CHAIN"
2357       include "COMMON.SBRIDGE"
2358       include "COMMON.INTERACT"
2359       include "COMMON.IOUNITS"
2360       include "COMMON.CLASSES"
2361       include "COMMON.VAR"
2362       include "COMMON.THERMAL"
2363       character*800 controlcard
2364       call card_concat(controlcard,.true.)
2365       call readi(controlcard,"NGRIDT",NGridT,200)
2366       call reada(controlcard,"DELTAT",deltaT,5.0d0)
2367       call reada(controlcard,"T0",GridT0,2.0d2)
2368       write (iout,*) "Parameters of thermal curves"
2369       write (iout,*) "NGRIDT",NGridT," DELTAT",deltaT," T0",GridT0
2370       return
2371       end
2372 c------------------------------------------------------------------------------
2373       subroutine daread_ene(iprot,istart_conf,iend_conf)
2374       implicit none
2375       include "DIMENSIONS"
2376       include "DIMENSIONS.ZSCOPT"
2377 #ifdef MPI
2378       include "mpif.h"
2379       include "COMMON.MPI"
2380 #endif
2381       include "COMMON.CHAIN"
2382       include "COMMON.CLASSES"
2383       include "COMMON.ENERGIES"
2384       include "COMMON.IOUNITS"
2385       include "COMMON.PROTFILES"
2386       include "COMMON.ALLPROT"
2387       include "COMMON.WEIGHTDER"
2388       include "COMMON.COMPAR"
2389       include "COMMON.VMCPAR"
2390       integer iprot,istart_conf,iend_conf
2391       integer i,ii,iii,j,k
2392 #ifdef DEBUG
2393       write (iout,*) "Calling DAREAD_ENE"
2394 #endif
2395 c      write (iout,*) "Reading: n_ene",n_ene
2396 c      call flush(iout)
2397 c      write (iout,*) "DAREAD_ENE",istart_conf,iend_conf
2398 c
2399 c Read conformations off a DA scratchfile.
2400 c
2401       do ii=istart_conf,iend_conf
2402         iii=list_conf(ii,iprot)
2403         i = ii - istart_conf + 1
2404         read(ientin,rec=iii) (enetb(i,j,iprot),j=1,n_ene),
2405      &     (enetb_orig(i,j,iprot),j=1,n_ene),
2406      &     (enetb_oorig(i,j,iprot),j=1,n_ene),
2407      &     eini(ii,iprot),entfac(ii,iprot),
2408      &     (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp),
2409      &     ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2410 #ifdef DEBUG
2411         write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot),
2412      &  entfac(ii,iprot),
2413         write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene)
2414         write (iout,'(18(1pe12.4))')
2415      &  ((nu(k,j,i,iprot)k=1,maxdimnat),j=1,natlike(iprot))
2416         call flush(iout)
2417 #endif
2418       enddo
2419
2420       return
2421       end
2422 c------------------------------------------------------------------------------
2423       subroutine dawrite_ene(iprot,istart_conf,iend_conf,unit_out)
2424       implicit none
2425       include "DIMENSIONS"
2426       include "DIMENSIONS.ZSCOPT"
2427 #ifdef MPI
2428       include "mpif.h"
2429       include "COMMON.MPI"
2430 #endif
2431       include "COMMON.CHAIN"
2432       include "COMMON.CLASSES"
2433       include "COMMON.ENERGIES"
2434       include "COMMON.IOUNITS"
2435       include "COMMON.PROTFILES"
2436       include "COMMON.ALLPROT"
2437       include "COMMON.WEIGHTDER"
2438       include "COMMON.VMCPAR"
2439       include "COMMON.COMPAR"
2440       integer iprot,istart_conf,iend_conf,unit_out
2441       integer i,ii,iii,j,k
2442 c      write (iout,*) "Writing: n_ene",n_ene
2443 c      call flush(iout)
2444 c      write (iout,*) "DAWRITE_ENE",istart_conf,iend_conf
2445 c
2446 c Write conformations to a DA scratchfile.
2447 c
2448       do ii=istart_conf,iend_conf
2449         iii=list_conf(ii,iprot)
2450         i = ii - istart_conf + 1
2451         write(unit_out,rec=iii) (enetb(i,j,iprot),j=1,n_ene),
2452      &     (enetb_orig(i,j,iprot),j=1,n_ene),
2453      &     (enetb_oorig(i,j,iprot),j=1,n_ene),
2454      &     eini(ii,iprot),entfac(ii,iprot),
2455      &     (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp),
2456      &     ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2457 #ifdef DEBUG
2458         write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot),
2459      &  entfac(ii,iprot)
2460         write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene)
2461         write (iout,'(18(1pe12.4))')
2462      &  ((nu(kj,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2463         call flush(iout)
2464 #endif
2465       enddo
2466
2467       return
2468       end
2469 c------------------------------------------------------------------------------
2470       subroutine daread_ccoords(iprot,istart_conf,iend_conf)
2471       implicit none
2472       include "DIMENSIONS"
2473       include "DIMENSIONS.ZSCOPT"
2474 #ifdef MPI
2475       include "mpif.h"
2476       include "COMMON.MPI"
2477 #endif
2478       include "COMMON.CHAIN"
2479       include "COMMON.CLASSES"
2480       include "COMMON.ENERGIES"
2481       include "COMMON.IOUNITS"
2482       include "COMMON.PROTFILES"
2483       include "COMMON.ALLPROT"
2484       include "COMMON.INTERACT"
2485       include "COMMON.VAR"
2486       include "COMMON.SBRIDGE"
2487       include "COMMON.GEO"
2488       include "COMMON.COMPAR"
2489       include "COMMON.VMCPAR"
2490       include "COMMON.WEIGHTDER"
2491       integer iprot,istart_conf,iend_conf
2492       integer i,j,k,ij,ii,iii
2493       integer len
2494       real*4 csingle(3,maxres2)
2495       character*16 form,acc
2496       character*32 nam
2497 c
2498 c Read conformations off a DA scratchfile.
2499 c
2500 #ifdef DEBUG
2501       write (iout,*) "DAREAD_COORDS"
2502       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2503       write (iout,*) "lenrec",lenrec(iprot)
2504       inquire(unit=ientin,name=nam,recl=len,form=form,access=acc)
2505       write (iout,*) "len=",len," form=",form," acc=",acc
2506       write (iout,*) "nam=",nam
2507       call flush(iout)
2508 #endif
2509       do ii=istart_conf,iend_conf
2510         iii=list_conf(ii,iprot)
2511         ij = ii - istart_conf + 1
2512 #ifdef DEBUG
2513         write (iout,*) "Reading binary file, record",iii," ii",ii
2514         call flush(iout)
2515 #endif
2516         read(icbase,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2517      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2518      &    nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot),
2519      &    entfac(ii,iprot),
2520      &    ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot))
2521         do i=1,2*nres
2522           do j=1,3
2523             c(j,i)=csingle(j,i)
2524           enddo
2525         enddo
2526 #ifdef DEBUG
2527         write (iout,*) "iprot",iprot," ii",ii
2528         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2529         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2530         write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot)
2531         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2532         call flush(iout)
2533 #endif
2534         call store_ccoords(iprot,ii-istart_conf+1)
2535       enddo
2536       return
2537       end
2538 c------------------------------------------------------------------------------
2539       subroutine dawrite_ccoords(iprot,istart_conf,iend_conf,unit_out)
2540       implicit none
2541       include "DIMENSIONS"
2542       include "DIMENSIONS.ZSCOPT"
2543 #ifdef MPI
2544       include "mpif.h"
2545       include "COMMON.MPI"
2546 #endif
2547       include "COMMON.CHAIN"
2548       include "COMMON.INTERACT"
2549       include "COMMON.CLASSES"
2550       include "COMMON.ENERGIES"
2551       include "COMMON.IOUNITS"
2552       include "COMMON.PROTFILES"
2553       include "COMMON.ALLPROT"
2554       include "COMMON.VAR"
2555       include "COMMON.SBRIDGE"
2556       include "COMMON.GEO"
2557       include "COMMON.COMPAR"
2558       include "COMMON.WEIGHTDER"
2559       include "COMMON.VMCPAR"
2560       real*4 csingle(3,maxres2)
2561       integer iprot,istart_conf,iend_conf
2562       integer i,j,k,ii,ij,iii,unit_out
2563       integer len
2564       character*16 form,acc
2565       character*32 nam
2566 c
2567 c Write conformations to a DA scratchfile.
2568 c
2569 #ifdef DEBUG
2570       write (iout,*) "DAWRITE_COORDS"
2571       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2572       write (iout,*) "lenrec",lenrec(iprot)
2573       inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
2574       write (iout,*) "len=",len," form=",form," acc=",acc
2575       write (iout,*) "nam=",nam
2576       call flush(iout)
2577 #endif
2578       do ii=istart_conf,iend_conf
2579         iii=list_conf(ii,iprot)
2580         ij = ii - istart_conf + 1
2581         call restore_ccoords(iprot,ii-istart_conf+1)
2582 #ifdef DEBUG
2583         write (iout,*) "Writing binary file, record",iii," ii",ii
2584         call flush(iout)
2585 #endif
2586         do i=1,2*nres
2587           do j=1,3
2588             csingle(j,i)=c(j,i)
2589           enddo
2590         enddo
2591         write(unit_out,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2592      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2593      &    nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot),
2594      &    entfac(ii,iprot),
2595      &    ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot))
2596 #ifdef DEBUG
2597         write (iout,*) "kbatch",kbatch(ii,iprot)," iscore",
2598      &    iscore(ii,0,iprot)
2599         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2600         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2601         write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot)
2602         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2603         call flush(iout)
2604 #endif
2605       enddo
2606       return
2607       end
2608 c------------------------------------------------------------------------------
2609       subroutine store_ccoords(iprot,ii)
2610       implicit none
2611       include "DIMENSIONS"
2612       include "DIMENSIONS.ZSCOPT"
2613       include "COMMON.VAR"
2614       include "COMMON.CHAIN"
2615       include "COMMON.ALLPROT"
2616       include "COMMON.IOUNITS"
2617       include "COMMON.GEO"
2618       include "COMMON.SBRIDGE"
2619       double precision thetnorm
2620       integer iprot,i,j,ii
2621       do i=1,nres_zs(iprot)
2622         do j=1,3
2623           c_zs(j,i,ii,iprot)=c(j,i)
2624         enddo
2625       enddo
2626       do i=nnt_zs(iprot),nct_zs(iprot)
2627         do j=1,3
2628           c_zs(j,i+nres,ii,iprot)=c(j,i+nres)
2629         enddo
2630       enddo
2631 c 5/7/02 AL: store sbridge info
2632       nss_zs(ii,iprot)=nss
2633       do i=1,nss
2634         ihpb_zs(i,ii,iprot)=ihpb(i)
2635         jhpb_zs(i,ii,iprot)=jhpb(i)
2636       enddo
2637       return
2638       end
2639 c------------------------------------------------------------------------------
2640       subroutine restore_ccoords(iprot,ii)
2641       implicit none
2642       include "DIMENSIONS"
2643       include "DIMENSIONS.ZSCOPT"
2644       include "COMMON.INTERACT"
2645       include "COMMON.VAR"
2646       include "COMMON.ALLPROT"
2647       include "COMMON.SBRIDGE"
2648       include "COMMON.CHAIN"
2649       include "COMMON.IOUNITS"
2650       integer iprot,i,j,ii
2651       do i=1,nres_zs(iprot)
2652         do j=1,3
2653           c(j,i)=c_zs(j,i,ii,iprot)
2654         enddo
2655       enddo
2656       do i=nnt_zs(iprot),nct_zs(iprot)
2657         do j=1,3
2658           c(j,i+nres)=c_zs(j,i+nres,ii,iprot)
2659         enddo
2660       enddo
2661 c 5/7/02 AL: restore sbridge info
2662       nss=nss_zs(ii,iprot)
2663       do i=1,nss
2664         ihpb(i)=ihpb_zs(i,ii,iprot)+nres
2665         jhpb(i)=jhpb_zs(i,ii,iprot)+nres
2666 c        forcon(i)=fbr
2667 c        dhpb(i)=dbr
2668       enddo
2669 #ifdef DEBUG
2670         write (iout,*) "restore_ccoords",ii
2671         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2672         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2673         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2674         call flush(iout)
2675 #endif
2676       return
2677       end
2678 c------------------------------------------------------------------------------
2679       subroutine card_concat(card,to_upper)
2680       implicit none
2681       include 'DIMENSIONS.ZSCOPT'
2682       include "COMMON.IOUNITS"
2683       character*(*) card
2684       character*80 karta,ucase
2685       logical to_upper
2686       integer ilen
2687       external ilen
2688       read (inp,'(a)') karta
2689       if (to_upper) karta=ucase(karta)
2690       card=' '
2691       do while (karta(80:80).eq.'&')
2692         card=card(:ilen(card)+1)//karta(:79)
2693         read (inp,'(a)') karta
2694         if (to_upper) karta=ucase(karta)
2695       enddo
2696       card=card(:ilen(card)+1)//karta
2697       return
2698       end
2699 c------------------------------------------------------------------------------
2700       subroutine readi(rekord,lancuch,wartosc,default)
2701       implicit none
2702       character*(*) rekord,lancuch
2703       integer wartosc,default
2704       integer ilen,iread
2705       external ilen
2706       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2707       if (iread.eq.0) then
2708         wartosc=default
2709         return
2710       endif
2711       iread=iread+ilen(lancuch)+1
2712       read (rekord(iread:),*) wartosc
2713       return
2714       end
2715 c----------------------------------------------------------------------------
2716       subroutine reada(rekord,lancuch,wartosc,default)
2717       implicit none
2718       character*(*) rekord,lancuch
2719       character*80 aux
2720       double precision wartosc,default
2721       integer ilen,iread
2722       external ilen
2723       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2724       if (iread.eq.0) then
2725         wartosc=default
2726         return
2727       endif
2728       iread=iread+ilen(lancuch)+1
2729       read (rekord(iread:),*) wartosc
2730       return
2731       end
2732 c----------------------------------------------------------------------------
2733       subroutine multreadi(rekord,lancuch,tablica,dim,default)
2734       implicit none
2735       integer dim,i
2736       integer tablica(dim),default
2737       character*(*) rekord,lancuch
2738       character*80 aux
2739       integer ilen,iread
2740       external ilen
2741       do i=1,dim
2742         tablica(i)=default
2743       enddo
2744       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2745       if (iread.eq.0) return
2746       iread=iread+ilen(lancuch)+1
2747       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2748    10 return
2749       end
2750 c----------------------------------------------------------------------------
2751       subroutine multreada(rekord,lancuch,tablica,dim,default)
2752       implicit none
2753       integer dim,i
2754       double precision tablica(dim),default
2755       character*(*) rekord,lancuch
2756       character*80 aux
2757       integer ilen,iread
2758       external ilen
2759       do i=1,dim
2760         tablica(i)=default
2761       enddo
2762       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2763       if (iread.eq.0) return
2764       iread=iread+ilen(lancuch)+1
2765       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2766    10 return
2767       end
2768 c----------------------------------------------------------------------------
2769       subroutine reads(rekord,lancuch,wartosc,default)
2770       implicit none
2771       character*(*) rekord,lancuch,wartosc,default
2772       character*80 aux
2773       integer ilen,lenlan,lenrec,iread,ireade
2774       external ilen
2775       logical iblnk
2776       external iblnk
2777       lenlan=ilen(lancuch)
2778       lenrec=ilen(rekord)
2779       iread=index(rekord,lancuch(:lenlan)//"=")
2780 c      print *,"rekord",rekord," lancuch",lancuch
2781 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2782       if (iread.eq.0) then
2783         wartosc=default
2784         return
2785       endif
2786       iread=iread+lenlan+1
2787 c      print *,"iread",iread
2788 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2789       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2790         iread=iread+1
2791 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2792       enddo
2793 c      print *,"iread",iread
2794       if (iread.gt.lenrec) then
2795          wartosc=default
2796         return
2797       endif
2798       ireade=iread+1
2799 c      print *,"ireade",ireade
2800       do while (ireade.lt.lenrec .and.
2801      &   .not.iblnk(rekord(ireade:ireade)))
2802         ireade=ireade+1
2803       enddo
2804       wartosc=rekord(iread:ireade)
2805       return
2806       end
2807 c----------------------------------------------------------------------------
2808       subroutine multreads(rekord,lancuch,tablica,dim,default)
2809       implicit none
2810       integer dim,i
2811       character*(*) rekord,lancuch,tablica(dim),default
2812       character*80 aux
2813       integer ilen,lenlan,lenrec,iread,ireade
2814       external ilen
2815       logical iblnk
2816       external iblnk
2817       do i=1,dim
2818         tablica(i)=default
2819       enddo
2820       lenlan=ilen(lancuch)
2821       lenrec=ilen(rekord)
2822       iread=index(rekord,lancuch(:lenlan)//"=")
2823 c      print *,"rekord",rekord," lancuch",lancuch
2824 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2825       if (iread.eq.0) return
2826       iread=iread+lenlan+1
2827       do i=1,dim
2828 c      print *,"iread",iread
2829 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2830       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2831         iread=iread+1
2832 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2833       enddo
2834 c      print *,"iread",iread
2835       if (iread.gt.lenrec) return
2836       ireade=iread+1
2837 c      print *,"ireade",ireade
2838       do while (ireade.lt.lenrec .and.
2839      &   .not.iblnk(rekord(ireade:ireade)))
2840         ireade=ireade+1
2841       enddo
2842       tablica(i)=rekord(iread:ireade)
2843       iread=ireade+1
2844       enddo
2845       end
2846 c----------------------------------------------------------------------------
2847       subroutine split_string(rekord,tablica,dim,nsub)
2848       implicit none
2849       integer dim,nsub,i,ii,ll,kk
2850       character*(*) tablica(dim)
2851       character*(*) rekord
2852       integer ilen
2853       external ilen
2854       do i=1,dim
2855         tablica(i)=" "
2856       enddo
2857       ii=1
2858       ll = ilen(rekord)
2859       nsub=0
2860       do i=1,dim
2861 C Find the start of term name
2862         kk = 0
2863         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
2864           ii = ii+1
2865         enddo
2866 C Parse the name into TABLICA(i) until blank found
2867         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
2868           kk = kk+1
2869           tablica(i)(kk:kk)=rekord(ii:ii)
2870           ii = ii+1
2871         enddo
2872         if (kk.gt.0) nsub=nsub+1
2873         if (ii.gt.ll) return
2874       enddo
2875       return
2876       end
2877 c-------------------------------------------------------------------------------
2878       integer function iroof(n,m)
2879       ii = n/m
2880       if (ii*m .lt. n) ii=ii+1
2881       iroof = ii
2882       return
2883       end