source/maxlik/src-Fmatch_safe/readrtns_MP.F

   1       subroutine read_general_data(*)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5 #ifdef MPI
   6       include "mpif.h"
   7       include "COMMON.MPI"
   8       integer ierror,kolor,klucz
   9 #endif
  10       include "COMMON.WEIGHTS"
  11       include "COMMON.NAMES"
  12       include "COMMON.VMCPAR"
  13       include "COMMON.TORSION"
  14       include "COMMON.INTERACT"
  15       include "COMMON.ENERGIES"
  16       include "COMMON.MINPAR"
  17       include "COMMON.IOUNITS"
  18       include "COMMON.TIME1"
  19       include "COMMON.PROTFILES"
  20       include "COMMON.CHAIN"
  21       include "COMMON.CLASSES"
  22       include "COMMON.THERMAL"
  23       include "COMMON.FFIELD"
  24       include "COMMON.OPTIM"
  25       include "COMMON.CONTROL"
  26       include "COMMON.SCCOR"
  27       include "COMMON.SPLITELE"
  28       character*800 controlcard
  29       integer i,j,k,l,ii,n_ene_found,ist1,iet1,ist2,iet2,ls,le
  30       integer ind,itype1,itype2,itypf,itypsc,itypp,itypt1,itypt2
  31       integer ilen,lenpot,lenpre
  32       external ilen
  33       character*4 liczba,liczba1
  34 #ifndef ISNAN
  35       external proc_proc
  36 #endif
  37       character*16 ucase
  38       character*16 key
  39       external ucase
  40       double precision xchuj,weitemp,weitemp_low,weitemp_up
  41       character*80 item(7)
  42       integer nitem
  43       integer rescode
  44
  45 c      write (iout,*) "Enter read_general_data"
  46 c      call flush(iout)
  47
  48       lenpot=ilen(pot)
  49       lenpre=ilen(prefix)
  50 c      do i=1,max_ene
  51 c        ename(i)=" "
  52 c      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       out_forces=index(controlcard,"OUT_FORCES").gt.0
 113       write (iout,*) "OUT_FORCES ",out_forces
 114       unres_pdb = index(controlcard,"UNRES_PDB").gt.0
 115 c      write (iout,*) "UNRES_PDB ",unres_pdb
 116 c Energy calculation settings section
 117       call readi(controlcard,'TORMODE',tor_mode,0)
 118 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 119        call reada(controlcard,'BOXX',boxxsize,100.0d0)
 120        call reada(controlcard,'BOXY',boxysize,100.0d0)
 121        call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 122 c      print *,me," checkpoint 2"
 123 c Cutoff range for interactions
 124       call reada(controlcard,"R_CUT",r_cut,15.0d0)
 125       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 126       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 127       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 128       if (lipthick.gt.0.0d0) then
 129        bordliptop=(boxzsize+lipthick)/2.0
 130        bordlipbot=bordliptop-lipthick
 131 C      endif
 132       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
 133      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 134       buflipbot=bordlipbot+lipbufthick
 135       bufliptop=bordliptop-lipbufthick
 136       if ((lipbufthick*2.0d0).gt.lipthick)
 137      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 138       endif
 139 c      write(iout,*) "bordliptop=",bordliptop
 140 c      write(iout,*) "bordlipbot=",bordlipbot
 141 c      write(iout,*) "bufliptop=",bufliptop
 142 c      write(iout,*) "buflipbot=",buflipbot
 143       call readi(controlcard,'SYM',symetr,1)
 144 c      print *,me," checkpoint 3"
 145 c      write (iout,*) "n_ene",n_ene
 146       call flush(iout)
 147 c Energy-term-weights section
 148       wname(4)="WCORRH"
 149 C Read third record: weights
 150       do i=1,max_ene
 151         ww(i)=ww0(i)
 152       enddo
 153 c      print *,me," checkpoint 4"
 154 C Read fourth record: masks
 155       call card_concat(controlcard,.true.)
 156 c      write (iout,*) "card2",controlcard
 157       do i=1,n_ene
 158         key = "MASK_"//wname(i)(2:ilen(wname(i)))
 159         call readi(controlcard,key(:ilen(key)),imask(i),0)
 160       enddo
 161 c      print *,me," checkpoint 5"
 162 C Read fifth record: lower limits of weights
 163       call card_concat(controlcard,.true.)
 164 c      write (iout,*) "card3",controlcard
 165       do i=1,n_ene
 166         key = "WLOW_"//wname(i)(2:ilen(wname(i)))
 167         call reada(controlcard,key(:ilen(key)),ww_low(i),ww(i))
 168       enddo
 169 C Read sixth record: upper limits of weights
 170       call card_concat(controlcard,.true.)
 171 c      write (iout,*) "card4",controlcard
 172       do i=1,n_ene
 173         key = "WUP_"//wname(i)(2:ilen(wname(i)))
 174         call reada(controlcard,key(:ilen(key)),ww_up(i),ww(i))
 175       enddo
 176 C Read seventh record: VMC parameters
 177       call card_concat(controlcard,.true.)
 178 c      write (iout,*) "card5",controlcard
 179       call readi(controlcard,"MODE",mode,3)
 180       call readi(controlcard,"NSCANCYCLE",nscancycle,3)
 181       call readi(controlcard,"MAXSTEP_SCAN",maxstep_scan,50)
 182       call reada(controlcard,"RSTEP_SCAN",step_scan,1.0d-1)
 183       call readi(controlcard,"READ_STAT",read_stat,3)
 184       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 185 c      init_ene = index(controlcard,"INIT_ENE").gt.0 .and. read_stat.gt.1
 186       init_ene = .true.
 187       call readi(controlcard,"NMCM",nmcm,0)
 188       call readi(controlcard,"MAXSCAN",maxscan,0)
 189       call readi(controlcard,"MAXMIN",maxmin,100)
 190       call readi(controlcard,"MAXFUN",maxfun,100)
 191       call reada(controlcard,"TOLF",tolf,1.0d-6)
 192       call reada(controlcard,"RTOLF",rtolf,1.0d-6)
 193       out_minim=0
 194       if (index(controlcard,"OUT_MINIM").gt.0) out_minim=iout
 195       print_ini=0
 196       if (index(controlcard,"PRINT_INI").gt.0) print_ini=1
 197       print_fin=0
 198       if (index(controlcard,"PRINT_FIN").gt.0) print_fin=1
 199       print_stat=0
 200       if (index(controlcard,"PRINT_STAT").gt.0) print_stat=1
 201       call reada(controlcard,"RSTIME",rstime,9.0d2)
 202       call reads(controlcard,"MINIMIZER",minimizer,"SUMSL")
 203       call readi(controlcard,"OPT_MODE",opt_mode,0)
 204       mod_other_params=index(controlcard,"OPTIMIZE_OTHER").gt.0
 205       if (read_stat.eq.0 .and. mod_other_params) then
 206         write (iout,*) "Error: only optimization of energy-term",
 207      &    " weights can be performed with READ_STAT=",read_stat
 208         call flush(iout)
 209         return1
 210       endif
 211       if (index(controlcard,"RESTART").gt.0) then
 212         restart_flag=.true.
 213       else
 214         restart_flag=.false.
 215       endif
 216 c      print *,me," checkpoint 6"
 217       return
 218       end
 219 c-------------------------------------------------------------------------------------------------
 220       subroutine read_pmf_data(*)
 221       implicit none
 222       include "DIMENSIONS"
 223       include "DIMENSIONS.ZSCOPT"
 224 #ifdef MPI
 225       include "mpif.h"
 226       include "COMMON.MPI"
 227       integer ierror,kolor,klucz
 228 #endif
 229       include "COMMON.WEIGHTS"
 230       include "COMMON.NAMES"
 231       include "COMMON.VMCPAR"
 232       include "COMMON.TORSION"
 233       include "COMMON.INTERACT"
 234       include "COMMON.ENERGIES"
 235       include "COMMON.MINPAR"
 236       include "COMMON.IOUNITS"
 237       include "COMMON.TIME1"
 238       include "COMMON.PROTFILES"
 239       include "COMMON.CHAIN"
 240       include "COMMON.CLASSES"
 241       include "COMMON.THERMAL"
 242       include "COMMON.FFIELD"
 243       include "COMMON.OPTIM"
 244       include "COMMON.CONTROL"
 245       include "COMMON.SCCOR"
 246       include "COMMON.SPLITELE"
 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       write (iout,*) "TORSION_PMF", TORSION_PMF," TURN_PMF ",TURN_PMF,
 270      &  " EELLO_PMF",EELLO_PMF
 271       call reada(controlcard,"WELLO_PMF",wello_pmf,1.0d0)
 272       call reada(controlcard,"WELLO_PMF_LOW",wello_pmf_low,0.0d0)
 273       call reada(controlcard,"WELLO_PMF_UP",wello_pmf_up,5.0d0)
 274       call reada(controlcard,"WTURN_PMF",wturn_pmf,1.0d0)
 275       call reada(controlcard,"WTURN_PMF_LOW",wturn_pmf_low,0.0d0)
 276       call reada(controlcard,"WTURN_PMF_UP",wturn_pmf_up,5.0d0)
 277       call reada(controlcard,"WPMF",wpmf,1.0d0)
 278       write (iout,*) "nloctyp",nloctyp
 279       call multreada(controlcard,"E0",e0(0,-nloctyp+1),
 280      &  (ntyp+1)*(2*nloctyp-1),0.0d0)
 281       call multreada(controlcard,"E0_LOW",e0_low(0,-nloctyp+1),
 282      &  (ntyp+1)*(2*nloctyp-1),
 283      &  -1.0d2)
 284       call multreada(controlcard,"E0_UP",e0_up(0,-nloctyp+1),
 285      &  (ntyp+1)*(2*nloctyp-1),
 286      &  1.0d2)
 287       write (iout,*) "WTURN_PMF",WTURN_PMF,WTURN_PMF_LOW,WTURN_PMF_UP
 288       write (iout,*) "WELLO_PMF",WELLO_PMF,WELLO_PMF_LOW,WELLO_PMF_UP
 289       write (iout,*) "E0"
 290       do i=0,2
 291         do j=-2,2
 292           write (iout,"(2i5,3f15.5)") i,j,e0(i,j),e0_low(i,j),e0_up(i,j)
 293         enddo
 294       enddo
 295 #endif
 296       return
 297       end
 298 c-----------------------------------------------------------------------
 299       subroutine read_optim_parm(*)
 300       implicit none
 301       include "DIMENSIONS"
 302       include "DIMENSIONS.ZSCOPT"
 303 #ifdef MPI
 304       include "mpif.h"
 305       include "COMMON.MPI"
 306       integer ierror,kolor,klucz
 307 #endif
 308       include "COMMON.WEIGHTS"
 309       include "COMMON.NAMES"
 310       include "COMMON.VMCPAR"
 311       include "COMMON.TORSION"
 312       include "COMMON.LOCAL"
 313       include "COMMON.INTERACT"
 314       include "COMMON.ENERGIES"
 315       include "COMMON.MINPAR"
 316       include "COMMON.IOUNITS"
 317       include "COMMON.TIME1"
 318       include "COMMON.PROTFILES"
 319       include "COMMON.CHAIN"
 320       include "COMMON.CLASSES"
 321       include "COMMON.THERMAL"
 322       include "COMMON.FFIELD"
 323       include "COMMON.OPTIM"
 324       include "COMMON.CONTROL"
 325       include "COMMON.SCCOR"
 326       character*800 controlcard
 327       character*4 reskind
 328       integer i,j,k,l,ii,n_ene_found,ist1,iet1,ist2,iet2,ls,le
 329       integer ind,ind1,ind2,itype1,itype2,itypf,itypsc,itypp,
 330      &   itypt1,itypt2,masktemp,iblock,iaux,itypa
 331       integer ilen,lenpot,lenpre
 332       external ilen
 333       double precision aux,vb_low,vb_up,vb_rel
 334       character*4 liczba,liczba1
 335 #ifndef ISNAN
 336       external proc_proc
 337 #endif
 338       character*16 ucase
 339       character*16 key
 340       external ucase
 341       double precision xchuj,weitemp,weitemp_low,weitemp_up
 342       character*80 item(7)
 343       character*3 typf,typa
 344       integer nitem
 345       integer rescode
 346       integer ind_shield /25/
 347
 348       lenpot=ilen(pot)
 349       lenpre=ilen(prefix)
 350       write (iout,*) "MOD_OTHER_PARAMS ",mod_other_params
 351       do iblock=1,2
 352       do itypt1=-ntyp,ntyp
 353         do itypt2=-ntyp,ntyp
 354           mask_tor(0,itypt1,itypt2,iblock)=0
 355           weitor(0,itypt1,itypt2,iblock)=1.0d0
 356           weitor_low(0,itypt1,itypt2,iblock)=1.0d0
 357           weitor_up(0,itypt1,itypt2,iblock)=1.0d0
 358         enddo
 359       enddo
 360       do itypt1=-ntyp,ntyp
 361         do itypt2=-ntyp,ntyp
 362           do l=0,3
 363             mask_tor(l,itypt1,itypt2,iblock)=0
 364             weitor(l,itypt1,itypt2,iblock)=1.0
 365             weitor_low(l,itypt1,itypt2,iblock)=1.0
 366             weitor_up(l,itypt1,itypt2,iblock)=1.0
 367           enddo
 368         enddo
 369       enddo
 370       enddo
 371 #ifdef DEBUG
 372           write (iout,*) "ntyp",ntyp
 373           do itypt1=-ntyp,ntyp
 374             do itypt2=-ntyp,ntyp
 375               write (iout,*) "itypt1",itypt1," itypt2",itypt2,
 376      &      "weitor",weitor(0,itypt1,itypt2,1),eitor(0,itypt1,itypt2,2)
 377             enddo
 378           enddo
 379 #endif
 380       if (mod_other_params) then
 381 c        write (*,*)
 382 c     &"Internal parameters of UNRES energy components will be optimized"
 383         call card_concat(controlcard,.true.)
 384         write (iout,*) "mod_side ",controlcard
 385         call flush(iout)
 386         mod_side = (index(controlcard,"MOD_SIDE").gt.0)
 387         if (mod_side) then
 388           nsingle_sc=0
 389           npair_sc=0
 390           call card_concat(controlcard,.true.)
 391           do while ( index(controlcard,"END").eq.0 )
 392             call split_string(controlcard,item(1),4,nitem)
 393             if (nitem.eq.1 .or. item(2)(:1).eq."*") then
 394               nsingle_sc=nsingle_sc+1
 395               ityp_ssc(nsingle_sc)=rescode(1,item(1),0)
 396               if (nitem.lt.3) then
 397                 epss_low(nsingle_sc)=0.02d0
 398               else
 399                 read (item(3),*) epss_low(nsingle_sc)
 400               endif
 401               if (nitem.lt.4) then
 402                 epss_up(nsingle_sc)=0.0d0
 403               else
 404                 read (item(4),*) epss_up(nsingle_sc)
 405               endif
 406             else
 407               npair_sc=npair_sc+1
 408               ityp_psc(1,npair_sc)=rescode(1,item(1),0)
 409               ityp_psc(2,npair_sc)=rescode(1,item(2),0)
 410               if (nitem.lt.3) then
 411                 epsp_low(npair_sc)=0.02d0
 412               else
 413                 read (item(3),*) epsp_low(npair_sc)
 414               endif
 415               if (nitem.lt.4) then
 416                 epsp_up(npair_sc)=0.0d0
 417               else
 418                 read (item(4),*) epsp_up(npair_sc)
 419               endif
 420             endif
 421             call card_concat(controlcard,.true.)
 422           enddo
 423           if (nsingle_sc+npair_sc.eq.0) mod_side=.false.
 424           call card_concat(controlcard,.true.)
 425         endif
 426         mod_side_other=
 427      &    index(controlcard,"MOD_SIDE_OTHER").gt.0
 428         write (iout,*) "mod_side_other ",controlcard,mod_side_other
 429         if (mod_side_other) then
 430           mod_side_other=.false.
 431           do i=1,ntyp
 432             do j=1,5
 433               mask_sigma(j,i)=0
 434             enddo
 435           enddo
 436           call card_concat(controlcard,.true.)
 437 c          write (iout,*) "mod_side_oter ",controlcard
 438           do while ( index(controlcard,"END").eq.0 )
 439            call reads(controlcard,"RESKIND",reskind,"   ")
 440            itypsc=rescode(1,reskind,0)
 441            if (itypsc.lt.1 .or. itypsc.gt.20) then
 442               write (iout,*) "Error in SC optimization data;",
 443      &          " SC type must not be dummy, type is" ,restyp," ",itypsc
 444               write (*,*) "Error in SC optimization data;",
 445      &          " SC type must not be dummy, type is" ,restyp," ",itypsc
 446               return1
 447            endif
 448 c Sigma0 appear in all potentials
 449            call readi(controlcard,"MASK_SIGMA0",mask_sigma(1,itypsc),0)
 450 c Anisotropy parameters only in BP, GB, and GBV potentials
 451            if (ipot.gt.2) then
 452            call readi(controlcard,"MASK_SIGMAII",mask_sigma(2,itypsc),0)
 453            call readi(controlcard,"MASK_CHIP",mask_sigma(3,itypsc),0)
 454            call readi(controlcard,"MASK_ALP",mask_sigma(4,itypsc),0)
 455            call reada(controlcard,"SIGMA0_LOW",sigma_low(1,itypsc),0.d0)
 456            call reada(controlcard,"SIGMA0_UP",sigma_up(1,itypsc),0.0d0)
 457            call reada(controlcard,"SIGMAII_LOW",sigma_low(2,itypsc),
 458      &        0.0d0)
 459            call reada(controlcard,"SIGMAII_UP",sigma_up(2,itypsc),0.0d0)
 460            call reada(controlcard,"CHIP_LOW",sigma_low(3,itypsc),0.d0)
 461            call reada(controlcard,"CHIP_UP",sigma_up(3,itypsc),0.0d0)
 462            call reada(controlcard,"ALP_LOW",sigma_low(4,itypsc),0.d0)
 463            call reada(controlcard,"ALP_UP",sigma_up(4,itypsc),0.0d0)
 464            endif
 465 c r0 only in LJ and LJK potentials
 466            if (ipot.eq.2 .or. ipot.eq.5) then
 467            call readi(controlcard,"MASK_R0",mask_sigma(5,itypsc),0)
 468            call reada(controlcard,"R0_LOW",sigma_low(5,itypsc),0.d0)
 469            call reada(controlcard,"R0_UP",sigma_up(5,itypsc),0.0d0)
 470            endif
 471            do k=1,5
 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,5
 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,5)
 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         mod_ang=.false.
 698         do i=0,nthetyp-1
 699           mask_ang(i)=0
 700         enddo
 701
 702         call card_concat(controlcard,.true.)
 703         do while (index(controlcard,'END').eq.0)
 704           mod_ang=.true.
 705           write (iout,'(a)') "angle: ",controlcard
 706           call reads(controlcard,"TYPE",typa,"   ")
 707           itypa=rescode(1,typa,0)
 708 c          write (iout,*) "TYPA ",typa," itypa",itypa
 709           if (itypa.eq.0 .or. itypa.gt.ntyp) then
 710             write (*,*) "Error specifying angle parms"
 711             stop
 712           endif
 713           itypa=itortyp(itypa)
 714           mask_ang(itypa)=1
 715           write (iout,*) "bend type",itypa
 716           call reada(controlcard,"VB_LOW",vb_low,-1.0d5)
 717           call reada(controlcard,"VB_UP",vb_up,1.0d5)
 718           call reada(controlcard,"VB_REL",vb_rel,0.0d0)
 719           write (iout,*) "VB_LOW",vb_low," VB_UP",vb_up," VB_REL",vb_rel
 720           do i=1,nbend_kcc_TB(itypa)
 721             if (vb_rel.gt.0) then
 722               write (iout,*) "v1bend_chyb=",v1bend_chyb(i,itypa)
 723               v1bend_low(i,itypa)=v1bend_chyb(i,itypa)*(1.0d0
 724      &           -dsign(vb_rel,v1bend_chyb(i,itypa)))
 725               v1bend_up(i,itypa)=v1bend_chyb(i,itypa)*(1.0d0
 726      &           +dsign(vb_rel,v1bend_chyb(i,itypa)))
 727             else
 728               v1bend_low(i,itypa)=vb_low
 729               v1bend_up(i,itypa)=vb_up
 730             endif
 731           enddo
 732           call card_concat(controlcard,.true.)
 733
 734         enddo
 735         call card_concat(controlcard,.true.)
 736
 737         write (iout,*) "Boundaries of angle potential coefficients"
 738         write (iout,'(3a)') "Index","      Low bound","       Up bound"
 739         do i=0,nthetyp-1
 740           if (mask_ang(i).eq.0) cycle
 741           write (iout,'(2a)') 'Type ',restyp(iloctyp(i))
 742           do j=1,nbend_kcc_TB(i)
 743           write (iout,'(i5,2f15.1)') j,v1bend_low(j,i),v1bend_up(j,i)
 744           enddo
 745         enddo
 746
 747         endif
 748
 749         write (iout,*) "mod_fourier ",controlcard
 750         call flush(iout)
 751         mod_fourier(nloctyp)=index(controlcard,"MOD_FOURIER")
 752 #ifdef NEWCORR
 753         if (mod_fourier(nloctyp).gt.0) then
 754           mod_fourier(nloctyp)=0
 755           do i=0,nloctyp-1
 756             mod_fourier(i)=0
 757             do ii=1,3
 758               do j=1,2
 759                 mask_bnew1(ii,j,i)=0
 760                 mask_bnew2(ii,j,i)=0
 761                 mask_ccnew(ii,j,i)=0
 762                 mask_ddnew(ii,j,i)=0
 763               enddo
 764             enddo
 765             do ii=1,2
 766               do j=1,2
 767                 do k=1,2
 768                   mask_eenew(ii,j,k,i)=0
 769                 enddo
 770               enddo
 771             enddo
 772             do ii=1,3
 773               mask_e0new(ii,i)=0
 774             enddo
 775           enddo
 776           call card_concat(controlcard,.true.)
 777           do while ( index(controlcard,"END").eq.0 )
 778 c            write(iout,*) controlcard
 779             call reads(controlcard,"TYPF",typf,"   ")
 780             itypf=rescode(1,typf,0)
 781 c            write (iout,*) "TYPF ",typf," itypf",itypf
 782             if (itypf.eq.0 .or. itypf.gt.ntyp) then
 783               write (*,*) "Error specifying fourier parms"
 784               stop
 785             endif
 786             itypf=itype2loc(itypf)
 787             write (iout,*) "local type",itypf
 788
 789             if (index(controlcard,"B1_LOW").gt.0) then
 790
 791             call readi(controlcard,"IND",ind,0)
 792             call readi(controlcard,"COEFF",ii,0)
 793             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 794               write (iout,*)
 795      &          "Error specifying B1, components undefined",ind,ii
 796               stop
 797             endif
 798             mask_bnew1(ii,ind,itypf)=1
 799             call reada(controlcard,"B1_LOW",bnew1_low(ii,ind,itypf),
 800      &        0.1d0)
 801             call reada(controlcard,"B1_UP",bnew1_up(ii,ind,itypf),
 802      &        0.0d0)
 803             mod_fourier(itypf)=mod_fourier(itypf)
 804      &        +mask_bnew1(ii,ind,itypf)
 805
 806             else if (index(controlcard,"B2_LOW").gt.0) then
 807
 808             mask_bnew2(ii,ind,itypf)=1
 809             call readi(controlcard,"IND",ind,0)
 810             call readi(controlcard,"COEFF",ii,0)
 811             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 812               write (iout,*)
 813      &          "Error specifying B2, components undefined",ind,ii
 814               stop
 815             endif
 816             call reada(controlcard,"B2_LOW",bnew2_low(ii,ind,itypf),
 817      &        0.1d0)
 818             call reada(controlcard,"B2_UP",bnew2_up(ii,ind,itypf),
 819      &        0.0d0)
 820             mod_fourier(itypf)=mod_fourier(itypf)
 821      &        +mask_bnew2(ii,ind,itypf)
 822
 823             else if (index(controlcard,"C_LOW").gt.0) then
 824
 825             call readi(controlcard,"IND",ind,0)
 826             call readi(controlcard,"COEFF",ii,0)
 827             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 828               write (iout,*)
 829      &          "Error specifying C, components undefined",ind,ii
 830               stop
 831             endif
 832             mask_ccnew(ii,ind,itypf)=1
 833             call reada(controlcard,"C_LOW",ccnew_low(ii,ind,itypf),.1d0)
 834             call reada(controlcard,"C_UP",ccnew_up(ii,ind,itypf),0.0d0)
 835             mod_fourier(itypf)=mod_fourier(itypf)
 836      &        +mask_ccnew(ii,ind,itypf)
 837
 838             else if (index(controlcard,"D_LOW").gt.0) then
 839
 840             call readi(controlcard,"IND",ind,0)
 841             call readi(controlcard,"COEFF",ii,0)
 842             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 843               write (iout,*)
 844      &          "Error specifying D, components undefined",ind,ii
 845               stop
 846             endif
 847             mask_ddnew(ii,ind,itypf)=1
 848             call reada(controlcard,"D_LOW",ddnew_low(ii,ind,itypf),
 849      &        0.1d0)
 850             call reada(controlcard,"D_UP",ddnew_up(ii,ind,itypf),
 851      &        0.0d0)
 852             mod_fourier(itypf)=mod_fourier(itypf)
 853      &        +mask_ddnew(ii,ind,itypf)
 854
 855             else if (index(controlcard,"E0_LOW").gt.0) then
 856
 857             call readi(controlcard,"COEFF",ii,0)
 858             if (ii.eq.0 .or. ii.gt.3) then
 859               write (iout,*)
 860      &          "Error specifying E0, components undefined",ii
 861               stop
 862             endif
 863             mask_e0new(ii,itypf)=1
 864             call reada(controlcard,"E0_LOW",e0new_low(ii,itypf),
 865      &        0.1d0)
 866             call reada(controlcard,"E0_UP",e0new_up(ii,itypf),
 867      &        0.0d0)
 868             mod_fourier(itypf)=mod_fourier(itypf)
 869      &        +mask_e0new(ii,itypf)
 870
 871             else if (index(controlcard,"E_LOW").gt.0) then
 872
 873             call readi(controlcard,"IND1",ind1,0)
 874             call readi(controlcard,"IND2",ind2,0)
 875             call readi(controlcard,"COEFF",ii,0)
 876             if (ind1.eq.0 .or. ind2.eq.0 .or. ii.eq.0 .or. ii.gt.2) then
 877               write (iout,*)
 878      &          "Error specifying E, components undefined",ind1,ind2,ii
 879               stop
 880             endif
 881             mask_eenew(ii,ind1,ind2,itypf)=1
 882             call reada(controlcard,"E_LOW",
 883      &        eenew_low(ii,ind1,ind2,itypf),0.1d0)
 884             call reada(controlcard,"E_UP",
 885      &        eenew_up(ii,ind1,ind2,itypf),0.0d0)
 886             mod_fourier(itypf)=mod_fourier(itypf)
 887      &        +mask_eenew(ii,ind1,ind2,itypf)
 888             endif
 889             call card_concat(controlcard,.true.)
 890           enddo
 891           call card_concat(controlcard,.true.)
 892           write (iout,*) "mod_fouriertor card ",controlcard
 893           mod_fouriertor(nloctyp)=index(controlcard,"MOD_FOURIERTOR")
 894           write (iout,*) "mod_fouriertor value",mod_fouriertor(nloctyp)
 895           write (2,*) "SPLIT_FOURIERTOR",SPLIT_FOURIERTOR,
 896      &      " tor_mode",tor_mode
 897           IF (SPLIT_FOURIERTOR .and. tor_mode.eq.2
 898      &        .and. mod_fouriertor(nloctyp).gt.0) THEN
 899           do i=0,nloctyp-1
 900             mod_fouriertor(i)=0
 901             do ii=1,3
 902               do j=1,2
 903                 mask_bnew1tor(ii,j,i)=0
 904                 mask_bnew2tor(ii,j,i)=0
 905                 mask_ccnewtor(ii,j,i)=0
 906                 mask_ddnewtor(ii,j,i)=0
 907               enddo
 908             enddo
 909             do ii=1,2
 910               do j=1,2
 911                 do k=1,2
 912                   mask_eenewtor(ii,j,k,i)=0
 913                 enddo
 914               enddo
 915             enddo
 916             do ii=1,3
 917               mask_e0newtor(ii,i)=0
 918             enddo
 919           enddo
 920           call card_concat(controlcard,.true.)
 921           do while ( index(controlcard,"END").eq.0 )
 922 c            write(iout,*) controlcard
 923             call reads(controlcard,"TYPF",typf,"   ")
 924             itypf=rescode(1,typf,0)
 925 c            write (iout,*) "TYPF ",typf," itypf",itypf
 926             if (itypf.eq.0 .or. itypf.gt.ntyp) then
 927               write (*,*) "Error specifying fourier parms"
 928               stop
 929             endif
 930             itypf=itype2loc(itypf)
 931             write (iout,*) "local type",itypf
 932
 933             if (index(controlcard,"B1_LOW").gt.0) then
 934
 935             call readi(controlcard,"IND",ind,0)
 936             call readi(controlcard,"COEFF",ii,0)
 937             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 938               write (iout,*)
 939      &          "Error specifying B1, components undefined",ind,ii
 940               stop
 941             endif
 942             mask_bnew1tor(ii,ind,itypf)=1
 943             call reada(controlcard,"B1_LOW",bnew1tor_low(ii,ind,itypf),
 944      &        0.1d0)
 945             call reada(controlcard,"B1_UP",bnew1tor_up(ii,ind,itypf),
 946      &        0.0d0)
 947             mod_fouriertor(itypf)=mod_fouriertor(itypf)
 948      &        +mask_bnew1tor(ii,ind,itypf)
 949
 950             else if (index(controlcard,"B2_LOW").gt.0) then
 951
 952             mask_bnew2tor(ii,ind,itypf)=1
 953             call readi(controlcard,"IND",ind,0)
 954             call readi(controlcard,"COEFF",ii,0)
 955             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 956               write (iout,*)
 957      &          "Error specifying B2, components undefined",ind,ii
 958               stop
 959             endif
 960             call reada(controlcard,"B2_LOW",bnew2tor_low(ii,ind,itypf),
 961      &        0.1d0)
 962             call reada(controlcard,"B2_UP",bnew2tor_up(ii,ind,itypf),
 963      &        0.0d0)
 964             mod_fouriertor(itypf)=mod_fouriertor(itypf)
 965      &        +mask_bnew2tor(ii,ind,itypf)
 966
 967             else if (index(controlcard,"C_LOW").gt.0) then
 968
 969             call readi(controlcard,"IND",ind,0)
 970             call readi(controlcard,"COEFF",ii,0)
 971             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 972               write (iout,*)
 973      &          "Error specifying C, components undefined",ind,ii
 974               stop
 975             endif
 976             mask_ccnewtor(ii,ind,itypf)=1
 977             call reada(controlcard,"C_LOW",ccnewtor_low(ii,ind,itypf),
 978      &       .1d0)
 979             call reada(controlcard,"C_UP",ccnewtor_up(ii,ind,itypf),
 980      &       0.0d0)
 981             mod_fouriertor(itypf)=mod_fouriertor(itypf)
 982      &        +mask_ccnewtor(ii,ind,itypf)
 983
 984             else if (index(controlcard,"D_LOW").gt.0) then
 985
 986             call readi(controlcard,"IND",ind,0)
 987             call readi(controlcard,"COEFF",ii,0)
 988             if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
 989               write (iout,*)
 990      &          "Error specifying D, components undefined",ind,ii
 991               stop
 992             endif
 993             mask_ddnewtor(ii,ind,itypf)=1
 994             call reada(controlcard,"D_LOW",ddnewtor_low(ii,ind,itypf),
 995      &        0.1d0)
 996             call reada(controlcard,"D_UP",ddnewtor_up(ii,ind,itypf),
 997      &        0.0d0)
 998             mod_fouriertor(itypf)=mod_fouriertor(itypf)
 999      &        +mask_ddnewtor(ii,ind,itypf)
1000
1001             else if (index(controlcard,"E0_LOW").gt.0) then
1002
1003             call readi(controlcard,"COEFF",ii,0)
1004             if (ii.eq.0 .or. ii.gt.3) then
1005               write (iout,*)
1006      &          "Error specifying E0, components undefined",ii
1007               stop
1008             endif
1009             mask_e0newtor(ii,itypf)=1
1010             call reada(controlcard,"E0_LOW",e0newtor_low(ii,itypf),
1011      &        0.1d0)
1012             call reada(controlcard,"E0_UP",e0newtor_up(ii,itypf),
1013      &        0.0d0)
1014             mod_fouriertor(itypf)=mod_fouriertor(itypf)
1015      &        +mask_e0newtor(ii,itypf)
1016
1017             else if (index(controlcard,"E_LOW").gt.0) then
1018
1019             call readi(controlcard,"IND1",ind1,0)
1020             call readi(controlcard,"IND2",ind2,0)
1021             call readi(controlcard,"COEFF",ii,0)
1022             if (ind1.eq.0 .or. ind2.eq.0 .or. ii.eq.0 .or. ii.gt.2) then
1023               write (iout,*)
1024      &          "Error specifying E, components undefined",ind1,ind2,ii
1025               stop
1026             endif
1027             mask_eenewtor(ii,ind1,ind2,itypf)=1
1028             call reada(controlcard,"E_LOW",
1029      &        eenewtor_low(ii,ind1,ind2,itypf),0.1d0)
1030             call reada(controlcard,"E_UP",
1031      &        eenewtor_up(ii,ind1,ind2,itypf),0.0d0)
1032             mod_fouriertor(itypf)=mod_fouriertor(itypf)
1033      &        +mask_eenewtor(ii,ind1,ind2,itypf)
1034             endif
1035             call card_concat(controlcard,.true.)
1036           enddo
1037           call card_concat(controlcard,.true.)
1038           write (2,*) "End read FOURIERTOR ",controlcard
1039
1040           ENDIF ! SPLIT_FOURIERTOR
1041
1042         endif
1043         do itypf=0,nloctyp-1
1044           write (iout,*) "itypf",itypf," mod_fourier",
1045      &    mod_fourier(itypf)
1046           mod_fourier(nloctyp)=mod_fourier(nloctyp)
1047      &      +mod_fourier(itypf)
1048         enddo
1049         write (iout,*) "mod_fourier all",mod_fourier(nloctyp)
1050         do itypf=0,nloctyp-1
1051           write (iout,*) "itypf",itypf," mod_fouriertor",
1052      &    mod_fouriertor(itypf)
1053           mod_fouriertor(nloctyp)=mod_fouriertor(nloctyp)
1054      &      +mod_fouriertor(itypf)
1055         enddo
1056         write (iout,*) "mod_fouriertor all",mod_fouriertor(nloctyp)
1057 #else
1058         if (mod_fourier(nloctyp).gt.0) then
1059           call card_concat(controlcard,.true.)
1060           do while ( index(controlcard,"END").eq.0 )
1061             call readi(controlcard,"TYPF",typf,"   ")
1062             itypf=rescode(1,typf,0)
1063             if (itypf.eq.0 .or. itypf.gt.ntyp) then
1064               write (*,*) "Error specifying fourier parms"
1065               stop
1066             endif
1067             itypf=itype2loc(itypf)
1068             call readi(controlcard,"IND",ind,0)
1069             call reada(controlcard,"B_LOW",b_low(ind,itypf),0.1d0)
1070             call reada(controlcard,"B_UP",b_up(ind,itypf),0.0d0)
1071             mask_fourier(ind,itypf)=1
1072             mod_fourier(itypf)=mod_fourier(itypf)
1073      &         +mask_fourier(ind,itypf)
1074             mod_fourier(nloctyp)=mod_fourier(nloctyp)
1075      &         +mask_fourier(ind,itypf)
1076             call card_concat(controlcard,.true.)
1077           enddo
1078           call card_concat(controlcard,.true.)
1079         endif
1080         do itypf=0,nloctyp-1
1081           write (iout,*) "itypf",itypf," mod_fourier",mod_fourier(itypf)
1082           mod_fourier(nloctyp)=mod_fourier(nloctyp)+mod_fourier(itypf)
1083         enddo
1084         write (iout,*) "mod_fourier all",mod_fourier(nloctyp)
1085 #endif
1086         do i=1,2
1087           do j=1,2
1088             do k=1,2
1089               mask_elec(i,j,k)=0
1090             enddo
1091           enddo
1092         enddo
1093         mod_elec=index(controlcard,"MOD_ELEC").gt.0
1094         if (mod_elec) then
1095           mod_elec=.false.
1096           call card_concat(controlcard,.true.)
1097           do while ( index(controlcard,"END").eq.0 )
1098 c            write (iout,*) "controlcard ",controlcard
1099             call readi(controlcard,"TYPE1",itype1,0)
1100             call readi(controlcard,"TYPE2",itype2,0)
1101 c            write (iout,*) "itype1",itype1," itype2",itype2
1102             if (itype1.eq.0 .or. itype1.gt.2 .or. itype2.eq.0
1103      &           .or. itype2.gt.2) then
1104               write (iout,*) "Error specifying elec parms"
1105               stop
1106             endif
1107             if (index(controlcard,"EPP").gt.0) then
1108               mod_elec=.true.
1109               mask_elec(itype1,itype2,1)=1
1110               mask_elec(itype2,itype1,1)=1
1111               call reada(controlcard,"EPP_LOW",epp_low(itype1,itype2),
1112      &          0.1d0)
1113               epp_low(itype2,itype1)=epp_low(itype1,itype2)
1114               call reada(controlcard,"EPP_UP",epp_up(itype1,itype2),
1115      &          0.0d0)
1116               epp_up(itype2,itype1)=epp_up(itype1,itype2)
1117             endif
1118             if (index(controlcard,"RPP").gt.0) then
1119               mod_elec=.true.
1120               mask_elec(itype1,itype2,2)=1
1121               mask_elec(itype2,itype1,2)=1
1122               call reada(controlcard,"RPP_LOW",rpp_low(itype1,itype2),
1123      &          0.1d0)
1124               rpp_low(itype2,itype1)=rpp_low(itype1,itype2)
1125               call reada(controlcard,"RPP_UP",rpp_up(itype1,itype2),
1126      &          0.0d0)
1127               rpp_up(itype2,itype1)=rpp_up(itype1,itype2)
1128             endif
1129             if (index(controlcard,"ELPP6").gt.0) then
1130               mod_elec=.true.
1131               mask_elec(itype1,itype2,3)=1
1132               mask_elec(itype2,itype1,3)=1
1133               call reada(controlcard,"ELPP6_LOW",
1134      &          elpp6_low(itype1,itype2),0.1d0)
1135               elpp6_low(itype2,itype1)=elpp6_low(itype1,itype2)
1136               call reada(controlcard,"ELPP6_UP",
1137      &          elpp6_up(itype1,itype2),0.0d0)
1138               elpp6_up(itype2,itype1)=elpp6_up(itype1,itype2)
1139             endif
1140             if (index(controlcard,"ELPP3").gt.0) then
1141               mod_elec=.true.
1142               mask_elec(itype1,itype2,4)=1
1143               mask_elec(itype2,itype1,4)=1
1144               call reada(controlcard,"ELPP3_LOW",
1145      &          elpp3_low(itype1,itype2),0.1d0)
1146               elpp3_low(itype2,itype1)=elpp3_low(itype1,itype2)
1147               call reada(controlcard,"ELPP3_UP",
1148      &          elpp3_up(itype1,itype2),0.0d0)
1149               elpp3_up(itype2,itype1)=elpp3_up(itype1,itype2)
1150             endif
1151             call card_concat(controlcard,.true.)
1152           enddo
1153           call card_concat(controlcard,.true.)
1154         endif
1155         do i=1,20
1156           do j=1,2
1157             do k=1,2
1158               mask_scp(i,j,k)=0
1159             enddo
1160           enddo
1161         enddo
1162         mod_scp=index(controlcard,"MOD_SCP").gt.0
1163         if (mod_scp) then
1164           mod_scp=.false.
1165           call card_concat(controlcard,.true.)
1166           do while (index(controlcard,"END").eq.0)
1167             if (index(controlcard,"EPSCP").gt.0) then
1168               mod_scp=.true.
1169               call readi(controlcard,"ITYPSC",itypsc,0)
1170               call readi(controlcard,"ITYPP",itypp,0)
1171               if (itypsc.gt.20 .or. itypp.eq.0 .or. itypp.gt.2) then
1172                 write (iout,*) "Error specifying scp parms"
1173                 stop
1174               endif
1175               mask_scp(itypsc,itypp,1)=1
1176               call reada(controlcard,"EPSCP_LOW",
1177      &          epscp_low(itypsc,itypp),0.1d0)
1178               call reada(controlcard,"EPSCP_UP",
1179      &          epscp_up(itypsc,itypp),0.0d0)
1180             endif
1181             if (index(controlcard,"RSCP").gt.0) then
1182               mod_scp=.true.
1183               call readi(controlcard,"ITYPSC",itypsc,0)
1184               call readi(controlcard,"ITYPP",itypp,0)
1185               mask_scp(itypsc,itypp,2)=1
1186               call readi(controlcard,"ITYPSC",itypsc,0)
1187               call readi(controlcard,"ITYPP",itypp,0)
1188               if (itypsc.gt.20 .or. itypp.eq.0 .or. itypp.gt.2) then
1189                 write (iout,*) "Error specifying scp parms"
1190                 stop
1191               endif
1192               call reada(controlcard,"RSCP_LOW",
1193      &          rscp_low(itypsc,itypp),0.1d0)
1194 c              write(iout,*)itypsc,itypp,rscp_low(itypsc,itypp)
1195               call reada(controlcard,"RSCP_UP",
1196      &          rscp_up(itypsc,itypp),0.0d0)
1197 c              write(iout,*)itypsc,itypp,rscp_up(itypsc,itypp)
1198             endif
1199             call card_concat(controlcard,.true.)
1200           enddo
1201         endif
1202         write (iout,*) "END ",controlcard
1203         call flush(iout)
1204         write (iout,*) "mod_fourier",mod_fourier(nloctyp),
1205      &   " mod_elec ",mod_elec," mod_scp ",mod_scp,
1206      &   " mod_sccor",mod_sccor," mod_ang ",mod_ang,
1207      &   " ind_shield", imask(ind_shield)
1208         mod_other_params=
1209      &   mod_fourier(nloctyp).gt.0 .or. mod_elec .or. mod_scp
1210      &     .or. mod_sccor .or. mod_ang .or. mod_side_other
1211      &     .or. imask(ind_shield).eq.1
1212         if (read_stat.lt.2. .and. mod_other_params) then
1213           write (iout,*)"Error - only weights and sidechain parameters",
1214      &     " can be optimized with READ_STAT=",read_stat
1215           call flush(iout)
1216           return1
1217         endif
1218         init_ene = init_ene .or. read_stat.eq.2
1219       endif
1220       write (iout,*) "End read: MOD_OTHER_PARAMS ",mod_other_params
1221 c      write (*,*) "MOD_SIDE ",mod_side," MOD_FOURIER",
1222 c     & mod_fourier(nloctyp),
1223 c     & " MOD_ELEC ",mod_elec," MOD_SCP ",mod_scp
1224       init_ene=init_ene .or. mod_other_params
1225 c      write (iout,*) "init_ene",init_ene
1226 c      write (iout,*)
1227       call flush(iout)
1228       return
1229       end
1230 c-----------------------------------------------------------------------------
1231       subroutine print_general_data(*)
1232       implicit none
1233       include "DIMENSIONS"
1234       include "DIMENSIONS.ZSCOPT"
1235 #ifdef MPI
1236       include "mpif.h"
1237       include "COMMON.MPI"
1238       integer ierror,kolor,klucz
1239 #endif
1240       include "COMMON.WEIGHTS"
1241       include "COMMON.NAMES"
1242       include "COMMON.VMCPAR"
1243       include "COMMON.TORSION"
1244       include "COMMON.LOCAL"
1245       include "COMMON.INTERACT"
1246       include "COMMON.ENERGIES"
1247       include "COMMON.MINPAR"
1248       include "COMMON.IOUNITS"
1249       include "COMMON.TIME1"
1250       include "COMMON.PROTFILES"
1251       include "COMMON.CHAIN"
1252       include "COMMON.CLASSES"
1253       include "COMMON.THERMAL"
1254       include "COMMON.FFIELD"
1255       include "COMMON.OPTIM"
1256       include "COMMON.CONTROL"
1257       include "COMMON.SCCOR"
1258       character*800 controlcard
1259       integer i,j,k,l,ii,n_ene_found,itypt,jtypt
1260       integer ind,itype1,itype2,itypf,itypsc,itypp
1261       integer ilen,lenpot,lenpre
1262       external ilen
1263       character*4 liczba,liczba1
1264       if (mode.eq.1) then
1265         write (iout,*) "Single-point target function evaluation"
1266       else if (mode.eq.2) then
1267         write (iout,*) "Grid search of the parameter space"
1268       else if (mode.eq.3) then
1269         write (iout,*) "Minimization of the target function"
1270       else
1271         write (iout,*) "Wrong MODE",mode,", should be 1, 2, or 3"
1272         call flush(iout)
1273         return1
1274       endif
1275       write (iout,*) "RESCALE_MODE",rescale_mode
1276 c      mod_other_params=index(controlcard,"OPTIMIZE_OTHER").gt.0
1277 c      if (read_stat.eq.0 .and. mod_other_params) then
1278 c        write (iout,*) "Error: only optimization of energy-term",
1279 c     &    " weights can be performed with READ_STAT=",read_stat
1280 c        call flush(iout)
1281 c        return1
1282 c      endif
1283       if (mode.eq.3) then
1284         write (iout,*) "Parameters of the SUMSL procedure:"
1285         write (iout,'(a,t15,i5)') "MAXMIN",maxmin
1286         write (iout,'(a,t15,i5)') "MAXFUN",maxfun
1287         write (iout,'(a,t15,e15.7)') "TOLF",tolf
1288         write (iout,'(a,t15,e15.7)') "RTOLF",rtolf
1289       endif
1290       if (mod_other_params) then
1291         write (iout,*)
1292      &"Internal parameters of UNRES energy components will be optimized"
1293 c        call card_concat(controlcard,.true.)
1294 c        mod_side = (index(controlcard,"MOD_SIDE").gt.0)
1295         if (mod_side) then
1296           write (iout,*) "SC epsilons to be optimized:"
1297           write (iout,*) "Single: eps(i,j)=eps(i,j)+(x(i)+x(j))/2"
1298           do i=1,nsingle_sc
1299             write (iout,*) restyp(ityp_ssc(i)),epss_low(i),epss_up(i)
1300           enddo
1301           write (iout,*) "Pairs: eps(i,j)=eps(i,j)+x(i,j):"
1302           do i=1,npair_sc
1303             write (iout,*) restyp(ityp_psc(1,i)),
1304      &       restyp(ityp_psc(2,i)),epsp_low(i),epsp_up(i)
1305           enddo
1306         endif
1307         if (mod_sccor) then
1308           write (iout,*)"Torsional weights/coefficients to be optimized"
1309           write(iout,'(a)') "TYP1 TYP2    WEIGHT",
1310      &     "   LOWER BOUND UPPER BOUND"
1311           do itypt=-nsccortyp,nsccortyp
1312             do jtypt=-nsccortyp,nsccortyp
1313               do l=1,3
1314               if (mask_tor(l,itypt,jtypt,1).gt.0) then
1315                 write(iout,'(3a4,3f10.5)')l,restyp(itypt),restyp(jtypt),
1316      &           weitor(l,itypt,jtypt,1),weitor_low(l,itypt,jtypt,1),
1317      &           weitor_up(l,itypt,jtypt,1)
1318               endif
1319               enddo
1320             enddo
1321           enddo
1322         endif
1323 c 7/8/17 AL: Optimization of the bending parameters
1324         if (mod_ang) then
1325           write (iout,*) "Boundaries of angle potential coefficients"
1326           write (iout,'(3a)') "Index","     Low bound","       Up bound"
1327           do i=0,nthetyp
1328             if (mask_ang(i).eq.0) cycle
1329             write (iout,'(2a)') 'Type ',restyp(iloctyp(i))
1330             do j=1,nbend_kcc_TB(i)
1331             write (iout,'(i5,2f15.3)') j,v1bend_low(j,i),v1bend_up(j,i)
1332             enddo
1333           enddo
1334         endif
1335         if (mod_fourier(nloctyp).gt.0) then
1336           write (iout,*) "Fourier coefficients to be optimized"
1337           do itypf=0,nloctyp-1
1338             if (mod_fourier(itypf).gt.0) then
1339               write (iout,'(3a,i2)') "Type ",restyp(iloctyp(itypf)),
1340      &           " number of coeffs",mod_fourier(itypf)
1341               write(iout,'(a)') "NAME  COEFF    LOWER BOUND UPPER BOUND"
1342 #ifdef NEWCORR
1343              do j=1,2
1344                do k=1,3
1345                  if (mask_bnew1(k,j,itypf).gt.0)
1346      &            write (iout,'(2hB1,2i2,f10.5)') k,j,bnew1(k,j,itypf),
1347      &             bnew1_low(k,j,itypf),bnew1_up(k,j,itypf)
1348                enddo
1349              enddo
1350              do j=1,2
1351                do k=1,3
1352                  if (mask_bnew2(k,j,itypf).gt.0)
1353      &            write (iout,'(2hB2,2i2,3f11.5)') k,j,bnew2(k,j,itypf),
1354      &             bnew2_low(k,j,itypf),bnew2_up(k,j,itypf)
1355                enddo
1356              enddo
1357              do j=1,2
1358                do k=1,3
1359                  if (mask_ccnew(k,j,itypf).gt.0)
1360      &            write (iout,'(2hCC,2i2,3f11.5)') k,j,ccnew(k,j,itypf),
1361      &             ccnew_low(k,j,itypf),ccnew_up(k,j,itypf)
1362                enddo
1363              enddo
1364              do j=1,2
1365                do k=1,3
1366                  if (mask_ddnew(k,j,itypf).gt.0)
1367      &            write (iout,'(2hDD,2i2,3f11.5)') k,j,ddnew(k,j,itypf),
1368      &             ddnew_low(k,j,itypf),ddnew_up(k,j,itypf)
1369                enddo
1370              enddo
1371              do j=1,2
1372                if (mask_e0new(j,itypf).gt.0)
1373      &            write (iout,'(2hE0,i2,3f11.5)') j,e0new(j,itypf),
1374      &             e0new_low(j,itypf),e0new_up(j,itypf)
1375              enddo
1376              do j=1,2
1377                do k=1,2
1378                  do l=1,2
1379                    if (mask_eenew(l,k,j,itypf).gt.0)
1380      &              write (iout,'(2hEE,3i2,3f11.5)')
1381      &               l,k,j,eenew(l,k,j,itypf),eenew_low(l,k,j,itypf),
1382      &               eenew_up(l,k,j,itypf)
1383                  enddo
1384                enddo
1385              enddo
1386 #else
1387               do i=2,13
1388                 if (mask_fourier(i,itypf).gt.0) then
1389                   write (iout,'(1hB,i2,3f11.5)')
1390      &            i,b(i,itypf),b_low(i,itypf),b_up(i,itypf)
1391                 endif
1392               enddo
1393 #endif
1394             endif
1395           enddo
1396         endif
1397         if (mod_elec) then
1398           write (iout,*)
1399           write (iout,*) "Electrostatic parameters to be optimized"
1400           do itype1=1,2
1401             do itype2=1,itype1
1402               if (mask_elec(itype1,itype2,1).gt.0)
1403      &        write (iout,'(2i3," EPP",f8.3," EPP_LOW",f8.3,
1404      &          " EPP_UP",f8.3)')
1405      &          itype1,itype2,epp(itype1,itype2),epp_low(itype1,itype2),
1406      &          epp_up(itype1,itype2)
1407               if (mask_elec(itype1,itype2,2).gt.0)
1408      &        write (iout,'(2i3," RPP",f8.3," RPP_LOW",f8.3,
1409      &          " RPP_UP",f8.3)')
1410      &          itype1,itype2,rpp(itype1,itype2),rpp_low(itype1,itype2),
1411      &          rpp_up(itype1,itype2)
1412               if (mask_elec(itype1,itype2,3).gt.0)
1413      &        write (iout,'(2i3," ELPP6",f8.3," ELPP6_LOW",f8.3,
1414      &          " ELPP6_UP",f8.3)')
1415      &          itype1,itype2,elpp6(itype1,itype2),
1416      &          elpp6_low(itype1,itype2),elpp6_up(itype1,itype2)
1417               if (mask_elec(itype1,itype2,4).gt.0)
1418      &        write (iout,'(2i3," ELPP3",f8.3," ELPP3_LOW",f8.3,
1419      &          " ELPP3_UP",f8.3)')
1420      &          itype1,itype2,elpp3(itype1,itype2),
1421      &          elpp3_low(itype1,itype2),elpp3_up(itype1,itype2)
1422             enddo
1423           enddo
1424         endif
1425         if (mod_scp) then
1426           mod_scp=.false.
1427           write (iout,*)
1428           write (iout,*) "SCP parameters to be optimized:"
1429           if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
1430      &        mask_scp(0,2,2) .gt. 0) then
1431             write (iout,*)
1432      & "SCP parameters are assumed to depend on peptide group type only"
1433             do j=1,2
1434               if (mask_scp(0,j,1).gt.0)
1435      &          write (iout,'(i3," EPSCP",f8.3," EPSCP_LOW",f8.3,
1436      &           " EPSCP_UP",f8.3)') j,eps_scp(1,j),epscp_low(0,j),
1437      &           epscp_up(0,j)
1438               if (mask_scp(0,j,2).gt.0)
1439      &          write (iout,'(i3," RSCP",f8.3," RSCP_LOW",f8.3,
1440      &           " RSCP_UP",f8.3)') j,rscp(1,j),rscp_low(0,j),
1441      &           rscp_up(0,j)
1442             enddo
1443           else
1444             do i=1,20
1445               do j=1,2
1446                 if (mask_scp(i,j,1).gt.0)
1447      &            write (iout,'(2i3," EPSCP",f8.3," EPSCP_LOW",f8.3,
1448      &             " EPSCP_UP",f8.3)') i,j,eps_scp(i,j),epscp_low(i,j),
1449      &             epscp_up(i,j)
1450                 if (mask_scp(i,j,2).gt.0)
1451      &            write (iout,'(2i3," RSCP",f8.3," RSCP_LOW",f8.3,
1452      &             " RSCP_UP",f8.3)') i,j,rscp(i,j),rscp_low(i,j),
1453      &             rscp_up(i,j)
1454               enddo
1455             enddo
1456           endif
1457         endif
1458       endif
1459       write (iout,*) "MOD_OTHER_PARAMS ",mod_other_params
1460       write (iout,*) "MOD_SIDE ",mod_side," MOD_FOURIER",
1461      & mod_fourier(nloctyp),
1462      & " MOD_ELEC ",mod_elec," MOD_SCP ",mod_scp," mod_ang",mod_ang
1463       init_ene=init_ene .or. mod_other_params
1464       write (iout,*) "init_ene",init_ene
1465       write (iout,*)
1466       call flush(iout)
1467       return
1468       end
1469 c-----------------------------------------------------------------------------
1470       subroutine read_protein_data(*)
1471       implicit none
1472       include "DIMENSIONS"
1473       include "DIMENSIONS.ZSCOPT"
1474 #ifdef MPI
1475       include "mpif.h"
1476       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1477       include "COMMON.MPI"
1478 #endif
1479       include "COMMON.CONTROL"
1480       include "COMMON.CHAIN"
1481       include "COMMON.CLASSES"
1482       include "COMMON.COMPAR"
1483       include "COMMON.ENERGIES"
1484       include "COMMON.IOUNITS"
1485       include "COMMON.PROTFILES"
1486       include "COMMON.PROTNAME"
1487       include "COMMON.VMCPAR"
1488       include "COMMON.OPTIM"
1489       include "COMMON.WEIGHTS"
1490       include "COMMON.NAMES"
1491       include "COMMON.ALLPROT"
1492       include "COMMON.THERMAL"
1493       include "COMMON.TIME1"
1494       include "COMMON.WEIGHTDER"
1495       character*64 nazwa,key
1496       character*16000 controlcard
1497       character*16000 all_protfiles
1498       integer i,j,k,kk,l,iprot,ii,if,ib,ibatch,nn,nn1,iww,maskcheck,
1499      &  il,inat,ilevel,weightread,jfrag,jclass,nfragm,iparm
1500       integer nrec,nlines,iscor
1501       double precision energ,wangnorm(maxT),wq(maxT),wrms(maxT),
1502      &  wrgy(maxT),wsign(maxT),wangnat(maxT),wqnat(maxT),wrmsnat(maxT),
1503      &  wrgynat(maxT),wcorangnorm(maxT),wcorq(maxT),
1504      &  wcorrms(maxT),wcorrgy(maxT),wcorsign(maxT),wcorangnat(maxT),
1505      &  wcorqnat(maxT),wcorrmsnat(maxT),wcorrgynat(maxT),
1506      &  angnormlow(maxT),qlow(maxT),rmslow(maxT),
1507      &  rgylow(maxT),signlow(maxT),angnatlow(maxT),
1508      &  qnatlow(maxT),rmsnatlow(maxT),rgynatlow(maxT),
1509      &  angcorlow(maxT),qcorlow(maxT),rmscorlow(maxT),rgycorlow(maxT),
1510      &  signcorlow(maxT),angcornatlow(maxT),qcornatlow(maxT),
1511      &  rmscornatlow(maxT),rgycornatlow(maxT),signcornatlow(maxT),
1512      &  angnormup(maxT),qup(maxT),rmsup(maxT),rgyup(maxT),signup(maxT),
1513      &  angnatup(maxT),qnatup(maxT),rmsnatup(maxT),rgynatup(maxT),
1514      &  angcorup(maxT),qcorup(maxT),rmscorup(maxT),rgycorup(maxT),
1515      &  signcorup(maxT),
1516      &  angcornatup(maxT),qcornatup(maxT),rmscornatup(maxT),
1517      &  rgycornatup(maxT),signcornatup(MaxT)
1518       integer ilen,iroof
1519       external ilen,iroof
1520       double precision ebia(maxprot),rjunk
1521       character*10 liczba
1522       character*64 zeros /
1523      &'0000000000000000000000000000000000000000000000000000000000000000'
1524      & /
1525       logical lerr
1526 c      print *,"Processor",me," calls read_protein_data"
1527
1528 C Read seventh record: general data of proteins used for calibration
1529       call card_concat(controlcard,.true.)
1530 c      write(2, *)controlcard
1531       call readi(controlcard,"NPROT",nprot,0)
1532       pdbref=index(controlcard,"PDBREF").gt.0
1533       print_refstr=index(controlcard,"PRINT_REFSTR").gt.0
1534       if (nprot.eq.0) then
1535         write(iout,*) "Error: at least one protein must be specified."
1536         call flush(iout)
1537         return1
1538       endif
1539       protname(0)="all"
1540       do iprot=1,nprot
1541         read (inp,'(a)') protname(iprot)
1542         write (iout,*)
1543         write (iout,*) "Reading data of protein",iprot," named ",
1544      &   protname(iprot)
1545         call card_concat(controlcard,.true.)
1546         maxlik(iprot)=index(controlcard,"MAXLIK").gt.0
1547         fmatch(iprot)=index(controlcard,"FMATCH").gt.0
1548         call readi(controlcard,"ISTART_MD",istart_MD(iprot),1)
1549         call readi(controlcard,"IEND_MD",iend_MD(iprot),maxstr)
1550         call readi(controlcard,"IFREQ_MD",ifreq_MD(iprot),1)
1551         call reada(controlcard,"WFORCE",wforce(iprot),1.0d0)
1552         mdbox(iprot)=index(controlcard,"MDBOX").gt.0
1553         write (iout,*) "Protein: ",
1554      &   protname(iprot)(:ilen(protname(iprot)))," maxlik ",
1555      &   maxlik(iprot)," fmatch",fmatch(iprot)," ifreq_MD:",
1556      &   ifreq_MD(iprot)," istart_MD",istart_MD(iprot)," iend_MD",
1557      &   iend_MD(iprot)," wforce",wforce(iprot),ifreq_MD(iprot)
1558         call reada(controlcard,"ENECUT_MIN",enecut_min(iprot),15.0d0)
1559         call reada(controlcard,"ENECUT_MAX",enecut_max(iprot),100.0d0)
1560         call reada(controlcard,"ENECUT",enecut(iprot),enecut_min(iprot))
1561         if (enecut(iprot).lt.enecut_min(iprot))
1562      &    enecut(iprot)=enecut_min(iprot)
1563         if (enecut_max(iprot).le.enecut_min(iprot))
1564      &    enecut_max(iprot)=2*enecut_min(iprot)
1565         write (iout,'(3(a,f9.1))') "ENECUT",enecut(iprot)," ENECUT_MIN",
1566      &   enecut_min(iprot)," ENECUT_MAX",enecut_max(iprot)
1567         call readi(controlcard,"HEFAC",hefac(iprot),50)
1568         call readi(controlcard,"HTFAC",htfac(iprot),50)
1569         call readi(controlcard,"HELOW",hemax_low(iprot),20)
1570         call readi(controlcard,"HTLOW",htmax_low(iprot),20)
1571         write (iout,*) "iprot",iprot,
1572      &   " hefac",hefac(iprot)," helow",hemax_low(iprot),
1573      &   " htfac",htfac(iprot)," htlow",htmax_low(iprot)
1574 c 7/27/2013 AL Read maximum likelihood data
1575
1576         if (maxlik(iprot)) then
1577
1578         call card_concat(controlcard,.true.)
1579         call readi(controlcard,"NBETA_L",nbeta(1,iprot),0)
1580         write (iout,*) "NBETA_L",nbeta(1,iprot)
1581         caonly(iprot)=index(controlcard,"CAONLY").gt.0
1582         sconly(iprot)=index(controlcard,"SCONLY").gt.0
1583         rmscomp(iprot)=index(controlcard,"RMSCOMP").gt.0
1584         anglecomp(iprot)=index(controlcard,"ANGLECOMP").gt.0
1585         sidecomp(iprot)=index(controlcard,"SIDECOMP").gt.0
1586         call reada(controlcard,"SIGMA",sigma2(iprot),4.0d0)
1587         call reada(controlcard,"SIGMAANG",sigmaang2(iprot),4.0d0)
1588         call reada(controlcard,"SIGMASIDE",sigmaside2(iprot),4.0d0)
1589         write (iout,*) "RMSCOMP",rmscomp(iprot)," SIGMA",sigma2(iprot),
1590      &     " CAONLY ",caonly(iprot)," SCONLY",sconly(iprot)
1591         write (iout,*) "ANGLECOMP",anglecomp(iprot),
1592      &     " SIGMAANG",sigmaang2(iprot)
1593         write (iout,*) "SIDECOMP",sidecomp(iprot),
1594      &     " SIGMASIDE",sigmaside2(iprot)
1595         do ib=1,nbeta(1,iprot)
1596           read(inp,*)betaT(ib,1,iprot),weilik(ib,iprot),
1597      &      pdbfile(ib,iprot)
1598           write (iout,*) i,betaT(ib,1,iprot),weilik(ib,iprot),
1599      &      pdbfile(ib,iprot)
1600         enddo
1601 c 7/27/2013 AL Read heat-capacity data
1602         call card_concat(controlcard,.true.)
1603         call readi(controlcard,"NBETA_CV",nbeta(2,iprot),0)
1604         call reada(controlcard,"WCV",wcv(iprot),1.0d0)
1605         call reada(controlcard,"BASE",heatbase(iprot),0.0d0)
1606         write (iout,*) "NBETA_CV",nbeta(2,iprot)," WCV",wcv(iprot)
1607         do ib=1,nbeta(2,iprot)
1608           read(inp,*) betaT(ib,2,iprot),target_cv(ib,iprot),
1609      &      weiCv(ib,iprot)
1610           weiCv(ib,iprot)=weiCv(ib,iprot)*wcv(iprot)
1611           write (iout,*) betaT(ib,2,iprot),target_cv(ib,iprot),
1612      &      weiCv(ib,iprot)
1613         enddo
1614         write (iout,*) "Experimental heat capacity curve"
1615         do ib=1,nbeta(2,iprot)
1616           write (iout,'(f6.2,2f10.5)') betaT(ib,2,iprot),
1617      &     target_cv(ib,iprot),weiCv(ib,iprot)
1618         enddo
1619         write (iout,'(a,f10.5)') "Baseline",heatbase(iprot)
1620         call flush(iout)
1621 c Conformation-dependent averages
1622         call card_concat(controlcard,.true.)
1623         call readi(controlcard,"NATLIKE",natlike(iprot),0)
1624         do i=1,natlike(iprot)
1625           call card_concat(controlcard,.true.)
1626           call reada(controlcard,"WNAT",wnat(i,iprot),1.0d0)
1627           call readi(controlcard,"NUMNAT",numnat(i,iprot),1)
1628           call readi(controlcard,"NATDIM",natdim(i,iprot),1)
1629           do ib=1,nbeta(i+2,iprot)
1630             read(inp,*)betaT(ib,i+2,iprot),(weinat(k,ib,i,iprot),
1631      &             nuexp(k,ib,i,iprot),k=1,natdim(i,iprot))
1632           enddo
1633         enddo
1634         do i=1,natlike(iprot)+2
1635           do j=1,nbeta(i,iprot)
1636             betaT(j,i,iprot)=1.0d0/(Rgas*betaT(j,i,iprot))
1637             write (2,*) "R i",i," j",j," beta",betaT(j,i,iprot)
1638           enddo
1639         enddo
1640         do iparm=1,nparmset
1641
1642 C Read names of files with the data for protein IPROT
1643         call card_concat(controlcard,.false.)
1644
1645         if (iparm.eq.myparm) then
1646         call split_string(controlcard,protfiles(1,iprot),
1647      &     maxfile_prot,nfile_prot(iprot))
1648 #ifdef DEBUG
1649         write(iout,*)"iprot",iprot," nfile",nfile_prot(iprot)
1650         write(iout,*)
1651      &     (protfiles(i,iprot),i=1,nfile_prot(iprot))
1652 #endif
1653         endif
1654
1655         enddo ! iparm
1656         endif ! maxlik(iprot)
1657
1658 c Read molecular information of the protein
1659         call molread_zs(iprot)
1660 c 3/31/04 AL Read the reference structures of protein iprot
1661 c        print *,"Calling read_ref_structure"
1662         if (maxlik(iprot)) call read_ref_structure(iprot,*11)
1663
1664         if (fmatch(iprot)) then
1665           call card_concat(controlcard,.false.)
1666           call reads(controlcard,"atomdef",fpdbfile(iprot),"inpcrd.pdb")
1667           call reads(controlcard,"mdcrd",fcoordfile(iprot),
1668      &      "equil2.mdpdb")
1669           call reads(controlcard,"mdforce",fforcefile(iprot),
1670      &      "equil2.mdfrc")
1671           write (iout,*) "protein:    ",
1672      &      protname(iprot)(:ilen(protname(iprot)))
1673           write (iout,*) "fpdbfile:   ",
1674      &      fpdbfile(iprot)(:ilen(fpdbfile(iprot)))
1675           write (iout,*) "fcoordfile: ",
1676      &      fcoordfile(iprot)(:ilen(fcoordfile(iprot)))
1677           write(iout,*)  "fforcefile: ",
1678      &      fforcefile(iprot)(:ilen(fforcefile(iprot)))
1679           write (iout,*) "Calling atom_partition"
1680           call atom_partition(iprot)
1681           write (iout,*) "After atom_partition"
1682 c          call read_allat_coord_forces(iprot)
1683 c          write (iout,*) "After forces"
1684         endif
1685 #ifdef DEBUG
1686         write (iout,*) "Protein",iprot," left READ_REF_STRUCTURE"
1687 #endif
1688       enddo ! iprot
1689       return
1690    11 return1
1691       end
1692 c-------------------------------------------------------------------------------
1693       subroutine read_database(*)
1694       implicit none
1695       include "DIMENSIONS"
1696       include "DIMENSIONS.ZSCOPT"
1697 #ifdef MPI
1698       include "mpif.h"
1699       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1700       include "COMMON.MPI"
1701 #endif
1702       include "COMMON.CHAIN"
1703       include "COMMON.INTERACT"
1704       include "COMMON.CLASSES"
1705       include "COMMON.ENERGIES"
1706       include "COMMON.IOUNITS"
1707       include "COMMON.PROTFILES"
1708       include "COMMON.PROTNAME"
1709       include "COMMON.VMCPAR"
1710       include "COMMON.NAMES"
1711       include "COMMON.ALLPROT"
1712       include "COMMON.WEIGHTS"
1713       include "COMMON.WEIGHTDER"
1714       include "COMMON.VAR"
1715       include "COMMON.SBRIDGE"
1716       include "COMMON.GEO"
1717       include "COMMON.COMPAR"
1718       include "COMMON.OPTIM"
1719       include "COMMON.FMATCH"
1720       character*128 nazwa
1721       character*16000 controlcard
1722       character*16000 all_protfiles
1723       character*4 liczba,liczba1
1724       integer i,j,jj,jjj,jj_old,icount,k,kk,l,iprot,ii,if,ib,ibatch,
1725      &  nn,nn1,inan
1726       integer ixdrf,iret,itmp
1727       integer nrec,nlines,iscor
1728       double precision energ,t_acq,tcpu
1729       integer ilen,iroof
1730       external ilen,iroof
1731       double precision rjunk
1732       integer ntot_all(0:maxprot,0:maxprocs-1)
1733       logical lerr
1734       double precision energia(0:max_ene),etot
1735       real*4 csingle(3,maxres2),reini,refree,rmsdev,prec
1736       integer Previous,Next
1737 c      print *,"Processor",me," calls read_protein_data"
1738 #ifdef MPI
1739       if (me.eq.master) then
1740         Previous=MPI_PROC_NULL
1741       else
1742         Previous=me-1
1743       endif
1744       if (me.eq.nprocs-1) then
1745         Next=MPI_PROC_NULL
1746       else
1747         Next=me+1
1748       endif
1749 #endif
1750 c Set the scratchfile names
1751       write (liczba,'(bz,i4.4)') me
1752       write (liczba1,'(bz,i4.4)') myparm
1753       do iprot=1,nprot
1754         write (iout,*) "READ_DATABASE: iprot",iprot,
1755      &   " maxlik ",maxlik(iprot)," fmatch ",fmatch(iprot)
1756         call flush(iout)
1757         if (maxlik(iprot)) then
1758 c 1/27/05 AL Change stored coordinates to single precision and don't store
1759 c         energy components in the binary databases.
1760         lenrec(iprot)=12*(nres_zs(iprot)+nct_zs(iprot)-nnt_zs(iprot)+1)
1761      &   +4*(2*nss_zs(1,iprot)+1)+8*natlike(iprot)*maxdimnat+16
1762 c 4/13/04 AL Add space for similarity measure
1763         lenrec_ene(iprot) = (2*nntyp+3*n_ene+2)*8
1764      &   +8*natlike(iprot)*maxdimnat
1765 #ifdef DEBUG
1766         call flush(iout)
1767         write (iout,*) "Protein i"," lenrec",lenrec(iprot)
1768         write (iout,*) "lenrec_ene",lenrec_ene(iprot)
1769         call flush(iout)
1770 #endif
1771         bprotfiles(iprot) = scratchdir(:ilen(scratchdir))//
1772      &       "/"//protname(iprot)(:ilen(protname(iprot)))//
1773      &       "_par"//liczba1//"_"//liczba//".xbin"
1774         benefiles(iprot) = scratchdir(:ilen(scratchdir))//
1775      &       "/"//protname(iprot)(:ilen(protname(iprot)))//
1776      &       "_par"//liczba1//"_"//liczba//".enebin"
1777 c        write (iout,*) "scratchfile ",
1778 c     &    bprotfiles(iprot)(:ilen(bprotfiles(iprot)))
1779         endif
1780         if (fmatch(iprot)) then
1781 c 1/27/05 AL Change stored coordinates to single precision and don't store
1782 c         energy components in the binary databases.
1783         write (iout,*) "nres_zs",nres_zs(iprot),"nnt_zs",nnt_zs(iprot),
1784      &    " nct_zs",nct_zs(iprot)," nsite",nsite(iprot)
1785         lenrec_MD(iprot)=12*(nres_zs(iprot)+nct_zs(iprot)-nnt_zs(iprot)
1786      &   +1+nsite(iprot))
1787      &   +4*(ns_zs(iprot)+1)
1788 c 4/13/04 AL Add space for similarity measure
1789         lenrec_forces(iprot) = maxres6*n_ene*8
1790 #ifdef DEBUG
1791         write (iout,*)
1792      &   "Protein ",protname(iprot)(:ilen(protname(iprot))),
1793      &   " lenrec_MD",lenrec_MD(iprot),
1794      &   " lenrec_forces",lenrec_forces(iprot)
1795         call flush(iout)
1796 #endif
1797         bprotfiles_MD(iprot) = scratchdir(:ilen(scratchdir))//
1798      &       "/"//protname(iprot)(:ilen(protname(iprot)))//
1799      &       "_par"//liczba1//"_"//liczba//".fxbin"
1800         bforcefiles(iprot) = scratchdir(:ilen(scratchdir))//
1801      &       "/"//protname(iprot)(:ilen(protname(iprot)))//
1802      &       "_par"//liczba1//"_"//liczba//".fcompbin"
1803         write (iout,*) "bprotfiles_MD ",
1804      &    bprotfiles_MD(iprot)(:ilen(bprotfiles_MD(iprot)))
1805         write (iout,*) "bforcefiles ",
1806      &    bforcefiles(iprot)(:ilen(bforcefiles(iprot)))
1807         call flush(iout)
1808         endif
1809       enddo
1810
1811       do i=0,nprot
1812         ntot(i)=0
1813         ntot_MD(i)=0
1814       enddo
1815       do iprot=1,nprot
1816         if (maxlik(iprot)) call read_MREMD_coord(iprot,Previous,Next)
1817         if (fmatch(iprot)) call read_allat_coord_forces(iprot)
1818       enddo ! iprot
1819       write(iout,*)"A total of",ntot(0)," MREMD conformations read."
1820       write(iout,*)"A total of",ntot_MD(0)," MD conformations read."
1821 #ifdef MPI
1822 c Check if everyone has the same number of MREMD conformations
1823       call MPI_Allgather(ntot(0),maxprot+1,MPI_INTEGER,
1824      &  ntot_all(0,0),maxprot+1,MPI_INTEGER,MPI_Comm_World,IERROR)
1825       lerr=.false.
1826       do i=0,nprocs-1
1827         if (i.ne.me) then
1828           do j=1,nprot
1829             if (.not.maxlik(iprot)) cycle
1830             if (ntot(j).ne.ntot_all(j,i)) then
1831               write (iout,*)
1832      &         "Number of MREMD conformations at processor",i,
1833      &         " for protein",j," differs from that at processor",me,
1834      &         ntot(j),ntot_all(j,i)
1835               lerr = .true.
1836             endif
1837           enddo
1838         endif
1839       enddo
1840       if (lerr) then
1841         write (iout,*)
1842         write (iout,*)"Number of MREMD conformations read by processors"
1843         write (iout,'(10x,7a10)') (protname(i),i=0,nprot)
1844         write (iout,*)
1845         do i=0,nprocs-1
1846           write (iout,'(8i10)') i,(ntot_all(j,i),j=0,nprot)
1847         enddo
1848         write (iout,*) "Calculation terminated."
1849         call flush(iout)
1850         return1
1851       endif
1852       return
1853 c Check if everyone has the same number of MD conformations
1854       call MPI_Allgather(ntot_MD(0),maxprot+1,MPI_INTEGER,
1855      &  ntot_all(0,0),maxprot+1,MPI_INTEGER,MPI_Comm_World,IERROR)
1856       lerr=.false.
1857       do i=0,nprocs-1
1858         if (i.ne.me) then
1859           do j=1,nprot
1860             if (.not.fmatch(iprot)) cycle
1861             if (ntot(j).ne.ntot_all(j,i)) then
1862               write (iout,*)"Number of MD conformations at processor",i,
1863      &         " for protein",j," differs from that at processor",me,
1864      &         ntot_MD(j),ntot_all(j,i)
1865               lerr = .true.
1866             endif
1867           enddo
1868         endif
1869       enddo
1870       if (lerr) then
1871         write (iout,*)
1872         write (iout,*) "Number of MD conformations read by processors"
1873         write (iout,'(10x,7a10)') (protname(i),i=0,nprot)
1874         write (iout,*)
1875         do i=0,nprocs-1
1876           write (iout,'(8i10)') i,(ntot_all(j,i),j=0,nprot)
1877         enddo
1878         write (iout,*) "Calculation terminated."
1879         call flush(iout)
1880         return1
1881       endif
1882       return
1883 #endif
1884  1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
1885       call flush(iout)
1886 #ifdef MPI
1887       call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1888 #endif
1889       return1
1890       end
1891 c------------------------------------------------------------------------
1892       subroutine read_MREMD_coord(iprot,previous,next)
1893       implicit none
1894       include "DIMENSIONS"
1895       include "DIMENSIONS.ZSCOPT"
1896 #ifdef MPI
1897       include "mpif.h"
1898       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1899       include "COMMON.MPI"
1900 #endif
1901       include "COMMON.CHAIN"
1902       include "COMMON.INTERACT"
1903       include "COMMON.CLASSES"
1904       include "COMMON.ENERGIES"
1905       include "COMMON.IOUNITS"
1906       include "COMMON.PROTFILES"
1907       include "COMMON.PROTNAME"
1908       include "COMMON.VMCPAR"
1909       include "COMMON.NAMES"
1910       include "COMMON.ALLPROT"
1911       include "COMMON.WEIGHTS"
1912       include "COMMON.WEIGHTDER"
1913       include "COMMON.VAR"
1914       include "COMMON.SBRIDGE"
1915       include "COMMON.GEO"
1916       include "COMMON.COMPAR"
1917       include "COMMON.OPTIM"
1918       character*128 nazwa
1919       character*16000 controlcard
1920       character*16000 all_protfiles
1921       character*4 liczba,liczba1
1922       integer i,j,jj,jjj,jj_old,icount,k,kk,l,iprot,ii,if,ib,ibatch,
1923      &  nn,nn1,inan
1924       integer ixdrf,iret,itmp
1925       integer nrec,nlines,iscor
1926       double precision energ,t_acq,tcpu
1927       integer ilen,iroof
1928       external ilen,iroof
1929       double precision rjunk
1930       integer ntot_all(0:maxprot,0:maxprocs-1)
1931       logical lerr
1932       double precision energia(0:max_ene),etot
1933       real*4 csingle(3,maxres2),reini,refree,rmsdev,prec
1934       integer Previous,Next
1935         jj_old=1
1936         call restore_molinfo(iprot)
1937         open (ientout,file=bprotfiles(iprot),status="unknown",
1938      &    form="unformatted",access="direct",recl=lenrec(iprot))
1939 c Change AL 12/30/2017
1940         if (.not.mod_other_params)
1941      &  open (istat,file=benefiles(iprot),status="unknown",
1942      &    form="unformatted",access="direct",recl=lenrec_ene(iprot))
1943 c Read conformations from binary DA files (one per batch) and write them to
1944 c a binary DA scratchfile.
1945         jj=0
1946         jjj=0
1947         write (liczba,'(bz,i4.4)') me
1948 #ifdef MPI
1949         IF (ME.EQ.MASTER) THEN
1950 c Only the master reads the database; it'll send it to the other procs
1951 c through a ring.
1952 #endif
1953           t_acq = tcpu()
1954           icount=0
1955           itmp=0
1956           do if=1,nfile_prot(iprot)
1957             nazwa=protfiles(if,iprot)
1958      &        (:ilen(protfiles(if,iprot)))//".cx"
1959 #ifdef DEBUG
1960             write (iout,*) "Opening file ",nazwa(:ilen(nazwa))
1961 #endif
1962 #if (defined(AIX) && !defined(JUBL))
1963             call xdrfopen_(ixdrf,nazwa, "r", iret)
1964 #else
1965             call xdrfopen(ixdrf,nazwa, "r", iret)
1966 #endif
1967             if (iret.eq.0) goto 1111
1968             do i=1,maxstr
1969
1970               prec=10000.0
1971 #if (defined(AIX) && !defined(JUBL))
1972               call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
1973               if (iret.eq.0) goto 102
1974               call xdrfint_(ixdrf, nss, iret)
1975               if (iret.eq.0) goto 102
1976               do j=1,nss
1977                 call xdrfint_(ixdrf, ihpb(j), iret)
1978                 if (iret.eq.0) goto 102
1979                 call xdrfint_(ixdrf, jhpb(j), iret)
1980                 if (iret.eq.0) goto 102
1981               enddo
1982               call xdrffloat_(ixdrf,reini,iret)
1983               if (iret.eq.0) goto 102
1984               call xdrffloat_(ixdrf,refree,iret)
1985               if (iret.eq.0) goto 102
1986               call xdrfint_(ixdrf,natlik,iret)
1987               if (iret.eq.0) goto 102
1988               do j=1,natlik
1989                 call xdrfint(ixdrf,natliktemp(j),iret)
1990                 if (iret.eq.0) goto 102
1991                 do k=1,natliktemp(j)
1992                   call xdrffloat(ixdrf,nutemp(k,j),iret)
1993                   if (iret.eq.0) goto 102
1994                 enddo
1995               enddo
1996 #else
1997 c               write (0,*) "me",me," iprot",iprot," i",i
1998               call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret)
1999               if (iret.eq.0) goto 102
2000               call xdrfint(ixdrf, nss, iret)
2001               if (iret.eq.0) goto 102
2002               do k=1,nss
2003                 call xdrfint(ixdrf, ihpb(k), iret)
2004                  if (iret.eq.0) goto 102
2005                 call xdrfint(ixdrf, jhpb(k), iret)
2006                 if (iret.eq.0) goto 102
2007               enddo
2008               call xdrffloat(ixdrf,reini,iret)
2009               if (iret.eq.0) goto 102
2010               call xdrffloat(ixdrf,refree,iret)
2011               if (iret.eq.0) goto 102
2012 #ifdef FINAL
2013               call xdrfint(ixdrf,natlik,iret)
2014               if (iret.eq.0) goto 102
2015               do j=1,natlik
2016                 call xdrfint(ixdrf,natliktemp(j),iret)
2017                 if (iret.eq.0) goto 102
2018                 do k=1,natliktemp(j)
2019                   call xdrffloat(ixdrf,nutemp(k,j),iret)
2020                   if (iret.eq.0) goto 102
2021                 enddo
2022               enddo
2023 #else
2024               natlik=0
2025               call xdrffloat(ixdrf,rmsdev,iret)
2026               if (iret.eq.0) goto 102
2027 c               write (2,*) "rmsdev",rmsdev
2028               call xdrfint(ixdrf,iscor,iret)
2029               if (iret.eq.0) goto 102
2030 c               write (2,*) "iscor",iscor
2031 #endif
2032 #endif
2033               eini(jj+1,iprot)=reini
2034               entfac(jj+1,iprot)=refree
2035               do k=1,nres
2036                 do l=1,3
2037                   c(l,k)=csingle(l,k)
2038                 enddo
2039               enddo
2040               do k=nnt,nct
2041                 do l=1,3
2042                   c(l,nres+k)=csingle(l,nres+k-nnt+1)
2043                 enddo
2044               enddo
2045 #ifdef DEBUG
2046               write (iout,'(5hREAD ,i5,2f15.4)')
2047      &         jj+1,eini(jj+1,iprot),entfac(jj+1,iprot)
2048               write (iout,*) "natlik",natlik
2049               do l=1,natlik
2050                 write (iout,*) "natliketemp",natliktemp(l)
2051                 write(iout,*) (nutemp(k,l),k=1,natliktemp(l))
2052               enddo
2053               write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
2054               write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
2055               call flush(iout)
2056 #endif
2057               call add_new_cconf(jjj,jj,jj_old,icount,iprot,
2058      &            Next)
2059             enddo
2060   102       continue
2061             write (iout,*) "Protein ",protname(iprot),
2062      &        i-1," conformations read from DA file ",
2063      &        nazwa(:ilen(nazwa))
2064             write (iout,*) jj," conformations read so far"
2065 #if (defined(AIX) && !defined(JUBL))
2066             call xdrfclose_(ixdrf,nazwa,iret)
2067 #else
2068             call xdrfclose(ixdrf,nazwa,iret)
2069 #endif
2070 c              print *,"file ",nazwa(:ilen(nazwa))," closed"
2071           enddo
2072 #ifdef MPI
2073 #ifdef DEBUG
2074           write (iout,*) "jj_old",jj_old," jj",jj
2075 #endif
2076           call write_and_send_cconf(icount,jj_old,jj,iprot,Next)
2077           if (icount.gt.0) call MPI_Send(0,1,MPI_INTEGER,Next,570,
2078      &           WHAM_COMM,IERROR)
2079           jj_old=jj+1
2080 #else
2081           call write_and_send_cconf(icount,jj_old,jj,iprot,Next)
2082 #endif
2083           t_acq = tcpu() - t_acq
2084           write (iout,*) "Processor",me," protein",iprot,
2085      &        " batch",ibatch," time for conformation read/send",t_acq
2086 #ifdef MPI
2087         ELSE
2088 c A worker gets the confs from the master and sends them to its neighbor
2089         t_acq = tcpu()
2090         call receive_and_pass_cconf(icount,jj_old,jj,iprot,
2091      &       Previous,Next)
2092         t_acq = tcpu() - t_acq
2093           write (iout,*) "Processor",me," protein",iprot,
2094      &         " batch",ibatch,
2095      &         " time for conformation receive/send",t_acq
2096         ENDIF
2097 #endif
2098         ntot(iprot)=jj
2099         write (iout,*) "Protein",
2100      &    protname(iprot)(:ilen(protname(iprot))),", ",ntot(iprot),
2101      &   " conformatons read ",ntot(iprot)," conformations written to ",
2102      &  bprotfiles(iprot)(:ilen(bprotfiles(iprot)))
2103         ntot(0) = ntot(0)+ntot(iprot)
2104         close(ientout)
2105         close(istat)
2106       return
2107  1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
2108       call flush(iout)
2109 #ifdef MPI
2110       call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2111 #endif
2112       end
2113 c------------------------------------------------------------------------------
2114       subroutine add_new_cconf(jjj,jj,jj_old,icount,iprot,Next)
2115       implicit none
2116       include "DIMENSIONS"
2117       include "DIMENSIONS.ZSCOPT"
2118       include "COMMON.CHAIN"
2119       include "COMMON.INTERACT"
2120       include "COMMON.LOCAL"
2121       include "COMMON.CLASSES"
2122       include "COMMON.ENERGIES"
2123       include "COMMON.IOUNITS"
2124       include "COMMON.PROTFILES"
2125       include "COMMON.PROTNAME"
2126       include "COMMON.VMCPAR"
2127       include "COMMON.WEIGHTS"
2128       include "COMMON.NAMES"
2129       include "COMMON.ALLPROT"
2130       include "COMMON.WEIGHTDER"
2131       include "COMMON.VAR"
2132       include "COMMON.SBRIDGE"
2133       include "COMMON.GEO"
2134       include "COMMON.COMPAR"
2135 #ifdef ISNAN
2136       logical isnan
2137 #endif
2138       integer i,j,jj,jjj,jj_old,icount,k,kk,l,iprot,ii,ib,ibatch,
2139      &  nn,nn1,inan,Next,itj
2140       double precision etot,energia(0:max_ene)
2141       jjj=jjj+1
2142 c      if (entfac(jj+1,iprot).gt.-10.0d0
2143 c     &     .or. entfac(jj+1,iprot).lt.-150.0d0) then
2144 c        write (iout,*) "Entropy factor out of range for conformation",
2145 c     &    jjj,entfac(jj+1,iprot),", conformation skipped."
2146 c        return
2147 c      endif
2148       call int_from_cart1(.false.)
2149       do j=nnt+1,nct
2150         if (vbld(j).lt.2.0d0 .or. vbld(j).gt.6.5d0) then
2151           write (iout,*) "nnt",nnt," nct",nct
2152           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j)
2153           write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
2154           write (iout,*) "The Cartesian geometry is:"
2155           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
2156           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
2157           write (iout,*) "The internal geometry is:"
2158           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2159           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2160           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2161           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2162           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2163           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2164           write (iout,*)
2165      &      "This conformation WILL NOT be added to the database."
2166           return
2167         endif
2168       enddo
2169       do j=nnt,nct
2170         itj=itype(j)
2171         if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
2172           write (iout,*) "nnt",nnt," nct",nct
2173           write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
2174           write (iout,*) "The Cartesian geometry is:"
2175           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
2176           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
2177           write (iout,*) "The internal geometry is:"
2178           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2179           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2180           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2181           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2182           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2183           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2184           write (iout,*)
2185      &      "This conformation WILL NOT be added to the database."
2186           return
2187         endif
2188       enddo
2189       do j=3,nres
2190         if (theta(j).le.0.0d0) then
2191           write (iout,*)
2192      &      "Zero theta angle(s) in conformation",ii
2193           write (iout,*) "The Cartesian geometry is:"
2194           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
2195           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
2196           write (iout,*) "The internal geometry is:"
2197           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2198           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2199           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2200           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2201           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2202           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2203           write (iout,*)
2204      &      "This conformation WILL NOT be added to the database."
2205           return
2206         endif
2207         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
2208       enddo
2209       if (.not. init_ene) then
2210       etot=0.0d0
2211       do j=1,n_ene
2212         etot=etot+ww(j)*enetb(icount+1,j,iprot)
2213       enddo
2214       inan=0
2215 c detecting NaNQ
2216 #ifdef ISNAN
2217 #ifdef AIX
2218       if (isnan(etot).ne.0) inan=1
2219 #else
2220       if (isnan(etot)) inan=1
2221 #endif
2222 #else
2223 #ifdef WINPGI
2224       idumm=proc_proc(etot,inan)
2225 #else
2226       call proc_proc(etot,inan)
2227 #endif
2228 #endif
2229       else
2230       inan=0
2231       endif
2232       if (inan.eq.1) then
2233         write (iout,*) "NaNs detected in some of the energy",
2234      &    " components for protein",iprot," batch",ibatch,
2235      &    " conformation",jjj
2236         write (iout,*) "etot",etot
2237         write (iout,*) "The Cartesian geometry is:"
2238         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2239         write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
2240         write (iout,*) "The internal geometry is:"
2241         write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2242         write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2243         write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2244         write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2245         write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
2246         write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2247         write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2248         write (iout,*) "The components of the energy are:"
2249         energia(0)=etot
2250         do k=1,n_ene
2251           energia(k)=enetb(jj+1,k,iprot)
2252         enddo
2253         call enerprint(energia(0))
2254         write (iout,*)
2255      &    "This conformation WILL NOT be added to the database."
2256         call flush(iout)
2257       else
2258         jj=jj+1
2259 #ifdef DEBUG
2260         write (iout,'(e15.5,16i5)') entfac(icount+1,iprot),
2261      &        iscore(icount+1,0,iprot),ibatch
2262         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2263         write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
2264         write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2265         write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2266         write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2267         write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2268         write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
2269         write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2270         write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2271         write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
2272         write (iout,'(5e15.5)') (enetb(icount+1,j,iprot),j=1,n_ene)
2273 c        write (iout,'(2e15.5)') (eneps(1,j,icount+1,iprot),
2274 c     &        eneps(2,j,icount+1,iprot),j=1,nntyp)
2275 #endif
2276 c        write (iout,*) "First nu in readrtms"
2277         call flush(iout)
2278         icount=icount+1
2279         list_conf(jj,iprot)=jj
2280         call store_cconf_from_file(jj,icount,iprot)
2281         do j=1,natlike(iprot)
2282           do k=1,natlik
2283             if (k.eq.numnat(j,iprot)) then
2284               do l=1,natdim(j,iprot)
2285                 nu(l,k,j,iprot)=nutemp(l,k)
2286               enddo
2287               exit
2288             endif
2289           enddo
2290         enddo
2291 c        write (iout,*) "End First nu in readrtms"
2292         call flush(iout)
2293         if (icount.eq.maxstr_proc) then
2294 #ifdef DEBUG
2295           write (iout,* ) "jj_old",jj_old," jj",jj
2296           write (iout,*) "Calling write_and_send_cconf"
2297           call flush(iout)
2298 #endif
2299           call write_and_send_cconf(icount,jj_old,jj,iprot,Next)
2300           jj_old=jj+1
2301 #ifdef DEBUG
2302           write (iout,*) "Exited write_and_send_cconf"
2303           call flush(iout)
2304 #endif
2305           icount=0
2306         endif
2307       endif
2308       return
2309       end
2310 c------------------------------------------------------------------------------
2311       subroutine store_cconf_from_file(jj,icount,iprot)
2312       implicit none
2313       include "DIMENSIONS"
2314       include "DIMENSIONS.ZSCOPT"
2315       include "COMMON.CHAIN"
2316       include "COMMON.SBRIDGE"
2317       include "COMMON.INTERACT"
2318       include "COMMON.IOUNITS"
2319       include "COMMON.CLASSES"
2320       include "COMMON.ALLPROT"
2321       include "COMMON.VAR"
2322       integer i,j,jj,icount,ibatch,iprot
2323 c Store the conformation that has been read in
2324       do i=1,2*nres
2325         do j=1,3
2326           c_zs(j,i,icount,iprot)=c(j,i)
2327         enddo
2328       enddo
2329       nss_zs(icount,iprot)=nss
2330       do i=1,nss
2331         ihpb_zs(i,icount,iprot)=ihpb(i)
2332         jhpb_zs(i,icount,iprot)=jhpb(i)
2333       enddo
2334       return
2335       end
2336 c------------------------------------------------------------------------------
2337       subroutine write_and_send_cconf(icount,jj_old,jj,iprot,Next)
2338       implicit none
2339       include "DIMENSIONS"
2340       include "DIMENSIONS.ZSCOPT"
2341 #ifdef MPI
2342       include "mpif.h"
2343       integer IERROR
2344       include "COMMON.MPI"
2345 #endif
2346       include "COMMON.WEIGHTS"
2347       include "COMMON.CHAIN"
2348       include "COMMON.SBRIDGE"
2349       include "COMMON.INTERACT"
2350       include "COMMON.IOUNITS"
2351       include "COMMON.CLASSES"
2352       include "COMMON.VAR"
2353       include "COMMON.ALLPROT"
2354       include "COMMON.ENERGIES"
2355       include "COMMON.WEIGHTDER"
2356       include "COMMON.OPTIM"
2357       include "COMMON.COMPAR"
2358       integer icount,jj_old,jj,Next,iprot
2359 c Write the structures to a scratch file
2360 #ifdef MPI
2361 c Master sends the portion of conformations that have been read in to the neighbor
2362 #ifdef DEBUG
2363       write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
2364       call flush(iout)
2365 #endif
2366       call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM,IERROR)
2367 #ifdef DEBUG
2368       write (iout,*) "Processor",me," Next",next," sent icount=",icount
2369       write (iout,*) "Processor",me," jj_old",jj_old," jj",jj
2370       call flush(iout)
2371 #endif
2372       if (icount.gt.0) then
2373       call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER,
2374      &    Next,571,WHAM_COMM,IERROR)
2375       call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2376      &    Next,572,WHAM_COMM,IERROR)
2377       if (.not. init_ene) call MPI_Send(enetb(jj_old,1,iprot),
2378      &    maxstr*n_ene,
2379      &    MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,IERROR)
2380       call MPI_Send(nu(1,1,jj_old,iprot),
2381      &    maxdimnat*natlike(iprot)*icount,
2382      &    MPI_DOUBLE_PRECISION,
2383      &    Next,577,WHAM_COMM,IERROR)
2384       call MPI_Send(eini(jj_old,iprot),icount,
2385      &    MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR)
2386       call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2387      &    Next,579,WHAM_COMM,IERROR)
2388       call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2,
2389      &    MPI_REAL,Next,580,WHAM_COMM,IERROR)
2390       if (.not. init_ene) call MPI_Send(eneps(1,1,1,iprot),
2391      &    2*icount*nntyp,
2392      &    MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR)
2393       endif
2394 #endif
2395       call flush(iout)
2396       call dawrite_ccoords(iprot,jj_old,jj,ientout)
2397 c Change AL 20/12/2017
2398       if (.not. mod_other_params)
2399      &call dawrite_ene(iprot,jj_old,jj,istat)
2400       return
2401       end
2402 c------------------------------------------------------------------------------
2403 #ifdef MPI
2404       subroutine receive_and_pass_cconf(icount,jj_old,jj,iprot,Previous,
2405      &  Next)
2406       implicit none
2407       include "DIMENSIONS"
2408       include "DIMENSIONS.ZSCOPT"
2409       include "mpif.h"
2410       integer IERROR,STATUS(MPI_STATUS_SIZE)
2411       include "COMMON.MPI"
2412       include "COMMON.CHAIN"
2413       include "COMMON.SBRIDGE"
2414       include "COMMON.INTERACT"
2415       include "COMMON.IOUNITS"
2416       include "COMMON.CLASSES"
2417       include "COMMON.ALLPROT"
2418       include "COMMON.ENERGIES"
2419       include "COMMON.VAR"
2420       include "COMMON.GEO"
2421       include "COMMON.WEIGHTS"
2422       include "COMMON.WEIGHTDER"
2423       include "COMMON.COMPAR"
2424       include "COMMON.OPTIM"
2425       integer i,j,k,icount,jj_old,jj,iprot,Previous,Next
2426       icount=1
2427 #ifdef DEBUG
2428       write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
2429       call flush(iout)
2430 #endif
2431       do while (icount.gt.0)
2432 c      call MPI_Probe(Previous,570,WHAM_COMM,STATUS,IERROR)
2433       call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,WHAM_COMM,
2434      &     STATUS,IERROR)
2435 #ifdef DEBUG
2436       write (iout,*)"Processor",me," previous",previous," icount",icount
2437       call flush(iout)
2438 #endif
2439       call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM,
2440      &     IERROR)
2441 #ifdef DEBUG
2442       write (iout,*) "Processor",me," icount",icount
2443       call flush(iout)
2444 #endif
2445       if (icount.eq.0) return
2446       call MPI_Recv(nss_zs(1,iprot),icount,MPI_INTEGER,
2447      &    Previous,571,WHAM_COMM,STATUS,IERROR)
2448       call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER,
2449      &  Next,571,WHAM_COMM,IERROR)
2450       call MPI_Recv(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2451      &    Previous,572,WHAM_COMM,STATUS,IERROR)
2452       call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2453      &  Next,572,WHAM_COMM,IERROR)
2454       if (.not. init_ene) then
2455         call MPI_Recv(enetb(jj_old,1,iprot),maxstr*n_ene,
2456      &   MPI_DOUBLE_PRECISION,Previous,576,WHAM_COMM,STATUS,IERROR)
2457         call MPI_Send(enetb(jj_old,1,iprot),maxstr*n_ene,
2458      &   MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,STATUS,IERROR)
2459       endif
2460       call MPI_Recv(nu(1,1,jj_old,iprot),
2461      &  maxdimnat*natlike(iprot)*icount,
2462      &  MPI_DOUBLE_PRECISION,
2463      &  Previous,577,WHAM_COMM,STATUS,IERROR)
2464       call MPI_Send(nu(1,1,jj_old,iprot),
2465      &  maxdimnat*natlike(iprot)*icount,
2466      &  MPI_DOUBLE_PRECISION,
2467      &  Next,577,WHAM_COMM,IERROR)
2468       call MPI_Recv(eini(jj_old,iprot),icount,
2469      &  MPI_DOUBLE_PRECISION,Previous,578,WHAM_COMM,STATUS,IERROR)
2470       call MPI_Send(eini(jj_old,iprot),icount,
2471      &  MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR)
2472       call MPI_Recv(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2473      &  Previous,579,WHAM_COMM,STATUS,IERROR)
2474       call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2475      &  Next,579,WHAM_COMM,IERROR)
2476       call MPI_Recv(c_zs(1,1,1,iprot),3*icount*maxres2,
2477      &  MPI_REAL,Previous,580,WHAM_COMM,STATUS,IERROR)
2478       call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2,
2479      &  MPI_REAL,Next,580,WHAM_COMM,IERROR)
2480       if (.not. init_ene) then
2481         call MPI_Recv(eneps(1,1,1,iprot),2*icount*nntyp,
2482      &   MPI_DOUBLE_PRECISION,Previous,582,WHAM_COMM,STATUS,IERROR)
2483         call MPI_Send(eneps(1,1,1,iprot),2*icount*nntyp,
2484      &   MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR)
2485       endif
2486       jj=jj_old+icount-1
2487       do i=jj_old,jj
2488         list_conf(i,iprot)=i
2489       enddo
2490       call dawrite_ccoords(iprot,jj_old,jj,ientout)
2491 c Change AL 12/20/2017
2492       if (.not. mod_other_params)
2493      &call dawrite_ene(iprot,jj_old,jj,istat)
2494       jj_old=jj+1
2495 #ifdef DEBUG
2496       write (iout,*) "Protein",iprot
2497       write (iout,*) "Processor",me," received",icount," conformations"
2498       do i=1,icount
2499         write (iout,'(8f10.4)') ((c_zs(l,k,i,iprot),l=1,3,k=1,nres)
2500         write (iout,'(8f10.4)')((c_zs(l,k,i+nres,iprot),l=1,3,k=nnt,nct)
2501         write (iout,'(16i5)') nss_zs(i,iprot),(ihpb_zs(k,i,iprot),
2502      &    jhpb_zs(k,i,iprot),k=1,nss_zs(i,iprot))
2503         write (iout,'(5e15.5)') (enetb(i,j,iprot),j=1,n_ene)
2504         write (iout,'(2e15.5)') (eneps(1,j,i,iprot),
2505      &        eneps(2,j,i,iprot),j=1,nntyp)
2506         write (iout,'(e15.5,16i5)') entfac(i),iscore(i,0,iprot),
2507      &        kbatch(i,iprot)
2508       enddo
2509 #endif
2510       enddo
2511       return
2512       end
2513 #endif
2514 c------------------------------------------------------------------------------
2515       subroutine read_thermal
2516       implicit none
2517       include "DIMENSIONS"
2518       include "DIMENSIONS.ZSCOPT"
2519       include "COMMON.CHAIN"
2520       include "COMMON.SBRIDGE"
2521       include "COMMON.INTERACT"
2522       include "COMMON.IOUNITS"
2523       include "COMMON.CLASSES"
2524       include "COMMON.VAR"
2525       include "COMMON.THERMAL"
2526       character*800 controlcard
2527       call card_concat(controlcard,.true.)
2528       call readi(controlcard,"NGRIDT",NGridT,200)
2529       call reada(controlcard,"DELTAT",deltaT,5.0d0)
2530       call reada(controlcard,"T0",GridT0,2.0d2)
2531       write (iout,*) "Parameters of thermal curves"
2532       write (iout,*) "NGRIDT",NGridT," DELTAT",deltaT," T0",GridT0
2533       return
2534       end
2535 c------------------------------------------------------------------------------
2536       subroutine daread_ene(iprot,istart_conf,iend_conf)
2537       implicit none
2538       include "DIMENSIONS"
2539       include "DIMENSIONS.ZSCOPT"
2540 #ifdef MPI
2541       include "mpif.h"
2542       include "COMMON.MPI"
2543 #endif
2544       include "COMMON.CHAIN"
2545       include "COMMON.CLASSES"
2546       include "COMMON.ENERGIES"
2547       include "COMMON.IOUNITS"
2548       include "COMMON.PROTFILES"
2549       include "COMMON.ALLPROT"
2550       include "COMMON.WEIGHTDER"
2551       include "COMMON.COMPAR"
2552       include "COMMON.VMCPAR"
2553       integer iprot,istart_conf,iend_conf
2554       integer i,ii,iii,j,k
2555 #ifdef DEBUG
2556       write (iout,*) "Calling DAREAD_ENE"
2557 #endif
2558 c      write (iout,*) "Reading: n_ene",n_ene
2559 c      call flush(iout)
2560 c      write (iout,*) "DAREAD_ENE",istart_conf,iend_conf
2561 c
2562 c Read conformations off a DA scratchfile.
2563 c
2564       do ii=istart_conf,iend_conf
2565         iii=list_conf(ii,iprot)
2566         i = ii - istart_conf + 1
2567         read(ientin,rec=iii) (enetb(i,j,iprot),j=1,n_ene),
2568      &     (enetb_orig(i,j,iprot),j=1,n_ene),
2569      &     (enetb_oorig(i,j,iprot),j=1,n_ene),
2570      &     eini(ii,iprot),entfac(ii,iprot),
2571      &     (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp),
2572      &     ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2573 #ifdef DEBUG
2574         write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot),
2575      &  entfac(ii,iprot),
2576         write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene)
2577         write (iout,'(18(1pe12.4))')
2578      &  ((nu(k,j,i,iprot)k=1,maxdimnat),j=1,natlike(iprot))
2579         call flush(iout)
2580 #endif
2581       enddo
2582
2583       return
2584       end
2585 c------------------------------------------------------------------------------
2586       subroutine dawrite_ene(iprot,istart_conf,iend_conf,unit_out)
2587       implicit none
2588       include "DIMENSIONS"
2589       include "DIMENSIONS.ZSCOPT"
2590 #ifdef MPI
2591       include "mpif.h"
2592       include "COMMON.MPI"
2593 #endif
2594       include "COMMON.CHAIN"
2595       include "COMMON.CLASSES"
2596       include "COMMON.ENERGIES"
2597       include "COMMON.IOUNITS"
2598       include "COMMON.PROTFILES"
2599       include "COMMON.ALLPROT"
2600       include "COMMON.WEIGHTDER"
2601       include "COMMON.VMCPAR"
2602       include "COMMON.COMPAR"
2603       integer iprot,istart_conf,iend_conf,unit_out
2604       integer i,ii,iii,j,k
2605 c      write (iout,*) "Writing: n_ene",n_ene
2606 c      call flush(iout)
2607 c      write (iout,*) "DAWRITE_ENE",istart_conf,iend_conf
2608 c
2609 c Write conformations to a DA scratchfile.
2610 c
2611       do ii=istart_conf,iend_conf
2612         iii=list_conf(ii,iprot)
2613         i = ii - istart_conf + 1
2614         write(unit_out,rec=iii) (enetb(i,j,iprot),j=1,n_ene),
2615      &     (enetb_orig(i,j,iprot),j=1,n_ene),
2616      &     (enetb_oorig(i,j,iprot),j=1,n_ene),
2617      &     eini(ii,iprot),entfac(ii,iprot),
2618      &     (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp),
2619      &     ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2620 #ifdef DEBUG
2621         write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot),
2622      &  entfac(ii,iprot)
2623         write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene)
2624         write (iout,'(18(1pe12.4))')
2625      &  ((nu(kj,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2626         call flush(iout)
2627 #endif
2628       enddo
2629
2630       return
2631       end
2632 c------------------------------------------------------------------------------
2633       subroutine daread_forces(iprot,istart_conf,iend_conf)
2634       implicit none
2635       include "DIMENSIONS"
2636       include "DIMENSIONS.ZSCOPT"
2637 #ifdef MPI
2638       include "mpif.h"
2639       include "COMMON.MPI"
2640 #endif
2641       include "COMMON.CHAIN"
2642       include "COMMON.CLASSES"
2643       include "COMMON.ENERGIES"
2644       include "COMMON.IOUNITS"
2645       include "COMMON.PROTFILES"
2646       include "COMMON.ALLPROT"
2647       include "COMMON.WEIGHTDER"
2648       include "COMMON.COMPAR"
2649       include "COMMON.VMCPAR"
2650       include "COMMON.FMATCH"
2651       include "COMMON.NAMES"
2652       include "COMMON.INTERACT"
2653       integer iprot,istart_conf,iend_conf
2654       integer i,ii,iii,j,jj,k,l
2655       integer kiri
2656       character*16 forma,acc
2657       character*128 nam
2658 #ifdef DEBUG
2659       write (iout,*) "Calling DAREAD_FORCES"
2660 #endif
2661 c      write (iout,*) "Reading: n_ene",n_ene
2662 c      call flush(iout)
2663 c      write (iout,*) "DAREAD_ENE",istart_conf,iend_conf
2664 c
2665 c Read conformations off a DA scratchfile.
2666 c
2667 #ifdef DEBUG
2668       inquire(unit=ientin,name=nam,recl=kiri,form=forma,access=acc)
2669       write (iout,*) "DAREAD_FORCES","len=",kiri," form=",forma,
2670      &  " acc=",acc
2671       write (iout,*) "nam=",nam
2672 #endif
2673       do ii=istart_conf,iend_conf
2674         iii=list_conf_MD(ii,iprot)
2675         i = ii - istart_conf + 1
2676         read(ientin,rec=iii) (((forctb(j,l,k,i,iprot),j=1,n_ene),l=1,3),
2677      &     k=1,nsite(iprot))
2678 #ifdef DEBUG
2679         write (iout,'(a,4i5)') "DAREAD_FORCES",ii,iii,i,nsite(iprot)
2680         write (iout,*) "CG force components"
2681         do k=1,n_ene
2682
2683           write (iout,'(a,i3,1x,a)') "Component",k,ename(k)
2684
2685           jj=nres
2686           do j=1,nres
2687             if (itype(j).eq.10 .or. itype(j).eq.ntyp1) then
2688               write (iout,'(a4,i5,3e15.5)')
2689      &        restyp(itype(j)),j,(forctb(k,l,j,i,iprot),l=1,3)
2690             else
2691               jj=jj+1
2692               write (iout,'(a4,i5,3e15.5,5x,3e15.5)')
2693      &        restyp(itype(j)),j,(forctb(k,l,j,i,iprot),l=1,3),
2694      &        (forctb(k,l,jj,i,iprot),l=1,3)
2695             endif
2696           enddo
2697
2698         enddo
2699         call flush(iout)
2700 #endif
2701       enddo
2702
2703       return
2704       end
2705 c------------------------------------------------------------------------------
2706       subroutine dawrite_forces(iprot,istart_conf,iend_conf,unit_out)
2707       implicit none
2708       include "DIMENSIONS"
2709       include "DIMENSIONS.ZSCOPT"
2710 #ifdef MPI
2711       include "mpif.h"
2712       include "COMMON.MPI"
2713 #endif
2714       include "COMMON.CHAIN"
2715       include "COMMON.CLASSES"
2716       include "COMMON.ENERGIES"
2717       include "COMMON.IOUNITS"
2718       include "COMMON.PROTFILES"
2719       include "COMMON.ALLPROT"
2720       include "COMMON.WEIGHTDER"
2721       include "COMMON.VMCPAR"
2722       include "COMMON.COMPAR"
2723       include "COMMON.FMATCH"
2724       include "COMMON.NAMES"
2725       include "COMMON.INTERACT"
2726       integer iprot,istart_conf,iend_conf,unit_out
2727       integer i,ii,iii,j,jj,k,l
2728       integer kiri
2729       character*16 forma,acc
2730       character*128 nam
2731 c      write (iout,*) "Writing: n_ene",n_ene
2732 c      call flush(iout)
2733 c      write (iout,*) "DAWRITE_ENE",istart_conf,iend_conf
2734 c
2735 c Write conformations to a DA scratchfile.
2736 c
2737 #ifdef DEBUG
2738       inquire(unit=unit_out,name=nam,recl=kiri,form=forma,access=acc)
2739       write (iout,*) "DAWRITE_FORCES","len=",kiri," form=",forma,
2740      &  " acc=",acc
2741       write (iout,*) "nam=",nam
2742       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2743       write (iout,*) "list",(list_conf_MD(k,iprot),
2744      & k=istart_conf,iend_conf)
2745 #endif
2746       do ii=istart_conf,iend_conf
2747         iii=list_conf_MD(ii,iprot)
2748         i = ii - istart_conf + 1
2749         write(unit_out,rec=iii) (((forctb(j,l,k,i,iprot),j=1,n_ene),
2750      &     l=1,3),k=1,nsite(iprot))
2751 #ifdef DEBUG
2752         write (iout,'(a,4i5)') "DAWRITE_FORCES",ii,iii,i,nsite(iprot)
2753         write (iout,*) "CG force components"
2754         do k=1,n_ene
2755
2756           write (iout,'(a,i3,1x,a)') "Component",k,ename(k)
2757
2758           jj=nres
2759           do j=1,nres
2760             if (itype(j).eq.10 .or. itype(j).eq.ntyp1) then
2761               write (iout,'(a4,i5,3e15.5)')
2762      &        restyp(itype(j)),j,(forctb(k,l,j,i,iprot),l=1,3)
2763             else
2764               jj=jj+1
2765               write (iout,'(a4,i5,3e15.5,5x,3e15.5)')
2766      &        restyp(itype(j)),j,(forctb(k,l,j,i,iprot),l=1,3),
2767      &        (forctb(k,l,jj,i,iprot),l=1,3)
2768             endif
2769           enddo
2770
2771         enddo
2772         call flush(iout)
2773 #endif
2774       enddo
2775
2776       return
2777       end
2778 c------------------------------------------------------------------------------
2779       subroutine daread_ccoords(iprot,istart_conf,iend_conf)
2780       implicit none
2781       include "DIMENSIONS"
2782       include "DIMENSIONS.ZSCOPT"
2783 #ifdef MPI
2784       include "mpif.h"
2785       include "COMMON.MPI"
2786 #endif
2787       include "COMMON.CHAIN"
2788       include "COMMON.CLASSES"
2789       include "COMMON.ENERGIES"
2790       include "COMMON.IOUNITS"
2791       include "COMMON.PROTFILES"
2792       include "COMMON.ALLPROT"
2793       include "COMMON.INTERACT"
2794       include "COMMON.VAR"
2795       include "COMMON.SBRIDGE"
2796       include "COMMON.GEO"
2797       include "COMMON.COMPAR"
2798       include "COMMON.VMCPAR"
2799       include "COMMON.WEIGHTDER"
2800       integer iprot,istart_conf,iend_conf
2801       integer i,j,k,ij,ii,iii
2802       integer len
2803       real*4 csingle(3,maxres2)
2804       character*16 form,acc
2805       character*32 nam
2806 c
2807 c Read conformations off a DA scratchfile.
2808 c
2809 #ifdef DEBUG
2810       write (iout,*) "DAREAD_COORDS"
2811       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2812       write (iout,*) "lenrec",lenrec(iprot)
2813       inquire(unit=ientin,name=nam,recl=len,form=form,access=acc)
2814       write (iout,*) "len=",len," form=",form," acc=",acc
2815       write (iout,*) "nam=",nam
2816       call flush(iout)
2817 #endif
2818       do ii=istart_conf,iend_conf
2819         iii=list_conf(ii,iprot)
2820         ij = ii - istart_conf + 1
2821 #ifdef DEBUG
2822         write (iout,*) "Reading binary file, record",iii," ii",ii
2823         call flush(iout)
2824 #endif
2825         read(icbase,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2826      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2827      &    nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot),
2828      &    entfac(ii,iprot),
2829      &    ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot))
2830         do i=1,2*nres
2831           do j=1,3
2832             c(j,i)=csingle(j,i)
2833           enddo
2834         enddo
2835 #ifdef DEBUG
2836         write (iout,*) "iprot",iprot," ii",ii
2837         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2838         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2839         write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot)
2840         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2841         call flush(iout)
2842 #endif
2843         call store_ccoords(iprot,ii-istart_conf+1)
2844       enddo
2845       return
2846       end
2847 c------------------------------------------------------------------------------
2848       subroutine dawrite_ccoords(iprot,istart_conf,iend_conf,unit_out)
2849       implicit none
2850       include "DIMENSIONS"
2851       include "DIMENSIONS.ZSCOPT"
2852 #ifdef MPI
2853       include "mpif.h"
2854       include "COMMON.MPI"
2855 #endif
2856       include "COMMON.CHAIN"
2857       include "COMMON.INTERACT"
2858       include "COMMON.CLASSES"
2859       include "COMMON.ENERGIES"
2860       include "COMMON.IOUNITS"
2861       include "COMMON.PROTFILES"
2862       include "COMMON.ALLPROT"
2863       include "COMMON.VAR"
2864       include "COMMON.SBRIDGE"
2865       include "COMMON.GEO"
2866       include "COMMON.COMPAR"
2867       include "COMMON.WEIGHTDER"
2868       include "COMMON.VMCPAR"
2869       real*4 csingle(3,maxres2)
2870       integer iprot,istart_conf,iend_conf
2871       integer i,j,k,ii,ij,iii,unit_out
2872       integer len
2873       character*16 form,acc
2874       character*32 nam
2875 c
2876 c Write conformations to a DA scratchfile.
2877 c
2878 #ifdef DEBUG
2879       write (iout,*) "DAWRITE_COORDS"
2880       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2881       write (iout,*) "lenrec",lenrec(iprot)
2882       inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
2883       write (iout,*) "len=",len," form=",form," acc=",acc
2884       write (iout,*) "nam=",nam
2885       call flush(iout)
2886 #endif
2887       do ii=istart_conf,iend_conf
2888         iii=list_conf(ii,iprot)
2889         ij = ii - istart_conf + 1
2890         call restore_ccoords(iprot,ii-istart_conf+1)
2891 #ifdef DEBUG
2892         write (iout,*) "Writing binary file, record",iii," ii",ii
2893         call flush(iout)
2894 #endif
2895         do i=1,2*nres
2896           do j=1,3
2897             csingle(j,i)=c(j,i)
2898           enddo
2899         enddo
2900         write(unit_out,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2901      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2902      &    nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot),
2903      &    entfac(ii,iprot),
2904      &    ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot))
2905 #ifdef DEBUG
2906         write (iout,*) "kbatch",kbatch(ii,iprot)," iscore",
2907      &    iscore(ii,0,iprot)
2908         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2909         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2910         write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot)
2911         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2912         call flush(iout)
2913 #endif
2914       enddo
2915       return
2916       end
2917 c------------------------------------------------------------------------------
2918       subroutine daread_MDcoords(iprot,istart_conf,iend_conf)
2919       implicit none
2920       include "DIMENSIONS"
2921       include "DIMENSIONS.ZSCOPT"
2922 #ifdef MPI
2923       include "mpif.h"
2924       include "COMMON.MPI"
2925 #endif
2926       include "COMMON.CHAIN"
2927       include "COMMON.CLASSES"
2928       include "COMMON.ENERGIES"
2929       include "COMMON.IOUNITS"
2930       include "COMMON.PROTFILES"
2931       include "COMMON.ALLPROT"
2932       include "COMMON.INTERACT"
2933       include "COMMON.VAR"
2934       include "COMMON.SBRIDGE"
2935       include "COMMON.GEO"
2936       include "COMMON.COMPAR"
2937       include "COMMON.VMCPAR"
2938       include "COMMON.WEIGHTDER"
2939       include "COMMON.FMATCH"
2940       integer iprot,istart_conf,iend_conf
2941       integer i,j,k,ij,ii,iii
2942       integer len
2943       real*4 csingle(3,maxres2),fsingle(3,maxres2)
2944       character*16 form,acc
2945       character*32 nam
2946 c
2947 c Read conformations off a DA scratchfile.
2948 c
2949 #ifdef DEBUG
2950       write (iout,*) "DAREAD_MDCOORDS"
2951       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2952       write (iout,*) "lenrec_MD",lenrec_MD(iprot)
2953       inquire(unit=ientin,name=nam,recl=len,form=form,access=acc)
2954       write (iout,*) "len=",len," form=",form," acc=",acc
2955       write (iout,*) "nam=",nam
2956       call flush(iout)
2957 #endif
2958       do ii=istart_conf,iend_conf
2959         iii=list_conf_MD(ii,iprot)
2960         ij = ii - istart_conf + 1
2961 #ifdef DEBUG
2962         write (iout,*) "Reading binary file, record",iii," ii",ii
2963         call flush(iout)
2964 #endif
2965         read(icbase,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2966      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2967      &    ((CGforce_MD(j,i,ii,iprot),j=1,3),i=1,nsite(iprot)),
2968      &    nss,(ihpb(i),jhpb(i),i=1,nss)
2969         do i=1,2*nres
2970           do j=1,3
2971             c(j,i)=csingle(j,i)
2972           enddo
2973         enddo
2974 #ifdef DEBUG
2975         write (iout,*) "iprot",iprot," ii",ii
2976         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2977         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2978         write (iout,'(8f10.5)') ((CGforce_MD(j,i,ii,iprot),j=1,3),i=1,
2979      &    nsite(iprot))
2980         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2981         call flush(iout)
2982 #endif
2983         call store_MDcoords(iprot,ii-istart_conf+1)
2984       enddo
2985       return
2986       end
2987 c------------------------------------------------------------------------------
2988       subroutine dawrite_MDcoords(iprot,istart_conf,iend_conf,unit_out)
2989       implicit none
2990       include "DIMENSIONS"
2991       include "DIMENSIONS.ZSCOPT"
2992 #ifdef MPI
2993       include "mpif.h"
2994       include "COMMON.MPI"
2995 #endif
2996       include "COMMON.CHAIN"
2997       include "COMMON.INTERACT"
2998       include "COMMON.CLASSES"
2999       include "COMMON.ENERGIES"
3000       include "COMMON.IOUNITS"
3001       include "COMMON.PROTFILES"
3002       include "COMMON.ALLPROT"
3003       include "COMMON.VAR"
3004       include "COMMON.SBRIDGE"
3005       include "COMMON.GEO"
3006       include "COMMON.COMPAR"
3007       include "COMMON.WEIGHTDER"
3008       include "COMMON.VMCPAR"
3009       include "COMMON.FMATCH"
3010       real*4 csingle(3,maxres2)
3011       integer iprot,istart_conf,iend_conf
3012       integer i,j,k,ii,ij,iii,unit_out
3013       integer len
3014       character*16 form,acc
3015       character*32 nam
3016 c
3017 c Write conformations to a DA scratchfile.
3018 c
3019 #ifdef DEBUG
3020       write (iout,*) "DAWRITE_MDCOORDS"
3021       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
3022       write (iout,*) "lenrec_MD",lenrec_MD(iprot)
3023       inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
3024       write (iout,*) "len=",len," form=",form," acc=",acc
3025       write (iout,*) "nam=",nam
3026       call flush(iout)
3027 #endif
3028       do ii=istart_conf,iend_conf
3029         iii=list_conf_MD(ii,iprot)
3030         ij = iii - istart_conf + 1
3031         call restore_MDcoords(iprot,ii-istart_conf+1)
3032 #ifdef DEBUG
3033         write (iout,*) "Writing binary file, record",iii," ii",ii,
3034      &  " ij",ij
3035         write (iout,*) "nss",nss
3036         call flush(iout)
3037 #endif
3038         do i=1,2*nres
3039           do j=1,3
3040             csingle(j,i)=c(j,i)
3041           enddo
3042         enddo
3043         write(unit_out,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
3044      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
3045      &    ((CGforce_MD(j,i,ij,iprot),j=1,3),i=1,nsite(iprot)),
3046      &    nss,(ihpb(i),jhpb(i),i=1,nss)
3047 #ifdef DEBUG
3048         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
3049         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
3050         write (iout,'(8f10.5)') ((CGforce_MD(j,i,ij,iprot),j=1,3),
3051      &    i=1,nres)
3052         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
3053         call flush(iout)
3054 #endif
3055       enddo
3056       return
3057       end
3058 c------------------------------------------------------------------------------
3059       subroutine store_ccoords(iprot,ii)
3060       implicit none
3061       include "DIMENSIONS"
3062       include "DIMENSIONS.ZSCOPT"
3063       include "COMMON.VAR"
3064       include "COMMON.CHAIN"
3065       include "COMMON.ALLPROT"
3066       include "COMMON.IOUNITS"
3067       include "COMMON.GEO"
3068       include "COMMON.SBRIDGE"
3069       double precision thetnorm
3070       integer iprot,i,j,ii
3071       do i=1,nres_zs(iprot)
3072         do j=1,3
3073           c_zs(j,i,ii,iprot)=c(j,i)
3074         enddo
3075       enddo
3076       do i=nnt_zs(iprot),nct_zs(iprot)
3077         do j=1,3
3078           c_zs(j,i+nres,ii,iprot)=c(j,i+nres)
3079         enddo
3080       enddo
3081 c 5/7/02 AL: store sbridge info
3082       nss_zs(ii,iprot)=nss
3083       do i=1,nss
3084         ihpb_zs(i,ii,iprot)=ihpb(i)
3085         jhpb_zs(i,ii,iprot)=jhpb(i)
3086       enddo
3087       return
3088       end
3089 c------------------------------------------------------------------------------
3090       subroutine restore_ccoords(iprot,ii)
3091       implicit none
3092       include "DIMENSIONS"
3093       include "DIMENSIONS.ZSCOPT"
3094       include "COMMON.INTERACT"
3095       include "COMMON.VAR"
3096       include "COMMON.ALLPROT"
3097       include "COMMON.SBRIDGE"
3098       include "COMMON.CHAIN"
3099       include "COMMON.IOUNITS"
3100       integer iprot,i,j,ii
3101       do i=1,nres_zs(iprot)
3102         do j=1,3
3103           c(j,i)=c_zs(j,i,ii,iprot)
3104         enddo
3105       enddo
3106       do i=nnt_zs(iprot),nct_zs(iprot)
3107         do j=1,3
3108           c(j,i+nres)=c_zs(j,i+nres,ii,iprot)
3109         enddo
3110       enddo
3111 c 5/7/02 AL: restore sbridge info
3112       nss=nss_zs(ii,iprot)
3113       do i=1,nss
3114         ihpb(i)=ihpb_zs(i,ii,iprot)+nres
3115         jhpb(i)=jhpb_zs(i,ii,iprot)+nres
3116 c        forcon(i)=fbr
3117 c        dhpb(i)=dbr
3118       enddo
3119 #ifdef DEBUG
3120         write (iout,*) "restore_ccoords",ii
3121         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
3122         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
3123         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
3124         call flush(iout)
3125 #endif
3126       return
3127       end
3128 c------------------------------------------------------------------------------
3129       subroutine store_MDcoords(iprot,ii)
3130       implicit none
3131       include "DIMENSIONS"
3132       include "DIMENSIONS.ZSCOPT"
3133       include "COMMON.VAR"
3134       include "COMMON.CHAIN"
3135       include "COMMON.ALLPROT"
3136       include "COMMON.IOUNITS"
3137       include "COMMON.GEO"
3138       include "COMMON.SBRIDGE"
3139       include "COMMON.FMATCH"
3140       double precision thetnorm
3141       integer iprot,i,j,ii
3142       do i=1,nres_zs(iprot)
3143         do j=1,3
3144           cg_MD(j,i,ii,iprot)=c(j,i)
3145         enddo
3146       enddo
3147       do i=nnt_zs(iprot),nct_zs(iprot)
3148         do j=1,3
3149           cg_MD(j,i+nres,ii,iprot)=c(j,i+nres)
3150         enddo
3151       enddo
3152 c 5/7/02 AL: store sbridge info
3153       nss_MD(ii,iprot)=nss
3154       do i=1,nss
3155         ihpb_MD(i,ii,iprot)=ihpb(i)
3156         jhpb_MD(i,ii,iprot)=jhpb(i)
3157       enddo
3158       return
3159       end
3160 c------------------------------------------------------------------------------
3161       subroutine restore_MDcoords(iprot,ii)
3162       implicit none
3163       include "DIMENSIONS"
3164       include "DIMENSIONS.ZSCOPT"
3165       include "COMMON.INTERACT"
3166       include "COMMON.VAR"
3167       include "COMMON.ALLPROT"
3168       include "COMMON.SBRIDGE"
3169       include "COMMON.CHAIN"
3170       include "COMMON.IOUNITS"
3171       include "COMMON.FMATCH"
3172       integer iprot,i,j,ii
3173       do i=1,nres_zs(iprot)
3174         do j=1,3
3175           c(j,i)=cg_MD(j,i,ii,iprot)
3176         enddo
3177       enddo
3178       do i=nnt_zs(iprot),nct_zs(iprot)
3179         do j=1,3
3180           c(j,i+nres)=cg_MD(j,i+nres,ii,iprot)
3181         enddo
3182       enddo
3183 c 5/7/02 AL: restore sbridge info
3184       nss=nss_MD(ii,iprot)
3185       do i=1,nss
3186         ihpb(i)=ihpb_MD(i,ii,iprot)+nres
3187         jhpb(i)=jhpb_MD(i,ii,iprot)+nres
3188 c        forcon(i)=fbr
3189 c        dhpb(i)=dbr
3190       enddo
3191 #ifdef DEBUG
3192         write (iout,*) "restore_ccoords",ii
3193         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
3194         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
3195         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
3196         call flush(iout)
3197 #endif
3198       return
3199       end
3200 c------------------------------------------------------------------------------
3201       subroutine card_concat(card,to_upper)
3202       implicit none
3203       include 'DIMENSIONS.ZSCOPT'
3204       include "COMMON.IOUNITS"
3205       character*(*) card
3206       character*80 karta,ucase
3207       logical to_upper
3208       integer ilen
3209       external ilen
3210       read (inp,'(a)') karta
3211       if (to_upper) karta=ucase(karta)
3212       card=' '
3213       do while (karta(80:80).eq.'&')
3214         card=card(:ilen(card)+1)//karta(:79)
3215         read (inp,'(a)') karta
3216         if (to_upper) karta=ucase(karta)
3217       enddo
3218       card=card(:ilen(card)+1)//karta
3219       return
3220       end
3221 c------------------------------------------------------------------------------
3222       subroutine readi(rekord,lancuch,wartosc,default)
3223       implicit none
3224       character*(*) rekord,lancuch
3225       integer wartosc,default
3226       integer ilen,iread
3227       external ilen
3228       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
3229       if (iread.eq.0) then
3230         wartosc=default
3231         return
3232       endif
3233       iread=iread+ilen(lancuch)+1
3234       read (rekord(iread:),*) wartosc
3235       return
3236       end
3237 c----------------------------------------------------------------------------
3238       subroutine reada(rekord,lancuch,wartosc,default)
3239       implicit none
3240       character*(*) rekord,lancuch
3241       character*80 aux
3242       double precision wartosc,default
3243       integer ilen,iread
3244       external ilen
3245       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
3246       if (iread.eq.0) then
3247         wartosc=default
3248         return
3249       endif
3250       iread=iread+ilen(lancuch)+1
3251       read (rekord(iread:),*) wartosc
3252       return
3253       end
3254 c----------------------------------------------------------------------------
3255       subroutine multreadi(rekord,lancuch,tablica,dim,default)
3256       implicit none
3257       integer dim,i
3258       integer tablica(dim),default
3259       character*(*) rekord,lancuch
3260       character*80 aux
3261       integer ilen,iread
3262       external ilen
3263       do i=1,dim
3264         tablica(i)=default
3265       enddo
3266       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
3267       if (iread.eq.0) return
3268       iread=iread+ilen(lancuch)+1
3269       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
3270    10 return
3271       end
3272 c----------------------------------------------------------------------------
3273       subroutine multreada(rekord,lancuch,tablica,dim,default)
3274       implicit none
3275       integer dim,i
3276       double precision tablica(dim),default
3277       character*(*) rekord,lancuch
3278       character*80 aux
3279       integer ilen,iread
3280       external ilen
3281       do i=1,dim
3282         tablica(i)=default
3283       enddo
3284       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
3285       if (iread.eq.0) return
3286       iread=iread+ilen(lancuch)+1
3287       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
3288    10 return
3289       end
3290 c----------------------------------------------------------------------------
3291       subroutine reads(rekord,lancuch,wartosc,default)
3292       implicit none
3293       character*(*) rekord,lancuch,wartosc,default
3294       character*80 aux
3295       integer ilen,lenlan,lenrec,iread,ireade
3296       external ilen
3297       logical iblnk
3298       external iblnk
3299       lenlan=ilen(lancuch)
3300       lenrec=ilen(rekord)
3301       iread=index(rekord,lancuch(:lenlan)//"=")
3302 c      print *,"rekord",rekord," lancuch",lancuch
3303 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3304       if (iread.eq.0) then
3305         wartosc=default
3306         return
3307       endif
3308       iread=iread+lenlan+1
3309 c      print *,"iread",iread
3310 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3311       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3312         iread=iread+1
3313 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3314       enddo
3315 c      print *,"iread",iread
3316       if (iread.gt.lenrec) then
3317          wartosc=default
3318         return
3319       endif
3320       ireade=iread+1
3321 c      print *,"ireade",ireade
3322       do while (ireade.lt.lenrec .and.
3323      &   .not.iblnk(rekord(ireade:ireade)))
3324         ireade=ireade+1
3325       enddo
3326       wartosc=rekord(iread:ireade)
3327       return
3328       end
3329 c----------------------------------------------------------------------------
3330       subroutine multreads(rekord,lancuch,tablica,dim,default)
3331       implicit none
3332       integer dim,i
3333       character*(*) rekord,lancuch,tablica(dim),default
3334       character*80 aux
3335       integer ilen,lenlan,lenrec,iread,ireade
3336       external ilen
3337       logical iblnk
3338       external iblnk
3339       do i=1,dim
3340         tablica(i)=default
3341       enddo
3342       lenlan=ilen(lancuch)
3343       lenrec=ilen(rekord)
3344       iread=index(rekord,lancuch(:lenlan)//"=")
3345 c      print *,"rekord",rekord," lancuch",lancuch
3346 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3347       if (iread.eq.0) return
3348       iread=iread+lenlan+1
3349       do i=1,dim
3350 c      print *,"iread",iread
3351 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3352       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3353         iread=iread+1
3354 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3355       enddo
3356 c      print *,"iread",iread
3357       if (iread.gt.lenrec) return
3358       ireade=iread+1
3359 c      print *,"ireade",ireade
3360       do while (ireade.lt.lenrec .and.
3361      &   .not.iblnk(rekord(ireade:ireade)))
3362         ireade=ireade+1
3363       enddo
3364       tablica(i)=rekord(iread:ireade)
3365       iread=ireade+1
3366       enddo
3367       end
3368 c----------------------------------------------------------------------------
3369       subroutine split_string(rekord,tablica,dim,nsub)
3370       implicit none
3371       integer dim,nsub,i,ii,ll,kk
3372       character*(*) tablica(dim)
3373       character*(*) rekord
3374       integer ilen
3375       external ilen
3376       do i=1,dim
3377         tablica(i)=" "
3378       enddo
3379       ii=1
3380       ll = ilen(rekord)
3381       nsub=0
3382       do i=1,dim
3383 C Find the start of term name
3384         kk = 0
3385         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
3386           ii = ii+1
3387         enddo
3388 C Parse the name into TABLICA(i) until blank found
3389         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
3390           kk = kk+1
3391           tablica(i)(kk:kk)=rekord(ii:ii)
3392           ii = ii+1
3393         enddo
3394         if (kk.gt.0) nsub=nsub+1
3395         if (ii.gt.ll) return
3396       enddo
3397       return
3398       end
3399 c-------------------------------------------------------------------------------
3400       integer function iroof(n,m)
3401       ii = n/m
3402       if (ii*m .lt. n) ii=ii+1
3403       iroof = ii
3404       return
3405       end