source/unres/src_MD-M/readrtns_CSA.F

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