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