source/unres/src_MD-M/readrtns_CSA.F

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