source/unres/src_MD/readrtns.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 csa      if (modecalc.eq.8) then
  32 csa       inquire (file="fort.40",exist=file_exist)
  33 csa       if (.not.file_exist) call csaread
  34 csa      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,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
  48      &     ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
  49       enddo
  50 #ifdef MPI
  51       endif
  52 #endif
  53 c      print *,"Processor",myrank," leaves READRTNS"
  54       return
  55       end
  56 C-------------------------------------------------------------------------------
  57       subroutine read_control
  58 C
  59 C Read contorl data
  60 C
  61       implicit real*8 (a-h,o-z)
  62       include 'DIMENSIONS'
  63 #ifdef MP
  64       include 'mpif.h'
  65       logical OKRandom, prng_restart
  66       real*8  r1
  67 #endif
  68       include 'COMMON.IOUNITS'
  69       include 'COMMON.TIME1'
  70       include 'COMMON.THREAD'
  71       include 'COMMON.SBRIDGE'
  72       include 'COMMON.CONTROL'
  73       include 'COMMON.MCM'
  74       include 'COMMON.MAP'
  75       include 'COMMON.HEADER'
  76 csa      include 'COMMON.CSA'
  77       include 'COMMON.CHAIN'
  78       include 'COMMON.MUCA'
  79       include 'COMMON.MD'
  80       include 'COMMON.FFIELD'
  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,'CONSTR_HOMOL',constr_homology,0)
 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            write(*,*) "CSA not supported in this version"
 181            stop
 182 csa        modecalc=8
 183 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 184 crc
 185 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 186 crc
 187 crc        modecalc=9
 188 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 189 cfmc        modecalc=10
 190       else if (index(controlcard,'SOFTREG').gt.0) then
 191         modecalc=11
 192       else if (index(controlcard,'CHECK_BOND').gt.0) then
 193         modecalc=-1
 194       else if (index(controlcard,'TEST').gt.0) then
 195         modecalc=-2
 196       else if (index(controlcard,'MD').gt.0) then
 197         modecalc=12
 198       else if (index(controlcard,'RE ').gt.0) then
 199         modecalc=14
 200       endif
 201
 202       lmuca=index(controlcard,'MUCA').gt.0
 203       call readi(controlcard,'MUCADYN',mucadyn,0)
 204       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 205       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 206      & then
 207        write (iout,*) 'MUCADYN=',mucadyn
 208        write (iout,*) 'MUCASMOOTH=',muca_smooth
 209       endif
 210
 211       iscode=index(controlcard,'ONE_LETTER')
 212       indphi=index(controlcard,'PHI')
 213       indback=index(controlcard,'BACK')
 214       iranconf=index(controlcard,'RAND_CONF')
 215       i2ndstr=index(controlcard,'USE_SEC_PRED')
 216       gradout=index(controlcard,'GRADOUT').gt.0
 217       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 218
 219       if(me.eq.king.or..not.out1file)
 220      & write (iout,'(2a)') diagmeth(kdiag),
 221      &  ' routine used to diagonalize matrices.'
 222       return
 223       end
 224 c--------------------------------------------------------------------------
 225       subroutine read_REMDpar
 226 C
 227 C Read REMD settings
 228 C
 229       implicit real*8 (a-h,o-z)
 230       include 'DIMENSIONS'
 231       include 'COMMON.IOUNITS'
 232       include 'COMMON.TIME1'
 233       include 'COMMON.MD'
 234 #ifndef LANG0
 235       include 'COMMON.LANGEVIN'
 236 #else
 237       include 'COMMON.LANGEVIN.lang0'
 238 #endif
 239       include 'COMMON.INTERACT'
 240       include 'COMMON.NAMES'
 241       include 'COMMON.GEO'
 242       include 'COMMON.REMD'
 243       include 'COMMON.CONTROL'
 244       include 'COMMON.SETUP'
 245       character*80 ucase
 246       character*320 controlcard
 247       character*3200 controlcard1
 248       integer iremd_m_total
 249
 250       if(me.eq.king.or..not.out1file)
 251      & write (iout,*) "REMD setup"
 252
 253       call card_concat(controlcard)
 254       call readi(controlcard,"NREP",nrep,3)
 255       call readi(controlcard,"NSTEX",nstex,1000)
 256       call reada(controlcard,"RETMIN",retmin,10.0d0)
 257       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 258       mremdsync=(index(controlcard,'SYNC').gt.0)
 259       call readi(controlcard,"NSYN",i_sync_step,100)
 260       restart1file=(index(controlcard,'REST1FILE').gt.0)
 261       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 262       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 263       if(max_cache_traj_use.gt.max_cache_traj)
 264      &           max_cache_traj_use=max_cache_traj
 265       if(me.eq.king.or..not.out1file) then
 266 cd       if (traj1file) then
 267 crc caching is in testing - NTWX is not ignored
 268 cd        write (iout,*) "NTWX value is ignored"
 269 cd        write (iout,*) "  trajectory is stored to one file by master"
 270 cd        write (iout,*) "  before exchange at NSTEX intervals"
 271 cd       endif
 272        write (iout,*) "NREP= ",nrep
 273        write (iout,*) "NSTEX= ",nstex
 274        write (iout,*) "SYNC= ",mremdsync
 275        write (iout,*) "NSYN= ",i_sync_step
 276        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 277       endif
 278
 279       t_exchange_only=(index(controlcard,'TONLY').gt.0)
 280       call readi(controlcard,"HREMD",hremd,0)
 281       if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
 282         write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
 283       endif
 284       if(usampl.and.hremd.gt.0) then
 285             write (iout,'(//a)')
 286      &      "========== ERROR: USAMPL and HREMD cannot be used together"
 287 #ifdef MPI
 288             call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 289 #endif
 290             stop
 291       endif
 292
 293
 294       remd_tlist=.false.
 295       if (index(controlcard,'TLIST').gt.0) then
 296          remd_tlist=.true.
 297          call card_concat(controlcard1)
 298          read(controlcard1,*) (remd_t(i),i=1,nrep)
 299          if(me.eq.king.or..not.out1file)
 300      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 301       endif
 302       remd_mlist=.false.
 303       if (index(controlcard,'MLIST').gt.0) then
 304          remd_mlist=.true.
 305          call card_concat(controlcard1)
 306          read(controlcard1,*) (remd_m(i),i=1,nrep)
 307          if(me.eq.king.or..not.out1file) then
 308           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 309           iremd_m_total=0
 310           do i=1,nrep
 311            iremd_m_total=iremd_m_total+remd_m(i)
 312           enddo
 313           if(hremd.gt.1)then
 314            write (iout,*) 'Total number of replicas ',
 315      &       iremd_m_total*hremd
 316           else
 317            write (iout,*) 'Total number of replicas ',iremd_m_total
 318           endif
 319          endif
 320       endif
 321       if(me.eq.king.or..not.out1file)
 322      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 323       return
 324       end
 325 c--------------------------------------------------------------------------
 326       subroutine read_MDpar
 327 C
 328 C Read MD settings
 329 C
 330       implicit real*8 (a-h,o-z)
 331       include 'DIMENSIONS'
 332       include 'COMMON.IOUNITS'
 333       include 'COMMON.TIME1'
 334       include 'COMMON.MD'
 335 #ifndef LANG0
 336       include 'COMMON.LANGEVIN'
 337 #else
 338       include 'COMMON.LANGEVIN.lang0'
 339 #endif
 340       include 'COMMON.INTERACT'
 341       include 'COMMON.NAMES'
 342       include 'COMMON.GEO'
 343       include 'COMMON.SETUP'
 344       include 'COMMON.CONTROL'
 345       include 'COMMON.SPLITELE'
 346       character*80 ucase
 347       character*320 controlcard
 348
 349       call card_concat(controlcard)
 350       call readi(controlcard,"NSTEP",n_timestep,1000000)
 351       call readi(controlcard,"NTWE",ntwe,100)
 352       call readi(controlcard,"NTWX",ntwx,1000)
 353       call reada(controlcard,"DT",d_time,1.0d-1)
 354       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 355       call reada(controlcard,"DAMAX",damax,1.0d1)
 356       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 357       call readi(controlcard,"LANG",lang,0)
 358       RESPA = index(controlcard,"RESPA") .gt. 0
 359       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 360       ntime_split0=ntime_split
 361       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 362       ntime_split0=ntime_split
 363       call reada(controlcard,"R_CUT",r_cut,2.0d0)
 364       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 365       rest = index(controlcard,"REST").gt.0
 366       tbf = index(controlcard,"TBF").gt.0
 367       call readi(controlcard,"HMC",hmc,0)
 368       tnp = index(controlcard,"NOSEPOINCARE99").gt.0
 369       tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
 370       tnh = index(controlcard,"NOSEHOOVER96").gt.0
 371       if (RESPA.and.tnh)then
 372         xiresp = index(controlcard,"XIRESP").gt.0
 373       endif
 374       call reada(controlcard,"Q_NP",Q_np,0.1d0)
 375       usampl = index(controlcard,"USAMPL").gt.0
 376       mdpdb = index(controlcard,"MDPDB").gt.0
 377       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 378       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 379       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 380       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 381       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 382       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 383       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 384       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 385       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 386       large = index(controlcard,"LARGE").gt.0
 387       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 388       rattle = index(controlcard,"RATTLE").gt.0
 389 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 390       nset=0
 391       if(usampl) then
 392         call read_fragments
 393       endif
 394
 395       if(me.eq.king.or..not.out1file) then
 396        write (iout,*)
 397        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 398        write (iout,*)
 399        write (iout,'(a)') "The units are:"
 400        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 401        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 402      &  " acceleration: angstrom/(48.9 fs)**2"
 403        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 404        write (iout,*)
 405        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 406        write (iout,'(a60,f10.5,a)')
 407      &  "Initial time step of numerical integration:",d_time,
 408      &  " natural units"
 409        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 410        if (RESPA) then
 411         write (iout,'(2a,i4,a)')
 412      &    "A-MTS algorithm used; initial time step for fast-varying",
 413      &    " short-range forces split into",ntime_split," steps."
 414         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 415      &   r_cut," lambda",rlamb
 416        endif
 417        write (iout,'(2a,f10.5)')
 418      &  "Maximum acceleration threshold to reduce the time step",
 419      &  "/increase split number:",damax
 420        write (iout,'(2a,f10.5)')
 421      &  "Maximum predicted energy drift to reduce the timestep",
 422      &  "/increase split number:",edriftmax
 423        write (iout,'(a60,f10.5)')
 424      & "Maximum velocity threshold to reduce velocities:",dvmax
 425        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 426        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 427        if (rattle) write (iout,'(a60)')
 428      &  "Rattle algorithm used to constrain the virtual bonds"
 429       endif
 430       reset_fricmat=1000
 431       if (lang.gt.0) then
 432         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 433         call reada(controlcard,"RWAT",rwat,1.4d0)
 434         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 435         surfarea=index(controlcard,"SURFAREA").gt.0
 436         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 437         if(me.eq.king.or..not.out1file)then
 438          write (iout,'(/a,$)') "Langevin dynamics calculation"
 439          if (lang.eq.1) then
 440           write (iout,'(a/)')
 441      &      " with direct integration of Langevin equations"
 442          else if (lang.eq.2) then
 443           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 444          else if (lang.eq.3) then
 445           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 446          else if (lang.eq.4) then
 447           write (iout,'(a/)') " in overdamped mode"
 448          else
 449           write (iout,'(//a,i5)')
 450      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 451           stop
 452          endif
 453          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 454          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 455          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 456          write (iout,'(a60,f10.5)')
 457      &   "Scaling factor of the friction forces:",scal_fric
 458          if (surfarea) write (iout,'(2a,i10,a)')
 459      &     "Friction coefficients will be scaled by solvent-accessible",
 460      &     " surface area every",reset_fricmat," steps."
 461         endif
 462 c Calculate friction coefficients and bounds of stochastic forces
 463         eta=6*pi*cPoise*etawat
 464         if(me.eq.king.or..not.out1file)
 465      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 466      &   ,eta
 467         gamp=scal_fric*(pstok+rwat)*eta
 468         stdfp=dsqrt(2*Rb*t_bath/d_time)
 469         do i=1,ntyp
 470           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 471           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 472         enddo
 473         if(me.eq.king.or..not.out1file)then
 474          write (iout,'(/2a/)')
 475      &   "Radii of site types and friction coefficients and std's of",
 476      &   " stochastic forces of fully exposed sites"
 477          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 478          do i=1,ntyp
 479           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 480      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 481          enddo
 482         endif
 483       else if (tbf) then
 484         if(me.eq.king.or..not.out1file)then
 485          write (iout,'(a)') "Berendsen bath calculation"
 486          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 487          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 488          if (reset_moment)
 489      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 490      &   count_reset_moment," steps"
 491          if (reset_vel)
 492      &    write (iout,'(a,i10,a)')
 493      &    "Velocities will be reset at random every",count_reset_vel,
 494      &   " steps"
 495         endif
 496       else if (tnp .or. tnp1 .or. tnh) then
 497         if (tnp .or. tnp1) then
 498            write (iout,'(a)') "Nose-Poincare bath calculation"
 499            if (tnp) write (iout,'(a)')
 500      & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
 501            if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
 502         else
 503            write (iout,'(a)') "Nose-Hoover bath calculation"
 504            write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
 505               nresn=1
 506               nyosh=1
 507               nnos=1
 508               do i=1,nnos
 509                qmass(i)=Q_np
 510                xlogs(i)=1.0
 511                vlogs(i)=0.0
 512               enddo
 513               do i=1,nyosh
 514                WDTI(i) = 1.0*d_time/nresn
 515                WDTI2(i)=WDTI(i)/2
 516                WDTI4(i)=WDTI(i)/4
 517                WDTI8(i)=WDTI(i)/8
 518               enddo
 519               if (RESPA) then
 520                if(xiresp) then
 521                  write (iout,'(a)') "NVT-XI-RESPA algorithm"
 522                else
 523                  write (iout,'(a)') "NVT-XO-RESPA algorithm"
 524                endif
 525                do i=1,nyosh
 526                 WDTIi(i) = 1.0*d_time/nresn/ntime_split
 527                 WDTIi2(i)=WDTIi(i)/2
 528                 WDTIi4(i)=WDTIi(i)/4
 529                 WDTIi8(i)=WDTIi(i)/8
 530                enddo
 531               endif
 532         endif
 533
 534         write (iout,'(a60,f10.5)') "Temperature:",t_bath
 535         write (iout,'(a60,f10.5)') "Q =",Q_np
 536         if (reset_moment)
 537      &  write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 538      &   count_reset_moment," steps"
 539         if (reset_vel)
 540      &    write (iout,'(a,i10,a)')
 541      &    "Velocities will be reset at random every",count_reset_vel,
 542      &   " steps"
 543
 544       else if (hmc.gt.0) then
 545          write (iout,'(a)') "Hybrid Monte Carlo calculation"
 546          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 547          write (iout,'(a60,i10)')
 548      &         "Number of MD steps between Metropolis tests:",hmc
 549
 550       else
 551         if(me.eq.king.or..not.out1file)
 552      &   write (iout,'(a31)') "Microcanonical mode calculation"
 553       endif
 554       if(me.eq.king.or..not.out1file)then
 555        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 556        if (usampl) then
 557           write(iout,*) "MD running with constraints."
 558           write(iout,*) "Equilibration time ", eq_time, " mtus."
 559           write(iout,*) "Constraining ", nfrag," fragments."
 560           write(iout,*) "Length of each fragment, weight and q0:"
 561           do iset=1,nset
 562            write (iout,*) "Set of restraints #",iset
 563            do i=1,nfrag
 564               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 565      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 566            enddo
 567            write(iout,*) "constraints between ", npair, "fragments."
 568            write(iout,*) "constraint pairs, weights and q0:"
 569            do i=1,npair
 570             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 571      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 572            enddo
 573            write(iout,*) "angle constraints within ", nfrag_back,
 574      &      "backbone fragments."
 575            write(iout,*) "fragment, weights:"
 576            do i=1,nfrag_back
 577             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 578      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 579      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 580            enddo
 581           enddo
 582         iset=mod(kolor,nset)+1
 583        endif
 584       endif
 585       if(me.eq.king.or..not.out1file)
 586      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 587       return
 588       end
 589 c------------------------------------------------------------------------------
 590       subroutine molread
 591 C
 592 C Read molecular data.
 593 C
 594       implicit real*8 (a-h,o-z)
 595       include 'DIMENSIONS'
 596 #ifdef MPI
 597       include 'mpif.h'
 598       integer error_msg
 599 #endif
 600       include 'COMMON.IOUNITS'
 601       include 'COMMON.GEO'
 602       include 'COMMON.VAR'
 603       include 'COMMON.INTERACT'
 604       include 'COMMON.LOCAL'
 605       include 'COMMON.NAMES'
 606       include 'COMMON.CHAIN'
 607       include 'COMMON.FFIELD'
 608       include 'COMMON.SBRIDGE'
 609       include 'COMMON.HEADER'
 610       include 'COMMON.CONTROL'
 611       include 'COMMON.DBASE'
 612       include 'COMMON.THREAD'
 613       include 'COMMON.CONTACTS'
 614       include 'COMMON.TORCNSTR'
 615       include 'COMMON.TIME1'
 616       include 'COMMON.BOUNDS'
 617       include 'COMMON.MD'
 618       include 'COMMON.REMD'
 619       include 'COMMON.SETUP'
 620       character*4 sequence(maxres)
 621       integer rescode
 622       double precision x(maxvar)
 623       character*256 pdbfile
 624       character*320 weightcard
 625       character*80 weightcard_t,ucase
 626       dimension itype_pdb(maxres)
 627       common /pizda/ itype_pdb
 628       logical seq_comp,fail
 629       double precision energia(0:n_ene)
 630
 631       integer ilen
 632       external ilen
 633 C
 634 C Body
 635 C
 636 C Read weights of the subsequent energy terms.
 637       if(hremd.gt.0) then
 638
 639        k=0
 640        do il=1,hremd
 641         do i=1,nrep
 642          do j=1,remd_m(i)
 643           i2set(k)=il
 644           k=k+1
 645          enddo
 646         enddo
 647        enddo
 648
 649        if(me.eq.king.or..not.out1file) then
 650         write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
 651         write (iout,*) 'Current weights for processor ',
 652      &                 me,' set ',i2set(me)
 653        endif
 654
 655        do i=1,hremd
 656          call card_concat(weightcard)
 657          call reada(weightcard,'WLONG',wlong,1.0D0)
 658          call reada(weightcard,'WSC',wsc,wlong)
 659          call reada(weightcard,'WSCP',wscp,wlong)
 660          call reada(weightcard,'WELEC',welec,1.0D0)
 661          call reada(weightcard,'WVDWPP',wvdwpp,welec)
 662          call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 663          call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 664          call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 665          call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 666          call reada(weightcard,'WTURN3',wturn3,1.0D0)
 667          call reada(weightcard,'WTURN4',wturn4,1.0D0)
 668          call reada(weightcard,'WTURN6',wturn6,1.0D0)
 669          call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 670          call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 671          call reada(weightcard,'WBOND',wbond,1.0D0)
 672          call reada(weightcard,'WTOR',wtor,1.0D0)
 673          call reada(weightcard,'WTORD',wtor_d,1.0D0)
 674          call reada(weightcard,'WANG',wang,1.0D0)
 675          call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 676          call reada(weightcard,'SCAL14',scal14,0.4D0)
 677          call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 678          call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 679          call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 680          call reada(weightcard,'TEMP0',temp0,300.0d0)
 681          if (index(weightcard,'SOFT').gt.0) ipot=6
 682 C 12/1/95 Added weight for the multi-body term WCORR
 683          call reada(weightcard,'WCORRH',wcorr,1.0D0)
 684          if (wcorr4.gt.0.0d0) wcorr=wcorr4
 685
 686          hweights(i,1)=wsc
 687          hweights(i,2)=wscp
 688          hweights(i,3)=welec
 689          hweights(i,4)=wcorr
 690          hweights(i,5)=wcorr5
 691          hweights(i,6)=wcorr6
 692          hweights(i,7)=wel_loc
 693          hweights(i,8)=wturn3
 694          hweights(i,9)=wturn4
 695          hweights(i,10)=wturn6
 696          hweights(i,11)=wang
 697          hweights(i,12)=wscloc
 698          hweights(i,13)=wtor
 699          hweights(i,14)=wtor_d
 700          hweights(i,15)=wstrain
 701          hweights(i,16)=wvdwpp
 702          hweights(i,17)=wbond
 703          hweights(i,18)=scal14
 704          hweights(i,21)=wsccor
 705
 706        enddo
 707
 708        do i=1,n_ene
 709          weights(i)=hweights(i2set(me),i)
 710        enddo
 711        wsc    =weights(1)
 712        wscp   =weights(2)
 713        welec  =weights(3)
 714        wcorr  =weights(4)
 715        wcorr5 =weights(5)
 716        wcorr6 =weights(6)
 717        wel_loc=weights(7)
 718        wturn3 =weights(8)
 719        wturn4 =weights(9)
 720        wturn6 =weights(10)
 721        wang   =weights(11)
 722        wscloc =weights(12)
 723        wtor   =weights(13)
 724        wtor_d =weights(14)
 725        wstrain=weights(15)
 726        wvdwpp =weights(16)
 727        wbond  =weights(17)
 728        scal14 =weights(18)
 729        wsccor =weights(21)
 730
 731
 732       else
 733        call card_concat(weightcard)
 734        call reada(weightcard,'WLONG',wlong,1.0D0)
 735        call reada(weightcard,'WSC',wsc,wlong)
 736        call reada(weightcard,'WSCP',wscp,wlong)
 737        call reada(weightcard,'WELEC',welec,1.0D0)
 738        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 739        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 740        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 741        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 742        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 743        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 744        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 745        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 746        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 747        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 748        call reada(weightcard,'WBOND',wbond,1.0D0)
 749        call reada(weightcard,'WTOR',wtor,1.0D0)
 750        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 751        call reada(weightcard,'WANG',wang,1.0D0)
 752        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 753        call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
 754        call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
 755        call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
 756        call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
 757        call reada(weightcard,'SCAL14',scal14,0.4D0)
 758        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 759        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 760        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 761        call reada(weightcard,'TEMP0',temp0,300.0d0)
 762        if (index(weightcard,'SOFT').gt.0) ipot=6
 763 C 12/1/95 Added weight for the multi-body term WCORR
 764        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 765        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 766        weights(1)=wsc
 767        weights(2)=wscp
 768        weights(3)=welec
 769        weights(4)=wcorr
 770        weights(5)=wcorr5
 771        weights(6)=wcorr6
 772        weights(7)=wel_loc
 773        weights(8)=wturn3
 774        weights(9)=wturn4
 775        weights(10)=wturn6
 776        weights(11)=wang
 777        weights(12)=wscloc
 778        weights(13)=wtor
 779        weights(14)=wtor_d
 780        weights(15)=wstrain
 781        weights(16)=wvdwpp
 782        weights(17)=wbond
 783        weights(18)=scal14
 784        weights(21)=wsccor
 785       endif
 786        weights(25)=wdfa_dist
 787        weights(26)=wdfa_tor
 788        weights(27)=wdfa_nei
 789        weights(28)=wdfa_beta
 790
 791       if(me.eq.king.or..not.out1file)
 792      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 793      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 794      &  wturn4,wturn6,
 795      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 796    10 format (/'Energy-term weights (unscaled):'//
 797      & 'WSCC=   ',f10.6,' (SC-SC)'/
 798      & 'WSCP=   ',f10.6,' (SC-p)'/
 799      & 'WELEC=  ',f10.6,' (p-p electr)'/
 800      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 801      & 'WBOND=  ',f10.6,' (stretching)'/
 802      & 'WANG=   ',f10.6,' (bending)'/
 803      & 'WSCLOC= ',f10.6,' (SC local)'/
 804      & 'WTOR=   ',f10.6,' (torsional)'/
 805      & 'WTORD=  ',f10.6,' (double torsional)'/
 806      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 807      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 808      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 809      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 810      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 811      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 812      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 813      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 814      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 815      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 816      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 817      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 818      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 819       if(me.eq.king.or..not.out1file)then
 820        if (wcorr4.gt.0.0d0) then
 821         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 822      &   'between contact pairs of peptide groups'
 823         write (iout,'(2(a,f5.3/))')
 824      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 825      &  'Range of quenching the correlation terms:',2*delt_corr
 826        else if (wcorr.gt.0.0d0) then
 827         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 828      &   'between contact pairs of peptide groups'
 829        endif
 830        write (iout,'(a,f8.3)')
 831      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 832        write (iout,'(a,f8.3)')
 833      &  'General scaling factor of SC-p interactions:',scalscp
 834       endif
 835       r0_corr=cutoff_corr-delt_corr
 836       do i=1,20
 837         aad(i,1)=scalscp*aad(i,1)
 838         aad(i,2)=scalscp*aad(i,2)
 839         bad(i,1)=scalscp*bad(i,1)
 840         bad(i,2)=scalscp*bad(i,2)
 841       enddo
 842       call rescale_weights(t_bath)
 843       if(me.eq.king.or..not.out1file)
 844      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 845      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 846      &  wturn4,wturn6,
 847      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 848    22 format (/'Energy-term weights (scaled):'//
 849      & 'WSCC=   ',f10.6,' (SC-SC)'/
 850      & 'WSCP=   ',f10.6,' (SC-p)'/
 851      & 'WELEC=  ',f10.6,' (p-p electr)'/
 852      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 853      & 'WBOND=  ',f10.6,' (stretching)'/
 854      & 'WANG=   ',f10.6,' (bending)'/
 855      & 'WSCLOC= ',f10.6,' (SC local)'/
 856      & 'WTOR=   ',f10.6,' (torsional)'/
 857      & 'WTORD=  ',f10.6,' (double torsional)'/
 858      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 859      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 860      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 861      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 862      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 863      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 864      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 865      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 866      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 867      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 868      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 869      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 870      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 871       if(me.eq.king.or..not.out1file)
 872      & write (iout,*) "Reference temperature for weights calculation:",
 873      &  temp0
 874       call reada(weightcard,"D0CM",d0cm,3.78d0)
 875       call reada(weightcard,"AKCM",akcm,15.1d0)
 876       call reada(weightcard,"AKTH",akth,11.0d0)
 877       call reada(weightcard,"AKCT",akct,12.0d0)
 878       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 879       call reada(weightcard,"V2SS",v2ss,7.61d0)
 880       call reada(weightcard,"V3SS",v3ss,13.7d0)
 881       call reada(weightcard,"EBR",ebr,-5.50D0)
 882       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 883       do i=1,maxres
 884         dyn_ss_mask(i)=.false.
 885       enddo
 886       do i=1,maxres-1
 887         do j=i+1,maxres
 888           dyn_ssbond_ij(i,j)=1.0d300
 889         enddo
 890       enddo
 891       call reada(weightcard,"HT",Ht,0.0D0)
 892       if (dyn_ss) then
 893         ss_depth=ebr/wsc-0.25*eps(1,1)
 894         Ht=Ht/wsc-0.25*eps(1,1)
 895         akcm=akcm*wstrain/wsc
 896         akth=akth*wstrain/wsc
 897         akct=akct*wstrain/wsc
 898         v1ss=v1ss*wstrain/wsc
 899         v2ss=v2ss*wstrain/wsc
 900         v3ss=v3ss*wstrain/wsc
 901       else
 902         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 903       endif
 904
 905       if(me.eq.king.or..not.out1file) then
 906        write (iout,*) "Parameters of the SS-bond potential:"
 907        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 908      & " AKCT",akct
 909        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 910        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 911        write (iout,*)" HT",Ht
 912        print *,'indpdb=',indpdb,' pdbref=',pdbref
 913       endif
 914       if (indpdb.gt.0 .or. pdbref) then
 915         read(inp,'(a)') pdbfile
 916         if(me.eq.king.or..not.out1file)
 917      &   write (iout,'(2a)') 'PDB data will be read from file ',
 918      &   pdbfile(:ilen(pdbfile))
 919         open(ipdbin,file=pdbfile,status='old',err=33)
 920         goto 34
 921   33    write (iout,'(a)') 'Error opening PDB file.'
 922         stop
 923   34    continue
 924 c        print *,'Begin reading pdb data'
 925         call readpdb
 926         do i=1,2*nres
 927           do j=1,3
 928             crefjlee(j,i)=c(j,i)
 929           enddo
 930         enddo
 931 c        print *,'Finished reading pdb data'
 932         if(me.eq.king.or..not.out1file)
 933      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 934      &   ' nstart_sup=',nstart_sup
 935         do i=1,nres
 936           itype_pdb(i)=itype(i)
 937         enddo
 938         close (ipdbin)
 939         nnt=nstart_sup
 940         nct=nstart_sup+nsup-1
 941         call contact(.false.,ncont_ref,icont_ref,co)
 942
 943         if (sideadd) then
 944 C Following 2 lines for diagnostics; comment out if not needed
 945          write (iout,*) "Before sideadd"
 946          call intout
 947          if(me.eq.king.or..not.out1file)
 948      &    write(iout,*)'Adding sidechains'
 949          maxsi=1000
 950          do i=2,nres-1
 951           iti=itype(i)
 952           if (iti.ne.10) then
 953             nsi=0
 954             fail=.true.
 955             do while (fail.and.nsi.le.maxsi)
 956               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 957               nsi=nsi+1
 958             enddo
 959             if(fail) write(iout,*)'Adding sidechain failed for res ',
 960      &              i,' after ',nsi,' trials'
 961           endif
 962          enddo
 963 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
 964          call chainbuild
 965         endif
 966 C Following 2 lines for diagnostics; comment out if not needed
 967 c        write (iout,*) "After sideadd"
 968 c        call intout
 969       endif
 970       if (indpdb.eq.0) then
 971 C Read sequence if not taken from the pdb file.
 972         read (inp,*) nres
 973 c        print *,'nres=',nres
 974         if (iscode.gt.0) then
 975           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 976         else
 977           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 978         endif
 979 C Convert sequence to numeric code
 980         do i=1,nres
 981           itype(i)=rescode(i,sequence(i),iscode)
 982         enddo
 983 C Assign initial virtual bond lengths
 984         do i=2,nres
 985           vbld(i)=vbl
 986           vbld_inv(i)=vblinv
 987         enddo
 988         do i=2,nres-1
 989           vbld(i+nres)=dsc(itype(i))
 990           vbld_inv(i+nres)=dsc_inv(itype(i))
 991 c          write (iout,*) "i",i," itype",itype(i),
 992 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 993         enddo
 994       endif
 995 c      print *,nres
 996 c      print '(20i4)',(itype(i),i=1,nres)
 997       do i=1,nres
 998 #ifdef PROCOR
 999         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
1000 #else
1001         if (itype(i).eq.21) then
1002 #endif
1003           itel(i)=0
1004 #ifdef PROCOR
1005         else if (itype(i+1).ne.20) then
1006 #else
1007         else if (itype(i).ne.20) then
1008 #endif
1009           itel(i)=1
1010         else
1011           itel(i)=2
1012         endif
1013       enddo
1014       if(me.eq.king.or..not.out1file)then
1015        write (iout,*) "ITEL"
1016        do i=1,nres-1
1017          write (iout,*) i,itype(i),itel(i)
1018        enddo
1019        print *,'Call Read_Bridge.'
1020       endif
1021       call read_bridge
1022 C 8/13/98 Set limits to generating the dihedral angles
1023       do i=1,nres
1024         phibound(1,i)=-pi
1025         phibound(2,i)=pi
1026       enddo
1027       read (inp,*) ndih_constr
1028       if (ndih_constr.gt.0) then
1029         read (inp,*) ftors
1030         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1031         if(me.eq.king.or..not.out1file)then
1032          write (iout,*)
1033      &   'There are',ndih_constr,' constraints on phi angles.'
1034          do i=1,ndih_constr
1035           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1036          enddo
1037         endif
1038         do i=1,ndih_constr
1039           phi0(i)=deg2rad*phi0(i)
1040           drange(i)=deg2rad*drange(i)
1041         enddo
1042         if(me.eq.king.or..not.out1file)
1043      &   write (iout,*) 'FTORS',ftors
1044         do i=1,ndih_constr
1045           ii = idih_constr(i)
1046           phibound(1,ii) = phi0(i)-drange(i)
1047           phibound(2,ii) = phi0(i)+drange(i)
1048         enddo
1049       endif
1050       nnt=1
1051 #ifdef MPI
1052       if (me.eq.king) then
1053 #endif
1054        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1055        do i=1,nres
1056          write (iout,'(a3,i5,2f10.1)')
1057      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1058        enddo
1059 #ifdef MP
1060       endif
1061 #endif
1062       nct=nres
1063 cd      print *,'NNT=',NNT,' NCT=',NCT
1064       if (itype(1).eq.21) nnt=2
1065       if (itype(nres).eq.21) nct=nct-1
1066
1067 C     Bartek:READ init_vars
1068 C     Initialize variables!
1069 C     Juyong:READ read_info
1070 C     READ fragment information!!
1071 C     both routines should be in dfa.F file!!
1072
1073       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1074      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1075        call init_dfa_vars
1076        print*, 'init_dfa_vars finished!'
1077        call read_dfa_info
1078        print*, 'read_dfa_info finished!'
1079       endif
1080 C
1081       if (pdbref) then
1082         if(me.eq.king.or..not.out1file)
1083      &   write (iout,'(a,i3)') 'nsup=',nsup
1084         nstart_seq=nnt
1085         if (nsup.le.(nct-nnt+1)) then
1086           do i=0,nct-nnt+1-nsup
1087             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1088               nstart_seq=nnt+i
1089               goto 111
1090             endif
1091           enddo
1092           write (iout,'(a)')
1093      &            'Error - sequences to be superposed do not match.'
1094           stop
1095         else
1096           do i=0,nsup-(nct-nnt+1)
1097             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1098      &      then
1099               nstart_sup=nstart_sup+i
1100               nsup=nct-nnt+1
1101               goto 111
1102             endif
1103           enddo
1104           write (iout,'(a)')
1105      &            'Error - sequences to be superposed do not match.'
1106         endif
1107   111   continue
1108         if (nsup.eq.0) nsup=nct-nnt
1109         if (nstart_sup.eq.0) nstart_sup=nnt
1110         if (nstart_seq.eq.0) nstart_seq=nnt
1111         if(me.eq.king.or..not.out1file)
1112      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1113      &                 ' nstart_seq=',nstart_seq
1114       endif
1115 c--- Zscore rms -------
1116       if (nz_start.eq.0) nz_start=nnt
1117       if (nz_end.eq.0 .and. nsup.gt.0) then
1118         nz_end=nnt+nsup-1
1119       else if (nz_end.eq.0) then
1120         nz_end=nct
1121       endif
1122       if(me.eq.king.or..not.out1file)then
1123        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1124        write (iout,*) 'IZ_SC=',iz_sc
1125       endif
1126 c----------------------
1127       call init_int_table
1128       if (refstr) then
1129         if (.not.pdbref) then
1130           call read_angles(inp,*38)
1131           goto 39
1132    38     write (iout,'(a)') 'Error reading reference structure.'
1133 #ifdef MPI
1134           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1135           stop 'Error reading reference structure'
1136 #endif
1137    39     call chainbuild
1138           call setup_var
1139 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1140           nstart_sup=nnt
1141           nstart_seq=nnt
1142           nsup=nct-nnt+1
1143           do i=1,2*nres
1144             do j=1,3
1145               cref(j,i)=c(j,i)
1146             enddo
1147           enddo
1148           call contact(.true.,ncont_ref,icont_ref,co)
1149         endif
1150         if(me.eq.king.or..not.out1file)
1151      &   write (iout,*) 'Contact order:',co
1152         if (pdbref) then
1153         if(me.eq.king.or..not.out1file)
1154      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1155         do i=1,ncont_ref
1156           do j=1,2
1157             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1158           enddo
1159           if(me.eq.king.or..not.out1file)
1160      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1161      &     icont_ref(1,i),' ',
1162      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1163         enddo
1164         endif
1165       endif
1166 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1167       if (constr_dist.gt.0) then
1168         call read_dist_constr
1169       endif
1170
1171
1172       if (constr_homology.gt.0) then
1173         call read_constr_homology
1174         if (indpdb.gt.0 .or. pdbref) then
1175           do i=1,2*nres
1176             do j=1,3
1177               c(j,i)=crefjlee(j,i)
1178               cref(j,i)=crefjlee(j,i)
1179             enddo
1180           enddo
1181         endif
1182       else
1183         homol_nset=0
1184       endif
1185
1186
1187       if (nhpb.gt.0) call hpb_partition
1188 c      write (iout,*) "After read_dist_constr nhpb",nhpb
1189 c      call flush(iout)
1190       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1191      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1192      &    modecalc.ne.10) then
1193 C If input structure hasn't been supplied from the PDB file read or generate
1194 C initial geometry.
1195         if (iranconf.eq.0 .and. .not. extconf) then
1196           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1197      &     write (iout,'(a)') 'Initial geometry will be read in.'
1198           if (read_cart) then
1199             read(inp,'(8f10.5)',end=36,err=36)
1200      &       ((c(l,k),l=1,3),k=1,nres),
1201      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1202             call int_from_cart1(.false.)
1203             do i=1,nres-1
1204               do j=1,3
1205                 dc(j,i)=c(j,i+1)-c(j,i)
1206                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1207               enddo
1208             enddo
1209             do i=nnt,nct
1210               if (itype(i).ne.10) then
1211                 do j=1,3
1212                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1213                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1214                 enddo
1215               endif
1216             enddo
1217             return
1218           else
1219             call read_angles(inp,*36)
1220           endif
1221           goto 37
1222    36     write (iout,'(a)') 'Error reading angle file.'
1223 #ifdef MPI
1224           call mpi_finalize( MPI_COMM_WORLD,IERR )
1225 #endif
1226           stop 'Error reading angle file.'
1227    37     continue
1228         else if (extconf) then
1229          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1230      &    write (iout,'(a)') 'Extended chain initial geometry.'
1231          do i=3,nres
1232           theta(i)=90d0*deg2rad
1233          enddo
1234          do i=4,nres
1235           phi(i)=180d0*deg2rad
1236          enddo
1237          do i=2,nres-1
1238           alph(i)=110d0*deg2rad
1239          enddo
1240          do i=2,nres-1
1241           omeg(i)=-120d0*deg2rad
1242          enddo
1243         else
1244           if(me.eq.king.or..not.out1file)
1245      &     write (iout,'(a)') 'Random-generated initial geometry.'
1246
1247
1248 #ifdef MPI
1249           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1250      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1251 #endif
1252             do itrial=1,100
1253               itmp=1
1254               call gen_rand_conf(itmp,*30)
1255               goto 40
1256    30         write (iout,*) 'Failed to generate random conformation',
1257      &          ', itrial=',itrial
1258               write (*,*) 'Processor:',me,
1259      &          ' Failed to generate random conformation',
1260      &          ' itrial=',itrial
1261               call intout
1262
1263 #ifdef AIX
1264               call flush_(iout)
1265 #else
1266               call flush(iout)
1267 #endif
1268             enddo
1269             write (iout,'(a,i3,a)') 'Processor:',me,
1270      &        ' error in generating random conformation.'
1271             write (*,'(a,i3,a)') 'Processor:',me,
1272      &        ' error in generating random conformation.'
1273             call flush(iout)
1274 #ifdef MPI
1275             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1276    40       continue
1277           endif
1278 #else
1279    40     continue
1280 #endif
1281         endif
1282       elseif (modecalc.eq.4) then
1283         read (inp,'(a)') intinname
1284         open (intin,file=intinname,status='old',err=333)
1285         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1286      &  write (iout,'(a)') 'intinname',intinname
1287         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1288         goto 334
1289   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1290 #ifdef MPI
1291         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1292 #endif
1293         stop 'Error opening angle file.'
1294   334   continue
1295
1296       endif
1297 C Generate distance constraints, if the PDB structure is to be regularized.
1298       if (nthread.gt.0) then
1299         call read_threadbase
1300       endif
1301       call setup_var
1302       if (me.eq.king .or. .not. out1file)
1303      & call intout
1304       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1305         write (iout,'(/a,i3,a)')
1306      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1307         write (iout,'(20i4)') (iss(i),i=1,ns)
1308        if (dyn_ss) then
1309           write(iout,*)"Running with dynamic disulfide-bond formation"
1310        else
1311         write (iout,'(/a/)') 'Pre-formed links are:'
1312         do i=1,nss
1313           i1=ihpb(i)-nres
1314           i2=jhpb(i)-nres
1315           it1=itype(i1)
1316           it2=itype(i2)
1317           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1318      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1319      &    ebr,forcon(i)
1320         enddo
1321         write (iout,'(a)')
1322        endif
1323       endif
1324       if (ns.gt.0.and.dyn_ss) then
1325           do i=nss+1,nhpb
1326             ihpb(i-nss)=ihpb(i)
1327             jhpb(i-nss)=jhpb(i)
1328             forcon(i-nss)=forcon(i)
1329             dhpb(i-nss)=dhpb(i)
1330           enddo
1331           nhpb=nhpb-nss
1332           nss=0
1333           call hpb_partition
1334           do i=1,ns
1335             dyn_ss_mask(iss(i))=.true.
1336           enddo
1337       endif
1338       if (i2ndstr.gt.0) call secstrp2dihc
1339 c      call geom_to_var(nvar,x)
1340 c      call etotal(energia(0))
1341 c      call enerprint(energia(0))
1342 c      call briefout(0,etot)
1343 c      stop
1344 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1345 cd    write (iout,'(a)') 'Variable list:'
1346 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1347 #ifdef MPI
1348       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1349      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1350      &  'Processor',myrank,': end reading molecular data.'
1351 #endif
1352       return
1353       end
1354 c--------------------------------------------------------------------------
1355       logical function seq_comp(itypea,itypeb,length)
1356       implicit none
1357       integer length,itypea(length),itypeb(length)
1358       integer i
1359       do i=1,length
1360         if (itypea(i).ne.itypeb(i)) then
1361           seq_comp=.false.
1362           return
1363         endif
1364       enddo
1365       seq_comp=.true.
1366       return
1367       end
1368 c-----------------------------------------------------------------------------
1369       subroutine read_bridge
1370 C Read information about disulfide bridges.
1371       implicit real*8 (a-h,o-z)
1372       include 'DIMENSIONS'
1373 #ifdef MPI
1374       include 'mpif.h'
1375 #endif
1376       include 'COMMON.IOUNITS'
1377       include 'COMMON.GEO'
1378       include 'COMMON.VAR'
1379       include 'COMMON.INTERACT'
1380       include 'COMMON.LOCAL'
1381       include 'COMMON.NAMES'
1382       include 'COMMON.CHAIN'
1383       include 'COMMON.FFIELD'
1384       include 'COMMON.SBRIDGE'
1385       include 'COMMON.HEADER'
1386       include 'COMMON.CONTROL'
1387       include 'COMMON.DBASE'
1388       include 'COMMON.THREAD'
1389       include 'COMMON.TIME1'
1390       include 'COMMON.SETUP'
1391 C Read bridging residues.
1392       read (inp,*) ns,(iss(i),i=1,ns)
1393       print *,'ns=',ns
1394       if(me.eq.king.or..not.out1file)
1395      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1396 C Check whether the specified bridging residues are cystines.
1397       do i=1,ns
1398         if (itype(iss(i)).ne.1) then
1399           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1400      &   'Do you REALLY think that the residue ',
1401      &    restyp(itype(iss(i))),i,
1402      &   ' can form a disulfide bridge?!!!'
1403           write (*,'(2a,i3,a)')
1404      &   'Do you REALLY think that the residue ',
1405      &    restyp(itype(iss(i))),i,
1406      &   ' can form a disulfide bridge?!!!'
1407 #ifdef MPI
1408          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1409          stop
1410 #endif
1411         endif
1412       enddo
1413 C Read preformed bridges.
1414       if (ns.gt.0) then
1415       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1416       if(fg_rank.eq.0)
1417      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1418       if (nss.gt.0) then
1419         nhpb=nss
1420 C Check if the residues involved in bridges are in the specified list of
1421 C bridging residues.
1422         do i=1,nss
1423           do j=1,i-1
1424             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1425      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1426               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1427      &      ' contains residues present in other pairs.'
1428               write (*,'(a,i3,a)') 'Disulfide pair',i,
1429      &      ' contains residues present in other pairs.'
1430 #ifdef MPI
1431               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1432               stop
1433 #endif
1434             endif
1435           enddo
1436           do j=1,ns
1437             if (ihpb(i).eq.iss(j)) goto 10
1438           enddo
1439           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1440    10     continue
1441           do j=1,ns
1442             if (jhpb(i).eq.iss(j)) goto 20
1443           enddo
1444           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1445    20     continue
1446           dhpb(i)=dbr
1447           forcon(i)=fbr
1448         enddo
1449         do i=1,nss
1450           ihpb(i)=ihpb(i)+nres
1451           jhpb(i)=jhpb(i)+nres
1452         enddo
1453       endif
1454       endif
1455       return
1456       end
1457 c----------------------------------------------------------------------------
1458       subroutine read_x(kanal,*)
1459       implicit real*8 (a-h,o-z)
1460       include 'DIMENSIONS'
1461       include 'COMMON.GEO'
1462       include 'COMMON.VAR'
1463       include 'COMMON.CHAIN'
1464       include 'COMMON.IOUNITS'
1465       include 'COMMON.CONTROL'
1466       include 'COMMON.LOCAL'
1467       include 'COMMON.INTERACT'
1468 c Read coordinates from input
1469 c
1470       read(kanal,'(8f10.5)',end=10,err=10)
1471      &  ((c(l,k),l=1,3),k=1,nres),
1472      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1473       do j=1,3
1474         c(j,nres+1)=c(j,1)
1475         c(j,2*nres)=c(j,nres)
1476       enddo
1477       call int_from_cart1(.false.)
1478       do i=1,nres-1
1479         do j=1,3
1480           dc(j,i)=c(j,i+1)-c(j,i)
1481           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1482         enddo
1483       enddo
1484       do i=nnt,nct
1485         if (itype(i).ne.10) then
1486           do j=1,3
1487             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1488             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1489           enddo
1490         endif
1491       enddo
1492
1493       return
1494    10 return1
1495       end
1496 c----------------------------------------------------------------------------
1497       subroutine read_threadbase
1498       implicit real*8 (a-h,o-z)
1499       include 'DIMENSIONS'
1500       include 'COMMON.IOUNITS'
1501       include 'COMMON.GEO'
1502       include 'COMMON.VAR'
1503       include 'COMMON.INTERACT'
1504       include 'COMMON.LOCAL'
1505       include 'COMMON.NAMES'
1506       include 'COMMON.CHAIN'
1507       include 'COMMON.FFIELD'
1508       include 'COMMON.SBRIDGE'
1509       include 'COMMON.HEADER'
1510       include 'COMMON.CONTROL'
1511       include 'COMMON.DBASE'
1512       include 'COMMON.THREAD'
1513       include 'COMMON.TIME1'
1514 C Read pattern database for threading.
1515       read (icbase,*) nseq
1516       do i=1,nseq
1517         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1518      &   nres_base(2,i),nres_base(3,i)
1519         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1520      &   nres_base(1,i))
1521 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1522 c    &   nres_base(2,i),nres_base(3,i)
1523 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1524 c    &   nres_base(1,i))
1525       enddo
1526       close (icbase)
1527       if (weidis.eq.0.0D0) weidis=0.1D0
1528       do i=nnt,nct
1529         do j=i+2,nct
1530           nhpb=nhpb+1
1531           ihpb(nhpb)=i
1532           jhpb(nhpb)=j
1533           forcon(nhpb)=weidis
1534         enddo
1535       enddo
1536       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1537       write (iout,'(a,i5)') 'nexcl: ',nexcl
1538       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1539       return
1540       end
1541 c------------------------------------------------------------------------------
1542       subroutine setup_var
1543       implicit real*8 (a-h,o-z)
1544       include 'DIMENSIONS'
1545       include 'COMMON.IOUNITS'
1546       include 'COMMON.GEO'
1547       include 'COMMON.VAR'
1548       include 'COMMON.INTERACT'
1549       include 'COMMON.LOCAL'
1550       include 'COMMON.NAMES'
1551       include 'COMMON.CHAIN'
1552       include 'COMMON.FFIELD'
1553       include 'COMMON.SBRIDGE'
1554       include 'COMMON.HEADER'
1555       include 'COMMON.CONTROL'
1556       include 'COMMON.DBASE'
1557       include 'COMMON.THREAD'
1558       include 'COMMON.TIME1'
1559 C Set up variable list.
1560       ntheta=nres-2
1561       nphi=nres-3
1562       nvar=ntheta+nphi
1563       nside=0
1564       do i=2,nres-1
1565         if (itype(i).ne.10) then
1566           nside=nside+1
1567           ialph(i,1)=nvar+nside
1568           ialph(nside,2)=i
1569         endif
1570       enddo
1571       if (indphi.gt.0) then
1572         nvar=nphi
1573       else if (indback.gt.0) then
1574         nvar=nphi+ntheta
1575       else
1576         nvar=nvar+2*nside
1577       endif
1578 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1579       return
1580       end
1581 c----------------------------------------------------------------------------
1582       subroutine gen_dist_constr
1583 C Generate CA distance constraints.
1584       implicit real*8 (a-h,o-z)
1585       include 'DIMENSIONS'
1586       include 'COMMON.IOUNITS'
1587       include 'COMMON.GEO'
1588       include 'COMMON.VAR'
1589       include 'COMMON.INTERACT'
1590       include 'COMMON.LOCAL'
1591       include 'COMMON.NAMES'
1592       include 'COMMON.CHAIN'
1593       include 'COMMON.FFIELD'
1594       include 'COMMON.SBRIDGE'
1595       include 'COMMON.HEADER'
1596       include 'COMMON.CONTROL'
1597       include 'COMMON.DBASE'
1598       include 'COMMON.THREAD'
1599       include 'COMMON.TIME1'
1600       dimension itype_pdb(maxres)
1601       common /pizda/ itype_pdb
1602       character*2 iden
1603 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1604 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1605 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1606 cd     & ' nsup',nsup
1607       do i=nstart_sup,nstart_sup+nsup-1
1608 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1609 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1610         do j=i+2,nstart_sup+nsup-1
1611           nhpb=nhpb+1
1612           ihpb(nhpb)=i+nstart_seq-nstart_sup
1613           jhpb(nhpb)=j+nstart_seq-nstart_sup
1614           forcon(nhpb)=weidis
1615           dhpb(nhpb)=dist(i,j)
1616         enddo
1617       enddo
1618 cd      write (iout,'(a)') 'Distance constraints:'
1619 cd      do i=nss+1,nhpb
1620 cd        ii=ihpb(i)
1621 cd        jj=jhpb(i)
1622 cd        iden='CA'
1623 cd        if (ii.gt.nres) then
1624 cd          iden='SC'
1625 cd          ii=ii-nres
1626 cd          jj=jj-nres
1627 cd        endif
1628 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1629 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1630 cd     &  dhpb(i),forcon(i)
1631 cd      enddo
1632       return
1633       end
1634 c----------------------------------------------------------------------------
1635       subroutine map_read
1636       implicit real*8 (a-h,o-z)
1637       include 'DIMENSIONS'
1638       include 'COMMON.MAP'
1639       include 'COMMON.IOUNITS'
1640       character*3 angid(4) /'THE','PHI','ALP','OME'/
1641       character*80 mapcard,ucase
1642       do imap=1,nmap
1643         read (inp,'(a)') mapcard
1644         mapcard=ucase(mapcard)
1645         if (index(mapcard,'PHI').gt.0) then
1646           kang(imap)=1
1647         else if (index(mapcard,'THE').gt.0) then
1648           kang(imap)=2
1649         else if (index(mapcard,'ALP').gt.0) then
1650           kang(imap)=3
1651         else if (index(mapcard,'OME').gt.0) then
1652           kang(imap)=4
1653         else
1654           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1655           stop 'Error - illegal variable spec in MAP card.'
1656         endif
1657         call readi (mapcard,'RES1',res1(imap),0)
1658         call readi (mapcard,'RES2',res2(imap),0)
1659         if (res1(imap).eq.0) then
1660           res1(imap)=res2(imap)
1661         else if (res2(imap).eq.0) then
1662           res2(imap)=res1(imap)
1663         endif
1664         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1665           write (iout,'(a)')
1666      &    'Error - illegal definition of variable group in MAP.'
1667           stop 'Error - illegal definition of variable group in MAP.'
1668         endif
1669         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1670         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1671         call readi(mapcard,'NSTEP',nstep(imap),0)
1672         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1673           write (iout,'(a)')
1674      &     'Illegal boundary and/or step size specification in MAP.'
1675           stop 'Illegal boundary and/or step size specification in MAP.'
1676         endif
1677       enddo ! imap
1678       return
1679       end
1680 c----------------------------------------------------------------------------
1681 csa      subroutine csaread
1682 csa      implicit real*8 (a-h,o-z)
1683 csa      include 'DIMENSIONS'
1684 csa      include 'COMMON.IOUNITS'
1685 csa      include 'COMMON.GEO'
1686 csa      include 'COMMON.CSA'
1687 csa      include 'COMMON.BANK'
1688 csa      include 'COMMON.CONTROL'
1689 csa      character*80 ucase
1690 csa      character*620 mcmcard
1691 csa      call card_concat(mcmcard)
1692 csa
1693 csa      call readi(mcmcard,'NCONF',nconf,50)
1694 csa      call readi(mcmcard,'NADD',nadd,0)
1695 csa      call readi(mcmcard,'JSTART',jstart,1)
1696 csa      call readi(mcmcard,'JEND',jend,1)
1697 csa      call readi(mcmcard,'NSTMAX',nstmax,500000)
1698 csa      call readi(mcmcard,'N0',n0,1)
1699 csa      call readi(mcmcard,'N1',n1,6)
1700 csa      call readi(mcmcard,'N2',n2,4)
1701 csa      call readi(mcmcard,'N3',n3,0)
1702 csa      call readi(mcmcard,'N4',n4,0)
1703 csa      call readi(mcmcard,'N5',n5,0)
1704 csa      call readi(mcmcard,'N6',n6,10)
1705 csa      call readi(mcmcard,'N7',n7,0)
1706 csa      call readi(mcmcard,'N8',n8,0)
1707 csa      call readi(mcmcard,'N9',n9,0)
1708 csa      call readi(mcmcard,'N14',n14,0)
1709 csa      call readi(mcmcard,'N15',n15,0)
1710 csa      call readi(mcmcard,'N16',n16,0)
1711 csa      call readi(mcmcard,'N17',n17,0)
1712 csa      call readi(mcmcard,'N18',n18,0)
1713 csa
1714 csa      vdisulf=(index(mcmcard,'DYNSS').gt.0)
1715 csa
1716 csa      call readi(mcmcard,'NDIFF',ndiff,2)
1717 csa      call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1718 csa      call readi(mcmcard,'IS1',is1,1)
1719 csa      call readi(mcmcard,'IS2',is2,8)
1720 csa      call readi(mcmcard,'NRAN0',nran0,4)
1721 csa      call readi(mcmcard,'NRAN1',nran1,2)
1722 csa      call readi(mcmcard,'IRR',irr,1)
1723 csa      call readi(mcmcard,'NSEED',nseed,20)
1724 csa      call readi(mcmcard,'NTOTAL',ntotal,10000)
1725 csa      call reada(mcmcard,'CUT1',cut1,2.0d0)
1726 csa      call reada(mcmcard,'CUT2',cut2,5.0d0)
1727 csa      call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1728 csa      call readi(mcmcard,'ICMAX',icmax,3)
1729 csa      call readi(mcmcard,'IRESTART',irestart,0)
1730 csac!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1731 csa      ntbankm=0
1732 csac!bankt
1733 csa      call reada(mcmcard,'DELE',dele,20.0d0)
1734 csa      call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1735 csa      call readi(mcmcard,'IREF',iref,0)
1736 csa      call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1737 csa      call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1738 csa      call readi(mcmcard,'NCONF_IN',nconf_in,0)
1739 csa      call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1740 csa      write (iout,*) "NCONF_IN",nconf_in
1741 csa      return
1742 csa      end
1743 c----------------------------------------------------------------------------
1744 cfmc      subroutine mcmfread
1745 cfmc      implicit real*8 (a-h,o-z)
1746 cfmc      include 'DIMENSIONS'
1747 cfmc      include 'COMMON.MCMF'
1748 cfmc      include 'COMMON.IOUNITS'
1749 cfmc      include 'COMMON.GEO'
1750 cfmc      character*80 ucase
1751 cfmc      character*620 mcmcard
1752 cfmc      call card_concat(mcmcard)
1753 cfmc
1754 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1755 cfmc      write(iout,*)'MAXRANT=',maxrant
1756 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1757 cfmc      write(iout,*)'MAXFAM=',maxfam
1758 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1759 cfmc      write(iout,*)'NNET1=',nnet1
1760 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1761 cfmc      write(iout,*)'NNET2=',nnet2
1762 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1763 cfmc      write(iout,*)'NNET3=',nnet3
1764 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1765 cfmc      write(iout,*)'ILASTT=',ilastt
1766 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1767 cfmc      write(iout,*)'MAXSTR=',maxstr
1768 cfmc      maxstr_f=maxstr/maxfam
1769 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1770 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1771 cfmc      write(iout,*)'NMCMF=',nmcmf
1772 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1773 cfmc      write(iout,*)'IFOCUS=',ifocus
1774 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1775 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1776 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1777 cfmc      write(iout,*)'INTPRT=',intprt
1778 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1779 cfmc      write(iout,*)'IPRT=',iprt
1780 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1781 cfmc      write(iout,*)'IMAXTR=',imaxtr
1782 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1783 cfmc      write(iout,*)'MAXEVEN=',maxeven
1784 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1785 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1786 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1787 cfmc      write(iout,*)'INIMIN=',inimin
1788 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1789 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1790 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1791 cfmc      write(iout,*)'NTHREAD=',nthread
1792 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1793 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1794 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1795 cfmc      write(iout,*)'MAXPERT=',maxpert
1796 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1797 cfmc      write(iout,*)'IRMSD=',irmsd
1798 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1799 cfmc      write(iout,*)'DENEMIN=',denemin
1800 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1801 cfmc      write(iout,*)'RCUT1S=',rcut1s
1802 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1803 cfmc      write(iout,*)'RCUT1E=',rcut1e
1804 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1805 cfmc      write(iout,*)'RCUT2S=',rcut2s
1806 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1807 cfmc      write(iout,*)'RCUT2E=',rcut2e
1808 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1809 cfmc      write(iout,*)'DPERT1=',d_pert1
1810 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1811 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1812 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1813 cfmc      write(iout,*)'DPERT2=',d_pert2
1814 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1815 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1816 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1817 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1818 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1819 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1820 cfmc      d_pert1=deg2rad*d_pert1
1821 cfmc      d_pert1a=deg2rad*d_pert1a
1822 cfmc      d_pert2=deg2rad*d_pert2
1823 cfmc      d_pert2a=deg2rad*d_pert2a
1824 cfmc      d_pert2b=deg2rad*d_pert2b
1825 cfmc      d_pert2c=deg2rad*d_pert2c
1826 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1827 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1828 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1829 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1830 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1831 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1832 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1833 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1834 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1835 cfmc      write(iout,*)'RCUTINI=',rcutini
1836 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1837 cfmc      write(iout,*)'GRAT=',grat
1838 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1839 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1840 cfmc
1841 cfmc      return
1842 cfmc      end
1843 c----------------------------------------------------------------------------
1844       subroutine mcmread
1845       implicit real*8 (a-h,o-z)
1846       include 'DIMENSIONS'
1847       include 'COMMON.MCM'
1848       include 'COMMON.MCE'
1849       include 'COMMON.IOUNITS'
1850       character*80 ucase
1851       character*320 mcmcard
1852       call card_concat(mcmcard)
1853       call readi(mcmcard,'MAXACC',maxacc,100)
1854       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1855       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1856       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1857       call readi(mcmcard,'MAXREPM',maxrepm,200)
1858       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1859       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1860       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1861       call reada(mcmcard,'E_UP',e_up,5.0D0)
1862       call reada(mcmcard,'DELTE',delte,0.1D0)
1863       call readi(mcmcard,'NSWEEP',nsweep,5)
1864       call readi(mcmcard,'NSTEPH',nsteph,0)
1865       call readi(mcmcard,'NSTEPC',nstepc,0)
1866       call reada(mcmcard,'TMIN',tmin,298.0D0)
1867       call reada(mcmcard,'TMAX',tmax,298.0D0)
1868       call readi(mcmcard,'NWINDOW',nwindow,0)
1869       call readi(mcmcard,'PRINT_MC',print_mc,0)
1870       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1871       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1872       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1873       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1874       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1875       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1876       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1877       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1878       if (nwindow.gt.0) then
1879         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1880         do i=1,nwindow
1881           winlen(i)=winend(i)-winstart(i)+1
1882         enddo
1883       endif
1884       if (tmax.lt.tmin) tmax=tmin
1885       if (tmax.eq.tmin) then
1886         nstepc=0
1887         nsteph=0
1888       endif
1889       if (nstepc.gt.0 .and. nsteph.gt.0) then
1890         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1891         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1892       endif
1893 C Probabilities of different move types
1894       sumpro_type(0)=0.0D0
1895       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1896       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1897       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1898       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1899       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1900       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1901       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1902       do i=1,MaxMoveType
1903         print *,'i',i,' sumprotype',sumpro_type(i)
1904         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1905         print *,'i',i,' sumprotype',sumpro_type(i)
1906       enddo
1907       return
1908       end
1909 c----------------------------------------------------------------------------
1910       subroutine read_minim
1911       implicit real*8 (a-h,o-z)
1912       include 'DIMENSIONS'
1913       include 'COMMON.MINIM'
1914       include 'COMMON.IOUNITS'
1915       character*80 ucase
1916       character*320 minimcard
1917       call card_concat(minimcard)
1918       call readi(minimcard,'MAXMIN',maxmin,2000)
1919       call readi(minimcard,'MAXFUN',maxfun,5000)
1920       call readi(minimcard,'MINMIN',minmin,maxmin)
1921       call readi(minimcard,'MINFUN',minfun,maxmin)
1922       call reada(minimcard,'TOLF',tolf,1.0D-2)
1923       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1924       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1925       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1926       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1927       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1928      &         'Options in energy minimization:'
1929       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1930      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1931      & 'MinMin:',MinMin,' MinFun:',MinFun,
1932      & ' TolF:',TolF,' RTolF:',RTolF
1933       return
1934       end
1935 c----------------------------------------------------------------------------
1936       subroutine read_angles(kanal,*)
1937       implicit real*8 (a-h,o-z)
1938       include 'DIMENSIONS'
1939       include 'COMMON.GEO'
1940       include 'COMMON.VAR'
1941       include 'COMMON.CHAIN'
1942       include 'COMMON.IOUNITS'
1943       include 'COMMON.CONTROL'
1944 c Read angles from input
1945 c
1946        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1947        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1948        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1949        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1950
1951        do i=1,nres
1952 c 9/7/01 avoid 180 deg valence angle
1953         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1954 c
1955         theta(i)=deg2rad*theta(i)
1956         phi(i)=deg2rad*phi(i)
1957         alph(i)=deg2rad*alph(i)
1958         omeg(i)=deg2rad*omeg(i)
1959        enddo
1960       return
1961    10 return1
1962       end
1963 c----------------------------------------------------------------------------
1964       subroutine reada(rekord,lancuch,wartosc,default)
1965       implicit none
1966       character*(*) rekord,lancuch
1967       double precision wartosc,default
1968       integer ilen,iread
1969       external ilen
1970       iread=index(rekord,lancuch)
1971       if (iread.eq.0) then
1972         wartosc=default
1973         return
1974       endif
1975       iread=iread+ilen(lancuch)+1
1976       read (rekord(iread:),*,err=10,end=10) wartosc
1977       return
1978   10  wartosc=default
1979       return
1980       end
1981 c----------------------------------------------------------------------------
1982       subroutine readi(rekord,lancuch,wartosc,default)
1983       implicit none
1984       character*(*) rekord,lancuch
1985       integer wartosc,default
1986       integer ilen,iread
1987       external ilen
1988       iread=index(rekord,lancuch)
1989       if (iread.eq.0) then
1990         wartosc=default
1991         return
1992       endif
1993       iread=iread+ilen(lancuch)+1
1994       read (rekord(iread:),*,err=10,end=10) wartosc
1995       return
1996   10  wartosc=default
1997       return
1998       end
1999 c----------------------------------------------------------------------------
2000       subroutine multreadi(rekord,lancuch,tablica,dim,default)
2001       implicit none
2002       integer dim,i
2003       integer tablica(dim),default
2004       character*(*) rekord,lancuch
2005       character*80 aux
2006       integer ilen,iread
2007       external ilen
2008       do i=1,dim
2009         tablica(i)=default
2010       enddo
2011       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2012       if (iread.eq.0) return
2013       iread=iread+ilen(lancuch)+1
2014       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2015    10 return
2016       end
2017 c----------------------------------------------------------------------------
2018       subroutine multreada(rekord,lancuch,tablica,dim,default)
2019       implicit none
2020       integer dim,i
2021       double precision tablica(dim),default
2022       character*(*) rekord,lancuch
2023       character*80 aux
2024       integer ilen,iread
2025       external ilen
2026       do i=1,dim
2027         tablica(i)=default
2028       enddo
2029       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2030       if (iread.eq.0) return
2031       iread=iread+ilen(lancuch)+1
2032       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2033    10 return
2034       end
2035 c----------------------------------------------------------------------------
2036       subroutine openunits
2037       implicit real*8 (a-h,o-z)
2038       include 'DIMENSIONS'
2039 #ifdef MPI
2040       include 'mpif.h'
2041       character*16 form,nodename
2042       integer nodelen
2043 #endif
2044       include 'COMMON.SETUP'
2045       include 'COMMON.IOUNITS'
2046       include 'COMMON.MD'
2047       include 'COMMON.CONTROL'
2048       integer lenpre,lenpot,ilen,lentmp
2049       external ilen
2050       character*3 out1file_text,ucase
2051       character*3 ll
2052       external ucase
2053 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2054       call getenv_loc("PREFIX",prefix)
2055       pref_orig = prefix
2056       call getenv_loc("POT",pot)
2057       call getenv_loc("DIRTMP",tmpdir)
2058       call getenv_loc("CURDIR",curdir)
2059       call getenv_loc("OUT1FILE",out1file_text)
2060 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2061       out1file_text=ucase(out1file_text)
2062       if (out1file_text(1:1).eq."Y") then
2063         out1file=.true.
2064       else
2065         out1file=fg_rank.gt.0
2066       endif
2067       lenpre=ilen(prefix)
2068       lenpot=ilen(pot)
2069       lentmp=ilen(tmpdir)
2070       if (lentmp.gt.0) then
2071           write (*,'(80(1h!))')
2072           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
2073           write (*,'(80(1h!))')
2074           write (*,*)"All output files will be on node /tmp directory."
2075 #ifdef MPI
2076         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2077         if (me.eq.king) then
2078           write (*,*) "The master node is ",nodename
2079         else if (fg_rank.eq.0) then
2080           write (*,*) "I am the CG slave node ",nodename
2081         else
2082           write (*,*) "I am the FG slave node ",nodename
2083         endif
2084 #endif
2085         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2086         lenpre = lentmp+lenpre+1
2087       endif
2088       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2089 C Get the names and open the input files
2090 #if defined(WINIFL) || defined(WINPGI)
2091       open(1,file=pref_orig(:ilen(pref_orig))//
2092      &  '.inp',status='old',readonly,shared)
2093        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2094 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2095 C Get parameter filenames and open the parameter files.
2096       call getenv_loc('BONDPAR',bondname)
2097       open (ibond,file=bondname,status='old',readonly,shared)
2098       call getenv_loc('THETPAR',thetname)
2099       open (ithep,file=thetname,status='old',readonly,shared)
2100 #ifndef CRYST_THETA
2101       call getenv_loc('THETPARPDB',thetname_pdb)
2102       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2103 #endif
2104       call getenv_loc('ROTPAR',rotname)
2105       open (irotam,file=rotname,status='old',readonly,shared)
2106 #ifndef CRYST_SC
2107       call getenv_loc('ROTPARPDB',rotname_pdb)
2108       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2109 #endif
2110       call getenv_loc('TORPAR',torname)
2111       open (itorp,file=torname,status='old',readonly,shared)
2112       call getenv_loc('TORDPAR',tordname)
2113       open (itordp,file=tordname,status='old',readonly,shared)
2114       call getenv_loc('FOURIER',fouriername)
2115       open (ifourier,file=fouriername,status='old',readonly,shared)
2116       call getenv_loc('ELEPAR',elename)
2117       open (ielep,file=elename,status='old',readonly,shared)
2118       call getenv_loc('SIDEPAR',sidename)
2119       open (isidep,file=sidename,status='old',readonly,shared)
2120 #elif (defined CRAY) || (defined AIX)
2121       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2122      &  action='read')
2123 c      print *,"Processor",myrank," opened file 1"
2124       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2125 c      print *,"Processor",myrank," opened file 9"
2126 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2127 C Get parameter filenames and open the parameter files.
2128       call getenv_loc('BONDPAR',bondname)
2129       open (ibond,file=bondname,status='old',action='read')
2130 c      print *,"Processor",myrank," opened file IBOND"
2131       call getenv_loc('THETPAR',thetname)
2132       open (ithep,file=thetname,status='old',action='read')
2133 c      print *,"Processor",myrank," opened file ITHEP"
2134 #ifndef CRYST_THETA
2135       call getenv_loc('THETPARPDB',thetname_pdb)
2136       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2137 #endif
2138       call getenv_loc('ROTPAR',rotname)
2139       open (irotam,file=rotname,status='old',action='read')
2140 c      print *,"Processor",myrank," opened file IROTAM"
2141 #ifndef CRYST_SC
2142       call getenv_loc('ROTPARPDB',rotname_pdb)
2143       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2144 #endif
2145       call getenv_loc('TORPAR',torname)
2146       open (itorp,file=torname,status='old',action='read')
2147 c      print *,"Processor",myrank," opened file ITORP"
2148       call getenv_loc('TORDPAR',tordname)
2149       open (itordp,file=tordname,status='old',action='read')
2150 c      print *,"Processor",myrank," opened file ITORDP"
2151       call getenv_loc('SCCORPAR',sccorname)
2152       open (isccor,file=sccorname,status='old',action='read')
2153 c      print *,"Processor",myrank," opened file ISCCOR"
2154       call getenv_loc('FOURIER',fouriername)
2155       open (ifourier,file=fouriername,status='old',action='read')
2156 c      print *,"Processor",myrank," opened file IFOURIER"
2157       call getenv_loc('ELEPAR',elename)
2158       open (ielep,file=elename,status='old',action='read')
2159 c      print *,"Processor",myrank," opened file IELEP"
2160       call getenv_loc('SIDEPAR',sidename)
2161       open (isidep,file=sidename,status='old',action='read')
2162 c      print *,"Processor",myrank," opened file ISIDEP"
2163 c      print *,"Processor",myrank," opened parameter files"
2164 #elif (defined G77)
2165       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2166       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2167 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2168 C Get parameter filenames and open the parameter files.
2169       call getenv_loc('BONDPAR',bondname)
2170       open (ibond,file=bondname,status='old')
2171       call getenv_loc('THETPAR',thetname)
2172       open (ithep,file=thetname,status='old')
2173 #ifndef CRYST_THETA
2174       call getenv_loc('THETPARPDB',thetname_pdb)
2175       open (ithep_pdb,file=thetname_pdb,status='old')
2176 #endif
2177       call getenv_loc('ROTPAR',rotname)
2178       open (irotam,file=rotname,status='old')
2179 #ifndef CRYST_SC
2180       call getenv_loc('ROTPARPDB',rotname_pdb)
2181       open (irotam_pdb,file=rotname_pdb,status='old')
2182 #endif
2183       call getenv_loc('TORPAR',torname)
2184       open (itorp,file=torname,status='old')
2185       call getenv_loc('TORDPAR',tordname)
2186       open (itordp,file=tordname,status='old')
2187       call getenv_loc('SCCORPAR',sccorname)
2188       open (isccor,file=sccorname,status='old')
2189       call getenv_loc('FOURIER',fouriername)
2190       open (ifourier,file=fouriername,status='old')
2191       call getenv_loc('ELEPAR',elename)
2192       open (ielep,file=elename,status='old')
2193       call getenv_loc('SIDEPAR',sidename)
2194       open (isidep,file=sidename,status='old')
2195 #else
2196       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2197      &action='read')
2198        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2199 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2200 C Get parameter filenames and open the parameter files.
2201       call getenv_loc('BONDPAR',bondname)
2202       open (ibond,file=bondname,status='old',action='read')
2203       call getenv_loc('THETPAR',thetname)
2204       open (ithep,file=thetname,status='old',action='read')
2205 #ifndef CRYST_THETA
2206       call getenv_loc('THETPARPDB',thetname_pdb)
2207       print *,"thetname_pdb ",thetname_pdb
2208       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2209       print *,ithep_pdb," opened"
2210 #endif
2211       call getenv_loc('ROTPAR',rotname)
2212       open (irotam,file=rotname,status='old',action='read')
2213 #ifndef CRYST_SC
2214       call getenv_loc('ROTPARPDB',rotname_pdb)
2215       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2216 #endif
2217       call getenv_loc('TORPAR',torname)
2218       open (itorp,file=torname,status='old',action='read')
2219       call getenv_loc('TORDPAR',tordname)
2220       open (itordp,file=tordname,status='old',action='read')
2221       call getenv_loc('SCCORPAR',sccorname)
2222       open (isccor,file=sccorname,status='old',action='read')
2223       call getenv_loc('FOURIER',fouriername)
2224       open (ifourier,file=fouriername,status='old',action='read')
2225       call getenv_loc('ELEPAR',elename)
2226       open (ielep,file=elename,status='old',action='read')
2227       call getenv_loc('SIDEPAR',sidename)
2228       open (isidep,file=sidename,status='old',action='read')
2229 #endif
2230 #ifndef OLDSCP
2231 C
2232 C 8/9/01 In the newest version SCp interaction constants are read from a file
2233 C Use -DOLDSCP to use hard-coded constants instead.
2234 C
2235       call getenv_loc('SCPPAR',scpname)
2236 #if defined(WINIFL) || defined(WINPGI)
2237       open (iscpp,file=scpname,status='old',readonly,shared)
2238 #elif (defined CRAY)  || (defined AIX)
2239       open (iscpp,file=scpname,status='old',action='read')
2240 #elif (defined G77)
2241       open (iscpp,file=scpname,status='old')
2242 #else
2243       open (iscpp,file=scpname,status='old',action='read')
2244 #endif
2245 #endif
2246       call getenv_loc('PATTERN',patname)
2247 #if defined(WINIFL) || defined(WINPGI)
2248       open (icbase,file=patname,status='old',readonly,shared)
2249 #elif (defined CRAY)  || (defined AIX)
2250       open (icbase,file=patname,status='old',action='read')
2251 #elif (defined G77)
2252       open (icbase,file=patname,status='old')
2253 #else
2254       open (icbase,file=patname,status='old',action='read')
2255 #endif
2256 #ifdef MPI
2257 C Open output file only for CG processes
2258 c      print *,"Processor",myrank," fg_rank",fg_rank
2259       if (fg_rank.eq.0) then
2260
2261       if (nodes.eq.1) then
2262         npos=3
2263       else
2264         npos = dlog10(dfloat(nodes-1))+1
2265       endif
2266       if (npos.lt.3) npos=3
2267       write (liczba,'(i1)') npos
2268       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2269      &  //')'
2270       write (liczba,form) me
2271       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2272      &  liczba(:ilen(liczba))
2273       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2274      &  //'.int'
2275       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2276      &  //'.pdb'
2277       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2278      &  liczba(:ilen(liczba))//'.mol2'
2279       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2280      &  liczba(:ilen(liczba))//'.stat'
2281       if (lentmp.gt.0)
2282      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2283      &      //liczba(:ilen(liczba))//'.stat')
2284       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2285      &  //'.rst'
2286       if(usampl) then
2287           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2288      & liczba(:ilen(liczba))//'.const'
2289       endif
2290
2291       endif
2292 #else
2293       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2294       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2295       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2296       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2297       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2298       if (lentmp.gt.0)
2299      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2300      &    //'.stat')
2301       rest2name=prefix(:ilen(prefix))//'.rst'
2302       if(usampl) then
2303          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2304       endif
2305 #endif
2306 #if defined(AIX) || defined(PGI)
2307       if (me.eq.king .or. .not. out1file)
2308      &   open(iout,file=outname,status='unknown')
2309 c#define DEBUG
2310 #ifdef DEBUG
2311       if (fg_rank.gt.0) then
2312         write (liczba,'(i3.3)') myrank/nfgtasks
2313         write (ll,'(bz,i3.3)') fg_rank
2314         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2315      &   status='unknown')
2316       endif
2317 #endif
2318 c#undef DEBUG
2319       if(me.eq.king) then
2320        open(igeom,file=intname,status='unknown',position='append')
2321        open(ipdb,file=pdbname,status='unknown')
2322        open(imol2,file=mol2name,status='unknown')
2323        open(istat,file=statname,status='unknown',position='append')
2324       else
2325 c1out       open(iout,file=outname,status='unknown')
2326       endif
2327 #else
2328       if (me.eq.king .or. .not.out1file)
2329      &    open(iout,file=outname,status='unknown')
2330 c#define DEBUG
2331 #ifdef DEBUG
2332       if (fg_rank.gt.0) then
2333         print "Processor",fg_rank," opening output file"
2334         write (liczba,'(i3.3)') myrank/nfgtasks
2335         write (ll,'(bz,i3.3)') fg_rank
2336         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2337      &   status='unknown')
2338       endif
2339 #endif
2340 c#undef DEBUG
2341       if(me.eq.king) then
2342        open(igeom,file=intname,status='unknown',access='append')
2343        open(ipdb,file=pdbname,status='unknown')
2344        open(imol2,file=mol2name,status='unknown')
2345        open(istat,file=statname,status='unknown',access='append')
2346       else
2347 c1out       open(iout,file=outname,status='unknown')
2348       endif
2349 #endif
2350 csa      csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2351 csa      csa_seed=prefix(:lenpre)//'.CSA.seed'
2352 csa      csa_history=prefix(:lenpre)//'.CSA.history'
2353 csa      csa_bank=prefix(:lenpre)//'.CSA.bank'
2354 csa      csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2355 csa      csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2356 csa      csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2357 csac!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2358 csa      csa_int=prefix(:lenpre)//'.int'
2359 csa      csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2360 csa      csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2361 csa      csa_in=prefix(:lenpre)//'.CSA.in'
2362 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2363 C Write file names
2364       if (me.eq.king)then
2365       write (iout,'(80(1h-))')
2366       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2367       write (iout,'(80(1h-))')
2368       write (iout,*) "Input file                      : ",
2369      &  pref_orig(:ilen(pref_orig))//'.inp'
2370       write (iout,*) "Output file                     : ",
2371      &  outname(:ilen(outname))
2372       write (iout,*)
2373       write (iout,*) "Sidechain potential file        : ",
2374      &  sidename(:ilen(sidename))
2375 #ifndef OLDSCP
2376       write (iout,*) "SCp potential file              : ",
2377      &  scpname(:ilen(scpname))
2378 #endif
2379       write (iout,*) "Electrostatic potential file    : ",
2380      &  elename(:ilen(elename))
2381       write (iout,*) "Cumulant coefficient file       : ",
2382      &  fouriername(:ilen(fouriername))
2383       write (iout,*) "Torsional parameter file        : ",
2384      &  torname(:ilen(torname))
2385       write (iout,*) "Double torsional parameter file : ",
2386      &  tordname(:ilen(tordname))
2387       write (iout,*) "SCCOR parameter file : ",
2388      &  sccorname(:ilen(sccorname))
2389       write (iout,*) "Bond & inertia constant file    : ",
2390      &  bondname(:ilen(bondname))
2391       write (iout,*) "Bending parameter file          : ",
2392      &  thetname(:ilen(thetname))
2393       write (iout,*) "Rotamer parameter file          : ",
2394      &  rotname(:ilen(rotname))
2395       write (iout,*) "Threading database              : ",
2396      &  patname(:ilen(patname))
2397       if (lentmp.ne.0)
2398      &write (iout,*)" DIRTMP                          : ",
2399      &  tmpdir(:lentmp)
2400       write (iout,'(80(1h-))')
2401       endif
2402       return
2403       end
2404 c----------------------------------------------------------------------------
2405       subroutine card_concat(card)
2406       implicit real*8 (a-h,o-z)
2407       include 'DIMENSIONS'
2408       include 'COMMON.IOUNITS'
2409       character*(*) card
2410       character*80 karta,ucase
2411       external ilen
2412       read (inp,'(a)') karta
2413       karta=ucase(karta)
2414       card=' '
2415       do while (karta(80:80).eq.'&')
2416         card=card(:ilen(card)+1)//karta(:79)
2417         read (inp,'(a)') karta
2418         karta=ucase(karta)
2419       enddo
2420       card=card(:ilen(card)+1)//karta
2421       return
2422       end
2423 c----------------------------------------------------------------------------------
2424       subroutine readrst
2425       implicit real*8 (a-h,o-z)
2426       include 'DIMENSIONS'
2427       include 'COMMON.CHAIN'
2428       include 'COMMON.IOUNITS'
2429       include 'COMMON.MD'
2430       include 'COMMON.CONTROL'
2431       open(irest2,file=rest2name,status='unknown')
2432       read(irest2,*) totT,EK,potE,totE,t_bath
2433       do i=1,2*nres
2434          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2435       enddo
2436       do i=1,2*nres
2437          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2438       enddo
2439       if(usampl.or.homol_nset.gt.1) then
2440              read (irest2,*) iset
2441       endif
2442       close(irest2)
2443       return
2444       end
2445 c---------------------------------------------------------------------------------
2446       subroutine read_fragments
2447       implicit real*8 (a-h,o-z)
2448       include 'DIMENSIONS'
2449 #ifdef MPI
2450       include 'mpif.h'
2451 #endif
2452       include 'COMMON.SETUP'
2453       include 'COMMON.CHAIN'
2454       include 'COMMON.IOUNITS'
2455       include 'COMMON.MD'
2456       include 'COMMON.CONTROL'
2457       read(inp,*) nset,nfrag,npair,nfrag_back
2458       if(me.eq.king.or..not.out1file)
2459      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2460      &  " nfrag_back",nfrag_back
2461       do iset=1,nset
2462          read(inp,*) mset(iset)
2463        do i=1,nfrag
2464          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2465      &     qinfrag(i,iset)
2466          if(me.eq.king.or..not.out1file)
2467      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2468      &     ifrag(2,i,iset), qinfrag(i,iset)
2469        enddo
2470        do i=1,npair
2471         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2472      &    qinpair(i,iset)
2473         if(me.eq.king.or..not.out1file)
2474      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2475      &    ipair(2,i,iset), qinpair(i,iset)
2476        enddo
2477        do i=1,nfrag_back
2478         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2479      &     wfrag_back(3,i,iset),
2480      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2481         if(me.eq.king.or..not.out1file)
2482      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2483      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2484        enddo
2485       enddo
2486       return
2487       end
2488 c-------------------------------------------------------------------------------
2489       subroutine read_dist_constr
2490       implicit real*8 (a-h,o-z)
2491       include 'DIMENSIONS'
2492 #ifdef MPI
2493       include 'mpif.h'
2494 #endif
2495       include 'COMMON.SETUP'
2496       include 'COMMON.CONTROL'
2497       include 'COMMON.CHAIN'
2498       include 'COMMON.IOUNITS'
2499       include 'COMMON.SBRIDGE'
2500       integer ifrag_(2,100),ipair_(2,100)
2501       double precision wfrag_(100),wpair_(100)
2502       character*500 controlcard
2503 c      write (iout,*) "Calling read_dist_constr"
2504 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2505 c      call flush(iout)
2506       call card_concat(controlcard)
2507       call readi(controlcard,"NFRAG",nfrag_,0)
2508       call readi(controlcard,"NPAIR",npair_,0)
2509       call readi(controlcard,"NDIST",ndist_,0)
2510       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2511       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2512       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2513       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2514       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2515 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2516 c      write (iout,*) "IFRAG"
2517 c      do i=1,nfrag_
2518 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2519 c      enddo
2520 c      write (iout,*) "IPAIR"
2521 c      do i=1,npair_
2522 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2523 c      enddo
2524       if (.not.refstr .and. nfrag.gt.0) then
2525         write (iout,*)
2526      &  "ERROR: no reference structure to compute distance restraints"
2527         write (iout,*)
2528      &  "Restraints must be specified explicitly (NDIST=number)"
2529         stop
2530       endif
2531       if (nfrag.lt.2 .and. npair.gt.0) then
2532         write (iout,*) "ERROR: Less than 2 fragments specified",
2533      &   " but distance restraints between pairs requested"
2534         stop
2535       endif
2536       call flush(iout)
2537       do i=1,nfrag_
2538         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2539         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2540      &    ifrag_(2,i)=nstart_sup+nsup-1
2541 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2542         call flush(iout)
2543         if (wfrag_(i).gt.0.0d0) then
2544         do j=ifrag_(1,i),ifrag_(2,i)-1
2545           do k=j+1,ifrag_(2,i)
2546 c            write (iout,*) "j",j," k",k
2547             ddjk=dist(j,k)
2548             if (constr_dist.eq.1) then
2549             nhpb=nhpb+1
2550             ihpb(nhpb)=j
2551             jhpb(nhpb)=k
2552               dhpb(nhpb)=ddjk
2553             forcon(nhpb)=wfrag_(i)
2554             else if (constr_dist.eq.2) then
2555               if (ddjk.le.dist_cut) then
2556                 nhpb=nhpb+1
2557                 ihpb(nhpb)=j
2558                 jhpb(nhpb)=k
2559                 dhpb(nhpb)=ddjk
2560                 forcon(nhpb)=wfrag_(i)
2561               endif
2562             else
2563               nhpb=nhpb+1
2564               ihpb(nhpb)=j
2565               jhpb(nhpb)=k
2566               dhpb(nhpb)=ddjk
2567               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2568             endif
2569 #ifdef MPI
2570             if (.not.out1file .or. me.eq.king)
2571      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2572      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2573 #else
2574             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2575      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2576 #endif
2577           enddo
2578         enddo
2579         endif
2580       enddo
2581       do i=1,npair_
2582         if (wpair_(i).gt.0.0d0) then
2583         ii = ipair_(1,i)
2584         jj = ipair_(2,i)
2585         if (ii.gt.jj) then
2586           itemp=ii
2587           ii=jj
2588           jj=itemp
2589         endif
2590         do j=ifrag_(1,ii),ifrag_(2,ii)
2591           do k=ifrag_(1,jj),ifrag_(2,jj)
2592             nhpb=nhpb+1
2593             ihpb(nhpb)=j
2594             jhpb(nhpb)=k
2595             forcon(nhpb)=wpair_(i)
2596             dhpb(nhpb)=dist(j,k)
2597 #ifdef MPI
2598             if (.not.out1file .or. me.eq.king)
2599      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2600      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2601 #else
2602             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2603      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2604 #endif
2605           enddo
2606         enddo
2607         endif
2608       enddo
2609       do i=1,ndist_
2610         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2611      &     ibecarb(i),forcon(nhpb+1)
2612         if (forcon(nhpb+1).gt.0.0d0) then
2613           nhpb=nhpb+1
2614           if (ibecarb(i).gt.0) then
2615             ihpb(i)=ihpb(i)+nres
2616             jhpb(i)=jhpb(i)+nres
2617           endif
2618           if (dhpb(nhpb).eq.0.0d0)
2619      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2620         endif
2621       enddo
2622 #ifdef MPI
2623       if (.not.out1file .or. me.eq.king) then
2624 #endif
2625       do i=1,nhpb
2626           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2627      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2628       enddo
2629       call flush(iout)
2630 #ifdef MPI
2631       endif
2632 #endif
2633       return
2634       end
2635 c-------------------------------------------------------------------------------
2636
2637       subroutine read_constr_homology
2638
2639       include 'DIMENSIONS'
2640 #ifdef MPI
2641       include 'mpif.h'
2642 #endif
2643       include 'COMMON.SETUP'
2644       include 'COMMON.CONTROL'
2645       include 'COMMON.CHAIN'
2646       include 'COMMON.IOUNITS'
2647       include 'COMMON.MD'
2648       include 'COMMON.GEO'
2649       include 'COMMON.INTERACT'
2650 c
2651 c For new homol impl
2652 c
2653       include 'COMMON.VAR'
2654 c
2655
2656 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2657 c    &                 dist_cut
2658 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2659 c    &    sigma_odl_temp(maxres,maxres,max_template)
2660       character*2 kic2
2661       character*24 model_ki_dist, model_ki_angle
2662       character*500 controlcard
2663       integer ki, i, j, k, l
2664       logical lprn /.true./
2665 c
2666 c     FP - Nov. 2014 Temporary specifications for new vars
2667 c
2668       double precision rescore_tmp,x12,y12,z12
2669       double precision, dimension (max_template,maxres) :: rescore
2670       character*24 pdbfile,tpl_k_rescore
2671 c -----------------------------------------------------------------
2672 c Reading multiple PDB ref structures and calculation of retraints
2673 c not using pre-computed ones stored in files model_ki_{dist,angle}
2674 c FP (Nov., 2014)
2675 c -----------------------------------------------------------------
2676 c
2677 c
2678 c Alternative: reading from input
2679       call card_concat(controlcard)
2680       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2681       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2682       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2683       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2684       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2685
2686       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2687       if (homol_nset.gt.1)then
2688          call card_concat(controlcard)
2689          read(controlcard,*) (waga_dist1(i),i=1,homol_nset)
2690          call card_concat(controlcard)
2691          read(controlcard,*) (waga_angle1(i),i=1,homol_nset)
2692          write(iout,*) "iset distance_weight angle_weight"
2693          do i=1,homol_nset
2694            write(iout,*) i,waga_dist1(i),waga_angle1(i)
2695          enddo
2696          iset=mod(kolor,homol_nset)+1
2697       else
2698        iset=1
2699 c       call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0)
2700 c       call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0)
2701       endif
2702
2703       write (iout,*) "nnt",nnt," nct",nct
2704       call flush(iout)
2705
2706
2707       lim_odl=0
2708       lim_dih=0
2709 c
2710 c  New
2711 c
2712       lim_theta=0
2713       lim_xx=0
2714 c
2715 c  Reading HM global scores (prob not required)
2716 c
2717 c     open (4,file="HMscore")
2718 c     do k=1,constr_homology
2719 c       read (4,*,end=521) hmscore_tmp
2720 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
2721 c       write(*,*) "Model", k, ":", hmscore(k)
2722 c     enddo
2723 c521  continue
2724
2725 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2726       do k=1,constr_homology
2727
2728         read(inp,'(a)') pdbfile
2729 c  Next stament causes error upon compilation (?)
2730 c       if(me.eq.king.or. .not. out1file)
2731 c         write (iout,'(2a)') 'PDB data will be read from file ',
2732 c    &   pdbfile(:ilen(pdbfile))
2733         open(ipdbin,file=pdbfile,status='old',err=33)
2734         goto 34
2735   33    write (iout,'(a)') 'Error opening PDB file.'
2736         stop
2737   34    continue
2738 c        print *,'Begin reading pdb data'
2739 c
2740 c Files containing res sim or local scores (former containing sigmas)
2741 c
2742
2743         write(kic2,'(bz,i2.2)') k
2744
2745         tpl_k_rescore="template"//kic2//".sco"
2746 c       tpl_k_sigma_odl="template"//kic2//".sigma_odl"
2747 c       tpl_k_sigma_dih="template"//kic2//".sigma_dih"
2748 c       tpl_k_sigma_theta="template"//kic2//".sigma_theta"
2749 c       tpl_k_sigma_d="template"//kic2//".sigma_d"
2750
2751         unres_pdb=.false.
2752         call readpdb
2753 c
2754 c     Distance restraints
2755 c
2756 c          ... --> odl(k,ii)
2757 C Copy the coordinates from reference coordinates (?)
2758         do i=1,2*nres
2759           do j=1,3
2760             c(j,i)=cref(j,i)
2761 c           write (iout,*) "c(",j,i,") =",c(j,i)
2762           enddo
2763         enddo
2764 c
2765 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2766 c
2767 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2768           open (ientin,file=tpl_k_rescore,status='old')
2769           do irec=1,maxdim ! loop for reading res sim
2770             if (irec.eq.1) then
2771                rescore(k,irec)=0.0d0
2772                goto 1301
2773             endif
2774             read (ientin,*,end=1401) rescore_tmp
2775 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2776             rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2777 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2778  1301     continue
2779           enddo
2780  1401   continue
2781           close (ientin)
2782 c         open (ientin,file=tpl_k_sigma_odl,status='old')
2783 c         do irec=1,maxdim ! loop for reading sigma_odl
2784 c            read (ientin,*,end=1401) i, j,
2785 c    &                                sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
2786 c            sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
2787 c    &       sigma_odl_temp(i+nnt-1,j+nnt-1,k)
2788 c         enddo
2789 c 1401   continue
2790 c         close (ientin)
2791         if (waga_dist.gt.0.0d0) then
2792           ii=0
2793           do i = nnt,nct-2 ! right? without parallel.
2794             do j=i+2,nct ! right?
2795 c         do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
2796 c           do j=i+2,nres ! ibid
2797 c         do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
2798 c           do j=i+2,nct ! ibid
2799               ii=ii+1
2800 c             write (iout,*) "k",k
2801 c             write (iout,*) "i",i," j",j," constr_homology",
2802 c    &                       constr_homology
2803               ires_homo(ii)=i
2804               jres_homo(ii)=j
2805 c
2806 c Attempt to replace dist(i,j) by its definition in ...
2807 c
2808               x12=c(1,i)-c(1,j)
2809               y12=c(2,i)-c(2,j)
2810               z12=c(3,i)-c(3,j)
2811               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2812               odl(k,ii)=distal
2813 c
2814 c             odl(k,ii)=dist(i,j)
2815 c             write (iout,*) "dist(",i,j,") =",dist(i,j)
2816 c             write (iout,*) "distal = ",distal
2817 c             write (iout,*) "odl(",k,ii,") =",odl(k,ii)
2818 c            write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2819 c    &                      "rescore(",k,j,") =",rescore(k,j)
2820 c
2821 c  Calculation of sigma from res sim
2822 c
2823 c             if (odl(k,ii).le.6.0d0) then
2824 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
2825 c    Other functional forms possible depending on odl(k,ii), eg.
2826 c
2827             if (odl(k,ii).le.dist_cut) then
2828               sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
2829 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
2830             else
2831               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
2832      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2833
2834 c   Following expr replaced by a positive exp argument
2835 c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2836 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
2837
2838 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
2839 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
2840             endif
2841 c
2842               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
2843 c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
2844 c
2845 c             sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
2846 c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of use of res. similarity
2847             enddo
2848 c           read (ientin,*) sigma_odl(k,ii) ! 1st variant
2849           enddo
2850 c         lim_odl=ii
2851 c         if (constr_homology.gt.0) call homology_partition
2852         endif
2853 c
2854 c     Theta, dihedral and SC retraints
2855 c
2856         if (waga_angle.gt.0.0d0) then
2857 c         open (ientin,file=tpl_k_sigma_dih,status='old')
2858 c         do irec=1,maxres-3 ! loop for reading sigma_dih
2859 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2860 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2861 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2862 c    &                            sigma_dih(k,i+nnt-1)
2863 c         enddo
2864 c1402   continue
2865 c         close (ientin)
2866           do i = nnt+3,nct ! right? without parallel.
2867 c         do i=1,nres ! alternative for bounds acc to readpdb?
2868 c         do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
2869 c         do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
2870             dih(k,i)=phiref(i) ! right?
2871 c           read (ientin,*) sigma_dih(k,i) ! original variant
2872 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
2873 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2874 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2875 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
2876 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
2877
2878             sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
2879      &                     rescore(k,i-2)+rescore(k,i-3)  !  right expression ?
2880 c
2881 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
2882 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2883 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2884 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
2885 c   Instead of res sim other local measure of b/b str reliability possible
2886             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2887 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2888             if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
2889 c           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
2890           enddo
2891         endif
2892
2893         if (waga_theta.gt.0.0d0) then
2894 c         open (ientin,file=tpl_k_sigma_theta,status='old')
2895 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2896 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2897 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2898 c    &                              sigma_theta(k,i+nnt-1)
2899 c         enddo
2900 c1403   continue
2901 c         close (ientin)
2902
2903           do i = nnt+2,nct ! right? without parallel.
2904 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
2905 c         do i=ithet_start,ithet_end ! with FG parallel.
2906              thetatpl(k,i)=thetaref(i)
2907 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2908 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
2909 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2910 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
2911 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
2912              sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
2913      &                        rescore(k,i-2) !  right expression ?
2914              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2915
2916 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2917 c                             rescore(k,i-2) !  right expression ?
2918 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2919              if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
2920           enddo
2921         endif
2922
2923         if (waga_d.gt.0.0d0) then
2924 c       open (ientin,file=tpl_k_sigma_d,status='old')
2925 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2926 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2927 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2928 c    &                          sigma_d(k,i+nnt-1)
2929 c         enddo
2930 c1404   continue
2931           close (ientin)
2932
2933           do i = nnt,nct ! right? without parallel.
2934 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
2935 c         do i=loc_start,loc_end ! with FG parallel.
2936              if (itype(i).eq.10) goto 1 ! right?
2937                xxtpl(k,i)=xxref(i)
2938                yytpl(k,i)=yyref(i)
2939                zztpl(k,i)=zzref(i)
2940 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2941 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2942 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2943 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
2944                sigma_d(k,i)=rescore(k,i) !  right expression ?
2945                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2946
2947 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
2948 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2949 c              read (ientin,*) sigma_d(k,i) ! 1st variant
2950                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
2951     1     continue
2952           enddo
2953         endif
2954         close(ientin)
2955       enddo
2956       if (waga_dist.gt.0.0d0) lim_odl=ii
2957       if (constr_homology.gt.0) call homology_partition
2958       if (constr_homology.gt.0) call init_int_table
2959       write (iout,*) "homology_partition: lim_theta= ",lim_theta,
2960      & "lim_xx=",lim_xx
2961 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2962 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2963 c
2964 c Print restraints
2965 c
2966       if (.not.lprn) return
2967       write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2968       write (iout,*) "Distance restraints from templates"
2969       do ii=1,lim_odl
2970        write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
2971      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
2972       enddo
2973       write (iout,*) "Dihedral angle restraints from templates"
2974       do i=nnt+3,lim_dih
2975         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
2976      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2977       enddo
2978       write (iout,*) "Virtual-bond angle restraints from templates"
2979       do i=nnt+2,lim_theta
2980         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
2981      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2982       enddo
2983       write (iout,*) "SC restraints from templates"
2984       do i=nnt,lim_xx
2985         write(iout,'(i5,10(4f8.2,4x))') i,
2986      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
2987      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
2988       enddo
2989
2990 c -----------------------------------------------------------------
2991       return
2992       end
2993 c----------------------------------------------------------------------
2994
2995 #ifdef WINIFL
2996       subroutine flush(iu)
2997       return
2998       end
2999 #endif
3000 #ifdef AIX
3001       subroutine flush(iu)
3002       call flush_(iu)
3003       return
3004       end
3005 #endif
3006
3007 c------------------------------------------------------------------------------
3008       subroutine copy_to_tmp(source)
3009       include "DIMENSIONS"
3010       include "COMMON.IOUNITS"
3011       character*(*) source
3012       character* 256 tmpfile
3013       integer ilen
3014       external ilen
3015       logical ex
3016       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3017       inquire(file=tmpfile,exist=ex)
3018       if (ex) then
3019         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3020      &   " to temporary directory..."
3021         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3022         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3023       endif
3024       return
3025       end
3026 c------------------------------------------------------------------------------
3027       subroutine move_from_tmp(source)
3028       include "DIMENSIONS"
3029       include "COMMON.IOUNITS"
3030       character*(*) source
3031       integer ilen
3032       external ilen
3033       write (*,*) "Moving ",source(:ilen(source)),
3034      & " from temporary directory to working directory"
3035       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3036       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3037       return
3038       end
3039 c------------------------------------------------------------------------------
3040       subroutine random_init(seed)
3041 C
3042 C Initialize random number generator
3043 C
3044       implicit real*8 (a-h,o-z)
3045       include 'DIMENSIONS'
3046 #ifdef AMD64
3047       integer*8 iseedi8
3048 #endif
3049 #ifdef MPI
3050       include 'mpif.h'
3051       logical OKRandom, prng_restart
3052       real*8  r1
3053       integer iseed_array(4)
3054 #endif
3055       include 'COMMON.IOUNITS'
3056       include 'COMMON.TIME1'
3057       include 'COMMON.THREAD'
3058       include 'COMMON.SBRIDGE'
3059       include 'COMMON.CONTROL'
3060       include 'COMMON.MCM'
3061       include 'COMMON.MAP'
3062       include 'COMMON.HEADER'
3063 csa      include 'COMMON.CSA'
3064       include 'COMMON.CHAIN'
3065       include 'COMMON.MUCA'
3066       include 'COMMON.MD'
3067       include 'COMMON.FFIELD'
3068       include 'COMMON.SETUP'
3069       iseed=-dint(dabs(seed))
3070       if (iseed.eq.0) then
3071         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3072      &    'Random seed undefined. The program will stop.'
3073         write (*,'(/80(1h*)/20x,a/80(1h*))')
3074      &    'Random seed undefined. The program will stop.'
3075 #ifdef MPI
3076         call mpi_finalize(mpi_comm_world,ierr)
3077 #endif
3078         stop 'Bad random seed.'
3079       endif
3080 #ifdef MPI
3081       if (fg_rank.eq.0) then
3082       seed=seed*(me+1)+1
3083 #ifdef AMD64
3084       iseedi8=dint(seed)
3085       if(me.eq.king .or. .not. out1file)
3086      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3087       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3088       OKRandom = prng_restart(me,iseedi8)
3089 #else
3090       do i=1,4
3091        tmp=65536.0d0**(4-i)
3092        iseed_array(i) = dint(seed/tmp)
3093        seed=seed-iseed_array(i)*tmp
3094       enddo
3095       if(me.eq.king .or. .not. out1file)
3096      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3097      &                 (iseed_array(i),i=1,4)
3098       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3099      &                 (iseed_array(i),i=1,4)
3100       OKRandom = prng_restart(me,iseed_array)
3101 #endif
3102       if (OKRandom) then
3103         r1=ran_number(0.0D0,1.0D0)
3104         if(me.eq.king .or. .not. out1file)
3105      &   write (iout,*) 'ran_num',r1
3106         if (r1.lt.0.0d0) OKRandom=.false.
3107       endif
3108       if (.not.OKRandom) then
3109         write (iout,*) 'PRNG IS NOT WORKING!!!'
3110         print *,'PRNG IS NOT WORKING!!!'
3111         if (me.eq.0) then
3112          call flush(iout)
3113          call mpi_abort(mpi_comm_world,error_msg,ierr)
3114          stop
3115         else
3116          write (iout,*) 'too many processors for parallel prng'
3117          write (*,*) 'too many processors for parallel prng'
3118          call flush(iout)
3119          stop
3120         endif
3121       endif
3122       endif
3123 #else
3124       call vrndst(iseed)
3125       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3126 #endif
3127       return
3128       end