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