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*320 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,20
 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       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 684       do i=1,maxres
 685         dyn_ss_mask(i)=.false.
 686       enddo
 687       do i=1,maxres-1
 688         do j=i+1,maxres
 689           dyn_ssbond_ij(i,j)=1.0d300
 690         enddo
 691       enddo
 692       call reada(weightcard,"HT",Ht,0.0D0)
 693       if (dyn_ss) then
 694         ss_depth=ebr/wsc-0.25*eps(1,1)
 695         Ht=Ht/wsc-0.25*eps(1,1)
 696         akcm=akcm*wstrain/wsc
 697         akth=akth*wstrain/wsc
 698         akct=akct*wstrain/wsc
 699         v1ss=v1ss*wstrain/wsc
 700         v2ss=v2ss*wstrain/wsc
 701         v3ss=v3ss*wstrain/wsc
 702       else
 703         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 704       endif
 705
 706       if(me.eq.king.or..not.out1file) then
 707        write (iout,*) "Parameters of the SS-bond potential:"
 708        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 709      & " AKCT",akct
 710        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 711        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 712        write (iout,*)" HT",Ht
 713        print *,'indpdb=',indpdb,' pdbref=',pdbref
 714       endif
 715       if (indpdb.gt.0 .or. pdbref) then
 716         read(inp,'(a)') pdbfile
 717         if(me.eq.king.or..not.out1file)
 718      &   write (iout,'(2a)') 'PDB data will be read from file ',
 719      &   pdbfile(:ilen(pdbfile))
 720         open(ipdbin,file=pdbfile,status='old',err=33)
 721         goto 34
 722   33    write (iout,'(a)') 'Error opening PDB file.'
 723         stop
 724   34    continue
 725 c        write (iout,*) 'Begin reading pdb data'
 726 c        call flush(iout)
 727         call readpdb
 728 c        write (iout,*) 'Finished reading pdb data'
 729 c        call flush(iout)
 730         if(me.eq.king.or..not.out1file)
 731      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 732      &   ' nstart_sup=',nstart_sup
 733         do i=1,nres
 734           itype_pdb(i)=itype(i)
 735         enddo
 736         close (ipdbin)
 737         nnt=nstart_sup
 738         nct=nstart_sup+nsup-1
 739         call contact(.false.,ncont_ref,icont_ref,co)
 740
 741         if (sideadd) then
 742          if(me.eq.king.or..not.out1file)
 743      &    write(iout,*)'Adding sidechains'
 744          maxsi=1000
 745          do i=2,nres-1
 746           iti=itype(i)
 747           if (iti.ne.10 .and. itype(i).ne.21) then
 748             nsi=0
 749             fail=.true.
 750             do while (fail.and.nsi.le.maxsi)
 751               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 752               nsi=nsi+1
 753             enddo
 754             if(fail) write(iout,*)'Adding sidechain failed for res ',
 755      &              i,' after ',nsi,' trials'
 756           endif
 757          enddo
 758         endif
 759       endif
 760       if (indpdb.eq.0) then
 761 C Read sequence if not taken from the pdb file.
 762         read (inp,*) nres
 763 c        print *,'nres=',nres
 764         if (iscode.gt.0) then
 765           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 766         else
 767           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 768         endif
 769 C Convert sequence to numeric code
 770         do i=1,nres
 771           itype(i)=rescode(i,sequence(i),iscode)
 772         enddo
 773 C Assign initial virtual bond lengths
 774         do i=2,nres
 775           vbld(i)=vbl
 776           vbld_inv(i)=vblinv
 777         enddo
 778         do i=2,nres-1
 779           vbld(i+nres)=dsc(itype(i))
 780           vbld_inv(i+nres)=dsc_inv(itype(i))
 781 c          write (iout,*) "i",i," itype",itype(i),
 782 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 783         enddo
 784       endif
 785 c      print *,nres
 786 c      print '(20i4)',(itype(i),i=1,nres)
 787       do i=1,nres
 788 #ifdef PROCOR
 789         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 790 #else
 791         if (itype(i).eq.21) then
 792 #endif
 793           itel(i)=0
 794 #ifdef PROCOR
 795         else if (itype(i+1).ne.20) then
 796 #else
 797         else if (itype(i).ne.20) then
 798 #endif
 799           itel(i)=1
 800         else
 801           itel(i)=2
 802         endif
 803       enddo
 804       if(me.eq.king.or..not.out1file)then
 805        write (iout,*) "ITEL"
 806        do i=1,nres-1
 807          write (iout,*) i,itype(i),itel(i)
 808        enddo
 809        print *,'Call Read_Bridge.'
 810       endif
 811       call read_bridge
 812 C 8/13/98 Set limits to generating the dihedral angles
 813       do i=1,nres
 814         phibound(1,i)=-pi
 815         phibound(2,i)=pi
 816       enddo
 817       read (inp,*) ndih_constr
 818       if (ndih_constr.gt.0) then
 819         read (inp,*) ftors
 820         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 821         if(me.eq.king.or..not.out1file)then
 822          write (iout,*)
 823      &   'There are',ndih_constr,' constraints on phi angles.'
 824          do i=1,ndih_constr
 825           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 826          enddo
 827         endif
 828         do i=1,ndih_constr
 829           phi0(i)=deg2rad*phi0(i)
 830           drange(i)=deg2rad*drange(i)
 831         enddo
 832         if(me.eq.king.or..not.out1file)
 833      &   write (iout,*) 'FTORS',ftors
 834         do i=1,ndih_constr
 835           ii = idih_constr(i)
 836           phibound(1,ii) = phi0(i)-drange(i)
 837           phibound(2,ii) = phi0(i)+drange(i)
 838         enddo
 839       endif
 840       nnt=1
 841 #ifdef MPI
 842       if (me.eq.king) then
 843 #endif
 844        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 845        do i=1,nres
 846          write (iout,'(a3,i5,2f10.1)')
 847      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 848        enddo
 849 #ifdef MP
 850       endif
 851 #endif
 852       nct=nres
 853 cd      print *,'NNT=',NNT,' NCT=',NCT
 854       if (itype(1).eq.21) nnt=2
 855       if (itype(nres).eq.21) nct=nct-1
 856       if (pdbref) then
 857         if(me.eq.king.or..not.out1file)
 858      &   write (iout,'(a,i3)') 'nsup=',nsup
 859         nstart_seq=nnt
 860         if (nsup.le.(nct-nnt+1)) then
 861           do i=0,nct-nnt+1-nsup
 862             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 863               nstart_seq=nnt+i
 864               goto 111
 865             endif
 866           enddo
 867           write (iout,'(a)')
 868      &            'Error - sequences to be superposed do not match.'
 869           stop
 870         else
 871           do i=0,nsup-(nct-nnt+1)
 872             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 873      &      then
 874               nstart_sup=nstart_sup+i
 875               nsup=nct-nnt+1
 876               goto 111
 877             endif
 878           enddo
 879           write (iout,'(a)')
 880      &            'Error - sequences to be superposed do not match.'
 881         endif
 882   111   continue
 883         if (nsup.eq.0) nsup=nct-nnt
 884         if (nstart_sup.eq.0) nstart_sup=nnt
 885         if (nstart_seq.eq.0) nstart_seq=nnt
 886         if(me.eq.king.or..not.out1file)
 887      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 888      &                 ' nstart_seq=',nstart_seq
 889       endif
 890 c--- Zscore rms -------
 891       if (nz_start.eq.0) nz_start=nnt
 892       if (nz_end.eq.0 .and. nsup.gt.0) then
 893         nz_end=nnt+nsup-1
 894       else if (nz_end.eq.0) then
 895         nz_end=nct
 896       endif
 897       if(me.eq.king.or..not.out1file)then
 898        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 899        write (iout,*) 'IZ_SC=',iz_sc
 900       endif
 901 c----------------------
 902       call init_int_table
 903       if (refstr) then
 904         if (.not.pdbref) then
 905           call read_angles(inp,*38)
 906           goto 39
 907    38     write (iout,'(a)') 'Error reading reference structure.'
 908 #ifdef MPI
 909           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 910           stop 'Error reading reference structure'
 911 #endif
 912    39     call chainbuild
 913           call setup_var
 914 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 915           nstart_sup=nnt
 916           nstart_seq=nnt
 917           nsup=nct-nnt+1
 918           kkk=1
 919           do i=1,2*nres
 920             do j=1,3
 921               cref(j,i,kkk)=c(j,i)
 922             enddo
 923           enddo
 924           call contact(.true.,ncont_ref,icont_ref,co)
 925         endif
 926 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 927         call flush(iout)
 928         if (constr_dist.gt.0) call read_dist_constr
 929         write (iout,*) "After read_dist_constr nhpb",nhpb
 930         call hpb_partition
 931         if(me.eq.king.or..not.out1file)
 932      &   write (iout,*) 'Contact order:',co
 933         if (pdbref) then
 934         if(me.eq.king.or..not.out1file)
 935      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 936         do i=1,ncont_ref
 937           do j=1,2
 938             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 939           enddo
 940           if(me.eq.king.or..not.out1file)
 941      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 942      &     icont_ref(1,i),' ',
 943      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 944         enddo
 945         endif
 946       endif
 947       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 948      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 949      &    modecalc.ne.10) then
 950 C If input structure hasn't been supplied from the PDB file read or generate
 951 C initial geometry.
 952         if (iranconf.eq.0 .and. .not. extconf) then
 953           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 954      &     write (iout,'(a)') 'Initial geometry will be read in.'
 955           if (read_cart) then
 956             read(inp,'(8f10.5)',end=36,err=36)
 957      &       ((c(l,k),l=1,3),k=1,nres),
 958      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 959             write (iout,*) "Exit READ_CART"
 960             write (iout,'(8f10.5)')
 961      &       ((c(l,k),l=1,3),k=1,nres),
 962      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 963             call int_from_cart1(.true.)
 964             write (iout,*) "Finish INT_TO_CART"
 965             do i=1,nres-1
 966               do j=1,3
 967                 dc(j,i)=c(j,i+1)-c(j,i)
 968                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 969               enddo
 970             enddo
 971             do i=nnt,nct
 972               if (itype(i).ne.10 .and. itype(i).ne.21) then
 973                 do j=1,3
 974                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 975                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 976                 enddo
 977               endif
 978             enddo
 979             return
 980           else
 981             call read_angles(inp,*36)
 982           endif
 983           goto 37
 984    36     write (iout,'(a)') 'Error reading angle file.'
 985 #ifdef MPI
 986           call mpi_finalize( MPI_COMM_WORLD,IERR )
 987 #endif
 988           stop 'Error reading angle file.'
 989    37     continue
 990         else if (extconf) then
 991          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 992      &    write (iout,'(a)') 'Extended chain initial geometry.'
 993          do i=3,nres
 994           theta(i)=90d0*deg2rad
 995          enddo
 996          do i=4,nres
 997           phi(i)=180d0*deg2rad
 998          enddo
 999          do i=2,nres-1
1000           alph(i)=110d0*deg2rad
1001          enddo
1002          do i=2,nres-1
1003           omeg(i)=-120d0*deg2rad
1004          enddo
1005         else
1006           if(me.eq.king.or..not.out1file)
1007      &     write (iout,'(a)') 'Random-generated initial geometry.'
1008
1009
1010 #ifdef MPI
1011           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1012      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1013 #endif
1014             do itrial=1,100
1015               itmp=1
1016               call gen_rand_conf(itmp,*30)
1017               goto 40
1018    30         write (iout,*) 'Failed to generate random conformation',
1019      &          ', itrial=',itrial
1020               write (*,*) 'Processor:',me,
1021      &          ' Failed to generate random conformation',
1022      &          ' itrial=',itrial
1023               call intout
1024
1025 #ifdef AIX
1026               call flush_(iout)
1027 #else
1028               call flush(iout)
1029 #endif
1030             enddo
1031             write (iout,'(a,i3,a)') 'Processor:',me,
1032      &        ' error in generating random conformation.'
1033             write (*,'(a,i3,a)') 'Processor:',me,
1034      &        ' error in generating random conformation.'
1035             call flush(iout)
1036 #ifdef MPI
1037             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1038    40       continue
1039           endif
1040 #else
1041           write (*,'(a)')
1042      &      ' error in generating random conformation.'
1043           stop
1044    40     continue
1045 #endif
1046         endif
1047       elseif (modecalc.eq.4) then
1048         read (inp,'(a)') intinname
1049         open (intin,file=intinname,status='old',err=333)
1050         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1051      &  write (iout,'(a)') 'intinname',intinname
1052         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1053         goto 334
1054   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1055 #ifdef MPI
1056         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1057 #endif
1058         stop 'Error opening angle file.'
1059   334   continue
1060
1061       endif
1062 C Generate distance constraints, if the PDB structure is to be regularized.
1063       if (nthread.gt.0) then
1064         call read_threadbase
1065       endif
1066       call setup_var
1067       if (me.eq.king .or. .not. out1file)
1068      & call intout
1069       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1070         write (iout,'(/a,i3,a)')
1071      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1072         write (iout,'(20i4)') (iss(i),i=1,ns)
1073        if (dyn_ss) then
1074           write(iout,*)"Running with dynamic disulfide-bond formation"
1075        else
1076         write (iout,'(/a/)') 'Pre-formed links are:'
1077         do i=1,nss
1078           i1=ihpb(i)-nres
1079           i2=jhpb(i)-nres
1080           it1=itype(i1)
1081           it2=itype(i2)
1082           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1083      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1084      &    ebr,forcon(i)
1085         enddo
1086         write (iout,'(a)')
1087        endif
1088       endif
1089       if (ns.gt.0.and.dyn_ss) then
1090           do i=nss+1,nhpb
1091             ihpb(i-nss)=ihpb(i)
1092             jhpb(i-nss)=jhpb(i)
1093             forcon(i-nss)=forcon(i)
1094             dhpb(i-nss)=dhpb(i)
1095           enddo
1096           nhpb=nhpb-nss
1097           nss=0
1098           call hpb_partition
1099           do i=1,ns
1100             dyn_ss_mask(iss(i))=.true.
1101           enddo
1102       endif
1103       if (i2ndstr.gt.0) call secstrp2dihc
1104 c      call geom_to_var(nvar,x)
1105 c      call etotal(energia(0))
1106 c      call enerprint(energia(0))
1107 c      call briefout(0,etot)
1108 c      stop
1109 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1110 cd    write (iout,'(a)') 'Variable list:'
1111 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1112 #ifdef MPI
1113       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1114      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1115      &  'Processor',myrank,': end reading molecular data.'
1116 #endif
1117       return
1118       end
1119 c--------------------------------------------------------------------------
1120       logical function seq_comp(itypea,itypeb,length)
1121       implicit none
1122       integer length,itypea(length),itypeb(length)
1123       integer i
1124       do i=1,length
1125         if (itypea(i).ne.itypeb(i)) then
1126           seq_comp=.false.
1127           return
1128         endif
1129       enddo
1130       seq_comp=.true.
1131       return
1132       end
1133 c-----------------------------------------------------------------------------
1134       subroutine read_bridge
1135 C Read information about disulfide bridges.
1136       implicit real*8 (a-h,o-z)
1137       include 'DIMENSIONS'
1138 #ifdef MPI
1139       include 'mpif.h'
1140 #endif
1141       include 'COMMON.IOUNITS'
1142       include 'COMMON.GEO'
1143       include 'COMMON.VAR'
1144       include 'COMMON.INTERACT'
1145       include 'COMMON.LOCAL'
1146       include 'COMMON.NAMES'
1147       include 'COMMON.CHAIN'
1148       include 'COMMON.FFIELD'
1149       include 'COMMON.SBRIDGE'
1150       include 'COMMON.HEADER'
1151       include 'COMMON.CONTROL'
1152       include 'COMMON.DBASE'
1153       include 'COMMON.THREAD'
1154       include 'COMMON.TIME1'
1155       include 'COMMON.SETUP'
1156 C Read bridging residues.
1157       read (inp,*) ns,(iss(i),i=1,ns)
1158       print *,'ns=',ns
1159       if(me.eq.king.or..not.out1file)
1160      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1161 C Check whether the specified bridging residues are cystines.
1162       do i=1,ns
1163         if (itype(iss(i)).ne.1) then
1164           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1165      &   'Do you REALLY think that the residue ',
1166      &    restyp(itype(iss(i))),i,
1167      &   ' can form a disulfide bridge?!!!'
1168           write (*,'(2a,i3,a)')
1169      &   'Do you REALLY think that the residue ',
1170      &    restyp(itype(iss(i))),i,
1171      &   ' can form a disulfide bridge?!!!'
1172 #ifdef MPI
1173          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1174          stop
1175 #endif
1176         endif
1177       enddo
1178 C Read preformed bridges.
1179       if (ns.gt.0) then
1180       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1181       if(fg_rank.eq.0)
1182      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1183       if (nss.gt.0) then
1184         nhpb=nss
1185 C Check if the residues involved in bridges are in the specified list of
1186 C bridging residues.
1187         do i=1,nss
1188           do j=1,i-1
1189             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1190      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1191               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1192      &      ' contains residues present in other pairs.'
1193               write (*,'(a,i3,a)') 'Disulfide pair',i,
1194      &      ' contains residues present in other pairs.'
1195 #ifdef MPI
1196               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1197               stop
1198 #endif
1199             endif
1200           enddo
1201           do j=1,ns
1202             if (ihpb(i).eq.iss(j)) goto 10
1203           enddo
1204           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1205    10     continue
1206           do j=1,ns
1207             if (jhpb(i).eq.iss(j)) goto 20
1208           enddo
1209           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1210    20     continue
1211           dhpb(i)=dbr
1212           forcon(i)=fbr
1213         enddo
1214         do i=1,nss
1215           ihpb(i)=ihpb(i)+nres
1216           jhpb(i)=jhpb(i)+nres
1217         enddo
1218       endif
1219       endif
1220       return
1221       end
1222 c----------------------------------------------------------------------------
1223       subroutine read_x(kanal,*)
1224       implicit real*8 (a-h,o-z)
1225       include 'DIMENSIONS'
1226       include 'COMMON.GEO'
1227       include 'COMMON.VAR'
1228       include 'COMMON.CHAIN'
1229       include 'COMMON.IOUNITS'
1230       include 'COMMON.CONTROL'
1231       include 'COMMON.LOCAL'
1232       include 'COMMON.INTERACT'
1233 c Read coordinates from input
1234 c
1235       read(kanal,'(8f10.5)',end=10,err=10)
1236      &  ((c(l,k),l=1,3),k=1,nres),
1237      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1238       do j=1,3
1239         c(j,nres+1)=c(j,1)
1240         c(j,2*nres)=c(j,nres)
1241       enddo
1242       call int_from_cart1(.false.)
1243       do i=1,nres-1
1244         do j=1,3
1245           dc(j,i)=c(j,i+1)-c(j,i)
1246           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1247         enddo
1248       enddo
1249       do i=nnt,nct
1250         if (itype(i).ne.10 .and. itype(i).ne.21) then
1251           do j=1,3
1252             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1253             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1254           enddo
1255         endif
1256       enddo
1257
1258       return
1259    10 return1
1260       end
1261 c----------------------------------------------------------------------------
1262       subroutine read_threadbase
1263       implicit real*8 (a-h,o-z)
1264       include 'DIMENSIONS'
1265       include 'COMMON.IOUNITS'
1266       include 'COMMON.GEO'
1267       include 'COMMON.VAR'
1268       include 'COMMON.INTERACT'
1269       include 'COMMON.LOCAL'
1270       include 'COMMON.NAMES'
1271       include 'COMMON.CHAIN'
1272       include 'COMMON.FFIELD'
1273       include 'COMMON.SBRIDGE'
1274       include 'COMMON.HEADER'
1275       include 'COMMON.CONTROL'
1276       include 'COMMON.DBASE'
1277       include 'COMMON.THREAD'
1278       include 'COMMON.TIME1'
1279 C Read pattern database for threading.
1280       read (icbase,*) nseq
1281       do i=1,nseq
1282         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1283      &   nres_base(2,i),nres_base(3,i)
1284         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1285      &   nres_base(1,i))
1286 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1287 c    &   nres_base(2,i),nres_base(3,i)
1288 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1289 c    &   nres_base(1,i))
1290       enddo
1291       close (icbase)
1292       if (weidis.eq.0.0D0) weidis=0.1D0
1293       do i=nnt,nct
1294         do j=i+2,nct
1295           nhpb=nhpb+1
1296           ihpb(nhpb)=i
1297           jhpb(nhpb)=j
1298           forcon(nhpb)=weidis
1299         enddo
1300       enddo
1301       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1302       write (iout,'(a,i5)') 'nexcl: ',nexcl
1303       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1304       return
1305       end
1306 c------------------------------------------------------------------------------
1307       subroutine setup_var
1308       implicit real*8 (a-h,o-z)
1309       include 'DIMENSIONS'
1310       include 'COMMON.IOUNITS'
1311       include 'COMMON.GEO'
1312       include 'COMMON.VAR'
1313       include 'COMMON.INTERACT'
1314       include 'COMMON.LOCAL'
1315       include 'COMMON.NAMES'
1316       include 'COMMON.CHAIN'
1317       include 'COMMON.FFIELD'
1318       include 'COMMON.SBRIDGE'
1319       include 'COMMON.HEADER'
1320       include 'COMMON.CONTROL'
1321       include 'COMMON.DBASE'
1322       include 'COMMON.THREAD'
1323       include 'COMMON.TIME1'
1324 C Set up variable list.
1325       ntheta=nres-2
1326       nphi=nres-3
1327       nvar=ntheta+nphi
1328       nside=0
1329       do i=2,nres-1
1330         if (itype(i).ne.10 .and. itype(i).ne.21) then
1331           nside=nside+1
1332           ialph(i,1)=nvar+nside
1333           ialph(nside,2)=i
1334         endif
1335       enddo
1336       if (indphi.gt.0) then
1337         nvar=nphi
1338       else if (indback.gt.0) then
1339         nvar=nphi+ntheta
1340       else
1341         nvar=nvar+2*nside
1342       endif
1343 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1344       return
1345       end
1346 c----------------------------------------------------------------------------
1347       subroutine gen_dist_constr
1348 C Generate CA distance constraints.
1349       implicit real*8 (a-h,o-z)
1350       include 'DIMENSIONS'
1351       include 'COMMON.IOUNITS'
1352       include 'COMMON.GEO'
1353       include 'COMMON.VAR'
1354       include 'COMMON.INTERACT'
1355       include 'COMMON.LOCAL'
1356       include 'COMMON.NAMES'
1357       include 'COMMON.CHAIN'
1358       include 'COMMON.FFIELD'
1359       include 'COMMON.SBRIDGE'
1360       include 'COMMON.HEADER'
1361       include 'COMMON.CONTROL'
1362       include 'COMMON.DBASE'
1363       include 'COMMON.THREAD'
1364       include 'COMMON.TIME1'
1365       dimension itype_pdb(maxres)
1366       common /pizda/ itype_pdb
1367       character*2 iden
1368 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1369 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1370 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1371 cd     & ' nsup',nsup
1372       do i=nstart_sup,nstart_sup+nsup-1
1373 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1374 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1375         do j=i+2,nstart_sup+nsup-1
1376           nhpb=nhpb+1
1377           ihpb(nhpb)=i+nstart_seq-nstart_sup
1378           jhpb(nhpb)=j+nstart_seq-nstart_sup
1379           forcon(nhpb)=weidis
1380           dhpb(nhpb)=dist(i,j)
1381         enddo
1382       enddo
1383 cd      write (iout,'(a)') 'Distance constraints:'
1384 cd      do i=nss+1,nhpb
1385 cd        ii=ihpb(i)
1386 cd        jj=jhpb(i)
1387 cd        iden='CA'
1388 cd        if (ii.gt.nres) then
1389 cd          iden='SC'
1390 cd          ii=ii-nres
1391 cd          jj=jj-nres
1392 cd        endif
1393 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1394 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1395 cd     &  dhpb(i),forcon(i)
1396 cd      enddo
1397       return
1398       end
1399 c----------------------------------------------------------------------------
1400       subroutine map_read
1401       implicit real*8 (a-h,o-z)
1402       include 'DIMENSIONS'
1403       include 'COMMON.MAP'
1404       include 'COMMON.IOUNITS'
1405       character*3 angid(4) /'THE','PHI','ALP','OME'/
1406       character*80 mapcard,ucase
1407       do imap=1,nmap
1408         read (inp,'(a)') mapcard
1409         mapcard=ucase(mapcard)
1410         if (index(mapcard,'PHI').gt.0) then
1411           kang(imap)=1
1412         else if (index(mapcard,'THE').gt.0) then
1413           kang(imap)=2
1414         else if (index(mapcard,'ALP').gt.0) then
1415           kang(imap)=3
1416         else if (index(mapcard,'OME').gt.0) then
1417           kang(imap)=4
1418         else
1419           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1420           stop 'Error - illegal variable spec in MAP card.'
1421         endif
1422         call readi (mapcard,'RES1',res1(imap),0)
1423         call readi (mapcard,'RES2',res2(imap),0)
1424         if (res1(imap).eq.0) then
1425           res1(imap)=res2(imap)
1426         else if (res2(imap).eq.0) then
1427           res2(imap)=res1(imap)
1428         endif
1429         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1430           write (iout,'(a)')
1431      &    'Error - illegal definition of variable group in MAP.'
1432           stop 'Error - illegal definition of variable group in MAP.'
1433         endif
1434         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1435         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1436         call readi(mapcard,'NSTEP',nstep(imap),0)
1437         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1438           write (iout,'(a)')
1439      &     'Illegal boundary and/or step size specification in MAP.'
1440           stop 'Illegal boundary and/or step size specification in MAP.'
1441         endif
1442       enddo ! imap
1443       return
1444       end
1445 c----------------------------------------------------------------------------
1446       subroutine csaread
1447       implicit real*8 (a-h,o-z)
1448       include 'DIMENSIONS'
1449       include 'COMMON.IOUNITS'
1450       include 'COMMON.GEO'
1451       include 'COMMON.CSA'
1452       include 'COMMON.BANK'
1453       include 'COMMON.CONTROL'
1454       character*80 ucase
1455       character*620 mcmcard
1456       call card_concat(mcmcard)
1457
1458       call readi(mcmcard,'NCONF',nconf,50)
1459       call readi(mcmcard,'NADD',nadd,0)
1460       call readi(mcmcard,'JSTART',jstart,1)
1461       call readi(mcmcard,'JEND',jend,1)
1462       call readi(mcmcard,'NSTMAX',nstmax,500000)
1463       call readi(mcmcard,'N0',n0,1)
1464       call readi(mcmcard,'N1',n1,6)
1465       call readi(mcmcard,'N2',n2,4)
1466       call readi(mcmcard,'N3',n3,0)
1467       call readi(mcmcard,'N4',n4,0)
1468       call readi(mcmcard,'N5',n5,0)
1469       call readi(mcmcard,'N6',n6,10)
1470       call readi(mcmcard,'N7',n7,0)
1471       call readi(mcmcard,'N8',n8,0)
1472       call readi(mcmcard,'N9',n9,0)
1473       call readi(mcmcard,'N14',n14,0)
1474       call readi(mcmcard,'N15',n15,0)
1475       call readi(mcmcard,'N16',n16,0)
1476       call readi(mcmcard,'N17',n17,0)
1477       call readi(mcmcard,'N18',n18,0)
1478
1479       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1480
1481       call readi(mcmcard,'NDIFF',ndiff,2)
1482       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1483       call readi(mcmcard,'IS1',is1,1)
1484       call readi(mcmcard,'IS2',is2,8)
1485       call readi(mcmcard,'NRAN0',nran0,4)
1486       call readi(mcmcard,'NRAN1',nran1,2)
1487       call readi(mcmcard,'IRR',irr,1)
1488       call readi(mcmcard,'NSEED',nseed,20)
1489       call readi(mcmcard,'NTOTAL',ntotal,10000)
1490       call reada(mcmcard,'CUT1',cut1,2.0d0)
1491       call reada(mcmcard,'CUT2',cut2,5.0d0)
1492       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1493       call readi(mcmcard,'ICMAX',icmax,3)
1494       call readi(mcmcard,'IRESTART',irestart,0)
1495 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1496       ntbankm=0
1497 c!bankt
1498       call reada(mcmcard,'DELE',dele,20.0d0)
1499       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1500       call readi(mcmcard,'IREF',iref,0)
1501       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1502       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1503       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1504       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1505       write (iout,*) "NCONF_IN",nconf_in
1506       return
1507       end
1508 c----------------------------------------------------------------------------
1509 cfmc      subroutine mcmfread
1510 cfmc      implicit real*8 (a-h,o-z)
1511 cfmc      include 'DIMENSIONS'
1512 cfmc      include 'COMMON.MCMF'
1513 cfmc      include 'COMMON.IOUNITS'
1514 cfmc      include 'COMMON.GEO'
1515 cfmc      character*80 ucase
1516 cfmc      character*620 mcmcard
1517 cfmc      call card_concat(mcmcard)
1518 cfmc
1519 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1520 cfmc      write(iout,*)'MAXRANT=',maxrant
1521 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1522 cfmc      write(iout,*)'MAXFAM=',maxfam
1523 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1524 cfmc      write(iout,*)'NNET1=',nnet1
1525 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1526 cfmc      write(iout,*)'NNET2=',nnet2
1527 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1528 cfmc      write(iout,*)'NNET3=',nnet3
1529 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1530 cfmc      write(iout,*)'ILASTT=',ilastt
1531 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1532 cfmc      write(iout,*)'MAXSTR=',maxstr
1533 cfmc      maxstr_f=maxstr/maxfam
1534 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1535 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1536 cfmc      write(iout,*)'NMCMF=',nmcmf
1537 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1538 cfmc      write(iout,*)'IFOCUS=',ifocus
1539 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1540 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1541 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1542 cfmc      write(iout,*)'INTPRT=',intprt
1543 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1544 cfmc      write(iout,*)'IPRT=',iprt
1545 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1546 cfmc      write(iout,*)'IMAXTR=',imaxtr
1547 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1548 cfmc      write(iout,*)'MAXEVEN=',maxeven
1549 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1550 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1551 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1552 cfmc      write(iout,*)'INIMIN=',inimin
1553 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1554 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1555 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1556 cfmc      write(iout,*)'NTHREAD=',nthread
1557 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1558 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1559 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1560 cfmc      write(iout,*)'MAXPERT=',maxpert
1561 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1562 cfmc      write(iout,*)'IRMSD=',irmsd
1563 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1564 cfmc      write(iout,*)'DENEMIN=',denemin
1565 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1566 cfmc      write(iout,*)'RCUT1S=',rcut1s
1567 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1568 cfmc      write(iout,*)'RCUT1E=',rcut1e
1569 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1570 cfmc      write(iout,*)'RCUT2S=',rcut2s
1571 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1572 cfmc      write(iout,*)'RCUT2E=',rcut2e
1573 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1574 cfmc      write(iout,*)'DPERT1=',d_pert1
1575 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1576 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1577 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1578 cfmc      write(iout,*)'DPERT2=',d_pert2
1579 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1580 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1581 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1582 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1583 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1584 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1585 cfmc      d_pert1=deg2rad*d_pert1
1586 cfmc      d_pert1a=deg2rad*d_pert1a
1587 cfmc      d_pert2=deg2rad*d_pert2
1588 cfmc      d_pert2a=deg2rad*d_pert2a
1589 cfmc      d_pert2b=deg2rad*d_pert2b
1590 cfmc      d_pert2c=deg2rad*d_pert2c
1591 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1592 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1593 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1594 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1595 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1596 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1597 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1598 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1599 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1600 cfmc      write(iout,*)'RCUTINI=',rcutini
1601 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1602 cfmc      write(iout,*)'GRAT=',grat
1603 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1604 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1605 cfmc
1606 cfmc      return
1607 cfmc      end
1608 c----------------------------------------------------------------------------
1609       subroutine mcmread
1610       implicit real*8 (a-h,o-z)
1611       include 'DIMENSIONS'
1612       include 'COMMON.MCM'
1613       include 'COMMON.MCE'
1614       include 'COMMON.IOUNITS'
1615       character*80 ucase
1616       character*320 mcmcard
1617       call card_concat(mcmcard)
1618       call readi(mcmcard,'MAXACC',maxacc,100)
1619       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1620       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1621       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1622       call readi(mcmcard,'MAXREPM',maxrepm,200)
1623       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1624       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1625       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1626       call reada(mcmcard,'E_UP',e_up,5.0D0)
1627       call reada(mcmcard,'DELTE',delte,0.1D0)
1628       call readi(mcmcard,'NSWEEP',nsweep,5)
1629       call readi(mcmcard,'NSTEPH',nsteph,0)
1630       call readi(mcmcard,'NSTEPC',nstepc,0)
1631       call reada(mcmcard,'TMIN',tmin,298.0D0)
1632       call reada(mcmcard,'TMAX',tmax,298.0D0)
1633       call readi(mcmcard,'NWINDOW',nwindow,0)
1634       call readi(mcmcard,'PRINT_MC',print_mc,0)
1635       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1636       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1637       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1638       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1639       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1640       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1641       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1642       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1643       if (nwindow.gt.0) then
1644         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1645         do i=1,nwindow
1646           winlen(i)=winend(i)-winstart(i)+1
1647         enddo
1648       endif
1649       if (tmax.lt.tmin) tmax=tmin
1650       if (tmax.eq.tmin) then
1651         nstepc=0
1652         nsteph=0
1653       endif
1654       if (nstepc.gt.0 .and. nsteph.gt.0) then
1655         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1656         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1657       endif
1658 C Probabilities of different move types
1659       sumpro_type(0)=0.0D0
1660       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1661       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1662       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1663       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1664       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1665       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1666       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1667       do i=1,MaxMoveType
1668         print *,'i',i,' sumprotype',sumpro_type(i)
1669         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1670         print *,'i',i,' sumprotype',sumpro_type(i)
1671       enddo
1672       return
1673       end
1674 c----------------------------------------------------------------------------
1675       subroutine read_minim
1676       implicit real*8 (a-h,o-z)
1677       include 'DIMENSIONS'
1678       include 'COMMON.MINIM'
1679       include 'COMMON.IOUNITS'
1680       character*80 ucase
1681       character*320 minimcard
1682       call card_concat(minimcard)
1683       call readi(minimcard,'MAXMIN',maxmin,2000)
1684       call readi(minimcard,'MAXFUN',maxfun,5000)
1685       call readi(minimcard,'MINMIN',minmin,maxmin)
1686       call readi(minimcard,'MINFUN',minfun,maxmin)
1687       call reada(minimcard,'TOLF',tolf,1.0D-2)
1688       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1689       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1690       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1691       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1692       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1693      &         'Options in energy minimization:'
1694       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1695      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1696      & 'MinMin:',MinMin,' MinFun:',MinFun,
1697      & ' TolF:',TolF,' RTolF:',RTolF
1698       return
1699       end
1700 c----------------------------------------------------------------------------
1701       subroutine read_angles(kanal,*)
1702       implicit real*8 (a-h,o-z)
1703       include 'DIMENSIONS'
1704       include 'COMMON.GEO'
1705       include 'COMMON.VAR'
1706       include 'COMMON.CHAIN'
1707       include 'COMMON.IOUNITS'
1708       include 'COMMON.CONTROL'
1709 c Read angles from input
1710 c
1711        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1712        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1713        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1714        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1715
1716        do i=1,nres
1717 c 9/7/01 avoid 180 deg valence angle
1718         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1719 c
1720         theta(i)=deg2rad*theta(i)
1721         phi(i)=deg2rad*phi(i)
1722         alph(i)=deg2rad*alph(i)
1723         omeg(i)=deg2rad*omeg(i)
1724        enddo
1725       return
1726    10 return1
1727       end
1728 c----------------------------------------------------------------------------
1729       subroutine reada(rekord,lancuch,wartosc,default)
1730       implicit none
1731       character*(*) rekord,lancuch
1732       double precision wartosc,default
1733       integer ilen,iread
1734       external ilen
1735       iread=index(rekord,lancuch)
1736       if (iread.eq.0) then
1737         wartosc=default
1738         return
1739       endif
1740       iread=iread+ilen(lancuch)+1
1741       read (rekord(iread:),*,err=10,end=10) wartosc
1742       return
1743   10  wartosc=default
1744       return
1745       end
1746 c----------------------------------------------------------------------------
1747       subroutine readi(rekord,lancuch,wartosc,default)
1748       implicit none
1749       character*(*) rekord,lancuch
1750       integer wartosc,default
1751       integer ilen,iread
1752       external ilen
1753       iread=index(rekord,lancuch)
1754       if (iread.eq.0) then
1755         wartosc=default
1756         return
1757       endif
1758       iread=iread+ilen(lancuch)+1
1759       read (rekord(iread:),*,err=10,end=10) wartosc
1760       return
1761   10  wartosc=default
1762       return
1763       end
1764 c----------------------------------------------------------------------------
1765       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1766       implicit none
1767       integer dim,i
1768       integer tablica(dim),default
1769       character*(*) rekord,lancuch
1770       character*80 aux
1771       integer ilen,iread
1772       external ilen
1773       do i=1,dim
1774         tablica(i)=default
1775       enddo
1776       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1777       if (iread.eq.0) return
1778       iread=iread+ilen(lancuch)+1
1779       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1780    10 return
1781       end
1782 c----------------------------------------------------------------------------
1783       subroutine multreada(rekord,lancuch,tablica,dim,default)
1784       implicit none
1785       integer dim,i
1786       double precision tablica(dim),default
1787       character*(*) rekord,lancuch
1788       character*80 aux
1789       integer ilen,iread
1790       external ilen
1791       do i=1,dim
1792         tablica(i)=default
1793       enddo
1794       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1795       if (iread.eq.0) return
1796       iread=iread+ilen(lancuch)+1
1797       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1798    10 return
1799       end
1800 c----------------------------------------------------------------------------
1801       subroutine openunits
1802       implicit real*8 (a-h,o-z)
1803       include 'DIMENSIONS'
1804 #ifdef MPI
1805       include 'mpif.h'
1806       character*16 form,nodename
1807       integer nodelen
1808 #endif
1809       include 'COMMON.SETUP'
1810       include 'COMMON.IOUNITS'
1811       include 'COMMON.MD'
1812       include 'COMMON.CONTROL'
1813       integer lenpre,lenpot,ilen,lentmp
1814       external ilen
1815       character*3 out1file_text,ucase
1816       character*3 ll
1817       external ucase
1818 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1819       call getenv_loc("PREFIX",prefix)
1820       pref_orig = prefix
1821       call getenv_loc("POT",pot)
1822       call getenv_loc("DIRTMP",tmpdir)
1823       call getenv_loc("CURDIR",curdir)
1824       call getenv_loc("OUT1FILE",out1file_text)
1825 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1826       out1file_text=ucase(out1file_text)
1827       if (out1file_text(1:1).eq."Y") then
1828         out1file=.true.
1829       else
1830         out1file=fg_rank.gt.0
1831       endif
1832       lenpre=ilen(prefix)
1833       lenpot=ilen(pot)
1834       lentmp=ilen(tmpdir)
1835       if (lentmp.gt.0) then
1836           write (*,'(80(1h!))')
1837           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1838           write (*,'(80(1h!))')
1839           write (*,*)"All output files will be on node /tmp directory."
1840 #ifdef MPI
1841         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1842         if (me.eq.king) then
1843           write (*,*) "The master node is ",nodename
1844         else if (fg_rank.eq.0) then
1845           write (*,*) "I am the CG slave node ",nodename
1846         else
1847           write (*,*) "I am the FG slave node ",nodename
1848         endif
1849 #endif
1850         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1851         lenpre = lentmp+lenpre+1
1852       endif
1853       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1854 C Get the names and open the input files
1855 #if defined(WINIFL) || defined(WINPGI)
1856       open(1,file=pref_orig(:ilen(pref_orig))//
1857      &  '.inp',status='old',readonly,shared)
1858        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1859 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1860 C Get parameter filenames and open the parameter files.
1861       call getenv_loc('BONDPAR',bondname)
1862       open (ibond,file=bondname,status='old',readonly,shared)
1863       call getenv_loc('THETPAR',thetname)
1864       open (ithep,file=thetname,status='old',readonly,shared)
1865 #ifndef CRYST_THETA
1866       call getenv_loc('THETPARPDB',thetname_pdb)
1867       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1868 #endif
1869       call getenv_loc('ROTPAR',rotname)
1870       open (irotam,file=rotname,status='old',readonly,shared)
1871 #ifndef CRYST_SC
1872       call getenv_loc('ROTPARPDB',rotname_pdb)
1873       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1874 #endif
1875       call getenv_loc('TORPAR',torname)
1876       open (itorp,file=torname,status='old',readonly,shared)
1877       call getenv_loc('TORDPAR',tordname)
1878       open (itordp,file=tordname,status='old',readonly,shared)
1879       call getenv_loc('FOURIER',fouriername)
1880       open (ifourier,file=fouriername,status='old',readonly,shared)
1881       call getenv_loc('ELEPAR',elename)
1882       open (ielep,file=elename,status='old',readonly,shared)
1883       call getenv_loc('SIDEPAR',sidename)
1884       open (isidep,file=sidename,status='old',readonly,shared)
1885 #elif (defined CRAY) || (defined AIX)
1886       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1887      &  action='read')
1888 c      print *,"Processor",myrank," opened file 1"
1889       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1890 c      print *,"Processor",myrank," opened file 9"
1891 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1892 C Get parameter filenames and open the parameter files.
1893       call getenv_loc('BONDPAR',bondname)
1894       open (ibond,file=bondname,status='old',action='read')
1895 c      print *,"Processor",myrank," opened file IBOND"
1896       call getenv_loc('THETPAR',thetname)
1897       open (ithep,file=thetname,status='old',action='read')
1898 c      print *,"Processor",myrank," opened file ITHEP"
1899 #ifndef CRYST_THETA
1900       call getenv_loc('THETPARPDB',thetname_pdb)
1901       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1902 #endif
1903       call getenv_loc('ROTPAR',rotname)
1904       open (irotam,file=rotname,status='old',action='read')
1905 c      print *,"Processor",myrank," opened file IROTAM"
1906 #ifndef CRYST_SC
1907       call getenv_loc('ROTPARPDB',rotname_pdb)
1908       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1909 #endif
1910       call getenv_loc('TORPAR',torname)
1911       open (itorp,file=torname,status='old',action='read')
1912 c      print *,"Processor",myrank," opened file ITORP"
1913       call getenv_loc('TORDPAR',tordname)
1914       open (itordp,file=tordname,status='old',action='read')
1915 c      print *,"Processor",myrank," opened file ITORDP"
1916       call getenv_loc('SCCORPAR',sccorname)
1917       open (isccor,file=sccorname,status='old',action='read')
1918 c      print *,"Processor",myrank," opened file ISCCOR"
1919       call getenv_loc('FOURIER',fouriername)
1920       open (ifourier,file=fouriername,status='old',action='read')
1921 c      print *,"Processor",myrank," opened file IFOURIER"
1922       call getenv_loc('ELEPAR',elename)
1923       open (ielep,file=elename,status='old',action='read')
1924 c      print *,"Processor",myrank," opened file IELEP"
1925       call getenv_loc('SIDEPAR',sidename)
1926       open (isidep,file=sidename,status='old',action='read')
1927 c      print *,"Processor",myrank," opened file ISIDEP"
1928 c      print *,"Processor",myrank," opened parameter files"
1929 #elif (defined G77)
1930       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1931       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1932 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1933 C Get parameter filenames and open the parameter files.
1934       call getenv_loc('BONDPAR',bondname)
1935       open (ibond,file=bondname,status='old')
1936       call getenv_loc('THETPAR',thetname)
1937       open (ithep,file=thetname,status='old')
1938 #ifndef CRYST_THETA
1939       call getenv_loc('THETPARPDB',thetname_pdb)
1940       open (ithep_pdb,file=thetname_pdb,status='old')
1941 #endif
1942       call getenv_loc('ROTPAR',rotname)
1943       open (irotam,file=rotname,status='old')
1944 #ifndef CRYST_SC
1945       call getenv_loc('ROTPARPDB',rotname_pdb)
1946       open (irotam_pdb,file=rotname_pdb,status='old')
1947 #endif
1948       call getenv_loc('TORPAR',torname)
1949       open (itorp,file=torname,status='old')
1950       call getenv_loc('TORDPAR',tordname)
1951       open (itordp,file=tordname,status='old')
1952       call getenv_loc('SCCORPAR',sccorname)
1953       open (isccor,file=sccorname,status='old')
1954       call getenv_loc('FOURIER',fouriername)
1955       open (ifourier,file=fouriername,status='old')
1956       call getenv_loc('ELEPAR',elename)
1957       open (ielep,file=elename,status='old')
1958       call getenv_loc('SIDEPAR',sidename)
1959       open (isidep,file=sidename,status='old')
1960 #else
1961       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1962      &  readonly)
1963        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1964 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1965 C Get parameter filenames and open the parameter files.
1966       call getenv_loc('BONDPAR',bondname)
1967       open (ibond,file=bondname,status='old',readonly)
1968       call getenv_loc('THETPAR',thetname)
1969       open (ithep,file=thetname,status='old',readonly)
1970 #ifndef CRYST_THETA
1971       call getenv_loc('THETPARPDB',thetname_pdb)
1972       print *,"thetname_pdb ",thetname_pdb
1973       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1974       print *,ithep_pdb," opened"
1975 #endif
1976       call getenv_loc('ROTPAR',rotname)
1977       open (irotam,file=rotname,status='old',readonly)
1978 #ifndef CRYST_SC
1979       call getenv_loc('ROTPARPDB',rotname_pdb)
1980       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1981 #endif
1982       call getenv_loc('TORPAR',torname)
1983       open (itorp,file=torname,status='old',readonly)
1984       call getenv_loc('TORDPAR',tordname)
1985       open (itordp,file=tordname,status='old',readonly)
1986       call getenv_loc('SCCORPAR',sccorname)
1987       open (isccor,file=sccorname,status='old',readonly)
1988       call getenv_loc('FOURIER',fouriername)
1989       open (ifourier,file=fouriername,status='old',readonly)
1990       call getenv_loc('ELEPAR',elename)
1991       open (ielep,file=elename,status='old',readonly)
1992       call getenv_loc('SIDEPAR',sidename)
1993       open (isidep,file=sidename,status='old',readonly)
1994 #endif
1995 #ifndef OLDSCP
1996 C
1997 C 8/9/01 In the newest version SCp interaction constants are read from a file
1998 C Use -DOLDSCP to use hard-coded constants instead.
1999 C
2000       call getenv_loc('SCPPAR',scpname)
2001 #if defined(WINIFL) || defined(WINPGI)
2002       open (iscpp,file=scpname,status='old',readonly,shared)
2003 #elif (defined CRAY)  || (defined AIX)
2004       open (iscpp,file=scpname,status='old',action='read')
2005 #elif (defined G77)
2006       open (iscpp,file=scpname,status='old')
2007 #else
2008       open (iscpp,file=scpname,status='old',readonly)
2009 #endif
2010 #endif
2011       call getenv_loc('PATTERN',patname)
2012 #if defined(WINIFL) || defined(WINPGI)
2013       open (icbase,file=patname,status='old',readonly,shared)
2014 #elif (defined CRAY)  || (defined AIX)
2015       open (icbase,file=patname,status='old',action='read')
2016 #elif (defined G77)
2017       open (icbase,file=patname,status='old')
2018 #else
2019       open (icbase,file=patname,status='old',readonly)
2020 #endif
2021 #ifdef MPI
2022 C Open output file only for CG processes
2023 c      print *,"Processor",myrank," fg_rank",fg_rank
2024       if (fg_rank.eq.0) then
2025
2026       if (nodes.eq.1) then
2027         npos=3
2028       else
2029         npos = dlog10(dfloat(nodes-1))+1
2030       endif
2031       if (npos.lt.3) npos=3
2032       write (liczba,'(i1)') npos
2033       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2034      &  //')'
2035       write (liczba,form) me
2036       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2037      &  liczba(:ilen(liczba))
2038       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2039      &  //'.int'
2040       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2041      &  //'.pdb'
2042       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2043      &  liczba(:ilen(liczba))//'.mol2'
2044       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2045      &  liczba(:ilen(liczba))//'.stat'
2046       if (lentmp.gt.0)
2047      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2048      &      //liczba(:ilen(liczba))//'.stat')
2049       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2050      &  //'.rst'
2051       if(usampl) then
2052           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2053      & liczba(:ilen(liczba))//'.const'
2054       endif
2055
2056       endif
2057 #else
2058       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2059       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2060       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2061       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2062       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2063       if (lentmp.gt.0)
2064      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2065      &    //'.stat')
2066       rest2name=prefix(:ilen(prefix))//'.rst'
2067       if(usampl) then
2068          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2069       endif
2070 #endif
2071 #if defined(AIX) || defined(PGI)
2072       if (me.eq.king .or. .not. out1file)
2073      &   open(iout,file=outname,status='unknown')
2074 #ifdef DEBUG
2075       if (fg_rank.gt.0) then
2076         write (liczba,'(i3.3)') myrank/nfgtasks
2077         write (ll,'(bz,i3.3)') fg_rank
2078         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2079      &   status='unknown')
2080       endif
2081 #endif
2082       if(me.eq.king) then
2083        open(igeom,file=intname,status='unknown',position='append')
2084        open(ipdb,file=pdbname,status='unknown')
2085        open(imol2,file=mol2name,status='unknown')
2086        open(istat,file=statname,status='unknown',position='append')
2087       else
2088 c1out       open(iout,file=outname,status='unknown')
2089       endif
2090 #else
2091       if (me.eq.king .or. .not.out1file)
2092      &    open(iout,file=outname,status='unknown')
2093 #ifdef DEBUG
2094       if (fg_rank.gt.0) then
2095         write (liczba,'(i3.3)') myrank/nfgtasks
2096         write (ll,'(bz,i3.3)') fg_rank
2097         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2098      &   status='unknown')
2099       endif
2100 #endif
2101       if(me.eq.king) then
2102        open(igeom,file=intname,status='unknown',access='append')
2103        open(ipdb,file=pdbname,status='unknown')
2104        open(imol2,file=mol2name,status='unknown')
2105        open(istat,file=statname,status='unknown',access='append')
2106       else
2107 c1out       open(iout,file=outname,status='unknown')
2108       endif
2109 #endif
2110       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2111       csa_seed=prefix(:lenpre)//'.CSA.seed'
2112       csa_history=prefix(:lenpre)//'.CSA.history'
2113       csa_bank=prefix(:lenpre)//'.CSA.bank'
2114       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2115       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2116       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2117 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2118       csa_int=prefix(:lenpre)//'.int'
2119       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2120       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2121       csa_in=prefix(:lenpre)//'.CSA.in'
2122 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2123 C Write file names
2124       if (me.eq.king)then
2125       write (iout,'(80(1h-))')
2126       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2127       write (iout,'(80(1h-))')
2128       write (iout,*) "Input file                      : ",
2129      &  pref_orig(:ilen(pref_orig))//'.inp'
2130       write (iout,*) "Output file                     : ",
2131      &  outname(:ilen(outname))
2132       write (iout,*)
2133       write (iout,*) "Sidechain potential file        : ",
2134      &  sidename(:ilen(sidename))
2135 #ifndef OLDSCP
2136       write (iout,*) "SCp potential file              : ",
2137      &  scpname(:ilen(scpname))
2138 #endif
2139       write (iout,*) "Electrostatic potential file    : ",
2140      &  elename(:ilen(elename))
2141       write (iout,*) "Cumulant coefficient file       : ",
2142      &  fouriername(:ilen(fouriername))
2143       write (iout,*) "Torsional parameter file        : ",
2144      &  torname(:ilen(torname))
2145       write (iout,*) "Double torsional parameter file : ",
2146      &  tordname(:ilen(tordname))
2147       write (iout,*) "SCCOR parameter file : ",
2148      &  sccorname(:ilen(sccorname))
2149       write (iout,*) "Bond & inertia constant file    : ",
2150      &  bondname(:ilen(bondname))
2151       write (iout,*) "Bending parameter file          : ",
2152      &  thetname(:ilen(thetname))
2153       write (iout,*) "Rotamer parameter file          : ",
2154      &  rotname(:ilen(rotname))
2155       write (iout,*) "Threading database              : ",
2156      &  patname(:ilen(patname))
2157       if (lentmp.ne.0)
2158      &write (iout,*)" DIRTMP                          : ",
2159      &  tmpdir(:lentmp)
2160       write (iout,'(80(1h-))')
2161       endif
2162       return
2163       end
2164 c----------------------------------------------------------------------------
2165       subroutine card_concat(card)
2166       implicit real*8 (a-h,o-z)
2167       include 'DIMENSIONS'
2168       include 'COMMON.IOUNITS'
2169       character*(*) card
2170       character*80 karta,ucase
2171       external ilen
2172       read (inp,'(a)') karta
2173       karta=ucase(karta)
2174       card=' '
2175       do while (karta(80:80).eq.'&')
2176         card=card(:ilen(card)+1)//karta(:79)
2177         read (inp,'(a)') karta
2178         karta=ucase(karta)
2179       enddo
2180       card=card(:ilen(card)+1)//karta
2181       return
2182       end
2183 c----------------------------------------------------------------------------------
2184       subroutine readrst
2185       implicit real*8 (a-h,o-z)
2186       include 'DIMENSIONS'
2187       include 'COMMON.CHAIN'
2188       include 'COMMON.IOUNITS'
2189       include 'COMMON.MD'
2190       open(irest2,file=rest2name,status='unknown')
2191       read(irest2,*) totT,EK,potE,totE,t_bath
2192       do i=1,2*nres
2193          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2194       enddo
2195       do i=1,2*nres
2196          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2197       enddo
2198       if(usampl) then
2199              read (irest2,*) iset
2200       endif
2201       close(irest2)
2202       return
2203       end
2204 c---------------------------------------------------------------------------------
2205       subroutine read_fragments
2206       implicit real*8 (a-h,o-z)
2207       include 'DIMENSIONS'
2208 #ifdef MPI
2209       include 'mpif.h'
2210 #endif
2211       include 'COMMON.SETUP'
2212       include 'COMMON.CHAIN'
2213       include 'COMMON.IOUNITS'
2214       include 'COMMON.MD'
2215       include 'COMMON.CONTROL'
2216       read(inp,*) nset,nfrag,npair,nfrag_back
2217       if(me.eq.king.or..not.out1file)
2218      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2219      &  " nfrag_back",nfrag_back
2220       do iset=1,nset
2221          read(inp,*) mset(iset)
2222        do i=1,nfrag
2223          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2224      &     qinfrag(i,iset)
2225          if(me.eq.king.or..not.out1file)
2226      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2227      &     ifrag(2,i,iset), qinfrag(i,iset)
2228        enddo
2229        do i=1,npair
2230         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2231      &    qinpair(i,iset)
2232         if(me.eq.king.or..not.out1file)
2233      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2234      &    ipair(2,i,iset), qinpair(i,iset)
2235        enddo
2236        do i=1,nfrag_back
2237         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2238      &     wfrag_back(3,i,iset),
2239      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2240         if(me.eq.king.or..not.out1file)
2241      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2242      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2243        enddo
2244       enddo
2245       return
2246       end
2247 c-------------------------------------------------------------------------------
2248       subroutine read_dist_constr
2249       implicit real*8 (a-h,o-z)
2250       include 'DIMENSIONS'
2251 #ifdef MPI
2252       include 'mpif.h'
2253 #endif
2254       include 'COMMON.SETUP'
2255       include 'COMMON.CONTROL'
2256       include 'COMMON.CHAIN'
2257       include 'COMMON.IOUNITS'
2258       include 'COMMON.SBRIDGE'
2259       integer ifrag_(2,100),ipair_(2,100)
2260       double precision wfrag_(100),wpair_(100)
2261       character*500 controlcard
2262 c      write (iout,*) "Calling read_dist_constr"
2263 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2264 c      call flush(iout)
2265       call card_concat(controlcard)
2266       call readi(controlcard,"NFRAG",nfrag_,0)
2267       call readi(controlcard,"NPAIR",npair_,0)
2268       call readi(controlcard,"NDIST",ndist_,0)
2269       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2270       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2271       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2272       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2273       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2274 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2275 c      write (iout,*) "IFRAG"
2276 c      do i=1,nfrag_
2277 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2278 c      enddo
2279 c      write (iout,*) "IPAIR"
2280 c      do i=1,npair_
2281 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2282 c      enddo
2283       call flush(iout)
2284       do i=1,nfrag_
2285         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2286         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2287      &    ifrag_(2,i)=nstart_sup+nsup-1
2288 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2289         call flush(iout)
2290         if (wfrag_(i).gt.0.0d0) then
2291         do j=ifrag_(1,i),ifrag_(2,i)-1
2292           do k=j+1,ifrag_(2,i)
2293 c            write (iout,*) "j",j," k",k
2294             ddjk=dist(j,k)
2295             if (constr_dist.eq.1) then
2296             nhpb=nhpb+1
2297             ihpb(nhpb)=j
2298             jhpb(nhpb)=k
2299               dhpb(nhpb)=ddjk
2300             forcon(nhpb)=wfrag_(i)
2301             else if (constr_dist.eq.2) then
2302               if (ddjk.le.dist_cut) then
2303                 nhpb=nhpb+1
2304                 ihpb(nhpb)=j
2305                 jhpb(nhpb)=k
2306                 dhpb(nhpb)=ddjk
2307                 forcon(nhpb)=wfrag_(i)
2308               endif
2309             else
2310               nhpb=nhpb+1
2311               ihpb(nhpb)=j
2312               jhpb(nhpb)=k
2313               dhpb(nhpb)=ddjk
2314               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2315             endif
2316 #ifdef MPI
2317             if (.not.out1file .or. me.eq.king)
2318      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2319      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2320 #else
2321             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2322      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2323 #endif
2324           enddo
2325         enddo
2326         endif
2327       enddo
2328       do i=1,npair_
2329         if (wpair_(i).gt.0.0d0) then
2330         ii = ipair_(1,i)
2331         jj = ipair_(2,i)
2332         if (ii.gt.jj) then
2333           itemp=ii
2334           ii=jj
2335           jj=itemp
2336         endif
2337         do j=ifrag_(1,ii),ifrag_(2,ii)
2338           do k=ifrag_(1,jj),ifrag_(2,jj)
2339             nhpb=nhpb+1
2340             ihpb(nhpb)=j
2341             jhpb(nhpb)=k
2342             forcon(nhpb)=wpair_(i)
2343             dhpb(nhpb)=dist(j,k)
2344 #ifdef MPI
2345             if (.not.out1file .or. me.eq.king)
2346      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2347      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2348 #else
2349             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2350      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2351 #endif
2352           enddo
2353         enddo
2354         endif
2355       enddo
2356       do i=1,ndist_
2357         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2358         if (forcon(nhpb+1).gt.0.0d0) then
2359           nhpb=nhpb+1
2360           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2361 #ifdef MPI
2362           if (.not.out1file .or. me.eq.king)
2363      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2364      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2365 #else
2366           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2367      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2368 #endif
2369         endif
2370       enddo
2371       call flush(iout)
2372       return
2373       end
2374 c-------------------------------------------------------------------------------
2375 #ifdef WINIFL
2376       subroutine flush(iu)
2377       return
2378       end
2379 #endif
2380 #ifdef AIX
2381       subroutine flush(iu)
2382       call flush_(iu)
2383       return
2384       end
2385 #endif
2386 c------------------------------------------------------------------------------
2387       subroutine copy_to_tmp(source)
2388       include "DIMENSIONS"
2389       include "COMMON.IOUNITS"
2390       character*(*) source
2391       character* 256 tmpfile
2392       integer ilen
2393       external ilen
2394       logical ex
2395       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2396       inquire(file=tmpfile,exist=ex)
2397       if (ex) then
2398         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2399      &   " to temporary directory..."
2400         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2401         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2402       endif
2403       return
2404       end
2405 c------------------------------------------------------------------------------
2406       subroutine move_from_tmp(source)
2407       include "DIMENSIONS"
2408       include "COMMON.IOUNITS"
2409       character*(*) source
2410       integer ilen
2411       external ilen
2412       write (*,*) "Moving ",source(:ilen(source)),
2413      & " from temporary directory to working directory"
2414       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2415       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2416       return
2417       end
2418 c------------------------------------------------------------------------------
2419       subroutine random_init(seed)
2420 C
2421 C Initialize random number generator
2422 C
2423       implicit real*8 (a-h,o-z)
2424       include 'DIMENSIONS'
2425 #ifdef AMD64
2426       integer*8 iseedi8
2427 #endif
2428 #ifdef MPI
2429       include 'mpif.h'
2430       logical OKRandom, prng_restart
2431       real*8  r1
2432       integer iseed_array(4)
2433 #endif
2434       include 'COMMON.IOUNITS'
2435       include 'COMMON.TIME1'
2436       include 'COMMON.THREAD'
2437       include 'COMMON.SBRIDGE'
2438       include 'COMMON.CONTROL'
2439       include 'COMMON.MCM'
2440       include 'COMMON.MAP'
2441       include 'COMMON.HEADER'
2442       include 'COMMON.CSA'
2443       include 'COMMON.CHAIN'
2444       include 'COMMON.MUCA'
2445       include 'COMMON.MD'
2446       include 'COMMON.FFIELD'
2447       include 'COMMON.SETUP'
2448       iseed=-dint(dabs(seed))
2449       if (iseed.eq.0) then
2450         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2451      &    'Random seed undefined. The program will stop.'
2452         write (*,'(/80(1h*)/20x,a/80(1h*))')
2453      &    'Random seed undefined. The program will stop.'
2454 #ifdef MPI
2455         call mpi_finalize(mpi_comm_world,ierr)
2456 #endif
2457         stop 'Bad random seed.'
2458       endif
2459 #ifdef MPI
2460       if (fg_rank.eq.0) then
2461       seed=seed*(me+1)+1
2462 #ifdef AMD64
2463       iseedi8=dint(seed)
2464       if(me.eq.king .or. .not. out1file)
2465      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2466       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2467       OKRandom = prng_restart(me,iseedi8)
2468 #else
2469       do i=1,4
2470        tmp=65536.0d0**(4-i)
2471        iseed_array(i) = dint(seed/tmp)
2472        seed=seed-iseed_array(i)*tmp
2473       enddo
2474       if(me.eq.king)
2475      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2476      &                 (iseed_array(i),i=1,4)
2477       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2478      &                 (iseed_array(i),i=1,4)
2479       OKRandom = prng_restart(me,iseed_array)
2480 #endif
2481       if (OKRandom) then
2482 c        r1 = prng_next(me)
2483         r1=ran_number(0.0D0,1.0D0)
2484         if(me.eq.king)
2485      &   write (iout,*) 'ran_num',r1
2486         if (r1.lt.0.0d0) OKRandom=.false.
2487       endif
2488       if (.not.OKRandom) then
2489         write (iout,*) 'PRNG IS NOT WORKING!!!'
2490         print *,'PRNG IS NOT WORKING!!!'
2491         if (me.eq.0) then
2492          call flush(iout)
2493          call mpi_abort(mpi_comm_world,error_msg,ierr)
2494          stop
2495         else
2496          write (iout,*) 'too many processors for parallel prng'
2497          write (*,*) 'too many processors for parallel prng'
2498          call flush(iout)
2499          stop
2500         endif
2501       endif
2502       endif
2503 #else
2504       call vrndst(iseed)
2505       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2506 #endif
2507       return
2508       end