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