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