source/unres/src_MD-M/readrtns_CSA.F

   1       subroutine readrtns
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       logical file_exist
  12 C Read force-field parameters except weights
  13       call parmread
  14 C Read job setup parameters
  15       call read_control
  16 C Read control parameters for energy minimzation if required
  17       if (minim) call read_minim
  18 C Read MCM control parameters if required
  19       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  20 C Read MD control parameters if reqjuired
  21       if (modecalc.eq.12) call read_MDpar
  22 C Read MREMD control parameters if required
  23       if (modecalc.eq.14) then
  24          call read_MDpar
  25          call read_REMDpar
  26       endif
  27 C Read MUCA control parameters if required
  28       if (lmuca) call read_muca
  29 C Read CSA control parameters if required (from fort.40 if exists
  30 C otherwise from general input file)
  31       if (modecalc.eq.8) then
  32        inquire (file="fort.40",exist=file_exist)
  33        if (.not.file_exist) call csaread
  34       endif
  35 cfmc      if (modecalc.eq.10) call mcmfread
  36 C Read molecule information, molecule geometry, energy-term weights, and
  37 C restraints if requested
  38       call molread
  39 C Print restraint information
  40 #ifdef MPI
  41       if (.not. out1file .or. me.eq.king) then
  42 #endif
  43       if (nhpb.gt.nss)
  44      &write (iout,'(a,i5,a)') "The following",nhpb-nss,
  45      & " distance constraints have been imposed"
  46       do i=nss+1,nhpb
  47         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
  48       enddo
  49 #ifdef MPI
  50       endif
  51 #endif
  52 c      print *,"Processor",myrank," leaves READRTNS"
  53       return
  54       end
  55 C-------------------------------------------------------------------------------
  56       subroutine read_control
  57 C
  58 C Read contorl data
  59 C
  60       implicit real*8 (a-h,o-z)
  61       include 'DIMENSIONS'
  62 #ifdef MP
  63       include 'mpif.h'
  64       logical OKRandom, prng_restart
  65       real*8  r1
  66 #endif
  67       include 'COMMON.IOUNITS'
  68       include 'COMMON.TIME1'
  69       include 'COMMON.THREAD'
  70       include 'COMMON.SBRIDGE'
  71       include 'COMMON.CONTROL'
  72       include 'COMMON.MCM'
  73       include 'COMMON.MAP'
  74       include 'COMMON.HEADER'
  75       include 'COMMON.CSA'
  76       include 'COMMON.CHAIN'
  77       include 'COMMON.MUCA'
  78       include 'COMMON.MD'
  79       include 'COMMON.FFIELD'
  80       include 'COMMON.INTERACT'
  81       include 'COMMON.SETUP'
  82       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  83       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  84       character*80 ucase
  85       character*320 controlcard
  86
  87       nglob_csa=0
  88       eglob_csa=1d99
  89       nmin_csa=0
  90       read (INP,'(a)') titel
  91       call card_concat(controlcard)
  92 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  93 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  94       call reada(controlcard,'SEED',seed,0.0D0)
  95       call random_init(seed)
  96 C Set up the time limit (caution! The time must be input in minutes!)
  97       read_cart=index(controlcard,'READ_CART').gt.0
  98       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  99       call readi(controlcard,'SYM',symetr,1)
 100       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
 101       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 102       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 103       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 104       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 105       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 106       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 107       call reada(controlcard,'DRMS',drms,0.1D0)
 108       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 109        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 110        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 111        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 112        write (iout,'(a,f10.1)')'DRMS    = ',drms
 113        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 114        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 115       endif
 116       call readi(controlcard,'NZ_START',nz_start,0)
 117       call readi(controlcard,'NZ_END',nz_end,0)
 118       call readi(controlcard,'IZ_SC',iz_sc,0)
 119       timlim=60.0D0*timlim
 120       safety = 60.0d0*safety
 121       timem=timlim
 122       modecalc=0
 123       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 124       minim=(index(controlcard,'MINIMIZE').gt.0)
 125       dccart=(index(controlcard,'CART').gt.0)
 126       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 127       overlapsc=.not.overlapsc
 128       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 129       searchsc=.not.searchsc
 130       sideadd=(index(controlcard,'SIDEADD').gt.0)
 131       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 132       outpdb=(index(controlcard,'PDBOUT').gt.0)
 133       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 134       pdbref=(index(controlcard,'PDBREF').gt.0)
 135       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 136       indpdb=index(controlcard,'PDBSTART')
 137       extconf=(index(controlcard,'EXTCONF').gt.0)
 138       call readi(controlcard,'IPRINT',iprint,0)
 139       call readi(controlcard,'MAXGEN',maxgen,10000)
 140       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 141       call readi(controlcard,"KDIAG",kdiag,0)
 142       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 143       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 144      & write (iout,*) "RESCALE_MODE",rescale_mode
 145       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 146       if (index(controlcard,'REGULAR').gt.0.0D0) then
 147         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 148         modecalc=1
 149         refstr=.true.
 150       endif
 151       if (index(controlcard,'CHECKGRAD').gt.0) then
 152         modecalc=5
 153         if (index(controlcard,'CART').gt.0) then
 154           icheckgrad=1
 155         elseif (index(controlcard,'CARINT').gt.0) then
 156           icheckgrad=2
 157         else
 158           icheckgrad=3
 159         endif
 160       elseif (index(controlcard,'THREAD').gt.0) then
 161         modecalc=2
 162         call readi(controlcard,'THREAD',nthread,0)
 163         if (nthread.gt.0) then
 164           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 165         else
 166           if (fg_rank.eq.0)
 167      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 168           stop 'Error termination in Read_Control.'
 169         endif
 170       else if (index(controlcard,'MCMA').gt.0) then
 171         modecalc=3
 172       else if (index(controlcard,'MCEE').gt.0) then
 173         modecalc=6
 174       else if (index(controlcard,'MULTCONF').gt.0) then
 175         modecalc=4
 176       else if (index(controlcard,'MAP').gt.0) then
 177         modecalc=7
 178         call readi(controlcard,'MAP',nmap,0)
 179       else if (index(controlcard,'CSA').gt.0) then
 180         modecalc=8
 181 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 182 crc
 183 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 184 crc
 185 crc        modecalc=9
 186 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 187 cfmc        modecalc=10
 188       else if (index(controlcard,'SOFTREG').gt.0) then
 189         modecalc=11
 190       else if (index(controlcard,'CHECK_BOND').gt.0) then
 191         modecalc=-1
 192       else if (index(controlcard,'TEST').gt.0) then
 193         modecalc=-2
 194       else if (index(controlcard,'MD').gt.0) then
 195         modecalc=12
 196       else if (index(controlcard,'RE ').gt.0) then
 197         modecalc=14
 198       endif
 199
 200       lmuca=index(controlcard,'MUCA').gt.0
 201       call readi(controlcard,'MUCADYN',mucadyn,0)
 202       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 203       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 204      & then
 205        write (iout,*) 'MUCADYN=',mucadyn
 206        write (iout,*) 'MUCASMOOTH=',muca_smooth
 207       endif
 208
 209       iscode=index(controlcard,'ONE_LETTER')
 210       indphi=index(controlcard,'PHI')
 211       indback=index(controlcard,'BACK')
 212       iranconf=index(controlcard,'RAND_CONF')
 213       i2ndstr=index(controlcard,'USE_SEC_PRED')
 214       gradout=index(controlcard,'GRADOUT').gt.0
 215       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 216       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 217       if (me.eq.king .or. .not.out1file )
 218      &  write (iout,*) "DISTCHAINMAX",distchainmax
 219
 220       if(me.eq.king.or..not.out1file)
 221      & write (iout,'(2a)') diagmeth(kdiag),
 222      &  ' routine used to diagonalize matrices.'
 223       return
 224       end
 225 c--------------------------------------------------------------------------
 226       subroutine read_REMDpar
 227 C
 228 C Read REMD settings
 229 C
 230       implicit real*8 (a-h,o-z)
 231       include 'DIMENSIONS'
 232       include 'COMMON.IOUNITS'
 233       include 'COMMON.TIME1'
 234       include 'COMMON.MD'
 235 #ifndef LANG0
 236       include 'COMMON.LANGEVIN'
 237 #else
 238       include 'COMMON.LANGEVIN.lang0'
 239 #endif
 240       include 'COMMON.INTERACT'
 241       include 'COMMON.NAMES'
 242       include 'COMMON.GEO'
 243       include 'COMMON.REMD'
 244       include 'COMMON.CONTROL'
 245       include 'COMMON.SETUP'
 246       character*80 ucase
 247       character*320 controlcard
 248       character*3200 controlcard1
 249       integer iremd_m_total
 250
 251       if(me.eq.king.or..not.out1file)
 252      & write (iout,*) "REMD setup"
 253
 254       call card_concat(controlcard)
 255       call readi(controlcard,"NREP",nrep,3)
 256       call readi(controlcard,"NSTEX",nstex,1000)
 257       call reada(controlcard,"RETMIN",retmin,10.0d0)
 258       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 259       mremdsync=(index(controlcard,'SYNC').gt.0)
 260       call readi(controlcard,"NSYN",i_sync_step,100)
 261       restart1file=(index(controlcard,'REST1FILE').gt.0)
 262       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 263       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 264       if(max_cache_traj_use.gt.max_cache_traj)
 265      &           max_cache_traj_use=max_cache_traj
 266       if(me.eq.king.or..not.out1file) then
 267 cd       if (traj1file) then
 268 crc caching is in testing - NTWX is not ignored
 269 cd        write (iout,*) "NTWX value is ignored"
 270 cd        write (iout,*) "  trajectory is stored to one file by master"
 271 cd        write (iout,*) "  before exchange at NSTEX intervals"
 272 cd       endif
 273        write (iout,*) "NREP= ",nrep
 274        write (iout,*) "NSTEX= ",nstex
 275        write (iout,*) "SYNC= ",mremdsync
 276        write (iout,*) "NSYN= ",i_sync_step
 277        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 278       endif
 279       remd_tlist=.false.
 280       if (index(controlcard,'TLIST').gt.0) then
 281          remd_tlist=.true.
 282          call card_concat(controlcard1)
 283          read(controlcard1,*) (remd_t(i),i=1,nrep)
 284          if(me.eq.king.or..not.out1file)
 285      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 286       endif
 287       remd_mlist=.false.
 288       if (index(controlcard,'MLIST').gt.0) then
 289          remd_mlist=.true.
 290          call card_concat(controlcard1)
 291          read(controlcard1,*) (remd_m(i),i=1,nrep)
 292          if(me.eq.king.or..not.out1file) then
 293           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 294           iremd_m_total=0
 295           do i=1,nrep
 296            iremd_m_total=iremd_m_total+remd_m(i)
 297           enddo
 298           write (iout,*) 'Total number of replicas ',iremd_m_total
 299          endif
 300       endif
 301       if(me.eq.king.or..not.out1file)
 302      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 303       return
 304       end
 305 c--------------------------------------------------------------------------
 306       subroutine read_MDpar
 307 C
 308 C Read MD settings
 309 C
 310       implicit real*8 (a-h,o-z)
 311       include 'DIMENSIONS'
 312       include 'COMMON.IOUNITS'
 313       include 'COMMON.TIME1'
 314       include 'COMMON.MD'
 315 #ifndef LANG0
 316       include 'COMMON.LANGEVIN'
 317 #else
 318       include 'COMMON.LANGEVIN.lang0'
 319 #endif
 320       include 'COMMON.INTERACT'
 321       include 'COMMON.NAMES'
 322       include 'COMMON.GEO'
 323       include 'COMMON.SETUP'
 324       include 'COMMON.CONTROL'
 325       include 'COMMON.SPLITELE'
 326       character*80 ucase
 327       character*320 controlcard
 328
 329       call card_concat(controlcard)
 330       call readi(controlcard,"NSTEP",n_timestep,1000000)
 331       call readi(controlcard,"NTWE",ntwe,100)
 332       call readi(controlcard,"NTWX",ntwx,1000)
 333       call reada(controlcard,"DT",d_time,1.0d-1)
 334       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 335       call reada(controlcard,"DAMAX",damax,1.0d1)
 336       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 337       call readi(controlcard,"LANG",lang,0)
 338       RESPA = index(controlcard,"RESPA") .gt. 0
 339       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 340       ntime_split0=ntime_split
 341       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 342       ntime_split0=ntime_split
 343       call reada(controlcard,"R_CUT",r_cut,2.0d0)
 344       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 345       rest = index(controlcard,"REST").gt.0
 346       tbf = index(controlcard,"TBF").gt.0
 347       usampl = index(controlcard,"USAMPL").gt.0
 348       mdpdb = index(controlcard,"MDPDB").gt.0
 349       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 350       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 351       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 352       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 353       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 354       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 355       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 356       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 357       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 358       large = index(controlcard,"LARGE").gt.0
 359       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 360       rattle = index(controlcard,"RATTLE").gt.0
 361 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 362       nset=0
 363       if(usampl) then
 364         call read_fragments
 365       endif
 366
 367       if(me.eq.king.or..not.out1file) then
 368        write (iout,*)
 369        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 370        write (iout,*)
 371        write (iout,'(a)') "The units are:"
 372        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 373        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 374      &  " acceleration: angstrom/(48.9 fs)**2"
 375        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 376        write (iout,*)
 377        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 378        write (iout,'(a60,f10.5,a)')
 379      &  "Initial time step of numerical integration:",d_time,
 380      &  " natural units"
 381        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 382        if (RESPA) then
 383         write (iout,'(2a,i4,a)')
 384      &    "A-MTS algorithm used; initial time step for fast-varying",
 385      &    " short-range forces split into",ntime_split," steps."
 386         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 387      &   r_cut," lambda",rlamb
 388        endif
 389        write (iout,'(2a,f10.5)')
 390      &  "Maximum acceleration threshold to reduce the time step",
 391      &  "/increase split number:",damax
 392        write (iout,'(2a,f10.5)')
 393      &  "Maximum predicted energy drift to reduce the timestep",
 394      &  "/increase split number:",edriftmax
 395        write (iout,'(a60,f10.5)')
 396      & "Maximum velocity threshold to reduce velocities:",dvmax
 397        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 398        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 399        if (rattle) write (iout,'(a60)')
 400      &  "Rattle algorithm used to constrain the virtual bonds"
 401       endif
 402       reset_fricmat=1000
 403       if (lang.gt.0) then
 404         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 405         call reada(controlcard,"RWAT",rwat,1.4d0)
 406         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 407         surfarea=index(controlcard,"SURFAREA").gt.0
 408         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 409         if(me.eq.king.or..not.out1file)then
 410          write (iout,'(/a,$)') "Langevin dynamics calculation"
 411          if (lang.eq.1) then
 412           write (iout,'(a/)')
 413      &      " with direct integration of Langevin equations"
 414          else if (lang.eq.2) then
 415           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 416          else if (lang.eq.3) then
 417           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 418          else if (lang.eq.4) then
 419           write (iout,'(a/)') " in overdamped mode"
 420          else
 421           write (iout,'(//a,i5)')
 422      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 423           stop
 424          endif
 425          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 426          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 427          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 428          write (iout,'(a60,f10.5)')
 429      &   "Scaling factor of the friction forces:",scal_fric
 430          if (surfarea) write (iout,'(2a,i10,a)')
 431      &     "Friction coefficients will be scaled by solvent-accessible",
 432      &     " surface area every",reset_fricmat," steps."
 433         endif
 434 c Calculate friction coefficients and bounds of stochastic forces
 435         eta=6*pi*cPoise*etawat
 436         if(me.eq.king.or..not.out1file)
 437      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 438      &   ,eta
 439         gamp=scal_fric*(pstok+rwat)*eta
 440         stdfp=dsqrt(2*Rb*t_bath/d_time)
 441         do i=1,ntyp
 442           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 443           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 444         enddo
 445         if(me.eq.king.or..not.out1file)then
 446          write (iout,'(/2a/)')
 447      &   "Radii of site types and friction coefficients and std's of",
 448      &   " stochastic forces of fully exposed sites"
 449          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 450          do i=1,ntyp
 451           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 452      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 453          enddo
 454         endif
 455       else if (tbf) then
 456         if(me.eq.king.or..not.out1file)then
 457          write (iout,'(a)') "Berendsen bath calculation"
 458          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 459          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 460          if (reset_moment)
 461      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 462      &   count_reset_moment," steps"
 463          if (reset_vel)
 464      &    write (iout,'(a,i10,a)')
 465      &    "Velocities will be reset at random every",count_reset_vel,
 466      &   " steps"
 467         endif
 468       else
 469         if(me.eq.king.or..not.out1file)
 470      &   write (iout,'(a31)') "Microcanonical mode calculation"
 471       endif
 472       if(me.eq.king.or..not.out1file)then
 473        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 474        if (usampl) then
 475           write(iout,*) "MD running with constraints."
 476           write(iout,*) "Equilibration time ", eq_time, " mtus."
 477           write(iout,*) "Constraining ", nfrag," fragments."
 478           write(iout,*) "Length of each fragment, weight and q0:"
 479           do iset=1,nset
 480            write (iout,*) "Set of restraints #",iset
 481            do i=1,nfrag
 482               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 483      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 484            enddo
 485            write(iout,*) "constraints between ", npair, "fragments."
 486            write(iout,*) "constraint pairs, weights and q0:"
 487            do i=1,npair
 488             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 489      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 490            enddo
 491            write(iout,*) "angle constraints within ", nfrag_back,
 492      &      "backbone fragments."
 493            write(iout,*) "fragment, weights:"
 494            do i=1,nfrag_back
 495             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 496      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 497      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 498            enddo
 499           enddo
 500         iset=mod(kolor,nset)+1
 501        endif
 502       endif
 503       if(me.eq.king.or..not.out1file)
 504      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 505       return
 506       end
 507 c------------------------------------------------------------------------------
 508       subroutine molread
 509 C
 510 C Read molecular data.
 511 C
 512       implicit real*8 (a-h,o-z)
 513       include 'DIMENSIONS'
 514 #ifdef MPI
 515       include 'mpif.h'
 516       integer error_msg
 517 #endif
 518       include 'COMMON.IOUNITS'
 519       include 'COMMON.GEO'
 520       include 'COMMON.VAR'
 521       include 'COMMON.INTERACT'
 522       include 'COMMON.LOCAL'
 523       include 'COMMON.NAMES'
 524       include 'COMMON.CHAIN'
 525       include 'COMMON.FFIELD'
 526       include 'COMMON.SBRIDGE'
 527       include 'COMMON.HEADER'
 528       include 'COMMON.CONTROL'
 529       include 'COMMON.DBASE'
 530       include 'COMMON.THREAD'
 531       include 'COMMON.CONTACTS'
 532       include 'COMMON.TORCNSTR'
 533       include 'COMMON.TIME1'
 534       include 'COMMON.BOUNDS'
 535       include 'COMMON.MD'
 536       include 'COMMON.SETUP'
 537       character*4 sequence(maxres)
 538       integer rescode
 539       double precision x(maxvar)
 540       character*256 pdbfile
 541       character*320 weightcard
 542       character*80 weightcard_t,ucase
 543       dimension itype_pdb(maxres)
 544       common /pizda/ itype_pdb
 545       logical seq_comp,fail
 546       double precision energia(0:n_ene)
 547       integer ilen
 548       external ilen
 549 C
 550 C Body
 551 C
 552 C Read weights of the subsequent energy terms.
 553       call card_concat(weightcard)
 554       call reada(weightcard,'WLONG',wlong,1.0D0)
 555       call reada(weightcard,'WSC',wsc,wlong)
 556       call reada(weightcard,'WSCP',wscp,wlong)
 557       call reada(weightcard,'WELEC',welec,1.0D0)
 558       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 559       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 560       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 561       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 562       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 563       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 564       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 565       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 566       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 567       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 568       call reada(weightcard,'WBOND',wbond,1.0D0)
 569       call reada(weightcard,'WTOR',wtor,1.0D0)
 570       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 571       call reada(weightcard,'WANG',wang,1.0D0)
 572       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 573       call reada(weightcard,'SCAL14',scal14,0.4D0)
 574       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 575       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 576       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 577       call reada(weightcard,'TEMP0',temp0,300.0d0)
 578       if (index(weightcard,'SOFT').gt.0) ipot=6
 579 C 12/1/95 Added weight for the multi-body term WCORR
 580       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 581       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 582       weights(1)=wsc
 583       weights(2)=wscp
 584       weights(3)=welec
 585       weights(4)=wcorr
 586       weights(5)=wcorr5
 587       weights(6)=wcorr6
 588       weights(7)=wel_loc
 589       weights(8)=wturn3
 590       weights(9)=wturn4
 591       weights(10)=wturn6
 592       weights(11)=wang
 593       weights(12)=wscloc
 594       weights(13)=wtor
 595       weights(14)=wtor_d
 596       weights(15)=wstrain
 597       weights(16)=wvdwpp
 598       weights(17)=wbond
 599       weights(18)=scal14
 600       weights(21)=wsccor
 601       if(me.eq.king.or..not.out1file)
 602      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 603      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 604      &  wturn4,wturn6
 605    10 format (/'Energy-term weights (unscaled):'//
 606      & 'WSCC=   ',f10.6,' (SC-SC)'/
 607      & 'WSCP=   ',f10.6,' (SC-p)'/
 608      & 'WELEC=  ',f10.6,' (p-p electr)'/
 609      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 610      & 'WBOND=  ',f10.6,' (stretching)'/
 611      & 'WANG=   ',f10.6,' (bending)'/
 612      & 'WSCLOC= ',f10.6,' (SC local)'/
 613      & 'WTOR=   ',f10.6,' (torsional)'/
 614      & 'WTORD=  ',f10.6,' (double torsional)'/
 615      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 616      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 617      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 618      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 619      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 620      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 621      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 622      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 623      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 624       if(me.eq.king.or..not.out1file)then
 625        if (wcorr4.gt.0.0d0) then
 626         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 627      &   'between contact pairs of peptide groups'
 628         write (iout,'(2(a,f5.3/))')
 629      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 630      &  'Range of quenching the correlation terms:',2*delt_corr
 631        else if (wcorr.gt.0.0d0) then
 632         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 633      &   'between contact pairs of peptide groups'
 634        endif
 635        write (iout,'(a,f8.3)')
 636      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 637        write (iout,'(a,f8.3)')
 638      &  'General scaling factor of SC-p interactions:',scalscp
 639       endif
 640       r0_corr=cutoff_corr-delt_corr
 641       do i=1,20
 642         aad(i,1)=scalscp*aad(i,1)
 643         aad(i,2)=scalscp*aad(i,2)
 644         bad(i,1)=scalscp*bad(i,1)
 645         bad(i,2)=scalscp*bad(i,2)
 646       enddo
 647       call rescale_weights(t_bath)
 648       if(me.eq.king.or..not.out1file)
 649      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 650      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 651      &  wturn4,wturn6
 652    22 format (/'Energy-term weights (scaled):'//
 653      & 'WSCC=   ',f10.6,' (SC-SC)'/
 654      & 'WSCP=   ',f10.6,' (SC-p)'/
 655      & 'WELEC=  ',f10.6,' (p-p electr)'/
 656      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 657      & 'WBOND=  ',f10.6,' (stretching)'/
 658      & 'WANG=   ',f10.6,' (bending)'/
 659      & 'WSCLOC= ',f10.6,' (SC local)'/
 660      & 'WTOR=   ',f10.6,' (torsional)'/
 661      & 'WTORD=  ',f10.6,' (double torsional)'/
 662      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 663      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 664      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 665      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 666      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 667      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 668      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 669      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 670      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 671       if(me.eq.king.or..not.out1file)
 672      & write (iout,*) "Reference temperature for weights calculation:",
 673      &  temp0
 674       call reada(weightcard,"D0CM",d0cm,3.78d0)
 675       call reada(weightcard,"AKCM",akcm,15.1d0)
 676       call reada(weightcard,"AKTH",akth,11.0d0)
 677       call reada(weightcard,"AKCT",akct,12.0d0)
 678       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 679       call reada(weightcard,"V2SS",v2ss,7.61d0)
 680       call reada(weightcard,"V3SS",v3ss,13.7d0)
 681       call reada(weightcard,"EBR",ebr,-5.50D0)
 682       if(me.eq.king.or..not.out1file) then
 683        write (iout,*) "Parameters of the SS-bond potential:"
 684        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 685      & " AKCT",akct
 686        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 687        write (iout,*) "EBR",ebr
 688 c       print *,'indpdb=',indpdb,' pdbref=',pdbref
 689       endif
 690       if (indpdb.gt.0 .or. pdbref) then
 691         read(inp,'(a)') pdbfile
 692         if(me.eq.king.or..not.out1file)
 693      &   write (iout,'(2a)') 'PDB data will be read from file ',
 694      &   pdbfile(:ilen(pdbfile))
 695         open(ipdbin,file=pdbfile,status='old',err=33)
 696         goto 34
 697   33    write (iout,'(a)') 'Error opening PDB file.'
 698         stop
 699   34    continue
 700 c        print *,'Begin reading pdb data'
 701         call readpdb
 702 c        print *,'Finished reading pdb data'
 703         if(me.eq.king.or..not.out1file)
 704      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 705      &   ' nstart_sup=',nstart_sup
 706         do i=1,nres
 707           itype_pdb(i)=itype(i)
 708         enddo
 709         close (ipdbin)
 710         nnt=nstart_sup
 711         nct=nstart_sup+nsup-1
 712         call contact(.false.,ncont_ref,icont_ref,co)
 713
 714         if (sideadd) then
 715          if(me.eq.king.or..not.out1file)
 716      &    write(iout,*)'Adding sidechains'
 717          maxsi=1000
 718          do i=2,nres-1
 719           iti=itype(i)
 720           if (iti.ne.10 .and. itype(i).ne.21) then
 721             nsi=0
 722             fail=.true.
 723             do while (fail.and.nsi.le.maxsi)
 724               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 725               nsi=nsi+1
 726             enddo
 727             if(fail) write(iout,*)'Adding sidechain failed for res ',
 728      &              i,' after ',nsi,' trials'
 729           endif
 730          enddo
 731         endif
 732       endif
 733       if (indpdb.eq.0) then
 734 C Read sequence if not taken from the pdb file.
 735         read (inp,*) nres
 736 c        print *,'nres=',nres
 737         if (iscode.gt.0) then
 738           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 739         else
 740           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 741         endif
 742 C Convert sequence to numeric code
 743         do i=1,nres
 744           itype(i)=rescode(i,sequence(i),iscode)
 745         enddo
 746 C Assign initial virtual bond lengths
 747         do i=2,nres
 748           vbld(i)=vbl
 749           vbld_inv(i)=vblinv
 750         enddo
 751         do i=2,nres-1
 752           vbld(i+nres)=dsc(itype(i))
 753           vbld_inv(i+nres)=dsc_inv(itype(i))
 754 c          write (iout,*) "i",i," itype",itype(i),
 755 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 756         enddo
 757       endif
 758 c      print *,nres
 759 c      print '(20i4)',(itype(i),i=1,nres)
 760       do i=1,nres
 761 #ifdef PROCOR
 762         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 763 #else
 764         if (itype(i).eq.21) then
 765 #endif
 766           itel(i)=0
 767 #ifdef PROCOR
 768         else if (itype(i+1).ne.20) then
 769 #else
 770         else if (itype(i).ne.20) then
 771 #endif
 772           itel(i)=1
 773         else
 774           itel(i)=2
 775         endif
 776       enddo
 777       if(me.eq.king.or..not.out1file)then
 778        write (iout,*) "ITEL"
 779        do i=1,nres-1
 780          write (iout,*) i,itype(i),itel(i)
 781        enddo
 782        print *,'Call Read_Bridge.'
 783       endif
 784       call read_bridge
 785 C 8/13/98 Set limits to generating the dihedral angles
 786       do i=1,nres
 787         phibound(1,i)=-pi
 788         phibound(2,i)=pi
 789       enddo
 790       read (inp,*) ndih_constr
 791       if (ndih_constr.gt.0) then
 792         read (inp,*) ftors
 793         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 794         if(me.eq.king.or..not.out1file)then
 795          write (iout,*)
 796      &   'There are',ndih_constr,' constraints on phi angles.'
 797          do i=1,ndih_constr
 798           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 799          enddo
 800         endif
 801         do i=1,ndih_constr
 802           phi0(i)=deg2rad*phi0(i)
 803           drange(i)=deg2rad*drange(i)
 804         enddo
 805         if(me.eq.king.or..not.out1file)
 806      &   write (iout,*) 'FTORS',ftors
 807         do i=1,ndih_constr
 808           ii = idih_constr(i)
 809           phibound(1,ii) = phi0(i)-drange(i)
 810           phibound(2,ii) = phi0(i)+drange(i)
 811         enddo
 812       endif
 813       nnt=1
 814 #ifdef MPI
 815       if (me.eq.king) then
 816 #endif
 817        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 818        do i=1,nres
 819          write (iout,'(a3,i5,2f10.1)')
 820      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 821        enddo
 822 #ifdef MP
 823       endif
 824 #endif
 825       nct=nres
 826 cd      print *,'NNT=',NNT,' NCT=',NCT
 827       if (itype(1).eq.21) nnt=2
 828       if (itype(nres).eq.21) nct=nct-1
 829       if (pdbref) then
 830         if(me.eq.king.or..not.out1file)
 831      &   write (iout,'(a,i3)') 'nsup=',nsup
 832         nstart_seq=nnt
 833         if (nsup.le.(nct-nnt+1)) then
 834           do i=0,nct-nnt+1-nsup
 835             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 836               nstart_seq=nnt+i
 837               goto 111
 838             endif
 839           enddo
 840           write (iout,'(a)')
 841      &            'Error - sequences to be superposed do not match.'
 842           stop
 843         else
 844           do i=0,nsup-(nct-nnt+1)
 845             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 846      &      then
 847               nstart_sup=nstart_sup+i
 848               nsup=nct-nnt+1
 849               goto 111
 850             endif
 851           enddo
 852           write (iout,'(a)')
 853      &            'Error - sequences to be superposed do not match.'
 854         endif
 855   111   continue
 856         if (nsup.eq.0) nsup=nct-nnt
 857         if (nstart_sup.eq.0) nstart_sup=nnt
 858         if (nstart_seq.eq.0) nstart_seq=nnt
 859         if(me.eq.king.or..not.out1file)
 860      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 861      &                 ' nstart_seq=',nstart_seq
 862       endif
 863 c--- Zscore rms -------
 864       if (nz_start.eq.0) nz_start=nnt
 865       if (nz_end.eq.0 .and. nsup.gt.0) then
 866         nz_end=nnt+nsup-1
 867       else if (nz_end.eq.0) then
 868         nz_end=nct
 869       endif
 870       if(me.eq.king.or..not.out1file)then
 871        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 872        write (iout,*) 'IZ_SC=',iz_sc
 873       endif
 874 c----------------------
 875       call init_int_table
 876       if (refstr) then
 877         if (.not.pdbref) then
 878           call read_angles(inp,*38)
 879           goto 39
 880    38     write (iout,'(a)') 'Error reading reference structure.'
 881 #ifdef MPI
 882           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 883           stop 'Error reading reference structure'
 884 #endif
 885    39     call chainbuild
 886           call setup_var
 887 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 888           nstart_sup=nnt
 889           nstart_seq=nnt
 890           nsup=nct-nnt+1
 891           kkk=1
 892           do i=1,2*nres
 893             do j=1,3
 894               cref(j,i,kkk)=c(j,i)
 895             enddo
 896           enddo
 897           call contact(.true.,ncont_ref,icont_ref,co)
 898         endif
 899 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 900         call flush(iout)
 901         if (constr_dist.gt.0) call read_dist_constr
 902         write (iout,*) "After read_dist_constr nhpb",nhpb
 903         call hpb_partition
 904         if(me.eq.king.or..not.out1file)
 905      &   write (iout,*) 'Contact order:',co
 906         if (pdbref) then
 907         if(me.eq.king.or..not.out1file)
 908      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 909         do i=1,ncont_ref
 910           do j=1,2
 911             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 912           enddo
 913           if(me.eq.king.or..not.out1file)
 914      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 915      &     icont_ref(1,i),' ',
 916      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 917         enddo
 918         endif
 919       endif
 920       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 921      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 922      &    modecalc.ne.10) then
 923 C If input structure hasn't been supplied from the PDB file read or generate
 924 C initial geometry.
 925         if (iranconf.eq.0 .and. .not. extconf) then
 926           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 927      &     write (iout,'(a)') 'Initial geometry will be read in.'
 928           if (read_cart) then
 929             read(inp,'(8f10.5)',end=36,err=36)
 930      &       ((c(l,k),l=1,3),k=1,nres),
 931      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 932             write (iout,*) "Exit READ_CART"
 933             write (iout,'(8f10.5)')
 934      &       ((c(l,k),l=1,3),k=1,nres),
 935      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 936             call int_from_cart1(.true.)
 937             write (iout,*) "Finish INT_TO_CART"
 938             do i=1,nres-1
 939               do j=1,3
 940                 dc(j,i)=c(j,i+1)-c(j,i)
 941                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 942               enddo
 943             enddo
 944             do i=nnt,nct
 945               if (itype(i).ne.10 .and. itype(i).ne.21) then
 946                 do j=1,3
 947                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 948                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 949                 enddo
 950               endif
 951             enddo
 952             return
 953           else
 954             call read_angles(inp,*36)
 955           endif
 956           goto 37
 957    36     write (iout,'(a)') 'Error reading angle file.'
 958 #ifdef MPI
 959           call mpi_finalize( MPI_COMM_WORLD,IERR )
 960 #endif
 961           stop 'Error reading angle file.'
 962    37     continue
 963         else if (extconf) then
 964          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 965      &    write (iout,'(a)') 'Extended chain initial geometry.'
 966          do i=3,nres
 967           theta(i)=90d0*deg2rad
 968          enddo
 969          do i=4,nres
 970           phi(i)=180d0*deg2rad
 971          enddo
 972          do i=2,nres-1
 973           alph(i)=110d0*deg2rad
 974          enddo
 975          do i=2,nres-1
 976           omeg(i)=-120d0*deg2rad
 977          enddo
 978         else
 979           if(me.eq.king.or..not.out1file)
 980      &     write (iout,'(a)') 'Random-generated initial geometry.'
 981
 982
 983 #ifdef MPI
 984           if (me.eq.king  .or. fg_rank.eq.0 .and. (
 985      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
 986 #endif
 987             do itrial=1,100
 988               itmp=1
 989               call gen_rand_conf(itmp,*30)
 990               goto 40
 991    30         write (iout,*) 'Failed to generate random conformation',
 992      &          ', itrial=',itrial
 993               write (*,*) 'Processor:',me,
 994      &          ' Failed to generate random conformation',
 995      &          ' itrial=',itrial
 996               call intout
 997
 998 #ifdef AIX
 999               call flush_(iout)
1000 #else
1001               call flush(iout)
1002 #endif
1003             enddo
1004             write (iout,'(a,i3,a)') 'Processor:',me,
1005      &        ' error in generating random conformation.'
1006             write (*,'(a,i3,a)') 'Processor:',me,
1007      &        ' error in generating random conformation.'
1008             call flush(iout)
1009 #ifdef MPI
1010             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1011    40       continue
1012           endif
1013 #else
1014           do itrial=1,100
1015             itmp=1
1016             call gen_rand_conf(itmp,*30)
1017             goto 40
1018    30       write (iout,*) 'Failed to generate random conformation',
1019      &        ', itrial=',itrial
1020             write (*,*) 'Failed to generate random conformation',
1021      &        ', itrial=',itrial
1022           enddo
1023           write (iout,'(a,i3,a)') 'Processor:',me,
1024      &      ' error in generating random conformation.'
1025           write (*,'(a,i3,a)') 'Processor:',me,
1026      &      ' error in generating random conformation.'
1027           stop
1028    40     continue
1029 #endif
1030         endif
1031       elseif (modecalc.eq.4) then
1032         read (inp,'(a)') intinname
1033         open (intin,file=intinname,status='old',err=333)
1034         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1035      &  write (iout,'(a)') 'intinname',intinname
1036         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1037         goto 334
1038   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1039 #ifdef MPI
1040         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1041 #endif
1042         stop 'Error opening angle file.'
1043   334   continue
1044
1045       endif
1046 C Generate distance constraints, if the PDB structure is to be regularized.
1047       if (nthread.gt.0) then
1048         call read_threadbase
1049       endif
1050       call setup_var
1051       if (me.eq.king .or. .not. out1file)
1052      & call intout
1053       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1054         write (iout,'(/a,i3,a)')
1055      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1056         write (iout,'(20i4)') (iss(i),i=1,ns)
1057         write (iout,'(/a/)') 'Pre-formed links are:'
1058         do i=1,nss
1059           i1=ihpb(i)-nres
1060           i2=jhpb(i)-nres
1061           it1=itype(i1)
1062           it2=itype(i2)
1063           if (me.eq.king.or..not.out1file)
1064      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1065      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1066      &    ebr,forcon(i)
1067         enddo
1068         write (iout,'(a)')
1069       endif
1070       if (i2ndstr.gt.0) call secstrp2dihc
1071 c      call geom_to_var(nvar,x)
1072 c      call etotal(energia(0))
1073 c      call enerprint(energia(0))
1074 c      call briefout(0,etot)
1075 c      stop
1076 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1077 cd    write (iout,'(a)') 'Variable list:'
1078 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1079 #ifdef MPI
1080       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1081      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1082      &  'Processor',myrank,': end reading molecular data.'
1083 #endif
1084       return
1085       end
1086 c--------------------------------------------------------------------------
1087       logical function seq_comp(itypea,itypeb,length)
1088       implicit none
1089       integer length,itypea(length),itypeb(length)
1090       integer i
1091       do i=1,length
1092         if (itypea(i).ne.itypeb(i)) then
1093           seq_comp=.false.
1094           return
1095         endif
1096       enddo
1097       seq_comp=.true.
1098       return
1099       end
1100 c-----------------------------------------------------------------------------
1101       subroutine read_bridge
1102 C Read information about disulfide bridges.
1103       implicit real*8 (a-h,o-z)
1104       include 'DIMENSIONS'
1105 #ifdef MPI
1106       include 'mpif.h'
1107 #endif
1108       include 'COMMON.IOUNITS'
1109       include 'COMMON.GEO'
1110       include 'COMMON.VAR'
1111       include 'COMMON.INTERACT'
1112       include 'COMMON.LOCAL'
1113       include 'COMMON.NAMES'
1114       include 'COMMON.CHAIN'
1115       include 'COMMON.FFIELD'
1116       include 'COMMON.SBRIDGE'
1117       include 'COMMON.HEADER'
1118       include 'COMMON.CONTROL'
1119       include 'COMMON.DBASE'
1120       include 'COMMON.THREAD'
1121       include 'COMMON.TIME1'
1122       include 'COMMON.SETUP'
1123 C Read bridging residues.
1124       read (inp,*) ns,(iss(i),i=1,ns)
1125 c      print *,'ns=',ns
1126       if(me.eq.king.or..not.out1file)
1127      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1128 C Check whether the specified bridging residues are cystines.
1129       do i=1,ns
1130         if (itype(iss(i)).ne.1) then
1131           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1132      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
1133      &   ' can form a disulfide bridge?!!!'
1134           write (*,'(2a,i3,a)')
1135      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
1136      &   ' can form a disulfide bridge?!!!'
1137 #ifdef MPI
1138          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1139          stop
1140 #endif
1141         endif
1142       enddo
1143 C Read preformed bridges.
1144       if (ns.gt.0) then
1145       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1146       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1147       if (nss.gt.0) then
1148         nhpb=nss
1149 C Check if the residues involved in bridges are in the specified list of
1150 C bridging residues.
1151         do i=1,nss
1152           do j=1,i-1
1153             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1154      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1155               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1156      &      ' contains residues present in other pairs.'
1157               write (*,'(a,i3,a)') 'Disulfide pair',i,
1158      &      ' contains residues present in other pairs.'
1159 #ifdef MPI
1160               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1161               stop
1162 #endif
1163             endif
1164           enddo
1165           do j=1,ns
1166             if (ihpb(i).eq.iss(j)) goto 10
1167           enddo
1168           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1169    10     continue
1170           do j=1,ns
1171             if (jhpb(i).eq.iss(j)) goto 20
1172           enddo
1173           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1174    20     continue
1175           dhpb(i)=dbr
1176           forcon(i)=fbr
1177         enddo
1178         do i=1,nss
1179           ihpb(i)=ihpb(i)+nres
1180           jhpb(i)=jhpb(i)+nres
1181         enddo
1182       endif
1183       endif
1184       return
1185       end
1186 c----------------------------------------------------------------------------
1187       subroutine read_x(kanal,*)
1188       implicit real*8 (a-h,o-z)
1189       include 'DIMENSIONS'
1190       include 'COMMON.GEO'
1191       include 'COMMON.VAR'
1192       include 'COMMON.CHAIN'
1193       include 'COMMON.IOUNITS'
1194       include 'COMMON.CONTROL'
1195       include 'COMMON.LOCAL'
1196       include 'COMMON.INTERACT'
1197 c Read coordinates from input
1198 c
1199       read(kanal,'(8f10.5)',end=10,err=10)
1200      &  ((c(l,k),l=1,3),k=1,nres),
1201      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1202       do j=1,3
1203         c(j,nres+1)=c(j,1)
1204         c(j,2*nres)=c(j,nres)
1205       enddo
1206       call int_from_cart1(.false.)
1207       do i=1,nres-1
1208         do j=1,3
1209           dc(j,i)=c(j,i+1)-c(j,i)
1210           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1211         enddo
1212       enddo
1213       do i=nnt,nct
1214         if (itype(i).ne.10 .and. itype(i).ne.21) then
1215           do j=1,3
1216             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1217             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1218           enddo
1219         endif
1220       enddo
1221
1222       return
1223    10 return1
1224       end
1225 c----------------------------------------------------------------------------
1226       subroutine read_threadbase
1227       implicit real*8 (a-h,o-z)
1228       include 'DIMENSIONS'
1229       include 'COMMON.IOUNITS'
1230       include 'COMMON.GEO'
1231       include 'COMMON.VAR'
1232       include 'COMMON.INTERACT'
1233       include 'COMMON.LOCAL'
1234       include 'COMMON.NAMES'
1235       include 'COMMON.CHAIN'
1236       include 'COMMON.FFIELD'
1237       include 'COMMON.SBRIDGE'
1238       include 'COMMON.HEADER'
1239       include 'COMMON.CONTROL'
1240       include 'COMMON.DBASE'
1241       include 'COMMON.THREAD'
1242       include 'COMMON.TIME1'
1243 C Read pattern database for threading.
1244       read (icbase,*) nseq
1245       do i=1,nseq
1246         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1247      &   nres_base(2,i),nres_base(3,i)
1248         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1249      &   nres_base(1,i))
1250 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1251 c    &   nres_base(2,i),nres_base(3,i)
1252 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1253 c    &   nres_base(1,i))
1254       enddo
1255       close (icbase)
1256       if (weidis.eq.0.0D0) weidis=0.1D0
1257       do i=nnt,nct
1258         do j=i+2,nct
1259           nhpb=nhpb+1
1260           ihpb(nhpb)=i
1261           jhpb(nhpb)=j
1262           forcon(nhpb)=weidis
1263         enddo
1264       enddo
1265       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1266       write (iout,'(a,i5)') 'nexcl: ',nexcl
1267       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1268       return
1269       end
1270 c------------------------------------------------------------------------------
1271       subroutine setup_var
1272       implicit real*8 (a-h,o-z)
1273       include 'DIMENSIONS'
1274       include 'COMMON.IOUNITS'
1275       include 'COMMON.GEO'
1276       include 'COMMON.VAR'
1277       include 'COMMON.INTERACT'
1278       include 'COMMON.LOCAL'
1279       include 'COMMON.NAMES'
1280       include 'COMMON.CHAIN'
1281       include 'COMMON.FFIELD'
1282       include 'COMMON.SBRIDGE'
1283       include 'COMMON.HEADER'
1284       include 'COMMON.CONTROL'
1285       include 'COMMON.DBASE'
1286       include 'COMMON.THREAD'
1287       include 'COMMON.TIME1'
1288 C Set up variable list.
1289       ntheta=nres-2
1290       nphi=nres-3
1291       nvar=ntheta+nphi
1292       nside=0
1293       do i=2,nres-1
1294         if (itype(i).ne.10 .and. itype(i).ne.21) then
1295           nside=nside+1
1296           ialph(i,1)=nvar+nside
1297           ialph(nside,2)=i
1298         endif
1299       enddo
1300       if (indphi.gt.0) then
1301         nvar=nphi
1302       else if (indback.gt.0) then
1303         nvar=nphi+ntheta
1304       else
1305         nvar=nvar+2*nside
1306       endif
1307 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1308       return
1309       end
1310 c----------------------------------------------------------------------------
1311       subroutine gen_dist_constr
1312 C Generate CA distance constraints.
1313       implicit real*8 (a-h,o-z)
1314       include 'DIMENSIONS'
1315       include 'COMMON.IOUNITS'
1316       include 'COMMON.GEO'
1317       include 'COMMON.VAR'
1318       include 'COMMON.INTERACT'
1319       include 'COMMON.LOCAL'
1320       include 'COMMON.NAMES'
1321       include 'COMMON.CHAIN'
1322       include 'COMMON.FFIELD'
1323       include 'COMMON.SBRIDGE'
1324       include 'COMMON.HEADER'
1325       include 'COMMON.CONTROL'
1326       include 'COMMON.DBASE'
1327       include 'COMMON.THREAD'
1328       include 'COMMON.TIME1'
1329       dimension itype_pdb(maxres)
1330       common /pizda/ itype_pdb
1331       character*2 iden
1332 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1333 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1334 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1335 cd     & ' nsup',nsup
1336       do i=nstart_sup,nstart_sup+nsup-1
1337 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1338 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1339         do j=i+2,nstart_sup+nsup-1
1340           nhpb=nhpb+1
1341           ihpb(nhpb)=i+nstart_seq-nstart_sup
1342           jhpb(nhpb)=j+nstart_seq-nstart_sup
1343           forcon(nhpb)=weidis
1344           dhpb(nhpb)=dist(i,j)
1345         enddo
1346       enddo
1347 cd      write (iout,'(a)') 'Distance constraints:'
1348 cd      do i=nss+1,nhpb
1349 cd        ii=ihpb(i)
1350 cd        jj=jhpb(i)
1351 cd        iden='CA'
1352 cd        if (ii.gt.nres) then
1353 cd          iden='SC'
1354 cd          ii=ii-nres
1355 cd          jj=jj-nres
1356 cd        endif
1357 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1358 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1359 cd     &  dhpb(i),forcon(i)
1360 cd      enddo
1361       return
1362       end
1363 c----------------------------------------------------------------------------
1364       subroutine map_read
1365       implicit real*8 (a-h,o-z)
1366       include 'DIMENSIONS'
1367       include 'COMMON.MAP'
1368       include 'COMMON.IOUNITS'
1369       character*3 angid(4) /'THE','PHI','ALP','OME'/
1370       character*80 mapcard,ucase
1371       do imap=1,nmap
1372         read (inp,'(a)') mapcard
1373         mapcard=ucase(mapcard)
1374         if (index(mapcard,'PHI').gt.0) then
1375           kang(imap)=1
1376         else if (index(mapcard,'THE').gt.0) then
1377           kang(imap)=2
1378         else if (index(mapcard,'ALP').gt.0) then
1379           kang(imap)=3
1380         else if (index(mapcard,'OME').gt.0) then
1381           kang(imap)=4
1382         else
1383           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1384           stop 'Error - illegal variable spec in MAP card.'
1385         endif
1386         call readi (mapcard,'RES1',res1(imap),0)
1387         call readi (mapcard,'RES2',res2(imap),0)
1388         if (res1(imap).eq.0) then
1389           res1(imap)=res2(imap)
1390         else if (res2(imap).eq.0) then
1391           res2(imap)=res1(imap)
1392         endif
1393         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1394           write (iout,'(a)')
1395      &    'Error - illegal definition of variable group in MAP.'
1396           stop 'Error - illegal definition of variable group in MAP.'
1397         endif
1398         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1399         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1400         call readi(mapcard,'NSTEP',nstep(imap),0)
1401         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1402           write (iout,'(a)')
1403      &     'Illegal boundary and/or step size specification in MAP.'
1404           stop 'Illegal boundary and/or step size specification in MAP.'
1405         endif
1406       enddo ! imap
1407       return
1408       end
1409 c----------------------------------------------------------------------------
1410       subroutine csaread
1411       implicit real*8 (a-h,o-z)
1412       include 'DIMENSIONS'
1413       include 'COMMON.IOUNITS'
1414       include 'COMMON.GEO'
1415       include 'COMMON.CSA'
1416       include 'COMMON.BANK'
1417       include 'COMMON.CONTROL'
1418       character*80 ucase
1419       character*620 mcmcard
1420       call card_concat(mcmcard)
1421
1422       call readi(mcmcard,'NCONF',nconf,50)
1423       call readi(mcmcard,'NADD',nadd,0)
1424       call readi(mcmcard,'JSTART',jstart,1)
1425       call readi(mcmcard,'JEND',jend,1)
1426       call readi(mcmcard,'NSTMAX',nstmax,500000)
1427       call readi(mcmcard,'N0',n0,1)
1428       call readi(mcmcard,'N1',n1,6)
1429       call readi(mcmcard,'N2',n2,4)
1430       call readi(mcmcard,'N3',n3,0)
1431       call readi(mcmcard,'N4',n4,0)
1432       call readi(mcmcard,'N5',n5,0)
1433       call readi(mcmcard,'N6',n6,10)
1434       call readi(mcmcard,'N7',n7,0)
1435       call readi(mcmcard,'N8',n8,0)
1436       call readi(mcmcard,'N9',n9,0)
1437       call readi(mcmcard,'N14',n14,0)
1438       call readi(mcmcard,'N15',n15,0)
1439       call readi(mcmcard,'N16',n16,0)
1440       call readi(mcmcard,'N17',n17,0)
1441       call readi(mcmcard,'N18',n18,0)
1442
1443       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1444
1445       call readi(mcmcard,'NDIFF',ndiff,2)
1446       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1447       call readi(mcmcard,'IS1',is1,1)
1448       call readi(mcmcard,'IS2',is2,8)
1449       call readi(mcmcard,'NRAN0',nran0,4)
1450       call readi(mcmcard,'NRAN1',nran1,2)
1451       call readi(mcmcard,'IRR',irr,1)
1452       call readi(mcmcard,'NSEED',nseed,20)
1453       call readi(mcmcard,'NTOTAL',ntotal,10000)
1454       call reada(mcmcard,'CUT1',cut1,2.0d0)
1455       call reada(mcmcard,'CUT2',cut2,5.0d0)
1456       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1457       call readi(mcmcard,'ICMAX',icmax,3)
1458       call readi(mcmcard,'IRESTART',irestart,0)
1459 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1460       ntbankm=0
1461 c!bankt
1462       call reada(mcmcard,'DELE',dele,20.0d0)
1463       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1464       call readi(mcmcard,'IREF',iref,0)
1465       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1466       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1467       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1468       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1469       write (iout,*) "NCONF_IN",nconf_in
1470       return
1471       end
1472 c----------------------------------------------------------------------------
1473 cfmc      subroutine mcmfread
1474 cfmc      implicit real*8 (a-h,o-z)
1475 cfmc      include 'DIMENSIONS'
1476 cfmc      include 'COMMON.MCMF'
1477 cfmc      include 'COMMON.IOUNITS'
1478 cfmc      include 'COMMON.GEO'
1479 cfmc      character*80 ucase
1480 cfmc      character*620 mcmcard
1481 cfmc      call card_concat(mcmcard)
1482 cfmc
1483 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1484 cfmc      write(iout,*)'MAXRANT=',maxrant
1485 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1486 cfmc      write(iout,*)'MAXFAM=',maxfam
1487 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1488 cfmc      write(iout,*)'NNET1=',nnet1
1489 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1490 cfmc      write(iout,*)'NNET2=',nnet2
1491 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1492 cfmc      write(iout,*)'NNET3=',nnet3
1493 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1494 cfmc      write(iout,*)'ILASTT=',ilastt
1495 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1496 cfmc      write(iout,*)'MAXSTR=',maxstr
1497 cfmc      maxstr_f=maxstr/maxfam
1498 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1499 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1500 cfmc      write(iout,*)'NMCMF=',nmcmf
1501 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1502 cfmc      write(iout,*)'IFOCUS=',ifocus
1503 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1504 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1505 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1506 cfmc      write(iout,*)'INTPRT=',intprt
1507 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1508 cfmc      write(iout,*)'IPRT=',iprt
1509 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1510 cfmc      write(iout,*)'IMAXTR=',imaxtr
1511 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1512 cfmc      write(iout,*)'MAXEVEN=',maxeven
1513 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1514 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1515 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1516 cfmc      write(iout,*)'INIMIN=',inimin
1517 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1518 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1519 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1520 cfmc      write(iout,*)'NTHREAD=',nthread
1521 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1522 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1523 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1524 cfmc      write(iout,*)'MAXPERT=',maxpert
1525 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1526 cfmc      write(iout,*)'IRMSD=',irmsd
1527 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1528 cfmc      write(iout,*)'DENEMIN=',denemin
1529 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1530 cfmc      write(iout,*)'RCUT1S=',rcut1s
1531 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1532 cfmc      write(iout,*)'RCUT1E=',rcut1e
1533 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1534 cfmc      write(iout,*)'RCUT2S=',rcut2s
1535 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1536 cfmc      write(iout,*)'RCUT2E=',rcut2e
1537 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1538 cfmc      write(iout,*)'DPERT1=',d_pert1
1539 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1540 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1541 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1542 cfmc      write(iout,*)'DPERT2=',d_pert2
1543 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1544 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1545 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1546 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1547 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1548 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1549 cfmc      d_pert1=deg2rad*d_pert1
1550 cfmc      d_pert1a=deg2rad*d_pert1a
1551 cfmc      d_pert2=deg2rad*d_pert2
1552 cfmc      d_pert2a=deg2rad*d_pert2a
1553 cfmc      d_pert2b=deg2rad*d_pert2b
1554 cfmc      d_pert2c=deg2rad*d_pert2c
1555 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1556 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1557 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1558 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1559 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1560 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1561 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1562 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1563 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1564 cfmc      write(iout,*)'RCUTINI=',rcutini
1565 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1566 cfmc      write(iout,*)'GRAT=',grat
1567 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1568 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1569 cfmc
1570 cfmc      return
1571 cfmc      end
1572 c----------------------------------------------------------------------------
1573       subroutine mcmread
1574       implicit real*8 (a-h,o-z)
1575       include 'DIMENSIONS'
1576       include 'COMMON.MCM'
1577       include 'COMMON.MCE'
1578       include 'COMMON.IOUNITS'
1579       character*80 ucase
1580       character*320 mcmcard
1581       call card_concat(mcmcard)
1582       call readi(mcmcard,'MAXACC',maxacc,100)
1583       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1584       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1585       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1586       call readi(mcmcard,'MAXREPM',maxrepm,200)
1587       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1588       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1589       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1590       call reada(mcmcard,'E_UP',e_up,5.0D0)
1591       call reada(mcmcard,'DELTE',delte,0.1D0)
1592       call readi(mcmcard,'NSWEEP',nsweep,5)
1593       call readi(mcmcard,'NSTEPH',nsteph,0)
1594       call readi(mcmcard,'NSTEPC',nstepc,0)
1595       call reada(mcmcard,'TMIN',tmin,298.0D0)
1596       call reada(mcmcard,'TMAX',tmax,298.0D0)
1597       call readi(mcmcard,'NWINDOW',nwindow,0)
1598       call readi(mcmcard,'PRINT_MC',print_mc,0)
1599       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1600       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1601       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1602       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1603       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1604       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1605       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1606       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1607       if (nwindow.gt.0) then
1608         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1609         do i=1,nwindow
1610           winlen(i)=winend(i)-winstart(i)+1
1611         enddo
1612       endif
1613       if (tmax.lt.tmin) tmax=tmin
1614       if (tmax.eq.tmin) then
1615         nstepc=0
1616         nsteph=0
1617       endif
1618       if (nstepc.gt.0 .and. nsteph.gt.0) then
1619         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1620         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1621       endif
1622 C Probabilities of different move types
1623       sumpro_type(0)=0.0D0
1624       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1625       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1626       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1627       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1628       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1629       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1630       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1631       do i=1,MaxMoveType
1632         print *,'i',i,' sumprotype',sumpro_type(i)
1633         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1634         print *,'i',i,' sumprotype',sumpro_type(i)
1635       enddo
1636       return
1637       end
1638 c----------------------------------------------------------------------------
1639       subroutine read_minim
1640       implicit real*8 (a-h,o-z)
1641       include 'DIMENSIONS'
1642       include 'COMMON.MINIM'
1643       include 'COMMON.IOUNITS'
1644       character*80 ucase
1645       character*320 minimcard
1646       call card_concat(minimcard)
1647       call readi(minimcard,'MAXMIN',maxmin,2000)
1648       call readi(minimcard,'MAXFUN',maxfun,5000)
1649       call readi(minimcard,'MINMIN',minmin,maxmin)
1650       call readi(minimcard,'MINFUN',minfun,maxmin)
1651       call reada(minimcard,'TOLF',tolf,1.0D-2)
1652       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1653       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1654       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1655       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1656       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1657      &         'Options in energy minimization:'
1658       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1659      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1660      & 'MinMin:',MinMin,' MinFun:',MinFun,
1661      & ' TolF:',TolF,' RTolF:',RTolF
1662       return
1663       end
1664 c----------------------------------------------------------------------------
1665       subroutine read_angles(kanal,*)
1666       implicit real*8 (a-h,o-z)
1667       include 'DIMENSIONS'
1668       include 'COMMON.GEO'
1669       include 'COMMON.VAR'
1670       include 'COMMON.CHAIN'
1671       include 'COMMON.IOUNITS'
1672       include 'COMMON.CONTROL'
1673 c Read angles from input
1674 c
1675        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1676        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1677        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1678        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1679
1680        do i=1,nres
1681 c 9/7/01 avoid 180 deg valence angle
1682         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1683 c
1684         theta(i)=deg2rad*theta(i)
1685         phi(i)=deg2rad*phi(i)
1686         alph(i)=deg2rad*alph(i)
1687         omeg(i)=deg2rad*omeg(i)
1688        enddo
1689       return
1690    10 return1
1691       end
1692 c----------------------------------------------------------------------------
1693       subroutine reada(rekord,lancuch,wartosc,default)
1694       implicit none
1695       character*(*) rekord,lancuch
1696       double precision wartosc,default
1697       integer ilen,iread
1698       external ilen
1699       iread=index(rekord,lancuch)
1700       if (iread.eq.0) then
1701         wartosc=default
1702         return
1703       endif
1704       iread=iread+ilen(lancuch)+1
1705       read (rekord(iread:),*,err=10,end=10) wartosc
1706       return
1707   10  wartosc=default
1708       return
1709       end
1710 c----------------------------------------------------------------------------
1711       subroutine readi(rekord,lancuch,wartosc,default)
1712       implicit none
1713       character*(*) rekord,lancuch
1714       integer wartosc,default
1715       integer ilen,iread
1716       external ilen
1717       iread=index(rekord,lancuch)
1718       if (iread.eq.0) then
1719         wartosc=default
1720         return
1721       endif
1722       iread=iread+ilen(lancuch)+1
1723       read (rekord(iread:),*,err=10,end=10) wartosc
1724       return
1725   10  wartosc=default
1726       return
1727       end
1728 c----------------------------------------------------------------------------
1729       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1730       implicit none
1731       integer dim,i
1732       integer tablica(dim),default
1733       character*(*) rekord,lancuch
1734       character*80 aux
1735       integer ilen,iread
1736       external ilen
1737       do i=1,dim
1738         tablica(i)=default
1739       enddo
1740       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1741       if (iread.eq.0) return
1742       iread=iread+ilen(lancuch)+1
1743       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1744    10 return
1745       end
1746 c----------------------------------------------------------------------------
1747       subroutine multreada(rekord,lancuch,tablica,dim,default)
1748       implicit none
1749       integer dim,i
1750       double precision tablica(dim),default
1751       character*(*) rekord,lancuch
1752       character*80 aux
1753       integer ilen,iread
1754       external ilen
1755       do i=1,dim
1756         tablica(i)=default
1757       enddo
1758       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1759       if (iread.eq.0) return
1760       iread=iread+ilen(lancuch)+1
1761       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1762    10 return
1763       end
1764 c----------------------------------------------------------------------------
1765       subroutine openunits
1766       implicit real*8 (a-h,o-z)
1767       include 'DIMENSIONS'
1768 #ifdef MPI
1769       include 'mpif.h'
1770       character*16 form,nodename
1771       integer nodelen
1772 #endif
1773       include 'COMMON.SETUP'
1774       include 'COMMON.IOUNITS'
1775       include 'COMMON.MD'
1776       include 'COMMON.CONTROL'
1777       integer lenpre,lenpot,ilen,lentmp
1778       external ilen
1779       character*3 out1file_text,ucase
1780       character*3 ll
1781       external ucase
1782 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1783       call getenv_loc("PREFIX",prefix)
1784       pref_orig = prefix
1785       call getenv_loc("POT",pot)
1786       call getenv_loc("DIRTMP",tmpdir)
1787       call getenv_loc("CURDIR",curdir)
1788       call getenv_loc("OUT1FILE",out1file_text)
1789 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1790       out1file_text=ucase(out1file_text)
1791       if (out1file_text(1:1).eq."Y") then
1792         out1file=.true.
1793       else
1794         out1file=fg_rank.gt.0
1795       endif
1796       lenpre=ilen(prefix)
1797       lenpot=ilen(pot)
1798       lentmp=ilen(tmpdir)
1799       if (lentmp.gt.0) then
1800           write (*,'(80(1h!))')
1801           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1802           write (*,'(80(1h!))')
1803           write (*,*)"All output files will be on node /tmp directory."
1804 #ifdef MPI
1805         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1806         if (me.eq.king) then
1807           write (*,*) "The master node is ",nodename
1808         else if (fg_rank.eq.0) then
1809           write (*,*) "I am the CG slave node ",nodename
1810         else
1811           write (*,*) "I am the FG slave node ",nodename
1812         endif
1813 #endif
1814         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1815         lenpre = lentmp+lenpre+1
1816       endif
1817       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1818 C Get the names and open the input files
1819 #if defined(WINIFL) || defined(WINPGI)
1820       open(1,file=pref_orig(:ilen(pref_orig))//
1821      &  '.inp',status='old',readonly,shared)
1822        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1823 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1824 C Get parameter filenames and open the parameter files.
1825       call getenv_loc('BONDPAR',bondname)
1826       open (ibond,file=bondname,status='old',readonly,shared)
1827       call getenv_loc('THETPAR',thetname)
1828       open (ithep,file=thetname,status='old',readonly,shared)
1829       call getenv_loc('ROTPAR',rotname)
1830       open (irotam,file=rotname,status='old',readonly,shared)
1831       call getenv_loc('TORPAR',torname)
1832       open (itorp,file=torname,status='old',readonly,shared)
1833       call getenv_loc('TORDPAR',tordname)
1834       open (itordp,file=tordname,status='old',readonly,shared)
1835       call getenv_loc('FOURIER',fouriername)
1836       open (ifourier,file=fouriername,status='old',readonly,shared)
1837       call getenv_loc('ELEPAR',elename)
1838       open (ielep,file=elename,status='old',readonly,shared)
1839       call getenv_loc('SIDEPAR',sidename)
1840       open (isidep,file=sidename,status='old',readonly,shared)
1841 #elif (defined CRAY) || (defined AIX)
1842       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1843      &  action='read')
1844 c      print *,"Processor",myrank," opened file 1"
1845       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1846 c      print *,"Processor",myrank," opened file 9"
1847 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1848 C Get parameter filenames and open the parameter files.
1849       call getenv_loc('BONDPAR',bondname)
1850       open (ibond,file=bondname,status='old',action='read')
1851 c      print *,"Processor",myrank," opened file IBOND"
1852       call getenv_loc('THETPAR',thetname)
1853       open (ithep,file=thetname,status='old',action='read')
1854 c      print *,"Processor",myrank," opened file ITHEP"
1855       call getenv_loc('ROTPAR',rotname)
1856       open (irotam,file=rotname,status='old',action='read')
1857 c      print *,"Processor",myrank," opened file IROTAM"
1858       call getenv_loc('TORPAR',torname)
1859       open (itorp,file=torname,status='old',action='read')
1860 c      print *,"Processor",myrank," opened file ITORP"
1861       call getenv_loc('TORDPAR',tordname)
1862       open (itordp,file=tordname,status='old',action='read')
1863 c      print *,"Processor",myrank," opened file ITORDP"
1864       call getenv_loc('SCCORPAR',sccorname)
1865       open (isccor,file=sccorname,status='old',action='read')
1866 c      print *,"Processor",myrank," opened file ISCCOR"
1867       call getenv_loc('FOURIER',fouriername)
1868       open (ifourier,file=fouriername,status='old',action='read')
1869 c      print *,"Processor",myrank," opened file IFOURIER"
1870       call getenv_loc('ELEPAR',elename)
1871       open (ielep,file=elename,status='old',action='read')
1872 c      print *,"Processor",myrank," opened file IELEP"
1873       call getenv_loc('SIDEPAR',sidename)
1874       open (isidep,file=sidename,status='old',action='read')
1875 c      print *,"Processor",myrank," opened file ISIDEP"
1876 c      print *,"Processor",myrank," opened parameter files"
1877 #elif (defined G77)
1878       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1879       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1880 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1881 C Get parameter filenames and open the parameter files.
1882       call getenv_loc('BONDPAR',bondname)
1883       open (ibond,file=bondname,status='old')
1884       call getenv_loc('THETPAR',thetname)
1885       open (ithep,file=thetname,status='old')
1886       call getenv_loc('ROTPAR',rotname)
1887       open (irotam,file=rotname,status='old')
1888       call getenv_loc('TORPAR',torname)
1889       open (itorp,file=torname,status='old')
1890       call getenv_loc('TORDPAR',tordname)
1891       open (itordp,file=tordname,status='old')
1892       call getenv_loc('SCCORPAR',sccorname)
1893       open (isccor,file=sccorname,status='old')
1894       call getenv_loc('FOURIER',fouriername)
1895       open (ifourier,file=fouriername,status='old')
1896       call getenv_loc('ELEPAR',elename)
1897       open (ielep,file=elename,status='old')
1898       call getenv_loc('SIDEPAR',sidename)
1899       open (isidep,file=sidename,status='old')
1900 #else
1901       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1902      &  readonly)
1903        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1904 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1905 C Get parameter filenames and open the parameter files.
1906       call getenv_loc('BONDPAR',bondname)
1907       open (ibond,file=bondname,status='old',readonly)
1908       call getenv_loc('THETPAR',thetname)
1909       open (ithep,file=thetname,status='old',readonly)
1910       call getenv_loc('ROTPAR',rotname)
1911       open (irotam,file=rotname,status='old',readonly)
1912       call getenv_loc('TORPAR',torname)
1913       open (itorp,file=torname,status='old',readonly)
1914       call getenv_loc('TORDPAR',tordname)
1915       open (itordp,file=tordname,status='old',readonly)
1916       call getenv_loc('SCCORPAR',sccorname)
1917       open (isccor,file=sccorname,status='old',readonly)
1918       call getenv_loc('FOURIER',fouriername)
1919       open (ifourier,file=fouriername,status='old',readonly)
1920       call getenv_loc('ELEPAR',elename)
1921       open (ielep,file=elename,status='old',readonly)
1922       call getenv_loc('SIDEPAR',sidename)
1923       open (isidep,file=sidename,status='old',readonly)
1924 #endif
1925 #ifndef OLDSCP
1926 C
1927 C 8/9/01 In the newest version SCp interaction constants are read from a file
1928 C Use -DOLDSCP to use hard-coded constants instead.
1929 C
1930       call getenv_loc('SCPPAR',scpname)
1931 #if defined(WINIFL) || defined(WINPGI)
1932       open (iscpp,file=scpname,status='old',readonly,shared)
1933 #elif (defined CRAY)  || (defined AIX)
1934       open (iscpp,file=scpname,status='old',action='read')
1935 #elif (defined G77)
1936       open (iscpp,file=scpname,status='old')
1937 #else
1938       open (iscpp,file=scpname,status='old',readonly)
1939 #endif
1940 #endif
1941       call getenv_loc('PATTERN',patname)
1942 #if defined(WINIFL) || defined(WINPGI)
1943       open (icbase,file=patname,status='old',readonly,shared)
1944 #elif (defined CRAY)  || (defined AIX)
1945       open (icbase,file=patname,status='old',action='read')
1946 #elif (defined G77)
1947       open (icbase,file=patname,status='old')
1948 #else
1949       open (icbase,file=patname,status='old',readonly)
1950 #endif
1951 #ifdef MPI
1952 C Open output file only for CG processes
1953 c      print *,"Processor",myrank," fg_rank",fg_rank
1954       if (fg_rank.eq.0) then
1955
1956       if (nodes.eq.1) then
1957         npos=3
1958       else
1959         npos = dlog10(dfloat(nodes-1))+1
1960       endif
1961       if (npos.lt.3) npos=3
1962       write (liczba,'(i1)') npos
1963       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1964      &  //')'
1965       write (liczba,form) me
1966       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1967      &  liczba(:ilen(liczba))
1968       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1969      &  //'.int'
1970       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1971      &  //'.pdb'
1972       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1973      &  liczba(:ilen(liczba))//'.mol2'
1974       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1975      &  liczba(:ilen(liczba))//'.stat'
1976       if (lentmp.gt.0)
1977      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1978      &      //liczba(:ilen(liczba))//'.stat')
1979       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1980      &  //'.rst'
1981       if(usampl) then
1982           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1983      & liczba(:ilen(liczba))//'.const'
1984       endif
1985
1986       endif
1987 #else
1988       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1989       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1990       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1991       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1992       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1993       if (lentmp.gt.0)
1994      &  call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
1995      &    //'.stat')
1996       rest2name=prefix(:ilen(prefix))//'.rst'
1997       if(usampl) then
1998          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1999       endif
2000 #endif
2001 #if defined(AIX) || defined(PGI)
2002       if (me.eq.king .or. .not. out1file)
2003      &   open(iout,file=outname,status='unknown')
2004 #ifdef DEBUG
2005       if (fg_rank.gt.0) then
2006         write (liczba,'(i3.3)') myrank/nfgtasks
2007         write (ll,'(bz,i3.3)') fg_rank
2008         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2009      &   status='unknown')
2010       endif
2011 #endif
2012       if(me.eq.king) then
2013        open(igeom,file=intname,status='unknown',position='append')
2014        open(ipdb,file=pdbname,status='unknown')
2015        open(imol2,file=mol2name,status='unknown')
2016        open(istat,file=statname,status='unknown',position='append')
2017       else
2018 c1out       open(iout,file=outname,status='unknown')
2019       endif
2020 #else
2021       if (me.eq.king .or. .not.out1file)
2022      &    open(iout,file=outname,status='unknown')
2023 #ifdef DEBUG
2024       if (fg_rank.gt.0) then
2025         write (liczba,'(i3.3)') myrank/nfgtasks
2026         write (ll,'(bz,i3.3)') fg_rank
2027         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2028      &   status='unknown')
2029       endif
2030 #endif
2031       if(me.eq.king) then
2032        open(igeom,file=intname,status='unknown',access='append')
2033        open(ipdb,file=pdbname,status='unknown')
2034        open(imol2,file=mol2name,status='unknown')
2035        open(istat,file=statname,status='unknown',access='append')
2036       else
2037 c1out       open(iout,file=outname,status='unknown')
2038       endif
2039 #endif
2040       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2041       csa_seed=prefix(:lenpre)//'.CSA.seed'
2042       csa_history=prefix(:lenpre)//'.CSA.history'
2043       csa_bank=prefix(:lenpre)//'.CSA.bank'
2044       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2045       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2046       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2047 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2048       csa_int=prefix(:lenpre)//'.int'
2049       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2050       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2051       csa_in=prefix(:lenpre)//'.CSA.in'
2052 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2053 C Write file names
2054       if (me.eq.king)then
2055       write (iout,'(80(1h-))')
2056       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2057       write (iout,'(80(1h-))')
2058       write (iout,*) "Input file                      : ",
2059      &  pref_orig(:ilen(pref_orig))//'.inp'
2060       write (iout,*) "Output file                     : ",
2061      &  outname(:ilen(outname))
2062       write (iout,*)
2063       write (iout,*) "Sidechain potential file        : ",
2064      &  sidename(:ilen(sidename))
2065 #ifndef OLDSCP
2066       write (iout,*) "SCp potential file              : ",
2067      &  scpname(:ilen(scpname))
2068 #endif
2069       write (iout,*) "Electrostatic potential file    : ",
2070      &  elename(:ilen(elename))
2071       write (iout,*) "Cumulant coefficient file       : ",
2072      &  fouriername(:ilen(fouriername))
2073       write (iout,*) "Torsional parameter file        : ",
2074      &  torname(:ilen(torname))
2075       write (iout,*) "Double torsional parameter file : ",
2076      &  tordname(:ilen(tordname))
2077       write (iout,*) "SCCOR parameter file : ",
2078      &  sccorname(:ilen(sccorname))
2079       write (iout,*) "Bond & inertia constant file    : ",
2080      &  bondname(:ilen(bondname))
2081       write (iout,*) "Bending parameter file          : ",
2082      &  thetname(:ilen(thetname))
2083       write (iout,*) "Rotamer parameter file          : ",
2084      &  rotname(:ilen(rotname))
2085       write (iout,*) "Threading database              : ",
2086      &  patname(:ilen(patname))
2087       if (lentmp.ne.0)
2088      &write (iout,*)" DIRTMP                          : ",
2089      &  tmpdir(:lentmp)
2090       write (iout,'(80(1h-))')
2091       endif
2092       return
2093       end
2094 c----------------------------------------------------------------------------
2095       subroutine card_concat(card)
2096       implicit real*8 (a-h,o-z)
2097       include 'DIMENSIONS'
2098       include 'COMMON.IOUNITS'
2099       character*(*) card
2100       character*80 karta,ucase
2101       external ilen
2102       read (inp,'(a)') karta
2103       karta=ucase(karta)
2104       card=' '
2105       do while (karta(80:80).eq.'&')
2106         card=card(:ilen(card)+1)//karta(:79)
2107         read (inp,'(a)') karta
2108         karta=ucase(karta)
2109       enddo
2110       card=card(:ilen(card)+1)//karta
2111       return
2112       end
2113 c----------------------------------------------------------------------------------
2114       subroutine readrst
2115       implicit real*8 (a-h,o-z)
2116       include 'DIMENSIONS'
2117       include 'COMMON.CHAIN'
2118       include 'COMMON.IOUNITS'
2119       include 'COMMON.MD'
2120       open(irest2,file=rest2name,status='unknown')
2121       read(irest2,*) totT,EK,potE,totE,t_bath
2122       do i=1,2*nres
2123          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2124       enddo
2125       do i=1,2*nres
2126          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2127       enddo
2128       if(usampl) then
2129              read (irest2,*) iset
2130       endif
2131       close(irest2)
2132       return
2133       end
2134 c---------------------------------------------------------------------------------
2135       subroutine read_fragments
2136       implicit real*8 (a-h,o-z)
2137       include 'DIMENSIONS'
2138 #ifdef MPI
2139       include 'mpif.h'
2140 #endif
2141       include 'COMMON.SETUP'
2142       include 'COMMON.CHAIN'
2143       include 'COMMON.IOUNITS'
2144       include 'COMMON.MD'
2145       include 'COMMON.CONTROL'
2146       read(inp,*) nset,nfrag,npair,nfrag_back
2147       if(me.eq.king.or..not.out1file)
2148      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2149      &  " nfrag_back",nfrag_back
2150       do iset=1,nset
2151          read(inp,*) mset(iset)
2152        do i=1,nfrag
2153          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2154      &     qinfrag(i,iset)
2155          if(me.eq.king.or..not.out1file)
2156      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2157      &     ifrag(2,i,iset), qinfrag(i,iset)
2158        enddo
2159        do i=1,npair
2160         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2161      &    qinpair(i,iset)
2162         if(me.eq.king.or..not.out1file)
2163      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2164      &    ipair(2,i,iset), qinpair(i,iset)
2165        enddo
2166        do i=1,nfrag_back
2167         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2168      &     wfrag_back(3,i,iset),
2169      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2170         if(me.eq.king.or..not.out1file)
2171      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2172      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2173        enddo
2174       enddo
2175       return
2176       end
2177 c-------------------------------------------------------------------------------
2178       subroutine read_dist_constr
2179       implicit real*8 (a-h,o-z)
2180       include 'DIMENSIONS'
2181 #ifdef MPI
2182       include 'mpif.h'
2183 #endif
2184       include 'COMMON.SETUP'
2185       include 'COMMON.CONTROL'
2186       include 'COMMON.CHAIN'
2187       include 'COMMON.IOUNITS'
2188       include 'COMMON.SBRIDGE'
2189       integer ifrag_(2,100),ipair_(2,100)
2190       double precision wfrag_(100),wpair_(100)
2191       character*500 controlcard
2192 c      write (iout,*) "Calling read_dist_constr"
2193 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2194 c      call flush(iout)
2195       call card_concat(controlcard)
2196       call readi(controlcard,"NFRAG",nfrag_,0)
2197       call readi(controlcard,"NPAIR",npair_,0)
2198       call readi(controlcard,"NDIST",ndist_,0)
2199       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2200       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2201       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2202       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2203       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2204 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2205 c      write (iout,*) "IFRAG"
2206 c      do i=1,nfrag_
2207 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2208 c      enddo
2209 c      write (iout,*) "IPAIR"
2210 c      do i=1,npair_
2211 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2212 c      enddo
2213       call flush(iout)
2214       do i=1,nfrag_
2215         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2216         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2217      &    ifrag_(2,i)=nstart_sup+nsup-1
2218 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2219         call flush(iout)
2220         if (wfrag_(i).gt.0.0d0) then
2221         do j=ifrag_(1,i),ifrag_(2,i)-1
2222           do k=j+1,ifrag_(2,i)
2223 c            write (iout,*) "j",j," k",k
2224             ddjk=dist(j,k)
2225             if (constr_dist.eq.1) then
2226               nhpb=nhpb+1
2227               ihpb(nhpb)=j
2228               jhpb(nhpb)=k
2229               dhpb(nhpb)=ddjk
2230               forcon(nhpb)=wfrag_(i)
2231             else if (constr_dist.eq.2) then
2232               if (ddjk.le.dist_cut) then
2233                 nhpb=nhpb+1
2234                 ihpb(nhpb)=j
2235                 jhpb(nhpb)=k
2236                 dhpb(nhpb)=ddjk
2237                 forcon(nhpb)=wfrag_(i)
2238               endif
2239             else
2240               nhpb=nhpb+1
2241               ihpb(nhpb)=j
2242               jhpb(nhpb)=k
2243               dhpb(nhpb)=ddjk
2244               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2245             endif
2246 #ifdef MPI
2247             if (.not.out1file .or. me.eq.king)
2248      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2249      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2250 #else
2251             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2252      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2253 #endif
2254           enddo
2255         enddo
2256         endif
2257       enddo
2258       do i=1,npair_
2259         if (wpair_(i).gt.0.0d0) then
2260         ii = ipair_(1,i)
2261         jj = ipair_(2,i)
2262         if (ii.gt.jj) then
2263           itemp=ii
2264           ii=jj
2265           jj=itemp
2266         endif
2267         do j=ifrag_(1,ii),ifrag_(2,ii)
2268           do k=ifrag_(1,jj),ifrag_(2,jj)
2269             nhpb=nhpb+1
2270             ihpb(nhpb)=j
2271             jhpb(nhpb)=k
2272             forcon(nhpb)=wpair_(i)
2273             dhpb(nhpb)=dist(j,k)
2274 #ifdef MPI
2275             if (.not.out1file .or. me.eq.king)
2276      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2277      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2278 #else
2279             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2280      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2281 #endif
2282           enddo
2283         enddo
2284         endif
2285       enddo
2286       do i=1,ndist_
2287         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2288         if (forcon(nhpb+1).gt.0.0d0) then
2289           nhpb=nhpb+1
2290           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2291 #ifdef MPI
2292           if (.not.out1file .or. me.eq.king)
2293      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2294      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2295 #else
2296           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2297      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2298 #endif
2299         endif
2300       enddo
2301       call flush(iout)
2302       return
2303       end
2304 c-------------------------------------------------------------------------------
2305 #ifdef WINIFL
2306       subroutine flush(iu)
2307       return
2308       end
2309 #endif
2310 #ifdef AIX
2311       subroutine flush(iu)
2312       call flush_(iu)
2313       return
2314       end
2315 #endif
2316 c------------------------------------------------------------------------------
2317       subroutine copy_to_tmp(source)
2318       include "DIMENSIONS"
2319       include "COMMON.IOUNITS"
2320       character*(*) source
2321       character* 256 tmpfile
2322       integer ilen
2323       external ilen
2324       logical ex
2325       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2326       inquire(file=tmpfile,exist=ex)
2327       if (ex) then
2328         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2329      &   " to temporary directory..."
2330         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2331         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2332       endif
2333       return
2334       end
2335 c------------------------------------------------------------------------------
2336       subroutine move_from_tmp(source)
2337       include "DIMENSIONS"
2338       include "COMMON.IOUNITS"
2339       character*(*) source
2340       integer ilen
2341       external ilen
2342       write (*,*) "Moving ",source(:ilen(source)),
2343      & " from temporary directory to working directory"
2344       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2345       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2346       return
2347       end
2348 c------------------------------------------------------------------------------
2349       subroutine random_init(seed)
2350 C
2351 C Initialize random number generator
2352 C
2353       implicit real*8 (a-h,o-z)
2354       include 'DIMENSIONS'
2355 #ifdef MPI
2356       include 'mpif.h'
2357       logical OKRandom, prng_restart
2358       real*8  r1
2359       integer iseed_array(4)
2360 #endif
2361       include 'COMMON.IOUNITS'
2362       include 'COMMON.TIME1'
2363       include 'COMMON.THREAD'
2364       include 'COMMON.SBRIDGE'
2365       include 'COMMON.CONTROL'
2366       include 'COMMON.MCM'
2367       include 'COMMON.MAP'
2368       include 'COMMON.HEADER'
2369       include 'COMMON.CSA'
2370       include 'COMMON.CHAIN'
2371       include 'COMMON.MUCA'
2372       include 'COMMON.MD'
2373       include 'COMMON.FFIELD'
2374       include 'COMMON.SETUP'
2375       iseed=-dint(dabs(seed))
2376       if (iseed.eq.0) then
2377         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2378      &    'Random seed undefined. The program will stop.'
2379         write (*,'(/80(1h*)/20x,a/80(1h*))')
2380      &    'Random seed undefined. The program will stop.'
2381 #ifdef MPI
2382         call mpi_finalize(mpi_comm_world,ierr)
2383 #endif
2384         stop 'Bad random seed.'
2385       endif
2386 #ifdef MPI
2387       if (fg_rank.eq.0) then
2388       seed=seed*(me+1)+1
2389 #ifdef AMD64
2390       if(me.eq.king)
2391      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2392       OKRandom = prng_restart(me,iseed)
2393 #else
2394       do i=1,4
2395        tmp=65536.0d0**(4-i)
2396        iseed_array(i) = dint(seed/tmp)
2397        seed=seed-iseed_array(i)*tmp
2398       enddo
2399       if(me.eq.king)
2400      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2401      &                 (iseed_array(i),i=1,4)
2402       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2403      &                 (iseed_array(i),i=1,4)
2404       OKRandom = prng_restart(me,iseed_array)
2405 #endif
2406       if (OKRandom) then
2407 c        r1 = prng_next(me)
2408          r1=ran_number(0.0D0,1.0D0)
2409         if(me.eq.king)
2410      &   write (iout,*) 'ran_num',r1
2411         if (r1.lt.0.0d0) OKRandom=.false.
2412       endif
2413       if (.not.OKRandom) then
2414         write (iout,*) 'PRNG IS NOT WORKING!!!'
2415         print *,'PRNG IS NOT WORKING!!!'
2416         if (me.eq.0) then
2417          call flush(iout)
2418          call mpi_abort(mpi_comm_world,error_msg,ierr)
2419          stop
2420         else
2421          write (iout,*) 'too many processors for parallel prng'
2422          write (*,*) 'too many processors for parallel prng'
2423          call flush(iout)
2424          stop
2425         endif
2426       endif
2427       endif
2428 #else
2429       call vrndst(iseed)
2430       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2431 #endif
2432       return
2433       end