source/unres/src_MD-M/readrtns_CSA.F

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