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 c      write(iout,*) "bordliptop=",bordliptop
 245 c      write(iout,*) "bordlipbot=",bordlipbot
 246 c      write(iout,*) "bufliptop=",bufliptop
 247 c      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 c        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 c from old chainbuild
1111 C
1112 C Define the origin and orientation of the coordinate system and locate the
1113 C first three CA's and SC(2).
1114 C
1115       call orig_frame
1116 *
1117 * Build the alpha-carbon chain.
1118 *
1119       do i=4,nres
1120         call locate_next_res(i)
1121       enddo
1122 C
1123 C First and last SC must coincide with the corresponding CA.
1124 C
1125       do j=1,3
1126         dc(j,nres+1)=0.0D0
1127         dc_norm(j,nres+1)=0.0D0
1128         dc(j,nres+nres)=0.0D0
1129         dc_norm(j,nres+nres)=0.0D0
1130         c(j,nres+1)=c(j,1)
1131         c(j,nres+nres)=c(j,nres)
1132       enddo
1133 C
1134 C Define the origin and orientation of the coordinate system and locate the
1135 C first three CA's and SC(2).
1136 C
1137       call orig_frame
1138 *
1139 * Build the alpha-carbon chain.
1140 *
1141       do i=4,nres
1142         call locate_next_res(i)
1143       enddo
1144 C
1145 C First and last SC must coincide with the corresponding CA.
1146 C
1147       do j=1,3
1148         dc(j,nres+1)=0.0D0
1149         dc_norm(j,nres+1)=0.0D0
1150         dc(j,nres+nres)=0.0D0
1151         dc_norm(j,nres+nres)=0.0D0
1152         c(j,nres+1)=c(j,1)
1153         c(j,nres+nres)=c(j,nres)
1154       enddo
1155
1156 c
1157         else
1158           if(me.eq.king.or..not.out1file)
1159      &     write (iout,'(a)') 'Random-generated initial geometry.'
1160
1161
1162 #ifdef MPI
1163           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1164      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1165 #endif
1166             do itrial=1,100
1167               itmp=1
1168               call gen_rand_conf(itmp,*30)
1169               goto 40
1170    30         write (iout,*) 'Failed to generate random conformation',
1171      &          ', itrial=',itrial
1172               write (*,*) 'Processor:',me,
1173      &          ' Failed to generate random conformation',
1174      &          ' itrial=',itrial
1175               call intout
1176
1177 #ifdef AIX
1178               call flush_(iout)
1179 #else
1180               call flush(iout)
1181 #endif
1182             enddo
1183             write (iout,'(a,i3,a)') 'Processor:',me,
1184      &        ' error in generating random conformation.'
1185             write (*,'(a,i3,a)') 'Processor:',me,
1186      &        ' error in generating random conformation.'
1187             call flush(iout)
1188 #ifdef MPI
1189             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1190    40       continue
1191           endif
1192 #else
1193           write (*,'(a)')
1194      &      ' error in generating random conformation.'
1195           stop
1196    40     continue
1197 #endif
1198         endif
1199       elseif (modecalc.eq.4) then
1200         read (inp,'(a)') intinname
1201         open (intin,file=intinname,status='old',err=333)
1202         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1203      &  write (iout,'(a)') 'intinname',intinname
1204         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1205         goto 334
1206   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1207 #ifdef MPI
1208         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1209 #endif
1210         stop 'Error opening angle file.'
1211   334   continue
1212
1213       endif
1214 C Generate distance constraints, if the PDB structure is to be regularized.
1215       if (nthread.gt.0) then
1216         call read_threadbase
1217       endif
1218       call setup_var
1219       if (me.eq.king .or. .not. out1file)
1220      & call intout
1221       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1222         write (iout,'(/a,i3,a)')
1223      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1224         write (iout,'(20i4)') (iss(i),i=1,ns)
1225        if (dyn_ss) then
1226           write(iout,*)"Running with dynamic disulfide-bond formation"
1227        else
1228         write (iout,'(/a/)') 'Pre-formed links are:'
1229         do i=1,nss
1230           i1=ihpb(i)-nres
1231           i2=jhpb(i)-nres
1232           it1=itype(i1)
1233           it2=itype(i2)
1234           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1235      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1236      &    ebr,forcon(i)
1237         enddo
1238         write (iout,'(a)')
1239        endif
1240       endif
1241       if (ns.gt.0.and.dyn_ss) then
1242           do i=nss+1,nhpb
1243             ihpb(i-nss)=ihpb(i)
1244             jhpb(i-nss)=jhpb(i)
1245             forcon(i-nss)=forcon(i)
1246             dhpb(i-nss)=dhpb(i)
1247           enddo
1248           nhpb=nhpb-nss
1249           nss=0
1250           call hpb_partition
1251           do i=1,ns
1252             dyn_ss_mask(iss(i))=.true.
1253           enddo
1254       endif
1255       if (i2ndstr.gt.0) call secstrp2dihc
1256 c      call geom_to_var(nvar,x)
1257 c      call etotal(energia(0))
1258 c      call enerprint(energia(0))
1259 c      call briefout(0,etot)
1260 c      stop
1261 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1262 cd    write (iout,'(a)') 'Variable list:'
1263 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1264 #ifdef MPI
1265       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1266      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1267      &  'Processor',myrank,': end reading molecular data.'
1268 #endif
1269       return
1270       end
1271 c--------------------------------------------------------------------------
1272       logical function seq_comp(itypea,itypeb,length)
1273       implicit none
1274       integer length,itypea(length),itypeb(length)
1275       integer i
1276       do i=1,length
1277         if (itypea(i).ne.itypeb(i)) then
1278           seq_comp=.false.
1279           return
1280         endif
1281       enddo
1282       seq_comp=.true.
1283       return
1284       end
1285 c-----------------------------------------------------------------------------
1286       subroutine read_bridge
1287 C Read information about disulfide bridges.
1288       implicit real*8 (a-h,o-z)
1289       include 'DIMENSIONS'
1290 #ifdef MPI
1291       include 'mpif.h'
1292 #endif
1293       include 'COMMON.IOUNITS'
1294       include 'COMMON.GEO'
1295       include 'COMMON.VAR'
1296       include 'COMMON.INTERACT'
1297       include 'COMMON.LOCAL'
1298       include 'COMMON.NAMES'
1299       include 'COMMON.CHAIN'
1300       include 'COMMON.FFIELD'
1301       include 'COMMON.SBRIDGE'
1302       include 'COMMON.HEADER'
1303       include 'COMMON.CONTROL'
1304       include 'COMMON.DBASE'
1305       include 'COMMON.THREAD'
1306       include 'COMMON.TIME1'
1307       include 'COMMON.SETUP'
1308 C Read bridging residues.
1309       read (inp,*) ns,(iss(i),i=1,ns)
1310       print *,'ns=',ns
1311       if(me.eq.king.or..not.out1file)
1312      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1313 C Check whether the specified bridging residues are cystines.
1314       do i=1,ns
1315         if (itype(iss(i)).ne.1) then
1316           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1317      &   'Do you REALLY think that the residue ',
1318      &    restyp(itype(iss(i))),i,
1319      &   ' can form a disulfide bridge?!!!'
1320           write (*,'(2a,i3,a)')
1321      &   'Do you REALLY think that the residue ',
1322      &    restyp(itype(iss(i))),i,
1323      &   ' can form a disulfide bridge?!!!'
1324 #ifdef MPI
1325          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1326          stop
1327 #endif
1328         endif
1329       enddo
1330 C Read preformed bridges.
1331       if (ns.gt.0) then
1332       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1333       if(fg_rank.eq.0)
1334      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1335       if (nss.gt.0) then
1336         nhpb=nss
1337 C Check if the residues involved in bridges are in the specified list of
1338 C bridging residues.
1339         do i=1,nss
1340           do j=1,i-1
1341             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1342      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1343               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1344      &      ' contains residues present in other pairs.'
1345               write (*,'(a,i3,a)') 'Disulfide pair',i,
1346      &      ' contains residues present in other pairs.'
1347 #ifdef MPI
1348               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1349               stop
1350 #endif
1351             endif
1352           enddo
1353           do j=1,ns
1354             if (ihpb(i).eq.iss(j)) goto 10
1355           enddo
1356           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1357    10     continue
1358           do j=1,ns
1359             if (jhpb(i).eq.iss(j)) goto 20
1360           enddo
1361           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1362    20     continue
1363           dhpb(i)=dbr
1364           forcon(i)=fbr
1365         enddo
1366         do i=1,nss
1367           ihpb(i)=ihpb(i)+nres
1368           jhpb(i)=jhpb(i)+nres
1369         enddo
1370       endif
1371       endif
1372       return
1373       end
1374 c----------------------------------------------------------------------------
1375       subroutine read_x(kanal,*)
1376       implicit real*8 (a-h,o-z)
1377       include 'DIMENSIONS'
1378       include 'COMMON.GEO'
1379       include 'COMMON.VAR'
1380       include 'COMMON.CHAIN'
1381       include 'COMMON.IOUNITS'
1382       include 'COMMON.CONTROL'
1383       include 'COMMON.LOCAL'
1384       include 'COMMON.INTERACT'
1385 c Read coordinates from input
1386 c
1387       read(kanal,'(8f10.5)',end=10,err=10)
1388      &  ((c(l,k),l=1,3),k=1,nres),
1389      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1390       do j=1,3
1391         c(j,nres+1)=c(j,1)
1392         c(j,2*nres)=c(j,nres)
1393       enddo
1394       call int_from_cart1(.false.)
1395       do i=1,nres-1
1396         do j=1,3
1397           dc(j,i)=c(j,i+1)-c(j,i)
1398           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1399         enddo
1400       enddo
1401       do i=nnt,nct
1402         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1403           do j=1,3
1404             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1405             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1406           enddo
1407         endif
1408       enddo
1409
1410       return
1411    10 return1
1412       end
1413 c----------------------------------------------------------------------------
1414       subroutine read_threadbase
1415       implicit real*8 (a-h,o-z)
1416       include 'DIMENSIONS'
1417       include 'COMMON.IOUNITS'
1418       include 'COMMON.GEO'
1419       include 'COMMON.VAR'
1420       include 'COMMON.INTERACT'
1421       include 'COMMON.LOCAL'
1422       include 'COMMON.NAMES'
1423       include 'COMMON.CHAIN'
1424       include 'COMMON.FFIELD'
1425       include 'COMMON.SBRIDGE'
1426       include 'COMMON.HEADER'
1427       include 'COMMON.CONTROL'
1428       include 'COMMON.DBASE'
1429       include 'COMMON.THREAD'
1430       include 'COMMON.TIME1'
1431 C Read pattern database for threading.
1432       read (icbase,*) nseq
1433       do i=1,nseq
1434         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1435      &   nres_base(2,i),nres_base(3,i)
1436         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1437      &   nres_base(1,i))
1438 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1439 c    &   nres_base(2,i),nres_base(3,i)
1440 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1441 c    &   nres_base(1,i))
1442       enddo
1443       close (icbase)
1444       if (weidis.eq.0.0D0) weidis=0.1D0
1445       do i=nnt,nct
1446         do j=i+2,nct
1447           nhpb=nhpb+1
1448           ihpb(nhpb)=i
1449           jhpb(nhpb)=j
1450           forcon(nhpb)=weidis
1451         enddo
1452       enddo
1453       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1454       write (iout,'(a,i5)') 'nexcl: ',nexcl
1455       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1456       return
1457       end
1458 c------------------------------------------------------------------------------
1459       subroutine setup_var
1460       implicit real*8 (a-h,o-z)
1461       include 'DIMENSIONS'
1462       include 'COMMON.IOUNITS'
1463       include 'COMMON.GEO'
1464       include 'COMMON.VAR'
1465       include 'COMMON.INTERACT'
1466       include 'COMMON.LOCAL'
1467       include 'COMMON.NAMES'
1468       include 'COMMON.CHAIN'
1469       include 'COMMON.FFIELD'
1470       include 'COMMON.SBRIDGE'
1471       include 'COMMON.HEADER'
1472       include 'COMMON.CONTROL'
1473       include 'COMMON.DBASE'
1474       include 'COMMON.THREAD'
1475       include 'COMMON.TIME1'
1476 C Set up variable list.
1477       ntheta=nres-2
1478       nphi=nres-3
1479       nvar=ntheta+nphi
1480       nside=0
1481       do i=2,nres-1
1482         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1483           nside=nside+1
1484           ialph(i,1)=nvar+nside
1485           ialph(nside,2)=i
1486         endif
1487       enddo
1488       if (indphi.gt.0) then
1489         nvar=nphi
1490       else if (indback.gt.0) then
1491         nvar=nphi+ntheta
1492       else
1493         nvar=nvar+2*nside
1494       endif
1495 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1496       return
1497       end
1498 c----------------------------------------------------------------------------
1499       subroutine gen_dist_constr
1500 C Generate CA distance constraints.
1501       implicit real*8 (a-h,o-z)
1502       include 'DIMENSIONS'
1503       include 'COMMON.IOUNITS'
1504       include 'COMMON.GEO'
1505       include 'COMMON.VAR'
1506       include 'COMMON.INTERACT'
1507       include 'COMMON.LOCAL'
1508       include 'COMMON.NAMES'
1509       include 'COMMON.CHAIN'
1510       include 'COMMON.FFIELD'
1511       include 'COMMON.SBRIDGE'
1512       include 'COMMON.HEADER'
1513       include 'COMMON.CONTROL'
1514       include 'COMMON.DBASE'
1515       include 'COMMON.THREAD'
1516       include 'COMMON.TIME1'
1517       dimension itype_pdb(maxres)
1518       common /pizda/ itype_pdb
1519       character*2 iden
1520 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1521 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1522 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1523 cd     & ' nsup',nsup
1524       do i=nstart_sup,nstart_sup+nsup-1
1525 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1526 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1527         do j=i+2,nstart_sup+nsup-1
1528           nhpb=nhpb+1
1529           ihpb(nhpb)=i+nstart_seq-nstart_sup
1530           jhpb(nhpb)=j+nstart_seq-nstart_sup
1531           forcon(nhpb)=weidis
1532           dhpb(nhpb)=dist(i,j)
1533         enddo
1534       enddo
1535 cd      write (iout,'(a)') 'Distance constraints:'
1536 cd      do i=nss+1,nhpb
1537 cd        ii=ihpb(i)
1538 cd        jj=jhpb(i)
1539 cd        iden='CA'
1540 cd        if (ii.gt.nres) then
1541 cd          iden='SC'
1542 cd          ii=ii-nres
1543 cd          jj=jj-nres
1544 cd        endif
1545 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1546 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1547 cd     &  dhpb(i),forcon(i)
1548 cd      enddo
1549       return
1550       end
1551 c----------------------------------------------------------------------------
1552       subroutine map_read
1553       implicit real*8 (a-h,o-z)
1554       include 'DIMENSIONS'
1555       include 'COMMON.MAP'
1556       include 'COMMON.IOUNITS'
1557       character*3 angid(4) /'THE','PHI','ALP','OME'/
1558       character*80 mapcard,ucase
1559       do imap=1,nmap
1560         read (inp,'(a)') mapcard
1561         mapcard=ucase(mapcard)
1562         if (index(mapcard,'PHI').gt.0) then
1563           kang(imap)=1
1564         else if (index(mapcard,'THE').gt.0) then
1565           kang(imap)=2
1566         else if (index(mapcard,'ALP').gt.0) then
1567           kang(imap)=3
1568         else if (index(mapcard,'OME').gt.0) then
1569           kang(imap)=4
1570         else
1571           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1572           stop 'Error - illegal variable spec in MAP card.'
1573         endif
1574         call readi (mapcard,'RES1',res1(imap),0)
1575         call readi (mapcard,'RES2',res2(imap),0)
1576         if (res1(imap).eq.0) then
1577           res1(imap)=res2(imap)
1578         else if (res2(imap).eq.0) then
1579           res2(imap)=res1(imap)
1580         endif
1581         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1582           write (iout,'(a)')
1583      &    'Error - illegal definition of variable group in MAP.'
1584           stop 'Error - illegal definition of variable group in MAP.'
1585         endif
1586         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1587         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1588         call readi(mapcard,'NSTEP',nstep(imap),0)
1589         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1590           write (iout,'(a)')
1591      &     'Illegal boundary and/or step size specification in MAP.'
1592           stop 'Illegal boundary and/or step size specification in MAP.'
1593         endif
1594       enddo ! imap
1595       return
1596       end
1597 c----------------------------------------------------------------------------
1598       subroutine csaread
1599       implicit real*8 (a-h,o-z)
1600       include 'DIMENSIONS'
1601       include 'COMMON.IOUNITS'
1602       include 'COMMON.GEO'
1603       include 'COMMON.CSA'
1604       include 'COMMON.BANK'
1605       include 'COMMON.CONTROL'
1606       character*80 ucase
1607       character*620 mcmcard
1608       call card_concat(mcmcard)
1609
1610       call readi(mcmcard,'NCONF',nconf,50)
1611       call readi(mcmcard,'NADD',nadd,0)
1612       call readi(mcmcard,'JSTART',jstart,1)
1613       call readi(mcmcard,'JEND',jend,1)
1614       call readi(mcmcard,'NSTMAX',nstmax,500000)
1615       call readi(mcmcard,'N0',n0,1)
1616       call readi(mcmcard,'N1',n1,6)
1617       call readi(mcmcard,'N2',n2,4)
1618       call readi(mcmcard,'N3',n3,0)
1619       call readi(mcmcard,'N4',n4,0)
1620       call readi(mcmcard,'N5',n5,0)
1621       call readi(mcmcard,'N6',n6,10)
1622       call readi(mcmcard,'N7',n7,0)
1623       call readi(mcmcard,'N8',n8,0)
1624       call readi(mcmcard,'N9',n9,0)
1625       call readi(mcmcard,'N14',n14,0)
1626       call readi(mcmcard,'N15',n15,0)
1627       call readi(mcmcard,'N16',n16,0)
1628       call readi(mcmcard,'N17',n17,0)
1629       call readi(mcmcard,'N18',n18,0)
1630
1631       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1632
1633       call readi(mcmcard,'NDIFF',ndiff,2)
1634       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1635       call readi(mcmcard,'IS1',is1,1)
1636       call readi(mcmcard,'IS2',is2,8)
1637       call readi(mcmcard,'NRAN0',nran0,4)
1638       call readi(mcmcard,'NRAN1',nran1,2)
1639       call readi(mcmcard,'IRR',irr,1)
1640       call readi(mcmcard,'NSEED',nseed,20)
1641       call readi(mcmcard,'NTOTAL',ntotal,10000)
1642       call reada(mcmcard,'CUT1',cut1,2.0d0)
1643       call reada(mcmcard,'CUT2',cut2,5.0d0)
1644       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1645       call readi(mcmcard,'ICMAX',icmax,3)
1646       call readi(mcmcard,'IRESTART',irestart,0)
1647 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1648       ntbankm=0
1649 c!bankt
1650       call reada(mcmcard,'DELE',dele,20.0d0)
1651       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1652       call readi(mcmcard,'IREF',iref,0)
1653       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1654       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1655       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1656       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1657       write (iout,*) "NCONF_IN",nconf_in
1658       return
1659       end
1660 c----------------------------------------------------------------------------
1661 cfmc      subroutine mcmfread
1662 cfmc      implicit real*8 (a-h,o-z)
1663 cfmc      include 'DIMENSIONS'
1664 cfmc      include 'COMMON.MCMF'
1665 cfmc      include 'COMMON.IOUNITS'
1666 cfmc      include 'COMMON.GEO'
1667 cfmc      character*80 ucase
1668 cfmc      character*620 mcmcard
1669 cfmc      call card_concat(mcmcard)
1670 cfmc
1671 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1672 cfmc      write(iout,*)'MAXRANT=',maxrant
1673 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1674 cfmc      write(iout,*)'MAXFAM=',maxfam
1675 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1676 cfmc      write(iout,*)'NNET1=',nnet1
1677 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1678 cfmc      write(iout,*)'NNET2=',nnet2
1679 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1680 cfmc      write(iout,*)'NNET3=',nnet3
1681 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1682 cfmc      write(iout,*)'ILASTT=',ilastt
1683 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1684 cfmc      write(iout,*)'MAXSTR=',maxstr
1685 cfmc      maxstr_f=maxstr/maxfam
1686 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1687 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1688 cfmc      write(iout,*)'NMCMF=',nmcmf
1689 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1690 cfmc      write(iout,*)'IFOCUS=',ifocus
1691 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1692 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1693 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1694 cfmc      write(iout,*)'INTPRT=',intprt
1695 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1696 cfmc      write(iout,*)'IPRT=',iprt
1697 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1698 cfmc      write(iout,*)'IMAXTR=',imaxtr
1699 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1700 cfmc      write(iout,*)'MAXEVEN=',maxeven
1701 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1702 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1703 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1704 cfmc      write(iout,*)'INIMIN=',inimin
1705 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1706 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1707 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1708 cfmc      write(iout,*)'NTHREAD=',nthread
1709 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1710 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1711 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1712 cfmc      write(iout,*)'MAXPERT=',maxpert
1713 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1714 cfmc      write(iout,*)'IRMSD=',irmsd
1715 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1716 cfmc      write(iout,*)'DENEMIN=',denemin
1717 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1718 cfmc      write(iout,*)'RCUT1S=',rcut1s
1719 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1720 cfmc      write(iout,*)'RCUT1E=',rcut1e
1721 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1722 cfmc      write(iout,*)'RCUT2S=',rcut2s
1723 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1724 cfmc      write(iout,*)'RCUT2E=',rcut2e
1725 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1726 cfmc      write(iout,*)'DPERT1=',d_pert1
1727 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1728 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1729 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1730 cfmc      write(iout,*)'DPERT2=',d_pert2
1731 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1732 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1733 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1734 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1735 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1736 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1737 cfmc      d_pert1=deg2rad*d_pert1
1738 cfmc      d_pert1a=deg2rad*d_pert1a
1739 cfmc      d_pert2=deg2rad*d_pert2
1740 cfmc      d_pert2a=deg2rad*d_pert2a
1741 cfmc      d_pert2b=deg2rad*d_pert2b
1742 cfmc      d_pert2c=deg2rad*d_pert2c
1743 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1744 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1745 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1746 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1747 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1748 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1749 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1750 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1751 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1752 cfmc      write(iout,*)'RCUTINI=',rcutini
1753 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1754 cfmc      write(iout,*)'GRAT=',grat
1755 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1756 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1757 cfmc
1758 cfmc      return
1759 cfmc      end
1760 c----------------------------------------------------------------------------
1761       subroutine mcmread
1762       implicit real*8 (a-h,o-z)
1763       include 'DIMENSIONS'
1764       include 'COMMON.MCM'
1765       include 'COMMON.MCE'
1766       include 'COMMON.IOUNITS'
1767       character*80 ucase
1768       character*320 mcmcard
1769       call card_concat(mcmcard)
1770       call readi(mcmcard,'MAXACC',maxacc,100)
1771       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1772       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1773       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1774       call readi(mcmcard,'MAXREPM',maxrepm,200)
1775       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1776       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1777       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1778       call reada(mcmcard,'E_UP',e_up,5.0D0)
1779       call reada(mcmcard,'DELTE',delte,0.1D0)
1780       call readi(mcmcard,'NSWEEP',nsweep,5)
1781       call readi(mcmcard,'NSTEPH',nsteph,0)
1782       call readi(mcmcard,'NSTEPC',nstepc,0)
1783       call reada(mcmcard,'TMIN',tmin,298.0D0)
1784       call reada(mcmcard,'TMAX',tmax,298.0D0)
1785       call readi(mcmcard,'NWINDOW',nwindow,0)
1786       call readi(mcmcard,'PRINT_MC',print_mc,0)
1787       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1788       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1789       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1790       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1791       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1792       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1793       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1794       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1795       if (nwindow.gt.0) then
1796         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1797         do i=1,nwindow
1798           winlen(i)=winend(i)-winstart(i)+1
1799         enddo
1800       endif
1801       if (tmax.lt.tmin) tmax=tmin
1802       if (tmax.eq.tmin) then
1803         nstepc=0
1804         nsteph=0
1805       endif
1806       if (nstepc.gt.0 .and. nsteph.gt.0) then
1807         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1808         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1809       endif
1810 C Probabilities of different move types
1811       sumpro_type(0)=0.0D0
1812       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1813       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1814       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1815       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1816       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1817       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1818       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1819       do i=1,MaxMoveType
1820         print *,'i',i,' sumprotype',sumpro_type(i)
1821         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1822         print *,'i',i,' sumprotype',sumpro_type(i)
1823       enddo
1824       return
1825       end
1826 c----------------------------------------------------------------------------
1827       subroutine read_minim
1828       implicit real*8 (a-h,o-z)
1829       include 'DIMENSIONS'
1830       include 'COMMON.MINIM'
1831       include 'COMMON.IOUNITS'
1832       include 'COMMON.CONTROL'
1833       include 'COMMON.SETUP'
1834       character*80 ucase
1835       character*320 minimcard
1836       call card_concat(minimcard)
1837       call readi(minimcard,'MAXMIN',maxmin,2000)
1838       call readi(minimcard,'MAXFUN',maxfun,5000)
1839       call readi(minimcard,'MINMIN',minmin,maxmin)
1840       call readi(minimcard,'MINFUN',minfun,maxmin)
1841       call reada(minimcard,'TOLF',tolf,1.0D-2)
1842       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1843       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1844       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1845       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1846 #ifdef MPI
1847       if (.not. out1file .or. me.eq.king) then
1848 #endif
1849       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1850      &         'Options in energy minimization:'
1851       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1852      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1853      & 'MinMin:',MinMin,' MinFun:',MinFun,
1854      & ' TolF:',TolF,' RTolF:',RTolF
1855 #ifdef MPI
1856       endif
1857 #endif
1858       return
1859       end
1860 c----------------------------------------------------------------------------
1861       subroutine read_angles(kanal,*)
1862       implicit real*8 (a-h,o-z)
1863       include 'DIMENSIONS'
1864       include 'COMMON.GEO'
1865       include 'COMMON.VAR'
1866       include 'COMMON.CHAIN'
1867       include 'COMMON.IOUNITS'
1868       include 'COMMON.CONTROL'
1869 c Read angles from input
1870 c
1871        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1872        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1873        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1874        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1875
1876        do i=1,nres
1877 c 9/7/01 avoid 180 deg valence angle
1878         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1879 c
1880         theta(i)=deg2rad*theta(i)
1881         phi(i)=deg2rad*phi(i)
1882         alph(i)=deg2rad*alph(i)
1883         omeg(i)=deg2rad*omeg(i)
1884        enddo
1885       return
1886    10 return1
1887       end
1888 c----------------------------------------------------------------------------
1889       subroutine reada(rekord,lancuch,wartosc,default)
1890       implicit none
1891       character*(*) rekord,lancuch
1892       double precision wartosc,default
1893       integer ilen,iread
1894       external ilen
1895       iread=index(rekord,lancuch)
1896       if (iread.eq.0) then
1897         wartosc=default
1898         return
1899       endif
1900       iread=iread+ilen(lancuch)+1
1901       read (rekord(iread:),*,err=10,end=10) wartosc
1902       return
1903   10  wartosc=default
1904       return
1905       end
1906 c----------------------------------------------------------------------------
1907       subroutine readi(rekord,lancuch,wartosc,default)
1908       implicit none
1909       character*(*) rekord,lancuch
1910       integer wartosc,default
1911       integer ilen,iread
1912       external ilen
1913       iread=index(rekord,lancuch)
1914       if (iread.eq.0) then
1915         wartosc=default
1916         return
1917       endif
1918       iread=iread+ilen(lancuch)+1
1919       read (rekord(iread:),*,err=10,end=10) wartosc
1920       return
1921   10  wartosc=default
1922       return
1923       end
1924 c----------------------------------------------------------------------------
1925       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1926       implicit none
1927       integer dim,i
1928       integer tablica(dim),default
1929       character*(*) rekord,lancuch
1930       character*80 aux
1931       integer ilen,iread
1932       external ilen
1933       do i=1,dim
1934         tablica(i)=default
1935       enddo
1936       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1937       if (iread.eq.0) return
1938       iread=iread+ilen(lancuch)+1
1939       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1940    10 return
1941       end
1942 c----------------------------------------------------------------------------
1943       subroutine multreada(rekord,lancuch,tablica,dim,default)
1944       implicit none
1945       integer dim,i
1946       double precision tablica(dim),default
1947       character*(*) rekord,lancuch
1948       character*80 aux
1949       integer ilen,iread
1950       external ilen
1951       do i=1,dim
1952         tablica(i)=default
1953       enddo
1954       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1955       if (iread.eq.0) return
1956       iread=iread+ilen(lancuch)+1
1957       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1958    10 return
1959       end
1960 c----------------------------------------------------------------------------
1961       subroutine openunits
1962       implicit real*8 (a-h,o-z)
1963       include 'DIMENSIONS'
1964 #ifdef MPI
1965       include 'mpif.h'
1966       character*16 form,nodename
1967       integer nodelen
1968 #endif
1969       include 'COMMON.SETUP'
1970       include 'COMMON.IOUNITS'
1971       include 'COMMON.MD'
1972       include 'COMMON.CONTROL'
1973       integer lenpre,lenpot,ilen,lentmp
1974       external ilen
1975       character*3 out1file_text,ucase
1976       character*3 ll
1977       external ucase
1978 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1979       call getenv_loc("PREFIX",prefix)
1980       pref_orig = prefix
1981       call getenv_loc("POT",pot)
1982       call getenv_loc("DIRTMP",tmpdir)
1983       call getenv_loc("CURDIR",curdir)
1984       call getenv_loc("OUT1FILE",out1file_text)
1985 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1986       out1file_text=ucase(out1file_text)
1987       if (out1file_text(1:1).eq."Y") then
1988         out1file=.true.
1989       else
1990         out1file=fg_rank.gt.0
1991       endif
1992       lenpre=ilen(prefix)
1993       lenpot=ilen(pot)
1994       lentmp=ilen(tmpdir)
1995       if (lentmp.gt.0) then
1996           write (*,'(80(1h!))')
1997           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1998           write (*,'(80(1h!))')
1999           write (*,*)"All output files will be on node /tmp directory."
2000 #ifdef MPI
2001         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2002         if (me.eq.king) then
2003           write (*,*) "The master node is ",nodename
2004         else if (fg_rank.eq.0) then
2005           write (*,*) "I am the CG slave node ",nodename
2006         else
2007           write (*,*) "I am the FG slave node ",nodename
2008         endif
2009 #endif
2010         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2011         lenpre = lentmp+lenpre+1
2012       endif
2013       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2014 C Get the names and open the input files
2015 #if defined(WINIFL) || defined(WINPGI)
2016       open(1,file=pref_orig(:ilen(pref_orig))//
2017      &  '.inp',status='old',readonly,shared)
2018        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2019 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2020 C Get parameter filenames and open the parameter files.
2021       call getenv_loc('BONDPAR',bondname)
2022       open (ibond,file=bondname,status='old',readonly,shared)
2023       call getenv_loc('THETPAR',thetname)
2024       open (ithep,file=thetname,status='old',readonly,shared)
2025 #ifndef CRYST_THETA
2026       call getenv_loc('THETPARPDB',thetname_pdb)
2027       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2028 #endif
2029       call getenv_loc('ROTPAR',rotname)
2030       open (irotam,file=rotname,status='old',readonly,shared)
2031 #ifndef CRYST_SC
2032       call getenv_loc('ROTPARPDB',rotname_pdb)
2033       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2034 #endif
2035       call getenv_loc('TORPAR',torname)
2036       open (itorp,file=torname,status='old',readonly,shared)
2037       call getenv_loc('TORDPAR',tordname)
2038       open (itordp,file=tordname,status='old',readonly,shared)
2039       call getenv_loc('FOURIER',fouriername)
2040       open (ifourier,file=fouriername,status='old',readonly,shared)
2041       call getenv_loc('ELEPAR',elename)
2042       open (ielep,file=elename,status='old',readonly,shared)
2043       call getenv_loc('SIDEPAR',sidename)
2044       open (isidep,file=sidename,status='old',readonly,shared)
2045 #elif (defined CRAY) || (defined AIX)
2046       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2047      &  action='read')
2048 c      print *,"Processor",myrank," opened file 1"
2049       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2050 c      print *,"Processor",myrank," opened file 9"
2051 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2052 C Get parameter filenames and open the parameter files.
2053       call getenv_loc('BONDPAR',bondname)
2054       open (ibond,file=bondname,status='old',action='read')
2055 c      print *,"Processor",myrank," opened file IBOND"
2056       call getenv_loc('THETPAR',thetname)
2057       open (ithep,file=thetname,status='old',action='read')
2058 c      print *,"Processor",myrank," opened file ITHEP"
2059 #ifndef CRYST_THETA
2060       call getenv_loc('THETPARPDB',thetname_pdb)
2061       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2062 #endif
2063       call getenv_loc('ROTPAR',rotname)
2064       open (irotam,file=rotname,status='old',action='read')
2065 c      print *,"Processor",myrank," opened file IROTAM"
2066 #ifndef CRYST_SC
2067       call getenv_loc('ROTPARPDB',rotname_pdb)
2068       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2069 #endif
2070       call getenv_loc('TORPAR',torname)
2071       open (itorp,file=torname,status='old',action='read')
2072 c      print *,"Processor",myrank," opened file ITORP"
2073       call getenv_loc('TORDPAR',tordname)
2074       open (itordp,file=tordname,status='old',action='read')
2075 c      print *,"Processor",myrank," opened file ITORDP"
2076       call getenv_loc('SCCORPAR',sccorname)
2077       open (isccor,file=sccorname,status='old',action='read')
2078 c      print *,"Processor",myrank," opened file ISCCOR"
2079       call getenv_loc('FOURIER',fouriername)
2080       open (ifourier,file=fouriername,status='old',action='read')
2081 c      print *,"Processor",myrank," opened file IFOURIER"
2082       call getenv_loc('ELEPAR',elename)
2083       open (ielep,file=elename,status='old',action='read')
2084 c      print *,"Processor",myrank," opened file IELEP"
2085       call getenv_loc('SIDEPAR',sidename)
2086       open (isidep,file=sidename,status='old',action='read')
2087 c      print *,"Processor",myrank," opened file ISIDEP"
2088 c      print *,"Processor",myrank," opened parameter files"
2089 #elif (defined G77)
2090       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2091       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2092 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2093 C Get parameter filenames and open the parameter files.
2094       call getenv_loc('BONDPAR',bondname)
2095       open (ibond,file=bondname,status='old')
2096       call getenv_loc('THETPAR',thetname)
2097       open (ithep,file=thetname,status='old')
2098 #ifndef CRYST_THETA
2099       call getenv_loc('THETPARPDB',thetname_pdb)
2100       open (ithep_pdb,file=thetname_pdb,status='old')
2101 #endif
2102       call getenv_loc('ROTPAR',rotname)
2103       open (irotam,file=rotname,status='old')
2104 #ifndef CRYST_SC
2105       call getenv_loc('ROTPARPDB',rotname_pdb)
2106       open (irotam_pdb,file=rotname_pdb,status='old')
2107 #endif
2108       call getenv_loc('TORPAR',torname)
2109       open (itorp,file=torname,status='old')
2110       call getenv_loc('TORDPAR',tordname)
2111       open (itordp,file=tordname,status='old')
2112       call getenv_loc('SCCORPAR',sccorname)
2113       open (isccor,file=sccorname,status='old')
2114       call getenv_loc('FOURIER',fouriername)
2115       open (ifourier,file=fouriername,status='old')
2116       call getenv_loc('ELEPAR',elename)
2117       open (ielep,file=elename,status='old')
2118       call getenv_loc('SIDEPAR',sidename)
2119       open (isidep,file=sidename,status='old')
2120 #else
2121       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2122      &  readonly)
2123        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2124 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2125 C Get parameter filenames and open the parameter files.
2126       call getenv_loc('BONDPAR',bondname)
2127       open (ibond,file=bondname,status='old',readonly)
2128       call getenv_loc('THETPAR',thetname)
2129       open (ithep,file=thetname,status='old',readonly)
2130       call getenv_loc('ROTPAR',rotname)
2131       open (irotam,file=rotname,status='old',readonly)
2132       call getenv_loc('TORPAR',torname)
2133       open (itorp,file=torname,status='old',readonly)
2134       call getenv_loc('TORDPAR',tordname)
2135       open (itordp,file=tordname,status='old',readonly)
2136       call getenv_loc('SCCORPAR',sccorname)
2137       open (isccor,file=sccorname,status='old',readonly)
2138 #ifndef CRYST_THETA
2139       call getenv_loc('THETPARPDB',thetname_pdb)
2140 c      print *,"thetname_pdb ",thetname_pdb
2141       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2142 c      print *,ithep_pdb," opened"
2143 #endif
2144       call getenv_loc('FOURIER',fouriername)
2145       open (ifourier,file=fouriername,status='old',readonly)
2146       call getenv_loc('ELEPAR',elename)
2147       open (ielep,file=elename,status='old',readonly)
2148       call getenv_loc('SIDEPAR',sidename)
2149       open (isidep,file=sidename,status='old',readonly)
2150       call getenv_loc('LIPTRANPAR',liptranname)
2151       open (iliptranpar,file=liptranname,status='old',action='read')
2152 #ifndef CRYST_SC
2153       call getenv_loc('ROTPARPDB',rotname_pdb)
2154       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2155 #endif
2156 #endif
2157 #ifndef OLDSCP
2158 C
2159 C 8/9/01 In the newest version SCp interaction constants are read from a file
2160 C Use -DOLDSCP to use hard-coded constants instead.
2161 C
2162       call getenv_loc('SCPPAR',scpname)
2163 #if defined(WINIFL) || defined(WINPGI)
2164       open (iscpp,file=scpname,status='old',readonly,shared)
2165 #elif (defined CRAY)  || (defined AIX)
2166       open (iscpp,file=scpname,status='old',action='read')
2167 #elif (defined G77)
2168       open (iscpp,file=scpname,status='old')
2169 #else
2170       open (iscpp,file=scpname,status='old',readonly)
2171 #endif
2172 #endif
2173       call getenv_loc('PATTERN',patname)
2174 #if defined(WINIFL) || defined(WINPGI)
2175       open (icbase,file=patname,status='old',readonly,shared)
2176 #elif (defined CRAY)  || (defined AIX)
2177       open (icbase,file=patname,status='old',action='read')
2178 #elif (defined G77)
2179       open (icbase,file=patname,status='old')
2180 #else
2181       open (icbase,file=patname,status='old',readonly)
2182 #endif
2183 #ifdef MPI
2184 C Open output file only for CG processes
2185 c      print *,"Processor",myrank," fg_rank",fg_rank
2186       if (fg_rank.eq.0) then
2187
2188       if (nodes.eq.1) then
2189         npos=3
2190       else
2191         npos = dlog10(dfloat(nodes-1))+1
2192       endif
2193       if (npos.lt.3) npos=3
2194       write (liczba,'(i1)') npos
2195       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2196      &  //')'
2197       write (liczba,form) me
2198       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2199      &  liczba(:ilen(liczba))
2200       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2201      &  //'.int'
2202       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2203      &  //'.pdb'
2204       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2205      &  liczba(:ilen(liczba))//'.mol2'
2206       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2207      &  liczba(:ilen(liczba))//'.stat'
2208       if (lentmp.gt.0)
2209      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2210      &      //liczba(:ilen(liczba))//'.stat')
2211       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2212      &  //'.rst'
2213       if(usampl) then
2214           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2215      & liczba(:ilen(liczba))//'.const'
2216       endif
2217
2218       endif
2219 #else
2220       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2221       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2222       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2223       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2224       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2225       if (lentmp.gt.0)
2226      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2227      &    //'.stat')
2228       rest2name=prefix(:ilen(prefix))//'.rst'
2229       if(usampl) then
2230          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2231       endif
2232 #endif
2233 #if defined(AIX) || defined(PGI)
2234       if (me.eq.king .or. .not. out1file)
2235      &   open(iout,file=outname,status='unknown')
2236 #ifdef DEBUG
2237       if (fg_rank.gt.0) then
2238         write (liczba,'(i3.3)') myrank/nfgtasks
2239         write (ll,'(bz,i3.3)') fg_rank
2240         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2241      &   status='unknown')
2242       endif
2243 #endif
2244       if(me.eq.king) then
2245        open(igeom,file=intname,status='unknown',position='append')
2246        open(ipdb,file=pdbname,status='unknown')
2247        open(imol2,file=mol2name,status='unknown')
2248        open(istat,file=statname,status='unknown',position='append')
2249       else
2250 c1out       open(iout,file=outname,status='unknown')
2251       endif
2252 #else
2253       if (me.eq.king .or. .not.out1file)
2254      &    open(iout,file=outname,status='unknown')
2255 #ifdef DEBUG
2256       if (fg_rank.gt.0) then
2257         write (liczba,'(i3.3)') myrank/nfgtasks
2258         write (ll,'(bz,i3.3)') fg_rank
2259         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2260      &   status='unknown')
2261       endif
2262 #endif
2263       if(me.eq.king) then
2264        open(igeom,file=intname,status='unknown',access='append')
2265        open(ipdb,file=pdbname,status='unknown')
2266        open(imol2,file=mol2name,status='unknown')
2267        open(istat,file=statname,status='unknown',access='append')
2268       else
2269 c1out       open(iout,file=outname,status='unknown')
2270       endif
2271 #endif
2272       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2273       csa_seed=prefix(:lenpre)//'.CSA.seed'
2274       csa_history=prefix(:lenpre)//'.CSA.history'
2275       csa_bank=prefix(:lenpre)//'.CSA.bank'
2276       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2277       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2278       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2279 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2280       csa_int=prefix(:lenpre)//'.int'
2281       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2282       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2283       csa_in=prefix(:lenpre)//'.CSA.in'
2284 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2285 C Write file names
2286       if (me.eq.king)then
2287       write (iout,'(80(1h-))')
2288       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2289       write (iout,'(80(1h-))')
2290       write (iout,*) "Input file                      : ",
2291      &  pref_orig(:ilen(pref_orig))//'.inp'
2292       write (iout,*) "Output file                     : ",
2293      &  outname(:ilen(outname))
2294       write (iout,*)
2295       write (iout,*) "Sidechain potential file        : ",
2296      &  sidename(:ilen(sidename))
2297 #ifndef OLDSCP
2298       write (iout,*) "SCp potential file              : ",
2299      &  scpname(:ilen(scpname))
2300 #endif
2301       write (iout,*) "Electrostatic potential file    : ",
2302      &  elename(:ilen(elename))
2303       write (iout,*) "Cumulant coefficient file       : ",
2304      &  fouriername(:ilen(fouriername))
2305       write (iout,*) "Torsional parameter file        : ",
2306      &  torname(:ilen(torname))
2307       write (iout,*) "Double torsional parameter file : ",
2308      &  tordname(:ilen(tordname))
2309       write (iout,*) "SCCOR parameter file : ",
2310      &  sccorname(:ilen(sccorname))
2311       write (iout,*) "Bond & inertia constant file    : ",
2312      &  bondname(:ilen(bondname))
2313       write (iout,*) "Bending parameter file          : ",
2314      &  thetname(:ilen(thetname))
2315       write (iout,*) "Rotamer parameter file          : ",
2316      &  rotname(:ilen(rotname))
2317       write (iout,*) "Thetpdb parameter file          : ",
2318      &  thetname_pdb(:ilen(thetname_pdb))
2319       write (iout,*) "Threading database              : ",
2320      &  patname(:ilen(patname))
2321       if (lentmp.ne.0)
2322      &write (iout,*)" DIRTMP                          : ",
2323      &  tmpdir(:lentmp)
2324       write (iout,'(80(1h-))')
2325       endif
2326       return
2327       end
2328 c----------------------------------------------------------------------------
2329       subroutine card_concat(card)
2330       implicit real*8 (a-h,o-z)
2331       include 'DIMENSIONS'
2332       include 'COMMON.IOUNITS'
2333       character*(*) card
2334       character*80 karta,ucase
2335       external ilen
2336       read (inp,'(a)') karta
2337       karta=ucase(karta)
2338       card=' '
2339       do while (karta(80:80).eq.'&')
2340         card=card(:ilen(card)+1)//karta(:79)
2341         read (inp,'(a)') karta
2342         karta=ucase(karta)
2343       enddo
2344       card=card(:ilen(card)+1)//karta
2345       return
2346       end
2347 c----------------------------------------------------------------------------------
2348       subroutine readrst
2349       implicit real*8 (a-h,o-z)
2350       include 'DIMENSIONS'
2351       include 'COMMON.CHAIN'
2352       include 'COMMON.IOUNITS'
2353       include 'COMMON.MD'
2354       include 'COMMON.CONTROL'
2355       open(irest2,file=rest2name,status='unknown')
2356       read(irest2,*) totT,EK,potE,totE,t_bath
2357       totTafm=totT
2358       do i=1,2*nres
2359          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2360       enddo
2361       do i=1,2*nres
2362          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2363       enddo
2364       if(usampl.or.homol_nset.gt.1) then
2365              read (irest2,*) iset
2366       endif
2367       close(irest2)
2368       return
2369       end
2370 c---------------------------------------------------------------------------------
2371       subroutine read_fragments
2372       implicit real*8 (a-h,o-z)
2373       include 'DIMENSIONS'
2374 #ifdef MPI
2375       include 'mpif.h'
2376 #endif
2377       include 'COMMON.SETUP'
2378       include 'COMMON.CHAIN'
2379       include 'COMMON.IOUNITS'
2380       include 'COMMON.MD'
2381       include 'COMMON.CONTROL'
2382       integer iset1
2383       read(inp,*) nset,nfrag,npair,nfrag_back
2384       if(me.eq.king.or..not.out1file)
2385      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2386      &  " nfrag_back",nfrag_back
2387       do iset1=1,nset
2388          read(inp,*) mset(iset1)
2389        do i=1,nfrag
2390          read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2391      &     qinfrag(i,iset1)
2392          if(me.eq.king.or..not.out1file)
2393      &    write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2394      &     ifrag(2,i,iset1), qinfrag(i,iset1)
2395        enddo
2396        do i=1,npair
2397         read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2398      &    qinpair(i,iset1)
2399         if(me.eq.king.or..not.out1file)
2400      &   write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2401      &    ipair(2,i,iset1), qinpair(i,iset1)
2402        enddo
2403        do i=1,nfrag_back
2404         read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2405      &     wfrag_back(3,i,iset1),
2406      &     ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2407         if(me.eq.king.or..not.out1file)
2408      &   write(iout,*) "A",i,wfrag_back(1,i,iset1),
2409      &   wfrag_back(2,i,iset1),
2410      &   wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2411      &   ifrag_back(2,i,iset1)
2412        enddo
2413       enddo
2414       return
2415       end
2416 C---------------------------------------------------------------------------
2417       subroutine read_afminp
2418             implicit real*8 (a-h,o-z)
2419       include 'DIMENSIONS'
2420 #ifdef MPI
2421       include 'mpif.h'
2422 #endif
2423       include 'COMMON.SETUP'
2424       include 'COMMON.CONTROL'
2425       include 'COMMON.CHAIN'
2426       include 'COMMON.IOUNITS'
2427       include 'COMMON.SBRIDGE'
2428       character*320 afmcard
2429 c      print *, "wchodze"
2430       call card_concat(afmcard)
2431       call readi(afmcard,"BEG",afmbeg,0)
2432       call readi(afmcard,"END",afmend,0)
2433       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2434       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2435       print *,'FORCE=' ,forceAFMconst
2436 CCCC NOW PROPERTIES FOR AFM
2437        distafminit=0.0d0
2438        do i=1,3
2439         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2440        enddo
2441         distafminit=dsqrt(distafminit)
2442         print *,'initdist',distafminit
2443       return
2444       end
2445
2446 c-------------------------------------------------------------------------------
2447       subroutine read_dist_constr
2448       implicit real*8 (a-h,o-z)
2449       include 'DIMENSIONS'
2450 #ifdef MPI
2451       include 'mpif.h'
2452 #endif
2453       include 'COMMON.SETUP'
2454       include 'COMMON.CONTROL'
2455       include 'COMMON.CHAIN'
2456       include 'COMMON.IOUNITS'
2457       include 'COMMON.SBRIDGE'
2458       integer ifrag_(2,100),ipair_(2,100)
2459       double precision wfrag_(100),wpair_(100)
2460       character*500 controlcard
2461 c      write (iout,*) "Calling read_dist_constr"
2462 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2463 c      call flush(iout)
2464       call card_concat(controlcard)
2465       call readi(controlcard,"NFRAG",nfrag_,0)
2466       call readi(controlcard,"NPAIR",npair_,0)
2467       call readi(controlcard,"NDIST",ndist_,0)
2468       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2469       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2470       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2471       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2472       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2473 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2474 c      write (iout,*) "IFRAG"
2475 c      do i=1,nfrag_
2476 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2477 c      enddo
2478 c      write (iout,*) "IPAIR"
2479 c      do i=1,npair_
2480 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2481 c      enddo
2482       if (.not.refstr .and. nfrag.gt.0) then
2483         write (iout,*)
2484      &  "ERROR: no reference structure to compute distance restraints"
2485         write (iout,*)
2486      &  "Restraints must be specified explicitly (NDIST=number)"
2487         stop
2488       endif
2489       if (nfrag.lt.2 .and. npair.gt.0) then
2490         write (iout,*) "ERROR: Less than 2 fragments specified",
2491      &   " but distance restraints between pairs requested"
2492         stop
2493       endif
2494       call flush(iout)
2495       do i=1,nfrag_
2496         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2497         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2498      &    ifrag_(2,i)=nstart_sup+nsup-1
2499 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2500         call flush(iout)
2501         if (wfrag_(i).gt.0.0d0) then
2502         do j=ifrag_(1,i),ifrag_(2,i)-1
2503           do k=j+1,ifrag_(2,i)
2504 c            write (iout,*) "j",j," k",k
2505             ddjk=dist(j,k)
2506             if (constr_dist.eq.1) then
2507             nhpb=nhpb+1
2508             ihpb(nhpb)=j
2509             jhpb(nhpb)=k
2510               dhpb(nhpb)=ddjk
2511             forcon(nhpb)=wfrag_(i)
2512             else if (constr_dist.eq.2) then
2513               if (ddjk.le.dist_cut) then
2514                 nhpb=nhpb+1
2515                 ihpb(nhpb)=j
2516                 jhpb(nhpb)=k
2517                 dhpb(nhpb)=ddjk
2518                 forcon(nhpb)=wfrag_(i)
2519               endif
2520             else
2521               nhpb=nhpb+1
2522               ihpb(nhpb)=j
2523               jhpb(nhpb)=k
2524               dhpb(nhpb)=ddjk
2525               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2526             endif
2527 #ifdef MPI
2528             if (.not.out1file .or. me.eq.king)
2529      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2530      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2531 #else
2532             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2533      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2534 #endif
2535           enddo
2536         enddo
2537         endif
2538       enddo
2539       do i=1,npair_
2540         if (wpair_(i).gt.0.0d0) then
2541         ii = ipair_(1,i)
2542         jj = ipair_(2,i)
2543         if (ii.gt.jj) then
2544           itemp=ii
2545           ii=jj
2546           jj=itemp
2547         endif
2548         do j=ifrag_(1,ii),ifrag_(2,ii)
2549           do k=ifrag_(1,jj),ifrag_(2,jj)
2550             nhpb=nhpb+1
2551             ihpb(nhpb)=j
2552             jhpb(nhpb)=k
2553             forcon(nhpb)=wpair_(i)
2554             dhpb(nhpb)=dist(j,k)
2555 #ifdef MPI
2556             if (.not.out1file .or. me.eq.king)
2557      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2558      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2559 #else
2560             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2561      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2562 #endif
2563           enddo
2564         enddo
2565         endif
2566       enddo
2567       do i=1,ndist_
2568         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2569         if (forcon(nhpb+1).gt.0.0d0) then
2570           nhpb=nhpb+1
2571           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2572 #ifdef MPI
2573           if (.not.out1file .or. me.eq.king)
2574      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2575      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2576 #else
2577           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2578      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2579 #endif
2580         endif
2581       enddo
2582       call flush(iout)
2583       return
2584       end
2585 c-------------------------------------------------------------------------------
2586
2587       subroutine read_constr_homology
2588
2589       include 'DIMENSIONS'
2590 #ifdef MPI
2591       include 'mpif.h'
2592 #endif
2593       include 'COMMON.SETUP'
2594       include 'COMMON.CONTROL'
2595       include 'COMMON.CHAIN'
2596       include 'COMMON.IOUNITS'
2597       include 'COMMON.MD'
2598       include 'COMMON.GEO'
2599       include 'COMMON.INTERACT'
2600       include 'COMMON.NAMES'
2601 c
2602 c For new homol impl
2603 c
2604       include 'COMMON.VAR'
2605 c
2606
2607 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2608 c    &                 dist_cut
2609 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2610 c    &    sigma_odl_temp(maxres,maxres,max_template)
2611       character*2 kic2
2612       character*24 model_ki_dist, model_ki_angle
2613       character*500 controlcard
2614       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2615       logical lprn /.true./
2616       integer ilen
2617       external ilen
2618       logical liiflag
2619 c
2620 c     FP - Nov. 2014 Temporary specifications for new vars
2621 c
2622       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2623      &                 rescore3_tmp
2624       double precision, dimension (max_template,maxres) :: rescore
2625       double precision, dimension (max_template,maxres) :: rescore2
2626       double precision, dimension (max_template,maxres) :: rescore3
2627       character*24 pdbfile,tpl_k_rescore
2628 c -----------------------------------------------------------------
2629 c Reading multiple PDB ref structures and calculation of retraints
2630 c not using pre-computed ones stored in files model_ki_{dist,angle}
2631 c FP (Nov., 2014)
2632 c -----------------------------------------------------------------
2633 c
2634 c
2635 c Alternative: reading from input
2636       call card_concat(controlcard)
2637       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2638       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2639       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2640       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2641       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2642       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2643       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2644       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2645       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2646       if(.not.read2sigma.and.start_from_model) then
2647           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
2648      &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2649           start_from_model=.false.
2650       endif
2651       if(start_from_model .and. (me.eq.king .or. .not. out1file))
2652      &    write(iout,*) 'START_FROM_MODELS is ON'
2653       if(start_from_model .and. rest) then
2654         if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2655          write(iout,*) 'START_FROM_MODELS is OFF'
2656          write(iout,*) 'remove restart keyword from input'
2657         endif
2658       endif
2659       if (homol_nset.gt.1)then
2660          call card_concat(controlcard)
2661          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2662          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2663           write(iout,*) "iset homology_weight "
2664           do i=1,homol_nset
2665            write(iout,*) i,waga_homology(i)
2666           enddo
2667          endif
2668          iset=mod(kolor,homol_nset)+1
2669       else
2670        iset=1
2671        waga_homology(1)=1.0
2672       endif
2673
2674 cd      write (iout,*) "nnt",nnt," nct",nct
2675 cd      call flush(iout)
2676
2677
2678       lim_odl=0
2679       lim_dih=0
2680 c
2681 c      write(iout,*) 'nnt=',nnt,'nct=',nct
2682 c
2683       do i = nnt,nct
2684         do k=1,constr_homology
2685          idomain(k,i)=0
2686         enddo
2687       enddo
2688
2689       ii=0
2690       do i = nnt,nct-2
2691         do j=i+2,nct
2692         ii=ii+1
2693         ii_in_use(ii)=0
2694         enddo
2695       enddo
2696
2697       do k=1,constr_homology
2698
2699         read(inp,'(a)') pdbfile
2700         if(me.eq.king .or. .not. out1file)
2701      &    write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2702      &  pdbfile(:ilen(pdbfile))
2703         open(ipdbin,file=pdbfile,status='old',err=33)
2704         goto 34
2705   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
2706      &  pdbfile(:ilen(pdbfile))
2707         stop
2708   34    continue
2709 c        print *,'Begin reading pdb data'
2710 c
2711 c Files containing res sim or local scores (former containing sigmas)
2712 c
2713
2714         write(kic2,'(bz,i2.2)') k
2715
2716         tpl_k_rescore="template"//kic2//".sco"
2717
2718         unres_pdb=.false.
2719         if (read2sigma) then
2720           call readpdb_template(k)
2721         else
2722           call readpdb
2723         endif
2724 c
2725 c     Distance restraints
2726 c
2727 c          ... --> odl(k,ii)
2728 C Copy the coordinates from reference coordinates (?)
2729         do i=1,2*nres
2730           do j=1,3
2731             c(j,i)=cref(j,i,1)
2732 c           write (iout,*) "c(",j,i,") =",c(j,i)
2733           enddo
2734         enddo
2735 c
2736 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2737 c
2738 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2739           open (ientin,file=tpl_k_rescore,status='old')
2740           if (nnt.gt.1) rescore(k,1)=0.0d0
2741           do irec=nnt,nct ! loop for reading res sim
2742             if (read2sigma) then
2743              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2744      &                                rescore3_tmp,idomain_tmp
2745              i_tmp=i_tmp+nnt-1
2746              idomain(k,i_tmp)=idomain_tmp
2747              rescore(k,i_tmp)=rescore_tmp
2748              rescore2(k,i_tmp)=rescore2_tmp
2749              rescore3(k,i_tmp)=rescore3_tmp
2750              if (.not. out1file .or. me.eq.king)
2751      &        write(iout,'(a7,i5,3f10.5,i5)') "rescore",
2752      &                      i_tmp,rescore2_tmp,rescore_tmp,
2753      &                                rescore3_tmp,idomain_tmp
2754             else
2755              idomain(k,irec)=1
2756              read (ientin,*,end=1401) rescore_tmp
2757
2758 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2759              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2760 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2761             endif
2762           enddo
2763  1401   continue
2764         close (ientin)
2765         if (waga_dist.ne.0.0d0) then
2766           ii=0
2767           do i = nnt,nct-2
2768             do j=i+2,nct
2769
2770               x12=c(1,i)-c(1,j)
2771               y12=c(2,i)-c(2,j)
2772               z12=c(3,i)-c(3,j)
2773               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2774 c              write (iout,*) k,i,j,distal,dist2_cut
2775
2776             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2777      &            .and. distal.le.dist2_cut ) then
2778
2779               ii=ii+1
2780               ii_in_use(ii)=1
2781               l_homo(k,ii)=.true.
2782
2783 c             write (iout,*) "k",k
2784 c             write (iout,*) "i",i," j",j," constr_homology",
2785 c    &                       constr_homology
2786               ires_homo(ii)=i
2787               jres_homo(ii)=j
2788               odl(k,ii)=distal
2789               if (read2sigma) then
2790                 sigma_odl(k,ii)=0
2791                 do ik=i,j
2792                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2793                 enddo
2794                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2795                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2796      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2797               else
2798                 if (odl(k,ii).le.dist_cut) then
2799                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2800                 else
2801 #ifdef OLDSIGMA
2802                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2803      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2804 #else
2805                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2806      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2807 #endif
2808                 endif
2809               endif
2810               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2811             else
2812               ii=ii+1
2813               l_homo(k,ii)=.false.
2814             endif
2815             enddo
2816           enddo
2817         lim_odl=ii
2818         endif
2819 c
2820 c     Theta, dihedral and SC retraints
2821 c
2822         if (waga_angle.gt.0.0d0) then
2823 c         open (ientin,file=tpl_k_sigma_dih,status='old')
2824 c         do irec=1,maxres-3 ! loop for reading sigma_dih
2825 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2826 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2827 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2828 c    &                            sigma_dih(k,i+nnt-1)
2829 c         enddo
2830 c1402   continue
2831 c         close (ientin)
2832           do i = nnt+3,nct
2833             if (idomain(k,i).eq.0) then
2834                sigma_dih(k,i)=0.0
2835                cycle
2836             endif
2837             dih(k,i)=phiref(i) ! right?
2838 c           read (ientin,*) sigma_dih(k,i) ! original variant
2839 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
2840 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2841 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2842 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
2843 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
2844
2845             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2846      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
2847 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2848 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
2849 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2850 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2851 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
2852 c   Instead of res sim other local measure of b/b str reliability possible
2853             if (sigma_dih(k,i).ne.0)
2854      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2855 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2856           enddo
2857           lim_dih=nct-nnt-2
2858         endif
2859
2860         if (waga_theta.gt.0.0d0) then
2861 c         open (ientin,file=tpl_k_sigma_theta,status='old')
2862 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2863 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2864 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2865 c    &                              sigma_theta(k,i+nnt-1)
2866 c         enddo
2867 c1403   continue
2868 c         close (ientin)
2869
2870           do i = nnt+2,nct ! right? without parallel.
2871 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
2872 c         do i=ithet_start,ithet_end ! with FG parallel.
2873              if (idomain(k,i).eq.0) then
2874               sigma_theta(k,i)=0.0
2875               cycle
2876              endif
2877              thetatpl(k,i)=thetaref(i)
2878 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2879 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
2880 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2881 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
2882 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
2883              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
2884      &                        rescore(k,i-2))/3.0
2885 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2886              if (sigma_theta(k,i).ne.0)
2887      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2888
2889 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2890 c                             rescore(k,i-2) !  right expression ?
2891 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2892           enddo
2893         endif
2894
2895         if (waga_d.gt.0.0d0) then
2896 c       open (ientin,file=tpl_k_sigma_d,status='old')
2897 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2898 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2899 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2900 c    &                          sigma_d(k,i+nnt-1)
2901 c         enddo
2902 c1404   continue
2903
2904           do i = nnt,nct ! right? without parallel.
2905 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
2906 c         do i=loc_start,loc_end ! with FG parallel.
2907                if (itype(i).eq.10) cycle
2908                if (idomain(k,i).eq.0 ) then
2909                   sigma_d(k,i)=0.0
2910                   cycle
2911                endif
2912                xxtpl(k,i)=xxref(i)
2913                yytpl(k,i)=yyref(i)
2914                zztpl(k,i)=zzref(i)
2915 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2916 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2917 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2918 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
2919                sigma_d(k,i)=rescore3(k,i) !  right expression ?
2920                if (sigma_d(k,i).ne.0)
2921      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2922
2923 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
2924 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2925 c              read (ientin,*) sigma_d(k,i) ! 1st variant
2926           enddo
2927         endif
2928       enddo
2929 c
2930 c remove distance restraints not used in any model from the list
2931 c shift data in all arrays
2932 c
2933       if (waga_dist.ne.0.0d0) then
2934         ii=0
2935         liiflag=.true.
2936         do i=nnt,nct-2
2937          do j=i+2,nct
2938           ii=ii+1
2939           if (ii_in_use(ii).eq.0.and.liiflag) then
2940             liiflag=.false.
2941             iistart=ii
2942           endif
2943           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
2944      &                   .not.liiflag.and.ii.eq.lim_odl) then
2945              if (ii.eq.lim_odl) then
2946               iishift=ii-iistart+1
2947              else
2948               iishift=ii-iistart
2949              endif
2950              liiflag=.true.
2951              do ki=iistart,lim_odl-iishift
2952               ires_homo(ki)=ires_homo(ki+iishift)
2953               jres_homo(ki)=jres_homo(ki+iishift)
2954               ii_in_use(ki)=ii_in_use(ki+iishift)
2955               do k=1,constr_homology
2956                odl(k,ki)=odl(k,ki+iishift)
2957                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
2958                l_homo(k,ki)=l_homo(k,ki+iishift)
2959               enddo
2960              enddo
2961              ii=ii-iishift
2962              lim_odl=lim_odl-iishift
2963           endif
2964          enddo
2965         enddo
2966       endif
2967       if (constr_homology.gt.0) call homology_partition
2968       if (constr_homology.gt.0) call init_int_table
2969 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2970 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2971 c
2972 c Print restraints
2973 c
2974       if (.not.lprn) return
2975 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2976       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2977        write (iout,*) "Distance restraints from templates"
2978        do ii=1,lim_odl
2979        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
2980      &  ii,ires_homo(ii),jres_homo(ii),
2981      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
2982      &  ki=1,constr_homology)
2983        enddo
2984        write (iout,*) "Dihedral angle restraints from templates"
2985        do i=nnt+3,nct
2986         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
2987      &      (rad2deg*dih(ki,i),
2988      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2989        enddo
2990        write (iout,*) "Virtual-bond angle restraints from templates"
2991        do i=nnt+2,nct
2992         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
2993      &      (rad2deg*thetatpl(ki,i),
2994      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2995        enddo
2996        write (iout,*) "SC restraints from templates"
2997        do i=nnt,nct
2998         write(iout,'(i5,100(4f8.2,4x))') i,
2999      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3000      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3001        enddo
3002       endif
3003 c -----------------------------------------------------------------
3004       return
3005       end
3006 c----------------------------------------------------------------------
3007
3008 #ifdef WINIFL
3009       subroutine flush(iu)
3010       return
3011       end
3012 #endif
3013 #ifdef AIX
3014       subroutine flush(iu)
3015       call flush_(iu)
3016       return
3017       end
3018 #endif
3019 c------------------------------------------------------------------------------
3020       subroutine copy_to_tmp(source)
3021       include "DIMENSIONS"
3022       include "COMMON.IOUNITS"
3023       character*(*) source
3024       character* 256 tmpfile
3025       integer ilen
3026       external ilen
3027       logical ex
3028       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3029       inquire(file=tmpfile,exist=ex)
3030       if (ex) then
3031         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3032      &   " to temporary directory..."
3033         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3034         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3035       endif
3036       return
3037       end
3038 c------------------------------------------------------------------------------
3039       subroutine move_from_tmp(source)
3040       include "DIMENSIONS"
3041       include "COMMON.IOUNITS"
3042       character*(*) source
3043       integer ilen
3044       external ilen
3045       write (*,*) "Moving ",source(:ilen(source)),
3046      & " from temporary directory to working directory"
3047       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3048       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3049       return
3050       end
3051 c------------------------------------------------------------------------------
3052       subroutine random_init(seed)
3053 C
3054 C Initialize random number generator
3055 C
3056       implicit real*8 (a-h,o-z)
3057       include 'DIMENSIONS'
3058 #ifdef AMD64
3059       integer*8 iseedi8
3060 #endif
3061 #ifdef MPI
3062       include 'mpif.h'
3063       logical OKRandom, prng_restart
3064       real*8  r1
3065       integer iseed_array(4)
3066 #endif
3067       include 'COMMON.IOUNITS'
3068       include 'COMMON.TIME1'
3069       include 'COMMON.THREAD'
3070       include 'COMMON.SBRIDGE'
3071       include 'COMMON.CONTROL'
3072       include 'COMMON.MCM'
3073       include 'COMMON.MAP'
3074       include 'COMMON.HEADER'
3075       include 'COMMON.CSA'
3076       include 'COMMON.CHAIN'
3077       include 'COMMON.MUCA'
3078       include 'COMMON.MD'
3079       include 'COMMON.FFIELD'
3080       include 'COMMON.SETUP'
3081       iseed=-dint(dabs(seed))
3082       if (iseed.eq.0) then
3083         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3084      &    'Random seed undefined. The program will stop.'
3085         write (*,'(/80(1h*)/20x,a/80(1h*))')
3086      &    'Random seed undefined. The program will stop.'
3087 #ifdef MPI
3088         call mpi_finalize(mpi_comm_world,ierr)
3089 #endif
3090         stop 'Bad random seed.'
3091       endif
3092 #ifdef MPI
3093       if (fg_rank.eq.0) then
3094       seed=seed*(me+1)+1
3095 #ifdef AMD64
3096       iseedi8=dint(seed)
3097       if(me.eq.king .or. .not. out1file)
3098      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3099       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3100       OKRandom = prng_restart(me,iseedi8)
3101 #else
3102       do i=1,4
3103        tmp=65536.0d0**(4-i)
3104        iseed_array(i) = dint(seed/tmp)
3105        seed=seed-iseed_array(i)*tmp
3106       enddo
3107       if(me.eq.king)
3108      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3109      &                 (iseed_array(i),i=1,4)
3110       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3111      &                 (iseed_array(i),i=1,4)
3112       OKRandom = prng_restart(me,iseed_array)
3113 #endif
3114       if (OKRandom) then
3115 c        r1 = prng_next(me)
3116         r1=ran_number(0.0D0,1.0D0)
3117         if(me.eq.king .or. .not. out1file)
3118      &   write (iout,*) 'ran_num',r1
3119         if (r1.lt.0.0d0) OKRandom=.false.
3120       endif
3121       if (.not.OKRandom) then
3122         write (iout,*) 'PRNG IS NOT WORKING!!!'
3123         print *,'PRNG IS NOT WORKING!!!'
3124         if (me.eq.0) then
3125          call flush(iout)
3126          call mpi_abort(mpi_comm_world,error_msg,ierr)
3127          stop
3128         else
3129          write (iout,*) 'too many processors for parallel prng'
3130          write (*,*) 'too many processors for parallel prng'
3131          call flush(iout)
3132          stop
3133         endif
3134       endif
3135       endif
3136 #else
3137       call vrndst(iseed)
3138       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3139 #endif
3140       return
3141       end