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