source/unres/src_MD-M/readrtns_CSA.F

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