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