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