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