source/unres/src_MD-M/readrtns_CSA.F

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