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