source/maxlik/src-Fmatch/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 c      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,.false.)
 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,.false.)
 409               ityp_psc(2,npair_sc)=rescode(1,item(2),0,.false.)
 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,.false.)
 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,.false.))
 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,.false.))
 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,.false.))
 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,.false.))
 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,.false.)
 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,.false.)
 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,.false.)
 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,.false.)
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 (itype(j-1).eq.ntyp1 .or. itype(j).eq.ntyp1) cycle
2151         if (vbld(j).lt.2.0d0 .or. vbld(j).gt.6.5d0) then
2152           write (iout,*) "nnt",nnt," nct",nct
2153           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j)
2154           write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
2155           write (iout,*) "The Cartesian geometry is:"
2156           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
2157           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
2158           write (iout,*) "The internal geometry is:"
2159           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2160           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2161           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2162           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2163           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2164           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2165           write (iout,*)
2166      &      "This conformation WILL NOT be added to the database."
2167           return
2168         endif
2169       enddo
2170       do j=nnt,nct
2171         itj=itype(j)
2172         if (itype(j).ne.10 .and. itype(j).ne.ntyp1
2173      &    .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
2174           write (iout,*) "nnt",nnt," nct",nct
2175           write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
2176           write (iout,*) "The Cartesian geometry is:"
2177           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
2178           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
2179           write (iout,*) "The internal geometry is:"
2180           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2181           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2182           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2183           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2184           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2185           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2186           write (iout,*)
2187      &      "This conformation WILL NOT be added to the database."
2188           return
2189         endif
2190       enddo
2191       do j=3,nres
2192         if (itype(j).eq.ntyp1 .or. itype(j-1).eq.ntyp1
2193      &       .or. itype(j-1).eq.ntyp1) cycle
2194         if (theta(j).le.0.0d0) then
2195           write (iout,*)
2196      &      "Zero theta angle(s) in conformation",ii
2197           write (iout,*) "The Cartesian geometry is:"
2198           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
2199           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
2200           write (iout,*) "The internal geometry is:"
2201           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2202           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2203           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2204           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2205           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2206           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2207           write (iout,*)
2208      &      "This conformation WILL NOT be added to the database."
2209           return
2210         endif
2211         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
2212       enddo
2213       if (.not. init_ene) then
2214       etot=0.0d0
2215       do j=1,n_ene
2216         etot=etot+ww(j)*enetb(icount+1,j,iprot)
2217       enddo
2218       inan=0
2219 c detecting NaNQ
2220 #ifdef ISNAN
2221 #ifdef AIX
2222       if (isnan(etot).ne.0) inan=1
2223 #else
2224       if (isnan(etot)) inan=1
2225 #endif
2226 #else
2227 #ifdef WINPGI
2228       idumm=proc_proc(etot,inan)
2229 #else
2230       call proc_proc(etot,inan)
2231 #endif
2232 #endif
2233       else
2234       inan=0
2235       endif
2236       if (inan.eq.1) then
2237         write (iout,*) "NaNs detected in some of the energy",
2238      &    " components for protein",iprot," batch",ibatch,
2239      &    " conformation",jjj
2240         write (iout,*) "etot",etot
2241         write (iout,*) "The Cartesian geometry is:"
2242         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2243         write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
2244         write (iout,*) "The internal geometry is:"
2245         write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2246         write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2247         write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2248         write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2249         write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
2250         write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2251         write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2252         write (iout,*) "The components of the energy are:"
2253         energia(0)=etot
2254         do k=1,n_ene
2255           energia(k)=enetb(jj+1,k,iprot)
2256         enddo
2257         call enerprint(energia(0))
2258         write (iout,*)
2259      &    "This conformation WILL NOT be added to the database."
2260         call flush(iout)
2261       else
2262         jj=jj+1
2263 #ifdef DEBUG
2264         write (iout,'(e15.5,16i5)') entfac(icount+1,iprot),
2265      &        iscore(icount+1,0,iprot),ibatch
2266         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2267         write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
2268         write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
2269         write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
2270         write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
2271         write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
2272         write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
2273         write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
2274         write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
2275         write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
2276         write (iout,'(5e15.5)') (enetb(icount+1,j,iprot),j=1,n_ene)
2277 c        write (iout,'(2e15.5)') (eneps(1,j,icount+1,iprot),
2278 c     &        eneps(2,j,icount+1,iprot),j=1,nntyp)
2279 #endif
2280 c        write (iout,*) "First nu in readrtms"
2281         call flush(iout)
2282         icount=icount+1
2283         list_conf(jj,iprot)=jj
2284         call store_cconf_from_file(jj,icount,iprot)
2285         do j=1,natlike(iprot)
2286           do k=1,natlik
2287             if (k.eq.numnat(j,iprot)) then
2288               do l=1,natdim(j,iprot)
2289                 nu(l,k,j,iprot)=nutemp(l,k)
2290               enddo
2291               exit
2292             endif
2293           enddo
2294         enddo
2295 c        write (iout,*) "End First nu in readrtms"
2296         call flush(iout)
2297         if (icount.eq.maxstr_proc) then
2298 #ifdef DEBUG
2299           write (iout,* ) "jj_old",jj_old," jj",jj
2300           write (iout,*) "Calling write_and_send_cconf"
2301           call flush(iout)
2302 #endif
2303           call write_and_send_cconf(icount,jj_old,jj,iprot,Next)
2304           jj_old=jj+1
2305 #ifdef DEBUG
2306           write (iout,*) "Exited write_and_send_cconf"
2307           call flush(iout)
2308 #endif
2309           icount=0
2310         endif
2311       endif
2312       return
2313       end
2314 c------------------------------------------------------------------------------
2315       subroutine store_cconf_from_file(jj,icount,iprot)
2316       implicit none
2317       include "DIMENSIONS"
2318       include "DIMENSIONS.ZSCOPT"
2319       include "COMMON.CHAIN"
2320       include "COMMON.SBRIDGE"
2321       include "COMMON.INTERACT"
2322       include "COMMON.IOUNITS"
2323       include "COMMON.CLASSES"
2324       include "COMMON.ALLPROT"
2325       include "COMMON.VAR"
2326       integer i,j,jj,icount,ibatch,iprot
2327 c Store the conformation that has been read in
2328       do i=1,2*nres
2329         do j=1,3
2330           c_zs(j,i,icount,iprot)=c(j,i)
2331         enddo
2332       enddo
2333       nss_zs(icount,iprot)=nss
2334       do i=1,nss
2335         ihpb_zs(i,icount,iprot)=ihpb(i)
2336         jhpb_zs(i,icount,iprot)=jhpb(i)
2337       enddo
2338       return
2339       end
2340 c------------------------------------------------------------------------------
2341       subroutine write_and_send_cconf(icount,jj_old,jj,iprot,Next)
2342       implicit none
2343       include "DIMENSIONS"
2344       include "DIMENSIONS.ZSCOPT"
2345 #ifdef MPI
2346       include "mpif.h"
2347       integer IERROR
2348       include "COMMON.MPI"
2349 #endif
2350       include "COMMON.WEIGHTS"
2351       include "COMMON.CHAIN"
2352       include "COMMON.SBRIDGE"
2353       include "COMMON.INTERACT"
2354       include "COMMON.IOUNITS"
2355       include "COMMON.CLASSES"
2356       include "COMMON.VAR"
2357       include "COMMON.ALLPROT"
2358       include "COMMON.ENERGIES"
2359       include "COMMON.WEIGHTDER"
2360       include "COMMON.OPTIM"
2361       include "COMMON.COMPAR"
2362       integer icount,jj_old,jj,Next,iprot
2363 c Write the structures to a scratch file
2364 #ifdef MPI
2365 c Master sends the portion of conformations that have been read in to the neighbor
2366 #ifdef DEBUG
2367       write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
2368       call flush(iout)
2369 #endif
2370       call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM,IERROR)
2371 #ifdef DEBUG
2372       write (iout,*) "Processor",me," Next",next," sent icount=",icount
2373       write (iout,*) "Processor",me," jj_old",jj_old," jj",jj
2374       call flush(iout)
2375 #endif
2376       if (icount.gt.0) then
2377       call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER,
2378      &    Next,571,WHAM_COMM,IERROR)
2379       call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2380      &    Next,572,WHAM_COMM,IERROR)
2381       if (.not. init_ene) call MPI_Send(enetb(jj_old,1,iprot),
2382      &    maxstr*n_ene,
2383      &    MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,IERROR)
2384       call MPI_Send(nu(1,1,jj_old,iprot),
2385      &    maxdimnat*natlike(iprot)*icount,
2386      &    MPI_DOUBLE_PRECISION,
2387      &    Next,577,WHAM_COMM,IERROR)
2388       call MPI_Send(eini(jj_old,iprot),icount,
2389      &    MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR)
2390       call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2391      &    Next,579,WHAM_COMM,IERROR)
2392       call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2,
2393      &    MPI_REAL,Next,580,WHAM_COMM,IERROR)
2394       if (.not. init_ene) call MPI_Send(eneps(1,1,1,iprot),
2395      &    2*icount*nntyp,
2396      &    MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR)
2397       endif
2398 #endif
2399       call flush(iout)
2400       call dawrite_ccoords(iprot,jj_old,jj,ientout)
2401 c Change AL 20/12/2017
2402       if (.not. mod_other_params)
2403      &call dawrite_ene(iprot,jj_old,jj,istat)
2404       return
2405       end
2406 c------------------------------------------------------------------------------
2407 #ifdef MPI
2408       subroutine receive_and_pass_cconf(icount,jj_old,jj,iprot,Previous,
2409      &  Next)
2410       implicit none
2411       include "DIMENSIONS"
2412       include "DIMENSIONS.ZSCOPT"
2413       include "mpif.h"
2414       integer IERROR,STATUS(MPI_STATUS_SIZE)
2415       include "COMMON.MPI"
2416       include "COMMON.CHAIN"
2417       include "COMMON.SBRIDGE"
2418       include "COMMON.INTERACT"
2419       include "COMMON.IOUNITS"
2420       include "COMMON.CLASSES"
2421       include "COMMON.ALLPROT"
2422       include "COMMON.ENERGIES"
2423       include "COMMON.VAR"
2424       include "COMMON.GEO"
2425       include "COMMON.WEIGHTS"
2426       include "COMMON.WEIGHTDER"
2427       include "COMMON.COMPAR"
2428       include "COMMON.OPTIM"
2429       integer i,j,k,icount,jj_old,jj,iprot,Previous,Next
2430       icount=1
2431 #ifdef DEBUG
2432       write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
2433       call flush(iout)
2434 #endif
2435       do while (icount.gt.0)
2436 c      call MPI_Probe(Previous,570,WHAM_COMM,STATUS,IERROR)
2437       call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,WHAM_COMM,
2438      &     STATUS,IERROR)
2439 #ifdef DEBUG
2440       write (iout,*)"Processor",me," previous",previous," icount",icount
2441       call flush(iout)
2442 #endif
2443       call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM,
2444      &     IERROR)
2445 #ifdef DEBUG
2446       write (iout,*) "Processor",me," icount",icount
2447       call flush(iout)
2448 #endif
2449       if (icount.eq.0) return
2450       call MPI_Recv(nss_zs(1,iprot),icount,MPI_INTEGER,
2451      &    Previous,571,WHAM_COMM,STATUS,IERROR)
2452       call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER,
2453      &  Next,571,WHAM_COMM,IERROR)
2454       call MPI_Recv(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2455      &    Previous,572,WHAM_COMM,STATUS,IERROR)
2456       call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2457      &  Next,572,WHAM_COMM,IERROR)
2458       if (.not. init_ene) then
2459         call MPI_Recv(enetb(jj_old,1,iprot),maxstr*n_ene,
2460      &   MPI_DOUBLE_PRECISION,Previous,576,WHAM_COMM,STATUS,IERROR)
2461         call MPI_Send(enetb(jj_old,1,iprot),maxstr*n_ene,
2462      &   MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,STATUS,IERROR)
2463       endif
2464       call MPI_Recv(nu(1,1,jj_old,iprot),
2465      &  maxdimnat*natlike(iprot)*icount,
2466      &  MPI_DOUBLE_PRECISION,
2467      &  Previous,577,WHAM_COMM,STATUS,IERROR)
2468       call MPI_Send(nu(1,1,jj_old,iprot),
2469      &  maxdimnat*natlike(iprot)*icount,
2470      &  MPI_DOUBLE_PRECISION,
2471      &  Next,577,WHAM_COMM,IERROR)
2472       call MPI_Recv(eini(jj_old,iprot),icount,
2473      &  MPI_DOUBLE_PRECISION,Previous,578,WHAM_COMM,STATUS,IERROR)
2474       call MPI_Send(eini(jj_old,iprot),icount,
2475      &  MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR)
2476       call MPI_Recv(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2477      &  Previous,579,WHAM_COMM,STATUS,IERROR)
2478       call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2479      &  Next,579,WHAM_COMM,IERROR)
2480       call MPI_Recv(c_zs(1,1,1,iprot),3*icount*maxres2,
2481      &  MPI_REAL,Previous,580,WHAM_COMM,STATUS,IERROR)
2482       call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2,
2483      &  MPI_REAL,Next,580,WHAM_COMM,IERROR)
2484       if (.not. init_ene) then
2485         call MPI_Recv(eneps(1,1,1,iprot),2*icount*nntyp,
2486      &   MPI_DOUBLE_PRECISION,Previous,582,WHAM_COMM,STATUS,IERROR)
2487         call MPI_Send(eneps(1,1,1,iprot),2*icount*nntyp,
2488      &   MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR)
2489       endif
2490       jj=jj_old+icount-1
2491       do i=jj_old,jj
2492         list_conf(i,iprot)=i
2493       enddo
2494       call dawrite_ccoords(iprot,jj_old,jj,ientout)
2495 c Change AL 12/20/2017
2496       if (.not. mod_other_params)
2497      &call dawrite_ene(iprot,jj_old,jj,istat)
2498       jj_old=jj+1
2499 #ifdef DEBUG
2500       write (iout,*) "Protein",iprot
2501       write (iout,*) "Processor",me," received",icount," conformations"
2502       do i=1,icount
2503         write (iout,'(8f10.4)') ((c_zs(l,k,i,iprot),l=1,3,k=1,nres)
2504         write (iout,'(8f10.4)')((c_zs(l,k,i+nres,iprot),l=1,3,k=nnt,nct)
2505         write (iout,'(16i5)') nss_zs(i,iprot),(ihpb_zs(k,i,iprot),
2506      &    jhpb_zs(k,i,iprot),k=1,nss_zs(i,iprot))
2507         write (iout,'(5e15.5)') (enetb(i,j,iprot),j=1,n_ene)
2508         write (iout,'(2e15.5)') (eneps(1,j,i,iprot),
2509      &        eneps(2,j,i,iprot),j=1,nntyp)
2510         write (iout,'(e15.5,16i5)') entfac(i),iscore(i,0,iprot),
2511      &        kbatch(i,iprot)
2512       enddo
2513 #endif
2514       enddo
2515       return
2516       end
2517 #endif
2518 c------------------------------------------------------------------------------
2519       subroutine read_thermal
2520       implicit none
2521       include "DIMENSIONS"
2522       include "DIMENSIONS.ZSCOPT"
2523       include "COMMON.CHAIN"
2524       include "COMMON.SBRIDGE"
2525       include "COMMON.INTERACT"
2526       include "COMMON.IOUNITS"
2527       include "COMMON.CLASSES"
2528       include "COMMON.VAR"
2529       include "COMMON.THERMAL"
2530       character*800 controlcard
2531       call card_concat(controlcard,.true.)
2532       call readi(controlcard,"NGRIDT",NGridT,200)
2533       call reada(controlcard,"DELTAT",deltaT,5.0d0)
2534       call reada(controlcard,"T0",GridT0,2.0d2)
2535       write (iout,*) "Parameters of thermal curves"
2536       write (iout,*) "NGRIDT",NGridT," DELTAT",deltaT," T0",GridT0
2537       return
2538       end
2539 c------------------------------------------------------------------------------
2540       subroutine daread_ene(iprot,istart_conf,iend_conf)
2541       implicit none
2542       include "DIMENSIONS"
2543       include "DIMENSIONS.ZSCOPT"
2544 #ifdef MPI
2545       include "mpif.h"
2546       include "COMMON.MPI"
2547 #endif
2548       include "COMMON.CHAIN"
2549       include "COMMON.CLASSES"
2550       include "COMMON.ENERGIES"
2551       include "COMMON.IOUNITS"
2552       include "COMMON.PROTFILES"
2553       include "COMMON.ALLPROT"
2554       include "COMMON.WEIGHTDER"
2555       include "COMMON.COMPAR"
2556       include "COMMON.VMCPAR"
2557       integer iprot,istart_conf,iend_conf
2558       integer i,ii,iii,j,k
2559 #ifdef DEBUG
2560       write (iout,*) "Calling DAREAD_ENE"
2561 #endif
2562 c      write (iout,*) "Reading: n_ene",n_ene
2563 c      call flush(iout)
2564 c      write (iout,*) "DAREAD_ENE",istart_conf,iend_conf
2565 c
2566 c Read conformations off a DA scratchfile.
2567 c
2568       do ii=istart_conf,iend_conf
2569         iii=list_conf(ii,iprot)
2570         i = ii - istart_conf + 1
2571         read(ientin,rec=iii) (enetb(i,j,iprot),j=1,n_ene),
2572      &     (enetb_orig(i,j,iprot),j=1,n_ene),
2573      &     (enetb_oorig(i,j,iprot),j=1,n_ene),
2574      &     eini(ii,iprot),entfac(ii,iprot),
2575      &     (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp),
2576      &     ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2577 #ifdef DEBUG
2578         write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot),
2579      &  entfac(ii,iprot),
2580         write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene)
2581         write (iout,'(18(1pe12.4))')
2582      &  ((nu(k,j,i,iprot)k=1,maxdimnat),j=1,natlike(iprot))
2583         call flush(iout)
2584 #endif
2585       enddo
2586
2587       return
2588       end
2589 c------------------------------------------------------------------------------
2590       subroutine dawrite_ene(iprot,istart_conf,iend_conf,unit_out)
2591       implicit none
2592       include "DIMENSIONS"
2593       include "DIMENSIONS.ZSCOPT"
2594 #ifdef MPI
2595       include "mpif.h"
2596       include "COMMON.MPI"
2597 #endif
2598       include "COMMON.CHAIN"
2599       include "COMMON.CLASSES"
2600       include "COMMON.ENERGIES"
2601       include "COMMON.IOUNITS"
2602       include "COMMON.PROTFILES"
2603       include "COMMON.ALLPROT"
2604       include "COMMON.WEIGHTDER"
2605       include "COMMON.VMCPAR"
2606       include "COMMON.COMPAR"
2607       integer iprot,istart_conf,iend_conf,unit_out
2608       integer i,ii,iii,j,k
2609 c      write (iout,*) "Writing: n_ene",n_ene
2610 c      call flush(iout)
2611 c      write (iout,*) "DAWRITE_ENE",istart_conf,iend_conf
2612 c
2613 c Write conformations to a DA scratchfile.
2614 c
2615       do ii=istart_conf,iend_conf
2616         iii=list_conf(ii,iprot)
2617         i = ii - istart_conf + 1
2618         write(unit_out,rec=iii) (enetb(i,j,iprot),j=1,n_ene),
2619      &     (enetb_orig(i,j,iprot),j=1,n_ene),
2620      &     (enetb_oorig(i,j,iprot),j=1,n_ene),
2621      &     eini(ii,iprot),entfac(ii,iprot),
2622      &     (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp),
2623      &     ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2624 #ifdef DEBUG
2625         write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot),
2626      &  entfac(ii,iprot)
2627         write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene)
2628         write (iout,'(18(1pe12.4))')
2629      &  ((nu(kj,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2630         call flush(iout)
2631 #endif
2632       enddo
2633
2634       return
2635       end
2636 c------------------------------------------------------------------------------
2637       subroutine daread_forces(iprot,istart_conf,iend_conf)
2638       implicit none
2639       include "DIMENSIONS"
2640       include "DIMENSIONS.ZSCOPT"
2641 #ifdef MPI
2642       include "mpif.h"
2643       include "COMMON.MPI"
2644 #endif
2645       include "COMMON.CHAIN"
2646       include "COMMON.CLASSES"
2647       include "COMMON.ENERGIES"
2648       include "COMMON.IOUNITS"
2649       include "COMMON.PROTFILES"
2650       include "COMMON.ALLPROT"
2651       include "COMMON.WEIGHTDER"
2652       include "COMMON.COMPAR"
2653       include "COMMON.VMCPAR"
2654       include "COMMON.FMATCH"
2655       include "COMMON.NAMES"
2656       include "COMMON.INTERACT"
2657       integer iprot,istart_conf,iend_conf
2658       integer i,ii,iii,j,jj,k,l
2659       integer kiri
2660       character*16 forma,acc
2661       character*128 nam
2662 #ifdef DEBUG
2663       write (iout,*) "Calling DAREAD_FORCES"
2664 #endif
2665 c      write (iout,*) "Reading: n_ene",n_ene
2666 c      call flush(iout)
2667 c      write (iout,*) "DAREAD_ENE",istart_conf,iend_conf
2668 c
2669 c Read conformations off a DA scratchfile.
2670 c
2671 #ifdef DEBUG
2672       inquire(unit=ientin,name=nam,recl=kiri,form=forma,access=acc)
2673       write (iout,*) "DAREAD_FORCES","len=",kiri," form=",forma,
2674      &  " acc=",acc
2675       write (iout,*) "nam=",nam
2676 #endif
2677       do ii=istart_conf,iend_conf
2678         iii=list_conf_MD(ii,iprot)
2679         i = ii - istart_conf + 1
2680         read(ientin,rec=iii) (((forctb(j,l,k,i,iprot),j=1,n_ene),l=1,3),
2681      &     k=1,nsite(iprot))
2682 #ifdef DEBUG
2683         write (iout,'(a,4i5)') "DAREAD_FORCES",ii,iii,i,nsite(iprot)
2684         write (iout,*) "CG force components"
2685         do k=1,n_ene
2686
2687           write (iout,'(a,i3,1x,a)') "Component",k,ename(k)
2688
2689           jj=nres
2690           do j=1,nres
2691             if (itype(j).eq.10 .or. itype(j).eq.ntyp1) then
2692               write (iout,'(a4,i5,3e15.5)')
2693      &        restyp(itype(j)),j,(forctb(k,l,j,i,iprot),l=1,3)
2694             else
2695               jj=jj+1
2696               write (iout,'(a4,i5,3e15.5,5x,3e15.5)')
2697      &        restyp(itype(j)),j,(forctb(k,l,j,i,iprot),l=1,3),
2698      &        (forctb(k,l,jj,i,iprot),l=1,3)
2699             endif
2700           enddo
2701
2702         enddo
2703         call flush(iout)
2704 #endif
2705       enddo
2706
2707       return
2708       end
2709 c------------------------------------------------------------------------------
2710       subroutine dawrite_forces(iprot,istart_conf,iend_conf,unit_out)
2711       implicit none
2712       include "DIMENSIONS"
2713       include "DIMENSIONS.ZSCOPT"
2714 #ifdef MPI
2715       include "mpif.h"
2716       include "COMMON.MPI"
2717 #endif
2718       include "COMMON.CHAIN"
2719       include "COMMON.CLASSES"
2720       include "COMMON.ENERGIES"
2721       include "COMMON.IOUNITS"
2722       include "COMMON.PROTFILES"
2723       include "COMMON.ALLPROT"
2724       include "COMMON.WEIGHTDER"
2725       include "COMMON.VMCPAR"
2726       include "COMMON.COMPAR"
2727       include "COMMON.FMATCH"
2728       include "COMMON.NAMES"
2729       include "COMMON.INTERACT"
2730       integer iprot,istart_conf,iend_conf,unit_out
2731       integer i,ii,iii,j,jj,k,l
2732       integer kiri
2733       character*16 forma,acc
2734       character*128 nam
2735 c      write (iout,*) "Writing: n_ene",n_ene
2736 c      call flush(iout)
2737 c      write (iout,*) "DAWRITE_ENE",istart_conf,iend_conf
2738 c
2739 c Write conformations to a DA scratchfile.
2740 c
2741 #ifdef DEBUG
2742       inquire(unit=unit_out,name=nam,recl=kiri,form=forma,access=acc)
2743       write (iout,*) "DAWRITE_FORCES","len=",kiri," form=",forma,
2744      &  " acc=",acc
2745       write (iout,*) "nam=",nam
2746       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2747       write (iout,*) "list",(list_conf_MD(k,iprot),
2748      & k=istart_conf,iend_conf)
2749 #endif
2750       do ii=istart_conf,iend_conf
2751         iii=list_conf_MD(ii,iprot)
2752         i = ii - istart_conf + 1
2753         write(unit_out,rec=iii) (((forctb(j,l,k,i,iprot),j=1,n_ene),
2754      &     l=1,3),k=1,nsite(iprot))
2755 #ifdef DEBUG
2756         write (iout,'(a,4i5)') "DAWRITE_FORCES",ii,iii,i,nsite(iprot)
2757         write (iout,*) "CG force components"
2758         do k=1,n_ene
2759
2760           write (iout,'(a,i3,1x,a)') "Component",k,ename(k)
2761
2762           jj=nres
2763           do j=1,nres
2764             if (itype(j).eq.10 .or. itype(j).eq.ntyp1) then
2765               write (iout,'(a4,i5,3e15.5)')
2766      &        restyp(itype(j)),j,(forctb(k,l,j,i,iprot),l=1,3)
2767             else
2768               jj=jj+1
2769               write (iout,'(a4,i5,3e15.5,5x,3e15.5)')
2770      &        restyp(itype(j)),j,(forctb(k,l,j,i,iprot),l=1,3),
2771      &        (forctb(k,l,jj,i,iprot),l=1,3)
2772             endif
2773           enddo
2774
2775         enddo
2776         call flush(iout)
2777 #endif
2778       enddo
2779
2780       return
2781       end
2782 c------------------------------------------------------------------------------
2783       subroutine daread_ccoords(iprot,istart_conf,iend_conf)
2784       implicit none
2785       include "DIMENSIONS"
2786       include "DIMENSIONS.ZSCOPT"
2787 #ifdef MPI
2788       include "mpif.h"
2789       include "COMMON.MPI"
2790 #endif
2791       include "COMMON.CHAIN"
2792       include "COMMON.CLASSES"
2793       include "COMMON.ENERGIES"
2794       include "COMMON.IOUNITS"
2795       include "COMMON.PROTFILES"
2796       include "COMMON.ALLPROT"
2797       include "COMMON.INTERACT"
2798       include "COMMON.VAR"
2799       include "COMMON.SBRIDGE"
2800       include "COMMON.GEO"
2801       include "COMMON.COMPAR"
2802       include "COMMON.VMCPAR"
2803       include "COMMON.WEIGHTDER"
2804       integer iprot,istart_conf,iend_conf
2805       integer i,j,k,ij,ii,iii
2806       integer len
2807       real*4 csingle(3,maxres2)
2808       character*16 form,acc
2809       character*32 nam
2810 c
2811 c Read conformations off a DA scratchfile.
2812 c
2813 #ifdef DEBUG
2814       write (iout,*) "DAREAD_COORDS"
2815       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2816       write (iout,*) "lenrec",lenrec(iprot)
2817       inquire(unit=ientin,name=nam,recl=len,form=form,access=acc)
2818       write (iout,*) "len=",len," form=",form," acc=",acc
2819       write (iout,*) "nam=",nam
2820       call flush(iout)
2821 #endif
2822       do ii=istart_conf,iend_conf
2823         iii=list_conf(ii,iprot)
2824         ij = ii - istart_conf + 1
2825 #ifdef DEBUG
2826         write (iout,*) "Reading binary file, record",iii," ii",ii
2827         call flush(iout)
2828 #endif
2829         read(icbase,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2830      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2831      &    nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot),
2832      &    entfac(ii,iprot),
2833      &    ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot))
2834         do i=1,2*nres
2835           do j=1,3
2836             c(j,i)=csingle(j,i)
2837           enddo
2838         enddo
2839 #ifdef DEBUG
2840         write (iout,*) "iprot",iprot," ii",ii
2841         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2842         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2843         write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot)
2844         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2845         call flush(iout)
2846 #endif
2847         call store_ccoords(iprot,ii-istart_conf+1)
2848       enddo
2849       return
2850       end
2851 c------------------------------------------------------------------------------
2852       subroutine dawrite_ccoords(iprot,istart_conf,iend_conf,unit_out)
2853       implicit none
2854       include "DIMENSIONS"
2855       include "DIMENSIONS.ZSCOPT"
2856 #ifdef MPI
2857       include "mpif.h"
2858       include "COMMON.MPI"
2859 #endif
2860       include "COMMON.CHAIN"
2861       include "COMMON.INTERACT"
2862       include "COMMON.CLASSES"
2863       include "COMMON.ENERGIES"
2864       include "COMMON.IOUNITS"
2865       include "COMMON.PROTFILES"
2866       include "COMMON.ALLPROT"
2867       include "COMMON.VAR"
2868       include "COMMON.SBRIDGE"
2869       include "COMMON.GEO"
2870       include "COMMON.COMPAR"
2871       include "COMMON.WEIGHTDER"
2872       include "COMMON.VMCPAR"
2873       real*4 csingle(3,maxres2)
2874       integer iprot,istart_conf,iend_conf
2875       integer i,j,k,ii,ij,iii,unit_out
2876       integer len
2877       character*16 form,acc
2878       character*32 nam
2879 c
2880 c Write conformations to a DA scratchfile.
2881 c
2882 #ifdef DEBUG
2883       write (iout,*) "DAWRITE_COORDS"
2884       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2885       write (iout,*) "lenrec",lenrec(iprot)
2886       inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
2887       write (iout,*) "len=",len," form=",form," acc=",acc
2888       write (iout,*) "nam=",nam
2889       call flush(iout)
2890 #endif
2891       do ii=istart_conf,iend_conf
2892         iii=list_conf(ii,iprot)
2893         ij = ii - istart_conf + 1
2894         call restore_ccoords(iprot,ii-istart_conf+1)
2895 #ifdef DEBUG
2896         write (iout,*) "Writing binary file, record",iii," ii",ii
2897         call flush(iout)
2898 #endif
2899         do i=1,2*nres
2900           do j=1,3
2901             csingle(j,i)=c(j,i)
2902           enddo
2903         enddo
2904         write(unit_out,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2905      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2906      &    nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot),
2907      &    entfac(ii,iprot),
2908      &    ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot))
2909 #ifdef DEBUG
2910         write (iout,*) "kbatch",kbatch(ii,iprot)," iscore",
2911      &    iscore(ii,0,iprot)
2912         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2913         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2914         write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot)
2915         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2916         call flush(iout)
2917 #endif
2918       enddo
2919       return
2920       end
2921 c------------------------------------------------------------------------------
2922       subroutine daread_MDcoords(iprot,istart_conf,iend_conf)
2923       implicit none
2924       include "DIMENSIONS"
2925       include "DIMENSIONS.ZSCOPT"
2926 #ifdef MPI
2927       include "mpif.h"
2928       include "COMMON.MPI"
2929 #endif
2930       include "COMMON.CHAIN"
2931       include "COMMON.CLASSES"
2932       include "COMMON.ENERGIES"
2933       include "COMMON.IOUNITS"
2934       include "COMMON.PROTFILES"
2935       include "COMMON.ALLPROT"
2936       include "COMMON.INTERACT"
2937       include "COMMON.VAR"
2938       include "COMMON.SBRIDGE"
2939       include "COMMON.GEO"
2940       include "COMMON.COMPAR"
2941       include "COMMON.VMCPAR"
2942       include "COMMON.WEIGHTDER"
2943       include "COMMON.FMATCH"
2944       integer iprot,istart_conf,iend_conf
2945       integer i,j,k,ij,ii,iii
2946       integer len
2947       real*4 csingle(3,maxres2),fsingle(3,maxres2)
2948       character*16 form,acc
2949       character*32 nam
2950 c
2951 c Read conformations off a DA scratchfile.
2952 c
2953 #ifdef DEBUG
2954       write (iout,*) "DAREAD_MDCOORDS"
2955       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2956       write (iout,*) "lenrec_MD",lenrec_MD(iprot)
2957       inquire(unit=ientin,name=nam,recl=len,form=form,access=acc)
2958       write (iout,*) "len=",len," form=",form," acc=",acc
2959       write (iout,*) "nam=",nam
2960       call flush(iout)
2961 #endif
2962       do ii=istart_conf,iend_conf
2963         iii=list_conf_MD(ii,iprot)
2964         ij = ii - istart_conf + 1
2965 #ifdef DEBUG
2966         write (iout,*) "Reading binary file, record",iii," ii",ii
2967         call flush(iout)
2968 #endif
2969         read(icbase,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2970      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2971      &    ((CGforce_MD(j,i,ii,iprot),j=1,3),i=1,nsite(iprot)),
2972      &    nss,(ihpb(i),jhpb(i),i=1,nss)
2973         do i=1,2*nres
2974           do j=1,3
2975             c(j,i)=csingle(j,i)
2976           enddo
2977         enddo
2978 #ifdef DEBUG
2979         write (iout,*) "iprot",iprot," ii",ii
2980         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2981         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2982         write (iout,'(8f10.5)') ((CGforce_MD(j,i,ii,iprot),j=1,3),i=1,
2983      &    nsite(iprot))
2984         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2985         call flush(iout)
2986 #endif
2987         call store_MDcoords(iprot,ii-istart_conf+1)
2988       enddo
2989       return
2990       end
2991 c------------------------------------------------------------------------------
2992       subroutine dawrite_MDcoords(iprot,istart_conf,iend_conf,unit_out)
2993       implicit none
2994       include "DIMENSIONS"
2995       include "DIMENSIONS.ZSCOPT"
2996 #ifdef MPI
2997       include "mpif.h"
2998       include "COMMON.MPI"
2999 #endif
3000       include "COMMON.CHAIN"
3001       include "COMMON.INTERACT"
3002       include "COMMON.CLASSES"
3003       include "COMMON.ENERGIES"
3004       include "COMMON.IOUNITS"
3005       include "COMMON.PROTFILES"
3006       include "COMMON.ALLPROT"
3007       include "COMMON.VAR"
3008       include "COMMON.SBRIDGE"
3009       include "COMMON.GEO"
3010       include "COMMON.COMPAR"
3011       include "COMMON.WEIGHTDER"
3012       include "COMMON.VMCPAR"
3013       include "COMMON.FMATCH"
3014       real*4 csingle(3,maxres2)
3015       integer iprot,istart_conf,iend_conf
3016       integer i,j,k,ii,ij,iii,unit_out
3017       integer len
3018       character*16 form,acc
3019       character*32 nam
3020 c
3021 c Write conformations to a DA scratchfile.
3022 c
3023 #ifdef DEBUG
3024       write (iout,*) "DAWRITE_MDCOORDS"
3025       write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
3026       write (iout,*) "lenrec_MD",lenrec_MD(iprot)
3027       inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
3028       write (iout,*) "len=",len," form=",form," acc=",acc
3029       write (iout,*) "nam=",nam
3030       call flush(iout)
3031 #endif
3032       do ii=istart_conf,iend_conf
3033         iii=list_conf_MD(ii,iprot)
3034         ij = iii - istart_conf + 1
3035         call restore_MDcoords(iprot,ii-istart_conf+1)
3036 #ifdef DEBUG
3037         write (iout,*) "Writing binary file, record",iii," ii",ii,
3038      &  " ij",ij
3039         write (iout,*) "nss",nss
3040         call flush(iout)
3041 #endif
3042         do i=1,2*nres
3043           do j=1,3
3044             csingle(j,i)=c(j,i)
3045           enddo
3046         enddo
3047         write(unit_out,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
3048      &    ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
3049      &    ((CGforce_MD(j,i,ij,iprot),j=1,3),i=1,nsite(iprot)),
3050      &    nss,(ihpb(i),jhpb(i),i=1,nss)
3051 #ifdef DEBUG
3052         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
3053         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
3054         write (iout,'(8f10.5)') ((CGforce_MD(j,i,ij,iprot),j=1,3),
3055      &    i=1,nres)
3056         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
3057         call flush(iout)
3058 #endif
3059       enddo
3060       return
3061       end
3062 c------------------------------------------------------------------------------
3063       subroutine store_ccoords(iprot,ii)
3064       implicit none
3065       include "DIMENSIONS"
3066       include "DIMENSIONS.ZSCOPT"
3067       include "COMMON.VAR"
3068       include "COMMON.CHAIN"
3069       include "COMMON.ALLPROT"
3070       include "COMMON.IOUNITS"
3071       include "COMMON.GEO"
3072       include "COMMON.SBRIDGE"
3073       double precision thetnorm
3074       integer iprot,i,j,ii
3075       do i=1,nres_zs(iprot)
3076         do j=1,3
3077           c_zs(j,i,ii,iprot)=c(j,i)
3078         enddo
3079       enddo
3080       do i=nnt_zs(iprot),nct_zs(iprot)
3081         do j=1,3
3082           c_zs(j,i+nres,ii,iprot)=c(j,i+nres)
3083         enddo
3084       enddo
3085 c 5/7/02 AL: store sbridge info
3086       nss_zs(ii,iprot)=nss
3087       do i=1,nss
3088         ihpb_zs(i,ii,iprot)=ihpb(i)
3089         jhpb_zs(i,ii,iprot)=jhpb(i)
3090       enddo
3091       return
3092       end
3093 c------------------------------------------------------------------------------
3094       subroutine restore_ccoords(iprot,ii)
3095       implicit none
3096       include "DIMENSIONS"
3097       include "DIMENSIONS.ZSCOPT"
3098       include "COMMON.INTERACT"
3099       include "COMMON.VAR"
3100       include "COMMON.ALLPROT"
3101       include "COMMON.SBRIDGE"
3102       include "COMMON.CHAIN"
3103       include "COMMON.IOUNITS"
3104       integer iprot,i,j,ii
3105       do i=1,nres_zs(iprot)
3106         do j=1,3
3107           c(j,i)=c_zs(j,i,ii,iprot)
3108         enddo
3109       enddo
3110       do i=nnt_zs(iprot),nct_zs(iprot)
3111         do j=1,3
3112           c(j,i+nres)=c_zs(j,i+nres,ii,iprot)
3113         enddo
3114       enddo
3115 c 5/7/02 AL: restore sbridge info
3116       nss=nss_zs(ii,iprot)
3117       do i=1,nss
3118         ihpb(i)=ihpb_zs(i,ii,iprot)+nres
3119         jhpb(i)=jhpb_zs(i,ii,iprot)+nres
3120 c        forcon(i)=fbr
3121 c        dhpb(i)=dbr
3122       enddo
3123 #ifdef DEBUG
3124         write (iout,*) "restore_ccoords",ii
3125         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
3126         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
3127         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
3128         call flush(iout)
3129 #endif
3130       return
3131       end
3132 c------------------------------------------------------------------------------
3133       subroutine store_MDcoords(iprot,ii)
3134       implicit none
3135       include "DIMENSIONS"
3136       include "DIMENSIONS.ZSCOPT"
3137       include "COMMON.VAR"
3138       include "COMMON.CHAIN"
3139       include "COMMON.ALLPROT"
3140       include "COMMON.IOUNITS"
3141       include "COMMON.GEO"
3142       include "COMMON.SBRIDGE"
3143       include "COMMON.FMATCH"
3144       double precision thetnorm
3145       integer iprot,i,j,ii
3146       do i=1,nres_zs(iprot)
3147         do j=1,3
3148           cg_MD(j,i,ii,iprot)=c(j,i)
3149         enddo
3150       enddo
3151       do i=nnt_zs(iprot),nct_zs(iprot)
3152         do j=1,3
3153           cg_MD(j,i+nres,ii,iprot)=c(j,i+nres)
3154         enddo
3155       enddo
3156 c 5/7/02 AL: store sbridge info
3157       nss_MD(ii,iprot)=nss
3158       do i=1,nss
3159         ihpb_MD(i,ii,iprot)=ihpb(i)
3160         jhpb_MD(i,ii,iprot)=jhpb(i)
3161       enddo
3162       return
3163       end
3164 c------------------------------------------------------------------------------
3165       subroutine restore_MDcoords(iprot,ii)
3166       implicit none
3167       include "DIMENSIONS"
3168       include "DIMENSIONS.ZSCOPT"
3169       include "COMMON.INTERACT"
3170       include "COMMON.VAR"
3171       include "COMMON.ALLPROT"
3172       include "COMMON.SBRIDGE"
3173       include "COMMON.CHAIN"
3174       include "COMMON.IOUNITS"
3175       include "COMMON.FMATCH"
3176       integer iprot,i,j,ii
3177       do i=1,nres_zs(iprot)
3178         do j=1,3
3179           c(j,i)=cg_MD(j,i,ii,iprot)
3180         enddo
3181       enddo
3182       do i=nnt_zs(iprot),nct_zs(iprot)
3183         do j=1,3
3184           c(j,i+nres)=cg_MD(j,i+nres,ii,iprot)
3185         enddo
3186       enddo
3187 c 5/7/02 AL: restore sbridge info
3188       nss=nss_MD(ii,iprot)
3189       do i=1,nss
3190         ihpb(i)=ihpb_MD(i,ii,iprot)+nres
3191         jhpb(i)=jhpb_MD(i,ii,iprot)+nres
3192 c        forcon(i)=fbr
3193 c        dhpb(i)=dbr
3194       enddo
3195 #ifdef DEBUG
3196         write (iout,*) "restore_ccoords",ii
3197         write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
3198         write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
3199         write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
3200         call flush(iout)
3201 #endif
3202       return
3203       end
3204 c------------------------------------------------------------------------------
3205       subroutine card_concat(card,to_upper)
3206       implicit none
3207       include 'DIMENSIONS.ZSCOPT'
3208       include "COMMON.IOUNITS"
3209       character*(*) card
3210       character*80 karta,ucase
3211       logical to_upper
3212       integer ilen
3213       external ilen
3214       read (inp,'(a)') karta
3215       if (to_upper) karta=ucase(karta)
3216       card=' '
3217       do while (karta(80:80).eq.'&')
3218         card=card(:ilen(card)+1)//karta(:79)
3219         read (inp,'(a)') karta
3220         if (to_upper) karta=ucase(karta)
3221       enddo
3222       card=card(:ilen(card)+1)//karta
3223       return
3224       end
3225 c------------------------------------------------------------------------------
3226       subroutine readi(rekord,lancuch,wartosc,default)
3227       implicit none
3228       character*(*) rekord,lancuch
3229       integer wartosc,default
3230       integer ilen,iread
3231       external ilen
3232       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
3233       if (iread.eq.0) then
3234         wartosc=default
3235         return
3236       endif
3237       iread=iread+ilen(lancuch)+1
3238       read (rekord(iread:),*) wartosc
3239       return
3240       end
3241 c----------------------------------------------------------------------------
3242       subroutine reada(rekord,lancuch,wartosc,default)
3243       implicit none
3244       character*(*) rekord,lancuch
3245       character*80 aux
3246       double precision wartosc,default
3247       integer ilen,iread
3248       external ilen
3249       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
3250       if (iread.eq.0) then
3251         wartosc=default
3252         return
3253       endif
3254       iread=iread+ilen(lancuch)+1
3255       read (rekord(iread:),*) wartosc
3256       return
3257       end
3258 c----------------------------------------------------------------------------
3259       subroutine multreadi(rekord,lancuch,tablica,dim,default)
3260       implicit none
3261       integer dim,i
3262       integer tablica(dim),default
3263       character*(*) rekord,lancuch
3264       character*80 aux
3265       integer ilen,iread
3266       external ilen
3267       do i=1,dim
3268         tablica(i)=default
3269       enddo
3270       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
3271       if (iread.eq.0) return
3272       iread=iread+ilen(lancuch)+1
3273       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
3274    10 return
3275       end
3276 c----------------------------------------------------------------------------
3277       subroutine multreada(rekord,lancuch,tablica,dim,default)
3278       implicit none
3279       integer dim,i
3280       double precision tablica(dim),default
3281       character*(*) rekord,lancuch
3282       character*80 aux
3283       integer ilen,iread
3284       external ilen
3285       do i=1,dim
3286         tablica(i)=default
3287       enddo
3288       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
3289       if (iread.eq.0) return
3290       iread=iread+ilen(lancuch)+1
3291       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
3292    10 return
3293       end
3294 c----------------------------------------------------------------------------
3295       subroutine reads(rekord,lancuch,wartosc,default)
3296       implicit none
3297       character*(*) rekord,lancuch,wartosc,default
3298       character*80 aux
3299       integer ilen,lenlan,lenrec,iread,ireade
3300       external ilen
3301       logical iblnk
3302       external iblnk
3303       lenlan=ilen(lancuch)
3304       lenrec=ilen(rekord)
3305       iread=index(rekord,lancuch(:lenlan)//"=")
3306 c      print *,"rekord",rekord," lancuch",lancuch
3307 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3308       if (iread.eq.0) then
3309         wartosc=default
3310         return
3311       endif
3312       iread=iread+lenlan+1
3313 c      print *,"iread",iread
3314 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3315       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3316         iread=iread+1
3317 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3318       enddo
3319 c      print *,"iread",iread
3320       if (iread.gt.lenrec) then
3321          wartosc=default
3322         return
3323       endif
3324       ireade=iread+1
3325 c      print *,"ireade",ireade
3326       do while (ireade.lt.lenrec .and.
3327      &   .not.iblnk(rekord(ireade:ireade)))
3328         ireade=ireade+1
3329       enddo
3330       wartosc=rekord(iread:ireade)
3331       return
3332       end
3333 c----------------------------------------------------------------------------
3334       subroutine multreads(rekord,lancuch,tablica,dim,default)
3335       implicit none
3336       integer dim,i
3337       character*(*) rekord,lancuch,tablica(dim),default
3338       character*80 aux
3339       integer ilen,lenlan,lenrec,iread,ireade
3340       external ilen
3341       logical iblnk
3342       external iblnk
3343       do i=1,dim
3344         tablica(i)=default
3345       enddo
3346       lenlan=ilen(lancuch)
3347       lenrec=ilen(rekord)
3348       iread=index(rekord,lancuch(:lenlan)//"=")
3349 c      print *,"rekord",rekord," lancuch",lancuch
3350 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3351       if (iread.eq.0) return
3352       iread=iread+lenlan+1
3353       do i=1,dim
3354 c      print *,"iread",iread
3355 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3356       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3357         iread=iread+1
3358 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3359       enddo
3360 c      print *,"iread",iread
3361       if (iread.gt.lenrec) return
3362       ireade=iread+1
3363 c      print *,"ireade",ireade
3364       do while (ireade.lt.lenrec .and.
3365      &   .not.iblnk(rekord(ireade:ireade)))
3366         ireade=ireade+1
3367       enddo
3368       tablica(i)=rekord(iread:ireade)
3369       iread=ireade+1
3370       enddo
3371       end
3372 c----------------------------------------------------------------------------
3373       subroutine split_string(rekord,tablica,dim,nsub)
3374       implicit none
3375       integer dim,nsub,i,ii,ll,kk
3376       character*(*) tablica(dim)
3377       character*(*) rekord
3378       integer ilen
3379       external ilen
3380       do i=1,dim
3381         tablica(i)=" "
3382       enddo
3383       ii=1
3384       ll = ilen(rekord)
3385       nsub=0
3386       do i=1,dim
3387 C Find the start of term name
3388         kk = 0
3389         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
3390           ii = ii+1
3391         enddo
3392 C Parse the name into TABLICA(i) until blank found
3393         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
3394           kk = kk+1
3395           tablica(i)(kk:kk)=rekord(ii:ii)
3396           ii = ii+1
3397         enddo
3398         if (kk.gt.0) nsub=nsub+1
3399         if (ii.gt.ll) return
3400       enddo
3401       return
3402       end
3403 c-------------------------------------------------------------------------------
3404       integer function iroof(n,m)
3405       ii = n/m
3406       if (ii*m .lt. n) ii=ii+1
3407       iroof = ii
3408       return
3409       end