source/unres/src_MD-M/readrtns_CSA.F

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