source/unres/src_MD-M/readrtns_CSA.F

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