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